最優(yōu)化方法課程設(shè)計-最優(yōu)化大作業(yè)-用優(yōu)化算法求解函數(shù)最值問題

1、最優(yōu)化方法大作業(yè)用優(yōu)化算法求解函數(shù)最值問題摘要最優(yōu)化（optimization）是應(yīng)用數(shù)學(xué)的重要研究領(lǐng)域它是研究在給定約束之下如何尋求某些因素（的量），以使某一（或某些）指標(biāo)達到最優(yōu)的一些學(xué)科的總稱。最優(yōu)化問題一般包括最小化問題和最大化問題，而最大化問題可以通過簡單的轉(zhuǎn)化使之成最最小化問題。最小化問題分為兩類，即約束最小化和無約束最小化問題。在此報告中，前兩個問題屬于無約束最小化問題的求解，報告中分別使用了“牛頓法”和“共軛梯度法”。后兩個問題屬于有約束最小化問題的求解，報告中分別用“外點法”和“內(nèi)點法”求解。雖然命名不一樣，其實質(zhì)都是構(gòu)造“懲罰函數(shù)”或者“障礙函數(shù)”，通過拉格朗日乘子法將有約

2、束問題轉(zhuǎn)化為無約束問題進行求解。再此報告中，“外點法”和“內(nèi)點法”分別用了直接求導(dǎo)和調(diào)用“牛頓法”來求解無約束優(yōu)化問題。在此實驗中，用“共軛梯度法”對“牛頓法”所解函數(shù)進行求解時出現(xiàn)錯誤，報告中另取一函數(shù)用“共軛梯度法”求解得到正確的結(jié)果。此實驗中所有的函數(shù)其理論值都是顯見的，分析計算結(jié)果可知程序正確，所求結(jié)果誤差處于可接受范圍內(nèi)。報告中對所用到的四種方法在其使用以前都有理論說明，對“外點法”中懲罰函數(shù)和“內(nèi)點法”中障礙函數(shù)的選擇也有相應(yīng)的說明，另外，對此次試驗中的收獲也在報告的三部分給出。本報告中所用程序代碼一律用MATLAB編寫?！娟P(guān)鍵字】函數(shù)最優(yōu)化牛頓法共軛梯度法內(nèi)點法外點法MATLAB

