续志明 ^1,2 , 刘匡宇 ¹ , 白宇 ¹ , 王书宁 ¹
1. 清华大学 自动化系, 信息科学与技术国家实验室, 北京 100084;
2. 空军工程大学 理学院, 西安 710051

收稿日期：2016-11-15
基金项目：国家自然科学基金资助项目（61473165，61134012）；国家“九七三”重点基础研究项目（2012CB720505）
作者简介：续志明(1983-), 男, 博士研究生
通信作者：王书宁, 教授, E-mail:swang@tsinghua.edu.cn

摘要：连续分片线性规划是一类应用广泛的重要规划，寻找连续分片线性规划的全局最优解是研究这类规划的重点和难点。该文研究的是一种对此类规划进行全局寻优的确定性启发式算法。由于此类规划问题可以转化为凸多面体上的凹优化问题进行求解，因此利用凹函数的上水平集的凸性，该文提出可以通过直接穿透目标函数上水平集在其等值面上进行搜索，以逃离当前局部最优解进行全局寻优。该方法中每次逃离的搜索方向都通过山形凹目标函数的顶点投影来确定，因此称为山顶投影穿山法。在数值实验中，将所提出的山顶投影穿山法与CPLEX以及绕山法进行了比较，结果表明该算法在计算速度与全局寻优能力上性能优越。
关键词：全局优化分片线性凹优化割平面穿山
Hill tunneling method via peak subpoints for continuous piecewise linear programming
XU Zhiming^1,2, LIU Kuangyu¹, BAI Yu¹, WANG Shuning¹
1.National Laboratory for Information Science and Technology, Department of Automation, Tsinghua University, Beijing 100084, China;
2.Science College, Air Force Engineering University, Xi'an 710051, China


Abstract: Continuous piecewise linear (CPWL) optimization is a widely used mathematical programming method, and searching for the globally optimal solution is the research emphasis. This paper presents a heuristic algorithm with determinacy for the global optimization of CPWL programming, which is referred to as the hill tunneling via peak subpoint (HTPS) algorithm. A CPWL programming problem can be transformed into an equivalent concave optimization problem over a polyhedron; thus, the current algorithm makes use of the concavity of the super-level sets of concave piecewise linear functions to escape a local optimum to search on the other side of its contour surface by cutting across the super-level set. Each searching path for the hill tunneling is established using the peak subpoint of the "hill" shaped concave objective function. Numerical tests show that this algorithm is more efficient and has better global search capability than CPLEX or the hill detouring (HD) method, which shows its superior performance.
Key words: global optimizationpiecewise linearconcave minimizationcutting plane methodhill tunneling
连续分片线性规划是一类重要的数学规划问题。一方面由于其在全局上的灵活性，连续分片线性问题可以逼近任意连续非线性问题，从而比线性规划用途更为广阔；另一方面由于其在局部上的线性，相比于一般的非线性问题，有其结构和解法上的优势，可以利用线性规划很快地找到局部最优解，这使得连续分片线性规划具有广泛的应用。
在连续分片线性规划中，若目标函数与约束条件均为凸的，则其局部最优解就是全局最优解，应用分片线性单纯形法^[1-3]等方法相对容易解决。但对一般的分片线性规划问题，寻找全局最优已被证明是一个NP-hard问题^[4]。
寻求连续分片线性规划的全局最优解，目前的常规做法是通过引入整数变量将问题转化为混合整数规划(mixed integer programming, MIP)来处理^[5]。另一类成为国内外****研究热点的方法是启发式算法。实际上考虑到问题的规模，这类算法具有更好的实用性。这方面的研究方向包括随机性的智能搜索算法^[6]和确定性的算法^[7-8]。智能搜索算法由于本身原理上的随机性，算法性能无法很好地得到保证，所以本文将注意力主要集中于确定性的启发式算法。
确定性算法中特别值得一提的是Huang等^[7]提出的绕山法(hill detouring, HD)。该方法提出了一种先离开再进入可行域、沿目标函数等值面进行搜索的全局寻优策略，效果良好。受此启发，本文设计了对连续分片线性规划进行全局寻优的山顶投影穿山法(hill tunneling via peak subpoint, HTPS)。这种方法从某一局部最优点出发，通过直接穿透(而非如绕山法那样绕开)目标函数的上水平集，以更快地逃离当前局部最优陷阱，并到达其等值面的另一侧进行搜索。
1 预备工作1.1 问题的等价转化连续分片线性规划问题的一般表示为如下的条件极小化问题：

