優(yōu)化方法上機大作業(yè)

1、優(yōu)化方法 上機大作業(yè)機械工程與材料能源學部能源與動力學院能源與環(huán)境工程 x0=0;1T;%初始值s0=-1;1T;%初始搜索方向c1=0.1;c2=0.5;a=0;b=inf;d=1;n=0;x1=x0+d*s0;g0=(x0(2)-x0(1)2)*x0(1)-2*(1-x0(1);(x0(2)-x0(1)2);g1=(x1(2)-x1(1)2)*x1(1)-2*(1-x1(1);(x1(2)-x1(1)2);f1=(x1(2)-x1(1)2)2+(1-x1(1)2;f0=(x0(2)-x0(1)2)2+(1-x0(1)2;while(f0-f1<-c1*d*g0'*s0)|(g

2、1'*s0<c2*g0'*s0)if (f0-f1)<(-c1*d*g0'*s0)b=d;d=(d+a)/2;x1=x0+d*s0;g0=(x0(2)-x0(1)2)*x0(1)-2*(1-x0(1);(x0(2)-x0(1)2);g1=(x1(2)-x1(1)2)*x1(1)-2*(1-x1(1);(x1(2)-x1(1)2);f1=(x1(2)-x1(1)2)2+(1-x1(1)2;f0=(x0(2)-x0(1)2)2+(1-x0(1)2;elseif (g1')*s0)<(c2*(g0')*s0)a=d;if(2*d<=(d

3、+b)/2)d=2*d;elsed=(d+b)/2;endx1=x0+d*s0;g0=(x0(2)-x0(1)2)*x0(1)-2*(1-x0(1);(x0(2)-x0(1)2);g1=(x1(2)-x1(1)2)*x1(1)-2*(1-x1(1);(x1(2)-x1(1)2);f1=(x1(2)-x1(1)2)2+(1-x1(1)2;f0=(x0(2)-x0(1)2)2+(1-x0(1)2;endendx1f1=(x1(2)-x1(1)2)2+(1-x1(1)2x1 = -0.0000 1.0000d = 1.1102e-016f1 = 2function f = fun( x )%UNTI

4、TLED3 Summary of this function goes here% Detailed explanation goes heref=x(1)2-2*x(1)*x(2)+2*x(2)2+x(3)2+x(4)2-x(2)*x(3)+2*x(1)+3*x(2)-x(3);Endfunction g = fun( x )%UNTITLED4 Summary of this function goes here% Detailed explanation goes here g=2 -2 0 0;-2 4 -1 0;0 -1 2 0;0 0 0 2*x+2;3;-1;0; endx0=0

5、;0;0;0; %初始值eps=1.0e-4; %精度g0=gfun(x0);s0=-g0;n=0;syms d1;while norm(g0)>eps if n<3g=gfun(x0+d1*s0);d= double(solve(s0'*g);x1=x0+d*s0;g1=gfun(x1);if norm(g1)<eps n=n+1; x0=x1; breakelses0=-g1+(norm(g1)2/norm(g0)2)*s0;x0=x1;g0=g1;endelseif n=3x0=x1;g0=gfun(x0);s0=-g0;n=0;endn=n+1;endx0nf

6、un(x0)x0 = -4 -3 -1 0n = 3ans =-8function f= fun3_1(x )%FUN3 Summary of this function goes here% Detailed explanation goes heref=x(1)+2*x(2)2+exp(x(1)2+x(2)2);endfunction g= gfun3_1(x)%GFUN3_1 Summary of this function goes here% Detailed explanation goes here g=1+2*x(1)*exp(x(1)2+x(2)2);4*x(2)+2*x(2

7、)*exp(x(1)2+x(2)2);end(1) 最速下降法x0=0;1;%初始值eps=1.0e-5;%精度n=0;g0=gfun3_1(x0);syms d1;while norm(g0)>=epss0=-g0;g=gfun3_1(x0+d1*s0);d= double(solve(s0'*g);x1=x0+d*s0;g1=gfun3_1(x1);if( norm(g1)<eps) n=n+1; x0=x1;break;else x0=x1;g0=gfun3_1(x0);endn=n+1;endf0=fun3_1(x0)x0nf0 = 0.7729x0 = -0.41

8、94 0.0000n = 6（2）牛頓法function g2 = hesse(x)%HESSE Summary of this function goes here% Detailed explanation goes here g2=2*exp(x(1)2+x(2)2)+4*x(1)2*exp(x(1)2+x(2)2),4*x(1)*x(2)*exp(x(1)2+x(2)2) 4*x(1)*x(2)*exp(x(1)2+x(2)2),4+2*exp(x(1)2+x(2)2)+4*x(2)2*exp(x(1)2+x(2)2); endx0=0;1;%初始值eps=1.0e-5;%精度n=0;

9、g0=gfun3_1(x0);g20=hesse(x0);while norm(g0)>=eps d=-g20g0; x1=x0+d;g1=gfun3_1(x1);if( norm(g1)<eps) n=n+1; x0=x1;break;else x0=x1;g0=gfun3_1(x0);endn=n+1;endf0=fun3_1(x0)x0nf0 = 0.7729x0 = -0.4194 0.0000n =35（3）BFGSx0=0;1;%初始值eps=1.0e-5;%精度n=0;g0=gfun3_1(x0);syms d1;h0=eye(2);while norm(g0)>

10、;=epss0=-h0*g0;g=gfun3_1(x0+d1*s0);d= double(solve(s0'*g);x1=x0+d*s0;g1=gfun3_1(x1);if( norm(g1)<eps) n=n+1; x0=x1;break;else h0=h0-(h0*(g1-g0)*(g1-g0)'*h0)/(g1-g0)'*h0*(g1-g0)+. (x1-x0)*(x1-x0)')/(x1-x0)'*(g1-g0)+. (g1-g0)'*h0*(g1-g0)*(x1-x0)*(x1-x0)');x0=x1;g0=gfun3_

11、1(x0);endn=n+1;endf0=fun3_1(x0)x0nf0 = 0.7729x0 = -0.4194 0.0000n = 4functionx,lamk,exitflag,output=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), i

12、ndex(i)=0; endendwhile (k<=kmax) Aee= ; if(ne>0), Aee=Ae; end for(j=1:ni) if(index(j)>0), Aee=Aee; Ai(j,:); end end gk=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); end if(y>=0) exitflag=0

13、; 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) if(index(i)=0)&(Ai(i,:)*dk<0) tm1=(bi(i)-Ai(i,:)*x)/(Ai(i,:)*dk); if(tm1<tm) tm=tm1; ti=i; end end end alpha=min(alpha,tm); x = x+al

14、pha*dk; if(tm<1), index(ti)=1; end end if(exitflag=0), break; end k=k+1;endoutput.fval=0.5*x'*H*x+c'*x;output.iter=k; 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);else x = -g

15、invH*c; lambda = zeros(m,1);endcallqpact.m文件function callqpact 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'x, lambda, exitflag,output=qpact(H,c,Ae,be,Ai,bi,x0)運行結(jié)果： callqpactx = 1.4000 1.7000lambda = 0.8000exitflag = 0output = fval: -6.4500 iter: 7

16、Function x,mu,lambda,output=multphr(fun,hf,gf,dfun,dhf,dgf,x0)function psi=mpsi(x,fun,hf,gf,dfun,dhf,dgf,mu,lambda,sigma)function he=h1(x)he=-x(1)2-x(2)2+25.0;function f=f1(x)f=4*x(1)-x(2)2-12;function gi=g1(x)gi=10*x(1)-x(1)2+10*x(2)-x(2)2-34;function dhe = dh1(x)dhe = -1*x(1), -1*x(2)'function

17、 dgi = dg1(x)dgi = 10-2*x(1), 10-2*x(2)'function g=df1(x)g = 4, -2.0*x(2)'function x,val,k=bfgs(fun,gfun,x0,varargin)maxk=500; sigma=2.0; eta=2.0; theta=0.8; k=0; ink=0; epsilon=1e-5; x=x0; he=feval(hf,x); gi=feval(gf,x);n=length(x); l=length(he); m=length(gi);mu=0.1*ones(l,1); lambda=0.1*on

18、es(m,1);btak=10; btaold=10; while(btak>epsilon & k<maxk) x,ival,ik=bfgs('mpsi','dmpsi',x0,fun,hf,gf,dfun,dhf,dgf,mu,lambda,sigma); ink=ink+ik; he=feval(hf,x); gi=feval(gf,x); btak=0.0;for (i=1:l), btak=btak+he(i)2; end for i=1:m temp=min(gi(i),lambda(i)/sigma); btak=btak+te

19、mp2; end btak=sqrt(btak); if btak>epsilon if(k>=2 & btak > theta*btaold) sigma=eta*sigma; end for (i=1:l), mu(i)=mu(i)-sigma*he(i); end for (i=1:m) lambda(i)=max(0.0,lambda(i)-sigma*gi(i); end end k=k+1; btaold=btak; x0=x;endf=feval(fun,x);output.fval=f;output.iter=k;output.inner_iter=i

20、nk;output.bta=btak;f=feval(fun,x); he=feval(hf,x); gi=feval(gf,x);l=length(he); m=length(gi);psi=f; s1=0.0;for(i=1:l) psi=psi-he(i)*mu(i); s1=s1+he(i)2;endpsi=psi+0.5*sigma*s1;s2=0.0;for(i=1:m) s3=max(0.0, lambda(i) - sigma*gi(i); s2=s2+s32-lambda(i)2;endpsi=psi+s2/(2.0*sigma);maxk=500; rho=0.55; si

21、gma1=0.4; epsilon1=1e-5; k=0; n=length(x0); Bk=eye(n); while(k<maxk) gk=feval(gfun,x0,varargin:); if(norm(gk)<epsilon1), break; end dk=-Bkgk; m=0; mk=0; while(m<20) newf=feval(fun,x0+rhom*dk,varargin:); oldf=feval(fun,x0,varargin:); if(newf<oldf+sigma1*rhom*gk'*dk) mk=m; break; end m

22、=m+1; end x=x0+rhomk*dk; sk=x-x0; yk=feval(gfun,x,varargin:)-gk; if(yk'*sk>0) Bk=Bk-(Bk*sk*sk'*Bk)/(sk'*Bk*sk)+(yk*yk')/(yk'*sk); end k=k+1; x0=x;endval=feval(fun,x0,varargin:);function dpsi=dmpsi(x,fun,hf,gf,dfun,dhf,dgf,mu,lambda,sigma)dpsi=feval(dfun,x);he=feval(hf,x); gi=f

23、eval(gf,x);dhe=feval(dhf,x); dgi=feval(dgf,x);l=length(he); m=length(gi);for(i=1:l) dpsi=dpsi+(sigma*he(i)-mu(i)*dhe(:,i);endfor(i=1:m) dpsi=dpsi+(sigma*gi(i)-lambda(i)*dgi(:,i);end>> x0=1,1'x,mu,lambda,output=multphr('f1','h1','g1','df1','dh1','

24、dg1',x0)x = 1.0013 4.8987mu = 2.0312lambda = 0.7545output = fval: -31.9923 iter: 5 inner_iter: 58 bta: 4.3187e-07function 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; eps

25、ilon2=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=si

26、gma; else sigma=1.0/(tau+2*deta); end im=0; while(im<=20) if(phi1(x+rhoim*dk,sigma)-phi1(x,sigma)<eta*rhoim*dphi1(x,sigma,dk) mk=im; break; end im=im+1; if(im=20), mk=10; end end alpha=rhomk; x1=x+alpha*dk; hk,gk=cons(x1); dfk=df1(x1); Ae,Ai=dcons(x1); Ak=Ae; Ai; lamu=pinv(Ak)'*dfk; if(l&g

27、t;0&m>0) 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; yk=dlax(x1,mu,lam)-dlax(x,mu,lam); if(sk'*yk>0.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

28、*zk'/(sk'*zk)-(Bk*sk)*(Bk*sk)'/(sk'*Bk*sk); x=x1; k=k+1;endval=f1(x);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(l0>0&m0>0) p=f+1.0/sigma*(norm(h,1)+nor

29、m(gn,1); endfunction dp=dphi1(x,sigma,d)df=df1(x); h,g=cons(x); gn=max(-g,0);l0=length(h); m0=length(g);if(l0=0), dp=df'*d-1.0/sigma*norm(gn,1); endif(m0=0), dp=df'*d-1.0/sigma*norm(h,1); endif(l0>0&m0>0) dp=df'*d-1.0/sigma*(norm(h,1)+norm(gn,1);endfunction l=la(x,mu,lam)f=f1(x

30、); h,g=cons(x);l0=lemgth(h); m0=length(g);if(l0=0), l=f-lam*g; endif(m0=0), l=f-mu'*h; endif(l0>0&m0>0) l=f-mu'*h-lam'*g; endfunction dl=dlax(x,mu,lam)df=df1(x); Ae,Ai=dcons(x); m1,m2=size(Ai); l1,l2=size(Ae);if(l1=0), dl=df-Ai'*lam; endif(m1=0), dl=df-Ae'*mu; endif(l1&

31、gt;0&m1>0), dl=df-Ae'*mu-Ai'*lam; endfunction f=f1(x)f=4*x(1)-x(2)2-12;function df=df1(x)df=4, -2*x(2)'function h,g=cons(x)h=25-x(1)2-x(2)2;g=x(1);x(2);function dh,dg=dcons(x)dh=-2*x(1), -2*x(2);dg=1 0; 0 1;Qpsubp.m文件function d,mu,lam,val,k=qpsubp(dfk,Bk,Ae,hk,Ai,gk) n=length(dfk)

32、; 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(

33、dh)<epsilon) break; end A=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=Ab; 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); 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

34、+1); du=dz(n+2:n+l+1); end i=0; while (i<=20) if(l>0&m>0) dh1=dah(ep+rhoi*de,d+rhoi*dd,mu+rhoi*du,lam+rhoi*dl,dfk,Bk,Ae,hk,Ai,gk); 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); en

35、d if(norm(dh1)<=(1-sigma*(1-gamma*ep0)*rhoi)*norm(dh) mk=i; break; end i=i+1; if(i=20), mk=10; end end alpha=rhomk; if(l>0&m>0) 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; end if(m=0) ep=ep+alpha*de; d=d+alp

36、ha*dd; mu=mu+alpha*du; end k=k+1; end val=0.5*d'*Bk*d+dfk'*d; function 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(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); end end if(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,:)