第三章 非線性規(guī)劃

1、第三章 非線性規(guī)劃一、非線性規(guī)劃模型的建立1.1 非線性規(guī)劃模型簡介如果目標函數(shù)或約束條件中包含非線性函數(shù)，就稱這種規(guī)劃問題為非線性規(guī)劃問題。一般說來，解非線性規(guī)劃要比解線性規(guī)劃問題困難得多。而且，也不像線性規(guī)劃有單純形法這一通用方法，非線性規(guī)劃目前還沒有適于各種問題的一般算法，各個方法都有自己特定的適用范圍。下面通過實例歸納出非線性規(guī)劃數(shù)學(xué)模型的一般形式，介紹有關(guān)非線性規(guī)劃的基本概念。例1 （投資決策問題）某企業(yè)有n個項目可供選擇投資，并且至少要對其中一個項目投資。已知該企業(yè)擁有總資金A元，投資于第i(i=1, ,n)個項目需花資金ai元，并預(yù)計可收益bi元。試選擇最佳投資方案。解: 設(shè)投資

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.最佳投資方案應(yīng)是投資額最小而總收益最大的方案，所以這個最佳投資決策問題歸結(jié)為總資金以及決策變量（取0或1）的限制條件下，極大化總收益和總投資之比。因此，其數(shù)學(xué)模型為：nbxiimaxQ=i=1niai=1xins.t.0<ai=1ixiA (3.1)xi(1

3、-xi)=0,i=1, ,n.上面例題是在一組等式或不等式的約束下，求一個函數(shù)的最大值（或最小值）問題，其中目標函數(shù)或約束條件中至少有一個非線性函數(shù)，這類問題稱之為非線性規(guī)劃問題，簡記為（NP）。1.2 非線性規(guī)劃模型的一般形式minf(x)s.t.hj(x)0,j=1, ,q (3.2)gi(x)=0,其中x=x1 Ti=1, ,p xn稱為模型（NP）的決策變量，f稱為目標函數(shù)，gi(i=1, ,p)和hj(j=1, ,q)稱為約束函數(shù)。另外，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)為二次規(guī)劃問題。其中目標函數(shù)為自變量x的二次函數(shù)，約束條件則全是線性函數(shù)。二、非線性規(guī)劃的求解方法2.1 非線性規(guī)劃解的概念定義3.1 把滿足問題（3.2）中條件的解XEn稱為可行解（或可行點），所有可行點的集合稱為可行集（或可行域）。記為D。即D=X|hj(x)0,gi(x)=0,xEn 問題(3.2)可簡記為 minf(x)。xD定義3.2 對于問題(3.2)，設(shè) X*D，若存在>0,使得對一切XD，且X-X*D，都有fX()<f(X)，則稱X*是f(x)在D上的局部極小值點*（局部最優(yōu)解）。特別地當XX*

5、時，若f(X局部極小值點（嚴格局部最優(yōu)解）。)<f(X)，則稱X是f(x)在D上的嚴格定義3.3 對于問題(3.2),設(shè)X*D,對任意的XD,都有f(X*)<f(X)則稱X*是。特別地當XX時，若f(Xf(x)在D上的全局極小值點（全局最優(yōu)解）則稱X是f(x)在D上的嚴格全局極小值點（嚴格全局最優(yōu)解）。*)<f(X)，2.2 非線性規(guī)劃的基本解法 2.2.1無約束極值問題求解方法無約束極值問題的求解方法主要有：一維搜索算法；最速下降法；牛頓法和擬牛頓法。 2.2.2 有約束問題求解方法 (1) 近似規(guī)劃法近似規(guī)劃法的基本思想是將非線性規(guī)劃問題中的目標函數(shù)和約束條件近似為線性函

