# Nonlinear model predictive control trajectory tracking for a fixed-wing UAV in a deep-stall and perching landing maneuver with guaranteed stability

Huu Thien Nguyen, Ionela Prodan, Fernando a C C Fontes

► To cite this version:

Huu Thien Nguyen, Ionela Prodan, Fernando a C C Fontes. Nonlinear model predictive control trajectory tracking for a fixed-wing UAV in a deep-stall and perching landing maneuver with guaranteed stability. European Journal of Control, 2026, 88, 10.1016/j.ejcon.2025.101442. hal-05497443

HAL Id: hal-05497443

https://hal.science/hal-05497443v1

Submitted on 6 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Nonlinear model predictive control trajectory tracking for a fixed-wing UAV in a deep-stall and perching landing maneuver with guaranteed stability

Huu Thien Nguyen<sup>a,*</sup>, Ionela Prodan<sup>a, b,1</sup>, Fernando A.C.C. Fontes<sup>a</sup>

## ARTICLE INFO

2000 MSC:

49N90

65K10

93D20

Keywords:

Deep-stall

Perching

Nonlinear model predictive control (NMPC)

Trajectory tracking

Trajectory optimization

Fixed-wing UAVs

## ABSTRACT

We address the planning and control problems of a fixed-wing unmanned aerial vehicle (FW UAV) in a deep-stall and perching landing maneuver. We use optimal control methods to generate a trajectory that satisfies not only the aircraft maneuver requirement but also the recursive feasibility demands of the tracking control method that we use subsequently. To track the generated trajectory, we developed a Nonlinear Model Predictive Control (NMPC) scheme. We then outline the steps to systematically derive the NMPC terminal costs and the terminal sets tailored to this problem and formally prove the stability of the NMPC tracking scheme. Simulation results for different initial conditions show the viability of the developed scheme for the deep-stall and perching landing maneuvers. Finally, we examine the relationships between airspeed, wing loading, elevator deflection angle, and flight path angle in steady flight, including deep-stall landing.

## 1. Introduction

Aerobatic flights are intentional maneuvers that incorporate an abrupt change in an aircraft's attitude, an abnormal attitude, or abnormal acceleration, which are not necessary for normal flight (FAA, 2024). Aerobatic flights are always intriguing since they show the capabilities of aircraft that are rarely demonstrated during their nominal operations, as observed in commercial flights or popular fixed-wing unmanned aerial vehicle (FW UAV) applications like package delivery, surveillance, mapping, etc. One category of aerobatic flights is *post-stall maneuvers*, executed when the aircraft is in a stall.

A stall happens when the FW UAV's angle of attack (AoA) exceeds the critical AoA—the value defined when the lift coefficient is maximum—causing the airflow around its wings to detach, leading to sudden loss of lift force and reduced controllability. However, during this time, the drag force increases, and by exploiting the stall region for FW UAVs, fascinating post-stall maneuvers can be executed, namely deep-stall landing (Mathisen et al., 2021), perching landing (Wüest et al., 2024), and post-stall turning (Basescu & Moore, 2020). Real-life applications of post-stall maneuvers include the Raven B UAV (AeroVironment, 2021) that belly lands in a deep-stall and breaks apart upon touching the

ground to dissipate the energy or a cobra maneuver that is executed in fighter jets (Phillips, 2004, Chapter 3). In this paper, we only focus on deep-stall landing, perching landing, and a combination of them.

A deep-stall occurs when the aircraft is in a stall and its AoA exceeds far beyond its critical AoA. Subsequently, as the airflow surrounding the wings stabilizes, the aircraft dives stably into the post-stall region (Taniguchi, 2008). During deep-stall, the large AoA leaves the aircraft with large altitude losses and limited control power (White & Cooper, 1965). From a different perspective, perching is a bio-inspired maneuver that also makes use of airflow separation and the large drag force, but at a higher AoA, which could be more than 90°, to land an aircraft at a sufficiently small airspeed (AliKhan et al., 2013; Wickenheiser & Garcia, 2008).

Optimization-based techniques have been employed to generate a deep-stall landing and a perching landing trajectory by solving an optimization problem subject to a set of input and state constraints. Concerning *perching*, in Wickenheiser and Garcia (2008), an optimization problem with two cost functions is solved for a morphing UAV: minimization of the starting distance, and the altitude overshoot. The final velocity of the UAV is constrained to be 5% of the initial velocity. Roberts et al. (2009) generate a perching trajectory by solving an open-

<footnote><sup>a</sup> SYSTEC-ISR ARISE, Faculdade de Engenharia, Universidade do Porto, Rua Dr. Roberto Frias, Porto, 4200-465, Portugal</footnote>

<footnote><sup>b</sup> Univ. Grenoble Alpes, Grenoble INP, LCIS, Valence, F-26000, France</footnote>

<footnote>* Corresponding author.
E-mail addresses: huu-thien.nguyen@ieee.org (H.T. Nguyen), ionela.prodan@lcis.grenoble-inp.fr (I. Prodan), faf@fe.up.pt (F.A.C.C. Fontes).</footnote>

<footnote n="1">¹ Institute of Engineering Univ. Grenoble Alpes</footnote>
loop optimal control problem (OCP) with nominal quadratic state cost, input cost, and a terminal state cost. In AliKhan et al. (2013), an OCP minimizing the total distance (by formulating the cost as an integration of airspeed over time) while constraining the final states to get the perching landing trajectory for a fixed-wing aircraft in three dimensions is considered. The terminal airspeed is constrained to be less than or equal to 30% of the initial airspeed. Venkateswara Rao et al. (2015) also propose an optimization problem to minimize the distance traveled while constraining the final airspeed to be less than 5% of the initial one, but only with the longitudinal dynamics. Waldock et al. (2018) compare various optimization-based perching trajectories with reinforcement learning-based perching trajectories. Feroskhan et al. (2020) constrain the final airspeed to be equal to 0.05 times the initial airspeed while also minimizing the total distance traveled in a 3D perching to generate the reference perching trajectory. Recently, Wüest et al. (2024) studies the swooping (diving followed by perching) maneuvers in birds and morphing drones by solving an OCP that minimizes the final kinetic energy and the distance traveled at a high AoA while constraining the vertical displacement to one meter downward from the initial height. Regarding optimization in deep-stall landing, Mathisen et al. (2021) utilize NMPC to precisely deep-stall land an FW UAV. Basescu et al. (2023) employ NMPC with learned aerodynamics to dive and stall their FW UAVs.

There exist other works dealing with the post-stall regime that do not focus on optimization. Kim and Park (2022) guide their FW UAV in the approaching phase and deep-stall landing phase by tracking the computer mouse with live aerial images from an onboard camera. Two PI controllers are used to track the reference heading and flight path angles. Mehanovic et al. (2017) design a FW UAV with special landing gear that they call legs, which help their drone perch on a vertical wall.

After generating an optimal post-stall trajectory by solving a constrained optimization problem, several control techniques with stability assured have been used to track it. Lyapunov-based control techniques have been used predominantly to guarantee stability for the **perching** maneuver. Feroskhan et al. (2020) design backstepping controllers and Feroskhan and Go (2018) propose an adaptive sliding mode controller to track a 3D perching trajectory. Sliding mode control is also used in Venkateswara Rao and Go (2014) to track a longitudinal perching trajectory. The main drawback of these tracking techniques is that the input and state constraints are only certified from the trajectory generation step and are not built into the controllers themselves. Therefore, large initial variations in the initial states could lead to impracticality in controlling the UAVs. Moreover, finding the Lyapunov functions associated with these techniques is also a challenging task. Time-varying linear quadratic regulator (LQR) controllers are employed in Moore et al. (2014) to track a longitudinal perching trajectory of a fixed-wing glider and in Basescu and Moore (2020) to track a post-stall turning trajectory of an aerobatic aircraft. However, LQR controllers cannot take into account the input constraints of the UAV. Model Predictive Control appears to be a suitable candidate for trajectory tracking problems with its ability to handle input and state constraints (Mayne, 2001). Systems with time-delays are not considered in this research (see e.g., Aslam et al. (2025, 2024), where time-delays are explicitly considered in the models). Concerning deep-stall, Mathisen et al. (2021) and Basescu et al. (2023) do not consider the stability perspective of their NMPC controllers.

Various works in the literature have solved the problem of **guaranteed stability in trajectory tracking with an NMPC controller**. Faulwasser and Findeisen (2011) design an NMPC scheme for an asymptotically constant trajectory. The terminal cost comprises the cost from the initial time to the boundary of the terminal set, and the second part is inside the terminal set with a local controller. Their terminal sets are time-varying level sets of Lyapunov functions. Alessandretti et al. (2013) propose an NMPC tracking scheme with asymptotic arbitrarily small tracking error. The size of their terminal set is dependent on the constraints of the systems. Köhler et al. (2020a) augment the results of the quasi-infinite

MPC in Chen and Allgöwer (1998) for trajectory tracking with exponential stability using time-varying terminal weighting matrices and time-varying terminal sets.

This paper extends the results presented in Nguyen et al. (2023), where we proposed a recovery technique for FW UAVs in a maneuver comprising a deep-stall and a perch on a recovery net. However, in that work, the stability aspect of the NMPC tracking scheme in a deep-stall with a perch landing maneuver was not considered. In this work, we address the stability of the NMPC tracking task, referring to this problem as *stable deep-stall and perching landing*. The contributions of this paper, compared to existing works in the literature, are as follows:

* We generate an optimal trajectory that combines a deep-stall segment with a perching segment and includes a stability phase at the end of the trajectory for a hybrid-axis longitudinal fixed-wing UAV dynamics by solving optimal control problems.

* We employ an NMPC tracking framework with time-varying terminal ingredients, namely terminal costs and terminal sets, to asymptotically track a deep-stall and perching landing trajectory for an FW UAV while ensuring full satisfaction of input and state constraints.

* We present the stability properties of the NMPC tracking scheme with detailed proofs, then show their implementation by simulating the deep-stall and perching maneuver from various initial positions.

* We analyze the relations between wing loading, airspeed, elevator deflection angle, and flight path angle in a steady flight, to validate the effects of wing loading on airspeed and flight path angle in a steady deep-stall.

**Notations:** Let $\mathbb{l}_n$ be an identity matrix of size $n$, $0_{m \times n}$ be a zero matrix of size $m \times n$, and $\mathbb{0}$ be a zero matrix of appropriate dimension. We use $\bar{x}$ and $\bar{x}$ to denote the upper-bound and lower-bound defined for the variable $x$, respectively. $\|x\|_p^2$ represents the weighted-P norm of a vector $x$, that is, $\|x\|_p^2 = x^\top P x$, while $\|x\|$ is the Euclidean norm. Let $S_+^n = (S_{++}^n)^\top$ be the vector space of $n \times n$ real symmetric positive semidefinite (positive definite) matrices. Let $\square$ be the reference value, while $x_t$ is the value of the time-varying function $x(\cdot)$ at the time instant $t$, i.e. $x_t = x(t)$. We denote the discrete variable at the time-step $k \in \mathbb{N}$ by the notation $x[k]$. The symbol $\oplus$ denotes the Minkowski sum, i.e. $A \oplus B = \{a+b \mid a \in A, b \in B\}$.

The outline of the paper is as follows. Section 2 describes the problem formulation and recalls the longitudinal FW UAV dynamics. Section 3 formulates the OCPs to obtain reference deep-stall and perching trajectories. Next, in Section 4, we delineate the NMPC tracking framework with stability guarantees. Section 5 shows the generated optimal deep-stall and perching maneuvers and demonstrates the stability of the NMPC tracking controller through simulations. Finally, Section 6 summarizes the paper and suggests future directions.

## 2. Prerequisites

### 2.1. Problem formulation

The problem formulation is as follows, which is also depicted in Fig. 1. The stable deep-stall and perching landing comprises two parts. In the first part, the FW UAV executes a deep-stall landing from a steady-level flight to quickly reduce its altitude, then performs a perching maneuver to perch on a recovery net that is located close to the ground. In the second part, which occurs after the recovery net, the FW UAV will recover from post-stall to steady-level flight. The end of the second part is only “imaginary” for the purposes of the NMPC tracking task. The two phases are carried out autonomously by using the results of constrained OCPs. Although the second phase is not physically realized, it is included to ensure that the prediction horizon of the NMPC tracking scheme is sufficiently long to guarantee the feasibility of the first NMPC iteration. In practice, the NMPC tracking task ends when the FW UAV touches the recovery net.
<figure><img src="image_2.png" /><figcaption>Fig. 1. Sketch of the problem.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2. Coordinate systems, forces, and moment acting on a longitudinal model fixed-wing UAV.</figcaption></figure>

## 2.2. Longitudinal fixed-wing UAV dynamics

Consider the longitudinal FW UAV model in Beard and McLain (2012), Nguyen et al. (2023). We define two reference frames: the inertial frame $I\{O^I, x^I, z^I\}$, which is fixed, pointing north and down, with origin $O^I$ on the ground, below the recovery net; the body frame $B\{O^B, x^B, z^B\}$ is attached to the center of mass of the UAV, which moves along with the UAV, as in Fig. 2.

Since this work focuses on a nominal stability problem, we disregard all the disturbances (including wind and gusts) to the UAV. The longitudinal FW UAV dynamical model in Nguyen et al. (2023, Eq. 1) is therefore unsuitable, as it does not explicitly include the state variables angle of attack $\alpha$ and airspeed $V_a$, the two most relevant variables to analyze a deep-stall and perching maneuver.

Therefore, we now employ the hybrid-axis longitudinal dynamics of an FW UAV without wind, as presented in Stengel (2004, Chapter 5)

$$ \dot{x} = V_a \cos(\gamma), \qquad (1a) $$

$$ \dot{z} = -V_a \sin(\gamma), \qquad (1b) $$

$$ m\dot{V}_a = -F_D - mg \sin(\gamma) + F_p \cos(\alpha), \qquad (1c) $$

$$ m\dot{\alpha} = m\omega_y - \frac{1}{V_a} [F_L - mg \cos(\gamma) + F_p \sin(\alpha)], \qquad (1d) $$

$$ \dot{\theta} = \omega_y, \qquad (1e) $$

$$ \dot{\omega}_y = M/I_{yy}. \qquad (1f) $$

Here, the variables x and z are the positions of the FW UAV in the inertial frame $I$; m and g are the UAV's mass and the gravitational acceleration, respectively; $V_a$ and $\alpha$ are the airspeed and the angle of attack, $\theta$ is the

pitch angle, while $\gamma$ is the flight path angle, satisfying the relation

$$ \theta = \alpha + \gamma. \qquad (2) $$

Aerodynamically, $F_L$ and $F_D$ are lift and drag forces, and M is the pitching moment around the y$^B$ axis, given by

$$ F_i = \frac{1}{2} \rho V_a^2 S \left[ C_i(\alpha) + C_{i_q} \frac{c}{2V_a} \omega_y + C_{i\delta_e} \delta_e \right], \qquad (3a) $$

$$ M = \frac{1}{2} \rho V_a^2 S c \left[ C_m(\alpha) + C_{m_q} \frac{c}{2V_a} \omega_y + C_{m\delta_e} \delta_e \right], \qquad (3b) $$

with the aerodynamic coefficients

$$ C_i(\alpha) = [1 - \sigma(\alpha)] C_i^{\text{pre}}(\alpha) + \sigma(\alpha) C_i^{\text{post}}(\alpha), \qquad (4a) $$

$$ C_m(\alpha) = C_{m0} + C_{m\alpha}, \qquad (4b) $$

here, the subscript $\square_i$, for $i \in \{L, D\}$, denotes “lift” and “drag”, $\rho$ is the air density, $S$ is the planform area, $c$ is the mean chord of the wing, $\omega_y$ is the pitch rate, $C_{i_q}$, $C_{i\delta_e}$, $C_{m_q}$, $C_{m\delta_e}$ are aerodynamic constants, and $\delta_e$ is the elevator deflection angle. In the pre-stall regime, we employ the typical linear and quadratic functions for lift and drag (Beard & McLain, 2012), while in the post-stall regime, the aerodynamics are derived from the “flat-plate model” of Cory and Tedrake (2008)

$$ C_L^{\text{pre}}(\alpha) = C_{L0} + C_{L\alpha}, \qquad C_L^{\text{post}}(\alpha) = \sin(2\alpha), \qquad (5a) $$

$$ C_D^{\text{pre}}(\alpha) = C_{Dp} + \frac{[C_L^{\text{pre}}(\alpha)]^2}{\pi e_O A R}, \qquad C_D^{\text{post}}(\alpha) = 2 \sin^2(\alpha). \qquad (5b) $$

In (5), $C_{L0}$ and $C_{L\alpha}$ are aerodynamic constants, $C_{Dp}$ is the parasitic drag, $e_O$ is the Oswald efficiency factor, $AR$ is the wing aspect ratio, the AoA $\alpha$ has the unit $[\text{rad}]$, and the $\sigma(\alpha)$ in (4a) is a sigmoid function used as a blending function between the two regimes

$$ \sigma(\alpha) = \frac{1}{1 + e^{-M(\alpha - \alpha_0)}}, \qquad (6) $$

where $\tilde{M}$ is the transition rate and $\alpha_0$ is the cut-off AoA.

Finally, the propulsive force $F_p$ is

$$ F_p = \frac{1}{2} \rho S_{\text{prop}} C_{\text{prop}} (k_{\text{motor}} \delta_t)^2 - V_a^2, \qquad (7) $$

where $S_{\text{prop}}$ and $C_{\text{prop}}$ are propeller's parameters, $k_{\text{motor}}$ is the motor constant, and $\delta_t \in [0, 1]$ is the pulse-width modulation of the propeller.

The hybrid-axis longitudinal dynamical model of an FW UAV (1) is written in the shortened form as

