2. ナップサック問題と列生成法 #

2.1. 背景 #

ナップサック問題には次の側面がある．

数理最適化問題として線形に定式化をする場合の典型的な考え方を提供してくれる．
数式が素直で基本的であるため，別の問題の文脈の中で，ナップサック問題としての構造を見出しやすい．

後者について，ここでは列生成法との関連を述べる．

2.1.1. 列生成法と連続緩和問題 #

列生成法そのものは線形最適化問題を直接に解くというよりは，連続緩和の文脈でよく用いられる．整数最適化問題 \(P_{\mathbb{Z}}\) のように，応用範囲が広く，しかし問題規模が大きくなりやすい場合が典型的である．この場合，問題 \(P_{\mathbb{Z}}\) を連続緩和した問題 \(P\) を考え，列生成法を適用してなるべく少ない変数 (または制約条件式) で求解することになる．

これから説明することは問題の定式化を違った側面からアプローチすることもあって，抽象性が高くそして多分に技巧的な面があろう．それ故に一つずつ整理して述べたい．

2.2. 設定 #

列生成法の典型的な応用例として，ビンパッキング問題を採り上げよう．問題設定が単純であるため，列生成法のアルゴリズムに集中でき，教育的でもあると考えられる．

2.2.1. ビンパッキング問題 #

2.2.1.1. straight forward な定式化 #

ビンパッキング問題は次のように容易に定式化できる．

集合・要素

アイテム \(i\) とアイテム集合 \(I\)

(1)#\[i\in I\]
ビン \(b\) とビン集合 \(B\)

(2)#\[b\in B\]

定数

アイテム \(i\) の容量

(3)#\[IVOL_i\]
ビン \(b\) の共通する容量

(4)#\[BVOL\]

変数

アイテム \(i\) をビン \(b\) に入れるかどうか．

(5)#\[z_{i,b} \in \mathbb{B}\]
ビン \(b\) を使用するかどうか．

(6)#\[use_b \in \mathbb{B}\]

制約条件

アイテム \(i\) は必ず何れかのビンに入れる．

(7)#\[\sum_{b\in B} z_{i,b} = 1, \, \forall i\in I\]
ビンへあふれずに入れる．

(8)#\[\sum_{i\in I} IVOL_i\cdot z_{i,b} \leq BVOL\cdot use_b, \, \forall b\in B\]

目的関数

ビンの使用数を最小化

(9)#\[\mathrm{Minimize:}~ \sum_{b\in B} use_b\]

以上の定式化は straight forward な定式化である．幾つかのコメントをしよう．

アイテム一つにつきビンを一つ割り当てることを考えれば，必要なビンの総数の上限値は，アイテムの総数 \(N=|I|\) に等しい．

(10)#\[|B| = |I| = N\]
変数の総数のオーダーは \(O(|B||I|+|I|)=O(N^2)\) である．アイテムが \(10^2\) 個あれば，\(O(10^4)\) となってそれなりに大きな規模になる．
\(|B|=|I|\) という見積もりは最悪の場合を想定したものであり，実際のアイテム集合がビンの容量にほぼ等しいものだけで構成されていることは特殊なケースである．つまり変数が多過ぎる可能性がある．殆どが \(0\) 値をとる active でない変数と予想される．
制約条件式の総数のオーダーは \(O(|I| + |B|)=O(N)\) であり，線形にスケールする．

これらのことから，ビンパッキング問題を上記の定式化で直接に求解するのは，比較的小さな規模が現実的であり，もう少し大きな規模となると耐えられないと考えられる．一方で無駄な変数が多そうであるから，列生成法が有用ではないかと考えられる．

2.2.1.2. パターンという観点でみた定式化 #

straight forward な定式化を忘れて，違うアプローチによる定式化を試みよう．それはパターンという観点を中心に問題の離散構造を覆いつくす定式化である．

