A new regularized and robust LSTM learning algorithm for industrial soft sensors

Abstract

Modern industrial processes often have nonlinearity, multivariate, time-delay, and measurement outliers, which make accurate data-driven modeling of key performance indicators difficult. To address these issues, a robust and regularized long short-term memory (LSTM) neural network for soft sensors in complex industrial processes was proposed. First, a conventional LSTM architecture was used as the basic model to deal with nonlinearity and time delay. Thereafter, a novel LSTM loss function that combines the excellent resistance to outliers of Huber M-Loss with the superior model reduction capability of ℓ₁ regularization was designed. Subsequently, a backpropagation through time training algorithm for the proposed LSTM was developed, including the chain derivative calculation and updating formulas. The adaptive moment estimation was applied to perform the gradient update, while the grid search and moving window cross-validation were used to find the optimal hyperparameters. Finally, nonlinear artificial datasets with time series and outliers, as well as an industrial dataset of a desulfurization process, were applied to investigate the performance of the proposed soft sensor. Simulation results show that the proposed algorithm outperforms other state-of-the-art soft sensors in terms of predictive accuracy and training time. The causal relationship of the data-driven soft sensor trained by the proposed algorithm is consistent with the field operation and chemical reactions of the desulfurization process.

Keywords

Soft sensor backpropagation through time huber M-Loss long short-term memory

1 Introduction

In modern industrial processes, the accurate and real-time monitoring of key performance indicators is necessary to ensure process safety, production efficiency, and product quality. However, some indicators, such as the slurry composition, fluid viscosity, and gas concentration, are difficult to measure directly using hardware sensors. Specific online analyzers are often expensive and time-consuming to measure. Soft sensors provide fast, low-cost, and convenient alternatives to estimate hard-to-measure variables with easy-to-measure variables, such as temperature, pressure, and flow rate [1].

Soft sensors are categorized into two classes: mechanism modeling [2] and data-driven modeling [3]. Mechanistic modeling approaches have been developed using theorems of physics or chemistry, such as mass balance, energy balance, and reaction principles. However, owing to the complexity and ambiguity of actual industrial processes, it is impractical to establish sufficiently accurate mechanistic models. Data-driven soft sensors use the historical data collected by process automation systems to establish mathematical models without much dependency on the mechanism knowledge of the processes and are therefore much easier to implement. Commonly used data-driven soft sensors consist of linear algorithms, such as principal component analysis [4], partial least squares [5], nonnegative garrote [6], and nonlinear algorithms, such as support vector machines [7] and neural networks (NN) [8].

Recently, NNs have increasingly been used in soft-sensor development owing to their powerful nonlinear mapping capabilities. Depending on whether there are loops in the structure, NNs are divided into two categories: feedforward and recurrent NNs. Feedforward NNs, such as the multilayer perceptron (MLP) [9], convolutional neural network [10] and stochastic configuration networks [11], pass the data in only one direction from the input layer to the output layer. Such a structure is ideally suitable for modeling processes that are not temporally dynamic and is inefficient for time-series processes [12]. Recurrent neural networks (RNNs) are more suitable for time-series processes because they provide bidirectional signal paths by introducing loops to the network architecture [13].

To overcome the gradient vanishing and gradient explosion problems of a standard RNN, long short-term memory (LSTM) and gated recurrent units (GRU) can be used to control the flow of information; gated units can be designed to forget previous useless information and store useful information in the memory cell for long periods and have therefore been widely studied in the field of soft sensors. For example, Yuan et al. [14] proposed a spatiotemporal attention-based LSTM network for key variable prediction in the hydrocracking of industrial crude oil. Liu et al. [15] designed an LSTM network with a spatiotemporal attention mechanism to forecast dissolved oxygen. In [16], an LSTM-based soft sensor was developed and applied to predict the oxygen content in the boiler flue gas. In [17], a novel integrated GRU stacked autoencoder for extracting features and learning dynamic relations was proposed and applied to the soft sensor for a debutanizer column. Ma et al. [18] designed a bidirectional GRU structure by combining future and historical information and applied it to the quality prediction of a hot rolling process. Although GRU has better computational efficiency than LSTM, it is less competent for tasks that require long-term memory and long-distance dependencies.

Vaswani et al. [19] proposed a Transformer model, which can capture time-series information by introducing multi-head self-attention mechanism and position encoding. Geng et al. [20] developed a dynamic soft-sensing model with a Transformer and gated convolutional neural network and applied it to practical industrial processes. However, the Transformer depended highly on the self-attention mechanism, whose computational complexity was quadratic to the length of the input sequence.

Owing to the increasing number of input variables, the structures of NN-based soft sensors have become increasingly complex, which brings many redundant neurons and weights and undermines the performance of soft sensors. Tibshirani [21] proposed an ℓ₁ regularization algorithm, called the least absolute shrinkage and selection operator (LASSO), to shrink the coefficients of multiple linear regression problems. Recently, the ℓ₁ regularization has been extended to multiple nonlinear regression models. Sun et al. introduced LASSO into MLP to perform variable selection and developed a soft sensor for practical crude oil processing [22]. Subsequently, Fan et al. [23] proposed a double LASSO for MLP, not only to shrink the input variables but also to simplify the hidden layer. Wu et al. [24] designed an input selection algorithm with LASSO for LSTM, which is referred to as LASSO-LSTM.

However, the aforementioned LSTM-based algorithms have two drawbacks. First, these algorithms use the mean square error (MSE) as the loss function for training and optimization owing to their excellent characteristics, such as convexity, smoothness, and low computational cost under the Gaussian assumption. In [25], the statistical results demonstrated that model estimation with the MSE loss function has weak resistance to outliers. Practical industrial processes are sometimes unstable, which inevitably introduces outliers in the measured data. To address this problem, the design of an effective robust loss function helps improve the robustness of soft sensors. There are various robust loss functions, such as Huber M-Loss [26], maximum correlation criterion [27], and least absolute deviation [28], implemented in linear or nonlinear models. Among them, Huber M-Loss, which combines squared loss with absolute loss, is quadratic for small errors but grows linearly for large errors. Specifically, the Huber M-Loss can adapt to the process with different error distributions by adjusting the tunable hyperparameter.

The second drawback is the computational cost of variable selection and structural optimization for LSTM. For example, the spatiotemporal attention mechanism needs to calculate the relevance between the input variables and output variables and then assign different weights to the input variable at each time step. LASSO-LSTM introduces the ℓ₁ penalty into a well-trained LSTM to create a quadratic optimization equation and then achieves variable selection by solving the equation. These algorithms are time consuming owing to the intrinsic complexity of the LSTM structure. Therefore, the development of fast training and easily implemented optimization methods for LSTM has become significant.

In this study, a new regularized and robust LSTM learning algorithm for industrial soft sensors is proposed. The main contributions of this study are as follows:

To improve the resistance to outliers and reduce the influence of redundant input variables, Huber M-Loss combined with ℓ₁ regularization was introduced into the loss function of LSTM.

A backpropagation through time (BPTT) training algorithm for the proposed LSTM was developed, including derivative calculation and updating formulas of network weights.

The superiority of the proposed soft sensor was demonstrated by artificial datasets and its application in an actual desulfurization process.

The feature importance analysis shows that the presented soft-sensing model has valuable interpretability that is consistent with chemical reactions and field operations.

The remainder of this paper is organized as follows. Section 2 provides a brief introduction to related works. In Section 3, the proposed algorithm is described in detail. In Section 4, the proposed algorithm is verified using artificial datasets and is compared with other approaches. In Section 5, the application of the proposed soft sensor to a practical desulfurization process and related discussions are presented. Finally, concluding remarks and future work are presented in Section 6.

2 Related work

In this study, all vectors are column vectors and are denoted by boldface lower-case letters. Matrices are represented in bold capital letters. For vector k (or matrix K ), k ^T ( K ^T) is the transpose vector (matrix). The time index is given as a subscript, for example, k _t is the vector at time t, and $k_{t}^{i}$ is the i^th member of vector k _t.

2.1 Long short-term memory neural network

LSTM can predict time series by retaining the state at each moment during the learning process. When a model is built, a sliding regression method must be introduced to handle the data. For the original time series $D = {{x_{1}^{T}, x_{2}^{T}, \dots, x_{L}^{T}}$ , {d₁, d₂, …, d_L }}, the sliding window m is used to construct a multi-step time series as ${(x_{1}^{T}, x_{2}^{T}, \dots, x_{m}^{T}), d_{p + m}}$ , ${(x_{2}^{T}, x_{3}^{T}, \dots, x_{m + 1}^{T}), d_{p + m + 1}} \dots {(x_{L - m - p + 1}^{T}, x_{L - m - p + 2}^{T}, \dots, x_{L - p}^{T}), d_{L}}$ , where L is the length of the original time series, x_t denotes the input data at the t^th moment, d_t denotes the response output at the t^th moment, p represents the prediction step, and m determines the memory capacity of LSTM. The structure of the LSTM used in this study is shown in Fig. 1.

Fig. 1

Unit structure of LSTM.

The three inputs are the current input $x_{t} \in ℝ^{s}$ , the previous memory unit $c_{t - 1} \in ℝ^{u}$ , and the previous external state $h_{t - 1} \in ℝ^{u}$ , where s and u are the numbers of input and hidden neurons, respectively. Each new time series is learned through three gate controllers f _t, i _t, and o _t, which determine the forgetting, preserving, and transmitting of information, respectively. The forward calculation procedure is expressed as follows:

$f_{t} = σ (W_{f} \cdot [h_{t - 1}, x_{t}] + b_{f})$ (1)

$i_{t} = σ (W_{i} \cdot [h_{t - 1}, x_{t}] + b_{i})$ (2)

$o_{t} = σ (W_{o} \cdot [h_{t - 1}, x_{t}] + b_{o})$ (3)

${\tilde{c}}_{t} = \tanh (W_{c} \cdot [h_{t - 1}, x_{t}] + b_{c})$ (4)

$c_{t} = f_{t} ⊙ c_{t - 1} + i_{t} ⊙ {\tilde{c}}_{t}$ (5)

$h_{t} = o_{t} ⊙ \tanh (c_{t})$ (6) Here, ⊙ represents the dot product, and σ (·) and tanh (·) are sigmoid and hyperbolic tangent activation functions, respectively. $W_{c}, W_{o}, W_{f}, W_{i} \in ℝ^{u \times (u + s)}$ and $b_{c}, b_{o}, b_{f}, b_{i} \in ℝ^{u}$ are the weight matrices and bias vectors, respectively. Considering the output h _t, the final estimate is obtained by

${\hat{d}}_{t} = w_{v}^{T} (o_{t} ⊙ \tanh (f_{t} ⊙ c_{t - 1} + i_{t} ⊙ {\tilde{c}}_{t})) + b_{v},$ (7) where $w_{v} \in ℝ^{u}$ and b_v are the final regression coefficient vector and bias, respectively.

2.2 ℓ₁ regularization

The LASSO for the multivariate linear regression model is formulated as:

$\hat{w} = \underset{w}{argmin} \sum_{\forall (x, d) \in [X, D]} {(d - w^{T} x)}^{2} + λ {| w |}_{1},$ (8) where w ^T = [w₁, w₂, …, w_q] is the model coefficient vector and ${| w |}_{1} = \sum_{i = 1}^{q} | w_{i} |$ . x = [x₁, x₂, …, x_q] ^T is the vector of the input variables and d is the output variable. $X \in ℝ^{N \times q}$ and $D \in ℝ^{N}$ are the input and output datasets, respectively, where N is the total number of samples. λ is the shrinkage parameter.

The first term in the formula is the residual loss that controls the approximation degree of the model, whereas the second term is the ℓ₁ penalty that regularizes the model coefficients. λ is a nonnegative tunable parameter that controls the shrinking strength. When λ equals zero, model (8) is equivalent to an ordinary least square estimate. As λ keeps increasing, the coefficients of some variables shrink to zero, and the variable selection or model reduction is achieved.

2.3 Robust estimation with Huber M-Loss

Huber M-Loss is robust and adaptive because it combines a quadratic loss with a linear loss, and its expression is given as