$\min \left\{ {{f_0}\left( z \right)\left| {{\rm{s}}{\rm{.t}}{\rm{.}}\;{f_i}\left( z \right) \le 0,i = 1,2, \cdots ,N} \right.} \right\}.$

(1)

其中f_i(z), i=0, 1, …, N均为连续分片线性函数。
根据文[8]，连续分片线性规划总有精确罚函数存在，即总存在某个数${\hat r}$＞0，使得对任意r≥${\hat r}$，式(1)与下面的无约束问题具有相同的局部最优解集：

$\min {F_r}\left( z \right) = {f_0}\left( z \right) + r\sum\limits_{i = 1}^N {\max \left\{ {0,{f_i}\left( z \right)} \right\}} .$

(2)

文[9]指出，广义链接超平面(generalized hinging hyperplanes, GHH)模型具有完全表示能力，所以式(2)中的F_r(z)可以表示成如下的GHH形式：

$p\left( z \right) = \sum\limits_{k = 1}^{\hat N} {\mathop {\max }\limits_{1 \le j \le {n_k}} \left\{ {\mathit{\boldsymbol{\hat a}}_{kj}^{\rm{T}}z + {{\hat b}_{kj}}} \right\}} - \sum\limits_{k = 1}^{\bar N} {\mathop {\max }\limits_{1 \le j \le {l_k}} \left\{ {\mathit{\boldsymbol{\bar a}}_{kj}^{\rm{T}}z + {{\bar b}_{kj}}} \right\}} .$

(3)

其中所有max函数都称为GHH模型的基函数。这样就可以将模型式(2)化为DC规划(different of convex functions)minp(z)的形式。
进一步地，通过引入新变量λ_k，minp(z)又可以等价转化为如下的条件极小化问题：

$\begin{array}{l}\min \;\;\;\sum\limits_{k = 1}^{\hat N} {{\lambda _k}} - \sum\limits_{k = 1}^{\bar N} {\mathop {\max }\limits_{1 \le j \le {l_k}} \left\{ {\mathit{\boldsymbol{\bar a}}_{kj}^{\rm{T}}z + {{\bar b}_{kj}}} \right\}} ,\\{\rm{s}}{\rm{.}}\;{\rm{t}}{\rm{.}}\;\;\;\mathit{\boldsymbol{\hat a}}_{kj}^{\rm{T}}z + {{\hat b}_{kj}} \le {\lambda _k},1 \le j \le {n_k},1 \le k \le \hat N.\end{array}$

(4)

式(4)的目标函数是一个凹函数，所以这是一个凸多面体上的凹分片线性优化问题。这意味着任意连续分片线性规划问题式(1)都可以等价地化为具有线性约束的凹分片线性优化问题：

$\min \left\{ {f\left( x \right) = \sum\limits_{i = 1}^M {\mathop {\min }\limits_{1 \le j \le {N_i}} \left\{ {\mathit{\boldsymbol{a}}_{ij}^{\rm{T}}x + {b_{ij}}} \right\}} \left| {{\rm{s}}{\rm{.t}}{\rm{.}}\;Cx} \right. \le g} \right\}.$

(5)

其中：$C \in {\mathbb{R}^{m \times n}},g \in {\mathbb{R}^m},{a_{ij}} \in {\mathbb{R}^n},{b_{ij}} \in \mathbb{R},\forall i,j$。为方便起见，可以将本问题的可行域记为Ω。
特别的，凹优化问题式(5)的目标函数在x∈R²时图像形似一座山，相关算法的命名灵感正来源于此。
本文之后的分析都将围绕问题式(5)展开。不失一般性，以下假定Ω的顶点均为非退化的，并且Ω是封闭的有界集。事实上，若Ω无界，则可根据问题的实际背景添加约束b_l≤x≤b_u，把Ω化为有界集。
1.2 若干定义与性质1) 定义。
任意给定一点$x \in {\mathbb{R}^n}$，式(5)的目标函数值为

$f\left( x \right) = \sum\limits_{i = 1}^M {\mathop {\min }\limits_{1 \le j \le {N_i}} \left\{ {\mathit{\boldsymbol{a}}_{ij}^{\rm{T}}x + {b_{ij}}} \right\}} .$

(6)

此时?1≤i≤M，可定义关联于点x的指标集Φ_xⁱ为

$\mathit{\Phi }_x^i = \left\{ {j\left| {\mathit{\boldsymbol{a}}_{ij}^{\rm{T}}x + {b_{ij}}} \right. = \mathop {\min }\limits_{1 \le j \le {N_i}} \left\{ {\mathit{\boldsymbol{a}}_{ij}^{\rm{T}}x + {b_{ij}}} \right\}} \right\}.$

(7)

进而定义在给定点x处的指标空间Φ_x为

${\mathit{\Phi }_x} = \mathit{\Phi }_x^1\mathit{\Phi }_x^2 \cdots \mathit{\Phi }_x^M.$

(8)

定义1?对任意一个给定点x∈Rⁿ，任取φ=(k₁，k₂，…，k_M)∈Φ_x，令a$a_x^\varphi = \sum\limits_{i = 1}^M {{a_{i{k_i}}}} ,b_x^\varphi = \sum\limits_{i = 1}^M {{b_{i{k_i}}}} $, 则可以得到一个线性函数：

$\mathit{\boldsymbol{a}}_x^{\varphi {\rm{T}}}x + b_x^\varphi .$

(9)

称此线性函数为式(5)的目标函数在给定点x处的一个线性子面，记为Ψ_x^φ。
定义2?对任意γ∈R，称集合{x|f(x)=γ}为式(5)的目标函数的γ-等值面，记为Γ_γ。
定义3?若对式(5)的某一等值面Γ_γ与某一子面Ψ：a^Tx+b，有Γ_γ∩{x|f(x)=a^Tx+b}≠?，则定义相应的等值子面为

${\mathit{\Gamma }_{\gamma ,\mathit{\Psi }}} = \left\{ {x\left| {f\left( x \right)} \right. = \gamma ,f\left( x \right) = {\mathit{\boldsymbol{a}}^{\rm{T}}}x + b} \right\}.$

(10)

除上述定义外，本文后面还会涉及到所谓γ-扩张的概念，根据文[10]其定义如下。
定义4?若f(x)是凹函数，对某给定点x，设标量γ满足γ≤f(x)，θ₀是一个足够大的正数。对方向$d \in {\mathbb{R}^n}\backslash \left\{ 0 \right\},$令$\theta = {\rm{min}}\left\{ {{\theta _0},{\rm{sup}}\left\{ {t:f\left( {x + td} \right) \geqslant \gamma } \right\}} \right\}$，则称

$x' = x + \theta d,$

(11)

为f从x向d的γ-扩张。
2) 性质。
对分片线性凹函数f(x)，如下2条性质成立。
性质1?若f(x)=a^Tx+b，?x∈D在某开集$D \subset {\mathbb{R}^n}$上成立，则有f(x)≤a^Tx+b，?$x \in {\mathbb{R}^n}$。
性质2?对任意$\hat x \in {\mathbb{R}^n}$，f(x)的上水平集$\left\{ {x \in {\mathbb{R}^n}|f\left( x \right) \geqslant f\left( {\hat x} \right)} \right\}$是一个凸多面体。
这2条性质已经在文[7]中被证明。
2 跳出局部最优解的穿山法思想对式(5)全局寻优时面临的主要挑战是：当获得了问题的某一局部最优解之后，如何跳出当前局部最优陷阱，以继续搜索其他可能存在的更优解。
绕山法在这方面进行了有益的尝试，其逃离当前局部最优陷阱的方法为：以问题在当前局部最小解x₀附近的某个等值子面上的不可行点x₁为起点，一直沿等值面绕行凹目标函数“山”，直到等值面再次进入原问题可行域，然后重新开始搜索一个更优的局部最小解，如图 1所示。

图 1 绕山法示意图

图选项

直观上，逃离局部最优点时，更有效率的做法应该是直接去到距离当前局部最优点较远的区域进行搜索，而非绕山法这样从当前局部最优点附近开始用“绕”的方式逐渐逃离进行搜索。基于这种考虑，本文设计了一种新的逃离局部最优解的策略——穿山，如图 2所示。

图 2 穿山法示意图

图选项

