具有超線性收斂性的非線性半定規(guī)劃的濾子算法

趙 奇, 張 燕

(江蘇科技大學(xué) 張家港校區(qū) 基礎(chǔ)教學(xué)部,張家港 215600)

非線性半定規(guī)劃問(wèn)題(nonlinear semidefinite programming,NLSDP)是目前優(yōu)化研究的一個(gè)熱點(diǎn),在工程學(xué)上，很多地方需要計(jì)算上述NLSDP問(wèn)題, 如穩(wěn)定控制反饋模型[1-2]、桁架設(shè)計(jì)優(yōu)化[3]、魯棒優(yōu)化問(wèn)題[4-5]等.此外,該問(wèn)題還是相關(guān)應(yīng)用問(wèn)題, 如線性矩陣不等式[6]、二階錐規(guī)劃問(wèn)題[7]的基礎(chǔ).關(guān)于NLSDP問(wèn)題, 已經(jīng)有了一些有效的算法,如文獻(xiàn)[8]中提出了懲罰/障礙方法, 解決了帶線性矩陣不等式的凸半定規(guī)劃問(wèn)題.文獻(xiàn)[9]中提出了增廣拉格朗日函數(shù)方法, 并設(shè)計(jì)了該算法的PENNON 程序來(lái)解決 NLSDP和雙線性矩陣不等式問(wèn)題.文獻(xiàn)[10]中提出了求解NLSDP的原始-對(duì)偶內(nèi)點(diǎn)算法, 并分析了算法的全局和局部收斂性.文獻(xiàn)[11]中將非線性規(guī)劃中的序列二次規(guī)劃方法推廣到NLSDP, 提出了序列二次半定規(guī)劃方法.文獻(xiàn)[12-14]中將非線性規(guī)劃的無(wú)懲罰型方法推廣到非線性半定規(guī)劃問(wèn)題上去, 提出了濾子算法和可行上界控制方法.文中主要討論一種針對(duì)NLSDP問(wèn)題的濾子算法.該方法的思想來(lái)源于文獻(xiàn)[15]中提出的解決非線性規(guī)劃問(wèn)題的線性搜索濾子算法,在非線性規(guī)劃中,已經(jīng)證明了該算法具有全局和局部收斂性質(zhì).文中分析非線性半定規(guī)劃濾子算法局部收斂性,在嚴(yán)格互補(bǔ)松弛條件和帶sigma項(xiàng)的強(qiáng)二階充分條件下, 證明了算法產(chǎn)生的序列具有局部超線性收斂性.文獻(xiàn)[16-17]中對(duì)半定規(guī)劃局部收斂性質(zhì)做了研究工作,然而,這些理論分析中默認(rèn)步長(zhǎng)因子為1.在實(shí)際計(jì)算中,步長(zhǎng)因子不可能總為1,因此Maratos效應(yīng)也會(huì)發(fā)生. Maratos效應(yīng)指出, 對(duì)某些效益函數(shù),超線性收斂步可能被拒絕,從而破壞算法的超線性收斂性,因此必須采用一些手段來(lái)克服它. 文中去掉了步長(zhǎng)因子為1的默認(rèn)假設(shè),采用一種針對(duì)非線性半定規(guī)劃的二階校正技術(shù), 利用該技術(shù)避免Maratos效應(yīng).需要指出的是由于半定約束的存在,以前非線性規(guī)劃的技巧不可能直接推廣.

求解二次半定規(guī)劃子問(wèn)題時(shí), 往往用一個(gè)正定陣直接替代海塞陣近似,計(jì)算效果并不好, 文中在求解二次半定規(guī)劃子問(wèn)題時(shí),從SDPT3軟件包[18]的計(jì)算原理出發(fā), 利用最優(yōu)性條件的比較給出了拉格朗日乘子的表達(dá)形式，并將其應(yīng)用于BFGS校正中. 最后, 通過(guò)一些實(shí)際問(wèn)題的數(shù)值實(shí)驗(yàn), 表明該算法是穩(wěn)定而有效的, 驗(yàn)證了二階校正步和BFGS公式的作用.

1 問(wèn)題與算法

