# Dynamical Model and Tracking Controller for a Novel Flapping Wing Drone Platform

Jossué Cariño Escobar, Lina Le-Guellec, Edwin Van Ruymbeke, Nicolas Marchand, Thomas Engels, Franck Ruffier

## ► To cite this version:

Jossué Cariño Escobar, Lina Le-Guellec, Edwin Van Ruymbeke, Nicolas Marchand, Thomas Engels, et al.. Dynamical Model and Tracking Controller for a Novel Flapping Wing Drone Platform. 2025 International Conference on Unmanned Aircraft Systems (ICUAS), May 2025, Charlotte, NC, United States. 10.1109/ICUAS65942.2025.11007922. [hal-05001130v2](https://hal.science/hal-05001130v2)

HAL Id: hal-05001130

https://hal.science/hal-05001130v2

Submitted on 17 Apr 2025

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 Attribution 4.0 International License
# Dynamical Model and Tracking Controller for a Novel Flapping Wing
Drone Platform*

Jossué Cariño Escobar¹, Lina Le-Guellec², Edwin Van Ruymbeke³,
Nicolas Marchand², Thomas Engels¹ and Franck Ruffier¹

**Abstract**— This work focuses on the design and control of a novel type of Flapping-Wing Micro Aerial Vehicle (FWMAV). The drone, known as the X-Fly, is a new under-actuated robotic platform that also has an inner control loop to stabilize its roll angle thanks to an onboard IMU. Such assistance makes the X-Fly easier to pilot. The under-actuation and the flapping oscillations make the modelling and the control of the X-Fly a challenging task.

A dynamical model is introduced that is able to take advan-
tage of the stabilized roll dynamics to separate the platform
into two almost independent sub-systems, one for the altitude
and another for the position on the x-y plane.

A trajectory tracking controller for the altitude and a
circular trajectory are then proposed and tested in order to
corroborate the validity of the presented model.

I. INTRODUCTION

Flapping-Wing Micro Aerial Vehicles (FWMAVs), also
called ornithopters, have started to gain significant attention
in the last decade. In contrast with other types of fixed-wing
Unmanned Aerial Vehicles (UAVs), FWMAVs have a higher
thrust-to-weight ratio, agility, and maneuverability [1]. They
are safer and lighter than similarly sized multi-rotor UAVs
while retaining hover-capability, high maneuverability, and
efficiency at small scales. FWMAVs have gained significant
attention in research, particularly in the fields of bio-robotic
and bio-mimetics. However, their complex dynamics pose a
substantial challenge and has become the focus of control
and mathematical modeling research.

Some of the early studies focus on the design of the phys-
ical platform, with a special focus on agility and hovering
capabilities, like [2]–[5]. Research into the altitude control
can be seen in works like [6], [7].

More recent advances have been made in the field of
robotics using FWMAVs, such as:

*J.C.E. support was provided by the ANR Bucolyc grant (ANR-23-CE51-0037) to F.R. Also, L.L-G. was supported by a doctoral fellowship grant from Univ. Grenoble Alpes and the French Ministry of Defense (AID - Agence Innovation Défense). J.C.E., L.L-G., and F.R. were also supported by Aix Marseille University and the CNRS (Life Science, Information Science, and Engineering and Science & technology Institutes). The facilities for the experimental tests has been mainly provided by ROBOTEX 2.0 (Grants ROBOTEX ANR-10-EQPX-44-01 and TIRREX ANR-21-ESRE-0015).

¹Jossué Cariño Escobar, Thomas Engels & Franck Ruffier are with Aix-Marseille Université, CNRS, ISM, Marseille, France {jossue.carino-escobar, thomas.engels} @univ-amu.fr, franck.ruffier@cnrs.fr

²Lina Le-Guellec and Nicolas Marchand are with Univ. Grenoble Alpes, CNRS, GIPSA-lab, Grenoble, France lina.le-guellec@grenoble-inp.fr, nicolas.marchand@grenoble-inp.fr

³Edwin Van Ruymbeke is with XTIM BionicBird, Marseille, France

• Obstacle avoidance in [8], [9].

• Trajectory tracking in [10]–[13].

This paper introduces a comprehensive investigation into
the design, modeling, and experimental validation of a novel
flapping wing drone system based on the patent [14]. It
makes use of two actuators. The wing motor is in charge of
the flapping frequency using a gear mechanism. The second
motor is attached to an aerodynamic control surface the acts
like a rudder, but it also tenses or relaxes the wings, altering
both the yaw and the roll dynamics with a single actuator.

<figure><img src="image_2.png" /><figcaption>Fig. 1. An image on the X-Fly commercial platform that was modelled and controlled.</figcaption></figure>

The innovation of this new commercially available vehicle
named **X-fly** is the inclusion of an inner control loop based
on an IMU. The embedded IMU helps to stabilize the roll
dynamics, in contrast with the previous flapping-wing drone
–Metafly– that was devoid of IMU. For this reason, the roll
dynamics was stabilized by the help of the Vicon motion
capture system in [12].

The complex dynamics of FWMAVs have been studied extensively in the literature in works like [15], [16]. There has even been an interest in black-box, as seen in [10], [12], [17], and grey-box system identification tools, as seen in [18], [19], to obtain the dynamical model of ornithopters. The main issue with this type of techniques is that the model is often an approximation around an operation region, and they might fail if they diverge too much from this operating region.

We propose a mathematical model representation that
extends the traditional fixed-wing approximation by con-
sidering a unique force direction in the body frame, as
seen in [20], aligning with the drone's specific flapping
mechanism. This model integrates the tail as a fixed control
surface, linking translational and altitude dynamics into an
underactuated system, unlike the work in [21] that do not consider the heading actuator dynamics into consideration to simplify the model. It is proven that one of the desired effects of the tail on the system is the stability of the pitch angle, which in turn allows the flapping to generate lift.

The introduction to the notation used throughout this work and the model itself is presented in Section II.

To manage the underactuated nature of our FWMAV model, we design a new controller based on Lyapunov formalism. This controller exhibits robustness by operating effectively from takeoff, without requiring the drone's state to be near the desired state initially (Section III).

While many previous works have presented controllers for FWMAVs theoretically and validated them through simulations, our approach emphasizes experimental corroboration. We conducted extensive experiments to validate both our mathematical model representation and the performance of our new controller in practical scenarios (Section IV).

This paper contributes to the field of flapping wing drones by proposing a novel mathematical model that captures the unique dynamics of our FWMAV system, presenting a robust controller that operates effectively from takeoff, and experimentally validating these contributions. The remainder of this paper is structured as follows: Section II details our mathematical model representation, Section III introduces our new controller design, Section IV presents our experimental corroboration, and Section V concludes the paper with future research directions.

## A. NOTATION

Throughout this work, scalar variables are indicated by a lowercase letter $a \in \mathbb{R}$, while $n$ length vectors are denoted using overlined lowercase letters $\bar{a} \in \mathbb{R}^n$. A 3D vector is displayed with an arrow $\vec{a} \in \mathbb{R}^3$, and the reference frame $\mathcal{A}$ where it was defined can also be joined like so $^A\bar{a}$. Unitary vectors are shown with a hat $\hat{a}$ and are usually obtained by normalizing its homologous vector $\vec{a}$.

Matrices are presented with uppercase letters $A \in \mathbb{R}^{n \times n}$. Rotation matrices from the reference frame $\mathcal{A}$ to the reference frame $\mathcal{B}$ are indicated as $R^{A \to B}$.

## II. SYSTEM MODEL

Let us define the system state vector as

$$ \bar{x} := \begin{bmatrix} x & y & z & \dot{x} & \dot{y} & \dot{z} & \phi & \theta & \psi & \omega_x & \omega_y & \omega_z \end{bmatrix}^T \\ := \begin{bmatrix} \dot{p}^T & \dot{\bar{p}}^T & \bar{\eta}^T & \bar{\omega}^T \end{bmatrix}^T \quad (1) $$

where $\bar{p} := \begin{bmatrix} x & y & z \end{bmatrix}^T$, $\mathcal{B}\bar{\omega} := \bar{\omega} := \begin{bmatrix} \omega_x & \omega_y & \omega_z \end{bmatrix}^T$ are the position defined in the inertial frame $\mathcal{I}$ and angular speed defined in the body frame $\mathcal{B}$, respectively. The attitude is denoted using Euler angles $\bar{\eta} := \begin{bmatrix} \phi & \theta & \psi \end{bmatrix}^T$ and $\dot{\bar{p}}$ is the translational speed in the inertial frame.

The control vector is defined as

$$ \bar{u} := \begin{bmatrix} u_m & u_r \end{bmatrix}^T \quad (2) $$

where $u_m \in \mathbb{R}_+$ is the mean forward force due to the flapping movement and is assumed to be proportional to the flapping frequency. The term $u_r \in \mathbb{R}$ is the angle of the rudder aerodynamic surface.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Free body diagram of the flapping bird drone with the position and orientation of the body frame $\mathcal{B}$ w.r.t. the inertial frame $\mathcal{I}$. The rudder is located on the $xz$ plane while the tail is perpendicular to it.</figcaption></figure>

The system model of the vehicle as shown in Figure 2 is obtained using the Newton-Euler formalism.

$$ m \ddot{\bar{p}} = R^{\mathcal{B} \to \mathcal{I}} \begin{bmatrix} \begin{pmatrix} u_m \\ 0 \\ 0 \end{pmatrix} + \bar{f}_t + \bar{f}_r \end{bmatrix} + m \bar{g} \quad (3) $$

$$ J \dot{\bar{\omega}} = \bar{\tau}_{\phi} + \bar{\tau}_t + \bar{\tau}_r - \bar{\omega} \times J \omega - D_{\bar{\omega}} \bar{\omega} \quad (4) $$

The vectors $\bar{f}_t, \bar{f}_r$ represents the aerodynamical forces derived from the tail and the rudder, respectively. The vector $\bar{g} := (0 \ 0 \ -g)$ is the gravity's effect $g \in \mathbb{R}$ defined on the inertial frame $\mathcal{I}$. The constants $m \in \mathbb{R}_+$, $J \in \mathbb{R}_+^{3 \times 3}$ are the vehicle's mass and inertia matrix, respectively. The inertia matrix is assumed to be diagonal $J := \text{diag} \ [J_x \ J_y \ J_z]$ with positive components $J \in \mathbb{R}_+$, even though the vehicle only has one symmetry axis. This assumption implies that the off-diagonal elements are relatively small compared to the diagonal ones. The drag effects in the attitude system are considered in the matrix $D_{\bar{\omega}} := \text{diag} \ [d_{\bar{\omega},x} \ d_{\bar{\omega},y} \ d_{\bar{\omega},z}] \in \mathbb{R}^{3 \times 3}$, which will also be considered a diagonal matrix. The torque $\bar{\tau}_{\phi}$ is the effect of an internal roll inner controller as shown in Figure 7.

The vector $\bar{\tau}_t$ is the vector that denotes the distance from the center of mass of the vehicle to the tail's aerodynamic surface center as shown in Figure 3. The length $\bar{r}_t$ is the distance from the center of gravity to the rudder's aerodynamical surface center, as shown in Figure 4, and is defined as:

$$ \bar{r}_t := \bar{r}_{r,1} - l_{r,2} \begin{pmatrix} \cos u_r \\ \sin u_r \\ 0 \end{pmatrix} \quad (5) $$

The tail force $\bar{f}_t$ is calculated as:

$$ \bar{f}_t = \frac{1}{2} \rho d_t \hat{n}_t (\hat{n}_t^T \mathcal{B} \dot{\bar{p}})^2 \operatorname{sgn} (\hat{n}_t^T \mathcal{B} \dot{\bar{p}}) \quad (6) \\ := c_t \hat{n}_t \hat{n}_t^T \mathcal{B} \dot{\bar{p}} $$
where the vector $\hat{n}_t := (\sin \theta_t \ 0 \ \cos \theta_t)^T$ is the normal to the tail's surface and depends on the tail's angle $\theta_t$ as seen on Figure 3. The constants $\rho, d_t \in \mathbb{R}_+$ depict the air's density and the tail's drag coefficient, respectively. The vector $\mathcal{B}\vec{p} := R^{\mathcal{I}\to\mathcal{B}}\vec{p}$ is the vehicle's speed in the body frame. Some of the parameters can be grouped into $c_t \in \mathbb{R}$. Let us note that $\vec{f}_t$ only has components in the body axes $\mathcal{B}x$ and $\mathcal{B}y$, defined as $f_{t,x}$ and $f_{t,z}$ respectively.

The rudder force $\vec{f}_r$ is obtained as:

$$ \vec{f}_r = \frac{1}{2} \rho d_r \hat{n}_r (\hat{n}_r^T \mathcal{B} \dot{\vec{p}})^2 \operatorname{sgn} (\hat{n}_r^T \mathcal{B} \dot{\vec{p}}) \quad (7) \\ := c_r \hat{n}_r \hat{n}_r^T \mathcal{B} \dot{\vec{p}} $$

where $d_r$ is the rudder's drag coefficient and $\hat{n}_r := (-\sin u_r \ \cos u_r \ 0)^T$. Some of the parameters can be grouped into $c_r \in \mathbb{R}$. Let us note that $\vec{f}_r$ only has components in the body axes $\mathcal{B}x$ and $\mathcal{B}y$, defined as $f_{r,x}$ and $f_{r,y}$ respectively.

The torques due to the tail's aerodynamic surface $\vec{\tau}_t$ and the rudder $\vec{\tau}_r$ depend on their corresponding aerodynamical forces as:

$$ \vec{\tau}_t := \vec{r}_t \times \vec{f}_t \quad (8) $$

$$ \vec{\tau}_r := \vec{r}_r \times \vec{f}_r \quad (9) $$

such that $\vec{\tau}_t$ only has components in the $\mathcal{B}y$ plane, and $\vec{\tau}_r$ only has components in the $\mathcal{B}z$ plane. These components are defined as $\tau_{t,y}$ and $\tau_{r,z}$ respectively.

### A. Attitude Description

The model presented in (3) and (4) can be simplified due to the presence of the roll inner controller $\vec{\tau}_\phi$ as seen in Figure 7. This controller guarantees that the dynamics of the roll state $\bar{x}_\phi := [\phi \ \omega_x] \to \bar{0}$ are stable. The stability of $\phi$ towards zero, and the fact it is defined in the body frame, implies that the rotation matrix $R^{\mathcal{I}\to\mathcal{B}}$ can be assumed to only consider the effects of the yaw $\psi$ and pitch $\theta$ angles as:

$$ \begin{align} R^{\mathcal{I}\to\mathcal{B}} &= R(\theta) R(\psi) := \begin{bmatrix} \cos\theta & 0 & \sin\theta \\ 0 & 1 & 0 \\ -\sin\theta & 0 & \cos\theta \end{bmatrix} \nonumber \\ &= \begin{bmatrix} \cos\psi & \sin\psi & 0 \\ -\sin\psi & \cos\psi & 0 \\ 0 & 0 & 1 \end{bmatrix} \nonumber \\ &= \begin{bmatrix} \cos\theta & \cos\psi & \cos\theta & \sin\psi & \sin\theta \\ -\sin\psi & \cos\psi & 0 \\ -\sin\theta & \cos\psi & -\sin\theta & \sin\psi & \cos\theta \end{bmatrix} \nonumber \\ \implies R^{\mathcal{B}\to\mathcal{I}} &= \begin{bmatrix} \cos\theta & \cos\psi & -\sin\psi & -\sin\theta & \cos\psi \\ \cos\theta & \sin\psi & \cos\psi & -\sin\theta & \sin\psi \\ \sin\theta & 0 & \cos\theta & 0 & \sin\theta \end{bmatrix} \tag{10} \end{align} $$

In turn, this allows us to separate the full model into two almost independent models based on the yaw and pitch dynamics, one for controlling the altitude of the system $z$, and another to control its position on the $xy$ plane.

Another consequence of this is that the attitude kinematics can be described as:

$$ \begin{pmatrix} \dot{\theta} \\ \dot{\psi} \end{pmatrix} = \begin{bmatrix} 0 & 1 & 0 \\ \sin\theta & 0 & \cos\theta \end{bmatrix} \vec{\omega} \quad (11) $$

Furthermore, it will be shown in Section II-B that the pitch dynamics are stabilizable.

### B. Altitude Dynamics

<figure><img src="image_4.png" /><figcaption>Fig. 3. Free body diagram of the flapping bird drone lateral movement in the x-z body axis.</figcaption></figure>

The altitude dynamics are depicted in Figure 3. They are defined in the intermediate frame defined by the yaw $\psi$ using $R(\psi)$. Its speed state can be described with respect to $\dot{\vec{p}}$ with a change in variables:

$$ \begin{align} \dot{\vec{p}}_{xz} &:= \begin{bmatrix} \dot{x}_{xz} \\ \dot{z} \end{bmatrix} = \begin{bmatrix} \cos \psi & \sin \psi & 0 \\ 0 & 0 & 1 \end{bmatrix} \dot{\vec{p}} \tag{12} \\ &:= E_{xz} \dot{\vec{p}} \nonumber \end{align} $$

where $\dot{x}_{xz} \in \mathbb{R}$ is the horizontal speed in the intermediate reference frame and its dynamics are derived from (3) and (11) as:

$$ \begin{align} m \ddot{\vec{p}}_{xz} &= E_{xz} \ddot{\vec{p}} + \begin{bmatrix} (\cos \psi \dot{y} - \sin \psi \dot{x}) \\ 0 \end{bmatrix} \dot{\psi} \nonumber \\ &= (u_m + f_{r,x} + f_{t,x}) \begin{bmatrix} \cos \theta \\ \sin \theta \end{bmatrix} + f_{t,z} \begin{bmatrix} -\sin \theta \\ \cos \theta \end{bmatrix} \nonumber \\ &\quad + \begin{bmatrix} (\cos \psi \dot{y} - \sin \psi \dot{x})(\omega_x \sin \theta + \omega_z \cos \theta) \\ -mg \end{bmatrix} \tag{13} \end{align} $$

which comes from the fact that

$$ E_{xz} R^{\mathcal{B}\to\mathcal{I}} = \begin{bmatrix} \cos\theta & 0 & -\sin\theta \\ \sin\theta & 0 & \cos\theta \end{bmatrix} \quad (14) $$

The model (13) can be further simplified for control purposes if one considers $\omega_x \to 0$ and if the effects of the aerodynamic forces $\vec{f}_t, \vec{f}_r$ are grouped into a single drag force as:
$$
\begin{align}
m \ddot{\vec{p}}_{xz} &\approx u_m \begin{bmatrix} \cos \theta \\ \sin \theta \end{bmatrix} - d_{xz} \dot{\vec{p}}_{xz} \nonumber \\
&\quad + \begin{bmatrix} (\cos \psi \dot{y} - \sin \psi \dot{x}) \omega_z \cos \theta \\ -mg \end{bmatrix} \tag{15} \\
&:= u_m \begin{bmatrix} \cos \theta \\ \sin \theta \end{bmatrix} - d_{xz} \dot{\vec{p}}_{xz} + \begin{bmatrix} \xi_{xz} \\ -mg \end{bmatrix} \nonumber
\end{align}
$$

where $d_{xz} \in \mathbb{R}_+$ is considered a positive definite constant drag coefficient. The term $\xi_{xz} \in \mathbb{R}$ are the effects of the translational dynamics on the altitude system and can be considered as a perturbation.

Assuming that the angular velocity $\omega_x \to 0$, the pitch dynamics can be obtained from equations (4), (8) and (11), and using the pitch difference $\Delta\theta := \theta - \theta_t$ as:

$$
\begin{align}
\dot{\theta} &= \omega_y \tag{16} \\
J_y \dot{\omega}_y &= \tau_{t,y} - d_{\vec{\omega},y} \omega_y := c_{\tau_t} \hat{n}_t^T \mathcal{B} \dot{\vec{p}} - d_{\vec{\omega},y} \omega_y \nonumber \\
&= c_{\tau_t} \hat{n}_t^T R^{\mathcal{I} \to \mathcal{B}} \dot{\vec{p}} - d_{\vec{\omega},y} \omega_y \nonumber \\
&= c_{\tau_t} \begin{bmatrix} (\sin \theta_t \cos \theta - \cos \theta_t \sin \theta) \cos \psi \\ (\sin \theta_t \cos \theta - \cos \theta_t \sin \theta) \sin \psi \\ \sin \theta_t \sin \theta + \cos \theta_t \cos \theta \end{bmatrix}^T \dot{\vec{p}} - d_{\vec{\omega},y} \omega_y \nonumber \\
&= c_{\tau_t} \begin{bmatrix} -\sin \Delta\theta & | & \cos \Delta\theta \end{bmatrix} E_{xz} \dot{\vec{p}} - d_{\vec{\omega},y} \omega_y \nonumber \\
&= c_{\tau_t} \begin{bmatrix} -\sin \Delta\theta & | & \cos \Delta\theta \end{bmatrix} \dot{\vec{p}}_{xz} - d_{\vec{\omega},y} \omega_y \tag{17}
\end{align}
$$

where the variable $c_{\tau_t} \in \mathbb{R}$ is a scalar.

Thus the altitude model is described by equations (15), (16) and (17).

### C. Translational Dynamics on the x-y Plane

The translational dynamics are depicted in Figure 4. They are defined in the $x-y$ plane on the inertial frame $\mathcal{I}$. Its dynamics are derived for the position $\vec{p}_{xy} := [x \ y]^T$ and (3) as:

$$
\begin{align}
m \ddot{\vec{p}}_{xy} &= (u_m + f_{r,x} + f_{t,x}) \cos \theta \begin{bmatrix} \cos \psi \\ \sin \psi \end{bmatrix} \nonumber \\
&\quad + f_{r,y} \begin{bmatrix} -\sin \psi \\ \cos \psi \end{bmatrix} - (f_{t,z}) \sin \theta \begin{bmatrix} \cos \psi \\ \sin \psi \end{bmatrix} \tag{18}
\end{align}
$$

A similar analysis can be made as in (15) to simplify the effect of the forces due to the aerodynamic surfaces:

$$
m \ddot{\vec{p}}_{xy} \approx u_m \cos \theta \begin{bmatrix} \cos \psi \\ \sin \psi \end{bmatrix} - d_{xy} \dot{\vec{p}}_{xy} \quad (19)
$$

where $d_{xy} \in \mathbb{R}_+$ is a positive definite constant drag coefficient.

Assuming that the angular velocity $\omega_x \to 0$, the yaw dynamics can be obtained from equations (4), (9) and (11) as:

$$
\dot{\psi} = \omega_z \cos \theta \quad (20)
$$

$$
J_z \dot{\omega}_z = \tau_{r,z} - d_{\vec{\omega},z} \omega_z := c_{\tau_r} \sin u_r - d_{\vec{\omega},z} \omega_z \quad (21)
$$

<figure><img src="image_5.png" /><figcaption>Fig. 4. Free body diagram of the flapping bird drone horizontal movement in the x-y inertial axis.</figcaption></figure>

where the variable $c_{\tau_r} \in \mathbb{R}$ is a scalar.

Thus the translational model is described by equations (19), (20) and (21).

## III. CONTROLLER DESIGN

### A. Altitude Controller

Let us define a sub-state based on the altitude error $z_e := z - z_d$, where $z_d \in \mathbb{R}$ is the desired altitude, as:

$$
\bar{x}_z := \left[ \int z_e dt \quad z_e \quad \dot{z} \right]^T \quad (22)
$$

$$
\begin{align}
\dot{\bar{x}}_z &= \begin{bmatrix} 0 & 1 & 0 \\ 0 & 0 & 1 \\ 0 & 0 & -d_{xz} \end{bmatrix} \bar{x}_z + \begin{bmatrix} 0 \\ 0 \\ \sin \theta \end{bmatrix} u_m - \begin{bmatrix} 0 \\ 0 \\ mg \end{bmatrix} \nonumber \\
&:= A_z \bar{x}_z + B_z u_m + G_z \tag{23}
\end{align}
$$

Another sub-state can be defined for the pitch as:

$$
\bar{x}_{\theta} := \begin{bmatrix} \Delta\theta & \omega_y \end{bmatrix} \quad (24)
$$

$$
\begin{align}
\dot{\bar{x}}_{\theta} &= \begin{bmatrix} 0 & 1 \\ 0 & -d_{\vec{\omega},y} \end{bmatrix} \bar{x}_{\theta} + \begin{bmatrix} 0 \\ c_{\tau_t}(\dot{z} \cos \Delta\theta - \dot{x}_{xz} \sin \Delta\theta) \end{bmatrix} \nonumber \\
&\approx \begin{bmatrix} 0 & 1 \\ -c_{\tau_t} \dot{x}_{xz} & -d_{\vec{\omega},y} \end{bmatrix} \bar{x}_{\theta} + \begin{bmatrix} 0 \\ c_{\tau_t} \dot{z} \cos \Delta\theta \end{bmatrix} \nonumber \\
&:= A_{\theta} \bar{x}_{\theta} + G_{\theta} \tag{25}
\end{align}
$$

**Theorem 1:** Let us define the altitude controller $u_m$ as:

$$ u_m := -K_z \bar{x}_z \quad (26) $$
where $K_z \in \mathbb{R}^{1 \times 3}$ is a state feedback gain matrix designed such that it solves the Lyapunov equation

$$ (A_z - B_z K_z + G_z)^T P_z + P_z (A_z - B_z K_z + G_z) = -Q_z \quad (27) $$

on the trajectories of system (23) for certain positive symmetric matrices $P_z, Q_z \in \mathbb{R}_+^{3 \times 3}$.

Let us further assume that the pitch angle is bounded as $-\pi/2 < \theta < \pi/2$, and the effects of the aerodynamical forces $c_{\tau_r}$ and the lateral speed $\dot{x}_{xz}$ are bounded as

$$ 0 \le \dot{x}_{xz} < l_{\dot{x}_{xz}} \quad (28) $$

$$ 0 \le c_{\tau_r} < l_{c_{\tau_r}} \quad (29) $$

where the positive scalars $l_{\dot{x}_{xz}}, l_{c_{\tau_r}} \in \mathbb{R}_+$ represent the bounds.

Then the system described by the dynamical equations (15), (16) and (17) is uniformly stable in the sense of Lyapunov.

*Proof:* The bounds (28) and (29) applied to the model (25) imply that a pair of matrices $P_\theta, Q_\theta \in \mathbb{R}^{2 \times 2}$ exist such that:

$$ A_{\theta}^{T} P_{\theta} + P_{\theta} A_{\theta} = -Q_{\theta} \quad (30) $$

$$ \left[ \begin{array}{cc} \|Q_z\| & \|P_\theta\| \\ \|P_\theta\| & \|Q_\theta\| \end{array} \right] > 0 \quad (31) $$

Let us propose the following Lyapunov candidate function:

$$ V_z := \bar{x}_z^T P_z \bar{x}_z + \bar{x}_\theta^T P_\theta \bar{x}_\theta \quad (32) $$

whose differentiation w.r.t. time is:

$$ \dot{V}_z = 2\bar{x}_z^T P_z \dot{\bar{x}}_z + 2\bar{x}_\theta^T P_\theta \dot{\bar{x}}_\theta \quad (33) $$

Substituting from (27) and (30) into (33) results in:

$$ \begin{aligned} \dot{V}_z &= -\bar{x}_z^T Q_z \dot{\bar{x}}_z - \bar{x}_\theta^T Q_\theta \dot{\bar{x}}_\theta + 2\bar{x}_\theta^T P_\theta G_\theta \\ &\le -\|\bar{x}_z\|^2 \|Q_z\| + 2 \|\bar{x}_z\| \|P_\theta\| \|\bar{x}_\theta\| - \|\bar{x}_\theta\|^2 \|Q_\theta\| \end{aligned} \quad (34) $$

Due to the condition (31), this implies that $\dot{V}_z \le 0$, thus concluding the proof.

For a fixed tail $\theta_t$, tail torque $\tau_t$ changes sign when forward speed increases

<figure><img src="image_6.png" /><figcaption>Fig. 5. A graphical representation of the stability of $\theta$ due to the tail and the lateral speed</figcaption></figure>

Conditions (28) and (29) imply that the ornithopter is moving forwards. A more intuitive representation of the stability of the pitch system can be seen in Figure 5.

## B. Position Controller

Let us define a radius in the $x-y$ inertial frame as $r := \|\vec{p}_{xy}\|$. The control objective requires that the vehicle remains at a desired constant radius $r_d \in \mathbb{R}_+$. A robust way to achieve this is if the system has similar behavior to a limit cycle, which is explained in Theorem 2.

*Theorem 2:* Let us define the following desired lateral speed $\dot{\vec{p}}_{xy,d}$ as:

$$ \dot{\vec{p}}_{xy,d} := \begin{pmatrix} -y \\ x \end{pmatrix} - \mu \vec{p}_{xy} r_e \quad (35) $$

where $r_e := r^2 - r_d^2$ is the squared radius error and $\mu \in \mathbb{R}_+$ is a positive constant gain that determines how fast the vehicle reaches the desired radius.

If the vehicle's lateral speed tracks (35), then there is an unstable spiral at the origin surrounded by a stable limit cycle which grows out of the origin at exactly $r_d$ as seen in Figure 6.

<figure><img src="image_7.png" /><figcaption>Fig. 6. Development of a limit cycle using the proposed modified Hopf bifurcation. The $\mu$ parameter affects the convergence rate.</figcaption></figure>

*Proof:* The dynamics of (35) are derived from a Hopf bifurcation, as seen in [22, p. 12.4]. The main difference is that the limit cycle is always located at $r_d$, and the parameter $\mu$ is instead used as a design gain. If the speed of the vehicle has the dynamics proposed in (35), then this implies that $r \to r_d$ as $t \to \infty$.

Equation (35) gives us a desired speed for the position controller to track. Its differentiation w.r.t. time is also required for the tracking scheme.

$$ \begin{aligned} \ddot{\vec{p}}_{xy,d} &= \begin{pmatrix} -\dot{y} \\ \dot{x} \end{pmatrix} - \mu \dot{\vec{p}}_{xy,d} r_e - \mu \vec{p}_{xy} \dot{r}_e \\ &= \begin{pmatrix} -\dot{y} \\ \dot{x} \end{pmatrix} - \mu [r_e - 2 \vec{p}_{xy} \vec{p}_{xy}^T] \dot{\vec{p}}_{xy,d} \end{aligned} \quad (36) $$

*Theorem 3:* Let us define the control law $u_r$ for the translational system as:

$$ u_r := -K_{xy} \bar{x}_{xy} \quad (37) $$

where $K_{xy} \in \mathbb{R}^{1 \times 3}$ is a positive gain matrix, and the vector $\bar{x}_{xy} \in \mathbb{R}^{3 \times 1}$ is defined as:
<figure><img src="image_8.png" /></figure>

$$ \bar{x}_{xy} := \begin{bmatrix} \int \dot{\bar{p}}_{xy,d}^T \begin{pmatrix} -\sin \psi \\ \cos \psi \\ \omega_z \end{pmatrix} dt \\ \dot{\bar{p}}_{xy,d}^T \begin{pmatrix} -\sin \psi \\ \cos \psi \\ \omega_z \end{pmatrix} \end{bmatrix} = \begin{bmatrix} \int \|\dot{\bar{p}}_{xy,d}\| \sin \Delta \psi dt \\ \|\dot{\bar{p}}_{xy,d}\| \omega_z \sin \Delta \psi \end{bmatrix} \quad (38) $$

where $\Delta \psi \in \mathbb{R}$ is the angle difference between the desired speed vector $\dot{\bar{p}}_{xy,d}$ and the heading $\psi$.

Furthermore, let us assume that $K_{xy}$ and a positive symmetric matrix $P_{xy} \in \mathbb{R}^{3 \times 3}$ are chosen to satisfy the following Lyapunov equation:

$$ (A_{xy} - B_{xy} K_{xy})^T P_{xy} + P_{xy} (A_{xy} - B_{xy} K_{xy}) = -Q_{xy} \quad (39) $$

where the matrices $A_{xy} \in \mathbb{R}^{3 \times 3}$ and $B_{xy} \in \mathbb{R}^{3 \times 1}$ are defined as:

$$ A_{xy} := \begin{bmatrix} 0 & 1 & 0 \\ 0 & 0 & \|\dot{\bar{p}}_{xy,d}\| \cos \Delta \psi \\ 0 & 0 & -d_{\bar{\omega},z} \end{bmatrix} \quad (40) $$

$$ B_{xy} := \begin{bmatrix} 0 \\ 0 \\ c_{\tau_r} \end{bmatrix} \quad (41) $$

Let us further assume that:

$$ \Delta \dot{\psi} \approx \omega_z \quad (42) $$

$$ \dot{\bar{p}}_{xy,d}^T \ddot{\bar{p}}_{xy,d} \approx 0 \quad (43) $$

$$ \sin u_r \approx u_r \quad (44) $$

The control law (37) stabilizes the system described in (19), (20) and (21) in the sense of Lyapunov, if one assumes that the altitude system is stable i.e. $\frac{d}{dt}(u_m \cos \theta) \to 0$.

*Proof:* Let us propose the following Lyapunov function:

$$ V_{xy} := \bar{x}_{xy}^T P_{xy} \bar{x}_{xy} \quad (45) $$

whose differentiation w.r.t. time is

$$ \dot{V}_{xy} = (A_{xy} - B_{xy} K_{xy})^T P_{xy} + P_{xy} (A_{xy} - B_{xy} K_{xy}) \quad (46) $$

where the matrices $A_{xy} \in \mathbb{R}^{3 \times 3}$ and $B_{xy} \in \mathbb{R}^{3 \times 1}$ are obtained from the differentiation of $\bar{x}_{xy}$, defined in (38), and the assumptions (42), (43), and (44) as:

$$ \begin{aligned} \dot{\bar{x}}_{xy} &= \begin{bmatrix} \|\dot{\bar{p}}_{xy,d}\| \sin \Delta \psi \\ \Delta \dot{\psi} \|\dot{\bar{p}}_{xy,d}\| \cos \Delta \psi + \frac{\dot{\bar{p}}_{xy,d}^T \ddot{\bar{p}}_{xy}}{\|\dot{\bar{p}}_{xy,d}\|} \sin \Delta \psi \\ c_{\tau_r} \sin u_r - d_{\bar{\omega},z} \omega_z \end{bmatrix} \\ \approx & \begin{bmatrix} 0 & 1 & 0 \\ 0 & 0 & \|\dot{\bar{p}}_{xy,d}\| \cos \Delta \psi \\ 0 & 0 & -d_{\bar{\omega},z} \end{bmatrix} \bar{x}_{xy} + \begin{bmatrix} 0 \\ 0 \\ c_{\tau_r} \end{bmatrix} u_r \end{aligned} \quad (47) $$

∴ $\dot{V}_{xy} = -Q_{xy}$ from (39), thus concluding the proof.

Assumptions (42) and (43) can be regulated by the parameter $\mu$ from (35).

<figure><img src="image_9.png" /><figcaption>Fig. 7. Block diagram of the control system.</figcaption></figure>

The block diagram of the control system, including the inner control loop and outer control loop, are shown in Figure 7.

## IV. RESULTS

### A. Simulation Results

The simulations were made using the dynamic equations (15), (16), (17), (19), (20) and (21). The parameters used for the simulation are shown in Table I.

<table><caption>TABLE I
PARAMETERS FOR THE SIMULATION</caption><tr><td>m =</td><td>1.12e-2 g</td><td>J =</td><td>diag [ 1e-1 1e-1 1e-1 ]</td></tr><tr><td>t<sub>f</sub> =</td><td>5e1 s</td><td>p̃<sub>0</sub> =</td><td>[ -r<sub>d</sub> -r<sub>d</sub> 0 ]<sup>T</sup> m</td></tr><tr><td>t<sub>0</sub> =</td><td>0 s</td><td>dot̃<sub>0</sub> =</td><td>[ 0 0 0 ]<sup>T</sup> m/s</td></tr><tr><td>r<sub>d</sub> =</td><td>1.5 m</td><td>η<sub>0</sub> =</td><td>[ 0 0 0 ]<sup>T</sup> rad</td></tr><tr><td>z<sub>d</sub> =</td><td>1.5 m</td><td>K<sub>z</sub> =</td><td>[ 1e0 3.5e1 1e2 ]</td></tr><tr><td>c<sub>τ<sub>t</sub></sub> =</td><td>-1e0</td><td>g̃ =</td><td>[ 0 0 -9.80665 ]<sup>T</sup> m/s<sup>2</sup></td></tr><tr><td>c<sub>τ<sub>r</sub></sub> =</td><td>-1e0</td><td>K<sub>xy</sub> =</td><td>[ 0.2 0.4 0.5 ]</td></tr></table>

<figure><img src="image_10.png" /><figcaption>Fig. 8. Simulation results for the altitude, speed and the control inputs</figcaption></figure>

Figure 8 shows that the altitude system is asymptotically stable towards $z_d$. As the altitude system is able to converge rapidly, the pitch remains relatively constant for the rest of the simulation. The ornithopter's angular speed $\omega_z$ stabilizes
towards a constant value, and this can be seen in the oscillating yaw angle. The Figure also shows how the controllers stabilize asymptotically towards constant values.

<figure><img src="image_11.png" /><figcaption>Fig. 9. Simulation results for the position in the x-y plane</figcaption></figure>

Figure 9 is a representation of the translational movements of the X-Fly. Oscillations can be seen at the start, but eventually the system stabilizes close to the desired radius, shown in red.

## B. Experimental Results

The experiments were performed in an indoor flight arena equipped with the Vicon motion capture system. The parameters used for the controllers are shown in Table II.

<table><caption>TABLE II
PARAMETERS FOR THE EXPERIMENTS</caption><tr><td>k<sub>z,i</sub> = 1e<sup>-2</sup></td><td>k<sub>z,p</sub> = 8e<sup>-1</sup></td><td>k<sub>z,d</sub> = 1e<sup>-1</sup></td></tr><tr><td>K<sub>xy,i</sub> = 1e<sup>-4</sup></td><td>K<sub>xy,Δψ</sub> = 5e<sup>-1</sup></td><td>K<sub>xy,ωz</sub> = 1e<sup>-4</sup></td></tr></table>

<figure><img src="image_12.png" /><figcaption>Fig. 10. Experimental results for the altitude, orientation and control inputs</figcaption></figure>

Figure 10 shows how the yaw of the drone changes as it tracks the circle trajectory. The roll and pitch angles are shown to be stable, thus experimentally corroborating the stability of the roll inner controller. The altitude controller

as depicted in Theorem 1 is also stable, but the oscillations due to the flapping movement, the unmodelled dynamics and other simplifications are more evident in this plot. There is a considerable altitude error. This is because the speed of the vehicle drastically changes with very close turns, which does not give enough time for the altitude controller to compensate in this case due to coupled dynamics between both systems. However, the controller is still able to keep the altitude stable.

Figure 10 also shows the control commands given to the X-Fly during the experiment. Unlike the commands shown on Figure 8, the experimental motor commands indicate a Pulse Width Modulation (PWM) percentage from 0 to 1. The negative values in $u_r$ indicate a direction change in the current and angle of the rudder. A rate limiter was implemented on the controller to prevent the value of $u_m$ to change more than 20% between consecutive discrete commands. This helped to prevent shattering effects due to the flapping, and to smooth out the control action.

<figure><img src="image_13.png" /><figcaption>Fig. 11. Experimental results for the position in the x-y plane</figcaption></figure>

Figure 11 presents the position of the artificial bird on the x-y plane. It can be seen that it oscillates around the desired radius, which corroborates the stability as predicted in Figure 9.

A short video showing the X-Fly using the presented control strategy can be found in: https://youtu.be/G1-csBxXKZU

## V. CONCLUSIONS

This work has introduced the modeling and analysis of a new FWMAV with a focus on control. It was proven that the proposed separation of the dynamical model into an altitude system and a translational system was possible as seen in equations (15), (16), (17), (19), (20) and (21).

Furthermore, a controller for $u_m$ was proposed in Theorem 1 that was able to maintain a constant altitude, but also prove and ensure the stability of the pitch dynamics. A similar procedure was made for the control input $u_r$ in Theorem 3. Both of these models and theorems were corroborated through simulations and experimental tests.

The presented model was able to give a more detailed description of the system. However, the tracking could be improved by means of an integral action on the translational system.
It was noted that during the experiments the performance of the controllers was greatly affected by the battery level. It would be worth considering how this information could be integrated into the model.

## ACKNOWLEDGMENTS

We thank Nicolas Baker from the content agency “Animal pensant” and CNRS Images for the videos. We also thank Abdoullah Ndoye for his kind help and for the discussions.

## REFERENCES

[1] A. Hammad and S. F. Armanini, “Landing and takeoff capabilities of bioinspired aerial vehicles: A review,” *Bioinspiration & Biomimetics*, vol. 19, no. 3, p. 031001, Mar. 2024.

[2] M. Keennon, K. Klingebiel, and H. Won, “Development of the nano hummingbird: A tailless flapping wing micro air vehicle,” in *50th AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition*.

[3] K. Y. Ma, P. Chirarattananon, S. B. Fuller, and R. J. Wood, “Controlled flight of a biologically inspired, insect-scale robot,” *Science*, vol. 340, no. 6132, pp. 603–607, 2013.

[4] M. Karásek, F. T. Muijres, C. D. Wagter, B. D. W. Remes, and G. C. H. E. de Croon, “A tailless aerial robotic flapper reveals that flies use torque coupling in rapid banked turns,” *Science*, vol. 361, no. 6407, pp. 1089–1094, 2018.

[5] A. Chen, B. Song, Z. Wang, D. Xue, and K. Liu, “A novel actuation strategy for an agile bioinspired fwing performing a morphing-coupled wingbeat pattern,” *IEEE Transactions on Robotics*, vol. 39, no. 1, pp. 452–469, 2023.

[6] S. S. Baek and R. S. Fearing, “Flight forces and altitude regulation of 12 gram i-bird,” in *2010 3rd IEEE RAS & EMBS International Conference on Biomedical Robotics and Biomechatronics*, 2010, pp. 454–460.

[7] W. He, H. Huang, Y. Chen, et al., “Development of an autonomous flapping-wing aerial vehicle,” *Science China Information Sciences*, vol. 60, no. 6, May 2017.

[8] G. C. H. E. de Croon, M. A. Groen, C. De Wagter, B. Remes, R. Ruijsink, and B. W. van Oudheusden, “Design, aerodynamics and autonomy of the delfly,” *Bioinspiration & Biomimetics*, vol. 7, no. 2, p. 025003, May 2012.

[9] C. De Wagter, S. Tijmons, B. D. W. Remes, and G. C. H. E. de Croon, “Autonomous flight of a 20-gram flapping wing mav with a 4-gram onboard stereo vision system,” in *2014 IEEE International Conference on Robotics and Automation (ICRA)*, 2014, pp. 4982–4987.

[10] J. Oh, S. Kim, B. Lee, S. Kim, and J. Suk, “System identification & attitude control of avian-type flyer with flight test,” in *2017 11th Asian Control Conference (ASCC)*, 2017, pp. 1677–1682.

[11] J. Hoff, U. Syed, A. Ramezani, and S. Hutchinson, “Trajectory planning for a bat-like flapping wing robot,” in *2019 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)*, 2019, pp. 6800–6805.