6、數(shù)，并限制變量的取值范圍，從而得到一個近似線性規(guī)劃問題。對近似線性規(guī)劃問題求解可得原問題的一個近似解，從這個近似解出發(fā)，重復(fù)以上步驟，產(chǎn)生一個由線性規(guī)劃最優(yōu)解組成的序列。這樣的序列往往收斂于非線性規(guī)劃問題的最有解。 (2) 罰函數(shù)法罰函數(shù)法的基本思想是通過構(gòu)造罰函數(shù)將非線性規(guī)劃問題的求解，轉(zhuǎn)化為求解一系列無約束極值問題，這類方法也稱為序列無約束最小化方法，簡記為SUMT (SequentialUnconstrained Minimization Technique)。該方法主要有兩種形式，一種叫外罰函數(shù)法，另一種叫內(nèi)罰函數(shù)法。三、非線性規(guī)劃的Matlab解法3.1有約束問題Matlab中非線性

7、規(guī)劃的數(shù)學(xué)模型寫成以下形式:minf(x)AxBAeqx=Beq (3.5)C(x)0Ceq(x)=0其中f(x)是標量函數(shù)，A,B,Aeq,Beq是相應(yīng)維數(shù)的矩陣和向量，C(x),Ceq(x)是非線性向量函數(shù)。Matlab中的命令是:X=FMINCON(FUN,X0,A,B,Aeq,Beq,LB,UB,NONLCON,OPTIONS)它的返回值是向量x，其中FUN是用M文件定義的函數(shù)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文件定義的非線性向量函數(shù)C(x),Ceq(x)；OPTIONS定義了優(yōu)化參數(shù)，可以使用Matlab缺省的參數(shù)設(shè)置。例2 求下列非線性規(guī)劃問題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)是一個標量函數(shù)。Matlab中的基本命令是X,FVAL=FMINUNC(FUN,X0,OPTIONS,P1,P2, .)它的返回值是向量x的值和函數(shù)的極小值。FUN 是一個M文件，當FUN只有一個返回值時，它的返回值是函數(shù)f(x)；當FUN有兩個返回值時，它的第二個返回值是f(x)的一階導(dǎo)

10、數(shù)行向量；當FUN有三個返回值時，它的第三個返回值是f(x)的二階導(dǎo)數(shù)陣（Hessian陣）。X0是向量x的初始值，OPTIONS是優(yōu)化參數(shù)，使用缺省參數(shù)時，OPTIONS為空矩陣。P1，P2是可以傳遞給FUN的一些參數(shù)。例3 求函數(shù)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求多元函數(shù)的極值也可以使用Matlab的命令X

11、,FVAL= FMINSEARCH(FUN,X0,OPTIONS,P1,P2,.)3.3 罰函數(shù)法利用罰函數(shù)法，可將非線性規(guī)劃問題的求解，轉(zhuǎn)化為求解一系列無約束極值問題 考慮如下問題：minf(x)gi(x)0, i=1, ,r,s.t. hi(x)0, i=1, ,s,k(x)=0, i=1, ,t.i取一個充分大的數(shù) M>0 ，構(gòu)造函數(shù)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為適當?shù)男邢騡r(x)kt(x)hs(x)量，Matlab中可以直接利用 max和 min函數(shù)。）則以增廣目標函數(shù)P(x,M)為目標函數(shù)的無約束極值問題minP(x,M)的最優(yōu)解x也是原問題的最優(yōu)解。例4 用罰函數(shù)法求例2中的非線性規(guī)劃問題 解 首先，編寫 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 二次規(guī)劃若某非線性規(guī)劃的目標函數(shù)為自變量x的二次函數(shù)，約束條件又全是線性的，就稱這種規(guī)劃為二次規(guī)劃。Matlab中二次規(guī)劃的數(shù)學(xué)模型可表述如下：min 12xHx+fTTxs.t. AxbAeqx=beqLBxUB (3.7)這里H是實對稱矩陣，f,b,beq,LB,UB是列向量，A,Aeq是相應(yīng)維數(shù)的矩陣。Matlab中求解二次規(guī)劃的命令是X,FVAL= QUADPROG(H,f,A,b,Aeq,beq,LB,UB,X0,OPTIONS)X的返回值

