クソコード製造機

数理最適化とかPythonとか

GurobiのLazy Constraintで制約を逐次追加する

Gurobi

組合せ最適化問題について、整数計画の定式化をしていると、表現したい条件についての制約が非常に多くなってしまうことがある。

このような場合には、必要な制約のみを逐次追加していく切除平面法が用いられることがある。

本記事では、それをGurobiのLazy Constraintという機能を用いて実装していく。

問題としては、無向グラフ G = (V, E)と辺の重み d_eについての最小全域木問題に対する、カットセット不等式を用いた定式化を考える。
また、問題例としては、重みをランダムに定めた格子グラフを生成する。
その生成コードは次の通り。

import networkx as nx
import random
import gurobipy as gp
from gurobipy import GRB

n, m = 5, 5
lb, ub = 1, 100

random.seed(42)
G = nx.grid_2d_graph(n, m)
V = n * m
nodes = list(G.nodes)

d = {}
for e in G.edges:
    d[e] = random.randint(lb, ub)

カットセット不等式を用いた定式化は、次のようになる。


\displaystyle \begin{aligned} &\mathrm{minimize} &&\sum_{e \in E} d_e x_e &\\ &\mathrm{s. t.} & &\sum_{e \in E} x_e = |V| - 1,& \\ &&& \sum_{e \in \delta(S)} x_e \geq 1 &\{\} \neq S \subsetneq V \\ &&& x_e \in \{0, 1\}, & e \in E \end{aligned}

これをGurobiで素直に実装すると次のようになる。

model = gp.Model()
x = model.addVars(G.edges, vtype=GRB.BINARY)

model.setObjective(x.prod(d))

model.addConstr(gp.quicksum(x) == V - 1)

for s in range(1, 1 << V - 1):
    S = [nodes[i] for i in range(V) if s & (1 << i)]
    cut_x = [x[(u, v)] for u, v in G.edges if (u in S and v not in S) or (v in S and u not in S)]
    x_cutsum = gp.quicksum(cut_x)
    model.addConstr(x_cutsum >= 1)

model.optimize()

カットセット不等式は頂点数の指数本になるため、頂点数が20~30程度で限界がくる。

よって、カットセット不等式を逐次追加していくことを考える。
方針は次のようにする。

  1. ある頂点 s \in Vを定める。
  2. 分枝限定法で暫定解が得られたときに、 S \subset Vを、得られた暫定解によって sから到達可能な頂点集合とする。
  3. V = Sならば、何もしない。そうでなければ、 Sに対応するカットセット不等式を追加する。

この実装は次のようになる。深さ優先探索で、到達可能な頂点集合を求めている。
model.cbLazyで制約を追加するためには、事前にパラメータ「LazyConstraints」を1に設定する必要があることに注意。

model = gp.Model()
model.setParam("LazyConstraints", 1)
x = model.addVars(G.edges, vtype=GRB.BINARY)

model.setObjective(x.prod(d))

model.addConstr(gp.quicksum(x) == V - 1)

def callback(model, where):
    if where == GRB.Callback.MIPSOL:
        x_cur = model.cbGetSolution(x)
        S = [(0, 0)]
        que = [(0, 0)]
        while len(que) > 0:
            v = que.pop()
            for u in G[v]:
                if u in S:
                    continue
                if ((u, v) in x.keys() and x_cur[(u, v)] >= 0.5) \
                        or ((v, u) in x.keys() and x_cur[(v, u)] >= 0.5):
                    S.append(u)
                    que.append(u)

        if len(S) < V:
            cut_x = [x[(u, v)] for u, v in G.edges if (u in S and v not in S) or (v in S and u not in S)]
            x_cutsum = gp.quicksum(cut_x)
            model.cbLazy(x_cutsum >= 1)


model.optimize(callback)

頂点数を64(8×8)に設定したときのGurobiのログの抜粋が以下のようになる。

ログを見ると、255本のカットセット不等式が追加されていることがわかる。

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0 1703.50000    0    4          - 1703.50000      -     -    0s
     0     0 1715.50000    0    4          - 1715.50000      -     -    0s
     0     2 1726.50000    0    8          - 1726.50000      -     -    0s
* 2910   462              13    1861.0000000 1831.00000  1.61%   2.3    0s
* 2923   462               9    1842.0000000 1831.00000  0.60%   2.3    0s
* 3068   464              11    1834.0000000 1833.00000  0.05%   2.3    0s
* 3541     0              12    1833.0000000 1833.00000  0.00%   2.4    0s

Cutting planes:
  Lazy constraints: 255

カットセット不等式は 2^{64} - 2本となるため、本質的に必要な制約が少ないことがわかる。
また、クラスカル法でもとめた重みと一致していることも確認できるため、最小全域木として必要な制約は追加できていることもわかる。

ただし、これでもMiller-Tucker-Zemlin制約を用いた定式化のような本質的に制約の本数が少ない定式化に比べると劣ってしまうため、定式化が難しい場合の最終手段ととらえるべきだろう。