曹小華 劉道凡 涂圣才

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

0 引 言

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

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

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

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

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

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

圖1 典型的物流運(yùn)輸節(jié)點(diǎn)

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

(1)

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

(2)

阻塞預(yù)警模型為

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

(3)

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

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

L=t-tcurr

(4)

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

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

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

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

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

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

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

Δti=ti+1-ti

(5)

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

定義t=ti時(shí)刻該輸入流的流速度為

(6)

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

2.1.2等間隔化

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

(7)

(8)

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

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

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

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

(9)

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

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

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

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

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

(10)

2.2.1趨勢(shì)項(xiàng)模型

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

g(t)=k0t+b0

(11)

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

(12)

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

(13)

將其寫成更簡(jiǎn)潔的向量形式

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

(14)

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

(15)

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

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

(16)

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

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

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

(17)

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

2.2.2季節(jié)周期性

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

(18)

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

(19)

則s(t)為

s(t)=x(t)β

(20)

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

3 實(shí)驗(yàn)與結(jié)果

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

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

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

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

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

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

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

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

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

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

(21)

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

(22)

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

(23)

3.3 實(shí)驗(yàn)結(jié)果對(duì)比分析

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

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

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

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

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

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

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

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

4 結(jié) 束 語(yǔ)

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