本文讨论如下带线性约束的全局优化问题,

$ \begin{equation} \begin{aligned} &\min \ f(x), \\ & \ {\rm s.t.} \ Ax\leq b, \end{aligned} \end{equation} $

(1.1)

其中$ f(x): R^n \rightarrow R $为非凸函数, $ A $为$ m\times n $阶矩阵, $ b = (b_1, b_2, b_3, \cdots, b_m)^{\rm T}\in R^m $.对约束问题(1.1), 称所有满足约束的点为可行点, 由所有可行点组成的集合$ X = \{x\in R^n:Ax\leq b\} $为可行域.一般说来, 问题(1.1)有多个极小值, 用传统方法求解只能得到局部最优, 而且对于全局最优解的判定条件至今都缺少实用的研究成果, 这都使得问题求解仍然面临着许多困难.

填充函数法是求解多极值最优化问题的常用算法之一, 其主要思想是:在求得问题的一个局部极小点后, 构造填充函数.通过极小化该填充函数, 寻找比当前局部极小点更优的点, 进而得到更优的极小值^[1-5].与其他方法相比, 该算法思想简单, 而且仅用到经典算法, 所以易于实现且效率较高, 同时也可以推广到其他非线性问题以及离散数学规划的求解中^[6].

滤子技术最早由Fletcher和Leyffer提出, 他们详细讨论了滤子作为代替罚函数的工具在局部优化算法中的一些应用^{[7, 8]}.之后滤子技术在局部优化问题的求解中被认为是一种更有效的方法, 因其良好的数值效果, 许多学者继续进行了一系列的相关研究^[9-11].

梯度投影算法自从被Rosen^[12]提出后就引起了广泛的注意和系统的研究^{[13, 14]}, 由于该方法简单、实际应用的数值效果好, 在一些更有效的近代算法中也继续沿用了它的基本思想^{[15, 16]}.本文为了优化填充函数算法, 将滤子技术和改进的梯度投影方法融入到全局优化算法中, 提出了基于梯度投影的广义滤子填充函数算法求解问题(1.1).

为方便起见, 下面做一些记号说明.

$ J = \{1, 2, \cdots, m\} $为指标集, $ A $的第$ j $行为$ a_j^{\rm T} $, 记$ c_j(x) = a_j ^{\rm T}x-b_j, j\in J $.

记约束违反度函数$ h(x) = \max (0, c_j(x), j\in J) $, $ A(x) = (a_j, j\in J_0(x)) $, 其中$ J_0(x) = \{j|c_j(x)\geq 0, j\in J\} $, 并且记$ \overline{J_0(x)} = J\setminus J_0(x) $.

$ L(P) $表示问题(1.1)的局部最优解集合, $ G(P) $表示问题(1.1)的全局最优解集合.

$ x^* $是$ f(x) $的一个稳定点, $ S_1 (x^* ) = \{x\; |\; f(x)\geq f(x^* ), x\in X\setminus\{x^* \} \} $称为高水平集, $ S_2(x^*) = \{x\; |\; f(x)<f(x^* ), x\in X\} $称为低水平集.

$ {\rm int} X $表示可行域$ X $的内点集合, $ \partial X $表示可行域$ X $的边界点集合.

考虑问题(1.1), 我们首先提出以下假设:

[A1] $ f(x) $只有有限多个极小值, 即存在直径充分大的闭箱$ \Omega $能包含所有极小值.

[A2] $ \nabla f(x) $在$ \Omega $上连续.

[A3] $ \forall x\in R^n $, $ \{a_j, j\in J_0(x)\} $为线性无关向量组.

由A1和A2可知, 原问题的全局解必在有界闭区域$ X\bigcap\Omega $上, 因此可以认为$ X $是有界闭集, 且算法中的$ x $总是取自$ \Omega $.

定义2.1   设$ x^* $是问题$ \mathrm{(1.1)} $当前的局部极小值, 称$ T(x, x^*) $是$ f(x) $在$ x^* $处的广义填充函数, 如果$ T(x, x^* ) $满足

(ⅰ) $ T(x, x^*) $在高水平集$ S_1(x^*) $上没有稳定点;

(ⅱ)如果$ x^\ast $不是问题$ \mathrm{(1.1)} $的一个全局极小值点, 即$ x^* \notin G(P) $, 那么存在点$ x_1^* \in S_2(x^*) $, 使得$ x_1^* $是$ T(x, x^*) $的极小值点.

下面给出求解问题(1.1)的广义填充函数.

设$ x^* $是问题(1.1)的一个局部极小点, 定义单参数广义填充函数

$ \begin{equation} T(x, x^*, r) = \dfrac{1-{\rm e}^{-\frac{1}{r^2}[f(x)-f(x^*)+r]}}{1+\|x-x^*\|}, \end{equation} $

(2.1)

其中

$ \begin{equation} 0<r<\min\{|f(x_1^*)-f(x_2^*)|:f(x_1^*)\neq f(x_2^*);x_1^*, x_2^*\in L(P)\}. \end{equation} $

(2.2)

根据[A2], 显然$ T(x, x^*, r) $在可行域上是连续可微的.下面来讨论$ T(x, x^*, r) $的性质.

定理2.1   对充分小的$ r>0 $, 有如下结论成立

