萤火虫算法在供应链库存优化设计中的应用

赵文丹，韩 雪

(沈阳化工大学信息工程学院，辽宁 沈阳 110142)

0 引言

自新冠疫情以来，在国内外产业链不稳定的挑战下，企业间的竞争已经逐步转化为供应链成本间的竞争[1]。从物流角度看，供需不匹配是影响供应链成本的主要原因。库存过多会导致成本增加。当库存过少时，供应链有断链的风险。因此，在物资储存和分配方面，特别是对于库存控制，仍需优化以保证成本最小化，同时维持高水平的客户满意度。对于多级库存控制，较多学者从数理模型、应用算法等角度来求解。李卓群[2]将库存成本作为评价供应链性能的指标，建立了复杂非线性供应链系统动力学模型。在联合库存管理的研究下，Yu W[3]将蚁群算法和模糊模型的基本思想应用于供应链库存优化模型。在需求随机条件下，李沁洪[4]建立了3种模式下的库存协同控制模型。陈圣峻[5]在随机需求下选择供应商和库存协同研究。杨欣[6]利用改进的非支配遗传算法(non-dominated sorting genetic algorithms,NSGA-Ⅱ)解决制造服务组合优化问题。

基于随机需求扰动的前提，赵川[7]从控制理论角度推导出库存传递函数建立库存模型，提出基于自抗扰控制的优化模型。群智能算法在求解工程问题时得到了广泛应用，如粒子群算法[8]、蚁群算法[9]和萤火虫算法(firefly algorithm,FA)[10]等。FA易陷入局部最优而且无法逃离。为了提高FA性能，2018年，L.Zhou[11]等提出了1种混合FA，将正交质心对立学习引入FA，并与基准函数进行了比较。赵嘉[12]对萤火虫种群采用双样本学习来寻找全局最优解。

上述改进算法虽然能保证种群以更加高效的方式搜索空间，但是在面对工程问题时没有摆脱局部最优。基于上述分析，本文从萤火虫物理意义角度改进吸引力项和随机项，既保证群体内相互独立，又采用跳出局部最优解的策略，从而提高算法收敛速度。为了验证本文算法性能，使用3组测试函数，将本文算法与一些最新算法进行比较，并运用到工程领域中，证明算法的实用性。

1 模型描述

1.1 三级供应链混合模型

流程工业下多产品多目标生产线建立了三级混合模型。供应链的仿真模型由多个原材料的供应商、工厂、分销商和不同优先级的零售商组成。为考虑现实生活中的加急订单，特设VIP客户模拟紧急订单。其中1个为带有优先级的零售商，遵从订货优先原则。三级供应链模型如图1所示。

图1 三级供应链模型Fig.1 Three-level supply chain model

如图1所示，供应链以客户下游的订单来拉动整个系统的运作。图1中，实线表示物料流，虚线表示订单流。工厂为有资源限制的节点。工厂一条生产线生产A、B、C这3种产品，每次只加工1种产品。一般操作如下。

由零售商总结顾客1 d之内下订单的总量为累积订单，并在截至日期之前比如第二天清点库存，按照时间顺序来满足客户的订单。如果产品有库存，零售商将打包配送给客户，完成订单；否则，将订单视为待完成的订单。由于待完成的订单未能满足顾客的需要，会导致客户满意度 (customer satisfaction， CS)下降。所以零售商需要向分销商进行订货，从而产生订单来满足下游的需要。如果零售商过多地向分销商下订单会导致运输成本 (transportation cost，TC)增加，并且过多的库存会导致储存成本 (store cost，SC)增加。这一过程在供应链系统中重复出现，一直出现到工厂节点。

制定一个系统的决策使整条供应链的利润最大，也就是在保证成本最小化的同时还要维持较高的CS。因此，本文提出了多产品模型来模拟此类型的问题，并制定相应的策略，用不同算法来比较各种策略的优劣。

1.2 指标

1.2.1 成本指标

本文建立的模型以供应链多级库存成本最小为目标函数，同时维持高水平的CS。

(1)

式中：T为供应链总成本;k为节点，可表示供应商、分销商、零售商的节点;m为客户的数量;P为单位生产成本;S为单位存储成本;A为单位运输成本;W为分销商向供应商订货的数量，即订货量;M为零售商向分销商订货的数量；U为客户向零售商订货的数量;τ为时间常数。

1.2.2 CS

CS为零售商完成客户订单的比值，具体如式(2)、式(3)所示。

(2)

式中：C为客户对于三种产品的满意度;U为客户向零售商节点的订货量;R为客户收到产品的收货量;N为零售商节点个数。

式中：m为客户的数量；λ为三种产品各自占总产品数目的比例。

2 算法的基本原理

2.1 标准FA

在FA中，每只萤火虫个体代表问题的可行解，全局优化的过程实际上是移动每只萤火虫位置的过程。每只萤火虫向着亮度更高的萤火虫飞去，同时亮度正比于萤火虫之间的吸引度。当2只萤火虫距离变远时，由于空气介质的存在，光被吸收，在光减弱的同时吸引力变小。对于1个维度为D，种群个数为N的个体，第i只萤火虫的位置为x1，x2,…,xi。萤火虫之间的亮度和吸引力分别为I和β，萤火虫的相对亮度和吸引力分别如式(4)、式(5)所示。

I=I0×e-γ×rij2

(4)

式中：γ为光强度吸引系数;rij为萤火虫i与萤火虫j之间的欧氏距离；当r=0时，I0为萤火虫的初始荧光亮度。

β=β0×e-γ×rij2

(5)

式中：β0为萤火虫的自我吸引法，由萤火虫之间的亮度决定。

算法的位置更新如式(6)、式(7)所示。

式中：xi、xj为2只萤火虫i、j在维度D中的空间位置。

xi,t+1=xit+β0×e-γ×rij2(xjt-xit)+α×(rand-0.5)

(7)

式中：α为步长因子;rand表示(0，1)之间的随机分布。

2.2 改进的萤火虫算法

2.2.1 萤火虫的吸引力项改进

为了便于研究FA的位置更新方程，将FA分为初始位置项xit、吸引力项Aij和随机项Bij这3项。

Aij=β0×e-γ×rij×rij

(8)

在一维空间中，2只萤火虫xi与xj之间的距离设置为[0，10]。

随着萤火虫之间距离的增加，吸引力项的值从0开始快速增加。当距离大于3时，萤火虫间的吸引力值下降。为了克服这个缺点，在吸引力项中引入了吸引力加强因子，实现全局搜索策略。吸引力项Aij更新如式(10)所示。

式中：K为吸引力加强因子；b、h为待定常数，取值分别为10、100。

式中：γ、β为待定常数，值分别为0.05和1。

2.2.2 萤火虫的随机项改进

当FA进行到后期时，由于萤火虫之间距离逐渐减小，算法趋于局部收敛。引入随机项可以防止FA陷入局部最优，提高FA的搜索效率。为此，本文引入了阻尼振动曲线式(12)。根据振动曲线的优点可知：当萤火虫间距离减小时，式(11)起主要作用；当萤火虫间距离变大时，式(13)起主要作用。重新构建随机模块，如式(13)所示。

Bij=α×(rand-0.5)

(11)

y=A×e-a1×t×sin(a×t)

(12)

Bij=g×k×Ai,j×(rand-0.5)+f×

sin[2π×(rand-0.5)]

(13)

式中：Bij为随机项；k、f、g、a、a1为待定参数，目的是平衡搜索效率和收敛行为。

初值为k=0.8、f=0.89、g=0.62、a=1、a1=1。

2.2.3 自适应光系数

式中：t为当前迭代次数；T为最大迭代次数；s为待定参数，初值设为0.1。

随着迭代次数的增加，萤火虫之间距离减小，介质吸收光变少，改进光吸收系数变为自适应系数，由0.1增大到2。较小的初始值和较快的增加速度，使得算法收敛能力加强。

2.3 莱维飞行与混沌变异

在本文中定义种群N=10，根据适应度取前70%定义为优解，剩下的称为劣解，概率p为(0,1)均匀的随机数。将优解和依据概率p选择的劣解共同执行改进后的萤火虫飞行策略，将(1-p)概率的劣解引入莱维飞行策略，以保证种群的多样性。其位置迭代式为：

