Physics-Informed Neural Networks for UAV System Estimation

1. Introduction

One of the main challenges in control theory consists in realizing a dynamic model of the physical system under study to understand and predict its behavior over time. This may lead to a highly complex mathematical description of the considered system, also requiring to yield dedicated experiments to estimate the unknown model parameters [1]. This is particularly true for complex systems, such as, for instance, Unmanned Aerial Vehicle (UAVs), since they are considered unstable multiple-input and multiple-output (MIMO) dynamics systems. When actual real data are not available (e.g., due to the lack of reliable knowledge about the system behavior or the lack of or inaccuracy of measurements), an efficient alternative is the so-called system identification, which relies on obtaining the mathematical models of a dynamic system from measured data of its input and output parameters. Two approaches can be adopted for system identification, namely, grey box model, also known as semi-physical model (i.e., when the constructed model still has a number of unknown free parameters that are estimated through system identification) and black box model (i.e., when no prior model is available) [2]. Regardless of the adopted model, system identification strategies consist in designing and conducting an identification experiment to collect data, selecting the structure of the dynamic system also specifying which parameters are to be identified, and fitting the model parameters to the obtained data; the overall quality of the resulting model is then evaluated through a validation procedure [1].

Research on systems identification has garnered considerable attention in recent decades, starting from 50s, and has become a vital discipline in various applications within the field of automatic control. This includes areas such as robotics, industrial processes, reduced-order modeling, and model testing.

In recent decades, as researchers and engineers have gained access to larger amounts of data, one of the most widely used methods for system identification has been the Least Squares Method. This approach can be considered the foundation of data-driven system modeling.

In [3], a classification of models identification is presented based on the applied methodologies. The classification includes Fuzzy Logic Theory, Genetic Algorithm, Neural Network, Swarm Intelligence Optimization Algorithms, Auxiliary Model Identification Method, Hierarchical Method, and Stochastic Theory [4].

Referring to the last method, in [5] crucial principles of realization theory were introduced, particularly emphasizing controllability and observability concepts. Their focus was on linear systems identification, with a special emphasis on establishing the minimal state-space representation to define the subspace within which the system dynamics evolves. They pioneered the concept of state-variable equations, which facilitates the realization of an external description through an equivalent internal description of a dynamical system. A second algorithm addressing the same problem was simultaneously presented by Kalman (see [6]), utilizing the principles of controllability and observability and necessitating linear algebra-type computations. State-space models are particularly well-suited for this approach as they lend themselves to linear algebra techniques, robust numerical simulation, and modern control design methods. A few years later, Ho and Kalman [5] approached the same problem from a new perspective. They demonstrated that the minimum realization problem is equivalent to a representation problem involving a sequence of real matrices known as Markov parameters (pulse response samples).

During the 90s, expanding upon the foundational work by Gilbert and Kalman, various methods were developed to identify the most observable and controllable subspace of a system based on given input-output (I/O) data. Some years later, at NASA, Juang devised an approach to concurrently ascertain a linear state-space model and the corresponding Kalman filter using noisy input-output measurements (see [7] for details). Referred to as the Observer/Kalman Identification Algorithm (OKID) and exclusively formulated in the time domain, this method computes the Markov parameters of a linear system. Subsequently, both the state-space model and a corresponding observer are determined simultaneously from these parameters. The Kalman Filter was therefore a reference point for identification in the case of linear systems and its parameters. Subsequently, the theory was extended to nonlinear systems through the Extended Kalman Filter (EKF) in various versions (see, e.g., the work in [8] for the polynomial one), starting from a series of measurements subject to noise, and was successful in any application that can be linked to control theory.

In this regard, it is worth noting that, an interesting comparison among modern methods for parameter estimation of aircrafts is presented in [9]. In particular, by using real flight data, the authors find that the performance of Filter Error Methods (FEMs), such as Extended Kalman Filter (EKF) and Unscented Kalman Filter (UKF), prevail over Equation Error Methods (EEMs) and Output Error Methods (OEMs), particularly in the presence of turbulence and noise.

