Sparse neural networks with skip-connections for identification of
  aluminum electrolysis cell
Neural networks are rapidly gaining interest in nonlinear system
identification due to the model's ability to capture complex input-output
relations directly from data. However, despite the flexibility of the approach,
there are still concerns about the safety of these models in this context, as
well as the need for large amounts of potentially expensive data. Aluminum
electrolysis is a highly nonlinear production process, and most of the data
must be sampled manually, making the sampling process expensive and infrequent.
In the case of infrequent measurements of state variables, the accuracy and
open-loop stability of the long-term predictions become highly important.
Standard neural networks struggle to provide stable long-term predictions with
limited training data. In this work, we investigate the effect of combining
concatenated skip-connections and the sparsity-promoting $\ell_1$
regularization on the open-loop stability and accuracy of forecasts with short,
medium, and long prediction horizons. The case study is conducted on a
high-dimensional and nonlinear simulator representing an aluminum electrolysis
cell's mass and energy balance. The proposed model structure contains
concatenated skip connections from the input layer and all intermittent layers
to the output layer, referred to as InputSkip. $\ell_1$ regularized InputSkip
is called sparse InputSkip. The results show that sparse InputSkip outperforms
dense and sparse standard feedforward neural networks and dense InputSkip
regarding open-loop stability and long-term predictive accuracy. The results
are significant when models are trained on datasets of all sizes (small,
medium, and large training sets) and for all prediction horizons (short,
medium, and long prediction horizons.)
Erlend Torje Berg Lundby, Haakon Robinsson, Adil Rasheed, Ivar Johan Halvorsen, Jan Tommy Gravdahl
# Sparse neural networks with skip-connections for nonlinear system identificationfootnoteinfo

###### Abstract

Data-driven models such as neural networks are being applied more and more to safety-critical applications, such as the modeling and control of cyber-physical systems. Despite the flexibility of the approach, there are still concerns about the safety of these models in this context, as well as the need for large amounts of potentially expensive data. In particular, when long-term predictions are needed or frequent measurements are not available, the open-loop stability of the model becomes important. However, it is difficult to make such guarantees for complex black-box models such as neural networks, and prior work has shown that model stability is indeed an issue. In this work, we consider an aluminum extraction process where measurements of the internal state of the reactor are time-consuming and expensive. We model the process using neural networks and investigate the role of including skip connections in the network architecture as well as using \(\ell_{1}\) regularization to induce sparse connection weights. We demonstrate that these measures can greatly improve both the accuracy and the stability of the models for datasets of varying sizes.

D +
Footnote †: footnote
2016). Skip-connections were originally proposed by He et al. (2016) as a way to circumvent this, by introducing a shorter path between the early layers and the output. They were not only found to enable the training of significantly deeper networks, but Li et al. (2017) also demonstrated that they may help improve training convergence.

In the field of dynamical systems and control, we often design a model with a purpose in mind, such as the design of a control system or state observer. Crucially, we are interested in the behavior and performance of the controlled system in terms of objectives such as energy efficiency or yield. This implies that the model does not need to be perfectly accurate for the entire state space, so long as the resulting closed-loop performance is sufficient (known as _identification for control_ (I4C)). If high-frequency measurements from the system are available, only the short-term behavior of the model is important, since any drift out of the operational space is quickly corrected. However, if measurements are rarely available, such as in the aluminum electrolysis process that we consider, the long-term model behavior and open-loop stability become much more important. Stable long-term predictions can be important for decision-making, meaning that a model with good long-term stability and accuracy is inherently important.

In this work, we investigate the effects of adding skip connections and \(\ell_{1}\) regularization on the accuracy and stability of these models for short, medium, and long horizons. We address the following questions:

* How do skip connections affect the stability and generalization error of neural networks trained on high-dimensional nonlinear dynamical systems?
* How does sparsity affect stability and generalization error for neural networks with skip connections that model nonlinear dynamics?
* How does the amount of training data affect neural networks with skip connections compared to neural networks without skip connections?

