The ARMA model is usually suitable for short-term forecasts of time series data, and is widely applied in business, economics, engineering and other areas (Box et al., 2008). The ARMA is usually formulated as the following equations: ρ ( B ) y t = τ ( B ) a t , (1) Where, ρ ( B ) = 1 − ρ 1 B − ρ 2 B 2 − ⋯ − ρ p B p , τ ( B ) = τ 0 + τ 1 B + τ 2 B 2 + ⋯ + τ q B q ,

y t is the output of the model at time t , and a t is random shocks, such as white noise of Gaussian distribution. B is defined as the backward shift operator, i.e., y t − 1 = B y t . When ρ ( B ) = 1 , the ARMA model can be degenerated into the MA model. Similarly, when τ ( B ) = 1 , the model will become AR model. It is worth noting that AR and MA models are both special cases of an ARMA model. [ ρ 1 ⋯ ρ p , τ 0 , τ 1 ⋯ τ q ] represents the weight value of the items corresponding to [ y t − 1 ⋯ y t − p , a t , a t − 1 ⋯ a t − q ] . Parameters, [ ρ 1 ⋯ ρ p , τ 0 , τ 1 ⋯ τ q ] , are unknown and need to be estimated by using the collected data. In this paper, RLS and RELS algorithms are used for parameter identifications. The RLS is used to minimize the cost function as follows: J ( θ ) = ∑ i = 1 t [ y i − φ i T θ ] 2 , (2) Where, φ i = [ y i − 1 , ⋯ , y i − p , a i , a i − 1 , ⋯ , a i − q ] T , θ = [ ρ 1 , ⋯ , ρ p , τ 0 , τ 1 , ⋯ , τ q ] T ,

The RLS algorithm for estimating parameter θ can be expressed as (Ding, 2010): θ ^ t = θ ^ t − 1 + P t φ t [ y t − φ t T θ ^ t − 1 ] , (3) P t = P t − 1 + P t − 1 φ t φ t T P t − 1 1 + φ t T P t − 1 φ t , (4) Where P 0 = p 0 I p + q + 1 , P i ( i = 1 , ⋯ , t ) ∈ R ( p + q + 1 ) × ( p + q + 1 ) are the covariance matrix and I p + q + 1 is an identity matrix of the order p + q + 1 . p 0 is assumed as a large positive number, e.g., p 0 = 10 6 . θ ^ t is the estimated value of θ at time t . On the basis of RLS, RELS additionally take into account innovation [ e t − 1 , ⋯ , e t − r ] . The innovation in this study indicates the difference between the real value and the predictive, where e = y − φ T θ . The general form of ARMA-RELS can be extended as shown in Eq. 5. ρ ( B ) y t = τ ( B ) a t + c ( B ) e t , (5) Where, c ( B ) = c 1 B + c 2 B 2 + ⋯ + c r B r ,

The sets of other parameters also need to be updated, φ i = [ y i − 1 , ⋯ , y i − p , a i , a i − 1 , ⋯ , a i − q , e i − 1 , ⋯ , e i − r ] T , θ = [ ρ 1 , ⋯ , ρ p , τ 0 , τ 1 , ⋯ , τ q , c 1 , ⋯ , c r ] T ,

P 0 = p 0 I p + q + r + 1 , P i ( i = 1 , ⋯ , t ) ∈ R ( p + q + r + 1 ) × ( p + q + r + 1 ) , while the recursive form that is used to estimate the parameter θ remains unchanged.

However, aforementioned ARMA models omit the interference of exogenous variables on the prediction results. The ARMAX method was proposed to tackle this nuisance. The ARMAX model takes into account not only the effects from the historical series of output itself, but additionally the effects of exogenous inputs. The general expressions for ARMAX-RLS and ARMAX-RELS can be presented as Eqs 6, 7, respectively. ρ ( B ) y t = τ ( B ) a t + σ ( B ) u t , (6) ρ ( B ) y t = τ ( B ) a t + σ ( B ) u t + c ( B ) e t , (7) Where, σ ( B ) = σ 0 + σ 1 B + σ 2 B 2 + ⋯ + σ j B j ,

