文章快速检索    
  同济大学学报(自然科学版)  2019, Vol. 47 Issue (12): 1831-1836.  DOI: 10.11908/j.issn.0253-374x.2019.12.020
0

引用本文  

牟刚, 孙永超, 徐承龙. 基于改进树搜索方法的库存系统优化[J]. 同济大学学报(自然科学版), 2019, 47(12): 1831-1836.   DOI: 10.11908/j.issn.0253-374x.2019.12.020
MU Gang, SUN Yongchao, XU Chenglong. Inventory System Optimization Based on Improved Tree Search Method[J]. Journal of Tongji University (Natural Science), 2019, 47(12): 1831-1836.   DOI: 10.11908/j.issn.0253-374x.2019.12.020

基金项目

中央高校基本科研业务费

第一作者

牟刚(1978—),男,博士生,主要研究方向为机器学习与随机优化.E-mail: mug@tongji.edu.cn

通信作者

徐承龙(1963—),男,教授,博士生导师,理学博士,主要研究方向为金融工程与风险管理.E-mail: clxu@tongji.edu.cn

文章历史

收稿日期:2018-12-10
基于改进树搜索方法的库存系统优化
牟刚 1, 孙永超 1, 徐承龙 2     
1. 同济大学 数学科学学院,上海 200092;
2. 上海财经大学 数学学院,上海 200433
摘要:提出一种算法对库存系统的库存策略进行优化,保障库存系统服务水平达到一定水平的同时最小化库存系统成本.将库存系统优化问题抽象为一个随机优化问题,结合克里金插值和蒙特卡罗树搜索求解这一随机优化问题,提高运算效率.将算法应用于真实算例中取得了很好的效果.
关键词供应链管理    库存系统    克里金插值    蒙特卡罗树搜索    
Inventory System Optimization Based on Improved Tree Search Method
MU Gang 1, SUN Yongchao 1, XU Chenglong 2     
1. School of Mathematical Sciences, Tongji University, Shanghai 200092, China;
2. School of Mathematics, Shanghai University of Finance and Economics, Shanghai 200433, China
Abstract: Propose an algorithm to optimize the strategy for an inventory system, minimize the storage cost with constraint of satisfying service level. Abstract the inventory system optimization problem into a stochastic optimization problem, combine Kriging regression and Monte Carlo tree search to solve the stochastic optimization problem efficiently. Apply this method into a real world example and get positive feedback.
Key words: supply chain management    inventory optimization    Kriging interpolation    Monte Carlo tree search    

库存系统优化是供应链管理中的一个重大问题[1].主要考虑两方面问题:一方面过多的库存积压会给企业带来巨额的资金占用成本;另一方面若出现库存短缺的情况,会降低企业的服务水平.库存系统优化希望通过优化方法找到最优库存策略,在保证整个系统服务水平的基础上,最小化总体成本.

一般用随机优化的模型表示库存系统优化问题,用模拟方法进行求解.对于一个全球的供应链的实际问题,需要考虑的因素比较多,供应链节点数量大大增加,优化问题的计算量也相应增加.目前的解决方法是使用随机切平面方法、遗传算法或者粒子群优化算法等来进行优化.但是对于高维问题,这些算法的速度都很慢,很难满足企业的实际需要.

对每个仓库来说,基本的库存策略有(r, Q)策略,(r, S)策略等[2].(r, Q)策略指对库存进行连续性检查,当库存小于订货点水平r时,通过向上一级仓库订货进行库存补充,每次订货的数量为Q[3];(r, S)策略指对库存进行连续性检查,当库存小于订货点水平r时,通过向上一级仓库订货进行库存补充,控制订货量使得在订货后,库存总量为S[4].本文的算法过程和算例都是在(r, Q)策略下进行的,也可以应用于其他策略.