ビンパッキング問題を今一度思い出そう．ビンにアイテムを入れる方法にはいろいろな方法がある．これらをすべて列挙するとき，その中の何れか一つ以上の方法が最適解を与えていることは (実行可能である限り) 間違いない．パターンとはここで列挙するところの複数のビンにアイテムを入れる方法のことである．

例えば次を満たすビン集合 \(B\) とアイテム集合 \(I\) があったとしよう．

(11)#\[\begin{split}BVOL &= 5, \\ I &= \{1,2,3\}, \\ (IVOL_1, IVOL_2, IVOL_3) &= (1, 2, 4)\end{split}\]

このときアイテムを選び出すパターンをすべて列挙すると，次のパターンだけある．

パターン番号	パターン	合計容量	一つのビンに入るかどうか
\(0\)	\(\emptyset\)	\(0\)	yes
\(1\)	\(\{1\}\)	\(1\)	yes
\(2\)	\(\{2\}\)	\(2\)	yes
\(3\)	\(\{3\}\)	\(4\)	yes
\(4\)	\(\{1,2\}\)	\(3\)	yes
\(5\)	\(\{1,3\}\)	\(5\)	yes
\(6\)	\(\{2,3\}\)	\(6\)	no
\(7\)	\(\{1,2,3\}\)	\(7\)	no

このうち，興味あるパターンはアイテムを使用し，与えられたビンの容量に収まるパターンである．つまりパターン番号 \(1,2,3,4,5\) である．今回はビンの容量と各アイテムの容量から，一部のパターンのみが許容されたが，それら容量の値によってはパターン番号 \(0\) を除いてすべて許容されうる．そして「パターン一つにつき，一つのビンが対応している」ということに注意しよう．

これらのことからビンパッキング問題について次の定式化を考え出すことができる．但し straight forward な定式化で既に定義したものは，繰り返しとなるため再掲していない．

集合・要素

ビン一つに収まるパターンのパターン番号 \(p\) とそれらからなるパターン番号集合 \(PATTERNS\)

(12)#\[p\in PATTERNS\]

定数

パターン番号 \(p\) のパターンがアイテム \(i\) を含んでいるかどうか，というパターン番号とアイテムの関連付け情報 \(ASSOC\)

(13)#\[ASSOC_{p,i}\]

変数

パターン番号 \(p\) を採用するかどうか．

(14)#\[x_p \in \mathbb{B}\]

制約条件

アイテム \(i\) は必ず何れか一つ以上のパターン番号 \(p\) のパターンで選ぶ．

(15)#\[\sum_{p\in PATTERNS} ASSOC_{p,i} x_p \geq 1, \, \forall i\in I\]

目的関数

パターン数の最小化 (ビンの使用数を最小化)

(16)#\[\mathrm{Minimize:}~ \sum_{p\in PATTERNS} x_p\]

注釈

制約条件式 (15) より，ある一つのアイテムが一つしかないにもかかわらず，解によっては複数のビンに入れたパターンが選び出されうる定式化になっている．この場合でもビンの総数は正しくカウントされているため，後処理でアイテムの重複さえ除去できれば，どのビンにどのアイテムを入れるかという解を得ることができる．この後処理はビンパッキング問題よりも難しいとも考えられないし，また運用によって許容される解かもしれない．定式化にはこの種の冗長性がしばしば伴う．

上記の定式化についてもオーダーを見積もろう．パターンの数はアイテムの総数 \(N=|I|\) に対して，パターンの総数のオーダーは \(O(2^N)\) と指数的である．

(17)#\[\sum_{l=1}^N {N\choose l} = \sum_{l=0}^N {N\choose l} - {N\choose 0} = (1+1)^N - 1 = 2^N - 1\]

つまり変数の総数のオーダーは指数的にスケールする．一方で制約条件式の総数のオーダーは \(O(|I|)=O(N)\) と線形にスケールする．今回の定式化もまた変数の総数が多い．また straight forward の定式化と同様に，無駄なパターンが数多くあるとも考えられる．よって再び列生成法が有用ではないかと考えられるわけである ¹ ¹変数の数が多い場合には疎な記述も検討すべきであるが，今回のビンパッキング問題では任意の入力データに対して，疎な構造を見出すのは難しい．．

