Python 数学建模——cvxpy 规划求解器
文章目录
- 前言
- cvxpy 介绍
- 核心步骤
- 代码实例
- 整数规划
- 非线性规划
前言
在数学建模的过程中，难免会遇到规划问题。特别是国赛 C 题，问题往往被描述为一个非线性的复杂规划问题，在各问中调整约束条件或者目标函数，从而得到各问的答案。
诚然，当一个规划问题过于复杂时，我们更偏向于使用遗传算法、模拟退火、粒子群算法等智能算法，这类算法具有比较高的灵活度，可以面对各种各样复杂的约束条件。但是，当约束条件没那么复杂的时候，直接使用 Python 内部的规划求解器可能是一个更好的选择。
cvxpy 介绍
参考文献：Solver Features — CVXPY 1.5 documentation
cvxpy 可以求解线性规划、整数规划和非线性规划，是一个高适配度的规划求解器。还有一些求解器：
scipy.optimize.minimize：线性规划/非线性规划，求解非线性规划时需要给出一个解向量的初值，然后求解得到局部最优解。scipy.optimize.linprog：线性规划。可以参考这篇博客：数学建模学习·线性规划的Python求解_python规划求解模块-CSDN博客cvxopt.solvers：线性规划/非线性规划。
核心步骤
import cvxpy as cp。- 通过
x = cp.Variable(shape)创建形状为shape的连续型变量，比如x = cp.Variable((11,45))创建变量 X 11 × 45 X_{11\times45} X11×45。设置参数integer=True创建整数变量。 - 通过
obj = cp.Minimize(expr)设置最小化的目标函数expr。 - 设置约束条件
cons。 - 通过
prob = cp.Problem(obj,cons)构建问题模型。 - 通过
prob.solve()求解问题。
有两个参数，求解器
solver和是否显示求解的详细信息verbose。
默认solver=None，会自动选择最佳求解器。常用求解器GLPK_MI（目前我安装的求解器中，只有它和CPLEX可以解决整数规划），还可以安装其他的求解器，见参考文献。
默认verbose=False，也就是求解的时候不会打印详细信息。（这些详细信息没有什么价值）
- 通过
prob.value打印目标函数的最优解，通过x.value打印最优解下对应变量的值。
代码实例
整数规划
min z = 40 x 1 + 90 x 2 s . t . { 9 x 1 + 7 x 2 ≤ 56 7 x 1 + 20 x 2 ≥ 70 x 1 , x 2 ≥ 0 ; x 1 , x 2 ∈ Z \begin{matrix}\min{}&z=40{{x}_{1}}+90{{x}_{2}}\\&\\\mathrm{s.t.}&\left \{{\begin{matrix}9{{x}_{1}}+7{{x}_{2}}\le 56\\7{{x}_{1}}+20{{x}_{2}}\ge 70\\{{x}_{1}},{{x}_{2}}\ge 0;{{x}_{1}},{{x}_{2}}\in \mathbb Z\end{matrix}}\right .\end{matrix} mins.t.z=40x1+90x2⎩ ⎨ ⎧9x1+7x2≤567x1+20x2≥70x1,x2≥0;x1,x2∈Z 对于上面这个整数规划问题，使用下面的代码求解：
import cvxpy as cp
import numpy as np
a = np.array([[9,7],[-7,-20]])
b = np.array([56,-70])
c = np.array([40,90])
x = cp.Variable(2,integer=True)
obj = cp.Minimize(c * x)
cons = [a@x<=b, x>=0]
prob = cp.Problem(obj,cons)
prob.solve(solver='GLPK_MI',verbose=False)
print('最优值为：',prob.value)
print('最优解为：\n',x.value)
注：一般
@用于矩阵相乘。*用于标矢相乘，或者矩阵按元素对应相乘之和。
求解的结果是：
Long-step dual simplex will be used
最优值为： 350.0
最优解为：
[2. 3.]
假设构造的未知数矩阵是
X
11
×
45
X_{11\times45}
X11×45，代码cons = [cp.sum(x,axis = 0) <= 1]可以表示按列求和
X
11
×
45
X_{11\times45}
X11×45 得到的
45
45
45 个结果都小于等于
1
1
1。axis=1按行求和。和np.sum用法一样，只是cp.sum是针对含有未知数的版本。
注意：cvxpy 库在进行规划求解前，会使用 DCP 规则验证目标函数的凸性，不满足 DCP 规则时会报错。只要不出现x*x这类变量和变量相乘的式子，一般不会有问题。
非线性规划
min z = x 1 2 + x 2 2 + 3 x 3 2 + 4 x 4 2 + 2 x 5 2 − 8 x 1 − 2 x 2 − 3 x 3 − x 4 − 2 x 5 s . t . { x 1 + x 2 + x 3 + x 4 + x 5 ≤ 100 x 1 + 2 x 2 + 2 x 3 + x 4 + 6 x 5 ≤ 800 2 x 1 + x 2 + 6 x 3 ≤ 200 x 3 + x 4 + 5 x 5 ≤ 200 0 ≤ x i ≤ 99 , x i ∈ Z ( i = 1 , 2 , ⋯ , 5 ) \begin{matrix}\min{}&z={{x}_{1}^{2}}+{{x}_{2}^{2}}+3{{x}_{3}^{2}}+4{{x}_{4}^{2}}+2{{x}_{5}^{2}}-8{{x}_{1}}-2{{x}_{2}}-3{{x}_{3}}-{{x}_{4}}-2{{x}_{5}}\\&\\\mathrm{s.t.}&\left \{{\begin{matrix}{{x}_{1}}+{{x}_{2}}+{{x}_{3}}+{{x}_{4}}+{{x}_{5}}\le 100\\{{x}_{1}}+2{{x}_{2}}+2{{x}_{3}}+{{x}_{4}}+6{{x}_{5}}\le 800\\2{{x}_{1}}+{{x}_{2}}+6{{x}_{3}}\le 200\\{{x}_{3}}+{{x}_{4}}+5{{x}_{5}}\le 200\\0\le {{x}_{i}}\le 99,{{x}_{i}}\in \mathbb Z(i=1,2,\cdots ,5)\end{matrix}}\right .\end{matrix} mins.t.z=x12+x22+3x32+4x42+2x52−8x1−2x2−3x3−x4−2x5⎩ ⎨ ⎧x1+x2+x3+x4+x5≤100x1+2x2+2x3+x4+6x5≤8002x1+x2+6x3≤200x3+x4+5x5≤2000≤xi≤99,xi∈Z(i=1,2,⋯,5) 可以使用下面的代码求解这个问题：
import cvxpy as cp
import numpy as np
c1=np.array([1, 1, 3, 4, 2])
c2=np.array([-8, -2, -3, -1, -2])
a=np. array ([ [1, 1, 1, 1, 1], [1, 2, 2, 1, 6], [2, 1, 6, 0, 0],
[0, 0, 1, 1, 5]])
b=np. array ([400, 800, 200, 200])
x=cp.Variable(5,integer=True)
obj=cp.Minimize(c1*x**2+c2*x)
con=[0<=x, x<=99, a*x<=b]
prob= cp.Problem(obj, con)
prob.solve()
print("最优值为：",prob.value)
print("最优解为： \n",x.value)
然而，并不是说能够解决所有的非线性规划。涉及到 log x \log x logx或者 2 / x 2/x 2/x 之类的非线性目标函数，用目前的求解器还是无法进行。