Robust Interval Prediction of Intermittent Demand for Spare Parts Based on Tensor Optimization

Abstract

Demand for spare parts, which is triggered by element failure, project schedule and reliability demand, etc., is a kind of sensing data to the aftermarket service of large manufacturing enterprises. Prediction of the demand for spare parts plays a crucial role in inventory management and lifecycle quality management for the aftermarket service of large-scale manufacturing enterprises. In real-life applications, however, demand for spare parts occurs randomly and fluctuates greatly, and the demand sequence shows obvious intermittent distribution characteristics. Additionally, due to factors such as reporting mistakes made by personnel or environmental changes, the actual data of the demand for spare parts are prone to abnormal variations. It is thus hard to capture the evolutional pattern of the demand for spare parts by traditional time series forecasting methods. The reliability of prediction results is also reduced. To address these concerns, this paper proposes a tensor optimization-based robust interval prediction method of intermittent time series for the aftersales demand for spare parts. First, using the advantages of tensor decomposition to effectively mine intrinsic information from raw data, a sequence-smoothing network based on tensor decomposition and a stacked autoencoder is proposed. Tucker decomposition is applied to the hidden features of the encoder, and the obtained core tensor is reconstructed through the decoder, thus allowing us to smooth outliers in the original demand sequence. An alternating optimization algorithm is further designed to find the optimal sequence feature representation and tensor decomposition factors for the extraction of the evolutionary trend of the intermittent series. Second, an adaptive interval prediction algorithm with a dynamic update mechanism is designed to obtain point prediction values and prediction intervals for the demand sequence, thereby improving the reliability of the forecast. The proposed method is validated using the actual aftersales data from a large engineering manufacturing enterprise in China. The experimental results demonstrate that, compared with typical time series prediction methods, the proposed method can effectively grab the evolutionary trend of various intermittent series and improve the accuracy of predictions made with small-sample intermittent series. Moreover, the proposed method provides a reliable elastic prediction interval when distortion occurs in the prediction results, offering a new solution for intelligent planning decisions related to spare parts in practical maintenance.

1. Introduction

In complicated equipment manufacturing enterprises such as shield tunneling, rail transportation, and wind energy, the cost of maintaining an inventory of spare parts generally accounts for more than 60% of inventory costs [1]. Qualified inventory optimization [2] and flexible scheduling of parts [3] are critical to improve the efficiency of aftermarket service in lifecycle product management [4]. Due to various uncertain factors, such as element failure, project schedule, safe inventory level, etc., there will be an inventory shortage of spare parts in warehouses which, in turn, triggers the demand for spare parts. The demand for spare parts can then serve as a kind of sensing data to monitor the maintenance efficiency and evaluate the aftermarket service quality. Accurate prediction of the demand for spare parts plays a supporting role in intelligent inventory optimization. However, in practical operations, parts planning is often linked to new online projects or associated with the unavailability of spare parts, resulting in sporadic demand for spare parts. The data distribution of the demand for spare parts therefore has intermittent characteristics. Precisely predicting the demand of intermittent time series is challenging in spare parts management for manufacturing enterprises.

Due to the intermittent characteristics of these demand data, predicting demand relies heavily on time series prediction. Currently, time series prediction methods can be divided into three categories [5]: (1) statistical methods (e.g., exponential smoothing [6] and moving average [7]); (2) machine learning methods (e.g., support vector regression (SVR) [8], random forests (RF) [9], and LightGBM (light gradient boosting machine) [10,11]); and (3) deep learning methods (e.g., recurrent neural networks (RNN) [12] and long short-term memory (LSTM) [13]). These methods are often applicable to time series with strong periodicity and apparent trends. However, there are some challenges involved in effectively extracting evolutionary patterns from time series with strong randomness, poor continuity, and, especially, small sample sizes, leading to low prediction accuracy for intermittent series. To solve this problem, Croston [14] improved the exponential smoothing algorithm by decomposing intermittent time series into zero-interval and demand sequences. The exponential smoothing algorithm was then applied to each sequence separately, and the results were weighted to improve prediction performance. Syntetos et al. [15] further improved the Croston algorithm and designed the Syntetos–Boylan approximation (SBA) method. The SBA method introduced a bias term $(1-\alpha /2)$ to mitigate the uncertainty of intermittent distributions. In addition, some studies proposed different metrics, such as the average demand interval (ADI) and the square of the coefficient of variation (CV2) [16], to explore intermittent characteristics to better extract intrinsic information from demand sequences [17]. Another typical approach is to perform hierarchical clustering on demand sequences [18], which divides the original sequences with weak overall patterns into multiple clusters with more significant patterns. Then, different regression algorithms can be applied to the clusters for prediction. Shi et al. [19] proposed the block Hankel tensor-autoregressive integrated moving average (BHT-ARIMA) model, which uses tensor decomposition [20] to extract the intrinsic correlations among multidimensional small-sample sequences. Although the aforementioned methods have achieved certain results, they still have some limitations. First, they are mostly based on the assumption that all sequence demands have predictive value and disregard the interference of abnormal values. In fact, in actual business, due to special events such as natural disasters, emergencies, market fluctuations, and other factors, some spare parts have abnormal demand patterns, which require manual analysis for recognition. Second, demand prediction results are random and unreliable. Operations still want to obtain trustworthy results regarding demand prediction, but the existing methods fail to make valid decisions when the prediction results are distorted.

There is also a similar concept, i.e., abnormal demand forecasting, which refers to the process of forecasting and analyzing abnormal demands that may occur in the future. Guo et al. [21] utilized the passenger flow characteristics extracted by SVR into LSTM to predict abnormal passenger flow. Liu et al. [22] combined statistical learning and linear regression to build a model of the relationship between price discounts and demand for medical consumables. Li et al. [23] proposed a multi-scale radial basis function (MSRBF) network to predict subway passenger flow under special events. Nguyen et al. [24] combined LSTM with SVM to detect outliers in a demand sequence, and then used an LSTM neural network to predict the occurrence of outliers. Li et al. [25] proposed a combination model of time series and regression analysis, plus dummy variables, to predict special factors related to passenger flow. These methods aim to predict the occurrence of abnormal but meaningful events in many fields. However, in the situation described by this paper, abnormal demands are commonly found due to reporting mistakes made by personnel or environmental changes. Such abnormal demands are harmful to the prediction of the intermittent time series, since the evolutionary pattern of the demand for spare parts is disturbed. There is no value in predicting abnormal demands if such demands have no predictability. This paper is, therefore, solely devoted to predicting normal demands against the disturbances induced by abnormal demands, instead of predicting the abnormal demands.

To solve the problems mentioned above, this paper employs a tensor decomposition technique to smooth the abnormal demands in demand sequences. Tensor Tucker decomposition decomposes a tensor into a set of factor matrices and one core tensor that can describe the intrinsic information from raw data. With the decomposed forms, tensor Tucker decomposition brings the potential to extract the evolutionary trend from intermittent time series. The concept of a prediction interval is further introduced to tackle the problem of high uncertainty and less reliability in the prediction results. The methodology is as follows: First, a tensor autoencoder network is constructed, which performs tensor Tucker decomposition on the output features extracted from the hidden layers of the autoencoder. Second, the core tensor is decoded and reconstructed with an alternating optimization strategy to obtain the optimal feature representation of intermittent time series. Third, an adaptive prediction interval (API) algorithm is developed with a dynamical update mechanism to obtain point prediction values and prediction intervals, thus improving the reliability of the prediction results. Finally, the performance of the proposed method is validated using a set of real-life aftersales data from a large engineering manufacturing enterprise from China.

The contributions of this paper can be summarized as follows:

(1) An intermittent series smoothing algorithm is proposed. By integrating tensor Tucker decomposition into a stacked autoencoder network with an alternately optimizing scheme, the proposed algorithm extracts the evolutionary trend of the intermittent time series under irregular noise interference. Compared with existing time series anomaly detection methods, the proposed algorithm does not require any pre-fixed detection thresholds, and can adaptively identify outliers in the series. It is highly suitable for smoothing the outliers in the intermittent time series. Moreover, the proposed algorithm is universally applicable, e.g., the stacked autoencoder can be easily replaced by other deep models;

(2) An adaptive prediction interval algorithm is designed. Different from the existing point prediction methods using a deterministic prediction model, this algorithm incorporates the prediction intervals with the point prediction, which can match the intermittent characteristics of demand sequence for spare parts. This algorithm provides an effective solution to address the uncertainty in the demand for spare parts. According to the literature survey, there is no related research about interval prediction, specifically for demand prediction.

2. Background

2.1. Multi-Way Delay Embedding Transform

The multi-way delay embedding transform (MDT) [26] is capable of embedding low-order data into a high-dimensional space, which can be used to construct Hankel matrices or block Hankel tensors. The tensor obtained by MDT possesses the characteristics of low rank, which can smooth the original data to make them easier to train. $v={\left(v_1,\dots ,v_L\right)}^T\in R^L$ denotes a vector that is transformed into a Hankel matrix with a delay of $\tau$ through MDT. This process is referred to Hankelization of the vector. The transform process is shown in Equation (1).

$${\mathcal{H}}_{\tau }\left(vs.\right):=\left(\begin{array}{cccc}{v}_1& v_2& \cdots & v_{L-\tau +1}\\ v_2& v_3& \cdots & v_{L-\tau +2}\\ ⋮& ⋮& \ddots & ⋮\\ v_{\tau }& v_{\tau +1}& \cdots & v_L\end{array}\right)\in R^{\tau \times \left(L-\tau +1\right)}\frac{-b\pm \sqrt{b^2-4ac}}{2a}$$

(1)

First, the duplication matrix $S\in {\left\{0,1\right\}}^{\tau \times \left(L-\tau +1\right)\times L}$ is constructed with a delay of $\tau$, as shown in Equation (2).

$$S^T={\left(\begin{array}{cccc}{I}_{\tau }& & & \\ {}_{}& I_{\tau }& & \\ & & \ddots & \\ & & & I_{\tau }\end{array}\right)}_{\tau \times \tau }T$$

(2)

The vector v is transformed into a Hankel matrix denoted by ${\mathcal{H}}_{\tau }\left(vs.\right)$. The duplication matrix S is essentially a linear transformation. The vectorization expansion is shown in Equation (3):

$$vec\left({\mathcal{H}}_{\tau }\left(vs.\right)\right)=Sv,Sv\in R^{\tau \times \left(L-\tau +1\right)}$$

(3)

where $vec(·)$ expands the matrix along the column direction. The Hankel matrix through delay embedding can be shown in Equation (4).

$$\begin{array}{c}{\mathcal{H}}_{\tau }\left(vs.\right)=fol{d}_{\left(L,\tau \right)}\left(Sv\right):=v_H,\hfill \\ fol{d}_{\left(L,\tau \right)}:R^{\tau \times \left(L-\tau +1\right)}\to R^{\tau \times \left(L-\tau +1\right)}\hfill \end{array}$$

(4)

where $fol{d}_{\left(L,\tau \right)}$ folds a vector into a matrix.

