👌

MIPソルバーによる容量制約付き配送計画問題の解法についての簡単なまとめ

2020/12/27に公開

容量制約付き配送計画問題とは

複数のトラックを使って荷物を集積所（デポ）から各店舗に配送するときに、どのトラックがどの店舗をどの順番で回るのが最適か、という問題を配送計画問題と呼びます。何をもって「最適」とするかは問題に依存しますが、移動費用や車両数とすることが多いようです。

各トラックには積載量の上限があるので「各トラックの運ぶ荷物は積載容量を超えない」という制約を課したものは特に容量制約付き配送計画問題（capacitated vehicle routing problem, CVRP）と呼ばれており、配送計画問題の中で最も基本的な問題です。

混合整数計画法による定式化と部分順回路除去

ここでは、デポと店舗をあわせて $N$ 箇所あるものとし、インデックス $i, j \in \{0, 1, \dots, N\}$ で表すことにします。また、デポを表すインデックスを $0$ とします。

$C_{ij}$ を $i$ から $j$ への移動費用、 $x_{ij}\in\{0, 1\}$ を店舗（またはデポ） $i$ から $j$ に直接移動する場合に $1$ となる変数とすると、移動費用の合計を最小化する目的関数は以下の様に書けます。

\min_{x}\quad \sum C_{ij}x_{ij},

各店舗に入っていくトラックの数と出ていくトラックの数は等しくなります。また、1つの店舗を複数のトラックが訪れる必要はない^[1]ことを考えると、 $x_{ij}$ には以下の制約が必要になることがわかります。

\quad \sum_{j}x_{ij} = \sum_{j}x_{ji} = 1, \quad i\neq 0

この制約がないと、トラックは途中で分裂したり消滅したりしてしまいます。

さらにトラックの台数 $M$ が決まっている場合は、デポに対して以下の制約が追加されます。

\quad \sum_{j}x_{0j} = \sum_{j}x_{j0} = M.

実は、ここまでの制約ではまだ足りません。前述の通り容量制約を考慮する必要があるのに加えて、以下のような部分巡回路が発生してしまう可能性を排除する必要があります。

部分順回路除去と変数と制約の数

この部分巡回路除去については、いくつかの方法が提案されています。定式化の方法によって追加の変数の数や制約数が異なるので、参考までにまとめてみました。

定式化	追加の変数の数	部分巡回路除去に関わる制約数
部分巡回路除去定式化	$0$	$2^{N-1}$
ポテンシャル定式化	$N$	$N(N-1)$
単品種フロー定式化	$N(N-1)$	$N$ ^[2]
多品種フロー定式化	$N^2(N-1)$	$N(N-1)$ ^[3]
集合分割定式化^[4]	$(N-1)!$	0

以降では、それぞれの定式化について具体的に見ていきます。

部分巡回路除去定式化

部分巡回路除去定式化^[5]では、追加の変数を用いずに部分巡回路除去制約と容量制約を実現します。

部分巡回路除去定式化の詳細

\sum_{i, j\in S}x_{ij} \leq |S| - V(S), \quad ^\forall S \subset \{1,2,\dots,N\}

ここで $S$ は（デポを除く）店舗集合の部分集合、 $V(S)$ は $S$ に含まれる店舗すべての需要を賄うのに必要なトラックの台数を表します。下図を見ると分かる通り、デポを除く任意の部分集合 $S$ について、エッジの総数 $\sum_{i, j\in S}x_{ij}$ がトラック1台で運ぶ場合は $|S| - 1$ 、トラック2台で運ぶ場合は $|S| - 2$ となることから導ける制約です^[6]。

この方法では、すべての部分集合 $S$ に対して制約が必要なので、制約の数は $2^{N-1}$ になります^[7]。そのため、 $N$ が大きくなると制約の数が急激に増えてしまいます。例えば $N=30$ のとき、 $2^{29} = 536870912$ 個の制約が必要となり、現実的ではありません。そこで、後述する必要な制約だけを逐次追加していく方法が提案されています。

ポテンシャル定式化

ポテンシャル定式化^[8]では、ポテンシャルという概念を各店舗に導入し、部分巡回路除去制約と容量制約を実現します。ポテンシャルは物理で出てくるポテンシャル（位置エネルギー）と同じ概念で、店舗を訪れる毎に減っていく（もしくは増えていく）なにか、と考えるとわかりやすいと思います。部分巡回路除去のためにはどんなポテンシャルを使ってもよいのですが、CVRPでは各トラックの累積の配送量を使うと容量制約も同時に課せるので一石二鳥です。

