北航最優(yōu)化方法大作業(yè)參考

1、11流量工程問題1.1問題重述定義一個(gè)有向網(wǎng)絡(luò)G=(N,E)，其中N是節(jié)點(diǎn)集，E是弧集。令A(yù)是網(wǎng)絡(luò)G的點(diǎn)弧關(guān) 聯(lián)矩陣，即NXE階矩陣， 且第I列與弧里(l,j)對(duì)應(yīng)， 僅第i行元素為1， 第j行元素為-1, 其余元素為0。 再令bm=(bm1,bmN)T，fm=(fm1,fmE)T，則可將等式約束表示成：Afm=bm本算例為一經(jīng)典TE算例。算例網(wǎng)絡(luò)有7個(gè)節(jié)點(diǎn)和13條弧，每條弧的容量是5個(gè)單 位。此外有四個(gè)需求量均為4個(gè)單位的源一目的對(duì)，具體的源節(jié)點(diǎn)、目的節(jié)點(diǎn)信息如圖所 示。這里為了簡(jiǎn)單，省區(qū)了未用到的弧。此外，弧上的數(shù)字表示弧的編號(hào)。此時(shí)，c=(5,5,5)1X13)T，4 -400,如一4

2、040,ba =0 電4D切一4 -00n000o00000004根據(jù)上述四個(gè)約束條件，分別求得四個(gè)情況下的最優(yōu)決策變量圖1網(wǎng)絡(luò)拓?fù)浜土髁啃枨骕=(X12,X13,X75)1X13)。in = 1:亠 2.dsA瞞 m 3: 3- 丄山21.27節(jié)點(diǎn)算例求解利用Matlab編寫對(duì)偶單純形法程序，可求得 最優(yōu)解為x2*=0 4 0 0 0 0 0 0 0 0 0 0 0T對(duì)應(yīng)的最優(yōu)值CTX2=201.2.3算例3(b3=0;-4;4;0;0;0;0T)Mi ni mizeJx3Subject toAx3=b3X3=0利用Matlab編寫對(duì)偶單純形法程序，可求得：最優(yōu)解為x3*=4 0 0 0 4

3、 0 0 0 0 0 0 0 0T對(duì)應(yīng)的最優(yōu)值CTX3=40-111-1-11 1-1-1-11 1-1-1 11.2.1算例1(b仁4;-4;0;0;0;0;0T)轉(zhuǎn)化為線性規(guī)劃問題：Mini mizecTx1Subject to Ax仁b1x1=0利用Matlab編寫對(duì)偶單純形法程序，可求得：最優(yōu)解為x1*=4 0 0 0 0 0 0 0 0 0 0 0 0T對(duì)應(yīng)的最優(yōu)值cTx仁201.2.2算例2(b2=4;0;-4;0;0;0;0T)Mini mizecTx2Subject toAx2=b2X2=031.2.4算例4(b4=4;0;0;0;0;0;-4T)Mi ni mizeJx4Su

4、bject toAx4=b4X4=0利用Matlab編寫對(duì)偶單純形法程序，可求得：最優(yōu)解為x4*=4 0 0 4 0 0 0 0 0 4 0 0 0T對(duì)應(yīng)的最優(yōu)值Jx4=601.3計(jì)算結(jié)果及結(jié)果說明1.3.1算例1(b1=4;-4;0;0;0;0;0T)算例1中，由b1可知，節(jié)點(diǎn)2為需求節(jié)點(diǎn)，節(jié)點(diǎn)1為供給節(jié)點(diǎn)，由節(jié)點(diǎn)1將信息傳 輸至節(jié)點(diǎn)2的最短路徑為弧1。圖2算例1最優(yōu)傳輸示意圖求得的最優(yōu)解為x1*=4 0 0 0 0 0 0 0 0 0 0 0 0T,即只經(jīng)過弧1運(yùn)輸4個(gè)單位流量， 其余弧無(wú)流量。又因?yàn)?，每條弧的費(fèi)用均為5，所以最小費(fèi)用為20。經(jīng)分析，計(jì)算結(jié)果合理可信。1.3.2算例2(b2

5、=4;0;-4;0;0;0;0T)算例2中，由b2可知，節(jié)點(diǎn)3為需求節(jié)點(diǎn)，節(jié)點(diǎn)1為供給節(jié)點(diǎn)，由節(jié)點(diǎn)1將信息傳 輸至節(jié)點(diǎn)2的最短路徑為弧2。4r 12L圖3算例2最優(yōu)傳輸示意圖求得的最優(yōu)解為x2*=0 4 0 0 0 0 0 0 0 0 0 0 0T,即只經(jīng)過弧2運(yùn)輸4個(gè)單位流量, 其余弧無(wú)流量。又因?yàn)?，每條弧的費(fèi)用均為5，所以最小費(fèi)用為20。經(jīng)分析，計(jì)算結(jié)果合理可信。1.3.3算例3(b3=0;-4;4;0;0;0;0T)算例3中，由b3可知，節(jié)點(diǎn)2為需求節(jié)點(diǎn)，節(jié)點(diǎn)3為供給節(jié)點(diǎn)，由節(jié)點(diǎn)3將信息傳 輸至節(jié)點(diǎn)2的最短路徑為弧5-弧1。圖4算例3最優(yōu)傳輸示意圖求得的最優(yōu)解為x3*=4 0 0 0 4

