Расположение объекта

Для заданного значения параметра масштабирования $$N$$предположим, что существуют следующие значения:

Эти средства находятся на отдельных целочисленных точках сетки между 1 и $$N$$ в $$$$$$$$направлениях x и y. Для того, чтобы объекты имели отдельные местоположения, требуется этот $$\mathrm{f+w+s\le 1}$$. В этом примере принимают $$N\mathrm{=}\mathrm{}$$20$$,\mathrm{f}\mathrm{=}\mathrm{}$$$$0,05\mathrm{,}\mathrm{w}\mathrm{}$$= 0,05 $$\mathrm{}\mathrm{и}$$s = 0,1.

Производство и распределение

Есть $$$$продукты P, изготовленные фабриками. Принять $$P\mathrm{=}\mathrm{}$$20.

Спрос на каждый продукт $$p$$ в торговой точке $$s$$ равен $$d(s,$$p). Потребность - это количество, которое может быть продано в интервале времени. Одним из ограничений модели является удовлетворение потребности, что означает, что система производит и распределяет именно те количества потребности.

На каждом заводе и на каждом складе существуют ограничения по мощности.

Предположим, что каждая торговая точка получает свои поставки только с одного склада. Отчасти проблема заключается в определении наиболее дешевого картирования торговых точек на склады.

Затраты

Стоимость транспортировки продукции с завода на склад и со склада на торговую точку зависит от расстояния между объектами и от конкретного продукта. Если $$\mathrm{dist}(a,$$b) - расстояние между $$\mathrm{объектами}$$ a $$и$$b, то стоимость транспортировки $$\mathrm{продукта}$$ p между этими объектами является расстоянием, умноженным на стоимость транспортировки $$\mathrm{tcost}$$(p):

Расстояние в этом примере - это расстояние сетки, также называемое $${}_{\mathrm{L1}}$$ расстоянием. Это сумма абсолютных разностей координат $$x$$ и $$y$$.

Стоимость изготовления единицы продукта $$p$$ на заводе $$f$$ составляет $$\mathrm{pcost}(f,$$p).

Проблема оптимизации

Учитывая набор местоположений объекта, а также требования и ограничения по мощности, найдите:

Эти количества должны обеспечивать удовлетворение потребности и минимизацию общих затрат. Также каждая торговая точка обязана получать всю свою продукцию ровно с одного склада.

Переменные и уравнения для задачи оптимизации

Управляющие переменные, то есть переменные, которые можно изменить при оптимизации:

Целевой функцией для минимизации является

Ограничения:

$$\underset{}{}\mathrm{\sum wx}(p,f,w)\mathrm{\le pcap}$$(f, p) (производительность завода).

$$\underset{}{}\mathrm{\sum fx}(p,f,\underset{)}{w}\mathrm{=\sum s}(d(s,p)$$⋅y (s, w)) (спрос удовлетворяется).

$$\underset{}{}\underset{}{}\frac{\mathrm{\sum p\sum sd}(s,}{p)поворот}(p)\cdot y(s,w)$$≤wcap (w) (вместимость склада).

$$\underset{}{}\mathrm{\sum wy}(ы,w\mathrm{)}$$ = 1 (каждая торговая точка связана с одним складом).

$$x(p,f,\mathrm{w}$$) ≥0 (неотрицательное производство).

$$y(s,w)\mathrm{}\mathrm{}$$ϵ{0,1} ($$$$двоичный y).

Переменные $$x$$ и $$y$$ отображаются в функциях цели и ограничения линейно. Поскольку $$y$$ ограничен целочисленными значениями, проблема заключается в смешанно-целочисленной линейной программе (MILP).

Создать случайную проблему: Расположение объекта

Задайте значения $$$$$$$$$$$$$$$$параметров N, f, w и s и создайте местоположения оборудования.

rng(1) % for reproducibility
N = 20; % N from 10 to 30 seems to work. Choose large values with caution.
N2 = N*N;
f = 0.05; % density of factories
w = 0.05; % density of warehouses
s = 0.1; % density of sales outlets

F = floor(f*N2); % number of factories
W = floor(w*N2); % number of warehouses
S = floor(s*N2); % number of sales outlets

xyloc = randperm(N2,F+W+S); % unique locations of facilities
[xloc,yloc] = ind2sub([N N],xyloc);

Конечно, нереально брать случайные места для объектов. Этот пример предназначен для демонстрации методов решения, а не способов создания хороших местоположений объекта.

Постройте график объектов. Объекты 1 - F - заводы, F + 1 - F + W - склады, а F + W + 1 - F + W + S - торговые точки.

h = figure;
plot(xloc(1:F),yloc(1:F),'rs',xloc(F+1:F+W),yloc(F+1:F+W),'k*',...
    xloc(F+W+1:F+W+S),yloc(F+W+1:F+W+S),'bo');
lgnd = legend('Factory','Warehouse','Sales outlet','Location','EastOutside');
lgnd.AutoUpdate = 'off';
xlim([0 N+1]);ylim([0 N+1])

Создание случайных мощностей, затрат и потребностей

Создание случайных производственных затрат, мощностей, оборотов и потребностей.

P = 20; % 20 products

% Production costs between 20 and 100
pcost = 80*rand(F,P) + 20;

% Production capacity between 500 and 1500 for each product/factory
pcap = 1000*rand(F,P) + 500;

% Warehouse capacity between P*400 and P*800 for each product/warehouse
wcap = P*400*rand(W,1) + P*400;

% Product turnover rate between 1 and 3 for each product
turn = 2*rand(1,P) + 1;

% Product transport cost per distance between 5 and 10 for each product
tcost = 5*rand(1,P) + 5;

% Product demand by sales outlet between 200 and 500 for each
% product/outlet
d = 300*rand(S,P) + 200;

Эти случайные потребности и возможности могут привести к неосуществимым проблемам. Другими словами, иногда спрос превышает производственные и складские ограничения. Если изменить некоторые параметры и получить неосуществимую проблему, во время решения вы получите exitflag -2.

Создание матриц целей и ограничений и векторов

Вектор целевой функции obj в intlincon состоит из коэффициентов переменных $$x(p,f,$$ w) $$иy($$s, w). Так что есть естественноP*F*W + S*W коэффициенты в obj.

Один из способов генерации коэффициентов - начать с P-by-F-by-W множество obj1 для $$$$коэффициентов x и S-by-W множество obj2 для $$$$коэффициентов y (s, w). Затем преобразовать эти массивы в два вектора и объединить их вobj путем вызова

obj = [obj1(:);obj2(:)];

obj1 = zeros(P,F,W); % Allocate arrays
obj2 = zeros(S,W);

На протяжении генерации векторов цели и ограничений и матриц мы генерируем $$\mathrm{массив\; (p},f,$$w) или $$$$массив (s, w), а затем преобразуем результат в вектор.

Чтобы начать генерирование входных данных, создайте массивы расстояний distfw(i,j) и distsw(i,j).

distfw = zeros(F,W); % Allocate matrix for factory-warehouse distances
for ii = 1:F
    for jj = 1:W
        distfw(ii,jj) = abs(xloc(ii) - xloc(F + jj)) + abs(yloc(ii) ...
            - yloc(F + jj));
    end
end

distsw = zeros(S,W); % Allocate matrix for sales outlet-warehouse distances
for ii = 1:S
    for jj = 1:W
        distsw(ii,jj) = abs(xloc(F + W + ii) - xloc(F + jj)) ...
            + abs(yloc(F + W + ii) - yloc(F + jj));
    end
end

Создание записей obj1 и obj2.

for ii = 1:P
    for jj = 1:F
        for kk = 1:W
            obj1(ii,jj,kk) = pcost(jj,ii) + tcost(ii)*distfw(jj,kk);
        end
    end
end

for ii = 1:S
    for jj = 1:W
        obj2(ii,jj) = distsw(ii,jj)*sum(d(ii,:).*tcost);
    end
end

Объединение записей в один вектор.

obj = [obj1(:);obj2(:)]; % obj is the objective function vector

Теперь создайте матрицы ограничений.

Ширина каждой матрицы линейных ограничений - это длина obj вектор.

matwid = length(obj);

Существует два типа линейных неравенств: ограничения производственных мощностей и ограничения складских мощностей.

Есть P*F ограничения производственного потенциала, и W ограничения емкости склада. Матрицы ограничений довольно разрежены, порядка 1% ненулевых, поэтому сохраните память, используя разреженные матрицы.

Aineq = spalloc(P*F + W,matwid,P*F*W + S*W); % Allocate sparse Aeq
bineq = zeros(P*F + W,1); % Allocate bineq as full

% Zero matrices of convenient sizes:
clearer1 = zeros(size(obj1));
clearer12 = clearer1(:);
clearer2 = zeros(size(obj2));
clearer22 = clearer2(:);

% First the production capacity constraints
counter = 1;
for ii = 1:F
    for jj = 1:P
        xtemp = clearer1;
        xtemp(jj,ii,:) = 1; % Sum over warehouses for each product and factory
        xtemp = sparse([xtemp(:);clearer22]); % Convert to sparse
        Aineq(counter,:) = xtemp'; % Fill in the row
        bineq(counter) = pcap(ii,jj);
        counter = counter + 1;
    end
end

% Now the warehouse capacity constraints
vj = zeros(S,1); % The multipliers 
for jj = 1:S
    vj(jj) = sum(d(jj,:)./turn); % A sum of P elements
end

for ii = 1:W
    xtemp = clearer2;
    xtemp(:,ii) = vj;
    xtemp = sparse([clearer12;xtemp(:)]); % Convert to sparse
    Aineq(counter,:) = xtemp'; % Fill in the row
    bineq(counter) = wcap(ii);
    counter = counter + 1;
