二次约束二次规划问题的二元均值松弛定界算法

2021-02-05田福平高岳林

合肥工业大学学报（自然科学版） 2021年1期

田福平, 高岳林,2, 孙滢

(1.北方民族大学数学与信息科学学院,宁夏银川 750021; 2.合肥工业大学数学学院,安徽合肥 230601; 3.宁夏科学计算与智能信息处理协同创新中心,宁夏银川 750021)

本文主要考虑的二次约束二次规划((quadratically constrained quadratic programming,QQP)问题形式如下:

其中:Qs=(qsij)n×n为实对称矩阵;cs∈Rn;ds∈R,s=0,1,2,…,N;i,j=1,2,…,n;A∈Rm×n;b∈R。QQP问题的可行域记为F,假设F是非空有界闭集。

一般二次约束二次规划问题是非线性规划(nonlinear programming,NLP)问题,它在很多领域都有广泛的应用,如无线通信、网络、雷达、信号处理等[1-2]。当QQP问题中的Q0,Q1,Q2,…,Qs均为半正定矩阵时,可以利用文献[3]中的二阶锥规划方法对其进行求解;另一方面,一些特殊的非凸QQP问题可以等价地转化为凸问题[4-6]。但是,一般二次约束二次规划问题通常是NP难问题,因此很难求到此类问题的全局最优解。

学者们提出了各种各样的用于求解非凸二次规划问题的分支定界算法[7-8]。同时,也有一些可以用来求解非凸QQP问题全局最优解的求解器,如SCIP、GloMIQO、Ipopt、ANTIGONE、BARON、COUENNE等。这些算法和求解器已经在某些类型的问题上有了十分良好的性能,如BARON、SCIP、COUENNE求解器对于求解具有稀疏参数的线性约束二次规划问题是非常有效的;文献[9]提出的算法可以有效地求解带有盒约束的二次规划问题;文献[7]提出的算法专门用于求解带有可分离约束的混合整数二次规划问题。但是对于一般的QQP问题,目前尚没有较好的求解方法。在利用分支定界算法求解极小化问题时最关键的是定下界,能否成功构造出最优值的下界序列对算法设计成功与否是至关重要的,定下界过程要求把下界问题转化为一系列尽可能简单并易于求解的松弛子问题。文献[10]将松弛方法分为线性松弛、凸包络松弛、拉格朗日对偶松弛、二阶锥规划松弛和半定规划松弛。在这些不同类型的松弛中,目前求解全局最优解时最常用的方法是线性松弛。

本文通过引入辅助乘积变量将原QQP问题等价地转化为NLP问题,基于二元均值不等式将等价非线性规划问题转化为松弛线性规划(relaxation linear pngramming,RLP)问题提出了一种新的剖分策略,给出了超矩形缩减的方法及其算法的描述,并证明了该算法的收敛性,通过数值实验验证了本文提出的算法是可行有效的。

1 基于辅助变量的等价非线性规划

本文通过求解以下2n个线性规划得到x的上下界:

(1)

(2)

由(1)式得最优值lj,j=1,2,…,n;由(2)式得最优值uj,j=1,2,…,n,l=(l1,l2,…,ln),u=(u1,u2,…,un),令超矩形H={x|l≤x≤u},显然超矩形H是包含可行域F的最小超矩形。

任取超矩形H的一个子超矩形Hk⊆H,原问题的子问题QQP(Hk)可以写成如下形式:

QQP(Hk):

通过引入一个新的变量Y=(yij)n×n=xxT=(xixj)n×n,可以将问题QQP(Hk)转化为以下等价问题:

其中,·表示2个矩阵的内积,即2个矩阵对应元素乘积的和。

因为Q0=(q0ij)n×n;Qs=(qsij)n×n,s=0,1,2,…,N,i、j=1,2,…,n；Y=(yij)n×n=xxT=(xixj)n×n,所以此问题NLP(Hk)可以展开成以下形式:

NLP(Hk):

引理1 问题QQP(Hk)与问题NLP(Hk)等价。

证明若(Y,x)是NLP(Hk)的可行解,则有：

从而x是QQP(Hk)的可行解,并且

f0(x)=xTQ0x+c0Tx=Q0·xxT+c0Tx=

Q0·Y+c0Tx=f0(Y,x)。

反之,如果x是QQP(Hk)的可行解,令Y=xxT,那么

Qs·Y+csTx=Qs·xxT+csTx=

xTQsx+csTx≤ds,s=1,2,…,N,

于是(Y,x)是YNLP(Hk)的可行解,并且

因此,问题QQP(Hk)与问题NLP(Hk)等价。

2 基于二元均值不等式的线性松弛

本文基于二元均值不等式ab≤(a2+b2)/2,结合线性函数来构造QQP(Hk)问题的松弛线性规划,该松弛线性规划为问题QQP(Hk)提供了一个最优值的下界,具体操作如下:

对于每一个i={1,2,…,n},计算

(3)

(4)

一方面,对于每一个i=1,2,…,n都有li≤xi≤ui,由xili≥0,xjlj≥0，以及二元均值不等式可知:

(5)

成立。整理不等式(5)有:

xixj≤(ljxi+lixj-lilj)+

[(xi-li)2+(xj-lj)2]/2≤

(6)

另一方面,设函数G(x)=x2,x∈R过点(l,l2)和(u,u2),得到关于x的直线函数:

U(x)=(l+u)(x-l)+l2=(l+u)x-lu。

当x=xi时,U(xi)=(li+ui)xiliui;当x=xj时,U(xj)=(lj+ui)xjljuj,并且G(xi)≤U(xi),G(xj)≤U(xj)。因此可得:

(7)

将(7)式带入(6)式，得

xixj≤[U(xi)+U(xj)+2ljxi+2lixj-

(2li+uj-lj)xj-liui-ljuj+(li-lj)2]/2。

因此,可以得到以下不等式:

yij=xixj≤[(2lj+ui-li)+

(2li+uj-lj)xj-liui-

ljuj+(li-lj)2]/2

(8)

基于二次均值不等式,结合函数的性质,可以得到松弛线性规划RLP(Hk)如下:

以下引理说明了问题QQP(Hk)、NLP(Hk)、RLP(Hk)的可行解以及最优解之间的关系。

引理2 若(Yk,xk)为问题RLP(Hk)的可行解,并且满足Yk=xk(xk)T,则xk为QQP(Hk)的可行解。

引理3 若(Yk,xk)为问题RLP(Hk)的最优解,并且满足Yk=xk(xk)T,则(Yk,xk)为问题NLP(Hk)的全局最优解,从而xk为QQP(Hk)的全局最优解。

证明设问题RLP(Hk)的最优值为α(Hk),问题QQP(Hk)的全局最优值为Z(Hk)。

因为(Yk,xk)为RLP(Hk)的最优解,所以α(Hk)≤Z(Hk),并且α(Hk)≤f0(Yk,xk)。因为yijk=xikxjk,所以由引理2知,xk∈F,从而Z(xk)≤f0(xk),并且f0(xk)=f0(Yk,xk),因此α(Hk)=Z(Hk)=f0(Hk),即xk为QQP(Hk)的全局最优解。

定义1 若(Yk,xk)为RLP(Hk)的解,其中xk∈F并且fs(xk)-fs(Yk,xk)≤ε,s=0,1,2,…,N,则xk称为问题QQP(Hk)的强ε-全局最优解。

定义2 若(Yk,xk)为RLP(Hk)的解,其中xk∈F,并且f0(xk)-f0(Yk,xk)≤ε,s=0,1,2,…,N的ε-全局最优解。

3 超矩形的剖分与缩减

3.1 新的分支方法

为了适当地提高算法的分支效率,本文采用了一种新的剖分方法,令Hk=[lk,uk]∈H表示当前所要剖分的超矩形,在给出具体的剖分方法之前,首先给出以下说明:对所有的变量xi,其中li≤xi≤ui,当xi∈Hk=[lk,uk]时,抛物线g(xi)=(xi)2与直线l(xi)=(lik+uik)xilikuik(i=1,2,…,n)围成的面积记为S(lik,uik),如图1所示。显然

图1 关于变量xi抛物线直线围成的面积

用xk表示问题QQP(Hk)的当前最优解,对应的松弛问题的最优解为(Yk,xk),显然,xk∈D∩Hk=[lk,uk],进行以下形式的剖分:

else找到第1个xidk∈argmaxω

end if;

通过x′与x″的连线或连线所在的平面将超矩形Hk剖分为Hk1=[lk1,uk1]和Hk2=[lk2,uk2]2个子超矩形,则这2个子超矩形分别为:

3.2 超矩形的缩减

为了提高算法的收敛速度,本文采用超矩形缩减技术[11],假设RLP中所有的线性不等式约束的展开形式为:

具体缩减方法描述如下:

令Ik∶={1,2,…,m};

fori=1,2,…,mdo begin

ifrLi>bithen 停止。问题RLP在Hk上没有可行解(Hk被删除);

else ifrLi

Ik∶=Ik-{i},QQP问题的第i个线性不等式约束被删除;

else

forj=1,2,…,ndo

ifaij>0 then

else ifaij<0 then

end if;

end do;

end if;

end do;

为了方便起见,本文把该方法产生的新的整超矩形仍然记为Hk,它是原来整超矩形的子集。

4 算法的描述和收敛性分析

首先,为了方便算法的叙述,假定算法迭代到第k步时,本文给出以下记号:F表示QQP问题的可行域;Q表示QQP问题的当前所有可行点的集合;Hk表示当前迭代步所要剖分的超矩形;Ω表示剪枝后剩下的所有超矩形的集合;U表示QQP问题当前全局最优值的上界;L表示QQP问题当前全局最优值的下界。

基于上述记号,下面给出求解QQP问题的分支定界算法具体描述:

(1) 初始化。通过 (3) 式、(4) 式构造关于x的初始超矩形H=H0=[l0,u0],在超矩形H0上求解RLP问题;如果其对应的最优解和最优值分别记为(Y0、x0)、L;那么L就是问题QQP当前全局最优值的一个下界,如果x0∈F,令Q=Q∪x0,初始上界为U=min{f(x):x∈Q},并找到QQP问题当前的一个最好的解,那么给定所需容忍度ε>0,Ω=H0;置k=1,转入步骤(2)。

(2) 终止准则。若U-L≤ε或Ω≠φ,则迭代终止,输出QQP问题当前全局最优解x*以及全局最优值f(x*);否则转步骤(3)。

(3) 选择规则。在Ω中选择L所对应的超矩形Hk,即L=L(Hk)，此时Ω=ΩHk。

Ω=Ω∪Hki,i={1,2}。

(5) 超矩形的缩减。运用3.2中超矩形的缩减技术,对剖分后的子超矩形进行缩减和删除,把缩减后的超矩形仍然记为Ω=Hki,i={1,2}。

(6) 剪枝规则。让Ω=Ω{H:L(H)≥U,H∈Ω}。

置k=k+1,转到步骤(2)。

为了证明收敛性,令

φ=[(2lj+ui-li)xi+(2li+uj-lj)xj-

liui-ljuj-(li-uj)2]/2。

定理1 令εi=ui-li,对每一个i∈{1,2,…,n}以及任意的x∈F,当εi→0时,都有|φ-yij|→0。

证明令不等式(8)的右侧[(2lj+ui-li)+(2li+uj-lj)xj-liui-ljuj+(li-lj)2]/2=φ,又因为(xili)(xjlj)≥0,所以有xixj≥ljxi+lixjlilj,因此

|φ-yij|≤|φ-ljxi-lixj+ujlj|=

liui-ljuj-(li-lj)2]-ljxi-lixj+