u t is the exogenous inputs set. This set of exogenous inputs is possibly multidimensional, whose dimensionality depends on the number of input variables selected. The method used in this research to determine the input variables is described in Section 2.3. On the other hand, the formulations of the two parameter estimation algorithms, RLS and RELS, for the ARMAX model are basically consistent with the ARMA model, respectively. However, mainly due to the effect of the additional u t , some parameters are updated accordingly on the original basis as follows: φ i = [ y i − 1 , ⋯ , y i − p , a i , a i − 1 , ⋯ , a i − q , u i , u i − 1 , ⋯ , u i − j , ] T , θ = [ ρ 1 , ⋯ , ρ p , τ 0 , τ 1 , ⋯ , τ q , σ 0 , σ 1 , ⋯ , σ j ] T ,

P 0 = p 0 I p + q + j + 2 , P i ( i = 1 , ⋯ , t ) ∈ R ( p + q + j + 2 ) × ( p + q + j + 2 ) for RLS and, φ i = [ y i − 1 , ⋯ , y i − p , a i , a i − 1 , ⋯ , a i − q , e i − 1 , ⋯ , e i − r , u i , u i − 1 , ⋯ , u i − j , ] T , θ = [ ρ 1 , ⋯ , ρ p , τ 0 , τ 1 , ⋯ , τ q , c 1 , ⋯ , c r , σ 0 , σ 1 , ⋯ , σ j ] T ,

P 0 = p 0 I p + q + r + j + 2 , P i ( i = 1 , ⋯ , t ) ∈ R ( p + q + r + j + 2 ) × ( p + q + r + j + 2 ) for RELS. Then, p , q , r and j are calculated by the Akaike Information Criterion method as shown in Eq. 8. A I C ( p , q , r , j ) = − 2 ⁡ ln ( L ) + 2 Q , (8) Where Q is the number of parameters, L is the likelihood function. The optimal order of model is the { p , q , r , j } value satisfying the minimum A I C . Accordingly, due to the introduction of lag r , iterative k-step ahead prediction can be formulated as Eqs 9, 10. ρ ^ ( B ) y ^ t + k | t = τ ^ ( B ) a t + k | t + σ ^ ( B ) u t + k | t + ∑ i = k r c ^ i e t + k − i , 1 ≤ k ≤ r , (9) ρ ^ ( B ) y ^ t + k | t = τ ^ ( B ) a t + k | t + σ ^ ( B ) u t + k | t , k ≥ r , (10) Where [ ρ ^ , τ ^ , σ ^ , c ^ i , y ^ t + k | t ] are the corresponding estimated values of [ ρ , τ , σ , c i , y t + k ] . If k ≤ 0 , y ^ t + k | t = y t + k , meaning that the estimated output of the prediction model is equal to the real value.

Both RLS and RELS algorithms are simple but powerful to estimate unknown parameters without needing to calculate matrix inversion during iterative learning. These make them suitable for online model identification. RELS is actually a direct extension of RLS, aiming to reduce the influence of colored noise by adding residuals to the information vector φ and the parameter vector θ . Compared to the standard least squares method, RLS and RELS algorithm improve the identification performance of time series model at the expense of higher computational complexity.

Time series data derived from real industrial processes usually exhibit strong nonlinearity and high dynamics, which renders the monitoring of such data unsuitable if using linear models. Therefore, nonlinear methods based on neural networks are highly recommended for modeling such dataset. The standard NARNN is formulated as follows: y ^ t = N ( y t − 1 , y t − 2 , ⋯ , y t − d y ) , (11) Where y ^ t is the estimation of the output by a specific neural network at the t moment, [ y t − 1 , y t − 2 , ⋯ , y t − d y ] is the time series dataset, d y is the maximum output-memory order and N ( ∙ ) means a specific neural network. The distinction between NARNN and ordinary neural networks for multi-step prediction is that several observed data have to be replaced by the estimate of the network, so that Eq. 11 can also be reformulated as follows: y ^ t = N ( y ^ t − 1 , y ^ t − 2 , ⋯ , y ^ t − s , y t − s − 1 , ⋯ , y t − d y ) , (12) Where [ y ^ t − 1 , y ^ t − 2 , ⋯ , y ^ t − s ] are estimates of the output over the time period from t − s to t − 1 , respectively, s is the number of delay steps for autoregression, and [ y t − s − 1 , ⋯ , y t − d y ] are the observations from time t − d y to time t − s − 1 . The structure of NARNN is presented as Figure 1A.

