设施选址问题主要研究如何确定设施的最优数目和最佳位置以及如何为顾客分配设施。选址决策属于战略层决策，建设投资大、周期长、回收慢，且一经选定就需长期运营，因此，合理选址就显得十分重要。经典的选址模型^[1]，如覆盖问题、p-中值问题、无容量限制的设施选址问题，都是基于空间垄断的假设，即将要建立的服务设施将提供一种独特的产品或服务，并且独占这部分市场。然而，在大量的现实问题中这一模型假设不再成立：在新建设施进入市场之前，通常已有一个或多个竞争设施存在；在新建设施进入市场之后，还会有其他竞争设施也尝试进入这个市场。因此，在实际的研究中需考虑对手的竞争。尤其是在加油站、快餐店、超市、快递门店等盈利性设施的选址中，考虑对手的竞争显得尤为重要^[2]。

竞争选址问题指的是市场上存在2个或2个以上的同类产品或服务的提供者，通过空间选址和设施设计达到最大化市场份额或最大化利润的目标。竞争选址问题明确地考虑到竞争对手已经(或将要)进入目标市场，而新建设施将与它们竞争既有的市场份额。在对客户进行服务的过程中，竞争双方争夺客户需求，增加自身市场占有率，且任一方的决策将对竞争双方均产生作用^[3]。作为竞争选址问题的一种，带预见性的竞争是指在新建设施进入市场之前，已经考虑到竞争对手将紧随其后进入市场。此类问题被描述为Stackelberg博弈^[4]，即双寡头垄断市场的主从博弈。在该博弈中，有两个行动者(领导者和跟随者)相继进行决策。领导者能预见到跟随者的最优反应，且跟随者能在自己决策之前观测到领导者的决策。Hakimi^[5]首次将这个模型引入到选址领域，领导者将跟随者的最优反应纳入到决策过程中，并采取最小最大策略来维护自己的利益。文献[5]提出2类问题：medianoid和centroid，其中，(r|X_p)-medianoid problem(RXMP)是指在领导者从候选设施集合X_p中确定p个设施点之后，跟随者如何确定自己将要新建的r个设施的选址，从而使自己的市场份额最大。而(r|p)-centroid problem(RPCP)指的是领导者在选址时已经知道跟随者将会在自己行动之后新建r个设施，领导者需确定自己的p个设施的选址，从而使自己的市场份额最大。

RPCP考虑的是一个离散的网络，客户需求只能被最近的设施捕获，领导者和跟随者的目标都是最大化自己的市场份额。RPCP通常都建模成为一个二元双层线性规划问题^[6-8](binary bilevel linear pro- gramming problem，BBLP)。Hakimi^[4]证明了RPCP是NP难问题，Noltemeier等^[9]进一步证明了其是 ∑ ​ p 2 难问题。这是因为验证领导者的解是否可行需求解一个NP难问题，即最优化跟随者的策略。为求解这个高度复杂的问题，本文采用了一种新型的建模方法，即将跟随者的策略用一种启发式方法(贪婪算法)代替，然后将这种策略的结果以约束条件的形式，集成到领导者的问题中，这样初始的双层规划问题就可以作为一个单层问题来求解。

1 RPCP的研究现状

令I和J分别表示客户和候选设施的集合，每个客户有确定的需求(或购买力)μ_i，d_ij表示客户i∈I到设施j∈J的距离(或运输成本、接受服务的惩罚成本)。对 ∀ i∈I，j、k∈J，当d_ij<d_ik时，令参数a_ijk=1，否则取0。引入下面的二元决策变量：

x_j=1，如果领导者开放设施j；

y_j=1，如果跟随者开放设施j；

z_i=1，如果客户i的需求被领导者的设施捕获。

当领导者的决策为x时，跟随者仍能捕获客户i的需求的设施集合为

J i ( x ) = { j ∈ J | d i j ＜ min l ∈ J d i l ｜ x l = 1 } 。

(1)

于是，RPCP可以建模为如下的二元双层线性规划模型^[6-8]

max ∑ ​ i ∈ I μ i z i 。

(2)

s.t.

∑ ​ j ∈ J x j = p ;

(3)

x j ∈ { 0 , 1 } ， j ∈ J 。

(4)

其中，z是跟随者优化问题的最优解：

max ∑ ​ i ∈ I μ i ( 1 − z i ) 。

