《数学建模简明教程--基于python》学习笔记-第六章-线性规划模型-课后习题解答

关于cvxpy的坑你要知道!!!

• 安装坑(求解器坑)
  • 见0.1
• 官方文档坑
  
  • 官方文档的另一种打开方式，这应该是你想要的:cvxpy官方文档基本函数目录
• 非负约束条件坑
  • 非负约束条件可以在定义决策变量时直接添加，也可以在约束条件中添加。
  • 但是整数的非负约束一定要记得在约束条件中添加！
• 目标函数（obj）位置坑
  • obj应该在你定义完相关问题数据和变量之后定义。
  • 约束条件也最好在目标函数（obj）定义好之后再定义。
    

0 准备工作

0.1 安装cvxpy库

# !pip install cvxpy
# !pip install cvxopt 

• 安装cvxpy的同时，建议安装cvxopt库。它提供GLPK_MI等求解器,可用于求解混合整数规划。
• 如果你使用jupyter notebook编程,安装该库后需要重启才生效！ 注意：要重启，重载是没有用的

0.2 导入库

import cvxpy as cp
import numpy as np
from numpy import array as ar
import matplotlib.pyplot as plt
# 相关设置
plt.rc('font', family='SimHei')  # 显示中文字体
plt.rc('axes',unicode_minus=False)  # 显示负号
np.set_printoptions(precision=3, suppress=True)  # 设置print()时的小数位数为3，并关闭科学计数法
np.random.seed(2022)  # 设置随机数种子

先纵览一下线性规划的题目。

1 混合饲料的比例问题

题1：

# 用矩阵形式定义模型的各个部分
c = ar([0.3, 0.18])           # 定义目标向量/价值系数
a = ar([[1,0], [1,1],[1,1]])  # 定义约束矩阵/技术系数
b = ar([100, 6000, 5000])     # 定义等式右端项，可以表示有限的资源

• 关于c,a,b的具体含义和由来请往后看到模型部分。

# 定义决策变量
x = cp.Variable(2, nonneg=True)  # nonneg：非负

• 通过调用cp库的Variable类来生成一个实例x
• "2"表示这是一个2维的向量，包括两个决策变量
• x[0]和x[1]分别表示动物和谷物饲料的数量,和c中的价值系数对应
• "nonneg=True"表示非负数
• 你可能还会想知道Variable除了nonneg还有哪些类型,请运行
  x.attributes

# 定义最小化目标函数
obj = cp.Minimize(c@x)  # @:矩阵乘法

# 定义约束条件（这里我用线性方程组的形式）
cons = [a[0]@x >= b[0],  # 动物饲料要大于100=1000*0.5*0.2
        a[1]@x <= b[1],  # 总的饲料需求量应该小于供应量 
        a[2]@x >= b[2]]  # 总的饲料需求量应该多于最低饲料需求量5000=1000*0.5

# 把这个obj和cons打包成一个Problem
prob = cp.Problem(obj, cons)

# 求解这个线性规划问题
prob.solve(solver='GLPK_MI')  # 求解问题 solver:求解器
# GLPK (GNU Linear Programming Kit:GNU线性编程工具)用于建立线性规划LP和混合型整数规划MIP问题的建模语言，并对模型进行最优化求解。 

# 打印求解结果  
print('解的类型',prob.status)
print('最优解x*为',x.value)
print('最优值z*为',prob.value)

解的类型 optimal
最优解x*为 [ 100. 4900.]
最优值z*为 912.0

• 也可以输入prob.solution查看以上信息。
• 注意x的非负约束可以在定义x时添加，也可以在约束条件中添加
  • 如果不添加，结果可能出现无界解-inf!!!
• 关于"solver:求解器"，'Variable变量类型"以及其他方面的简单介绍，可参考b站cvxpy专栏的《cvxpy 常用功能汇总》
  
  

2 最小化仓库租借费用

题2:

2.1 问题分析

首先可以考察一下每个月的可租赁期限是否都一样，这就有点像强化学习里考察智能体在每个状态的动作空间是否一样，这里显然是不一样的。