ポテンシャル定式化の詳細

\begin{aligned} u_i + D_j - Q(1 - x_{ij}) &\leq u_{j}, \quad j \neq 0\\ D_i \leq u_i &\leq Q. \end{aligned}

ここで、 $u_i$ がポテンシャルで、店舗 $i$ を出発する時までに運んだ配送量を表します。また、 $D_{i}$ は店舗 $i$ の需要、 $Q$ はトラックの容量です。 $x_{ij} = 1$ 、つまりトラックが $i$ から $j$ に直接行く時、 $u_{i} + D_{j} \leq u_j$ となり、 $j$ までの配送量は、 $i$ までの配送量と $j$ の需要をたしたものを表していることがわかります。逆に $x_{ij} = 0$ のとき、 $u_i + D_j - Q \leq u_{j}$ となりますが、 $D_i \leq u_i \leq Q$ より、これは常に成り立つので、追加の制約になはりません。

さて、この定式化は下界が非常に悪いことが知られています^[9]。そこで、いろいろな不等式（例えば $x_{ij} = 1$ なら $x_{ji} = 1$ など、問題の構造を考えるとわかる条件）をつかってギリギリを攻めていくことを考えます。この操作のことを強化と呼びます。強化した制約式は、~~まったくよくわからない~~前述のものほどは解釈性がよくありませんが、求解時間を短くなる可能性があります。

\begin{aligned} &u_{i} + D_{j} - Q(1 - x_{ij}) + (Q - D_i - D_j)x_{ji} \leq u_j & j\neq 0\\ &u_i \geq D_i + \sum_{j}D_jx_{ji} \\ &u_i \leq Q - \sum_{j}D_{j}x_{ij} \\ &u_{i} \leq Q - (Q - D_i)x_{0i} \\ &u_{i} \leq Q - (Q - \max_{j\neq i} D_j - D_i)x_{0i} - \sum_{j\neq i}D_{j}x_{ij} \\ &D_{i} \leq u_i \leq Q. \end{aligned}

単品種フロー定式化

単品種フロー定式化では、各店舗ではなく店舗と店舗をつなぐエッジに条件付けすることで、巡回路除去制約を実現します。ポテンシャル定式化では各店舗までの累積の配送量に注目したのに対して、単品種フロー定式化では、各エッジの上を通る荷物の積載量に注目します。

単品種フロー定式化の詳細

\begin{aligned} \sum_{j}f_{ji} - \sum_{j}f_{ij} &= D_{i}, \\ 0 \leq f_{ij} &\leq Qx_{ij}. \end{aligned}

ここで、 $f_{ij}$ がエッジ $(i, j)$ の上を通る荷物の積載量です。1つめの式は店舗 $i$ で需要 $D_{i}$ だけ積荷が降ろされることを、2つめの式は、トラックが通らければ荷物も通らないことを表します。

ポテンシャル定式化と同様、多少強化することができて、

\begin{aligned} \sum_{j}f_{ji} - \sum_{j}f_{ij} &= D_{i}, \\ 0 \leq f_{ij} &\leq (Q - D_i)x_{ij}. \\ \end{aligned}

となります。

D_{j}x_{ij} \leq f_{ij}

もいえるはずですが、実験してみた限りではあまり速くはならないようです。

多品種フロー定式化

多品種フロー定式化では、単品種フロー定式化と同様に各エッジを通る荷物の積載量を考えますが、「どの店舗むけの荷物か」の情報も利用します。

多品種フロー定式化の詳細

\begin{aligned} \sum_{j}f^{k}_{ji} - \sum_{j}f^{k}_{ij} &= 0,& k\neq i,\ i\neq 0\\ \sum_{j}f^{i}_{ji} - \sum_{j}f^{i}_{ij} &= D_{i}, \\ \sum_{j}f^{k}_{j0} - \sum_{j}f^{k}_{0j} &= -D_{k}, \\ \sum_{k}f^{k}_{ij} &\leq Qx_{ij}, \\ 0 \leq f^{k}_{ij} &\leq D_{k}. \\ \end{aligned}

ここで $f^{k}_{ij}$ はエッジ $(i, j)$ の上を通る店舗 $k$ 向けの荷物の積載量です。
こちらも多少強化することができて、

