GLPKでより定式化が難しい問題を解いてみた

前回は定式化が簡単な問題を扱ったが、今回はより定式化が難しい問題を解くことで、GLPK（というより混合整数計画問題）の応用可能性の高さを示したい。
前回の内容を前提にしているので、前回をまだ読んでいない方はそちらからどうぞ。

解きたい問題

前回のコメントにあったこの問題を解く。

３枚中２枚はカンストできるときにできるだけ合成しつつカンスト余剰分を減らしたい・・・とか考えたくなっちゃいますね。

3枚全部はカンストできないという状況を作るために素材カードの枚数は減らしておく。

カード名	経験値
合成元カード1	223066
合成元カード2	301999
合成元カード3	225993
素材カード1	249157
素材カード2	348192
素材カード3	317660
素材カード4	282072
素材カード5	449708
素材カード6	87473
素材カード7	278373

問題の定式化

問題は2段階の構造を持っている。

1. カンストさせるカードの枚数は多い方がいい
2. カンストさせるカードの枚数が同じなら、余剰経験値は少ない方がいい

前回は2.の問題だけ考えていたようなものであるが、今回は1.も考えなければならない。
そこで、1.の問題に対応する変数として、新たに0-1変数$y_i$を導入することにする。
合成元カード$i$をカンストさせる場合は$y_i=0$、カンストさせない場合は$y_i=1$とする。（直感的な定義とは逆だがこうした方が後で式が簡単になる）

また、あるトリックを成立させるために、何でもいいので巨大な値を持つ定数として$M_1, M_2$を導入することにする。

最小化する目的関数は次の式である。

M_1 \sum_{i=1}^3 y_i + \sum_{i=1}^3 \max\left(b_i + \sum_{j=1}^{10} e_j x_{i, j} - c,\,0\right)

$M_1$を使うことで問題の優先順位を表現している。
なお、カンストさせない場合余剰経験値は0になるということを表現するためにmax関数が現れているが、これは後で消すことができるので今はこのままにしておく。

制約条件として以下のものがある。

\begin{align}
&\forall i \in \{n \in \mathbb{N}|1 \leq n \leq 3\}.\,\forall j \in \{n \in \mathbb{N}|1 \leq n \leq 10\}.\,x_{i, j} \in \{0, 1\} \\
&\forall i \in \{n \in \mathbb{N}|1 \leq n \leq 3\}.\,y_{i} \in \{0, 1\} \\
&\forall j \in \{n \in \mathbb{N}|1 \leq n \leq 10\}.\, \sum_{i=1}^3 x_{i, j} \leq 1 \\
&\forall i \in \{n \in \mathbb{N}|1 \leq n \leq 3\}.\, b_i + \sum_{j=1}^{10} e_j x_{i, j} - c \geq -M_2 y_i\\
\end{align}

最後の制約条件を説明すると、合成元カード$i$をカンストさせる場合、$y_i=0$にできるが、カンストさせられない場合、制約条件を満たすために$y_i=1$にせざるを得ないというトリックになっている。

このように巨大な定数を使う方法はbig-M法と呼ばれ、多くの論理関係を定式化できるという利点もあるが、数値計算が不安定になりやすいという欠点もある。
そのため、巨大な定数はなるべく使わない方が良く、使う場合でも制約式の意味を変えない範囲でぎりぎりまで低い値を使うのが良いとされる。（参考: ここまで解ける整数計画）

max関数を消すためさらにもう一捻りする。
連続変数$z_i$を導入する。

最小化する目的関数は次の式である。

M_1 \sum_{i=1}^3 y_i + \sum_{i=1}^3 z_i

制約条件として以下のものがある。

\begin{align}
&\forall i \in \{n \in \mathbb{N}|1 \leq n \leq 3\}.\,\forall j \in \{n \in \mathbb{N}|1 \leq n \leq 10\}.\,x_{i, j} \in \{0, 1\} \\
&\forall i \in \{n \in \mathbb{N}|1 \leq n \leq 3\}.\,y_{i} \in \{0, 1\} \\
&\forall j \in \{n \in \mathbb{N}|1 \leq n \leq 10\}.\, \sum_{i=1}^3 x_{i, j} \leq 1 \\
&\forall i \in \{n \in \mathbb{N}|1 \leq n \leq 3\}.\, b_i + \sum_{j=1}^{10} e_j x_{i, j} - c \geq -M_2 y_i\\
&\forall i \in \{n \in \mathbb{N}|1 \leq n \leq 3\}.\, b_i + \sum_{j=1}^{10} e_j x_{i, j} - c \leq z_i\\
&\forall i \in \{n \in \mathbb{N}|1 \leq n \leq 3\}.\, 0 \leq z_i\\
\end{align}

