Physics-guided neural networks for feedforward control with
  input-to-state stability guarantees
The increasing demand on precision and throughput within high-precision
mechatronics industries requires a new generation of feedforward controllers
with higher accuracy than existing, physics-based feedforward controllers. As
neural networks are universal approximators, they can in principle yield
feedforward controllers with a higher accuracy, but suffer from bad
extrapolation outside the training data set, which makes them unsafe for
implementation in industry. Motivated by this, we develop a novel
physics-guided neural network (PGNN) architecture that structurally merges a
physics-based layer and a black-box neural layer in a single model. The
parameters of the two layers are simultaneously identified, while a novel
regularization cost function is used to prevent competition among layers and to
preserve consistency of the physics-based parameters. Moreover, in order to
ensure stability of PGNN feedforward controllers, we develop sufficient
conditions for analyzing or imposing (during training) input-to-state stability
of PGNNs, based on novel, less conservative Lipschitz bounds for neural
networks. The developed PGNN feedforward control framework is validated on a
real-life, high-precision industrial linear motor used in lithography machines,
where it reaches a factor 2 improvement with respect to physics-based
mass-friction feedforward and it significantly outperforms alternative neural
network based feedforward controllers.
Max Bolderman, Hans Butler, Sjirk Koekebakker, Eelco van Horssen, Ramidin Kamidi, Theresa Spaan-Burke, Nard Strijbosch, Mircea Lazar
# Physics-guided neural networks for feedforward control

###### Abstract

Currently, there is an increasing interest in merging physics-based methods and artificial intelligence to push performance of feedforward controllers for high-precision back-reformarchics beyond what is achievable with linear feedforward control. In this paper, we develop a systematic design procedure for feedforward control using physics-guided neural networks (PGNNs) that can handle nonlinear and unknown dynamics. PGNNs effectively merge physics-based and NN-based models, and thereby result in nonlinear feedforward controllers with higher performance and the same reliability as classical, linear feedforward controllers. In particular, conditions are presented to validate (after training) and impose (before training) input-to-state stability (ISS) of PGNN feedforward controllers. The developed PGNN feedforward control framework is validated on a real-life, high-precision industrial linear motor used in lithography machines, where it reaches a factor \(2\) improvement with respect to conventional mass-friction feedforward.

 Feedforward control, neural networks, nonlinear system identification, high-precision back-reformarchics, linear motors.

## I Introduction

The field of high-precision methodronics requires continuously innovating control methods to facilitate the ever-increasing demands on both throughput as well as accuracy. For example, wafer scanners in lithography machines used for semiconductor manufacturing [1] require sub-nanometer position accuracy at velocities and accelerations exceeding \(1\ \frac{m}{s}\) and \(30\ \frac{m}{s^{2}}\), respectively, see [2, Chapter 9]. On a similar note, the ability to increase the throughput while decreasing the position error can allow for the use of components that are manufactured with larger tolerances and thereby improving the cost effectiveness when no further accuracy and throughput improvements are required. This is an objective in, for example, the manufacturing industry of printing applications using relatively low-cost stepping motors [3].

Feedforward control is a dominant actor in achieving this high position accuracy, while feedback control predominantly concerns closed-loop stability and disturbance rejection [4]. Robustness against varying references is achieved by generating the feedforward input by passing the reference through a model of the inverse system dynamics, see, e.g., [5, 6, 7]. As a result, performance of these so-called inversion-based feedforward controllers is limited by the accuracy of the model of the inverse dynamics [8].

Conventional feedforward control methods employ parametric models of the inverse system derived from physical knowledge, which are often linear [5, 7], but can also be linear-in-the-parameters (LIP) [6, 9] or nonlinear [10, 11]. Advantages of using physics-based models are the interpretable nature of the models with robust performance due to the reliance on physical knowledge, as well as a systematic stability analysis based on linear system theory. Moreover, methods exist for constructing stable inverses of linear nonminimum phase systems, i.e., systems that have an unstable inverse, using approximate solutions or non-causal filtering [12, 13]. Additionally, the training of linear and LIP models, i.e., offline optimization to find the model parameters for which the model fits the data, is typically a convex optimization problem. On the downside however, using models derived from physical knowledge yields structural model errors due to the intrinsic simplifications when deriving such models, e.g., when neglecting parasitic effects present in real-life applications such as electromagnetic distortions and complex friction forces [14]. Therefore, more general model structures that can learn nonlinear and unknown parasitic effects are needed for improving performance of feedforward controllers.

Black-box neural network (NN) models are a good candidate because of their universal approximation capabilities [15], and have been already used in feedforward control in, e.g., [16, 17]. The efficient backpropagation algorithm [18] enables training using large amounts of data, which is required for high-precision back-reformarchic systems due to their high sampling rates (\(>\ \mathcal{O}(10^{3})\ Hz\)). Currently, however, the use of NN-based feedforward controllers in practical applications remains limited, see, e.g., [19, 20] for experimental results. This is explained by the fact that using NNs comes with more difficult design challenges compared to linear physics-based models. For example, NNs require a tedious hyperparameter tuning [21]; the NN training can get stuck in local minima; and NNs suffer from poor extrapolation outside the training data [22]. Moreover, validating stability of NN-based feedforward controllers is non-trivial [23, 24], and the linear stable inversion tools do not apply for the nonlinear NN models.
The above analysis suggests that a systematic method for effectively merging physics-based models and NNs would enable feedforward controllers with higher performance for mechatronic systems with complex nonlinear dynamics due to the presence of parasitic effects. In this paper we develop a generalized feedforward control design methodology based on physics-guided neural networks (PGNNs), which effectively merge physics-based and NN-based models within a common model class. PGNNs, as described in this paper, were originally developed and employed in feedforward control design in [25], and can be regarded as a generalization of the parallel linear-NN model that has been used in nonlinear system identification, see, e.g., [26, Chapter 21]. The main contribution of this paper is the development of a generalized feedforward control design method based on PGNNs, with the following key features:

1. _Regularized training_ of PGNNs based on identified physical model parameters, to avoid competition among the NN and physical layer in the PGNN;
2. _Regularized training_ of PGNNs based on the physical model output to promote physical model compliance, and enhance robustness to non-training data;
3. _Input-to-state stability (ISS) guarantees_ for nonlinear PGNN-based feedforward controllers;
4. _Experimental validation_ on a real-life high-precision industrial linear motor for lithography machines.

**Remark I.1**: _Compared to previous conference papers [25, 27], we present the following original contributions in this journal paper: \((i)\) methods for regularized training, including the NN weights, and methods for optimal selection of the regularization parameter; \((ii)\) graceful degradation of the PGNN feedforward controller when it is operated on conditions that were not present in the training data; \((iii)\) ISS guarantees for PGNN feedforward controllers, which, in combination with stable inversion methods for linear systems, enable the design of a stable nonlinear feedforward controller for nonminimum phase systems; \((iv)\) novel, real-life experimental results._

The remainder of this paper is organized as follows: Section II introduces inversion-based feedforward control, followed by the problem statement in Section III. The methodology for PGNN feedforward controller design with regularized training and graceful degradation is presented in Section IV. Methods for testing and imposing ISS of PGNN feedforward controllers are developed in Section V. Section VI demonstrates effectiveness of the PGNN feedforward controller on a real-life coreless linear motor and a nonminimum phase simulation example. The main conclusions are summarized in Section VII.

For streamlining exposition of the results, in this paper all proofs are reported in Appendices.

## II Inversion-based feedforward control

### _Feedback-feedforward control architecture_

Fig. 1 displays a standard feedback-feedforward control scheme, where \(u(t)\) is the control input, and \(y(t)\) the system output, with time \(t\in\mathbb{R}_{\geq 0}\). A mechatronic motion system typically has a dominant linear part, denoted by \(G(s)\) with \(s\) the Laplace variable, and possible higher order modes \(\sum_{i}G_{i}(s)\). Additionally, \(g(\cdot)\) represents nonlinear parasitic effects that typically arise from manufacturing tolerances. The goal is to minimize the tracking error \(e(k):=r(k)-y(k)\), where \(r(k)\) is the reference, and \(y(k)\) is the measured at discrete-time indices \(k\in\mathbb{Z}_{\geq 0}\).

The control input \(u(t)\) is computed at discrete-time instances \(k\) according to

\[u(k)=u_{\text{fb}}(k)+u_{\text{ff}}(k), \tag{1}\]

where \(u_{\text{fb}}(k)\) is the feedback, and \(u_{\text{ff}}(k)\) the feedforward input. The zero-order-hold (ZOH) in Fig. 1 is a discrete-to-continous (D2C) operator, which lets \(u(t)=u(k)\) for \(t\in[kT_{s},(k+1)T_{s})\), with sampling time \(T_{s}\in\mathbb{R}_{>0}\). The feedback input is of the form

\[u_{\text{fb}}(k)=C(q)e(k), \tag{2}\]

where \(q\) is the forward-shift operator, e.g., \(e(k)=qe(k-1)\), and \(C(q)\) a rational transfer function. The design of the feedforward controller is presented in the following sections.

**Remark II.1**: _In this paper we focus on the tracking problem, i.e., minimizing the output tracking error \(e(k)\) with respect to a desired reference \(r(k)\). Nevertheless, the methods proposed in this work can be extended to reject known disturbances acting on the closed-loop system._

### _Linear feedforward controller design_

As an introductory example to feedforward controller design, we consider a moving mass experiencing viscous friction, i.e.,

\[m\ddot{y}(t)=u(t)-f_{v}\dot{y}(t),\quad t\in\mathbb{R}_{\geq 0}. \tag{3}\]

In (3), \(m\in\mathbb{R}_{>0}\) is the mass, \(f_{v}\in\mathbb{R}_{>0}\) the viscous friction coefficient, and \(u(t)\) and \(y(t)\) are respectively the force input and position output. Application of the Laplace transform to (3) gives

\[Y(s)=G(s)U(s):=\frac{1}{ms^{2}+f_{v}s}U(s). \tag{4}\]

The moving mass in (4) is a simple case of the general system displayed in Fig. 1, i.e., without higher order dynamics \(\sum_{i}G_{i}(s)\) and nonlinear parasitic effects \(g(\cdot)\).

In order to derive a discrete-time feedforward controller, we discretize the dynamics (4) using ZOH, which gives

