python优化包简介
以下的内容简要介绍了qpsolver库、cvxopt库以及ortools。以下是将会用到的引用代码。
import numpy as np
from qpsolvers import solve_qp
import cvxopt
from cvxopt import matrix,solvers
qpsolvers库中的solve_qp和cvxopt能够解优化问题,但是前者能够兼容numpy,后者函数参数只能使用matrix类型,否则会报错。
qpsolvers
qpsolvers中集成了许多优化算法,不仅有cvxopt,还有cvxpy、dense、osqp等等(但是似乎只能使用cvxopt?)
目前摸索出来的结论是这个只能用来进行普通的求解,假如想像matlab设置什么选项,没找到入口(但是代码里好像有?有空可以再找找)
这是qpsolvers的文档介绍
solver_qp
solver_qp是python中解决具有线性约束的二次规划的问题,数学建模中往往遇到的问题也是这类型的问题,此时就可以用这个函数去求解。
min
x
1
2
x
T
H
x
+
f
T
x
s
u
c
h
t
h
a
t
{
A
x
≤
b
,
A
e
q
x
=
b
e
q
,
l
b
≤
x
≤
u
b
\min_{x}{\frac{1}{2} x^T H x + f^T x } \\ such \ that \begin{cases} A x \leq b, \\ Aeq x = beq, \\ lb \leq x \leq ub \end{cases}
xmin21xTHx+fTxsuch that⎩
⎨
⎧Ax≤b,Aeqx=beq,lb≤x≤ub
函数
ans = solve_qp(H, f, solver)
ans = solve_qp(H, f, A, b, Aeq, beq, lb, ub, solver)
ans = solve_qp(H, f, A, b, solver, initvals, sym_proj, verbose)
- H:一般要求是对称矩阵,当不对称时需要加一个sym_proj=True的参数(好像没要求正定?)
- f:向量
- solver:要求只能选择qpsolvers.available_solvers里的选项(输出出来的话其实只有一个cvxopt,也就是只能且必须写这个)
- initvals:设置初始解,是一个字典,有x,s,y,z四个key。(在coneprog文件中)
- initvals[‘x’]是一个(n,1)的密集矩阵。
- initvals[‘s’]是一个(K,1)的密集矩阵,表示相对于圆锥C严格为正的向量。
- initvals[‘y’]是一个(p,1)的密集矩阵。
- initvals[‘z’]是一个(K,1)的密集矩阵,表示相对于圆锥C严格为正的向量。
eg.1
H=np.array([[1.,-1.],[-1.,2.]])
f=np.array([[-2.],[-6.]]).reshape((2,))
x = solve_qp(H,f,solver='cvxopt')
print("QP solution: x = {}".format(x))
A=np.array([[1.,1.],[-1.,2.],[2.,1.]])
b=np.array([[2.],[2.],[3.]]).reshape((3,))
#print(qpsolvers.available_solvers)
x = solve_qp(H, f, A,b, solver="cvxopt") #这里solver必须写,而且只能写这个,写其他会报错。
print("QP solution: x = {}".format(x))
寄!
没找到设置精度、显示过程等选项,设置初值这个会报错,因为初值是一个字典,但不知道初值应该加到哪一项(后面试用cvxopt好像找到设置初值的用法了,这个与cvxopt应该也类似,但是我没试过)
x0 = np.array([1,2])
x = solve_qp(H, f, A, b, solver='cvxopt', initvals=x0) # 会报错
cvxopt
cvxopt也是一个求解器,其实它好像已经被上面的qpsolvers集成进去了,奈何上面那个不给力,一点也不好用,文档还不多,就只能用这个了。这个的矩阵不支持ndarray,但是可以很方便地转换。也支持不少求解器,能够解决许多问题,如最小二乘、方程组求解、线性规划、二次规划、非线性凸优化等等问题。
还有,这个的matrix好像跟np数组排列方式不同,np数组是先行后列,这个是先列后行,需要注意!
这是cvxopt的官方文档
输出结果
输出结果有许多,基本都有’status’, ‘x’, ‘s’, ‘z’, ‘y’,
‘primal objective’, ‘dual objective’, ‘gap’, ‘relative gap’,
‘primal infeasibility’, ‘dual infeasibility’, ‘primal slack’,‘dual slack’
线性规划
线性规划可以使用cvxopt.solvers.lp,求解的问题形式大致如下:
m
i
n
i
m
i
z
e
c
T
x
s
u
b
j
e
c
t
t
o
G
x
≤
h
A
x
=
b
.
\begin{aligned} & minimize & c^{T}x\\ & subject \ to & Gx \leq h \\ & &Ax=b. \end{aligned}
minimizesubject tocTxGx≤hAx=b.
函数:cvxopt.solver.lp(c, G, h[, A, b[, solver[, primalstart…]]]])
- c,G,h,A,b 数据类型必须都是’d’,不能是’i’(整型)
- solver:可选None(conelp),‘glpk’,‘mosek’ (后面这俩是外部求解器,若想使用需先安装)
- primalstart:初始解,当solver为None时可选。是一个字典。
- x:是一个长度为n的向量,数据类型为’d’;
- s:是一个长度为m的向量,数据类型为’d’;
- dualstart:好像也是初始解,dual starting point。也是一个字典。
- y:是一个长度为p的向量,数据类型为’d’;
- z:是一个长度为m的向量,数据类型为’d’;
注意 当solver选项为None时,该函数使用前需满足以下条件:
- n ≥ \geq ≥ 1
- Rank(A) = p
- Rank([G; A]) = n
输出
x,s,y,z是原始和对偶最优条件的近似解,满足:
G
x
+
s
=
h
,
A
x
=
b
G
′
z
+
A
′
y
+
c
=
0
s
≥
0
,
z
≥
0
s
′
z
=
0
\begin{aligned} & Gx + s = h , Ax = b \\ & G'z + A'y + c = 0 \\ & s \geq 0 , z \geq 0 \\ & s'z = 0 \end{aligned}
Gx+s=h,Ax=bG′z+A′y+c=0s≥0,z≥0s′z=0
eg.2
m
i
n
i
m
i
z
e
−
4
x
1
−
5
x
2
s
u
b
j
e
c
t
t
o
2
x
1
+
x
2
≤
3
x
1
+
2
x
2
≤
3
x
1
≥
0
,
x
2
≥
0.
\begin{aligned} & minimize & -4x_1-5x_2\\ & subject \ to & 2x_1+x_2 \leq 3 \\ & & x_1+2x_2 \leq 3 \\ & & x_1 \geq 0 \ ,\ x_2 \geq 0. \end{aligned}
minimizesubject to−4x1−5x22x1+x2≤3x1+2x2≤3x1≥0 , x2≥0.
c = matrix([-4.,-5.]) #注意这里数后必须加点,否则它的类型不对会报错
print(c.size)
G = matrix([[2.,1.,-1.,0.],[1.,2.,0.,-1.]])
h = matrix([3.,3.,0.,0.])
print(G.typecode,h.typecode) #dtype
sol = solvers.lp(c,G,h)
print(sol['x'])
# 加上初始解的情况
x0 = matrix([10.,-1.])
s0 = matrix([1.,1.,0.7,1.]) #为啥还需要有s? 这里s每项都必须是正数,0也不行
primal = {'x':x0,'s':s0}
sol = solvers.lp(c,G,h, primalstart=primal)
print(sol['x'])
二次规划
二次规划可以使用cvxopt.solvers.qp,求解的问题形式大致如下:
m
i
n
i
m
i
z
e
1
2
x
T
P
x
+
q
T
x
s
u
b
j
e
c
t
t
o
G
x
≤
h
A
x
=
b
.
\begin{aligned} & minimize & \frac{1}{2}x^{T}Px+q^{T}x\\ & subject \ to & Gx \leq h \\ & &Ax=b. \end{aligned}
minimizesubject to21xTPx+qTxGx≤hAx=b.
函数:cvxopt.solver.qp(P, q[, G, h[, A, b[, solver[, initvals]]]])
- c,G,h,A,b 数据类型必须都是’d’,不能是’i’(整型)
- solver:可选None(conelp),‘mosek’
- primalstart:初始解,当solver为None时可选。是一个字典。
- x:是一个长度为n的向量,数据类型为’d’;
- s:是一个长度为m的向量,数据类型为’d’;
- dualstart:好像也是初始解,dual starting point。也是一个字典。
- y:是一个长度为p的向量,数据类型为’d’;
- z:是一个长度为m的向量,数据类型为’d’;
注意 该函数使用前需满足以下条件:
- P的下三角部分必须有值,而且P必须是正定的
输出
- status:有四个可能值
- ‘optimal’
- ‘unkonwn’
- ‘primal infeasible’:意味着找到了原始不可行性的证明,'x’和’s’为None,'z’和’y’是近似解。
- ‘dual infeasible’:意味着找到了双重不可行性的证明,'z’和’y’为None,'x’和’s’是近似解。
- x,s,y,z是原始和对偶最优条件的近似解,满足:
G x + s = h , A x = b P x + G ′ z + A ′ y + q = 0 s ≥ 0 , z ≥ 0 s ′ z = 0 Gx + s = h , Ax = b \\ Px + G'z + A'y + q = 0 \\ s \geq 0 , z \geq 0 \\ s'z = 0 Gx+s=h,Ax=bPx+G′z+A′y+q=0s≥0,z≥0s′z=0 - primal_objective:优化目标的最终值
- dual_objective:对偶问题的最优值
eg.3
m
i
n
i
m
i
z
e
2
x
1
2
−
4
x
2
2
+
6
x
1
x
2
−
3
x
1
−
5
x
2
s
u
b
j
e
c
t
t
o
2
x
1
+
x
2
≤
10
,
x
2
≥
−
1
x
1
−
3
x
2
=
5.
\begin{aligned} & minimize & 2x_1^2-4x_2^2+6x_1x_2-3x_1-5x_2\\ & subject \ to & 2x_1+x_2 \leq 10\ ,\ x_2 \geq -1 \\ & &x_1-3x_2 = 5. \end{aligned}
minimizesubject to2x12−4x22+6x1x2−3x1−5x22x1+x2≤10 , x2≥−1x1−3x2=5.
经过手动推导,我们发现这就是一个二次函数求最小值的问题,大概极值点在
x
2
=
−
76
44
x_2=-\frac{76}{44}
x2=−4476处,由于
x
2
≥
−
1
x_2 \geq -1
x2≥−1,所以最后结果应该是
x
2
=
−
1
,
x
1
=
2
x_2=-1\ ,\ x_1=2
x2=−1 , x1=2
# P = matrix([[2.,3.],[3.,-4.]])
# print(P) # 因为P是对称矩阵,所以行列互换没啥影响
P = matrix([[2., 3.], [3., -4.]]).T
q = matrix([-3.,-5.])
# G = matrix([[2.,1.],[0.,-1.]])
# print(G)
G = matrix([[2., 1.], [0., -1.]]).T
print(G)
h = matrix([10.,1.])
A = matrix([1.,-3.]).T
print(A.size)
b = matrix([5.])
sol = solvers.qp(P,q,G,h,A,b)
print(sol['x'])
print(sol['primal objective'])
非线性凸优化
非线性凸优化可以使用cvxopt.solvers.cp,求解的问题形式大致如下:
m
i
n
i
m
i
z
e
f
0
(
x
)
s
u
b
j
e
c
t
t
o
f
k
(
x
)
,
k
=
1
,
.
.
.
,
m
G
x
≤
h
A
x
=
b
.
\begin{aligned} & minimize & f_0(x)\\ & subject \ to & f_k(x), k=1,...,m \\ & &Gx \leq h \\ & &Ax=b. \end{aligned}
minimizesubject tof0(x)fk(x),k=1,...,mGx≤hAx=b.
函数:cvxopt.solver.qp(F[, G, h[, dims[, A, b[, kktsolver]]]])
- F是一个函数,要求能够返回如下值
- F会返回一个元组(mnl, x0),mnl是非线性不等式约束的数量,x0是定义域内的一个点。
- F(x)也会返回一个元组(f, Df),f是(mnl+1,1)数据类型为’d’的矩阵。Df是f在x点处的微分。假如x点不在f定义域内,返回None或(None,None)。
- F(x,z)会返回3个值,前两个为f, Df,最后一个是H。z也是’d’的(mnl+1,1)的矩阵,H是(n,n)的矩阵,代表海森矩阵。(但是Hij不知道具体代表什么~)
- G,h,A,b 数据类型必须都是’d’,不能是’i’(整型)
注意 该函数使用前需满足以下条件:
- 函数 f k f_k fk是凸函数且可二次微分
- 剩下那个条件说的不是人话,自己有兴趣可以看文档介绍,应该问题不大。
输出
- status:有四个可能值
- ‘optimal’:'x’为原始最优解,'snl’和’sl’是非线性和线性不等式约束中的对应松弛。'y’的解释也有,但是没看懂。
- ‘unkonwn’
其他
还有其他一些参数可调,能够决定是否显示求解过程,调整最大迭代代数,绝对精度,相对精度,容差等,具体可见这里 (其实上面的函数也有额外参数可调,不过没写出来)
eg.4
等式约束解析定心问题定义为:
m
i
n
i
m
i
z
e
−
∑
i
=
1
m
log
x
i
s
u
b
j
e
c
t
t
o
A
x
=
b
\begin{aligned} & minimize & -\sum_{i=1}^{m}{\log x_i}\\ & subject \ to & Ax=b \end{aligned}
minimizesubject to−i=1∑mlogxiAx=b
下面定义的函数acent解决了这个问题,假设它是可解决的。
from cvxopt import spdiag, log
A = matrix([1., -3.]).T
b = matrix([5.])
m,n = A.size
def F(x=None, z=None):
if x is None: return 0, matrix(1.0, (n, 1)) #?这里非线性约束没有,不应该是0.0吗,但是那样会报错
if min(x) <= 0.0: return None
f = -sum(log(x))
Df = -(x ** -1).T
if z is None: return f, Df
H = spdiag(z[0] * x ** -2)
return f, Df, H
# solvers.options['show_progress'] = False
sol = solvers.cp(F, A=A, b=b)
print(sol['x']) #这个结果最后会超过迭代次数停止
ortools
之所以把它单独列出来,是因为它不只是能够优化凸函数,而且有优化非凸函数的算法。
OR-Tools是一个用于优化的开源软件套件,用于解决车辆路径、流程、整数和线性规划以及约束编程等世界上最棘手的问题。
当前ortools提供的优化器包括:
- 约束规划
- 线性与混合整数规划
- 路径规划
- 调度规划
- 网络规划
- 装箱(对于三维装箱问题,需要用约束规划来求解)
相关文档
线性规划
相比上面的函数,这里添加求解的约束会直观很多,但是未来有可能遇到许多约束的情况,如何避免手动一个一个添加是一个问题。
编程套路
- 加载一个线性求解器的wrapper
- 实例化求解器
- 创建变量,为变量设置上下限
- 定义约束
- 定义目标函数
- 调用求解方法
- 输出求解结果
#1.加载一个线性求解器的wrapper
from ortools.linear_solver import pywraplp
#2.实例化求解器,这里使用的是Google自研的GLOP线性规划求解器
solver = pywraplp.Solver.CreateSolver('GLOP')
# 3.创建变量,为变量设置上下限
x = solver.NumVar(0, solver.infinity(), 'x')
y = solver.NumVar(0, solver.infinity(), 'y')
print('Number of variables =', solver.NumVariables())
#4.定义约束
# Constraint 0: x + 2y <= 14.
solver.Add(x + 2 * y <= 14.0)
# Constraint 1: 3x - y >= 0.
solver.Add(3 * x - y >= 0.0)
# Constraint 2: x - y <= 2.
solver.Add(x - y <= 2.0)
print('Number of constraints =', solver.NumConstraints())
#5.定义目标函数
solver.Maximize(3 * x + 4 * y)
#6.调用求解方法
status = solver.Solve()
#7.输出求解结果
if status == pywraplp.Solver.OPTIMAL:
print('Solution:')
print('Objective value =', solver.Objective().Value())
print('x =', x.solution_value())
print('y =', y.solution_value())
else:
print('The problem does not have an optimal solution.')
print('\nAdvanced usage:')
print('Problem solved in %f milliseconds' % solver.wall_time())
print('Problem solved in %d iterations' % solver.iterations())
上面第二个博客里对ortools介绍挺全的,这里就不再赘述了。