第三章非线性规划.docx

1、第三章非线性规划第三章非线性规划 1非线性规划1.1非线性规划的实例与定义如果目标函数或约束条件中包含非线性函数，就称这种规划问题为非线性规划问 题。一般说来，解非线性规划要比解线性规划问题困难得多。而且，也不象线性规划有 单纯形法这一通用方法， 非线性规划目前还没有适于各种问题的一般算法， 各个方法都有自己特定的适用范围。下面通过实例归纳出非线性规划数学模型的一般形式， 介绍有关非线性规划的基本概念。例1 (投资决策问题)某企业有 n个项目可供选择投资，并且至少要对其中一个项目投资。已知该企业拥有总资金 A元，投资于第i(i 1, ,n)个项目需花资金ai元,并预计可收益b元。试选择最佳投资

2、方案。解设投资决策变量为且总的投资金额不能超过总资金 A，故有限制条件n0 ai xi Ai 1另外，由于xi(i 1, ,n)只取值0或1，所以还有xi(1 xi) 0, i 1, ,n.最佳投资方案应是投资额最小而总收益最大的方案， 所以这个最佳投资决策问题归结为总资金以及决策变量(取 0或1)的限制条件下，极大化总收益和总投资之比。因此，其数学模型为：ni 1 max Q n aiXii 1ns.t. 0 aixi Ai 1Xi(1 Xi) 0, i 1, ,n.上面例题是在一组等式或不等式的约束下，求一个函数的最大值(或最小值)问题，其中目标函数或约束条件中至少有一个非线性函数， 这类

