はじめに

数理最適化には興味があってこのブログでもいくつか記事を書いてきたが、つまみ食いばかりして理解があいまいな部分もあるので、もう少しちゃんと基本から学ぶため以下の書籍を読み進めている。

作者:梅谷俊治
発売日: 2020/10/26
メディア: 単行本

ということで今回は、この書籍の2章に書かれている線形計画問題を解くアルゴリズム『2段階単体法』のまとめとPythonによる実装を記す。

線形計画問題

解きたい問題

目的関数が線形関数で、制約条件が全て線形の等式もしくは不等式である最適化問題を線形計画問題という。線形計画問題は全て以下の標準形として定式化できる。

$\displaystyle \begin{aligned} \max_{{\bf x}} &\ \ {\bf c}^T {\bf x} \\ {\rm s.t.} &\ \ {\bf A}{\bf x} \leq {\bf b} \\ &\ \ {\bf x} \geq {\bf 0} \end{aligned} \tag{1}$

最適化問題の変数が $n$ 個、非負制約を除く制約条件が $m$ 個であれば、 ${\bf x} \in \mathbb{R}^n, \ {\bf A} \in \mathbb{R}^{m \times n}, \ {\bf b} \in \mathbb{R}^m$ である。ここで、スラック変数 ${\bf s} \in \mathbb{R}^m$ を導入し不等式制約を等式制約に変換する。

$\displaystyle \begin{aligned} \max_{{\bf x}, \ {\bf s}} &\ \ {\bf c}^T {\bf x} \\ {\rm s.t.} &\ \ {\bf A}{\bf x} + {\bf s} = {\bf b} \\ &\ \ {\bf x} \geq {\bf 0} \\ &\ \ {\bf s} \geq {\bf 0} \end{aligned} \tag{2}$

この問題（2）の解は元の問題と同じ解が得られる。

最適解

線形計画問題（1）の解の存在性は、 ${\bf c}, \ {\bf A}, \ {\bf b}$ の値に応じて以下の3パターンに分かれる。

制約を満たす解が存在しない（実行不能）
目的関数を無限に大きくできる（非有界）
最適解がある

そして最適解があるパターンの場合、その最適解は不等式制約で形成される $n$ 次元凸多面体上の頂点上にある。従って最適化する上では凸多面体上の頂点のみを調べればよい。これを効率よく探索するのが単体法である。

また、凸多面体上の頂点ということは式（1）において、非負制約を含む $n + m$ 個の制約条件のうち少なくとも $n$ 個は等式を満たすことになる。さらに、これを式（2）で考えると等式を満たすということは対応する変数もしくはスラック変数が $0$ になることを意味する*1。つまり、最適解 ${\bf x}^{*}, \ {\bf s}^{*}$ の合わせて $n + m$ 個の要素のうち少なくとも $n$ 個は $0$ になる。

アルゴリズム

単体法の考え方

細かい解説は上記の書籍に譲るとして、ここでは大まかな考え方を述べる。まずは ${\bf x} = {\bf 0}$ が制約条件を満たすとして話を進める。この前提は ${\bf b}$ に負値が含まれていると成り立たないが、その場合の対処法は後述する。とりあえず ${\bf b}$ に負値は含まれないとする。

単体法では最適解候補である凸多面体上の1つの頂点から始め（上記の前提が満たされているなら、原点 ${\bf x} = {\bf 0}$ も頂点の1つなので原点を初期値としてよい）、目的関数が改善する隣接する頂点に移動するという操作を改善できなくなるまで繰り返す。先述の通り、頂点ということは $n$ 個の制約条件が等式を満たすことになるが、隣接する頂点に移動するということは、この等式を満たす制約を1つ入れ替える操作に対応する。

この操作を数式を使って記述する。まず式（2）を少し記号などを書き換える。

$\displaystyle \begin{aligned} \max_{{\bf x}_B, \ {\bf x}_N} &\ \ {\bf c}_N^T {\bf x}_N + {\bf c}_B^T {\bf x}_B \\ {\rm s.t.} &\ \ {\bf N}{\bf x}_N + {\bf B}{\bf x}_B = {\bf b} \\ &\ \ {\bf x}_B \geq {\bf 0} \\ &\ \ {\bf x}_N \geq {\bf 0} \end{aligned} \tag{3}$

