基于改进布谷鸟算法的物流配送中心选址研究

近年来，许多学者针对当下配送中心选址问题进行研究，主要包含单目标选址及多目标选址问题等[3]。陆丽丹，曹陆铖通过引入蚁群算法来对物流中心选址优化进行研究[4];刘善球，樊兵鹏将最小成本作为目标，运用遗传算法来解决快递物流配送中心选址问题[5];张聪等运用改进狮群算法以物流服务响应概率为约束条件研究了冷链配送中心选址问题[6]。但以上方法存在收敛速度、求解速度较慢，全局寻优能力较差等问题。

针对以上问题，本文提出了一种改进布谷鸟算法用于解决物流配送中心的选址问题，即将布谷鸟算法和粒子群算法进行有效组合，提高搜索效率，并引入自适应参数来避免算法过早收敛，增强全局寻优能力。实例仿真表明，本文的改进布谷鸟算法具有可行性、适用性，为企业降低物流成本，提高库存周转次数提供了参考依据。

混合整数规划模型(MIP)一般用于解决多设施选址规划问题，具有使用简单且最终解符合实际需求等优点[7]。本文所研究的物流配送中心选址问题是由生产制造部、物流配送中心、地级市批发商组成的三层级供应链模型，如图1所示。其主要分为生产制造部配送至物流配送中心及物流配送中心配送至地级市批发商两阶段，本文研究内容可视为带多约束条件的最优化混合整数规划模型。

在生产制造部和地级市批发商的地理位置以及需求量确定的情况下，这里只考虑物流配送中心的设置，在满足产销量平衡的条件下，将总成本降到最低。因此将物流配送中心选址问题描述为如下目标函数：

minC=Σi=1mΣj=1nΡijRijΤj+Σj=1nΣk=1lQjkRjkΤj+Σj=1nUjWjΤj         (1)" role="presentation" style="padding: 0px; margin: 0px; display: inline; line-height: normal; text-indent: 0px; text-align: left; letter-spacing: normal; word-spacing: normal; overflow-wrap: normal; text-wrap: nowrap; float: none; direction: ltr; max-width: none; max-height: none; min-width: 0px; min-height: 0px; border: 0px; position: relative;">C=Σi=1mΣj=1nPijRijTj+Σj=1nΣk=1lQjkRjkTj+Σj=1nUjWjTj         (1)$?=\underset{?=1}{\overset{?}{{\displaystyle \text{Σ}}}}\underset{?=1}{\overset{?}{{\displaystyle \text{Σ}}}}?{}_{??}?{}_{??}?{}_{?}+\underset{?=1}{\overset{?}{{\displaystyle \text{Σ}}}}\underset{?=1}{\overset{?}{{\displaystyle \text{Σ}}}}?{}_{??}?{}_{??}?{}_{?}+\underset{?=1}{\overset{?}{{\displaystyle \text{Σ}}}}?{}_{?}?{}_{?}?{}_{?}(1)$

其中，C为成本函数，m、n、l、分别表示生产制造部、物流配送中心、地级市批发商的总数，Pij和Rij分别表示第i个生产制造部至第j个物流配送中心的配送量和配送成本，Qjk和Rjk分别表示第j个物流配送中心至第k个地级市批发商的配送量和配送成本，Uj和Wj分别表示第j个物流配送中心的基础建设成本和货物存储能力，Tj为0-1变量，Tj=0表示第j个备选方案未被选中建立物流配送中心，反之Tj=1为被选中。

所述目标函数所需满足的约束条件包括：

a.每个生产制造部制造的快消品数量等于该生产制造部的制造能力：

Σj=1nΡij=Ai,i=1,2,⋯,m" role="presentation" style="padding: 0px; display: inline; line-height: normal; text-indent: 0px; text-align: left; letter-spacing: normal; word-spacing: normal; overflow-wrap: normal; text-wrap: nowrap; float: none; direction: ltr; max-width: none; max-height: none; min-width: 0px; min-height: 0px; border: 0px; position: relative; margin: 0px 30px !important;">Σj=1nPij=Ai,i=1,2,?,m$\underset{?=1}{\overset{?}{{\displaystyle \text{Σ}}}}?{}_{??}=?{}_{?},?=1,2,?,?$ (2)