在一月份，你可以采取租赁期限为1-4个月的合同，在二月份，你可以采取租赁期限为1-3个月的合同，而没必要租4个月，因为就这题目而言我们只考虑这四个月,三月和四月也是一样的道理。

办理时间\租赁期限	1月	2月	3月	4月
1月	✔	✔	✔	✔
2月	✔	✔	✔	×
3月	✔	✔	×	×
4月	✔	×	×	×

2.2 模型假设

（1）租金的计量采取收付实现制；
（2）在签订合同时交钱 ；

2.3 符号说明

$x_{ij}$: 表示在第$i$个月租了期限为$j$个月的仓库；

2.4 模型建立

根据符号说明，相应的将x定义成一个矩阵，i表示行，j表示列。

x = cp.Variable((4,4),nonneg=True)  
c = np.array([[28]*4, [45]*3+[10000], [60]*2+[10000]*2, [73]+[10000]*3]).T  # 10000表示很大的数

• cp.Variable有很多和np.array一样的函数，但是像完全像用np.array一样用cp.Variable难免会出错，但是我在百度等无果，也看不懂官方文档的组织架构，只能首先就是大胆尝试，还未遇到什么问题，
• 但你输入cp.Variable的实例时不像np.array会返回具体形式，而是会返回一个expression，看不到里面,嗨。

• 为了弄个明白，于是我决定在pycharm里直接查看源码，看个究竟。结果在看到这里时，它又提示我去看官方文档 cvxpy基本函数目录，简直了，这就是我要找的啊！想必是东方那回归本源的神秘力量在指引我吧~~呜呜…

obj = cp.Minimize(cp.sum(cp.multiply(c,x)))  # cp.multiply表示对应元素相乘

• obj的位置请注意一定要放在相关变量的后面，否则
  一、如果你没定义相关变量，报错了,这还算运气好的；
  二、坏的是它用了你前面或其他不知道哪里的变量，那你就可能求解得到一个inf，反复运行还找不出问题QAQ

cons = [cp.sum(x[0,:]) >= 1500,
       sum([cp.sum(x[1,:]), sum(x[0,1:])]) >= 1000,  
       sum([cp.sum(x[2,:]), sum(x[0,2:]), sum(x[1,1:])]) >= 1500,
       sum([cp.sum(x[3,:]), sum(x[0,3:]), sum(x[1,2:]), sum(x[2,1:])]) >= 1200] 
prob = cp.Problem(obj, cons)  # 生成问题实例

prob.solve(solver='GLPK_MI')  #  求解问题
print('最优解x*为\n',x.value)
print('最优值z*为',prob.value)

最优解x*为
 [[ 300.   -0.   -0. 1200.]
 [  -0.   -0.   -0.   -0.]
 [ 300.   -0.   -0.   -0.]
 [  -0.   -0.   -0.   -0.]]
最优值z*为 104400.0

3 护士值班问题

题3：

• 模型可参考书本p125。

# 问题数据
c = ar([0.3, 0.18])  # 定义目标向量
a = ar([[1,0], [1,1],[1,1]])  # 定义约束矩阵
b = ar([60, 70, 60, 50, 20, 30])  # 定义右端项 通常表示有限的资源

# 问题建模
x = cp.Variable(6, integer=True)  # 定义三个决策变量 （生成一个实例）
obj = cp.Minimize(cp.sum(x))
cons = [x[0]  >= b[0],
        x[1] + x[0] >= b[1],
        x[2] + x[1] >= b[2],
        x[3] + x[2] >= b[3],
        x[4] + x[3] >= b[4],
        x[5] + x[4] >= b[5],
        x>=0]
prob = cp.Problem(obj, cons)

# 问题求解
prob.solve(solver='GLPK_MI')  # 求解问题的 solver求解器
print('最优解的状态为:',prob.status)
print('最优解为:',x.value)
print('最优值为:',prob.value)

最优解的状态为: optimal
最优解为: [60. 10. 50.  0. 30.  0.]
最优值为: 150.0

