日々の流転

λ. 「異様に難しいドラクエの謎解き」を整数線形計画問題として解く

「異様に難しいドラクエの謎解き」をモデル検査器で解く の続き。 元の「人生を振り返る - 異様に難しいドラクエの謎解き」に「押す順番はどうでも良くて、四隅のボタンを押す回数の問題でしたね」とあって、この洞察を使えば整数線形計画問題として定式化して、もっとスマートに解けるなぁ、と思った。

整数線形計画問題による定式化

x1, x2, x3, x4 をそれぞれ左上、右上、右下、左下の石を押す回数とすると、最短の手順を求める問題は、以下の最適化問題として定式化できる。

Minimize
  x1 + x2 + x3 + x4
Subject To
  3 + x1 + x2 + x3 ≡ 0 mod 4  (左の像の条件)
  3 + x2 + x3 + x4 ≡ 0 mod 4  (上の像の条件)
      x3 + x4 + x1 ≡ 0 mod 4  (右の像の条件)
      x4 + x1 + x2 ≡ 0 mod 4  (下の像の条件)
  x1, x2, x3, x4 ≧ 0
  x1, x2, x3, x4 ∈ Z

追加的な変数を導入してmodを取り除けば、これは整数線形計画問題になる。その結果の問題をLPファイル形式にすると以下ようになる(baroque_tower.lp)。

Minimize
  x1 + x2 + x3 + x4
Subject To
  c1: x1 + x2 + x3 - 4 n1 = -3
  c2: x2 + x3 + x4 - 4 n2 = -3
  c3: x3 + x4 + x1 - 4 n3 = 0
  c4: x4 + x1 + x2 - 4 n4 = 0
Bounds
  0 <= x1
  0 <= x2
  0 <= x3
  0 <= x4
  0 <= n1
  0 <= n2
  0 <= n3
  0 <= n4
General
  x1 x2 x3 x4 n1 n2 n3 n4
End

MIPソルバで解いてみる

これをMIPソルバのlpsolveを使って解かせてみる。

> lp_solve.exe -rxli xli_CPLEX.dll baroque_tower.lp
set_XLI: Successfully loaded 'xli_CPLEX.dll'

Value of objective function: 6.00000000

Actual values of the variables:
x1                              1
x2                              2
x3                              2
x4                              1
n1                              2
n2                              2
n3                              1
n4                              1 

ということで、(x1,x2,x3,x4)=(1,2,2,1)という最適解が得られ、これは左上を1回、右上と右下を2回ずつ、左下を一回という、前回得られたものと同じ解が得られる。

実は最初はGLPKでも解こうと思ったけれど、GLPKだと、有界でない整数変数のせいで無限にブランチして止まらないっぽい。各変数に適当な上界を加えれば解ける。 あと、自分が以前に書いたオモチャのMIPソルバでも楽勝で解けた。

最適解の一意性

最適解の一意性については、MIPソルバを使うよりも、論理結合子を自由に使えるSMTソルバの方が便利なので、SMTソルバのZ3で確認してみることにする。

以下の内容(baroque_tower.smt2)をZ3に与える（オンラインでも Z3 online at RiSE4Fun で試せる）ことで、6手の他の解が存在しないことを確認できる。 詳しくは以下のコメント中にて。

; baroque_tower.smt2
(set-option :produce-models true)
(set-logic QF_LIA)

; 変数の宣言
(declare-fun x1 () Int)
(declare-fun x2 () Int)
(declare-fun x3 () Int)
(declare-fun x4 () Int)

; 前述の制約条件の追加
(assert (>= x1 0))
(assert (>= x2 0))
(assert (>= x3 0))
(assert (>= x4 0))
(assert (= (mod (+ 3 x1 x2 x3) 4) 0))
(assert (= (mod (+ 3 x2 x3 x4) 4) 0))
(assert (= (mod (+ x3 x4 x1) 4) 0))
(assert (= (mod (+ x4 x1 x2) 4) 0))

;; まずは適当に解を求めてみる
(check-sat)
; sat と出力され、解が存在することが分かるので、解を表示させる
(get-value (x1 x2 x3 x4))
; ((x1 1)
;  (x2 4946)
;  (x3 4990)
;  (x4 9749))

;; 既に最適解での目的関数値が6と分かっているので、その条件を加えて再度解いてみる  
(assert (= (+ x1 x2 x3 x4) 6))
(check-sat)
; 今度もsatと出力され、解が存在することが分かるので、 解を表示させる
(get-value (x1 x2 x3 x4))
; ((x1 1)
;  (x2 2)
;  (x3 2)
;  (x4 1))

