车辆路径问题（vehicle routing problem，VRP）由Dantzig等^[1]首次提出，是物流交通领域的核心问题，也是一类经典的组合优化问题^[2]. 包含产品成本、多车场、多车型在内的多约束车辆路径问题（rich vehicle routing problem，RVRP）是传统VRP问题的延伸，对于多中心分布式企业来说，在进行车辆调度时，不仅要考虑所服务的客户需求，还要综合考虑多个车场的位置、产品成本以及不同类型车辆的属性，在求解难度上较普通VRP问题大大增加，目前仍缺乏成熟且有效的求解算法.

由于实际问题的研究角度不同，不同学者对传统VRP问题添加了各种特定约束，衍生出许多新问题. 针对多车场多车型车辆路径问题^[3]（multi-depot heterogeneous vehicle routing problem，MDHVRP），Salhi等^[4]研究附带拣货作业的多车型车辆路径问题，并提出基于集合划分的启发式算法对其进行求解；Sundar等^[5]针对带有燃料限制的多车场多车型问题，设计分支切割算法，但该算法只适用于小规模情形；Luo等^[6]在混合蛙跳算法中引入极值优化算法对多车场带时间窗的车辆路径问题进行研究；Dayarian等^[7]以牛奶采集系统为例，提出分支定价法对多属性车辆路径问题进行研究；Calvet等^[8]考虑不同仓库配送对需求变化的潜在影响，设计统计学习技术与元启发式框架相结合的混合方法；de Oliveira等^[9]将MDVRP问题分解成多个经典VRP子问题，提出具有可变长度基因和搜索算子的并行进化算法，通过算例验证所提算法的有效性；Calvet等^[10]研究车队数量有限即仓库容量有限的多车场问题，提出结合蒙特卡罗模拟和元启发式算法的模拟启发式框架；杨烨^[11]对带时间窗的多车型满载车辆调度问题进行研究，并设计了“扫描-插入-遗传”启发式算法，但是其所得解很依赖扫描情况，往往无法得到全局最优解；叶勇等^[12]考虑配送中心的动态启用性，设计了改进狼群算法；马宇红等^[13]引入司机工资，包括基本工资和加班费，以最小配送费为目标，对多车场多车型车辆调度问题进行研究，但其在算法验证过程中，案例规模较小，并且将多车型简化为了单车型，降低了求解难度，使得结果不具备充足的信服力；刘丽姣^[14]将碳排放引入多车场多车型车辆路径问题中，提出单循环结构综合学习细菌觅食优化算法，通过单峰函数和多峰函数对其进行有效性验证，并证明多车场多车型较单车场单车型的优势；何东东等^[15]建立带时间窗的多车型绿色车辆路径模型，并设计了改进的禁忌搜索算法；Baradaran等^[16]以实际车辆配送系统为基础，研究具有多个附带优先级的硬时间窗的多车型车辆路径问题，提出3种多目标模型，并利用二进制人工蜂群算法进行求解验证；杨翔等^[17]针对模糊需求下的开放式多车场车辆路径问题，采用两阶段禁忌搜索算法进行研究.

多年来，学者们对MDHVRP问题进行了各式各样拓展性的研究，然而较少有学者研究不同产品成本对多车场车辆调度的影响. 2016年，鲁建厦等^[18]首次将产品成本考虑到多车场问题中去，不同地区用人成本、原材料成本不同，往往导致其产品成本不同，产品成本低的配送中心理应服务更多的客户，这也是很多大型跨国企业将加工厂建在发展中国家的原因，主要是为了利用发展中国家较为低廉的劳动力、原材料，从而在激烈的国际竞争中获得优势，赢取更大的利润，这也证明了研究产品成本对调度影响的重要性与必要性. 然而，鲁建厦等^[18]的初步研究并未证明考虑产品成本对优化调度的必要性，为此，本研究在其基础之上，进一步增加模型复杂度，以最小化总成本为目标，提出考虑产品成本、多车场、多车型在内的多约束车辆路径问题（rich vehicle routing problem，RVRP）模型. 设计改进的混合蛙跳算法（improved hybrid shuffled frog leaping algorithm，IHSFLA）进行求解，并通过实例验证模型和算法的有效性.

RVRP问题可以描述为：某调度中心在其服务区内拥有多个分布在不同地点的配送中心，每个配送中心拥有不同类型的车辆，其产品的成本也互不相同；该调度中心要为其所在服务区内的客户进行配送服务，要求在满足客户的需求条件下，安排合理的车辆配送方案使产品总成本及车辆配送成本之和最小.

