曹小华 刘道凡 涂圣才
(武汉理工大学物流工程学院 武汉 430063)
生产物流自动化系统可以实现不同场地、工序或设备之间物料等的传送.既定生产物流系统中各环节本身及彼此间的依赖关系在短期内难以调整.通常各环节都遵循生产计划,工作进程合理有序,但仍会出现异常情况,如数据错误、设备故障、人为失误或业务变动.张俊[1]提出了装配生产过程的实时状态模型,应用RFID(射频识别)技术进行实时数据采集,实现了电子看板的实时生产数据可视化.赵浩然等[2]分析了数字孪生车间与三维可视化监控之间的关系,提出一种三维可视化监控模式和实时数据驱动的虚拟车间运行模式.Cao等[3]建立了RFID数据的计算模型以获取多属性生产物流的状态数据,提出了一种相似度模型表征状态特征,结合聚类方法挖掘出生产物流中的异常.
关于生产物流的异常问题目前多采用实时监控的方案,即基于既定事实进行状态判别.对于生产物流运输路径的阻塞异常,尚未形成提前预警的成熟方案,实现阻塞提前预警可避免因既定异常而导致经济损失.在复杂设备的故障诊断领域,利用有效的时间序列预测模型进行故障预警是一种可行的方案[4].
关于时序预测的研究主要有:①时域分析法,如ARIMA(求和自回归滑动平均)和ARCH(自回归条件异方差)[5-6];②频域分析法,常见如傅里叶变换;③基于机器学习的方法,如BP神经网络[7]和LSTM(长短时记忆神经网络)[8]等.不同算法性能各异,同时大多数算法的预测精度随步长增加而显著下降.
为了实现对物流运输节点阻塞异常的预警,建立了描述运输节点阻塞的预警模型,并提出一种基于Prophet框架的流速度预测方法,基于实际生产物流系统中采集的不同情况下的数据,实现了对运输路径的流入/流出速度的有效预测,验证了该预测方法用于阻塞异常预警的有效性.
一个典型的物流运输节点见图1,其具有n个输入流和m个输出流,缓冲容量为b.在运输节点处,定义输入流的流入速度vin;在单位时间内流入该节点的产品数量,定义输出流的流出速度vout;在单位时间内流出该节点的产品数量.
图1 典型的物流运输节点
多流入多流出节点和输入/输出流的典型实例是生产物流系统中的缓存区和自动化输送线.假设当前时刻为tcurr,缓冲容量为bcurr,则在任意时刻t(t≥tcurr),运输节点的净流入量为
(1)
将式按等时间间隔Δt离散化,得到离散化式.
(2)
阻塞预警模型为
Δ(t)=bcurr-Q(t)-bthre
(3)
式中:bthre为设定的缓冲容量阈值.当Δ(t)<0时发生阻塞.
在当前时刻tcurr若能预测tcurr~t时间段内的流入与流出速度,即可预算出累积运输量,再根据式对未来t时刻的节点阻塞状态进行提前预警.定义阻塞预警提前期为
L=t-tcurr
(4)
在阻塞预警模型的基础上,将提出基于Prophet框架的流入/流出速度的预测方法,该预测方法的输入是tcurr时刻前流入/流出速度的历史序列,输出是阻塞预警提前期L时段内流入/流出速度的预测序列.此后用流速度v来统一表示vin和vout.
在实际生产过程中,RFID等设备会记录下每个产品流入/流出节点的时间点,见表1.表1前两列是2018年11月14日某物流公司的某运输节点处每一批(20个)产品流入该节点的时间点.但是从表1前两列无法直接得到适用于Prophet的流速度序列.
表1 每批产品流入节点的时间点及流速度序列
2.1.1数据转换
为了得到流速度序列,把相邻两批产品流入节点的时间差定义为时间间隔Δti,即
Δti=ti+1-ti
(5)
式中:ti为第i批产品流入节点的时间点,s/(20个).
定义t=ti时刻该输入流的流速度为
(6)
式中:v(ti)的单位为(20个)/min.得到流速度v的非等间隔的时间序列{v(ti),i=0,1,…,n},见表1第三列.
2.1.2等间隔化
适用于非等间隔序列的直接建模方法较少且不易实现[9],采用分段三次Hermite插值法得到等间隔数据.将上节得到的非等间隔的流速度序列{v(ti)}去除掉绝对时间信息,得到序列{u(t)},即u(t)=v(t+t1),见图2中样本节点.对于数据点集{(xk,yk),k=1,2,…,n},在相邻两个节点所构成的每个子段[xk,xk+1]上定义插值函数H3(x),即插值点x所对应的值为
(7)
(8)
端点x1和xn处的一阶导数值则取为相邻区段一阶差商的3倍.对于数据点集{(t,u(t)),t=0,t2-t1,…,tn-t1},从t=0开始以步长h=60 s作插值点进行插值,得到结果见图2.将等间隔数据点集恢复绝对时间信息得到等间隔化流速度序列{v(t),t=t1,60+t1,…,60l+t1},即v(t)=H3(t-t1),见图3.
由图3可知,流速度序列呈上升趋势并伴有波动,这可能是多种因素导致的,提出一种基于Prophet框架的流速度预测模型,该方法只考虑流速度自身随时间的变化.Prophet算法的基本思想是拟合历史数据的相关特性以学习出数据的走向.根据Cramer分解定理,时间序列可表示为[10]
y(t)=g(t)+s(t)+h(t)+εt
(9)
式中:y(t)为时间序列在时刻t的观测值;趋势项g(t)为序列的总体变化趋势;周期项s(t)为季节项,体现出序列的周期循环性质;节假日项h(t)为有规律的突发事件效应;高斯白噪声εt具有随机平稳特征,已无可提取的信息.
图2 分段三次Hermite等间隔化插值结果
图3 等间隔化的流速度序列
结合实际情况,剔除节假日项h(t),将y(t)表示为
y(t)=g(t)+s(t)+εt
(10)
2.2.1趋势项模型
Prophet算法中趋势项的核心函数是逻辑回归函数和分段线性函数.这里采用基于分段线性函数的趋势项模型.对于线性函数式,将斜率(增长率)项k0与偏置项b0稍加改造得到分段线性函数式.
g(t)=k0t+b0
(11)
g(t)=k(t)·t+b(t)
(12)
斜率项k(t)处于变化状态且未知待求,因此需从历史数据中检测出斜率突变的点.考虑在序列样本的时间跨度T内设置m个变点sj(j=1,2,…,m),构成变点向量S.默认在80%样本数据量内均匀设置25个变点.在变点sj处斜率的增量定义为δj,构成斜率增量列向量δ∈Rm.于是,斜率函数
(13)
将其写成更简洁的向量形式
k(t)=k0+a(t)Tδ
(14)
式中:列向量a(t)∈{0,1}m的元素aj(t)满足
(15)
考虑到变点造成的函数值的非连续性,需对偏置项b(t)做出调整,调整后的偏置项为
b(t)=b0+a(t)Tγ
(16)
式中:列向量γ={γj(j=1,2,…,m)}m,并且γj=-sjδj.
于是基于分段线性函数的趋势项模型改造为
g(t)=[k0+a(t)Tδ]t+[b0+a(t)Tγ]
(17)
参数δj是未知的,为防止g(t)发生过拟合现象,需加入稀疏先验对模型进行L1正则化,因此δ~Laplace(0,τ).τ的默认值为0.05,其值较大时模型可适应较大的增长趋势,但同时也将增加预测值的不确定性.
2.2.2季节周期性
考虑周期项s(t),将周期项s(t)用有限阶的离散傅里叶级数近似表示,即
(18)
式中:P为序列的周期;N为展开的阶数;β=[a1,b1,…,aN,bN]T为待估计的参数向量.如果令
(19)
则s(t)为
s(t)=x(t)β
(20)
同参数向量δ类似,需对模型进行L2正则化,因此β~Normal(0,σ2).σ值的大小体现出其对周期性的调节强度.序列展开阶数N相当于对周期项进行低通滤波,会产生截断误差,N越大,越能适应更复杂的周期性,但有可能使模型过拟合,可根据AIC等准则对模型进行参数的自动选择.
预处理后的数据划分成训练数据集(时间跨度为T)以及测试数据集(时间跨度为L).将训练数据输入预测模型得到时间跨度为L的预测值,并将Prophet模型的预测结果与LSTM和SARIMA(季节求和自回归滑动平均)的预测结果进行统计对比.
预处理后的序列的时间粒度为1 min,共201个数据,前181个为训练数据,后20个为测试数据.预测模型的参数设置见表2,采用最大后验估计(MAP)实现参数估计与区间估计.
表2 预测模型的参数设置
将训练数据输入Prophet模型,提取出流速度序列的组成成分见图4.图4a)为序列的总体趋势走向,图4b)为序列的周期特征,将两个成分应用加法得到图4c)的组合模型,如中间曲线所示,并给出了置信度为0.95下的序列估计区间.图4c)中折线显示了变点位置及增长率变化情况.
图4 流速度序列训练数据集的分解效果
根据从训练数据中学习到的模型待估参数和序列特征,得到图5的预测结果.图5为预测序列的值及置信度为0.95下的序列区间估计值,以及测试序列的实际值.
图5 模型的预测结果
均方根误差(eRMSE).预测值与实际值差值的向量2范数.
(21)
平均绝对误差(eMAE).预测值与实际值差值的向量1范数.
(22)
平均绝对百分比误差(eMAPE),为绝对误差与实际值之比的平均值,即
(23)
SARIMA的基本思想是对含有趋势性和周期性的序列进行d阶差分与S步差分,若能得到平稳序列则建立自回归滑动平均模型,并经白噪声检验、模型定阶与参数估计、模型检验后进行预测,再还原到原始序列的预测值.LSTM实现多步预测的基本思想是先从序列首值开始提取一定长度的序列,并按一定比例切分成输入数据(即训练步长)与输出数据(即标签),随后逐步滚动,提取出多组输入数据与标签训练LSTM网络,通过LSTM网络中的遗忘门、输入门、输出门及记忆单元来调节网络记忆状态,最后经迭代预测得到与测试序列长度相当的预测序列[11].
对于训练数据,经由一系列建模过程最终确定SARIMA模型为SARIMA(1,0,1)×(1,0,1)10后进行动态预测.LSTM模型的训练步长与预测步长分别设为40和20,每个LSTM单元的隐含层神经元数设为100,激活函数采用线性整流函数(ReLU),优化器采用Adam算法,损失函数采用均方误差(MSE),训练次数epochs=100.三种算法的预测结果对比见图6.
图6 三种算法的预测结果对比
3种算法的性能指标见表3.由图6和表3可知,LSTM模型相较于其余两者并未获得较好的预测效果;SARIMA与Prophet模型均表现出不错的效果,各统计指标均不超过1%.从模型结构及建模过程的复杂度来看,Prophet模型最简洁,LSTM模型最为复杂,SARIMA的建模及模型解释需相当深厚的专业知识,后两者都不太适用于物流流速度的预测.
表3 3种算法的性能指标
SARIMA等传统模型对平稳序列及较弱的非线性序列的分析效果显著,但是对于存在阶跃突变的序列,建模过程变得复杂甚至不可控,预测效果往往不理想.图7为存在阶跃突变的某等间隔流速度序列,Prophet的预测效果优于SARIMA,SARIMA中较优的SARIMA(1,0,0)×(0,0,0)20的预测结果近似于直线,意义不大,说明SARIMA模型不太适用于分段差分平稳序列的预测.
图7 序列存在突变的预测结果
在测试集上,Prophet方法的eRMSE、eMAE和eMAPE分别为0.063 8,0.048 1和0.042 5,而SARIMA方法的各项指标分别为0.097 4,0.079 9和0.666,均高于Prophet方法.物流生产过程中由于任务订单变化、生产线突发汇聚或分流,均会出现流速度序列阶跃等现象,Prophet模型具有突变点检测的功能,因此对此类干扰有较强的适应能力,可以获得较为满意的预测效果.
针对生产物流阻塞异常的问题,建立了表征阻塞状态的阻塞预警模型,通过预测出各运输路径的运输速度及累积运输量以实现对该物流运输节点的阻塞预警.在对原始物流数据进行数据转换及等间隔化的基础上,提出了基于Prophet框架的物流运输节点的流速度预测方法,实验结果表明该预测方法性能优良,实现了对运输路径的流速度的有效预测,并可以适应流速度突变所带来的影响,性能优于经典的SARIMA和LSTM预测方法.由于所提出的预测模型仅考虑流速度变量随时间的变化,但是实际生产过程中情况复杂多变,计划将其他因素引入预测模型,实现更准确地预警.