Фабрика, склад, модель выделения продаж: основанный на решателе

В этом примере показано, как настроить и решить смешано-целочисленную задачу линейного программирования. Проблема состоит в том, чтобы найти оптимальные уровни производства и распределения среди набора фабрик, складов и торговых точек. Для подхода, основанного на проблеме смотрите Фабрику, Склад, Модель Выделения Продаж: основанный на проблеме.

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

Местоположения средства

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

  • fN2 фабрики

  • wN2 склады

  • sN2 торговые точки

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

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

Существуют P продукты сделаны фабриками. Взять P=20.

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

Существуют полные ограничения на каждую фабрику и каждый склад.

  • Производство продукта p на фабрике f меньше pcap(f,p).

  • Способность склада w wcap(w).

  • Сумма продукта p это может быть транспортировано со склада w к торговой точке во временном интервале меньше turn(p)*wcap(w), где turn(p) текучесть кадров продукта p.

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

Затраты

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

dist(a,b)*tcost(p).

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

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

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

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

  • Производственный уровень каждого продукта на каждой фабрике

  • Расписание распределения для продуктов от фабрик до складов

  • Расписание распределения для продуктов от складов до торговых точек

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

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

Контрольные переменные, означая тех, которых можно изменить в оптимизации,

  • x(p,f,w) = сумма продукта p это транспортируется от фабрики f в склад w

  • y(s,w) = бинарное переменное принимающее значение 1, когда торговая точка s сопоставлен со складом w

Целевая функция, чтобы минимизировать

fpwx(p,f,w)(pcost(f,p)+tcost(p)dist(f,w))

+swp(d(s,p)tcost(p)dist(s,w)y(s,w)).

Ограничения

wx(p,f,w)pcap(f,p) (способность фабрики).

fx(p,f,w)=s(d(s,p)y(s,w)) (требование соблюдается).

psd(s,p)turn(p)y(s,w)wcap(w) (способность склада).

wy(s,w)=1 (каждая торговая точка сопоставляет в один склад).

x(p,f,w)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])

Figure contains an axes. The axes contains 3 objects of type line. These objects represent Factory, Warehouse, Sales outlet.

Сгенерируйте случайные мощности, затраты и требования

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

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);

Figure Optimization Plot Function contains an axes. The axes with title Best objective: 3.0952e+07, Relative gap: 0. contains 4 objects of type line. These objects represent Root LB, Cuts LB, Heuristics UB, New Solution.

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')

Figure contains an axes. The axes with title Mapping of sales outlets to warehouses contains 43 objects of type line. These objects represent Factory, Warehouse, Sales outlet.

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

Похожие темы