双层规划是一类具有上下两层递阶结构的系统决策优化问题, 其中上层问题的目标函数和约束条件不仅与上层决策变量有关, 而且还依赖于下层问题的最优解, 而下层问题的最优解又受上层决策变量的影响. 双层规划已被应用到诸多领域中, 比如交通运输^[1]、资源分配^[2]、价格问题、供应链管理^[3]和生产计划^[3]等. 尽管双层规划问题的应用广泛, 但对双层规划问题的求解仍然是非常困难的, 即使是线性双层规划问题, 也是强$ NP- $难问题, 比如Jeroslow ^[4]、Bard ^[5]、Hansen ^[6]、Vicente等人^[7], 给出了相应的证明. 另一个造成双层规划问题求解困难的原因是非凸性, 即使最简单的双层线性规划, 一般来说也是一个非凸优化问题, 尽管能找到问题的解, 通常只能是局部最优解而非全局最优解. 虽然双层规划问题求解困难, 但到目前为止, 仍有不少算法被用来求解该问题, 比如极点搜索算法^[8]、分枝定界算法^{[9, 10]}、罚函数算法^[11]、下降方法、互补旋转算法、非数值优化算法等. 其中绝大多数算法都只能求解双层规划问题的局部最优解, 而本文提出了一种填充函数法来求解双层规划问题的全局最优解.

填充函数法最早是由西安交通大学的葛仁溥教授^[12]首次提出的. 该算法的基本思想是: 先将目标函数局部极小化, 找到原始问题的一个局部最优解, 然后在找到的这个解处构造一个填充函数, 再对此填充函数进行局部极小化, 使得填充函数的局部极小点处的原始问题的目标函数值比之前的局部极小值小, 多次进行上述局部极小化和填充过程, 即可找到全局最优解. 该方法是求解优化问题的一种有效的确定性全局优化算法. 在利用该方法求解优化问题时, 关键是要构造一个填充函数以保证算法在实现过程中, 避免陷入局部最优, 从而找到更好的局部极小点.

填充函数法最初是被用来求解无约束优化问题的^[13--16], 后来又被扩展到用于求解其他优化问题, 比如, 非线性全局优化问题^[17]、等式约束优化问题^{[18, 19]}、非线性方程组等^[20]. 本文对填充函数的应用进行扩展, 用来求解双层规划问题.

本文余下的内容可分为如下几个部分: 第二部分, 双层规划问题的转化, 利用下层问题的KKT条件将双层规划转化为单层约束规划问题; 第三部分, 结合罚函数法, 构造了一种新的填充函数, 并探讨了它的性质; 第四部分, 基于构造的新的填充函数, 提出了一种求解双层规划问题的填充函数法; 第五部分, 对两个例子进行了数值实验; 最后一部分是总结.

本文考虑上、下两层都带有等式和不等式约束的双层规划问题, 其数学模型可描述为:

$ \begin{equation} \begin{array}{ll} \mathop {\min}\limits_x F(x, y) \quad \hbox{s.t.}\\ G_i(x, y)\le 0, i=1, 2, \cdots, m, \quad H_j(x, y)=0, j=1, 2, \cdots, n, \quad x\in X, \\ \mathop {\min}\limits_y f(x, y)\quad \hbox{s.t.}\\ g_a(x, y)\le 0, a=1, 2, \cdots, s, \quad h_b(x, y)=0, b=1, 2, \cdots, t, \quad y\in Y. \end{array} \end{equation} $

(2.1)

其中$ x\in \mathbb{R}^n $是上层问题的决策变量, $ y\in \mathbb{R}^n $是下层问题的决策变量, $ F(x, y)\; , G_i(x, y)\; , H_j(x, y) $分别是上层的目标函数和约束函数, $ f(x, y)\; , g_a(x, y)\; , h_b(x, y) $分别是下层的目标函数和约束函数. 对问题(2.1), 本文做如下假设.