3、问题称之为非线性规划问题，简记为(NP。可概括为一般形式min f (x)s.t. hj(x) 0, j 1, ,q (NP)gi(x) 0, i 1, , p其中x Xi Xn T称为模型(NF3)的决策变量，f称为目标函数，gi (i 1, , p)和 hj( j 1, ,q) 称为约束函数。另外， gi(x) 0 (i 1, ,p) 称为等式约束，hj(x) 0 (j 1, ,q) 称为不等式约束。对于一个实际问题，在把它归结成非线性规划问题时，一般要注意如下几点：( i )确定供选方案：首先要收集同问题有关的资料和数据，在全面熟悉问题的基 础上，确认什么是问题的可供选择的方案，并用一组

4、变量来表示它们。( ii )提出追求目标：经过资料分析，根据实际需要和可能，提出要追求极小化 或极大化的目标。并且，运用各种科学和技术原理，把它表示成数学关系式。(iii )给出价值标准：在提出要追求的目标之后，要确立所考虑目标的“好”或 “坏”的价值标准，并用某种数量形式来描述它。( iv )寻求限制条件：由于所追求的目标一般都要在一定的条件下取得极小化或 极大化效果， 因此还需要寻找出问题的所有限制条件， 这些条件通常用变量之间的一些 不等式或等式来表示。1.2线性规划与非线性规划的区别如果线性规划的最优解存在， 其最优解只能在其可行域的边界上达到 (特别是可行 域的顶点上达到) ；而非线

5、性规划的最优解(如果最优解存在)则可能在其可行域的任 意一点达到。1.3非线性规划的 Matlab 解法Matlab 中非线性规划的数学模型写成以下形式min f (x)Ax BAeq x BeqC(x) 0 ,Ceq(x) 0其中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 定义了线性约束

6、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 求下列非线性规划问题22min f(x) x1 x2 82x1 x2 02x1 x2 2 0x1,x2 0.(i )编写M文件fun 1.mfun ctio n f=fun 1(x);f=x(1F2+x(2)

7、A2+8;和M文件fun2.mfun ctio n g,h=fu n2(x);g=-x(1)A2+x(2);h=-x(1)-x(2)A2+2; % 等式约束(ii )在Matlab的命令窗口依次输入opti on s=optimset;x,y=fmi ncon( fun1 ,ran d(2,1),zeros(2,1), .fun2 , opti ons)就可以求得当x, 1,x2 1时，最小值y 10。1.4求解非线性规划的基本迭代格式记(NP)的可行域为K。若x* K ,并且f (x ) f (x), x K则称x*是(NP)的整体最优解，f(X*)是(NP)的整体最优值。如果有f (x )

8、 f (x), x K, 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)的严格局部最优值。由于线性规划的目标函数为线性函数，可行域为凸集，因而求出的最优解就是整 个可行域上的全局最优解。 非线性规划却不然，有时求出的某个解虽是一部分可行域上的极值点，但并不一定是整个可行域上的全局最优解。对于

9、非线性规划模型(NP)，可以采用迭代方法求它的最优解。迭代方法的基本思 想是：从一个选定的初始点 x Rn出发，按照某一特定的迭代规则产生一个点列xk，使得当xk是有穷点列时，其最后一个点是 (NP)的最优解；当xk是无穷点列时，它有极限点，并且其极限点是 (NP)的最优解。设xk Rn是某迭代方法的第k轮迭代点，xk 1 Rn是第k 1轮迭代点，记x x tkP ( 1)这里tk R1, pk Rn，|pk| 1，显然Pk是由点xk与点xk 1确定的方向。式(1 )就 是求解非线性规划模型(NP)的基本迭代格式。通常，我们把基本迭代格式(1)中的pk称为第k轮搜索方向，tk为沿pk方向的 步

10、长，使用迭代方法求解(NP)的关键在于，如何构造每一轮的搜索方向和确定适当的步 长。设x Rn,p 0 ,若存在 0，使f(x tp) f(x), t (0,)，称向量p是f在点x处的下降方向。设x Rn, p 0 ,若存在t 0，使x tp K ,称向量p是点x处关于K的可行方向。一个向量p，若既是函数f在点x处的下降方向，又是该点关于区域 K的可行方 向，则称之为函数 f在点x处关于K的可行下降方向。现在，我们给出用基本迭代格式(1)求解(NP)的一般步骤如下：0选取初始点x0，令k: 0。k1 构造搜索方向，依照一定规则，构造 f在点x处关于K的可行下降方向作为搜索方向pk。k k2寻求

11、搜索步长。以x为起点沿搜索方向 p寻求适当的步长tk，使目标函数值 有某种意义的下降。3求出下一个迭代点。按迭代格式( 1)求出xk1k kx tkP。若xk 1已满足某种终止条件，停止迭代。k 1 k4以xk 1代替xk，回到1。步。1.5凸函数、凸规划设f (x)为定义在 n维欧氏空间 E(n)中某个凸集 R上的函数，若对任何实数(0 1)以及R中的任意两点x和x(2)，恒有f ( x(1 )x(2) f(x)(1 )f (x(2)则称f (x)为定义在R上的凸函数。若对每一个 (0 1)和x(1) x(2) R恒有f ( x(1) (1 )x(2) f(x)(1 )f (x(2)则称f(

12、x)为定义在R上的严格凸函数。考虑非线性规划min f (x)x RR x|gj(x) 0,j 1,2, ,1假定其中f (x)为凸函数，gj(x)( j 1,2, ,l)为凸函数，这样的非线性规划称为凸规划。可以证明，凸规划的可行域为凸集，其局部最优解即为全局最优解，而且其最优 解的集合形成一个凸集。当凸规划的目标函数 f (x)为严格凸函数时，其最优解必定唯一(假定最优解存在)。由此可见，凸规划是一类比较简单而又具有重要理论意义的非 线性规划。无约束问题2.1一维搜索方法当用迭代法求函数的极小点时， 常常用到一维搜索，即沿某一已知方向求目标函数的极小点。一维搜索的方法很多， 常用的有：(1

13、)试探法(“成功一失败”，斐波那契法, 0.618法等)；插值法(抛物线插值法，三次插值法等) ；(3)微积分中的求根法(切线法，二分法等)。考虑一维极小化问题min f(t) (2)a t b若f (t)是a,b区间上的下单峰函数，我们介绍通过不断地缩短 a,b的长度，来搜索得(2)的近似最优解的两个方法。为了缩短区间a,b，逐步搜索得(2)的最优解t*的近似值，我们可以采用以下 途径：在a,b中任取两个关于a, b是对称的点t1和t2 (不妨设t2 t1，并把它们叫 做搜索点)，计算f(tj和f(t2)并比较它们的大小。对于单峰函数，若 f(t2) f (ti)，则必有t* a,ti，因而

14、a,ti是缩短了的单峰区间；若f(tj f (t2)，则有t* t2,b，故t2,b是缩短了的单峰区间；若 f(t2) f (ti)，则a,ti和t2,b都是缩短了的单峰。因此通过两个搜索点处目标函数值大小的比较， 总可以获得缩短了的单峰区间。对于新的单峰区间重复上述做法，显然又可获得更短的单峰区间。如此进行， 在单峰区间缩短到充分小时，我们可以取最后的搜索点作为( 2)最优解的近似值。应该按照怎样的规则来选取探索点，使给定的单峰区间的长度能尽快地缩短？2.1.1Fib on acci 法若数列 Fn满足关系：Fo Fi 1Fn Fn 2 Fn 1, n 2,3,则称Fn为Fib on acc

15、i数列，Fn称为第n个Fibo nacci数，称相邻两个 Fibo nacci数之为 Fibonacci 分数。当用斐波那契法以 n个探索点来缩短某一区间时，区间长度的第一次缩短率为弘，其后各次分别为 弘，也，1。由此，若t1和t2(t2 t1)是单峰区间a,bFn Fn Fn 2 F2中第1个和第2个探索点的话，那么应有比例关系如果要求经过一系列探索点搜索之后，使最后的探索点和最优解之间的距离不超 过精度 0,这就要求最后区间的长度不超过 ，即据此，我们应按照预先给定的精度 ，确定使(4)成立的最小整数n作为搜索次数,直到进行到第n个探索点时停止。用上述不断缩短函数 f(t)的单峰区间a,b

16、的办法，来求得问题(2)的近似解,是Kiefer(1953年)提出的，叫做 Fin bo nacci法，具体步骤如下：索次数n。2 k 1,aa,b b0，计算最初两个搜索点，按(3)计算t1和t2。3 while kn1f1f (t1), f2f(t2)iff1f2F(n1 k)心a t2;t2匕恙a(ba)F(nk)elseF(n 1k) /b忧t1t2 ;t2b(ab)F(nk)endk k 1end4当进行至k n 1时，1t1 t2 -(a b)2这就无法借比较函数值 f (tj和f(t2)的大小确定最终区间，为此，取1t2 (a b)21ti a (- )(b a)2其中 为任意小

17、的数。在t-和t2这两点中，以函数值较小者为近似极小点， 相应的函数值为近似极小值。并得最终区间 a,t1或t2,b。由上述分析可知，斐波那契法使用对称搜索的方法， 逐步缩短所考察的区间，它能以尽量少的函数求值次数，达到预定的某一缩短率。例3试用斐波那契法求函数 f(t) t2 t 2的近似极小点，要求缩短后的区间不大于区间1,3的0.08倍。程序留作习题。2.1.20.618 法若 0 ,满足比例关系1145 1称之为黄金分割数，其值为 0.61803398872黄金分割数.。Fn 11 FnFn 1Fn和Fib on acci分数之间有着重要的关系，它们是，n为偶数，Fn 1n ， n为奇

18、数。F n 1Fn 1lim n F n现用不变的区间缩短率 0.618，代替斐波那契法每次不同的缩短率， 就得到了黄金分割法(0.618法)。这个方法可以看成是斐波那契法的近似，实现起来比较容易，效果 也相当好，因而易于为人们所接受。用0.618法求解，从第2个探索点开始每增加一个探索点作一轮迭代以后，原单峰区间要缩短0.618倍。计算n个探索点的函数值可以把原区间 a0,b0连续缩短n 1次，因为每次的缩短率均为 ，故最后的区间长度为时，可用下式计算探索点个数 n :(bo ao)这就是说，当已知缩短的相对精度为n 1当然，也可以不预先计算探索点的数目 n ,而在计算过程中逐次加以判断，看

19、是否已满足了提出的精度要求。0.618法是一种等速对称进行试探的方法，每次的探索点均取在区间长度的 0.618倍和0.382倍处。2.2二次插值法对极小化问题(2)，当f(t)在a,b上连续时，可以考虑用多项式插值来进行一 维搜索。它的基本思想是：在搜索区间中，不断用低次(通常不超过三次)多项式来近 似目标函数，并逐步用插值多项式的极小点来逼近( 2)的最优解。2.3无约束极值问题的解法无约束极值问题可表述为min f (x), x E(n) (5)求解问题(5)的迭代法大体上分为两种：一是用到函数的一阶导数或二阶导数， 称为解析法。另一是仅用到函数值，称为直接法。2.3.1解析法2.3.1.

20、1梯度法(最速下降法)对基本迭代格式k 1 k 丄 kX X tkP ( 6)k k我们总是考虑从点 X出发沿哪一个方向p，使目标函数f下降得最快。微积分的知识 告诉我们，点X k的负梯度方向k k、P f(x )，是从点Xk出发使f下降最快的方向。为此，称负梯度方向 f(xk)为f在点Xk处的最速下降方向。按基本迭代格式(6)，每一轮从点Xk出发沿最速下降方向 f(xk)作一维搜索，来建立求解无约束极值问题的方法，称之为最速下降法。这个方法的特点是，每轮的搜索方向都是目标函数在当前点下降最快的方向。同时，用f(xk) 0或f(xk)| 作为停止条件。其具体步骤如下：1 选取初始数据。选取初始

21、点 X0，给定终止误差，令 k: 0。2求梯度向量。计算 f(xk)，若| f(xk)| ，停止迭代，输出xk。否则，进行3。构造负梯度方向。取pk f(xk).进行一维搜索。求tk，使得f (xk tkpk) min f (xk tpk)t 0令xk 1例4xk tk pk, k : k 1,转 2 。用最速下降法求解无约束非线性规划问题min f (x) x； 25x；10 6。其中x (X1,X2)T，要求选取初始点 X0 (2,2)t，终止误差解：(i) f (x) (2x1,50x2)T编写M文件detaf.m如下function f,df=detaf(x);f=x(1)A2+25*

22、x(2)A2;df(1)=2*x(1); df(2)=50*x(2);(ii )编写M文件zuisu.mclcx=2;2;f0,g=detaf(x);while norm(g)0.000001p=-g/norm(g);t=1.0;f=detaf(x+t*p);while ff0t=t/2;f=detaf(x+t*p);endx=x+t*pf0,g=detaf(x)end2.3.1.2Newton 法处的二次逼近式考虑目标函数f在点xf(x) Q(x) f(xk)f (xk)T(x xk) 扣 xk)T 2f (xk)(x xk)假定Hesse阵编写M文件nwfun.m如下： function

23、f,df,d2f=nwfun( x);f=x(1)A4+25*x(2)A4+x(1)A2*x(2)A2; df(1)=4*x(1)A3+2*x(1)*x(2)A2; df(2)=100*x (2) A3+2*x(1)A2*x(2); d2f(1,1)=12*x(1)A2+2*x(2)A2; d2f(1,2)=4*x(1)*x(2); d2f(2,1)=d2f(1,2); d2f(2,2)=300*x (2)A2+4*x(1)*x(2);(ii)编写M文件：clcx=2;2;f0,g1,g2=nwfun(x)while norm(g1)0.00001 %dead loop,for i=1:3p=

24、-i nv(g2)*g1,p=p/no rm(p)t=1.0,f=detaf(x+t*p)while ff0t=t/2,f=detaf(x+t*p),end x=x+t*pf0,g1,g2=nwfun(x)end如果目标函数是非二次函数，一般地说，用 Newt on法通过有限轮迭代并不能保证可求得其最优解。Newton法的优点是收敛速度快；缺点是有时不好用而需采取改进措施，此外，当 维数较高时，计算 2f(xk) 1的工作量很大。2.3.1.3 变尺度法变尺度法(Variable Metric Algorithm )是近20多年来发展起来的，它不仅是求解无约束极值问题非常有效的算法， 而且也已

25、被推广用来求解约束极值问题。 由于它既避免了计算二阶导数矩阵及其求逆过程， 又比梯度法的收敛速度快，特别是对高维问题具有显著的优越性，因而使变尺度法获得了很高的声誉。 下面我们就来简要地介绍一种变尺度法一DFP法的基本原理及其计算过程。这一方法首先由 Davidon在1959年提出，后经Fletcher和Powell加以改进。我们已经知道，牛顿法的搜索方向是 2f(xk)导数矩阵2f(xk)及其逆阵，我们设法构造另一个矩阵， 的逆阵2f(xk) 1，这一类方法也称拟牛顿法(下面研究如何构造这样的近似矩阵，并将它记为 H (k)。我们要求：每一步都能以现有的信息来确定下一个搜索方向； 每做一次选

26、代，目标函数值均有所下降； 这些近似矩阵最后应收敛于解点处的当f (x)是二次函数时， f (x )1 f(xk)，为了不计算二阶用它来逼近二阶导数矩阵Quasi-Newton Method)。Hesse阵的逆阵。k k 1其Hesse阵为常数阵A，任两点x和x 处的梯度之差为k k 1f(x ) A(xk 1 kx x A对于非二次函数，仿照二次函数的情形, 矩阵H(k 1)满足关系式k 1 kx x这就是常说的拟若令G(k)kx则式(7)变为kx1 f(xk1)f(xk)要求其H(k 1) f(xk 1) f (xk)Newto n 条件。xkf(xk 1) f(xk)1 xkHesse阵

27、的逆阵的第k 1次近似(7)(8)H(k 1) G(k)_ ( 9)1已知，用下式求 H化1)(设H和H 1)均为对称正定阵)； _ H(k) 其中 H(k)称为第现假定H(k)H(k1) IH(k)k次校正矩阵。显然,xk (H(k)_ ( 10)H (k 1应满足拟Newton条件(9)，即要求H (k) G(k)由此可以设想,H(k) G (k)H (k)的一种比较简单的形式是 xk(Q(k) )T H(k) G(k)(W(k)Txk H(k) G(k)(11)H(k)其中Q(k)和W(k)为两个待定列向量。将式(12)中的 H(k)代入(11)，得xk(Q(k)T G(k) H(k)

28、G(k)(W(k)T G(k) 这说明，应使xk H(k) G(k) (Q(k)T G(k)(W(k)T G(k) 1Q(k)W(k) kH(k) G(k)由式(13)得k( xk)T G(k) k( G(k)TH(k) G(k) 1若(xk)T G(k)和(G(k)TH(k)G(k)不等于零，则有上述矩阵称为尺度矩阵。通常，我们取第一个尺度矩阵 H (0)为单位阵，以后的尺度矩阵按式(18)逐步形成。可以证明：(i)当xk不是极小点且 H(k)正定时，式(17)右端两项的分母不为零，从而可 按式(18)产生下一个尺度矩阵 H (k 1)；(ii )若H(k)为对称正定阵，则由式(18)产生的

29、H (k 1)也是对称正定阵； iii )由此推出DFP法的搜索方向为下降方向。现将DFP变尺度法的计算步骤总结如下。1 给定初始点X0及梯度允许误差 0。2。若 f(x0) ，则X0即为近似极小点，停止迭代，否则，转向下一步。3。令H(0) I (单位矩阵)，0 (0) 0 p H f(x )在p0方向进行一维搜索，确定最佳步长0 0 0 0 min f (x p ) f(x 0P )如此可得下一个近似点1 0 0x x p4 一般地，设已得到近似点 xk，算出 f(Xk)，若I f(xk)| _则Xk即为所求的近似解，停止迭代；否则，计算 H(k):H(k 1) G(k 1)(G(k 1)T h (k 1)(G(k 1)T H(k 1) G(k 1)并令pk H (k) f(xk)，在pk方向上进行一维搜索，得 k，从而可得下一个近似点k 1 k kX X kP5若Xk 1满足精度要求，则 Xk 1即为所求的近似解，否则，转回 4。，直到求出某点满足精度要求为止。232直接法在无约束