# Sparse deep neural networks for modeling aluminum1 electrolysis dynamics

Erlend Torje Berg Lundby11footnotemark: 1<EMAIL_ADDRESS>Adil
Rasheed<EMAIL_ADDRESS>Jan Tommy Gravdahl<EMAIL_ADDRESS>Ivar Johan Halvorsen<EMAIL_ADDRESS>Norwegian University of
Science and Technology, Department of Engineering Cybernetics, Trondheim,
Norway  SINTEF Digital, Department of Mathematics and Cybernetics, Trondheim,
Norway

###### Abstract

Artificial neural networks have a broad array of applications today due to
their high degree of flexibility and ability to model nonlinear functions from
data. However, the trustworthiness of neural networks is limited due to their
black-box nature, their poor ability to generalize from small datasets, and
their inconsistent convergence during training. Aluminum electrolysis is a
complex nonlinear process with many interrelated sub-processes. Artificial
neural networks can potentially be well suited for modeling the aluminum
electrolysis process, but the safety-critical nature of this process requires
trustworthy models. In this work, sparse neural networks are trained to model
the system dynamics of an aluminum electrolysis simulator. The sparse model
structure has a significantly reduction in model complexity compared to a
corresponding dense neural network. We argue that this makes the model more
interpretable. Furthermore, the empirical study shows that the sparse models
generalize better from small training sets than dense neural networks.
Moreover, training an ensemble of sparse neural networks with different
parameter initializations show that the models converge to similar model
structures with similar learned input features.

###### keywords:

Aluminum electrolysis , Sparse neural networks , Data-driven modeling ,
Nonlinear dynamics , Ordinary differential equations

††journal: Arxiv

##  1 Introduction/ motivation

Aluminum is extracted from the electrolytic Hall-Héroult process. In this
process, alumina $(Al_{2}O_{3})$ is dissolved in cryolite $(Na_{3}AlF_{6})$
and then reduced to aluminum. The reaction is driven by a line current that is
sent through the electrolytic bath [1]. The dynamics of the aluminum
electrolysis process contain nonlinearities and interrelated sub-processes
which are challenging to model. Accurate dynamic models are crucial for
optimizing product quality and energy consumption. Furthermore, the harsh
environment of aluminum electrolysis cells requires extra effort to ensure
safe operation. To this end, research efforts have focused on modeling the
Hall-Héroult process accurately. Most of the models developed until recently
falls in the first principle physics-based models (PBM). For example, authors
in [2] developed a model of mass and energy balance based on the first law of
thermodynamics. The resulting model includes, among other things, a complete
control volume analysis, an extensive material balance, a 3D finite element
model (FEM) for modeling resistance in the cell lining and shell, and
Computational Fluid Dynamics (CFD) model for calculating gas velocity
streamlines. In [3], a multi-scale, multi-physics modeling framework including
magneto-hydrodynamics, bubble flow, thermal convection, melting and
solidification based on a set of chemical reactions was developed. These kinds
of PBMs, although highly interpretable, cannot capture all the complex
dynamics of the aluminum electrolysis due to the fact that they have numerous
assumptions, incomplete physics, and uncertainties in the input parameters.
For example, the magneto-hydrodynamic phenomena or the reactivity and species
concentration distribution are phenomena that are difficult to model [4]. Thus
these first-principle model’s predictions might deviate significantly from
observation.

In recent year owing to the availability of copious amount of data, cheap
computational resources, and major advancement in the algorithms, Data-Driven
Modeling (DDM) is coming up as an attractive alternative to PBM. The DDM has
the potential to accurately model even the poorly understood complex phenomena
(mentioned earlier) directly from the data. As a result, a broad array of
scientific communities explore their applicability in many engineering
applications. Some examples are the use of neural networks to model simulated
dynamics of a pressurized water nuclear reactor [5], identification of the
dynamics of the production and purification process of bioethanol [6],
prediction of chemical reactions [7]. Authors in [8] demonstrated convincingly
how DDM outperforms PBM in the absence of the full understanding of physics.
However, despite their advantages they suffer from certain shortcomings [9];
they are difficult to interpret, difficult to generalize to solve previously
unseen problems, and unstable to train. These are important shortcomings to
overcome before the models can be used in high stake applications. We discuss
each of these briefly.

Interpretability: This can be defined as the ability of a model to express
itself in human interpretable form [10]. A simple model like linear regression
having very few trainable parameters can be a good example of an interpretable
model. However, a DNN, which can be seen as a complex version of linear
regression, and constituting of millions of trainable parameters can be
extremely difficult or almost impossible to interpret. Reducing the complexity
of DNN models can lead to more intepretable models.

Generalizaibility: This refers to the models’ ability to predict outcome
values for unseen data from the same distribution as the training data. Highly
complex models are prone to overfitting, meaning that they do not generalize
to areas of the input space that training data do not cover. DNN models have
shown promising results in prediction accuracy by processing large amounts of
data. However, in many complex physical systems, such as aluminum
electrolysis, the cost, and challenges of data acquisition limit the training
data to partial information about the system dynamics. As a result, many DNN
models fail to generalize from limited data, and will provide inaccurate
predictions in regions poorly represented by the training data. Reducing the
complexity of the DNN can increase the generalizability of the model.

Stability: The training of DNN requires solving an optimization problem where
a cost function dependent on thousands or even million of parameters is
minimized. Depending upon the complexity of the problems multiple minima can
be encountered during the optimization, and the task of converging to the
global minima gets increasingly difficult with the number of parameters. By
stability in the context of the current work we mean that the optimization
leads to a reasonable minima which even if not global is close to it. To the
contrary an unstable DNN will be the one where the optimization process gets
stuck in a bad local minima.

Combining PBM and DDM in a hybrid modeling approach to address these
shortcomings is emerging for a wide range of scientific applications [9]. In
[8], the equations known from PBM are augmented by a corrective source term
generated by a DNN. In [11], the estimate calculated by a PBM is subtracted
from coarse measurements of the metal height in an aluminum electrolysis
simulator to find residual measurements. These residual measurements are used
by a compressed sensing algorithm to model the unknown disturbance signal of
the dynamics. In [12], it is proposed a novel hybrid modeling approach called
Physics Informed Neural Network (PINN). PINN utilizes known first principle
knowledge expressed in partial differential equations and their corresponding
boundary conditions to regularize a neural network that is trained to estimate
certain variables. Authors of [13] introduces a method called Physics Guided
Machine Learning (PGML). The main idea is to augment simplified theories
relevant to the dynamics of the system that is modeled with the learning
process of the resulting hybrid DNN model. The simplified theories are
considered as features in the resulting model. These features are added to
intermediate layers in the neural network architecture rather than to the
input layer in order to emphasize their importance. As most of the examples
above illustrates, neural networks are often highly important in hybrid
models. While hybrid models themselves addresses the shortcomings mentioned
above, DNNs utilized in hybrid models still suffer from them. This article
focuses on addressing interpretability, generalizability and training
stability of DNNs used in a purely DDM manner. This can also highly benefit
hybrid models.

In the aluminum electrolysis process, DNNs have been applied to predict
essential variables that are difficult to measure continuously. In [14] and
[15] neural networks were used to predict the temperature and the
concentration of alumina, aluminum fluoride, and calcium fluoride in the bath.
In [16], neural networks were used to model the carbon anode properties as a
function of physical and chemical properties. Despite the cited works and many
more available in the existing literature, issues related to interpretability,
generalizability and training stability of DNN models in the context of the
aluminum electrolysis process remains unaddressed.

To this end, in the current article, we apply the sparsity promoting
$\ell_{1}$ regularizer on the weights of each layer in DNN models that predict
the states in the aluminum electrolysis. The regularization pushes most of the
weights to zero, thus reducing the complexity drastically. We demonstrate that
the resulting sparse models have increased interpretability, generalizability
and training stability.

The article is structured as follows. Section 2 presents the relevant theory
for the work in the case study. Section 3 presents the method applied in the
paper and the experimental setup of the simulator for data generation. In
section 4, the results are presented and discussed. Finally, in section 5,
conclusions are given, and potential future work is presented.

##  2 Theory

###  2.1 Deep neural networks (DNN)

A DNN is a supervised machine learning algorithm [17] that can be denoted by

$\mathbf{y}=\mathbf{\hat{f}}(\mathbf{x};\bm{\theta}),$ (1)

where $\mathbf{y}\in\mathbb{R}^{s}$ is the output vector of the network model
and $s$ is the length of the output vector of the model.
$\mathbf{x}\in\mathbb{R}^{d}$ is the input vector to the network model, and
$d$ is the input dimension. Here, $\bm{\theta}\in\mathbb{R}^{p}$ denotes all
trainable parameters in the network model. $p$ is the number of parameters.
Each layer j+1 operates on the output vector from the previous layer
$\mathbf{Z}^{j}\in R^{L_{j}}$ and outputs a vector
$\mathbf{Z}^{j+1}\in\mathbb{R}^{L_{j+1}}$:

$\mathbf{Z}^{j+1}=\sigma(\mathbf{W}^{j+1}\mathbf{Z}^{j}+\mathbf{b}^{j+1}).$
(2)

$\mathbf{W}^{j+1}\in\mathbb{R}^{L_{j+1}\times L_{j}}$ is called the weight
matrix, and $\mathbf{b}^{j+1}\in\mathbb{R}^{L_{j+1}}$ is the called the bias
vector of layer $j+1$.
$\bm{\theta}=\\{\bm{\theta^{1}},\;...,\;\bm{\theta}^{j+1},\;...\\}$, and
$\bm{\theta}^{j+1}=\\{\mathbf{W}^{j+1},\;\mathbf{b}^{j+1}\\}$. $\sigma$ is a
non-linear activation function. That is,

$\displaystyle\sigma(\mathbf{W}^{j+1}\mathbf{Z}^{j}+\mathbf{b}^{j+1})=(\sigma(W^{j+1}_{1}\mathbf{Z}^{j}+b^{j+1}_{1}),\;...,$
$\displaystyle\;\sigma(W^{j+1}_{i}\mathbf{Z}^{j}+b^{j+1}_{i}),\;...,\;\sigma(W^{j+1}_{L_{j+1}}\mathbf{Z}^{j}+b^{j+1}_{L_{j+1}}))^{T}.$
(3)

$W^{j+1}_{i}$ are the row vectors of the weight matrix $\mathbf{W}_{j+1}$ and
$b^{j+1}_{i},\;i=1,\;...,\;L_{j+1}$ are the entries of the bias vector
$\mathbf{b}^{j+1}$. Thus, the activation function calculates the output of
each neuron in layer $j+1$ as a nonlinear function of the weighted sum of
outputs from the neurons in the previous layer plus a bias. Each neuron
outputs one value, and the weight in the consecutive layer determines the
importance of the output of each neuron in the current layer. The nonlinear
activation function $\sigma$ can, for example, be the sigmoid function,
hyperbolic tangent function or the binary step function to mention a few. For
the last decade or so, the popularity of the piece-wise linear (PWL)
activation function Rectified Linear Unit (ReLU) has grown exponentially. This
is in part due to its computational simplicity, representational sparsity and
non-vanishing gradients. The ReLU activation function is given by:

$\sigma(z)=max\\{0,\;z\\}.$ (4)

ReLU is the only activation function used in the current work.

###  2.2 Sparse neural networks and regularization

Dense neural networks are often overparameterized models, meaning that they
have more parameters than can be estimated from the data and thus often suffer
from overfitting. In [18], it is shown empirically that randomly initialized
dense neural networks contain subnetworks that can improve generalization
compared to the dense networks. These subnetworks, characterized by
significantly fewer non-zero trainable parameters than their dense
counterpart, are called sparse neural networks. Their utility can further be
seen in terms of increased computational performance for inference and
training, and increased storage and energy efficiency. Typically large-scale
models that require millions to billions of parameters and arithmetic
operations can highly benefit from such sparsification. To conclude
sparsification of complex models will lead to simpler models which are
relatively easier to interpret, generalize, and train.

Among the methods that can be used to sparsity a complex network,
regularization techniques are the most popular ones. In regularization,
penalty terms $R(\mathbf{w})$ defined on the weights are added to the cost
function $C$:

$C(\mathbf{x}_{i},\mathbf{y}_{i},\bm{\theta})=L(\mathbf{y}_{i},\mathbf{\hat{f}}(\mathbf{x};\bm{\theta}))+\lambda
R(\mathbf{w}).$ (5)

The vector $\bm{\theta}=\\{\mathbf{w},\;\mathbf{b}\\}$ denotes the adaptable
parameters, namely the weights $\mathbf{w}$ and biases $\mathbf{b}$ in the
network model. Furthermore, $\\{(\mathbf{x}_{i},\mathbf{y}_{i})\\}_{i=1}^{N}$
is the training data, $L(\cdot,\cdot)$ is a proper loss function and the
positive scalar coefficient $\lambda$ is a hyperparameter to weight the terms
$L(\cdot,\cdot)$ and $R(\cdot)$. The standard choice of loss function
$L(\cdot,\cdot)$ for regression tasks is the Mean Squared Error (MSE). In the
training process, the cost function is minimized to find optimal values of the
parameters:

$\displaystyle\bm{\theta}^{*}=\underset{\bm{\theta}}{\mathrm{argmin}}\left\\{\frac{1}{N}\sum_{i=1}^{N}C(\mathbf{x}_{i},\mathbf{y}_{i},\bm{\theta})\right\\}.$
(6)

The most intuitive sparsity promoting regularizer is the $\ell_{0}$ norm,
often referred to as the sparsity norm:

$R_{\ell_{0}}(\mathbf{w})=||\mathbf{w}||_{0}=\sum_{i}\begin{cases}1&w_{i}\neq
0,\\\ 0&w_{i}=0.\end{cases}$ (7)