${\bf x}_N = {\bf x}, \ {\bf x}_B = {\bf s}, \ {\bf c}_N = {\bf c}, \ {\bf c}_B = {\bf 0}, \ {\bf N} = {\bf A}, \ {\bf B} = {\bf I}$ とすれば式（2）になる。 ${\bf x}_N$ は等式を満たすとした制約に対応する変数を表す。先述の通り等式を満たす制約に対応する変数は $0$ になるので、 ${\bf x}_N$ は ${\bf 0}$ に固定される。この ${\bf x}_N$ に割り当てれた（つまり $0$ に固定する）変数を非基底変数、 ${\bf x}_B$ に割り当てられた変数を基底変数という。また、 ${\bf N}$ を非基底行列、 ${\bf B}$ を基底行列と呼ぶ。単体法は ${\bf x}_N$ を ${\bf 0}$ にするのが本当に最適なのかを調べ、最適でない（まだ大きくできる）なら非基底変数と基底変数を1つ入れ替える。それに伴って対応する ${\bf N}$ と ${\bf B}$ の列を交換し、対応する ${\bf c}_N$ と ${\bf c}_B$ の要素を交換する。これが凸多面体上の頂点を移動する操作に当たる。

例えば、初期値が原点から始まるとすれば、最初は非基底変数には元の最適化問題に含まれる変数 ${\bf x}$ が割り当てられており、基底変数にはスラック変数 ${\bf s}$ が割り当てられている。そして、この状態が最適でなければ基底変数、非基底変数の中から変数を1つずつ適切に選び交換する。また、それらに対応する基底行列、非基底行列の列および係数の要素を交換する。この時の操作のイメージはこんな感じ。

f:id:YamagenSakam:20201229094122p:plain

では、どうやって現在の状態が最適か否かを判断するか？入れ替える基底変数、非基底変数はどう選ぶか？というのが次のお話。

交換する変数の選択と停止、非有界の判定

まず ${\bf B}$ が正則だとすると式（3）における等式制約より、

$\displaystyle {\bf x}_B = {\bf B}^{-1} {\bf b} - {\bf B}^{-1} {\bf N} {\bf x}_N\tag{4}$

が成り立つ。これを式（3）の目的関数に代入すると、

$\displaystyle {\bf c}_B^T {\bf B}^{-1} {\bf b} + {\bar{\bf c}_N} ^T {\bf x}_N \tag{5}$

が得られる。ここで、

$\displaystyle {\bar{\bf c}_N} = {\bf c}_N - {\bf N}^T {\bf B}^{-T} {\bf c}_B \tag{6}$

である。式（5）, （6）より次のことがわかる。

${\bar{\bf c}_N} \leq {\bf 0}$ なら ${\bf x}_N = {\bf 0}$ が最適でこれ以上大きくできない。つまり、この時にアルゴリズムは終了。
${\bar{\bf c}_N}$ に正の要素が含まれるなら、その要素に対応する ${\bf x}_N$ の要素を $0$ より大きな値にすると改善できる。従ってその変数は基底変数にする方がよい。

これで、アルゴリズムの終了条件と交換するべき非基底変数がわかった。では次に交換の対象となる基底変数はどう選ぶべきか？今、選んだ非基底変数を $x_k$ とすると、当然 $x_k$ は大きくした方がよい。一方で、式（4）より、 $x_k$ を大きくしすぎて ${\bf x}_B$ の要素が負になってはいけない（非負制約があるため）。 ${\bf B}^{-1} {\bf N}$ の $x_k$ に対応する列を ${\bar{\bf a}}_k$ とすると、

$\displaystyle {\bf B}^{-1} {\bf b} - x_k {\bar{\bf a}}_k \geq {\bf 0} \tag{7}$

を満たすように $x_k$ を定める必要がある。ここで、 ${\bar{\bf a}}_k \leq {\bf 0}$ ならば $x_k$ をいくらでも大きくできる。つまり、非有界ということになる。

一方、 ${\bar{\bf a}}_k$ の中に正の要素が $m_{+}$ 個含まれる場合、その正の要素の1つを ${\bar a}_{ik}^+$ 、それに対応する ${\bar{\bf b}} = {\bf B}^{-1} {\bf b}$ の要素を ${\bar b}_{i}^+$ とすると、

$\displaystyle 0 \leq x_k \leq \frac{\bar{b}_{i}^+}{\bar{a}_{ik}^+} \tag{8}$

が全ての $i = 1 \cdots m_{+}$ について成り立っている必要がある。実行可能解を移動してる限り式（4）と ${\bf x}_N = {\bf 0}$ より ${\bar{\bf b}} \geq {\bf 0}$ なので式（8）を満たす $x_k$ は存在し、 $x_k$ は以下の値にまで大きくできる。