该问题基于如下假设：1）调度中心指挥多个配送中心，每个配送中心拥有足够多数量的不同型号的车辆；2）每辆车最多只被调度一次，车辆从配送中心出发，完成任务后返回到原配送中心；3）车辆车型不同，承载量不同；4）车辆车型不同，车辆的固定使用成本和行驶单位距离成本不同；5）每一辆车至少可以服务一个客户，即每个客户的需求量不会超过其车辆的最大承载量；6）每一个客户只能由一辆车服务一次；7）每个配送中心产品的成本互不相同.

RVRP问题可以用加权图G=（V，E）来表示，节点集V={C $ \cup D $}，C={1，2，···，N}为包含N个客户的客户集，D={N+1，N+2，···，N+H}为包含m个车场共H辆车的车辆集，E={（i，j）|i， $ j\in V $， $i \ne j $}为节点i和j之间边的集合. 目标是在满足约束的条件下，确定每辆车的顾客服务顺序，使得产品成本及车辆配送成本之和最小.

决策变量为

(1) $ {x}_{ij}^{h}=\left\{\begin{array}{l}1,\;\;{\text{从节点}}i{\text{到节点}}j{\text{由车辆}}h{\text{配送}},i\ne j;\\ 0,\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\text{其他}}.\end{array}\right. $

根据以上条件分析，建立RVRP数学模型如下：

(2) $\begin{split} {{Z}} =& {\rm{min}}\;\sum\limits_{h = 1}^H {\sum\limits_{i = 1}^{N + H} {\sum\limits_{j = 1}^N {{r_j}} } } {\alpha _h}x_{ij}^h + \sum\limits_{h = 1}^H {\sum\limits_{i = N + 1}^{N + H} {\sum\limits_{j = 1}^{N + H} {{C_h}} } } x_{ij}^h + \\ &\sum\limits_{h = 1}^H {\sum\limits_{i = 1}^{N + H} {\sum\limits_{j = 1}^{N + H} {{\beta _h}} } } {d_{ij}}x_{ij}^h. \end{split}$

约束条件如下：

(3) $ \sum\limits_{i = 1}^{N + H} {\sum\limits_{j = 1}^N {{r_j}} } x_{ij}^h \leqslant {Q_h};\;{{h}} = {\rm{1}},{\rm{2}}, \cdots ,{{H}}, $

(4) $ \sum\limits_{h = 1}^H {\sum\limits_{i = N + 1}^{N + H} {\sum\limits_{j = 1}^{N + H} {x_{ij}^h} } } \leqslant H, $

(5) $ \sum\limits_{h = 1}^H {\sum\limits_{i = 1}^{N + H} {x_{ij}^h} } = 1;\;{{j}} = {\rm{1}},{\rm{2}}, \cdots ,{{N}}, $

(6) $ \sum\limits_{h = 1}^H {\sum\limits_{j = 1}^{N + H} {x_{ij}^h} } = 1;\;{{i}} = {\rm{1}},{\rm{2}}, \cdots ,{{N}}, $

(7) $ \sum\limits_{i = 1}^{N + H} {x_{i(N + h)}^h} = \sum\limits_{j = 1}^{N + H} {x_{(N + h)j}^h} \leqslant1;\;{{h}} = {\rm{1}},{\rm{2}}, \cdots ,{{H}}, $

(8) $ \sum\limits_{i = N + 1}^{N + H} {\sum\limits_{j = N + 1}^{N + H} {x_{ij}^h} } = 0;\;{{h}} = {\rm{1}},{\rm{2}}, \cdots ,{{H}}. $

式中： $ {r}_{j} $为客户j的需求量， $ {d}_{ij} $为节点i到节点j之间的距离， $ \alpha $_h为车辆h所代表的配送中心产品的单位生产成本，C_h为车辆h的固定使用成本， $ {\beta }_{h} $为车辆h的单位距离行驶成本， $ {Q}_{h} $为车辆h的最大承载量. 式（2）表示总成本，第1部分为产品成本，第2部分为配送车辆的固定成本，第3部分为配送车辆的距离成本；式（3）表示每辆车的单次配送总量不能超过其最大承载量；式（4）表示调度中心使用车辆数不超过其各配送中心可用车辆总数；式（5）、（6）表示一个客户只能由一辆车服务一次；式（7）表示车辆从配送中心出发完成任务后返回原配送中心；式（8）表示车辆不能从一个配送中心行驶到另一个配送中心.

