floptで単体法(Simplex)実装

本稿では、floptモデリングツールで定義した線形計画問題に対し単体法を実行して解の取得を行います。Wikipediaにある問題を例に実装してみましょう。

主問題の定義

まず、Variable3つを用いて目的関数と、制約を作成します。

import flopt

x = flopt.Variable("x", lowBound=0)
y = flopt.Variable("y", lowBound=0)
z = flopt.Variable("z", lowBound=0)

prob = flopt.Problem(name="Primary problem", sense="Minimize")
prob += -2 * x - 3 * y - 4 * z
prob +=  3 * x + 2 * y +     z <= 10
prob +=  2 * x + 5 * y + 3 * z <= 15

問題の表示

prob.show()
>>> Name: Primary problem
>>>   Type         : Problem
>>>   sense        : Minimize
>>>   objective    : -2*x-(3*y)-(4*z)
>>>   #constraints : 2
>>>   #variables   : 3 (Continuous 3)
>>> 
>>>   C 0, name None, 3*x+2*y+z-10 <= 0
>>>   C 1, name None, 2*x+5*y+3*z-15 <= 0
>>>
>>>   V 0, name y, Continuous 0 <= y <= None
>>>   V 1, name z, Continuous 0 <= z <= None
>>>   V 2, name x, Continuous 0 <= x <= None

単体法(Simplex)

線形計画問題において、その最適解があるとすれば、少なくともそのうちの一つは実行可能領域の凸多面体の頂点に必ず存在するという性質があります。そのため、頂点を一つ一つ見ていけば良いですが、制約や変数が多いと頂点の個数自体かなりの数になります。単体法は頂点をめぐる効率の良い順序を提供し、全ての頂点を調べ上げる必要をなくしてくれるアルゴリズムになっています。

等式制約に変換

まず、実装する単体法のために不等式制約をスラック変数を用いて等式制約に変換します。.toEq()メソッドを用いることで、floptが自動でスラック変数の追加を行ってくれます。

prob = prob.toEq()
prob.show()
>>> Name: Primary problem
>>>   Type         : Problem
>>>   sense        : Minimize
>>>   objective    : -2*x-(3*y)-(4*z)
>>>   #constraints : 2
>>>   #variables   : 5 (Continuous 5)
>>> 
>>>   C 0, name None, 3*x+2*y+z-10+__0_slack == 0
>>>   C 1, name None, 2*x+5*y+3*z-15+__1_slack == 0
>>> 
>>>   V 0, name __1_slack, Continuous 0 <= __1_slack <= None
>>>   V 1, name __0_slack, Continuous 0 <= __0_slack <= None
>>>   V 2, name y, Continuous 0 <= y <= None
>>>   V 3, name z, Continuous 0 <= z <= None
>>>   V 4, name x, Continuous 0 <= x <= None

線形計画の形式に変換

flopt.convertを用いて線形計画の形式データとして、下記の形にベクトルや行列データを取り出します。lp.Aやlp.cでそれぞれのデータにアクセスできます。 前述の操作により制約はすべて等式に変換されたので、A, bのみに制約の情報が格納されています。

lp = flopt.convert.LpStructure.fromFlopt(prob)
lp.show()
>>> LpStructure
>>> obj  c.T.dot(x) + C
>>> s.t. Gx <= h
>>>      Ax == b
>>>      lb <= x <= ub
>>> 
>>> ...

下記では最大化の場合の単体法を実装するので最小化の場合は係数ベクトルcの符号を反転させます。

c, A, b = lp.c, lp.A, lp.b
if prob.sense.lower() == "minimize":
    c = -c

基底変数と非基底変数の設定

基底変数(basis)と非基底変数(nonbasis)はインデックスで管理します。ここでは簡単に変数の前からn-m個を非基底変数、それ以外を基底変数とします、ここでmは制約の個数、nは変数の個数です。非基底変数は値を0に固定する変数ですね、その上で等式制約を満たすような基底変数の値、基底解を初期解として設定します。

