《优化方法 上机大作业》由会员分享，可在线阅读，更多相关《优化方法 上机大作业（20页珍藏版）》请在金锄头文库上搜索。
1、优化方法 上机大作业能源与动力工程学院 2012 级 优化方法优化方法 上机大作业上机大作业机械工程与材料能源学部机械工程与材料能源学部能源与动力学院能源与动力学院能源与环境工程能源与环境工程联系方式：联系方式：15326178310优化方法 上机大作业能源与动力工程学院 2012 级 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-x
2、1(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-f1epsif n=eps s0=-g0; g=gfun3_1(x0+d1*s0); d= double(solve(s0*g); x1=x0+d*s0; g1=gfun3_1(x1); if( norm(g1)=epsd=-g20g0;x1=x0+d; g1=gfun3_1(x1); if( norm(g1)=eps s0=-h0*g0; g=gfun3_1(x0+d1*s0); d= double(solve(s
3、0*g); x1=x0+d*s0; g1=gfun3_1(x1); if( norm(g1)bi(i)+epsilon), index(i)=0; end end while (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;f
4、or(i=1:ni)if(index(i) break;endendendk=k+1;优化方法 上机大作业能源与动力工程学院 2012 级 else exitflag=1;alpha=1.0; tm=1.0;for(i=1:ni)if(index(i)=0)lambda = pinv(Ae*ginvH*Ae)*rb;x = ginvH*(Ae*lambda-c); elsex = -ginvH*c;lambda = zeros(m,1); end callqpact.m 文件 function callqpact H=2 0; 0 2; c=-2 -5; Ae= ; be= ; Ai=1 -2
5、; -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) 运行结果：callqpact x =1.4000优化方法 上机大作业能源与动力工程学院 2012 级 1.7000 lambda =0.8000 exitflag =0 output = fval: -6.4500iter: 7Function x,mu,lambda,output=multphr(fun,hf,gf,dfun,dhf,dgf,x0) function psi=mpsi(x,fu
6、n,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 dgi = dg1(x) dgi = 10-2*x(1), 10-2*x(2); function g=df1(x) g = 4, -2.0*x(2); f
7、unction 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; 优化方法 上机大作业能源与动力工程学院 2012 级 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*ones(m,1); btak=10; btaold=10; while(btakepsilon endfor
8、(i=1:l), mu(i)=mu(i)-sigma*he(i); endfor (i=1:m)lambda(i)=max(0.0,lambda(i)-sigma*gi(i);endendk=k+1;btaold=btak;x0=x; end f=feval(fun,x); output.fval=f; output.iter=k; output.inner_iter=ink; 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; f
9、or(i=1:l)psi=psi-he(i)*mu(i);s1=s1+he(i)2; end psi=psi+0.5*sigma*s1;优化方法 上机大作业能源与动力工程学院 2012 级 s2=0.0; for(i=1:m)s3=max(0.0, lambda(i) - sigma*gi(i);s2=s2+s32-lambda(i)2; end psi=psi+s2/(2.0*sigma); maxk=500; rho=0.55; sigma1=0.4; epsilon1=1e-5; k=0; n=length(x0); Bk=eye(n); while(k0)Bk=Bk-(Bk*sk*sk
10、*Bk)/(sk*Bk*sk)+(yk*yk)/(yk*sk);endk=k+1; x0=x; end val=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=feval(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
11、(:,i); end for(i=1:m)dpsi=dpsi+(sigma*gi(i)-lambda(i)*dgi(:,i); end优化方法 上机大作业能源与动力工程学院 2012 级 x0=1,1; x,mu,lambda,output=multphr(f1,h1,g1,df1,dh1, dg1,x0)x =1.00134.8987 mu =2.0312 lambda =0.7545 output = fval: -31.9923iter: 5inner_iter: 58bta: 4.3187e-07function x,mu,lam,val,k=sqpm(x0,mu0,lam0) max
12、k=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(k0 lam=lamu(l+1:l+m);endif(l=0), mu=; lam=lamu; endif(m=0), mu=lamu; lam=; endsk=a
13、lpha*dk; yk=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; end val=f1(x); function p=phi1(x,sigma) f=f1(x); h,g=cons(x); gn=max(-g,0); l0=length(h); m0=leng
14、th(g);优化方法 上机大作业能源与动力工程学院 2012 级 if(l0=0), p=f+1.0/sigma*norm(gn,1); end if(m0=0), p=f+1.0/sigma*norm(h,1); end if(l00 end function 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); end if(m0=0), dp=df*d-1.0/sigma*norm(h,1);
15、 end if(l00 end function l=la(x,mu,lam) f=f1(x); h,g=cons(x); l0=lemgth(h); m0=length(g); if(l0=0), l=f-lam*g; end if(m0=0), l=f-mu*h; end if(l00 end function 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; end if(m1=0), dl=df-Ae*mu; end if(l10 end function 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,
