Model Predictive Control Technique for Ducted Fan Aerial Vehicles Using Physics-Informed Machine Learning

Abstract

This paper proposes a model predictive control (MPC) approach for ducted fan aerial robots using physics-informed machine learning (ML), where the task is to fully exploit the capabilities of the predictive control design with an accurate dynamic model by means of a hybrid modeling technique. For this purpose, an indigenously developed ducted fan miniature aerial vehicle with adequate flying capabilities is used. The physics-informed dynamical model is derived offline by considering the forces and moments acting on the platform. On the basis of the physics-informed model, a data-driven ML approach called adaptive sparse identification of nonlinear dynamics is utilized for model identification, estimation, and correction online. Thereafter, an MPC-based optimization problem is computed by updating the physics-informed states with the physics-informed ML model at each step, yielding an effective control performance. Closed-loop stability and recursive feasibility are ensured under sufficient conditions. Finally, a simulation study is conducted to concisely corroborate the efficacy of the presented framework.

1. Introduction

1.1. Literature Review

The ducted fan aerial vehicles (DFAVs) are a type of vertical take-off and landing platforms that have the innate capabilities of a helicopter and fixed-wing aircraft [1,2]. These airborne vehicles are used in a wide number of applications in different areas such as transportation, inspection, aerial surveillance, and manipulation [3,4]. Due to their annular fuselage mechanism called duct [5], these aerial robots are known for their safety in congested, cluttered, and hazardous environments [6,7]. Having vector thrust capabilities [8], DFAVs also have several advantages over other configurations, such as fixed-wing platforms, helicopters, open-rotor, and shrouded quadrotors [9]. However, complex flow distribution makes it arduous to model these aircraft accurately [10], consequently making the flight control design more challenging.

Over the years, several flight control systems based on different strategies have been presented to control these vehicles [9]. In this regard, model predictive control (MPC) has proved itself as an advanced control approach that can address many issues associated with aerial applications [11]. One of the simplest techniques is linear MPC, which has been employed for the purpose of attitude control for these aerial platforms [12,13]. An improved unified control framework based on linear MPC for position and attitude control with Kalman filter as a disturbance observer (DOB) has also been studied and validated [14]. Moreover, robust MPC with a DOB and adaptive MPC with two DOBs have been presented to compensate for the effect of model inaccuracies and disturbances [15]. For achieving similar tasks with some improvements compared to robust MPC [15], composite disturbance rejection approaches based on MPC, referred to as MPC-based compound flight control (MPC-CFC), have been proposed [16,17].

In the aforementioned works, some problems can be solved to improve the control performance of future MPC-driven control systems.

• The performance of MPC depends on the accurate model of the system dynamics [18,19], which is challenging to obtain. This problem is even more prevalent in flight control, where varying flight conditions, disturbances, faults, and model uncertainties may lead to ineffective control performance.
• In the disturbance rejection MPC framework based on DOB (e.g., [16,17]), the dynamics of the disturbance profile cannot be faster than the dynamics of DOB, which implies that this type of composite method is not suitable to estimate and compensate for fast time-varying disturbances.

There are two types of dominant techniques for dynamical system modeling [20]: (1) Physics-informed dynamics, where the equations of motion are derived from governing physical rules, and (2) machine learning (ML) based data-driven system modeling [21]. To construct an end-to-end dynamical model, these techniques are employed separately. There is a growing potential for integrating data-driven techniques based on ML in the aerospace industry [22]. One of the commonly used ML techniques is the neural network (NN) [23,24,25,26]. Nevertheless, NN-based methods need massive training data, which may be uninterpretable and difficult to include constraints. These reasons limit their utilization for online system identification [27]. One alternative is to use sparse identification of nonlinear dynamics (SINDy) developed in [28], which has shown effective control and predictive capabilities with MPC in the presence of noise and is more suitable for low and medium data than the NN technique. Moreover, it has demonstrated strong parametric robustness, effective training, and execution time [27]. While the SINDy algorithm is an effective strategy for multiple learning, an adaptive SINDy technique is more suitable that can efficiently update and correct the dynamics online in a repetitive fashion [29,30]. However, it is challenging to seamlessly capture the correct model with ML techniques alone to satisfy the conservation laws and constraints. On the other hand, physics-informed dynamical models can capture constraints and conservation laws. However, they are either too simple or ultra-complex that may become too expensive to handle online due to high computational costs, particularly for varying system dynamics for MPC. In this aspect, a physics-informed ML hybrid-modeling technique can be employed to seamlessly incorporate physical models and data in all circumstances involving high-dimensional, physically understood, and uncertain contexts [31,32,33,34,35,36].

1.2. Contributions

Motivated by the above discussion, this article presents an MPC-based control strategy for DFAVs via a hybrid modeling approach called physics-informed ML. First, the physics-informed model (nominal system) for the aerial robot is constructed by considering the forces and moments offline. Moreover, it is considered that the nominal system is not affected by any disturbances or uncertainties. Thereafter, a data-driven approach called adaptive SINDy is employed for estimating and compensating for any parametric changes online during the flight. The optimization problem is solved online, where in each step, the physics-informed system is updated by the physics-informed ML. This way, the control action applied to the aerial vehicle is optimal with respect to the current state. In 

Table 1. A concise comparison between existing MPC-based control methods and the proposed technique for DFAVs (order by year).

References	Features 1	Advantages	Disadvantages
Linear MPC [12,13]	AMS	Simple design, effective	Lack adequate tracking and robustness.
		computations.
1. Disturbance rejection RMPC 2,	ACDMPS	Considered time-delays, discrete time optimization problem	1. RMPC 2: lack effective control performance,
2. Disturbance rejection adaptive			2. Adaptive MPC: design intricacy may not be effective, feasibility analysis is not available for both schemes.
Observer-based MPC [14]	ACDMPS	Effective computations, easy real-time implementation, considered time delays.	Ineffective control performance, recursive easibility cannot be established for the entire flight.
Compound RMPC 2 [16,17]	ACDMPS	-	Inadequate performance, suitable if the DOB’s dynamics is faster than disturbance dynamics.

¹ A: attitude tracking, C: composite technique, D: DOB, M: physics-informed modeling, P: position tracking, S: simulation study. ² RMPC: Robust MPC.

, the advantages and disadvantages of existing techniques over the proposed scheme are provided. The contributions of this manuscript are expressed in the following manner:

• To fully exploit the potential of MPC, an online-based predictive control algorithm using physics-informed ML without the utilization of DOB is presented, unlike [15,16,17]. By employing the methodology in such a way, the inherent issues encountered in the disturbance rejection MPC framework [16,17] can be solved.
• For efficient response and to avoid any usual computational complexity problem, only the data-driven part in the hybrid modeling is determined online for model correction. To further enhance the computational efficiency of the developed control algorithm, the physics-informed model is also updated by the physics-informed ML model in each step while solving the optimization problem.
• Unlike [12,13,14,15], theoretical properties such as recursive feasibility and closed-loop stability with constraints satisfaction under sufficient conditions are derived.
• In contrast to the existing disturbance rejection framework [16], the designed approach demonstrates effective performance. Furthermore, the constructed approach can be implemented in other robotics systems to attain similar goals.

1.3. Organization

In Section 2, the problem formulation is constructed in the physics-informed modeling and physics-informed ML model subsections. Furthermore, control objectives are defined, and a few assumptions with preliminaries are established for smooth control operation. The control approach includes MPC development and its feasibility and stability proofs and is provided in Section 3. Section 4 illustrates the numerical implementation that includes comparative analysis with some discussion. Section 5 concludes the article.

1.4. Notation

$\mathbb{N}$ and $\mathbb{R}$ denote all non-negative integers and real space, respectively. $\parallel (.)\parallel \triangleq \sqrt{{(.)}^{\mathrm{T}}(.)}$ is the Euclidean norm. $\mathrm{diag}\{{(.)}_{11},{(.)}_{22},...,{(.)}_n\}$ denotes the diagonal matrix with entries ${(.)}_{11},\cdots ,{(.)}_n\in \mathbb{R}$. Q-weighted norm is referred to as ${\parallel (.)\parallel }_Q\triangleq \sqrt{{(.)}^{\mathrm{T}}\mathit{Q}(.)}$, where $\mathit{Q}$ denotes a positive definite matrix. Consider any given system state/input as ∘, $\circ \left(\tau |t_l\right)$ means ∘ at $\tau$ predicted at $t_l$. Superscript ${\circ }^{*}$ is utilized to define state/input after optimization. Accents $\tilde{\circ }$ and $\widehat{\circ }$ represent the nominal state/nominal control input and estimated state/input, respectively. In the entire manuscript, physics-informed dynamics and physics-informed ML are referred to as nominal and real systems, and these terms are used interchangeably. A bold font style is employed to represent vectors and matrices.

2. Problem Formulation

The hybrid modeling procedure is divided into two phases, i.e., physics-informed modeling (offline) and data-driven scheme (online), which is based on lots of physics and less data, shown in Figure 1.

2.1. Physics-Informed Modelling

In Figure 2, the ducted fan aerial robot is represented by two coordinate frames, i.e., body-fixed frame ${\mathcal{B}}_f\triangleq \{{\mathcal{B}}_o,{\mathcal{B}}_x,{\mathcal{B}}_y,{\mathcal{B}}_z\}$ and inertial frame ${\mathcal{I}}_f\triangleq \{{\mathcal{I}}_o,{\mathcal{I}}_x,{\mathcal{I}}_y,{\mathcal{I}}_z\}$. ${\mathit{\xi }}^{\mathcal{I}}={[{\xi }_x^{\mathcal{I}},{\xi }_y^{\mathcal{I}},{\xi }_z^{\mathcal{I}}]}^{\mathrm{T}}\in {\mathbb{R}}^3$, ${\mathit{V}}_{\mathcal{I}}={[V_x^{\mathcal{I}},V_y^{\mathcal{I}},V_z^{\mathcal{I}}]}^{\mathrm{T}}\in {\mathbb{R}}^3$, ${\mathit{V}}_{\mathcal{B}}={[V_x^{\mathcal{B}},V_y^{\mathcal{B}},V_z^{\mathcal{B}}]}^{\mathrm{T}}\in {\mathbb{R}}^3$ and ${\mathit{\omega }}_{\mathcal{B}}={[p,q,r]}^{\mathrm{T}}\in {\mathbb{R}}^3$ are position, linear velocity in ${\mathcal{I}}_f$, linear velocity in ${\mathcal{B}}_f$, and angular velocity, respectively. ${\Delta }_{wind}={\left[{\Delta }_x^{\mathcal{B}}\left(t\right),{\Delta }_y^{\mathcal{B}}\left(t\right),{\Delta }_z^{\mathcal{B}}\left(t\right)\right]}^{\mathrm{T}}$ are disturbances acting on the flying robot. $\mathsf{\Theta }={[\varphi ,\theta ,\psi ]}^T\in {\mathbb{R}}^3$ denotes the attitude of the aircraft, where $\varphi$, $\theta$ and $\psi$ are roll, pitch, and yaw angle, respectively. From these rotation angles, the rotational matrix ${\mathit{T}}_{\mathcal{B}\mathcal{I}}\in SO\left(3\right)$ from ${\mathcal{B}}_f$ and ${\mathcal{I}}_f$ can be defined as (see among others, e.g., [9]):

