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

1、实用最优化方法Matlab程序设计程序如下：function lambda = 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;while de lambdak=(lambdak+a)/2; xk=x0 + lambdak*sk; fk=f(xk); e=f0-fk; d=-c1*lambdak*g0*sk;endlambdakfunction g=grad(x)g=zeros(2,1);

2、g(1)=-4*x(1)*(x(2)-x(1)2)-2*(1-x(1);g(2)=2*(x(2)-x(1)2);function fa=f(x)fa=(x(2)-x(1)2)2+(1-x(1)2;在命令窗口中输入x0=0;1; s0=-1;1;nonexact(x0,s0)输出结果为：程序如下：function x_star = cong(x0,eps)g0 = grad(x0);res0 = norm(g0);resk=res0;k = 0;xk =x0;while reskeps k=k+1; if k=1 sk=-grad(xk); else sk=-grad(xk)+resk2/res

3、02*s0; end lambdak=step(xk,sk); x0=xk; xk=x0 + lambdak*sk; res0=resk; resk=norm(grad(xk); s0=sk; fprintf(-the%d-th iteration, the residual is %f, the lambdak is %fnn, k ,resk,lambdak);endx_star=xk;disp(the optimal solution is);x_star%sub functionsfunction g=grad(x)g=zeros(4,1);g(1)=2*x(1)-2*x(2)+2;g

4、(2)=4*x(2)-2*x(1)-x(3)+3;g(3)=2*x(3)-x(2)-1;g(4)=2*x(4);function lambda=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,

5、eps)输出结果为- 当= 时程序如下：function x_star = dfph(x0,eps)g0 = grad(x0);res0 = norm(g0);res=res0;k = 0;xk =x0;while reseps k=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-th iteration, the residual is %f, the lambdak is %fnn, k

6、 ,res,lambdak);endx_star=xk;disp(the optimal solution is);x_starfunction lambda = 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;while de lambdak=(lambdak+a)/2; xk=x0 + lambdak*s0; fk=f(xk); e=f0-fk; d=-c1*lambdak*g0*s0;endla

7、mbda=lambdak;function g=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);function fa=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)输出结果为：当时程序如下：function x_star = newton1(x0,eps)g0 = grad(x0);res0 = norm(g0);res=res0;k = 0;xk =x0;while

8、 reseps k=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-th iteration, the residual is %f, the lambdak is %fnn, k ,res,lambdak);endx_star=xk;disp(the optimal solution is);x_starfunction lambda = nonexact(x0,s0)g0 = grad(x

9、0);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;while de lambdak=(lambdak+a)/2; xk=x0 + lambdak*s0; fk=f(xk); e=f0-fk; d=-c1*lambdak*g0*s0;endlambda=lambdak;function g=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)

10、2);function fa=f(x)fa=x(1)+2*x(2)2+exp(x(1)2+x(2)2);function G=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公式时：程序如下：fun

