随着烟草行业的持续发展，烟草制品的消费需求逐渐向更多元化和分散化的趋势发展. 在某些地广人稀、零售户分布不均衡的地区，由于地理位置和空间环境的限制，人工进行配送路径规划存在货物无法及时送达、配送成本极高的问题，造成零售户经营受影响以及物流中心成本过高的问题. 为了满足偏远地区分布不均衡的客户的需求，研究卷烟的配送路径. 提高配送效率、降低配送成本是卷烟配送过程中亟待解决的问题. 在地广人稀、烟草零售户分布不均的特定环境下，配送距离从1公里至数百公里不等，若按照现有的配送方案一一为客户进行配送，配送距离较长、配送次数较多，导致效率低下、配送成本增加. 现有配送主要通过货车公路运输，会导致汽车尾气排放增加. 因此，本研究引入快递行业的自提点模式，以期降低物流中心配送成本. 不仅给予零售户选择的权利、提高了配送的灵活性，也能降低配送成本、减少对环境的污染. 因此，本研究开展了针对地广人稀、零售户分布不均衡区域卷烟自提点选址及配送路径规划的研究. 自提点选址决定了客户从何处取货，而路径规划则决定了如何高效地将货物送到各个自提点. 两者结合研究有助于大幅度地降低运输成本.

近年来，车辆路径规划（vehicle routing problem，VRP）在卷烟配送领域的应用研究取得了显著进展. 在VRP问题的求解方面，与传统精确算法相比，启发式算法通过局部优化、邻域搜索和概率性跳出局部最优等机制，能够有效兼顾解的质量与计算效率. Archetti等^[1]针对分体配送车辆路径问题，提出变体的分体配送车辆路径，改进了启发式算法的精确性；范文兵等^[2]针对启发式算法收敛速度慢、配送成本高的问题，改进了传统的启发式规划算法，以跳出局部最优，在求解过程中取得了良好的效果；方文婷等^[3]针对蚁群算法初始阶段由于信息素不足导致收敛速度慢的问题，将A*算法与蚁群算法相结合，利用A*算法的全局收敛性和蚁群算法的正反馈性构造了一种混合蚁群算法. TAN等^[4]提出混合多目标进化算法，该算法在进化搜索中融入了多种启发式算法以进行局部开发，并引入帕累托最优性的概念来解决VRPTW中的多目标优化问题. 方昕^[5]提出的引入惯性权重的新型粒子群算法在都市圈多站点配送中表现优异，但其节点数太少. 上述研究多聚焦于零售户分布密集或半密集区域如城市商圈，对地广人稀、零售户分布高度不均衡区域的适配性不足. 快递行业通过引入自提点来降低物流中心的配送成本^[6]. 为了降低偏远区域配送成本，烟草公司也可借鉴该方式. 李谊等^[7]提出基于多视角聚类的配送中心选址模型，但其聚焦于宏观物流网络，未细化至零售户优先级；Zhu等^[8]的ZK-means算法与Lin等^[9]的密度峰值改进方法，虽然提升了聚类稳定性，但其距离度量仍以物理空间维度为主，未融合实际配送中的经济指标. 在卷烟实际配送过程中，还应考虑零售户的实际经济效益、历史经营情况. 因此，需要符合卷烟配送实际情况的自提点选址方法. 此外，尽管部分学者利用改进启发式算法求解卷烟配送问题，但其优化目标仍以单一成本或时间窗约束为主，缺乏对多目标协同优化的系统考量. 例如，Archetti等^[1]的分体配送模型虽能降低车辆空载率，却未考虑时间窗违约与零售户优先级惩罚成本的耦合效应；Li等^[10]与张平莉^[11]的研究在客户优先级建模上取得进展，但其优先级定义多依赖静态指标如订单量，未结合零售户的历史经营动态与区域经济特征，可能会导致稀疏分布场景下的配送灵活性不足. 上述局限性导致现有选址方案在稀疏区域易出现自提点覆盖不足或运营成本过高的问题. 多目标路径规划研究近年来逐步深入，但在地广人稀场景下的协同优化方面仍面临挑战.

现有研究方法在地广人稀、零售户分布不均衡的烟草物流场景中存在系统性局限. 选址问题的传统数学模型为P-中值数学模型^[12],依赖均匀分布假设，与本研究的地广人稀“孤岛式”零售网络会产生冲突. K-means类梯度算法是解决选址问题的经典算法，但其单纯依赖物理距离的聚类方式更易忽视低密度区域的经济效益与优先级需求，造成选址决策偏差. VRP的数学模型是混合整数规划模型，节点数增长会增加求解难度. 在偏远区域长距离的实际配送中，单次服务零售户数量增加会显著增加车辆运输成本，且在分布不均衡的场景下，高优先级零售户与低优先级零售户在时间窗要求上可能会存在冲突.

本研究以最小化车辆总成本、优先级成本以及自提点总成本为目标，构建带时间窗的多目标卷烟选址-配送路径规划模型. 主要创新点如下. 通过物流水平得分改进距离函数，提出经济K均值聚类（economic K-means，EKM）算法，引入自提点补贴机制，构建自提点选址模型，结合多阶段聚类与补贴机制实现经济效益与历史经营数据的综合优化，得到合理的卷烟自提点位置. 引入零售户优先级的概念，提出优先级模拟退火（priority simulated annealing，PSA）算法，将时间窗约束转化为优先级惩罚成本的软性目标，构建兼顾优先级成本最小化与车辆总成本最小化的多目标VRP模型. 最后，通过对Z地市区域进行卷烟配送路径规划，验证卷烟带时间窗的多目标选址-配送路径规划模型的有效性，并通过对比实验验证算法的优越性.

