Model Updating for Nonlinear Systems
with Stability Guarantees

Farhad Ghanipoor^1, Carlos Murguia^1, Peyman Mohajerin Esfahani^2, Nathan van de Wouw^1

1. Department of Mechanical Engineering, Eindhoven University of Technology, The Netherlands
2. Delft Center for Systems and Control, Delft University of Technology, The Netherlands

To improve the predictive capacity of system models in the input-output sense, this paper presents a framework for model updating via learning of modeling uncertainties in locally (and thus also in globally) Lipschitz nonlinear systems. First, we introduce a method to extend an existing known model with an uncertainty model so that stability of the extended model is guaranteed in the sense of set invariance and input-to-state stability. To achieve this, we provide two tractable semi-definite programs. These programs allow obtaining optimal uncertainty model parameters for both locally and globally Lipschitz nonlinear models, given uncertainty and state trajectories. Subsequently, in order to extract this data from the available input-output trajectories, we introduce a filter that incorporates an approximated internal model of the uncertainty and asymptotically estimates uncertainty and state realizations. This filter is also synthesized using semi-definite programs with guaranteed robustness with respect to uncertainty model mismatches, disturbances, and noise. Numerical simulations for a large data-set of a roll plane model of a vehicle illustrate the effectiveness and practicality of the proposed methodology in improving model accuracy, while guaranteeing stability.

§ INTRODUCTION

Dynamical systems modeling has been a key problem in many engineering and scientific fields, such as biology, physics, chemistry, and transportation. When modeling dynamical systems, it is of key importance to use well-established principles of physics and known prior system properties (e.g., stability and set-invariance) [1]. However, for many complex practical systems, we only tend to have partial knowledge of the physics governing their dynamics [2]. Even in cases where accurate physics-based models are established, as, e.g., in robotics, there exist inevitable (both parametric and non-parametric) uncertainties that impact the model's predictive accuracy.

For a class of nonlinear dynamical systems, this paper assumes a prior system model is available (however, results here are also applicable when no prior model is available). It focuses on learning models for uncertainties while guaranteeing the stability of extended models (prior models plus uncertainty representations), given available input-output data. This problem contrasts with black box modeling approaches, such as Neural Networks (NNs) or Gaussian Processes (GPs), because we incorporate prior relationships derived from first principles into both the modeling and learning framework. Moreover, this problem differs from typical grey-box identification problems because, in our case, a prior model with known parameters is available. Off-the-shelf grey-box system identification methods can cope with a subset of the problem under discussion, where no prior model is available. The problem considered here enables the characterization of hybrid system representations (referred here as the extended model) comprising both prior physics-based and uncertainty learned models. In what follows, a comparison of our approach with some related existing literature on hybrid modeling (based on both first-principles and data) is provided.

Existing Literature: Our approach is fundamentally different from existing Physics-Informed (PI) learning techniques where standard black box models (such as NNs or GPs) are trained constrained to satisfy physics-based relations [3, 4, 5]. Yazdani et al. in [3] use this technique to construct the so-called Physics Informed Neural Networks (PINNs). They constrain a known physics-based model during the training of a NN-based model (i.e., penalize loss function for the mismatch between the physics-based model and the NN as a soft constraint) and incorporate physics knowledge in the structure of the NN. Although this method provides a NN as a system model (as well as parameters for physics-based models), it does not give a closed-form expression for the uncertainty in known physics-based models (i.e., it does not account for modeling mismatches due to unmodeled physics). Furthermore, parameters of both physics-based models and NNs are learned simultaneously, which increases the computational burden.

The approach proposed in this paper also differs from the so-called Sparse Identification of Nonlinear Dynamics (SINDy) scheme [6]. SINDy assumes full knowledge of system states and their time derivatives. Then, based on known physics-based models and known variables (i.e., system states and their derivatives) a library of functions is generated that can be incorporated in dynamical models to account for uncertainty. To select the active functions in models, sparse identification algorithms are exploited. This approach has demonstrated accurate performance in sparse model identification of complex nonlinear systems [7, 8, 9]. However, SINDy not only requires full-state measurement but also requires the derivative of states to be known. Although the state derivatives can be approximated numerically if the complete state is known, most numerical methods are noise sensitive. Furthermore, the requirement of full-state measurements is a strong assumption for most dynamical systems. In our work, we do not require measurements of the full-state and its time derivative. The proposed algorithms need input-output data only.

Our approach augments a known physics-based model by a black-box model used as a correction term, see, e.g., [10, 11] and references therein. Such generic approach is also taken by Quaghebeur et al. in [12], who add an NN model to a known physics-based model with unknown parameters. This approach allows maintaining the basic structure of the model that comes from first principles, which improves interpretability. However, simulating the hybrid model at each iteration during the training process is necessary. This approach is evidently more computationally intensive compared to our proposed method, which eliminates the requirement of simulating the model in every iteration. Moreover, the main drawback of this method is that it assumes that the initial state of the dynamic system is known or at least it requires measuring all the states (full-state measurement) of the true dynamical system. This assumption is dropped in the proposed method considered here.

Another important advantage of the proposed method is guaranteeing the stability of the extended nonlinear model (i.e., the model consisting of the known physics-based model and the uncertainty model). The identification of stable models has (mainly) been widely studied in the context of discrete LTI systems [13, 14]. For further results on uncertainty learning for LTI systems, refer to [15]. However, the identification of nonlinear stable models is still under study, mainly focusing on identifying the complete dynamics in a black box fashion. For instance, kernel- or Koopman-based methods that enforce some form of model stability during the learning process are proposed in [16, 17] for the autonomous case and for the non-autonomous case in [18]. Moreover, there are results that aim to enforce model stability in other types of models, such as recurrent equilibrium network models [19] and Lur’e-type models [20, 21]. However, it is important to note that none of these methods can be directly applied to the specific problem we are addressing here, given the prior known model.

In this paper, we propose a framework for model updating via learning modeling uncertainties in (physics-based) models applicable to locally Lipschitz nonlinear systems. We first focus on learning uncertainty models, assuming that some realizations of input, estimated uncertainty, and estimated state are given. During uncertainty learning, we guarantee that trajectories of the extended model (i.e., known prior model plus uncertainty model) belong to a given invariant set for locally Lipschitz extended models, or ensure input-to-state stability for globally Lipschitz extended models. This is achieved by formulating the problem as a constrained supervised learning problem.

One key challenge in this problem involves the introduction of stability constraints, which is tackled using Lyapunov-based tools. The stability criteria usually result in an optimization problem that is non-convex. We address this challenge by proposing two different approaches for both locally and globally Lipschitz models:

* Cost Modification: The first approach involves a change of variables, which leads to the rewriting of cost function (Theorems <ref> and <ref> for locally and globally Lipschitz models, respectively).

* Constraint Modification: This approach introduces a sufficient (convex) condition to satisfy the stability constraint (Theorems <ref> and <ref> for locally and globally Lipschitz models, respectively).

For the sake of completeness and comparison, we also provide a sequential algorithm that alternates the use of some variables in the optimization problem to convexity the program. However, we show in the numerical section that (as it is known in related literature [23]) that initializing this algorithms is challenging, which further strengthens the importance of the  results provided here. We referred to the above mentioned sequential algorithm as the method of Sequential Convex Programming (SCP).

After addressing the challenge of non-convexity, the paper proceeds to discuss the practical implementation of the framework. It outlines a method for estimating uncertainty and state trajectories using input-output data and the known prior model (Proposition <ref>). In this context, uncertainty is considered as an unknown input affecting the system dynamics and the estimation of uncertainty and state trajectories is achieved using robust state and unknown input observers [22]. For a schematic overview of the proposed methodology, Model Updating for Nonlinear Systems with Stability guarantees (MUNSyS), see Figure <ref>.

Overview of the MUNSyS methodology.

In summary, the main contributions of this paper are as follows:

* Practically Applicable Framework:  The proposed framework for learning modeling uncertainties in locally (and globally) Lipschitz nonlinear models only requires input-output data. Another important aspect for practical applicability is the lower computational cost compared to existing methods. This is achieved by eliminating the requirement to simulate the model in every iteration.

* Stability Guarantee of Extended Models:  The proposed framework guarantees the following: 1. For locally Lipschitz extended models, trajectories of the extended model (i.e., known prior model plus uncertainty model) belong to a given invariant set; 2. For globally Lipschitz models, it ensures Input-to-State Stability (ISS).

* Convex Approximations for Non-Convex Programs: Two distinct approaches for both locally and globally Lipschitz models, namely cost and constraint modification, are proposed to offer tractable convex approximations of formulated non-convex optimization problems for uncertainty model learning, while ensuring stability guarantees.

This paper generalizes the preliminary results published in the conference paper [15]. In comparison to [15], which is applicable to LTI systems only, we present results for a broader class of systems, considering both locally and globally Lipschitz nonlinear systems. The results in [15] are a subclass of the results provided for globally Lipschitz nonlinear systems (see Theorem <ref>).

The set of nonnegative real numbers is represented by the symbol $\mathbb{R}^{+}$. The identity matrix of size $n \times n$ is denoted as $I_n$ or simply $I$ when the context specifies $n$. Similarly, matrices of dimensions $n \times m$ comprising only zeros are denoted as ${0}_{n \times m}$ or ${0}$ when the dimensions are clear. The first and second time-derivatives of a vector $x$ are denoted as $\dot{x}$ and $\ddot{x}$, respectively. For the $r^{th}$-order time-derivatives of vector $x$, the notation $x^{(r)}$ is employed. A positive definite matrix is symbolized by $X \succ 0$, and positive semi-definite matrices are indicated by $X \succeq 0$. Similarly, $X \prec 0$ is used for negative definite matrices, and $X \preceq 0$ for negative semi-definite matrices. The Hadamard (element-wise) power  of $n \in \mathbb{N}$ for a matrix $X$ is denoted by $X^{\circ n}$. The same notation is used for vectors. The $exp(x)$ shows the element-wise exponential function. The notation $(x_1, \ldots , x_n)$ signifies the column vector composed of elements $x_1, \ldots, x_n$, and this notation is also used when the components $x_i$ are vectors. Both Euclidean norm and the matrix norm induced by Euclidean norm are represented by the notation $||\cdot||$. The infinity-norm is denoted as $||\cdot||_\infty$. We employ the notation $\mathcal{L}_2(0, T)$ (or simply $\mathcal{L}_2$) to represent vector-valued functions $z: [0, T] \to \mathbb{R}^{k}$ satisfying $||z(t)||_{\mathcal{L}_2}^2 := \int_{0}^{T}||z(t)||^{2} dt < \infty$.

For a vector-valued signal $z$ defined for all $t \geq 0$, the $\mathcal{L}_\infty$ norm is denoted as $||z||_{\mathcal{L}_\infty} := \sup_{t \geq 0} ||z(t)||$. For a differentiable function $V: \mathbb{R}^{n} \to \mathbb{R}$, the row-vector of partial derivatives is denoted as $\frac{\partial V}{\partial x}$, and $\dot{V}(x)$ denotes the total derivative of $V(x)$ with respect to time (i.e., $\frac{\partial V}{\partial x} \frac{dx}{dt}$). The trace of a matrix $W$ is denoted as $tr(W)$. A continuous function $\bar \alpha:[0, a) \rightarrow[0, \infty)$ is said to belong to class $\mathcal{K}$ if it is strictly increasing and $\bar \alpha(0) = 0$. A continuous function $\bar \beta:[0, a) \times[0, \infty) \rightarrow[0, \infty)$ is said to belong to class $\mathcal{KL}$ if, for each fixed $s$, the mapping $\bar \beta(r,s)$ belongs to class  $\mathcal{K}$ with respect to $r$ and, for each fixed $r$, the mapping $\bar \beta(r,s)$ is decreasing with respect to $s$ and $\bar \beta(r, s) \rightarrow 0$ as $s \rightarrow \infty$. Time dependencies are often omitted for simplicity in notation.

§ PROBLEM FORMULATION

Consider the nonlinear system

