数学实验9结23李会平.docx

1、数学实验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

2、,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,clear allx0=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,maxfu

3、n,20000);x,fv,ef,out,lag,grad,hess=fmincon(yetifun1,x0,A1,b1,v1,yetifun2,opt)运行后得到结果：x = 0.0000 -0.0000 0.0000 0 100.0000 100.0000fv = -400.0000ef = 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，运行

4、x的值对初值极其敏感，为此编写如下代码：%9.4(2),A的需求量增长为600tclc,clear allx0=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=

5、fmincon(yetifun1,x0,A1,b1,v1,yetifun2,opt) 得到的结果是：x = 300 0 300 0 0 0fv = -600ef = 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.0000 0.0000 50.0000 1

6、50.0000 0fv = -750ef = 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*x

7、5=x4*x2;得到结果如下所示：Global optimal solution found. Objective value: -400.0000 Objective bound: -400.0000 Infeasibilities: 0.000000 Extended solver steps: 2 Total solver iterations: 140 Model Class: NLP Total variables: 6 Nonlinear variables: 4 Integer variables: 0 Total constraints: 9 Nonlinear constra

8、ints: 1 Total nonzeros: 28 Nonlinear nonzeros: 4Variable Value Reduced CostX1 0.000000 0.000000X2 0.000000 4.000000X3 0.000000 0.000000X4 0.000000 2.000000X5 100.0000 0.000000X6 100.0000 0.000000Row Slack or Surplus Dual Price1 -400.0000 -1.0000002 0.000000 200.00003 0.000000 600.00004 500.0000 0.00

9、00005 400.0000 0.0000006 400.0000 0.0000007 100.0000 0.0000008 0.000000 2.0000009 0.000000 0.2000000E-01分析：得到了全局最优解，结果与matlab在（1）问中得到的局部最优解是一致的。（2）A的需求量变成600时候，仿照上例中lingo的代码，将100改成600，仍然可以得到全局最优解，如下所示：Variable Value Reduced CostX1 300.0000 0.000000X2 0.000000 2.000000X3 300.0000 0.000000X4 0.000000

10、6.000000X5 0.000000 0.000000X6 0.000000 0.000000Row Slack or Surplus Dual Price1 -600.0000 -1.0000002 0.000000 400.00003 0.000000 1000.0004 200.0000 0.0000005 500.0000 0.0000006 200.0000 0.0000007 0.000000 1.0000008 200.0000 0.0000009 0.000000 0.1333333E-01分析：应该用300t甲和300t丙去生产A，最大利润为60万，这与matlab得到的局

11、部最优解一样，验证可知该解在可行域内。在寻找全局最优解时候用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;运

12、行后得到的全局最优解如下：Variable Value Reduced CostX1 0.000000 1.000000X2 0.000000 0.000000X3 0.000000 0.000000X4 50.00000 0.000000X5 150.0000 0.000000X6 0.000000 0.5000000Row Slack or Surplus Dual Price1 -750.0000 -1.0000002 0.000000 200.00003 0.000000 350.00004 450.0000 0.0000005 350.0000 0.0000006 500.0000

13、0.0000007 100.0000 0.0000008 0.000000 3.750000 9 0.000000 0.2000000E-01分析：这与matlab的局部最优解一样，但是总的来说，用lingo求解更加方便，因为lingo往往可以直接给出问题的全局最优解且编程实现更直观。9.8投资组合问题一、模型建立1.本题为投资组合模型的具体应用，要找出投资方案，首先要得到A、B、C三种股票的期望收益率Er（平均值可以作为它的无偏估计），然后是三种股票的方差和协方差，对于股票i,j他们的方差及协方差在概率统计中都有相应的无偏估计形式，即以上i，j为股票类别，k为各个时期，n为总时期数。特别的，

