數學模組化學習（1） 線性規劃

使用numpy+scipy編寫

經典線性規劃

經典線性規劃可被描述為下列問題:

m i n   c T x s u c h   t h a t   A u b x ≤ b u b A e q x = b e q l ≤ x ≤ u min\ c^Tx\\ such\ that\ A_{ub}x\leq b_{ub}\\ A_{eq}x=b_{eq}\\ l\leq x\leq u min cTxsuch that Aub​x≤bub​Aeq​x=beq​l≤x≤u

其中x，c，l，u都是n維向量，A是（m，n）矩陣，b是m維向量，這些式子限定了一個經典線性規劃的條件

在scipy.optimize.linprog中，其分别對應于以下參數表達：

import numpy as np
from scipy.optimize import linprog
c = np.array([-29.0, -45.0, 0.0, 0.0])
A_ub = np.array([[1.0, -1.0, -3.0, 0.0],
                 [-2.0, 3.0, 7.0, -3.0]])
b_ub = np.array([5.0, -10.0])
A_eq = np.array([[2.0, 8.0, 1.0, 0.0],
               [4.0, 4.0, 0.0, 1.0]])
b_eq = np.array([60.0, 60.0])
x0_bounds = (0, None)
x1_bounds = (0, 5.0)
x2_bounds = (-np.inf, 0.5)  # +/- np.inf can be used instead of None
 x3_bounds = (-3.0, None)
bounds = [x0_bounds, x1_bounds, x2_bounds, x3_bounds]
result = linprog(c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq, bounds=bounds)
           

在Optimize包函數中計算的結果一般傳回

OptimizeResult

，其結果整體如下圖，可以使用類似字典的操作方法取值：

典型可轉化為線性規劃問題的問題

絕對值求職

對于問題：

m i n   ∣ x 1 ∣ + . . . + ∣ x n ∣ s . t .   A x ≤ b min\ |x_1|+...+|x_n|\\ s.t.\ Ax\leq b min ∣x1​∣+...+∣xn​∣s.t. Ax≤b

可以令 u i = x i + ∣ x i ∣ 2 , v i = ∣ x i ∣ − x i 2 u_i=\frac{x_i+|x_i|}{2},v_i=\frac{|x_i|-x_i}{2} ui​=2xi​+∣xi​∣​,vi​=2∣xi​∣−xi​​，并由此有：

x i = u i − v i , ∣ x i ∣ = u i + v i x_i = u_i-v_i,|x_i|=u_i+v_i xi​=ui​−vi​,∣xi​∣=ui​+vi​

上述的線性規劃可以表示為：

m i n ∑ i = 1 n ( u i + v i ) s . t .   { A ( u − v ) ≤ b u , v ≥ 0 min\sum_{i=1}^n (u_i+v_i)\\ s.t.\ \begin{cases} A(u-v)\leq b \\ u,v\geq 0 \end{cases} mini=1∑n​(ui​+vi​)s.t. {A(u−v)≤bu,v≥0​

矩陣化表示，其中關于規劃的限制條件使用矩陣表示：

[ A , − A ] [ u v ] ≤ b \begin{bmatrix} A,-A \end{bmatrix} \begin{bmatrix} u\\v \end{bmatrix}\leq b [A,−A​][uv​]≤b

最終的期望：

[ c , c ] [ u v ] \begin{bmatrix} c,c \end{bmatrix} \begin{bmatrix} u\\v \end{bmatrix} [c,c​][uv​]

可以看出，實際上将問題轉化為了求 [ u v ] \begin{bmatrix}u\\v\end{bmatrix} [uv​]的值，是以對于如下問題：

可以轉換為以下函數：

def solve_abs_linear_prog():
    c = np.array([1, 2, 3, 4])

    A_ub = np.array([[1, -1, -1, 1],
                     [1, -1, 1, -3],
                     [1, -1, -2, 3]])
    b_ub = np.array([-2, -1, -0.5])
    res = optimize.linprog(c=np.array(c.tolist() + c.tolist()),
                           A_ub=np.column_stack((A_ub, -A_ub)),
                           b_ub=b_ub,
                           bounds=[(0, None)] * 8)
    print(res.x)
    print(round(np.array(c.tolist() + c.tolist()) @ res.x, 2))
           

結果（不得不吐槽python的數值精準度…）：