2.3. パターンと列生成法 #

パターンという観点は列生成法と相性が良い．というのも，列生成法では次の価格付け問題を解く必要があったが，そこにナップサック問題としての構造を見出すことができるからである．

(18)#\[\min_{j\in J\setminus J_k} (\vec{\gamma})_j = \min_{j\in J\setminus J_k} \left((\vec{c})_j - \sum_{i\in I}A_{ji}(\vec{\omega}_k^*)_i\right)\]

パターンを用いたビンパッキング問題の定式化の場合，考える主問題の標準形 \(P(\vec{b},A,\vec{c};I,J)\) および上記の問題の諸量は次のように対応する．但し \(PATTERNS_k(\subset PATTERNS)\) は列生成法での \(k\) 反復目のパターン番号集合である．

(19)#\[\begin{split}i \in I &\to i\in I, \\ (\vec{b})_i &\to 1, \\ j\in J &\to p\in PATTERNS, \\ J_k &\to PATTERNS_k, \\ (\vec{c})_j &\to 1, \\ A_{ji} &\to ASSOC_{p,i}\end{split}\]

よって考えるべき価格付け問題はパターン番号 \(p\) を通した \(ASSOC_{p,i}\) の自由度のみが変量とわかる．これから目的関数については最小化問題から次の最大化問題として焼き直すことができる．

(20)#\[\max_{p\in PATTERNS\setminus PATTERNS_k} \sum_{i\in I}ASSOC_{p,i}(\vec{\omega}_k^*)_i\]

ここまで整理したが，もう少し \(ASSOC_{p,i}\) について着目しよう． \(ASSOC_{p,i}\) とは「ビンに収まるパターン集合の中で，アイテム \(i\) がパターン番号 \(p\) のパターンに属しているかどうか」を表した．今，これを変数として扱いたいので，このために次の変数を定義しよう．

変数: アイテム \(i\) を使用するかどうか．

(21)#\[y_i\in \mathbb{B}\]

これはパターンとはアイテム集合の各アイテムを選び出すかということと等価であることを用いている．

(22)#\[\mathrm{パターン番号}~p \leftrightarrow \mathrm{パターン} \leftrightarrow \{y_i\}_{i\in I}\]

例えば (11) のインスタンスの場合，パターンは \((y_1,y_2,y_3)\) の \(0,1\) の並びと下表のように一対一対応する．

パターン番号	パターン	\((y_1,y_2,y_3)\)
\(0\)	\(\emptyset\)	\((0,0,0)\)
\(1\)	\(\{1\}\)	\((1,0,0)\)
\(2\)	\(\{2\}\)	\((0,1,0)\)
\(3\)	\(\{3\}\)	\((0,0,1)\)
\(4\)	\(\{1,2\}\)	\((1,1,0)\)
\(5\)	\(\{1,3\}\)	\((1,0,1)\)
\(6\)	\(\{2,3\}\)	\((0,1,1)\)
\(7\)	\(\{1,2,3\}\)	\((1,1,1)\)

これら対応の中で問題で扱うパターン番号集合 \(P\) は \(P=\{1,2,3,4,5\}\) であったが，このために必要な \(\{y_i\}_{i=1}^3\) は \(IVOL_1\cdot y_1 + IVOL_1\cdot y_1 + IVOL_1\cdot y_3 \leq BVOL=5\) を満たす \(\{y_i\}_{i=1}^3\) のことである．これは数理最適化問題の文脈では変数に課せられる制約条件に他ならない．

以上から，故に次のナップサック問題に価格付け問題を変形できたとわかる．

(23)#\[\begin{split}& \mathrm{Maximize:}~ \sum_{i\in I}(\vec{\omega}_k^*)_i\cdot y_i \\ & \mathrm{subject~to:}~ \\ & \sum_{i\in I} IVOL_i\cdot y_i \leq BVOL\end{split}\]

