大步長(zhǎng)路徑跟蹤內(nèi)點(diǎn)新算法

周廣付， 姚奕榮， 王筱莉

(上海大學(xué)理學(xué)院，上海200444)

在過(guò)去的20年里，內(nèi)點(diǎn)法的研究一直是非線性規(guī)劃及最優(yōu)化領(lǐng)域最引人注目的熱點(diǎn)［1-3］.內(nèi)點(diǎn)法已廣泛應(yīng)用于經(jīng)濟(jì)金融、工程控制、技術(shù)物理、物流配送、計(jì)算機(jī)科學(xué)及生物工程等領(lǐng)域，但是在實(shí)際的理論研究、測(cè)試及收斂性證明中也遇到了一些問(wèn)題：①對(duì)于初始點(diǎn)的選取，內(nèi)點(diǎn)法的核心思想是從問(wèn)題的可行域中的某一點(diǎn)出發(fā)，沿著中心路徑進(jìn)行搜索，最后到達(dá)問(wèn)題的最優(yōu)解.但是，對(duì)于一些具有多個(gè)約束的大規(guī)模問(wèn)題，一個(gè)初始可行點(diǎn)的選取是很困難的.在線性規(guī)劃問(wèn)題中，可采用一些非可行內(nèi)點(diǎn)法來(lái)克服這一困難［4-9］，但這些方法對(duì)于一般的非線性規(guī)劃問(wèn)題不能得出令人滿意的結(jié)果.在錐優(yōu)化模型中，文獻(xiàn)［10］提出了把自對(duì)偶嵌入技術(shù)推廣到隨機(jī)動(dòng)態(tài)規(guī)劃及更一般的錐優(yōu)化問(wèn)題中，從而使初始化難題得到解決.本工作通過(guò)引入一個(gè)輔助變量并在迭代過(guò)程中使該輔助變量的值逐步趨近于零，較好地解決了路徑跟蹤內(nèi)點(diǎn)算法初始點(diǎn)選取的難題.②對(duì)于收斂性的證明，路徑跟蹤內(nèi)點(diǎn)算法作為原對(duì)偶內(nèi)點(diǎn)算法的延伸和擴(kuò)展，憑借其良好的全局收斂性和較低的時(shí)間復(fù)雜度，一直受到人們的廣泛關(guān)注.但對(duì)迭代方向正交性的要求，已成為路徑跟蹤內(nèi)點(diǎn)算法運(yùn)用于求解非線性規(guī)劃問(wèn)題的一道障礙.通過(guò)對(duì)算法的研究和測(cè)試，本工作提出了一個(gè)新的關(guān)系不等式，并證明了該算法的全局收斂性.因此，本工作使路徑跟蹤內(nèi)點(diǎn)法的應(yīng)用得以拓展.

1 KKT條件

首先，考慮如下約束優(yōu)化問(wèn)題：

式中，f，ɡi：Rn→R(i=1，2，…，m)是二次連續(xù)可微的.

引入變量x0，則問(wèn)題(P)轉(zhuǎn)化為

定義該問(wèn)題的拉格朗日函數(shù)為

式中，

對(duì)KKT條件引入松弛變量s=(s1，s2，…，sm)T，可得到式(4)的等價(jià)形式為

式中，s，z≥0，S=diag(s1，s2，…，sm).

在實(shí)際求解過(guò)程中，通常采用Newton法來(lái)求解上述KKT條件等式，那么Newton方程的解即為Newton迭代方向Δw=(Δx-，Δs，Δy，Δz)T，有

2 路徑跟蹤內(nèi)點(diǎn)算法

中心路徑在內(nèi)點(diǎn)算法中扮演著非常重要的角色，它是由一類嚴(yán)格可行點(diǎn)組成的弧形軌跡.若有參數(shù)τ＞0，則中心路徑中的點(diǎn)可通過(guò)下列方程式進(jìn)行求解：

