SQP算法

以下为两个m函数，分别是qpsubp和sqpm，其中qpsubp用于求解二次规划子问题，sqpm是用基于拉格朗日函数Hesse矩阵的SQP方法求解约束优化问题。具体的输入输出值意义已在程序中说明，具体是使用方法也在两段代码的末尾给出。将以下代码除去示例部分，分别保存为两个.m文件即可在matlab中使用下面给出两端代码供大家参考和交流：代码qpsub：function[d,mu,lam,val,k]=qpsubp(dfk,Bk,Ae,hk,Ai,gk)%功能:求解二次规划子问题%输入:dfk是xk处的梯度，Bk是第k次近似Hesse矩阵，Ae,hk线性等式约束有关的参数%Ai,gk是线性不等式约束的有关参数%输出:d,val分别是最优解和最优值,mu,lam是乘子向量,k是迭代次数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)epsilon)break;endA=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(l0&m0)de=dz(1);dd=dz(2:n+1);du=dz(n+2:n+l+1);dl=dz(n+l+2:n+l+m+1);endif(l==0)de=dz(1);dd=dz(2:n+1);dl=dz(n+2:n+m+1);endif(m==0)de=dz(1);dd=dz(2:n+1);du=dz(n+2:n+l+1);endi=0;while(i=20)if(l0&m0)dh1=dah(ep+rho^i*de,d+rho^i*dd,mu+rho^i*du,lam+rho^i*dl,dfk,Bk,Ae,hk,Ai,gk);endif(l==0)dh1=dah(ep+rho^i*de,d+rho^i*dd,mu,lam+rho^i*dl,dfk,Bk,Ae,hk,Ai,gk);endif(m==0)dh1=dah(ep+rho^i*de,d+rho^i*dd,mu+rho^i*du,lam,dfk,Bk,Ae,hk,Ai,gk);endif(norm(dh1)=(1-sigma*(1-gamma*ep0)*rho^i)*norm(dh))mk=i;break;endi=i+1;if(i==20)mk=10;endendalpha=rho^mk;if(l0&m0)ep=ep+alpha*de;d=d+alpha*dd;mu=mu+alpha*du;lam=lam+alpha*dl;endif(l==0)ep=ep+alpha*de;d=d+alpha*dd;lam=lam+alpha*dl;endif(m==0)ep=ep+alpha*de;d=d+alpha*dd;mu=mu+alpha*du;endk=k+1;endval=f1(d);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%functionp=phi(ep,a,b)p=a+b-sqrt(a^2+b^2+2*ep^2);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%functiondh=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(l0&m0)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);endendif(l==0)dh(2:n+1)=Bk*d-Ai'*lam+dfk;for(i=1:m)dh(n+1+i)=phi(ep,lam(i),gk(i)+Ai(i,:)*d);endendif(m==0)dh(2:n+1)=Bk*d-Ae'*mu+dfk;dh(n+2:n+l+1)=hk+Ae*d;enddh=dh(:);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%functionbet=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%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%functionA=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(l0&m0)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];endif(l==0)A=[1,zeros(1,n),zeros(1,m);zeros(n,1),Bk,-Ai';v1,dd2*Ai,dd1];endif(m==0)A=[1,zeros(1,n),zeros(1,l);zeros(n,1),Bk,-Ae';zeros(l,1),Ae,zeros(l,l)];end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%例解二次规划问题minf(x)=x1^2+x1x2+2x2^2-6x1-2x2-12x3s.t.x1+x2+x3-2=0x1-2x2+3=0x1,x2,x3=0%%%%%%%%%目标函数f(x)%%%%%%%%%%%%%%%%functionf=f1(x)%f1.mf=x(1)^2+x(1)*x(2)+2*x(2)^2-6*x(1)-2*x(2)-12*x(3);在matlab命令窗口依次输入下列命令:dfk=[-6-2-12]';Bk=[210;140;000];Ae=[111];hk=[-2]';Ai=[1-20;100;010;001];gk=[3000]';[d,mu,lam,val,k]=qpsubp(dfk,Bk,Ae,hk,Ai,gk)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%代码sqpmfunction[x,mu,lam,val,k]=sqpm(x0,mu0,lam0)%功能:用基于拉格朗日函数Hesse矩阵的SQP方法求解约束优化问题:%minf(x)s.t.h_i(x)=0,i=1,...,l%输入:x0是初始点，mu0,lam0是乘子向量的初始值%输出:x,mu分别是近似最优点和相应的乘子%val是最优值,mh是约束函数的模,k是迭代次数maxk=100;%最大迭代次数n=length(x0);l=length(mu0);m=length(lam0);rho=0.5;eta=0.1;B0=eye(n);x=x0;mu=mu0;lam=lam0;Bk=B0;sigma=0.8;epsilon1=1e-6;epsilon2=1e-5;[hk,gk]=cons(x);dfk=df1(x);[Ae,Ai]=dcons(x);Ak=[Ae;Ai];k=0;while(kmaxk)[dk,mu,lam]=qpsubp(dfk,Bk,Ae,hk,Ai,gk);%求解子问题mp1=norm(hk,1)+norm(max(-gk,0),1);if(norm(dk,1)epsilon1)&(mp1epsilon2)break;end%检验终止准则deta=0.05;%罚参数更新tau=max(norm(mu,inf),norm(lam,inf));if(sigma*(tau+deta)1)sigma=sigma;elsesigma=1.0/(tau+2*deta);endim=0;%Armijo搜索while(im=20)temp=eta*rho^im*dphi1(x,sigma,dk);if(phi1(x+rho^im*dk,sigma)-phi1(x,sigma)temp)mk=im;break;endim=im+1;if(im==20)mk=10;endendalpha=rho^mk;x1=x+alpha*dk;[hk,gk]=cons(x1);dfk=df1(x1);[Ae,Ai]=dcons(x1);Ak=[Ae;Ai];lamu=pinv(Ak)'*dfk;%计算最小二乘乘子if(l0&m0)mu=lamu(1:l);lam=lamu(l+1:l+m);endif(l==0)mu=[];lam=lamu;endif(m==0)mu=lamu;lam=[];endsk=alpha*dk;%更新矩阵Bkyk=dlax(x1,mu,lam)-dlax(x,mu,lam);if(sk'*yk0.2*sk'*Bk*sk)theta=1;elsetheta=0.8*sk'*Bk*sk/(sk'*Bk*sk-sk'*yk);endzk=theta*yk+(1-theta)*Bk*sk;Bk=Bk+zk*zk'/(sk'*zk)-(Bk*sk)*(Bk*sk)'/(sk'*Bk*sk);x=x1;k=k+1;endval=f1(x);%p=phi1(x,sigma)%dd=norm(dk)%%%%%%%%%l1精确价值函数%%%%%%%%%%%%functionp=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);endif(m0==0)p=f+1.0/sigma*norm(h,1);endif(l00&m00)p=f+1.0/sigma*(norm(h,1)+norm(gn,1));end%%%%%%%%%价值函数的方向导数%%%%%%%%functiondp=dphi1(x,sigma,d)df=df1(x);[h,g]=cons(x);gn=max(-g,0);l0=le