(5)

s.t.

∑ ​ j ∈ J y j = r ；

(6)

1 − z i ≤ ∑ ​ j ∈ J i ( x ) y j , ∀ i ∈ I ;

(7)

y j ∈ { 0 , 1 } , 0 ≤ z i ≤ 1 ， ∀ i ∈ I , j ∈ J 。

(8)

上下层问题的目标函数(2)和(5)分别最大化领导者和跟随者的市场份额。这是一个零和博弈，所以约束(5)等价于 min ∑ ​ i ∈ I μ i z i ，即最小化领导者的市场份额。约束(3)和(6)分别限制领导者和跟随者开放p个、r个设施。约束(7)确定了服务每个客户i的设施。约束(4)和(8)设置了变量的取值范围，变量z可以松弛为连续类型，因为在最优解中会自动取0或者1。

下文中，将双层规划模型(2)~(8)称为RPCP，对给定的领导者决策x，为验证是否可行，需要精确求解跟随者问题(5)~(8)，记为RPCP-O。

近年来，很多算法(包括启发式算法和精确算法)都对“离散选址问题标杆库”^①中的算例进行了测试，规模可以达到100个客户、100个候选设施、p=r=20。启发式算法运用了p-中值问题的邻域操作(如概率交换)，在很小的迭代次数内，可以达到较高的精度；但是每一次迭代的复杂度仍然很高，因为需要精确求解RPCP-O来评估领导者的邻域解是否可行。截止到目前，比较有效的启发式方法是Davydov等^[7]提出的可变邻域搜索算法和随机禁忌搜索算法，Biesinger等^[8]设计的混合遗传算法，Alekseeva等^[10]提出的混合文化基因算法，以及Davy- dov^[11]设计的拉格朗日松弛-禁忌搜索算法。启发式算法虽可以在较短时间内找到较好的结果，但无法判断和最优解的间隙有多大。

精确算法，包括Alekseeva等^[12]提出的迭代精确法、Roboredo和Pessoa^[13]提出的分支切割法(RP-B&C)以及Alekseeva和Kochetov^[6]提出的改进迭代精确法(modified iterative exact method，MEM)，可以保证全局最优，但时间复杂度太高(当p=r=10时，求解时间超过10 h^[13])。Alekseeva等^[12]提出的迭代精确法是基于一个单层的二元优化模型，有着指数级的约束条件和决策变量；每次迭代中，都添加跟随者新的选址方案，然后再次求解此模型。Roboredo和Pessoa^[13]提出的RP-B&C与迭代精确法类似，但只有多项式级的决策变量，并且指数级的约束条件可以升维成更强的不等式。Alekseeva和Koche- tov^[5]提出的MEM改进了迭代精确法，只有多项式级的决策变量，并运用了文献[13]提出的强有效不等式。

由于启发式算法无法判断找到的可行解与最优解的间隙，而精确算法求解双层模型的时间复杂度又太高。因此，本文将RPCP建立成一个单层模型，用约束条件来表示跟随者的贪婪策略，从而避免精确求解跟随者的优化模型。数值实验表明所提出的模型框架对小规模和中等规模算例均有效；基于这个框架，下一步能设计相应的精确算法，从而更有效地求解这一类问题。

2 跟随者采用贪婪策略的单层模型 2.1 决策变量

假设给定领导者的设施，跟随者采取贪婪策略，逐个开放r个设施，记这一新模型为RPCP-G。建模方法如下。决策变量x、z的定义如前文所示，而对变量y重定义，并引入新的二元变量：

y_jt=1，如果跟随者在贪婪算法的第t次迭代中开放设施j；

s_jt=1，如果客户i的需求，被跟随者在第t次迭代中开放的设施所捕获；

v_ijt=1，如果客户i在贪婪算法的第t次迭代中，仍然可以被设施j捕获。

Daskin在最大覆盖问题中介绍了贪婪算法^[1]。引入到本文中，该算法的步骤如下。在初始阶段，给定领导者的决策向量x，更新变量v的值。在第t次迭代中，若领导者的所有设施都比j到客户i的距离更远，那么j就可以捕获到i(v_ijt=1)；否则不能捕获(v_ijt=0)。对每个关闭的设施，计算若将其开放可以捕获到的总需求量( ∑ ​ i ∈ I μ i v i j t )，然后跟随者选择开放可以捕获最大需求的设施(j)。被这个设施捕获到(s_it=1)的客户，在后续迭代中不再被考虑( v i j t ′ = 0 , ∀ t ′ ＞ t )。重复以上步骤，直到跟随者已经选择了r个设施。