本工作稱條件(9)為修正的KKT條件.條件(9)與條件(7)之間唯一的差別在于，條件(9)在最后一行中引入了參數(shù)τ，即把條件(7)第三行中的互補(bǔ)條件換成了要求sizi乘積值對(duì)任意下標(biāo)i均相等的條件.由條件(9)可定義中心路徑C={(，sτ，yτ，zτ)| τ＞0}.當(dāng)τ趨近于0時(shí)，條件(9)等價(jià)于條件(7);當(dāng)τ下降到0時(shí)，C收斂于某一點(diǎn)，那么該點(diǎn)必為原問(wèn)題的最優(yōu)點(diǎn).因此，沿著中心路徑搜索的方法提供了一條找到最優(yōu)解的途徑，沿著該途徑不僅可以保證s和z各個(gè)分量都嚴(yán)格大于0，還可以使所有sizi的乘積值幾乎以相同的速率下降到0.

當(dāng)sizi的各個(gè)分量均等于σμ時(shí)，Newton迭代方向(Δ，Δs，Δy，Δz)指向點(diǎn)(，sσμ，yσμ，zσμ)∈C;反之，由式(8)求得的迭代方向則直接指向滿足KKT條件(7)的點(diǎn).

下面，給出式(10)可解的充分條件.

引理1 如果矩陣H+ΔgT(x)S-1ZΔg(x)是正定的，則矩陣N(w)非奇異.

證明 考慮方程式

式中，(vx0，vx，vs，vy，vz)∈R1×Rn×Rm×R1×Rm，則有

由假設(shè)知，vx=0，因此，vx0=0，vs=0，vy=0，vz=0.由此，可求得(Δ，Δs，Δy，Δz)T.證畢.

路徑跟蹤算法要求每步產(chǎn)生的迭代點(diǎn)均落在一個(gè)包含中心路徑C的區(qū)域內(nèi)，且沿著這條中心路徑可搜索至原問(wèn)題的最優(yōu)解.為了避免迭代點(diǎn)過(guò)于接近非負(fù)區(qū)域的邊界，本工作要求每次迭代中的搜索方向必須使下一步產(chǎn)生的點(diǎn)列更加接近最優(yōu)點(diǎn).

最優(yōu)性算法的一個(gè)要點(diǎn)是如何在搜索空間中衡量滿足條件的點(diǎn)，因此，在路徑跟蹤算法中，對(duì)偶參數(shù)μ擔(dān)當(dāng)了這個(gè)重要的角色.當(dāng)k→∞時(shí)，μk下降至0，因此，迭代點(diǎn)(k，sk，yk，zk)逐步滿足KKT條件(7).

現(xiàn)定義單邊∞范數(shù)區(qū)域N-∞(γ)如下：

式中，γ∈(0，1］，F(xiàn)0={(，s，y，z)|s＞0}.若存在一點(diǎn)落在N-∞(γ)中，則對(duì)每個(gè)分量sizi必定大于μ與γ的乘積.這個(gè)約束是非常松的，當(dāng)γ→0時(shí)，N-∞(γ)就包含可行域F={(，s，y，z)|s≥0}中的大部分點(diǎn).本工作取γ=10-3.

下面介紹的大步長(zhǎng)路徑跟蹤算法，由于采用了一個(gè)較為廣泛的區(qū)域N-∞(γ)(γ→0)，因此，該算法具有良好的實(shí)用性.通過(guò)式(10)可求得搜索方向，步長(zhǎng)αk則取保證迭代點(diǎn)落于N-∞(γ)內(nèi)的最大值.

定義

大步長(zhǎng)路徑跟蹤內(nèi)點(diǎn)算法步驟如下：

(1)給定 ε＞0，γ，σ∈(0，1)，取 x0及＞max(ɡi(x0))，s0=0，y0和z0;

(2)如果‖r(wk)‖≤ε，停止;

(4)取αk作為α在區(qū)間［0，1］中的最大值，且有((αk)，sk(αk)，yk(αk)，zk(αk))∈N-∞(γ);

(6)令k∶=k+1，轉(zhuǎn)步驟(1).

3 收斂性分析

下面部分分析上述算法的收斂性.首先，介紹引理2［4］，并用它來(lái)證明引理3.引理3給出了ΔiΔzi(i=1，2，…，n)向量乘積的邊界.定理1給出了αk的一個(gè)下確界以及迭代過(guò)程中測(cè)度μ下降量的估計(jì)值，由此即可得算法的全局收斂性.為了證明算法的收斂性，給出如下假設(shè)：

(1)矩陣H+ΔɡT(x)S-1ZΔɡ(x)正定;

(2)不等式0≤ΔsTΔz≤(1-σ)sTz成立.

引理2 假設(shè)u，v是任意2個(gè)n維向量，且有uTv≥0，那么

式中，U=diag(u1，u2，…，un)，V=diag(v1，v2，…，vn).

證明 由假設(shè)(2)，易得

把式(10)的最后一行兩邊分別乘以(SZ)-1/2，并令D=S1/2Z-1/2，可得

因?yàn)?/p>

