混合整数計画法におけるテクニック例(Python-MIP)

はじめに

前回のブログ投稿で、PuLPとPython-MIPを使用して混合整数計画法(MIP)でTSPの問題を解いた。 少ない頂点(30頂点くらい)でないと短時間で探索を終えて最適解を出すことは難しかったが、 今回はいくつかのアプローチでそれよりも多い頂点で短時間で解を導出します。 以下のアプローチを使用しました。

• 制約の追加(TSPの最適解に含まれる可能性の高い辺である、コストの小さい辺のどれか使用する制約を加える)
• 近似解を初期解に加えて探索する
• Lazy Constraintsを使う

制約の追加

前回より制約を加えて探索する範囲を狭めることで高速化を行いました。 今回追加した制約はツアーの隣接頂点は頂点間距離の小さい上位の辺を必ず使用するという制約です。 TSPの最適解では遠くの頂点同士をつなぐということは滅多にありません。 ランダムに頂点の位置を決めたグラフであれば、ある頂点を端点とする上位20%の低コストの辺を使用することを制約に加えれば大体のケースで最適解を出せると思います。 ランダムな頂点でなくとも制約を弱めて上位50%とすればほぼどんなグラフでも最適解を出せるでしょう。 しかし完全グラフなどでは問題ありませんが、辺が少ない疎なグラフで頂点の場所が偏っているとのコーナーケースでは避けた方がいいと思います。

制約を加える具体的なコードは以下になります。 各頂点について辺のコストが小さい順にソートしてその上位のうちの1辺が使用されるので、 辺を使用するかどうかのバイナリ変数の和を1にする制約を加えていけばいいです。 コードでは範囲はrange(0, (n+1)//2)となっているため上位50%を使用します。

  for i in range(n):
    s = []
    for j in range(n):
        if i!=j:
            s.append((dists[i][j], j))
    s.sort()
    model += mip.xsum(x[i][s[k][1]] for k in range(0, (n+1)//2)) == 1
  

50頂点の完全グラフで比較を行いました。 グラフの構成の仕方は前回と同じくランダムに頂点のx,y座標を設定し、頂点間距離はユークリッド距離に小数点以下を除去したものを使用しました。 完全無向グラフとなります。

• 制約追加なし: 1000sec~(最適解は出たが1000secで終了せず)
• 上位50%のみ使用: 457sec
• 上位20%のみ使用: 143sec

いずれも同じ解が得られました。 たしかに制約の追加でより短時間で解を得ることができました。 ただ強い制約を加える方がより効果的な高速化につながりますが、 最適解の導出を妨げるような強い制約を加えないように注意が必要かと思います。

近似解を初期解に加える

分枝限定法では解が見つかるとそれを枝狩りに使用できますので、 より質の良い初期解であればあるほどさらに高速化します。 コードは以下を加えるだけです。 initsolの配列がTSP近似解のルートになります。 この近似解はシンプルな2-optとsimulated annealingを用いて導き出しました。 頂点数によってはsimulated annealingを何回か繰り返せば最適解が出てくるとは思いますが、それが最適解である保証はできません。

    model.start = [(x[initsol[i]][initsol[i+1]], 1) for i in range(n)]

それなりの最適解に近い解を使用しましたが、最適解を導き出すのは問題の規模によって時間がかかりました。 50頂点ではうまくいきましたが、100頂点TSPの導出は短時間ではできませんでした。

Lazy Constraints

Lazy Constraintsは最初は目的の制約より弱い制約で問題を解き、解が見つかるたびにその解が最終的な解の制約を満たさない場合に その許容解を禁止するような制約を後から課す手法です。 制約式が増えすぎて重くなってしまう場合や、解を見てから制約の計画を立てたいなどというときに用いるとよいと思います。 今回はsubtourを禁止する制約

$v_{i}-\left( n+1\right) x_{ij}\geq v_{j}-n$

を取り除きます。これは制約式が頂点の数に対して2乗に増えていくのでこれを抑えるためです。 subtourの制約が解かれるため解にはいくつかのsubtourが含まれている可能性がありますが、 subtourを有する解が見つかるたびにそのsubtourを禁止して除去する制約を加えます。 以下が加える制約になります。Cはsubtour頂点集合です。

$\sum _{\left( i,j\right) \in C}x_{ij}\leq \left| C\right| -1\forall C$

subtour自体は解からBFSのように探索すれば容易に取得可能です。 制約を加えるにはPython-MIPではmip.ConstrsGeneratorを使用します。 model.optimaize()の前にmodel.lazy_constrs_generator = SubTourLazyGenerator(x)の行を追加しLazyGeneratorを作成します。 SubTourLazyGeneratorは以下のコードになります。 公式のgithubのexampleにLazy Constraintsの使用例がありましたが、無駄な操作をして効率が悪い部分があるので高速化しました。 またLazy Constraintsを使用する場合は model.threads=-1とすると実行環境のせいかプログラムが落ちてしまったので、model.threads=2としました。

    class SubTourLazyGenerator(mip.ConstrsGenerator):
    def __init__(self, xv):
        self._x = xv

    def generate_constrs(self, model: mip.Model, depth: int = 0, npass: int = 0):
        x_, N, cp = model.translate(self._x), range(len(self._x)), mip.CutPool()
        g = [[j for j in N if x_[i][j].x >= 0.99] for i in N]

        vis = set()
        for v in N:
            if v in vis:
                continue
            tour = [v]
            vis.add(v)
            dq = deque()
            dq.append(v) 
            while dq:
                cur = dq.popleft()
                for nx in g[cur]:
                    if nx not in vis:
                        tour.append(nx)
                        dq.append(nx)
                        vis.add(nx)
            
            if len(tour) != len(N):
                AS = [(i, j) for (i, j) in product(tour, tour) if i != j]
                cut = mip.xsum(x_[i][j] for (i, j) in AS) <= len(tour) - 1
                cp.add(cut)

        for cut in cp.cuts:
            model += cut
    

lazy constraintsを使用したこのコードは100頂点のTSPを1000sec以内で解くことができました。 subtourの総数は頂点数に対して指数関数的ですが、 こちらが思った以上に効率的に制約の追加ができているようです。

まとめ

TSPを混合整数計画法で効率的に解くために制約の追加や近似解の導入、lazy constraintsを使用しました。 さわってみた感じではTSPでもそうですが他の問題を解く場合でも、問題の性質をよく把握して制約を考えることが必要であることが分かりました。 単純にソルバーに送り込むだけでは複雑な問題は安易に解けないということですね。 (TSPの場合はMIPを使用しなくともLKHといったもっと効率的な解法はありますが)