Deep active learning for nonlinear system identification

TL;DR

This paper introduces a novel deep active learning framework for nonlinear system identification, combining global exploration and local optimal control to efficiently acquire informative data and improve neural network modeling of dynamical systems.

Contribution

It formulates a static deep active learning approach integrating system dynamics exploration and optimal control for efficient data acquisition in nonlinear system identification.

Findings

Outperforms standard data acquisition methods in simulated case studies.

Uses ensemble neural network variance as an information measure.

Combines global exploration with local optimal control for data collection.

Abstract

The exploding research interest for neural networks in modeling nonlinear dynamical systems is largely explained by the networks' capacity to model complex input-output relations directly from data. However, they typically need vast training data before they can be put to any good use. The data generation process for dynamical systems can be an expensive endeavor both in terms of time and resources. Active learning addresses this shortcoming by acquiring the most informative data, thereby reducing the need to collect enormous datasets. What makes the current work unique is integrating the deep active learning framework into nonlinear system identification. We formulate a general static deep active learning acquisition problem for nonlinear system identification. This is enabled by exploring system dynamics locally in different regions of the input space to obtain a simulated dataset covering the broader input space. This simulated dataset can be used in a static deep active learning acquisition scheme referred to as global explorations. The global exploration acquires a batch of initial states corresponding to the most informative state-action trajectories according to a batch acquisition function. The local exploration solves an optimal control problem, finding the control trajectory that maximizes some measure of information. After a batch of informative initial states is acquired, a new round of local explorations from the initial states in the batch is conducted to obtain a set of corresponding control trajectories that are to be applied on the system dynamics to get data from the system. Information measures used in the acquisition scheme are derived from the predictive variance of an ensemble of neural networks. The novel method outperforms standard data acquisition methods used for system identification of nonlinear dynamical systems in the case study performed on simulated data.

Tables

Table 1. Table 1 : Values and units of the parameters of the vessel.

Constant	Value	Unit	Constant	Value	Unit
$m_{11}$	2389.657	kg	$m_{22}$	2533.911	kg
$m_{23}$	62.386	kg	$m_{32}$	28.141	kg
$m_{33}$	5068.910	$\mathrm{kg}\cdot m^2$	$X_{v_1}$	-27.632	$\mathrm{kg}\cdot s^{-1}$
$X_{\left|v_1\right|v_1}$	-110.064	$\mathrm{kg}\cdot s^{-1}$	$X_{v_1{v}_1{v}_1}$	-13.965	$\mathrm{kg}\cdot s^{-1}$
$Y_{v_2}$	-52.947	$\mathrm{kg}\cdot s^{-1}$	$Y_{\left|v_2\right|v_2}$	-116.486	$\mathrm{kg}\cdot s^{-1}$
$Y_{v_2{v}_2{v}_2}$	-24.313	$\mathrm{kg}\cdot s^{-1}$	$Y_{\left|r\right|v_2}$	-1540.383	$\mathrm{kg}\cdot s^{-1}$
$Y_r$	24.732	$\mathrm{kg}\cdot s^{-1}$	$Y_{\left|v_2\right|r}$	572.141	$\mathrm{kg}\cdot s^{-1}$
$Y_{\left|r\right|r}$	-115.457	$\mathrm{kg}\cdot s^{-1}$	$N_{v_2}$	3.5241	$\mathrm{kg}\cdot s^{-1}$
$N_{\left|v_2\right|v_2}$	-0.832	$\mathrm{kg}\cdot s^{-1}$	$N_{\left|r\right|v_2}$	336.827	$\mathrm{kg}\cdot s^{-1}$
$N_r$	-122.860	$\mathrm{kg}\cdot s^{-1}$	$N_{\left|r\right|r}$	-874.428	$\mathrm{kg}\cdot s^{-1}$
$N_{rrr}$	0.000	$\mathrm{kg}\cdot s^{-1}$	$N_{\left|v_2\right|r}$	-121.957	$\mathrm{kg}\cdot s^{-1}$

Table 2. Table 2 : Initial condition intervals for the states 𝐱 𝐱 \mathbf{x}

Variable	Initial condition interval
$x_1$	$[-0.2,\mathrm{\hspace{0.33em}1.4}]$
$x_2$	$[-0.2,\mathrm{\hspace{0.33em}0.2}]$
$x_3$	$[-0.2,\mathrm{\hspace{0.33em}0.2}]$

Table 3. Table 3 : Reference intervals for the states 𝐱 r ​ e ​ f subscript 𝐱 𝑟 𝑒 𝑓 \mathbf{x}_{ref}

Variable	Reference values interval
$x_1$	$[-0.3,\mathrm{\hspace{0.33em}1.3}]$
$x_2$	$[-0.3,\mathrm{\hspace{0.33em}0.3}]$
$x_3$	$[-0.3,\mathrm{\hspace{0.33em}0.3}]$

Equations

$$\dot{\mathbf{\eta }}$$

$$\mathbf{M}\dot{\mathbf{\nu }}+\mathbf{C}(\mathbf{\nu })\mathbf{\nu }+\mathbf{D}(\mathbf{\nu })\mathbf{\nu }$$

$$\mathbf{M}=\left[\begin{array}{ccc}{m}_{11}& 0& 0\\ 0& m_{22}& m_{23}\\ 0& m_{32}& m_{33}\end{array}\right],$$

$$\begin{array}{rl}{c}_{13}(\mathbf{\nu })& =-m_{22}{v}_2-m_{23}r\\ c_{23}(\mathbf{\nu })& =m_{11}{v}_1\\ c_{31}(\mathbf{\nu })& =-c_{13}(\mathbf{\nu })\\ c_{32}(\mathbf{\nu })& =-c_{23}(\mathbf{\nu })\end{array}\begin{array}{rl}{d}_{11}(\mathbf{\nu })& =-X_{v_1}-X_{|v_1|v_1}|v_1|-X_{v_1{v}_1{v}_1}{v}_1^2\\ d_{22}(\mathbf{\nu })& =-Y_{v_2}-Y_{|v_2|v_2}|v_2|-Y_{|r|v_2}|r|-Y_{v_2{v}_2{v}_2}{v}_2^2\\ d_{23}(\mathbf{\nu })& =-Y_r-Y_{|v_2|r}|v_2|-Y_{|r|r}|r|\\ d_{32}(\mathbf{\nu })& =-N_{v_2}-N_{|v_2|v_2}|v_2|-N_{|r|v_2}|r|\\ d_{33}(\mathbf{\nu })& =-N_r-N_{|v_2|r}|v_2|-N_{|r|r}|r|-N_{rrr}{r}^2\end{array}.$$

$$\begin{array}{rl}\dot{\mathbf{\eta }}& =\mathbf{R}(\psi )\mathbf{v}=\left[\begin{array}{ccc}\cos \psi & -\sin \psi & 0\\ \sin \psi & \cos \psi & 0\\ 0& 0& 1\end{array}\right]\left[\begin{array}{c}{v}_1\\ v_2\\ r\end{array}\right]\end{array},$$

$$\dot{\mathbf{\nu }}={\mathbf{M}}^{-1}\left(\mathbf{\tau }-\mathbf{C}(\mathbf{\nu })\mathbf{\nu }-\mathbf{D}(\mathbf{\nu })\mathbf{\nu }\right)$$

$$\dot{\mathbf{x}}=\mathbf{f}(\mathbf{x},\mathbf{u})\mathrm{with}\mathbf{x}=\mathbf{\nu }\text{and}\mathbf{u}=\mathbf{\tau },$$

$$\dot{\mathbf{x}}=\hat{\mathbf{f}}(\mathbf{x},\mathbf{u};\mathbf{\theta }),$$

$${\mathbf{\theta }}^{\ast }=\underset{\mathbf{\theta }}{\mathrm{argmin}}\left\{\frac{1}{N}\sum _{i=1}^{N}C({\mathbf{x}}_i,{\mathbf{u}}_i,{\dot{\mathbf{x}}}_i,\mathbf{\theta })\right\}.$$

$$C({\mathbf{x}}_i,{\mathbf{u}}_i,{\dot{\mathbf{x}}}_i,\mathbf{\theta })=J({\dot{\mathbf{x}}}_i,\hat{\mathbf{f}}({\mathbf{x}}_i,{\mathbf{u}}_i;\mathbf{\theta }))+\lambda R(\mathbf{\theta }).$$

$$J({\dot{\mathbf{x}}}_i,\hat{\mathbf{f}}({\mathbf{x}}_i,{\mathbf{u}}_i;\mathbf{\theta }))=({\dot{\mathbf{x}}}_i-\hat{\mathbf{f}}({\mathbf{x}}_i,{\mathbf{u}}_i;\mathbf{\theta }){)}^2.$$

$${}^j{\hat{\mathbf{x}}}_{t+1}{=}^j{\hat{\mathbf{x}}}_t+{\hat{\mathbf{f}}}_j{(}^j{\hat{\mathbf{x}}}_t,{\mathbf{u}}_t;\mathbf{\theta })\Delta T,$$

$${\mathbf{\mu }}_{t+1}=\frac{1}{M}\sum _{j=1}^M{}^j{\hat{\mathbf{x}}}_{t+1},$$

$${\mathbf{\sigma }}_{t+1}^2=\frac{1}{n}\sum _j^n{(}^j{\hat{\mathbf{x}}}_{t+1}-{\mathbf{\mu }}_{t+1}{)}^2.$$

$${\mathcal{B}}^{\ast }=\underset{\mathcal{B}\subseteq \mathcal{U}}{\mathrm{argmax}}{a}_{batch}(\mathcal{B},\hat{\mathbf{f}}(\mathcal{L})),$$

$$\{{\mathbf{u}}_t{\}}_{t=0}^{T-1,\ast }$$

$$s.t.{\hat{\mathbf{x}}}_{t+1}$$

$$H(y)=-{\int }_{\mathcal{Y}}p(y)log(p(y))dy,$$

$$H_{Gaussian}(y)=\frac{1}{2}log(2\pi exp({\sigma }^2(y)),$$

$${\mathcal{B}}^{\ast }$$

$${\left\{\{{\mathbf{u}}_t{\}}_{t=0}^{T-1,\ast },\{{\hat{\mathbf{x}}}_t{\}}_{t=1}^T\right\}}_{i=1}^b\subseteq \mathcal{X}\mathcal{U}$$

$${\mathbf{X}}_i=\left[\begin{array}{cccccc}{x}_1(0)& x_2(0)& x_3(0)& u_1(0)& u_2(0)& u_3(0)\\ x_1(1)& x_2(1)& x_3(1)& u_1(1)& u_2(1)& u_3(1)\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ x_1(t)& x_2(t)& x_3(t)& u_1(t)& u_2(t)& u_3(t)\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ x_1(T-1)& x_2(T-1)& x_3(T-1)& u_1(T-1)& u_2(T-1)& u_3(T-1)\\ x_1(T)& x_2(T)& x_3(T)& nan& nan& nan\end{array}\right]$$

$${\mathbf{Y}}_i=\left[\begin{array}{ccc}\frac{x_1(1)-x_1(0)}{\Delta T}& \frac{x_2(1)-x_2(0)}{\Delta T}& \frac{x_3(1)-x_3(0)}{\Delta T}\\ \frac{x_1(2)-x_1(1)}{\Delta T}& \frac{x_2(2)-x_2(1)}{\Delta T}& \frac{x_3(2)-x_3(1)}{\Delta T}\\ ⋮& ⋮& ⋮\\ \frac{x_1(T)-x_1(T-1)}{\Delta T}& \frac{x_2(T)-x_2(T-1)}{\Delta T}& \frac{x_3(T)-x_3(T-1)}{\Delta T}\end{array}\right]$$

$${\mathcal{S}}_i=[{\mathbf{X}}_i^{{}^{\prime }},{\mathbf{Y}}_i].$$

$${\mathcal{S}}_{batch}=[{\mathcal{S}}_1^T,{\mathcal{S}}_2^T,...,{\mathcal{S}}_b^T{]}^T,$$

$$a_{global}\left(\mathcal{B},{\left(\{{\mathbf{u}}_t{\}}_{t=0}^{T-1,\ast },\{{\hat{\mathbf{x}}}_t{\}}_{t=1}^T\right)}_{i=1}^b\right)=\sum _{i=1}^b\sum _{t=1}^{T}H({[}^j{\hat{\mathbf{x}}}_t{]}_{j=1}^M),$$

$$H({[}^j{\hat{\mathbf{x}}}_t{]}_{j=1}^M)\approx {\mathbf{\sigma }}_t^2.$$

$$a_{global}\left(\mathcal{B},{\left(\{{\mathbf{u}}_t{\}}_{t=0}^{T-1,\ast },\{{\hat{\mathbf{x}}}_t{\}}_{t=1}^T\right)}_{i=1}^b\right)=\sum _{i=1}^b\sum _{t=1}^T{\mathbf{\sigma }}_t^2\odot {\mathbf{s}}_{inv},$$

$$a_{local}=\sum _{t=1}^T{\mathbf{\sigma }}_t^2\odot {\mathbf{s}}_{inv}.$$

$$\frac{d\hat{\mathbf{x}}(t_i)}{dt}=\{\begin{array}{ll}\hat{\mathbf{f}}(\hat{\mathbf{x}}(t_i),\mathbf{u}(t_i)),& \text{if }\mathbf{x}(t_i)\text{ is not measured,}\\ \hat{\mathbf{f}}(\mathbf{x}(t_i),\mathbf{u}(t_i)),& \text{if }\mathbf{x}(t_i)\text{ is measured.}\end{array}$$

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

$$\text{AN-RFMSE}=\frac{1}{p}\sum _{i=1}^p\frac{1}{n}\sum _{j=1}^n{\left(\frac{{\hat{x}}_i(t_j)-x_i(t_j)}{std(x_i)}\right)}^2,$$

Full text

††journal: Arxiv

1 Introduction

Modeling dynamical systems is a cornerstone of most engineering applications where the system states of the processes change over time. Prediction models can be utilized in different ways. Models can be used in a control system setting to design safe and optimal control systems. Furthermore, if the prediction models exhibit a high degree of accuracy, they can be utilized to forecast system states over extended periods of time. Accurate simulations of the physical behavior across longer horizons can enhance our understanding of the underlying physical process and support stakeholders in informed decision making. Physics-Based Models (PBMs) are widely utilized in a range of engineering and scientific applications. PBMs are mathematical models derived from fundamental physical principles. These principles are fundamental laws that govern a particular aspect of the natural world. Moreover, these principles or fundamental laws describe the behavior of an observable phenomenon. Our understanding of these physical phenomena is primarily attained through the examination and analysis of observed occurrences of the phenomenon in question. On one side, this means that PBMs possess a high degree of interpretability, and demonstrate generalizability when modeling assumptions are sound. On the other side, the deductive nature of this modeling approach is highly biased, and potentially ignorant to unobserved or unknown physical phenomena. Increased access to abundant data, cheap computational resources, and many achievements and improvements in the Machine Learning (ML) society has created an enormous interest in a variety Data-Driven Models (DDMs) in many engineering fields including material science [1], biomechanics [2], production of biofuels [3], reservoir modeling in oil and gas industry [4], aluminum extraction [5], bioengineering [6], drug discovery [7] and more [8]. DDMs can model an underlying process directly from input-output data. In the scientific community of systems and control, the art of building mathematical models from observed input-output data is called system identification [9]. The research of identifying linear dynamics started in the late 1950’s [10]. Since then, a pool of theories covering linear system identification has been developed. Linear models typically require less data and are structured and well-behaved. Thus a linear model is preferred if it can approximate a system with satisfactory accuracy. However, more complex model structures are typically required to obtain highly accurate process models of complex nonlinear systems. The inescapable nature of a wilde range of real world problems and processes is a core motivation for the broad, and ever evolving field of nonlinear system identification. Determining the model structure for nonlinear problems is a fundamental problem in nonlinear system identification. Using Neural Networks (NNs) structures stands out as highly interesting, due to their remarkable abilities to model complex nonlinear phenomena. This has motivated researchers to investigate NNs as process models in various dynamical systems. In [11], NNs were used to identify the dynamics of a pressurized water nuclear reactor, authors of [12] used NNs to identify the dynamics of a purification process of bioethanol, in [13] a Deep Learning (DL) model was used to predict chemical reactions, and in [14], NNs were used to identify the dynamics of a Quadcopter. Unfortunately, NNs typically require large amounts of diverse data, which is expensive to generate and usually unavailable from dynamical systems. In [15], this issue is addressed by inducing sparsity in the NN models and shows that this reduces the data requirements for obtaining NNs with desirable accuracy. In addition to sparse NNs, [16] show that skip-connections can also contribute to increasing model accuracy when data size is limited. However, large amounts of training data is still needed to achieve satisfactory model accuracy. This may preclude the use of NNs in applications where data acquisition is expensive.

Active Learning (AL) aims to maximize model accuracy with a minimal amount of data by acquiring the most informative training data [17]. AL has been utilized in many fields, including image recognition [18], text classification [19], and object detection [20] to mention a few. In these scenarios, large unlabeled datasets are available, and the AL algorithm aims to choose the most informative samples among the unlabeled data. The goal is then to reduce the costly labeling process which is performed by human domain experts. In the context of dynamical systems, labeling output data may not incur significant expenses. However, data obtained from a processes under closed-loop feedback control often lacks the necessary information to identify NN models with acceptable performance, including accuracy and generalizability to operational regions of the input space, known as the state-action space. Obtaining informative datasets of a dynamical processes is costly due to for example interruption of the production or operation, and expenses of measuring certain states. Moreover, exciting the dynamical systems to regions with high model uncertainty can induce unforeseen and potentially severe incidents to the process and its surroundings. While the safety-critical nature of using NNs in controlled processes is a major research challenge itself, approached by for example techniques within reachability analysis [21, 22], we limit the scope of this work to the informativeness of the sampled data, addressing challenges related to costs of sampling large datasets in dynamical systems.

Introducing AL methods to experimental design for system identification introduces additional challenges to the AL problem. That is, most AL methods address static acquisition problems where, in principle, any location in the input space is directly accessible, or a dataset is sampled in advance. For dynamical systems, on the other hand, reaching a desired location in the state-action space requires system excitation through control inputs. The topic of optimal excitation has been addressed in the research field known as optimal experiment design [23]. In light of this, optimal experiment design can be seen as a subfield of AL, or they can be seen as related research topics. Anyways, AL, which originates from the computer science community provides a wide range of information-theoretic approaches as well as a well-defined learning framework, providing great inspiration to researchers working with system identification.

Authors of [24] propose AL strategies for the identification of a Gaussian Process (GP) model inspired by information theoretic measures. The most promising work they propose suggests optimizing a sequence of control inputs that maximize the predictive differential entropy along the state trajectory; a method outperforming state-of-the-art experimental design methods. The work of identifying GP models was extended in [25], to also include global explorations. The global search for initial states is done by exploring the informativeness of short trajectories from candidate initial states. When an informative initial state is acquired, the local exploration maximizes the predictive entropy along the state trajectory as in [24]. AL is also applied to acquire data that efficiently identify linear models by solving an Optimal Control Problem (OCP) that maximize the minimal eigenvalues of covariates of states [26]. An active learning approach to identify a restricted class of nonlinear dynamical modes whose state transitions depend linearly on a known feature embedding of state-action pairs was investigated in [27]. However, the research on active learning for system identification of NN models is, to the best of the authors’ knowledge, highly limited.

To that end, we extend the work of AL used to acquire the most informative data for system identification to NNs. That is,

In equation (18), we formalize a general Batch Mode Deep Active Learning (BMDAL) acquisition scheme for dynamical system identification referred to in this work as global exploration. The scheme is on the static deep active learning batch acquisition form presented in equation (14) from [17]. The static acquisition scheme for dynamical systems is enabled through local explorations, obtaining a set of simulated informative candidate state-action trajectories distributed around the state-action space. The novel DeepAL scheme builds upon the AL scheme presented in [25] that iteratively searches for the single most informative state-action trajectory for identifying a GP model through local and global explorations.

• 2.

The novel formulation of the static BMDAL acquisition for dynamical systems is utilized in a novel framework presented in Fig. 3. Ensembles of NNs are used to produce uncertainty estimates that assess the informativeness of single state-action trajectories.

• 3.

The general nature of the proposed BMDAL formulation allows for a wide range of query strategies to be applied. This is demonstrated by using two different query strategies in the case study.

The AL algorithm is showcased on the simulated dynamics of a 3 Degree Of Freedom (DOF) surface vessel with three states and three control inputs, yielding an input space of six dimensions. The simulator represents the dynamics of the MilliAmpere ferry[28], which is an experimental platform owed by NTNU. The simulation model is presented in Section 2.1.

2 Theory

2.1 Physics based simulator

We use the standard 3-Degrees of Freedom (3-DOF) model of a marine craft, which is a simplified model of a real vessel. The state of the vessel is described by the pose vector $\bm{\eta}=[x\,y\,\psi]^{\top}\in\mathbb{R}^{3}$ and the velocity vector $\bm{\nu}=[v_{1}\,v_{2}\,r]^{\top}\in\mathbb{R}^{3}$. The pose vector describes the position and orientation of the vessel in the North-East-Down (NED) frame with $x$ and $y$ being the position in the North and East directions, respectively, and $\psi$ being the heading angle. The velocity vector describes the velocity of the vessel in the body frame with $u$ and $v$ being the velocity in the surge and sway directions, respectively, and $r$ being the yaw rate. The model is formulated as a nonlinear system of ordinary differential equations (ODEs) as follows (see [29] for details):

[TABLE]

where $\mathbf{R}(\psi)\in\mathrm{SO}(3)$ is the rotation matrix from the body frame to the NED frame. The mass matrix $\mathbf{M}\in\mathbb{R}^{3\times 3}$, the Coriolis matrix $\mathbf{C}(\bm{\nu})\in\mathbb{R}^{3\times 3}$, and the damping matrix $\mathbf{D}(\bm{\nu})\in\mathbb{R}^{3\times 3}$ are given by [28] and express the inertia and Coriolis matrices as:

[TABLE]

where the elements $c_{ij}(\bm{\nu})$ and $d_{ij}(\bm{\nu})$ are given by

[TABLE]

The constant coefficients in Equation (3) are summarized in Table 1. The kinematics of the vessel is given by

[TABLE]

and its dynamics are governed by

[TABLE]

We are interested in learning the dynamics of the vessel in response to given forces and moments, which we consider as control input. In the remainder of the paper we use the following notation:

[TABLE]

i.e. $\mathbf{x}\in\mathbb{R}^{3}$ is the state vector and $\mathbf{u}\in\mathbb{R}^{3}$ is the control input, with dynamics described by Equation (5).

2.2 Data-driven modeling ODE’s using Neural Networks

Data-Driven Models (DDM)s can approximate an underlying process with data from the process. The underlying process can be represented by a nonlinear mapping from input to output space. Deep feed-forward neural networks known as Multi-Layer Perceptrons (MLP) have high flexibility and, thereby ability to model nonlinear mappings. In this work, the MLP is used to model Eq. Equation (6). The trainable parameters of the NN are denoted by $\bm{\theta}\in\mathbb{R}^{p}$ where $p$ is the number of parameters, and the network is denoted as

[TABLE]

where $\dot{\mathbf{x}}$ is the vector of time derivatives of the states. Optimal model parameters $\bm{\theta}^{*}$ are found according to the optimization of a cost function $C(\cdot)$:

[TABLE]

$N$ is the size of the training set ${\cal D}=\{([\mathbf{x}_{i},\mathbf{u}_{i}],\dot{\mathbf{x}}_{i})\}_{i=1}^{N}$, and cost function $C(\cdot)$ is given by:

[TABLE]

$J(\cdot)$ is a proper loss function. For regression tasks, the loss function is typically chosen to be the Mean Squared Error (MSE):

[TABLE]

$R(\bm{\theta})$ is a regularization term on the parameters. In this work, the $\ell_{1}$ regularization is used due to its sparsity promoting nature. $\lambda$ is a tunable scalar multiplied with the regularization term $R(\bm{\theta})$. The error of the cost function is propagated backward from the output layer, through the hidden layers, and to the input layer in the so-called back-propagation scheme.

2.3 Ensembles of neural networks

Ensemble learning includes methods that combine multiple models in making predictions. The main premise in ensemble learning is that errors made by individual models are likely to be compensated by other models such that the overall ensemble prediction on average improves prediction accuracy over individual models [30]. Deep ensembles, consisting of multiple Deep Neural Networks (DNNs), have gained significant attention in recent years due to their improved accuracy, uncertainty estimation, and robustness to out-of-distribution data. There are two well known methods of training ensembles of NNs; by bootstrapping, where the ensemble methods are trained on different bootstrap samples of the dataset, and by multiple random seeds, where the model parameters of the members are intialized with different random parameters and then trained on the entire dataset. While the bootstrapping method has been found to hurt performance of the NNs, using random initializations turn out to be a promising approach [31, 32]. In [33], the success of random initialization in deep ensembles is explained by the fact that this method explore diverse modes of the function space. That is, Bayesian neural networks, which do not perform as good as deep ensembles, only explore the proximity of one single mode of the function space. Ensembles can provide both point estimates and uncertainty estimates. The point estimate can be calculated as an average of the predictions. Consider an ensemble of $M$ NNs with different parameter initialization. When forecasting several timesteps ahead without feedback from measurements, we let each ensemble member estimate individual state forecast trajectories. The state forecast ${}^{j}\hat{\mathbf{x}}_{t+1}$ at time $t+1$ provided by an ensemble member $\hat{\mathbf{f}}_{j}$ is given by forward Euler integration:

[TABLE]

where $\Delta T$ is the time step. Then, the average prediction of $M$ predictions calculated by individual NN ensemble members at timestep $t+1$ is given by:

[TABLE]

Uncertainty-based strategies are utilized in the proposed AL method. Therefore, the main motivation of using deep ensembles in this work is to obtain an uncertainty estimate. NN predictions have two sources of uncertainty, namely model uncertainty (also known as epistemic uncertainty), and data uncertainty (also known as aleatoric uncertainty). Usually, these types of uncertainties are modeled separately. In this study, the data is sampled from a fully observed process without process disturbances or measurement noise. Thus, only model uncertainty is considered here. The model uncertainty is caused by shortcomings in the model. This includes errors in the training procedure such as bad training hyperparameters (learning rate, batch size, regularization etc.), insufficient model structure, or lack of information in the data [34]. In general, there are four different types of methods of estimating the uncertainty of a NN based on whether the NNs are deterministic or stochastic, and whether a single NN or multiple NNs are used to estimate the uncertainty. These are single deterministic methods, Bayesian methods, ensemble methods and test-time augmentation methods [34]. Ensemble methods have proven to be attractive in quantifying uncertainty of NNs. Ensemble methods are in several works compared to Bayesian methods. In [35], it is argued that ensemble based methods preform better than Bayesian Monte-Carlo Dropout approximations in DeepAL due to more calibrated predictive uncertainties. Both [36] and [37] came to the same conclusions, particularly under dataset shift. A simple way to quantify the uncertainty of these predictions is to calculate the empirical variance of NN predicitions for each output of the network

[TABLE]

Here, ${}^{j}\hat{\mathbf{x}}_{t+1}$ is the NN prediction made by ensemble member $j$ of $\mathbf{x}_{t+1}$, and $boldsymbol{\mu}_{t+1}$ is the mean ensemble prediction at time step $t+1$.

2.4 Deep active learning

Deep Active Learning (DeepAL) has emerged as a combined approach between DL and AL, addressing DL specific challenges within AL. This mainly includes dealing with over-confident uncertainty estimates of NN predictions, efficiently data-acquisition of data batches rather than the traditional AL one-by-one query method, and the joint optimization of the NN model and AL algorithm [17]. The majority of research on DeepAL focuses on static acquisition problems. Static acquisition problems refer to scenarios where data is already available, and any point in the input space can be acquired directly. Examples of such problems are visual data processing such as image classification [38] and object detection [39], Natural Language Processing (NLP) such as machine translation [40], text classification [41] and semantic analysis [42]. The static acquisition problem imply that there exists an unlabeled dataset $\mathcal{U}=\{\mathcal{Z}\}$ with $c$ input samples $\mathcal{Z}=\{\mathbf{z}_{1},\;\mathbf{z}_{2},\;...,\;\mathbf{z}_{c}\}$. The goal of DeepAL in the static acquisition problem is to acquire as few as possible of the unlabeled data in $\mathcal{U}$ for labeling by choosing the most informative samples. That includes designing a query strategy $Q$, $\mathcal{U}\xrightarrow{Q}\mathcal{L}$, where $\mathcal{L}=\{\mathcal{Z},\mathcal{Y}\}$ is a labeled dataset, and $\mathcal{Y}$ are labels corresponding to inputs $\mathcal{Z}$. The query strategy can be expressed in terms of an acquisition function $a_{batch}$ which acquires a batch $\mathcal{B}^{*}=\{z_{1}^{*},\;z_{2}^{*},\;...,\;z_{b}^{*}\}$ of samples to be labeled. The batch based query called BMDAL is the foundation of DeepAL. The DeepAL scheme is an iterative acquisition scheme, and one acquisition step is generally defined by:

[TABLE]

Here $\hat{\mathbf{f}}(\mathcal{\mathcal{L}})$ is the NN. $\mathcal{L}$ is the labeled data up until the given acquisition step, and the notation $\hat{\mathbf{f}}(\mathcal{\mathcal{L}})$ indicate that the NN is trained on this data. The acquisition function $a_{batch}$ is in general a function of the NN that is trained on the currently acuired data since the informativeness of new samples can be evaluated using this NN.

The acquisition function $a_{batch}$ defines the query strategy of the AL scheme. There exists a range of different query strategies in AL. Here we will shortly describe the strategies relevant to the case study. Uncertainty-based strategy is one of the most popular strategies in AL. The strategy aims to select samples in which the model predictions are most uncertain about. Uncertainty-based AL methods are typically computationally efficient and easy to implement. Moreover, these methods typically provides highly informative samples. One of the most utilized uncertainty-based methods calculates the predictive entropy $H[\mathbf{y}|\mathbf{x},\mathcal{L}]$ for a given sample $\mathbf{x}$. However, there are some concerns about applying uncertainty-based sampling strategies in BMDAL. Acquiring a batch of the most informative samples using an uncertainty measure can lead to a batch of very similar samples. Moreover, strategies of this type are often focused on examples close to a decision boundary, making it vulnerable to adverserial attacks [17]. Hence, a Hybrid strategy is often preferred, accounting for diversity in the sampled data. A method called Diverse Mini-Batch Active Learning (DMBAL) [43] adds informativeness to the optimization of a K-means algorithm in the weights of each candidate sample. In the DMBAL algorithm, informative estimates obtained by some informative measures are assigned as weights to the corresponding candidate samples. In each acquisition step, a batch $\mathcal{B}$ of $b$ samples closest to the centroids of the weighted K-means algorithm is added to the training set.

2.5 Deep Active learning in Dynamical systems

Although data from a dynamical system may be readily available from production or operation, it often provides limited information and is not well-suited for the purpose of system identification. Due to the physical nature of dynamical systems, arbitrary points in the state-action space cannot be directly accessed. In order to sample given data points from the state-action space, the dynamics must be excited by control inputs. This is a dynamic acquisition problem. In an attempt to maximize the information contained in this sampling process, an OCP can be defined over a finite horizon, maximizing some measure of information. This is here referred to as the local exploration. As the name indicates, this optimization is only efficient for shorter horizons, limiting the method to explore the dynamics in the proximity of the initial state. However, when identifying the input-output mapping of a nonlinear dynamical system, the entire operational window of the system must be explored.

Assuming a set of input data is already available, AL offers a robust approach for selecting the most informative data points from the input space. In the current acquisition problem for dynamical systems, we do not have access to a pre-sampled dataset. However, by locally exploring different parts of the input space, a set of simulated state-action trajectories can be obtained. With an available simulated dataset, a static AL acquisition problem for dynamical systems can be formulated. This is referred to as global exploration. The global exploration acquires the batch of initial states corresponding to the batch of state-action trajectories that maximizes a global batch acquisition function. Following the acquisition of a set of initial states through global exploration, a subsequent round of local exploration is conducted for each state in the batch. This local exploration entails a longer optimization horizon compared to the initial search conducted for all candidates during the global exploration. This is because the computational complexity of the OCP increases significantly with the horizon, making it necessary to restrict the horizon to a relatively short length when optimizing for all candidates prior to the global exploration. When control trajectories are obtained from the final local explorations these trajectories are applied to the real system from the corresponding acquired initial states. This is done under the assumption that the system is driven to each initial state using a specific control law. As the system evolves under the applied control sequences, data on the system states is collected.

2.5.1 Local exploration

Data sampled from a dynamical system should be properly excited by a control signal to obtain informative data that can be used for system identification. Local exploration refers to the dynamic AL acquisition problem of finding a control trajectory that informatively excites the system. Given an initial state $\mathbf{x}_{0}$ from where the dynamical system is excited, the local exploration can be formulated as an open loop finite horizon OCP , which yields a sequence of control inputs $\{\mathbf{u}\}_{k=0}^{T-1,\;*}=\{\mathbf{u}_{0},\;\mathbf{u}_{1},\;...,\;\mathbf{u}_{T-1}\}^{*}$ . In the context of active learning, the objective function is an acquisition function $a_{local}$ that measures the informativeness of the sequence of forecasted states $\{\hat{\mathbf{x}}_{t}\}_{t=1}^{T}=\{\hat{\mathbf{x}}_{1},\;\hat{\mathbf{x}}_{2},\;...,\;\hat{\mathbf{x}}_{T}\}$ given the candidate sequence of control inputs $\{\mathbf{u}_{k}\}_{k=0}^{T-1}=\{\mathbf{u}_{0},\;\mathbf{u}_{1},\;...,\;\mathbf{u}_{T-1}\}$ and an initial state $\mathbf{x}_{0}$:

[TABLE]

where $\hat{\mathbf{x}}_{0}=\mathbf{x}_{0}$. The standard strategy in Model Predictive Control (MPC) formulation is to only apply the first control input in the sequence and then solve the OCP again for each consecutive timestep until the end of the horizon. This scheme requires $T-1$ optimizations to obtain one sequence of control inputs, and is therefore computational expensive. An alternative that is computationally feasible is to optimize for the entire control sequence one time and apply the control sequence obtained from that one solution of the OCP. The authors of [24] developed an active learning scheme for a GP model. They suggested to maximize the sum of differential entropy of the GP model predictions over the control horizon of $T$ steps, such that $a_{local}=\sum_{i=0}^{T-1}H[\hat{\mathbf{f}}(\hat{\mathbf{x}}_{i},\mathbf{u}_{i};\bm{\theta})]$. The differential entropy of variable $y$ is defined by [44]

[TABLE]

where $p(y)$ is the probability density function. In this case, the probability density function represent the distribution over the predictions. If the variable $y$ is Gaussian distributed, the differential entropy is given by

[TABLE]

where $\sigma^{2}$ is the variance of the given prediction.

2.5.2 Global exploration

Exciting the system dynamics is essential to obtain informative data from a dynamical system. The local exploration formulated as an OCP in equation (15) provides a sound basis for exiting the problem locally. However, when the goal is to obtain the most informative data from the entire input space, solely depending on the optimization in equation (15) is inefficient. That is, the computational complexity increases drastically with optimization horizon. This put restrictions on how long the optimization horizon can be, and therefore also the area that a optimized state-action trajectory can span. Moreover, the uncertainty of state forecasts typically increases with each time step. This is highly relevant in the local exploration formulation since the corresponding OCP typically aims to maximize some uncertainty measure. With high levels of uncertainty, the actual states are likely to deviate from predicted states after longer horizons. Thus, the efficacy of local explorations as defined above is typically limited to exploring dynamics locally. Authors in [25] suggest partitioning the search problem into global and local explorations for actively learning a GP model. Building upon the work in [25], equation (18) provides a general formulation of the DeepAL optimization problem for dynamical systems, acquiring an optimal batch rather than single initial states. The global exploration consider a set $\mathcal{X}=\{\mathbf{x}_{0,\;1},\;\mathbf{x}_{0,\;2},\;...,\;\mathbf{x}_{0,\;c}\}$ of $c$ candidate initial states. For each of the candidate initial states $\mathbf{x}_{0,\;i}$, an optimal control trajectory $(\mathbf{u_{t}}_{t=0}^{T-1,\;*})_{i}$ is obtained by optimizing the OCP in equation (15). With an initial condition $\mathbf{x}_{0,\;i}$ and the obtained control trajectory $(\mathbf{u_{t}}_{t=0}^{T-1,\;*})_{i}$, the corresponding forecasted state trajectory $\left(\{\hat{\mathbf{x}}_{t}\}_{t=1}^{T}\right)_{i}$ is estimated by the model. One acquisition step of the global exploration is generally described in the following DeepAL optimization formulation:

[TABLE]

where $\mathcal{B}=\{\mathbf{x}_{0,\;1},\;\mathbf{x}_{0,\;2},\;...,\;\mathbf{x}_{0,\;b}\}$ is a candidate batch of initial conditions, and $\left\{\{\mathbf{u}_{t}\}_{t=0}^{T-1,\;*},\{\hat{\mathbf{x}}_{t}\}_{t=1}^{T}\right\}_{i=1}^{b}$ is the corresponding batch of simulated state-action trajectories. $\mathcal{XU}=\left\{\left\{\mathbf{u}_{t}\}_{t=0}^{T-1,\;*},\{\hat{\mathbf{x}}_{t}\}_{t=1}^{T}\right\}_{1},\;\left\{\mathbf{u}_{t}\}_{t=0}^{T-1,\;*},\{\hat{\mathbf{x}}_{t}\}_{t=1}^{T}\right\}_{2},\;...,\;\left\{\mathbf{u}_{t}\}_{t=0}^{T-1,\;*},\{\hat{\mathbf{x}}_{t}\}_{t=1}^{T}\right\}_{c}\right\}$ is the set of simulated candidate state-action trajectories. The acquired a batch $\mathcal{B}^{*}=\{\mathbf{x}_{0,\;1}^{*},\;\mathbf{x}_{0,\;2}^{*},\;...,\;\mathbf{x}_{0,\;b}^{*}\}$ of initial conditions corresponds to the batch simulated state-action trajectories $\left(\{\mathbf{u}_{t}\}_{t=0}^{T-1,\;*},\{\hat{\mathbf{x}}_{t}\}_{t=1}^{T}\right)_{i=1}^{b}$ that maximize some global batch acquisition function $a_{global}$. Since simulated state-action trajectories are already sampled from in a local exploration scheme, the global exploration becomes a static acquisition problem on the form of the standard DeepAL scheme presented in equation (14).

