实验7无约束优化

1、实验10 无约束优化 实验目的 1 1  掌握用matlab优化工具箱的基本用法，对不同算法进行初步分析、比较。2 2  练习用无约束优化方法建立和求解实际问题模型（包括非线性最小二乘拟合）。实验内容第2题：取不同的初值计算下列非线性规划，尽可能求出所有局部极小点，进而找出全局极小点，并对不同算法（搜索方向、步长搜索、数值梯度与分析梯度等）的结果进行分析、比较。（4），。其中，。问题分析：首先用数学方法计算出真实值， 将，代入，化z的表达式为：可见分母的值越小，z越小。当x1，x2的值在2.54之间时可能取最小值。为获得直观认识，我们再画出该函数的三维图

2、形和等高线图。x1,x2=meshgrid(2.5:0.1:4.5,2.5:0.1:4.5);c1=0.7;c2=0.73;y=-1./(x1-4).2+(x2-4).2+c1)-1./(x1-2.5).2+(x2-3.8).2+c2);mesh(x1,x2,y);pause; %画三维网格图contour(x1,x2,y,20) %画等高线图（20条）可以看出，z存在两个最小值。 对z求导：      matlab 程序设计及结果：下面用不同的算法进行计算：建立fun2.m文件：function y=fun2(x,c1,c2)y=-1./(x

3、(1)-4).2+(x(2)-4).2+c1)-1./(x(1)-2.5).2+(x(2)-3.8).2+c2);用数值方法计算函数的梯度：% comparing different algorithms: without using gradient vectorc1=0.7;c2=0.73;format short ex0=-1.9,2;'-case1: bfgs, hybrid 2,3 poly-'fopt=optimset('largescale','off', 'maxfunevals',1000);x1,v1,exit

4、1,out1=fminunc(fun2,x0,fopt,c1,c2)pause'-case2: dfp, hybrid 2,3 poly-'fopt=optimset('tolfun',1e-8,'tolx',1e-8,'largescale','off','hessupdate','dfp', 'maxfunevals',1000);x2,v2,exit2,out2=fminunc(fun2,x0,fopt,c1,c2)pause'-case3: steep

5、, hybrid 2,3 poly-'fopt=optimset('tolfun',1e-8,'tolx',1e-8,'largescale','off','hessupdate','steepdesc', 'maxfunevals',1000);x3,v3,exit3,out3=fminunc(fun2,x0,fopt,c1,c2)pause'-case4: bfgs, 3rd poly-'fopt=optimset('tolfun',1e-

6、8,'tolx',1e-8,'largescale','off','linesearchtype','cubicpoly', 'maxfunevals',1000);x4,v4,exit4,out4=fminunc(fun2,x0,fopt,c1,c2)pause'-case5: dfp, 3rd poly-'fopt=optimset('tolfun',1e-8,'tolx',1e-8,'largescale','off

7、9;,'hessupdate','dfp','linesearchtype','cubicpoly', 'maxfunevals',1000);x5,v5,exit5,out5=fminunc(fun2,x0,fopt,c1,c2)pause'-case6: steep, 3rd poly-'fopt=optimset('tolfun',1e-8,'tolx',1e-8,'largescale','off','hessupdat

8、e','steepdesc','linesearchtype','cubicpoly', 'maxfunevals',1000);x6,v6,exit6,out6=fminunc(fun2,x0,fopt,c1,c2) pause '+ results of solutions +'solutions=x1;x2;x3;x4;x5;x6''+ results of funvalues +'funvalues=v1,v2,v3,v4,v5,v6'+ funccount +

9、9;iterations=out1.funccount,out2.funccount,out3.funccount,out4.funccount,out5.funccount,out6.funccount得到结果为：搜索方向步长搜索最优解最优值迭代次数bfgs混合二、三次插值3.9060e+000 3.9875e+000-1.7752e+00054dfp混合二、三次插值3.9060e+000 3.9875e+000-1.7752e+00084最速下降混合二、三次插值3.9060e+000 3.9875e+000-1.7752e+000111bfgs三次插值3.9060e+000 3.9875e

10、+000-1.7752e+00057dfp三次插值3.9060e+000 3.9875e+000-1.7752e+00084最速下降三次插值3.9060e+000 3.9875e+000-1.7752e+000111用分析方法计算函数的梯度：建立grad.m文件给出fun2函数的梯度： function f,g=grad(x,c1,c2)f=-1./(x(1)-4).2+(x(2)-4).2+c1)-1./(x(1)-2.5).2+(x(2)-3.8).2+c2);% compute the function value at xif nargout > 1 % fun cal

