工場、倉庫、販売店割り当てモデル: 問題ベース

ライブスクリプトを開く

この例では、混合整数線形計画問題を設定および解決する方法を説明します。問題は、一連の工場、倉庫、販売店間の生産および流通の最適レベルを求めることです。ソルバーベースのアプローチについては、工場、倉庫、販売店割り当てモデル: ソルバーベースを参照してください。

この例では、まず工場、倉庫、販売店のランダムな位置を生成します。スケーリングパラメーター $N$ を自由に変更して、生産および流通施設が配置されている両方のグリッドのサイズをスケーリングすると共に、グリッドエリアの各タイプの施設の密度が $N$ に依存しないように、これらの施設の数をスケーリングします。

施設の位置

特定のスケーリングパラメーターの値 $N$ には以下があるとします。

$⌊ f N^{2} ⌋$ 個の工場
$⌊ w N^{2} ⌋$ 個の倉庫
$⌊ s N^{2} ⌋$ 個の販売店

これらの施設は、 $x$ および $y$ 方向で 1 と $N$ の間にある別個の整数グリッド点にあります。施設が別個の位置をもつためには、 $f + w + s \leq 1$ である必要があります。この例では、 $N = 20$ 、 $f = 0.05$ 、 $w = 0.05$ 、 $s = 0.1$ とします。

生産および流通

工場では、 $P$ 個の製品が製造されています。 $P = 20$ とします。

販売店 $s$ における各製品 $p$ の需要は $d (s, p)$ です。需要とは、ある期間において販売可能な数量です。このモデルの制約の 1 つは需要を満たすことです。つまり、システムが需要に厳密に見合う数量を生産および供給することです。

各工場および倉庫には能力の制約があります。

工場 $f$ における製品 $p$ の生産は $p c a p (f, p)$ 未満です。
倉庫 $w$ の容量は $w c a p (w)$ です。
特定の期間に倉庫 $w$ から販売店に輸送できる製品 $p$ の数量は $t u r n (p) * w c a p (w)$ 未満です。ここで、 $t u r n (p)$ は製品 $p$ の回転率です。

各販売店は 1 つの倉庫からのみ供給を受けるものとします。問題の一部として、倉庫に対し販売店を最も低コストに配置することがあります。

コスト

工場から倉庫および倉庫から販売店までの製品輸送コストは、施設間の距離および製品ごとに異なります。 $d i s t (a, b)$ が施設 $a$ と $b$ の間の距離であるとすると、これらの施設間で製品 $p$ を出荷するコストは輸送コスト $t c o s t (p)$ に距離を乗算した値になります。

$d i s t (a, b) * t c o s t (p) .$

この例における距離は $L_{1}$ 距離とも呼ばれるグリッド距離です。これは $x$ 座標と $y$ 座標の差の絶対値の総和です。

工場 $f$ で製品 $p$ を 1 単位生産するコストは $p c o s t (f, p)$ です。

最適化問題

一連の施設の位置と需要および能力の制約を所与として、以下を求めます。

各工場における各製品の生産レベル
工場から倉庫への製品の出荷スケジュール
倉庫から販売店への製品の出荷スケジュール

これらの数量は、需要が満たされ、総コストが最小化されるものでなければなりません。また、各販売店はすべての製品を 1 つの倉庫からのみ受け取る必要があります。

最適化問題の変数と方程式

制御変数 (最適化において変更可能な変数) は以下のとおりです。

$x (p, f, w)$ = 工場 $f$ から倉庫 $w$ に輸送する製品 $p$ の数量
$y (s, w)$ = 販売店 $s$ が倉庫 $w$ と関連付けられているときに値が 1 になるバイナリ変数

最小化する目的関数は以下のとおりです。

$\sum_{f} \sum_{p} \sum_{w} x (p, f, w) \cdot (p c o s t (f, p) + t c o s t (p) \cdot d i s t (f, w))$

$+ \sum_{s} \sum_{w} \sum_{p} (d (s, p) \cdot t c o s t (p) \cdot d i s t (s, w) \cdot y (s, w)) .$

