使用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 Aubx≤bubAeqx=beql≤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的数值精准度…):