\begin{equation} \label{eq:sys}
\left\{\begin{aligned}
\dot{x}_s &=A x_s+ B_u u+ S_g g(V_g x_s, u) + S_\eta \eta(V_\eta x_s, u) + B_\omega \omega,  \\
y_s&= C x_s + D_\nu \nu,
\end{aligned}\right.
\end{equation}

where $x_s \in {\mathbb{R}^{n}}$, $y_s \in {\mathbb{R}^{{m}}}$, and $u \in {\mathbb{R}^{{l}}}$ are system state, measured output and known input vectors, respectively. Function $g: \mathbb{R}^{n_{v_g}} \times \mathbb{R}^{l} \to \mathbb{R}^{n_g}$ is a known nonlinear vector field, and function $\eta: \mathbb{R}^{n_{v_\eta}} \times \mathbb{R}^{l} \to \mathbb{R}^{n_\eta}$ denotes unknown modeling uncertainty. Signals $\omega: \mathbb{R}^{+}  \to \mathbb{R}^{n_\omega}$ and $\nu: \mathbb{R}^{+}  \to \mathbb{R}^{m_\nu}$ are unknown bounded disturbances; the former with unknown frequency range and the latter typically with high-frequency content (e.g., related to measurement noise). Known matrices $(A, B_u, S_g, V_g, S_\eta, V_\eta, B_\omega, C, D_\nu)$ are of appropriate dimensions, with ${n},{m}, l, n_{v_g}, n_g, n_{v_\eta}, n_{\eta}, n_{\omega}, m_\nu \in \mathbb{N}$. The matrices $S_g$ and $S_\eta$ specify the equations explicitly incorporating the nonlinearity $g$ and uncertainty $\eta$, while the matrices $V_g$ and $V_\eta$ identify the states influencing the nonlinearity and uncertainty, respectively. Moreover, without loss of generality, we assume that zero is an equilibrium point of the system for $u = 0$.

In the following, we state the problem we aim to solve at a high abstraction level.

(Uncertainty Learning with Stability Guarantees) Consider system (<ref>) with known input and output signals, $u$ and $y_s$. We aim to learn a data-based model for the uncertainty function $\eta(\cdot)$ so that the extended model composed of the known part of (<ref>) and the learned uncertainty model are “stable". The objective is to construct a more accurate system model (at least applicable to trajectories close to the training data set).

We first assume the availability of a labeled data-set containing input, estimated uncertainty, and estimated state realizations. A method to obtained this data from input-output trajectories is given in Section <ref>. In what follows, we outline the problem settings for both locally and globally Lipschitz systems.

§.§ Problem Settings

We consider two model classes. For both model classes, the extended model (for the system in (<ref>)) is of the form:

\begin{equation} \label{eq:model}
\begin{aligned}
\dot{x} &= A x+ B_u u + S_g g(V_g x, u) + S_{\eta_l} \eta_l(V_\eta x, u), \\
\eta_l(V_\eta x,u) &:= \Theta_l V_\eta x+B_l  u +  \Theta_n h(V_\eta x,u),
\end{aligned}
\end{equation}

where $x \in {\mathbb{R}^{n}}$ is model state and function $\eta_l: \mathbb{R}^{n_{v_\eta}} \times \mathbb{R}^{l} \to \mathbb{R}^{n_{\eta_l}}$ is the uncertainty model that is parameterized by $\Theta_l, B_l$ and $\Theta_n$. Function $h: \mathbb{R}^{n_{v_\eta}} \times \mathbb{R}^{l} \to \mathbb{R}^{n_{h}}$ is a given nonlinear vector filed. This function contains the vector of basis functions that serve as candidates for the nonlinearities in the uncertainty model. Design matrices are collected as $\theta = (\Theta_l, B_l, \Theta_n)$ and have appropriate dimensions, $n_{\eta_l}, n_h \in \mathbb{N}$. Matrix $S_{\eta_l}$, similar to $S_\eta$ in (<ref>), indicates the explicit presence of the uncertainty model $\eta_l$ in the right-hand side and may differ from $S_{\eta}$.

§.§.§ Cost Function

Recall that, for now, we assume that estimated uncertainty and state realizations are available. Define the data vector corresponding to the $i$-th sample in time as

\begin{equation*}
\begin{aligned}
d_i := {\left[\begin{array}{cccc}
\hat{x}_i^\top V_\eta^\top & u_i^\top  & \hat{\eta}_i^\top  & h(V_\eta \hat{x}_i, u_i)^\top
\end{array}\right]^\top},
\end{aligned}
\end{equation*}

where $\hat{x}_i, u_i,$ and $\hat{\eta}_i$ represent the state estimate, input, and uncertainty estimate, respectively. Vector $h(V_\eta \hat{x}_i, u_i)$ corresponds to the evaluated nonlinearity in (<ref>) at the given realizations (of state estimation and input). Given $N$ samples of the data vector defined above, define the data matrix $D$ as follows:

\begin{equation} \label{eq:data_matrix}
D := \sum_{i = 1}^{N} d_i d_i^\top.
\end{equation}

Further, consider the error vector between the uncertainty model and its (given) estimate as

\begin{equation} \label{eq:T}
\begin{aligned}
e_i := \eta_l(V_\eta \hat{x}_i, u_i)  - \hat{\eta}_i = T d_i, \quad \text{where}  \quad T := {\left[\begin{array}{cccc}
\Theta_l  & B_l & -I & 	\Theta_n
\end{array}\right]}.
\end{aligned}
\end{equation}

Next, we introduce the following quadratic cost function to be minimized for the identification of $\Theta_l$, $B_l$, and $\Theta_n$:

\begin{equation} \label{eq:quadratic_cost}
\begin{aligned}
J := \sum_{i = 1}^{N} e_i^\top e_i = \sum_{i = 1}^{N} d_i^\top T^\top T d_i.
\end{aligned}
\end{equation}

In the following sections, for two model classes, we provide stability constraints and formulate the constrained supervised learning of uncertainty models as optimization problems.

§.§.§ Locally Lipschitz Model Class

We first consider locally Lipschitz nonlinearities in (<ref>), for both the basis functions $h(\cdot)$ and the known vector field $g(\cdot)$. Note that given bounded estimated state and input trajectories, we can always find ellipsoids that bound these trajectories. Hereafter, the system state and input trajectories are referred to as state and input sets, respectively.

Note that the set of all trajectories that can be generated by system in (<ref>) is not known a priori. We only have finite data realizations of the estimated state in response to some input trajectories. Given this, we seek models (<ref>) that generate state trajectories “close" to the available data set in response to inputs “close" to the known input trajectories. To this end, we embed estimated state data and input trajectories in some known ellipsoidal sets, $\mathcal{E}_{sys}$ and $\mathcal{E}_{u}$, respectively. This embedding can be obtained efficiently, e.g., exploiting results in <cit.>.

Having these ellipsoids $\mathcal{E}_{sys}$ and $\mathcal{E}_{u}$, to induce closeness of trajectories between system training data and model trajectories, we seek to enforce during learning that all trajectories generated by the model, in response to all input trajectories contained in $\mathcal{E}_{u}$, are contained in some known ellipsoidal set $\mathcal{E}_{inv}$. If we enforce the latter, and $\mathcal{E}_{inv} \subseteq \mathcal{E}_{sys}$, then we can guarantee that the model produces trajectories close to the training date set (close in the sense of set inclusion within $\mathcal{E}_{sys}$), see Figure <ref>.

Let the ellipsoidal set containing the state estimates training data be of the form:

\begin{equation} \label{eq:system_set}
\mathcal{E}_{sys} := \left\{x_s \mid x_s^{\top} F x_s \leq 1\right\},
\end{equation}

with $F \succeq 0$ of appropriate dimensions, and the input set of the form:

\begin{equation} \label{eq:input_set}
\mathcal{E}_{u} := \left\{u \mid u^{\top} U u \leq 1\right\},
\end{equation}

with $U \succeq 0$ of appropriate dimension.

Illustration of different (ellipsoidal) sets for the system state $\mathcal{E}_{sys}$, input $\mathcal{E}_{u}$, and model state $\mathcal{E}_{inv}$.

Because we seek to restrict model trajectories to the sets $\mathcal{E}_{u}$ and $\mathcal{E}_{sys}$, and consider models in (<ref>) with Lipschitz nonlinearities, we only require them to be Lipschitz within these sets. We formally state this in the following assumption.

(Locally Lipschitz Nonlinearities) The functions $g(V_g x, u)$ and $h(V_\eta x, u)$ in (<ref>) are locally Lipschitz in $\mathcal{E}_{u}$ and $\mathcal{E}_{sys}$, i.e., there exist known positive constants $l_{g_x}, l_{g_u}, l_{h_x}$ and $l_{h_u}$ satisfying

\begin{equation} \label{eq:locally_lipschitz}
\begin{aligned}
\|g(V_g x_1, u_1)-g(V_g x_2, u_2)\| &\leq  l_{g_x} \|V_g (x_1-x_2)\|+  l_{g_u}   \|u_1-u_2\|\,\\[2 mm]
\|h(V_\eta x_1, u_1)-h(V_\eta x_2, u_2)\| &\leq  l_{h_x} \|V_\eta (x_1-x_2)\|\
+  l_{h_u}   \|u_1-u_2\|\\
\end{aligned}
\end{equation}

for all $x_1,x_2 \in  \mathcal{E}_{sys} \subset \mathbb{R}^{n}$ and all ${u_1, u_2} \in \mathcal{E}_{u} \subset \mathbb{R}^{l}$.

In what follows, we formulate conditions to guarantee that model trajectories of (<ref>), with locally Lipschitz nonlinearities, belong, in forward time, to an ellipsoidal invariant set:

\begin{equation} \label{eq:model_set}
\mathcal{E}_{inv} := \left\{x \mid x^{\top} P x \leq 1\right\},
\end{equation}

guaranteeing $\mathcal{E}_{inv} \subseteq \mathcal{E}_{sys}$, for some positive definite matrix $P$. Conditions to ensure the latter can be formulated through Lyapunov-based stability tools and the S-procedure <cit.>. For the sake of readability, these conditions are derived in Appendix <ref>. There, it is shown that if the following conditions hold, the ellipsoid $\mathcal{E}_{inv}$ is a forward invariant set for (<ref>) and $\mathcal{E}_{inv} \subseteq \mathcal{E}_{sys}$:

\begin{align}
\Delta+\beta P  & P (B_u+S_{\eta_l} B_l)  & 0 \\
\star & (l_{g_u} + \bar l_{h_u}) I - \alpha U & 0 \\
\star & \star & (\alpha - \beta) I
\end{array}\right] \preceq 0, \label{eq:delta_locally_condition} \\
& \left[\begin{array}{cc}
\gamma     F -  P & 0 \\
* & - \gamma +1
\end{array}\right] \preceq 0, \label{eq:elipsoid_subset} \\
&l_{h_x} \|  \Theta_n  \|  \leq  \bar l_{h_x},   \label{eq:theta_n_bound}
\end{align}

where the involved components in the above equation are defined as

\begin{equation} \label{eq:delta_globally}
\begin{aligned}
\Delta &:= A^\top P+P A+  V_\eta ^\top \Theta_l^\top S_{\eta_l}^\top  P+ P S_{\eta_l} \Theta_l V_\eta + (l_{g_x} +l_{g_u}) P S_g S_g^\top P + (\bar l_{h_x} + \bar l_{h_u}) P S_{\eta_l} S_{\eta_l}^\top P \\
&+l_{g_x} V_g^\top V_g + \bar l_{h_x} V_\eta^\top V_\eta,
\end{aligned}
\end{equation}

\begin{equation} \label{eq:l_h_bar_u}
\begin{aligned}
\bar l_{h_u} &:= \bar l_{h_x} \frac{l_{h_u}}{l_{h_x}},
\end{aligned}
\end{equation}

and $\bar l_{h_x}, \alpha, \beta, \gamma$ are adjustable parameters. In fact, $\bar{l}_{h_x}$ characterizes the size of a set within which $\Theta_n$ is enforced to reside.

Note that the condition in (<ref>) is not convex in the invariant set shape matrix $P$ and uncertainty model parameters $\theta = (\Theta_l, B_l, \Theta_n)$, resulting in a non-convex optimization problem for uncertainty learning with invariance guarantees. Now, we can state the non-convex optimization problem for the locally Lipschitz model class, aiming to find a tractable convex solution later.

(Locally Lipschitz Model Class)
Consider a given data-set $D$ of input, estimated uncertainty, and state realizations and let Assumption <ref> be satisfied. Further consider given ellipsoidal sets $\mathcal{E}_{sys}$ and $\mathcal{E}_{u}$, as introduced in (<ref>) and (<ref>), respectively, containing the available input and state estimates data sets. Find the optimal parameters $\theta = (\Theta_l, B_l, \Theta_n)$ of the uncertainty model $\eta_l(\cdot)$ in (<ref>) (with locally Lipschitz nonlinearities) that minimizes the cost function $J$ in (<ref>) constrained to the ellipsoid $\mathcal{E}_{inv}$ in (<ref>) being a forward invariant set for the extended system model in (<ref>). In other words, solve the non-convex optimization problem:
P,Θ_l, B_l,Θ_n, α, γ

(<ref>)   and   P≻0.

In what follows, we formulate an analogue problem for the globally Lipschitz model class (that is assuming that the nonlinearities in (<ref>) are globally Lipschitz).

(Comparison of Locally and Globally Lipschitz Model Classes)
The main drawback of the globally Lipschitz model class is that it covers a smaller class of nonlinear systems compared to the locally Lipschitz class. However, we can enforce a stronger stability property (input-to-state stability [24]) for globally Lipschitz nonlinearities. Notably, for this model class, knowledge of the system state and input sets is no longer required. Furthermore, in the derivation of stability conditions for globally Lipschitz case, the need for a sufficient approximation (see S-procedure tools in Appendix <ref>) is eliminated. This elimination relaxes the conservatism typically introduced by the sufficient condition. Moreover, the absence of this condition reduces the number of tuning parameters required for optimizing the globally Lipschitz class. The above arguments motivate to consider both model classes.

§.§.§ Globally Lipschitz Model Class

For this model class, the nonlinearities in (<ref>) are assumed globally Lipschitz.

(Globally Lipschitz Nonlinearities) The functions $g(V_g x, u)$ and $h(V_\eta x, u)$ in (<ref>) are globally Lipschitz, i.e., there exist known positive constants $l_{g_x}, l_{g_u}, l_{h_x}$ and $l_{h_u}$ satisfying

