python优化算法包

最新推荐文章于 2024-04-24 13:35:29 发布

F77_CSDN

最新推荐文章于 2024-04-24 13:35:29 发布

阅读量2.9k

点赞数 2

分类专栏：数学建模文章标签： python

本文链接：https://blog.csdn.net/F77_CSDN/article/details/127017532

版权

数学建模专栏收录该内容

1 篇文章 0 订阅

订阅专栏

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}{\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}$

函数
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，求解的问题形式大致如下：
$\begin{aligned} & minimize & c^{T}x\\ & subject \ to & Gx \leq h \\ & &Ax=b. \end{aligned}$

函数：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是原始和对偶最优条件的近似解，满足：
$\begin{aligned} & Gx + s = h , Ax = b \\ & G'z + A'y + c = 0 \\ & s \geq 0 , z \geq 0 \\ & s'z = 0 \end{aligned}$

eg.2
$\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}$

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，求解的问题形式大致如下：
$\begin{aligned} & minimize & \frac{1}{2}x^{T}Px+q^{T}x\\ & subject \ to & Gx \leq h \\ & &Ax=b. \end{aligned}$

函数：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是原始和对偶最优条件的近似解，满足：
$\\ Px + G'z + A'y + q = 0 \\ s \geq 0 , z \geq 0 \\ s'z = 0$
primal_objective:优化目标的最终值
dual_objective:对偶问题的最优值

eg.3
$\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}$
经过手动推导，我们发现这就是一个二次函数求最小值的问题，大概极值点在 $x_2=-\frac{76}{44}$ 处，由于 $x_2 \geq -1$ ,所以最后结果应该是 $x_2=-1\ ,\ x_1=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，求解的问题形式大致如下：
$\begin{aligned} & minimize & f_0(x)\\ & subject \ to & f_k(x), k=1,...,m \\ & &Gx \leq h \\ & &Ax=b. \end{aligned}$
函数：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$ 是凸函数且可二次微分
剩下那个条件说的不是人话，自己有兴趣可以看文档介绍，应该问题不大。

输出

status:有四个可能值
- ‘optimal’:'x’为原始最优解,'snl’和’sl’是非线性和线性不等式约束中的对应松弛。'y’的解释也有，但是没看懂。
- ‘unkonwn’

其他
还有其他一些参数可调，能够决定是否显示求解过程，调整最大迭代代数，绝对精度，相对精度，容差等，具体可见这里（其实上面的函数也有额外参数可调，不过没写出来）

eg.4
等式约束解析定心问题定义为:
$\begin{aligned} & minimize & -\sum_{i=1}^{m}{\log x_i}\\ & subject \ to & Ax=b \end{aligned}$
下面定义的函数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提供的优化器包括：

约束规划
线性与混合整数规划
路径规划
调度规划
网络规划
装箱(对于三维装箱问题，需要用约束规划来求解)

相关文档

官方文档
重载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())