More recently, the use of many machine learning tools for control problems is gaining attention. The technology of artificial neural networks has undergone significant development since the late 20th century. Over time, neural networks have found applications in various fields, particularly for cyber-physical systems, with notable contributions in intelligent control, nonlinear optimization, computer vision, biomedical engineering, robotics, and system identification. Neural network identification leverages the structure of a nonlinear system, enabling the simulation of input-output relationships through the network’s ability to approximate any nonlinear mapping. It employs the self-adaptation and self-learning capabilities of neural networks to implement simple learning algorithms in engineering, ultimately achieving the forward or inverse model of the system through training. This method exhibits the following characteristics:

• No need to establish the model structure of the actual system, as the neural network itself serves as a model for network identification.
• It is capable of identifying any linear or nonlinear model.
• The neural network not only serves as a model but is also an actual system achievable through physics.

However, it comes with certain limitations:

• Local minimum problem.
• Lengthy training time and slow learning speed.
• Difficulty in extracting ideal training samples.
• Challenges in optimizing the network structure.
• Difficulty in completely solving the convergence problem theoretically for the neural network algorithm.

To address these limitations, numerous researchers have proposed improved neural network identification methods, resulting in favorable identification outcomes (see, e.g., [10,11]).

In the above described context, the main contribution of this paper is to exploit the power in the management and use of a limited amount of data available data in Neural Networks with the physical laws of the flight dynamics of a quadcopter. In particular, we use Physics-Informed Neural Networks (PINNs), also known as Theory-Trained Neural Networks, to estimate the drone model. In fact, in the last decade, aerial robots, and particularly small UAVs and drones, have witnessed a continuously increasing use in a wide range of services (see, e.g., [12]). As a consequence, the identification and control (and hence the system parameters estimation) of these systems is of paramount importance from different point of view (see [13,14]).

We conclude this section highlighting that the use of PINNs for system identification is recently reaching relevant interest in the research community. As an example, the work in [15] proposes the use of PINNs for non-linear system identification for power system dynamics. In particular, the final aim of this contribution is to discover the frequency dynamics of future power systems. To assess the performance of the presented model and validate their proposal, the authors compare the estimator against the UKF. PINNS are also used in [16] to discover the ordinary differential equation (ODE) of a rotor system from noise measurements and assess the healthy/faulty machine condition. The validation is performed on a test bench. System identification of structural systems with a multiphysics damping model through PINNs is presented in [17], while in [18] they are adopted to estimate motion and identify system parameters of a moored buoy under different sea states. To the best of authors knowledge, a first attempt in the use of PINNs for quadrotor dynamical modeling is proposed in [19], where a comparison with other existing approaches is performed, showing that PINNs outperform linearized mathematical models and classical black-box approaches based on Deep Neural Networks (DNNs) in terms of test error, while also exhibiting better capability of underlying existing relationships between parameters. In particular, in [19] the orientation estimation of a quadrotor is dealt with, while in the present paper we focus on all the state variables, therefore also the spatial ones.

The reminder of the paper is organized as follows. In Section 2 the drone dynamics model is recalled. Section 3 illustrates the two methodologies chosen to compare system model identification: a quick recall of the Extended Kalman Filter and then the idea on which PINNs are based is described. Section 4 shows the simulation and comparison results between the methods mentioned above. Some final considerations and future research ideas conclude the paper.

2. Mathematical Model of a Quadrotor

The quadrotor analyzed in this study is composed of a rigid structure with four rotors (for more details, refer to [20,21]). The propellers produce a force $F_i=b·{\omega }_{p,i}^2$, which is proportional to their angular velocity ${\omega }_{p,i}$, for $i=1,2,3,4$. Propellers 1 and 3 spin counterclockwise, while propellers 2 and 4 spin clockwise. In the following discussion, the quadrotor’s orientation will be described using Euler angles. The frames $RC(O,e_1,e_2,e_3)$ and $R\mathrm{\Gamma }(\mathrm{\Omega },{\epsilon }_1,{\epsilon }_2,{\epsilon }_3)$) represent the reference systems fixed to the Earth and the quadrotor, respectively, with $\mathrm{\Omega }$ located at the quadrotor’s center of mass (see Figure 1).

