周 迅，张惠珍

（上海理工大学管理学院，上海 200093）

0 引言

设施选址与车辆路径是供应链网络中的两个关键问题，对公司运营和盈利有着重大影响。因此，设施选址问题［1］（Facility Location Problem，FLP）与车辆路径问题［2］（Vehicle Routing Problem，VRP）引起了众多研究者关注。随着物流系统的日益复杂，选址和运输之间的耦合越来越高，储运之间关系日益密切。在解决许多实际问题时，应综合考虑两者关系进行决策，从而实现整体优化，由此形成具有容量的设施选址问题（Capacitated Location Routing Problem，CLRP）［3］，该问题是当前物流研究热点。近年来，越来越多的企业将主要精力投入到核心业务中以获得竞争优势，同时将企业物流活动外包给第三方物流供应商。通常服务车辆由第三方物流公司拥有，在完成客户服务之后，这些车辆返回到第三方物流公司，而不是企业配送中心或仓库，因此带来由第三方物流配送的开放式选址路径问题（Open Capacitated Location Routing Problem，OCLRP）［4］。CLRP 与相应OCLRP 区别如图1 所示。图中实线部分给出了OCLRP 问题解决方案，每辆车在为最后一个客户结束服务后无需回到仓库；图中实线和虚线整体构成了相应CL⁃RP 解决问题的方案，车辆在为最后一个客户结束服务后需返回仓库。从图1 可以看出，CLRP 与OCLRP 的主要区别为OCLRP 中的每条路径都是哈密顿路径，而CLRP 中的每条路径都是哈密顿回路。

Fig.1 The difference between CLRP and corresponding OCLRP图1 CLRP 与相应OCLRP 之间的区别

在实际配送中，OCLRP 会面临各种限制条件。比如，很多客户对交货期都有一定要求。配送服务最好提前在约定的时间窗口内完成，否则会因违反时间窗口约束而受到处罚。此外，在很多情况下，车辆不仅需要为客户进行送货服务，还有可能在客户所在位置进行取货，将货物运回收货点或第三方公司车辆中心。因此，本文将OCLRP 与同时取送货及软时间窗［5］相结合，提出了带软时间窗的同时取送货［6］开放式选址路径问题（Open Capacitated Loca⁃tion Routing Problem with Soft Time Windows and Simultane⁃ous Pickup-Delivery，OCLRP-STWSPD）。OCLRP-STWSPD问题集FLP 和VRP 两个NP-hard 问题为一体，并结合物流配送的实际情况增加了软时间窗和同时取送货限制，不仅复杂度更高，且更贴近现实。

当采用精确算法求解大规模组合优化NP-hard 问题时，通常难以得到最优解，且求解效率不高，而智能优化算法具有在较短的时间内对大规模优化问题提供近似最优解的优点。智能优化算法常被用来求解此类问题，如蚁群算法［7］、遗传算法［8］、模拟退火算法［9］、粒子群算法［10］等。烟花算法［11］（Fireworks Algorithm，FWA）是2010 年Tan &Zhu 受烟花爆炸产生火花这一自然现象启发而提出的一种新型智能优化算法，具有寻优效果好、稳定性强、效率高的优点。近年来，很多改进的FWA 算法已被提出，如反向烟花算法［12］、混沌烟花算法［13］、二进制烟花算法［14］等，同时FWA 也被广泛应用于云计算［15］、WEB 服务［16］、置换流水车间［17］等多个领域。鉴于OCLRP-STWSPD 问题的复杂性及FWA 良好的优化性能，本文对FWA 加以改进并用于求解OCLRP-STWSPD 问题。

本文根据OCLRP-STWSPD 问题特征，提出新的解表示方法，并对标准的FWA 加以改进，提出求解OCLRP-ST⁃WSPD 的离散烟花算法。实验证明该算法可有效解决目前第三方物流配送中同时取送货和配送时间约束问题，对节约企业配送成本、提高客户满意度具有重要意义，同时又可进一步拓展烟花算法应用领域。

1 带软时间窗同时取送货开放式选址路径问题

1.1 问题描述

OCLRP-STWSPD 是在传统OCLRP 问题上加入同时取送货和软时间窗限制，车辆需在客户能够接受的时间内进行配送或取货，否则将产生惩罚成本。OCLRP-STWSPD 依据客户点和候选仓库点空间位置与服务需求关系，确定开放仓库与配送路径，由开放仓库派出车辆为客户提供配送和取货服务，使物流系统总成本达到最小，并要求：①车辆对客户除了提供配送服务，还可能从客户处取走货物（客户货物量已知）；②在配送过程中每位客户有且仅能被1 辆车访问1 次，且客户配送需求或取货需求需得到满足，车辆在配送过程中任意时刻的装载量不能超过车辆最大装载量；③分配给任一开放仓库的客户总需求不能大于该仓库容量；④每辆车为最后1 个客户完成服务后不返回仓库，而是返回收货点或第三方公司车辆中心；⑤在配送过程中，每个客户有自己接受服务的时间窗限制，车辆可在客户预先给定的服务时间之外为其服务，即车辆可早到或晚到，但违反时间窗的服务会产生一定的惩罚成本；⑥若一位客户同时具有配送和取货需求，则服务顺序应该为先卸货后取货。

1.2 数学模型

模型主要参数设置包括：Vp为第三方物流公司的停车场集合；Vd为候选仓库集合，Vd={1,23,…,m}，其中m为候选仓库数量；Vc为客户点集合，Vc={m+1,m+2,…,m+n}，其中n为客户数量；V为停车场、候选仓库和客户点集合，V=Vp∪Vd∪Vc；K为配送车辆集合；Oi为仓库i开放成本；di为客户i配送量；bi为客户i取货量；Qi为仓库i容量；P为每一辆车最大装载量，假设所有配送车辆相同；cij为节点i到节点j的距离；e为单位距离的运输成本；w为每辆车固定使用成本；ETi为客户i接受服务的最早开始时间；LTi为客户i接受服务的最晚开始时间；Ti为配送车辆到达节点i的时间；si为客户i的服务时间；tij为从节点i到节点j的时间；pe为提前到达的单位时间惩罚成本；pl为延迟服务的单位时间惩罚成本；aik为车辆k在离开节点i时的装载量。

决策变量设置如下：

ujk=用于消除第k条配送路线中子回路的辅助变量

在上述参数和决策变量设置基础上，OCLRP-STWSPD数学模型可构建为：

其中，目标函数（1）表示最小化总成本，包括仓库开放成本、车辆的固定启用成本、行驶成本和惩罚成本；约束（2）确保每个客户只由1 辆车提供1 次服务；约束条件（3）表示车辆在任何时刻的载重量都不超过其最大装载量；约束条件（4）表示车辆在离开仓库时的载重量，即配送路线上所有客户配送量之和；式（5）为车辆在离开各个客户时的装载量的等式约束；约束条件（6）表示分配给仓库的客户总需求不能超过仓库的总容量；约束（7）保证每条配送路线连续性；约束（8）表明任意两个仓库之间没有连接；约束（9）确保配送车辆为最后1 位客户完成服务后不返回仓库；约束（10）说明每个客户都必须被一个开放的仓库提供服务；约束（11）表示两个决策变量之间的关系；约束（12）表示车辆在行驶过程中的时间窗约束；约束（13）用于消除子回路；约束（14）-（16）表示3 组不同的0-1 决策变量；约束（17）为消除子回路所用的辅助变量。

2 OCLRP-STWSPD 离散烟花算法求解

2.1 传统烟花算法

烟花算法［11］是一种受烟花爆炸现象启发，通过模拟烟花爆炸过程实现的群体智能优化算法。为了保持种群多样性，烟花算法设计了两个爆炸搜索过程。烟花算法将每一个烟花和烟花爆炸产生的每一个火花视为一个可行解，每一次烟花在不同爆炸半径内产生不同数量的火花。通过调整爆炸半径和火花数量，实现全局搜索与局部搜索平衡。该算法由爆炸算子、变异算子和选择策略3 部分组成。

2.1.1 爆炸算子

假设研究问题是最小化目标函数f(x)，F为解的可行域。首先在可行域F上初始化生成N个烟花Xi(i=1,2,3,…,N)，并评估初始烟花函数f(Xi)。烟花Xi通过爆炸算子在半径Ai上产生Si个火花，每个烟花爆炸半径Ai与产生的爆炸火花数Si的计算公式分别为：

其中，ymin为当前烟花种群对应的目标函数最小值，ymax为目标函数最大值；ε 是一个很小的常数，以防止分母为0。A0和M0为常数，分别用以调整每个烟花爆炸产生的半径大小和爆炸火花总数量。为了区分烟花优劣，函数值越小的烟花在小半径内产生的火花越多，相应局部搜索能力越强［15］；相反，函数值较高的烟花在较大范围内产生相对较少的火花，具有一定的全局搜索能力。

为控制火花数量过多或过少，引入式（18）所示的机制控制每个烟花产生的火花数量Si。

其中，a和b为两个常数，round为四舍五入函数。

2.1.2 变异算子

在变异算子中，通常采用高斯变异产生变异火花加强烟花种群多样性。高斯变异火花产生过程为：首先，在烟花种群中随机选择1 个烟花Xi;然后，随机选择Xi中的若干分量，并按照公式（19）进行变异操作［10］，重复M1次，即可产生M1个变异火花。

其中，r为服从均值和标准差均为1的高斯分布随机数。

2.1.3 选择策略

为了将当前烟花种群中的优秀信息传递到下一代种群中，烟花算法会在候选种群H中选择N个个体作为下一代烟花。其中，候选种群H由烟花、爆炸火花和高斯变异火花组成［10］。选择过程为：首先，采用精英策略保留最优秀的烟花，即函数值最小解；其次，使用轮盘赌策略从候选集合中选择N-1 个个体。其中每个个体被选中的概率为：

2.2 离散烟花算法

传统烟花算法最初用于连续域函数优化，算法步骤中的爆炸算子和变异算子仅适用于连续域空间搜索。由于离散优化问题与连续优化问题的区别，传统烟花算法不能直接应用于OCLRP-STWSPD 和其他离散优化问题求解。因此，本文依据OCLRP-STWSPD 问题特征和烟花算法搜索机制，对爆炸算子和变异算子重新定义，并针对OCLRPSTWSPD 提出新的解表示方法，在算法搜索过程中加入自适应策略，进而提出求解OCLRP-STWSPD 问题的离散烟花算法（Discrete Fireworks Algorithm，DFWA）。

2.2.1 初始解构建

（1）编码方式。假设有m个候选仓库和n个客户，则初始解由Vd={1,2,3,…,m}集合内的元素和Vc={m+1,m+2,…,m+n}中的元素以及一些辅助数字0 组成。其中，辅助数字0 用来分割路径，辅助数字0 的数量应该等于或小于表示大于或等于x的最小整数）。分配给每个仓库的客户位于它与下一个不同的仓库之间，如果在一个仓库之后没有客户，则表示该候选仓库没有开放。

图2 给出1 个10 位客户和3 个备选仓库的例子。如图2 所示，仓库3 未开放，仓库1、2 开放。其中，从仓库1 出发的两条配送路线为1 →9 →10→13 和1 →12→4；从仓库2 出发的两条配送路线为2→8 →11→5 和2 →6→7。

Fig.2 Solution representation图2 解的表示

（2）种群初始化。为了产生好的初始种群，本文采用基于贪婪思想的策略构建初始解。

步骤1：设Per为Vd中m个仓库的随机排列。将Per中的第i个元素表示为Peri，令i=1。

步骤2：设Cus为未被分配客户的集合，C(Peri)表示未被分配客户中需求不超过仓库Peri剩余能力的客户集合。首先，对未分配的客户进行分配，并从Peri仓库开始一条新的路线，将C(Peri)中距离Peri仓库最近的客户分配给Peri仓库；其次，将C(Peri)中距离上一个已分配客户最近的客户分配到当前路线中，如果该客户需求量大于车辆剩余容量，则将该客户从该路线中删除，并继续从C(Peri)中寻找距离上一个客户最近且需求量满足车辆剩余容量的客户。若不存在客户满足条件，则从Peri开始一条新的路线，并重复上述操作，直到Peri仓库剩余容量无法服务任何新客户为止。

步骤3：如果Cus不为空，则令i=i+1，并重复步骤2。否则，将使用2.1.1 中的编码方式对当前序列进行编码作为一个初始解。

步骤4：如果初始解的个数小于预先设定的烟花种群大小N，返回步骤1；否则，输出初始种群。

2.2.2 爆炸算子

