数学实验第5次作业-无约束优化

实用

无约束优化

-实验目的

1掌握MATLAB优化工具箱的基本用法，对不同算法进行初步分析、比较。

2练习用无约束优化方法建立和求解实际问题的模型（包括最小二乘拟合）o

二实验内容

1取不同的初值计算下列非线性规划，尽可能求出所有局部极小点，进而找出全局极小点,

并对不同算法（搜索方向、搜索步长、数值梯度与分析梯度等）的结果进行分析、比较。

minz=(无逐2)2(1——%!—%2(1—%i)S]

理论计算部分：

根据多元微积分中求解函数的（局部）极值点的方法，首先要计算原函数的梯度以及

Hessiari矩阵。对于原函数

22

Z=(%1%2)(1—-一%2(1—%1)平，X6JR

首先计算梯度

dxx

dz

dx2.

22

丁=2X2(XI%2)(1-X1)[1-Xi-%2(1-41)S]-2(1-X1)(%1X2)[1-Xi-X2(l-*1)5]

OX\

422

+2(5X2(1-Xj)-1)(X1X2)(1-%1),

dz7

2522

记=2X1(X1X2)(1-%1)(1-%1-%2(1-%1)S]-2(1-X1)(X1X2)(1-Xi),

然后计算Hessian矩阵

d2zd2z

标dxrdx2

d2zd2z

距

此处由于计算所得式子非常复杂，所以不再继续计算。理论上可以通过求解非线性方程组

PZ（X）=0确定驻点，然后将每个驻点分别代入Kz中判断矩阵是否正定，即可以确定该驻

点是否为局部最优解。但是，通过上面的计算可以得知，要求解的非线性方程组较为复杂，

及时解除解，要再带入Hessian矩阵判断正定也比较困难，所以此处不进行解析求解。

文档

实用

观察原式，可以简单的看出，原式一定大于0,所以与=0或1,电=0是全局极小点,

通过将最后一个式子因式分解

52242

[1-Xj-X2(l-X1)]=(1-XX)[1-X2(l-^l)].

还可以看出，(2,1)也是一个局部极小点。

用MATLAB解决问题：

对于目标函数是2维以下的情况,不妨先绘制出可能出现局部最优解的区域的函数图象,

这样有利于直观的判断函数在这个区域的大致情况。

在MATLAB中输入以下命令：

»[x,y]=meshgrid(-1:0.01:2,-1:0.01:2);

»z=(x.*y)."2.*(l-x).-2.*(l-x-y.*5).-2;

»mesh(x,y,z>xlabel('),ylabelCy*),zlabelCz*);

»figure

»contour(x,y,z,50)

画出原图像如下:

画出等高线如下:

文档

实用

从图像上可以看出，由于所选择的范围太大，导致图像中右边大面积的点的值因为太小，所

以相对于较大的值来讲全部被认为是0,使得无法确定零点。

缩小选取的范围，输入命令如下：

»[x,y]=meshgrid(-0.1:0.01:1.1,-0.1:0.01:1.1);

»z=(x.*y).-2.*(l-x).-2.*(l-x-y.-5).*2;

»mesh(x,y,z),xlabelCx*),ylabelCy*),zlabelCz);

»figure;

»contour(x,y,z,300)

画出图像如下:

画出等高线如下:

文档

实用

对于点(2,1),输入命令如下：

»[x,y]=meshgrid(l.9:0.001:2.h0.9:0.001:1.1);

»z=(x.*y).*2.♦(1-x).*2.*(l-x-y.-5).*2;

»mesh(x,y,z),xlabelCx*),ylabel('y')，zlabel('z');

»figure;

»contour(x,y,z,300)

画出图像如下:

画出等高线如下:

文档

f

文档

实用

X=

-0.9203-0.0000

9.8419e-009

exit=

1

out=

iterations:13

funcCoxmt:51

stepsize:1

firstorderopt:6.7403e-004

algorithm:medium-scale:Quasi-Nevtonlinesearch,

message:[1x471char]

然后改变初值进行计算，输出结果如下:

初值最优解函数值迭代次数函数调用次数

[-1.-1](0.0000,-0.5000)016

[-2,-2](-0.0003,6.2061e-0072272

-1.2880)

[0.1,0.2](0.0000,0.1719)1.1184e-018624

