線形計画問題を単体法（シンプレックス法）によって解くコードを Pythonで実装してみた

現在執筆中..

線形計画問題
単体法
単体法の実装
おまけ
- 線形計画問題の変形

線形計画問題

今回解く線形計画問題は不等式標準形とする．なお，任意の線形計画問題は不等式標準形に変換可能なので，このコードを適用させるには不等式標準形に変換（→おまけ）してから適用させる必要がある．ただし，任意の線形計画問題は不等式標準形に変形できるという事実があるので結構使えると思う．
不等式標準形の線形計画問題
$\displaystyle{\text 条件} A\bf{x}\ge\bf{b}, \bf{x}\ge \bf{0}{\text のもとで}\bf{c}^T\bf{x}{\text を最小にする}\bf{x}{\text を求める.}$
また，このような線形計画問題を以下のように表すことにする．
最小化: $\bf{c}^T\bf{x}$
条件: $A\bf{x}\ge\bf{b}, \bf{x}\ge\bf{0}$

具体的な例としては下の3次元空間での線形計画問題を取り扱う．
例1
最小化: $-x_1+2x_2-x_3$
条件: $-2x_1+2x_2\ge-6\\ -x_1-x_2+2x_3\ge-3\\ -x_1-x_2-x_3\ge-3\\x_1, x_2, x_3\ge0$
この場合
$\displaystyle{A=\left(\begin{array}{ccc} -2&2&0 \\-1&-1&2 \\-1&-1&-1\end{array}\right)}$ , $\displaystyle{\bf{b}=\left(\begin{array}{c} -6 \\-3 \\-3\end{array}\right)}$ , $\displaystyle{\bf{c}=\left(\begin{array}{c} -1 \\2 \\-1\end{array}\right)}$ , $\displaystyle{\bf{x}=\left(\begin{array}{c} x_1 \\x_2 \\x_3\end{array}\right)}$
とすれば，上の
最小化: $\bf{c}^T\bf{x}$
条件: $A\bf{x}\ge\bf{b}, \bf{x}\ge\bf{0}$
という形で表す事ができるとわかる．とにかく条件を満たす $x_1, x_2, x_3$ で目標関数を最小にするものを求めるということを行う．そのような線形計画問題のための解法の1つが単体法（シンプレックス法）である．というわけで最初に単体法の解説を簡単に書いておいて，その後pythonのコードに焼き直す．

単体法

単体法の

単体法の実装

とにかく上のアルゴリズムを Pythonでコーディングすると↓のようになる．

from copy import deepcopy
import numpy as np
import traceback