The inverse transform of multi-way delay embedding for vectors can convert the data from a high-dimensional space to a lower-dimensional target space. This can be calculated using Equations (5) and (6).

$${\mathcal{H}}_{\tau }^{-1}\left(V_H\right)=S^{†}vec\left(V_H\right)$$

(5)

$$S^{†}:={\left(S^{T}S\right)}^{-1}{S}^T$$

(6)

where † is the Moore–Penrose inverse matrix [27].

2.2. Tensor Decomposition

Tensor decomposition aims to decompose high-order tensor data into low-rank matrices or vectors [28]. It is commonly applied in data compression, dimensionality reduction, feature extraction, etc. Tensor Tucker decomposition decomposes an Nth-order tensor $\chi \in {\mathbb{R}}^{I_1\times I_2\times \cdots \times I_N}$ into the product of a core tensor ${\varsigma }_t\in {\mathbb{R}}^{J_1\times J_2\times \cdots \times J_N}$ and N factor matrices $U^{\left(n\right)}\in {\mathbb{R}}^{I_n\times J_n}$, as shown in Equation (7). The factor matrices obtained from Tucker decomposition represent the principal components of the tensor’s modular expansion, while the core tensor captures the correlations between these components.

$$\chi =\varsigma \times {}_1{U}^{\left(1\right)}\times {}_2{U}^{\left(2\right)}\cdots \times {}_N{U}^{\left(N\right)}$$

(7)

where $\varsigma \times {}_1{U}^{\left(n\right)}$ is the n-mode product of the modular (n) expansion of tensor S and matrix $U^{\left(n\right)}\in {\mathbb{R}}^{I_n\times J_n}$:

$$\begin{array}{c}{\left[\varsigma \times U^{\left(n\right)}\right]}_{j_1\cdots j_{n-1}{i}_n{j}_{n+1}\cdots j_N}={\displaystyle \sum _{j_n=1}^{J_n}}{g}_{j_1\cdots j_{n-1}{i}_n{j}_{n+1}\cdots j_N}{u}_{i_n{j}_n}\hfill \\ \varsigma \times U^{\left(n\right)}\in {\mathbb{R}}^{J_1\times J_2\times \cdots \times J_N}\hfill \end{array}$$

(8)

With the above equation, one data point in the tensor can be expanded by Equation (9) as:

$$x_{i_1{i}_2\cdots i_N}=\sum _{j_1,j_2,\cdots ,j_N}^{}{g}_{j_1\cdots j_N}{u}_{i_1{j}_1}^{\left(1\right)}{u}_{i_2{j}_2}^{\left(2\right)}\cdots u_{i_3{j}_3}^{\left(3\right)}$$

(9)

Figure 1 illustrates a three-order tensor decomposed using Tucker decomposition, resulting in the product of a smaller core tensor and three factor matrices.

2.3. LightGBM

LightGBM, developed by Microsoft engineers [10], is a gradient-boosting framework based on decision trees. By utilizing feature parallelism during its training process, LightGBM assigns discrete features to multiple bins and constructs decision trees using histogram algorithms, which provides a quick and efficient training mechanism for LightGBM. Additionally, LightGBM employs a sparse feature algorithm to significantly reduce memory consumption, which is suitable for training with extremely large-scale data. Moreover, LightGBM utilizes a leaf-wise tree growth strategy, which leads to faster convergence and higher accuracy compared to traditional level-wise strategies.

3. Proposed Method

In this section, a tensor optimization-based robust interval prediction method of intermittent time series is presented. This method consists of two parts: (1) a sequence smoothing network based on tensor decomposition and a stacked autoencoder, which aims to smooth anomalous demand values in the original sequence and to extract the evolutionary trend of demand as well; and (2) an adaptive prediction interval algorithm, which aims to construct a reliable prediction interval to avoid the oversupply risk or inventory shortage caused by inaccurate predictions. In this method, the role of tensor Tucker decomposition is: (1) extracting core tensors from the demand sequence for representing the evolutionary trend; and (2) smoothing the outliers in the sequence. The negative interference of anomalous demands can then be effectively reduced under an unsupervised mode. Moreover, training a LightGBM model with core tensors can better mine the trend information from the demand sequence. Tensor Tucker decomposition is believed suitable for intermittent time series forecasting with small samples.

The flowchart of the proposed method is shown in Figure 2. First, the hidden features are extracted from the original data using a stacked autoencoder, then tensor Tucker decomposition is performed on the hidden features to obtain the core tensors. An alternating optimization scheme is designed to obtain the optimal core tensors by alternately updating the autoencoder parameters and tensor factor matrices. Second, an adaptive interval prediction algorithm is constructed. The interval is calculated using the predicted values and prediction residuals from LightGBM estimators. Finally, a dynamic update mechanism is used to adjust the width of prediction interval. The detailed implementation will be presented as follows.

3.1. Sequence Smoothing Network

Denote the demand sequence by $X=\left\{x_1,x_2,\dots ,x_m\right\}\in R^{m\times n}$, where m indicates the sample number and n represents the time dimension. The sequence can be encoded and mapped into a hidden layer, as shown in Equation (10).

$$h=f(Wx+b)$$

(10)

where h represents the hidden layer features, $f(·)$ is the activation function of the encoding layer, and $W\in {\mathbb{R}}^r$ and $b\in {\mathbb{R}}^r$ are the weight matrix and bias vector of the encoding layer, respectively. The obtained feature set is denoted by $H=\{h_1^{},h_2^{},\cdots ,h_m^{}\}\in R^{m\times l}$, where l represents the dimension of the deep features.

