第三章非线性规划 §1非线性规划 1.1非线性规划的实例与定义 如果目标函数或约束条件中包含非线性函数,就称这种规划问题为非线性规划问 题。一般说来,解非线性规划要比解线性规划问题困难得多。而且,也不象线性规划有 单纯形法这一通用方法,非线性规划目前还没有适于各种问题的一般算法,各个方法都 有自己特定的适用范围。 下面通过实例归纳出非线性规划数学模型的一般形式,介绍有关非线性规划的基本 概念 例1(投资决策问题)某企业有n个项目可供选择投资,并且至少要对其中一个 项目投资。已知该企业拥有总资金A元,投资于第i(i=1,…,n)个项目需花资金a1元 并预计可收益b元。试选择最佳投资方案。 解设投资决策变量为 1,决定投资第个项目 ,i=1 0,决定不投资第个项目 则投资总额为∑ax1,投资总收益为∑bx1。因为该公司至少要对一个项目投资,并 且总的投资金额不能超过总资金A,故有限制条件 0<>ax.≤A 另外,由于x,(i=1,…,n)只取值0或1,所以还有 x,(1-x1)=0,i=1, 最佳投资方案应是投资额最小而总收益最大的方案,所以这个最佳投资决策问题归 结为总资金以及决策变量(取0或1)的限制条件下,极大化总收益和总投资之比。因 此,其数学模型为: ∑bx max O=i=l t.0 ax≤A x,(1-x;)=0,i=1,…,n 上面例题是在一组等式或不等式的约束下,求一个函数的最大值(或最小值)问 题,其中目标函数或约束条件中至少有一个非线性函数,这类问题称之为非线性规划问 题,简记为(NP)。可概括为一般形式 min f(x) sLth,(x)≤0,j=1,…q (NP) 8(x)=0,i=1,…,p
-19- 第三章 非线性规划 §1 非线性规划 1.1 非线性规划的实例与定义 如果目标函数或约束条件中包含非线性函数，就称这种规划问题为非线性规划问 题。一般说来，解非线性规划要比解线性规划问题困难得多。而且，也不象线性规划有 单纯形法这一通用方法，非线性规划目前还没有适于各种问题的一般算法，各个方法都 有自己特定的适用范围。 下面通过实例归纳出非线性规划数学模型的一般形式，介绍有关非线性规划的基本 概念。 例 1 （投资决策问题）某企业有 n 个项目可供选择投资，并且至少要对其中一个 项目投资。已知该企业拥有总资金 A 元，投资于第 i(i = 1, ,n) 个项目需花资金 i a 元， 并预计可收益 i b 元。试选择最佳投资方案。 解 设投资决策变量为 = ， 决定不投资第 个项目 决定投资第 个项目 i i xi 0 1, ，i =1, ,n ， 则投资总额为 = n i i i a x 1 ，投资总收益为 = n i i i b x 1 。因为该公司至少要对一个项目投资，并 且总的投资金额不能超过总资金 A ，故有限制条件 = n i ai xi A 1 0 另外，由于 x (i 1, ,n) i = 只取值 0 或 1，所以还有 x (1 x ) 0, i 1, ,n. i − i = = 最佳投资方案应是投资额最小而总收益最大的方案，所以这个最佳投资决策问题归 结为总资金以及决策变量（取 0 或 1）的限制条件下，极大化总收益和总投资之比。因 此，其数学模型为： = = = n i i i n i i i a x b x Q 1 1 max s.t. = n i ai xi A 1 0 x (1 x ) 0, i 1, ,n. i − i = = 上面例题是在一组等式或不等式的约束下，求一个函数的最大值（或最小值）问 题，其中目标函数或约束条件中至少有一个非线性函数，这类问题称之为非线性规划问 题，简记为（NP）。可概括为一般形式 min f (x) s.t. hj (x) 0, j = 1, , q (NP) gi (x) = 0, i = 1, , p
其中x=[x1…xn称为模型(N)的决策变量,f称为目标函数,g(i=1…,p) 和h,(j=1,…,q)称为约束函数。另外,g,(x)=0(=1…,p)称为等式约束, h(x)≤0(=1…,q)称为不等式约束。 对于一个实际问题,在把它归结成非线性规划问题时,一般要注意如下几点: (i)确定供选方案:首先要收集同问题有关的资料和数据,在全面熟悉问题的基 础上,确认什么是问题的可供选择的方案,并用一组变量来表示它们。 (ⅱ)提出追求目标:经过资料分析,根据实际需要和可能,提出要追求极小化 或极大化的目标。并且,运用各种科学和技术原理,把它表示成数学关系式。 (i)给出价值标准:在提出要追求的目标之后,要确立所考虑目标的“好”或 坏”的价值标准,并用某种数量形式来描述它 (iv)寻求限制条件:由于所追求的目标一般都要在一定的条件下取得极小化或 极大化效果,因此还需要寻找出问题的所有限制条件,这些条件通常用变量之间的一些 不等式或等式来表示 1.2线性规划与非线性规划的区别 如果线性规划的最优解存在,其最优解只能在其可行域的边界上达到(特别是可行 域的顶点上达到);而非线性规划的最优解(如果最优解存在)则可能在其可行域的任 意一点达到 1.3非线性规划的 Matlab解法 Matlab中非线性规划的数学模型写成以下形式 min f(x) Ax< B Aeg·x=Beq C(x)≤0 Ceq(x)=0 其中f(x)是标量函数,A,B,Aeg,Beq是相应维数的矩阵和向量,C(x).Ceq(x)是非 线性向量函数。 Matlab中的命令是 X=FMINCON(FUN, XO, A, B, Aeg, Beg, LB, UB, NONLCON, OPTIONS) 它的返回值是向量x,其中FUN是用M文件定义的函数f(x):X0是x的初始值 A,B,Aeq,Beq定义了线性约束A*X≤B,Aeq*X=Beq,如果没有等式约束,则 A=[],B=[],Aeq=[],Beq=[]:LB和UB是变量x的下界和上界,如果上界和下界没有约 束,则LB=[],UB=[],如果x无下界,则LB=-inf,如果x无上界,则UB=inf; NONLCON 是用M文件定义的非线性向量函数C(x),Ceq(x); OPTIONS定义了优化参数,可以使 用 Matlab缺省的参数设置。 例2求下列非线性规划问题 min f(x)=x,+x,+8 x-x2≥0 ≥0
-20- 其中 T n x [x x ] = 1 称为模型（NP）的决策变量， f 称为目标函数， i g (i = 1, , p) 和 h ( j 1, ,q) j = 称为约束函数。另外， gi (x) = 0 (i = 1, , p) 称为等式约束， hj (x) 0 ( j = 1, ,q) 称为不等式约束。 对于一个实际问题，在把它归结成非线性规划问题时，一般要注意如下几点： （i）确定供选方案：首先要收集同问题有关的资料和数据，在全面熟悉问题的基 础上，确认什么是问题的可供选择的方案，并用一组变量来表示它们。 （ii）提出追求目标：经过资料分析，根据实际需要和可能，提出要追求极小化 或极大化的目标。并且，运用各种科学和技术原理，把它表示成数学关系式。 （iii）给出价值标准：在提出要追求的目标之后，要确立所考虑目标的“好”或 “坏”的价值标准，并用某种数量形式来描述它。 （iv）寻求限制条件：由于所追求的目标一般都要在一定的条件下取得极小化或 极大化效果，因此还需要寻找出问题的所有限制条件，这些条件通常用变量之间的一些 不等式或等式来表示。 1.2 线性规划与非线性规划的区别 如果线性规划的最优解存在，其最优解只能在其可行域的边界上达到（特别是可行 域的顶点上达到）；而非线性规划的最优解（如果最优解存在）则可能在其可行域的任 意一点达到。 1.3 非线性规划的 Matlab 解法 Matlab 中非线性规划的数学模型写成以下形式 min f (x) = = ( ) 0 ( ) 0 Ceq x C x Aeq x Beq Ax B , 其中 f (x) 是标量函数， A, B, Aeq, Beq 是相应维数的矩阵和向量， C(x),Ceq(x) 是非 线性向量函数。 Matlab 中的命令是 X=FMINCON(FUN,X0,A,B,Aeq,Beq,LB,UB,NONLCON,OPTIONS) 它的返回值是向量 x ，其中 FUN 是用 M 文件定义的函数 f (x) ；X0 是 x 的初始值； A,B,Aeq,Beq 定义了线性约束 A* X B, Aeq * X = Beq ，如果没有等式约束，则 A=[],B=[],Aeq=[],Beq=[]；LB 和 UB 是变量 x 的下界和上界，如果上界和下界没有约 束，则 LB=[]，UB=[]，如果 x 无下界，则 LB=-inf，如果 x 无上界，则 UB=inf；NONLCON 是用 M 文件定义的非线性向量函数 C(x),Ceq(x) ；OPTIONS 定义了优化参数，可以使 用 Matlab 缺省的参数设置。 例 2 求下列非线性规划问题 − − + = − = + + , 0. 2 0 0 min ( ) 8 1 2 2 1 2 2 2 1 2 2 2 1 x x x x x x f x x x
(i)编写M文件fun1.m function f=funl(x) f=x(1)^2+x(2)^2+8 和M文件fun2.m function [g, h]=fun2(x)i g=-x(1)^2+x(2); h=-x(1)-x(2)^2+2;号等式约束 (i)在 Matlab的命令窗口依次输入 [x,y]= fmincon("fun1',rand(2,1),[],[],[],[], zeros(2,1),[], options 就可以求得当x1=L,x2=1时,最小值y=10 14求解非线性规划的基本迭代格式 记(NP)的可行域为K 若x∈K,并且 f(x)≤f(x),Vx∈K 则称x是(NP)的整体最优解,f(x)是(NP)的整体最优值。如果有 f(x)0,使
-21- （i）编写 M 文件 fun1.m function f=fun1(x); f=x(1)^2+x(2)^2+8; 和M文件fun2.m function [g,h]=fun2(x); g=-x(1)^2+x(2); h=-x(1)-x(2)^2+2; %等式约束 （ii）在 Matlab 的命令窗口依次输入 options=optimset; [x,y]=fmincon('fun1',rand(2,1),[],[],[],[],zeros(2,1),[], ... 'fun2', options) 就可以求得当 x1 =1, x2 =1 时，最小值 y = 10 。 1.4 求解非线性规划的基本迭代格式 记（NP）的可行域为 K 。 若 x K * ，并且 f (x ) f (x), xK * 则称 * x 是（NP）的整体最优解， ( ) * f x 是(NP)的整体最优值。如果有 * * f (x ) f (x), xK, x x 则称 * x 是（NP）的严格整体最优解， ( ) * f x 是(NP)的严格整体最优值。 若 x K * ，并且存在 * x 的邻域 ( ) * N x ，使 f (x ) f (x), x N (x ) K * * ， 则称 * x 是（NP）的局部最优解， ( ) * f x 是(NP)的局部最优值。如果有 f (x ) f (x), x N (x ) K * * 则称 * x 是（NP）的严格局部最优解， ( ) * f x 是(NP)的严格局部最优值。 由于线性规划的目标函数为线性函数，可行域为凸集，因而求出的最优解就是整 个可行域上的全局最优解。非线性规划却不然，有时求出的某个解虽是一部分可行域上 的极值点，但并不一定是整个可行域上的全局最优解。 对于非线性规划模型(NP)，可以采用迭代方法求它的最优解。迭代方法的基本思 想是：从一个选定的初始点 n x R 出发，按照某一特定的迭代规则产生一个点列 { } k x ，使得当 { } k x 是有穷点列时，其最后一个点是(NP)的最优解；当 { } k x 是无穷点 列时，它有极限点，并且其极限点是(NP)的最优解。 设 k n x R 是某迭代方法的第 k 轮迭代点， k n x R +1 是第 k +1 轮迭代点，记 k k k k x = x + t p +1 （1） 这里 , , 1 1 = k n k t k R p R p ，显然 k p 是由点 k x 与点 k +1 x 确定的方向。式（1）就 是求解非线性规划模型(NP)的基本迭代格式。 通常，我们把基本迭代格式（1）中的 k p 称为第 k 轮搜索方向， k t 为沿 k p 方向的 步长，使用迭代方法求解(NP)的关键在于，如何构造每一轮的搜索方向和确定适当的步 长。 设 x R , p 0 n ，若存在 0 ，使
f(x+l)0,使 x+tp∈K 称向量p是点x处关于K的可行方向。 个向量p,若既是函数∫在点x处的下降方向,又是该点关于区域K的可行方 向,则称之为函数∫在点x处关于K的可行下降方向 现在,我们给出用基本迭代格式(1)求解(NP)的一般步骤如下: 0°选取初始点x,令k:=0。 1°构造搜索方向,依照一定规则,构造∫在点x处关于K的可行下降方向作为 搜索方向p 2°寻求搜索步长。以x为起点沿搜索方向p寻求适当的步长t,使目标函数值 有某种意义的下降 3°求出下一个迭代点。按迭代格式(1)求出 若x*已满足某种终止条件,停止迭代。 4°以x+代替x,回到1°步。 1.5凸函数、凸规划 设f(x)为定义在n维欧氏空间E中某个凸集R上的函数,若对任何实数 a(0<a<1)以及R中的任意两点x和x2),恒有 f(a(+(1-a)x2)≤f(x)+(1-a)(x2) 则称∫(x)为定义在R上的凸函数 若对每一个a(0<a<1)和x≠x(2)∈R恒有 f(ax"+(1-a)x2)<f(x)+(1-a)f(x2) 则称∫(x)为定义在R上的严格凸函数。 考虑非线性规划 fo R={x|g,(x)≤0,j=1,2,…,l} 假定其中f(x)为凸函数,g1(x)j=1,2,…D为凸函数,这样的非线性规划称为 凸规划 可以证明,凸规划的可行域为凸集,其局部最优解即为全局最优解,而且其最优 解的集合形成一个凸集。当凸规划的目标函数∫(x)为严格凸函数时,其最优解必定唯 (假定最优解存在)。由此可见,凸规划是一类比较简单而又具有重要理论意义的非 线性规划 §2无约束问题 21一维搜索方法 当用迭代法求函数的极小点时,常常用到一维搜索,即沿某一已知方向求目标函数
-22- f (x + tp) f (x), t (0, )， 称向量 p 是 f 在点 x 处的下降方向。 设 x R , p 0 n ，若存在 t 0 ，使 x +tpK ， 称向量 p 是点 x 处关于 K 的可行方向。 一个向量 p ，若既是函数 f 在点 x 处的下降方向，又是该点关于区域 K 的可行方 向，则称之为函数 f 在点 x 处关于 K 的可行下降方向。 现在，我们给出用基本迭代格式（1）求解(NP)的一般步骤如下： 0°选取初始点 0 x ，令 k := 0。 1° 构造搜索方向，依照一定规则，构造 f 在点 k x 处关于 K 的可行下降方向作为 搜索方向 k p 。 2° 寻求搜索步长。以 k x 为起点沿搜索方向 k p 寻求适当的步长 k t ，使目标函数值 有某种意义的下降。 3°求出下一个迭代点。按迭代格式（1）求出 k k k k x = x + t p +1 。 若 k +1 x 已满足某种终止条件，停止迭代。 4°以 k +1 x 代替 k x ，回到 1°步。 1.5 凸函数、凸规划 设 f (x) 为定义在 n 维欧氏空间 (n) E 中某个凸集 R 上的函数，若对任何实数 (0 1) 以及 R 中的任意两点 (1) x 和 (2) x ，恒有 ( (1 ) ) ( ) (1 ) ( ) (1) (2) (1) (2) f x + − x f x + − f x 则称 f (x) 为定义在 R 上的凸函数。 若对每一个 (0 1) 和 x x R (1) (2) 恒有 ( (1 ) ) ( ) (1 ) ( ) (1) (2) (1) (2) f x + − x f x + − f x 则称 f (x) 为定义在 R 上的严格凸函数。 考虑非线性规划 = = { | ( ) 0, 1,2, , } min ( ) R x g x j l f x j x R 假定其中 f (x) 为凸函数， g (x)( j 1,2, ,l) j = 为凸函数，这样的非线性规划称为 凸规划。 可以证明，凸规划的可行域为凸集，其局部最优解即为全局最优解，而且其最优 解的集合形成一个凸集。当凸规划的目标函数 f (x) 为严格凸函数时，其最优解必定唯 一（假定最优解存在）。由此可见，凸规划是一类比较简单而又具有重要理论意义的非 线性规划。 §2 无约束问题 2.1 一维搜索方法 当用迭代法求函数的极小点时，常常用到一维搜索，即沿某一已知方向求目标函数
的极小点。一维搜索的方法很多,常用的有:(1)试探法(“成功一失败”,斐波那契法 0618法等):插值法(抛物线插值法,三次插值法等):(3)微积分中的求根法(切线 法,二分法等)。 考虑一维极小化问题 min f(o 若f(m)是[a,b区间上的下单峰函数,我们介绍通过不断地缩短[a,6的长度,来 搜索得(2)的近似最优解的两个方法。 为了缩短区间[ab],逐步搜索得(2)的最优解r的近似值,我们可以采用以下 途径:在[a,b]中任取两个关于[a,b]是对称的点和t2(不妨设l20,这就要求最后区间的长度不超过δ,即 (4) 据此,我们应按照预先给定的精度dδ,确定使(4)成立的最小整数n作为搜索次数
-23- 的极小点。一维搜索的方法很多，常用的有：（1）试探法（“成功—失败”，斐波那契法， 0.618 法等）；插值法（抛物线插值法，三次插值法等）；（3）微积分中的求根法（切线 法，二分法等）。 考虑一维极小化问题 min f (t) atb （2） 若 f (t) 是 [a,b] 区间上的下单峰函数，我们介绍通过不断地缩短 [a,b] 的长度，来 搜索得（2）的近似最优解的两个方法。 为了缩短区间 [a,b] ，逐步搜索得（2）的最优解 * t 的近似值，我们可以采用以下 途径：在 [a,b] 中任取两个关于 [a,b] 是对称的点 1 t 和 2 t （不妨设 2 1 t t ，并把它们叫 做搜索点），计算 ( ) 1 f t 和 ( ) 2 f t 并比较它们的大小。对于单峰函数，若 ( ) ( ) 2 1 f t f t ， 则必有 [ , ] 1 * t a t ，因而 [ , ] 1 a t 是缩短了的 单峰区 间；若 ( ) ( ) 1 2 f t f t ，则有 [ , ] 2 * t t b ，故 [ , ] t 2 b 是缩短了的单峰区间；若 ( ) ( ) 2 1 f t = f t ，则 [ , ] 1 a t 和 [ , ] t 2 b 都是 缩短了的单峰。因此通过两个搜索点处目标函数值大小的比较，总可以获得缩短了的单 峰区间。对于新的单峰区间重复上述做法，显然又可获得更短的单峰区间。如此进行， 在单峰区间缩短到充分小时，我们可以取最后的搜索点作为（2）最优解的近似值。 应该按照怎样的规则来选取探索点，使给定的单峰区间的长度能尽快地缩短？ 2.1.1 Fibonacci 法 若数列{ Fn }满足关系： F0 = F1 = 1 , 2,3, , Fn = Fn−2 + Fn−1 n = 则称 { } Fn 为 Fibonacci 数列， Fn 称为第 n 个 Fibonacci 数，称相邻两个 Fibonacci 数之 比 n n F F −1 为 Fibonacci 分数。 当用斐波那契法以 n 个探索点来缩短某一区间时，区间长度的第一次缩短率为 n n F F −1 ，其后各次分别为 2 1 2 1 3 , , , F F F F F F n n n n − − − 。由此，若 1 t 和 ( ) 2 2 1 t t t 是单峰区间 [a,b] 中第 1 个和第 2 个探索点的话，那么应有比例关系 n n F F b a t 1 a −1 = − − ， n n F F b a t 2 a −2 = − − 从而 ( ) 1 1 b a F F t a n n = + − − ， ( ) 2 2 b a F F t a n n = + − − （3） 它们关于 [a,b] 确是对称的点。 如果要求经过一系列探索点搜索之后，使最后的探索点和最优解之间的距离不超 过精度 0 ，这就要求最后区间的长度不超过 ，即 − Fn b a （4） 据此，我们应按照预先给定的精度 ，确定使（4）成立的最小整数 n 作为搜索次数
直到进行到第n个探索点时停止。 用上述不断缩短函数f(1)的单峰区间[a,b]的办法,来求得问题(2)的近似解, 是 Kiefer((1953年)提出的,叫做 Fibonacci法,具体步骤如下 选取初始数据,确定单峰区间[ao3b],给出搜索精度δ>0,由(4)确定搜 索次数n 2°k=1,a=ao,b=b,计算最初两个搜索点,按(3)计算和t2 3° while k0,满足比例关系 称之为黄金分割数,其值为o=√5-1 0.6180339887 黄金分割数O和 Fibonacci分数之间有着重要的关系,它们是
-24- 直到进行到第 n 个探索点时停止。 用上述不断缩短函数 f (t) 的单峰区间 [a,b] 的办法，来求得问题（2）的近似解， 是 Kiefer(1953 年)提出的，叫做 Finbonacci 法，具体步骤如下： 1° 选取初始数据，确定单峰区间 [ , ] a0 b0 ，给出搜索精度 0 ，由（4）确定搜 索次数 n。 2° 0 0 k =1,a = a ,b = b ，计算最初两个搜索点，按（3）计算 1 t 和 2 t 。 3°while k n −1 ( ), ( ) 1 1 2 2 f = f t f = f t if 1 2 f f ( ) ( ) ( 1 ) ; ; 2 2 1 1 b a F n k F n k a t t t t a − − − − = = = + else ( ) ( ) ( 1 ) ; ; 1 1 2 2 a b F n k F n k b t t t t b − − − − = = = + end k = k +1 end 4° 当进行至 k = n −1 时， ( ) 2 1 t 1 = t 2 = a + b 这就无法借比较函数值 ( ) 1 f t 和 ( ) 2 f t 的大小确定最终区间，为此，取 = + + − = + )( ) 2 1 ( ( ) 2 1 1 2 t a b a t a b 其中 为任意小的数。在 1 t 和 2 t 这两点中，以函数值较小者为近似极小点，相应的函数 值为近似极小值。并得最终区间 [ , ] 1 a t 或 [ , ] t 2 b 。 由上述分析可知，斐波那契法使用对称搜索的方法，逐步缩短所考察的区间，它能 以尽量少的函数求值次数，达到预定的某一缩短率。 例 3 试用斐波那契法求函数 ( ) 2 2 f t = t −t + 的近似极小点，要求缩短后的区间 不大于区间 [−1,3] 的 0.08 倍。 程序留作习题。 2.1.2 0.618 法 若 0 ，满足比例关系 − = 1 1 称之为黄金分割数，其值为 0.6180339887 2 5 1 = − = 。 黄金分割数 和 Fibonacci 分数之间有着重要的关系，它们是
F 2n<0< n为偶数 FO、、F n为奇数。 2°c=lim 现用不变的区间缩短率0.618,代替斐波那契法每次不同的缩短率,就得到了黄金 分割法(0.618法)。这个方法可以看成是斐波那契法的近似,实现起来比较容易,效果 也相当好,因而易于为人们所接受。 用0618法求解,从第2个探索点开始每增加一个探索点作一轮迭代以后,原单 峰区间要缩短0618倍。计算n个探索点的函数值可以把原区间[a,b]连续缩短n-1 次,因为每次的缩短率均为,故最后的区间长度为 b-a0) 这就是说,当已知缩短的相对精度为δ时,可用下式计算探索点个数n <6 当然,也可以不预先计算探索点的数目n,而在计算过程中逐次加以判断,看是否已满 足了提出的精度要求 0.618法是一种等速对称进行试探的方法,每次的探索点均取在区间长度的0.618 倍和0.382倍处 22二次插值法 对极小化问题(2),当∫()在[a,b]上连续时,可以考虑用多项式插值来进行一 维搜索。它的基本思想是:在搜索区间中,不断用低次(通常不超过三次)多项式来近 似目标函数,并逐步用插值多项式的极小点来逼近(2)的最优解。 23无约束极值问题的解法 无约束极值问题可表述为 mnf(x),x∈E(n (5) 求解问题(5)的迭代法大体上分为两种:一是用到函数的一阶导数或二阶导数, 称为解析法。另一是仅用到函数值,称为直接法。 2.3.1解析法 2.3.1.1梯度法(最速下降法) 对基本迭代格式 x +t,p 我们总是考虑从点x出发沿哪一个方向p,使目标函数∫下降得最快。微积分的知识 告诉我们,点x的负梯度方向 P 是从点x出发使∫下降最快的方向。为此,称负梯度方向-V(x2)为∫在点x处的 最速下降方向。 按基本迭代格式(6),每一轮从点x*出发沿最速下降方向-Vf(x)作一维搜索 来建立求解无约束极值问题的方法,称之为最速下降法
-25- 1° 1 1 − − n n n n F F F F ，n 为偶数， 1 1 − − n n n n F F F F ，n 为奇数。 2° n n n F F 1 lim − → = 。 现用不变的区间缩短率 0.618，代替斐波那契法每次不同的缩短率，就得到了黄金 分割法（0.618 法）。这个方法可以看成是斐波那契法的近似，实现起来比较容易，效果 也相当好，因而易于为人们所接受。 用 0.618 法求解，从第 2 个探索点开始每增加一个探索点作一轮迭代以后，原单 峰区间要缩短 0.618 倍。计算 n 个探索点的函数值可以把原区间 [ , ] a0 b0 连续缩短 n −1 次，因为每次的缩短率均为 ，故最后的区间长度为 1 0 0 ( ) − − n b a 这就是说，当已知缩短的相对精度为 时，可用下式计算探索点个数 n ： n−1 当然，也可以不预先计算探索点的数目 n ，而在计算过程中逐次加以判断，看是否已满 足了提出的精度要求。 0.618 法是一种等速对称进行试探的方法，每次的探索点均取在区间长度的 0.618 倍和 0.382 倍处。 2.2 二次插值法 对极小化问题（2），当 f (t) 在 [a,b] 上连续时，可以考虑用多项式插值来进行一 维搜索。它的基本思想是：在搜索区间中，不断用低次（通常不超过三次）多项式来近 似目标函数，并逐步用插值多项式的极小点来逼近（2）的最优解。 2.3 无约束极值问题的解法 无约束极值问题可表述为 ( ) min ( ), n f x x E （5） 求解问题（5）的迭代法大体上分为两种：一是用到函数的一阶导数或二阶导数， 称为解析法。另一是仅用到函数值，称为直接法。 2.3.1 解析法 2.3.1.1 梯度法（最速下降法） 对基本迭代格式 k k k k x = x + t p +1 （6） 我们总是考虑从点 k x 出发沿哪一个方向 k p ，使目标函数 f 下降得最快。微积分的知识 告诉我们，点 k x 的负梯度方向 ( ) k k p = −f x ， 是从点 k x 出发使 f 下降最快的方向。为此，称负梯度方向 ( ) k −f x 为 f 在点 k x 处的 最速下降方向。 按基本迭代格式（6），每一轮从点 k x 出发沿最速下降方向 ( ) k −f x 作一维搜索， 来建立求解无约束极值问题的方法，称之为最速下降法
这个方法的特点是,每轮的搜索方向都是目标函数在当前点下降最快的方向。同 时,用Vf(x)=0或Ⅳ(x2)≤E作为停止条件。其具体步骤如下 1°选取初始数据。选取初始点x0,给定终止误差,令k:=0。 2求梯度向量。计算V(x2),若x)≤E,停止达代,输出x2。否则, 3°构造负梯度方向。取 4°进行一维搜索。求l,使得 f(x+tip)=mn f(x"+tp) 令 k:=k+1,转 例4用最速下降法求解无约束非线性规划问题 min f(x)=xi+25x2 其中x=(x1x2),要求选取初始点x"=(2,2),终止误差E=10。 解:(i)Vf(x)=(2x1,50x2) 编写M文件 deaf. m如下 function [f, dfl=deaf(x)i f=x(1)^2+25+x(2)^2 df(1)=2*x(1) df(2)=50+x(2); (i)编写M文件 Zulu.m x=[2;2] f0, g]=deaf(x)i while norm(g)>0.000001 p=-g/ norm(g) t=l 0; f=deaf (x+t*p)i while f>fo t=t/2; f=deaf (x+t*p)i =X+t大 f0, g]=deaf(x) 23.1.2 Newton法 考虑目标函数∫在点x处的二次逼近式 f(x)Q(x)=f(x)+Vf(x)(x-x)+-(x-x)'Vf(x)(x-x) 假定 Hesse阵
-26- 这个方法的特点是，每轮的搜索方向都是目标函数在当前点下降最快的方向。同 时，用 ( ) = 0 k f x 或 ( ) k f x 作为停止条件。其具体步骤如下： 1°选取初始数据。选取初始点 0 x ，给定终止误差，令 k := 0。 2°求梯度向量。计算 ( ) k f x ， 若 ( ) k f x ，停止迭代，输出 k x 。否则， 进行 3°。 3°构造负梯度方向。取 ( ) k k p = −f x . 4° 进行一维搜索。求 k t ，使得 ( ) min ( ) 0 k k t k k k f x + t p = f x + tp 令 , : 1, 1 = + = + + x x t p k k k k k k 转 2°。 例 4 用最速下降法求解无约束非线性规划问题 2 2 2 min f (x) = x1 + 25x 其中 T x (x , x ) = 1 2 ，要求选取初始点 T x (2,2) 0 = ，终止误差 6 10− = 。 解：（i） T f (x) (2x ,50x ) = 1 2 编写 M 文件 detaf.m 如下 function [f,df]=detaf(x); f=x(1)^2+25*x(2)^2; df(1)=2*x(1); df(2)=50*x(2); （ii）编写M文件zuisu.m clc x=[2;2]; [f0,g]=detaf(x); while norm(g)>0.000001 p=-g'/norm(g); t=1.0;f=detaf(x+t*p); while f>f0 t=t/2;f=detaf(x+t*p); end x=x+t*p [f0,g]=detaf(x) end 2.3.1.2 Newton 法 考虑目标函数 f 在点 k x 处的二次逼近式 ( ) ( )( ) 2 1 ( ) ( ) ( ) ( ) ( ) k k T k k T 2 k k f x Q x = f x + f x x − x + x − x f x x − x 假定 Hesse 阵
af(x) af(x) ax ax f(x) a(x) 正定 由于V2f(x4)正定,函数Q的稳定点x4+是Q(x)的最小点。为求此最小点,令 VQ(x*+)=Vf(x)+V2f(x*)x-x)=0, 即可解得 x=x-lV f(x Vf(x) 对照基本迭代格式(1),可知从点x出发沿搜索方向 p=-lVf(x] Vf(x) 并取步长t4=1即可得Q(x)的最小点x。通常,把方向p叫做从点x2出发的 Newton方向。从一初始点开始,每一轮从当前迭代点出发,沿 Newton方向并取步长 为1的求解方法,称之为 Newton法。其具体步骤如下: 1°选取初始数据。选取初始点x°,给定终止误差E>0,令k=0 2°求梯度向量。计算V(x3),若Vf(x)≤E,停止迭代,输出x。否则, 进行3° 3°构造 Newton方向。计算[v2f(x2),取 p=-[Vf(x] Vf(x) 4°求下一迭代点。令x+=x+p4,k=k+1,转2° 例5用 Newton法求解 min f(x)=x1+25x 2+xixi 选取x°=(2,2),E=10°。 解:(i)V(x)=[4x1+2x1x2100x2+2x2x2 Vf 2 4x1x2300x2+4x2 编写M文件 nwfun m如下: function [f, df, d2f]=nwfun (x)i df(1)=4+x(1)^3+2*x(1)*x(2)^2 df(2)=100*x(2)^3+2+x(1)^2*x(2); d2f(1,1)=12*x(1)^2+2+x(2)^2 d2f(1,2)=4+x(1)*x(2) (2,1)=d2f(1,2) d2f(2,2)=300*x(2)^2+4+x(1)+x(2) (i)编写M文件: x=[2;2]
-27- = 2 2 1 2 1 2 2 1 2 2 ( ) ( ) ( ) ( ) ( ) n k n k n k k k x f x x x f x x x f x x f x f x 正定。 由于 ( ) 2 k f x 正定，函数 Q 的稳定点 k +1 x 是 Q(x) 的最小点。为求此最小点，令 ( ) ( ) ( )( ) 0 1 2 1 = + − = k + k k k+ k Q x f x f x x x ， 即可解得 [ ( )] ( ) k 1 k 2 k 1 k x = x − f x f x + − . 对照基本迭代格式（1），可知从点 k x 出发沿搜索方向。 [ ( )] ( ) k 2 k 1 k p = − f x f x − 并取步长 t k = 1 即可得 Q(x) 的最小点 k +1 x 。通常，把方向 k p 叫做从点 k x 出发的 Newton 方向。从一初始点开始，每一轮从当前迭代点出发，沿 Newton 方向并取步长 为 1 的求解方法，称之为 Newton 法。其具体步骤如下： 1°选取初始数据。选取初始点 0 x ，给定终止误差 0 ，令 k := 0。 2°求梯度向量。计算 ( ) k f x ，若 ( ) k f x ，停止迭代，输出 k x 。否则， 进行 3°。 3°构造 Newton 方向。计算 2 1 [ ( )]− k f x ，取 [ ( )] ( ) k 2 k 1 k p = − f x f x − . 4° 求下一迭代点。令 , : 1 1 = + = + + x x p k k k k k ，转 2°。 例 5 用 Newton 法求解， 2 2 2 1 4 2 4 min f (x) = x1 + 25x + x x 选取 T x (2,2) 0 = ， 6 10− = 。 解：（i） T f (x) [4x 2x x 100x 2x x ] 2 2 1 3 2 2 1 2 3 = 1 + + + + = 2 1 2 1 2 2 1 2 2 2 2 2 1 4 300 4 12 2 4 x x x x x x x x f 编写 M 文件 nwfun.m 如下： function [f,df,d2f]=nwfun(x); f=x(1)^4+25*x(2)^4+x(1)^2*x(2)^2; df(1)=4*x(1)^3+2*x(1)*x(2)^2; df(2)=100*x(2)^3+2*x(1)^2*x(2); d2f(1,1)=12*x(1)^2+2*x(2)^2; d2f(1,2)=4*x(1)*x(2); d2f(2,1)=d2f(1,2); d2f(2,2)=300*x(2)^2+4*x(1)*x(2); （ii）编写 M 文件： clc x=[2;2];
f0, gl, g2]=nw fun(x) while norm(g1)>0.00001 号dead1oop,fori=1:3 p=-inv(g2)*gl, p=p/norm(p) t=l 0, f=detaf(x+t*p while f>fO t=t/2, f=detaf (x+t*p) end [fo, gl, g2]=nw fun(x) 如果目标函数是非二次函数,一般地说,用 Newton法通过有限轮迭代并不能保证 可求得其最优解。 Newton法的优点是收敛速度快;缺点是有时不妤用而需采取改进措施,此外,当 维数较高时,计算-[V2f(x2)的工作量很大。 2.3.1.3变尺度法 变尺度法( Variable Metric Algorithm)是近20多年来发展起来的,它不仅是求解 无约束极值问题非常有效的算法,而且也已被推广用来求解约束极值问题。由于它既避 免了计算二阶导数矩阵及其求逆过程,又比梯度法的收敛速度快,特别是对高维问题具 有显著的优越性,因而使变尺度法获得了很高的声誉。下面我们就来简要地介绍一种变 尺度法一DFP法的基本原理及其计算过程。这一方法首先由 Davidon在1959年提出, 后经 Fletcher和 Powell加以改进。 我们已经知道,牛顿法的搜索方向是-V2f(x)vf(x2),为了不计算二阶 导数矩阵[Vf(x)及其逆阵,我们设法构造另一个矩阵,用它来逼近二阶导数矩阵 的逆阵V2f(x),这一类方法也称拟牛顿法( Quasi-Newton Method)。 下面研究如何构造这样的近似矩阵,并将它记为H)。我们要求:每一步都能以 现有的信息来确定下一个搜索方向;每做一次选代,目标函数值均有所下降;这些近似 矩阵最后应收敛于解点处的Hese阵的逆阵。 当f(x)是二次函数时,其Hese阵为常数阵A,任两点x和x处的梯度之差为 VfO Vf(x=A( 或 x+l-x=A-IVf(x+)-Vf(x') 对于非二次函数,仿照二次函数的情形,要求其 Hesse阵的逆阵的第k+1次近似 矩阵H(+满足关系式 x*+l-x=H+IVf(x+)-Vf(x') (7) 这就是常说的拟 Newton条件。 若令 △G(=Vf(x4)-Vf(x2) (8) 则式(7)变为 Axk=H(+△G 现假定厅已知,用下式求H(+1(设F()和H均为对称正定阵
-28- [f0,g1,g2]=nwfun(x) while norm(g1)>0.00001 %dead loop,for i=1:3 p=-inv(g2)*g1',p=p/norm(p) t=1.0,f=detaf(x+t*p) while f>f0 t=t/2,f=detaf(x+t*p), end x=x+t*p [f0,g1,g2]=nwfun(x) end 如果目标函数是非二次函数，一般地说，用 Newton 法通过有限轮迭代并不能保证 可求得其最优解。 Newton 法的优点是收敛速度快；缺点是有时不好用而需采取改进措施，此外，当 维数较高时，计算 2 1 [ ( )]− − k f x 的工作量很大。 2.3.1.3 变尺度法 变尺度法（Variable Metric Algorithm）是近 20 多年来发展起来的，它不仅是求解 无约束极值问题非常有效的算法，而且也已被推广用来求解约束极值问题。由于它既避 免了计算二阶导数矩阵及其求逆过程，又比梯度法的收敛速度快，特别是对高维问题具 有显著的优越性，因而使变尺度法获得了很高的声誉。下面我们就来简要地介绍一种变 尺度法—DFP 法的基本原理及其计算过程。这一方法首先由 Davidon 在 1959 年提出， 后经 Fletcher 和 Powell 加以改进。 我们已经知道，牛顿法的搜索方向是 [ ( )] ( ) 2 k 1 k − f x f x − ，为了不计算二阶 导数矩阵 [ ( )] 2 k f x 及其逆阵，我们设法构造另一个矩阵，用它来逼近二阶导数矩阵 的逆阵 2 1 [ ( )]− k f x ，这一类方法也称拟牛顿法（Quasi-Newton Method）。 下面研究如何构造这样的近似矩阵，并将它记为 (k ) H 。我们要求：每一步都能以 现有的信息来确定下一个搜索方向；每做一次选代，目标函数值均有所下降；这些近似 矩阵最后应收敛于解点处的 Hesse 阵的逆阵。 当 f (x) 是二次函数时，其 Hesse 阵为常数阵 A ，任两点 k x 和 k +1 x 处的梯度之差为 ( ) ( ) ( ) k 1 k k 1 k f x −f x = A x − x + + 或 [ ( ) ( )] k 1 k 1 k 1 k x − x = A f x −f x + − + 对于非二次函数，仿照二次函数的情形，要求其 Hesse 阵的逆阵的第 k +1 次近似 矩阵 (k +1) H 满足关系式 [ ( ) ( )] k 1 k (k 1) k 1 k x − x = H f x −f x + + + （7） 这就是常说的拟 Newton 条件。 若令 = − = − + + k k k k k k x x x G f x f x 1 ( ) 1 ( ) ( ) （8） 则式（7）变为 k (k 1) (k ) x = H G + ， （9） 现假定 (k ) H 已知，用下式求 (k +1) H （设 (k ) H 和 (k +1) H 均为对称正定阵）；