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

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

3、點(diǎn)法 外點(diǎn)法 MATLAB一，問(wèn)題描述1， 分別用共軛梯度發(fā)法和牛頓法來(lái)求解一下優(yōu)化問(wèn)題2， 分別用外點(diǎn)法和內(nèi)點(diǎn)發(fā)求解一下優(yōu)化問(wèn)題二、問(wèn)題求解1.1 用牛頓法求解 1.1.1問(wèn)題分析：取步長(zhǎng)為1而沿著牛頓方向迭代的方法稱(chēng)為牛頓法，在牛頓法中，初始點(diǎn)的取值隨意，在以后的每次迭代中，直到終止條件成立時(shí)停止。1.1.2 問(wèn)題求解注：本程序編程語(yǔ)言為MATLAB，終止條件為，初始取值為10 10 10 10M文件（求解函數(shù)）如下：function s=newton1(f,c,eps)%c 是初值，eps為允許誤差值 if nargin=2 eps=1.0e-16;elseif nargin<1

4、error('') % returnendsyms x1 x2 x3 x4x=x1 x2 x3 x4.'grad = jacobian(f).'hesse = jacobian(grad);a=grad;b=hesse;i=1;gradk=subs(a,x1 x2 x3 x4,c(1) c(2) c(3) c(4);hessek=subs(b,x1 x2 x3 x4,c(1) c(2) c(3) c(4); pk=-1*(hessekgradk); x=tihuan(c);while norm(gradk)>=eps x=x+pk; gradk=subs(

5、a,x1 x2 x3 x4,x(1) x(2) x(3) x(4); hessek=subs(b,x1 x2 x3 x4,x(1) x(2) x(3) x(4); pk=-hessekgradk; i=i+1;enddisp('the times of iteration is:')disp(i)disp('The grad is:')disp(gradk)disp('and the result is:')x=x.'disp(x)return“tihuan”子函數(shù):function x=tihuan(x)x(1)=x(1);x(2)=x

6、(2);x(3)=x(3);x(4)=x(4);end調(diào)用方式如下：syms x1 x2 x3 x4f=(x1+10*x2)2+5*(x3-x4)2+(x2-2*x3)4+10*(x1-x4)4;c=10 10 10 10'%初始值newton1(f,c,eps);1.1.3 計(jì)算結(jié)果如下： 由上述結(jié)果可知，當(dāng)?shù)螖?shù)達(dá)到47次時(shí)滿(mǎn)足終止條件，此時(shí)x為1.0e-005 * -0.1111 0.0111 0.0095 0.0095,顯然，此題的理論解為0 0 0 0,分析上述結(jié)果，與理論解的誤差處于可接受范圍之內(nèi)。求解完成。1.2 用共軛梯度法求解函數(shù) 用共軛梯度法求解上述函數(shù)的程序代碼

7、如下：1.2.1問(wèn)題分析：取，當(dāng)搜索到時(shí)，共軛方向，此時(shí)，與A共軛，用右乘上式得，由得，若不滿(mǎn)足條件，進(jìn)行下一次迭代。1.1.2 問(wèn)題求解注：程序所用語(yǔ)言為MATLAB，精度為syms x1 x2 x3 x4 t0 t1f=(x1+10*x2)2+5*(x3-x4)2+(x2-2*x3)4+10*(x1-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(

8、a,x1 x2 x3 x4,c(1) c(2) c(3) c(4); x=tihuan(c); p0=0;while norm(gradk)>=eps p0=-gradk; y=x; x=x+t0*p0; k=0; gradk=subs(a,x1 x2 x3 x4,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,x1 x2 x3 x4,x(1) x(2) x(3) x(4); while

9、norm(gradk)>=eps if k+1=n gradk2=subs(a,x1 x2 x3 x4,x(1) x(2) x(3) x(4); gradk1=subs(a,x1 x2 x3 x4,y(1) y(2) y(3) y(4); lamda=norm(gradk2).2/norm(gradk1).2; p0=-gradk2+lamda*p0; k=k+1; else k=0; p0=-subs(a,x1 x2 x3 x4,x(1) x(2) x(3) x(4); end clear y; y=x; x=x+t1*p0; gradk=subs(f,x1 x2 x3 x4,x(1)

10、 x(2) x(3) x(4); m=solve(gradk); t1=real(m); t1=eval(t1(1); x=x+t1*p0;x=eval(x);clear m;clear t1;syms t1 gradk=subs(a,x1 x2 x3 x4,x(1) x(2) x(3) x(4); end disp(x.') return;end disp(x.')此程序?yàn)橐怀醪匠绦?，假設(shè)初值為10;10;10;10，則第一次運(yùn)算得t0=0.0011,lamda=0.9291,迭代后的x=NaN?，F(xiàn)用共軛梯度法求解另一函數(shù) 對(duì)上述程序稍加改動(dòng)來(lái)求解本題的代碼如下：注：程序所用

11、語(yǔ)言為MATLAB，精度為function s=gongegrad2(f,c,eps)%c 是初值，eps為允許誤差值 if nargin=2 %eps=1.0e-16;elseif nargin<1 error('') returnendticsyms x1 x2 t0 t1grad1 = diff(f,x1);grad2=diff(f,x2);grad=grad1;grad2;a=grad;i=1;n=40;gradk=subs(a,x1 x2,c(1) c(2); x=tihuan(c); p0=0;while norm(gradk)>=eps p0=-gra

12、dk; y=x; x=x+t0*p0; k=0; gradk=subs(f,x1 x2,x(1) x(2); w=solve(gradk); t0=real(w); t0=eval(t0); t0=t0(1); x=y+t0*p0; gradk=subs(a,x1 x2,x(1) x(2); while norm(gradk)>=eps if k+1=n gradk2=subs(a,x1 x2,x(1) x(2); gradk1=subs(a,x1 x2,y(1) y(2); lamda=norm(gradk2)2/norm(gradk1)2; p0=-gradk2+lamda*p0;

13、k=k+1; else k=0; p0=-subs(a,x1 x2,x(1) x(2); end clear y; y=x; x=x+t1*p0; gradk=subs(f,x1 x2,x(1) x(2); m=solve(gradk); t1=real(m); t1=eval(t1(1); x=y+t1*p0; clear m;clear t1;syms t1 gradk=subs(a,x1 x2,x(1) x(2); end disp(sprintf('the last point we want is %f %f',x(1),x(2); disp(sprintf('

14、;the times used to recursion is %f',k); disp(sprintf('the function value is %f',x(1)2+25*x(2)2); toc return;end disp(sprintf('the last point we want is %f %f',x(1),x(2); disp(sprintf('the times used to recursion is %f',k); disp(sprintf('the function value is %f',x

15、(1)2+25*x(2)2); toc“tihuan”子函數(shù)為：function x=tihuan(x)v,g=size(x);for i=1:vx(i)=x(i);end程序調(diào)用方式為：clear allclcsyms x1 x2 t0 t1f=x12+25*x22;c=2;2;%初值gongegrad2(f,c,eps)程序結(jié)果如下：由上述結(jié)果知，用共軛梯度法對(duì)上述函數(shù)求解需要迭代三次得到最優(yōu)解0，此時(shí)x為0 0;符合理論分析的結(jié)果，求解完成。2.1 用外點(diǎn)法法求解函數(shù) 2.1.1 問(wèn)題分析 外點(diǎn)法為序列無(wú)約束最優(yōu)化方法，其基本思想為將條件函數(shù)吸收到目標(biāo)函數(shù)中進(jìn)行求解。其在數(shù)學(xué)上的直觀(guān)理解

16、是拉格朗日乘子法：，為總代價(jià)，為價(jià)格，為罰款。即在經(jīng)濟(jì)學(xué)上總代價(jià)為價(jià)格和罰款的和。此時(shí) 稱(chēng)為增廣目標(biāo)函數(shù)，通常取 2.1.2 問(wèn)題求解 兩種方法求解程序如下： 2.1.2.1 注：程序所用語(yǔ)言為MATLAB，終止條件為，程序中無(wú)約束優(yōu)化部分通過(guò)求導(dǎo)實(shí)現(xiàn)。 M文件如下： ticclc%c 是初值，eps為允許誤差值 if nargin=1 eps=1.0e-16;elseif nargin<1 error('') returnendsyms mx1,x2=solve('3*x12+2*m*(x1+x2-1)=0','2*x2+2*m*(x1+x2-1

17、)=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>eps t=10*t;k=k+1; x1,x2=solve('3*x12+2*m*(x1+x2-1)=0','2*x2+2*m*(x1+x2-1)=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);end

18、x1=eval(x1);x2=eval(x2);f1=x1(1)3+x2(1)2;f2=x1(2)3+x2(2)2;if f1<f2disp(sprintf('the final x is %f %f',x1(1),x2(1);disp(sprintf('the final function value is %f',f1);else disp(sprintf('the final x is %f %f',x1(2),x2(2);disp(sprintf('the final function value is %f',f2

19、); enddisp(sprintf('the times used to recursion is %f',k)toc 調(diào)用方式如下 syms x1 x2 mT=x13+x22+m*(x1+x2-1)2;waidian(T,eps);function s=waidian(T,eps)實(shí)驗(yàn)結(jié)果如下由上述結(jié)果可知當(dāng)?shù)?7次時(shí)能夠達(dá)到終止條件，此時(shí), 函數(shù)得到最優(yōu)解0.368870.求解完成。2.1.2.2，程序中無(wú)約束優(yōu)化部分通過(guò)調(diào)用“牛頓法”完成代碼如下；ticclear all;close all;clcsyms x1 x2m=1;x0=2;2;eps=1.0e-6;T=x

20、13+x22+m*(x1+x2-1)2;c=10;k=1;while 1s=newton1(T,x0,eps);%調(diào)用newton法x1=s(1);x2=s(2);if m*(x1+x2-1)2>eps m=c*m;k=k+1;syms x1 x2 T=x13+x22+m*(x1+x2-1)2;else disp(sprintf('the final result is %f %f',x1,x2); disp(sprintf('the function value is %f',x13+x22); disp(sprintf('the times u

21、sed to recursion is %f',k); break;endendtoc實(shí)驗(yàn)結(jié)果如下： 由上述結(jié)果可知當(dāng)?shù)?次時(shí)能夠達(dá)到終止條件，此時(shí), 函數(shù)得到最優(yōu)解0.368869.求解完成。2.2 用內(nèi)點(diǎn)法法求解函數(shù) 2.2.1 問(wèn)題分析 同外點(diǎn)法，內(nèi)點(diǎn)法也為序列無(wú)約束最優(yōu)化方法，其基本思想為將條件函數(shù)吸收到目標(biāo)函數(shù)中進(jìn)行求解。 稱(chēng)為障礙函數(shù)、圍墻函數(shù)或者懲罰函數(shù)。通常取 若在上一次迭代中x處于可解區(qū)域外部，則“圍墻”不夠高，在下一次迭代時(shí)加高“圍墻函數(shù)”。2.2.2 問(wèn)題求解 求解程序如下： 注：程序所用語(yǔ)言為MATLAB，本程序出于對(duì)運(yùn)行時(shí)間和算法簡(jiǎn)單的考慮，選用終止條件為，

22、此時(shí)可以認(rèn)為當(dāng)x接近求解邊界時(shí)（圍墻）為無(wú)窮大，以確保函數(shù)的自變量迭代發(fā)生在可解區(qū)域內(nèi)部2.2.2.1，程序中無(wú)約束優(yōu)化部分通過(guò)求導(dǎo)實(shí)現(xiàn)。 ticclear all;close all;clc;eps=106;syms w x1 x2I=x13+x22+w*(1/(x1+x2-1);x10=2;x20=2;x1,x2=solve('3*x12-w/(x1+x2-1)2=0','2*x2-w/(x1+x2-1)2=0');c=10;k=1;x1=limit(x1,w,c);x2=limit(x2,w,c);x1=eval(x1);x2=eval(x2);x1=re

23、al(x1(1);x2=real(x2(1);bound=(1/(x1+x2-1);while bound<=eps c=c/10;k=k+1;x1,x2=solve('3*x12-w/(x1+x2-1)2=0','2*x2-w/(x1+x2-1)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('the final x is %f %f

24、',x1,x2);disp(sprintf('the final function value is %f',x1(1)3+x2(1)2);disp(sprintf('the times used to recursion is %f',k)toc 實(shí)驗(yàn)結(jié)果如下： 由上述結(jié)果可知當(dāng)?shù)?5次時(shí)能夠達(dá)到終止條件，此時(shí), 函數(shù)得到最優(yōu)解0.368870.求解完成，結(jié)果與“外點(diǎn)法”相同。2.2.2.2，程序中無(wú)約束優(yōu)化部分通過(guò)調(diào)用“牛頓法”完成（）代碼如下；ticclear all;close all;clcsyms x1 x2w=10;x0=2;2;eps=1.0e-6;I=x13+x22+w*(1/(x1+x2-1)2);c=0.3;k=1;while 1s=newton1(I,x0,eps);%調(diào)用newton法x1=s(1);x2=s(2);if sqrt(x1-x0(1)2+(x2-x0(2)2)>eps w=c*w