最優(yōu)化課程設(shè)計(jì)--共軛梯度法算法分析與實(shí)現(xiàn)

1、<p><b>　　最優(yōu)化方法課程設(shè)計(jì)</b></p><p><b>　?。ㄔO(shè)計(jì)程序）</b></p><p>　　題目 共軛梯度法算法分析與實(shí)現(xiàn) </p><p>　　班級(jí) / 學(xué)號(hào) 14140101/2011041401011 </p><p>　　學(xué) 生 姓 名

2、 黃中武 </p><p>　　指 導(dǎo) 教 師 王吉波 王微微 </p><p>　　課 程 設(shè) 計(jì) 任 務(wù) 書</p><p>　　課 程 名 稱 最優(yōu)化方法課程設(shè)計(jì) </p><p>　　院（系） 理學(xué)院 專業(yè) 信息與計(jì)算科學(xué)

3、 </p><p>　　課程設(shè)計(jì)題目 共軛梯度法算法分析與實(shí)現(xiàn) </p><p>　　課程設(shè)計(jì)時(shí)間: 2014 年 6月 16日至 2014 年 6月 27日</p><p>　　課程設(shè)計(jì)的要求及內(nèi)容：</p><p><b>　　[要求]</b></p><p>　　學(xué)習(xí)態(tài)度要

4、認(rèn)真，要積極參與課程設(shè)計(jì)，鍛煉獨(dú)立思考能力；</p><p>　　嚴(yán)格遵守上機(jī)時(shí)間安排；</p><p>　　按照MATLAB編程訓(xùn)練的任務(wù)要求來(lái)編寫程序；</p><p>　　根據(jù)任務(wù)書來(lái)完成課程設(shè)計(jì)論文；</p><p>　　報(bào)告書寫格式要求按照沈陽(yáng)航空航天大學(xué)“課程設(shè)計(jì)報(bào)告撰寫規(guī)范”；</p><p>　　報(bào)告上

5、交時(shí)間：課程設(shè)計(jì)結(jié)束時(shí)上交報(bào)告；</p><p>　　嚴(yán)禁抄襲行為，一旦發(fā)現(xiàn)，課程設(shè)計(jì)成績(jī)?yōu)椴患案瘛?lt;/p><p>　　運(yùn)用共軛梯度法求解無(wú)約束最優(yōu)化問題</p><p>　　要求：1）了解求解無(wú)約束最優(yōu)化問題的共軛梯度法；</p><p>　　2）繪出程序流程圖；</p><p>　　3）編寫求解無(wú)約束最優(yōu)化問題的

6、共軛梯度法MATLAB程序；</p><p>　　4）利用編寫文件求解某無(wú)約束最優(yōu)化問題；</p><p><b>　　5）給出程序注釋。</b></p><p>　　指導(dǎo)教師 年 月 日</p><p>　　負(fù)責(zé)教師

7、 年 月 日</p><p>　　學(xué)生簽字 年 月 日</p><p><b>　　沈陽(yáng)航空航天大學(xué)</b></p><p>　　課 程 設(shè) 計(jì) 成 績(jī) 評(píng) 定 單</p><p>　　課 程 名 稱 最優(yōu)化理論與算法課程設(shè)計(jì)

8、 </p><p>　　院（系） 理學(xué)院 專業(yè) 信息與計(jì)算科學(xué) </p><p>　　課程設(shè)計(jì)題目 共軛梯度法算法分析與實(shí)現(xiàn) </p><p>　　學(xué)號(hào) 2011041401011 姓名 黃中武 </p><p>&l

9、t;b>　　指導(dǎo)教師評(píng)語(yǔ)：</b></p><p>　　課程設(shè)計(jì)成績(jī) </p><p>　　指導(dǎo)教師簽字 </p><p><b>　　年 月 日</b></p><p><b>　　目

10、錄</b></p><p><b>　　一、正文1</b></p><p><b>　　二、總結(jié)2</b></p><p><b>　　參考文獻(xiàn)3</b></p><p><b>　　附 錄4</b></p><p&

11、gt;<b>　　一、正文</b></p><p>　　一 無(wú)約束最優(yōu)化問題的共軛梯度法</p><p>　　共軛梯度法最初是由Hesteness和Stiefel于1952年為求解線形方程組而提出的。后來(lái)，人們把這種方法用于求解無(wú)約束最優(yōu)化問題，使之成為一種重要的最優(yōu)化方法。</p><p>　　下面，重點(diǎn)介紹Fletcher-Reeves共軛

