追記：参考リンク

はじめに

組合せ最適化(いわゆる混合整数最適化)では、さまざまなルールをシンプルに記述することができます。
ここでは、パズルを例題に、いくつかのテクニックを簡単にご紹介します。
これらのテクニックのいくつかは、Pythonによる部分もあります。数理最適化のモデル化においてPythonは相性がよいといえるでしょう。

環境構築

サクッと試したい場合：Binderというサービスを使うことにより、ブラウザだけで試せます。詳しくは、無料JupyterサービスのBinderの紹介をご覧ください。
きちんと試したい場合：Anacondaを導入後、下記を実行してください。

pip install pulp ortoolpy unionfind more_itertools

準備

モデラーとしては、PuLPを用います。(PuLPについて)
後述のコード例では、下記を省略します。

python

%matplotlib inline
import numpy as np, matplotlib.pyplot as plt
from itertools import groupby
from pulp import *
from unionfind import unionfind
from ortoolpy import addvar, addvars, addbinvar, addbinvars
m = LpProblem()

メソッド	説明
groupby	キーが同じものをグループ化する
unionfind	素集合データ構造用のクラス
addvar	変数を1つ作成
addvars	変数のリストを作成
addbinvar	0-1変数を1つ作成
addbinvars	0-1変数のリストを作成
LpProblem	PuLPの数理最適化モデル

対象パズル

SaitoTsutomu/OptForPuzzleに用意しました。

テクニック

変数をnp.arrayで作成

効用：演算の高速化、多彩なアクセス
対象パズル：「数独」など
例：

python

v = np.array(addbinvars(9, 9, 9)) # (1)
m += lpSum(v[i,:,k]) == 1 # (2)
m += lpSum(v[i:i+3,j:j+3,k]) == 1 # (3)
w = addvars(9, 9) # (4)
m += lpDot(range(9), v[i,j])+1 == r[i][j] # (4)

(1)では、9x9x9の多次元配列の0-1変数を作成しています。それぞれの次元は、行、列、数を表しています。
(2)では、$i$行目に数字$k+1$が1つだけを意味しています。
(3)では、左上が$(i,j)$となる3x3の領域で、数字$k+1$が1つだけを意味しています。
また、この$v$を使うと、(4)のように、マス$(i,j)$の数字を変数$w_{ij}$のように表現できます。

結果をnp.vectorizeで取得

効用：演算の高速化、多彩なアクセス
対象パズル：「数独」など(可視化は「黒どこ」など)
例：

python

r = np.vectorize(value)(v).astype(int) # (1)
plt.imshow(1-r, cmap='gray', interpolation='none') # (2)

変数をnp.arrayで作成し最適化を解いた結果は、上記の(1)のように結果$r$を取得すると便利です。こうすると、高速に簡単に結果を得ることができ、さらにnumpyの多彩な機能を使って処理を続けられます。
2次元の白黒の結果を図として確認する場合は、(2)のようにmatplotlibを使うと簡単に可視化できます。

なお、変数をpandasのDataFrameのSeriesで管理しているときは、applyで同様にできます。

すべて同じ

効用：効率的なモデル化
対象パズル：「ペイントエリア」など
例：

python

for vi, vj in zip(v, v[1:]):
    m += vi == vj

変数の1次元配列$v$のすべての要素が同じ値を取らないといけない場合、上記のように記述できます。
(変数そのものを置き換える方法もあります)

周りの数

効用：効率的なモデル化
対象パズル：「クリーク」「シャカシャカ」「数コロ」「のりのり」「ペイントエリア」「ボンバーパズル」

マスの変数$v[i,j]$の上下左右の変数の和を使いたいとき、下記のように$w$を使うことによりができます。
例：

python

u = np.zeros((nh+2, nw+2), dtype=object)
v = u[1:-1,1:-1] = np.array(addbinvars(nh, nw))
w = u[:-2,1:-1]+u[2:,1:-1]+u[1:-1,:-2]+u[1:-1,2:]

これは、vの周りに1周多い配列$u$を用意し、うまくスライスを使った例です。

