Python数学实验与建模课后习题第5章解析

5.1 求下列线性规划的解

m a x z = 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 ⩽ 400 , 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 , i = 1 , 2 , 3 , 4 ; x 5 ⩾ − 10. max z=8x_1-2x_2+3x_3-x_4-2x_5, \\s.t. \left\{ \begin{array}{l} x_1+x_2+x_3+x_4+x_5\leqslant400, \\x_1+2x_2+2x_3+x_4+6x_5\leqslant800, \\2x_1+x_2+6x_3\leqslant200, \\x_3+x_4+5x_5\leqslant200, 0\leqslant x_i\leqslant99,i=1,2,3,4; \\x_5\geqslant-10. \end{array} \right. maxz=8x1−2x2+3x3−x4−2x5,s.t.⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧x1+x2+x3+x4+x5⩽400,x1+2x2+2x3+x4+6x5⩽800,2x1+x2+6x3⩽200,x3+x4+5x5⩽200,0⩽xi⩽99,i=1,2,3,4;x5⩾−10.

#scipy.optimize
from scipy.optimize import linprog
c=[-8,2,-3,1,2]
A=[[1,1,1,1,1],[1,2,2,1,6],[2,1,6,0,0],[0,0,1,1,5]]
b=[400,800,200,200]
bounds=((0,99),(0,99),(0,99),(0,99),(-10,None))
res=linprog(c,A,b,None,None,bounds)
print("最优解:",-res.fun)
print("最优解为:",res.x)

#cvxopt.solvers
import numpy as np
from cvxopt import matrix,solvers
c=matrix([-8.,2,-3,1,2])
A=matrix([[1.,1,1,1,1],[1,2,2,1,6],[2,1,6,0,0],[0,0,1,1,5],
[-1,0,0,0,0],[0,-1,0,0,0],[0,0,-1,0,0],[0,0,0,-1,0],[0,0,0,0,-1],
[1,0,0,0,0],[0,1,0,0,0],[0,0,1,0,0],[0,0,0,1,0]]).T
b=matrix([400.,800,200,200,0,0,0,0,10,99,99,99,99])
sol=solvers.lp(c,A,b)
print("最优解为:\n",sol['x'])
print("最优值为:",-sol['primal objective'])
#如有要求x非负之类的条件，不建议使用cvxopt.solvers，A矩阵构造复杂

#cvxpy
import cvxpy as cp
from matplotlib.pyplot import *
rc('text',usetex='True');rc('font',size=16)
x=cp.Variable((5,1))
y=cp.Variable(1)
obj=cp.Minimize([-8,2,-3,1,2]@x)
a1=np.array([[1,1,1,1,1],[1,2,2,1,6],[2,1,6,0,0],[0,0,1,1,5]])
a2=np.array([[400],[800],[200],[200]])
con=[a1@x<=a2,x[:-1,0]>=0,x[:-1,0]<=99,x[-1,0]>=-10]
prob=cp.Problem(obj,con)
prob.solve()
print("最优值为:",-prob.value)
print("最优解为:\n",x.value)

得解： x 1 = 99 , x 2 = 0 , x 3 = 0.3333 , x 4 = 0 , x 5 = − 10 , z 的最大值为 813 x_1=99,x_2=0,x_3=0.3333,x_4=0,x_5=-10,z的最大值为813 x1=99,x2=0,x3=0.3333,x4=0,x5=−10,z的最大值为813

5.2 求5.4节模型三的解。

模型如下：
m i n { s x 5 − ( 1 − s ) ( 0.05 x 0 + 0.27 x 1 + 0.19 x 2 + 0.185 x 3 + 0.185 x 4 ) } , s . t . { 0.025 x 1 − x 5 ⩽ 0 0.015 x 2 − x 5 ⩽ 0 0.055 x 3 − x 5 ⩽ 0 0.026 x 4 − x 5 ⩽ 0 x 0 + 1.01 x 1 + 1.02 x 2 + 1.045 x 3 + 1.065 x 4 = 1 x i ⩾ 0 , i = 0 , 1 , 2 , 3 , 4 , 5 min \{sx_5-(1-s)(0.05x_0+0.27x_1+0.19x_2+0.185x_3+0.185x_4)\},\\ s.t.\left\{ \begin{array}{l} 0.025x_1-x_5\leqslant0\\ 0.015x_2-x_5\leqslant0\\ 0.055x_3-x_5\leqslant0\\ 0.026x_4-x_5\leqslant0\\ x_0+1.01x_1+1.02x_2+1.045x_3+1.065x_4=1\\ x_i\geqslant0,i=0,1,2,3,4,5 \end{array} \right. min{sx5−(1−s)(0.05x0+0.27x1+0.19x2+0.185x3+0.185x4)},s.t.⎩⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎧0.025x1−x5⩽00.015x2−x5⩽00.055x3−x5⩽00.026x4−x5⩽0x0+1.01x1+1.02x2+1.045x3+1.065x4=1xi⩾0,i=0,1,2,3,4,5
其中 x 5 x_5 x5为 max ⁡ 1 ⩽ i ⩽ n { q i x i } \max_{1\leqslant i\leqslant n}\{q_ix_i\} max1⩽i⩽n{qixi}

#linprog
from scipy.optimize import *
from matplotlib.pyplot import *
rc('text',usetex=True);rc('font',size=16)
A=[[0,0.025,0,0,0,-1],[0,0,0.15,0,0,-1],[0,0,0,0.55,0,-1],[0,0,0,0,0.026,-1]]
b=[0,0,0,0]
Aeq=[[1,1.01,1.02,1.045,1.065,0]]
beq=[1]
bound=((0,None),(0,None),(0,None),(0,None),(0,None),(0,None))
s=0;ss=[];aa=[]
while s<=1:c=[-(1-s)*0.05,-(1-s)*0.27,-(1-s)*0.19,-(1-s)*0.185,-(1-s)*0.185,s]res=linprog(c,A,b,Aeq,beq)ss.append(s);aa.append(-res.fun)s=s+0.01
plot(ss,aa,'r.')
show()

#cvxpy
import cvxpy as cp
import numpy as np
from matplotlib.pyplot import *
rc('text',usetex=True);rc('font',size=16)
x=cp.Variable(6)
a1=np.array([0.025,0.015,0.055,0.026])
a2=np.array([0.05,0.27,0.19,0.185,0.185])
a3=np.array([1,1.01,1.02,1.045,1.065])
con=[cp.multiply(a1,x[1:5])-x[5]<=0,a3@x[:-1]==1,x>=0]
s=0;ss=[];aa=[]
while s<=1:obj=cp.Minimize(s*x[-1]-(1-s)*a2@x[:-1])prob=cp.Problem(obj,con)prob.solve(solver='GLPK_MI')ss.append(s);aa.append(-prob.value)s=s+0.01
plot(ss,aa,'r.')
show()

5.3 某股民决定对6家公司的股票进行投资，根据对这6家公司的了解，估计了这6家公司股票的明年预期收益和这6种股票收益的协方差矩阵的数据见表5.6要获得至少25%的预期收益，最小风险是多少？

表5.6 公司股票明年预期收益和收益的协方差矩阵数据

股票	收益率/%	公司1	公司2	协方差公司3	公司4	公司5	公司6
公司1	20	0.032	0.005	0.03	-0.031	-0.027	0.01
公司2	42	0.005	0.1	0.085	-0.07	-0.05	0.02
公司3	100	0.03	0.085	0.333	-0.11	-0.02	0.042
公司4	50	-0.031	-0.07	-0.11	0.125	0.05	-0.06
公司5	46	-0.027	-0.05	-0.02	0.05	0.065	-0.02
公司6	30	0.01	0.02	0.042	-0.06	-0.02	0.08

建立模型：
m i n { x 6 } s . t . { 0.032 x 0 − x 6 ⩽ 0 , 0.1 x 1 − x 6 ⩽ 0 , 0.333 x 2 − x 6 ⩽ 0 , 0.125 x 3 − x 6 ⩽ 0 , 0.065 x 4 − x 6 ⩽ 0 , 0.08 x 5 − x 6 ⩽ 0 , x 0 + x 1 + x 2 + x 3 + x 4 + x 5 = 1 , − 20 x 0 − 42 x 1 − 100 x 2 − 50 x 3 − 46 x 4 − 30 x 5 ⩽ − 25 , 0 ⩽ x i ⩽ 1 , i = 0 , 1 , ⋯ , 5 , x 6 ⩾ 0. min\{x_6\}\\ s.t.\left\{ \begin{array}{l} 0.032x_0-x_6\leqslant 0,\\ 0.1x_1-x_6\leqslant 0,\\ 0.333x_2-x_6\leqslant 0,\\ 0.125x_3-x_6\leqslant 0,\\ 0.065x_4-x_6\leqslant 0,\\ 0.08x_5-x_6\leqslant 0,\\ x_0+x_1+x_2+x_3+x_4+x_5=1,\\ -20x_0-42x_1-100x_2-50x_3-46x_4-30x_5\leqslant -25,\\ 0\leqslant x_i\leqslant 1,i=0,1,\cdots ,5,x_6\geqslant 0. \end{array} \right. min{x6}s.t.⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧0.032x0−x6⩽0,0.1x1−x6⩽0,0.333x2−x6⩽0,0.125x3−x6⩽0,0.065x4−x6⩽0,0.08x5−x6⩽0,x0+x1+x2+x3+x4+x5=1,−20x0−42x1−100x2−50x3−46x4−30x5⩽−25,0⩽xi⩽1,i=0,1,⋯,5,x6⩾0.

#cvxpy
import cvxpy as cp
import numpy as np
x=cp.Variable(7)
a1=np.array([1,1,1,1,1,1])
a2=np.array([20,42,100,50,46,30])
a3=np.array([0.032,0.1,0.333,0.125,0.065,0.08])
obj=cp.Minimize(x[6])
con=[cp.multiply(x[:6],a3)-x[6]<=0,a1@x[:-1]==1,a2@x[:-1]>=25,x[:-1]>=0,x[:-1]<=1,x[-1]>=0]
prob=cp.Problem(obj,con)
prob.solve(solver='GLPK_MI')
print('最小风险：',prob.value)

#linprog
from scipy.optimize import *
c=[0,0,0,0,0,0,1]
A=[[-20,-42,-100,-50,-46,-30,0],[0.032,0,0,0,0,0,-1],[0,0.1,0,0,0,0,-1],
[0,0,0.333,0,0,0,-1],[0,0,0,0.125,0,0,-1],
[0,0,0,0,0.065,0,-1],[0,0,0,0,0,0.08,-1]]
b=[25,0,0,0,0,0,0]
Aeq=[[1,1,1,1,1,1,0]]
beq=[1]
bound=((0,1),(0,1),(0,1),(0,1),(0,1),(0,1),(0,None))
res=linprog(c,A,b,Aeq,beq,bound)
print("最小风险:",res.fun)

5.4 某糖果厂用原料A,B,C加工成三种不同牌号的糖果甲,乙,丙。已知各种牌号糖果中A,B,C含量、原料成本、各种原料的每月限制用量，三种牌号糖果的单位加工费及售价，如表5.7所示。问该厂生产这三种牌号糖果各多少千克，才能使其获利最大。试建立这个问题的线性规划模型。

表5.7 成产规格等数据

原料	甲	乙	丙	原料成本/(元/kg)	每月限制用量/kg
A	⩾ 60 % \geqslant60\% ⩾60%	⩾ 30 % \geqslant30\% ⩾30%		2	2000
B				1.5	2500
C	⩽ 20 % \leqslant20\% ⩽20%	⩽ 50 % \leqslant50\% ⩽50%	⩽ 60 % \leqslant60\% ⩽60%	1	1200
加工费/(元/kg)	0.5	0.4	0.3
售价/(元/kg)	3.4	2.85	2.25

设

	甲	乙	丙
A	x 1 x_1 x1	x 2 x_2 x2	x 3 x_3 x3
B	x 4 x_4 x4	x 5 x_5 x5	x 6 x_6 x6
C	x 7 x_7 x7	x 8 x_8 x8	x 9 x_9 x9

建立模型：
m i n z = − 0.9 x 1 − 0.45 x 2 + 0.05 x 3 − 1.4 x 4 − 0.95 x 5 − 0.45 x 6 − 1.9 x 7 − 1.45 x 8 − 0.95 x 9 s . t . { 0.6 ( x 1 + x 4 + x 7 ) ⩽ x 1 0.3 ( x 2 + x 5 + x 8 ) ⩽ x 2 x 7 ⩽ 0.2 ( x 1 + x 4 + x 7 ) x 8 ⩽ 0.5 ( x 2 + x 5 + x 8 ) x 9 ⩽ 0.6 ( x 3 + x 6 + x 9 ) x 1 + x 2 + x 3 ⩽ 2000 x 4 + x 5 + x 6 ⩽ 2500 x 7 + x 8 + x 9 ⩽ 1200 min z=-0.9x_1-0.45x_2+0.05x_3-1.4x_4-0.95x_5-0.45x_6-1.9x_7-1.45x_8-0.95x_9 \\s.t.\left\{ \begin{array}{l} 0.6(x_1+x_4+x_7)\leqslant x_1\\ 0.3(x_2+x_5+x_8)\leqslant x_2\\ x_7\leqslant 0.2(x_1+x_4+x_7)\\ x_8\leqslant 0.5(x_2+x_5+x_8)\\ x_9\leqslant 0.6(x_3+x_6+x_9)\\ x_1+x_2+x_3\leqslant 2000\\ x_4+x_5+x_6\leqslant 2500\\ x_7+x_8+x_9\leqslant 1200 \end{array} \right. minz=−0.9x1−0.45x2+0.05x3−1.4x4−0.95x5−0.45x6−1.9x7−1.45x8−0.95x9s.t.⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧0.6(x1+x4+x7)⩽x10.3(x2+x5+x8)⩽x2x7⩽0.2(x1+x4+x7)x8⩽0.5(x2+x5+x8)x9⩽0.6(x3+x6+x9)x1+x2+x3⩽2000x4+x5+x6⩽2500x7+x8+x9⩽1200

#linprog
from scipy.optimize import *
import numpy as np
c=[-0.9,-0.45,0.05,-1.4,-0.95,-0.45,-1.9,-1.45,-0.95]
A=[[-0.4,0,0,0.6,0,0,0.6,0,0],[0,-0.7,0,0,0.3,0,0,0.3,0],[-0.2,0,0,-0.2,0,0,0.8,0,0],[0,-0.5,0,0,-0.5,0,0,0.5,0],[0,0,-0.6,0,0,-0.6,0,0,0.4],[1,1,1,0,0,0,0,0,0],[0,0,0,1,1,1,0,0,0],[0,0,0,0,0,0,1,1,1]]
b=[0,0,0,0,0,2000,2500,1200]
bound=((0,2000),(0,2000),(0,2000),(0,2500),(0,2500),(0,2500),(0,1200),(0,1200),(0,1200))
res=linprog(c,A,b,None,None,bound)
a=np.array(np.reshape(res.x,(3,3)))
print("应生产糖果\n甲:%.0f千克"%(a[0,0]+a[1,0]+a[2,0]))
print("乙:%.0f千克"%(a[0,1]+a[1,1]+a[2,1]))
print("丙:%.0f千克"%(a[0,2]+a[1,2]+a[2,2]))

5.5 试求多目标线性规划问题

m a x z 1 = 100 x 1 + 90 x 2 + 80 x 3 + 70 x 4 m i n z 2 = 3 x 2 + 2 x 4 s . t . { x 1 + x 2 ⩾ 30 , x 3 + x 4 ⩾ 30 , 3 x 1 + 2 x 3 ⩽ 120 , 3 x 2 + 2 x 4 ⩽ 48 , x i ⩾ 0 , i = 1 , 2 , 3 , 4. maxz_1=100x_1+90x_2+80x_3+70x_4 \\min z_2=3x_2+2x_4 \\s.t. \left\{ \begin{array}{l} x_1+x_2\geqslant30, \\x_3+x_4\geqslant30, \\3x_1+2x_3\leqslant120, \\3x_2+2x_4\leqslant48, \\x_i\geqslant0,i=1,2,3,4. \end{array} \right. maxz1=100x1+90x2+80x3+70x4minz2=3x2+2x4s.t.⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧x1+x2⩾30,x3+x4⩾30,3x1+2x3⩽120,3x2+2x4⩽48,xi⩾0,i=1,2,3,4.

#把max放进约束条件,不妨假定小于等于10000
#linprog
from scipy.optimize import linprog
c=[0,0,3,2]
A=[[-1,-1,0,0],[0,0,-1,-1],[3,0,2,0],[0,3,0,2],[100,90,80,70]]
b=[[-30],[-30],[120],[48],[10000]]
LB=[0]*len(c)
UB=[None]*len(c)
bound=tuple(zip(LB,UB))
res=linprog(c,A,b,None,None,bound)
print("目标函数的最小解:",res.fun)
print("最优解为:",res.x)