14、当i=j时候，得到的结果为方差计算方法。在lingo中求期望收益率和协方差矩阵时候，需要通过以上公式变成计算，但是在matlab中直接提供了cov命令可以得到几个向量的协方差矩阵，大大简化了运算。2.本体的目标即为在给定收益率情况下，让风险最小。以原题要求为例即st:其中，而 ，cov为向量的协方差矩阵。（二）模型求解1.matlab实现（1）在matlab中编程如下所示：采用fmincon函数求解时：1）目标函数：function f=gupiaofun1(x,cov)f=x*cov*x;end2）命令文件：clc,clear alla1=1.300 1.103 1.216 0.954 0.

15、929 1.056 1.038 1.089 1.090 1.083 1.035 1.176;a2=a1-ones(1,12);b1=1.225 1.290 1.216 0.728 1.144 1.107 1.321 1.305 1.195 1.390 0.928 1.715;b2=b1-ones(1,12);c1=1.149 1.260 1.419 0.922 1.169 0.965 1.133 1.732 1.021 1.131 1.006 1.908;c2=c1-ones(1,12);cov=cov(a2,b2,c2)%将a1,b1,c1转化为列向量，用cov命令求协方差矩阵Er=mean

16、(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.5301 0.3564 0.1135fv = 0.0224ef = 1分析：结果是收敛的，且三种股票的持股比例分别为

17、53.01%，35.64%，11.35%，如果采用quadprog命令求解，则命令文件如下所示：clc,clear alla1=1.300 1.103 1.216 0.954 0.929 1.056 1.038 1.089 1.090 1.083 1.035 1.176;a2=a1-ones(1,12);b1=1.225 1.290 1.216 0.728 1.144 1.107 1.321 1.305 1.195 1.390 0.928 1.715;b2=b1-ones(1,12);c1=1.149 1.260 1.419 0.922 1.169 0.965 1.133 1.732 1.02

18、1 1.131 1.006 1.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)得到结果为：

19、x = 0.5301 0.3564 0.1135fv = 0.0224ef = 1分析：得到结果与用fmincon函数的结果一样，但是一般来说针对二次规划的问题，用quadprog函数往往会更好一点。问题一：如果年收益率在10%到100%之间变化，投资组合和风险如何变化；解答：针对不同的收益率，直接将原来的约束的右边换成相应的约束即可，但是考虑到Er=（0.0891，0.2137，0.2346），即最高收益率为23.46%，所以更高的收益率期望已不可能实现，所以以下只分析10%到24%区间的情况，变成如下所示：%期望收益率在10%到24%之间浮动的时候，去0.2个百分点为1个间隔分析clc,c

20、lear all a1=1.300 1.103 1.216 0.954 0.929 1.056 1.038 1.089 1.090 1.083 1.035 1.176;a2=a1-ones(1,12);b1=1.225 1.290 1.216 0.728 1.144 1.107 1.321 1.305 1.195 1.390 0.928 1.715;b2=b1-ones(1,12);c1=1.149 1.260 1.419 0.922 1.169 0.965 1.133 1.732 1.021 1.131 1.006 1.908;c2=c1-ones(1,12);cov=cov(a2,b2,c

21、2); 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);for i=0.1:0.002:0.24k=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;

