无约束最优化问题的基本研究毕业论文

1、 . . . 关于无约束最优化问题求解的基本研究关于无约束最优化问题求解的基本研究 . . . 摘要摘要无约束最优化计算方法是数值计算领域中十分活跃的研究课题之一，快速的求解无约束最优化问题，除了自身的重要性以外，还体现在它也构成一些约束最优化问题的子问题.因此，对于无约束最优化问题，如何快速有效的求解一直是优化工作者十分关心的事.论文研究求解无约束最优化问题的几种主要的导数法，并且讨论了这些方法的优缺点以与每种方法的适用围.同事论文分别对每种方法给出了具体实例，并对例子进行了 matlab 软件实现关键词：关键词：无约束最优化; 导数法 ;极值 ; 精确度 AbstractAbstractU

2、nconstrained optimization numerical calculation method is very active in the field of research, one of the most rapidly solving unconstrained optimization problems, in addition to its importance, is also reflected in some of the constraints that it also constitutes a sub-problem of optimization prob

3、lems. Therefore, for unconstrained optimization problems, how fast and effective solution has been optimized workers very concerned about. Thesis for solving unconstrained optimization problems several major derivative method, and discusses the advantages and disadvantages of these methods as well a

4、s the scope of application of each method. Colleagues papers for each method were specific examples are given, and examples of the matlab softwareKeyword：Unconstrained optimization Derivative method Extremum Accuracy目录目录摘要摘要- - 1 1 - -ABSTRACTABSTRACT- - 2 2 - -第一章绪论第一章绪论- - 4 4 - -1.11.1 研究背景与意义研究背

5、景与意义- 4 - . . . - 2 - / 391.21.2 问题阐述与简介问题阐述与简介- 4 -第二章无约束问题的极值条件第二章无约束问题的极值条件- - 6 6 - -2.1.2.1. 无约束极值问题无约束极值问题- 6 -2.22.2 必要条件必要条件- 6 -2.32.3 二阶充分条件二阶充分条件- 8 -2.42.4 充要条件充要条件- 8 -第三章求解无约束最优化的几种主要方法第三章求解无约束最优化的几种主要方法- - 1010 - -3.13.1 最速下降法最速下降法- 10 -3.23.2 牛顿法牛顿法- 15 -3.33.3 修正牛顿法修正牛顿法- 19 -3.43.4

6、 共轭梯度法共轭梯度法- 23 -3.53.5 变尺度法变尺度法- 26 -结束语结束语- - 3737 - -参考文献参考文献- - 3838 - -致致- - 3939 - -第一章第一章 绪论绪论1.11.1 研究背景与意义研究背景与意义追求最优化目标是人类共同的理想，最优化就是从众多可能方案中选出最佳方案，以达到最优目标.最优化理论和算法是在第二次世界大战后迅速发展起来的一门新兴的应用数学分支，它是一门应用性很强的年轻学科.虽然最优化可以追朔到很古老的极值问题，但是直到 1947 年 Dantzig 提出一般线性规划问题 . . . - 3 - / 39的单纯形法之后，它才成为一门独立

7、的学科.近三、四十年来随着现代科技的发展和电子计算机的广泛应用，进一步推动了最优化的迅猛发展与其理论和算法的研究.现在最优化理论已广泛应用与生产、管理、军事国防、政府决策、交通运输、经济规划等方面.无约束最优化计算方法不仅本身有着不少实际应用，而且与约束最优化计算方法有着紧密的联系：一方面有些处理无约束最优化问题的方法能直接推广应用于约束最优化问题；另一方面，还可以把一些约束最优化问题转化为无约束最优化问题来处理.因此从这个意义上讲，无约束最优化计算方法也是处理约束最优化问题的基本方法.研究求解无约束最优化问题的有关理论和算法，在近几十年来迅速发展并且日趋成熟.随着计算机的发展和普遍应用，作为

8、一种有效的最优化方法无约束最优化方法在工程设计、管理优化、系统分析等方面的应用日益开拓，愈来愈受到应用部门的重视，所以研究无约束最优化问题的计算方法是意义重大的1.21.2 问题阐述与简介问题阐述与简介 无约束法指寻求 元实函数在整个维向量空间上的最优值点n f Xnn的方法.这类方法的意义在于：虽然实用规划问题大多是有约束的，但许多约束最优化方法可将有约束问题转化为若干无约束问题来求解. 无约束最优化方法大多是逐次一维搜索的迭代算法.这类迭代算法可分为两类.一类不涉与导数，只用到函数值，称为直接法.另一类需要用目标函数的导函数，称为解析法.这些迭代算法的基本思想是：在一个近似点处选定一个有利