The quadrotor’s position in $RC$ is represented by $p={(x,y,z)}^T$, while its orientation is defined by the Euler angles $\alpha ={(\varphi ,\theta ,\psi )}^T$, where $\varphi \in [-\pi /2,\pi /2)$, $\theta \in (-\pi /2,\pi /2)$, and $\psi \in [-\pi ,\pi )$ correspond to the roll, pitch, and yaw angles, respectively. The 3–2–1 sequence is utilized here, as in [22]. Furthermore, $v={(v_x,v_y,v_z)}^T$ and $\omega ={({\omega }_1,{\omega }_2,{\omega }_3)}^T$ denote the linear and angular velocities of the quadrotor’s center of mass, expressed in $RC$ and $R\mathrm{\Gamma }$, respectively. The quadrotor’s translational dynamics are expressed in $RC$, while its rotational dynamics are expressed in $R\mathrm{\Gamma }$:

$$\begin{array}{cc}\hfill \dot{p}& =v\hfill \\ \hfill \dot{v}& ={\displaystyle \frac{1}{m}}\left(f_p+f_g\right)\hfill \\ \hfill \dot{\alpha }& =M\left(\alpha \right)\omega \hfill \\ \hfill \dot{\omega }& =J^{-1}\left(-\tilde{\omega }J\omega +{\tau }_p+{\tau }_{gy}\right)\hfill \end{array}$$

(1)

where m represents the mass of the quadrotor, $f_p=\mathcal{R}\left(\alpha \right){\phi }_p$ is the force exerted by the propellers in $RC$ (${\phi }_p$ is the same force, but expressed in $R\mathrm{\Gamma }$), and J (a positive definite symmetric matrix in ${\mathbb{R}}^{3\times 3}$, expressed in $R\mathrm{\Gamma }$) is the quadrotor’s inertia matrix, and

$$\tilde{\omega }=\left(\begin{array}{ccc}0& -{\omega }_3& {\omega }_2\\ {\omega }_3& 0& -{\omega }_1\\ -{\omega }_2& {\omega }_1& 0\end{array}\right)$$

(2)

is the so-called dyadic representation of $\omega$. Furthermore,

$${\phi }_p=\left(\begin{array}{c}0\\ 0\\ u_p\end{array}\right),u_p=\sum _{i=1}^4{F}_i,f_g=\left(\begin{array}{c}0\\ 0\\ -mg\end{array}\right)$$

(3)

$${\tau }_p=\left(\begin{array}{c}{\tau }_1\\ {\tau }_2\\ {\tau }_3\end{array}\right)=\left(\begin{array}{c}\ell b({\omega }_{p,2}^2-{\omega }_{p,4}^2)\\ \ell b({\omega }_{p,3}^2-{\omega }_{p,1}^2)\\ c({\omega }_{p,1}^2-{\omega }_{p,2}^2-{\omega }_{p,3}^2-{\omega }_{p,4}^2)\end{array}\right)$$

(4)

are the input forces (as shown in 3) and moments (as shown in Section 2) generated by the propellers (expressed in $R\mathrm{\Gamma }$), where ℓ is the distance from the center of mass $CG$ to the rotor shaft, and b (with units [b]=${\displaystyle \frac{N,s^2}{ra{d}^2}}$) and c (with units [c]=${\displaystyle \frac{N,s^2,m}{ra{d}^2}}$) are the thrust and drag coefficients, respectively. It is evident that:

$$\begin{array}{c}\left(\begin{array}{c}{F}_1\\ F_2\\ F_3\\ F_4\end{array}\right)={\left(\begin{array}{cccc}1& 1& 1& 1\\ 0& \ell & 0& -\ell \\ -\ell & 0& \ell & 0\\ c& -c& c& -c\end{array}\right)}^{-1}\left(\begin{array}{c}{u}_p\\ {\tau }_1\\ {\tau }_2\\ {\tau }_3\end{array}\right),\\ \left(\begin{array}{c}{\omega }_{p,1}^2\\ {\omega }_{p,2}^2\\ {\omega }_{p,3}^2\\ {\omega }_{p,4}^2\end{array}\right)=\left(\begin{array}{cccc}\frac{1}{4b}& 0& \frac{-1}{2b\ell }& \frac{1}{4c}\\ \frac{1}{4b}& \frac{1}{2b\ell }& 0& \frac{-1}{4c}\\ \frac{1}{4b}& 0& \frac{1}{2b\ell }& \frac{1}{4c}\\ \frac{1}{4b}& \frac{-1}{2b\ell }& 0& \frac{-1}{4c}\end{array}\right)\left(\begin{array}{c}{u}_p\\ {\tau }_1\\ {\tau }_2\\ {\tau }_3\end{array}\right),\\ {\omega }_{p,i}=\sqrt{{\displaystyle \frac{F_i}{b}}},0\le F_i\le F_{i,max},0\le {\omega }_{p,i}\le {\omega }_{p,i,max}\end{array}$$

