序列二次规划法.pdf
- 文档编号:14650451
- 上传时间:2023-06-25
- 格式:PDF
- 页数:24
- 大小:557.82KB
序列二次规划法.pdf
《序列二次规划法.pdf》由会员分享,可在线阅读,更多相关《序列二次规划法.pdf(24页珍藏版)》请在冰点文库上搜索。
序列二次规划算法序列二次规划算法1序列二次规划法简介序列二次规划法简介非线性规划问题是目标函数或约束条件中包含非线性函数的规划问题。
一般说来,解非线性规划要比解线性规划问题困难得多。
而且,也不像线性规划有单纯形法这一通用方法,非线性规划目前还没有适于各种问题的一般算法,各个方法都有自己特定的适用范围。
当目标函数和约束条件具有良好的分析性质时,人们喜欢用间接法分析与求解约束最优化问题。
利用间接法求解最优化问题的途径一般有两种,另一种途径是在可行域内使目标函数下降的迭代点法,如可行点法。
另一种是利用目标函数和约束条件构造增广目标函数,借此将约束最优化问题转化为无约束最优化问题,然后利用求解无约束最优化问题的方法间接求解新目标函数的局部最优解或稳定点,如人们所熟悉的惩罚函数法,乘子法和序列二次规划等。
序列二次规划算法是目前公认的求解约束非线性优化问题的最有效方法之一,与其它优化算法相比,其最突出的优点是收敛性好、计算效率高、边界搜索能力强,受到了广泛重视及应用。
但其迭代过程中的每一步都需要求解一个或多个二次规划子问题。
一般地,由于二次规划子问题的求解难以利用原问题的稀疏性、对称性等良好特性,随着问题规模的扩大,其计算工作量和所需存贮量是非常大的。
因此,目前的序列二次规划方法一般只适用于中小型问题。
另外,由于大型二次规划问题的求解通常使用迭代法,所需解的精度越高,花费的时间就越多,稳定性也就越差,相对线性方程组的求解理论来说,二次规划的求解是不完善的。
2序列二次规划的研究序列二次规划的研究最优化理论及方法是一个具有广泛应用背景的研究领域。
它研究诸如从众多的方案中选出最优方案等问题,常见的各种模型如线性规划,二次规划,非线性规划,多目标规划等。
最优化理论及方法已经在经济计划,工程设计,生产管理,交通运输等领域得到了广泛的应用,吸引了很多学者的研究。
一般非线性约束的数学规划问题是:
min().()0(1,2,.,)()0(1,2,.,)uvfXstgXuphXvm其中nXR是决策变量,()fX是目标函数,()ugX,()vhX分别是不等式约束函数和等式约束函数。
在求解上述问题的诸多算法中序列二次规划(SQP)是最有效的一类方法之一由于它具有超线性收敛速度,对它的研究一直是非线性规划的一个热点,许多学者对它进行了研究和改进序列二次规划(SQP)算法最早由Wilson(1963)提出。
Wilson提出了Newton-SQP方法。
60年代末和70年代初,随着解无约束优化问题的拟Newton法的发展,拟Newton-SQP算法的研究引起了学者们的高度重视。
Polomares-Mangasarian(1976)提出了一类拟Newton-SQP方法。
在该类方法中,他采用对Lagrange函数的整体Hesse矩阵(即既对x又对Lagrange乘子的Hesse矩阵)的近似进行修正的方法来修正。
Han(1976,1977)证明了PSB-SQP和BFGS-SQP方法的局部收敛性。
而且,若采用精确罚函数进行搜索,则算法具有全局收敛性。
自此,SQP算法的研究受到了广泛的重视。
迄今为止,SQP算法的研究工作已取得了丰硕的成果。
SQP算法的主要优点之一是它的全局收敛性和局部超线性收敛性。
为使算法具有全局收敛性,通常要求子问题(QP)中的矩阵H对称正定。
使得(QP)产生的方向S是某个效益函数的下降方向。
另一方面,若矩阵H对称正定,则子问题(QP)是一个严格凸二次规划。
此时,该问题有唯一解。
而且,其求解也相对容易。
因此,研究矩阵H的对称正定性是一项非常重要的工作。
并己引起了学者们的重视。
Powell(1976)提出了一个修正的BFGS公式。
该公式可保证产生的矩阵是对称正定的。
数值计算的结果表明利用该修正公式的SQP算法的数值效果较好。
然而该算法的局部收敛性质尚不清楚。
Coleman和Conn(1984)提出了一种既约Hessian-SQP方法。
即用一系列的正定矩阵去逼近解*,*S处的既约Hessian矩阵。
Noeedal和Overton(1985)提出了单边及双边的投影Hessian方法。
在单边投影Hessian方法中证明了局部1-步超线性收敛性而在双边投影Hessian方法中证明了在某种假设条件下的2-步超线性收敛性质。
但他的BFGS修正公式是在满足一定的修正准则才进行修正。
Xie和Byrd(1999)在RHSQP方法的基础上加上线性搜索得到全局收敛性算法。
他并且把Nocedal和Overton的修正法则改为了一种新的修正法则一正曲率法则。
这样能得到更好的数值结果。
并且证明了当利用两个常用的效益函数,Flether效益函数和1l精确罚函数进行线性搜索时,算法具有全局收敛性和R-线性收敛性。
R.Fontecilla(1988)提出了一种2-步超线性收敛算法。
与Coleman和Conn不同之处在于他将步长犷分为互相垂直的两个向量:
即水平方向的步长kh和垂直方向的步长kv。
且仅利用水平方向的步长kh对H进行修正。
与coleman和Conn不同之处在于对kS进行分解时他使用了一个正交投影算子而不是标准正交基,比较容易满足连续性条件。
比前者适用的范围更广。
而且,他把这一方法用在了几种不同的乘子修正公式中。
在二阶充分条件下,证明了该算法是2-步超线性收敛的。
若在每一次迭代中再计算一次约束函数的值,则算法为1-步超线性收敛。
但没有全局收敛性研究。
3序列二次规划算法推导过程序列二次规划算法推导过程序列二次规划(SQP)算法是将复杂的非线性约束最优化问题转化为比较简单的二次规划(QP)问题求解的算法。
所谓二次规划问题就是目标函数为二次函数,约束函数为线性函数的最优化问题。
二次规划问题是最简单的非线性约束最优化问题。
3.1序列二次规划算法思想序列二次规划算法思想非线性约束最优化问题:
min().()0(1,2,.,)()0(1,2,.,)uvfXstgXuphXvm(1-1)利用泰勒展开把非线性约束问题式(1-1)的目标函数在迭代点Xk简化成二次函数,把约束函数简化成线性函数后得到的就是如下的二次规划问题:
T2TTT1min()()()2.()()0(1,2,.,)()()0(1,2,.,)kkkkkkkkuukkkvvfXXXfXXXfXXXstgXXXgXuphXXXhXvm(1-2)此问题是原约束最优化问题的近似问题,但其解不一定是原问题的可行点。
为此,令kSXX将上述二次规划问题变成关于变量的S的问题,即T2TTT1min()()()2.()()0(1,2,.,)()()0(1,2,.,)kkkkuukkvvfXSfXSfXSstgXSgXuphXShXvm(1-3)令2Teq12T12Teq12T12()()(),(),.,()(),(),.,()(),(),.,()(),(),.,()kkkkkmkkkpkkkmkkkpHfXCfXAhXhXhXAgXgXgXBhXhXhXBgXgXgX(1-4)将式(1-4)变成二次规划问题的一般形式,即TTeqeq1min2.SHSCSstASBASB(1-5)求解此二次规划问题,将其最优解*S作为原问题的下一个搜索方向kS,并在该方向上进行原约束问题目标函数的约束一维搜索,就可以得到原约束问题的一个近似解1kX。
反复这一过程,就可以求得原问题的最优解。
上述思想得以实现的关键在于如何计算函数的二阶导数矩阵H,如何求解式(1-5)所示的二次规划问题。
3.2二阶导数矩阵的计算二阶导数矩阵的计算二阶导数矩阵的近似计算可以利用拟牛顿(变尺度)法中变尺度矩阵计算的DFP公式TT1TTkkkkkkkkkkkkkXXHqqHHHqXqHq(1-6)或BFGS公式1TTTTTTT1kkkkkkkkkkkkkkkkkkkHHXXXqXXqHqHqqXqHXq(1-7)3.3二次规划问题的求解二次规划问题的求解二次规划问题式(1-7)的求解分为以下两种情况。
(1)等式约束二次规划问题TTeqeq1min()2.fXSHSCSstASB(1-8)其拉格朗日函数为TTTeqeq1min(S,)()2LSHSCSASB由多元函数的极值条件(S,)0L得Teqeqeq00HSCAASB写成矩阵形式,即Teqeqeq0CHASBA(1-9)式(1-9)其实就是以T,S为变量的线性方程组,而且变量数和方程数都等于nm。
由线性代数知,此方程要么无解,要么有惟一解。
如果有解,利用消元变换可以方便的地求出该方程的惟一解,记作11T,kkS。
根据k-t条件,若此解中的乘子向量1k不全为零,则1kS就是等式约束二次规划问题式(1-8)的最优解*S,即1*kSS。
(2)一般约束二次规划问题对于一般约束下的二次规划问题式(1-5),在不等式约束条件中找出迭代点kX的起作用的约束,将等式约束和起作用的约束组成新的约束条件,构成新的等式约束问题:
TT11min()2.=knijjjiEIjfXSHSCSstasb(1-10)其中,E代表等式约束下的集合,kI代表不等式约束中起作用约束的下标集合。
此式即式(1-8),可以用同样的方法求解。
在求得式(1-10)的解11T,kkS之后,根据k-t条件,若解中对应原等式约束条件的乘子不全为零,对应起作用约束条件的乘子不小于零,则1kS就是所求一般约束二次规划问题式(1-5)的最优解*S。
综上所述,在迭代点kX上先进行矩阵kH的变更,在构造和求解相应的二次规划子问题,并该子问题最优解*S作为下一次迭代的搜索方向kS。
然后在该方向上对原非线性最优化问题目标函数进行约束一维搜索,得到下一个迭代点1kX,并判断收敛精度是否满足。
重复上述过程,直到迭代点1kX最终满足终止准则,得到原非线性最约束问题的最优解*X为止。
这种算法称为二次规划法,它是目前求解非线性约束最优化问题的常用方法,简称SQP法。
序列二次规划算法的迭代步奏如下:
1给定初始点0X、收敛精度,令0HI(单位矩阵),置0k.2在点kX简化原问题为二次规划问题式(1-10)。
3求解二次规划问题,并令*kSS。
4在方向kS上对原问题目标函数进行约束一维搜索,得点1kX。
5终止判断,若1kX满足给定的精度的终止准则,则令1*kXX,1*()kffX,输出最优解,终止计算,否则转6。
6按式(1-6)或(1-7)修正1kH,令1kk,转2继续迭代。
4MATLAB程序程序%*%*functiond,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是迭代次数%功能:
求解二次规划子问题%输入:
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;%k为迭代次数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)0&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;%mk=0;mm=0;while(mm0&m0)dh1=dah(ep+rhoi*de,d+rhoi*dd,mu+rhoi*du,lam+rhoi*dl,dfk,Bk,Ae,hk,Ai,gk);endif(l=0)dh1=dah(ep+rhoi*de,d+rhoi*dd,mu,lam+rhoi*dl,dfk,Bk,Ae,hk,Ai,gk);endif(m=0)dh1=dah(ep+rhoi*de,d+rhoi*dd,mu+rhoi*du,lam,dfk,Bk,Ae,hk,Ai,gk);endif(norm(dh1)0&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;end%*functionp=phi(ep,a,b)p=a+b-sqrt(a2+b2+2*ep2);%*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;fori=1:
mdh(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;fori=1:
mdh(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);%*functiondd1,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*ep2);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%*%*functionx,mu,lam,val,k=sqpm(x0,mu0,lam0)%功能:
用基于拉格朗日函数Hesse阵的SQP方法求解约束优化问题:
%minf(x)s.t.h_i(x)=0,i=1,.,l.%输入:
x0是初始点,mu0是乘子向量的初始值%输出:
x,mu分别是近似最优点及相应的乘子%val是最优值,mh是约束函数的模,k是迭代次数。
maxk=1000;%最大迭代次数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)if(phi1(x+rhoim*dk,sigma)-phi1(x,sigma)0&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);%*%*%*精确价值函数*%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=length(h);m0=length(g);if(l0=0),dp=df*d-1.0/sigma*norm(gn,1);endif(m0=0),dp=df*d-1.0/sigma*norm(h,1);endif(l00&m00)dp=df*d-1.0/sigma*(norm(h,1)+norm(gn,1);end%*%*拉格朗日函数*%functionl=la(x,mu,lam)f=f1(x);%调用目标函数文件h,g=cons(x);%调用约束函数文件l0=lemgth(h);m0=length(g);if(l0=0),l=f-lam*g;endif(m0=0),l=f-mu*h;endif(l0.0&m0.0)l=f-mu*h-lam*g;end%*%*拉格朗日函数的梯度*%functiondl=dlax(x,mu,lam)df=df1(x);%调用目标函数梯度文件Ae,Ai=dcons(x);%调用约束函数Jacobi矩阵文件m1,m2=size(Ai);l1,l2=size(Ae);if(l1=0),dl=df-Ai*lam;endif(m1=0),dl=df-Ae*mu;endif(l10&m10),dl=df-Ae*mu-Ai*lam;end%*functionf=f1(x)f=;%(目标函数)%*functiondf=df1(x)df=;%目标函数关于各个变量的导数,中间用“;”隔开%*functionh,g=cons(x)h=;%输入等式约束函数,中间用“;”隔开g=;%输入不等式约束函数,中间用“;”隔开%*functiondh,dg=dcons(x)dg=;%等式约束函数关于各个变量的导数,中间用“,”隔开dh=;%不等式约束函数关于各个变量的导数,中间用“,”隔开%*%在MATLAB的命令窗口输入以下命令x0=;%各个变量的初值,中间用“空格”隔开;mu0=;%等式约束函数的个数,有几个就有几个“0”,中间用“空格”隔开lam0=;%不等式约束函数的个数,有几个就有几个“0”,中间用“空格”隔开x,mu,k=sqpm(x0,mu0,lam0)%*%*5实例检验和计算实例检验和计算对于上述程序,我们首先应该验证其正确性,才能用来解决问题。
例1:
根据例1的要求,改动相应的f1.m,df1.m,cons.m,dcons.m文件,在命令窗口输入相应的命令,结果如下图所示:
对于上述结果,我们可以用相应的MATLAB的优化函数来解决,看所得结果是否相同。
编写MATLAB目标函数文件functionc,ce=cxmcon(x)ce=pi*x
(1)*x
(2)+pi*x
(1)2-150;c=;非线性约束函数返回变量分别是c和ce两个量,其中,前者为不等式约束的数学描述,后者为非线性等式约束,如果某个约束不存在,则应该将其赋值为空矩阵。
观察可知,例1只有非线性等式约束,没有非线性不等式约束,则A,B,Aeq,Beq,都将为空矩阵了。
另外,应该给出搜索的初值,向量x0=3,3,因此,可以调用fmincon()函数求解此约束最优化问题。
在命令窗口输入以下内容:
f=(x)-pi*x
(1)2*x
(2);ff=optimset;ff.LargeScale=off;ff.TolX=1e-15;ff.TolFun=1e-20;ff.TolCon=1e-20;x0=1;1;xm=0;0;xM=;A=;B=;Aeq=;Beq=;x,f_opt,c,d=fmincon(f,x0,A,B,Aeq,Beq,xm,xM,cxmcon)程序运行后的优化结果:
例2:
根据例2的要求,改动相应的f1.m,df1.m,cons.m,dcons.m文件,在命令窗口输入相应的命令,结果如下图所示:
调用相应的MATLAB的优化函数来解决,看所得结果是否相同。
编写MATLAB目标函数文件:
functionc,ce=ccon(x)ce=;c=-x
(1)2+6*x
(1)-4*x
(2)+11;x
(1)*x
(2)-3*x
(2)-exp(x
(1)-3)+1;x(
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 序列 二次 规划