The $\ell_{0}$ norm counts the number of nonzero weights. Unfortunately, the
$\ell_{0}$ norm has several drawbacks that make it less suitable for
optimization. The $\ell_{0}$ norm is nondifferentiable and sensitive to
measurement noise. Furthermore, in terms of computational complexity, the
problem of $\ell_{0}$ norm is shown to be NP-hard[19]. The $\ell_{1}$ norm is
a convex relaxation of the $\ell_{0}$ norm, and is given by:

$R_{\ell_{1}}(\mathbf{w})=||\mathbf{w}||_{1}=\sum_{i}|w_{i}|.$ (8)

Due to its geometrical characteristics, $\ell_{1}$ minimization is sparsity
promoting. However, the $\ell_{1}$ norm usually does not reduce the weights to
zero but rather to very small magnitudes. Thus, magnitude pruning can be
applied after $\ell_{1}$ minimization to achieve true sparse models.

Figure 1: Illustration of $\ell_{1}$ regularization. $\beta_{1}$ and
$\beta_{2}$ are the model parameters. $\hat{\beta}$ is the MSE estimate. The
ellipses show the contours of the error from the MSE estimate. The blue
diamond illustrates the $\ell_{1}$ constraints.

Fig. 1 illustrates how $\ell_{1}$ regularization constraints the parameters to
a sparse solution. In this case $\beta_{1}$ is pushed towards zero.

###  2.3 Region bounds for PWA neural networks

The complexity of neural networks with Piecewise Affine (PWA) activation
functions such as ReLU can be analyzed by looking at how the network
partitions the models’ input space to an exponential number of linear response
regions [20, 21]. For each region in the input space, the PWA neural network
has a linear response for the output. Authors in [22] present asymptotic upper
and lower bounds for maximum number of regions for a DNN with ReLU activation:

$\displaystyle\begin{aligned}
\textrm{Lower}&:\Omega\left(\left(\frac{n}{d}\right)^{(L-1)d}n^{d}\right),\\\
\textrm{Upper}&:\mathcal{O}\left(n^{dL}\right).\\\ \end{aligned}$ (9)

$d$ is the input dimension, $L$ is the number of hidden layers, and $n$ is the
number of activation functions or neurons in each layer. The bounds in Eq. (9)
are valid for networks with the same number of neurons in each layer. The
bonds for networks with an uneven number of neurons show similar exponential
results and are thus not included for convenience. Eq. (9) illustrates the
exponential relation between the input dimension, number of neurons, and the
depth of the network. For realistic amounts of data sampled from a physical
system, the number of linear regions that a relatively small dense neural
network partition the input space into exceeds the sampled data by several
orders of magnitude. Thus, in order to generalize to larger areas of the
models’ input space, the number of regions needs to be reduced drastically.
This motivates sparsifying the model.

###  2.4 Simulation model

The simulation model used in this article is based on the mass and energy
balance of an aluminum electrolysis cell. The derivation of the model is found
in the A. In this section, the simulation model is put in a state space form,
where constants, system states, and control inputs are denoted by by symbols
$k_{i}$, $x_{i}$, or $u_{i}$ respectively. The simulation model can be
expressed as a nonlinear system of ODE’s with $8$ states
$\mathbf{x}\in\mathbb{R}^{8}$ and $5$ inputs $\mathbf{u}\in\mathbb{R}^{5}$ on
the form:

$\dot{\mathbf{x}}=f(\mathbf{x},\mathbf{u}),$ (10)

where $\dot{\mathbf{x}}\in\mathbb{R}^{8}$ is the time derivative of the states
$\mathbf{x}$, and $f(\mathbf{x},\mathbf{u})$ is a nonlinear function. Table 1
presents the mapping from physical quantities to variables.

Table 1: Table of states and inputs Variable | Physical meaning | Unit
---|---|---
$x_{1}$ | mass side ledge | $kg$
$x_{2}$ | mass $Al_{2}O3$ | $kg$
$x_{3}$ | mass $ALF_{3}$ | $kg$
$x_{4}$ | mass $Na_{3}AlF_{6}$ | $kg$
$x_{5}$ | mass metal | $kg$
$x_{6}$ | temperature bath | ${}^{\circ}C$
$x_{7}$ | temperature side ledge | ${}^{\circ}C$
$x_{8}$ | temperature wall | ${}^{\circ}C$
$u_{1}$ |  $Al_{2}O_{3}$ feed | $kg$
$u_{2}$ | Line current | $kA$
$u_{3}$ |  $AlF_{3}$ feed | $kg$
$u_{4}$ | Metal tapping | $kg$
$u_{5}$ | Anode-cathode distance | $cm$

The nonlinear functions in Eqs. 11) - (15 that partly describes the dynamics
of the system states are defined in advance of presenting the system dynamics
in order to simplify the expressions in Eq. 10:

$\displaystyle g_{1}$
$\displaystyle=\,991.2+1.12pr_{x_{3}}-0.13pr_{x_{3}}^{2.2}+0.061pr_{x_{3}}^{1.5}$
$\displaystyle\,-\frac{7.93pr_{x_{2}}}{1+pr_{x_{3}}(0.094-0.0017pr_{x_{3}}-0.0023pr_{x_{2}})}$
(11) $\displaystyle g_{2}$
$\displaystyle=\,exp(2.496-\frac{2068.4}{273+x_{6}}-2.07c_{x_{2}})$ (12)
$\displaystyle g_{3}$ $\displaystyle=\,0.531+6.958\cdot 10^{-7}u_{2}-2.51\cdot
10^{-12}u_{2}^{2}$ $\displaystyle\,+3.1\cdot
10^{-18}u_{2}^{3}+\frac{0.43-0.14(pr_{x_{2}}-pr_{x_{2,crit}})}{1+7.35(pr_{x_{2}}-pr_{x_{2,crit}})}$
(13) $\displaystyle g_{4}$ $\displaystyle=\,\frac{0.5517+3.8168\cdot
10^{-6}u_{2}}{1+8.271\cdot 10^{-6}u_{2}}$ (14) $\displaystyle g_{5}$
$\displaystyle=\,\frac{3.8168\cdot 10^{-6}\cdot g_{3}\cdot g_{4}\cdot
u_{2}}{g_{2}(1-g_{3})}$ (15)

$g_{1}$ is the liquidus temperature $T_{liq}$, $g_{2}$ is the electrical
conductivity $\kappa$ from Eq. (71), $g_{3}$ is the bubble coverage $\phi$
from Eq. 75, $g_{4}$ is the bubble thickness $d_{bub}$ from Eq. (74) and
$g_{5}$ is the bubble voltage drop $U_{bub}$. Notice that:

$\displaystyle c_{x_{2}}$ $\displaystyle=\frac{x_{2}}{x_{2}+x_{3}+x_{4}}$
$\displaystyle c_{x_{3}}$ $\displaystyle=\frac{x_{3}}{x_{2}+x_{3}+x_{4}}$
$\displaystyle pr_{x_{2}}$ $\displaystyle=100\cdot c_{x_{2}}$ $\displaystyle
pr_{x_{3}}$ $\displaystyle=100\cdot c_{x_{3}}$

With the nonlinear functions $g_{1},\;...\;g_{5}$ described, the state space
equations in Eq. 10 is described by the following:

$\displaystyle\dot{x}_{1}$
$\displaystyle=\,\frac{k_{1}(g_{1}-x_{7})}{x_{1}k_{0}}-k_{2}(x_{6}-g_{1})$
(17) $\displaystyle\dot{x}_{2}$ $\displaystyle=\,u_{1}-k_{3}u_{2}$ (18)
$\displaystyle\dot{x}_{3}$ $\displaystyle=\,u_{3}-k_{4}u_{1}$ (19)
$\displaystyle\dot{x}_{4}$
$\displaystyle=\,-\left(\frac{k_{1}(g_{1}-x_{7})}{x_{1}k_{0}}-k_{2}(x_{6}-g_{1})\right)+k_{5}u_{1}$
(20) $\displaystyle\dot{x}_{5}$ $\displaystyle=\,k_{6}u_{2}-u_{4}$ (21)
$\displaystyle\dot{x}_{6}$
$\displaystyle=\,\frac{\alpha}{x_{2}+x_{3}+x_{4}}\biggl{[}u_{2}\left(g_{5}+\frac{u_{2}u_{5}}{2620g_{2}}\right)$
$\displaystyle\,\;-k_{9}\frac{x_{6}-x_{7}}{k_{10}+k_{11}k_{0}x_{1}}$
$\displaystyle\,\;-\left(k_{7}(x_{6}-g_{1})^{2}-k_{8}\frac{(x_{6}-g_{1})(g_{1}-x_{7})}{k_{0}x_{1}}\right)\biggr{]}$
(22) $\displaystyle\dot{x}_{7}$
$\displaystyle=\,\frac{\beta}{x_{1}}\biggl{[}-\biggl{(}\quad
k_{12}(x_{6}-g_{1})(g_{1}-x_{7})$
$\displaystyle\,\;-k_{13}\frac{(g_{1}-x_{7})^{2}}{k_{0}x_{1}}\biggr{)}$
$\displaystyle\,+\frac{k_{9}(g_{1}-x_{7})}{k_{15}k_{0}x_{1}}-\frac{k_{9}(x_{7}-x_{8})}{k_{14}+k_{15}k_{0}x_{1}}\biggr{]}$
(23) $\displaystyle\dot{x}_{8}$
$\displaystyle=\,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)$ (24)

##  3 Method and experimental setup

###  3.1 Sparsity promoting regularization

In this article, sparse DNN models are utilized to predict state variables in
an aluminum electrolysis simulator. All the weight matrices in a DNN model are
$\ell_{1}$ regularized to impose sparsity on the model. Fig. 2 illustrates how
weights are enumerated according to their input and output nodes. Layer $j$
has $i$ nodes and layer $(j+1)$ has $r$ nodes. Layer $j=0$ corresponds to the
input layer, and consist of measured or estimated states $\mathbf{x}(t)$ and
control inputs $\mathbf{u}(t)$ at time step $t$. The output layer consist of
the estimated time derivatives of the states $\mathbf{\dot{x}}(t)$ at time
step $t$.

Figure 2: Enumerated weights according to their input and output nodes between
layer $j$ and $(j+1)$

The weight matrix $\mathbf{W}_{j+1}$ maps the weights from layer $j$ to layer
$j+1$. $\mathbf{W}_{j+1}$ is arranged as follows:

$\mathbf{W}_{j+1}=\begin{bmatrix}w_{11}&w_{12}&\;...\;&w_{1i}\\\
w_{21}&w_{22}&\;...\;&w_{2i}\\\ \vdots&\ddots&\\\
w_{r1}&w_{r2}&\ldots&w_{ri}\end{bmatrix}.$ (25)

Regularization terms $R_{\ell_{1},j+1}$ are defined for each weight matrix
$\mathbf{W}_{j+1}$ :

$R_{\ell_{1},j+1}=\sum_{i,k\in W_{j+1}}|w_{i,k}|.$ (26)

where $w_{i,k}$ are the model weights from layer $j$ to $j+1$, or
equivalently, the elements of $\mathbf{W}_{j+1}$. The regularization terms for
each layer are added to the cost function:

$\displaystyle\mathbf{w}^{*}$ $\displaystyle=\
\underset{\mathbf{w}}{\mathrm{argmin}}\biggl{\\{}\frac{1}{N}\sum_{i=1}^{N}(\mathbf{y_{i}}-\mathbf{f(x_{i}}))^{2}+\lambda_{1}R_{\ell_{1},1}...$
(27) $\displaystyle\
\quad+\lambda_{j+1}R_{\ell_{1},j+1}+...+\lambda_{L}R_{\ell_{1},L}\biggr{\\}},$

where $\frac{1}{N}\sum_{i=1}^{N}(\mathbf{y_{i}}-\mathbf{f(x_{i}}))^{2}$ is the
MSE, and $\lambda_{j+1},j=0,\;...,3$ is a layer-specific hyperparameter that
determines how the weights in $\mathbf{W}_{j+1}$ are penalized.

###  3.2 Experimental setup and data generation

Table 2: Constants in the simulator Constant | Physical meaning | Nummeric value
---|---|---
$k_{0}$ | $1/(\rho_{sl}A_{sl})$ | $2\cdot 10^{-5}$
$k_{1}$ | $2k_{sl}A_{sl}/\Delta_{fus}H_{cry}$ | $7.5\cdot 10^{-4}$
$k_{2}$ | $h_{bath-sl}A_{sl}/\Delta_{fus}H_{cry}$ | $0.18$
$k_{3}$ | $0.002\frac{M_{Al_{2}O_{3}}\cdot CE}{z\cdot F}$ | $1.7\cdot 10^{-7}$
$k_{4}$ | $C_{Na_{2}O}\frac{4M_{AlF_{3}}}{3M_{Na_{2}O}}$ | $0.036$
$k_{5}$ | $C_{Na_{2}O}\frac{2M_{cry}}{3M_{Na_{2}O}}$ | $0.03$
$k_{6}$ | $0.002\frac{M_{Al}\cdot CE}{z\cdot F}$ | $4.43\cdot 10^{-8}$
$k_{7}$ | $k_{2}\cdot c_{p_{cry,\;liq}}$ | $338$
$k_{8}$ | $k_{1}\cdot c_{p_{cry,\;liq}}$ | $1.41$
$k_{9}$ | $A_{sl}$ | $17.92$
$k_{10}$ | $1/h_{bath-sl}$ | $0.00083$
$k_{11}$ | $1/(2k_{sl})$ | $0.2$
$k_{12}$ | $k_{2}\cdot c_{p_{cry,\;s}}$ | $237.5$
$k_{13}$ | $k_{1}\cdot c_{p_{cry,\;s}}$ | $0.99$
$k_{14}$ | $x_{wall}/(2k_{wall})$ | $0.0077$
$k_{15}$ | $1/(2k_{sl})$ | $0.2$
$k_{16}$ | $T_{0}$ | $35$
$k_{17}$ | $1/(m_{wall}c_{p,\;wall})$ | $5.8\cdot 10^{-7}$
$k_{18}$ | $1/h_{wall-0}$ | $0.04$
$\alpha$ | $1/c_{p_{bath,\;liq}}$ | $5.66\cdot 10^{-4}$
$\beta$ | $1/c_{p_{cry,\;sol}}$ | $7.58\cdot 10^{-4}$