(5)

with $i=1,2,3,4$, where $F_{i,max}$ and ${\omega }_{p,i,max}$ represent the maximum forces and angular velocities for each propeller, constrained by physical limitations. Additionally, $f_g$ in (3) refers to the gravitational force, expressed in $RC$. Vectors expressed in $R\mathrm{\Gamma }$ are converted into vectors in $RC$ using the rotation matrix

$$\mathcal{R}\left(\alpha \right)=\left(\begin{array}{ccc}{c}_{\theta }{c}_{\psi }& s_{\varphi }{s}_{\theta }{c}_{\psi }-c_{\varphi }{s}_{\psi }& c_{\varphi }{s}_{\theta }{c}_{\psi }+s_{\varphi }{s}_{\psi }\\ c_{\theta }{s}_{\psi }& s_{\varphi }{s}_{\theta }{s}_{\psi }+c_{\varphi }{c}_{\psi }& c_{\varphi }{s}_{\theta }{s}_{\psi }-s_{\varphi }{c}_{\psi }\\ -s_{\theta }& s_{\varphi }{c}_{\theta }& c_{\varphi }{c}_{\theta }\end{array}\right)$$

(6)

where $c_{★}=cos(★)$, $s_{★}=sin(★)$, $★=\varphi ,\theta ,\psi$. The angular velocity dynamics are expressed using the following matrix:

$$M\left(\alpha \right)=\left(\begin{array}{ccc}1& s_{\varphi }{\mathrm{tg}}_{\theta }& c_{\varphi }{\mathrm{tg}}_{\theta }\\ 0& c_{\varphi }& -s_{\varphi }\\ 0& s_{\varphi }{\mathrm{sc}}_{\theta }& c_{\varphi }{\mathrm{sc}}_{\theta }\end{array}\right)$$

where ${\mathrm{tg}}_{★}=tan(★)$ and ${\mathrm{sc}}_{★}=sec(★)$ with $★=\varphi ,\theta ,\psi$.Assuming small angles for $\varphi$ (roll) and $\theta$ (pitch), which is reasonable for a quadrotor performing non-aggressive maneuvers, this matrix can be approximated by the identity matrix, i.e.,$M\left(\alpha \right)\simeq I_{3x3}$[23]. Therolling torque ${\tau }_1$ is generated by the forces $F_2$ and $F_4$, while the pitching torque ${\tau }_2$ is generated by the forces $F_1$ and $F_3$. According to Newton’s third law, the propellers exert a yawing torque ${\tau }_3$ on the quadrotor body in the direction opposite to the propeller rotation. Furthermore, the gyroscopic torque arising from the propeller rotations is given by

$${\tau }_{gy}=\sum _{i=1}^4{(-1)}^i{J}_{p,i}{\omega }_{p,i}\tilde{\omega }{ϵ}_3$$

(7)

where $J_{p,i}$, for $i=1,2,3,4$, represents the moment of inertia of the $i_{\mathrm{th}}$ motor and propeller about its axis of rotation. Finally, the gyroscopic torque can also be expressed as

$${\tau }_{gy}^{\circ }=\sum _{i=1}^4{(-1)}^i{J}_p^{\circ }{\omega }_{p,i}\tilde{\omega }{ϵ}_3=J_p^{\circ }{\omega }_p\left(\begin{array}{c}{\omega }_2\\ -{\omega }_1\\ 0\end{array}\right)$$

(8)

where ${\omega }_p=-{\omega }_{p,1}+{\omega }_{p,2}-{\omega }_{p,3}+{\omega }_{p,4}$ is referred to as the rotor relative speed. Given these conditions, the mathematical model (1) of the quadrotor can be rewritten as:

