二次約束二次規(guī)劃問題的二元均值松弛定界算法

田福平, 高岳林,2, 孫 瀅

(1.北方民族大學(xué) 數(shù)學(xué)與信息科學(xué)學(xué)院,寧夏 銀川 750021; 2.合肥工業(yè)大學(xué) 數(shù)學(xué)學(xué)院,安徽 合肥 230601; 3.寧夏科學(xué)計(jì)算與智能信息處理協(xié)同創(chuàng)新中心,寧夏 銀川 750021)

本文主要考慮的二次約束二次規(guī)劃((quadratically constrained quadratic programming,QQP)問題形式如下:

其中:Qs=(qsij)n×n為實(shí)對稱矩陣;cs∈Rn;ds∈R,s=0,1,2,…,N;i,j=1,2,…,n;A∈Rm×n;b∈R。QQP問題的可行域記為F,假設(shè)F是非空有界閉集。

一般二次約束二次規(guī)劃問題是非線性規(guī)劃(nonlinear programming,NLP)問題,它在很多領(lǐng)域都有廣泛的應(yīng)用,如無線通信、網(wǎng)絡(luò)、雷達(dá)、信號處理等[1-2]。當(dāng)QQP問題中的Q0,Q1,Q2,…,Qs均為半正定矩陣時,可以利用文獻(xiàn)[3]中的二階錐規(guī)劃方法對其進(jìn)行求解;另一方面,一些特殊的非凸QQP問題可以等價地轉(zhuǎn)化為凸問題[4-6]。但是,一般二次約束二次規(guī)劃問題通常是NP難問題,因此很難求到此類問題的全局最優(yōu)解。

學(xué)者們提出了各種各樣的用于求解非凸二次規(guī)劃問題的分支定界算法[7-8]。同時,也有一些可以用來求解非凸QQP問題全局最優(yōu)解的求解器,如SCIP、GloMIQO、Ipopt、ANTIGONE、BARON、COUENNE等。這些算法和求解器已經(jīng)在某些類型的問題上有了十分良好的性能,如BARON、SCIP、COUENNE求解器對于求解具有稀疏參數(shù)的線性約束二次規(guī)劃問題是非常有效的;文獻(xiàn)[9]提出的算法可以有效地求解帶有盒約束的二次規(guī)劃問題;文獻(xiàn)[7]提出的算法專門用于求解帶有可分離約束的混合整數(shù)二次規(guī)劃問題。但是對于一般的QQP問題,目前尚沒有較好的求解方法。在利用分支定界算法求解極小化問題時最關(guān)鍵的是定下界,能否成功構(gòu)造出最優(yōu)值的下界序列對算法設(shè)計(jì)成功與否是至關(guān)重要的,定下界過程要求把下界問題轉(zhuǎn)化為一系列盡可能簡單并易于求解的松弛子問題。文獻(xiàn)[10]將松弛方法分為線性松弛、凸包絡(luò)松弛、拉格朗日對偶松弛、二階錐規(guī)劃松弛和半定規(guī)劃松弛。在這些不同類型的松弛中,目前求解全局最優(yōu)解時最常用的方法是線性松弛。

本文通過引入輔助乘積變量將原QQP問題等價地轉(zhuǎn)化為NLP問題,基于二元均值不等式將等價非線性規(guī)劃問題轉(zhuǎn)化為松弛線性規(guī)劃(relaxation linear pngramming,RLP)問題提出了一種新的剖分策略,給出了超矩形縮減的方法及其算法的描述,并證明了該算法的收斂性,通過數(shù)值實(shí)驗(yàn)驗(yàn)證了本文提出的算法是可行有效的。

1 基于輔助變量的等價非線性規(guī)劃

本文通過求解以下2n個線性規(guī)劃得到x的上下界:

(1)

(2)

由(1)式得最優(yōu)值lj,j=1,2,…,n;由(2)式得最優(yōu)值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)轉(zhuǎn)化為以下等價問題:

其中,·表示2個矩陣的內(nèi)積,即2個矩陣對應(yīng)元素乘積的和。