In order to adequately represent the temporal information between samples in the feature set H, an MDT with the operations of multi-linear duplication and multi-way folding is employed to transform the original sequence to a three-order tensor. This process can also reduce noise disturbance. Denote by $H=\{h_1^{},h_2^{},\cdots ,h_m^{}\}\in R^{m\times l}$ the tensor of $\mathbb{X}=\{{\mathcal{X}}_1,\dots ,{\mathcal{X}}_{m-\tau +1}\}\subseteq {\mathbb{R}}^{l\times \tau \times (m-\tau +1)}$, then the MDT for $\mathbb{X}$ can be defined as:

$${\mathcal{X}}_i={\mathcal{H}}_{\tau }\left(H\right)=Fol{d}_{(m,\tau )}(H\times {}_1{S}_1\times \cdots \times {}_{m-\tau +1}{S}_{m-\tau +1})$$

(11)

where $\tau$ and m represent the time window size and sample length, respectively, and S is a duplication matrix. In this paper, $\tau$ is set to 6.

With the tensors constructed by MDT, the Tucker decomposition technique is introduced to obtain the three-order core tensor $\mathcal{G}=\{g_1^1,\dots ,g_{{}_{m-\tau +1}}^l\}\in R^{l\times \tau \times (m-\tau +1)}$ that represents the essential information of the sequence:

$$\begin{array}{c}\mathcal{G}=\mathbb{X}\times {}_1{U^{\left(1\right)}}^{\mathrm{T}}\times {}_2{U^{\left(2\right)}}^{\mathrm{T}}\times \cdots \times {}_V{U^{\left(v\right)}}^{\mathrm{T}}\hfill \\ s.t.{U^{\left(v\right)}}^{\mathrm{T}}{U}^{\left(v\right)}=I,v=1,\cdots ,V\hfill \end{array}$$

(12)

where ${\left\{U^{\left(v\right)}\right\}}_{v=1}^V$ is the projection matrix and can maximally preserve the temporal continuity between core tensors, where the superscript v represents the tensor dimension, and $m-\tau +1$ represents the sample length after reconstruction. Equation (12) means that the decomposition result consists of a core tensor and a series of factor matrices. Obviously, the core tensor $\mathcal{G}$ contains the intrinsic information of an evolutionary trend. By optimizing ${\left\{U^{\left(v\right)}\right\}}_{v=1}^V$, the inherent correlation between feature sequences can be sufficiently captured, while noise interference can also be reduced. The loss of tensor reconstruction can be calculated as:

$$\begin{array}{c}{\mathcal{L}}_{Tensor}=min{∥\mathbb{X}-\widehat{\mathbb{X}}∥}_F^2\hfill \end{array}$$

(13)

where $\mathbb{X}$ is the original tensor, and $\mathbb{X}=\mathcal{G}\times {}_1{U}^{\left(1\right)}\times {}_2{U}^{\left(2\right)}\times \cdots \times {}_H{U}^{\left(H\right)}$ is the tensor reconstructed from the factor matrix ${\left\{U^{\left(v\right)}\right\}}_{v=1}^V$.

To implement the decoding operation, it is necessary to transform the core tensor $\mathcal{G}$ back to the original input space. Here, the inverse MDT transform [6] is applied to $\mathcal{G}$ by reversing the transformation along the time dimension to obtain a second-order core tensor ${\mathcal{G}}^{\prime }=\{g_1^{\prime },\dots ,g_m^{\prime }\}\in R^{m\times l}$, as shown in Equation (14). The core tensor ${\mathcal{G}}^{\prime }$ is then used as the input of decoding to update the network.

$${\mathcal{G}}^{\prime }=\mathcal{H}\left(\mathcal{G}\right)=Unfol{d}_{(m,\tau )}\left(\mathcal{G}\right)\times {}_1{S_1}^{†}\times \cdots \times {}_{m-\tau +1}{S}_{m-\tau +1}^{†}$$

(14)

where † is the Moore–Penrose pseudo-inverse.

Through decoding ${\mathcal{G}}^{\prime }$, the reconstructed data $\widehat{x}$ can be obtained as follows:

$$\widehat{x}=f^{*}(W^{*}{g}^{\prime }+b^{*})$$

(15)

where $f^{*}(·)$ is the activation function of the decoding layer, $W^{*}\in {\mathbb{R}}^{n\times r}$ is the weight matrix of the decoding layer, and $b^{*}\in {\mathbb{R}}^n$ is the bias vector of the decoding layer. Consequently, the reconstruction loss of the autoencoder can be calculated as:

$$\underset{}{{\mathcal{L}}_{AE}=}\frac{1}{m}\sum _{i=1}^m\frac{1}{2}{∥\widehat{x}-x∥}^2+\frac{\lambda }{2}({{∥W∥}_F}^2+{{∥W^{*}∥}_F}^2)$$

(16)

where $\lambda$ is the weight decay parameter and $∥·∥$ is the Frobenius norm.

Based on the aforementioned analysis, the whole loss function is:

$${\mathcal{L}}_{loss}=\sum _{i=1}^M{\mathcal{L}}_{AE}+\eta {\mathcal{L}}_{Tensor}$$

(17)

Minimizing Equation (17) can smooth anomalous demands in the demand sequence and extract the evolutionary trend of the demand for spare parts. The key idea of this process is to reduce the significant deviations of anomalous demands, making them suitable for intermittent sequence anomaly detection. The crucial aspect of this process lies in utilizing the core tensor to represent the evolutionary trend. As Equation (12) indicates, ${\left\{U^{\left(v\right)}\right\}}_{v=1}^V$ is randomly initialized. Therefore, the optimization of tensor decomposition in minimizing Equation (17) is required to obtain the optimal representation of the core tensors.