(ⅰ)函数$ T(x, x^*, r) $在集合$ S_1(x^*)\setminus \partial X $上没有稳定点;

(ⅱ)当$ x-x^* $与$ \{a_j, j\in J_0 (x)\} $线性无关时, 函数$ T(x, x^*, r) $在集合$ S_1 (x^* )\bigcap \partial X $上没有稳定点.

证 (ⅰ) $ \forall x\in S_1(x^*)\setminus \partial X $, 有

$ \begin{eqnarray*} &&\dfrac{(x-x^*)^{\rm T}}{\|x-x^*\|}\nabla T(x, x^*, r)\\ & = &\dfrac{(x-x^*)^{\rm T}}{\|x-x^*\|} \frac{ -{\rm e}^{-\frac{1}{r^2}[f(x)-f(x^*)+r]}\big[-\frac{1}{r^2}\nabla f(x)\big] (1+\|x-x^*\|)-\big\{1-{\rm e}^{-\frac{1}{r^2}[f(x)-f(x^*)+r]}\big\}\frac{x-x^*}{\|x-x^*\|}}{(1+\|x-x^*\|)^2}\\ & = & \frac{{\rm e}^{-\frac{1}{r^2}[f(x)-f(x^*)+r]}}{r^2(1+\|x-x^*\|)}\; \bigg\{\frac{(x-x^*)^{\rm T}}{\|x-x^*\|}\nabla f(x) -r^2{\rm e}^{ \frac{1}{r^2}[f(x)-f(x^*)+r]}T(x, x^*, r)\bigg\}. \end{eqnarray*} $

根据A2, $ \big|\frac{(x-x^*)^{\rm T}}{\|x-x^*\|}\nabla f(x)\big|\leq \|\nabla f(x)\| $在可行域内是有界的.而由$ x\in S_1(x^*)\setminus\partial X $可知$ f(x)\geq f(x^\ast) $, 则$ T(x, x^*, r)>0 $.所以当$ r>0 $充分小时, 有

$ \begin{equation} \frac{(x-x^*)^{\rm T}}{\|x-x^*\|}\nabla T(x, x^*, r)<0, \end{equation} $

(2.3)

即函数$ T(x, x^*, r) $在集合$ S_1(x^*)\setminus \partial X $上没有稳定点.

(ⅱ)假设$ \exists x_L\in S_1 (x^* )\bigcap \partial X $是$ T(x, x^*, r) $的稳定点, 则有$ \Lambda(x_L ) = (\lambda_1, \lambda_2, \cdots, \lambda_m)^{\rm T} $, 使得

$ \begin{equation} \nabla T(x_L, x^*, r) = \sum\limits_{j = 1}^m \lambda_j a_j \end{equation} $

(2.4)

成立, 其中$ \lambda_j\geq 0, \lambda_j c_j (x_L ) = 0, j\in J $.而由(ⅰ)的证明可知

$ \begin{equation} \nabla T(x_L, x^*, r) = \dfrac{{\rm e}^{-\frac{1}{r^2} [f(x_L )-f(x^* )+r] }}{r^2 (1+\|x_L-x^* \|)}\nabla f(x_L )-\dfrac{T(x_L, x^*, r)}{1+\|x_L-x^*\|}\dfrac{x_L-x^*}{\|x_L-x^*\|} \end{equation} $

(2.5)

且当$ r>0 $充分小时, 显然上式右端第一项趋于$ 0 $.则由(2.4), (2.5)式知

$ x_L-x^* = -\sum\limits_{j = 1}^m \dfrac{\lambda_j \|x_L-x^*\|(1+\|x_L-x^*\|)}{T(x_L, x^*, r)}a_j, $

这与已知条件矛盾, 故函数$ T(x, x^*, r) $在集合$ S_1 (x^* )\bigcap X $上没有稳定点.

由(2.3)式, 显然有

推论2.1   $ x-x^* $是函数$ T(x, x^*, r) $在$ S_1(x^*)\setminus\partial X $处的一个严格下降方向.

定理2.2  如果$ x^* \in L(P) $, 但是$ x^* \notin G(P) $, 那么$ T(x, x^*, r) $在$ S_2(x^*) $上至少有一个极小值点.

证   由于$ G(P) $非空和(2.2)式的定义, 则必有$ x_G \in G(P) $, 使得$ T(x_G, x^*, r)<0 $.而$ T(x, x^*, r) $在$ X $上连续, 所以必存在函数的最小值点$ \bar{x} \in X $, 使得

$ T(\bar{x}, x^*, r)\leq T(x_G, x^*, r)<0. $

也就是说$ 1-{\rm e}^{-\frac{1}{r^2}[f(x)-f(x^*)+r]}<0 $, 即$ f(\bar x)< f(x^*) $, $ \bar{x} \in S_2(x^*) $.

推论2.2   如果$ x\in G(P) $, 则对充分小的$ r>0 $, 函数$ T(x, x^*, r) $在可行域内部$ {\rm int} X $没有稳定点.

由定理2.1和定理2.2可知, $ T(x, x^*, r) $是一个广义填充函数.

滤子最初定义为由两个相互竞争的函数$ \phi(x) $和$ \psi(x) $组成的数对集合, 记为$ (\phi(x), \psi(x)) $.为了之后讨论方便, 本小节将给定的一阶连续可微函数$ Q(x) $作为目标函数, 即考虑如下问题