# m is the number of equal constraints
# n is the number of variables
m, n = A.shape

# create a solution container
x = np.zeros(n)

# create indexes of basis and nonbasis variables
indexes = list(range(n))
basis = indexes[n - m :]
nonbasis = indexes[: n - m]
x[basis] = np.linalg.inv(A[:, basis]).dot(b)
x

Step1: 被約費用の計算

基底行列(B)と非基底行列(N)を制約行列Aから取り出して、それを元に被約費用(reduced_cost)の計算を行います。reduced_cost[i] > 0であれば非基底変数nonbasis[i]を増加させることで、目的関数を増加させることができるため、基底解は最適解ではありません。逆に全てのreduced_costの要素が0以下であれば非基底変数の変動による目的関数の増加は見込めないので、今の基底解が最適解がであることが分かります。

# set basis and nonbasis matrix
B, N = A[:, basis], A[:, nonbasis]

# calculate reduced_cost
B_inv = np.linalg.inv(B)
b_bar = B_inv.dot(b)
y = B_inv.T.dot(c[basis])
reduced_cost = c[nonbasis] - N.T.dot(y)
print("\treduced_cost", reduced_cost)
if np.all(reduced_cost <= 0):
    print("Optimal found !!")
    break

Step2: 基底変数と非基底変数の入れ替え

reduced_costが正であるような非基底変数の中で最もインデックスが小さいもの(nonbasis[k])を選択し、その非基底変数を増加させます。これを最小添字規則と呼びます。一方で、reduced_costの要素が大きい非基底変数を増加させるのは最大係数規則なのですが、巡回が生じてアルゴリズムが終了しない場合があるため、今回は最小添字規則での実装を行っています。また非基底変数を最大どれほど増加できるかは等式制約から自然と導くことができます。

# select nonbasis pivot k
k = np.where(reduced_cost > 0)[0][0].item()

# calculate update value of nonbasis pivot k
a_bar = B_inv.dot(N[:, k])
if np.all(a_bar <= 0):
    print("Unbounded")
    break
tmp = np.where(a_bar > 0, a_bar, 1.0)
theta_ = np.where(a_bar > 0, b_bar / tmp, np.inf)
theta = np.min(theta_)

# select basis pivot l
l = np.argmin(theta_)
print("\tk", k, ", theta", theta)

# update solution x
x[nonbasis[k]] = theta
x[basis] = b_bar - theta * a_bar

# swap basis[l] and nonbasis[k]
basis.append(nonbasis.pop(k))
nonbasis.append(basis.pop(l))
basis.sort()
nonbasis.sort()

ここで, a_barの要素が全て負であれば、nonbasis[k]を制約を満たし目的関数を増加させながら、無限に増加させることができるため、問題が非有界であることが分かります。

あとは、このStepが終了したら、Step1に戻り、最適解が見つかるか、非有界であることが見つかるまでループを繰り返すだけです。

補助問題

前述した単体法で、最初に非基底変数を変数の前から順に決めましたが、その際に実行可能解になるような基底解が見つからない場合があります。その場合単体法がうまく動作しません。例えば次のような問題だと、上記のように非基底変数を決めるとxとyが非基底変数となり、これらを0に固定してしまうと、実行可能解が見つかりません。

import flopt

x = flopt.Variable("x", lowBound=0)
y = flopt.Variable("y", lowBound=0)
z = flopt.Variable("z", lowBound=0)

prob = flopt.Problem(name="Primary problem", sense="Minimize")
prob += - 2 * x - 3 * y
prob +=   3 * x + 2 * y == 10
prob += - 2 * x - 5 * y  <= -6

そこでまず、単体法を実行する前に実行可能解を見つけるための問題を解きます。これを補助問題と呼び、補助問題とその後実行する単体法を合わせて二段階単体法と呼びます。補助問題ではそれぞれの制約の実行可能度合いを表す補助変数を用いることで、実行可能解が存在するか、存在するならどのような値かを見つけます。

不等式制約に変換