因?yàn)镼0=(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,結(jié)合線性函數(shù)來構(gòu)造QQP(Hk)問題的松弛線性規(guī)劃,該松弛線性規(guī)劃為問題QQP(Hk)提供了一個最優(yōu)值的下界,具體操作如下:

對于每一個i={1,2,…,n},計(jì)算

(3)

(4)

一方面,對于每一個i=1,2,…,n都有l(wèi)i≤xi≤ui,由xili≥0,xjlj≥0，以及二元均值不等式可知:

(5)

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

xixj≤(ljxi+lixj-lilj)+

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

(6)

另一方面,設(shè)函數(shù)G(x)=x2,x∈R過點(diǎn)(l,l2)和(u,u2),得到關(guān)于x的直線函數(shù):

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

當(dāng)x=xi時,U(xi)=(li+ui)xiliui;當(dāng)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)

基于二次均值不等式,結(jié)合函數(shù)的性質(zhì),可以得到松弛線性規(guī)劃RLP(Hk)如下:

以下引理說明了問題QQP(Hk)、NLP(Hk)、RLP(Hk)的可行解以及最優(yōu)解之間的關(guān)系。

引理2 若(Yk,xk)為問題RLP(Hk)的可行解,并且滿足Yk=xk(xk)T,則xk為QQP(Hk)的可行解。

引理3 若(Yk,xk)為問題RLP(Hk)的最優(yōu)解,并且滿足Yk=xk(xk)T,則(Yk,xk)為問題NLP(Hk)的全局最優(yōu)解,從而xk為QQP(Hk)的全局最優(yōu)解。

證明設(shè)問題RLP(Hk)的最優(yōu)值為α(Hk),問題QQP(Hk)的全局最優(yōu)值為Z(Hk)。

因?yàn)?Yk,xk)為RLP(Hk)的最優(yōu)解,所以α(Hk)≤Z(Hk),并且α(Hk)≤f0(Yk,xk)。因?yàn)閥ijk=xikxjk,所以由引理2知,xk∈F,從而Z(xk)≤f0(xk),并且f0(xk)=f0(Yk,xk),因此α(Hk)=Z(Hk)=f0(Hk),即xk為QQP(Hk)的全局最優(yōu)解。

定義1 若(Yk,xk)為RLP(Hk)的解,其中xk∈F并且fs(xk)-fs(Yk,xk)≤ε,s=0,1,2,…,N,則xk稱為問題QQP(Hk)的強(qiáng)ε-全局最優(yōu)解。

定義2 若(Yk,xk)為RLP(Hk)的解,其中xk∈F,并且f0(xk)-f0(Yk,xk)≤ε,s=0,1,2,…,N的ε-全局最優(yōu)解。

3 超矩形的剖分與縮減

3.1 新的分支方法

為了適當(dāng)?shù)靥岣咚惴ǖ姆种?本文采用了一種新的剖分方法,令Hk=[lk,uk]∈H表示當(dāng)前所要剖分的超矩形,在給出具體的剖分方法之前,首先給出以下說明:對所有的變量xi,其中l(wèi)i≤xi≤ui,當(dāng)xi∈Hk=[lk,uk]時,拋物線g(xi)=(xi)2與直線l(xi)=(lik+uik)xilikuik(i=1,2,…,n)圍成的面積記為S(lik,uik),如圖1所示。顯然

圖1 關(guān)于變量xi拋物線直線圍成的面積

用xk表示問題QQP(Hk)的當(dāng)前最優(yōu)解,對應(yīng)的松弛問題的最優(yōu)解為(Yk,xk),顯然,xk∈D∩Hk=[lk,uk],進(jìn)行以下形式的剖分:

else找到第1個xidk∈argmaxω

end if;

通過x′與x″的連線或連線所在的平面將超矩形Hk剖分為Hk1=[lk1,uk1]和Hk2=[lk2,uk2]2個子超矩形,則這2個子超矩形分別為:

3.2 超矩形的縮減

