AERO-LQG: Aerial-Enabled Robust Optimization for LQG-Based Quadrotor Flight Controller

TL;DR

AERO-LQG introduces an evolutionary strategy-based framework for tuning LQG parameters in quadrotor control, significantly improving energy efficiency and robustness in mission-specific flight modes.

Contribution

It presents a novel aerial-enabled robust optimization method for LQG tuning using evolutionary strategies, addressing the challenge of weight selection in optimal control.

Findings

Performance improvements of several tens of percent in simulated hover mode.

Demonstrates the effectiveness of evolutionary strategies for control parameter optimization.

Potential for enhanced energy efficiency and robustness in quadrotor flight control.

Abstract

Quadrotors are indispensable in civilian, industrial, and military domains, undertaking complex, high-precision tasks once reserved for specialized systems. Across all contexts, energy efficiency remains a critical constraint: quadrotors must reconcile the high power demands of agility with the minimal consumption required for extended endurance. Meeting this trade-off calls for mode-specific optimization frameworks that adapt to diverse mission profiles. At their core lie optimal control policies defining error functions whose minimization yields robust, mission-tailored performance. While solutions are straightforward for fixed weight matrices, selecting those weights is a far greater challenge-lacking analytical guidance and thus relying on exhaustive or stochastic search. This interdependence can be framed as a bi-level optimization problem, with the outer loop determining weights a priori. This work introduces an aerial-enabled robust optimization for LQG tuning (AERO-LQG), a framework employing evolutionary strategy to fine-tune LQG weighting parameters. Applied to the linearized hovering mode of quadrotor flight, AERO-LQG achieves performance gains of several tens of percent, underscoring its potential for enabling high-performance, energy-efficient quadrotor control. The project is available at GitHub.

Tables

Table 1. TABLE I : Comparative analysis of weight optimization methods, ordered top-to-bottom by their novelty and precision.

Symbol	Optimization method	Approach	Pros	Cons	Complexity
MT	Manual tuning	Trial-and-error	Trivial implementation	Inefficient (exhaustive)	$O\left(nf\right)$
BR	Bryson’s rule [28]	Heuristic	Intuitive; Fast initial guess	Inherently suboptimal	$O\left(f\right)$
PS	Particle swarm [29]	Swarm-based	Robust against local minima	Premature convergence	$O\left(knf\right)$
GA	Genetic algorithm [30]	Evolutionary	Global search; gradient-free	Computationally intensive	$O\left(knf\right)$
BY	Bayesian [31]	Surrogate-based	Sample-efficient	Scalability issues	$O\left(kf+k^3\right)$
CMA	Covariance matrix adaptation [32]	Simulation-based	Robust in rugged search spaces	Slower convergence	$O\left(knf\right)$

Table 2. TABLE II : Comparison of optimizers by final-time errors.

	Position [m]		Orientation [deg]		Effort [N$\cdot$s]
	$\Vert {𝒆}_{{𝝃}^{ℐ}}\Vert$	$\Vert {𝜺}_{{𝝃}^{ℐ}}\Vert$	$\Vert {𝒆}_{{𝜼}^{ℐ}}\Vert$	$\Vert {𝜺}_{{𝜼}^{ℐ}}\Vert$	${\int }_0^T𝒖\left(t\right)dt$
MT	$\to \infty$	$\to \infty$	$\to \infty$	$\to \infty$	$\to \infty$
BR	0.48	1.12	0.57	12.3	$>{10}^3$
PS	0.22	0.77	0.18	6.13	293.8
GA	0.03	0.09	0.22	2.05	119.1
BY	0.31	0.27	0.47	3.21	143.1
CMA	0.01	0.05	0.17	1.74	98.6

Equations

$$\dot{\mathbf{x}}(t)=\mathbf{f}(\mathbf{x}(t),\mathbf{u}(t)),$$

$$\dot{\mathbf{x}}(t)=\mathbf{f}({\mathbf{x}}_e,{\mathbf{u}}_e)=0.$$