制約は次のとおりです。

$\sum_{w} x (p, f, w) \leq p c a p (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)}{t u r n (p)} \cdot y (s, w) \leq w c a p (w)$ (倉庫の容量)。

$\sum_{w} y (s, w) = 1$ (各販売店が 1 つの倉庫に関連付けられている)。

$x (p, f, w) \geq 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;

これらのランダムな需要と能力によって実行不可能な問題につながることがあります。つまり、需要が生産と倉庫の能力を超える場合があります。パラメーターを変更して実行不可能な問題が発生した場合、解を生成する際に終了フラグ -2 が表示されます。

変数および制約の生成

問題の指定を開始するには、距離の配列 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

最適化問題の変数を作成します。x は生産、つまり、次元が P x F x W の連続変数を表します。y は販売店の倉庫へのバイナリ割り当て、つまり、S x W の変数を表します。

x = optimvar('x',P,F,W,'LowerBound',0);
y = optimvar('y',S,W,'Type','integer','LowerBound',0,'UpperBound',1);

ここで制約を作成します。最初の制約は生産能力の制約です。

capconstr = sum(x,3) <= pcap';

次の制約は各販売店における需要の充足です。

demconstr = squeeze(sum(x,2)) == d'*y;

各倉庫には能力の制約があります。

warecap = sum(diag(1./turn)*(d'*y),1) <= wcap';

最後に、各販売店が正しく 1 つの倉庫に結びついている必要があります。

salesware = sum(y,2) == ones(S,1);

問題および目的の作成

最適化問題を作成します。

factoryprob = optimproblem;

目的関数には 3 つの部分があります。最初の部分は、生産コストの和です。

objfun1 = sum(sum(sum(x,3).*(pcost'),2),1);

2 番目の部分は、工場から倉庫への輸送コストの和です。

objfun2 = 0;
for p = 1:P
    objfun2 = objfun2 + tcost(p)*sum(sum(squeeze(x(p,:,:)).*distfw));
end

3 番目の部分は、倉庫から販売店への輸送コストの和です。

r = sum(distsw.*y,2); % r is a length s vector
v = d*(tcost(:));
objfun3 = sum(v.*r);

最小化する目的関数は、この 3 つの部分の和です。

factoryprob.Objective = objfun1 + objfun2 + objfun3;

制約を問題に含めます。

factoryprob.Constraints.capconstr = capconstr;
factoryprob.Constraints.demconstr = demconstr;
factoryprob.Constraints.warecap = warecap;
factoryprob.Constraints.salesware = salesware;

問題を解く

何百もの出力の線が表示されないように、反復表示をオフにします。解決進行状況を監視するためのプロット関数を含めます。プロット関数を含めるために、'legacy' アルゴリズムを指定します。

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

ソルバーを呼び出して解を求めます。

[sol,fval,exitflag,output] = solve(factoryprob,'options',opts);

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

解の検証

終了フラグと解の実行不可能性を調べます。

exitflag

exitflag = 
    OptimalSolution

infeas1 = max(max(infeasibility(capconstr,sol)))

infeas1 = 2.3647e-11

infeas2 = max(max(infeasibility(demconstr,sol)))

infeas2 = 2.2737e-13

infeas3 = max(infeasibility(warecap,sol))

infeas3 = 0

infeas4 = max(infeasibility(salesware,sol))

infeas4 = 8.8818e-16

厳密に整数値になるように解の y 部分を丸めます。これらの変数が厳密に整数ではない可能性がある理由を理解するには、一部の “整数” 解は整数ではないを参照してください。

sol.y = round(sol.y); % get integer solutions

各倉庫に関連付けられている販売店数はいくつでしょうか。この場合、一部の倉庫は関連付けられた店舗がゼロであることに注意してください。これは、倉庫が最適解で使用されていないことを示します。

outlets = sum(sol.y,1)

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(sol.y(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')

緑色の線がない黒い * は使用されない倉庫を表します。