实用最优化方法编程大作业

1、实用最优化方法编程大作业班 级： 姓 名： 指导老师： 学 号： 大连理工大学2015年11月27日版本号：MATLAB 7.11.0 (R2010b)【文件名 WP.m】function x = WP(f,x0,var,s0,eps)clcsyms x1 x2 ;c1=0.1;c2=0.5;b=inf;lambda=1;if nargin=4 eps=1.0e-6;endgradf=jacobian(f,var); g0=subs(gradf,var,x0);f0=subs(f,var,x0);gs=g0*s0'a=0;j=0;while j<1000 new_x=x0+lam

2、bda*s0; new_f=subs(f,var,new_x); left=f0-new_f; new_g=subs(gradf,var,new_x); new_gs=new_g*s0' right=-1*c1*lambda*gs; if left<right %不满足第一个条件 j=j+1; b=lambda; lambda=0.5*(lambda+a); else %满足第一个条件 left2=new_gs; right2=c2*g0; if left2<right2 %不满足第二个条件 j=j+1; a=lambda; lambda=min(2*lambda,0.5*

3、(lambda+b); else x=lambda; break; end endend在Command Window 输入：syms x1 x2WP(100*(x2-x12)2+(1-x1)2,-1,1,x1 x2,1 1)程序运行后可得出结果：ans=0.0039【文件名minGETD.m】function x,minf = minGETD(f,x0,var,eps)clcsyms x1 x2 x3format long;n=length(x0);if nargin=3 eps=1.0e-3;endx0=x0'syms lambda;gradf=jacobian(f,var);g=

4、subs(gradf,var,x0);s=-g'k=0;while 1 tol=norm(double(g); if tol<=eps x=x0; break; end x1=x0+lambda*s; f1=subs(f,var,x1); dy1=diff(f1); lambda0=solve(dy1); x1=x0+lambda0*s; g1=subs(gradf,var,x1) tol=norm(double(g1) if tol<=eps x=x1; break; end if k+1=n x0=x1; continue; else miu=dot(g1,g1)/do

5、t(g,g) s=-g1'+miu*s; k=k+1; x0=x1; endendx在Command Window 输入：syms x1 x2 x3x0=0 0 0;var=x1 x2 x3;f=x12-2*x1*x2+2*x22+x32-x2*x3+2*x1+3*x2-x3;minGETD(f,x0,var,eps)程序运行后可得出结果：x=-236894/59711,-178563/59711,-59465/59711可认为最终解为-4,-3,-1。【文件名Excise.m】function x,minf = Excise_3(f,x0,var,method,eps) % meth