[3,2](1.8782,1.6809)8.2624e-011939

[2,0.9](2.0248,0.9066)1.1459e-014530

从上表中可以得到以下结论，极小值为0,但是极小点不唯一。当初值选择的不同时，最优

解会向%=0或1或犯=0或者点(2,1)收敛。但是初值不同，收敛的速度不同，初值越

接近最优解，收敛越快。

文档

实用

下面用不同的搜索方向和步长搜索进行计算

在命令栏中输入以下命令（用三种搜索方向（BFGS,DFP和最速下降法）以及两种步长搜索

（混合二、三次插值和三次插值））：

»Xcomparingdifferentalgorithas:withoutusinggradientvector

formatshorte

x0=[0.I,0.2];

optl=optiMsetCLargeScale,,*offzMaxFunEvals*,1000);%casel:拟牛顿法BFGS公式，昆台二次和三次多项式插值

[xl,vl,exit1,out1]=fminunc(®work1,xO,opt1);

fopt=optinset(opt1,*HessUpdate*,*dfp*);%case2:拟牛顿法DFP公式，混合二次和三次多顼式插值

[x2,v2,exit2,out2]=fninunc(ork1,xO,fopt);

fopt=optlaset(opt1,'HessUpdate*,'steepdesc*);%case3:最速下降法，混合二次和三次多项式插值

[x3,v3,exit3,out3]=f»inunc(dtrorkl,xO,fopt);

opt2=optiBset(opt1,'LineSearchType?,'cubicpoly,);%caseL拟牛顿法BFGS公式，三次多项式插值

[x4,v4,exit4,out4]=fninunc(dvork1,xO,opt2);

fopt=optlaset(opt2,,HessUpdate,Jdfp');%case5:拟牛顿法DFP公式，混合二次和三次多项式插值

[x5,v5,exit5,out5]=fMinunc(dworkl,xO,fopt);

fopt=optiaset(opt2,,HessUpdate*,*steepdesc,);%case6:最速下降法，混合二次和三次多项式插值

[x6,v6,exit6,out6]=f»inunc(dworkl,x0,fopt);

»solutions=[x1;x2;x3;x4;x5;x6];%输出结果

funvalues=[vl;v2;v3;v4;v5;v6].

iterations=[out1.iterations;out2.iterations;out3.iterations;out4.iterations,out5.iterations;out6.iterations];

funcCount=[out1.funcCount;out2.funcCount;out3.funcCount;out4.funcCount;out5.funcCount;out6.funcCount],

[solutions,funvalues,iterations,funcCount]

分析梯度只需将命令中的

optl=optimset^LargeScale\off','MaxFunEvals'^OOO)

改为

optl=optimset^LargeScale','off','MaxFunEvals\lOOO'GradObj\'on')

即可。同时函数文件需要作如下改动(此处式子过长不适宜截图)：

function[yzg]=workl(x)

y=(x(l)*x(2))A2*(l-x(l))A2*(l-x(l)-x(2)*(l-x(l))A5)A2;

ifnargout>l

g(l)=2*x(2)*(x(l)*x(2))*(l-x(l))A2*(1-x(1)-x(2)*(1-x(1))A5)A2-2*(1-x

(1))*(x(l)*x(2))A2*(l-x(l)-x(2)*(l-x(l))A5)A2+2*(5*x(2)*(l-x(l))A4-l

)*(x(l)*x(2))A2*(l-x(l))A2;

g(2)=2*x(l)*(x(l)*x(2))*(l-x(l))A2*(1-x(1)-x(2)*(l-x(l))A5)A2-2*(1-x

(1))A5*(x(l)*x(2))A2*(l-x(l))A2;

end

文档

实用

然后在命令栏中输入以下命令:

»KcoMparingdifferentalgorithms:withoutusinggradientvector

foraatshorte

x0=[0.1,0.2];

optl=optinsetCLargeScale*,'off^,'MaxFunEvals,,lOOO/GradObj'Jon');Kcasel:捌牛顿法BFGS公式，混合二次和三次多项式插值

[xl,vl,exit1,outl]=faxnunc<dworkl,xO,opt1);

fopt=optinset(opt1,*HessUpdate*,*dfp/);!kase2:旗牛籁法DFP公式，混台二次和三次多项式插值

