大步长路径跟踪内点新算法

周广付， 姚奕荣， 王筱莉

(上海大学理学院，上海200444)

在过去的20年里，内点法的研究一直是非线性规划及最优化领域最引人注目的热点［1-3］.内点法已广泛应用于经济金融、工程控制、技术物理、物流配送、计算机科学及生物工程等领域，但是在实际的理论研究、测试及收敛性证明中也遇到了一些问题：①对于初始点的选取，内点法的核心思想是从问题的可行域中的某一点出发，沿着中心路径进行搜索，最后到达问题的最优解.但是，对于一些具有多个约束的大规模问题，一个初始可行点的选取是很困难的.在线性规划问题中，可采用一些非可行内点法来克服这一困难［4-9］，但这些方法对于一般的非线性规划问题不能得出令人满意的结果.在锥优化模型中，文献［10］提出了把自对偶嵌入技术推广到随机动态规划及更一般的锥优化问题中，从而使初始化难题得到解决.本工作通过引入一个辅助变量并在迭代过程中使该辅助变量的值逐步趋近于零，较好地解决了路径跟踪内点算法初始点选取的难题.②对于收敛性的证明，路径跟踪内点算法作为原对偶内点算法的延伸和扩展，凭借其良好的全局收敛性和较低的时间复杂度，一直受到人们的广泛关注.但对迭代方向正交性的要求，已成为路径跟踪内点算法运用于求解非线性规划问题的一道障碍.通过对算法的研究和测试，本工作提出了一个新的关系不等式，并证明了该算法的全局收敛性.因此，本工作使路径跟踪内点法的应用得以拓展.

1 KKT条件

首先，考虑如下约束优化问题：

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

引入变量x0，则问题(P)转化为

定义该问题的拉格朗日函数为

式中，

对KKT条件引入松弛变量s=(s1，s2，…，sm)T，可得到式(4)的等价形式为

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

在实际求解过程中，通常采用Newton法来求解上述KKT条件等式，那么Newton方程的解即为Newton迭代方向Δw=(Δx-，Δs，Δy，Δz)T，有

2 路径跟踪内点算法

中心路径在内点算法中扮演着非常重要的角色，它是由一类严格可行点组成的弧形轨迹.若有参数τ＞0，则中心路径中的点可通过下列方程式进行求解：

本工作称条件(9)为修正的KKT条件.条件(9)与条件(7)之间唯一的差别在于，条件(9)在最后一行中引入了参数τ，即把条件(7)第三行中的互补条件换成了要求sizi乘积值对任意下标i均相等的条件.由条件(9)可定义中心路径C={(，sτ，yτ，zτ)| τ＞0}.当τ趋近于0时，条件(9)等价于条件(7);当τ下降到0时，C收敛于某一点，那么该点必为原问题的最优点.因此，沿着中心路径搜索的方法提供了一条找到最优解的途径，沿着该途径不仅可以保证s和z各个分量都严格大于0，还可以使所有sizi的乘积值几乎以相同的速率下降到0.

当sizi的各个分量均等于σμ时，Newton迭代方向(Δ，Δs，Δy，Δz)指向点(，sσμ，yσμ，zσμ)∈C;反之，由式(8)求得的迭代方向则直接指向满足KKT条件(7)的点.

下面，给出式(10)可解的充分条件.

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

证明 考虑方程式

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

由假设知，vx=0，因此，vx0=0，vs=0，vy=0，vz=0.由此，可求得(Δ，Δs，Δy，Δz)T.证毕.

路径跟踪算法要求每步产生的迭代点均落在一个包含中心路径C的区域内，且沿着这条中心路径可搜索至原问题的最优解.为了避免迭代点过于接近非负区域的边界，本工作要求每次迭代中的搜索方向必须使下一步产生的点列更加接近最优点.

最优性算法的一个要点是如何在搜索空间中衡量满足条件的点，因此，在路径跟踪算法中，对偶参数μ担当了这个重要的角色.当k→∞时，μk下降至0，因此，迭代点(k，sk，yk，zk)逐步满足KKT条件(7).

现定义单边∞范数区域N-∞(γ)如下：

式中，γ∈(0，1］，F0={(，s，y，z)|s＞0}.若存在一点落在N-∞(γ)中，则对每个分量sizi必定大于μ与γ的乘积.这个约束是非常松的，当γ→0时，N-∞(γ)就包含可行域F={(，s，y，z)|s≥0}中的大部分点.本工作取γ=10-3.

下面介绍的大步长路径跟踪算法，由于采用了一个较为广泛的区域N-∞(γ)(γ→0)，因此，该算法具有良好的实用性.通过式(10)可求得搜索方向，步长αk则取保证迭代点落于N-∞(γ)内的最大值.

定义

大步长路径跟踪内点算法步骤如下：

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

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

(4)取αk作为α在区间［0，1］中的最大值，且有((αk)，sk(αk)，yk(αk)，zk(αk))∈N-∞(γ);

(6)令k∶=k+1，转步骤(1).

3 收敛性分析

下面部分分析上述算法的收敛性.首先，介绍引理2［4］，并用它来证明引理3.引理3给出了ΔiΔzi(i=1，2，…，n)向量乘积的边界.定理1给出了αk的一个下确界以及迭代过程中测度μ下降量的估计值，由此即可得算法的全局收敛性.为了证明算法的收敛性，给出如下假设：

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

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

引理2 假设u，v是任意2个n维向量，且有uTv≥0，那么

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

证明 由假设(2)，易得

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

因为

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

利用展开平方欧氏范数和关系式 sTz=nμ，eTe=n，有

证毕.

定理1 由假设(1)和(2)，给定参数γ，σ∈(0，1)，则存在一个常数δ∈(0，1)，使得

对所有的k≥0成立.

由于αk是［0，1］区间内α的最大值，那么αk有下确界

由引理3可知，对任意的i=1，2，…，n，有

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

整理上式，可得

接下来，估计在第k步迭代中测度μ的下降量.由式(10)，(16)，(17)的最后一行以及假设，可得

由此可知，α(1-α)是一个关于α的二次凹函数，并且，在任意给定的区间范围内，函数的最小值可在区间的端点处取到.所以，在最后一式中引入替代估计参数

4 数值试验

本工作测试了3个数值算例，结果表明，所提出的大步长路径跟踪内点算法是可行的.

例1

例2

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

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

［1］ BOY，QINGX，FENGG 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.