$ \begin{equation} \begin{aligned} &\min \ Q(x), \\ & \ {\rm s.t.} \ Ax\leq b. \end{aligned} \end{equation} $

(3.1)

取$ \phi(x) $为$ Q(x) $, $ \psi(x) $为$ h(x) $, 下面给出"支配"的相关概念.

定义3.1   点$ x^k $支配点$ x^l $当且仅当$ Q(x^k)\leq Q(x^l) $且$ h(x^k)\leq h(x^l) $.

定义3.2   滤子是由一列互不支配的数对组成, 即若点$ x^k $和$ x^l $都在滤子中, 那么$ Q(x_k)\leq Q(x_l) $和$ h(x_k)\leq h(x_l) $两者必有之一不成立.

此外, 为了保证算法的下降性, 在本文中如果说点$ x_{k+1} $"被滤子接受"当且仅当存在滤子中的点$ x_l $, 使得下面两个不等式

$ \begin{eqnarray} &&Q(x_{k+1})< Q(x_l)-\beta h(x_l), \end{eqnarray} $

(3.2)

$ \begin{eqnarray} &&h(x_{k+1})< (1-\eta)h(x_l) \end{eqnarray} $

(3.3)

之一成立, 其中$ \beta, \eta \in (0, 1) $.

从以上概念可以看出, 滤子能够作为一个评判标准来决定对当前试探步(沿下降方向寻找合适步长时假设的迭代点)是否接受.但是, 以(3.2)和(3.3)式作为滤子集的接受标准可能导致算法收敛到一个可行的非稳定点.因此, 本文在下降搜索时又另外采用了其他准则.

对当前迭代点$ x_k $, 记

$ \begin{eqnarray*} &&U(x_k) = (u_j(x_k), j\in J_0(x))^{\rm T} = -(A(x_k)^{\rm T} A(x_k))^{-1} A(x_k)^{\rm T} \nabla Q(x_k), \\ &&P(x_k) = I-A(x_k)(A(x_k)^{\rm T} A(x_k))^{-1}A(x_k)^{\rm T}, \end{eqnarray*} $

其中$ I $为$ n $阶单位矩阵, 称$ P(x_k) $为$ x_k $处的投影矩阵, 并记$ x_k $处的搜索方向为

$ \begin{equation} d_k = \begin{cases}-P(x_k)\nabla Q(x_k), &\mbox{当}\; \; x_k\in X, \\ -P(x_k)\nabla Q(x_k)+\rho_k B(x_k)^{\rm T}\omega, &\mbox{当}\; \; x_k\notin X, \end{cases} \end{equation} $

(3.4)

其中

$ \begin{eqnarray*} &&\rho_k = \dfrac{\nabla Q(x_k)^{\rm T} P(x_k)\nabla Q(x_k)+h(x_k)}{2|U(x_k)^{\rm T}\omega|+1}, \\ &&B(x_k) = (A(x_k)^{\rm T} A(x_k))^{-1}A(x_k)^{\rm T}, \; \; \omega = (\omega_j, j\in J_0(x_k))^{\rm T}, \; \; \omega_j = -1. \end{eqnarray*} $

对于当前点的试探步长$ \alpha_k $, 本文要求一个足够下降量标准(转换条件)

$ \begin{equation} m_k(\alpha_k)<0\ \ , \ \ [-m_k(\alpha_k)]^{s_1}(\alpha_k)^{1-{s_1}}>\delta_1[h(x_k)]^{s_2} \end{equation} $

(3.5)

成立.固定参数为$ \delta_1>0 $, $ s_2>1 $, $ s_1>2s_2 $, 其中$ m_k(\alpha_k): = \alpha_k\nabla Q(x_k)^{\rm T} d_k. $

若转换条件(3.5)成立, 则可以不再受滤子接受准则约束, 只要求目标函数的Armijo条件

$ \begin{equation} Q(x_k(\alpha_k))\leq Q(x_k)+\delta_2 m_k(\alpha_k) \end{equation} $

(3.6)

成立.其中$ \delta_2\in(0, \frac{1}{2}) $是一个固定的常数.不难看出, 若对某一试探步长$ \alpha_k $, 转换条件(3.5)成立, 而Armijo条件(3.6)不成立, 那么对于再小一点的步长, 转换条件也将有可能不再成立.

有时算法在搜索过程中目标函数逐步下降但迭代点却离可行域越来越远, 此时就需要进入可行性恢复阶段.下面简述下可行性恢复的具体过程.

设点$ x\in\Omega $但$ x\notin X $, 给定一个固定的参数$ \varepsilon>0 $, 不妨设有$ c_{i_j}(x)> \varepsilon $, $ j = 1, 2, \cdots, k $, $ c_{i_j}(x)\leq\varepsilon $, $ j = k+1, \cdots, m $, 即存在$ k $个约束不在可接受的范围内, 需要进行可行性恢复, 则求解问题

$ \begin{equation} \begin{aligned} & \min \ h(x), \\ & {\rm s.t.} \ \ c_{i_j}(x)\leq\varepsilon, j = k+1, \cdots, m. \end{aligned} \end{equation} $

