实验9_非线性规划_陈雨_2010012199.pdf
实验实验 9 非线性非线性规划规划 陈雨陈雨 基应基应 02 2010012199 【实验目的】掌握用 MATLAB 优化工具箱和 LINGO 解非线性规划的方法;练习建立实际问题的非线性规划模型。【实验内容】题目 4 4.1 数学模型 我们先约定如下记号:1:a 生产 A 所用甲原料 2:a 生产 A 所用乙原料 3:a 生产 A 所用丙原料 1:b 生产 B 所用甲原料 2:b 生产 B 所用乙原料 3:b 生产 B 所用丙原料:a A 产品产量:b B 产品产量 1:c 原料甲进货量 2:c 原料乙进货量 3:c 原料丙进货量 其中,1,2,3iiicab i 3311,iiiiaa bb 根据题目,可以列出如下约束条件与目标函数:原料甲乙先混合,则有:1 22 1a ba b 产品 A、B 的含硫量约束:1231233%1%2%2.5%3%1%2%1.5%aaaabbbb 原料甲、乙、丙的供应量约束:123,500c c c 产品 A、B 最大市场需求量约束:100,200ab 非负性约束:,0,1,2,3iia bi 目标函数即净利润:12391561610zabccc 4.2 模型优化 4.2.1 第一小问 我们直接利用 MATLAB 优化工具箱函数 fmincon,注意我们需要将规划问题化为函数输入格式,详细见函数代码 exp904.m:注:由于 word 格式限制,上述代码中几处换行在 MATLAB 执行时需做调整。function exp904 clear all;clc;x0=0;0;100;0;200;0;A1=(eye(3),zeros(3)+zeros(3),eye(3);b1=500*ones(3,1);v1=zeros(6,1);x,fv,ef,out=fmincon(objfun,x0,A1,b1,v1,confun)function f=objfun(x)f=-(9*(x(1)+x(2)+x(3)+15*(x(4)+x(5)+x(6)-6*(x(1)+x(4)-16*(x(2)+x(5)-10*(x(3)+x(6);end function c1,c2=confun(x)c1=0.005,-0.015,-0.005,0,0,0*x;0,0,0,0.015,-0.005,0.005*x;1,1,1,0,0,0*x-100;0,0,0,1,1,1*x-200;c2=x(1)*x(5)-x(2)*x(4);end end 注:此处123123,Txa a a b b b 根据 MATLAB 的部分输出结果,我们得到结论:公司最优进货量为,甲 0t,乙 100t,丙 100t;并全部生产 200t 产品 B;最大净利润为 400千元(40 万元)。x=-0.0000 0 0.0000 -0.0000 100.0000 100.0000 fv=-400 ef=1 out=iterations:8 funcCount:56 lssteplength:1 stepsize:1.0958e-19 algorithm:medium-scale:SQP,Quasi-Newton,line-search firstorderopt:4.4409e-16 constrviolation:4.3831e-20 message:1x783 char 4.2.2 第二小问 我们需要改变约束条件,另外由于规划对初值敏感,我们还需改变初值 x0,x0=300;0;300;0;0;0 程序如下:function exp90402 clear all;clc;x0=300;0;300;0;0;0;A1=(eye(3),zeros(3)+zeros(3),eye(3);b1=500*ones(3,1);v1=zeros(6,1);x,fv,ef,out=fmincon(objfun,x0,A1,b1,v1,confun)function f=objfun(x)f=-(9*(x(1)+x(2)+x(3)+15*(x(4)+x(5)+x(6)-.6*(x(1)+x(4)-.16*(x(2)+x(5)-10*(x(3)+x(6);end function c1,c2=confun(x)c1=0.005,-0.015,-0.005,0,0,0*x;0,0,0,.0.015,-0.005,0.005*x;.1,1,1,0,0,0*x-600;0,0,0,1,1,1*x-200;c2=x(1)*x(5)-x(2)*x(4);end end 此时公司的最优决策为:购入原料甲 300t,原料乙 0t,原料丙 300t;全部生产产品 A 600t;净利润为 600 千元(60 万元)。x=300 0 300 0 0 0 fv=-600 ef=1 out=iterations:1 funcCount:7 lssteplength:1 stepsize:0 algorithm:medium-scale:SQP,Quasi-Newton,line-search firstorderopt:8.8818e-16 constrviolation:0 message:1x783 char 4.2.3 第三小问 4.2.3.1 情况(1)function exp9040301 clear all;clc;x0=0;0;100;0;200;0;A1=(eye(3),zeros(3)+zeros(3),eye(3);b1=500*ones(3,1);v1=zeros(6,1);x,fv,ef,out=fmincon(objfun,x0,A1,b1,v1,confun)function f=objfun(x)f=-(9*(x(1)+x(2)+x(3)+15*(x(4)+x(5)+x(6)-.6*(x(1)+x(4)-.13*(x(2)+x(5)-10*(x(3)+x(6);end function c1,c2=confun(x)c1=0.005,-0.015,-0.005,0,0,0*x;0,0,0,.0.015,-0.005,0.005*x;.1,1,1,0,0,0*x-100;0,0,0,1,1,1*x-200;c2=x(1)*x(5)-x(2)*x(4);end end 公司的最优决策为,购入原料甲 50t,原料乙 150t,原料丙 0t;全部生产 B 产品 200t;最大净利润为 750 千元(75 万元)。x=0 0 0 50.0000 150.0000 -0.0000 fv=-750.0000 ef=1 out=iterations:4 funcCount:28 lssteplength:1 stepsize:4.2019e-14 algorithm:medium-scale:SQP,Quasi-Newton,line-search firstorderopt:2.1316e-13 constrviolation:5.6843e-14 message:1x783 char 4.2.3.2 情况(2)注:选取以上初值有另一组解,同样得到最大净利润。function exp9040302 clear all;clc;x0=300;0;300;0;0;0;A1=(eye(3),zeros(3)+zeros(3),eye(3);b1=500*ones(3,1);v1=zeros(6,1);x,fv,ef,out=fmincon(objfun,x0,A1,b1,v1,confun)function f=objfun(x)f=-(9*(x(1)+x(2)+x(3)+15*(x(4)+x(5)+x(6)-.6*(x(1)+x(4)-.13*(x(2)+x(5)-10*(x(3)+x(6);end function c1,c2=confun(x)c1=0.005,-0.015,-0.005,0,0,0*x;0,0,0,.0.015,-0.005,0.005*x;.1,1,1,0,0,0*x-600;0,0,0,1,1,1*x-200;c2=x(1)*x(5)-x(2)*x(4);end end function exp9040302 clear all;clc;x0=0;0;0;0;0;0;A1=(eye(3),zeros(3)+zeros(3),eye(3);b1=500*ones(3,1);v1=zeros(6,1);x,fv,ef,out=fmincon(objfun,x0,A1,b1,v1,confun)function f=objfun(x)f=-(9*(x(1)+x(2)+x(3)+15*(x(4)+x(5)+x(6)-.6*(x(1)+x(4)-.13*(x(2)+x(5)-10*(x(3)+x(6);end function c1,c2=confun(x)c1=0.005,-0.015,-0.005,0,0,0*x;0,0,0,.0.015,-0.005,0.005*x;.1,1,1,0,0,0*x-600;0,0,0,1,1,1*x-200;c2=x(1)*x(5)-x(2)*x(4);end end 公司的最优决策为,购入原料甲 450t,原料乙 150t,原料丙 0t;全部生产 A 产品 600t;最大净利润为 750 千元(75 万元)。x=450.0000 150.0000 0 -0.0000 0.0000 -0.0000 fv=-750.0000 ef=4 out=iterations:3 funcCount:21 lssteplength:1 stepsize:9.8012e-07 algorithm:medium-scale:SQP,Quasi-Newton,line-search firstorderopt:1.5921e-06 constrviolation:9.0039e-10 message:1x762 char 公司的最优决策为,购入原料甲 50t,原料乙 150t,原料丙 0t;全部生产 B 产品 200t;最大净利润为 750 千元(75 万元)。x=0.0000 0 0.0000 50.0000 150.0000 0.0000 fv=-750 ef=1 out=iterations:9 funcCount:63 lssteplength:1 stepsize:1.9013e-18 algorithm:medium-scale:SQP,Quasi-Newton,line-search firstorderopt:4.4409e-16 constrviolation:1.2705e-19 message:1x783 char 题目 8 8.1 数学模型 记号约定如下:,:a b c 分别投资于 A,B,C 三支股票的资金比率:r 三支股票的收益向量:E 投资组合期望收益:Var 投资组合收益方差(风险)我们可以列出约束条件和目标函数:1,0abca b c 123()()()15%EaE rbE rcE r 123222123121323()()()()2(,)2(,)2(,)zVarVar arbrcra Var rb Var rc Var rabCov r racCov r rbcCov r r 8.2 模型优化 8.2.1 期望年收益至少为 15%,我们只需最小化投资组合的风险。clear all;clc;A=1.300 1.225 1.149;1.103 1.290 1.260;1.216 1.216 1.419;0.954 0.728 0.922;0.929 1.144 1.169;1.056 1.107 0.965;1.038 1.321 1.133;1.089 1.305 1.732;1.090 1.195 1.021;1.083 1.390 1.131;1.035 0.928 1.006;1.176 1.715 1.908;Corr=corrcoef(A);r=sum(A)./12-1;Std=std(A);H=2*cov(A);A1=-r;b1=-0.15;A2=ones(1,3);b2=1;v1=zeros(3,1);x0=0,0,1;x,fv,ef,out=quadprog(H,A1,b1,A2,b2,v1,x0)三种股票的资金投入分别为:53.0%,35.6%,11.4%。8.2.2 第一小问 x=0.5301 0.3564 0.1135 fv=0.0224 ef=1 out=iterations:3 constrviolation:0 algorithm:active-set message:Optimization terminated.firstorderopt:1.3878e-17 cgiterations:注:红线代表 A,蓝线代表 B,绿线代表 C 值得注意的点在收益率为 21.9%左右,此时投资组合的结构发生变化。超过该点时,A 不再包含于投资组合中,仅仅持有 B 与 C。另外风险在该点后会加速上升。8.2.3 第二小问 增加一种投资方式,亦即增加一个决策变量。clear all;clc;A=1.300 1.225 1.149;1.103 1.290 1.260;1.216 1.216 1.419;0.954 0.728 0.922;0.929 1.144 1.169;1.056 1.107 0.965;1.038 1.321 1.133;1.089 1.305 1.732;1.090 1.195 1.021;1.083 1.390 1.131;1.035 0.928 1.006;1.176 1.715 1.908;r=sum(A)./12-1;H=2*cov(A);A1=-r;A2=ones(1,3);b2=1;v1=zeros(3,1);x0=0,0,1;for i=100:234 b1=-0.001*i;x(:,i),fv(:,i)=quadprog(H,A1,b1,A2,b2,v1,);re(i)=0.001*i;end fv=fv.0.5;subplot(1,2,1),plot(re,x(1,:),r,re,x(2,:),b,re,x(3,:),g);subplot(1,2,2),plot(re,fv),grid,xlabel(return),ylabel(risk);clear all;clc;A=1.300 1.225 1.149;1.103 1.290 1.260;1.216 1.216 1.419;0.954 0.728 0.922;0.929 1.144 1.169;1.056 1.107 0.965;1.038 1.321 1.133;1.089 1.305 1.732;1.090 1.195 1.021;1.083 1.390 1.131;1.035 0.928 1.006;1.176 1.715 1.908;r=sum(A)./12-1;r=r,0.05;H=2*cov(A);H=H,zeros(3,1);zeros(1,4);A1=-r;b1=-0.15;A2=ones(1,4);b2=1;v1=zeros(4,1);x0=0,0,1,0;x,fv,ef,out=quadprog(H,A1,b1,A2,b2,v1,x0)购买 A,B,C 与国库券的比率为:8.7%,42.9%,14.3%,34.1%。8.2.4 第三小问 考虑将 C 出售,另购入适量的 A 与 B。考虑交易费用,我们需对约束条件做些修改,详细见如下程序:x=0.0869 0.4285 0.1434 0.3412 fv=0.0208 ef=1 out=iterations:2 constrviolation:1.1102e-16 algorithm:active-set message:Optimization terminated.firstorderopt:2.3096e-17 cgiterations:clear all;clc;A=1.300 1.225 1.149;1.103 1.290 1.260;1.216 1.216 1.419;0.954 0.728 0.922;0.929 1.144 1.169;1.056 1.107 0.965;1.038 1.321 1.133;1.089 1.305 1.732;1.090 1.195 1.021;1.083 1.390 1.131;1.035 0.928 1.006;1.176 1.715 1.908;r=sum(A)./12-1;H=2*cov(A);A1=-(r(1)-0.01),-(r(2)-0.01),-(r(3)+0.01);b1=-0.15-0.15*0.01+0.01*0.85;A2=1.01,1.01,0.99;b2=1-0.01*0.15+0.01*0.85;v1=0.5,0.35,0;x0=0,1,0;x,fv,ef,out=quadprog(H,A1,b1,A2,b2,v1,x0)我们应该出售部分 C,购买 A,而 B 不变。最后的比率为:52.65%,35%,12.3%。x=0.5265 0.3500 0.1230 fv=0.0226 ef=1 out=iterations:2 constrviolation:0 algorithm:active-set message:Optimization terminated.firstorderopt:2.7756e-17 cgiterations: