曹小華 劉道凡 涂圣才

(武漢理工大學(xué)物流工程學(xué)院 武漢 430063)

0 引 言

生產(chǎn)物流自動化系統(tǒng)可以實現(xiàn)不同場地、工序或設(shè)備之間物料等的傳送.既定生產(chǎn)物流系統(tǒng)中各環(huán)節(jié)本身及彼此間的依賴關(guān)系在短期內(nèi)難以調(diào)整.通常各環(huán)節(jié)都遵循生產(chǎn)計劃，工作進(jìn)程合理有序，但仍會出現(xiàn)異常情況，如數(shù)據(jù)錯誤、設(shè)備故障、人為失誤或業(yè)務(wù)變動.張俊[1]提出了裝配生產(chǎn)過程的實時狀態(tài)模型，應(yīng)用RFID(射頻識別)技術(shù)進(jìn)行實時數(shù)據(jù)采集，實現(xiàn)了電子看板的實時生產(chǎn)數(shù)據(jù)可視化.趙浩然等[2]分析了數(shù)字孿生車間與三維可視化監(jiān)控之間的關(guān)系，提出一種三維可視化監(jiān)控模式和實時數(shù)據(jù)驅(qū)動的虛擬車間運行模式.Cao等[3]建立了RFID數(shù)據(jù)的計算模型以獲取多屬性生產(chǎn)物流的狀態(tài)數(shù)據(jù)，提出了一種相似度模型表征狀態(tài)特征，結(jié)合聚類方法挖掘出生產(chǎn)物流中的異常.

關(guān)于生產(chǎn)物流的異常問題目前多采用實時監(jiān)控的方案，即基于既定事實進(jìn)行狀態(tài)判別.對于生產(chǎn)物流運輸路徑的阻塞異常，尚未形成提前預(yù)警的成熟方案，實現(xiàn)阻塞提前預(yù)警可避免因既定異常而導(dǎo)致經(jīng)濟(jì)損失.在復(fù)雜設(shè)備的故障診斷領(lǐng)域，利用有效的時間序列預(yù)測模型進(jìn)行故障預(yù)警是一種可行的方案[4].

關(guān)于時序預(yù)測的研究主要有：①時域分析法，如ARIMA(求和自回歸滑動平均)和ARCH(自回歸條件異方差)[5-6]；②頻域分析法，常見如傅里葉變換；③基于機(jī)器學(xué)習(xí)的方法，如BP神經(jīng)網(wǎng)絡(luò)[7]和LSTM(長短時記憶神經(jīng)網(wǎng)絡(luò))[8]等.不同算法性能各異，同時大多數(shù)算法的預(yù)測精度隨步長增加而顯著下降.

為了實現(xiàn)對物流運輸節(jié)點阻塞異常的預(yù)警，建立了描述運輸節(jié)點阻塞的預(yù)警模型，并提出一種基于Prophet框架的流速度預(yù)測方法，基于實際生產(chǎn)物流系統(tǒng)中采集的不同情況下的數(shù)據(jù)，實現(xiàn)了對運輸路徑的流入/流出速度的有效預(yù)測，驗證了該預(yù)測方法用于阻塞異常預(yù)警的有效性.

1 物流運輸節(jié)點的阻塞預(yù)警模型

一個典型的物流運輸節(jié)點見圖1，其具有n個輸入流和m個輸出流，緩沖容量為b.在運輸節(jié)點處，定義輸入流的流入速度vin；在單位時間內(nèi)流入該節(jié)點的產(chǎn)品數(shù)量，定義輸出流的流出速度vout；在單位時間內(nèi)流出該節(jié)點的產(chǎn)品數(shù)量.

圖1 典型的物流運輸節(jié)點

多流入多流出節(jié)點和輸入/輸出流的典型實例是生產(chǎn)物流系統(tǒng)中的緩存區(qū)和自動化輸送線.假設(shè)當(dāng)前時刻為tcurr，緩沖容量為bcurr，則在任意時刻t(t≥tcurr)，運輸節(jié)點的凈流入量為

(1)

將式按等時間間隔Δt離散化，得到離散化式.

(2)

阻塞預(yù)警模型為

Δ(t)=bcurr-Q(t)-bthre

(3)

式中：bthre為設(shè)定的緩沖容量閾值.當(dāng)Δ(t)<0時發(fā)生阻塞.

在當(dāng)前時刻tcurr若能預(yù)測tcurr～t時間段內(nèi)的流入與流出速度，即可預(yù)算出累積運輸量，再根據(jù)式對未來t時刻的節(jié)點阻塞狀態(tài)進(jìn)行提前預(yù)警.定義阻塞預(yù)警提前期為

L=t-tcurr

(4)

2 流速度預(yù)測方法

在阻塞預(yù)警模型的基礎(chǔ)上，將提出基于Prophet框架的流入/流出速度的預(yù)測方法，該預(yù)測方法的輸入是tcurr時刻前流入/流出速度的歷史序列，輸出是阻塞預(yù)警提前期L時段內(nèi)流入/流出速度的預(yù)測序列.此后用流速度v來統(tǒng)一表示vin和vout.

2.1 數(shù)據(jù)預(yù)處理

在實際生產(chǎn)過程中，RFID等設(shè)備會記錄下每個產(chǎn)品流入/流出節(jié)點的時間點，見表1.表1前兩列是2018年11月14日某物流公司的某運輸節(jié)點處每一批(20個)產(chǎn)品流入該節(jié)點的時間點.但是從表1前兩列無法直接得到適用于Prophet的流速度序列.

表1 每批產(chǎn)品流入節(jié)點的時間點及流速度序列

2.1.1數(shù)據(jù)轉(zhuǎn)換

為了得到流速度序列，把相鄰兩批產(chǎn)品流入節(jié)點的時間差定義為時間間隔Δti，即

Δti=ti+1-ti

(5)

式中：ti為第i批產(chǎn)品流入節(jié)點的時間點，s/(20個).

定義t=ti時刻該輸入流的流速度為

(6)

式中：v(ti)的單位為(20個)/min.得到流速度v的非等間隔的時間序列{v(ti),i=0,1,…,n}，見表1第三列.

2.1.2等間隔化

適用于非等間隔序列的直接建模方法較少且不易實現(xiàn)[9]，采用分段三次Hermite插值法得到等間隔數(shù)據(jù).將上節(jié)得到的非等間隔的流速度序列{v(ti)}去除掉絕對時間信息，得到序列{u(t)}，即u(t)=v(t+t1)，見圖2中樣本節(jié)點.對于數(shù)據(jù)點集{(xk,yk),k=1,2,…,n}，在相鄰兩個節(jié)點所構(gòu)成的每個子段[xk,xk+1]上定義插值函數(shù)H3(x)，即插值點x所對應(yīng)的值為

(7)

(8)

端點x1和xn處的一階導(dǎo)數(shù)值則取為相鄰區(qū)段一階差商的3倍.對于數(shù)據(jù)點集{(t,u(t)),t=0,t2-t1,…,tn-t1}，從t=0開始以步長h=60 s作插值點進(jìn)行插值，得到結(jié)果見圖2.將等間隔數(shù)據(jù)點集恢復(fù)絕對時間信息得到等間隔化流速度序列{v(t),t=t1,60+t1,…,60l+t1}，即v(t)=H3(t-t1)，見圖3.

2.2 基于Prophet框架的流速度預(yù)測模型

由圖3可知，流速度序列呈上升趨勢并伴有波動，這可能是多種因素導(dǎo)致的，提出一種基于Prophet框架的流速度預(yù)測模型，該方法只考慮流速度自身隨時間的變化.Prophet算法的基本思想是擬合歷史數(shù)據(jù)的相關(guān)特性以學(xué)習(xí)出數(shù)據(jù)的走向.根據(jù)Cramer分解定理，時間序列可表示為[10]

y(t)=g(t)+s(t)+h(t)+εt

(9)

式中：y(t)為時間序列在時刻t的觀測值；趨勢項g(t)為序列的總體變化趨勢；周期項s(t)為季節(jié)項，體現(xiàn)出序列的周期循環(huán)性質(zhì)；節(jié)假日項h(t)為有規(guī)律的突發(fā)事件效應(yīng)；高斯白噪聲εt具有隨機(jī)平穩(wěn)特征，已無可提取的信息.

圖2 分段三次Hermite等間隔化插值結(jié)果

圖3 等間隔化的流速度序列

結(jié)合實際情況，剔除節(jié)假日項h(t)，將y(t)表示為

y(t)=g(t)+s(t)+εt

(10)

2.2.1趨勢項模型

Prophet算法中趨勢項的核心函數(shù)是邏輯回歸函數(shù)和分段線性函數(shù).這里采用基于分段線性函數(shù)的趨勢項模型.對于線性函數(shù)式，將斜率(增長率)項k0與偏置項b0稍加改造得到分段線性函數(shù)式.

g(t)=k0t+b0

(11)

g(t)=k(t)·t+b(t)

(12)

