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\times45}$ 。设置参数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打印最优解下对应变量的值。

代码实例

整数规划

$\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}$ 对于上面这个整数规划问题，使用下面的代码求解：

import cvxpy as cp
import numpy as npa = 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\times45}$ ，代码cons = [cp.sum(x,axis = 0) <= 1]可以表示按列求和 $X_{11\times45}$ 得到的 $45$ 个结果都小于等于 $1$ 。axis=1按行求和。和np.sum用法一样，只是cp.sum是针对含有未知数的版本。
注意：cvxpy 库在进行规划求解前，会使用 DCP 规则验证目标函数的凸性，不满足 DCP 规则时会报错。只要不出现x*x这类变量和变量相乘的式子，一般不会有问题。

非线性规划

$\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}$ 可以使用下面的代码求解这个问题：

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)