[12] A. Ndoye, J. J. Castillo-Zamora, S. Samorah-Laki, R. Miot, E. Van Ruymbeke, and F. Ruffier, “Vector field aided trajectory tracking by a 10-gram flapping-wing micro aerial vehicle,” in *2023 IEEE International Conference on Robotics and Automation (ICRA)*, 2023, pp. 5379–5385.

[13] C. Qian, R. Chen, P. Shen, Y. Fang, J. Yan, and T. Li, “Trajectory generation and tracking control for flapping wing robot three-dimensional flight,” *IEEE/ASME Transactions on Mechatronics*, pp. 1–13, 2024.

[14] E. Van Ruymbeke, *Flying toy able to move by the flapping of wings*, US Patent 8,382,546, Feb. 2013.

[15] A. Bakhtiari, S. E. Haghighi, and A. Maghsoudpour, “Modeling and control of a flapping wing robot,” *Proceedings of the Institution of Mechanical Engineers, Part K: Journal of Multi-body Dynamics*, vol. 233, no. 1, pp. 174–181, 2019.

[16] A. Banazadeh and N. Taymourtash, “Adaptive attitude and position control of an insect-like flapping wing air vehicle,” *Nonlinear Dynamics*, vol. 85, no. 1, pp. 47–66, Jul. 2016.

