数学实验——无约束优化.docx

数学实验报告
实验6 无约束优化

实验6 无约束优化
分1 黄浩 2011011743
一、 实验目的
1. 掌握用MATLAB 优化工具箱的基本用法，对不同算法进行初步分析、比较。
2. 练习用无约束优化方法建立和求解实际问题模型（包括非线性最小二乘拟合）。
二、 实验内容
1. 《数学实验》第二版（问题2.1）
问题叙述：
取不同的初值计算非线性规划：
minz=x1x22(1-x1)2[1-x1-x2(1-x1)5]2
尽可能求出所有局部极小点，进而找出全局极小点，并对不同算法（搜索方向、步长搜索、数值梯度与分析梯度等）的结果进行分析、比较。
实验过程：
首先绘制这个函数的三维图形以及等高线（程序见四.1），结果如下：

通过观察这两幅图，可以得到，x2确定时，x1越负，函数值越大，x1确定时，x2绝对值越大，函数值越大。但对于x1正向偏离0的情况，并没有很好的反映，于是扩大绘图范围，做出下图（程序见四.2）：

由上面两幅图可见，方程像是一个四角被捏起的花布，而且z的最小值为0。因此只要求解该方程的零点，即得到了方程的局部极小点，且若将原方程变形为：
z=x1x22(1-x1)4[1-x2(1-x1)4]2
我们容易发现，该方程的零点为：x2=0或x1=0或x1=1或x2(1-x1)4=1
在求解零点之前，先针对一个零点，不妨用x1=1，x2=1，分析不同算法的优劣。
在matlab的无约束优化中，可以使用fminumc和fminsearch两种函数，搜索方向的算法有BFGS公式、DFP公式和最速下降法三种（书中还提到的Gill-Murray公式在matlab中已经不再使用），步长的一维搜索有混合二次三次多项式插值和三次多项式插值两种方法，另外，在求解函数梯度是也有数值方法和分析方法两种。在对上述四类算法因素进行分析时，我们采用控制变量法，每次只保持一种或两种算法因素改变，分析它的精度及效率。
（一） 分析fminumc与fminsearch两种方法的精度及效率
选择初值为x1=0.8,x2=0.8，使用fminunc和fminsearch的默认算法及控制参数，输出结果如下（程序见四.3、四.4）：

