線形計画問題

線形計画問題の解法の一つである「シンプレックス法」をPythonで実装してみました。

def LP(c, A, b, xb_idx):
    xn_idx = np.arange(A.shape[1])
    xn_idx = np.delete(xn_idx, xb_idx)

    while True:
        B = A[:, xb_idx]
        N = A[:, xn_idx]
        B_inv = np.linalg.inv(B) # Bの逆行列
        b_bar = np.dot(B_inv, b) # bバー。(xB, xN) = (bバー, 0)
        if np.any(b_bar < 0) :
            print('initial basic solution isn\'t solution!')
            return 0
        cB = c[xb_idx]
        cN = c[xn_idx]
        pi = np.dot(B_inv.T, cB) #　シンプレックス乗数
        rcf = cN - np.dot(N.T, pi) # 相対コスト係数
        
        if np.all(rcf >= 0) :
            solution = np.zeros_like(c, dtype=float)
            solution[xb_idx] = b_bar
            return solution
        
        for i in range(len(rcf)):
            if rcf[i] < 0 :
                choose = xn_idx[i]
                choose_idx = i
                break
            
        y = np.dot(B_inv, A[:, choose])
        y = y.reshape(len(y), 1)
        
        if np.all(y <= 0) :
            print('Unbounded')
            return 0
        
        theta_opt = b_bar / y
        
        theta = np.min(theta_opt)
        dropout = np.argmin(theta_opt)
        
        temp = xn_idx[choose_idx]
        xn_idx[choose_idx] = xb_idx[dropout]
        xb_idx[dropout] = temp
        xb_idx = sorted(xb_idx)
        xn_idx = sorted(xn_idx)

具体的にはこう使います。

\begin{align}
&minimize &-2x_1-3x_2 \\
&\text{subject to} &4x_1+3x_2+x_3&=120\\
&&x_1+2x_2+x_4&=60\\
&&x_i \geq 0 &&\forall{i = 1, 2, 3, 4}
\end{align}

こんな時は

A = np.array([[4,3,1,0],[1,2,0,1]])
b = np.array([[120,60]]).T
c = np.array([[-2,-3,0,0]]).T
xb_idx = [2,3]
x = LP(c, A, b, xb_idx)
print('Optimum solution')
print(x)
print('objective function -> ',np.round(np.dot(c.T, x)[0][0]))
print('constraints -> ', np.all(np.round(np.dot(A, x) - b, 7) == 0))

出力は

Optimum solution
[[12.]
[24.]
[ 0.]
[ 0.]]
objective function -> -96.0
constraints -> True