$$\begin{array}{cc}\hfill \dot{x}\left(t\right)& =v_x\left(t\right),\dot{y}\left(t\right)=v_y\left(t\right),\dot{z}\left(t\right)=v_z\left(t\right)\hfill \\ \hfill {\dot{v}}_x\left(t\right)& =\left(c_{\varphi \left(t\right)}{s}_{\theta \left(t\right)}{c}_{\psi \left(t\right)}+s_{\varphi \left(t\right)}{s}_{\psi \left(t\right)}\right){\displaystyle \frac{u_p\left(t\right)}{m}}\hfill \\ \hfill {\dot{v}}_y\left(t\right)& =\left(c_{\varphi \left(t\right)}{s}_{\theta \left(t\right)}{s}_{\psi \left(t\right)}-s_{\varphi \left(t\right)}{c}_{\psi \left(t\right)}\right){\displaystyle \frac{u_p\left(t\right)}{m}}\hfill \\ \hfill {\dot{v}}_z\left(t\right)& =-g+c_{\varphi \left(t\right)}{c}_{\theta \left(t\right)}{\displaystyle \frac{u_p\left(t\right)}{m}}\hfill \\ \hfill \dot{\varphi }\left(t\right)& ={\omega }_1\left(t\right),\dot{\theta }\left(t\right)={\omega }_2\left(t\right),\dot{\psi }\left(t\right)={\omega }_3\left(t\right)\hfill \\ \hfill {\dot{\omega }}_1\left(t\right)& ={\displaystyle \frac{J_2-J_3}{J_1}}{\omega }_2\left(t\right){\omega }_3\left(t\right)+{\displaystyle \frac{J_p}{J_1}}{\omega }_p\left(t\right){\omega }_2\left(t\right)+{\displaystyle \frac{1}{J_1}}{\tau }_1\left(t\right)\hfill \\ \hfill {\dot{\omega }}_2\left(t\right)& ={\displaystyle \frac{J_3-J_1}{J_2}}{\omega }_1\left(t\right){\omega }_3\left(t\right)-{\displaystyle \frac{J_p}{J_2}}{\omega }_p\left(t\right){\omega }_1\left(t\right)+{\displaystyle \frac{1}{J_2}}{\tau }_2\left(t\right)\hfill \\ \hfill {\dot{\omega }}_3\left(t\right)& ={\displaystyle \frac{J_1-J_2}{J_3}}{\omega }_1\left(t\right){\omega }_2\left(t\right)+{\displaystyle \frac{1}{J_3}}{\tau }_3\left(t\right).\hfill \end{array}$$

(9)

where the state space vector is $\left[xyz{v}_x{v}_y{v}_z\varphi \theta \psi {\omega }_1\omega 2{\omega }_3\right]$, the control input are $\left[u_p{\tau }_1{\tau }_2{\tau }_3\right]$. Given the aims of the research, it is assumed that we have a limited amount of data available in the form input-output (state). Moreover, the system identification methods that will be presented in Section 3 are in discrete form, Thus, the quadrotor system equations (9) will be implemented in discrete form by keeping the functions constant over each time interval $[t_i,t_{i+1})$, where $t_{i+1}-t_i=\mathsf{\Delta }t$ for $i=0,\dots ,i_{max}-1$, and $t_{i_{max}}=t_f$ represents the final mission time. This discrete-time formulation allows for accurate modeling within each interval, enhancing the overall performance of the UAV system. The digitalization is obtained by a first-order discretization.

3. System Identification Methods

3.1. Extended Kalman Filter

The EKF ([24]) is an adaptation of the standard Kalman filter designed for use when the system and/or measurement models are nonlinear. The approach underlying the extended Kalman filter is based on the following procedure. Given the system type:

$$\begin{array}{c}\hfill s_k=f(s_{k-1},{\mu }_{k-1})+n_{1,k}\\ \hfill z_k=h\left(x_k\right)+n_{2,k}\end{array}$$

(10)

where $s_k$ is the state space vector, $z_k$ is the output and ${\mu }_{k-1}$ is the control input, it approximates the nonlinear functions of the system $x_k$ and $z_k$, through a Taylor series expansion stopped at the first order around the current estimate, thus making the system linear. The noise vectors $n_{1,k}$ and $n_{2,k}$ represent disturbances affecting the state and measurements, respectively, and are assumed to be uncorrelated with each other. Given a random variable s we want to know the probability density of the variable y obtained from the transformation of the variable s. Moments up to first order are studied.