本文选择贪婪策略作为跟随者响应的原因主要有2点。首先，贪婪算法广泛运用于设施选址问题中^[1]，定理3将证明当给定领导者决策，跟随者的最佳策略(模型RPCP-O的解)和贪婪策略(模型RPCP-G的解)相同时，那么就找到了原问题的最优解。其次，这种方法让跟随者的决策以多项式级的约束条件嵌入到领导者问题中成为可能。

① Discrete Location Problems benchmark library，见http://www.math.nsc.ru/ AP/benchmarks/Competitive/p_med_comp_eng.html

2.2 单层模型

max ∑ ​ i ∈ I μ i z i 。

(2)

s.t.

∑ ​ j ∈ J x j = p ；

(3)

∑ ​ j ∈ J y j t = 1 ， ∀ t = 1 , … , r ;

(9)

x j + ∑ ​ t = 1 r y j t ≤ 1 ， ∀ j ∈ J ；

(10)

z i + ∑ ​ t = 1 r s i t = 1 , ∀ i ∈ I ;

(11)

1 − v i j 1 ≤ p − ∑ ​ k ∈ J a i j k x k ， ∀ i ∈ I , j ∈ J ；

(12)

v i j 1 + x k ≤ 1 + a i j k , ∀ i ∈ I , j 、 k ∈ J ;

(13)

v i j t ≥ v i j 1 − ∑ ​ t ′ = 1 t − 1 s i t ′ ， ∀ i ∈ I , j ∈ J , t = 2 , ⋯ , r ；

(14)

v i j t ≤ 1 − ∑ ​ t ′ = 1 t − 1 s i t ′ , ∀ i ∈ I , j ∈ J , t = 2 , ⋯ , r ；

(15)

v i j t ≤ v i j , t − 1 ， ∀ i ∈ I , j ∈ J , t = 2 , ⋯ , r ；

(16)

∑ ​ i ∈ I μ i v i j t ≤ ∑ ​ i ∈ I μ i s i t , ∀ j ∈ J , t = 1 , ⋯ , r ;

(17)

s i t ≥ y i t + v i j t − 1 ， i ∈ I , j ∈ J , t = 1 , ⋯ , r ；

(18)

s i t ≤ 1 − y i t + v i j t , i ∈ I , j ∈ J , t = 1 , ⋯ , r ;

(19)

x j , y j t ∈ { 0 , 1 } ， j ∈ J , t = 1 , ⋯ , r ；

(20)

0 ≤ z i , s i t , v i j t ≤ 1 , i ∈ I , j ∈ J , t = 1 , ⋯ , r 。

(21)

约束(9)强调了跟随者在每次迭代中需要开放一个设施。约束(10)指每个设施只能被两个决策者开放一次。约束(11)指当算法终止时，每个客户只能被一个决策者捕获到(没有服务距离的限制)。

领导者先行决策，当跟随者行动“可以被捕获”可以通过下面的3个规则来进行说明。

1) 在t=1次迭代中，对任意客户i和任意设施j，若领导者的所有设施到i的距离都比j要远，那么j可以捕获到客户i，即只有当 ∑ ​ k ∈ J a i j k x k = p (对∀k∈J，当xk=1时，都有a_ijk=1)时，才有v_ij1=1，如约束(12)所示；反之，若领导者至少有一个设施比j近，跟随者就不能捕获到客户i，即当存在x_k=1且a_ijk=0时，有v_ij1=0，如约束(13)所示。

2) 在t=2, …, r次迭代中，若客户i最初可以被设施j捕获( v i j 1 = 1 )，且至今仍未被跟随者捕获( ∑ ​ t ′ = 1 t − 1 s i t ′ = 0 )，那么在第t次迭代中，仍然可以被捕获，即v_ijt=1，如约束(14)所示；该客户一旦被捕获( ∑ ​ t ′ = 1 t − 1 s i t ′ = 1 )，那么就不能再被捕获到，即v_ijt=0，如约束(15)所示。

3) 不能被给定设施捕获的客户，将一直不能被捕获，如约束(16)所示。

