生産現場で役に立つかもしれない最適化手法1【GLPK】

[<--目次へ](https://qiita.com/ki073/items/508b3c20eb03078703dd)
## はじめに
生産現場などでありそうな課題を取り上げて、整数計画法を使って考えてみました。GLPKのインストールや使い方は上のリンクを参考にしてください。

## 最適化の条件
* 製造した製品をパレットに入れて出荷したいが、1つのパレットにはupperLimitQuantity個(ここでは48)の製品を載せることができる。可能な限り最大数まで製品をパレットに載せたい。ただしパレット重量の上限を超えない範囲とする。
* 同じ製品はできるだけ同じパレットに載せたい


## プログラム

GLPKで動くプログラムです。

1. パレット積載数の上限upperLimitQuantityでdefautの値を設定
<details><summary>個々のパレットでdefautとは違う値を設定するときには`maxQuantityInPallet[3]=20;`のようにパレット番号を指定して変更してください（変数名に注意）。
ただし、積載数量が Pallet1 >= Pallet2 >= Pallet3 ...のようにパレット番号が1つ前より大きくならないように制約条件を設定していますので注意してください。</summary></details>
1. パレット積載数の下限を設定
<details><summary>`minQuantityInPallet[パレット番号]`でパレットごとの最低積載数を設定できます。defaultは0となっています。積載数を固定するには、最大積載数と最低積載数を同じ値にすることでできます。</summary></details>
1. パレット積載重量の上限の設定
<details><summary>データブロック中に必ず書き込んでください.
例:`param maxPalletWeight default 1000;`
個々のパレットで別々の値を設定することもできます。方法は「データセクションの書き方」を参考にしてください。製品の重量を考慮しない場合は、weightを0にすると多少は速くなるかもしれません.(後述)</summary></details>
1. 使用パレット数
<details><summary>製品数から最低限必要なパレット数を求めます。それ以上のパレットを使う場合には`minNPallets`を設定するとで増やすことができます。</summary></details>
1. 積載する製品情報
<details><summary>下記のプログラムを参考にしてください。data2.txtは重量を考慮しない場合ですので参考にしてください。</summary></details>

```GLPK:pallet04.mod
# Version 0.4.2
set Product; # 製品名を読み込む
param weight{Product} default 0; # 1個あたりの重量を読み込む
param quantity{Product} integer; # 製品の個数を読み込む
param totalProductQuantity := sum{pd in Product}quantity[pd]; # 製品総数を計算
param upperLimitQuantity integer; # パレットの最大積載数を読み込む
param requireNPallets := ceil((totalProductQuantity)/upperLimitQuantity); # パレット必要数を計算
param minNPallets integer, default 1; # 上のパレット必要数より大きいとそのパレット数を使う。通常は1でよい。

printf "*********************\n";
printf "%d個の製品があり パレット当たり %d 個積載すると", totalProductQuantity, upperLimitQuantity;
printf "最低 %d パレット必要です。\n", requireNPallets;
printf "内訳  %d 個積載できるパレットが %d ", upperLimitQuantity, totalProductQuantity div upperLimitQuantity;
printf (if totalProductQuantity mod upperLimitQuantity>0 then "と %d 個積載できるパレットが 1 " else "") & "必要です\n", totalProductQuantity mod upperLimitQuantity;
printf "*********************\n";

param nPallets := max(requireNPallets, minNPallets); # 使用パレット数
set Pallet := 1..nPallets; # パレット番号をつける
param maxQuantityInPallet{Pallet} default upperLimitQuantity; # 各パレットに載せることができる最大個数
param maxPalletWeight{Pallet}; # パレットに載せることができる最大重量を読み込む
param minQuantityInPallet{Pallet} integer default 0; # 各パレットの最低積載数

# パレット別の製品積載数量
var nProductInPallet{pd in Product, Pallet} integer >=0, <=min(upperLimitQuantity, quantity[pd]);
# パレット別の製品が積載されているかを示すバイナリ変数
var containProductInPallet{Product, Pallet} binary;

# 製品の総数を合わす
s.t. constraintTotalQuantity{pd in Product}: sum{pt in Pallet}nProductInPallet[pd, pt]=quantity[pd];
# パレットの最大積載数以下に
s.t. constrainPalletMaxQuantity{pt in Pallet}: sum{pd in Product}nProductInPallet[pd, pt]<=maxQuantityInPallet[pt];
# 積載数はパレット番号が大きくなると少なく
s.t. constrainQuantity{pt in 1..nPallets-1}: sum{pd in Product}nProductInPallet[pd,pt]>=sum{pd in Product}nProductInPallet[pd,pt+1];
# パレットに製品が含まれるか
s.t. containPallet1{pd in Product, pt in Pallet}: containProductInPallet[pd,pt]<=nProductInPallet[pd,pt];
s.t. containPallet2{pd in Product, pt in Pallet}: containProductInPallet[pd,pt]*maxQuantityInPallet[pt]>=nProductInPallet[pd,pt];
# パレット総重量が限度を超えないように
s.t. constrainPalletMaxWeight{pt in Pallet}: sum{pd in Product: weight[pd]>0}weight[pd]*nProductInPallet[pd,pt]<=maxPalletWeight[pt];
# パレットの最低積載数以上に
s.t. constrainPalletMinQuantity{pt in Pallet:  minQuantityInPallet[pt]>0}: sum{pd in Product}nProductInPallet[pd, pt]>=minQuantityInPallet[pt];

minimize minimizeSumContain: sum{pd in Product, pt in Pallet}containProductInPallet[pd,pt];

solve;

# 出力
printf "\nパレット別積載製品\n";
for{pt in Pallet}{
    printf "Pallet %d  積載数 %d", pt, sum{pd in Product}nProductInPallet[pd,pt];
    printf (if sum{pd in Product}weight[pd]*nProductInPallet[pd,pt]>0 then "    重量 %f <= %f\n" else "\n"), sum{pd in Product}weight[pd]*nProductInPallet[pd,pt], maxPalletWeight[pt];
    printf{pd in Product : nProductInPallet[pd,pt]>0} "    %s   %d\n", pd,  nProductInPallet[pd,pt];
}
printf "\n製品別積載パレット\n";
for{pd in Product}{
    printf "%s  ", pd;
    printf{pt in Pallet : nProductInPallet[pd,pt]>0} "   Pallet %d: %d", pt  ,  nProductInPallet[pd,pt];
    printf "\n";
}

# Excelへ画面上でのコピペ用のタブ区切りデータ
printf "\nFor Excel\n";
printf{pt in Pallet} "\tPallet %s", pt; printf "\n";
for{pd in Product}{
    printf "%s", pd;
    printf{pt in Pallet} (if nProductInPallet[pd, pt]>0 then "\t%s" else "\t"),  nProductInPallet[pd, pt];
    printf "\n";
}
printf "合計";
printf{pt in Pallet} "\t%d", sum{pd in Product}nProductInPallet[pd, pt];
printf "\n";
```

(2021.12.27) パレットの最低積載数が反映されていなかったので修正しています。

データ部分を分けていますが、上のプログラム部分と繋いで1つのファイルにして実行してください。

```GLPK:data1.txt
data;
# パレットに載せることができる最大個数
param upperLimitQuantity := 48;
# パレットに載せることができる最大重量
param maxPalletWeight default 520;
# パレットの最大積載数を変更
param maxQuantityInPallet[8] := 24;
# 製品名 個数 1個あたりの重量
param : Product : quantity weight :=
"製品1"   14  7.632
"製品2"    5  9.648
"製品3"   58  10.08
"製品4"   19  11.28
"製品5"   22  9.24
"製品6"   16  11.28
"製品7"   16  7.56
"製品8"   17  7.224
"製品9"    8  5.6112
"製品10"  38  10.815
"製品11"   1  7.287
"製品12"  10  7.344
"製品13"   4  6.384
"製品14"   0  10.332
"製品15"  10  10.332
"製品16"  20  6.8832
"製品17"  32  12.4416
"製品18"  26  6.48
"製品19"  10  10.896
"製品20"   1  8.256
"製品21"  32  7.6293
;
end;
```

重量を考慮しない場合です。上と比べて`param : Product : quantity:=`
でweightを抜いています。

```GLPK:data2.txt
data;
# パレットに載せることができる最大個数
param upperLimitQuantity := 48;
# パレットに載せることができる最大重量
param maxPalletWeight default 1000;
# パレットの最大積載数を変更
param maxQuantityInPallet[5] := 32;
# 製品名 個数
param : Product : quantity:=
"品番E" 16
"品番B" 7
"品番I" 4
"品番G" 10
"品番Q" 10
"品番A" 1
"品番H" 10
"品番K" 2
"品番S" 30
"品番D" 54
"品番L" 32
"品番N" 14
"品番F" 0
"品番R" 16
"品番P" 10
"品番M" 3
"品番O" 5
;
end;
```

(2021.12.24 追加)

###なかなか結果が出ない場合
　+ 83065: mip =   1.700000000e+01 >=   1.600000000e+01   5.9% (548; 16715)
のような表示がでるがなかなか結果が出力されないことがあります。目的関数(minimizeの行の式）の計算結果が17と16との間になることが分かったということです。最適解の16になる可能性があるが、許容解(17)が得られていますので、オプション `--tmlim 120` の追加で秒数を指定して計算を打ち切り、許容解を出力することができます。

###細かな条件の設定

データブロック(data;より後)に以下の設定をすることで、条件を設定できます。

1. パレット数を増やす
`param minNPallets := 7;`
のようにパレット数を指定することで数を増やすことができます。プログラムを動かすと最初の方に必要パレット数が出力されますのでそれを参考に設定してください。

1. パレットへの最大積載数と最小積載数を設定
`maxQuantityInPallet[6] := 25;`
`minQuantityInPallet[6] := 20;`
のように最大積載数と最小積載数を設定でます。[]の中はパレット番号です。
ただし積載数は、パレット番号n>=パレット番号n+1 のように番号が大きくなるに従って同じか減少するように制約をかけていますので、番号の一番大きいパレットのみに設定するのが良いと思います。

1. 最大積載重量を設定
`param maxPalletWeight[6] :=1000;`
パレットごとに最大積載重量を変更できます。

###パレットに載せる個数を8の倍数に

コメント欄にあるように、パレットに載せる個数を8の倍数にしたいとのことですが、ここまでのプログラムでは考慮されていません。データブロックに
`param maxQuantityInPallet[5] := 32;`
のように、番号の一番大きなパレットに載せる最大個数を8の倍数を指定することで、条件を満たすことがあります。この方法は速度低下もなく実用的なのですが、明示的に8の倍数に制約する方法を以下に示します。
solve;の上に次のプログラムを追加してみてください。
(2021.12.24 変更)

```GLPK:pallet04_ext
set checkDiv8 := 1..nPallets; # 割り切れるかチェックする範囲
var quotient{pt in checkDiv8} integer, >=0, <=maxQuantityInPallet[pt] div 8; # 商
var divisible{checkDiv8} binary; # 8で割り切れると1、余りがあると0になる変数
s.t. div8_1{pt in checkDiv8}: sum{pd in Product}nProductInPallet[pd,pt]>=quotient[pt]*8;
s.t. div8_2{pt in checkDiv8}: sum{pd in Product}nProductInPallet[pd,pt]<=quotient[pt]*8+7*(1-divisible[pt]);
s.t. div8_3: sum{pt in checkDiv8}divisible[pt]>=card(checkDiv8)-1; # 割り切れないパレット数が1つ以下
```

速度低下はそうなさそうなので、全部のパレットを対象に8の倍数かチェックするように変更しました。1つのパレットを除いて8の倍数になるように制約条件追加しています。
速度低下が問題になるようでしたら、次のように変更すると番号の大きい方から2つのパレットだけをチェクするようになります。
`checkDiv8 := nPallets-1..nPallets; # 割り切れるかチェックする範囲`