IF制約

効用：条件によって成り立つ場合を表現
対象パズル：「カナオレ」など
変数$x,y$があったとき、「$x == 1$」ならば「$y \le a$」としたいときは、十分に大きな数$M$を使用して、下記のように書くことができます。
例：

python

m += y <= a + M(1-x)

また、「$x==1$」ならば「$y=a$」としたいときは、下記のように書くことができます。
例：

python

m += y <= a + M(1-x)
m += y >= a - M(1-x)

連結制約

効用：簡単なモデル化
対象パズル：「黒どこ」など

「黒どこ」のように「全ての白マスがつながっていること」という制約があります。このような制約を数式で表そうとすると、結構大変だったりします。そこで、一旦無視して、接続していなかったら、その解を禁止して解きなおすやり方があります。

例：連結性の確認

python

from unionfind import unionfind
r = np.array(
    [[0,0,0],
     [1,1,1],
     [0,0,0]])
print(unionfind.isconnected(1-r))
r[1][1] = 0
print(unionfind.isconnected(1-r))
>>>
False
True

結果の$r$は、0:白、1:黒を表しているとします。unionfind.isconnected(1-r) で白が連結しているかを確認できます。
これを踏まえると、下記のようにして連結になるまで、モデルを変更していくことができます。

python

while True:
    m.solve()
    r = np.vectorize(value)(v).astype(int)
    if unionfind.isconnected(1-r):
        break
    m += lpSum(v[r==1]) <= r.sum() - 1

上記のコードはシンプルですが、パズルによっては繰り返しが長くなることがあります。その場合は、「結果の全黒マスを禁止」するのではなく、黒マスで分断された白のエリアごとに「その白を囲っている黒マスを禁止」することにより、効率よく解くことができます。

エリア(部屋や国)ごとの処理

効用：シンプルな入力とプログラム
対象パズル：「因子の部屋」「カントリーロード」「さとがえり」「スターバトル」「タイルペイント」「チョコナ」「のりのり」「へやわけ」「ペイントエリア」
例：

python

data = """\
AABBB
AABCC
ADDDC
DDECC
EEEEC""".split()
v = np.array(addbinvars(nn, nn))
for _, d in groupby(sorted(zip(''.join(data), v.flat)), lambda x:x[0]):
    m += lpSum(c[1] for c in d) == 1

2次元のマスの表が上記 data のように、部屋に分かれています。このとき、「各部屋の中で黒マスは1つ」という条件は、groupbyを使うと上記のようにシンプルに書けます。

候補から選択

効用：簡単なモデル化
対象パズル：「カナオレ]「さとがえり]「四角に切れ]「ぬりかべ]「ののぐらむ]「フィルマット]「ブロックパズル]「ホタルビーム]

2つの0-1変数$x,y$があったとき、多くてもどちらか1方だけ成り立つ(=1)ようにしたければ、下記のようにできます。

例：

python

m += x + y <= 1

パズルのルールを数式で表すことが困難な場合、ルールを満たす組合せを列挙し、その中から選ばせることで簡単にモデル化できることがあります。このような定式化は典型問題の集合分割問題にあたります。

「カナオレ」の候補の作成関数を見てみましょう。
例：

python

def makecand(lst, y, x, d, l, p, u):
    yy, xx = y+[0, -1, 0, 1][d], x+[-1, 0, 1, 0][d] # (1)
    if 0 <= yy < nh and 0 <= xx < nw and (yy,xx) not in u:
        if p == l - 1:
            lst.append(u + [(yy, xx)])
            return
        for k in range(4):
            makecand(lst, yy, xx, k, l, p + 1, u + [(yy,xx)])

引数のdが方向になっています。(1)で 1マス追加し、pが長さが達したらlstに追加して終了。そうでなければ、4方向に探索を繰り返します。
「カナオレ」をこのような方法以外でモデル化するのは、困難でしょう。
AMPLなどのモデリング言語では、柔軟な記述ができません。このようにモデル化でPythonを使うのは、非常に有用といえます。

以上