li)(xi-li)|+|(uj-lj)(xj-lj)|

又因为εi→0,所以|φ-yij|→0。

定理2 令εi=ui-li,对于每一个i∈{1,2,…,n}以及任意的x∈F,当εi→0时,都有maxS(li,ui)→0。

证明

maxS(li,ui)=

因为ui-li→0,所以maxS(li,ui)→0,从而对于每一个i都有S(li,ui)→0。

定理4 设QQP问题的可行域F是n维的,所要剖分的关于变量x的超矩形是n维的,且超矩形的剖分是耗尽的,若算法在有限步k终止,则xk必是问题的全局最优解,或者所产生的无穷可行点列{xk}的任一聚点都是QQP问题的全局最优解。

证明(1) 若算法是有限步终止且算法终止时迭代了k步,则根据终止准则可知Uk-Lk≤ε,即

f(xk)-Lk≤ε

(9)

假设此时的全局最优解为x*,可知：

Uk=f(xk)≥f(x*)≥Lk

(10)

因此,将 (9)式、(10) 式联立可得:

f(xk)+ε≥f(x*)+ε≥Lk+ε≥f(xk)

(11)

故(1)得证。

(2) 如果算法的迭代是无限的,QQP问题的可行解序列为{xk},对应的线性松弛问题的可行解序列为{(xk,Yk)}。