$${\mathit{T}}_{\mathcal{B}\mathcal{I}}=\left[\begin{array}{ccc}{C}_{\psi }{C}_{\theta }& C_{\psi }{S}_{\theta }{S}_{\varphi }-S_{\psi }{C}_{\varphi }& C_{\psi }{S}_{\theta }{C}_{\varphi }+S_{\psi }{S}_{\varphi }\\ S_{\psi }{C}_{\theta }& S_{\psi }{S}_{\theta }{S}_{\varphi }+C_{\psi }{C}_{\varphi }& S_{\psi }{S}_{\theta }{C}_{\varphi }-C_{\psi }{S}_{\varphi }\\ -S_{\theta }& C_{\theta }{S}_{\varphi }& C_{\theta }{C}_{\varphi }\end{array}\right],$$

(1)

where $C_{(*)}=cos(*)$ and $S_{(*)}=sin(*)$, with $(*)$ denotes any angle. The rotational kinematics of the airborne platform is described in the following form:

$$\dot{\mathsf{\Theta }}=\mathsf{\Omega }{\mathit{\omega }}_{\mathcal{B}},\mathsf{\Omega }=\left[\begin{array}{ccc}{C}_{\theta }& 0& -S_{\theta }\\ S_{\theta }{S}_{\varphi }/C_{\varphi }& 1& C_{\theta }{S}_{\varphi }/C_{\varphi }\\ S_{\theta }/C_{\varphi }& 0& C_{\varphi }{C}_{\varphi }\end{array}\right].$$

(2)

Other flight components in Figure 2 are expressed as follows:

$$\begin{array}{cc}\hfill V_{wind}& =\sqrt{{\left(V_x^{\mathcal{B}}-{\Delta }_x^{\mathcal{B}}\right)}^2+{\left(V_y^{\mathcal{B}}-{\Delta }_y^{\mathcal{B}}\right)}^2+{\left(V_z^{\mathcal{B}}-{\Delta }_z^{\mathcal{B}}\right)}^2},\hfill \\ \hfill \alpha & =-arccos\left(\frac{V_z^{\mathcal{B}}-{\Delta }_z^{\mathcal{B}}}{V_{wind}}\right),0\le \alpha \le \pi ,\hfill \\ \hfill \beta & =arctan\left(\frac{V_y^{\mathcal{B}}-{\Delta }_y^{\mathcal{B}}}{V_x^{\mathcal{B}}-{\Delta }_x^{\mathcal{B}}}\right),-\frac{\pi }{2}\le \beta \le \frac{\pi }{2},\hfill \end{array}$$

(3)

where $V_{wind}$, $\alpha$, and $\beta$ are denoted as airspeed, angle of attack, and side-slip angle, respectively. The translational and rotational dynamics of the VTOL aircraft can be formulated as