[x2,v2,exit2,out2]=fBinunc(®»rorkl,xO,fopt);

fopt=opti»set(opt1,'HessUpdate'Jsteepdesc*);%case3：最龙下法，混合二次和三次多项式插值

[x3,v3,exit3,out3]=f»munc(9vorkl,xO,fopt);

opt2=opti»set(opt1,'LmeSearchType,/cubicpolyJ);\casZ:椒牛顿法BFGS公式，三次多项式插值

[x4,v4,exit4,out4]=fMinunc(®workl,xO,opt2);

fopt=optunset(opt2,,HessUpdate*/dfp');肌ase5:拟牛顿法DFP公式，混合二次和三次多顼式插值

[x5,v5,exit5,out5]=fainxmc(®workl,xO,fopt);

fopt=optiMset(opt2,,HessUpdate,,*steepdesc,);务caseb:最漫下降法，温含二次和三次多项式临值

[x6,v6,exit6,out6]=fBinuncOworkbxO,fopt);

»solutions=[xl;x2;x3;x4;x5;x6];%输出结果

funvalues=[vl；v2,v3,v4;v5,v6];

iterations=[out1.iterations;out2.iterations;out3.iterations,out4.iterations.out5.iterations,out6.iterations],

funcCount=[out1.funcCount,out2.fxmcCount,out3.funcCount;out4.funcCount,out5.funcCount,out6.funcCoxmt];

[solutions,funvalues,iterations,funcCount]

输出结果如下表所示:

数值梯度

搜索方向步长搜索最优解X最优解y函数值迭代次数调用次数

BFGS-3.7957e-0071.7187e-0012.9188e-015521

混合二、

DFP-1.9022e-0071.7109e-0017.2776e-016521

三次插值

最速下降1.6709e-0051.6971e-0015.5436e-01218111

BFGS-3.7957e-0071.7187e-0012.9188e-015521

DFP三次插值-1.9022e-0071.7109e-0017.2776e-016521

最速下降1.6709e-0051.6971e-0015.5436e-01218111

分析梯度

搜索方向步长搜索最优解X最优解y函数值迭代次数调用次数

BFGS-1.4776e-0051.7567e-0014.5788e-01257

混合二、

DFP-5.7471e-0061.7434e-0016.8437e-01357

三次插值

最速下降1.6238e-0051.7418e-0015.4550e-0122245

BFGS-1.4776e-0051.7567e-0014.5788e-01257

DFP三次插值-5.7471e-0061.7434e-0016.8437e-01357

最速下降1.6238e-0051.7418e-0015.4550e-0122245

文档

实用

分析上表，可得到以下结果：

1函数的极小值为0,分别在打=0或1或必=0或者点(2,1)处取得。

2当初值不同时，对于得到的最优解、迭代的次数以及函数的调用次数都会有影响。而且

当初值不同时，得到的解不同，且得到的都是局部最优解。

3从表中可以看出，使用拟牛顿法的BFGS公式或者。FP公式时，需要的迭代次数相差不

多，但是使用最速下降法是需要的迭代系数相对就会较多，函数调用次数也会相应增多。

4分析梯度与数值梯度相比，迭代次数相差不大，相差较多的是函数调用次数。而且在本

题中，从数值来看，数值梯度相较于分析梯度要较好一些。

下面用自己实现的最速下降法和牛顿法来求解

最速下降法

在命令栏中输入以下内容：

X1(1)=0.5;

x2(l)=0.5;

i=l：

Ewhile1

a=xl(i);

b=x2(i);

p=2.*a.*(b.'2).*((l-a).,2).♦(1-a-b.*((l-a).*5).*2)..

-2.♦(1-a).*((a.*b).*2).*(1-a-b.*((l-a).*5)."2)..

+2.*(a.*b).'2.*(1-a-b.*(l-a)."5).♦(5.*b.*(a.-4)-1);

q=2.*b.*(a.*2).*((l-a)."2).*((1-a-b.*((l-a).*5)).*2)..

-2.*((a.*b)."2).*((l-a)."2).*(1-a-b.*((l~a).*5))...,

.*((l-a).*5);

i=i+l;

xl(i)=a-p;

x2(i)=b-q;

if((abs(xl(i)-xl(i-1))<=le-6)44(abs(x2(i)-x2(i-l))<=le-6))

break

end

-end

其中前面的p与q是手动求出的梯度，xl(l)=0.5,久2(1)=0.5为初值。。

输出迭代的结果，%!=0.9614,x2=0,2428,迭代次数为16166次。

文档

实用

2有一组数据&,%)(i=l,2,3……33),其中右=10(i-l),%由下表给出。现要求用

这组数据拟合函数

x-Xst

/(x,t)=Xj+x2e~^+x3e

中的参数x,初值可选为(0.5,1.5,-1,0,0.01,0.02),用GN和LM两种方法求解。对

力作一扰动，即%+e“勺为(-0.05,0.05)内的随机数，观察并分析迭代收敛是否会变慢。

i%i%i%

10.844120.718230.478

20.908130.685240.467

30.932140.658250.457

40.936150.628260.448

50.925160.603270.438

60.908170.580280.431

70.881180.558290.424

80.850190.538300.420

90.818200.522310.414

100.784210.506320.411

110.751220.490330.406

初步解决：

首先编制函数的M文件:

functiony=vork2(x,t,c)

y=x(l)+x(2)*exp(-x(4)*t)+x(3)*exp(-x(5)*t)-c;

然后在命令栏中输入以下命令：

»x0=[0.5,1.5,-1,0.Ob0.02];

»i=[l:l:33];

»t=10*(i-l);

»c=[0,8440.9080.9320.9360.9250.9080.8810.8500.8180.7840.751...

0.7180.6850.6580.6280.6030.5800.5580.5380.5220.5060.490...

0.4780.4670.4570.4480.4380.4310.4240.4200.4140.4110.406];

将所有参数输入到MATLAB中。

先用LM法计算，输入以下命令:

文档

实用

»optl=optimsetCLargeScale*,*off*,*MaxFunEvals*,1000);

»[xl,norml,resl,exit1,outl]=lsqnonlinCwork2*,xO,[],[],opt1,t,c)

此时，先不使用大规模算法，输出结果如下:

X1=

0.37541.9358-1.46470.01290.0221

norml=

5.4649e-005

resl=

Columns1through9

0.0026-0.0045-0.00100.00110.00290.00020.00050.0006-0.0006

Columns10through18

-0.0006-0.0013-0.00100.0007-0.00170.00080.0003-0.00010.0005

Columns19through27

0.0011-0.0006-0.00050.00120.0003-0.0003-0.0006-0.00080.0010

Columns28through33

0.00070.0012-0.00060.0003-0.0013-0.0003

exit1=

3

out1=

iterations:35

funcCount:230

stepsize:2.9676bo05

cgiterations:[]

firstorderopt:4.9405e-008

algorithm:'Levenberg-Marquardt

message:[1x111char]

然后再在命令栏中输入以下命令:

文档

实用

»optl=optinsetCLargeScale*,*on',*MaxFxmEvals*,1000);

»[xl,norml,resl,exit1,outl]=lsqnonlin(，work2,,xO,[],[],opt1,t,c)

此时使用大规模算法，输出结果如下：

X1=

0.37541.9358-1.46470.01290.0221

nornl=

5.4649e-005

resl=

Columns1through9

0.0026-0.0045-0.00100.00110.00290.00020.00050.0006-0.0006

Colunns10through18

-0.0006-0.0013-0.00100.0007-0.00170.00080.0003-0.00010.0005

Coluinns19through27

0.0011-0.0006-0.00050.00120.0003-0.0003-0.0006-0.00080.0010

Columns28through33

0.00070.0012-0.00060.0003-0.0013-0.0003

exit1=

1

out1=

firstorderopt:2.2762e-008

iterations:8

funcCount:54

cgiterations:0

algorithm:,large-scale:trust-regionreflectiveNewton'

message:[1x137char]

从两个输出结果可以看出，当使用了大规模算法之后，迭代次数和函数调用次数要明显少于

没有使用大规模算法的时候。但是二者的解一样。

再用GN法计算，输入以下命令，同样的，先不使用大规模算法：

文档

实用

»optl=optimsetCLargeScale,/off,JMaxFunEvals',1000);

»opt2=optimset(opt1/LevenbergMarquardt，,?off');

»[x2,norm2,res2,exit2,out2]=lsqnonlinCvork2,,xO,[],[],opt2,t,c)

输出结果如下:

x2=

0.37541.9358-1.46470.01290.0221

nor»2=

5.4649e-005

res2=

Columns1through9

0.0026-0.0045-0.00100.00110.00290.00020.00050.0006-0.0006

Columns10through18

-0.0006-0.0013-0.00100.0007-0.00170.00080.0003-0.00010.0005

Columns19through27

0.0011-0.0006-0.00050.00120.0003-0.0003-0.0006-0.00080.0010

Columns28through33

0.00070.0012-0.00060.0003-0.0013-0.0003

exit2=

out2=

iterations:9

funcCount:81

stepsize:0.9999

cgiterations:[]

firstorderopt:[]

algorithm:,medium-scale:Gauss-Nerton,line-search,

message:[1x147char]

然后再在命令栏中输入以下命令：

文档

实用

»optl=optimset('LargeScale','on,,sMaxFunEvals*,1000);

»opt2=optiinset(opt1,'LevenbergMarquardt1,'off');

»[x2,norm2,res2,exit2,out2]=lsqnonlin(,work2,,xO,[],[],opt2,t,c)

此时使用大规模算法，输出结果如下：

x2=

0.37541.9358-1.46470.01290.0221

nom2=

5.4649bo05

res2=

Columns1through9

0.0026-0.0045-0.00100.00110.00290.00020.00050.0006-0.0006

Colunns10through18

-0.0006-0.0013-0.00100.0007-0.00170.00080.0003-0.00010.0005

Columns19through27

0.0011-0.0006-0.00050.00120.0003-0.0003-0.0006-0.00080.0010

Colunns28through33

0.00070.0012-0.00060.0003-0.0013-0.0003

exit2=

1

out2=

firstorderopt:2.2762e-008

iterations:8

funcCount:54

cgiterations:0

algorithm:,large-scale:trust-regionreflectiveNevton/

message:[1x137char]

从输出结果可以看出，使用了大规模算法之后，迭代次数以及函数调用次数确实减少了很多。

但是最终算出的结果是一样的。最终算出的结果如下表所示：

文档

实用

0.375411.93585一1.464690.012870.02212

迭代次数的比较：

不使用大规模算法使用大规模算法

LM法358

GN法98

可以看出，LM法的收敛速度不如GN法快。

下面研究当％发生扰动之后对于迭代收敛速度的影响。(以下全部不使用大规模算法)

首先是当e»的所有元素均是相等的时候。

e,=-0.4,在命令栏中输入以下内容：

法：

»cl=c-0.04;

»optl=optinsetCLargeScale*,*off*,*MaxFunEvals*,1000);

»[xl,norml,resl,exit1,outl]=lsqnonlinCwork2*,xO,[],[],opt1,t,cl)

输出结果如下：

xl=

0.33541.9358-1.46470.01290.0221

out1=

iterations:35

funcCount:230

stepsize:2.7338e-005

cgiterations:[]

firstorderopt:4.2791e-008

algorithm:JLevenberg-Marquardt,

message:[1x111char]

GN法:

文档

实用

»cl=c-O.04;

»optl=optimsetCLargeScale"Joff*,'MaxFunEvals,,1000);

»opt2=optimset(opt1,LevenbergMarquardt,,off');

»[x2,norm2,res2,exit2,out2]=lsqnonlinCwork2*,xO,[],[],opt2,t,cl)

输出结果如下：

x2=

0.33541.9358-1.46470.01290.0221

out2=

iterations:9

funcCount:81

stepsize:1.0000

cgiterations:[]

firstorderopt:[]

algorithm:medium-scale:Gauss-Newton,line-search'

message:[1x147char]

=-0,1,在命令栏中输入以下内容：

法：

»c2=c-0.01;

»optl=optimset('LargeScale'Joff',MaxFunEvals,,1000);

»[xl,norml,resl,exit1,outl]=lsqnonlinCvork2*,xO,[],[],opt1,t,c2)

输出结果如下:

文档

实用

xl=

0.36541.9358-1.46470.01290.0221

out1=

iterations:35

funcCount:230

stepsize:2.5547e-005

cgiterations:[]

firstorderopt:3.8092e-008

algorithm:,Levenberg-Marquardt,

message:[lxll1char]

GN法:

»c2=c-0.01;

»optl=optimsetCLargeScale,,'off'MaxFunEvals',1000);

»opt2=optimset(opt1,'LevenbergMarquardtJJoff');

»[x2,norm2,res2,exit2,out2]=lsqnonlinCwork2J,xO,[],[],opt2,t,c2)

