第三章 非线性规划 §1 非线性规划 非线性规划的实例与定义 如果目标函数或约束条件中包含非线性函数，就称这种规划问题为非线性规划问题。一般说来，解非线性规划要比解线性规划问题困难得多。而且，也不象线性规划有单纯形法这一通用方法，非线性规划目前还没有适于各种问题的一般算法，各个方法都有自己特定的适用范围。 下面通过实例归纳出非线性规划数学模型的一般形式，介绍有关非线性规划的基本概念。 例1 （投资决策问题）某企业有n个项目可供选择投资，并且至少要对其中一个项目投资。已知该企业拥有总资金A元，投资于第i(i=1,L,n)个项目需花资金a元，i并预计可收益b元。试选择最佳投资方案。 i解 设投资决策变量为 ⎧1,决定投资第i个项目 x=，i=1,L,n， ⎨i0，决定不投资第i个项目⎩nn则投资总额为ax，投资总收益为bx。因为该公司至少要对一个项目投资，并∑ii∑iii=1i=1且总的投资金额不能超过总资金A，故有限制条件 n 0<ax≤A ∑iii=1另外，由于x(i=1,L,n)只取值0或1，所以还有 i x(1−x)=0,i=1,L,n. ii最佳投资方案应是投资额最小而总收益最大的方案，所以这个最佳投资决策问题归结为总资金以及决策变量（取0或1）的限制条件下，极大化总收益和总投资之比。因此，其数学模型为： nbx∑iii=1maxQ= nax∑iii=. 0<ax≤A ∑iii=1 x(1−x)=0,i=1,L,n. ii上面例题是在一组等式或不等式的约束下，求一个函数的最大值（或最小值）问题，其中至少有一个非线性函数，这类问题称之为非线性规划问题。可概括为一般形式 minf(x) (x)≤0,j=1,L,q (NP) j g(x)=0,i=1,L,p i -32-
T其中x=[xLx]称为模型（NP）的决策变量，f称为目标函数，g(i=1,L,p)1ni和h(j=1,L,q)称为约束函数。另外，g(x)=0 (i=1,L,p)称为等式约束，jih(x)≤0 (j=1,L,q)称为不等式的约束。 j对于一个实际问题，在把它归结成非线性规划问题时，一般要注意如下几点： （i）确定供选方案：首先要收集同问题有关的资料和数据，在全面熟悉问题的基础上，确认什么是问题的可供选择的方案，并用一组变量来表示它们。 （ii）提出追求目标：经过资料分析，根据实际需要和可能，提出要追求极小化或极大化的目标。并且，运用各种科学和技术原理，把它表示成数学关系式。 （iii）给出价值标准：在提出要追求的目标之后，要确立所考虑目标的“好”或“坏”的价值标准，并用某种数量形式来描述它。 （iv）寻求限制条件：由于所追求的目标一般都要在一定的条件下取得极小化或极大化效果，因此还需要寻找出问题的所有限制条件，这些条件通常用变量之间的一些不等式或等式来表示。 线性规划与非线性规划的区别 如果线性规划的最优解存在，其最优解只能在其可行域的边界上达到（特别是可行域的顶点上达到）；而非线性规划的最优解（如果最优解存在）则可能在其可行域的任意一点达到。 非线性规划的Matlab解法 Matlab中非线性规划的数学模型写成以下形式 minf(x) Ax≤B⎧⎪Aeq⋅x=Beq⎪ , ⎨C(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定义了线性约束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 求下列非线性规划 222min f(x)=x+x+x+8 . x−x+x≥0 12323 x+x+x≤20 1232 −x−x+2=0 12 -33-
2 x+2x=3 23 x,x,x≥0 123解 （i）编写M文件定义目标函数 function f=fun1(x); f=sum(x.^2)+8; （ii）编写M文件定义非线性约束条件 function [g,h]=fun2(x); g=[-x(1)^2+x(2)-x(3)^2 x(1)+x(2)^2+x(3)^3-20]; %非线性不等式约束 h=[-x(1)-x(2)^2+2 x(2)+2*x(3)^2-3]; %非线性等式约束 （iii）编写主程序文件如下： options=optimset('largescale','off'); [x,y]=fmincon('fun1',rand(3,1),[],[],[],[],zeros(3,1),[], ... 'fun2', options) 就可以求得当x=,x=,x=时，最小值y=。 求解非线性规划的基本迭代格式 记（NP）的可行域为K。 *若x∈K，并且 * f(x)≤f(x),∀x∈K **则称x是（NP）的整体最优解，f(x)是(NP)的整体最优值。如果有 **f(x)<f(x),∀x∈K,x≠x **则称x是（NP）的严格整体最优解，f(x)是(NP)的严格整体最优值。 ***若x∈K，并且存在x的邻域N(x)，使 δ** f(x)≤f(x),∀x∈N(x)IK， δ**则称x是（NP）的局部最优解，f(x)是(NP)的局部最优值。如果有 **f(x)<f(x),∀x∈N(x)IK δ**则称x是（NP）的严格局部最优解，f(x)是(NP)的严格局部最优值。 由于线性规划的目标函数为线性函数，可行域为凸集，因而求出的最优解就是整个可行域上的全局最优解。非线性规划却不然，有时求出的某个解虽是一部分可行域上的极值点，但却并不一定是整个可行域上的全局最优解。 对于非线性规划模型(NP)，可以采用迭代方法求它的最优解。迭代方法的基本思0n想是：从一个选定的初始点x∈R出发，按照某一特定的迭代规则产生一个点列kkk{x}，使得当{x}是有穷点列时，其最后一个点是(NP)的最优解；当{x}是无穷点列时，它有极限点，并且其极限点是(NP)的最优解。 knk+1n设x∈R是某迭代方法的第k轮迭代点，x∈R是第k+1轮迭代点，记 k+1kk x=x+tp （1） k1knkkkk+1这里t∈R,p∈R,p=1，并且p的方向是从点x向着点x的方向。式（1）k就是求解非线性规划模型(NP)的基本迭代格式。 -34-
kk通常，我们把基本迭代格式（1）中的p称为第k轮搜索方向，t为沿p方向的k步长，使用迭代方法求解(NP)的关键在于，如何构造每一轮的搜索方向和确定适当的步长。 n设x∈R,p≠0，若存在δ>0，使 f(x+tp)<f(x),∀t∈(0,δ)， 称向量p是f在点x处的下降方向。 n设x∈R,p≠0，若存在t>0，使 x+tp∈K， 称向量p是点x处关于K的可行方向。 一个向量p，若既是函数f在点x处的下降方向，又是该点关于区域K的可行方向，称之为函数f在点x处关于K的可行下降方向。 现在，我们给出用基本迭代格式（1）求解(NP)的一般步骤如下： 00° 选取初始点x，令k:=0。 k1° 构造搜索方向，依照一定规划，构造f在点x处关于K的可行下降方向作为k搜索方向p。 kk2° 寻求搜索步长。以x为起点沿搜索方向p寻求适当的步长t，使目标函数值k有某种意义的下降。 3° 求出下一个迭代点。按迭代格式（1）求出 k+1kk x=x+tp。 kk+1若x已满足某种终止条件，停止迭代。 k+1k 4° 以x代替x，回到1°步。 凸函数、凸规划 (n)设f(x)为定义在n维欧氏空间E中某个凸集R上的函数，若对任何实数(1)(2)α(0<α<1)以及R中的任意两点x和x，恒有 (1)(2)(1)(2) f(αx+(1−α)x)≤αf(x)+(1−α)f(x) 则称f(x)为定义在R上的凸函数。 (1)(2)若对每一个α(0<α<1)和x≠x∈R恒有 (1)(2)(1)(2) f(αx+(1−α)x)<αf(x)+(1−α)f(x) 则称f(x)为定义在R上的严格凸函数。 考虑非线性规划 min f(x)⎧⎪x∈R ⎨R={x|g(x)≤0,j=1,2,L,l}⎪j⎩假定其中f(x)为凸函数，g(x)(j=1,2,L,l)为凸函数，这样的非线性规划称为j凸规划。 可以证明，凸规划的可行域为凸集，其局部最优解即为全局最优解，而且其最优解的集合形成一个凸集。当凸规划的目标函数f(x)为严格凸函数时，其最优解必定唯一（假定最优解存在）。由此可见，凸规划是一类比较简单而又具有重要理论意义的非 -35-
线性规划。 §2 无约束问题 一维搜索方法 当用迭代法求函数的极小点时，常常用到一维搜索，即沿某一已知方向求目标函数的极小点。一维搜索的方法很多，常用的有：（1）试探法（“成功—失败”，斐波那契法，法等）；（2）插值法（抛物线插值法，三次插值法等）；（3）微积分中的求根法（切线法，二分法等）。 考虑一维极小化问题 minf(t) （2） a≤t≤b若f(t)是[a,b]区间上的下单峰函数，我们介绍通过不断地缩短[a,b]的长度，来搜索得（2）的近似最优解的两个方法。 *为了缩短区间[a,b]，逐步搜索得（2）的最优解t的近似值，我们可以采用以下途径：在[a,b]中任取两个关于[a,b]是对称的点t和t（不妨设t<t，并把它们叫1221做搜索点），计算f(t)和f(t)并比较它们的大小。对于单峰函数，若f(t)<f(t)，1221*则必有t∈[a,t]，因而[a,t]是缩短了的单峰区间；若f(t)<f(t)，则有1112*t∈[t,b]，故[t,b]是缩短了的单峰区间；若f(t)=f(t)，则[a,t]和[t,b]都是222112缩短了的单峰。因此通过两个搜索点处目标函数值大小的比较，总可以获得缩短了的单峰区间。对于新的单峰区间重复上述做法，显然又可获得更短的单峰区间。如此进行，在单峰区间缩短到充分小时，我们可以取最后的搜索点作为（2）最优解的近似值。 应该按照怎样的规则来选取探索点，使给定的单峰区间的长度能尽快地缩短？ Fibonacci法 如用F表示计算n个函数值能缩短为单位长区间的最大原区间长度，可推出F满nn足关系 F=F=1 01 F=F+F,n=2,3,L, nn−2n−1数列{F}称为Fibonacci数列，F称为第n个Fibonacci数，相邻两个Fibonacci数之nnFn−1比称为Fibonacci分数。 Fn 当用斐波那契法以n个探索点来缩短某一区间时，区间长度的第一次缩短率为FFFFn−1n−2n−31，其后各次分别为,,L,。由此，若t和t(t<t)是单峰区间[a,b]1221FFFFnn−1n−22中第1个和第2个探索点的话，那么应有比例关系 t−aFt−aF1n−12n−2=， = b−aFb−aFnn从而 FFn−1n−2 t=a+(b−a)，t=a+(b−a) （3） 12FFnn它们关于[a,b]确是对称的点。 -36-
如果要求经过一系列探索点搜索之后，使最后的探索点和最优解之间的距离不超过精度δ>0，这就要求最后区间的长度不超过δ，即 b−a ≤δ （4） Fn据此，我们应按照预先给定的精度δ，确定使（4）成立的最小整数n作为搜索次数，直到进行到第n个探索点时停止。 用上述不断缩短函数f(t)的单峰区间[a,b]的办法，来求得问题（2）的近似解，是Kiefer(1953年)提出的，叫做Finbonacci法，具体步骤如下： o1 选取初始数据，确定单峰区间[a,b]，给出搜索精度δ>0，由（4）确定搜00索次数n。 o2 k=1,a=a,b=b，计算最初两个搜索点，按（3）计算t和t。 0012o3 while k<n−1 f=f(t),f=f(t) 1122 if f<f 12F(n−1−k) a=t;t=t;t=a+(b−a) 2211F(n−k) else F(n−1−k) b=t;t=t;t=b+(a−b) 1122F(n−k)end k=k+1 end o4 当进行至k=n−1时， 1 t=t=(a+b) 122这就无法借比较函数值f(t)和f(t)的大小确定最终区间，为此，取 121⎧t=(a+b)2⎪⎪2 ⎨1⎪t=a+(+ε)(b−a)1⎪⎩2其中ε为任意小的数。在t和t这两点中，以函数值较小者为近似极小点，相应的函数12值为近似极小值。并得最终区间[a,t]或[t,b]。 12 由上述分析可知，斐波那契法使用对称搜索的方法，逐步缩短所考察的区间，它能以尽量少的函数求值次数，达到预定的某一缩短率。 2例3 试用斐波那契法求函数f(t)=t−t+2的近似极小点，要求缩短后的区间不大于区间[−1,3]的倍。 程序留作习题。 法 若ω>0，满足比例关系 -37-
ω1−ω = 1ω5−1称之为黄金分割数，其值为ω==。 2黄金分割数ω和Fibonacci分数之间有着重要的关系 Fn−1ω=lim。 n→∞Fn现用不变的区间缩短率，代替斐波那契法每次不同的缩短率，就得到了黄金分割法（法）。这个方法可以看成是斐波那契法的近似，实现起来比较容易，效果也相当好，因而易于为人们所接受。 用法求解，从第2个探索点开始每增加一个探索点作一轮迭代以后，原单峰区间要缩短倍。计算n个探索点的函数值可以把原区间[a,b]连续缩短n−100次，因为每次的缩短率均为μ，故最后的区间长度为 n−1 (b−a)μ 00这就是说，当已知缩短的相对精度为δ时，可用下式计算探索点个数n： n−1μ≤δ 当然，也可以不预先计算探索点的数目n，而在计算过程中逐次加以判断，看是否已满足了提出的精度要求。 法是一种等速对称进行试探的方法，每次的探索点均取在区间长度的倍和倍处。 二次插值法 对极小化问题（2），当f(t)在[a,b]上连续时，可以考虑用多项式插值来进行一维搜索。它的基本思想是：在搜索区间中，不断用低次（通常不超过三次）多项式来近似目标函数，并逐步用插值多项式的极小点来逼近（2）的最优解。 无约束极值问题的解法 无约束极值问题可表述为 (n)min f(x),x∈E （5） 求解问题（5）的迭代法大体上分为两点：一是用到函数的一阶导数或二阶导数，称为解析法。另一是仅用到函数值，称为直接法。 解析法 梯度法（最速下降法） 对基本迭代格式 k+1kkx=x+tp （6） kkk我们总是考虑从点x出发沿哪一个方向p，使目标函数f下降得最快。微积分的知识k告诉我们，点x的负梯度方向 kk p=−∇f(x)， kkk是从点x出发使f下降最快的方向。为此，称负梯度方向−∇f(x)为f在点x处的最速下降方向。 -38-
kk按基本迭代格式（6），每一轮从点x出发沿最速下降方向−∇f(x)作一维搜索，来建立求解无约束极值问题的方法，称之为最速下降法。 这个方法的特点是，每轮的搜索方向都是目标函数在当前点下降最快的方向。同kk时，用∇f(x)=0或∇f(x)≤ε作为停止条件。其具体步骤如下： 01°选取初始数据。选取初始点x，给定终止误差，令k:=0。 kkk2°求梯度向量。计算∇f(x)， 若∇f(x)≤ε，停止迭代，输出x。否则，进行3°。 3° 构造负梯度方向。取 kkp=−∇f(x). 4° 进行一维搜索。求t，使得 kkkkk f(x+tp)=minf(x+tp) kt≥0k+1kk令x=x+tp,k:=k+1,转2°。 k例4 用最速下降法求解无约束非线性规划问题 22 minf(x)=x+25x 12T0T其中x=(x,x)，要求选取初始点x=(2,2)。 12T解 （i）∇f(x)=(2x,50x) 12编写M文件，定义函数f(x)及其梯度列向量如下 function [f,df]=detaf(x); f=x(1)^2+25*x(2)^2; df=[2*x(1) 50*x(2)]; （ii）编写主程序文件如下： clc x=[2;2]; [f0,g]=detaf(x); while norm(g)> p=-g/norm(g); t=;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 x,f0 Newton法 k考虑目标函数f在点x处的二次逼近式 1kkTkkT2kkf(x)≈Q(x)=f(x)+∇f(x)(x−x) +(x−x)∇f(x)(x−x) 2假定Hesse阵 -39-
2k2k⎡∂f(x)∂⎤f(x)L⎢2⎥∂x∂x∂x11n⎢⎥2k∇f(x)=MLM ⎢⎥2k2k∂f(x)∂f(x)⎢⎥L2⎢⎥∂x∂x∂xn1n⎣⎦正定。 2kk+1由于∇f(x)正定，函数Q的驻点x是Q(x)的极小点。为求此极小点，令 k+1k2kk+1k∇Q(x)=∇f(x)+∇f(x)(x−x)=0， 即可解得 k+1k2k−1kx=x−[∇f(x)]∇f(x). k对照基本迭代格式（1），可知从点x出发沿搜索方向。 k2k−1k p=−[∇f(x)]∇f(x) k+1kk并取步长t=1即可得Q(x)的最小点x。通常，把方向p叫做从点x出发的kNewton方向。从一初始点开始，每一轮从当前迭代点出发，沿Newton方向并取步长为1的求解方法，称之为Newton法。其具体步骤如下： 01°选取初始数据。选取初始点x，给定终止误差ε>0，令k:=0。 kkk2°求梯度向量。计算∇f(x)，若∇f(x)≤ε，停止迭代，输出x。否则，进行3°。 2k−13°构造Newton方向。计算[∇f(x)]，取 k2k−1k p=−[∇f(x)]∇f(x). k+1kk4° 求下一迭代点。令x=x+p,k:=k+1，转2°。 例5 用Newton法求解， 4422 minf(x)=x+25x+xx 12120T选取x=(2,2)。 3232T解：（i）∇f(x)=[4x+2xx100x+2xx] 11221222⎡⎤12x+2x4xx21212∇f= ⎢⎥224xx300x+2x⎣1221⎦（ii）编写M文件如下： function [f,df,d2f]=nwfun(x); f=x(1)^4+25*x(2)^4+x(1)^2*x(2)^2; df=[4*x(1)^3+2*x(1)*x(2)^2;100*x(2)^3+2*x(1)^2*x(2)]; d2f=[2*x(1)^2+2*x(2)^2,4*x(1)*x(2) 4*x(1)*x(2),300*x(2)^2+2*x(1)^2]; （III）编写主程序文件如下： clc x=[2;2]; [f0,g1,g2]=nwfun(x); while norm(g1)> p=-inv(g2)*g1; x=x+p; -40-
[f0,g1,g2]=nwfun(x); end x, f0 如果目标函数是非二次函数，一般地说，用Newton法通过有限轮迭代并不能保证可求得其最优解。 为了提高计算精度，我们在迭代时可以采用变步长计算上述问题，编写主程序文件example5_2如下： clc,clear x=[2;2]; [f0,g1,g2]=nwfun(x); while norm(g1)> p=-inv(g2)*g1;p=p/norm(p); t=;f=nwfun(x+t*p); while f>f0 t=t/2;f=nwfun(x+t*p); end x=x+t*p; [f0,g1,g2]=nwfun(x); end x,f0 Newton法的优点是收敛速度快；缺点是有时不好用而需采取改进措施，此外，当2k−1维数较高时，计算−[∇f(x)]的工作量很大。 变尺度法 变尺度法（Variable Metric Algorithm）是近20多年来发展起来的，它不仅是求解无约束极值问题非常有效的算法，而且也已被推广用来求解约束极值问题。由于它既避免了计算二阶导数矩阵及其求逆过程，又比梯度法的收敛速度快，特别是对高维问题具有显著的优越性，因而使变尺度法获得了很高的声誉。下面我们就来简要地介绍一种变尺度法—DFP法的基本原理及其计算过程。这一方法首先由Davidon在1959年提出，后经Fletcher和Powell加以改进。 2k−1k 我们已经知道，牛顿法的搜索方向是−[∇f(x)]∇f(x)，为了不计算二阶2k导数矩阵[∇f(x)]及其逆阵，我们设法构造另一个矩阵，用它来逼近二阶导数矩阵2k−1的逆阵[∇f(x)]，这一类方法也称拟牛顿法（Quasi-Newton Method）。 (k) 下面研究如何构造这样的近似矩阵，并将它记为H。我们要求：每一步都能以现有的信息来确定下一个搜索方向；每做一次选代，目标函数值均有所下降；这些近似矩阵最后应收敛于解点处的Hesse阵的逆阵。 kk+1当f(x)是二次函数时，其Hesse阵为常数阵A，任两点x和x处的梯度之差为 k+1kk+1k ∇f(x)−∇f(x)=A(x−x) 或 k+1k−1k+1kx−x=A[∇f(x)−∇f(x)] 对于非二次函数，仿照二次函数的情形，要求其Hesse阵的逆阵的第k+1次近似(k+1)矩阵H满足关系式 k+1k(k+1)k+1kx−x=H[∇f(x)−∇f(x)] （7） -41-
这就是常说的拟Newton条件。 若令 (k)k+1k⎧ΔG=∇f(x)−∇f(x) （8） ⎨kk+1kΔx=x−x⎩则式（7）变为 k(k+1)(k) Δx=HΔG， （9） (k)(k+1)(k)(k+1)现假定H已知，用下式求H（设H和H均为对称正定阵）； (k+1)(k)(k)H=H+ΔH （10） (k)(k+1)其中ΔH称为第k次校正矩阵。显然，H应满足拟Newton条件（9），即要求 k(k)(k)(k)Δx=(H+ΔH)ΔG 或 (k)(k)k(k)(k)ΔHΔG=Δx−HΔG （11） (k)由此可以设想， ΔH的一种比较简单的形式是 (k)k(k)T(k)(k)(k)TΔH=Δx(Q)−HΔG(W) (12) (k)(k)其中Q和W为两个待定列向量。 (k)将式（12）中的ΔH代入（11），得 k(k)T(k)(k)(k)(k)T(k)k(k)(k)Δx(Q)ΔG−HΔG(W)ΔG=Δx−HΔG 这说明，应使 (k)T(k)(k)T(k)(Q)ΔG=(W)ΔG=1 （13） (k)考虑到ΔH应为对称阵，最简单的办法就是取 (k)k⎧Q=ηΔx⎪k （14） ⎨(k)(k)(k)⎪W=ξHΔG⎩k由式（13）得 kT(k)(k)T(k)(k)η(Δx)ΔG=ξ(ΔG)HΔG=1 （15） kkkT(k)(k)T(k)(k)若(Δx)ΔG和(ΔG)HΔG不等于零，则有 11⎧η==kkT(k)(k)Tk⎪(Δx)ΔG(ΔG)Δx⎪ （16） ⎨1⎪ξ=k(k)T(k)(k)⎪(ΔG)HΔG⎩于是，得校正矩阵 kkT(k)(k)(k)T(k)Δx(Δx)HΔG(G)ΔH(k)ΔH=− (17) (k)Tk(k)T(k)(k)(ΔG)Δx(ΔG)HΔG从而得到 kkT(k)(k)(k)T(k)Δx(Δx)HΔG(G)ΔH(k+1)(k)H=H+− （18） (k)Tk(k)T(k)(k)(ΔG)Δx(ΔG)HΔG(0)上述矩阵称为尺度矩阵。通常，我们取第一个尺度矩阵H为单位阵，以后的尺度矩 -42-
阵按式（18）逐步形成。可以证明： k(k)（i）当x不是极小点且H正定时，式（17）右端两项的分母不为零，从而可(k+1)按式（18）产生下一个尺度矩阵H； (k)(k+1)（ii）若H为对称正定阵，则由式（18）产生的H也是对称正定阵； （iii）由此推出DFP法的搜索方向为下降方向。 现将DFP变尺度法的计算步骤总结如下。 01°给定初始点x及梯度允许误差ε>0。 002°若∇f(x)≤ε，则x即为近似极小点，停止迭代，否则，转向下一步。 3°令 (0)H=I（单位矩阵）， 0(0)0p=−H∇f(x) 0在p方向进行一维搜索，确定最佳步长λ： 00000 minf(x+λp)=f(x+λp) 0λ如此可得下一个近似点 100 x=x+λp 0kk 4°一般地，设已得到近似点x，算出∇f(x)，若 0 ∇f(x)≤ε k(k)则x即为所求的近似解，停止迭代；否则，计算H： k−1k−1T(k−1)(k−1)(k−1)T(k−1)Δx(Δx)HΔG(G)ΔH(k)(k−1)H=H+ − (k−1)Tk−1(k−1)T(k−1)(k−1)(ΔG)Δx(ΔG)HΔGk(k)kk并令p=−H∇f(x)，在p方向上进行一维搜索，得λ，从而可得下一个近似点 kk+1kk x=x+λp kk+1k+15°若x满足精度要求，则x即为所求的近似解，否则，转回4°，直到求出某点满足精度要求为止。 直接法 在无约束非线性规划方法中，遇到问题的目标函数不可导或导函数的解析式难以表示时，人们一般需要使用直接搜索方法。同时，由于这些方法一般都比较直观和易于理解，因而在实际应用中常为人们所采用。下面我们介绍Powell方法。 这个方法主要由所谓基本搜索、加速搜索和调整搜索方向三部分组成，具体步骤如下： 01° 选取初始数据。选取初始点x，n个线性无关初始方向，组成初搜索方向组01n−1{p,p,L,p}。给定终止误差ε>0，令k:=0。 0k01n−12°进行基本搜索。令y:=x，依次沿{p,p,L,p}中的方向进行一维搜12n索。对应地得到辅助迭代点y,y,L,y，即 jj−1j−1 y=y+tp j−1j−1j−1j−1j−1f(y+tp)=minf(y+tp),j=1,L,n j−1t≥0 -43-
nn0nk+1n3°构造加速方向。令p=y−y，若p≤ε，停止迭代，输出x=y。否则进行4°。 4°确定调整方向。按下式 m−1mj−1jf(y)−f(y)=max{f(y)−f(y)|1≤j≤n} 找出m。若 0nn0m−1mf(y)−2f(y)+f(2y−y)<2[f(y)−f(y)] 成立，进行5°。否则，进行6°。 5°调整搜索方向组。令 k+1nnnnnnx=y+tp:f(y+tp)=minf(y+tp). nnt≥0同时，令 01n+10m−1m+1n−1n {p,p,L,p}:={p,L,p,p,L,p,p}， k+1k:=k+1，转2°。 k+1n 6°不调整搜索方向组。令x:=y,k:=k+1，转2°。 Matlab求无约束极值问题 在Matlab工具箱中，用于求解无约束极值问题的函数有fminunc和fminsearch，用法介绍如下。 求函数的极小值 minf(x) , x其中x可以为标量或向量。 Matlab中fminunc的基本命令是 [X,FVAL]=FMINUNC(FUN,X0,OPTIONS,P1,P2, ...) 其中的返回值X是所求得的极小点，FVAL是函数的极小值，其它返回值的含义参见相关的帮助。FUN 是一个M文件，当FUN只有一个返回值时，它的返回值是函数f(x)；当FUN有两个返回值时，它的第二个返回值是f(x)的梯度向量；当FUN有三个返回值时，它的第三个返回值是f(x)的二阶导数阵（Hessian阵）。X0是向量x的初始值，OPTIONS是优化参数，可以使用缺省参数。P1，P2是可以传递给FUN的一些参数。 222例6 求函数f(x)=100(x−x)+(1−x)的最小值。 211解：编写M文件如下： function [f,g]=fun2(x); f=100*(x(2)-x(1)^2)^2+(1-x(1))^2; g=[-400*x(1)*(x(2)-x(1)^2)-2*(1-x(1));200*(x(2)-x(1)^2)]; 编写主函数文件如下： options = optimset('GradObj','on'); [x,y]=fminunc('fun2',rand(1,2),options) 即可求得函数的极小值。 在求极值时，也可以利用二阶导数，编写M文件如下： function [f,df,d2f]=fun3(x); f=100*(x(2)-x(1)^2)^2+(1-x(1))^2; df=[-400*x(1)*(x(2)-x(1)^2)-2*(1-x(1));200*(x(2)-x(1)^2)]; d2f=[-400*x(2)+1200*x(1)^2+2,-400*x(1) -400*x(1),200]; -44-
编写主函数文件如下： options = optimset('GradObj','on','Hessian','on'); [x,y]=fminunc('fun3',rand(1,2),options) 即可求得函数的极小值。 求多元函数的极值也可以使用Matlab的fminsearch命令，其使用格式为： [X,FVAL,EXITFLAG,OUTPUT]＝FMINSEARCH(FUN,X0,OPTIONS,P1,P2,...) 例7 求函数f(x)=sin(x)+3取最小值时的x值。 解 编写f(x)的M文件如下： function f=fun4(x); f=sin(x)+3; 编写主函数文件如下： x0=2; [x,y]=fminsearch(@fun4,x0) 即求得在初值2附近的极小点及极小值。 §3 约束极值问题 带有约束条件的极值问题称为约束极值问题，也叫规划问题。 求解约束极值问题要比求解无约束极值问题困难得多。为了简化其优化工作，可采用以下方法：将约束问题化为无约束问题；将非线性规划问题化为线性规划问题，以及能将复杂问题变换为较简单问题的其它方法。 库恩—塔克条件是非线性规划领域中最重要的理论成果之一，是确定某点为最优点的必要条件，但一般说它并不是充分条件（对于凸规划，它既是最优点存在的必要条件，同时也是充分条件）。 二次规划 若某非线性规划的目标函数为自变量x的二次函数，约束条件又全是线性的，就称这种规划为二次规划。 Matlab中二次规划的数学模型可表述如下： 1TT min xHx+fx, 2Ax≤b⎧. ⎨Aeq⋅x=beq⎩这里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 求解二次规划 22⎧ min f(x)=2x−4xx+4x−6x−3x 112212⎪ x+x≤3⎪12 ⎨ 4x+x≤9⎪12⎪ x,x≥0⎩12解 编写如下程序： -45-
h=[4,-4;-4,8]; f=[-6;-3]; a=[1,1;4,1]; b=[3;9]; [x,value]=quadprog(h,f,a,b,[],[],zeros(2,1)) 求得 ⎡⎤ x= , Min f(x)=−。 ⎢⎥⎣⎦ 罚函数法 利用罚函数法，可将非线性规划问题的求解，转化为求解一系列无约束极值问题，因而也称这种方法为序列无约束最小化技术，简记为 SUMT (Sequential Unconstrained Minization Technique)。 罚函数法求解非线性规划问题的思想是，利用问题中的约束函数作出适当的罚函数，由此构造出带参数的增广目标函数，把问题转化为无约束非线性规划问题。主要有两种形式，一种叫外罚函数法，另一种叫内罚函数法，下面介绍外罚函数法。 考虑问题： minf(x) ⎧g(x)≤0, i=1,L,r,i⎪. h(x)≥0,j=1,L,s, ⎨j⎪k(x)=0,m=1,L,t ⎩m取一个充分大的数 M>0 ，构造函数 rstP(x,M)=f(x)+Mmax(g(x),0)−Mmin(h(x),0)+M|k(x)| ∑∑∑iiii=1i=1i=1⎛G(x)⎞⎛H(x)⎞⎛⎞⎛⎞⎜⎟⎜⎟（或P(x,M)=f(x)+Msummax⎜⎟−Msummin⎜⎟+M||K(x)|| ⎜⎟⎜⎟⎜⎟⎜⎟00⎝⎠⎝⎠⎝⎠⎝⎠这里G(x)=[g(x)Lg(x)]，H(x)=[h(x)Lh(x)]， 1r1sK(x)=k(x)Lk(t)，Matlab中可以直接利用 max、min和sum函数。）则以1t增广目标函数P(x,M)为目标函数的无约束极值问题 minP(x,M) 的最优解x也是原问题的最优解。 例9 求下列非线性规划 22⎧min f(x)=x+x+812⎪2 x−x≥0⎪12 ⎨2 −x−x+2=0⎪12⎪ x,x≥0.⎩12解 (i)编写 M 文件 function 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)+... -46-
M*abs(-x(1)-x(2)^2+2); 或者是利用Matlab的求矩阵的极小值和极大值函数编写如下： 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); 我们也可以修改罚函数的定义，编写如下： function g=test(x); M=50000; f=x(1)^2+x(2)^2+8; g=f-M*min(min(x),0)-M*min(x(1)^2-x(2),0)+M*(-x(1)-x(2)^2+2)^2; (ii)在Matlab命令窗口输入 [x,y]=fminunc('test',rand(2,1)) 即可求得问题的解。 Matlab求约束极值问题 在Matlab优化工具箱中，用于求解约束最优化问题的函数有：fminbnd、fmincon、quadprog、fseminf、fminimax，上面我们已经介绍了函数fmincon和quadprog。 fminbnd函数 求单变量非线性函数在区间上的极小值 minf(x) , x∈[x,x] 12xMatlab的命令为 [X,FVAL] = FMINBND(FUN,x1,x2,OPTIONS)， 它的返回值是极小点x和函数的极小值。这里fun 是用M文件定义的函数或Matlab中的单变量数学函数。 2例10 求函数 f(x)=(x−3)−1, x∈[0,5] 的最小值。 解 编写M文件 function f=fun5(x); f=(x-3)^2-1; 在Matlab的命令窗口输入 [x,y]=fminbnd('fun5',0,5) 即可求得极小点和极小值。 fseminf函数 求 min{F(x)|C(x)≤0,Ceq(x)=0,PHI(x,w)≤0} xA*x≤B⎧. ⎨Aeq*x=Beq⎩其中C(x),Ceq(x),PHI(x,w)都是向量函数；w是附加的向量变量，w的每个分量都限定在某个区间内。 上述问题的Matlab命令格式为 X=FSEMINF(FUN,X0,NTHETA,SEMINFCON,A,B,Aeq,Beq) 其中FUN用于定义目标函数F(x)；X0为x的初始值；NTHETA是半无穷约束PHI(x,w)的个数；函数SEMINFCON 用于定义非线性不等式约束C(x)，非线性等 -47-
式约束Ceq(x)和半无穷约束PHI(x,w)的每一个分量函数，函数SEMINFCON有两个输入参量X 和 S，S是推荐的取样步长，也许不被使用。 222例11 求函数f(x)=(x−)+(x−)+(x−)取最小值时的x值,123约束为： 12K(x,w)=sin(wx)cos(wx)−(w−50)−sin(wx)−x≤1 1111121133100012K(x,w)=sin(wx)cos(wx)−(w−50)−sin(wx)−x≤1 222221223310001≤w≤100，1≤w≤100 12解 （1）编写M文件定义目标函数如下： function f=fun6(x,s); f=sum(().^2); （2）编写M文件定义约束条件如下： function [c,ceq,k1,k2,s]=fun7(x,s); c=[];ceq=[]; if isnan(s(1,1)) s=[,0; 0]; end %取样值 w1=1:s(1,1):100; w2=1:s(2,1):100; %半无穷约束 k1=sin(w1*x(1)).*cos(w1*x(2))-1/1000*(w1-50).^2-sin(w1*x(3))-x(3)-1; k2=sin(w2*x(2)).*cos(w2*x(1))-1/1000*(w2-50).^2-sin(w2*x(3))-x(3)-1; %画出半无穷约束的图形 plot(w1,k1,'-',w2,k2,'+'); （3）调用函数fseminf 在Matlab的命令窗口输入 [x,y]=fseminf(@fun6,rand(3,1),2,@fun7) 即可。 fminimax函数 ⎧⎫求解minmaxF(x) ⎨⎬xF⎩i⎭A*x≤b⎧⎪Aeq*x=Beq⎪⎪. C(x)≤0 ⎨⎪Ceq(x)=0⎪⎪LB≤x≤UB⎩其中F(x)={F(x),L,F(x)}。 1m上述问题的Matlab命令为 X=FMINIMAX(FUN,X0,A,B,Aeq,Beq,LB,UB,NONLCON) -48-
例12 求函数族{f(x),f(x),f(x),f(x),f(x)}取极大极小值时的x值。其中： 1234522⎧f(x)=2x+x−48x−40x+30411212⎪22f(x)=−x−3x⎪212⎪f(x)=x+3x−18 ⎨312⎪f(x)=−x−x412⎪⎪f(x)=x+x−8512⎩解 （1）编写M文件定义向量函数如下： function f=fun8(x); f=[2*x(1)^2+x(2)^2-48*x(1)-40*x(2)+304 -x(1)^2-3*x(2)^2 x(1)+3*x(2)-18 -x(1)-x(2) x(1)+x(2)-8]; （2）调用函数fminimax [x,y]=fminimax(@fun8,rand(2,1)) 利用梯度求解约束优化问题 x221例13 已知函数f(x)=e(4x+2x+4xx+2x+1)，且满足非线性约束： 12122xx−x−x≤−⎧1212 ⎨xx≥−10⎩12求minf(x) x分析：当使用梯度求解上述问题时，效率更高并且结果更准确。 题目中目标函数的梯度为： x221⎡⎤e(4x+2x+4xx+8x+6x+1)121212 ⎢⎥x1e(4x+4x+2)⎢⎥⎣12⎦解 （1）编写M文件定义目标函数及梯度函数： function [f,df]=fun9(x); f=exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+2*x(2)+1); df=[exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+8*x(1)+6*x(2)+1);exp(x(1))*(4*x(2)+4*x(1)+2)]; （2）编写M文件定义约束条件及约束条件的梯度函数： function [c,ceq,dc,dceq]=fun10(x); c=[x(1)*x(2)-x(1)-x(2)+;-x(1)*x(2)-10]; dc=[x(2)-1,-x(2);x(1)-1,-x(1)]; ceq=[];dceq=[]; （3）调用函数fmincon，编写主函数文件如下： %采用标准算法 options=optimset('largescale','off'); %采用梯度 options=optimset(options,'GradObj','on','GradConstr','on'); [x,y]=fmincon(@fun9,rand(2,1),[],[],[],[],[],[],@fun10,options) -49-
Matlab优化工具箱的用户图形界面解法 Matlab优化工具箱中的optimtool命令提供了优化问题的用户图形界面解法。optimtool可应用到所有优化问题的求解，计算结果可以输出到Matlab工作空间中。 图1 优化问题用户图形界面解法示意图 例14 用optimtool重新求解例1。 利用例1已经定义好的函数fun1和fun2。在Matlab命令窗口运行optimtool，就打开图形界面，如图1所示，填入有关的参数，未填入的参数取值为空或者为默认值，然后用鼠标点一下“start”按钮，就得到求解结果，再使用“file”菜单下的“Export to Workspace…”选项，把计算结果输出到Matlab工作空间中去。 §4 飞行管理问题 在约10，000m高空的某边长160km的正方形区域内，经常有若干架飞机作水平飞行。区域内每架飞机的位置和速度向量均由计算机记录其数据，以便进行飞行管理。当一架欲进入该区域的飞机到达区域边缘时，记录其数据后，要立即计算并判断是否会与区域内的飞机发生碰撞。如果会碰撞，则应计算如何调整各架（包括新进入的）飞机飞行的方向角，以避免碰撞。现假定条件如下： 1）不碰撞的标准为任意两架飞机的距离大于8km; 2）飞机飞行方向角调整的幅度不应超过30度； 3）所有飞机飞行速度均为每小时800km; 4）进入该区域的飞机在到达区域边缘时，与区域内飞机的距离应在60km以上； 5）最多需考虑6架飞机； 6）不必考虑飞机离开此区域后的状况。 请你对这个避免碰撞的飞行管理问题建立数学模型，列出计算步骤，对以下数据进行计算（方向角误差不超过度），要求飞机飞行方向角调整的幅度尽量小。 设该区域4个顶点的座标为(0,0)，(160,0)，(160,160)，(0,160)。记录数据见表1。 -50-
表1 飞行记录数据 飞机编号 横座标x 纵座标y 方向角（度） 1150140243285852363 150 155 4145501595130150230新进入0 0 52 注：方向角指飞行方向与x轴正向的夹角。 为方便以后的讨论，我们引进如下记号： D为飞行管理区域的边长； Ω为飞行管理区域，取直角坐标系使其为[0,D]×[0,D]； a为飞机飞行速度，a=800km/h； 00(x,y)为第i架飞机的初始位置； ii(x(t),y(t))为第i架飞机在t时刻的位置； ii00θ为第i架飞机的原飞行方向角，即飞行方向与x轴夹角，0≤θ<2π； iiππΔθ为第i架飞机的方向角调整，−≤Δθ≤； ii660θ=θ+Δθ为第i架飞机调整后的飞行方向角。 模型一 根据相对运动的观点在考察两架飞机i和j的飞行时，可以将飞机i视为不动而飞机j以相对速度 v=v−v=(acosθ−acosθ,asinθ−asinθ) （19） ijjijiji相对于飞机i运动，对（19）式进行适当的化约可得 θ−θθ+θθ+θjijiji v=2asin(−sin,cos) 222θ−θθ+θθ+θππjijiji =2asin(cos(+),cos(+)) （20） 22222θ+θπij不妨设θ≥θ，此时相对飞行方向角为β=+，见图2。 jiij22 图2 相对飞行方向角 由于两架飞机的初始距离为 002002 r(0)=(x−x)+(y−y) （21） ijijij -51-
80 α=arcsin （22） ijr(0)ij则只要当相对飞行方向角β满足 ij00 α<β<2π−α （23） ijijij时，两架飞机不可能碰撞（见图2）。 0 记β为调整前第j架飞机相对于第i架飞机的相对速度（矢量）与这两架飞机连ij线（从j指向i的矢量）的夹角（以连线矢量为基准，逆时针方向为正，顺时针方向为负）。则由式（23）知，两架飞机不碰撞的条件为 100 β+(Δθ+Δθ)>α （24） ijijij2其中 0 β＝相对速度v的幅角－从n指向m的连线矢量的幅角 mnmniθiθnme−e =arg (x+iy)−(x+iy)mmnn0（注意β表达式中的i表示虚数单位）这里我们利用复数的幅角，可以很方便地计算mn0角度β（m,n=1,2,L,6）。 mn本问题中的优化目标函数可以有不同的形式：如使所有飞机的最大调整量最小；所有飞机的调整量绝对值之和最小等。这里以所有飞机的调整量绝对值之和最小为目标函数，可以得到如下的数学规划模型： 6 minΔθ ∑ii=. β+(Δθ+Δθ)>α，i,j=1,2,L,6，i≠j ijijij2o Δθ≤30， i=1,2,L,6 i利用如下的程序： clc,clear x0=[150 85 150 145 130 0]; y0=[140 85 155 50 150 0]; q=[243 236 159 230 52]; xy0=[x0; y0]; d0=dist(xy0); %求矩阵各个列向量之间的距离 d0(find(d0==0))=inf; a0=asind(8./d0) %以度为单位的反函数 xy1=x0+i*y0 xy2=exp(i*q*pi/180) for m=1:6 for n=1:6 if n~=m b0(m,n)=angle((xy2(n)-xy2(m))/(xy1(m)-xy1(n))); end -52-
end end b0=b0*180/pi; dlmwrite('',a0,'delimiter', '\t','newline','PC'); fid=fopen('','a'); fwrite(fid,'~','char'); %往纯文本文件中写LINGO数据的分割符 dlmwrite('',b0,'delimiter', '\t','newline','PC','-append','roffset', 1) 0求得α的值如表2所示。 ij0表2 α的值 ij 1 2 3 4 5 6 0 1 0 2 3 4 5 0 6 0 求得β的值如表3所示。 ij0表3 β的值 ij 1 2 3 4 5 6 0 1 0 9 2 0 3 0 4 0 5 9 0 6 上述飞行管理的数学规划模型可如下输入LINGO求解： model: sets: plane/1..6/:delta; link(plane,plane):alpha,beta; endsets data: alpha=@file(''); !需要在alpha的数据后面加上分隔符"~"; beta=@file(''); enddata min=@sum(plane:@abs(delta)); @for(plane:@bnd(-30,delta,30)); @for(link(i,j)|i#ne#j:@abs(beta(i,j)+*delta(i)+*delta(j))>alpha(i,j)); end ooo求得的最优解为Δθ=，Δθ=−，Δθ=，其它调整356角度为0。 模型二 -53-
两架飞机i,j不发生碰撞的条件为 22(x(t)−x(t))+(y(t)−y(t))>64 （25） ijij1≤i≤5，i+1≤j≤6，0≤t≤min{T,T} ij其中T,T分别表示第i,j架飞机飞出正方形区域边界的时刻。这里 ij00x(t)=x+atcosθ，y(t)=y+atsinθ，i=1,2,L,n； iiiiiiπ0θ=θ+Δθ，|Δθ|≤，i=1,2,L,n； iiii6 下面我们把约束条件（25）加强为对所有的时间t都成立，记 ~222~~l=(x(t)−x(t))+(y(t)−y(t))−64=a(i,j)t+b(i,j)t+c(i,j) i,jijijθ−θij22~其中a(i,j)=4asin， 2~b(i,j)=2a[(x(0)−x(0))(cosθ−cosθ)+(y(0)−y(0))(sinθ−sinθ)] ijijijij22~c(i,j)=(x(0)−x(0))+(y(0)−y(0))−64 ijij则两架i,j飞机不碰撞的条件是 ~2~~Δ(i,j)=b(i,j)−4a(i,j)c(i,j)<0 （26） 这样我们建立如下的非线性规划模型 62 (Δθ) ∑ii=1 . Δ(i,j)<0，1≤i≤5， i+1≤j≤6 π Δθ≤，i=1,2,L,6 i6习 题 三 1. 用最速下降法（梯度法）求函数： 22 f(x)=4x+6x−2x−2xx−2x 1211220T的极大点。给定初始点x=(1,1)。 2. 试用牛顿法求解： 1minf(x)=− 22x+x+212(0)T取初始点x=(4,0)，并将采用变步长和采用固定步长λ=时的情形做比较。 3. 某工厂向用户提供发动机，按合同规定，其交货数量和日期是：第一季度末交40台，第二季末交60台，第三季末交80台。工厂的最大生产能力为每季100台，每2季的生产费用是f(x)=50x+（元），此处x为该季生产发动机的台数。若工厂生产的多，多余的发动机可移到下季向用户交货，这样，工厂就需支付存贮费，每台发动机每季的存贮费为4元。问该厂每季应生产多少台发动机，才能既满足交货合同，又使工厂所花费的费用最少（假定第一季度开始时发动机无存货）。 4. 用Matlab的非线性规划命令fmincon求解飞行管理问题的模型二。 5. 用罚函数法求解飞行管理问题的模型二。 -54-
6. 求下列问题的解 22 maxf(x)=2x+3x+3x+x+x 1122322 . x+2x+x+2x+x≤10 1122322 x+x+x+x−x≤50 112232 2x+x+2x+x≤40 11232 x+x=2 13 x+2x≥1 12 x≥0,x,x不约束 1237． 图3 影院剖面示意图 图3为影院的剖面示意图，座位的满意程度主要取决于视角α和仰角β。视角α是观众眼睛到屏幕上、下边缘视线的夹角，α越大越好；仰角β是观众眼睛到屏幕上边缘视线与水平线的夹角，β太大使人的头部过分上仰，引起不舒服感，一般要求β不超过o30。 记影院屏幕高h，上边缘距地面高H，地板线倾角θ，第一排和最后一排座位与屏幕水平距离分别为d和D，观众平均座高为c（指眼睛到地面的距离）。已知参数h=，H=5，d=，D=19，c=（单位：m）。 o（1）地板线倾角θ=10，问最佳座位在什么地方？ o（2）求地板线倾角（一般不超过20），使所有观众的平均满意程度最大。 （3）地板线设计成什么形状可以进一步提高观众的满意程度。 -55-