9、搜索方向，沿这个方向进行一维寻查，得出新的近似点.然后对新点施行同样手续，如此反复迭代，直到满足预定的精度要求为止.根据搜索方向的取法不同，可以有各种算法.属于直接型的算法有交替方向法（又称坐标轮换法） 、模式搜索法、旋转方向法、鲍威尔共轭方向法和单纯形加速法等.属于解析型的算法有：梯度法：又称最速下降法.这是早期的解析法，收敛速度较慢.牛顿法：收敛速度快，但不稳定，计算也较困难.共轭梯度法：收敛较快，效果较好.变尺度法：这是一类效率较高的方法.其中达维登-弗莱彻-鲍威尔变尺度法，简称 DFP 法，是最常用的方法.本文主要研究无约束最优化问题中主要的几种解析法的算法理论，并对各个方法进行了举例

10、分析和 matlab 软件实现. . . . - 4 - / 39第二章第二章 无约束问题的极值条件无约束问题的极值条件2.1.2.1. 无约束极值问题无约束极值问题考虑非线性规划问题 （2.1） nminRfXX其中是定义在上的实函数，这个问题是求在维欧式空间的极小 f XnR f Xn . . . - 5 - / 39点，称为无约束极值问题，这是一个古典的极值问题.2.22.2 必要条件必要条件为研究函数的极值条件，先介绍一个定理 f X定理定理 2.12.1 设函数 在点可微，如果存在方向，使，则 f XXd 0fX d存在，使得对每个，有.00, ffXdX证明证明 函数在的一阶 Ta

11、ylor 展开式为fXdX =+d+TfffXdXXd （2.2） =+TffdXXd其中当时，.00d由于，当充分小时，在（1.3.2）式中 0fX d d+0,TfXd因此存在，使得时，有00, +0TfdXd从而由（2.2）式得出. ffXdX利用上述定理可以证明局部极小点的一阶必要条件.定理定理 2.22.2 设函数在点可微，若时局部最小点，则梯度 f XXX =fX0证明证明 用反证法，设，令方向，则有 0fX d=fX- 2=0TTffffXdXXX-根据定理 2.1，必存在，使得当时，成立00, . . . - 6 - / 39 ffXdX这与是局部极小点矛盾X下面，利用函数的

12、Hesse 矩阵，给出局部极小点的二阶必要条件 f X定理定理 2.32.3 设函数在点处二次可微，若是局部极小点，则梯度 f XXX，并且 Hesse 矩阵半正定 0fX 2fX证明证明 定理 1.3.2 已经证明，现在只需证明 Hesse 矩阵半正 =fX0 2fX定.设是任意一个维向量，由于在处二次可微，且，则有dn f XX fX0， 2221=+2TfffXdXdX dd经移项整理，得到 （2.3） 222212TfffdXdXdX d由于是局部极小点，当充分小时，必有x ffXdX因此由（2.3）式推得， 20TfdX d即是半正定的. 2fX2.32.3 二阶充分条件二阶充分条件

13、 下面给出局部极小点的二阶充分条件定理定理 2.42.4 设函数在点处可微，若梯度，且 Hesse 矩阵 f XX fX0正定，则是局部极小点. 2fXX . . . - 7 - / 39证明证明 由于在的二阶 Taylor 展开式为 ,fX0 f XX （2.4） 2212TfffXXXXXXXXX设的最小特征值为，由于正定，必有 2fXmin0 2fX 22minTfXXXXXXX从而由（1.3.4）式得出 22min212ffXXXXXXXX当时，因此存在的邻域，当XX220XXXXX,NX时，即是的局部极小点,NXX ffXXX f X2.42.4 充要条件充要条件前面的几个定理分别给

14、出无约束极值的必要条件和充分条件，这些条件都不是充分必要条件，而且利用这些条件只能研究局部极小点.下面在函数凸性的假设下，给出全局极小点的充分必要条件定理定理 2.52.5 设是定义在上的可微凸函数，则为全局极小点 f XnnXX的充分必要条件是梯度 fX0证明证明 必要性是显然的，若是全局极小点，自然是局部极小点，根据定理X1.3.2，必有. fX0 现在证明充分性，设，则对任意的，有 fX0nX，由于是可微的凸函数，则有 0TfXXX f X， =TffffXXXXXX即是全局极小点X 在上述定理中，如果是严格凸函数，则全局极小值是唯一的 f X 上面介绍的几个极值条件，是针对极小化问题给

15、出的，对于极大化问题，可 . . . - 8 - / 39以给出类似的定理例例 2.12.1 利用极值条件解下列问题 22221121min12deffxxxxX先求驻点.由于， 3111422fxxx222fxx令，即 fx03114220 xx，220 x 解此方程组，得到驻点.12,1,0TTx xX再利用极值条件判断是否为极小点，由于目标函数的 Hesse 矩阵X， 221122002xfX由此可知. 210002fX显然为正定矩阵，根据定理 1.3.4，驻点是局部最小点 2fX1,0TX第三章第三章 求解无约束最优化的几种主要方法求解无约束最优化的几种主要方法3.13.1 最速下降法

