数学建模算法与应用 线性规划（cvxpy包）
数学建模算法与应用 线性规划（使用cvxpy包）
说明
使用python中cvxpy库完成《数学建模算法与应用》中课后习题
因为本人也是初学者，若代码有错误还请各位指出
cvxpy库的使用
- 这里以简单例子来使用cvxpy库
maxz=4x1+3x2max\ z = 4x_1 + 3x_2 max z=4x1+3x2
s.t.={2x1+x2≤10x1+x2≤8x2≤7xi≥0s.t.= \begin{cases} 2x_1 + x_2 \leq 10\\ x_1 + x_2 \leq 8\\ x_2 \leq 7\\ x_i \geq 0 \end{cases} s.t.=⎩⎪⎪⎪⎨⎪⎪⎪⎧2x1+x2≤10x1+x2≤8x2≤7xi≥0
import cvxpy as cp
import numpy as np
# 决策变量
n = 2
# cvxpy使用cp.Variable(n,intger = True)中的intger参数来规定x变量为整数
x = cp.Variable(n, integer=True)
# 约束1
A1 = np.array([[2,1],[1,1],[0,1]])
b1 = np.array([10,8,7])
# 约束2
A2 = np.array([[1,0],[0,1]])
b2 = np.array([0,0])
# 目标函数
c = np.array([4, 3])
# 将最大化目标变为最小化
c = -c
# 定义问题，添加约束条件
prob = cp.Problem(cp.Minimize(c.T @ x),[A1 @ x <= b1, A2 @ x >= b2])
# 求解
ans = -round(prob.solve(solver = cp.ECOS_BB))
# 输出结果
print("目标函数最大值:", ans)
# 对x向量各元素取整数后再输出
value = []
for item in x.value:value.append(round(item))
print(value)
输出：
目标函数最大值: 26
[2, 6]
习题
习题1.1
maxz=3x1−x2−x3max \ z = 3x_1-x_2-x_3 max z=3x1−x2−x3
s.t.={x1−2x2+x3≤11−4x1+x2+2x3≥3−2x1+x3=1xi≥0s.t.= \begin{cases} x_1 - 2x_2 + x_3\leq 11\\ -4x_1 + x_2 + 2x_3 \geq 3\\ -2x_1 + x_3 = 1\\ x_i \geq 0 \end{cases} s.t.=⎩⎪⎪⎪⎨⎪⎪⎪⎧x1−2x2+x3≤11−4x1+x2+2x3≥3−2x1+x3=1xi≥0
import cvxpy as cp
import numpy as np
# 决策变量
n = 3
x = cp.Variable(n, integer=False)# 约束1
A1 = np.array([1, -2, 1])
b1 = np.array([11])# 约束2
A2 = np.array([-2, 0, 1])
b2 = np.array([1])# 约束3
A3 = np.array([[-4, 1, 2],[1, 0, 0],[0, 1, 0],[0, 0, 1]])
b3 = [3, 0, 0, 0]# 目标函数
c = np.array([3, -1, -1])
c = -c
# 定义问题，添加约束条件
prob = cp.Problem(cp.Minimize(c.T @ x),[A1 @ x <= b1, A2 @ x == b2, A3 @ x >= b3])# 求解
ans = -round(prob.solve(solver=cp.ECOS_BB))# 输出结果
print("目标函数最大值:", ans)
# 对x向量各元素取整数后再输出
print(x.value)
输出：
目标函数最大值: 2
[4. 1. 9.]
习题1.2
minz=∣x1∣+2∣x2∣+3∣x3∣+4∣x4∣min \ z = |x_1| + 2|x_2| + 3|x_3| + 4|x_4| min z=∣x1∣+2∣x2∣+3∣x3∣+4∣x4∣
s.t.={x1−x2−x3+x4=0x1−x2+x3−3x4=1x1−x2−2x3+3x4=−0.5s.t.= \begin{cases} x_1 - x_2 - x_3 + x_4 = 0\\ x_1 - x_2 + x_3 - 3x_4 = 1\\ x_1 - x_2 -2x_3 + 3x_4 = -0.5\\ \end{cases} s.t.=⎩⎪⎨⎪⎧x1−x2−x3+x4=0x1−x2+x3−3x4=1x1−x2−2x3+3x4=−0.5
- 做变量变换ui=xi+∣xi∣2u_{i} = \frac{x_{i}+|x_i|}{2}ui=2xi+∣xi∣,vi=∣xi∣−xi2v_{i} = \frac{|x_i|-x_i}{2}vi=2∣xi∣−xi，将模型变换为线性规划模型
mincT(u+v)min \ c^T(u+v) min cT(u+v)
s.t.={A(u−v)≤bu,v≥0s.t.= \begin{cases} A(u-v) \leq b\\ u,v \geq 0 \end{cases} s.t.={A(u−v)≤bu,v≥0
import cvxpy as cp
import numpy as np
# 决策变量
n = 4
u = cp.Variable(n, integer=False)
v = cp.Variable(n, integer=False)# 约束1
A1 = np.array([[1, -1, -1, 1],[1, -1, 1, -3],[1, -1, -2, 3]])
b1 = np.array([0, 1, -0.5])# 约束2
A2 = np.ones((4, 4))
for i in range(A2.shape[0]):for j in range(A2.shape[1]):if i == j:passelse:A2[i, j] = A2[i, j]*0b2 = np.array([0, 0, 0, 0])# 目标函数
c = np.arange(1, 5, 1)# 定义问题，添加约束条件
prob = cp.Problem(cp.Minimize(c.T @ (u+v)),[A1 @ (u-v) == b1, A2 @ u >= b2, A2 @ v >= b2])# 求解
ans = prob.solve(solver=cp.ECOS_BB)# 输出结果
print("目标函数最小值:", ans)
# 对x向量各元素取整数后再输出
print(u.value - v.value)print(A2)
输出：
目标函数最小值: 1.25
[ 2.50000000e-01 -2.97320702e-11 2.28691875e-11 -2.50000000e-01]
习题1.3
某厂生产三种产品Ⅰ、Ⅱ、Ⅲ。每种产品要经过AAA、BBB两道工序加工。设该厂有两种规格的设备能完成AAA工序，它们以A1,A2A_1,A_2A1,A2表示；有三种规格的设备能完成B工序，它们以B1B_1B1、B2B_2B2、B3B_3B3表示。产品Ⅰ可在AAA、BBB任何一种规格设备上加工。产品Ⅱ可在任何规格的AAA设备上加工，但完成BBB工序时，只能在B1B_{1}B1设备上加工；产品Ⅲ只能在A2A_2A2和B2B_2B2设备上加工。已知在各种机床设备的单件工时、原材料费、产品销售单价、各种设备有效台时以及满负荷操作时机床设备的费用如表所示，求最优的生产计划，使该厂利润最大。
算法设计
对产品Ⅰ来说，设以A1,A2A_1,A_2A1,A2完成AAA工序的产品分别为x1x_1x1、x2x_2x2件，转入BBB工序时，以B1B_1B1、B2B_2B2、B3B_3B3完成BBB工序的产品分别为x3x_3x3、x4x_4x4、x5x_5x5件；对产品Ⅱ来说，设以以A1,A2A_1,A_2A1,A2完成AAA工序的产品分别为x6x_6x6、x7x_7x7件，转入BBB工序时，以B1B_1B1完成BBB工序的产品为x8x_8x8件；对产品Ⅲ来说，设以A2A_2A2完成AAA工序的产品为x9x_9x9，则以B2B_2B2完成BBB工序的产品也为x9x_9x9件。
maxz=(1.25−0.25)(x1+x2)+(2−0.35)x8+(2.8−0.5)x9−300600(5x1+10x6)−32110000(7x2+9x7+12x9)−2504000(6x3+8x8)−7837000(4x4+11x9)−200400×7x5\begin{aligned} max \ z &= (1.25-0.25)(x_1 + x_2) + (2 - 0.35) x_8 + (2.8 - 0.5 ) x_9\\ &-\frac{300}{600}(5x_1 + 10x_6)-\frac{321}{10000}(7x_2 + 9x_7 + 12x_9)\\ &-\frac{250}{4000}(6x_3 + 8x_8)-\frac{783}{7000}(4x_4 + 11x_9)-\frac{200}{400} \times 7x_5 \end{aligned} max z=(1.25−0.25)(x1+x2)+(2−0.35)x8+(2.8−0.5)x9−600300(5x1+10x6)−10000321(7x2+9x7+12x9)−4000250(6x3+8x8)−7000783(4x4+11x9)−400200×7x5
s.t.={5x1+10x6≤60007x2+9x7+12x9≤100006x3+8x8≤40004x4+11x9≤70007x5≤4000x1+x2=x3+x4+x5x6+x7=x8xi≥0,i=1,2,...,9s.t.= \begin{cases} 5x_1 + 10x_6 \leq 6000\\ 7x_2 +9 x_7 + 12x_9 \leq 10000\\ 6x_3 + 8x_8 \leq 4000\\ 4x_4 + 11x_9 \leq 7000\\ 7x_5 \leq 4000\\ x_1 + x_2 = x_3 + x_4 + x_5\\ x_6 + x_7 = x_8\\ x_i \geq 0,i = 1,2,...,9 \end{cases} s.t.=⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧5x1+10x6≤60007x2+9x7+12x9≤100006x3+8x8≤40004x4+11x9≤70007x5≤4000x1+x2=x3+x4+x5x6+x7=x8xi≥0,i=1,2,...,9
# 决策变量
n = 9
x = cp.Variable(n, integer=True)# 约束1
A1 = np.array([[5,0,0,0,0,10,0,0,0],[0,7,0,0,0,0,9,0,12],[0,0,6,0,0,0,0,8,0],[0,0,0,4,0,0,0,0,11],[0,0,0,0,7,0,0,0,0]])
b1 = np.array([6000,10000,4000,7000,4000])# 约束2
A2 = np.ones((9, 9))
for i in range(A2.shape[0]):for j in range(A2.shape[1]):if i == j:passelse:A2[i, j] = A2[i, j]*0b2 = np.array([0,0,0,0,0,0,0,0,0])#约束3
A3=np.array([[1,1,-1,-1,-1,0,0,0,0],[0,0,0,0,0,1,1,-1,0]])
b3=np.array([0,0])# 定义问题，添加约束条件
prob = cp.Problem(cp.Maximize(x[0] + x[1] + 1.65 * x[7] + 2.3 * x[8] - 0.05 * (5 * x[0] + 10 * x[5])-0.0321 * (7 * x[1] + 9 * x[6] + 12 * x[8])-25/400 * (6 * x[2] + 8 * x[7])-783/7000 * (4 * x[3] + 11 * x[8])-0.35 * x[4]),[A1 @ x <= b1, A2 @ x >= b2, A3 @ x == b3])# 求解
ans = prob.solve(solver=cp.ECOS_BB)# 输出结果
print("目标函数最大值:", ans)
# 对x向量各元素取整数后再输出
print(x.value)
输出：
目标函数最大值: 1146.4142000003724
[ 1.20000000e+03 2.30000000e+02 1.73068974e-09 8.59000000e+025.71000000e+02 -1.28932697e-09 5.00000000e+02 5.00000000e+023.24000000e+02]
习题1.4
一架货机有三个货舱：前舱、中舱和后舱。三个货舱所能装载的货物的最大质量和体积有限制。为了飞机的平衡，三个货舱装载的货物质量必须与其最大的容许量成正比例。现有四类货物用该货机进行装运，货物的规格以及装运后获得的利润如表所示。
每种货物可以无限细分；
每种货物可以分布在一个或者多个货舱内；
不同的货物可以放在同一个货舱内，并且可以保证不留空隙。
如何装运，使货机飞行利润最大
算法设计
用i=1,2,3,4i = 1,2,3,4i=1,2,3,4分别表示货物1、货物2、货物3、货物4；j=1,2,3j = 1,2,3j=1,2,3分别表示前舱、中舱和后舱。设xij(i=1,2,3,4;j=1,2,3)x_{ij}(i = 1,2,3,4;j = 1,2,3)xij(i=1,2,3,4;j=1,2,3)表示第iii种货物装在第jjj个货舱内的质量，wj,vj(j=1,2,3)w_{j},v_{j}(j = 1,2,3)wj,vj(j=1,2,3)分别表示第jjj个舱的质量限制和体积限制，ai,bi,ci(i=1,2,3,4)a_{i},b_{i},c_{i}(i = 1,2,3,4)ai,bi,ci(i=1,2,3,4)分别表示可以运输的第iii种货物的质量，单位质量所占的空间和单位货物的利润。
maxz=∑i=14ci∑j=13xijmax \ z = \sum^{4}_{i = 1}c_{i}\sum^{3}_{j = 1}x_{ij} max z=i=1∑4cij=1∑3xij
s.t.={∑j=13xij≤ai,i=1,2,3,4∑i=14xij≤wi,j=1,2,3∑i=14bixij≤vi,i=1,2,3∑i=14xi110=∑i=14xi216=∑i=14xi38xij≥0,i=1,2,3,4;j=1,2,3s.t.= \begin{cases} \sum\limits^{3}_{j = 1} x_{ij}\leq a_{i},\ \ i = 1,2,3,4\\ \sum\limits^{4}_{i = 1} x_{ij}\leq w_{i},\ \ j = 1,2,3\\ \sum\limits^{4}_{i = 1} b_{i}x_{ij}\leq v_{i},\ \ i = 1,2,3\\ \frac{\sum\limits^{4}_{i = 1}x_{i1}}{10} = \frac{\sum\limits^{4}_{i = 1}x_{i2}}{16} = \frac{\sum\limits^{4}_{i = 1}x_{i3}}{8}\\ x_{ij} \geq 0,\ \ i = 1,2,3,4;j = 1,2,3 \end{cases} s.t.=⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧j=1∑3xij≤ai, i=1,2,3,4i=1∑4xij≤wi, j=1,2,3i=1∑4bixij≤vi, i=1,2,310i=1∑4xi1=16i=1∑4xi2=8i=1∑4xi3xij≥0, i=1,2,3,4;j=1,2,3
此次在变量设定中使用了
nonneg参数规定xijx_{ij}xij非负，可以不用写最后一条约束涉及到变量求和的地方使用
cvxpy包中内置sum函数，即cp.sum()更多内置函数参考cvxpy官方文档
# 决策变量
n = (4,3)
# nonneg参数，变量是否为非负
x = cp.Variable(n,nonneg = True)# 约束1
b1 = np.array([18,15,23,12])# 约束2
b2 = np.array([10,16,8])# 约束3
A3 = np.array([480,650,580,390])
b3 = np.array([6800,8700,5300])
# 定义问题，添加约束条件
c = np.array([3100,3800,3500,2850])
prob = cp.Problem(cp.Maximize(c @ cp.sum(x,axis = 1)),[cp.sum(x,axis = 1) <= b1,cp.sum(x,axis = 0) <= b2,A3 @ x <= b3,cp.sum(x[:,0])/10 == cp.sum(x[:,1])/16,cp.sum(x[:,1])/16 == cp.sum(x[:,2])/8])# 求解
ans = prob.solve(solver=cp.ECOS_BB)# 输出结果
print("目标函数最大值:", ans)
# 对x向量各元素取整数后再输出
print(np.sum(x.value,axis = 1))
输出：
目标函数最大值: 121515.7894845719
[3.19620414e-08 1.50000000e+01 1.59473684e+01 3.05263157e+00]
习题1.5
某部门在今后五年内考虑给下列项目投资，已知：
项目A，从第一年到第四年每年年初需要投资，并于次年末回收本利115%；
项目B，从第三年初需要投资，到第五年末能回收本利125%，但规定最大投资额不超过4万元；
项目C，第二年初需要投资，到第五年末能回收本利140%，但规定最大投资额不超过3万元；
项目D，五年内每年初可购买公债，于当年末归还，并加利息6%。
该部门现有资金10万元，问它应如何确定给这些项目每年的投资额，使到第五年年末拥有的资金的本利总额最大。
算法设计
用j=1,2,3,4j = 1,2,3,4j=1,2,3,4,分别表示项目AAA、BBB、CCC、DDD，用xij(i=1,2,3,4,5)x_{ij}(i = 1,2,3,4,5)xij(i=1,2,3,4,5)分别表示第i年年初给项目AAA、BBB、CCC、DDD的投资额。
maxz=1.15x41+1.40x23+1.25x32+1.06x54max \ z = 1.15x_{41} + 1.40x_{23}+ 1.25x_{32}+1.06x_{54} max z=1.15x41+1.40x23+1.25x32+1.06x54
s.t.={x11+x14=100000x21+x23+x24=1.06x14x31+x32+x34=1.15x11+1.06x24x41+x44=1.15x21+1.06x34x54=1.15x31+1.06x44x32≤40000x23≤30000xij≥0,i=1,2,3,4,5;j=1,2,3,4s.t.= \begin{cases} x_{11} + x_{14} = 100000\\ x_{21} + x_{23} + x_{24} = 1.06x_{14}\\ x_{31} + x_{32} + x_{34} = 1.15x_{11} + 1.06x_{24}\\ x_{41} + x{44} = 1.15x_{21} + 1.06x_{34}\\ x_{54} = 1.15x_{31} + 1.06x_{44}\\ x_{32} \leq 40000\\ x_{23} \leq 30000\\ x_{ij} \geq 0,\ \ i = 1,2,3,4,5;j = 1,2,3,4 \end{cases} s.t.=⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧x11+x14=100000x21+x23+x24=1.06x14x31+x32+x34=1.15x11+1.06x24x41+x44=1.15x21+1.06x34x54=1.15x31+1.06x44x32≤40000x23≤30000xij≥0, i=1,2,3,4,5;j=1,2,3,4
- 换了一个求解器，GLOP
# 决策变量
n = (5,4)
# nonneg参数，变量是否为非负
x = cp.Variable(n,nonneg = True)prob = cp.Problem(cp.Maximize(1.15 * x[3,0] + 1.40 * x[1,2] + 1.25 * x[2,1] + 1.06 * x[4,3]),[x[0,0] + x[0,3] == 100000,x[1,0] + x[1,2] + x[1,3] == 1.06 * x[0,3],x[2,0] + x[2,1] + x[2,3] == 1.15 * x[0,0] + 1.06 * x[1,3],x[3,0] + x[3,3] == 1.15*x[1,0] + 1.06*x[2,3],x[4,3] == 1.15*x[2,0] + 1.06 * x[3,3],x[2,1] <= 40000,x[1,2] <= 30000])# 求解
ans = prob.solve(solver=cp.GLOP)# 输出结果
print("目标函数最大值:", ans)
# 对x向量各元素取整数后再输出
print(x.value)
输出：
目标函数最大值: 143750.0
[[34782.60869565 -0. -0. 65217.39130435][39130.43478261 -0. 30000. -0. ][ -0. 40000. -0. -0. ][45000. -0. -0. -0. ][ -0. -0. -0. -0. ]]
习题1.6
食品厂用三种原料生产两种糖果，糖果的成分要求和销售价如表所示。该厂根据订单至少需要生产600kg高级奶糖、800kg水果糖，为求最大利润，试建立线性规划模型求解。
算法设计
用i=1,2i = 1,2i=1,2分别表示高级奶糖和水果糖，用j=1,2,3j = 1,2,3j=1,2,3分别表示原料AAA、BBB、CCC。设xij(i=1,2;j=1,2,3)x_{ij}(i = 1,2;j = 1,2,3)xij(i=1,2;j=1,2,3)表示生产第iii种糖用的第jjj种原料的量，aia_{i}ai表示第iii种糖果的需求量，bjb_{j}bj表示第jjj种原料的可供量。
maxz=4x11+12x12+16x13−5x21+3x22+7x23max \ z = 4x_{11} + 12x_{12} + 16x_{13}-5x_{21} + 3x_{22} + 7x_{23} max z=4x11+12x12+16x13−5x21+3x22+7x23
s.t.={∑j=13xij≥ai,i=1,2∑i=12xij≤bj,j=1,2,3x11≥50%(x11+x12+x13)x12≥25%(x11+x12+x13)x13≤10%(x11+x12+x13)x21≤40%(x21+x22+x23)x22≤40%(x21+x22+x23)x23≥15(x21+x22+x23)xij≥0,i=1,2;j=1,2,3s.t.= \begin{cases} \sum\limits^{3}_{j = 1}x_{ij} \geq a_{i}, \ i = 1,2\\ \sum\limits^{2}_{i = 1}x_{ij} \leq b_{j}, \ j = 1,2,3\\ x_{11} \geq 50\%(x_{11} + x_{12} + x_{13})\\ x_{12} \geq 25\%(x_{11} + x_{12} + x_{13})\\ x_{13} \leq 10\%(x_{11} + x_{12} + x_{13})\\ x_{21} \leq 40\%(x_{21} + x_{22} + x_{23})\\ x_{22} \leq 40\%(x_{21} + x_{22} + x_{23})\\ x_{23} \geq 15(x_{21} + x_{22} + x_{23})\\ x_{ij} \geq 0,\ \ i = 1,2;j = 1,2,3 \end{cases} s.t.=⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧j=1∑3xij≥ai, i=1,2i=1∑2xij≤bj, j=1,2,3x11≥50%(x11+x12+x13)x12≥25%(x11+x12+x13)x13≤10%(x11+x12+x13)x21≤40%(x21+x22+x23)x22≤40%(x21+x22+x23)x23≥15(x21+x22+x23)xij≥0, i=1,2;j=1,2,3
# 决策变量
n = (2,3)
# nonneg参数，变量是否为非负
x = cp.Variable(n,nonneg = True)
# 约束1
b1 = [600,800]# 约束2
b2 = [600,750,625]
# cp.sum()中axis = 1按行求和
# cp.sum()中axis = 0按列求和
prob = cp.Problem(cp.Maximize(np.array([4,12,16]) @ x[0,:] + np.array([-5,3,7]) @ x[1,:]),[cp.sum(x,axis = 1) >= b1,cp.sum(x,axis = 0) <= b2,x[0,0] >= 0.5 * cp.sum(x[0,:],axis = 0),x[0,1] >= 0.25 * cp.sum(x[0,:],axis = 0),x[0,2] <= 0.1 * cp.sum(x[0,:],axis = 0),x[1,0] <= 0.4 * cp.sum(x[1,:],axis = 0),x[1,1] <= 0.4 * cp.sum(x[1,:],axis = 0),x[1,2] >= 0.15 * cp.sum(x[1,:],axis = 0)])# 求解
ans = prob.solve(solver=cp.GLOP)# 输出结果
print("目标函数最大值:", ans)
# 对x向量各元素取整数后再输出
print(x.value)
输出：
目标函数最大值: 14200.0
[[600. 575. -0.][ -0. 175. 625.]]
习题1.7
求解下列线性规划问题，其中矩阵A=(aij)100×150A = (a_{ij})_{100 \times 150}A=(aij)100×150中的元素aija_{ij}aij为[0,10][0,10][0,10]上的随机整数。
maxvmax \ v max v
s.t.={∑i=1100aijxi≥v,j=1,2,...,150∑i=1100xi=100xi≥0,i=1,2,...,100s.t.= \begin{cases} \sum\limits^{100}_{i = 1}a_{ij}x_{i} \geq v, \ j = 1,2,...,150\\ \sum\limits^{100}_{i = 1}x_{i} = 100\\ x_{i} \geq 0,\ i = 1,2,...,100 \end{cases} s.t.=⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧i=1∑100aijxi≥v, j=1,2,...,150i=1∑100xi=100xi≥0, i=1,2,...,100
#生成随机数种子，使结果固定
np.random.seed(42)
a = np.random.randint(1, 11, size=(100,150))
# 决策变量
n = np.array([100,1])
# nonneg参数，变量是否为非负
x = cp.Variable(n[0],nonneg = True)
v = cp.Variable(n[1])# cp.sum()中axis = 1按行求和
# cp.sum()中axis = 0按列求和
prob = cp.Problem(cp.Maximize(v),[a.T @ x >= v,cp.sum(x) == 100])# 求解
ans = prob.solve(solver=cp.GLOP)# 输出结果
print("目标函数最大值:", ans)
# 对x向量各元素取整数后再输出
print(x.value)
输出：
目标函数最大值: 538.2560346542277
[-0. -0. 0.32868907 1.83981903 1.60824063 1.429344122.35357916 -0. 3.38452136 -0. 0.00595007 1.6272227-0. 0.56945258 -0. 2.54100266 1.43375965 0.863279972.11702508 2.40072583 -0. 2.52730273 -0. -0.1.13784709 -0. 0.62438124 -0. 0.11470483 2.13261487-0. 2.78707508 2.7390843 1.22599971 1.36865498 4.551308193.26987533 -0. -0. -0. 4.26209229 1.43778780.20195965 0.63892766 0.13196033 -0. -0. 1.318078780.75303873 3.36903244 -0. 0.81368899 0.2527292 1.944166551.58858316 0.21577316 1.7780698 -0. -0. 1.771560090.22421819 -0. 0.94448177 -0. -0. -0.1.13899834 3.02491754 2.22561055 3.26464305 -0. 1.389167722.06832119 -0. 2.6891509 -0. 0.42024825 1.69604933-0. 0.03469465 -0. 2.85642199 2.70167707 -0.1.01923824 1.87565017 0.98045756 -0. 0.50588775 1.73655473-0. -0. -0. 0.06062595 -0. 1.469786141.08026003 1.13403002 -0. -0. ]
习题1.8
设xn+1=max1≤i≤n{qixi}x_{n+1} = \mathop{max}\limits_{1 \leq i \leq n}\{q_{i}x_{i}\}xn+1=1≤i≤nmax{qixi}，则模型二可以线性化为
minxn+1min \ x_{n+1} min xn+1
s.t.={qixi≤xn+1,i=1,2,...,n∑i=0n(ri−pi)xi≥kM∑i=0n(1+pi)xi=Mxi≥0,i=0,1,...,ns.t.= \begin{cases} q_{i}x_{i} \leq x_{n+1},\ i = 1,2,...,n\\ \sum\limits^{n}_{i = 0}(r_i - p_i)x_i \geq k_{M}\\ \sum\limits^{n}_{i = 0}(1 + p_i)x_i = M\\ x_{i} \geq 0,\ i = 0,1,...,n \end{cases} s.t.=⎩⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎧qixi≤xn+1, i=1,2,...,ni=0∑n(ri−pi)xi≥kMi=0∑n(1+pi)xi=Mxi≥0, i=0,1,...,n
# 决策变量
n = 6
# nonneg参数，变量是否为非负
x = cp.Variable(n,nonneg = True)# 约束1
A1 = np.array([[0,0.025,0,0,0,-1],[0,0,0.015,0,0,-1],[0,0,0,0.055,0,-1],[0,0,0,0,0.026,-1],[-0.05,-0.27,-0.19,-0.185,-0.185,0]])# 约束2
A2 = np.array([1,1.01,1.02,1.045,1.065,0])
b2 = np.array([1])# 每轮结果
ANS = []
for k in np.arange(0.05,0.255,0.005):b1 = np.array([0,0,0,0,-k])prob = cp.Problem(cp.Minimize(x[5]),[A1 @ x <= b1,A2 @ x == b2])# 求解ans = prob.solve(solver=cp.GLOP)ANS.append([ans])if k >= 0.205 and k<= 0.210:print(x.value)
输出：
[-0. 0.30887428 0.51479046 0.1403974 0.01524453 0.00772186]
绘制图形观察
import matplotlib as mpl
import matplotlib.pyplot as plt#中文显示问题
plt.rcParams['font.sans-serif']=['SimHei']
plt.rcParams['axes.unicode_minus'] = False# notebook嵌入图片
%matplotlib inline
# 提高分辨率
%config InlineBackend.figure_format='retina'# 忽略警告
import warnings
warnings.filterwarnings('ignore')plt.figure(dpi = 600,figsize = (5,5))
plt.plot(np.arange(0.05,0.255,0.005),ANS)
plt.xlabel("Profit")
plt.ylabel("Risk")
plt.title("Profit and Risk chart")
plt.savefig("Profit and Risk chart")
结语
最后这里提供给大家源代码文件
数学建模算法与应用 线性规划（cvxpy包）相关推荐
- 个人数学建模算法库之线性规划模型
线性规划模型 模型描述 Matlab代码 基于求解器 基于问题 Python代码 Scipy库 Plup库 模型描述 m i n z = f T x ( 1 ) s . t . { A x ≤ b ( ...
- 司守奎《数学建模算法与应用》课后习题：线性规划
系列目录 司守奎<数学建模算法与应用>课后习题:线性规划 模拟退火算法解决旅行商问题详解 遗传算法解决旅行商问题 MATLAB实现 基于模拟退火优化的投影寻踪评价法求解供货能力评价模型 目 ...
- 数学建模算法与应用学习day1——线性规划问题整数规划问题
以下内容来自司守奎编写的数学建模算法与应用学习,主要是记录自己的学习历程,转载还请标明出处! 一.线性规划 知识点 1.1线性规划问题 1.1.2线性规划解的概念 f = [-2 ; -3 ; 5]; ...
- 数学建模算法学习笔记 已完结
这是为了准备国赛突击学习的模型算法,我在原有的基础上加上自己的理解虽然不知道对不对,就是为了记录下自己学的模型他究竟是个什么东西,语言通俗,但是极不准确,只适合做一个大概的了解,建议大家详细的还是要看 ...
- 数学建模算法：支持向量机_从零开始的算法：支持向量机
数学建模算法:支持向量机 从零开始的算法 (Algorithms From Scratch) A popular algorithm that is capable of performing lin ...
- matlab中x从0到5不含0,关于MATLAB的数学建模算法学习笔记
关于MATLAB的数学建模算法学习笔记 目录 线性规划中应用: (3) 非线性规划: (3) 指派问题;投资问题:(0-1问题) (3) 1)应用fmincon命令语句 (3) 2)应用指令函数:bi ...
- python数学建模（二）线性规划2.实战(思路清晰\过程完整、详细)
文章目录 (一)简单陈述本文章的内容 (二)线性规划例题(实战) 2.1 实战题目 2.2 符号规定和基本假设 2.3 模型的分析 2.4 模型的建立 2.5 模型一的求解和分析 2.5.1 (代码) ...
- LL1分析构造法_数学建模算法--最优赋权法（含代码）
数学建模算法--最优赋权法(含代码) 作者:郑铿城 本次介绍数学建模和科研写作的方法--最优赋权法最优赋权法经常用于分析评价类问题,从该算法的名称就可以看到,该算法首先要体现"最优" ...
- python dendrogram_【聚类分析】《数学建模算法与应用》第十章 多元分析 第一节 聚类分析 python实现...
第十章 多元分析 第一节 聚类分析 介绍 这里是司守奎教授的<数学建模算法与应用>全书案例代码python实现,欢迎加入此项目将其案例代码用python实现 GitHub项目地址:Math ...
最新文章
- Jenkins 2.16.3默认没有Launch agent via Java Web Start，如何配置使用
- python画图代码turtle-使用Python的turtle模块画图的方法
- DDD - 聚合与聚合根_如何理解 Respository与DAO
- java 阻塞 socket_java socket非阻塞I/O
- Elipse 、Idea配置 Java-Code-Formatter
- Linux学习总结（55）——Linux 运维常用脚本
- Linux常用最基础命令总结
- 基于javaweb的员工绩效考核系统
- YAML语法详细总结
- 基于BIND实现智能DNS解析
- 腾讯企业邮箱申请步骤
- 【微信JSSDK】PHP版微信录音文件下载
- Microsoft 解决方案框架版本 3.0 概述（MSF3.0）
- CSS中使盒子移动方法总结
- Java微信公众平台开发(一)--接入微信公众平台
- 普元 AppServer 6.5 将springboot应用部署到应用服务器，上传文件时报错：Caused by: org.springframework.web.multipart.Multipar
- 梦兴阁给我生活带来的改变
- 形态等位点对迭代次数的贡献
- 矩阵的初等变换的应用
- Docker compose - 最开始的version 字段是什么，为什么要写这个字段