實用最優(yōu)化方法Matlab程序設計

1、實用最優(yōu)化方法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利用乘子法求解下面的約束優(yōu)化問題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先寫出求解二次規(guī)劃子問題的程序程序如下：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序列二次規(guī)劃程序為：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