The NARNN model is primarily concerned with historical series of the target variables as shown in Eq. 12. The information carried by exogenous inflow data is ignored in this modeling process, and then NARXNN model is proposed to make use of this information. The NARXNN can be formulated as follows: y ^ t = N ( y ^ t − 1 , y ^ t − 2 , ⋯ , y ^ t − s , y t − s − 1 , ⋯ , y t − d y , x t − 1 , ⋯ , x t − d x ) , (13) Where [ x t − 1 , ⋯ , x t − d x ] is a matrix consisting of exogenous input variables, the dimension of x depends on the quantity of exogenous input variables, and d x is the maximum delay index of exogenous input variables. The structure of NARXNN is different from the NARNN slightly, mainly with the addition of several extra inputs. The structure of NARXNN is presented as Figure 1B.

The two neural networks, NARNN and NARXNN, update weights in each layer by using the Bayesian regularization backpropagation algorithm (MacKay, 1992). Training samples are shown as following set, D = { ( x ¯ 1 , y ¯ 1 ) , ( x ¯ 2 , y ¯ 2 ) , ⋯ , ( x ¯ N , y ¯ N ) } , where x ¯ i ( i = 1,2 , ⋯ , N ) and y ¯ i ( i = 1,2 , ⋯ , N ) represent the input and output of the neural network, respectively. Given a neural network, called M , let g ( x ; w , M ) be the response of network M with respect to the input x , and w denotes the weight of network. The optimal parameters can be achieved by minimizing the quadratic cost function: E D = ∑ i = 1 N [ g ( x ¯ i ; w , M ) − y ¯ i ] 2 . (14)

The objective function is extended from F = E D to F = α E D + β E w to prevent the overfitting. The regularization term E w is denoted as: E w = ∑ n ω n 2 . (15)

Note that α and β are unknown parameters of the objective function F . When α is larger, the accuracy of the model to the training samples is enhanced, and similarly when β is larger, the generalization ability of the model is enhanced. The Bayesian regularization generally treats the network weights as random variables and the detailed methods for estimating the values of weights w , α and β can be found in (Dan Foresee and Hagan, 1997).

NARNN and NARXNN are variants of ARMA together with neural network, combing both the dynamic and recurrent properties. Both methods do not require strict stationarity of the target time series. On the other hand, it is should be noted that NARXNN needs more reasonable computational cost (Cadenas et al., 2016).

Feature selection is of vital importance to improve the performance of the model, especially whose predictions depend on a number of extrinsic inputs to some extent. Excellent choices of inputs not only help to provide accurate results, but also speed up calculations and reduce the number of sensor installations, all of which lead to operational cost savings. The SVM-RFE combined with clustering method proposed in this study ranks the features of the continuous processes based on backward elimination.

The first thing worth noting is that the support vector machine classifier achieves the distinction between two classes by searching for the optimal hyperplane in a high-dimensional space (Rakotomamonjy, 2003). For the binary classification problem with training data set { X , Y } , where X ∈ R n are the features and Y ∈ { − 1,1 } , there exists a hyperplane or decision function of the following form. f ( X ) = 〈 w , Φ ( X ) 〉 + b , (16) Where Φ ( X ) refers to the mapping relationship from features X to the high dimensional space. The parameters ( w , b ) are determined by minimizing the weights and the distance of each misclassified data to the hyperplane, before which the features X need to be normalized. The optimization problem can be written as: min w , ξ 1 2 ‖ w ‖ 2 + C ∑ k = 1 m ξ k 2 , (17) y k f ( x k ) ≥ 1 − ξ k s . t   y k ∈ Y ,   x k ∈ X   f o r ∀ k Where C is used as a penalty factor to weigh the importance of misclassification. SVM-RFE compares the impact of different remaining subsets on the classification by backward elimination of features, with the aim of preserving the subset of features that are most beneficial to the classification. The ranking of features is achieved through multiple iterations of elimination in the above work until the remaining feature subset is empty, and the criterion of elimination at each step for a given feature i can be expressed as: R c ( i ) = | ∇ ‖ w ( i ) ‖ 2 | , (18) Where ∇ denotes the weight difference between the previous subset and the one whose i t h feature is eliminated. The feature that minimizes R c will be removed after one round of loops, which means that the remaining feature subset has the least difference in classification performance from the feature subset containing the removed feature. It is worth mentioning that SVM-RFE was widely used in feature selection for binary classification problems, while it was rarely used in continuous classification problems.

