最優(yōu)化大作業(yè)

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

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

3、點法 外點法 MATLAB一，問題描述1， 分別用共軛梯度發(fā)法和牛頓法來求解一下優(yōu)化問題2， 分別用外點法和內(nèi)點發(fā)求解一下優(yōu)化問題二、問題求解1.1 用牛頓法求解 1.1.1問題分析：取步長為1而沿著牛頓方向迭代的方法稱為牛頓法，在牛頓法中，初始點的取值隨意，在以后的每次迭代中，直到終止條件成立時停止。1.1.2 問題求解注：本程序編程語言為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 計算結果如下： 由上述結果可知，當?shù)螖?shù)達到47次時滿足終止條件，此時x為1.0e-005 * -0.1111 0.0111 0.0095 0.0095,顯然，此題的理論解為0 0 0 0,分析上述結果，與理論解的誤差處于可接受范圍之內(nèi)。求解完成。1.2 用共軛梯度法求解函數(shù) 用共軛梯度法求解上述函數(shù)的程序代碼

7、如下：1.2.1問題分析：取，當搜索到時，共軛方向，此時，與A共軛，用右乘上式得，由得，若不滿足條件，進行下一次迭代。1.1.2 問題求解注：程序所用語言為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.')此程序為一初步程序，假設初值為10;10;10;10，則第一次運算得t0=0.0011,lamda=0.9291,迭代后的x=NaN?，F(xiàn)用共軛梯度法求解另一函數(shù) 對上述程序稍加改動來求解本題的代碼如下：注：程序所用

11、語言為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)程序結果如下：由上述結果知，用共軛梯度法對上述函數(shù)求解需要迭代三次得到最優(yōu)解0，此時x為0 0;符合理論分析的結果，求解完成。2.1 用外點法法求解函數(shù) 2.1.1 問題分析 外點法為序列無約束最優(yōu)化方法，其基本思想為將條件函數(shù)吸收到目標函數(shù)中進行求解。其在數(shù)學上的直觀理解

16、是拉格朗日乘子法：，為總代價，為價格，為罰款。即在經(jīng)濟學上總代價為價格和罰款的和。此時 稱為增廣目標函數(shù)，通常取 2.1.2 問題求解 兩種方法求解程序如下： 2.1.2.1 注：程序所用語言為MATLAB，終止條件為，程序中無約束優(yōu)化部分通過求導實現(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ù)?7次時能夠達到終止條件，此時, 函數(shù)得到最優(yōu)解0.368870.求解完成。2.1.2.2，程序中無約束優(yōu)化部分通過調(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ù)?次時能夠達到終止條件，此時, 函數(shù)得到最優(yōu)解0.368869.求解完成。2.2 用內(nèi)點法法求解函數(shù) 2.2.1 問題分析 同外點法，內(nèi)點法也為序列無約束最優(yōu)化方法，其基本思想為將條件函數(shù)吸收到目標函數(shù)中進行求解。 稱為障礙函數(shù)、圍墻函數(shù)或者懲罰函數(shù)。通常取 若在上一次迭代中x處于可解區(qū)域外部，則“圍墻”不夠高，在下一次迭代時加高“圍墻函數(shù)”。2.2.2 問題求解 求解程序如下： 注：程序所用語言為MATLAB，本程序出于對運行時間和算法簡單的考慮，選用終止條件為，

22、此時可以認為當x接近求解邊界時（圍墻）為無窮大，以確保函數(shù)的自變量迭代發(fā)生在可解區(qū)域內(nèi)部2.2.2.1，程序中無約束優(yōu)化部分通過求導實現(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ù)?5次時能夠達到終止條件，此時, 函數(shù)得到最優(yōu)解0.368870.求解完成，結果與“外點法”相同。2.2.2.2，程序中無約束優(yōu)化部分通過調(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