.toIneq()メソッドを用いることで等式制約を二つの不等式制約に分解できます。これを用いて、以降の記述を簡単にするために全ての制約を不等式制約に変換しておきます。

# create auxiliary problem
ineq_prob = prob.toIneq()
ineq_prob.name = "Auxiliary problem"
ineq_prob.sense = "minimize"

ineq_prob.show()
>>> Name: Auxiliary problem
>>>   Type         : Problem
>>>   sense        : minimize
>>>   objective    : -2*x-(3*y)
>>>   #constraints : 3
>>>   #variables   : 2 (Continuous 2)
>>> 
>>>   C 0, name None, 3*x+2*y-10 <= 0
>>>   C 1, name None, -(3*x+2*y-10) <= 0
>>>   C 2, name None, -2*x-(5*y)+6 <= 0

補助変数の追加

そして二つの正の値を取る補助変数(artificial, artificial)を追加します。補助問題の目的関数はartificialとして、全ての不等式制約の左辺からartificialを差し引き、artificial = _artificialという制約を追加して準備終了です。この目的関数が不等式制約の乖離の度合いを表しており、もし目的関数値が0で実行可能解(すなわち、最適解)が見つかれば、主問題の制約を全て満たす解も見つかっていることになります。その主問題の実行可能解をもとに単体法を実行する、という流れです。

# add artificial variables
artificial = flopt.Variable("artificial", lowBound=0)
_artificial = flopt.Variable("_artificial", lowBound=0)

# set objective function ( = _artificial )
ineq_prob += _artificial

# add artificial in inequality constraints 
for const in ineq_prob.constraints:
    const.expression -= artificial

# add a constraint
ineq_prob += artificial == _artificial
ineq_prob.show()
>>> Name: Auxiliary problem
>>>   Type         : Problem
>>>   sense        : minimize
>>>   objective    : _artificial+0
>>>   #constraints : 4
>>>   #variables   : 4 (Continuous 4)
>>> 
>>>   C 0, name None, 3*x+2*y-10-artificial <= 0
>>>   C 1, name None, -(3*x+2*y-10)-artificial <= 0
>>>   C 2, name None, -2*x-(5*y)+6-artificial <= 0
>>>   C 3, name None, artificial-_artificial == 0

また、この問題は全ての制約を等式制約に変換すれば、先ほど実装した単体法を適用して解を求めることができます。（一石二鳥）。

解の保存

単体法による求解後、得られた解を変数に保存することで、その後の単体法へ値を引き継ぎます。

for var, value in zip(lp.x, x):
    var.setValue(value)

実行可能解を元に単体法を実行

これは簡単で、flopt.Value(lp.x)で先ほど保存した変数の値を取り出したのち、それを解x変数に格納し、値が0の変数を非基底変数とするだけです。個数だけ気をつければOK.

x = flopt.Value(lp.x).astype(np.float64)
nonbasis = np.where(x == 0)[0].tolist()
basis = np.where(x != 0)[0].tolist()
while len(nonbasis) > n - m:
    basis.append(nonbasis.pop())
assert len(nonbasis) == n - m

単体法コードのまとめ

これまでのコードのまとめです。純粋な単体法を補助問題のどちらも扱えるコードです。

import flopt
import flopt.convert

import enum
import numpy as np


class SolverStatus(enum.Enum):
    Optimal = enum.auto()
    Unbounded = enum.auto()
    Feasible = enum.auto()
    Infeasible = enum.auto()