$\displaystyle x_k = \min_i{\frac{\bar{b}_i^+}{\bar{a}_{ik}^+}} \tag{9}$

そして、この式（9）を最小にする $i$ に対応する基底変数は $0$ になる。つまり非基底変数になる。よって、交換対象の基底変数はこの $i$ に対応する変数を選ぶことになる。

交換する非基底変数、基底変数の選び方をまとめると、

非基底変数からは式（6）の正の要素に対応する変数を選ぶ。なければその時が最適であるためアルゴリズム停止。
基底変数からは式（9）を最小にする $i$ に対応する変数を選ぶ。 ${\bar{\bf a}}_k \leq {\bf 0}$ ならば非有界。

初期化（2段階単体法）

ここまでは、 ${\bf x} = {\bf 0}$ が制約を満たすとしこれを初期値として話を進めてきたが、 ${\bf b}$ に負値が含まれる場合 ${\bf x} = {\bf 0}$ は制約を満たさない。そもそも制約を満たす解などなく実行不能である可能性もある。そういことを考慮すると、まず実行可能解が存在するかを判断し、存在するなら適切な初期値を見つける必要がある。そのために以下の補助問題を考える。*2

$\displaystyle \begin{aligned} \min_{{\bf x}, \ x_a} &\ \ x_a \\ {\rm s.t.} &\ \ {\bf A}{\bf x} - x_a {\bf I} \leq {\bf b} \\ &\ \ {\bf x} \geq {\bf 0} \\ & \ \ x_a \geq 0 \end{aligned} \tag{10}$

追加した変数 $x_a$ は人工変数と呼ばれ、式（10）はこの $x_a$ を大きくすれば制約を満たせるので実行可能解は存在する。補助問題は最適化問題（1）が解を持つためには、どれくらい制約条件を修正するべきか？を求める問題とみなせる。従って最適解が $x_a^{*}=0$ であれば制約条件を修正する必要がない、つまり制約条件を満たす実行可能な解が存在する。そして、この時の ${\bf x}^{*}$ は制約条件を満たすのでこれを初期値として本来解きたい問題を解くようにすればよい。逆に最適解で $x_a^{*}>0$ であれば実行不能である。

実装としては、補助問題を先述の単体法の手続きで解いて、最適解 $x_a^{*}=0$ となった時の ${\bf B}, \ {\bf N}, \ {\bf c}_B, \ {\bf c}_N$ から $x_a$ に対応する列および要素を抜いたものを初期値としてもう一度単体法を行えばよい。このようにまず補助問題を解いて初期値を求め、その後解きたい問題を解く手続きを2段階単体法という。

ただし、 ${\bf b}$ には負値が含まれる場合、この補助問題についても原点 ${\bf x} = {\bf 0}, x_a = 0$ は制約を満たさない。そこで、 ${\bf b}$ の最小要素に対応するスラック変数を非基底変数、人工変数を基底変数とする入れ替えを行う。こうすれば実行可能になるのでこれを補助問題の初期値とすればよい。

処理の流れまとめ

ということで、2段階単体法の処理をまとめると以下のようになる。

2段階単体法の流れ

式（10）の補助問題を単体法で解く。
最適値が0でなければ実行不能として終了。
最適値が0になった場合、その時の非基底変数、基底変数の割り当てに基づき、 ${\bf B}, \ {\bf N}, \ {\bf c}_B, \ {\bf c}_N$ を並び替え。
これを初期値として単体法で式（1）の問題を解く。

そして単体法本体部分はこんな感じ。

単体法の流れ

${\bar{\bf b}} = {\bf B}^{-1} {\bf b}$ と式（6）で ${\bar {\bf c}}_N$ を求める。
${\bar {\bf c}}_N \leq {\bf 0}$ であれば最適解 ${\bf x}_B = {\bar{\bf b}}, {\bf x}_N = {\bf 0}$ 、最適値 ${\bf c}_B^{T} {\bf x}_B$ として終了。
そうでなければ、 ${\bar {\bf c}}_N$ の正の要素 $c_k$ に対応する非基底変数 $x_k$ を選択。
${\bf B}^{-1} {\bf N}$ の $x_k$ に対応する列 ${\bar{\bf a}}_k$ が ${\bar{\bf a}}_k \leq {\bf 0}$ であれば非有界として終了。
そうでなければ、式（9）の $i$ に対応する基底変数を選択。
選択した非基底、基底変数およびそれらに対応する ${\bf N}, {\bf B}$ の列、 ${\bf c}_N, {\bf c}_B$ の要素を交換。
1.に戻る。