(3.7)

对于进入可行性恢复阶段的判定, 本文采用如下准则.

如果当前迭代步$ \alpha_{k} $足够大, 转换条件(3.5)对于某个$ \alpha\leq\alpha_k $是成立的, 那么仍有可能存在某个小一点的步长能被Armijo条件(3.6)接受, 此时就不用进入可行性恢复阶段.故由(3.5)式可知, 若$ \nabla Q(x_k)^{\rm T}d_k<0 $, 只要$ \alpha_k>\dfrac{\delta_1[h(x_k)]^{s_2}}{[-\nabla Q(x_k)^{\rm T}d_k]^{s_1}} $, 就不进入可行性恢复阶段.

然而, 当转换条件(3.5)不成立时, 考虑如下两个线性近似

$ \begin{eqnarray*} &&Q(x_k+\alpha_k d_k) = Q(x_k)+\alpha_k\nabla Q(x_k)^{\rm T} d_k+o(\|\alpha_{k}d_k\|), \\ &&h(x_k+\alpha_k d_k) = \max\{0, c_j(x_k)+\alpha_k a_j^{\rm T} d_k+o(\|\alpha_{k}d_k\|), j\in\overline{J_0(x)}\}. \end{eqnarray*} $

当$ a_j^{\rm T}d_k<0 $时, 若$ \alpha_k\leq -\dfrac{\eta c_j(x_k)}{a_j^{\rm T} d_k} $, 则无法满足约束违反度函数的足够下降量.类似地, 当$ \nabla Q(x_k)^{\rm T} d_k<0 $时, 若$ \alpha_k\leq -\dfrac{\beta Q(x_k)}{\nabla Q(x_k)^{\rm T} d_k} $, 目标函数的足够下降量也不再满足.即当$ a_j^Td_k<0 $时, 对目标函数的下降存在一个最小试探步长

$ \begin{equation} \alpha_k^{\min}: = \theta \cdot \begin{cases} \min\bigg\{\dfrac{\delta_1[h(x_k)]^{s_2}}{[-\nabla Q(x_k)^{\rm T}d_k]^{s_1}}, -\dfrac{\beta Q(x_k)}{\nabla Q(x_k)^{\rm T} d_k}, -\dfrac{\eta c_j(x_k)}{a_j^{\rm T} d_k}\bigg\}, &\mbox{若}\; \; \nabla Q(x_k)^{\rm T}d_k<0, \\ -\dfrac{\eta c_j(x_k)}{a_j^{\rm T} d_k}, &\mbox{否则}. \end{cases} \end{equation} $

(3.8)

当步长$ \alpha_k<\alpha_k^{\min} $时, 算法进入可行性恢复阶段, 其中$ \theta \in (0, 1] $为一常数, 用来补充线性化时被省略的高阶项, 避免不必要的可行性恢复.

当$ a_j^{\rm T}d_k = 0 $时, 最小试探步长即为

$ \begin{equation} \alpha_k^{\min}: = \theta \cdot \min\bigg\{\dfrac{\delta_1[h(x_k)]^{s_2}}{[-\nabla Q(x_k)^{\rm T}d_k]^{s_1}}, -\dfrac{\beta Q(x_k)}{\nabla Q(x_k)^{\rm T} d_k}\bigg\}. \end{equation} $

(3.9)

本文研究的是带线性约束的全局优化问题, 同时从目标函数、填充函数和约束违反度函数三个角度考虑, 所以将一般滤子推广到三维, 构成滤子$ (f(x), T(x, x^*), h(x)) $.即

(ⅰ)点$ x_k $支配点$ x_l $当且仅当$ f(x_k)\leq f(x_l) $, $ T(x_k, x^*)\leq T(x_l, x^*) $和$ h(x_k)\leq h(x_l) $同时成立.

(ⅱ)若点$ x_k $和点$ x_l $都在滤子中, 那么在$ f(x_k)\leq f(x_l) $, $ T(x_k, x^*)\leq T(x_l, x^*) $和$ h(x_k)\leq h(x_l) $中至少有一个不等式不成立.

在本文中, 用$ F_k $表示点$ \{x_l|l\leq k\} $的集合, 也就是说$ (f(x_l), T(x_l, x^*), h(x_l)) $是当前滤子中的元素, $ F_k $为当前滤子.另外, $ |F| $表示滤子$ F $所包含的元素个数.

此外, 如果说点$ x_{k+1} $"被滤子$ F_k $接受"当且仅当存在点$ x_l\in F_k $, 使得下面三个不等式

$ \begin{eqnarray} &&f(x_{k+1})< f(x_l)-\beta_1 h(x_l), \end{eqnarray} $

(4.1)

$ \begin{eqnarray} &&T(x_{k+1}, x^*)< T(x_l, x^*)-\beta_2 h(x_l), \end{eqnarray} $

(4.2)

$ \begin{eqnarray} &&h(x_{k+1})< (1-\eta)h(x_l) \end{eqnarray} $

(4.3)

之一成立, 其中$ \beta_1, \beta_2, \eta \in(0, 1) $.

定义4.1   将数组$ (f(x), T(x, x^*), h(x)) $加入到当前滤子中, 并把被$ (f(x), T(x, x^*), h(x)) $支配的所有数组移出滤子的过程称为更新滤子.即