$$\begin{array}{c}\hfill s\sim N(\widehat{s},{\sigma }_s^2)=\widehat{s}+{\delta }_s\mathrm{with}{\delta }_s\sim N(0,{\sigma }_s^2)\end{array}$$

(11)

$$\begin{array}{c}\hfill y=f\left(s\right)=f(\widehat{s}+{\delta }_s)\end{array}$$

(12)

whose development is $y=f\left(\widehat{s}\right)+\lambda f{\delta }_s+{\displaystyle \frac{1}{2}}{\mathsf{\Delta }}^{2}f{\delta }_s^2$. Then, $J_s^h={\left[{\displaystyle \frac{\partial h}{\partial s}}\right]}_{s={\widehat{s}}_{k|k}}$ is defined as the Jacobian matrix of the function h with respect to s and $J_s^f={\left[{\displaystyle \frac{\partial f}{\partial s}}\right]}_{s={\widehat{s}}_{k|k}}$ is the Jacobian matrix of f respect to s. Moreover, the conditional probability function is assumed to be Gaussian. Using these approximations the system turns out to be linear, so the Kalman filter can be applied, obtaining the following recursive equations (respectively referred to the prediction and the correction):

$$\begin{array}{c}\hfill {\widehat{s}}_{k|k-1}=f(s_{k|k-1},{\mu }_k)\\ \hfill P_{k|k-1}=J_s^f{P}_{k-1|k-1}{J}_s^{f^T}+Q\end{array}$$

(13)

$$\begin{array}{c}\hfill n_{2,k}=z_k-h\left({\widehat{s}}_{k|k-1}\right)\\ \hfill S_k=J_s^h{P}_{k|k-1}{J}_s^{h^T}+R_k\\ \hfill K_k=P_{k|k-1}{J}_s^{h^T}{S}_k^{-1}\\ \hfill s_{k|k}=s_{k|k-1}+K_k{n}_{2,k}\\ \hfill P_{k|k}=(I-K_k{J}_s^h)P_{k|k-1}\end{array}$$

(14)

where the matrices Q and R have the same properties as in the linear Kalman filter. However, there are no assurances regarding the quality of the estimates produced, and the extended Kalman filter is highly sensitive to the accuracy of the initial estimates (see [24]).

3.2. Physics-Informed Neural Networks

PINNs are introduced to integrate physical laws, typically described by Ordinary Differential Equations (ODEs), into Deep Neural Networks (DNNs) (see [25]). This approach trains DNNs in a supervised fashion to comply with given physical laws, enabling the automatic discovery of data-driven solutions for ODEs, as demonstrated in the application under consideration (Equation (9)). The core idea behind PINNs is to incorporate the differential equation into the loss function, as illustrated in Figure 2, enhancing the network’s robustness and facilitating accurate approximations even in data-scarce scenarios. In this study, we focus on parametrized and nonlinear partial differential equations in the general form:

$${\partial }_t\widehat{y}+N\left[\widehat{y};\lambda \right]=0.t\in \left[0,T\right]$$

(15)

where $\widehat{y}(t,x)$ represents the latent (hidden) solution or the state of the dynamic system, and $N\left[·\right]$ denotes a nonlinear differential operator parametrized by $\lambda$. We define $l(t,x)$ as the expression on the left-hand-side of Equation (15), i.e.:

$$l:=\partial \widehat{y}+N\left[\widehat{y}\right]$$

(16)

and then, we continue by modeling $\widehat{y}(t,x)$ using a deep neural network. In this context, $\widehat{y}$ serves as the output of a layered architecture neural network, denoted by $l_w\left(t\right)$, where $\widehat{y}=l_w\left(t\right)$, let $l_w$ denotes the mapping function learned by a deep network with adaptive weights w. In this context, the neural network is expected to learn the solution of a specified ODE as a function of continuous time t.

By using automatic differentiation and the chain rule, we can derive a neural network representing $\widehat{y}(t,x)$. Importantly, this network shares the same parameters as the one representing $\widehat{y}(t,x)$, but may employ different activation functions due to the influence of the differential operator N. Assuming an autonomous system, we train a neural network $\widehat{y}\left(t\right)$ by optimizing its shared parameters with those of $\widehat{y}(t,x)$ and $l(t,x)$. Our goal is to minimize a Mean Squared Error (MSE) cost function.