実装

以下が実装した2段階単体法のpythonコード。

import numpy as np
import copy

def swapNB(cB, cN, B, N, Bidx, Nidx, swap_pos_B, swap_pos_N):
    # idx swap
    tmp = Nidx[swap_pos_N]
    Nidx[swap_pos_N] = Bidx[swap_pos_B]
    Bidx[swap_pos_B] = tmp
    
    # N, B swap
    tmp = np.copy(N[:, swap_pos_N])
    N[:, swap_pos_N] = B[:, swap_pos_B]
    B[:, swap_pos_B] = tmp

    #cN, cB, swap
    tmp = cN[swap_pos_N]
    cN[swap_pos_N] = cB[swap_pos_B]
    cB[swap_pos_B] = tmp   

# 単体法本体
def Simplex(cB_init, cN_init, B_init, N_init, Bidx_init, Nidx_init, b):
    Bidx = copy.copy(Bidx_init)
    Nidx = copy.copy(Nidx_init)
    cB = copy.copy(cB_init)
    cN = copy.copy(cN_init)
    B = copy.copy(B_init)
    N = copy.copy(N_init)
    while(1):
        invB = np.linalg.inv(B)
        b_ = invB @ b
        y = invB.T @ cB
        cN_ = cN - N.T @ y
        # 全て0以下なら停止
        if (cN_ <= 0).all() == True:
            xB = b_
            return cB @ xB, xB, cB, cN, B, N, Bidx, Nidx
        # cN_の正の要素の中で一番最初の要素に
        # 対応する変数交換する非基底変数として選択
        k = np.where(cN_ > 0)[0][0]
        
        ak_ = invB @ N[:, k]
        # 全て0以下なら非有界
        if((ak_ <= 0).all() == True):
            raise Exception('unbounded error')
        # ak_の中で、b_i/a_k_iが最小になる要素iに
        # 対応する変数を交換する基底変数として選択
        tmp = np.ones(b_.shape) * np.inf
        tmp[ak_ > 0] = b_[ak_ > 0] / ak_[ak_ > 0]
        l = np.argmin(tmp)
        swapNB(cB, cN, B, N, Bidx, Nidx, l, k)

# 初期値を探す（2段階単体法の1段目）
def SearchFeasibleDicstionary(cB_init, cN_init, B_init, N_init, 
                              Bidx_init, Nidx_init, b):
    # bに負が含まれる場合、実行可能な辞書を探す
    bmin = np.min(b)
    if bmin < 0:
        Bidx_aux = copy.copy(Bidx_init)
        Nidx = copy.copy(Nidx_init)
        cB_aux = copy.copy(cB_init)
        cN = copy.copy(cN_init)
        B_aux = copy.copy(B_init)
        N = copy.copy(N_init)
        
        # 人工変数のindex
        artificial_var_idx = cB_init.shape[0] + cN_init.shape[0]

        # 人工変数を非基底変数にする
        Nidx_aux = copy.copy(Nidx)
        Nidx_aux.append(artificial_var_idx)

        # 人工変数に対応する要素と列を追加
        cN_aux = np.zeros(cN.shape[0] + 1)
        # min x_a -> max -x_a　という最適化問題を解くため係数は-1
        cN_aux[cN.shape[0]] = -1
        N_aux = np.concatenate([N, np.ones((N.shape[0], 1)) * -1], axis = 1)

        # bの最小値に対応する基底変数と人工変数を入れ替え
        k = np.argmin(b)
        swapNB(cB_aux, cN_aux, B_aux, N_aux, 
                Bidx_aux, Nidx_aux, k, N_aux.shape[1] - 1)

        # 補助問題を解いて実行可能解を見つける
        z, _, res_cB, res_cN, res_B, res_N, res_Bidx, res_Nidx = \
         Simplex(cB_aux, cN_aux, B_aux, N_aux, Bidx_aux, Nidx_aux, b)
        
        # 最適解が負なら実行不能
        if z < 0:
            raise Exception('infeasible error')  
        # 得られた辞書のidxから人工変数を削除
        if artificial_var_idx in res_Nidx:
            # 人工変数が非基底変数にあるならそのまま削除
            res_Nidx.remove(artificial_var_idx)
        else:
            #退化して最適解の基底変数に人工変数が含まれる場合は
            # 非基底変数1つと交換する
            r = res_Bidx.index(artificial_var_idx)
            for i in range(len(res_Nidx)):
                swapNB(res_cB, res_cN, res_B, res_N, res_Bidx, res_Nidx, r, i)
                # 基底行列が正則か確かめる
                if(np.linalg.matrix_rank(res_B) == res_B.shape[0]):
                    res_Nidx.remove(artificial_var_idx)
                    break
                # 正則でなければ元に戻して、交換するindexを変えてもう一度
                swapNB(res_cB, res_cN, res_B, res_N, res_Bidx, res_Nidx, i, r)

        # 得られた基底、非基底変数の割り当てに基づいて
        # 基底、非基底行列、係数を再構成
        res_cNB = np.concatenate([cN_init, cB_init])[res_Nidx + res_Bidx]
        res_cN = res_cNB[:cN_init.shape[0]]
        res_cB = res_cNB[cN_init.shape[0]:]
        res_NB = np.concatenate([N_init, B_init],
                                        axis = 1)[:, res_Nidx + res_Bidx]
        res_N = res_NB[:, :cN_init.shape[0]]
        res_B = res_NB[:, cN_init.shape[0]:]
        return res_Bidx, res_Nidx, res_cB, res_cN, res_B, res_N
        
    # 負がなければそのままreturn
    return Bidx_init, Nidx_init, cB_init, cN_init, B_init, N_init

