数学实验9结23李会平.docx
《数学实验9结23李会平.docx》由会员分享，可在线阅读，更多相关《数学实验9结23李会平.docx（24页珍藏版）》请在冰豆网上搜索。
数学实验9结23李会平
实验九非线性规划
【实验目的】
1.掌握用matlab优化工具箱和lingo解非线性规划的方法；
2.练习建立实际问题的非线性规划模型
【实验内容】
9.4生产决策问题
一、模型建立
根据题意假设分别有
吨（t）的甲乙丙原料用于生产A产品，有
吨（t）的甲乙丙原料用于生产B产品。
由于甲乙先混合，再和丙混合，所以在A,B产品中，甲乙的比例相同的，本题需要在满足原料和需求量约束的条件下最大化理论，即
st:
整理以后变成：
st:
二、模型求解
1.用matlab实现
（1）在matlab中通过fmincon函数实现的编程如下：
目标函数：
function[f]=yetifun1（x）
c=[-3,7,1,-9,1,-5];
f=c*x';
end
约束函数：
function[c,ceq]=yetifun2（x）
c=0;
ceq=x
（1）*x（5）-x（4）*x
（2）;
end
主函数：
%9.4
（1）问题
clc,clearall
x0=[20,20,20,20,30,30];
A1=[0.5,-1.5,-0.5,0,0,0;...
0,0,0,1.5,-0.5,0.5;...
eye（3）,eye（3）;...
1,1,1,0,0,0;...
0,0,0,1,1,1];
b1=[0,0,500,500,500,100,200];
v1=zeros（1,6）;
opt=optimset（'largescale','off','maxiter',3000,'maxfun',20000）;
[x,fv,ef,out,lag,grad,hess]=fmincon（@yetifun1,x0,A1,b1,[],[],v1,[],@yetifun2,opt）
运行后得到结果：
x=0.0000-0.00000.00000100.0000100.0000
fv=-400.0000
ef=1
另外，ef=1说明结果是收敛的，据此可以决定生产决策，用100t乙和丙去生产B产品，最后得到的利润为400k即40万。
（2）A的需求量增长为600t时候，重新进行规划，只需要将右端项量的b1由b1=[0,0,500,500,500,100,200]改成[0,0,500,500,500,600,200]，运行x的值对初值极其敏感，为此编写如下代码：
%9.4
（2）,A的需求量增长为600t
clc,clearall
x0=[300,0,300,0,0,0];
A1=[0.005,-0.015,-0.005,0,0,0;...
0,0,0,0.015,-0.005,0.005;...
eye（3）,eye（3）;...
1,1,1,0,0,0;...
0,0,0,1,1,1];
b1=[0,0,500,500,500,600,200];
v1=zeros（1,6）';
opt=optimset（'largescale','off','maxiter',3000000,'maxfun',2000000）;
[x,fv,ef,out,lag,grad,hess]=fmincon（@yetifun1,x0,A1,b1,[],[],v1,[],@yetifun2,opt）
得到的结果是：
x=3000300000
fv=-600
ef=1
分析：
用300t的甲和300t的丙生产A是最优的，当然此处的最优结果是通过不断变换初值得到的结果，相比之下，用lingo寻找结果则好得多。
（3）如果B的价格下降为13千元/t，此时将原目标函数
变为：
对费用向量的系数进行调整以后即目标函数的代码变为：
function[f]=yetifun1（x）
c=[-3,4,1,-9,-2,-5];
f=c*x';
end
得到的结果为：
x=0.0000-0.00000.000050.0000150.00000
fv=-750
ef=1
分析：
生产决策变成将50t甲和150t乙用于生产B产品，最大利润是75万。
2.lingo实现
matlab难以给出问题的全局最优解，而lingo在这方面做得更好，所以对以上的问题均通过lingo实现与matlab做对比：
（1）lingo实现的代码相对简单：
min=-3*x1+7*x2+x3-9*x4+x5-5*x6;
0.005*x1-0.015*x2-0.005*x3<=0;
0.015*x4-0.005*x5+0.005*x6<=0;
x1+x4<=500;
x2+x5<=500;
x3+x6<=500;
x1+x2+x3<=100;
x4+x5+x6<=200;
x1*x5=x4*x2;
得到结果如下所示：
Globaloptimalsolutionfound.
Objectivevalue:
-400.0000
Objectivebound:
-400.0000
Infeasibilities:
0.000000
Extendedsolversteps:
2
Totalsolveriterations:
140
ModelClass:
NLP
Totalvariables:
6
Nonlinearvariables:
4
Integervariables:
0
Totalconstraints:
9
Nonlinearconstraints:
1
Totalnonzeros:
28
Nonlinearnonzeros:
4
VariableValueReducedCost
X10.0000000.000000
X20.0000004.000000
X30.0000000.000000
X40.0000002.000000
X5100.00000.000000
X6100.00000.000000
RowSlackorSurplusDualPrice
1-400.0000-1.000000
20.000000200.0000
30.000000600.0000
4500.00000.000000
5400.00000.000000
6400.00000.000000
7100.00000.000000
80.0000002.000000
90.0000000.2000000E-01
分析：
得到了全局最优解，结果与matlab在
（1）问中得到的局部最优解是一致的。
（2）A的需求量变成600时候，仿照上例中lingo的代码，将100改成600，仍然可以得到全局最优解，如下所示：
VariableValueReducedCost
X1300.00000.000000
X20.0000002.000000
X3300.00000.000000
X40.0000006.000000
X50.0000000.000000
X60.0000000.000000
RowSlackorSurplusDualPrice
1-600.0000-1.000000
20.000000400.0000
30.0000001000.000
4200.00000.000000
5500.00000.000000
6200.00000.000000
70.0000001.000000
8200.00000.000000
90.0000000.1333333E-01
分析：
应该用300t甲和300t丙去生产A，最大利润为60万，这与matlab得到的局部最优解一样，验证可知该解在可行域内。
在寻找全局最优解时候用lingo更方便。
同时看到，市场需求的变化对生产决策有着重要的影响，本题中，因为对A的需求大幅增加，所以生产商将放弃B的生产，而只生产A。
（3）B的价格将为13时候，此时的lingo程序也只需将目标函数变化一下，代码为：
min=-3*x1+4*x2+x3-9*x4-2*x5-5*x6;
0.005*x1-0.015*x2-0.005*x3<=0;
0.015*x4-0.005*x5+0.005*x6<=0;
x1+x4<=500;
x2+x5<=500;
x3+x6<=500;
x1+x2+x3<=600;
x4+x5+x6<=200;
x1*x5=x4*x2;
运行后得到的全局最优解如下：
VariableValueReducedCost
X10.0000001.000000
X20.0000000.000000
X30.0000000.000000
X450.000000.000000
X5150.00000.000000
X60.0000000.5000000
RowSlackorSurplusDualPrice
1-750.0000-1.000000
20.000000200.0000
30.000000350.0000
4450.00000.000000
5350.00000.000000
6500.00000.000000
7100.00000.000000
80.0000003.750000
90.0000000.2000000E-01
分析：
这与matlab的局部最优解一样，但是总的来说，用lingo求解更加方便，因为lingo往往可以直接给出问题的全局最优解且编程实现更直观。
9.8投资组合问题
一、模型建立
1.本题为投资组合模型的具体应用，要找出投资方案，首先要得到A、B、C三种股票的期望收益率Er（平均值可以作为它的无偏估计），然后是三种股票的方差和协方差，对于股票i,j他们的方差及协方差在概率统计中都有相应的无偏估计形式，即
以上i，j为股票类别，k为各个时期，n为总时期数。
特别的，当i=j时候，得到的结果为方差计算方法。
在lingo中求期望收益率和协方差矩阵时候，需要通过以上公式变成计算，但是在matlab中直接提供了cov命令可以得到几个向量的协方差矩阵，大大简化了运算。
2.本体的目标即为在给定收益率情况下，让风险最小。
以原题要求为例即
st:
其中
，而
‘，cov为
向量的协方差矩阵。
（二）模型求解
1.matlab实现
（1）在matlab中编程如下所示：
①采用fmincon函数求解时：
1）目标函数：
functionf=gupiaofun1（x,cov）
f=x*cov*x';
end
2）命令文件：
clc,clearall
a1=[1.3001.1031.2160.9540.9291.0561.0381.0891.0901.0831.0351.176];
a2=a1-ones（1,12）;
b1=[1.2251.2901.2160.7281.1441.1071.3211.3051.1951.3900.9281.715];
b2=b1-ones（1,12）;
c1=[1.1491.2601.4190.9221.1690.9651.1331.7321.0211.1311.0061.908];
c2=c1-ones（1,12）;
cov=cov（[a2',b2',c2']）%将a1,b1,c1转化为列向量，用cov命令求协方差矩阵
Er=[mean（a2）,mean（b2）,mean（c2）]';
%期望年收益至少达到15%时候的投资，采用fmincon求解
x0=[1/3,1/3,1/3];
opt=optimset（'largescale','off','maxiter',3000,'maxfun',20000）;
[x,fv,ef,out,lag,grad,hess]=fmincon（@gupiaofun1,x0,-Er',-0.15,ones（1,3）,1,zeros（1,3）,ones（1,3）,[],opt,cov）
3）得到的结果如下所示：
x=0.53010.35640.1135
fv=0.0224
ef=1
分析：
结果是收敛的，且三种股票的持股比例分别为53.01%，35.64%，11.35%，
②如果采用quadprog命令求解，则命令文件如下所示：
clc,clearall
a1=[1.3001.1031.2160.9540.9291.0561.0381.0891.0901.0831.0351.176];
a2=a1-ones（1,12）;
b1=[1.2251.2901.2160.7281.1441.1071.3211.3051.1951.3900.9281.715];
b2=b1-ones（1,12）;
c1=[1.1491.2601.4190.9221.1690.9651.1331.7321.0211.1311.0061.908];
c2=c1-ones（1,12）;
cov=cov（[a2',b2',c2']）
Er=[mean（a2）,mean（b2）,mean（c2）]';
%用二次规划的方法求解
H=2*cov;%由于方差是x*cov*x’,所以在系数为0.5情况下，H2倍
x0=[1/3,1/3,1/3];
opt=optimset（'largescale','off','maxiter',3000,'maxfun',20000）;
[x,fv,ef,out,lag]=quadprog（H,zeros（1,3）,-Er',-0.15,ones（1,3）,1,zeros（1,3）,ones（1,3）,x0,opt）
得到结果为：
x=0.53010.35640.1135
fv=0.0224
ef=1
分析：
得到结果与用fmincon函数的结果一样，但是一般来说针对二次规划的问题，用quadprog函数往往会更好一点。
问题一：
如果年收益率在10%到100%之间变化，投资组合和风险如何变化；
解答：
针对不同的收益率，直接将原来的约束
的右边换成相应的约束即可，但是考虑到Er=（0.0891，0.2137，0.2346），即最高收益率为23.46%，所以更高的收益率期望已不可能实现，所以以下只分析10%到24%区间的情况，变成如下所示：
%期望收益率在10%到24%之间浮动的时候，去0.2个百分点为1个间隔分析
clc,clearall
a1=[1.3001.1031.2160.9540.9291.0561.0381.0891.0901.0831.0351.176];
a2=a1-ones（1,12）;
b1=[1.2251.2901.2160.7281.1441.1071.3211.3051.1951.3900.9281.715];
b2=b1-ones（1,12）;
c1=[1.1491.2601.4190.9221.1690.9651.1331.7321.0211.1311.0061.908];
c2=c1-ones（1,12）;
cov=cov（[a2',b2',c2']）;
Er=[mean（a2）,mean（b2）,mean（c2）]';
H=2*cov;
x0=[1/3,1/3,1/3];
opt=optimset（'largescale','off','maxiter',3000,'maxfun',20000）;
fori=0.1:
0.002:
0.24
k=round（（i-0.1）/0.002+1）;%表示在进行第几次运算
T（k,1）=i;
S（k,1）=i;
[x,fv,ef,out,lag]=quadprog（H,zeros（1,3）,-Er',-i,ones（1,3）,1,zeros（1,3）,ones（1,3）,x0,opt）;
T（k,2:
4）=x;
S（k,2）=fv;
end
T
plot（S（:
1）,S（:
2））
得到的结果为：
T=
0.10000.91450.07250.0129
0.10200.89920.08390.0170
0.10400.88380.09520.0210
0.10600.86840.10660.0250
0.10800.85300.11790.0290
0.11000.83770.12930.0330
0.11200.82230.14070.0371
0.11400.80690.15200.0411
0.11600.79150.16340.0451
0.11800.77610.17470.0491
0.12000.76080.18610.0532
0.12200.74540.19740.0572
0.12400.73000.20880.0612
0.12600.71460.22010.0652
0.12800.69930.23150.0692
0.13000.68390.24290.0733
0.13200.66850.25420.0773
0.13400.65310.26560.0813
0.13600.63770.27690.0853
0.13800.62240.28830.0894
0.14000.60700.29960.0934
0.14200.59160.31100.0974
0.14400.57620.32230.1014
0.14600.56080.33370.1055
0.14800.54550.34510.1095
0.15000.53010.35640.1135
0.15200.51470.36780.1175
0.15400.49930.37910.1215
0.15600.48400.39050.1256
0.15800.46860.40180.1296
0.16000.45320.41320.1336
0.16200.43780.42450.1376
0.16400.42240.43590.1417
0.16600.40710.44730.1457
0.16800.39170.45860.1497
0.17000.37630.47000.1537
0.17200.36090.48130.1578
0.17400.34560.49270.1618
0.17600.33020.50400.1658
0.17800.31480.51540.1698
0.18000.29940.52670.1738
0.18200.28400.53810.1779
0.18400.26870.54940.1819
0.18600.25330.56080.1859
0.18800.23790.57220.1899
0.19000.22250.58350.1940
0.19200.20720.59490.1980
0.19400.19180.60620.2020
0.19600.17640.61760.2060
0.19800.16100.62890.2100
0.20000.14560.64030.2141
0.20200.13030.65160.2181
0.20400.11490.66300.2221
0.20600.09950.67440.2261
0.20800.08410.68570.2302
0.21000.06870.69710.2342
0.21200.05340.70840.2382
0.21400.03800.71980.2422
0.21600.02260.73110.2463
0.21800.00720.74250.2503
0.220000.69720.3028
0.222000.60160.3984
0.22400.00000.50600.4940
0.226000.41040.5896
0.228000.31470.6853
0.230000.21910.7809
0.232000.12350.8765
0.234000.02790.9721
0.2360-0.0030-0.00001.0030
0.2380-0.0072-0.00001.0072
0.2400-0.0114-0.00001.0114
分析：
上面显示了在不同期望收益率情况下需要进行的投资比例的变化，注意到21.80%到22.00%之间发生了投资组合种类的变化，所以将该区间细分求出更精确的T矩阵如下如下所示：
T=
0.21800.00720.74250.2503
0.21810.00650.74310.2505
0.21820.00570.74360.2507
0.21830.00490.74420.2509
0.21840.00420.74480.2511
0.21850.00340.74530.2513
0.21860.00260.74590.2515
0.21870.00190.74650.2517
0.21880.00110.74700.2519
0.21890.00030.74760.2521
0.219000.74500.2550
0.219100.74020.2598
0.219200.73550.2645
0.219300.73070.2693
0.219400.72590.2741
0.219500.72110.2789
0.219600.71630.2837
0.219700.71160.2884
0.219800.70680.2932
0.219900.70200.2980
0.220000.69720.3028
分析：
所以当期望收益率达到21.90%时候就应该只投资B、C两种股票，这是因为B、C两种股票的收益率更高的原因，为了得到更大收益率，必须选择风险更高的B、C两种股票。
可以画出最小风险随着期望
升级会员 