4 销售代理点覆盖范围问题

题4:

• 如果一个销售代理点可以覆盖周围多个区域的话，那就有一点冰汽时代的味道了~

4.1 问题分析

题目要求在尽可能多覆盖学生的情况下，确定一个最优的销售代理点建立方案。各个区域间有相邻和不相邻，这种关系可以抽像成一个邻接矩阵。并且这是一个连通图，不存在孤立点。此外各个区域的人数有多有少。

4.2 符号说明

$A$: 邻接矩阵;
$x$: 决策变量向量,$x_i=0$或$1$。0表示该区域未被覆盖，1表示被覆盖；
$c$: 价值系数向量,$c_i$表示该区域的学生数；
$\theta$: 零向量；

4.3 建模建立

该问题的目标函数是覆盖的大学生总数 $z=c^{T}X$。
决策变量$x_i$还应该满足以下约束条件:

首先,一共有2个销售点，那就意味着可以覆盖4个区域，因此有${\sum }_{i=1}^7{x}_i=4$。

其次，一个销售代理点可以向邻接的一个区域销售，因为是一对一的关系，所以区域和代理点的身份可以看作是可以互换的。

所以当有一个区域被销售覆盖时，那这个区域周围必然有另一个区域也被销售覆盖。也就是覆盖区域应该两个两个连着出现。 即

$x-Ax\le \theta$

综上所述，该问题的线性规划模型为

$$minz=c^{T}X$$

