大连理工优化方法 增广拉格朗日方法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<maxk)
gk=feval(dmpsi,x0,fun,hf,gfun,dfun,dhf,dgf,mu,lambda,sigma);
if(norm(gk)<epsilon1)
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<oldf+sigma1*rho^m*gk'*dk)
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<maxk)
[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。