穿山策略是从当前局部最优点x₀出发，沿指向目标函数“山”的山顶在自变量空间的投影点P的方向，直接穿透目标函数上水平集，快速到达较远处的目标函数“山”的另一侧开始搜索。这种策略的逃离方向相当于从山顶下方穿过，可以利用山体最大程度屏蔽当前局部最优解的影响, 所以称这种方法为山顶投影穿山法。这种方式有可能快速摆脱当前局部最优解的影响，从而实现加速逃离的目的。
3 山顶投影穿山法3.1 算法总体框架由于有很多成熟的下降算法^[7]可用来求式(5)的局部极小解，因此对算法框架的描述从获得了式(5)的某一局部最小顶点x₀开始。
1) 初始化当前迭代轮数k=1，当前局部最小值γ=f(x₀)，以及Ω₀←Ω，x₁←x₀。
2) 在运算时间超出给定阈值T前迭代如下操作：
a) 从x_k出发进行γ-扩张，得到n个扩张点x₁^k，x₂^k，…，x_n^k；
b) 用x₁^k，x₂^k，…，x_n^k确定一个超平面H_k，以H_k作为割平面，将Ω_k－1割去一部分不含优于x₀的可行点的区域，得到一个新的更小的可行域Ω_k={x|C^kx≤g^k}；
c) 若Ω_k=?，x₀即全局最优，算法终止，否则求目标函数在Ω_k内的山顶投影点P_k；
d) 从x_k出发沿方向P_k-x_k穿山，得到目标函数在山另一侧的某γ-等值子面Γ^k；
e) 检验Γ^k或其附近的等值子面是否重新进入了可行域；
f) 若已重新进入，则跳出迭代，令Ω←Ω_k，并从重新进入可行域处再次开始求解局部最优顶点，否则计算Γ_f(x0)上在Γ^k及其附近离Ω_k最近的点x_k+1；
g) k←k+1。
3.2 算法实现细节与技术分析1) 第一轮迭代时的γ-扩张。
在第一轮迭代时k=1，此时γ-扩张的起点x₁就是当前局部最优顶点x₀，方向分别为从x₀指向其在可行域上各相邻顶点的方向。换句话说，第一轮迭代时的γ-扩张，就是从x₁=x₀出发，沿Ω₀=Ω过x₁的各棱进行直线搜索，直到与目标函数等值面Γ_f(x0)相交，各交点即为第一轮迭代的扩张点。
在x₀非退化的情形下，会有n条棱满足条件，所得的扩张点x_i¹，i=1，2，…，n使得各x_i¹-x₁均仿射独立，从而x₁与x₁¹，x₂¹，…，x_n¹可定义一个n维单纯形S₁。
下面给出第一轮迭代n个扩张方向的具体求法。
由于x₀是Ω的一个(非退化的)顶点，所以问题式(5)的约束条件里会有n个不相关的在x₀处起作用(取“=”)，记这n个约束为

$\mathit{\boldsymbol{B}}x \le \mathit{\boldsymbol{\bar b}}.$

(12)

则x₀=B^-1b，其中B=(b₁，b₂，…，b_n)^T∈$\mathbb{R}$^n×n，b=(b₁，b₂，…，b_n)^T∈$\mathbb{R}$ⁿ。
命题1?对式(12)中的B，若B^-1=(d₁，d₂，…，d_n)，则-d₁，-d₂，…，-d_n就是从x₀出发沿棱指向其在可行域的n个相邻顶点的方向。
证明：由BB^-1=I可得，对任意1≤i，j≤n及i≠j，有b_j^Td_i=0，b_i^Td_i=1。由于b_j是超平面{x|b_j^Tx=b_j}的法向量，因此d_i//{x|b_j^Tx=b_j}对任意i≠j成立。从而

${\mathit{\boldsymbol{d}}_i}//\bigcap\limits_{\begin{array}{*{20}{c}}{j = 1}\\{j \ne 1}\end{array}}^n {\left\{ {x\left| {\mathit{\boldsymbol{b}}_j^{\rm{T}}x = {{\bar b}_j}} \right.} \right\}} ,\;\;\forall 1 \le i \le n.$

(13)

这说明d_i平行于各超平面{x|b_j^Tx=b_j}，j≠i所交成的棱。
另一方面，对任意t>0，有

$\begin{array}{*{20}{c}}{\mathit{\boldsymbol{b}}_i^{\rm{T}}\left( {{x_0} + t{\mathit{\boldsymbol{d}}_i}} \right) = \mathit{\boldsymbol{b}}_i^{\rm{T}}{x_0} + t\mathit{\boldsymbol{b}}_i^{\rm{T}}{\mathit{\boldsymbol{d}}_i} = \mathit{\boldsymbol{b}}_i^{\rm{T}}{x_0} + t > }\\{\mathit{\boldsymbol{b}}_i^{\rm{T}}{x_0} = {{\bar b}_i}.}\end{array}$

(14)

即x₀+td_i不满足约束式(12)，从x₀出发的方向d_i是指向可行域外的。这说明-d_i是从x₀出发沿棱指向可行域内的，即指向可行域内x₀的相邻顶点的方向。这样就得到了第一轮γ-扩张的n个方向。
2) 割平面的引入。
第k轮迭代中，γ-扩张的起点x_k与所得的n个扩张点x₁^k，x₂^k，…，x_n^k可以确定一个单纯形S_k。由于Ω_k-1与S_k均为凸集，再结合目标函数f(x)的凹性，易得Ω_k-1∩S_k内已没有比当前局部最优点x₀更优的点。
启发算法可以在每一轮迭代中将可行域割掉一部分而不影响寻求原问题的全局最优。具体操作中是通过添加割平面来实现这个目的的。
作为n维单纯形S_k的n个顶点，x₁^k，x₂^k，…，x_n^k必然可以唯一确定一个n-1维超平面H_k：c_k^Tx=g_k，这就是第k轮迭代所产生的割平面。
显然x_k作为γ-扩张的起点，与集合Ω_k-1\S_k分别在H_k的两侧。根据割平面的方程在Ω_k-1中添加新约束，即可得到一个新的更小的可行域：

$\begin{array}{l}{\mathit{\Omega }_k} = \left\{ {x\left| {\begin{array}{*{20}{c}}{{C^{k - 1}}x \le {g^{k - 1}}}\\{\left[ {{\rm{sign}}\left( {\mathit{\boldsymbol{c}}_k^{\rm{T}}{x_k} - {g_k}} \right)} \right]\left( {c_k^{\rm{T}}{x_k} - {g_k}} \right) \le 0}\end{array}} \right.} \right\} = \\\;\;\;\;\;\;\;\left\{ {x\left| {{C^k}x \le {g^k}} \right.} \right\}.\end{array}$

(15)

其中sign(·)为符号函数。
添加割平面H_k后，若出现Ω_k=?，则表明整个可行域内已没有比当前局部最优点x₀更优的点了。此时算法终止，x₀即为原问题的全局最优解。
3) 山顶投影点的计算。
目标函数f(x)在Ω_k内的山顶投影点P_k，即为如下的凸优化问题的解：

$\begin{array}{l}\max \;\;\;f\left( x \right) = \sum\limits_{i = 1}^M {\mathop {\min }\limits_{1 \le j \le {N_i}} \left\{ {\mathit{\boldsymbol{a}}_{ij}^{\rm{T}}x + {b_{ij}}} \right\}} ;\\{\rm{s}}{\rm{.}}\;{\rm{t}}{\rm{.}}\;\;\;\;{C^k}x \le {g^k}.\end{array}$

(16)

问题式(16)可以转化为线性规划问题：

$\begin{array}{l}\min \;\;\;\;\sum\limits_{i = 1}^M {{\lambda _i}} ;\\{\rm{s}}{\rm{.}}\;{\rm{t}}{\rm{.}}\;\;\;\;{C^k}x \le {g^k},\\\;\;\;\;\;\;\;\;\; - a_{ij}^{\rm{T}}x - {b_{ij}} \le {\lambda _i},\\\;\;\;\;\;\;\;\;\;1 \le i \le M,1 \le j \le {N_i}.\end{array}$

(17)

若式(17)的解为${\left[ {{{\mathit{\boldsymbol{\tilde x}}}^{\rm{T}}},{{\tilde \lambda }_1},{{\tilde \lambda }_2}, \cdots ,{{\tilde \lambda }_M}} \right]^{\rm{T}}}$，则第k轮迭代的山顶投影点即为${P_k} = \tilde x$。
4) 穿山操作与线性子面的确定。
在第k轮迭代获得了山顶投影点P_k后，从x_k出发沿方向P_k-x_k进行γ-扩张，这就是第k轮迭代中的穿山操作。
第k轮穿山所得的扩张点x′_k位于目标函数等值面Γ_γ上。任取某一φ∈Φ_x′k，通过式(9)得到线性子面Ψ_x′^φk，以此作为第k轮穿山得到的线性子面Ψ_k，并进而获得相应的等值子面，作为第k轮穿山得到的等值子面Γ^k。即令Ψ_k←Ψ_x′^φk，并相应的令Γ^k←Γ_γ，Ψk，其中记Ψ_k为

${\mathit{\Psi }_k}:a_k^{\rm{T}}x + {b_k}.$

(18)

5) 对穿山所得等值子面的检验。
穿山操作后，需对本轮所得的等值子面Γ^k进行检验，以判断Γ_γ在此附近是否重新进入了当前可行域Ω_k，即判断在Γ^k或其附近是否有Γ_γ∩Ω_k≠?。
对Γ^k的检验通过解如下线性规划问题来实现：

$\begin{array}{l}\min \;\;\;\;s;\\{\rm{s}}{\rm{.}}\;{\rm{t}}{\rm{.}}\;\;\;\;\;{C^k}x \le {g^k} \le S,\\\;\;\;\;\;\;\;\;\;\;\mathit{\boldsymbol{a}}_k^{\rm{T}} + {b_k} \le f\left( {{x_k}} \right).\end{array}$

(19)

其中：s∈R，S=(s₁，s₂，…，s_n)^T；C^k、g^k来自于式(15)中描述Ω_k的方程组；a^T_kx+b_k为式(18)中Ψ_k的定义方程。
设式(19)的最优解为$\left( {\hat s,\hat x} \right)$，则${\hat s}$度量了Γ^k与Ω_k在某种意义上的距离，如图 3所示。

图 3 对等势子面的检查示意图

图选项

式(19)的计算结果有以下3种可能。
a) 检查结果一：$\hat s > 0,\mathit{\boldsymbol{a}}_k^{\rm{T}}\hat x + {b_k} = f\left( {{x_k}} \right)$。
条件$\hat s > 0$表示Γ^k∩Ω_k=?，条件$\mathit{\boldsymbol{a}}_k^{\rm{T}}\hat x + {b_k} = f\left( {{x_k}} \right)$表示${\hat x}$∈Γ_f(xk)，如图 3所示。此时应以${\hat x}$作为下一轮γ-扩张的起始点x_k+1。
b) 检查结果二：$\hat s > 0,\mathit{\boldsymbol{a}}_k^{\rm{T}}\hat x + {b_k} ＜ f\left( {{x_k}} \right)$。
这种情况表明${\hat x}$?Γ_γ，且根据性质1易知有$f\left( {\hat x} \right) ＜ f\left( {{x_k}} \right)$。
此时需要将Γ^k更新为一个其附近的距离Ω_k更近的等势子面，并重新进行检验。具体来说，以目标函数f(x)在${\hat x}$处的一个线性子面${\psi _{\hat x}}:a_{\hat x}^{\rm{T}}\hat x + {b_{\hat x}}$作为新的Ψ_k，并重新定义相应的新Γ^k，然后对此新Γ^k重新求解式(19)。由于目标函数上水平集{x|f(x)≥f(x_k)}是凸集，因此这种更新操作必然可以得到一个更小的最优${\hat s}$。
重复这种更新与检查，直到出现${\hat s}$≤0或${\hat s}$>0，$a_k^{\rm{T}}\hat x + {b_k} = f\left( {{x_k}} \right)$。
c) 检查结果三：${\hat s}$≤0。
此时Γ^k∩Ω_k≠?，${\hat x}$即为交点，可以以${\hat x}$为起点重新开始寻找局部最优。
6) k≥2时的γ-扩张。
当k≥2时，若第k-1轮迭代中对等值子面的检验出现检查结果一，则要开始第k轮的迭代。此时，取第k-1轮迭代时对等值子面进行检查所得的${\hat x}$，作为第k轮迭代中γ-扩张的起点x_k；而γ-扩张的方向将通过下述方式来确定。
在非退化假设下，式(19)的约束条件在最优解$\left( {\hat s,\hat x} \right)$处有n+1个起作用(取到“=”)。由于检查结果一意味着$a_k^{\rm{T}}\hat x + {b_k} = f\left( {{x_k}} \right)$，因此跟s有关的剩余约束不等式里会有n个在$\left( {\hat s,\hat x} \right)$处取“=”，记这n个约束为