\[y(k)=G(q)u(k):=\frac{b_{1}q^{-1}+b_{2}q^{-2}}{1+a_{1}q^{-1}+a_{2}q^{-2}}u(k), \tag{5}\]

where \(a_{1},a_{2},b_{1},b_{2}\in\mathbb{R}\) are the discrete-time transfer function coefficients. Suppose that the feedforward yields perfect tracking, i.e., \(e(k)=r(k)-y(k)=0\), then (2) gives \(u_{\text{fb}}(k)=0\), such that \(u(k)=u_{\text{ff}}(k)\) from (1). Correspondingly, from (5) we compute the feedforward controller
as

\[u_{\text{ff}}(k)=G^{-1}(q)r(k)=\frac{1+a_{1}q^{-1}+a_{2}q^{-2}}{b_{1 }q^{-1}+b_{2}q^{-2}}r(k) \tag{6}\] \[=\frac{1}{b_{1}}r(k+1)+\frac{a_{1}}{b_{1}}r(k)+\frac{a_{2}}{b_{2}} r(k-1)-\frac{b_{2}}{b_{1}}u_{\text{ff}}(k-1).\]

The linear feedforward controller design directly extends to systems with higher order dynamics. Then, (5) becomes

\[y(k)=G(q)u(k):=\frac{q^{-n_{k}}\sum_{i=1}^{n_{k}}b_{i}q^{-i}}{1+ \sum_{i=1}^{n_{a}}a_{i}q^{-i}}u(k), \tag{7}\]

with \(n_{a}\in\mathbb{Z}_{\geq 0}\), \(n_{b}\in\mathbb{Z}_{>0}\) the order of the dynamics, and \(n_{k}\in\mathbb{Z}_{\geq 0}\) the number of pure input delays. We denote \(a_{0}=1\) and, following the same reasoning as for (6), obtain the feedforward controller corresponding to (7) as

\[u_{\text{ff}}(k)=\sum_{i=0}^{n_{a}}\frac{a_{i}}{b_{1}}r(k+n_{k}+1-i)-\sum_{i= 2}^{n_{b}}\frac{b_{i}}{b_{1}}u_{\text{ff}}(k+1-i). \tag{8}\]

In order to implement the feedforward controller (8) we require two assumptions:

1. _Reference preview:_ future reference values up to \(r(k+n_{k}+1)\) are known at time instant \(k\);
2. _Stable inverse dynamics:_\(G^{-1}(q)\) is stable, i.e., has poles within the unit circle, or, equivalently, \(G(q)\) is minimum phase. This requires that the zeros of \(\sum_{i=1}^{n_{b}}b_{i}q^{n_{b}-i}\) are within the unit circle.

**Remark II.2**: _When \(G(q)\) is nonminimum phase, it is common practice to either approximate \(G(q)\) using a minimum phase transfer function or to use non-causal filtering to obtain a stable exact inverse of \(G(q)\), see, e.g., [12] for an overview._

### _Nonlinear feedforward controller design_

The linear feedforward control framework can be extended to deal with the nonlinearities \(g(\cdot)\), by starting from a nonlinear correspondent of (7), i.e., a nonlinear autoregressive exogeneous (NARX) representation, which is given as

\[y(k)=h\big{(}[y(k-1),\ldots,y(k-n_{a}), \tag{9}\] \[u(k-n_{k}-1),\ldots,u(k-n_{k}-n_{b})]^{T}\big{)},\]

where \(h:\mathbb{R}^{n_{a}+n_{b}}\rightarrow\mathbb{R}\) is a nonlinear function that describes the complete dynamics. We assume that there exists an inverse representation of (9), such that, with a slight abuse of notation, we can write

\[u(k)=h^{-1}\big{(}\phi(k)\big{)}, \tag{10}\]

where \(\phi(k):=[y(k+n_{k}+1),\ldots,y(k+n_{k}-n_{a}+1),u(k-1),\ldots,u(k-n_{b}+1)]^{T}\). Following the same reasoning as for the linear feedforward (5), we obtain the ideal nonlinear feedforward controller as

\[u_{\text{ff}}(k)=h^{-1}\big{(}\phi_{\text{ff}}(k)\big{)}, \tag{11}\]

with \(\phi_{\text{ff}}(k):=[r(k+n_{k}+1),\ldots,r(k+n_{k}-n_{a}+1),u_{\text{ff}}(k-1 ),\ldots,u_{\text{ff}}(k-n_{b}+1)]^{T}\).

In general, the function \(h^{-1}\) in (11) is unknown. Moreover, the precision demanded by industry exceeds manufacturing tolerances, which implies that a machine-specific \(h^{-1}\) needs to be found. For these reasons, a systematic data-based approach for finding a model of the inverse system dynamics \(h^{-1}\) would be desirable.

### _Identification for feedforward control_

We discuss the three main aspects to perform a direct identification of the inverse dynamics (10): the data set, the model class, and the identification criterion.

_Data set:_ we consider the availability of a data set that is generated on the system displayed in Fig. 1, i.e., satisfying dynamics (9). As a result, with a slight abuse of notation, we have

\[Z^{N}=\{\phi_{0},u_{0},\ldots,\phi_{N-1},u_{N-1}\}, \tag{12}\]

where \(\phi_{i}:=\phi(i)\) and \(u_{i}:=u(i)\) for \(i\in\{0,\ldots,N-1\}\) in (10) during the data generating experiment, where \(N\in\mathbb{Z}_{>0}\) is the number of data points.

Fig. 1: Feedback–feedforward control scheme for a system that is described by a nominal linear model \(G(s)\), possible higher order dynamics \(\sum_{i}G_{i}(s)\), and nonlinear parasitic effects \(g(\cdot)\) which can depend on derivatives of \(u\) and \(y\).
_Model class:_ we require a model parametrization of the inverse dynamics.

**Definition II.1**: _A model parametrization of the inverse dynamics is given as_

\[\hat{u}\big{(}\theta,\phi(k)\big{)}=f\big{(}\theta,\phi(k)\big{)}, \tag{13}\]

_where \(\hat{u}\big{(}\theta,\phi(k)\big{)}\) is the prediction of the input \(u(k)\), \(\theta\in\mathbb{R}^{n_{\theta}}\) denotes the vector of parameters with \(n_{\theta}\in\mathbb{Z}_{>0}\), and \(f:\mathbb{R}^{n_{\theta}}\times\mathbb{R}^{n_{a}+n_{b}}\to\mathbb{R}\) is a user-defined function._

_Identification criterion:_ the identification criterion defines the best choice of parameters \(\theta\) to have the model (13) fit the inverse system dynamics (10). Typically, the identification criterion aims to minimize a cost function, i.e.,

\[\hat{\theta}=\arg\min_{\theta}V(\theta,Z^{N}), \tag{14}\]

such as the mean-squared error (MSE)

\[V(\theta,Z^{N})=V_{\text{MSE}}(\theta,Z^{N}):=\frac{1}{N}\sum_{i=0}^{N-1} \big{(}u_{i}-\hat{u}(\theta,\phi_{i})\big{)}^{2}. \tag{15}\]

Similar to (8) and (11), the feedforward controller is obtained by computing the input \(u(k)\) that yields the output \(y(k)=r(k)\) for the identified model, i.e., (13) with \(\theta=\hat{\theta}\) from (14), such that

\[u_{\text{ff}}(k)=\hat{u}\big{(}\hat{\theta},\phi_{\text{ff}}(k)\big{)}=f\big{(} \hat{\theta},\phi_{\text{ff}}(k)\big{)}. \tag{16}\]

The model parametrization (13) is a crucial choice made by the user. Conventional methods employ a model that is derived from physical knowledge of the system, such that

\[f\big{(}\theta,\phi(k)\big{)}=f_{\text{phy}}\big{(}\theta_{\text{phy}},\phi(k )\big{)}, \tag{17}\]

where \(\theta_{\text{phy}}\in\mathbb{R}^{n_{\theta_{\text{phy}}}}\) are the physical parameters. For reasons that become apparent in Section IV, we denote \(\theta_{\text{phy}}^{*}\) as the parameters of the physical model (17) identified according to the MSE identification criterion (14), (15). As an alternative to physics-based models (17), black-box NNs can be used to parametrize the inverse dynamics (10), i.e.,

\[\hat{u}\big{(}\theta_{\text{NN}},\phi(k)\big{)} =f_{\text{NN}}\big{(}\theta_{\text{NN}},\phi(k)\big{)}\] \[=W_{L+1}\alpha_{L}\Big{(}\ldots\alpha_{1}\big{(}W_{1}\phi(k)+B_{1 }\big{)}\Big{)}+B_{L+1}, \tag{18}\]

where \(\alpha_{l}:\mathbb{R}^{n_{l}}\to\mathbb{R}^{n_{l}}\) denotes the aggregation of activation functions with \(n_{l}\in\mathbb{Z}_{>0}\) the number of neurons in layer \(l\in\{0,\ldots,L\}\), and \(L\in\mathbb{Z}_{>0}\) the number of hidden layers. The parameters \(\theta_{\text{NN}}:=[\text{col}(W_{1})^{T},B_{1}^{T},\ldots,\text{col}(W_{L+1 })^{T},B_{L+1}^{T}]^{T}\) are the concatenation of all weights \(W_{l}\in\mathbb{R}^{n_{l}\times n_{l-1}}\) and biases \(B_{l}\in\mathbb{R}^{n_{l}}\), where \(\text{col}(W_{l})\) stacks the columns of \(W_{l}\). A NN model (18) can approximate any continuous nonlinear mapping on a bounded set with arbitrary accuracy, provided that it consists of a sufficient number of neurons and/or hidden layers [15]. This concludes the required preliminary knowledge to define the problem considered in this paper.

## III Problem statement

To formally characterize the inherent limitation of a physical model (17), we define the corresponding structural model errors.

**Definition III.1**: _Given the physical model (17), we define the structural model error as_

\[g\big{(}\phi(k)\big{)}:=h^{-1}\big{(}\phi(k)\big{)}-f_{\text{phy}}\big{(} \theta_{\text{phy}}^{*},\phi(k)\big{)}. \tag{19}\]

_Notice that Fig. 1 displays the structural model error \(g\big{(}\phi(k)\big{)}\) for the situation in which the forward dynamics are considered and the physical model is linear, i.e., \(f_{\text{phy}}\big{(}\theta_{\text{phy}},\phi(k)\big{)}=\theta_{\text{phy}}^{T} \phi(k)\). The inverse dynamics (10) can be rewritten into_

\[u(k)=f_{\text{phy}}\big{(}\theta_{\text{phy}}^{*},\phi(k)\big{)}+g\big{(}\phi( k)\big{)}. \tag{20}\]

In order to derive conditions for which the feedforward controller (11) achieves zero tracking error, we define a set of operating conditions \(\mathcal{R}\) as follows.

**Definition III.2**: _Denote \(\Phi_{\text{ff}}\subseteq\mathbb{R}^{n_{a}+n_{b}}\) as all regressor points \(\phi_{\text{ff}}(k)\) supplied to the feedforward controller (16) for all references \(r(k)\) and all \(k\). Then, the set of operating conditions \(\mathcal{R}\) is defined as_

\[\mathcal{R}:=\Phi_{\text{ff}}\cup\{\phi_{0},...,\phi_{N-1}\}. \tag{21}\]

_As an example, it is possible to obtain the operating conditions \(\mathcal{R}\) for the moving mass (3) by considering maxima and minima on the position, velocity, and acceleration._

By comparing the inverse dynamics (10) with the feed-forward input (11), we obtain the following condition for perfect tracking

\[u_{\text{ff}}(k)=f\big{(}\hat{\theta},\phi_{\text{ff}}(k)\big{)}=h^{-1}\big{(} \phi_{\text{ff}}(k)\big{)},\quad\forall\ \phi_{\text{ff}}(k)\in\mathcal{R}. \tag{22}\]

Consequently, due to the structural model errors (20), condition (22) cannot be satisfied when using a physical model (17). On the other hand, when using a NN model (18), it is possible to find a sufficient number of hidden layers and neurons per layer in combination with a set of parameters \(\theta_{\text{NN}}^{*}\) such that (22) is satisfied. In practice however, the tedious hyperparameter tuning and non-convex optimization can return suboptimal identified parameters \(\hat{\theta}_{\text{NN}}\) that do not even reach performance of the physics-based feedforward. Additionally, dangerous situations can arise when the data set \(Z^{N}\) does not cover \(\mathcal{R}\) sufficiently well, e.g., as a result of limitations on either time or safety there can be no training data available for all relevant operating conditions. Finally, there is only limited research available for validating stability of NN-based feedforward controllers, see, e.g., [23, 24], and no tools to enforce stability a priori.

To solve the above issues of NN-based feedforward controllers, in this paper, we develop a generalized feedforward control design methodology based on PGNNs. To this end, we introduce the formal definition of a PGNN model class.

**Definition III.3**: _The PGNN model parameterization of the inverse dynamics is given as_

\[\hat{u}\big{(}\theta,\phi(k)\big{)}=f_{\text{phy}}\big{(}\theta_{\text{phy}}, \phi(k)\big{)}+f_{\text{NN}}\big{(}\theta_{\text{NN}},T\big{(}\phi(k)\big{)} \big{)}, \tag{23}\]
where \(\theta=[\theta_{\text{NN}}^{T},\theta_{\text{phy}}^{T}]^{T}\) are the PGNN parameters, and \(T:\mathbb{R}^{n_{a}+n_{b}}\rightarrow\mathbb{R}^{n_{0}}\) is an input transformation, with \(n_{0}\in\mathbb{Z}_{>0}\) the number of NN inputs.

**Remark III.1**: _The input transformation \(T\) can be used to improve the numerical properties of the PGNN training, as well as to create physically relevant inputs. As a practical example, consider a moving mass and let \(\phi_{i}=[y_{i+1},y_{i}]^{T}\) with \(y_{i}\) the measured position. Then, choosing \(T(\phi_{i})=\left[\frac{y_{i+1}-y_{i}}{T_{s}},y_{i}\right]^{T}\) with \(T_{s}\in\mathbb{R}_{>0}\) the sampling time, can improve numerical properties of the training when \(T_{s}\) is small, such that \(y_{i+1}\approx y_{i}\). As a second example, let \(y_{i}\) be the measured rotation of a rotary motor. Then, choosing \(T(\phi_{i})=\left[\frac{y_{i+1}-y_{i}}{T_{s}},\text{mod}(y_{i})\right]^{T}\) with \(\text{mod}(y_{i})\) the remainder after division of \(y_{i}\) by \(2\pi\) can improve the training convergence and extrapolation in terms of the rotation \(y_{i}\), since \(T(\cdot)\) imposes the rotational reproducible behaviour._

**Remark III.2**: _The PGNN (23) imposes physical knowledge by embedding a known physical model within a NN. Note that this is different from physics-informed neural networks (PINNs) [28, 29], which impose physical knowledge via a regularization term in the cost function that penalizes deviations of the NN output with respect to the output of the available physical model. The use of PINNs for feedforward control was explored in, e.g., [30]. It is also worth to point out that, conforming to this definition, the PGNNs developed in [31] are PINNs, i.e., the physical model is not embedded within a NN therein._

## IV Physics-guided neural networks for feedforward control

In this section, we introduce a regularized cost function for training the PGNN (24) to avoid competition between the physical and NN layers. We demonstrate that, under suitable assumptions on the model class, training data, and training convergence, the PGNN recovers the inverse system dynamics. Moreover, we consider the situation in which these assumptions are violated and introduce methods to ensure that \(i)\) the PGNN improves with respect to using only the physical model in terms of data fit via an optimized parameter initialization, and \(ii)\) we impose a form of graceful degradation for the PGNN when operated at conditions that were not present in the training data.

For ease of presentation, in the remainder of this paper we consider \(T\big{(}\phi(k)\big{)}=\phi(k)\), such that the PGNN (23) reduces to

\[\hat{u}\big{(}\theta,\phi(k)\big{)}=f_{\text{phy}}\big{(}\theta_{\text{phy}}, \phi(k)\big{)}+f_{\text{NN}}\big{(}\theta_{\text{NN}},\phi(k)\big{)}. \tag{24}\]

### _Regularized PGNN training and inverse system identification_

The PGNN (24) is generally overparameterized, i.e., the NN can identify parts of the physical model which results in a parameter drift during training. Therefore, we employ a regularized cost function

\[V(\theta,Z^{N})=V_{\text{MSE}}(\theta,Z^{N})+\lambda V_{\text{reg}}(\theta), \tag{25}\]

with \(\lambda>0\), and the regularization cost is

\[V_{\text{reg}}(\theta):=\left\|\begin{bmatrix}\Lambda_{\text{NN}}&0\\ 0&\Lambda_{\text{phy}}\end{bmatrix}\left(\theta-\begin{bmatrix}0\\ \theta_{\text{phy}}^{*}\end{bmatrix}\right)\right\|_{2}^{2}. \tag{26}\]

In (26), \(\Lambda_{\text{NN}},\Lambda_{\text{phy}}\) are matrices that define the relative importance of the different regularization terms. Note that, \(\Lambda_{\text{NN}}\) relates to the standard \(\mathcal{L}_{2}\) regularization for the network weights and biases [26, Chapter 7], while \(\Lambda_{\text{phy}}\) solves the overparameterization of the PGNN model (24).

**Remark IV.1**: _The hyperparameter \(\lambda\) in (25) relates the importance of the data fit with the regularization costs, and can be tuned using, e.g., the L-curve method. The L-curve displays the data fit (\(x\)-axis) and regularization cost (\(y\)-axis) for varying values of \(\lambda\). The optimal choice for \(\lambda\) is the point where the curvature is largest [32]. An example is shown in Fig. 3 with \(\Lambda_{\text{phy}}=\Lambda_{\text{NN}}=I\) using data that is generated on a coreless linear motor, which is discussed in Section VI. The optimal choice is \(\lambda=0.07\). Although the differences in data fit \(V_{\text{MSE}}\) in Fig. 3 seem minor, the minimal data fit is limited by the noise that is present in the data set. Correspondingly, the effect of \(\lambda\) is enlarged when the noise is subtracted._

Fig. 3: L–curve obtained by training the PGNN (24) according to (14), (25) with \(\Lambda_{\text{phy}}=\Lambda_{\text{NN}}=I\) for \(20\) values of \(\lambda\) logarithmically spaced in \([10^{-5},10]\), and optimal choice \(\lambda=0.07\) (red circle).

Fig. 2: Schematic overview of the physics–guided neural network.
The PGNN (24) recovers the inverse dynamics (10) when the following assumptions are satisfied.

**Assumption IV.1**: _There exists a \(\theta_{\text{NN}}^{\ast}\) such that \(f_{\text{NN}}\big{(}\theta_{\text{NN}}^{\ast},\phi(k)\big{)}=g\big{(}\phi(k)\big{)}\) for all \(\phi(k)\in\mathcal{R}\)._

**Assumption IV.2**: _For \(\theta_{\text{NN}}^{A}\neq\theta_{\text{NN}}^{B}\) with \(f_{\text{NN}}\big{(}\theta_{\text{NN}}^{A},\phi(k)\big{)}\neq f_{\text{NN}} \big{(}\theta_{\text{NN}}^{B},\phi(k)\big{)}\) for some \(\phi(k)\in\mathcal{R}\), it holds that_

\[\frac{1}{N}\sum_{i=0}^{N-1}\big{(}f_{\text{NN}}(\theta_{\text{NN}}^{A},\phi_{i })-f_{\text{NN}}(\theta_{\text{NN}}^{B},\phi_{i})\big{)}^{2}>0. \tag{27}\]

**Assumption IV.3**: _The minimization of (25) over \(\theta\) yields a global optimum._

Assumption IV.1 dictates that there should exist a choice of parameters for which the PGNN (24) recovers the original inverse system (10), i.e., the system should be in the model set. Assumption IV.2 relates to the persistence of excitation.

**Proposition IV.1** (PGNN consistency): _Consider the PGNN (24) that is used to identify the inverse dynamics (10) according to identification criterion (14) with cost function (25) using \(\Lambda_{\text{NN}}=0\). Suppose that Assumptions IV.1, IV.2, and IV.3 hold. Then, the identified PGNN parameters satisfy \(\hat{\theta}=[\hat{\theta}_{\text{phy}}^{T},\hat{\theta}_{\text{NN}}^{T}]^{T}= [\theta_{\text{phy}}^{\ast^{T}},\theta_{\text{NN}}^{\ast^{T}}]^{T}\)._

