实用最优化方法Matlab程序设计

实用最优化方法Matlab程序设计程序如下：functionlambda=nonexact(x0,s0)g0=grad(x0);f0=f(x0);a=0;c1=0.1;c2=0.5;lambdak=1;sk=s0;d=-c1*lambdak*g0'*sk;xk=x0+lambdak*sk;f1=f(xk);e=f0-f1;whiledelambdak=(lambdak+a)/2;xk=x0+lambdak*sk;fk=f(xk);e=f0-fk;d=-c1*lambdak*g0'*sk;endlambdakfunctiong=grad(x)g=zeros(2,1);g(1)=-4*x(1)*(x(2)-x(1)^2)-2*(1-x(1));g(2)=2*(x(2)-x(1)^2);functionfa=f(x)fa=(x(2)-x(1)^2)^2+(1-x(1))^2;在命令窗口中输入x0=[0;1];s0=[-1;1];nonexact(x0,s0)输出结果为：程序如下：functionx_star=cong(x0,eps)g0=grad(x0);res0=norm(g0);resk=res0;k=0;xk=x0;whilereskepsk=k+1;ifk==1sk=-grad(xk);elsesk=-grad(xk)+resk^2/res0^2*s0;endlambdak=step(xk,sk);x0=xk;xk=x0+lambdak*sk;res0=resk;resk=norm(grad(xk));s0=sk;fprintf('------the%d-thiteration,theresidualis%f,thelambdakis%f\n\n',k,resk,lambdak);endx_star=xk;disp('theoptimalsolutionis');x_star%subfunctionsfunctiong=grad(x)g=zeros(4,1);g(1)=2*x(1)-2*x(2)+2;g(2)=4*x(2)-2*x(1)-x(3)+3;g(3)=2*x(3)-x(2)-1;g(4)=2*x(4);functionlambda=step(x,s)a=2*x(1)*s(1)-2*x(2)*s(1)-2*x(1)*s(2)+4*x(2)*s(2)+2*s(3)*x(3)+2*x(4)*s(4)-s(2)*x(3)-s(3)*x(2)+2*s(1)+3*s(2)-s(3);b=2*s(1)^2+4*s(2)^2-4*s(1)*s(2)+2*s(3)^2+2*s(4)^2-2*s(2)*s(3);lambda=-a/b;在命令窗口中输入x0=[0;0;0;0];eps=1e-6;cong(x0,eps)输出结果为-当kH=2I时程序如下：functionx_star=dfph(x0,eps)g0=grad(x0);res0=norm(g0);res=res0;k=0;xk=x0;whileresepsk=k+1;H=eye(length(x0));sk=-H*grad(xk);lambdak=nonexact(xk,sk);x0=xk;xk=x0+lambdak*sk;res=norm(grad(xk));fprintf('------the%d-thiteration,theresidualis%f,thelambdakis%f\n\n',k,res,lambdak);endx_star=xk;disp('theoptimalsolutionis');x_starfunctionlambda=nonexact(x0,s0)g0=grad(x0);f0=f(x0);a=0;c1=0.1;c2=0.5;lambdak=1;d=-c1*lambdak*g0'*s0;xk=x0+lambdak*s0;f1=f(xk);e=f0-f1;whiledelambdak=(lambdak+a)/2;xk=x0+lambdak*s0;fk=f(xk);e=f0-fk;d=-c1*lambdak*g0'*s0;endlambda=lambdak;functiong=grad(x)g=zeros(2,1);g(1)=2*x(1)*exp(x(1)^2+x(2)^2)+1;g(2)=4*x(2)+2*x(2)*exp(x(1)^2+x(2)^2);functionfa=f(x)fa=x(1)+2*x(2)^2+exp(x(1)^2+x(2)^2);在命令窗口中输入x0=[0;1];eps=1e-3;dfph(x0,eps)输出结果为：当12()kkHfx时程序如下：functionx_star=newton1(x0,eps)g0=grad(x0);res0=norm(g0);res=res0;k=0;xk=x0;whileresepsk=k+1;H=inv(GRAND(xk));sk=-H*grad(xk);lambdak=nonexact(xk,sk);x0=xk;xk=x0+lambdak*sk;res=norm(grad(xk));fprintf('------the%d-thiteration,theresidualis%f,thelambdakis%f\n\n',k,res,lambdak);endx_star=xk;disp('theoptimalsolutionis');x_starfunctionlambda=nonexact(x0,s0)g0=grad(x0);f0=f(x0);a=0;c1=0.1;c2=0.5;lambdak=1;d=-c1*lambdak*g0'*s0;xk=x0+lambdak*s0;f1=f(xk);e=f0-f1;whiledelambdak=(lambdak+a)/2;xk=x0+lambdak*s0;fk=f(xk);e=f0-fk;d=-c1*lambdak*g0'*s0;endlambda=lambdak;functiong=grad(x)g=zeros(2,1);g(1)=2*x(1)*exp(x(1)^2+x(2)^2)+1;g(2)=4*x(2)+2*x(2)*exp(x(1)^2+x(2)^2);functionfa=f(x)fa=x(1)+2*x(2)^2+exp(x(1)^2+x(2)^2);functionG=GRAND(x)G=zeros(2,2);G(1)=2*exp(x(1)^2+x(2)^2)+4*x(1)^2*exp(x(1)^2+x(2)^2);G(2)=4*x(1)*x(2)*exp(x(1)^2+x(2)^2);G(3)=4*x(1)*x(2)*exp(x(1)^2+x(2)^2);G(4)=4+(2+4*x(2)^2)*exp(x(1)^2+x(2)^2);在命令窗口中输入x0=[0;1];eps=1e-3;newton1(x0,eps)输出结果为：取kH为BFGS公式时：程序如下：functionx_star=bfgs1(x0,eps)g0=grad(x0);res0=norm(g0);res=res0;k=0;xk=x0;whileresepsk=k+1;ifk==1H=eye(length(x0));elseH0=H;mu=1+dg'*H0*dg/(dx'*dg);H=H0+(mu*dx*dx'-H0*dg*dx'-dx*dg'*H0)/(dx'*dg);endsk=-H*grad(xk);lambdak=nonexact(xk,sk);x0=xk;xk=x0+lambdak*sk;dx=xk-x0;dg=grad(xk)-grad(x0);res=norm(grad(xk));fprintf('------the%d-thiteration,theresidualis%f,thelambdakis%f\n\n',k,res,lambdak);endx_star=xk;fm=f(xk);disp('theoptimalsolutionis');x_starfmfunctionlambda=nonexact(x0,s0)g0=grad(x0);f0=f(x0);a=0;c1=0.1;c2=0.5;lambdak=1;d=-c1*lambdak*g0'*s0;xk=x0+lambdak*s0;f1=f(xk);e=f0-f1;whiledelambdak=(lambdak+a)/2;xk=x0+lambdak*s0;fk=f(xk);e=f0-fk;d=-c1*lambdak*g0'*s0;endlambda=lambdak;functiong=grad(x)g=zeros(2,1);g(1)=2*x(1)*exp(x(1)^2+x(2)^2)+1;g(2)=4*x(2)+2*x(2)*exp(x(1)^2+x(2)^2);functionfa=f(x)fa=x(1)+2*x(2)^2+exp(x(1)^2+x(2)^2);在命令窗口中输入x0=[0;1];eps=1e-3;bfgs1(x0,eps)输出结果为：程序如下：functionx=qpact(H,c,Ae,be,Ai,bi,x0)epsilon=1.0e-9;err=1.0e-6;k=0;x=x0;n=length(x);kmax=1.0e3;ne=length(be);ni=length(bi);lamk=zeros(ne+ni,1);index=ones(ni,1);for(i=1:ni)if(Ai(i,:)*xbi(i)+epsilon)index(i)=0;endendwhile(kkmax)%求解子问题%Aee=[];if(ne0)Aee=Ae;endfor(j=1:ni)if(index(j)0)Aee=[Aee;Ai(j,:)];endendgk=H*x+c;[m1,n1]=size(Aee);[dk,lamk]=qsubp(H,gk,Aee,zeros(m1,1));if(norm(dk)=err)y=0.0;if(length(lamk)ne)[y,jk]=min(lamk(ne+1:length(lamk)));endif(y=0)exitflag=0;elseexitflag=1;for(i=1:ni)if(index(i)&(ne+sum(index(1:i)))==jk)index(i)=0;break;endendendk=k+1;elseexitflag=1;alpha=1.0;%求步长%tm=1.0;for(i=1:ni)if((index(i)==0)&(Ai(i,:)*dk0))tm1=(bi(i)-Ai(i,:)*x)/(Ai(i,:)*dk);if(tm1tm)tm=tm1;ti=i;endendendalpha=min(alpha,tm);x=x+alpha*dk;if(tm1)%修正有效集%index(ti)=1;endendif(exitflag==0)break;endk=k+1;endfn=f(x);x_star=x;x_starfn%求解子问题%function[x,lambda]=qsubp(H,c,Ae,be)ginvH=pinv(H);[m,n]=size(Ae);if(m0)rb=Ae*ginvH*c+be;lambda=pinv(Ae*ginvH*Ae')*rb;x=ginvH*(Ae'*lambda-c);elsex=-ginvH*c;lambda=0;endfunctionfn=f(x)fn=(x(1)-1)^2+(x(2)-2.5)^2;在命令窗口中输入H=[20;02];c=[-2;-5];Ae=[];be=[];Ai=[1-2;-1-2;-12;10;01];bi=[-2;-6;-2;0;0];x0=[0;0];qpact(H,c,Ae,be,Ai,bi