$$ \dot{\mathbf{x}} = f_{ha}(\mathbf{x}, \mathbf{u}), \qquad (8) $$

with $\mathbf{x} \triangleq [x, z, V_a, \alpha, \theta, \omega_y]^\mathrm{T} \in \mathcal{X} \subseteq \mathbb{R}^2 \times \mathbb{R}_{>0} \times \mathbb{R}^3$ and $\mathbf{u} \triangleq [\delta_e, \delta_t]^\mathrm{T} \in \mathcal{U} \subseteq \mathbb{R} \times [0, 1].$
The number of states is $n_s = 6$ and the number of control inputs is $n_u = 2$.

We now present the definitions of a steady flight and a steady-level flight from the dynamical perspective.

**Definition 1 (Steady flight).** The FW UAV is in steady flight when its states satisfy

$$ \dot{V}_a = 0, V_a = V_{a_r} > 0, \dot{\alpha} = 0, \dot{\theta} = 0, \dot{\omega}_y = 0. \quad (9) $$

**Definition 2 (Steady-level flight).** The FW UAV is in steady-level flight when it is in steady flight with the flight path angle equal to zero (i.e., $\gamma = 0$).

### 3. Trajectory generation for the stable deep-stall and perching maneuver

We decompose the reference trajectory into two parts. The first part (i) is when the FW UAV deep-stalls and perches on the recovery net, and the second part (ii) is when the UAV recovers from the post-stall region to fly in steady-level flight after the net. Note that in real-life scenarios, the tracking will stop when the FW UAV reaches the net. The reason for the second part is threefold. First (ii-1), because we are dealing with a stability problem, the UAV must be in a steady state when time goes to infinity, and a steady-level flight is a steady state for FW UAVs. Second (ii-2), by having a recovery reference trajectory with an “imaginary” steady-level flight trajectory at the end, the prediction horizon length could be arbitrarily long to satisfy the recursive feasibility property of a stable NMPC scheme, i.e., the first NMPC tracking iteration is guaranteed to give a feasible optimal solution. Third (ii-3), if there was no recovery net, the UAV could execute another nominal landing procedure low to the ground (starting from the starting position of the “imaginary” steady-level flight), including a flare maneuver followed by a touchdown.

We generate the landing trajectory by solving two consecutive OCPs, the first is the deep-stall and perching phase (the landing phase) in the interval $[t_0, t_f]$, and the other for the post-stall to steady-level flight (the stability phase) in the interval $[t_f, t_s]$. We then concatenate the solutions for the desired trajectory.

The landing phase trajectory, which comprises the deep-stall and perching maneuvers, is generated by solving an OCP:

$$ \min_{u(\cdot)} J_1(\mathbf{x}, \mathbf{u}) = \int_{t_0}^{t_f} \ell(\mathbf{x}(t), \mathbf{u}(t))dt, \quad (10) $$

$$ \text{subject to: } \dot{x}(t) = f_{ha}(x(t), u(t)), \ t \in [t_0, t_f] \quad (11a) $$

$$ \mathbf{u}(t) \in \mathcal{V}_g, \mathbf{x}(t) \in \mathcal{X}_{g_1}, \ t \in [t_0, t_f] \quad (11b) $$

$$ \mathbf{x}(t_0) = \mathbf{x}_0, \quad (11c) $$

$$ \mathbf{x}(t_f) \in \mathcal{X}_{g_1}, \quad (11d) $$

Then, the stability phase is obtained by solving the second OCP:

$$ \min_{u(\cdot)} J_2(\mathbf{x}, \mathbf{u}) = \int_{t_f}^{t_0} \ell(\mathbf{x}(t), \mathbf{u}(t))dt, \quad (12) $$

$$ \text{subject to: } \dot{x}(t) = f_{ha}(x(t), u(t)), \ t \in [t_f, t_s] \quad (13a) $$

$$ \mathbf{u}(t) \in \mathcal{V}_g, \mathbf{x}(t) \in \mathcal{X}_{g_2}, \ t \in [t_f, t_s] \quad (13b) $$

$$ \mathbf{x}(t_f) = \mathbf{x}_f, \quad (13c) $$

$$ \mathbf{x}(t_s) \in \mathcal{X}_{g_2}, \quad (13d) $$

In (11a) and (13a), $f_{ha}(\cdot)$ is the longitudinal dynamics as in (8). In (10)–(13), the timestamps $t_0, t_f$, and $t_s$ are the initial time, the terminal time, and the stability time, respectively; the sets $\mathcal{V} \subset \mathbb{R}^{n_u}, \mathcal{X}_{g_1} \subset \mathbb{R}^{n_s}, \mathcal{X}_{g_2} \subset \mathbb{R}^{n_s}, \mathcal{X}_{f_{g_1}} \subset \mathbb{R}^{n_s}, \mathcal{X}_{f_{g_2}} \subset \mathbb{R}^{n_s}$ are the input constraint set, state constraint set in the landing phase, state constraint set in the stability phase, terminal state constraint set in the landing phase, and terminal state constraint set in the stability phase, accordingly. They satisfy the conditions $\mathcal{X}_{f_{g_1}} \subset \mathcal{X}_{g_1}$ and $\mathcal{X}_{f_{g_2}} \subset \mathcal{X}_{g_2}$. The state $\mathbf{x}_0$ in (11c) is the steady-level flight state before the UAV starts executing the landing. The state $\mathbf{x}_f$ in (13c)

is $x(t_f)$ in (11d), that is, $x(t_f) = \mathbf{x}_f$, which is the terminal state of the UAV when it perches on the recovery net, obtained by solving (10)–(11). We also define the state $x(t_s) \triangleq \mathbf{x}_s$ in (13d) to be the state of the steady-level flight at the end of the stability phase.

Finally, in (10) and (12), we use a common running cost function $\ell(\mathbf{x}(t), \mathbf{u}(t))$ defined as

$$ \ell(\mathbf{x}(t), \mathbf{u}(t)) = \|\delta_t\|_{P_{\delta t}}^2 + \|\Delta\mathbf{u}\|_{P_{\Delta U}}^2. \quad (14) $$

Here, we want to minimize the thrust $\delta_t$ and the input deviation $\Delta\mathbf{u}$, with $P_{\delta t}$ and $P_{\Delta U}$ are a weighting number and a matrix. The input deviation $\Delta\mathbf{u}$ will be discussed in detail later in Section 5.1.

Let $(\mathbf{x}_r, \mathbf{u}_r)$ be the solution of the OCPs (10)–(13) by concatenating the solutions of the OCPs (10)–(11) and (12)–(13), when starting at a specific state $\mathbf{x}_0 \in \mathcal{X}_0 \subset \mathbb{R}^6$ at $t_0$. This is our reference for the tracking task. When the FW UAV starts at different initial positions $\mathbf{x}_0 \in \mathcal{X}_0$, we generate different trajectory-control input pairs $(\mathbf{x}_r, \mathbf{u}_r)$. Let us define the reference state set $\mathcal{X}_t$ to be a set of functions from $[t_0, t_s]$ to $\mathbb{R}^6$ that is a solution of the OCPs (10)–(13) for some $\mathbf{x}_0 \in \mathcal{X}_0$.

### 4. NMPC tracking framework for the stable deep-stall and perching landing

After obtaining a reference trajectory, we design a stable NMPC tracking scheme. Now, let us define the error state $\mathbf{x}_{e_t} \triangleq \mathbf{x}_t - \mathbf{x}_r$, and the error input $\mathbf{u}_{e_t} \triangleq \mathbf{u}_t - \mathbf{u}_{r_t}$. Linearizing the dynamics (1) around the reference trajectory $(\mathbf{x}_r, \mathbf{u}_r)$ using a first-order Taylor approximation gives the time-varying error dynamics

$$ \begin{aligned} \dot{\mathbf{x}}_{e_t} &= \frac{\partial f_{ha}}{\partial \mathbf{x}}(\mathbf{x}_{r_t}, \mathbf{u}_{r_t})\mathbf{x}_{e_t} + \frac{\partial f_{ha}}{\partial \mathbf{u}}(\mathbf{x}_{r_t}, \mathbf{u}_{r_t})\mathbf{u}_{e_t} + \Phi_t(\mathbf{x}_{r_t}, \mathbf{u}_{r_t}, \mathbf{x}_{e_t}, \mathbf{u}_{e_t}, \mu) \\ &= A_t \mathbf{x}_{e_t} + B_t \mathbf{u}_{e_t} + \Phi_t(\mathbf{x}_{r_t}, \mathbf{u}_{r_t}, \mathbf{x}_{e_t}, \mathbf{u}_{e_t}, \mu) \\ &\triangleq g_{ha}(t, \mathbf{x}_{e_t}, \mathbf{u}_{e_t}). \end{aligned} \quad (15) $$

Here, the remainder term $\Phi_t(\mathbf{x}_{r_t}, \mathbf{u}_{r_t}, \mathbf{x}_{e_t}, \mathbf{u}_{e_t}, \mu) \in \mathbb{R}^6$ has its $i$-th row component calculated by (Courant & John, 1989, Section 1.7.c.)

$$ \begin{aligned} \Phi_{t_f}(\mathbf{x}_{r_f}, \mathbf{u}_{r_f}, \mathbf{x}_{e_f}, \mathbf{u}_{e_f}, \mu) &= \frac{1}{2} \mathbf{x}_{e_f}^\top \mathbf{H}_{xx} f_{ha_l}(\mathbf{x}_{\mu}, \mathbf{u}_{\mu}) \mathbf{x}_{e_f} \\ &+ \frac{1}{2} \mathbf{u}_{e_f}^\top \mathbf{H}_{uu} f_{ha_l}(\mathbf{x}_{\mu}, \mathbf{u}_{\mu}) \mathbf{u}_{e_f} + \mathbf{x}_{e_f}^\top \mathbf{H}_{xu} f_{ha_l}(\mathbf{x}_{\mu}, \mathbf{u}_{\mu}) \mathbf{u}_{e_f}, \end{aligned} \quad (16) $$

where $(\mathbf{x}_{\mu}, \mathbf{u}_{\mu}) = (\mathbf{x}_{r_f} + \mu \mathbf{x}_{e_f}, \mathbf{u}_{r_f} + \mu \mathbf{u}_{e_f})$, for some $0 < \mu < 1$, $\mathbf{H}_{xx} f_{ha_l}(\mathbf{x}_{\mu}, \mathbf{u}_{\mu}) \in \mathbb{R}^{6 \times 6}$, $\mathbf{H}_{uu} f_{ha_l}(\mathbf{x}_{\mu}, \mathbf{u}_{\mu}) \in \mathbb{R}^{2 \times 2}$, and $\mathbf{H}_{xu} f_{ha_l}(\mathbf{x}_{\mu}, \mathbf{u}_{\mu}) \in \mathbb{R}^{6 \times 2}$ are Hessian matrices of the $i$-th row scalar component of the dynamics $f_{ha}$ as in (8), that is, $f_{ha_l}: \mathcal{X} \times \mathcal{U} \to \mathbb{R}$. Using the time-varying error dynamics (15), we formulate the NMPC trajectory tracking problem.

At the time instant $t$, with the error state $\mathbf{x}_{e_f}$, we solve an open-loop OCP (Grüne & Pannek, 2017):

$$ \min_{\bar{\mathbf{u}}(\cdot)} J \left( t, \mathbf{x}_{e_f}, \bar{\mathbf{x}}_{e}(\cdot), \bar{\mathbf{u}}_{e}(\cdot) \right) = \int_{t_f}^{t+T_p} L(\tau, \bar{\mathbf{x}}_{e}(\tau), \bar{\mathbf{u}}_{e}(\tau)) d\tau + F(t+T_p, \bar{\mathbf{x}}_{e}(t+T_p)), \quad (17) $$

$$ \text{subject to: } \dot{\bar{\mathbf{x}}}_{e}(\tau) = g_{ha}(\tau, \bar{\mathbf{x}}_{e}(\tau), \bar{\mathbf{u}}_{e}(\tau)), \quad \text{a.e. } \tau \in [t, t+T_p], \quad (18a) $$

$$ \bar{\mathbf{u}}_{e}(\tau) + \mathbf{u}_{r}(\tau) \in \mathcal{V}_{\text{MPC}}, \quad \text{a.e. } \tau \in [t, t+T_p], \quad (18b) $$

$$ \bar{\mathbf{x}}_{e}(\tau) + \mathbf{x}_{r}(\tau) \in \mathcal{X}_{\text{MPC}}, \quad \forall \tau \in [t, t+T_p], \quad (18c) $$

$$ \bar{\mathbf{x}}_{e}(t) = \mathbf{x}_{e_f}, \quad (18d) $$

$$ \bar{\mathbf{x}}_{e}(t+T_p) + \mathbf{x}_{r}(t+T_p) \in \mathcal{X}_{\text{MPC}}(t+T_p). \quad (18e) $$

In (17)–(18), the dynamical model $g_{ha}(\tau, \cdot, \cdot)$ as in (18a) is the time-varying error dynamics as in (15); the variable $\bar{\mathbf{x}}_{e}(\tau) = \mathbf{x}_{e}(\tau; t, \mathbf{x}_{e_f}, \bar{\mathbf{u}}_{e}(\cdot))$ is the predicted error state at the time $\tau$, starting from the error state $\mathbf{x}_{e_f}$ at $t$, under the predicted control input $\bar{\mathbf{u}}_{e}(\cdot)$, which is the optimal solution to the OCP. The variables $(\mathbf{x}_r, \mathbf{u}_r)$ are the optimal reference trajectory, which are the solutions of the OCPs (10)–(13). The set $\mathcal{V}_{\text{MPC}} \subset \mathbb{R}^{n_u}$ in (18b) is the input constraint set, the set $\mathcal{X}_{\text{MPC}} \subset \mathbb{R}^{n_s}$ in (18c) is the state constraint set, and the set $\mathcal{X}_{\text{MPC}}(t) \subset \mathbb{R}^{n_s}$ in (18e) is the terminal state
constraint set, while $T_p$ in (17) is the prediction horizon. Let $\mathbf{u}_e^*(\cdot)$ be the solution of (17)–(18). The first part of the solution, $\mathbf{u}_e^*(\tau)$, $\tau \in [t, t + \delta]$, is applied to the system, then the time shifts $t \leftarrow t + \delta$, and the problem (17)–(18) is solved repeatedly. The MPC law $\mathbf{u}_{\mathrm{e} \mathrm{MPC}}^*(\cdot)$ is constructed by the concatenation of the solutions of (17)–(18) that are applied to the system (Fontes, 2001; Grüne & Pannek, 2017).

We define respectively the time-varying stage cost $L$ and terminal cost $F$ in (17) as follows

$$L(\tau, \mathbf{x}_e(\tau), \mathbf{u}_e(\tau)) = \|\mathbf{x}_e(\tau)\|_Q^2 + \|\mathbf{u}_e(\tau)\|_R^2, \quad (19a)$$

$$F(t + T_p, \mathbf{x}_e(t + T_p)) = \|\mathbf{x}_e(t + T_p)\|_{P_t + T_p}^2. \quad (19b)$$

In (19a), we have weighting matrices $Q \in \mathbb{S}_{++}^6$ and $R \in \mathbb{S}_+^2$, while the time-varying matrix $P_t \in \mathbb{S}_{++}^6$ as in (19b) that we will define later in Section 4.1$^2$.

The terminal error state is constrained to be inside a time-varying terminal set $\mathcal{X}_{\mathrm{fMPC}}(t)$ defined as a terminal ellipsoid

$$\mathcal{X}_{\mathrm{fMPC}}(t) = \mathcal{E}_t(\mathbf{x}_{e_t}) \triangleq \{\mathbf{x}_{e_t} \in \mathbb{R}^6 | \mathbf{x}_{e_t}^\top P_t \mathbf{x}_{e_t} \le Y_t, Y_t > 0\}. \quad (20)$$

**Remark 1.** Note that, in this NMPC formulation, we have a time-varying terminal set $\mathcal{X}_{\mathrm{fMPC}}(t)$ due to the time-varying terminal weighting matrix $P_t$ and time-varying $Y_t$.

We now present the sufficient conditions for the stability of the NMPC scheme (17)–(18).

**Assumption 1.** The following conditions hold for all $t \ge t_0$:

(A) There is no model mismatch and there is no disturbance acting on the system;

(B) The input constraint set $\mathcal{V}_{\mathrm{MPC}}$ is compact and the state constraint set $\mathcal{X}_{\mathrm{MPC}}$ is closed;

(C) The time-varying terminal cost $F(t, \mathbf{x}_e(t))$ is continuously differentiable;

(D) There exists a $\mathcal{K}$-function $\kappa(\cdot)$ such that

$$\kappa(\|\mathbf{x}_e(t)\|) \le L(t, \mathbf{x}_e(t), \mathbf{u}_e(t)); \quad (21)$$

(E) The time-varying terminal ellipsoid $\mathcal{X}_{\mathrm{fMPC}}(t)$ satisfies $\mathcal{X}_{\mathrm{fMPC}}(t) \oplus \mathcal{X}_t \in \mathcal{X}_{\mathrm{MPC}}$, $\mathcal{X}_{\mathrm{fMPC}}(t)$ is a closed set, and the equilibrium state $0 \in \mathcal{X}_{\mathrm{fMPC}}(t)$;

(F) There exists a time-varying local controller $\mathbf{u}_{e\mathrm{LOC}}(t, \mathbf{x}_{e_t})$ satisfies $\mathbf{u}_{e\mathrm{LOC}}(t, \mathbf{x}_{e_t}) + \mathbf{u}_{P_t} \in \mathcal{V}_{\mathrm{MPC}}$, for all $\mathbf{x}_{e_t} \in \mathcal{X}_{\mathrm{fMPC}}(t)$;

(G) The time-varying terminal set $\mathcal{X}_{\mathrm{fMPC}}(t)$ is forward invariant under the local controller $\mathbf{u}_{e\mathrm{LOC}}(\tau, \mathbf{x}_{e_t})$, i.e.,

$$\forall \mathbf{x}_{e_t} \in \mathcal{X}_{\mathrm{fMPC}}(t), \mathbf{x}_e(\tau; t, \mathbf{x}_{e_t}, \mathbf{u}_{e\mathrm{LOC}}(\tau, \mathbf{x}_{e_t}) \in \mathcal{X}_{\mathrm{fMPC}}(\tau), \text{ for all } \tau \ge t; \quad (22)$$

(H) The terminal cost is a local Lyapunov function: for all $\tau \ge t$ and for all $\mathbf{x}_e(\tau; t, \mathbf{x}_{e_t}, \mathbf{u}_{e\mathrm{LOC}}(\tau, \mathbf{x}_{e_t})) \in \mathcal{X}_{\mathrm{fMPC}}(\tau)$:

$$\frac{d F(\tau, \mathbf{x}_e(\tau))}{dt} \le -L(\tau, \mathbf{x}_e(\tau), \mathbf{u}_{e\mathrm{LOC}}(\tau, \mathbf{x}_{e_t})). \quad (23)$$

**Theorem 1.** Let the Assumption 1 hold, then the NMPC scheme (17)–(18) achieves the recursive feasibility and the closed-loop error dynamics $\dot{x}_e(t) = g_{ha}(t, \mathbf{x}_e(t), \mathbf{u}_e(t))$ are asymptotically stable.

**Proof.** The result is a modification of Faulwasser (2013, Corollary 3.2), augmented from a time-varying NMPC stabilization problem (Faulwasser, 2013, Theorem 2.1) with its detailed proof in Faulwasser (2013, Appendix A) and Fontes (2001). Here, we briefly present the sketch of the proof, where the difference here is that we are considering the solution of the NMPC scheme as the input error, not just the input.

**Recursive feasibility:**

Let $\mathbf{u}_e^*(\cdot) : [t, t + T_p] \to \mathcal{V}_{\mathrm{MPC}}$ be the solution (optimal predicted input) of the OCP (17)–(18) at the time instant $t$. The equivalent optimal predicted state is $\mathbf{x}_e^*(\tau) = \mathbf{x}_e(\tau; t, \mathbf{x}_{e_t}, \mathbf{u}_e^*(\cdot))$, for $\tau \in [t, t + T_p]$. For $\tau \in [t, t + \delta]$,

applying the optimal predicted input to the system while assuming there is no model mismatch and no external disturbance, as in Assumption 1(A), the real state at the time instant $t + T_p$, which is the optimal predicted state at $t + T_p$, $\mathbf{x}_e(t + T_p) = \mathbf{x}_e^*(t + T_p) = \mathbf{x}_e(t + T_p; t, \mathbf{x}_{e_t}, \mathbf{u}_e^*(\cdot))$, is inside the terminal set $\mathcal{X}_{\mathrm{fMPC}}(t + T_p)$.

Now, we solve the OCP (17)–(18) at the time instant $t + \delta$, an admissible (not necessarily optimal) input is

$$\tilde{\mathbf{u}}_e(\tau) = \begin{cases} \mathbf{u}_e^*(\tau), & \text{for } \tau \in [t + \delta, t + T_p] \\ \mathbf{u}_{e\mathrm{LOC}}(\tau, \mathbf{x}_{e_t}) & \text{choose} \\ \mathbf{u}_{e\mathrm{LOC}}(t + T_p, \mathbf{x}_{e_t + T_p}) & \text{for } \tau \in [t + T_p, t + T_p + \delta]. \end{cases} \quad (24)$$

Using Assumption 1(G):

$$\begin{align} \forall \mathbf{x}_{e_t + T_p} \in \mathcal{X}_{\mathrm{fMPC}}(t + T_p), \nonumber \\ \mathbf{x}_e(\tau; t + T_p, \mathbf{x}_{e_t + T_p}, \mathbf{u}_{e\mathrm{LOC}}(\tau, \mathbf{x}_{e_t}) \in \mathcal{X}_{\mathrm{fMPC}}(\tau), \text{ for all } \tau \ge t + T_p, \tag{25} \end{align}$$

hence, $\mathbf{x}_e(t + T_p + \delta; t + T_p, \mathbf{x}_{e_t + T_p}, \mathbf{u}_{e\mathrm{LOC}}(t + T_p, \mathbf{x}_{e_t + T_p}) \in \mathcal{X}_{\mathrm{fMPC}}(t + T_p + \delta)$. The admissible state equivalent to the feasible input (24) is

$$\bar{\mathbf{x}}_e(\tau) = \begin{cases} \mathbf{x}_e^*(\tau), & \text{for } \tau \in [t + \delta, t + T_p] \\ \mathbf{x}_e(\tau; t + T_p, \mathbf{x}_{e_t + T_p}, \mathbf{u}_{e\mathrm{LOC}}(\tau, \mathbf{x}_{e_t}) & \text{for } \tau \in [t + T_p, t + T_p + \delta]. \end{cases} \quad (26)$$

Using induction for other time instants along the reference trajectory, we obtain the proof of recursive feasibility.

**Asymptotic stability:**

Let $t_k = t_0 + k\delta$, for $k \in \mathbb{N}$, and consider the cost functional

$$J(t_k, \mathbf{x}_{e_t k}, \bar{\mathbf{x}}_{e_t k}(\cdot), \bar{\mathbf{u}}_{e_t k}(\cdot)) = \int_{t_k}^{t_k + T_p} L(\tau, \bar{\mathbf{x}}_{e_t k}(\tau), \mathbf{u}_{e_t k}^*(\tau)) d\tau + F(t_k + T_p, \bar{\mathbf{x}}_{e_t k}^*(t_k + T_p)). \quad (27)$$

Let us define the optimal value function of the OCP (17)–(18) for $t \in [t_k, t_{k+1}]$ as

$$V_{t_k}(t, \mathbf{x}_{e_t}) \triangleq \int_t^{t_k + T_p} L(\tau, \mathbf{x}_{e_t}^*(\tau), \mathbf{u}_{e_t}^*(\tau)) d\tau + F(t_k + T_p, \mathbf{x}_{e_t}^*(t_k + T_p)). \quad (28)$$

The optimal value function when $t = t_k$ is

$$\begin{align} V_{t_k}(t, \mathbf{x}_{e_t k}) &= \int_{t_k}^{t_k + T_p} L(\tau, \mathbf{x}_{e_t k}^*(\tau), \mathbf{u}_{e_t k}^*(\tau)) d\tau + F(t_k + T_p, \mathbf{x}_{e_t k}^*(t_k + T_p)) \nonumber \\ &= \underbrace{\int_{t_k}^{t} L(\tau, \mathbf{x}_{e_t k}^*(\tau), \mathbf{u}_{e_t k}^*(\tau)) d\tau}_{\ge 0} + \int_{t}^{t_k + T_p} L(\tau, \mathbf{x}_{e_t k}^*(\tau), \mathbf{u}_{e_t k}^*(\tau)) d\tau \nonumber \\ &\quad + F(t_k + T_p, \mathbf{x}_{e_t k}^*(t_k + T_p)) \ge V_{t_k}(t, \mathbf{x}_{e_t}). \tag{29} \end{align}$$

Therefore, $V_{t_k}(t, \mathbf{x}_{e_t})$ is decreasing for $t \in [t_k, t_k + \delta]$ and $\forall k \in \mathbb{N}$. Next, let us consider the time instant $t_{k+1} = t_k + \delta$. We have

$$V_{t_{k+1}}(t_{k+1}, \mathbf{x}_{e_{t_{k+1}}}) \le J(t_{k+1}, \mathbf{x}_{e_{t_{k+1}}}, \bar{\mathbf{x}}_{e_{t_{k+1}}}(\cdot), \bar{\mathbf{u}}_{e_{t_{k+1}}}(\cdot)), \quad (30)$$

in which

$$\begin{align} J(t_{k+1}, \mathbf{x}_{e_{t_{k+1}}}, \bar{\mathbf{x}}_{e_{t_{k+1}}}(\cdot), \bar{\mathbf{u}}_{e_{t_{k+1}}}(\cdot)) & \nonumber \\ &= \int_{t_{k+\delta}}^{t_{k+\delta+T_p}} L(\tau, \bar{\mathbf{x}}_{e_{t_{k+1}}}(\tau), \bar{\mathbf{u}}_{e_{t_{k+1}}}(\tau)) d\tau + F(t_k + \delta + T_p, \bar{\mathbf{x}}_{e_{t_{k+1}}}(t_k + \delta + T_p)) \nonumber \\ &= \int_{t_{k+\delta}}^{t_{k+T_p}} L(\tau, \bar{\mathbf{x}}_{e_{t_{k+1}}}(\tau), \bar{\mathbf{u}}_{e_{t_{k+1}}}(\tau)) d\tau + \int_{t_{k+T_p}}^{t_{k+\delta+T_p}} L(\tau, \bar{\mathbf{x}}_{e_{t_{k+1}}}(\tau), \bar{\mathbf{u}}_{e_{t_{k+1}}}(\tau)) d\tau \nonumber \\ &\quad + F(t_k + \delta + T_p, \bar{\mathbf{x}}_{e_{t_{k+1}}}(t_k + \delta + T_p)). \tag{31} \end{align}$$

For $\tau \in [t_k + \delta, t_k + \delta + T_p]$, $\mathbf{u}_{e_t k}^*(\tau) = \bar{\mathbf{u}}_{e_{t_{k+1}}}(\tau)$, as in (24), and $\bar{\mathbf{x}}_{e_{t_{k+1}}}(\tau; t_{k+1}, \mathbf{x}_{e_{t_{k+1}}}, \bar{\mathbf{u}}_{e_{t_{k+1}}}(\cdot)) = \mathbf{x}_{e_t}^*(\tau; t_k, \mathbf{x}_{e_t}, \mathbf{u}_e^*(\cdot))$, as in (26), so

<footnote n="2">2 We will employ the subscript $\square$ for all time-varying matrices for readability.</footnote>
that

$$
\begin{align}
& V_{t_{k+1}}(t_{k+1}, \mathbf{x}_{e_{t_{k+1}}}) - V_{t_k}(t_k, \mathbf{x}_{e_{t_k}}) \nonumber \\
& \le \int_{t_k + T_p}^{t_{k+\delta} + T_p} L(\tau, \bar{\mathbf{x}}_{e_{t_{k+1}}}(\tau), \bar{\mathbf{u}}_{e_{t_{k+1}}}(\tau)) d\tau + F(t_k + \delta + T_p, \bar{\mathbf{x}}_{e_{t_{k+1}}}(t_k + \delta + T_p)) \nonumber \\
& \quad - F(t_k + T_p, \bar{\mathbf{x}}_{e_{t_k}}(t_k + T_p)) - \int_{t_k}^{t_{k+\delta}} L(\tau, \mathbf{x}_{e_{t_k}}^*(\tau), \mathbf{u}_{e_{t_k}}^*(\tau)) d\tau \nonumber \\
& \le - \int_{t_k}^{t_{k+\delta}} L(\tau, \mathbf{x}_{e_{t_k}}^*(\tau), \mathbf{u}_{e_{t_k}}^*(\tau)) d\tau \le - \int_{t_k}^{t_{k+\delta}} \kappa(\|\mathbf{x}_e(\tau; t_k, \mathbf{x}_{e_{t_k}}, \mathbf{u}_e^*(\cdot))\|) d\tau. \tag{32}
\end{align}
$$

The first three terms in the second inequality in (32) cancel out due to Assumption 1(H).
Let us define the MPC value function

$$
V^{\delta}(t, \mathbf{x}_{e_{t_k}}) \triangleq V_{t_k}(t_k, \mathbf{x}_{e_{t_k}}) - \int_{t_k}^{t} L(\tau, \mathbf{x}_{e_{t_k}}^{*}(\tau), \mathbf{u}_{e_{t_k}}^{*}(\tau)) d\tau \quad (33)
$$

with $t_k = \max_k \{t_k \mid t_k \le t\}$. Using the inequality (32) iteratively, we ob-
tain

$$
\begin{align}
V^{\delta}(t_k, \mathbf{x}_{e_{t_k}}) - V^{\delta}(t_0, \mathbf{x}_{e_{t_0}}) &= V_{t_k}(t_k, \mathbf{x}_{e_{t_k}}) - V_{t_0}(t_0, \mathbf{x}_{e_{t_0}}) \nonumber \\
&\le - \sum_{i=0}^{k} \int_{t_i}^{t_{i+1}} \kappa(\|\mathbf{x}_e(\tau; t_i, \mathbf{x}_{e_{i_i}}, \mathbf{u}_e^*(\cdot))\|) d\tau, \tag{34}
\end{align}
$$

which leads to

$$
V^{\delta}(t_k, \mathbf{x}_{e_{t_k}}) - V^{\delta}(t_0, \mathbf{x}_{e_{t_0}}) \le - \int_{t_0}^{t_k} \kappa (\|\mathbf{x}_e(\tau; t_0, \mathbf{x}_{e_{t_0}}, \mathbf{u}_{e_{\text{MPC}}}^*(\cdot))\|) d\tau, \quad (35)
$$

hence

$$
V^{\delta}(t, \mathbf{x}_{e_t}) - V^{\delta}(t_0, \mathbf{x}_{e_{t_0}}) \le - \int_{t_0}^{t} \kappa (\|\mathbf{x}_e(\tau; t_0, \mathbf{x}_{e_{t_0}}, \mathbf{u}_{e_{\text{MPC}}}^*(\cdot))\|) d\tau, \quad \forall t \ge t_0. \quad (36)
$$

Since $V^{\delta}(t_0, \mathbf{x}_{e_{t_0}}) < \infty$ and $V^{\delta}(t, \mathbf{x}_{e_t}) > 0$, then $V^{\delta}(t, \mathbf{x}_{e_t}) < \infty$, we obtain

$$
\int_{t_0}^{t} \kappa (\| \mathbf{x}_e(\tau; t_0, \mathbf{x}_{e_{t_0}}, \mathbf{u}_{e_{\text{MPC}}}^*(\cdot)) \|) d\tau < \infty, \quad \forall t \ge t_0. \quad (37)
$$

Finally, applying the following Lemma, we can achieve the asymptotic
stability.

**Lemma 1.** *Michalska and Vinter (1994, Lemma 4)* Let $h(\cdot) : \mathbb{R}^{n_s} \rightarrow \mathbb{R}_{\ge 0}$ be a continuous, positive definite function and $x$ be an absolutely continuous function on $\mathbb{R}$. If $\|x(\cdot)\| \in L^{\infty}(0, \infty) < \infty$, $\|\dot{x}(\cdot)\| \in L^{\infty}(0, \infty) < \infty$, and $\lim_{t \to \infty} \int_0^t h(x(\tau)) d\tau < \infty$, then $x(t) \to 0$ as $t \to \infty$.

Therefore, $\lim_{t \to \infty} \| \mathbf{x}_e(t; t_0, \mathbf{x}_{e_{t_0}}, \mathbf{u}_{e_{\text{MPC}}}^*(\cdot)) \| = 0$. $\square$

After obtaining the reference trajectory, we briefly present steps to
design terminal ingredients, namely the terminal cost $F(t, \mathbf{x}_{e_t})$ as in (19b)
and the terminal set $\mathcal{X}_{f_{\text{MPC}}}(t)$ as in (20), for our NMPC tracking scheme
with guaranteed stability for the FW UAV in the deep-stall and perchling
maneuver.

The condition 1(B) can be satisfied with an appropriate choice of
the constraint sets. The condition 1(C) is satisfied with the terminal cost
defined as in (19b). The condition 1(D) is satisfied with the stage cost
$L(\cdot, \cdot)$ defined as in (19a):

$$
L(t, \mathbf{x}_e(t), \mathbf{u}_e(t)) = \|\mathbf{x}_e(t)\|_Q^2 + \|\mathbf{u}_e(t)\|_R^2 \ge \|\mathbf{x}_e(t)\|_Q^2 \ge \lambda_{\min}(Q) \|\mathbf{x}_e(t)\|^2. \quad (38)
$$

Now, we move forward to design the terminal cost and terminal
set that satisfy the conditions 1(E)–1(H). First, the local controller
$\mathbf{u}_{e_{\text{LOC}}}(t, \mathbf{x}_{e_t})$ as in 1(F)–1(H) is chosen as a feedback controller

$$
\mathbf{u}_{e_{\text{LOC}}}(t, \mathbf{x}_{e_t}) \triangleq K_t \mathbf{x}_{e_t}, \quad (39)
$$

with $K_t \in \mathbb{R}^{n_s \times n_u}$ to be defined. In Section 4.1, we solve an optimization problem and utilize a linear matrix inequality (LMI) obtained from the

local Lyapunov function to find the time-varying matrices $K_t$ as in (39) and the time-varying positive definite matrices $P_t$ as in (20) satisfying the condition 1(H), which implies the fulfillment of the condition 1(G). Next, in Section 4.2, we find the time-varying variable $\Upsilon_{1t} > 0$ so that the time-varying terminal ellipsoid satisfies the conditions 1(E) and 1(F). Finally, due to the model mismatch between the real dynamics and the linearized dynamics as in (15), we have to re-check the condition (23) in Section 4.3 to find the variable $\Upsilon_t \in (0, \Upsilon_{1t}]$ as in (20).

4.1. Stability proof (find P<sub>t</sub>)

We first formulate a time-varying ellipsoid

$$
\mathcal{X}_{f_{\text{MPC}}}^{1}(t) = \mathcal{E}_{t}^{1}(\mathbf{x}_{e}) \triangleq \{\mathbf{x}_{e} \in \mathbb{R}^{6} \mid \mathbf{x}_{e}^{\top} P_{t} \mathbf{x}_{e} \le \Upsilon_{1t}\}, \quad (40)
$$

with time-varying $\Upsilon_{1t} > 0$. This ellipsoid will **temporarily** serve as the terminal set. We now proceed to find the time-varying terminal matrices $P_t$ and the time-varying variable $\Upsilon_{1t} > 0$. The design of the matrices $P_t$ is based on the work of Köhler et al. (2020a), by solving an optimization problem containing an LMI.

**Proposition 1.** *The terminal weighting matrices $P_t$ as in (19b), (40) and the feedback gain matrices $K_t$ in the local controller $\mathbf{u}_{e_{\text{LOC}}}(t, \mathbf{x}_{e_t}) = K_t \mathbf{x}_{e_t}$ as in Assumptions 1(F)–1(H), (39) are obtained from the optimization problem:*

$$
\min_{X_t, Y_t, X_{\min}} - \log \det(X_{\min}), \quad (41a)
$$

$$
\begin{equation}
\begin{split}
\text{s.t.} \quad \left[ \begin{array}{@{}cccc@{}}
((A_t + \epsilon I_n)X_t + B_t Y_t)^\top & & & \\
Q^{1/2} X_t & * & * & \\
R^{1/2} Y_t & & * & \\
& \ddots & \ddots & \\
& & \ddots & \\
& & & -\mathbf{1}_{n_s} \\
& & & \mathbf{0} \\
& & & -\mathbf{1}_{n_u}
\end{array} \right] \\
&\le 0, \quad \text{for } t \in [t_0, t_s];
\end{split}
\tag{41b}
\end{equation}
$$

$$
0 < X_{\min} \le X_t, \quad \text{for } t \in [t_0, t_s]; \tag{41c}
$$

$$
\frac{d X_{t_s}}{dt} = 0, \tag{41d}
$$

with $P_t = X_t^{-1}$, $K_t = Y_t P_t$, and $*$ represents the equivalent transpose compo-
nents of a matrix.

**Proof.** First, let us recall the condition of the local Lyapunov condition as in 1(H) for $\tau \ge t \ge t_0$:

$$
\frac{d F(\tau, \mathbf{x}_{e_t})}{dt} \le -L(\tau, \mathbf{x}_e(\tau), \mathbf{u}_{e_{\text{LOC}}}(\tau, \mathbf{x}_{e_t})). \quad (42)
$$

Plugging the local controller as in (39),

$$
\mathbf{u}_{e_t} = \mathbf{u}_{e_{\text{LOC}}}(t, \mathbf{x}_{e_t}) = K_t \mathbf{x}_{e_t}, \quad (43)
$$

into the linearized dynamics (15) gives

$$
\begin{align}
\dot{\mathbf{x}}_{e_t} &= A_t \mathbf{x}_{e_t} + B_t K_t \mathbf{x}_{e_t} + \Phi_t (\mathbf{x}_{r_t}, \mathbf{u}_{r_t}, \mathbf{x}_{e_t}, K_t \mathbf{x}_{e_t}, \mu) \nonumber \\
&= (A_t + B_t K_t) \mathbf{x}_{e_t} + \Phi_t (\mathbf{x}_{r_t}, \mathbf{u}_{r_t}, \mathbf{x}_{e_t}, K_t \mathbf{x}_{e_t}, \mu). \tag{44}
\end{align}
$$

Taking the time derivative of the terminal cost (19b) gives

$$
\begin{align}
\dot{F}(t, \mathbf{x}_{e_t}) &= \dot{\mathbf{x}}_{e_t}^{\top} P_t \mathbf{x}_{e_t} + \mathbf{x}_{e_t}^{\top} \dot{P}_t \mathbf{x}_{e_t} + \mathbf{x}_e^{\top} P_t \mathbf{x}_e \nonumber \\
&= \mathbf{x}_{e_t}^{\top} \left[ (A_t + B_t K_t)^{\top} P_t + P_t (A_t + B_t K_t) \right] \mathbf{x}_{e_t} \nonumber \\
&\quad + 2 \mathbf{x}_{e_t}^{\top} P_t \Phi_t (\mathbf{x}_{r_t}, \mathbf{u}_{r_t}, \mathbf{x}_{e_t}, K_t \mathbf{x}_{e_t}, \mu) + \mathbf{x}_{e_t}^{\top} \dot{P}_t \mathbf{x}_{e_t}, \tag{45}
\end{align}
$$

here, we employ the fact that the term $\Phi_t^\top(\mathbf{x}_{e_t}) P_t \mathbf{x}_{e_t}$ is a scalar, hence,

$$
\begin{align}
\Phi_t^\top(\mathbf{x}_{r_t}, \mathbf{u}_{r_t}, \mathbf{x}_{e_t}, K_t \mathbf{x}_{e_t}, \mu) P_t \mathbf{x}_{e_t} &= [\Phi_t^\top(\mathbf{x}_{r_t}, \mathbf{u}_{r_t}, \mathbf{x}_{e_t}, K_t \mathbf{x}_{e_t}, \mu) P_t \mathbf{x}_{e_t}]^\top \nonumber \\
&= \mathbf{x}_{e_t}^\top P_t \Phi_t(\mathbf{x}_{r_t}, \mathbf{u}_{r_t}, \mathbf{x}_{e_t}, K_t \mathbf{x}_{e_t}, \mu). \tag{46}
\end{align}
$$

The quadratic term $\mathbf{x}_{e_t}^\top P_t \Phi_t(\mathbf{x}_{r_t}, \mathbf{u}_{r_t}, \mathbf{x}_{e_t}, K_t \mathbf{x}_{e_t}, \mu)$ is bounded from above as

$$
\begin{align}
& \mathbf{x}_{e_t}^\top P_t \Phi_t(\mathbf{x}_{r_t}, \mathbf{u}_{r_t}, \mathbf{x}_{e_t}, K_t \mathbf{x}_{e_t}, \mu) \le \|\mathbf{x}_{e_t}^\top P_t\| \|\Phi_t(\mathbf{x}_{r_t}, \mathbf{u}_{r_t}, \mathbf{x}_{e_t}, K_t \mathbf{x}_{e_t}, \mu)\| \nonumber \\
& \le \|P_t\| \left[ \sup \frac{\Phi_t(\mathbf{x}_{r_t}, \mathbf{u}_{r_t}, \mathbf{x}_{e_t}, K_t \mathbf{x}_{e_t}, \mu)}{\|\mathbf{x}_{e_t}\|} \right] \|\mathbf{x}_{e_t}\|^2 \le \epsilon \mathbf{x}_{e_t}^\top P_t \mathbf{x}_{e_t}, \tag{46}
\end{align}
$$
in which

$$ \sup \frac{\Phi_t(\mathbf{x}_{t_r}, \mathbf{u}_{t_r}, \mathbf{x}_{e_t}, K_t \mathbf{x}_{e_t}, \mu)}{\|\mathbf{x}_{e_t}\|} \le \epsilon \frac{\lambda_{\min}(P_t)}{\|P_t\|} \quad \text{and} \quad \lambda_{\min}(P_t) \|\mathbf{x}_{e_t}\|^2 \le \mathbf{x}_{e_t}^\top P_t \mathbf{x}_{e_t}, \qquad (47) $$

which leads to

$$ \Phi_t(\mathbf{x}_{t_r}, \mathbf{u}_{t_r}, \mathbf{x}_{e_t}, K_t \mathbf{x}_{e_t}, \mu) \le \epsilon \mathbf{x}_{e_t}, \qquad (48) $$

With an upper bound for the remainder term, $\Phi_t(\mathbf{x}_{t_r}, \mathbf{u}_{t_r}, \mathbf{x}_{e_t}, K_t \mathbf{x}_{e_t}, \mu) \le \epsilon \mathbf{x}_{e_t}$, the time derivative of the terminal cost (45) is now bounded from above as

$$ \begin{aligned} \dot{F}(t, \mathbf{x}_{e_t}) &\le \mathbf{x}_{e_t}^\top (A_t + B_t K_t)^\top P_t + P_t (A_t + B_t K_t)] \mathbf{x}_{e_t} + 2 \mathbf{x}_e^\top P_t \epsilon \mathbf{x}_{e_t} + \mathbf{x}_e^\top \dot{P}_t \mathbf{x}_{e_t} \\ &= \mathbf{x}_{e_t}^\top [(A_t + B_t K_t + \epsilon \mathbb{1}_{n_s})^\top P_t + P_t (A_t + B_t K_t + \epsilon \mathbb{1}_{n_s})] \mathbf{x}_{e_t} + \mathbf{x}_e^\top \dot{P}_t \mathbf{x}_{e_t}, \end{aligned} \qquad (49) $$

where $n_s$ is the number of states. Therefore, the Lyapunov condition (42) is satisfied if and only if

$$ (A_t + B_t K_t + \epsilon \mathbb{1}_{n_s})^\top P_t + P_t (A_t + B_t K_t + \epsilon \mathbb{1}_{n_s}) + d P_t / dt \le -(Q + K_t^\top R K_t). \qquad (50) $$

Using the changes of variables $X_t = P_t^{-1}$ and $Y_t = K_t X_t$, (i.e., $K_t = Y_t X_t^{-1} = Y_t P_t$, hence, $K_t^\top = X_t^{-1} Y_t^\top$), the inequality (50) turns into

$$ (A_t + B_t Y_t X_t^{-1} + \epsilon \mathbb{1}_{n_s})^\top X_t^{-1} + X_t^{-1} (A_t + B_t Y_t X_t^{-1} + \epsilon \mathbb{1}_{n_s}) + \frac{d X_t^{-1}}{dt} + (Q + X_t^{-1} Y_t^\top R Y_t X_t^{-1}) \le 0. \qquad (51) $$

After some algebraic manipulations, taking into account the facts that $X_t^{-1} X_t = X_t X_t^{-1} = \mathbb{1}_{n_s}$, $X_t^\top = X_t$, and from the time derivative of $X_t X_t^{-1} = \mathbb{1}_{n_s}$, $\frac{d X_t^{-1}}{dt} = -X_t^{-1} \frac{d X_t}{dt} X_t^{-1}$, hence, $X_t^\top \frac{d X_t^{-1}}{dt} X_t = -\frac{d X_t}{dt}$, we obtain

$$ ((A_t + \epsilon \mathbb{1}_{n_s})X_t + B_t Y_t)^\top + ((A_t + \epsilon \mathbb{1}_{n_s})X_t + B_t Y_t) - \frac{d X_t}{dt} - \begin{bmatrix} Q^{1/2} X_t \\ R^{1/2} Y_t \end{bmatrix}^\top \begin{bmatrix} -\mathbb{1}_{n_s} & 0 \\ 0 & -\mathbb{1}_{n_u} \end{bmatrix} \begin{bmatrix} Q^{1/2} X_t \\ R^{1/2} Y_t \end{bmatrix} \le 0, \qquad (52) $$

here, $n_s$ is the number of states.

Finally, to find the terminal weighting matrices $P_t$ as in (19b) and (40), we solve the optimization problem (41). The values of the feedback gain matrices $K_t$ as in (43) are followed by $K_t = Y_t P_t$. The reasons for the problem formulation (41) are presented next.

The LMI (41b) is obtained directly from the nonlinear equality constraint (52).

Now, consider $n$ ellipsoids, $\mathcal{E}_1 = \{\mathbf{x} | \mathbf{x}^\top X^{-1} \mathbf{x} \le 1\}$, $\mathcal{E}_2 = \{\mathbf{x} | \mathbf{x}^\top X^{-1} \mathbf{x} \le 1\}$, ..., $\mathcal{E}_n = \{\mathbf{x} | \mathbf{x}^\top X^{-1} \mathbf{x} \le 1\}$, and the ellipsoid $\mathcal{E}_{\min} = \{\mathbf{x} | \mathbf{x}^\top X^{-1} \mathbf{x} \le 1\}$, with the symmetric positive definite matrices $X_1, X_2, \dots, X_{\min} \in \mathbb{S}_{+\min}^n$. We want there to exist the ellipsoid $\mathcal{E}_{\min}$ inscribing in the ellipsoids $\mathcal{E}_1, \mathcal{E}_2, \dots, \mathcal{E}_n$, that is, $\mathcal{E}_{\min} \subseteq \mathcal{E}_1, \mathcal{E}_{\min} \subseteq \mathcal{E}_2, \dots, \mathcal{E}_{\min} \subseteq \mathcal{E}_n$, which are satisfied if and only if $0 < X_{\min} \le X_1, 0 < X_{\min} \le X_2, \dots, 0 < X_{\min} \le X_n$ (Boyd & Vandenberghe, 2004, Example 2.18). This explains the inequality constraint (41c).

Next, in Grötschel et al. (1988, Chapter 3), the volume of an ellipsoid centered at $\mathbf{x}_0$,

$$ \mathcal{E} = \{\mathbf{x} \in \mathbb{R}^{n_s} | (\mathbf{x} - \mathbf{x}_0)^\top P(\mathbf{x} - \mathbf{x}_0) \le 1\}, \qquad (53) $$

is calculated as:

$$ \text{vol}(\mathcal{E}) = \frac{\pi^{n_s/2}}{\Gamma(n_s/2 + 1)} \times \frac{1}{\sqrt{\det(P)}}, \qquad (54) $$

where $\Gamma(n_s) = \int_0^\infty e^{-t} t^{n_s-1} dt$, for $n_s > 0$.

We would also want to have the largest ellipsoid $\mathcal{E}_{\min}$ possible, which means the smallest det($P_{\min}$) (from (54)), hence, the objective is to minimize $\text{logdet}(P_{\min}) = \text{logdet}(X_{\min}^{-1}) = -\text{logdet}(X_{\min})$ (c.f. Boyd & Vandenberghe, 2004, Section 8.4.2). This explains the objective function (41a).

Finally, when $t \ge t_s$, the matrices $A_t$ and $B_t$ are constant (since the FW UAV is in steady-level flight), we want the time derivative of the matrix $P_t$ to be zero, so we enforce the condition (41d). The proof is complete. □

After obtaining the symmetric positive definite terminal matrix $P_t$, we conclude that the set $\mathcal{E}_t^1(\mathbf{x}_e) = \{\mathbf{x}_e \in \mathbb{R}^6 | \mathbf{x}_e^\top P_t \mathbf{x}_e \le Y_{1,t}\}$ is forward invariant. Next, we move forward to find the value of $Y_{1,t}$, such that the ellipsoid $\mathcal{E}_t^1(\mathbf{x}_e)$ satisfies the condition 1(E).

## 4.2. Constraint satisfaction (find $Y_{1,t}$)

Consider all the terminal error states that belong to the ellipsoid (20)

$$ \begin{aligned} \mathcal{E}_t^1(\mathbf{x}_{e_t}) &= \{\mathbf{x}_{e_t} \in \mathbb{R}^{n_s} | \mathbf{x}_{e_t}^\top P_t \mathbf{x}_{e_t} \le Y_{1,t}\} \\ &= \{\mathbf{x}_{e_t} \in \mathbb{R}^{n_s} | \mathbf{x}_{e_t} = \sqrt{Y_{1,t}} P_t^{-1/2} \mathbf{v}_t, \|\mathbf{v}_t\| \le 1, \mathbf{v}_t \in \mathbb{R}^{n_s}\}. \end{aligned} \qquad (55) $$

We want to find the largest $\mathcal{E}_t^1(\mathbf{x}_{e_t})$, which is analogous to finding the largest $Y_{1,t}$, such that $\forall \mathbf{x}_{e_t} \in \mathcal{E}_t^1$, $\mathbf{u}_t = K_t \mathbf{x}_{e_t} + \mathbf{u}_{r_t} \in \mathcal{V}_{\text{MPC}}$, and $\mathbf{x}_t = \mathbf{x}_{e_t} + \mathbf{x}_{r_t} \in \mathcal{X}_{\text{MPC}}$, with the state and input constraint sets written as polytopes

$$ \mathcal{X}_{\text{MPC}} \triangleq \{\mathbf{x} \in \mathbb{R}^{n_s} | \mathbf{a}_x^\top \mathbf{x} \le \mathbf{b}_x\}, \qquad (56a) $$

$$ \mathcal{V}_{\text{MPC}} \triangleq \{\mathbf{u} \in \mathbb{R}^{n_u} | \mathbf{a}_u^\top \mathbf{u} \le \mathbf{b}_u\}. \qquad (56b) $$

Plugging $\mathbf{x}_t$, $\mathbf{u}_t$, and $\mathbf{x}_{e_t}$ (55) into (56) yields the following constraints

$$ \sup_{\|\mathbf{v}_t\| \le 1} \begin{cases} \mathbf{a}_x^\top [\mathbf{x}_{t_r} + \sqrt{Y_{1,t}} P_t^{-1/2} \mathbf{v}_t] \le \mathbf{b}_x, \\ \mathbf{a}_u^\top [\mathbf{u}_{t_r} + K_t(\sqrt{Y_{1,t}} P_t^{-1/2} \mathbf{v}_t)] \le \mathbf{b}_u. \end{cases} \qquad (57) $$

Therefore, the problem of finding the largest ellipsoid $\mathcal{E}_t^1(\mathbf{x}_{e_t})$ is equivalent to an optimization problem:

$$ \min_{Y_{1,t}} -Y_{1,t}, \qquad (58a) $$

$$ \text{s.t.} \quad \left\{ \begin{array}{ll} \| \sqrt{Y_{1,t}} P_t^{-1/2} \mathbf{a}_x \| + \mathbf{a}_x^\top \mathbf{x}_{t_r} & \le \mathbf{b}_x, \\ \| \sqrt{Y_{1,t}} P_t^{-1/2} K_t^\top \mathbf{a}_u \| + \mathbf{a}_u^\top \mathbf{u}_{t_r} & \le \mathbf{b}_u, \end{array} \right. \qquad (58b) $$

for $t \in [t_0, t_s]$.

## 4.3. Derivation of the time-varying terminal sets (find $Y_t$)

Noticing that with the ellipsoid $\mathcal{E}_t^1(\mathbf{x}_{e_t})$, we have not yet considered the inequality (48), the remainder term is bounded from above by: $\Phi_t(\mathbf{x}_{e_t}) \le \epsilon \mathbf{x}_{e_t}$, or from (44):

$$ g_{ha}(t, \mathbf{x}_{e_t}, K_t \mathbf{x}_{e_t}) - (A_t + B_t K_t) \mathbf{x}_{e_t} \le \epsilon \mathbf{x}_{e_t}. \qquad (59) $$

We now derive our time-varying terminal set ellipsoid $\mathcal{X}_{\text{MPC}}(t) = \mathcal{E}_t(\mathbf{x}_{e_t}) = \{\mathbf{x}_{e_t} \in \mathbb{R}^6 | \mathbf{x}_{e_t}^\top P_t \mathbf{x}_{e_t} \le Y_t\}$ as in (20) by reducing $Y_{1,t}$ to $Y_t$ such that the local stability condition (59) holds $\forall \mathbf{x}_{e_t} \in \mathcal{E}_t(\mathbf{x}_{e_t})$ with a chosen $\epsilon > 0$ and $(A_t, B_t)$ for $t \in [t_0, t_s]$. Therefore, all conditions 1(B)-1(H) in Assumption 1 are satisfied, and the stability of the closed-loop system is guaranteed in Theorem 1.

We present in Algorithm 1 the steps to design an NMPC scheme for trajectory tracking with stability ensured. We also show in Algorithm 2 the bisection algorithm to find the value of $Y_t$ as in (20).

# 5. Simulation results for the stable deep-stall and perching problem

## 5.1. Simulation parameters

For simulations, we use the Aerosonde FW UAV whose parameters are presented in Burston et al. (2014), Beard and McLain (2012, Appendix E.2). The aerodynamic coefficients of the FW UAV are plotted in
**Algorithm 1** NMPC design for trajectory tracking (adapted from Chen and Allgöwer (1998), Köhler et al. (2020a)).

**data** The dynamics $f_{ha}(\cdot)$ in (8), the reference trajectory $(\mathbf{x}_{r_t}, \mathbf{u}_{r_t})$ from solving (10)–(13), the time-varying error dynamics $g_{ha}(\cdot)$ as in (15)

1: Calculate $A_t = \frac{\partial f_{ha}}{\partial x}(\mathbf{x}_{r_t}, \mathbf{u}_{r_t}) \in \mathbb{R}^{6 \times 6}$, $B_t = \frac{\partial f_{ha}}{\partial u}(\mathbf{x}_{r_t}, \mathbf{u}_{r_t}) \in \mathbb{R}^{6 \times 2}$ as in (15)

2: Choose $Q \in \mathbb{S}_{++}^6$, $R \in \mathbb{S}_+^2$ in (19a)

3: Choose a positive $\epsilon > 0$ as in (48)

4: Solve the optimization problem (41) for $P_t$ and the feedback gain matrix $K_t$

5: Calculate $A_{K_t} = A_t + B_t K_t$ as in (44)

6: Find the largest $\Upsilon_{1_t} > 0$ by solving the optimization problem (58)

7: Construct the time-varying terminal invariant set $\mathcal{X}_{f_{MPC}}(t)$

$$ \mathcal{X}_{f_{MPC}}(t) = \{ \mathbf{x}_{e_t} \in \mathbb{R}^6 | \mathbf{x}_{e_t}^\top P_t \mathbf{x}_{e_t} \le \Upsilon_t \} \quad (60) $$

with $\Upsilon_t \in (0, \Upsilon_{1_t}]$ is a solution of an optimization problem:

$$ \Lambda = \max_{\mathbf{x}_{e_t}} \{ \mathbf{x}_{e_t}^\top P_t [ g_{ha}(t, \mathbf{x}_{e_t}, K_t \mathbf{x}_{e_t}) - A_{K_t} \mathbf{x}_{e_t} ] - \epsilon \mathbf{x}_{e_t}^\top P_t \mathbf{x}_{e_t} \}, \quad (61a) $$

$$ \text{s.t.} \quad \mathbf{u}_t = K_t \mathbf{x}_{e_t} + \mathbf{u}_{r_t} \in \mathcal{U}_{\text{MPC}} \text{ for all } \mathbf{x}_{e_t} \in \mathcal{X}_{f_{\text{MPC}}}(t), \quad (61b) $$

$$ \mathbf{x}_t = \mathbf{x}_{e_t} + \mathbf{x}_{r_t} \in \mathcal{X}_{\text{MPC}} \text{ for all } \mathbf{x}_{e_t} \in \mathcal{X}_{f_{\text{MPC}}}(t), \quad (61c) $$

$$ \mathbf{x}_{e_t}^\top P_t \mathbf{x}_{e_t} \le \Upsilon_t, \quad (61d) $$

$$ \Lambda \le 0. \quad (61e) $$

result The time-varying terminal sets $\mathcal{X}_{f_{MPC}}(t)$.

**Algorithm 2** Bisection algorithm to find $\Upsilon_t$ for the trajectory tracking problem (adapted from Köhler et al. (2020a), Müller and Allgöwer (2023)).

**data** the time-varying error dynamics $g_{ha}(\cdot)$ as in (15), $\mathbf{x}_{r_t}, \mathbf{u}_{r_t}, \epsilon, K_t, A_{K_t}, P_t$, and $\Upsilon_{1_t}$

1: **for** $k = 1$ : RefPoints **do**

2: $t_k \leftarrow t_0 + k\delta$

3: $\Upsilon_{1_t} \leftarrow \Upsilon_{1_{t_k}}$

4: $\bar{\Upsilon} \leftarrow \Upsilon_{1_{t_k}}$

5: $\underline{\Upsilon} \leftarrow 0$

6: $g(\mathbf{x}_e) \triangleq \mathbf{x}_e^\top P_{t_k} [ g_{ha}(t_k, \mathbf{x}_{e_{t_k}}, K_{t_k} \mathbf{x}_{e_{t_k}}) - A_{K_{t_k}} \mathbf{x}_e ] - \epsilon \mathbf{x}_e^\top P_{t_k} \mathbf{x}_e$

7: **for** $i = 1$ : Iteration **do**

8: $g_{\text{temp}} \leftarrow +\infty$

9: **for** $j = 1$ : RandomPoints **do**

10: Generate a random vector $\mathbf{x}_{\text{rand}}$ with $\|\mathbf{x}_{\text{rand}}\| = 1$

11: Formulate a temporary ellipsoid $\mathcal{E}_t^k(\mathbf{x}_e) = \{\mathbf{x}_e \in \mathbb{R}^{n_s} | \mathbf{x}_e^\top P_{t_k} \mathbf{x}_e \le \Upsilon_{1_{t_k}}\}$

12: Generate a test vector $\mathbf{x}_{\text{test}} = \sqrt{\Upsilon_{1_{t_k}}} P_{t_k}^{-1/2} \mathbf{x}_{\text{rand}}$, satisfying $\mathbf{x}_{\text{test}} \in \mathcal{E}_t^k$

13: Calculate $g_{\text{obj}} = -g(\mathbf{x}_{\text{test}})$

14: $g_{\text{temp}} = \min(g_{\text{temp}}, g_{\text{obj}})$

15: **end for**

16: $g_{\text{temp}} = -g_{\text{temp}}$

17: **if** $g_{\text{temp}} \ge 0$ **then**

18: $\bar{\Upsilon} \leftarrow \Upsilon_{1_{t_k}}$

19: **else**

20: $\underline{\Upsilon} \leftarrow \Upsilon_{1_{t_k}}$

21: **end if**

22: $\Upsilon_{1_{t_k}} \leftarrow (\bar{\Upsilon} + \underline{\Upsilon})/2$

23: **end for**

24: $\Upsilon_{1_{t_k}} \leftarrow \underline{\Upsilon}$

25: **end for**

26: **result** $\Upsilon_{1_{t_k}}$

<figure><img src="image_4.png" /><figcaption>Fig. 3. Aerodynamic coefficients.</figcaption></figure>

Fig. 3. The critical AoA ($\alpha_c$) is 24.07°. The vertical position of recovery net is $z_{\text{net}} \in [-4.7, -1.7]$ (m).

Regarding the generation task, the state and input constraint sets in the generation task are

$$ \mathcal{X}_g = \mathcal{X}_{g_1} \cup \mathcal{X}_{g_2}, \quad (62a) $$

$$ \mathcal{V}_g = \{ \mathbf{u} \in \mathbb{R}^2 | [-40^\circ, 0.1]^\top \le [\delta_e, \delta_t]^\top \le [40^\circ, 0.9]^\top \}, \quad (62b) $$

We choose the initial time $t_0$, the terminal time $t_f$, and the stability time $t_s$ as in (10)–(13) to be 0s, 24s, and 28s, respectively. In the stage cost (14), the weighting terms are chosen as $P_{\delta_t} = 1000$ and $P_{\Delta U} = 4000\delta_2$. We choose the airspeed for the steady-level flight at the beginning of the landing phase to be $V_{a_0} = 25$ (m/s) and the airspeed for the steady-level flight at the end of the stability phase is $V_{a_s} = 16$ (m/s). We solve the OCPs (10)–(11) with 121 initial positions $(x, -z)$ inside a square $[-290, -270]$ (m) $\times$ $[190, 210]$ (m).

Concerning the tracking task, the state and input constraint sets in the NMPC tracking task are

$$ \mathcal{X}_{\text{MPC}} = \mathcal{X}_{\text{MPC}_1} \cup \mathcal{X}_{\text{MPC}_2}, \quad (63a) $$

$$ \mathcal{V}_{\text{MPC}} = \{ \mathbf{u} \in \mathbb{R}^2 | [-60^\circ, 0]^\top \le [\delta_e, \delta_t]^\top \le [60^\circ, 1]^\top \}. \quad (63b) $$

We split the MPC state constraint set $\mathcal{X}_{\text{MPC}}$ into two sets: $\mathcal{X}_{\text{MPC}_1}$ before the net, and $\mathcal{X}_{\text{MPC}_2}$ after the net.

When $x < x_{\text{net}}$, we have the following constraints

$$ \mathcal{X}_{\text{MPC}_1} \triangleq \left\{ \begin{array}{l} \bar{x}_{\text{MPC}} \le x \le \bar{x}_{\text{MPC}} \\ \bar{z}_{\text{net}} \le -z \le a_{\bar{x}} x + \bar{z}_{\text{net}} \\ \bar{V}_{a_{\text{MPC}}} \le V_a \le a_{\bar{V}_a} x + \bar{V}_a \text{MPC} \\ \underline{\alpha}_{\text{MPC}} \le \alpha \le \bar{\alpha}_{\text{MPC}} \\ \underline{\theta}_{\text{MPC}} \le \theta \le \bar{\theta}_{\text{MPC}} \\ \underline{\omega}_y \le \omega_y \le \bar{\omega}_y \text{MPC} \end{array} \right. , \quad (64) $$

here, the variable $a_{\bar{V}_a}^x$ is the minimum (algebraic value) slope of the curve $V_a(x)$ before the net obtained from the generation task.
<table><caption>**Table 1**
Comparison of various bounds for the elevator and stabilator deflections in the literature.</caption><thead><tr><th>Work</th><th>Elevator/stabilator range</th><th>Deep-stall or Perching</th><th>2D or 3D</th></tr></thead><tbody><tr><td>Kim and Park (2022)</td><td>δ<sub>e</sub> = -40°</td><td>Deep-stall</td><td>3D experiments</td></tr><tr><td>Mathisen et al. (2021)</td><td>δ<sub>e</sub> ∈ [-40°, 40°]</td><td>Deep-stall</td><td>2D simulations</td></tr><tr><td>Cunis et al. (2020)</td><td>δ<sub>e</sub> = -60°</td><td>Deep-stall</td><td>2D and 3D analysis</td></tr><tr><td>Feroskhan et al. (2020)</td><td>δ<sub>e</sub> ∈ [-30°, 30°]</td><td>Perching</td><td>2D simulations</td></tr><tr><td>Venkateswara Rao et al. (2015)</td><td>δ<sub>e</sub> ∈ [-30°, 30°]</td><td>Perching</td><td>2D simulations</td></tr><tr><td>Alikhan and Go (2015)</td><td>δ<sub>e</sub> ∈ [-40°, 40°]</td><td>Perching</td><td>3D simulations</td></tr><tr><td>Pointner et al. (2011)</td><td>δ<sub>e</sub> = -45°</td><td>Deep-stall</td><td>2D simulations</td></tr><tr><td>Taniguchi (2008)</td><td>δ<sub>e</sub> = -55°</td><td>Deep-stall</td><td>2D simulations</td></tr><tr><td></td><td>δ<sub>s</sub> = -35°</td><td></td><td>3D experiments</td></tr><tr><td>Sim (1984)</td><td>δ<sub>s</sub> ∈ [-30°, -40°, -50°, -60°]</td><td>Deep-stall</td><td>3D experiments</td></tr></tbody></table>

<table><caption>**Table 2**
State and input constraints for the stable deep-stall and perching landing.</caption><thead><tr><th rowspan="2">Parameters</th><th colspan="4">Generation task</th><th colspan="2">NMPC trajectory tracking</th></tr><tr><th>t<sub>0</sub></th><th>t<sub>0</sub> &lt; t &lt; t<sub>f</sub></th><th>t<sub>f</sub></th><th>t<sub>f</sub> &lt; t &lt; t<sub>s</sub></th><th>t<sub>s</sub></th><th>t<sub>0</sub></th><th>t<sub>0</sub> &lt; t &le; t<sub>s</sub></th></tr></thead><tbody><tr><th>State constraint set</th><th>χ<sub>s1</sub></th><th>χ<sub>s1</sub> &subset; χ<sub>s2</sub></th><th>χ<sub>s2</sub></th><th>χ<sub>s2</sub> &subset; χ<sub>s3</sub></th><th>χ<sub>s3</sub></th><th>χ<sub>MPC</sub></th><th></th></tr><tr><th>x [m]</th><td>x<sub>0</sub></td><td>x<sub>0</sub> &le; x(t) &le; 0</td><td>0 &le; x(t<sub>f</sub>) &le; 10<sup>-4</sup></td><td>0 &le; x(t) &le; 100</td><td>59.5 &le; x(t<sub>s</sub>) &le; 60.5</td><td>x<sub>0</sub></td><td></td></tr><tr><th>z [m]</th><td>z<sub>0</sub></td><td>-250 &le; z(t) &le; -2</td><td>-4.2 &le; z(t<sub>f</sub>) &le; -2.2</td><td></td><td>z<sub>s</sub></td><td>z<sub>0</sub></td><td></td></tr><tr><th>V<sub>a</sub> [m/s]</th><td>V<sub>a0</sub></td><td>2 &le; V<sub>a</sub>(t) &le; 28</td><td>2 &le; V<sub>a0</sub>(t) &le; 7</td><td>2 &le; V<sub>a</sub>(t) &le; 28</td><td>V<sub>a0</sub></td><td>V<sub>a0</sub></td><td>χ<sub>MPC1</sub> &cup; χ<sub>MPC2</sub></td></tr><tr><th>α [deg]</th><td>α<sub>0</sub></td><td>0° &le; α(t) &le; 100°</td><td>90° &le; α(t<sub>f</sub>) &le; 100°</td><td>0° &le; α(t) &le; 100°</td><td>α<sub>s</sub></td><td>α<sub>0</sub></td><td></td></tr><tr><th>θ [deg]</th><td>θ<sub>0</sub></td><td>0° &le; θ(t) &le; 120°</td><td>0° &le; θ(t<sub>f</sub>) &le; 90°</td><td>0° &le; θ(t) &le; 120°</td><td>θ<sub>s</sub></td><td>θ<sub>0</sub></td><td></td></tr><tr><th>ω<sub>y</sub> [rad/s]</th><td>0</td><td>-1.5 &le; ω<sub>y</sub>(t) &le; 1.5</td><td></td><td></td><td>0</td><td>ω<sub>y0</sub></td><td></td></tr><tr><th>u(t) &isin; U</th><td>u<sub>0</sub></td><td>[-40°, 0.1]<sup>T</sup> &le; [δ<sub>e</sub>, δ<sub>s</sub>]<sup>T</sup> &le; [40°, 0.9]<sup>T</sup></td><td></td><td></td><td></td><td>[-60°, 0]<sup>T</sup> &le; [δ<sub>e</sub>, δ<sub>s</sub>]<sup>T</sup> &le; [60°, 1]<sup>T</sup></td><td></td></tr></tbody></table>

The constraints (64) are transformed into a polytope as

$$ \left[ \begin{matrix} -1 & 0 & 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 \\ -a_{-z_x} & -1 & 0 & 0 & 0 & 0 \\ 0 & 0 & -1 & 0 & 0 & 0 \\ -a_{-V_a^x} & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & -1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & -1 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 \end{matrix} \right] \begin{bmatrix} x \\ z \\ V_a \\ \alpha \\ \theta \\ \omega_y \end{bmatrix} \le \begin{bmatrix} -\frac{x_{-\text{MPC}}}{x_{\text{MPC}}} \\ -\frac{z_{-\text{net}}}{z_{\text{net}}} \\ \bar{z}_{\text{net}} \\ -\frac{V_a_{\text{MPC}}}{V_a_{\text{MPC}}} \\ \frac{-\alpha_{\text{MPC}}}{\alpha_{\text{MPC}}} \\ -\frac{\theta_{\text{MPC}}}{\theta_{\text{MPC}}} \\ -\frac{\omega_y}{\omega_y_{\text{MPC}}} \\ -\frac{\bar{\omega}_{y_{\text{MPC}}}}{\bar{\omega}_{y_{\text{MPC}}}} \end{bmatrix} \quad (65) $$

When $x \ge x_{\text{net}}$, we have the following constraints

$$ \chi_{\text{MPC2}} \triangleq \left\{ \begin{array}{l} \frac{x_{-\text{MPC}}^{+}}{x_{\text{MPC}}} \le x \le \bar{x}_{\text{MPC}}^{+} \\ \bar{z}_{\text{net}} \le -z \le \bar{z}_{\text{net}} \\ \frac{V_a}{V_{a_{\text{MPC}}}} \le V_a \le a_{-V_a^x}^{+} x^x + \bar{V}_{a_{\text{MPC}}} \\ \frac{\alpha}{\alpha_{\text{MPC}}} \le \alpha \le \bar{\alpha}_{\text{MPC}} \\ \theta_{\text{MPC}} \le \theta \le \bar{\theta}_{\text{MPC}} \\ \omega_y \le \omega_y \le \bar{\omega}_{y_{\text{MPC}}} \end{array} \right. , \quad (66) $$

here, the variable $a_{-V_a^x}^+$ is the maximum (algebraic value) slope of the curve $V_a(x)$ after the net obtained from the generation task.

The constraints (66) are transformed into another polytope as

$$ \left[ \begin{matrix} -1 & 0 & 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & -1 & 0 & 0 & 0 & 0 \\ 0 & 0 & -1 & 0 & 0 & 0 \\ -a_{-V_a^x}^+ & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & -1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & -1 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & -1 \\ 0 & 0 & 0 & 0 & 0 & 1 \end{matrix} \right] \begin{bmatrix} x \\ z \\ V_a \\ \alpha \\ \theta \\ \omega_y \end{bmatrix} \le \begin{bmatrix} -\frac{x_{-\text{MPC}}^{+}}{x_{\text{MPC}}} \\ -\frac{z_{-\text{net}}}{z_{\text{net}}} \\ \bar{z}_{\text{net}} \\ -\frac{V_a}{V_{a_{\text{MPC}}}} \\ \frac{-\alpha_{\text{MPC}}}{\alpha_{\text{MPC}}} \\ -\frac{\theta_{\text{MPC}}}{\theta_{\text{MPC}}} \\ -\frac{\omega_y}{\omega_y_{\text{MPC}}} \\ -\frac{\bar{\omega}_{y_{\text{MPC}}}}{\bar{\omega}_{y_{\text{MPC}}}} \end{bmatrix} \quad (67) $$

For the tracking task, we track the reference trajectory generated in Section 5.2 that starts from $(x, -z) = (-280, 200)$. Furthermore, to test the stability of the NMPC tracking scheme (17)–(18), we solve the tracking problem 100 times, with the initial error state of the FW UAV taking an arbitrary value from the set $\{x_e \in \mathbb{R}^6 | x_e - x_e \le x_e \le x_e + x_e\}$ with $x_e = [1, 1, 2, 5^\circ, 5^\circ, \pi/16]^\text{T}$, which is inside the MPC state constraint set discussed in Tables 2–3. The weighting matrices $Q$ and $R$ in (19a) are chosen to be constants as $Q = \text{diag}\{10, 10, 10, 10, 10, 10\}$ and $R = \text{diag}\{1, 1\}$. The upper bound $\epsilon$ as in (48) is chosen to be 0.1. All the state and input constraints are gathered in Table 2, with the specific bounds for the states in the NMPC tracking as in (64)–(66) are presented in Table 3. We choose the terminal $-z$ of the FW UAV in the generation task (at the beginning of the "imaginary" steady-level flight) to be the middle height of the net.

**Remark 2.** We list in Table 1 the bounds and minimum values for the elevator deflection angle $\delta_e$ and the stabilator (tail's horizontal stabilizer) deflection angle $\delta_s$ in post-stall maneuver-related works in the literature. For manned aerial vehicles, in the 1960s, NASA did some experiments on a modified SGS 1–36 sailplane with the minimum stabilator angle down to -60° (Sim, 1984). There exist works with a lower bound for the elevator deflection down to -40° (Alikhan & Go, 2015;
<table><caption>**Table 3**
Upper and lower bounds in the NMPC tracking for the stable deep-stall and perching landing.</caption><thead><tr><th>x̄<sub>MPC</sub> [m]</th><th>x̄<sup>+</sup><sub>MPC</sub> [m]</th><th>x̄<sup>-</sup><sub>MPC</sub> [m]</th><th>a<sup>-</sup><sub>z<sub>x</sub></sub></th><th>a<sup>+</sup><sub>v<sub>x</sub></sub></th><th>a<sup>-</sup><sub>v<sub>x</sub></sub></th></tr></thead><tbody><tr><td>x<sub>0</sub> - 50</td><td>4.0629</td><td>-5.6762</td><td>x<sub>f</sub> + 50</td><td>-1.1568</td><td>-0.6839</td><td>1.3213</td></tr><tr><td>z̄<sub>net</sub> [m]</td><td>z<sub>net</sub> [m]</td><td>V̄<sub>aMPC</sub> [m/s]</td><td>V<sub>aMPC</sub> [m/s]</td><td>θ̄<sub>MPC</sub> [deg]</td><td>θ̄<sub>MPC</sub> [deg]</td><td>ω<sub>γMPC</sub> [rad/s]</td></tr><tr><td>4.7</td><td>1.7</td><td>7.5</td><td>0.01</td><td>150°</td><td>-30°</td><td>π/2</td></tr></tbody></table>

Kim & Park, 2022; Mathisen et al., 2021), and the lowest is -60° in Cunis et al. (2020). This strengthens our choices for the input constraint in the generation task, as in (62b) and relaxing them in the tracking task as in (63b).

We solve the LMI (41) with the YALMIP toolbox (Lofberg, 2004) using the MOSEK solver (MOSEK, 2024) in MATLAB 2023a. The strict positive definite constraint $0 < X_{\min} \le X_t$ in (41c) is transformed into inequality constraints as $X_{\min} \le \mathbb{I}_6$ and $X_t \le \mathbb{I}_6$, which is due to the fact that $X_t \le \mathbb{I}_6 \Rightarrow X_t^{-1} \ge \mathbb{I}_6^{-1} \Rightarrow P_t \ge \mathbb{I}_6 > 0$. To avoid numerical errors, we also added a small positive definite matrix as $X_{\min} + 10^{-8}\mathbb{I}_6 \le X_t$. All the OCPs are transformed into NLPs, then are solved using the “Direct Multiple Shooting” method (Bock & Plitt, 1984), using the CasADi toolbox (Andersson et al., 2019) with solver IPOPT (Wächter & Biegler, 2006), in MATLAB R2023a. The simulations are carried out on a lab computer, with an AMD Ryzen 5 2600 6-core processor, 3.4GHz, 12 CPUs, and 16GB of RAM, running Windows 10. The time interval $[t_0, t_f]$ is divided into $N$ arcs, the time interval $[t_f, t_s]$ is split into $N'$ arcs, and the stability time $t_s$ is chosen such that

$$ \frac{t_f - t_0}{N} = \frac{t_s - t_f}{N'} = \delta. \qquad (68) $$

The equivalent time instants are

$$ \text{For } t \in [t_0, t_f]: \quad t_{k1} = t_0 + k_1 \delta, \text{ for } k_1 \in \{0, 1, \dots, N\}; \qquad (69a) $$

$$ \text{For } t \in [t_f, t_s]: \quad t_{k2} = t_f + k_2 \delta, \text{ for } k_2 \in \{0, 1, \dots, N'\}. \qquad (69b) $$

The dynamical model is discretized using the Runge-Kutta 4th-order method, with step $\delta = 0.1$.

The cost function $J_1(\mathbf{x}, \mathbf{u}) = \int_{t_0}^{t_f} \ell(\mathbf{x}(t), \mathbf{u}(t))dt$ as in (10) is discretized into a Riemann sum

$$ J_{1d} = \left( \sum_{k_1=0}^{N-1} \|\delta_t[k_1]\|_{P_{\delta t}}^2 + \sum_{k_1=0}^{N-2} \|\mathbf{u}[k_1+1] - \mathbf{u}[k_1]\|_{P_{\delta U}}^2 \right) \delta. \qquad (70) $$

with an additional input constraint regarding the steady-level flight

$$ \mathbf{u}[k_1 = 0] = \mathbf{u}_0, \qquad (71) $$

where $\mathbf{u}_0$ is the steady-level flight input before the landing phase. In addition, the cost function $J_2(\mathbf{x}, \mathbf{u}) = \int_{t_f}^{t_s} \ell(\mathbf{x}(t), \mathbf{u}(t))dt$ as in (12) is discretized into a Riemann sum

$$ J_{2d} = \left( \sum_{k_2=0}^{N'-1} \|\delta_t[k_2]\|_{P_{\delta t}}^2 + \sum_{k_2=0}^{N'-2} \|\mathbf{u}[k_2+1] - \mathbf{u}[k_2]\|_{P_{\delta U}}^2 \right) \delta, \qquad (72) $$

with an additional input constraint

$$ \mathbf{u}[k_2 = 0] = \mathbf{u}_r[N-1], \qquad (73) $$

where $\mathbf{u}_r[N-1]$ is the optimal solution obtained by solving the OCP (10)-(11) in the landing phase.

Regarding the prediction horizon $T_p$ in (17), we start with a two-step (0.2s) prediction horizon, and incrementally increase it until the first NMPC iteration gives a solution (feasible). Since the NMPC scheme has the recursive feasibility property, other NMPC iterations afterward also give solutions with state and input constraints guarantees.

## 5.2. Results of the generation task for the stable deep-stall and perching problem

Out of the 121 initial positions $(x, -z)$ inside a square $[-290, -270] \times [190, 210]$ $(m)$, only 5 initial positions give infeasible solutions, which are on the top left corner and could be seen in Fig. 4.

Fig. 5 shows zoomed-in segments of the generated trajectories with $x \in [-50, 5]$ $(m)$, focusing on the region before the net, while Fig. 6 gathers the zoomed-in segments for $x \in [-10, 60]$ $(m)$, concentrating on the region behind the net. All trajectories reach the upper bound of $-z$ in the set $X_{g2}$, and the FW UAV eventually reaches steady-level flight around $x = 59.5$ $(m)$.

The bundle of generated airspeed $V_a$, angle of attack $\alpha$, pitch angle $\theta$, and pitch rate $\omega_y$ are presented in Figs. 7–11, with a zoom-in view for the angle of attack for $x \in [-20, 20]$ $(m)$ in Fig. 9. Figs. 12–13 collect the bundle of elevator deflection $\delta_e$ and pulse-width modulation $\delta_t$. There exists a stable deep-stall phase with constant states and inputs for all trajectories. All the generated states and inputs satisfy the state and input constraint sets (2).

We store the values of the matrices $P_t$, $K_t$, $X_t$, $X_{\min}$, $Y_t$, and the time-varying variables $Y_{1t}$ and $Y_t$ online at https://doi.org/10.6084/m9.fishshare.27186912.v3.

Fig. 14 shows the base-10 logarithmic scale of $Y_{1t}$ of all the ellipsoids $E_t^1(x_{e_t}) = \{x_e \in \mathbb{R}^6 | x_{e_t}^T P_t x_{e_t} \le Y_{1t}\}$ as in (55), satisfying the polytopic constraints (58b), along the generated optimal landing trajectory.

Next, we apply Algorithm 2 to find $Y_t$ and plot its base-10 logarithmic scale along the reference trajectory in Fig. 15.

## 5.3. Results of the tracking task for the stable deep-stall and perching problem

We present the NMPC tracking simulation results in Figs. 16–31.

Fig. 16 gathers all the prediction horizon values $T_p$ of the 100 tracking simulations, which are obtained by incrementally increasing the prediction horizon from 2 steps until the first NMPC iteration is solvable. Fig. 17 depicts the mean solving time of each NMPC iteration for all 100 scenarios, with the grand mean solving time of 1.1105s.

Figs. 18–25 focus on the trajectory tracking results. We choose not to plot the terminal ellipsoids $E_t(x_{e_t})$ as in (20) along the simulated trajectories for two reasons: (i) they are merely fictitious sets without any true physical significance, and (ii) they are too small to observe. All the trajectories eventually converge to the reference trajectory at the horizontal position $x$ around $-220$ $(m)$.

We present the zoomed-in segments of the tracking results for $x \le -200$ $(m)$ in Figs. 26–31. It is clear that by tightening the input constraints in the generation task and relaxing them in the tracking task, the FW UAV has control authority to track the nominal reference trajectory, although starting from random initial states inside the state constraint set.

To visualize the tracking performance more and to quantify it, we introduce the absolute error and the root mean square (RMS) error. The absolute errors for 100 tracking cases at each instant $t_k = t_0 + k\delta$, for $k \in \{0, 1, \dots, 280\}$, are calculated by

$$ iAE[k] = \|\mathbf{x}[k] - \mathbf{x}_r[k]\|, \quad \text{for } i \in \{1, 2, \dots, 100\}. \qquad (74) $$
<figure><img src="image_5.png" /><figcaption>Fig. 4. The bundle of generated trajectories with feasible and infeasible initial positions.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. The bundle of generated trajectories before the net.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6. The bundle of generated trajectories after the net.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7. The bundle of generated airspeed $V_a$.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Fig. 8. The bundle of the generated angle of attack $\alpha$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9. The zoomed-in view of the generated bundle angle of attack $\alpha$ for $x \in [-20, 20] (m)$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10. The bundle of generated pitch angle $\theta$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 11. The bundle of generated pitch rate $\omega_y$.</figcaption></figure>

The RMS errors of all 100 tracking scenarios are calculated using the equation

$$ \dot{i}\text{RMSE} = \sqrt{\frac{1}{281} \sum_{k=0}^{280} \|\mathbf{x}[k] - \mathbf{x}_r[k]\|^2}, \quad \text{for } i \in \{1, 2, \dots, 100\}. \quad (75) $$

The absolute errors for the 100 tracking cases, as in (74), are presented in Fig. 32. We can see that in all cases, the absolute errors become less than or equal to 0.01 when $x$ surpasses -238.9233 (m). The 1% settling time (i.e., when $\dot{i}AE \le 0.01$) is gathered in Fig. 33, with the mean settling time being 2.5030 (s). Fig. 34 depicts the RMS errors as in (75) for 100 scenarios, with the mean RMS error equal to 0.1548.
<figure><img src="image_13.png" /><figcaption>Fig. 12. The bundle of generated elevator deflection $\delta_e$.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 13. The bundle of generated pulse-width modulation $\delta_t$.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 14. The value of $\Upsilon_t$, along the reference trajectory.</figcaption></figure>

The detailed simulation results are stored online at https://doi.org/10.6084/m9 figshare.30258457.v1.

**Remark 3.** We emphasize that the linearization (15) is around the reference trajectory, which includes “abnormal” states, for example, large AoA (up to 100°) and large pitch (up to more than 90°), which are not used in commercial aviation. However, these values are just nominal states of the FW UAV in the sense of mathematics and dynamics while executing the deep-stall and perching maneuver. The error dynamics around these states (which belong to the reference trajectory) are still small enough to apply the NMPC tracking controller (17)–(18).

**Remark 4.** For dynamic missions, we identify two possible problems. First, when the sampling time might be crucial, we can employ the *mesh refinement* technique (Paiva & Fontes, 2019) to adaptively change the sampling time, i.e., shorter before and after the net to capture the nonlinear dynamics of the FW UAV in perching, and longer in the stable deep-stall phase to reduce the computational effort. Second, if the FW UAV has different initial states from the bundle of generated reference trajectories, or there is model mismatch, one should not generate new trajectories online with the provided NMPC tracking scheme in this

paper, because this entails finding new $P_t$ and $\Upsilon_t$. Instead, the recommended approach is to increase the prediction horizon to achieve the recursive feasibility property.

**Remark 5.** Concerning the real-life implementation of the proposed NMPC tracking scheme, we propose the following refinements: (i) the computational time can be reduced for embedded hardware by using the C code generation feature of CasADi (Andersson et al., 2019), or the acados toolbox (Verschueren et al., 2022); (ii) although we do not consider explicitly the robustness aspect of the tracking scheme in this paper, the NMPC controller has the inherent robustness property (Magni & Sepulchre, 1997) which can handle the uncertainties that may arise. Moreover, if the FW UAV encounters disturbances in its practical operation, one should not generate a new reference trajectory and solve for the terminal ingredients of the NMPC tracking scheme again, but just increase the prediction horizon such that the NMPC iteration gives a feasible solution; (iii) regarding the constraint tightening in the generation task, interested readers could consult the work of Köhler et al. (2020b, Appendix B), where the authors present a systematic approach
<figure><img src="image_16.png" /><figcaption>Fig. 15. The value of $Y_t$ along the reference trajectory.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 16. The prediction horizons in each scenario.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 17. The mean calculating time for an NMPC iteration in each scenario.</figcaption></figure>

to selecting the constraint tightening with bounded disturbances for nonlinear systems in the discrete-time domain.

### 5.4. The relations between wing loading, airspeed, elevator deflection angle, and flight path angle in steady flight

This section investigates the relations between wing loading $W_L$, airspeed $V_a$, elevator deflection angle $\delta_e$, and flight path angle $\gamma$ in steady flight.

The wing loading $W_L$ is defined by the mass of the UAV divided by its wing's area as

$$W_L = \frac{m}{bc} \qquad (76)$$

Low wing loading increases the aircraft's agility, reduces turning radius in a coordinated turn (Phillips, 2004, Section 3.9), and induces low landing speed (Phillips, 2004, Section 3.8). Aircraft with low wing loading also have high maneuverability and require a shorter landing distance (Brandt et al., 2004). However, low wing loading aircraft will be

more susceptible to disturbances (Nelson, 1998, Chapter 6). This could be mitigated with a low aspect ratio (Anderson, 2016, Chapter 5).

We gather in Table 4 the physical parameters (mass, wing area, wing loading, and aspect ratio) of the aircraft used in post-stall maneuver-related works in the literature. The majority of the aircraft have their wing loading less than or equal to 10 (kg/m²).

Now, considering the dynamics in steady flight with only four states $\mathbf{x}_{st} \triangleq [V_a, \alpha, \theta, \omega_y]^\top$ and two inputs $\mathbf{u}_{st} \triangleq [\delta_e, \delta_t]^\top$, satisfying the ODE $\dot{\mathbf{x}}_{st} = f_{st}(\mathbf{x}_{st}, \mathbf{u}_{st})$, where $f_{st}(\cdot)$ gathers the right-hand side of the following dynamics

$$\dot{V}_a = -F_D/m - g \sin(\gamma) + F_p \cos(\alpha)/m, \qquad (77a)$$

$$\dot{\alpha} = \omega_y - \frac{1}{mV_a} [F_L/m - g \cos(\gamma) + F_p \sin(\alpha)/m], \qquad (77b)$$

$$\dot{\theta} = \omega_y, \qquad (77c)$$

$$\dot{\omega}_y = M/I_{yy}. \qquad (77d)$$

First, we want to know the relations between the wing loading, airspeed, and elevator deflection angle in steady-level flight. To observe the
<figure><img src="image_19.png" /><figcaption>Fig. 18. The bundle of NMPC from various initial positions.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 19. The zoomed-in view of the positions (x, -z) before and after the net.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Fig. 20. The bundle of airspeed $V_a$ from NMPC tracking.</figcaption></figure>

effect of wing loading on steady-level flight airspeed, we change the size of the wing by modifying the wingspan b and the chord c simultaneously, while maintaining the wing aspect ratio $AR = b^2/S = b/c$:

$$ \frac{m}{b'c'} = \epsilon_{W_L} \frac{m}{bc} \quad \text{and} \quad \frac{b'}{c'} = \frac{b}{c}. \qquad (78) $$

where $\epsilon_{W_L}$ is a constant that will be varied to changed the aspect ratio, $b'$ and $c'$ are the modified wingspan and chord. The aspect ratio directly influences the induced drag as in (5b), so we want to isolate the effect of the aspect ratio and only focus on the effect of wing loading in our analysis.

By solving (78), we obtain a way to modify the shape of the wing for our test

$$ b' = \frac{b}{\sqrt{\epsilon_{W_L}}}, \quad c' = \frac{c}{\sqrt{\epsilon_{W_L}}}. \qquad (79) $$

Then, we solve an optimization problem to get the minimum airspeed in steady-level flight:

$$ V_a = \arg \min_{V_a} V_a^2, \qquad (80) $$

$$ \text{subject to: } \dot{\mathbf{x}}_{st} = f_{st}(\mathbf{x}_{st}, \mathbf{u}_{st}), \qquad (81a) $$

$$ \dot{V}_a = 0, \dot{\alpha} = 0, \dot{\omega}_y = 0, \gamma = 0, \omega_y = 0, W_L = W_{L,SL}, \qquad (81b) $$

$$ -10^{\circ} \le \alpha \le 110^{\circ}, \quad 0.01 \le V_a \le 30, \qquad (81c) $$
<figure><img src="image_22.png" /><figcaption>Fig. 21. The bundle of angle of attack $\alpha$ from NMPC tracking.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Fig. 22. The bundle of pitch angle $\theta$ from NMPC tracking.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Fig. 23. The bundle of pitch rate $\omega_y$ from NMPC tracking.</figcaption></figure>

$$ \delta_e \le \delta_e \le -\delta_e, \quad 0 \le \delta_t \le 1. \qquad (81d) $$

We vary the lower bounds on the elevator $\delta_e \in \{-30^\circ, -40^\circ, -50^\circ, -60^\circ\}$. The initial guesses for the states and inputs are $\mathbf{x}_{\text{guess}_1} = [10, 0.001, 0.001, 0]^\mathrm{T}$ and $\mathbf{u}_{\text{guess}_1} = [\delta_e, 0.35]^\mathrm{T}$. The wing loading $\mathbf{W}_{L,SL}$ is chosen by sweeping various values of $\epsilon_{W_L}$ in $\{0.1, 0.2, \dots, 0.9, 1\}$.

The steady-level flight states and inputs are obtained by solving the set of equations: $V_a = V_{a_r}, \omega_y = 0, \gamma = 0, \dot{V}_a = 0, \dot{\alpha} = 0$, for $\delta_e, \delta_t, \theta$, and $\alpha$.

We gather in Tables 5 and 6 the wing loading, states, and control inputs of the FW UAV obtained by solving the optimization problem (80)–(81), for $\delta_e = -40^\circ$, and $\delta_e = -60^\circ$, respectively. For the same FW UAV, the lower the wing loading, the lower elevator bound $\delta_e = -60^\circ$

gives smaller steady-level airspeed at a higher pitch angle (and AoA), while consuming more energy (larger input PWM $\delta_t$), compared to the lower elevator bound $\delta_e = -40^\circ$.

We also present in Fig. 35 the optimal airspeed $V_a$ obtained from solving (80)–(81) with respect to various wing loading values $W_L$ and lower bound of the elevator deflections $\delta_e$.

From Fig. 35, two immediate comments are made. First, with the same lower bound on elevator deflection, the smaller the wing loading, the smaller the steady-level airspeed. Second, at the same wing loading, the smaller the lower bound on the elevator deflection, the smaller the steady-level airspeed that the FW UAV achieves.

Based on the mathematical aspect of the steady-level flight, we infer that the steady deep-stall state is just another set of states and in-
<figure><img src="image_25.png" /><figcaption>Fig. 24. The bundle of elevator deflection $\delta_e$ from NMPC tracking.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Fig. 25. The bundle of thrust PWM $\delta_t$ from NMPC tracking.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>Fig. 26. The zoomed-in view of the airspeed bundle $V_a$ from NMPC tracking.</figcaption></figure>

puts satisfying a set of equations. Therefore, in a second test, we want to know the relation between the airspeed, wing loading, and elevator deflection angle in deep-stall, by solving an optimization problem to get the minimum airspeed in steady deep-stall when the elevator is in $-40^{\circ} \le \delta_e \le -31^{\circ}$:

$$V_a = \arg \min_{V_a} V_a^2, \qquad (82)$$

$$\text{subject to: } \dot{\mathbf{x}}_{st} = f_{st}(\mathbf{x}_{st}, \mathbf{u}_{st}), \qquad (83a)$$

$$\dot{V}_a = 0, \dot{\alpha} = 0, \dot{\omega}_y = 0, \gamma = \gamma_{st}, \omega_y = 0, \qquad (83b)$$

$$-10^{\circ} \le \alpha \le 110^{\circ}, \qquad (83c)$$

$$-40^{\circ} \le \delta_e \le -31^{\circ}, 0 \le \delta_t \le 1. \qquad (83d)$$

with $\gamma_{st} \in \{-50^{\circ}, -45^{\circ}, \dots, 0^{\circ}\}$, and $\epsilon_{W_L}$ in $\{0.25, 0.5, 0.75, 1\}$. The guesses for the states and inputs are $\mathbf{x}_{\text{guess}_2} = [10, 0.001, 0.001, 0]^\top$ and $\mathbf{u}_{\text{guess}_1} = [-40^{\circ}, 0.35]^\top$.

We show in Fig. 36 the optimal airspeed as the solution of the optimization problem (82)-(83) with various wing loading and lower elevator bounds. For the same wing loading, the steeper the deep-stall landing (the smaller the flight path angle), the larger the airspeed. Moreover, the smaller the wing loading, the smaller the airspeed in a stable deep-stall.

At $\gamma_{st} = -50^{\circ}$, the wing loading $W_L = 6.14$ and $W_L = 12.27$ could not give a feasible solution at $\delta_e$ around $-40^{\circ}$. In other words, with the lower elevator bound $\delta_e = -40^{\circ}$, the wing loading $W_L = 6.14$ and $W_L = 12.27$, the FW UAV could not execute a deep-stall maneuver with a steady flight path angle $\gamma_{st} = -50^{\circ}$.
<figure><img src="image_28.png" /><figcaption>Fig. 27. The zoomed-in view of the angle of attack bundle $\alpha$ from NMPC tracking.</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>Fig. 28. The zoomed-in view of the pitch angle bundle $\theta$ from NMPC tracking.</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>Fig. 29. The zoomed-in view of the pitch rate bundle $\omega_y$ from NMPC tracking.</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>Fig. 30. The zoomed-in view of the elevator deflection bundle $\delta_e$ from NMPC tracking.</figcaption></figure>
<figure><img src="image_32.png" /><figcaption>Fig. 31. The zoomed-in view of the thrust PWM bundle $\delta_t$ from NMPC tracking.</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>Fig. 32. The absolute errors from NMPC tracking.</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>Fig. 33. The 1% settling time from NMPC tracking.</figcaption></figure>

<figure><img src="image_35.png" /><figcaption>Fig. 34. The RMS errors from NMPC tracking.</figcaption></figure>
<table><caption>**Table 4**
Comparison of wing loading for different UAVs in the literature.</caption><thead><tr><th>Work</th><th>UAV mass [kg]</th><th>Wing area [m²]</th><th>Wing loading [kg/m²]</th><th>AR [-]</th></tr></thead><tbody><tr><td rowspan="2">Basescu et al. (2023)</td><td>3.94</td><td>0.4630</td><td>8.5097</td><td>5.0164</td></tr><tr><td>4.24</td><td>0.4017</td><td>10.5551</td><td>5.7819</td></tr><tr><td>Kim and Park (2022)</td><td>2.5</td><td>0.705</td><td>3.546</td><td>14.524</td></tr><tr><td>Nguyen et al. (2023), Mathisen et al. (2021)</td><td>13.5</td><td>0.55</td><td>24.55</td><td>15.24</td></tr><tr><td>Cunis et al. (2020)</td><td>1.55</td><td>0.277</td><td>5.596</td><td>9.540</td></tr><tr><td>Feroskhan et al. (2020)</td><td>0.9</td><td>0.3</td><td>3</td><td>4.8</td></tr><tr><td>Venkateswara Rao and Yanhua (2019)</td><td>1.2</td><td>0.25</td><td>4.8</td><td>4</td></tr><tr><td>Venkateswara Rao and Go (2016), Alikhan and Go (2015)</td><td>0.0635</td><td>0.0347</td><td>1.8300</td><td>5</td></tr><tr><td>Venkateswara Rao et al. (2015), Venkateswara Rao and Go (2014)</td><td>0.8</td><td>0.25</td><td>3.2</td><td>4</td></tr><tr><td>Puopolo and Jacob (2015)</td><td>0.88</td><td>0.2961</td><td>2.9720</td><td>6.7142</td></tr><tr><td>Pointner et al. (2011)</td><td>1.5</td><td>0.28</td><td>5.3571</td><td>4.1935</td></tr><tr><td>Taniguchi (2008)</td><td>0.231</td><td>0.159</td><td>1.453</td><td>12</td></tr></tbody></table>

<table><caption>**Table 5**
Varying steady-level flight with $\delta_e \in [-40°, 40°]$.</caption><thead><tr><th>ε<sub>W<sub>L</sub></sub></th><th>W<sub>L</sub></th><th>b</th><th>c</th><th>V<sub>a</sub></th><th>α [deg]</th><th>θ [deg]</th><th>δ<sub>e</sub> [deg]</th><th>δ<sub>t</sub></th></tr></thead><tbody><tr><td>0.1</td><td>2.45</td><td>9.16</td><td>0.60</td><td>3.85</td><td>49.11</td><td>49.11</td><td>-40</td><td>0.33</td></tr><tr><td>0.2</td><td>4.91</td><td>6.47</td><td>0.42</td><td>5.45</td><td>49.11</td><td>49.11</td><td>-40</td><td>0.34</td></tr><tr><td>0.3</td><td>7.36</td><td>5.29</td><td>0.35</td><td>6.67</td><td>49.11</td><td>49.11</td><td>-40</td><td>0.34</td></tr><tr><td>0.4</td><td>9.82</td><td>4.58</td><td>0.30</td><td>7.70</td><td>49.11</td><td>49.11</td><td>-40</td><td>0.34</td></tr><tr><td>0.5</td><td>12.27</td><td>4.09</td><td>0.27</td><td>8.61</td><td>49.11</td><td>49.11</td><td>-40</td><td>0.35</td></tr><tr><td>0.6</td><td>14.73</td><td>3.74</td><td>0.25</td><td>9.43</td><td>49.11</td><td>49.11</td><td>-40</td><td>0.35</td></tr><tr><td>0.7</td><td>17.18</td><td>3.46</td><td>0.23</td><td>10.19</td><td>49.11</td><td>49.11</td><td>-40</td><td>0.35</td></tr><tr><td>0.8</td><td>19.64</td><td>3.24</td><td>0.21</td><td>10.89</td><td>49.11</td><td>49.11</td><td>-40</td><td>0.36</td></tr><tr><td>0.9</td><td>22.09</td><td>3.05</td><td>0.20</td><td>11.55</td><td>49.11</td><td>49.11</td><td>-40</td><td>0.36</td></tr><tr><td>1</td><td>24.55</td><td>2.9</td><td>0.19</td><td>12.18</td><td>49.11</td><td>49.11</td><td>-40</td><td>0.36</td></tr></tbody></table>

<table><caption>**Table 6**
Varying steady-level flight with $\delta_e \in [-60°, 60°]$.</caption><thead><tr><th>ε<sub>W<sub>L</sub></sub></th><th>W<sub>L</sub></th><th>b</th><th>c</th><th>V<sub>a</sub></th><th>α [deg]</th><th>θ [deg]</th><th>δ<sub>e</sub> [deg]</th><th>δ<sub>t</sub></th></tr></thead><tbody><tr><td>0.1</td><td>2.45</td><td>9.16</td><td>0.60</td><td>2.17</td><td>75.42</td><td>75.42</td><td>-60</td><td>0.39</td></tr><tr><td>0.2</td><td>4.91</td><td>6.47</td><td>0.42</td><td>3.07</td><td>75.42</td><td>75.42</td><td>-60</td><td>0.39</td></tr><tr><td>0.3</td><td>7.36</td><td>5.29</td><td>0.35</td><td>3.76</td><td>75.42</td><td>75.42</td><td>-60</td><td>0.39</td></tr><tr><td>0.4</td><td>9.82</td><td>4.58</td><td>0.30</td><td>4.34</td><td>75.42</td><td>75.42</td><td>-60</td><td>0.39</td></tr><tr><td>0.5</td><td>12.27</td><td>4.09</td><td>0.27</td><td>4.85</td><td>75.42</td><td>75.42</td><td>-60</td><td>0.39</td></tr><tr><td>0.6</td><td>14.73</td><td>3.74</td><td>0.25</td><td>5.31</td><td>75.42</td><td>75.42</td><td>-60</td><td>0.39</td></tr><tr><td>0.7</td><td>17.18</td><td>3.46</td><td>0.23</td><td>5.74</td><td>75.42</td><td>75.42</td><td>-60</td><td>0.39</td></tr><tr><td>0.8</td><td>19.64</td><td>3.24</td><td>0.21</td><td>6.14</td><td>75.42</td><td>75.42</td><td>-60</td><td>0.39</td></tr><tr><td>0.9</td><td>22.09</td><td>3.05</td><td>0.20</td><td>6.51</td><td>75.42</td><td>75.42</td><td>-60</td><td>0.39</td></tr><tr><td>1</td><td>24.55</td><td>2.9</td><td>0.19</td><td>6.86</td><td>75.42</td><td>75.42</td><td>-60</td><td>0.39</td></tr></tbody></table>

<figure><img src="image_36.png" /><figcaption>**Fig. 35.** Steady airspeed $V_a$ with various wing loading values $W_L$ and lower bounds for the elevator $\delta_e$.</figcaption></figure>
<figure><img src="image_37.png" /><figcaption>Fig. 36. Steady airspeed $V_a$ with various flight path angle $\gamma$ and wing loading values $W_L$, at $\delta_e = -40^\circ$.</figcaption></figure>

<figure><img src="image_38.png" /><figcaption>Fig. 37. Quantitative tracking results with varying upper remainder bound $\epsilon$.</figcaption></figure>
From this analysis, we conclude that compared to just belly landing a FW UAV in deep-stall, for example, as in Mathisen et al. (2021), Taniguchi (2008), our deep-stall and perching maneuver can reduce the terminal airspeed of the FW UAV when it touches the recovery net close to the ground. Moreover, for different “desired” steady deep-stall flight path angle in the generation task, we will have different values of reference airspeed, which lead to changes in the state constraints $\mathcal{X}_{\text{MPC}}$. More specifically, different values of steady flight path angle cause changes in the term $a_{zx}^-$ in the state constraint $z_{\text{net}}^- \le -z \le a_{zx}^-$ and $\bar{z}_{\text{net}}^-$ as in (64); different airspeed profiles yield changes in the term $a_{V_{ax}}^-$ in the state constraint $V_{a_{\text{MPC}}} \le V_a \le a_{V_{ax}}^-$ and $\bar{V}_{a_{\text{MPC}}}^-$ as in (64), and in the term $a_{V_{ax}}^+$ in the state constraint $V_{a_{\text{MPC}}} \le V_a \le a_{V_{ax}}^+$ and $\bar{V}_{a_{\text{MPC}}}^-$ as in (66). These changes also induce the change in the terminal set $\mathcal{X}_{f_{\text{MPC}}}$ as in (20), following the changes in the optimization problem (58).

### 5.5. Tracking performance with different remainder bounds

To examine the tracking performance with different remainder bounds, we run the simulations with the same 100 initial conditions, with additional various values of $\epsilon$, more specifically, $\epsilon = 0.01$, $\epsilon = 1$, and $\epsilon = 10$, then gather the results in Fig. 37. We see that the difference in RMS errors is insignificant (Fig. 37b). Going from $\epsilon = 0.01$ to $\epsilon = 1$, the prediction horizons decrease (Fig. 37a), which leads to a decrease in the mean calculating time (Fig. 37c). The tradeoff of this is the increase in the 1% settling time in Fig. 37d. When $\epsilon = 10$, the mean 1% settling time is shortest (Fig. 37d); however, there are four cases when the mean calculating time is more excessive (Fig. 37c), more than 3.5s for each NMPC iteration. The detailed box plots can be viewed online at https://doi.org/10.6084/m9.fishshare.30258340.v2. There, we also have a figure depicting the volume of the terminal ellipsoids and the time-varying variable $Y_i$ for each $\epsilon$.

## 6. Conclusions and future work

We successfully demonstrate a deep-stall and perching landing maneuver with stability ensured for a fixed-wing UAV. This is done in two steps: first, generating an optimal landing trajectory and then tracking it from various random initial states. The landing trajectory is found by solving two consecutive OCPs. The asymptotic tracking task is then handled by an NMPC scheme with time-varying terminal costs and the time-varying ellipsoid sets as terminal sets. We also analyze the effect of wing loading on steady airspeed and flight path angle in level flight and deep-stall. We conclude that the terminal perching airspeed is strongly dependent on the wing loading of the FW UAV.

Future work will revolve around the real-life implementation of this maneuver, namely: i) tackling the robustness problem to find the relations between the bounded disturbances and the aforementioned tightening; ii) designing more efficient computational implementations for online tracking control by enlarging the terminal ellipsoids, for example, by tuning the gains of a feedback linearization controller which serves as the local controller inside the terminal set (Nguyen et al., 2024); iii) applying the stable deep-stall and perching landing on different FW UAVs with various airframes (mass, wing loading, chord, wing span, aspect ratio, wing shape), to validate the feasibility, versatility, and effectiveness of the technique.

### CRediT authorship contribution statement

**Huu Thien Nguyen:** Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Formal analysis, Data curation, Project administration, Funding acquisition, Conceptualization; **Ionela Prodan:** Writing – review & editing, Supervision, Project administration, Funding acquisition; **Fernando A.C.C. Fontes:** Writing – review & editing, Validation, Supervision, Resources, Project administration, Funding acquisition, Conceptualization.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgments

The authors are grateful for the contributions of Professor Fernando Lobo Pereira at the initial stages of this work. This research is supported by FCT/MCTES(PIDDAC), through projects 2020.07959.BD (https://doi.org/10.54499/2020.07959.BD), 2022.02320.PTDC-KEFCODE (https://doi.org/10.54499/2022.02320.PTDC), and 2022.02801.PTDC-UPWIND-ATOL (https://doi.org/10.54499/2022.02801.PTDC). The work of I. Prodan has benefited from a French government grant managed by the Agence Nationale de la Recherche under the France 2030 program, reference ANR-23-IACL-0006; the La Région Auvergne-RhôneAlpes, Pack Ambition Recherche 2021 - PlanMAV, reference RECPLAMALCIR; International Research Booster Univ. Grenoble Alpes, reference C900IDX24V30.

### References

* AeroEnvironment Inc. (2021). Raven® B Group 1 UAS | Surveillance and Reconnaissance Drone. https://www.avinc.com/uaas/raven.

* Alessandretti, A., Aguiar, A. P., & Jones, C. N. (2013). Trajectory-tracking and path-following controllers for constrained underactuated vehicles using model predictive control. In *2013 European control conference (ECC)* (pp. 1371–1376). Zurich: IEEE. https://doi.org/10.23919/ECC.2013.6669717

* Alikhan, M., & Go, T. H. (2015). Aerodynamic modeling and optimization of sideslip perching maneuver. In *AIAA Atmospheric flight mechanics conference*. Kissimmee, Florida: American Institute of Aeronautics and Astronautics. https://doi.org/10.2514/6.2015-0235

* Alikhan, M., Peyada, N. K., & Go, T. H. (2013). Flight dynamics and optimization of three-dimensional perching maneuver. *Journal of Guidance, Control, and Dynamics*, 36(6), 1791–1797. https://doi.org/10.2514/1.58894

* Anderson, J. D. (2016). Introduction to Flight. (Eighth ed.). New York, NY: McGraw-Hill Education.

* Andersson, J. A. E., Gillis, J., Horn, G., Rawlings, J. B., & Diehl, M. (2019). CasADI: A software framework for nonlinear optimization and optimal control. *Mathematical Programming Computation*, 11(1), 1–36. https://doi.org/10.1007/s12532-018-0139-4

* Aslam, M. S., Bilal, H., & Vasilakos, A. V. (2025). Self-triggered scheme design for Takagi-Sugeno fuzzy model based on mismatch premise variable with time-varying delay. *IEEE Transactions on Automation Science and Engineering*, 22, 15536–15548. https://doi.org/10.1109/TASE.2025.3570089

* Aslam, M. S., Shamrooz, S., & Bilal, H. (2024). Fuzzy PD-sliding mode control design for networked system with time delays. *European Journal of Control*, 78, 101020. https://doi.org/10.1016/j.ejcon.2024.101020

* Basescu, M., & Moore, J. (2020). Direct NMPC for post-stall motion planning with fixed-wing UAVs. In *2020 IEEE International conference on robotics and automation (ICRA)* (pp. 9592–9598). Paris, France: IEEE. https://doi.org/10.1109/ICRA40945.2020.9196724

* Basescu, M., Yeh, B., Scheuer, L., Wolfe, K., & Moore, J. (2023). Precision post-stall landing using NMPC with learned aerodynamics. *IEEE Robotics and Automation Letters*, 8(5), 3031–3038. https://doi.org/10.1109/LRA.2023.3264738

* Beard, R. W., & McLain, T. W. (2012). Small Unmanned Aircraft: Theory and Practice. Princeton, N.J: Princeton University Press.

* Bock, H. G., & Plitt, K. J. (1984). A multiple shooting algorithm for direct solution of optimal control problems. *IFAC Proceedings Volumes*, 17(2), 1603–1608. https://doi.org/10.1016/S1474-6670(17)61205-9

* Boyd, S. P., & Vandenberghe, L. (2004). Convex Optimization. Cambridge, UK; New York: Cambridge University Press.

* Brandt, S. A., Bertin, J. J., Stiles, R. J., & Whitford, R. (2004). Introduction to Aeronautics: A Design Perspective, Second Edition. Reston, VA: American Institute of Aeronautics and Astronautics. https://doi.org/10.2514/4.862007

* Burston, M. T., Sabatini, R., Clothier, R., Gardi, A., & Ramasamy, S. (2014). Reverse engineering of a fixed wing unmanned aircraft 6-DoF model for navigation and guidance applications. *Applied Mechanics and Materials*, 629, 164–169. https://doi.org/10.4028/www.scientific.net/AMM.629.164

* Chen, H., & Allgöwer, F. (1998). A quasi-infinite horizon nonlinear model predictive control scheme with guaranteed stability. *Automatica*, 34(10), 1205–1217. https://doi.org/10.1016/S0005-1098(98)00073-9

* Cory, R., & Tedrake, R. (2008). Experiments in fixed-wing UAV perching. In *AIAA Guidance, navigation and control conference and exhibit*. Honolulu, Hawaii: AIAA. https://doi.org/10.2514/6.2008-7256

* Courant, R., & John, F. (1989). Introduction to Calculus and Analysis: Volume II. New York, NY: Springer New York. https://doi.org/10.1007/978-1-4613-8958-3

* Cunis, T., Condomines, J.-P., Burlion, L., & la Cour-Harbo, A. (2020). Dynamic stability analysis of aircraft flight in deep stall. *Journal of Aircraft*, 57(1), 143–155. https://doi.org/10.2514/1.C035455
FAA (2024). Aerobatic flight. https://www.ecfr.gov/current/title-14/part-91/section-91-303.

Faulwasser, T. (2013). Optimization-Based Solutions to Constrained Trajectory-Tracking and Path-Following Problems. Ph.D. thesis. Magdeburg, Universität, Diss., 2012.

Faulwasser, T., & Findeisen, R. (2011). A model predictive control approach to trajectory tracking problems via time-varying level sets of lyapunov functions. In *IEEE Conference on decision and control and european control conference* (pp. 3381–3386). Orlando, FL, USA: IEEE. https://doi.org/10.1109/CDCT.2011.6160492

Feroskhan, M., & Go, T. H. (2016). Dynamics of sideslip perching maneuver under dynamic stall influence. *Aerospace Science and Technology*, **50**, 220–233. https://doi.org/10.1016/j.ast.2016.01.005

Feroskhan, M., & Go, T. H. (2018). Performance of sideslip perching maneuver with an adaptive gain control feature. *Aerospace Science and Technology*, **78**, 648–660. https://doi.org/10.1016/j.ast.2018.05.026

Feroskhan, M., Zheng, Z., & Go, T. H. (2020). Solutions to planar aircraft perching problem utilizing sideslip maneuvering. *Journal of Aerospace Engineering*, **33**(6), 04020066. https://doi.org/10.1061/(ASCE)AS.1943-5525.0001189

Fontes, F. A. C. C. (2001). A general framework to design stabilizing nonlinear model predictive controllers. *Systems & Control Letters*, **42**(2), 127–143. https://doi.org/10.1016/S0167-6911(00)00084-0

Grötschel, M., Lovász, L., & Schrijver, A. (1988). Geometric Algorithms and Combinatorial Optimization (vol. 2). Algorithms and Combinatorics. Berlin, Heidelberg: Springer Berlin Heidelberg. https://doi.org/10.1007/978-3-642-97881-4

Grüne, L., & Pannek, J. (2017). Nonlinear Model Predictive Control. Communications and Control Engineering. Cham: Springer International Publishing. https://doi.org/10.1007/978-3-319-46024-6

Kim, D., & Park, S. (2022). Vision-assisted deep stall landing for a fixed-wing UAV. *Journal of Field Robotics*, . https://doi.org/10.1002/rob.22100

Köhler, J., Müller, M. A., & Allgöwer, F. (2020a). A nonlinear model predictive control framework using reference generic terminal ingredients. *IEEE Transactions on Automatic Control*, **65**(8), 3576–3583. https://doi.org/10.1109/TAC.2019.2949350

Köhler, J., Müller, M. A., & Allgöwer, F. (2020b). A nonlinear model predictive control framework using reference generic terminal ingredients – extended version. *IEEE Transactions On Automatic Control*, **65**(8), 3576–3583. https://doi.org/10.1109/TAC.2019.2949350

Lofberg, J. (2004). YALMIP: A toolbox for modeling and optimization in MATLAB. In *2004 IEEE International conference on robotics and automation* (pp. 284–289). Taipei, Taiwan: IEEE. https://doi.org/10.1109/CACSD.2004.1393890

Magni, L., & Sepulchre, R. (1997). Stability margins of nonlinear receding-horizon control via inverse optimality. *Systems & Control Letters*, **32**(4), 241–245. https://doi.org/10.1016/S0167-6911(97)00079-0

Mathisen, S., Gryte, K., Gros, S., & Johansen, T. A. (2021). Precision deep-stall landing of fixed-wing UAVs using nonlinear model predictive control. *Journal of Intelligent & Robotic Systems*, **101**(1), 24. https://doi.org/10.1007/s10846-020-01264-3

Mayne, D. Q. (2001). Control of constrained dynamic systems. *European Journal of Control*, **7**(2-3), 87–99. https://doi.org/10.3166/ejc.7.8.7-99

Mehanovic, D., Bass, J., Courteau, T., Rancourt, D., & Lussier Desbiens, A. (2017). Autonomous thrust-assisted perching of a fixed-wing UAV on vertical surfaces. In M. Mangan, M. Cutkosky, A. Mura, P. F. M. J. Verschure, T. Prescott, & N. Lepora (Eds.), *Biomimetic and biohybrid systems* (pp. 302–314). Cham: Springer International Publishing (vol. 10384). https://doi.org/10.1007/978-3-319-63537-8_26

Michalska, H., & Vinter, R. B. (1994). Nonlinear stabilization using discontinuous moving-horizon control. *IMA Journal of Mathematical Control and Information*, **11**(4), 321–340. https://doi.org/10.1093/imamci/11.4.321

Moore, J., Cory, R., & Tedrake, R. (2014). Robust post-stall perching with a simple fixed-wing glider using LQR-Trees. *Bioinspiration & Biomimetics*, **9**(2), 025013. https://doi.org/10.1088/1748-3182/9/2/025013

MOSEK (2024). The MOSEK optimization toolbox for MATLAB manual. Version 10.2. http://docs.mosek.com/latest/toolbox/index.html.

Müller, M., & Allgöwer, F. (2023). EECI Course on nonlinear and data-based model predictive control.

Nelson, R. C. (1998). Flight Stability and Automatic Control. (2nd ed.). Boston, Massachusetts: WCB/McGraw Hill.

Nguyen, H. T., Nguyen, N. T., & Prodan, I. (2024). Notes on the terminal region enlargement of a stabilizing NMPC design for a multicopter. *Automatica*, **159**, 111375. https://doi.org/10.1016/j.automatica.2023.111375

Nguyen, H. T., Prodan, I., & Fontes, F. A. C. C. (2023). Trajectory optimization and NMPC tracking for a fixed-wing UAV in deep stall with perch landing*. In *2023 European control conference (ECC)* (pp. 1–7). https://doi.org/10.23919/ECC57647.2023.10178188

Paiva, L. T., & Fontes, F. A. C. C. (2019). Sampled-data model predictive control: adaptive time-mesh refinement algorithms and guarantees of stability. *Discrete & Continuous Dynamical Systems - B*, **24**(5), 2335–2364. https://doi.org/10.3934/cdcsb.2019098

Phillips, W. F. (2004). Mechanics of Flight. Hoboken, N.J.: Wiley.

Pointner, W., Kotsis, G., Langthaler, P., Naderhirn, M., & Sense, A. (2011). Using formal methods to verify safe deep stall landing of a MAV. In *2011 IEEE/AIAA 30th digital avionics systems conference* (p. 10). https://doi.org/10.1109/DASC.2011.6096086

Puopolo, M., & Jacob, J. D. (2015). Model for longitudinal perch maneuvers of a fixed-wing unmanned aerial vehicle. *Journal of Aircraft*, **52**(6), 2021–2031. https://doi.org/10.2514/1.C033136

Roberts, J. W., Cory, R., & Tedrake, R. (2009). On the controllability of fixed-wing perching. In *2009 American control conference* (pp. 2018–2023). St. Louis, MO, USA: IEEE. https://doi.org/10.1109/ACC.2009.5160526

Sim, A. (1984). Flight characteristics of a manned, low-speed, controlled deep stallvehicle. In *11th atmospheric flight mechanics conference*. Seattle, WA, U.S.A.: American Institute of Aeronautics and Astronautics. https://doi.org/10.2514/6.1984-2074

Stengel, R. F. (2004). Flight Dynamics. Princeton, NJ: Princeton University Press.

Taniguchi, H. (2008). Analysis of deepstall landing for UAV. In *26th congress of international council of the aeronautical sciences. anchorage, alaska, USA*. https://www.icas.org/ICAS_ARCHIVE/ICAS2008/PAPERS/408.PDF.

Venkateswara Rao, D. M. K. K., & Go, T. H. (2014). Optimization, stability analysis, and trajectory tracking of perching maneuvers. *Journal of Guidance, Control, and Dynamics*, **37**(3), 879–888. https://doi.org/10.2514/1.60787

Venkateswara Rao, D. M. K. K., Tang, H., & Go, T. H. (2015). A parametric study of fixed-wing aircraft perching maneuvers. *Aerospace Science and Technology*, **42**, 459–469. https://doi.org/10.1016/j.ast.2015.02.006

Venkateswara Rao, D. M. K. K., & Yanhua, W. (2019). A vision-based navigation system for perching aircraft. *Journal of Intelligent & Robotic Systems*, **95**(2), 555–566. https://doi.org/10.1007/s10846-018-0807-7

Verschueren, R., Frison, G., Kouzoupis, D., Frey, J., Duijkeren, N. V., Zanelli, A., Novoselnik, B., Albin, T., Quirynen, R., & Diehl, M. (2022). Acados—a modular open-source framework for fast embedded optimal control. *Mathematical Programming Computation*, **14**(1), 147–183. https://doi.org/10.1007/s12532-021-00208-8

Wächter, A., & Biegler, L. T. (2006). On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming. *Mathematical Programming*, **106**(1), 25–57. https://doi.org/10.1007/s10107-004-0559-y

Waldock, A., Greatwood, C., Salama, F., & Richardson, T. (2018). Learning to perform a perched landing on the ground using deep reinforcement learning. *Journal of Intelligent & Robotic Systems*, **92**(3-4), 685–704. https://doi.org/10.1007/s10846-017-0696-1

White, M. D., & Cooper, G. E. (1965). Simulator studies of the deep stall. In *NASA Conference on aircraft operating problems: a compilation of the papers presented* (pp. 101–111). https://ntrs.nasa.gov/citations/20000011980.

Wickenheiser, A. M., & Garcia, E. (2008). Optimization of perching maneuvers through vehicle morphing. *Journal of Guidance, Control, and Dynamics*, **31**(4), 815–823. https://doi.org/10.2514/1.33819

Wüest, V., Jeger, S., Feroskhan, M., Ajanic, E., Bergonti, F., & Floreano, D. (2024). Agile perching maneuvers in birds and morphing-wing drones. *Nature Communications*, **15**(1), 8330. https://doi.org/10.1038/s41467-024-52369-4