Data for the aluminium electrolysis process is generated by integrating the
non-linear ODEs given by Eqs. (17) - (24) with a set of chosen initial values
for the state variables $\mathbf{x}(t_{0})$, and fourth-order Runge-Kutta
(RK4) algorithm. The initial conditions of each variable $x_{i}$ for each
simulation are randomly chosen from a given interval of possible initial
conditions:

Table 3: Initial conditions for system variables Variable | Initial condition interval
---|---
$x_{1}$ | $[3260,\;3260]$
$c_{x_{2}}$ | $[0.02,\;0.03]$
$c_{x_{3}}$ | $[0.10,\;0.12]$
$x_{4}$ | $[13500,\;14000]$
$x_{5}$ | $[9950,\;10000]$
$x_{6}$ | $[975,\;975]$
$x_{7}$ | $[816,\;816]$
$x_{8}$ | $[580,\;580]$

Table 3 shows the intervals for initial conditions of each variable. For
$x_{2}$ and $x_{3}$, concentrations $c_{x_{2}}$ and $c_{x_{3}}$ are given.

Data-driven models depend on a high degree of variation in the training data
to be reliable and valid in a large area of the input space. The input signal
determines how the system is excited and thus what data is available for
modeling and parameter estimation. Operational data from a controlled, stable
process is generally characterized by a low degree of variation. Even large
amounts of data sampled over a long period from a controlled process can not
guarantee that the variation in the training data i large enough to ensure
that the trained model generalizes to unseen data. Random excitations are
added to the input signals to increase the variation in the sampled data. The
intuition is that the random excitation will push the dynamics out of the
standard operating condition so that variation in the training data increases.
In general, each control input $i$ is given by:

$u_{i}=\textrm{Deterministic term + Random term}.$ (28)

The control inputs $u_{1},\;u_{3}$ and $u_{4}$ are impulses. The random term
is zero for these control inputs when the deterministic term is zero. The
deterministic term is a proportional controller. The control inputs $u_{2}$
and $u_{5}$ are always nonzero. These control inputs have constant
deterministic terms and a random term that changes periodically. The random
term stays constant for a certain period $\Delta T_{rand}$ before changing to
a new randomly determined constant. Choosing the period $\Delta T_{rand}$ is a
matter of balancing different objectives. On the one hand, it is desirable to
choose a large period $\Delta T_{rand}$ so that the system can stabilize and
evolve under the given conditions to reveal the system dynamics under the
given conditions. On the other hand, it is desirable to test the systems under
many different operational conditions. By empirically testing different
periods $\Delta T_{rand}$, and seeing how the dynamics evolve in simulation,
it turns out that setting $\Delta T_{rand}=30\Delta T$ is a fair compromise
between the two.

Table 4: Control functions Input | Deterministic term | Random term interval | $\Delta T_{rand}$
---|---|---|---
$u_{1}$ | $3e4(0.023-c_{x_{2}})$ | $[-2.0,2.0]$ | $\Delta T$
$u_{2}$ | $14e3$ | $[-7e3,7e3]$ | $30\cdot\Delta T$
$u_{3}$ | $13e3(0.105-c_{x_{3}})$ | $[-0.5,\;0.5]$ | $\Delta T$
$u_{4}$ | $2(x_{5}-10e3)$ | $[-2.0,\;2.0]$ | $\Delta T$
$u_{5}$ | $0.05$ | $[-0.015,\;0.015]$ | $30\cdot\Delta T$

Table 4 gives the numerical values of the deterministic term of the control
input, the interval of values for the random terms, and the duration $\Delta
T_{rand}$ of how long the random term is constant before either becoming zero
($u_{1},u_{3},u_{4}$) or changing to a new randomly chosen value
($u_{2},u_{5}$). One simulation $i$ with a given set of initial conditions is
simulated for $1000$ time steps, and each time step $\Delta T=30s$. The
simulation generate the data matrix as in Eq. (29):

$\mathbf{X}=\leavevmode\resizebox{160.43727pt}{}{
$\begin{bmatrix}x_{1}(0)&x_{2}(0)&\;...\;&x_{8}(0)&u_{1}(0)&\;...\;&u_{5}(0)\\\
\\\ x_{1}(1)&x_{2}(1)&\;...\;&x_{8}(1)&u_{1}(1)&\;...\;&u_{5}(1)\\\
\vdots&\vdots&\ddots&\vdots&\vdots&\ddots&\vdots\\\
x_{1}(k)&\ldots&x_{i}(k)&\ldots&u_{1}(k)&\ldots&u_{5}(k)\\\
\vdots&\vdots&\ddots&\vdots&\vdots&\ddots&\vdots\\\
x_{1}(999)&x_{2}(999)&\;...\;&x_{8}(999)&u_{1}(999)&\;...\;&u_{5}(999)\\\
\end{bmatrix}.$ }$ (29)

The number $k$ within the parenthesis of variable $i$ indicate the time step
for when $x_{i}(k)$ is sampled. The target values are then calculated as

$\mathbf{Y}=\begin{bmatrix}\frac{x_{1}(1)-x_{1}(0)}{\Delta
T}&\;...\;&\frac{x_{8}(1)-x_{8}(0)}{\Delta T}\\\ \\\
\frac{x_{1}(2)-x_{1}(1)}{\Delta T}&\;...\;&\frac{x_{8}(2)-x_{8}(1)}{\Delta
T}\\\ \vdots&\ddots&\vdots\\\ \frac{x_{1}(k)-x_{1}(k-1)}{\Delta
T}&\;...\;&\frac{x_{8}(k)-x_{8}(k-1)}{\Delta T}\\\ \vdots&\ddots&\vdots\\\
\frac{x_{1}(999)-x_{1}(998)}{\Delta
T}&\;...\;&\frac{x_{8}(999)-x_{8}(998)}{\Delta T}\\\ \end{bmatrix}.$

(30)

The training set consist of smaller training sets $\mathcal{S}_{k}$:

$\mathcal{S}_{train}=\left\\{\mathcal{S}_{1},\;\mathcal{S}_{2},\;...,\;\mathcal{S}_{k},\;....,\;\mathcal{S}_{n}\right].$
(31)

Here, $n$ is the number of simulated time-series $\mathbf{X}$. Each training
set $\mathcal{S}_{k}$ from simulation $k$ are put in pairs for regression:

$\mathcal{S}_{k}=[\mathbf{X},\mathbf{Y}]=\begin{bmatrix}[\mathbf{x}(1),\mathbf{u}(1)],&\mathbf{y}(1)\\\
\vdots&\vdots\\\ [\mathbf{x}(k),\mathbf{u}(k)],&\mathbf{y}(k)\\\
\vdots&\vdots\\\ [\mathbf{x}(999),\mathbf{u}(999)],&\mathbf{y}(999)\\\
\end{bmatrix}.$ (32)

The number of time series simulations $n$ varies in the experiments to
evaluate the model performance as a function of training size. The test set
also consists of several time series simulations generated in the same way
described above. The test set is given by:

$\mathcal{S}_{test}=\\{\\{X_{1}\\},\;\\{X_{2}\\},\;...,\;\\{X_{p}\\}\\}.$ (33)

$\mathcal{S}_{test}(i)=\\{X_{i}\\}=\\{([\mathbf{x}_{k},\;\mathbf{u}_{k}]\\}_{k=1}^{1000}$
is one simulated time series, and $\\{\mathbf{x}_{k}\\}_{k=1}^{1000}$ is being
forcasted by the models. $p=20$ is the number of simulated time series in the
test set. In all experiments 20 models of each dense and sparse networks with
different initialization are trained on a training set and evaluated on a test
set.

###  3.3 Performance metrics

In aluminum electrolysis, data is generally sampled at rare instants during
the operation. Thus, a model needs to accurately forecast several time steps
without feedback from measurements to ensure safe and optimal operation.
Therefore, the models’ capability to estimate the states $\mathbf{x}$ over a
given time horizon without measurement feedback becomes an important measure
of performance. The initial conditions $\mathbf{x}(t_{0})$ are given to the
models. Then the consecutive $n$ time steps of the states are estimated
$\\{\mathbf{\hat{x}}(t_{1}),\;...,\;\mathbf{\hat{x}}(t_{n})\\}$. This is
called a rolling forecast. The model estimates the time derivatives of the
states $d\mathbf{\hat{x}_{i}}/dt$ based on the current control inputs
$\mathbf{u}(t_{i})$ and initial conditions $\mathbf{x}_{0}(t)$ if $t=t_{0}$,
or the estimate of the current state variables $\mathbf{\hat{x}}(t_{i})$ if
$t>t_{0}$:

$\frac{d\mathbf{\hat{x}}(t_{i})}{dt}=\begin{cases}\hat{f}(\mathbf{\hat{x}}(t_{i}),\;\mathbf{u}(t_{i})),&\text{if
}t_{i}>t_{0}\\\ \hat{f}(\mathbf{x}_{0}(t_{i}),\;\mathbf{u}(t_{i})),&\text{if
}t_{i}=t_{0}\end{cases}$ (34)

Then, the next state estimate $\mathbf{x}(t_{i+1})$ is calculated as

$\mathbf{\hat{x}}(t_{i+1})=\mathbf{\hat{x}}(t_{i})+\frac{d\mathbf{\hat{x}}(t_{i})}{dt}\cdot\Delta
T.$ (35)

The rolling forecast can be computed for each of the states $x_{i}$ for one
set of test trajectories $\mathcal{S}_{test}$. However, presenting the rolling
forecast of multiple test sets would render the interpretation difficult. By
introducing a measure called Average Normalized Rolling Forecast Mean Squared
Error (AN-RFMSE) that compresses the information about model performance, the
models can easily be evaluated on a large number of test sets. The AN-RFMSE is
a scalar defined as:

$\textrm{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})}{std(x_{i})}\right)^{2},$
(36)

where $\hat{x}_{i}(t_{j})$ is the model estimate of the simulated state
variable $x_{i}$ at time step $t_{j}$, $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 the normalized
rolling forecast MSE is averaged over. Hence, for every model $\hat{f}_{j}$
and every test set time series $\mathcal{S}_{test}(i)$, there is a
corresponding AN-RFMSE. This generates a matrix, where each row represents
individual model instances, and every column represents one test set
simulation $\mathcal{S}_{test}(i)$. Each entry in the matrix is the AN-RFMSE
for a given model instance and a given $\mathcal{S}_{test}(i)$. The matrix is
given by:

$\leavevmode\resizebox{173.44534pt}{}{ $\textrm{AN-
RFMSE}_{mat}=\begin{bmatrix}\textrm{AN-RFMSE}_{11}&\ldots&\textrm{AN-
RFMSE}_{1n}\\\ \vdots&\vdots&\vdots\\\ \textrm{AN-
RFMSE}_{k1}&\ldots&\textrm{AN-RFMSE}_{kn}\end{bmatrix}$}.$ (37)

$k$ is the number of models, and $n$ is the number of time series simulations
in the test set $\mathcal{S}_{test}$. There are two $\textrm{AN-RFMSE}_{mat}$
matrices, one for sparse models and one for dense models. Averaging over all
columns at each row, that is, averaging over all test set time series for each
model instance, generates a vector

$\overline{\textrm{AN-
RFMSE}}_{vec}=\begin{bmatrix}\frac{1}{n}\sum_{i_{i}}^{n}\textrm{AN-
RFMSE}_{1i}\\\ \vdots\\\ \frac{1}{n}\sum_{i_{i}}^{n}\textrm{AN-
RFMSE}_{ki}\end{bmatrix}=\begin{bmatrix}\overline{\textrm{AN-RFMSE}_{1}}\\\
\vdots\\\ \overline{\textrm{AN-RFMSE}_{k}}\end{bmatrix}.$

(38)

The elements of $\overline{\textrm{AN-RFMSE}}_{vec}$ is the average AN-RFMSE
over all test set time series for every model instance.

##  4 Results and discussion

In this section, we present and discuss the main findings of the work. In
doing so, we will analyze the results from the perspective of
interpretability, generalizability and training stability.

###  4.1 Interpretability perspective

As discussed earlier, the interpretability of a model is the key to its
acceptability in high-stake applications like the aluminum extraction process
considered here. Unfortunately, highly complex dense neural networks having
thousands to millions of parameters which were the starting point for the
modeling here are almost impossible to interpret. Refer to the Fig. 3 that
shows the model structure of a dense DNN model learned for the generated data.

Figure 3: Model structure of each output function $\\{f_{1},..,f_{8}\\}$ for
one of the trained dense neural network models. Blue circles represent
neurons, inputs and outputs in the model. $X_{i}$ represent system variable
$i$ in the input layer, $U_{i}$ represents control input $i$ in the input
layer, and $Z_{i}$ represents latent variable $i$ in the layer it is
visualized. The directed edges indicate weights in the model.

Fortunately, through the regularization it was possible to significantly
reduce the model complexity resulting in a drastically reduced number of
trainable parameters (see the Figs. 4-11). It can be argued that the reduced
model complexity of sparse neural networks increases the model
interpretability. With domain knowledge about the aluminum electrolysis
process, the sparse models can be evaluated as we will do in the remainder of
this section.

The results related to the interpretability aspect is presented in the form of
model structure plots which can be used to explain the input-output mapping of
the models. If the model structures are very sensitive to the initialization,
then there interpretation will not make sense therefore, 100 DNNs with
different initialization are trained independently and their common trends are
emphasized in the discussions. We now present each of the model outputs
$\\{\hat{f}_{1}(\mathbf{x},\mathbf{u}),\;...,\;\hat{f}_{8}(\mathbf{x},\mathbf{u})\\}$.
It is worth mentioning that these outputs are estimates of each of the time
derivatives of the states $\\{\dot{x}_{1},\;...,\;\dot{x}_{8}\\}$
respectively.