$$MSE=MS{E}_{\widehat{y}}+\gamma MS{E}_l.$$

(17)

where:

$$MS{E}_{\widehat{y}}={\displaystyle \frac{1}{N_{\widehat{y}}}}\sum _{i=1}^{N_{\widehat{y}}}{\displaystyle \frac{1}{N_t}}\sum _{j=1}^{N_t}{|{\widehat{y}}_i\left(t^j\right)-{\widehat{y}}_{i,ref}^j|}^2$$

(18)

$$MS{E}_l={\displaystyle \frac{1}{N_{\widehat{y}}}}\sum _{i=1}^{N_{\widehat{y}}}{\displaystyle \frac{1}{N_l}}\sum _{j=1}^{N_l}{\left|l\left({\widehat{y}}_i\left(t^k\right)\right)\right|}^2$$

(19)

being $0\le \gamma \le 1$ a hyper-parameter that should reflect how confident we are in the physical constraints of our system, $N_t$ represent the total number of training data samples, $N_l$ the number of collocation points, and $N_{\widehat{y}}$ the number of outputs produced by the neural network. For each output i, we denote the network’s prediction as ${\widehat{y}}_i\left(·\right)$. Given a data pair $(t^j,{\widehat{y}}_{i,ref}^j)$, where j indexes the pair and ${\widehat{y}}_{i,ref}^j$ is the desired output, we can compare it to the network’s prediction ${\widehat{y}}_i\left(·\right)$. In particular, the outputs of the system equations (9) in the application considered are enclosed in the following vector $\widehat{y}=\left[xyz\varphi \theta \psi \right]$.

The initial loss term $MS{E}_{\widehat{y}}$ is associated with the conventional regression cost function applied to the acquired training data ${(t^j,{\widehat{y}}_{i,\mathrm{ref}}^j)}_{j=1}^{N_t}$, typically used to establish the boundary (initial or terminal) conditions of ODEs during their solution. The second loss term, $MS{E}_l$, penalizes deviations in the behavior of $\widehat{y}\left(t\right)$ as measured by $l\left(\widehat{y}\right)$. This ensures that the solution adheres to the required physical properties, as defined by $l\left(\widehat{y}\right)$, at a specific set of randomly chosen collocation points ${\left\{t^k\right\}}_{k=1}^{N_l}$. The experimental findings indicate a substantial reduction in the needed training data size $N_t$ to learn specific dynamical behaviors. This reduction is attributed to the a priori information incorporated from $MS{E}_l$. Considering the assumed representation of the differential equation of the physical system as $l\left(\widehat{y}\right)=0$, the term $MS{E}_l$ serves as an indicator of how effectively the PINN conforms to the solution of the physical model. The physics-informed method outlined in this study offers a unified framework that combines a pre-existing theoretical model, potentially approximate, with measured data from processes. This framework is intended to address shortcomings in the theoretical model or increase the effectiveness of sample data in process modeling.

4. Simulation Results

In this part, we execute some computer simulations to validate the proposed approach and theoretical findings.

4.1. PINNs Hyperparameters Tuning

Hyperparameter tuning involves selecting the best values for a neural network’s hyperparameters, which are parameters set before training that significantly impact the model’s performance. This process is crucial for enhancing the accuracy and efficiency of the model, helping to achieve optimal results. For instance, the learning rate controls the speed at which the model learns. If this value is too high or too low, the model may not fit the data effectively. Therefore, finding the right combination of hyperparameters is key to ensuring good model performance.

There are various methods for hyperparameter tuning such as grid search, random search, and Bayesian optimization. Grid search involves defining a range of hyperparameter values and systematically evaluating the model for every possible combination. In contrast, random search selects random combinations of hyperparameter values to evaluate the model, which can be more efficient since it doesn’t require testing all combinations. Bayesian optimization, a more sophisticated approach, uses a probabilistic model that relates hyperparameters to performance metrics, helping to predict which hyperparameter values are likely to improve the model’s performance. A broad review of tuning algorithms is available in [26], which also explores a genetic approach. The fundamental hyperparameters to adjust include the number of layers, the number of neurons per layer, the learning rate, and the batch size for the neural network. For Physics-Informed Neural Networks (PINNs), there is an additional hyperparameter called lambda, which controls the balance between the data and the model. Future research will likely focus on a deeper theoretical investigation of this specific parameter.

