GLPK
- GNU Linear Programming Kitの略
- 線形計画問題、整数計画問題、混合整数計画問題のソルバー
- この種の最適化問題のソルバーとしては商用ソフトのCPLEXなどが高性能でよく使われている[1]
- GLPKはCPLEXなどと比べると性能はかなり落ちるが、無料という利点がある
インストール
Ubuntu:
sudo apt-get install glpk-utils
CentOS:
sudo yum install glpk-utils
さっそく使ってみる。
解きたい問題
カードを合成することでカードを強化していくタイプのゲームをプレイしていた。
システムは次のようなものであった。
- 各カードはそれぞれ経験値を持っており、経験値が多いカードほど強い。
- 経験値は合成によってのみ増やすことができる。カードAにカードBを合成することでカードBは消え、カードAの経験値にカードBの経験値が加算される。
- 経験値は一定の値に達するとカンストする。
今手持ちのカードが13枚あり、それぞれ次のような経験値を持っているとする。
- 223066
- 301999
- 225993
- 249157
- 348192
- 317660
- 282072
- 449708
- 87473
- 278373
- 310193
- 207917
- 170856
そして、経験値は1000000でカンストするとする。
最初の3枚を限界まで強化したいが、カンスト分を上回って無駄になる経験値ができるだけ出ないようにしたい。
今少なくとも手持ちの経験値で最初の3枚を限界まで強化できることは既に分かっているとすると、どう合成すれば無駄が無いだろうか。
問題の定式化
最初の3枚を合成元カード1, 2, 3と呼ぶことにする。
他の10枚を素材カード1, 2, …, 10と呼ぶことにする。
合成元カード1, 2, 3の経験値をそれぞれ$b_1, b_2, b_3$とおく。
素材カード1, 2, …, 10の経験値をそれぞれ$e_1, e_2, \ldots, e_{10}$とおく。
合成元カード$i$に素材カード$j$を合成するかどうかを$x_{i, j}$で表す。
$x_{i, j}=1$なら合成し、$x_{i, j}=0$なら合成しない。
経験値のカンスト値を$c$とおく。
最小化する目的関数は次の式である。
\sum_{i=1}^3 \left(b_i + \sum_{j=1}^{10} e_j x_{i, j} - c\right)
制約条件として以下のものがある。
\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\}.\, b_i + \sum_{j=1}^{10} e_j x_{i, j} \geq c\\
&\forall j \in \{n \in \mathbb{N}|1 \leq n \leq 10\}.\, \sum_{i=1}^3 x_{i, j} \leq 1 \\
\end{align}
最後の制約条件は、素材カードはどれか一枚の合成元カードにしか合成できないことを示している。
問題をGLPKの設定ファイルとして書く
GLPKで問題を解く際は2つのファイルを用意する。
一つはモデルファイルと呼ばれるもので、もう一つはデータファイルと呼ばれるものである。
データファイルには$b_i$や$e_j$や$c$のような定数を書き、モデルファイルにはそれ以外の、定数に依存しない問題の定義を書く。
このようにすることで、定数部分のみが異なる同種の問題に対して、データファイルを書き換えるだけで問題を解くことができるようになっているのである。
モデルファイル
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;
var x{I, J} binary;
minimize WASTED_EXP: sum{i in I}(b[i] + sum{j in J}(e[j] * x[i, j]) - c);
s.t. COUNTER_STOP{i in I}: b[i] + sum{j in J}(e[j] * x[i, j]) >= c;
s.t. UNIQUE{j in J}: sum{i in I}(x[i, j]) <= 1;
end;
paramはデータファイルで定数として与える値である。
配列の値をデータファイルから与える場合は上記のように、setでインデックス集合を作るような工夫が必要になる。
varは最適化問題で最適化したいパラメータである。
binaryを後ろに付けることでその値が0か1のどちらかをとるという制約を表せる。
minimizeは最小化する目的関数を表し、s.t.は制約条件を表す。
どちらもWASTED_EXPやCOUNTER_STOPのような適当な名前を付けることができる。
データファイル
efficient_fusion.dat
param i_max := 3;
param b :=
1 223066
2 301999
3 225993
;
param j_max := 10;
param e :=
1 249157
2 348192
3 317660
4 282072
5 449708
6 87473
7 278373
8 310193
9 207917
10 170856
;
param c := 1000000;
end;
こちらは見たまんまなので多分解説しなくても大丈夫だと思う。
GLPKで問題を解く
このようなコマンドを打つ。
glpsol -m efficient_fusion.mod -d efficient_fusion.dat -o efficient_fusion.out
すると標準出力に途中経過などの情報が出力される。
GLPSOL: GLPK LP/MIP Solver 4.40
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...
17 lines were read
Reading data section from efficient_fusion.dat...
21 lines were read
Generating WASTED_EXP...
Generating COUNTER_STOP...
Generating UNIQUE...
Model has been successfully generated
glp_mpl_build_prob: row WASTED_EXP; constant term -2248942 ignored
ipp_basic_tech: 1 row(s) and 0 column(s) removed
ipp_reduce_bnds: 1 pass(es) made, 0 bound(s) reduced
ipp_basic_tech: 0 row(s) and 0 column(s) removed
ipp_reduce_coef: 1 pass(es) made, 0 coefficient(s) reduced
glp_intopt: presolved MIP has 13 rows, 30 columns, 60 non-zeros
glp_intopt: 30 integer columns, all of which are binary
Scaling...
A: min|aij| = 1.000e+00 max|aij| = 4.497e+05 ratio = 4.497e+05
GM: min|aij| = 9.501e-01 max|aij| = 1.052e+00 ratio = 1.108e+00
EQ: min|aij| = 9.027e-01 max|aij| = 1.000e+00 ratio = 1.108e+00
2N: min|aij| = 5.000e-01 max|aij| = 1.000e+00 ratio = 2.000e+00
Constructing initial basis...
Size of triangular part = 13
Solving LP relaxation...
0: obj = -2.248942000e+06 infeas = 8.579e+00 (0)
* 15: obj = 0.000000000e+00 infeas = 0.000e+00 (0)
OPTIMAL SOLUTION FOUND
Integer optimization begins...
+ 15: mip = not found yet >= -inf (1; 0)
+ 95: >>>>> 1.705870000e+05 >= 0.000000000e+00 100.0% (23; 1)
+ 127: >>>>> 1.349990000e+05 >= 0.000000000e+00 100.0% (29; 5)
+ 307: >>>>> 1.160290000e+05 >= 0.000000000e+00 100.0% (65; 20)
+ 337: >>>>> 1.044670000e+05 >= 0.000000000e+00 100.0% (65; 30)
+ 2060: >>>>> 8.681300000e+04 >= 0.000000000e+00 100.0% (230; 315)
+ 2486: >>>>> 8.311400000e+04 >= 0.000000000e+00 100.0% (229; 470)
+ 5092: >>>>> 4.752600000e+04 >= 8.640000000e+02 98.2% (330; 1127)
+ 6239: mip = 4.752600000e+04 >= tree is empty 0.0% (0; 2951)
INTEGER OPTIMAL SOLUTION FOUND
Time used: 0.2 secs
Memory used: 0.5 Mb (540895 bytes)
Writing MIP solution to `efficient_fusion.out'...
結果はefficient_fusion.outに出力されている。
efficient_fusion.out
Problem: efficient_fusion
Rows: 14
Columns: 30 (30 integer, 30 binary)
Non-zeros: 90
Status: INTEGER OPTIMAL
Objective: WASTED_EXP = 47526 (MINimum)
No. Row name Activity Lower bound Upper bound
------ ------------ ------------- ------------- -------------
1 WASTED_EXP 2.29647e+06
2 COUNTER_STOP[1]
800182 776934
3 COUNTER_STOP[2]
698865 698001
4 COUNTER_STOP[3]
797421 774007
5 UNIQUE[1] 1 1
6 UNIQUE[2] 1 1
7 UNIQUE[3] 0 1
8 UNIQUE[4] 1 1
9 UNIQUE[5] 1 1
10 UNIQUE[6] 0 1
11 UNIQUE[7] 1 1
12 UNIQUE[8] 1 1
13 UNIQUE[9] 1 1
14 UNIQUE[10] 1 1
No. Column name Activity Lower bound Upper bound
------ ------------ ------------- ------------- -------------
1 x[1,1] * 0 0 1
2 x[1,2] * 0 0 1
3 x[1,3] * 0 0 1
4 x[1,4] * 1 0 1
5 x[1,5] * 0 0 1
6 x[1,6] * 0 0 1
7 x[1,7] * 0 0 1
8 x[1,8] * 1 0 1
9 x[1,9] * 1 0 1
10 x[1,10] * 0 0 1
11 x[2,1] * 1 0 1
12 x[2,2] * 0 0 1
13 x[2,3] * 0 0 1
14 x[2,4] * 0 0 1
15 x[2,5] * 1 0 1
16 x[2,6] * 0 0 1
17 x[2,7] * 0 0 1
18 x[2,8] * 0 0 1
19 x[2,9] * 0 0 1
20 x[2,10] * 0 0 1
21 x[3,1] * 0 0 1
22 x[3,2] * 1 0 1
23 x[3,3] * 0 0 1
24 x[3,4] * 0 0 1
25 x[3,5] * 0 0 1
26 x[3,6] * 0 0 1
27 x[3,7] * 1 0 1
28 x[3,8] * 0 0 1
29 x[3,9] * 0 0 1
30 x[3,10] * 1 0 1
Integer feasibility conditions:
KKT.PE: max.abs.err = 0.00e+00 on row 0
max.rel.err = 0.00e+00 on row 0
High quality
KKT.PB: max.abs.err = 0.00e+00 on row 0
max.rel.err = 0.00e+00 on row 0
High quality
End of output
重要な部分を抜粋。
Status: INTEGER OPTIMAL
これは最適解が見つかったことを示している。
No. Column name Activity Lower bound Upper bound
------ ------------ ------------- ------------- -------------
1 x[1,1] * 0 0 1
2 x[1,2] * 1 0 1
見つかった解はActivityの欄に書かれている。
この場合、$x_{1,1}$は0、$x_{1,2}$は1が見つかった解だということである。
解答
| カード名 | 経験値 |
|---|---|
| 合成元カード1 | 223066 |
| 合成元カード2 | 301999 |
| 合成元カード3 | 225993 |
| 素材カード1 | 249157 |
| 素材カード2 | 348192 |
| 素材カード3 | 317660 |
| 素材カード4 | 282072 |
| 素材カード5 | 449708 |
| 素材カード6 | 87473 |
| 素材カード7 | 278373 |
| 素材カード8 | 310193 |
| 素材カード9 | 207917 |
| 素材カード10 | 170856 |
- 合成元カード1に素材カード4, 8, 9を合成する
- 合成元カード2に素材カード1, 5を合成する
- 合成元カード3に素材カード2, 7, 10を合成する
まとめ
日常のちょっとした問題も線形計画問題や整数計画問題に帰着することがあるので、ソルバーを使えると便利かもしれない。
コメントによる補足
他の方から頂いた補足。
あと、せっかくMathProgで書いているのでしたら、解を表示するところもやると面白いですよ。
モデルファイルのend;の前にこんな感じに書けば表示できるはず・・・（最近書いてないのでうろ覚え）solve; for{j in J}{ printf "%d -> %d\n", j, sum{i in I}i * x[i][j]; }制約式のUNIQUEがあるので掛け算してsumを取ることで選択した値が出てくるのがポイントですね。
試してみる。
efficient_fusion.mod
（途中略）
solve;
for {j in J} {
printf "%d -> %d\n", j, sum{i in I}(i * x[i, j]);
}
end;
（途中略）
INTEGER OPTIMAL SOLUTION FOUND
Time used: 0.2 secs
Memory used: 0.5 Mb (540895 bytes)
1 -> 2
2 -> 3
3 -> 0
4 -> 1
5 -> 2
6 -> 0
7 -> 3
8 -> 1
9 -> 1
10 -> 3
Model has been successfully processed
Writing MIP solution to `efficient_fusion.out'...
情報ありがとうございました。