Table 5: Frequency of learned features for each output $\\{\hat{f}_{1},\;...,\;\hat{f}_{8}\\}$. Each column $i$ correspond to an output $\hat{f}_{i}$ of the neural network. Each row element $j$ correspond to one of the features $\\{x_{1},\;x_{2},\;...,\;x_{8},\;u_{1},\;u_{2},\;...,\;u_{5}\\}$. The value of the table element $(i,j)$ is the percent of how many out of one hundred models of the output $f_{i}$ that feature $j$ occurs.  Feature | Output functions
---|---
$\hat{f}_{1}$ | $\hat{f}_{2}$ | $\hat{f}_{3}$ | $\hat{f}_{4}$ | $\hat{f}_{5}$ | $\hat{f}_{6}$ | $\hat{f}_{7}$ | $\hat{f}_{8}$
$x_{1}$ | $86$ | $1$ | $1$ | $86$ | $2$ | $100$ | $99$ | $100$
$x_{2}$ | $100$ | $2$ | $2$ | $100$ | $1$ | $100$ | $100$ | $100$
$x_{3}$ | $100$ | $2$ | $2$ | $100$ | $0$ | $93$ | $100$ | $87$
$x_{4}$ | $100$ | $0$ | $0$ | $100$ | $0$ | $87$ | $100$ | $87$
$x_{5}$ | $2$ | $0$ | $0$ | $2$ | $0$ | $2$ | $2$ | $2$
$x_{6}$ | $100$ | $2$ | $2$ | $100$ | $1$ | $100$ | $100$ | $87$
$x_{7}$ | $22$ | $1$ | $1$ | $21$ | $1$ | $7$ | $89$ | $89$
$x_{8}$ | $100$ | $1$ | $1$ | $100$ | $2$ | $87$ | $100$ | $100$
$u_{1}$ | $100$ | $100$ | $100$ | $100$ | $0$ | $100$ | $100$ | $87$
$u_{2}$ | $4$ | $0$ | $0$ | $4$ | $1$ | $100$ | $18$ | $18$
$u_{3}$ | $2$ | $0$ | $100$ | $2$ | $0$ | $2$ | $4$ | $4$
$u_{4}$ | $2$ | $0$ | $0$ | $2$ | $100$ | $3$ | $3$ | $3$
$u_{5}$ | $3$ | $0$ | $0$ | $3$ | $1$ | $100$ | $18$ | $18$

####  4.1.1 Model output $\hat{f}_{1}$

The simulation model for the first output $f_{1}$ defined in Eq. (17) is a
function of the features $\\{x_{1},x_{2},x_{3},x_{4},x_{6},x_{7}\\}$. $f_{1}$
can further be divided into three sums
$f_{1}=h_{1}(x_{1},x_{2},x_{3},x_{4},x_{7})+h_{2}(x_{6})+h_{3}(x_{2},x_{3},x_{4})$,
where $h_{1}=k_{1}\frac{g_{1}(x_{2},x_{3},x_{4})-x_{7}}{k_{0}x_{1}}$,
$h_{2}=-k_{2}x_{6}$ and $h_{3}=g_{1}(x_{2},x_{3},x_{4})$. $g_{1}$ is a
nonlinear function defined in Eq. (11).

Figure 4: Most common learned structures for
$\hat{f}_{1}(\mathbf{x},\mathbf{u})$.

Fig. 4 shows the three most common learned structures of the first output of
the neural network $\hat{f}_{1}(\mathbf{x},\mathbf{u})$. In total, these
structures account for $86\%$ of all learned structures of $\hat{f}_{1}$. The
top structure forms the resulting structure for $63\%$ of the models. It is a
function of seven inputs
$\hat{f}_{1}=f(x_{1},x_{2},x_{3},x_{4},x_{6},x_{8},u_{1})$. All input features
are connected to the same neuron in the first layer. Moreover, it is only one
neuron in each hidden layer. The upper bound in Eq. (9) states that this model
structure only has $n^{dL}=1^{7\cdot 3}=1$ region. This is equivalent to
stating that the model collapses to one linear model. The middle and the
bottom structures of Fig. 4 has more than one neuron in the hidden layers.
However, the structures can be divided into two disconnected subnetworks since
the hidden neurons are not connected before they are added in the output
layer. Hence, also these models collapse to linear models with a single
region. This means that the neural networks do not capture the nonlinear
dynamics in the simulation model. All structures of Fig. 4 are erroneously
including $x_{8}$ and $u_{1}$ in their feature basis. Besides, $x_{7}$ which
is present in the simulation model $f_{1}$ is not found by the top and bottom
structures in Fig. 4. In fact, $x_{7}$ is found only in $22\%$ of the models
$\hat{f}_{1}$ according to Table 5. The exact cause of this erroneous feature
selection is not trivial. However, $x_{8}$ which is the wall temperature
correlates highly with the side ledge temperature $x_{7}$. Thus, $x_{8}$ can
possible have been learned as a feature of $f_{1}$ instead of $x_{7}$.
Moreover, the alumina feed $u_{1}$ on the other hand affects the time
derivative of the mass of alumina $\dot{x}_{2}$, the time derivative of mass
of aluminum fluoride $\dot{x}_{3}$ and the time derivative of cryolite
$\dot{x}_{4}$ directly. All these variables affect $\dot{x}_{1}$ through the
liquidus temperature $g_{1}(x_{2},x_{3},x_{4})$. To understand how this might
cause the learning algorithm to find $u_{1}$ as a feature of $\hat{f}_{1}$,
consider the following: let $u_{1}$ be zero until time $t$. Then,
$\\{x_{2},x_{3},x_{4}\\}$ will be updated due to $u_{1}$ at the next sampled
time step $t+1$. However, the fourth-order Runge Kutta solver splits the
sampling interval into 4 smaller intervals
$\\{t+0.25,\;t+0.5,\;t+0.75,\;t+1\\}$ and solve the ODE eqautions at all these
time steps. Thus, the state variables $\\{x_{2},x_{3},x_{4}\\}$ are updated
already at time $t+0.25$. Since $\dot{x}_{1}$ is depending on these variables,
$x_{1}$ will be updated at $t+0.5$. Therefore, at time $t+1$, when data is
sampled, $x_{1}$ would also be changed. Hence, the learning algorithm finds
$u_{1}$ to affect the time derivative $\dot{x}_{1}$. This could might have
been solved by shortening the sampling interval. Furthermore, $x_{1}$ not
included as a feature in $14\%$ of the models. This might be a combination of
parameter initialization and that $x_{1}$ is multiplied by the small constant
$k_{0}=2\cdot 10^{-5}$.

Figure 5: Most common learned structure for
$\hat{f}_{2}(\mathbf{x},\mathbf{u})$. $97\%$ of $f_{2}$ ends up with this
structure.

####  4.1.2 Model output $\hat{f}_{2}$

Fig. 5 shows the most common learned structure among the models $\hat{f}_{2}$
that models the time derivative of alumina $f_{2}=\dot{x}_{2}$. $97\%$ of the
structures end up as the structure in Fig. 5. The simulation model $f_{2}$ in
Eq. (18) is a linear model dependent on $\\{u_{1},u_{2}\\}$. The learned
models $\hat{f}_{2}$ only finds $u_{1}$ as the relevant feature. The reason
for this might be that $u_{2}$ is proportional to a very small constant
$k_{3}=1.7\cdot 10^{-7}$. Variations in the line current $u_{2}$ are not big
enough to be of significance for the learning algorithm. The dynamics caused
by $u_{2}$ is instead captured in a bias in the models.

Figure 6: Most common learned structures for
$\hat{f}_{3}(\mathbf{x},\mathbf{u})$.

####  4.1.3 Model output $\hat{f}_{3}$

Fig. 6 shows the four most common learned structures of $\hat{f}_{3}$.
$\hat{f}_{3}$ models the time derivative of the aluminum fluoride mass
$\dot{x}_{3}$ in the cell. The simulation model $\dot{x}_{3}=f_{3}$ in Eq.
(19) is a linear model depending on the features $\\{u_{1},u_{3}\\}$, and in
all structures in Fig. 6, only these two features are found. As shown in Table
5, these features are found in $100\%$ of the trained models. The structures
found are mainly linear models, however in the second and forth structures,
there are some weights that connects the features in intermediate layers.

Figure 7: Most common learned structures for
$\hat{f}_{4}(\mathbf{x},\mathbf{u})$.

####  4.1.4 Model output $\hat{f}_{4}$

Fig. 7 show the four most common learned structures among models $\hat{f}_{4}$
that models the mass rate of liquid cryolite $Na_{3}AlF_{6}$ in the bath,
namely $\dot{x}_{4}$. $\dot{x}_{4}$ is simulated by the simulation model
$\dot{x}_{4}=f_{4}$ in Eq. (19). $f_{4}$ consist of the features
$\\{x_{1},x_{2},x_{3},x_{4},x_{6},x_{7},u_{1}\\}$. Table 5 show that
$\\{x_{2},x_{3},x_{4},x_{6},x_{8},u_{1}\\}$ are included in $100\%$ of the
learned models, $x_{1}$ is included in $86\%$ of the models and $x_{7}$ is
included in only $21\%$ of the models. As for $\hat{f}_{1}$, $x_{8}$ is
erroneously included in the basis of the model. This might be explained by
that $x_{7}$ and $x_{8}$ highly correlates, and that the wrong feature is
included. The structures in Fig. 7 are all forming linear models. However, the
simulation model $\dot{x}_{4}=f_{4}$ is partly nonlinear. Thus, the
approximation $\hat{f}_{4}$ oversimplifies the dynamics. This might be caused
by a high weighting of the sparse regularisation term. If the loss function is
less penaltized, it is room for more weights in the model and therefore also
more nonlinearities.

####  4.1.5 Model output $\hat{f}_{5}$

Figure 8: Most common learned structures for
$\hat{f}_{5}(\mathbf{x},\mathbf{u})$.

Fig. 8 show the most common learned structure among the models $\hat{f}_{5}$
and include $98\%$ of the learned model structures. $\hat{f}_{5}$ is modeling
the mass rate of produced aluminum in the cell $\dot{x}_{5}$. The $x_{5}$ time
series are produced by the simulation model $\dot{x}_{5}=f_{5}$ in Eq. 21.
$f_{5}$ is a linear model dependent on the features $\\{u_{2},u_{4}\\}$.
However, most of the model structures are only depending on $u_{4}$. This can
be caused by the fact that $u_{2}$ is proportional to a very small constant
$k_{6}=4.43\cdot 10^{-8}$. Thus, variations in $u_{2}$ might not be large
enough for the learning algorithm to find $u_{2}$ significant as a basis for
$\hat{f}_{5}$.

####  4.1.6 Model output $\hat{f}_{6}$

Figure 9: Most common learned structures for
$\hat{f}_{6}(\mathbf{x},\mathbf{u})$.

Fig. 9 show the most common model structures of $\hat{f}_{6}$. $\hat{f}_{6}$
models the bath temperature time derivative $\dot{x}_{6}$. The bath
temperature $x_{6}$ is simulated by the ODE in Eq. (22). It is a nonlinear
equation depending on $\\{x_{1},x_{2},x_{3},x_{4},x_{6},x_{7},u_{2},u_{5}\\}$.
The most common structure, learned by $57\%$ of the models $\hat{f}_{6}$ is
illustrated in the top plot of Fig. 9. This structure has the basis
$\\{x_{1},x_{2},x_{3},x_{4},x_{6},x_{8},u_{1},u_{2},u_{5}\\}$. Hence, it finds
$u_{1}$ and $x_{8}$, which is erroneously found in many of the structures
above. A possible explanation for this trend is given above. The structure has
two neurons in the first layer, one with the basis
$\\{x_{1},x_{2},x_{6},x_{8},u_{1},u_{2},u_{5}\\}$ and one with the basis
$\\{x_{1},x_{2},x_{3},x_{4},x_{6},x_{8},u_{1}\\}$. Since the model collapses
to a linear model, all terms are summed in the end. Thus, the separation of
the basis is thus of little importance. The second plot from above in Fig. 9
is the second most common structure, and accounts for $11\%$ of the models
$\hat{f}_{6}$. It has the same feature basis as the most common structure, but
they are arranged differently in the first layer. However, since both model
structures are linear models, this arrangement is of minor importance. The
third structure in Fig. 9 which account for $7\%$ of the structures of
$\hat{f}_{6}$ has the basis $\\{x_{1},x_{2},x_{6},u_{1},u_{2},u_{5}\\}$.
Compared to the other structures, $\\{x_{3},x_{4},x_{8}\\}$ is not present.
Since it only happens rarely, this is maybe partly caused by bad parameter
initialization. The fourth most common structure, which accounts for $6\%$ of
the structures, has the same feature basis as the first and the second most
common structures. This structure is plotted in the bottom of Fig. 9. While
all other structures in Fig. 9 have one linear response region, the fourth
most common model structure models some nonlinearities. That is, neurons in
the first hidden layer are connected in the second hidden layer. Hence, the
input space must be divided into several linear response regions for this
structure.

####  4.1.7 Model output $\hat{f}_{7}$

Figure 10: Most common learned structures for
$\hat{f}_{7}(\mathbf{x},\mathbf{u})$.

Fig. 10 show the two most common structures for the models $\hat{f}_{7}$.
$\hat{f}_{7}$ models the time derivative of the side ledge temperature
$\dot{x}_{7}$. $\dot{x}_{7}=f_{7}$ is simulated by the ODE in Eq. (23).
$\dot{x}_{7}$ is depending on the feature basis
$\\{x_{1},x_{2},x_{3},x_{4},x_{6},x_{7},x_{8}\\}$. Table 5 states that the
basis $\\{x_{2},x_{3},x_{4},x_{6},x_{8},u_{1}\\}$ is present for $100\%$ of
the models, $x_{1}$ is present for $99\%$ of the models and $x_{7}$ is present
in $89\%$ of the models. The top plot in Fig. 10 show the structure that
accounts for $63\%$ of the models. The bottom plot account for $10\%$ of the
model structures of $\hat{f}_{7}$. These two structures collapse to linear
models, and have the same feature basis
$\\{x_{1},x_{2},x_{3},x_{4},x_{6},x_{7},x_{8},u_{1}\\}$, but have minor
differences in how weights are connected between input layer and first hidden
layer. $u_{1}$ is also for this model output erroneously found as a basis, and
a possible explanation is mentioned above.

