第三章 非线性规划

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

2、决策变量为xi=n1,决定投资第i个项目i个项目n0，决定不投资第，i=1, ,n则投资总额为aixi，投资总收益为bixi。因为该公司至少要对一个项目投资，并且总i=1i=1的投资金额不能超过总资金A，故有限制条件n0<ai=1ixiA另外，由于xi(i=1, ,n)只取值0或1，所以还有xi(1-xi)=0,i=1, ,n.最佳投资方案应是投资额最小而总收益最大的方案，所以这个最佳投资决策问题归结为总资金以及决策变量（取0或1）的限制条件下，极大化总收益和总投资之比。因此，其数学模型为：nbxiimaxQ=i=1niai=1xins.t.0<ai=1ixiA (3.1)xi(1

3、-xi)=0,i=1, ,n.上面例题是在一组等式或不等式的约束下，求一个函数的最大值（或最小值）问题，其中目标函数或约束条件中至少有一个非线性函数，这类问题称之为非线性规划问题，简记为（NP）。1.2 非线性规划模型的一般形式minf(x)s.t.hj(x)0,j=1, ,q (3.2)gi(x)=0,其中x=x1 Ti=1, ,p xn称为模型（NP）的决策变量，f称为目标函数，gi(i=1, ,p)和hj(j=1, ,q)称为约束函数。另外，gi(x)=0 (i=1, ,p)称为等式约束，hj(x)0(j=1, ,q)称为不等式约束。特别称min f(x),xE (3.3) n为无约束极

4、值问题。称 min 12xHx+fTTx, s.t. Axb (3.4)为二次规划问题。其中目标函数为自变量x的二次函数，约束条件则全是线性函数。二、非线性规划的求解方法2.1 非线性规划解的概念定义3.1 把满足问题（3.2）中条件的解XEn称为可行解（或可行点），所有可行点的集合称为可行集（或可行域）。记为D。即D=X|hj(x)0,gi(x)=0,xEn 问题(3.2)可简记为 minf(x)。xD定义3.2 对于问题(3.2)，设 X*D，若存在>0,使得对一切XD，且X-X*D，都有fX()<f(X)，则称X*是f(x)在D上的局部极小值点*（局部最优解）。特别地当XX*

5、时，若f(X局部极小值点（严格局部最优解）。)<f(X)，则称X是f(x)在D上的严格定义3.3 对于问题(3.2),设X*D,对任意的XD,都有f(X*)<f(X)则称X*是。特别地当XX时，若f(Xf(x)在D上的全局极小值点（全局最优解）则称X是f(x)在D上的严格全局极小值点（严格全局最优解）。*)<f(X)，2.2 非线性规划的基本解法 2.2.1无约束极值问题求解方法无约束极值问题的求解方法主要有：一维搜索算法；最速下降法；牛顿法和拟牛顿法。 2.2.2 有约束问题求解方法 (1) 近似规划法近似规划法的基本思想是将非线性规划问题中的目标函数和约束条件近似为线性函

6、数，并限制变量的取值范围，从而得到一个近似线性规划问题。对近似线性规划问题求解可得原问题的一个近似解，从这个近似解出发，重复以上步骤，产生一个由线性规划最优解组成的序列。这样的序列往往收敛于非线性规划问题的最有解。 (2) 罚函数法罚函数法的基本思想是通过构造罚函数将非线性规划问题的求解，转化为求解一系列无约束极值问题，这类方法也称为序列无约束最小化方法，简记为SUMT (SequentialUnconstrained Minimization Technique)。该方法主要有两种形式，一种叫外罚函数法，另一种叫内罚函数法。三、非线性规划的Matlab解法3.1有约束问题Matlab中非线性