class Simplex():
    def __init__(self, A, b, c):
        self.A = A
        self.b = b
        self.c = c
        self.m, self.n = A.shape  # m: 方程式の数, n: 変数の次元

        # 辞書の定義
        self.dic = np.zeros((self.m + 1, self.n + 1), dtype='float64')
        for i in range(self.m):
            for j in range(self.n):
                self.dic[i + 1, j + 1] += A[i, j]
        for j in range(self.n):
            self.dic[0, j + 1] += c[j]
        for i in range(self.m):
            self.dic[i + 1, 0] += -b[i]
        self.x_n = np.zeros(self.n, dtype='float64')  # 非基底変数
        self.x_m = np.zeros(self.m, dtype='float64')  # 基底変数
        self.x_n_index = np.arange(0, self.n, dtype='int64')
        self.x_m_index = np.arange(self.n, self.m + self.n, dtype='int64')

    def run(self):
        # 原点が許容解であるか確認
        for j in range(self.n):
            self.x_n[j] = 0.0
        self.x_m = self.update_basis()  # 非基底変数が原点のときの基底変数を代入
        is_feasible = True
        for x in self.x_m:
            is_feasible = is_feasible and (x >= 0.0)
        if not is_feasible:
            try:  # 原点が非許容である場合エラー出力
                raise Exception('The origin is not feasible.')
            except BaseException:
                print(traceback.format_exc())

        z_sol = self.loop()

        x_sol = [0 for _ in range(self.n)]
        x_concat = np.concatenate([self.x_n, self.x_m])
        index_concat = np.concatenate([self.x_n_index, self.x_m_index])
        for j, x in zip(index_concat, x_concat):
            if j < self.n:
                x_sol[j] = x

        return z_sol, np.array(x_sol)

    def loop(self):
        while True:
            self.show_dic()
            # 非基底変数の選択
            j_pivot = -1
            for j, a in enumerate(self.dic[0, 1:]):
                if a < 0:
                    j_pivot = j + 1
                    break
            if j_pivot == -1:  # 負である係数が存在しない場合終了
                return self.get_cost()

            # 基底変数の選択
            base_candidate = {}
            for i in range(1, self.m + 1):
                if self.dic[i, j_pivot] < 0:
                    base_candidate[i] = self.dic[i, 0] \
                        / abs(self.dic[i, j_pivot])
            if len(base_candidate) == 0:
                # 基底の候補がない場合はこの線形計画問題は非有界
                raise Exception('This problem is unbounded.')
            d_min = min(base_candidate.values())
            i_min_list = [kv[0] for kv in base_candidate.items()
                          if kv[1] == d_min]
            # 基底の候補から添え字の最小となる基底を選択
            i_pivot = self.m + self.n
            index_pivot = self.m + self.n
            for i_min in i_min_list:
                if index_pivot >= self.x_m_index[i_min - 1]:
                    i_pivot = i_min
                    index_pivot = self.x_m_index[i_min - 1]

            # 選択された非基底変数を許容性を保存しつつ最大限まで増加させる
            self.x_n[j_pivot - 1] = d_min
            self.x_m[i_pivot - 1] = 0.0

            # ピボット演算でdic, x_m, x_n, x_m_index, x_n_indexを更新
            self.pivot(i_pivot, j_pivot)

    def pivot(self, i_pivot, j_pivot):
        self.x_m[i_pivot - 1], self.x_n[j_pivot - 1] = \
            self.x_n[j_pivot - 1], self.x_m[i_pivot - 1]
        self.x_m_index[i_pivot - 1], self.x_n_index[j_pivot - 1] = \
            self.x_n_index[j_pivot - 1], self.x_m_index[i_pivot - 1]

        dic_prev = deepcopy(self.dic)
        for i in range(self.m + 1):
            if i == i_pivot:
                c = -1.0 / dic_prev[i_pivot, j_pivot]
                self.dic[i, :] = dic_prev[i, :] * c
                self.dic[i_pivot, j_pivot] = -c
            else:
                c = dic_prev[i, j_pivot] / dic_prev[i_pivot, j_pivot]
                self.dic[i, :] = dic_prev[i, :] - dic_prev[i_pivot, :] * c
                self.dic[i, j_pivot] = c

    def update_basis(self):
        return np.dot(self.dic[1:, 1:], self.x_n) + self.dic[1:, 0]

    def get_cost(self):
        return np.dot(self.dic[0, 1:], self.x_n) + self.dic[0, 0]

    def show_dic(self):
        print('    |          |', end='')
        for x in self.x_n_index:
            print('%10.2f ' % (x), end='')
        print('')
        print('-' * (4 + 11 + 11 * self.n))

        print('cost|%10.2f|' % (self.dic[0, 0]), end='')
        for d in self.dic[0, 1:]:
            print('%10.2f ' % (d), end='')
        print('')

        for i in range(1, self.m + 1):
            print('%4d|%10.2f|' % (self.x_m_index[i - 1], self.dic[i, 0]), end='')
            for d in self.dic[i, 1:]:
                print('%10.2f ' % (d), end='')
            print('')
        print('')


# 実際に実行してみる
if __name__ == "__main__":
    A = np.array([[-2, -2, 1], [-2, 0, -4], [4, -3, 1]], dtype='float64')
    b = np.array([4, 4, 1], dtype='float64')
    c = np.array([-2, -1, -1], dtype='float64')
    print(Simplex(A, b, c).run())

まずはSimplexクラスのコンストラクタの説明から．

    def __init__(self, A, b, c):
        self.A = A
        self.b = b
        self.c = c
        self.m, self.n = A.shape  # m: 方程式の数, n: 変数の次元

とりあえず不等式標準形の各定数をメンバ変数として持っておく（結局使うことはなかったが...）．そして，行列Aの形状から方程式の数（不等式の数）： m，変数の次元：n，をメンバ変数に格納する．次に，辞書を定義する．

        # 辞書の定義
        self.dic = np.zeros((self.m + 1, self.n + 1), dtype='float64')
        for i in range(self.m):
            for j in range(self.n):
                self.dic[i + 1, j + 1] += A[i, j]
        for j in range(self.n):
            self.dic[0, j + 1] += c[j]
        for i in range(self.m):
            self.dic[i + 1, 0] += -b[i]
        self.x_n = np.zeros(self.n, dtype='float64')  # 非基底変数
        self.x_m = np.zeros(self.m, dtype='float64')  # 基底変数
        self.x_n_index = np.arange(0, self.n, dtype='int64')
        self.x_m_index = np.arange(self.n, self.m + self.n, dtype='int64')