2.4. 列生成法の実装 #

PySIMPLE でビンパッキング問題を列生成法によって解く実装をしよう．列生成法が汎用的な枠組みでもあることから，共通化可能な部分についても確認しておこう．

2.4.1. 設計概要 #

振り返りのために，列生成法の手続きを次にまとめておこう．

\(k=1\) とする．
主問題 \(P\) のうち，適当な数の変数を選び出し，\(J_k\) を定義する．
限定主問題 \(P|_k\) を求解して，実行可能解 \(\vec{x}_k^{\star}\) を構成する．
限定双対問題 \(D|_k\) を求解して，最適解 \(\vec{\omega}_k^*\) を求める．
最適解 \(\vec{\omega}_k^*\) に対して価格付け問題に等価なナップサック問題 (23) を解いて，主問題 \(P\) の被約費用が非負かどうか決定する．非負ならば \(\vec{x}_k^{\star}\) を最適解とする．負ならば \(k\to k+1\) および次式で \(J_k\) を更新する．手続き 3 へ．

(24)#\[\begin{split}\bar{J} &= \operatorname*{argmin}_{j\in J} (\vec{\gamma})_j, \\ J_k &\leftarrow J_k \cup \bar{J}\end{split}\]

列生成法は汎用的な枠組みであるが，そのうち共通しない部分は次の設計である．

パターン番号集合の初期集合 \(PATTERNS_1\) の生成
価格付け問題の解き方

場合によっては他の部分も問題に応じた実装が必要になろう．

今回，パターン番号集合の初期集合 \(PATTERNS_1\) としては，一つのアイテムを一つのビンに割り当てる次の集合を考えよう．

(25)#\[PATTERNS_1 = \bigcup_{i\in I} \{i\}\]

上記を除いた部分については共通化が可能であるから，列生成法の処理を次のスーパークラスとして実装をまとめることにしよう．

2.4.2. スーパークラスの実装例 #

PySIMPLE による列生成法を司るスーパークラスの実装例 (column_generator.py) として次を挙げる．

PySIMPLE Sample Code (column_generator.py)

# -*- coding: utf-8 -*-
from abc import ABCMeta, abstractmethod
from pysimple import *