####  4.1.8 Model output $\hat{f}_{8}$

Figure 11: Most common learned structures for
$\hat{f}_{8}(\mathbf{x},\mathbf{u})$.

Fig. 11 show the two most common model structures for the last model output
$\hat{f}_{8}$. $\hat{f}_{8}$ models the time derivative of the wall
temperature $\dot{x}_{8}=f_{8}$, which is simulated in Eq. (24). $\dot{x}_{8}$
depends on the feature basis $\\{x_{1},x_{7},x_{8}\\}$. However, the most
common learned structure for $\hat{f}_{8}$ has the basis
$\\{x_{1},x_{2},x_{3},x_{4},x_{6},x_{7},x_{8},u_{1}\\}$. This is the exact
same structure as the most common learned structure for $\hat{f}_{7}$.
Therefore, a possible explanation is that $\hat{f}_{8}$ adapts the same
parameters as $\hat{f}_{7}$ in some cases as they highly correlates. The
bottom plot in Fig. 11 show the second most common learned structure of the
model output $\hat{f}_{8}$. This structure is learned by $11\%$ of the models
$\hat{f}_{8}$. The feature basis for this structure is
$\\{x_{1},x_{2},x_{7},x_{8}\\}$, and reminds more of the actual basis. In this
structure, there is only one erroneous learned feature, namely $x_{2}$. Fig.
4-11 and Table 5 show that the sparse learning is quite consistent in finding
the same feature basis and structure with similar characteristics. However,
some differences that affect the models are present.

It is clear that doing a similar analysis for models in Fig. 3 is impossible
as all interconnections make the models a black box. On average, $93\%$ of the
weights in the sparse DNN models are pruned. For the outputs
$\hat{f}_{1},\;\hat{f}_{4},\;\hat{f}_{6},\;\hat{f}_{7}$ and $\hat{f}_{8}$,
approximately $40\%$ of the input features are pruned of the model structures.
For $\hat{f}_{2},\;\hat{f}_{3}$ and $\hat{f}_{5}$, $85-95\%$ of the input
features are pruned. For all outputs of the sparse DNN models, around
$85-95\%$ of the neurons are pruned at each layer. In a neural network, the
number of matrix operations only decreases if neurons are pruned. That is,
removing a neuron in layer $j$ is equivalent to removing a row in weight
matrix $\mathbf{W}_{j}$ and a column in weight matrix $\mathbf{W}_{j+1}$. The
dense models have a compact model structure, where most of the weights are
nonzero. The dense DNN models in the case study have the shapes 13-15-14-12-8.
The first number is the number of features, the second, third, and fourth
numbers are the numbers of neurons in hidden layers, and the last number is
the number of outputs. This shape gives $669$ matrix operations in a forward
pass. An average sparse DNN, has the shape $13-6-6-6-8$. This gives $198$
matrix operations. Thus, the number of matrix operations in the forward pass
of a sparse DNN model is reduced by approximately $70\%$.

###  4.2 Generalizability perspective

This section focuses on the models’ performance on test data in terms of
accuracy and uncertainty. Furthermore, we also investigate the impact of the
training data quantity and prediction horizon on the performance measured in
terms of $\overline{\textrm{AN-RFMSE}}$.

####  4.2.1 Comparison of sparse and dense rolling forecast

Fig. 12 and Fig. 13 show the performance of 20 sparse and 20 dense DNN models
with different parameter initialization forecasting the state variables
$\mathbf{x}$ in one of the time series in the test set
$\mathcal{S}_{test}(i)=\\{\mathbf{X}_{i}\\}$ as defined in Eq. (33). The
models are trained on a data set
$\mathcal{S}_{train}=\left\\{\\{\mathbf{X}_{1},\;\mathbf{Y}_{1}\\},\;\\{\mathbf{X}_{2},\;\mathbf{Y}_{2}\\},\;...,\;\\{\mathbf{X}_{10},\;\mathbf{Y}_{10}\\}\right\\}$
consisting of ten time series $\\{\mathbf{X}_{1},\;..,\;\mathbf{X}_{10}\\}$
with 999 time steps each.

TruthMean sparsecontrol inputstd. dev. sparse

Figure 12: Sparse rolling forecast of state variables $\\{x_{1},..,x_{8}\\}$
at each time instant. The true values of $\mathbf{x}$ and control inputs
$\mathbf{u}$ are taken from one simulated set of test set trajectories
$\mathbf{X}_{i}\in\mathbf{X}_{test}$. The orange dotted line shows the average
of 20 forecasts calculated by 20 sparse neural network models with different
parameter initialization. The orange band shows the standard deviation of the
same 20 forecasts calculated by 20 sparse models.

TruthMean densecontrol inputstd. dev. dense

Figure 13:  Dense rolling forecast of state variables $\\{x_{1},..,x_{8}\\}$
at each time instant. The true values of $\mathbf{x}$ and control inputs
$\mathbf{u}$ are taken from one simulated set of test set trajectories
$\mathbf{X}_{i}\in\mathbf{X}_{test}$. The orange dashed line shows the average
of 20 forecasts calculated by 20 dense neural network models with different
parameter initialization. The orange band shows the standard deviation of the
same 20 forecasts calculated by 20 dense models. The plots are divided into
two regions on the x-axis with two different scales on the y-axis. The first
region from time step zero until approximately time step 400 has the y-scale
on the left side of the plots. The second region, from time step 500 to time
step 700, has the y-scale on the right side of the plots. The first region
shows the average and the standard deviation before the average forecast
standard deviation drifts exponentially from the true values. The second
region shows the exponential drift of the mean and standard deviation of the
forecasts.

Fig. 12 and Fig. 13 indicate that the forecasts of sparse and dense models are
showing similar performance for the first time steps after they are given the
initial conditions. However, while the forecasts calculated by sparse models
show a consistently slow drift from the simulated values of $\mathbf{x}$, the
mean and standard deviation of forecasts calculated by dense models suddenly
drifts exponentially. The narrow banded standard deviation of sparse neural
networks can indicate that these models converge to models with similar
characteristics during training despite different parameter initialization.
Furthermore, the consistently slow drift between the sparse DNN model forecast
of $\mathbf{x}$ and the true values of $\mathbf{x}$ indicate that the sparse
models are generalizing better as they are showing good forecasting
capabilities in a broader region than the dense DNN models. Fig. 12 and Fig.
13 show some interesting results that indicate better generalization of sparse
DNN models than dense DNN models and that the convergence of model parameters
for sparse DNN models are more robust to random initialization than dense DNN
models are.

####  4.2.2 Impact of the training data quantity and prediction horizon

Fig. 14 shows the median, maximum and minimum elements of the
$\overline{\textrm{AN-RFMSE}}_{vec}$ vector.

Median $\overline{\textrm{AN-RFMSE}}$ sparseMedian $\overline{\textrm{AN-
RFMSE}}$ dense

(a) Median, maximum and minimum $\overline{\textrm{AN-RFMSE}}$ for 200
timesteps of forecasting.

(b) Median, maximum and minimum $\overline{\textrm{AN-RFMSE}}$ for 300
timesteps of forecasting.

(c) Median, maximum and minimum $\overline{\textrm{AN-RFMSE}}$ for 500
timesteps of forecasting.

(d) Median, maximum and minimum $\overline{\textrm{AN-RFMSE}}$ for 1000
timesteps of forecasting.

Figure 14: Median, maximum and minimum $\overline{\textrm{AN-RFMSE}}$. There
are five groups of models of both sparse and dense DNNs trained on data sets
with different sizes. For each group of models, there is a corresponding
colored bar indicating median values and an error bar indicating maximum and
minimum values of the vector $\overline{\textrm{AN-RFMSE}}_{vec}$. Moreover,
the size of the training sets are indicated on the x-axis of the subplots. The
blue bar shows the median $\overline{\textrm{AN-RFMSE}}$ among 20 sparse DNN
models over 20 test sets. The orange bar shows the median
$\overline{\textrm{AN-RFMSE}}$ among 20 dense DNN models over 20 test sets.
For each subfigure, the $\overline{\textrm{AN-RFMSE}}$ are calculated for a
given number of timesteps reported in the captions of each subfigure (14(a) \-
14(d)). Notice the logarithmic scale of the y-axis.

Fig. 14 contains a good amount of information about the model performance of
dense and sparse DNN models. Figs. 14(a) to 14(d) report median and extreme
values of $\overline{\textrm{AN-RFMSE}}_{vec}$ over four different time
horizons for five groups of models trained on five data sets with different
sizes. Hence, the results in Fig. 14 show how dense and sparse DNN models
perform with varying amounts of training data over varying time horizons.
Figs. 14(a) to 14(d) show that the groups of sparse models trained on data
sets with varying size show similar results, both in terms of median
$\overline{\textrm{AN-RFMSE}}$ and extreme values. However, as Fig. 14(d)
shows, there seems to be a small trend that groups with more training data
perform slightly better over longer forecasting horizons. Furthermore, the
band between the minimum and maximum values of $\overline{\textrm{AN-RFMSE}}$
is overall relatively small for all groups of models and all forecast horizons
for sparse models. The converging behavior of the performance of groups of
sparse models as a function of the amount of data in the training set
indicates that only small amounts of data are required to gain significance
for the model parameters. While the sparse models show stable performance
across groups of models with different amounts of training data and slowly
increasing values of $\overline{\textrm{AN-RFMSE}}$ proportional to the length
of the forecasting horizon, the same cannot be said about the performance of
the dense models. When considering the dense models, Figs. 14(a) to 14(d)
indicate that there is a trend where both median, minimum and maximum values
of $\overline{\textrm{AN-RFMSE}}$ decreases significantly as sizes of training
set decreases. This expected trend indicates that the performance improves
with increasing dataset size. However, the trend is not consistent for all
groups of dense DNN models for all forecasting horizons. Furthermore, the
maximum values of $\overline{\textrm{AN-RFMSE}}$ for groups of dense DNN
models for longer forecasting horizons such as in Fig. 14(c) and Fig. 14(d)
show that $\overline{\textrm{AN-RFMSE}}$ exponentially increases for some of
the models within the groups. This indicates that the dense DNN models are
likely to have some input regions where the model output is not sound. If the
model estimate enters a region with poorly modeled dynamics, the model
estimate might drift exponentially. For short-term prediction, that is in Fig.
14(a) and Fig. 14(b), the trend is that dense models show better performance
for median and minimum values, especially within the groups with large
training sets. This may be because dense models have more flexibility in terms
of more parameters. Furthermore, the sparse regularizers defined in Section
2.2 are guiding the structure of the network. Sparse regularizers are based on
assumptions about the model structure that might ignore some of the dynamics
in the system. However, for longer forecasting horizons (Fig. 14(c) and Fig.
14(d)), sparse models are always showing better performance than dense models
in terms of median $\overline{\textrm{AN-RFMSE}}$. This is a typical example
of a bias-variance trade-off. For all forecasting horizons and within all
groups of training set sizes, sparse models are always showing a smaller
maximum value of $\overline{\textrm{AN-RFMSE}}$.

###  4.3 Training stability perspective

Sparsification on a large ensemble of neural networks gives similar sparse
structures. This has been shown in the structure plots in the Figs.4-11.
Furthermore, Table 5 show that sparse models to a large extent finds the same
feature basis for each of the model outputs
$\\{\hat{f}_{1},\;...,\;\hat{f}_{8}\\}$.

Moreover, Fig. 14 clearly indicates for all forecasting horizons that the
uncertainty bounds for dense models are much larger than those for sparse
models. For the longer horizons, some of the dense models tend to blow up.
This is probably due to that the model enters a region of the input space
where it overfits. This can be seen as poor generalization to that specific
area. However, for shorter prediction horizons such as in Fig.14(a), the large
bound between minimum and maximum $\overline{\textrm{AN-RFMSE}}$ is not of an
order that indicates a blow up. Still, this uncertainty bounds are larger for
dense models than for sparse models. This indicates that the sparse models are
more likely to converge to a similar minimum than the dense models. The
results presented in Fig. 14 also indicates that the sparse models are likely
to converge already for small amounts of data. This is not the case for dense
models.

##  5 Conclusions and future work

This article presents a sparse neural network model that approximates a set of
nonlinear ODEs based on time series data sampled from the system variables
included in the ODEs. The set of nonlinear ODEs represents an aluminum
electrolysis simulator based on the mass and energy balance of the process.
This includes nonlinear and interrelated models of electrochemical and thermal
subprocesses. The sparsity in the model is achieved by imposing sparsity using
$\ell_{1}$ regularization on the weights. The main takeaway from the work can
be concluded as follows:

  * 1.

$\ell_{1}$ regularization drastically reduces the number of parameters in the
DNN. In our case we witnessed a 93% reduction in the parameters compared to
the corresponding dense neural network.

  * 2.

The sparse neural network was more interpretable using the domain knowledge of
the aluminium electrolysis process. In contrast the dense neural networks were
completely black-box.

  * 3.

Sparse neural networks were consistently more stable compared to the dense
counterparts. This was reflected in the model uncertainty estimates based on a
large ensemble of models.

  * 4.

For short prediction horizon, the median accuracy of dense models are better
than for sparse models. However, the uncertainty bounds of dense models are
much larger.

  * 5.

For longer prediction horizons, sparse models outperforms dense models both in
terms of higher median accuracy and lower uncertainty bounds.