其中，Ai表示第i个生产制造部的制造能力。

b.每个地级市批发商接收到的快消品总量要满足该批发商的需求量：

Σj=1nQjk=Bk,k=1,2,⋯,l" role="presentation" style="padding: 0px; display: inline; line-height: normal; text-indent: 0px; text-align: left; letter-spacing: normal; word-spacing: normal; overflow-wrap: normal; text-wrap: nowrap; float: none; direction: ltr; max-width: none; max-height: none; min-width: 0px; min-height: 0px; border: 0px; position: relative; margin: 0px 30px !important;">Σj=1nQjk=Bk,k=1,2,?,l$\underset{?=1}{\overset{?}{{\displaystyle \text{Σ}}}}?{}_{??}=?{}_{?},?=1,2,?,?$ (3)

其中，Bk表示第k个地级市批发商的需求总量。

c.每个物流配送中心的快消品进出量要持平：

Σi=1mΡij=Σk=1lQjk,j=1,2,⋯,n" role="presentation" style="padding: 0px; display: inline; line-height: normal; text-indent: 0px; text-align: left; letter-spacing: normal; word-spacing: normal; overflow-wrap: normal; text-wrap: nowrap; float: none; direction: ltr; max-width: none; max-height: none; min-width: 0px; min-height: 0px; border: 0px; position: relative; margin: 0px 30px !important;">Σi=1mPij=Σk=1lQjk,j=1,2,?,n$\underset{?=1}{\overset{?}{{\displaystyle \text{Σ}}}}?{}_{??}=\underset{?=1}{\overset{?}{{\displaystyle \text{Σ}}}}?{}_{??},?=1,2,?,?$ (4)

d.从各生产制造部运至物流配送中心的快消品总量不能超过其建设容量：

Σi=1mΡij≤Wj,j=1,2,⋯,n" role="presentation" style="padding: 0px; display: inline; line-height: normal; text-indent: 0px; text-align: left; letter-spacing: normal; word-spacing: normal; overflow-wrap: normal; text-wrap: nowrap; float: none; direction: ltr; max-width: none; max-height: none; min-width: 0px; min-height: 0px; border: 0px; position: relative; margin: 0px 30px !important;">Σi=1mPij≤Wj,j=1,2,?,n$\underset{?=1}{\overset{?}{{\displaystyle \text{Σ}}}}?{}_{??}\le ?{}_{?},?=1,2,?,?$ (5)

其中，Wj表示第j个物流配送中心的存储能力。

Pij≥0,Qjk≥0,i=1,2,…,m

j=1,2,…,n k=1,2,…,l (6)

e.所有物流配送中心容量之和等于总的需求量：

Σj=1nWj=U,j=1,2,⋯,n" role="presentation" style="padding: 0px; display: inline; line-height: normal; text-indent: 0px; text-align: left; letter-spacing: normal; word-spacing: normal; overflow-wrap: normal; text-wrap: nowrap; float: none; direction: ltr; max-width: none; max-height: none; min-width: 0px; min-height: 0px; border: 0px; position: relative; margin: 0px 30px !important;">Σj=1nWj=U,j=1,2,?,n$\underset{?=1}{\overset{?}{{\displaystyle \text{Σ}}}}?{}_{?}=?,?=1,2,?,?$ (7)

f.所有地级市批发商所需规格集都由物流配送中心整车配送。

布谷鸟搜索算法(Cuckoo Search, CS)是一种启发式算法，最早由著名学者Xin-She Yang与Suash Deb提出，其算法的主要灵感来自于自然界中布谷鸟的繁育方式[8]。其核心思想主要基于两个策略：布谷鸟巢寄生繁育方式以及Lévy飞行机制。