まず，辞書の形状はm+1×n+1の行列であるので，そのような形状ですべての要素が0である行列self.dicを定義して，それにA, b, cの数値を代入して行っている．また，非基底変数と基底変数を宣言する．単体法はそのアルゴリズムを通して非基底変数と基底変数をスワップすることを何回も行う．したがって，現在の非基底変数または基底変数の何番目の変数が $x_1$ や $x_2$ に対応するのか記憶しておかなければならない．そこで，k番目（0オリジン）の非基底変数に対応する変数が $x_i$ であるとすると，self.x_n_index[k]=iというふうに，非基底変数に対してはself.x_n_index，基底変数に対してはself.x_m_indexに対応する変数の添字を記憶しておく．コンストラクタは以上で終了．
メソッドrun()にて単体法のアルゴリズムを実行する．

    def run(self):
        # 原点が許容解であるか確認
        for j in range(self.n):
            self.x_n[j] = 0.0
        self.x_m = self.update_basis()  # 非基底変数が原点のときの基底変数を代入
        is_feasible = True
        for x in self.x_m:
            is_feasible = is_feasible and (x >= 0.0)
        if not is_feasible:
            try:  # 原点が非許容である場合エラー出力
                raise Exception('The origin is not feasible.')
            except BaseException:
                print(traceback.format_exc())

        z_sol = self.loop()

ここでは最初に原点が許容解であるか確認する．（原点が許容解ではないならば二段階単体法を用いる必要がある．）まずは，全ての非基底変数に0を代入し，それと辞書を用いて基底変数の値を求める．そして，基底変数（今は1度もピボットしていないのでスラック変数）が負であるならば原点は非許容であるので，例外を吐いて終了する．許容解であるならば単体法のループ部分（loop()メソッド，戻り値：目的関数）に入る．なお，このコードでは辞書と非基底変数から基底変数，目的関数を求めるためにそれぞれ，update_basis(), get_cost()というメソッドを定義した．

    def update_basis(self):
        return np.dot(self.dic[1:, 1:], self.x_n) + self.dic[1:, 0]

    def get_cost(self):
        return np.dot(self.dic[0, 1:], self.x_n) + self.dic[0, 0]

さて，単体法の心臓部分のループ．なお，計算の過程がわかるように1ループごとに辞書の内容をメソッドshow_dic()にて出力する．

   def loop(self):
        while True:
            self.show_dic()
            # 非基底変数の選択
            j_pivot = -1
            for j, a in enumerate(self.dic[0, 1:]):
                if a < 0:
                    j_pivot = j + 1
                    break
            if j_pivot == -1:  # 負である係数が存在しない場合終了
                return self.get_cost()

            # 基底変数の選択
            base_candidate = {}
            for i in range(1, self.m + 1):
                if self.dic[i, j_pivot] < 0:
                    base_candidate[i] = self.dic[i, 0] \
                        / abs(self.dic[i, j_pivot])
            if len(base_candidate) == 0:
                # 基底の候補がない場合はこの線形計画問題は非有界
                raise Exception('This problem is unbounded.')
            d_min = min(base_candidate.values())
            i_min_list = [kv[0] for kv in base_candidate.items()
                          if kv[1] == d_min]
            # 基底の候補から添え字の最小となる基底を選択
            i_pivot = self.m + self.n
            index_pivot = self.m + self.n
            for i_min in i_min_list:
                if index_pivot >= self.x_m_index[i_min - 1]:
                    i_pivot = i_min
                    index_pivot = self.x_m_index[i_min - 1]

            # 選択された非基底変数を許容性を保存しつつ最大限まで増加させる
            self.x_n[j_pivot - 1] = d_min
            self.x_m[i_pivot - 1] = 0.0

            # ピボット演算でdic, x_m, x_n, x_m_index, x_n_indexを更新
            self.pivot(i_pivot, j_pivot)