针对OCLRP-STWSPD 离散特征，本文采用交换操作实现爆炸算子［5］，并将爆炸半径Ai定义为交换次数。具体操作为对烟花Xi，随机选取Xi上的Ai对客户，并对选中位置上的元素进行交换，得到新的序列。

图3 为爆炸算子具体操作。如图3 所示，当烟花Xi爆炸半径Ai为2 时，对烟花序列Xi选择两对客户进行交换操作，选择5 和7、10 和12 进行交换操作形成新序列。

Fig.3 Illustration of generating a spark with Ai=2图3 为2 时的爆炸算子操作

2.2.3 变异算子

（1）变异算子实现。DFWA 变异算子是通过插入和逆转两种搜索机制［5］实现的。变异火花通过这两种操作产生，两种操作概率均为0.5。

生成随机数rand，若rand小于0.5，则使用插入操作：随机选取Xi的两个位置a、b，将b位置上的元素插入到a位置元素的前面；若rand大于等于0.5，则使用逆转操作：随机选取Xi的两个位置a、b，逆转a、b之间的元素排列。

图4 和图5 分别给出了变异算子的两种实现方法。图4 为插入操作，将元素9 插入到元素7 前面产生一个新的序列；图5 为逆转操作，将元素5 和元素13 之间的所有元素逆转其排列顺序生成1 个新的序列。

Fig.4 Illustration of the insertion move图4 插入操作

Fig.5 Illustration of the inverse move图5 逆转操作

（2）自适应策略。基本思想为：如果当前代搜索到更优解，则认为当前烟花种群具有较强的寻优能力；否则，说明可能已经陷入局部最优。基于该思想，算法利用变异算子中的参数q改变种群寻优能力。若当前代搜索到更优解，则减小q以加强局部搜索；若当前最优解没有得到改进，则增大q以提高接受较差解的概率，以便更快跳出局部最优解。

在离散优化过程中，往往很难在1 次迭代中得到较好的解。因此，当烟花种群落入局部极值时，往往需多次迭代才能跳出局部收敛。鉴于此，本文设定如果连续10 代不改进最优值，则增大q值。公式（22）-（23）分别用来增大和减小q值。其中，r和h均为常数，分别为增长因子和衰退因子。

2.3 算法步骤

基于上述操作，求解OCLRP-STWSPD 的DFWA 具体步骤为（见图6）：

步骤1：初始化。首先，初始化烟花种群大小N、爆炸因子A0、爆炸火花数M0、变异火花数M1、控制参数q、增长因子r、衰退因子h、最大迭代次数T；然后使用贪婪策略初始化烟花种群。

步骤2：对烟花种群进行爆炸算子和变异算子操作，产生新的爆炸火花和变异火花。

步骤3：从候选种群（烟花，爆炸火花，变异火花）中选择N个烟花作为下一代种群。

步骤4：若当前代搜索过程产生更优解，则减小q；若当前代解连续10 代没有得到改进，则增大q值。

步骤5：重复步骤2、3、4、5，直到达到最佳迭代次数T，输出最优解。

Fig.6 The flow of DFWA algorithm图6 DFWA 算法流程

2.4 时间复杂度分析

假设初始烟花种群大小为N，n名客户和m个仓库，爆炸算子产生M0个爆炸火花，变异算子产生M1个变异火花。DFWA 在构建初始解阶段需生成N个初始解，且对每个仓库进行客户分配排序，其时间复杂度为O(N×m×n)；在爆炸算子阶段，N个烟花产生M0个爆炸火花，其时间复杂度为O(N×M0)；在变异算子阶段产生了个M1变异火花，其时间复杂度为O(M1)。在选择阶段需从候选种群中保留最优个体，涉及到所有候选种群个体排序，其时间复杂度为O((N+M0+M1)log2(N+M0+M1))。综合各个步骤，DF⁃WA 时间复杂度为O(N×m×n)+O(N×M0) +O(M1)+O((N+M0+M1)log2(N+M0+M1))=O(N×m×n)+O(N×M0) +O((N+M0+M1)log2(N+M0+M1))。

3 实验分析

3.1 参数设置与灵敏度分析

参数设置可能影响DFWA 性能。因此，本文进行灵敏度分析以确定最佳参数设置。实验在Windows10 系统下的MATLAB R2016a 环境进行。测试参数值设置如下：