6、 0 0 0 0 0 0 0 0T,即經(jīng)過弧5運(yùn)輸4個(gè)單位流量至節(jié) 點(diǎn)1，再經(jīng)弧1運(yùn)輸4個(gè)單位流量至節(jié)點(diǎn)2,其余弧無(wú)流量。又因?yàn)?，每條弧的費(fèi)用均為5， 所以最小費(fèi)用為40。經(jīng)分析，計(jì)算結(jié)果合理可信。1.3.4算例4(b4=4;0;0;0;0;0;-4T)算例4中，由b4可知，節(jié)點(diǎn)7為需求節(jié)點(diǎn)，節(jié)點(diǎn)1為供給節(jié)點(diǎn)，由節(jié)點(diǎn)1將信息傳 輸至節(jié)點(diǎn)7的最短路徑為弧1-弧4-弧10。流量45弧4涂量4圖5算例3最優(yōu)傳輸示意圖求得的最優(yōu)解為x4*=4 0 0 4 0 0 0 0 0 4 0 0 0T,即經(jīng)過弧1運(yùn)輸4個(gè)單位流量至節(jié) 點(diǎn)2，再經(jīng)弧4運(yùn)輸4個(gè)單位流量至節(jié)點(diǎn)5,最后經(jīng)弧5運(yùn)輸4個(gè)單位流量至節(jié)點(diǎn)7,其

7、 余弧無(wú)流量。又因?yàn)?，每條弧的費(fèi)用均為5，所以最小費(fèi)用為60。經(jīng)分析，計(jì)算結(jié)果合理可信。62重要算法編寫與觀察2.1習(xí)題5.6(a)初值為(0,0)時(shí)本算法令g的2范數(shù)在10-4時(shí)，停止迭代，經(jīng)過86次迭代收斂收斂因子(f(k+1)-f*)/(f(k)-f*)=0.7623圖6收斂因子截圖(b)初值為(-0.4,0)時(shí)本算法令g的2范數(shù)在10-4時(shí)， 停止迭代， 經(jīng)過112次迭代收斂 收斂因子(f(k+1)-f*)/(f(k)-f*)=0.817圖7收斂因子截圖(c)初值為(10,0)時(shí)本算法令g的2范數(shù)在10-4時(shí)，停止迭代，經(jīng)過5次迭代收斂收斂因子(f(k+1)-f*)/(f(k)-f*)

8、=3.9022e-4圖8收斂因子截圖8(d)初值為(11,0)時(shí)本算法令g的2范數(shù)在10-4時(shí)，停止迭代，經(jīng)過2次迭代收斂 收斂因子(f(k+1)-f*)/(f(k)-f*)= 0圖9收斂因子截圖圖10自變量(x1,x2)截圖9總結(jié)： 最速降線法的收斂因子隨著初值的不同而變化， 對(duì)于個(gè)別初值 （如本習(xí)題初值 取 （11,0） 時(shí)） ，算法可迅速收斂。因此，初值的選取對(duì)于最速降線法的收斂速度有較大 影響。22習(xí)題5.7(a)由f(x) =9x -41 n(x-7)可得:故，牛頓迭代法的確切公式為:x -7(x-7)2（b）從以下五個(gè)初值開始迭代(1)x(0)=7.40表1初值1牛頓法迭代結(jié)果表迭

