实验9非线性规划

实验 9

【实验目的】

非线性规划

1. 了解非线性优化问题求解的基本思想，掌握用 MATLAB 优化工具箱解非线性规划的基本 方法； 2. 练习建立实际问题的非线性规划模型。 【实验内容】 【题目 3】 对问题

2 min{100(x2  x1 )2  (1  x1 )2  90(x 4  x3 )2  (1  x3 )2  10.1 1  x2 )2 [( 2

    (1  x 4 )2 ]  19.8 x 2  1)(x 4  1)}增加以下条件，并分别取初值(3, 1, 3, 1)和 (

(3, , , )，求解非线性规划： 131

（1）  10 

x i  10 ；

（2） 10  x i  10,x 1x 2  x 1  x 2  1.5  0,x 1x 2  10  0, 100  x 1x 2x 3x 4  100 ；  （3）  10  xi  10,x1x 2  x1  x 2  1.5  0,x1x 2  10  0,x1  x 2  0,x1x 2x3x 4  16 。 再试取不同的初值或用分析梯度计算，比较计算结果，你能从中得到什么启示？ 1.1 目标函数的 M 文件的编写

2 设 f(x )  100(x2  x1 )2  (1  x1 )2  90 x 4  x3 )2  (1  x3 )2  10.1 1  x2 )2 ( [( 2

 (1  x 4 )2 ]  19.8 x 2  1)(x 4  1)}。 (

现在需要 f(x)的梯度，下面利用 matlab 里的 diff 命令求解梯度，建立文件 tidu.m

symsx1x2x3x4; z=100*(x2-x1^2)^2+(1-x1)^2+90*(x4-x3^2)^2+(1-x3)^2+10.1*((1-x2)^2+(1x4)^2)+19.8*(x2-1)*(x4-1); zx1=diff(z,x1);zx2=diff(z,x2);zx3=diff(z,x3);zx4=diff(z,x4); [zx1; zx2; zx3; zx4 ]

输出结果如下 ans = 2*x1 - 400*x1*(x2 - x1^2) - 2 - 200*x1^2 + (1101*x2)/5 + (99*x4)/5 - 40 2*x3 - 360*x3*(x4 - x3^2) - 2 - 180*x3^2 + (99*x2)/5 + (1001*x4)/5 - 40

由于可能用到黑塞矩阵，下面用 diff 命令求函数 f(x)的黑塞矩阵，建立文件 hessian.m

zx1x1=diff(zx1,x1);zx1x2=diff(zx1,x2);zx1x3=diff(zx1,x3);zx1x4=diff(z x1,x4); zx2x2=diff(zx2,x2);zx2x3=diff(zx2,x3);zx2x4=diff(zx2,x4); zx3x3=diff(zx3,x3);zx3x4=diff(zx3,x4); zx4x4=diff(zx4,x4); [zx1x1,zx1x2,zx1x3,zx1x4;zx1x2,zx2x2,zx2x3,zx2x4;zx1x3,zx2x3,zx3x3,zx 3x4;zx1x4,zx2x4,zx3x4,zx4x4;]

输出结果如下： ans = [ 1200*x1^2 - 400*x2 + 2, -400*x1, 0, 0] [ -400*x1, 1101/5, 0, 99/5] [ 0, 0, 1080*x3^2 - 360*x4 + 2, -360*x3] [ 0, 99/5, -360*x3, 1001/5] 下面编写目标函数的 M 文件 fun3.m：

f=100*(x(2)-x(1)^2)^2+(1-x(1))^2+90*(x(4)-x(3)^2)^2+(1-x(3))^2+10.1*( (1-x(2))^2+(1-x(4))^2)+19.8*(x(2)-1)*(x(4)-1); ifnargout>1 g(1)=2*x(1) - 400*x(1)*(x(2) - x(1)^2) - 2; g(2)=- 200*x(1)^2 + (1101*x(2))/5 + (99*x(4))/5 - 40 g(3)=2*x(3) - 360*x(3)*(x(4) - x(3)^2) - 2 g(4)=- 180*x(3)^2 + (99*x(2))/5 + (1001*x(4))/5 - 40 end ifnargout>2 H=[ 1200*x(1)^2 - 400*x(2) + 2, -400*x(1), -400*(x)1, 0, 0, end end 99/5, 1101/5, 0, -360*x3, 0, 99/5 1001/5]; 0

0, 1080*x3^2 - 360*x4 + 2, -360*x3

1.2 第一小问的求解 由于第一小问的约束条件只有上下限约束：  10  束的 M 文件。 又由于第一小问的约束条件只有上下限约束， 除了一般的中等规模算法之外， 还可使用 大规模算法。 下面利用 matlab 的 fmincon 命令

求解，分别使用中等规模的数值、分析算法以及大规 模的分析梯度算法求解。建立文件 fmincon31.m。

x0=[-3,-1,-3,-1];

x i  10 ，故无需编写非线性约

%中等规模算法（数值算法） opt1=optimset('largescale','off','MaxIter',3000,'MaxFun',20000); [x11,fv11,ef11,out11]=fmincon(@fun3,x0,[],[],[],[],[-10,-10,-10,-10], [10,10,10,10],[],opt1); %中等规模算法（分析算法） opt1=optimset('largescale','off','MaxIter',3000,'MaxFun',20000); opt2=optimset(opt1,'GradObj','on','GradCon','on','DerivativeCheck','o n'); [x12,fv12,ef12,out12]=fmincon(@fun3,x0,[],[],[],[],[-10,-10,-10,-10], [10,10,10,10],[],opt2); %大规模算法（分析算法） opt1=optimset('largescale','on','MaxIter',3000,'MaxFun',20000); opt2=optimset(opt1,'GradObj','on','GradCon','on','Hessian','on','Deri vativeCheck','on'); [x13,fv13,ef13,out13]=fmincon(@fun3,x0,[],[],[],[],[-10,-10,-10,-10], [10,10,10,10],[],opt2); % results of solutions solutions=[x11;x12;x13]; funvalues=[fv11;fv12;fv13]; iterations=[out11.funcCount;out12.funcCount;out13.funcCount]; rr1=[solutions,funvalues,iterations]

按照上述以上程序，分别设置初值为[3,1,3,1]，[-2,-2,-2,-2]，[2,2,2,2]，并将结果整理得 到如下表格： 初值 [-3,-1, -3,-1] [3,1, 3,1] [-2,-2, -2,-2] [2,2, 2,2] 算法 中等规模 大规模 中等规模 大规模 中等规模 大规模 中等规模 大规模 梯度 数值 分析 分析 数值 分析 分析 数值 分析 分析 数值 分析 分析 最优解 x1 1.0002 1.0002 0.9984 1.0001 1.0002 1.0011 1.0000 1.0000 0.9997 1.0001 1.0001 1.0010 最优解 X2 1.0004 1.0004 0.9968 1.0002 1.0004 1.0023 1.0000 1.0000 0.9995 1.0002 1.0002 1.0020 最优解 X3 0.9998 0.9998 1.0017 1.0000 0.9998 0.9987 1.0000 1.0000 1.0002 0.9999 0.9999 0.9990 最优解 X4 0.9996 0.9996 1.0034 1.0000 0.9997 0.9975 1.0000 1.0000 1.0003 0.9998 0.9998 0.9979 最优值 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0 0 0.0000 0.0000 0.0000 0.0000 目标函数 调用次数 136.0000 67.0000 204.0000 368.0000 160.0000 37.0000 12.0000 9.0000 302.0000 102.0000 54.0000 49.0000

对上表中的数据进行分析，可以得到以下结论： （1）使用中等规模算法时，若给出分析梯度，但是目标函数的调用次数会明显减少。

故对本问题，给出分析梯度对结果有积极影响。 （2）对于本题，初值取[-3,-1,-3,-1]和[-2,-2,-2,-2]时，使用中等规模算法的迭代次数和目 标函数使用次数明显小于使用大规模算法的，但是初值取[3,1,3,1]和[2,2,2,2]时，使用中等规 模算法的迭代次数和目标函数使用次数大于使用大规模算法的。 （3）对于本题，使用大规模算法得到的最优解和最优值均不如使用大规模算法得到的 结果好。 （4）对于本题，初值的改变并不影响最终的结果。 （5） 对第一小问，最优解为 x1= x2= x3= x4=1，最优值为 f(x)=0。 1.3 第二小问的求解 由于第二小问的约束条件不是线

性的，故需要编写约束函数的 M 文件。 1.3.1 非线性约束函数的 M 文件的编写 第二小文的非线性约束条件有：

x1x 2  x1  x 2  1.5  0,

x1x 2  10  0,

 100  x 1x 2x 3x 4  100

上下界约束条件有：

 10  x i  10,

编写如下文件 con32.m

function [c,ceq,g,geq]=con32(x) c=[x(1)*x(2)-x(1)-x(2)+1.5;-10-x(1)*x(2);x(1)*x(2)*x(3)*x(4)-100;-100 -x(1)*x(2)*x(3)*x(4)]; ceq=0; ifnargout>2 g=[x(2)-1,-x(2),x(2)*x(3)*x(4),-x(2)*x(3)*x(4);x(1)-1,-x(1),x(1)*x(3) *x(4),-x(1)*x(3)*x(4);0,0,x(1)*x(2)*x(4),-x(1)*x(2)*x(4);0,0,x(1)*x(2 )*x(3),-x(1)*x(2)*x(3)]; geq=[0;0;0;0]; end end

1.3.2 用 fmincon 命令求解 编写如下程序 fmincon32.m

x0=[-3,-1,-3,-1]; %-------------数值梯度------------opt1=optimset('largescale','off','MaxIter',3000,'MaxFun',20000); [x1,fv1,ef1,out1]=fmincon(@fun3,x0,[],[],[],[],[-10,-10,-10,-10],[10,

10,10,10],@con32,opt1); [c11,c12]=con32(x); %--------------·分析梯度-------------opt1=optimset('largescale','off','MaxIter',3000,'MaxFun',20000); opt2=optimset(opt1,'GradObj','on','GradCon','on','DerivativeCheck','o n'); [x2,fv2,ef2,out2]=fmincon(@fun3,x0,[],[],[],[],[-10,-10,-10,-10],[10, 10,10,10],@con32,opt2); [c21,c22]=con32(x); %--------- results of solutions-------solutions=[x1;x2]; funvalues=[fv1;fv2]; iterations=[out1.funcCount;out2.funcCount]; rr1=[solutions,funvalues,iterations]

按照上述方法，将初值分别设为[3,1,3,1]，[-2,-2,-2,-2]，[2,2,2,2]，并将输出结果进行整 理，得到如下表格： 初值 [-3,-1, -3,-1] [3,1, 3,1] [-2,-2, -2,-2] [2,2, 2,2] 梯度 数值 分析 数值 分析 数值 分析 数值 分析 最优解 x1 最优解 X2 最优解 X3 最优解 X4 最优值 目标函数 调用次数

1.4200 1.4200 1.4206 1.4206 1.4200 1.4200 -1.1082 -1.1082

-0.1904 -0.1904 -0.1888 -0.1888 -0.1904 -0.1904 1.2372 1.2372

1.4661 1.4661 -1.4394 -1.4394 1.4660 1.4660 0.8803 0.8803

2.1511 2.1511 2.0813 2.0813 2.1511 2.1511 0.7743 0.7742

487.9879 487.9879 493.7994 493.7994 487.9877 487.9876 4.4898 4.4898

97 52 111 60 164 96 282 122

通过对数据的分析，可以得到以下结论： （1）使用中等规模算法（SQP）时，若给出分析梯度，目标函数的调用次数会明显减 少。故对本问题，给出分析梯度对结果有积极影响。 （2）无论采用是不是分析梯度，得到的最优解和最优值都是一样的。 （3）对于本题，选择的初值不同，得到的结果也不同。 1.4 第三小问的求解 1.4.1 非线性约束函数的 M 文件的编写 第二小问的非线性约束条件有：

x1x 2  x1  x 2  1.5  0,

x1x 2  10  0,

x1x 2x 3x 4  16

上下界约束条件有：

 10  x i  10,

线性约束有：

x 1  x 2  0,

编写如下文件con33.m

function [c,ceq,g,geq]=con33(x) c=[x(1)*x(2)-x(1)-x(2)+1.5;-10-x(1)*x(2)]; ceq=x(1)*x(2)*x(3)*x(4)-16; ifnargout>2 g=[x(2)-1,-x(2);x(1)-1,-x(1);0,0;0,0]; geq=[x(2)*x(3)*x(4);x(1)*x(3)*x(4);x(1)*x(2)*

82.8559 ef = 5 out = iterations: 199 funcCount: 3552 lssteplength: 1 stepsize: 1.1586e-005 algorithm: 'medium-scale: SQP, Quasi-Newton, line-search' firstorderopt: 13.5346 constrviolation: 3.3885e-020 message: [1x843 char] 将运行结果进行整理，得到如下表格：

工地(i) C12 C12

1 0 3

2 5 0

3 0 4

4 0 7

5 0 6

6 11 0

料场 (7.2500,7.2500) (3.3565, 5.6060)

总吨公里数为82.8559，比用原临时料场时的总吨公里数减小了53.3441。 2.4 程序改进 这个问题的等式约束简单，可以直接融入目标函数中去。 修改目标函数 M 文件如下。

function f=xuanzhifun1(x) f=x(5)*((x(1)-1.25)^2+(x(3)-1.25)^2)^0.5 ..., +x(6)*((x(1)-8.75)^2+(x(3)-0.75)^2)^0.5 ..., +x(7)*((x(1)-0.5)^2+(x(3)-4.75)^2)^0.5 ..., +x(8)*((x(1)-5.75)^2+(x(3)-5)^2)^0.5 ..., +x(9)*((x(1)-3)^2+(x(3)-6.5)^2)^0.5 ..., +x(10)*((x(1)-7.25)^2+(x(3)-7.25)^2)^0.5 ..., +(3-x(5))*((x(2)-1.25)^2+(x(4)-1.25)^2)^0.5 ..., +(5-x(6))*((x(2)-8.75)^2+(x(4)-0.75)^2)^0.5 ..., 21 +(4-x(7))*((x(2)-0.5)^2+(x(4)-4.75)^2)^0.5 ..., +(7-x(8))*((x(2)-5.75)^2+(x(4)-5)^2)^0.5 ..., +(6-x(9))*((x(2)-3)^2+(x(4)-6.5)^2)^0.5 ..., +(11-x(10))*((x(2)-7.25)^2+(x(4)-7.25)^2)^0.5; end

这样， 目标函数自变量的维数下降了不少， 可能对程序的运行和问题的解决有积极影响。 相应的，运行程序修改如下。

x0=[5,2,1,7,0,0,0,0,0,0]; A1=[0,0,0,0,1,1,1,1,1,1;0,0,0,0,-1,-1,-1,-1,-1,-1]; b1=[20,-16]; v1=[-inf,-inf, -inf, -inf,0,0,0,0,0,0]; v2=[inf, inf, inf, inf, 3,5,4,7,6,11]; opt1=optimset('largescale','off','MaxIter',1000,'MaxFun',20000); [x,fv,ef,out]=fmincon(@xuanzhifun1,x0,A1,b1,[],[],v1,v2,[],opt1)

输出结果如下： x= 3.2250 7.2500 5.6523 6.0000 0 fv= 82.8256 ef = 5 out = iterations: 35 funcCount: 439 7.2500 3.0000 0 4.0000 7.0000

lssteplength: 1 stepsize: 1.9703e-005 algorithm: 'medium-scale: SQP, Quasi-Newton, line-search' firstorderopt: 4.8991 constrviolation: 0 message: [1x843 char] 总吨公里数为82.8256， 比程序改进前算出来的结果82.8559还要略小一些， 比用原临时 料场时的总吨公里数减小了53.3744。 另外，程序改进之后，迭代次数和目标函数调用次数均下降了许多，说明这个改进对问 题的解决有积极影响。

实验 9

【实验目的】

非线性规划

1. 了解非线性优化问题求解的基本思想，掌握用 MATLAB 优化工具箱解非线性规划的基本 方法； 2. 练习建立实际问题的非线性规划模型。 【实验内容】 【题目 3】 对问题

2 min{100(x2  x1 )2  (1  x1 )2  90(x 4  x3 )2  (1  x3 )2  10.1 1  x2 )2 [( 2

    (1  x 4 )2 ]  19.8 x 2  1)(x 4  1)}增加以下条件，并分别取初值(3, 1, 3, 1)和 (

(3, , , )，求解非线性规划： 131

（1）  10 

x i  10 ；

（2） 10  x i  10,x 1x 2  x 1  x 2  1.5  0,x 1x 2  10  0, 100  x 1x 2x 3x 4  100 ；  （3）  10  xi  10,x1x 2  x1  x 2  1.5  0,x1x 2  10  0,x1  x 2  0,x1x 2x3x 4  16 。 再试取不同的初值或用分析梯度计算，比较计算结果，你能从中得到什么启示？ 1.1 目标函数的 M 文件的编写

2 设 f(x )  100(x2  x1 )2  (1  x1 )2  90 x 4  x3 )2  (1  x3 )2  10.1 1  x2 )2 ( [( 2

 (1  x 4 )2 ]  19.8 x 2  1)(x 4  1)}。 (

现在需要 f(x)的梯度，下面利用 matlab 里的 diff 命令求解梯度，建立文件 tidu.m

symsx1x2x3x4; z=100*(x2-x1^2)^2+(1-x1)^2+90*(x4-x3^2)^2+(1-x3)^2+10.1*((1-x2)^2+(1x4)^2)+19.8*(x2-1)*(x4-1); zx1=diff(z,x1);zx2=diff(z,x2);zx3=diff(z,x3);zx4=diff(z,x4); [zx1; zx2; zx3; zx4 ]

输出结果如下 ans = 2*x1 - 400*x1*(x2 - x1^2) - 2 - 200*x1^2 + (1101*x2)/5 + (99*x4)/5 - 40 2*x3 - 360*x3*(x4 - x3^2) - 2 - 180*x3^2 + (99*x2)/5 + (1001*x4)/5 - 40

由于可能用到黑塞矩阵，下面用 diff 命令求函数 f(x)的黑塞矩阵，建立文件 hessian.m

zx1x1=diff(zx1,x1);zx1x2=diff(zx1,x2);zx1x3=diff(zx1,x3);zx1x4=diff(z x1,x4); zx2x2=diff(zx2,x2);zx2x3=diff(zx2,x3);zx2x4=diff(zx2,x4); zx3x3=diff(zx3,x3);zx3x4=diff(zx3,x4); zx4x4=diff(zx4,x4); [zx1x1,zx1x2,zx1x3,zx1x4;zx1x2,zx2x2,zx2x3,zx2x4;zx1x3,zx2x3,zx3x3,zx 3x4;zx1x4,zx2x4,zx3x4,zx4x4;]

输出结果如下： ans = [ 1200*x1^2 - 400*x2 + 2, -400*x1, 0, 0] [ -400*x1, 1101/5, 0, 99/5] [ 0, 0, 1080*x3^2 - 360*x4 + 2, -360*x3] [ 0, 99/5, -360*x3, 1001/5] 下面编写目标函数的 M 文件 fun3.m：

f=100*(x(2)-x(1)^2)^2+(1-x(1))^2+90*(x(4)-x(3)^2)^2+(1-x(3))^2+10.1*( (1-x(2))^2+(1-x(4))^2)+19.8*(x(2)-1)*(x(4)-1); ifnargout>1 g(1)=2*x(1) - 400*x(1)*(x(2) - x(1)^2) - 2; g(2)=- 200*x(1)^2 + (1101*x(2))/5 + (99*x(4))/5 - 40 g(3)=2*x(3) - 360*x(3)*(x(4) - x(3)^2) - 2 g(4)=- 180*x(3)^2 + (99*x(2))/5 + (1001*x(4))/5 - 40 end ifnargout>2 H=[ 1200*x(1)^2 - 400*x(2) + 2, -400*x(1), -400*(x)1, 0, 0, end end 99/5, 1101/5, 0, -360*x3, 0, 99/5 1001/5]; 0

0, 1080*x3^2 - 360*x4 + 2, -360*x3

1.2 第一小问的求解 由于第一小问的约束条件只有上下限约束：  10  束的 M 文件。 又由于第一小问的约束条件只有上下限约束， 除了一般的中等规模算法之外， 还可使用 大规模算法。 下面利用 matlab 的 fmincon 命令

求解，分别使用中等规模的数值、分析算法以及大规 模的分析梯度算法求解。建立文件 fmincon31.m。

x0=[-3,-1,-3,-1];

x i  10 ，故无需编写非线性约

%中等规模算法（数值算法） opt1=optimset('largescale','off','MaxIter',3000,'MaxFun',20000); [x11,fv11,ef11,out11]=fmincon(@fun3,x0,[],[],[],[],[-10,-10,-10,-10], [10,10,10,10],[],opt1); %中等规模算法（分析算法） opt1=optimset('largescale','off','MaxIter',3000,'MaxFun',20000); opt2=optimset(opt1,'GradObj','on','GradCon','on','DerivativeCheck','o n'); [x12,fv12,ef12,out12]=fmincon(@fun3,x0,[],[],[],[],[-10,-10,-10,-10], [10,10,10,10],[],opt2); %大规模算法（分析算法） opt1=optimset('largescale','on','MaxIter',3000,'MaxFun',20000); opt2=optimset(opt1,'GradObj','on','GradCon','on','Hessian','on','Deri vativeCheck','on'); [x13,fv13,ef13,out13]=fmincon(@fun3,x0,[],[],[],[],[-10,-10,-10,-10], [10,10,10,10],[],opt2); % results of solutions solutions=[x11;x12;x13]; funvalues=[fv11;fv12;fv13]; iterations=[out11.funcCount;out12.funcCount;out13.funcCount]; rr1=[solutions,funvalues,iterations]

按照上述以上程序，分别设置初值为[3,1,3,1]，[-2,-2,-2,-2]，[2,2,2,2]，并将结果整理得 到如下表格： 初值 [-3,-1, -3,-1] [3,1, 3,1] [-2,-2, -2,-2] [2,2, 2,2] 算法 中等规模 大规模 中等规模 大规模 中等规模 大规模 中等规模 大规模 梯度 数值 分析 分析 数值 分析 分析 数值 分析 分析 数值 分析 分析 最优解 x1 1.0002 1.0002 0.9984 1.0001 1.0002 1.0011 1.0000 1.0000 0.9997 1.0001 1.0001 1.0010 最优解 X2 1.0004 1.0004 0.9968 1.0002 1.0004 1.0023 1.0000 1.0000 0.9995 1.0002 1.0002 1.0020 最优解 X3 0.9998 0.9998 1.0017 1.0000 0.9998 0.9987 1.0000 1.0000 1.0002 0.9999 0.9999 0.9990 最优解 X4 0.9996 0.9996 1.0034 1.0000 0.9997 0.9975 1.0000 1.0000 1.0003 0.9998 0.9998 0.9979 最优值 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0 0 0.0000 0.0000 0.0000 0.0000 目标函数 调用次数 136.0000 67.0000 204.0000 368.0000 160.0000 37.0000 12.0000 9.0000 302.0000 102.0000 54.0000 49.0000

对上表中的数据进行分析，可以得到以下结论： （1）使用中等规模算法时，若给出分析梯度，但是目标函数的调用次数会明显减少。

故对本问题，给出分析梯度对结果有积极影响。 （2）对于本题，初值取[-3,-1,-3,-1]和[-2,-2,-2,-2]时，使用中等规模算法的迭代次数和目 标函数使用次数明显小于使用大规模算法的，但是初值取[3,1,3,1]和[2,2,2,2]时，使用中等规 模算法的迭代次数和目标函数使用次数大于使用大规模算法的。 （3）对于本题，使用大规模算法得到的最优解和最优值均不如使用大规模算法得到的 结果好。 （4）对于本题，初值的改变并不影响最终的结果。 （5） 对第一小问，最优解为 x1= x2= x3= x4=1，最优值为 f(x)=0。 1.3 第二小问的求解 由于第二小问的约束条件不是线

性的，故需要编写约束函数的 M 文件。 1.3.1 非线性约束函数的 M 文件的编写 第二小文的非线性约束条件有：

x1x 2  x1  x 2  1.5  0,

x1x 2  10  0,

 100  x 1x 2x 3x 4  100

上下界约束条件有：

 10  x i  10,

编写如下文件 con32.m

function [c,ceq,g,geq]=con32(x) c=[x(1)*x(2)-x(1)-x(2)+1.5;-10-x(1)*x(2);x(1)*x(2)*x(3)*x(4)-100;-100 -x(1)*x(2)*x(3)*x(4)]; ceq=0; ifnargout>2 g=[x(2)-1,-x(2),x(2)*x(3)*x(4),-x(2)*x(3)*x(4);x(1)-1,-x(1),x(1)*x(3) *x(4),-x(1)*x(3)*x(4);0,0,x(1)*x(2)*x(4),-x(1)*x(2)*x(4);0,0,x(1)*x(2 )*x(3),-x(1)*x(2)*x(3)]; geq=[0;0;0;0]; end end

1.3.2 用 fmincon 命令求解 编写如下程序 fmincon32.m

x0=[-3,-1,-3,-1]; %-------------数值梯度------------opt1=optimset('largescale','off','MaxIter',3000,'MaxFun',20000); [x1,fv1,ef1,out1]=fmincon(@fun3,x0,[],[],[],[],[-10,-10,-10,-10],[10,

10,10,10],@con32,opt1); [c11,c12]=con32(x); %--------------·分析梯度-------------opt1=optimset('largescale','off','MaxIter',3000,'MaxFun',20000); opt2=optimset(opt1,'GradObj','on','GradCon','on','DerivativeCheck','o n'); [x2,fv2,ef2,out2]=fmincon(@fun3,x0,[],[],[],[],[-10,-10,-10,-10],[10, 10,10,10],@con32,opt2); [c21,c22]=con32(x); %--------- results of solutions-------solutions=[x1;x2]; funvalues=[fv1;fv2]; iterations=[out1.funcCount;out2.funcCount]; rr1=[solutions,funvalues,iterations]

按照上述方法，将初值分别设为[3,1,3,1]，[-2,-2,-2,-2]，[2,2,2,2]，并将输出结果进行整 理，得到如下表格： 初值 [-3,-1, -3,-1] [3,1, 3,1] [-2,-2, -2,-2] [2,2, 2,2] 梯度 数值 分析 数值 分析 数值 分析 数值 分析 最优解 x1 最优解 X2 最优解 X3 最优解 X4 最优值 目标函数 调用次数

1.4200 1.4200 1.4206 1.4206 1.4200 1.4200 -1.1082 -1.1082

-0.1904 -0.1904 -0.1888 -0.1888 -0.1904 -0.1904 1.2372 1.2372

1.4661 1.4661 -1.4394 -1.4394 1.4660 1.4660 0.8803 0.8803

2.1511 2.1511 2.0813 2.0813 2.1511 2.1511 0.7743 0.7742

487.9879 487.9879 493.7994 493.7994 487.9877 487.9876 4.4898 4.4898

97 52 111 60 164 96 282 122

通过对数据的分析，可以得到以下结论： （1）使用中等规模算法（SQP）时，若给出分析梯度，目标函数的调用次数会明显减 少。故对本问题，给出分析梯度对结果有积极影响。 （2）无论采用是不是分析梯度，得到的最优解和最优值都是一样的。 （3）对于本题，选择的初值不同，得到的结果也不同。 1.4 第三小问的求解 1.4.1 非线性约束函数的 M 文件的编写 第二小问的非线性约束条件有：

x1x 2  x1  x 2  1.5  0,

x1x 2  10  0,

x1x 2x 3x 4  16

上下界约束条件有：

 10  x i  10,

线性约束有：

x 1  x 2  0,

编写如下文件con33.m

function [c,ceq,g,geq]=con33(x) c=[x(1)*x(2)-x(1)-x(2)+1.5;-10-x(1)*x(2)]; ceq=x(1)*x(2)*x(3)*x(4)-16; ifnargout>2 g=[x(2)-1,-x(2);x(1)-1,-x(1);0,0;0,0]; geq=[x(2)*x(3)*x(4);x(1)*x(3)*x(4);x(1)*x(2)*

82.8559 ef = 5 out = iterations: 199 funcCount: 3552 lssteplength: 1 stepsize: 1.1586e-005 algorithm: 'medium-scale: SQP, Quasi-Newton, line-search' firstorderopt: 13.5346 constrviolation: 3.3885e-020 message: [1x843 char] 将运行结果进行整理，得到如下表格：

工地(i) C12 C12

1 0 3

2 5 0

3 0 4

4 0 7

5 0 6

6 11 0

料场 (7.2500,7.2500) (3.3565, 5.6060)

总吨公里数为82.8559，比用原临时料场时的总吨公里数减小了53.3441。 2.4 程序改进 这个问题的等式约束简单，可以直接融入目标函数中去。 修改目标函数 M 文件如下。

function f=xuanzhifun1(x) f=x(5)*((x(1)-1.25)^2+(x(3)-1.25)^2)^0.5 ..., +x(6)*((x(1)-8.75)^2+(x(3)-0.75)^2)^0.5 ..., +x(7)*((x(1)-0.5)^2+(x(3)-4.75)^2)^0.5 ..., +x(8)*((x(1)-5.75)^2+(x(3)-5)^2)^0.5 ..., +x(9)*((x(1)-3)^2+(x(3)-6.5)^2)^0.5 ..., +x(10)*((x(1)-7.25)^2+(x(3)-7.25)^2)^0.5 ..., +(3-x(5))*((x(2)-1.25)^2+(x(4)-1.25)^2)^0.5 ..., +(5-x(6))*((x(2)-8.75)^2+(x(4)-0.75)^2)^0.5 ..., 21 +(4-x(7))*((x(2)-0.5)^2+(x(4)-4.75)^2)^0.5 ..., +(7-x(8))*((x(2)-5.75)^2+(x(4)-5)^2)^0.5 ..., +(6-x(9))*((x(2)-3)^2+(x(4)-6.5)^2)^0.5 ..., +(11-x(10))*((x(2)-7.25)^2+(x(4)-7.25)^2)^0.5; end

这样， 目标函数自变量的维数下降了不少， 可能对程序的运行和问题的解决有积极影响。 相应的，运行程序修改如下。

x0=[5,2,1,7,0,0,0,0,0,0]; A1=[0,0,0,0,1,1,1,1,1,1;0,0,0,0,-1,-1,-1,-1,-1,-1]; b1=[20,-16]; v1=[-inf,-inf, -inf, -inf,0,0,0,0,0,0]; v2=[inf, inf, inf, inf, 3,5,4,7,6,11]; opt1=optimset('largescale','off','MaxIter',1000,'MaxFun',20000); [x,fv,ef,out]=fmincon(@xuanzhifun1,x0,A1,b1,[],[],v1,v2,[],opt1)

输出结果如下： x= 3.2250 7.2500 5.6523 6.0000 0 fv= 82.8256 ef = 5 out = iterations: 35 funcCount: 439 7.2500 3.0000 0 4.0000 7.0000

lssteplength: 1 stepsize: 1.9703e-005 algorithm: 'medium-scale: SQP, Quasi-Newton, line-search' firstorderopt: 4.8991 constrviolation: 0 message: [1x843 char] 总吨公里数为82.8256， 比程序改进前算出来的结果82.8559还要略小一些， 比用原临时 料场时的总吨公里数减小了53.3744。 另外，程序改进之后，迭代次数和目标函数调用次数均下降了许多，说明这个改进对问 题的解决有积极影响。