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

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

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

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

Есть $$P$$ продукты, произведенная заводами. Бери $$P=20$$.

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

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

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

Затраты

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

Расстояние в этом примере является расстоянием сетки, также известным как $$L_1$$ расстояние. Это сумма абсолютного различия в $$x$$ координаты и $$y$$ координаты.

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

Задача оптимизации

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

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

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

Переменные управления, означающие те, которые вы можете изменить в оптимизации,

Целевая функция для минимизации

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

$$\sum _{w}x(p,f,w)\le pcap(f,p)$$ (мощность завода).

$$\sum _{f}x(p,f,w)=\sum _s(d(s,p)\cdot y(s,w))$$ (спрос удовлетворяется).

$$\sum _p\sum _s\frac{d(s,p)}{turn(p)}\cdot y(s,w)\le wcap(w)$$ (вместимость склада).

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

$$x(p,f,w)\ge 0$$ (неотрицательная продукция).

$$y(s,w)ϵ\{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);

На протяжении генерации целевых и ограничительных векторов и матриц мы генерируем $$(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')

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