matlab與數(shù)學(xué)建模-第03章非線性

/第三章§1例1（投資決策問題）某企業(yè)有n個(gè)項(xiàng)目可供選擇投資，并且至少要對(duì)其中一個(gè) A元，投資于第i(i1,L,n)個(gè)項(xiàng)目需花 并預(yù)計(jì)可收益bi元。試選擇最佳投資方案。解設(shè)投資決策變量為

決定投資第i

，i1,L,n0，決定不投資第i個(gè)項(xiàng) 則投資總額為aixi，投資總收益為bixi。因?yàn)樵摴局辽僖獙?duì)一個(gè)項(xiàng)目投資，i in0aixiixi(1xi)0,i1,L,結(jié)為總以及決策變量（取0或1）的限制條件下，極大化總收益和總投資之比。因此，nbinnmaxQnais.t.0aixiixi(1xi)0,i1,L,minf( hj(x)gi(x)

j1,L,qi1,L,p

xx1Lxn]T稱為模型（NP）f稱為目標(biāo)函數(shù)，gi(i1,Lp)hjj1,Lq)稱為約束函數(shù)。另外，gix)0(i1,L,p)稱為等式約束，hjx)0(j1,Lq)稱為不等式的約束。非線性規(guī)劃的中非線性規(guī)劃的數(shù)學(xué)模型寫成以下形式minf(AxAeqxC(x) A,B,Aeq,Beq定義了線性約束A*XB,Aeq*XBeq，如果沒有線性約束，則A=[],B=[],Aeq=[],Beq=[]；LB和UB是變量x的下界和上界，如果上界和下界沒有約束，則LB=[]，UB=[]，如x無(wú)下界，則LB各分量都為-inf，如果x無(wú)上UB定義了優(yōu)化參數(shù)，可以使用缺省的參數(shù)設(shè)置。 minf(x)x2x2x2 13 x2x2x213xx2x3 2x1x223x22x23x1,x2,x3解（i）Mfun1.m編寫M文件fun2.mx(1)+x(2)^2+x(3)^3-20];%非線性不等式約束;%編寫主程序文件example2.m如下：'fun2',options)記（NP）的可行域?yàn)镵。x*Kf(x*)f( x則稱x*是（NP）的整體最優(yōu)解，f(x*)是(NP)的整體最優(yōu)值。如果有f(x*)f( xK,xx*Kx*N(x*，使f(x*)f( xN(x*)IK則稱x*是（NP）的局部最優(yōu)解，f(x*)是(NP)的局部最優(yōu)值。如果有f(x*)f( xN(x*)I則稱x*是（NP）的嚴(yán)格局部最優(yōu)解，f(x*)是(NP)的嚴(yán)格局部最想是：從一個(gè)選定的初始點(diǎn)x0Rn{xk}，使得當(dāng){xk是有窮點(diǎn)列時(shí)，其最后一個(gè)點(diǎn)是(NP)的最優(yōu)解；當(dāng){xk是無(wú)窮點(diǎn)xkRn是某迭代方法的第kxk1Rn是第k1xk1xktk 這里tkR1pkRn,pk1pkxkxk1的方向。式（1）通常，把基本迭代格式（1）中的pk稱為第k輪搜索方向，tk為沿pk方向的xRnp0，若存在0fxtpfxt(0,，稱向量pfx處的下降方向。xRnp0，若存在t0xtpK向，稱之為函數(shù)fx處關(guān)于K的可行下降方向。現(xiàn)在，給出用基本迭代格式（1）求解(NP)的一般步驟如下0°x0，令k:01°fxkK的可行下降方向作為搜索方向pk。2xkpk尋求適當(dāng)?shù)牟介L(zhǎng)tk，使目標(biāo)函數(shù)值3°求出下一個(gè)迭代點(diǎn)。按迭代格式（1）xk1xktkpk若xk1已滿足某種終止條件，停fx為定義在n維歐氏空間E(n)中某個(gè)凸集R上的函數(shù)，若對(duì)任何實(shí)數(shù)(01)R中的任意x(1)x(2)f(x(1)(1)x(2))f(x(1))(1)f(x(2)若對(duì)每一個(gè)(01)x(1x(2Rf(x(1)(1)x(2))f(x(1))(1)f(x(2)

f(R{x|gj(x)0,j1,2,L,fx為凸函數(shù)，gjxj1,2,L,l為凸函數(shù)，這樣的非線性規(guī)劃稱為 （2（at

f 若f(t)是[a,b]區(qū)間上的下單峰函數(shù)，介紹通過不斷地縮短[a,b]的長(zhǎng)度，來(lái)為了縮短區(qū)間[a,b]，逐步搜索得（2）的最優(yōu)解t*的近似值，可以采用以下途徑：在[ab中任取兩個(gè)關(guān)于[ab是對(duì)稱的點(diǎn)t1和t2（不妨設(shè)t2t1則必有t*[at]，因而[at]是縮短了的單峰區(qū)間；若f(t)f(t)，則有t*t*[tb

f(t)f(t，則[at

和[tb FibonacciF0F1FnFn2Fn1 n數(shù)列{FnFibonacciFn稱為第nFibonacciFibonacci數(shù)之Fn1Fibonacci分?jǐn)?shù)。FnFn1，其后各次分別為Fn2Fn3,LF1。由此，若t和t(tt是單峰區(qū)間[a

2 t1aFn1，t2aFn2ba ba taFn1(ba)

aFn2(b

F F 它們關(guān)于[a,b]確是對(duì)稱的點(diǎn)過精度0，這就要求最后區(qū)間的長(zhǎng)度不超過，即ba

應(yīng)按照預(yù)先給定的精度，確定使（4）成立的最小整數(shù)n作為搜索次數(shù)，直到進(jìn)行到第n個(gè)探索點(diǎn)時(shí)停止。f(t的單峰區(qū)間[ab的辦法，來(lái)求得問題（2）的近似解，是Kiefer(1953年)，叫做Finbonacci法，具體步驟如下：1o選取初始數(shù)據(jù)，確定單峰區(qū)間[a0b0，給出搜索精度0，由（4）確定搜索次數(shù)n。2ok1,aabb，計(jì)算最初兩個(gè)搜索點(diǎn)，按（3）計(jì)算t和t 3owhileknf1f(t1),f2f(t2iff1at;tt;taF(n1k)(b

2 1

F(nkF(n1k

bt;tt;tb (a1 2 F(nkkk4o當(dāng)進(jìn)行至kn1t1

1(ab) t22(ata(

)(b t1和t2這兩點(diǎn)中，以函數(shù)值較小者為近似極小點(diǎn)，相應(yīng)的函數(shù)值為近似極小值。并得最終區(qū)間[a,t1或[t2b]。3f(tt2t2的近似極小點(diǎn)，要求縮短后的區(qū)間不大于區(qū)間[1,3]的0.08倍。 0.618若01 稱之為黃金分割數(shù)，其值為

510.6180339887L2黃金分割數(shù)FibonaccilimFn1n0.6182個(gè)探索點(diǎn)開始每增加一個(gè)探索點(diǎn)作一輪迭代以后，原單峰區(qū)間0.618倍。計(jì)算n個(gè)探索點(diǎn)的函數(shù)值可以把原區(qū)間[a0,b0]連續(xù)n1次，因?yàn)槊看蔚目s短率均為，故最后的區(qū)間長(zhǎng)度為

)n這就是說，當(dāng)已知縮短的相對(duì)精度為時(shí)，可用下式計(jì)算探索點(diǎn)個(gè)數(shù)nn10.618minf( xE(n xk1xktk 告訴，點(diǎn)xk的負(fù)梯度方向pkf(xk)xkf下降最快的方向。為此，稱負(fù)梯度方向fxkfxk處的時(shí)，用fxk0fxk)1°選取初始數(shù)據(jù)。選取初始點(diǎn)x0，給定終止誤差，令k:02°求梯度向量。計(jì)算fxkf(xk3°pkf(xk)4°進(jìn)行一維搜索。求tkf(xktpk)minf(xktpk

txk1xktkpkk:k1,2例 minf(x)x2 xx1x2)Tx02,2)T。解（i）fx)(2x150x2TM文件detaf.mf(xwhilenorm(g)>0.000001whilef>f0Newton考慮目標(biāo)函數(shù)f在點(diǎn)xk處的二次近f(x)Q(x)

f(xk)f(xk)T(xxk)1(xxk)T2f(xk)(xxk)2f(xk 2f(xk)

x2f(xk) M

n2f(xk

)2f(xk)xx x

由于2fxk)正定，函數(shù)Q的駐點(diǎn)xk1是Qx)的極小點(diǎn)。為求Qxk1fxk2fxkxk1xk0，xk1xk[2f(xk)]1f(xk)對(duì)照基本迭代格式（1）可知從點(diǎn)xk出發(fā)沿搜索方向。pk[2f(xk)]1f(xktk1即可得Qx的最小點(diǎn)xk1。通常，把方向pk叫做從點(diǎn)xk出發(fā)的NewtonNewton方向并取步長(zhǎng)為1的求解方法，稱之為Newton法。其具體步驟如下：1x0，給定終止誤差0，令k:02°求梯度向量。計(jì)算fxk

f(xk

3°構(gòu)Newton方向。計(jì)算[2fxk)]1pk[2f(xk)]1f(xk)4°求下一迭代點(diǎn)。令xk1xkpk,k:k1，轉(zhuǎn)2°。例5 用Newton法求解，minf(x)x425x4x2x1x02,2)T

11 100x32x2x112x2

124x2

12f 300x2122x2

1 1（ii）編寫M文件nwfun.m如下：function[f,df,d2f]=nwfun(x);x,whilef>f0Newton法的優(yōu)點(diǎn)是收斂速度快；缺點(diǎn)是有時(shí)不好用而需采取改進(jìn)措施，此外，當(dāng)維數(shù)較高時(shí)，計(jì)算[2f(xk)]1的工作量很大。變尺度法（VariableMetricAlgorithm）20多年來(lái)發(fā)展起來(lái)的，它不僅是求解尺度法—DFPDavidon1959年提出，后經(jīng)Fletcher和Powell加以改進(jìn)。已經(jīng)知道，法的搜索方向是2fxk)]1fxk)，為了不計(jì)算二階導(dǎo)數(shù)矩陣[2f(xk)]及其逆陣，設(shè)法構(gòu)造另一個(gè)矩陣，用它來(lái)近二階導(dǎo)數(shù)矩陣的逆陣[2f(xk)]1，這一類方法也稱擬法（Quasi-NewtonMethod。下面研究如何構(gòu)造這樣的近似矩陣，并將它記H(k)。要求：每一步都能以矩陣最后應(yīng)收斂于解點(diǎn)處的Hesse陣的逆陣。當(dāng)fx)是二次函Hesse陣為常數(shù)陣Axkxk1處的梯度之f(xk1)f(xk)A(xk1xk或xk1xkA1[f(xk1)f(xkH(k1)滿足xk1xkH(k1)[f(xk1)f(xk G(k)f(xk1)f(xkxkxk1

xkH(k1)G(k) 現(xiàn)假H(k)已知，用下H(k1)（設(shè)H(k)H(k1)均為對(duì)稱正定陣H(k1)H(k)H(k 其中H(k)k次校正矩陣。顯H(k1)Newton條件（9）即要xk(H(k)H(k))G(k或

H(k)G(k)xkH(k)G(kH(k)的一種比較簡(jiǎn)單的

H(k)xk(Q(k))TH(k)G(k)(W(k)其中Q(k)和W(k)為兩個(gè)待將式（12）中的H(k)代入（11）xk(Q(k))TG(k)H(k)G(k)(W(k))TG(k)xkH(k)G(k(Q(k))TG(k)(W(k))TG(k)考慮到H(k)應(yīng)為對(duì)稱陣，最簡(jiǎn)單的辦Q(k)kW(k)H(k)G(k

k(xk)TG(k)k(G(k))TH(k)G(k)若(xk)TGk)和(Gk)THkGk)

(xk)TG(k

(G(k))T

(G(k))TH(k)G(k

H(k)

xk(xk)T

H(k)G(k)(G(k))TH(k(G(k))TH(k)G(k

H(k

H(k

xk(xk)T(G(k))T

H(k)G(k)(G(k))TH(k(G(k))TH(k)G(k

上述矩陣稱為尺度矩陣。通常，取第一個(gè)尺度矩陣H 為單位陣，以后的尺度)若H(k)為對(duì)稱正定陣，則由式（18）產(chǎn)生的H(k1)現(xiàn)將DFP 1°給定初始點(diǎn)x0及梯度允許誤差0。2fx0

H(0)I（單位矩陣p0H(0)f(x0p0方向進(jìn)行一維搜索，確定最佳步長(zhǎng)0minf(x0p0)x1x00

f(x00p04xk，算出fxkf(x0)則xk即為所求的近似解，停止迭代；否則，計(jì)算H(k)H(k

H(k

xk1(xk1 (G(k1))Txk

H(k1)G(k1)(G(k1))TH(k(G(k1))TH(k1)G(kpkH(k)fxkpk方向上進(jìn)行一維搜索，得kkxk1xkk5°若xk1滿足精度要求xk1即為所求的近似解，否則，轉(zhuǎn)回4°，直到求理解，因而在實(shí)際應(yīng)用中常為人們所采用。下面介紹Powell方法。°{p01pn1}。給定終止誤差0，令k:02y0:xk，依次沿{p01pn1中的方向進(jìn)行一維搜索。對(duì)應(yīng)地得到輔助迭代點(diǎn)y1,y2,L,yn，即jyjyj1 pjf(yj1

pj1)

f(yj1tp j1,L, /3°pnyny0，若

xk1yn找出m

f(ym1)f(ym)max{f(yj1)f(yj)|1jf(y0)2f(yn)f(2yny0)2[f(ym1)f(ymxk1yn

pn:f(yntpn)minf(yntpn) t{p0,p1pn1}k1:{p0,L,pm1,pm1,L,pn1,pn}k:k12 minf(x)x其中x可以為標(biāo)量或中fminunc的基本[X,FVAL]=FMINUNC(FUN,X0,OPTIONS,P1,P2,FUNfxFUN有三個(gè)返回例6 求函數(shù)f(x)100(xx2)2(1x)2的最小值。 解：編寫Mfun2.m如下：function[f,g]=fun2(x);g=[-400*x(1)*(x(2)-x(1)^2)-2*(1-x(1));200*(x(2)-編寫主函數(shù)文件example6.m如下：options=optimset('GradObj','on');function[f,df,d2f]=fun3(x); -編寫主函數(shù)文件example62.moptions=optimset('GradObj','on','Hessian','on'); 例7 求函數(shù)f(x)sin(x)3取最小值時(shí)的x值。解編寫f(x)的M文件fun4.m如下：functionf=fun4(x);編寫主函數(shù)文件example7.m如下： 中二次規(guī)劃的數(shù)學(xué)模型可表述如下：

1xTHxfTx,AxAeqx中求解二次規(guī)劃令是[X,FVAL]= 看在指令中運(yùn)行helpqurog后的幫助。8minf(x)2x24xx4x26x

x2

1 4x1x2 x1,x2解h=[4,-4;-f=[-6;-x x

Minfx)11.0250SUMT(SequentialUnconstrainedMinizationTechnique)。minf(gi(x)0,i1,L,s.t.hj(x)0,j1,L,k (x)0,m1,L,k取一個(gè)充分大的數(shù)M0 P(x,M)f(x)Mmax(gi(x),0)Mmin(hi(x),0)M|ki(x)

G(x)

H(x)

（或P(x,M)f(x)Msummax Msummin M||K(x) 這里G(x)g1 gr(x)，H(x)h1 hsK(x)k1(x) kt(t)， 增廣目標(biāo)函數(shù)PxM為目標(biāo)函數(shù)的無(wú)約束極值問題例9求下列非線性規(guī)劃minf(x)x2x2 x2x 2 x1

2 x,x 解(i)Mtest.mfunctiong=test(x); g=f-M*sum(min([x';zeros(1,2)]))-M*min(x(1)^2-M*abs(-x(1)- fminbndminf(x)x

x[x1,x2它的返回值是極小點(diǎn)x和函數(shù)的極小值。這里fun是用M文件定義的函數(shù)或 10求函數(shù)f(x)x3)21,x[0,5]的最小值。解編寫Mfun5.mfunctionf=(x-3)^2- fseminf函數(shù)min{F(x)|C(x)0,Ceq(x)0,PHI(x,w)x

A*xAeq*x其中FUN用于定義目標(biāo)函數(shù)F(x；X0x的初始值；NTHETA是半無(wú)窮約束PHI(xwSEMINFCON用于定義非線性不等式約束C(x)輸入?yún)⒘縓和S，S是的取樣步長(zhǎng)，也許不被使用。例11 求函數(shù)f(x)(x10.5)2(x20.5)2(x30.5)2取最小值時(shí)的x值,K(x,w)sin(wx)cos(wx) (w50)2sin(wx)x 1 1 1 K(x,w)sin(w2x2)cos(w2x) (w50)2sin(wx)x

2 1w1100，1w2解（1）Mfun6.mfunctionf=fun6(x,s);編寫Mfun7.mifs=[0.2,0;0.2%半無(wú)窮k1=sin(w1*x(1)).*cos(w1*x(2))-1/1000*(w1-50).^2-sin(w1*x(3))-x(3)-k2=sin(w2*x(2)).*cos(w2*x(1))-1/1000*(w2-50).^2-sin(w2*x(3))-x(3)- fminimax函數(shù)x A*xAeq*x

C(x)Ceq(x)F(xF1x),LFm(x)}。上述問題的命令為例 求函數(shù)族{f1(x),f2(x),f3(x),f4(x),f5(x)}取極大極小值時(shí)的x值其中f(x)2x2x248x40x f(x)x23x f3(x)x13x2f4(x)x1f5(x)x1x2解（1）Mfun8.mfunctionf=[2*x(1)^2+x(2)^2-48*x(1)--x(1)^2-x(1)+3*x(2)--x(1)-x(1)+x(2)-（2）調(diào)用函數(shù)fminimax例 已知函數(shù)f(x)ex1(4x22x24xx2x1)，且滿足非線性約束1x1x2x1x2xx1求minfx

1 ex1(4x22x24xx8x6x ex1(4x1

1 解（1）Mfun9.mfunction[f,df]=fun9(x);編寫Mfun10.mfunction[c,ceq,dc,dceq]=fun10(x);優(yōu)化工具箱中的optimtool 1優(yōu)化問題用戶圖形界面解法示意例 利用例1已經(jīng)定義好的函數(shù)fun1和fun2。在命令窗口運(yùn)行optimtool，就打后用鼠標(biāo)點(diǎn)一下“start”按鈕，就得到求解結(jié)果，再使用“file”菜單下的“ExporttoWorkspace…”選項(xiàng)，把計(jì)算結(jié)果輸出到工作空間中去。§410，000m160km的正方形區(qū)域內(nèi)，經(jīng)常有若干架飛機(jī)作水平行計(jì)算（方向角誤差不超過0.01度，要求飛機(jī)飛行方向角調(diào)整的幅度盡量小。1飛行記錄數(shù)1345D為飛行管理區(qū)域的邊長(zhǎng)；為飛行管理區(qū)域，取直角坐標(biāo)系使其為[0,D][0,D](x0y0為第i架飛機(jī)的初始 (xi(tyi(t為第i架飛機(jī)在tii為第i

0

為第

根據(jù)相對(duì)運(yùn)動(dòng)的觀點(diǎn)在兩架飛機(jī)i和j的飛行時(shí)，可以將飛機(jī)i視為不動(dòng)而飛機(jī)j以相對(duì)速度vijvjvi(acosjacosi,asinjasiniv2asinji(sinji,cosji

2a

ji

),

ji

ii不妨設(shè)ji，此時(shí)相對(duì)飛行方向角為ij 2

22相對(duì)飛行方向(x0x0)2(x0x0)2(y0y0 0 則只要當(dāng)相對(duì)飛行方向角ij02 2記0為調(diào)整前第j架飛機(jī)相對(duì)于第i架飛機(jī)的相對(duì)速度（矢量）與這兩架線（j指向i的矢量）的夾角（以連線矢量為基準(zhǔn)，逆時(shí)針方向?yàn)檎槙r(shí)針方向?yàn)?1()

0mn＝相對(duì)速度vmn的幅角－從nm的連線矢量的幅0

ein（注意0表達(dá)式中的i表示虛數(shù)單位）這里 角度0（m,n1,2,L,6。6 01()0，i,j1,2,L,6，i i 130y0=[1408515550150q=[243236220.5159230xy0=[x0; form=1:6forif %往純文本文件中寫LINGO數(shù)據(jù)dlmwrite('txt1.txt',b0,'delimiter','\t','newline','PC','-append','roffset',求得0的值如2所示20的12345