\begin{equation} \label{eq:globally_lipschitz}
\begin{aligned}
\|g(V_g x_1, u_1)-g(V_g x_2, u_2)\| &\leq  l_{g_x} \|V_g (x_1-x_2)\|
+  l_{g_u}   \|u_1-u_2\|,\\[2 mm]
\|h(V_\eta x_1, u_1)-h(V_\eta x_2, u_2)\| &\leq l_{h_x} \|V_\eta (x_1-x_2)\|
+ l_{h_u}   \|u_1-u_2\|,\\
\end{aligned}
\end{equation}

for all $x_1,x_2 \in {\mathbb{R}^{n}}$ and all ${u_1, u_2} \in {\mathbb{R}^{l}}$.

In what follows, we formulate conditions to ensure Input-to-State Stability (ISS) of the extended model in (<ref>) satisfying Assumption <ref>. In the following definition, we introduce the notion of ISS for the extended model [24].

(Input-to-State Stability) The extended model (<ref>) is said to be ISS with respect to the input $u(t)$ if there exist a class $\mathcal{K} \mathcal{L}$ function $\bar \beta(\cdot)$ and a class $\mathcal{K}$ function $\bar \alpha(\cdot)$ such that for any initial state $x(t_0)$ and any bounded input $u(t)$, the solution $x(t)$ of (<ref>) exists for all finite $t \geq t_{0}$ and satisfies

\begin{equation}\label{ISS_def}
\|x(t)\| \leq \bar \beta\left(\left\|x\left(t_{0}\right)\right\|, t-t_{0}\right) + \bar \alpha( \sup _{t_0 \leq \tau \leq t}\|u (\tau)\| ).
\end{equation}

For readability, ISS conditions for the model in (<ref>) satisfying Assumption <ref> are derived in Appendix <ref>. It is shown that if the condition (<ref>) together with the following condition hold,

\begin{equation} \label{eq:delta_condition}
\begin{aligned}
&\Delta \prec 0,  \\
\end{aligned}
\end{equation}

where $\Delta$ is as defined in (<ref>), then, the extended model in (<ref>) with globally Lipschitz nonlinearities is ISS with respect to the input $u(t)$.

Now, we can state the non-convex optimization problem for the globally Lipschitz model class, aiming to find a tractable convex approximation later.

(Globally Lipschitz Model Class) Consider a given data-set $D$ of input, estimated uncertainty, and state realizations and let Assumption <ref> be satisfied. Find the optimal parameters $\theta = (\Theta_l, B_l, \Theta_n)$ of the uncertainty model $\eta_l(\cdot)$ in (<ref>) (with globally Lipschitz nonlinearities) that minimizes the cost function $J$ in (<ref>), such that the extended system model in (<ref>) is ISS with respect to input $u(t)$. In other words, solve the non-convex optimization problem:
P,Θ_l, B_l,Θ_n

(<ref>), (<ref>), P≻0.

The challenge in Problem Settings <ref> and <ref> is that stability constraints are not convex, resulting in non-convex optimization problems. To tackle this challenge, we provide two approximate convex solutions for both problem settings using two distinct approaches: first, cost modification approach and second, constraint modification approach. Moreover, as an alternative solution, we provide a procedure to solve the original non-convex problems via sequential convex programming.

§ COST MODIFICATION APPROACH

We convexify the stability constraints by a change of variables and rewrite the cost function in (<ref>) in them. First, we provide the approximate solution for Problem Setting <ref>, followed by the approximate solution for Problem Setting <ref>.

§.§ Locally Lipschitz Model Class

The following theorem formalizes the associated convex optimization problem obtained via the cost modification approach as an approximation to the non-convex optimization problem in (<ref>).

(Stable Locally Lipschitz Model Learning with Modified Cost)
Consider system (<ref>), a given data-set $D$ of input, estimated uncertainty, and state realizations. In addition, consider given ellipsoidal sets $\mathcal{E}_{sys}$ and $\mathcal{E}_{u}$, as introduced in (<ref>) and (<ref>) with shape matrices $F \succeq 0$ and $U \succeq 0$, respectively. Further consider the extended system model in (<ref>), under Assumption <ref> with Lipschitz constants $l_{g_u}, l_{g_x}, l_{h_u}$ and $l_{h_x}$. Consider the following convex program:
\beginmini!|s|[2]
P,S, R, Z, W, α, γ

M_11     M_12      0     M_14     M_15

⋆    M_22     0     0     0

⋆    ⋆     M_33    0     0

⋆    ⋆    ⋆    -I     0

⋆    ⋆    ⋆    ⋆    -I
l̅_h_x  I     l_h_x Z^⊤
l̅_h_x  (2  μ_1  P - μ_1^2 I)
]	≽0
γF -  P     0
- γ+1
2 μ_2 P     T̃D̃^⊤    μ_2 I
I      0
*     W
P≻0,   α, γ≥0

\endmini!
with the involved matrices defined as follows:

\begin{equation} \label{eq:delta_lmi_components} \tag{16d}
\begin{aligned}
M_{1} &:= A^\top P+P A+  V_\eta^\top S^\top + S V_\eta  +l_{g_x} V_g^\top V_g + \bar l_{h_x} V_\eta^\top V_\eta,  &M_{11} &:=  M_{1} + \beta P,\\
M_{12} &:=	P B_u+R,
&M_{22} &:= (l_{g_u} + \bar l_{h_u}) I - \alpha U, \\
M_{14} &:=\sqrt{l_{g_x} +l_{g_u}}  P S_g,
&M_{15} &:=\sqrt{\bar l_{h_x} + \bar l_{h_u}} P, \\
M_{33} &:=		 (\alpha - \beta) I,
&\tilde{T} &:= {\left[\begin{array}{cccc}
S  & R & -P & Z
\end{array}\right]},
\end{aligned}
\end{equation}

for given positive scalars $\bar l_{h_x}, \beta, \mu_1, \mu_2$, and $\bar l_{h_u}$ as defined in (<ref>), where $\tilde{D}$ is the Cholesky decomposition of the data matrix $D$ in (<ref>) (i.e., $D = \tilde{D}^\top \tilde{D}$), and the remaining matrices are the known parts of the system dynamics in (<ref>). Denote part of the optimizers of (<ref>) as $P^\star$, $S^\star$, $R^\star,Z^\star$, and $W^\star$. Then, the following parameters of the extended model (<ref>), $S_{\eta_l} = I, \Theta_l = \Theta_l^\star = P^{\star^{-1}} S^\star ,B_l = B_l^\star = P^{\star^{-1}} R^\star$ and $\Theta_n = \Theta_n^\star = P^{\star^{-1}} Z^\star$ guarantee that the ellipsoid $\mathcal{E}_{inv}$ in (<ref>) with $P = P^{\star}$ is a forward invariant set for the extended system model in (<ref>) and that $\mathcal{E}_{inv} \subseteq \mathcal{E}_{sys}$. In addition, it holds that the cost $J$ of (<ref>) satisfies $J \leq tr(W^\star)$. As such, (<ref>) represents an approximate convexified version of the problem in (<ref>).

The proof can be found in Appendix <ref>.

The scalar parameters $\bar l_{h_x}, \beta, \mu_1, \mu_2$, and $\bar l_{h_u}$  in Theorem <ref> are tuned for the minimal feasible cost $tr(W)$ by a line search.

§.§ Globally Lipschitz Model Class

The subsequent theorem formalizes the convex optimization problem using the cost modification approach, serving as an approximation to the non-convex optimization problem in (<ref>).
(Stable Globally Lipschitz Model Learning with Modified Cost)
Consider system (<ref>), a given data-set $D$ of input, estimated uncertainty, and state realizations. In addition, consider the extended system model in (<ref>), under Assumption <ref> with Lipschitz constants $l_{g_u}, l_{g_x}, l_{h_u}$ and $l_{h_x}$. Consider the following convex program:
\beginmini!|s|[2]
P,S, R, Z, W

M_1      M_14     M_15

⋆    -I     0

⋆    ⋆    -I
l̅_h_x  I     l_h_x Z^⊤
l̅_h_x  (2  μ_1  P - μ_1^2 I)
2 μ_2 P     T̃D̃^⊤    μ_2 I
I      0
*     W
]≽0,   P ≻0

\endmini!
with $M_{1}, M_{14}, M_{15}, \tilde{T}$ as defined in (<ref>), for given positive scalars $\bar l_{h_x}, \mu_1,$ and $\mu_2$, where $\tilde{D}$ is the Cholesky decomposition of the data matrix $D$ in (<ref>), and the remaining matrices are the known parts of the system dynamics in (<ref>). Denote part of the optimizers of (<ref>) as $P^\star$, $S^\star$, $R^\star,Z^\star$, and $W^\star$. Then, the following parameters of the extended model (<ref>), $S_{\eta_l} = I, \Theta_l = \Theta_l^\star = P^{\star^{-1}} S^\star ,B_l = B_l^\star = P^{\star^{-1}} R^\star$ and $\Theta_n = \Theta_n^\star = P^{\star^{-1}} Z^\star$ guarantee that the extended model in (<ref>) is ISS with respect to input $u(t)$. In addition, it holds that the cost $J$ of (<ref>) satisfies $J \leq tr(W^\star)$. As such, (<ref>) represents an approximate convexified version of the problem in (<ref>). As a special case, when the Lipschitz constants are zero (no nonlinearity) the conditions in (<ref>) and (<ref>) reduces to

\begin{equation} \label{eq:sdp_quadratic_cost_globally_linear} \tag{17c}
A^\top P+P A+  V_\eta^\top S^\top + S V_\eta \prec 0.
\end{equation}

The proof follows the line of reasoning of the proof of Theorem <ref> and it is omitted here due to space constraints.

(Surrogate Convex Optimizations with Modified Cost)
We remark that the semi-definite programs in (<ref>) and (<ref>) are not equivalent to the non-convex optimization problems in (<ref>) and (<ref>), respectively (i.e., these are convex approximations). This approximation is due to setting $S_{\eta_l} = I$ and using two sufficient conditions (two lower bounds) in the derivation of LMI conditions in each of the convex programs (see details at the end of the proof of Theorem <ref>). Although by letting $S_{\eta_l} = I$, we do not use the known structure of uncertainty, this makes the problem tractable. Note that here, we do not use knowledge of uncertainty structure.

Next, we follow a different approach to formulate alternative surrogate (approximate) convex optimization problems for non-convex optimization problems in Problem Settings <ref> and <ref>.

§ CONSTRAINT MODIFICATION APPROACH

In this section, instead of changing the model-related variables in the stability constraints, we formulate sufficient conditions (upper bounds) for the stability constraints which are linear in all the optimization parameters. This convexifies the optimization problems in (<ref>) and (<ref>). Similar to previous section, first, we provide the approximate solution for Problem Setting <ref>, followed by the approximate solution for Problem Setting <ref>.

§.§ Locally Lipschitz Model Class

The following theorem provides the convex approximation of the non-convex optimization problem in Problem Setting <ref> using the constraint modification approach.

(Stable Locally Lipschitz Model Learning with Modified Constraint)
Consider system (<ref>), a given data-set $D$ of input, estimated uncertainty, and state realizations. In addition, consider given ellipsoidal sets $\mathcal{E}_{sys}$ and $\mathcal{E}_{u}$, as introduced in (<ref>) and (<ref>) with characteristic matrices $F \succeq 0$ and $U \succeq 0$, respectively, and the extended system model in (<ref>), under Assumption <ref> with Lipschitz constants $l_{g_u}, l_{g_x}, l_{h_u}$ and $l_{h_x}$. Consider the following convex program:
\beginmini!|s|[2]
Q,Θ_l, B_l, Θ_n, W, α, γ

N_11      N_12     0    N_14      N_15     N_16

⋆    N_22     0     0     0     0

⋆     ⋆    N_33     0     0     0

⋆    ⋆     ⋆     -2 γ̅I     0     0

⋆    ⋆     ⋆    ⋆     -I     0

⋆    ⋆     ⋆     ⋆     ⋆    -I
l̅_h_x I      l_h_x Θ_n^⊤

⋆    l̅_h_x I
] ≽0
γF -  2 μ_3 I + μ_3^2 Q     0
- γ+1
tr(T D T^⊤)≤tr(W),   Q≻0,   α, γ≥0

\endmini!
with the involved matrices defined as follows:

\begin{equation} \label{eq:delta_lmi_components_constraint_app} \tag{18e}
\begin{aligned}
N_{1} &:=  A Q +Q A^\top+ (l_{g_x} +l_{g_u}) S_g S_g^\top + (\bar l_{h_x}+ \bar l_{h_u}) S_{\eta_l} S_{\eta_l}^\top, &N_{11}& :=  N_{1} + \beta Q,\\
N_{12} &:=	B_u+S_{\eta} B_l,
&N_{14} &:= S_\eta \Theta_l V_\eta +\bar \gamma Q, \\
N_{15} &:= \sqrt{l_{g_x}} Q V_g^\top,
&N_{16} &:= \sqrt{\bar l_{h_x}} Q V_\eta^\top, \\
N_{22} &:= (l_{g_u} + \bar l_{h_u}) I - \alpha U,
&N_{33} &:=		 (\alpha - \beta) I,
\end{aligned}
\end{equation}