3 Method and setup

In this section, the experimental setup, as well as the methods used in the case study is presented. The data is generated by integrating the nonlinear ODEs in equation (6) with a set of initial values for the states $\mathbf{x}_{0}$ using the fourth-order Runge-Kutta (RK4) numerical integration algorithm. In the DeepAL method presented in this work, a batch of initial states are chosen from a set of candidate states according to the optimization in equation (18). The query strategies defined by the global acquisition function $a_{global}$ are described in Section 3.1. The control trajectories that excite the system dynamics are acquired in the local exploration scheme defined in equation (15). The local acquisition function $a_{local}$ in this scheme is an uncertainty-based strategy also described in detail in Section 3.1. The benchmark method chooses the set of initial conditions randomly. Moreover, the control input trajectory from each initial state are chosen according to the Amplitude modulated Pseudo Random Binary Signal (APRBS) used to identify nonlinear dynamics with NNs in works like [45], [16] and [46]. In each loop of the AL scheme, a batch of $b=10$ time series $\{\mathbf{X}_{1},\;\mathbf{X}_{2},\;...,\;\mathbf{X}_{i},\;...,\;\mathbf{X}_{b}\}$ are obtained. The $i^{\prime}th$ time series is obtained by simulating the dynamics from an initial condition $\mathbf{x}_{0}=[x_{1}(0),\;x_{2}(0),\;x_{3}(0)]$ over a horizon of $T=15$ with timesteps $\Delta T=0.1sec$. This yields a time series $\mathbf{X}_{i}$:

[TABLE]

Hence, the control inputs $\mathbf{u}(t)$ are defined until timestep $t=T-1$. That is, at the last step there is no need for a control input since there is no next state to be calculated. For each time series $i$ in the batch, the output label $\mathbf{Y}_{i}$ for training is calculated by the forward Euler formula:

[TABLE]

We define a new matrix $\mathbf{X}_{i}^{{}^{\prime}}$ that contains all but the last row of the $i^{\prime}th$ time series $\mathbf{X}_{i}$. Then $\mathbf{X}_{i}^{{}^{\prime}}$ and $\mathbf{Y}_{i}$ are paired as inputs and outputs:

[TABLE]

This is done for all simulations in the batch. Then the input-output pairs are stacked:

[TABLE]

before added to the training data $\mathcal{D}_{train}$. Before the training is conducted, the inputs and outputs in the training set $\mathcal{D}_{train}$ are normalized, shuffled, and put in mini-batches for training. In each loop of the learning scheme, an ensemble of $M=10$ NNs are trained on all training data sampled up until that time.

3.1 Novel DeepAL scheeme for dynamical systems

The DeepAL acquisition scheme comprises a global exploration scheme and a local exploration scheme. The global exploration scheme will for each acquisition step in the AL loop choose a batch of $b$ initial states $\mathcal{B}^{*}=\{\mathbf{x}_{0,\;1}^{*},\;\mathbf{x}_{0,\;2}^{*},\;...,\;\mathbf{x}_{0,\;b}^{*}\}$ among a set of $c$ candidates $\mathcal{X}=\{\mathbf{x}_{0,\;1},\;\mathbf{x}_{0,\;2},\;...,\;\mathbf{x}_{0,\;c}\}$ according to the AL optimization problem defined in equation (18). For each of the candidate initial states $\mathbf{x}_{0,\;i}$ a state-action trajectory is obtained according to the local exploration in equation (15). The query strategy is defined by the global batch acquisition function $a_{global}$, which quantifies the informativeness of batches of state-action trajectories corresponding to initial state candidates. A simple uncertainty-based acquisition function that sum of predictive entropies along all candidate trajectories is given by:

[TABLE]

where $[^{j}\hat{\mathbf{x}}_{t}]_{j=1}^{M}$ is the set of ensemble forecasts at timestep $t$. Assuming that the ensemble predictions are approximately normally distributed around the mean prediction $\bm{\mu}_{t}$ given in equation (12), and that predicted states are uncorrelated, maximizing the entropy will become approximately the same as to maximizing the empirical variance given in equation (13):

[TABLE]

In order to scale the optimization problem according to the magnitude of states in the state vector, the empirical variance of state $k$, $\sigma_{t,\;k}$, $k\in\{1,\;2,\;3\}$ can be divided by the standard deviation of the $k^{\prime}th$ state based on the currently sampled dataset. defining the vector $\mathbf{s}_{inv}=[\frac{1}{std_{1}},\;\frac{1}{std_{2}},\;\frac{1}{std_{3}}]^{T}$, where $std_{k}$ is the standard deviation of the $k^{\prime}th$ state $x_{k}$, the resulting acquisition function can be defined as:

[TABLE]

where $\odot$ is the Hadamard product operator that takes the element-wise multiplication of the two vectors. The resulting acquisition function is purely uncertainty based and does not take into account the similarity between samples. Hybrid acquisition strategies takes into account both uncertainty of individual samples as well as the similarities between samples in a candidate batch $\mathcal{B}$. An intuitive hybrid acquisition method DMBAL adds informativeness to the optimization of a weighted K-means algorithm. That is, the algorithm acquires the closest sample to each of the $b$ centroids found by a weighted K-means, where the weight is some informative measure. A simple adaption of the algorithm to the problem at hand is given in Algorithm 1:

The method reuses the uncertainty measure defined in equation (25). In addition, the method aims to add diversity of samples by comparing the similarities of the candidate initial conditions in the modified DMBAL method. The $c$ candidates in $\mathcal{X}=\{\mathbf{x}_{0,\;1},\mathbf{x}_{0,\;2},\;...,\mathbf{x}_{0,c}\}$ are at each acquisition step uniformly sampled from the intervals given in Table 2:

The local exploration scheme obtains a sequence of control inputs by optimizing the dynamic acquisition problem formulated as an OCP in equation (15) from a given initial state. The local acquisition function $a_{local}$ used in the optimization defined in equation (15) is the same as the uncertainty-based global acquisition function defined in equation( 25), but only for a single initial state and the corresponding trajectory. That is, the local acquisition function is:

[TABLE]

A schematic illustration of the novel DeepAL scheme is presented in Fig. 3

The OCP solved to generate control sequence is solved in the optimization framework CasADi [47], while the design and training of DNNs used in the optimization is done using the DL framework PyTorch [48]. The ML-CasADi package developed by authors of [49] is used to integrate the two frameworks.

3.2 Performance metrics

We focus on the forecast error for several steps ahead in a so called rolling forecast to measure the generalization error of the NNs. That is, given initial conditions $\mathbf{x}_{0}=\mathbf{x}(t_{0})$ and a sequence of $n$ control inputs $\{\mathbf{u}(t_{0}),\;\mathbf{u}(t_{1}),\;...,\;\mathbf{u}(t_{n-1})\}$, a NN forecast the consecutive $n$ timesteps $\{\hat{\mathbf{x}}(t_{1}),\;\hat{\mathbf{x}}(t_{2}),\;...,\;\hat{\mathbf{x}}(t_{n})\}$. The NN $\hat{\mathbf{f}}$ predict the time derivatives of the states at timestep $i$, that is $d\mathbf{x}(t_{i})/dt$, based on the current state prediction $\hat{\mathbf{x}}(t_{i})$ and control input $\mathbf{u}(t_{i})$:

[TABLE]

Hence, the model either uses the previous calculated forecast of the state $\hat{\mathbf{x}}(t_{i})$ or the true state $\mathbf{x}(t_{i})$ depending on whether the state is measured or not. Here, we assume no measurement noise, making the measured and true state the same. In the case study, the states are measured at every $50$ timestep, meaning that the model forecast for this horizon before it is corrected by measurements. After estimating the time derivative at timestep $i$, the forecast at the next timestep is calculated by the forward Euler formula:

[TABLE]

The rolling forecast accuracy measure used in for example [46], [15] and [16] called Average Normalized Rolling Forecast Mean Squared Error (AN-RFMSE) is a scalar defined by

[TABLE]

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{D}_{train}$, $p=3$ is the number of state variables and $n$ is the number of time steps the normalized rolling forecast MSE is averaged over.

3.3 Test set generation

The utilization of DNNs in modeling dynamical systems is driven by their capability to represent intricate relationships with a high degree of accuracy. When proper measures are taken to address safety considerations, they have the potential to enhance the optimality of MPC. As a result, evaluating the sampling strategies on a testset generated by using an optimal control policy is a subject of significant interest. The MPC used when generating the testset solves an OCP that minimize a quadratic cost function:

[TABLE]

$\mathbf{x}_{ref,\;k}$ is the desired reference signal, and $\mathbf{Q}$ and $\mathbf{R}$ are weighting matrices. The subscript $k$ indicates the value of the variable at timestep $k$. The function $\mathbf{f}(\cdot,\cdot)$ is the simulation model itself. Given an initial condition, the optimization problem in equation Equation (30) is solved for $n$ steps. Both the sequence of states $\{\mathbf{x}_{k}\}_{k=0}^{n-1}$ and control inputs $\{\mathbf{u}_{k}\}_{k=0}^{n-1}$ are decision variables in the optimization, and can be extracted from the solution of the optimization. The testset consist of $100$ time series with initial conditions uniformly sampled from the intervals in Table 2.