3.2. Alternating Optimization Scheme

Minimization of Equation (17) includes the optimization of ${\mathcal{L}}_{AE}$ and ${\mathcal{L}}_{Tensor}$. The ${\mathcal{L}}_{AE}$ can be solved using a stochastic gradient descent (SGD) strategy [29]. However, SGD cannot be directly applied to tensor decomposition. Therefore, this paper adopts an alternating optimization strategy: fix ${\left\{U^{\left(v\right)}\right\}}_{v=1}^V$, and update $W,W^{*}$; then fix the updated $W,W^{*}$, and update ${\left\{U^{\left(v\right)}\right\}}_{v=1}^V$. These two steps are performed alternately until convergence. It should be noted that, since the number of tensor optimization iterations is typically smaller than the number of $W,W^{*}$ updates, updating ${\left\{U^{\left(v\right)}\right\}}_{v=1}^V$ is set to stop when the difference between two consecutive tensor optimizations of ${\left\{U^{\left(v\right)}\right\}}_{v=1}^V$ is less than a specific threshold. $W,W^{*}$ is updated until the model converges. With the initialized W and $W^{*}$, the specific optimization process is as follows:

(1) Fix ${\left\{U^{\left(v\right)}\right\}}_{v=1}^V$, and update $W,W^{*}$;

1. Encoding stage: the partial derivatives of ${\mathcal{L}}_{AE}$ with respect to the parameters are:

$$\left\{\begin{array}{c}\frac{\partial {\mathcal{L}}_{AE}}{\partial W}=\frac{2}{M}{\displaystyle \sum _{m=1}^M}({\widehat{x}}^{\left(m\right)}-x^{\left(m\right)})·\frac{\partial {({\widehat{x}}^{\left(m\right)}-x^{\left(m\right)})}^T}{\partial W}+2\lambda \frac{\partial ({{∥W∥}_F}^2+{{∥W^{*}∥}_F}^2)}{\partial W}\hfill \\ \frac{\partial {\mathcal{L}}_{AE}}{\partial b}=\frac{2}{M}{\displaystyle \sum _{m=1}^M}({\widehat{x}}^{\left(m\right)}-x^{\left(m\right)})·\frac{\partial {({\widehat{x}}^{\left(m\right)}-x^{\left(m\right)})}^T}{\partial b}\hfill \end{array}\right.$$

(18)

The partial derivatives corresponding to the error term on each sample can be obtained by:

$$\left\{\begin{array}{c}\frac{\partial ({\widehat{x}}^{\left(m\right)}-x^{\left(m\right)})}{\partial W}={\displaystyle \sum _{m=1}^M}\frac{\partial {\left({\widehat{x}}^{\left(m\right)}\right)}^T}{\partial W}=f^{\prime }·diag\left(X^{\left(m\right)}\right)\in {\mathbb{R}}^{1\times r}\hfill \\ \frac{\partial ({\widehat{x}}^{\left(m\right)}-x^{\left(m\right)})}{\partial b}={\displaystyle \sum _{m=1}^M}\frac{\partial {\left({\widehat{x}}^{\left(m\right)}\right)}^T}{\partial b}=f^{\prime }·1_r\in {\mathbb{R}}^{1\times 1}\hfill \end{array}\right.$$

(19)

where “·” represents the dot product operator, $diag(·)$ is the diagonal matrix, $1_r$ is the unit column vector of size r, and $f^{\prime }$ is the derivative of the activation function;

2. Decoding stage: the partial derivatives of ${\mathcal{L}}_{AE}$ with respect to the parameters are:

$$\left\{\begin{array}{c}\frac{\partial {\mathcal{L}}_{AE}}{\partial W^{*}}=\frac{2}{M}{\displaystyle \sum _{m=1}^M}({\widehat{x}}^{\left(m\right)}-x^{\left(m\right)})·\frac{\partial {({\widehat{x}}^{\left(m\right)}-x^{\left(m\right)})}^T}{\partial W^{*}}+2\lambda \frac{\partial ({{∥W∥}_F}^2+{{∥W^{*}∥}_F}^2)}{\partial W^{*}}\hfill \\ \frac{\partial {\mathcal{L}}_{AE}}{\partial b^{*}}=\frac{2}{M}{\displaystyle \sum _{m=1}^M}({\widehat{x}}^{\left(m\right)}-x^{\left(m\right)})·\frac{\partial {({\widehat{x}}^{\left(m\right)}-x^{\left(m\right)})}^T}{\partial b^{*}}\hfill \end{array}\right.$$

(20)

For easy analysis, the error propagation term (i.e., the derivative of the error term with respect to the hidden layer output) is briefly denoted by:

$${\delta }_{H^{\left(m\right)}}=\left(\frac{\partial ({\widehat{x}}^{\left(m\right)}-x^{\left(m\right)})}{\partial H^{\left(m\right)}}\right)$$

(21)

Furthermore, following the chain rule, we have:

$$\left\{\begin{array}{c}\frac{\partial ({\widehat{x}}^{\left(m\right)}-x^{\left(m\right)})}{\partial W^{*}}={\delta }_{H^{\left(m\right)}}·\frac{\partial H^{\left(m\right)}}{\partial W^{*}}={\delta }_{H^{\left(m\right)}}·{(f^{\prime }\odot diag\left(H^{\left(m\right)}\right))}^T\hfill \\ \frac{\partial ({\widehat{x}}^{\left(m\right)}-x^{\left(m\right)})}{\partial b^{*}}={\delta }_{H^{\left(m\right)}}·\frac{\partial H^{\left(m\right)}}{\partial b^{*}}={\delta }_{H^{\left(m\right)}}·{(f^{\prime }\odot 1_r)}^T\hfill \end{array}\right.$$

(22)

where “·” and $diag(·)$ have the same meaning as in the encoding stage;

3. Substitute Equations (19) and (22) into Equations (18) and (20), respectively, and the model parameters can be updated as:

$$\left\{\begin{array}{c}{W}^{(l+1)}=W^{\left(l\right)}-\xi ·\frac{\partial {\mathcal{L}}_{AE}}{\partial W}\left|{}_{W=W^{\left(l\right)}}\right.\hfill \\ b^{(l+1)}=b^{\left(l\right)}-\xi ·\frac{\partial {\mathcal{L}}_{AE}}{\partial b}\left|{}_{b=b^{\left(l\right)}}\right.\hfill \\ {W^{*}}^{(l+1)}={W^{*}}^{\left(l\right)}-\xi ·\frac{\partial {\mathcal{L}}_{AE}}{\partial W^{*}}\left|{}_{W^{*}={W^{*}}^{\left(l\right)}}\right.\hfill \\ {b^{*}}^{(l+1)}={b^{*}}^{\left(l\right)}-\xi ·\frac{\partial {\mathcal{L}}_{AE}}{\partial b^{*}}\left|{}_{b^{*}={b^{*}}^{\left(l\right)}}\right.\hfill \end{array}\right.$$

(23)

where $\xi$ is the learning rate.

(2) Fix $W,W^{*}$ and update ${\left\{U^{\left(v\right)}\right\}}_{v=1}^V$.

Updating ${\left\{U^{\left(v\right)}\right\}}_{v=1}^V$ can be achieved by minimizing the tensor reconstruction error as shown in Equation (13). Since ${\left\{U^{\left(v\right)}\right\}}_{v=1}^V$ is an orthogonal matrix, Equation (13) can be rewritten as follows:

$$\begin{array}{c}{∥\mathbb{X}-\widehat{\mathbb{X}}∥}_F^2={∥vec\left(\mathbb{X}\right)-(U^{\left(V\right)}\otimes U^{(V-1)}\otimes \cdots \otimes U^{\left(1\right)})·vec\left(\mathcal{G}\right)∥}_F^2\hfill \\ ={∥vec\left(\mathbb{X}\right)-(U^{\left(V\right)}\otimes U^{(V-1)}\otimes \cdots \otimes U^{\left(1\right)})·{(U^{\left(V\right)}\otimes U^{(V-1)}\otimes \cdots \otimes U^{\left(1\right)})}^T·vec\left(\mathbb{X}\right)∥}_F^2\hfill \\ ={∥vec\left(\mathbb{X}\right)∥}_F^2-{∥U^T·vec\left(\mathbb{X}\right)∥}_F^2\hfill \end{array}$$

(24)

where $U^T=U^{\left(H\right)}\otimes U^{(H-1)}\otimes \cdots \otimes U^{\left(1\right)}$. To minimize Equation (24), we should maximize the following equation:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {∥U^T·vec\left(\mathbb{X}\right)∥}_F^2=\left\{\begin{array}{c}{∥U^{{\left(1\right)}^T}·A^{\left(1\right)}·(U^{\left(H\right)}\otimes U^{(H-1)}\otimes \cdots \otimes U^{\left(2\right)})∥}_F^2\hfill \\ {∥U^{{\left(2\right)}^T}·A^{\left(2\right)}·(U^{\left(H\right)}\otimes U^{(H-1)}\otimes \cdots \otimes U^{\left(1\right)})∥}_F^2\hfill \\ \cdots \hfill \\ {∥U^{{\left(H\right)}^T}·A^{\left(H\right)}·(U^{(H-1)}\otimes U^{(H-2)}\otimes \cdots \otimes U^{\left(1\right)})∥}_F^2\hfill \end{array}\right.\end{array}\end{array}$$

(25)

where $A^{\left(i\right)}$ is the unfolding matrix of the tensor $\mathbb{X}$ along the i-th dimension. The alternating least squares method [30] can be used to solve it, where each factor matrix in different directions is fixed sequentially to find the least squares solution. This process is shown in Equation (26):

$$U_{k+1}^{\left(v\right)}=\underset{{\left\{U^{\left(v\right)}\right\}}_{v=1}^V}{argmin}{∥\widehat{\mathbb{X}}-\mathcal{G}\times {}_M{U}_k^{{}^{\left(V\right)}}\times {}_{M-1}\cdots \times {}_1{U}_k^{{}^{\left(1\right)}}∥}_F^2$$

(26)

3.3. Adaptive Prediction Intervals (API)

After obtaining the optimal feature representation, this paper designs the API algorithm to generate the point prediction values and reliable prediction intervals. The details of the API algorithm are presented as follows. The API algorithm consists of two stages: training and prediction. At the training stage, a fixed number of LightGBM estimators are fitted from a subset of training data. Then, the predicted value from all the LightGBM estimators is aggregated using the leave-one-out (LOO) strategy to generate LOO prediction factors and residuals. At the prediction stage, the API algorithm calculates the LOO predicted value by averaging the prediction factors from each test sample. With the predicted values, the center of the prediction interval is determined. The prediction intervals are then established using the LOO residuals. The width of prediction interval is updated using a dynamic updating strategy. The specific implementation steps are as follows:

First, a LightGBM model f is trained using the training samples ${\left\{(x_t,y_t)\right\}}_{t=1}^T$, and the prediction interval at the t-th time step is calculated as:

$${\widehat{C}}_t^{\alpha }=[{\widehat{f}}_{-t}\left(x_t\right)+q_{\widehat{\beta }},{\widehat{f}}_{-t}\left(x_t\right)+q_{(1-\alpha +\widehat{\beta })}]$$

(27)

where $\alpha$ is the significance factor, ${\widehat{f}}_{-t}$ represents the t-th estimator of f, $q_{\widehat{\beta }}$ is the $\widehat{\beta }$-quantile on ${\left\{{\widehat{\varsigma }}_i\right\}}_{i=t-1}^{t-T}$, and $q_{(1-\alpha +\widehat{\beta })}$ is the $(1-\alpha +\widehat{\beta })$-quantile on ${\left\{{\widehat{\varsigma }}_i\right\}}_{i=t-1}^{t-T}$. The LOO prediction residual ${\widehat{\varsigma }}_i$ and $\widehat{\beta }$ are defined as follows:

$${\widehat{\varsigma }}_i=y_i-{\widehat{f}}_{-t}\left(x_t\right)$$

(28)

$$\widehat{\beta }=\underset{\beta \in [0,\alpha ]}{argmin}({\widehat{f}}_{-t}\left(x_t\right)+q_{(1-\alpha +\widehat{\beta })}-{\widehat{f}}_{-t}\left(x_t\right)+q_{\widehat{\beta }})$$

(29)

Second, the interval center is ${\widehat{f}}_{-t}\left(x_t\right)$, and the interval width is the difference between the $(1-\alpha +\widehat{\beta })$ and $\widehat{\beta }$-quantiles on the past T residuals.

However, the obtained interval values cannot adequately address the problem of large fluctuations of demands for different spare parts. This section proposes an adaptive update mechanism to improve the reliability of the prediction interval. The specific step is as follows: First, each demand sequence is divided into a zero-interval sequence (i.e., the gap between two subsequent demands occur) and a non-zero sequence (i.e., the real demand values). The squared coefficient of variation is then calculated for the two sequences, which are denoted by $c{v}_1$ and $c{v}_2$, respectively. A larger coefficient indicates greater sequence fluctuation, thus requiring a larger interval width for prediction, and vice versa. Second, initialize the interval width parameter $\alpha =0.1$ and update $\alpha$ by:

$$\alpha =\alpha +c{v}_1+c{v}_2$$

(30)

The interval update mechanism, as shown in Equation (30), can improve the rationality and reliability of the prediction interval for demand sequences with different volatility and intermittency characteristics. Obviously, different values of interval widths will directly affect the decision of inventory management, e.g., spare part coverage rate. The spare part coverage rate, reflected by the interval coverage rate in this paper, is defined as the ratio of the number of demands covered by the interval to the total number of demands in the sequence. The mechanism begins by setting an initial interval with a larger width, and then reduces the width to meet the fluctuations in the demand sequence. The update process is stopped when the interval width reaches a certain threshold. In this study, the threshold of the interval width is determined just by the interval coverage rate. The setting of the coverage rate relies heavily on business logic. Different kinds of maintenance tasks or enterprises have different requirement for the setting. In this experiment, we received help from the maintenance engineers from our cooperative enterprise and set the lower limit of coverage rate to be 60%. We also observe that, with this threshold, the prediction results become stable, which indicates that the threshold runs well.

4. Validation Results and Analysis

4.1. Experimental Settings

The experimental data for validation are the real-life demand data of spare parts from a large engineering manufacturing enterprise from China. The data set contains 75 sequences of different spare parts, in which each sequence covers 34 months from November 2018 to August 2021. The training data are the demand values in the first 33 months, while the data of the last month are for the test. The enterprise usually has one month in advance to make inventory plans and implement the allocation. Therefore, in this experiment we mainly focus on the prediction value of the last month in the whole sequence. To visualize the intermittent characteristic of the demand data, we randomly select two spare parts and illustrate their demand sequences in Figure 3.

For a fair comparison, we introduce six representative methods of time series forecasting in this experiment. These six methods are applicable for different kinds of time series. It is worth noting that BHT-ARIMA, which also adopts tensor decomposition and joint optimization strategy, can be viewed as the SOTA method for intermittent time series forecasting. See

Table 1. Six prediction methods for comparison.

Method Name	Type	Implementation
Croston [14]	Intermittent time series forecasting	Exponential smoothing
BHT-ARIMA [19]		Joint optimization with ARIMA and tensor decomposition
ARIMA [31]	One-dimensional time series forecasting	Autoregressive modeling with moving average
Random Forest [9]	Multidimensional time series forecasting	Ensemble prediction
LightGBM [10]		Decision trees using histogram algorithm
LSTM [13]	Temporal deep learning method	Modeling long-term dependencies with memory units

for detailed information.

In the experiment, the parameters $\alpha$ and $\beta$ in the Croston method are set in the range of [0.13, 0.26], and the step size is set to 0.07. The parameters p and q of BHT-ARIMA are set to 1 and 3, respectively. For ARIMA, the parameters p and q are set to 2 and 3, respectively. For Random Forest, a grid search strategy is adopted to select the optimal parameters, with the parameters determined as max_depth = 5, n_estimators = 30, learning = 0.05, and num_leaves = 20. The LSTM hidden layer is set to 20, and the learning rate is 0.001.

4.2. Result Analysis

In this experiment, the demand data of the first 33 months are used as training and the demand data of the 34th month is for the test. The predicted values and prediction intervals obtained by the proposed method are shown in Figure 4.

As suggested by the enterprise’s engineers, a prediction result within a range of ±30% around the real values can be regarded correct. From Figure 4, our method works well when dealing with high-demand data. The reason is that, besides the unsupervised feature extraction capability of the autoencoder network, the core tensors can represent the evolutionary trend of the raw demand data well. More importantly, when the prediction results are heavily biased, the prediction intervals obtained by the proposed method can effectively cover the true value. It is clear that the prediction results by the proposed method can improve the decision-making ability of enterprises in inventory management.

Figure 5 shows the comparison of different prediction methods on the demand sequences with different degrees of intermittent distribution.

From Figure 5, the zero-interval of the No.1437 sequence is relatively stable, so the Croston method can achieve good results. However, the Croston method performs worse than the other methods on the sequences No.0392 and No.0052 with large demand fluctuations. It indicates that Croston is only suitable for the intermittent series with a stationary zero-interval. The other conventional methods also have similar effects. These methods all have certain limitations and are only applicable to certain types of intermittent time series. Our method outperforms the other methods. The key part of our method is the introduction of tensor decomposition for smoothing anomalous demands, which is beneficial for the effective extraction of evolutionary trends from intermittent time series. With no surprise, our method shows superiority in the forecasting of intermittent time series with different distribution characteristics.

We further evaluate the effectiveness of the designed prediction interval. As stated above, a reliable prediction interval is able to provide a reference for inventory management, even if the point prediction is heavily biased or distorted. Figure 6 shows the prediction interval and predicted value of the demand sequences No. 0392, 0602, 0348, and 3082.

From Figure 6, the point prediction values of the sequences No. 0392 and No. 0602 are more accurate than the others. The prediction intervals also cover the real values, while the interval range remains small. On the contrary, the demand number in the sequences No. 0348 and No. 3082 is smaller, while the deviation of point prediction results is rather large. The prediction interval can cover the real value, so the enterprise can obtain precise information for the spare parts in this interval and make a more reliable plan on inventory management.

For numerical comparison, we introduce three commonly-used error metrics: MAE, RMSE, and RMSSE. The prediction errors obtained by different methods are listed in

Table 2. Numerical comparison of different methods in terms of three error metrics on all 75 demand sequences.

Algorithm	MAE	RMSE	RMSSE
Random Forest	1.85	2.77	0.79
ARIMA	1.99	4.40	0.73
Croston	1.77	2.79	0.84
LightGBM	1.85	2.80	0.78
LSTM	1.74	3.20	0.91
BHT-ARIMA	1.67	2.73	0.71
Our method	1.64	2.57	0.58

. Our method obtains the lowest prediction error in terms of all three metrics. With these results, our method is believed to provide a new reliable solution for enterprises to realize inventory management and parts scheduling.

The mainstream method in this field is the demand prediction method based on deep learning models. One advantage of the proposed method, compared to the current mainstream deep learning model, is introduction of tensor decomposition. By applying tensor decomposition, the core tensor can be extracted from the demand sequences. This effectively mitigates the impact of outliers in the demand sequences and diminishes their interference on the prediction results. Additionally, by utilizing LightGBM prediction with the core tensor, the proposed method can better capture the demand trend information in the demand sequences. The proposed method is then more suitable than deep learning methods for predicting intermittent time series with small samples.

We must point out the potential disadvantages of tensor Tucker decomposition. Actually, the process of tensor decomposition is rather computationally expensive. In the experiment, each round of tensor decomposition needs about 0.6 s, while the alternating optimization needs 3.07 s on average. The total cost of the proposed method is 30 s per round. In contrast, LSTM takes an average of 2.15 s per round, and the other shallow models like Croston, Random Forest, and ARIMA needs much less time. The high cost raises the potential risks of the proposed method for applications that require high real-time performance.

5. Conclusions

In this paper, a new tensor optimization-based robust interval prediction method of intermittent time series is proposed to forecast the demand for spare parts. With the introduction of tensor decomposition, the proposed method can smooth the anomalous demand while preserving the intrinsic information of evolutionary trend from the demand sequences. Moreover, to tackle the distorted or biased prediction results of intermittent time series, the API algorithm is designed to transform traditional point prediction to adaptive interval prediction. The proposed method is able to provide a trustworthy interval for the prediction results, and can accurately reflect the uncertainty of the prediction results.

This study is fundamental to develop the effect of inventory management. As a typical extended application, the proposed method is critical to update the current safety stock model to a dynamic version by integrating the evolutionary trend of demand for spare parts. For instance, the prediction results from this paper can adjust the upper bound of safety stock to match the fluctuations in the actual demands. We believe this operation can decrease the cost level while keeping maintenance.

Another interesting issue is the correlation among the demand sequences for different spare parts. The demand for two spare parts is probably correlated due to many factors such as project cycle, climatic influence, failure probability, etc. As proven by many prior works, the correlation information between two time series is valuable. Especially for intermittent time series prediction, analyzing the correlation among some sequences is believed to improve the predictability, which is critical for intermittent time series. The joint prediction of two or more sequences is able to enhance the prediction performance in terms of accuracy and numerical stability. We plan to explore the utilization of correlation information in the future work. A feasible implementation is running adaptive clustering algorithms with profile coefficients to determine the proper sequence clusters before the prediction, which is expected to improve the predictability. We can further adopt multi-output regression or multi-task learning algorithms to achieve joint prediction on the sequences in the same cluster. We plan to explore the utilization of correlation information in the future work.

In our future work, transfer learning will be studied for the prediction of demand for spare parts. In actual enterprises, the demand data are usually insufficient, especially for newly deployed equipment. Considering the distinction between different types of equipment, we plan to build a transfer learning model to incorporate the evolutionary information of the demand from available equipment. The reliability of the prediction interval is also an interesting work. It is necessary to find a reliable method to converge the prediction interval into a prediction point with a high confidence level.