end

Существует два типа линейных ограничений равенства: ограничение, что спрос удовлетворяется, и ограничение, что каждая торговая точка соответствует одному складу.

Aeq = spalloc(P*W + S,matwid,P*W*(F+S) + S*W); % Allocate as sparse
beq = zeros(P*W + S,1); % Allocate vectors as full

counter = 1;
% Demand is satisfied:
for ii = 1:P
    for jj = 1:W
        xtemp = clearer1;
        xtemp(ii,:,jj) = 1;
        xtemp2 = clearer2;
        xtemp2(:,jj) = -d(:,ii);
        xtemp = sparse([xtemp(:);xtemp2(:)]'); % Change to sparse row
        Aeq(counter,:) = xtemp; % Fill in row
        counter = counter + 1;
    end
end

% Only one warehouse for each sales outlet:
for ii = 1:S
    xtemp = clearer2;
    xtemp(ii,:) = 1;
    xtemp = sparse([clearer12;xtemp(:)]'); % Change to sparse row
    Aeq(counter,:) = xtemp; % Fill in row
    beq(counter) = 1;
    counter = counter + 1;
end

Связанные ограничения и целочисленные переменные

Целочисленные переменные: length(obj1) + 1 до конца.

intcon = P*F*W+1:length(obj);

Верхние границы - от length(obj1) + 1 до конца тоже.

lb = zeros(length(obj),1);
ub = Inf(length(obj),1);
ub(P*F*W+1:end) = 1;

Отключите итеративное отображение, чтобы не получать сотни строк выходных данных. Включите функцию графика для контроля выполнения решения.

opts = optimoptions('intlinprog','Display','off','PlotFcn',@optimplotmilp);

Решить проблему

Сгенерированы все входные данные решателя. Вызовите решатель, чтобы найти решение.

[solution,fval,exitflag,output] = intlinprog(obj,intcon,...
                                     Aineq,bineq,Aeq,beq,lb,ub,opts);

if isempty(solution) % If the problem is infeasible or you stopped early with no solution
    disp('intlinprog did not return a solution.')
    return % Stop the script because there is nothing to examine
end

Анализ решения

Решение возможно в пределах заданных допусков.

exitflag

exitflag = 1

infeas1 = max(Aineq*solution - bineq)

infeas1 = 8.2991e-12

infeas2 = norm(Aeq*solution - beq,Inf)

infeas2 = 1.6428e-11

Проверьте, что целочисленные компоненты действительно целые числа, или достаточно близки, чтобы было разумно округлить их. Чтобы понять, почему эти переменные могут быть не совсем целыми числами, см. Некоторые «целочисленные» решения не целые числа.

diffint = norm(solution(intcon) - round(solution(intcon)),Inf)

diffint = 1.1990e-13

Некоторые целочисленные переменные не являются целыми числами, но все очень близки. Так что округляем целочисленные переменные.

solution(intcon) = round(solution(intcon));

Проверьте выполнимость округленного решения и изменение значения целевой функции.

infeas1 = max(Aineq*solution - bineq)

infeas1 = 8.2991e-12

infeas2 = norm(Aeq*solution - beq,Inf)

infeas2 = 5.8435e-11

diffrounding = norm(fval - obj(:)'*solution,Inf)

diffrounding = 1.8626e-08

Округление решения не существенно изменило его осуществимость.

Проще всего изучить решение, изменив его исходные размеры.

solution1 = solution(1:P*F*W); % The continuous variables
solution2 = solution(intcon); % The integer variables
solution1 = reshape(solution1,P,F,W);
solution2 = reshape(solution2,S,W);

Например, сколько торговых точек связано с каждым складом? Обратите внимание, что в этом случае некоторые склады имеют 0 связанных торговых точек, что означает, что склады не используются в оптимальном решении.

outlets = sum(solution2,1) % Sum over the sales outlets

outlets = 1×20

     3     0     3     2     2     2     3     2     3     1     1     0     0     3     4     3     2     3     2     1

Постройте график соединения между каждой торговой точкой и ее складом.

figure(h);
hold on
for ii = 1:S
    jj = find(solution2(ii,:)); % Index of warehouse associated with ii
    xsales = xloc(F+W+ii); ysales = yloc(F+W+ii);
    xwarehouse = xloc(F+jj); ywarehouse = yloc(F+jj);
    if rand(1) < .5 % Draw y direction first half the time
        plot([xsales,xsales,xwarehouse],[ysales,ywarehouse,ywarehouse],'g--')
    else % Draw x direction first the rest of the time
        plot([xsales,xwarehouse,xwarehouse],[ysales,ysales,ywarehouse],'g--')
    end
end
hold off

title('Mapping of sales outlets to warehouses')

Черный * без зеленых линий представляет неиспользуемые склады.