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

上机大作业II定义目标函数funfunctionf=fun(x)x1=x(1);x2=x(2);f=4*x1-x2^2-12;定义目标函数梯度函数dfunfunctionf=dfun(x)x2=x(2);f=[4;-2*x2];定义等式约束函数hffunctionqua=hf(x)qua=25-x(1)^2-x(2)^2;定义等式约束函数梯度函数dhffunctionqua=dhf(x)qua=[-2*x(1);-2*x(2)];定义不等式约束函数gfunfunctioninq=gfun(x)inq=10*x(1)-x(1)^2+10*x(2)-x(2)^2-34;定义不等式约束梯度数dgffunctioninq=dgf(x)inq=[10-2*x(1);10-2*x(2)];定义增广拉格朗日函数mpsifunctionpsi=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;fori=1:lpsi=psi-he(i)*mu(i);s1=s1+he(i)^2;endpsi=psi+0.5*sigma*s1;s2=0.0;fori=1:ms3=max(0.0,lambda(i)-sigma*gi(i));s2=s2+s3^2-lambda(i)^2;endpsi=psi+s2/(2.0*sigma);定义增广拉格朗日函数梯度函数dmpsifunctiondpsi=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);fori=1:ldpsi=dpsi+(sigma*he(i)-mu(i))*dhe(:,i);endfori=1:mdpsi=dpsi+(sigma*gi(i)-lambda(i))*dgi(:,i);end定义BFGS法函数函数bfgsfunction[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(kmaxk)gk=feval(dmpsi,x0,fun,hf,gfun,dfun,dhf,dgf,mu,lambda,sigma);if(norm(gk)epsilon1)break;enddk=-Bk\gk;m=0;mk=0;while(m20)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(newfoldf+sigma1*rho^m*gk'*dk)mk=m;break;endm=m+1;endx=x0+rho^mk*dk;sk=x-x0;yk=feval(dmpsi,x,fun,hf,gfun,dfun,dhf,dgf,mu,lambda,sigma)-gk;if(yk'*sk0)Bk=Bk-((Bk*sk)*sk'*Bk)/(sk'*Bk*sk)+(yk*yk')/(yk'*sk);endk=k+1;x0=x;endval=feval(mpsi,x0,fun,hf,gfun,dfun,dhf,dgf,mu,lambda,sigma);定义增广拉格朗日乘子法函数multphrfunctionanswer=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(btakepsilon&&kmaxk)[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;fori=1:lbtak=btak+he(i)^2;endfori=1:mtemp=min(gi(i),lambda(i)/sigma);btak=btak+temp^2;endbtak=sqrt(btak);ifbtakepsilonif(k=2&&btaktheta*btaold)sigma=eta*sigma;endfori=1:lmu(i)=mu(i)-sigma*he(i);endfori=1:mlambda(i)=max(0.0,lambda(i)-sigma*gi(i));endendk=k+1;btaold=btak;x0=x;endf=feval(fun,x);xfmulambdak运行求解x0=[0;0]x0=00multphr('fun','hf','gfun','dfun','dhf','dgf',x0)x=1.001281489564374.89871784708758f=-31.9923105871169mu=1.01559644571312lambda=0.754451167977228k=4