Dantzig-Wolfe分解と列生成法

この記事はJij Inc. Advent Calendar 2023の10日目の記事です。
株式会社Jij の松山です。

この記事では列生成法の話をしようと思います。
列生成法は、変数の規模が大きい問題を解くための手法の一つです。宮本先生の「はじめての列生成法」^[1]や梅谷先生の「資材切出し問題と列生成法」など良い説明がネットで見つかると思いますので、ここでは、具体例ではなく実務に役に立たない教科書的な部分を説明してみようかなと思います。
この記事の内容は主に二つの教科書^[2][3]を参考にしています。

Dantzig-Wolfe分解

ここでは、まず初めに次のようなブロック構造を内部に持つ整数計画問題を考えてみます。

$$max⁡c1⋅x1+c2⋅x2s.t.D1x1+D2x2≤d1B1x1≤b1B2x2≤b2x1∈Z+n1, x2∈Z+n2 \begin{align} &\max \quad \bm{c}_1\cdot\bm{x}_1 + \bm{c}_2\cdot\bm{x}_2 \\ &\mathrm{ s.t.}\quad D_1\bm{x}_1 + D_2\bm{x}_2\leq \bm{d}_1\\ &\hspace{8mm} B_1\bm{x}_1 \hspace{13mm} \leq \bm{b}_1\\ &\hspace{8mm} \hspace{13mm} B_2\bm{x}_2 \leq \bm{b}_2\\ &\notag\hspace{8mm}\bm{x}_1\in \Z^{n_1}_+, \ \bm{x}_2\in \Z^{n_2}_+ \end{align}$$

$x1\bm{x}_1$は$n1n_1$個の変数の組であり、$x2\bm{x}_2$は$n2n_2$個の変数の組であるとします。また、$Z+\Z_+$は正の整数の集合を表すとします。制約条件(2)は$mm$本の制約式、制約条件(3),(4)はそれぞれ$m1,m2m_1,m_2$本の制約式で書かれるとします。
この問題では、制約条件(3)では変数の組$x1\bm{x}_1$だけに依存し、制約条件(4)では変数の組$x2\bm{x}_2$だけに依存しています。しかし、制約条件(2)があるためにこの二つの変数の組が分離できなくなってしまっています。
ここでブロック構造とは、線形制約全体を一つの制約行列で書いたときに、(3),(4)の係数がブロック行列になっているからです。

このような構造の問題に対して、用いることのできる手法の一つがDantzig-Wolfe分解です。
Dantzig-Wolfe分解を使うことで、問題を元の問題よりも小さい子問題と元の問題よりも制約の少ない主問題へと分けることができます。
今回は二つのブロックがある問題ですが、複数あっても同じように扱うことができます。

まず初めに、制約条件(3),(4)で作られる実行可能領域$X1,X2X_1,X_2$をそれぞれ、

$$X1={x∈Rn1:B1x≤b1,x1∈Z+n1}, X2={x∈Rn2:B1x≤b1,x2∈Z+n2} X_1 = \left\{\bm{x} \in \R^{n_1} : B_1\bm{x} \leq \bm{b}_1, \bm{x}_1\in \Z^{n_1}_+\right\},\ X_2 = \left\{\bm{x} \in \R^{n_2} : B_1\bm{x} \leq \bm{b}_1,\bm{x}_2\in \Z^{n_2}_+\right\}$$

とかいておきます。すると、元の問題は、

$$max⁡∑j=12cj⋅xjs.t.∑j=12Djxj≤d1xj∈Xj∀j={1,2}\begin{align} &\max \quad \sum_{j=1}^2\bm{c}_j\cdot\bm{x}_j \\ &\mathrm{ s.t.}\quad \sum_{j=1}^2D_j\bm{x}_j \leq \bm{d}_1\\ &\notag\hspace{8mm} x_j \in X_j \quad \forall j =\{1,2\} \end{align}$$

と書くことができます。これは、元の問題を書き換えただけですが、これを見ると、制約条件(2)-(4)を同時に満たすものを見つけるのではなく、先に$X1X_1$,$X2X_2$に属する解を見つけておき、その中から制約条件(6)を満たし、目的関数を最大化するような解を見つければ良いのではないかと考えることができます。