在给定的库存策略下,库存系统仍然存在多个随机变量,包括每个仓库每天收到的订单量以及货物在两个仓库间交货时间等,因此库存系统的成本和服务水平是有一定随机性的黑箱函数,最终优化问题的目标函数和限制条件是随机黑箱函数的期望值,用多次的场景模拟的平均值对其期望进行估计,问题的维度取决于仓库的数量.利用克里金插值可以显著减少运算量,克里金插值最早被应用在统计地质学[5].Ankenman等[6]在传统克里金方法的基础上,考虑了样本点处取值的随机波动,构造了随机克里金方法.克里金方法还可以用于求解优化问题,Jones等[7]提出的高效全局优化方法(efficient global optimization, EGO).利用克里金方法,通过每次求解预期改进(expectation improvement, EI)函数最大的问题,迭代求解全局最优问题;Huang等[8]在Jones等[7]的基础上,使用了不同的预期改进函数,使求解更快.

蒙特卡罗树搜索是一种在给定的领域中通过在决策空间进行随机抽样,从而生成搜索树来寻找最佳决策的方法,广泛应用于人工智能领域中.其中最为代表的是阿尔法围棋程序[9]在5局3胜制的比赛中以4:1赢得了韩国棋手李世石九段.本文对蒙特卡罗树搜索方法进行了改进,并和克里金插值方法相结合,用于求解多级库存系统的优化问题.在搜索空间内利用克里金插值,得到目标函数在整个搜索空间的插值函数形式,然后将搜索空间分成多个子空间并通过蒙特卡罗树搜索选择最优的子空间,重复进行这一过程,最终求得最优解.本文将蒙特卡罗树搜索方法与克里金方法相结合解决黑箱函数的优化问题.

1 问题描述 1.1 多级库存系统

分布式多级库存系统是一个树形结构,每个节点只有一个上级节点,但是可能有多个下级节点,没有上级节点的称为根节点,没有下级节点的称为叶节点,每个节点表示一个仓库.假设根节点仓库总保持有充足库存.除根节点仓库外,每天每个仓库都会收到客户的订单,订单量服从某一随机分布.

图 1所示的是本文考虑的分布式库存系统的一个案例结构,其中U(a, b)表示(a, b)间的均匀分布, 除根节点仓库外共有4个仓库,编号分别为1,2,3,4.假设每个仓库每天收到的订单量分别服从正态分布N(1 000, 500),N(600, 400),N(500, 300),N(400, 300).4个仓库都使用(r, Q)库存策略.每天每个仓库都会检查库存水平,当第i个仓库的库存小于ri时,会向上级仓库发出订货订单,订货量为Qi.货物从上级节点运送到下级节点的运输时间服从均匀分布U(10, 15),U(4, 8),U(5, 10),U(4, 10).策略下每个仓库的库存水平变化方式如图 2所示.

图 1 多级库存系统的结构 Fig.1 The structure of multi-echelon inventory system
图 2 (r, Q)策略下库存水平变化情况 Fig.2 Dynamic of inventory system under (r, Q) strategy

图 2中,现有库存(on-hand inventory, OHI, 以SOHI表示)指仓库中现有货物量,缺货(Backorder, 以SBO表示)指当现有库存不能满足订单量而不能发货的数量,在途库存(on-order inventory, OI,以SOI表示)指已经向上级仓库订货而尚未到达的货物量,库存状况(inventory position, IP,以SIP表示)定义为现有库存减缺货量加在途库存,即:SIP=SOHI-SBO+SOI

发货时间(leading time, LD, 以TLD表示)指从向上级仓库下订单到该笔订单送达的时间,发货时间有如下表达式:

$ {T_{{\rm{LD}}}} = \left\{ \begin{array}{l} {T_{{\rm{PT}}}}\;\;\;\;\;\;\;\;\;\;上级仓库有足够的{\mathit{S}_{{\rm{OHI}}}}\\ {T_{{\rm{PT}}}} + T_{{\rm{LD}}}^\prime \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;上级仓库没有足够的{\mathit{S}_{{\rm{OHI}}}} \end{array} \right. $

其中:T′LD为上级仓库的发货时间; 运输时间(processing time, PT,以TPT表示)是货物在两个节点间运输的时间.当上级节点的SOHI大于下级节点的订单量时可以立即发货,TLD就是TPT;否则,需要等待上级节点补充库存.

在求解最优库存策略时,目标函数为某段时间内的总成本,限制条件为这段时间内的服务水平大于某一常数.总成本分为三部分,每天每单位在运输中将要运输至第i个节点的货物产生COi的在途成本,每天每单位在第i个节点仓库中的货物产生CHi的储存成本,每天第i个节点每单位的缺货量产生CSi的缺货成本;同时用及时供应率POTIF描述服务水平,为能立即发货的订单总量除以总订单量.

1.2 样本均值估计

x =(r1, r2, r3, r4, Q1, Q2, Q3, Q4), 库存系统的输入包含客户订单量和订单的运输时间等随机变量,记这些随机变量为w,总成本为C(x, w),及时供应率为F(x, w), 供应率的限制为Fmin.需要求解的优化问题为

$ \begin{array}{l} \;\;\mathop {\min }\limits_\mathit{\boldsymbol{x}} \phi (\mathit{\boldsymbol{x}}) = {E_w}[C(\mathit{\boldsymbol{x}}, w)]\\ {\rm{s}}{\rm{.t}}{\rm{.}}\;\psi (\mathit{\boldsymbol{x}}) = {E_w}[F(\mathit{\boldsymbol{x}}, w)] \cdots {F_{\min }} \end{array} $

这个随机优化问题的目标函数和限制条件都是某个随机变量的期望,且这个期望没有显式表达式,因此需要利用n次场景模拟下C(x, wi)和F(x, wi)的取值的均值ϕn(x)和ψn(x)代替期望,其中, i=1, …, n.这种方法成为样本均值估计(sample average approximation, SAA),Kim等[10]证明了利用这种方法求解优化问题时的收敛性和收敛速度.

1.3 基于代理的模拟

对于给定x,某一个场景wi下的成本和服务水平C(x, wi)和F(x, wi)仍然没有显示表达式,这是因为无法直接通过xwi直接得到观测时间内节点间的订单和运输情况,必须对于整个库存系统的运行过程进行模拟,才能根据库存策略判断订单和运输的产生,进而得到成本和服务水平.基于代理的模拟可以方便地模拟库存系统的动态过程,Chu等[3],Rahmandad等[11]都使用这一方法对库存系统进行模拟.

本文利用客户代理,订单代理,运输代理,仓库代理模拟整个库存系统的运行.每天客户代理产生随机的订单量,订单代理将信息从下级节点传到上级节点,运输代理将货物从上级节点传到下级节点,仓库代理根据订单代理和运输代理更新库存状况,并根据库存策略判断是否需要产生新的订单和运输.在观测时间0到T内循环这一过程,通过这4类代理的协同作用,模拟整个库存系统的运行情况,并最终返回所需的两个函数C(x, wi)和F(x, wi)的数值.

由于需要在0到T时间内循环模拟,同时对每个x需要进行n次这样的场景模拟,因此为得到某一库存策略x对应的期望成本和服务水平,需要与nT成正比的运算量.为了减小运算量,使用克里金插值.

2 克里金插值

克里金插值法广泛应用于统计地质学中,也被用于全局优化算法中[12].最普通的克里金方法假设,响应面也就是要估计的函数有如下形式:

$ Y(\mathit{\boldsymbol{x}}) = f{(\mathit{\boldsymbol{x}})^{\rm{T}}}\mathit{\boldsymbol{\beta }} + M(\mathit{\boldsymbol{x}})$

其中,M是从某一高斯随机场的一次抽样,高斯随机场满足条件E[M(X)]=0,E[M2(x)] < ∞,且当xx′的取值相近时,M(x)和M(x′)的取值也趋于相近,M(x)和M(x′)的协方差可以写成τ2R(|x - x′|; θ),其中τ2为随机场的方差,R(|x - x′|; θ)为xx′之间的变异系数,θ是这个函数的参数,变异系数函数需要满足R(0;θ)=1,且当x和x′的距离趋于无穷时R(|x-x′|; θ)→0.常用的变异系数的形式有指数型:R(d; θ)=e-θdθ>0, 高斯型:R(d; θ)=e-θd2, θ>0等,本文在计算时使用高斯型的变异系数函数.f(x)是一组已知的基函数,β为系数,通常可以通过最小二乘法求得,在实际计算中,经常取f(x)Tβ=β0[4, 6].

假设已知x1, …, xm处的函数值Y(x1), …, Y(xm),记Y=(Y(x1), …, Y(xm))T,克里金插值通过Y(x1), …, Y(xm)的线性组合估计未知点x0处的函数值Y(x0),即:

$ \widehat Y({x_0}) = \mathit{\boldsymbol{\lambda }}{({x_0})^{\rm{T}}}\mathit{\boldsymbol{Y}} + {\lambda _0}({x_0}). $

其中: λ (x0)为m维向量; λ0(x0)为常数,通过最小化估计的均方误差选取最优的线性组合的系数λ(x0)和λ0(x0),最终$\widehat Y$(x0)有如下形式[12]

$ \widehat Y({x_0}) = {\beta _0} + {\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_M}{({x_0}, \cdot)^{\rm{T}}}{\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}_M}^{ - 1}(Y - {\beta _0}{\bf{1}}{_m}) $