given positive scalars $\bar l_{h_x}, \beta, \mu_3, \bar \gamma$, and $\bar l_{h_u}$ as defined in (<ref>). The matrices $D, T$ are defined in (<ref>) and (<ref>), respectively, and the remaining matrices are the known parts of the system dynamics in (<ref>). Denote part of the optimizers of (<ref>) as $Q^\star$, $\Theta_l^\star$, $B_l^\star$, $\Theta_n^\star$, and $W^\star$. Then, the following parameters of the extended model (<ref>), $S_{\eta_l} = S_\eta, \Theta_l = \Theta_l^\star,B_l = B_l^\star$ and $\Theta_n = \Theta_n^\star$ guarantee that the ellipsoid $\mathcal{E}_{inv}$ as defined in (<ref>) with $P = Q^{{\star}^{-1}}$ is a forward invariant set for the extended system model in (<ref>). In addition, it holds that the cost $J$ of (<ref>) satisfies $J \leq tr(W^\star)$. As such (<ref>) represents an approximate convexified version of the problem in (<ref>).

The proof can be found in Appendix <ref>.

§.§ Globally Lipschitz Model Class

The subsequent theorem presents a formalization of the convex optimization problem derived through the constraint modification approach, serving as an approximation to the non-convex optimization problem in (<ref>).

(Stable Globally Lipschitz Model Learning with Modified Constraint)
Consider system (<ref>), a given data-set $D$ of input, estimated uncertainty, and state realizations. In addition, consider the extended system model in (<ref>), under Assumption <ref> with Lipschitz constants $l_{g_u}, l_{g_x}, l_{h_u}$ and $l_{h_x}$. Consider the following convex program:
\beginmini!|s|[2]
Q, Θ_l, B_l, Θ_n, W

N_1      N_14     N_15     N_16

⋆    -2 γ̅I     0      0

⋆    ⋆    -I     0

⋆    ⋆     ⋆    -I
l̅_h_x I      l_h_x Θ_n^⊤

⋆    l̅_h_x I
tr(T D T^⊤)≤tr(W),  Q≻0

\endmini!
with $N_{1}, N_{14}, N_{15}, N_{16}$ as defined in (<ref>), given positive scalars $\bar l_{h_x}$ and $\bar \gamma$. The matrices $D, T$ are defined in (<ref>) and (<ref>), respectively, and the remaining matrices are the known parts of the system dynamics in (<ref>). Denote part of the optimizers of (<ref>) as $Q^\star$, $\Theta_l^\star$, $B_l^\star$, $\Theta_n^\star$, and $W^\star$. Then, the following parameters of the extended model (<ref>), $S_{\eta_l} = S_\eta, \Theta_l = \Theta_l^\star,B_l = B_l^\star$ and $\Theta_n = \Theta_n^\star$ guarantee that the extended model in (<ref>) is ISS with respect to input $u(t)$. In addition, it holds that the cost $J$ of (<ref>) satisfies $J \leq tr(W^\star)$. As such (<ref>) represents an approximate convexified problem of the problem in (<ref>). As a special case, when the Lipschitz constants are zero (no nonlinearity) the conditions in (<ref>) and (<ref>) reduces to

\begin{equation} \label{eq:sdp_quadratic_constraint_globally_linear}
\left[\begin{array}{cccc} \tag{19c}
A Q +Q A^\top   & S_\eta \Theta_l +\bar \gamma Q\\
\star & -2 \bar \gamma I\\
\end{array}\right] \prec 0.
\end{equation}

The proof follows the line of reasoning of the proof of Theorem <ref> and is omitted for the sake of brevity.

(Surrogate Convex Optimizations with Modified Constraint)
The semi-definite programs in (<ref>) and (<ref>) are convex approximations of the non-convex optimization problems in (<ref>) and (<ref>), respectively. Note that, a disadvantage of the constraint modification approach compared to the cost modification approach is that to ensure the feasibility of the semi-definite problems in (<ref>) and (<ref>), the known $A$ matrix of the system in (<ref>) has to be Hurwitz. Moreover, in deriving a convex approximation for the stability constraint, a non-dissipative term is substituted by a term that can potentially be dissipative (using Young's inequality). Thus, this solution may be more conservative compared to the cost modification approach (see for example, results in Section <ref>).  On the other hand, unlike the cost modification approach, Theorems <ref> and <ref> use the knowledge of uncertainty structure by setting $S_{\eta_l} = S_\eta$, which is potentially beneficial.

In the following section, an alternative solution for solving the non-convex problems in Problem Settings <ref> and <ref> is provided using a procedure based on the technique of sequential convex programming.

§ SEQUENTIAL CONVEX PROGRAM APPROACH

Another approach that can provide solution for the non-convex optimization problems in (<ref>) and (<ref>) is Sequential Convex Programming (SCP). This method typically involves fixing some parameters and solving a series of convex sub-problems iteratively [23]. For the given non-convex optimization problems, we can fix Lyapunov function parameter $P$ or uncertainty model parameters $\theta = (\Theta_l, B_l, \Theta_n)$. The former transforms the non-convex optimization problem(s) in (<ref>) (or in (<ref>)) to feasibility problem(s) and the latter to convex optimization sub-problem(s). Then, this transformed problems can be solved iteratively until the convergence of the cost(s) of the convex optimization sub-problem(s).

Next, an SCP algorithm is presented as a sub-optimal solution for the non-convex optimization problem(s) in (<ref>) (or in (<ref>)). Steps 1 and 2 in this algorithm can be swapped to create a dual algorithm. In the loop below, initialization for $\theta$ is necessary, whereas for the dual loop, initialization for $P$ is required. SCP algorithm consisting of three steps:

* Fix $\theta$ and find a feasible solution for the transformed feasibility problems in (<ref>) (or (<ref>));

* Fix $P$ as the feasible solution from Step 2 and solve the transformed sub-convex problem(s) in (<ref>) (or (<ref>)).

* Redo Steps 1 and 2 until the cost of the transformed sub-convex problem(s) in Step 2 converge. Denote part of the optimizers of the transformed sub-convex problem(s) with $\theta = \theta^{\star}$ that provides the uncertainty model parameters. These uncertainty parameters guarantee stability of the extended model in (<ref>).

(SCP and Convex Approximations Comparison)
The effectiveness of this SCP method is recognized to rely on the initialization [23]. For example, in the simulation results in Section <ref>, no other initialization was found except the solution of the convex approximation optimizations provided in Sections <ref> and <ref>, highlighting the significance of the proposed convex programs. Moreover, depending on the optimization problem and initialization, the SCP method might require many iterations for convergence, making this method expensive in terms of computations. On the other hand, the drawbacks of the convex approximation programs are that they require some tuning and, due to conservatism in their derivation, their solutions might be more sub-optimal (i.e., have a larger cost) than the SCP method.

In the results so far, we assumed that state and uncertainty realizations are available, which is, in practice typically not the case. In what follows, we present a solution for uncertainty and state estimation based on only input and output data.

§ UNCERTAINTY AND STATE ESTIMATION

First, we formulate the uncertainty and state estimation problem before providing a solution for that problem. Consider system in (<ref>) and the required assumptions below to ensure that the problem is well-posed as is common in the existing literature [25, 26, 27].

The following assumptions stand throughout this section:

* State and Input Boundedness:  The state $x_s(t)$ and input $u(t)$ remain bounded over any finite time interval.

* $\mathcal{C}^r$ Uncertainty Vector: The uncertainty function $\eta(V_\eta x_s(t), u(t))$ in (<ref>) is $r$ times differentiable with respect to time, i.e., the total time derivatives $\eta^{(1)}(V_\eta x_s(t),u(t))$,$\eta^{(2)}(V_\eta x_s(t),u(t))$, ..., $\eta^{(r)}(V_\eta x_s(t),u(t))$ exist and are continuous.

* Bounded Disturbances: The disturbance vectors $\omega(t)$ and $\nu(t)$ in (<ref>) are bounded on any finite time interval, and $\nu(t)$ is once differentiable in $t$, i.e., $\dot{\nu}(t)$ exists, and is continuous and bounded over any finite time interval.

The following filter is designed for uncertainty and state estimation:

\begin{equation} 		\label{eq:filter}
\left\{\begin{aligned}
\dot{z} =&f(z,u,y_s;\psi), \\
\hat{\eta} =&\phi_1(z,y_s;\psi), \\
\hat{x}_s =&\phi_2(z,y_s;\psi),
\end{aligned}\right.
\end{equation}

where $z \in {\mathbb{R}^{n_z}}$ is the internal state of the filter with ${n_z} \in \mathbb{N}$. Functions $f: \mathbb{R}^{n_z} \times \mathbb{R}^{l} \times \mathbb{R}^{m} \to \mathbb{R}^{n_z}$, $\phi_{1}: \mathbb{R}^{n_z} \times \mathbb{R}^{m} \to \mathbb{R}^{n_\eta}$, and $\phi_{2}: \mathbb{R}^{n_z} \times \mathbb{R}^{m} \to \mathbb{R}^{n}$ characterize the filter structure, $\psi$ denotes design parameters. Define $\hat{x}_d := (\hat{\eta}, \hat{x}_s)$ (representing the estimate of both the uncertainty and the state) and its estimation error as

\begin{equation} \label{eq:ed}
e_d :=\hat{x}_d - x_d,
\end{equation}

where ${x}_d := ({\eta}, {x}_s)$. Later, it will be shown that, for the uncertainty-state estimator design proposed in Section <ref>, the estimation error dynamics, explicated in Equation (<ref>) in Section <ref>, exhibits $(\omega,  \eta^{(r)},\nu, \dot{\nu})$ as a perturbation input. Having that in mind, we can state the problem we aim to solve at a high abstraction level.

(Uncertainty and State Estimation - Abstract Level) Consider the system (<ref>) with known input and output signals, $u$ and $y_s$, and the uncertainty-state estimator filter (<ref>). For given $r$, design the filter parameters $\psi$ such that the following properties are guaranteed:

1) Stability: The estimation error dynamics is input-to-state stable with respect to the perturbation input $(\omega,  \eta^{(r)},\nu, \dot{\nu})$;

2) Disturbance Attenuation: For $\nu = 0$, the $\mathcal{L}_2$-gain from $(\omega,  \eta^{(r)})$ to $e_d$ is bounded by some known $c_1 > 0$, for $t \geq 0$;

3) Noise Rejection: For $ (\omega,  \eta^{(r)})= 0$, the $\mathcal{L}_2-\mathcal{L}_\infty$ induced gain from $(\nu, \dot{\nu})$ to $e_d$ is bounded by some known $c_2 > 0$, for $t \geq 0$.

In what follows, we present a synthesis procedure to design the filter dynamics so that the three properties stated in Sub-Problem <ref> are realized.

§.§ Ultra Local Uncertainty Representation

First some preliminaries are discussed, which are required to design the estimator filter (<ref>). Considering that the uncertainty $\eta(x_s(t), u(t))$ in (<ref>) is an implicit function of time, for all $x_s(t)$ and $u(t)$, we can write an entry-wise $r$-th order Taylor time-polynomial approximation at time $t$ of $\eta$ as $\bar \eta = a_0 + a_1 t+ \dots + a_{r-1} t^{r-1}$ with coefficients $a_i \in {\mathbb{R}^{n_\eta}}, i= 0, \dots, r-1$. This model can be written in state-space form as

