基于多精度模型和粒子群算法的(s,S)库存优化

朱晨波,王 星,王伍娜

(浙江工业大学 管理学院,浙江 杭州 310023)

库存决策一直是企业管理决策中的重要一环,各类企业都面临着复杂的库存决策问题[1-4]。自从Arrow等[5]于20世纪50年代提出了多阶段定期盘点的库存模型以来,(s,S)库存系统的研究一直为学者们所关注。很多文献证明了(s,S)库存策略的最优性,比如,在各阶段的需求是独立同分布的情况下,Scarf[6]和Iglehart[7]证明了(s,S)库存策略的最优性;在各阶段的需求是时序相关的情况下,Sethi等[8]用马氏需求来构建需求的时序相关性,以最小化总库存成本为目标,不仅证明了(s,S)策略对于需求时序相关的、具有固定订货成本的库存系统而言是最优的,而且证明了(s,S)策略在库存容量和服务水平约束、非平稳无限期的系统等情况下也是最优的。

在(s,S)库存策略的研究领域,优化库存策略参数(计算最优的s和S)一直是一个重要的课题。针对需求时序相关的(s,S)库存系统的研究,由于其系统较为复杂且库存策略参数较多,一般需要设计启发式算法优化该系统的库存策略。笔者为一类需求时序相关的(s,S)库存系统的策略优化问题设计了改进型粒子群优化算法,通过引入多精度模型,并基于次序变换和优化抽样的多精度优化方法,提高了传统的粒子群优化算法在计算最优的库存策略参数时的计算效率。

1 研究背景

研究库存策略参数优化的文献主要可以分为基于需求、成本特性设计的优化算法和一般的启发式算法两个方面。在基于需求、成本特性设计的优化算法方面,一般而言,其所研究的库存系统的需求、成本都具有较多特定的结构和假设。Liu等[9]和Liu等[10]所研究的库存系统都假设需求的到达是一个泊松过程,据此分别推导出库存控制策略参数的最优解析解。Song等[11]研究的库存系统的需求是一个马尔可夫驱动的泊松过程,假设其相关的成本函数是线性的,在此假设下设计了库存控制策略参数的优化算法。Walid等[12]研究了与Song等[11]类似的问题,针对需求为马尔可夫驱动的泊松过程,提出了更有效的算法来优化库存控制策略的参数。Feng等[13]研究了多产品、需求服从泊松分布的连续盘点的库存系统,首先用马氏决策过程的方法计算最优策略,但仅限于小规模的问题;然后根据最优策略的结构,设计了启发式策略;最后提供了计算该策略参数的算法。在使用一般的启发式算法方面,Tsou[14]使用类电磁机制算法和粒子群优化算法来优化多目标的连续盘点的库存模型。Taleizadeh等[15]使用模糊仿真和遗传算法来优化具有模糊需求和多产品的库存系统。Saracoglu等[16]使用遗传算法优化多产品多阶段的(Q,r)库存模型。与基于需求、成本特性设计的优化算法相比,一般的启发式算法适用于优化假设较少、结构较复杂和一般化的库存系统,其缺点是因为没有考虑目标系统的结构特性以及通常需要结合仿真方法,从而导致算法效率较低。

近年来,基于次序变换和优化抽样的多精度优化方法(Multi-fidelity optimization with ordinal transformation and optimal sampling,MO2TOS)常被用来提升启发式算法的计算效率。Xu等[17]首次提出了MO2TOS方法,该方法引入了次序变换和优化抽样的思想,利用低精度模型对原有的高精度模型进行辅助优化。之后,MO2TOS方法被用于研究多种类型的问题,Zhao等[18]在随机客户订单的调度问题研究中,为了提高候选解的评估效率,建造了一个低精度模型来选择性能较好的解进行高精度模型仿真,利用低精度模型的特点,开发了一种有次序变换与优化抽样的仿真优化算法。郝春锋等[19]应用MO2TOS算法框架求解包含很多不确定性的MRO调度问题。Chiuc等[20]将次序变换的思想与进化算法相结合,提出了一种新的算法,以提高大规模确定性优化问题的计算效率。Li等[21]首次在元优化系统中使用多精度策略,提出了一种基于多精度策略的元优化方法来加速参数优化。Shu等[22]提出了一种多精度贝叶斯优化方法,该方法采用分层克里格法构建多精度元模型,将标准的贝叶斯方法扩展到多精度优化,利用低精度模型的信息进一步降低优化成本,结果表明该方法能够以较低的计算成本获得全局最优解。