In practical industrial processes, the variation in feature values often leads to continuous variation in the output. This begs a question that needs to be addressed. When each output point is treated directly as a separate label, the volume of the feature set corresponding to each label is so small that each feature set does not have the ability to characterize a specific label. Therefore, such continuous processes cannot be directly classified as a multi-label feature classification problem. In order to rank the features for this type of data, this study proposed SVM-RFE combined with clustering algorithm, as shown in Figure 2.

The successive outputs are first clustered to form new classes, and thus similar outputs can be grouped into homogeneous classes to enhance the differences between the new classes, e.g., { y 1 , y 2 , ⋯ , y n } → { L 1 , L 2 , ⋯ , L m } ,   n ≫ m . The features corresponding to the same type of label are also merged into the same group, e.g., { x 1 , x 2 , ⋯ , x n } → { X 1 , X 2 , ⋯ , X m } ,   n ≫ m . In this way, the original problem is successfully transformed into a multi-label feature classification problem. When SVM is used for multi-classification problems, it is usually transformed into a series of binary classification problems that are handled separately and then summarized for the final result. The new data set { ( X 1 , L 1 ) , ( X 2 , L 2 ) , ⋯ , ( X m , L m ) } is partitioned using the so-called the one-versus-all method, and subsequently the features are individually ranked using SVM-RFE. The final summary of the ranking results for each group is R m × v as shown in Eq. 19. R = ( R 1 ⋮ R m ) = ( r 11 ⋯ r 1 v ⋮ ⋱ ⋮ r m 1 ⋯ r m v ) , (19) Where R i ( i = 1 , ⋯ , m ) represents the ranking of the features that distinguish { X i , L i } from the rest of dataset using SVM-RFE and the order, r i 1 to r i v , is decreasing according to the importance of the v optional features. For the ranking matrix R , this research considers that the importance index of the same column is consistent, which unifies the weights assigned to the same column. Then the final ranking of all features, { r 1 , r 2 , ⋯ , r v } , can be obtained by counting the frequency of each feature in each column.

In the final step, the model is trained by sequentially increasing the number of exogenous inputs to the model, depending on the importance of the features. Then, the appropriate training set is obtained by comparing the model test results under the AIC criterion. SVM-RFE combined with clustering algorithm migrates the feature ranking method for binary classification problems to a new application scenario and solves the problem of feature ranking for continuous processes. This approach implements feature selection while keeping the original data of the features intact and visually explaining the input variables selection.

In this study, MSE and R are used as the performance evaluation metrics of the model, which can be calculated as follows: M S E = ∑ t = 1 n ( y ( t ) − y e t ( t ) ) 2 n , (20) R = ∑ t = 1 n ( y ( t ) − y ¯ ) ( y e t ( t ) − y ¯ e t ) ∑ t = 1 n ( y ( t ) − y ¯ ) 2 ∑ t = 1 n ( y e t ( t ) − y ¯ e t ) 2 , (21) Where y and y e t refer to the true values of the system output and the estimated values of the prediction model, respectively; y ¯ and y ¯ e t are the mean values of y and y e t . n is the total number of evaluation samples. The smaller the MSE, the smaller the error of the model. R is in the range of [ 0,1 ] , the closer it is to 1, the better the performance of model.