Insight Edgeのデータサイエンティストのki_ieです。今日は最近勉強した数理最適化関連の技術を紹介します。

はじめに

本記事では、混合整数計画(MIP)問題に対する列生成法の発展的な使い方を紹介します。列生成法は巨大な線形計画(LP)問題を解くための厳密解法です。列生成法の有名な応用先として、多数の子問題と少数の子問題をまたぐ制約とに分解できるMIP問題のヒューリスティック解法があります。本記事では、このヒューリスティック解法を典型的なユースケースから少し拡張して利用する(具体的には、子問題に連続変数が登場するような場合でも利用できるように拡張する)方法を考えます。

本記事の構成は以下の通りです。

列生成法 : 列生成法について簡単に説明します。
列生成法のMIP問題への応用 - 典型的なケース : 列生成法のMIP問題への典型的な応用例を紹介します。
列生成法のMIP問題への応用 - 発展編 : 子問題に連続変数が登場するようなMIP問題での列生成法を紹介します。
プログラム実行例 : 簡単なサンプルプログラムを添付します。

列生成法

列生成法とは、大規模なLP問題を解くための技法です。この章では、列生成法の手続きと、その基本的な性質について簡単に説明します。

対象とするLP問題

次のLP問題 $(P)$ を考えます。

主問題 $(P)$
目的関数: max. $c^{⊤} x$
制約: $A x \leq b$ , $x \geq 0$

後ほど解析に利用するので、ここで $(P)$ の双対問題 $(D)$ も定義しておきましょう。

双対問題 $(D)$
目的関数: min. $b^{⊤} y$
制約: $A^{⊤} y \geq c, y \geq 0$

通常、LP問題はLPソルバーで高速に解くことができるのため求解に苦労することはないですが、 Aの列数があまりにも多く、Aの列を事前にすべて列挙しつくすことができない という場合状況が変わります。このような場合に用いられるのが列生成法です。

列生成法の手順

列生成法では「ほとんどの $A$ の列(と対応する変数)を無視して問題を解いても、最適解は得られるだろう」と山を張って、必要最小限の $A$ の列を列挙することで、効率的に $(P)$ を解こうとします。

