大学数学实验七_无约束优化
7(10)无约束最优化问题
无约束最优化问题
三,极值的充分条件
定理2 充分条件) 定理2 (充分条件) 设函数 z = f ( x , y )在点( x0 , y0 ) 的某邻域内连续, 有一阶及二阶连续偏导数, 的某邻域内连续 有一阶及二阶连续偏导数 又 f x ( x0 , y0 ) = 0, f y ( x0 , y0 ) = 0, 令 fxx ( x0 , y0 ) = A, fxy ( x0 , y0 ) = B, f yy ( x0 , y0 ) = C,
18
无约束最优化问题
作业
习题7.10 (112页 习题7.10 (112页) (A)2. 3.(2) 6. (B) 1. 2. 6.
19
�
一元函数 f ( x , y0 ) 在点 x0 处取得有极小值 处取得有极小值, 表示动点 P ( x , y ) ∈ U ( P0 , δ ),且 P ( x , y )沿直线
17
无约束最优化问题
y = y0上, 并沿该直线 即沿平行于 轴的正负 并沿该直线(即沿平行于 即沿平行于Ox轴的正负
方向)趋向于 方向 趋向于P0 ( x0 , y0 )时, f ( x, y) > f ( x0 , y0 ). 它们的关系是: 它们的关系是 取得极大(小 值 f ( x , y ) 在点 ( x0 , y0 ) 取得极大 小)值 f ( x0 , y )和f ( x , y0 )分别在 y0点和x0点 取得极大(小 值 取得极大 小)值.
下半个圆锥面
x
点取极大值. 也是最大值). 在(0,0)点取极大值 (也是最大值 点取极大值 也是最大值 马鞍面
z
O
y
O
x
y
4
无约束最优化问题
数学实验——无约束优化
实验6无约束优化分1黄浩43实验目的1. 掌握用MATLAB优化工具箱的基本用法，对不同算法进行初步分析、比较2. 练习用无约束优化方法建立和求解实际问题模型（包括非线性最小二乘拟合）。
二、实验内容1. 《数学实验》第二版（问题2.1）问题叙述：取不同的初值计算非线性规划：尽可能求出所有局部极小点，进而找出全局极小点，并对不同算法（搜索方向、步长搜索、数值梯度与分析梯度等）的结果进行分析、比较。
实验过程：首先绘制这个函数的三维图形以及等高线（程序见四.1），结果如下:s M tn0-”19 A8 A1G \ 5 -14 -13 \ 2 A1通过观察这两幅图，可以得到，x2确定时，x1越负，函数值越大，x1确定时，x2绝对值越大，函数值越大。
但对于x1正向偏离0的情况，并没有很好的反映，于是扩大绘图范围，做出下图（程序见四.2）：-1 -10由上面两幅图可见，方程像是一个四角被捏起的花布，而且z的最小值为0< 因此只要求解该方程的零点，即得到了方程的局部极小点，且若将原方程变形为:我们容易发现，该方程的零点为：x2=0或x1=0或x1=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=1, z=0,我们便可以比较出不同算法的精度。
无约束优化
无约束优化一、实验目的1、了解无约束优化的基本算法；2、掌握Matlab优化工具箱的基本用法；3、掌握用Matlab求解无约束优化实际问题。
二、实验要求能够掌握Matlab优化工具箱中fminunc，fminearch，lqnonlin，lqcurvefit的基本用法，能够对控制参数进行设置，能够对不同算法进行选择和比较。
fminunc为无约束优化提供了大型优化和中型优化算法.由option中的参数LargeScale控制：LargeScale=’on’(默认值),使用大型算法LargeScale=’off’,使用中型算法fminunc为中型优化算法的搜索方向提供了3种算法，由option中的参数HeUpdate控制：HeUpdate=’bfg’（默认值），拟牛顿法的BFGS公式；HeUpdate=’dfp’，拟牛顿法的DFP公式；HeUpdate=’teepdec’，最速下降法fminunc为中型优化算法的步长一维搜索提供了两种算法，由option 中参数LineSearchType控制：LineSearchType=’quadcubic’(缺省值)，混合的二次和三次多项式插值；LineSearchType=’cubicpoly’，三次多项式插搜索步长的算法选择（lqnonlin，lqcurvefit）LevenbergMarquardt=‘off’(GN法)LevenbergMarquardt=‘on’(LM法，缺省值)22+2某2+4某1某2+2某2+1)e某1例minf(某)=(4某11、编写M-文件fun1.m:functionf=fun1(某)f=e某p(某(1))某(4某某(1)^2+2某某(2)^2+4某某(1)某某(2)+2某某(2)+1);2、输入M文件myprg3.m如下:某0=[-1,1];某=fminunc('fun1',某0)y=fun1(某)三、实验内容1.求下列函数的极小值点：222f(某)=某1+4某2+9某3-2某1+18某22g(某)=某1+32某2-2某1某2+某1-2某22某2y22、求解min(+)对(a,b)的不同取值如(1,1)和(9,1)，及不同算法(搜索方ab向、步长搜索、数值梯度与分析梯度等)的结果进行分析、比较。
无约束常用优化方法
步长 ，作前进（或后退）试探．如试探成功（目
标函数值有所减小），则按步长序列
，加
大步长（注意每次加大步长都是由初始点算起），直
至试探失败（目标函数值比前一次的有所增加）时，
则取其前一次的步长作为沿这个坐标轴方向搜索的最
优步长，并计算出该方向上的终止点，而后以这个终
止点为始点再进行下一坐标轴方向的搜索，并重复上
处
显然 是二次函数，并且还是正定二次函数，所以 是凸函数且存在唯一全局极小点．为求此极小点，令
即可解得
即
（5.9）
对照基本迭代公式,易知，式（5.9）中的搜索方向
步长因子
方向
是直指点 处近似二次函数
的极小点的方向．此时称此方向为从点 出发的
Newton方向．从初始点开始，每一轮从当前迭代点出发，
沿Newton方向并取步长 的算法称为Newton法．
另外，共轭梯度法不要求精确的直线搜 索．但是，不精确的直线搜索可能导致迭代 出来的向量不再共轭，从而降低方法的效 能．克服的办法是，重设初始点，即把经过 n次迭代得到的Xn作为初始点重新迭代．
五、坐标轮换法
在坐标轮换法中，沿各个坐标轴方向进行一维搜索
时，常选用最优步长法或加速步长法．加速步长法从
初始点出发，沿搜索（坐标轴）方向先取一个较小的
三、共轭方向法
1、概念
通常，我们把从任意点
出发，依次沿某组共轭
方向进行一维搜索的求解最优化问题的方法，叫做共
轭方向法．
2、特点
• 一般地，在n维空间中可以找出n个互相共轭的方向，对于n元正 定二次函数，从任意初始点出发，顺次沿这n个共轭方向最多作n 次直线搜索就可以求得目标函数的极小点．这就是共轭方向法的 算法形成的基本思想．
实验7 无约束优化
数学实验作业（第八周）郭明钊 2012011880 化21一、原子位置问题1、 问题分析：题目中给出了各个原子之间的距离关系，所要求的就是每个原子在平面直角坐标系之中的具体位置。
则可以假设第i 个原子的坐标为(,)i i x y ，且假设第一个原子的位置为坐标原点，即(0,0)，则问题所求就转化为了使得2222[()()]i ji j ij ij x x y y d -+--∑达到最小值时的解。
这样问题就转化为了无约束优化：2222min [()()]i j i j ij ij x x y y d -+--∑，在这里，建立数组x ，其中有50个数，()i x i 为奇数为第i 个原子的横坐标，()i x i 为偶数为第i 个原子的纵坐标。
2、 使用matlab 中的lsqnonlin 函数实现：首先建立函数m 文件function y=distance(x,d)y(1)=(x(2*4-1))^2+(x(2*4))^2-d(1)^2;y(2)=(x(2*12-1))^2+(x(2*12))^2-d(2)^2;y(3)=(x(2*13-1))^2+(x(2*13))^2-d(3)^2;y(4)=(x(2*17-1))^2+(x(2*17))^2-d(4)^2;y(5)=(x(2*21-1))^2+(x(2*21))^2-d(5)^2;y(6)=(x(2*5-1)-x(2*2-1))^2+(x(2*5)-x(2*2))^2-d(6)^2;y(7)=(x(2*16-1)-x(2*2-1))^2+(x(2*16)-x(2*2))^2-d(7)^2;y(8)=(x(2*17-1)-x(2*2-1))^2+(x(2*17)-x(2*2))^2-d(8)^2;y(9)=(x(2*25-1)-x(2*2-1))^2+(x(2*25)-x(2*2))^2-d(9)^2;y(10)=(x(2*5-1)-x(2*3-1))^2+(x(2*5)-x(2*3))^2-d(10)^2;y(11)=(x(2*20-1)-x(2*3-1))^2+(x(2*20)-x(2*3))^2-d(11)^2;y(12)=(x(2*21-1)-x(2*3-1))^2+(x(2*21)-x(2*3))^2-d(12)^2;y(13)=(x(2*24-1)-x(2*3-1))^2+(x(2*24)-x(2*3))^2-d(13)^2;y(14)=(x(2*5-1)-x(2*4-1))^2+(x(2*5)-x(2*4))^2-d(14)^2;y(15)=(x(2*12-1)-x(2*4-1))^2+(x(2*12)-x(2*4))^2-d(15)^2;y(16)=(x(2*24-1)-x(2*4-1))^2+(x(2*24)-x(2*4))^2-d(16)^2;y(17)=(x(2*8-1)-x(2*6-1))^2+(x(2*8)-x(2*6))^2-d(17)^2;y(18)=(x(2*13-1)-x(2*6-1))^2+(x(2*13)-x(2*6))^2-d(18)^2;y(19)=(x(2*19-1)-x(2*6-1))^2+(x(2*19)-x(2*6))^2-d(19)^2;y(20)=(x(2*25-1)-x(2*6-1))^2+(x(2*25)-x(2*6))^2-d(20)^2;y(21)=(x(2*8-1)-x(2*7-1))^2+(x(2*8)-x(2*7))^2-d(21)^2;y(22)=(x(2*14-1)-x(2*7-1))^2+(x(2*14)-x(2*7))^2-d(22)^2;y(23)=(x(2*16-1)-x(2*7-1))^2+(x(2*16)-x(2*7))^2-d(23)^2;y(24)=(x(2*20-1)-x(2*7-1))^2+(x(2*20)-x(2*7))^2-d(24)^2;y(25)=(x(2*21-1)-x(2*7-1))^2+(x(2*21)-x(2*7))^2-d(25)^2;y(26)=(x(2*14-1)-x(2*8-1))^2+(x(2*14)-x(2*8))^2-d(26)^2;y(27)=(x(2*18-1)-x(2*8-1))^2+(x(2*18)-x(2*8))^2-d(27)^2;y(28)=(x(2*13-1)-x(2*9-1))^2+(x(2*13)-x(2*9))^2-d(28)^2;y(29)=(x(2*15-1)-x(2*9-1))^2+(x(2*15)-x(2*9))^2-d(29)^2;y(30)=(x(2*22-1)-x(2*9-1))^2+(x(2*22)-x(2*9))^2-d(30)^2;y(31)=(x(2*11-1)-x(2*10-1))^2+(x(2*11)-x(2*10))^2-d(31)^2;y(32)=(x(2*13-1)-x(2*10-1))^2+(x(2*13)-x(2*10))^2-d(32)^2;y(33)=(x(2*19-1)-x(2*10-1))^2+(x(2*19)-x(2*10))^2-d(33)^2;y(34)=(x(2*20-1)-x(2*10-1))^2+(x(2*20)-x(2*10))^2-d(34)^2;y(35)=(x(2*22-1)-x(2*10-1))^2+(x(2*22)-x(2*10))^2-d(35)^2;y(36)=(x(2*18-1)-x(2*11-1))^2+(x(2*18)-x(2*11))^2-d(36)^2;y(37)=(x(2*25-1)-x(2*11-1))^2+(x(2*25)-x(2*11))^2-d(37)^2;y(38)=(x(2*15-1)-x(2*12-1))^2+(x(2*15)-x(2*12))^2-d(38)^2;y(39)=(x(2*17-1)-x(2*12-1))^2+(x(2*17)-x(2*12))^2-d(39)^2;y(40)=(x(2*15-1)-x(2*13-1))^2+(x(2*15)-x(2*13))^2-d(40)^2;y(41)=(x(2*19-1)-x(2*13-1))^2+(x(2*19)-x(2*13))^2-d(41)^2;y(42)=(x(2*15-1)-x(2*14-1))^2+(x(2*15)-x(2*14))^2-d(42)^2;y(43)=(x(2*16-1)-x(2*14-1))^2+(x(2*16)-x(2*14))^2-d(43)^2;y(44)=(x(2*20-1)-x(2*16-1))^2+(x(2*20)-x(2*16))^2-d(44)^2;y(45)=(x(2*23-1)-x(2*16-1))^2+(x(2*23)-x(2*16))^2-d(45)^2;y(46)=(x(2*18-1)-x(2*17-1))^2+(x(2*18)-x(2*17))^2-d(46)^2;y(47)=(x(2*19-1)-x(2*17-1))^2+(x(2*19)-x(2*17))^2-d(47)^2;y(48)=(x(2*20-1)-x(2*19-1))^2+(x(2*20)-x(2*19))^2-d(48)^2;y(49)=(x(2*23-1)-x(2*19-1))^2+(x(2*23)-x(2*19))^2-d(49)^2;y(50)=(x(2*24-1)-x(2*19-1))^2+(x(2*24)-x(2*19))^2-d(50)^2;y(51)=(x(2*23-1)-x(2*21-1))^2+(x(2*23)-x(2*21))^2-d(51)^2;y(52)=(x(2*23-1)-x(2*22-1))^2+(x(2*23)-x(2*22))^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.3746 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.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,47)]; %设初值[x,norms,res]=lsqnonlin(@distance,x0,[],[],[],d)a=reshape(x,2,25)'b=a(:,1)';c=a(:,2)';plot(b,c,'*') %在坐标系中显示出各个原子的位置3、结果如下：误差平方和：norms = 0.1625误差向量res =8.4224e-002 -5.6716e-002 5.4597e-0025.6860e-003 -1.1847e-001 5.8619e-0023.8879e-002 1.1819e-001 2.2805e-0023.8320e-002 -3.4772e-003 -3.9762e-0021.1452e-002 5.3897e-002 -4.3139e-0022.1071e-002 -4.4450e-002 9.6968e-003-5.6455e-002 4.0990e-002 2.7157e-0021.4843e-001 1.5864e-002 -3.7006e-003-1.3585e-001 2.6711e-002 -3.8390e-0039.0472e-003 2.9573e-002 -4.3175e-003-2.5224e-003 -4.7600e-004 1.2934e-002-2.4690e-002 7.6172e-003 2.5536e-003-5.2482e-003 7.3814e-002 -3.1225e-003-6.2153e-002 4.1666e-002 1.3685e-0016.4955e-002 2.3253e-003 -6.2132e-002-1.5985e-003 -5.7399e-002 -1.5199e-0021.3400e-002 -1.8521e-002 1.2003e-0013.1653e-003各个原子的坐标数值（第一列为横坐标，第二列为纵坐标）ans =0 00.2496 1.37301.6630 1.60100.7720 0.64120.7953 1.17011.2807 0.94921.4662 0.83651.2996 0.50230.3616 0.50871.1538 0.82271.1655 1.3985-0.3685 -0.03120.2287 0.81570.9857 0.85610.5856 0.44400.6025 1.9132-0.2599 1.35380.5372 0.16430.7983 1.36701.1261 1.01080.8162 0.91311.6189 0.70121.0248 0.15481.4666 0.46970.8870 1.0702 显示在坐标系中：当改变初值时x0=[zeros(1,2),ones(1,48)]; %设初值得norms =0.2298ans =0 00.7773 1.54251.2131 1.73990.8442 0.52390.6526 0.98631.4017 1.14390.7255 0.70991.0215 0.8220-0.0902 0.93201.2946 1.03960.7993 1.3167-0.4498 0.26080.4560 0.67820.6244 1.06130.1326 1.11251.2374 1.98050.5434 1.25021.2808 0.02161.3079 0.47031.0750 0.94510.2150 1.25031.0575 0.51940.1391 0.52860.3273 1.01521.0902 0.9464当换成其他初值时，答案也会明显不同。
matlab实验7 无约束优化
热动71 马千里 970669实验七 无约束优化实验目的1. 掌握MATLAB 优化工具箱的基本用法，对不同算法进行初步分析、比较。
2. 练习实际问题的非线性最小二乘拟合。
实验内容3. 求解)12424(min 22122211++++x x x x x e x ，初值(-1,1)，对不同算法的结果进行分析、比较。
解：编制函数fun.mfunction y=fun(x)y=exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+2*x(2)+1);经过实验，用不同的算法fminu 得到的结果都是(0.5, -1)，所不同的是迭代的次数有不同。
以上结果表明，MA TLAB 的缺省值对于此题优化效果较好，迭代次数最少。
6.《中国统计年鉴(1995)》给出表8的数据，试据此拟合生产函数中的参数，如何看待用线性最小二乘法和非线性最小二乘法拟合的结果。
解：根据Gobb-Douglas 生产函数，用Q, K, L 分别表示产值、资金、劳动力。
βαL aK L)Q(K,=a. 考虑线性最小二乘拟合。
对生产函数取对数L K a Q ln ln ln ln βα++=设R=[1 lnK lnL] x=[lna α β]’ y=[lnQ]R=1.0000 -1.3988 1.57231.0000 -1.0829 1.60691.0000 -0.9556 1.63481.0000 -0.8389 1.66361.0000 -0.5987 1.69261.0000 -0.4951 1.71071.0000 -0.4394 1.73591.0000 -0.2854 1.76401.0000 -0.0371 1.78221.0000 0.4053 1.7954 1.0000 0.6389 1.8160则Rx=y 代表一组超定方程组。
在MATLAB 下用x=R\y 得到最小二乘意义下的解 x=-3.23600.62082.3728即a=exp(-3.2360)=0.0393 α=0.6208 β=2.3728b.非线性模型function f=pp1(c)Q=[0.7171 0.8964 1.0202 1.1962 1.4928 1.6909 1.8531 2.1618 2.6635 3.45154.5006];K=[0.2469 0.3386 0.3846 0.4322 0.5495 0.6095 0.6444 0.7517 0.9636 1.49981.8944];L=[4.8179 4.9873 5.1282 5.2783 5.4334 5.5329 5.6740 5.8360 5.9432 6.02206.1470];f=Q-c(1)*k.^c(2).*L.^c(3);c0=[1 1 1]c=leastsq(‘pp1’, c0)y=sum(pp1(c).*pp1(c)) 计算误差平方和得到c= 0.0357 0.6300 2.4290y= 0.0394即a=0.0357 α=0.6300 β=2.4290再对线性模型计算误差平方和以做比较y=sum(pp1([0.0393 0.6208 2.3728]).* pp1([0.0393 0.6208 2.3728]))得到y=0.0441将以上结果列表如下是在不同意义下的最小二乘解。
无约束优化
数学实验七无约束优化化21 张冶5.某分子由25个原子组成，并且已经通过实验测量得到了其中某些原子对之间的距离（假设在平面结构上讨论），确定每个原子的位置关系。
解：由已知得，找到25个原子的相对坐标，使其满足题中所给数据，使其误差最小。
令第一个点的坐标为（0,0），其余点的坐标为(x i,y j), 则只需找出余下24个点的坐标使∑[(x i−x j)2+(y i−y j)2−d ij2]2i,j最小即可。
设x=[x2 y2 x3 y3 x4 y4 ······x12 y12······x24 y24]编写M文件如下：function y=yuanzi(x,d)y(1)=(x(2*4-3))^2+(x(2*4-2))^2-d(1)^2;y(2)=(x(2*12-3))^2+(x(2*12-2))^2-d(2)^2;y(3)=(x(2*13-3))^2+(x(2*13-2))^2-d(3)^2;y(4)=(x(2*17-3))^2+(x(2*17-2))^2-d(4)^2;y(5)=(x(2*21-3))^2+(x(2*21-2))^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-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*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)-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(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*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))^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)=(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*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-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.75110.4439 0.8363 0.3208 0.1574 1.2736 0.5781 0.9254 0.6401 0.2467 0.4727 1.3840 0.4366 1.03071.3904 0.5725 0.7660 0.4394 1.0952 1.0422 1.8255 1.4325 1.0851 0.4995 1.2277 1.1271 0.70600.8052];x0=ones(1,48);[x,norms]=lsqnonlin(@yuanzi,x0,[],[],[],d); x=[0,0,x] ;D=reshape(x,2,25)'normsfor i=1:25plot(D(i,1),D(i,2),'*')hold on;end得到如下结果：D =0 01.7141 -0.04411.8719 -1.06420.8171 -0.55351.2606 -0.33781.5533 -0.56431.0677 -0.34581.0960 -0.73060.1433 -0.25541.4056 -0.53960.8949 -0.23200.0899 0.55620.5029 -0.70131.1012 -0.28760.5994 -0.29362.1994 -0.41771.3519 -0.07831.4262 -1.50811.2498 -1.17381.1860 -0.67640.7593 -1.02011.1465 -0.92950.3792 -0.40761.3248 -0.05561.2114 -0.5142norms =0.3008画出分布图：并且经过检验，发现当初值x0不同时，可能得到不同的结果。
Exp07_实验7_无约束优化
k k T k k k k
T
(x ) f
k T
k
3.2 Broyden-Fletcher-Goldfarb-Shanno(BFGS)公式
G
k
f
k
(f
k T
k T
)
(f
) x
k T
k
k
G x (x ) G (x ) G x
k
k
k
k T k
k
k T
k
k
H
k
(1
(f ) H f (x ) f
x
k 1
x [ F ( x )]
k
k
1
F (x )
k
F ( x) f ( x)
3 拟Newton方法 目的 思路
x
k 1 k
不计算Hessian阵，克服病态、不正定、计 算复杂等缺陷，同时保持收敛较快的优点 回顾解方程组 F(x)=0的拟牛顿法
k 1
x [ F ( x )]
2 Newton方法 将f(xk+1)在xk点作泰勒展开至二阶项，用d替代dk
f (x
k 1
) f (x d ) f (x ) f
k k
T
( x )d
k
1 2
d f ( x )d
T 2 k
求d使f(xk+1)极小右端对d导数为0 f ( x k ) 2 f ( x k ) d 0 牛顿方程 牛顿方向 迭代格式
k 1
x
x d
k k
k
f (x
k 1
) f (x )
k
无约束优化
1. 用黑塞矩阵判断驻点的性质 X 已知函数f(X)的驻点 ，可以利用驻点处的黑塞矩阵来 判断驻点的性质： （1）若 2 f ( X ) 是正定的，则驻点 X 是极小点； 2 （2）若 f ( X ) 是负定的，则驻点 X 是极大点； （3）若 2 f ( X ) 是不定的，则驻点 X 不是极值点； 2 （4）若 f ( X ) 是半定的，则驻点 X 可能是极值点； 也可能不是极值点，视高阶导数情况而定.
4 1 3 2 2 3 2 1 2
2 3
2 4 x13 2 x1 x2 x3 2 2 f ( X ) 6 x2 x1 4 x3 6 x3 4 x2 2 x1 x3
12 x12 2 x2 2 x1 2 x3 2 f ( X ) 2 x1 12 x2 4 2 x3 4 6 2 x 1
无约束优化
标准形式：
其中
X R
min f X n
X R
f : R n R1
X R
max f X = min [ f X ] n n
数学预备知识
1 梯度
定义1 设f(X)是定义在n维欧氏空间的
R n 上的
可微函数，我们称
(
为函数 f(X)在点X处的梯度，记为▽ f(X). 1，梯度方向是函数在点X处增长最快的方向， 负梯度方向是函数在点X处下降最快的方向. 2，梯度的模是函数沿这一方向的变化率.
2 极值点的必要条件和充分条件
n X R 定义1 对于问题（1），设 是任一给定点，P是
是非零向量，若存在一个数 0，使得对于任意 (0, ) 都有 f ( X P) f ( X ) ,则称P是f(X)在 X 处的下降方向.
7无约束最优化的解析法
第七章 无约束最优化的解析法本章主要内容：最速下降法及其收敛性与收敛速度 Newton 切线法及其收敛性与收敛速度 阻尼Newton 法 共轭梯度法及其收敛性 变度量法、最小二乘法教学目的及要求：掌握最速下降法并理解其收敛性与收敛速度，掌握Newton 切线法并理解其收敛性与收敛速度，了解阻尼Newton 法；掌握共轭梯度法并理解其收敛性；了解变度量法、最小二乘法。
教学重点：最速下降法．教学难点：变度量法．教学方法：启发式．教学手段：多媒体演示、演讲与板书相结合．教学时间：6学时．教学内容：§7.1 最速下降法考虑无约束最优化问题m i n ()f x ， （7.1.1）其中:n f R R →具有一阶连续偏导数．算法7-1（最速下降法）Step1 选取初始数据．选取初始点0x ，给定允许误差0ε>，令0k =． Step2 检查是否满足终止准则．计算()k f x ∇，若()k f x ε∇<，迭代终止，k x 为问题（7.1.1）的近似最优解；否则，转Step3．Step3 进行一维搜索．取()k k d f x =-∇，求k λ和1k x +，使得()min ()k k k k k f x d f x d λλλ≥+=+， 1k k k k x x d λ+=+．令:1k k =+，返回Step2．特别地，考虑1m i n ()2T T f x x Qx b x c =++， （7.1.2） 其中,n n n x R Q R ⨯∈∈为正定矩阵，,n b R c R ∈∈．设第k 次迭代点为k x ，从点k x 出发沿()k f x -∇作一维搜索，得1()k k k k x x f x λ+=-∇，其中k λ为最优步长．根据定理 6.1.1，有1()()0T k k f x f x +∇∇=．而(),n f x Q x b x R∇=+∀∈， 所以1()()()k k k k f x f x Q f x λ+∇=∇-∇，从而(()())()0T k k k k f x Q f x f x λ∇-∇∇=，而Q 正定，即()()0T k k f x Q f x ∇∇>，故由上式解出()()()()T k k k T k k f x f x f x Q f x λ∇∇=∇∇， （7.1.3） 于是1()()()()()T k k k k k T k k f x f x x x f x f x Q f x +∇∇=-∇∇∇， （7.1.4） 这是最速下降法用于问题（7.1.2）的迭代公式．例1 用最速下降法求解问题2212min ()4f x x x =+， （7.1.5）其中12(,)T x x x =．取初始点(0)(1,1)T x =，允许误差0.1ε=．解 问题（7.1.5）中的f 是正定二次函数，且800,,0020Q b c ⎛⎫⎛⎫=== ⎪ ⎪⎝⎭⎝⎭．f 在点12(,)T x x x =处的梯度12()(8,2)T f x x x ∇=．第一次迭代：令搜索方向(0)(0)()(8,2)T d f x =-∇=--，(0)d ε==>，从点(0)x 出发沿(0)d 作一维搜索，由（7.1.3）式和（7.1.4）式有0680.130769520λ==， (1)(1,1)0.130769(8,2)(0.046152,0.738462)T T T x =+--=-．第二次迭代：令(1)(1)()(0.369216, 1.476924)T d f x =-∇=-，(1) 1.522375d ε==>，从点(1)x 出发沿(1)d 作一维搜索，按（7.1.4）式得(2)(0.101537,0.147682)T x =．第三次迭代：令(2)(2)()(0.812296,0.295364)T d f x =-∇=--，(2)0.864329d ε==>，按（7.1.4）式求得(3)(0.009747,0.107217)T x =-．第四次迭代：令(3)(3)()(0.077976,0.214434)T d f x =-∇=-，(3)0.228171d ε==>，按（7.1.4）式求得(4)(0.019126,0.027816)T x =．第五次迭代：令(4)(4)()(0.153008,0.055632)T d f x =-∇=--，(4)0.162807d ε==>，按（7.1.4）式求得(5)(0.001835,0.020195)T x =-．此时，(5)()f x ε∇=<，已满足精度要求，故得问题（7.1.5）的近似最优解(5)(0.001835,0.020195)T x =-．实际上问题（7.1.5）的最优解为(0,0)T x =．定理7.1.1 设:n f R R →具有一阶连续偏导数，0n x R ∈，记0()f x α=，假定水平集(,)S f α有界，令{}k x 是由最速下降法求解问题（7.1.1）产生的点列，则（1）当{}k x 是有穷点列时，其最后一个点是f 的平稳点；（2）当{}k x 是无穷点列时，它必有极限点，并且任一极限点都是f 的平稳点．定理7.1.2 设:n f R R →具有二阶连续偏导数，由最速下降法解问题（7.1.1）产生的点列{}k x 收敛于x ．若存在0ε>和0M m >>，使得当x x ε-<时，有222(),T n m y y f x y M y y R ≤∇≤∀∈， （7.1.7） 则{}k x 线性收敛于x ．§7.2 Newton 法21()()()()()()()()2T T k k k k k k f x x f x f x x x x x f x x x ϕ≈=+∇-+-∇-． 令 ()0x ϕ∇=，即2()()()0k k k f x f x x x ∇+∇-=，解之得211[()]()k k k k x x f x f x -+=-∇∇．算法7-2（Newton 法）Step1 选取初始数据．选取初始点0x ，给定允许误差0ε>，令0k =． Step2 检查是否满足终止准则．计算()k f x ∇，若()k f x ε∇<，迭代终止，k x 为问题（7.1.1）的近似最优解；否则，转Step3．Step3 构造Newton 方向．计算21[()]k f x -∇，取21[()]()k k k d f x f x -=-∇∇． Step4 求下一个迭代点．令1k k k x x d +=+，:1k k =+，返回Step2．例2 用Newton 法求解问题（7.1.5），仍取初始点(0)(1,1)T x =，允许误差0.1ε=．解 (0)()(8,2)T f x ∇=，2(0)80()02f x ⎛⎫∇= ⎪⎝⎭，故2(0)11/80[()]01/2f x -⎛⎫∇= ⎪⎝⎭，(0)2(0)1(0)1[()]()1d f x f x -⎛⎫=-∇∇=- ⎪⎝⎭，(1)(0)(0)(1,1)(1,1)(0,0)T T T x x d =+=-=．由于 (1)()00.1f x ∇=<，迭代结束，得(1)x 为问题（7.1.5）的最优解．定理7.2.1 设:n f R R →具有三阶连续偏导数，,()0n x R f x ∈∇=，若存在0ε>和0m >，使得当x x ε-≤时，有22(),T n m y y f x y y R ≤∇∀∈， （7.2.2） 则当初始点0x 充分接近x 时，由Newton 法解问题（7.1.1）产生的点列{}k x 收敛于x ，并有二阶收敛速度．算法7-3（阻尼Newton 法）Step1 选取初始数据．选取初始点0x ，给定允许误差0ε>，令0k =． Step2 检查是否满足终止准则．计算()k f x ∇，若()k f x ε∇<，迭代终止，k x 为问题（7.1.1）的近似最优解；否则，转Step3．Step3 构造Newton 方向．计算21[()]k f x -∇，取21[()]()k k k d f x f x -=-∇∇． Step4 进行一维搜索．求k λ和1k x +，使得()min ()k k k k k f x d f x d λλλ≥+=+， 1k k k k x x d λ+=+．令:1k k =+，返回Step2．例3 用阻尼Newton 法求解下面问题：222121min ()(1)2()f x x x x =-+-，（7.2.6） 其中12(,)T x x x =．取初始点(0)(0,0)T x =，允许误差0.1ε=．解 第一次迭代：212112212(1)8()()4()x x x x f x x x ⎛⎫----∇= ⎪-⎝⎭， 22212111168()28()84x x x x f x x ⎛⎫--+-∇= ⎪-⎝⎭， 故 (0)(0)()(2,0),()2T f x f x ε∇=-∇=>，2(0)2(0)1201/20(),[()]0401/4f x f x -⎛⎫⎛⎫∇=∇= ⎪ ⎪⎝⎭⎝⎭．于是，Newton 方向 (0)2(0)1(0)[()]()(1,0)T d f x f x -=-∇∇=，从(0)x 出发沿(0)d 作一维搜索，即求(0)(0)2400min ()min(1)2f x d λλλλλ≥≥+=-+的最优解，得到012λ=．令 (1)(0)(0)0(1/2,0)T x x d λ=+=．(1)(1)()(0,1),()1T f x f x ε∇=-∇=>．第二次迭代：2(1)2(1)184111(),[()]48124f x f x --⎛⎫⎛⎫∇=∇= ⎪ ⎪-⎝⎭⎝⎭， (1)2(1)1(1)[()]()(1/4,1/2)T d f x f x -=-∇∇=．从(1)x 出发沿(1)d 作一维搜索，即求(1)(1)24001min ()min [8(2)(2)]128f x d λλλλλ≥≥+=-+- 的最优解，得到12λ=．令(2)(1)(1)1(1,1)T x x d λ=+=．(1)(1)()(0,1),()1T f x f x ε∇=-∇=>．此时，(2)(2)()(0,0),()0T f x f x ε∇=∇=<，得问题（7.2.6）的最优解为(2)(1,1)T x =，这是惟一的最优解．定理7.2.2 设:n f R R →具有二阶连续偏导数，0n x R ∈，记0()f x α=，假定水平集(,)S f α有界，并且对一切2,()n x R f x ∈∇正定．若{}k x 是由阻尼Newton法求解问题（7.1.1）产生的点列，则（1）当{}k x 是有穷点列时，其最后一个点是f 的唯一极小点；（2）当{}k x 是无穷点列时，它必收敛于f 的唯一极小点．§7.3 共轭梯度法最速下降法和Newton 法是最基本的无约束最优化方法，它们的特性各异：前者计算量较小而收敛速度慢；后者虽然收敛速度快，但需要计算目标函数的Hesse 矩阵及其逆矩阵，故计算量大．本节介绍一类无需计算二阶导数并且收敛速度快的方法．定义 设n n Q R ⨯∈为正定矩阵．若n R 中的向量组011,,,m d d d - 满足0,,0,1,,1,T i j d Qd i j m i j =∀=-≠ ，则称011,,,m d d d - 是Q 共轭的．定理7.3.1 设n n Q R ⨯∈是正定矩阵，n R 中非零向量组011,,,m d d d - 是Q 共轭的，则这m 个向量线性无关．定理7.3.2 设n p R ∈，011,,,n d d d - 是n R 中线性无关的向量组，若p 与每个i d 都正交，则0p =．考虑正定二次函数的无约束最优化问题1min ()2T T f x x Qx b x c =++， （7.3.3） 其中n n Q R ⨯∈为正定矩阵，,n b R c R ∈∈．定理7.3.3 设n n Q R ⨯∈为正定矩阵，011,,,n d d d - 是n R 中一组Q 共轭的非零向量．对于问题（7.3.3），若从任意点0n x R ∈出发依次沿011,,,n d d d - 进行一维搜索，则至多经过n 次迭代可得问题（7.3.3）的最优解．算法7-4（共轭方向法）给定一个正定矩阵n n Q R ⨯∈．Step1 选取初始数据．选取初始点0x ，给定允许误差0ε>． Step2 选取初始搜索方向．计算0()f x ∇，求出0d ，使00()0T f x d ∇<，令0k =． Step3 检查是否满足终止准则．若()k f x ε∇<，迭代终止；否则，转Step4． Step4 进行一维搜索．求k λ和1k x +，使得()min ()k k k k k f x d f x d λλλ≥+=+， 1k k k k x x d λ+=+．Step5 选取搜索方向．求1k d +使10,0,1,,T k j d Qd j k +== ，令:1k k =+，返回Step3．如果用共轭方向法求解正定二次函数的无约束最优化问题1min ()2T T f x x Qx b x c =++， （7.3.3） 其中n n Q R ⨯∈为正定矩阵，,n b R c R ∈∈（此时算法中的正定矩阵应与二次函数的正定矩阵一致），那么容易推出迭代公式为1()T k k k k k T k kf x d x x d d Qd +∇=-． （7.3.10） 对于求解问题（7.1.1），我们还有如下一些方法．Fletcher-Reeves （FR ）公式：212()()k k k f x f x α+∇=∇；Dixon-Myers （DM ）公式：11()()()T k k k T k k f x f x d f x α++∇∇=-∇； Polak-Ribiere-Polyak （PRP ）公式：11()[()()]()()T k k k k T k k f x f x f x f x f x α++∇∇-∇=∇∇． 算法7-5（FR 共轭梯度法）Step1 选取初始数据．选取初始点0x ，给定允许误差0ε>． Step2 检查是否满足终止准则．计算0()f x ∇，若0()f x ε∇<，迭代终止，0x为（7.1.1）的近似最优解；否则，转Step3．Step3 构造初始搜索方向．计算00(),0d f x k =-∇=．Step4 进行一维搜索．求k λ和1k x +，使得()min ()k k k k k f x d f x d λλλ≥+=+， 1k k k k x x d λ+=+．Step5 检查是否满足终止准则．计算1()k f x +∇，若1()k f x ε+∇<，迭代终止，1k x +为（7.1.1）的近似最优解；否则，转Step6．Step6 检查迭代次数．若1k n +=，令0:n x x =，返回Step3；否则，转Step7．Step7 构造共轭方向．用FR 公式取11()k k k k d f x d α++=-∇+，212()()k k k f x f x α+∇=∇，令:1k k =+，返回Step4． 注意，如果算法7-4的Step7中k α的形式改为DM 公式或PRP 公式，则分别得到DM 共轭梯度法和PRP 共轭梯度法．定理7.3.5 设:n f R R →具有一阶连续偏导数，0n x R ∈，记0()f x α=，并假设水平集(,)S f α有界．若{}k x 是由共轭梯度法（包括任何一种仅仅与算法7-5中k α的形式不同的共轭梯度法）解问题（7.1.1）所产生的点列，则（1）当{}k x 是有穷点列时，其最后一个点是f 的平稳点；（2）当{}k x 是无穷点列时，它必有极限点，并且任一极限点都是f 的平稳点．可以证明：共轭梯度法产生的点列{}k x 是n 步二阶收敛的，即2lim sup k n k k x x q x x +→∞-≤<∞-．例4 用FR 共轭梯度法求解问题（7.2.6）222121min ()(1)2()f x x x x =-+-，仍取初始点(0)(0,0)T x =，允许误差0.1ε=．222121min ()(1)2()f x x x x =-+-解 因为212112212(1)8()()4()x x x x f x x x ⎛⎫----∇= ⎪-⎝⎭， 所以(0)(0)()(2,0),()2T f x f x ε∇=-∇=>．令(0)(0)()(2,0)T d f x =-∇=，从(0)x 出发，沿(0)d 进行一维搜索，得(1)(0)(0)001/4,(1/2,0)T x x d λλ==+=．从而(1)(1)()(0,1),()1T f x f x ε∇=-∇=>．由FR 公式有2(1)02(0)()14()f x f x α∇==∇， 因此，新的搜索方向为 (1)(1)(0)0()(1/2,1)T d f x d α=-∇+=．从(1)x 出发，沿(1)d 进行一维搜索，得(2)(1)(1)111,(1,1)T x x d λλ==+=．此时(2)(2)()(0,0),()0T f x f x ε∇=∇=<．得问题（7.2.6）的最优解为(2)(1,1)T x =．§7.4 变度量法前面介绍的最速下降法和阻尼Newton 法，它们的迭代公式可以统一为1,(),k k k k k k k x x d d G f x λ+=+⎧⎨=-∇⎩ （7.4.1）其中,n n k k G R λ⨯∈是从点k x 出发沿k d 进行一维搜索的最优步长．当k n G I =（n 阶单位矩阵）时，（7.4.1）式即为最速下降法的迭代公式；当21[()]k k G f x -=∇时，（7.4.1）式就是阻尼Newton 法的迭代公式．因此，如果能够使k G 的选取既不需要计算Hesse 矩阵及其逆矩阵，又能很好地近似于21[()]k f x -∇，则由（7.4.1）式确定的迭代算法将会保持Newton 法收敛速度快的优点，同时又具有计算简单的特性．设问题（7.1.1）中目标函数f 具有二阶连续偏导数，且2()k f x ∇正定．为使由（7.4.1）中第2式确定的搜索方向是f 在点k x 处的下降方向，根据定理1.2.3，应当要求()0T k k f x d ∇<，或即()()0T k k k f x G f x ∇∇>，所以我们应要求（7.4.1）式中的k G 是正定矩阵．设在第1k +次迭代后得到1k x +，将f 在点1k x +处作Taylor 展开，取二阶近似，得21111111()()()()()()()2TT k k k k k k f x f x f x x x x x f x x x ++++++≈+∇-+-∇-， 对上式两边求梯度，有2111()()()()k k k f x f x f x x x +++∇≈∇+∇-，令k x x =，得到2111()()()()k k k k k f x f x f x x x +++∇-∇≈∇-，即21111[()][()()]k k k k k f x f x f x x x -+++∇∇-∇≈-．易知，当f 为正定二次函数时，上式成为等式，即 21111[()][()()]k k k k k f x f x f x x x -+++∇∇-∇=-． （7.4.2）因为具有正定Hesse 矩阵的函数在极小点附近可用二次函数很好地近似，所以如果我们迫使1k G +满足类似于（7.4.2）式的关系式，即令111[()()]k k k k k G f x f x x x +++∇-∇=-， （7.4.3）则1k G +就可以很好地近似于211[()]k f x -+∇．因此称（7.4.3）式为拟Newton 条件．为方便起见，记1k k k x x x +∆=-，1()()k k k g f x f x +∆=∇-∇，1k k k G G G +∆=-，并称k G ∆为校正矩阵，则拟Newton 条件可以写成k k k k k G g x G g ∆∆=∆-∆． （7.4.4）综上所述，我们得到如下的一类算法．算法7-6（拟Newton 法）Step1 选取初始数据．选取初始点0n x R ∈和初始矩阵0G ，要求0G 为正定矩阵（可取0n G I =），给定允许误差0ε>，令0k =．Step2 检查是否满足终止准则．计算()k f x ∇，若()k f x ε∇<，迭代终止，k x 为问题（7.1.1）的近似最优解；否则，转Step3．Step3 构造搜索方向．令()k k k d G f x =-∇．Step4 进行一维搜索．求k λ和1k x +，使得()min ()k k k k k f x d f x d λλλ≥+=+， 1k k k k x x d λ+=+．Step5 产生校正矩阵．求出满足（7.4.4）的校正矩阵k G ∆，令1,:1k k k G G G k k +=+∆=+，返回Step2．算法7-7（DFP 法）Step1 选取初始数据．选取初始点0x 和初始矩阵0n G I =，给定允许误差0ε>．Step2 检查是否满足终止准则．计算0()f x ∇，若0()f x ε∇<，迭代终止，0x 为（7.1.1）的近似最优解；否则，转Step3．Step3 构造初始DFP 方向．取00()d f x =-∇，令0k =．Step4 进行一维搜索．求出k λ和1k x +，使得01()min (),.k k k k k k k k k f x d f x d x x d λλλλ≥++=+⎧⎪⎨=+⎪⎩ Step5 检查是否满足终止准则．计算1()k f x +∇，若1()k f x ε+∇<，迭代终止，1k x +为（7.1.1）的近似最优解；否则，转Step6．Step6 检查迭代次数．若1k n +=，令0:n x x =，返回Step3；否则，转Step7．Step7 构造DFP 方向．用DFP 公式1T T k k k k k k k k T T k k k k kx x G g g G G G x g g G g +∆∆∆∆=+-∆∆∆∆，算出1k G +，取111()k k k d G f x +++=-∇，令:1k k =+，返回Step4．定理7.4.4 设:n f R R →具有一阶连续偏导数，0n x R ∈，记0()f x α=，并假设水平集(,)S f α有界．若{}k x 是由DFP 法求解问题（7.1.1）产生的点列，则（1）当{}k x 是有穷点列时，其最后一个点是f 的平稳点；（2）当{}k x 是无穷点列时，它必有极限点，并且其任一极限点都是f 的平稳点．可以证明DFP 法具有超线性收敛性．算法7-8（BFGS 法）Step1 选取初始数据．选取初始点0x 和初始矩阵0n G I =，给定允许误差0ε>．Step2 检查是否满足终止准则．计算0()f x ∇，若0()f x ε∇<，迭代终止，0x 为（7.1.1）的近似最优解；否则，转Step3．Step3 构造初始BFGS 方向．取00()d f x =-∇，令0k =．Step4 进行一维搜索．求出k λ和1k x +，使得01()min (),.k k k k k k k k k f x d f x d x x d λλλλ≥++=+⎧⎪⎨=+⎪⎩ Step5 检查是否满足终止准则．计算1()k f x +∇，若1()k f x ε+∇<，迭代终止，1k x +为（7.1.1）的近似最优解；否则，转Step6．Step6 检查迭代次数．若1k n +=，令0:n x x =，返回Step3；否则，转Step7． Step7 构造BFGS 方向．用BFGS 公式111()T T T T k k k k k k k k k k k k k T T T k k k k k kx x g G g G G x g G G g x x g x g x g +⎛⎫∆∆∆∆=++-∆∆+∆∆ ⎪∆∆∆∆∆∆⎝⎭， 算出1k G +，取111()k k k d G f x +++=-∇，令:1k k =+，返回Step4．可以证明BFGS 法具有超线性收敛性．§7.5 最小二乘法在实际应用中，我们经常遇到目标函数为若干个函数的平方和的最优化问题： 21min ()()mi i s x f x ==∑， （7.5.1）其中n x R ∈，一般假设m n ≥，这类问题称为最小二乘问题．当每个()i f x 都是线性函数时，问题（7.5.1）称为线性最小二乘问题，否则，称为非线性最小二乘问题．由于最小二乘问题相对于一般无约束最优化问题而言具有特殊形式，因此除能运用本章前面介绍的一般求解方法外，还应有更为简便有效的方法．一、线性最小二乘法当()i f x 为线性函数时，即1(),1,2,,ni ij j i j f x a x b i m ==-=∑ ，问题（7.5.1）就成为线性最小二乘问题．如令11111,n m mn m a a b A b a a b ⎛⎫⎛⎫ ⎪ ⎪== ⎪ ⎪ ⎪ ⎪⎝⎭⎝⎭，12()((),(),,())T m f x f x f x f x = ，则 2()()()()()T T s x f x f x Ax b Ax b Ax b ==--=-， 从而问题（7.5.1）可表示为2min ()s x Ax b =-． （7.5.2） 因为()()()2T T T T T s x Ax b Ax b x A Ax b Ax b b =--=-+，所以()s x 为二次函数，且2()22,()2T T T s x A Ax A b s x A A ∇=-∇=．显然，对一切n y R ∈，均有20T T y A Ay Ay =≥，即知，T A A 是半正定矩阵，从而由定理2.3.7知，()s x 是n R 上的凸函数．定理7.5.1 n x R ∈是问题（7.5.2）的最优解的充要条件是x 满足方程组T T A Ax A b =． （7.5.3） 容易知道，矩阵T A A 正定的充要条件是对于一切非零向量n y R ∈，有20T T y A Ay Ay =>．记1212(,,,),(,,,)T n n y y y y A p p p == ，则上式等价于10nj j j Ay p y ==≠∑． （7.5.4）而（7.5.4）式成立的充要条件是向量组12,,,n p p p 线性无关，这又等价于rank A n =．从而T A A 为正定矩阵的充要条件是rank A n =，换句话说，()s x 为正定二次函数的充要条件是rank A n =．定理7.5.2 若rank A n =，则1()T T x A A A b -= （7.5.5）为问题（7.5.2）的唯一最优解．显然，方程组Ax b =有解的充分必要条件是问题（7.5.2）的最优值min ()0s x =．例5 求解线性最小二乘问题 2min Ax b -， （7.5.6）其中31223,3141A b ⎛⎫⎛⎫ ⎪ ⎪=-=- ⎪ ⎪ ⎪ ⎪--⎝⎭⎝⎭．解 因为11472671,()726714350T T A A A A --⎛⎫⎛⎫== ⎪ ⎪-⎝⎭⎝⎭，所以，由公式（7.5.5），求得问题（7.5.6）的最优解22673213/14137141343/103501x ⎛⎫-⎛⎫⎛⎫⎛⎫ ⎪=-= ⎪⎪ ⎪ ⎪-⎝⎭⎝⎭⎝⎭ ⎪-⎝⎭． 由于问题（7.5.6）的最优值为()()11.454285710T Ax b Ax b --=≠．因此，方程组12121232,233,41x x x x x x +=⎧⎪-=-⎨⎪-+=-⎩无解．二、Gauss-Newton 法现在讨论非线性最小二乘问题2m i n ()()s x f x =， （7.5.7） 其中12,()((),(),,()),()n T m i x R f x f x f x f x f x ∈= 不全是线性函数，且假定()i f x 具有一阶连续偏导数，1,2,,i m = ．算法7-9（阻尼Gauss-Newton 法）Step1 选取初始数据．选取初始点0x ，给定允许误差0ε>，令0k =． Step2 检查是否满足终止准则．计算()k f x 及Jacobi 矩阵()k f x ∇，若()()T k k f x f x ε∇<，迭代终止，k x 为问题（7.5.7）的近似极小点；否则，转Step3．Step3 构造Gauss-Newton 方向．计算1[()()]T k k f x f x -∇∇，取1[()()]()()T T k k k k k d f x f x f x f x -=-∇∇∇．Step4 进行一维搜索．求k λ和1k x +，使得()min ()k k k k k s x d s x d λλλ≥+=+， 1k k k k x x d λ+=+．令:1k k =+，返回Step2．同阻尼Newton 法一样，阻尼Gauss-Newton 法具有全局收敛性．三、Levenberg-Marquardt 法算法7-10（LM 法）Step1 选取初始数据．选取初始点0x ，给定初始参数00α>，放大因子1β>及允许误差0ε>，令0k =．Step2 求目标函数值．计算()k f x 及()k s x ． Step3 求Jacobi 矩阵．计算()k f x ∇． Step4 检查是否满足终止准则．若()()T k k f x f x ε∇<，迭代终止，k x 为问题（7.5.7）的近似最优解；否则，转Step5．Step5 构造LM 方向．计算1[()()]T k k k n f x f x I α-∇∇+，令1[()()]()()T T k k k k n k k d f x f x I f x f x α-=-∇∇+∇．Step6 检查目标函数是否下降．计算()k k f x d +和()k k s x d +， 若()()k k k s x d s x +<，转Step8；否则，转Step7．Step7 放大参数．令:k k αβα=，返回Step5．Step8 缩小参数．令11,/,:1k k k k k x x d k k ααβ++=+==+，返回Step2． 初始参数0α和放大因子β应取适当数值，例如，根据经验可取00.01,10αβ==．可以证明，LM 法具有全局收敛性．。
无约束优化与约束优化
淮海工学院实验报告书课程名称：数学实验实验名称：无约束优化与约束优化班级数学091姓名：耿萍学号：090911107 日期：2012.5.18 地点数学实验室指导教师：曹卫平成绩：数理科学系1.实验目的：（1）学习建立最优化数学模型，并通过MATLAB求解。
（2）掌握MATLAB优化工具箱的基本用法，对不同算法进行初步分析、比较。
（3）练习实际问题的非线性最小二乘拟合。
（4）掌握用MATLAB优化工具包解线性规划和非线性规划。
（5）练习建立实际问题的线性规划和非线性规划模型。
2.实验内容：（1）某工厂生产甲乙两种产品，一件甲用A原料1公斤，B原料5公斤；一件已用A原料2公斤，B原料4公斤。
现有A原料20公斤，B 原料70公斤。
甲乙产品每件售价分别为20元和30元。
问如何安排生产使收入最大？当A原料或B原料，以每公斤2元买进B原料。
问如何安排计划使收入最大？当A原料或B原料增加1公斤时对收入影响多大。
又若可用有限资金以每公斤6元买进A原料，以每公斤2元买进原料。
问如何安排计划使利润最大？解释得到的结果。
（2）炼油厂将A,B,C三种原油加工成甲乙丙三种汽油。
一桶原油加工成一桶汽油的费用为4元，每天至多能加工汽油14000桶。
原油的买入价、买入量、辛烷值、硫含量，及汽油的卖入价、需求量、辛烷值、硫含量由下表给出。
问如何安排生产计划，在满足需求的条件下使利润最大？一般来说，做广告可以增加销售，估计一天向一种汽油投入一元广告费，可以使这种汽油日销售量增加10桶。
问如何安排生产和广告计划使利润最大化？原油类别买入价（元/桶）买入量（桶/天）辛烷值（%）硫含量（%）A 45 <=5000 12 0.5B 35 <=5000 6 2.0C 25 <=5000 8 3.0汽油类别卖出价（元/桶）需求量（桶/天）辛烷值（%）硫含量（%）甲70 3000 >=10 <=1.0乙60 2000 >=8 <=2.0丙50 1000 >=6 <=1.0（3）根据实验1表2中给出的美国人口数据，用非线性最小二乘法拟合阻滞增长模型中的参数，注意只拟合2个参数r，xm，与拟合3个参数r，xm，x0有何区别，并与实验1用线性最小二乘法拟合的结果进行比较。
数学实验 7：无约束优化
实验 7：无约束优化习题2：取不同的初值计算下列非线性规划，尽可能求出所有局部极小点，进而找出全局极小点，并对不同算法的结果进行分析比较：c a a c a a x x x x z T T 222111)()(1)()(1min +---+---=，2R ∈x其中)73.0,7.0(),(21=c c ，T a )4,4(1=，Ta )8.3,5.2(2=。
1． 程序设计1)构建题目中的方程function [f,g]=func(x) a1=[4,4]'; a2=[2.5,3.8]'; c=[0.7,0.73];f=-1/((x-a1)'*(x-a1)+c(1))-1/((x-a2)'*(x-a2)+c(2));%自定义目标函数的梯度函数 if nargout>1g(1)=2*(x(1)-a1(1))/((x(1)-a1(1))^2+(x(2)-a1(2))^2+c(1))^2+2*(x(1)-a2(1))/((x(1)-a2(1))^2+(x(2)-a2(2))^2+c(2))^2g(2)=2*(x(2)-a1(2))/((x(1)-a1(1))^2+(x(2)-a1(2))^2+c(1))^2+2*(x(2)-a2(2))/((x(1)-a2(1))^2+(x(2)-a2(2))^2+c(2))^2 end2)进行非线性规划并比较不同算法的不同opt1=optimset('LargeScale','off','MaxFunEvals',1000,'TolFun',1e-8,'To lX',1e-8,’GradObj ’,’off ’);%关闭大规模算法，规定最大调用次数，以及计算精度，拟牛顿法DFGS 公式,是否采用自设梯度%函数可调opt2=optimset(opt1,'HessUpdate','dfp');%拟牛顿法BFP 公式 opt3=optimset(opt1,'HessUpdate','gillmurray'); %拟牛顿法GM 公式 opt4=optimset(opt1,'HessUpdate','steepdesc'); %最速下降法opt5=optimset(opt1,'lineSearchType',' cubicpoly');%采用三次多项式插值 opt6=optimset(opt5,'HessUpdate','dfp'); opt7=optimset(opt5,'HessUpdate','gillmurray'); opt8=optimset(opt5,'HessUpdate','steepdesc');x0=[3,3]';%设置不同初值，比较结果[x1,v1,exit1,out1]=fminunc(@func,x0,opt1) %进行非线性规划[x2,v2,exit2,out2]=fminunc(@func,x0,opt2)[x3,v3,exit3,out3]=fminunc(@func,x0,opt3)[x4,v4,exit4,out4]=fminunc(@func,x0,opt4)[x5,v5,exit5,out5]=fminunc(@func,x0,opt5)[x6,v6,exit6,out6]=fminunc(@func,x0,opt6)[x7,v7,exit7,out7]=fminunc(@func,x0,opt7)[x8,v8,exit8,out8]=fminunc(@func,x0,opt8)3)绘制原函数的三维图形，找到全局极小点n1=2;n2=5; %绘图范围m=300; %取点个数x1=linspace(n1,n2,m);x2=linspace(n1,n2,m);[X,Y]=meshgrid(x1,x2);F=0;for i=1:mfor j=1:mF(i,j)=func([x1(i),x2(j)]');endendmesh(X,Y,F)xlabel('X1');ylabel('X2');zlabel('Func');pause;contour(x1,x2,F,50); %绘制等高线50条xlabel('X1');ylabel('X2');2．运行结果及分析可以看到采用BFGS与GM算法的结果完全相同，事实上，在运行时MatLab会输出警告，GM方法采用的应该还是默认的BFGS算法。
大学数学实验七_无约束优化1
0.4228 0.0000 1.1436e-011
表 1.2 分析方法计算梯度得到的结果
迭代 次数
4 4 18 4 4 18
目标函数引 用次数
6 6 37 6 6 37
比较分析方法计算梯度得到的结果和用数值方法计算梯度得到的结果，发现都能比较准 确的得到最优解 x2=0。但是，用分析方法计算梯度时，目标函数引用次数明显减小（减小 到原来的 1/3 左右），说明用分析方法计算梯度比较节省资源。
12
输出的图像如下。
图 1.5 局部极小点分布
从上图中可以看出，局部极小点集中在
或
或
三条直线上，以及第二
象限中的一条形似双曲线一支的一条曲线上。放大上图，在这条曲线上随机取 10 个点，输
出它们的最优值以判断它们是不是全局最小解。
x1
x2
f(x)
1
-0.9519
0.0689
8.1051e-010
2
三、对不同算法的结果进行分析比较（分析方法计算梯度） 首先计算 f(x)的梯度。
梯度 ：
而
修改函数文件 fun21，将上面计算的梯度增加进去。
8
将 GradObj 参数设置为 on，具体程序与之前所示一样，修改部分如下。
程序运行后得到的结果整理在下表中。
情况
1 2 3 4 5 6
搜索方向
BFGS DFP 最速下降 BFGS DFP 最速下降
最优解 步长搜索
x1
最优解 x2
最优值
混合二三 次插值
0.4196 0.4197 0.4228
0.0000 0.0000 0.0000
5.8944e-015 5.5970e-015 1.1436e-011
无约束优化
第八节 Powell法（方向加速法）
Powell法是利用共轭方向能够加速收敛旳性质 所形成旳一种搜索算法。
一、共轭方向旳生成
二、基本算法
三、改善旳算法
在鲍维尔基本算法中，每一轮迭代都用连结始点和终点 所产生出旳搜索方向去替代原来向量组中旳第一种向量， 而不论它旳“好坏”。
gk1 gk 与 d k 旳共轭方向 d j 正交。
图4-9 共轭梯度法旳几何阐明
第六节变尺度法
变尺度法旳基本思想：
前面讨论旳梯度法和牛顿法，它们旳迭代公式能够看作下列 公式旳特例。
xk1 xk k Hf xk
变尺度法是对牛顿法旳修正，它不是计算二阶导数旳矩阵和 它旳逆矩阵，而是设法构造一种对称正定矩阵H来替代Hesse 矩阵旳逆矩阵。并在迭代过程中，使其逐渐逼近H-1 。
优化设计追求目旳函数值最小，若搜索方向取该点旳负梯度 方向，使函数值在该点附近旳范围内下降最快。
按此规律不断走步，形成下列迭代算法：
xk1 xk akf xk
以负梯度方向为搜索方向，所以称最速下降法或梯度法。
搜索方向拟定为负梯度方向，还需拟定步长因子ak
即求一维搜索旳最佳步长，既有
f
xk1 0
f xk 2 f xk xk1 xk 0
xk1 xk 2 f xk 1 f xk
这是多元函数求极值旳牛顿法迭代公式。
例4-2 用牛顿法求 f x x12 25x xk k 2 f
数值法
能够处理复杂函数及没有数学体现式 旳优化设计问题
xk1 xk ak d k
搜索方向问题是无约束优化措施旳关键。
实验七(无约束优化)
一、实验目的1、掌握MATLAB优化工具箱的基本用法，对不同算法进行初步分析、比较。
2、练习用无约束优化方法建立和求解实际问题的模型（包括最小二乘拟合）。
二、实验容项目一：某分子由25个原子组成，并且已经通过实验测量得到了其中某些原子队之间的距离（假设在平面结构上讨论），如表所示。
请你确定每个原子的位置关系。
原子对距离原子对距离原子对距离原子对距离(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.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.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问题分析：每个原子的位置都是未知的，在坐标系中只有相对的位置参数，不妨固定原子1的坐标为（0，0）。
实验7-无约束优化
'---case6:gillmurray,3rdpoly-----'
fopt=optimset(opt2,'HessUpdate','steepdesc');
[x6,v6,exit6,out6]=fminunc('fun72',x0,fopt)
2.看z坐标可知此时取的自变量范围太宽，应缩小范围寻找最小值
取[x,y]=meshgrid(-0.1:0.01:1.1,-0.1:0.01:1.1)（x为以上x1,y为以上x2）
无论从图像上看还是从表达式上看，都很容易得知全局极小点不可能出现在x,y很大的位置。
由等高线可以看出最小值在x1=1或0或x2=0的附近。
opt2=optimset(opt1,'LineSearchType','cubicpoly');
[x4,v4,exit4,out4]=fminunc('fun72',x0,opt2)
'---case5:dfp,3rdpoly----'
fopt=optimset(opt2,'HessUpdate','dfp');
3.尝试取不同初值找局部极小点，取x0=[-0.5:0.05:1.5,-0.5:0.05:1.5]
看这些不同初值对应的局部极小点位置特点
4.对不同算法做对比
5.用分析方法计算梯度再求解
【结果分析1】
1. 直观认识函数，三维视图
由图已经看不出最优解的位置，因为范围取太大而使图像太扁，下面缩小了两自变量范围，再作图得：
清华大学数学实验_实验7无约束优化1
实验7无约束优化生医0 王言2010013212【实验目的】1．掌握用MATLAB优化工具箱的基本用法，对不同算法进行初步分析、比较。
2．练习用无约束优化方法建立和求解实际问题模型（包括最小二乘拟合）。
【实验内容】5. 某分子由25个原子组成，并且已经通过实验测量得到了其中某些原子对之间的距离（假设在平面结构上讨论），如表7.8所示。
请你确定每个原子的位置关系。
表7.8解：分析与建模： 不妨设第i 点坐标为，其中第一个点坐标为，问题即为求达到最小值的解，则问题转化为无约束优化：，其中：Matlab 代码如下：M 文件：function f=location(x,d);x(1)=0; x(26)=0;f(1)=(x(4))^2+(x(29))^2-d(1)^2; f(2)=(x(12))^2+(x(37))^2-d(2)^2; f(3)=(x(13))^2+(x(38))^2-d(3)^2; f(4)=(x(17))^2+(x(42))^2-d(4)^2; f(5)=(x(21))^2+(x(36))^2-d(5)^2;f(6)=(x(5)-x(2))^2+(x(30)-x(27))^2-d(6)^2; f(7)=(x(16)-x(2))^2+(x(41)-x(27))^2-d(7)^2; f(8)=(x(17)-x(2))^2+(x(42)-x(27))^2-d(8)^2; f(9)=(x(25)-x(2))^2+(x(50)-x(27))^2-d(9)^2; f(10)=(x(5)-x(3))^2+(x(30)-x(28))^2-d(10)^2; f(11)=(x(20)-x(3))^2+(x(45)-x(28))^2-d(11)^2; f(12)=(x(21)-x(3))^2+(x(46)-x(28))^2-d(12)^2; f(13)=(x(24)-x(3))^2+(x(49)-x(28))^2-d(13)^2;),(i i y x )0,0(),(11=y x ∑--+-ji ij j i j id y y x x,2222])()[(∑--+-ji ij j i j i d y y x x ,2222])()[(min ⎥⎦⎤⎢⎣⎡=252432252432y y y y x x x x xf(14)=(x(5)-x(4))^2+(x(30)-x(29))^2-d(14)^2;f(15)=(x(12)-x(4))^2+(x(37)-x(29))^2-d(15)^2;f(16)=(x(24)-x(4))^2+(x(49)-x(29))^2-d(16)^2;f(17)=(x(8)-x(6))^2+(x(33)-x(31))^2-d(17)^2;f(18)=(x(13)-x(6))^2+(x(38)-x(31))^2-d(18)^2;f(19)=(x(19)-x(6))^2+(x(44)-x(31))^2-d(19)^2;f(20)=(x(25)-x(6))^2+(x(50)-x(31))^2-d(20)^2;f(21)=(x(8)-x(7))^2+(x(33)-x(32))^2-d(21)^2;f(22)=(x(14)-x(7))^2+(x(39)-x(32))^2-d(22)^2;f(23)=(x(16)-x(7))^2+(x(41)-x(32))^2-d(23)^2;f(24)=(x(20)-x(7))^2+(x(45)-x(32))^2-d(24)^2;f(25)=(x(21)-x(7))^2+(x(46)-x(32))^2-d(25)^2;f(26)=(x(14)-x(8))^2+(x(39)-x(33))^2-d(26)^2;f(27)=(x(18)-x(8))^2+(x(43)-x(33))^2-d(27)^2;f(28)=(x(13)-x(9))^2+(x(38)-x(34))^2-d(28)^2;f(29)=(x(15)-x(9))^2+(x(40)-x(34))^2-d(29)^2;f(30)=(x(22)-x(9))^2+(x(47)-x(34))^2-d(30)^2;f(31)=(x(11)-x(10))^2+(x(36)-x(35))^2-d(31)^2;f(32)=(x(13)-x(10))^2+(x(38)-x(35))^2-d(32)^2;f(33)=(x(19)-x(10))^2+(x(44)-x(35))^2-d(33)^2;f(34)=(x(20)-x(10))^2+(x(45)-x(35))^2-d(34)^2;f(35)=(x(22)-x(10))^2+(x(47)-x(35))^2-d(35)^2;f(36)=(x(18)-x(11))^2+(x(43)-x(36))^2-d(36)^2;f(37)=(x(25)-x(11))^2+(x(50)-x(36))^2-d(37)^2;f(38)=(x(15)-x(12))^2+(x(40)-x(37))^2-d(38)^2;f(39)=(x(17)-x(12))^2+(x(42)-x(37))^2-d(39)^2;f(40)=(x(15)-x(13))^2+(x(40)-x(38))^2-d(40)^2;f(41)=(x(19)-x(13))^2+(x(44)-x(38))^2-d(41)^2;f(42)=(x(15)-x(14))^2+(x(40)-x(39))^2-d(42)^2;f(43)=(x(16)-x(14))^2+(x(41)-x(39))^2-d(43)^2;f(44)=(x(20)-x(16))^2+(x(45)-x(41))^2-d(44)^2;f(45)=(x(23)-x(16))^2+(x(48)-x(41))^2-d(45)^2;f(46)=(x(18)-x(17))^2+(x(43)-x(42))^2-d(46)^2;f(47)=(x(19)-x(17))^2+(x(44)-x(42))^2-d(47)^2;f(48)=(x(20)-x(19))^2+(x(45)-x(44))^2-d(48)^2;f(49)=(x(23)-x(19))^2+(x(48)-x(44))^2-d(49)^2;f(50)=(x(24)-x(19))^2+(x(49)-x(44))^2-d(50)^2;f(51)=(x(23)-x(21))^2+(x(48)-x(46))^2-d(51)^2;f(52)=(x(23)-x(22))^2+(x(48)-x(47))^2-d(52)^2;end主文件：d=[0.9607 0.4399 0.8143 1.3765 1.2722 0.5294 0.6144 0.3766 0.6893 0.94880.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.92540.6401 0.2467 0.4727 1.3840 0.4366 1.0307 1.3904 0.5725 0.7660 0.4394 1.09521.0422 1.8255 1.4325 1.0851 0.4995 1.2277 1.1271 0.7060 0.8052]';x0=[zeros(1,1),ones(1,24),zeros(1,1),ones(1,24)];[x,norm1]=lsqnonlin(@location,x0,[],[],[],d);xnorm1运行结果x =Columns 1 through 100 0.1925 1.3539 0.8210 0.6910 1.11740.7470 0.9049 0.4348 1.0297Columns 11 through 200.4536 -0.4604 0.7045 0.5812 0.3521 0.12710.4361 1.5726 1.3993 0.9407Columns 21 through 300.2642 1.5103 0.9707 0.3404 0.8626 00.9330 1.6449 0.5483 0.9793Columns 31 through 401.3858 1.2540 0.9459 0.6425 1.2754 1.24390.2204 0.4145 1.2862 0.8651Columns 41 through 500.3007 1.2927 0.4275 0.7643 0.9601 1.83231.3226 1.9180 1.0981 1.1166norm1 =0.0316改变初值：x0=[zeros(1,2),ones(1,22),zeros(1,2),ones(1,22)];结果是：x =Columns 1 through 100 0.8285 0.5479 0.8712 0.8802 1.5054 0.7397 1.0346 0.0074 1.3647Columns 11 through 201.0327 -0.4290 0.4778 0.6123 0.1921 1.4165 0.6782 1.3906 1.2647 1.0951Columns 21 through 300.0732 1.0151 0.1085 0.7128 1.4644 1.00001.5112 0.0944 0.5051 0.9830Columns 31 through 400.8822 0.7628 0.6989 1.0033 0.7893 1.2701 0.3074 0.6494 0.9350 1.1299Columns 41 through 501.6768 1.1750 -0.0659 0.2455 0.6879 1.1010 0.3047 0.4028 1.2196 1.2407Columns 51 through 521.0000 1.0000norm1 =0.2330这一组norm1较大，因此第一组初值较好。
常用的无约束优化方法
一系列沿坐标轴方向的一维搜索问题来求解。在每一 次迭代中，只改变 n 个变量中的一个，其余变量固定 不动，因此常称为单变量法或变量交错法或降维法
#
4.1 坐标轮换法-迭代步长的确定
(1) 最优步长
在沿坐标轴方向的搜索中，利用一维优化方法来确定沿该方向 上具有最小目标函数值的步长，即：
基本方向组为：
S1(
k
)
,S(k 2) Nhomakorabea,
各方向最优步长： a1(k ) , a2(k ) ,
,
S
(k n
)
, an(k )
新生方向：
Sk
x(k) n
x(k) 0
a S (k ) (k ) 11
a2(k ) S2(k )
an(k ) Sn(k )
若在优化搜索过程中出现1(k) =0(或近似等于0)，则方向 Sk
F2 < F3 F2 F3
其中： △ 是在第k 环方向组中，依次沿各方向优化搜索函数值下
降量的最大值，即Sm(k) 方向函数下降量最大
#
鲍威尔修正算法的方向淘汰
x
(k 2
)
S (k) 3
S
( 2
k
)
x1(k )
S1(k )
x(k) m
x(k) m 1
函数下降量△
S (k) n
x(k) 0
(F1)
S (k)
x(k) n
(F2)
第四章 常用的无约束优化 方法
王桂从
无约束优化问题的数学模型
min F(x)
x [x1 x2 xn ] Rn
求上述问题最优解(x*, F*)的方法称为无约束优化方法 无约束优化方法理论研究开展的比较早，构成的优 化方法已很多，也比较成熟。使用无约束优化方法 ，不仅可以直接求无约束优化设计问题的最优解， 而且通过对无约束优化方法的研究给约束优化方法 建立明确的概念及提供良好的基础，某些优化设计 方法就是先把优化设计问题转化为无约束问题后， 再直接用无约束优化方法求解。
无约束优化
牛顿法
牛顿法的基本思想就是利用二次函数来代替原目标函数， 以二次函数的极小点作为原目 标函数的极小点的近似，并逐步逼近该点。 设一般目标函数 f ( x ) 具有一、 二阶偏导数， 此时， 在 X k 处做泰勒展开并取值二次项，
得：
1 f ( X ) ( X ) f ( X k ) f ( X k )T ( X X k ) ( X X k )T 2 f ( X k )T ( X X k ) 2
为搜索方向，式中：
F ( X k ) F ( X ) , x1
k
F ( X k ) F ( X k ) ,.... x2 xn
T
第 k 1 次得到的新点：
X k 1 X k k F ( X k )
一般步长 1 常采用沿负梯度方向做一维搜索：
习 题 一 ： 试 用 梯 度 法 求 目 标 函 数 F ( X ) ( x1 1) ( x2 1) 的 极 小 值 ， 设 初 始 点
2
2
X 0 [0, 0]T ， 0.01 。
习题二：试用梯度法求目标函数 F ( X ) x1 25 x2 的最优解。初始点 X 0 [2, 2]T ，迭代精 度 0.05 。
其步
阻尼牛顿法一般步骤： 1 2 选取初始数据。选取初始点X 0 , 给定允许误差 0，令k 0 检验是否满足终止准则。计算f ( X k ), 若 f ( X k ) , 迭代终 止，X k 为问题的近似最优解；否则，转Step3. 3
2 k 构造牛顿方向。计算 f ( X ) ，取 1
0 0.130769
T X1 （ 11 ， ） 0.13.769(8, 2)T ( 0.046152, 0.738462)T 第二次迭代：
08实验7无约束优化2
f(x)=0, x=fsolve (‘f’,x0) xRn x=fminunc(‘f’,x0) x=lsqnonlin (‘f’,x0) X=fminsearch(‘f’,x0) x=lsqcurvefit(‘f’,x0)
1、fminbnd
fminunc fminsearch
调用格式为： [x, fv, fe, out]=fminbnd(@f, x0, opt, p1, p2,….) 例1、求解
问 题 模 型 基本命令的用法 .M文件
min f(x), xR
无约束极小 (非线性规划)
x=fmin(‘f’,x0)
function f=f(x) function f=f(x)
function r=F(x) function yy=F(x,t)
min f(x), xRn x=fminu(‘f’,x0) Ｘ=fmins(‘f’,x0)
t c(t) 0.25 19.21 0.5 18.15 1 15.36 1.5 14.10 2 12.89 3 9.32 4 7.45 6 5.24 8 3.01
t=[.25 .5 1 1.5 2 3 4 6 8]; c=[19.21 18.15 15.36 14.1 12.89 9.32 7.45 5.24 3.01];
f ( x k ) ( f x'1 , f x'2 , , f x'n )T
f (1x* ) 0 Gk 1d k 1 f ( xk )或xk H k 1f k
2.2 搜索步长的确定 一维优化问题
min f ( x k a d k )
a
其极小值点应为方程 (d k )T f ( x k a d k ) 0 的精确解. 搜索方法有:
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
迭代次数：4 目标函数调用次数：18
5、搜索方向：DFP；搜索步长：三次插值
输出结果如下 6
x5 =
0.4185 0.0000
v5 =
5.9198e-015
exit5 =
1
out5 =
iterations: 4
funcCount: 18
stepsize: 1
firstorderopt: 2.2075e-008
funcCount: 18
stepsize: 1
firstorderopt: 2.2209e-008
algorithm:
'medium-scale:
Quasi-Newton line search'
message: [1x468 char]
最优解 x1=0.4185 最优解 x2=0.0000
最优值=5.9926e-015
输出结果如下。
x6 =
0.4189 0.0000
v6 =
9.6391e-012
exit6 =
1
out6 =
iterations: 18
funcCount: 111
stepsize: 10
firstorderopt: 8.7861e-007
algorithm:
'medium-scale:
Quasi-Newton line search'
i
i
i
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
大学数学实验七 无约束优化 实验报告
【实验目的】 1、掌握 MATLAB 优化工具箱的基本用法，对不同算法进行初步分析、比较。 2、练习用无约束优化方法建立和求解实际问题的模型（包括非线性最小二
乘拟合）。 【实验内容】
2 取不同的初值计算下列非线性规划，尽可能求出所有局部极小点，进而找出全局极小点， 并对不同算法（搜索方向、搜索步长、数值梯度与分析梯度等）的结果进行分析、比较。 （1）
-0.8468
0.08596
1.9026e-010
3
-0.7043
0.1185
2.9597e-009
4
-0.5373
0.179
3.4633e-009
5
-0.4258
0.242
5.7723e-010
6
-0.3158
0.3336
3.6726e-011
7
-0.2145
0.4596
1.0024e-010
8
-0.1114
上面的程序中，x 和 y 分别表示 和 ，z 表示 f(x1, x2)的值。自变量的范围均取在[-2, 2] 上，根据输出的图像进一步缩小范围。
输出的三维图像如下。
1
图 1.1 三维图像
图 1.2 等高线
第一幅图中，最大的函数值的数量级已到 107，说明画图范围取得过大，函数值为 0 的 电在该图中均呈现一平面，故需要缩小自变量的范围再绘图。修改程序如下。
ቤተ መጻሕፍቲ ባይዱ
表 1.1 数值方法计算梯度得到的结果
迭代 次数
4 4 18 4 4 18
目标函数引 用次数
18 18 111 18 18 111
从上表可以看出，无论选择什么搜索方向和搜索步长，都能得到比较准确的结果（x2=0）。 但是从最优值、迭代次数和目标函数引用次数 4 个参数来看，BFGS 和 DFP 明显比最速下 降得到的结果要好。搜索步长的选择对这个问题没有影响。
message: [1x468 char]
最优解 x1=0.4189 最优解 x2=0.0000
最优值=9.6391e-012
1 表示收敛
迭代次数：18 目标函数调用次数：111
7
程序运行结果整理在下表中。
情况
1 2 3 4 5 6
搜索方向
BFGS DFP 最速下降 BFGS DFP 最速下降
最优解 步长搜索
v3 =
9.6391e-012
exit3 =
1
out3 =
iterations: 18
funcCount: 111
stepsize: 10
firstorderopt: 8.7861e-007
algorithm:
'medium-scale:
Quasi-Newton line search'
message: [1x468 char]
Quasi-Newton line search' message: [1x468 char]
最优解 x1=0.4185 最优解 x2=0.0000
最优值=5.9198e-015
1 表示收敛
迭代次数：4 目标函数调用次数：18
3、搜索方向：最速下降法；搜索步长：混合二三次插值
输出结果如下。
x3 =
0.4189 0.0000
10
牛顿法的迭代公式如下：
初值依旧取
。
编写计算程序如下。
11
这个程序中设置的停止迭代的条件也是
，认为此时得到的解是局部
极小点。输出的结果是(1.0000, 0.0000)。虽然根据理论分析可知这点是全局极小点，可是由
于初值给的是(0.4, 0.2)，输出的结果和用其他方法得到的结果非常不同，可能是由于 f(x)的
最优解 x1=0.4185 最优解 x2=0.0000
最优值=5.9926e-015
1 表示收敛
迭代次数：4 目标函数调用次数：18
2、搜索方向：DFP；搜索步长：混合二三次插值
4
输出结果如下。
x2 = 0.4185 0.0000
v2 = 5.9198e-015
exit2 = 1
out2 = iterations: 4 funcCount: 18 stepsize: 1 firstorderopt: 2.2075e-008 algorithm: 'medium-scale:
x1
最优解 x2
最优值
混合二三 次插值
0.4185 0.4185 0.4189
0.0000 0.0000 0.0000
5.9926e-015 5.9198e-015 9.6391e-012
0.4185 三次插值 0.4185
0.0000 5.9926e-015 0.0000 5.9198e-015
0.4189 0.0000 9.6391e-012
输出结果如下。
x1 = 0.4185 0.0000
v1 = 5.9926e-015
exit1 = 1
out1 = iterations: 4 funcCount: 18 stepsize: 1 firstorderopt: 2.2209e-008 algorithm: 'medium-scale:
Quasi-Newton line search' message: [1x468 char]
设
函数 是一个二变量的、乘积形式的函数，而且可以预见，其梯度向量和 Hessian 矩
阵的表达式比较复杂，故若通过求梯度的方法求解本题会十分麻烦。
从这个函数表达式可以看出
恒
故若某个 能使
，应该就是最优解。从表达式容易观察得，当
或
或
时，
。
下面用 MATLAB 求解这个问题。 一、输出三维图像和等高线，直观观察最优解所在位置 首先，先让 MATLAB 输出 的三维图像，直观地观察最优解的大致范围。程序如下。
四、自编程序实现用最速下降法和牛顿法进行计算
最速下降法的迭代公式如下：
初值依旧取
。之前也已经求出了 。编写计算程序如下。
9
上面的程序中设置的停止迭代的条件是
，认为此时得到的解是局部
极小点。输出的结果是(0.4194, 0.0000)，和用 MATLAB 自带命令求解出来的结果很接近。
观察其收敛过程（部分）：
利用之前绘三维图像和等高线的程序，在
和
上画出三
维图像和等高线。
输出结果如下。
图 1.6 部分局部极小点分布三维图
图 1.7 部分局部极小点分布等高线
从上图可以看出，在第二象限的确存在一条由局部极小点构成的线。根据 f(x)的表达式： 当
时， 也=0。故上式也是本题的一个最优解。 14
下面在图 1.5 的基础上，绘出 输出图像如下。
0.6555
1.2704e-010
9
-0.04218
0.8477
1.3768e-012
10
-0.02059
0.9218
3.6449e-012
表 1.3 第二象限中形似双曲线一支的曲线上的局部最小点对应的函数值
13
观察上表数据，发现这些点的最优解都已经十分靠近 0，很有可能也是全局极小点。需
要利用 MATLAB 进行进一步判断。
0.4228 0.0000 1.1436e-011
表 1.2 分析方法计算梯度得到的结果
迭代 次数