《实验7 无约束优化.doc》由会员分享，可在线阅读，更多相关《实验7 无约束优化.doc（8页珍藏版）》请在三一文库上搜索。
1、实验10 无约束优化实验目的1 1 掌握用MATLAB优化工具箱的基本用法，对不同算法进行初步分析、比较。2 2 练习用无约束优化方法建立和求解实际问题模型（包括非线性最小二乘拟合）。实验内容第2题：取不同的初值计算下列非线性规划，尽可能求出所有局部极小点，进而找出全局极小点，并对不同算法（搜索方向、步长搜索、数值梯度与分析梯度等）的结果进行分析、比较。（4），。其中，。问题分析：首先用数学方法计算出真实值， 将，代入，化z的表达式为：可见分母的值越小，z越小。当x1，x2的值在2.54之间时可能取最小值。为获得直观认识，我们再画出该函数的三维图形和等高线图。x1,x2=meshgrid(2.
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(1)-4).2+(x(2)-4).2+c1)-1./(x(1)-2.5).2+(x(2)-3.8).2+c2
3、);用数值方法计算函数的梯度：% 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,exit1,out1=fminunc(fun2,x0,fopt,c1,c2)pause-case2: dfp, hybrid 2,3 poly-fopt=optimset(TolFun,1e-8,T
4、olX,1e-8,LargeScale,off,HessUpdate,dfp, MaxFunEvals,1000);x2,v2,exit2,out2=fminunc(fun2,x0,fopt,c1,c2)pause-case3: steep, 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
5、-fopt=optimset(TolFun,1e-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,HessUpdate,dfp,LineSearchType,cubicpoly, MaxFunEvals,1000);x5,v5,exit5,out5=fminun
6、c(fun2,x0,fopt,c1,c2)pause-case6: steep, 3rd poly-fopt=optimset(TolFun,1e-8,TolX,1e-8,LargeScale,off,HessUpdate,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 +funvalue
7、s=v1,v2,v3,v4,v5,v6+ funcCount +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+000
8、111BFGS三次插值3.9060e+000 3.9875e+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
9、nargout 1 % fun called 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
10、algorithms: using gradient 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=optimset(LargeScale,off,HessUpdate,dfp, MaxFunEvals,100
11、0,GradObj,on );x2,v2,exit2,out2=fminunc(grad,x0,fopt,c1,c2)pause-case3: steep, hybrid 2,3 poly-fopt=optimset(LargeScale,off,HessUpdate,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,cubicpol
12、y, MaxFunEvals,1000,GradObj,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,out5=fminunc(grad,x0,fopt,c1,c2)pause-case6: steep, 3rd poly-fopt=optimset(LargeScal
13、e,off,HessUpdate,steepdesc,LineSearchType,cubicpoly, MaxFunEvals,1000,GradObj,on );x6,v6,exit6,out6=fminunc(grad,x0,fopt,c1,c2)pause + results of solutions +solutions=x1;x2;x3;x4;x5;x6funvalues=v1,v2,v3,v4,v5,v6iterations=out1.funcCount,out2.funcCount,out3.funcCount,out4.funcCount,out5.funcCount,out
14、6.funct得到结果为：搜索方向步长搜索最优解最优值迭代次数BFGS混合二、三次插值3.9060e+000 3.9875e+000-1.7752e+0009DFP混合二、三次插值3.9060e+000 3.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.7222
15、e+00012改变x0的值，比如（3，3）分别运行以上两程序，采用数值梯度算法，得到结果：搜索方向步长搜索最优解最优值迭代次数BFGS混合二、三次插值2.6071e+000 3.8143e+000-1.7222e+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最速
16、下降三次插值2.6071e+000 3.8143e+000-1.7222e+00069采用分析梯度算法，得到结果：搜索方向步长搜索最优解最优值迭代次数BFGS混合二、三次插值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+00012BFGS三次插值2.6071e+000 3.8143e+000-1.7222e+0009DFP三次插值2.6071e+000 3.8143e+000-1.7222
17、e+0009最速下降三次插值2.6071e+000 3.8143e+000-1.7222e+00012 结果分析：从上面表格中的数据可以得知对于函数的求解无论用数值法还是分析法，无论用什么搜索方向，什么步长搜索都可以得到最优解和最优值。但是从上面数据明显可以看出用分析法求解的迭代次数要小于数值法的迭代次数。对于函数z最速下降法的迭代次数最多，这是因为最速下降法使用负梯度作为搜索方向，只利用到一阶导数项，收敛阶数只为1，因此尽管开始时收敛较快，但是接近最优解时，收敛就变慢了。而BFGS和DFP算法均利用到了二阶导数项，收敛阶数较高，因此收敛速率较快。第5题：某分子由25个原子组成，并且已经通过实
18、验测量得到了其中某些原子队之间的距离（假设在平面结构上讨论），如表7.8所示。请你确定每个原子的位置关系。表7.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
19、)0.5781(20,16)1.0422(5,2)0.5294(19,6)0.6810(13,10)0.9254(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.
20、4366(24,19)1.1271(21,3)1.1090(21,7)0.7511(15,12)1.0307(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*1
21、1)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-
22、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(2
23、*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)
24、-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(
25、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*1
26、1-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)2
27、-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)=(
28、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*1
29、7-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-2
30、)-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.3
31、904,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.5229 1.0616 -0.1718 0.718
32、5 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所示
33、），记中心室和吸收室的容积分别为V，V 1，而t时刻的血液浓度分别为；中心室的排除速率为k，吸收速率为（这里k和分别是中心室和吸收室血药浓度变化率与浓度本身的比例系数）。设时刻口服剂量为d的药物，容易写出吸收室的血药浓度的微分方程为中心室血药浓度的变化率由两部分组成：与c成正比的排除（比例系数k）；与成正比的吸收（比例系数）。再考虑到中心室和吸收室的容积分别为V，V 1，得到的微分方程为由以上两个微分方程不难解出中心室血药浓度。在制定给药方案时必须知道这种药物的3个参数，实际中通常通过实验数据确定。设时刻口服一定剂量的药物，表7.11是实验数据，请由此确定。表7.11t0.0830.1670.
34、250.500.751.01.510.921.127.336.435.538.434.8t2.253.04.06.08.010.012.024.223.615.78.28.32.21.8问题分析及模型假设：由及b=d/V得：Matlab 程序设计及结果：下面采用lsqnonlin命令，利用控制参数进行计算。输入程序：function f=problem8(x,t,c)f=x(1)*x(2)*(exp(-x(3)*t)-exp(-x(2)*t)/(x(2)-x(3)-c;x0=5,3,2;%初值%x0=8,3,1t=0.083 0.167 0.25 0.50 0.75 1.0 . 1.5 2.2
35、5 3.0 4.0 6.0 8.0 10.0 12.0;c=10.9 21.1 27.3 36.4 35.5 38.4 . 34.8 24.2 23.6 15.7 8.2 8.3 2.2 1.8;%第一个控制参数，设定不采用大规模算法，设定函数的最大调用次数opt1=optimset(Largescale,off,MaxfunEvals,800);%第二个控制参数，搜索方向使用LM法opt2=optimset(opt1,LevenbergMarquardt,on);x,norm,res=lsqnonlin(exe81fun,x0,opt2,t,c)改变初值，结果相同：x = 46.8275 3.6212 0.2803norm = 34.2317res = Columns 1 through 9 1.1090 -0.3885 -0.5052 -0.5820 2.2767 -1.4071 -1.6859 2.8020 -1.7058 Columns 10 through 14 0.8438 1.2452 -2.9075 0.8787 -0.0423所以，所求参数值分别为：b=46.8275， k1=3.6212， k=0.2803。收获与体会：本次实验通过实践基本掌握了非线性无约束优化和最小二乘优化问题的方法，在尝试和探索过程中，可以观察到，初值的设定对最终的优化结果的影响很大。