卷烟配送，即车辆从配送中心出发，按预设路线完成零售户送货后再返程的闭环作业. 卷烟公司基于需求量将零售户划分为不同的等级. 高等级客户因需求量大需优先配送以降低时效损失. 针对配送过程中配送效率低、成本高的问题，以及根据Z地市区域现有状况，构建自提点选址模型和带时间窗的多目标车辆配送路径规划模型. 问题描述如下.

（a）基本信息：零售户坐标、需求量、订单价格、车辆等信息均已知.

（b）距离函数：计算物流水平得分，得出新的距离函数.

（c）自提点位置：提出EKM聚类算法进行多次聚类，得到最优自提点位置，提出自提点补贴机制.

（d）配送方式：1-N类型，即一个配送中心为多个目标零售户点以及自提点提供配送服务.

（e）惩罚成本：在整体配送过程中，如果未按照优先级进行配送，就会产生惩罚成本，惩罚成本用时间窗来定义. 本研究的时间窗类型为软时间窗，若车辆未在规定时间内到达，则会产生相应的早到或晚到的优先级惩罚成本.

为了方便分析和研究，做出如下假设.

（a）配送车辆对每个顾客点只访问一次，且最终返回配送中心；

（b）配送车辆和自提车辆的型号均相同，并且车辆在运输过程中的速度始终保持不变，忽略运输过程中的堵车、天气恶劣的情况；

（c）车辆固定成本周期内不变，车辆每公里运输费用已知且保持不变；

（d）零售户群可选择送货上门或用户自提服务，同一零售户仅能选择一种服务方式，不存在自提的零售户临时改为要求送货上门的情况；

（e）零售户优先级的评定仅考虑需求量，不受其他因素影响；

（f）自提补贴仅针对距离自提点特定范围内的顾客生效，大于指定配送范围时无补贴.

为了得到合理的自提点位置，通过物流水平得分对每个零售户计算其经济效益. 其物流水平指标涵盖需求量、经理评分及信用等级等. 采用熵权法进行客观赋权：基于指标数据离散度，通过信息熵计算初始权重并修正，形成兼顾数据分布特性的综合权重体系. 具体表达式如下.

1）利用极差法，得到归一化判断矩阵. 令$ {X}_{ij} $、$ {Y}_{ij} $分别表示第$ i $个零售户第$ j $个正向指标和逆向指标，$ i=1{,}2,\cdots ,n;j=1{,}2,\cdots ,m $，n为零售户数量，m为指标个数.

正向指标、逆向指标表达式分别如下：

(1)$ {X}_{ij}=\frac{{X}_{ij}-\min\; {X}_{ij}}{\max \;{X}_{ij}-\min \;{X}_{ij}}, $

(2)$ {Y}_{ij}=\frac{\max\; {X}_{ij}-{X}_{ij}}{\max\; {X}_{ij}-\min\; {X}_{ij}}. $

2）计算第$ i $个零售户第$ j $个指标的比重$ {P}_{ij} $：

(3)$ {P}_{ij}={{Y}_{ij}}\Big/{\displaystyle\sum \limits_{i=1}^{n}{Y}_{ij}}. $

3）通过$ {P}_{ij} $计算第$ j $个指标的信息熵值$ {e}_{j} $：

(4)$ {e}_{j}=-\frac{1}{\ln \, n}\sum \limits_{i=1}^{n}{P}_{ij}\ln\; {P}_{ij}. $

4）得到第$ j $个指标的信息熵冗余度$ {d}_{j} $：

(5)$ {d}_{j}=1-{e}_{j}. $

5）计算第$ j $个指标权重$ {W}_{j} $：

(6)$ {W}_{j}={{d}_{j}}\Big/{\displaystyle\sum \limits_{j=1}^{m}{d}_{j}} .$

6）根据各项指标权重$ {W}_{j} $计算各个零售户$ i $的物流水平得分$ {Z}_{i} $：

(7)$ {Z}_{i}=\sum \limits_{j=1}^{m}{P}_{ij}{W}_{j}. $

由于零售户的级别存在差异，优先服务高需求客户可提升整体满意度. 若配送时未按照零售户优先级的顺序配送，会导致提前送达或超时配送，则会产生一定的优先级惩罚成本. 单位时间优先级惩罚成本是指车辆未在限定时间内送到零售户点，每个单位时间会产生额外的优先级惩罚成本.

假定单位时间优先级惩罚成本由分段函数表示，表达式如下：