综上所述,针对结构较复杂、一般化的库存系统,通常结合仿真方法和启发式算法来优化库存策略,但算法的效率较低。笔者的研究将多精度模型以及MO2TOS方法整合到粒子群优化算法中,以提升该算法的效率,并用来优化需求时序相关的(s,S)库存策略。

2 (s,S)库存系统的多精度模型描述

笔者研究由马尔可夫过程驱动的平稳需求的(s,S)库存系统,假设不存在订货提前期,且每个周期的缺货是需要完全回补的,库存系统的成本包括固定和变动的订货成本、库存持有成本以及缺货成本。假设Dn为周期n的需求,而{Jn:n≥1}为驱动{Dn:n≥0}的马尔可夫链,且{Jn:n≥1}在需求状态空间Ω上不可约并且遍历,其转移概率矩阵为P(pij)(i,j∈Ω),稳态概率为πj(j∈Ω)。进一步地,假设{Jn:n≥1}是平稳的,即P(Jn=j)=πj(n=1,2,…),假设(s,S)库存策略依赖于需求状态,令(sj,Sj)为与需求状态j∈Ω相对应的参数,并且qj=Sj-sj。令Xn为周期n期初的库存水平,根据(s,S)库存策略,如果库存水平低于sJn,那么需要下订单使库存水平回到SJn,否则不需要下订单,库存水平将减去当期需求从而继续下降。

当需求服从位相型分布(如指数分布、爱尔朗分布等)时,Hu等[23]提出了两种数值方法可以快速计算上述库存系统的库存水平和库存成本。但是,当需求服从正态分布或截断的正态分布时,上述数值方法不再适用,只能使用仿真方法计算库存水平和库存成本,这将耗费大量的计算时间。把需求服从指数分布的库存系统称为低精度模型,把需求服从截断的正态分布的库存系统称为高精度模型。基于500个备选解使用数值方法计算低精度模型的库存成本的值,按照库存成本从小到大的次序为这些备选解排序,具体如图1的虚线所示。运用仿真方法计算低精度模型每一个备选解所对应的高精度模型的库存成本,并一一对应地展现在同一张图上,具体如图1的实线所示。

图1 次序变换后低精度与高精度模型的库存成本对比Fig.1 The comparison ofinventory costs between low-fidelity and high-fidelity models after ordinal transformation

由图1可知:与高精度模型相比,低精度模型虽然库存成本值有很大的误差,但是计算出来的备选解的库存成本大小次序基本上与高精度模型一致。针对低精度模型,与库存水平相关的成本指标如下：

1) 阶段n订货量的表达式为

(1)

2) 阶段n是否订货的表达式为

(2)

3) 阶段n的期初库存量表达式为

Vn=max{Xn,0}

(3)

4) 阶段n的期初缺货量表达式为

Wn=max{-Xn,0}

(4)

G({sj},{qj};j∈Ω)=KE[δ(U)]+
cE[U]+hE[V]+rE[W]

(5)

该长期平均成本表达式可以用库存水平的各阶矩的线性组合来表示,对于服从位相型分布的需求,库存水平的各阶矩可通过Hu等[23]所提出来的数值方法快速精确地计算得到,因此低精度模型的目标函数值(长期平均库存成本)也可以被很快地计算出来。需要指出的是,因为可变采购成本的值是不变的,所以进行分析的成本为