输出结果如下：

x2=

0.36541.9358-1.46470.01290.0221

out2=

iterations:9

funcCount:81

stepsize:0.9998

cgiterations:[]

firstorderopt:[]

algorithm:medium-scale:Gauss-Newton,line-search'

message:[1x147char]

=0.1,在命令栏中输入以下内容:

LM法，输出结果如下：

文档

实用

0.38541.9358-1.46470.01290.0221

out1=

iterations:35

funcCount:230

stepsize:3.0194e-005

cgiterations:[]

firstorderopt:5.0934e-008

algorithm:*Levenberg-Marquardt?

message:[1x111char]

GN法，输出结果如下:

0.38541.9358-1.46470.01290.0221

out2=

iterations:9

funcCount:81

stepsize:0.9999

cgiterations:[]

firstorderopt:[]

algorithm:*medium-scale:Gauss-Nevton,line-search*

message:[1x147char]

6=0.4,在命令栏中输入以下内容:

法，输出结果如下：

xl=

0.41541.9358-1.46470.01290.0221

文档

实用

out1=

iterations:35

funcCount:230

stepsize:2.4177e-005

cgiterations:[]

firstorderopt:3.4802e-008

algorithm:Levenberg-Marquardt1

message:[1x111char]

GN法,输出结果如下：