布谷鸟巢寄生繁育方式主要在于宿主鸟巢的选择上，当布谷鸟进入繁殖期，就会去寻找和自己繁殖期、育雏期相近的宿主鸟。在自然进化的发展过程中，布谷鸟鸟蛋的颜色及形状都会和宿主鸟相似，将其下的鸟蛋寄生于其他鸟类的巢穴中，让宿主鸟帮自己孵化鸟蛋。在一般情况下，一个宿主的鸟巢中布谷鸟只会产下一枚鸟蛋，以避免被宿主鸟发现。

在自然界中昆虫或者动物会通过随机或者拟随机的方式进行觅食，他们的飞行轨迹遵循着幂率分布，呈现出Lévy飞行的特征[9]。Lévy飞行是一种看似随机却有规律可循的轨迹形式，步长较大的长距离轨迹和步长较小的短距离轨迹交替出现，如图2所示。在搜索算法中通过Lévy飞行这一方式，在迭代过程中的初期通过大步长来进行巢穴搜索，增加全局的搜索能力，避免出现局部最优且过早收敛的现象[10];在算法迭代的后期，通过较小的步长进行更加精细的搜索从而达到全局最优。

本文采用的是基于物流配送中心编号的编码机制，即从k个候选位置中选择出n个最终地址，这样可以构成编码为：

Z={a1,a2,…,ai,…,aj,…,an},i,j∈{1,2,…,n},i≠j⇒ai≠aj,n≤k (8)

其中，Z表示物流配送中心选址方案。

以上的一组编码表示为一种物流配送中心的选址方案。在算法迭代过程中会根据这样的编码方案进行迭代计算。

首先，通过设置改进算法的参数，随机生成NP个个体，这里个体即为一组编码，因算法采用的是自适应步长因子及发现概率，故此不进行设置。将迭代的最大次数设置为T,作为判断迭代是否结束的终止条件。

进而执行布谷鸟算法产生新解，布谷鸟算法中的新解更新公式为：

C→it+1=C→it+λstep(t)⊕L(δ)" role="presentation" style="padding: 0px; display: inline; line-height: normal; text-indent: 0px; text-align: left; letter-spacing: normal; word-spacing: normal; overflow-wrap: normal; text-wrap: nowrap; float: none; direction: ltr; max-width: none; max-height: none; min-width: 0px; min-height: 0px; border: 0px; position: relative; margin: 0px 30px !important;">C? t+1i=C? ti+λstep(t)⊕L(δ)$\overrightarrow{?}{}_{?}^{?+1}=\overrightarrow{?}{}_{?}^{?}+?{}_{\text{s}\text{t}\text{e}\text{p}}(?)\oplus ?(?)$ (9)

其中，λstep(t)表示步长因子，t表示迭代次数，步长因子的设置与初始化问题的规模有关，C→it" role="presentation" style="padding: 0px; margin: 0px; display: inline; line-height: normal; text-indent: 0px; text-align: left; letter-spacing: normal; word-spacing: normal; overflow-wrap: normal; text-wrap: nowrap; float: none; direction: ltr; max-width: none; max-height: none; min-width: 0px; min-height: 0px; border: 0px; position: relative;">C? ti$\overrightarrow{?}{}_{?}^{?}$表示当前布谷鸟个体寄居的巢穴，当前的寄居巢穴为下一轮迭代寻找新巢穴的关键因素。L(δ)表示为Lévy飞行公式。

再通过粒子群算法对个体产生的新解为Ρ→it+1" role="presentation" style="padding: 0px; margin: 0px; display: inline; line-height: normal; text-indent: 0px; text-align: left; letter-spacing: normal; word-spacing: normal; overflow-wrap: normal; text-wrap: nowrap; float: none; direction: ltr; max-width: none; max-height: none; min-width: 0px; min-height: 0px; border: 0px; position: relative;">P? t+1i$\overrightarrow{?}{}_{?}^{?+1}$,其中粒子群算法中获取新解的计算公式如下：

