基于改进树搜索方法的库存系统优化

牟 刚, 孙永超, 徐承龙

(1.同济大学 数学科学学院,上海 200092;2.上海财经大学 数学学院,上海 200433)

库存系统优化是供应链管理中的一个重大问题[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 多级库存系统的结构

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

图2 (r,Q)策略下库存水平变化情况

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

在求解最优库存策略时,目标函数为某段时间内的总成本,限制条件为这段时间内的服务水平大于某一常数.总成本分为三部分,每天每单位在运输中将要运输至第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.需要求解的优化问题为

s.t.ψ(x)=Ew[F(x,w)]…Fmin

这个随机优化问题的目标函数和限制条件都是某个随机变量的期望,且这个期望没有显式表达式,因此需要利用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)仍然没有显示表达式,这是因为无法直接通过x和wi直接得到观测时间内节点间的订单和运输情况,必须对于整个库存系统的运行过程进行模拟,才能根据库存策略判断订单和运输的产生,进而得到成本和服务水平.基于代理的模拟可以方便地模拟库存系统的动态过程,Chu等[3],Rahmandad等[11]都使用这一方法对库存系统进行模拟.

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

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

2 克里金插值

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

Y(x)=f(x)Tβ+M(x)

其中,M是从某一高斯随机场的一次抽样,高斯随机场满足条件E[M(X)]=0,E[M2(x)]<∞,且当x和x′的取值相近时,M(x)和M(x′)的取值也趋于相近,M(x)和M(x′)的协方差可以写成τ2R(|x-x′|;θ),其中τ2为随机场的方差,R(|x-x′|;θ)为x和x′之间的变异系数,θ是这个函数的参数,变异系数函数需要满足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),即:

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

3 蒙特卡罗树搜索

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

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

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

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

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

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

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

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

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

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

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

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

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

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

5 实际案例的数值结果

图3 蒙特卡罗树搜索与克里金相结合的算法

Fig.3 Monte Carlo tree simulation algorithm with Kriging

表1 数值计算结果

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

图4 最优值与搜索时间的关系

图5 限制条件和搜索时间的关系

图6 最优解标准差和搜索时间的关系

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

表2 Chu等[3]算法的数值计算结果

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

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

6 结论

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

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

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

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