x2=

0.41541.9358-1.46470.01290.0221

out2=

iterations:9

funcCount:81

stepsize:0.9994

cgiterations:[]

firstorderopt:[]

algorithm:'medium-scale:Gauss-Newton,line-search'

message:[1x147char]

通过以上结果可以得到下面的表格：

结果的比较

尤1尤2九3无4工5

-法0.33541.9358-1.46470.01290.0221

0.04GN法0.33541.9358-1.46470.01290.0221

-法0.36541.9358-1.46470.01290.0221

0.01GN法0.36541.9358-1.46470.01290.0221

LM法0.38541.9358-1.46470.01290.0221

0.01

GN法0.38541.9358-1.46470.01290.0221

0.04LM法0.41541.9358-1.46470.01290.0221

文档

实用

GN法0.41541.9358-1.46470.01290.0221

迭代速度的比较

迭代次数

通过以上的比较可以看出，当力有一定的小扰动时㈤W(-0.05,0.05)),首先发生变

化的是Xi,从表中可以看出，至始至终在变化的只有勺，而其他的解都没有发生变化。其次,

迭代次数与函数调用次数没有发生变化，无论扰动有多大，迭代次数和函数调用次数都还是

保持原来的值不变，而且依然是LM法的迭代次数要大于GN法的迭代次数。