$ \begin{eqnarray} F_{k+1}& = &F_k\cup\{x_{k+1}\}\setminus \{x_l\in F_k |f(x_{k+1})\\ &\leq& f(x_l), T(x_{k+1}, x^*)\leq T(x_l, x^*), h(x_{k+1})\leq h(x_l)\}. \end{eqnarray} $

(4.4)

根据三维滤子的定义可知, 如果初始滤子集非空, 那么在任意点处的滤子集非空, 即$ |F|\geq 1 $.在提出广义滤子填充函数算法之前, 我们先给出试探步长$ \alpha_k $的迭代算法.

回溯线搜索算法

步骤0  令$ \alpha_k = 1 $.

步骤1  若$ \alpha_k<\alpha_k^{\min} $时, 转步骤5;否则, 计算下一迭代$ x_k(\alpha_k) = x_k+\alpha_k d_k $.

步骤2  若滤子中存在数组能支配$ x_k(\alpha_k) $, 则拒绝试探步, 转步骤4;否则, 进入步骤3.

步骤3  检验足够下降量

3.1情况I  若$ x_k\in X $时, (3.5)和(3.6)式同时成立, 并且$ x_k(\alpha_k)\in X $, 则$ \alpha_k $为所求步长; 否则, 转步骤4.

3.2情况II   若$ x_k\notin X $且(3.5)式成立时, (3.6)式成立, 则$ \alpha_k $为所求步长; 否则, 转步骤4.

3.3情况III   若$ x_k\notin X $且(3.5)式不成立时, 试探点能被滤子接受, 则$ \alpha_k $为所求步长; 否则, 转步骤4.

步骤4  选择新的试探步长$ \alpha_k = \dfrac{1}{2}\alpha_k $, 转步骤2.

步骤5  可行性恢复:求解问题(3.7), 并更新滤子集.

本算法的目的是对当前迭代点$ x_k $和方向$ d_k $, 寻找合适的下降步长.当$ x_k $不需要进行可行性恢复且下一步迭代$ x_k(\alpha_k) $不被滤子支配时, 试探步长$ \alpha_k $还需要满足足够下降量或者滤子接受准则.算法中对足够下降量的检验作了详细的情况分析, 确保当前点无论是否在可行域内一定存在可被接受的步长.

下面给出主算法.

算法

步骤0  任取$ x_0\in R^n $, 给定精度参数$ \varepsilon>0 $, 邻域半径$ \delta>0 $, $ r_0, r, G, N $为变量$ x $的维数.令$ k: = 0 $, 若当前点没有填充函数值, 则令其分量$ T = Z $, $ Z $充分大.初始化滤子集$ F_0 = {(f(x_0), Z, h(x_0))} $.

步骤1  取$ Q(x) $为$ f(x) $, 若有$ d_k = -P(x_k)\nabla f(x_k ) = 0 $, 当$ U(x_k)\geq 0 $和$ h(x_k) = 0 $同时成立, 令$ x^* = x_k $, 转步骤5;当存在$ u_j(x_k)<0 $, $ j\in J_0(x) $时, 转步骤4;否则, 进入步骤2.

步骤2  回溯线搜索, 若因可行性恢复跳出, 转步骤1;否则, 进入步骤3.

步骤3  令$ x_{k+1} = x_k(\alpha_k) $, $ k\leftarrow k+1 $, 更新滤子, 转步骤1.

步骤4  令$ J_0(x_k)\setminus\{j\} $, 重新构造$ A(x_k) $, 转步骤1.

步骤5  在$ x^* $处构造填充函数$ T(x, x^*, r) $, 令$ S = 1 $.

步骤6  若$ S>2N $, 则令$ r = \dfrac{r}{10} $, 转步骤12;否则, 令$ k: = 0 $, 取$ x_k\in O(x^*, \delta)\setminus\{x^*\} $.

步骤7  取$ Q(x) $为$ T(x) $, 若$ d_k = 0 $, 令$ F = \{x^*\} $, $ S = S+1 $, 转步骤6;否则进入步骤8.

步骤8  回溯线搜索, 若因可行性恢复跳出, 转步骤11;否则, 进入步骤9.

步骤9  令$ x_{k+1} = x_k(\alpha _k) $, $ k\leftarrow k+1 $, 更新滤子.

步骤10  若$ |F| = 1 $, 则令$ k: = 0 $, 转步骤1;否则, 进入步骤7.

步骤11  若$ |F|>G $, 则令$ F = \{x^*\} $, $ S = S+1 $, 转步骤6;否则, 转步骤7.

步骤12  若$ r<r_0 $, 则停止, $ x^* $为全局最优; 否则, 进入步骤5.

下面讨论该算法的相关性质.

引理4.1   矩阵$ P(x) $是半正定矩阵, 且满足$ P(x)^2 = P(x) $.

定理4.1   (ⅰ)当$ x_k\in X $且非KKT点时, 由算法得到的$ d_k $满足

$ \nabla Q(x_k)^{\rm T} d_k<0, \ \ \ a_j^{\rm T} d_k = 0, j\in J_0(x_k); $

(ⅱ)当$ x_k\notin X $时, 由算法得到的$ d_k $满足$ a_j^{\rm T} d_k<0, j\in J_0(x_k) $.