22、endTplot(S(:,1),S(:,2)得到的结果为：T = 0.1000 0.9145 0.0725 0.0129 0.1020 0.8992 0.0839 0.0170 0.1040 0.8838 0.0952 0.0210 0.1060 0.8684 0.1066 0.0250 0.1080 0.8530 0.1179 0.0290 0.1100 0.8377 0.1293 0.0330 0.1120 0.8223 0.1407 0.0371 0.1140 0.8069 0.1520 0.0411 0.1160 0.7915 0.1634 0.0451 0.1180 0.7761 0

23、.1747 0.0491 0.1200 0.7608 0.1861 0.0532 0.1220 0.7454 0.1974 0.0572 0.1240 0.7300 0.2088 0.0612 0.1260 0.7146 0.2201 0.0652 0.1280 0.6993 0.2315 0.0692 0.1300 0.6839 0.2429 0.0733 0.1320 0.6685 0.2542 0.0773 0.1340 0.6531 0.2656 0.0813 0.1360 0.6377 0.2769 0.0853 0.1380 0.6224 0.2883 0.0894 0.1400

24、0.6070 0.2996 0.0934 0.1420 0.5916 0.3110 0.0974 0.1440 0.5762 0.3223 0.1014 0.1460 0.5608 0.3337 0.1055 0.1480 0.5455 0.3451 0.1095 0.1500 0.5301 0.3564 0.1135 0.1520 0.5147 0.3678 0.1175 0.1540 0.4993 0.3791 0.1215 0.1560 0.4840 0.3905 0.1256 0.1580 0.4686 0.4018 0.1296 0.1600 0.4532 0.4132 0.1336

25、 0.1620 0.4378 0.4245 0.1376 0.1640 0.4224 0.4359 0.1417 0.1660 0.4071 0.4473 0.1457 0.1680 0.3917 0.4586 0.1497 0.1700 0.3763 0.4700 0.1537 0.1720 0.3609 0.4813 0.1578 0.1740 0.3456 0.4927 0.1618 0.1760 0.3302 0.5040 0.1658 0.1780 0.3148 0.5154 0.1698 0.1800 0.2994 0.5267 0.1738 0.1820 0.2840 0.538

26、1 0.1779 0.1840 0.2687 0.5494 0.1819 0.1860 0.2533 0.5608 0.1859 0.1880 0.2379 0.5722 0.1899 0.1900 0.2225 0.5835 0.1940 0.1920 0.2072 0.5949 0.1980 0.1940 0.1918 0.6062 0.2020 0.1960 0.1764 0.6176 0.2060 0.1980 0.1610 0.6289 0.2100 0.2000 0.1456 0.6403 0.2141 0.2020 0.1303 0.6516 0.2181 0.2040 0.11

27、49 0.6630 0.2221 0.2060 0.0995 0.6744 0.2261 0.2080 0.0841 0.6857 0.2302 0.2100 0.0687 0.6971 0.2342 0.2120 0.0534 0.7084 0.2382 0.2140 0.0380 0.7198 0.2422 0.2160 0.0226 0.7311 0.2463 0.2180 0.0072 0.7425 0.2503 0.2200 0 0.6972 0.3028 0.2220 0 0.6016 0.3984 0.2240 0.0000 0.5060 0.4940 0.2260 0 0.41

28、04 0.5896 0.2280 0 0.3147 0.6853 0.2300 0 0.2191 0.7809 0.2320 0 0.1235 0.8765 0.2340 0 0.0279 0.9721 0.2360 -0.0030 -0.0000 1.0030 0.2380 -0.0072 -0.0000 1.0072 0.2400 -0.0114 -0.0000 1.0114分析：上面显示了在不同期望收益率情况下需要进行的投资比例的变化，注意到21.80%到22.00%之间发生了投资组合种类的变化，所以将该区间细分求出更精确的T矩阵如下如下所示：T = 0.2180 0.0072 0.74

29、25 0.2503 0.2181 0.0065 0.7431 0.2505 0.2182 0.0057 0.7436 0.2507 0.2183 0.0049 0.7442 0.2509 0.2184 0.0042 0.7448 0.2511 0.2185 0.0034 0.7453 0.2513 0.2186 0.0026 0.7459 0.2515 0.2187 0.0019 0.7465 0.2517 0.2188 0.0011 0.7470 0.2519 0.2189 0.0003 0.7476 0.2521 0.2190 0 0.7450 0.2550 0.2191 0 0.7402

30、 0.2598 0.2192 0 0.7355 0.2645 0.2193 0 0.7307 0.2693 0.2194 0 0.7259 0.2741 0.2195 0 0.7211 0.2789 0.2196 0 0.7163 0.2837 0.2197 0 0.7116 0.2884 0.2198 0 0.7068 0.2932 0.2199 0 0.7020 0.2980 0.2200 0 0.6972 0.3028分析：所以当期望收益率达到21.90%时候就应该只投资B、C两种股票，这是因为B、C两种股票的收益率更高的原因，为了得到更大收益率，必须选择风险更高的B、C两种股票。可以画出最小风险随着期望