【線形計画法②】シンプレックス法でトラックの積載率が最大になる条件を求めてみた
貨物が３種類以上のトラック積載率最大化問題は一次方程式では解けない
トラックには容積と重量という２つの制約があり、積載量が限られます。
そして、貨物の比重によってどちらの制約に引っかかるかが決まります。
布団のような軽い貨物ばかりを積めば、荷台は満杯でも重量は半分くらい余裕があるかもしれません。
逆に、水のように重たい貨物ばかりを積めば、重量ギリギリまで積んでも荷台は半分くらい余ってしまうかもしれません。
どちらの場合もトラックの使い方としてはムダがあり、軽い貨物と重たい貨物をミックスして、容積としても重量としても１００％に近い積載率にするのが効率の良い積み方です。
そのために軽い貨物と重たい貨物をどの比率でミックスすれば良いかについては、下記の記事で解説しています。
簡単な一次方程式で解く方法ですが、この方法は貨物の種類が２種類の場合にしか使えません。
「経験と勘」でも何とかなります。
ところが、これが４種類くらいになってくるとそうはいきません。
無数にある組み合わせから最適解を見つけるのは「経験と勘」では困難です。
今回は「経験と勘」では絶対に無理な、貨物の種類がいくつあっても解ける方法を紹介します。
そのために、OR（オペレーションズリサーチ）の一手法である線形計画法を使います。
前回の記事で線形計画法を掃き出し法で解く方法を紹介しました。
今回は変数が多いため、より効率良く解くための手法であるシンプレックス法を使って解いていきます。
業務改善／ロジスティクスエンジニアリングの求人数が多いのはここだ！
積載率最大化問題を線形計画法に定式化する
次のような状況を想定します。
――――――――――――――――――――――――
４種類の貨物があります。
貨物A
１箱の重量：12kg
１箱の容積：0.05m3
貨物B
１箱の重量：10kg
１箱の容積：0.01m3
貨物C
１箱の重量：8kg
１箱の容積：0.08m3
貨物D
１箱の重量：13kg
１箱の容積：0.02m3
これらの貨物を最大積載容積50m3、最大積載重量25tの大型トラックで運ぶ時、容積ベースでも重量ベースでも積載率を最大にするためには、どの貨物をどれだけ積めば良いでしょうか？
但し、諸事情から次の条件を満たさなくてはいけません。
- 貨物Aと貨物Bの合計重量は8t以内
- 貨物Cと貨物Dの合計容積は25m3以内
――――――――――――――――――――――――
たった４種類の貨物しかありませんが、この問題をすぐに解ける人はまずいないでしょう。
しかし線形計画法を使えば、４０種類でも４００種類でも解けます。
それではまず線形計画法の問題に定式化していきます。
貨物A～Dの積載箱数をそれぞれx、y、z、wとすると、最大積載重量（２５ｔ）と最大積載容積（５０m3）の制約から次の不等式が成り立ちます。
12x + 10y + 8z + 13w ≦ 25000 ・・・①
0.05x + 0.01y + 0.08z + 0.02w ≦ 50 ・・・②
また、おまけの２つの制約条件は次のような不等式で表せます。
12x + 10y ≦ 8000 ・・・③
0.08z + 0.02w ≦ 25 ・・・④
そして重量ベース、容積ベースそれぞれの積載率は、次のように計算できます。
重量ベースの積載率＝ (12x + 10y + 8z + 13w) / 25000
容積ベースの積載率＝ (0.05x + 0.01y + 0.08z + 0.02w) / 50
そして今回はこれら２つの積載率を同時に最大化したいので、２つの平均を目的関数とします。
目的関数＝{(12x + 10y + 8z + 13w) / 25000 + (0.05x + 0.01y + 0.08z + 0.02w) / 50} / 2 ・・・⑤
これで線形計画法の問題に定式化できました。
連立方程式にする
次に、この問題を連立方程式の形にします。
そのために、非負の変数s1～s4を使って式①～④を次のように書き替えます。（s1からs4をスラック変数といいます）
12x + 10y + 8z + 13w + s1 = 25000 ・・・①’
0.05x + 0.01y + 0.08z + 0.02w + s2 = 50 ・・・②’
12x + 10y +s3 = 8000 ・・・③’
0.08z + 0.02w +s4 = 25 ・・・④’
また目的関数の値をmと置くと、式⑤は次のような式で表すことができます。
m = {(12x + 10y + 8z + 13w) / 25000 + (0.05x + 0.01y + 0.08z + 0.02w) / 50} / 2
⇒
m – (12/25000+0.05/50)/2*x – (10/25000+0.01/50)/2*y – (8/25000+0.08/50)/2*z – (13/25000+0.02/50)/2*w = 0 ・・・⑤’
これで式①’～⑤’による連立方程式ができました。
ここでの注意点は、
「変数が３つあるなら式は最低でも３つないと解は定まらない」
と中学で習ったことを覚えている人もいると思いますが、変数の数より式の数の方が多くないと解が一意に定まらないということです。
今回の問題では、変数が９個あるのに式は５式しかないので、解は一意に決まりません。
無数にある解の中から、目的関数の値が最大になるような変数を見つけるところに難しさがあります。
拡大係数行列にする
このような難しさはありますが、普通の連立方程式を解くように掃き出し法が使えますので、まずは５つの式から拡大係数行列を作成します。
定数項を右辺に集めて、変数ごとに列を揃えて書くと次のようになります。
ここでa1～a4はそれぞれ式⑤’中のx、y、z、wの係数です。
制約条件で囲まれる図形に沿って最適解を探すのがシンプレックス法
後は、この行列を単位行列が部分的にできるように行基本変形していくのですが、よく見ると下記のように既に単位行列ができています。
（列を並び替えていますが、行基本変形と同じように列基本変形をしても連立方程式の解は変わりません）
従って、
S1 = 25000
S2 = 5000
S3 = 8000
S4 = 2500
m = 0
は、この連立方程式の無数にある解の中の１つになっています。
そして、残りの変数x、y、z、wを０にすると、制約条件で囲まれる図形の頂点の１つの座標になります。（基底解といいます）
つまり、(s1, s2, s3, s4, m, x, y, z, w) = (25000, 5000, 8000, 2500, 0, 0, 0, 0, 0)の点は、制約条件で囲まれる図形の１つの頂点になっているということです。
では頂点は何個あるのかというと、この問題では１２６個あります。
なぜでしょうか？
次の拡大係数行列を見て下さい。
単位行列になっている部分は、列ごとに見ると、各変数が単位ベクトルになっています。
このように単位ベクトルになっている変数のことを基底変数というのですが、９個の変数のうちどの５個を基底変数にするかの組み合わせの数は9C5通り（＝１２６）あります。
組み合わせの数の計算方法については、こちらを参照して下さい。
ここで基底変数の解は、拡大係数行列の最右列に計算される値になります。
一方、非基底変数はこの行列では求まらないのですが、これらをすべてゼロとすると制約条件で囲まれる図形の頂点になります。
従って、この問題では合計１２６個の頂点があり、その中からm（目的関数の値）が最大になる頂点を探すという最適化問題になります。
ですので、拡大係数行列を何百回も行基本変形して１２６通りの解を求めてもいいのですが、それでは余りにも時間がかかりすぎます。
そこでムダな解はスキップして、意味のある解だけを選択的に調べるためのアルゴリズムがシンプレックス法です。
このアルゴリズムは１９４７年に開発され、現在でも様々なシーンで使われている優れものです。
このアルゴリズムは簡単にいうと、掃き出し法を行う順番を決めるためのものです。
掃き出し法は、解を求めたい変数の列が単位ベクトルになるように行基本変形することでした。
単位ベクトルにするということは、その列のいずれかの成分を１にして、残りはすべて０にすることです。
シンプレックス法は、どの列（つまりどの変数）のどの行（つまりその列のどの成分）を１にすべきかということを順番に教えてくれます。
シンプレックス法をExcelで実行する
以下で具体的に手順を見てみましょう。
まず最初に最下行を見ます。
この中で一番小さい成分を探します。（ステップ①）
すると-0.00096です。
この成分は４列目にあります。
次に４列目の成分で最右列の成分を割ります。
例えば１行４列の成分は８ですが、これで１行１０列の成分２５０００を割るのです。
すると３１２５になります。
これを４列目のすべての成分について行うと、３行目はゼロ割りになってしまい解がありませんが、残りは３１２５、６２５、３１２.５と求まります。
この中で最小値は３１２.５で、それは４行目にあります。（ステップ②）
ステップ①で４列目が、ステップ②で４行目が選ばれました。
この場合、４行４列の成分である０.０８を１になるように行基本変形しなさいというのがシンプレックス法が教えてくれる順番です。
実際にこの列の掃き出しを終了すると、次のようになります。
そうしたら今やった①から③までのステップを繰り返します。
その結果、今度は２列目を掃き出し処理することになります。
これを最下行の成分がすべて非負になるまで繰り返します。
今回の例では４回繰り返したところで下記のように終了しました。
答えの見方は次の通りです。
s1=750、x=447.4、y=263.2、w=1250の時、目的関数mは最大値0.985になるということです。
この中でs1は不等式を等式で表すために導入したスラック変数なので無視すると、商品Aを４４７.４箱、商品Bを２６３.２箱、商品Dを１２５０箱積めば、重量ベースと容積ベースの積載率の平均値が９８.５％となり、これが最大積載率になります。
２つの積載率を分けて計算すると重量ベースで９７％、容積ベースで１００％ですので、ほぼ満載できていることになります。
- 貨物Aと貨物Bの合計重量は8t以内
- 貨物Cと貨物Dの合計容積は25m3以内
という面倒な制約条件付きでここまでの積載率になるようにするのは「経験と勘」では難しいはずで、線形計画法の有用性が分かります。
またシンプレックス法を使うことで、少ない回数の行基本変形で最適解に辿り着くことができました。
１２６通りの組み合わせの中からピンポイントに効率良く計算した結果です。
このように、まともに解いていたら計算量が膨大になってしまう線形計画法を、少ない計算量で効率良くためのアルゴリズムがシンプレックス法です。
ちなみに、このシンプレックス法はExcelのソルバーでも解くことができます。
これを使えば、今まで述べてきたような知識がなくても線形計画法を解くことができます。
Udemyの関連講座
Optimization with Excel: Operations Research without Coding
Optimization with Gurobi, CBC, IPOPT. Linear programming, nonlinear, genetic algorithm. Using Excel, without coding
Assignment and Transportation Problem Operations Research 01
The Hungarian Assignment Problem and Transportation Problem
Optimization with Python: Solve Operations Research Problems
Solve optimization problems with CPLEX, Gurobi, Pyomo… using linear programming, nonlinear, evolutionary algorithms…