12、梯度法，簡(jiǎn)稱FR法。</p><p>　　共軛梯度法的基本思想是把共軛性與最速下降法相結(jié)合，利用已知點(diǎn)處的梯度構(gòu)造一組共軛方向，并沿這組方向進(jìn)行搜索，求出目標(biāo)函數(shù)的極小點(diǎn)。根據(jù)共軛梯度方向的基本性質(zhì)，這種方法具有二次終止性。首先討論對(duì)于二次函數(shù)的共軛梯度法，然后再把這種方法推廣到極小化一般函數(shù)的情形?？紤]問題</p><p>　　其中A是對(duì)稱正定矩陣，c是常數(shù)。</p>&l

13、t;p><b>　　具體求解方法如下：</b></p><p>　　首先，任意給定一個(gè)初始點(diǎn)，計(jì)算出目變函數(shù)在這點(diǎn)的梯度，若=0，則停止計(jì)算，否則，令</p><p>　　沿方向搜索，得到點(diǎn)，計(jì)算在處的梯度，若，則利用和構(gòu)造第二個(gè)搜索方向，再沿搜索。</p><p>　　一般地，若已知點(diǎn)和搜索方向，則從出發(fā)，沿進(jìn)行搜索，得到</p&

14、gt;<p><b>　　其中步長(zhǎng)滿足</b></p><p>　　此時(shí)可求出的顯式表達(dá)。令</p><p><b>　　求得極小點(diǎn)，令</b></p><p>　　根據(jù)二次函數(shù)的梯度表達(dá)式，即</p><p><b>　?。?0.3.16）</b></p&

15、gt;<p>　　由（10.3.16）式得到</p><p><b>　　（10.3.17）</b></p><p>　　計(jì)算在處的梯度，若，則停止計(jì)算；否則，用和構(gòu)造下一個(gè)搜索方向，并使和關(guān)于共軛。按此設(shè)想，令 </p><p><b>　　（10.3.18）</b></p><p&g

16、t;<b>　　上式兩端左乘，并令</b></p><p><b>　　由此得到</b></p><p><b>　　（10.3.19）</b></p><p>　　再?gòu)某霭l(fā)，沿方向搜索</p><p>　　綜上分析，在第一個(gè)搜索方向取負(fù)梯度的前提下，重復(fù)使用公式（10.3.1

17、4），（10.3.17），（10.3.18）和（10.3.19），就能伴隨計(jì)算點(diǎn)的增加，構(gòu)造出一組搜索方向。</p><p><b>　　二 程序流程圖</b></p><p><b>　　N</b></p><p><b>　　Y </b></p><p>　　圖一 共軛梯度

18、法程序流程圖</p><p>　　三 共軛梯度法的MATLAB程序</p><p>　　function [x,val,k]=frcg(fun,gfun,x0)</p><p>　　maxk=5000;</p><p>　　rho=0.6;sigma=0.4;</p><p>　　k=0;epsilon=1e-7;&l

19、t;/p><p>　　n=length(x0);</p><p>　　while (k<maxk)</p><p>　　g=feval(gfun,x0);</p><p>　　itern=k-(n+1)*floor(k/(n+1));</p><p>　　itern=itern+1;</p><p

20、>　　if(itern==1)</p><p><b>　　d=-g;</b></p><p><b>　　else</b></p><p>　　beta=(g'*g)/(g0'*g0);</p><p>　　d=-g+beta*d0;</p><p>

21、<b>　　gd=g'*d;</b></p><p><b>　　if(gd>=0)</b></p><p><b>　　d=-g;</b></p><p><b>　　end</b></p><p><b>　　end</b

22、></p><p>　　if (norm(g)<epsilon)</p><p><b>　　break;</b></p><p><b>　　end</b></p><p><b>　　m=0;</b></p><p><b>　