7、规划的数学模型写成以下形式:minf(x)AxBAeqx=Beq (3.5)C(x)0Ceq(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*XB,Aeq*X=Beq，如果没有等式约束，则A=,B=,Aeq=,Beq=；LB和UB是变量x的下界和上界，如果上界和下界没有约束，则LB=，

8、UB=，如果x无下界，则LB=-inf，如果x无上界，则UB=inf；NONLCON是用M文件定义的非线性向量函数C(x),Ceq(x)；OPTIONS定义了优化参数，可以使用Matlab缺省的参数设置。例2 求下列非线性规划问题min f(x)=x1+x2+8x12-x202-x1-x2+2=0 (3.6)x,x0.1222解 （i）编写M文件fun1.mfunction f=fun1(x);f=x(1)2+x(2)2+8;和M文件fun2.mfunction g,h=fun2(x);g=-x(1)2+x(2);h=-x(1)-x(2)2+2; %等式约束（ii）在Matlab的命令窗口依次

9、输入options=optimset;x,y=fmincon('fun1',rand(2,1),zeros(2,1), . 'fun2', options)x =11y =103.2 无约束问题Matlab中无约束极值问题写成以下形式minf(x) x其中x是一个向量，f(x)是一个标量函数。Matlab中的基本命令是X,FVAL=FMINUNC(FUN,X0,OPTIONS,P1,P2, .)它的返回值是向量x的值和函数的极小值。FUN 是一个M文件，当FUN只有一个返回值时，它的返回值是函数f(x)；当FUN有两个返回值时，它的第二个返回值是f(x)的一阶导

10、数行向量；当FUN有三个返回值时，它的第三个返回值是f(x)的二阶导数阵（Hessian阵）。X0是向量x的初始值，OPTIONS是优化参数，使用缺省参数时，OPTIONS为空矩阵。P1，P2是可以传递给FUN的一些参数。例3 求函数f(x)=100(x2-x1)+(1-x1)的最小值。解：编写M文件fun2.m如下：function f=fun3(x);f=100*(x(2)-x(1)2)2+(1-x(1)2;在Matlab命令窗口输入5 222x=fminunc('fun2',rand(1,2) x =1.0000 1.0000求多元函数的极值也可以使用Matlab的命令X

11、,FVAL= FMINSEARCH(FUN,X0,OPTIONS,P1,P2,.)3.3 罚函数法利用罚函数法，可将非线性规划问题的求解，转化为求解一系列无约束极值问题 考虑如下问题：minf(x)gi(x)0, i=1, ,r,s.t. hi(x)0, i=1, ,s,k(x)=0, i=1, ,t.i取一个充分大的数 M>0 ，构造函数rstP(x,M)=f(x)+Mmax(i=1gi(x),0)-Mmin(i=1hi(x),0)+M|ki=1i(x)|（或P(x,M)=f(x)+M1max(G(x),0)+M2min(H(x),0)+M3|K(x)|g1(x)h1(x)k1(x)这

12、里 G(x)= ，H(x)= ，K(x)= ，M1,M2,M3为适当的行向gr(x)kt(x)hs(x)量，Matlab中可以直接利用 max和 min函数。）则以增广目标函数P(x,M)为目标函数的无约束极值问题minP(x,M)的最优解x也是原问题的最优解。例4 用罚函数法求例2中的非线性规划问题 解 首先，编写 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+2);其次，在Matlab命令

13、窗口输入x,y=fminunc('test',rand(2,1)x =0.82140.4447y =4.9050e+0043.4 二次规划若某非线性规划的目标函数为自变量x的二次函数，约束条件又全是线性的，就称这种规划为二次规划。Matlab中二次规划的数学模型可表述如下：min 12xHx+fTTxs.t. AxbAeqx=beqLBxUB (3.7)这里H是实对称矩阵，f,b,beq,LB,UB是列向量，A,Aeq是相应维数的矩阵。Matlab中求解二次规划的命令是X,FVAL= QUADPROG(H,f,A,b,Aeq,beq,LB,UB,X0,OPTIONS)X的返回值

14、是向量x，FVAL的返回值是目标函数在X处的值。例5 求解二次规划 min f(x)=2x12-4x1x2+4x22-6x1-3x2 x1+x23 (3.8) 4x1+x29 x,x012解 在Matlab命令窗口输入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)Warning: Large-scale method does not currently solve this problem formulation,switching to medium-scale method.> In q

15、uadprog at 242Optimization terminated. x =1.9500 1.0500 value =-11.0250四、非线性规划建模实例4.1 选址问题某公司有6个建筑工地要开工,每个工地的位置(用平面坐标系a,b表示,距离单位:km)及水泥日用量d(单位:t)由下表给出.现计划建两个临时料场,日储量各有20t.假设从料场到工地之间均有直线道路相连.试确定料场的位置,使各料场对各建筑工地的运输量与路程乘积之和为最小。工地位置(a,b)及水泥日用量d4.1.1 模型建立记工地的位置为(ai,bi),水泥的日用量为di,i=1,6；料场的位置为(xj,yj),日储量为e

16、j,j=1,2；从料场j向工地i的运送量为Xij。在各工地用量必须满足和各料场运送量不超过日储量的条件下,使总的吨千米数最小,这是非线性规划问题。可建立数学模型如下：62minf=i=1j=1Xijxj-ai+yj-bi222Xij=di,j=16s.t. Xijej,i=1Xij0i=1,2 ,6j=1,2 (3.9)4.1.2 模型求解利用MATLAB编程（留给读者）求解得：两个料场的坐标分别为：(6.3875,4.3943),(5.7511,7.1867), 总的吨千米数最小为105.4626。由料场A,B向6个工地运料方案见下表：4.2 飞行管理问题在约10，000m高空的某边长160

17、km的正方形区域内，经常有若干架飞机作水平飞行。区域内每架飞机的位置和速度向量均由计算机记录其数据，以便进行飞行管理。当一架欲进入该区域的飞机到达区域边缘时，记录其数据后，要立即计算并判断是否会与区域内的飞机发生碰撞。如果会碰撞，则应计算如何调整各架（包括新进入的）飞机飞行的方向角，以避免碰撞。现假定条件如下：1）不碰撞的标准为任意两架飞机的距离大于8km; 2）飞机飞行方向角调整的幅度不应超过30度； 3）所有飞机飞行速度均为每小时800km;4）进入该区域的飞机在到达区域边缘时，与区域内飞机的距离应在60km以上； 5）最多需考虑6架飞机；6）不必考虑飞机离开此区域后的状况。请你对这个避免

18、碰撞的飞行管理问题建立数学模型，列出计算步骤，对以下数据进行计算（方向角误差不超过0.01度），要求飞机飞行方向角调整的幅度尽量小。设该区域4个顶点的座标为(0,0)，(160,0)，(160,160)，(0,160)。记录数据为： 飞机编号 横座标x 纵座标y 方向角（度） 1 150 140 243 2 85 85 236 3 150 155 220.54 145 50 159 5 130 150 230 新进入 0 0 52 注：方向角指飞行方向与x轴正向的夹角。4.2.1 问题分析首先我们来研究两架飞机不碰撞的条件：记两架飞机的初始位置为：(xi0,yi0),(x0,yj)。j时刻t飞

19、机的位置为：xi(t)=xi+vtcosi,0y(t)=y+vtsini,ii(i=i,j)两架飞机的距离rij=(xi(t)-xj(t)222+(yi(t)-yj(t)22=v(cosi-cosj)+(sini-j)ti0j0i22+2v(x-x)(cosi-cosj)+(y-y)(sini-sinj)t+(xi-xj)+(yi-yj),220j引入记号aij=v(cosi-cosj)+(sini-j),bij=2v(x-x)(cosi-cosj)+(y-y)(sini-sinj),0i0j0i0j222则两架飞机的距离可表示为：rij(t)=aijt+bijt+rij(0)222因此两架飞

20、机不碰撞的条件是：rij(t)=aijt+bijt+rij(0)>64。由于不必考虑在区域外的碰撞，所以两架飞机都在区域中的时间为tij=min(ti,tj) 其中他ti为第i架飞机在区域内的时间：22200D-xi0D-yiyi3,if0i<,tanior<i<2,-tani,00vcosi2D-xi2D-xi00D-x0D-yiD-yii,if0<i,taniori<,-tani,00vsin2D-x2xiiiti= 000-xi,if<,-tanD-yior<3,tanyi,iiii00vcos2xi2xii000yiyi33-yivsin

21、,if<i2,tanix0or2i<2,-taniD-x0iii4.2.2 模型建立记第i架飞机的初始方向角为i，调整后的方向角为i=i+i，其中i为调整N角度。则目标函数为总的调整量：f=i=1i结合前面的不碰撞条件可得下述非线性规划模型：Nminf=i,i=12s.t.rij(t)>64,t<tij,i,j=1, ,N,ij ,i=1, ,Ni64.2.3 模型求解N首先将目标函数改为：f=i=12i,其次考虑到区域对角线的长度为2D，从而任一架飞机在区域内停留的时间不会超过tm=22Dv，所以约束条件中的t<tij可修改为：drij(t)dt22t<t

22、m。注意到rij(t)是t的二次函数，可以利用222=0，求得t=-bij2aij。若0<t<tij，则代人rij(t)=aijt+bijt+rij(0)即可求得rij(t)。下面是求解飞行管理问题的Matlab程序： 先写两个函数文件如下：function f=air1(x) %目标函数f=x*x'function c ceq=air2(x) %非线性约束函数x0=150 85 150 145 130 0;y0=140 85 155 59 159 0;alpha0=243 236 220.5 159 230 52*pi/180;v=800;D=160;co=cos(alpha0+x);si=sin(alpha0+x);tm=sqrt(2)*D/v;for i=2:6for j=1:i-1a(i,j)=v2*(co(i)-co(j)2+(si(i)-si(j)2); b(i,j)=2*v*(x0(i)-x0(j)*(co(i)-co(j)+(y0(i)-y0(j)*(si(i)-si(j);r(i,j)=(x