def simplex(prob, auxiliary=False, initial_values=False):

    assert all(var.type() == flopt.VarContinuous for var in prob.getVariables()), "all variables must be continuous"
    assert all(var.getLb() == 0 for var in prob.getVariables()), "lower bounds of all variables must be 0"

    if auxiliary:
        assert all(
            var.name not in {"artificial", "_artificial"} for var in prob.getVariables()
        )

        # create auxiliary problem
        ineq_prob = prob.toIneq()
        ineq_prob.name = "Auxiliary problem"
        ineq_prob.sense = "minimize"

        # add artificial variables
        artificial = flopt.Variable("artificial", lowBound=0)
        _artificial = flopt.Variable("_artificial", lowBound=0)
        
        # set objective function ( = _artificial )
        ineq_prob += _artificial

        # add artificial in inequality constraints 
        for const in ineq_prob.constraints:
            const.expression -= artificial
            
        # add a constraint
        ineq_prob += artificial == _artificial

        prob = ineq_prob.toEq()
        variables = [_artificial, artificial] + sorted(
            prob.getVariables() - {artificial, _artificial},
            key=lambda v: ("__" in v.name, v.name),
        )

    else:

        prob = prob.toEq()
        variables = None

    # conver LP structure with all equal constraints
    lp = flopt.convert.LpStructure.fromFlopt(prob, x=variables)
    prob.show()

    # objective c^T x
    # constraints A x = b
    c, A, b = lp.c, lp.A, lp.b
    if prob.sense.lower() == "minimize":
        c = -c

    # m is the number of equal constraints
    # n is the number of variables
    m, n = A.shape

    # x is a solution
    x = np.zeros(n)

    # create indexes of basis and nonbasis variables
    indexes = list(range(n))

    if auxiliary:

        basis = indexes[:1] + indexes[n - m + 1 :]
        nonbasis = indexes[1 : n - m + 1]

        k = np.argmin(lp.b[:-1]) + 1
        basis.append(nonbasis.pop(0))
        nonbasis.append(basis.pop(k))
        basis.sort()
        nonbasis.sort()

        x[basis] = np.linalg.inv(A[:, basis]).dot(b)

    elif initial_values:

        x = flopt.Value(lp.x).astype(np.float64)
        nonbasis = np.where(x == 0)[0].tolist()
        basis = np.where(x != 0)[0].tolist()
        while len(nonbasis) > n - m:
            basis.append(nonbasis.pop())
        assert len(nonbasis) == n - m

    else:

        basis = indexes[n - m :]
        nonbasis = indexes[: n - m]
        x[basis] = np.linalg.inv(A[:, basis]).dot(b)

    status = None
    n_iter = 0
    while True:
        print("iter", n_iter, "x", x, "basis", basis, "nonbasis", nonbasis, "obj", c.dot(x))

        # set basis and nonbasis matrix
        B, N = A[:, basis], A[:, nonbasis]

        # calculate reduced_cost
        B_inv = np.linalg.inv(B)
        b_bar = B_inv.dot(b)
        y = B_inv.T.dot(c[basis])
        reduced_cost = c[nonbasis] - N.T.dot(y)
        print("\treduced_cost", reduced_cost)
        if np.all(reduced_cost <= 0):
            if auxiliary:
                if lp.c.dot(x) <= 0:
                    status = SolverStatus.Feasible
                else:
                    status = SolverStatus.Infeasible
            else:
                status = SolverStatus.Optimal
            break

        # select nonbasis pivot k
        k = np.where(reduced_cost > 0)[0][0].item()

        # calculate update value of nonbasis pivot k
        a_bar = B_inv.dot(N[:, k])
        if np.all(a_bar <= 0):
            assert not auxiliary
            status = SolverStatus.Unbounded
            break
        tmp = np.where(a_bar > 0, a_bar, 1.0)
        theta_ = np.where(a_bar > 0, b_bar / tmp, np.inf)
        theta = np.min(theta_)

        # select basis pivot l
        l = np.argmin(theta_)
        print("\tk", k, ", theta", theta)

        # update solution x
        x[nonbasis[k]] = theta
        x[basis] = b_bar - theta * a_bar

        # swap basis[l] and nonbasis[k]
        basis.append(nonbasis.pop(k))
        nonbasis.append(basis.pop(l))
        basis.sort()
        nonbasis.sort()

        n_iter += 1

    for var, value in zip(lp.x, x):
        var.setValue(value)

    return status

2段階最適化のコード