贪婪策略是通过约束条件(17)实现的。在每一次迭代中，跟随者捕获到的总需求，至少等于单个设施所能捕获到的最大需求。在第t次迭代中，若跟随者开放了设施j(即y_jt=1)，当v_ijt=1 时，客户i可以被捕获(即s_it=1)，如约束(18)所示；反之，当v_ijt=0时，客户i不能被捕获(即s_it=0 )，如约束(19)所示。

最后，约束条件(20)和(21)定义了变量的特性，z、s和v可以松弛为连续变量，将在定理1进行证明。

3 模型的分析

定理1 (变量松弛原理)。当变量z、s和v被松弛为连续类型后，在最优解中仍一定取整数值。

证明：在任意可行解中，约束(20)将变量x和y限制为取值0或1。对任意的客户i和设施j，若领导者开放了设施k(即x_k=1)，且k比j到客户i的距离更近(即a_ijk=0)，那么由约束(13)，一定有v_ij1=0；反之，若领导者开放的所有设施，都比j到i的距离更远，即当 ∑ ​ k ∈ J a i j k x k = p 时，j在第一次迭代中，可以捕获到客户i，那么由约束(12)，一定有v_ij1=1。因此，对 ∀ i, j，变量v_ij1一定取0或者1。

约束条件(9)表示跟随者在每次迭代中，都选择一个设施。不失一般性，设第t次迭代中开放的设施为j即y_jt=1，若y_ijt=1，约束(18)保证了s_it=1；若y_ijt=0，约束(19)保证了s_it=0。即当y_ijt取值0或1时，s_it也取0或者1。

当t>1时，由约束条件(16)，有 v i j t ≤ v i j , t − 1 ≤ 1 ，因此，当 v i j , t − 1 = 0 时，一定有 v i j t = 0 ；反之若 v i j , t − 1 =1 ，当 s i , t − 1 = 0 时，约束(14)暗示有v_ijt=1，当 s i , t − 1 = 1 时，约束(15)暗示有v_ijt=0。因此，当 v i j , t − 1 和 s i , t − 1 取0或者1时，v_ijt也一定取0或1。

从二元变量v_ij1开始，在贪婪算法中，迭代地运用上面的结论，可以保证对 ∀ i, j, t，变量v_ijt、s_it一定取值为0或1。

最后，目标函数(2)是最大化变量z_i的线性函数，约束(11)有 z i = 1 − ∑ ​ t = 1 r s i t ，因此保证了z_i也一定取值0或者1。

定理2 　(模型正确性原理)。模型RPCP-G的解，对应于跟随者的贪婪策略。

证明：在每一次迭代t中，跟随者都会开放能够捕获到最多客户需求的设施。即 y j t = 1 ⇒ j = argmax j ′ ∑ ​ ​ i ∈ I μ i v i j t 。当捕获最大客户需求的设施唯一时，则反之也成立。否则，一定存在某个 j ′ ≠ j 使得 ∑ ​ i ∈ I μ i v i j ′ t ＞ ∑ ​ i ∈ I μ i v i j t 。于是根据式(17)有 ∑ ​ ​ i ∈ I μ i s i t ＞ ∑ ​ i ∈ I μ i v i j t 。然而，若y_jt=1，根据式(18)和式(19)有s_it=v_ijt，于是 ∑ ​ i ∈ I μ i s i t = ∑ ​ i ∈ I μ i v i j t ，产生矛盾。

另外，在每次迭代中，开放的设施j将捕获所有能捕获到的客户，并且已经被捕获的客户在后续迭代中将不能再次被捕获。在这种策略中，跟随者将开放r个设施，其结果构成了贪婪算法的解。