In this study, the random search method was employed for hyperparameter tuning.

4.2. Model and Performance Comparison

To assess the effectiveness of the proposed Physics-Informed Neural Networks method and to make a comparison with the Extended Kalman Filter-based model identification approach outlined in Section 3, a simulation of the quadrotor system described in Section 2 was conducted using Simulink from MATLAB^®. We want to underline that we have available measured data in limited quantities that link the input and output (corresponding to the state). This makes the problem difficult and a hybrid approach between the use of neural networks and impositions of physical laws represents an excellent compromise between results obtained and computational load. In Table 1 the values of the known and unknown model parameters are reported.

An important factor to take into account and which is a strong point of PINNs is that for them every system parameter is assumed to be unknown, while for the EKF it is assumed that some trivial parameters, such as the mass and the distance of the motors from the center of mass, are known. Furthermore, at the simulation level, noise is added to all measured available data in order to make the simulation more realistic. The considered network is comprised of four layers, with each hidden layer featuring 80 neurons activated by the hyperbolic tangent function.

To evaluate the performance of two proposed system estimators in some missions, various metrics including Mean Absolute Error (MAE), Mean Square Error (MSE), Integral Squared Error (ISE), Integral Absolute Error (IAE), and Integral Time–weighted Absolute Error (ITAE) have been utilized as comparative indicators for system identification. The performance indices are defined for the output state variables $\widehat{y}$ (the same reasoning applies to other variables) concerning the reference $y_{ref}$ that needs to be tracked across N sampling instances. They are defined as follows:

$$MAE={\displaystyle \frac{{\sum }_{i=1}^N[y_{ref,i}-y_i]}{N}},MSE={\displaystyle \frac{1}{N}}\sum _{i=1}^N{(y_{ref,i}-y_i)}^2$$

$$ISE=\int {(y_{ref}-y)}^{2}dt,IAE=\int \mid y_{ref}-y\mid dt$$

$$ITAE=\int (t\mid y_{ref}-y\mid )dt$$

where $y_{ref,i}$ and $y_i$ respectively are the values of the variables $y_{ref}$ and y at sampling time i, and ISE, IAE, and ITAE integrate over time. Furthermore, the gains in the EKF and the architecture, some parameters of the neural network as well as a weight on the two cost functionals considered to calculate MSE in the PINNs assume an important value for evaluating the results. Some preliminary findings are reported in Table 2, Table 3, Table 4, Table 5, Table 6 and Table 7, illustrating the disparity between the EKF and PINNs estimators in each performance metric of the monitored variables for the x position (Table 2), y position (Table 3), z position (Table 4), $\varphi$ roll angle (Table 5), $\theta$ pitch angle (Table 6), and $\psi$ yaw angle (Table 7).

Figure 3, Figure 4, Figure 5, Figure 6 and Figure 7 show the trends of the Extended Kalman Filter and PINNs to approximate the real system respectively for position x, yz, and the angles $\varphi$, $\theta$, and $\psi$. Both methods will achieve excellent results considering the type of maneuver presented here, the disturbances present. But it is evident that the PINNs outperform the EKF, with the additional strength of requiring a reduced computational effort, despite they require training prior to their use or during the simulation.

5. Conclusions

In this research work we investigated the praticatability of employing a physics-informed machine learning technique, specifically physics-informed neural networks (PINNs), for nonlinear system model identification in the context of a quadrotor when input-output available real data are limited. Our discussion centered around the potential of PINNs to replace intricate nonlinear dynamics with a more computationally efficient approximation. Leveraging automatic differentiation, PINN approximations enable a cost-effective and straightforward computation of derivatives of the state. The technique is compared on experimental data with the extended Kalman filter, which according to previous literature is the method that has always obtained the best performance. The simulation comparison is performed on both spatial and orientation variables and shows excellent performances from both techniques, but PINNs are better both at the error level with respect to real data and at the computational level, but with a greater preliminary workload, at a higher training level.

Future research extensions will involve analyzing the results concerning a comparative computational analysis and variations in some identification parameters like $\gamma$ to understand the performance moving from an hybrid weighted approach data driven-model based.