基于梯度投影的广义滤子填充函数方法

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

(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的维数较高,程序的总运行时间相对长了很多.