基于梯度投影的廣義濾子填充函數方法

張慧雯,王 薇,李 民,徐以汎

(1.華東理工大學數學系,上海 200237)

(2.復旦大學管理學院,上海 200433)

1 引言與假設

本文討論如下帶線性約束的全局優化問題,

其中f(x):Rn→R為非凸函數,A為m×n階矩陣,b=(b1,b2,b3,···,bm)T∈Rm.對約束問題(1.1),稱所有滿足約束的點為可行點,由所有可行點組成的集合X={x∈Rn:Ax≤b}為可行域.一般說來,問題(1.1)有多個極小值,用傳統方法求解只能得到局部最優,而且對于全局最優解的判定條件至今都缺少實用的研究成果,這都使得問題求解仍然面臨著許多困難.

填充函數法是求解多極值最優化問題的常用算法之一,其主要思想是:在求得問題的一個局部極小點后,構造填充函數.通過極小化該填充函數,尋找比當前局部極小點更優的點,進而得到更優的極小值[1?5].與其他方法相比,該算法思想簡單,而且僅用到經典算法,所以易于實現且效率較高,同時也可以推廣到其他非線性問題以及離散數學規劃的求解中[6].

濾子技術最早由Fletcher和Leyffer提出,他們詳細討論了濾子作為代替罰函數的工具在局部優化算法中的一些應用[7,8].之后濾子技術在局部優化問題的求解中被認為是一種更有效的方法,因其良好的數值效果,許多學者繼續進行了一系列的相關研究[9?11].

梯度投影算法自從被Rosen[12]提出后就引起了廣泛的注意和系統的研究[13,14],由于該方法簡單、實際應用的數值效果好,在一些更有效的近代算法中也繼續沿用了它的基本思想[15,16].本文為了優化填充函數算法,將濾子技術和改進的梯度投影方法融入到全局優化算法中,提出了基于梯度投影的廣義濾子填充函數算法求解問題(1.1).

為方便起見,下面做一些記號說明.

J={1,2,···,m}為指標集,A的第j行為aTj,記cj(x)=aTjx?bj,j∈J.

記約束違反度函數h(x)=max(0,cj(x),j∈J),A(x)=(aj,j∈J0(x)),其中J0(x)={j|cj(x)≥0,j∈J},并且記.

L(P)表示問題(1.1)的局部最優解集合,G(P)表示問題(1.1)的全局最優解集合.