We make the following contributions:

* We perform a black box system identification of an aluminum electrolysis cell using different NN architectures.
* We demonstrate that the accuracy and open-loop stability of the resulting models is greatly improved by using \(\ell_{1}\) weight regularization and incorporating skip connections into the architecture.
* This advantage is consistent across datasets of varying sizes.

## 2 Theory

### Physics-based model for aluminum extraction

We evaluate NNs for nonlinear system identification by first training them on synthetic data generated from a known physics-based models (PBM). The model used in this work describes the internal dynamics of an aluminum electrolysis cell based on the Hall-Heroult process. Figure 1 shows a diagram of the electrolysis cell. Traditional PBMs of such systems are generally constructed by studying the mass/energy balance of the chemical reactions. Lundby et al. (2022) presents a more detailed exposition of the model that we use in this work. The system is described by a set of ordinary differential equations (ODE):

\[\dot{\mathbf{x}}=\mathbf{f}(\mathbf{x},\mathbf{u}), \tag{1}\]

where \(\mathbf{x}\in\mathbb{R}^{8}\) and \(\mathbf{u}\in\mathbb{R}^{5}\) represent the time-varying states and inputs of the system respectively. The full set of equations are:

\[\dot{x}_{1} =\frac{k_{1}(g_{1}-x_{7})}{x_{1}k_{0}}-k_{2}(x_{6}-g_{1}) \tag{2a}\] \[\dot{x}_{2} =\,u_{1}-k_{3}u_{2}\] (2b) \[\dot{x}_{3} =\,u_{3}-k_{4}u_{1}\] (2c) \[\dot{x}_{4} =\,-\frac{k_{1}(g_{1}-x_{7})}{x_{1}k_{0}}+k_{2}(x_{6}-g_{1})+k_{ 5}u_{1}\] (2d) \[\dot{x}_{5} =\,k_{6}u_{2}-u_{4}\] (2e) \[\dot{x}_{6} =\,\frac{\alpha}{x_{2}+x_{3}+x_{4}}\Bigg{[}u_{2}g_{5}+\frac{u_{ 2}^{2}u_{5}}{2620g_{2}}-k_{7}(x_{6}-g_{1})^{2}\] (2f) \[+k_{8}\frac{(x_{6}-g_{1})(g_{1}-x_{7})}{k_{0}x_{1}}-k_{9}\frac{x_ {6}-x_{7}}{k_{10}+k_{11}k_{0}x_{1}}\Bigg{]}\] \[\dot{x}_{7} =\,\frac{\beta}{x_{1}}\Bigg{[}\frac{k_{9}(g_{1}-x_{7})}{k_{15}k_{ 0}x_{1}}-k_{12}(x_{6}-g_{1})(g_{1}-x_{7})\] (2g) \[+\,\frac{k_{13}(g_{1}-x_{7})^{2}}{k_{0}x_{1}}\,-\frac{x_{7}-x_{8}} {k_{14}+k_{15}k_{0}x_{1}}\Bigg{]}\] \[\dot{x}_{8} =\,k_{17}k_{9}\left(\frac{x_{7}-x_{8}}{k_{14}+k_{15}k_{0}\cdot x _{1}}-\,\frac{x_{8}-k_{16}}{k_{14}+k_{18}}\right), \tag{2h}\]

where the intrinsic properties \(g_{i}\) of the bath mixture are given as:

Figure 1: Schematic of the setup
\[g_{1} =991.2+112c_{x_{3}}+61c_{x_{3}}^{1.5}-3265.5c_{x_{3}}^{2.2} \tag{3a}\] \[-\frac{793c_{x_{2}}}{-23c_{x_{2}}c_{x_{3}}-17c_{x_{3}}^{2}+9.36c_{x_ {3}}+1}\] \[g_{2} =\exp\,\left(2.496-\frac{2068.4}{273+x_{6}}-2.07c_{x_{2}}\right)\] (3b) \[g_{3} =0.531+3.06\cdot 10^{-18}u_{1}^{3}-2.51\cdot 10^{-12}u_{1}^{2}\] (3c) \[+6.96\cdot 10^{-7}u_{1}-\frac{14.37(c_{x_{2}}-c_{x_{2},crit})-0.431 }{735.3(c_{x_{2}}-c_{x_{2},crit})+1}\] \[g_{4} =\frac{0.5517+3.8168\cdot 10^{-6}u_{2}}{1+8.271\cdot 10^{-6}u_{2}}\] (3d) \[g_{5} =\frac{3.8168\cdot 10^{-6}g_{3}g_{4}u_{2}}{g_{2}(1-g_{3})}. \tag{3e}\]

See Table 1 for a description of these quantities.

The values of these constants can be found in Lundby et al. (2022). The dynamics of the system are relatively slow. The control inputs \(u_{1},\ u_{3}\) and \(u_{4}\) are therefore well modeled as impulses that represent discrete events involving the addition or removal of substances. This results in step changes in the linear states \(x_{2},x_{3},x_{5}\), which act as accumulator states for the mass of the corresponding substance (see Table 1). The control inputs \(u_{2}\) and \(u_{5}\) are piecewise constant, and always nonzero. The inputs \(\mathbf{u}\) are determined by a simple proportional controller \(\boldsymbol{\pi}(\mathbf{x})\). The simulation model is derived in Lundby et al. (2022), and we refer to that article for further details.

### Deep neural network with skip connections

A NN with \(L\) layers can be compactly written as an alternating composition of affine transformations \(\mathbf{Wz}+\mathbf{b}\) and nonlinear activation functions \(\boldsymbol{\sigma}:\mathbb{R}^{n}\mapsto\mathbb{R}^{n}\):

\[\hat{\mathbf{f}}(\mathbf{z})=\hat{\mathbf{f}}_{L}\circ\cdots\circ\hat{ \mathbf{f}}_{2}\circ\hat{\mathbf{f}}_{1} \tag{4}\]

where the activation function \(\boldsymbol{\sigma}_{i}\), weight matrix \(\mathbf{W}_{i}\), and bias vector \(\mathbf{b}_{i}\) correspond to the \(i\)th layer of the network. The universal approximation property of NNs makes them very attractive as a flexible model class when a lot of data is available. The representation capacity is generally understood to increase with both the depth and the width (the number of neurons in each layer), although early attempts to train very deep networks found them challenging to optimize using backpropagation due to the vanishing gradients problem. One of the major developments that enabled researchers to train deep NNs with many layers is the _skip connection_. A skip connection is simply an additional inter-layer connection that bypasses some of the layers of the network. This provides alternate pathways through which the loss can be backpropagated to the early layers of the NN, which helps mitigate the issues of vanishing and exploding gradients, which were major hurdles to training deeper models. In this work, we utilize a modified DenseNet architecture as proposed by Huang et al. (2017), where the outputs of earlier layers are concatenated to all the consecutive layers. We simplify the structure such that the model only contains skip connections from the input layer to all consecutive layers. We call this architecture InputSkip, which has reduced complexity compared to DenseNet. This design is motivated by the fact that the output of each layer (including the final output) becomes a sum of both a linear and a nonlinear transformation of the initial input \(\mathbf{x}\). Hence, the skip connections from the input layer to consecutive layers facilitate the reuse of the input features for modeling different linear and nonlinear relationships more independently of each other.

## 3 Method and setup

In this section, we present all the details of data generation and its preprocessing, and the methods that are required to reproduce the work. The steps can be briefly summarized as follows:

* Use Equation (2) with random initial conditions to generate 140 trajectories with 5000 timesteps each. Set aside 40 for training and 100 for testing. Construct 3 datasets by selecting 10,20 and 40 trajectories respectively.
* For each model class and dataset, train 10 instances on the training data.
* Repeat all experiments with \(\ell_{1}\) regularization, see loss function in Equation (5).
* Use trained models to generate predicted trajectories along the test set and compare them to the 100 test trajectories.

### Data generation

Equation (2) was discretized using the RK4 scheme with a fixed timestep \(h=10\,\mathrm{s}\) and numerically integrated on the interval \([0,5000h]\). We used uniformly randomly sampled initial conditions from the intervals shown in Table 2 to generate 140 unique trajectories. We set aside 40 trajectories for training and 100 of the trajectories as a test set. The 40 training trajectories were used to create 3 datasets of varying sizes (small, medium, large), namely 10, 20, and 40 trajectories. In total, the datasets contained 50000, 100000, and 200000 individual data points respectively.

Equation (2) also depends on the input signal \(\mathbf{u}\). In practice, this is given by a deterministic control policy \(\mathbf{u}=\boldsymbol{\pi}(\mathbf{x})\) that stabilizes the system and keeps the state \(\mathbf{x}\) within some region of the state space that is suitable for safe operation. We found that this was insufficient to successfully train our models, because the controlled trajectories showed very little variation after some time, despite having different initial conditions. This lack of diversity in the dataset resulted in models that could not generalize to unseen states, a situation that frequently arose during evaluation. To inject more variety into the data and sample states \(\mathbf{x}\) outside of the standard operational area, we used a stochastic controller

\[\boldsymbol{\pi}_{s}(\mathbf{x})=\boldsymbol{\pi}(\mathbf{x})+\mathbf{r}(t)\]

that introduced random perturbations \(\mathbf{r}(t)\) to the input. These perturbations were sampled using the Amplitude-modulated Pseudo-Random Binary Signal (APRBS) method proposed by Winter and Breitsamter (2018) for nonlinear system identification.

In system identification it is typical to optimize the model to estimate the function \(\hat{\mathbf{x}}=\mathbf{f}(\mathbf{x},\mathbf{u})\). However, this is not feasible for Equation (2) because the inputs \(\mathbf{u}\) are not differentiable. Instead, we discretize the trajectories using the forward Euler difference and use this as the regression variable:

\[\mathbf{y}_{k}=\frac{\mathbf{x}_{k+1}-\mathbf{x}_{k}}{h}\]

The datasets are then constructed as sets of the pairs \(([\mathbf{x}_{k},\mathbf{u}_{k}],\mathbf{y}_{k})\).

### Training setup

We optimize the models by minimizing the following loss function using stochastic gradient descent:

\[\mathbf{J}_{\theta}=\frac{1}{|\mathcal{B}|}\sum_{i\in\mathcal{B}}(\mathbf{y}_{ i}-\mathbf{\hat{f}}(\mathbf{x}_{i},\mathbf{u}_{i}))^{2}+\lambda\sum_{j=1}^{L}| \mathbf{W}_{j}| \tag{5}\]

where \(\mathcal{B}\) is a _batch_ of randomly sampled subset of indices from the dataset, \(L\) is the number of layers of the NN, and \(\lambda\) is the regularization parameter. This loss function is the sum of the mean squared error (MSE) of the model \(\hat{\mathbf{f}}\) with respect to the regression variables \(\mathbf{y}\), and the \(\ell_{1}\) norm of the connection weight matrices \(\mathbf{W}_{i}\) in all layers. We used a batch size of \(|\mathcal{B}|=128\). We used the popular ADAM solver proposed by Kingma and Ba (2014) with default parameters to minimize Equation (5).

### Evaluation of model accuracy

As previously mentioned, we are interested in evaluating the long-term predictive accuracy of the models. Starting from a given initial condition \(\mathbf{x}(t_{0})\), the model \(\mathbf{\hat{f}}(\mathbf{x},\mathbf{u})\) is used to generate an estimated trajectory using the recurrence:

\[\hat{\mathbf{x}}_{k+1}=\hat{\mathbf{x}}_{k}+h\,\hat{\mathbf{f}}(\hat{\mathbf{ x}}_{k},\mathbf{u}_{k}) \tag{6}\]

where \(\hat{\mathbf{x}}_{0}=\mathbf{x}_{0}\). Note that the input signal \(\mathbf{u}_{k}\) is replayed directly from the test trajectory. Borrowing a term from the field of time-series analysis, we refer to this as a _rolling forecast_. To evaluate the accuracy of a model over multiple trajectories, we define the Average Normalized Rolling Forecast Mean Squared Error (AN-RFMSE):

\[\text{AN-RFMSE}=\frac{1}{p}\sum_{i=1}^{p}\frac{1}{n}\sum_{j=1}^{n}\left(\frac{ \hat{x}_{i}(t_{j})-x_{i}(t_{j})}{\text{std}(x_{i})}\right)^{2}, \tag{7}\]

where \(\hat{x}_{i}(t_{j})\) is the model estimate of the simulated state variable \(x_{i}\) at time step \(t_{j}\), \(\text{std}(x_{i})\) is the standard deviation of variable \(x_{i}\) in the training set \(\mathcal{S}_{train}\), \(p=8\) is the number of state variables and \(n\) is the number of time steps being averaged over.

### Evaluation of model stability

A symptom of model instability is that its predictions can _blow-up_, which is characterized by a rapid (often exponential) increase in prediction error. More precisely, we say that a blow-up occurs when the normalized mean absolute error for all system states exceeds three (this corresponds to standard deviations). We detect this as follows:

\[\max_{j<n}\left[\frac{1}{p}\sum_{i=1}^{p}\left(\frac{|\hat{x}_{i}(t_{j})-x_{i} (t_{j})|}{\text{std}(x_{i})}\right)\right]>3 \tag{8}\]

where \(p=8\) is again the number of state variables and \(n\) is the number of time steps to consider. This is a conservative estimate. However, this does not lead to any significant underestimation of the number of blow-ups. This is because once a model starts to drift rapidly, it very quickly exceeds the normal error of three standard deviations.

\begin{table}
\begin{tabular}{l|l} \hline Variable & Initial condition interval \\ \hline \(x_{1}\) & [2060, 4460] \\ \(c_{x_{2}}\) & [0.02, 0.05] \\ \(c_{x_{3}}\) & [0.09, 0.13] \\ \(x_{4}\) & [11500, 16000] \\ \(x_{5}\) & [9550, 10600] \\ \(x_{6}\) & [940, 990] \\ \(x_{7}\) & [790, 850] \\ \(x_{8}\) & [555, 610] \\ \hline \end{tabular}
\end{table}
Table 2:
## 4 Results and Discussions

We characterize the different model classes (PlainDense, PlainSparse, InputSkipDense, InputSkipSparse) by estimating their blow-up frequencies and their rolling forecast mean squared error (RFMSE) on the validation data. The blow-up frequency is an interesting measure since it can indicate how stable the model is in practice.

We perform a Monte Carlo analysis by training 10 instances of each model class and evaluating these on 100 trajectories randomly generated using the true model, yielding 1000 data points for each model class. We repeat the experiments for 3 different dataset sizes to study the data efficiency of the models.

Figure 2 presents the total number of blow-ups recorded within each model class after \(100h\), \(2000h\), and \(5000h\) (short, medium, and long term respectively). For simplicity, blow-ups were detected by thresholding the computed variance of a predicted trajectory and manually inspected. It is clear that for short time horizons all the models exhibit robust behavior independently of the size of the training datasets. However, for medium and long time horizons, PlainDense, PlainSparse, and InputSkipDense architectures exhibit a significant number of blow-ups and therefore instability. Figure 1(a) - 1(c) show that PlainDense is generally the most unstable, with up to 67% of all trajectories resulting in a blow-up. For the smallest amount of training data (Figure 1(a)) PlainSparse and InputSkipDense have similar blow-up frequencies. For larger datasets, the PlainSparse architecture shows significantly better stability than both PlainDense and InputSkipDense. InputSkipDense and PlainDense both show better stability with increasing amounts of training data in terms of fewer blow-ups. However, both these dense models still suffer from significant amounts of blow-ups.