6、od表算法，1时为最速下降法，2时为牛顿法，3为BFGSclcticsyms x1 x2format long;if nargin = 4 eps = 1.0e-6;endn=length(x0);H0=eye(2);x0=x0'gradf=jacobian(f,var);g0=subs(gradf,var,x0)'s=-H0*g0;k=0;j=0; % 检查初始点是否为最优点tol=norm(double(g0);if tol<=eps x=x0; minf=subs(f,var,x); return;end while 1 lambda0=WP(f,x0',v

7、ar,s'); x1=x0+lambda0*s; f1=subs(f,var,x1); g1=subs(gradf,var,x1)' tol=norm(double(g1); if tol<=eps x=x1; break; end if k+1=n x0=x1; H0=eye(2); s=-subs(gradf,var,x0)' k=0; continue; else switch method case 1 % 最速下降法 H1=eye(n); case 2 % 牛顿法 Hesse=jacobian(gradf,var) H1=inv(subs(Hesse,v

8、ar,x1); case 3 % BFGS d_x=x1-x0; d_g=g1-g0; v=d_x/(d_x'*d_g)-H0*d_g/(d_g'*H0*d_g); H1=H0+(d_x*d_x')/(d_x'*d_g)-(H0*d_g)*(H0*d_g)')/. (d_g'*H0*d_g) +(d_g'*H0*d_g)*v*v' end j=j+1; s=-H1*g1; k=k+1; x0=x1; g0=g1; H0=H1; endendx % 最优解minf=subs(f,var,x) % 最优值j % 迭代次数toc % 运

9、行时间format short;在Command Window 输入：syms x1 x2Excise_3(x1+2*x22+exp(x12+x22),1 0,x1 x2,1)程序运行后可得出结果（最速下降法，method=1）：x=-0.419364591338560,0minf=0.772914478780059j=10Elapsed time is 2.328044 seconds.仅改变输入变量method的值得到牛顿法（method=2），程序运行后可得出结果：x=-0.419365009463280,0minf=0.772914478780027j=3Elapsed time is

10、 0.876141 seconds.仅改变输入变量method的值得到BFGS法（method=3），程序运行后可得出结果：x=-0.419364824015761,0minf=0.772914478779971j=8Elapsed time is 3.731816 seconds.可以看出，仅就习题3来说，牛顿法迭代次数（3次）最少，运行时间（0.87s）最短。【文件名Excise_4.m】function x,minf=Excise_4(x0)clcformat long;syms x1 x2;x=x1 x2'f=x12+2*x22-2*x1-6*x2-2*x1*x2;G=2 -2

11、;-2 4;g=-2 -6'b=1 2 0 0'AA=1/2 1/2; -1 2; -1 0; 0 -1;A=AA; A1=; A0=; b0=; b1=; for i=1:length(A(:,1) if A(i,:)*x0=b(i) A10=A(i,:); A1=A1;A10; b1=b1;b(i); else A00=A(i,:); A0=A0;A00; b0=b0;b(i); endend while 1 if size(A1)=0 0 ZA=G; Zg=-(G*x0+g); dl=inv(ZA)*Zg; d=dl; lambda=dl(3:length(dl),:);

12、 else ZA=G A1'A1 zeros(length(A1(:,1); Zg=-(G*x0+g);zeros(length(A1(:,1),1); dl=inv(ZA)*Zg; d=dl(1:2,:); lambda=dl(3:length(dl),:); end if norm(d)<1.0e-6 d=0 0' end lmin=min(lambda); if d=zeros(2,1) if lmin>=0 break; else a=find(lambda=lmin); A0=A0;A1(a,:); b0=b0;b1(a); A1(a,:)=; b1(a)

13、=; end else x0ba=x0+d; newb=A0*x0ba-b0; bmax=max(newb); if bmax<=0 x0=x0ba; else A2=; b2=; for i=1:length(A0(:,1) if A0(i,:)*d>0 A2=A2;A0(i,:); b2=b2;b0(i); end end alfav=(A2*x0-b2)./(A2*d); alfa=min(-alfav); x0=x0+alfa.*d; end A=AA; A1=; A0=; b0=; b1=; for i=1:length(A(:,1) if A(i,:)*x0=b(i)

14、A10=A(i,:); A1=A1;A10; b1=b1;b(i); else A00=A(i,:); A0=A0;A00; b0=b0;b(i); end end endendx=x0minf=subs(f,findsym(f),x0)习题要求取两种点：（1）可行域内部点（取x0=0 0）；(2)可行域边界点（取x0=2/3 4/3）。在Command Window 输入：Excise_4(0 0)程序运行后可得出结果：x=0.8,1.2minf=-7.2同理，在Command Window 输入：Excise_4(2/3 4/3)可得相同结果。【文件名Excise_5.m （脚本）】cle

15、arclcformat long;syms x1 x2 lambda miu1 miu2 miu3eps=1.0e-6;c=8; miu0=0 0 0' lambda=0;x=x1 x2'miu=miu1 miu2 miu3'f=4*x1-x22-12; h=25-x12-x22;g1=10*x1-x12-x22+10*x2-34;g2=x1; g3=x2;zg=g1 g2 g3'comp=0 0 0'x0=0 0'while 1 zg0=subs(zg,findsym(zg),x0) y=Lagrange(x,c,miu0,lambda,x0)

16、 x1=BFGS(x0,x,c,miu0,lambda) h1=subs(h,findsym(h),x1) zg1=subs(zg,findsym(zg),x1) b=min(zg1,-1/c*miu0) if h12+(norm(b)2<eps2 break; else lambda=lambda+c*h1 miu0=min(comp,miu0+c*zg1) x0=x1 endendminx=x1minf=subs(f,findsym(f),x1)【文件名Lagrange.m (计算增广的Lagrange表达式)】function y=Lagrange(x,c,miu,lambda,x

17、0);f1=4*x(1)-x(2)2-12; f=f1 0'h1=25-x(1)2-x(2)2; h=h1 0'g1=10*x(1)-x(1)2-x(2)2+10*x(2)-34;g2=x(1); g3=x(2);zg=g1 g2 g3'comp=0 0 0'y=f(1)+lambda*h(1)+c/2*(h(1)2);zg0=subs(zg,findsym(zg),x0);for i=1:3 if miu(i)+c*zg0(i)<0 y=y+1/(2*c)*(miu(i)+c*zg(i)2-miu(i)2); else y=y-1/(2*c)*(miu(

18、i)2); endend【文件名newWP.m（精简第一题WP.m）】function a=newWP(x,x0,s,c,miu,lambda)eps=1.0e-6;c1=0.1;c2=0.5;a=0;b=inf;lambda0=1;while 1 f=Lagrange(x,c,miu,lambda,x0); gf=jacobian(f, x'); f0=subs(f,findsym(f),x0); gf0=subs(gf,findsym(gf),x0)' x1=x0+lambda0*s; f1=subs(f,findsym(f),x1); gf1=subs(gf,findsy

19、m(gf),x1); left1=f0-f1; right1=-c1*lambda0*(gf0'*s); left2=gf1*s; right2=c2*(gf0'*s); if left1>=right1 if left2>=right2 break; else a=lambda0; lambda1=min(2*lambda0,(lambda0+b)/2); end else b=lambda0; lambda1=(lambda0+a)/2; end lambda0=lambda1;enda=lambda0;【文件名BFGS.m (精简习题3的method=3为适应乘子法的BFGS方法)】function minx=BFGS(x0,x,c,miu,lambda)s