11、led with 2 output arguments g(1)=2.*(x(1)-4)./(x(1)-4).2+(x(2)-4).2+c1).2)+2.*(x(1)-2.5)./(x(1)-2.5).2+(x(2)-3.8).2+c2).2); % compute the gradient evaluated at x g(2)=2.*(x(2)-4)./(x(1)-4).2+(x(2)-4).2+c1).2)+2.*(x(2)-3.8)./(x(1)-2.5).2+(x(2)-3.8).2+c2).2); endcomparing different algorithms: using g

12、radient vectorformat short ex0=-1.9,2;c1=0.7;c2=0.73;'-case1: bfgs, hybrid 2,3 poly-'fopt=optimset('largescale','off', 'maxfunevals',1000,'gradobj','on' );x1,v1,exit1,out1=fminunc(grad,x0,fopt,c1,c2)pause'-case2: dfp, hybrid 2,3 poly-'fopt=opti

13、mset('largescale','off','hessupdate','dfp', 'maxfunevals',1000,'gradobj','on' );x2,v2,exit2,out2=fminunc(grad,x0,fopt,c1,c2)pause'-case3: steep, hybrid 2,3 poly-'fopt=optimset('largescale','off','hessupdate','

14、;steepdesc', 'maxfunevals',1000,'gradobj','on' );x3,v3,exit3,out3=fminunc(grad,x0,fopt,c1,c2)pause'-case4: bfgs, 3rd poly-'fopt=optimset('largescale','off','linesearchtype','cubicpoly', 'maxfunevals',1000,'gradobj',&

15、#39;on' );x4,v4,exit4,out4=fminunc(grad,x0,fopt,c1,c2)pause'-case5: dfp, 3rd poly-'fopt=optimset('largescale','off','hessupdate','dfp','linesearchtype','cubicpoly', 'maxfunevals',1000,'gradobj','on' );x5,v5,exit5,out

16、5=fminunc(grad,x0,fopt,c1,c2)pause'-case6: steep, 3rd poly-'fopt=optimset('largescale','off','hessupdate','steepdesc','linesearchtype','cubicpoly', 'maxfunevals',1000,'gradobj','on' );x6,v6,exit6,out6=fminunc(grad,x0,fop

17、t,c1,c2)pause '+ results of solutions +'solutions=x1;x2;x3;x4;x5;x6'funvalues=v1,v2,v3,v4,v5,v6iterations=out1.funccount,out2.funccount,out3.funccount,out4.funccount,out5.funccount,out6.funct得到结果为：搜索方向步长搜索最优解最优值迭代次数bfgs混合二、三次插值3.9060e+000 3.9875e+000-1.7752e+0009dfp混合二、三次插值3.9060e+000 3.

18、9875e+000-1.7752e+0009最速下降混合二、三次插值3.9060e+000 3.9875e+000-1.7222e+00012bfgs三次插值3.9060e+000 3.9875e+000-1.7222e+0009dfp三次插值3.9060e+000 3.9875e+000-1.7222e+0009最速下降三次插值3.9060e+000 3.9875e+000-1.7222e+00012改变x0的值，比如（3，3）分别运行以上两程序，采用数值梯度算法，得到结果：搜索方向步长搜索最优解最优值迭代次数bfgs混合二、三次插值2.6071e+000 3.8143e+000-1.722

19、2e+00027dfp混合二、三次插值2.6071e+000 3.8143e+000-1.7222e+00027最速下降混合二、三次插值2.6071e+000 3.8143e+000-1.7222e+00069bfgs三次插值2.6071e+000 3.8143e+000-1.7222e+00030dfp三次插值2.6071e+000 3.8143e+000-1.7222e+00027最速下降三次插值2.6071e+000 3.8143e+000-1.7222e+00069 采用分析梯度算法，得到结果：搜索方向步长搜索最优解最优值迭代次数bfgs混合二、三次插值2.6071e+000

20、 3.8143e+000-1.7222e+0009dfp混合二、三次插值2.6071e+000 3.8143e+000-1.7222e+0009最速下降混合二、三次插值2.6071e+000 3.8143e+000-1.7222e+00012bfgs三次插值2.6071e+000 3.8143e+000-1.7222e+0009dfp三次插值2.6071e+000 3.8143e+000-1.7222e+0009最速下降三次插值2.6071e+000 3.8143e+000-1.7222e+00012      结果分析：从上面表格中的数据可以得知对于函

21、数的求解无论用数值法还是分析法，无论用什么搜索方向，什么步长搜索都可以得到最优解和最优值。但是从上面数据明显可以看出用分析法求解的迭代次数要小于数值法的迭代次数。对于函数z最速下降法的迭代次数最多，这是因为最速下降法使用负梯度作为搜索方向，只利用到一阶导数项，收敛阶数只为1，因此尽管开始时收敛较快，但是接近最优解时，收敛就变慢了。而bfgs和dfp算法均利用到了二阶导数项，收敛阶数较高，因此收敛速率较快。   第5题：某分子由25个原子组成，并且已经通过实验测量得到了其中某些原子队之间的距离（假设在平面结构上讨论），如表7.8所示。请你确定每个原子的位置关系。表7

22、.8原子对距离原子对距离原子对距离原子对距离（4，1）0.9607(5,4)0.4758(18,8)0.8363(15,13)0.5725（12,1）0.4399(12,4)1.3402(13,9)0.3208(19,13)0.7660(13,1)0.8143(24,4)0.7006(15,9)0.1574(15,14)0.4394(17,1)1.3765(8,6)0.4945(22,9)1.2736(16,14)1.0952(21,1)1.2722(13,6)1.0559(11,10)0.5781(20,16)1.0422(5,2)0.5294(19,6)0.6810(13,10)0.925

23、4(23,16)1.8255(16,2)0.6144(25,6)0.3587(19,10)0.6401(18,17)1.4325(17,2)0.3766(8,7)0.3351(20,10)0.2467(19,17)1.0851(25,2)0.6893(14,7)0.2878(22,10)0.4727(20,19)0.4995(5,3)0.9488(16,7)1.1346(18,11)1.3840(23,19)1.2277(20,3)0.8000(20,7)0.3870(25,11)0.4366(24,19)1.1271(21,3)1.1090(21,7)0.7511(15,12)1.0307(

24、23,21)0.7060(24,3)1.1432(14,8)0.4439(17,12)1.3904(23,22)0.8025       问题分析及模型假设： 据题意，不妨设第i点坐标为，其中第一个点坐标为，问题为求达到最小值的解，问题转化为无约束优化：，其中：      matlab 程序设计及结果：使用优化工具箱中的sqnonlin函数程序如下： function f=distance(x,d);y(1)=(x(2*3-1)2+(x(2*3)2-d(1)2;y(2)=(x(2*11-1)2+(x(2*

25、11)2-d(2)2;y(3)=(x(2*12-1)2+(x(2*12)2-d(3)2;y(4)=(x(2*16-1)2+(x(2*16)2-d(4)2;y(5)=(x(2*20-1)2+(x(2*20)2-d(5)2;y(6)=(x(2*5-3)-x(2*2-3)2+(x(2*5-2)-x(2*2-2)2-d(6)2;y(7)=(x(2*16-3)-x(2*2-3)2+(x(2*16-2)-x(2*2-2)2-d(7)2;y(8)=(x(2*17-3)-x(2*2-3)2+(x(2*17-2)-x(2*2-2)2-d(8)2;y(9)=(x(2*25-3)-x(2*2-3)2+(x(2*25

26、-2)-x(2*2-2)2-d(9)2;y(10)=(x(2*5-3)-x(2*3-3)2+(x(2*5-2)-x(2*3-2)2-d(10)2;y(11)=(x(2*20-3)-x(2*3-3)2+(x(2*20-2)-x(2*3-2)2-d(11)2;y(12)=(x(2*21-3)-x(2*3-3)2+(x(2*21-2)-x(2*3-2)2-d(12)2;y(13)=(x(2*24-3)-x(2*3-3)2+(x(2*24-2)-x(2*3-2)2-d(13)2;y(14)=(x(2*5-3)-x(2*4-3)2+(x(2*5-2)-x(2*4-2)2-d(14)2;y(15)=(x(

27、2*12-3)-x(2*4-3)2+(x(2*12-2)-x(2*4-2)2-d(15)2;y(16)=(x(2*24-3)-x(2*4-3)2+(x(2*24-2)-x(2*4-2)2-d(16)2;y(17)=(x(2*8-3)-x(2*6-3)2+(x(2*8-2)-x(2*6-2)2-d(17)2;y(18)=(x(2*13-3)-x(2*6-3)2+(x(2*13-2)-x(2*6-2)2-d(18)2;y(19)=(x(2*19-3)-x(2*6-3)2+(x(2*19-2)-x(2*6-2)2-d(19)2;y(20)=(x(2*25-3)-x(2*6-3)2+(x(2*25-2

28、)-x(2*6-2)2-d(20)2;y(21)=(x(2*8-3)-x(2*7-3)2+(x(2*8-2)-x(2*7-2)2-d(21)2;y(22)=(x(2*14-3)-x(2*7-3)2+(x(2*14-2)-x(2*7-2)2-d(22)2;y(23)=(x(2*16-3)-x(2*7-3)2+(x(2*16-2)-x(2*7-2)2-d(23)2;y(24)=(x(2*20-3)-x(2*7-3)2+(x(2*20-2)-x(2*7-2)2-d(24)2;y(25)=(x(2*21-3)-x(2*7-3)2+(x(2*21-2)-x(2*7-2)2-d(25)2;y(26)=(x

29、(2*14-3)-x(2*8-3)2+(x(2*14-2)-x(2*8-2)2-d(26)2;y(27)=(x(2*18-3)-x(2*8-3)2+(x(2*18-2)-x(2*8-2)2-d(27)2;y(28)=(x(2*13-3)-x(2*9-3)2+(x(2*13-2)-x(2*9-2)2-d(28)2;y(29)=(x(2*15-3)-x(2*9-3)2+(x(2*15-2)-x(2*9-2)2-d(29)2;y(30)=(x(2*22-3)-x(2*9-3)2+(x(2*22-2)-x(2*9-2)2-d(30)2;y(31)=(x(2*11-3)-x(2*10-3)2+(x(2*

30、11-2)-x(2*10-2)2-d(31)2;y(32)=(x(2*13-3)-x(2*10-3)2+(x(2*13-2)-x(2*10-2)2-d(32)2;y(33)=(x(2*19-3)-x(2*10-3)2+(x(2*19-2)-x(2*10-2)2-d(33)2;y(34)=(x(2*20-3)-x(2*10-3)2+(x(2*20-2)-x(2*10-2)2-d(34)2;y(35)=(x(2*22-3)-x(2*10-3)2+(x(2*22-2)-x(2*10-2)2-d(35)2;y(36)=(x(2*18-3)-x(2*11-3)2+(x(2*18-2)-x(2*11-2)

31、2-d(36)2;y(37)=(x(2*25-3)-x(2*11-3)2+(x(2*25-2)-x(2*11-2)2-d(37)2;y(38)=(x(2*15-3)-x(2*12-3)2+(x(2*15-2)-x(2*12-2)2-d(38)2;y(39)=(x(2*17-3)-x(2*12-3)2+(x(2*17-2)-x(2*12-2)2-d(39)2;y(40)=(x(2*15-3)-x(2*13-3)2+(x(2*15-2)-x(2*13-2)2-d(40)2;y(41)=(x(2*19-3)-x(2*13-3)2+(x(2*19-2)-x(2*13-2)2-d(41)2;y(42)=

32、(x(2*15-3)-x(2*14-3)2+(x(2*15-2)-x(2*14-2)2-d(42)2;y(43)=(x(2*16-3)-x(2*14-3)2+(x(2*16-2)-x(2*14-2)2-d(43)2;y(44)=(x(2*20-3)-x(2*16-3)2+(x(2*20-2)-x(2*16-2)2-d(44)2;y(45)=(x(2*23-3)-x(2*16-3)2+(x(2*23-2)-x(2*16-2)2-d(45)2;y(46)=(x(2*18-3)-x(2*17-3)2+(x(2*18-2)-x(2*17-2)2-d(46)2;y(47)=(x(2*19-3)-x(2*

33、17-3)2+(x(2*19-2)-x(2*17-2)2-d(47)2;y(48)=(x(2*20-3)-x(2*19-3)2+(x(2*20-2)-x(2*19-2)2-d(48)2;y(49)=(x(2*23-3)-x(2*19-3)2+(x(2*23-2)-x(2*19-2)2-d(49)2;y(50)=(x(2*24-3)-x(2*19-3)2+(x(2*24-2)-x(2*19-2)2-d(50)2;y(51)=(x(2*23-3)-x(2*21-3)2+(x(2*23-2)-x(2*21-2)2-d(51)2;y(52)=(x(2*23-3)-x(2*22-3)2+(x(2*23-

34、2)-x(2*22-2)2-d(52)2;接着运行：d=0.9607,0.4399,0.8143,1.3765,1.2722,0.5294,0.6144,0.3766,0.6893,0.9488,0.8000,1.1090,1.1432,0.4758,. 1.3402,0.7006,0.4945,1.0559,0.6810,0.3587,0.3351,0.2878,1.1346,0.3870,0.7511,0.4439,0.8363,0.3208,. 0.1574,1.2736,0.5781,0.9254,0.6401,0.2467,0.4727,1.3840,0.4366,1.0307,1.

35、3904,0.5725,0.7660,0.4394,. 1.0952,1.0422,1.8255,1.4325,1.0851,0.4995,1.2277,1.1271,0.7060,0.8052'x0=zeros(1,3),ones(1,45);x,norms=lsqnonlin(distance,x0,d);reshape(x,2,24)'结果如下： ans =  1.3355 0.0202 0.7213 -0.4567 0.6289 0.7704 1.0009 0.4489 0.6989 0.7247 0.4403 0.2500 0.2567 0.522

36、9 1.0616 -0.1718 0.7185 0.5979 0.6186 0.0171 0.1948 -0.4782 0.7421 -0.3334 0.6894 0.3821 1.0491 0.1071 1.2177 -0.5765 1.3140 0.3613 0.2867 1.3600 0.2281 0.2281 0.7393 0.3475 1.1384 0.5637 0.2603 0.8126 0.9425 1.2295 1.3042 0.5325 0.7727 0.4098 改变初值，会有不同的答案。  第8题：口服给药的方式相当于先有一个将药物从肠胃吸收入血液的过程，这个过程可简化为在药物进入中心室之前有一个吸收室（如图7.5所示），记中心室和吸收室的容积分别为v，v 1，而t时刻的血液浓度分别为；中心室的排除速率为k，吸收速率为（这里k和分别是中心室和吸收室血药浓度变化率与浓度本身的比例系数）。设时刻口服剂量为d的药物，容易写出吸收室的血药浓度的微分方程为中心室血药