$$\dot{\mathbf{x}}(t)\approx \mathbf{A}({\mathbf{x}}_e,{\mathbf{u}}_e)(\mathbf{x}(t)-{\mathbf{x}}_e)+\mathbf{B}({\mathbf{x}}_e,{\mathbf{u}}_e)(\mathbf{u}(t)-{\mathbf{u}}_e),$$

$$\delta \mathbf{x}(t)=\mathbf{x}(t)-{\mathbf{x}}_e\text{and}\delta \mathbf{u}(t)=\mathbf{u}(t)-{\mathbf{u}}_e.$$

$$\displaystyle\mathbf{A}(\bm{x}_{e},\bm{u}_{e})=\frac{\partial\bm{f}}{\partial\bm{x}}\bigg{|}_{(\bm{x}_{e},\bm{u}_{e})}\,\text{and}\ \,\mathbf{B}(\bm{x}_{e},\bm{u}_{e})=\frac{\partial\bm{f}}{\partial\bm{u}}\bigg{|}_{(\bm{x}_{e},\bm{u}_{e})},$$

$$\dot{\mathbf{x}}(t)$$

$$\mathbf{y}(t)$$

$$\mathbf{w}(t)$$

$$\mathbf{v}(t)$$

$$\dot{\mathbf{P}}(t)=\mathbf{AP}(t)+\mathbf{P}(t){\mathbf{A}}^{\top }+\mathbf{W}.$$

$$\left[\begin{array}{c}\mathbf{x}(t)\\ \mathbf{y}(t)\end{array}\right]\sim \mathcal{N}\left(\left[\begin{array}{c}\mathbf{\mu }(t)\\ \mathbf{C}\mathbf{\mu }(t)\end{array}\right],\left[\begin{array}{cc}\mathbf{P}(t)& \mathbf{P}(t){\mathbf{C}}^{\top }\\ \mathbf{CP}(t)& \mathbf{CP}(t){\mathbf{C}}^{\top }+\mathbf{V}\end{array}\right]\right).$$

$$\dot{\mathbf{\mu }}(t)$$

$$\underset{\mathbf{u}(t)}{\min }{\mathcal{J}}_{\text{in}}\triangleq \underset{T\to \infty }{\lim }\mathbb{E}\left[\frac{1}{T}{\int }_0^T(\Vert \mathbf{x}(t){\Vert }_{\mathbf{Q}}^2+\Vert \mathbf{u}(t){\Vert }_{\mathbf{R}}^2)\mathrm{d}t\right],$$

$$\mathbf{e}(t)$$

$$\hat{\mathbf{\epsilon }}(t)$$

$$\left[\begin{array}{c}\dot{\mathbf{x}}\\ \dot{\mathbf{e}}\end{array}\right]=\left[\begin{array}{cc}\mathbf{A}-\mathbf{BK}& \mathbf{BK}\\ 0& \mathbf{A}-\mathbf{LC}\end{array}\right]\left[\begin{array}{c}\mathbf{x}\\ \mathbf{e}\end{array}\right]+\left[\begin{array}{cc}\mathbf{I}& 0\\ \mathbf{I}& -\mathbf{L}\end{array}\right]\left[\begin{array}{c}\mathbf{w}\\ \mathbf{v}\end{array}\right],$$

$${\mathcal{R}}_{\mathcal{B}}^{\mathcal{I}}:=\left({\mathbf{e}}_x^{\mathcal{B}},{\mathbf{e}}_y^{\mathcal{B}},{\mathbf{e}}_z^{\mathcal{B}}\right)\in SO(3),$$

$$\left(\begin{array}{c}{\mathbf{f}}^{\mathcal{B}}\\ {\mathbf{\tau }}^{\mathcal{B}}\end{array}\right)=\left(\begin{array}{cc}m{\mathbf{I}}_3& 0_3\\ 0_3& \mathbf{J}\end{array}\right)\left(\begin{array}{c}{\dot{\mathbf{v}}}^{\mathcal{B}}\\ {\dot{\mathbf{\omega }}}^{\mathcal{B}}\end{array}\right)+\left(\begin{array}{c}{\mathbf{\omega }}^{\mathcal{B}}\times m{\mathbf{v}}^{\mathcal{B}}\\ {\mathbf{\omega }}^{\mathcal{B}}\times \mathbf{J}{\mathbf{\omega }}^{\mathcal{B}}\end{array}\right),$$

