GLPKのプログラムをMiniSatで動かす

<--目次へ
GLPKでは少し複雑になると最適解までなかなか到達しないことがよくあります。ソルバーの中ではかなり遅めですが、GLPKからMiniSatを使うことで、劇的に速度が上がります。
生産現場で役に立つかもしれない最適化手法1
のプログラムをMiniSatで解く方法を説明します。
GLPKにはMiniSatが含まれていませんので、インストールしておくことも忘れずに。

プログラムの変更

制約式(s.t.)中の変数、パラメータを整数型に

オプション--minisatを付け加えて　glpsolを実行したらわかりますが、制約式に浮動小数点のパラメータが含まれているとエラーになります。
例のプログラムでは、制約式中に浮動小数点パラメータweightとmaxProductWeightが使われていますので、これを整数に変更します。有効ケタ数を確保するためにここでは1000倍したweight_intとmaxProductWeight_intに代入をして使用しています。

目的関数(minimize, maximize)を制約式に

目的関数があっても無視されてエラーが出ませんので注意してください。これを制約式に変えます。元がminimizeだったので、適当な値以下になる制約式にします。ここでは最後のパレット以外は最大数maxQuantityInPalletまで積載した時の値quantitiyLastPalletを計算してそれを上限に使っています。もしそれが満たされない場合には「UNSATISFIABLE」を出力して停止します。その場合は条件を満たすまで、手で1つずつ大きくしていきます
(実際にはallowFreeSpaceを1つずつ大きくしていくことで実現しています)
ちょっと手間がかかりますが、計算速度が劇的に上がりますのであまり気にならないと思います。

変更後のプログラムです。--minisatオプションを加えることでエラーなく実行できます。
ここに示したデータでは
param allowFreeSpace := 1;
でないと条件を満たせませんので、使われる時にはデータファイルを変更してください。

pallet04.mod

# Version 0.4
set Product; # 製品名を読み込む
param maxQuantityInPallet integer; # パレットに載せることができる最大個数を読み込む
param quantity{Product} integer; # 製品の個数を読み込む
param weight{Product}; # 1個あたりの重量を読み込む
param maxProductWeight; # パレットに載せることができる最大重量を読み込む
param nPallet := ceil((sum{p in Product}quantity[p])/maxQuantityInPallet); # パレットの必要数を計算する
set Pallet := 1..nPallet; # パレット番号
# 追加
param allowFreeSpace; # 最後のパレット以外に許容される空き個数を読み込む
param weight_int{p in Product}:= round(weight[p]*1000); # 整数に変換 ここでは1000倍してグラムに
param maxProductWeight_int:= round(maxProductWeight*1000); # 同じく整数に
param quantitiyLastPallet := sum{p in Product}quantity[p]-maxQuantityInPallet*(nPallet-1); # 全部詰まっている時の最後のパレットの製品数
# ここまで

var assignProductInPallet{Product, Pallet} binary; # 製品のパレットへの割り当てを入れるバイナリ変数

# 製品を1つのパレットに割り当てる
s.t. restrictOnePallett{p in Product}: sum{q in Pallet}assignProductInPallet[p,q]==1;
# 最後のパレット以外は載せることができる最大個数以下に制限
s.t. restrictMaxQuantity{q in 1..nPallet-1}: sum{p in Product}assignProductInPallet[p,q]*quantity[p]<=maxQuantityInPallet;
# パレットに載せることができる最大重量以下に（変更あり)
s.t. limitPalletWeight{q in Pallet}: sum{p in Product}assignProductInPallet[p,q]*quantity[p]*weight_int[p]<=maxProductWeight_int;
# 最後のパレットの積載個数を制限 (変更あり)
s.t. restrictLastPattel: sum{p in Product}assignProductInPallet[p, nPallet]*quantity[p]<=quantitiyLastPallet+allowFreeSpace;

solve;

# 出力
for{q in Pallet}{
	# 最大積載数になっていないパレットには印をつける
    printf "Pallet %d %s", q, if sum{p in Product}assignProductInPallet[p,q]*quantity[p]!=maxQuantityInPallet then "******" else "";
    printf "(Quantitiy: %d, Weight: %f kg)\n", sum{p in Product}assignProductInPallet[p,q]*quantity[p], sum{p in Product}assignProductInPallet[p,q]*quantity[p]*weight[p];
    printf{p in Product : assignProductInPallet[p,q]==1 and quantity[p]!=0} "    %s   %d\n", p, quantity[p];
}
# Excelへ画面上でのコピペ用のタブ区切りデータ
printf "\nFor Excel\n";
printf{q in Pallet} "\tPallet %s", q; printf "\n";
for{p in Product}{
    printf "%s", p;
    printf{q in Pallet} "\t%s", if assignProductInPallet[p,q]*quantity[p]!=0 then quantity[p] else "";
    printf "\n";
}
printf "合計";
printf{q in Pallet} "\t%d", sum{p in Product}assignProductInPallet[p,q]*quantity[p];
printf "\n";

データセクション(下記)にはparam allowFreeSpaceが追加されています。

data2a.dat

data;
# パレットに載せることができる最大個数
param maxQuantityInPallet := 48;
# パレットに載せることができる最大重量
param maxProductWeight := 1000;
# 最後のパレット以外に許容される空き個数 通常は0で UNSATISFIABLEで停止すると1ずつ増やす
param allowFreeSpace := 0;
# 製品名 個数 1個あたりの重量kg
param : Product : quantity weight :=
"製品1"  14  7.632
"製品2"  10  9.648
"製品3"  39  10.08
"製品4"  22  11.28
"製品6"  20  11.28
"製品7"  20  7.56
"製品8"  8  7.224
"製品9"  8  5.6112
"製品10"  36  10.815
"製品12"  6  7.344
"製品13"  4  6.384
"製品14"  0  10.332
"製品15"  4  10.332
"製品16"  10  6.8832
"製品17"  10  12.4416
"製品18"  12  6.48
"製品19"  16  10.896
"製品21"  19  7.6293
;
end;