class ColumnGenerator(metaclass=ABCMeta):
    iternum = property(lambda self: self.__iternum)
    patternnum = property(lambda self: self.__pn0 + self.__iternum)

    @abstractmethod
    def create_init_pattern(self, *, pattern, b=1, c=1):
        """set self.pattern, self.b, self.c"""
        self.pattern = pattern
        self.b = b
        self.c = c
        self.__iternum = 0
        self.__pn0 = len(pattern)

    # type は solve_pricing_problem の変数 type と同じ
    def solve_primal_problem_for_selection_of_pattern(self, *, type=bin, silent=True):  # 主問題
        """生成されたパターンから b を満たす組合せを選択する"""
        pi = Element(value=self.pattern)
        p = Element(set=pi(0).set)
        b = Parameter(index=pi(1), value=self.b)
        c = Parameter(index=pi(0), value=self.c)
        ASSOC = Parameter(index=pi, value=self.pattern)  # パターン p にアイテム i を含むか (0-1)．
        x = Variable(index=p, type=type, lb=0)  # パターン p を選択するか (type=bin)，幾つ選択するか (type=int)．

        prb = Problem()
        prb += Sum(ASSOC[pi]*x[pi(0)], pi(0)) >= b[pi(1)], 'cons'  # アイテム i は何れかのパターンで選ぶ．
        prb += Sum(c[p]*x[p])
        prb.solve(silent=silent)
        assert prb.status == NuoptStatus.OPTIMAL
        return x[x[p].val>0].val, prb['cons'].dual

    def solve_dual_problem_for_selection_of_pattern(self, *, silent=True):  # 双対問題
        """solve_primal_problem_for_selection_of_pattern の双対問題"""
        pi = Element(value=self.pattern)
        i = Element(set=pi(1).set)
        b = Parameter(index=pi(1), value=self.b)
        c = Parameter(index=pi(0), value=self.c)
        ASSOC = Parameter(index=pi, value=self.pattern)
        omega = Variable(index=pi(1), lb=0, name='ω')

        prb = Problem(type=max)
        prb += Sum(ASSOC[pi]*omega[pi(1)], pi(1)) <= c[pi(0)]
        prb += Sum(b[i]*omega[i])
        prb.solve(silent=silent)
        assert prb.status == NuoptStatus.OPTIMAL
        return omega.val

    @abstractmethod
    def solve_pricing_problem(self, dualval, *, silent=True):
        """組合せを 1 つ生成する"""
        pass  # return y.val, prb.objective.val

    def update_pattern(self, yval):
        """update self.pattern, self.c"""
        self.pattern.update({(self.patternnum,i): v for (i,), v in yval.items()})

    def solve(self, *, eps=0, maxiter=100):
        # 初期パターンを生成する．
        self.create_init_pattern()

        objval = eps + 1e-6  # eps より大きい適当な値
        while objval > eps and self.iternum < maxiter:
            self.__iternum += 1

            # 限定主問題 `P|_k` を求解して，実行可能解 `\vec{x}_k^{\star}` を構成する．
            # 限定双対問題 `D|_k` を求解して，最適解 `\vec{\omega}_k^*` を求める．
            _, dualval = self.solve_primal_problem_for_selection_of_pattern(type=float)  # 主問題の双対変数を使う．
            # dualval = self.solve_dual_problem_for_selection_of_pattern()  # 双対問題を使う．

            # 価格付け問題を求解して，パターン番号集合を更新する．
            yval, objval = self.solve_pricing_problem(dualval)
            # print('update:', dict(yval.items()))
            self.update_pattern(yval)
            print(f'{self.iternum}: {objval}', flush=True)

        # 実行可能解を取得する．
        yval, _ = self.solve_primal_problem_for_selection_of_pattern()
        return yval

上記のコードの中で，初期集合の生成部分 create_init_pattern と価格付け問題の求解部分 solve_pricing_problem は，個別の問題に応じて設計するサブクラスに具体的な実装を記述することになる．ここでは初期集合でパターンの検証用に全列挙を設定することもできるよう，幾つかの配慮がなされている．

    @abstractmethod
    def create_init_pattern(self, *, pattern, b=1, c=1):
        """set self.pattern, self.b, self.c"""
        self.pattern = pattern
        self.b = b
        self.c = c
        self.__iternum = 0
        self.__pn0 = len(pattern)

    @abstractmethod
    def solve_pricing_problem(self, dualval, *, silent=True):
        """組合せを 1 つ生成する"""
        pass  # return y.val, prb.objective.val

限定主問題または限定双対問題の求解はそれぞれ solve_primal_problem_for_selection_of_pattern と solve_dual_problem_for_selection_of_pattern が担い，強く共通化されている部分である．

    def solve_primal_problem_for_selection_of_pattern(self, *, type=bin, silent=True):  # 主問題
        """生成されたパターンから b を満たす組合せを選択する"""
        pi = Element(value=self.pattern)
        p = Element(set=pi(0).set)
        b = Parameter(index=pi(1), value=self.b)
        c = Parameter(index=pi(0), value=self.c)
        ASSOC = Parameter(index=pi, value=self.pattern)  # パターン p にアイテム i を含むか (0-1)．
        x = Variable(index=p, type=type, lb=0)  # パターン p を選択するか (type=bin)，幾つ選択するか (type=int)．

        prb = Problem()
        prb += Sum(ASSOC[pi]*x[pi(0)], pi(0)) >= b[pi(1)], 'cons'  # アイテム i は何れかのパターンで選ぶ．
        prb += Sum(c[p]*x[p])
        prb.solve(silent=silent)
        assert prb.status == NuoptStatus.OPTIMAL
        return x[x[p].val>0].val, prb['cons'].dual

    def solve_dual_problem_for_selection_of_pattern(self, *, silent=True):  # 双対問題
        """solve_primal_problem_for_selection_of_pattern の双対問題"""
        pi = Element(value=self.pattern)
        i = Element(set=pi(1).set)
        b = Parameter(index=pi(1), value=self.b)
        c = Parameter(index=pi(0), value=self.c)
        ASSOC = Parameter(index=pi, value=self.pattern)
        omega = Variable(index=pi(1), lb=0, name='ω')

        prb = Problem(type=max)
        prb += Sum(ASSOC[pi]*omega[pi(1)], pi(1)) <= c[pi(0)]
        prb += Sum(b[i]*omega[i])
        prb.solve(silent=silent)
        assert prb.status == NuoptStatus.OPTIMAL
        return omega.val