令u=D-1Δs，v=DΔz，由引理2，得

利用展開平方歐氏范數(shù)和關(guān)系式 sTz=nμ，eTe=n，有

證畢.

定理1 由假設(shè)(1)和(2)，給定參數(shù)γ，σ∈(0，1)，則存在一個(gè)常數(shù)δ∈(0，1)，使得

對(duì)所有的k≥0成立.

由于αk是［0，1］區(qū)間內(nèi)α的最大值，那么αk有下確界

由引理3可知，對(duì)任意的i=1，2，…，n，有

累加方程SkΔzk+ZkΔsk=-SkZke+σμke(由式(10)的最后一行求得)的n個(gè)分量，利用式(17)和μk，μk(α)的定義，即得

整理上式，可得

接下來(lái)，估計(jì)在第k步迭代中測(cè)度μ的下降量.由式(10)，(16)，(17)的最后一行以及假設(shè)，可得

由此可知，α(1-α)是一個(gè)關(guān)于α的二次凹函數(shù)，并且，在任意給定的區(qū)間范圍內(nèi)，函數(shù)的最小值可在區(qū)間的端點(diǎn)處取到.所以，在最后一式中引入替代估計(jì)參數(shù)

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

本工作測(cè)試了3個(gè)數(shù)值算例，結(jié)果表明，所提出的大步長(zhǎng)路徑跟蹤內(nèi)點(diǎn)算法是可行的.

例1

例2

表1 算法迭代9步Table 1 Computations 9 iterations

表2 算法迭代6步Table 2 Computations 6 iterations

［1］ BOY，QINGX，F(xiàn)ENGG C.On the complexity of a combined homotopy interior method for convex programming［J］.Journal of Computational and Applied Mathematics，2007，200：32-46.

［2］ GEORG S.Path-following and augmented lagrangian methods for contact problems in linear elasticity［J］.Journal of Computational and Applied Mathematics，2007，203：533-547.

［3］ RENATOD C M，TAKASHIT.A strong bound on the integral of the central path curvature and its relationship with the iteration-complexity of primal-dual path-following LP algorithms［J］.Mathematical Programming，2008，115：105-149.

［4］ JORGEN，WRIGHTS J.Numerical optimization［M］.Boston：Springer，1999：60-112.

［5］ JIMB，SONGX.A non-interior predictor-corrector path following algorithm for the monotone linear complementarity problem ［J］. Mathematical Programming，2000，87：113-130.

［6］ BURKEJ，XUS.Complexity of a noninterior pathfollowing method for the linear complementarity problem［J］.Journal of Optimization Theory and Applications，2002，112：53-76.

［7］ FREUNDR M.A potential-function reduction algorithm for solving a linear program directly from an infeasible warm start［J］.Mathematical Programming，1991，52：441-466.

［8］ FORSGRENA.On warm starts for interior methods［M］.Boston：Springer，2006：51-66.

［9］ CORNELISR，TAMáST，JEAN-PHILIIPE V.Interior point methods for linear optimization［M］.Boston：Springer，2006：55-76.

［10］ ZHANGS Z.A new self-dual embedding method for convex programming ［J］. Journal of Global Optimization，2004，29：479-496.