证    (ⅰ)显然, 当$ x_k\in X $且非KKT点时, 有

$ \nabla Q(x_k)^{\rm T} d_k = -\nabla Q(x_k)^{\rm T} P(x)\nabla Q(x) = -\|P(x)\nabla Q(x)\|^2<0. $

另外

$ A(x_k)^{\rm T} d_k = -A(x_k)^{\rm T} (I-A(x) (A(x)^{\rm T} A(x))^{-1} A(x)^{\rm T})\nabla Q(x_k) = 0, $

即$ a_j^T d_k = 0 $, $ j\in J_0(x_k) $, 结论得证.

(ⅱ)当$ x_k\notin X $时,

$ A(x_k)^{\rm T} d_k = -A(x_k)^{\rm T} P(x_k)\nabla Q(x_k)+\rho_k A(x_k)^{\rm T} B(x_k)^{\rm T}\omega = \rho_k\omega<0, $

即$ a_j^{\rm T} d_k<0 $, $ j\in J_0(x_k) $, 结论得证.

定理4.2  回溯线搜索有限终止.

证   (ⅰ)若$ x_k\in X $时, $ \alpha_k^{\min} = 0 $.由定理4.1知在非KKT点处, $ \nabla Q(x_k)^T d_k<0 $, 则必有$ m_k(\alpha_k)<0 $, 并且$ [-m_k(\alpha_k)]^{s_1}(\alpha_k )^{1-s_1}>0 = \delta_1 [h(x_k)]^{s_2} $, 即(3.5)式成立.

另一方面, 当$ \alpha_k $足够小时, 有

$ \begin{eqnarray*} Q(x_k(\alpha_k))-Q(x_k) & = &Q(x_k)+\alpha_k\nabla Q(x_k)^T d_k+o(\alpha_k d_k)-Q(x_k)\\ & = &\alpha_k\nabla Q(x_k)^T d_k+o(\alpha_k d_k)\leq \delta_2 \alpha_k \nabla Q(x_k)^T d_k, \end{eqnarray*} $

即(3.6)式成立.因此若$ x_k\in X $时, 算法必能找到合适的$ \alpha_k $被接受.

(ⅱ)若$ x_k\notin X $时, $ \alpha_k^{\min}>0 $, 此时要么在步骤3产生一个可被接受的新迭代, 要么进入可行性恢复阶段.

综上, 算法必能找到合适的$ \alpha_k $, 回溯线搜索总会有限终止.

为了证明算法性质, 本文给出如下假设

[A4]存在$ M_m>0 $, 使得$ |m_k(\alpha)|\leq M_m \alpha $.

定理4.3   若假设都成立, 则有

$ \begin{equation} \lim\limits_{k\rightarrow\infty}h(x_k) = 0. \end{equation} $

(4.5)

证 (ⅰ)若滤子集仅在有限步扩大, 假设$ K $, 当迭代$ k>K $时, 滤子集不再扩充.则由算法可知, 对所有$ k\geq K $, (3.5)和(3.6)式都成立, 则

$ \delta_1 [h(x_k)]^{s_2}<[-m_k (\alpha_k)]^{s_1} \alpha_k^{1-s_1}\leq M_m^{s_1} \alpha_k. $

故$ \alpha_k>\dfrac{\delta_1 [h(x_k)]^{s_2}}{M_m^{s_1}} $.由$ 1-\dfrac{1}{s_1}>0 $, 有

$ Q(x_{k+1})-Q(x_k)\leq \delta_2 m_k(\alpha_k)<-\delta_2 \delta_1^{\frac{1}{s_1}}(\alpha_k)^{1-\frac{1}{s_1}}[h(x_k)]^{\frac{s_2}{s_1}}<-\delta_1 \delta_2 M_m^{1-{s_1}} [h(x_k)]^{s_2}. $

故对$ i = 1, 2, \cdots $,

$ Q(x_{K+i}) = Q(x_K)+\sum\limits_{k = K}^{K+i-1}[Q(x_{k+1})-Q(x_k)]<Q(x_K)-\delta_1\delta_2 M_m^{1-{s_1}}\sum\limits_{k = K}^{K+i-1}[h(x_k)]^{s_2}. $

由$ Q(x_{K+i}) $下有界可知, (4.5)式成立.

(ⅱ)若滤子集一直扩大, 下降搜索时所得序列$ \{x_{K_i}\} $, 令$ \{x_{k_i}\} $是其子序列, 滤子集在$ k_i\in \{K_i\} $扩大, 对所有$ i $, 不妨设存在常数$ C_1\in R $和$ C_2>0 $, 使得$ Q(x_{k_i} )\geq C_1, \; h(x_{k_i} )\leq C_2. $

由定理4.1和可行性恢复可知, 必有$ \lim\limits_{i\rightarrow\infty}h(x_{k_i}) = 0 $成立, 则(4.5)式成立.

下面验证滤子的相关性质.

定理4.4   设$ x_k \in S_1(x^*) $, 当前滤子集为$ F_k $.若$ r>0 $充分小, 则一定存在$ \alpha >0 $, 使得$ x_{k+1} = x_k+\alpha d $被滤子$ F_k $接受.