\begin{aligned} \sum_{j}f^{k}_{ji} - \sum_{j}f^{k}_{ij} &= 0,& k\neq i,\ i\neq 0\\ \sum_{j}f^{i}_{ji} - \sum_{j}f^{i}_{ij} &= D_{i}, \\ \sum_{j}f^{k}_{j0} - \sum_{j}f^{k}_{0j} &= -D_{k}, \\ D_{j}x_{ij} \leq \sum_{k}f^{k}_{ij} &\leq (Q - D_i)x_{ij}, \\ 0 \leq f^{k}_{ij} &\leq D_{k}. \end{aligned}

となります。

集合分割定式化

集合分割定式化は少し特殊な定式化です。ルート（一つのトラックの走る経路）候補をすべて列挙しておいて、その中から適切な組合せを見つけ出します。

当然、どうやってルートを列挙するか、という問題に出くわします。一般の配送計画問題では色々な制約があるために、ルート候補数が列挙可能な数に限られていることがあり、そういった場合に有効なアプローチのようです。また、ルート候補数がそこまで絞り込めない場合は、良さそうなルートの候補を徐々に追加していくアプローチを取ります。この方法を列生成法と呼びます。また、分枝限定法の中に列生成法を組み込む分枝価格法も提案されています。

配送計画問題に対する列生成法のもう少し詳しい説明については「はじめての配送計画の列生成法」を御覧ください。

集合分割定式化の詳細

\begin{aligned} \min \sum_{r}C_{ir}z_{r} & \\ \sum_{r}A_{ir}z_{r} &= 1, \\ z_{r} &\in \{0, 1\}. \end{aligned}

あるいは少し条件を緩めて^[10]

\begin{aligned} \min \sum_{r}C_{ir}z_{r} & \\ \sum_{r}A_{ir}z_{r} &\geq 1, \\ z_{r} &\in \{0, 1\}. \end{aligned}

とします。ここで $z_{r}$ が、ルート $r$ を使う場合に $1$ となる変数で、 $A_{ir}$ はルート $r$ に店舗 $i$ が含まれている場合に $1$ となる定数です。ルートを事前に列挙しているため、他の定式化に出てきた $x_{ij}$ は出てきません。

多項式サイズのアルゴリズムについての補足

部分巡回路除去定式化と集合分割定式化を除く、変数と制約の数が高々多項式サイズの定式化については、Aksen et al. (2018) で詳しく比較されているので、興味なる方はそちらもご覧ください。問題にもよるので一概には言えませんが、単品種フロー定式化に車両数の上限と下限を加えると高速に解ける傾向があるようです。

試しに標準的なベンチマークセットであるCVRPLibの set_A/A-n32-k5 を使って求解時間を比較してみました。厳密な評価をしたわけではありませんが、定式化によって求解時間が大きく異なることはわかるかと思います。

定式化	求解時間
ポテンシャル定式化	> 10min
単品種フロー定式化	41.7sec
多品種フロー定式化	> 10min

部分巡回路除去制約の逐次追加と分枝カット法

さて、部分巡回路除去定式化において必要な制約だけを逐次的に追加していく最も簡単な方法は、とりあえず部分巡回路除去制約を無視して最適解を求めて、部分巡回路ができたら、その巡回路を除去する制約を追加していくことです^[11]。ここでは、この方法を切除平面法（簡便法）と呼ぶことにします。

切除平面法（簡便法）をPuLPで実装すると、以下のようになります。前述のとおり、必要がなくなるまで制約を追加しては求解を繰り返します。

最適化の実行

model = build_model(n_nodes, dists, n_vehicles)
while True:
    assert model.solve() == pulp.LpSolutionOptimal
    # 制約を追加
    if add_cuts(model, demands, capacity) == 0:
        break

build_model では、部分巡回路除去制約を無視したモデルを構築しています。

部分巡回路を無視したモデルの構築

def build_model(n_nodes, dists, n_vehicles):
    nodes = list(range(n_nodes))

    model = pulp.LpProblem()
    x = np.array(pulp.LpVariable.matrix('x', (nodes, nodes), cat=pulp.LpBinary))
    # xの対角成分を 0 に固定する
    np.vectorize(lambda x, a: x.bounds(a, a)).fix_value(np.diagonal(x), 0)

    model.setObjective(pulp.lpSum(dists*x))

    for i in nodes[1:]:
        model.addConstraint(pulp.lpSum(x[i, :]) == 1)
        model.addConstraint(pulp.lpSum(x[:, i]) == 1)  
    model.addConstraint(pulp.lpSum(x[0, :]) == n_vehicles)
    model.addConstraint(pulp.lpSum(x[:, 0]) == n_vehicles)

    model._x = x
    return model