G({sj},{qj};j∈Ω)=
KE[δ(U)]+hE[V]+rE[W]

(6)

3 基本算法

3.1 粒子群优化算法

粒子群优化(Partical swarm optimization,PSO)算法是Kennedy等[24]于1995年提出来的,是一种典型的群体智能优化算法。之后,很多学者对该算法进行了各种不同的改进,如Ratnaweera等[25]、Li[26]、Li等[27]、胡旺等[28]、邱飞岳等[29]以及李浩君等[30]。在该算法中,粒子具有速度和位置两个属性。每个粒子都被看作是D维空间中的一个点,也是目标函数的一个候选解。第i个粒子的位置表示为xi=(xi1,xi2,…,xiD),其中xid∈[xmin,xmax],d=1,2,…,D,xmin和xmax分别为xid所能取到的最小值和最大值。第i个粒子的速度表示为vi=(vi1,vi2,…,viD),其中vid∈[vmin,vmax],d=1,2,…,D,vmin和vmax分别为速度的最小值和最大值。每个粒子的适应度值由预定义的研究问题的目标函数值所决定,每个粒子会记忆在整个迭代过程中自身曾找到的最优位置(即个体最优解),把第i个粒子所找到的最优位置记为pbesti,把整个种群所发现的全局最佳位置(全局最优解)记为gbest(gbest是所有pbesti中的最佳值)。在每次迭代中,每个粒子会更新自己的速度和位置,其更新迭代表达式为

vid(t+1)=w×vid(t)+c1×rand()×(pbestid(t)-
xid(t))+c2×rand()×(gbestd(t)-xid(t))

(7)

xid(t+1)=xid(t)+vid(t+1)

(8)

式中:i为第i个粒子,i=1,2,…,N,N为粒子群中的粒子总数;d为粒子的第d个维度,d=1,2,…,D,D为搜索空间的总维度;t为当前的迭代次数;c1和c2为学习因子;rand()为(0,1)之间的随机数;vid(t)为第t次迭代,第i个粒子在第d维上的速度;xid(t)为第t次迭代,第i个粒子在第d维上的位置;pbestid(t)为到第t次迭代为止第i个粒子在第d维上的最优位置;gbestd(t)为到第t次迭代为止整个种群找到的在第d维上的最优位置;w为惯性权重,用于调节对解空间的搜索范围。对w进行迭代,计算式为

w=wmin+(wmax-wmin)×(T-t)/T

(9)

式中:wmin为w的最小值;wmax为w的最大值;T为最大迭代次数;t为当前迭代次数。同时考虑精度与最大迭代次数来为粒子群算法设定终止条件,其中精度的计算式为

(10)

式中:yi为在第i次迭代时得到的最优适应度值;yi-1为在第i-1次迭代时得到的最优适应度值,当连续一定的次数,精度都小于某个阈值,那么迭代停止。另外,设定最大迭代次数,即当程序迭代到了最大次数就停止运行。

3.2 MO2TOS方法

MO2TOS方法由次序变换(Ordinal transformation,OT)和优化抽样(Optimal sampling,OS)两部分组成。OT方法首先利用低精度模型对所有的解进行评估得到每个解的适应度值；然后根据适应度值的大小进行从小到大排序,将原解空间转换为一维有序空间。OS由某种计算量分配方法,确定每个组应被分配到的抽样数目,反复使用OS,直到计算量预算用完为止。MO2TOS方法流程如图2所示。

图2 MO2TOS方法流程图Fig.2 The flow chart of the MO2TOS method

用MO2TOS方法挑选粒子群时,在次序变换阶段,首先基于所有的备选解计算低精度模型的适应度值,按照适应度值从小到大的次序为这些备选解排序(以库存成本最小为目标,故适应度值越小越好)；其次根据排序把这些备选解分成k组,并从每组当中抽样出n0个备选解,从而得到kn0个初始的备选解；再次进入优化抽样阶段,只要备选解的数量小于N(N为设定的粒子群中的粒子数量),就增加Δ个抽样数；最后根据一定的优化抽样规则确定每个组应被分配到的抽样数的增加数目,直至备选解(粒子)的数量达到N为止。其中,笔者采用的计算量分配方法是Chen等[31]在2000年提出的最优计算量分配(Optimal computing budget allocation,OCBA)方法。

