线性规划模型形式
线性规划模型的一般形式为
max
(
min
)
z
=
∑
j
=
1
n
c
j
x
j
;
(
1
)
s.t.
{
∑
j
=
1
n
a
i
j
x
j
≤
(
≥
,
=
)
b
i
,
i
=
1
,
2
,
⋯
,
m
,
x
j
≥
0
,
j
=
1
,
2
,
⋯
,
n
.
(
2
)
\begin{gathered} \max (\min ) z=\sum_{j=1}^{n} c_{j} x_{j} ; \quad (1) \\ \text { s.t. } \begin{cases} \sum_{j=1}^{n} a_{i j} x_{j} \leq(\geq,=) b_{i}, i=1,2, \cdots, m, \\ x_{j} \geq 0, \quad j=1,2, \cdots, n . \end{cases} \quad (2) \end{gathered}
max(min)z=j=1∑ncjxj;(1) s.t. {∑j=1naijxj≤(≥,=)bi,i=1,2,⋯,m,xj≥0,j=1,2,⋯,n.(2)
或矩阵形式
max
(
min
)
z
=
c
T
x
;
(
1
)
s.t.
{
A
x
≤
(
≥
,
=
)
b
,
x
≥
0.
(
2
)
\begin{gathered} &\max (\min ) z=c^{T} x; \quad (1) \\ &\text { s.t. }\left\{\begin{array}{l} A x \leq(\geq,=) b, \\ x \geq 0. \end{array}\right. \quad (2) \end{gathered}
max(min)z=cTx;(1) s.t. {Ax≤(≥,=)b,x≥0.(2)
其中,式 (1) 称为目标函数,式 (2) 称为约束条件。
线性规划问题的Python求解
1. scipy.optimize.linprog
标准型式:
min
z
=
c
T
x
s.t.
{
A
⋅
x
≤
b
A
e
q
⋅
x
=
b
e
q
L
b
≤
x
≤
U
b
\begin{gathered} \min \quad z=c^{T} x \\ \text { s.t. }\left\{\begin{array}{l} A \cdot x \leq b \\ A e q \cdot x=b e q \\ L b \leq x \leq U b \end{array}\right. \end{gathered}
minz=cTx s.t. ⎩⎨⎧A⋅x≤bAeq⋅x=beqLb≤x≤Ub
调用方式:
from scipy.optimize import linprog
res=linprog(c, A=None, b=None, Aeq=None, beq=None, bounds=None, method='simplex')
c, A, b, Aeq, beq 分别为上述标准型式中的目标向量, 不等号约束, 与等号约束.
bounds 是决策向量的下界和上界所组成的
n
n
n 个元素的元组. 如
−
∞
<
x
1
<
+
∞
- \infty < {x_1} < + \infty
−∞<x1<+∞,
−
3
⩽
x
2
<
+
∞
- 3 \leqslant {x_2} < + \infty
−3⩽x2<+∞, 则 bounds=((None,None),(-3,None))
# Lb, Ub 转化为 bounds
bounds=tuple(zip(Lb, Ub))
默认每个决策变量下界为 0 0 0, 上界为 + ∞ + \infty +∞
返回值:
print(res.fun) #最优值 (目标函数最小值)
print(res.x ) #最优解
2. cvxopt.solvers
标准型:
min
z
=
c
T
x
s.t.
{
A
⋅
x
≤
b
A
e
q
⋅
x
=
b
e
q
\begin{gathered} &\min \quad z=c^{T} x \\ &\text { s.t. }\left\{\begin{array}{l} A \cdot x \leq b \\ A e q \cdot x=b e q \end{array}\right. \end{gathered}
minz=cTx s.t. {A⋅x≤bAeq⋅x=beq
调用方式:
import numpy as np
from cvxopt import matrix, solvers
sol=solvers.lp(c,A,b,Aeq,beq)
#c,A,b,Aeq,beq 为 cvxopt.matrix 矩阵
'''
可用 cvxopt.matrix(array,dims) 把array按照dims重新排成矩阵
省略dims: 如果array为np.array, 则为其原本形式;
如果array为list, 则认为 list 中用逗号分隔开的为一"列"
'''
返回值:
print(sol['x']) #最优解
print(sol['primal objective']) #最优值
注:
- 在程序中虽然没有直接使用 NumPy 库中的函数,也要加载
- 数据如果全部为整型数据,也必须写成浮点型数据
3. cvxpy
调用方式:
import cvxpy as cp
import numpy as np
import pandas as pd
x=cp.Variable(shape=(), name=None, var_id=None, **kwargs)
#shape: 表示形状, 可以使用元组表示矩阵
#name: 变量名字, 可以使用字符串
#例如 x=cp.Variable((m,n)) # m*n 的矩阵
#或 x=cp.Variable(shape=(3,3), name='cov', symmetric=True)
'''数域
boolean 布尔型, integer 整数型
x=cp.Variable(shape=(1),boolean=True)
y=cp.Variable(shape=(1),integer=True)
neg 负数, nonneg 非负
x=cp.Variable(shape=(1),nonneg=True)
pos 正数, nonpos 非正
x=cp.Variable(shape=(1),nonpos=True)
complex 复数, imag 虚数
'''
obj=cp.Minimize(f) #构造目标函数
#可选 cp.Minimize, cp.Maximize
'''
cvxpy 的函数
使用 * 进行矩阵标量和向量标量乘法。
使用 @ 进行矩阵-矩阵和矩阵-向量乘法。
使用 cp.multiply 进行元素乘法
'''
con=[] #构造约束条件的列表
prob=cp.Problem(obj,con) #构造模型
prob.solve(solver='GLPK_MI',verbose=True) #求解模型, verbose=True 表明显示具体求解信息
#显示已安装的求解器: cp.installed_solvers()
返回值:
print(prob.value) # 最优值
print(x.value) # 最优解
示例:
min
∑
i
=
1
m
∑
j
=
1
n
c
i
j
x
i
j
,
s.t.
{
∑
j
=
1
n
x
i
j
≤
b
⃗
,
i
=
1
,
2
,
⋯
,
m
,
∑
i
=
1
m
x
i
j
=
b
e
q
⃗
,
j
=
1
,
2
,
⋯
,
n
,
L
b
≤
X
≤
U
b
,
X
=
(
x
i
j
)
,
i
=
1
,
2
,
⋯
,
m
;
j
=
1
,
2
,
⋯
,
n
.
\begin{gathered} &\min \quad \sum_{i=1}^{m} \sum_{j=1}^{n} c_{i j} x_{i j}, \\ &\text { s.t. }\left\{\begin{array}{l} \displaystyle{\sum_{j=1}^{n} x_{i j} \leq \vec{b}, \quad i=1,2, \cdots, m,} \\ \displaystyle{\sum_{i=1}^{m} x_{i j}= \vec{beq}, \quad j=1,2, \cdots, n,} \\ \displaystyle{\mathbf{Lb} \leq \mathbf{X} \leq \mathbf{Ub}, \quad \mathbf{X}=(x_{ij}), i=1,2, \cdots, m ; j=1,2, \cdots, n .} \end{array}\right. \end{gathered}
mini=1∑mj=1∑ncijxij, s.t. ⎩⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎧j=1∑nxij≤b,i=1,2,⋯,m,i=1∑mxij=beq,j=1,2,⋯,n,Lb≤X≤Ub,X=(xij),i=1,2,⋯,m;j=1,2,⋯,n.
import cvxpy as cp
import numpy as np
import pandas as pd
x=cp.Variable((m,n)) #构造m*n维变量
obj=cp.Minimize(cp.sum(cp.multiply(c,x))) #构造目标函数
con=[cp.sum(x,axis=1,keepdims=True)<=b,cp.sum(x,axis=0,keepdims=True)==beq,Lb<=x<=Ub] #构造约束条件
prob=cp.Problem(obj,con) #构造模型
prob.solve(solver='GLPK_MI',verbose=True) #求解模型
print(prob.value) # 最优值
print(x.value) # 最优解