次に、$XjX_j$の凸包$conv(Xj)\mathrm{conv}(X_j)$を考えます。凸包の内部の点$x\bm{x}$は、凸包$conv(Xj)\mathrm{conv}(X_j)$の端点 ${vk}k∈K\{v_k\}_{k\in K}$の凸結合と端射線 ${rh}h∈H\{r_h\}_{h\in H}$の錐結合

$$conv(Xj)={xj∈Rnj:xj=∑k∈Kλjkvjk+∑h∈Hμjhrjh, ∑k∈Kλjk=1, λjk∈R+, μjh∈R+}\mathrm{conv}(X_j) =\left\{\bm{x}_j \in \R^{n_j} : \bm{x}_j = \sum_{k\in K} \lambda_{jk} \bm{v}_{jk} + \sum_{h\in H} \mu_{jh} \bm{r}_{jh},\ \sum_{k\in K} \lambda_{jk} = 1,\ \lambda_{jk}\in \R_+,\ \mu_{jh}\in \R_+ \right\}$$

で書くことができます。これは、通常の線形計画で出てくる多面体を制約条件という内部の点が満たすべき条件で表現したものではなく、多面体の端点や端射線の線型結合で表現したものになっています。
端点や端射線は聞き馴染みのないものかもしれませんが、大まかに多面体の頂点が端点であり、多面体の非有界方向へのへりが端射線に対応します。ここで非有界方向が出てくるのは、一般の線形計画問題の実行可能解が必ずしも有界であるとは限らないためです。
注意としては、ここで表現される多面体は、整数条件も考えている、つまり格子点に対する凸包を考えているので、一般には線形計画$Bix≤biB_i\bm{x} \leq \bm{b}_i$で表現される凸包とは異なります。
これを元の問題に代入してあげると問題を次のように書き直すことができます。

$$max⁡∑j=12(∑k∈Kλjkcj⋅vjk+∑h∈Hμjhcj⋅rjh)s.t.∑j=12(∑k∈KDjλjkvjk+∑h∈HDjμjhrjh)≤d1∑k∈Kλjk=1∀j={1,2}∑k∈K∑k∈Kλjkvjk+∑h∈Hμjhrjh∈Z∀j={1,2}λjk≥0, μjk≥0 \begin{align} &\max \quad \sum_{j=1}^2 \left( \sum_{k\in K} \lambda_{jk} \bm{c}_j\cdot\bm{v}_{jk} + \sum_{h\in H} \mu_{jh} \bm{c}_j\cdot\bm{r}_{jh} \right) \\ &\mathrm{ s.t.}\quad \sum_{j=1}^2\left(\sum_{k\in K} D_j\lambda_{jk} \bm{v}_{jk} + \sum_{h\in H} D_j\mu_{jh} \bm{r}_{jh}\right) \leq \bm{d}_1\\ &\hspace{8mm}\sum_{k\in K} \lambda_{jk} = 1\quad \forall j =\{1,2\}\\ &\hspace{8mm}\sum_{k\in K} \sum_{k\in K} \lambda_{jk} \bm{v}_{jk} + \sum_{h\in H} \mu_{jh} \bm{r}_{jh}\in \Z\quad \forall j=\{1,2\}\\ &\notag\hspace{8mm} \lambda_{jk}\geq 0,\ \mu_{jk} \geq 0\\ \end{align}$$

このように元の問題に対して、$XjX_j$が作る凸包の点を代入して得られたこの問題をDantzig-Wolfe Reformulationと呼びます。制約条件(10)は$conv(Xj)\mathrm{conv}(X_j)$の内部の点が整数条件を満たすようにするために加えており、元の問題の整数条件に対応します。この問題は、$XjX_j$それぞれを解いて得られた解をいかにうまく組み合わせて、制約条件を満たす最適解を見つけるかという問題になっています。

さて、ここまでの流れを整理しておくと、まず初めに、$conv(Xj)\mathrm{conv}(X_j)$の全ての端点と端射線を列挙し、その次に(7)-(10)で記述される問題を解くという流れになっています。このようなやり方には明らかに無理があります。なぜなら、端点と端射線を全て列挙するというのはあまり現実的ではありませんし、そもそも全てを列挙する労力で元の問題を解けてしまえそうです。

