2023年优化方法上机大作业.doc

优化措施 上机大作业机械工程与材料能源学部能源与动力学院能源与环境工程 联络方式：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)|(g1*s0c2*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*deps if n3g=gfun(x0+d1*s0);d= double(solve(s0*g);x1=x0+d*s0;g1=gfun(x1);if norm(g1)=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 d=-g20g0; x1=x0+d;g1=gfun3_1(x1);if( norm(g1)=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)bi(i)+epsilon), index(i)=0; endendwhile (k0), 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)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) 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 = 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)运行成果： callqpactx = 1.4000 1.7000lambda = 0.8000exitflag = 0output = fval: -6.4500 iter: 7Function 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 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*ones(m,1);btak=10; btaold=10; while(btakepsilon & kepsilon 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=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;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; sigma1=0.4; epsilon1=1e-5; k=0; n=length(x0); Bk=eye(n); while(kmaxk) gk=feval(gfun,x0,varargin:); if(norm(gk)epsilon1), break; end dk=-Bkgk; m=0; mk=0; while(m20) newf=feval(fun,x0+rhom*dk,varargin:); oldf=feval(fun,x0,varargin:); if(newf0) 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=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(:,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,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; 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,inf); if(sigma*(tau+deta)1) sigma=sigma; else sigma=1.0/(tau+2*deta); end im=0; 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; yk=dlax(x1,mu,lam)-dlax(x,mu,lam); if(sk*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);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); 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(l00&m00) dp=df*d-1.0/sigma*(norm(h,1)+norm(gn,1);endfunction 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; endif(m0=0), l=f-mu*h; endif(l00&m00) 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(l10&m10), 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); 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(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; while (i0&m0) 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); 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; end if(m=0) ep=ep+alpha*de; d=d+alpha*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(l0&m0) 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,:)*d); end end if(m=0) dh(2:n+1)=Bk*d-Ae*mu+dfk; dh(n+2:n+l+1)=hk+Ae*d; end dh=dh(:); function bet=beta(ep,d,mu,lam,dfk,Bk,Ae,hk,Ai,gk,gamma) dh=dah(ep,d,mu,lam,dfk,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; end function A=JacobiH(ep,d,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