多约束车辆路径问题属于典型的NP难题^[19]，在有效时间内，很难利用精确算法进行计算，因此通常采用启发式算法对其进行规划求解. 混合蛙跳算法（shuffled frog leaping algorithm，SFLA）由Eusuff等^[20]提出，模仿一群青蛙在觅食时的信息交互行为. 该算法结合模因算法和粒子群算法的优点，概念简单、调整参数少，具有较强的全局搜索能力，已成功在车辆路径规划^[21]、柔性车间调度^[22]领域内运用. 在混合蛙跳算法中，种群由一群青蛙组成，每一只青蛙代表一种可能的解决方案，将这些青蛙按照某种规则分到不同的族群内，青蛙在各自的族群内部进行交流演化，在一定时间后，所有族群的青蛙会聚集在一起，进行融合，然后再行划分、进化，循环以往，直至满足收敛条件.

为了加快算法的收敛速度，避免陷入局部最优解，本研究提出基于聚类的改进混合蛙跳算法来求解RVRP问题. 采用四标准聚类算法^[23]对所有客户进行基于配送中心的聚类；通过客户邻近矩阵产生初始解的方式构造初始种群，将初始种群分为不同族群，各族群内部进行交流、进化；将所有族群进行混合、重排，再次进行分组，循环以往，直至满足收敛条件.

算法通用框架如下所示.

1 输入初始参数

2 对客户进行四标准聚类分析

3 根据聚类结果以及邻近矩阵构造初始解

4 对初始解计算适应度进行升序排列，并按规则分 到N_f个族群内

5 记录全局最优解P_g

6 while 不满足全局收敛条件

7 for i=1∶m

8 While 不满足局部交流次数

9 在族群内部产生子群

10 记录子群内局部最优解P_b，最劣解P_w