$${\mathbf{f}}^{\mathcal{B}}=\sum _{i=1}^4{\mathbf{f}}_{z,i}^{\mathcal{B}}-m{\text{g}}^{\mathcal{B}}.$$

$${\mathbf{f}}_{z,i}^{\mathcal{B}}=f_{z,i}{\text{e}}_z^{\mathcal{B}}=\left(k_T{\Omega }_i^2\right){\text{e}}_z^{\mathcal{B}},$$

$${\mathbf{\tau }}_{\varphi }^{\mathcal{B}}=\ell \left({\mathbf{f}}_{z,4}^{\mathcal{B}}-{\mathbf{f}}_{z,2}^{\mathcal{B}}\right){\text{e}}_x^{\mathcal{B}},$$

$${\mathbf{\tau }}_{\theta }^{\mathcal{B}}=\ell \left({\mathbf{f}}_{z,3}^{\mathcal{B}}-{\mathbf{f}}_{z,1}^{\mathcal{B}}\right){\text{e}}_y^{\mathcal{B}}.$$

$${\mathbf{\tau }}_{\psi }^{\mathcal{B}}=\sum _{i=1}^4{\tau }_{M,i}=(-1{)}^i\left(k_M{\Omega }_i^2\right){\text{e}}_z^{\mathcal{B}}.$$

$${\mathbf{\xi }}^{\mathcal{I}}={\left(x,y,z\right)}^{\top },$$

$${\dot{\mathbf{\xi }}}^{\mathcal{I}}={\mathcal{R}}_{\mathcal{B}}^{\mathcal{I}}{\mathbf{v}}^{\mathcal{B}}.$$

$${\mathbf{\eta }}^{\mathcal{I}}=(\varphi ,\theta ,\psi {)}^{\top },$$

$${\dot{\mathbf{\eta }}}^{\mathcal{I}}={\text{W}}_{\mathbf{\eta }}{\mathbf{\omega }}^{\mathcal{B}}=\left(\begin{array}{ccc}1& {\text{s}}_{\varphi }{\text{t}}_{\theta }& {\text{c}}_{\varphi }{\text{t}}_{\theta }\\ 0& {\text{c}}_{\varphi }& -{\text{s}}_{\varphi }\\ 0& {\text{s}}_{\varphi }/{\text{c}}_{\theta }& {\text{c}}_{\varphi }/{\text{c}}_{\theta }\end{array}\right)\left(\begin{array}{c}p\\ q\\ r\end{array}\right),$$

$$\mathbf{u}(\Omega )=\left(\begin{array}{c}{\mathbf{u}}_z\\ {\mathbf{u}}_{\varphi }\\ {\mathbf{u}}_{\theta }\\ {\mathbf{u}}_{\psi }\end{array}\right)=\left(\begin{array}{c}{f}_z{\text{e}}_z^{\mathcal{B}}\\ {\tau }_{\varphi }{\text{e}}_x^{\mathcal{B}}\\ {\tau }_{\theta }{\text{e}}_y^{\mathcal{B}}\\ {\tau }_{\psi }{\text{e}}_z^{\mathcal{B}}\end{array}\right)\in {\mathbb{R}}^4.$$

$$\mathbf{x}={\left(\begin{array}{cccc}({\mathbf{\xi }}^{\mathcal{I}}{)}^{\top }& ({\mathbf{v}}^{\mathcal{B}}{)}^{\top }& ({\mathbf{\eta }}^{\mathcal{I}}{)}^{\top }& ({\mathbf{\omega }}^{\mathcal{B}}{)}^{\top }\end{array}\right)}^{\top }\in {\mathbb{R}}^{12},$$