16、最速下降法 最速下降法又称为梯度法，是 1847 年由著名数学家 Cauchy 给出的，它是解析法中最古老的一种，其他解析方法或是它的变形，或是受它的启发而得到的，因此它是最优化方法的基础.作为一种基本的算法，他在最优化方法中占有重要 . . . - 9 - / 39地位.3.1.13.1.1 最速下降法的算法原理最速下降法的算法原理最速下降法的搜索法向是目标函数的负梯度方向，最速下降法从目标函数的负梯度方向一直前进，直到到达目标函数的最低点.已知目标函数在点的梯度为：( )kX( )( )( )( )12.TkkkknffffxxxXXXX当求目标函数的最小点时，由于函数沿负梯度方向下降最快

17、，故在点的探( )kX索方向应取该点的负梯度方向，即( )( )( )kkkff XSX显然，为单位向量.这样第次迭代计算所得的新点为( )kS1k ( )( )(1)( )( )( )( )( )kkkkkkkkffXXXSXX负梯度仅给出了最优化方向，而没有给出步长的大小，所以可能有各种各样的最速下降的过程，它们依赖于的大小.( )( )kkfX3.1.23.1.2 步长步长的两种取法：的两种取法： k一种方法是任意给定一个初始步长，使满足条件：( )( )( )( )()()kkkkffXSX另外一种方法是沿负梯度方向做一维探索，以求解一维最优化问题的最优步长，即对目标函数极小，以得到最

18、优步长：( )( )( )( )( )0min()()kkkkkffXSXS以此最优步长作为由点出发沿该点的负梯度方向探索的步长.( )kX( )k这种方法的迭代计算的收敛性，可用以下三式中的任一式或二式作为准则来进行判断：( )1( )(1)2( )( )(1)3kkkkkkffffXXXXXX . . . - 10 - / 393.1.33.1.3 最速下降法算法步骤最速下降法算法步骤用最速下降法求无约束多维极值问题的算法步骤如下：min( ),nfX X（1）取初始点，精度，令 0X00k （2） 计算搜索方向，其中表示函数在点( )( )()kkf VX( )()kfX( )f X处的