考慮如下形式的非線性半定規(guī)劃問(wèn)題:

(1)

問(wèn)題(1)中f:Rn→R,h:Rn→Rp,G:Rn→Sm都是光滑函數(shù).Sm表示m×m實(shí)對(duì)稱矩陣的線性空間. 文中矩陣內(nèi)積定義為〈A,B〉=trace(AB),A,B∈Sm,Sm表示負(fù)半定偏序，即A?B當(dāng)且僅當(dāng)A-B為負(fù)半定矩陣.

在當(dāng)前迭代步, 求解如下形式的二次半定規(guī)劃子問(wèn)題:

(2)

式中：gk=f(xk)；Bk為含有問(wèn)題(1) 的二階導(dǎo)數(shù)信息的對(duì)稱矩陣；Dh(x)為h(x)的Jacobi矩陣,線性算子DG(x)定義如下:

?p∈Rn

其伴隨算子DG(x)T定義如下:

二階算子D2G(x):Rn×Rn→Sm定義如下:

再定義線性化預(yù)測(cè)下降量：

和約束違反度函數(shù)：

v(x)=‖h(x)‖1+(λ1(G(x)))+

式中，λ1(A)為矩陣A的最大特征值, 且規(guī)定(α)+=max{0,α}.

如果下面的條件成立, 即

mk(α)<0, [-mk(α)]sf[α]1-sf>δ[v(xk)]sv

(3)

則步長(zhǎng)因子α滿足條件:

f(xk+αdk)-f(xk)≤ηfmk(α)

(4)

式中，sv>1,sf≥1,ηf∈(0,0.5),稱滿足條件(3～4)的迭代為f-型迭代. 否則,要求α沿著搜索方向dk使得約束違反度改善或目標(biāo)函數(shù)下降兩者之一成立, 即下面兩個(gè)條件：

v(xk+αdk)≤(1-γv)v(xk)，或者,

f(xk+αdk)≤f(xk)-γfv(xk)

(5)

至少有一個(gè)成立,式中γv,γf∈(0,1).