;; 更に、今見つかった最適解を除外する制約を加えて、再度解いてみる
(assert (not (and (= x1 1) (= x2 2) (= x3 2) (= x4 1))))
(check-sat)
; 今度はunsatと出力され、従って今除外された以外の解がないことが分かる

関連リンク

• 『「異様に難しいドラクエの謎解き」をモデル検査器で解く』代わりに数学で解く - Elliptium

ひょっとしたら関係するかもしれない書籍など

λ. ゴモリーの小数カットの導出メモ

Outline of an algorithm for integer solutions to linear programs でのゴモリーカット(Gomory's fractional cut)が、割と天下りに定義して、そのあと必要な性質を満たすことを証明しているけれど、これがちょっと気に入らなかったので、演繹的に導出してみた。(An algorithm for the mixed integer problem のゴモリーの混合整数カットの方ではないので注意)

準備

N を(0を含む)自然数の集合とし、(R, ≤_R), (Z, ≤_Z) をそれぞれ実数と整数の順序集合とする。 また、 I : (Z, ≤_Z) → (R, ≤_R) を順序集合間の通常の埋め込みとするが、たいてい省略する。

⌈-⌉, ⌊-⌋ : (R, ≤_R) → (Z, ≤_Z) をそれぞれ正の無限大方向への切り上げ(ceiling)と負の無限大方向への切り捨て(floor)とする。

随伴、モナド、コモナド、双対性

⌈-⌉ と ⌊-⌋ はそれぞれ I の左随伴と右随伴になっている。すなわち、以下が成り立つ。

• ∀a∈R, b∈Z. ⌈a⌉ ≤_Z b ⇔ a ≤_R I(b)
• ∀a∈Z, b∈R. I(a) ≤_R b ⇔ a ≤_Z ⌊b⌋

したがって、 I∘⌈-⌉, I∘⌊-⌋ : (R, ≤_R) → (R, ≤_R) はそれぞれ (R, ≤_R) 上のモナド、コモナドとなっており、以下が成り立つ。

• ⌊a⌋ ≤ a ≤ ⌈a⌉
• ⌊⌊a⌋⌋ = ⌊a⌋
• ⌈a⌉ = ⌈⌈a⌉⌉

また、 ⌈-⌉∘I, ⌊-⌋∘I : (Z, ≤_Z) → (Z, ≤_Z) はそれぞれ (Z, ≤_Z) 上のコモナド、モナドとなっており、こっちは更に強く、以下が成り立っている。

• ⌈I(a)⌉ = a for all a∈Z
• a = ⌊I(a)⌋ for all a∈Z

また、(-1)×- : (R, ≤_R) → (R,≥_R) によって、⌈-⌉と⌊-⌋は双対になっており、 ⌈a⌉ = - ⌊- a⌋ が成り立つ。

⌈-⌉ と ⌊-⌋ の分配性

a≤⌈a⌉ および b≤⌈b⌉ より a+b ≤_R ⌈a⌉+⌈b⌉ = I(⌈a⌉+⌈b⌉) で、随伴性より ⌈a+b⌉ ≤_Z ⌈a⌉+⌈b⌉ 。 双対性より ⌊a+b⌋ ≥_Z ⌊a⌋+⌊b⌋ 。

さらに、帰納法により ∀n∈N. ⌈n×a⌉ ≤_Z n×⌈a⌉ もいえる。 同じく、双対性より ∀n∈N. ⌊n×a⌋ ≥_Z n× ⌊a⌋

ゴモリーの小数カットの導出

x_i, y ∈ N が

• Σ_i a_i x_i + y = b …… (1)

を満たすとする。すると、

• Σ_i a_i x_i + y ≥ b … (2)
• Σ_i a_i x_i + y ≤ b … (3)

(3)の両辺に⌊-⌋を適用すると、⌊-⌋ の単調性より

• ⌊Σ_i a_i x_i + y⌋ ≤ ⌊b⌋ … (4)

⌊-⌋ の分配性と ⌊y⌋=y より

• ⌊Σ_i a_i x_i + y⌋ ≥ Σ_i ⌊a_i x_i⌋ + ⌊y⌋ ≥ Σ_i ⌊a_i⌋ x_i + ⌊y⌋ = Σ_i ⌊a_i⌋ x_i + y … (5)

(4), (5)に推移律を適用して

• Σ_i ⌊a_i⌋ x_i + y ≤ ⌊b⌋ … (6)

両辺に -1 をかけて

• Σ_i (- ⌊a_i⌋) x_i - y ≥ -⌊b⌋ … (7)

(2) と (7) の両辺を足し合わせて、

• Σ_i (a_i - ⌊a_i⌋) x_i ≥ b - ⌊b⌋

■

λ. 強双対性を使って線形計画を制約充足問題に還元する

Linear and Discrete Optmization の双対性に関する部分(参考: 第五週目メモ)を復習していて気になったのだけど、強双対性の条件 c^T x^* = b^T y^* って線形式で書けてしまっているのよね。 ということは主問題と双対問題の制約条件と強双対性の連言をとっても線形(QF_LRA)の範囲に収まっていて、これの解(実行可能解)を求めれば、目的関数に関する最適化をしなくても最適解が求まることになるはず。 ということで試してみた。

主問題

簡単な問題として、Linear and Discrete Optmization で最初の例として使われていた、以下の例を考える。

MAXIMIZE 
  100 x1 + 125 x2
SUBJECT TO
  3 x1 + 6 x2 <= 30
  8 x1 + 4 x2 <= 44
BOUNDS
  0 <= x1 <= 5
  0 <= x2 <= 4
END

最適解は (x1,x2)=(4,3) で、目的関数値は775。

双対問題

この問題の双対問題は以下のようになる。

MINIMIZE
  30 y1 + 44 y2 + 5 y3 + 4 y4
SUBJECT TO
  3 y1 + 8 y2 + y3 >= 100
  6 y1 + 4 y2 + y4 >= 125
BOUNDS
  0 <= y1
  0 <= y2
  0 <= y3
  0 <= y4
END

最適解は (y1,y2,y3,y4) = (50/3, 25/4, 0, 0) で、目的関数値は775。

両方を合わせたもの

MAXIMIZE 
  0 x1
SUBJECT TO
  3 x1 + 6 x2 <= 30
  8 x1 + 4 x2 <= 44
  3 y1 + 8 y2 + y3 >= 100
  6 y1 + 4 y2 + y4 >= 125
  100 x1 + 125 x2 - 30 y1 - 44 y2 - 5 y3 - 4 y4 = 0
BOUNDS
  0 <= x1 <= 5
  0 <= x2 <= 4
  0 <= y1
  0 <= y2
  0 <= y3
  0 <= y4
END

(x1,x2,y1,y2,y3,y4) = (4, 3, 50/3, 25/4, 0, 0) はこの問題の(実行可能)解で、このときの目的関数値はどちらも775。

Z3の場合

最適化機能を持たない充足可能性ソルバの例として、SMTソルバのZ3で試してみる（オンラインでも Z3 online at RiSE4Fun で試せる）。

(declare-fun x1 () Real)
(declare-fun x2 () Real)
(define-fun obj-primal () Real
  (+ (* 100.0 x1) (* 125.0 x2)))
(assert (<= (+ (* 3.0 x1) (* 6.0 x2)) 30.0))
(assert (<= (+ (* 8.0 x1) (* 4.0 x2)) 44.0))
(assert (and (<= 0.0 x1) (<= x1 5.0)))
(assert (and (<= 0.0 x2) (<= x2 4.0)))

(declare-fun y1 () Real)
(declare-fun y2 () Real)
(declare-fun y3 () Real)
(declare-fun y4 () Real)
(define-fun obj-dual () Real
  (+ (* 30.0 y1) (* 44.0 y2) (* 5.0 y3) (* 4.0 y4)))
(assert (>= (+ (* 3.0 y1) (* 8.0 y2) y3) 100.0))
(assert (>= (+ (* 6.0 y1) (* 4.0 y2) y4) 125.0))
(assert (<= 0.0 y1))
(assert (<= 0.0 y2))
(assert (<= 0.0 y3))
(assert (<= 0.0 y4))

(assert (= obj-primal obj-dual))

(check-sat) ; => sat

(get-value (obj-primal x1 x2))
; => ((obj-primal 775.0) (x1 4.0) (x2 3.0))

(get-value (obj-dual y1 y2 y3 y4))
; => ((obj-dual 775.0) (y1 (/ 50.0 3.0))
;     (y2 (/ 25.0 4.0)) (y3 0.0) (y4 0.0))

感想とか

整数上の場合だと「解が一つ見つかった後に、目的関数値がそれよりも真に良いということを表す制約を追加して、再度解き直す」ということを繰り返すことで、有界な場合には最適化を行えるのに対して、実数上だと有界であってもZeno状態みたいなことに陥って、最適解に到達できない場合があるはず。 なので、実数の場合には充足可能性ソルバでは最適化問題は解けないと思っていたけれど、こういうやり方をすれば解けるんだなぁ。 ということで面白かった。

まあ、浮動小数点数で計算するような場合には、誤差蓄積で解が見つからなくなるような場合もあるような気はするけれど。