$\mathit{\boldsymbol{B}}x - \mathit{\boldsymbol{\tilde g}} \le \mathit{\boldsymbol{\tilde S}}.$

(20)

类似命题1的分析，若记B^-1=(d₁，d₂，…，d_n)，则-d₁，-d₂，…，-d_n分别平行于可行域Ω_k-1的距Γ^k-1最近的顶点指向其n个相邻顶点的方向。从而-d₁，-d₂，…，-d_n就是进行本轮γ-扩张的n个方向。
4 数值试验本部分将通过仿真试验对比各算法的全局寻优性能，试验环境为Core i3 3.30 GHz CPU和8 GB内存，HTPS与HD使用商业软件为Matlab 2010a，求解等价MIP使用商业软件CPLEX 12.4。试验中作为比较对象的其他文献方法，其算例结果均是由本文再现了其算法程序后在此计算环境下仿真所得，而非引自相关文献。
4.1 试验方案记待测问题集P由若干子集P_n，m组成，每个P_n，m包含20个随机产生的形如下式的问题，每个问题N=n且M=m：

$\min \;\;\;\sum\limits_{i = 1}^M {{\beta _i}\mathop {\min }\limits_{1 \le j \le {N_i}} \left\{ {\mathit{\boldsymbol{a}}_{ij}^{\rm{T}}x + {b_{ij}}} \right\}} .$

(21)

其中：β_i∈{1，－1}，a_ij∈$\mathbb{R}$^N，b_ij∈$\mathbb{R}$，?i，j，并设各N_i的选取上界为L，即N_i≤L。为方便起见，不失一般性，将所有问题都限定在x∈[0，1]^N上，并令L=3。
同时记待测算法集S={HTPS，HD，MIP}，试验中分别考察S中3类算法在待测问题集P上的求解性能。其中MIP表示将问题式(21)转化为等价的MIP问题^[5]进行求解，HD和HTPS分别表示将式(21)转化为等价的凹优化问题式(5)后用HD及HTPS算法进行求解。就具体的凹转化而言，?p∈P，首先将待测问题的目标函数看做式(3)中的p(z)，然后以节1.1的方式将式(21)化为式(4)的形式，再给其添加约束x∈[0，1]^N，这样待测问题就转化为了式(5)的形式。
试验中，首先用一个N=2、M=30的实例验证穿山法的可行性及展示其具体求解过程，然后再考察和比较以S中3类算法分别求解问题集P中各类问题的性能。试验中参数N从{3，5，8，10，30，50，80，100}中取值，参数M从{10，30，50，80，100，120，150，180，200}中取值。所有试验的时间上限阈值都设定为T=500 s。
4.2 穿山法演示验证为验证穿山法的可行性，首先以前述方式生成一个N=2、M=30的定义在x∈[0，1]²上的具体算例，如图 4所示，然后演示用HTPS算法对其求解的具体过程。

图 4 演示验证算例目标函数

图选项

图 5在自变量空间上展示了穿山法求解此算例的具体过程。从初始可行点(0，0)出发，寻求局部最优得到x₁，并陷入其中。此时通过穿山迭代，跳出此局部最优陷阱，到达与x₁具有相同函数值的非局部最小点x₂。再次寻求局部最优得到x₃，然后穿山到达与x₃具有相同函数值的x₄。最终从x₄出发搜索得到了全局最优解x₅。

图 5 穿山法求解过程演示

图选项

这个二维实例上的测试验证了穿山法的效果，并直观展示了其跳出局部最优陷阱的能力。
4.3 各算法仿真试验与结果分析1) 运算速度的比较。
通过计算以算法s∈S求解P_n，m内各问题的平均运算时间t_s，来考察其运算速度，如图 6和7所示。

图 6 各算法在N=5时的平均运算时间

图选项

图 7 各算法在M=50时的平均运算时间

图选项