文中濾子集合定義為一個(gè)集合F?[0,+∞)×R;該集合為一個(gè)拒絕域,包含所有的在k個(gè)迭代步不能被接受的(v,f)點(diǎn)對(duì).規(guī)定一個(gè)嘗試步xk(α)=xk+αdk能夠被濾子集合接受時(shí),其對(duì)應(yīng)的(v,f)點(diǎn)對(duì)不在拒絕域中,即(v(xk(α)),f(xk(α))?Fk.

在開(kāi)始時(shí)，定義F0={(v,f)∈R2:v>vmax},其中vmax>v(x0).在整個(gè)計(jì)算過(guò)程中, 當(dāng)新的迭代步xk+1被接受后,如果不是f-型迭代,則濾子集合采用下面的公式進(jìn)行更新:

Fk+1:=Fk∪{(v,f)∈R2:v≥(1-γv)v(xk),f≥f(xk)-γfv(xk)}

(6)

否則, 濾子集合保持不動(dòng),Fk+1:=Fk.

如果α太小,則認(rèn)為算法做的線性搜索失敗,定義臨界步長(zhǎng)因子：

s.t.h(xk+dk)+Dh(xk)d=0,

(7)

下面給出帶二階校正步的濾子算法：

步驟2:計(jì)算搜索方向.解子問(wèn)題(2)來(lái)得到搜索方向dk.若(2)不相容,轉(zhuǎn)步驟8.

步驟3: 檢查中止條件.若‖dk‖≤ε,則停止.

步驟4: 線性搜索過(guò)程:

4.1: 令l:=0,αk,0:=1.

4.3: 若(v(xk(αk,l)),f(xk(αk,l))∈Fk, 則轉(zhuǎn)步驟4.5.

4.4: 若式(3)成立,但式(4)不成立,則轉(zhuǎn)步驟4.5;若(3)和(4)都成立,則轉(zhuǎn)步驟5.

若式(3)不成立,式(5)中的條件至少有一個(gè)成立,則轉(zhuǎn)步驟5;否則,若式(3)不成立,式(5)中的條件都不成立, 則轉(zhuǎn)步驟4.5.

4.8:αk,l+1∈[T1αk,l,T2αk,l],l:=l+1,轉(zhuǎn)步驟4.2.

步驟5: 更新迭代步:αk=αk,l,xk+1=xk+αkdk.

步驟6: 用式(6)更新濾子集合.

步驟7:k:=k+1.計(jì)算fk,hk,gk,Gk,Bk,vk,轉(zhuǎn)步驟2.

步驟8: 可行性恢復(fù)階段.計(jì)算新的迭代點(diǎn)xk+1，使得約束違反度v(x)滿足v(xk+1)≤(1-γv)·v(xk).

更新濾子集合：

Fk+1:=Fk∪{(v,f)∈R2:v≥(1-γv)v(xk),f≥f(xk)-γfv(xk)},k:=k+1,轉(zhuǎn)步驟2.

2 假設(shè)條件與算法收斂性結(jié)果

假設(shè)A:

A1: 函數(shù)f(x),h(x),G(x)是二次連續(xù)可微的.

A2: 存在緊凸集合Ω，使得對(duì)任意k,xk,xk+dk∈Ω.

A3: 存在兩個(gè)常數(shù)0

A4: MFCQ約束規(guī)格[12]在任意可行點(diǎn)處成立.

定義下面的濾子更新指標(biāo)集合A={k|Fk+1≠Fk},|A|表示其中元素個(gè)數(shù),在假設(shè)條件A下,有下面的全局收斂性結(jié)果, 該定理證明的思想方法和文獻(xiàn)[15]中的相似.

定理1:在假設(shè)條件A下,考慮由算法生成的序列{xk},則下面情況之一成立:

(1) |A|<+∞, {xk}的每個(gè)聚點(diǎn)都是問(wèn)題(1)的KKT點(diǎn).

(2) |A|=+∞, {xk}k∈A的每個(gè)聚點(diǎn)都是都是問(wèn)題(1)的KKT點(diǎn).

局部收斂性的分析除假設(shè)A外,還需要下面的假設(shè)B.其中B2,B3,B4詳細(xì)解釋參考文獻(xiàn)[19], B5類似于文獻(xiàn)[20]中假設(shè).

假設(shè)B:

B1:xk→x*,其中x*是問(wèn)題(1)的KKT點(diǎn).

B2: 非退化條件在(x*,λ*,Y*)成立,其中(λ*,Y*)是相關(guān)的拉格郎日乘子.

B3: 嚴(yán)格互補(bǔ)條件成立,即

rank(G(x*))=r

rank(Y*)=m-r

B4: 帶sigma項(xiàng)的二階充分條件成立.

B5:‖(Wk-Bk)dk‖=o(‖dk‖),

下面給出局部收斂性的主要結(jié)果:

首先, 給出二階校正步的一些性質(zhì).這些性質(zhì)的證明已經(jīng)在文獻(xiàn)[19]中給出.

引理2:在假設(shè)條件A和B下,有下面的結(jié)果成立：

利用上面的結(jié)果和子問(wèn)題(2)的最優(yōu)性條件,可以推出下述結(jié)論.

引理3:在假設(shè)條件A和B下, 對(duì)充分大的k,若xk+1=xk+dk,則：

為了證明濾子算法的局部收斂性, 需要構(gòu)造下面的一個(gè)罰函數(shù)作為證明工具, 即：

Pσ(x)=f(x)+σv(x)

和這個(gè)罰函數(shù)在Pσ(xk+d)的二次近似模型：

引入的罰函數(shù)僅僅作為證明的工具,而在算法中并沒(méi)用到.同樣利用引理2,和泰勒展開(kāi)式, 可以得到下述結(jié)果.

引理4:在假設(shè)條件A,B下,若max{tr(Yk+1),‖λk+1‖∞}≤σ<+∞,則k充分大時(shí)，

(8)

根據(jù)子問(wèn)題(2)的KKT條件有：

gk+Bkdk+Dh(xk)Tλk+1+DG(xk)TYk+1=0

(9)

把式(9)兩邊和dk作內(nèi)積, 結(jié)合dk的定義可得

(10)

注意到Y(jié)k+10,將下面的不等式：

(11)

這樣從式(8,11),結(jié)合假設(shè)A4,就得到：

(12)

因此第一個(gè)結(jié)果得證.

從文獻(xiàn)[19] 中 引理 4.9 中的(4.41) 可知, 對(duì)充分大的k有

(13)

利用式(11,13)、引理2,以及假設(shè)A4可知, 對(duì)充分大的k有

(14)

(15)

因此, 從式(8，12，14)可得：

再結(jié)合式(15)，可知：

利用引理2,3,4的結(jié)果,結(jié)合文獻(xiàn)[20]中引理4.5,4.6和定理4.7的證明,可以得到:

引理5:在假設(shè)條件A和B下,對(duì)充分大的k,總有：

引理5意味著對(duì)充分大的k,全步長(zhǎng)(步長(zhǎng)因子為1)或者全步長(zhǎng)帶二階校正步總能被接受.下面引理的證明和文獻(xiàn)[17]中定理3.3的證明類似.

引理6:在假設(shè)條件A和B下, 對(duì)充分大的k,有下面的關(guān)系成立:

綜合引理5和6的結(jié)果,可以推出下面的局部收斂結(jié)論:

定理8:子問(wèn)題(2)的乘子可以用下式得到:

3 數(shù)值試驗(yàn)

為了測(cè)試算法的效果,在Intel(R) Core(TM) i5-3450, 16GB RAM 配置的電腦上進(jìn)行數(shù)值試驗(yàn). 用Matlab 2012a編寫(xiě)程序,子問(wèn)題的求解采用

SDPT3軟件包.Bk按照下面的公式進(jìn)行校正：

qk=xL(xk+1,λk+1,Yk+1)-xL(xk,λk+1,Yk+1)

算法中使用下面的參數(shù), 精度ε要求和文獻(xiàn)[19]中的一樣.

B0=I,γv=0.05,β1=0.9,β2=0.75,ν=2.1,δ=0.01,sf=5,sv=4.5,γα=0.8,ηf=0.001,T1=T2=0.6,ε=1.0e-4

測(cè)試問(wèn)題1: 控制力反饋問(wèn)題[13，22]

式中，QF=CTFTFC+I,AF=A+BFC.為方便和文獻(xiàn)[20]中方法比較,初始點(diǎn)選擇L=I,F=-0.1E,其中I為單位矩陣, 而E為元素都為1的矩陣.計(jì)算結(jié)果見(jiàn)表1.

表1 SOFP問(wèn)題測(cè)試結(jié)果Table 1 Results of SOFP

與文獻(xiàn)[20]中無(wú)懲罰型方法相比，該算法總體效果較好，但考慮到濾子算法需要儲(chǔ)存濾子集合的信息, 因此總體說(shuō)來(lái)兩種方法各有優(yōu)劣.

通過(guò)改變初始點(diǎn)進(jìn)行數(shù)值試驗(yàn), 對(duì)AC4模型取初始點(diǎn)為L(zhǎng)=I,F=O時(shí),可以看出二階校正步的效果.提高ε精度要求,分別對(duì)調(diào)用二階校正步和不用二階校正步做數(shù)值試驗(yàn),明顯前者計(jì)算效果較好,結(jié)果如表2.

表2 二階校正結(jié)果Table 2 Results of SOC

測(cè)試問(wèn)題2:近端串?dāng)_問(wèn)題[21]

式中：ri,ai為區(qū)間[0,1]上隨機(jī)數(shù)；參數(shù)P取1；初始變量Xii=ti=1,i=1,2,…，n.該測(cè)試問(wèn)題說(shuō)明定理8還原乘子后所得BFGS公式的效果. 該結(jié)果比較如圖1.(a)使用了BFGS校正公式來(lái)更新，(b)用單位矩陣計(jì)算.明顯看出使用了BFGS校正公式來(lái)更新后, 用很少的迭代次數(shù), 函數(shù)曲線即趨向于穩(wěn)定，而用單位矩陣則需要較多的迭代次數(shù), 明顯前者計(jì)算效果好.

圖1 使用BFGS公式的效果Fig.1 Results of using BFGS