$$\left(\begin{array}{c}{\dot{\mathbf{\xi }}}^{\mathcal{I}}\\ {\dot{\mathbf{v}}}^{\mathcal{B}}\\ {\dot{\mathbf{\eta }}}^{\mathcal{I}}\\ {\dot{\mathbf{\omega }}}^{\mathcal{B}}\end{array}\right)=\left(\begin{array}{c}({\mathcal{R}}_{\mathcal{B}}^{\mathcal{I}}){\mathbf{v}}^{\mathcal{B}}\\ m^{-1}{\mathbf{u}}_z-{\mathbf{\omega }}^{\mathcal{B}}\times {\mathbf{v}}^{\mathcal{B}}-\text{g}{\text{e}}_z^{\mathcal{B}}\\ ({\text{W}}_{\mathbf{\eta }}){\mathbf{\omega }}^{\mathcal{B}}\\ {\mathbf{J}}^{-1}\left(({\mathbf{u}}_{\varphi }+{\mathbf{u}}_{\theta }+{\mathbf{u}}_{\psi })-{\mathbf{\omega }}^{\mathcal{B}}\times \mathbf{J}{\mathbf{\omega }}^{\mathcal{B}}\right)\end{array}\right).$$

$$\frac{\partial \mathbf{f}}{\partial \mathbf{x}}=\left(\begin{array}{cccc}{0}_3& {\mathbf{I}}_3+(\mathbf{\eta }{)}_{\times }& -({\mathbf{v}}^{\mathcal{B}}{)}_{\times }& 0_3\\ 0_3& -({\mathbf{\omega }}^{\mathcal{B}}{)}_{\times }& -({\text{g}}^{\mathcal{I}}{)}_{\times }& ({\mathbf{v}}^{\mathcal{B}}{)}_{\times }\\ 0_3& 0_3& 0_3& {\mathbf{I}}_3\\ 0_3& 0_3& 0_3& {\nabla }_{\mathbf{\omega }}({\dot{\mathbf{\omega }}}^{\mathcal{B}})\end{array}\right),$$

$${\nabla }_{\mathbf{\omega }}({\dot{\mathbf{\omega }}}^{\mathcal{B}})={\mathbf{J}}^{-1}\left((\mathbf{J}{\mathbf{\omega }}^{\mathcal{B}}{)}_{\times }-{\mathbf{\omega }}^{\mathcal{B}}\times \mathbf{J}\right).$$

Full text

I Introduction

Unmanned aerial vehicles (UAVs), particularly quadrotors, have rapidly evolved from research prototypes into practical tools across a broad range of applications. Their versatility and agility enable complex missions—such as surveillance [1], inspection [2], delivery [3], and search and rescue [4]—making them valuable assets in both civilian and military domains [5]. From a dynamical systems perspective, quadrotor flight modes span a wide stability spectrum—from near-equilibrium hovering to highly nonlinear, aggressive maneuvers. As a result, no single control strategy can achieve optimal performance across all regimes. Transitions between control modes are therefore common, especially as energy constraints impose competing demands: controllers must either deliver agility for high-intensity missions [6], or conversely, limit maneuverability to conserve energy and extend flight time [7, 8].

Monolithic control lacks such flexibility; adaptive scheduling with task-dependent adjustments is needed [9, 10, 11]. Frequent mode switching, however, complicates state estimation, which remains fragile to the quadrotor’s nonlinear and underactuated dynamics—despite advances in hardware [12], inertial sensing [13], and AI-driven perception [14]. To mitigate these issues, linear quadratic Gaussian (LQG) control remains a staple in aerial robotics, valued for its theoretical rigor and optimality under specific assumptions [15].

Many stabilized flight modes can be approximated as linear systems, allowing precise LQG control, provided operating conditions stay near the linearization point and modeling uncertainties are managed [16]. Quadrotors are no exception—their hovering mode, in particular, can be effectively linearized to facilitate both control and state estimation. However, while the quadratic form of the objective cost function is inherently convex for fixed weights, tuning these values constitutes an additional outer optimization task—highly non-convex due to the coupled estimator–controller dynamics [17]. This raises a key question:

Could bi-level optimization redefine the limits of LQG tuning for quadrotor stability ?

