实用最优化方法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 delambdak=(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);g(1)=

2、-4*x(1)*(x(2)-x(1)人2)-2*(1-x(1);g(2)=2*(x(2)-x (1)睥function fa=f(x)fa=(x(2)-x(1)A2)A2+(1-x(1)A2;在命令窗 口中输入 x0=0;1; s0=-1;1;nonexact(x0,s0)输出结果为： kO= 0;1; s0=-1;1;nonesact(kO sO) lambdak =1.1102e-016程序如下：function x_star = cong(x0,eps)g0 = grad(x0);res0 = norm(g0);resk=res0;k = 0;xk =x0;while reskepsk=

3、k+1;if k=1sk=-grad(xk);elsesk=-grad(xk)+reskA2/res0A2*s0;endlambdak=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=g

4、rad(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);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);lam

5、bda=-a/b;在命令窗口中输入 x0=0;0;0;0;eps=1e-6;cong(x0,eps)输出结果为 k0= 0;0;0;0;eps=le-6;cong(kOj eps)thel thit erat iorijtheresidual is3. 674235,the lambdakis0.500000the2 thit erat iorijtheresidual is1. 896238,the lambdakis0.482759the3 thiteration,theresidual is0. 000000,the lambdakis0.690476the optimal soluti

6、on isK_star =-4.0000-3.0000-1. 0000 ans =-4.0000-3.0000-1. 00000当H k = 2时程序如下：function x_star = dfph(x0,eps)g0 = grad(x0);res0 = norm(g0);res=res0;k = 0;xk =x0;while resepsk=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 iter

7、ation, 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(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 delambdak=(lambdak+a)/2;xk=x0 + lambdak*s0

8、;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)A2+x(2)A2)+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)A2+exp(x(1)A2+x(2)A2);在命令窗口中输入 x0=0;1;eps=1e-3; dfph(x0,eps)输出结果为： kO= LO;1J ;eps=le-3;dfph(kO, eps)thel thitera

9、tion,theresidual is1.320364,the lambdakis 0.125000the2 thit erat ioiijtheresidual is0. 664208,the lambdakis 0.250000the3 thit erat iorijtheresidual is0. 344799,the lambdakis 0.250000the4 thit erat iorijtheresidual is0. 195525,the lambdakis 0.250000the5 thit erat iorijtheresidual is0. 115409,the lamb

10、dakis 0.250000the6 thiteration,theresidual is0.068684,the lambdakis 0.250000the? thiteration,theresidual is0. 040933,the lambdakis 0.250000the8 thiteration,theresidual is0. 024400,the lambdakis 0.250000the9 thit erat iorijtheresidual is0. 014546,the lambdakis 0.250000the 10 th it erat iorijthe resid

11、ualis 0.008671,the lambdakis0.250000thel 1 th it erat iorijthe residualis 0. 005169,the lambdakis0.250000the 12 th it erat iorijthe residualis 0. 003082,the lambdakis0.250000the 13 th it erat iorijthe residualis 0. 001837,the lambdakis0.250000thel4 th it erat iorijthe residualis 0. 001095,the lambda

12、kis0.250000the 15 th it erat iorijthe residualis 0. 000653,the lambdakis0.250000the optimal solution iss_star =-0.4194-0.0001fm 二0.7729当 % =松 f (xk ) -1 时 程序如下：function x_star = newton1(x0,eps)g0 = grad(x0);res0 = norm(g0);res=res0;k = 0;xk =x0;while resepsk=k+1;H=inv(GRAND(xk);sk=-H*grad(xk);lambda

13、k = 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(x0);f0=f(x0);a=0;c1=0.1;c2=0.5;lambdak=1;d=-c1*lambdak*g0*s0;

14、xk=x0 + lambdak*s0;f1=f(xk);e=f0-f1;while delambdak=(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)A2+x(2)A2);function fa=f(x)fa=x(1)+2*x(2)A2+exp(x(1)A2+x(2)A2);functi

15、on G=GRAND(x)G=zeros(2,2);G(1)=2*exp(x(1)A2+x(2)A2)+4*x(1)A2*exp(x(1)A2+x(2)A2);G(2)=4*x(1)*x(2)*exp(x(1)A2+x(2)A2);G(3)=4*x(1)*x(2)*exp(x(1)A2+x(2)A2);G(4)=4+(2+4*x(2)A2)*exp(x(1)A2+x(2)A2);在命令窗口中输入 x0=0;1;eps=1e-3;newton1(x0,eps)输出结果为： z0= LO ; 1J ; eps= 1 e-3 ; nevrt onl (迥 eps)thel thiteration,t

16、heresidual is 3.650135,the lambdak is1.000000the2 thiteration,theresidual is 0.550720,the lambdak is1.000000the3 thit erat ion3theresidual is 0.009899the lambdak is1.000000the4 thit erat ioiijtheresidual is 0.000024the lambdak is1.000000the optimal solution isK_star =-0.41940.0000fm =0.7729取Hk为BFGS公

