昨日は数式への落とし込みをやった。
今日はそれを実際にGLPKで解いてみる。
GLPKのダウンロードとインストール
まずはGLPKのダウンロードとインストールから。
これは簡単で、GLPKのサイトからtarボールをダウンロードしてきて、makeするだけ。
詳細は省略。
ちゃんとインストール出来ると、glpsolというコマンドが使えるようになる。
モデルの記述
ソルバーに問題を解かせるためには、まず問題を数式で表現してやらないといけない。
これはモデルの記述と呼ばれて、特定のモデル記述言語を使って行うことが多い。
GLPKの場合、何個かのモデル記述言語に対応していて、AMPL（A Mathematical Programming Language）というモデル記述言語（とほぼ同じモデル記述言語）を使うことが出来る。
書いたモデルが、以下。
#---------------------------------------- # 20150715.mod #---------------------------------------- # 以下の問題を解く。 # # 次の枠内に # 数字がいくつずつあるかを # 括弧内に数字で記入しなさい。 # +--------------------+ # | 2015/07/15 | # | 0 ( ) 5 ( ) | # | 1 ( ) 6 ( ) | # | 2 ( ) 7 ( ) | # | 3 ( ) 8 ( ) | # | 4 ( ) 9 ( ) | # +--------------------+ #---------------------------------------- #---------------------------------------- # モデル #---------------------------------------- #------------------------------ # 集合の定義 #------------------------------ # 対象の数字 set target; # 各桁の数字 set nil; set index := target union nil; # 桁を示す数字 set digit; # 制約C,D用 set l_digit; # 制約D用 set zero; set nonzero := index diff zero; #------------------------------ # 変数の定義 #------------------------------ # その数字が入れば1,そうでなければ0とする var x{target, index, digit} binary; #------------------------------ # パラメータの定義 #------------------------------ # その数字の引数での数 param arg{target}; #------------------------------ # 目的関数 #------------------------------ # ダミー minimize obj: 0; #------------------------------ # 制約式 #------------------------------ # 自己言及の条件を満たす subject to A {k in target}: 1 * (sum {i in target} i * x[k,i,1]) + 10 * (sum {i in target} i * x[k,i,2]) = sum {j in target, n in digit} x[j,k,n] + arg[k]; # 各桁に入る数字は高々1つ subject to B {k in target, n in digit}: sum {i in index} x[k,i,n] = 1; # n桁目に数字が入っていないなら、 # (n+1)桁目にも数字が入っていない subject to C {k in target, n in l_digit}: x[k,-1,n+1] >= x[k,-1,n]; # 数字が入っている桁で一番大きい桁の数字は0でない subject to D {k in target, n in l_digit}: x[k,-1,n+1] <= sum {i in nonzero} x[k,i,n]; # N桁目には数字が入っていない subject to E {k in target}: x[k,-1,3] = 1; # 1桁目には数字が絶対入っている subject to F {k in target}: x[k,-1,1] = 0; #---------------------------------------- # データ #---------------------------------------- data; #------------------------------ # 集合 #------------------------------ # 各桁の数字 set nil := -1; # 制約E用 set zero := 0; # 対象の数字 set target := 0 1 2 3 4 5 6 7 8 9; # 桁を示す数字 set digit := 1 2 3; # 制約B,E用 set l_digit := 1 2; #------------------------------ # パラメータ #------------------------------ # 数字列 param arg := 0 3 1 3 2 2 3 1 4 1 5 3 6 1 7 2 8 1 9 1 ;
詳細は説明しないけど、雰囲気は伝わるんじゃないかな。
glpsolの実行
さて、モデルが書けたら、あとはGLPKのglpsolコマンドを実行して、問題を解かせる。
$ glpsol -m 20150715.mod -o /dev/stdout > 20150715.ans
-mオプションでモデルファイルを指定し、-oオプションで答えの出力先を指定する。
ただ、glpsolは答え以外にもいろんな情報を出力するので、上記のようにそれらをまとめて標準出力にしてしまって、リダイレクトして書き出したファイルを見た方が、情報は見やすいかもしれない。
書き出した内容は、以下のような感じ。
GLPSOL: GLPK LP/MIP Solver, v4.55
Parameter(s) specified in the command line:
-m 20150715.mod -o /dev/stdout
Reading model section from 20150715.mod...
Reading data section from 20150715.mod...
20150715.mod:141: warning: unexpected end of file; missing end statement inserted
141 lines were read
Generating obj...
Generating A...
Generating B...
Generating C...
Generating D...
Generating E...
Generating F...
Model has been successfully generated
GLPK Integer Optimizer, v4.55
101 rows, 330 columns, 1071 non-zeros
330 integer variables, all of which are binary
Preprocessing...
12 hidden covering inequaliti(es) were detected
50 rows, 132 columns, 467 non-zeros
132 integer variables, all of which are binary
Scaling...
A: min|aij| = 1.000e+00 max|aij| = 2.000e+01 ratio = 2.000e+01
GM: min|aij| = 2.659e-01 max|aij| = 3.761e+00 ratio = 1.414e+01
EQ: min|aij| = 7.071e-02 max|aij| = 1.000e+00 ratio = 1.414e+01
2N: min|aij| = 6.250e-02 max|aij| = 1.250e+00 ratio = 2.000e+01
Constructing initial basis...
Size of triangular part is 50
Solving LP relaxation...
GLPK Simplex Optimizer, v4.55
50 rows, 132 columns, 467 non-zeros
0: obj = 0.000000000e+00 infeas = 1.879e+01 (0)
* 27: obj = 0.000000000e+00 infeas = 8.885e-16 (0)
OPTIMAL LP SOLUTION FOUND
Integer optimization begins...
+ 27: mip = not found yet >= -inf (1; 0)
+ 921: >>>>> 0.000000000e+00 >= 0.000000000e+00 0.0% (17; 90)
+ 921: mip = 0.000000000e+00 >= tree is empty 0.0% (0; 161)
INTEGER OPTIMAL SOLUTION FOUND
Time used: 0.0 secs
Memory used: 0.5 Mb (511021 bytes)
Writing MIP solution to '/dev/stdout'...
Problem: 20150715
Rows: 101
Columns: 330 (330 integer, 330 binary)
Non-zeros: 1071
Status: INTEGER OPTIMAL
Objective: obj = 0 (MINimum)
No. Row name Activity Lower bound Upper bound
------ ------------ ------------- ------------- -------------
1 obj 0
2 A[0] 3 3 =
3 A[1] 3 3 =
4 A[2] 2 2 =
5 A[3] 1 1 =
6 A[4] 1 1 =
7 A[5] 3 3 =
8 A[6] 1 1 =
9 A[7] 2 2 =
10 A[8] 1 1 =
11 A[9] 1 1 =
12 B[0,1] 1 1 =
13 B[0,2] 1 1 =
〜省略〜
80 D[9,1] 0 -0
81 D[9,2] 0 -0
82 E[0] 1 1 =
83 E[1] 1 1 =
84 E[2] 1 1 =
85 E[3] 1 1 =
86 E[4] 1 1 =
87 E[5] 1 1 =
88 E[6] 1 1 =
89 E[7] 1 1 =
90 E[8] 1 1 =
91 E[9] 1 1 =
92 F[0] 0 -0 =
93 F[1] 0 -0 =
94 F[2] 0 -0 =
95 F[3] 0 -0 =
96 F[4] 0 -0 =
97 F[5] 0 -0 =
98 F[6] 0 -0 =
99 F[7] 0 -0 =
100 F[8] 0 -0 =
101 F[9] 0 -0 =
No. Column name Activity Lower bound Upper bound
------ ------------ ------------- ------------- -------------
1 x[0,0,1] * 0 0 1
2 x[0,1,1] * 0 0 1
3 x[0,2,1] * 0 0 1
4 x[0,3,1] * 1 0 1
5 x[0,4,1] * 0 0 1
6 x[0,5,1] * 0 0 1
7 x[0,6,1] * 0 0 1
8 x[0,7,1] * 0 0 1
9 x[0,8,1] * 0 0 1
10 x[0,9,1] * 0 0 1
11 x[0,0,2] * 0 0 1
12 x[0,1,2] * 0 0 1
13 x[0,2,2] * 0 0 1
14 x[0,3,2] * 0 0 1
〜省略〜
325 x[8,-1,1] * 0 0 1
326 x[8,-1,2] * 1 0 1
327 x[8,-1,3] * 1 0 1
328 x[9,-1,1] * 0 0 1
329 x[9,-1,2] * 1 0 1
330 x[9,-1,3] * 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
この中で、
INTEGER OPTIMAL SOLUTION FOUND Time used: 0.0 secs Memory used: 0.5 Mb (511021 bytes)
という部分から、最適解が見つかって、それを見つけるのに0.0秒かかって（ほとんどかかってない！）、メモリを0.5MB使ったということが分かる。
そして、答えはx[k,i,n]がズラッと並んでいる部分のActivityの列を見ればいい。
ここで1になっている変数を拾い上げてみると、
- k=0
- x[0,3,1]
- x[0,-1,2]
- x[0,-1,3]
- k=1
- x[1,5,1]
- x[1,-1,2]
- x[1,-1,3]
- k=2
- x[2,6,1]
- x[2,-1,2]
- x[2,-1,3]
- k=3
- x[3,2,1]
- x[3,-1,2]
- x[3,-1,3]
- k=4
- x[4,2,1]
- x[4,-1,2]
- x[4,-1,3]
- k=5
- x[5,4,1]
- x[5,-1,2]
- x[5,-1,3]
- k=6
- x[6,2,1]
- x[6,-1,2]
- x[6,-1,3]
- k=7
- x[7,2,1]
- x[7,-1,2]
- x[7,-1,3]
- k=8
- x[8,1,1]
- x[8,-1,2]
- x[8,-1,3]
- k=9
- x[9,1,1]
- x[9,-1,2]
- x[9,-1,3]
という感じになっている。
これを解釈すると、答えは次のようになる。
+--------------------+ | 2015/07/15 | | 0 ( 3) 5 ( 4) | | 1 ( 5) 6 ( 2) | | 2 ( 6) 7 ( 2) | | 3 ( 2) 8 ( 1) | | 4 ( 2) 9 ( 1) | +--------------------+
一つ一つ確認していくと、ちゃんと正解が得られていることが分かると思う。
このように、パズルを最適化問題として定式化することが出来るし、ちゃんと定式化してやれば、ソルバーを使ってサクッと解くことが出来る。
そして、小さい問題（といっても制約が約100個、変数が約300個あるのだけど）なら、フリーソフトのGLPKでも簡単に解くことが出来る。
ということで、もっと最適化問題が一般的になればいいのになぁ、と。
今日はここまで！