Ρ→it+1=Ρ→it+φi(Ρ→it-Ρ→kt)" role="presentation" style="padding: 0px; display: inline; line-height: normal; text-indent: 0px; text-align: left; letter-spacing: normal; word-spacing: normal; overflow-wrap: normal; text-wrap: nowrap; float: none; direction: ltr; max-width: none; max-height: none; min-width: 0px; min-height: 0px; border: 0px; position: relative; margin: 0px 30px !important;">P? t+1i=P? ti+φi(P? ti−P? tk)$\overrightarrow{?}{}_{?}^{?+1}=\overrightarrow{?}{}_{?}^{?}+?{}_{?}\left(\overrightarrow{?}{}_{?}^{?}-\overrightarrow{?}{}_{?}^{?}\right)$ (10)

其中，φi代表区间[-1,1]上的随机数，且k≠i,通过粒子群算法产生的新的可能解将与原来的解进行比较，基于贪婪选择策略对解进行排序，并保留较好的解Ρ→it+1" role="presentation" style="padding: 0px; margin: 0px; display: inline; line-height: normal; text-indent: 0px; text-align: left; letter-spacing: normal; word-spacing: normal; overflow-wrap: normal; text-wrap: nowrap; float: none; direction: ltr; max-width: none; max-height: none; min-width: 0px; min-height: 0px; border: 0px; position: relative;">P? t+1i$\overrightarrow{?}{}_{?}^{?+1}$。

通过混合布谷鸟算法和粒子群算法产生的新解表示为Μ→it+1" role="presentation" style="padding: 0px; margin: 0px; display: inline; line-height: normal; text-indent: 0px; text-align: left; letter-spacing: normal; word-spacing: normal; overflow-wrap: normal; text-wrap: nowrap; float: none; direction: ltr; max-width: none; max-height: none; min-width: 0px; min-height: 0px; border: 0px; position: relative;">M? t+1i$\overrightarrow{\text{Μ}}{}_{\text{i}}^{\text{t}+1}$,其公式如下：

Μ→it+1=σ×Ρ→it+1+(1-σ)×C→it+1" role="presentation" style="padding: 0px; display: inline; line-height: normal; text-indent: 0px; text-align: left; letter-spacing: normal; word-spacing: normal; overflow-wrap: normal; text-wrap: nowrap; float: none; direction: ltr; max-width: none; max-height: none; min-width: 0px; min-height: 0px; border: 0px; position: relative; margin: 0px 30px !important;">M? t+1i=σ×P? t+1i+(1−σ)×C? t+1i$\overrightarrow{?}{}_{?}^{?+1}=?\times \overrightarrow{?}{}_{?}^{?+1}+(1-?)\times \overrightarrow{?}{}_{?}^{?+1}$ (11)

其中，σ表示权重系数，σ∈[0,1]。

在传统的布谷鸟搜索算法中，布谷鸟鸟蛋的被发现概率一般设置为常数值0.25,这就易造成在算法迭代过程中会出现因概率固定而导致较优位置和较劣位置被等概率替换的现象。类似地，传统布谷鸟算法的步长因子也通常为常数，因此，算法迭代后期会影响整体的收敛速度及收敛进度，故本文将其设置为自适应的步长因子，其公式分别如下：

θ(t)=θmax-θmax-θminΤ⋅t;" role="presentation" style="padding: 0px; display: inline; line-height: normal; text-indent: 0px; text-align: left; letter-spacing: normal; word-spacing: normal; overflow-wrap: normal; text-wrap: nowrap; float: none; direction: ltr; max-width: none; max-height: none; min-width: 0px; min-height: 0px; border: 0px; position: relative; margin: 0px 30px !important;">θ(t)=θmax−θmax−θminT⋅t;$?(?)=?{}_{\max }-\frac{?{}_{\max }-?{}_{\min }}{?}\cdot ?;$

λstep(t)=λmax-λmax-λminΤ⋅t" role="presentation" style="padding: 0px; display: inline; line-height: normal; text-indent: 0px; text-align: left; letter-spacing: normal; word-spacing: normal; overflow-wrap: normal; text-wrap: nowrap; float: none; direction: ltr; max-width: none; max-height: none; min-width: 0px; min-height: 0px; border: 0px; position: relative; margin: 0px 30px !important;">λstep(t)=λmax−λmax−λminT⋅t$\text{λ}{}_{\text{s}\text{t}\text{e}\text{p}}(?)=\text{λ}{}_{\max }-\frac{\text{λ}{}_{\max }-\text{λ}{}_{\min }}{?}\cdot ?$ (12)