Each of the time series are generated by ten optimizations of equation Equation (30) each with a horizon of $50$ timesteps. The final state $\mathbf{x}_{n}$ in one optimization is then the initial state of the next time series, such that the $i^{\prime}th$ time series in the testset can be written as:

[TABLE]

where $n=50$. The value of the references $\mathbf{x}_{ref}$ are uniformly drawn from Table 3 are constant for the optimization horizon. Hence each time series has ten different references over $500$ timesteps. Hence, the testset can be written as:

[TABLE]

4 Results and discussion

The case study presented in this section investigate the efficacy of global and local explorations compared to benchmark random sampling methods. Moreover, the study presents the effect of the global hybrid strategy DMBAL for different values of the prefilter hyperparameter $\beta$, where the special case of $\beta=1$ can be considered as an uncertainty based acquisition strategy.

4.1 Information based and random sampling

In order to quantify the efficacy of global and local explorations in the proposed DeepAL sampling scheme compared to benchmark random sampling methods, we define three fundamental data acquisition schemes. All schemes use both a global-, and a local sampling method. The three fundamental methods are based on either an information theoretic approach or a random sampling strategy for both local and global exploration. Fig. 4 shows a schematic presentation of how the different schemes are derived from the two sampling strategies.

The fundamental sampling schemes are namely Global Informative Local Informative(GI-LI), Global Random Local Informative(GR-LI), and Global Random Local Random(GR-LR). The GI-LI method is described in Section 3.1, and is using information theoretic sampling strategies both locally and globally. The global exploration in GI-LI acquires initial conditions from where to conduct the local explorations. The global exploration is using the local exploration method to assign measures of informativeness to the candidate initial conditions. The GR-LI method is globally random, and a batch of initial states is acquired by uniformly sampling within the interval of states presented in Table 2. The local exploration method of GR-LI is the exact same as the local sampling method used GI-LI. Th GR-LR method uses the same random global strategy as GR-LI, and uses the APRBS method to excite the dynamics locally.

Fig. 5(d) presents the performance of the three sampling schemes for after all data acquisition steps. The error bound show the $25-75$ percentile of AN-RFMSE values. The upper bound is particularly interesting since it gives an intuition about the models ability to generalize to a broader set of the test set trajectories. The results show that the GI-LI method outperforms the GR-LI and the GR-LR methods, both in terms of higher mean accuracy as well as significantly lower values for the $75$ percentile of AN-RFMSE values. Moreover, the globally random, locally informative GR-LI method shows better performance in terms of higher mean accuracy an lower $75$ percentile AN-RFMSE values than the purely random $GR-LR$ method. However, the superiority of GR-LI method over the GR-LR method is not as significant as the superiority of GI-LI over the two others, indicating that the globally informative step is of major importance.

Fig. 6(f) shows the mean and uncertainty bounds of ensemble forecasts trained on data sampled with GI-LI and GR-LR methods. The plots illustrate how the GI-LI might provide data that gives improved mean predictions as well as narrower and better calibrated uncertainty bounds.

4.2 Uncertainty based and hybrid global strategy

Fig. 7 show the performance of GI-LI methods using the global DMBAL method with different prefilter parameter $\beta$ for different stages of the AL scheme. The algorithm only considers the $b\cdot\beta$ samples with highest informativeness score according to a given informativeness measure, where a batch $\mathcal{B}^{*}$ of $b=10$ initial conditions are chosen. Hence, $\beta=1$ corresponds to simply picking the $b$ samples with the highest informativeness score. Since the chosen informative measure is uncertainty based, using $\beta=1$ means that the global method is uncertainty based. Using $\beta>1$ means that the global algorithm takes into account diversity by choosing the samples closest to centroid centers of a k-means algorithm, which means that the overall method is a hybrid AL strategy, combining uncertainty measures and diversity measures. Fig 7 shows that the average performance in all plots is approximately the same for all choices of $\beta$. In Fig. 7(a) and 7(a) showing performance for models trained on data acquired after respectively $10$ to $30$ AL loops and $30$ to $40$ AL loops, it is difficult to conclude any significant differences, other than that the method using $\beta=10$ has higher extremum AN-RFMSE values, as well as more stretched distribution towards higher AN-RFMSE values, indicating that lower choices of $\beta$ give better results. Fig 7(c)-7(f) show the performance of the method for different $\beta$ values for later stages of the AL scheme. The most evident result is the tendency that the upper extreme values of AN-RFMSE have a minimum for $\beta=3$, and that both lower and higher values of $\beta$ give higher extreme AN-RFMSE values. Moreover, it seems like the $95$ percentile also is the lowest for $\beta=3$. Moreover, again, using high values of $\beta$ turns out to give worse result as $\beta=10$ show the worst performance. The overall results from the conducted case study indicate that using the DMBAL algorithm in the global acquisition scheme is of minor importance compared to acquiring a batch of samples with the highest uncertainty scores. However, choosing the right values of the prefilter hyperparameter $\beta$, the modeling errors in extreme cases can potentially be reduced. The limited value of the hybrid method used in this study might be explained by that it only compares the similarity of initial states rather than similarities between candidate trajectories. Moreover, the shortcomings of uncertainty based methods are known to be more present if the batch size is greater. With the currently choice of batch size $b=10$ trajectories being sampled at each acquisition step, the power of hybrid strategies might not be present.

5 Conclusions and future work

The main conclusions from the work can be itemized as follows:

The globally random, locally informative based GR-LI strategy show slightly better results than a globally random, locally random strategy in terms of mean accuracy an lower values of the $75$ procentile of AN-RFMSE values, indicating better generalization. However, the novel GI-LI DeepAL scheme significantly outperform GR-LI and GR-LR schemes both in terms of mean accuracy and $75$ procentile values of AN-RFMSE. This indicates that global explorations are of major importance with respect to achieving higher accuracy and generalization.

• 2.

The DMBAL approach, which emphasizes diversity in the selection of samples, might reduce the upper bound of extreme error values, provided that the prefilter hyperparameter is carefully chosen. This method is compared to simply selecting the top $b$ samples based on some uncertainty measure, when sampling a batch of initial conditions globally. However, the DMBAL approach does not exhibit significant improvements over the global uncertainty based method beyond this in the given case study. The method only compares the similarity between initial states of the candidate trajectories, rather than similarity between the whole candidate trajectories. This leaves out potentially important information.

The novel DeepAL framework is flexible and allows for a range of AL aqcuisition strategies. Conducting a comparative study including different AL acquisition functions in the framework would be highly interesting as it could increase our knowledge about efficient sampling of dynamical systems. Global hybrid strategies that can consider similarity of entire state trajectories are of particular interest since the currently tested hybrid strategy that only compares initial conditions seems to be of limited value.

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.