$$\left\{\begin{array}{c}{\ddot{\mathit{\xi }}}_{\mathcal{I}}={\dot{\mathit{V}}}_{\mathcal{I}}={\mathit{T}}_{\mathcal{B}\mathcal{I}}{\mathit{F}}_{\mathcal{B}}/m+{\mathit{G}}_a,\hfill \\ {\dot{\mathit{\omega }}}_{\mathcal{B}}={\mathit{I}}^{-1}\left({\mathit{M}}_{\mathcal{B}}+\mathit{I}{\mathit{\omega }}_{\mathcal{B}}\times {\mathit{\omega }}_{\mathcal{B}}\right),\hfill \end{array}\right.$$

(4)

where m, ${\mathit{F}}_{\mathcal{B}}$, ${\mathit{M}}_{\mathcal{B}}$, and $\mathit{I}$ are the mass, total force, total moments, and diagonal matrix consisting of the moment of inertia, respectively. Moreover, ${\mathit{G}}_a={\left[\begin{array}{ccc}0\hfill & 0\hfill & g\hfill \end{array}\right]}^{\mathrm{T}}$, and g is the gravitational acceleration. Velocity of the airflow exhausted from the fan ($V_{af}$) can be expressed as:

$$V_{af}=-\frac{V_z^{\mathcal{B}}-{\Delta }_z^{\mathcal{B}}}{2}+\sqrt{{\left(\frac{V_z^{\mathcal{B}}-{\Delta }_z^{\mathcal{B}}}{2}\right)}^2+\frac{T}{2\rho A_{fd}}},$$

(5)

where $T=A_{\omega T}{\omega }_r^2$ is the thrust, and ${\omega }_r$ and $A_{\omega T}$ denote the rotation speed and thrust coefficient, respectively. Moreover, $\rho$ and $A_{fd}$ are the free air density and fan disk area, respectively. Next, we describe the local airspeed on left-wing ($V_{\ell w}$) and right-wing ($V_{rw}$) with local $\alpha$ on left-wing (${\alpha }_{{\ell }_w}$) and right-wing (${\alpha }_{r_w}$) are given as follows:

$$\begin{array}{cc}\hfill V_{{\ell }_w}& =\sqrt{{\left(V_x^{\mathcal{B}}+r{\ell }_w-{\Delta }_x\right)}^2+{\left(V_z^{\mathcal{B}}-{\Delta }_z\right)}^2},\hfill \\ \hfill V_{r_w}& =\sqrt{{\left(V_x^{\mathcal{B}}-r{\ell }_w-{\Delta }_x\right)}^2+{\left(V_z^{\mathcal{B}}-{\Delta }_z\right)}^2},\hfill \\ \hfill {\alpha }_{{\ell }_w}& =C_{{\alpha }_{l_w}}=-\left(V_z^{\mathcal{B}}-{\Delta }_z\right)/V_{{\ell }_w},0\le {\alpha }_{{\ell }_w}\le \pi ,\hfill \\ \hfill {\alpha }_{r_w}& =C_{{\alpha }_{r_w}}=-\left(V_z^{\mathcal{B}}-{\Delta }_z\right)/V_{r_w},0\le {\alpha }_{r_w}\le \pi ,\hfill \end{array}$$

(6)

where ${\ell }_w$ represents the lever arm from the ${\mathcal{B}}_z$ axis to the aerodynamic center of the left and right/left wing. Generally, an aircraft is affected by four main forces during flight (lift, weight, thrust, and drag), depicted in Figure 2c. In the current scenario, the ${\mathit{F}}_{\mathcal{B}}$ can be defined as follows:

$$\begin{array}{cc}\hfill {\mathit{F}}_{\mathcal{B}}& =\underset{\mathrm{Thrust} \mathrm{force}}{\underbrace{\left[\begin{array}{c}0\\ 0\\ -A_{\omega T}{\omega }_r^2\end{array}\right]}}+\underset{\mathrm{Duct}\text{-}\mathrm{body} \mathrm{force}}{\underbrace{V_{wind}^2\left[\begin{array}{c}-\left(A_{L_d}{C}_{\alpha }+A_{D_d}{S}_{\alpha }\right)C_{\beta }\\ -\left(A_{L_d}{C}_{\alpha }+A_{D_d}{S}_{\alpha }\right)S_{\beta }\\ -A_{L_d}{S}_{\alpha }+A_{D_d}{C}_{\alpha }\end{array}\right]}}\hfill \\ & +\underset{\mathrm{Wing} \mathrm{force}}{\underbrace{\left[\begin{array}{c}-\left(A_{L_w}{V}_{{\ell }_w}^2{C}_{{\alpha }_{{\ell }_w}}+A_{L_w}{V}_{r_w}^2{C}_{{\alpha }_{r_w}}+A_{D_w}{V}_{{\ell }_w}^2{S}_{{\alpha }_{{\ell }_w}}+A_{D_w}{V}_{r_w}^2{S}_{{\alpha }_{r_w}}\right)\\ 0\\ -A_{L_w}{V}_{{\ell }_w}^2{S}_{{\alpha }_{{\ell }_w}}-A_{L_w}{V}_{r_w}^2{S}_{{\alpha }_{r_w}}+A_{D_w}{V}_{{\ell }_w}^2{C}_{{\alpha }_{{\ell }_w}}+A_{D_w}{V}_{r_w}^2{C}_{{\alpha }_{r_w}}\end{array}\right]}}\hfill \\ \hfill & +\underset{\mathrm{Momentum} \mathrm{drag}}{\underbrace{\rho A_{fd}{V}_{af}{\zeta }_{ratio}\left[\begin{array}{c}{V}_x^{\mathcal{B}}-{\Delta }_x^{\mathcal{B}}\\ V_y^{\mathcal{B}}-{\Delta }_y^{\mathcal{B}}\\ 0\end{array}\right]}},\hfill \end{array}$$

(7)

where $A_{L_d}$ and $A_{D_d}$ are the lift and drag coefficients, respectively. Moreover, $A_{D_w}$ and $A_{L_w}$ represent the drag and lift coefficients of the half-wing section of the aircraft, respectively. ${\zeta }_{ratio}$ denotes the damping ratio and is defined as

$${\zeta }_{ratio}=\frac{V_{sx}-V_{sx}^{{}^{\prime }}}{V_{sx}}.$$

(8)

The total moment is expressed in the following manner:

$$\begin{array}{cc}\hfill {\mathit{M}}_{\mathcal{B}}& =\underset{\mathrm{Control} \mathrm{surfaces}}{\underbrace{A_{cs}{V}_{af}^2\left[\begin{array}{c}-l_1{\delta }_1+l_1{\delta }_3\\ -l_1{\delta }_2+l_1{\delta }_4\\ l_2{\delta }_1+l_2{\delta }_2+l_2{\delta }_3+l_2{\delta }_4\end{array}\right]}}+\underset{\mathrm{Fan} \mathrm{torque}}{\underbrace{\left[\begin{array}{c}0\\ 0\\ k_{\omega s}{\omega }_r^2\end{array}\right]}}+\underset{\mathrm{Gyroscopic} \mathrm{effect}}{\underbrace{I_f{\omega }_r\left[\begin{array}{c}-q\\ p\\ 0\end{array}\right]}}\hfill \\ & +\underset{\mathrm{Aerodynamic} \mathrm{pitching} \mathrm{moment}}{\underbrace{{ϵ}_{drag}\rho A_{fd}{V}_{af}{\zeta }_{ratio}\left[\begin{array}{c}{V}_x^{\mathcal{B}}-{\Delta }_x^{\mathcal{B}}\\ V_y^{\mathcal{B}}-{\Delta }_y^{\mathcal{B}}\\ 0\end{array}\right]+{ϵ}_D{V}_{wind}^2\left[\begin{array}{c}-\left(A_{L_d}{C}_{\alpha }+A_{D_d}{S}_{\alpha }\right)C_{\beta }\\ -\left(A_{L_d}{C}_{\alpha }+A_{D_d}{S}_{\alpha }\right)S_{\beta }\\ -A_{L_d}{S}_{\alpha }+A_{D_d}{C}_{\alpha }\end{array}\right]}}\hfill \\ & +\underset{\mathrm{Wing} \mathrm{moment}}{\underbrace{\left[\begin{array}{c}0\\ A_M{V}_{{\ell }_w}^2+A_M{V}_{{\ell }_w}^2\\ -{\ell }_w{A}_{L_w}{V}_{{\ell }_w}^2+{\ell }_w{A}_{L_w}{V}_{r_w}^2\end{array}\right]}}+\underset{\mathrm{Anti}-\mathrm{rotation}\mathrm{torque}}{\underbrace{A_{ar}{V}_{af}^2}},\hfill \end{array}$$

(9)

where $k_{\omega s}$, ${ϵ}_{drag}$, ${ϵ}_D$, $I_f$, $A_{ar}$, $l_1$, $l_2$, $A_{cs}$, ${\delta }_i$, and $A_M$ represent the rotation speed with respect to fan torque coefficient, lever arms of momentum drag, duct body force, fan inertia, constant coefficient, lever arms related to roll/pitch-axis, and yaw-axis, deflection angle to force coefficient, the deflection angle of i-th control surface, and moment coefficient of the half-wing section of the vehicle, respectively. Consider the drone’s system defined as follows:

$${\dot{\mathit{x}}}_s=f_s\left({\mathit{x}}_s\left(t\right),{\mathit{u}}_s\left(t\right);{\mathit{\eta }}_s\right),{\mathit{x}}_s\left(0\right)={\mathit{x}}_0,$$

(10)

where ${\mathit{x}}_s={\left[{\mathit{\xi }}_{\mathcal{I}}{\mathit{V}}_{\mathcal{I}}\mathsf{\Theta }{\mathit{\omega }}_{\mathcal{B}}\right]}^{\mathrm{T}}\in {\mathbb{R}}^{12}$ denotes the system states. $f_s$ and ${\mathit{\eta }}_s$ are the dynamics and the aircraft’s parameters. It is considered that the system (10) is directly affected by disturbances and uncertainties. ${\mathit{u}}_s\left(t\right)={\left[T{\mathit{u}}_{s2}^{\mathrm{T}}\right]}^{\mathrm{T}}\in {\mathbb{R}}^4$ represents control input, where ${\mathit{u}}_{s2}$ denotes the regulated moment input and is provided as follows:

$${\mathit{u}}_{s2}=\left[\begin{array}{c}{u}_{11}\\ u_{12}\\ u_{13}\end{array}\right]=\left[\begin{array}{c}-l_1{\delta }_1+l_1{\delta }_3\\ -l_1{\delta }_2+l_1{\delta }_4\\ l_2{\delta }_1+l_2{\delta }_2+l_2{\delta }_3+l_2{\delta }_4\end{array}\right].$$

Next, by considering only the physics-informed model of the aircraft as a nominal system:

$${\dot{\tilde{\mathit{x}}}}_s\left(t\right)=f_s\left({\tilde{\mathit{x}}}_s\left(t\right),{\tilde{\mathit{u}}}_s\left(t\right)\right).$$

(11)

The discrepancy modeling problem is derived from the difference between the quantity of interest (i.e., forces and moments) from a physics-informed model ${\nu }_m\left(t\right)$ and estimated value ${\nu }_o\left(t\right)$. Hence, discrepancy $\Delta \nu \left(t\right)$ is expressed as:

$$\Delta \nu \left(t\right)={\nu }_o\left(t\right)-{\nu }_m\left(t\right).$$

(12)

Under ideal conditions, the estimated value and physics-informed model may have the same magnitude in (12), i.e., $\Delta \nu \left(t\right)\approx 0$. Nonetheless, such an outcome in actual flight may not be possible due to several factors.

2.2. Data-Driven ML-Adaptive SINDy

Inspired by the existing work on learning from discrepancy model [20] and integrating it with adaptive SINDy by [30] for compensating for any abrupt changes in the system dynamics during different flight conditions of the DFAV, a hybrid modeling approach based on physics-informed dynamics and adaptive SINDy is constructed, shown in Figure 3. The sudden changes in the adaptive SINDy involve only addition, deletion, and modifications. This is useful as it enables us to effectively determine the new and missing terms with less data than to identify the entire model from scratch. There are three basic kinds of model changes in the adaptive SINDy:

• Modification: If the whole model is unchanged except for model parameters, then least square regression will be employed on the known model to find new parameters;
• Deletion: If a few terms are removed, then sparse regression can be utilized on the sparse coefficients to identify which terms have been taken out;
• Addition: To find the model error, the SINDy regression will identify the sparsest combination of inactive terms.

2.2.1. Baseline Model

By using the standard procedure of SINDy, we measure the required number of snapshots $m_s$ of ${\mathit{x}}_s$ and ${\mathit{u}}_s$ and rearrange them into two data matrices ${\mathit{X}}_s$ and ${\mathit{U}}_s$:

$$\begin{array}{cc}\hfill {\mathit{X}}_s=& \left[\begin{array}{cccc}{\xi }_x^{\mathcal{I}}\left(t_1\right)& {\xi }_y^{\mathcal{I}}\left(t_1\right)& \cdots & \psi \left(t_1\right)\\ ⋮& ⋮& \ddots & ⋮\\ {\xi }_x^{\mathcal{I}}\left(t_i\right)& {\xi }_y^{\mathcal{I}}\left(t_i\right)& \cdots & \psi \left(t_i\right)\\ ⋮& ⋮& \ddots & ⋮\\ {\xi }_x^{\mathcal{I}}\left(t_m\right)& {\xi }_y^{\mathcal{I}}\left(t_m\right)& \cdots & \psi \left(t_m\right)\end{array}\right],\hfill \\ \hfill {\mathit{U}}_s=& \left[\begin{array}{cccc}T\left(t_1\right)& u_{11}\left(t_1\right)& \cdots & u_{13}\left(t_1\right)\\ ⋮& ⋮& \ddots & ⋮\\ T\left(t_i\right)& u_{11}\left(t_i\right)& \cdots & u_{13}\left(t_i\right)\\ ⋮& ⋮& \ddots & ⋮\\ T\left(t_m\right)& u_{11}\left(t_m\right)& \cdots & u_{13}\left(t_m\right)\end{array}\right],\hfill \end{array}$$

(13)

where $t_i$ is the sampling time. Simplify Equation (13) and transform it into the following form:

$$\begin{array}{cc}\hfill {\mathit{X}}_s& ={\left[{\mathit{x}}_{s1}\left(t_1\right){\mathit{x}}_{s2}\left(t_2\right){\mathit{x}}_{s3}\left(t_3\right)...{\mathit{x}}_{ms}\left(t_m\right)\right]}^{\mathrm{T}},\hfill \\ \hfill {\mathit{U}}_s& ={\left[{\mathit{u}}_{s1}\left(t_1\right){\mathit{u}}_{s2}\left(t_2\right){\mathit{u}}_{s3}\left(t_3\right)...{\mathit{u}}_{ms}\left(t_m\right)\right]}^{\mathrm{T}},\hfill \end{array}$$

(14)

where ${\mathit{u}}_{ms}$ and ${\mathit{x}}_{ms}$ are control and state vector at the mth sampling time. The value of the aircraft’s state is calculated by the numerical differentiation.

$$\begin{array}{cc}\hfill {\dot{\mathit{X}}}_s& ={\left[{\dot{\mathit{x}}}_{s1}\left(t_1\right){\dot{\mathit{x}}}_{s2}\left(t_2\right){\dot{\mathit{x}}}_{s3}\left(t_3\right)...{\dot{\mathit{x}}}_{ms}\left(t_m\right)\right]}^{\mathrm{T}},\hfill \\ & =\left[\begin{array}{cccc}{\dot{\xi }}_x^{\mathcal{I}}\left(t_1\right)& {\dot{\xi }}_y^{\mathcal{I}}\left(t_1\right)& \cdots & \dot{\psi }\left(t_1\right)\\ ⋮& ⋮& \ddots & ⋮\\ {\dot{\xi }}_x^{\mathcal{I}}\left(t_i\right)& {\dot{\xi }}_y^{\mathcal{I}}\left(t_i\right)& \cdots & \dot{\psi }\left(t_i\right)\\ ⋮& ⋮& \ddots & ⋮\\ {\dot{\xi }}_x^{\mathcal{I}}\left(t_m\right)& {\dot{\xi }}_y^{\mathcal{I}}\left(t_m\right)& \cdots & \dot{\psi }\left(t_m\right)\end{array}\right].\hfill \end{array}$$

(15)

The vehicle model depends primarily on the physics-informed modeling method with less concentration on the ML technique as it is mainly useful for any discrepancy and abrupt changes in the parametric values. Therefore, a similar assumption to [29] is required, i.e., either the obtained data are prefiltered through relevant frameworks [37,38] or it contains no noise interference. Next, a candidate library function $\Xi ({\mathit{X}}_s,{\mathit{U}}_s)$ is developed using the data matrices ${\mathit{X}}_s$ and ${\mathit{U}}_s$, in which selected function is arbitrary, and is either a trigonometric or polynomial function. In general, polynomial functions are a relatively basic choice. Therefore, the complexity of the library can be increased by including trigonometric functions. Nevertheless, both functions and physics-informed knowledge of the vehicle will be utilized to construct the library function:

$$\Xi ({\mathit{X}}_s,{\mathit{U}}_s)=\left[{\mathit{1}}^{\mathrm{T}},{\mathit{X}}_s^{\mathrm{T}},{({\mathit{X}}_s\otimes {\mathit{X}}_s)}^{\mathrm{T}},{({\mathit{X}}_s\otimes {\mathit{U}}_s)}^{\mathrm{T}},\cdots ,S_{{\mathit{X}}_s^{\mathrm{T}}},S_{{\mathit{U}}_s^{\mathrm{T}}},\cdots \right],$$

(16)

where ${\mathit{X}}_s\otimes {\mathit{X}}_s$ and ${\mathit{X}}_s\otimes {\mathit{U}}_s$ are defined as:

$$\begin{array}{cc}\hfill & {\mathit{X}}_s\otimes {\mathit{X}}_s=\left[\begin{array}{cccccc}{\xi }_x^{\mathcal{I}}{\xi }_x^{\mathcal{I}}\left(t_1\right)& {\xi }_x^{\mathcal{I}}{\xi }_y^{\mathcal{I}}\left(t_1\right)& \cdots & {\xi }_y^{\mathcal{I}}{\xi }_y^{\mathcal{I}}\left(t_1\right)& \cdots & {\psi }^2\left(t_1\right)\\ {\xi }_x^{\mathcal{I}}{\xi }_x^{\mathcal{I}}\left(t_2\right)& {\xi }_x^{\mathcal{I}}{\xi }_y^{\mathcal{I}}\left(t_2\right)& \cdots & {\xi }_y^{\mathcal{I}}{\xi }_y^{\mathcal{I}}\left(t_2\right)& \cdots & {\psi }^2\left(t_2\right)\\ ⋮& ⋮& \cdots & ⋮& \cdots & ⋮\\ {\xi }_x^{\mathcal{I}}{\xi }_x^{\mathcal{I}}\left(t_m\right)& {\xi }_x^{\mathcal{I}}{\xi }_y^{\mathcal{I}}\left(t_m\right)& \cdots & {\xi }_y^{\mathcal{I}}{\xi }_y^{\mathcal{I}}\left(t_m\right)& \cdots & {\psi }^2\left(t_m\right)\end{array}\right],\hfill \\ \hfill & {\mathit{X}}_s\otimes {\mathit{U}}_s=\left[\begin{array}{cccccc}{\xi }_x^{\mathcal{I}}T\left(t_1\right)& {\xi }_x^{\mathcal{I}}{u}_{11}\left(t_1\right)& \cdots & {\xi }_y^{\mathcal{I}}{u}_{11}\left(t_1\right)& \cdots & \psi u_{13}\left(t_1\right)\\ {\xi }_x^{\mathcal{I}}T\left(t_2\right)& {\xi }_x^{\mathcal{I}}{u}_{11}\left(t_2\right)& \cdots & {\xi }_y^{\mathcal{I}}{u}_{11}\left(t_2\right)& \cdots & \psi u_{13}\left(t_2\right)\\ ⋮& ⋮& \cdots & ⋮& \cdots & ⋮\\ {\xi }_x^{\mathcal{I}}T\left(t_m\right)& {\xi }_x^{\mathcal{I}}{u}_{11}\left(t_m\right)& \cdots & {\xi }_y^{\mathcal{I}}{u}_{11}\left(t_m\right)& \cdots & \psi u_{13}\left(t_m\right)\end{array}\right].\hfill \end{array}$$

(17)

The dynamics of the vehicle model only depend on the few nonlinear terms in practice. The sparse model of $f_k(.)$ in system (10) can be set as follows:

$$f_k({\mathit{X}}_s,{\mathit{U}}_s,{\chi }_k)=\sum _{j=0}^{p_s}{\chi }_{kj}{\mathsf{\Lambda }}_j({\mathit{x}}_s,{\mathit{u}}_s),k=1,2,3,\cdots ,12,$$

(18)

where $p_s+1$ represents the number of candidate functions. ${\chi }_{kj}$ and ${\mathsf{\Lambda }}_j$ are the candidate functions of the j-th column in $\Xi ({\mathit{X}}_s,{\mathit{U}}_s)$ and the weighted coefficient of ${\mathsf{\Lambda }}_j({\mathit{x}}_s,{\mathit{u}}_s)$ related to the k-th state, i.e., ${\mathit{\chi }}_k={[{\chi }_{k0},{\chi }_{k1},{\chi }_{k2},\cdots ,{\chi }_{k{p}_s}]}^{\mathrm{T}}$, respectively. In the third step, sparse optimization is employed to determine a sparse model:

$${\mathit{\chi }}_k=\underset{{\mathit{\chi }}_k}{argmin}{∥{\dot{\mathit{X}}}_{ks}-{\mathit{\chi }}_k\Xi {({\mathit{X}}_s,{\mathit{U}}_s)}^{\mathrm{T}}∥}_2+{\lambda }_k{∥{\mathit{\chi }}_k∥}_1,$$

(19)

where ${\dot{\mathit{X}}}_{ks}$ and ${\lambda }_k$ denote the k-th row of $\dot{\mathit{X}}$ and sparsity-promoting hyper-parameter, respectively. The theoretical results for convergences of SINDy framework have been provided in [39]. The computed value of ${\mathit{\chi }}_k$ is stored in sparse matrix $\mathsf{\Psi }$:

$$\mathsf{\Psi }=\left[\begin{array}{cccccc}{\chi }_{01}& {\chi }_{02}& \cdots & {\chi }_{0k}& \cdots & {\chi }_{0n}\\ {\chi }_{11}& {\chi }_{12}& \cdots & {\chi }_{1k}& \cdots & {\chi }_{1n}\\ ⋮& ⋮& \ddots & ⋮& \ddots & ⋮\\ {\chi }_{p1}& {\chi }_{p2}& \cdots & {\chi }_{pk}& \cdots & {\chi }_{pn}\end{array}\right]=\left[{\mathit{\chi }}_1{\mathit{\chi }}_2\cdots {\mathit{\chi }}_k\cdots {\mathit{\chi }}_{12}\right].$$

(20)

The baseline model of the DFAV’s dynamics can be written in the following expression:

$${\mathit{X}}_s\approx \Xi ({\mathit{X}}_s,{\mathit{U}}_s)\mathsf{\Psi }.$$

(21)

To summarize the entire process, the baseline model is determined, and a grid search is employed to identify the optimal hyper-parameter selection. In grid-search, all the combinations of hyper-parameters are evaluated, and the best available set is chosen. This process is executed one time to find the best set of hyper-parameters for future updates. The baseline model can be expressed as sparse coefficients in ${\mathit{\chi }}_0$.

2.2.2. Estimation of Model Divergence

Due to various factors (e.g., flight condition, nonlinearity, turbulence, actuator faults, etc.), the model parameters may vary at different times, resulting in divergence problems. Hence, a prediction framework that can quickly estimate any variation in the vehicle model is needed. In this work, the classical-predictor–correction approach is utilized to detect the variation caused by the aforementioned factors [40]. The predictor step is executed over a time $\overline{\tau }$ during the interval $t,t+\overline{\tau }$ using a valid model at time t. The divergence of the predictor ${\widehat{\mathit{x}}}_s$ and measured state ${\mathit{x}}_s$ is calculated at $t+\tau$ as $\parallel {\mathit{x}}_{div}\parallel =\parallel {\widehat{\mathit{x}}}_s(t+\tau )-{\mathit{x}}_s(t+\tau )\parallel$. The main theme is to determine when the model and the related state measurements diverge faster as opposed to the predicted ones by the dynamics of the system. Among many techniques, the divergence of the trajectory can be expressed by the Lyapunov exponent:

$$\lambda =\underset{\tau \to \infty }{lim}\underset{{\mathit{x}}_{div}\left(t_0\right)\to 0}{lim}\frac{〈log\left(\frac{{\mathit{x}}_{div}\left(t_0+\tau \right)}{{\mathit{x}}_{div}\left(t_0\right)}\right)〉}{\tau },$$

(22)

and its inverse sets of fastest time scale [30]. Since the information about the dynamical system for the prediction step, the Lyapunov exponent is computed via tangent space. The requirement for an effective detection necessitates computing tolerance ${\lambda }_{div}$ and determining the divergence time. If the prediction horizon by the local Lyapunov exponent and provided time scale are inconsistent, a divergence issue will arise between the measurements and the model. In the start, when the predicted and measured values deviate, then the divergence horizon is defined as:

$$T_{div}\left(t\right)=\underset{t_{div}}{argmin}\parallel {\widehat{\mathit{x}}}_s(t+t_{div})-{\mathit{x}}_s(t+t_{div})\parallel >∥{\lambda }_{div}∥.$$

(23)

Hence, it is assumed that, at time t, the real system and prediction model diverge if the model value and measured time scale differ:

$${\overline{\lambda }}_{div}\left(t\right)>\frac{log\left({\lambda }_{div}\right)-log(\parallel {\mathit{x}}_{div}\left(t\right)\parallel )}{T_{div}\left(t\right)},$$

(24)

where ${\lambda }_{div}$ is the maximum value obtained from (22) during the time scale $[t,t+T_{div}]$, & ${\overline{\lambda }}_{div}={\langle {\lambda }_{div}\left(t^{{}^{\prime }}\right)\rangle }_{t^{{}^{\prime }}\in [t,t+T_{div}]}$ with ${\lambda }_{div}\left(t\right)=\max ({\lambda }_{di{v}_1},{\lambda }_{di{v}_2},{\lambda }_{di{v}_3},\cdots ,{\lambda }_{di{v}_{12}})$.

2.2.3. Adaptive Model Recovery

The following process is implemented for the quick recovery of model parameters. In the beginning, new data onto the existing sparse data ${\mathit{\chi }}_0$ are regressed to determine varying parameters. Thereafter, the deletion of terms is identified by executing the sparse regression on the columns of $\Xi$ that in ${\mathit{\chi }}_0$ corresponds to the nonzero rows. In the presence of a residual error, a sparse model in the inactive columns of $\Xi$ that corresponds to zero rows in ${\mathit{\chi }}_0$ is fit for this error. Through this process, new terms may be introduced. This procedure is only performed when there is a divergence issue present.

2.3. Control Objective

By considering the physics-informed model as a nominal system ${\tilde{\mathit{x}}}_s$ and real system ${\mathit{x}}_s$ obtained as a result of the physics-informed ML procedure, suppose the aerial vehicle is flying with speed $\upsilon$, and the plane follows a reference trajectory ${\mathit{x}}_{\mathrm{ref}}$ in the presence of different factors (e.g., disturbances, uncertainties, actuator faults). From a control perspective, the tracking error of both control input ${\mathit{e}}_u={\mathit{u}}_s-{\mathit{u}}_{\mathrm{ref}}$ and aircraft’s state ${\mathit{e}}_x={\mathit{x}}_s-{\mathit{x}}_{\mathrm{ref}}$ needs to be minimized. Similarly, ${\tilde{\mathit{e}}}_x$ and ${\tilde{\mathit{e}}}_u$ represent nominal state and control input tracking error, respectively. Hence, the cost function is designed as follows:

$$\begin{array}{cc}\hfill J\left({\tilde{\mathit{e}}}_x\left(t\right),{\tilde{\mathit{e}}}_u\left(t\right)\right)& ={\int }_t^{t+t_p}\underset{\mathrm{Stage} \cos \mathrm{t}}{\underbrace{{\mathcal{S}}_c\left({\tilde{\mathit{e}}}_x\left(\tau \right|t),{\tilde{\mathit{e}}}_u\left(\tau \right|t)\right)d\tau }}+\underset{\mathrm{Terminal} \mathrm{penalty}}{\underbrace{{\mathcal{T}}_p\left({\tilde{\mathit{e}}}_x(t+t_p|t)\right)}},\hfill \\ \hfill J\left({\tilde{\mathit{e}}}_x\left(t\right),{\tilde{\mathit{e}}}_u\left(t\right)\right)& ={\int }_t^{t+t_p}\underset{\mathrm{Stage} \mathrm{cost}}{\underbrace{\left({∥{\tilde{\mathit{e}}}_x\left(\tau \right|t)∥}_{\mathit{Q}}^2+{∥{\tilde{\mathit{e}}}_u\left(\tau \right|t)∥}_{\mathit{P}}^2\right)d\tau }}+\underset{\mathrm{Terminal} \mathrm{penalty}}{\underbrace{{∥{\tilde{\mathit{e}}}_x(t+t_p|t)∥}_{\mathit{R}}^2}},\hfill \end{array}$$

(25)

where $t_p=n_s{t}_i$ represents the prediction horizon with $n_s$ denoting the predictive step. Moreover, $\mathit{P}$, $\mathit{Q}$, and $\mathit{R}$ are the positive definitive weighting matrices in (25).

For the smooth and effective functioning of controller design, some preliminaries are needed before proceeding to the main control framework. Recall the nominal system (11), where only physics-informed modeling is considered.

Assumption 1.

It is assumed that $f_s\left({\tilde{\mathit{x}}}_s\left(t\right),{\tilde{\mathit{u}}}_s\left(t\right)\right)$ is Lipschitz continuous with ${\tilde{\mathit{u}}}_s\in {\mathbb{U}}_{\mathrm{MPC}}$ is locally Lipschitz with a constant $\mathcal{L}$ in ${\tilde{\mathit{x}}}_s\left(t\right)$. Hence, the following condition $∥f_s\left({\tilde{\mathit{x}}}_{s1},{\tilde{\mathit{u}}}_s\right)-f_s\left({\tilde{\mathit{x}}}_{s2},{\tilde{\mathit{u}}}_s\right)∥$$⩽\mathcal{L}∥{\tilde{\mathit{x}}}_{s1}-{\tilde{\mathit{x}}}_{s2}∥$ holds.

Assumption 2.

We assume that nominal system (11) is known and differentiable with all the initial states available. Thus, one can use Jacobian linearization at the origin, which is assumed to be stabilizable. With this assumption, a linear state feedback control law ${\mathit{u}}_s=\mathit{K}{\tilde{\mathit{e}}}_x$ can be obtained in such a form that is asymptotically stable if no disturbances are present.

To further elaborate on the nature of the control system, Assumption A2 from [16] with an additional condition can be adopted in the following assumption:

Assumption 3.

(i) For the nominal tracking error, the terminal controller and terminal region$\mathbb{T}$are such that, if${\tilde{\mathit{e}}}_x\left(t_l+t_p|t_l\right)\in \mathbb{T}$, then by applying${\tilde{\mathit{u}}}_s\left(\tau |t_{l+1}\right)=\mathit{K}{\tilde{\mathit{e}}}_x\left(\tau |t_{l+1}\right)$to the vehicle over the time interval$\tau \in \left[t_l+t_p,t_{l+1}+t_p\right)$, it satisfies the following expressions of the locally stabilizing controller:

$${\tilde{\mathit{e}}}_x\left(\tau |t_l\right)\in \mathbb{T},{\tilde{\mathit{u}}}_s=\mathit{K}{\tilde{\mathit{e}}}_x\in {\mathbb{U}}_{\mathrm{MPC}}.$$

(26)

(ii) the condition related to stage and terminal cost is provided as follows:

$${\dot{\mathcal{T}}}_p\left({\tilde{\mathit{e}}}_x\left(\tau |t_l\right)\right)+{\mathcal{S}}_c\left({\tilde{\mathit{e}}}_x\left(\tau |t_l\right),{\tilde{\mathit{e}}}_u\left(\tau |t_l\right)\right)\le 0.$$

(27)

Remark 1.

Assumption 3 in (26) is a standard assumption employed in the MPC community [41]. However, in the current scenario, its usage is only restricted to an ideal situation, and the linear state feedback control law in Assumption 2 is only available if the tracking error lies in the terminal region, which is not always the case in the control system for an aircraft.

Several techniques to determine the terminal region have been developed over the years (e.g., [42,43,44]). Inspired by [16], a terminal region is defined as:

$$\mathbb{T}=\left\{{\tilde{\mathit{e}}}_x:{\tilde{\mathcal{F}}}_1|{\tilde{\mathit{e}}}_{\xi }|+{\tilde{\mathcal{F}}}_2|{\tilde{\mathit{e}}}_V|+{\tilde{\mathcal{F}}}_3|{\tilde{\mathit{e}}}_{\mathsf{\Theta }}|+{\tilde{\mathcal{F}}}_4\left|{\tilde{\mathit{e}}}_{\omega }\right|<1-\varkappa \right\},$$

(28)

where $\varkappa$ can be appropriately selected based on the control design. $\tilde{\mathcal{F}}$ and ${\tilde{\mathit{e}}}_{(•)}$ are the feedback gains and state error between reference and actual trajectory with parameters satisfying $p_{\iota }{q}_{\iota }\le 0.25$, ${\tilde{\mathcal{F}}}_{\iota }\in \left(\frac{1\pm \sqrt{1-4p_{\iota }{q}_{\iota }}}{p_{\iota }{q}_{\iota }}\right)$, $\iota \in \{1,2,3,...,12\}$.

Furthermore, disturbances are also present in the system. Therefore, an assumption is needed to prove the recursive feasibility and stability in the subsequent section.

Assumption 4.

It is assumed that turbulence ${\Delta }_{\mathrm{ext}}$ may also be present along with ${\Delta }_{wind}$. Hence, the dynamics of total disturbances acting on the aircraft are defined as:

$$\begin{array}{cc}\hfill {\dot{\Delta }}_{\mathrm{total}}& ={\Delta }_{\mathrm{ext}}+{\Delta }_{wind},\hfill \\ \hfill & ={\left[{\Delta }_v{\Delta }_{\omega }\right]}^{\mathrm{T}}+{\Delta }_{wind},\hfill \end{array}$$

(29)

where ${\Delta }_v$ and ${\Delta }_{\omega }$ denote the turbulence acting on the translational and rotational components of the vehicle. It is assumed that the total disturbances are bounded, i.e., $\parallel {\Delta }_{\mathrm{total}}\parallel \le \gimel$, which can be fast-time varying unlike [16].

Furthermore, crucial definitions related to stability proof are given as follows:

Definition 1.

[45] The tracking error system is input-to-state stable (ISS), if there exists a $\mathcal{K}\mathcal{L}$ function $\mathsf{\Gamma }(·,·):{\mathbb{R}}_{\ge 0}\times {\mathbb{R}}_{\ge 0}\to \mathbb{R}$ & $\mathcal{K}$ function $\mathsf{{\rm Y}}(·)$ such that $t\ge 0$, it satisfies that

$$∥{\mathit{e}}_x\left(t\right)∥\le \mathsf{\Gamma }\left(∥{\mathit{e}}_x\left(t_0\right)∥,t\right)+\mathsf{{\rm Y}}(\gimel ).$$

(30)

Definition 2.

[45] A function$V(.)$is referred to as ISS-Lyapunov function for tracking error system, if there exist${\mathcal{K}}_{\infty }$functions${\beth }_i(.)$, and $\mathcal{K}$function$\daleth (.)$such that$\mathsf{\forall }{\mathit{e}}_x\in {\mathbb{R}}^2$

$${\beth }_1\left(∥{\mathit{e}}_x\left(t_l\right)∥\right)\le V\left({\mathit{e}}_x\left(t_l\right)\right)\le {\beth }_2\left(∥{\mathit{e}}_x\left(t_l\right)∥\right),$$

(31)

$$V\left({\mathit{e}}_x\left(t_{l+1}\right)\right)-V\left({\mathit{e}}_x\left(t_l\right)\right)\le -{\beth }_3\left(∥{\mathit{e}}_x\left(t_l\right)∥\right)+\daleth (\gimel ).$$

(32)

If the Definitions 1 and 2 satisfy and no disturbances are acting on the aerial vehicle, then the tracking error also vanishes [46].

3. Control Framework

This section explains the presented control approach, illustrated in Figure 4. Reference trajectory $\mathcal{P}$ is employed to generate a particular flight envelop in different scenarios during the time interval $t\in [0,10]$, which is expressed as:

$$\mathcal{P}\left(t\right)={\left[{\mathit{\xi }}_{\mathrm{ref}},{\mathit{V}}_{\mathrm{ref}},{\mathit{\eta }}_{\mathrm{ref}},{\mathit{\omega }}_{\mathrm{ref}}\right]}^{\mathrm{T}}.$$

(33)

The optimization control problem is solved at a time sequence $\left\{t_l|l\in \mathbb{N},t_{l+1}-t_l=t_i\right\}$:

Problem 1.

$$\begin{array}{cc}\hfill \underset{{\tilde{\mathit{u}}}_s\left(\tau |t_l\right)}{\min }& J\left({\tilde{\mathit{e}}}_x\left(t\right),{\tilde{\mathit{e}}}_u\left(t\right)\right),\hfill \end{array}$$

(34a)

$$\begin{array}{cc}\hfill \mathrm{subject}\mathrm{to}& {\tilde{\mathit{x}}}_s\left(t_l|t_l\right)={\mathit{x}}_s\left(t_l\right),\hfill \end{array}$$

(34b)

$$\begin{array}{cc}\hfill & {\dot{\tilde{\mathit{x}}}}_s\left(\tau |t_l\right)=f_s\left({\tilde{\mathit{x}}}_s\left(\tau |t_l\right),{\tilde{\mathit{u}}}_s\left(\tau |t_l\right)\right),\hfill \end{array}$$

(34c)

$$\begin{array}{cc}\hfill & {\tilde{\mathit{u}}}_s\left(\tau |t_l\right)\in {\mathbb{U}}_{\mathrm{MPC}},\hfill \end{array}$$

(34d)

$$\begin{array}{cc}\hfill & {\tilde{\mathit{x}}}_s\left(\tau |t_l\right)\in \mathbb{X},\hfill \end{array}$$

(34e)

$$\begin{array}{cc}\hfill & ∥{\tilde{\mathit{e}}}_x\left(\tau |t_l\right)∥\le \frac{t_p\times \varpi }{\tau -t_l},\hfill \end{array}$$

(34f)

$$\begin{array}{cc}\hfill & {\tilde{\mathit{e}}}_x\left(t_l+t_p|t_l\right)\in {\mathbb{T}}_{\vartheta },\hfill \end{array}$$

(34g)

where ${\mathbb{T}}_{\vartheta }=\left\{{\tilde{\mathit{e}}}_x:∥{\tilde{\mathit{e}}}_x∥\le \vartheta \right\}$ is a robust terminal region, and $\varpi$ is designed as follows:

$$\varpi =\frac{\upsilon -\varkappa }{{({\mathcal{F}}_1^2+{\mathcal{F}}_2^2+{\mathcal{F}}_3^2+{\mathcal{F}}_4^2)}^{1/2}},\vartheta <\varpi ,$$

(35)

Remark 2.

For Problem 1, the nominal system (physics-informed dynamics) is updated by the actual state (physics-informed ML model) at each step. Consequently, the provided optimization problem must be solved online. Nevertheless, such an updating framework produces an optimal state with respect to the current one, and since the disturbances are not considered in Problem 1, the designed optimization problem consumes less computational resources.

Problem 1 determines the minimizing sequence over the interval $[t_l,t_l+t_p)$:

$${\tilde{\mathit{u}}}_s^{*}\left(t_l\right)=\left\{{\tilde{\mathit{u}}}_s^{*}\left(t_l|t_l\right),{\tilde{\mathit{u}}}_s^{*}\left(t_{l+1}|t_l\right),\dots ,{\tilde{\mathit{u}}}_s^{*}\left(t_{l+N}|t_l\right)\right\}.$$

(36)

Hence, the drone’s control signal over interval $\left[t_l,t_{l+1}\right)$ is defined as:

$${\mathit{u}}_s={\tilde{\mathit{u}}}_s^{*}\left(t_l|t_l\right),$$

(37)

where ${\tilde{\mathit{u}}}_s^{*}\left(t_l|t_l\right)$ is the initial control action of ${\tilde{\mathit{u}}}_s^{*}\left(t_l\right)$.

A systematic procedure to apply the control action to the aircraft is provided in Algorithm 1. Nonetheless, if we assume that unknown disturbances are present, a tracking error between the reference and actual trajectory may exist. Furthermore, in order to satisfy the state constraints in the presence of the disturbances, the error on the upper bound between the actual and predicted nominal state is given as follows:

$$\begin{array}{cc}\hfill & ∥{\mathit{x}}_s\left(t_{l+1}\right)-{\tilde{\mathit{x}}}_s^{*}\left(t_{l+1}|t_l\right)∥\hfill \\ & =\parallel {\mathit{x}}_s\left(t_l\right)+{\int }_{t_l}^{t_{l+1}}{f}_s\left({\mathit{x}}_s\left(\tau \right),{\tilde{\mathit{u}}}_s\left(\tau \right)\right)d\tau -{\tilde{\mathit{x}}}_s^{*}\left(t_l|t_l\right)\hfill \\ & -{\int }_{t_l}^{t_{l+1}}{f}_s\left({\tilde{\mathit{x}}}_s^{*}\left(\tau |t_l\right),{\tilde{\mathit{u}}}_s^{*}\left(t_l|t_l\right)\right)d\tau \parallel ,\hfill \\ & \le \gimel t_i+\upsilon {\int }_{t_l}^{t_{l+1}}∥{\mathit{x}}_s\left(\tau \right)-{\tilde{\mathit{x}}}_s^{*}\left(\tau |t_l\right)∥d\tau ,\hfill \\ & \le \gimel t_i{e}^{\upsilon t_i}.\hfill \end{array}$$

(38)

Algorithm 1: MPC-based control using physics-informed ML

Theorem 1.

Suppose the Assumptions 2 and 3 hold, and the Problem 1 is feasible at initial time $t_0$. Then, Problem 1 is feasible for all intervals under Algorithm 1, if the following conditions satisfy:

$$\gimel \le \frac{e^{-\upsilon t_p}}{t_i}(\varpi -\vartheta ),\tilde{\mathcal{F}}{t}_i\ge ln\left(\varpi \right)-ln\left(\vartheta \right),$$

(39)

where$\tilde{\mathcal{F}}=min\left\{{\tilde{\mathcal{F}}}_i\right\}$.

$$\vartheta \ge \frac{\varpi (t_p-t_i)}{t_p}.$$

(40)

Proof. 

A feasible control sequence ${\tilde{\mathit{u}}}_s\left(\tau |t_{l+1}\right)$ at $t_{l+1}$ is provided as follows:

$${\tilde{\mathit{u}}}_s\left(\tau |t_{l+1}\right)=\left\{\begin{array}{c}{\tilde{\mathit{u}}}_s^{*}\left(\tau |t_l\right),\tau \in \left[t_{l+1},t_l+t_p\right),\hfill \\ \mathit{K}{\tilde{\mathit{e}}}_x\left(\tau |t_l\right),\tau \in \left[t_l+t_p,t_{l+1}+t_p\right).\hfill \end{array}\right.$$

(41)

As the nominal state is updated by the actual state, i.e., ${\tilde{\mathit{x}}}_s\left(t_{l+1}|t_{l+1}\right)={\mathit{x}}_s\left(t_{l+1}\right)$, and by considering the time interval $\tau \in \left[t_{l+1},t_l+t_p\right)$, we obtain

$$\begin{array}{cc}\hfill & ∥{\tilde{\mathit{x}}}_s\left(\tau |t_{l+1}\right)-{\tilde{\mathit{x}}}_s^{*}\left(\tau |t_l\right)∥\hfill \\ & =\parallel {\mathit{x}}_s\left(t_{l+1}\right)+{\int }_{t_{l+1}}^{\tau }{f}_s\left({\tilde{\mathit{x}}}_s\left(s|t_{l+1}\right),{\tilde{\mathit{u}}}_s^{*}\left(s|t_l\right)\right)ds\hfill \\ & -{\tilde{\mathit{x}}}_s^{*}\left(t_{l+1}|t_l\right)-{\int }_{t_{l+1}}^{\tau }{f}_s\left({\tilde{\mathit{x}}}_s^{*}\left(s|t_l\right),{\tilde{\mathit{u}}}_s^{*}\left(s|t_l\right)\right)ds\parallel .\hfill \end{array}$$

(42)

Utilizing Grönwall–Bellman inequality, we achieve:

$$∥{\tilde{\mathit{x}}}_s\left(\tau |t_{l+1}\right)-{\tilde{\mathit{x}}}_s^{*}\left(\tau |t_l\right)∥\le \gimel t_i{e}^{\upsilon \left(\tau -t_{l+1}+t_i\right)}.$$

(43)

Employing triangle inequality and substituting $t_l+t_p$ in (43), we have

$$∥{\tilde{\mathit{e}}}_x\left(t_l+t_p|t_{l+1}\right)∥\le ∥{\tilde{\mathit{e}}}_x^{*}\left(t_l+t_p|t_l\right)∥+\gimel t_i{e}^{\upsilon t_i}.$$

(44)

Since the ensuing expressions hold, $∥{\tilde{\mathit{e}}}_x^{*}\left(t_l+t_p|t_l\right)∥\le \vartheta$ and $\gimel \le \frac{e^{-\upsilon t_p}}{t_i}(\varpi -\vartheta )$. Therefore, one can obtain:

$$∥{\tilde{\mathit{e}}}_x\left(t_l+t_p|t_{l+1}\right)∥\le \varpi ,$$

(45)

which implies ${\tilde{\mathit{e}}}_x\left(t_l+t_p|t_{l+1}\right)\in \mathbb{T}$. Next, consider the control input $\mathit{u}\left(\tau \right)=\mathit{K}{\tilde{\mathit{e}}}_x\in {\mathbb{U}}_{\mathrm{MPC}}$ is directly applied to the drone during $\tau \in \left[t_l+t_p,t_{l+1}+t_p\right)$.

$$\begin{array}{cc}\hfill \frac{d}{d\tau }{∥{\tilde{\mathit{e}}}_x\left(\tau \right|t_{l+1})∥}^2& =-2({\tilde{\mathcal{F}}}_1{\tilde{e}}_{\xi }^2+{\tilde{\mathcal{F}}}_2{\tilde{e}}_V^2+{\tilde{\mathcal{F}}}_3{\tilde{e}}_{\mathsf{\Theta }}^2+{\tilde{\mathcal{F}}}_4{\tilde{e}}_{\omega }^2),\hfill \\ \hfill & \le -2\tilde{\mathcal{F}}{∥{\tilde{\mathit{e}}}_x\left(\tau \right|t_{l+1})∥}^2.\hfill \end{array}$$

(46)

By using the comparison principle,

$$∥{\tilde{\mathit{e}}}_x\left(t_{l+1}+t_p|t_{l+1}\right)∥\le ∥{\tilde{\mathit{e}}}_x\left(t_l+t_p|t_{l+1}\right)∥e^{-\widetilde{\mathcal{F}{t}_i}}.$$

Through $\tilde{\mathcal{F}}{t}_i\ge ln\left(\varpi \right)-ln\left(\vartheta \right)$:

$$∥{\tilde{\mathit{e}}}_x\left(t_{l+1}+t_p|t_{l+1}\right)∥\le \vartheta .$$

(47)

This comprehensively ensures that ${\tilde{\mathit{u}}}_s\left(\tau |t_{l+1}\right)$ is able to force the nominal tracking error into the terminal region ${\mathbb{T}}_{\vartheta }$, which satisfies the constraint (34g).

To prove the constraint satisfaction of (34f), we initially consider the time interval $\tau \in [t_{l+1},t+t_p)$. From (43), it follows that:

$$∥{\tilde{\mathit{e}}}_x\left(\tau |t_{l+1}\right)∥\le ∥{\tilde{\mathit{e}}}_x^{*}\left(\tau |t_l\right)∥+\gimel t_i{e}^{\upsilon t_i}.$$

(48)

Because of (39) and $∥{\tilde{\mathit{e}}}_x\left(\tau |t_l\right)∥\le \frac{t_p}{\tau -t_l}$, we obtain

$$∥{\tilde{\mathit{e}}}_x\left(\tau |t_{l+1}\right)∥\le \frac{t_p\times \varpi }{\tau -t_l}+(\varpi -\vartheta ).$$

(49)

From (40), we have

$$\varpi -\vartheta \le \frac{t_i\varpi }{t_p-t_i}\le \frac{t_i\varpi t_p}{(\tau -t_{l+1})(\tau -t_l)}.$$

(50)

Substituting (50) into (49), we achieve

$$∥{\tilde{\mathit{e}}}_x\left(\tau |t_{l+1}\right)∥\le \frac{t_p\times \varpi }{\tau -t_{l+1}},$$

(51)

which comprehensively proves that that state constraint (34f) is satisfied. Next, we consider the time interval $\tau \in [t_k+t_p,t_{l+1}+t_p)$. Through the proof of (45) and (47), and since $\frac{\varpi t_i}{\tau -t_{l+1}}\ge \varpi$, $∥{\tilde{\mathit{e}}}_x\left(\tau |t_{l+1}\right)∥\le \frac{t_p\times \varpi }{\tau -t_{l+1}}$ is satisfied over $\tau \in [t_k+t_p,t_{l+1}+t_p)$. This completes the proof of constraint satisfaction of (34f). The aforementioned mathematical analysis conclusively proves that the designed candidate solution (41) is a feasible solution. This establishes that Problem 1 is feasible for the entire time and concludes the feasibility proof. □

Theorem 2.

Supposing that all the conditions in Theorem 1 meet, then the aircraft will satisfy constraints in the closed loop control operation.

Proof. 

For the entire time interval $t\ge 0$, the constraints satisfied in Theorem 1, then the actual system will ensure constraint satisfaction in the close loop scenario despite the presence of disturbances. □

The following theorem provides the stability analysis of the proposed technique.

Theorem 3.

Suppose the drone is controlled through (37) and by Algorithm 1 with all conditions holding in Theorems 1 and 2, then the tracking error system is ISS if

$$\underline{q}{\vartheta }^2>\frac{1}{2}\gimel e^{\upsilon t_p}(\varpi +\vartheta )+\frac{{\overline{q}}^2{\gimel }^2{t}_i}{2\upsilon }(e^{2\upsilon t_p}-e^{2\upsilon t_i})+\frac{2{\overline{q}}^2\gimel \varpi }{\sqrt{2}\upsilon }{(\frac{t_p^2}{t_i}-t_p)}^{\frac{1}{2}}{(e^{2\upsilon t_p}-e^{2\upsilon t_i})}^{\frac{1}{2}},$$

(52)

where$\underline{q}=min\left\{q_k\right\}$and$\overline{q}=max\left\{q_k\right\}$with$k=\{1,2,3,\cdots ,12\}$.

Proof. 

Choose a Lyapunov function

$$V\left({\mathit{e}}_x\left(t_l\right)\right)=J\left({\tilde{\mathit{e}}}_x^{*}\left(t_l\right),{\tilde{\mathit{e}}}_u^{*}\left(t_l\right)\right).$$

(53)

Through the Riemann integral principle, a constant $r_1$ exists such that the subsequent expression satisfies $V\left({\mathit{e}}_x\right)\ge r_1{\mathcal{S}}_c\left({\mathit{e}}_x,{\mathit{e}}_u\right)\triangleq {\beth }_1\left(∥{\mathit{e}}_x∥\right)$. From (27), $V\left({\mathit{e}}_x\left(t_l\right)\right)\le {\mathcal{T}}_p\left({\mathit{e}}_x\left(t_l\right)\right)+{\mathcal{T}}_p\left({\mathit{e}}_x\left(t_l+t_p\mid t_l\right)\right),\forall {\mathit{e}}_x\in {\mathbb{T}}_{\vartheta }$. Due to ${\mathcal{T}}_p(.)$ in the terminal region with respect to time. Hence, $V\left({\mathit{e}}_x\left(t_l\right)\right)\le 2g\left({\mathit{e}}_x\left(t_l\right)\right),\forall {\mathit{e}}_x\in {\mathbb{T}}_{\vartheta }$. As the origin lies in the ${\mathbb{T}}_{\vartheta }$, & $2g\left({\mathit{e}}_x\left(t_l\right)\right)\le \vartheta$$\forall {\mathit{e}}_x\in {\mathbb{T}}_{\vartheta }$, it satisfies that $2{\mathcal{T}}_p\left({\mathit{e}}_x\left(t_l\right)\right)\ge {\mathbb{T}}_{\vartheta }$. Because of the feasibility of Problem 1, an upper bound $r_2>\vartheta$ for $V\left({\mathit{e}}_x\left(t_l\right)\right)$ exists. Thus, ${\beth }_2\left(∥{\mathit{e}}_x\left(t_l\right)∥\right)=\frac{r_2}{\vartheta }{\mathcal{T}}_p\left({\mathit{e}}_x\left(t_l\right)\right)$ is a ${\mathcal{K}}_{\infty }$ function such that ${\beth }_2\left(∥{\mathit{e}}_x\left(t_l\right)∥\right)\ge r_2$, thereby satisfying ${\beth }_2\left(∥{\mathit{e}}_x\left(t_l\right)∥\right)\ge V\left({\mathit{e}}_x\left(t_l\right)\right)$. This confirms the existence of ${\mathcal{K}}_{\infty }$ functions ${\beth }_1$, ${\beth }_2$ satisfying (31).

Lyapunov function at $t_{l+1}$ with the actual system

$$V\left({\mathit{e}}_x\left(t_{l+1}\right)\right)=J\left({\tilde{\mathit{e}}}_x\left(t_{l+1}\right),{\tilde{\mathit{e}}}_u\left(t_{l+1}\right)\right),$$

(54)

so

$$\begin{array}{cc}\hfill \Delta V& =V\left({\mathit{e}}_x\left(t_{l+1}\right)\right)-V\left({\mathit{e}}_x\left(t_l\right)\right)\hfill \\ & \le J\left({\tilde{\mathit{e}}}_x\left(t_{l+1}\right),{\tilde{\mathit{e}}}_u\left(t_{l+1}\right)\right)-J\left({\tilde{\mathit{e}}}_x^{*}\left(t_l\right),{\tilde{\mathit{e}}}_u^{*}\left(t_l\right)\right)\hfill \\ & \triangleq \Delta V_1+\Delta V_2+\Delta V_3,\hfill \end{array}$$

(55)

where

$$\begin{array}{cc}\hfill \Delta V_1& ={\int }_{t_{l+1}}^{t_l+t_p}\left({∥{\tilde{\mathit{e}}}_x\left(\tau |t_{l+1}\right)∥}_Q^2-{∥{\tilde{\mathit{e}}}_x^{*}\left(\tau |t_l\right)∥}_Q^2\right)d\tau ,\hfill \\ \hfill \Delta V_2& ={\int }_{t_l+t_p}^{t_{l+1}+t_p}\left({∥{\tilde{\mathit{e}}}_x\left(\tau |t_{l+1}\right)∥}_Q^2+{∥{\tilde{\mathit{e}}}_u\left(\tau |t_{l+1}\right)∥}_P^2\right)d\tau \hfill \\ & +{∥{\tilde{\mathit{e}}}_x\left(t_{l+1}+t_p|t_{l+1}\right)∥}_R^2-{∥{\tilde{\mathit{e}}}_x^{*}\left(t_l+t_p|t_l\right)∥}_R^2,\hfill \\ \hfill \Delta V_3& =-{\int }_{t_l}^{t_{l+1}}\left({∥{\tilde{\mathit{e}}}_x^{*}\left(\tau |t_l\right)∥}_Q^2+{∥{\tilde{\mathit{e}}}_u^{*}\left(\tau |t_l\right)∥}_P^2\right)d\tau .\hfill \end{array}$$

Let $\Delta V_1$

$$\begin{array}{cc}\hfill \Delta V_1& \le {\int }_{t_{l+1}}^{t_l+t_p}\left({∥{\tilde{\mathit{e}}}_x\left(\tau |t_{l+1}\right)-{\tilde{\mathit{e}}}_x^{*}\left(\tau |t_l\right)∥}_Q\right)\left({∥{\tilde{\mathit{e}}}_x\left(\tau |t_{l+1}\right)∥}_Q+{∥{\tilde{\mathit{e}}}_x^{*}\left(\tau |t_l\right)∥}_Q\right)d\tau \hfill \\ \hfill & \le {\int }_{t_{l+1}}^{t_l+t_p}\left[{\overline{q}}^2\gimel t_i{e}^{\upsilon \left(\tau +t_i-t_{l+1}\right)}\times \left(2∥{\tilde{\mathit{e}}}_x^{*}\left(\tau |t_l\right)∥+\gimel t_i{e}^{\upsilon \left(\tau +t_i-t_{l+1}\right)}\right)\right]d\tau \hfill \\ \hfill & ={\int }_{t_{l+1}}^{t_l+t_p}\left[2{\overline{q}}^2{e}^{\upsilon \left(\tau +t_i-t_{l+1}\right)}∥{\tilde{\mathit{e}}}_x^{*}\left(\tau |t_l\right)∥+{\overline{q}}^2{e}^{2\upsilon \left(\tau +t_i-t_{l+1}\right)}\right]d\tau \hfill \\ & \le {\int }_{t_{l+1}}^{t_l+t_p}2\overline{q}{e}^{\upsilon \left(\tau +t_i-t_{l+1}\right)}∥{\tilde{\mathit{e}}}_x^{*}\left(\tau |t_l\right)∥d\tau +\frac{{\overline{q}}^2}{2\upsilon }\left(e^{2\upsilon t_p}-e^{2\upsilon t_i}\right).\hfill \end{array}$$

(56)

Applying Hölder inequality to the first term,

$$\Delta V_1\le {\left({\int }_{t_{l+1}}^{t_l+t_p}{∥{\tilde{\mathit{e}}}_x^{*}\left(\tau |t_l\right)∥}^{2}d\tau \right)}^{\frac{1}{2}}\frac{2{\overline{q}}^2\gimel t_i}{\sqrt{2}\upsilon }{(e^{2\upsilon t_p}-e^{2\upsilon t_i})}^{\frac{1}{2}}+\frac{{\overline{q}}^2{\gimel }^2{t}_i^2}{2\upsilon }\left(e^{2\upsilon t_p}-e^{2\upsilon t_i}\right).$$

(57)

Suppose $\Delta V_2$

$$\begin{array}{cc}\hfill \Delta V_2=& {\int }_{t_l+T}^{t_{l+1}+t_p}{∥{\tilde{\mathit{e}}}_x\left(\tau |t_{l+1}\right)∥}_Q^2+{∥{\tilde{\mathit{e}}}_u\left(\tau |t_{l+1}\right)∥}_P^{2}d\tau +{∥{\tilde{\mathit{e}}}_x\left(t_{l+1}+t_p|t_{l+1}\right)∥}_R^2\hfill \\ \hfill & -{∥{\tilde{\mathit{e}}}_x^{*}\left(t_l+t_p|t_l\right)∥}_R^2+{∥{\tilde{\mathit{e}}}_x\left(t_l+t_p|t_{l+1}\right)∥}_R^2-{∥{\tilde{\mathit{e}}}_x\left(t_l+t_p|t_{l+1}\right)∥}_R^2.\hfill \end{array}$$

(58)

Integrating (27) from $t_l+t_p$ into $t_{l+1}+t_p$ and substituting it into (58)

$$\begin{array}{cc}\hfill \Delta V_2\le & {∥{\tilde{\mathit{e}}}_x\left(t_l+t_p|t_{l+1}\right)∥}_R^2-{∥{\tilde{\mathit{e}}}_x^{*}\left(t_l+t_p|t_l\right)∥}_R^2\hfill \\ \hfill \le & \left(\frac{1}{2}∥{\tilde{\mathit{e}}}_x\left(t_l+t_p|t_{l+1}\right)-{\tilde{\mathit{e}}}_x^{*}\left(t_l+t_p|t_l\right)∥\right)\hfill \\ & \times \left(∥{\tilde{\mathit{e}}}_x\left(t_l+t_p|t_{l+1}\right)∥+∥{\tilde{\mathit{e}}}_x^{*}\left(t_l+t_p|t_l\right)∥\right)\hfill \\ \hfill \le & \frac{1}{2}\daleth \beth e^{\upsilon t_p}(\vartheta +\varpi ).\hfill \end{array}$$

(59)

Consider $\Delta V_3$

$$\begin{array}{cc}\hfill \Delta V_3& <-{\int }_{t_l}^{t_{l+1}}{∥{\tilde{\mathit{e}}}_x^{*}\left(\tau |t_l\right)∥}_Q^{2}d\tau ,\hfill \\ \hfill & \le -\underline{q}{t}_i{\vartheta }^2.\hfill \end{array}$$

(60)

By the addition of (57), (59) and (60), $\Delta V\triangleq \Delta V_1+\Delta V_2+\Delta V_3$ satisfies

$$\begin{array}{cc}\hfill \Delta V<& -\underline{q}{t}_i{\vartheta }^2+\frac{1}{2}\gimel t_i{e}^{\upsilon t_p}(\vartheta +\varpi )+\frac{2{\overline{q}}^2\gimel t_i}{\sqrt{2}\upsilon }{\left(e^{2\upsilon t_p}-e^{2\upsilon t_i}\right)}^{\frac{1}{2}}\hfill \\ & +\frac{{\overline{q}}^2{\gimel }^2{t}_i^2}{2\upsilon }\left(e^{2\upsilon t_p}-e^{\upsilon t_i}\right).\hfill \end{array}$$

(61)

From the condition (52), $\Delta V<0$ satisfies. When ${\mathit{e}}_x\in {\mathbb{T}}_{\vartheta }$, re-assume $\Delta V_1$ and $\Delta V_3$:

$$\begin{array}{cc}\hfill \Delta V_1& \le {\int }_{t_{l+1}}^{t_l+t_p}2{\overline{q}}^2\gimel t_i\vartheta (e^{\upsilon \tau -t_{l+2}}d\tau +\frac{{\overline{q}}^2{\gimel }^2{t}_i^2}{2\upsilon }\left(e^{2\upsilon t_p}-e^{2\upsilon t_i}\right)\hfill \\ \hfill & =\frac{2{\overline{q}}^2\beth t_i\vartheta }{\upsilon }(e^{\upsilon t_p}-e^{\upsilon t_i})+\frac{\overline{q}{\overline{q}}^2{\beth }^2{t}_i^2}{2\upsilon }(e^{\upsilon t_p}-e^{\upsilon t_i}).\hfill \end{array}$$

(62)

Because of the decreasing property of ${∥{\tilde{\mathit{e}}}_x\left(\tau |t_l\right)∥}_Q^2$ in ${\mathbb{T}}_{\vartheta }$, then $\Delta V_3$ follows

$$\Delta V_3\le -\underline{q}{t}_i{∥{\tilde{\mathit{e}}}_x^{*}\left(t_{l+1}|t_l\right)∥}^2.$$

(63)

From (48), we obtain ${∥{\mathit{e}}_x\left(t_{l+1}\right)∥}^2\le {∥{\mathit{e}}_x^{*}\left(t_{l+1}|t_l\right)∥}^2+{\gimel }^2{t}_i^2{e}^{2\upsilon t_i}+2\vartheta \gimel t_i^2{e}^{\vartheta t_i}$. As a result,

$$\Delta V_3\le -\underline{q}{t}_i{∥{\tilde{\mathit{e}}}_x\left(t_{l+1}\right)∥}^2\underline{q}{\gimel }^2{t}_i^3{e}^{2\vartheta \gimel }+2\underline{q}\vartheta \gimel t_i^2{e}^{\vartheta t_i}.$$

Consequently, it satisfies that

$$\Delta V\le -\underline{q}{t}_i{∥{\tilde{\mathit{e}}}_x\left(t_{l+1}\right)∥}^2+\daleth (\gimel )$$

(64)

where

$\daleth (\gimel )=\frac{2{\overline{q}}^2\beth t_i\vartheta }{\upsilon }(e^{\upsilon t_p}-e^{\upsilon t_i})+\frac{{\overline{q}}^2{\beth }^2{t}_i^2}{2\upsilon }(e^{\upsilon t_p}-e^{\upsilon t_i})+\frac{1}{2}\beth \daleth e^{\upsilon t_p}(\vartheta +\varpi )+\underline{q}{\beth }^2{t}_i^3{e}^{2\upsilon t_i}+2\underline{q}\vartheta \beth t_i^2{e}^{\upsilon t_i}$ is a $\mathcal{K}$ function w.r.t. ℶ. Hence, the stability is comprehensively proved. □

4. Comparative Analysis

4.1. Simulation Studies

The parameters of the aerial robot are provided as follows: $m=1.6$ kg, $I_{xx}=I_{yy}=0.025$ kg.m${}^2$, $I_{zz}=0.005$ kg.m${}^2$, $l_1=0.17$ m, $l_2=0.06$ m, $g=9.8$ m/s${}^2$. 30% parametric uncertainty for the entire flight is considered during the simulation study. To validate the reliability of the controller, a high-speed $\upsilon =20$ m/s is selected for the comparative analysis. In the cost function, $\mathit{R}=\mathrm{diag}\left\{0.10.10.10.111\right\}$, and the values of $q_i$ and $p_i$ are 1 and 0.25, respectively. In Problem 1, the numerical values of $t_p$, $t_i$, $\vartheta$, ${\mathcal{F}}_i$, $\varpi$ and $\varkappa$ are chosen as 10 s, 0.05 s, 0.04, 4, 0.08 and 0.2, respectively. $\lambda$ is selected as 0.02, 0.03, 0.04, 0.042, 0.05, 0.051, 0.06, 0.07, 0.08, 0.085, 0.09, and 0.095. The constraints are $-1\le {\mathit{\xi }}^{\mathcal{I}}\le 1$, $-2.5\le {\mathit{V}}_{\mathcal{I}}\le 2.5$, $0\le T\le 4$, $-\pi /4\le \mathsf{\Theta }\le \pi /4$, and $-2\le {\mathit{\omega }}_{\mathcal{B}}\le 2$.

In order to corroborate the control performance, different trajectories are defined that show that the DFAV has the properties of fixed-wing aircraft and helicopters. Suppose the drone is flying with an angular velocity of 2 rad/s. In the first trajectory (65), the aircraft takes off, and after attaining a certain altitude, it performs a bank-turn motion:

$$\begin{array}{cc}\hfill {\xi }_x& =sin\left(2t\right),\hfill \\ \hfill {\xi }_y& =5cos\left(2t\right),\hfill \\ \hfill {\xi }_z& =50-e^{-3t}.\hfill \end{array}$$

(65)

In the second instance, we validate the control performance with a helical 3D trajectory:

$$\begin{array}{cc}\hfill {\xi }_x& =sin\left(2t\right),\hfill \\ \hfill {\xi }_y& =5cos\left(2t\right),\hfill \\ \hfill {\xi }_z& =50t-cos(-2t).\hfill \end{array}$$

(66)

Moreover, since the aerial robot is flying with high speed $\upsilon$, we validate the performance with a straight and level flight (SLF), a kind of horizontal flight phase of DFAVs [14]. During this flight, the drone must keep a constant altitude (e.g., 40 m) during the entire time, while the aircraft should fly at the desired angles, i.e., $\alpha =$5${}^{\circ }$ and $\theta =$5${}^{\circ }$. In this case, some initial conditions are required to be considered, e.g., ${\xi }_x\left(0\right)=0$ m, ${\xi }_z\left(0\right)=39$ m, $\alpha \left(0\right)=4.8^{\circ }$, and $\theta \left(0\right)=4.8^{\circ }$.

Three different cases are considered:

• Case (1): Performance in the presence of parametric uncertainties;
• Case (2): Performance in the presence of parametric uncertainties and disturbances;
• Case (3): Performance in the presence of faults, model uncertainties and disturbances.

Before proceeding to the aforementioned case studies, the correlation between flight performance and endurable wind speed is provided. In the previous work based on this aerial vehicle [47], aircraft was flown at various speeds under different weather conditions. From these real-time experiments, it was found that the maximum endurable wind speed is approximately 8 m/s. Therefore, it can be concluded that effective flight performance can be achieved if the wind speed is within the maximum endurable wind speed limit. However, the flight performance degrades as the wind speed passes the maximum endurable wind speed threshold. If the aircraft keeps flying over this maximum endurable wind speed, it encounters a structural failure, or a severe fault may occur. Nonetheless, the severity of this structural failure or fault depends on the magnitude and the total duration of the wind gust it encounters. To keep in view the contributions and main objective of this article, we deal with the control performance only under three case studies, under calm weather (Case 1), and wind gust is present (Case 2). However, it is assumed that it may not cause any structural faults because the magnitude of the wind gust is mostly under 8 m/s during the entire flight. Lastly, it is considered in the third case study that the partial fault exists associated with the unavailability of the dynamical information of the wings whose probability of occurrence is less than 3%. This occurrence level is uncommon, but to precisely verify the efficacy of the control approach, the worst-case situation (less than 3%) shall be analyzed in the third case study.

4.1.1. Case(1): Trajectory Tracking Response under Model Uncertainties

A comparative analysis of both strategies under model uncertainties is shown in Figure 5. In Figure 5a, the simulation results based on 3D trajectory (65) and (66) are illustrated. Furthermore, the tracking response under model inaccuracies derived from SLF flight is depicted in Figure 5b. Both controllers are primarily constructed on the principle of nominal MPC with similar attributes, so they exhibit adequate tracking performance. Nevertheless, the proposed technique demonstrates effective control performance with fast response and rapid convergence time.

4.1.2. Case(2): Trajectory Tracking Response under Uncertainties and Disturbances

Recall disturbance model in (29): the continuous Dryden wind turbulence model MIL-HDBK-1797B is employed using Simulink with default parameters chosen. Due to some minor limitations in this turbulence model, wind components are added, depicted in Figure 6. These wind components are composed of sinusoidal waveform ${\Delta }_{\mathrm{wind}}={[sin\left(8t\right)sin\left(6.01t\right)sin\left(4.02t\right)]}^{\mathrm{T}}$. The trajectory tracking performance under model uncertainties and disturbances are shown in Figure 7. Note that the MPC-CFC performance, as opposed to the developed approach, deteriorates despite having a dedicated DOB mechanism. This conclusively indicates that DOB in MPC-CFC [16,17] is only suitable for slow time-varying disturbances.

4.1.3. Case(3): Tracking Response in the Presence of Model Uncertainties, Disturbances, and Fault

In the third case study, suppose a fault occurs during the initial few seconds of the simulations due to the presence of a strong wind gust, considered in (29), and this fault is considered as a partial malfunction in one of the wings of the aircraft, whose role is depicted in Figure 2, is now either unavailable or unknown in this case. The tracking performance subject to disturbances, uncertainties, and fault is shown in Figure 8. More specifically, this fault has an adverse effect while controlling the attitude of the aerial platform illustrated in Figure 8b. In this aspect, as adaptive SINDY is equipped to handle and compensate for model divergence and correction, unlike DOB developed offline in [16,17], the presented approach by exploiting the full potential of MPC reveals effective control performance compared to the MPC-CFC technique [16,17].

4.1.4. Performance Analysis Based on Tracking Error and Computations

To analyze the effectiveness of both techniques, mean absolute error (MAE) is utilized and is defined as follows:

$$\mathrm{MAE}=\sum _{j=1}^{n_{\mathrm{obs}}}\left|e\left(t_j\right)\right|/n_{\mathrm{obs}},$$

(67)

where $e\left(t_j\right)$ and $n_{\mathrm{obs}}$ are the tracking error in each time step and the sum of a number of observations, respectively. In

Table 2. Performance analysis based on tracking error (MAE).

	Scenario	MPC-CFC	Proposed		Scenario	MPC-CFC	Proposed
${\xi }_y$(m)	Case(1)-Figure 5a	1.5752	0.0016	${\xi }_z$(m)	Case(1)-Figure 5b	0.1711	0.1424
	Case(1)-Figure 5a	0.1569	0.0157		Case(1)-Figure 5a	0.0011	1.09$\times {10}^{-6}$
	Case(2)-Figure 7a	2.3279	0.0141		Case(2)-Figure 7b	0.0999	0.0987
	Case(2)-Figure 7a	0.2357	0.0235		Case(2)-Figure 7a	0.0016	9.82$\times {10}^{-6}$
	Case(3)-Figure 8a	2.8687	0.0157		Case(3)-Figure 8b	0.122	0.0782
	Case(3)-Figure 8a	0.4729	0.0314		Case(3)-Figure 8a	0.002	1.091$\times {10}^{-5}$
${\xi }_x$(m)	Case(1)-Figure 5a	0.03066	0.0003	$\theta$(${}^{°}$)	Case(1)-Figure 5b	0.0482	0.0378
	Case(1)-Figure 5a	0.0325	0.0033		Case(2)-Figure 7b	0.0376	0.0349
	Case(2)-Figure 7a	0.4429	0.0029		Case(3)-Figure 8b	0.0748	0.0284
	Case(2)-Figure 7a	0.0486	0.049	$\alpha$(${}^{°}$)	Case(1)-Figure 5b	0.0466	0.0460
	Case(3)-Figure 8a	0.5446	0.0033		Case(2)-Figure 7b	0.0371	0.0344
	Case(3)-Figure 8a	0.0966	0.0065		Case(3)-Figure 8b	0.069	0.0389

, a comparative analysis is illustrated. The trajectory error suggests that the designed approach shows effective tracking performance.

The numerical simulations are executed on Windows-driven Intel-based dual-core processor (2.2 GHz) with 8 GB of RAM. The optimization Problem 1 is written in MATLAB 2021a via CasADi [48], 3.5.5 version, and computed by NLP solver IPOPT [49], 3.12.3 version. The computational cost of these schemes is determined by the time required to compute the control action during the trajectory tracking response, depicted in Figure 9. Based on these findings, the average normalized computational time for the entire simulation study for the developed algorithm and MPC-CFC strategy is 18.62 ms and 19.51 ms, respectively.

4.2. Discussion

Different case studies have shown that the developed algorithm demonstrates an efficient tracking performance compared to the MPC-CFC technique. It is primarily due to the hybrid modeling approach integrated with the control mechanism that enables the MPC-based controller to perform effectively despite uncertainties, disturbances, and faults. This also confirms our earlier assertion that physics-informed dynamics and DOB in the disturbance rejection MPC approaches (e.g., [16,17]) limited the capability of MPC due to inadequate modeling and ineffective estimation. Nevertheless, the computational complexity may be affected if a large amount of data are utilized online rather than with less data. With regard to the ML technique related to adaptive SINDy, this work has not addressed any shortcomings concerning this method (e.g., noise sensitivity and long-term memory), as they are beyond the scope of this article. Regarding flight robustness, the proposed framework is designed based on a more suitable selection of the terminal region, which yields adequate stability margins as opposed to [16,17]. In spite of all the aforementioned benefits, the presented MPC-based technique using a physics-informed ML model still ensures sub-optimal performance.

5. Conclusions and Future Directions

This paper has presented an MPC-based control design for the ducted fan aircraft utilizing a hybrid modeling scheme known as physics-informed ML. In the beginning, the physics-informed model was derived offline from the drone with sufficient capabilities. Thereafter, an online-based data-driven modeling technique is integrated with physics-informed dynamics to determine the actual model. Afterwards, an MPC-based control algorithm was developed by updating the physics-informal dynamics (nominal system) with real states. ISS stability and recursive feasibility were proven under adequate conditions. Finally, simulations were conducted under three different scenarios starting from challenging (case 1), worse (case 2), and worst-case situation (case 3), which revealed that the developed framework exhibits effective control performance.

For future studies, the worst-case scenario can be utilized as a baseline to effectively construct a more refined fault-tolerant controller. Nevertheless, it is noteworthy that the developed framework is not designed as a fault-tolerant controller. Moreover, the worst-case scenario is only verified under a partial fault occurrence. Furthermore, the designed strategy also ensures sub-optimal performance like [16,17] for aerial applications.

In the end, several future directions can be pursued to improve the efficacy of the proposed strategy:

• Physics-informed modeling: the designed model is basically derived from the principle of the Newton–Euler method. Other mathematical models, such as the Lagrange-Euler approach, can be employed;
• Data-driven ML: the developed ML scheme was inspired by adaptive SINDy [30], where different ways of improving the ML model’s capability can be explored. Moreover, further investigation is required to find an effective technique to enhance computational efficiency with more data and less physics;
• Real-time implementation: future work will involve real-time testing of the presented framework. For this purpose, a more suitable solver can be used for code generation especially developed for real-time embedded optimization.