其中，θ(t)表示布谷鸟算法的发现概率，λstep(t)表示布谷鸟算法的步长因子，θmax和θmin分别表示发现概率的上下界限，λmax和λmin分别表示步长因子的最大值和最小值，T表示改进算法的最大迭代次数，t表示当前的迭代次数。

对宿主巢穴实现扰动，并更新解。对于基本的布谷鸟算法来说，主要是通过Lévy飞行来实现寄居巢穴的更新，但是祖先种群与子代种群之间没有关联性，故这里引入适应度权重因子来衡量寄生巢穴的优劣性。针对各地级市批发商的销量构造客户价值权重因子，如下：

Wi=|fbest-fifworst-fbest|" role="presentation" style="padding: 0px; display: inline; line-height: normal; text-indent: 0px; text-align: left; letter-spacing: normal; word-spacing: normal; overflow-wrap: normal; text-wrap: nowrap; float: none; direction: ltr; max-width: none; max-height: none; min-width: 0px; min-height: 0px; border: 0px; position: relative; margin: 0px 30px !important;">Wi=∣∣fbest−fifworst−fbest∣∣$?{}_{?}=\left|\frac{?{}_{????}-?{}_{?}}{?{}_{?????}-?{}_{????}}\right|$ (13)

其中，Wi表示第i个寄居巢穴的适应度权重因子，fi表示第i个布谷鸟巢穴的适应度值，fworst和fbest分别表示布谷鸟寄居巢穴的适应度值的上下界限。

在布谷鸟算法实现随机Lévy飞行算法过程中，根据个体筑巢的关系，增加巢穴之间的适应度权重因子，实现布谷鸟巢穴的更新，最终迭代最优位置公式如下：

S→it+1=S→it+λstep(t)⊕L(δ)⊕(S→it-S→bestt)+η(S→it-S→jt)(Wit-Wjt)         (14)" role="presentation" style="padding: 0px; display: inline; line-height: normal; text-indent: 0px; text-align: left; letter-spacing: normal; word-spacing: normal; overflow-wrap: normal; text-wrap: nowrap; float: none; direction: ltr; max-width: none; max-height: none; min-width: 0px; min-height: 0px; border: 0px; position: relative; margin: 0px 30px !important;">S? t+1i=S? ti+λstep(t)⊕L(δ)⊕(S? ti−S? tbest)+η(S? ti−S? tj)(Wti−Wtj)         (14)$\overrightarrow{?}{}_{?}^{?+1}=\overrightarrow{?}{}_{?}^{?}+?{}_{????}(?)\oplus ?(?)\oplus \left(\overrightarrow{?}{}_{?}^{?}-\overrightarrow{?}{}_{????}^{?}\right)+?\left(\overrightarrow{?}{}_{?}^{?}-\overrightarrow{?}{}_{?}^{?}\right)\left(?{}_{?}^{?}-?{}_{?}^{?}\right)(14)$

其中，λstep(t)表示在执行Lévy飞行过程中第t次的步长因子，L(δ)表示Lévy随机飞行公式，S→bestt" role="presentation" style="padding: 0px; margin: 0px; display: inline; line-height: normal; text-indent: 0px; text-align: left; letter-spacing: normal; word-spacing: normal; overflow-wrap: normal; text-wrap: nowrap; float: none; direction: ltr; max-width: none; max-height: none; min-width: 0px; min-height: 0px; border: 0px; position: relative;">S? tbest$\overrightarrow{?}{}_{\text{b}\text{e}\text{s}\text{t}}^{?}$表示求解过程中在第t次迭代中的最优位置，η为随机数，η∈[0,1],Wit" role="presentation" style="padding: 0px; margin: 0px; display: inline; line-height: normal; text-indent: 0px; text-align: left; letter-spacing: normal; word-spacing: normal; overflow-wrap: normal; text-wrap: nowrap; float: none; direction: ltr; max-width: none; max-height: none; min-width: 0px; min-height: 0px; border: 0px; position: relative;">ti${}_{?}^{?}$和Wjt" role="presentation" style="padding: 0px; margin: 0px; display: inline; line-height: normal; text-indent: 0px; text-align: left; letter-spacing: normal; word-spacing: normal; overflow-wrap: normal; text-wrap: nowrap; float: none; direction: ltr; max-width: none; max-height: none; min-width: 0px; min-height: 0px; border: 0px; position: relative;">tj${}_{?}^{?}$表示第i个、第j个布谷鸟寄生巢穴的适应度权重因子。

