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

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

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

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

1 引言與假設

本文討論如下帶線性約束的全局優(yōu)化問題,

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

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

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

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

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

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

記約束違反度函數(shù)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)的局部最優(yōu)解集合,G(P)表示問題(1.1)的全局最優(yōu)解集合.

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

intX表示可行域X的內(nèi)點集合,?X表示可行域X的邊界點集合.

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

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

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

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

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

2 廣義填充函數(shù)

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

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

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

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

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

其中

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

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

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

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

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

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

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

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

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

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

由(2.3)式,顯然有

推論2.1x?x?是函數(shù)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上連續(xù),所以必存在函數(shù)的最小值點,使得

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

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

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

3 濾子和梯度投影

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

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

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

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

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

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

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

對當前迭代點xk,記

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

其中

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

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

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

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

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

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

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

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

4 基于梯度投影的廣義濾子填充函數(shù)算法及其性質(zhì)

本文研究的是帶線性約束的全局優(yōu)化問題,同時從目標函數(shù)、填充函數(shù)和約束違反度函數(shù)三個角度考慮,所以將一般濾子推廣到三維,構(gòu)成濾子(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所包含的元素個數(shù).

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

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

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

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

回溯線搜索算法

步驟0令αk=1.

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

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

步驟3檢驗足夠下降量

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

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

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

步驟4選擇新的試探步長,轉(zhuǎn)步驟2.

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

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

下面給出主算法.

算法

步驟0任取x0∈Rn,給定精度參數(shù)ε>0,鄰域半徑δ>0,r0,r,G,N為變量x的維數(shù).令k:=0,若當前點沒有填充函數(shù)值,則令其分量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,轉(zhuǎn)步驟5;當存在uj(xk)<0,j∈J0(x)時,轉(zhuǎn)步驟4;否則,進入步驟2.

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

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

步驟4令J0(xk){j},重新構(gòu)造A(xk),轉(zhuǎn)步驟1.

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

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

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

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

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

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

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

步驟12若r

下面討論該算法的相關性質(zhì).

引理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,回溯線搜索總會有限終止.

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

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

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

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

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

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

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

下面驗證濾子的相關性質(zhì).

定理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)的下降方向.根據(jù)算法的下降性以及濾子的定義,必有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,結(jié)論顯然成立.否則,對所有,根據(jù)濾子的定義,點x?肯定在Fk內(nèi).而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?),有.

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

根據(jù)定理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 數(shù)值結(jié)果

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

算法的數(shù)值結(jié)果在表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秒.由于迭代得的局部極小點比較多,在本文中只列出部分迭代結(jié)果.

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