11、ction x_star = bfgs1(x0,eps)g0 = grad(x0);res0 = norm(g0);res=res0;k = 0;xk =x0;while reseps k=k+1; if k=1 H=eye(length(x0); else H0=H; mu=1+dg*H0*dg/(dx*dg); H=H0+(mu*dx*dx-H0*dg*dx-dx*dg*H0)/(dx*dg); end sk=-H*grad(xk); lambdak = nonexact(xk,sk); x0=xk; xk=x0 + lambdak*sk; dx=xk-x0; dg=grad(xk)-gr

12、ad(x0); res=norm(grad(xk); fprintf(-the%d-th iteration, the residual is %f, the lambdak is %fnn, k ,res,lambdak);endx_star=xk;fm=f(xk);disp(the optimal solution is);x_starfmfunction lambda = 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

13、(xk);e=f0-f1;while de lambdak=(lambdak+a)/2; xk=x0 + lambdak*s0; fk=f(xk); e=f0-fk; d=-c1*lambdak*g0*s0;endlambda=lambdak;function g=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);function fa=f(x)fa=x(1)+2*x(2)2+exp(x(1)2+x(2)2);在命令窗口中输入x0=0;1;eps=1e-3;bfgs1(x

14、0,eps)输出结果为：程序如下：function x=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 (k0) Aee=Ae; end for(j=1:ni) if(index(j)0) Aee=Aee; Ai(j,:);

15、 end end gk=H*x+c; m1,n1 = size(Aee); dk,lamk=qsubp(H,gk,Aee,zeros(m1,1); if(norm(dk)ne) y,jk=min(lamk(ne+1:length(lamk); end if(y=0) exitflag=0; else exitflag=1; for(i=1:ni) if(index(i)&(ne+sum(index(1:i)=jk) index(i)=0; break; end end end k=k+1; else exitflag=1; alpha=1.0; %求步长% tm=1.0; for(i=1:ni

16、) if(index(i)=0)&(Ai(i,:)*dk0) tm1=(bi(i)-Ai(i,:)*x)/(Ai(i,:)*dk); if(tm1tm) tm=tm1; ti=i; end end end alpha=min(alpha,tm); x=x+alpha*dk; if(tm0) rb=Ae*ginvH*c + be; lambda=pinv(Ae*ginvH*Ae)*rb; x=ginvH*(Ae*lambda-c);else x=-ginvH*c; lambda=0;endfunction fn=f(x)fn=(x(1)-1)2+(x(2)-2.5)2;在命令窗口中输入H=2 0

17、;0 2;c=-2;-5;Ae= ; be= ;Ai=1 -2;-1 -2;-1 2;1 0;0 1;bi=-2;-6;-2;0 ;0;x0=0;0;qpact(H,c,Ae,be,Ai,bi,x0)输出结果为：程序如下：function x_star = lagrange(x0,epsn)v=1;u=0;m=h(u,x0);while mepsn if (u+4*d(x0)0 u=u+4*d(x0); else u=0; end v=v+4*e(x0); x0=bfgs(v,u,x0); m=h(u,x0);endx_tar=x0;x_tarfn=fn(x0);fnfunction star

18、 = bfgs(v0,u0,x0)eps=1e-4;g0 = grad(v0,u0,x0);res0 = norm(g0);res=res0;k = 0;xk =x0;while reseps k=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;function lambda = nonexact(v0,u0,x0,s0)g0 = grad(v0,u0,x0);f0=f(v0,u0,x0);a=0;c1=

19、0.1;c2=0.5;lambdak=1;p=-c1*lambdak*g0*s0;xk=x0 + lambdak*s0;f1=f(v0,u0,xk);q=f0-f1;while pq lambdak=(lambdak+a)/2; xk=x0 + lambdak*s0; fk=f(v0,u0,xk); q=f0-fk; p=-c1*lambdak*g0*s0;endlambda=lambdak;function g=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)

20、-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);function f=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-u2);function h=h(u,x)if (10*x(1)-x(1)2+10*x(2)-x(2)2-34)(-u/4)

21、 h=(25-x(1)2-x(2)2)2+(10*x(1)-x(1)2+10*x(2)-x(2)2-34)2;else h=(25-x(1)2-x(2)2)2+(-u/4)2;endfunction d=d(x)d=10*x(1)-x(1)2+10*x(2)-x(2)2-34;function e=e(x)e=25-x(1)2-x(2)2;function fn=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=

22、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,

23、Ai,gk); if(norm(dh)0&m0) de=dz(1); dd=dz(2:n+1); du=dz(n+2:n+l+1); dl=dz(n+l+2:n+l+m+1); end if(l=0) de=dz(1); dd=dz(2:n+1); dl=dz(n+2:n+m+1); end if(m=0) de=dz(1); dd=dz(2:n+1); du=dz(n+2:n+l+1); end i=0; mk=0; while (mk0&m0) dh1=dah(ep+rhoi*de,d+rhoi*dd,mu+rhoi*du,lam+rhoi*dl,dfk,Bk,Ae,hk,Ai,gk);

24、end if(l=0) dh1=dah(ep+rhoi*de,d+rhoi*dd,mu,lam+rhoi*dl,dfk,Bk,Ae,hk,Ai,gk); end if(m=0) dh1=dah(ep+rhoi*de,d+rhoi*dd,mu+rhoi*du,lam,dfk,Bk,Ae,hk,Ai,gk); end if(norm(dh1)0&m0) ep=ep+alpha*de; d=d+alpha*dd; mu=mu+alpha*du; lam=lam+alpha*dl; end if(l=0) ep=ep+alpha*de; d=d+alpha*dd; lam=lam+alpha*dl;

25、end if(m=0) ep=ep+alpha*de; d=d+alpha*dd; mu=mu+alpha*du; end k=k+1;endfunction p=phi(ep,a,b)p=a+b-sqrt(a2+b2+2*ep2);function dh=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(l0&m0) dh(2:n+1)=Bk*d-Ae*mu-Ai*lam+dfk; dh(n+2:n+l+1)=hk+Ae*d;

26、 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(:);function bet=beta(ep,d,mu,lam,dfk,Bk,Ae,hk,Ai,gk,gamma)dh=dah(ep,d,mu,lam,df

27、k,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*ep2); dd1(i,i)=1-lam(i)/fm; dd2(i,i)=1-(gk(i)+Ai(i,:)*d)/fm; v1(i)=-2*ep/fm;endfunction A=JacobiH(ep,d

28、,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(l0&m0) 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), z

29、eros(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

30、;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(kmaxk) dk,mu,lam=qpsubp(dfk,Bk,Ae,hk,Ai,gk); %求解子问题 mp1=norm(hk,1)+norm(max(-gk,0),1); if(norm(dk,1)epsilon1)&(mp1epsilon2) break; end %检验终止准则 deta=0.05; %罚参数更新 tau=max(norm(mu,inf),norm(lam,i

31、nf); if(sigma*(tau+deta)1) sigma=sigma; else sigma=1.0/(tau+2*deta); end im=0; %Armijo搜索 while(im=20) if(phi1(x+rhoim*dk,sigma)-phi1(x,sigma)0&m0) mu=lamu(1:l); lam=lamu(l+1:l+m); end if(l=0) mu=; lam=lamu; end if(m=0) mu=lamu; lam=; end sk=alpha*dk; %更新矩阵Bk yk=dlax(x1,mu,lam)-dlax(x,mu,lam); if(sk*

32、yk0.2*sk*Bk*sk) theta=1; else theta=0.8*sk*Bk*sk/(sk*Bk*sk-sk*yk); end zk=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精确价值函数%function p=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/sigma*norm(h,1);endif(l00&m00) p=f+1.0/sigma*(norm(h,1)+norm(gn,1);end%价值函数的方向导数%function dp=