x?是f(x)的一個穩定點,S1(x?)={x|f(x)≥f(x?),x∈X{x?}}稱為高水平集,S2(x?)={x|f(x)

intX表示可行域X的內點集合,?X表示可行域X的邊界點集合.

考慮問題(1.1),我們首先提出以下假設:

[A1]f(x)只有有限多個極小值,即存在直徑充分大的閉箱?能包含所有極小值.

[A2]?f(x)在?上連續.

[A3]?x∈Rn,{aj,j∈J0(x)}為線性無關向量組.T

由A1和A2可知,原問題的全局解必在有界閉區域X?上,因此可以認為X是有界閉集,且算法中的x總是取自?.

2 廣義填充函數

定義2.1設x?是問題(1.1)當前的局部極小值,稱T(x,x?)是f(x)在x?處的廣義填充函數,如果T(x,x?)滿足

(i)T(x,x?)在高水平集S1(x?)上沒有穩定點;

(ii)如果x?不是問題(1.1)的一個全局極小值點,即x?/∈G(P),那么存在點,使得是T(x,x?)的極小值點.

下面給出求解問題(1.1)的廣義填充函數.

設x?是問題(1.1)的一個局部極小點,定義單參數廣義填充函數

其中

根據[A2],顯然T(x,x?,r)在可行域上是連續可微的.下面來討論T(x,x?,r)的性質.

定理2.1對充分小的r>0,有如下結論成立

(i)函數T(x,x?,r)在集合S1(x?)?X上沒有穩定點;

(ii)當x?x?與{aj,j∈J0(x)}線性無關時,函數T(x,x?,r)在集合S1(x?)T?X上沒有穩定點.

證(i)?x∈S1(x?)?X,有

根據A2,f(x)≥f(x?),則T(x,x?,r)>0.所以當r>0充分小時,有在可行域內是有界的.而由x∈S1(x?)?X可知

即函數T(x,x?,r)在集合S1(x?)?X上沒有穩定點.

成立,其中λj≥0,λjcj(xL)=0,j∈J.而由(i)的證明可知

且當r>0充分小時,顯然上式右端第一項趨于0.則由(2.4),(2.5)式知

這與已知條件矛盾,故函數T(x,x?,r)在集合上沒有穩定點.

由(2.3)式,顯然有

推論2.1x?x?是函數T(x,x?,r)在S1(x?)?X處的一個嚴格下降方向.

定理2.2如果x?∈L(P),但是x?G(P),那么T(x,x?,r)在S2(x?)上至少有一個極小值點.

證由于G(P)非空和(2.2)式的定義,則必有xG∈G(P),使得T(xG,x?,r)<0.而T(x,x?,r)在X上連續,所以必存在函數的最小值點,使得

也就是說1?e?r12[f(x)?f(x?)+r]<0,即.

推論2.2如果x∈G(P),則對充分小的r>0,函數T(x,x?,r)在可行域內部intX沒有穩定點.

由定理2.1和定理2.2可知,T(x,x?,r)是一個廣義填充函數.

3 濾子和梯度投影

濾子最初定義為由兩個相互競爭的函數φ(x)和ψ(x)組成的數對集合,記為(φ(x),ψ(x)).為了之后討論方便,本小節將給定的一階連續可微函數Q(x)作為目標函數,即考慮如下問題

取φ(x)為Q(x),ψ(x)為h(x),下面給出“支配”的相關概念.

定義3.1點xk支配點xl當且僅當Q(xk)≤Q(xl)且h(xk)≤h(xl).

定義3.2濾子是由一列互不支配的數對組成,即若點xk和xl都在濾子中,那么Q(xk)≤Q(xl)和h(xk)≤h(xl)兩者必有之一不成立.

此外,為了保證算法的下降性,在本文中如果說點xk+1“被濾子接受”當且僅當存在濾子中的點xl,使得下面兩個不等式

之一成立,其中β,η∈(0,1).

從以上概念可以看出,濾子能夠作為一個評判標準來決定對當前試探步(沿下降方向尋找合適步長時假設的迭代點)是否接受.但是,以(3.2)和(3.3)式作為濾子集的接受標準可能導致算法收斂到一個可行的非穩定點.因此,本文在下降搜索時又另外采用了其他準則.

對當前迭代點xk,記

其中I為n階單位矩陣,稱P(xk)為xk處的投影矩陣,并記xk處的搜索方向為

其中

對于當前點的試探步長αk,本文要求一個足夠下降量標準(轉換條件)

成立.固定參數為δ1>0,s2>1,s1>2s2,其中mk(αk):=αk?Q(xk)Tdk.

若轉換條件(3.5)成立,則可以不再受濾子接受準則約束,只要求目標函數的Armijo條件

有時算法在搜索過程中目標函數逐步下降但迭代點卻離可行域越來越遠,此時就需要進入可行性恢復階段.下面簡述下可行性恢復的具體過程.

設點x∈?但,給定一個固定的參數ε>0,不妨設有,即存在k個約束不在可接受的范圍內,需要進行可行性恢復,則求解問題

對于進入可行性恢復階段的判定,本文采用如下準則.

如果當前迭代步αk足夠大,轉換條件(3.5)對于某個α≤αk是成立的,那么仍有可能存在某個小一點的步長能被Armijo條件(3.6)接受,此時就不用進入可行性恢復階段.故由(3.5)式可知,若?Q(xk)Tdk<0,只要,就不進入可行性恢復階段.

然而,當轉換條件(3.5)不成立時,考慮如下兩個線性近似

4 基于梯度投影的廣義濾子填充函數算法及其性質

本文研究的是帶線性約束的全局優化問題,同時從目標函數、填充函數和約束違反度函數三個角度考慮,所以將一般濾子推廣到三維,構成濾子(f(x),T(x,x?),h(x)).即

(i)點xk支配點xl當且僅當f(xk)≤f(xl),T(xk,x?)≤T(xl,x?)和h(xk)≤h(xl)同時成立.

(ii)若點xk和點xl都在濾子中,那么在f(xk)≤f(xl),T(xk,x?)≤T(xl,x?)和h(xk)≤h(xl)中至少有一個不等式不成立.

在本文中,用Fk表示點{xl|l≤k}的集合,也就是說(f(xl),T(xl,x?),h(xl))是當前濾子中的元素,Fk為當前濾子.另外,|F|表示濾子F所包含的元素個數.

此外,如果說點xk+1“被濾子Fk接受”當且僅當存在點xl∈Fk,使得下面三個不等式

之一成立,其中β1,β2,η ∈(0,1).

定義4.1將數組(f(x),T(x,x?),h(x))加入到當前濾子中,并把被(f(x),T(x,x?),h(x))支配的所有數組移出濾子的過程稱為更新濾子.即

根據三維濾子的定義可知,如果初始濾子集非空,那么在任意點處的濾子集非空,即|F|≥1.在提出廣義濾子填充函數算法之前,我們先給出試探步長αk的迭代算法.

回溯線搜索算法

步驟0令αk=1.

步驟1若時,轉步驟5;否則,計算下一迭代xk(αk)=xk+αkdk.

步驟2若濾子中存在數組能支配xk(αk),則拒絕試探步,轉步驟4;否則,進入步驟3.

步驟3檢驗足夠下降量

3.1情況I若xk∈X時,(3.5)和(3.6)式同時成立,并且xk(αk)∈X,則αk為所求步長;否則,轉步驟4.

3.2情況II若xkX且(3.5)式成立時,(3.6)式成立,則αk為所求步長;否則,轉步驟4.

3.3情況III若xkX且(3.5)式不成立時,試探點能被濾子接受,則αk為所求步長;否則,轉步驟4.

步驟4選擇新的試探步長,轉步驟2.

步驟5可行性恢復:求解問題(3.7),并更新濾子集.

本算法的目的是對當前迭代點xk和方向dk,尋找合適的下降步長.當xk不需要進行可行性恢復且下一步迭代xk(αk)不被濾子支配時,試探步長αk還需要滿足足夠下降量或者濾子接受準則.算法中對足夠下降量的檢驗作了詳細的情況分析,確保當前點無論是否在可行域內一定存在可被接受的步長.

下面給出主算法.

算法

步驟0任取x0∈Rn,給定精度參數ε>0,鄰域半徑δ>0,r0,r,G,N為變量x的維數.令k:=0,若當前點沒有填充函數值,則令其分量T=Z,Z充分大.初始化濾子集F0=(f(x0),Z,h(x0)).

步驟1取Q(x)為f(x),若有dk=?P(xk)?f(xk)=0,當U(xk)≥0和h(xk)=0同時成立,令x?=xk,轉步驟5;當存在uj(xk)<0,j∈J0(x)時,轉步驟4;否則,進入步驟2.

步驟2回溯線搜索,若因可行性恢復跳出,轉步驟1;否則,進入步驟3.

步驟3令xk+1=xk(αk),k←k+1,更新濾子,轉步驟1.

步驟4令J0(xk){j},重新構造A(xk),轉步驟1.

步驟5在x?處構造填充函數T(x,x?,r),令S=1.

步驟6若S>2N,則令,轉步驟 12;否則,令k:=0,取xk∈O(x?,δ){x?}.

步驟7取Q(x)為T(x),若dk=0,令F={x?},S=S+1,轉步驟6;否則進入步驟8.

步驟8回溯線搜索,若因可行性恢復跳出,轉步驟11;否則,進入步驟9.

步驟9令xk+1=xk(αk),k←k+1,更新濾子.

步驟10若|F|=1,則令k:=0,轉步驟1;否則,進入步驟7.

步驟11若|F|>G,則令F={x?},S=S+1,轉步驟6;否則,轉步驟7.

步驟12若r

下面討論該算法的相關性質.

引理4.1矩陣P(x)是半正定矩陣,且滿足P(x)2=P(x).

定理4.1(i)當xk∈X且非KKT點時,由算法得到的dk滿足

(ii)當xkX時,由算法得到的dk滿足.

證(i)顯然,當xk∈X且非KKT點時,有

另外

(ii)當xkX時,

定理4.2回溯線搜索有限終止.

證(i)若xk∈X時,.由定理4.1知在非KKT點處,?Q(xk)Tdk<0,則必有mk(αk)<0,并且,即 (3.5) 式成立.

另一方面,當αk足夠小時,有

即(3.6)式成立.因此若xk∈X時,算法必能找到合適的αk被接受.

綜上,算法必能找到合適的αk,回溯線搜索總會有限終止.

為了證明算法性質,本文給出如下假設

[A4]存在Mm>0,使得|mk(α)|≤Mmα.

定理4.3若假設都成立,則有

證(i)若濾子集僅在有限步擴大,假設K,當迭代k>K時,濾子集不再擴充.則由算法可知,對所有k≥K,(3.5)和(3.6)式都成立,則

故對i=1,2,···,

由Q(xK+i)下有界可知,(4.5)式成立.

(ii)若濾子集一直擴大,下降搜索時所得序列{xKi},令{xki}是其子序列,濾子集在ki∈{Ki}擴大,對所有i,不妨設存在常數C1∈R和C2>0,使得Q(xki)≥C1,h(xki)≤C2.

下面驗證濾子的相關性質.

定理4.4設xk∈S1(x?),當前濾子集為Fk.若r>0充分小,則一定存在α>0,使得xk+1=xk+αd被濾子Fk接受.

證由定理4.1可知,當xk∈X且不是問題(3.1)的KKT點時,由算法得到的dk滿足

所以當xk∈S1(x?)時,算法中的迭代方向一定是f(xk)或T(xk,x?,r)的下降方向.根據算法的下降性以及濾子的定義,必有xl∈Fk,使得f(xk+1)

下面討論濾子集Fk和低水平集S2(x?)的關系.

定理4.5設xk處的濾子集為Fk={xl|1≤l≤k}.如果點xk+1被Fk接受,并且支配所有xl(?xl∈Fk),則xk+1∈S2(x?).

證若xk+1,則對,xk+1必然無法支配,故xk+1∈X.若存在,因為xk+1支配xl,結論顯然成立.否則,對所有,根據濾子的定義,點x?肯定在Fk內.而xk+1被Fk接受且支配x?,由不等式(4.1)–(4.3)可知f(xk+1)

定理4.6設xk處的濾子集為Fk={xl|1≤l≤k}?S1(x?),且S2(x?).則必存在點xN∈S2(x?)使得?xl∈Fk,xN支配xl.

證首先,?xN∈S2(x?),有.

另外,根據S1(x?)和S2(x?)的定義,xl和xN均在可行域內,則h(xl)=h(xN)=0,所以xN支配xl.

根據定理4.6,如果濾子集里只有一個元素,且該元素支配x?,則算法進入到x?的低水平集里.

定理4.7設xk處的濾子集為Fk={xl|1≤l≤k}?S2(x?),進一步搜索得x??∈L(P),則x??必能支配Fk中所有的點.

證由得

則當r充分小時,T(x??,x?,r)

5 數值結果

本節主要驗證基于梯度投影的廣義濾子填充函數算法的有效性.這些算例都是在Matlab2014b運行環境下運行,處理器為Intel(R)Core(TM)i5-2410M CPU@2.30GHZ 2.30GHZ,系統類型為64位操作系統,基于x64的處理器.算法的參數設置如下:s1=2.5,s2=1.2,δ1=δ2=β1=β2=η=10?6,r0=1,r=10?3,G=500,δ=10?3.

算法的數值結果在表1–4中列出,其中各個符號定義如下.

k:第k次求解得局部極小點;:進行第k次局部極小點迭代的初始點;:第k個局部極小點;:第k個局部極小值;Fk:在第k次達到局部極小點時的濾子.

算例5.1

全局極小值為f(x∞)=0.42196,總運行時間為8.513秒.

算例5.2

表1:算例5.1

表2:算例5.2

全局極小值為f(x∞)=?147.26943,總運行時間為4.566秒.

算例5.3

表3:算例5.3

表4:算例5.4

算例5.4

取在可行域外的初始點

得到全局極小點

全局極小值為f(x∞)=0.55285,總運行時間為409.516秒.由于迭代得的局部極小點比較多,在本文中只列出部分迭代結果.

以上的數值結果說明了基于梯度投影的廣義濾子填充函數算法的有效性.算例5.1與算例5.2的最優解在內部,其中算例5.1首先搜索到邊界上的KKT點,再通過填充函數迭代到低水平集,并最終找到了目標函數的穩定點.算例5.2首先在可行域內找到穩定點,再通過兩階段迭代找到最優解.算例5.3的最優解在邊界上,其搜索情況與算例5.1類似.算例5.4的情況較為復雜,首先從可行域外搜索到了KKT點,然后在多次迭代后達到最優解.前三個算例由于規模較小,耗時均較短.而算例5.4的維數較高,程序的總運行時間相對長了很多.