While these sparse models show promising results within interpretability and
generalizability, there is still a high potential for improvement. There is a
desire to increase prediction accuracy and decrease the bias of the sparse
models. This might be addressed by investigating other sparsity structures at
different layers that better compromise the bias-variance trade-off. One
possible direction is to inject simplified theories known from first principle
into the neural network to possibly increase accuracy.

##  6 Acknowledgments

This work was supported by the industry partners Borregaard, Elkem, Hydro,
Yara and the Research Council of Norway through the project TAPI: Towards
Autonomy in Process Industries, project number 294544.

## References

  * [1] K. Grotheim, H. Kvande, Introduction to Aluminium Electrolysis-Understanding the Hall-Heroult Process, Aluminium-Verlag, Dusseldorf, Germany, 1993.
  * [2] V. Gusberti, D. S. Severo, B. J. Welch, M. Skyllas-Kazacos, Modeling the mass and energy balance of different aluminium smelting cell technologies, in: Light Metals 2012, Springer, 2012, pp. 929–934.
  * [3] K. E. Einarsrud, I. Eick, W. Bai, Y. Feng, J. Hua, P. J. Witt, Towards a coupled multi-scale, multi-physics simulation framework for aluminium electrolysis, Applied Mathematical Modelling 44 (2017) 3–24.
  * [4] P. Mandin, R. Wüthrich, H. Roustan, Industrial aluminium production: the hall-heroult process modelling, ECS Transactions 19 (26) (2009) 1.
  * [5] A. Naimi, J. Deng, A. Abdulrahman, V. Vajpayee, V. Becerra, N. Bausch, Dynamic neural network-based system identification of a pressurized water reactor, 2020, p. 100 – 104.
  * [6] E. M. Rentería-Vargas, C. J. Zúñiga Aguilar, J. Y. Rumbo Morales, F. D. J. S. Vázquez, M. De-La-Torre, J. A. Cervantes, E. S. Bustos, M. Calixto Rodríguez, Neural network-based identification of a psa process for production and purification of bioethanol, IEEE Access 10 (2022) 27771–27782.
  * [7] D. Fooshee, A. Mood, E. Gutman, M. Tavakoli, G. Urban, F. Liu, N. Huynh, D. Van Vranken, P. Baldi, Deep learning for chemical reaction prediction, Mol. Syst. Des. Eng. 3 (2018) 442–452.
  * [8] S. S. Blakseth, A. Rasheed, T. Kvamsdal, O. San, Deep neural network enabled corrective source term approach to hybrid analysis and modeling, Neural Networks 146 (2022) 181–199.
  * [9] O. San, A. Rasheed, T. Kvamsdal, Hybrid analysis and modeling, eclecticism, and multifidelity computing toward digital twin revolution (2021).  arXiv:2103.14629.
  * [10] Y. Zhang, P. Tiňo, A. Leonardis, K. Tang, A survey on neural network interpretability (2021).  arXiv:2012.14261.
  * [11] E. T. B. Lundby, A. Rasheed, J. T. Gravdahl, I. J. Halvorsen, A novel hybrid analysis and modeling approach applied to aluminum electrolysis process, Journal of Process Control 105 (2021) 62–77.
  * [12] M. Raissi, P. Perdikaris, G. Karniadakis, Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations, Journal of Computational Physics 378 (2019) 686–707.
  * [13] S. Pawar, O. San, B. Aksoylu, A. Rasheed, T. Kvamsdal, Physics guided machine learning using simplified theories, Physics of Fluids 33 (1) (2021) 011701.
  * [14] P. R. Chermont, F. M. Soares, R. C. De Oliveira, Simulations on the bath chemistry variables using neural networks, in: TMS Light Metals, Vol. 2016-January, 2016.
  * [15] A. M. F. de Souza, F. M. Soares, M. A. G. de Castro, N. F. Nagem, A. H. d. J. Bitencourt, C. d. M. Affonso, R. C. L. de Oliveira, Soft sensors in the primary aluminum production process based on neural networks using clustering methods, Sensors 19 (23) (2019).
  * [16] D. Bhattacharyay, D. Kocaefe, Y. Kocaefe, B. Morais, An artificial neural network model for predicting the co 2 reactivity of carbon anodes used in the primary aluminum production, Neural computing and Applications 28 (3) (2017) 553–563.
  * [17] I. Goodfellow, Y. Bengio, A. Courville, Deep learning, MIT press, 2016.
  * [18] J. Frankle, M. Carbin, The lottery ticket hypothesis: Finding sparse, trainable neural networks, arXiv preprint arXiv:1803.03635 (2018).
  * [19] B. K. Natarajan, Sparse approximate solutions to linear systems, SIAM journal on computing 24 (2) (1995) 227–234.
  * [20] G. F. Montufar, R. Pascanu, K. Cho, Y. Bengio, On the number of linear regions of deep neural networks, in: Z. Ghahramani, M. Welling, C. Cortes, N. Lawrence, K. Q. Weinberger (Eds.), Advances in Neural Information Processing Systems, Vol. 27, Curran Associates, Inc., 2014.
  * [21] R. Pascanu, G. Montufar, Y. Bengio, On the number of response regions of deep feed forward networks with piece-wise linear activations, ICLR, 2014.
  * [22] T. Serra, C. Tjandraatmadja, S. Ramalingam, Bounding and counting linear regions of deep neural networks, CoRR abs/1711.02114 (2017).  arXiv:1711.02114.
  * [23] 14 - general physical properties, in: W. Gale, T. Totemeier (Eds.), Smithells Metals Reference Book (Eighth Edition), eighth edition Edition, Butterworth-Heinemann, Oxford, 2004, pp. 14–1–14–45.
  * [24] I. Dincer, C. Zamfirescu, Chapter 1 - fundamentals of thermodynamics, in: I. Dincer, C. Zamfirescu (Eds.), Advanced Power Generation Systems, Elsevier, Boston, 2014, pp. 1–53.
  * [25] S. Skogestad, Chemical and energy process engineering, CRC press, 2008.
  * [26] A dynamic model for the energy balance of an electrolysis cell, Chemical Engineering Research and Design 74 (8) (1996) 913–933.
  * [27] T. Drengstig, On process model representation and alf3 dynamics of aluminum electrolysis cells, Ph.D. thesis, Ph. D. thesis, Norwegian Univ. of Science and Tech.(NUST) (1997).
  * [28] S. W. Jessen, Mathematical modeling of a hall héroult aluminium reduction cell, Master’s thesis, Technical University of Denmark, DTU, DK-2800 Kgs. Lyngby, Denmark (2008).
  * [29] T. M. Hyde, B. J. Welch, The gas under anodes in aluminium smelting cells. part i: Measuring and modelling bubble resistance under horizontally oriented electrodes, LIGHT METALS-WARRENDALE- (1997) 333–340.
  * [30] A. Solheim, S. Rolseth, E. Skybakmoen, L. Støen, Å. Sterten, T. Støre, Liquidus temperature and alumina solubility in the system na 3 alf 6-alf 3-lif-caf 2-mgf 2, in: Essential Readings in Light Metals: Aluminum Reduction Technology, Volume 2, John Wiley & Sons, 2013, pp. 73–82.

##  Appendix A Simulation model

In this section we will follow a purely physics-based approach to deriving the
equations. At appropriate places we will highlight the challenges and
assumption

The dynamical system simulated in this study is generated by the set of
ordinary differential equations (ODE’s) in Eqs. (17) - (24). This system of
equations is derived from a simplified model of an aluminum electrolysis cell.
This model comprises simplified energy and mass balance of the electrolysis
cell. The model consider an energy balance based on sideways heat transfer,
energy transfer between side ledge and bath due to melting and freezing of
cryolite, and energy input as a function of electrical resistance in the
electrolyte and voltage drop due to bubbles. The mass balance includes mass
transfer between side ledge and bath, input of $Al_{2}O_{3}$ and $AlF_{3}$,
production of metal and consumption of the raw material $Al_{2}O_{3}$ and
tapping of metal from the electrolysis cell. In Table 1, the system states and
inputs are described. The purpose of the simulation model is not to mimic the
exact dynamics of an aluminum electrolysis cell, but rather to generate
nonlinear dynamics similar to what occurs in a real aluminum electrolysis.

###  A.1 Heat capacity

Heat capacity is a measure of the amount of thermal energy a body of a certain
material can store for a given temperature and volume and is given by the
definition [23]:

$C=\frac{\delta q}{\delta T}.$ (39)

$C[J/^{\circ}C]$ is the heat capacity, $\delta q[J]$ is an infinitesimal heat
quantity and $\delta T[^{\circ}C]$. Specific heat capacity $c_{p}$ is heat
capacity at constant pressure per unit of mass:

$c_{p}=\left(\frac{dh}{dT}\right)_{p},$ (40)