假设2.1   假设$ f(x, y) $, $ g_a(x, y), a=1, 2, \cdots, s. $和$ h_b(x, y), b=1, 2, \cdots, t. $是光滑函数, $ f $和$ g_a(x, y), a=1, 2, \cdots, s. $是凸函数, $ h_b(x, y), b=1, 2, \cdots, t. $是关于$ y $的线性函数.

由假设2.1可知, 利用下层问题的KKT最优性条件, 可将双层规划问题等价转化为连续的单层约束优化问题.

假设2.2  假设问题(2.1)的局部极小点存在且个数有限. 设下层问题的拉格朗日函数为:

$ \begin{equation} L(x, y, \lambda, \mu)=f(x, y)+\sum\limits_{a=1}^s \lambda_a g_a(x, y)+\sum\limits_{b=1}^t \mu_b h_b(x, y). \end{equation} $

(2.2)

则下层问题的KKT最优性条件为:

$ \begin{equation} \begin{array}{ll} \frac{\partial f}{\partial y}+\sum\limits_{a=1}^s \lambda_a\frac{\partial g_a}{\partial y}+\sum\limits_{b=1}^t \mu_b\frac{\partial h_b}{\partial y}=0, \\ g_a(x, y)\le 0, a=1, 2, \cdots, s, \\ \lambda_a\ge 0, a=1, 2, \cdots, s, \\ \lambda_a g_a(x, y)=0, a=1, 2, \cdots, s, \\ h_b(x, y)=0, b=1, 2, \cdots, t, \\ \mu_b\ne 0, b=1, 2, \cdots, t, \\ y\in Y.\\ \end{array} \end{equation} $

(2.3)

那么双层规划问题(2.1)可转化为如下的单层约束规划问题:

$ \begin{equation} \begin{array}{ll} \mathop {\min}\limits_{x, y} F(x, y) \quad \hbox{s.t.}\\ G_i(x, y)\le 0, i=1, 2, \cdots, m, \quad H_j(x, y)=0, j=1, 2, \cdots, n, \\ \frac{\partial f}{\partial y}+\sum\limits_{a=1}^s \lambda_a\frac{\partial g_a}{\partial y}+\sum\limits_{b=1}^t \mu_b\frac{\partial h_b}{\partial y}=0, \\ g_a(x, y)\le 0, a=1, 2, \cdots, s, \\ \lambda_a\ge 0, a=1, 2, \cdots, s, \\ \lambda_a g_a(x, y)=0, a=1, 2, \cdots, s, \\ h_b(x, y)=0, b=1, 2, \cdots, t, \\ \mu_b\ne 0, b=1, 2, \cdots, t, \\ x\in X, \quad y\in Y.\\ \end{array} \end{equation} $

(2.4)

令$ z:=(x, y, \lambda, \mu) $, $ Z:=\big\{z=(x, y, \lambda, \mu)\mid x\in X, y\in Y, \lambda_a\ge 0, \mu_b\ne 0\big\} $, $ L(z)=\lambda_a g_a(x, y) $, $ \wp(z)= \frac{\partial f}{\partial y}+ \sum\limits_{a=1}^s \lambda_a \frac{\partial g_a}{\partial y}+ \sum\limits_{b=1}^t \mu_b \frac{\partial h_b}{\partial y} $. 则单层约束规划问题(2.4)可化为

$ \begin{equation} \begin{array}{ll} \min F(z) \quad \hbox{s.t.}\\ G_i(z)\le 0, i=1, 2, \cdots, m, \quad H_j(z)=0, j=1, 2, \cdots, n, \\ \wp(z)=0, \\ g_a(z)\le 0, a=1, 2, \cdots, s, \quad h_b(z)=0, b=1, 2, \cdots, t, \\ L(z)=0, \\ z\in Z.\\ \end{array} \end{equation} $

(2.5)

此时问题(2.5)等价于问题(2.1), 通过求解问题(2.5)的最优解就可以求解原始问题(2.1)的最优解.

定义3.1 ^[20]    函数$ P(z, z^*):Z\to \mathbb{R} $被称为$ F(z) $在局部极小点$ z^* $处的填充函数, 如果它满足如下性质:

(1) $ z^* $是$ P(z, z^*) $在$ Z $上的一个严格局部极大点;

(2) 对任意的$ z\in S_1 $, 有$ \nabla P(z, z^* ) \ne 0 $, 其中$ S_1=\big\{z|F(z)\ge F(z^*), z\in Z|\{z^*\}\big\} $;

(3) 若$ z^* $不是全局极小点, 即$ F(z) $存在另一个局部极小点$ z^{**}\in {\rm int}(Z) $, 满足$ F(z^{**})<F(z^*) $, 则$ P(z, z^*) $存在一个极小点$ z'\in N(z^{**}, \delta)\subset Z $. 其中$ N(z^{**}, \delta) $是点$ z^{**} $处的$ \delta $邻域.

在下文中, 介绍满足定义3.1的填充函数$ F_f(z, z^*, q, u) $. 首先, 本文设计了两个不同的连续函数$ \phi_r(t) $和$ \varphi_r(t) $. 函数$ \phi_r(t) $具有如下的性质: 当$ t\ge0 $时, $ \phi_r(t)=1 $; 当$ t\leq-r $时, $ \phi_r(t)=0 $; 另外它是单调递增的. $ \varphi_r(t) $具有如下的性质: 当$ t\ge r $时, $ \varphi_r(t)=t+r $; 当$ t\leq 0 $时, $ \varphi_r(t)=0 $; 另外它是单调递增的.

对于任何给定的$ r>0 $,