最小化する目的関数にmax関数が現れている場合は、このように各オペランドに対応した制約条件を追加することでmax関数を消すことができる。

問題をGLPKの設定ファイルとして書く

efficient_fusion.mod

param i_max;
set I := 1..i_max;
param b{I};

param j_max;
set J := 1..j_max;
param e{J};

param c;

param M1;
param M2;

var x{I, J} binary;
var y{I} binary;
var z{I} >= 0;

minimize OBJECTIVE_FUNC: M1 * sum{i in I}(y[i]) + sum{i in I}(z[i]);
s.t. UNIQUE{j in J}: sum{i in I}(x[i, j]) <= 1;
s.t. COUNTER_STOP_IF_POSSIBLE{i in I}: b[i] + sum{j in J}(e[j] * x[i, j]) - c >= -M2 * y[i];
s.t. MAX{i in I}: b[i] + sum{j in J}(e[j] * x[i, j]) - c <= z[i];

solve;

for {j in J} {
    printf "%d -> %d\n", j, sum{i in I}(i * x[i, j]) * (1 - y[i]);
}

end;

efficient_fusion.dat

param i_max := 3;
param b :=
1 223066
2 301999
3 225993
;
param j_max := 7;
param e :=
1 249157
2 348192
3 317660
4 282072
5 449708
6 87473
7 278373
;
param c := 1000000;
param M1 := 10000000;
param M2 := 10000000;
end;

GLPKで問題を解く

glpsol -m efficient_fusion.mod -d efficient_fusion.dat -o efficient_fusion.out

GLPSOL: GLPK LP/MIP Solver, v4.52
Parameter(s) specified in the command line:
 -m efficient_fusion.mod -d efficient_fusion.dat -o efficient_fusion.out
Reading model section from efficient_fusion.mod...
29 lines were read
Reading data section from efficient_fusion.dat...
20 lines were read
Generating OBJECTIVE_FUNC...
Generating UNIQUE...
Generating COUNTER_STOP_IF_POSSIBLE...
Generating MAX...
Model has been successfully generated
GLPK Integer Optimizer, v4.52
14 rows, 27 columns, 75 non-zeros
24 integer variables, all of which are binary
Preprocessing...
3 constraint coefficient(s) were reduced
13 rows, 27 columns, 69 non-zeros
24 integer variables, all of which are binary
Scaling...
 A: min|aij| =  1.000e+00  max|aij| =  7.769e+05  ratio =  7.769e+05
GM: min|aij| =  9.565e-01  max|aij| =  1.045e+00  ratio =  1.093e+00
EQ: min|aij| =  9.565e-01  max|aij| =  1.000e+00  ratio =  1.045e+00
2N: min|aij| =  5.000e-01  max|aij| =  1.076e+00  ratio =  2.152e+00
Constructing initial basis...
Size of triangular part is 13
Solving LP relaxation...
GLPK Simplex Optimizer, v4.52
13 rows, 27 columns, 69 non-zeros
      0: obj =   0.000000000e+00  infeas =  1.373e+02 (0)
*    13: obj =   1.183598591e+07  infeas =  0.000e+00 (0)
*    21: obj =   3.041532485e+06  infeas =  0.000e+00 (0)
OPTIMAL LP SOLUTION FOUND
Integer optimization begins...
+    21: mip =     not found yet >=              -inf        (1; 0)
+    39: >>>>>   1.017962200e+07 >=   1.000000000e+07   1.8% (8; 0)
+    52: >>>>>   1.009214900e+07 >=   1.000000000e+07   0.9% (9; 2)
+   137: >>>>>   1.008799200e+07 >=   1.000000000e+07   0.9% (16; 35)
+   226: >>>>>   1.005656100e+07 >=   1.000000000e+07   0.6% (14; 91)
+   242: mip =   1.005656100e+07 >=     tree is empty   0.0% (0; 155)
INTEGER OPTIMAL SOLUTION FOUND
Time used:   0.0 secs
Memory used: 0.2 Mb (181060 bytes)
1 -> 3
2 -> 1
3 -> 0
4 -> 3
5 -> 1
6 -> 0
7 -> 3
Model has been successfully processed
Writing MIP solution to `efficient_fusion.out'...

解答

• 合成元カード1に素材カード2, 5を合成する
• 合成元カード3に素材カード1, 4, 7を合成する

まとめ

GLPKで解ける問題は意外と多い。