最初に非基底変数の選択．選択できない場合，現在の値が最適解であるのでループ終了．次に基底変数の選択．ここで選択できない場合は線形計画問題が非有界であるため最適解はない．ループ終了．さて，無事に非基底変数・基底変数のペアが選択されたら次は基底変数が0になるまで非基底を増加させる．それ以上増加させたら基底がマイナスになって非許容となるためそこで止める．そしてピボットで諸々の値を更新する．
ピボットの部分はこんな感じ．

    def pivot(self, i_pivot, j_pivot):
        self.x_m[i_pivot - 1], self.x_n[j_pivot - 1] = \
            self.x_n[j_pivot - 1], self.x_m[i_pivot - 1]
        self.x_m_index[i_pivot - 1], self.x_n_index[j_pivot - 1] = \
            self.x_n_index[j_pivot - 1], self.x_m_index[i_pivot - 1]

        dic_prev = deepcopy(self.dic)
        for i in range(self.m + 1):
            if i == i_pivot:
                c = -1.0 / dic_prev[i_pivot, j_pivot]
                self.dic[i, :] = dic_prev[i, :] * c
                self.dic[i_pivot, j_pivot] = -c
            else:
                c = dic_prev[i, j_pivot] / dic_prev[i_pivot, j_pivot]
                self.dic[i, :] = dic_prev[i, :] - dic_prev[i_pivot, :] * c
                self.dic[i, j_pivot] = c

さて，ループが終わったらrun()メソッドに戻ろう．

        z_sol = self.loop()

        x_sol = [0 for _ in range(self.n)]
        x_concat = np.concatenate([self.x_n, self.x_m])
        index_concat = np.concatenate([self.x_n_index, self.x_m_index])
        for j, x in zip(index_concat, x_concat):
            if j < self.n:
                x_sol[j] = x

        return z_sol, np.array(x_sol)

ループが無事終了したら変数z_solに最適値が帰ってくるはずだ．そしてz_solを与える解をx_solに代入してそれらを返そう．x_solの求め方は基底・非基底に対応する添字を小さい方からn個取り出して．その値をx_solに格納して終了！（投げやり）
さて，上の例で示した線形計画問題を解いてみる．
例1
最小化: $-x_1+2x_2-x_3$
条件: $-2x_1+2x_2\ge-6\\ -x_1-x_2+2x_3\ge-3\\ -x_1-x_2-x_3\ge-3\\x_1, x_2, x_3\ge0$

if __name__ == "__main__":
    A = np.array([[-2, 2, 0], [-1, -1, 2], [-1, -1, -1]], dtype='float64')
    b = np.array([-6, -3, -3], dtype='float64')
    c = np.array([-1, 2, -1], dtype='float64')
    print(Simplex(A, b, c).run())

すると標準出力は

    |          |      0.00       1.00       2.00 
------------------------------------------------
cost|      0.00|     -1.00       2.00      -1.00 
   3|      6.00|     -2.00       2.00       0.00 
   4|      3.00|     -1.00      -1.00       2.00 
   5|      3.00|     -1.00      -1.00      -1.00 

    |          |      3.00       1.00       2.00 
------------------------------------------------
cost|     -3.00|      0.50       1.00      -1.00 
   0|      3.00|     -0.50       1.00       0.00 
   4|      0.00|      0.50      -2.00       2.00 
   5|      0.00|      0.50      -2.00      -1.00 

    |          |      3.00       1.00       5.00 
------------------------------------------------
cost|     -3.00|      0.00       3.00       1.00 
   0|      3.00|     -0.50       1.00      -0.00 
   4|      0.00|      1.50      -6.00      -2.00 
   2|      0.00|      0.50      -2.00      -1.00 

(-3.0, array([3., 0., 0.]))

こんな感じ．辞書の更新がされて最後に目的関数-3.0, 最適解 $(x_1, x_2, x_3)=(3, 0, 0)$ ということでいい感じに計算できてましたということでおしまい．

おまけ

線形計画問題の変形

等式標準形から，不等式標準形への変形を行う．
不等式標準形
最小化: $\bf{c}^T\bf{x}$
条件: $A\bf{x}\ge\bf{b}, \bf{x}\ge\bf{0}$
等式標準形
最小化: $\bf{\tilde{c}}^T\bf{\tilde{x}}$
条件: $A\bf{\tilde{x}}=\bf{\tilde{b}}, \bf{\tilde{x}}\ge\bf{0}$
等式標準系から不等式標準系への変形には次の式変形を用いる．