To answer this question, we develop an aerial-enabled robust optimization for LQG tuning (AERO-LQG), tailored to modern quadrotors. As illustrated in Fig. 1, an outer evolutionary loop refines the weighting matrices, while an inner loop computes the LQG gains, together forming a nested architecture. Evaluated across diverse metrics, AERO-LQG delivers unprecedented gains—reducing tracking error by over 55% and extending flight endurance by nearly 8%.

The paper unfolds as follows: Section II presents the problem setup and technical foundations; Section III outlines the proposed approach; Section IV discusses the findings; and Section V concludes the study.

II Problem Setup

Quadrotor behavior is modeled as a continuous-time dynamical system, defined by the following differential equation

[TABLE]

where $\bm{x}(t)\in\mathbb{R}^{n}$ is the state vector, $\bm{u}(t)\in\mathbb{R}^{k}$ is the control input vector, and $\bm{f}:\mathbb{R}^{n}\times\mathbb{R}^{k}\to\mathbb{R}^{n}$ is a smooth function that defines its dynamics. The solution to this system describes the trajectory of the state $\bm{x}(t)$, starting from an initial condition $\bm{x}(0)$ and evolving under the influence of the input $\bm{u}(t)$.

Most aerial flight modes such as take-off, uncoordinated turns, maneuvers, and landing, require the system states to change dramatically, i.e. $\dot{\bm{x}}(t)\neq\mathbf{0}$, thereby causing $\bm{f}$ to evolve nonlinearly in time. That said, certain operating points, if controlled properly, enable maintaining the system at an equilibrium, denoted by subscript $e$, with their states demonstrating negligible change [18], such that

[TABLE]

The equilibrium behavior is assessed using local linearity, approximating the nonlinear dynamics $\bm{f}$ with a linear model near the equilibrium point $(\bm{x}_{e},\bm{u}_{e})$. General fixed-wing configurations permit linear approximation of several trimmed (stable) flight modes—such as straight-and-level flight, coordinated turns, and steady climbs or descents—thanks to the gliding capability of their large wing area. Rotorcraft, by contrast, depend entirely on rapid, coordinated motor corrections, yielding predominantly nonlinear dynamics [17, 19]. Consequently, the only flight mode that consistently maintains equilibrium under small perturbations is hovering (center), as depicted in the finite state machine in Fig. 2.

II-A Equilibrium Point

Based on linear systems theory, steady-state dynamics allow the nonlinear system (1) to be linearized around the equilibrium point $(\bm{x}_{e},\bm{u}_{e})$ through a first-order Taylor expansion [20]. Since by definition, quadrotor hovering balances the net forces and net moments into zero, the following approximation holds around the hovering equilibrium point

[TABLE]

where $\bm{x}(t)$ is the state vector, $\bm{u}(t)$ the control inputs, and the deviations from equilibrium are expressed as

[TABLE]

Assuming linear time-invariant (LTI) dynamics in the vicinity of the equilibrium point, the system Jacobians with respect to the state and input are defined as

[TABLE]

reducing the error dynamics to the standard LTI form

[TABLE]

where matrix $\mathbf{C}$ extracts the measurable states from $\bm{x}$. Model discrepancies and unmodeled dynamics are represented by zero-mean white Gaussian noise processes, with process covariance $\mathbf{W}$ and measurement covariance $\mathbf{V}$

[TABLE]

These propagate the open-loop uncertainty, whose state covariance ${\mathbf{P}}(t)$ evolves by the continuous-time Lyapunov equation

[TABLE]

At any fixed time $t$, the joint distribution of the states and outputs can thus be approximated as Gaussian, namely

[TABLE]

where $\bm{\mu}(t)$ is the nominal noise-free state trajectory, initialized as $\bm{\mu}(t_{0})=\mathbb{E}[\bm{x}(t_{0})]$, with time evolution given by

[TABLE]

II-B Inner-Level Optimization