$ρ (d, \hat{d}) = {\begin{matrix} \frac{1}{2} {(\hat{d} - d)}^{2}, | \hat{d} - d | ⩽ δ \\ δ | \hat{d} - d | - \frac{1}{2} δ^{2}, | \hat{d} - d | > δ \end{matrix}$ (9)

Here, parameter δ controls the degree to which Huber loss is quadratic or linear. The quadratic loss is based on squared errors and punishes larger errors more, which makes it effective when the errors are normally distributed. However, if outliers exist in the data, the square part of the function results in a biased estimation of the model. By measuring the absolute differences, the linear loss is less sensitive to outliers in the data. Huber M-Loss incorporates quadratic loss and linear loss in one expression, which is robust to outliers and differentiable everywhere. The hyperparameter δ that regulates the degree of quadratic or linear, can be determined through a data-driven cross-validation approach.

3 Proposed methodology

In this study, a new soft-sensing algorithm based on robust and regularized LSTM is proposed. A detailed description of the proposed algorithm is provided below.

3.1 LSTM with Huber M-Loss and ℓ₁ regularization

Compared to the standard LSTM, there are two major revisions to the loss function of the proposed LSTM. First, to improve the resistance to outliers, the MSE loss of the conventional LSTM is replaced by the Huber M-Loss. By adjusting the value of δ in the Huber M-Loss, the LSTM can adapt to various processes with different error distributions. Second, the ℓ₁ penalty is added to the loss function to reduce the model complexity. As λ increases, the values of the weight matrices continuously shrink until some of them are forced to zero, implying that some redundant units are eliminated from the LSTM. Hence, it reduces the possibility of overfitting and improves the generalization performance.

Combining Equation (8) with Equation (9), the loss function is designed as

$\begin{matrix} l (d_{t}, {\hat{d}}_{t}) = ρ (d_{t}, {\hat{d}}_{t}) + λ {| W |}_{1} \\ {| W |}_{1} = {| W_{f} |}_{1} + {| W_{i} |}_{1} + {| W_{o} |}_{1} + | W_{c} |_{1}, \end{matrix}$ (10) where λ controls the shrinkage strength of the weight matrices, and δ determines the quadratic or linear degree of Huber M-Loss. The determination of the optimal combination of δ and λ is discussed in Section 3.2.

In this study, a BPTT training algorithm is developed for the proposed LSTM to minimize $l (d_{t}, {\hat{d}}_{t})$ . First, for the output gradient,

$\frac{l (d_{t}, {\hat{d}}_{t})}{\partial {\hat{d}}_{t}} = {\begin{matrix} {\hat{d}}_{t} - d_{t}, - δ ⩽ {\hat{d}}_{t} - d_{t} ⩽ δ \\ - δ, {\hat{d}}_{t} - d_{t} < - δ \\ δ, {\hat{d}}_{t} - d_{t} > δ \end{matrix}$ (11) Then, the external state of the LSTM at time t is h _t, and the error term at time t is defined as $δ_{t} \underset{=}{def} \frac{l (d_{t}, {\hat{d}}_{t})}{\partial h_{t}}$ ,

$\frac{l (d_{t}, {\hat{d}}_{t})}{\partial h_{t}} = \frac{l (d_{t}, {\hat{d}}_{t})}{\partial {\hat{d}}_{t}} \frac{\partial (w_{v}^{T} h_{t} + b_{v})}{\partial h_{t}} = w_{v} \frac{l (d_{t}, {\hat{d}}_{t})}{\partial {\hat{d}}_{t}}$ (12) The four components f _t, i _t, ${\tilde{c}}_{t}$ , o _t, and their corresponding error terms, are defined as follows:

$θ_{f, t} = W_{f} [h_{t - 1}, x_{t}] + b_{f, t} = W_{fh} h_{t - 1} + W_{fx} x_{t} + b_{f, t}$ (13) ∥ $θ_{i, t} = W_{i} [h_{t - 1}, x_{t}] + b_{i, t} = W_{ih} h_{t - 1} + W_{ix} x_{t} + b_{i, t}$ (14) ∥ $θ_{\tilde{c}, t} = W_{c} [h_{t - 1}, x_{t}] + b_{c, t} = W_{ch} h_{t - 1} + W_{cx} x_{t} + b_{c, t}$ (15) ∥ $θ_{o, t} = W_{o} [h_{t - 1}, x_{t}] + b_{o, t} = W_{oh} h_{t - 1} + W_{ox} x_{t} + b_{o, t},$ (16)

where $δ_{f, t} \overset{def}{=} \frac{l (d_{t}, {\hat{d}}_{t})}{\partial θ_{f, t}}$ , $δ_{i, t} \overset{def}{=} \frac{l (d_{t}, {\hat{d}}_{t})}{\partial θ_{i, t}}$ , $δ_{\tilde{c}, t} \overset{def}{=} \frac{l (d_{t}, {\hat{d}}_{t})}{\partial θ_{\tilde{c}, t}}$ , $δ_{o, t} \overset{def}{=} \frac{l (d_{t}, {\hat{d}}_{t})}{\partial θ_{o, t}}$ .

The gradient of the relevant parameters can be obtained by the chain derivative rule:

${\begin{matrix} \frac{l (d_{t}, {\hat{d}}_{t})}{\partial W_{oh}} = \sum_{j = 1}^{t} (δ_{o, j} h_{j - 1}^{T} + λ sign (W_{oh, j})) \\ \frac{l (d_{t}, {\hat{d}}_{t})}{\partial W_{ox}} = \sum_{j = 1}^{t} (δ_{o, j} x_{j}^{T} + λ sign (W_{ox, j})) \\ \frac{l (d_{t}, {\hat{d}}_{t})}{\partial b_{o}} = \sum_{j = 1}^{t} δ_{o, j} \end{matrix}$ (17)

${\begin{matrix} \frac{l (d_{t}, {\hat{d}}_{t})}{\partial W_{fh}} = \sum_{j = 1}^{t} (δ_{f, j} h_{j - 1}^{T} + λ sign (W_{fh, j})) \\ \frac{l (d_{t}, {\hat{d}}_{t})}{\partial W_{fx}} = \sum_{j = 1}^{t} (δ_{f, j} x_{j}^{T} + λ sign (W_{fx, j})) \\ \frac{l (d_{t}, {\hat{d}}_{t})}{\partial b_{f}} = \sum_{j = 1}^{t} δ_{f, j} \end{matrix}$ (18)

${\begin{matrix} \frac{l (d_{t}, {\hat{d}}_{t})}{\partial W_{ih}} = \sum_{j = 1}^{t} (δ_{i, j} h_{j - 1}^{T} + λ sign (W_{ih, j})) \\ \frac{l (d_{t}, {\hat{d}}_{t})}{\partial W_{ix}} = \sum_{j = 1}^{t} (δ_{i, j} x_{j}^{T} + λ sign (W_{ix, j})) \\ \frac{l (d_{t}, {\hat{d}}_{t})}{\partial b_{i}} = \sum_{j = 1}^{t} δ_{i, j} \end{matrix}$ (19)

${\begin{matrix} \frac{l (d_{t}, {\hat{d}}_{t})}{\partial W_{ch}} = \sum_{j = 1}^{t} (δ_{\tilde{c}, j} h_{j - 1}^{T} + λ sign (W_{ch, j})) \\ \frac{l (d_{t}, {\hat{d}}_{t})}{\partial W_{cx}} = \sum_{j = 1}^{t} (δ_{\tilde{c}, j} x_{j}^{T} + λ sign (W_{cx, j})) \\ \frac{l (d_{t}, {\hat{d}}_{t})}{\partial b_{c}} = \sum_{j = 1}^{t} δ_{\tilde{c}, j} \end{matrix}$ (20)

The detailed derivation process of the weight parameters of the proposed LSTM is shown in Appendix A. In practical applications, soft sensors are usually trained and executed in industrial computers, which may have different integrated development environments (IDEs). The BPTT formulas can be easily implemented in different IDEs for industrial applications. After completing the derivation of the BPTT algorithm for the proposed LSTM, the adaptive moment estimation (Adam) of [29] is implemented to perform the gradient update. Adam is an efficient optimization approach for the large-scale weight estimation of NNs and can provide fast, stable, and lower error convergence for the proposed model.

3.2 Hyperparameters optimization

Hyperparameters are essential for neural network learning because they have a significant impact on model performance. There are two categories of hyperparameters in the proposed algorithm. The first category comprises the regular training hyperparameters of the LSTM, such as the learning rate, number of hidden layer neurons, maximum epochs, and size m in the sliding window, which can be determined by experience and some trials. The second category comprises the advanced hyperparameters, including δ in Huber M-Loss and λ in ℓ₁ regularization, which are tuned by grid search and moving window cross-validation (MWCV) [30].

Grid search: In this study, the process is performed in three steps: (1) Set the bound domain of the hyperparameters. The bound domain of δ is set to [0 . 1, 1], and the bound domain of λ is set to [0, λ_ub], where λ_ub is sufficiently large to make the model underfit. (2) A list of possible hyperparameter combinations is generated from the defined domains. (3) Find the best hyperparameter combination through two loops, in which MWCV is used to evaluate the average loss function in the inner loop, and the outer loop ensures that every possible hyperparameter combination is traversed.

MWCV: In this study, we construct multiple windows that move synchronously along the sub-dataset with a step size of 1 to split a validation set Z_val = Z_k+1, and consider the union of the previous sub-datasets as the training dataset Z_train = [Z₁, Z₂ … Z_k]. Under the given combination, MWCV is repeated K - 1 times to obtain an average validation loss. After trying every possible combination of δ and λ, the combination with the smallest validation loss is considered the optimal [δ, λ].

Combining the above methods, the pseudocode for determining the optimal [δ, λ] is presented in Algorithm 1.

Algorithm 1 Pseudo-code of finding the optimal [δ, λ]
Input: Z ={ X _train, d _train } that is constructed from the startup
1. Generate possible hyperparameter combinations of δ and λ within predefined domains;
2. Separate Z into K sub-datasets [ Z ₁, Z ₂, …, Z _K];
3. for each possible combination: // Grid search
4. for k = 1 up to K - 1: // MWCV
5. Train a robust and regularized LSTM with Z _train = [ Z ₁, Z ₂, …, Z _k];
6. Evaluate the Huber M-Loss of the trained LSTM with Z _val = Z _k+1;
7. end for
8. Calculate average Huber M-Loss for the current combination;
9. end for
10. Output: the optimal combination [δ, λ] with the minimum Huber M-Loss.

3.3 Computational flow of proposed algorithm

The overall computational flow of the proposed algorithm is as follows: First, a multi-step time-series dataset is constructed and divided into a training set and a testing set. Second, the loss function of LSTM is reformulated by introducing Huber M-Loss and ℓ₁ regularization into the standard LSTM. Subsequently, the optimal [δ, λ] is determined using Algorithm 1. Finally, BPTT and Adam are used to train the proposed LSTM with the designed loss function. The pseudocode of the proposed ℓ₁-HM-LSTM is presented in Algorithm 2.

Algorithm 2 Pseudo-code of the ℓ₁-HM-LSTM
Input: the original time series dataset $D = {x_{1}^{T}, x_{2}^{T}, \dots, x_{L}^{T}}, {d_{1}, d_{2}, \dots, d_{L}}$
1. Construct a multi-step time series dataset for LSTM: D ′ = { X ₁, X ₂, …, X _n-p } , {d_p+m, d_p+m+1, … , d_L } where $X_{i} = {x_{i}^{T}, x_{i + 1}^{T}, \dots, x_{i + m - 1}^{T}}$ ;
2. Divide the dataset D ’ into { X _train, d _train} and { X _test, d _test};
3. Reformulate the LSTM loss function with Huber M-Loss and ℓ₁ regularization;
4. Obtain the optimal hyperparameter combination [δ^, λ^] via Algorithm 1;
5. Initialize all parameters W and b of robust LSTM model;
6. for i = 1 up to Epoch: // The training process of ℓ₁- HM-LSTM
7. Calculate the output of LSTM unit at each time step by solving Equations (1)–(7), and obtain the predicted value;
8. Compute the loss function l with Equation (10) under [δ^, λ^];
9. Use the BPTT algorithm of Equations (17)–(20) to compute the parametric gradients of the prediction error of the ℓ₁-HM-LSTM;
10. Update W and b with Adam optimizer to minimize the loss l;
11. end for
Output: soft sensor model of ℓ₁-HM-LSTM

4 Simulation results on artificial datasets

4.1 Experimental setup

In this study, the proposed algorithm was simulated on numerical cases and a practical industrial dataset and compared with LASSO-MLP [22], LSTM, Transformer, LASSO-LSTM [24], STA-LSTM [15] and ℓ₁-HM-GRU. To ensure fairness, all the LSTM-based soft-sensing algorithms had the same network structure and training methods, while the hyperparameters of LASSO-MLP were determined using cross-validation. The optimal hyperparameters of different algorithms involved in the experiments are shown in Appendix B.

The algorithms were simulated in the same environment, which had a Ryzen 7 5800 H (3.2 GHz) CPU with 16GB RAM, a Windows 11 operating system, and Python 3.6 programming software. The following statistical indices were used to evaluate the algorithms’ performance.

1) Root mean square error (RMSE):

$RMSE = \sqrt{\frac{1}{N_{test}} \sum_{i = 1}^{N_{test}} {(d_{i} - {\hat{d}}_{i})}^{2}}$ (21) 2) Mean absolute error (MAE):

$MAE = \frac{1}{N_{test}} \sum_{i = 1}^{N_{test}} | d_{i} - {\hat{d}}_{i} |$ (22) 3) Coefficient of determination (R²):

$R^{2} = 1 - \frac{\sum_{i = 1}^{N_{test}} {(d_{i} - {\hat{d}}_{i})}^{2}}{\sum_{i = 1}^{N_{test}} {(d_{i} - {\bar{d}}_{i})}^{2}},$ (23) where d_i and ${\hat{d}}_{i}$ are the real and predicted values, respectively, of the ith test sample. ${\bar{d}}_{i}$ represents the average of all the predicted output values, and N_test is the number of test samples.

4.2 Data generation

The artificial datasets were generated based on the Friedman dataset [31], which has five relevant input variables and 2000 samples. The output variable was derived as:

$\begin{matrix} y_{t} = & 10 \sin (π x_{t}^{1} x_{t}^{2}) + 20 {(x_{t}^{3} - 0 . 5)}^{2} + 10 \\ (x_{t - 1}^{4} + x_{t - 2}^{4} + x_{t - 3}^{4} + x_{t - 4}^{4} + x_{t - 5}^{4}) + 5 x_{t}^{5} + ɛ \end{matrix}$ (24)

Here, the independent variable x is randomly sampled from a uniform distribution (0,1) and ɛ denotes the Gaussian noise. After that, redundant variables and outliers were added to test the performance of the algorithm.

The detailed procedure of the dataset generation is as follows. First, the dataset of five relevant input variables (x¹, x², x³, x⁴, x⁵) was generated as X ₁, and the delay characteristic was added to the variable x⁴ with five moments [ $x_{t - 1}^{4}, x_{t - 2}^{4}, x_{t - 3}^{4}, x_{t - 4}^{4}, x_{t - 5}^{4}$ ]. Then, the dataset of the output variable Y , was obtained by Equation (24). Furthermore, the dataset of irrelevant variables [x⁶, x⁷, x⁸, …] was randomly generated as a redundant dataset X ₂, and let X ={ X ₁ ; X ₂ }. For adding outliers, it randomly selected samples from X and Y with a given contaminate ratio η and added the non-Gaussian differences Δx and Δy to the selected samples. The pseudocode of the data generation is presented in Algorithm 3.

Algorithm 3 Pseudo-code of data generation
Input: the number of irrelevant variables ℏ, the contamination ratio η;
1. Generate the dataset consisting of five relevant input variables as X ₁ ;
2. The delay characteristic is added to the variable x⁴ with five moments;
3. The output variable is calculated by Equation (24), construct the output dataset Y ;
4. Generate the dataset X ₂ consisting of ℏ irrelevant variables;
5. Let X ={ X ₁ ; X ₂ } ;
6. Outlier generation ( If needed ):
7. Randomly select samples from X and Y with a given contaminate ratio η;
8. Add the non-Gaussian differences Δx and Δy to the selected samples;
9. End if
Output: Dataset { X , Y }

4.3 Artificial dataset of 10 input variables with normal distribution

Five irrelevant input variables have been added to the test case, and the data are normally distributed. Table 1 shows the mean, standard deviation, and worst and best experimental results for each algorithm over 10 runs. The comparative results demonstrate that ℓ₁-HM-LSTM is better than the other algorithms for all accuracy indices. Meanwhile, the MLP-based model, as a typical feedforward neural network, is not suitable for modeling time series with different durations.

Table 1
Simulation results on the artificial dataset of 10 input variables with normal distribution

Algorithms R ² RMSE MAE Time

Mean Std Worst Best Mean Std Worst Best Mean Std Worst Best (s)

LASSO-MLP 0.234 0.016 0.197 0.260 6.828 0.074 6.993 6.714 5.531 0.039 5.616 5.469 2.3

LSTM 0.861 0.005 0.851 0.870 2.900 0.056 3.015 2.813 2.283 0.061 2.383 2.256 22.5

Transformer 0.873 0.006 0.864 0.884 2.778 0.069 2.871 2.656 2.176 0.069 2.281 2.099 33.7

LASSO-LSTM 0.877 0.009 0.854 0.887 2.724 0.099 2.975 2.622 2.174 0.098 2.439 2.078 231.4

STA-LSTM 0.887 0.002 0.884 0.891 2.622 0.024 2.654 2.574 2.038 0.019 2.066 2.004 203.1

ℓ₁-HM-GRU 0.972 0.002 0.969 0.975 1.282 0.048 1.361 1.218 1.015 0.035 1.074 0.971 25.3

ℓ₁-HM-LSTM 0.974 0.001 0.971 0.976 1.240 0.030 1.307 1.199 0.982 0.025 1.046 0.951 27.1

Algorithms	R ²	RMSE	MAE	Time
LASSO-MLP	0.234	0.016	0.197	0.260	6.828	0.074	6.993	6.714	5.531	0.039	5.616	5.469	2.3
LSTM	0.861	0.005	0.851	0.870	2.900	0.056	3.015	2.813	2.283	0.061	2.383	2.256	22.5
Transformer	0.873	0.006	0.864	0.884	2.778	0.069	2.871	2.656	2.176	0.069	2.281	2.099	33.7
LASSO-LSTM	0.877	0.009	0.854	0.887	2.724	0.099	2.975	2.622	2.174	0.098	2.439	2.078	231.4
STA-LSTM	0.887	0.002	0.884	0.891	2.622	0.024	2.654	2.574	2.038	0.019	2.066	2.004	203.1
ℓ₁-HM-GRU	0.972	0.002	0.969	0.975	1.282	0.048	1.361	1.218	1.015	0.035	1.074	0.971	25.3
ℓ₁-HM-LSTM	0.974	0.001	0.971	0.976	1.240	0.030	1.307	1.199	0.982	0.025	1.046	0.951	27.1

Moreover, the time cost of completing a single model training after all the hyperparameters are determined is compared. Notable is the fact that the time cost of the ℓ₁-HM-LSTM is close to that of LSTM and much smaller than other LSTM-based models. Because the proposed algorithm is a conventional LSTM with a modified loss function, the training cost does not increase significantly. However, both LASSO-LSTM and STA-LSTM are additional optimizations of the well-trained LSTM model, so their training costs increase significantly.

To demonstrate the effect of ℓ₁ regularization and the interpretability of the trained model, a mean impact value (MIV) analysis for all input variables was performed on the conventional LSTM and the ℓ₁-HM-LSTM. MIV is an indicator that measures the importance of the input variables to the output in complex black box models [32]. The computational steps are as follows: First, for every input variable of a well-trained model, the data of the corresponding column in the input dataset increases and decreases by 10% to obtain two modified datasets. Second, the two modified datasets are input into the model to obtain their corresponding outputs. Third, the absolute value of the difference in the outputs is calculated to quantify the effects of the input variable on the output, represented as Δ_i, i = 1, 2, …, s. After calculating the MIV of every input variable, an MIV vector Δ = (Δ₁, Δ₂, …, Δ_s) is obtained. The contribution ratio of variable i to the output can be measured as

$δ_{i} = \frac{Δ_{i}}{\sum_{j = 1}^{s} Δ_{j}} \times 100 %$ (25)

The contribution ratio assessment was performed with the testing data on the best models of the LSTM and the ℓ₁-HM-LSTM over 10 runs. The percentage diagrams of the contribution ratios of all variables are shown in Fig. 2.

Fig. 2

(a) Feature importance distribution of the LSTM model (b) Feature importance distribution of the ℓ₁-HM-LSTM model.

First, the relevant variables, that is, x₁ ∼ x₅, have larger shares of contribution ratio than irrelevant variables in both diagrams, which is consistent with Equation (24). Second, the contribution ratio of variable x₄ is significantly larger than that of the other variables because it has multiple time delays to the output, as shown in Equation (24). The results show that LSTM can effectively deal with short- and long-term dependencies between the input and output. Third, for the conventional LSTM, the total contribution ratio of irrelevant variables is so significant that it significantly deteriorates model performance. For the ℓ₁-HM-LSTM model, the contribution ratio of irrelevant variables is compressed to a negligible scale.

In addition, to verify the potential relationship between the optimal network structure and the feature importance analysis, the sum of the absolute values of the weights connecting each input variable in the optimal ℓ₁-HM-LSTM are calculated. Then, the percentage of each input variable to all variables is calculated according to their respective sums and shown in Fig. 3. The top five largest values are variables 4, 3, 1, 2, and 5, which is consistent with the feature importance analysis of Fig. 2. (b). The result further demonstrates the interpretability of the soft sensor model with our approach.

Fig. 3

Feature importance distribution of ℓ₁-HM-LSTM model with weights analysis.

To explore the model reduction ability of the ℓ₁ regularization, the shrinkage on the weights connecting a relevant variable and an irrelevant variable is presented. For the test case, the number of weights connecting one input variable is 120. The frequency distributions of the weights connecting the relevant variable x₄ and the irrelevant variable x₈ of the LSTM and the ℓ₁-HM-LSTM are shown in Figs. 4 and 5, respectively. First, the weights are normally distributed in the LSTM model, while most of the weights are shrunk to very small values close to zero in the ℓ₁-HM-LSTM model. The results demonstrate the powerful model reduction ability of the ℓ₁ regularization. Second, the weights connecting the irrelevant variable x₈ are shrunk much more completely than those of variable x₄, illustrating the shrinking precision of the ℓ₁ regularization.

Fig. 4

(a) Frequency distribution of the weights connecting x₄ before sparsity (b) Frequency distribution of the weights connecting x₄ after sparsity.

Fig. 5

(a) Frequency distribution of the weights connecting x₈ before sparsity (b) Frequency distribution of the weights connecting x₈ after sparsity.

4.4 Artificial dataset of 10 input variables with outliers

The artificial dataset consists of 10 input variables and 10% of the samples are augmented with outliers to demonstrate the robustness of the proposed ℓ₁-HM-LSTM. The simulation results for the modified dataset are presented in Table 2.

Table 2
Simulation results on the artificial dataset of 10 input variables with 10% outliers

Algorithms R ² RMSE MAE Time

Mean Std Worst Best Mean Std Worst Best Mean Std Worst Best (s)

LASSO-MLP 0.161 0.009 0.148 0.178 6.800 0.038 6.851 6.731 5.458 0.053 5.523 5.382 2.18

LSTM 0.648 0.045 0.563 0.712 4.619 0.294 5.154 4.182 3.633 0.244 4.099 3.269 22.7

Transformer 0.786 0.010 0.766 0.800 3.608 0.084 3.768 3.488 2.900 0.075 3.009 2.794 32.6

LASSO-LSTM 0.818 0.012 0.800 0.838 3.163 0.111 3.317 2.986 2.515 0.097 2.650 2.390 233.6

STA-LSTM 0.866 0.006 0.854 0.876 2.728 0.062 2.847 2.623 2.151 0.065 2.287 2.068 211.4

ℓ₁-HM-GRU 0.965 0.003 0.959 0.970 1.450 0.081 1.564 1.336 1.129 0.067 1.230 1.029 25.1

ℓ₁-HM-LSTM 0.969 0.001 0.968 0.973 1.352 0.044 1.417 1.271 1.064 0.032 1.114 1.000 26.9

Algorithms	R ²	RMSE	MAE	Time
LASSO-MLP	0.161	0.009	0.148	0.178	6.800	0.038	6.851	6.731	5.458	0.053	5.523	5.382	2.18
LSTM	0.648	0.045	0.563	0.712	4.619	0.294	5.154	4.182	3.633	0.244	4.099	3.269	22.7
Transformer	0.786	0.010	0.766	0.800	3.608	0.084	3.768	3.488	2.900	0.075	3.009	2.794	32.6
LASSO-LSTM	0.818	0.012	0.800	0.838	3.163	0.111	3.317	2.986	2.515	0.097	2.650	2.390	233.6
STA-LSTM	0.866	0.006	0.854	0.876	2.728	0.062	2.847	2.623	2.151	0.065	2.287	2.068	211.4
ℓ₁-HM-GRU	0.965	0.003	0.959	0.970	1.450	0.081	1.564	1.336	1.129	0.067	1.230	1.029	25.1
ℓ₁-HM-LSTM	0.969	0.001	0.968	0.973	1.352	0.044	1.417	1.271	1.064	0.032	1.114	1.000	26.9

Compared to the results on the artificial dataset with normal distribution, the performance of all algorithms deteriorates. The ℓ₁-HM-LSTM shows better performance than the other algorithms, with small performance degradation.

4.5 Artificial datasets of 50 input variables with outliers

To further verify the robustness of the proposed algorithm when the number of outliers and redundancy increase, larger artificial datasets of 50 input variables with different degrees of outliers were constructed. The contamination ratio η ranged from 0.1 to 0.6. The statistical results with η = 0.1 are shown in Table 3. It can be seen that LSTM suffers from overfitting owing to the existence of redundant input variables and outliers, while the proposed ℓ₁-HM-LSTM shows significant advantages over the other algorithms.

Table 3
Simulation results on the artificial dataset of 50 input variables with 10% outliers

Algorithms R ² RMSE MAE Time

Mean Std Worst Best Mean Std Worst Best Mean Std Worst Best (s)

LASSO-MLP 0.113 0.018 0.083 0.145 6.747 0.071 6.861 6.626 5.444 0.061 5.550 5.332 39.7

LSTM 0.181 0.031 0.145 0.249 6.558 0.126 6.704 6.280 5.211 0.113 5.415 5.046 27.3

Transformer 0.529 0.021 0.493 0.559 4.971 0.110 5.158 4.819 3.976 0.093 4.121 3.778 101.5

LASSO-LSTM 0.599 0.059 0.447 0.679 4.525 0.327 5.328 4.058 3.635 0.293 4.333 3.243 263.3

STA-LSTM 0.835 0.003 0.827 0.840 2.916 0.032 2.983 2.873 2.295 0.024 2.336 2.272 215.9

ℓ₁-HM-GRU 0.878 0.007 0.869 0.893 2.527 0.081 2.619 2.367 2.063 0.064 2.136 1.939 28.7

ℓ₁-HM-LSTM 0.911 0.008 0.902 0.931 2.152 0.101 2.261 1.895 1.710 0.089 1.798 1.479 31.6

Algorithms	R ²	RMSE	MAE	Time
LASSO-MLP	0.113	0.018	0.083	0.145	6.747	0.071	6.861	6.626	5.444	0.061	5.550	5.332	39.7
LSTM	0.181	0.031	0.145	0.249	6.558	0.126	6.704	6.280	5.211	0.113	5.415	5.046	27.3
Transformer	0.529	0.021	0.493	0.559	4.971	0.110	5.158	4.819	3.976	0.093	4.121	3.778	101.5
LASSO-LSTM	0.599	0.059	0.447	0.679	4.525	0.327	5.328	4.058	3.635	0.293	4.333	3.243	263.3
STA-LSTM	0.835	0.003	0.827	0.840	2.916	0.032	2.983	2.873	2.295	0.024	2.336	2.272	215.9
ℓ₁-HM-GRU	0.878	0.007	0.869	0.893	2.527	0.081	2.619	2.367	2.063	0.064	2.136	1.939	28.7
ℓ₁-HM-LSTM	0.911	0.008	0.902	0.931	2.152	0.101	2.261	1.895	1.710	0.089	1.798	1.479	31.6

To demonstrate how the hyperparameters δ and λ affect the performance of the proposed algorithm, a dataset of 50 input variables and η = 0 . 3 was taken. After the MWCV process,λ = 0 . 0015 and δ = 0 . 4 were determined as the optimal hyperparameter combination of the ℓ₁-HM-LSTM. Then, the model performance was compared by regulating one hyperparameter while holding the other constant to illustrate the influence of a single hyperparameter.

First, let λ = 0 . 0015; the RMSE with different δ is given in Fig. 6. A small δ means that the Huber loss function favors MAE and therefore makes the model estimation biased, while a large δ means that the function prefers MSE and makes the model less sensitive to outliers. Second, let δ = 0 . 4; the RMSE with different λ is presented in Fig. 7. Obviously, the smallest RMSE appears at λ = 0 . 0015, where the redundant nodes are well shrunk, and an optimal model is obtained. If λ is very small, it means there is little shrinkage on the redundant nodes, which makes the model performance inadequate. As the λ continues to increase beyond the optimal value, more and more valuable nodes are removed, which makes the model performance deteriorate successively.

Fig. 6

The validation RMSE of ℓ₁-HM-LSTM with different δ when λ = 0 . 0015.

Fig. 7

The validation RMSE of ℓ₁-HM-LSTM with different λ when δ = 0 . 4.

The comparative results with different values of η are shown in Fig. 8, which demonstrates that ℓ₁-HM-LSTM has the best performance in each evaluation index. Especially, when η > 0 . 4, the R² of the other algorithms other than ℓ₁-HM-GRU is below 0.7, which means that the model performance of these algorithms began to collapse. Besides, the ℓ₁-HM-GRU is obviously worse than ℓ₁-HM-LSTM in prediction accuracy. Since the GRU has only one update gates to control information updating, it is less effective than LSTM in handling complex sequence modeling tasks with long-term dependencies. Although both ℓ₁ regularization and STA technology can reduce the model complexity and improve the model performance, the LSTM with MSE loss has insufficient resistance to outliers. The proposed Huber M-Loss with ℓ₁ regularization not only simplifies the LSTM model, but has sufficient resistance to outliers, and therefore produces more robust models.

Fig. 8

Performance comparison with different η.

5 Application to the SO₂ estimation of a desulfurization process

5.1 Process description

In this section, the proposed algorithm is applied to predict the SO₂ content in the gas emitted from the desulfurization process of an actual thermal power plant. In the plant, a double absorption tower structure with limestone-gypsum wet flue gas desulfurization (FGD) technology is adopted. A schematic of the process is shown in Fig. 9.

Fig. 9

Process flow chart of double-tower desulfurization system.

The process flow is briefly described as follows: First, the lump limestone is ground and crushed into powder using a wet ball mill and mixed with water to generate limestone slurry, which is then further treated and entered the oxidation reaction pool in the absorption tower to form a circulating slurry. Second, the circulating slurry is sent into the spray layer by a circulating water pump, which is atomized in the spray layer and fully in contact with the flue gas to produce gypsum slurry. Finally, the gypsum slurry is dewatered by the gypsum treatment system to form gypsum by-products. The moisture of the purified flue gas is removed by the absorption tower demister, and the dust is removed by the electrostatic precipitator before it is discharged into the atmosphere. The principle of the chemical reaction is as follows.

$\begin{matrix} CaC O_{3} + S O_{2} + \frac{1}{2} O_{2} + 2 H_{2} O \\ ⇄ CaS O_{4} \cdot 2 H_{2} O + C O_{2} . \end{matrix}$ (26)

In this system, the SO₂ concentration in the flue gas emitted from the secondary absorption tower is a key performance indicator related to environmental standards and energy consumption. However, this indicator is unstable, and an online analyzer requires regular maintenance. It is necessary to design soft sensors to provide calibration information and improve the system’s reliability. In addition, the indicator can be used as the output to design a model predictive control system to improve absorption efficiency and reduce energy consumption.

However, there are two challenges when performing data-driven modeling of the process. First, the double-tower structure increases the duration and complexity of the process, which leads to a considerable time delay between the input and output, as well as model redundancy. Second, owing to process disturbances, instrument deviations, and personnel operation errors, there are some outliers in the measured data. Therefore, it is of great significance to develop a soft sensor with long-term memory, resistance to outliers, and reduced model redundancy.

5.2 Predicting the SO₂ concentration from the emitting flue gas

The dataset collected from the FGD process consists of 53 variables, of which SO₂ is the response variable. After data analysis and preprocessing, a candidate input variable pool consisting of 35 variables was determined, as shown in Table 4. A total of 2676 samples were collected at one-minute sampling intervals, the first 80% of which were taken for training, and the others were taken for testing.

Table 4
Candidate input variables of desulfurization system

No. Name No. Name

1 Generator power 19 #1 tower exit temp. 3

2 Furnace raw FG SO₂ content 20 Estimated value of DE.

3 #1 tower exit FG SO₂ content 21 Chimney inlet FG flow rate

4 #1 tower gypsum slurry pH 22 Furnace raw FG flow rate

5 Limestone slurry flow rate 23 Chimney inlet FG press.

6 #2 tower gypsum slurry pH 24 Furnace raw FG soot content

7 #2 tower feed flow rate 25 Furnace raw FG NO_x content

8 #1 tower #1 CSP current 26 Furnace raw FG O₂ content

9 #1 tower #2 CSP current 27 Estimated median of DE

10 #1 tower #3 CSP current 28 #2 tower exit temp. 1

11 #1 tower #4 CSP current 29 #2 tower exit temp. 2

12 #2 tower #1 CSP current 30 #2 tower exit press.

13 #2 tower #2 CSP current 31 #1 tower exit FG O₂ content

14 #1 tower inlet FG temp. 32 #1 tower liquid level

15 Total air volume 33 #2 tower liquid level

16 Total coal 34 Slurry discharge pump current

17 #1 tower exit FG temp. 1 35 Slurry cyclone inlet press.

18 #1 tower exit FG temp. 2

No.	Name	No.	Name
1	Generator power	19	#1 tower exit temp. 3
2	Furnace raw FG SO₂ content	20	Estimated value of DE.
3	#1 tower exit FG SO₂ content	21	Chimney inlet FG flow rate
4	#1 tower gypsum slurry pH	22	Furnace raw FG flow rate
5	Limestone slurry flow rate	23	Chimney inlet FG press.
6	#2 tower gypsum slurry pH	24	Furnace raw FG soot content
7	#2 tower feed flow rate	25	Furnace raw FG NO_x content
8	#1 tower #1 CSP current	26	Furnace raw FG O₂ content
9	#1 tower #2 CSP current	27	Estimated median of DE
10	#1 tower #3 CSP current	28	#2 tower exit temp. 1
11	#1 tower #4 CSP current	29	#2 tower exit temp. 2
12	#2 tower #1 CSP current	30	#2 tower exit press.
13	#2 tower #2 CSP current	31	#1 tower exit FG O₂ content
14	#1 tower inlet FG temp.	32	#1 tower liquid level
15	Total air volume	33	#2 tower liquid level
16	Total coal	34	Slurry discharge pump current
17	#1 tower exit FG temp. 1	35	Slurry cyclone inlet press.
18	#1 tower exit FG temp. 2

*FG = flue-gas; CSP = circulating slurry pump; DE = desulfurization efficiency.

In this study, the dataset was taken without any modification, and the anomaly detection toolkit (ADTK) [33] was used to detect the presence of outliers from the dataset. Figure 10 shows the detected outliers of the studied dataset with ADTK criteria, in which 71 samples were identified as outliers.

Fig. 10

Distribution of measured SO₂ concentration.

Table 5 shows the statistical results of the SO₂ prediction with different soft sensors. First, it can be observed that the conventional LSTM model performs poorly owing to excessive irrelevant variables and outliers. Second, the LASSO-MLP shows unsatisfactory performance owing to the time delays with unknown durations caused by the design of the double-tower structure and the complex reaction mechanism. The feedforward structure of MLP is inefficient for the modeling of time-series processes. Third, although the Transformer model using the attention mechanism can capture global key features, the encoder is easily disturbed by outliers and has poor robustness. Fourth, although LASSO-LSTM and STA-LSTM can capture the time-series characteristics of the process and simplify the model, their performance is less favorable. The proposed ℓ₁-HM-LSTM is superior to other algorithms in all criteria, indicating that model redundancy and outliers in the measured data are handled well. Our approach has better generalization performance and higher robustness in industrial processes with outliers. The best predictive results with ℓ₁-HM-LSTM and different algorithms’ error distributions are shown in Figs. 11 and 12, respectively. The fitting curve demonstrates that the proposed algorithm can accurately follow the dynamics of the target variable. The box-plot error distribution of our approach is the densest, which illustrates the superiority of our algorithm in actual industrial applications.

Fig. 11

Fitting curve of SO₂ concentration with ℓ₁-HM-LSTM.

Fig. 12

Error box plots of different algorithms.

Table 5

Simulation results on the SO₂ prediction

Algorithms	R ²				RMSE				MAE				Time
	Mean	Std	Worst	Best	Mean	Std	Worst	Best	Mean	Std	Worst	Best	(s)
LASSO-MLP	0.496	0.048	0.436	0.579	3.204	0.155	3.395	2.933	2.610	0.119	2.760	2.379	17.3
LSTM	0.444	0.039	0.402	0.514	3.370	0.121	3.498	3.152	2.663	0.131	2.798	2.386	35.1
Transformer	0.725	0.010	0.706	0.740	2.371	0.046	2.450	2.302	1.757	0.078	1.900	1.656	74.4
LASSO-LSTM	0.732	0.012	0.717	0.752	2.337	0.052	2.404	2.250	1.853	0.050	1.947	1.764	345.6
STA-LSTM	0.777	0.007	0.756	0.782	2.133	0.035	2.234	2.109	1.743	0.037	1.853	1.719	218.5
ℓ₁-HM-GRU	0.822	0.010	0.807	0.835	1.905	0.055	1.969	1.834	1.551	0.053	1.631	1.480	40.2
ℓ₁-HM-LSTM	0.850	0.009	0.830	0.864	1.747	0.054	1.860	1.667	1.340	0.054	1.460	1.279	42.7

5.3 Feature importance analysis

To demonstrate the interpretability of the soft-sensor model trained with our approach, a feature importance analysis was performed on the conventional LSTM and the ℓ₁-HM-LSTM, with the percentage diagrams of the contribution ratio shown in Fig. 13. In the percentage diagram of the conventional LSTM, the contribution ratio of each variable is relatively even. The relevant input variables do not contribute significantly more than the irrelevant variables in the LSTM model. In contrast, the distribution of the contribution ratio of the ℓ₁-HM-LSTM is much more explicit.

Fig. 13

Feature importance distribution of the LSTM model and the ℓ₁-HM-LSTM model.

In the ℓ₁-HM-LSTM model, the first variable that prevails is variable 3, and the outlet flue gas SO₂ content of absorption tower #1. It is highly correlated with the output. The second highest contributing variable is variable 12, the current of #1 circulating slurry pump of absorption tower #2. In the field operation, this variable is used to control the amount of slurry sprayed in absorption tower #2, which is important for the SO₂ absorption. Moreover, variables 2, 5, and 8 are the other three most influential variables. Variable 2 is the inlet flue gas SO₂ content of absorption tower #1, which has an inevitable but smaller impact on the output than variable 3. Variable 5 is the flow rate of the limestone slurry in the absorption tower, which appears to influence the SO₂ absorption of the system, as shown in Equation (26). Variable 8, that is, the #1 circulating slurry pump of absorption tower #1, is used to control the amount of slurry sprayed in absorption tower #1.

There are four circulating pumps in absorption tower #1 and three pumps in absorption tower #2. Under normal field conditions, the “3 + 2” operation mode is adopted, in which three pumps in tower #1 and two in tower #2 are kept running, and the others are idle. Sometimes, the “2 + 1” operation mode is taken to save energy when the load on the power plant is low at night. The #1 pumps of both towers run as manipulated variables to control the SO₂ absorption. As it is closer to the process outlet, pump #1 of tower #2 contributes more to the output variable than pump #1 of tower #1, as shown in Fig. 8. In summary, the feature importance distribution presented by the ℓ₁-HM-LSTM is consistent with the chemical reaction and field operation of the plant.

5.4 Ablation study

To further illustrate the contribution of each component in the ℓ₁-HM-LSTM model, ablation experiments were conducted on the artificial dataset of 50 input variables with outliers and η = 0.3 as well as the FGD process dataset. The experimental results are shown in Fig. 14. It can be observed that ℓ₁ regularization and Huber M-Loss are effective for the experimental cases, and removing each component leads to a decline in the model performance. By comparison, removing the ℓ₁ regularization has a more significant impact on performance than removing the Huber M-Loss does, indicating that ℓ₁ regularization with model parsimony contributes more to prediction accuracy when there are excessive redundant variables and outliers.

Fig. 14

Ablation study results in terms of (a) the R², (b) the RMSE and (c) the MAE on the artificial dataset and the FGD dataset.

6 Conclusion

Considering that the soft sensing of practical industrial processes often experiences multivariate, time series, and outliers, a robust LSTM with Huber M-Loss and ℓ₁ regularization was proposed. The proposed algorithm takes advantage of the adaptability of the Huber M-Loss to outliers and applies ℓ₁ regularization to achieve input variable selection and model reduction. The weight-updating formulas for the designed LSTM were derived using the BPTT algorithm. The hyperparameters of Huber M-Loss and ℓ₁ regularization were determined using the MWCV method. Numerical datasets with different dimensions and outliers were used to test the performance. The comparative results with other approaches demonstrate the excellent robustness and rapidity of the proposed algorithm. Furthermore, the proposed algorithm was implemented for the SO₂ estimation of an actual desulfurization process. Our approach presents sufficiently accurate models that can provide a solid foundation for the optimization and control system improvement of the process.

Regarding future research, our approach can be improved in a few ways. First, the proposed algorithm is limited to regression problems, and the algorithm can be extended to classification problems. Moreover, the proposed algorithm can only train offline models, which is very effective for the process with stable operating conditions. However, the model performance degrades when concept drift occurs or operating conditions change. Therefore, another direction for future research is to develop online learning or incremental learning techniques for such algorithms.

Footnotes

Acknowledgments

The authors would like to thank the editors and reviewers for the time and effort they spent reviewing this paper as well as for their detailed and constructive comments for the paper improvement in terms of presentation and quality.

Funding

This work is funded in part by the Open Foundation of State Key Laboratory of Process Automation in Mining & Metallurgy (BGRIMM-KZSKL-2021-07), the Shandong Provincial Natural Science Foundation of China (ZR2021MF022).

Appendix A

In this section, the detailed derivation of the gradient update equations for the proposed ℓ₁-HM-LSTM is provided as follows. At time t - 1, the error term of BPTT is defined as,

(A.1) $A . 1 δ_{t - 1}^{T} = \frac{l (d_{t}, {\hat{d}}_{t})}{\partial h_{t - 1}} = \frac{l (d_{t}, {\hat{d}}_{t})}{\partial h_{t}} \frac{\partial h_{t}}{\partial h_{t - 1}} = δ_{t}^{T} \frac{\partial h_{t}}{\partial h_{t - 1}}$ According to Equations (1)–(4), f _t, i _t, ${\tilde{c}}_{t}$ , and o _t are functions of h _t-1. Following the chain derivation rule,

(A.2) $\begin{matrix} δ_{t}^{T} \frac{\partial h_{t}}{\partial h_{t - 1}} \\ = δ_{t}^{T} \frac{\partial h_{t}}{\partial o_{t}} \frac{\partial o_{t}}{\partial θ_{o, t}} \frac{\partial θ_{o, t}}{\partial h_{t - 1}} + δ_{t}^{T} \frac{\partial h_{t}}{\partial c_{t}} \frac{\partial c_{t}}{\partial f_{t}} \frac{\partial f_{t}}{\partial θ_{f, t}} \frac{\partial θ_{f, t}}{\partial h_{t - 1}} + \\ δ_{t}^{T} \frac{\partial h_{t}}{\partial c_{t}} \frac{\partial c_{t}}{\partial i_{t}} \frac{\partial i_{t}}{\partial θ_{i, t}} \frac{\partial θ_{i, t}}{\partial h_{t - 1}} + δ_{t}^{T} \frac{\partial h_{t}}{\partial c_{t}} \frac{\partial c_{t}}{\partial {\tilde{c}}_{t}} \frac{\partial {\tilde{c}}_{t}}{\partial θ_{\tilde{c}, t}} \frac{\partial θ_{\tilde{c}, t}}{\partial h_{t - 1}} \\ = δ_{o, t}^{T} \frac{\partial θ_{o, t}}{\partial h_{t - 1}} + δ_{f, t}^{T} \frac{\partial θ_{f, t}}{\partial h_{t - 1}} + δ_{i, t}^{T} \frac{\partial θ_{i, t}}{\partial h_{t - 1}} + δ_{\tilde{c}, t}^{T} \frac{\partial θ_{\tilde{c}, t}}{\partial h_{t - 1}} \end{matrix}$ Next, each partial derivative of Equation (A.2) is calculated. Considering the output gate o _t as an example, according to Equation (6), the partial derivative $\frac{\partial h_{t}}{\partial o_{t}}$ can be obtained as

(A.3) $\begin{matrix} \frac{\partial h_{t}}{\partial o_{t}} ≜ & (\begin{matrix} \frac{\partial h_{t}^{1}}{\partial o_{t}^{1}} & \dots & \frac{\partial h_{t}^{1}}{\partial o_{t}^{u}} \\ ⋮ & ⋱ & ⋮ \\ \frac{\partial h_{t}^{u}}{\partial o_{t}^{1}} & \dots & \frac{\partial h_{t}^{u}}{\partial o_{t}^{u}} \end{matrix}) \\ = (\begin{matrix} tanh (c_{t}^{1}) & \dots & 0 \\ ⋮ & ⋱ & ⋮ \\ 0 & \dots & tanh (c_{t}^{u}) \end{matrix}) = Π^{(tanh (c_{t}))}, \end{matrix}$ where Πs^(tanh(c_t)) = diag (tanh(c_t)). Similarly, because o _t = σ ( θ_o,t), we obtain $\frac{\partial o_{t}}{\partial θ_{o, t}} = Π^{(o_{t} ⊙ (1 - o_{t}))}$ .

From Equation (16), we obtain

(A.4) $\begin{matrix} \frac{\partial θ_{o, t}}{\partial h_{t - 1}} = & (\begin{matrix} \frac{\partial θ_{o, t}^{1}}{\partial h_{t - 1}^{1}} & \dots & \frac{\partial θ_{o, t}^{1}}{\partial h_{t - 1}^{u}} \\ ⋮ & ⋱ & ⋮ \\ \frac{\partial θ_{o, t}^{u}}{\partial h_{t - 1}^{1}} & \dots & \frac{\partial θ_{o, t}^{u}}{\partial h_{t - 1}^{u}} \end{matrix}) \\ = (\begin{matrix} w_{11} & \dots & w_{1 u} \\ ⋮ & ⋱ & ⋮ \\ w_{u 1} & \dots & w_{uu} \end{matrix}) = W_{oh} \end{matrix}$ By computing all the partial derivatives in Equation (A.2), we obtain

(A.5) $\begin{matrix} δ_{t - 1}^{T} \\ = δ_{o, t}^{T} \frac{\partial θ_{o, t}}{\partial h_{t - 1}} + δ_{f, t}^{T} \frac{\partial θ_{f, t}}{\partial h_{t - 1}} + δ_{i, t}^{T} \frac{\partial θ_{i, t}}{\partial h_{t - 1}} + δ_{\tilde{c}, t}^{T} \frac{\partial θ_{\tilde{c}, t}}{\partial h_{t - 1}} \\ = δ_{o, t}^{T} W_{oh} + δ_{f, t}^{T} W_{fh} + δ_{i, t}^{T} W_{ih} + δ_{\tilde{c}, t}^{T} W_{ch} \end{matrix}$ According to the definitions, δ_o,t, δ_f,t, δ_i,t, and $δ_{\tilde{c}, t}$ are expressed as

(A.6) $δ_{o, t}^{T} = δ_{t}^{T} ⊙ \tanh {(c_{t})}^{T} ⊙ o_{t}^{T} ⊙ {(1 - o_{t})}^{T}$

(A.7) $\begin{matrix} δ_{f, t}^{T} & = δ_{t}^{T} ⊙ {[o_{t} ⊙ (1 - \tanh {(c_{t})}^{2})]}^{T} \\ ⊙ c_{t - 1}^{T} ⊙ [f_{t} ⊙ (1 - f_{t})]^{T} \end{matrix}$

(A.8) $\begin{matrix} δ_{i, t}^{T} & = δ_{t}^{T} ⊙ {[o_{t} ⊙ (1 - \tanh {(c_{t})}^{2})]}^{T} \\ ⊙ {\tilde{c}}_{t}^{T} ⊙ [i_{t} ⊙ (1 - i_{t})]^{T} \end{matrix}$

(A.9) $\begin{matrix} δ_{\tilde{c}, t}^{T} & = δ_{t}^{T} ⊙ {[o_{t} ⊙ (1 - \tanh {(c_{t})}^{2})]}^{T} \\ ⊙ i_{t}^{T} ⊙ {(1 - {\tilde{c}}_{t}^{2})}^{T} \end{matrix}$ Equations (A.6)–(A.9) are the backpropagation of errors over time at one moment. Substituting them into Equation (A.5), the error term δ_t-1 is obtained as

(A.10) $\begin{matrix} δ_{t - 1}^{T} = δ_{t}^{T} ⊙ tanh {(c_{t})}^{T} ⊙ o_{t}^{T} ⊙ {(1 - o_{t})}^{T} W_{o h} + δ_{t}^{T} ⊙ \\ {[o_{t} ⊙ (1 - \tanh {(c_{t})}^{2})]}^{T} \\ ⊙ c_{t - 1}^{T} ⊙ [f_{t} ⊙ (1 - f_{t})]^{T} W_{f h} + \\ δ_{t}^{T} ⊙ {[o_{t} ⊙ (1 - \tanh {(c_{t})}^{2})]}^{T} \\ ⊙ {\tilde{c}}_{t}^{T} ⊙ [i_{t} ⊙ (1 - i_{t})]^{T} W_{i h} + \\ δ_{t}^{T} ⊙ {[o_{t} ⊙ (1 - \tanh {(c_{t})}^{2})]}^{T} ⊙ i_{t}^{T} \\ ⊙ {(1 - {\tilde{c}}_{t}^{2})}^{T} W_{c h} \end{matrix}$

Note that we only provide the variations applied to the output gate of the LSTM network, whereas the variations for the remaining gates can be derived in a similar manner. The weight gradient for W _oh is the sum of the gradients at each moment.

At time t, the gradient of W _oh,t can be expressed as

(A.11) $\begin{matrix} \frac{l (d_{t}, {\hat{d}}_{t})}{\partial W_{oh, t}} & = \frac{l (d_{t}, {\hat{d}}_{t})}{\partial θ_{o, t}} \frac{\partial θ_{o, t}}{\partial W_{oh, t}} + λ \frac{\partial {| W |}_{1}}{\partial W_{oh, t}} \\ = δ_{o, t} h_{t - 1}^{T} + λ sign (W_{oh, t}), \end{matrix}$ where

(A.12) $sign (x) = {\begin{matrix} - 1, x < 0 \\ 0, x = 0 \\ 1 x > 0 \end{matrix}$ Then, the final gradient W _oh is obtained as

(A.13) $\frac{l (d_{t}, {\hat{d}}_{t})}{\partial W_{oh}} = \sum_{j = 1}^{t} (δ_{o, j} h_{j - 1}^{T} + λ sign (W_{oh, j}))$ For the weight gradient of W _ox

(A.14) $\begin{matrix} \frac{l (d_{t}, {\hat{d}}_{t})}{\partial W_{ox, t}} & = \frac{l (d_{t}, {\hat{d}}_{t})}{\partial θ_{o, t}} \frac{\partial θ_{o, t}}{\partial W_{ox, t}} + λ \frac{\partial {| W |}_{1}}{\partial W_{ox, t}} \\ = δ_{o, t} x_{t}^{T} + λ sign (W_{ox, t}) \end{matrix}$

(A.15) $\frac{l (d_{t}, {\hat{d}}_{t})}{\partial W_{ox}} = \sum_{j = 1}^{t} (δ_{o, j} x_{j}^{T} + λ sign (W_{ox, j}))$ For the bias gradient of b _o

(A.16) $\frac{l (d_{t}, {\hat{d}}_{t})}{\partial b_{o, t}} = \frac{l (d_{t}, {\hat{d}}_{t})}{\partial θ_{o, t}} \frac{\partial θ_{o, t}}{\partial b_{o, t}} = δ_{o, t}$

(A.17) $\frac{l (d_{t}, {\hat{d}}_{t})}{\partial b_{o}} = \sum_{j = 1}^{t} δ_{o, j}$ Similar to W _oh, W _ox and b _o, the gradient update equations for W _fh, W _fx, b _f, W _ih, W _ix, b _i, W _ch, W _cx and b _c can also be derived.

Appendix B

Table B4

Optimal hyperparameters of different algorithms on the SO₂ prediction

Hyperparameters	LASSO-MLP	Transformer	LSTM	LASSO-LSTM	STA-LSTM	ℓ₁-HM-GRU	ℓ₁-HM-LSTM
Input dimension	35	35	35	35	35	35	35
Output dimension	1	1	1	1	1	1	1
Number of hidden neurons	40	40	40	40	40	40	40
λ	0.5	–	–	0.04		0.003	0.002
δ	–	–	–	–	–	0.4	0.5
Initial learning rate	0.001	0.001	0.001	0.001	0.001	0.001	0.001
Sliding window m	–	8	8	8	8	8	8
Max epochs	1000	1000	1000	1000	1000	1000	1000

References

Sun

and Ge

, A survey on deep learning for data-driven soft sensors, IEEE Transactions on Industrial Informatics 17(9) (2021), 5853–5866.

Wang

, Bi

, Zhao

, Zhou

, Kim

and Xu

, A study on thermal stress deformation using analytical methods based on the temperature distribution of storage material in a rotary air-preheater, Applied Thermal Engineering 29(11–12) (2009), 2350–2357.

Ren

, Wang

, Laili

and Zhang

, A data-driven self-supervised lstm-deepfm model for industrial soft sensor, IEEE Transactions on Industrial Informatics 18(9) (2021), 5859–5869.

Yin

, Li

, Sun

, Feng

, Yan

and Ma

, Mixed kernel principal component weighted regression based on just-in-time learning for soft sensor modeling, Measurement Science and Technology 33(1) (2021), 015102.

Wang

, Wang

J.S.

, Wang

S.Y.

, Xing

and Li

X.D.

, ANFIS soft sensing model of SMB chromatographic separation process based on new adaptive population evolution particle swarm optimization algorithm, Journal of Intelligent & Fuzzy Systems (2021), 1–26.

Pan

C.C.

, Bai

, Yang

G.K.

, Wong

D.S.H.

and Jang

S.S.

, An inferential modeling method using enumerative PLS based nonnegative garrote regression, Journal of Process Control 22(9) (2012), 1637–1646.

Lai

, Li

and Pan

, A combined model based on feature selection and support vector machine for PM2. 5 prediction, Journal of Intelligent & Fuzzy Systems 40(5) (2021), 10099–10113.

Wang

, Jia

Q.S.

, Zhou

, Bi

, Qiao

and Abusorrah

, Artificial neural networks for water quality soft-sensing in wastewater treatment: a review, Artificial Intelligence Review 55(1) (2022), 565–587.

Guo

, Guo

, Xu

, Xia

and Sun

, MLP neural network-based regional logistics demand prediction, Neural Computing and Applications 33(9) (2021), 3939–3952.

10.

Zheng

, Tian

, Jiang

and Yang

, Layer-wise learning based stochastic gradient descent method for the optimization of deep convolutional neural network, Journal of Intelligent & Fuzzy Systems 37(4) (2019), 5641–5654.

11.

Wang

and Li

, Stochastic configuration networks: fundamentals and algorithms, IEEE Transactions on Cybernetic 47(10) (2017), 3466–3479.

12.

Rahimilarki

, Gao

, Jin

and Zhang

, Convolutional neural network fault classification based on time-series analysis for benchmark wind turbine machine, Renewable Energy 185 (2022), 916–931.

13.

Selva

Kumar

Santhose

and Girirajkumar

, Z-source inverter fed induction motor drive control using particle swarm optimization recurrent neural network, Journal of Intelligent & Fuzzy Systems 28(6) (2015), 2749–2760.

14.

Yuan

, Li

, Shardt

Y.A.

, Wang

and Yang

, Deep learning with spatiotemporal attention-based LSTM for industrial soft sensor model development, IEEE Transactions on Industrial Electronics 68(5) (2020), 4404–4414.

15.

Liu

Y.Q.

, Zhang

, Song

L.H.

and Chen

Y.Y.

, Attention-based recurrent neural networks for accurate short-term and long-term dissolved oxygen prediction, Computers and Electronics in Agriculture 165 (2019), 104964.

16.

Pan

, Su

, Huang

and Wang

, LSTM-based soft sensor design for oxygen content of flue gas in coal-fired power plant, Transactions of the Institute of Measurement and Control 43(1) (2021), 78–87.

17.

Tian

, Xu

, Zhu

Q.X.

and He

Y.L.

, Novel stacked input-enhanced supervised autoencoder integrated with gated recurrent unit for soft sensing, IEEE Transactions on Instrumentation and Measurement 71 (2022) 1–9.

18.

, Wang

and Peng

, A novel bidirectional gated recurrent unit-based soft sensor modeling framework for quality prediction in manufacturing processes, IEEE Sensors Journal 22(19) (2022), 1861–18619.

19.

Vaswani

et al. Attention is all you need., Advances in Neural Information Processing Systems (2017), 30.

20.

Geng

, Chen

, Meng

and Han

, Novel transformer based on gated convolutional neural network for dynamic soft sensor modeling of industrial processes, IEEE Transactions on Industrial Informatics 18(3) (2021), 1521–1529.

21.

Tibshirani

, Regression shrinkage and selection via the LASSO, Journal of the Royal Statistical Society: Series B (Methodological) 58(1) (1996), 267–288.

22.

Sun

, Huang

S.H.

, Wong

D.S.H.

and Jang

S.S.

, Design and application of a variable selection method for multilayer perceptron neural network with LASSO, IEEE Transactions on Neural Networks and Learning Systems 28(6) (2016), 1386–1396.

23.

Fan

, Tao

, Zheng

and Jang

S.S.

, A data-driven soft sensor based on multilayer perceptron neural network with a double LASSO approach, IEEE Transactions on Instrumentation and Measurement 69(7) (2019), 3972–3979.

24.

, Yu

, Xu

, Cao

and Sun

, Nonlinear dynamic soft-sensing modeling of NOx emission of a selective catalytic reduction denitration system, IEEE Transactions on Instrumentation and Measurement 71 (2022), 1–11.

25.

, Nsofor

G.C.

and Song

, A comparative analysis of predictive data mining techniques, International Journal of Rapid Manufacturing 1(2) (2009), 150–172.

26.

Yang

et al. A hybrid robust system considering outliers for electric load series forecasting, Applied Intelligence 52(2) (2022), 1630–1652.

27.

, Huang

and Wang

, Robust stochastic configuration networks with maximum correntropy criterion for uncertain data regression, Information Sciences 473 (2019), 3–86.

28.

Chen

, Li

, Yan

, Guo

and Liu

, Least absolute deviation-based robust support vector regression, Knowledge-Based Systems 131 (2017), 183–194.

29.

Kingma

D.P.

and Ba

, Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980 (2014).

30.

Bergmeir

and Benítez

J.M.

, On the use of cross-validation for time series predictor evaluation, Information Sciences 191 (2012), 192–213.

31.

Friedman

J.H.

, Multivariate adaptive regression splines, The Annals of Statistics 19(1) (1991), 1–67.

32.

Zheng

, Lan

, Zhang

, Kainz

and Chen

, Prediction of MRI RF exposure for implantable plate devices using artificial neural network, IEEE Transactions on Electromagnetic Compatibility 62(3) (2019), 673–681.

33.

Anomaly Detection Toolkit, User Guide. [Online]. Available: https://arundo-adtk.readthedocs-hosted.com/en/stable/userguide.html.

A new regularized and robust LSTM learning algorithm for industrial soft sensors

Abstract

Keywords

1 Introduction

2 Related work

2.1 Long short-term memory neural network

3.1 LSTM with Huber M-Loss and ℓ1 regularization

3.3 Computational flow of proposed algorithm

4 Simulation results on artificial datasets

4.1 Experimental setup

5.1 Process description

Footnotes

Acknowledgments

Funding

Appendix A

Appendix B

References

3.1 LSTM with Huber M-Loss and ℓ₁ regularization