斜率項k(t)處于變化狀態(tài)且未知待求，因此需從歷史數(shù)據(jù)中檢測出斜率突變的點.考慮在序列樣本的時間跨度T內(nèi)設(shè)置m個變點sj(j=1,2,…,m)，構(gòu)成變點向量S.默認(rèn)在80%樣本數(shù)據(jù)量內(nèi)均勻設(shè)置25個變點.在變點sj處斜率的增量定義為δj，構(gòu)成斜率增量列向量δ∈Rm.于是，斜率函數(shù)

(13)

將其寫成更簡潔的向量形式

k(t)=k0+a(t)Tδ

(14)

式中：列向量a(t)∈{0,1}m的元素aj(t)滿足

(15)

考慮到變點造成的函數(shù)值的非連續(xù)性，需對偏置項b(t)做出調(diào)整，調(diào)整后的偏置項為

b(t)=b0+a(t)Tγ

(16)

式中：列向量γ={γj(j=1,2,…,m)}m，并且γj=-sjδj.

于是基于分段線性函數(shù)的趨勢項模型改造為

g(t)=[k0+a(t)Tδ]t+[b0+a(t)Tγ]

(17)

參數(shù)δj是未知的，為防止g(t)發(fā)生過擬合現(xiàn)象，需加入稀疏先驗對模型進(jìn)行L1正則化，因此δ～Laplace(0,τ).τ的默認(rèn)值為0.05，其值較大時模型可適應(yīng)較大的增長趨勢，但同時也將增加預(yù)測值的不確定性.

2.2.2季節(jié)周期性

考慮周期項s(t)，將周期項s(t)用有限階的離散傅里葉級數(shù)近似表示，即

(18)

式中：P為序列的周期；N為展開的階數(shù)；β=[a1,b1,…,aN,bN]T為待估計的參數(shù)向量.如果令

(19)

則s(t)為

s(t)=x(t)β

(20)

同參數(shù)向量δ類似，需對模型進(jìn)行L2正則化，因此β～Normal(0,σ2).σ值的大小體現(xiàn)出其對周期性的調(diào)節(jié)強(qiáng)度.序列展開階數(shù)N相當(dāng)于對周期項進(jìn)行低通濾波，會產(chǎn)生截斷誤差，N越大，越能適應(yīng)更復(fù)雜的周期性，但有可能使模型過擬合，可根據(jù)AIC等準(zhǔn)則對模型進(jìn)行參數(shù)的自動選擇.

3 實驗與結(jié)果

預(yù)處理后的數(shù)據(jù)劃分成訓(xùn)練數(shù)據(jù)集(時間跨度為T)以及測試數(shù)據(jù)集(時間跨度為L).將訓(xùn)練數(shù)據(jù)輸入預(yù)測模型得到時間跨度為L的預(yù)測值，并將Prophet模型的預(yù)測結(jié)果與LSTM和SARIMA(季節(jié)求和自回歸滑動平均)的預(yù)測結(jié)果進(jìn)行統(tǒng)計對比.

3.1 模型設(shè)置與預(yù)測結(jié)果

預(yù)處理后的序列的時間粒度為1 min，共201個數(shù)據(jù)，前181個為訓(xùn)練數(shù)據(jù)，后20個為測試數(shù)據(jù).預(yù)測模型的參數(shù)設(shè)置見表2，采用最大后驗估計(MAP)實現(xiàn)參數(shù)估計與區(qū)間估計.

表2 預(yù)測模型的參數(shù)設(shè)置

將訓(xùn)練數(shù)據(jù)輸入Prophet模型，提取出流速度序列的組成成分見圖4.圖4a)為序列的總體趨勢走向，圖4b)為序列的周期特征，將兩個成分應(yīng)用加法得到圖4c)的組合模型，如中間曲線所示，并給出了置信度為0.95下的序列估計區(qū)間.圖4c)中折線顯示了變點位置及增長率變化情況.

圖4 流速度序列訓(xùn)練數(shù)據(jù)集的分解效果

根據(jù)從訓(xùn)練數(shù)據(jù)中學(xué)習(xí)到的模型待估參數(shù)和序列特征，得到圖5的預(yù)測結(jié)果.圖5為預(yù)測序列的值及置信度為0.95下的序列區(qū)間估計值，以及測試序列的實際值.

圖5 模型的預(yù)測結(jié)果

3.2 性能度量指標(biāo)

均方根誤差(eRMSE).預(yù)測值與實際值差值的向量2范數(shù).

(21)

平均絕對誤差(eMAE).預(yù)測值與實際值差值的向量1范數(shù).

(22)

平均絕對百分比誤差(eMAPE)，為絕對誤差與實際值之比的平均值，即

(23)

3.3 實驗結(jié)果對比分析