每种参数组合独立运行10 次，并使用平均目标函数值与平均运行时间为评价指标。当研究其中一个参数与算法性能之间关系时，其他参数保持不变。

3.1.1 烟花种群大小N

如图7 所示，实线与虚线曲线分别表示在不同参数N下的目标函数值和计算时间。随着N的增大，解空间也随之变大，因此得到优秀解的概率也会增加。从图1 可以看出，N对DFWA 的影响较大，当N增大时，目标值会不断减小，计算时间也会相应提高。当N提高到1.25（n+m）后，解的质量提升速度明显下降，而计算时间依然会增加，因此综合解质量和计算时间，本文选取N=1.25()n+m作为烟花种群大小。

Fig.7 Sensitivity analysis of parameter N图7 参数N 的灵敏度分析

3.1.2 爆炸因子A0

在DFWA 中，A0用于调整爆炸半径，以控制产生每个火花所需交换次数。当A0较小时，交换次数较小，交换产生的火花不能得到充分改善；但是当A0较大时，也可能陷入无意义迭代，增加计算时间。如图8 所示，当A0=2N时，计算结果较好，相应计算时间也相对合理。

Fig.8 Sensitivity analysis of parameter A0图8 参数A0的灵敏度分析

3.1.3 爆炸火花数M0

爆炸火花数M0越少，对邻域搜索力度越小，得到最优解的可能性越小。但是，较多的爆炸火花可能与穷举法非常相似，会大幅增加计算时间。如图9 所示，当M0=(n+m)时，DFWA 寻优效果和计算时间较为合理。因此，本文将爆炸火花数量设置为M0=(n+m)。

Fig.9 Sensitivity analysis of parameter M0图9 参数M0的灵敏度分析

3.1.4 变异火花数M1

变异火花数M1和爆炸火花数M0效果相似。变异火花可提高种群多样性，可有效防止算法陷入局部最优解。M1值越小，其搜索范围越小；M1值越大，其搜索范围越大，种群多样性也越高，但M1取值较大时，相应计算时间也越高。如图10 所示，当M1=(n+m)时，可较好平衡解质量和计算时间。

Fig.10 Sensitivity analysis of parameter M1图10 参数M1的灵敏度分析

3.1.5 控制参数q

在DFWA 中，参数q用来调节接受较差解的概率。当q较小时，接受差解的概率也较小，可加速收敛速度，同时也可降低种群多样性，加大陷入局部最优解的可能性；当q较大时，接受差解的概率较大，种群多样性提高的同时也降低了收敛速度。图11 给出了在不同q值大小下DFWA 的表现，可以看到当q值较小或较大时，DFWA 寻优效果并不理想。当q=1 时，DFWA 可较好平衡收敛速度与寻优效果。

Fig.11 Sensitivity analysis of parameter q图11 参数q 的灵敏度分析

3.2 算例分析与结果比较

目前尚未有标准的实验数据，为了验证本文OCLRPSTWSPD 模型与DFWA 算法有效性与可行性，对文献［19］和［20］给出的CLRP 算例加以改造进行实验。算例中各节点坐标、仓库费用和容量、车辆最大容量、客户配送量以及各项固定成本费用均采用文献［19］、［20］所给算例原始数据，客户取货量采用均匀分布的方式生成，均匀分布上下限为对应算例中配送量最小值和最大值，而客户时间窗上下限ETi和LTi同样采用均匀分布的方式生成，均匀分布区间分别为［0，10］和［ETi+30,ETi+60]。DFWA 参数设置为参数设置与灵敏度分析中找到的最优参数，此外，增长因子r=1.1、衰退因子h=0.98、最大迭代次数T=500。

第一组数据客户数为20 或50，候选仓库数均为5，单位运输成本为100，每辆车固定租赁成本为1 000，单位惩罚成本Pe和Pl分别为50 和80，车辆行驶速度为3；第二组数据客户数量范围从8 到50，候选仓库数为2 或5，单位运输成本为1，单位惩罚成本Pe和Pl分别设置为5 和8，车辆行驶速度为5。表1 和表2 分别给出了DFWA 对第一组算例中每个算例独立求解10 次的结果，包括最小成本值、最小成本值出现次数、平均成本值、平均运行时间以评估算法性能。如表1 所示，算法求解20-5 的算例时，均可达到最小成本值，平均运行时间最长达10.33s。对于50-5 的算例，算法可在10 次运行中至少2 次达到最小成本值，平均运行时间最长达32.93s。表2 给出了DFWA 对第二组算例中每个算例独立求解10 次的结果。对于第二组算例，算法在10 次运行中至少3 次可达到最小成本值，且计算时间不高于31.63s。