Having identified the hovering operating point as locally linear under small perturbations enables the application of linear control theory. In this framework, solving the Riccati equations yields optimal controller ($\mathbf{K}$) and estimator ($\mathbf{L}$) gains, whereby the Kalman filter provides state estimates to the linear quadratic regulator (LQR). These gains are obtained by minimizing a quadratic cost over a finite time horizon $T$,

[TABLE]

subject to (11), with a-priori specified $\mathbf{Q}\succeq 0$ and $\mathbf{R}\succ 0$ weighting deviations in $\bm{x}(t)$ and $\bm{u}(t)$, respectively [21, 16]. By defining the estimation and the control errors as

[TABLE]

the plant-estimator dynamics admits the block-state form of

[TABLE]

highlighting the separation principle, with an upper-triangular form that decouples the estimator dynamics [22, 23, 24].

II-C System Dynamics

We begin by outlining the assumptions and coordinate frames relevant to hovering, a near-equilibrium condition that justifies standard modeling practices [25, 26]: (i) the quadrotor is a rigid symmetric body; (ii) thrust and moments scale with the square of rotor speeds; (iii) aerodynamic drag is neglected due to low speed; (iv) the center of mass coincides with the body-frame origin. The orientation of the body frame $\mathcal{B}$ relative to the inertial frame $\mathcal{I}$ is given by the rotation matrix

[TABLE]

which is constructed through a sequence of three consecutive rotations: $\mathcal{R}_{z}(\psi)$ about the $\mathbf{e}_{z}^{\mathcal{I}}$ axis, $\mathcal{R}_{y}(\theta)$ about the $\mathbf{e}_{y}^{\mathcal{I}}$ axis, and $\mathcal{R}_{x}(\phi)$ about the $\mathbf{e}_{x}^{\mathcal{I}}$ axis, where $\{\phi,\theta\}\in\left(-\frac{\pi}{2},\frac{\pi}{2}\right)$ and $\psi\in\left(-\pi,\pi\right)$ correspond to the roll, pitch, and yaw angles, respectively. Building on the Newton-Euler formalism [27], the coupled translational and rotational dynamics are expressed as

[TABLE]

where, $m$ is the body mass; $\bm{J}$, $\bm{I}_{3}$, and $\bm{0}_{3}$ are $3\times 3$ inertia tensor, identity, and zero matrices. Forces $\bm{f}$ and moments $\bm{\tau}$, expressed in the body frame, govern the time derivatives of $\bm{v}^{\mathcal{B}}$ and $\bm{\omega}^{\mathcal{B}}$. In fixed-point hover, near-zero body velocities make drag negligible, so the net body-frame force includes only gravity $\mathbf{g}^{\mathcal{B}}=\text{g}\,\mathbf{e}_{z}^{\mathcal{B}}$ and vertical rotor thrusts $\bm{f}_{z,i}^{\mathcal{B}}$, yielding

[TABLE]

Each rotor spins along the body-frame $z$-axis, generating thrust by the $i$-th propeller modeled as

[TABLE]

where $k_{T}$ is a lumped coefficient. Attitude stabilization is achieved through differential thrust between opposing rotors spaced by $\ell$. The rolling moment arises from the thrust difference between motors #4 (left) and #2 (right), given by

[TABLE]

while the pitch moment is controlled by the longitudinal thrust difference between motors #1 (front) and #3 (rear), as

[TABLE]

Finally, yaw arises from asymmetry in rotor drag torques, i.e., when $\sum_{i}\tau_{M,i}\neq 0$. Defining $k_{M}$ as the torque coefficient and positive yaw as counter-clockwise about $\textbf{e}_{z}^{\mathcal{B}}$, the yaw torque is

[TABLE]

Fig. 3 depicts the quadrotor free-body diagram: airframe (thick brown), reference axes (black), rotor thrust vectors (green), control inputs (azure), and rotor spin actuation (orange); gravity and inertial displacement are shown as dashed guides.

II-D System Kinematics

The platform’s trajectory is derived by time-integrating its dynamics and applying a coordinate transformation to an inertial reference frame, letting the position vector be [15]

[TABLE]

with the corresponding body velocities related by

[TABLE]

The Euler angles vector in the inertial frame is denoted as