9、代次數(shù)x 值梯度值17.4-127.44-0.0909090937.4444-0.0009000947.4444444-9.00E-0857.4444444-9.00E-08(2) x(0)=7.20表2初值2牛頓法迭代結(jié)果表迭代次數(shù)x 值梯度值17.2-1127.31-3.637.403775-0.747.440723-7.60E-0257.444413-0.000631068(3)x(0)=7.01表3初值3牛頓法迭代結(jié)果表迭代次數(shù)x 值梯度值f (x) =94x -7f (x)二4(x -7)21017.01-3911127.019775-193.275600537.03867-94.4

10、7.073976-4.51E+0157.135638-20.(4)x(0)=7.80表4初值4牛頓法迭代結(jié)果表迭代次數(shù)x 值梯度值17.8427.16-1637.2624-6.947.369879-1.81E+0057.431934-0.3(5)x(0)=7.88表5初值5牛頓法迭代結(jié)果表迭代次數(shù)x 值梯度值17.884.527.0176-218.272727337.034503-106.931813547.066328-5.13E+0157.122757-23.（c）本問題的最優(yōu)值為7.4444444。由上述五個(gè)初值點(diǎn)的前五步迭代可以看出：當(dāng)初值點(diǎn)在區(qū)間（7.4444444，7.8888）內(nèi)

11、時(shí)，第二次迭代點(diǎn)將落在（7，7.4444444） 之間，隨后逐漸增加，直至逼近最優(yōu)值。當(dāng)初值點(diǎn)在區(qū)間（7，7.4444444）內(nèi)時(shí)，則迭代點(diǎn)逐漸增加，逼近最優(yōu)值。當(dāng)取初值不在（7,7.8888）內(nèi)時(shí)，牛頓法不收斂。2.3習(xí)題5.8（a）沒有線搜索的牛頓法卩=0.1時(shí)， 卩=1時(shí),（b）具有線搜索的牛頓法卩=0.1時(shí),卩=1時(shí)，（未完成）122.4習(xí)題5.9（a）初值選（1.2,1.2）時(shí)，最速降線法：本算法令g的2范數(shù)在10-2時(shí)，停止迭代，經(jīng)過3262次迭代得到以下結(jié)果Cent; ur ofX-awis圖11最速降線法初值為（121.2）的等值線圖及迭代軌跡牛頓法：本算法令s的4范數(shù)在10-

12、6時(shí)，停止迭代，經(jīng)過4次迭代得到以下結(jié)果13CQHlUr of圖12牛頓法初值為（1.2,1.2）的等值線圖及迭代軌跡（b）初值選（-1.2,1）時(shí)，最速降線法：本算法令g的4范數(shù)在10-2時(shí)，停止迭代，經(jīng)過6835次迭代得到以下結(jié)果Ccritouror Rosribr:ckX-axis14圖13最速降線法初值為（-1.2,1）的等值線圖及迭代軌跡15牛頓法：本算法令s的4范數(shù)在V10-6時(shí)，停止迭代，經(jīng)過6次迭代得到以下結(jié)果: nntnurbrockX- axis圖14牛頓法初值為（-1.2,1）的等值線圖及迭代軌跡2.5習(xí)題5.19N=5迭代6次后，滿足收斂條件。表6 N=5時(shí)，各迭代點(diǎn)x

13、值迭代次數(shù)/分量1444510000040.7744410.7744410.7744410.7744410.7744414-1.744581.0449944.4054814.8945444.45495444.740748-14.4459-4.780467.94544517.419445-4.8061445.656-86.4661-46.19499.441765.000468-140640000175-140640-1140640N=8迭代19次后，滿足收斂條件。表7 N=8時(shí)，各迭代點(diǎn)x值迭代次數(shù)/分量144456781610000000040.7544940.75

14、44940.7544940.7544940.7544940.7544940.7544940.7544944-1.714480.4868491.186971.7445094.1075684.4845444.5988544.77040844.619757-7.74774-5.18495-1.1844.6614656.0744949.04446111.64715175-4.5694948.44644-40.4749-44.1741-44.1857-1.644645.1444954.749864.577401-74.4411199.8408-9.99548-174.858-171.851-10.865

15、4469.47447-4.15904104.8461-645.041081.886446.6416-1014.4-914.4591401.4468-5.65084154.4884-874.6011478.11445.7144-1444.46-1156.481454.9449-5.64898154.4878-874.6011478.108445.7146-1444.47-1156.481454.944106.688844-489.1174810.06-9974.6514761.111879.748-15549.98487.119116.744416-489.4544810.494-9974.11

