实用最优化方法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;whiled>elambdak=(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;whileresk>epsk=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)输出结果为-当=时程序如下：functionx_star=dfph(x0,eps)g0=grad(x0);res0=norm(g0);res=res0;k=0;xk=x0;whileres>epsk=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;whiled>elambdak=(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)输出结果为：当时程序如下：functionx_star=newton1(x0,eps)g0=grad(x0);res0=norm(g0);res=res0;k=0;xk=x0;whileres>epsk=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;whiled>elambdak=(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)输出结果为：取为BFGS公式时：程序如下：functionx_star=bfgs1(x0,eps)g0=grad(x0);res0=norm(g0);res=res0;k=0;xk=x0;whileres>epsk=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;whiled>elambdak=(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,:)*x>bi(i)+epsilon)index(i)=0;endendwhile(k<kmax)%求解子问题%Aee=[];if(ne>0)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,:)*dk<0))tm1=(bi(i)-Ai(i,:)*x)/(Ai(i,:)*dk);if(tm1<tm)tm=tm1;ti=i;endendendalpha=min(alpha,tm);x=x+alpha*dk;if(tm<1)%修正有效集%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(m>0)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,x0)输出结果为：程序如下：functionx_star=lagrange(x0,epsn)v=1;u=0;m=h(u,x0);whilem>epsnif(u+4*d(x0))>0u=u+4*d(x0);elseu=0;endv=v+4*e(x0);x0=bfgs(v,u,x0);m=h(u,x0);endx_tar=x0;x_tarfn=fn(x0);fnfunctionstar=bfgs(v0,u0,x0)eps=1e-4;g0=grad(v0,u0,x0);res0=norm(g0);res=res0;k=0;xk=x0;whileres>epsk=k+1;sk=-grad(v0,u0,x0);lambdak=nonexact(v0,u0,xk,sk);x0=xk;xk=x0+lambdak*sk;res=norm(grad(v0,u0,xk));endstar=xk;functionlambda=nonexact(v0,u0,x0,s0)g0=grad(v0,u0,x0);f0=f(v0,u0,x0);a=0;c1=0.1;c2=0.5;lambdak=1;p=-c1*lambdak*g0'*s0;xk=x0+lambdak*s0;f1=f(v0,u0,xk);q=f0-f1;whilep>qlambdak=(lambdak+a)/2;xk=x0+lambdak*s0;fk=f(v0,u0,xk);q=f0-fk;p=-c1*lambdak*g0'*s0;endlambda=lambdak;functiong=grad(v,u,x)g=zeros(2,1);g(1)=4-2*x(1)*v-8*x(1)*(25-x(1)^2-x(2)^2)+(u+4*(10*x(1)-x(1)^2+10*x(2)-x(2)^2-34))*(10-2*x(1));g(2)=-2*x(2)-2*x(2)*v-8*x(2)*(25-x(1)^2-x(2)^2)+(u+4*(10*x(1)-x(1)^2+10*x(2)-x(2)^2-34))*(10-2*x(2));functionf=f(v,u,x)f=4*x(1)-x(2)^2-12+v*(25-x(1)^2-x(2)^2)+2*(25-x(1)^2-x(2)^2)^2+0.125*((u+4*(10*x(1)-x(1)^2+10*x(2)-x(2)^2-34))^2-u^2);functionh=h(u,x)if(10*x(1)-x(1)^2+10*x(2)-x(2)^2-34)>(-u/4)h=(25-x(1)^2-x(2)^2)^2+(10*x(1)-x(1)^2+10*x(2)-x(2)^2-34)^2;elseh=(25-x(1)^2-x(2)^2)^2+(-u/4)^2;endfunctiond=d(x)d=10*x(1)-x(1)^2+10*x(2)-x(2)^2-34;functione=e(x)e=25-x(1)^2-x(2)^2;functionfn=fn(x)fn=4*x(1)-x(2)^2-12;在命令窗口中输入x0=[1;1];epsn=1e-4;lagrange(x0,epsn)输出结果为：先写出求解二次规划子问题的程序程序如下：function[d,mu,lam,val,k]=qpsubp(dfk,Bk,Ae,hk,Ai,gk)n=length(dfk);l=length(hk);m=length(gk);gamma=0.05;epsilon=1.0e-6;rho=0.5;sigma=0.2;ep0=0.05;mu0=0.05*zeros(l,1);lam0=0.05*zeros(m,1);d0=ones(n,1);u0=[ep0;zeros(n+l+m,1)];z0=[ep0;d0;mu0;lam0,];k=0;z=z0;ep=ep0;d=d0;mu=mu0;lam=lam0;while(k<=150)dh=dah(ep,d,mu,lam,dfk,Bk,Ae,hk,Ai,gk);if(norm(dh)<epsilon)break;endA=JacobiH(ep,d,mu,lam,dfk,Bk,Ae,hk,Ai,gk);b=beta(ep,d,mu,lam,dfk,Bk,Ae,hk,Ai,gk,gamma)*u0-dh;dz=A\b;if(l>0&m>0)de=dz(1);dd=dz(2:n+1);du=dz(n+2:n+l+1);dl=dz(n+l+2:n+l+m+1);endif(l==0)de=dz(1);dd=dz(2:n+1);dl=dz(n+2:n+m+1);endif(m==0)de=dz(1);dd=dz(2:n+1);du=dz(n+2:n+l+1);endi=0;mk=0;while(mk<=20)if(l>0&m>0)dh1=dah(ep+rho^i*de,d+rho^i*dd,mu+rho^i*du,lam+rho^i*dl,dfk,Bk,Ae,hk,Ai,gk);endif(l==0)dh1=dah(ep+rho^i*de,d+rho^i*dd,mu,lam+rho^i*dl,dfk,Bk,Ae,hk,Ai,gk);endif(m==0)dh1=dah(ep+rho^i*de,d+rho^i*dd,mu+rho^i*du,lam,dfk,Bk,Ae,hk,Ai,gk);endif(norm(dh1)<=(1-sigma*(1-gamma*ep0)*rho^i)*norm(dh))mk=i;break;endi=i+1;if(i==20)mk=10;endendalpha=rho^mk;if(l>0&m>0)ep=ep+alpha*de;d=d+alpha*dd;mu=mu+alpha*du;lam=lam+alpha*dl;endif(l==0)ep=ep+alpha*de;d=d+alpha*dd;lam=lam+alpha*dl;endif(m==0)ep=ep+alpha*de;d=d+alpha*dd;mu=mu+alpha*du;endk=k+1;endfunctionp=phi(ep,a,b)p=a+b-sqrt(a^2+b^2+2*ep^2);functiondh=dah(ep,d,mu,lam,dfk,Bk,Ae,hk,Ai,gk)n=length(dfk);l=length(hk);m=length(gk);dh=zeros(n+l+m+1,1);dh(1)=ep;if(l>0&m>0)dh(2:n+1)=Bk*d-Ae'*mu-Ai'*lam+dfk;dh(n+2:n+l+1)=hk+Ae*d;for(i=1:m)dh(n+l+1+i)=phi(ep,lam(i),gk(i)+Ai(i,:)*d);endendif(l==0)dh(2:n+1)=Bk*d-Ai'*lam+dfk;for(i=1:m)dh(n+1+i)=phi(ep,lam(i),gk(i)+Ai(i,:)*d);endendif(m==0)dh(2:n+1)=Bk*d-Ae'*mu+dfk;dh(n+2:n+l+1)=hk+Ae*d;enddh=dh(:);functionbet=beta(ep,d,mu,lam,dfk,Bk,Ae,hk,Ai,gk,gamma)dh=dah(ep,d,mu,lam,dfk,Bk,Ae,hk,Ai,gk);bet=gamma*norm(dh)*min(1,norm(dh));function[dd1,dd2,v1]=ddv(ep,d,lam,Ai,gk)m=length(gk);dd1=zeros(m,m);dd2=zeros(m,m);v1=zeros(m,1);for(i=1:m)fm=sqrt(lam(i)^2+(gk(i)+Ai(i,:)*d)^2+2*ep^2);dd1(i,i)=1-lam(i)/fm;dd2(i,i)=1-(gk(i)+Ai(i,:)*d)/fm;v1(i)=-2*ep/fm;endfunctionA=JacobiH(ep,d,mu,lam,dfk,Bk,Ae,hk,Ai,gk)n=length(dfk);l=length(hk);m=length(gk);A=zeros(n+l+m+1,n+l+m+1);[dd1,dd2,v1]=ddv(ep,d,lam,Ai,gk);if(l>0&m>0)A=[1,zeros(1,n),zeros(1,l),zeros(1,m);zeros(n,1),Bk,-Ae',-Ai';zeros(l,1),Ae,zeros(l,l),zeros(l,m);v1,dd2*Ai,zeros(m,l),dd1];endif(l==0)A=[1,zeros(1,n),zeros(1,m);zeros(n,1),Bk,-Ai';v1,dd2*Ai,dd1];endif(m==0)A=[1,zeros(1,n),zeros(1,l);zeros(n,1),Bk,-Ae';zeros(l,1),Ae,zeros(l,l)];end序列二次规划程序为：function[x,mu,lam,val,k]=sqpm(x0,mu0,lam0)maxk=100;%最大迭代次数n=length(x0);l=length(mu0);m=length(lam0);rho=0.5;eta=0.1;B0=eye(n);x=x0;mu=mu0;lam=lam0;Bk=B0;sigma=0.8;epsilon1=1e-6;epsilon2=1e-5;[hk,gk]=cons(x);dfk=df1(x);[Ae,Ai]=dcons(x);Ak=[Ae;Ai];k=0;while(k<maxk)[dk,mu,lam]=qpsubp(dfk,Bk,Ae,hk,Ai,gk);%求解子问题mp1=norm(hk,1)+norm(max(-gk,0),1);if(norm(dk,1)<epsilon1)&(mp1<epsilon2)break;end%检验终止准那么deta=0.05;%罚参数更新tau=max(norm(mu,inf),norm(lam,inf));if(sigma*(tau+deta)<1)sigma=sigma;elsesigma=1.0/(tau+2*deta);endim=0;%Armijo搜索while(im<=20)if(phi1(x+rho^im*dk,sigma)-phi1(x,sigma)<eta*rho^im*dphi1(x,sigma,dk))mk=im;break;endim=im+1;if(im==20)mk=10;endendalpha=rho^mk;x1=x+alpha*dk;[hk,gk]=cons(x1);dfk=df1(x1);[Ae,Ai]=dcons(x1);Ak=[Ae;Ai];lamu=pinv(Ak)'*dfk;%计算最小二乘乘子if(l>0&m>0)mu=lamu(1:l);lam=lamu(l+1:l+m);endif(l==0)mu=[];lam=lamu;endif(m==0)mu=lamu;lam=[];endsk=alpha*dk;%更新矩阵Bkyk=dlax(x1,mu,lam)-dlax(x,mu,lam);if(sk'*yk>0.2*sk'*Bk*sk)theta=1;elsetheta=0.8*sk'*Bk*sk/(sk'*Bk*sk-sk'*yk);endzk=theta*yk+(1-theta)*Bk*sk;Bk=Bk+zk*zk'/(sk'*zk)-(Bk*sk)*(Bk*sk)'/(sk'*Bk*sk);x=x1;k=k+1;endval=f1(x);%%%%l1精确价值函数%%%%%functionp=phi1(x,sigma)f=f1(x);[h,g]=cons(x);gn=max(-g,0);l0=length(h);m0=length(g);if(l0==0)p=f+1.0/sigma*norm(gn,1);endif(m0==0)p=f+1.0/