可以看出，在问题规模较小时3种算法的平均运算时间差别不大。随着问题规模的增长，t_MIP迅速大幅增长，这是由于求解MIP的分支定界法就本质而言也属于一种穷举算法。同时，t_HD与t_HTPS增长缓慢，并且在N不大时二者基本处于同一量级。
2) 优化效果的比较。
为使对各算法优化效果的评价在问题的不同规模下具有一致性，参考文[11]引入几个相对性能指标。以s∈S求解p∈P时，基于所解得的最优值v_p，s的函数E_p，s=v_p，s-min{v_p，s^*：s^*∈S}+1，分别定义相对优化率PR和优胜率SR如下：

${\rm{P}}{{\rm{R}}_{p,s}} = \frac{{{E_{p,s}}}}{{\min \left\{ {{E_{p,{s^ * }}}:{s^ * } \in S} \right\}}},$

(22)

${\rm{S}}{{\rm{R}}_s} = \frac{{{\rm{card}}\left\{ {p \in P:{E_{p,s}} \le {E_{p,{s^ * }}},\forall {s^ * } \in S} \right\}}}{{{\rm{card}}P}}.$

(23)

其中“card”是用来求集合的势的函数。
图 8和9通过计算SR值比较了各算法对不同规模问题的相对优化性能。

图 8 各算法在N=5时的SR值

图选项

图 9 各算法在M=50时的SR值

图选项

可以看出，除N≤10时MIP的优化效果略微优于HD外，其余维数下MIP的优化效果最差，而HTPS在所有维数下都是性能最好的。
5 结论本文针对连续分片线性规划的全局寻优问题，设计了山顶投影穿山法，从某一局部最优点出发，直接穿透目标函数上水平集，到达其等值面的另一侧进行搜索，以实现快速逃离当前局部最优陷阱的目的。本文设计的穿山法只用到了线性规划与直线搜索等成熟技术，并且在过程中通过添加割平面进一步提高了搜索效率。在数值实验中，与CPLEX和绕山法相比，穿山法在面对各种规模的问题时，在计算速度与全局寻优能力上都性能优越。

参考文献

[1]	Fourer R. A simplex algorithm for piecewise-linear programming I:Derivation and proof[J]. Mathematical Programming, 1985, 33(2): 204–233. DOI:10.1007/BF01582246
[2]	Fourer R. A simplex algorithm for piecewise-linear programming Ⅱ:Finiteness, feasibility and degeneracy[J]. Mathematical Programming, 1988, 41(1-3): 281–315. DOI:10.1007/BF01580769
[3]	Fourer R. A simplex algorithm for piecewise-linear programming Ⅲ:Computational analysis and applications[J]. Mathematical Programming, 1992, 53(1-3): 213–235. DOI:10.1007/BF01585703
[4]	Keha A B, De Farias J I R, Nemhauser G L. Abranch-and-cut algorithm without binary variables for nonconvex piecewise pinear optimization[J]. Operations Research, 2006, 54(5): 847–858. DOI:10.1287/opre.1060.0277
[5]	Vielma J P, Ahmed S, Nemhauser G. Mixed-integer models for nonseparable piecewise linear optimization:Unifying framework and extensions[J]. Operations Research, 2010, 58(2): 303–315. DOI:10.1287/opre.1090.0721
[6]	Xi X, Xu J, Mu X, et al. Continuous piecewise linear programming via concave optimization and genetic algorithm[C]//2012 IEEE 51st Annual Conference on Decision and Control (CDC). Maui, HI, USA:IEEE, 2012:2509-2514.
[7]	Huang X, Xu J, Mu X, et al. The hill detouring method for minimizing hinging hyperplanes functions[J]. Computers & Operations Research, 2012, 39(7): 1763–1770.
[8]	Huang X, Xu J, Wang S. Exact penalty and optimality condition for nonseparable continuous piecewise linear programming[J]. Journal of Optimization Theory and Applications, 2012, 155(1): 145–164. DOI:10.1007/s10957-012-0032-7
[9]	Wang S, Sun X. Generalization of hinging hyperplanes[J]. Information Theory, IEEE Transactions on, 2005, 51(12): 4425–4431. DOI:10.1109/TIT.2005.859246
[10]	Horst R, Tuy H. Global Optimization:Deterministic Approaches[M]. Berlin: Springer, 1996.
[11]	Dolan E D, More J J. Benchmarking optimization software with performance profiles[J]. Mathematical Programming, Series B, 2002, 91(2): 201–213. DOI:10.1007/s101070100263