最优解x1
最优解x2
最优值z
迭代次数
z调用次数
fminunc
9.9400e-001
7.9238e-001
8.0435e-010
12
39
fminsearch
1.0000e+000
5.8842e-001
9.5246e-057
67
123
因为精确解为x1=1，z=0，我们便可以比较出不同算法的精度。
由上表可见，fminunc的迭代次数和目标函数的调用次数较少，效率较高，但最优值的精度较低，与之相反的是，fminsearch的迭代次数和目标函数的调用次数较多，效率比较低，然而精度很高。
（二） 比较不同搜索方向和不同步长的算法精度及效率
使用fminunc函数，选择初值为x1=0.6,x2=0.8，分别用BFGS公式、DFP公式和最速下降法三种搜索方向，混合二次三次多项式插值和三次多项式插值两种步长算法，共2*3=6种方式进行运算，结果如下（程序见四.5）：
搜索方向
步长算法
最优解x1
最优解x2
最优值z
迭代次数
z调用次数
BFGS
混合二次三次多项式插值
1.0072e+000
7.6715e-001
1.5634e-009
7
27
DFP
1.0074e+000
7.5404e-001
1.7336e-009
5
21
最速下降法
9.9290e-001
7.2793e-001
1.3266e-009
39
255
BFGS
三次多项式插值
1.0072e+000
7.6715e-001
1.5634e-009
7
27
DFP
1.0074e+000
7.5404e-001
1.7336e-009
5
21
最速下降法
9.9290e-001
7.2793e-001
1.3266e-009
39
255
显然，精确解应为x1=1，z=0.由上表可知，计算效率：DFP>BFGS>最速下降法，而且DFP和BFGS的迭代次数近似，最速下降法的迭代次数远远大于前两者。这是因为最速下降法的收敛阶数为1，开始时收敛非常快，但当接近最优解时，收敛就变慢了。而前两种算法均利用到了二阶导数项，收敛阶数较高，因此总体上收敛速率较快，计算效率更高。对于解的精度，三者近似相同，没有很大的区别。
此外，通过对比可以看出，对于混合二次三次多项式插值和三次多项式插值两种步长算法，他们的最优解以及迭代次数都是一样的。但事实上不是二者没有区别，我通过翻阅matlab的warning消息，得知LineSearchType的控制参数已经不存在了，因而matlab只使用了某一种步长算法。
（三） 比较数值方法和分析方法的精度及效率
首先编写含梯度的函数M文件（程序见四.6），然后选择初值x1=0.6，x2=0.8，分别使用数值方法和分析方法进行计算（程序见四.7），结果如下：
搜索方向
算法
最优解x1
最优解x2
最优值
迭代次数
z调用次数
BFGS
数值方法
1.00720E+00
7.67150E-01
1.56340E-09
7
27
DFP
1.00740E+00
7.54040E-01
1.73360E-09
5
21
最速下降法
9.92900E-01
7.27930E-01
1.32660E-09
39
255
BFGS
分析方法
1.00720E+00
7.67150E-01
1.56340E-09
7
9
DFP
1.00740E+00
7.54040E-01
1.73360E-09
5
7
最速下降法
9.92900E-01
7.27930E-01
1.32660E-09
39
85
由上表可知，分析方法的目标函数调用次数少于数值方法，在进行大规模计算的时候，分析方法不失为数值方法的一种替代。但从本例来看，分析方法对于提高精度的效果不明显，而且对迭代次数也没有减少。
在进行完上述算法的比较之后，我们开始逐个检验之前得到的极小值点，选定fminunc函数和BFGS搜索方向进行分析, 取不同的初值（程序见四.8），结果如下：
初值x01
初值x02
期望解X1
期望解X2
最优解x1
最优解x2
最优值
0.3
0.3
0
0
2.2336e-001
2.8437e-006
1.4678e-013
0.2
-2
0
---
-7.5875e-009
-1.3095e+000
5.2657e-016
1.2
4
1
---
1.0026e+000
3.9964e+000
7.2181e-010
0.3
18
0.5
16
1.0088e+000
1.7980e+001
1.9765e-006
0.45
15.6
0.5
16
4.9682e-001
1.5599e+001
2.6331e-014
在第四步实验（即上表第四行）中，作者本来希望初值（0.3,18）能收敛到期望值（0.5,16），然而却最终得到的最优解为（1,18），虽然也为极值点，但搜索方向却没有按照预期进行。于是进行了第五步的实验，将初值与期望值的偏差进一步缩小，结果获得了期望的最优解。
由上表的最优值，我们可以判定，当x2=0或x1=0或x1=1或x2(1-x1)4=1时，目标函数值为0，又因为目标函数是非负连续的，因而零点即为局部极小值，也为全局极小值。
但这只是局部极小值的充分条件，仅凭此并不能找到所有的局部极小值。然而，通过对图像的观察，可以直观地判断出零点和极小值点是等价的，至于更分析上的探讨，由于本人数学能力有限，无法做出全面的证明。
得出结论：
该函数的极小值点为：x2=0或x1=0或x1=1或x2(1-x1)4=1，极小值点即为零点，因而极小值为0.
对于不同算法的优劣，有如下四点结论：
a) fminunc的迭代次数和目标函数的调用次数较少，效率较高，但最优值的精度较低，与之相反的是，fminsearch的迭代次数和目标函数的调用次数较多，效率比较低，然而精度很高。
b) 对于DFP、BFGS和最速下降法这三种搜索方向，计算效率：DFP>BFGS>最速下降法，而三者的计算精度在数量级上是近似的。
c) 对于混合二次三次多项式插值和三次多项式插值两种步长算法，由于matlab中已经没有了步长的设置参数，因此无法进行比较。
d) 对于数值方法和分析方法计算梯度，结论是分析方法可以大幅减少目标函数的调用次数，但对于提高精度的贡献很小，迭代次数也没有区别。