其中: ΣM(x0, ·)为x1, …, xmx0的协方差向量; ΣMx1, …, xm之间的协方差矩阵; 1mm维全1向量.

本文将克里金插值应用在在库存管理问题中,首先选取m个自变量的取值x1, …, xm,模拟求得这m个点处成本函数和服务水平函数C(xi, wi, j)和F(xi, wi, j),其中wi, j表示第i个自变量取值的第j个场景模拟.用不同场景下的均值ϕn(x)和φn(x)作为对期望ϕ(x)和φ(x)的估计,再对这两个函数利用克里金插值,得到任一点x0处的取值$\widehat \phi $(x)和$\widehat \varphi $(x).使用拉丁超立方的方法选择初始m个点的取值[14].

3 蒙特卡罗树搜索

蒙特卡罗树搜索[13]分成4个迭代的步骤:

(1) 选择:从一个根节点出发,持续使用子节点的树搜索策略来扩展树结构.一个节点是可扩展的,如果该节点处于非终极状态并有未访问的子节点.

(2) 扩展:根据可行的动作将一个或者多个子节点加入树结构.

(3) 模拟:在新的节点根据默认策略来进行模拟,生成奖励(效用)函数的结果.

(4) 回溯:将模拟的结果回溯到选择的节点中,并更新对应节点相应的统计值.

