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

实用最优化方法编程大作业班级：姓名：指导老师：学号：大连理工大学2015年11月27日1版本号：MATLAB7.11.0(R2010b)【文件名WP.m】functionx=WP(f,x0,var,s0,eps)clcsymsx1x2;c1=0.1;c2=0.5;b=inf;lambda=1;ifnargin==4eps=1.0e-6;endgradf=jacobian(f,var);g0=subs(gradf,var,x0);f0=subs(f,var,x0);gs=g0*s0';a=0;j=0;whilej1000new_x=x0+lambda*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;ifleftright%不满足第一个条件j=j+1;b=lambda;lambda=0.5*(lambda+a);else%满足第一个条件left2=new_gs;right2=c2*g0;ifleft2right2%不满足第二个条件j=j+1;a=lambda;2lambda=min(2*lambda,0.5*(lambda+b));elsex=lambda;break;endendend在CommandWindow输入：symsx1x2WP(100*(x2-x1^2)^2+(1-x1)^2,[-1,1],[x1x2],[11])程序运行后可得出结果：ans=0.0039【文件名minGETD.m】function[x,minf]=minGETD(f,x0,var,eps)clcsymsx1x2x3formatlong;n=length(x0);ifnargin==3eps=1.0e-3;endx0=x0';symslambda;gradf=jacobian(f,var);g=subs(gradf,var,x0);s=-g';k=0;while1tol=norm(double(g));iftol=epsx=x0;break;3endx1=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))iftol=epsx=x1;break;endifk+1==nx0=x1;continue;elsemiu=dot(g1,g1)/dot(g,g)s=-g1'+miu*s;k=k+1;x0=x1;endendx在CommandWindow输入：symsx1x2x3x0=[000];var=[x1x2x3];f=x1^2-2*x1*x2+2*x2^2+x3^2-x2*x3+2*x1+3*x2-x3;minGETD(f,x0,var,eps)程序运行后可得出结果：x=[-236894/59711,-178563/59711,-59465/59711]可认为最终解为[-4,-3,-1]。4【文件名Excise.m】function[x,minf]=Excise_3(f,x0,var,method,eps)%method表算法，1时为最速下降法，2时为牛顿法，3为BFGSclcticsymsx1x2formatlong;ifnargin==4eps=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));iftol=epsx=x0;minf=subs(f,var,x);return;endwhile1lambda0=WP(f,x0',var,s');x1=x0+lambda0*s;5f1=subs(f,var,x1);g1=subs(gradf,var,x1)';tol=norm(double(g1));iftol=epsx=x1;break;endifk+1==nx0=x1;H0=eye(2);s=-subs(gradf,var,x0)';k=0;continue;elseswitchmethodcase1%最速下降法H1=eye(n);case2%牛顿法Hesse=jacobian(gradf,var)H1=inv(subs(Hesse,var,x1));case3%BFGSd_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'endj=j+1;s=-H1*g1;k=k+1;x0=x1;g0=g1;H0=H1;endendx%最优解minf=subs(f,var,x)%最优值j%迭代次数toc%运行时间formatshort;在CommandWindow输入：symsx1x26Excise_3(x1+2*x2^2+exp(x1^2+x2^2),[10],[x1x2],1)程序运行后可得出结果（最速下降法，method=1）：x=[-0.419364591338560,0]minf=0.772914478780059j=10Elapsedtimeis2.328044seconds.仅改变输入变量method的值得到牛顿法（method=2），程序运行后可得出结果：x=[-0.419365009463280,0]minf=0.772914478780027j=3Elapsedtimeis0.876141seconds.仅改变输入变量method的值得到BFGS法（method=3），程序运行后可得出结果：x=[-0.419364824015761,0]minf=0.772914478779971j=8Elapsedtimeis3.731816seconds.可以看出，仅就习题3来说，牛顿法迭代次数（3次）最少，运行时间（0.87s）最短。【文件名Excise_4.m】function[x,minf]=Excise_4(x0)clcformatlong;symsx1x2;x=[x1x2]';f=x1^2+2*x2^2-2*x1-6*x2-2*x1*x2;G=[2-2;-24];g=[-2-6]';b=[1200]';AA=[1/21/2;-12;-10;0-1];A=AA;A1=[];A0=[];b0=[];b1=[];7fori=1:length(A(:,1))ifA(i,:)*x0==b(i)A10=A(i,:);A1=[A1;A10];b1=[b1;b(i)];elseA00=A(i,:);A0=[A0;A00];b0=[b0;b(i)];endendwhile1ifsize(A1)==[00]ZA=G;Zg=-(G*x0+g);dl=inv(ZA)*Zg;d=dl;lambda=dl(3:length(dl),:);elseZA=[GA1';A1zeros(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),:);endifnorm(d)1.0e-6d=[00]';endlmin=min(lambda);ifd==zeros(2,1)iflmin=0break;elsea=find(lambda==lmin);A0=[A0;A1(a,:)];b0=[b0;b1(a)];A1(a,:)=[];b1(a)=[];endelsex0ba=x0+d;newb=A0*x0ba-b0;8bmax=max(newb);ifbmax=0x0=x0ba;elseA2=[];b2=[];fori=1:length(A0(:,1))ifA0(i,:)*d0A2=[A2;A0(i,:)];b2=[b2;b0(i)];endendalfav=(A2*x0-b2)./(A2*d);alfa=min(-alfav);x0=x0+alfa.*d;endA=AA;A1=[];A0=[];b0=[];b1=[];fori=1:length(A(:,1))ifA(i,:)*x0==b(i)A10=A(i,:);A1=[A1;A10];b1=[b1;b(i)];elseA00=A(i,:);A0=[A0;A00];b0=[b0;b(i)];endendendendx=x0minf=subs(f,findsym(f),x0)习题要求取两种点：（1）可行域内部点（取x0=[00]’）；(2)可行域边界点（取x0=[2/34/3]’）。在CommandWindow输入：Excise_4([00]’)程序运行后可得出结果：x=[0.8,1.2]minf=-7.2同理，在CommandWindow输入：Excise_4([2/34/3]’)可得相同结果。9【文件名Excise_5.m（脚本）】clearclcformatlong;symsx1x2lambdamiu1miu2miu3eps=1.0e-6;c=8;miu0=[000]';lambda=0;x=[x1x2]';miu=[miu1miu2miu3]';f=4*x1-x2^2-12;h=25-x1^2-x2^2;g1=10*x1-x1^2-x2^2+10*x2-34;g2=x1;g3=x2;zg=[g1g2g3]';comp=[000]';x0=[00]';while1zg0=subs(zg,findsym(zg),x0)y=Lagrange(x,c,miu0,lambda,x0)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)ifh1^2+(norm(b))^2eps^2break;elselambda=lambda+c*h1miu0=min(comp,miu0+c*zg1)x0=x1end10endminx=x1minf=subs(f,findsym(f),x1)【文件名Lagrange.m(计算增广的Lagrange表达式)】functiony=Lagrange(x,c,miu,lambda,x0);f1=4*x(1)-x(2)^2-12;f=[f10]';h1=25-x(1)^2-x(2)^2;h=[h10]';g1=10*x(1)-x(1)^2-x(2)^2+10*x(2)-34;g2=x(1);g3=x(2);zg=[g1g2g3]';comp=[000]';y=f(1)+lambda*h(1)+c/2*(h(1)^2);zg0=subs(zg,findsym(zg),x0);fori=1:3ifmiu(i)+c*zg0(i)0y=y+1/(2*c)*((miu(i)+c*zg(i))^2-miu(i)^2);elsey=y-1/(2*c)*(miu(i)^2);endend【文件名newWP.m（精简第一题WP.m）】functiona=newWP(x,x0,s,c,miu,lambda)eps=1.0e-6;c1=0.1;c2=0.5;a=0;b=inf;lambda0=1;while1f=Lagrange(x,c,miu,lambda,x0);gf=jacobian(f,x');f0=subs(f,findsym(f),x0);gf0=subs(gf,findsym(gf),x0)';x1=x0