SARIMA的基本思想是對含有趨勢性和周期性的序列進(jìn)行d階差分與S步差分，若能得到平穩(wěn)序列則建立自回歸滑動平均模型，并經(jīng)白噪聲檢驗、模型定階與參數(shù)估計、模型檢驗后進(jìn)行預(yù)測，再還原到原始序列的預(yù)測值.LSTM實現(xiàn)多步預(yù)測的基本思想是先從序列首值開始提取一定長度的序列，并按一定比例切分成輸入數(shù)據(jù)(即訓(xùn)練步長)與輸出數(shù)據(jù)(即標(biāo)簽)，隨后逐步滾動，提取出多組輸入數(shù)據(jù)與標(biāo)簽訓(xùn)練LSTM網(wǎng)絡(luò)，通過LSTM網(wǎng)絡(luò)中的遺忘門、輸入門、輸出門及記憶單元來調(diào)節(jié)網(wǎng)絡(luò)記憶狀態(tài)，最后經(jīng)迭代預(yù)測得到與測試序列長度相當(dāng)?shù)念A(yù)測序列[11].

對于訓(xùn)練數(shù)據(jù)，經(jīng)由一系列建模過程最終確定SARIMA模型為SARIMA(1,0,1)×(1,0,1)10后進(jìn)行動態(tài)預(yù)測.LSTM模型的訓(xùn)練步長與預(yù)測步長分別設(shè)為40和20，每個LSTM單元的隱含層神經(jīng)元數(shù)設(shè)為100，激活函數(shù)采用線性整流函數(shù)(ReLU)，優(yōu)化器采用Adam算法，損失函數(shù)采用均方誤差(MSE)，訓(xùn)練次數(shù)epochs=100.三種算法的預(yù)測結(jié)果對比見圖6.

圖6 三種算法的預(yù)測結(jié)果對比

3種算法的性能指標(biāo)見表3.由圖6和表3可知，LSTM模型相較于其余兩者并未獲得較好的預(yù)測效果；SARIMA與Prophet模型均表現(xiàn)出不錯的效果，各統(tǒng)計指標(biāo)均不超過1%.從模型結(jié)構(gòu)及建模過程的復(fù)雜度來看，Prophet模型最簡潔，LSTM模型最為復(fù)雜，SARIMA的建模及模型解釋需相當(dāng)深厚的專業(yè)知識，后兩者都不太適用于物流流速度的預(yù)測.

表3 3種算法的性能指標(biāo)

SARIMA等傳統(tǒng)模型對平穩(wěn)序列及較弱的非線性序列的分析效果顯著，但是對于存在階躍突變的序列，建模過程變得復(fù)雜甚至不可控，預(yù)測效果往往不理想.圖7為存在階躍突變的某等間隔流速度序列，Prophet的預(yù)測效果優(yōu)于SARIMA，SARIMA中較優(yōu)的SARIMA(1,0,0)×(0,0,0)20的預(yù)測結(jié)果近似于直線，意義不大，說明SARIMA模型不太適用于分段差分平穩(wěn)序列的預(yù)測.

圖7 序列存在突變的預(yù)測結(jié)果

在測試集上，Prophet方法的eRMSE、eMAE和eMAPE分別為0.063 8，0.048 1和0.042 5，而SARIMA方法的各項指標(biāo)分別為0.097 4，0.079 9和0.666，均高于Prophet方法.物流生產(chǎn)過程中由于任務(wù)訂單變化、生產(chǎn)線突發(fā)匯聚或分流，均會出現(xiàn)流速度序列階躍等現(xiàn)象，Prophet模型具有突變點檢測的功能，因此對此類干擾有較強(qiáng)的適應(yīng)能力，可以獲得較為滿意的預(yù)測效果.

4 結(jié) 束 語

針對生產(chǎn)物流阻塞異常的問題，建立了表征阻塞狀態(tài)的阻塞預(yù)警模型，通過預(yù)測出各運輸路徑的運輸速度及累積運輸量以實現(xiàn)對該物流運輸節(jié)點的阻塞預(yù)警.在對原始物流數(shù)據(jù)進(jìn)行數(shù)據(jù)轉(zhuǎn)換及等間隔化的基礎(chǔ)上，提出了基于Prophet框架的物流運輸節(jié)點的流速度預(yù)測方法，實驗結(jié)果表明該預(yù)測方法性能優(yōu)良，實現(xiàn)了對運輸路徑的流速度的有效預(yù)測，并可以適應(yīng)流速度突變所帶來的影響，性能優(yōu)于經(jīng)典的SARIMA和LSTM預(yù)測方法.由于所提出的預(yù)測模型僅考慮流速度變量隨時間的變化，但是實際生產(chǎn)過程中情況復(fù)雜多變，計劃將其他因素引入預(yù)測模型，實現(xiàn)更準(zhǔn)確地預(yù)警.