In comparison, almost no blow-ups are recorded when using the InputSkipSparse architecture, even for the small training dataset. In Figure 2, the orange bars corresponding to the blow-up frequency of InputSkipSparse models are not visible for any of the training sets due to the significantly lower number of blow-ups. For InputSkipSparse models trained on the smallest dataset, only 3 out of 1000 possible blow-ups were reported for the longest horizon. Apart from that, no blow-ups were reported for the InputSkipSparse models.

Only a few blow-ups were recorded after \(5000h\) in the medium term.

Figure 3 presents a violin plot of the accuracy of each model class, expressed in terms of RFMSE over different time horizons. Only the plot for the smallest dataset (50000 points) is shown, due to the results being very similar. A larger width of the violin indicates a higher density of that given RFMSE value, while the error bars show the minimum and maximum recorded RFMSE values. The model estimates that blew up (see Figure 2) are not included. In this way, we estimate the generalization performance of the models only within their regions of stability. Note that the violin plots for model classes with many blow-ups are made using fewer samples, and can be seen as slightly "cherry-picked". Nonetheless, the InputSkipSparse architecture consistently yields more accurate results, up to an order of magnitude better than the others in the long term.

## 5 Conclusion and Future Work

In this work, we compared the performance of two different model structures trained both with and without sparsity promoting \(\ell_{1}\) regularization. The two model types are standard Multi-Layer Perceptrons (MLP), and a more specialized architecture that includes skip connections from the input layer to all consecutive layers. This yields four different model structures, which we call PlainDense, PlainSparse, InputSkipDense, and InputSkipSparse. The main conclusions of the article are as follows:

* NNs with skip connections are more stable for predictions over long time horizons compared to standard MLPs. Furthermore, the accuracy of NNs with skip

Figure 2: Divergence plot: Number of trajectories that blow-up over different time horizons. The total number of trajectories is 1000, so the values can be read as a permille.
connections is consistently higher for all forecasting horizons.
* The application of sparsity-promoting \(\ell_{1}\) regularization significantly improves the stability of both the standard MLP and InputSkip architectures. This improvement was more apparent for models with the InputSkip architecture.
* The InputSkipSparse showed satisfactory stability characteristics even when the amount of training data was restricted. This suggests that this architecture is more suitable for system identification tasks than the standard MLP structure.

The case study shows that both sparsity-promoting regularization and skip connections can result in more stable NN models for system identification tasks while requiring less data, as well as improving their multi-step generalization for both short, medium, and long prediction horizons. Despite the encouraging performance of the sparse-skip networks, we can not guarantee similar performance for noisy data, as we have only investigated the use of synthetic data devoid of any noise. However, such a study will be an interesting line of future work. This case study also has relevance beyond the current setup. In more realistic situations, we often have a partial understanding of the system we wish to model (see Equation (2)), and only wish to use data-driven methods to correct a PBM when it disagrees with the observations (e.g. due to a faulty assumption). As shown in Robinson et al. (2022), combining PBMs and data-driven methods in this way also has the potential to inject instability into the system. Finding new ways to improve or guarantee out-of-sample behavior for data-driven methods is therefore of paramount importance to improve the safety of such systems.

## Acknowledgements

This work was supported by the industry partners Borgregaard, Elkem, Hydro, Yara and the Research Council of Norway through the projects TAPI: Towards Autonomy in Process Industries (grant no. 294544) and EXAIGON: Explainable AI systems for gradual industry adoption (grant no. 304843)