\begin{equation} \label{eq:fault_model}
\left\{\begin{aligned}
\dot{\bar \zeta}_j &= \bar \zeta_{j+1},  \qquad 0 < j < r, \\
\dot{\bar \zeta}_{r} &= 0,\\
\bar \eta &= \bar \zeta_1,
\end{aligned}\right.
\end{equation}

where $\bar \zeta_j \in {\mathbb{R}^{n_\eta}}$. Clearly, in the above model, we have $\bar \eta^{(r)} = {0}$, which might not be true for actual uncertainty signal $\eta$. Under Assumption <ref>, the actual internal state-space representation of $\eta$ is as follows:

\begin{equation} \label{eq:fault_system}
\left\{\begin{aligned}
\dot{\zeta}_j &= \zeta_{j+1},  \qquad 0 < j < r, \\
\dot{\zeta}_{r} &= \eta^{(r)},\\
\eta &= \zeta_1,
\end{aligned}\right.
\end{equation}

where $\zeta_j \in {\mathbb{R}^{n_\eta}}$. Clearly, the accuracy of the approximate model (<ref>) increases as $\eta^{(r)}$ goes to zero (entry-wise). In the following, to design the uncertainty-state estimator we augment the system state, $x_s(t)$, with the states of the actual uncertainty internal state $\zeta_j(t), j \in\left\{1, \ldots, r\right\}$, and augment the system dynamics in (<ref>) with (<ref>). We then design a filter (observer) for the augmented system to simultaneously estimate $x_s$ and $\zeta_j$ using model (<ref>). We remark that proper selection of the number of the uncertainty derivatives, $r$, added to the approximated model (<ref>) (and (<ref>)) is problem-dependent, see [28] for discussion on selection of $r$.

§.§ Augmented Dynamics

Based on the internal representation of the uncertainty in (<ref>), define the augmented state $x_a:= (x_s, \zeta_1,\zeta_2,\ldots,  \zeta_{r})$, and rewrite the augmented dynamics using (<ref>) and (<ref>) as

\begin{equation}		\label{eq:augmented_system}
\begin{aligned}
\left\{\begin{aligned}
\dot{x}_{a} &=A_{a} x_{a}+B_{u_a} u+ S_{g_a} g(V_{g_a} x_{a}, u) +B_{\omega_a} \omega_{a},\\
y_s &= C_{a} x_{a} +  D_{\nu} \nu,
\end{aligned}\right.
\end{aligned}
\end{equation}

where the involved matrices in the above equation are defined as

\begin{equation} 		\label{eq:augmented_matrices}
\begin{aligned}
A_a &:=
\begin{bmatrix}
A & S_\eta & {0} \\
{0} & {0} & I_{d_n} \\
{0} & {0} & {0}
\end{bmatrix}, \quad
B_{u_a} :=	\left[\begin{array}{cc}
B_u \\
\end{array}\right],  \quad
S_{g_a} := \left[\begin{array}{cc}
S_g^\top &
\end{array}\right]^\top, \quad
V_{g_a} := [\begin{array}[]{ll}V_g &0 \end{array}],\\
B_{\omega_a} &:= \left[\begin{array}{cc}
B_{\omega} & 0 \\
0 & 0 \\
0 &  I_{n_\eta}
\end{array}\right], \quad
\omega_{a} :=	\left[\begin{array}{cc}
\omega \\
\eta^{{(r)}}
\end{array}\right], \quad
C_{a} := \left[\begin{array}{ll}	C	&0 \end{array}\right], \quad
V_{\eta_a} := [\begin{array}[]{ll}V_\eta &0 \end{array}]
\end{aligned}
\end{equation}

with $d_n := (r-1)n_{\eta}$.

§.§ Uncertainty-State Estimator

In this section, considering the uncertainty-state estimator general structure in (<ref>), inspired from observer-based approaches, we design $f(\cdot)$ and $\phi_i(\cdot), i=1,2$, as follows:

\begin{equation}
\begin{aligned} 		\label{eq:observer_dynamics}
f(z,u,y;\psi) &= N z+G u_a+L y_s +M S_{g_a} {g}(V_{g_a}\hat{x}_a+ H\color{black}(y_s - C_{a} \hat{x}_a),u), \\
\phi_i (z,y;\psi) =& \bar{C}_i (z-E y_s), \\
\end{aligned}
\end{equation}

with $\hat{x}_a = z-E y_s$, filter state $z \in {\mathbb{R}^{n_z}}$, $n_z = n+r n_{\eta}$,

\begin{equation*}
\bar{C}_1 := {\left[\begin{array}{ccc}
{0} & I_{n_\eta} & {0}
\end{array}\right]}, \quad
\bar{C}_2 := {\left[\begin{array}{cc}
I_n & {0}
\end{array}\right]},
\end{equation*}

and matrices $(N,G,L)$ defined as

\begin{equation}		\label{eq:observer_matrices}
\begin{aligned}
N &:=M A_a-K C_a,
&M&:=I+E C_a, \\[1mm]
G&:=M B_a,
&L&:=K(I+C_a E)-M A_a E.
\end{aligned}
\end{equation}

Matrices $E.K$ and $H$ are filter gains to be designed which can be collected as $\psi = (E,K,H)$. Note that according to (<ref>), the part of the augmented state, $x_a$, that we use to reconstruct uncertainty and state signals is $\bar{C}_a x_{a}$ with

\begin{equation} 	\label{eq:c_bar_a}
\bar{C}_a := {\left[\begin{array}{cc}
\bar{C}_1^\top &
\bar{C}_2^\top
\end{array}\right]^\top}.
\end{equation}

In the following section, we analyze the estimator error dynamics.

§.§ Uncertainty-State Estimator Error Dynamics

Consider the augmented state estimate $\hat{x}_a$ and let us define estimation error as

\begin{equation*}
e:=\hat{x}_a-x_a=z-x_a-E y_s=z-M x_a-E D_\nu \nu.
\end{equation*}

Then, given the algebraic relations in (<ref>), the estimation error dynamics can be written as

\begin{equation}	\label{eq:error_daymics_only}
\left\{\begin{aligned}
\dot{e} &= N e + M S_{g_a} \delta g - M B_{\omega_a} \omega_{a} + B_{\nu_a} \nu_a, \\
e_d &= \bar{C}_a e,
\end{aligned}\right.
\end{equation}

where the involved components in the above equation are defined as

\begin{equation} \label{eq:error_dyamincs_related_matrices}
\begin{aligned}
\delta g &:= g(V_{g_a}\hat{x}_{a} + H (y_s - C_{a} \hat{x}_{a}), u)-g(V_{g_a} {x}_{a}, u), \\
B_{\nu_a} &:= {\left[\begin{array}{cc}
KD_\nu  & -E D_\nu
\end{array}\right]}, \qquad \nu_a := {\left[\begin{array}{cc}
\nu^\top(t) & \dot{\nu}^\top(t)
\end{array}\right]^\top}.
\end{aligned}
\end{equation}

Now, we can restate Problem <ref> in a more formal way.

(Uncertainty and State Estimation) Consider the system (<ref>) with known input and output signals, $u$ and $y_s$. Furthermore, consider the internal uncertainty dynamics (<ref>), its Taylor approximation (<ref>), the augmented dynamics (<ref>), the uncertainty-state estimator (<ref>) with functions defined in (<ref>), and let Assumption <ref> be satisfied. Design the estimator gain matrices $\psi = (E,K, H)$ such that the following properties are guaranteed:

1) Stability: There exist a class $\mathcal{K} \mathcal{L}$ function $\bar \beta(\cdot)$ and a class $\mathcal{K}$ function $\bar \alpha(\cdot)$ such that for any initial estimation error $e(t_0)$ and any bounded input $\bar \omega_{a} :=(\omega , \eta^{(r)}, \nu, \dot{\nu})$, the solution $e(t)$ of (<ref>) exists for all $t \geq t_{0}$ and satisfies

\begin{equation*}\label{ISS_def}
\|e(t)\| \leq \bar \beta\left(\left\|e\left(t_{0}\right)\right\|, t-t_{0}\right) + \bar \alpha( \sup _{t_0 \leq \tau \leq t}\|\bar \omega_a (\tau)\| )
\end{equation*}

2) Disturbance Attenuation: for $\nu=\dot{\nu}=0$, it holds that

\begin{equation} \label{eq:J1}
\tilde{J}_1(\psi) := \sup_{(\omega , \eta)} \frac{\|e_d\|_{\mathcal{L}_2}}{\|(\omega , \eta^{(r)})\|_{\mathcal{L}_2}}
\end{equation}

is bounded by some known $c_1 > 0$;

3) Noise Rejection: for $\omega=\eta^{(r)}=0$, it holds that

\begin{equation} \label{eq:J2}
\tilde{J}_2(\psi) := \sup_{\nu} \frac{\|e_d\|_{\mathcal{L}_2}}{\|(\nu, \dot{\nu})\|_{\mathcal{L}_\infty}}
\end{equation}

is bounded by some known $c_2 > 0$.

The essence of Problem <ref> is to find an uncertainty-state estimator that, firstly,  ensures a bounded estimation error, $e(t)$, for any input (input to state boundedness); when $\nu = 0$ the $\mathcal{L}_2$-gain of the mapping from $(\omega , \eta^{(r)})$ to $e_d$ (the desired estimation error) is upper bounded by some $c_1>0$; and when $\omega=\eta^{(r)} = 0$, the $\mathcal{L}_2-\mathcal{L}_\infty$ induced gain from $(\nu, \dot{\nu})$ to $e_d$ is upper bounded by some $c_2>0$; and, for $\omega_{a} = {0}$ and $\nu_{a} = {0}$, $e(t)$ goes to zero asymptotically (internal stability).

§.§ Uncertainty-State Estimator Design

In the following proposition, we provide the solution of Problem <ref> as a semi-definite problem, where we seek to minimize the $\mathcal{L}_2$-gain from $\omega_a$ to $e_d$ for an acceptable upper bound on the $\mathcal{L}_2-\mathcal{L}_\infty$ induced norm from $\nu_a$ to $e_d$ (there exist a trade-off between these two norms, see [22]). Moreover, we add the Input-to-State Stability (ISS) constraint with respect to filter error dynamics input $(\omega_a, \nu_a)$ to this program to enforce that stability of the resulting estimation filter.

(Estimator Design [22])
Consider the system (<ref>), the augmented dynamics (<ref>), the uncertainty-state estimator (<ref>) with $f(\cdot)$ and $\phi(\cdot)$ as defined in (<ref>), and the corresponding estimation error dynamics (<ref>), under Assumption <ref> with Lipschitz constant $l_{g_x}$. Further, consider the following convex program:

\begin{equation}
\begin{array}{cl}
\min \limits_{\Pi, \bar R, \bar Q, H, \rho, \sigma}  & \rho \\
\text{s.t.} & 						\vspace{1 mm}
\left[\begin{array}{cc}
X_{11} & X_{12} \\
* & -I
\end{array}\right] \prec 0,\\
& \left[\begin{array}{ccc}
L_{11} & -(\Pi + \bar R C_{a}) B_{\omega_a} & X_{12} \\
* & -{\rho} a I  & 0 \\
* & * & -I
\end{array}\right]\preceq 0, \\[7mm]
X_{11} & H_{12} & 0 & X_{12}\\
* & -b^2 I  & T_\nu^\top H^\top & 0\\
* & * & -I & 0\\
* & * & * & -I
\end{array}\right]\preceq0, \\
\Pi	 & \bar{C}_a^\top \\
* & \sigma I
\end{array}\right] \succeq 0, \quad \Pi \succ 0, \quad \sigma \leq \sigma_{max}\\
\end{array}
\label{eq:mimization}
\end{equation}

with the involved matrices defined as follows:

\begin{equation*}
\begin{aligned}
X_{11} &:= A_{a}^{\top}  \Pi +A_{a}^{\top} C_{a}^{\top} \bar R^{\top}-C_{a}^\top \bar Q^\top + \Pi A_{a} +\bar R C_{a} A_{a}  - \bar Q C_{a} + l_{g_x} (V_{g_a}^\top V_{g_a} -V_{g_a}^\top  H  C_{a} - C_{a}^\top  H^\top  V_{g_a}), \\
X_{12} &:= {\left[\begin{array}{c}
\sqrt{2 l_{g_x}}\big((\Pi+\bar R C_{a})S_{g_a}\big)^\top  \\  \sqrt{l_{g_x}}(H C_{a})^\top
\end{array}\right]^\top}, \\
L_{11} &:=	X_{11}+ a \bar{C}_a^\top \bar{C}_a,\quad
H_{12} := [\begin{array}[]{ll}  \bar Q D_\nu  & - \bar R D_\nu \end{array}], \quad   T_\nu  := [\begin{array}[]{ll}  D_\nu & 0 \end{array}],
\end{aligned}
\end{equation*}

given scalars $a, \sigma_{max} > 0$, $b$, matrix $\bar{C}_a$ as in (<ref>), and the remaining matrices in (<ref>) and (<ref>). Denote the optimizers as $\Pi^\star$, $\bar R^\star$, $\bar Q^\star$, $H^\star$, $\rho^\star$ and $\sigma^\star$. Then, the following parameters of (<ref>), $\psi= \psi^\star = \{E^\star = \Pi^{\star^{-1}} \bar R^\star ,K^\star = \Pi^{\star^{-1}} \bar Q^\star,H^\star\}$ guarantees the following:

* The estimation error dynamics in (<ref>) is ISS with respect to perturbation input $(\omega_a, \nu_a)$;

* $\tilde{J}_1(\cdot)$ in (<ref>) is upper bounded by $\sqrt{ \rho^\star}$, i.e., the $\mathcal{L}_2$-gain of (<ref>) with $\nu = 0$ from $\omega_a = (\omega , \eta^{(r)})$ to $e_d$ is upper bounded by $\sqrt{ \rho^\star}$;

* $\tilde{J}_2(\cdot)$ in (<ref>) is upper bounded by $\sqrt{b \sigma^\star}$, i.e., the $\mathcal{L}_2-\mathcal{L}_\infty$ induced norm of (<ref>) with $\omega_a= 0$ from $\nu_a = (\nu, \dot{\nu})$ to $e_d$ is upper bounded by $\sqrt{b \sigma^\star}$.

Note that the solution we have provided for uncertainty-state estimation is applicable to globally Lipschitz nonlinear systems (i.e., Assumption <ref> is satisfied). However, for locally Lipschitz nonlinear systems, a similar unknown input-state estimator proposed in [29] can be used.

§ CASE STUDY

In this section, we evaluate the proposed methods via a case study. We assess the obtained model accuracy for the considered class of nonlinear systems in the sense of output accuracy. The case study is a four-degree-of-freedom roll plane maneuver of a vehicle with nonlinear suspension, as considered in [30, 31]. See Figure <ref> for a schematic of this system. The equations of motion and the system parameters are given in Appendix <ref>. The nonlinearity of the system is rooted from nonlinear stiffness and damping between the roll bar and both tires. As available measured outputs, we have the relative displacement and relative velocity between the roll bar and both tires (i.e.,  $y_{s_1} = q_1 - q_3, y_{s_2} = q_2 - q_4, y_{s_3} = \dot{q}_1 - \dot{q}_3, y_{s_4} = \dot{q}_2 - \dot{q}_4$, see Figure <ref>). The system inputs are the position of the tires ($u_1$ and $u_2$ in Figure <ref>). The training and test data sets are generated by providing different system inputs, as described in Section <ref>.

Roll plane system schematic.

