第三章 非线性规划

《第三章 非线性规划》由会员分享，可在线阅读，更多相关《第三章 非线性规划（21页珍藏版）》请在装配图网上搜索。

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

2、i解 设投资决策变量为i = 1,，n，11,决定投资第i个项目i o,决定不投资第i个项目则投资总额为 ax，投资总收益为工bx。因为该公司至少要对一个项目投资，并i ii ii=1i=1且总的投资金额不能超过总资金A，故有限制条件0 乙 a x Aiii =1另外，由于x （i = 1,n）只取值0或1,所以还有ix （1 一 x ） = 0, i = 1, ,n.ii最佳投资方案应是投资额最小而总收益最大的方案，所以这个最佳投资决策问题归 结为总资金以及决策变量（取0或1）的限制条件下，极大化总收益和总投资之比。因 此，其数学模型为： b xmax Q =ii-4=1 a xiii=1s

3、.t. 0 a x Aiii =1x （1 一 x ） = 0, i = 1,n.ii面例题是在一组等式或不等式的约束下，求一个函数的最大值（或最小值）问题，其中至少有一个非线性函数，这类问题称之为非线性规划问题。可概括为一般形式s.t.其中x = x1(NP)min f （ x）h .(x) 0, j = 1,qjg (x) = 0, i = 1,pix T称为模型（NP）的决策变量，f称为目标函数，g （i = h,P） ni和h.(j = 1,q)称为约束函数。另外，g.(X) = 0 (i = 1,p)称为等式约束， jih.(x) 0 (j = 1,q)称为不等式的约束。对于一个实际

4、问题，在把它归结成非线性规划问题时，一般要注意如下几点：(i) 确定供选方案：首先要收集同问题有关的资料和数据，在全面熟悉问题的基 础上，确认什么是问题的可供选择的方案，并用一组变量来表示它们。(ii) 提出追求目标：经过资料分析，根据实际需要和可能，提出要追求极小化 或极大化的目标。并且，运用各种科学和技术原理，把它表示成数学关系式。(iii) 给出价值标准：在提出要追求的目标之后，要确立所考虑目标的“好”或 “坏”的价值标准，并用某种数量形式来描述它。(iv) 寻求限制条件：由于所追求的目标一般都要在一定的条件下取得极小化或 极大化效果，因此还需要寻找出问题的所有限制条件，这些条件通常用变

5、量之间的一些 不等式或等式来表示。1.2 线性规划与非线性规划的区别 如果线性规划的最优解存在，其最优解只能在其可行域的边界上达到(特别是可行域的顶点上达到)；而非线性规划的最优解(如果最优解存在)则可能在其可行域的任 意一点达到。1.3 非线性规划的 Matlab 解法Matlab 中非线性规划的数学模型写成以下形式min f ( x)Ax BAeq - x = BeqVC (x) 0Ceq( x) = 0其中f (x)是标量函数，A,B, Aeq,Beq是相应维数的矩阵和向量，C(x), Ceq(x)是非 线性向量函数。Matlab 中的命令是 X=FMINCON(FUN,X0,A,B,A

6、eq,Beq,LB,UB,NONLCON,OPTIONS)它的返回值是向量x，其中FUN是用M文件定义的函数f (x) ； X0是x的初始值； A,B,Aeq,Beq定义了线性约束A * X 0V 1 2x x 2 + 2 = 01 2x , x 0.12(i)编写M文件funl.mfunction f=fun1(x);f=x(1)入2+x(2)入2+8;和M文件fun2.mfunctiong,h=fun2(x);g=-x(1)入2+x(2);h=-x(1)-x(2)入2+2; %等式约束(ii)在Matlab的命令窗口依次输入 options=optimset(largescale,off)

7、;x,y=fmincon(fun1,rand(2,1),zeros(2,1), . fun2, options)就可以求得当x = 1, x = 1时，最小值y = 10。121.4 求解非线性规划的基本迭代格式记(NP)的可行域为K。若x* eK，并且f (x*) f (x), Vx e K则称x*是(NP)的整体最优解，f(x*)是(NP)的整体最优值。如果有f (x*) f (x), Vx e K,x 丰 x*则称x*是(NP)的严格整体最优解，f (x*)是(NP)的严格整体最优值。若x* e K，并且存在x*的邻域N&(x*)，使f (x*) f (x), Vx e N (x*) A

