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

2021-02-05 09:30:20田福平高岳林

合肥工業大學學報(自然科學版) 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