优化方法 上机大作业

优化方法 上机大作业

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级