For the above-mentioned system, we have (unknown) parametric uncertainty in the linear stiffness between the roll bar and tires (i.e., uncertainty in the linear parameters $\tilde{k}_1, \tilde{k}_2$ shown in Figure <ref>), and the nonlinear stiffnesses (i.e., $\tilde{k}_{1_n}, \tilde{k}_{2_n}$ shown in Figure <ref>) are completely unknown. The problem is to find a data-based model for the mentioned uncertainties and guarantee stability of the extended model (i.e., known model together with uncertainty model). The proposed methods (three algorithms) presume that a (training) data-set of inputs, estimated uncertainties, and estimated state realizations is given (using uncertainty and state estimators discussed in Section <ref>). The results of uncertainty and state estimation for one training data-set are provided in Appendix <ref>.

Following the proposed method for this problem, since the nonlinear stiffness is unknown, first we assume that we have basis functions as given by $h(\cdot)$ in (<ref>) that are representative of the unknown nonlinearity (the proposed method is evaluated for a range of basis functions in Section <ref>).

§.§ Training and Test Data-Sets

The training and test data-sets are generated via different system inputs. System inputs comprise combinations of sinusoidals at different frequency ranges (to cover different range of inputs). This type of multi-sine signal is usual for system identification of the considered case study <cit.>. The system input $ u = (u_1, u_2)$ can be written as

\begin{equation*} \label{eq:input_sim}
\begin{aligned}
u_1 = \frac{\max {\tilde \alpha_i}}{\sum_{i=1}^{\tilde n} \tilde \alpha_i} \sum_{i=1}^{\tilde n} \tilde \alpha_i \sin \left(\tilde \omega_i t+\tilde \phi_i\right), \quad t \in[0, t_f]  \quad
u_2 = 0,
\end{aligned}
\end{equation*}

where the parameters

$$\left(\tilde \alpha_i\right)_{i=1}^{\tilde n}, \left(\tilde \omega_i\right)_{i=1}^{\tilde n}, \left(\tilde \phi_i\right)_{i=1}^{\tilde n},$$

and ${\tilde n}$ are randomly drawn from uniform distributions in the intervals $[0.01,0.1]$ $m, [0.6 \pi,3 \pi]$ $rad/s , [0,0.94\pi]$ $rad,$ and $[2,10]$, respectively. The parameter $t_f$ indicates the final time of the simulation.

Five training data sets and 1000 test data sets of inputs have been drawn from the aforementioned system input distributions.

Histogram for comparison of system outputs and the outputs of different models with cubic basis function (i.e., $h = (V_\eta x)^{\circ 3}$) for the test data.

§.§ Model Accuracy

Next, we compare the accuracy of the trained uncertainty models via the proposed methods. We compare the results against the true system model used to generate the data. We fit three extended models using the known physics-based model and the learned uncertainty models via the proposed methods described in Sections <ref> and <ref>. Additionally, we examine a scenario where no stability guarantees for the extended model is considered (i.e., only by minimizing the cost function in (<ref>)). Thus, for the test data-set, we compare the system output with the outputs of three extended models and a model that is not extended and neglects uncertainties (i.e., only consist physics-based model).

Given the large number of test data-sets (1000 test data-sets), this comparison is depicted via the histogram shown in Figure <ref>. The horizontal axis is a metric for the output error computed via $\sum_{t=0}^{t_f = 20} \| y_s(t) - C x(t) \|$. The histogram presents the distribution of output errors for the test data by dividing it into intervals and displaying the frequency of occurrences within each interval. In the results shown in Figure <ref>, the basis function representing the true system is a cubic function (i.e., $h = (V_\eta x)^{\circ 3} $ in (<ref>)). For results considering basis functions other than the cubic function, see Section <ref>. From Figure <ref>, one can conclude that both learning strategies (the results in blue and magenta) substantially enhance model quality compared to the model without a learned uncertainty model (green). Note that the SCP method is initialized with the solution of the approximated convex program. Other initializations, such as setting the parameters of the uncertainty model $\theta = (\Theta_l, B_l, \Theta_n)$ to zero or fixing the invariance set shape matrix $P$ as a constant multiplied by the identity matrix, did not yield a feasible solution. This further highlights the significance of the proposed convex approximation.

(The Price of Stability)
Imposing the stability constraint turns out to deteriorate the expected output errors (compare the average of the red histogram with the blue and magenta histograms in Figure <ref>). The root cause of this can be due to the fact that the stability constraints make the resulting optimization program non-convex. Namely, we have to resort to a sub-optimal solution among all possible models considered (the magenta and blue histograms in Figure <ref>), whereas ignoring this condition allows us to stay in a convex setting where we can find the global optimal model (the red histogram). It is, however, worth noting that this better performance at the average level may come at the cost of some rare instances in the tail of the histogram (see the zoom window in Figure <ref>). In this light, one can see this average performance degradation as “the price of stability" to ensure that any model we estimate meets our prior stability assumption.

§.§ Effect of Basis Functions

Figures <ref> and <ref>, similar to Fig <ref>, illustrate the histogram for the models output errors. Unlike Fig <ref>, the basis functions (i.e., $h(\cdot)$ in (<ref>)) are not same as the true nonlinearity in the system. Now, the basis functions are $h = [(V_\eta x)^{\circ 2^\top} , (V_\eta x)^{\circ 3^\top}]^\top$ and $h = [(V_\eta x)^{\circ 2^\top}, exp(V_\eta x)^\top, (V_\eta x)^{\circ 3^\top}]^\top$ in Figures <ref> and <ref>, respectively. The results in these two figures indicate that the selection of the basis function does not significantly affect the accuracy of the extended model, especially for the results with stability guarantees. That is, having a model that uses the “wrong" basis function is still better that using the model that does not account for uncertainty (at least in the considered use case).

Histogram for comparison of system outputs and the outputs of different models with basis function $h = [(V_\eta x)^{\circ 2^\top} , (V_\eta x)^{\circ 3^\top}]^\top$ for the test data.

Histogram for comparison of system outputs and the outputs of different models with correct basis function $h = [(V_\eta x)^{\circ 2^\top}, exp(V_\eta x)^\top, (V_\eta x)^{\circ 3^\top}]^\top$ for the test data.

§ CONCLUSION

This paper has proposed a framework for model updating via learning modeling uncertainties in locally (and globally) Lipschitz nonlinear models. Moreover, we have ensured the stability of the extended model, which consists of a prior known model and the learned uncertainty model. The proposed framework includes two steps. First, we have assumed that uncertainty and state estimates are known. Under this assumption, we have introduced two semi-definite programs: Theorems <ref> and <ref> for locally (and Theorems <ref> and <ref> for globally) Lipschitz nonlinear models, to learn uncertainty models while ensuring the stability of the extended model (invariant set and ISS for locally and globally Lipschitz nonlinear models, respectively). Second, we have proposed a filter, designed using the semi-definite program outlined in Proposition <ref>, to estimate uncertainty and state, based on the known prior model and input-output data. Simulations for a large data-set of a roll plane model of a vehicle have demonstrated the performance and potential of the proposed approach.

§ APPENDIX

§ INVARIANT SET CONDITIONS FOR MODEL (<REF>)

Let $W(x) = x^\top P x$ with $P \succ 0$ be a Lyapunov function candidate. We need to derive the derivative of the Lyapunov function candidate along the trajectories of the extended system model (<ref>). Before that, first, we impose the condition $\mathcal{E}_{inv} \subseteq \mathcal{E}_{sys}$ on the model invariant set and the system state set as defined in (<ref>) and (<ref>), respectively. By doing so, later, the Lipschitz property of the nonlinearities of extended system model (<ref>) can be used in deriving a tractable condition for $\dot{W}~\leq~0$. It is shown in <cit.> that the condition $\mathcal{E}_{inv} \subseteq \mathcal{E}_{sys}$ can be written as (<ref>) by applying $\mathcal{S}$-procedure.

From (<ref>) and the Lipschitz conditions for the known nonlinearity $g(\cdot)$ and the uncertainty model nonlinearity $h(\cdot)$ in (<ref>) (together with that zero is an equilibrium point of system for $u=0$), it follows that

\begin{equation}
\begin{aligned} \dot{W}(x) = &x^\top \left(A^\top P+P A+  V_\eta ^\top \Theta_l^\top S_{\eta_l}^\top  P+ P S_{\eta_l} \Theta_l V_\eta \right) x+2 x^\top P ( S_g g( V_g x,u)+S_{\eta_l} \Theta_n h( V_\eta x,u)) \\
& +2 x^\top P (B_u +S_{\eta_l} B_l)u\\
\leq &x^\top \left(A^\top P+P A+  V_\eta ^\top \Theta_l^\top S_{\eta_l}^\top  P+ P S_{\eta_l} \Theta_l V_\eta \right) x+2 \|x^\top P S_g g( V_g x,u) \|\\
&+ 2 \|x^\top P S_{\eta_l} \Theta_n h( V_\eta x,u)\| +2 x^\top P (B_u +S_{\eta_l} B_l)u\\
\leq &x^\top \left(A^\top P+P A+  V_\eta ^\top \Theta_l^\top S_{\eta_l}^\top  P+ P S_{\eta_l} \Theta_l V_\eta \right) x +2 \|x^\top P S_g \| \| g( V_g x,u) \|\\
&+ 2 \|x^\top P S_{\eta_l} \| \| \Theta_n h (V_\eta x,u)\| +2 x^\top P (B_u +S_{\eta_l} B_l)u\\
\leq &x^\top \left(A^\top P+P A+  V_\eta ^\top \Theta_l^\top S_{\eta_l}^\top  P+ P S_{\eta_l} \Theta_l V_\eta \right) x +2 \|x^\top P S_g \|  (l_{g_x} \|V_g x\| + l_{g_u} \|u\|) \\
&+ 2 \|x^\top P S_{\eta_l} \| \| \Theta_n \| (l_{h_x} \|V_\eta x\| + l_{h_u} \|u\|) +2 x^\top P (B_u +S_{\eta_l} B_l)u.
\label{eq:lyapanov_ap_1}
\end{aligned}
\end{equation}

Now, by imposing the condition on parameters of the nonlinear part of uncertainty model given in (<ref>) and then a specific form of the Cauchy–Schwarz inequality an upper bound for the right-hand-side of (<ref>) can be found as follows:

\begin{equation}
\begin{aligned}
\dot{W}(x)\leq &x^\top \left(A^\top P+P A+  V_\eta ^\top \Theta_l^\top S_{\eta_l}^\top  P+ P S_{\eta_l} \Theta_l V_\eta \right) x+2 \|x^\top P S_g \|  (l_{g_x} \|V_g x\| + l_{g_u} \|u\|) \\
&+ 2 \|x^\top P S_{\eta_l} \| (\bar l_{h_x} \|V_\eta x\| + \bar l_{h_u} \|u\|) +2 x^\top P (B_u +S_{\eta_l} B_l)u\\
\leq &x^\top \left(A^\top P+P A+  V_\eta ^\top \Theta_l^\top S_{\eta_l}^\top  P+ P S_{\eta_l} \Theta_l V_\eta \right) x+ (l_{g_x} +l_{g_u}) \|x^\top P S_g \|^{2}+ l_{g_x}\|V_g x\|^{2} + l_{g_u} \|u\|^{2} \\
& + (\bar l_{h_x} +\bar l_{h_u}) \|x^\top P S_{\eta_l} \|^{2}+ \bar l_{h_x} \|V_\eta x\|^{2} + \bar l_{h_u} \|u\|^{2} +2 x^\top P (B_u +S_{\eta_l} B_l)u\\
= &x^\top\big(A^\top P+P A+  V_\eta ^\top \Theta_l^\top S_{\eta_l}^\top  P+ P S_{\eta_l} \Theta_l V_\eta + (l_{g_x} +l_{g_u}) P S_g S_g^\top P + (\bar l_{h_x} + \bar l_{h_u}) P S_{\eta_l} S_{\eta_l}^\top P \\
&+l_{g_x} V_g^\top V_g + \bar l_{h_x} V_\eta^\top V_\eta \big) x + (l_{g_u} + \bar l_{h_u}) u^\top u  +2 x^\top P (B_u +S_{\eta_l} B_l)u  \\
= & x^\top \Delta x + 2 x^\top P (B_u+S_{\eta_l} B_l) u +(l_{g_u} + \bar l_{h_u}) u^\top u
\end{aligned}
\label{eq:lyapanov_ap}
\end{equation}

with $\bar l_{h_u}$ as defined in (<ref>) and $\Delta$ as defined in (<ref>). Therefore, from (<ref>), it follows that

\begin{equation}
\begin{aligned} \dot{W}(x) \leq& x^\top \Delta x + 2 x^\top P (B_u+S_{\eta_l} B_l) u +(l_{g_u} + \bar l_{h_u}) u^\top u.
\end{aligned}
\label{eq:lyapanov}
\end{equation}

If we can find a $P$ such that $\dot{W} \leq 0$ whenever $u \in \mathcal{E}_{u}$ (as defined in (<ref>)) and $W \geq 1$ along the trajectories of (<ref>); Then, the ellipsoid $\mathcal{E}_{inv}$ (as defined in (<ref>)) is an invariant set of (<ref>) [32]. To this end, consider the vector $\left[x^\top u^\top 1 \right]^\top$. Using (<ref>), the condition $\dot{W} \leq 0$ can be stated as