為了提高算法的收斂速度,本文采用超矩形縮減技術(shù)[11],假設(shè)RLP中所有的線性不等式約束的展開形式為:

具體縮減方法描述如下:

令I(lǐng)k∶={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;

為了方便起見,本文把該方法產(chǎn)生的新的整超矩形仍然記為Hk,它是原來整超矩形的子集。

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

首先,為了方便算法的敘述,假定算法迭代到第k步時,本文給出以下記號:F表示QQP問題的可行域;Q表示QQP問題的當(dāng)前所有可行點(diǎn)的集合;Hk表示當(dāng)前迭代步所要剖分的超矩形;Ω表示剪枝后剩下的所有超矩形的集合;U表示QQP問題當(dāng)前全局最優(yōu)值的上界;L表示QQP問題當(dāng)前全局最優(yōu)值的下界。

基于上述記號,下面給出求解QQP問題的分支定界算法具體描述:

(1) 初始化。通過 (3) 式、(4) 式構(gòu)造關(guān)于x的初始超矩形H=H0=[l0,u0],在超矩形H0上求解RLP問題;如果其對應(yīng)的最優(yōu)解和最優(yōu)值分別記為(Y0、x0)、L;那么L就是問題QQP當(dāng)前全局最優(yōu)值的一個下界,如果x0∈F,令Q=Q∪x0,初始上界為U=min{f(x):x∈Q},并找到QQP問題當(dāng)前的一個最好的解,那么給定所需容忍度ε>0,Ω=H0;置k=1,轉(zhuǎn)入步驟(2)。

(2) 終止準(zhǔn)則。若U-L≤ε或Ω≠φ,則迭代終止,輸出QQP問題當(dāng)前全局最優(yōu)解x*以及全局最優(yōu)值f(x*);否則轉(zhuǎn)步驟(3)。

(3) 選擇規(guī)則。在Ω中選擇L所對應(yīng)的超矩形Hk,即L=L(Hk)，此時Ω=ΩHk。

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

(5) 超矩形的縮減。運(yùn)用3.2中超矩形的縮減技術(shù),對剖分后的子超矩形進(jìn)行縮減和刪除,把縮減后的超矩形仍然記為Ω=Hki,i={1,2}。

(6) 剪枝規(guī)則。讓Ω=Ω{H:L(H)≥U,H∈Ω}。

置k=k+1,轉(zhuǎn)到步驟(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,當(dāng)εi→0時,都有|φ-yij|→0。

證明令不等式(8)的右側(cè)[(2lj+ui-li)+(2li+uj-lj)xj-liui-ljuj+(li-lj)2]/2=φ,又因?yàn)?xili)(xjlj)≥0,所以有xixj≥ljxi+lixjlilj,因此

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

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

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

又因?yàn)棣舏→0,所以|φ-yij|→0。

定理2 令εi=ui-li,對于每一個i∈{1,2,…,n}以及任意的x∈F,當(dāng)εi→0時,都有maxS(li,ui)→0。

證明

maxS(li,ui)=

因?yàn)閡i-li→0,所以maxS(li,ui)→0,從而對于每一個i都有S(li,ui)→0。

定理4 設(shè)QQP問題的可行域F是n維的,所要剖分的關(guān)于變量x的超矩形是n維的,且超矩形的剖分是耗盡的,若算法在有限步k終止,則xk必是問題的全局最優(yōu)解,或者所產(chǎn)生的無窮可行點(diǎn)列{xk}的任一聚點(diǎn)都是QQP問題的全局最優(yōu)解。

證明(1) 若算法是有限步終止且算法終止時迭代了k步,則根據(jù)終止準(zhǔn)則可知Uk-Lk≤ε,即

f(xk)-Lk≤ε

(9)

假設(shè)此時的全局最優(yōu)解為x*,可知：

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

(10)

因此,將 (9)式、(10) 式聯(lián)立可得:

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

(11)

故(1)得證。

(2) 如果算法的迭代是無限的,QQP問題的可行解序列為{xk},對應(yīng)的線性松弛問題的可行解序列為{(xk,Yk)}。