# max c.T @ x
# s.t. A @ x <= b
#      x >= 0
def SolveLP(c, A, b):
    cB = np.zeros(A.shape[0])
    cN = np.copy(c)
    N = np.copy(A)
    B = np.eye(A.shape[0])
    # 実際の変数のインデックス(0~n-1)
    actual_variable_idx = list(range(c.shape[0]))
    # スラック変数のインデックス(n~n+m-1)
    slack_variable_idx = list(range(c.shape[0], c.shape[0] + A.shape[0]))
    Nidx = copy.copy(actual_variable_idx)
    Bidx = copy.copy(slack_variable_idx)
   
    # 2段階単体法
    new_Bidx, new_Nidx, new_cB, new_cN, new_B, new_N = \
        SearchFeasibleDicstionary(cB, cN, B, N, Bidx, Nidx, b)
    z, xB_opt, cB_opt, cN_opt, B_opt, N_opt, Bidx_opt, Nidx_opt = \
        Simplex(new_cB, new_cN, new_B, new_N, new_Bidx, new_Nidx, b)
    
    # 得られた解から実際の変数（スラック変数以外の変数）を並び替えて出力
    x_opt = np.concatenate([xB_opt, np.zeros(cN_opt.shape[0])])
    x_actual_opt = x_opt[np.argsort(Bidx_opt + Nidx_opt)][actual_variable_idx]
    return z, x_actual_opt

実装は上述の処理を割と素直に書いた。現在どの変数が非基底変数でどの変数が基底変数にであるかについては、変数の添え字をNidx, Bidxというリストで持つことで管理している。 numpyは行優先にもかかわらず列を取り出す操作が多々あるので、効率的には行列を転置して行を取り出す実装にした方がよかったかもしれない*3。もう少し効率面で改善できそうだったり数値計算的に問題がありそうな箇所はあるが、まぁとりあえず書籍に載ってる演習問題を入力してみて、正しい答えを返すのでOKとする。

上述のソースコードに加えて検証用コード、検証に使ったテストケースを書いたテキストファイルを以下のgithub リポジトリに上げている。

github.com

1つこの実装方法でハマったのは、1段階目の補助問題を解く際に退化*4して人工変数が基底変数の中に含まれる最適値 $0$ の最適解が得られるケースがあったこと。補助問題の最適値が $0$ なのでもちろん実行可能なのだが、こうなると人工変数に対応する列を削除する際、基底行列から削除してしまい基底行列が正方でなくなってしまう。こうなった場合の対処として、人工変数といずれかの非基底変数を入れ替える処置を取っている。補助問題の定義より最適解が $0$ なら人工変数も $0$ になっているはずなので、同じく $0$ となる非基底変数を入れ替えても問題ない。ただし、入れ替えた結果の基底行列 ${\bf B}$ が正則でなくなってしまうこともある。そのため、1つずつ入れ替えを試して正則になったものを採用するようにした（この問題に対するもう少しスマートな解決方法はないものか）。