23、　mk=0;</b></p><p>　　while (m<20)</p><p>　　if (feval(fun,x0+rho^m*d)<feval(fun,x0)+sigma*rho^m*g'*d)</p><p><b>　　mk=m;</b></p><p><b>　　

24、break;</b></p><p><b>　　end</b></p><p><b>　　m=m+1;</b></p><p><b>　　end</b></p><p>　　x0=x0+rho^mk*d;</p><p>　　val=f

25、eval(fun,x0);</p><p><b>　　g0=g;</b></p><p><b>　　d0=d;</b></p><p><b>　　k=k+1;</b></p><p><b>　　end</b></p><p>

26、;<b>　　x=x0;</b></p><p>　　val=feval(fun,x);</p><p>　　function f=fun(x)</p><p>　　f=100*(x(1)^2-x(2))^2+(x(1)-1)^2;</p><p>　　function g=gfun(x)</p><

27、p>　　g=[400*x(1)*(x(1)^2-x(2))+2*(x(1)-1),-200*(x(1)^2-x(2))]'</p><p>　　四 利用所編程序求解實(shí)際問題</p><p>　　例 </p><p><b>　　該問題有精確解</b></p><p>　　利用所編程序求解此問

28、題，在命令窗口輸入：</p><p><b>　　運(yùn)行結(jié)果為：</b></p><p><b>　　五 程序注釋</b></p><p>　　function [x,val,k]=frcg(fun,gfun,x0)</p><p>　　%用FR共軛梯度法求解無(wú)約束問題：min f(x)</p&

29、gt;<p>　　%輸入：x0時(shí)初始點(diǎn)，fun，gfun分別是目標(biāo)函數(shù)和梯度</p><p>　　%輸出：x，val分別是近似最優(yōu)點(diǎn)和最優(yōu)值，k是迭代次數(shù)。</p><p>　　maxk=5000;%最大迭代次數(shù)</p><p>　　rho=0.6;sigma=0.4;</p><p>　　k=0;epsilon=1e-7;&l

30、t;/p><p>　　n=length(x0);</p><p>　　while (k<maxk)</p><p>　　g=feval(gfun,x0);%計(jì)算梯度</p><p>　　itern=k-(n+1)*floor(k/(n+1));</p><p>　　itern=itern+1;%計(jì)算搜索方向</

31、p><p>　　if(itern==1)</p><p><b>　　d=-g;</b></p><p><b>　　else</b></p><p>　　beta=(g'*g)/(g0'*g0);</p><p>　　d=-g+beta*d0;</p>

32、;<p><b>　　gd=g'*d;</b></p><p><b>　　if(gd>=0)</b></p><p><b>　　d=-g;</b></p><p><b>　　end</b></p><p><b>

33、;　　end</b></p><p>　　if (norm(g)<epsilon)%檢驗(yàn)終止條件</p><p><b>　　break;</b></p><p><b>　　end</b></p><p><b>　　m=0;</b></p>

34、<p><b>　　mk=0;</b></p><p>　　while (m<20)</p><p>　　if (feval(fun,x0+rho^m*d)<feval(fun,x0)+sigma*rho^m*g'*d)</p><p><b>　　mk=m;</b></p>&

35、lt;p><b>　　break;</b></p><p><b>　　end</b></p><p><b>　　m=m+1;</b></p><p><b>　　end</b></p><p>　　x0=x0+rho^mk*d;</p>

36、;<p>　　val=feval(fun,x0);</p><p><b>　　g0=g;</b></p><p><b>　　d0=d;</b></p><p><b>　　k=k+1;</b></p><p><b>　　end</b>&

37、lt;/p><p><b>　　x=x0;</b></p><p>　　val=feval(fun,x);</p><p>　　function f=fun(x)</p><p>　　f=100*(x(1)^2-x(2))^2+(x(1)-1)^2;%需要求解的函數(shù)</p><p>　　functio

38、n g=gfun(x)</p><p>　　g=[400*x(1)*(x(1)^2-x(2))+2*(x(1)-1),-200*(x(1)^2-x(2))]'%梯度</p><p><b>　　二、總結(jié)</b></p><p>　　不知不覺，一個(gè)學(xué)期就這樣過(guò)去了，隨著期末尾聲的到來(lái)，《最優(yōu)化》這門課也已經(jīng)結(jié)課了，接著而來(lái)的是為期兩個(gè)星期

39、的課程設(shè)計(jì)。</p><p>　　作為信息與計(jì)算科學(xué)系的一員，我們?cè)谡n堂上接觸最多的就是理論部分，剛好課程設(shè)計(jì)可以給我們驗(yàn)證理論，聯(lián)系實(shí)際的機(jī)會(huì)，所以，這次課程設(shè)計(jì)，我一直都在認(rèn)真完成。在拿到自己的課程設(shè)計(jì)題目的時(shí)候，我對(duì)它進(jìn)行了一定時(shí)間的審題，題目中要求我用共軛梯度法求解無(wú)約束優(yōu)化問題， 所以我拿到課本，仔細(xì)地看了共軛梯度法的原理與例題求解，以便幫助我更好的理解共軛梯度法的應(yīng)用。</p><

40、p>　　我記得，再一次研究生面試中，有的老師問過(guò)考生：你在你的大學(xué)生涯中學(xué)到了什么？如果讓我回答，我會(huì)說(shuō)：我學(xué)到了一種發(fā)現(xiàn)問題與解決問題的方法，這種方法是在我的學(xué)習(xí)中慢慢形成的一種潛移默化的思想，關(guān)于學(xué)到的內(nèi)容，我們不會(huì)一直保留太清晰的記憶，也不會(huì)在以后的生活中用到太多，就像買菜不會(huì)用到微積分，不會(huì)求導(dǎo)，但是，數(shù)學(xué)方法和思想給我們思考問題所帶來(lái)的方法論卻是大有裨益。</p><p>　　這次課程設(shè)計(jì)，我發(fā)現(xiàn)

41、我的很多不足，比如在MATLAB上的應(yīng)用上，我還不熟練，這與我在平時(shí)的練習(xí)中沒有加以重視有關(guān)，而且在《最優(yōu)化》的很多知識(shí)點(diǎn)上，自己還沒有很好地掌握。</p><p>　　最后 ，我要感謝我們的任課教師王吉波老師，感謝他在給我們上課中的帶來(lái)的幽默感和知識(shí)的傳授。在此，祝愿老師身體健康，工作順利。</p><p>　　完稿日期: 2014 年 6月 21日</p><p

42、><b>　　參考文獻(xiàn)</b></p><p>　　[1] 陳寶林.《最優(yōu)化理論與算法》[M].北京,清華大學(xué)出版社,2013.</p><p>　　[2] 劉衛(wèi)國(guó).《MATLAB程序設(shè)計(jì)教程》[M].北京,中國(guó)水利水電出版社,2010.</p><p>　　[3] 馬昌鳳.《最優(yōu)化方法及其MATLAB程序設(shè)計(jì)》[M].北京,科學(xué)出版社,

43、2010.</p><p><b>　　附 錄</b></p><p><b>　　源程序:</b></p><p>　　function [x,val,k]=frcg(fun,gfun,x0)</p><p>　　maxk=5000;</p><p>　　rho=0.6;s

44、igma=0.4;</p><p>　　k=0;epsilon=1e-7;</p><p>　　n=length(x0);</p><p>　　while (k<maxk)</p><p>　　g=feval(gfun,x0);</p><p>　　itern=k-(n+1)*floor(k/(n+1));<

45、;/p><p>　　itern=itern+1;</p><p>　　if(itern==1)</p><p><b>　　d=-g;</b></p><p><b>　　else</b></p><p>　　beta=(g'*g)/(g0'*g0);</p

46、><p>　　d=-g+beta*d0;</p><p><b>　　gd=g'*d;</b></p><p><b>　　if(gd>=0)</b></p><p><b>　　d=-g;</b></p><p><b>　　end

47、</b></p><p><b>　　end</b></p><p>　　if (norm(g)<epsilon)</p><p><b>　　break;</b></p><p><b>　　end</b></p><p><b

48、>　　m=0;</b></p><p><b>　　mk=0;</b></p><p>　　while (m<20)</p><p>　　if (feval(fun,x0+rho^m*d)<feval(fun,x0)+sigma*rho^m*g'*d)</p><p><b&g

49、t;　　mk=m;</b></p><p><b>　　break;</b></p><p><b>　　end</b></p><p><b>　　m=m+1;</b></p><p><b>　　end</b></p><

50、;p>　　x0=x0+rho^mk*d;</p><p>　　val=feval(fun,x0);</p><p><b>　　g0=g;</b></p><p><b>　　d0=d;</b></p><p><b>　　k=k+1;</b></p>&l

51、t;p><b>　　end</b></p><p><b>　　x=x0;</b></p><p>　　val=feval(fun,x);</p><p>　　function f=fun(x)</p><p>　　f=100*(x(1)^2-x(2))^2+(x(1)-1)^2;</p