Xi,t+1=Xit+β0×e-γ×rij2×(Xjt-Xit)+α×

sign(rand-0.5)⊗levy(β)

(15)

式中：Xit为t时刻萤火虫所在位置;α为步长控制因子，本文取值为 1;sign()为符号函数；⊗为内积运算。

式中：v和u是服从方差为1、均值为0的标准正态分布的随机变量。

式中：β为常数，本文取β=1.5；Γ为伽马函数。

在算法迭代过程中，对萤火虫局部最优解进行检测。若某个体连续5代未发生变化，说明个体达到局部最优，易陷入局部最优解，则执行以下混沌变异过程。

式中：d为步长;xmax和xmin分别为搜索空间范围的边界值。

xit=xmin+μ×(xmax-xmin)×d×(1-d)

(19)

式中：μ为混沌控制参量。

本文中μ=4，而且种群中的每个粒子的位置都有(1-p)的概率进行混沌变异。这进一步增强了全局的搜索能力，扩大了搜索空间和种群的多样性。

2.4 改进萤火虫算法的流程

改进萤火虫算法(improve firefly algorithm,IFA)流程如图2所示。

图2 IFA流程图Fig.2 IFA flowchart

3 算例分析

3.1 基准函数测试

仿真试验环境是Matlab2018b，采取3个基准测试函数来验证IFA的优化性能，与FA进行对比。其中F1、F2为单峰函数，F3为多峰函数，取种群数目N=25，迭代次数T=1 000，维度D=10，运行20次。其中，基本测试函数F1为：

式中：x在F1中范围为[-10,10]。

测试函数F2为：

式中：x在F2中范围为[-30,30]。

测试函数F3为：

(22)

式中：x在F3中范围为[-32,32]。

FA与IFA的最优适应度与平均适应度对比如表1所示。

表1 FA与IFA的最优适应度与平均适应度对比Tab.1 Optimal fitness and average fitness compared with FA and IFA

由表1可知，IFA比FA具有更好的适应度值、更优的平均值。由此说明改进后的算法真实有效。

以测试函数F1为例，算法收敛速度对比如图3所示。

图3 算法收敛速度对比图Fig.3 Comparison of algorithm convergence speed

由图3可知，IFA 在70代收敛；FA在90代收敛；和声搜索(harmony search,HS)算法在120代收敛；人工蜂群(artifical bee colomy,ABC)算法在130代收敛。由此可以看出，IFA收敛速度最快。迭代后期可以看出:虽然HS算法、ABC算法、FA都已收敛，但是精度不高;IFA具有更快的速度收敛而且精度最高。

3.2 工程算例

在Matlab2018b中，用正弦信号模拟客户3种产品的订单，按照控制率零售商向上级订货。控制率为：

U1=λ×(R-P)

(23)

式中：U1为零售商节点向分销商节点订货的件数；P为零售商库存；λ为零售商向上级订货的比例。

式中：U为客户订货量；K1、K2为λ对应的参数。

按照控制率零售商向分销商订货，分销商向供应商订货，计算供应链总成本，用IFA进行例证，用ABC算法与HS算法进行对照仿真试验。不同算法成本对比如图4所示。

图4 不同算法成本对比图Fig.4 Cost comparison of different algorithms

将IFA用于实际供应链三级库存模型中，HS算法易陷入局部最优解，ABC算法没有达到最优解还降低了CS。各算法对应指标如表2所示。由表2可知，IFA能够避免陷入局部最优解，同时成本减少10.1%，保持CS 98%以上。

表2 各算法对应指标Tab.2 Corresponding indexes of each algorithm

4 结论

本文提出了IFA，引入吸引力增强因子和自适应光系数来增加算法全局搜索速度。针对FA在寻优过程中容易陷入局部最优解的问题，采用莱维飞行和混沌变异，在增加种群多样性的同时避免局部最优。数值分析证明，改进后的算法真实有效，应用于供应链库存模型中可减少库存成本98.9万元、提高CS至98.4%。该结果证明了IFA在供应链库存控制上是更优的解决方案 。