See Appendix I.

From Assumption IV.1 and (20), it is observed that the identified PGNN recovers the inverse dynamics for all operating conditions, i.e.,

\[f_{\text{phy}}\big{(}\hat{\theta}_{\text{phy}},\phi(k)\big{)}+f_{\text{NN}} \big{(}\hat{\theta}_{\text{NN}},\phi(k)\big{)}=h^{-1}\big{(}\phi(k)\big{)},\ \forall\ \phi(k)\in\mathcal{R}, \tag{28}\]

such that the zero tracking error condition (22) is satisfied.

**Remark IV.2**: _A parallel linear-NN model structure was also employed for feedforward control design in [33]. Therein, an alternative regularization method based on orthogonal projection was developed to avoid the competition between the linear and NN layers under the assumption that \(g(\cdot)\) is identically zero outside \(\mathcal{R}\)._

### _Optimized PGNN parameter selection_

When Assumption IV.1 or IV.3 is violated, recovering the inverse dynamics is not guaranteed, or the training does not find the optimal parameters. The training, i.e., the optimization over \(\theta\) of (25), typically employs a nonlinear optimization scheme, that attains parameter estimates \(\theta^{(j)}\) where \(j\in\{0,\ldots,J\}\), with \(J\in\mathbb{Z}_{\geq 0}\) the number of solver iterations. The estimated parameter vector is chosen as

\[\hat{\theta}=\theta^{(j)},\quad j=\text{arg}\min_{j\in\{0,\ldots,J\}}V\big{(} \theta^{(j)},Z^{N}\big{)}. \tag{29}\]

In what follows, we will derive a specific choice for the parameters \(\theta^{(j)}\) which achieve a smaller value of the cost function (25) compared to using only the physical model. For simplicity, we consider a LIP physical model (17), i.e., \(f_{\text{phy}}\big{(}\theta_{\text{phy}},\phi(k)\big{)}=\theta_{\text{phy}}^{T} T_{\text{phy}}\big{(}\phi(k)\big{)}\), and rewrite the PGNN (24) according to

\[\hat{u}\big{(}\theta,\phi(k)\big{)}=\theta_{\text{L}}^{T}\phi_{\text{L}}\big{(} \theta_{\text{NL}},\phi(k)\big{)}:=\theta_{\text{L}}^{T}\begin{bmatrix}\alpha_ {l}\big{(}\theta_{\text{NL}},\phi(k)\big{)}\\ 1\\ T_{\text{phy}}\big{(}\phi(k)\big{)}\end{bmatrix}, \tag{30}\]

where \(\theta_{\text{L}}:=[\text{col}(W_{L+1})^{T},B_{L+1},\theta_{\text{phy}}^{T}]^{T}\) are the parameters in which the PGNN (30) is linear, and \(\theta_{\text{NL}}:=[\text{col}(W_{1})^{T},B_{1}^{T},\ldots,\text{col}(W_{L})^ {T},B_{L}^{T}]^{T}\), such that \(\theta=[\theta_{\text{NL}}^{T},\theta_{\text{L}}^{T}]^{T}\). From (30), we observe that an equivalent of the physical model (17) is obtained by selecting \(\theta_{\text{L}}=\overline{\theta}_{\text{L}}:=[0,0,\theta_{\text{phy}}^{\ast ^{T}}]^{T}\). We rewrite the regularization term (26) into

\[V_{\text{reg}}(\theta)=\left\|\begin{bmatrix}\Lambda_{\text{NL}}&\Lambda\\ \Lambda^{T}&\Lambda_{\text{L}}\end{bmatrix}\left(\begin{bmatrix}\theta_{\text{ NL}}\\ \theta_{\text{L}}\end{bmatrix}-\begin{bmatrix}0\\ \overline{\theta}_{\text{L}}\end{bmatrix}\right)\right\|_{2}^{2}, \tag{31}\]

where \(\Lambda\), \(\Lambda_{\text{NL}}^{T}=\Lambda_{\text{NL}}\), and \(\Lambda_{\text{L}}^{T}=\Lambda_{\text{L}}\) are obtained by selecting rows and columns of \(\Lambda_{\text{NN}}\) and \(\Lambda_{\text{phy}}\) accordingly. Then, given any parameter set \(\theta_{\text{NL}}^{(j)}\), choose \(\theta_{\text{L}}^{(j)}\) according to

\[\begin{split}\theta_{\text{L}}^{(j)}=& M(\theta_{\text{NL}}^{(j)})^{-1}\Bigg{(}\frac{1}{N}\sum_{i=0}^{N-1}u_{i}\phi_{ \text{L}}(\theta_{\text{NL}}^{(j)},\phi_{i})\\ &+\lambda\Big{(}\big{(}\Lambda_{\text{L}}^{2}+\Lambda^{T}\Lambda\big{)} \overline{\theta}_{\text{L}}-\big{(}\Lambda^{T}\Lambda_{\text{NL}}+\Lambda_{ \text{L}}\Lambda^{T}\big{)}\theta_{\text{NL}}^{(j)}\Big{)}\Bigg{)},\end{split} \tag{32}\]

where \(M(\theta_{\text{NL}}^{(j)})\) is given as

\[\begin{split} M(\theta_{\text{NL}}^{(j)}):=& \frac{1}{N}\sum_{i=0}^{N-1}\phi_{\text{L}}(\theta_{\text{NL}}^{(j)},\phi_{i}) \phi_{\text{L}}(\theta_{\text{NL}}^{(j)},\phi_{i})^{T}\\ &+\lambda\big{(}\Lambda_{\text{L}}^{2}+\Lambda^{T}\Lambda\big{)}. \end{split} \tag{33}\]

Note that (32) is valid only if \(M(\theta_{\text{NL}}^{(j)})\) is nonsingular, which yields a unique solution. Nonsingularity of \(M(\theta_{\text{NL}}^{(j)})\) is also referred to as the persistence of excitation.

**Proposition IV.2** (PGNN improves over physics): _Consider the PGNN (30) with \(\theta_{\text{NL}}^{(j)}\) given, e.g., initialized randomly for \(j=0\) or attained during training for \(j\neq 0\). Suppose that \(M(\theta_{\text{NL}}^{(j)})\) is nonsingular, and choose \(\theta_{\text{L}}^{(j)}\) according to (32). Then, for the cost function (25), we have_

\[V\left(\begin{bmatrix}\theta_{\text{NL}}^{(j)}\\ \theta_{\text{L}}^{(j)}\end{bmatrix},Z^{N}\right)\leq V\left(\begin{bmatrix} \theta_{\text{NL}}^{(j)}\\ \overline{\theta}_{\text{L}}\end{bmatrix},Z^{N}\right), \tag{34}\]

_with strict inequality if and only if_

\[\begin{split}\frac{1}{N}\sum_{i=0}^{N-1}\Big{(}u_{i}\phi_{ \text{L}}&\big{(}\theta_{\text{NL}}^{(j)},\phi_{i}\big{)}-\phi_{\text{L}} \big{(}\theta_{\text{NL}}^{(j)},\phi_{i}\big{)}\phi_{\text{L}}\big{(}\theta_{ \text{NL}}^{(j)},\phi_{i}\big{)}^{T}\overline{\theta}_{\text{L}}\Big{)}\\ &-\lambda(\Lambda^{T}\Lambda_{\text{NL}}+\Lambda_{\text{L}}\Lambda^{T} )\theta_{\text{NL}}^{(j)}\neq 0.\end{split} \tag{35}\]

_Moreover, if \(\Lambda_{\text{NN}}\) is such that \(\Lambda=0\), i.e., the cross terms between \(\theta_{\text{NL}}\) and \([W_{L+1}^{T},B_{L+1}]^{T}\) are not penalized, the PGNN achieves a better data fit compared to the physical model, i.e.,_

\[V_{\text{MSE}}\left(\begin{bmatrix}\theta_{\text{NL}}^{(j)}\\ \theta_{\text{L}}^{(j)}\end{bmatrix},Z^{N}\right)\leq V_{\text{MSE}}(\theta_{ \text{phy}}^{\ast},Z^{N}), \tag{36}\]

_with strict inequality if (35) holds._

See Appendix II.