3、，問題描述1，分別用共軛梯度發(fā)法和牛頓法來求解一下優(yōu)化問題minf(x)=(x+10x+5(x1232，分別用外點法和內(nèi)點發(fā)求解一下優(yōu)化問題minx3+x212s.tx+x-1>012二、問題求解1.1用牛頓法求解min+10x)2+5(x1231.1.1問題分析：取步長為1而沿著牛頓方向迭代的方法稱為牛頓法，在牛頓法中，初始點的取值隨意，在以后的每次迭代中，x=x-t2f(x)LVf(x)，直到終止條件k+1kkk成立時停止。1.1.2問題求解注：本程序編程語言為MATLAB，終止條件為|VfCjl2<10-16，初始取值為10101010M文件(求解函數(shù))如下：function

4、s=newton1(f,c,eps)%c是初值，eps為允許誤差值ifnargin=2eps=1.0e-16;elseifnargin<1error('')%returnendsymsx1x2x3x4x=x1x2x3x4.'grad=jacobian(f).'hesse=jacobian(grad);a=grad;b=hesse;i=1;gradk=subs(a,x1x2x3x4,c(1)c(2)c(3)c(4);hessek=subs(b,x1x2x3x4,c(1)c(2)c(3)c(4);pk=-1*(hessekgradk);x=tihuan(c);

5、whilenorm(gradk)>=epsx=x+pk;gradk=subs(a,x1x2x3x4,x(1)x(2)x(3)x(4);hessek=subs(b,x1x2x3x4,x(1)x(2)x(3)x(4);pk=-hessekgradk;i=i+1;enddisp('thetimesofiterationis:')disp(i)disp('Thegradis:')disp(gradk)disp('andtheresultis:')x=x.'disp(x)return“tihuan”子函數(shù)：functionx=tihuan(x

6、)x(1)=x(1);x(2)=x(2);x(3)=x(3);x(4)=x(4);end調(diào)用方式如下：symsx1x2x3x4f=(xl+10*x2廠2+5*(x3-x4廠2+(x2-2*x3廠4+10*(xl-x4廠4;c=10101010'%初始值newton1(f,c,eps);1.1.3計算結(jié)果如下：>Tnlawtjn.：ilthetiiiLmufiterationis:47Thegrajdis:1.0e-016*-0.7031-0.DOOC0.DOOCa7031anthe"esu_tis：1.Oe-005*廣1由上述結(jié)果可知，當(dāng)?shù)螖?shù)達到47次時滿足終止條件

7、，此時x為1.0e-005*-0.11110.01110.00950.0095,顯然，此題的理論解為0000,分析上述結(jié)果，與理論解的誤差處于可接受范圍之內(nèi)。求解完成。1.2用共軛梯度法求解函數(shù)14minf(x)=(x+10x)2+5(x一x)2+(x一2x1+10(x一x112342用共軛梯度法求解上述函數(shù)的程序代碼如下：1.2.1問題分析：xk+1時，共軛方向取p=-Vf(x)，當(dāng)搜索到00pk+i=Vf(x)+Xp,k=0,1,，n2，此時,k+1kkpAp=Vf(x)Ap+XpTApk+1kk+1kkkkVf(x)TApp與pA共軛，用Ap右乘上式k+1kk，由PT+1Apk二0屮&q

8、uot;=0,1,.,n2，若不滿足條件，進行下一次迭代。pTAppk1.1.2問題求解注：程序所用語言為MATLAB，精度為eps二10-16symsx1x2x3x4t0tlf=(xl+10*x2廠2+5*(x3-x4廠2+(x2-2*x3廠4+10*(xl-x4廠4;c=10;10;10;10;grad1=diff(f,x1);grad2=diff(f,x2);grad3=diff(f,x3);grad4=diff(f,x4);grad=grad1;grad2;grad3;grad4;a=grad;i=1;n=40;gradk=subs(a,x1x2x3x4,c(1)c(2)c(3)c(4

9、);x=tihuan(c);p0=0;whilenorm(gradk)>=epsp0=-gradk;y=x;x=x+t0*p0;k=0;gradk=subs(a,x1x2x3x4,x(1)x(2)x(3)x(4);w=solve(gradk(1)+gradk(2)+gradk(3)+gradk(4);t0=real(w);t0=eval(t0);t0=t0(1);x=y+t0*p0;gradk=subs(a,x1x2x3x4,x(1)x(2)x(3)x(4);whilenorm(gradk)>=epsifk+1=ngradk2=subs(a,x1x2x3x4,x(1)x(2)x(3

10、)x(4);gradk1=subs(a,x1x2x3x4,y(1)y(2)y(3)y(4);lamda二norm(gradk2)."2/norm(gradkl)."2;p0=-gradk2+lamda*p0;k=k+1;elsek=0;p0=-subs(a,x1x2x3x4,x(1)x(2)x(3)x(4);endcleary;y=x;x=x+t1*p0;gradk=subs(f,x1x2x3x4,x(1)x(2)x(3)x(4);m=solve(gradk);t1=real(m);t1=eval(t1(1);x=x+t1*p0;x=eval(x);clearm;clear

11、t1;symst1gradk=subs(a,x1x2x3x4,x(1)x(2)x(3)x(4);enddisp(x.')return;enddisp(x.')此程序為一初步程序，假設(shè)初值為10;10;10;10，則第一次運算得t0=0.0011,lamda=0.9291,迭代后的x=NaN?，F(xiàn)用共軛梯度法求解另一函數(shù)minf(x)=x2+25x212對上述程序稍加改動來求解本題的代碼如下：注：程序所用語言為MATLAB，精度為eps二10-16functions=gongegrad2(f,c,eps)%c是初值，eps為允許誤差值ifnargin=2%eps=1.0e-16;e

12、lseifnargin<1error('')returnendticsymsx1x2t0t1grad1=diff(f,x1);grad2=diff(f,x2);grad=grad1;grad2;a=grad;i=1;n=40;gradk=subs(a,x1x2,c(1)c(2);x=tihuan(c);p0=0;whilenorm(gradk)>=epsp0=-gradk;x=x+t0*p0;k=0;gradk=subs(f,x1x2,x(1)x(2);w=solve(gradk);t0=real(w);t0=eval(t0);t0=t0(1);x=y+t0*p0;

13、gradk=subs(a,x1x2,x(1)x(2);whilenorm(gradk)>=epsifk+1=ngradk2=subs(a,x1x2,x(1)x(2);gradk1=subs(a,x1x2,y(1)y(2);lamda二norm(gradk2)"2/norm(gradkl)"2;p0=-gradk2+lamda*p0;k=k+1;elsek=0;p0=-subs(a,x1x2,x(1)x(2);endcleary;y=x;x=x+t1*p0;gradk=subs(f,x1x2,x(1)x(2);m=solve(gradk);t1=real(m);t1=e

14、val(t1(1);x=y+t1*p0;clearm;cleart1;symst1gradk=subs(a,x1x2,x(1)x(2);enddisp(sprintf('thelastpointwewantis%f%f',x(1),x(2);disp(sprintf('thetimesusedtorecursionis%f',k);disp(sprintf('thefunctionvalueis%f',x(l)"2+25*x(2)"2);tocreturn;enddisp(sprintf('thelastpointwe

15、wantis%f%f',x(1),x(2);disp(sprintf('thetimesusedtorecursionis%f',k);disp(sprintf('thefunctionvalueis%f',x(l)"2+25*x(2)"2);toctihuan”子函數(shù)為：functionx=tihuan(x)v,g=size(x);fori=1:vx(i)=x(i);end程序調(diào)用方式為：clearallclcsymsx1x2t0t1f=x2+25*x2"2;c=2;2;%初值gongegrad2(f,c,eps)程序結(jié)果

16、如下：CommandWindowthe1astpointwewantis0.0000000.000000thetimesusedtorecursionis3.000000thefunctionvalueis0.000000Elapsedtimeis0.273592seconds.由上述結(jié)果知，用共軛梯度法對上述函數(shù)求解需要迭代三次得到最優(yōu)解0此時x為00;符合理論分析的結(jié)果，求解完成。2.1用外點法法求解函數(shù)Immx3+x2112Is.tx+x-1>0122.1.1問題分析外點法為序列無約束最優(yōu)化方法，其基本思想為將條件函數(shù)吸收到目標(biāo)函數(shù)中進行求解。其在數(shù)學(xué)上的直觀理解是拉格朗日乘子法：

17、minTCx;M=min|fO+M£min(0,g(x)1,minTx;M為總代價，fG)為價格，M£Imin(0,g0)1為罰款。即在經(jīng)濟學(xué)上總代價為價格和罰款的和。ii=1此時tnin(0,g(x)1=|Q)i當(dāng)g(x)>0,當(dāng)g(jc0,(=1m)ii稱TX;M=f(x)+M£g+(x)為增廣目標(biāo)函數(shù)，通常取ii=1當(dāng)g(x)>0當(dāng)g(x)<0i2.1.2問題求解兩種方法求解程序如下：2.1.2.1注：程序所用語言為MATLAB，終止條件為-g(xR10，程序中i無約束優(yōu)化部分通過求導(dǎo)實現(xiàn)。M文件如下：ticclc%c是初值，eps為允許誤

18、差值ifnargin=1eps=1.0e-16;elseifnargin<1error('')returnendsymsmxl,x2=solve('3*x2+2*m*(xl+x2-l)=0','2*x2+2*m*(xl+x2-l)=0');t=1;k=1;x1=limit(x1,m,t);x2=limit(x2,m,t);bound=max(eval(x1(1)+x2(1)-1,eval(x1(2)+x2(2)-1);while-bound>epst=10*t;k=k+1;xl,x2=solve('3*x2+2*m*(xl+x

19、2-l)=0','2*x2+2*m*(xl+x2-l)=0');x1=limit(x1,m,t);x2=limit(x2,m,t);bound=max(eval(x1(1)+x2(1)-1,eval(x1(2)+x2(2)-1);endx1=eval(x1);x2=eval(x2);fl=xl(l)"3+x2(l)"2;f2=xl(2廠3+x2(2廠2;iff1<f2disp(sprintf('thefinalxis%f%f',x1(1),x2(1);disp(sprintf('thefinalfunctionvalue

20、is%f',f1);elsedisp(sprintf('thefinalxis%f%f',x1(2),x2(2);disp(sprintf('thefinalfunctionvalueis%f',f2);enddisp(sprintf('thetimesusedtorecursionis%f',k)toc調(diào)用方式如下symsx1x2mT=x3+x2"2+m*(xl+x2-l)"2;waidian(T,eps);functions=waidian(T,eps)實驗結(jié)果如下*口CommandWindowthefinalxi

21、sEC.5485840.451416tliefinalfuncticnvalusis0.368870the:imesusedtorecursionis17.000000E1apsedtimeis0.468610seconds.»由上述結(jié)果可知當(dāng)?shù)?7次時能夠達到終止條件，此時,x=0.5485841x=0.4514162函數(shù)得到最優(yōu)解0.368870.求解完成。2.1.2.2，程序中無約束優(yōu)化部分通過調(diào)用“牛頓法”完成代碼如下；ticclearall;closeall;clcsymsx1x2m=1;x0=2;2;eps=1.0e-6;T=x1入3+x2入2+m*(x1+x2-1)入

22、2;c=10;k=1;while1s=newton1(T,xO,eps);%調(diào)用newton法x1=s(1);x2=s(2);ifm*(x1+x2-1)入2epsm=c*m;k=k+1;symsx1x2T=x1入3+x2入2+m*(x1+x2-1)入2;elsedisp(sprintf('thefinalresultis%f%f',x1,x2);disp(sprintf('thefunctionvalueis%f',x1入3+x2入2);disp(sprintf('thetimesusedtorecursionis%f',k);break;end

23、endtoc實驗結(jié)果如下：CommandWindowthefinalresultis0.5485840.451416thefunctionvalueis0.368869thetimesusedtorecursionis7000000由上述結(jié)果可知當(dāng)?shù)?次時能夠達到終止條件，此時,Jx=0.548584|X=0.4514162函數(shù)得到最優(yōu)解0.368869.求解完成。2.2用內(nèi)點法法求解函數(shù)Iminx3+x2112Is.tx+x-1>0122.2.1問題分析同外點法，內(nèi)點法也為序列無約束最優(yōu)化方法，其基本思想為將條件函數(shù)吸收到目標(biāo)函數(shù)中進行求解。I(x;w)=f(x)+wB(x)C為邊界

24、時B(x')為+s,x不為邊界時B(x)tf(x)稱B(x)為障礙函數(shù)、圍墻函數(shù)或者懲罰函數(shù)。通常取B(x)=乏ln(g(x貶或者藝占i匸igi(x丿匸igi)2若在上一次迭代中x處于可解區(qū)域外部，則“圍墻”不夠高，在下一次迭代時加高“圍墻函數(shù)”。2.2.2問題求解求解程序如下：注：程序所用語言為MATLAB，本程序出于對運行時間和算法簡單的考慮，選用B(x)=蘭一)i=1'i終止條件為g(x106，此時可以認(rèn)為當(dāng)x接近求解邊界時B(x)(圍墻)為i無窮大，以確保函數(shù)的自變量迭代發(fā)生在可解區(qū)域內(nèi)部2.2.2.1，程序中無約束優(yōu)化部分通過求導(dǎo)實現(xiàn)。ticclearall;clos

25、eall;clc;eps=10飛；symswx1x2I=x3+x2'2+w*(l/(xl+x2-l);x10=2;x20=2;xl,x2=solve('3*x2-w/(xl+x2-l)"2=0','2*x2-w/(xl+x2-l)"2=0');c=10;k=1;x1=limit(x1,w,c);x2=limit(x2,w,c);x1=eval(x1);x2=eval(x2);x1=real(x1(1);x2=real(x2(1);bound=(1/(x1+x2-1);whilebound<=epsc=c/10;k=k+1;xl,

26、x2=solve('3*x2-w/(xl+x2-l)"2=0','2*x2-w/(xl+x2-l)"2=0');x1=limit(x1,w,c);x2=limit(x2,w,c);x1=eval(x1);x2=eval(x2);x1=real(x1(1);x2=real(x2(1);bound=(1/(x1+x2-1);enddisp(sprintf('thefinalxis%f%f',x1,x2);disp(sprintf('thefinalfunctionvalueis%f',xl(l)"3+x2

27、(l)"2);disp(sprintf('thetimesusedtorecursionis%f',k)toc實驗結(jié)果如下：CommandWindowthefinalxis0.5485840.451416thefinalfunctio門valueis0.368870thetimesusedtorecursionis15.000000Elapsedtimeis2.150889seconds.»由上述結(jié)果可知當(dāng)?shù)?5次時能夠達到終止條件，此時,Jx=0.548584|X=0.4514162函數(shù)得到最優(yōu)解0.368870.求解完成，結(jié)果與“外點法”相同。2.2.2.2，程序中無約束優(yōu)化部分通過調(diào)用“牛頓法”完成(c=,收斂準(zhǔn)則為|xk-xk<eps)代碼如下；ticclearall;closeall;clcsymsx1x2w=10;x0=2;2;eps=1.0e-6;I=x1入3+x2入2+w*(1/(x1+x2-1)入2);c=0.3;k=1;while1s=newton1(I,x0,eps);%調(diào)用newton法x1=s(1);x2=s(2);ifsqrt(x1-x0(1)入2+(x2-x0(2)入2)epsw=c*w;k=k+1;x0(1)=x1;x0(2)=x2;