4 基于粒子群算法的(s,S)库存策略优化

针对需求时序相关的(s,S)库存系统,提出了3种基于粒子群优化的算法对库存策略参数进行优化,分别为传统粒子群优化算法(算法1)、考虑多精度模型的粒子群优化算法(算法2)以及考虑MO2TOS和多精度模型的粒子群优化算法(算法3)。

4.1 传统粒子群优化算法

在传统粒子群优化算法(算法1)中仅使用高精度模型来计算粒子的适应度值,算法1的实现步骤为

初始化步骤:随机地初始化N个粒子,令t=0,评估每个粒子的适应度值,确定pbest,pb,gbest,gb。其中pb为每个粒子自身的最优适应度值(评估pbest得到的最优适应度值);gb为整个粒子群的全局最优适应度值(评估gbest得到的最优适应度值)。

迭代步骤为

步骤1根据式(9)计算w。

步骤2根据式(7,8)调整粒子的速度和位置,更新pbest,pb,gbest,gb。

步骤3根据式(10)计算ε。如果ε≤precision,那么criterion1=criterion1+1;否则criterion1=0。其中precision为给定的精度数值；criterion为精度连续在某一阈值下的最大迭代次数；criterion1为精度连续在某一阈值下的累计迭代次数。

步骤4如果criterion1≥criterion或者t>maxgen,停止迭代,输出返回值gbest,gb;否则,令t=t+1,并转到步骤1。其中maxgen为最大迭代次数。

4.2 考虑多精度模型的粒子群优化算法

在多精度模型的粒子群优化算法(算法2)中,首先基于低精度模型进行优化,然后基于高精度模型进行优化。结合低精度与高精度的多精度优化方法可以使算法加速收敛到高精度模型下的最优解,从而以较低的计算成本获得全局最优解。具体地,首先,基于低精度模型使用粒子群优化算法对库存系统进行优化,其优化过程与算法1相同,并得到一组pbest和gbest;然后,基于高精度模型计算低精度模型得到这些粒子的适应度值,从而得到另一组更新了适应度值的“初始粒子群”;最后,再次使用粒子群优化算法对高精度模型进行优化,并得到最终的最优解。算法2的实现步骤为

1) 第1次初始化步骤。随机初始化N个带有速度和位置属性的粒子,令t=0,用低精度模型评估每个粒子的适应度值,确定pbest,pb,gbest,gb。

2) 第1次迭代步骤。使用算法1的迭代步骤。

3) 第2次初始化步骤。基于高精度模型计算上述迭代步骤结束后得到粒子群的pbest和gbest的适应度值,重新得到一组更新了适应度值的“初始粒子群”。

4) 第2次迭代步骤。使用算法1的迭代步骤。

4.3 考虑MO2TOS和多精度模型的粒子群优化算法

在算法1和算法2中,初始粒子群的生成是随机的。在考虑MO2TOS和多精度模型的粒子群优化的算法(算法3)中,使用MO2TOS方法对初始粒子进行择优选择。

算法3的实现步骤为

1) 第1次初始化步骤。

步骤1随机初始化M个仅有位置属性的粒子,令t=0,i=0,用低精度模型评估每个粒子,得到适应度值。其中M为进行OT的粒子总数。

步骤2根据适应度值对M个初始化粒子进行从小到大排序,均等分为k个组,并从每组中随机选取n0个粒子。

步骤3用高精度模型评估被选择粒子的适应度值,并计算出每组当前选中的粒子的适应度值的均值和方差。