现在将蒙特卡罗树搜索的方法用于随机优化问题.问题可以定义为, 对于一个N维超立方体Ω={xi∈[xli, xui],其中i=1, …, Nxi为第i维变量,xlixui分别为第i维变量的上下界,找出一个N维整数向量X,使效用函数f(X)的值最小(或最大).

相应的蒙特卡罗树搜索的算法如下:

(1) 选择:每次扩展后选择效用函数最优的子空间.

(2) 扩展:对选择的节点进行扩展,在选择的节点下生成2N个子节点.

(3) 模拟:在子空间内部进行均匀抽样并计算相应的效用函数.

(4) 回溯:将抽样计算的结果回溯到父节点,并计算相应的平均值和方差.

通过蒙特卡罗树搜索,在相应的计算资源的预算下,选择最优的子超立方体,然后将相应的子超立方体为根节点,重复进行蒙特卡罗树搜索,直至收敛到一个点为止.

4 蒙特卡罗方法和克里金方法的结合

针对多级库存系统问题,首先在原始空间中进行抽样, 再进行克里金回归,得出相应的成本函数和约束函数的插值函数形式$\widehat \phi $(x)和$\widehat \psi $(x).从而构造出效用函数如下:

$ {f_1} = \widehat \phi (\mathit{\boldsymbol{x}}) + p \cdot {I_{{\rm{sign}}}}({\rm{max}}({F_{{\rm{min}}}} - \widehat \psi (\mathit{\boldsymbol{x}}), 0)) $

其中:p为惩罚项系数; Fmin为约束条件的阀值; Isign(x)为示性函数.在进行蒙特卡罗树搜索的过程中可以直接使用构造出的效用函数进行计算,避免了反复的模拟抽样,大大提高了计算的效率.

当选择了一个新的子超立方体时,将在子超立方体里构造新的效用函数fi(x),i=2, …, M,保证计算的精度.最终的算法如图 3所示.图中, UCT表示搜索树的上置信确界.

图 3 蒙特卡罗树搜索与克里金相结合的算法 Fig.3 Monte Carlo tree simulation algorithm with Kriging
5 实际案例的数值结果

数值实验应用在图 1所示的库存系统结构中,随机参数的分布如图 1,且各个随机变量相互独立.Fmin=0.975,CO(1)=100,CO(i)=500,i=2, …, 4,CH(i)=1,i=1, …, 4,CS(1)=1 000,CS(i)=2 000, i=2, …, 4, T=120, 克里金差值样本点个数m=50,每个样本点模拟的场景数n=100.初始的搜索空间为x(i)∈[xl(i), xu(i)],当i=1, 2, 3, 5, 6, 7时,xl(i)=4 096,xu(i)=20 479,当i=4, 8时,xl(i)=32 768,xl(i)=65 536.搜索时间的范围从10~100 s.为了表明算法关于搜索时间的收敛效果,对每个搜索时间,运行十次改进的树搜索算法,比较不同搜索时间下所得最优解的均值和方差,具体结果见表 1.

下载CSV 表 1 数值计算结果 Tab.1 Numerical result

表 1中,ϕ(x)表示最优成本的均值, $\overline \psi $(x)表示及时供应率的均值;通过最优解的二范数的标准差s($\left\| X \right\|$)体现解的稳定性.通过表 1中的数据可知,利用蒙特卡罗树搜索和克里金差值结合的方法,所求得的解基本可以满足限制条件,同时随着搜索时间的增长,得到的最优解也越来越小,图 4图 5分别是通过搜索算法得到的最优解处的成本和及时供应率与搜索时间的关系.

图 4 最优值与搜索时间的关系 Fig.4 Relationship of ϕ(x) and search time
图 5 限制条件和搜索时间的关系 Fig.5 Relationship of $\overline \psi $(x) and search time

通过图 4图 5可以看出,随着搜索时间的增长,求得的最优值的大小不断下降,同时及时供应率也基本满足限制条件.图 6反映了所得最优解的稳定性和搜索时间的关系.

图 6 最优解标准差和搜索时间的关系 Fig.6 Relationship of s($\left\| x \right\|$) and search time

通过图 6可以看出,随着搜索时间的增长,最优解二范数的标准差也呈下降趋势,从一定程度上可以表明,所得最优解更加稳定.综合图 4~图 6,可以体现出算法关于搜索时间的收敛情况.针对图 5图 6的突变点,认为这是由于目标函数不够光滑和算法本身的随机性造成的.目标函数是使用随机样本均值代替的,在做克里金插值和蒙特卡罗树搜索时也有随机取样的步骤,这个也会造成曲线不够光滑的现象.为了进一步体现本文中提出的方法的优势,使用Chu等[3]所提出的算法对同一问题进行了计算,具体结果见表 2.

下载CSV 表 2 Chu等[3]算法的数值计算结果 Tab.2 Numerical result of the method in Chu, et al[3]

表 2列出了使用Chu等[3]所提出的算法对同一问题进行十次运算的结果.从结果中可以看出,这种方法所需时间为2 000 s左右,所得到的最优值和及时供应率在1010和0.95左右浮动.因此本文所提出的方法可以使用更少的时间得到更好的运算结果.

本文的方法应用到了某500强企业全球供应链的库存系统优化的工作中.该企业为制药企业,需要在保证全球病人医药供给的前提下,通过控制变量,降低药品的库存成本.由于实际问题中,业务流程比较复杂,黑箱函数需要对复杂的业务流程按照预先设定的时间粒度进行模拟,相应的计算量非常大.基于随机切平面方法的基准方法需要耗时平均约30 min,而使用本文改进的方法可以在100 s得到优化的结果,大大提高了计算的效率,从而帮助企业更好的进行实时决策.

6 结论

由于库存系统的复杂性和随机性,在解决库存系统高维优化问题时,需要求解的优化问题的目标函数和限制条件都是一个黑箱函数的期望形式,没有显示表达.因此对库存系统成本和服务水平的取值和优化需要耗费很大的运算量.如果使用基于传统的梯度下降方法,在耗费巨大运算量的同时,也很容易使求解陷入局部极小值点.

本文提出了一种结合蒙特卡罗树搜索法和克里金插值法的优化算法,可以高效地求解库存系统的优化问题,在满足系统服务水平大于某一阈值的情况下最小化总成本.蒙特卡罗树搜索是一种被广泛用于马尔可夫决策过程的方法,本文通过对搜索空间的不断切割,将优化问题转化为一个马尔可夫决策过程,并利用树搜索方法求解,使最终得到的结果收敛于全局最优解.

在求解时通过克里金插值可以减小运算量,提高搜索效率.在搜索过程中,在每一步搜索的子空间上利用克里金插值的显示公式得到对应的成本和服务水平的取值,使问题可以更高效地求解.

通过数值计算可以看出,改进的树搜索方法关于搜索时间收敛,随着搜索时间的增长,最优值越来越小,且求得的解越来越稳定.

参考文献
[1]
WASSICK J M, AGARWAL A, AKIYA N, et al. Addressing the operational challenges in the development, manufacture, and supply of advanced materials and performance products[J]. Computers & Chemical Engineering, 2012, 47: 157
[2]
黄晟.基于集中控制策略的汽车维修备件分布式多级库存研究[D].上海: 东华大学, 2013.
HUANG Sheng. A study on distributed multi-Level inventory system of automobile maintenance parts in centralized control background. Shanghai: Donghua University, 2013. http://cdmd.cnki.com.cn/Article/CDMD-10255-1013161940.htm
[3]
CHU Y, YOU F, WASSICK J M, et al. Simulation-based optimization framework for multi-echelon inventory systems under uncertainty[J]. Computers & Chemical Engineering, 2015, 73: 1
[4]
BILES W E, KLEIJNEN J P C, BEERS W C M V, et al. Kriging metamodeling in constrained simulation optimization: an explorative study[C]// Conference on Winter Simulation: 40 Years! the Best Is Yet to Come.[S.l.]: IEEE Press, 2007: 355-362.
[5]
李晓军, 李培楠, 朱合华, 等. 基于贝叶斯克里金的地下空间多源数据建模[J]. 同济大学学报(自然科学版), 2014, 42(3): 406
LI Xiaojun, LI Peinan, ZHU Hehua, et al. Geomodeling with integration of multi source data by Bayesian Kriging in underground space[J]. Journal of Tongji University (Natural Science), 2014, 42(3): 406 DOI:10.3969/j.issn.0253-374x.2014.03.013
[6]
ANKENMAN B, NELSON B L, STAUM J. Stochastic Kriging for simulation metamodeling[J]. Operations Research, 2010, 58(2): 371
[7]
JONES D R, SCHONLAU M, WELCH W J. Efficient global optimization of expensive black-box functions[J]. Journal of Global Optimization, 1998, 13(4): 455 DOI:10.1023/A:1008306431147
[8]
HUANG D, ALLEN T T, NOTZ W I, et al. Global optimization of stochastic black-box systems via sequential Kriging meta-models[J]. Journal of Global Optimization, 2006, 34(3): 441 DOI:10.1007/s10898-005-2454-3
[9]
SILVER D, HUANG A, MADDISON C J, et al. Mastering the game of Go with deep neural networks and tree search[J]. Nature, 2016, 529(7587): 484 DOI:10.1038/nature16961
[10]
KIM S, PASUPATHY R, HENDERSON S G. A guide to sample average approximation—handbook of simulation optimization[M]. New York: Springer, 2015
[11]
RAHMANDAD H, STERMAN J. Heterogeneity and network structure in the dynamics of diffusion: comparing agent-based and differential equation models[J]. Management Science, 2008, 54(5): 998 DOI:10.1287/mnsc.1070.0787
[12]
邹林君.基于Kriging模型的全局优化方法研究[D].武汉: 华中科技大学, 2011.
ZOU Linjun. A search of global optimization method based on kriging model[D]. Wuhan: Huazhong University of Science and Technology, 2011.
[13]
BROWNE C B, POWLEY E, WHITEHOUSE D, et al. A survey of Monte Carlo tree search methods[J]. IEEE Transactions on Computational Intelligence & Ai in Games, 2012, 4(1): 1