11 利用P_b更新最劣解得到 $\scriptsize{P}_{{\rm{w}}}^{'}$

12 if f（ $\scriptsize{P}_{{\rm{w}}}^{'}$）>f（P_w）

13 P_w= $\scriptsize {P}_{{\rm{w}}}^{'}$

14 else

15 利用P_g更新最劣解得到 $\scriptsize {P}_{{\rm{w}}}^{'}$

16 if f（ $\scriptsize {P}_{{\rm{w}}}^{'}$）>f（P_w）

17 P_w= $\scriptsize{P}_{{\rm{w}}}^{'}$

18 else

19 随机产生一个解替代P_w

20 end if

21 end if

22 end while

23 记录族群内局部最优解P_b

24 while 不满足邻域搜索次数时

25 对P_b进行邻域搜索产生 $\scriptsize{P}_{{\rm{b}}}^{'}$

26 if f（ $\scriptsize {P}_{{\rm{b}}}^{'}$）>f（P_b）

27 P_b= $\scriptsize{P}_{{\rm{b}}}^{'}$

28 else

29 以概率选择是否接受新解

30 end if

31 end while

32 end for

33 对m个族群进行混洗融合，然后进行升序分配

34 end while

35 输出最优解P_g，f（P_g）

SFLA的初始模型是用来解决连续问题的，而车辆路径问题属于多维离散性问题，因此必须改变其编码方式使其可以解决离散性问题.

传统的VRP问题多采用一段式的自然数编码方式，即在一段自然数编码中，包含了所有车辆的行驶路线，如由1个配送中心6个客户组成的配送网络，其可能存在的编码方式为0132040650，其中0代表配送中心，数字1~6代表客户编号，对其进行解码为0-1-3-2-0，0-4-0，0-6-5-0，即6个客户分3辆车进行配送服务，然而传统的编码方式在解决多车场多车型问题时，普遍存在表达复杂、难以进行交叉、容易产生非法解等问题，因此，本研究采用多染色体编码方式^[24]. 多染色体编码方式就是在进行编码表达时，采用多条染色体的形式，对应到车辆路径问题中来，就是每一辆车都对应一条染色体编码，所有车辆的染色体编码组成一个个体，这样就避免了传统编码方式中难以表达车场、车型的问题，同时，因为每辆车都有自己的编码，所以在不同个体或者不同车辆进行信息交流时，更加方便，且不易产生非法解.

对于一个包含了m个车场、共H辆车、N个客户的RVRP问题，用编号1~N表示客户，N+1~N+H分别表示来自不同车场的共H辆车，分别对每辆车进行编码，每辆车服务客户的顺序称为该车的基因链. 如图1所示，例如2个车场包含2种类型的车共4辆，要为服务区内的9个客户进行服务，分别用数字1~9表示客户，10~13表示来自不同车场的车辆，10~11表示车场1中2种型号的车辆，12~13表示车场2中2种型号的车辆，分别对4辆车进行编码，则其可能存在的编码方式为[10，1，6，4；11，2，7；12，8，5；13，3，9]，分别表示车辆10从车场1出发，依次服务客户1、6、4，最后回到车场1；车辆11从车场1出发，依次服务客户2、7，最后回到车场1；车辆12从车场2出发，依次服务客户8、5，最后回到车场2；车辆13从车场2出发，依次服务客户3、9，最后回到车场2.

初始解的优劣对算法至关重要，为了提高初始解的优异性，加快算法搜索速度，采取先聚类后根据邻近矩阵生成初始解的方法. 在生成初始解之前，首先对客户进行聚类分析，将客户分给不同的车场，聚类分析可以有效降低计算的规模与复杂度，大幅度提高算法的收敛速度.

在文献[13]提出的三标准聚类算法的基础上，添加距离中位值的衡量标准，采用四标准聚类算法对所有客户进行聚类，并分配给对应的车场. 四标准指的是中位距离、平均距离、平均距离方差、最小距离，可以更好地反映客户所在地区的情况. 4种计算方式如图2所示。该方法比常用的K-means聚类方法更有效^[23]，分配结果更加合理，具体流程如下所示.

1 初始化， $ {C}_{{\rm{all}}}={C} $//C为所有客户集

2 while not isempty $ \left({C}_{{\rm{all}}}\right) $

3 每个车场及其所要服务的客户为一个聚类

4 计算 $ {C}_{{\rm{all}}} $[1]与各个聚类中所有点的中位距离

5 If dis_med $ ({C}_{{\rm{all}}} $[1]，2 $ )- $ dis_med $ ({C}_{{\rm{all}}} $[1]，1 $ ) \geqslant$ $ 0.1\times $dis_med $ ({C}_{{\rm{all}}} $[1]，2 $ ) $ // dis_med $ ({C}_{{\rm{all}}} $[1]，n $ ) $为客户 $ {C}_{{\rm{all}}} $[1]到第n近聚类的中位距离

6 把客户 $ {C}_{{\rm{all}}} $[1]分给距离它平均距离最短的聚类

7 ${C}_{{\rm{all}}}={C}_{{\rm{all}}}[2:{\rm{end}}]$

8 else

9 计算 $ {C}_{{\rm{all}}} $[1]与各个聚类的平均距离

10 if dis_avg $ ({C}_{{\rm{all}}} $[1]，2 $ )- $ dis_avg $ ({C}_{{\rm{all}}} $[1]，1 $ ) \geqslant$ $ 0.1\times $dis_avg $ ({C}_{{\rm{all}}} $[1]，2 $ ) $ // dis_avg $ ({C}_{{\rm{all}}} $[1]，n $ ) $为客户 $ {C}_{{\rm{all}}} $[1]到第n近聚类的平均距离

11 把客户 $ {C}_{{\rm{all}}} $[1]分给距离它平均距离最短的聚类

12 ${C}_{{\rm{all}}}={C}_{{\rm{all}}}[2:{\rm{end}}]$

13 else

14 计算 $ {C}_{{\rm{all}}} $[1]到各个聚类的平均距离方差

15 if dis_var $ ({C}_{{\rm{all}}}\left[1\right],\;{{m}})\leqslant $0.4// dis_var $ ({C}_{{\rm{all}}}\left[1\right],\;{{m}}) $为客户 $ {C}_{{\rm{all}}}\left[1\right] $与第m个聚类的平均距离方差

16 把客户 $ {C}_{{\rm{all}}} $[1]分给第m个聚类

17 ${C}_{{\rm{all}}}={C}_{{\rm{all}}}[2:{\rm{end}}]$

18 else

19 计算 $ {C}_{{\rm{all}}}\left[1\right] $与每个聚类中最近客户的距离

20 将 $ {C}_{{\rm{all}}}\left[1\right] $分给与聚类中最近客户距离最小的聚类

21 end if

22 end if

23 end if

24 end while

邻近矩阵用来描述2个点的邻近关系^[25]，通过两点之间的距离构造邻近矩阵，表达式为

(9) $ \lambda \left(i\leftarrow j\right)={\left({\lambda }_{ij}\right)}_{N\times (N+m)};\;i\in C,j\in C\cup M. $

式中：M为包含m个车场的车场集，M={1，2，···，m}； $ {\lambda }_{{{i}}j} $为客户i和客户（或车场）j的邻近程度， $ {\lambda }_{{{i}}j}\in ${0，1，2，3，···，N+m−1}，如 $ {\lambda }_{12} $=3表示客户2与客户1邻近度为3，即在所有客户和车场中，客户2距离客户1第3近，当i= j时，定义 $ {\lambda }_{ij} $= 0.

定义概率选择公式为

(10) $ P\left(i,j\right)=\frac{1/{\lambda }_{ij}}{\displaystyle\sum\limits_{j = 1}^{{C_k} + 1} {1/{\lambda _{ij}}}} ;\;i\in {C}_{k},j\in {C}_{k}\cup \left\{k\right\},i\ne j. $

式中： $ {C}_{k} $表示为车场k所要服务的客户，当车辆服务客户i之后，根据概率公式选择将要服务的下一个客户j，客户j距离客户i越近，则被选中的几率越高，为了保证种群的多样性，其他距离较远的客户依然有一定的概率被选中，然后依次选取客户j的下一个客户，直至该车辆饱和，无法承担更多客户为止.

构造初始解的步骤如下所示.

1 初始化，capacity(h)=0， $ {R}_{h}=\left\{{h}\right\} $， $ {h}_{c}=H $ //capacity(h)为车辆h的当前装载量， $ {R}_{h} $为车辆h的配送路线， $ {h}_{c} $为车场k的车辆集

2 for k=1∶m

3 while not is empty( $ {C}_{k} $)

4 随机在 $ {h}_{c} $中选择一辆车h

5 $ {h}_{c} $= $ {h}_{c}\backslash h $// 表示将车辆h从车辆h_c中移除

6 根据概率公式在 $ {C}_{k} $中选择一个距离车场k较近的客户i分给车h

7 $ {C}_{k}={C}_{k}\backslash i $，将i加到 $ {R}_{h} $尾部

8 While 车辆h没有满载时

9 根据概率公式选择该车的下一个服务客户j

10 if ${\rm{capacity}}\left( {{h}} \right) + {r_j} \leqslant {Q_h}$

11 $ {C}_{k}={C}_{k}\backslash j $，将j加到 $ {R}_{h} $尾部

12 else

13 break

14 end if

15 end while

16 end while

17 end for

如图3所示，将初始种群按照适应度大小进行升序排列，位列第1的青蛙被分配给第1个族群，位列第2的青蛙被分配给第2个族群，以次类推，直至将所有青蛙分别分配到N_f个族群中，每个族群中包含S_f个青蛙.

在传统的SFLA中，青蛙的学习策略为向族群内最优个体P_b跳跃或者向整个青蛙种群内最优个体P_g跳跃，然而这种进化策略并没有充分利用到其余青蛙个体的信息，无法进行全面有效的搜索，导致算法易提前陷入局部最优解. 为了充分利用族群内部的有效信息，增加算法的搜索宽度，受小生境的启发，在族群内部引入子群的概念，在大小为n的族群内部随机选择s只青蛙组成一个子群（s<n），子群内最优个体为P_b，最差个体为P_w，青蛙种群最优个体为P_g，青蛙在跳跃过程中，信息交流方式借鉴了遗传算法中的交叉模式，同时，为了提高信息交换的质量，保证信息的有效性，车辆只与同车场的车辆进行信息交流，其具体跳跃步骤如下.

1）在族群中随机选择s个青蛙（s<n）形成一个子群，子群内最优个体为P_b，最差个体为P_w，青蛙种群最优个体为P_g；

2）P_w向P_b跳跃，具体过程如图4所示.

a）在P_b中随机选择某个车场中某辆车h的基因链，记为b，在P_w中也选择相同车辆的基因链，记为w；

b）分别将b、w这2条基因链中共有的客户基因存放到基因库F_bw中，而将其独有的基因分别存放在基因库F_b和F_w中；

c）交换b、w这2条基因链的客户基因；

d）遍历P_w中基因链w以外的其他客户基因，删除所有与F_b共享的基因，以此类推，对P_b个体采取同样的操作；

e）随机选择F_w中的一个客户基因，并将其插入到与车辆h相同车场的车辆的基因链中，若超载，则随机插入到该车场另一辆车的基因链中，如果该客户基因无法被该车场的任意一辆车进行服务，则随机选择该车场的某个客户并将其加入到F_w中，然后继续对该客户进行分配，直至将该客户分配出去. 对F_w中每个客户均执行上述操作，直至F_w为空. 循环若干次，如果仍然有客户无法被分配出去，则跳出循环，保持个体P_w不变；以此类推，对P_b个体采取同样的操作；

f）对其余m−1个车场均执行以上操作；

3）P_b跳跃过后的青蛙命名为 $ {P}_{{\rm{b}}}^{'} $，计算其适应度f( $ {P}_{{\rm{b}}}^{'} $)，与f(P_b)进行比较，如果f( $ {P}_{{\rm{b}}}^{'} $)>f(P_b)，用 $ {P}_{{\rm{b}}}^{'} $替换P_b，否则保持P_b不变；

4）P_w跳跃过后的青蛙命名为 $ {P}_{{\rm{w}}}^{'} $，计算其适应度f( $ {P}_{{\rm{w}}}^{'} $)，与f(P_w)进行比较，如果f( $ {P}_{{\rm{w}}}^{'} $)>f(P_w)，用 $ {P}_{{\rm{w}}}^{'} $替换P_w，否则将步骤2）中的P_b替换成P_g重新进行跳跃操作，再次比较跳跃过后的f( $ {P}_{{\rm{w}}}^{'} $)与f(P_w)，如果跳跃后的适应度有所提高，即f( $ {P}_{{\rm{w}}}^{'} $)>f(P_w)，用 $ {P}_{{\rm{w}}}^{'} $替换P_w，否则，随机生成一只青蛙替换P_w.

聚类分析并不能完全准确地将客户分配给对应的车场，导致每个车场内部所形成的车辆服务路线从整体来看并不是最优的，因此，车场之间必须进行信息交流，对部分客户进行重新分配，从而提升找到最优解的概率.

同构造邻近矩阵一样，构造远离矩阵，远离矩阵用来描述客户与车场的远离程度，通过客户与车场之间的距离构造远离矩阵：

(11) $ S\left( {x \leftarrow y} \right) = {({S_{xy}})_{m \times {m_x}}};x = 1,\; \cdots ,\;m;\;y = 1, \;\cdots ,\;{m_x}. $

式中： $ {S}_{xy}$为客户y和车场x的远离程度， $ {S}_{xy}\in ${1，2，3，···， $ {m}_{x} $−1}，如 $ {S}_{12} $=3表示客户2与车场1的远离度为3，即在车场1所服务客户中，客户2距离车场1第3远； $ {m}_{x} $为车场x所要服务的客户个数.

定义概率选择公式为

(12) $ {{P}}\left( {{{x}},{{y}}} \right) = \frac{{1/s\left( {x,y} \right)}}{\displaystyle\sum\limits_{y = 1}^{{m_x}} {1/s\left( {x,y} \right)} };x = 1,\; \cdots ,\;m;\;y = 1, \;\cdots ,\;{m_x}. $

根据概率公式，选择距离车场x较远的客户y，进行车场之间的信息交流. 在车场在进行信息交流时，主要有以下2种方式：1）客户互换. 分属不同车场的2个客户互换车场，如图5所示，根据式（12）选择距离车场A较远的客户y和距离车场B较远的客户z，互换2个客户. 2）客户搬迁. 将原本属于某个车场的客户重新分配给另一个车场进行服务，如图6所示，根据式（12）选择距离车场A较远的客户y，然后根据式（10）将客户y分配给距离其较近的车场B进行服务.

邻域搜索的具体步骤如下.

1 输入：青蛙P

2 for c=1∶m

3 while 不满足搜索次数

4 if rand (0,1)<0.5

5 根据概率公式（12）选择较远客户y

6 根据概率公式（10）将客户y分配给较近车场，形成新解 $\small {P}^{'}$

7 if f（ $\scriptsize {P}^{'}$）>f（P）

8 $\scriptsize P={P}^{'}$

9 else if 满足搜索次数

10 以模拟退火的接受准则决定是否接受新解

11 if 接受新解

12 保持P不变，用 $\scriptsize {P}^{'}$替换族群中最差青蛙

13 end if

14 end if

15 else

16 根据概率公式（12）选择较远客户y

17 选择另一车场，根据概率公式（12）选择其较远客户z

18 将y和z互换车场，形成新解 $\scriptsize{P}^{'}$

19 if f（ $\scriptsize {P}^{'}$）>f（P）

20 $\scriptsize P={P}^{'}$

21 else if 满足搜索次数

22 以模拟退火的接受准则决定是否接受新解

23 if 接受新解

24 保持P不变，用 $\scriptsize {P}^{'}$替换族群中最差青蛙

25 end if

26 end if

27 end if

28 end while

29 end for

为了充分发挥较劣解中的优良信息，采取模拟退火中的接受准则来一定程度上地接受邻域搜索过程中所产生的较劣解，提高算法逃离局部最优解的能力. 最后算法接受的解的表达式如下：

(13) $ {P_{{\rm{new}}}} = \left\{ {\begin{array}{*{20}{l}} {P',\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;f\left( {P'} \right) < f\left( {{P}} \right)}\\ \;\;\;\;\;\;\;{{\text{或}}\;\min\; \left( {1,\exp\; \left( {\dfrac{{ - \left( {{{f}}\left( {P'} \right) - {{f}}\left( {{P}} \right)} \right)}}{{{T_g}}}} \right)} \right) \geqslant \gamma ;}\\ {P,\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\text{其他}}.} \end{array}} \right. $

式中：P为邻域搜索前的解， $ {P}^{'} $为邻域搜索后的解， $ f\left(P\right) $为个体P的适应度， $ {T}_{g} $为第g次迭代时的温度， $ \gamma $为[0，1.0]的随机数.

将所有族群内的青蛙进行混合，删除其中完全相同的青蛙，只保留一只在种群中，然后随机生成若干只青蛙，使蛙群的总数量保持不变.

在IHSFLA中，将每一代种群看作一种状态，种群的进化则可以看作状态间的转移，由蛙跳策略可知，当前种群的状态仅和上一代种群状态有关. 设总体状态空间为G，则在算法迭代过程中，每一代种群G（t）对应马尔科夫链模型中的一种状态，种群更新则对应马尔科夫链模型中不同状态的转移. 由混合蛙跳算法的收敛性^[26]可知，当族群内出现最佳个体时，子代个体将不断靠近最优个体：

(14) $\mathop {{\rm{lim}}}\limits_{{{t}} \to \infty } P\left\{ {{X_{\rm{b}}} \in G(t)} \right\} = 1.$

式中：X_b为全局最优解。 式（14）表示IHSFLA将最终收敛于全局最优.

为了充分验证本研究算法的有效性，对MDHVRP问题以及RVRP问题分别进行求解验算，利用MATLAB R2014b仿真实验平台进行2组实验测试，测试环境为Intel（R）Core（TM）i7-4712MQ CPU 2.30 GHz，8 GB RAM.

参照文献[27]的例子进行求解验算，具体描述如下：有一个包含了3个车场、3种类型共10辆车的配送中心要为其服务区内的30位客户进行相关配送服务，要求安排合理的车辆配送路线，使得所有车辆的总成本最小.

根据文献[28]设置改进混合蛙跳算法参数如下. 蛙群规模F=300，族群数N_f=20，族群中的青蛙数目S_f=15，每个族群的局部搜索次数N_s=10，邻域搜索次数N_n=5，子群规模S_z=12，种群最大迭代次数G=300，初始温度θ=1000 °C，降温速率q=0.9.

根据文献[27]设置对比算法参数如下. 多染色体遗传算法：种群规模Pop_size=400；个体间交叉概率P_c1=0.5；个体内交叉概率P_c2=0.1；变异概率P_m=0.03. 单染色体遗传算法：种群规模Pop_size=100；交叉概率P_c=0.6；变异概率P_m=0.2，精英选择比率P_s=0.1.

利用Matlab软件对该问题进行计算求解，车辆最优配送方案如表1所示，最优路径如图7所示，配送总费用为8 019元，车辆平均利用率为79%，算法的进化曲线如图8所示，在第142代开始收敛. 图中，I为迭代次数，C为总费用. 韩胜军^[27]利用多染色体遗传算法对该问题进行求解，最优路径如图9所示，配送总费用为8306元，算法在第591代收敛. 对比得出，本研究IHSFLA求解质量更高，同时在收敛速度上较多染色体遗传算法更快. 将算法运行100次，将其结果与文献[27]结果比较，如表2所示. 可以看出，与多染色体遗传算法相比，解的平均质量提高了11%，与单染色体遗传算法相比，解的平均质量提高了42%，与粒子群算法相比，解的平均质量提高了33%. IHSFLA 不仅有着更高的求解质量，而且算法的波动性更小，稳定性更高，充分说明了IHSFLA在求解此类问题上的优异性.

某石油企业有A、B、C、D共4个配送中心，产品成本互不相同，共有13辆大型、中型、小型3种不同类型的车辆，现有50位客户需要进行石油配送服务，其具体信息如表3、4所示，要求安排合理的车辆及其配送的行驶路线，使产品总成本及所有车辆配送成本之和最小，使企业利润最大化.

为了验证本研究初始解构造方法的有效性，先后用随机初始解、先聚类后随机方式初始解以及先聚类后邻近矩阵初始解3种不同的初始解构造方式对该问题分别进行求解，确定蛙群规模总数量F=400，族群数N_f=20，族群中的青蛙数目S_f=20，每个族群的局部搜索次数N_s=10，邻域搜索次数N_n=5，子群规模S_z=16，种群最大迭代次数G=400，初始温度θ=1000 °C，降温速率q=0.9.

3种方式产生初始解的平均分布情况如图10所示.可以看出，随机方式产生初始解的成本普遍约为22500，而先聚类后随机方式产生初始解的成本约为15000，较随机方式产生的初始解质量已经有了较大的提升，而先聚类后邻近矩阵产生初始解的方式产生的初始解更优，成本约为12500. 可以看出聚类对初始解质量的提升，同时，构造邻近矩阵的方式会进一步优化初始解，而好的初始解会有效加快算法的迭代速度，减少算法的运行时间.

运用3种方式分别对该问题进行多次求解，其平均进化曲线如图11所示.随机产生初始解的方式因为初始解质量较差，收敛速度较慢，在400次的迭代过程中，并没有进行收敛，而先聚类后产生初始解的方式明显更优，其中，先聚类后利用邻近矩阵产生初始解的方式所产生的初始解质量最高，收敛速度最快，在约第250代已经收敛，最优解为10 743. 最优配送方案如表5所示，车辆平均利用率为90%，最优配送路线如图12所示，各车辆配送路线基本没有交叉重叠的部分，也间接说明该配送方案的质量较高.

为了验证所研究问题的实际意义，在不考虑产品成本的情况下，对该问题进行求解，在重复求解100次后，平均总成本为11467元，而在考虑产品成本时，平均总成本为10860元，相对于不考虑产品成本的最优调度方案，考虑成本的最优调度方案总成本减少了6%，占产品成本的13%.

如图13所示为不考虑产品成本的最优调度方案. 当不考虑产品成本时，配送成本为5937元，总成本为10979元；当考虑产品成本时，最优调度方案如图12所示，配送成本为6002元，总成本为10743元. 对比可知，不考虑产品成本的最优调度方案虽然比考虑产品成本的最优调度方案配送成本略少，但是，总成本却多了236元，占产品成本的5%，也就是说，考虑产品成本的调度方案更符合公司的利益要求，可以取得更大的经济效益. 同时，随着客户订单的增多，调度次数的增加，考虑产品成本的调度方案的优越性会体现的更加明显，说明所研究问题具有较强的现实意义.

为了充分验证所提算法的有效性以及优越性，将其测试结果与传统蛙跳算法（SFLA）、单循环结构综合学习细菌觅食优化算法^[14]（SRCLBFO）、增强蚁群算法^[29]（EACO）进行比较分析，SFLA参数与本研究所提算法的一致，SRCLBFO、EACO参数如下. SRCLBFO^[14]：菌落数S=30，趋向操作次数N_c= 5，复制操作数N_re= 5，迁徙操作数N_s= 5，惯性因子ω= 0.9，学习概率P_c= 0.1，加速因子c₁= 1.2，c₂= 0.6. EACO^[29]：信息素挥发系数控制因子γ=1，信息素挥发系数初始值ρ₀=0.9，残留信息素重要程度α= 1.25，启发信息素重要程度β= 2.5，初始化信息素浓度参数P_m= 1.1，信息素增量常数W= 500.

将4种算法分别运行100次，求其最优解的平均值、标准差以及算法所用时间，其中算法时间以最优解连续100代没有变化或者迭代次数达到1000来计算. 输出其1000次的平均进化曲线，各算法比较结果如表6所示，平均进化曲线如图14所示.

由表6可以看出，相较于其他3种算法，本研究所提出的改进混合蛙跳算法均值和标准差都是最优的，表明所提算法在有着较好寻优能力的同时，又有着较好的稳定性. 同时，由图14可以看出，SFLA和SRCLBFO由于初始解较差，迭代速度较慢，耗时较长，在第600代后才慢慢进行收敛，且结果相对较差，而EACO虽然有着较强的进化能力，收敛速度较快，用时最短，但是通过2层聚类算法简化问题较易陷入局部最优解. 综合来看，本研究所提IHSFLA算法进化速度较快，用时较短，又有较强的全局搜索能力，可以有效求解此类问题.

（1）针对多中心分布式企业存在的不同配送中心产品成本不同的现实情况，以车辆固定成本、行驶成本、产品成本等综合成本为基础，建立以总成本最小化为目标的多约束车辆路径问题模型.

（2）为了求解该问题模型，设计改进的混合蛙跳算法. 该算法利用四标准聚类方法对客户进行聚类，缩小求解规模；利用邻近矩阵产生质量更高的初始解，加快进化速度；在青蛙族群内部设计子群，弥补了传统蛙跳算法只依靠局部最优解和全局最优解进行进化的缺点，充分利用了族群内部青蛙的有效信息；设计先内后外的交流模式，使车辆之间的信息交流更有方向性；利用远离矩阵进行邻域搜索，在增加算法搜索宽度的同时，使搜索更有针对性，更具有效性；引用模拟退火接受准则，充分利用较劣解中的有效信息，增加算法逃离局部最优解的能力. 算例对比分析表明，改进的混合蛙跳算法收敛速度快、搜索能力强，具有较好的精度与稳定性.

（3）通过对比分析考虑产品成本前、后的调度方案，发现考虑产品成本的调度方案较不考虑之前总成本平均减少6%，占产品总成本的13%，证明考虑产品成本在车辆调度中的必要性，为多中心分布式企业提供了更好的调度思路，更有利于其实现利润的最大化.

节能减排是最近社会一直倡导的，为了响应国家号召，发挥企业的社会价值，同时结合企业自身的效益，未来将会对考虑碳交易成本的车辆调度问题作进一步的研究.