\begin{equation} \label{eq:vdot_cond}
E_1 = \left[\begin{array}{ccc}
\Delta  & P (B_u+S_{\eta_l} B_l)  & 0 \\
\star & (l_{g_u} + \bar l_{h_u}) I & 0 \\
\star & \star & 0
\end{array}\right] \preceq 0.
\end{equation}

Similarly, the condition (<ref>) can be written as

\begin{equation} \label{eq:u_cond}
F_1 =\left[\begin{array}{ccc}
0  & 0 & 0 \\
\star & U & 0 \\
\star & \star & -1
\end{array}\right] \preceq 0.
\end{equation}

Finally, the condition $W(x) \geq 1$ can be restated as

\begin{equation} \label{eq:v_cond}
G_1 =\left[\begin{array}{ccc}
-P  & 0 & 0 \\
\star & 0 & 0 \\
\star & \star & 1
\end{array}\right] \preceq 0.
\end{equation}

Therefore, for constructing an invariant set of trajectories, we require (<ref>) ($\dot{W} \leq 0$) to hold when (<ref>) and (<ref>) hold. To satisfy these conditions, we can apply $\mathcal{S}$-procedure <cit.>, which states that there should exist non-negative constants $\alpha$ and $\beta$ such that the following inequality holds:

\begin{equation}
E_1 - \alpha F_1 - \beta G_1 \preceq 0.
\end{equation}

The condition above can be written as (<ref>). Note that in the derivation of (<ref>), we have imposed the condition on Lipschitz set condition in (<ref>) and the norm constraint on parameters of nonlinear part of uncertainty model in (<ref>). Therefore, the condition (<ref>) together with (<ref>), and (<ref>) imply that the ellipsoid $\mathcal{E}_{inv}$ (as defined in (<ref>)) is an invariant set of the extended model (<ref>).

§ ISS CONDITIONS FOR MODEL (<REF>)

Let $W(x):={x}^\top P {x}$ be an ISS Lyapunov function candidate. Similar to the proof in Appendix <ref>, given the globally Lipschitz property of the model nonlinearities in (<ref>) and the imposed condition on parameters of nonlinear part of uncertainty model in (<ref>), the time-derivative of the candidate ISS Lyapunov function, $\dot{W}$, can be written in the form of (<ref>). The rest of the proof follows the line of reasoning of the proof of Proposition 1 in [22], which for brevity is omitted here. It is shown in [22] for a Lyapunov derivative inequality similar to (<ref>) that, in order to conclude ISS property of a dynamics, $\Delta$ has to be negative definite. Note that in the derivation of $\Delta$ in (<ref>), we have imposed the condition on the parameters of nonlinear part of uncertainty model in (<ref>). Therefore, the condition $\Delta \prec 0$ together with (<ref>) imply that the extended model (<ref>) is ISS with respect to input $u$.

§ PROOF OF THEOREM <REF>

We show that the non-convex optimization problem in (<ref>) can be convexified as in (<ref>). First, to make the problem tractable, we set the $S_{\eta_l} = I$ in extended system model (<ref>). Then, we convexify the non-convex conditions in (<ref>) as follows.

Consider the stability constraint (<ref>). By applying changes of variables as $S := P \Theta_l$, $R := P B_l$ and the Schur complement, the stability constraint in (<ref>) is equivalent to (<ref>). Favorably, by the introduced change of variables, the stability constraint becomes convex (an LMI). However, since the cost $J$ in (<ref>) is a function of $\Theta_l = P^{-1} S$ and $B_l = P^{-1} R$, the cost is not convex in $S, R$ and $P$, after the changes of variables. Therefore, we have to convexify the cost to arrive at a convex optimization problem formulation. For scalar cost function $J$ in (<ref>), we have the following equality:

J = \sum_{i = 1}^{N} d_i^\top T^\top T d_i = \sum_{i = 1}^{N} tr(d_i^\top T^\top T d_i).

Due to cyclic property of the trace operator, the above cost can be written as

J = tr(T^\top T \sum_{i = 1}^{N} d_i d_i^\top).

By the definition of the matrix $D$ in (<ref>) and the cyclic property of the trace, we have

J = tr(T D T^\top).

Now, we can write the epigraph form of the optimization problem in (<ref>) as follows:
P,S, R, Θ_n, W, α, γ

(<ref>),(<ref>), (<ref>)
P≻0, α, γ≥0
tr(T D T^⊤)≤tr(W).

Due to monotonicity of trace, the last constraint in (<ref>) can be transformed to the following constraint:

\begin{equation*}
T D T^\top \preceq W.
\end{equation*}

Note that even with the application of the above transformation, the optimization remains equivalent to (<ref>) <cit.>. By applying the Schur complement to the above inequality, we have

\begin{equation*}
\left[\begin{array}{cc}
W & T\tilde{D}^\top \\
\star & I
\end{array}\right] \succeq 0.
\end{equation*}

Note that by construction, the data matrix $D$ is always symmetric and positive semi-definite. Therefore, its Cholesky decomposition (i.e., $D = \tilde{D}^\top \tilde{D}$) always exists. By applying the congruence transformation of $\text{diag}(P,I)$ to the above inequality, we obtain the following equivalent inequality:

\begin{equation} \label{eq:GDG}
\left[\begin{array}{cc}
P W P &  \tilde{T} \tilde{D}^\top \\
\star & I
\end{array}\right] \succeq 0.
\end{equation}

Note that in the above inequality to introduce $\tilde{T} := P T = {\left[\begin{array}{cccc}
S  & R & -P & P \Theta_n
\end{array}\right]}$, the change of variable as $Z := P \Theta_n$ is applied.  Now, by substituting the lower bound $2\mu_2 P-\mu_2^2 W^{-1}$, with given positive scalar $\mu_2$, for $P W P$ and applying the Schur complement, the LMI in (<ref>) is obtained.

Consider another non-convex condition (upper bound on the parameters of nonlinear part of uncertainty model) as (<ref>). Based on <cit.>, the condition (<ref>) can be written as

\begin{equation} \label{eq:theta_n_bound_matrix_form}
\left[\begin{array}{cc}
\bar l_{h_x} I &  l_{h_x} \Theta_n^\top \\
\star & \bar l_{h_x} I
\end{array}\right] \succeq 0.
\end{equation}

By applying the congruence transformation of $\text{diag}(I,P)$ to the above inequality, we obtain the following equivalent inequality

\begin{equation*}
\left[\begin{array}{cc}
\bar l_{h_x} I &  l_{h_x} Z^\top \\
\star & \bar l_{h_x} P^2
\end{array}\right] \succeq 0.
\end{equation*}

Now, by substituting the lower bound $2\mu_1 P-\mu_1^2 I$, with given positive scalar $\mu_1$, for $P^2$, the LMI in (<ref>) is obtained.

In conclusion, instead of the non-convex optimization problem in (<ref>), we provide an approximation of that problem in the form of the semi-definite program in (<ref>). We remark that the approximation arises from initially setting $S_{\eta_l} = I$ at the beginning of the proof. Additionally, we use the lower bound of $2\mu_2 P-\mu_2^2 W^{-1}$ for $P W P$ in the derivation of the LMI in (<ref>), and the lower bound of $2\mu_1 P-\mu_1^2 I$ for $P^2$ in the derivation of the LMI in (<ref>).

§ PROOF OF THEOREM <REF>

Here, we follow the same line of reasoning as in the proof of Theorem <ref> by showing that the non-convex optimization problem in (<ref>) can be convexified as in (<ref>). First, we set $S_{\eta_l} = S_\eta$ in extended model (<ref>) (i.e., here we use the knowledge of uncertainty structure). Then, we convexify the non-convex conditions as follows.

Consider the constraint (<ref>). By applying the congruence transformation of $\text{diag}(P^{-1},I,I)$ to (<ref>), we obtain the following equivalent inequality:

\begin{equation} \label{eq:Q_form_delta_condition_locally}
\begin{aligned}
\left[\begin{array}{ccc}
Q \Delta Q +\beta Q  & B_u+S_{\eta} B_l & 0 \\
\star & (l_{g_u} + \bar l_{h_u}) I - \alpha U & 0 \\
\star & \star & \alpha - \beta
\end{array}\right] \prec 0,
\end{aligned}
\end{equation}

where $Q := P^{-1}$. Using (<ref>), the term $Q \Delta Q+\beta Q $ in the above inequality can be written as follows:

\begin{equation} \label{eq:Q_form_delta_with_beta}
\begin{aligned}
Q \Delta Q+\beta Q & = A Q +Q A^\top+ Q V_\eta ^\top \Theta_l^\top S_{\eta}^\top + S_{\eta} \Theta_l V_\eta Q + (l_{g_x} +l_{g_u}) S_g S_g^\top   \\
&+ (\bar l_{h_x}+ \bar l_{h_u}) S_{\eta_l} S_{\eta_l}^\top +l_{g_x} Q V_g^\top V_g Q + \bar l_{h_x} Q  V_\eta^\top V_\eta Q.
\end{aligned}
\end{equation}

Using Young's inequality $\big(X^\top Y+Y^\top X \preceq \frac{1}{2}(X+\bar Z Y)^\top \bar Z^{-1} (X+\bar Z Y)$, with a symmetric positive definite matrix $\bar Z\big)$, we can find the following upper bound for (<ref>):

\begin{equation*}
\begin{aligned}
Q \Delta Q+\beta Q &\leq A Q +Q A^\top +  \frac{1}{2}(V_\eta ^\top \Theta_l^\top S_\eta^\top +\bar Z Q)^\top \bar Z^{-1} (V_\eta ^\top \Theta_l^\top S_\eta^\top +\bar Z Q)  \\
& + (l_{g_x} +l_{g_u}) S_g S_g^\top + (\bar l_{h_x} + \bar l_{h_u}) S_{\eta_l} S_{\eta_l}^\top +l_{g_x} Q V_g^\top V_g Q+ \bar l_{h_x} Q  V_\eta^\top V_\eta Q.
\end{aligned}
\end{equation*}

Now, using the above upper bound for $Q \Delta Q+\beta Q$ and the Schur complement, the following sufficient condition can be found for the inequality in (<ref>):

\begin{equation*}
\left[\begin{array}{cccccc}
N_{11}  & N_{12} & 0 &S_\eta \Theta_l V_\eta +Q \bar Z  &  N_{15} & N_{16}\\
\star   & N_{22} & 0 & 0 & 0 & 0 \\
\star   &  \star & N_{33} & 0 & 0 & 0\\
\star  & \star  &  \star  &  -2 \bar Z & 0 & 0 \\
\star  & \star  &  \star  & \star &  -I & 0 	\\
\star  & \star  &  \star  &  \star &  \star & -I
\end{array}\right] \prec 0
\end{equation*}

with $N_{11}, N_{12}, N_{15}, N_{16}, N_{22}$ and $N_{33}$ as defined in (<ref>). Given nonlinear term $Q\bar Z$ in the above inequality we select $\bar Z = \bar \gamma I$, with a positive scalar $\bar \gamma$. Therefore, the above inequality can be written as (<ref>). The non-convex condition (upper bound on the parameters of nonlinear part of uncertainty model) given in (<ref>) can be written as (<ref>). For the details see derivation of (<ref>) in proof of Theorem <ref>.

The Lipschitz set condition (<ref>), in terms of $Q = P^{-1}$ can be written as

\begin{equation*}
\left[\begin{array}{cc}
\gamma   F -  Q^{-1} & 0 \\
* & - \gamma +1
\end{array}\right] \preceq 0.
\end{equation*}

Now, by substituting the upper bound $\gamma F -  2 \mu_3 I + \mu_3^2 Q$, with given positive scalar $\mu_3$, for $\gamma F -  Q^{-1}$, the LMI in (<ref>) is obtained. Furthermore, the quadratic cost can be treated as constraint (<ref>) by writing the epigraph form of the optimization problem (see (<ref>) for the epigraph form in proof of Theorem <ref>). Thus, instead of the non-convex optimization problem in (<ref>), we provide the approximate semi-definite program in (<ref>).  We remark that the approximation arises from using a sufficient condition (Young's inequality together with selection of $\bar Z = \bar \gamma I$) in the derivation of the LMI in (<ref>). Additionally, we use the upper bound $\gamma F -  2 \mu_3 I + \mu_3^2 Q$ for $\gamma F -  Q^{-1}$ in the derivation of the LMI in (<ref>).

§ ROLL PLANE SYSTEM DESCRIPTION

The equations of motion of the roll plane system are given in <cit.>. Here, we provide the state space representation. By defining the state vector $x_s = [x_{s_1}, x_{s_2},\ldots,x_{s_7}, x_{s_8}]^\top := [q_1, q_2, q_3, q_4,\dot{q}_1,\ldots, \dot{q}_4]^\top$, where $q_i$ and $\dot{q}_i$ are the displacements (as shown in Figure <ref>) and the corresponding velocities, the system dynamics can be described in the form of (<ref>), where we have

\begin{equation*}
\begin{aligned}
A &=\left[\begin{array}{cc}
0    & I \\
-\tilde{M}^{-1}\tilde{K} & -\tilde{M}^{-1}\tilde{C} \end{array}\right], \quad
B_u =\left[\begin{array}{c}
0 \\
\tilde{M}^{-1} \tilde{K}_u
\end{array}\right], \quad
S_g = S_\eta =\left[\begin{array}{c}
0 \\
- \tilde{M}^{-1} \tilde{S}
\end{array}\right], \\
g(\cdot) &= 0.2 \tilde{c}_{1,n} tanh(V_g x_s) \quad
D_\nu = I, \quad B_\omega= 0, \quad
\eta(\cdot) = \delta \tilde{k}_1 V_{\eta} x_s + \tilde{k}_{1_n} (V_{\eta} x_s)^{\circ 3}, \\
\tilde{C} & 0 \\
0 & \tilde{C}
\end{array}\right], \quad
\tilde{C}=\left[\begin{array}{cccc}
1 & 0 & -1 & 0 \\
0 & 1 & 0 & -1 \\
\end{array}\right] \quad
V_g =\left[\begin{array}{cc}
0 & 10 \tilde{C}
\end{array}\right], \quad
V_\eta =\left[\begin{array}{cc}
\tilde{C}  & 0
\end{array}\right],
\end{aligned}
\end{equation*}