ここから次のステップに進めるアイディアは、全てを列挙するのではなく考慮する価値がありそうな解だけを列挙するということです。そして、考慮する価値がありそうな解というものを効率よく探そうというのが、列生成法です。

さて、列生成法の下準備として、(7)-(10)から$xx$が整数である条件に対応する制約条件(10)を抜いてしまい、元の問題の連続緩和問題を考えます。

$$max⁡∑j=12(∑k∈Kλjkcj⋅vjk+∑h∈Hμjhcj⋅rjh)s.t.∑j=12(∑k∈KDjλjkvjk+∑h∈HDjμjhrjh)≤d1∑k∈Kλjk=1∀j={1,2}λjk≥0, μjk≥0 \begin{align} &\max \quad \sum_{j=1}^2 \left( \sum_{k\in K} \lambda_{jk} \bm{c}_j\cdot\bm{v}_{jk} + \sum_{h\in H} \mu_{jh} \bm{c}_j\cdot\bm{r}_{jh} \right) \\ &\mathrm{ s.t.}\quad \sum_{j=1}^2\left(\sum_{k\in K} D_j\lambda_{jk} \bm{v}_{jk} + \sum_{h\in H} D_j\mu_{jh} \bm{r}_{jh}\right) \leq \bm{d}_1\\ &\hspace{8mm}\sum_{k\in K} \lambda_{jk} = 1\quad \forall j =\{1,2\}\\ &\notag\hspace{8mm} \lambda_{jk}\geq 0,\ \mu_{jk} \geq 0\\ \end{align}$$

このような緩和問題をDantzig-Wolfe Relaxationと呼びます。
この問題をここではLinearized Master Problem(LMP)と呼ぶことにします。

列生成法

先ほども述べたように列生成法の基本的なアイディアは全てを最初から列挙するのではなく、最初は少ない数だけ列挙しておいて、徐々に必要な解を追加していくという方法です。
列生成法の基本的なアイディアは単体法から来ています。単体法には、reduced cost $c~\tilde{c}$という、非基底変数を増加させたときの目的関数の改善量を表す量が存在します。単体法では、このreduced costが正のものに対応する変数を増加させることで、目的関数を改善することができます。列生成法では、最大のreduced costを与える解を追加する価値がある解であると考え、最大のreduced costを与える解を見つける問題を子問題として解いていくことで、解を追加していきます。

より詳細にこれについてみていきたいと思います。
まず、列挙する数を制限した$K′⊂KK'\subset K$, $H′⊂HH'\subset H$というsubsetを考え、考慮するものをこのsubsetだけに制限した問題を考えます。