PySIMPLE では主問題を求解後に，対応する双対問題についても参照できるため，これらどちらかを求解するだけでよい．

列生成法で行う逐次的求解およびパターン番号集合の更新処理は solve が担う．

PySIMPLE Sample Code

    def solve(self, *, eps=0, maxiter=100):
        # 初期パターンを生成する．
        self.create_init_pattern()

        objval = eps + 1e-6  # eps より大きい適当な値
        while objval > eps and self.iternum < maxiter:
            self.__iternum += 1

            # 限定主問題 `P|_k` を求解して，実行可能解 `\vec{x}_k^{\star}` を構成する．
            # 限定双対問題 `D|_k` を求解して，最適解 `\vec{\omega}_k^*` を求める．
            _, dualval = self.solve_primal_problem_for_selection_of_pattern(type=float)  # 主問題の双対変数を使う．
            # dualval = self.solve_dual_problem_for_selection_of_pattern()  # 双対問題を使う．

            # 価格付け問題を求解して，パターン番号集合を更新する．
            yval, objval = self.solve_pricing_problem(dualval)
            # print('update:', dict(yval.items()))
            self.update_pattern(yval)
            print(f'{self.iternum}: {objval}', flush=True)

        # 実行可能解を取得する．
        yval, _ = self.solve_primal_problem_for_selection_of_pattern()
        return yval

最適性条件の判定は while 文の終了条件として織り込まれている．

2.4.3. サブクラスの実装例 #

ビンパッキング問題に特化したサブクラスの実装例 (bin_packing.py) が次である．

PySIMPLE Sample Code (bin_packing.py)

# -*- coding: utf-8 -*-
from itertools import combinations
from random import sample, randint, seed; seed(0)

from pysimple import *
from column_generator import ColumnGenerator

class BinPacking(ColumnGenerator):
    def __init__(self, BVOL, IVOL):
        self.BVOL = BVOL  # ビンの容量
        self.IVOL = IVOL  # アイテムの容量

    def create_init_pattern(self):
        pattern = dict.fromkeys(enumerate(self.IVOL, 1), 1)
        super().create_init_pattern(pattern=pattern)

    def solve_pricing_problem(self, dualval, *, silent=True):
        """ビンに収まるアイテムの組合せを 1 つ生成する"""
        i = Element(value=self.IVOL.keys())  # アイテム
        IVOL = Parameter(index=i, value=self.IVOL)  # アイテム i の容量
        DUALVAL = Parameter(index=i, value=dualval, name='ω*')
        y = BinaryVariable(index=i)  # アイテム i を使用するかどうか

        prb = Problem(type=max)
        prb += Sum(IVOL[i]*y[i]) <= self.BVOL, 'アイテムをビンに収める'
        prb += Sum(DUALVAL[i]*y[i])
        prb.solve(silent=silent)
        return y[y[i].val>0].val, prb.objective.val

    def visualize(self, zval):
        for (p,), val in zval.items():
            lentimes = [(len_, self.pattern[p,t]) for t, len_ in self.IVOL.items() if (p, t) in self.pattern]
            pstr = '+'.join(f'{l}{(f"*{times}","")[times==1]}' for l, times in lentimes)
            line = '+'+''.join(('-'*(l-1)+'+')*times for l, times in lentimes)
            print(f'pattern{p:2}({pstr:^20}){val:2}枚: {line}')