步骤4根据OCBA方法,计算分配给每组的粒子数Nj,并在第j组中随机选取Nj个粒子(j=1,2,…,k)。

步骤5i=i+1,如果i<(N-n0*k)/b,转到步骤3,否则,通过MO2TOS方法得到N个仅有位置属性的粒子,并初始化每个粒子的速度,通过低精度模型评估粒子的适应度值,确定pbest,pb,gbest,gb,进入第一次迭代步骤。其中b为每一次分配的粒子数。

2) 第1次迭代步骤。使用算法1的迭代步骤。

3) 第2次初始化步骤。基于高精度模型计算上述迭代步骤结束后得到粒子群的pbest和gbest的适应度值,重新得到一组更新了适应度值的“初始粒子群”。

4) 第2次迭代步骤。使用算法1的迭代步骤。

5 数值算例

设计算例来测试笔者所提的算法,算法的编程使用C++语言,并在1.60～2.11 GHz、Intel i5-10210U处理器、16 G的笔记本电脑上测试。设置了12种不同的参数组合进行实验,包括两组不同的概率转移矩阵P1,P2和稳态概率π1,π2,以及两组不同的指数分布参数λ1,λ2,具体赋值为

另外,令N=50,D=5,maxgen=100,k=5,b=10,ε=0.002,criterion=5,c1=2,c2=2,n0=5。在算例中,高精度模型的适应度值是经过1 000万次仿真得到的,运行一次大约需要5～6 s,而低精度模型的适应度值是使用Hu等[23]所提出的数值算法计算得到的,运行一次平均只需要0.08 s。数值算例的结果如表1～4所示。

表1 3种算法的结果比较(π1,P1,λ1)Table 1 The comparison of results of the three algorithms (π1, P1, λ1)

表2 3种算法的结果比较(π1,P1,λ2)Table 2 The comparison of results of the three algorithms (π1, P1, λ2)

表3 3种算法的结果比较(π2,P2,λ1)Table 3 The comparison of results of the three algorithms (π2, P2, λ1)

表4 3种算法的结果比较(π2,P2,λ2)Table 3 The comparison of results of the three algorithms (π2, P2, λ2)

由表1～4可知:无论在哪种参数组合下,3种算法得到的最低库存成本都很接近,最终都能收敛到最优解。在算法的迭代次数上,算法1的迭代次数是算法2或算法3的2到3倍,算法2和算法3的计算效率均比算法1要高。但是,算法2和算法3的迭代次数在不同的情况下互有优劣,总体上不分伯仲,可见在粒子群算法的初始化阶段引入MO2TOS方法去改进初始粒子群的质量并不能提升算法2的计算效率,说明粒子群算法初始解的优劣选择对于算法的收敛速度没有显著影响。

对(π1,P1,λ1)和(π1,P1,λ2)两组参数设置下的算法收敛情况进行绘图,如图3,4所示。由图3,4可知:算法2和算法3收敛速度很快,迭代5次或者更少就能找到最优解;而算法1的收敛速度就慢一些。

图3 3种算法的收敛情况比较(π1,P1,λ1)Fig.3 The comparison of the convergency of the three algorithms (π1, P1, λ1)

图4 3种算法的收敛情况比较(π1,P1,λ2)Fig.4 The comparison of the convergency of the three algorithms (π1, P1, λ2)

6 结 论

针对一类需求时序相关的(s,S)库存系统的策略优化问题,把需求服从截断的正态分布的库存系统称为高精度模型,把需求服从指数分布的库存系统称为低精度模型,其中高精度模型为研究的目标模型,低精度模型为优化计算时所用的辅助模型,提出了两种改进型粒子群优化算法。其中一种算法将多精度模型和粒子群优化算法相结合,形成混合粒子群优化算法;另一种算法在前一种算法的基础上还在粒子群初始化过程中引入了MO2TOS方法。数值算例结果表明:相较于传统的粒子群优化算法,笔者提出的改进型粒子群优化算法具有明显更高的计算效率。