然后比较的元素全部是随机数的情况

LM法在命令栏中输入以下内容：

»e=-0.05+0.l*rand(l,33);

»cl=c+e;

J

»optl=optimset('LargeScale*,off*9*MaxFunEvals*91000);

»[xl,norml,resl,exit1,outl]=lsqnonlin('work2,,xO,[],[],opt1,t,cl)

输出结果如下：

文档

实用

X1=

0.30331.0629-0.48490.00780.0287

norml=

0.0288

resl=

Columns1through9

0.0058-0.02610.0448-0.0380-0.01090.03950.0232-0.0010-0.0401

Colunns10through18

-0.03740.0459-0.0321-0.02650.0203-0.00730.05900.0309-0.0182

Columns19through27

-0.0061-0.02570.00310.0647-0.0198-0.0314-0.0090-0.0200-0.0198

Coluxins28through33

0.0119-0.01710.0264-0.02930.03350.0069

exit1=

3

out1=

iterations:32

funcCount:212

stepsize:2.2798e-004

cgiterations:[]

firstorderopt:2.5547e-005

algorithm:'Levenberg-MarquardtJ

message:[1x111char]

GN法在命令栏中输入以下内容：

»opt2=optimset(opt1,'LevenbergMarquardt7,off');

»[x2,nonn2,res2,exit2,out2]=lsqnonlinCwork2,,xO,[],[],opt2,t,cl)

输出结果如下：

文档

实用

x2=

0.30331.0629-0.4850