14、是向量x，F(xiàn)VAL的返回值是目標函數(shù)在X處的值。例5 求解二次規(guī)劃 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四、非線性規(guī)劃建模實例4.1 選址問題某公司有6個建筑工地要開工,每個工地的位置(用平面坐標系a,b表示,距離單位:km)及水泥日用量d(單位:t)由下表給出.現(xiàn)計劃建兩個臨時料場,日儲量各有20t.假設(shè)從料場到工地之間均有直線道路相連.試確定料場的位置,使各料場對各建筑工地的運輸量與路程乘積之和為最小。工地位置(a,b)及水泥日用量d4.1.1 模型建立記工地的位置為(ai,bi),水泥的日用量為di,i=1,6；料場的位置為(xj,yj),日儲量為e

16、j,j=1,2；從料場j向工地i的運送量為Xij。在各工地用量必須滿足和各料場運送量不超過日儲量的條件下,使總的噸千米數(shù)最小,這是非線性規(guī)劃問題。可建立數(shù)學(xué)模型如下：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), 總的噸千米數(shù)最小為105.4626。由料場A,B向6個工地運料方案見下表：4.2 飛行管理問題在約10，000m高空的某邊長160

17、km的正方形區(qū)域內(nèi)，經(jīng)常有若干架飛機作水平飛行。區(qū)域內(nèi)每架飛機的位置和速度向量均由計算機記錄其數(shù)據(jù)，以便進行飛行管理。當一架欲進入該區(qū)域的飛機到達區(qū)域邊緣時，記錄其數(shù)據(jù)后，要立即計算并判斷是否會與區(qū)域內(nèi)的飛機發(fā)生碰撞。如果會碰撞，則應(yīng)計算如何調(diào)整各架（包括新進入的）飛機飛行的方向角，以避免碰撞?，F(xiàn)假定條件如下：1）不碰撞的標準為任意兩架飛機的距離大于8km; 2）飛機飛行方向角調(diào)整的幅度不應(yīng)超過30度； 3）所有飛機飛行速度均為每小時800km;4）進入該區(qū)域的飛機在到達區(qū)域邊緣時，與區(qū)域內(nèi)飛機的距離應(yīng)在60km以上； 5）最多需考慮6架飛機；6）不必考慮飛機離開此區(qū)域后的狀況。請你對這個避免

18、碰撞的飛行管理問題建立數(shù)學(xué)模型，列出計算步驟，對以下數(shù)據(jù)進行計算（方向角誤差不超過0.01度），要求飛機飛行方向角調(diào)整的幅度盡量小。設(shè)該區(qū)域4個頂點的座標為(0,0)，(160,0)，(160,160)，(0,160)。記錄數(shù)據(jù)為： 飛機編號 橫座標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。由于不必考慮在區(qū)域外的碰撞，所以兩架飛機都在區(qū)域中的時間為tij=min(ti,tj) 其中他ti為第i架飛機在區(qū)域內(nè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，調(diào)整后的方向角為i=i+i，其中i為調(diào)整N角度。則目標函數(shù)為總的調(diào)整量：f=i=1i結(jié)合前面的不碰撞條件可得下述非線性規(guī)劃模型：Nminf=i,i=12s.t.rij(t)>64,t<tij,i,j=1, ,N,ij ,i=1, ,Ni64.2.3 模型求解N首先將目標函數(shù)改為：f=i=12i,其次考慮到區(qū)域?qū)蔷€的長度為2D，從而任一架飛機在區(qū)域內(nèi)停留的時間不會超過tm=22Dv，所以約束條件中的t<tij可修改為：drij(t)dt22t<t

22、m。注意到rij(t)是t的二次函數(shù)，可以利用222=0，求得t=-bij2aij。若0<t<tij，則代人rij(t)=aijt+bijt+rij(0)即可求得rij(t)。下面是求解飛行管理問題的Matlab程序： 先寫兩個函數(shù)文件如下：function f=air1(x) %目標函數(shù)f=x*x'function c ceq=air2(x) %非線性約束函數(shù)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