[TABLE]

with its rate of change related to the body angular rates $(p,q,r)$ through the transformation matrix $\textbf{{W}}_{\bm{\eta}}$, such that

[TABLE]

where ’t’ is the tangent function. For derivation, thrust and torques are idealized as direct control input functions, hence

[TABLE]

Finally, by incorporating the generalized coordinates into the system model, the augmented state vector is defined as

[TABLE]

with the nonlinear state-space model in (1) evolves as

[TABLE]

III Proposed Approach

We begin by analyzing the linearized dynamics about the hover equilibrium, outlining the outer evolutionary LQG optimization process, and lastly presenting our proposed implementation amongst several state-of-the-art tuning methods.

III-A Hovering Linearization

Since hovering occurs near equilibrium with small state variables, local behavior can be approximated using small-angle assumptions. This simplifies (26) to $\dot{\bm{\eta}}^{\mathcal{I}}\approx\bm{\omega}^{\mathcal{B}}$ and $\mathcal{R}_{\mathcal{B}}^{\mathcal{I}}\approx\bm{I}_{3}+(\bm{\eta})_{\times}$, respectively, where $(\cdot)_{\times}$ denotes the skew-symmetric operator. With these approximations, the Jacobian of the nonlinear dynamics in (29) is computed as

[TABLE]

with the corresponding sub-Jacobian of $\dot{\bm{\omega}}^{\mathcal{B}}$ given by

[TABLE]

Assuming a hovering equilibrium defined by a state vector at a fixed point and a nominal thrust control input

[TABLE]

substituting this operating point into the state Jacobian yields the linearized system dynamics model

[TABLE]

where higher-order centrifugal terms and their Coriolis couplings vanish. Likewise, linearizing the control input around the thrust required to balance gravity gives

[TABLE]

Within the state estimation process, these linearized dynamics hold locally around the hover equilibrium $(\bm{x}_{e},\bm{u}_{e})$, enabling state propagation via the state-space model (6).

III-B Outer-Level Optimization

Gradient-based methods, while effective in high-dimensional settings, are unsuitable here since analytical gradients are unavailable: the cost depends implicitly on closed-loop dynamics. The LQG cost (12) is set a priori by $\mathbf{Q}$ and $\mathbf{R}$, whose iterative tuning forms an inner feedback loop within the outer optimizer. As hovering prioritizes position holding, followed by attitude and yaw stabilization, the outer cost is defined as

[TABLE]

where $\lambda\ll 1$ penalizes the translational–rotational error trade-off. For an $n$-dimensional state-space model (6), the search spans $2n^{2}$ weight parameters, placing the bi-level tuning firmly in the domain of black-box optimization, where weights are iteratively refined without explicit gradients or access to internal system structure.

Remark 1: To enforce positive definiteness, both weight matrices are parameterized via Cholesky factorization, $\mathbf{Q}=\bm{L}_{\text{Q}}\bm{L}_{\text{Q}}^{\top}$, and $\mathbf{R}=\bm{L}_{\text{R}}\bm{L}_{\text{R}}^{\top}$, with $\bm{L}$ being lower triangular, effectively reducing the parameter count to $n(n+1)$.

Algorithm 1 formalizes the proposed framework in a high-level, optimizer-agnostic form, with the inner LQG control loop nested within the outer weight-optimization process.

Table I lists the six optimization techniques used as benchmarks in this study, distinguished by their key characteristics and ranked in ascending order of complexity. To assess the effectiveness of our covariance matrix adaptation (CMA) implementation—highlighted in green at the bottom—complexity is expressed in terms of number of parameters ($n$), iterations ($k$), and cost per function evaluation ($f$).

IV Analysis and Results

To explore the intricacies of LQG tuning in quadrotor hovering, AERO-LQG results are organized into three subsections.

IV-A Cost landscape

