(一)简单陈述本文章的内容
python建模会持续更新,用途是只作为个人笔记。我博客中的所有资料都可通过我提供的链接永久获取,希望大家一起相互促进,相互努力。
本章博客介绍用线性规划进行真题实战。思路清晰、易懂,过程详细,如有疑问或者错误可以和我沟通哈。
本文所有代码、资料获取链接:https://pan.baidu.com/s/10qPE4pNtGKpCppCNTZ6Inw
提取码:nvyu
(二)线性规划例题(实战)
上一篇我们讲了python数学建模(一)线性规划,现在我们来一道真题实战。让我们学会怎么合理运用之前所学习过的求解线性规划问题的库和模块。
2.1 实战题目
市场上有n种资产 s i s_i si(i = 1,2,···,n)可以选择,现用数额为 M 的子昂当大的资金作为一个时期的投资。这 n 种资产在这一时期内购买 s i s_i si 的平均收益为 r i r_i ri ,风险损失率为 q i q_i qi ,投资越分散,总的风险越少,总体风险可用投资的 s i s_i si 中最大的一个风险来度量。
购买 s i s_i si 时要交付交易费,费率为 p i p_i pi ,当购买额不超过给定值 u i u_i ui 时,交易费按照购买 u i u_i ui 计算。另外,假定同期银行存款利率时 r 0 r_0 r0,既无交易费也无风险费( r 0 r_0 r0 = 5%)。
已知 n = 4 时相关数据如下表所列:
s i s_i si | r i r_i ri/% | q i q_i qi/% | p i p_i pi/% | u i u_i ui/元 |
---|---|---|---|---|
s 1 s_1 s1 | 28 | 2.5 | 1 | 103 |
s 2 s_2 s2 | 21 | 1.5 | 2 | 198 |
s 3 s_3 s3 | 23 | 5.5 | 4.5 | 52 |
s 4 s_4 s4 | 25 | 2.6 | 6.5 | 40 |
试给该公司涉及一种投资组合方案,即用给定资金 M ,有选择地购买若干种资产或存银行生息,使净收益尽可能大,使总体风险尽可能小。
2.2 符号规定和基本假设
- 符号规定
- s i s_i si 表示第 i 种投资项目,如股票、债券等, i = 0,1,···,n,其中 s 0 s_0 s0指存入银行;
- r i r_i ri, p i p_i pi, q i q_i qi 分别表示 s i s_i si 的平均收益率、交易费率、风险损失率, i = 0,1,···,n,其中 p 0 p_0 p0=0, q 0 q_0 q0=0;
- u i u_i ui 表示 s i s_i si 的交易定额, i = 0,1,···,n;
- x i x_i xi 表示投资项目 s i s_i si 的资金,i = 0,1,···,n;
- a 表示投资风险度;
- Q 表示总体收益。
- 基本假设
- 投资数额 M 相当大,为了便于计算,假设 M=1;
- 投资越分散,总的风险越少;
- 总体风险用投资项目 s i s_i si 中最大的一个风险来度量;
- n + 1 种资产 s i s_i si 之间时相互独立的;
- 在投资的这一时期内, r i r_i ri, p i p_i pi, q i q_i qi为定值,不受意外因素影响;
- 净收益和总体风险只受 r i r_i ri, p i p_i pi, q i q_i qi影响,不受其他因素影响。
2.3 模型的分析
总体风险用所投资的 s i s_i si 最大的一个风险来衡量,即:
max{ q i q_i qi x i x_i xi|i = 1,2,···,n }.
购买 s i s_i si (i=1,2,···,n) 所付交易是一个分段函数,即:
交 易 费 = { p i x i , x i > u i , p i u i , 0 < x i < u i , 交易费 = \begin{cases} p_i x_i , x_i \gt u_i,\\ p_i u_i , 0\lt x_i \lt u_i,&\\ \end{cases} 交易费={
pixi,xi>ui,piui,0<xi<ui,
而题目所给定的定值 u i u_i ui (单元: 元) 相对总投资 M 是很少的, p i p_i pi u i u_i ui更小,这样购买 s i s_i si 的净收益可以简化为( r i r_i ri - p i p_i pi) x i x_i xi。要使净收益尽可能大,总体风险尽可能小,这是一个多目标规划模型。
目标函数为:
{ m a x ∑ i = 0 n ( r i − p i ) x i , m i n { m a x ( 0 ≤ x i ≤ u i ) q i x i . \begin{cases} max \sum_{i=0}^n (r_i - p_i)x_i,\\ min \begin{cases} max_ {(0\le x_i \le u_i)} {q_ix_i}. \end{cases} \end{cases} {
max∑i=0n(ri−pi)xi,min{
max(0≤xi≤ui)qixi.
约束条件为:
{ ∑ i = 0 n ( i + p i ) x i = M , x i ≥ 0 , i = 1 , 2 , ⋅ ⋅ ⋅ , n . \begin{cases} \sum_{i=0}^n (i + p_i)x_i = M, \\ x_i \ge 0 ,&i=1,2,···,n.\\ \end{cases} {
∑i=0n(i+pi)xi=M,xi≥0,i=1,2,⋅⋅⋅,n.
2.4 模型的建立
模型简化:
在实际投资中,投资者承受风险的程度不一样,若给定风险的一个界限 a,使得最大风险 ( q i x i M ) ≤ a \left(\frac{q_ix_i}{M}\right) \le a (Mqixi)≤a ,可找到相应的投资方案。这样把多目标规划变成一个目标的线性规划。
模型一:固定风险水平水平、优化收益
max ∑ i = 0 n ( r i − p i ) x i , \sum_{i=0}^n(r_i - p_i)x_i, ∑i=0n(ri−pi)xi,
{ ( q i x i M ) ≤ a , i = 1 , 2 , ⋅ ⋅ ⋅ , n , ∑ i = 0 n ( 1 + p i ) x i = M , x i ≥ 0 , i = 1 , 2 , ⋅ ⋅ ⋅ , n . \begin{cases} \left(\frac{q_ix_i}{M}\right) \le a,&i=1,2,···,n,\\ \sum_{i=0}^n(1 + p_i)x_i = M,&x_i \ge0,i=1,2,···,n.\\ \end{cases} {
(Mqixi)≤a,∑i=0n(1+pi)xi=M,i=1,2,⋅⋅⋅,n,xi≥0,i=1,2,⋅⋅⋅,n.
带入数据如下:
min f = [-0.05,-0.27,-0.19,-0.185,-0.185]·[ x 0 , x 1 , x 2 , x 3 , x 4 x_0,x_1,x_2,x_3,x_4 x0,x1,x2,x3,x4] T ^T T
{ x 0 + 1.01 x 1 + 1.02 x 2 + 1.045 x 3 + 1.065 x 4 = 1 , 0.025 x 1 ≤ a , 0.015 x 2 ≤ a , 0.055 x 3 ≤ a , 0.026 x 4 ≤ a , x i ≤ 0 , i = 0 , 1 , ⋅ ⋅ ⋅ , 4. \begin{cases} x_0+1.01x_1+1.02x_2+1.045x_3+1.065x_4=1,\\ 0.025x_1 \le a,\\ 0.015x_2 \le a,\\ 0.055x_3 \le a,\\ 0.026x_4 \le a,\\ x_i \le 0, i=0,1,···,4. \end{cases} ⎩⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎧x0+1.01x1+1.02x2+1.045x3+1.065x4=1,0.025x1≤a,0.015x2≤a,0.055x3≤a,0.026x4≤a,xi≤0,i=0,1,⋅⋅⋅,4.
模型二:固定盈利水平、极小化风险
min ( m a x 1 ≤ i ≤ n max_{1 \le i \le n} max1≤i≤n( q i x i q_ix_i qixi))
s . t . { ∑ i = 0 n ( r i − p i ) x i ≥ k , ∑ i = 0 n ( 1 + p i ) x i = M , x i ≥ 0 , i = 0 , 1 , ⋅ ⋅ ⋅ , n . s.t. \begin{cases} \sum_{i=0}^n(r_i-p_i)x_i \ge k,\\ \sum_{i=0}^n(1+p_i)x_i = M,\\ x_i \ge 0, i = 0,1,···,n. \end{cases} s.t.⎩⎪⎨⎪⎧∑i=0n(ri−pi)xi≥k,∑i=0n(1+pi)xi=M,xi≥0,i=0,1,⋅⋅⋅,n.
设
x n + 1 x_{n+1} xn+1 = m a x 1 ≤ i ≤ n max_{1 \le i \le n} max1≤i≤n( q i x i q_ix_i qixi)则模型可以线性化为:
min x n + 1 x_{n+1} xn+1,
s . t . { q i x i ≥ x n + 1 , i = 1 , 2 , ⋅ ⋅ ⋅ , n , ∑ i = 0 n ( r i − p i ) x i ≥ k , ∑ i = 0 n ( 1 + p i ) x i = M , x i ≥ 0 , i = 0 , 1 , ⋅ ⋅ ⋅ , n . s.t. \begin{cases} \ q_ix_i \ge x_{n+1},&i = 1,2,···,n,\\ \sum_{i=0}^n(r_i-p_i)x_i \ge k,\\ \sum_{i=0}^n(1+p_i)x_i = M,\\ x_i \ge 0, i = 0,1,···,n. \end{cases} s.t.⎩⎪⎪⎪⎨⎪⎪⎪⎧ qixi≥xn+1,∑i=0n(ri−pi)xi≥k,∑i=0n(1+pi)xi=M,xi≥0,i=0,1,⋅⋅⋅,n.i=1,2,⋅⋅⋅,n,
带入数据如下:
min x 5 x_5 x5
s . t . { 0.025 x 1 − x 5 ≤ 0 , 0.015 x 1 − x 5 ≤ 0 , 0.055 x 3 − x 5 ≤ 0 , 0.026 x 4 − x 5 ≤ 0 , − 0.005 x 0 − 0.27 x 1 − 0.19 x 2 − 0.185 x 3 − 0.185 x 4 ≤ − k , 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 , ⋅ ⋅ ⋅ , 5. s.t. \begin{cases} 0.025x_1-x_5 \le0,\\ 0.015x_1-x_5 \le0,\\ 0.055x_3-x_5 \le0,\\ 0.026x_4-x_5 \le0,\\ -0.005x_0-0.27x_1-0.19x_2-0.185x_3-0.185x_4 \le -k,\\ x_0+1.01x_1+1.02x_2+1.045x_3+1.065x_4 = 1.\\ x_i \ge0,i = 0,1,···,5. \end{cases} s.t.⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧0.025x1−x5≤0,0.015x1−x5≤0,0.055x3−x5≤0,0.026x4−x5≤0,−0.005x0−0.27x1−0.19x2−0.185x3−0.185x4≤−k,x0+1.01x1+1.02x2+1.045x3+1.065x4=1.xi≥0,i=0,1,⋅⋅⋅,5.
模型三:两个目标函数加权求和
min w( m a x 1 ≤ i ≤ n max_{1 \le i \le n} max1≤i≤n ( q i x i ) (q_ix_i) (qixi) ) - (1 - w) ∑ i = 0 n ( r i − p i ) x i \sum_{i=0}^n(r_i - p_i)x_i ∑i=0n(ri−pi)xi,
s . t . { ∑ i = 0 n ( 1 + p i ) x i = m , x i ≥ 0 , i = 0 , 2 , ⋅ ⋅ ⋅ , n . s.t. \begin{cases} \sum_{i=0}^n(1+p_i)x_i = m,\\ x_i \ge 0,i = 0,2,···,n.\\ \end{cases} s.t.{
∑i=0n(1+pi)xi=m,xi≥0,i=0,2,⋅⋅⋅,n.
带入数据如下:
设
x n + 1 x_{n+1} xn+1 = m a x 1 ≤ i ≤ n max_{1 \le i \le n} max1≤i≤n( q i x i q_ix_i qixi)则模型可以线性化:
min (w x n + 1 x_{n+1} xn+1 - (1-w)[-0.05,-0.27,-0.19,-0.185,-0.185]·[ x 0 , x 1 , x 2 , x 3 , x 4 x_0,x_1,x_2,x_3,x_4 x0,x1,x2,x3,x4] T ^T T)
s . t . { 0.025 x 1 − x 5 ≤ 0 , 0.015 x 1 − 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 , ⋅ ⋅ ⋅ , 5. s.t. \begin{cases} 0.025x_1-x_5 \le0,\\ 0.015x_1-x_5 \le0,\\ 0.055x_3-x_5 \le0,\\ 0.026x_4-x_5 \le0,\\ x_0+1.01x_1+1.02x_2+1.045x_3+1.065x_4 = 1.\\ x_i \ge0,i = 0,1,···,5. \end{cases} s.t.⎩⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎧0.025x1−x5≤0,0.015x1−x5≤0,0.055x3−x5≤0,0.026x4−x5≤0,x0+1.01x1+1.02x2+1.045x3+1.065x4=1.xi≥0,i=0,1,⋅⋅⋅,5.
2.5 模型一的求解和分析
2.5.1 (代码)求解模型一
由于 a 是任意给定的风险度,所以不同的投资者有不同的风险度。所以定从
a = 0 开始,以步长 Δ \Delta Δa = 0.001 进行循环搜索,程序取下:
from matplotlib import pyplot as plt
from numpy import ones, diag, c_, zeros
from scipy.optimize import linprog
c = [-0.05,-0.27,-0.19,-0.185,-0.185]
A = c_[zeros(4),diag([0.025,0.015,0.055,0.026])]
Aeq = [[1,1.01,1.02,1.045,1.065]]; beq = [1]
a = 0; aa = []; ss = []
while a<0.05:
b = ones(4)*a
res = linprog(c,A,b,Aeq,beq)
x = res.x; Q = -res.fun
aa.append(a); ss.append(Q)
a = a + 0.001
#设置图形大小
plt.figure(figsize=(20,8),dpi=80)
#plt.plot(aa,ss,'r*')
plt.scatter(aa,ss,color="#DB7093")
plt.xlabel('a'); plt.ylabel('Q',rotation=90)
plt.savefig("模块一.png")
plt.show()
2.5.2 模型一(结果)分析
由上图可知:
(1)风险大,收入也会大。
(2)当投资越分散时,投资者承担的风险越小
,这与题意一致。冒险的投资者会出现集中投资的情况,保守的投资者则尽量分散投资。
(3)在 a = 0.006 附近有一个转折点,在这一点左边,风险增加很少时,李瑞增长很快。在这一点右边,风险增加很大时,利润增长很缓慢,所以对于风险和收益没有特殊偏好的投资者来说,应该选择曲线的转折点作为最优投资组合,大约是 a = 0.6%,Q = 20%,所对应投资方案为:
所对应投资方案为:
风险度 a = 0.006 ; 收益 Q = 0.2019;
$x_0$=0, $x_1$=0.24, $x_2$=0.4, $x_3$=0.1091, $x_4$=0.2212;
2.6 模型二的求解和分析
2.6.1 (代码)求解模型二
from matplotlib import pyplot as plt
import numpy as np
import cvxpy as cp
x = cp.Variable(6,pos=True) # 决策向量6个,且都为正数
obj = cp.Minimize(x[5]) # 优化得目标函数(构造目标函数)
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])
k = 0.05; kk = []; ss = []
while k<0.27:
con = [cp.multiply(a1, x[1:5])-x[5]<=0, a2@x[:-1]>=k, a3@x[:-1]==1]
prob = cp.Problem(obj,con) #把目标函数与约束函数传进Problem函数种
prob.solve(solver='GLPK_MI') #返回优化值
#print("最优值为:\n", prob.value)
#print("最优解为:\n", x.value)
kk.append(k); ss.append(prob.value)
k = k+0.005
#设置图形大小
plt.figure(figsize=(20,8),dpi=80)
plt.plot(kk,ss,'r*')
#plt.scatter(kk,ss,color="#DB7093")
plt.xlabel('k')
plt.ylabel('R',rotation=90)
plt.savefig("模块二.png")
plt.show()
2.6.2 模型二(结果)分析
- 当收益 k = 0.21 时, x 0 x_0 x0 = 0, x 1 x_1 x1 = 0.3089, x 2 x_2 x2 = 0.5148, x 3 x_3 x3 = 0.1404, x 4 x_4 x4 = 0.0152;
- 风险最小为 0.0077.
2.7 模型三的求解和分析
2…1 (代码)求解模型三
import numpy as np
import cvxpy as cp
import pylab as plt
from matplotlib import font_manager
my_font = font_manager.FontProperties(fname="C:/windows/Fonts/simsun.ttc")
x = cp.Variable(6,pos=True)
r = np.array([0.05,0.28,0.21,0.23,0.25])
p = np.array([0,0.01,0.02,0.045,0.065])
q = np.array([0,0.025,0.015,0.055,0.026])
def LP(w):
V = [] #风险初始化
Q = [] #收益初始化
X = [] #最优解的初始化
con = [(1+p) @ x[:-1] == 10000] # 这里M=10000和1时,输出没有区别,就是一个比例问题。
for i in range (1,5):
con.append(q[i] * x[i] <= x[5])
for i in range(len(w)):
obj = cp.Minimize(w[i] * x[5] - (1-w[i]) * ((r-p) @ x[:-1]))
prob = cp.Problem(obj, con)
prob.solve(solver='GLPK_MI')
print("最优解为:\n", x.value)
print("最优值为:\n", prob.value)
xx = x.value #提出所有决策变量的取值
V.append(max(q * xx[:-1]))
Q.append((r-p) @ xx[:-1])
X.append(xx)
print('w=',w)
print('V=',np.round(V,2))
print('Q=',np.round(Q,2))
plt.figure()
plt.plot(V,Q,'*-')
plt.grid('on')
plt.xlabel('风险(元)',fontproperties=my_font)
plt.ylabel('收益(元)',fontproperties=my_font)
return x
w1 = np.arange(0,1.1,0.1)
LP(w1); print("-----------------------")
w2 = np.array([0.766,0.767,0.810,0.811,0.824,0.825,0.962,0.963,1.0])
X = LP(w2)
print(X[-3])
plt.show()
OUT:
(问什么我要把这么多结果输出呢?当然是下面分析是能用上的啦
)
Long-step dual simplex will be used
最优解为:
[ 0. 9900.99009901 -0. -0. -0.
247.52475248]
最优值为:
-2673.2673267326736
Long-step dual simplex will be used
最优解为:
[ 0. 9900.99009901 -0. -0. -0.
247.52475248]
最优值为:
-2381.188118811882
Long-step dual simplex will be used
最优解为:
[ 0. 9900.99009901 -0. -0. -0.
247.52475248]
最优值为:
-2089.1089108910896
Long-step dual simplex will be used
最优解为:
[ 0. 9900.99009901 -0. -0. -0.
247.52475248]
最优值为:
-1797.0297029702972
Long-step dual simplex will be used
最优解为:
[ 0. 9900.99009901 -0. -0. -0.
247.52475248]
最优值为:
-1504.9504950495052
Long-step dual simplex will be used
最优解为:
[ 0. 9900.99009901 -0. -0. -0.
247.52475248]
最优值为:
-1212.871287128713
Long-step dual simplex will be used
最优解为:
[ 0. 9900.99009901 -0. -0. -0.
247.52475248]
最优值为:
-920.7920792079208
Long-step dual simplex will be used
最优解为:
[ 1.93267624e-12 9.90099010e+03 0.00000000e+00 -0.00000000e+00
-0.00000000e+00 2.47524752e+02]
最优值为:
-628.7128712871287
Long-step dual simplex will be used
最优解为:
[ 0. 3690.03690037 6150.06150062 -0. -0.
92.25092251]
最优值为:
-359.16359163591625
Long-step dual simplex will be used
最优解为:
[ 0. 2375.83953945 3959.73256575 1079.92706339 2284.46109563
59.39598849]
最优值为:
-148.16737761865005
Long-step dual simplex will be used
最优解为:
[10000. -0. -0. -0. -0. -0.]
最优值为:
-0.0
w= [0. 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1. ]
V= [247.52 247.52 247.52 247.52 247.52 247.52 247.52 247.52 92.25 59.4
0. ]
Q= [2673.27 2673.27 2673.27 2673.27 2673.27 2673.27 2673.27 2673.27 2164.82
2016.24 500. ]
-----------------------
Long-step dual simplex will be used
最优解为:
[ 1.93267624e-12 9.90099010e+03 0.00000000e+00 -0.00000000e+00
-0.00000000e+00 2.47524752e+02]
最优值为:
-435.94059405940595
Long-step dual simplex will be used
最优解为:
[ 0. 3690.03690037 6150.06150062 -0. -0.
92.25092251]
最优值为:
-433.6469864698647
Long-step dual simplex will be used
最优解为:
[ 0. 3690.03690037 6150.06150062 -0. -0.
92.25092251]
最优值为:
-336.5928659286592
Long-step dual simplex will be used
最优解为:
[ 2.27373675e-13 3.13971743e+03 5.23286238e+03 1.42714428e+03
-0.00000000e+00 7.84929356e+01]
最优值为:
-334.37419723133996
Long-step dual simplex will be used
最优解为:
[ 2.27373675e-13 3.13971743e+03 5.23286238e+03 1.42714428e+03
-0.00000000e+00 7.84929356e+01]
最优值为:
-305.9759288330718
Long-step dual simplex will be used
最优解为:
[ 0. 2375.83953945 3959.73256575 1079.92706339 2284.46109563
59.39598849]
最优值为:
-303.83990219737484
Long-step dual simplex will be used
最优解为:
[4.54747351e-13 2.37583954e+03 3.95973257e+03 1.07992706e+03
2.28446110e+03 5.93959885e+01]
最优值为:
-19.47809063357112
Long-step dual simplex will be used
最优解为:
[10000. -0. -0. -0. -0. -0.]
最优值为:
-18.500000000000018
Long-step dual simplex will be used
最优解为:
[10000. -0. -0. -0. -0. -0.]
最优值为:
-0.0
w= [0.766 0.767 0.81 0.811 0.824 0.825 0.962 0.963 1. ]
V= [247.52 92.25 92.25 78.49 78.49 59.4 59.4 0. 0. ]
Q= [2673.27 2164.82 2164.82 2105.99 2105.99 2016.24 2016.24 500. 500. ]
2.6.2 模型三(结果)分析
2…6.2.1 分析输出结果一
可以得到当w取不同值时风险和收益的计算结果如下所示:
表一:
w | 0.0 | 0.1 | 0.2 | 0.3 | 0.4 | 0.5 | 0.6 | 0.7 | 0.8 | 0.9 | 1.0 |
---|---|---|---|---|---|---|---|---|---|---|---|
风险 | 247.52 | 247.52 | 247.52 | 247.52 | 247.52 | 247.52 | 247.52 | 247.52 | 92.25 | 59.4 | 0 |
收益 | 2673.27 | 2673.27 | 2673.27 | 2673.27 | 2673.27 | 2673.27 | 2673.27 | 2673.27 | 2164.82 | 2016.24 | 500 |
从以上数据可以看出,当投资偏好系数 w ≤ \le ≤ 0.7,所对应的风险收益和风险均达到最大值。此时,收益为 2673.27 元,风险为 247.52 元,此时应将全部资金均用来购买资产 s 1 s_1 s1;当 w 由 0.7 增加到 1.0 时,收益和风险均呈现下降趋势,特别,当 w = 1.0 时,收益和风险达到最小值,收益为 500 元,风险为 0 元,此时应将所有资金全部存入银行。
解答疑惑:上面的分析为啥选择方案
s 1 s_1 s1呢?:
答:这就是上面第一种输出结果的最优方案:也就是:
min w( m a x 1 ≤ i ≤ n max_{1 \le i \le n} max1≤i≤n ( q i x i ) (q_ix_i) (qixi) ) - (1 - w) ∑ i = 0 n ( r i − p i ) x i \sum_{i=0}^n(r_i - p_i)x_i ∑i=0n(ri−pi)xi等于0时,即最优值为( -0.0)时,又因为 x i x_i xi 表示投资项目 s i s_i si的资金,即为最优解 [10000. -0. -0. -0. -0. -0.]。这就是上面代码输出结果为那么多,我还要把它复制到博文中。
最优解为:
[10000. -0. -0. -0. -0. -0.]
最优值为:
-0.0
w= [0. 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1. ]
V= [247.52 247.52 247.52 247.52 247.52 247.52 247.52 247.52 92.25 59.4
0. ]
Q= [2673.27 2673.27 2673.27 2673.27 2673.27 2673.27 2673.27 2673.27 2164.82
2016.24 500. ]
2.6.1.2 分析输出结果二
为了更好地描述收益和风险的对应关,可将 w 取值进一步细化,重新计算的部分数据如下表二所示。
w | 0.766 | 0.767 | 0.810 | 0.811 | 0.824 | 0.825 | 0.962 | 0.963 | 1.0 |
---|---|---|---|---|---|---|---|---|---|
风险 | 247.52 | 92.25 | 95.25 | 78.49 | 78.49 | 59.4 | 59.4 | 0 | 0 |
收益 | 2673.27 | 2164.82 | 2164.82 | 2105.99 | 2105.99 | 2016.24 | 2016.24 | 500 | 500 |
从第二幅图可以看出,投资得收益越大,风险越大。投资者可以更具自己对风险喜好得不同,选择合适得投资方案。曲线的拐点坐标约为(59.4,2016.24),此时对应得投资方案时购买资产 s 1 s_1 s1、 s 2 s_2 s2、 s 3 s_3 s3、 s 4 s_4 s4的资金分别为2375.84元,3959.73元,1079.93元、2284.46元。存入银行的资金为 0 元,这对于风险和收益没有明显偏好的投资者是一个比较合适的选择。
这个投资方案怎么分析出来的:
(曲线的拐点坐标约为(59.4,2016.24),此时对应得投资方案时购买资产 s 1 s_1 s1、 s 2 s_2 s2、 s 3 s_3 s3、 s 4 s_4 s4的资金分别为2375.84元,3959.73元,1079.93元、2284.46元。)
答:因为我们在前面设
x n + 1 x_{n+1} xn+1 = m a x 1 ≤ i ≤ n max_{1 \le i \le n} max1≤i≤n( q i x i q_ix_i qixi)使模型可以线性化:
,又因为 m a x 1 ≤ i ≤ n max_{1 \le i \le n} max1≤i≤n( q i x i q_ix_i qixi)表示最大风险,n = 4,所以 x 5 x_5 x5 = m a x 1 ≤ i ≤ 5 max_{1 \le i \le 5} max1≤i≤5( q i x i q_ix_i qixi) =59.4(元)时,对应的是输出中
最优解为(为了便于理解我在下面结果加上对应的符号):
[ x 0 x_0 x0 = 0. , x 1 x_1 x1 = 2375.83953945 , x 2 x_2 x2 =3959.73256575 , x 3 x_3 x3 =1079.92706339 , x 4 x_4 x4 =2284.46109563 , x 5 x_5 x5 = 59.39598849]
最优解为:
[ 0. 2375.83953945 3959.73256575 1079.92706339 2284.46109563
59.39598849]
最优值为:
-303.83990219737484
(三)尝试用cvxopt.solvers 模块解题
提出问题:对于本文2.6.1 (代码)求解模型二
用的是cvxpy库
方法求解的,我们是否可以用cvxopt.solvers 模块或者pulp库
求解呢!
下面我们看一下用cvxopt.solvers 模块
来求解模型二,代码取下:
import numpy as np
from matplotlib import pyplot as plt
from cvxopt import matrix,solvers
c = matrix([0.,0,0,0,0])
A = matrix([[0.025,0,0,0,1],[0.015,0,0,0,1],[0.055,0,0,0,1],[0.026,0,0,0,1],[-0.05,-0.27,-0.19,-0.185,-0.185]]).T # matrix([2.,1]):记住matrix中第一个元素不许是float型
k = 0.05;kk = []
ss = []
while k < 0.27:
b = matrix([0., 0, 0, 0, -k])
Aeq = matrix([1., 1.01, 1.02, 1.045, 1.065], (1, 5)) # Aeq为行向量
beq = matrix(1.)
sol = solvers.lp(c, A, b, Aeq, beq)
kk.append(k); ss.append(sol['primal objective'])
k = k + 0.005
#设置图形大小
plt.figure(figsize=(20,8),dpi=80)
#plt.plot(kk,ss,'r*')
plt.scatter(kk,ss,color="#DB7093")
plt.xlabel('k')
plt.ylabel('R',rotation=90)
plt.savefig("cvxopt.solvers求解模块二.png")
plt.show()
结果\问题分析:
与之前cvxpy库
求解模型二得到的图形完全不一样,本人尝试过pupl库求解的输出图形和cvxopt.solvers 模块一样。为什么呢!上一篇博文python数学建模(一)线性规划.1我们介绍了四种求解线性规划的问题的方法(Scipy、pulp、cvxopt.solvers 、cvxpy),难道cvxpy库求解模型二是错误的吗?
答:显然不是,而是这里不能用Scipy、pulp、cvxopt.solvers 三种方法求解模型,因为cvxpy库适用于凸优化问题方面的线性问题,而Scipy库、pupl库和solver模块是用于标准线性规划问题的,所以cvxpy库与Scipy库、pupl库、solver模块是有很大区别的。(所以当遇到求解凸优化的题用cvxpy库,标准线性问题可以用Scipy库、pupl库、solver模块。
)
(四)修改上一篇博文中的一个题目
由于本人的失误在上一篇博文python数学建模(一)线性规划.1 中目录 (3.2.3 cvxopt.solvers 模块求解)中题目的代码复制成其他题目的代码了,但是两个提都属于cvxopt.solvers 模块求解
问题,现在都给大家补充上。
4.1 补充题一
min z = − 4 x 1 − 5 x 2 -4x_1 - 5x_2 −4x1−5x2,
s . t . { 2 x 1 + x 2 ≤ 3 , x 1 + 2 x 2 ≤ 3 , x 1 ≥ 0 , x 2 ≥ 0. s.t. \begin{cases} 2x_1 + x_2 \le 3,\\ x_1 + 2x_2 \le 3,&\\ x_1 \ge 0, x_2 \ge 0. \end{cases} s.t.⎩⎪⎨⎪⎧2x1+x2≤3,x1+2x2≤3,x1≥0,x2≥0.
import numpy as np
from cvxopt import matrix,solvers
c = matrix([-4.,-5]); A = matrix([[2.,1],[1,2],[-1,0],[0,-1]]).T
b = matrix([3.,3,0,0]); sol = solvers.lp(c,A,b)
print("最优解为:\n",sol['x'])
print("最优值为:\n",sol['primal objective'])
OUT
pcost dcost gap pres dres k/t
0: -8.1000e+00 -1.8300e+01 4e+00 0e+00 8e-01 1e+00
1: -8.8055e+00 -9.4357e+00 2e-01 2e-16 4e-02 3e-02
2: -8.9981e+00 -9.0049e+00 2e-03 6e-16 5e-04 4e-04
3: -9.0000e+00 -9.0000e+00 2e-05 1e-16 5e-06 4e-06
4: -9.0000e+00 -9.0000e+00 2e-07 1e-16 5e-08 4e-08
Optimal solution found.
最优解为:
[ 1.00e+00]
[ 1.00e+00]
最优值为:
-8.999999811406719
4.2 补充题一
min z = − 2 x 1 − 5 x 2 -2x_1 - 5x_2 −2x1−5x2,
s . t . { − x 1 + x 2 ≤ 1 , x 1 + x 2 ≥ 2 , x 1 − 2 x 2 ≤ 4 , x 2 ≥ 0 , x 1 + x 2 = 3.5. s.t. \begin{cases} -x_1 + x_2 \le 1,\\ x_1 + x_2 \ge 2,&\\ x_1 -2 x_2 \le 4,\\ x_2 \ge 0,\\ x_1 + x_2 = 3.5.\\ \end{cases} s.t.⎩⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎧−x1+x2≤1,x1+x2≥2,x1−2x2≤4,x2≥0,x1+x2=3.5.
# 线性规划 使用cvxopy.solvers模块求解
import numpy
from cvxopt import matrix,solvers
c = matrix([2.,1]); A = matrix([[-1.,1],[-1,-1],[1,-2],[0,-1]]).T # matrix([2.,1]):记住matrix中第一个元素不许是float型
b = matrix([1.,-2,4,0]); Aeq = matrix([1.,2],(1,2)) # Aeq为行向量
beq = matrix(3.5); sol = solvers.lp(c,A,b,Aeq,beq)
print("最优解为:\n",sol['x'])
print("最优值为:\n",sol['primal objective'])
OUT:
D:\ProgramData\Anaconda3\python.exe C:/Users/Administrator/PycharmProjects/pythonProject/数学建模/线性规划/线性规划3.cvxopy.solvers.py
pcost dcost gap pres dres k/t
0: 5.5556e+00 1.2222e+00 1e+01 0e+00 2e+00 1e+00
1: 4.6038e+00 3.7995e+00 2e+00 1e-16 4e-01 2e-01
2: 2.5229e+00 2.4639e+00 2e-01 2e-16 4e-02 4e-02
3: 2.5002e+00 2.4997e+00 2e-03 1e-16 4e-04 4e-04
4: 2.5000e+00 2.5000e+00 2e-05 4e-16 4e-06 4e-06
5: 2.5000e+00 2.5000e+00 2e-07 1e-16 4e-08 4e-08
Optimal solution found.
最优解为:
[ 5.00e-01]
[ 1.50e+00]
最优值为:
2.500000024611047