\begin{equation*}
\begin{aligned}
\tilde{M} &= \left[\begin{array}{cccc}
\frac{\tilde{m}}{2} & \frac{\tilde{m}}{2} & 0 & 0 \\
\frac{-\tilde{I}}{\tilde{L}} & \frac{\tilde{I}}{\tilde{L}} & 0 & 0 \\
0 & 0 & \tilde{m}_{t_1} & 0 \\
0 & 0 & 0 & \tilde{m}_{t_2}
\end{array}\right], \quad
\tilde{K} = \left[\begin{array}{cccc}
\tilde{k}_1 & \tilde{k}_2 & -\tilde{k}_1 & -\tilde{k}_2 \\
\frac{-\tilde{L}}{2}\tilde{k}_1 & \frac{\tilde{L}}{2}\tilde{k}_2 & \frac{\tilde{L}}{2}\tilde{k}_1 & \frac{-\tilde{L}}{2}\tilde{k}_2 \\
-\tilde{k}_1 & 0 & \tilde{k}_1+\tilde{k}_{t_1} & 0 \\
0 & -\tilde{k}_2 & 0 & \tilde{k}_2+\tilde{k}_{t_2}
\end{array}\right], \\
\tilde{C} &= \left[\begin{array}{cccc}
\tilde{c}_1 & \tilde{c}_2 & -\tilde{c}_1 & -\tilde{c}_2 \\
\frac{-\tilde{L}}{2}\tilde{c}_1 & \frac{\tilde{L}}{2}\tilde{c}_2 & \frac{\tilde{L}}{2}\tilde{c}_1 & \frac{-\tilde{L}}{2}\tilde{c}_2 \\
-\tilde{c}_1 & 0 & \tilde{c}_1 & 0 \\
0 & -\tilde{c}_2 & 0 & \tilde{c}_2
\end{array}\right], \quad
\tilde{K}_u =\left[\begin{array}{cc}
0 & 0\\
0 & 0 \\
\tilde{k}_{t_1} & 0 \\
0 & \tilde{k}_{t_2}  \\
\end{array}\right], \quad
\tilde{S}=\left[\begin{array}{cc}
1 & 1\\
\frac{-\tilde{L}}{2} & \frac{\tilde{L}}{2} \\
-1 & 0 \\
0 & -1  \\
\end{array}\right].
\end{aligned}
\end{equation*}

The parameter values are listed in Table <ref>.

Parameter values of roll plane system.

Parameter                                                               Value      Unit
$\tilde{m}$                   $580$                                                                            $\mathrm{~kg}$

$\tilde{m}_{t_1}, \tilde{m}_{t_2}$                              $36.26$                                                                  $\mathrm{~kg}$

$\tilde{I}$                              $63.3316$                                                                  $\mathrm{~kg} \cdot \mathrm{m}^2$

$\tilde{L}$                              $1.524$                                                                  $\mathrm{~m}$

$\tilde{c}_1, \tilde{c}_2$                              $710.70$                                                                  $\frac{\mathrm{N} \cdot \mathrm{s}}{\mathrm{m}}$

$\tilde{c}_{1,n}, \tilde{c}_{2,n}$                      $0.71$                                                                  $\frac{\mathrm{N} \cdot \mathrm{s}}{\mathrm{m}}$

$\tilde{k}_1, \tilde{k}_2$                              $19357.2$                                                                  $\frac{\mathrm{N}}{\mathrm{m}}$

$\tilde{k}_{t_1}, \tilde{k}_{t_2}$                              $96319.76$                                                                  $\frac{\mathrm{N}}{\mathrm{m}}$

$\tilde{k}_{1_n}, \tilde{k}_{2_n}$                              $15000$                                                                  $\frac{\mathrm{N}}{\mathrm{m}^3}$

$\delta \tilde{k}_1, \delta \tilde{k}_2$                              $5807.2$                                                                  $\frac{\mathrm{N}}{\mathrm{m}}$

Uncertainty entries and their estimates  for a randomly selected test data.

Deriving terms of the uncertainty $\eta(\cdot)$, i.e., $V_\eta x_s$ and their estimates, i.e., $V_\eta \hat{x}_s$, for a randomly selected test data.

§ RESULTS FOR UNCERTAINTY AND STATE ESTIMATION

For one test data set in the considered case study (i.e., the roll plane system), the entries of uncertainty $\eta$ in (<ref>) and their estimates are depicted in Figure <ref>. Moreover, Figure <ref> illustrates the deriving terms of the uncertainty $\eta(\cdot)$, i.e., $V_\eta x_s$ and their estimations, i.e., $V_\eta \hat{x}_s$. As demonstrated, the proposed uncertainty and state estimation method yield highly accurate results.

§ ACKNOWLEDGEMENTS

This publication is part of the project Digital Twin project
4.3 with project number P18-03 of the research programme
Perspectief which is (mainly) financed by the Dutch Research
Council (NWO).

[1]
Ian R Manchester, Max Revay, and Ruigang Wang.
Contraction-based methods for stable identification and robust
machine learning: a tutorial.
In 2021 60th IEEE Conference on Decision and Control (CDC),
pages 2955–2962. IEEE, 2021.

[2]
Amirhassan Abbasi and C Nataraj.
Physics-informed machine learning for uncertainty reduction in time
response reconstruction of a dynamic system.
IEEE Internet Computing, 26(4):35–44, 2022.

[3]
Alireza Yazdani, Lu Lu, Maziar Raissi, and George Em Karniadakis.
Systems biology informed deep learning for inferring parameters and
hidden dynamics.
PLoS computational biology, 16(11):e1007575, 2020.

[4]
Mitchell Daneker, Zhen Zhang, George Em Karniadakis, and Lu Lu.
Systems biology: Identifiability analysis and parameter
identification via systems-biology-informed neural networks.
In Computational Modeling of Signaling Networks, pages 87–105.
Springer, 2023.

[5]
Karan Taneja, Xiaolong He, QiZhi He, Xinlun Zhao, Yun-An Lin, Kenneth J Loh,
and Jiun-Shyan Chen.
A feature-encoded physics-informed parameter identification neural
network for musculoskeletal systems.
Journal of biomechanical engineering, 144(12):121006, 2022.

[6]
Steven L Brunton, Joshua L Proctor, and J Nathan Kutz.
Discovering governing equations from data by sparse identification of
nonlinear dynamical systems.
Proceedings of the national academy of sciences,
113(15):3932–3937, 2016.

[7]
Kathleen Champion, Bethany Lusch, J Nathan Kutz, and Steven L Brunton.
Data-driven discovery of coordinates and governing equations.
Proceedings of the National Academy of Sciences,
116(45):22445–22451, 2019.

[8]
Kathleen Champion, Peng Zheng, Aleksandr Y Aravkin, Steven L Brunton, and
J Nathan Kutz.
A unified sparse optimization framework to learn parsimonious
physics-informed models from data.
IEEE Access, 8:169259–169271, 2020.

[9]
Jean-Christophe Loiseau and Steven L Brunton.
Constrained sparse galerkin regression.
Journal of Fluid Mechanics, 838:42–67, 2018.

[10]
Mariane Yvonne Schneider, Ward Quaghebeur, Sina Borzooei, Andreas Froemelt,
Feiyi Li, Ramesh Saagi, Matthew J Wade, Jun-Jie Zhu, and Elena Torfs.
Hybrid modelling of water resource recovery facilities: status and
Water Science and Technology, 85(9):2503–2524, 2022.

[11]
William Bradley, Jinhyeun Kim, Zachary Kilwein, Logan Blakely, Michael
Eydenberg, Jordan Jalvin, Carl Laird, and Fani Boukouvala.
Perspectives on the integration between first-principles and
data-driven modeling.
Computers & Chemical Engineering, 166:107898, 2022.

[12]
Ward Quaghebeur, Ingmar Nopens, and Bernard De Baets.
Incorporating unmodeled dynamics into first-principles models through
machine learning.
IEEE Access, 9:22014–22022, 2021.

[13]
Seth L Lacy and Dennis S Bernstein.
Subspace identification with guaranteed stability using constrained
IEEE Transactions on automatic control, 48(7):1259–1263, 2003.

[14]
Loris Di Natale, Muhammad Zakwan, Philipp Heer, Giancarlo Ferrari Trecate, and
Colin N Jones.
Simba: System identification methods leveraging backpropagation.
arXiv preprint arXiv:2311.13889, 2023.

[15]
Farhad Ghanipoor, Carlos Murguia, Peyman Mohajerin Esfahani, and Nathan van de
Uncertainty learning for lti systems with stability guarantees.
available at arXiv:2310.20568, 2023.

[16]
Mohammad Khosravi and Roy S Smith.
Nonlinear system identification with prior knowledge on the region of
IEEE Control Systems Letters, 5(3):1091–1096, 2020.

[17]
Mohammad Khosravi.
Representer theorem for learning koopman operators.
IEEE Transactions on Automatic Control, 2023.

[18]
MF Shakib, Roland Tóth, AY Pogromsky, A Pavlov, and Nathan van de Wouw.
Kernel-based learning of stable nonlinear state-space models.
In 2023 62nd IEEE Conference on Decision and Control (CDC),
pages 2897–2902. IEEE, 2023.

[19]
Max Revay, Ruigang Wang, and Ian R Manchester.
Recurrent equilibrium networks: Unconstrained learning of stable and
robust dynamical models.
In 2021 60th IEEE Conference on Decision and Control (CDC),
pages 2282–2287. IEEE, 2021.

[20]
MF Shakib, AY Pogromsky, Alexey Pavlov, and Nathan van de Wouw.
Fast identification of continuous-time lur’e-type systems with
stability certification.
IFAC-PapersOnLine, 52(16):227–232, 2019.

[21]
Max Revay, Ruigang Wang, and Ian R Manchester.
Recurrent equilibrium networks: Flexible dynamic models with
guaranteed stability and robustness.
IEEE Transactions on Automatic Control, 2023.

[22]
Farhad Ghanipoor, Carlos Murguia, Peyman Mohajerin Esfahani, and Nathan van de
Robust fault estimators for nonlinear systems: An ultra-local model
available at arXiv:2305.14036, 2023.

[23]
Stephen Boyd, Laurent El Ghaoui, Eric Feron, and Venkataramanan Balakrishnan.
Linear matrix inequalities in system and control theory.
SIAM, 1994.

[24]
Eduardo D Sontag et al.
On the input-to-state stability property.
Eur. J. Control, 1(1):24–36, 1995.

[25]
Eduardo D Sontag.
Input to state stability: Basic concepts and results.
In Nonlinear and optimal control theory, pages 163–220.
Springer, 2008.

[26]
Jian Han, Xiuhua Liu, Xianwen Gao, and Xinjiang Wei.
Intermediate observer-based robust distributed fault estimation for
nonlinear multiagent systems with directed graphs.
IEEE Transactions on Industrial Informatics, 16(12):7426–7436,

[27]
B. Jiang, M. Staroswiecki, and V. Cocquempot.
Fault accommodation for nonlinear dynamic systems.
IEEE Transactions on Automatic Control, 51:1578–1583, 2006.

[28]
Farhad Ghanipoor, Carlos Murguia, Peyman Mohajerin Esfahani, and Nathan van de
Ultra local nonlinear unknown input observers for robust fault
In 2022 IEEE 61st Conference on Decision and Control (CDC),
pages 918–923. IEEE, 2022.

[29]
Farhad Ghanipoor, Carlos Murguia, Peyman Mohajerin Esfahani, and Nathan van de
Linear fault estimators for nonlinear systems: An ultra-local model
IFAC-PapersOnLine, 56(2):11693–11698, 2023.

[30]
Bas M Kessels, Rob HB Fey, and Nathan van de Wouw.
Real-time parameter updating for nonlinear digital twins using
inverse mapping models and transient-based features.
Nonlinear Dynamics, 111(11):10255–10285, 2023.

[31]
Emmanuel Blanchard.
Polynomial chaos approaches to parameter estimation and control
design for mechanical systems with uncertain parameters.
PhD thesis, Virginia Tech, 2010.

[32]
Yankai Lin, Michelle S Chong, and Carlos Murguia.
Plug-and-play secondary control for safety of lti systems under
available at arXiv:2212.00593, 2022.

[33]
Ryan James Caverly and James Richard Forbes.
Lmi properties and applications in systems, stability, and control
available at arXiv:1903.08599, 2019.