8、 K，8则称x*是(NP)的局部最优解，f(x*)是(NP)的局部最优值。如果有f (x*) 0，使f (元 + tp) 0，使X + tp g K ,称向量p是点X处关于K的可行方向。一个向量P ,若既是函数f在点X处的下降方向，又是该点关于区域K的可行方 向，称之为函数f在点X处关于K的可行下降方向。现在，我们给出用基本迭代格式(1)求解(NP)的一般步骤如下：0选取初始点x0，令k := 0。1构造搜索方向，依照一定规划，构造f在点xk处关于K的可行下降方向作为 搜索方向 pk 。2寻求搜索步长。以xk为起点沿搜索方向pk寻求适当的步长t，使目标函数值 k 有某种意义的下降。3 求出下一

9、个迭代点。按迭代格式(1)求出Xk+1 = Xk + t pk。k若 Xk+1 已满足某种终止条件，停止迭代。4以Xk+1代替Xk，回到1步。1.5 凸函数、凸规划设f (X)为定义在n维欧氏空间E (n)中某个凸集R上的函数，若对任何实数 a (0 a 1)以及R中的任意两点x和x，恒有f (ax(1)+ (1 -a)x(2) af (x)+ (1 -a)f (x) 则称f ( x )为定义在R上的凸函数。若对每一个a(0a 1)和x(1)丰xg R恒有 f(ax(1)+(1-a)x(2)af(x(1)+(1-a)f(x(2) 则称f ( x )为定义在R上的严格凸函数。考虑非线性规划min

10、 f ( x ) xgR、R = x I g (x) 0, j = 1,2,1假定其中f (X)为凸函数，g.(X)(j二1,2,1)为凸函数，这样的非线性规划称为 凸规划。可以证明，凸规划的可行域为凸集，其局部最优解即为全局最优解，而且其最优 解的集合形成一个凸集。当凸规划的目标函数f (X)为严格凸函数时，其最优解必定唯 一(假定最优解存在)。由此可见，凸规划是一类比较简单而又具有重要理论意义的非 线性规划。2 无约束问题2.1 一维搜索方法 当用迭代法求函数的极小点时，常常用到一维搜索，即沿某一已知方向求目标函数 的极小点。一维搜索的方法很多，常用的有：(1)试探法(“成功失败”，斐波那

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

12、t ，因而a, t 是缩短了的单峰区间；若f (t ) f (t )，则有1 1 1 2t* g t ,b，故t ,b是缩短了的单峰区间；若f (t ) = f (t )，则a,t 和t ,b都是 2 2 2 1 1 2 缩短了的单峰。因此通过两个搜索点处目标函数值大小的比较，总可以获得缩短了的单 峰区间。对于新的单峰区间重复上述做法，显然又可获得更短的单峰区间。如此进行， 在单峰区间缩短到充分小时，我们可以取最后的搜索点作为(2)最优解的近似值。应该按照怎样的规则来选取探索点，使给定的单峰区间的长度能尽快地缩短？2.1.1 Fibonacci 法如用F表示计算n个函数值能缩短为单位长区间的最

13、大原区间长度，可推出F满 nn 足关系F 二 F 二 101F = F + F , n = 2,3, ,nn2n 1数列F 称为Fibonacci数列，F称为第n个Fibonacci数，相邻两个Fibonacci数之nnF比n-1称为Fibonacci分数。Fn当用斐波那契法以n个探索点来缩短某一区间时，区间长度的第一次缩短率为FF F FpT1，其后各次分别为一n-1, n-3，,一1。由此，若t和t (t 0，这就要求最后区间的长度不超过，即ba 0，由(4)确定搜 00索次数n。2 k 1,a a ,b b，计算最初两个搜索点，按(3)计算t和t。0 0 1 23 while k n 1

14、 f (t2)FF(T (b -。)f1 f(t1), f21 1 2 if f 0，满足比例关系 1 一1 5 一 1称之为黄金分割数，其值为 2 o,6180339887。厶黄金分割数和Fibonacci分数之间有着重要的关系，它们是FF1n-1 W W L，FFnn-1F2O =lim。Fn T8 人n 为偶数n 为奇数。n现用不变的区间缩短率 0.618，代替斐波那契法每次不同的缩短率，就得到了黄金分割法（0.618 法）。这个方法可以看成是斐波那契法的近似，实现起来比较容易，效果 也相当好，因而易于为人们所接受。用 0.618 法求解，从第 2 个探索点开始每增加一个探索点作一轮迭代

15、以后，原单峰区间要缩短0.618倍。计算n个探索点的函数值可以把原区间a ,b 连续缩短n -1 00次，因为每次的缩短率均为卩，故最后的区间长度为（b - a ）卩 n-100这就是说，当已知缩短的相对精度为时，可用下式计算探索点个数n : 卩n-1 当然，也可以不预先计算探索点的数目n，而在计算过程中逐次加以判断，看是否已满 足了提出的精度要求。0.618 法是一种等速对称进行试探的方法，每次的探索点均取在区间长度的 0.618 倍和 0.382倍处。2.2 二次插值法对极小化问题（2），当f （t）在a,b上连续时，可以考虑用多项式插值来进行一 维搜索。它的基本思想是:在搜索区间中，不断

16、用低次（通常不超过三次）多项式来近 似目标函数，并逐步用插值多项式的极小点来逼近（2）的最优解。2.3 无约束极值问题的解法 无约束极值问题可表述为min f （x）, x g E（n）（5）求解问题（5）的迭代法大体上分为两点:一是用到函数的一阶导数或二阶导数， 称为解析法。另一是仅用到函数值，称为直接法。2.3.1 解析法2.3.1.1 梯度法（最速下降法） 对基本迭代格式Xk+1 = Xk + t pk（6）k我们总是考虑从点Xk出发沿哪一个方向pk，使目标函数f下降得最快。微积分的知识 告诉我们，点 xk 的负梯度方向Pk = -f （Xk ），是从点Xk出发使f下降最快的方向。为此，

17、称负梯度方向-Yf （Xk ）为f在点xk处的 最速下降方向。按基本迭代格式（6）,每一轮从点Xk出发沿最速下降方向（Xk）作一维搜索， 来建立求解无约束极值问题的方法，称之为最速下降法。这个方法的特点是，每轮的搜索方向都是目标函数在当前点下降最快的方向。同 时，用P 3 ) = 0或|Vf (xk) 作为停止条件。其具体步骤如下：1 选取初始数据。选取初始点X0，给定终止误差，令k := 0。2求梯度向量。计算Vf (xk)， 若 |Vf (xk )| 0令 xk+1 = xk + tpk , k := k +1,转 2。k例 4 用最速下降法求解无约束非线性规划问题min f (x) =

18、x 2 + 25 x 212其中x二(x ,x )T，要求选取初始点x0二(2,2)T。12解: (i) Vf (x)二(2x ,50x )t12编写 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.mclcx=2;2;f0,g=detaf(x);while norm ( g ) 0.000001p=-g/norm(g);t=1.0;f=detaf(x+t*p);whileff0t=t/2;f=detaf(x+t*p);endx=x+t*pf0,

19、g=detaf(x)end2.3.1.2 Newton 法考虑目标函数f在点xk处的二次逼近式f ( x)Q ( x) = f ( xk ) + Vf ( xk )t ( x 一 xk ) + 2( x 一 xk )t V 2 f ( xk )( x 一 xk )假定 Hesse 阵Qx212f (Xk )x QxV2 f (xk)=1nQ2 f (xk)Qx Qxn1Q2 f (xk)Qx 2n正定。由于V2/(xk)正定，函数Q的驻点xk+i是Q(x)的极小点。为求此极小点，令VQ(xk+1) = Vf (xk ) + V2 f (xk )(xk+1 一 xk ) = 0 ,即可解得xk+

20、1 = xk - V2 f (xk )-1 Vf (xk ).对照基本迭代格式(1),可知从点xk出发沿搜索方向。pk 二-V 2 f (xk )-1 Vf (xk)并取步长t二1即可得Q(x)的最小点Xk+1。通常，把方向pk叫做从点Xk出发的 kNewton 方向。从一初始点开始，每一轮从当前迭代点出发，沿 Newton 方向并取步长 为 1 的求解方法，称之为 Newton 法。其具体步骤如下：1选取初始数据。选取初始点x0，给定终止误差 0，令k := 0。2求梯度向量。计算Vf (xk)，若|Vf (xk )| 0.00001p=-inv(g2)*g1x=x+pf0,g1,g2=nw

21、fun(x)end如果目标函数是非二次函数，一般地说，用Newton法通过有限轮迭代并不能保证 可求得其最优解。为了提高计算精度，我们在迭代时可以采用变步长计算上述问题，程序如下： clcx=2;2;f0,g1,g2=nwfun(x)while norm(g1)0.00001p=-inv(g2)*g1,p=p/norm(p) t=1.0,f=nwfun(x+t*p) while ff0t=t/2,f=nwfun(x+t*p),endx=x+t*pf0,g1,g2=nwfun(x)endNewton 法的优点是收敛速度快；缺点是有时不好用而需采取改进措施，此外，当 维数较高时，计算-V 2 f

22、(xk)-1的工作量很大。2.3.1.3 变尺度法变尺度法(Variable Metric Algorithm)是近20多年来发展起来的，它不仅是求解 无约束极值问题非常有效的算法，而且也已被推广用来求解约束极值问题。由于它既避 免了计算二阶导数矩阵及其求逆过程，又比梯度法的收敛速度快，特别是对高维问题具 有显著的优越性，因而使变尺度法获得了很高的声誉。下面我们就来简要地介绍一种变 尺度法一DFP法的基本原理及其计算过程。这一方法首先由Davidon在1959年提出， 后经 Fletcher 和 Powell 加以改进。我们已经知道，牛顿法的搜索方向是 -V2 f(xk)-1Vf(xk)，导数

23、矩阵V2f (xk)及其逆阵，我们设法构造另一个矩阵，用它来逼近二阶导数矩阵 的逆阵V2 f (xk )-1，这一类方法也称拟牛顿法(Quasi-Newton Method)。下面研究如何构造这样的近似矩阵，并将它记为H(k)。我们要求：每一步都能以 现有的信息来确定下一个搜索方向；每做一次选代，目标函数值均有所下降；这些近似 矩阵最后应收敛于解点处的Hesse阵的逆阵。当/(x)是二次函数时，其Hesse阵为常数阵A，任两点xk和xk+1处的梯度之差为Vf ( xk+1) -Vf ( xk ) = A( xk+1 一 xk )或xk+1 一 xk = A-1Vf ( xk+1) -Vf (

24、xk )对于非二次函数，仿照二次函数的情形，要求其Hesse阵的逆阵的第k +1次近似 矩阵H (k+1)满足关系式xk +1 - xk = H (k +DVf ( xk+1) -Vf ( xk )(7)这就是常说的拟Newton条件。若令AG( k)=Vf (xk+i)-Vf (xk)Axk = Xk+1 Xk8)则式(7)变为 _Axk = H(k+i)AG(k),(9)现假定H(k)已知，用下式求H(k+i)(设H(k)和H(k+i)均为对称正定阵)；H (k+1) = H (k) + AH (k)(10)其中AH(k)称为第k次校正矩阵。显然，H(k+1)应满足拟Newton条件(9)

25、,即要求Axk = ( H (k) + AH (k )AG (k)或AH (k) AG (k) = Axk H (k) AG (k) 由此可以设想，AH(k)的一种比较简单的形式是AH (k) = Axk (Q (k) T 一 H (k) AG (k )(W (k) T其中Q (k)和W (k)为两个待定列向量。将式(12)中的AH(k)代入(11),得Axk (Q (k) T AG (k) 一 H (k) AG(k )(W (k) T AG (k) = Axk 一 H (k) AG (k)这说明，应使(Q(k)TAG(k) = (W(k)TAG(k) =1考虑到AH (k)应为对称阵，最简单

26、的办法就是取 JQ(k)=耳 AxkW (k) =g H (k) AG (k)k由式(13)得_耳(Axk )t AG (k) =g ( AG (k )tH (k) AG(k) = 1 k若(Axk ) T AG (k)和(AG (k) T H (k) AG (k)不等于零，则有1 = 1(Axk )T AG(k)( AG(k)T Axk1(AG (k) tH (k) AG (k)于是，得校正矩阵_Axk ( Axk ) TH (k) AG (k )(G (k) T AH (k)AH (k)=一(AG (k) T Axk( AG (k) tH (k) AG(k)11)(12)13)14 )(

27、15 )( 16 )(17)从而得到H (k + 1) = H (k) +Axk (Axk )T(AG(k)TAxkH(k)AG(k)(G(k)TAH(k)(AG (k) tH (k) AG (k)上述矩阵称为尺度矩阵。通常，我们取第一个尺度矩阵H (0)为单位阵，以后的尺度矩 阵按式(18)逐步形成。可以证明：(i)当xk不是极小点且H(k)正定时，式(17)右端两项的分母不为零，从而可 按式(18)产生下一个尺度矩阵H(k+1)；_(ii) 若H(k)为对称正定阵，则由式(18)产生的H(k+1)也是对称正定阵；(iii) 由此推出DFP法的搜索方向为下降方向。现将DFP变尺度法的计算步骤

28、总结如下。1给定初始点x0及梯度允许误差S 0。32若|Vf(xo)|，则xo即为近似极小点，停止迭代，否则，转向下一步。令H(0)= I (单位矩阵), po 二-H(0)Vf (xo)在p0方向进行一维搜索确定最佳步长九0：min九f (x0 + Xp0)= f (x0 + 九 p0)0如此可得下一个近似点x 1 = x 0 + X p 004 一般地，设已得到近似点xk，算出Vf (xk)，若Vf ( x 0)| 0，令k ：= 0。2进行基本搜索。令y0 ：= xk，依次沿p0,p1,pn-1中的方向进行一维搜索。对应地得到辅助迭代点y1,y2,yn，即yj = yj-1 +1 pj-

29、1j-1f (yj-1 +1 pj-1) = min f (yj-1 + tpj-1)j = 1,n输出 xk+1 = yn。丿Tt 01|，停止迭代，3构造加速方向。令pn = yn -y0，若p”|8否则进行 4。4确定调整方向。按下式f (ymT) - f (ym)二 max f (yj-1) - f (yj)11 j n 找出m。若f (yo)-2f (yn) + f (2yn - yo) 0p0, p1, p n + 1 ：= p 0,，pm-1, p m +1, , p n -1, pn，k+1k := k +1，转 2。6不调整搜索方向组。令xk+1 := yn,k := k +

30、1，转2。2.4 Matlab 求无约束极值问题在Matlab工具箱中，用于求解无约束极值问题的函数有fminunc和fminsearch,用 法介绍如下。求函数的极小值min f (x),x其中 x 可以为标量或向量。Matlab 中 fminunc 的基本命令是X,FVAL,EXITFLAG,OUTPUT,GRAD,HESSIAN=FMINUNC(FUN,X0,OPTIONS,P1,P2, .)其中的返回值X是所求得的极小点，FVAL是函数的极小值，其它返回值的含义参见相 关的帮助。FUN是一个M文件，当FUN只有一个返回值时，它的返回值是函数f (x)； 当FUN有两个返回值时，它的第二

31、个返回值是f (x)的梯度向量；当FUN有三个返回 值时，它的第三个返回值是f (x)的二阶导数阵(Hessian阵)。X0是向量X的初始值， OPTIONS是优化参数，可以使用确省参数。Pl, P2是可以传递给FUN的一些参数。例6求函数f (x) = 100(x - x 2)2 + (1 - x )2的最小值。2 1 1解：编写 M 文件 fun2.m 如下：function f,g=fun2(x);f=100*(x(2)_x(1)W2+(1_x(1)入2; g=400*x(1)*(x(2)x(1)入2)2*(1x(1);200*(x(2)x(1)入2);在Matlab命令窗口输入opti

32、ons = optimset(GradObj,on);fminunc(fun2,rand(1,2),options)即可求得函数的极小值。在求极值时，也可以利用上二阶导数，编写M文件fun3.m如下：function f,df,d2f=fun3(x);f=100*(x(2)-x(1)A2)A2+(1-x(1)A2;df=-400*x(l)*(x(2)-x(l)人2)-2*(l-x(l);200*(x(2)-x(l)人2); d2f=-400*x(2)+1200*x(l)A2+2,-400*x(l)-400*x(1),200;在Matlab命令窗口输入options = optimset(Gra

33、dObj,on,Hessian,on); fminunc(fun3,rand(1,2),options)即可求得函数的极小值。求多元函数的极值也可以使用 Matlab 的 fminsearch 命令，其使用格式为： X,FVAL,EXITFLAGOUTPUT=FMINSEARCH(FUN,XO,OPTIONS,P1,P2,.) 例7求函数f (x) = sin(x) + 3取最小值时的x值。解 编写f (x)的M文件fun4.m如下：function f=fun4(x);f=sin(x)+3;在命令窗口输入x0=2;x,y=fminsearch(fun4,x0)即求得在初值 2 附近的极小点及

34、极小值。3 约束极值问题 带有约束条件的极值问题称为约束极值问题，也叫规划问题。 求解约束极值问题要比求解无约束极值问题困难得多。为了简化其优化工作，可采 用以下方法：将约束问题化为无约束问题；将非线性规划问题化为线性规划问题，以及 能将复杂问题变换为较简单问题的其它方法。库恩塔克条件是非线性规划领域中最重要的理论成果之一，是确定某点为最优点 的必要条件，但一般说它并不是充分条件(对于凸规划，它既是最优点存在的必要条件， 同时也是充分条件)。3.1 二次规划若某非线性规划的目标函数为自变量x的二次函数，约束条件又全是线性的，就称 这种规划为二次规划。Matlab 中二次规划的数学模型可表述如下

35、：min 2 xtHx + f tx,s.t. Ax b这里h是实对称矩阵，f,b是列向量，a是相应维数的矩阵。Matlab中求解二次规划的命令是X,FVAL= QUADPROG(H,f,A,b,Aeq,beq,LB,UB,X0,OPTIONS)X的返回值是向量x，FVAL的返回值是目标函数在X处的值。(具体细节可以参看在 Matlab 指令中运行 help quadprog 后的帮助)。例 8 求解二次规划min f (x) = 2x 2-4x x + 4x 2-6x -3x1 1 2 2 1 2x + x 3 1 24 x + x 012解 编写如下程序：h=4,-4;-4,8;f=-6;

36、-3;a=1,1;4,1;b=3;9;x,value=quadprog(h,f,a,b,zeros(2,1) 求得1.9500x 二,Min f (x) = 11.0250。1.05003.2 罚函数法 利用罚函数法，可将非线性规划问题的求解，转化为求解一系列无约束极值问题 因而也称这种方法为序列无约束最小化技术，简记为 SUMT (Sequential Unconstrained Minization Technique)。罚函数法求解非线性规划问题的思想是，利用问题中的约束函数作出适当的罚函 数，由此构造出带参数的增广目标函数，把问题转化为无约束非线性规划问题。主要有 两种形式，一种叫外罚

37、函数法，另一种叫内罚函数法，下面介绍外罚函数法。考虑(NL)问题：min f ( x)g (x) 0, j = 1,s,jk (x) = 0, m = 1,tm取一个充分大的数M 0，构造函数P(x, M) = f (x) + M 工 max(g (x),0) 一 M 工 min(h (x),0) +M 工 | k (x) |iiii =1i =1i =1G (x)k0(或 P(x, M) = f (x) + Msuml maxii=1(H (x)、一 Msum min k 0(x )1min+ M II K(x)11、 丿丿k这里 G(x) = Lg (x)g (x)J, H (x) = h

38、 (x)h1r1sK(x) = k (x) k (t)，Matlab中可以直接利用max、1t增广目标函数P (x, M)为目标函数的无约束极值问题丿丿和 sum 函数。)则以minP(x,M )的最优解x也是原问题的最优解。例 9 求下列非线性规划min f (x) = x 2 + x 2 + 81 2x 2 一 x 0 0.12解(i)编写M文件test.mfunction g=test(x);M=50000;f=x(1)入2+x(2)入2+8;g=f-M*min(x(1),0)-M*min(x(2),0)-M*min(x(1)入2-x(2),0). +M*abs(-x(1)-x(2)入2

39、+2);或者是利用Matlab的求矩阵的极小值和极大值函数编写test.m如下： function g=test(x); M=50000;f=x(1)入2+x(2)入2+8; g=f-M*sum(min(x;zeros(1,2)-M*min(x(1)入2-x(2),0).+M*abs(-x(1)-x(2)入2+2);我们也可以修改罚函数的定义，编写test.m如下： function g=test(x);M=50000;f=x(l)人2+x(2)人2+8; g=f-M*min(min(x),O)-M*min(x(l)人2-x(2),0)+M*abs(-x(l)-x(2)人2+2);(ii)在M

40、atlab命令窗口输入 x,y=fminunc(test,rand(2,l) 即可求得问题的解。3.3 Matlab求约束极值问题在 Matlab 优化工具箱中，用于求解约束最优化问题的函数有： fminbnd、fmincon、 quadprog、fseminf、fminimax, 上面我们已经介绍了函数 fmincon 和 quadprog。3.3.1 fminbnd 函数 求单变量非线性函数在区间上的极小值min f (x), x e x , x xl 2Matlab 的命令为X,FVAL = FMINBND(FUN,xl,x2,OPTIONS)，它的返回值是极小点x和函数的极小值。这里f

41、un是用M文件定义的函数或Matlab中 的单变量数学函数。例10求函数f (x)二(x - 3)2 -1, x e 0,5的最小值。解 编写 M 文件 fun5.m function f=fun5(x); f=(x-3)A2-1;在Matlab的命令窗口输入 x,y=fminbnd(fun5,0,5) 即可求得极小点和极小值。3.3.2 fseminf 函数求min F(x) I C(x) 0, Ceq(x) = 0,PHI(x, w) 0xf A * x Bs.t.彳Aeq * x 二 Beq其中C(x),Ceq(x),PHI(x, w)都是向量函数；w是附加的向量变量，w的每个分量都 限

42、定在某个区间内。上述问题的 Matlab 命令格式为 X=FSEMINF(FUN,X0,NTHETA,SEMINFCON,A,B,Aeq,Beq) 其中FUN用于定义目标函数F(x) ； X0为x的初始值；NTHETA是半无穷约束PHI(x, w)的个数；函数SEMINFCON用于定义非线性不等式约束C(x)，非线性等 式约束Ceq(x)和半无穷约束PHI(x, w)的每一个分量函数，函数SEMINFCON有两个输入参量X和S, S是推荐的取样步长，也许不被使用。例 11 求函数 f (x) = (x 一0.5)2 + (x - 0.5)2 + (x - 0.5)2 取最小值时的 x 值,12

43、3 约束为：K (x, w ) = sin(w x )cos(w x ) 一(w 一 50)2 一 sin(w x ) 一 x 11 1 1 1 1 2 1000 1 1 3 3K (x,w ) = sin(w x ) cos(w x )一2 2 2 2 2 1(w 一 50)2 一 sin(w x ) 一 x 11000 2 2 3 31 w 100， 1 w 10012解 (1 )编写 M 文件 fun6.m 定义目标函数如下： function f=fun6(x,s);f=sum(x-0.5).人2);(2) 编写 M 文件 fun7.m 定义约束条件如下： function c,ceq

44、,k1,k2,s=fun7(x,s);c=;ceq=;if isnan(s(1,1) s=0.2,0;0.2 0;end %取样值 w1=1:s(1,1):100; w2=1:s(2,1):100;%半无穷约束kl=sin(wl*x(l).*cos(wl*x(2)-l/1000*(wl-50).人2-sin(wl*x(3)-x(3)-l; k2=sin(w2*x(2).*cos(w2*x(l)-l/1000*(w2-50).人2-sin(w2*x(3)-x(3)-l; %画出半无穷约束的图形plot(wl,kl,-,w2,k2,+);(3) 调用函数fseminf在 Matlab 的命令窗口输

45、入 x,y=fseminf(fun6,rand(3,l),2,fun7) 即可。3.3.3 fminimax 函数 、求解 min xFiA * x bAeq * x = Beqs.t. C(x) 0Ceq (x) = 0LB x UB其中 F(x) = F (x), F (x)。lm 上述问题的 Matlab 命令为 X=FMINIMAX(FUN,X0,A,B,Aeq,Beq,LB,UB,NONLCON) 例 12求函数族f (x),f (x),f (x),f (x), f (x)取极大极小值时的x值。其中: l2345f (x) = 2x2 + x2 一 48x 一 40x + 3041

46、12 1 2f (x) = 一x 2 一 3 x 22 1 2 f (x) = x + 3x 一 183 12f (x) = 一x - x4 12f (x) = x + x 一 8I 512解 (1 )编写 M 文件 fun8.m 定义向量函数如下 function f=fun8(x);f=2*x(l)人2+x(2)人2-48*x(l)-40*x(2)+304 -x(1)A2-3*x(2)A2x(1)+3*x(2)-18 -x(1)-x(2) x(1)+x(2)-8;(2)调用函数 fminimax x,y=fminimax(fun8,rand(2,1)3.3.4 利用梯度求解约束优化问题例1

47、3已知函数f (x) = exi(4x2 + 2x2 + 4x x + 2x +1)，且满足非线性约束:1 2 1 2 2f x x x x 1012求 min f ( x)x 分析:当使用梯度求解上述问题时，效率更高并且结果更准确。 题目中目标函数的梯度为:ex】(4x2 + 2x2 + 4x x + 8x + 6x +1)1 2 1 2 1 2ex】(4x + 4x + 2)12解 (1)编写 M 文件 fun9.m 定义目标函数及梯度函数: function f,df=fun9(x);f=exp(x(1)*(4*x(1)A2+2*x(2)A2+4*x(1)*x(2)+2*x(2)+1);

48、df=exp(x(1)*(4*x(1)A2+2*x(2)A2+4*x(1)*x(2)+8*x(1)+6*x(2)+1);exp(x(1)*(4*x(2) +4*x(1)+2);(2)编写M文件fun10.m定义约束条件及约束条件的梯度函数： function c,ceq,dc,dceq=fun10(x);c=x(1)*x(2)-x(1)-x(2)+1.5;-x(1)*x(2)-10;dc=x(2)-1,-x(2);x(1)-1,-x(1);ceq=;dceq=;( 3)调用函数 fmincon%采用标准算法 options=optimset(largescale,off);%采用梯度optio

49、ns=optimset(options,GradObj,on,GradConstr,on);x,y=fmincon(fun9,rand(2,1),fun10,options)4 飞行管理问题在约 10，000m 高空的某边长 160km 的正方形区域内，经常有若干架飞机作水平 飞行。区域内每架飞机的位置和速度向量均由计算机记录其数据，以便进行飞行管理。 当一架欲进入该区域的飞机到达区域边缘时，记录其数据后，要立即计算并判断是否会 与区域内的飞机发生碰撞。如果会碰撞，则应计算如何调整各架（包括新进入的）飞机 飞行的方向角，以避免碰撞。现假定条件如下：1）不碰撞的标准为任意两架飞机的距离大于8km;2）飞机飞行方向角调整的幅度不应超过30 度；3）所有飞机飞行速度均为每小时 800km;4）进入该区域的飞机在到达区域边缘时，与区域内飞机的距离应在60km以上；5）最多需考虑6架飞机；6）不必考虑飞机离开此区域后的状况。 请你对这个避免碰撞的飞行管理问题建立数学模型，列出计算步骤，对以下数据进行计算（