$$s.t.\{\begin{array}{l}{\sum }_{i=1}^7{x}_i=4,\\ x-Ax\le \theta ,\\ x_i=0或1,\end{array}$$

4.4 模型求解

全部计算的python程序如下:

1.问题数据

根据上图中的标号构建邻接矩阵。

a = np.array([[0,1,1,0,0,0,0],  # 构建邻接矩阵
             [1,0,1,1,0,0,0],
             [1,1,0,1,0,0,0],
             [0,1,1,0,1,1,1],
             [0,1,0,1,0,1,1],
             [0,0,0,1,1,0,1],
             [0,0,0,1,0,1,0]])
c = ar([34,29,42,21,56,18,71]) # 各地学生人数

2.问题构建

x = cp.Variable(7, boolean=True)  # boolean: 布尔型01变量
obj = cp.Maximize(x@s)
cons = [cp.sum(x) == 4,
        x - a @ x <= 0]
prob = cp.Problem(obj, cons)

3.问题求解

prob.solve(solver='GLPK_MI')  # 求解问题的 solver求解器
print('解的类型',prob.status)
print('最优解为x*:',x.value)
print('最优值为z*:',prob.value)

解的类型 optimal
最优解为x*: [0. 0. 1. 1. 1. 0. 1.]
最优值为z*: 190.0

5 客车接送问题

题5：

符号说明
$a_{ij}$: 表示让j号客车前往i号地区的费用,i=1,2,3,4,5,j=1,2,…8；
$x_{ij}$: $x_{ij}$=0或1，1表示让j号客车前往i号地区；

问题数据

a = ar([[125, 100, 95,128,108, 97,113, 130],
        [87,121,99,96,85,89,118,97],
        [108,116,110,89,93,87,92,129],
        [116,120,86,118,100,96,85,105],
        [92,87,90,92,102,88,99,88]])

问题构建

x = cp.Variable((5,8), boolean=True)  # x等于0表示
obj = cp.Minimize(cp.sum(cp.multiply(x,a)))
cons= [cp.sum(x, axis = 0) <= 1,  # 按列求和，表示每一辆车只能去一个地区，或者一个都不去。
       cp.sum(x) == 5,  # 一共只派遣五辆车
       cp.sum(x, axis = 1) ==1]  # 按行求和，表示每个地区要有一辆汽车，不多不少
prob = cp.Problem(obj, cons)

问题求解

prob.solve(solver='GLPK_MI')  # 求解问题的 solver求解器
print('解的类型',prob.status)
print('最优解为x*:',x.value)
print('最优值为z*:',prob.value)

解的类型 optimal
最优解为x*: [[0. 0. 1. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 1. 0. 0. 0.]
 [0. 0. 0. 0. 0. 1. 0. 0.]
 [0. 0. 0. 0. 0. 0. 1. 0.]
 [0. 1. 0. 0. 0. 0. 0. 0.]]
最优值为z*: 439.0

6 钢管下料问题

• 与书本例题一摸一样，所以在这不做多余的介绍。

# d是7种切割模式对应的切割19米长钢管得到的各种长度短钢管数量，以及各种模式余料长度（最后一列）
d = ar([4,0,0,3,
3,1,0,1,
2,0,1,3,
1,2,0,3,
1,1,1,1,
0,3,0,1,
0,0,2,3]).reshape(-1,4)
a = d[:,:3].T  # 各种模式切得的4，6，8米长度的钢管数量
r = d[:,-1]   # 各种模式的余料长度

两种目标函数都可以，二选一。

x = cp.Variable(7, integer=True)
b = np.array([50,20,15])  # 每种钢管需求数量
obj1 = cp.Minimize(r@x)  # 目标函数一：让余料和最小  
con = [x>=0, a@x>=b]
prob1 = cp.Problem(obj1, con)
prob1.solve(solver='GLPK_MI')

print('解的类型',prob.status)
print('最优解为x*:',x.value)
print('最优值为z*:',prob1.value)
print('余料的长度：',r@x.value)
print('三种短钢管的数量为',a@x.value)

解的类型 optimal
最优解为x*: [ 0. 12.  0.  0. 15.  0.  0.]
最优值为z*: 27.0
余料的长度： 27.0
三种短钢管的数量为 [51. 27. 15.]

obj2 = cp.Minimize(sum(x))  # 目标函数2：让钢管数量最少 
prob2 = cp.Problem(obj2, con)
prob2.solve(solver='GLPK_MI')

print('解的类型',prob.status)
print('最优解为x*:',x.value)
print('最优值为z*:',prob2.value)
print('余料的长度：',r@x.value)
print('三种短钢管的数量为',a@x.value)

解的类型 optimal
最优解为x*: [ 5.  5.  0.  0. 15.  0.  0.]
最优值为z*: 25.0
余料的长度： 35.0
三种短钢管的数量为 [50. 20. 15.]

7 1998年全国大学生数学建模竞赛A题（考虑投资阈值）

图7:

• 投资阈值费用： 购买一种金融资产最少要付的交易费，不投不用。
• 简化的模型（没有考虑投资阈值）在书本中已经给出，而这题要我们考虑阈值，所以我们在原来的基础上修改模型。

7.1 符号说明

$p$: 费率 ；
$r$: 平均收益率；
$q$: 风险损失率；
$x$: 决策变量，$x_i$表示银行存款以及每种金融资产投资的数量,i=0,1,2,3,4, 注意$x_5$表示的是总体投资风险；

用以下符号来实现交易费的三段函数线性化，

$u$: 需要缴纳的4种金融资产的费用向量;
$m$: 已知的4种金融资产的投资阈值费用向量;
$v$: 表示4种金融资产的投资阈值费用向量且分量可为0（这取决于是否投资该项资产）;
$n$: 01变量,用于判断$v_i$是否应该为0;
$y$: 01变量,用于选择$v_i$和$p_i{x}_i$的较大者;
$Ui$: 第i种资产的投资阈值；

问题数据

r = ar([0.05,0.28,0.21,0.23,0.25])  # 平均收益率  
q = ar([0,0.025,0.015,0.055,0.026])  # 风险损失率 
p = ar([0,0.01,0.02,0.045,0.065])  # 费率
U = ar([103,198,52,40])  # 投资阈值 当投资额度大于Ui时按照pixi算费用
m = np.multiply(p[1:], U)  # 计算投资阈值费用

注意：第一列表示银行存款，所以他的风险损失率和费率都为0。

7.2 模型建立

在购买$S_i$资产时所产生的交易费用为

$$f_i(x_i)\{\begin{array}{ll}0,& x_i=0,\\ m_i,& 0<x_i<U_i,\\ p_i{x}_i,& x_i\ge U_i,\end{array}$$

首先将这个三段函数约束转化为$u_i=max(v_i,p_i{x}_i)$和$v_i=0$或者$m_i$（$v_i$当且仅当$p_i{x}_i=0$时可以等于0），
然后继续线性化$u_i$和$v_i$这两个部分即可实现这个三段函数的线性化。则模型线性化为

$$min(w\underset{0\le i\le 4}{\max }(q_i{x}_i)-(1-w)(\sum _{i=0}^4{r}_i{x}_i-\sum _{i=0}^4{c}_i))$$

$$s.t.\{\begin{array}{ll}{q}_i{x}_i\le x_5,& i=0,1,...,4.\\ {\sum }_{i=0}^4{x}_i+{\sum }_{i=0}^3{u}_i=10000& \\ x_i\ge 0,& i=0,1,...,4.\\ x_i\le 10000n,& i=1,2,3,4.\\ n=0或1& \\ m_i-m_i(1-n_i)\le v_i\le m_i-m_i(1-n_i),& i=1,2,3,4\\ u_i\le 10000y_i+v_i,& i=1,2,3,4,\\ u_i\ge v_i,& i=1,2,3,4,\\ p_i{x}_i\le u_i\le p_i{x}_i+10000(1-y_i),& i=1,2,3,4,\\ y=0或1& \end{array}$$

模型的前三行差不多就是原来的简化模型，后面我通过4,5,6三行约束线性化$v_i$,通过最后4行约束线性化$u_i$。

线性化方法可参考该博文或书本

7.3 模型求解

# 根据构建的模型定义相关变量 
# 注意 这里用向量的形式表示各变量
x = cp.Variable(6, nonneg=True)
u = cp.Variable(4, nonneg=True)
n = cp.Variable(4, boolean=True)
v = cp.Variable(4, nonneg=True)
y = cp.Variable(4, boolean=True)

# 构建目标函数
obj = cp.Minimize(0.5 * x[5] - 0.5 * (r@x[:-1]-cp.sum(u)))

这里一定要注意obj的位置要放在相关变量赋值后面，不然会用你上一次变量生成obj这就会出现问题。我一直显示-inf，一直到换了obj的位置才得以纠正。

# 模型构建
cons = [cp.sum(x[:-1]) + cp.sum(u) == 10000,
        cp.multiply(q,x[0:5])<=x[5],  
        x[1:5] <= 10000*n,
        v <= m - cp.multiply(m, 1-n),  # 这里用向量的形式表示部分模型
        v >= m - cp.multiply(m, 1-n),
        u <= 10000*y + v,
        u >= v,
        u >= cp.multiply(p[1:5],x[1:5]),
        u <= cp.multiply(p[1:5],x[1:5]) + 10000*(1-y)]
prob = cp.Problem(obj, cons)

# 模型求解
prob.solve(solver='GLPK_MI')
print('解的类型',prob.status)
print('最优解为x*:',x.value)
print('最优值为z*:',obj.value)

解的类型 optimal
最优解为x*: [   0.    9900.99    -0.      -0.      -0.     247.525]
最优值为z*: -1212.871287128713

7.4 对不同w(风险偏好系数)做进一步分析

revenue = []
risk = []
objs = []
for w in np.linspace(0,1):
    obj = cp.Minimize(w * x[5] - (1-w) * (r@x[:-1]-cp.sum(u))) 
    cons = [cp.sum(x[:-1]) + cp.sum(u) == 10000,
        cp.multiply(q,x[0:5])<=x[5],  
        x[1:5] <= 10000*n,
        v <= m - cp.multiply(m, 1-n),  # 这里用向量的形式表示部分模型
        v >= m - cp.multiply(m, 1-n),
        u <= 10000*y + v,
        u >= v,
        u >= cp.multiply(p[1:5],x[1:5]),
        u <= cp.multiply(p[1:5],x[1:5]) + 10000*(1-y)]
    prob = cp.Problem(obj, cons)
    prob.solve(solver='GLPK_MI')
    objs.append(obj.value)
    risk.append(x.value[5])
    revenue.append(r@x.value[:5])

# 把risk风险作为自变量， 收益和风险偏好系数作为因变量
# 生成用于绘图的变量
x = risk
y1 = np.linspace(0,1)  # w:风险偏好系数
y2 = revenue

fig = plt.figure(figsize=(10, 4))
ax = fig.add_subplot(1,1,1)
# 绘制左边 
h1 =ax.scatter(x, y1, facecolors="none", edgecolors="b",label='风险偏好系数w')
ax.set_ylabel('风险偏好系数w')
ax.set_title("风险与收益对应关系图",y=-0.25)  
ax.set_xlabel('风险损失（元）')
plt.grid()
ax.legend(bbox_to_anchor=(1, 0.9)) # (1,1)右上角
# 绘制右边
ax2 = ax.twinx()
h2 = ax2.scatter(x,y2,c='r',label='收益（元）')
ax2.set_ylabel('收益（元）') 
# # 相关设置
ax2.spines['left'].set_color('b')
ax2.spines['right'].set_color('r')
ax.tick_params(axis='y', colors='b')
ax2.tick_params(axis='y', colors='r')
ax.yaxis.label.set_color('b')
ax2.yaxis.label.set_color('r')
ax2.legend()
ax2.legend(bbox_to_anchor=(0.18, 0.9)) # (1,1)右上角

plt.plot(risk,objs)
plt.text(x=10,y=100,s='目标函数值曲线')

Text(10, 100, '目标函数值曲线')

从上图可以看到随着风险损失的上升，收益在持续上升，而风险偏好系数在下降。此外，二者都呈现间断的跳跃形式，而目标函数值则呈现一个相对缓和的变化，其实就是w自身的变化导致的。

np.diff(risk).argmin()  # np.diff(risk)用后面的减去前面的 找出减少数最小的，也就是落差最大的地方

37

print(risk[37],risk[38])  # 绘图时先绘制92 
print(y1[37],y1[38])  # y1是风险偏好系数w 0~1的列表

247.52475247524754 92.25092250922509
0.7551020408163265 0.7755102040816326

这里如果把风险偏好系数当作自变量，风险损失当成因变量，会发现当风险偏好系数小于0.76时，风险损失发生了暴涨，质变！

学习感慨/困惑

1.最近在管理运筹学和数学建模中都正好在讲线性规划，于是我先完成了第六章的习题。我已经学习了单纯形法，大M法，两阶段法，对偶单纯形法，灵敏度分析，还自学了动态规划。不得不说各种方法熟能生巧，但是我依稀记得刚入门线性规划时自己怎么算都算不对的情形！在python里，我不知道它求解器用的都是什么方法，所以它于我而言相当于是黑箱。另外也不得不感慨问题的千奇百怪，光线性规划就已经4种解法了，还有那么多规划，那么多解法，只能感叹前途漫漫了。
答：运筹学学习路线。

2.数学建模和人工智能的联系和区别？
答：我想可从二者定义出发加以解释。

3.当求解问题出现inf时，如何高效找到问题所在？目前已知的有看约束条件，看目标函数，看函数用法，看变量定义，但效率都很低就是，因为我不知道直接原因，而返回的结果也只有inf，不像报错还会指出是你哪行代码的问题，所以就会有很多要一一排查的东西！！！

4.多目标规划里当两个目标的大小不是一个量级时权重$w$是否要做调整？怎么调整？

5.在用cvxpy库处理大型线性规划问题（涉及上万行乘以上万行的矩阵时）时发生内存不足的现象！怎么办？
答：可考虑使用gurobi库，它相对cvxpy而言，cvxpy在运行时似乎会老老实实地生成大型的数组，而gurobi好像要聪明很多，运行时占用的内存小很多，很厉害。