Condition (35) is obtained by rewriting \(\theta_{\text{L}}^
between the output of the final hidden layer \(\alpha_{L}(\cdot)\) and the structural model errors present in the data \(g(\phi_{i})\).

**Remark IV.3**: _Proposition IV.2 can be directly extended for PGNNs with physical models that are not LIP. To see this, rewrite the physical model_

\[f_{\text{phy}}\big{(}\theta_{\text{phy}},\phi(k)\big{)}=\theta_{\text{L-phy}}^{T} T_{\text{phy}}\big{(}\theta_{\text{NL-phy}},\phi(k)\big{)}, \tag{37}\]

_and choose \(\theta_{\text{L}}^{(j)}=[\text{col}(W_{L+1}^{(j)})^{T},B_{L+1}^{(j)},{\theta_{ \text{L-phy}}^{(j)}}^{T}]^{T}\) according to (32) using \(\theta_{\text{NL}}^{(j)}=[\text{col}(W_{1}^{(j)})^{T},B_{1}^{(j)}{}^{T},..., \text{col}(W_{L}^{(j)})^{T},B_{L}^{(j)}{}^{T},\theta_{\text{NL-phy}}^{T}]^{T}\)._

In practice, the optimized parameter selection is used after training, i.e., update \(\theta_{\text{L}}^{(j)}\) for \(j=J\), during training for each \(j=\{0,...,J\}\), or as an initialization for \(j=0\). Note that, from (29) and Proposition IV.2, if (35) holds, we have that

\[V(\hat{\theta},Z^{N})\leq V(\theta^{(0)},Z^{N})<V(\overline{\theta},Z^{N}), \tag{38}\]

when (32) is used for initialization of \(\theta_{\text{L}}^{(0)}\).

### _Enhancing PGNN robustness via regularized training_

In general, there is no systematic method to validate Assumption IV.2 in practice. For this reason, [34] gives as a user guideline to sample the complete domain of interest, i.e., the operating conditions \(\mathcal{R}\). For reasons of time and safety, it can be infeasible to design an experiment that achieves this. Correspondingly, we require a form of graceful degradation for situations in which the PGNN is operated on conditions that were not present in the training data. Since there is no data available for these conditions, the physical model is the only source of information. Therefore, we aim to have the PGNN comply with the physical model via a regularization term in the cost function, such that (25) becomes

\[V(\theta,Z^{N})=V_{\text{MSE}}(\theta,Z^{N})+\lambda V_{\text{reg}}(\theta)+ \gamma V_{\text{phy}}(\theta,Z^{E}). \tag{39}\]

In (39), \(\gamma\in\mathbb{R}_{>0}\) is a regularization parameter, and the regularization cost is

\[V_{\text{phy}}(\theta,Z^{E}):=\frac{1}{E}\sum_{i=0}^{E-1}\big{(}f_{\text{phy} }(\theta_{\text{phy}}^{*},\phi_{i}^{E})-\hat{u}(\theta,\phi_{i}^{E})\big{)}^{2}. \tag{40}\]

The set \(Z^{E}=\{\phi_{0}^{E},\ldots,\phi_{E-1}^{E}\}\) contains user-defined regressor points for which we aim to have the PGNN comply with the physical model, e.g., due to the absence of reliable data in \(Z^{N}\). Thereby, the aim is to have \(Z^{N}\cup Z^{E}\) cover the operating conditions \(\mathcal{R}\) up to high accuracy, which is done automatically by Algorithm 1. In Algorithm 1, \(C(\zeta,Z^{N},Z^{E})\) is an objective function that specifies the goal of the optimization. For example, choosing \(\phi_{i}^{E}\) to maximize the minimum squared Euclidean distance with respect to all regressor points, is achieved by choosing

\[C(\zeta,Z^{N},Z^{E}):=\min_{\phi\in Z^{N},Z^{E}}\left\|\phi-\zeta\right\|_{2}^ {2}. \tag{41}\]

Algorithm 1 iterates until a stopping criterion is met, e.g., by fixing a maximum number of points, or a minimum threshold \(\epsilon\) for the objective function (41).

**Remark IV.4**: _The results of Proposition IV.2 can be extended directly to the cost function (39) by appropriately revising the computation of \(\theta_{\text{L}}^{(j)}\) in (32), i.e., compute the least squares solution of (39) instead of (25)._

Consider a mechanical system as in (3), where the operating conditions \(\mathcal{R}\) are represented by a maximum position of \(0.15\)\(m\) and velocity of \(0.2\)\(\frac{m}{s}\), and neglect the acceleration for the sake of simplicity. Since it was deemed unsafe to travel the full stroke during the data generating experiment, the training data \(Z^{N}\) describes \(\mathcal{R}\) only in a limited range of the position, see the blue crosses in Fig. 4. Consequently, application of Algorithm 1 generates the set \(Z^{E}\), see the orange crosses in Fig. 4, such that \(Z^{N}\cup Z^{E}\) cover the complete operating conditions \(\mathcal{R}\).

```
Initialize\(Z^{E}=\{\}\), \(i=0\), while\(i<E-1\)\(\wedge\)\(C\big{(}\phi(i-1),Z^{N},Z^{E}\big{)}>\epsilon\)do \(\phi_{i}^{E}=\text{arg}\max_{\zeta\in\mathcal{R}}C\big{(}\zeta,Z^{N},Z^{E} \big{)}\), \(Z^{E}=Z^{E}\cup\phi_{i}^{E}\), \(i=i+1\). endwhile
```

**Algorithm 1** Design algorithm for \(Z^{E}\).

## V ISS of PGNN-based feedforward controllers

Stability of a feedforward controller (16) is a prerequisite for safe operation of the closed-loop system, i.e., for bounded reference values \(r(k)\) the feedforward input \(u_{\text{ff}}(k)\) should remain bounded. For linear feedforward controllers, stability is tested by ensuring that the poles of the feedforward controller transfer function are within the unit circle. However, the nonlinear NN in the PGNN feedforward (16), (24) complicates the stability assessment. For the sake of presentation, we consider a PGNN (24) with linear physical model, i.e., \(f_{\text{phy}}\big{(}\theta_{\text{phy}},\phi(k)\big{)}=\theta_{\text{phy}}^{T} \phi(k)\), and highlight some further generalizations. Next, to facilitate the stability analysis of PGNN-based feedforward controllers,

Fig. 4: Two–dimensional illustration of the regressor points in the CLM data set \(Z^{N}\) discussed in Section VI, and the regressor points \(Z^{E}\) designed following Algorithm 1.
we define

\[\hat{a}:=[\hat{a}_{1},\ldots,\hat{a}_{n_{a}+1}]^{T},\quad\hat{b}:=[\hat{b}_{1}, \ldots,\hat{b}_{n_{b}-1}]^{T}, \tag{42}\]

such that \(\hat{\theta}_{\text{phy}}=[\hat{a}^{T},\hat{b}^{T}]^{T}\). Moreover, we define

\[\phi_{r}(k) :=[r(k+n_{k}+1),\ldots,r(k+n_{k}-n_{a}+1)]^{T},\] \[\phi_{u_{\text{ff}}}(k) :=[u_{\text{ff}}(k-1),\ldots,u_{\text{ff}}(k-n_{b}+1)]^{T}, \tag{43}\]

such that \(\phi_{\text{ff}}(k)=[\phi_{r}(k)^{T},\phi_{u_{\text{ff}}}(k)^{T}]^{T}\). Note that, for a linear feedforward controller (8), stability is determined only by \(\hat{b}\). Then, by choosing \(x(k):=\phi_{u_{\text{ff}}}(k)\), we obtain a state-space representation of the PGNN (24) as

\[x(k+1)= A(\hat{b})x(k) \tag{44}\] \[+B\left(\hat{a}^{T}\phi_{r}(k)+f_{\text{NN}}\left(\hat{\theta}_{ \text{NN}},\begin{bmatrix}\phi_{r}(k)\\ x(k)\end{bmatrix}\right)\right),\]

where \(B=[1,O]^{T}\in\mathbb{R}^{n_{b}-1}\), \(A(\hat{b})=\begin{bmatrix}\hat{b}^{T}\\ I\quad O\end{bmatrix}\in\mathbb{R}^{(n_{b}-1)\times(n_{b}-1)}\), and \(\phi_{r}(k)\) is the external input.

Because the state-space PGNN feedforward (44) has \(\phi_{r}(k)\) as an exogeneous input, it makes sense to pursue ISS guarantees [35], which is a stronger property than nominal stability and guarantees boundedness of the state \(x(k)\) for any bounded reference \(\phi_{r}(k)\). First, we provide a general ISS result for the PGNN state-space representation (44). Afterwards, we show how ISS of the resulting PGNN feedforward controllers can be imposed a priori during the PGNN training by adjusting the identification criterion.

We recall the following notions: a function \(\kappa:\mathbb{R}_{\geq 0}\rightarrow\mathbb{R}_{\geq 0}\) is a \(\mathcal{K}\)-function if it is continuous, strictly increasing and \(\kappa(0)=0\). It is a \(\mathcal{K}_{\infty}\)-function if it is a \(\mathcal{K}\)-function and \(\kappa(s)\rightarrow\infty\) as \(s\rightarrow\infty\). A function \(\eta:\mathbb{R}_{\geq 0}\times\mathbb{R}_{\geq 0}\) is a \(\mathcal{KL}\)-function if, for each fixed \(k\geq 0\), the function \(\eta(\cdot,k)\) is a \(\mathcal{K}\)-function and for each fixed \(s\geq 0\), the function \(\eta(s,\cdot)\) is decreasing and \(\eta(s,k)\to 0\) as \(k\rightarrow\infty\).

**Definition V.1** (Iss [35]): _A system \(x(k+1)=f\big{(}x(k),\phi_{r}(k+1)\big{)}\) is ISS if there exist a \(\mathcal{KL}\)-function \(\eta:\mathbb{R}_{\geq 0}\times\mathbb{R}_{\geq 0}\rightarrow\mathbb{R}_{ \geq 0}\) and a \(\mathcal{K}\)-function \(\kappa\) such that, for each finite input sequence \(\phi_{r}\) and each initial state \(x(0)=\xi\in\mathbb{R}^{n_{b}-1}\), it holds that_

\[|x(k)|\leq\eta\big{(}|\xi|,k\big{)}+\kappa\big{(}\|\phi_{r}\|\big{)} \tag{45}\]

_for each \(k\in\mathbb{Z}_{+}\)._

**Remark V.1**: _If the state-space PGNN representation (44) is ISS with respect to the external input \(\phi_{r}(k+1)\), we have that:_

1. \(x(k)\) _remains bounded for bounded_ \(\phi_{r}(k)\) _and, consequently,_ \(u_{\text{ff}}(k)\) _remains bounded;_
2. \(x(k)\to 0\) _for_ \(\phi_{r}(k)\to 0\) _and, consequently,_ \(u_{\text{ff}}(k)\to 0\)_._

In order to analyze ISS of the PGNN feedforward controller (44), we will use a continuous ISS-Lyapunov function, which implies ISS as shown in [35].

**Definition V.2** (Iss-Lyapunov [35]): _A continuous function \(V:\mathbb{R}^{n}\rightarrow\mathbb{R}_{\geq 0}\) is called an ISS-Lyapunov function for a discrete-time state-space system \(x(k+1)=f\big{(}x(k),\phi_{r}(k)\big{)}\) if the following conditions hold:_

1. _There exists_ \(\mathcal{K}_{\infty}\)_-functions_ \(\kappa_{1}\)_,_ \(\kappa_{2}\) _such that_ \[\kappa_{1}\big{(}|x(k)|\big{)}\leq V\big{(}x(k)\big{)}\leq\kappa_{2}\big{(}|x( k)|\big{)},\forall\quad x(k)\in\mathbb{R}^{n_{b}-1}.\] (46)
2. _There exists a_ \(\mathcal{K}_{\infty}\)_-function_ \(\kappa_{3}\) _and a_ \(\mathcal{K}\)_-function_ \(\sigma\)_, such that_ \[V\big{(}x(k+1)\big{)}-V\big{(}x(k)\big{)}\leq-\kappa_{3}(|x(k)|)+\sigma(|\phi_ {r}(k+1)|),\] (47) \[\forall\;x(k)\in\mathbb{R}^{n_{b}-1},\;\forall\;\phi_{r}(k)\in \mathbb{R}^{n_{b}+1}.\]

**Remark V.2**: _For all commonly applied activation functions, the NN (18) is globally Lipschitz, i.e., there exists an \(L:=[L_{r}^{T},L_{u_{\text{ff}}}^{T}]^{T}\) such that_

\[\left|f_{\text{NN}}\big{(}\hat{\theta}_{\text{NN}},\phi_{\text{ ff}}^{A}(k)\big{)}-f_{\text{NN}}\big{(}\hat{\theta}_{\text{NN}},\phi_{\text{ff}}^{B}(k) \big{)}\right| \tag{48}\] \[\leq L^{T}\left|\phi_{\text{ff}}^{A}(k)-\phi_{\text{ff}}^{B}(k) \right|.\]

_Using backpropagation, it is possible to find values for \(L\), e.g.,_

\[L^{T} =\max_{\phi_{\text{ff}}(k)}\frac{\partial f_{\text{NN}}\big{(} \hat{\theta}_{\text{NN}},\phi_{\text{ff}}(k)\big{)}}{\partial\phi_{\text{ff}}(k)} \tag{49}\] \[=\max_{\phi_{\text{ff}}(k)}\hat{W}_{l+1}\text{diag}(\alpha_{l}^{ \prime})\ldots\hat{W}_{2}\text{diag}(\alpha_{1}^{\prime})\hat{W}_{1}\leq\Pi_{i=1} ^{l+1}|\hat{W}_{i}|,\]

_where \(\alpha_{l}^{\prime}=\frac{\partial\alpha_{l}(x)}{\partial x}\big{|}_{x=x_{l}}\), with \(x_{l}\) the input to NN layer \(l\). The upperbound in (49) is obtained by substitution of \(\text{diag}(\alpha_{l}^{\prime})=I\), which holds for \(\tanh\), ReLU, and several other activation functions._

**Assumption V.1**: _The origin \(x(k)=0\) is an equilibrium for the unexcited PGNN (44), i.e., \(\phi_{r}(k)=0\), such that (44) gives_

\[f_{\text{NN}}\big{(}\hat{\theta}_{\text{NN}},0\big{)}=0. \tag{50}\]

_Note that, we can introduce a coordinate transformation, e.g., \(\zeta(k)=x(k)+\varepsilon\) to satisfy Assumption V.1._

**Assumption V.2**: _There exists a \(P\succ 0\) such that \(Q:=P-A(\hat{b})^{T}PA(\hat{b})\succ 0\)._

Assumption V.2 requires that the parameters \(\hat{b}\) are such that the physical model is stable, i.e., \(A(\hat{b})\) is a Schur matrix. A \((P,Q)\) pair satisfying Assumption V.2 is typically obtained by choosing \(Q\succ 0\) and solving the discrete-time Lyapunov equation to obtain \(P\succ 0\). We denote \(\lambda_{\text{min}}(Q)\) and \(\lambda_{\text{max}}(Q)\) as the smallest and largest eigenvalue of \(Q\), respectively.

**Theorem V.1** (PGNN feedforward ISS): _Consider the PGNN feedforward (16), (24) with linear physical model, and its state-space representation (44). Suppose that Assumptions V.1 and V.2 hold. Let \((P,Q)\) satisfy Assumption V.2, and define_

\[c_{\beta}:=B^{T}P\left(I+\frac{1}{\beta\lambda_{\text{min}}(Q)}AA^{T}P\right)B. \tag{51}\]

_Then, if there exists a \(\beta>0\) such that_

\[L_{u_{\text{ff}}}^{T}L_{u_{\text{ff}}}<\frac{(1-\beta)\lambda_{\text{min}}(Q)}{c_ {\beta}}, \tag{52}\]
the PGNN feedforward state-space representation (44) is ISS.

See Appendix III.

**Remark V.3**: _Most NN stability research is conducted for systems operating in closed-loop with a NN-based feedback controller, see, e.g., [24, 36] which provide local stability conditions via linearization and integral quadratic constraints, respectively. A global ISS result for a black box NN with no physical model embedded is presented in [23]._

**Remark V.4**: _The value of \(\beta\) for which the right-hand side in (52) is maximal, is_

\[\begin{split}\beta=&\frac{1}{\lambda_{\text{ min}}(Q)B^{T}PB}\Big{(}-B^{T} PAA^{T}PB+\\ &\sqrt{B^{T}P\big{(}\lambda_{\text{min}}(Q)I+AA^{T}P\big{)}BB^{T} PAA^{T}PB}\Big{)}.\end{split} \tag{53}\]

Note that, since \(B\) is a column, the square root and division are scalar operations. Eq. (53) is obtained by setting the derivative w.r.t. \(\beta\) of the right hand side of (52) equal to zero, and choosing the option for which \(\beta>0\).

**Remark V.5**: _A similar result can be obtained for a PGNN with a general nonlinear physical model, when a quadratic Lyapunov function is available for the physical model._

The ISS condition (52) in Theorem V.1 can be validated _after_ training by using \(\beta\) in (53), the upperbound \(L_{u_{\text{ff}}}=\begin{bmatrix}O,&I\end{bmatrix}L\) in (49) for some \((P,Q)\) pair.

Recall that \(\theta=[\theta_{\text{NN}}^{T},\theta_{\text{phy}}^{T}]^{T}\), \(\theta_{\text{phy}}=[a^{T},b^{T}]^{T}\) for the linear physical model, and \(\theta_{\text{phy}}^{*}\) are the parameters \(\theta_{\text{phy}}\) identified using only the physical model (17) with MSE cost function (15). With the aim to ensure _before_ training that the PGNN is ISS, we fix \(\hat{b}=b^{*}\), and constrain the network weights to satisfy the ISS condition (52).

**Lemma V.1** (Training imposed ISS - stable physical model): _Consider the PGNN feedforward controller with linear physical model, such that it admits a state-space representation of the form (44). Suppose that Assumption V.1 holds, that a \((P,Q)\) pair satisfying Assumption V.2 is available, and choose \(\beta\) as in (53). Define the set_

\[\begin{split}\Theta&:=\Bigg{\{}\theta\in\mathbb{R}^{n _{\theta}}\;\bigg{|}\;\big{(}b=b^{*}\big{)}\;\wedge\\ &\qquad\left(\left\|\;\big{(}\Pi_{l=1}^{L}|W_{l}|\big{)}\;\begin{bmatrix} 0\\ I\end{bmatrix}\;\right\|_{2}^{2}<\frac{(1-\beta)\lambda_{\text{min}}(Q)}{c_{ \beta}}\right)\Bigg{\}},\end{split} \tag{54}\]

_and train the PGNN according to identification criterion_

\[\hat{\theta}=\text{arg}\min_{\theta\in\Theta}\;V\big{(}\theta,Z^{N}\big{)}. \tag{55}\]

_Then, the training returns an ISS PGNN feedforward controller (16), (24)._

See Appendix IV.

Assumption V.2 is violated when \(A(\hat{b})\) is not Schur, i.e., the system is nonminimum phase. There are two main approaches to obtain a stable feedforward controller for linear systems, see, e.g., [12]:

1. Non-causal feedforward controller design;
2. Stable approximate inversion techniques, such as ZPETC, ZMETC, NPZ-ignore.

In the remainder of this section, we demonstrate approaches to assess and impose ISS of PGNN feedforward controllers for nonminimum phase systems.

The nonminimum phase behaviour can be dealt with by extending the preview window of the PGNN, i.e., we adjust the PGNN (24) to

\[\begin{split}\hat{u}\big{(}\theta,\tilde{\phi}(k)\big{)}=[a^{T},b^{T}]\begin{bmatrix}\tilde{\phi}_{y}(k)\\ \tilde{\phi}_{u}(k)\end{bmatrix}+f_{\text{NN}}\left(\theta_{\text{NN}}, \begin{bmatrix}\tilde{\phi}_{y}(k)\\ \tilde{\phi}_{u}(k)\end{bmatrix}\right),\\ \tilde{\phi}_{y}(k):=[r(k+n_{k}+n_{\text{pw}}+1),...,r(k+n_{k}-n_{a}+1)]^{T},\\ \tilde{\phi}_{u}(k):=[u(k-1),...,u(k-n_{b}+n_{\text{us}}+1)]^{T}, \end{split} \tag{56}\]

where \(n_{\text{us}}\in\mathbb{Z}_{>0}\) are the number of unstable eigenvalues of \(A(\hat{b})\) of the original PGNN (24), and \(n_{\text{pw}}\in\mathbb{Z}_{>0}\) is the number of extended preview samples. The non-causal PGNN (56) is then trained following Lemma V.1.

**Remark V.6**: _In order to satisfy Assumption V.2, the parameters \(b^{*}\) should be such that \(A(b^{*})\) is Schur. This is achieved either by choosing \(n_{\text{pw}}\) large enough, such that the identification of \(\theta_{\text{phy}}^{*}=[a^{*T},b^{*T}]^{T}\) returns stable parameters \(b^{*}\), or by fixing \(b^{*}\) to retain the stable eigenvalues of the original \(A(b^{*})\)._

As an alternative approach, it is possible to employ linear stable approximation inversion tools on the linear part of the dynamics. In order to see this, we rewrite the state-space PGNN (44) into

\[u_{\text{ff}}(k)=G(\hat{b},q^{-1})\left(\hat{a}^{T}\phi_{r}(k)+f_{\text{NN}} \left(\hat{\theta},\begin{bmatrix}\phi_{r}(k)\\ \phi_{u_{\text{ff}}}(k)\end{bmatrix}\right)\right), \tag{57}\]

where \(G(\hat{b},q^{-1}):=\frac{1}{1-\sum_{i=1}^{n_{b}}b_{i}q^{-i}}\). Let \(\tilde{G}(\tilde{b},q^{-1}):=\frac{\sum_{i=0}^{n_{b}}\tilde{b}_{i}^{T}q^{-i}}{1- \sum_{i=1}^{n_{b}^{\text{ff}}}b_{i}^{T}q^{-i}}\) be a stable approximate of \(G(\hat{b},q^{-1})\) obtained via a user-preferred method which gives \(\tilde{n}_{b}^{n},\tilde{n}_{b}^{n}\in\mathbb{Z}_{\geq 0}\) and \(\tilde{b}_{i}^{n},\tilde{b}_{i}^{d}\in\mathbb{R}\). Then, a PGNN with stable physical model is obtained from (57) as

\[u_{\text{ff}}(k)=\tilde{G}(\tilde{b},q^{-1})\left(\hat{a}^{T}\phi_{r}(k)+f_{ \text{NN}}\left(\hat{\theta},\begin{bmatrix}\phi_{r}(k)\\ \phi_{u_{\text{ff}}}(k)\end{bmatrix}\right)\right). \tag{58}\]

If the stable approximation method creates terms with \(q^{d}\), \(d\in\mathbb{Z}_{>0}\) in the numerator of \(G(\tilde{b},q^{-1})\) in (58), \(\phi_{u_{\text{ff}}}(k)\) in the NN can only contain feedforward inputs up until \(u_{\text{ff}}(k-d-1)\) in order to compute the feedforward according (58).

ISS of the stable approximated PGNN (58) is imposed by training the unstable PGNN (57) according a revised version of Lemma V.1. To do so, we rewrite the stably inverted PGNN (58) into state-space form using \(x(k):=[u_{\text{ff}}(k),...,u_{\text{ff}}(k-n_{b}-\tilde{n}_{b}^{n}+1)]^{T}\), and use the triangular inequality with Lipschitz bound (48) to construct \(L_{r+}\) and \(L_{x}\) such that

\[\left|\sum_{i=0}^{\tilde{n}_{b}^{n}}\tilde{b}_{i}^{n}f_{\text{NN}}\big{(}\hat{ \theta}_{\text{NN}},\phi_{\text{ff}}(k-i)\big{)}\right|\leq[L_{r+}^{T},L_{x}^{T}] \begin{bmatrix}\phi_{r+}(k)\\ x(k)\end{bmatrix}, \tag{59}\]
with \(\phi_{r+}(k)=[r(k+n_{k}+1),...,r(k+n_{k}-n_{a}-\tilde{n}_{b}^{n}+1)]^{T}\).

**Lemma V.2** (Training imposed ISS - unstable physical model): _Consider the PGNN (24) with linear physical model. Suppose that Assumption V.1 holds, and choose \(\beta\) as in (53). Define_

\[\begin{split}\Theta:=\Bigg{\{}&\theta\in\mathbb{R} ^{n_{\theta}}\Biggm{|}\big{(}b=b^{*}\big{)}\\ &\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad
3. _PGNN-based feedforward with graceful degradation_ with NN dimensions \(n_{1}=16\), \(l=1\) (18), using physical model (61) identified according to (14), (39) with \(\lambda=0.07\), \(\Lambda_{\text{phy}}=\Lambda_{\text{NN}}=I\), and \(\gamma=0.1\).

_Results:_ Fig. 7 the tracking error \(e(k)\) resulting achieved by the different feedforward controllers for \(r_{1}(k)\) in Fig. 6. In order to test robustness of the PGNN, we compute the mean-absolute error (MAE)

\[MAE=\frac{1}{N}\sum_{k=0}^{N-1}|e(k)|, \tag{63}\]

for two references \(r_{1}(k)\) and \(r_{2}(k)\) in Fig. 6 using different velocities \(\max\big{(}\dot{r}_{i}(k)\big{)}=n\,0.025\ \frac{m}{s}\), \(i=\{1,2\}\), \(n=\{1,...,6\}\). Note that, since \(\max\big{(}r_{2}(k)\big{)}=0.15\ m\), the feedforward requires to extrapolate for \(r_{2}(k)\), see Fig. 4. Fig. 8 displays the resulting MAE values, and Table I lists the MAE for \(r_{1}(k)\).

For \(r_{1}(k)\), both the PGNNs outperform the physics-based feedforward in terms of the MAE for all tested velocities. A significant drop in performance is observed when using the PGNN with \(\gamma=0\) on reference \(r_{2}(k)\), which is caused by the PGNN extrapolating badly for \(r(k)>0.1\ m\). The physics-based regularization term (40) prevents this, as seen by the MAE of the PGNN with \(\gamma=0.1\).

### _Nonminimum phase rotating-translating mass_

_System dynamics:_ we consider a translating-rotating mass with force input \(u(k)\) and position output \(y(k)\) at opposite sides of the centre of mass, see Fig. 9. The dynamics are given as

\[\begin{split} y(t)=&\left(\frac{1}{ms^{2}+f_{v}s}- \frac{l_{y}^{2}}{Ms^{2}+2l_{x}ds+2l_{x}k}\right)\\ &\left(u(t)-g\big{(}y(t)\big{)}\right),\\ g\big{(}y(t)\big{)}=& c\sin\left(\frac{2\pi}{l_{m}}y(t )\right),\end{split} \tag{64}\]

where \(l_{x},l_{y}\in\mathbb{R}_{\geq 0}\) are the width and height of the mass \(m\in\mathbb{R}_{>0}\), \(M=\frac{1}{5}m(l_{x}^{2}+l_{y}^{2})\) is the moment of intertia, \(f_{v}\in\mathbb{R}_{>0}\) the viscous friction, and \(d,k\in\mathbb{R}_{>0}\) the damping and spring constant counteracting rotation at both ends of the mass. The function \(g\big{(}y(t)\big{)}\) is the cogging force and is assumed _unknown_, with \(l_{m}\in\mathbb{R}_{>0}\) the magnet pole pitch and \(c\in\mathbb{R}_{>0}\) the cogging magnitude. Parameter values are listed in Table II. The system (64) is controlled in closed-loop by a ZOH-discretized version of

\[C(s)=5\cdot 10^{3}\frac{s+4\pi}{s+20\pi}, \tag{65}\]

which achieves a \(1.22\ Hz\) bandwidth.

_Training data generation:_ data is generated in closed-loop by sampling \(u(k)\) and \(y(k)\) at a frequency of \(1\ kHz\), while exciting the system with the reference \(r(k)\) in the top window of Fig. 10, and adding a white noise with variance \(50\ N^{2}\) to the input \(u(k)\).

_Feedforward controllers:_ ZOH discretization of the transfer function in (64) gives \(n_{a}=n_{b}=4\), \(n_{k}=0\). The identified feedforward controllers are _unstable_, due to the nonminimum phase transfer function in the system (64). We test the following situations:

1. _No feedforward_, i.e., \(u_{\text{ff}}(k)=0\);
2. _Physics-based feedforward_ with ZPETC stable in

\begin{table}
\begin{tabular}{|c|c|c|c|c|c|c|} \hline \(m\) & \(l_{x},l_{y}\) & \(M\) & \(f_{v}\) & \(k\) & \(d\) & \(l_{m}\) & \(c\) \\ \hline \hline \(20\) & \(1\) & \(\frac{40}{3}\) & \(50\) & \(25\cdot 10^{3}\) & \(\frac{575}{3}\) & \(0.05\) & \(1\) \\ \hline \(kg\) & \(m\) & \(kpm^{2}\) & \(\frac{kg}{s}\) & \(\frac{kg}{s}\) & \(\frac{kg}{s}\) & \(m\) & \(\frac{kgm}{s}\) \\ \hline \end{tabular}
\end{table} TABLE II: Parameter values of the rotating–translating mass displayed in Fig. 9.

Fig. 6: References \(r_{1}(k)\) and \(r_{2}(k)\) (top window), and the generated feedforward input for \(r_{1}(k)\) using PGNNs with \(\lambda=0\) (middle window) and \(\lambda=0.07\) (bottom window).

Fig. 7: Tracking error for reference \(r_{1}(k)\) resulting from the different feedforward controllers.
version, i.e., (16), (17) with linear physical model \(f_{\text{phy}}\big{(}b_{\text{phy}},\phi(k)\big{)}=\theta_{\text{phy}}^{T}\phi(k)\) and parameters identified according to (14) with MSE cost function (15);
3. _PGNN-based feedforward_ with NN dimensions \(n_{1}=16\) and \(l=1\), using ZPETC stable inversion, i.e., (58) identified according to Lemma 5.2 which guarantees that the feedforward control input remains bounded;
4. _Physics-based feedforward_ with extended preview, i.e., (16), (17) with linear physical model \(f_{\text{phy}}\big{(}\theta_{\text{phy}},\tilde{\phi}(k)\big{)}=\theta_{\text{ phy}}^{T}\tilde{\phi}(k)\), \(n_{\text{pw}}=20\), and parameters identified according to (14) with MSE cost function (15);
5. _PGNN-based feedforward_ with NN dimensions \(n_{1}=16\) and \(l=1\) and extended preview window, i.e., (56) with \(n_{\text{pw}}=20\) identified according to Lemma 5.1.

We choose \(\tilde{Q}=Q=I\) to find \(\Theta\) in (54) for the extended preview PGNN (56) or in (60) for the stably inverted PGNN (58). Consider Lemma 5.2 as an example. We find \(P\) by solving \(A(\tilde{b}^{d})^{T}PA(\tilde{b}^{d})-P+\tilde{Q}=0\), which gives \(\beta=0.4982\) from (53). Consequently, we use \(\tilde{b}^{n}\) in (59) obtained from ZPETC approximation and combine it with (60), to obtain \(\|L_{x}\|_{2}^{2}\leq\|403L_{\text{\tiny{up}}}\|_{2}^{2}<1.96\cdot 10^{-4}\), which directly translates into NN parameter constraints when using (49).

_Results:_ Fig. 10 shows the tracking error resulting from the aforementioned feedforward controllers for the reference shown in the top window. It is clear that the PGNN manages to significantly outperform the physics-based feedforward controller in terms of tracking error, as is confirmed by the mean-squared error (MSE) values listed in Table III.

## VII Conclusions

In this paper, we presented a generalized framework for nonlinear feedforward control using physics-guided neural networks. It was demonstrated that, under suitable assumptions, the PGNNs can recover the inverse system dynamics exactly. Moreover, it was shown that the use of physical knowledge within the PGNN framework enabled tools to ensure improved training convergence with respect to the physical model, as well as providing a form of graceful degradation. A systematic method for analyzing input-to-state stability of PGNN-based feedforward controllers was provided, based on which the identification criterion was revised to impose ISS before training. The PGNN feedforward showed superior performance compared to physics-based feedforward controllers on a real-life industrial linear motor and a nonminimum phase simulation example. This shows that indeed, the developed PGNN-based feedforward controller design methodology offers a systematic and reliable solution for increasing performance of mechatronic systems with unknown, complex nonlinearities and parasitic effects.

Fig. 8: MAE of the tracking errors for \(r_{1}(k)\) and \(r_{2}(k)\) displayed in Fig. 7 with different velocities \(\max\big{(}\dot{r}(k)\big{)}\).

Fig. 10: Tracking error resulting from different feedforward controllers.

\begin{table}
\begin{tabular}{|l|l|l|l|l|l|} \cline{2-7} \multicolumn{1}{c|}{**No FF**} & \multicolumn{1}{c|}{**Un.-2Perc**} & \multicolumn{1}{c|}{**pons-2Perc**} & \multicolumn{1}{c|}{**Un.-PW**} & \multicolumn{1}{c|}{**pons-PW**} \\ \hline \hline
**MSE** & \(4.50\cdot 10^{-5}\) & \(1.44\cdot 10^{-7}\) & \(2.10\cdot 10^{-8}\) & \(1.59\cdot 10^{-7}\) & \(4.74\cdot 10^{-4}\) \\ \hline \end{tabular}
\end{table} TABLE III: MSE in \([m^{2}]\) of the tracking errors in Fig. 10.

Fig. 9: Rotating–translating mass with actuation and sensing on opposite sides of the centre of mass.
## Appendix I Proof of Proposition 4.1

Substitution of the inverse dynamics (20) and PGNN (24) into the cost function (25) with \(\Lambda_{\text{NN}}=0\) gives

\[\frac{1}{N} \sum_{i=0}^{N-1}\Big{(}f_{\text{phy}}(\theta_{\text{phy}}^{*},\phi _{i})+g(\phi_{i})-f_{\text{phy}}(\theta_{\text{phy}},\phi_{i}) \tag{66}\] \[-f_{\text{NN}}(\theta_{\text{NN}},\phi_{i})\Big{)}^{2}+\lambda \left\|\begin{bmatrix}0&0\\ 0&\Lambda_{\text{phy}}\end{bmatrix}\left(\theta-\begin{bmatrix}0\\ \theta_{\text{phy}}^{*}\end{bmatrix}\right)\right\|_{2}^{2}.\]

Both terms in (66) are non-negative, such that the global minimum is attained for \(\hat{\theta}_{\text{phy}}=\theta_{\text{phy}}^{*}\) (regularization term) and \(\hat{\theta}_{\text{NN}}=\theta_{\text{NN}}^{*}\) (MSE term, after substitution of \(\theta_{\text{phy}}=\theta_{\text{phy}}^{*}\)).

## Appendix II Proof of Proposition 4.2

Proof:: The proof of (34) follows directly by observing that \(\theta_{\text{L}}^{(j)}\) in (32) is the least squares solution of (14), (25) for the PGNN (30) given \(\theta_{\text{NL}}^{(j)}\). Since \(M(\theta_{\text{NL}}^{(j)})\) is non-singular, \(\theta_{\text{L}}^{(j)}\) is unique, such that (34) holds with strict inequality if and only if \(\theta_{\text{L}}^{(j)}\neq\overline{\theta}_{\text{L}}\). Observe that, \(\theta_{\text{L}}^{(j)}-\overline{\theta}_{\text{L}}=\theta_{\text{L}}^{(j)}- M\big{(}\theta_{\text{NL}}^{(j)}\big{)}^{-1}M\big{(}\theta_{\text{NL}}^{(j)} \big{)}\overline{\theta}_{\text{L}}\neq 0\) gives condition (35) after substitution of (32) and (33) and using nonsingularity of \(M(\theta_{\text{NL}}^{(j)})\).

Secondly, (36) follows by observing that choosing \(\theta_{\text{L}}^{(j)}\neq\overline{\theta}_{\text{L}}\) cannot decrease \(V_{\text{reg}}\) when \(\Lambda=0\). Correspondingly, (34) states that \(V_{\text{MSE}}\) must decrease if (35) holds, such that

\[V_{\text{MSE}}\left(\begin{bmatrix}\theta_{\text{NL}}^{(j)}\\ \theta_{\text{L}}^{(j)}\end{bmatrix},Z^{N}\right) \leq V\left(\begin{bmatrix}\theta_{\text{NL}}^{(j)}\\ \overline{\theta}_{\text{L}}\end{bmatrix},Z^{N}\right) \tag{67}\] \[=V_{\text{MSE}}(\theta_{\text{phy}}^{*},Z^{N}),\]

and with strict inequality if (35) holds.

## Appendix III Proof of Theorem 5.1

Proof:: The proof follows by showing that \(V\big{(}x(k)\big{)}=x(k)^{T}Px(k)\) is an ISS-Lyapunov function as in Definition 5.2. Condition (46) is satisfied with \(\kappa_{1}\big{(}x(k)\big{)}=\lambda_{\text{min}}(P)x(k)^{T}x(k)\) and \(\kappa_{2}\big{(}x(k)\big{)}=\lambda_{\text{max}}(P)x(k)^{T}x(k)\). We compute the difference \(V\big{(}x(k+1)\big{)}-V\big{(}x(k)\big{)}\) to obtain

\[V\big{(}x(k+1)\big{)}-V\big{(}x(k)\big{)}=-x(k)^{T}Qx(k) \tag{68}\] \[\qquad\qquad+2C_{1}B^{T}PAx(k)+B^{T}PBC_{1}^{2},\] \[C_{1}:=\hat{a}^{T}\phi_{r}(k+1)+f_{\text{NN}}\left(\hat{\theta}_ {\text{NN}},\begin{bmatrix}\phi_{r}(k+1)\\ x(k)\end{bmatrix}\right).\]

For a term \(2p^{T}q\) and a \(\varepsilon\in\mathbb{R}_{>0}\), we can complete the squares as:

\[2p^{T}q=\varepsilon p^{T}p-\varepsilon(p-\frac{1}{\varepsilon}q)^{T}(p-\frac{ 1}{\varepsilon}q)+\frac{1}{\varepsilon}q^{T}q\leq\varepsilon p^{T}p+\frac{1}{ \varepsilon}q^{T}q. \tag{69}\]

We complete the squares (69) of \(2C_{1}B^{T}PAx(k)\) in (68) using \(\varepsilon=\beta\lambda_{\text{min}}(Q)\), \(\beta\in\mathbb{R}_{>0}\), \(p=x(k)\), and \(q=A^{T}PBC_{1}\) to obtain

\[V\big{(}x(k+1)\big{)}-V\big{(}x(k)\big{)}\leq- (1-\beta)\lambda_{\text{min}}(Q)x(k)^{T}x(k) \tag{70}\] \[+c_{\beta}C_{1}^{2},\]

with \(c_{\beta}:=B^{T}PB\left(I+\frac{1}{\beta\lambda_{\text{min}}(Q)}AA^{T}P\right)B\). Similarly, by substituting \(C_{1}\) in (70) and completing the squares (69) for \(2\hat{a}^{T}\phi_{r}(k)f_{\text{NN}}\) using \(\varepsilon=\beta_{1}\), \(\beta_{1}\in\mathbb{R}_{>0}\), \(p=f_{\text{NN}}\) and \(q=\hat{a}^{T}\phi_{r}(k+1)\), we obtain

\[V\big{(}x(k+1)\big{)}-V\big{(}x(k)\big{)} \leq-(1-\beta)\lambda_{\text{min}}(Q)x(k)^{T}x(k) \tag{71}\] \[+c_{\beta}(1+\beta_{1})f_{\text{NN}}\left(\hat{\theta}_{\text{NN} },\begin{bmatrix}\phi_{r}(k+1)\\ x(k)\end{bmatrix}\right)^{2}\] \[+c_{\beta}(1+\frac{1}{\beta_{1}})\big{(}\hat{a}^{T}\phi_{r}(k+1) \big{)}^{2},\]

with \(\beta_{1}\in\mathbb{R}_{>0}\). Finally, substitution of the Lipschitz condition (48) with \(\phi_{r}^{B}(k)=0\) and (50), and completing the squares (69) for \(2\hat{L}_{uq}^{T}x(k)L_{r}^{T}\phi_{r}(k+1)\) using \(\varepsilon=\beta_{2}\), \(\beta_{2}\in\mathbb{R}_{>0}\), \(p=L_{uq}^{T}x(k)\), and \(q=L_{r}^{T}\phi_{r}(k+1)\), gives

\[V\big{(}x(k+1)\big{)}-V\big{(}x(k)\big{)} \leq-x(k)^{T}\Big{(}(1-\beta)\lambda_{\text{min}}(Q)-c_{\beta}. \tag{72}\] \[(1+\beta_{1})(1+\beta_{2})L_{uq}L_{u_{\sharp}}^{T}\Big{)}x(k)+ \phi_{r}(k+1)^{T}c_{\beta}\cdot\] \[\Big{(}(1+\beta_{1})(1+\frac{1}{\beta_{2}})L_{r}L_{r}^{T}+(1+\frac {1}{\beta_{1}})\hat{a}\hat{a}^{T}\Big{)}\phi_{r}(k+1)\] \[=:-\kappa_{3}\big{(}|x(k)|\big{)}+\sigma\big{(}|\phi_{r}(k+1)| \big{)}.\]

It is clear that \(\kappa_{3}\big{(}|x(k)|\big{)}\) is a \(\mathcal{K}_{\infty}\)-function if the matrix between \(x(k)^{T}\) and \(x(k)\) is positive definite, which reduces to the scalar condition (52) when using the Cauchy-Schwartz inequality, i.e., \(\big{(}L_{uq}^{T}x(k)\big{)}^{2}=x(k)^{T}L_{uq}L_{uq}^{T}x(k)\leq L_{uq}^{T}L_{uq} x(k)^{T}x(k)\), and choosing \(\beta_{1}\), \(\beta_{2}\) arbitrary small.

## Appendix IV Proof of Lemma 5.1

Proof:: The proof follows directly from Proposition 5.1 and observing that limiting \(\theta\in\Theta\) ensures that (52) is satisfied by using the NN Lipschitz bound in (49).

## Appendix V Proof of Lemma 5.2

Proof:: The proof follows by rewriting the stable inverted PGNN (58) into state-space form, and following the proof of Proposition 5.1 with the Lipschitz bound in (59). A \((P,\tilde{Q})\) pair exists, since the stable inversion of \(G(\tilde{b},q^{-1})\) gives a Schur \(A(\tilde{b})\).
