大连理工优化方法-增广拉格朗日方法MATLAB程序

上机大作业II

定义目标函数fun

function f=fun(x)

x1=x(1);

x2=x(2);

f=4*x1-x2^2-12;

定义目标函数梯度函数dfun

function f=dfun(x)

x2=x(2);

f=[4;-2*x2];

定义等式约束函数hf

function qua=hf(x)

qua=25-x(1)^2-x(2)^2;

定义等式约束函数梯度函数dhf

function qua=dhf(x)

qua=[-2*x(1);-2*x(2)];

定义不等式约束函数gfun

function inq=gfun(x)

inq=10*x(1)-x(1)^2+10*x(2)-x(2)^2-34;

定义不等式约束梯度数dgf

function inq=dgf(x)

inq=[10-2*x(1);10-2*x(2)];

定义增广拉格朗日函数mpsi

function psi=mpsi(x,fun,hf,gfun,dfun,dhf,dgf,mu,lambda,sigma) f=feval(fun,x);

he=feval(hf,x);

gi=feval(gfun,x);

l=length(he);

m=length(gi);

psi=f;

s1=0;

for i=1:l

psi=psi-he(i)*mu(i);

s1=s1+he(i)^2;

end

psi=psi+0.5*sigma*s1;

s2=0.0;

for i=1:m

s3=max(0.0, lambda(i) - sigma*gi(i));

s2=s2+s3^2-lambda(i)^2;

end

psi=psi+s2/(2.0*sigma);

定义增广拉格朗日函数梯度函数dmpsi

function dpsi=dmpsi(x,fun,hf,gfun,dfun,dhf,dgf,mu,lambda,sigma)

dpsi=feval(dfun,x);

he=feval(hf,x);

gi=feval(gfun,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);

end

for i=1:m

dpsi=dpsi+(sigma*gi(i)-lambda(i))*dgi(:,i);

end

定义BFGS法函数函数bfgs

function [x,val,k]=bfgs(mpsi,dmpsi,x0,fun,hf,gfun,dfun,dhf,dgf,mu,lambda,sigma) maxk=1000;

rho=0.5;

sigma1=0.4;

epsilon1=1e-4;

k=0;

n=length(x0);

Bk=eye(n);

while(k

gk=feval(dmpsi,x0,fun,hf,gfun,dfun,dhf,dgf,mu,lambda,sigma);

if(norm(gk)

break;

end

dk=-Bk\gk;

m=0;

mk=0;

while(m<20)

newf=feval(mpsi,x0+rho^m*dk,fun,hf,gfun,dfun,dhf,dgf,mu,lambda,sigma); oldf=feval(mpsi,x0,fun,hf,gfun,dfun,dhf,dgf,mu,lambda,sigma);

if(newf

mk=m;

break;

end

m=m+1;

end

x=x0+rho^mk*dk;

sk=x-x0;

yk=feval(dmpsi,x,fun,hf,gfun,dfun,dhf,dgf,mu,lambda,sigma)-gk;

if(yk'*sk>0)

Bk=Bk-((Bk*sk)*sk'*Bk)/(sk'*Bk*sk)+(yk*yk')/(yk'*sk);

end

k=k+1;

x0=x;

end

val=feval(mpsi,x0,fun,hf,gfun,dfun,dhf,dgf,mu,lambda,sigma);

定义增广拉格朗日乘子法函数multphr

function answer=multphr(fun,hf,gfun,dfun,dhf,dgf,x0)

maxk=5000;

sigma=2.0;

eta=2.0;

theta=0.8;

k=0;

ink=0;

epsilon=1e-4;

x=x0;

he=feval(hf,x);

gi=feval(gfun,x);

l=length(he);

m=length(gi);

mu=0.1*ones(l,1);

lambda=0.1*ones(m,1);

btak=10;

btaold=10;

while(btak>epsilon&&k

[x,v,ik]=bfgs('mpsi','dmpsi',x0,fun,hf,gfun,dfun,dhf,dgf,mu,lambda,sigma); ink=ink+ik;

he=feval(hf,x);

gi=feval(gfun,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+temp^2;

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;

end

f=feval(fun,x);

x

f

mu

lambda

k

运行求解

>> x0=[0;0]

x0 =

>> multphr('fun','hf','gfun','dfun','dhf','dgf',x0) x =

1.00128148956437

4.89871784708758

f =

-31.9923105871169

mu =

1.01559644571312

lambda =

0.754451167977228

k =

4