def two_phase_simplex(prob):
    print("------------------------------------------------------------------")
    print("  Phase 1 -- Solve auxiliary problem to find a feasible solution")
    print("------------------------------------------------------------------")
    status = simplex(prob, auxiliary=True)
    if status == SolverStatus.Infeasible:
        print("This problem is infeasible")
        return None
    else:
        print("This problem is feasible")
        print("---------------------------------------------------------------")
        print("  Phase 2 -- Solve primary problem")
        print("---------------------------------------------------------------")
        prob.show()
        status = simplex(prob, initial_values=True)
        if status == SolverStatus.Unbounded:
            print("This problem is unbounded")
            return None
        else:
            print("Optimal found!!")
            print("\toptimal objective value = ", prob.getObjectiveValue())
            for var in prob.getVariables():
                print("\t", var.name, "=", var.value())

先ほどの問題を実行してみると、こんな感じ。

import flopt

x = flopt.Variable("x", lowBound=0)
y = flopt.Variable("y", lowBound=0)
z = flopt.Variable("z", lowBound=0)

prob = flopt.Problem(name="Primary problem", sense="Minimize")
prob += -2 * x - 3 * y
prob += 3 * x + 2 * y == 10
prob += - 2 * x - 5 * y  <= -6

two_phase_simplex(prob)

Name: Primary problem
  Type         : Problem
  sense        : Minimize
  objective    : -2*x-(3*y)
  #constraints : 2
  #variables   : 2 (Continuous 2)

  C 0, name None, 3*x+2*y-10 == 0
  C 1, name None, -2*x-(5*y)+6 <= 0

------------------------------------------------------------------
  Phase 1 -- Solve auxiliary problem to find a feasible solution
------------------------------------------------------------------
Name: Auxiliary problem
  Type         : Problem
  sense        : minimize
  objective    : _artificial+0
  #constraints : 4
  #variables   : 7 (Continuous 7)

  C 0, name None, 3*x+2*y-10-artificial+__50_slack == 0
  C 1, name None, -(3*x+2*y-10)-artificial+__51_slack == 0
  C 2, name None, -2*x-(5*y)+6-artificial+__52_slack == 0
  C 3, name None, artificial-_artificial == 0

iter 0 x [10. 10.  0.  0. 20.  0.  4.] basis [0, 1, 4, 6] nonbasis [2, 3, 5] obj -10.0
    reduced_cost [ 3.  2. -1.]
    k 0 , theta 3.3333333333333335
iter 1 x [0.         0.         3.33333333 0.         0.         0.
 0.66666667] basis [1, 2, 4, 6] nonbasis [0, 3, 5] obj 0.0
    reduced_cost [-1. -0. -0.]
This problem is feasible
---------------------------------------------------------------
  Phase 2 -- Solve primary problem
---------------------------------------------------------------
Name: Primary problem
  Type         : Problem
  sense        : Minimize
  objective    : -2*x-(3*y)
  #constraints : 2
  #variables   : 2 (Continuous 2)

  C 0, name None, 3*x+2*y-10 == 0
  C 1, name None, -2*x-(5*y)+6 <= 0

Name: Primary problem
  Type         : Problem
  sense        : Minimize
  objective    : -2*x-(3*y)
  #constraints : 2
  #variables   : 3 (Continuous 3)

  C 0, name None, 3*x+2*y-10 == 0
  C 1, name None, -2*x-(5*y)+6+__53_slack == 0

iter 0 x [3.33333333 0.         0.        ] basis [0, 2] nonbasis [1] obj 6.666666666666667
    reduced_cost [1.66666667]
    k 0 , theta 5.0
iter 1 x [ 0.  5. 19.] basis [1, 2] nonbasis [0] obj 15.0
    reduced_cost [-2.5]
Optimal found!!
    optimal objective value =  -15.0
     x = 0.0
     y = 5.0

最適解として(x, y, z) = (0.0, 5.0)が見つかっていますね。

まとめ

本稿ではfloptでモデル化された線形計画問題から情報を取り出して、二段階単体法を実行するコードを紹介しました。