2. 《数学实验》第二版（问题6）
问题描述：
有一组数据(ti,,yi)(i=1,2,…,33)，其中ti=10(i-1)，yi由表中给出。现在用这组数据拟合函数：f(x,t)=x1+x2e-x4t+x3e-x5t
其中参数x，初值可选为(0.5,1.5,-1,0.01,0.02)，用GN和LM两种方法求解。对yi作一扰动，即yi+ei，ei为(-0.05,0.05)内的随机数，观察并分析迭代收敛是否会变慢。
i
yi
i
yi
i
yi
1
0.844
12
0.718
23
0.478
2
0.908
13
0.685
24
0.467
3
0.932
14
0.658
25
0.457
4
0.936
15
0.628
26
0.448
5
0.925
16
0.603
27
0.438
6
0.908
17
0.580
28
0.431
7
0.881
18
0.558
29
0.424
8
0.850
19
0.538
30
0.420
9
0.818
20
0.522
31
0.414
10
0.784
21
0.506
32
0.411
11
0.751
22
0.490
33
0.406
实验过程：
在matlab中，实现最小二乘拟合有lsqnonlin和lsqcurvefit，鉴于后者的代码更清晰简洁，因而我使用了lsqcurvefit来进行拟合。
又因为在此版本的matlab中，已经不能设置LevenbergMarquardt的控制参数，即不能手动设置非线性最小二乘拟合的下降方向，因此下面我只给出了LM法的拟合结果（程序见四.9、四.10）：
x1
x2
x3
x4
x5
迭代数
z调用次数
norm
0.37541
1.9358
-1.4647
0.01287
0.02212
8
54
5.4649e-005
因此，对于函数f(x,t)=x1+x2e-x4t+x3e-x5t，五个参数的取值分别为0.37541、1.9358、-1.4647、0.01287、0.02212时，与上述数据点的差别最小，拟合度最高。误差平方和为5.4649e-005
然后，对yi作一扰动，即yi+ei，ei为(-0.05,0.05)内的随机数，并重复四次（程序见四.11）所得结果如下：
序号
x1
x2
x3
x4
x5
迭代数
z调用次数
norm
1
0.401
2.399
-1.996
0.015
0.023
9
60
2.2390e-002
2
0.342
1.307
-0.817
0.010
0.025
8
54
2.2705e-002
3
0.407
25.920
-25.500
0.018
0.018
50
306
3.0588e-002
4
0.367
1.540
-1.104
0.012
0.026
11
72
2.1883e-002
由上表可以看出，扰动后输出的误差平方和norm数量级在0.01左右，远远大于未扰动之前的误差。因此，当一个随机小扰动加到y上之后，新的数据点一般会产生较大的偏离，致使数据点与函数的吻合性变差了。不仅如此，当加入扰动后，迭代数与目标函数调用次数也有或多或少的增加，这说明收敛速度变慢了，这也是由于数据点拟合度不高的必然结果。
得出结论：
对于函数f(x,t)=x1+x2e-x4t+x3e-x5t，五个参数的取值分别为0.37541、1.9358、-1.4647、0.01287、0.02212时，拟合度最高。而且在加入扰动后，拟合度和精度大大下降，收敛速度也有所变慢。
三、 实验总结
本次实验是利用matlab优化工具箱进行无约束优化，并进行多元函数极小值的分析及非线性含参量函数的最小二乘拟合。这一部分无论在基本原理上还是代码上都比之前的实验更有难度，如果说之前的实验是灰箱的话，那么我认为本次实验完全是“半黑”箱，各种搜索方向和步长算法完全是照搬代码，而对其内部的原理不甚清楚，这就导致在得出实验结果之后，只能进行表面上的分析而不能从内部原理来解释各种算法的优劣。
当然，对于今后可能的工程应用而言，数学的工具性高于它的理论性，通过这次实验，我也基本掌握了函数极值的基本求解思路及非线性最小二乘的方法，对算法和精度控制的有关参数也能够熟练应用，达到了基本要求。
四、 程序清单
1. 第一题——绘制第一幅三维图和等高线
[x1,x2]=meshgrid(-2:0.001:1,-1:0.001:1);
z=(x1.*x2).^2.*(1-x1).^2.*(1-x1-x2.*(1-x1).^5).^2;
mesh(x1,x2,z)
pause;
contour(x1,x2,z,120)

2. 第一题——绘制第二幅三维图和等高线
[x1,x2]=meshgrid(-8:0.01:10,-5:0.01:5);
z=(x1.*x2).^2.*(1-x1).^2.*(1-x1-x2.*(1-x1).^5).^2;
mesh(x1,x2,z)
pause;
contour(x1,x2,z,220)

3. 第一题——设置无约束优化的函数
function z=profun1(x)
z=(x(1)*x(2))^2*(1-x(1))^2*(1-x(1)-x(2)*(1-x(1))^5)^2;
end

4. 第一题——比较fminunc和fminsearch的优劣
x0=[0.8,0.8];
format shorte
[x1,z1,ef1,out1]=fminunc(@profun1,x0)
[x2,z2,ef2,out2]=fminsearch(@profun1,x0)