根据适应度值来进行优化选择，通过混合前代种群和后代种群来进行合并排序，按照适应度的值从优到劣排序，选择最优的NP个巢穴继续进行后续的迭代。在一轮迭代后，根据终止条件判断结束迭代，终止条件包括迭代结果是否已经收敛、是否达到了迭代次数、迭代得到的最优解与前面得到的最优解平均值的差值是否已经小于设定的阈值等。终止条件已满足，则将该算法的最优解输出，否则继续进行解的迭代更新。在完成迭代后收敛结果，获得物流配送中心选址方案。对最优解中的编码方案进行解码，得到物流中心的选址方案，即从k个候选物流中心方案中选择出n个作为最后的物流配送中心，使得这n个地址能够满足总成本最小、生产制造商和地级市批发商的产销平衡等条件。

根据上述数学模型及算法分析，这里设置算例来评估改进布谷鸟算法在求解物流配送中心的可行性及性能。因产业升级，某大型快消品公司现需要新建物流配送中心为40个地级市批发商提供服务，在满足产销平衡的条件下达到总成本最低，其中新建的物流配送中心从上述的40个地级市中选取。

对于目标函数，假设其所需的约束条件已经满足，此实验案例仅研究改进布谷鸟算法来求解物流配送中心在40个地级市批发商之间的设置方案。生产制造部到物流配送中心、物流配送中心到地级市批发商每单位的运输成本为1元，物流配送中心的单位建造成本为10元；迭代次数设置为200,迭代次数的上限为500。将改进布谷鸟算法的发现概率的上下限分别设置为θmax=0.5,θmin=0.1,步长因子的上下限分别为λmax=0.015和λmin=0.005。

本文所涉及的地级市批发商序号、坐标以及需求量如表1所示，生产制造部位置坐标如表2所示。

将上文所述的参数设置及相关数据代入到改进布谷鸟算法当中，通过Python代码进行算法的实现，最终得到的物流配送中心选址方案如图3所示。最后产生4个物流配送中心为40个地级市批发商提供物流配送服务；物流配送中心同生产制造部之间的关系如图4所示。

此外，为验证在求解物流配送中心选址方案上，改进布谷鸟算法具有一定的性能及优势，将其与传统布谷鸟算法、粒子群算法、人工蜂群算法对于同一实例进行求解对比，结果如表3所示。

由上述对比可知，相较于其他算法，改进布谷鸟算法得到的选址方案总费用在运行时间上有一定优势，能够较快收敛于最优解，提高了求解速度。同时，对于最终的总费用，改进布谷鸟算法比传统布谷鸟算法节省了192,219.41元，节省比例为7.5%;比粒子群算法节省了79,988.5元，节省比例为3.35%;比人工蜂群算法节省了146,008.36元，节省比例为5.8%。此外，在库存周转次数上，改进布谷鸟算法较为显著地提高了库存周转次数，减少了库存积压占用，提高了企业供应链的灵活性。由此可知，改进布谷鸟算法在求解物流配送中心选址问题上是一种有效的算法，并具有一定的可靠性和适用性。

本文研究分析了企业物流配送中心选址问题，提出一种改进布谷鸟搜索算法，基于传统布谷鸟算法在求解问题易早熟、易陷入局部最优这一现状，通过融合粒子群优化算法及引入自适应参数，有效提高了在求解物流配送中心选址问题的效率，降低了求解的时间。通过实例仿真，相较于传统启发式算法，改进的布谷鸟算法在具有较好的全局搜索能力及较短的计算时间的同时，能够较大程度地降低配送成本，优化物流配送效率，提高库存周转次数。