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

5.1 求下列线性规划的解

$$\begin{gathered} maxz=8x_1-2x_2+3x_3-x_4-2x_5, \\ s.t.\{\begin{array}{l}{x}_1+x_2+x_3+x_4+x_5⩽400,\\ x_1+2x_2+2x_3+x_4+6x_5⩽800,\\ 2x_1+x_2+6x_3⩽200,\\ x_3+x_4+5x_5⩽200,0⩽x_i⩽99,i=1,2,3,4;\\ x_5⩾-10.\end{array} \end{gathered}$$

#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$

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$为 max ⁡ 1 ⩽ i ⩽ n { q i x i } \max_{1\leqslant i\leqslant n}\{q_ix_i\} max1⩽i⩽n​{qi​xi​}

#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\%$$	$$⩾30\%$$		2	2000
B				1.5	2500
C	$$⩽20\%$$	$$⩽50\%$$	$$⩽60\%$$	1	1200
加工费/(元/kg)	0.5	0.4	0.3
售价/(元/kg)	3.4	2.85	2.25

设

	甲	乙	丙
A	$$x_1$$	$$x_2$$	$$x_3$$
B	$$x_4$$	$$x_5$$	$$x_6$$
C	$$x_7$$	$$x_8$$	$$x_9$$

建立模型：
$$\begin{gathered} minz=-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.\{\begin{array}{l}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\end{array} \end{gathered}$$

#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 试求多目标线性规划问题

$$\begin{gathered} max{z}_1=100x_1+90x_2+80x_3+70x_4 \\ min{z}_2=3x_2+2x_4 \\ s.t.\{\begin{array}{l}{x}_1+x_2⩾30,\\ x_3+x_4⩾30,\\ 3x_1+2x_3⩽120,\\ 3x_2+2x_4⩽48,\\ x_i⩾0,i=1,2,3,4.\end{array} \end{gathered}$$

#把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)