(8)$ p_{i}'=\begin{cases} {c}_{0},\;{f}_{0}\leqslant {h}_{i}\leqslant {f}_{1};\\ {c}_{1},\;{f}_{1}\leqslant {h}_{i}\leqslant {f}_{2};\\ {c}_{2},\;{f}_{2}\leqslant {h}_{i}\leqslant {f}_{3};\\ {c}_{3},\;{f}_{3}\leqslant {h}_{i}\leqslant {f}_{4};\\ {c}_{4},\;{f}_{4}\leqslant {h}_{i}\leqslant {f}_{5};\\ {c}_{5},\;{f}_{5}\leqslant {h}_{i}\leqslant {f}_{6}.\\ \end{cases} q_{i}'=\begin{cases} c_{0}',\;f_{0}'\leqslant {h}_{i}\leqslant f_{1}';\\ c_{1}',\;f_{1}'\leqslant {h}_{i}\leqslant f_{2}';\\ c_{2}',\;f_{2}'\leqslant {h}_{i}\leqslant f_{3}';\\ c_{3}',\;f_{3}'\leqslant {h}_{i}\leqslant f_{4}';\\ c_{4}',\;f_{4}'\leqslant {h}_{i}\leqslant f_{5}';\\ c_{5}',\;f_{5}'\leqslant {h}_{i}\leqslant f_{6}'.\\ \end{cases} $

$ \mathrm{式中}：{c}_{y} $和$ {{{c}_{y}}}' $分别表示历史和现有评定等级为$ y\;(y=0{,}1,2{,}3,4{,}5) $的零售户的单位时间优先级惩罚成本，$ {h}_{i} $表示零售户$ i $的需求量，$ {f}_{0} $~$ {f}_{6} $、$ f_{0}' $~$ f_{6}' $为常数.

零售户优先级级别由零售户过去的评级和现有的评级共同构成，具体表达式如下：

(9)$ {G}_{i}=\alpha p_{i}'+\beta q_{i}', $

(10)$ \alpha +\beta =1. $

式中：$ p_{i}' $和$ q_{i}' $表示零售户历史和现有的单位时间优先级惩罚成本；$ {G}_{i} $代表第$ i $个零售户的单位时间优先级惩罚成本；$ \alpha $和$ \beta $分别代表零售户历史等级所占的权重系数和零售户现有等级所占的权重系数，总和为1.

为了进一步激励商家选择自提^[6]，提出补贴机制. 补贴机制是指在零售户选择自提的情况下，通过提供一定金额的补贴，鼓励商家积极选择自提，从而减少物流成本并提升整体效率. 本研究补贴机制根据卷烟企业历史订单和配送情况，与企业讨论后开展试点工作，补贴金额以客户订购金额的一定比例为基础，且规定补贴金额不超过订购金额的0.6%. 由于卷烟零售户的订购卷烟量都存在一定限额，故在自提点范围内用户可以选择步行或者骑电瓶车前往自提. 根据实际情况，零售户一般一周订购一次，一次平均50条烟，订购金额约0.5万元至1.0万元，补贴金额约为零售户购买金额的0.3%~0.6%. 针对距离自提点$ {D}_{\max} $范围内的零售户提供一定的自提补贴. 本研究$ {D}_{\max} $设定为3 km. 补贴成本表达式如下：