根据步骤(4)、步骤(5),随着x的超矩形的不断加(x*)细,则有:

Lk=f0(xk,Yk)≤f0(x*)≤f0(xk)=Uk,

k=1,2,…

(12)

因为下界序列{Lk=f0(xk,Yk)}是单调不减且有界的，上界序列{Uk=f0(xk)}是单调不增且有界,所以原问题的上下界都是收敛的。对(12)式两边取极限可得:

(13)

L=f(x*)=U

(14)

因此,序列{xk}的每一个聚点x*都是问题QQP的全局最优解。

5 数值实验

本文给出几个算例来证明本文算法的有效性。算法所有测试过程均用Matlab9.0.0.341360(R2015a),在Inter(R)Core (TM)i7-6700,CPU@3.40 GHz,16 GB内存,64位Windows7操作系统的计算机上运行，并且文中所有线性规划均用对偶单纯形法求解。

算例1[12-13]

算例2[14]

算例3[14]

算例4[12]

算例5[15-16]

算例6[15,17]

算例7[12,18]

算例8[16]

算例9[19]

算例10[20]

算例11

首先生成对称矩阵Q0、Qs,其中Q0和Qs(s=1,2,…,N)的每一个元素均在[-1,1]上随机生成,本文利用特征值分解可以得到Q0=PTD0P及Qs=PTDsP,从而得到对角矩阵D0和Ds,其中D0的前r个分量在[-10,0]上随机生成；其余对角分量在[0,10]上随机生成；Ds的对角分量在[1,100]上随机生成；c0为n维零向量；cs为在[-100,100]上随机生成的n维向量；ds为在[1,50]上随机生成的数。