$$max⁡∑j=12(∑k∈K′λjkcj⋅vjk+∑h∈H′μjhcj⋅rjh)s.t.∑j=12(∑k∈K′Djλjkvjk+∑h∈H′Djμjhrjh)≤d1∑k∈K′λjk=1∀j={1,2}λjk≥0, μjk≥0 \begin{align} &\max \quad \sum_{j=1}^2 \left( \sum_{k\in K'} \lambda_{jk} \bm{c}_j\cdot\bm{v}_{jk} + \sum_{h\in H'} \mu_{jh} \bm{c}_j\cdot\bm{r}_{jh} \right) \\ &\mathrm{ s.t.}\quad \sum_{j=1}^2\left(\sum_{k\in K'} D_j\lambda_{jk} \bm{v}_{jk} + \sum_{h\in H'} D_j\mu_{jh} \bm{r}_{jh}\right) \leq \bm{d}_1\\ &\hspace{8mm}\sum_{k\in K'} \lambda_{jk} = 1\quad \forall j =\{1,2\}\\ &\notag\hspace{8mm} \lambda_{jk}\geq 0,\ \mu_{jk} \geq 0\\ \end{align}$$

この問題をRestricted Mater Problem(RMP)と呼びます。制約条件(15)に対応する双対変数を$πi\pi_i$, 制約条件(16)に対応する双対変数を$σj\sigma_j$とします。このRMPの実行可能解は、必ずLMPの実行可能解になります。なぜならば、 $λjk=0, (k∈K∖K′), μjh=0, h∈H∖H′\lambda_{jk} = 0,\ (k \in K\setminus K'),\ \mu_{jh} = 0,\ h \in H\setminus H'$とおくことで、LMPの実行可能解にすることができるからです。

さて、この線形計画問題を解いて、最適解$(λ∗,μ∗)(\bm{\lambda}^*, \bm{\mu}^*)$及び、双対変数 $(π∗,σ∗)(\bm{\pi}^*, \bm{\sigma}^*)$を得たとします。
これを活用して、reduced costを求める子問題を考えます。ここで、各子問題で求めるreduced costは

$$ζj=max⁡xj∈Xj(cj−π⋅D)⋅xk−σj,j={1,2}\begin{align}\zeta_j = \max_{\bm{x}_j\in X_j} (\bm{c}_j - \bm{\pi}\cdot D)\cdot \bm{x}_k- \sigma_j,\quad j = \{1,2\}\end{align}$$

この子問題をpricing problemと呼ぶこともあります。この問題はそれぞれ独立しているので、それぞれで解くことができます。この子問題を解き、追加すべき解を見つけていきます。
さて、どのように追加すべき解を選択するかをみてみると、

1. もし、(17)が非有界であれば、$(cj−π⋅D)⋅rh>0(\bm{c}_j - \bm{\pi}\cdot D)\cdot \bm{r}_h>0$が存在したことになるので、この端射線を$H′H'$に追加します。
2. もし、(17)が非有界ではなく、$ζj>0\zeta_j>0$であれば、$(cj−π⋅D)⋅vk−σj>0(\bm{c}_j - \bm{\pi}\cdot D)\cdot \bm{v}_k- \sigma_j > 0$となる端点が存在したことになるので、この端点を$K′K'$に追加します。
3. もし、(17)が非有界ではなく、$ζj≤0\zeta_j\leq 0$であれば、追加すべきものがないことになるので、ここでアルゴリズムを停止します。
   となります。

以上が基本的な列生成法の流れになります。
なので、列生成法の流れとしては、

1. RMPを解くことで、双対変数 $(π∗,σ∗)(\bm{\pi}^*, \bm{\sigma}^*)$を得る。
2. 最大reduced costを求める子問題$ζj=max⁡xj∈Xj(cj−π⋅D)⋅xk−σj\zeta_j = \max_{\bm{x}_j\in X_j} (\bm{c}_j - \bm{\pi}\cdot D)\cdot \bm{x}_k- \sigma_j$を解き、新たに集合$K′,H′K',H'$に追加すべき解を見つける。
3. $ζj≤0\zeta_j\leq 0$であれば、アルゴリズムを停止し、そうでなければ、解をRMPに追加し、1に戻る。
   となります。

結局何を解いたのか

さて、ここまで長い道のりを経て列生成法を説明してきました。
では、この列生成法を用いることで、元々解きたかった問題(1)-(4)は解けたのでしょうか？お気づきの方もいると思いますが、解けていません。結局、この道のりで解けたのは、元の問題を線形緩和したLMPだけです。そのため、必ずしも、ここで得た解が整数性を満たしているとは限りません。
ですが、適当に丸めるとか、最終的に得られた$K′,H′K',H'$を使って(7)-(10)を解き直すなどを行うことで、問題によっては比較的良い解を得ることができます。また、厳密解が欲しい場合には、列生成法と分枝限定法を組み合わせた分枝価格法を用いることもできます。

最後に

＼Rustエンジニア・数理最適化エンジニア募集中！／
株式会社Jijでは、数学や物理学のバックグラウンドを活かし、量子計算と数理最適化のフロンティアで活躍するRustエンジニア、数理最適化エンジニアを募集しています！
詳細は下記のリンクからご覧ください。皆さんのご応募をお待ちしております！
Rustエンジニア: https://open.talentio.com/r/1/c/j-ij.com/pages/51062
数理最適化エンジニア: https://open.talentio.com/r/1/c/j-ij.com/pages/75132

脚注

1. 宮本 裕一郎, はじめての列生成法, オペレーションズ・リサーチ, 57 (2012), 198-204.

2. L. A. Wolsey, Integer programming. John Wiley & Sons, 2020.

3. M. Conforti et al. Integer programming. Springer International Publishing, 2014.