应用数学MATHEMATICA APPLICATA 2021,34(1):240-252一种新的二次约束二次规划问题的分支定界算法黄小利1,2,高岳林1,2,谢金宵1,2,谷剑峰1,2(1.北方民族大学数学与信息科学学院,宁夏银川750021;2.宁夏科学计算与智能信息处理协同创新中心,宁夏银川750021)摘要：本文为了获得二次约束二次规划(QCQP)问题的全局最优解,提出一种新的参数化线性松弛分支定界算法.该算法利用参数化线性松弛技术,得到(QCQP)的全局最小值的下界,并利用区域缩减技术以最大限度地删除不可行区域,加快该算法的收敛速度.数值实验表明,本文提出的算法是有效并且可行的.关键词：二次约束二次规划;全局优化;分支定界;参数化线性松弛;区域缩减中图分类号：O221.2AMS(2000)主题分类：90C20;90C30;90C57文献标识码：A 文章编号：1001-9847(2021)01-0240-131.引言本文主要考虑以下形式的二次约束二次规划问题:(QCQP):min f 0(x )=n ∑j =1n ∑q =1a 0jq x j x q +n ∑j =1a 0j x j +ε0,s.t.f i (x )=n ∑j =1n ∑q =1a ijq x j x q +n ∑j =1a i j x j +εi ≤0,i =1,...,M x ∈X 0={x ∈R n |Ax ≤b,l 0≤x ≤u 0},这里f i (x ),i =0,1,...,M 是非凸二次函数.(a i jq )n ×n 是n ×n 阶实对称矩阵,a ij ,εi ∈R ,i =0,1,...,M,j =1,...,n,q =1,...,n,A ∈R m ×n ,b ∈R m ,l 0=(l 01,...,l 0n )T ,u 0=(u 01,...,u 0n )T ,X 0为非空有界闭集,记H=[l 0,u 0],H 0是超矩形.非凸二次约束二次规划(QCQP)问题在实际生活中应用范围广泛[1−4],并且通常是全局优化问题,但在解决过程中往往存在许多局部最优解而不是全局最优解,且(QCQP)问题是一个NP-hard 问题[5],这就使得在理论和计算方面存在巨大的挑战.因此,寻找一种有效的算法全局求解(QCQP)问题是十分有必要的.从(QCQP)问题的研究历程来看,已经有许多算法都适用于求解全局(QCQP)问题,现如今分支定界算法已成为求解该问题最常用的工具之一.根据分支定界算法的框架,基于对(QCQP)问题的松弛,在分支定界树的每个节点求解松弛子问题的下界,得到一个高质量的下界对原问题起着至关重要的作用.目前的(QCQP)松弛方法建立在多种松弛技术上,包括∗收稿日期：2020-07-08基金项目：国家自然科学基金项目(11161001;61561001);宁夏高等教育一流学科建设基金(NXYLXK 2017B09);北方民族大学研究生创新项目(YCX20103);北方民族大学重大专项(ZDZX201901)作者简介：黄小利,女,汉族,陕西人,研究方向:最优化理论与方法.通讯作者：高岳林.第1期黄小利等：一种新的二次约束二次规划问题的分支定界算法241基于凹凸包络的线性规划松弛[6−7],拉格朗日对偶松弛[6,8],基于二阶锥规划(SOCP)的松弛[9],基于半定规划(SDP)的松弛[10].在求解盒约束非凸二次规划时,文[11]在有限分支中使用半定规划松弛技术.在求解带线性约束的非凸二次规划问题(LCQP)时,文[12]基于一阶KKT条件的有限分支和多面体半定松弛求解(LCQP).在求解带球和线性约束的非凸二次规划问题时,文[13]利用了半定规划松弛逼近非凸二次规划问题.在求解带一个非凸二次约束的二次规划问题时,文[14]将问题转化为一个无非凸约束的参数二次规划问题,通过迭代求解并在每次迭代中更新参数.在求解带多个二次约束的非凸二次规划问题时,文[15]给出了非凸二次约束二次规划(QCQP)问题的新的凸松弛方法,将每个非凸约束分解,并松弛为两个二阶锥约束,将这些非凸约束与线性约束的乘积线性化,以构造一些新的有效约束;文[16]基于特征值分解的分支定界算法来求解(QCQP)问题,在特征值分解中,非凸分量由负特征值和相应的特征向量表示,提出了一种基于半定松弛的分支定界算法来求解(QCQP);文[17]用一种线性化方法将初始非凸规划问题化为一系列线性松弛规划问题;文[8]将(QCQP)问题转化为拉格朗日对偶优化问题,并在更新拉格朗日乘子的同时连续求解子问题;文[18-20]提出了一种参数线性化方法来生成(QCQP)的参数线性松弛规划,通过对可行域和目标函数的线性松弛,将初始非凸问题化为一系列线性规划问题;文[21]利用重构线性化技术在分支树内通过连续线性化来估计所有二次项;文[22]基于二次约束线性规划的半定松弛,引入了一种特殊的惩罚方法,得到了所谓的条件拟凸松弛,证明了条件拟凸松弛比标准半定松弛能提供更紧的界.非凸(QCQP)问题的困难源于二次项的非凸分量,松弛过程中只需对目标和约束函数中的二次项的非凸分量部分进行操作.本文给出了求解非凸二次约束二次规划(QCQP)问题的全局优化算法.该算法结合分支定界操作和区域缩减技术,通过对非凸二次项进行新的参数化线性松弛,使用线性松弛技术将非凸(QCQP)问题转化为线性松弛规划问题,提出的线性化技术在不增加新的变量和约束的情况下,将线性松弛规划问题嵌入到分支定界操作中.为了加快算法的收敛速度,在分支和定界过程中采用了区域缩减技术,为当前不存在全局最优解的区域提供了理论可能性,可以显著减少松弛间隙.最后,数值实验结果表明,该方法能较好地解决所有存在于全局最优解中的(QCQP)问题.本文的组成结构如下:第一部分提出了二次约束二次规划问题,并且概述了在不同的约束条件下求解二次规划问题的相关算法.第二部分提出了新的参数化线性松弛技术从而得到含参数的线性松弛规划问题.第三部分为了加快算法收敛速度从而提出区域缩减技术.第四部分将第二部分得到的参数化线性松弛规划,以及第三部分的区域缩减技术嵌套在分支定界框架内,得到参数化线性松弛分支定界算法以求解非凸二次约束二次规划(QCQP).第五部分通过实例证明算法的有效性和可行性.第六部分得出相关结论.2.新的参数化线性松弛技术为了求解原问题(QCQP),本文采用一种新的参数化线性松弛技术得到原问题和其子问题的下界函数.令H={x=(x1,x2,...,x n)T∈R n,l≤x≤u}⊆H0,其中l=(l1,...,l n)T,u= (u1,...,u n)T.与文[18]类似,我们定义ρ=(ρjq)n×n是一个对称参数矩阵且ρjq∈{0,1}.对于任何的x∈H,且j∈{1,2,...,n},q∈{1,2,...,n},j=q,ρjq∈{0,1}给出以下定义:x j(ρjq)=l j+ρjq(u j−l j);x j(1−ρjq)=l j+(1−ρjq)(u j−l j);x q(ρjq)=l q+ρjq(u q−l q);x q(1−ρjq)=l q+(1−ρjq)(u q−l q).对于非凸二次函数的非线性项,接下来我们提出新的线性松弛技术,定义如下:ϕ(x j,x q)=x j x q;ϕ(x j,x q)=14{[x j(ρjq)+3x q(ρjq)+x q(1−ρjq)−x j(1−ρjq)]x j242应用数学2021+[3x j(ρjq)+x q(ρjq)−x q(1−ρjq)+x j(1−ρjq)]x q−[x j(ρjq)+x q(ρjq)]2+[x j(ρjq)−x q(1−ρjq)][x j(1−ρjq)−x q(ρjq)]};ϕ(x j,x q)=14{[3x q(1−ρjq)−x j(ρjq)+x q(ρjq)+x j(1−ρjq)]x j+[3x j(ρjq)−x q(1−ρjq)+x q(ρjq)+x j(1−ρjq)]x q−[x j(ρjq) +x q(ρjq)][x j(1−ρjq)+x q(1−ρjq)]+[x j(ρjq)−x q(1−ρjq]2};ϑ(x j,x q)=ϕ(x j,x q)−ϕ(x j,x q), ϑ(x j,x q)=ϕ(x j,x q)−ϕ(x j,x q).对于任何的x∈H,且j∈{1,2,...,n},q∈{1,2,...,n},j=q,ρjj∈{0,1},ϕ(x j)=x2j,ϕ(x j)=2x j(ρjj)x j−x2j(ρjj);ϕ(x j)=(x j(ρjj)+x j(1−ρjj))x j−x j(ρjj)x j(1−ρjj);ϑ(x j)=ϕ(x j)−ϕ(x j); ϑ(x j)=ϕ(x j)−ϕ(x j);显然,x j(0)=l j,x j(1)=u j,x q(0)=l q,x q(1)=u q.定理1对于任意的x∈H⊆H0,有(a)对任意x j∈[l j,u j],j∈{1,2,...,n},有ϕ(x j)≤ϕ(x j)≤ϕ(x j);(b)对任意j∈{1,2,...,n},q∈{1,2,...,n},j=q,有ϕ(x j,x q)≤ϕ(x j,x q)≤ϕ(x j,x q);(c)当∥l−u∥→0时,得到ϑ(x j), ϑ(x j),ϑ(x j,x q), ϑ(x j,x q)→0.证(a)在每一个区间[l j,u j],j∈{1,2,...,n}对单变量函数ϕ(x j)=x2j进行线性缩放可得ϕ(x j)的上下界逼近函数,进而ϕ(x j)−ϕ(x j)=x2j−(2x j(ρjj)x j−x2j(ρjj))=(x j−x j(ρjj))2≥0,ϕ(x j)−ϕ(x j)=[x j(ρjj)+x j(1−ρjj)]x j−x j(ρjj)x j(1−ρjj)−x2j≥0.故ϕ(x j)≤ϕ(x j)≤ϕ(x j);(b)把(x j+x q)和(x j−x q)分别看成一个单变量,则(x j+x q)2和(x j−x q)2分别是在区间[l j+l q,u j+u q]和[l j−u q,u j−l q]上关于(x j+x q)以及(x j−x q)的凸函数,则由(a)(x j+x q)2≥2[x j(ρjq)+x q(ρjq)](x j+x q)−[x j(ρjq)+x q(ρjq)]2;(2.1)(x j+x q)2≤[x j(ρjq)+x q(ρjq)+x j(1−ρjq)+x q(1−ρjq)](x j+x q)−[x j(ρjq)+x q(ρjq)][x j(1−ρjq)+x q(1−ρjq)];(2.2) (x j−x q)2≥2[x j(ρjq)−x q(1−ρjq)](x j−x q)−[x j(ρjq)−x q(1−ρjq)]2;(2.3)(x j−x q)2≤[x j(ρjq)−x q(1−ρjq)+x j(1−ρjq)−x q(ρjq)](x j−x q)−[x j(ρjq)−x q(1−ρjq)][x j(1−ρjq)−x q(ρjq)].(2.4)结合式(2.1)-(2.4),ϕ(x j,x q)=14[(x j+x q)2−(x j−x q)2]≥14{2[x j(ρjq)+x q(ρjq)](x j+x q)−[x j(ρjq)+x q(ρjq)]2−[x j(ρjq)−x q(1−ρjq) +x j(1−ρjq)−x q(ρjq)](x j−x q)+[x j(ρjq)−x q(1−ρjq)][x j(1−ρjq)−x q(ρjq)]}=14{[x j(ρjq)+3x q(ρjq)+x q(1−ρjq)−x j(1−ρjq)]x j+[3x j(ρjq)+x q(ρjq)−x q(1−ρjq)+x j(1−ρjq)]x q−[x j(ρjq)+x q(ρjq)]2+[x j(ρjq)−x q(1−ρjq)][x j(1−ρjq)−x q(ρjq)]}第1期黄小利等：一种新的二次约束二次规划问题的分支定界算法243 =ϕ(x j,x q);ϕ(x j,x q)=14[(x j+x q)2−(x j−x q)2]≤14{[x j(ρjq)+x q(ρjq)+x j(1−ρjq)+x q(1−ρjq)](x j+x q)−[x j(ρjq)+x q(ρjq)][x j(1−ρjq)+x q(1−ρjq)]−2[x j(ρjq)−x q(1−ρjq)](x j−x q)+[x j(ρjq)−x q(1−ρjq)]2}=14{[3x q(1−ρjq)−x j(ρjq)+x q(ρjq)+x j(1−ρjq)]x j+[3x j(ρjq)−x q(1−ρjq)+x q(ρjq)+x j(1−ρjq)]x q−[x j(ρjq)+x q(ρjq)][x j(1−ρjq)+x q(1−ρjq)]+[x j(ρjq)−x q(1−ρjq]2}=ϕ(x j,x q).因此ϕ(x j,x q)≤ϕ(x j,x q)≤ϕ(x j,x q).(c)由于ϑ(x j)=ϕ(x j)−ϕ(x j)=[x j−x j(ρjj)]2是在区间[l j,u j]上的凸函数,所以ϑ(x j)在l j或者u j可以取得最大值.即max x j∈[l j,u j]ϑ(x j)=0;maxx j∈[l j,u j]ϑ(x j)=(l j−u j)2.(2.5)由 ϑ(x j)=ϕ(x j)−ϕ(x j)=[x j(ρjj)+x j(1−ρjj)]x j−x j(ρjj)x j(1−ρjj)−x2j可得, ϑ(x j)在l i+u j2取得最大值,即max x j∈[l j,u j] ϑ(xj)=(l j−u j)24.(2.6)由(2.5)和(2.6)式可知,当|l j−u j|→0时,maxx j∈[l j,u j]ϑ(x j)→0;(2.7)max x j∈[l j,u j] ϑ(xj)→0.(2.8)令ϑ(x j+x q)=ϕ(x j+x q)−ϕ(x j+x q), ϑ(x j+x q)=ϕ(x j+x q)−ϕ(x j+x q),ϑ(x j−x q)=ϕ(x j−x q)−ϕ(x j−x q), ϑ(x j−x q)=ϕ(x j−x q)−ϕ(x j−x q),使用同样的方法,当∥l−u∥→0时,max(x j+x q)∈[l j+l q,u j+u q]ϑ(x j+x q)→0;(2.9)max(x j+x q)∈[l j+l q,u j+u q] ϑ(xj+x q)→0;(2.10)max(x j−x q)∈[l j−u q,u j−l q]ϑ(x j−x q)→0;(2.11)max(x j−x q)∈[l j−u q,u j−l q] ϑ(xj−x q)→0.(2.12)因为ϑ(x j,x q)=ϕ(x j,x q)−ϕ(x j,x q)=14[(x j+x q)2−(x j−x q)2]−14{2[x j(ρjq)+x q(ρjq)](x j+x q)−[x j(ρjq)+x q(ρjq)]2−[x j(ρjq)−x q(1−ρjq)+x j(1−ρjq)−x q(ρjq)](x j−x q)−[x j(ρjq)−x q(1−ρjq)][x j(1−ρjq)−x q(ρjq)]}=14{(x j+x q)2−2[x j(ρjq)+x q(ρjq)](x j+x q)+[x j(ρjq)+x q(ρjq)]2}244应用数学2021+14{[x j (ρjq )−x q (1−ρjq )+x j (1−ρjq )−x q (ρjq )](x j −x q )−[x j (ρjq )−x q (1−ρjq )][x j (1−ρjq )−x q (ρjq )]−(x j +x q )2}≤14max (x j +x q )∈[l j +l q ,u j +u q ]ϑ(x j +x q )+14max (x j −x q )∈[l j −u q ,u j −l q ]ϑ(x j −x q ).根据(2.9)及(2.12)式得,当∥l −u ∥→0时,有ϑ(x j ,x q )→0.和以上证明方法类似,当∥l −u ∥→0时,我们可以得到 ϑ(x j ,x q )→0.证毕.对于任意x j ∈[l j ,u j ],j ∈{1,...,n },我们得到f l i (x ),i =0,1,...,M.即f i (x )=n ∑j =1n ∑q =1,j =qa i jq x j x q +n ∑j =1a i jj x 2j +n ∑j =1a i j x j +εi=n ∑j =1,a i jq >0n ∑q =1,j =q a i jq x j x q+n ∑j =1,a i jq <0n ∑q =1,j =q a i jq x j x q+n ∑j =1,a i jj >0a i jj x 2j +n ∑j =1,a i jj <0a i jj x 2j +n ∑j =1a i j x j +εi≥n ∑j =1,a i jq >0n ∑q =1,j =q a i jq ϕ(x j ,x q )+n ∑j =1,a i jq <0n ∑q =1,j =q a i jq ϕ(x j ,x q )n ∑j =1,a i jj >0a i jj ϕ(x j )+n ∑j =1,a i jj <0a i jj ϕ(x j )+n ∑j =1a i j x j +εi=f l i (x ).(2.13)通过以上讨论,我们得到原问题(QCQP)在子矩形H 上的线性松弛规划问题如下:(LRP):min f l 0(x ),s.t.f l i (x )≤0,i =1,...,M,x ∈X ={x ∈R n |Ax ≤b,l ≤x ≤u }.定理2对于任意的x ∈H =[l,u ]⊂H 0,有(a)f i (x )≥f l i (x ),i =0,1,...,M ;(b)当∥l −u ∥→0时,f i (x )−f l i (x )→0,i =0,1,...,M.证(a)由(2.13)式我们很容易得到f i (x )≥f l i (x ),i =0,1,...,M ;(b)由f l i (x )的定义和结论(a),对每一个i =0,1,...,M,我们有f i (x )−f li (x )=n ∑j =1n ∑q =1,j =qa i jq x j x q +n ∑j =1a i jj x 2j +n ∑j =1a i j x j +εi−[n ∑j =1a i j x j+n ∑j =1,a i jj >0a i jj ϕ(x j )+n ∑j =1,a i jj <0a i jj ϕ(x j )+n ∑j =1,a i jq >0n ∑q =1,j =q a i jq ϕ(x j ,x q )+n ∑j =1,a i jq <0n ∑q =1,j =q a i jq ϕ(x j ,x q )+εi ]=n ∑j =1,a i jj >0a i jj (x 2j −ϕ(x j ))−n ∑j =1,a i jj <0a i jj (ϕ(x j )−x 2j )第1期黄小利等：一种新的二次约束二次规划问题的分支定界算法245+n∑j=1,a ijq>0n∑q=1,j=qa ijq(x j x q−ϕ(x j,x q))−n∑j=1,a ijq<0n∑q=1,j=qa ijq(ϕ(x j,x q)−x j x q).由定理1的(c)知,当∥l−u∥→0时,ϑ(x j), ϑ(x j),ϑ(x j,x q), ϑ(x j,x q)→0.因此对于每一个i=0,1,...,M,当∥l−u∥→0时,有f i(x)−f li(x)→0,i=0,1,...,M.证明完毕.由定理2可以知道,对每一个子矩形H来讲,线性松弛规划问题(LRP)的目标函数和约束条件函数是小于等于原问题(QCQP)的目标函数和约束条件函数的,也就是说定理2保证了问题(LRP)的最优值为原问题(QCQP)能够提供可靠的下界,当∥l−u∥→0时,线性松弛规划问题(LRP)限逼近原问题(QCQP),这说明了算法是全局收敛的.3.缩减技术为了加快算法的收敛速度和计算速度,在这一部分,基于文[17-18],提出适用于本算法的区域缩减技术.当算法迭代到第k步时我们需要判断在矩形区域H上原问题(QCQP)是否存在全局最优解.此处缩减技术的前提是未删掉在矩形区域H上原问题(QCQP)的全局最优解,进而对区域进行缩减,具体见定理3和定理4.不失一般性,假设线性松弛规划问题LRP中的每一个线性函数φi(x),i=0,1,...,M,M+ 1,...,M+m,可以表示为φi(x)=n∑j=1αij x j+ξi,i=0,1,...,M,M+1,...,M+m.设UB k是算法迭代到第k步原问题(QCQP)当前已知的最好上界,并令µk s =UB k−(n∑j=1,j=smin{αk0jx j+ξk}),s=1,...,n,δk is =b ki−(n∑j=1,j=smin{αkijl kj,αkiju kj}+ξk i),i=1,...,M,M+1,...,M+m,s=1,...,n.定理3对于任何的子矩形H k⊆H⊆H0且H k=[l k,u k].如果存在s∈{1,2,...,n}使得µs<min{αk0s l ks,αk0su ks},则在子矩形H k上原问题(QCQP)没有全局最优解;否则,如果存在s∈{1,2,...,n}使得当αk0s >0且µks<αk0su ks,则原问题(QCQP)在子矩形H k1=(H k1j)n×1上没有全局最优解;当αk0s <0且µks<αk0sl ks,则原问题(QCQP)在子矩形H k2=(H k2j)n×1上没有全局最优解.其中,H k1j =H kj,j=s,j=1,...,n,(µksαk0s,u ks]∩H ks,j=s;H k2j=H kj,j=s,j=1,...,n,[l ks,µksαk0s)∩H ks,j=s.证如果存在s∈{1,2,...,n}使得µks <min{αk0sl ks,αk0su ks},则UB k−n∑j=1,j=nmin{αk0jx j+ξk}<min{αk0s l k s,αk0s u k s},即UB k<n∑j=1,j=nmin{αk0jx j+ξk}+min{αk0s l k s,αk0s u k s}≤n∑j=1αk0jx j+ξk=φk(x).246应用数学2021因此,原问题(QCQP)在H k上没有全局最优解.接下来记x的第s个分量为x s,因为x s∈(µk sαk0s ,u ks],所以有µk sαk0s <x ks≤u k s.当αk0s>0时,µk s<αk0s x k s,对于一切的x∈H k1,由µk s的定义和上面的不等式,有UB k−n∑j=1,j=smin{αk0jx j+ξk}<αk0s x k s,也就是UB k<n∑j=1,j=smin{αk0jx j+ξk}+αk0s x k s≤n∑j=1αk0jx j+ξk=φk(x).因此,对于任意x∈H k1,线性松弛函数φ0(x)要大于原问题(QCQP)的上界UB k,所以在H k1上不存在原问题(QCQP)的全局最优解.类似地,对于所有的x∈H k2,存在s∈{1,2,...,n}使得αk0s <0且µks<αk0sl ks时,在H k2这个子矩形上也不存在原问题(QCQP)的全局最优解.证毕.定理4对于任何的子矩形H k⊆H0且H k=[l k,u k].如果δk is<min{αk is l k s,αk is u k s},i= 1,...,M,M+1,...,M+m,则原问题在H k上没有全局最优解.否则,若存在s∈{1,2,...,n}使得αkis >0且δkis<αkisu ks,i=1,...,M,M+1,...,M+m,则原问题在子矩形H k3上没有全局最优解;若存在s∈{1,2,...,n}使得αkis <0且δkis<αkisl ks,则原问题在子矩形H k4上没有全局最优解.其中,H k3j =H kjj=s,j=1,...,n,(δkisαkis,u ks]∩H ksj=s;H k4j=H kjj=s,j=1,...,n,[l ks,δkisαkis]∩H ksj=s;证使用与证明定理3类似的方法,我们可以得到此定理成立.由定理3和定理4,采用缩减技术删掉不含有原问题(QCQP)全局最优解的区域.令H= (H j)n×1是H0任一子矩形,其中H j=[l j,u j].得到缩减规则如下:1)如果µs<min{αk0s l ks,αk0su ks},则H=∅;否则,如果αk0s>0且µk s<αk0s u k s,则令u k j=min{µk sαk0s ,u kj};H=(H j)n×1,H j=[l j,u j](j=1,...,n);如果α0s<0且µk s<α0s l k s,则令l k j=max{l kj ,δk isαkis},H=(H j)n×1,H j=[l j,u j](j=1,...,n).2)对每一个i=1,...,M,M+1,...,M+m,如果δkis <min{αkisl ks,αkisu ks},则H=∅;否则,如果αkis >0且δkis<αkisu ks则令u kj=min{δk isαkis,u kj},H=(H j)n×1,H j=[l j,u j](j=1,...,n);如果αkis <0且δkis<αkisl ks,则令l kj=max{l kj,δk isαkis},H=(H j)n×1,H j=[l j,u j](j=1,...,n).4.分支定界算法及其收敛性分析为了得到原问题(QCQP)的全局最优解,提出了参数化线性松弛分支定界算法.在矩形H⊆H0被剖分后的子集上求解对应的线性松弛规划问题.分支过程中,需要依据一定的剖分规则将初始超矩形H0分成两个新的子矩形.本文所采用的是标准矩形二分方法.考虑任一通过H={x∈R n|l j≤x j≤u j,j=1,...,n}⊆H0所确定的子节点问题.分支规则如下所示:1)令p=arg max{u j−l j|j=1,...,n};2)令ωp=l j+u j2;3)设新的子矩形区域H1={x∈R n|l j≤x j≤u j,j=q,l q≤x q≤ωp};H2={x∈R n|l j≤x j≤u j,j=q,ωp≤x q≤u q}.第1期黄小利等：一种新的二次约束二次规划问题的分支定界算法247通过以上分支规则,将矩形区域H剖分成了两个子矩形H1和H2.为了方便起见,我们把缩减后产生新的矩形区域仍然记为H k,它是初始超矩形区域的子集.记LB k是线性松弛规划问题(LRP)在子矩形H k上的最优值,x k=x(H k)是相对应的最优解.结合前面所提到的线性松弛规划问题,区域缩减技术以及分支规则,为求解问题(QCQP)提出分支定界算法,步骤如下:步0(初始化)设置误差精度ϵ>0,初始迭代次数k=0,Q表示剪支后所有矩形的集合,初始的可行集为W=∅,初始上界UB0=+∞.在初始矩形H0上求解线性规划问题LRP(H0),如果不可行,则原问题(QCQP)无解;否则将求得的最优值和最优解分别记为LB0和x0.如果x0对(QCQP)可行,置可行集为W= W∪{x0},则初始上界UB0=f0(x0),Q=Q∪{H0},置k=1,转步1;步1(终止准则)如果UB k−LB k<ϵ,则迭代终止,x k是问题(QCQP)的ϵ-全局最优解.否则,转步2;步2(矩形区域的分支)在Q中选择出当前最小下界所对应的矩形H k,通过本文所提出的分支规则将矩形区域H k剖分为两个新的子矩形H k1和H k2,用H k表示新剖分的矩形集合;步3(矩形区域的缩减)对每一个子矩形H kt∈H k,t=0,1,2,利用本文所给出的区域缩减技术分别对子矩形进行缩减;步4(剪支规则)令Q=Q\{H k:LB k−UB k>ε,H k∈Q},转步5;步5(更新上界)如果W=∅,那么UB k保持不变;如果W=∅,则更新上界,UB k= min{f0(x):x∈W},当前最好的可行点为x∗:=arg min x∈W f0(x);步6(更新下界)如果Q=∅,那么LB k保持不变;如果Q=∅,那么LB k=min{LRP: H k∈Q},置k=k+1,返回步1继续执行.为了证明本文的分支定界算法是全局收敛的,故给出下面的定理5,以保证当∥l−u∥→0时线性松弛规划问题(LRP)地逼近于原问题(QCQP).定理5假设(QCQP)问题的可行域是非空的,且矩形的剖分是耗尽的,则以下两个结论成立:(a)如果算法迭代在有限步终止,则(QCQP)问题的全局最优解将在算法终止时得到;(b)如果算法迭代在有限步不能终止,并且产生无穷多个可行点的序列{x k},则{x k}的任何一个聚点x∗都是问题(QCQP)的全局最优解.证(a)如果算法迭代在有限步终止,那么我们假设算法是在第k次迭代终止,k≥0.根据算法的步1得到UB k−LB k≤ϵ,表示原问题(QCQP)存在可行解并设为x∗,使得UB k=f0(x∗),同时有f0(x∗)−LB k≤ϵ.(4.1)假设v是原问题(QCQP)的全局最优值,由下界的计算过程表明LB k≤v.(4.2)因为x∗是原问题(QCQP)的可行解,有f0(x∗)≥v.(4.3)联立式(4.1)-(4.3)得v≤f0(x∗)≤LB k+ϵ≤v+ϵ.(4.4)所以,x∗是问题(QCQP)的ϵ-全局最优解.(b)如果算法的迭代是无限的,则通过求解一系列问题(LRP)可产生原问题的可行解序列{x k}∈H0.根据分支定界算法的结构,随着可行域的不断剖分加细我们有(x k)≤f0(x∗)≤f(x k)=UB k,k=1,2,···.(4.5) LB k=f l248应用数学2021因为下界序列{LB k=f l(x k)}是单调不减且有上界,且上界序列{UB k=f0(x k)}是单调不增且有下界,故原问题的上下界都是收敛的.对式(4.5)两边取极限可得lim k→∞f l(x k)≤f0(x∗)≤limk→∞f0(x k).(4.6)于是,令LB=limk→∞f l(x k),UB=limk→∞f0(x k),那么式(4.6)将转化为LB≤f0(x∗)≤UB.(4.7)不失一般性,假设超矩形序列{H k=[l k,u k]}满足x k∈H k以及H k+1⊆H k.在算法进行迭代过程中,最终矩形的剖分是耗尽的,即∞∩k=1H k={x k},根据函数f0(x)的连续性,我们有LB≤f0(x∗)≤limk→∞f0(x k)=UB.(4.8)所以,序列{x k}的每一个聚点x∗都是问题(QCQP)的全局最优解.5.数值实验给出以下几个算例证明本文算法的有效性.本文算法中所有的线性规划问题均选用对偶单纯形法求解,算法所有测试过程均用MATLAB9.2.0.538062(R2017a)在Inter(R)Core (TM)i5-8250U,***********,8GB内存,64位Windows10操作系统的计算机上运行.算例1[17−18]:min−x21+x1x2+x22+x1−2x2 s.t.x1+x2≤6,−2x21+x22+2x1+x2≤−4,1≤x1,x2≤6.算例2[17−18]:min x21+x22s.t.0.3x1x2≥1,2≤x1≤5,1≤x2≤3.算例3[17−18]:{min6x21+4x22+5x1x2,s.t.−6x1x2≤−48,0≤x1,x2≤10算例4[17]:min x1x2−2x1+x2+1s.t.8x22−6x1−16x2≤−11,−x22+3x1+2x2≤7,1≤x1≤2.5,1≤x2≤2.225.算例5[17]:min x1s.t.14x1+12x2−116x21−116x22≤1, 114x21+114x22−37x1−37x2≤−1, 1≤x1≤5.5,1≤x2≤5.5第1期黄小利等：一种新的二次约束二次规划问题的分支定界算法249算例6[17]:min−4x2+(x1−1)2+x22−10x23s.t.x21+x22+x23≤2,(x1−2)2+x22+x23≤2,2−√2≤x1≤√2,0≤x2,x3≤√2.算例7[21]:min f0(x)=12⟨x,Q0x⟩+⟨d0x⟩s.t.f i(x)=12⟨x,Q i x⟩+⟨d i x⟩≤b i,i=1,2,...,m, 0≤x j≤10,j=1,2,...,n,其中Q0∈R n×n的每个元素在区间[0,1]随机产生,Q i∈R n×n(i=1,2,...,m)的每个元素在区间[-1,0]随机产生;d0∈R n的每个元素在区间[0,1]随机产生;d i∈R n(i=1,2,...,m)的每个元素在区间[-1,0]随机产生;每一个b i∈R(i=1,2,...,m)在区间[-300,-90]产生.算例8[17−18]:min−n∑i=1x2is.t.n∑i=1x i≤j,j∈{1,2,...,n}, x i≥0,i∈{1,2,...,n}.表5.1算例1-6利用本文算法产生的数值结果E Refs Matrixρx∗f0(x∗)Iters Time(s) 1本文(00;00)(5.0,1.0)-16.0130.361252 1[17]-(5.0,1.0)-16.0130.37129 1[18](00;00)(5.0,1.0)-16.0140.482本文(01;10)(2.00000,1.66667) 6.77778250.73218 2[17]-(2.000534,1.666430) 6.779125285660.59 2[18](1,0;0,1)(2.000000,1.6666668) 6.77777838013.26 3本文(01;10)(2.555705,3.130251)118.3836180249.83 3[17]-(2.560885,3.123981)118.386615823321.796 3[18](00;00)(2.555765,3.1301777)118.38367173334.82 4本文(10;01)(2.0,1.0)-1.010.072806 4[18](01;10)(2.0,1.0)-1.02 1.5267 5本文(01;10)(1.177125;2.177124) 1.177125138 3.895[17]-(1.177391,2.176737) 1.17739191 2.845[18](01;10)(1.177125,2.177124) 1.177125131 3.766本文(110;111;011)(1.0,0.315864,0.924541)-10.7114433048.436 6[17]-(1.0,0.181981,0.983302)-11.363636153043.4782 6[18](000;000;000)(1.0,0.088388,0.994369)-10.23343694 2.79250应用数学2021表5.2算例7利用本文算法产生的数值结果m n Iters time(s)m n Iters Time(s)5584 2.173******** 6.825173105121 3.430311572787.80387330594 2.72547758202856.29689350563 1.936146595969166.40530760567 2.16155251013957430.4662805148 4.579592106234 6.701242100576 2.49782210780422.848295200537 1.42267110107331215.5171365005138 3.66477310208124123.487181表5.3算例8利用本文算法产生的数值结果n Refs Optimal value MatrixρIters Time(s)5本文-25(0)5×510.155[17]-25-115 3.82565[18]-25(0)5×5186 6.3310本文-100(0)10×1010.310[17]-100-47016.40610[18]-100(0)10×1072627.820本文-400(0)20×2010.6320[17]-400-193074.6547620[18]-400(0)20×202961140.2630本文-900(0)30×3010.8530[18]-900(0)30×306541245.939450本文-1500(0)50×501 1.38750[18]-1500(0)50×5019176931.588100本文-10000(0)100×1001 3.345200本文-40000(0)200×2001 6.9545为方便起见,表5.1-5.3的具体符号表示如下:E:算例编号;Refs:文献;Matrixρ:对称参数矩阵;f0(x∗):最优值;x∗:最优解;m:二次约束的个数;n:目标函数的维数;Iters:平均迭代次数;Time(s):平均运行时间,算法运行过程中的精度为10−6.算例1-6的计算结果如表5.1所示,利用文[17-18]和本算法进行比较,从总体上看本文要优于文[17-18]的计算效果,而从算例3和算例5的结果来看,本文计算时间要比文[18]长,但比文[17]用时短,证明了我们提出的算法是可行有效的.为了进一步证明我们的算法能够解决带有多个二次约束的二次规划问题,对算例7进行了如表5.2结果所示的中高维度的伪随机试验,在表5.2中,我们把参数矩阵ρ∈R n×n固定成ρ=(0)n×n.当约束条件个数固定,维数逐渐增大的情况下,迭代次数没有明显的增加,且运行时间较少;当维数固定,二次约束的个数不断增加时,迭代次数增大且运行时间越长.表5.3中,列出了我们的算法和文[17-18]中算法的数值结果,参数矩阵ρ的元素是全为0的对称矩阵,通过比较数值结果,不管从迭代次数还是从CPU运行时间,都可以看出我们的算法是优于文[17-18]的,随着维数的增加,运行时间在增加,但迭代次数只有1次,说明本算法对算例8的维数大小不是很敏感.数值结果表明我们的算法是有效的.6.结论针对一般的二次约束二次规划问题,本文基于参数化方法求解关于二次项以及双变量的第1期黄小利等：一种新的二次约束二次规划问题的分支定界算法251线性松弛问题.为了加快算法的收敛速度,利用参数化线性松弛分支定界算法,以及结合区域缩减技术,解决了二次约束二次规划问题.通过剖分矩形区域并解决参数化线性子问题,提出的算法最终收敛到原问题(QCQP)的全局最优解,在第五部分数值实验中,验证了在求解二次约束二次规划问题时本文算法是有效可行的.参考文献：[1]LOIOLA E M,ABREU N M M D,BOAVENTURA-NETTO P O,et al.A survey for the quadraticassignment problem[J].European Journal of Operational Research,2007,176(2):657-690.[2]KONAR A,SIDITOPOULOS N,MEHANNA O.Parametric frugal sensing of power spectra formoving average models.IEEE Transactions on Signal Processing.,2015,63(5):1073-1085.[3]SIDITOPOULOS N D,DAVIDSON T N,LUO Z Q.Transmit beamforming for physical-layer mul-ticasting[J].IEEE Transactions on Signal Processing,2006,54(6):2239-2251.[4]MATSKANI E,SIDITOPOULOS N D,LUO Z Q,et al.Convex approximation techniques for jointmultiuser downlink beamforming and admission Control[J].IEEE Transactions on Wireless Commu-nications,2008,7(7):2682-2693.[5]PARDALOS P M,VAVASIS S A.Quadratic programming with one negative eigenvalue is NP-hard[J].Journal of Global Optimization,1991,1(1):15-22[6]DU D Z,PARDALOS P M,WU W L.Reformulation-linearization methods for global optimization[J].Encyclopedia of Optimization,2000:2182-2186.[7]GAO Y,WEI F.A new bound-and-reduce approach of nonconvex quadratic programming problem-s[J].Applied Mathematics and Computation,2015,250:298-308.[8]YAMADA S,TAKEDA A.Successive lagrangian relaxation algorithm for nonconvex quadratic opti-mization[J].Journal of Global Optimization,2018,71(2):313-339.[9]KIM S,KOJIMA M.Second order cone programming relaxation of nonconvex quadratic optimizationproblems[J].Optimization Methods and Software,2001,15(3-4):201-224.[10]LUO Z,MA W,SO A M,YE Y,ZHANG S.Semidefifinite relaxation of quadratic optimizationproblems[J].IEEE Signal Processing Magazine,2010,27(3):20-34.[11]BURER S,VANDENBUSSCHE D.Globally solving box-constrained nonconvex quadratic programswith semidefinite-basedfinite branch-and-bound[J].Computational Optimization and Applications, 2009,43(2):181-195.[12]CHEN J,BURER S.Globally solving nonconvex quadratic programming problems via completelypositive programming[J].Mathematical Programming Computation,2012,4(1):33-52.[13]TSENG P.Further results on approximating nonconvex quadratic optimization by semidefinite pro-gramming relaxation[J].SIAM Journal on Optimization,2001,14(1):268-283.[14]KEYANPOUR M,OSMANPOUR N.On solving quadratically constrained quadratic programmingproblem with one non-convex constraint[J].OPSEARCH,2018,55(2):320-336.[15]JIANG R,LI D.Second order cone constrained convex relaxations for nonconvex quadratically con-strained quadratic programming[J].Journal of Global Optimization,2019,75(2):461-494.[16]LU C,DENG Z,ZHOU J,et al.A sensitive-eigenvector based global algorithm for quadraticallyconstrained quadratic programming[J].Journal of Global Optimization,2019,73(2):371-388. [17]JIAO H,CHEN Y.A global optimization algorithm for generalized quadratic programming[J].Journalof Applied Mathematics,2013,2013:1-9.[18]JIAO H,LIU S,LU N.A parametric linear relaxation algorithm for globally solving nonconvexquadratic programming[J].Applied Mathematics and Computation,2015,250:973-985.[19]QU S J,ZHANG K C,JI Y.A global optimization algorithm using parametric linearization relax-ation[J].Applied Mathematics and Computation,2007,186(1):763-771.252应用数学2021[20]WU H,ZHANG K C.A new accelerating method for global non-convex quadratic optimization withnon-convex quadratic constraints[J].Applied Mathematics and Computation,2008,197(2):810-818.[21]AUDET C,HANSEN P,JAUMARD B,et al.A branch and cut algorithm for nonconvex quadraticallyconstrained quadratic programming[J].Mathematical Programming,2000,87(1):131-152.[22]LUO H,BAI X,PENG J.Enhancing semidefinite relaxation for quadratically constrained quadraticprogramming via penalty methods[J].Journal of Optimization Theory and Applications,2019,180(3): 964-992.A New Branch and Bound Algorithm for QuadraticallyConstrained Quadratic Programming ProblemHUANG Xiaoli1,2,GAO Yuelin1,2,XIE Jinxiao1,2,GU Jianfeng1,2(1.School of Mathematics and Information Science,North Minzu University,Yinchuan 750021,China;2.Ningxia Scientific Computing and Intelligent Information ProcessingCo-Innovation Center,Yinchuan750021,China)Abstract:In this paper,in order to obtain globally optimal solution of the quadratically con-strained quadratic programming problem(QCQP),a new parametric linearized relaxation branch and bound algorithm is proposed.Our algorithm utilize the parametric linearized relaxation technique to get the lower bound of the global minimum value of(QCQP),using the region-reduction-technique to delete the infeasible region as much as possible in which can accelerate the convergence speed of the algorithm. The numerical experiments show that the proposed algorithm is effective and feasible.Key words:Quadratically constrained quadratic programming;Global optimization;Branch and bound;Parametric linearized relaxation;Region-reduction。