add_cuts では、NetworkX を利用して巡回路を抽出し、必要な制約を追加しています。

制約の追加

def add_cuts(model, demands, capacity):
    # NetworkX のグラフに変換
    graph = nx.from_numpy_array(
        np.vectorize(lambda x: x.value())(model._x)
    )
    # 各トラックのルートと部分巡回路を取得
    graph.remove_node(0)
    tours = list(nx.connected_components(graph))

    n_cuts = 0
    for tour in tours:
        tour = list(tour)
        if len(tour) <= 1:
            continue

        # エッジ数と、必要なトラックの台数を計算
        n_edges = pulp.lpSum(model._x[tour, :][:, tour])
        v = math.ceil(demands[tour].sum()/capacity)

        constraint = n_edges <= len(tour) - v
        if not constraint.valid():
            model.addConstraint(constraint)
            n_cuts += 1
    return n_cuts

set_A/A-n32-k5 の結果は以下のとおりです。

定式化	求解時間
切除平面法（簡便法）	> 10min

実は切除平面法（簡便法）は、巡回セールスマン問題に対しては効率的な方法であることが知られています。そこで、店舗数や需要量はそのままで積載容量をかえながら実験してみると以下のようになりました。（set_A/A-n32-k5 の店舗の需要量の合計は440です。）

積載容量	台数	求解時間
450	1	1.45sec
400	2	1.43sec
350	2	1.64sec
300	2	1.91sec
250	2	1.90sec
200	3	3.71sec
150	3	> 10min
100	5	> 10min

必要なトラック台数が1台のときは巡回セールスマン問題とおなじなので、非常に高速に動作します。一方で必要なトラック台数が増えると求解しなおす回数が増え、急激に時間がかかってしまう傾向があるようです。ちなみに、厳密解法ではなく近似解法であると捉え、上界と下界との差が一定以下となったら計算を打ち切るようにすると以下のようになりました。近似解法であると割り切れば、条件をかえずとも現実的な時間で求解できるようです。

ギャップ	距離	距離/最適値	求解時間
5%	792	1.006	113s
10%	823	1.045	32s
20%	912	1.16	47s

さて、上記の実装では、新しい制約を追加するたびに、最適化問題をはじめから解き直していました。商用ソルバーの Gurobi を使うと、求解途中にコールバックを使って制約を追加できます。この方法は分子カット法と呼ばれています。コードは以下の様になります ^[12]。

部分巡回路を無視したモデルの構築

def build_model(n_nodes, dists, demands, capacity, n_vehicles):
    nodes = list(range(n_nodes))

    model = gp.Model('cvrp.cutset_gurobi')
    
    x = {}
    for i in nodes:
        for j in nodes:
            if i == j:
                continue
            x[i, j] = model.addVar(vtype=gp.GRB.BINARY, name=f'x_{i}_{j}')

    model.setObjective(gp.quicksum(dists[k]*v for k, v in x.items()), gp.GRB.MINIMIZE)

    for i in nodes[1:]:
        model.addConstr(gp.quicksum(x[i, j] for j in nodes if i != j) == 1)
        model.addConstr(gp.quicksum(x[j, i] for j in nodes if i != j) == 1)
    model.addConstr(gp.quicksum(x[0, j] for j in nodes[1:]) == n_vehicles)
    model.addConstr(gp.quicksum(x[j, 0] for j in nodes[1:]) == n_vehicles)

    model._x = x
    model._capacity = capacity
    model._demands = demands
    return model

コールバック関数

def callback(model, where):
    # 新たな解が求まった時に、必要な制約を追加する
    if where != gp.GRB.Callback.MIPSOL:
        return

    edges = []
    for (i, j) in model._x:
        if model.cbGetSolution(model._x[i, j]) > 0.5:
            edges.append((i, j))
    graph = nx.from_edgelist(edges, create_using=nx.DiGraph)
    for tour in nx.simple_cycles(graph):
        dept_contained = 0 in tour
        if dept_contained:
            tour.remove(0)

        # 必要な車両台数とエッジ数を計算
        n_vehicles = math.ceil(model._demands[tour].sum()/model._capacity)
        n_edges = gp.quicksum(
            model._x[i, j]
            for i in tour
            for j in tour if i != j
        )
        if len(tour) >= 2 and (not dept_contained or n_vehicles > 1):
            # コールバック関数内では cbLazy を使って制約を追加できる
            model.cbLazy(n_edges <= len(tour) - n_vehicles)