首先精确求解模型RPCP-G，此时跟随者采用贪婪策略，令X^G和Y^G(X^G)分别表示领导者和跟随者的解；然后给定X^G，精确求解模型RPCP-O，得到跟随者的解为Y^*(X^G)。假设2种情况下，领导者捕获的客户需求分别是f(X^G, Y^G(X^G)和f(X^G, Y^*(X^G)。又设当精确求解双层规划模型RPCP时，它们的解分别是X^*和Y^*(X^*)。最优解表示为f(X^*, Y^*(X^*)。一般的，有下面的结论成立。

定理3 　(模型最优性条件)。 f ( X G , Y * ( X G ) ) ≤ f ( X * , Y * ( X * ) ) ≤ f ( X G , Y G ( X G ) ) ，即给定领导者的决策X^G，模型RPCP-G和RPCP-O的最优值分别提供了原问题的上界和下界。

证明：解(X^G, Y^*(X^G)是原问题的可行解，因为当领导者采取决策X^G时，跟随者实际会做出的决策即是F^*(X^G)。原问题是最大化问题，即最大化领导者的市场份额，显然任意可行解提供了一个有效下界，于是有 f ( X G , Y * ( X G ) ) ≤ f ( X * , Y * ( X * ) ) 。

RPCP是一个零和博弈，因此当跟随者的策略不是最优时，领导者的市场份额将增加，于是有 f ( X * , Y * ( X * ) ) ≤ f ( X * , Y G ( X * ) ) 。又因X^G是模型RPCP-G的最优解，所以有 f ( X * , Y G ( X * ) ) ≤ f ( X G , Y G ( X G ) ，联立有 f ( X * , Y * ( X * ) ) ≤ f ( X G ， Y G ( X G ) ) 。

因此，当模型RPCP-G和RPCP-O的目标函数值相等时，就找到了原问题的最优解。

4 算例测试

这一节进行数值实验来验证所建立的单层模型的有效性，模型采用优化求解器Gurobi 6.5进行求解，通过API由C#调用在VS2012平台上编程实现。实验在一台CPU为Intel Core Dual i3 2.53 GHz、内存为3 GB的笔记本上运行。

选取“离散选址问题标杆库”^①中的RPCP算例，其中客户集合和候选设施集合相同(即I=J)，在7000×7000的欧氏平面中随机选取，客户i和设施j之间的距离d_ij用欧氏距离表示，客户需求有2种情况：μ_i~U(1, 200), ∀ i∈I和μ_i=1, ∀ i∈I。每个算例的节点数量都是|I|=|J|=100，为减小规模，本文只选取前20个或者30个节点进行测试，这些节点同时作为客户和候选设施。

固定p=r=5，表1和表2分别统计了情形μ_i~U(1, 200)和μ_i=1的计算结果。对每个算例，列“总需求”表示了所考虑客户的需求之和，即总市场份额；用Gurobi求解器将模型RPCP-G计算到最优解，计算时间表示在列“时间/s”中，以s为单位，最优目标函数值表示在列“上界”中；然后固定领导者的决策，将模型RPCP-O求到最优，目标函数值表示在列“下界”中；上界相对于下界的间隙，表示在列“间隙/%”中。

从表1中可以看出，当规模较小(|I|=|J|=20)时，全部20个算例中，可以找到17个最优解；对其他3个算例，最优间隙也非常小。当规模较大(|I|=|J|=30)时，全部20个算例中，也可以找到8个最优解。由此可以推断出，本文所建立的单层模型对小规模算例非常有效，这是因为当规模较小时，跟随者可选的设施数量较少，贪婪策略能以较大概率逼近最优策略。比较2种情形下的计算时间，可以观察到，表2所列的μ_i=1情形下的计算时间，约为表1所列的μ_i~U(1, 200)情形下计算时间的2倍。这是因为当客户的需求都相等时，捕获客户的差异较小，搜索时间会更长。当|I|=|J|=30时，实际节点数为60，部分算例的计算时间已经超过1 h。对更大规模的算例，模型RPCP-G将不能用求解器直接求解，需要设计相应的列生成或者分支切割算法，这将是下一步研究课题。

5 结论

RPCP是一类经典的设施选址问题，常建模成为一个双层线性规划问题。启发式算法无法保证最优性，精确算法的时间复杂度太高。本文的一个主要贡献是用约束条件表示跟随者的响应，并集成到领导者的问题中。这样可以把双层规划问题简化成单层问题。本文假设跟随者采取贪婪策略，给定领导者的决策，求解模型RPCP-G和RPCP-O可以分别得到原问题的上界和下界。理论推导和数值实验证明了模型可以为领导者提供近似最优解，实验也表明计算负荷随问题规模显著增加。下一步研究重点是脱离MIP求解器、开发适用于此单层模型的有效算法。本文研究的问题保持了竞争选址领域的共性，未来可以考虑不同的顾客选择行为，不同的跟随者的设施数量。还可以考虑在其他应用场景下引入竞争的概念，如电子商务环境下配送中心的选址^[14]及不确定环境下的设施选址^[15]等。