优化方法 上机大作业
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; end
val=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); end if(m0==0), p=f+1.0/sigma*norm(h,1); end if(l0>0&m0>0)
p=f+1.0/sigma*(norm(h,1)+norm(gn,1)); 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); end if(l0>0&m0>0)
dp=df'*d-1.0/sigma*(norm(h,1)+norm(gn,1)); 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(l0>0&m0>0)
l=f-mu'*h-lam'*g; 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(l1>0&m1>0), dl=df-Ae'*mu-Ai'*lam; 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)
能源与动力工程学院 2012级
优化方法 上机大作业
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) 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=A\\b; 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+1); du=dz(n+2:n+l+1); end i=0; while (i<=20) if(l>0&m>0) dh1=dah(ep+rho^i*de,d+rho^i*dd,mu+rho^i*du,lam+rho^i*dl,dfk,Bk,Ae,hk,Ai,gk); end if(l==0) dh1=dah(ep+rho^i*de,d+rho^i*dd,mu,lam+rho^i*dl,dfk,Bk,Ae,hk,Ai,gk); end 能源与动力工程学院 2012级 优化方法 上机大作业 if(m==0) dh1=dah(ep+rho^i*de,d+rho^i*dd,mu+rho^i*du,lam,dfk,Bk,Ae,hk,Ai,gk); end if(norm(dh1)<=(1-sigma*(1-gamma*ep0)*rho^i)*norm(dh)) mk=i; break; end i=i+1; if(i==20), mk=10; end end alpha=rho^mk; 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+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(a^2+b^2+2*ep^2); 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) 能源与动力工程学院 2012级 优化方法 上机大作业 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*ep^2); 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(l>0&m>0) A=[1, zeros(1,n), zeros(1,l), zeros(1,m); zeros(n,1), Bk, -Ae', - Ai'; zeros(l,1), Ae, zeros(l,l), zeros(l,m) ; v1, dd2*Ai, zeros(m,l), dd1]; end if(l==0) A=[1, zeros(1,n), zeros(1,m); zeros(n,1), Bk, -Ai'; v1, dd2*Ai, dd1]; end if(m==0) A=[1, zeros(1,n), zeros(1,l); zeros(n,1), Bk, -Ae'; zeros(l,1), Ae, zeros(l,l), ]; 能源与动力工程学院 2012级 优化方法 上机大作业 end > x0=[3 3]'; mu0=[0]'; lam0=[0 0]'; [x,mu,lam,val,k]=sqpm(x0,mu0,lam0) x = -0.0000 5.0000 mu = 1.0000 lam = 4.0000 -0.0000 val = -37.0000 k = 5 c=[3;1;1]; A=[2 1 1;1 -1 -1]; b=[2;-1]; vlb=[0;0;0]; [x,fval]=linprog(c,A,b,[],[],vlb) Optimization terminated. x = 0.0000 0.5000 0.5000 fval = 能源与动力工程学院 2012级
优化方法 上机大作业