最適化計算の実行

    model = build_model(n_nodes, dists, demands, capacity, n_vehicles)
    # コールバックで制約を追加するには LazyConstraints オプションを設定する必要がある
    model.setParam('LazyConstraints', 1)
    # コールバックを指定して最適化
    model.optimize(callback=callback)

結果は以下のとおりです。切除平面法（簡便法）と違い、かなり高速に動作するようです。

定式化	求解時間
分枝カット法	47.6sec

集合分割定式化と列生成法

前述の通り、集合分割定式化では変数（1台のトラックのルート）を列挙する必要がありますが、ルートの数は $(N-1)!$ となり、実際に列挙するのは不可能です。そこで、システマティックに候補を列挙する方法が考案されていますが、その一つが列生成法です。列生成法では、一部の変数（ルート）しか考慮に入れていないような限定的な集合分割定式化を考え、その線形緩和問題の双対においては制約を追加することがルート候補を追加することに対応しているという性質を利用します。
具体的な手順は以下のとおりで、

初期ルートを設定
候補が生成されなくなくなるまでルートを生成する
1. 双対問題を解いて子問題を生成
2. 子問題を解いてルート候補を生成
ルート候補をもとに最適解を求める

というステップで最適化を行います。コードは以下の様になります。

# モデル
dual_model = build_dual_model(n_nodes, dists, n_vehicles)
slave_model = build_slave_model(n_nodes, demands, capacity, dists, n_vehicles)

# 初期ルートを設定
routes = []
for i in range(1, n_nodes):
    route = [0, i]
    routes.append(route)
    update_dual_model(dual_model, route)

# ルートが生成されなくなるまで繰り返す
while True:
    dual_model.optimize()
    y_opt = [y.x for y in dual_model._y.values()]
    route = generate_route(slave_model, y_opt)
    if route is None:
        break
    routes.append(route)
    update_dual_model(dual_model, route)

primal_model = build_primal_model(n_nodes, routes, dists, n_vehicles)
primal_model.optimize()

双対問題は以下のようになります。

def build_dual_model(n_nodes, dists, n_vehicles):
    nodes = list(range(n_nodes))
    model = gp.Model('cvrp.cgen_gurobi.master')
    
    y = {}
    for i in nodes:
        y[i] = model.addVar(lb=0, vtype=gp.GRB.CONTINUOUS, name=f'y_{i}')

    model.setObjective(gp.quicksum(y.values()) + (n_vehicles - 1)*y[0], gp.GRB.MAXIMIZE)

    model._y = y
    model._dists = dists
    return model

双対問題へのルート候補の追加

def update_dual_model(model, route):
    model.addConstr(gp.quicksum(model._y[i] for i in route) <= cost(route, model._dists))
    model.update()

子問題はに負の距離を許した最短経路問題になります。子問題はループの中で何回も解き直すことになるので、以下に高速に求解するかがポイントとなりますが、今回は分枝カット法で実装しています。

def build_slave_model(n_nodes, demands, capacity, dists, n_vehicles):
    nodes = list(range(n_nodes))

    model = gp.Model()
    model.setParam('LazyConstraints', 1)
    
    x = {}
    for i in nodes:
        for j in nodes:
            if i == j:
                continue
            x[i, j] = model.addVar(vtype=gp.GRB.BINARY, name=f'x_{i}_{j}')
    w = {}
    for i in nodes[1:]:
        w[i] = model.addVar(vtype=gp.GRB.BINARY, name=f'w_{i}')

    for i in nodes[1:]:
        model.addConstr(gp.quicksum(x[i, j] for j in nodes if i != j) == w[i])
        model.addConstr(gp.quicksum(x[j, i] for j in nodes if i != j) == w[i])
    model.addConstr(gp.quicksum(x[0, j] for j in nodes[1:]) == 1)
    model.addConstr(gp.quicksum(x[j, 0] for j in nodes[1:]) == 1)

    model.addConstr(gp.quicksum(demands[i]*w[i] for i in w) <= capacity)

    model._x = x
    model._w = w
    model._dists = dists
    model._capacity = capacity
    model._demands = demands
    model._n_vehicles = n_vehicles
    return model