def create_random_data(*, BVOL, N):  # BVOL: ビンの容量, N: アイテム種類数
    """テスト用のランダムデータを生成する"""
    # アイテム集合
    I = [f'item{i}' for i in range(1, N+1)]
    # 容量重複がない IVOL[i] を dict 形式でランダム作成
    IVOL = dict(zip(I, sample(range(1, BVOL//2+1), N)))
    return BVOL, IVOL

def solve_at_once(BVOL, IVOL, *, silent=True):
    """一度に解く straight forward な定式化"""
    i = Element(value=IVOL.keys())  # アイテム
    b = Element(value=range(len(i.set)))  # ビン (数は上限のアイテム数とする)
    IVOL = Parameter(index=i, value=IVOL)  # アイテム i の容量
    BVOL = Parameter(value=BVOL) # ビンの容量
    z = BinaryVariable(index=(i,b))  # アイテム i をビン b に入れるかどうか
    use = BinaryVariable(index=b)  # ビン b を使用するかどうか

    prb = Problem()
    prb += Sum(z[i,b], b) == 1 , 'いずれかのビンに入れる'
    prb += Sum(IVOL[i]*z[i,b], i) <= BVOL*use[b],  'アイテムをビンに収める'
    prb += Sum(use[b])  # ビンの総数

    prb.solve(silent=silent)
    assert prb.status == NuoptStatus.OPTIMAL
    print(use[use[b].val>0].val)  # 使ったビン
    print(z[z[i,b].val>0].val)  # 選ばれたアイテムとビンの組合せ

def solve_with_enumerate_all_patterns(BVOL, IVOL):
    """全列挙で解く定式化"""
    def enumerate_all_patterns(BVOL, IVOL):
        """容量 BVOL 以下になる IVOL の組合せを全列挙する"""
        v2i = {v: i for i, v in IVOL.items()}  # IVOL -> item
        patterns, p = [], 1
        for i in range(len(v2i)):
            for comb in combinations(v2i, i+1):
                if sum(comb) <= BVOL:
                    patterns.extend((p, v2i[v]) for v in comb)
                    p += 1
        return dict.fromkeys(patterns, 1)

    bp = BinPacking(BVOL, IVOL)
    pattern = enumerate_all_patterns(BVOL, IVOL)
    ColumnGenerator.create_init_pattern(bp, pattern=pattern)
    zval, _ = bp.solve_primal_problem_for_selection_of_pattern()
    bp.visualize(zval)

def solve_with_column_generation_method(BVOL, IVOL):
    """列生成法で解く定式化"""
    bp = BinPacking(BVOL, IVOL)
    zval = bp.solve(eps=1)
    bp.visualize(zval)


if __name__ == '__main__':
    BVOL, IVOL = create_random_data(BVOL=70, N=20)

    # 1. 素直な定式化で解く場合
    # solve_at_once(BVOL, IVOL)

    # 2. 全列挙で解く場合 (BVOL=70, N=20 で 6700 パターン)
    # solve_with_enumerate_all_patterns(BVOL, IVOL)

    # 3. 列生成で解く場合
    solve_with_column_generation_method(BVOL, IVOL)

設計概要で述べた「初期集合の生成」と「価格付け問題の求解」の実装は次のとおり．

    def create_init_pattern(self):
        pattern = dict.fromkeys(enumerate(self.IVOL, 1), 1)
        super().create_init_pattern(pattern=pattern)

    def solve_pricing_problem(self, dualval, *, silent=True):
        """ビンに収まるアイテムの組合せを 1 つ生成する"""
        i = Element(value=self.IVOL.keys())  # アイテム
        IVOL = Parameter(index=i, value=self.IVOL)  # アイテム i の容量
        DUALVAL = Parameter(index=i, value=dualval, name='ω*')
        y = BinaryVariable(index=i)  # アイテム i を使用するかどうか

        prb = Problem(type=max)
        prb += Sum(IVOL[i]*y[i]) <= self.BVOL, 'アイテムをビンに収める'
        prb += Sum(DUALVAL[i]*y[i])
        prb.solve(silent=silent)
        return y[y[i].val>0].val, prb.objective.val

bin_packing.py では検証のために次の三つの関数が実装されている．

solve_at_once (straight forward な定式化による求解)
solve_with_enumerate_all_patterns (最初に全列挙して求解)
solve_with_column_generation_method (列生成法での求解)

def solve_at_once(BVOL, IVOL, *, silent=True):
    """一度に解く straight forward な定式化"""
    i = Element(value=IVOL.keys())  # アイテム
    b = Element(value=range(len(i.set)))  # ビン (数は上限のアイテム数とする)
    IVOL = Parameter(index=i, value=IVOL)  # アイテム i の容量
    BVOL = Parameter(value=BVOL) # ビンの容量
    z = BinaryVariable(index=(i,b))  # アイテム i をビン b に入れるかどうか
    use = BinaryVariable(index=b)  # ビン b を使用するかどうか

    prb = Problem()
    prb += Sum(z[i,b], b) == 1 , 'いずれかのビンに入れる'
    prb += Sum(IVOL[i]*z[i,b], i) <= BVOL*use[b],  'アイテムをビンに収める'
    prb += Sum(use[b])  # ビンの総数

    prb.solve(silent=silent)
    assert prb.status == NuoptStatus.OPTIMAL
    print(use[use[b].val>0].val)  # 使ったビン
    print(z[z[i,b].val>0].val)  # 選ばれたアイテムとビンの組合せ

def solve_with_enumerate_all_patterns(BVOL, IVOL):
    """全列挙で解く定式化"""
    def enumerate_all_patterns(BVOL, IVOL):
        """容量 BVOL 以下になる IVOL の組合せを全列挙する"""
        v2i = {v: i for i, v in IVOL.items()}  # IVOL -> item
        patterns, p = [], 1
        for i in range(len(v2i)):
            for comb in combinations(v2i, i+1):
                if sum(comb) <= BVOL:
                    patterns.extend((p, v2i[v]) for v in comb)
                    p += 1
        return dict.fromkeys(patterns, 1)

    bp = BinPacking(BVOL, IVOL)
    pattern = enumerate_all_patterns(BVOL, IVOL)
    ColumnGenerator.create_init_pattern(bp, pattern=pattern)
    zval, _ = bp.solve_primal_problem_for_selection_of_pattern()
    bp.visualize(zval)

def solve_with_column_generation_method(BVOL, IVOL):
    """列生成法で解く定式化"""
    bp = BinPacking(BVOL, IVOL)
    zval = bp.solve(eps=1)
    bp.visualize(zval)

2.4.4. 求解 #

BVOL=200, N=100 とアイテム数が \(100\) でそれぞれの求解方法を試すと， straight forward な定式化では数分必要となるが，列生成法での求解は数秒で終えることができた．全列挙は言わずもがなである．

技法集

列生成法という技法

用語集

リソース

PySIMPLE 列生成法: https://www.msi.co.jp/solution/nuopt/docs/pysimple/sample/column_generator.html

参考文献

[1]

田中俊二. 大規模組合せ最適化問題に対する数理アプローチの基礎. 計測と制御, 56(12):967–972, 2017. URL: https://doi.org/10.11499/sicejl.56.967.

[2]

宮本裕一郎. はじめての列生成法. オペレーションズ・リサーチ = Communications of the Operations Research Society of Japan : 経営の科学, 57(4):198–204, 2012. URL: https://orsj.org/wp-content/corsj/or57-4/or57_4_198.pdf.

[3]

久保幹夫, 田村明久, and 松井知己. 応用数理計画ハンドブック. 朝倉書店, 2002.

引用書式

BibTeX

ナップサック問題と列生成法

目次