$ \begin{equation} \phi_r(t) = \left\{ {\begin{array}{*{20}c} {1, } & {t \ge 0, } \\ {\log_{1+r^3 }\ (-2(t+r)^3+3r(t+r)^2+1), } & { - r < t < 0, } \\ {0, } & {t \le - r.} \\ \end{array}} \right. \end{equation} $

(3.1)

$ \begin{equation} \varphi_r(t) = \left\{ {\begin{array}{*{20}c} {t+r, } & {t \ge r, } \\ { (t+r)\log_{1+r^3 }\ (-2t^3+3rt^2+1), } & { 0 < t < r, } \\ {0, } & {t \le 0.} \\ \end{array}} \right. \end{equation} $

(3.2)

则有

$ \begin{equation} \phi_r^\prime(t) = \left\{ {\begin{array}{*{20}c} {0, } & {t \ge 0, } \\ { \frac{-6(t+r)^2+6r(t+r)}{{-2(t+r)^3+3r(t+r)^2+1}}, } & { - r < t < 0, } \\ {0, } & {t \le - r.} \\ \end{array}} \right. \end{equation} $

(3.3)

$ \begin{equation} \varphi_r^\prime(t) = \left\{ {\begin{array}{*{20}c} {1, } & {t \ge r, } \\ {\log_{1+r^3 }\ (-2t^3+3rt^2+1)+(t+r) \frac{-6t^2+6rt}{{-2t^3+3rt^2+1}} , } & { 0 < t < r, } \\ {0, } & {t \le 0.} \\ \end{array}} \right. \end{equation} $

(3.4)

显然$ \phi_r(t) $和$ \varphi_r(t) $在$ \mathbb{R} $上是连续可微的. 当$ r\rightarrow0^+ $时, $ \varphi_r(t)\rightarrow max(t, 0) $. 为此, 本文构造了一个填充函数:

$ \begin{equation} \begin{split} F_f(z, z^*, q, u)=&\frac{1}{{||z - z^* ||^2+1}}\phi_q\bigg (F(z)-F(z^*)+\frac{1}{2}u\bigg(\sum\limits_{i=1}^m \varphi_\frac{1}{u}^2(G_i(z))\\ &+\sum\limits_{a=1}^s \varphi_\frac{1}{u}^2(g_a(z))+\sum\limits_{j=1}^n H_j^2(z)+ \sum\limits_{b=1}^t h_b^2(z)+L^2(z)+\wp^2(z)\bigg)\bigg). \end{split} \end{equation} $

(3.5)

其中$ q>0, u>0 $是两个参数, $ u>0 $是惩罚参数, $ \phi_r(t) $和$ \varphi_r(t) $由式(3.1)、(3.2) 定义.

以下定理表明$ F_f(z, z^*, q, u) $是满足定义3.1的一类填充函数.

定理3.1   假设$ q>0 $充分小, 且$ u>0 $足够大, 则$ z^* $是$ F_f(z, z^*, q, u) $在$ Z $上的一个全局极大点. 证明: 由式(3.5)可得,

$ \begin{equation} \begin{split} F_f(z^*, z^*, q, u)=&\frac{1}{{||z^* - z^* ||^2+1}}\phi_q\bigg(F(z^*)-F(z^*)+\frac{1}{2}u\bigg(\sum\limits_{i=1}^m \varphi_\frac{1}{u}^2(G_i(z^*))\\ &+\sum\limits_{a=1}^s \varphi_\frac{1}{u}^2(g_a(z^*))+\sum\limits_{j=1}^n H_j^2(z^*)+ \sum\limits_{b=1}^t h_b^2(z^*)+L^2(z^*)+\wp^2(z^*)\bigg)\bigg). \end{split} \end{equation} $

(3.6)

又根据$ \phi_r(t) $的性质知, 当$ t\ge0 $时, $ \phi_r(t)=1 $, 有$ F_f(z, z^*, q, u)\leq F_f(z^*, z^*, q, u) $. 因此$ z^* $是$ F_f(z, z^*, q, u) $在$ Z $上的一个全局极大点.

定理3.1表明函数$ F_f(z, z^*, q, u) $满足定义3.1的条件(1). 一般情况下, 要寻求函数的全局极大点是非常困难的事情. 然而, 此处构造的函数$ F_f(z, z^*, q, u) $的全局极大点非常容易求得, 就为$ z^* $.

定理3.2   假设$ q>0 $充分小, 且$ u>0 $足够大, 则对于任何$ z\in S_1\setminus\{z^*\} $有$ \nabla F_f(z, z^*, q, u)\ne 0 $, 其中$ S_1=\big\{z\in S \mid F(z)\ge F(z^*)\big\} $.

证明: 假设$ z\in S_1 $, 则$ F(z)>F(z^*) $, 得出

$ F_f(z, z^*, q, u)=\frac{1}{{||z - z^* ||^2+1}}, \; \; \; \; \; \nabla F_f(z, z^*, q, u)=-\frac{2(z-z^*)}{\big({||z - z^* ||^2+1}\big)^2}. $

又$ z\; \ne z^* $则$ \nabla F_f(z, z^*, q, u)\ne 0 $.

定理3.2表明函数$ F_f(z, z^*, q, u) $满足定义3.1的条件(2).

定理3.3   假设$ q>0 $充分小, 且$ u>0 $足够大, 若$ z^* $不是全局极小点, 即$ F(z) $存在另一个极小点$ z^{**}\in int(Z) $, 使得$ F(z^{**})<F(z^*) $, 则当$ 0<q<\big(F(z^*)-F(z^{**})\big) $时, $ F_f(z, z^*, q, u) $一定存在一个极小点$ z' $, 使得$ z'\in N(z^{**}, \delta)\subset Z $, 其中$ N(z^{**}, \delta) $是点$ z^{**} $处的$ \delta $邻域.

证明: 因$ F(z) $是连续的, 则一定存在点$ \bar{z}^* $, $ z'\in N(z^{**}, \delta) $, 满足$ F(\bar{z}^*)=F(z^*) $, $ F(\bar{z}^*)>F(z')>F(z^{**}) $, 其中$ z'=\bar{z}^*+\varepsilon(\bar{z}^*-z^*) $. 当$ 0<q<\big(F(z^*)-F(z^{**})\big) $时, 令$ F(z')\leq F(z^*)- q $, 则$ F_f(z', z^*, q, u)=0 $. 又$ F_f(z, z^*, q, u)\ge 0 $对所有的$ z\in Z $都成立. 所以$ z' $是$ F_f(z, z^*, q, u) $的极小点, 则结论成立.

定理3.3表明函数$ F_f(z, z^*, q, u) $满足定义3.1的条件(3).

本节基于上一节构造的填充函数$ F_f(z, z^*, q, u) $, 提出了一种求解双层规划问题的填充函数法, 下面介绍相应的求解问题(2.5)的填充函数算法4.1.

步0   选择充分小的正数$ q $, $ \lambda_L $和充分大的正数$ u $, 选择一个正的整数$ K $和方向$ e_i $, $ i=1, \ldots, K $. 选择一个初始点$ z_0\in Z $. 置$ k := 0 $;

步1   从$ z_k $开始, 通过局部搜索方法找到问题(2.5)的一个局部极小点$ z_k^* $;

步2   令

$ \begin{equation} \begin{split} F_f(z, z_k^*, q, u)=&\frac{1}{{||z - z_k^* ||^2+1}}\phi_q\bigg(F(z)-F(z_k^*)+\frac{1}{2}u\bigg(\sum\limits_{i=1}^m \varphi_\frac{1}{u}^2(G_i(z))\\ &+\sum\limits_{a=1}^s \varphi_\frac{1}{u}^2(g_a(z))+\sum\limits_{j=1}^n H_j^2(z)+ \sum\limits_{b=1}^t h_b^2(z)+L^2(z)+\wp^2(z)\bigg)\bigg). \end{split} \end{equation} $

(4.1)

其中$ \phi_r(t) $和$ \varphi_r(t) $由式(3.1)、(3.2) 定义, 置$ l=1 $和$ \lambda=1 $;

步3

(a) 若$ l\le K $, 转(b); 否则, 转步5;

(b) 若$ \lambda\ge \lambda_L $, 则置$ w_k^l:=z_k^*+\lambda e_l $, 转(c); 否则, 置$ l:=l+1 $, $ \lambda=1 $, 转(a);

(c) 若$ w_k^l\in Z $, 转(d); 否则, 置$ \lambda:= \frac{\lambda}{2} $, 转(b);

(d) 若$ f(w_k^l)<f(z_k^*) $, 则置$ z_{k+1}:=w_k^l $, $ k:=k+1 $, 转步1; 否则, 转步4;

步4   从$ w_k^l $出发, 搜索如下填充函数问题的局部极小点:

$ \begin{equation} \mathop {\min }\limits_{z \in Z } F_f(z, z_k^* , q, u), \end{equation} $

(4.2)

设$ w_k^* $是问题(4.2)一个已获得的局部极小点, 若它满足$ f(w_k^*)<f(z_k^*) $, 则置$ z_{k+1}:=w_k^* $, $ k:=k+1 $, 转步1; 否则, 置$ l=l+1 $, $ \lambda=1 $, 转步3(a);

步5   令$ z_s=z^*_k $, 停止.

算法4.1的主要思想:

该算法按以下方法选取步0中的方向$ e_i $, $ i=1, \ldots, K $.例如, 当$ n=2 $时, 取$ K=6n $, 方向$ e_i $被选为:

$ e_i = \left( {\cos \frac{{2\pi (i - 1)}}{K}, \sin \frac{{2\pi (i - 1)}}{K}} \right), i = 1, \ldots , K. $

当$ n\ge 3 $时, 取$ K=2n $, $ i=1, \ldots, n $. $ e_i $中的第$ i $个分量为1, 其他分量为0; $ i=n+1, \ldots, K $. $ e_i $中的第$ i $个分量为-1, 其他分量为0.

在步0中选取一个合适的初始点, 从这个初始点开始, 利用局部极小化方法, 可以获得问题(2.5)的一个局部极小点$ z_k^* $. 接下来的任务是确定$ z_k^* $是全局极小点或者找到比$ F(z_k^*) $函数值更小的局部极小点.

步3是选取合适的初始点, 去极小化填充函数问题(4.2).若方向集$ \{e_i, i=1, \ldots, K\} $足够大, 当无论怎么选取初始点时, 问题(4.2)的局部极小点$ w_k^* $均在$ Z $的边界处获得, 且$ F_f(w_k^*, z_k^*, q, u)>0 $, 则称步1中找到的局部极小点$ z_k^* $为问题(2.5)的近似全局极小点. 当$ q $取得足够小时, 该点即被认为是问题(2.5)的全局极小点.

定理4.1   算法4.1产生的序列$ \{z_k^*\} $收敛于问题(2.5)的全局极小点.

证明：由假设2.2可知, 问题(2.1)的局部极小点存在且个数有限, 同时算法4.1中产生的迭代序列$ \{z_k^*\} $满足$ F(z_{k+1}^*)<F(z_k^*) $, 所以$ \{z_k^*\} $是一个使得函数值下降的序列, 由此可知, 算法4.1产生的序列$ \{z_k^*\} $收敛于问题(2.5)的全局极小点.

本节应用填充函数算法来求解两个数值例子, 以揭示算法的可行性, 并与文献中的数值结果进行比较. 所有的数值实验均在Matlab2018a中执行. 算法使用Matlab2018a工具箱中的fmincon函数来局部极小化问题. 在整个数值实验中, 算法中的参数为: $ q=10^{-6} $, $ u= 4^5 $, $ \lambda_L= \frac{1}{2^5} $.

例1见文献[22].

$ \begin{equation} \begin{array}{ll} \mathop {\min}\limits_x F=x^2+y^2-16x-5xy \quad \hbox{s.t.}\quad 0\leq x\leq 20, \\ \mathop {\min}\limits_y f=-y\quad \hbox{s.t.}\quad x+y-20\leq 0, \quad 0\leq y\leq 10. \end{array} \end{equation} $

(5.1)

该问题的全局最优解为$ (x, y)=(11.1429, 8.8571)^T $, 其对应的目标函数值为$ -469.1429 $. 例$ 1 $的数值计算结果见表 1. 由表 1可知, 通过代入不同的初始点, 算法经过一次迭代, 即可找到全局最优解. 数值结果表明该算法是可行的.

例2见文献[24].

$ \begin{equation} \begin{array}{ll} \mathop {\min}\limits_{x_1, x_2} F=\begin{pmatrix}x_1&x_2\end{pmatrix}\begin{pmatrix}2&3\\4&1\end{pmatrix} \begin{pmatrix}y_1\\y_2\end{pmatrix} \quad \hbox{s.t.}\quad x_1+x_2=1, \quad x_1, x_2\ge 0, \\ \mathop {\min}\limits_{y_1, y_2}f=\begin{pmatrix}x_1&x_2\end{pmatrix}\begin{pmatrix}-1&-4\\-3&-2\end{pmatrix} \begin{pmatrix}y_1\\y_2\end{pmatrix} \quad \hbox{s.t.}\quad y_1+y_2=1, \quad y_1, y_2\ge 0.\\ \end{array} \end{equation} $

(5.2)

该问题的全局最优解为$ (x_1, x_2, y_1, y_2)=(0.0048, 0.0576, 0.1057, 0.8943)^T $, 其对应的目标函数值为0.0899. 例$ 2 $的数值计算结果见表 2. 由表 2可知, 通过六次迭代, 即可找到全局最优解. 数值结果表明该算法是可行的.

本文根据定义3.1, 结合罚函数提出了一种新的填充函数并探讨了它的性质. 通过下层问题的KKT条件将双层规划转化为单层约束规划问题; 最后基于构造的新的填充函数, 提出了一种求解双层规划问题的填充函数法, 并对两个例子进行了数值实验, 结果表明填充函数算法是可行的. 但在算法中, 由于初始点选取的不同, 就会出现找到最优解需要的迭代次数不同, 简言之, 算法对于初值的选取太过于依赖, 这在后续的研究中会进一步改进.