5. 第一题——比较不同搜索方向和不同步长的算法精度及效率
x0=[0.6,0.8];
format shorte
opt=optimset('LargeScale','off','MaxFunEvals',1000);
opt1=optimset(opt,'HessUpdate','bfgs');
opt2=optimset(opt,'HessUpdate','dfp');
opt3=optimset(opt,'HessUpdate','steepdesc');
opt4=optimset(opt,'HessUpdate','bfgs','LineSearchType','cubicpoly');
opt5=optimset(opt,'HessUpdate','dfp','LineSearchType','cubicpoly');
opt6=optimset(opt,'HessUpdate','steepdesc','LineSearchType','cubicpoly');
[x,fv,exit,out]=fminunc('profun1',x0,opt1)
pause
[x,fv,exit,out]=fminunc('profun1',x0,opt2)
pause
[x,fv,exit,out]=fminunc('profun1',x0,opt3)
pause
[x,fv,exit,out]=fminunc('profun1',x0,opt4)
pause
[x,fv,exit,out]=fminunc('profun1',x0,opt5)
pause
[x,fv,exit,out]=fminunc('profun1',x0,opt6)
pause

6. 第一题——设置分析方法的函数
function [z,g]=profun11(x);
z=(x(1)*x(2))^2*(1-x(1))^2*(1-x(1)-x(2)*(1-x(1))^5)^2;
if nargout>1
g(1)=2*x(1)*x(2)^2*(1-x(1)-x(2)*(1-x(1))^5)*((2*x(1)^2-3*x(1)+1)*(1-x(1)-x(2)*(1-x(1))^5)+(x(1)^3-2*x(1)^2+x(1))*(5*x(2)*(1-x(1))^4-1));
g(2)=2*x(1)^2*(1-x(1))^2*x(2)*(1-x(1)-x(2)*(1-x(1))^5)^2-2*x(1)^2*(1-x(1))^7*x(2)^2*(1-x(1)-x(2)*(1-x(1))^5);
end

7. 第一题——比较数值方法和分析方法计算梯度
x0=[0.6,0.8];
format shorte
opt=optimset('LargeScale','off','MaxFunEvals',1000);
opt1=optimset(opt,'HessUpdate','bfgs');
opt2=optimset(opt,'HessUpdate','dfp');
opt3=optimset(opt,'HessUpdate','steepdesc');
opt=optimset('LargeScale','off','MaxFunEvals',1000,'GradObj','on');
opt4=optimset(opt,'HessUpdate','bfgs');
opt5=optimset(opt,'HessUpdate','dfp');
opt6=optimset(opt,'HessUpdate','steepdesc');
[x1,fv1,exit,out1]=fminunc('profun1',x0,opt1)
[x2,fv2,exit,out2]=fminunc('profun1',x0,opt2)
[x3,fv3,exit,out3]=fminunc('profun1',x0,opt3)
[x4,fv4,exit,out4]=fminunc('profun11',x0,opt4)
[x5,fv5,exit,out5]=fminunc('profun11',x0,opt5)
[x6,fv6,exit,out6]=fminunc('profun11',x0,opt6)
x=[x1;x2;x3;x4;x5;x6]
fv=[fv1;fv2;fv3;fv4;fv5;fv6]
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]

8. 第一题——选定fminunc函数和BFGS搜索方向, 取不同的初值以寻找极小值点
x0=[10,6]; %x0可以变化%
format shorte
opt=optimset('LargeScale','off','MaxFunEvals',1000,'HessUpdate','bfgs');
[x1,z1,ef1,out1]=fminunc(@profun1,x0,opt)

9. 第二题——设置非线性最小二乘拟合的函数
function f=profun2(x,t)
f=x(1)+x(2)*exp(-x(4)*t)+x(3)*exp(-x(5)*t);
end

10. 第二题——无扰动时，最小二乘拟合
i=[0:32];
t=10*i;
y=[.844 .908 .932 .936 .925 .908 .881 .850 .818 .784 .751 .718 .685 .658 .628 .603 .580 .558 .538 .522 .506 .490 .478 .467 .457 .448 .438 .431 .424 .420 .414 .411 .406];
x0=[0.5,1.5,-1,0.01,0.02];
[x,norm,res,ef,out,lam,jac]=lsqcurvefit('profun2',x0,t,y)

11. 第二题——有扰动时，最小二乘拟合
i=[0:32];
t=10*i;
e=0.1*rand(1,33)-0.05;
y=[.844 .908 .932 .936 .925 .908 .881 .850 .818 .784 .751 .718 .685 .658 .628 .603 .580 .558 .538 .522 .506 .490 .478 .467 .457 .448 .438 .431 .424 .420 .414 .411 .406];
x0=[0.5,1.5,-1,0.01,0.02];
y=y+e;
[x,norm,res,ef,out,lam,jac]=lsqcurvefit(@profun2,x0,t,y)

11