When applied to a 12-dimensional state space rich in cross-coupled dynamics, the LQG quadratic cost (12) yields a highly non-convex landscape over the gain parameters. In such settings, gradient-based optimizers are prone to becoming trapped in inferior local minima, with little prospect of escape. To demonstrate, we compute the cost of a 10-second standardized gust trajectory, assuming isotropic 3D error and using L2 norms for $\mathbf{Q}_{\xi}$ and $\mathbf{R}_{\xi}$ for clarity. As shown in Fig. 4, the scalar surface is riddled with local minima, exposing the shortcomings of gradient-based techniques.

This, in turn, accentuates the challenge of proper LQG tuning: large regions of the cost surface correspond to poor controller performance, with log-scale costs on the order of $\mathcal{J}_{\text{out}}\gtrsim 10^{5}$. Yet, precise calibration of $\mathbf{Q}_{\xi}$ and $\mathbf{R}_{\xi}$ norms reduces the cost by nearly seven orders of magnitude, achieving near-global optimality and substantially improved performance.

IV-B LQG inherent instability

Fig. 5 illustrates the onset of instability in quadrotor states over a 10-second standardized gust trajectory [7], with the 3D error norms from (13) and (14) on the left, and all four control inputs (27) on the right. Even under careful manual tuning, the system diverges. In delicate hovering, LQG’s narrow stability margins make it highly sensitive to small disturbances or accumulated noise; once error assumptions are violated, divergence grows unchecked, and the actuators respond with rapid but ineffective corrections, oscillating between their mechanical limits.

This instability arises from the tight coupling between the estimator and controller, where even small deviations can degrade the Kalman filter, initiating a feedback loop that amplifies errors. Unlike LQR, which acts directly on measured states, LQG relies on state estimates, making it inherently more sensitive to estimation inaccuracies and tuning imperfections.

IV-C Performance Evaluation

Building on the need for robust, global optimization methods that withstand large perturbations, this section: (i) presents our CMA implementation, (ii) analyzes its performance, and (iii) benchmarks it against several state-of-the-art approaches.

Fig. 6 illustrates CMA’s iterative refinement across five stages, guiding the search distribution from random initialization (left) to convergence with the predefined criteria (right).

Normalized for convenience, both $\mathbf{Q}$ and $\mathbf{R}$ progressively concentrate their weight along the diagonal, while off-diagonal entries diminish—indicating negligible cross-correlation.

Fig. 7 follows the same top-to-bottom error analysis—position, velocity, orientation, and angular velocity—under the finalized parameters. In contrast to Fig. 5, where saturated controllers failed to restore stability, the trajectories here recover immediately from the same unit disturbance. A residual drift persists, visible as sawtooth patterns in the insets; however, the actuators remain within limits, achieving stabilization with significantly lower control effort (rightmost column).

Table II provides supporting numbers: MT at the top exhibits complete divergence ($\infty$), illustrating why manual tuning falls short next to automated algorithms. In contrast, the optimizers below perform far better, with CMA surpassing the leading GA by reducing estimation and control errors by at least 33% and 55%, respectively. Fig. 8 shows these differences across eight flight metrics from [7], benchmarking all candidate optimizers (excluding MT) with color-coded distinction.

As shown, CMA (green) attains optimality in five of eight criteria, excelling in precision and energy efficiency but lagging in agility and responsiveness. Overall, evolutionary strategies (PS, GA, CMA) consistently outperform heuristic methods (MT, BR) and Bayesian optimization by avoiding premature convergence in dense local minima. These advantages provide tangible benefits, particularly in resource-constrained applications where extended flight endurance is critical.

V Discussion and Conclusion

This work introduced AERO-LQG, a systematic weight-tuning framework for quadrotor hovering that addresses the narrow stability margins and highly non-convex cost landscapes characteristic of such systems. By decoupling the inner control loop from the outer stochastic optimizer, weight initialization is delegated to a powerful evolutionary strategy that remains robust to complex cost landscapes. Benchmarking across eight flight metrics showed the consistent superiority of our CMA implementation in key hovering capabilities—energy efficiency, flight endurance, and disturbance rejection—a result particularly significant for resource-constrained platforms. While demonstrated on an aerial system, the approach generalizes to other robots, offering resource-aware, high-performance control for real-world demands.