为了更进一步验证DFWA 性能，本文选择文献［18］中的OCLRP 算例进行直接计算并与文献［18］中的算法进行比较，使用DFWA 运行每组算例10 次后取其中最优值，其比较结果如表3 所示。

表3 提供了CPLEX 优化器、SA［18］、DFWA 在18 组算例上的表现。从求解结果方面看，CPLEX 求解器可以在16 组算例上达到最优值，而SA（模拟退火算法）和DFWA 均可以在17 组算例上取得最优值；从求解效率上来看，CPLEX 优化器在小规模且客户数不多于22 的算例上效率更高，甚至在一些算例上可以在1s 以内得到优化结果。但在大规模的算例上，CPLEX 求解效率较低，某些大规模算例在4 小时内仍未求得最优解。而SA 和DFWA 在大规模算例上的求解效率明显更高，从时间上来看，DFWA 在13 组算例上的计算时间少于SA，且计算时间大幅度下降。综合来看，DF⁃WA 在求解此类问题时，能够在较短时间内求解出令人满意的解，且具有寻优效果好、求解效率高的优点。

Table 1 The OCLRP-STWSPD results of the first set of data表1 第一组数据的OCLRP-STWSPD 结果

Table 2 The OCLRP-STWSPD results of the second set of data表2 第二组数据的OCLRP-STWSPD 结果

3.3 OCLRP-STWSPD 与CLRP-STWSPD 比较

为证明OCLRP-STWSPD 模型在现实生活中的适用性，本部分对上述数据集在相同时间窗和相同参数设置下，进行带软时间窗的同时取送货有容量选址路径问题（Capacitated Location Routing Problem with Soft Time Win⁃dows and Simultaneous Pickup-Delivery，CLRP-STWSPD）实验，并将CLRP-STWSPD 运行结果与上述OCLRP-STWSPD结果进行对比。

表4 给出了运用DFWA 对OCLRP-STWSPD 和OCL⁃RP-STWSPD 在所给算例上的求解结果、差值（差值=CL⁃RP-STWSPD－ OCLRP-STWSPD）和Gap值（Gap=。

Table 3 Comparison of DFWA and other algorithms表3 DFWA 与其他算法的比较结果

Table 4 Comparison of solution results of OCLRP-STWSPD and CLRP-STWSPD表4 OCLRP-STWSPD 与CLRP-STWSPD 求解结果对比

从表4 可知，在相同时间窗下，OCLRP-STWSPD 成本值远小于CLRPSTW-STWSPD，最大Gap达到30.20%，可节省高达15 270 单位的成本。这是由于OCLRP-STWSPD 物流系统外包给第三方公司，为所有客户完成服务后的车辆无需返回仓库以节省开支。这些结果进一步证明了OCL⁃RP-STWSPD 在现实生活中的适用性，同时也表明如果公司将配送活动外包给第三方，可节省大量配送成本。

4 结语

本文提出了一个应用背景更符合实际的带软时间窗同时取送货开放式选址路径问题，并以最小化仓库开放成本、配送成本、固定车辆车本、惩罚成本之和为目标建立数学模型。针对传统烟花算法求解离散问题的局限性，基于交叉逆转等操作重新定义爆炸算子与变异算子以维持种群多样性，并加入自适应策略避免过早陷入局部最优解，提高算法搜索能力。同时，在算法中加入新的解表示方式分割求解路径。最后通过算例求解与结果比较，证明了DFWA 在求解OCLRP-STWSPD 问题上有效性和可行性，可很好地解决目前第三方物流配送中面临的现实问题，对企业节约配送成本与提高客户满意度具有重要意义。下一步将集中从3 个方面开展更深入研究：①考虑更多的约束条件，比如考虑交通路况的OCLRP、带模糊需求的OCLRP、多目标的OCLRP；②设计其他启发式算法解决OCLRP 问题；③采用DFWA 解决FLP 与VRP 其他变体。