[17] S. F. Armanini, C. C. de Visser, and G. de Croon, “Black-box lti modelling of flapping-wing micro aerial vehicle dynamics,” in *AIAA Atmospheric Flight Mechanics Conference*. 2015.

[18] H. Gim, S. Kim, J. Suk, and S. Cho, “Longitudinal system identification of ornithopter with automated flight tests**this research was supported by a grant to bio-mimetic robot research center funded by defense acquisition program administration (ud130070id).,” *IFAC-PapersOnLine*, vol. 49, no. 17, pp. 194–199, 2016, 20th IFAC Symposium on Automatic Control in AerospaceACA 2016, ISSN: 2405-8963.

[19] J. Grauer, E. Ulrich, J. Hubbard, D. Pines, and J. S. Humbert, “Testing and system identification of an ornithopter in longitudinal flight,” *Journal of Aircraft*, vol. 48, no. 2, pp. 660–667, 2011.

[20] A. Ndoye, S. S. Laki, R. Miot, E. van Ruymbeke, and F. Ruffier, “Ornithopter’s intra-flapping body pitch highly depends on wingbeat frequency,” *Computer Methods in Biomechanics and Biomedical Engineering*, vol. 24, no. Sup 1, S101–S103, 2021.

[21] W. He, X. Mu, L. Zhang, and Y. Zou, “Modeling and trajectory tracking control for flapping-wing micro aerial vehicles,” *IEEE/CAA Journal of Automatica Sinica*, vol. 8, no. 1, pp. 148–156, 2021.

[22] D. Jordan and P. Smith, *Nonlinear ordinary differential equations: an introduction for scientists and engineers*. Oxford University Press, 2007.