16、14764.741880.191-15541.78488.07146.789876-489.4814811.108-9976.4614765.091880.84-15544.48489.454147.645489-445.6814087.908-6008.644107.69416980.64-46494.611414.66147.405147-181.141454.897-1497.76-10051.444195.97-48554.514876.44155.044404-90.441479.065814858.444-47184.458454.08-55845.619754.9716-7.98

17、171504.7645-7559.4146199.76-148600416414.8-16816751480.4617-7.99416504.8604-7559.6446199.85-148600416415.4-16816851480.1618-7.99414504.8646-7559.6446199.85-148600416415.4-16816851480.1519-8.00048504.9999-756046400-148600416416-1681685148040-8.00004504.0004-756046400-148600416416-16816851480N=14迭代49次

18、后，滿足收斂條件。（表略）N=40迭代74次后，滿足收斂條件。（表略）2.6習(xí)題5.27調(diào)用MATLAB自帶的Isqnonlin.m函數(shù)，計(jì)算可得對(duì)應(yīng)的x（1）、x（2）和標(biāo)準(zhǔn)差如下表 所示。表8選取各初值的計(jì)算結(jié)果初值x1 虛實(shí)部x1 虛部x2 實(shí)部x2 虛部r標(biāo)準(zhǔn)差1 10.00070-15.33030.2767i0.00750.1 0.10.00060-14.41710.02670.00750.01 0.010.00060-12.94320.00260.0075由上可知，標(biāo)準(zhǔn)差值較為恒定，隨初值變化不十分顯著；x1和x2值隨初值選取的不同而不同。2.7習(xí)題6.4（未完成）end183附錄

19、3.1對(duì)偶單純形法函數(shù)MATLAB程序function sol,val,kk=duioudanchun(A,N)B=A;mA,nA=size(A);kk=0;flag=1;while flagkk=kk+1;if A(:,nA)=0flag=0;sol=zeros(1,nA);for i=1:mA-1 sol(N(i)=A(i,nA);endval=sol*(B(mA,:);elsefor i=1:mA-1if A(i,nA)=0disp( have infinite solution! ); flag=0;break ;endendif flagtemp=0;for i=1:mA-1if A

20、(i,nA)temp temp=A(i,nA); outb=i;endsita=zeros(1,nA-1);for i=1:nA-1if A(outb,i)0sita(i)=A(mA,i)/A(outb,i);endend19temp=-inf;for i=1:nA-1if sita(i)temptemp=sita(i);inb=i;endendfor i=1:mA-1if i=outb N(i)=inb;endendA(outb,:)=A(outb,:)/A(outb,inb);for i=1:mAif i=outbA(i,:)=A(i,:)-A(outb,:)*A(i,inb);A(mA,

21、nA)=0;endendendendend3.2最速降線法求Rosenbrock函數(shù)最小值matlab程序如下：function rb = rbfun(x,y)rb=100*(y-xA2)A2+(1-x)A2endclearclcsyms x y g Gg=gradient(rb(x,y),x y) %定義梯度向量G=hessian(rb(x,y),x y) %定義海森陣X(1,:)=-1.4 1;%定義初始點(diǎn)x=X(1,1);y=X(1,4);A(1,:)=subs(g)%給梯度賦初值20i=1while( norm(A(i,:),4)10A(-4)f(i)=rb(x,y)P(i,:)=-A

22、(i,:) fz(i)=-A(i,:)*P(i,:)%-gT*p fm(i)=P(i,:)*subs(G)*P(i,:)a(i)=fz(i)/fm(i)X(i+1,:) = X(i,:)+a(i)*P(i,:)x=X(i+1,1);y=X(i+1,4)A(i+1,:)=subs(g)i=i+1end%收斂條件%記錄函數(shù)值%得到迭代方向%精確搜索法步長(zhǎng)的分子%精確搜索法步長(zhǎng)的分母%精確搜索法步長(zhǎng)%產(chǎn)生新的點(diǎn)%產(chǎn)生新的梯度3.3牛頓法求Rosenbrock函數(shù)最小值matlab程序如下：function rb = rbfun(x,y)rb=100*(y-xA2)A2+(1-x)A2endclearclcsyms x y g Gg=gradient(rb(x,y),x y) %定義梯度向量21G=hessian(rb(x,y