(11)$ {S}_{5}=\left\{\begin{array}{ll}{w}_{2}{D}_{kj}-{C}_{{\mathrm{s}}}{D}_{kj},&{D}_{kj}\leqslant {D}_{\max};\\ 0 ,& {D}_{kj}> {D}_{\max}.\end{array}\right.$

式中：$ k $表示配送车辆编号，$ j $表示目标零售户的编号， $ {w}_{2} $表示车辆$ k $的单位距离行驶成本，$ {D}_{kj} $表示自提点到零售户的距离，$ {C}_{{\mathrm{s}}} $表示单位公里自提点提供的补贴金额.

建立自提点总成本最低的卷烟选址模型，考虑自提点建设及运营成本、自提点补贴成本，提出EKM聚类算法求解，得到最优自提点位置. 数学模型表达式如下：

(12)$ {\mathrm{Min}}\;{S}_{0}=\sum \limits_{e=1}^{b}\sum \limits_{j=1}^{N_{\text{{c}}}}\sum \limits_{l=1}^{a}{\zeta }_{ejl}d({S _{jl}},{H}_{el}). $

(13)$\begin{split} &{\mathrm{s.t.}} \\ &\displaystyle\sum \limits_{l=1}^{a}{\zeta }_{ejl}=1,\;\forall e\in \left\{1{,}2,\cdots ,b\right\}, \; \forall j\in \left\{1{,}2,\cdots ,{N}_{{\mathrm{c}}}\right\};\end{split}$

(14)$ \begin{split}& {H}_{el}={\displaystyle\sum \limits_{j=1}^{{N}_{{\mathrm{c}}}}{\zeta }_{ejl}{S }_{jl}}\Big/{\displaystyle\sum \limits_{j=1}^{{N}_{{\mathrm{c}}}}{\zeta }_{ejl}},\\&\;\forall e\in \left\{1{,}2,\cdots ,b\right\},\;\forall l\in \left\{1{,}2,\cdots ,a\right\}; \end{split}$

(15)$ {M}_{jp}=\left\{\begin{array}{ll}1,& \displaystyle\sum \limits_{l=1}^{a}{\zeta }_{ejl}d({S}_{il},{H}_{el})\geqslant g,\;\forall j\in \left\{1{,}2,\cdots ,g\right\};\\ 0,& 其他.\end{array}\right. $

式(12)为目标函数$ {S}_{0} $，表示通过调整聚类中心来最小化数据点与其所属聚类中心之间的距离的平方和最小；式(13)为数据点分配约束，确保每个数据点都只能被分配到一个聚类中心；式(14)为中心更新约束，确保每个阶段的聚类中心是通过在该阶段分配给聚类的数据点的加权平均计算得到的; 式(15)表示只有当某个数据点距离自提点$ g $公里以外时，相应的二进制变量$ {M}_{jp} $才等于1. 其中，$ {S }_{jl} $表述数据点$ ；{H}_{el} $表示所属聚类中心$ ；d({S }_{jl},{H}_{el}) $表示两者之间的欧氏距离的平方和$ ；{N}_{{\mathrm{c}}} $为表示包含配送中心与配送节点的集合；$ a $为聚类中心数量；$ b $为阶段数量；$ {\zeta }_{ejl} $为0-1变量，$ {\zeta }_{ejl}=1 $表示第$ e $次迭代时，零售户$ j $归属于第$ l $个聚类中心；$ i $和$ j $表示目标零售户的编号.

K-means聚类算法是一种无监督学习算法，有收敛速度较快、聚类效果较好的优点. 但传统K-means聚类算法按距离的相似性进行聚类，例如欧氏距离、曼哈顿距离，其仅仅是物理空间上的距离，未考虑零售户实际经营情况带来的影响. 因此，提出经济效益距离函数，表达式如下：

(16)$ D\left({x}_{i},{x}_{j}\right)={{{{D}_{ij}}}^{2}}\Big/{({Z}_{i}{Z}_{j})}. $

式中：$ {{{D}_{ij}}}^{2} $为空间指标；$ {Z}_{i} $和$ {Z}_{j} $分别代表任意2个零售户的得分情况，零售户间相互独立、互不影响.

如图1所示为融合经济效益指标的EKM算法流程图. 针对分散的零售户点，通过EKM聚类算法进行区域划分，通过多阶段聚类得到每个区域最合理的自提点位置，作为后续配送过程中待配送点之一. 零售户可就近选择自提，减少无法及时送达的情况.

轮廓系数是常用的聚类评估指标，用于衡量聚类的紧密度和分离度. EKM聚类算法应当能够生成更加紧密且更好分离的聚类结果. 因此，通过引入轮廓系数进行评估. 轮廓系数的取值范围为[−1.0, 1.0]，其计算方法如下.

1）对于每一个样本$ i $，计算同一个簇中其他样本与$ i $之间的平均距离，记作$ a(i) $.

2）对于每一个样本$ a(i) $，计算其最近的其他簇的所有样本与$ i $之间的平均距离，记作$ b(i) $.

3）样本$ i $的轮廓系数$ s(i) $定义为

(17)$ s(i)=\frac{b(i)-a(i)}{\max\; \left\{a(i),b(i)\right\}} .$

整个数据集最终的轮廓系数是通过计算所有样本轮廓系数的均值得到的. 轮廓系数越大越好，接近1.0表示较好的聚类效果，而接近−1.0表示较差的聚类效果.

Iris（鸢尾花）数据集是用于测试和评价算法性能的经典数据集，分别应用于K-means算法与EKM算法，得到轮廓系数对比结果，如表1所示. 由表1可知，随着聚类簇数的增加，EKM聚类算法的轮廓系数下降较缓，且在终点处系数较高，验证了算法改进的有效性. 当算法生成更接近1.0的轮廓系数时，表明簇之间更紧密且分离更好，即算法的改进是有效的. 算法对比见表1.

建立考虑零售户优先级成本、车辆总成本以及自提点总成本最低的带时间窗的多目标卷烟配送路径规划模型.

1）目标函数.

本研究根据不同零售户和自提点的需求量设定优先等级，并将其转化为单位时间优先级惩罚成本，同时基于惩罚时间窗，得到最终的优先级惩罚成本. 总成本由车辆的固定成本$ {S}_{1} $、车辆运输成本$ {S}_{2} $、优先级惩罚成本$ {S}_{3} $、自提点建设成本及运营成本$ {S}_{4} $、自提点补贴成本$ {S}_{5} $构成.

车辆固定成本表达式如下：

(18)$ {S}_{1}=\sum \limits_{i=1}^{{N}_{{\mathrm{c}}}}\sum \limits_{k=1}^{K}{w}_{1}{x}_{oik}. $

式中：$ {w}_{1} $表示每辆车的发车成本；K表示车辆集合；$ {x}_{oik} $为0-1变量，当车辆k从配送中心到需要配送的节点$ i $时为1.

车辆运输成本表达式如下：

(19)$ {S}_{2}=\sum \limits_{i=0}^{{N}_{{\mathrm{c}}}}\sum \limits_{j=0}^{{N}_{{\mathrm{c}}}}\sum \limits_{k=1}^{K}{d}_{ij}{w}_{2}{x}_{ijk}. $

式中：$ {x}_{ijk} $为0-1变量，当车辆k从零售户$ i $配送到$ j $时为1；$ {d}_{ij} $为零售户$ i $与$ j $之间的距离.

优先级惩罚成本表达式如下：

(20)$ {S}_{3}=\sum \limits_{i=1}^{{N}_{{\mathrm{c}}}}\sum \limits_{k=1}^{K}{x}_{ijk}\left\{\max \left[0,\left({d}_{ij}/v-{\sigma }_{i}\right)\right]\times{G}_{i}\right\} .$

式中：$ v $为车速，$ {\sigma }_{i} $为零售户$ i $可接受的车辆最晚到达时间.

自提点建设成本及运营成本表达式如下：

(21)$ {S}_{4}=X{C}_{{\mathrm{c}}}+X{C}_{0}{Z}_{ip}. $

式中：$ X $为自提点数量；$ {C}_{{\mathrm{c}}} $为自提点建设成本；$ {C}_{0} $为自提点运营成本；$ {Z}_{ip} $为0-1变量，当零售户$ i $被自提点$ p $服务时取1.

自提点补贴成本表达式如下：

(22)$ {S}_{5}=\sum \limits_{i=1}^{{N}_{{\mathrm{c}}}}\sum \limits_{p=1}^{{N}_{{\mathrm{p}}}}{Z}_{ip}\left\{\max \left[0,{w}_{2}{D}_{kj}-{C}_{\text{s}}{D}_{kj}\right]\right\}. $

式中：$ {N}_{{\mathrm{p}}} $为自提点集合.

2）数学模型.

(23)$ {F}_{2}=\min \;({S}_{1}+{S}_{2}+{S}_{3}+{S}_{4}+{S}_{5}). $

(24)$ {\mathrm{s.t.}} \;\; \;\;{\sum }_{p\in {{N}_{{\mathrm{p}}}}}{Z}_{ip}\leqslant 1,\;\forall i\in \left\{1{,}2,\cdots ,{N}_{{\mathrm{c}}}\right\}; $

(25)$ {\sum }_{p\in {{N}_{{\mathrm{p}}}}}{W}_{p0}\leqslant {\sum }_{p\in {{N}_{{\mathrm{p}}}}}{V}_{p}; $

(26)$ {C}_{{\mathrm{s}}}\geqslant 0; $

(27)$ {C}_{{\mathrm{s}}}\leqslant {w}_{2}{d}_{ij},\;\forall j\in \left\{1{,}2,\cdots ,{N}_{{\mathrm{c}}}\right\}; $

(28)$ {w}_{2}{D}_{kj}\leqslant {w}_{2}{d}_{ij},\;\forall i\in \left\{1{,}2,\cdots ,{N}_{{\mathrm{c}}}\right\},\;j\in \left\{1{,}2,\cdots ,{N}_{{\mathrm{c}}}\right\}; $

(29)$ \sum \limits_{k=1}^{K}{y}_{ik}=1,\;\forall i\in \left\{1{,}2,\cdots ,{N}_{{\mathrm{c}}}\right\};$

(30)$ \sum \limits_{k=1}^{K}\sum \limits_{i=0}^{{N}_{{\mathrm{c}}}}{x}_{ijk}=1,\;\forall j\in \left\{1{,}2,\cdots ,{N}_{{\mathrm{c}}}\right\};$

(31)$ \sum \limits_{k=1}^{K}\sum \limits_{i=1}^{{N}_{{\mathrm{c}}}}{y}_{ik}=n;$

(32)$ \sum \limits_{i=1}^{{N}_{{\mathrm{c}}}}{y}_{ik}{h}_{i}\leqslant Q,\;k\in \left\{1{,}2,\cdots ,K\right\} ;$

(33)$ \sum \limits_{j=1}^{{N}_{{\mathrm{c}}}}{x}_{0jk}=1,\;k\in \left\{1{,}2,\cdots ,K\right\} ;$

(34)$ \sum \limits_{i=1}^{{N}_{{\mathrm{c}}}}{x}_{0ik}=1,\;k\in \left\{1{,}2,\cdots ,K\right\}; $

(35)$ {\sum }_{i\in {{N}_{{\mathrm{c}}}}}{\sum }_{j\in {{N}_{{\mathrm{c}}}}}{x}_{ijk}\leqslant |N|-1,\;k\in \left\{1{,}2,\cdots ,K\right\} ;$

(36)$ {\varepsilon }_{i}\leqslant {T}_{i}\leqslant {\sigma }_{i}-{\theta }_{i},\;\forall i\in \left\{1{,}2,\cdots ,{N}_{{\mathrm{c}}}\right\} ;$

(37)$ {T}_{j}\geqslant {T}_{i}+{\theta }_{i}+{d}_{ij}/v,\;\forall i\in \left\{1{,}2,\cdots ,{N}_{{\mathrm{c}}}\right\},\;j\in \left\{1{,}2,\cdots ,{N}_{{\mathrm{c}}}\right\}; $

(38)$ \begin{split} &{x}_{ijk},{y}_{ik}\in \{0{,}1\},\;\forall i\in \left\{1{,}2,\cdots ,{N}_{{\mathrm{c}}}\right\},\\&\;j\in \left\{1{,}2,\cdots ,{N}_{{\mathrm{c}}}\right\},\;k\in \left\{1{,}2,\cdots ,K\right\}\end{split} .$

式(23)为目标函数，$ {F}_{2} $表示总成本最小；式(24)表示每个零售户最多只被一个自提点服务；式(25)表示每个零售户最多由配送中心访问一次，其中$ {W}_{p0}$、${V}_{p}$均为0-1变量，当自提点被配送中心服务时，$ {W}_{p0}$为0，当自提点p开放运营时，${V}_{p} $为0；式(26)为非负性约束，表示补贴金额始终是非负的；式(27)表示补贴金额不能超过物流中心配送到户的成本；式(28)表示自提成本与补贴金额之和不能超过物流中心总的配送到户成本；式(29)、 (30)表示每个零售户只能用一辆车服务；式(31)表示车辆遍历所有节点；式(32)表示每条线路上的零售户需求量不能超过车载能力；式(33)、(34)表示所有任务的起点是配送中心，终点也是配送中心；式(35)表示子回路消除约束；式(36) 表示车辆到达开始服务零售户$ i $的时间晚于最早到达时间$ {\varepsilon }_{i} $，早于最晚到达时间$ {\sigma }_{i} $与服务零售户$ i $的时间差，其中$ {T}_{i} $为到达开始服务零售户$ i $的时间，$ {\theta }_{i} $为服务零售户$ i $需要的时间；式(37)表示保证了车辆在时间窗内服务每个零售户；式(38)表示车辆$ k $从零售户$ i $行驶到$ j $，则$ {x}_{ijk} $=1，否则为0. 若零售户$ i $的需求由车辆$ k $来完成，则$ {y}_{ik} $=1，否则为0.

模拟退火算法（simulated annealing，SA）是基于Monte-Carlo迭代实现全局最优的一种算法^[13]. 本研究对SA进行如下改进：

1）路径构造. 对于n个零售户点，直接产生n个1~n间互不重复的自然数排列，表示优先级配送顺序，并依约束划分配送路径.

2）增加记忆函数. 记忆函数可帮助算法更好地记忆过去的状态，并根据这些信息来指导搜索过程，从而提高收敛速度和搜索效率，减少算法的运行时间. 增加记忆函数可以防止概率性的最优解错失，因此在算法运行中，不断地会产生并接受新解. 记忆函数可以基于过去的状态和性能来更新当前状态的选择. 其更新规则可以根据过去状态的性能和记忆程度来进行. 假设每次迭代都会得到一个新的候选状态，就可根据以下规则更新记忆函数.

(a)计算候选状态$ {s}' $的性能.

(b)如果$ {s}' $的性能优于当前状态$ M(s) $，则将$ M ({s}') $更新为一个较大的值，以表明该状态的重要性增加了.

(c)将当前状态$ s $的记忆程度也进行更新，以保留过去状态的信息. 这可以通过将$ M(s) $更新为一个介于当前记忆程度和新的性能值之间的值来实现，以平衡新信息的获取和过去信息的保留. 记忆函数中更新记忆函数$ M({s}') $的规则、更新当前状态$ M(s) $记忆程度的规则分别如下：

(39)$ M({s}')=\max \;(M({s}'),{\mathrm{per}}({s}')), $

(40)$ M(s)=\eta M(s)+(1-\eta ){\mathrm{per}}(s). $

式中： $ 1-\eta $为介于0和1.0之间的衰减系数，用于控制过去信息的保留程度.$ M(s) $和$ M({s}') $分别表示状态$ s $和状态$ {s}' $的重要性. 这个记忆程度可以被视为状态$ s $和状态$ {s}' $在搜索中的优先级. 高的$ M(s) $和$ M({s}') $意味着状态$ s $和$ {s}' $在搜索中更加重要，有可能更频繁地被选取或者更加影响搜索方向. 反之，则在搜索中的重要性较低，可能被更少地考虑或者被更快地淘汰.$ {\mathrm{per}}(s) $和$ {\mathrm{per}}({s}') $分别表示状态$ s $和$ {s}' $的性能或评价函数值. 这个评价函数通常用来衡量状态$ s $和$ {s}' $在解空间中的优劣程度. 在不同的问题领域中，评价函数$ {\mathrm{per}}(s) $可以有不同的定义，通常取决于具体的问题和优化目标，在本问题中代表最小化成本. $ {s}' $与之前的状态$ s $不同，$ {s}' $是由某种移动或转换操作产生的候选解. $ {\mathrm{per}}({s}') $表示问题的性质和优化目标，在本问题中代表优先级成本指标.

在SA中，通过评价函数的值来判断新状态$ {s}' $是否比当前状态$ s $更优秀. 如果$ {\mathrm{per}}({s}') $较小，则意味着新状态更好，算法有可能接受它作为当前状态，从而继续探索. 否则，算法可能以一定的概率接受差于当前状态的新状态，以避免陷入局部最优解. 通过引入记忆函数，SA可以利用过去的信息，指导正确的搜索过程，加快收敛速度，最终达到最优解.

3）增加判断函数. 引入判断函数可以增强新状态的随机性，同时也有利于提高算法的效率. 在初始状态的处理过程中，引入判断函数到不同种群. 针对小规模，采取随机全排列的方式；针对大规模，可以协作不同的变化算法. 判断函数通过优化Metropolis接受准则来实现：

（a）计算新状态$ {g}' $的性能表现.

（b）计算当前状态$ g $与新状态$ {g}' $之间的差异，通常表示为能量差或成本差.

(41)$ {\Delta }E={\mathrm{mani}}\;({g}')-{\mathrm{mani}}\;(g) .$

（c）根据$ {\Delta }E $的值来判断是否接受新状态$ {g}' $.

(42)$ \left.\begin{array}{ll}\text{如果Δ}E \lt 0,& \text{接受}{{{g}}'};\\ \text{如果}\Delta {E}\geqslant 0,& \text{以概率}{\mathrm{e}}^{-\Delta {E}/{T}}\text{接受}{{{g}}'}.\end{array} \right\}$

通过增加判断函数，可以更灵活地在搜索过程中进行状态接受与拒绝的决策，从而更有效地探索搜索空间并逐步收敛到全局最优解.

4）确定迭代步长. 本研究采用固定的迭代步长，即针对每个温度均迭代相同的步数.

采用整数编码的形式对路线进行编码. 算法采用该编码方式能直观地体现零售户的优先级别，如示例：0表示配送中心，如有2个零售户点，1个自提点，1辆车，零售户1的优先级惩罚成本为50，自提点z2的优先级惩罚成本为80，零售户3的优先级惩罚成本为35，则解可表示为0-z2-1-3-0.

初始解$ {p}_{0} $的生成步骤如下.

1）以配送中心为起点和终点，随即生成初始路径，分别计算所有零售户和自提点的优先级惩罚成本以规定其配送顺序.

2）计算每条路径的总距离，保证总距离不超过车辆的最远行驶距离.

3）判断每条路径所需配送的所有零售户以及自提点的总需求重量与车辆剩余装载量的比值$ \beta $. 若$ \beta $小于等于1，则该条路径零售户的总需求可以被满足；若$ \beta $大于1，则须拆分某些零售户的需求.

4）基于步骤1）的零售户配送序列，在考虑车辆行驶距离和载重量的基础上，依次计算每个零售户和自提点插入该路径前后所产生的成本，同时，分别计算所有零售户和自提点的优先级惩罚成本，最终选取成本最小的位置. 若零售户没有合适的插入位置，就建立一条新的路径.

5）重复以上步骤直至全部零售户的需求都被满足，增加记忆函数，基于过去状态的性能和记忆程度来进行，在搜索过程中引导状态的选择，最终输出最优解P，得到最优函数值f（P₀）.

选择3种邻域构造方法，分别对路径内邻域和路径间邻域进行构造. 路径内邻域构造是指仅对单条路径的零售户需求节点进行变换，路径内邻域构造分为路径内2-opt法、路径内点插入法和路径内点交换法.

1）路径内2-opt法：随机选择同一路径上任意2个零售户点，将2个配送点之间的所有配送点进行逆序，如图2(a)所示.

2）路径内点插入法：随机删除某条路径上的某个配送点，将该配送点插入该路径的其他任何一个位置，如图2(b)所示.

3）路径内点交换法：随机选取同一路径上的任意2个配送点，将这2个配送点互换位置，如图2(c)所示.

路径间的邻域构造是指对2条路径的节点进行变换. 通过扩大邻域搜索空间来提高解的质量，对可行解进行局部优化. 路径间邻域构造方法同路径内邻域构造相同.

1）路径间2-opt法：随机选取任意2条路径中任意2个零售户点，将2个零售户点之间的所有点进行逆序，如图3(a)所示.

2）路径间点插入法：随机删除某条路径上的某个零售户点，将该零售户点插入另一条路径的任何一个位置，如图3(b)所示.

3）路径间点交换法：随机选取任意2条路径中任意2个零售户点，将2个零售户点互换位置，如图3(c)所示.

通过编码、初始解生成和邻域构造，车辆从配送中心出发，依据各零售户的优先级惩罚成本，可以得到多条配送路线，具体如下：配送中心-自提点/零售户-配送中心.

PSA算法流程如图4所示.

本研究提出的PSA算法在传统SA框架基础上进行了如下改进. 1）记忆函数增强全局搜索：引入记忆函数式(39)、(40)，通过衰减系数平衡历史状态与当前解的性能差异. 该机制在迭代中记录优质解的特征，避免因随机扰动丢失潜在最优解. 2）多维度邻域构造策略：设计了3种邻域生成规则，路径内和路径间的2-opt优化与点插入和点交换. PSA兼顾收敛速度与解质量，解的表示直观且计算效率较高，在迭代过程中产生的解大多为可行解，算法收敛速度加快.

在VRPTW问题中，通常选用Solomon算例进行数值实验. Solomon算例中的R类数据集的点的坐标位置是随机生成的；C类数据集的点的坐标位置有明显的聚簇；而RC类数据集中点则同时具有随机和聚簇的特点，数据点分布较不均衡，符合文章研究问题的特点. 初始温度和温度下降系数是SA算法中须测试的关键参数. 故设置初始温度和下降速率的正交实验，每个参数设计3个因子级别，其值如表2所示. 其中，$ {T}_{0} $与$ {V}_{0} $分别为初始温度和下降速率. 对于每个参数组合，重复运行50次，将总成本作为评估参数组合性能的指标. 可以看出，初始温度为500，温度下降速率为0.995的参数组合效果最好.

以最小化总成本为目标，选取50个顾客规模，利用 PSA、SA、TS、GA这4种算法针对6个算例进行求解. 4种算法的求解性能结果如表3所示. 其中，RC101、RC102、RC103、RC201、RC202、RC203均为RC类算例，C表示最优成本，R表示最短路径. PAS的成本最低且运行时间最少，故PSA改进是有效的.

对Z地市区域内所有配送点的经纬度坐标、卷烟配送时间、订单需求量等数据进行规范化处理. 文中所有的经维度均指东经、北纬. 根据经纬度-公里转换规则，将零售户坐标经度和纬度均转化为公里数，纬度差1度或经度差1度对应的实际距离均约为111 km，例如经度39.84°和纬度98.57°可以表示为经度4422.12 km、纬度10941.7 km. 如表4所示为一个示例. 其中，X与Y分别表示经度(×111 km)、纬度(×111 km).

在案例参数中，自提点建设成本为5000元/个，自提点运营成本为2000元/月，自提点补贴成本为10元/km. 每辆车发车成本为200元，限载量为2000 kg，车辆单位配送成本为1.12元/km. 相关参数设置如下：初始温度T₀=500，终止温度T_S=1，温度下降速率r=0.997，最大迭代次数I=2069. 算法采用MATLAB2022b软件编程实现，计算机运行环境为WIN10系统，CPU为Intel(R) Core(TM)i7-7700k，内存为16 GB.

文章以Z地市区域所负责的200个点为案例进行研究，验证本算法模型的实际应用情况. 零售户现实的坐标分布见图5.

1）自提点选址：基于分布不均衡的卷烟零售户，依据EKM聚类算法进行多阶段聚类，分为区域1~5，得到13个卷烟自提点. 其中，区域1有2个自提点、区域2 有2个自提点，区域3有3个自提点，区域4有3个自提点，区域5有3个自提点. 卷烟自提点选址结果如图6所示，坐标及物流水平得分如表5所示.

2）筛选待配送零售户点：设置距离阈值为3 km，筛选出在阈值范围内选择配送的零售户、阈值范围外的零售户，得到待配送到户的零售户点共29个. 因此，待配送的零售户和自提点共42个，如表6所示. Z地市区域的历史优先级成本和现有优先级成本如图7、8所示. 其中，HQ、NQ、PC分别表示历史需求量、现有需求量与优先级成本.

3）最优配送路线图：待配送点编号表示为1、2、3、$\cdots $、41、42，自提点编号表示为z1、z2、$\cdots $、z13. 最优配送路线图如图9所示，共由4辆车配送.

通过50次随机重复实验^[14]，得到实验结果如下. 最优配送方案如表7所示，其中配送路线1需要配送4个自提点和10个零售户点，最优公里数为141.95 km，路径最优成本为158元，如图10(a)所示. 配送路线2需要配送2个自提点和16个零售户点，最优公里数为147.32 km，路径最优成本为165元，如图10(b)所示. 配送路线3需要配送3个自提点，最优公里数为202.68 km，路径最优成本为227元，如图10(c)所示. 配送路线4需要配送4个自提点和3个零售户点，最优公里数为725.00 km，路径最优成本为812元，如图10(d)所示. 总公里数为1216.95 km，总成本共计为1362元. Z地市区域原有的配送模式主要依赖人工分配车辆，按固定配送路线配送，短板明显，难以根据订单量、零售户数、配送距离等实时调整路线方案. 以案例中的200个零售户点为例，原有配送方案的总公里数为1502.65 km，总成本为1698元，所需车辆数为8. 通过实验表明，与传统人工配送方案相比，本研究提出的PSA算法通过优先级动态调整机制与多目标路径优化，显著降低了运输距离与车辆使用数量，从而实现了总成本节约近20%. 一是缩短了运输距离，PSA算法通过零售户优先级排序与路径内/间邻域搜索，有效避免了车辆绕行，运输距离较原方案减少了19.0%. 二是优化了车辆数量，原有方案依赖固定路线分配，而PSA算法通过整合高优先级零售户的配送任务，在相同载重约束下进一步减小了车辆空驶率，最终仅需4辆车即可完成全部任务，减少了车辆固定成本.

禁忌搜索算法（TS）的局部搜索能力较强，且不易陷入循环，广泛应用于路径规划问题^[15]. 因此，针对带时间窗的多目标选址-配送路径规划模型，对SA、TS与PSA进行对比. PSA算法得出总公里数为1216.95 km，总成本为1362元. SA算法求解结果如表8所示，总公里数为2109.27 km，总成本为2452元. TS算法最优路径结果如表9所示，总公里数为1649.30 km，总成本为1847 元. 3种算法实验所得的最优路径R、最优成本均值C、平均收敛迭代次数G及算法平均运行时间T的结果如表10所示，并根据算法迭代过程绘制了3种算法的平均收敛过程图，如图11所示. 其中，G表示迭代次数，C表示成本. 由图11可知，TS算法收敛较快，平均迭代次数为1450次，但对初始解依赖较强，难以得到最优解，算法最优路径为1694.30 km，最优均值为1847元，算法平均运行时间为388.0 s. SA算法收敛较快，平均迭代次数为1150次，但算法运行成本较高且时间较长, 算法最优路径为2109.27 km，最优均值为2362元，算法平均运行时间为410.4 s. 而本研究提出的PSA算法，虽然平均迭代次数为1500次，稍微长于前两者，但增强了局部搜索能力，其求解精度更高，稳定性更好，求得的成本更低，算法最优路径为1216.95 km，最优均值为1362元，算法平均运行时间为341.0 s，从而验证了PSA算法的有效性.

在地广人稀、分布不均衡区域的配送过程中，存在配送成本高、效率低、配送路径规划不合理的问题，因此，本研究提出自提点补贴机制鼓励零售户自提，提出零售户优先级定义，构造优先级惩罚成本，建立考虑零售户优先级成本、车辆总成本及自提点总成本最低的带时间窗的多目标卷烟选址-配送路径规划模型，并提出EKM-PSA算法进行求解，得到最优自提点位置和最优路径规划. 对实验案例进行50次随机重复实验验证了所建立的带时间窗的多目标选址-配送路径规划模型是可行有效的. 但是本研究缺乏对临时插单、路网实时变化的考虑，未来将针对动态路网下的配送进一步研究.