具体的な手順を見てみましょう。列生成法では、最初に $A$ の列の一部 $A^{'}$ を列挙し、 $(P)$ の実行可能な部分問題 $(P^{'})$ を作ります。ここで、 $x^{'}$ は $(P)$ において $A^{'}$ の列に対応する変数のベクトルで、 $c^{'}$ は同様に対応する目的関数の重みです。( $(P^{'})$ が実行可能になるような $A^{'}$ の列挙は簡単にできると想定します。)

部分問題 $(P^{'})$
目的関数: max. $c^{' ⊤} x^{'}$
制約: $A^{'} x^{'} \leq b$ , $x^{'} \geq 0$

ここで $(P^{'})$ の最適解を $x_{*}^{'}$ とします。 $x_{*}^{'}$ は $(P)$ の実行可能解ですが、これがもし $(P)$ の最適解であればうれしいですね！また、もしそうでなければ、適当な列を $A^{'}$ に加えて部分問題の最適解の目的関数値を改善することにしましょう。

最適性の確認

まず $x_{*}^{'}$ の $(P)$ での最適性を確認しましょう。解析のために、 $(P^{'})$ の双対問題 $(D^{'})$ を次のように定義し、この最適解を $y_{*}^{'}$ とします。

双対問題 $(D^{'})$
目的関数: min. $b^{⊤} y$
制約: $A^{' ⊤} y \geq c^{'}, y \geq 0$

このとき、 $x_{*}^{'}, y_{*}^{'}$ の性質は以下のように整理できます。

$x_{*}^{'}$ は $(P)$ の実行可能解
$x_{*}^{'}$ の $(P)$ での目的関数値と、 $y_{*}^{'}$ の $(D)$ での目的関数値は一致する。
(というのも、 $x_{*}^{'}$ の $(P)$ での目的関数値 = $x_{*}^{'}$ の $(P^{'})$ での目的関数値 = $y_{*}^{'}$ の $(D^{'})$ での目的関数値 = $y_{*}^{'}$ の $(D)$ での目的関数値、であるから)

ここに「3. $y_{*}^{'}$ は $(D)$ の実行可能解」が加われば主実行可能解と双対実行可能解の目的関数値が一致していることになり、 $x_{*}^{'}$ は $(P)$ の最適解であることがわかりハッピーエンドです。 $y_{*}^{'}$ が $(D)$ の実行可能解であることを確認するためには、 $A^{⊤} y \geq c$ を確認すればOKです。ここで $A = [A^{'} | A^{″}], c = [c^{'} | c^{″}]$ と分解してこの条件を再度確認すると、 $[\begin{matrix} A^{' ⊤} \\ A^{″ ⊤} \end{matrix}] y \geq [\begin{matrix} c^{'} \\ c^{″} \end{matrix}]$ です。 $A^{' ⊤} y \geq c^{'}$ は $y_{*}^{'}$ が $(D^{'})$ の実行可能解であることから担保されているので、下段の $A^{″ ⊤} y \geq c^{″}$ が確認できればOKです。

何らかの方法で $A^{″ ⊤} y \geq c^{″}$ が証明できれば最適性の証明は終わりです。一般には、 $A$ の列をすべて列挙することが現実的でないという仮定の下ではこれを愚直に証明することは難しいでしょう。しかし、すこし唐突ですが 何らかの効率の良い方法により $c^{″} - A^{″ ⊤} y$ の最大値を求めることできれば (← まだ腑に落ちないで大丈夫です！)、最適性の証明が可能となります。つまり、 $c^{″} - A^{″ ⊤} y$ の最大値が $0$ 以下の場合、 $x_{*}^{'}$ が $(P)$ の最適解です。列生成法は、 $c^{″} - A^{″ ⊤} y$ の最大値を求められるような特殊な状況で用いられるアルゴリズムなのです。

列の生成

$c^{″} - A^{″ ⊤} y$ の最大値が $0$ より大の場合はどうでしょうか。この場合、列生成法のアルゴリズムはこの最大値を与える列を $A^{'}$ に加えて $(P^{'}), (D^{'})$ を解き、また最適性の確認に戻ります。 $A$ の列数は有限なのでこの手続きはいずれ停止します。実際には、 $A$ のごく一部の列を加えるだけでこの手続きが終了することが期待されます。

最小限の列生成法の解説はこの通りですが、

Aの列数があまりにも多く、Aの列を事前にすべて列挙しつくすことができない
何らかの効率の良い方法により $c^{″} - A^{″ ⊤} y$ の最大値を求めることができる

という条件を満たす具体的なケースが想像しづらく、そんなケースが実在するのか怪しく見えてしまいますね。これらの条件が実際に満たされるケースは列生成法のMIP問題への応用を見ると理解しやすいので、次の章で考えましょう。

列生成法のMIP問題への応用 - 典型的なケース

問題の定義

列生成法のMIP問題への応用を見ていきましょう。次の問題を考えます。

多数のタンクを部屋に収納したいと考えている。タンクにはそれぞれサイズが定まっている。
部屋は、同じサイズのものが複数与えられている。
与えられた部屋に収納可能な範囲で「タンクのサイズの総和」を最大化したい。

これは、数理最適化問題としては次のように記述されます。

タンク配置問題1 $(P_{1})$
パラメタ
$I$ : タンクの集合。 $| I | = m$
$J = {0, 1, \dots, n - 1}$ : 部屋の集合。 $| J | = n$
$s_{i} \in N_{> 0} (i \in I)$ : タンクのサイズ。
$S \in N_{> 0}$ : 各部屋の容量 (すべての部屋で共通とする)
変数
$x_{i j} \in {0, 1} (i \in I, j \in J)$ : タンク $i$ を部屋 $j$ に設置するか否か
目的関数
すべての部屋に置かれたタンクのサイズの和を最大化したい
max $\sum_{i \in I, j \in J} s_{i} x_{i j}$
制約
タンクは一度しか使えない
$\sum_{j \in J} x_{i j} \leq 1 (i \in I)$
部屋に置けるタンクのサイズの和には上限がある
$\sum_{i \in I} s_{i} x_{i j} \leq S (j \in J)$

$(P_{1})$ は、変数、制約の数という意味では効率よく表現されたMIP問題で「 $A$ の列数があまりにも多く、 $A$ の列を事前にすべて列挙しつくすことができない」ような問題ではありません。このままMIPソルバーに問題なく入力が可能です。しかし、このような多数の子問題の解を張り合わせて全体の解とするようなMIP問題では、上手に変形して列生成法を利用することで効率的に(理論保証のない)実行可能解を得られるケースがあります。

問題の変形

列生成法を利用できる形に問題を変形しましょう。まず、すべての部屋 $j \in J$ について、その部屋に入れられるタンクの組合せ全ての集合 $Q (j)$ を定義します。この集合は、問題規模に応じて指数的に大きくなります。これらの組合せ $q \in Q (j), j \in J$ について、部屋 $j$ の組合せ $q$ におけるタンクサイズの和 $s_{j q}^{'}$ とします。また部屋 $j$ の組合せ $q$ にタンク $i$ が含まれるかどうかを $ξ_{i j q}$ とします。これらを問題のパラメタだと思うと、部屋 $j$ で組合せ $q$ を採用するか否かの0/1変数 $X_{j q}$ を用いて、 $(P_{1})$ は次のように変形できます。

タンク配置問題1 - 変形版 $(M P_{1})$
パラメタ
$I$ : タンクの集合。 $| I | = m$
$J = {0, 1, . . ., n - 1}$ : 部屋の集合。 $| J | = n$
$Q (j) (j \in J)$ : 部屋 $j$ に入れられるタンクの組合せの集合
$s_{j q}^{'} (j \in J, q \in Q (j))$ : 部屋 $j$ の組合せ $q$ におけるタンクサイズ和
$ξ_{i j q} (j \in J, q \in Q (j))$ : 部屋 $j$ の組合せ $q$ にタンク $i$ が含まれるか。
変数
$X_{j q} \in {0, 1} (j \in J, q \in Q (j))$ : 部屋 $j$ で組合せ $q$ を採用するか
目的関数
すべての部屋に置かれたタンクのサイズの和を最大化したい
max $\sum_{j \in J, q \in Q (j)} s_{j q}^{'} X_{j q}$
制約
タンクは一度しか使えない
$\sum_{j \in J, q \in Q (j)} X_{j q} ξ_{i j q} \leq 1 (i \in I)$
組合せ選択
$\sum_{q \in Q (j)} X_{j q} \leq 1 (j \in J)$

列生成法の利用

問題 $(M P_{1})$ には $\sum_{j \in J} | Q (j) |$ 個の変数・列が登場するため、問題全体を記述しきってからこの問題をソルバーに解かせるというのは現実的ではありませんが、列生成法を用いて有用であろう少数の変数・列を列挙をすることで、良い実行可能界が算出できます列生成法を利用した $(M P_{1})$ の解き方は次のような手順です。背景には、 $(M P_{1}^{'})$ に登場する変数・列をだけを考慮しても、問題 $(M P_{1})$ の良い実行可能解が得られるだろうという期待があります。

列生成法を利用した $(M P_{1})$ の解き方:

$(M P_{1})$ の線形緩和問題 $(M L P_{1})$ の最適解を求める。ここで列生成法を利用する。
$(M P_{1})$ のうち前段の列生成法中に利用した変数・列だけに注目しが部分問題 $(M P_{1}^{'})$ を定義し、これをMIPソルバーで解く。

では列生成法を実行するために必要だった「 $c^{″} - A^{″ ⊤} y$ の最大値を求める」という部分はどのように実装すれば良いのでしょうか。 $(M L P_{1})$ の部分問題 $(M L P_{1}^{'})$ において $j \in J$ について列挙されている列の集合が $Q^{'} (j) (Q^{'} (j) \subseteq Q (j))$ であるとしましょう。このとき、「 $c^{″} - A^{″ ⊤} y$ の最大値を求める」というのは、以下の最適化問題を解くことにほかなりません。ここで $π_{i} (i \in I), μ_{j} (j \in J)$ はそれぞれ制約1, 2 に対応する $(M L P_{1}^{'})$ における双対最適解の値です。

子問題1 $(S P_{1})$
変数
$j \in J$
$q \in Q (j)$
目的関数
max. $s_{j q}^{'} - \sum_{i \in I} ξ_{i j q} π_{i} - μ_{j}$

この問題は、 $j$ を固定して $s_{j q}^{'}, ξ_{i j q}$ の定義に立ち返ると次の整数計画(IP)問題として記述できます。(補足: $(S P_{1})$ の定義では、本来は変数を $j \in J, q \in Q (j) ∖ Q^{'} (j)$ ととるのが $A^{″}$ の定義と一貫した方法です。しかし子問題の最適化結果が既知の列を返すことは列生成の終了判定条件が満たされる時だけなので、 $j \in J, q \in Q (j)$ で代用して問題ありません。またMIPで記述する際に $Q (j) ∖ Q^{'} (j)$ の取り扱いが難くなるので、このような変数の取り方を採用しています。)

子問題1' $(S P_{1}^{'} (j))$
変数
$x_{i} \in {0, 1} (i \in I)$ : タンク $i$ を組合せに入れるか否か。
目的関数
max $\sum_{i \in I} s_{i} x_{i} - \sum_{i \in I} π_{i} - μ_{j}$
制約
タンクは一度しか使えない
$\sum_{j \in J} x_{i} \leq 1 (i \in I)$
部屋に置けるタンクのサイズの和には上限がある
$\sum_{i \in I} s_{i} x_{i} \leq S$

$(S P_{1}^{'} (j))$ はナップサック問題なので、( $S$ が大きすぎなければ) 高速に解くことが可能です。これを利用して、列生成法の「 $c^{″} - A^{″ ⊤} y$ の最大値を求める」という手続きを実行できます。

問題 $(M P_{1})$ に対する列生成法を利用したヒューリスティックは、まとめると以下のような流れになります。

$(M P_{1})$ の線形緩和問題 $(M L P_{1})$ の最適解を求める。ここで列生成法を利用する。
- 初期解生成: 各部屋 $j$ に「一つもタンクを含まない組合せ」のみからなる集合を $Q^{'} (j)$ として与えればOK。一つもタンクを置かないというのは実行可能解であるため。
- 列の生成: $(S P_{1}^{'} (j))$ を利用して終了判定および列の追加を繰り返す。
$(M L P_{1})$ の記述に利用した変数・列だけを利用して、問題 $(M P_{1})$ の部分問題 $(M P_{1}^{'})$ をMIPソルバーで解く。

この問題に限らず、子問題が切り出せるMIP問題では、上記の流れのヒューリスティック解法が利用できます。具体的な問題によって精度の良し悪しがあったり、初期解生成に一工夫を加えたりという手間が必要なこともありますが、実装や個別の問題における数学的な整理がそこまで難しくないため、列生成法は子問題に分解可能な難しいMIP問題を解くためのツールとして頻繁に用いられています。

列生成法のMIP問題への応用 - 発展編

ここまで、子問題がIP問題であるようなMIP問題に列生成法が有効であることを確認しました。では、子問題の構造を持ちながら子問題が連続変数を含むMIP問題になっているケースに列生成法を適用するとき、どのような注意・変更が必要なのでしょうか。実際の問題を題材に考えていきましょう。

問題の定義

前述の問題を拡張した次の問題を考えます。

多数のタンクを部屋に収納し、それらのタンクに水を入れて保管したいと考えている。タンクにはそれぞれ自重とサイズが与えられている。
部屋は、同じサイズのものが $n$ 個与えられている。
与えられた部屋に収納可能な範囲で「タンクに入れる水の量」を最大化したい。
部屋は1次元的に等間隔で並んでおり、その座標は $0, 1, 2, \dots, n - 1$ である。
各部屋の重さは収納したタンクと水の重さの合計であり、建物の構造上、全体の重心が両端の部屋の中央(座標で $(n - 1) / 2$ )にある必要がある。

これは、数理最適化問題としては次のように記述されます。

タンク配置問題2 $(P_{2})$
パラメタ - $I$ : タンクの集合。 $| I | = m$ - $J = {0, 1, \dots, n - 1}$ : 部屋の集合。 $| J | = n$ - $s_{i} \in N_{> 0} (i \in I)$ : タンクのサイズ。 - $w_{i} \in R_{> 0} (i \in I)$ : タンクの重さ。 - $S \in N_{> 0}$ : 各部屋の容量 (すべての部屋で共通とする)
変数
$x_{i j} \in {0, 1} (i \in I, j \in J)$ : タンク $i$ を部屋 $j$ に設置するか否か
$y_{j} \geq 0 (j \in J)$ : 部屋 $j$ のタンクに合計どれだけ水を入れるか。
目的関数
タンクに注ぐ水の量の和を最大化したい
max $\sum_{j \in J} y_{j}$
制約
タンクは一度しか使えない
$\sum_{j \in J} x_{i j} \leq 1 (i \in I)$
部屋に置けるタンクのサイズの和には上限がある
$\sum_{i \in I} s_{i} x_{i j} \leq S (j \in J)$
部屋のタンクに入れる水の量の上限は、置かれたタンクのサイズの和によって決まる
$y_{j} \leq \sum_{i \in I} s_{i} x_{i j} (j \in J)$
重心が中央にあること
$\sum_{j \in J} W_{j} (j - c) = 0$
ただし $c = (n - 1) / 2$ , $W_{j} = y_{j} + \sum_{i \in I} w_{i} x_{i j}$

問題の変形

この問題も子問題の構造を持ちます。列生成法の適用を考えるため、子問題を列挙するような形に変形してみましょう。

タンク配置問題2 - 変形版1 $(M P_{2.1})$
パラメタ
$I$ : タンクの集合。 $| I | = m$
$J = {0, 1, . . ., n - 1}$ : 部屋の集合。 $| J | = n$
$Q (j) (j \in J)$ : 部屋 $j$ に入れられるタンクと水量の組合せの集合
$s_{j q}^{'} (j \in J, q \in Q (j))$ : 部屋 $j$ の組合せ $q$ におけるタンクサイズ和
$w_{j q}^{'} (j \in J, q \in Q (j))$ : 同、タンクの重さ。
$ξ_{i j q} (j \in J, q \in Q (j))$ : 部屋 $j$ の組合せ $q$ にタンク $i$ が含まれるか。
$λ_{j q} (j \in J, q \in Q (j))$ : 部屋 $j$ の組合せ $q$ で、タンクにどれだけの水をいれるか。
変数
$X_{j q} \in {0, 1} (j \in J, q \in Q (j))$ : 部屋 $j$ で組合せ $q$ を採用するか
目的関数
タンクに入った水の量の最大化
max $\sum_{j \in J, q \in Q (j)} λ_{j q} X_{j q}$
制約
タンクは一度しか使えない
$\sum_{j \in J, q \in Q (j)} X_{j q} ξ_{i j q} \leq 1 (i \in I)$
組合せ選択
$\sum_{q \in Q (j)} X_{j q} \leq 1 (j \in J)$
重心が中央にあること
$\sum_{j \in J} W_{j} (j - c) = 0$
ただし $c = (n - 1) / 2$ , $W_{j} = \sum_{q \in Q (j)} (w_{j q}^{'} + λ_{j q}) X_{j q}$

$(M P_{2.1})$ 正しい変形になっていますが、 $Q (j)$ が無限集合になっている点が気になります。また列生成法でそのごく一部 $Q^{'} (j) (j \in J)$ を列挙した後、本来連続変数で自由に(計算コスト低く)最適化できる $y_{j}$ が $λ_{j q}$ で(いくつかの飛び飛びの値に)固定されてしまっているのも気になりますね。今回は制約3.のように等式制約が入っているので、これは特に都合が悪いです。そこで、これらの課題を回避するために次のような定式化を考えます。

タンク配置問題2 - 変形版2 $(M P_{2.2})$
パラメタ
$I$ : タンクの集合。 $| I | = m$
$J = {0, 1, . . ., n - 1}$ : 部屋の集合。 $| J | = n$
$R (j) (j \in J)$ : 部屋 $j$ に入れられるタンクの組合せの集合
$s_{j r}^{'} (j \in J, r \in R (j))$ : 部屋 $j$ の組合せ $r$ におけるタンクサイズ和
$w_{j r}^{'} (j \in J, r \in R (j))$ : 同、タンクの重さ。
$ξ_{i j r} (j \in J, r \in R (j))$ : 部屋 $j$ の組合せ $r$ にタンク $i$ が含まれるか。
変数
$X_{j r} \in {0, 1} (j \in J, r \in R (j))$ : 部屋 $j$ で組合せ $r$ を採用するか
$Y_{j r} \in [0, 1] (j \in J, r \in R (j))$ : 部屋 $j$ で組合せ $r$ を採用したときに、何割の水を入れるか。
目的関数
タンクに入った水の量の最大化 max $\sum_{j \in J, r \in R (j)} s_{j r}^{'} Y_{j r}$
制約
タンクは一度しか使えない
$\sum_{j \in J, r \in R (j)} X_{j r} ξ_{i j r} \leq 1 (i \in I)$
組合せ選択
$\sum_{r \in R (j)} X_{j r} \leq 1 (j \in J)$
Yの制限
$Y_{j r} \leq X_{j r} (j \in J, r \in R (j))$
重心が中央にあること
$\sum_{j \in J} W_{j} (j - c) = 0$
ただし $c = (n - 1) / 2$ , $W_{j} = \sum_{r \in R (j)} w_{j r}^{'} X_{j r} + s_{j r}^{'} Y_{j r}$

この定式化では $R (j)$ は有限集合です。また、これの部分集合 $R^{'} (j) (j \in J)$ を列挙した後でも、連続変数の自由度は $Y_{j r}$ の形で保たれています。しかし今度は $R (j)$ の一要素 $r$ に対して、2つの列( $X_{i r}$ , $Y_{j r}$ )が対応してしまい、素朴に列生成を行うことができません。

しかし実は $(M P_{2.2})$ の定式化を用いると、列生成的なヒューリスティックを組み立てることができます。

列生成法の利用

$(M P_{2.2})$ を用いたヒューリスティックの手順は以下の通りです:

$(M P_{2.2})$ の線形緩和問題 $(M L P_{2.2})$ の最適解を求める。ここで列生成法を利用する。
- 初期解生成 :
  - $R^{'} (j)$ : 各部屋 $j$ に「一つもタンクを含まない組合せ」のみからなる集合を $R^{'} (j)$ として与えればOK。
- 終了判定と列生成
  - $R (j) (j \in J)$ の部分集合 $R^{'} (j)$ が列挙されているとする。これに対応する $(M L P_{2.2})$ の部分問題 $(M L P_{2.2}^{'})$ を解く。
  - 子問題 $(S P_{2}^{'})$ を利用して終了判定および列の追加を繰り返す。
    - 最適性が証明された場合、 2. に移行
    - 最適性が証明されなかった場合
      - $R^{'} (j)$ にタンクの組合せを追加する。
部分MIP問題の求解
- $(M P_{2.2})$ の $Q^{'} (j)$ に対応する部分問題 $(M P_{2.2}^{'})$ をMIPソルバーで解く

ここで、子問題 $(S P_{2}^{'})$ は典型的なケースで利用した論法をそのまま利用して、以下の最適化問題で定義します。ここで $π_{i} (i \in I), μ_{j} (j \in J), γ$ はそれぞれ制約1, 2, 3に対応する $(M L P_{2.1}^{'})$ の双対最適解における値です。

子問題2' $(S P_{2}^{'} (j))$
変数
$x_{i} \in {0, 1} (i \in I)$ : タンク $i$ を組合せに入れるか否か。
$y$ : どれだけの水を入れるか。
目的関数
max $\sum_{i \in I} y_{i} - \sum_{i \in I} π_{i} - μ_{j} - (w_{i} x_{i} + y) (j - c) γ$
制約
タンクは一度しか使えない
$\sum_{j \in J} x_{i} \leq 1 (i \in I)$
部屋に置けるタンクのサイズの和には上限がある
$\sum_{i \in I} s_{i} x_{i} \leq S$

有限回で手続きが終了するのか、 $j$ について $(S P_{2}^{'} (j))$ の目的関数値が 0 以下となったときに $(M L P_{2.2})$ が厳密に解けているのか、という部分の確認は少し複雑なので省略します。いずれも、 $(M P_{2.1})$ の $Q (j)$ の要素を子問題において連続変数が(整数変数を固定したときの)端点になっている解に絞った問題 $(M P V_{2})$ を考え、前述の手続きを $(M P V_{2})$ の線形緩和問題 $(M L P V_{2})$ (これは $(M L P_{2.2})$ と等価になります) に対する列生成法の手続きとして解釈することで、証明ができます。

この手続きの「終了判定と列生成」は、 $Q^{'} (j)$ に列を追加する際に「同じ整数変数値の列全て」を追加することで、有限回 (最大でも $\sum_{j \in J} | Q (j) |$ 回) で停止するようになっています。また、部分MIP問題の求解のステップでは $(M P_{2.2})$ のモデルを用いて子問題の連続変数を自由にコントロールして、等式制約に問題なく対処できる手続きになっています。

ここまで、子問題に連続変数が登場するようなMIP問題でも列生成法を使ったヒューリスティック解法の設計ができる場合がある、というご紹介でした。最後に $(P_{2})$ に対するナイーブな解法・列生成法的な解法のサンプルプログラムと実行結果を紹介して、この記事の結びとします。

プログラム実行例

準備

import numpy as np
import pulp  # MIPソルバーCBC へのラッパー https://github.com/coin-or/pulp
import itertools
from dataclasses import dataclass

@dataclass
class Params:  
    # (P_2) の問題パラメタをまとめるためのクラス。s や w の値はランダム生成するようにする。
    # 初期化に利用するメタパラメタ
    m: int = 5
    n: int = 10    
    S: int = 8
    
    # ランダムシード
    rand_seed: int = 0
    
    # 初期化後に生成するパラメタ
    I: range = None
    J: range = None
    s: np.ndarray = None
    w: np.ndarray = None
    
    def __post_init__(self):
        np.random.seed(self.rand_seed)
        self.I = range(self.m)
        self.J = range(self.n)
        self.s = np.random.randint(1, self.S, size=(self.m,))
        self.w = (1/2) * self.s

def rename_constraint(name: str, c: pulp.LpConstraint) -> pulp.LpConstraint:
    # pulp の制約に名前をつけてそのまま返す。デバッグ用。
    c.setName(name)
    return c

ナイーブな方法

def generate_problem(params: Params):    
    prob = pulp.LpProblem(sense=pulp.LpMaximize)
    
    ## 変数定義
    x = {
        (i, j): 
        pulp.LpVariable(name=f'x_{i}_{j}', cat=pulp.LpBinary)
        for i in params.I
        for j in params.J
    }
    y = {
        j: 
        pulp.LpVariable(name=f'y_{j}',lowBound=0, upBound=None, cat=pulp.LpContinuous)
        for j in params.J
    }
    
    ## 目的関数
    max_obj = pulp.lpSum(y.values())
    
    
    ## 補助的な線形関数
    # 線形関数を補助的に定義
    部屋内合計タンクサイズ = {
        j: pulp.lpSum(params.s[i] * x[i, j] for i in params.I)
        for j in params.J
    }
    
    ## 制約
    # タンクは一度しか使えない
    タンク利用_constr = {
        i: pulp.lpSum(x[i, j] for j in params.J) <= 1
        for i in params.I
    }
 
    # 部屋に置けるタンクのサイズの和には上限がある
    部屋内合計タンクサイズ_上限_constr = {
        j: 部屋内合計タンクサイズ[j] <= params.S
        for j in params.J
    }
 
    # 部屋のタンクに入れる水の量の上限は、置かれたタンクのサイズの和によって決まる
    部屋内水量_上限_constr = {
        j: y[j] <= 部屋内合計タンクサイズ[j]
        for j in params.J
    }
     
    # 部屋の重心が中心にあること
    c = (params.n-1)/2
    W = {j: pulp.lpSum(params.w[i] * x[i, j] for i in params.I) + y[j] for j in params.J}
    重心_constr = (pulp.lpSum( W[j] * (j-c) for j in params.J) == 0)
    

    # 目的関数・制約を最適化問題に追加
    prob.setObjective(max_obj)
    
    for c in itertools.chain(
        タンク利用_constr.values(),
        部屋内合計タンクサイズ_上限_constr.values(),
        部屋内水量_上限_constr.values(),
        [重心_constr]
    ):
        prob.addConstraint(c)
    
    return prob, x, y

my_params = Params(m=100, n=30, S=20, rand_seed=0)  # パラメタ生成
my_prob, *_ = generate_problem(my_params)  # 問題生成
# optimization_status = my_prob.solve(pulp.PULP_CBC_CMD(timeLimit=10, msg=1))  # 求解
optimization_status = my_prob.solve(pulp.PULP_CBC_CMD(timeLimit=10, msg=0))  # (ブログ整形用にmsg=0に設定)

結果は次の通り。目的関数値562.88の実行可能解が得られているが、上界値(Upper bound)は600とまだ最適化の余地がある。

Welcome to the CBC MILP Solver 
Version: 2.10.3 
Build Date: Dec 15 2019 

(中略)

Result - Stopped on time limit

Objective value:                562.87931034
Upper bound:                    600.000
Gap:                            -0.06
Enumerated nodes:               8030
Total iterations:               84038
Time (CPU seconds):             9.93
Time (Wallclock seconds):       10.00

Option for printingOptions changed from normal to all
Total time (CPU seconds):       9.94   (Wallclock seconds):       10.01

列生成法的な方法

def solve_colgen(params: Params, time_limit=10, max_loop=10, eps=0.0001):
    # タンクの組合せとその情報
    ## 変数定義 
    X = [[pulp.LpVariable(name=f'X_{j}_{0}', cat=pulp.LpBinary)] for j in params.J]  
    # 0/1 離散変数。 X[j][r] は j についてパターン r のタンクの置き方を採用することを意味。
    Y = [[pulp.LpVariable(name=f'Y_{j}_{0}', cat=pulp.LpContinuous, lowBound=0.0, upBound=1.0)] for j in params.J] 
    # 0-1 連続変数。 Y[j][r] は j についてパターン r のタンクに、容量に対してどれだけの割合の(0~1)水を入れるか。
    
    ## 生成された列(変数)のデータ
    ss = [[0.0] for j in params.J]  # ss[j][r] は部屋 j のパターン r のタンクの置き方における、水の容量。
    ww = [[0.0] for j in params.J]  # ww[j][r] は (同) タンク自体の重さ
    xi = {(i,j):[0.0] for j in params.J for i in params.I}  # xi[i,j][r] は (同) において、タンク i を利用するか否か (0/1)
    
    def solve_master_prob():
        master_prob =  pulp.LpProblem(sense=pulp.LpMaximize)

        # 目的関数
        max_obj = pulp.lpSum([
            ss[j][r] * Y[j][r]
            for j in params.J
            for r in range(len(X[j]))
        ])

        # 制約
        ## パターンの選択
        パターン選択_constr = {
            j: 
            rename_constraint(
                f'パターン選択_constr_{j}', 
                pulp.lpSum(X[j][r] for r in range(len(X[j]))) <= 1 
            )
            for j in params.J
        }
        
        ## Yの制限
        Yの制限_constr = {
            (j, r): 
            rename_constraint(
                f'Yの制限_constr_{j}_{r}',
                Y[j][r] <= X[j][r]
            )
            for j in params.J
            for r in range(len(X[j]))
        }        
        
        ## タンクは一度しか使えない
        タンク利用_constr = {
            i:
            rename_constraint(
                f'タンク利用_constr_{i}', 
                pulp.lpSum(
                    xi[i,j][r] * X[j][r] 
                    for j in params.J 
                    for r in range(len(X[j]))
                ) <= 1
            )
            for i in params.I
        }
        ## 部屋の重心が中心にあること
        c = (params.n-1)/2
        W = {
            j: 
            pulp.lpSum(ww[j][r] * X[j][r] + ss[j][r] * Y[j][r] for r in range(len(X[j])))
            for j in params.J
        }
        重心_constr = rename_constraint(
            '重心_constr',
            (pulp.lpSum( W[j] * (j-c) for j in params.J) == 0)
        )
       
    
        # 目的関数・制約を最適化問題に追加
        master_prob.setObjective(max_obj)
        
        for c in itertools.chain(
            パターン選択_constr.values(),
            Yの制限_constr.values(),
            タンク利用_constr.values(),
            [重心_constr],
        ):
            master_prob.addConstraint(c)                
        master_prob.solve(pulp.PULP_CBC_CMD(mip=0, msg=0))
        
        パターン選択_constr_dual = {r: v.pi for r,v in パターン選択_constr.items()}
        Yの制限_constr_dual = {r: v.pi for r,v in Yの制限_constr.items()}
        タンク利用_constr_dual = {r: v.pi for r,v in タンク利用_constr.items()}
        重心_constr_dual = 重心_constr.pi
        
        return master_prob, パターン選択_constr_dual, Yの制限_constr_dual, タンク利用_constr_dual, 重心_constr_dual,
    
    # 子問題の変数定義
    x = {
        (i, j): 
        pulp.LpVariable(name=f'x_{i}_{j}', cat=pulp.LpBinary)
        for i in params.I
        for j in params.J
    }

    y = {
        j: 
        pulp.LpVariable(name=f'y_{j}',lowBound=0, upBound=None, cat=pulp.LpContinuous)
        for j in params.J
    }
    
    def solve_sub_prob(パターン選択_constr_dual, Yの制限_constr_dual, タンク利用_constr_dual, 重心_constr_dual):
        # 陽に部屋ごとに問題を分離して解くことが自然だが、
        # 問題が明らかに分解可能な形なので、MIPソルバーが自動的に分解して求解することを期待して全ての j について
        # まとめた問題をソルバーに投げている。
        sub_prob = pulp.LpProblem(sense=pulp.LpMaximize, name='Sub')
        
        # 目的関数
        c = (params.n-1)/2
        W = {
            j: pulp.lpSum(params.w[i] * x[i, j] for i in params.I) + y[j] 
            for j in params.J
        }
        
        max_obj_部屋別 = {
            j: (
                y[j]  # 元の目的関数
                - pulp.lpSum(
                    タンク利用_constr_dual[i] * x[i,j]
                    for i in params.I 
                )  # タンク利用制約の双対変数由来
                - 重心_constr_dual * W[j] * (j-c)  # 重心制約の双対変数由来
                - パターン選択_constr_dual[j] # パターン選択制約の双対変数由来  TODO これでいいのか？
            )
            for j in params.J
        }
        max_obj = pulp.lpSum(max_obj_部屋別.values())
        
        # 線形関数を補助的に定義
        部屋内合計タンクサイズ = {
            j: pulp.lpSum(params.s[i] * x[i, j] for i in params.I)
            for j in params.J
        }
        
        # 部屋に置けるタンクのサイズの和には上限がある
        部屋内合計タンクサイズ_上限_constr = {
            j: rename_constraint(f'部屋内合計タンクサイズ_上限_constr_{j}', 部屋内合計タンクサイズ[j] <= params.S)
            for j in params.J
        }

        # 部屋のタンクに入れる水の量の上限は、置かれたタンクのサイズの和によって決まる
        部屋内水量_上限_constr = {
            j: rename_constraint(f'部屋内水量_上限_constr_{j}', y[j] == 部屋内合計タンクサイズ[j])
            # j: rename_constraint(f'部屋内水量_上限_constr_{j}', y[j] <= 部屋内合計タンクサイズ[j])
            # TODO DEV
            for j in params.J
        }    

        # 目的関数・制約を最適化問題に追加
        sub_prob.setObjective(max_obj)
        
        for c in itertools.chain(
            部屋内合計タンクサイズ_上限_constr.values(),
            部屋内水量_上限_constr.values(),
        ):
            sub_prob.addConstraint(c)

        sub_prob.solve(pulp.PULP_CBC_CMD(timeLimit=time_limit, msg=0))
        
        # マスター問題の update
        for j in params.J:
            if max_obj_部屋別[j].value() > eps:
                r_ub = len(X[j])
                X[j].append(
                    pulp.LpVariable(name=f'X_{j}_{r_ub}', cat=pulp.LpBinary)
                )
                Y[j].append(
                    pulp.LpVariable(name=f'Y_{j}_{r_ub}', cat=pulp.LpContinuous, lowBound=0.0, upBound=1.0)
                )
                ss[j].append(部屋内合計タンクサイズ[j].value())
                ww[j].append(pulp.lpSum(params.w[i] * x[i, j] for i in params.I).value())
                for i in params.I:
                    xi[i,j].append(x[i,j].value())
                
        return sub_prob, max_obj_部屋別      
    

    def get_naive_problem_variable_values():        
        _x = {
            (i, j): pulp.LpAffineExpression()
            for i in params.I
            for j in params.J
        }
        for i, j, in itertools.product(params.I, params.J):                
            for r in range(len(X[j])):
                _x[i, j] += (xi[i,j][r] * X[j][r])
        
        _y = {
            j: pulp.LpAffineExpression()
            for j in params.J
        }
        for j in params.J:
            for r in range(len(X[j])):
                _y[j] += (ss[j][r] * Y[j][r])

        _x_vals = {k: v.value() for k, v in _x.items()}
        _y_vals = {k: v.value() for k, v in _y.items()}
            
        return _x_vals, _y_vals, _x, _y
        
    
    for _i in range(max_loop):
        print(f'\niter: {_i} ----------------')
        # master_prob, タンク利用_constr_dual, 重心_constr_dual = solve_master_prob()
        master_prob, パターン選択_constr_dual, Yの制限_constr_dual, タンク利用_constr_dual, 重心_constr_dual = solve_master_prob()
        print(f'  MASTER_LP_STATUS:')        
        print(f'    Result: {pulp.LpStatus[master_prob.status]}')
        print(f'    MasterLP_Obj: {master_prob.objective.value()}')
        # print(f'    Duals:')
        # print(f'      パターン選択_constr_dual {list(パターン選択_constr_dual.values())}')
        # # print(f'      {Yの制限_constr_dual=}')
        # print(f'      タンク利用_constr_dual {list(タンク利用_constr_dual.values())}')
        # print(f'      重心_constr_dual {重心_constr_dual}')
        # print('***')
        # print(master_prob)
        # for j in params.J:
        #     for r in range(len(X[j])):
        #         if X[j][r].value() > 0.0:
        #             print(f'X[{j}][{r}] = {X[j][r].value()}')
        #         if Y[j][r].value() > 0.0:
        #             print(f'Y[{j}][{r}] = {Y[j][r].value()}')

        sub_prob, max_obj_部屋別 = solve_sub_prob(パターン選択_constr_dual, Yの制限_constr_dual, タンク利用_constr_dual, 重心_constr_dual)
        SubprobObj_max = max(e.value() for e in max_obj_部屋別.values())
        print(f'  SUBPROB_MIP_STATUS:')
        print(f'    Result: {pulp.LpStatus[sub_prob.status]}')
        print(f'    SubprobObj: {SubprobObj_max}')                
        # for j in params.J:
        #     print(f'      {j} : {max_obj_部屋別[j].value()}')        
        # print('***')
        # print(sub_prob)
        # for i in params.I:
        #     for j in params.J:
        #         if x[i,j].value() > 0.0:
        #             print(f'x[{i},{j}] = {x[i,j].value()}')
        # for j in params.J:
        #     if y[j].value() > 0.0:
        #         print(f'y[{j}] = {y[j].value()}')
        # 終了判定
        if SubprobObj_max <= eps:
            break
    
    # 生成された列を使ってMIP部分問題を解き、最終的な結果とする
    master_prob, パターン選択_constr_dual, Yの制限_constr_dual, タンク利用_constr_dual, 重心_constr_dual = solve_master_prob()        
    master_prob.solve(pulp.PULP_CBC_CMD(timeLimit=time_limit, msg=0))        
    print(f'MASTER_MIP_STATUS:')        
    print(f'  MasterMIP_Obj: {master_prob.objective.value()}')       
    
    # 解の検証
    x_vals, y_vals, *_ = get_naive_problem_variable_values()
    prob_for_validation, x_v, y_v = generate_problem(params)
    # ナイーブなモデルの変数を固定して、valid となるかを確認。
    for k, v in x_v.items():
        v.setInitialValue(x_vals[k])
    for k, v in y_v.items():
        v.setInitialValue(y_vals[k])
    print(f'  ValidationMIP_Obj: {prob_for_validation.objective.value()}')       
    print(f'  ValidationMIP_Valid: {prob_for_validation.valid(eps=eps)}')

solve_colgen(my_params, time_limit=5, max_loop=100)

iter: 0 ----------------
  MASTER_LP_STATUS:
    Result: Optimal
    MasterLP_Obj: None
  SUBPROB_MIP_STATUS:
    Result: Optimal
    SubprobObj: 20.0

iter: 1 ----------------
  MASTER_LP_STATUS:
    Result: Optimal
    MasterLP_Obj: 199.36507944
  SUBPROB_MIP_STATUS:
    Result: Optimal
    SubprobObj: 29.206349135000004

iter: 2 ----------------
  MASTER_LP_STATUS:
    Result: Optimal
    MasterLP_Obj: 260.0
  SUBPROB_MIP_STATUS:
    Result: Optimal
    SubprobObj: 20.0

iter: 3 ----------------
  MASTER_LP_STATUS:
    Result: Optimal
    MasterLP_Obj: 327.69230796
  SUBPROB_MIP_STATUS:
    Result: Optimal
    SubprobObj: 20.0

iter: 4 ----------------
  MASTER_LP_STATUS:
    Result: Optimal
    MasterLP_Obj: 408.96969695999996
  SUBPROB_MIP_STATUS:
    Result: Optimal
    SubprobObj: 23.0303030435

iter: 5 ----------------
  MASTER_LP_STATUS:
    Result: Optimal
    MasterLP_Obj: 456.8902041020001
  SUBPROB_MIP_STATUS:
    Result: Optimal
    SubprobObj: 20.0

iter: 6 ----------------
  MASTER_LP_STATUS:
    Result: Optimal
    MasterLP_Obj: 483.91531642000007
  SUBPROB_MIP_STATUS:
    Result: Optimal
    SubprobObj: 23.099094506

iter: 7 ----------------
  MASTER_LP_STATUS:
    Result: Optimal
    MasterLP_Obj: 502.4378701800001
  SUBPROB_MIP_STATUS:
    Result: Optimal
    SubprobObj: 20.0

iter: 8 ----------------
  MASTER_LP_STATUS:
    Result: Optimal
    MasterLP_Obj: 532.7289676199999
  SUBPROB_MIP_STATUS:
    Result: Optimal
    SubprobObj: 20.905528775

iter: 9 ----------------
  MASTER_LP_STATUS:
    Result: Optimal
    MasterLP_Obj: 535.1733378319999
  SUBPROB_MIP_STATUS:
    Result: Optimal
    SubprobObj: 20.0

iter: 10 ----------------
  MASTER_LP_STATUS:
    Result: Optimal
    MasterLP_Obj: 566.835042106
  SUBPROB_MIP_STATUS:
    Result: Optimal
    SubprobObj: 20.0

iter: 11 ----------------
  MASTER_LP_STATUS:
    Result: Optimal
    MasterLP_Obj: 589.4684529115998
  SUBPROB_MIP_STATUS:
    Result: Optimal
    SubprobObj: 20.0

iter: 12 ----------------
  MASTER_LP_STATUS:
    Result: Optimal
    MasterLP_Obj: 600.0000004290002
  SUBPROB_MIP_STATUS:
    Result: Optimal
    SubprobObj: 5.149999989795262e-07
MASTER_MIP_STATUS:
  MasterMIP_Obj: 578.9655172
  ValidationMIP_Obj: 578.9655172
  ValidationMIP_Valid: True

比較

ナイーブな方法での目的関数値が562.9、列生成法的なアプローチでは579.0ということで、一例だけの比較に意味はないですが少しだけ良い結果が得られています。

Insight Edge Tech Blog

子問題に連続変数が現れる問題における列生成法的アプローチ

はじめに

列生成法

対象とするLP問題

列生成法の手順

最適性の確認

列の生成

列生成法のMIP問題への応用 - 典型的なケース

問題の定義

問題の変形

列生成法の利用

列生成法のMIP問題への応用 - 発展編

問題の定義

問題の変形

列生成法の利用

プログラム実行例

準備

ナイーブな方法

列生成法的な方法

比較

参考文献