17、式时：程序如下：function x_star = bfgs1(x0,eps) g0 = grad(x0);res0 = norm(g0);res=res0;k = 0;xk =x0;while resepsk=k+1;if k=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

18、(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(xk);e

19、=f0-f1;while delambdak=(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)A2+x(2)A2);function fa=f(x)fa=x(1)+2*x(2)A2+exp(x(1)A2+x(2)A2);在命令窗口中输入 x0=0;1;eps=1e-3;bfgs1(x0,e

20、ps)输出结果为： kO= 0;1 ;eps=le-3;bfgsl (kOj eps)thel thiterat iorijtheresidual is1.320364,the lamb dakis0.125000the2 thiterat iorijtheresidual is0. 812021,the lamb dakis0.500000the3 thiterat iorijtheresidual is0.125024,the lamb dakis1.000000the4 thiterat iorijtheresidual is0.013504,the lamb dakis1.000000

21、the5 thiterat iorijtheresidual is0. 000282,the lamb dakis1.000000the optimal solution isK_star =-0.4194-0.0000 fm 二0.7729程序如下：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

22、(Ai(i,:)*xbi(i)+epsilon)index(i)=0;endendwhile (k0)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)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) in

23、dex(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(tm0)rb=Ae*ginvH*c + be;lambda=pinv(Ae*ginvH*Ae)*rb;x=ginvH*(Ae*lambda-c);elsex=-ginvH*c;lam

24、bda=0;endfunction fn=f(x)fn=(x(1)-1)A2+(x(2)-2.5)A2;在命令窗口中输入H=2 0;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)输出结果为: H=2 0;0 2;c= -2;-5;Ae= ; be=;Ai=l -2;-l -2;-l 2; 1 0;0 1 bi= -2;-6;-2;0 ;0;kO=O;O;qpact (Hj Cj Ae bej Aij bi, kO)z_star =1.4OOO1.

25、7000 fn 二0.SOOO利用乘子法求解下面的约束优化问题min 4xi 12s.t. 25 xf x| = 0程序如下程序如下:function x_star = lagrange(x0,epsn)10 xi x： + 10 x2 x项34 2 0 xi X2 2 v=1;u=0;m=h(u,x0);while mepsnif (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);fnfunction star = bfgs(v0,u0,x0)e

26、ps=1e-4;g0 = grad(v0,u0,x0);res0 = norm(g0);res=res0;k = 0;xk =x0;while resepsk=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=0.1;c2=0.5;lambdak=1;p=-

27、c1*lambdak*g0*s0;xk=x0 + lambdak*s0;f1=f(v0,u0,xk);q=f0-f1;while pqlambdak=(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)-x(2)人2- 34)*(10-2*x(1);g(

28、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)A2-12+v*(25-x(1)A2-x(2)A2)+2*(25-x(1)A2-x(2)A2)A2+0.125*(u+4*(10*x(1)-x(1)A2+10*x(2)-x(2)A2-34)A2-uA2);function h=h(u,x)if (10*x(1)-x(1)A2+10*x(2)-x(2)A2-34)(-u/4)h=(25-x

29、(1)A2-x(2)A2)A2+(10*x(1)-x(1)A2+10*x(2)-x(2)A2-34)A2;elseh=(25-x(1)A2-x(2)A2)A2+(-u/4)A2;endfunction d=d(x)d=10*x(1)-x(1)A2+10*x(2)-x(2)A2-34;function e=e(x)e=25-x(1)A2-x(2)A2;function fn=fn(x)fn=4*x(1)-x(2)A2-12;在命令窗 口中输入 x0=1;1;epsn=1e-4; lagrange(x0,epsn)输出结果为： zO=Ll:U ;epsn=le-4; lagrange(kO eps

30、n)s_tar =0.97804.9034fn 二-32.1319先写出求解二次规划子问题的程序程序如下：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=

31、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)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);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 (mk0&m0)dh1=dah(ep+rhoAi

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

33、=lam+alpha*dl; end if(m=0) ep=ep+alpha*de; mu=mu+alpha*du; end k=k+1;endfunction p=phi(ep,a,b)p=a+b-sqrt(aA2+bA2+2*epA2);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*

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

35、,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)A2+(gk(i)+Ai(i,:)*d)A2+2*epA2);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,mu,lam,

36、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), v1,Ae,dd2*Ai,zeros(l,l),zeros(m,l),zeros(l,m) ;dd1;endif(l=0)A=1,zeros(1,n),zeros(1,m);zeros(n,1),Bk

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

38、;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,inf);if(sigma*(tau+deta)1)sigma=sigma;elsesigma=1.0/(tau+

39、2*deta);endim=0; %Armijo 搜索while(im=20)if(phi1(x+rhoAim*dk,sigma)-phi1(x,sigma)0&m0)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*yk0.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精确价值函数%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