证   由定理4.1可知, 当$ x_k\in X $且不是问题(3.1)的KKT点时, 由算法得到的$ d_k $满足

$ \nabla Q(x_k)^{\rm T} d_k<0;\ a_j^{\rm T} d_k = 0, \ j\in J_0(x), $

所以当$ x_k\in S_1(x^*) $时, 算法中的迭代方向一定是$ f(x_k) $或$ T(x_k, x^*, r) $的下降方向.根据算法的下降性以及滤子的定义, 必有$ x_l\in F_k $, 使得$ f(x_{k+1})< f(x_l) $或$ T(x_{k+1}, x^*)< T(x_l, x^*) $成立, 即$ x_{k+1} $能被滤子$ F_k $接受.

下面讨论滤子集$ F_k $和低水平集$ S_2(x^*) $的关系.

定理4.5   设$ x_k $处的滤子集为$ F_k = \{x_l\; |\; 1\leq l\leq k\} $.如果点$ x_{k+1} $被$ F_k $接受, 并且支配所有$ x_l\; (\forall x_l \in F_k) $, 则$ x_{k+1}\in S_2(x^*) $.

证   若$ x_{k+1}\notin X $, 则对$ x_l\in F_k\bigcap X $, $ x_{k+1} $必然无法支配, 故$ x_{k+1}\in X $.若存在$ x_l\in F_k\bigcap S_2 (x^* ) $, 因为$ x_{k+1} $支配$ x_l $, 结论显然成立.否则, 对所有$ x_l\in F_k\bigcap(S_1 (x^* )\bigcup\{x^* \}) $, 根据滤子的定义, 点$ x^* $肯定在$ F_k $内.而$ x_{k+1} $被$ F_k $接受且支配$ x^* $, 由不等式(4.1)-(4.3)可知$ f(x_{k+1} )<f(x^* ) $, 即$ x_{k+1}\in S_2 (x^* ) $.

定理4.6   设$ x_k $处的滤子集为$ F_k = \{x_l\; |\; 1\leq l\leq k\}\subseteq S_1(x^*) $, 且$ S_2(x^*)\neq \emptyset $.则必存在点$ x_N \in S_2(x^*) $使得$ \forall x_l \in F_k $, $ x_N $支配$ x_l $.

证   首先, $ \forall x_N \in S_2(x^*) $, 有$ f(x_N)< f(x_l)( \forall \; x_l \in F_k\bigcap S_1(x^*) $).

又由于(2.2)式和$ S_2 (x^* ) \neq \emptyset $, 一定存在$ x_N \in S_2(x^*) $, 使得$ {\rm e}^{-\frac{1}{r^2} [f(x_N )-f(x^* )+r] }>1 $, 即$ T(x_N, x^*, r)<0 $.因为$ x_l\in F_k\bigcap S_1 (x^* ) $, 显然$ T(x_l, x^*, r)>0 $, 则有$ T(x_N, x^*, r)< T(x_l, x^*, r), \forall \; x_l \in F_k\bigcap S_1(x^*) $.

另外, 根据$ S_1 (x^* ) $和$ S_2 (x^* ) $的定义, $ x_l $和$ x_N $均在可行域内, 则$ h(x_l ) = h(x_N ) = 0 $, 所以$ x_N $支配$ x_l $.

根据定理4.6, 如果滤子集里只有一个元素, 且该元素支配$ x^* $, 则算法进入到$ x^* $的低水平集里.

定理4.7   设$ x_k $处的滤子集为$ F_k = \{x_l\; |\; 1\leq l\leq k\}\subseteq S_2 (x^* ) $, 进一步搜索得$ x^{**}\in L(P) $, 则$ x^{**} $必能支配$ F_k $中所有的点.

证   由$ f(x^{**} )<f(x_l )(\forall x_l\in F_k\bigcap S_2 (x^* )) $得

$ \begin{equation} {\rm e}^{-\frac{1}{r^2} [f(x^{**} )-f(x^* )+r] }>{\rm e}^{-\frac{1}{r^2} [f(x_l )-f(x^* )+r] }, \end{equation} $

(4.6)

则当$ r $充分小时, $ T(x^{**}, x^*, r)<T(x_l, x^*, r) $必能成立, 并且$ h(x^{**}) = h(x_l) = 0 $, 所以$ x^{**} $必支配$ x_l $.

本节主要验证基于梯度投影的广义滤子填充函数算法的有效性.这些算例都是在Matlab2014b运行环境下运行, 处理器为Intel(R) Core(TM) i5-2410M CPU @ 2.30GHZ 2.30GHZ, 系统类型为64位操作系统, 基于x64的处理器.算法的参数设置如下: $ s_1 = 2.5 $, $ s_2 = 1.2 $, $ \delta_1 = \delta_2 = \beta_1 = \beta_2 = \eta = 10^{-6} $, $ r_0 = 1 $, $ r = 10^{-3} $, $ G = 500 $, $ \delta = 10^{-3} $.

算法的数值结果在表 1-4中列出, 其中各个符号定义如下.

$ k $:第$ k $次求解得局部极小点; $ x_k^0 $:进行第$ k $次局部极小点迭代的初始点; $ x_k^* $:第$ k $个局部极小点; $ f(x_k^*) $:第$ k $个局部极小值; $ F_k $:在第$ k $次达到局部极小点时的滤子.