def slave_callback(model, where):
    if where != gp.GRB.Callback.MIPSOL:
        return

    edges = []
    for (i, j) in model._x:
        if model.cbGetSolution(model._x[i, j]) > 0.5:
            edges.append((i, j))
    graph = nx.from_edgelist(edges, create_using=nx.DiGraph)
    for tour in nx.simple_cycles(graph):
        dept_contained = 0 in tour
        if dept_contained:
            tour.remove(0)

        n_edges = gp.quicksum(
            model._x[i, j]
            for i in tour
            for j in tour if i != j
        )
        if len(tour) >= 2:
            model.cbLazy(n_edges <= len(tour) - 1)

子問題を使って使ってルート候補を生成します。

def generate_route(slave_model, y_opt, timelimit):
    start = time.time()
    cost = gp.quicksum(slave_model._dists[i, j]*slave_model._x[i, j] for i, j in slave_model._x)
    objective =  cost - gp.quicksum(y_opt[i]*slave_model._w[i] for i in slave_model._w) - slave_model._n_vehicles*y_opt[0]

    slave_model.setObjective(objective)
    slave_model.update()
    slave_model.optimize(callback=slave_callback)
    
    edges = []
    for (i, j), x in slave_model._x.items():
        if x.x >= 0.5:
            edges.append((i, j))
    graph = nx.from_edgelist(edges, create_using=nx.DiGraph)
    return next(nx.simple_cycles(graph))

主問題では生成されたルート候補のなかから最適な組合せを探します。

def build_primal_model(n_nodes, routes, dists, n_vehicles):
    model = gp.Model()
    a = defaultdict(list)
    for r, route in enumerate(routes):
        for i in route:
            a[i].append(r)

    z = {}
    for r, route in enumerate(routes):
        z[r] = model.addVar(vtype=gp.GRB.BINARY, name=f'z_{r}')

    model.setObjective(gp.quicksum(cost(route, dists)*z[r] for r in z))

    for i in range(n_nodes):
        model.addConstr(
            gp.quicksum(z[r] for r in a[i]) >= 1
        )

    model.addConstr(gp.quicksum(z.values()) == n_vehicles)

    model._z = z
    return model

おなじデータセットで実験した結果は以下のとおりです。列生成法はあくまで近似解法である点に注意が必要です。
子問題を138回も解いており、以下に子問題を高速に解くかがポイントとなりそうです。

定式化	距離	距離/最適値	ルート皇甫嵩	求解時間
列生成法	1098	1.36	138	10min

まとめ

容量制約付き配送計画問題（Capacitated Vehicle Routing Problem, CVRP）の混合整数計画法を用いた定式化について整理しました。
特に制約や変数の数が非常に大きい場合について実験をおこない、定式化によって求解時間が大きく変わることがわかったかと思います。

脚注

基本的な CVRP では、1店舗の需要がトラックの容量を超えることはない様になっています ↩︎
$f_{ij}$ と $x_{ij}$ の関係式（つなぎ制約）を含めると $N^2$ ↩︎
$f^{k}_{ij}$ と $x_{ij}$ の関係式（つなぎ制約）を含めると $2N(N-1)$ ↩︎
$x_{ij}$ すら使っておらず、他の定式化と比べてかなり特殊な定式化方法ですが、変数の数や成約数を比較するため、敢えてここに記載しています ↩︎
2添字運搬車移動定式化とも呼ばれるようです ↩︎
より直接的に、 $S$ から外に向かうエッジ（上図の細い線、カットセットと呼びます）の本数がトラックの台数x2 となることを利用する方法もあります。 ↩︎
正確には $|S|=1$ の場合は次数制約と等価のため少し減ります ↩︎
あるいは Miller-Tucker-Zemlin 定式化と呼ばれるようです ↩︎
混合整数計画法では、最適解の上限と下限を狭めていくような計算を繰り返し、上限と下限が一致したら最適解であるとする方法が一般的ですが、この定式化では下限がなかなかあがらないことが知られています ↩︎
集合分割問題ではなく集合被覆問題と呼ばれる形式になります。 ↩︎
詳しくは「あたらしい数理最適化: Python言語とGurobiで解く」を御覧ください。 ↩︎
PuLP(cbc)では記述できないので、Gurobiのコードになっています ↩︎