where $c_{p}\;[J/(kg^{\circ}C]$ is the specific heat capacity, $h\;[J/kg]$ is
specific enthalpy and $T\;[^{\circ}C]$ is temperature. The subscript $p$
indicates constant pressure. In the process of aluminum electrolysis, the
pressure can be assumed to be constant at $p=1\;[atm]$.

###  A.2 Energy and mass balance

The first law of thermodynamics known as the energy conservation principle
states the following [24]:

$\frac{dE_{i}}{dt}=\dot{E}_{in,\;i}-\dot{E}_{out,\;i}.$ (41)

$\frac{dE_{i}}{dt}\;[W]$ is the change of energy of species $i$ in the system,
$\dot{E}_{in,\;i}\;[W]$ is the energy input rate and $\dot{E}_{out,\;i}\;[W]$
is the energy output rate of species $i$ the system. System is here used
synonymous to control volume. The energy of the system can be transferred
through heat, work or through the energy associated with the mass crossing the
system boundary. This can be expressed as follows:

$\displaystyle\dot{E}_{in,\;i}$ $\displaystyle=$
$\displaystyle\dot{Q}_{in,\;i}+\dot{W}_{in,\;i}+\dot{m}_{in,\;i}e_{in,\;i}$
(42) $\displaystyle\dot{E}_{out,\;i}$ $\displaystyle=$
$\displaystyle\dot{Q}_{out,\;i}+\dot{W}_{out,\;i}+\dot{m}_{out,\;i}e_{out,\;i}$
(43)

where $\dot{Q}_{in,\;i}\;[W]$ and $\dot{Q}_{out,\;i}\;[W]$ are the rates of
heat in and out of the system and $\dot{W}_{in,\;i}\;[W]$ and
$\dot{W}_{out,\;i}\;[W]$ is the rate of work generated on the system.
$\dot{m}_{in,\;i}\;[kg/s]$ and $\dot{m}_{out,\;i}\;[kg/s]$ is the mass rate
into the system and out of the system respectively, whereas
$e_{in,\;i}\;[J/kg]$ and $e_{out,\;i}\;[J/kg]$ is the specific energy of the
mass entering and leaving the system. The specific energy can be formulated
as:

$e=u+\frac{1}{2}v^{2}+gz,$ (44)

where $u\;[J/kg]$ is the specific internal energy, $\frac{1}{2}v^{2}\;[J/kg]$
is the specific energy related to velocity $v\;[m/s]$, and $gz\;[J/kg]$ is the
specific energy related to elevation difference $z\;[m]$.

The change of system energy $\frac{dE_{i}}{dt}=\dot{E}_{i}$ can be written as:

$\frac{dE_{i}}{dt}=\frac{d(m_{i}e_{i})}{dt},$ (45)

where $m\;[kg]$ is the mass of the system and $e_{i}\;[J/kg]$ is the specific
energy of the system. Since the relevant control volumes are related to an
aluminum electrolysis, it is reasonable to neglect the terms
$\frac{1}{2}v^{2}$ and $gz$. In this work, $\dot{Q}=Q_{in,\;i}-Q_{out,\;i}$ is
defined as positive when net heat is provided to the system, and
$W=W_{in,\;i}-W_{out,\;i}$ is positive when work is added to the system. Thus,
the resulting energy equation can be formulated as:

$\frac{d(m_{i}u_{i})}{dt}=m\frac{du_{i}}{dt}+u_{i}\frac{dm_{i}}{dt}=\dot{m}_{in,\;i}u_{in,\;i}-\dot{m}_{out,\;i}u_{out,\;i}+\dot{Q}+\dot{W}.$
(46)

Work $W\;[J]$ is organized transfer of energy. $W$ can be divided into several
types of work [25]:

$W=W_{flow}+W_{\Delta V}+W_{s}+W_{el}+W_{other}.$ (47)

$W_{flow}$ is the work associated with the volume displacements of streams
that enter and exit the system, $W_{\Delta V}$ is the work associated with
changes of the volume of the system, $W_{s}$ is the mechanical work supplied
using movable machinery, $W_{el}$ is the electrochemical work supplied when
the system is connected to an external electric circuit. $W_{other}$ is the
sum of other types of work, for example if surface areas changes or
electromagnetic work. For an aluminum electrolysis, $W_{\Delta
V}=W_{s}=W_{other}\approx 0$. $W_{flow}$ is given by:

$W_{flow}=pV,$ (48)

where $p$ is pressure and $V$ is volume. Enthalpy $H\;[J]$ is given by:

$H=U+pV.$ (49)

where $U=m\cdot u\;[J]$. Furthermore, $H=h\cdot m$. Thus:

$m\frac{du_{i}}{dt}+u_{i}\frac{dm_{i}}{dt}=\dot{m}_{in,\;i}h_{in,\;i}-\dot{m}_{out,\;i}h_{out,\;i}+\dot{Q}+\dot{W}_{el}.$
(50)

Assuming that the flow work is neglectable compared to the other quantities in
the energy equation for aluminum electrolysis gives that $H\approx U$. Recall
that $c_{p}=\frac{dh}{dT}$. Hence:

$\frac{du_{i}}{dt}\approx\frac{dh_{i}}{dt}=\frac{\partial h_{i}}{\partial
T_{i}}\frac{dT_{i}}{dt}=c_{p_{i}}\frac{dT_{i}}{dt},$ (51)

where $\frac{dT}{dt}=\dot{T_{i}}\;[^{\circ}C]$ is the temperature derivative
with respect to time. This yields:

$\frac{dT_{i}}{dt}=\frac{1}{m_{i}c_{p_{i}}}\left(\dot{m}_{in,\;i}h_{in,\;i}-\dot{m}_{out,\;i}h_{out,\;i}+\dot{Q}+\dot{W_{el}}-u_{i}\frac{dm_{i}}{dt}\right)$
(52)

The mass rate equation can be formulated as:

$\frac{dm_{i}}{dt}=\dot{m}_{in}-\dot{m}_{out}+\sum_{j=1}^{n_{j}}r_{i,j},$ (53)

where $r_{i,j}$ is the reaction rate of species $i$ being produced or consumed
in reaction $j$. Assuming that the contents of the control volume to be
perfectly mixed and assuming that the flow work is neglectable gives:

$u_{i}\approx h_{i}=h_{out,\;i}.$ (54)

Hence, the resulting temperature specific energy equation for component $i$ in
a control volume is given by:

$\frac{dT_{i}}{dt}=\frac{1}{m_{i}c_{p_{i}}}\left(\dot{m}_{in,i}(h_{in,i}-h_{out,i})+\dot{Q}+\dot{W}_{el}-h_{out,\;i}\sum_{j=1}^{n_{j}}r_{i,\;j}.\right)$
(55)

The latter equation states that the time derivative of the temperature in the
control volume is dependent on the composition of species in the control
volume. It is assumed that the temperature is equal for all components in the
control volume. Furthermore, it is assumed that there is a common heat loss
$\dot{Q}$ from a control volume to other control volumes, and that electrical
power $\dot{W}_{el}$ is performed on the whole control volume instead of on
different components in the control volume. When different components are
mixed in a control volume, the enthalpy of this mix is more complex than
adding the enthalpy of individual species. However, the complexity of mixed
enthalpy is left out of this simulation model. The heat capacity of a mix of
components in a control volume $c_{p_{cv}}$ is simplified to be constant
despite of that $c_{p_{cv}}$ varies with composition and temperature in the
control volume. The values for different species and control volumes are taken
from [26] and [27]. Thus, the simplified simulation equation for the
temperature derivative in a control volume is given by:

$\displaystyle\frac{dT_{cv}}{dt}$
$\displaystyle=\,\frac{1}{m_{cv}c_{p_{cv}}}\left(\left[\sum_{i=1}^{n_{i}}\dot{m}_{in,i}(h_{in,i}-h_{out,i})\right]+\dot{Q}+\dot{W}_{el}\right.$
$\displaystyle\,\left.-\sum_{i=1}^{n_{i}}\left[h_{out,\;i}\sum_{j=1}^{n_{j}}r_{i,\;j}\right]\right)$
(56)

where $m_{cp}$ is the sum of masses in the control volume.

###  A.3 Heat transfer

Heat transfer $\dot{Q}[W]$ will from this point be referred to as $Q$, meaning
that the dot is omitted. the In the process of aluminum electrolysis, the two
most important principles for heat transfer are convection and conduction.
Conduction is heat transfer through molecular motion within a solid material.
The expression for conduction is given by

$Q=-k\cdot A\cdot\frac{\partial T}{\partial x}.$ (57)

$Q\;[W]$ is heat transferred, $A\;[m^{2}]$ is the area the heat is transferred
through $\frac{\partial T}{\partial x}\;[^{\circ}C/m]$ is the temperature
gradient in the direction $x$ that the heat is transferred, and
$k\;[W/(m^{\circ}C)]$ is the thermal conductivity, a material dependent
proportionality constant. For a fixed cross-section area, the one dimensional
steady state heat flow through a wall of thickness $x\;[m]$ from $x=0$ with
temperature $T_{1}$ to $x=1$ with temperature $T_{1}$ integrates to:

$Q=k\cdot A\cdot\frac{T_{1}-T_{2}}{x},$ (58)

where $T_{1}>T_{2}$. Thermal conductive resistance for a plane wall can be
extracted from the latter expression:

$R_{cond}=\frac{x}{k\cdot A},$ (59)

where $R_{cond}[^{\circ}C/W]$ is the thermal resistance, $x[m]$ is the
thickness of the solid material in the direction heat is transferred, and $k$
and $A$ are as mentioned above. Thermal conductive analysis is analogous to an
electrical circuit, where the temperature difference is analogous to the
potential difference $V$, the heat flow is analogous to the electrical current
$I$ and thermal resistance is analogous to electrical resistance $R_{el}$.
Convection is the heat transfer through the mass motion of a fluid. Heat
transfer between a surface at temperature $T_{s}$ and a fluid at a bulk
temperature $T_{f}$ is due to convection. Convection can be formulated as:

$Q=h\cdot A\cdot(T_{s}-T_{f}),$ (60)

where $A[m^{2}]$ is the contact surface between a solid surface and the
liquid, the the heat transfer coefficient $h\;[W/(m^{2})^{\circ}C]$ is the
proportionality constant between the heat flux and the thermodynamic driving
force for the flow of heat, i.e. the temperature difference $(T_{s}-T_{f})$.

Thermal resistance can be defined for a fluid $R_{conv}$, and is given by:

$R_{conv}=\frac{1}{h\cdot A}.$ (61)

As for electrical circuits, thermal resistances can be coupled in series, and
the reciprocal of the total resistance equals the sum of reciprocals of
individual resistances:

$\frac{1}{R_{tot}}=\sum_{i=1}^{N}\frac{1}{R_{i}}$ (62)

Eq. 62 together with the assumption of stationary heat transfer makes it
possible to calculate the heat transfer from one edge to the other through
several resistors in series as the temperature difference between the edges
divided by the sum of reciprocals of individual resistors:

$Q=\frac{T_{1}-T_{N+1}}{\sum_{i=1}^{N}R_{i}},$ (63)

where $T_{1}>T_{2}>...>T_{N+1}$ is temperature and $R_{i}$ are the resistors.
In the simulation model, the heat transfer is assumed to be piecewise
stationary, meaning that the heat transfer is assumed to be constant from the
middle of one control volume to the middle of the adjacent control volume.
However, the heat transfer is not assumed to be stationary through several
control volumes. Thus, there are separate energy balances for each control
volume. Heat transfer is only considered through the side walls of the
electrolysis cell.

Figure 15: Convection and conduction trough several materials

Convective heat transfer $Q_{bath-liq}$ from the bath to the surface of the
side ledge

$Q_{bath-liq}=h_{bath-sl}A_{sl}(T_{bath}-T_{liq}).$ (64)

Conductive transfer $Q_{liq-sl}$ from surface of side ledge to the center of
the side ledge

$Q_{liq-sl}=\frac{2k_{sl}A_{sl}(T_{liq}-T_{sl})}{x_{sl}}.$ (65)

Conductive heat transfer $Q_{sl-wall}$ from center of side ledge to the center
of the side wall

$Q_{sl-
wall}=\frac{A_{sl}(T_{sl}-T_{wall})}{(x_{wall}/2k_{wall})+(x_{sl}/2k_{sl})}$
(66)

The heat transfer from the middle of the wall to the ambient $Q_{wall-0}$
consists of conductive heat transfer from the middle of the wall to the
surface of the wall, and the convection from the surface of the wall to the
ambient air

$Q_{wall-0}=\frac{A_{sl}(T_{wall}-T_{0})}{(1/h_{wall-0})+(x_{wall}/2k_{wall})}.$
(67)

###  A.4 Electrochemical power

Electrochemical power $\dot{W}_{el}\;[W]$, from now referred to as $P_{el}$ is
the amount of energy transferred to a system from a electrical circuit and is
defined as:

$P_{el}=U_{cell}\cdot I_{line},$ (68)

where $U_{cell}\;[V]$ is the applied cell voltage and $I_{line}\;[A]$ is the
line current sent through the electrolyte. The cell voltage is composed of
three different types of voltage contributions. these are the decomposition
voltage, which is the theoretical minimum potential for the decomposition of
alumina, overvoltage, meaning the excess voltage due to electrode polarization
and ohmic voltage drops, due to resistance of various sections in the cell
[28]. These contributions can be divided into smaller contributions caused by
different effects in different parts of the cells. To make the mathematical
expression in the resulting nonlinear simulation model less comprehensive,
only ohmic voltage drop contributions are included. These are:

  * 1.

Electrolyte voltage drop $U_{el}\;[V]$

  * 2.

Bubble voltage drop $U_{bub}\;[V]$

The voltage drop over the electrolyte is due to the electrical resistivity of
the electrolyte. Assuming uniform current density, the resistance of the
elecrolyte is given by:

$R_{el}=\frac{1}{\kappa}\frac{d}{A}.$ (69)

$R_{el}\;[\Omega]$ is the electrical resistance, $\kappa\;[1/(\Omega m)]$ is
electrical conductivity, $d\;[m]$ is the interpolar distance and $A\;[m^{2}]$
is the total surface of the anodes. The expression for electrical conductivity
is given by [1]:

$\displaystyle\kappa$
$\displaystyle=\,exp\biggl{(}2.0156-\frac{2068.4}{T_{bath}+273}+0.4349\cdot
BR-2.07C_{al_{2}O_{3}}-0.5C_{CaF_{2}}$
$\displaystyle\,-1.66C_{MgF_{2}}+1.78C_{LiF}+0.77C_{Li_{3}AlF_{6}}\biggr{)}.$
(70)

$T_{bath}\;[^{\circ}C]$ is the temperature of the electrolyte, $BR\;[-]$ is
the bath ratio, while $C_{x}[-]$ is the concentration of substance $x$. $BR$
is assumed to be constant at $1.2$, $C_{MgF_{2}}=0.01$, $C_{CaF_{2}}=0.05$,
$C_{LiF}=0$ and $C_{Li_{3}AlF_{6}}=0$. Thus, $\kappa$ can be simplified to:

$\kappa=exp\biggl{(}2.496-\frac{2068.4}{T_{bath}+273}-2.07C_{Al_{2}O_{3}}\biggr{)}.$
(71)

The voltage drop due to resistance in the electrolyte is given by:

$U_{EL}=R_{EL}\cdot I_{line}$ (72)

Gas accumulation beneath the anode surface which reduces the cross-sectional
area of the electrolyte in that zone. Thus, the effective resistivity
increases and causes the so called bubble voltage drop $U_{bub}$ [29]:

$U_{bub}=\frac{d_{bub}\cdot j_{A}}{\kappa}\frac{\phi}{1-\phi}.$ (73)

$j_{A}\;[A/cm^{2}]$ is the anode current density $d_{bub}\;[cm]$ in the bubble
layer thickness and $\phi\;[-]$ is the bubble coverage as a fraction of the
anode:

$d_{bub}=\frac{0.5517+j_{A}}{1+2.167j_{A}},$ (74)

and

$\displaystyle\phi$
$\displaystyle=\,0.509+0.1823j_{A}-0.1723j_{A}^{2}+0.05504j_{A}^{3}$
$\displaystyle\,+\frac{0.4322-0.3781BR}{1-1.637BR}+\frac{0.431-0.1437(x_{Al_{2}O3}-x_{Al_{2}O3}^{AE})}{1+7.353(x_{Al_{2}O3}-x_{Al_{2}O3}^{AE})}$
(75)

$x_{Al_{2}O_{3}}\;[-]$ is the weight percent of alumina in the bath and
$x_{Al_{2}O_{3}}^{AE}\;[-]$ is the weight percent of alumina at where the
anode effect occurs, in this case it is assumed that
$x_{Al_{2}O_{3}}^{AE}=2.0$. Since the simulation model is simplified to only
include contributions from $U_{bub}$ and $U_{EL}$, the total applied cell
voltage in the simulation model is given by:

$U_{cell}=U_{EL}+U_{bub}.$ (76)

###  A.5 Mass rates

The substances considered in the simulation model are alumina ($Al_{2}O_{3}$),
aluminum fluoride ($AlF_{3}$) and cryolite ($Na_{3}AlF_{6}$) in the bath and
liquid aluminum ($Al$) in the metal pad below the bath. Aluminum is extracted
from alumina, which is dissolved in the electrolytic bath. In addition to
alumina, carbon anodes are consumed in the net reaction, producing molten
aluminum and carbon dioxide gas ($CO_{2}$):

$2Al_{2}O_{3}+3C\rightarrow 3CO_{2}+4Al.$ (77)

This reaction occurs at a rate according to Farraday’s law for aluminum
electrolysis [27]:

$r_{al}=\frac{CE\cdot I_{line}}{z\cdot F},$ (78)

where $r_{al}\;[mol/s]$ is the reaction rate of the primary reaction of the
Hall-Héroult process presented in Eq. 77, $z=12$ is the number of electrones
in the reaction, $F=96486.7[(A\cdot s)/mol]$ is the Faraday constant, and
$CE[-]$ is the current efficiency, assumed constant at $CE=0.95$. The fed
alumina contains several impurities [27]. In the simulation model derived used
in this work, only sodium oxide ($Na_{2}O$) is considered as additions in the
feeded alumina. The content of $Na_{2}O$ in alumina react as:

$3Na_{2}O+4AlF_{3}\rightarrow 2Na_{3}AlF_{6}+Al_{2}O_{3}.$ (79)

$3$ Mol $Na_{2}O$ reacts with $4$ Mol of $AlF_{3}$ and produces $2$ Mol of
$Na_{3}AlF_{6}$ and $1$ Mol of $Al_{2}O_{3}$. The reaction rate of the latter
reaction $r_{bath}\;[kmol/s]$ can be formulated as:

$r_{bath}=\frac{C_{Na_{2}O}}{3M_{Na_{2}O}}u_{Al_{2}O_{3}},$ (80)

where $C_{Na_{2}O}\;[-]$ is the weight percent of $Na_{2}O$ in the alumina
feed, $M_{Na_{2}O}\;[g/mol]$ is the molar mass of $Na_{2}O$ and
$u_{Al_{2}O_{3}}\;kg/s$ is the rate of alumina feed. The reaction in Eq. 79
affects the mass balance of both $AlF_{3}$, $Na_{3}AlF_{6}$ and $Al_{2}O_{3}$.
Therefore, $r_{bath}$ is included in the mass balance equations of all these
species. The general mass rate Eq. 53 is used in the derivation the mass rate
of side ledge, cryolite, alumina, aluminum fluoride and metal. The control
volumes in which there are nonzero mass rates are the bath/electrolyte, the
metal pad and the side ledge. The mass rates in the electrolyte are:

$\dot{m}_{Al_{2}O_{3}}=(1-C_{Na_{2}O})u_{Al_{2}O_{3}}-\frac{2}{1000}r_{al}M_{Al_{2}O_{3}}+r_{bath}M_{Al_{2}O_{3}}.$
(81)

$\dot{m}_{Al_{2}O_{3}}\;[kg/s]$ is the mass rate of alumina and
$M_{Al_{2}O_{3}}\;[g/mol]$ is the molar mass of alumina.
$\frac{2}{1000}r_{al}M_{Al_{2}O_{3}}\;[kg/s]$ is the reaction rate of alumina
produced due to the reaction in Eq. 77 and $r_{bath}M_{Al_{2}O_{3}}\;[kg/s]$
is the reaction rate of alumina due to the reaction in Eq. 79. The mass rate
of $AlF_{3}$ is given by:

$\dot{m}_{AlF_{3}}=u_{AlF_{3}}-4r_{bath}M_{AlF_{3}},$ (82)

where $\dot{m}_{AlF_{3}}\;[kg/s]$ is the mass rate of aluminum fluoride,
$u_{AlF_{3}}\;[kg/s]$ is the input rate of aluminum flouride and
$4r_{bath}M_{AlF_{3}}\;[kg/s]$ is the reaction rate of produced aluminum
fluoride from the reaction in equation (79). The mass rate of cryolite in the
bath is given by:

$\dot{m}_{cry}=w_{fus}+2r_{bath}M_{cry}.$ (83)

$\dot{m}_{cry}\;[kg/s]$ is the mass rate of cryolite in the electrolyte,
$2r_{bath}M_{cry}\;[kg/s]$ is the reaction rate of produced cryolite due to
the reaction in Eq. 79 and $w_{fus}\;[kg/s]$ is the mass rate of cryolite
transferred between the side ledge and the bath. $w_{fus}$ is given by:

$w_{fus}=\frac{Q_{bath-liq}-Q_{liq-sl}}{\Delta_{fus}H_{cry}}.$ (84)

$Q_{bath-liq}$ and $Q_{liq-sl}$ is given in Eqs. 64 and (65) respectivelym and
$\Delta_{fus}H_{cry}$ is the heat of fusion for cryolite, i.e. the amount of
energy required to melt one $kg$ of cryolite. The heat of fusion for cryolite
at $1000^{\circ}C$ is $\Delta_{fus}H_{cry}=119495\;[J/kg]$ [27], and is
assumed to be constant in the simulation model. The side ledge is necessary to
withstand the highly corrosive molten cryolite in the oven. The side ledge
consists of frozen cryolite [27]. The mass rate of side ledge is therefore the
transfer of cryolite between the electrolyte and side ledge due to melting and
freezing:

$\dot{m}_{sl}=-w_{fus}.$ (85)

$\dot{m}_{sl}\;[kg/s]$ is the mass rate of side ledge, and $w_{fus}\;[kg/s]$
is given above. The mass rate of aluminum is given by:

$\dot{m}_{Al}=\frac{2}{1000}r_{al}M_{Al}-u_{tap}.$ (86)

$\dot{m}_{Al}\;[kg/s]$ is the mass rate of aluminum,
$\frac{2}{1000}r_{al}M_{Al}\;[kg/s]$ is the reaction rate of produced aluminum
due to the reaction in equation (77), and $u_{tap}\;[kg/s]$ is the control
input of tapping metal from the oven.

###  A.6 Temperature derivatives

Eq. (56) is used to calculate the temperature derivatives in the electrolyte,
side ledge and side wall. As menitoned above, $\dot{Q}=Q$ and
$\dot{W}_{el}=P_{el}$. In the electrolyte, the energy is transferred in and
out of the control volume in many different ways. Heat $Q_{bath-sl}$ is
tranferred through convection from the bath to the side ledge surface
($Q_{bath-liq}$) and from the side ledge surface to the center of the side
ledge with conductive heat $Q_{liq-sl}$. The resulting heat transfer can be
formulated as:

$Q_{bath-sl}=\frac{T_{bath}-T_{sl}}{(x_{sl}/2k_{sl}A_{sl})+1/(h_{bath-
sl}A_{sl})}.$ (87)

Energy is transferred through mass transfer in several ways. Energy needed to
heat and melt substances fed as control input $u$ is given by:

$\displaystyle E_{u}$ $\displaystyle=\
\Delta_{fus}H_{Al_{2}O_{3}}u_{Al_{2}O_{3}}+\Delta_{fus}H_{AlF_{3}}u_{AlF_{3}}$
$\displaystyle\
+(T_{bath}-T_{in})(\bar{c}_{p_{Al_{2}O_{3}}}u_{Al_{2}O_{3}}+\bar{c}_{p_{AlF_{3}}}u_{AlF_{3}}).$
(88)

$\Delta_{fus}H_{i}\;[J/kg]$ is the specific heat of fusion for substance $i$,
and $\bar{c}_{p_{i}}\;[J/(^{\circ}Ckg)]$ is the average heat capacity from
$T_{in}$ to $T_{bath}$. Energy is also transferred through mass transfer
between the electrolyte and the side ledge. When side ledge (frozen cryolite)
melts into the bath, energy is required both to heat and melt the frozen
cryolite. The energy required to heat the frozen cryolite is:

$E_{tc,\;liq}=w_{fus}c_{p,cry,liq}(T_{bath}-T_{liq}).$ (89)

$w_{fus}\;[kg/s]$ is as mentioned above the mass rate of cryolite between bath
and side ledge, and is positive when cryolite melts and is transferred to the
bath. $c_{p,cry,liq}\;[J/(^{\circ}Ckg)]$ is the heat capacity of molten
cryolite $T_{bath}\;[^{\circ}C]$ is the bath temperature and
$T_{liq}\;[^{\circ}C]$ is the liquidus temperature at which cryolite melts and
freezes. $E_{tc,\;liq}\;[W]$ is the energy required to heat molten cryolite
from liquidus temperature to bath temperature. The subscript $tc$ stands for
temperature change and the subscript $liq$ indicates that the substance is
liquid. When $E_{tc,liq}$ is positive, it is becauce $W_{fus}$ is positive,
indicating that cryolite is melting. Thus, when $E_{tc,liq}$ is negative,
cryolite is freezing. The energy needed to melt frozen cryolite is given by:

$E_{sc}=w_{fus}\Delta_{fus}H_{cry}.$ (90)

$E_{sc}\;[W]$ is the energy required to melt the mass of frozen electrolyte.
The subscript $sc$ stands for state change, meaning that it transitions
between solid and liquid. $\Delta_{fus}H_{cry}\;[J/kg]$ is the heat of fusion
for cryolite as mentioned above. When $E_{sc}$ is positive, cryolite is melted
and energy is required whereas when $E_{sc}$ is negative, cryolite freezes,
and energy is released. Energy is required for the primary reaction in Eq. 77
since it is endothermic:

$E_{r_{al}}=r_{al}\Delta H_{r_{al}}.$ (91)

$E_{r_{al}}\;[W]$ is the amount of energy required for the reaction to take
place, $r_{al}[mol/s]$ is as mentioned above the reaction rate of the primary
reaction in the process, and $\Delta H_{r_{al}}=2197582\;[J/mol]$ is the
enthalpy of reaction for (77). Since the reaction in equation (79) is
exothermic, this reaction releases energy to its surroundings. This is given
by:

$E_{r_{bath}}=1000r_{bath}\Delta H_{r_{bath}}.$ (92)

$E_{r_{bath}}\;[W]$ is the energy released due to the reaction in (79)
$r_{bath}\;[kmol/s]=1000[mol/s]$ is the reaction rate of (79) and $\Delta
H_{r_{bath}}=-993283\;[J/mol]$ is the entalpy of reation for (79). The time
derivative for the bath temperature is given by:

$\dot{T}_{bath}=\frac{P_{el}-Q_{bath-
sl}-E_{u}-E_{tc,\;liq}-E_{sc}-E_{r_{al}}-E_{r_{bath}}}{m_{bath}c_{p_{bath,\;liq}}}.$
(93)

$m_{bath}\;[kg]$ is the mass of the liquid bath and
$c_{p_{bath,\;liq}}\;[J/(^{\circ}C\cdot kg)]$ is the enthalpy of the liquid
bath. For the side ledge control volume, the energy transfer is through
melting and freezing of cryolite and heat transfer. The energy transfer
related to melting and freezing of cryolite is given by:

$E_{tc\;sol}=w_{fus}c_{p_{cry,\;sol}}(T_{liq}-T_{sl}).$ (94)

$E_{tc\;sol}\;[W]$ is the energy required into the side ledge when frozen
cryolite heats from side ledge temperature $T_{sl}$ to liquidus temperature
$T_{liq}$. $w_{fus}\;[kg/s]$ is given above and
$c_{p_{cry,\;sol}}\;[J/(^{\circ}C\cdot kg)]$ is the heat capacity of solid
cryolite. The time derivative of the side ledge temperature is given by:

$\dot{T}_{sl}=\frac{Q_{liq-sl}-Q_{sl-
wall}-E_{tc\;sol}}{m_{sl}c_{p_{cry,\;sol}}}$ (95)

$\dot{T}_{sl}\;[^{\circ}C/s]$ is the temperature change in the side ledge,
$Q_{liq-sl}\;[W]$ and $Q_{sl-wall}\;[W]$ is the heat in and out of the side
ledge respectively. $m_{sl}\;[kg]$ is the mass of the side ledge, and
$c_{p_{cry,\;sol}}\;[J/(^{\circ}C\cdot kg)]$ is the heat capacity of solid
cryolite. There is no mass transfer through the wall. Therefore, the only
energy transfer through this control volume is through heat transfer. The time
derivative of the wall temperature $\dot{T}_{wall}\;[\leavevmode\nobreak\
\circ C/s]$ is given by:

$\dot{T}_{wall}=\frac{Q_{sl-wall}-Q_{wall-0}}{m_{wall}c_{p_{wall}}}.$ (96)

$Q_{sl-wall}\;[W]$ is the heat from the side ledge to the wall,
$Q_{wall-0}\;[W]$ is the heat from the wall to the ambient, $m_{wall}\;[kg]$
is the mass og the wall and $c_{p_{wall}}\;[J/(^{\circ}C\cdot kg)]$ is the
heat capacity of the wall.

###  A.7 Liquidus temperature

In [30], the liquidus temperature $T_{liq}$ was determined for primary
crystalization of cryolite $(Na_{3}AlF_{6})$ in a system consisting of the
bath components $Na_{3}AlF_{6}-AlF_{3}-LiF-CaF_{2}-MgF_{2}-KF$. The liquidus
temperature was determined by thermal analysis in a vertical alumina tube
furnace under argon atmosphere. An empirical cure was fitted, which is valid
from temperatures $1011^{\circ}C$ to approximately $800^{\circ}C$. The curve
is given by:

$\displaystyle T_{liq}$ $\displaystyle=\
1011+0.50[AlF_{3}]-0.13[AlF_{3}]^{2.2}$ $\displaystyle\
-\frac{3.45[CaF_{2}]}{1+0.0173[CaF_{2}]}$ $\displaystyle\
+0.124[CaF_{2}][AlF_{3}]-0.00542\left([CaF_{2}][AlF_{3}]\right)^{1.5}$
$\displaystyle\
-\frac{7.93[Al_{2}O_{3}]}{1+0.0936[Al_{2}O_{3}]-0.0017[Al_{2}O_{3}]^{2}-0.0023[AlF_{3}][Al_{2}O_{3}]}$
$\displaystyle\ -\frac{8.90[LiF]}{1+0.0047[LiF]+0.0010[AlF_{3}]^{2}}$
$\displaystyle\ -3.95[MgF_{2}]-3.95[KF].$ (97)

$[x]$ denote the weight-$\%$ of component $x$. In the simulator, it is assumed
that the following components are constant at values $[MgF_{2}]=1\%$,
$[CaF_{2}]=5\%$, $[KF]=[LiF]=0\%$. This yields:

$\displaystyle T_{liq}$ $\displaystyle=\
991.2+1.12[AlF_{3}]-0.13[AlF_{3}]^{2.2}+0.061[AlF_{3}]^{1.5}$ $\displaystyle\
-\frac{7.93[Al_{2}O_{3}]}{1+0.0936[AlF_{3}]-0.0017[AlF_{3}]^{2}-0.0023[AlF_{3}][Al_{2}O_{3}]}.$
(98)

###  A.8 Further simplifications in the simulation model

In addition to assumptions and simplifications accounted for in the article,
some additional simplifications are made in the simulation model. The reason
for this is to simplify the analytical expression in the ODE’s used to
simulate the dynamics of the simulation model. The ODE’s will still describe a
complex nonlinear system, but comparing predictive models with the simulation
models, and thus analysing the performance of the novel predictive models will
be clearer. From the expression for $\dot{T}_{bath}$ in (93), the terms
$E_{u},\;E_{sc},\;E_{r_{Al}}$ and $E_{r_{bath}}$ are omitted. Thus, in the
simulation model, the expression for $\dot{T}_{bath}$ is given by:

$\dot{T}_{bath,\;sim}=\frac{P_{el}-Q_{bath-
sl}-E_{tc,\;liq}}{m_{bath}c_{p_{bath,\;liq}}}.$ (99)

This neglects some essential physical effects in the process. This is
justified by the argument that the main purpose of this work is to evaluate
data driven models on a complex nonlinear system, rather than simulating the
dynamics of an aluminum electrolysis cell as good as possible.