算例5.1

$ \begin{equation*} \begin{aligned} \min & \ \ f(x) = x_1^2 + x_2^2 - \cos (18x_1) - \cos (18x_2), \ \\ {\rm s.t.}\ & x_1+x_2\leq -2, \; \; x_1-5x_2\leq 3.5, \; \; -3\leq x_i\leq 2, i = 1, 2. \end{aligned} \end{equation*} $

取在可行域外的初始点$ x_1^0 = (0.25397, 1.82675)^{\rm T}, $得到全局极小点

$ x^\infty = (-1.38766, -0.69384)^{\rm T}, $

全局极小值为$ f(x^\infty) = 0.42196 $, 总运行时间为8.513秒.

算例5.2

$ \begin{equation*} \begin{aligned} \min\ & \ \ f(x) = \big(\sum\limits_{i = 1}^5 i\cos{[(i+1) x_1+i]}\big)\big(\sum\limits_{i = 1}^5 i\cos{[(i+1) x_2+i]}\big)\\ \ & \ \ \ \ \ \ \ \ \ \ \ +[(x_1+1.42513)^2+(x_2+0.80032)^2], \\ {\rm s.t.}\ & \ \ \ 10x_1+5x_2\leq -10, \; \; 5x_1-10x_2\leq -10, \; \; -10\leq x_i\leq 10, i = 1, 2. \end{aligned} \end{equation*} $

取在可行域外的初始点$ x_1^0 = (-2.37284, 4.58849){\rm ^T} $, 得到全局极小点

$ x^\infty = (-7.70562, -0.80032)^{\rm ^T}, $

全局极小值为$ f(x^\infty ) = -147.26943 $, 总运行时间为4.566秒.

算例5.3

$ \begin{equation*} \begin{aligned} \min\ & \ \ f(x) = -25(x_1-2)^2-(x_2-2)^2 -(x_3-1)^2-(x_4-4)^2-(x_5-1)^2-(x_6-4)^2, \\ {\rm s.t.}\ & \ \ \ x_3+x_4\geq 4, \; \; x_5+x_6\geq 4, \; \; x_1-3x_2\leq 2, \; \; -x_1+x_2\leq 2, \; \; x_1+x_2\leq 6, \; \; x_1+x_2\geq 2, \\& \ \ \ 0 \leq x_1 \leq 6, \; \; 0 \leq x_2 \leq 8, \; \; 1 \leq x_3 \leq 5, \; \; 0 \leq x_4 \leq 6, \; \; 1 \leq x_5 \leq 5, \; \; 0 \leq x_6 \leq 10. \end{aligned} \end{equation*} $

取在可行域外的初始点$ x_1^0 = (3.28329, 0.83175, 0.89576, 1.54505, 5.04430, 1.52569)^{\rm ^T} $, 得到全局极小点$ x^\infty = (5, 1, 5, 0, 5, 10)^{\rm ^T} $, 全局极小值为$ f(x^\infty ) = -310 $, 总运行时间为15.939秒.

算例5.4

$ \begin{equation*} \begin{aligned} \min\ & \ \ f(x) = \sum\limits_{i = 1}^{20} \big[x_i^2 - \dfrac{1}{10}\cos(5\pi x_i)\big], \\ {\rm s.t.}\ & \ \ \ x_i+x_{i+1}\geq 0.5, i = 1, 2, \cdots, 19, \; \; -1\leq x_i\leq 1, i = 1, 2, \cdots, 20. \end{aligned} \end{equation*} $

取在可行域外的初始点

$ \begin{eqnarray*} x_1^0& = &(0.9058, 0.1270, 0.9134, 0.6324, 0.8147, 0.0975, 0.2785, 0.5469, 0.9575, 0.9649, \\ &&0.1576, 0.9706, 0.9572, 0.4854, 0.8003, 0.1419, 0.4218, 0.9157, 0.7922, 0.9595){\rm ^T}, \end{eqnarray*} $

得到全局极小点

$ \begin{eqnarray*} x^\infty& = &(0.4291, 0.0709, 0.4291, 0.0709, 0.4291, 0.0709, 0.4291, 0.0709, 0.4291, 0.0709, \\ &&0.4291, 0.0709, 0.4291, 0.0709, 0.4291, 0.0709, 0.4291, 0.0709, 0.4291, 0.0709)^{\rm ^T}, \end{eqnarray*} $

全局极小值为$ f(x^\infty ) = 0.55285 $, 总运行时间为409.516秒.由于迭代得的局部极小点比较多, 在本文中只列出部分迭代结果.

以上的数值结果说明了基于梯度投影的广义滤子填充函数算法的有效性.算例5.1与算例5.2的最优解在内部, 其中算例5.1首先搜索到边界上的KKT点, 再通过填充函数迭代到低水平集, 并最终找到了目标函数的稳定点.算例5.2首先在可行域内找到稳定点, 再通过两阶段迭代找到最优解.算例5.3的最优解在边界上, 其搜索情况与算例5.1类似.算例5.4的情况较为复杂, 首先从可行域外搜索到了KKT点, 然后在多次迭代后达到最优解.前三个算例由于规模较小, 耗时均较短.而算例5.4的维数较高, 程序的总运行时间相对长了很多.