算例1～算例10的计算结果见表1所列。表1中,x*为原问题的最优解;f(x*)为目标函数的最优值;I为迭代次数;ε为容忍度。

表1 算例1～算例10的计算结果

从表1可以看出,对于算例2、3，从迭代次数上看本文算法均优于其他算法,计算结果与其他算法相同;针对算例4、10，从计算结果可知本文算法分别要优于文献[12,20]算法,迭代次数大于文献[12,20]; 对于算例1,本文算法计算结果与其他算法相同,迭代次数大于其他算法;针对算例6、9,本文算法的迭代次数明显小于其他算法,计算结果与其他算法相当;最后针对算例5、7、8,从迭代次数和计算结果上看本文的算法要弱于其他算法。

以上对算例1～算例10 的计算结果进行分析，说明了本文算法是可行有效的。

为了进一步说明本文算法的有效性,针对算例11,进行了若干的随机实验,结果见表2所列,表2中，AI为平均迭代次数;At为迭代终止时CPU平均运行时间;m为二次约束的个数;n为目标函数的维数;r为目标函数矩阵中负特征值的个数。

表2中的容忍度均为5×10-4。通过表2的计算结果可以看出,当n和r固定时,随着二次约束m个数的增多,平均迭代次数AI和迭代时间At都在增加;若n和m固定,当r