19、梯度；( )kX（3） 若，则停止计算；否则，从出发，沿进行一维搜索，( )kV( )kX( )kV即求，使得.此处的一维搜索k( )( )( )( )0()min)kkkkkffXV(XV可以用黄金分割法等算法（4） 令，转步骤（2）.(1)( )( ),1kkkkkkXXV 例例 3.13.1 试用最速下降法求目标函数的极小值，设初始点 2212( )4fxxX；收敛要求.(0)22TX2110解：原函数的梯度，在点的梯 1212( )(28TTfffxxxxXX)X(0)X度为.(0)416TfX梯度的模为22(0)(0)(0)2212()()41616.492fffxxXXX梯度的负方

20、向为(0)14160.2430.97016.492TT S(1)1(1)(0)(0)(1)220.24320.24320.97020.970 xx XXS(0)(0)22(20.243 )4(20.970 )fXS(0)(0)2(20.243 )( 0.243)8(20.970 )( 0.970)7.64516.492dfdXS令，求出，算得(0)(0)0dfdXS*2.157 . . . - 11 - / 39(1)(1)1.4760.0923 ,2.9520.738TTfXX梯度的模为22(1)2.9520.7383.043f X根据收敛准则，故未达到要求，应继续探索.(1)213.043

21、10fX下一步探索放向为(1)12.9520.7380.9700.2433.043TT S(2)1(2)(1)(1)(2)21.4760.9701.4760.9700.09230.2430.09230.243xxXXS(1)(1)22(1.4760.970 )4( 0.09230.243 )fXS，得到(1)(1)2.3543.0430dfdXS*1.293(2)(2)0.2220.222 ,0.4441.776TTfXX(2)211.83110fX未达到收敛要求，所以还应继续探索，下一步探索方向为(2)10.4441.7760.2420.9701.831TT S(3)1(3)(2)(2)(3

22、)20.2220.2420.2220.970 xxXXS(2)(2)22(1)(1)(0.2220.242 )4(0.2220.970 )7.6441.8300fdfdXSXS得到：*0.239(3)(3)0.1640.0098 ,0.3280.0784TTfxx(3)210.33710fx继续探索，当探索到点时，(7)0.00160.000096Tx，达到预定的收敛要求，因而可认为为最(7)210.003210fx*(7)xx优点，而为极小值. *226(0.0016)4( 0.000096)2.596 100fx . . . - 12 - / 393.1.43.1.4 最速下降法的最速下降

23、法的 matlabmatlab 实现实现首先建立 M 文件：function x,val,k=grad(fun,gfun,x0)%功能：用最速下降法求解无约束问题： minf（x）%输入：x0是初始点，fun，gfun分别是目标函数和梯度%输出：x,val分别是近似最优点和最优值，k是迭代次数.maxk=5000; %最大迭代次数rho=0.5; sigma=0.4;k=0; eps=1e-7;while(kmaxk) g=feval(gfun,x0);%计算梯度 d=-g; %计算搜索方向 if(norm(d)eps),break;end m=0;mk=0; while(m20) %Armi

24、jo搜索 if(feval(fun,x0+rhom*d)feval(fun,x0)+sigma*rhom*g*d) mk=m;break; end m=m+1; end x0=x0+rhomk*d; k=k+1;endx=x0;val=feval(fun,x0);然后建立fun，gfun的m文件function f=fun(x)f=x(1)2+4*x(2)2;function g=gfun(x)g=2*x(1),8*x(2);求解结果为 x = 0 0val = 0 . . . - 13 - / 39k = 23.1.53.1.5 最速下降法的锯齿现象最速下降法的锯齿现象 最速下降法对初始点的

25、选择要求不高，每一轮迭代工作量较少，它可以很快的从初始点达到极小点附近.但在接近极小点时，最速下降法却会出现锯齿现象，收敛速度很慢，因为对一般二元函数，在极小点附近可用极小点的二阶泰勒多项式来近似，而后者为凸函数时，他的等值线是一族共心椭圆，特别是当椭圆比较扁平时，最速下降法的收敛速度越慢. 至于最速下降法出现锯齿现象的原因，可以作如下粗略解释：由在的泰勒展式： f X kX =kkkfff kXXXP =kkfPX为了出从出发沿方向的极小点，令 kX kP =0,TkkkkfXPP则有. 10kkffXX即方向与方向正交.这表明迭代产生的点列11kkf PX =kkfPX所循路径是“之”字行

26、的.当接近极小点时，每次迭代移动的步长 kX kX*X很小，这样就呈现出锯齿现象，影响了收敛速度.因此常常将梯度法与其他方法结合起来使用.3.1.63.1.6 最速下降法的说明最速下降法的说明最速下降法的优点是算法简单,每次迭代计算量小,占用存量小,即使从一个不好的初始点出发,往往也能收敛到局部极小点.但它有一个严重缺点就是收敛速度慢.沿负梯度方向函数下降很快的说法,容易使人们产生一种错觉,认为这一定是最理想的搜索方向,沿该方向搜索时收敛速度应该很快,然而事实证明,梯度法的收敛速度并不快.特别对等值线(面)有狭长深谷形状的函数,收敛速度更慢.其原因是由于每次迭代后下一次搜索方向总是与前一次搜索

27、方向相互垂直,如此 . . . - 14 - / 39继续下去就产生所谓的锯齿现象.即从直观上看,在远离极小点的地方每次迭代可能使目标函数有较大的下降,但是在接近极小点的地方,由于锯齿现象,从而导致每次迭代行进距离缩短,因而收敛速度不快.3.23.2 牛顿法牛顿法3.2.13.2.1 牛顿法算法思想与原理牛顿法算法思想与原理如前面所提到的，最速下降法在最初几步迭代中函数值下降很快外，总的说来下降的并不快，且愈接近极值点下降的愈慢.因此，应寻找使目标函数下降更快的方法.牛顿法就是一种收敛很快的方法，其基本思路是利用二次曲线来逐点近似原目标函数，以二次曲线的极小值点来近似原目标函数的极小值点并逐渐

28、逼近改点.一维目标函数 在点逼近用的二次曲线（即泰勒二次多项式）为( )f X( )kX( )( )( )( )( )( )21()()()()()()2kkkkkkfffXXXXXXXX此二次函数的极小点可由求得.( )()0kX对于维问题，为目标函数在点逼近用的二次曲线为：nn( )f X( )kX( )( )( )( )( )2( )( )1()()() . .().2kkkkkTkkfff XXXXXXXXXX令式中的，则上式可改写为：Hessian2( )( )()()kkfHXX( )( )( )( )( )( )( )1()()() . .().2( )kkkkkTkkffHf

29、XXXXXXXXXXX当时可求得二次曲线的极值点，且当且仅当改点处的矩阵( )0X( )XHessian为正定时有极小值点.由上式得：( )( )( )( )()()kkkfH XXXXX令，则( )0X( )( )( )()()0kkkfHXXXX若为可逆矩阵，将上式等号两边左乘，则得 ( )()kH X1( )()kHX1( )( )( )()()0kkknHfIXXXX整理后得 . . . - 15 - / 391( )( )( )()()kkkHfXXXX当目标函数是二次函数时，牛顿法变得极为简单、有效，这时是一个()f X( )()kH X常数矩阵，式变成精确( )( )( )( )

30、( )( )( )1()()() . .().2( )kkkkkTkkffHf XXXXXXXXXXX表达式，而利用式作一次迭代计算所得的就是最优点1( )( )( )()()kkkHfXXXXX.在一般情况下不一定是二次函数，则不能一步就能求出极小值，即极小值点不*X( )f X在方向上，但由于在点附近函数与是近似的，所1( )( )()()kkHfXX( )kX( )X( )f X以这个方向可以作为近似方向，可以用式求出点作为1( )( )( )()()kkkHfXXXXX一个逼近点.这时式可改成牛顿法的一般迭代公式：(1)kX1( )( )( )()kkkHfXXX(X1(1)( )(

31、)( )()()kkkkHfXXXX式中称为牛顿方向，通过这种迭代，逐次向极小值点逼近.1( )( )()()kkHfXX*X牛顿法是基于多元函数的泰勒展开而来的，它将作为探索1( )( )()()kkHfXX方向，因此它的迭代公式可直接写出来：1(1)( )( )( )()()kkkkHfXXXX3.2.23.2.2 牛顿法算法步骤牛顿法算法步骤（1） 给定初始点，与精度，令；(0)x00k （2） 若，停止，极小点为，否则转步骤（3） ；( )()kfX( )kX（3） 计算，令；12( )()kfX1( )( )( )()()kkkHf SXX（4） 令，转步骤（2）.(1)( )( )

32、kkkXXS1kk . . . - 16 - / 39例例 3.23.2 试用牛顿法求目标函数的极小点.22121212( )10460fxxx xxxX解解：取，则( )00TkX112( )21221010()244kfxxxffxxxX( )222112( )2( )22221221()()12kkkXXffxx xHfffx xx XX1(1)( )( )( )()()0211081012463kkkkXXH Xf X 即为最小点(1)86k X*X3.2.33.2.3 牛顿法的牛顿法的 matlabmatlab 实现实现首先建立 M 文件 newton.mfunction x,fva

33、l,iterations=newton(fun,x0,tol)%这是一个用Newton法求解无约束优化问题的matlab程序%第一个参数fun是一个包含函数,梯度,hessian矩阵的M文件%比如：% function f,grad,hessian=fun_grad_hess(x)% f=x(1)3+(x(2)-1)2;% grad=3*x(1)2;2*x(2)-2;% hessian=6*x(1) 0;0 2;% end%调用上面函数格式为:x,.=newton_armijo(fun_grad_hess,.)%第二个参数x0是迭代的初始点(列向量)%第三个参数tol是求的结果的精度if na

34、rgin=2 tol=1e-6;elseif nargin=tol p=-inv(H)*grad; x=x+p; fval,grad,H=fun(x); count=count+1;endif nargout=3 iterations=count;end然后编写待求函数和梯度，hessian矩阵组成的M文件fun_grad_hess1.m如下：function f,grad,hessian=fun_grad_hess1(x)f=x(1)2+x(2)2-x(1)*x(2)-10*x(1)-4*x(2)+60;grad=2*x(1)-x(2)-10;2*x(2)-x(1)-4;hessian=2

35、-1;-1 2;end最后在工作窗口调用函数newton.m如下：x,fval,lamla=newton(fun_grad_hess1,0;0,1e-6)所得结果如下：x = 8.0000 6.0000fval = 8.0000lamla = 13.2.43.2.4 牛顿法的说明牛顿法的说明从上例看到,用 Newton 法求解,只经一轮迭代就得到最优解.这并不偶然,由Newton 方向的构造知,对于正定二次函数,Newton 方向就是指向其极小点的方向.因此，用 Newton 法解目标函数为正定二次函数的无约束最优化问题，只需一次迭代即可得到最优解.对于目标函数不是二次函数的无约束最优化问题,

36、一般地说, . . . - 18 - / 39用 Newton 法通过有限轮迭代并不能保证可求得最优解.但因目标函数在最优解附近近似于二次函数,因此当先取接近于最优解的初始点用 Newton 法求解时,其收敛速度一般是快的.可证在初始点离最优解不远的条件下,Newton 法是二次收敛的.但初始点远离最优解时,此法并不一定收敛Newton 法具有二次收敛的优点,但它存在下面四个严重的缺点: 虽每次迭代不会使目标函数 f(X)上升,但不能保证 f(X)下降.当 Hesse 矩阵非正定时,Newton 法的搜索将会失败. 对初始点要求严格.一般要求比较接近或有利于接近极值点,而这在实际计算中是比较难

37、办的. 因搜索方向要求 Hesse 矩阵的逆,在某迭代时可能求不出此方向.若目标函数Hesse 矩阵为奇异,则不有构成牛顿方向,无法迭代牛顿方向构造困难,计算相当复杂,除了求梯度以外还需计算 Hesse 矩阵与其逆矩阵,占用机器存相当大.3.33.3 修正牛顿法修正牛顿法3.3.13.3.1 修正牛顿法的基本原理修正牛顿法的基本原理为了克服 Newton 法的缺点,人们保留选取 Newton 方向作为搜索方向,摒弃其步长恒取 1,而用一维搜索确定最优步长,由此产生的算法称为修正 Newton 法.3.3.23.3.2 修正牛顿法的迭代步骤修正牛顿法的迭代步骤（1） 给定初始点，与精度，令；(0

38、)X00k （2） 若，停止，极小点为，否则转步骤（3） ；( )()kfX( )kX（3） 计算，令；12( )()kfX1( )( )( )()()kkkHf SXX（4） 用一维搜索法求，使得，令( )( )( )( )0()min()kkkkkffXSXS，转步骤（2）.(1)( )( )kkkkXXS1kk例例 3.33.3 用修正牛顿法求解下列问题： ， 22121122min22fxxxx xxX 初始点 ， . 000 X610 . . . - 19 - / 39 解解 计算与 . fX 2fX , , 1212142122xxfxx X 24222fX故有 ， ， 011fX

39、 024222fX 1021122112f X因而牛顿方向为 10002111122311122ff SXX从出发，沿牛顿方向 作一维搜索，令步长变量为，记最优步长为， 0X 0S0则有 ， 001033022 XS故 22002333222225542f XS令 0055022fXS则有 . . . - 20 - / 39，01故 1000110133022 XXS计算： 1fX ， 131412023012122f X故有 ，停止计算，输出，即为极小点. 1fX 1X 1*132XX3.3.3 修正你顿法的 matlab 实现首先建立 m 文件：function x,val,k=revis

40、enm(fun,gfun,Hess,x0)%功能：用修正牛顿法求解无约束问题：minf(x)%输入：x0是初始点，fun,gfun,Hess分别是求% 目标函数值，梯度，Hesse矩阵的函数%输出： x,val分别是近似最优点和最优值，k是迭代次数.n=length(x0);maxk=150;rho=0.55;sigma=0.4;tau=0.0;k=0; epsilion=1e-5;while(kmaxk) gk=feval(gfun,x0); %功能：计算梯度 muk=norm(gk)(1+tau); Gk=feval(Hess,x0); %功能：计算Hess矩阵 Ak=Gk+muk*eye

41、(n); dk=-Akgk;%解方程组Gk*dk=-gk,计算搜索方向 if(norm(gk)epsilion),break;end %检验终止准则 m=0;mk=0; while(m20) %用Armijo搜索步长 if(feval(fun,x0+rhom*dk)feval(fun,x0)+sigma*rhom*gk*dk) mk=m;break; End然后建立fun，gfun和Hess矩阵的m文件如下：function f=fun(x)f=x(1)-x(2)+2*x(1)2+2*x(1)*x(2)+x(2)2;function g=gfun(x) . . . - 21 - / 39g=1

42、+4*x(1)+2*x(2),-1+2*x(1)+2*x(2);function He=Hess(x)n=length(x);He=zeros(n,n);He=4,2;2,2;然后在工作窗口输入：x0=0,0; x,val,k=revisenm(fun,gfun,Hess,x0)所得结果如下x = -1.0000 1.5000val = -1.2500k = 63.3.43.3.4 修正牛顿法的说明修正牛顿法的说明修正 Newton 法克服了 Newton 法的缺点.特别是,当迭代点接近于最优解时,此法具有收敛速度快的优点,对初始点的选择要求不严.但是,修正 Newton 法仍需要计算目标函数

43、的 Hesse 矩阵和逆矩阵,所以求解的计算量和存贮量均很大.另外,当目标函数的 Hesse 矩阵在某点处出现奇异时,迭代将无法进行.因此修正Newton 法仍有相当的局限性3.43.4 共轭共轭梯度法梯度法 首先介绍共轭方向的概念： 定义定义 3.4.13.4.1 设是维欧式空间中的两个向量，即，若有（即, x yn,nx y0Tx y） ，就称与是两个正交的向量，又设是一个阶对称正定矩阵，若有,0 x yxyAn . . . - 22 - / 39，则称向量与关于共轭正交，简称关于共轭0Tx AyxyAA定义定义 3.4.23.4.2 设一组非零向量，为阶对称正定阵，若下 12,nnppp

44、An试成立：.称向量组关于共 0, ,1,2,Tijij i jnpAp 12,npppA轭.也称它们为一组共轭方向（或称为的个共轭方向）.AAn共轭梯度法最早是由 Hestenes 和 Stiefle（1952）提出来的，用于解正定系数矩阵的线性方程组，在这个基础上，Fletcher 和 Reeves （1964）首先提出了解非线性最优化问题的共轭梯度法.由于共轭梯度法不需要矩阵存储，且有较快的收敛速度和二次终止性等优点，现在共轭梯度法已经广泛地应用于实际问题中.共轭梯度法是一个典型的共轭方向法，它的每一个搜索方向是互相共轭的，而这些搜索方向仅仅是负梯度方向与上一次迭代的搜索方向的组合，因d

45、此，存储量少，计算方便3.4.13.4.1 共轭梯度法的基本思想与原理共轭梯度法的基本思想与原理 T1min2TfbCXX AXX其中为阶对称正定阵，是常数An,nnxbC由上是正定二次函数， f X fbXAX故有 . 1kkkkff XXAP0,1,1kn又设以迭代点沿搜索方向进行一维搜索时采用最佳一维搜索求步长因子kXkP，即k10minkkkkkkkkkffXXPXPXP . . . - 23 - / 39进行多次迭代，检验若 ，则就是最优解.否则就继续进行迭代. 0kfXkX3.4.23.4.2 共轭梯度法的迭代步骤共轭梯度法的迭代步骤已知目标函数,终止限 fX0（1） 选取初始点,

46、给定终止限0X0（2） 求初始梯度.计算 ,若停止迭代输出; 否则转0fX0fX0X（3）（3）构造初始搜索方向.取,令,转(4). 00Pf X0k （4）进行一维搜索.求使得,令 kt0minkkkkKtftftXPXP，转（5）1kKkktXXP（5）求梯度向量.计算，若，停止迭代输出；否则1kfX1kfX1kX转 （6）（6）检验迭代次数.若，令，转（3） ；否则转（7）1kn 0=kXX（7）构造共轭方向.取，令， 212kkkffXX11kkkkf PXP1kk转（4）例：例： 用共轭梯度法求，其中 2212min25fxxX6120,2,2,10TTx xXX解解 0001000

47、00244242 1001002100tfttt PXXXP 令，则 22min2425 2 100minttt 500032100160dttdt210110201.919883.8397514.767320.0200370.0014720.003070.1535910016ftffXXXX . . . - 24 - / 39新的搜索方向1003.8397543.845650.0014720.153591000.00615f 1PXP有 2111121211.919883.8456529.5799714.7673000.499230.003070.00615tdtfttttdtXXPXP下一

48、迭代点21111.919883.8456500.499230.003070.006150t XXP由于62010fX 停止迭代输出所求得*20,0TXX3.4.33.4.3 共轭梯度法的共轭梯度法的 matlabmatlab 实现实现首先建立如下 m 文件：function x,val,k=frcg(fun,gfun,x0)%功能：用共轭梯度法求解无约束问题：minf（x）%输入：x0是初始点，fun,gfun分别是目标函数和梯度%输出：x,val分别是近似最优点和最优值，k是迭代次数.maxk=5000; %最大迭代次数rho=0.6;sigma=0.4;k=0; eps=1e-3;n=le

49、ngth(x0);while(k=0.0) d=-g; end end if(norm(g)eps),break;end %检验终止条件 m=0;mk=0; while(m20) %Armijo搜索 if(feval(fun,x0+rhom*d) x0=2,2; x,val,k=frcgfun,gfun,x03.4.43.4.4 共轭梯度法的说明共轭梯度法的说明 实际上,可以把共轭梯度法看作是最速下降法的一种改进.当令时,就变0k为最速下降法.共轭梯度法由于不涉与矩阵,仅仅有向量运算,因而存储量小,适合于维数较高的优化问题.另外,共轭梯度法可以不要求精确的直线搜索.但是,不精确的直线搜索可能导

50、致迭代出来的向量不再共轭,从而降低方法的效能.克服的办法是,重设初始点,即把经过次迭代得到的作为初始点重新迭代.1n1nX计算实践指出,用比用重作初始点要好.1nXnX3.53.5 变尺度法变尺度法变尺度法是在牛顿法的基础上发展起来的，它和梯度法亦有密切关系.变尺度法避免了 Newton 法在每次迭代都要计算目标函数的 Hesse 矩阵和它的逆矩阵而导致随问题的维数增加计算量迅速增加.3.5.13.5.1 变尺度法的基本原理变尺度法的基本原理 . . . - 26 - / 39在 Newton 法中，基本迭代公式1kkkktXXP其中，1kt ，21()()kkkff PXX于是有 （1）11

51、,0,1,2kkkkkXXG g其中是初始点，和分别是目标函数在点的梯度和 Hesse 矩阵.0XkgkG f XkX为了消除这个迭代公式中的 Hesse 逆矩阵，可用某种近似矩阵1kG来替换它，即构造一个矩阵序列去逼近 Hesse 逆矩阵序列kkHH X kH1kG此时上式变为1kkkkXXH g事实上，式中无非是确定了第次迭代的搜索方向，为了取得更大kkk PH gk的灵活性，我们考虑更一般的的迭代公式 （2 2）1kkkkktXXH g其中步长因子通过从出发沿作直线搜索来确定.式（2）是代表ktkXkkk PH g很长的一类迭代公式.例如，当（单位矩阵）时，它变为最速下降法的迭kHI代公

52、式.为使确实与近似并且有容易计算的特点，必须对附加某些条kH1kGkH件：第一， 为保证迭代公式具有下降性质，要求中的每一个矩阵都是 kH对称正定的，理由是，为使搜索方向是下降方向，只要kkk PH g0TTkkkkk g PgH g . . . - 27 - / 39成立即可，即0Tkkkg H g成立.当对称正定时，此公式必然成立，从而保证式（2）具有下降性质.kH第二，要求之间的迭代具有简单形式.显然，kH （3）1kkkHHE是最简单的形式了.其中称为校正矩阵，式（3）称为校正公式.kE第三，必须满足拟 Newton 条件.即： （4）111()()kkkkkHggXX为了书写方便也记

53、1kkkygg1kkkSXX于是拟 Newton 条件可写为 (5)1kkkHyS有式（3）和（5）知，必须满足kE()kkkkHEyS 或 （6）kkkkkE ySH y3.5.23.5.2 DFPDFP 算法的基本思想和原理算法的基本思想和原理DFP 校正是第一个拟牛顿校正是 1959 年由 Davidon 提出的后经 Fletcher和 Powell 改进故名之为 DFP 算法它也是求解无约束优化问题最有效的算法之一.DFP 算法基本原理考虑如下形式的校正公式 (7)1TTkkkkkkkkHHU UV V其中kkVU ,是特定n维向量，kk,是待定常数.这时，校正矩阵是 . . . -

54、28 - / 39TTkkkkkkkEU UV V现在来确定kE根据拟 Newton 条件，必须满足（6） ，于是有kETTkKkkkkkkkkU UV VySH y或.TTkkkkkkkkkkkyU UV V ySH y满足这个方程的待定向量和有无穷多种取法，下面是其中的一种：kUkVTkkkkkU U ySTkkkkkk V V yH y注意到和都是数量，不妨取TkkU yTkkV y，kkUSkkkVH y同时定出，.1kTkkS y1kTkkk y H y将这两式代回（7）得 （8）1TTkkkkkkkkTTkkkkkS SH y y HHHS yy H y这就是 DFP 校正公式.3

55、.5.33.5.3 DFPDFP 法的迭代步骤法的迭代步骤已知目标函数与其梯度，问题的维数，终止限. f X g Xn（1） 选定初始点.计算，.00()ffX00()ggX（2） 置，.0HI00,0k Pg（3） 作直线搜索；计算，.1(,)kkklsXXP11()kkffX11()kkggX . . . - 29 - / 39（4） 判别终止准则是否满足：若满足，则打印，结束；否转11(,)kkfX（5）.（5） 若，则置，转（2） ；否则转kn010101,kkkffXXgg（6）.（6） 计算，1kkkSXX，1kkkygg，1TTkkkkkkkkTTkkkkkS SH y y HH

56、HS yy H y，111kkk PHg 置，转（3）.1kk例例 3.5.13.5.1 用 DFP 法求 22120min4,1,1TxxX解解：取时，DFP 法与最速下降法有一样的的第一代迭代点0HI 由例 2.1 可知 ， 12220,808xgxQX0011120.738461.47692,=180.046160.36923 XgXg，0100100.261541.046160.523088.36923SXXygg0000001000000TTTTS SH y y HHHS yy H y 因为 00000008.89236,70.31762TTTS yy H yy yT00000000

57、0.068400.273610.273614.37778,0.273611.094454.377870.04401TTS SH y y Hy y . . . - 30 - / 39 所以 11.003800.031490.031490.12697H搜索方向为1111.414960.09340 PH g直线搜索2110.738461.9416,0.046160.09340ttXXP由知110dftdtXP0.49423t 所以 200 X20gX是极小点2X3.5.43.5.4 DFPDFP 法的法的 matlabmatlab 实现实现 首先建立 M 文件为：function x,f,tdc=d

58、fp(x0,A,fun,gfun)%功能：用 DFP 算法求解无约束问题：minf(x),x0是初始点.%输出：x,f 分别是近似最优点和最优值，计算运行循环次数.%A 为所求问题 Hesse 矩阵，fun 为问题函数，gfun 为梯度函数.ep=1e-6;%ep 为终止限精度k=1;n=length(x0);Hk=eye(n);X(:,1)=x0;G(:,1)=gfun(X(:,1);X(:,k+1)=X(:,k)-G(:,1)*G(:,1)/(G(:,1)*A*G(:,1)*G(:,1);%一步最速下降法G(:,k+1)=gfun(X(:,k+1);tdc=0;%计算迭代此数while(n

59、orm(G(:,k+1)ep) if k=n+1 x0=X(:,k+1); g0=gfun(x0); Hk=eye(n);k=1; pk=-g0; X(:,k)=x0; G(:,k)=g0; else sk=X(:,k+1)-X(:,k); yk=G(:,k+1)-G(:,k); . . . - 31 - / 39 Hk=Hk+sk*sk/(sk*yk)-Hk*yk*yk*Hk/(yk*Hk*yk); pk=Hk*G(:,k+1); k=k+1; end X(:,k+1)=X(:,k)+(pk*pk/(pk*A*pk)*pk; G(:,k+1)=gfun(X(:,k+1); tdc=tdc+1

60、;endx=X(:,k+1);f=fun(x);对于例3.6.1输入命令： x0=1,1; A=2 0;0 8; fun=inline(x(1)2+4*x(2)2); gfun=inline(2*x(1);8*x(2); x,f,tdc=dfp(x0,A,fun,gfun)即可得到结果3.5.53.5.5 BFGSBFGS 法的基本思想和原理法的基本思想和原理另一个有效和著名的变尺度算法是 Broyden, Fletcher(1970), Goldfarb(1969)和 Shanno(1970)共同研究的结果,因而叫做 BFGS 法. 其基本原理为： 考虑校正公式 ，1TTTTTkkkkkkk