根據(jù)步驟(4)、步驟(5),隨著x的超矩形的不斷加(x*)細(xì),則有:

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

k=1,2,…

(12)

因?yàn)橄陆缧蛄衶Lk=f0(xk,Yk)}是單調(diào)不減且有界的，上界序列{Uk=f0(xk)}是單調(diào)不增且有界,所以原問題的上下界都是收斂的。對(12)式兩邊取極限可得:

(13)

L=f(x*)=U

(14)

因此,序列{xk}的每一個聚點(diǎn)x*都是問題QQP的全局最優(yōu)解。

5 數(shù)值實(shí)驗(yàn)

本文給出幾個算例來證明本文算法的有效性。算法所有測試過程均用Matlab9.0.0.341360(R2015a),在Inter(R)Core (TM)i7-6700,CPU@3.40 GHz,16 GB內(nèi)存,64位Windows7操作系統(tǒng)的計(jì)算機(jī)上運(yùn)行，并且文中所有線性規(guī)劃均用對偶單純形法求解。

算例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]上隨機(jī)生成,本文利用特征值分解可以得到Q0=PTD0P及Qs=PTDsP,從而得到對角矩陣D0和Ds,其中D0的前r個分量在[-10,0]上隨機(jī)生成；其余對角分量在[0,10]上隨機(jī)生成；Ds的對角分量在[1,100]上隨機(jī)生成；c0為n維零向量；cs為在[-100,100]上隨機(jī)生成的n維向量；ds為在[1,50]上隨機(jī)生成的數(shù)。

算例1～算例10的計(jì)算結(jié)果見表1所列。表1中,x*為原問題的最優(yōu)解;f(x*)為目標(biāo)函數(shù)的最優(yōu)值;I為迭代次數(shù);ε為容忍度。

表1 算例1～算例10的計(jì)算結(jié)果

從表1可以看出,對于算例2、3，從迭代次數(shù)上看本文算法均優(yōu)于其他算法,計(jì)算結(jié)果與其他算法相同;針對算例4、10，從計(jì)算結(jié)果可知本文算法分別要優(yōu)于文獻(xiàn)[12,20]算法,迭代次數(shù)大于文獻(xiàn)[12,20]; 對于算例1,本文算法計(jì)算結(jié)果與其他算法相同,迭代次數(shù)大于其他算法;針對算例6、9,本文算法的迭代次數(shù)明顯小于其他算法,計(jì)算結(jié)果與其他算法相當(dāng);最后針對算例5、7、8,從迭代次數(shù)和計(jì)算結(jié)果上看本文的算法要弱于其他算法。

以上對算例1～算例10 的計(jì)算結(jié)果進(jìn)行分析，說明了本文算法是可行有效的。

為了進(jìn)一步說明本文算法的有效性,針對算例11,進(jìn)行了若干的隨機(jī)實(shí)驗(yàn),結(jié)果見表2所列,表2中，AI為平均迭代次數(shù);At為迭代終止時CPU平均運(yùn)行時間;m為二次約束的個數(shù);n為目標(biāo)函數(shù)的維數(shù);r為目標(biāo)函數(shù)矩陣中負(fù)特征值的個數(shù)。

表2中的容忍度均為5×10-4。通過表2的計(jì)算結(jié)果可以看出,當(dāng)n和r固定時,隨著二次約束m個數(shù)的增多,平均迭代次數(shù)AI和迭代時間At都在增加;若n和m固定,當(dāng)r

表2 隨機(jī)算例11的計(jì)算結(jié)果

6 結(jié) 論

針對二次約束二次規(guī)劃問題,本文利用二元均值不等式的性質(zhì)將等價問題中的非線性約束函數(shù)進(jìn)行了線性松弛，并結(jié)合所提出的剖分技術(shù)以及分支定界框架構(gòu)造了新的分支定界算法。數(shù)值實(shí)驗(yàn)結(jié)果驗(yàn)證了本文算法的有效性和可行性。