Evaluation of Quadrotor Flip-Robustness and Attitude Recovery Under Directional Wind Disturbances: A Comparative Study of Metaheuristic-Optimized Controllers

1. Introduction

Quadrotor unmanned aerial vehicles (UAVs) are widely used in inspection, mapping, monitoring, transportation, and autonomous-flight research because of their vertical take-off and landing, hovering, and maneuverability. Their nonlinear, multivariable, underactuated, and open-loop unstable dynamics, however, strongly couple translational and rotational motion and make closed-loop performance sensitive to modeling uncertainty, aerodynamic loads, actuator limits, and external disturbances [1,2,3]. Consequently, meaningful controller assessment should extend beyond nominal hover and isolated step responses.

No single control architecture provides the best combination of tracking accuracy, robustness, control effort, and implementation complexity under all conditions. Comparative studies report different rankings depending on the plant, trajectory, constraints, and performance metrics [4,5,6]. Fair comparison therefore requires a common nonlinear plant, identical actuator limits and disturbance scenarios, and consistent evaluation criteria.

PID control remains attractive because of its simplicity and low computational cost, although its performance depends strongly on gain selection and operating conditions. Robust-adaptive and intelligent PID variants [7,8], as well as GA- and PSO-based tuning approaches [9,10], have therefore been investigated; PSO has also been combined with sliding-mode tuning and disturbance-observer compensation [11]. LQR provides systematic state-feedback design and has demonstrated effective quadrotor stabilization and trajectory tracking [4,12], but its linearized design requires validation on the nonlinear constrained plant. Robust $H_{\infty }$ formulations explicitly address disturbance attenuation and have been applied to cascaded quadrotor control and wind rejection [13,14].

SMC is particularly attractive for rejecting matched uncertainties and disturbances. Quadrotor studies include foundational backstepping and sliding-mode designs [15], fast-terminal manifolds [16], finite- and fixed-time disturbance observers [17,18], experimentally evaluated disturbance-observer SMC [19], dynamic sliding modes [20], adaptive integral-terminal designs [21], saturation-aware adaptive SMC [22], high-order disturbance observers [23], and super-twisting or finite-time formulations [24,25]. However, differences in models, constraints, objectives, and disturbance definitions make direct cross-study ranking of PID, LQR, $H_{\infty }$, and SMC difficult.

Two comparative studies are especially relevant. Al-Qadasi et al. [26] compared GA- and PSO-tuned PID, LQR, and $H_{\infty }$ attitude controllers under nominal response, wind gusts, and actuator saturation, but excluded SMC, full six-degree-of-freedom trajectory tracking, and directional robustness analysis. Rinaldi et al. [6] compared PID, LQR, model predictive control, feedback linearization, and SMC using a common nonlinear Newton–Euler model, but did not examine cross-controller metaheuristic tuning or a complete directional disturbance sweep.

Wind disturbances can produce tracking errors, attitude excursions, increased control activity, and reduced recovery margin. Existing approaches include extended-state observers [27], state-dependent sliding gains [28], fixed-time adaptive observers under turbulent wind [18], fractional-order SMC [29], quaternion-based robust control [30], and higher-generation sliding-mode trajectory control [31]. Most evaluations, however, consider prescribed gusts, turbulence realizations, or only a few disturbance directions, providing limited information about the azimuths at which attitude, rate, recovery, or actuator-saturation limits are first violated.

This directional limitation is important because disturbance loads and their induced moments vary with orientation and application point. A controller that performs well for one disturbance direction may therefore experience substantially different actuator demands in another. Moreover, resistance to a finite translational gust does not necessarily establish attitude robustness against a directly applied disturbance moment. These mechanisms should therefore be evaluated separately while identifying the constraint responsible for any failure.

Optimization introduces an additional comparison issue. RCGA and PSO can tune nonlinear controller parameters without analytical gradients, but their outcomes depend on the objective function, parameter bounds, evaluation budget, and stochastic realization. Existing applications demonstrate the effectiveness of both methods [9,10,11,26], but do not establish universal optimizer superiority. The present study therefore employs budget-matched single runs and interprets the resulting optimizer comparison strictly as within-study evidence.

Accordingly, three gaps remain: a unified comparison of fundamentally different controllers on one constrained nonlinear plant; evaluation under sustained coupled three-dimensional tracking and finite gusts; and a systematic directional robustness protocol capable of identifying failure sectors and binding actuator constraints. This study addresses these gaps by evaluating PID, LQR, $H_{\infty }$, and SMC controllers tuned using RCGA and PSO under common plant dynamics, objectives, constraints, and simulation conditions. The evaluation includes nominal step tracking, a finite-duration translational gust, wind-disturbed helical tracking, and an attitude-only $0^{\circ }$–${360}^{\circ }$ force-azimuth sweep at ${15}^{\circ }$ resolution. In the directional study, the applied force acts through an aerodynamic center-of-pressure offset to generate a body-frame disturbance moment, thereby separating direct attitude robustness from the translational-gust scenarios.

The main contributions of this study are summarized as follows:

1.

A nonlinear six-degree-of-freedom Newton–Euler quadrotor model, based on an established comparative framework [6], is used within a common cascaded architecture for all controller evaluations.

2.

PID, LQR, $H_{\infty }$, and SMC are optimized and evaluated under identical plant dynamics, objective functions, actuator constraints, and simulation conditions to enable a consistent cross-controller comparison.

3.

RCGA and PSO are compared using identical objective functions, parameter bounds, plant models, and computational budgets. The resulting comparison is interpreted as controlled within-study evidence rather than a statistical claim of universal optimizer superiority.

4.

The optimized controllers are evaluated during sustained three-dimensional helical-trajectory tracking to assess coupled translational and rotational performance beyond conventional hover and isolated step responses.

5.

A constant-magnitude body-frame force acting through an aerodynamic center-of-pressure offset is swept over the complete ${360}^{\circ }$ azimuthal range at ${15}^{\circ }$ resolution to quantify directional sensitivity, stable coverage, failure sectors, recovery behavior, angular-rate response, and actuator saturation.

6.

Translational-gust and direct-moment scenarios are evaluated separately, allowing tracking degradation, attitude excursion, recovery, control effort, actuator saturation, and the mechanism responsible for any protocol-defined failure to be assessed without conflating the two disturbance mechanisms.

The remainder of this paper is organized as follows. Section 2 presents the modelling, controller design, optimization, stability, and evaluation methods. Section 3 reports the nominal, finite-disturbance, helical-tracking, and directional-robustness results. Section 4 discusses the principal findings, related work, practical implications, and limitations. Finally, Section 5 presents the conclusions and future research directions.

2. Materials and Methods

This section presents the common mathematical and numerical framework used for quadrotor modelling, controller design, optimization, stability analysis, simulation benchmarking, and directional-robustness evaluation.

2.1. Quadrotor System Modelling and Problem Formulation

A common nonlinear quadrotor model is used for controller design, optimization, and comparative evaluation. Analytical claims are restricted to the stated assumptions, while simulation results are interpreted as numerical evidence.

2.1.1. Reference Frames, Notation, and Modelling Assumptions

Two right-handed frames are employed: the North–East–Down (NED) inertial frame ${\mathcal{F}}_I=\{O_I;X_I,Y_I,Z_I\}$ and the body-fixed frame ${\mathcal{F}}_B=\{O_B;X_B,Y_B,Z_B\}$ attached to the centre of mass. The $X_B$, $Y_B$, and $Z_B$ axes point forward, starboard, and downward, respectively, consistent with standard Newton–Euler quadrotor modelling [4,5,7]. Under NED, altitude is $h=-Z$.

The quadrotor uses a plus (+) configuration with arm length ℓ, as illustrated in Figure 1. The virtual control inputs are

$$\begin{array}{cc}\hfill U_1& =T_1+T_2+T_3+T_4,\hfill \end{array}$$

(1a)

$$\begin{array}{cc}\hfill U_2& =T_2-T_4,\hfill \end{array}$$

(1b)

$$\begin{array}{cc}\hfill U_3& =T_1-T_3,\hfill \end{array}$$

(1c)

$$\begin{array}{cc}\hfill U_4& =k_T(T_1-T_2+T_3-T_4),\hfill \end{array}$$

(1d)

where $U_1$ is collective thrust, $U_2$ and $U_3$ are roll/pitch differential-thrust inputs, and $U_4$ is the yaw moment.

The model adopts the following assumptions:

(A1) 

The vehicle is a rigid body of constant mass m, with its centre of mass at $O_B$.

(A2) 

The airframe is symmetric, with

$$J=diag(I_x,I_y,I_z),I_x=I_y.$$

(2)

(A3) 

Rotor thrust and reaction torque satisfy

$$T_i=b{\Omega }_i^2,Q_i=d{\Omega }_i^2,$$

(3)

while blade flapping, ground effect, induced-flow transients, and rotor interactions are neglected.

(A4) 

Aerodynamic damping is approximated by

$$F_d=-K_t\dot{r},{\tau }_d=-K_r\omega ,$$

(4)

with

$$K_t=diag(k_x,k_y,k_z)>0,K_r=diag(k_p,k_q,k_r)>0,$$

(5)

consistent with comparative quadrotor models [6,26].

(A5) 

Attitude is represented by ZYX Euler angles

$$\eta ={[\varphi ,\theta ,\psi ]}^T,$$

(6)

with $\left|\theta \right|<\pi /2$.

(A6) 

A continuous-time full-state-feedback model is assumed; sensor noise, estimation errors, delays, discretization, and battery dynamics are neglected.

(A7) 

Environmental disturbances are bounded, piecewise-continuous forces applied at a fixed aerodynamic reference point, as defined in Section 2.2.

Principal symbols used throughout the formulation are summarized in Table 1.

2.1.2. Attitude Representation and Kinematics

The ZYX body-to-inertial rotation matrix is [7,10]

$$R_B^I=\left[\begin{array}{ccc}cos\psi cos\theta & cos\psi sin\theta sin\varphi -sin\psi cos\varphi & cos\psi sin\theta cos\varphi +sin\psi sin\varphi \\ sin\psi cos\theta & sin\psi sin\theta sin\varphi +cos\psi cos\varphi & sin\psi sin\theta cos\varphi -cos\psi sin\varphi \\ -sin\theta & cos\theta sin\varphi & cos\theta cos\varphi \end{array}\right].$$

(7)

The translational and attitude kinematics are

$$\dot{r}=v,v=R_B^I{v}^B,$$

(8)

and

$$\dot{\eta }=W\left(\eta \right)\omega ,$$

(9)

where

$$W\left(\eta \right)=\left[\begin{array}{ccc}1& sin\varphi tan\theta & cos\varphi tan\theta \\ 0& cos\varphi & -sin\varphi \\ 0& {\displaystyle \frac{sin\varphi }{cos\theta }}& {\displaystyle \frac{cos\varphi }{cos\theta }}\end{array}\right].$$

(10)

The transformation is singular at $\theta =\pm \pi /2$. At hover, $W\left(0\right)=I_3$, giving the local approximation $\dot{\eta }\approx \omega$ used for linear controller design; the nonlinear relation is retained in simulation.

2.1.3. Translational Dynamics

Newton’s second law in ${\mathcal{F}}_I$ gives

$$m\ddot{r}=mg{e}_3-U_1{R}_B^I{e}_3-K_t\dot{r}+F_w^I,$$

(11)

with

$$F_w^I=R_B^I{F}_w^B.$$

(12)

The component dynamics are

$$\begin{array}{cc}\hfill \ddot{X}& =-{\displaystyle \frac{U_1}{m}}(cos\psi sin\theta cos\varphi +sin\psi sin\varphi )-{\displaystyle \frac{k_x}{m}}\dot{X}+{\displaystyle \frac{F_{w,x}^I}{m}},\hfill \end{array}$$

(13a)

$$\begin{array}{cc}\hfill \ddot{Y}& =-{\displaystyle \frac{U_1}{m}}(sin\psi sin\theta cos\varphi -cos\psi sin\varphi )-{\displaystyle \frac{k_y}{m}}\dot{Y}+{\displaystyle \frac{F_{w,y}^I}{m}},\hfill \end{array}$$

(13b)

$$\begin{array}{cc}\hfill \ddot{Z}& =g-{\displaystyle \frac{U_1}{m}}cos\theta cos\varphi -{\displaystyle \frac{k_z}{m}}\dot{Z}+{\displaystyle \frac{F_{w,z}^I}{m}}.\hfill \end{array}$$

(13c)

These equations constitute the nonlinear translational model [4,5,8]. Horizontal motion is generated indirectly by tilting the collective-thrust vector, motivating the cascaded position/attitude architecture used in this study.

2.1.4. Rotational Dynamics

Including rotor gyroscopic effects, aerodynamic damping, control moments, and external disturbances, Euler’s rigid-body equation is

$$J\dot{\omega }+\omega \times (J\omega +h_r)={\tau }_c-K_r\omega +{\tau }_w^B,$$

(14)

where

$$h_r=J_r{\Omega }_r{e}_3,$$

(15)

$${\Omega }_r={\Omega }_1-{\Omega }_2+{\Omega }_3-{\Omega }_4,$$

(16)

and

$${\tau }_c={[\ell U_2,\ell U_3,U_4]}^T.$$

(17)

The component dynamics are

$$\begin{array}{cc}\hfill \dot{p}& ={\displaystyle \frac{I_y-I_z}{I_x}}qr-{\displaystyle \frac{J_r}{I_x}}q{\Omega }_r-{\displaystyle \frac{k_p}{I_x}}p+{\displaystyle \frac{\ell }{I_x}}{U}_2+{\displaystyle \frac{{\tau }_{w,x}^B}{I_x}},\hfill \end{array}$$

(18a)

$$\begin{array}{cc}\hfill \dot{q}& ={\displaystyle \frac{I_z-I_x}{I_y}}pr+{\displaystyle \frac{J_r}{I_y}}p{\Omega }_r-{\displaystyle \frac{k_q}{I_y}}q+{\displaystyle \frac{\ell }{I_y}}{U}_3+{\displaystyle \frac{{\tau }_{w,y}^B}{I_y}},\hfill \end{array}$$

(18b)

$$\begin{array}{cc}\hfill \dot{r}& ={\displaystyle \frac{I_x-I_y}{I_z}}pq-{\displaystyle \frac{k_r}{I_z}}r+{\displaystyle \frac{U_4}{I_z}}+{\displaystyle \frac{{\tau }_{w,z}^B}{I_z}}.\hfill \end{array}$$

(18c)

These standard nonlinear rotational equations [4,5,8,9] retain inertial coupling and rotor gyroscopic effects.

2.1.5. Rotor Aerodynamics and Control Allocation

The mapping between squared rotor speeds and virtual inputs is

$$\left[\begin{array}{c}{U}_1\\ U_2\\ U_3\\ U_4\end{array}\right]=\left[\begin{array}{cccc}b& b& b& b\\ 0& b& 0& -b\\ b& 0& -b& 0\\ d& -d& d& -d\end{array}\right]\left[\begin{array}{c}{\Omega }_1^2\\ {\Omega }_2^2\\ {\Omega }_3^2\\ {\Omega }_4^2\end{array}\right].$$

(19)

Equivalently,

$$\begin{array}{cc}\hfill U_1& =b({\Omega }_1^2+{\Omega }_2^2+{\Omega }_3^2+{\Omega }_4^2),\hfill \end{array}$$

(20a)

$$\begin{array}{cc}\hfill U_2& =b({\Omega }_2^2-{\Omega }_4^2),\hfill \end{array}$$

(20b)

$$\begin{array}{cc}\hfill U_3& =b({\Omega }_1^2-{\Omega }_3^2),\hfill \end{array}$$

(20c)

$$\begin{array}{cc}\hfill U_4& =d({\Omega }_1^2-{\Omega }_2^2+{\Omega }_3^2-{\Omega }_4^2).\hfill \end{array}$$

(20d)

Defining the allocation matrix as $\Gamma$,

$$\left[\begin{array}{c}{\Omega }_1^2\\ {\Omega }_2^2\\ {\Omega }_3^2\\ {\Omega }_4^2\end{array}\right]={\Gamma }^{-1}U,$$

(21)

with physically admissible rotor speeds

$${\Omega }_i=\sqrt{max\left(0,{\left[{\Gamma }^{-1}U\right]}_i\right)},i=1,\dots ,4.$$

(22)

The reaction-torque-to-thrust ratio is

$$\kappa ={\displaystyle \frac{d}{b}}.$$

(23)

For the adopted parameters, $\kappa =0.0444$ m, indicating lower direct yaw authority than the roll/pitch authority generated through the arm length $\ell =0.2$ m.

2.1.6. Compact Nonlinear State-Space Representation

Define

$$x={[r^T,v^T,{\eta }^T,{\omega }^T]}^T\in {\mathbb{R}}^{12},u={[U_1,U_2,U_3,U_4]}^T\in {\mathbb{R}}^4,$$

(24)

and

$$w\left(t\right)={[{\left(F_w^B\right)}^T,{\left({\tau }_w^B\right)}^T]}^T\in {\mathbb{R}}^6.$$

(25)

The constrained nonlinear plant is

$$\dot{x}=f\left(x\right)+G\left(x\right)sat\left(u\right)+D\left(x\right)w\left(t\right),y=x,$$

(26)

where

$$f\left(x\right)=\left[\begin{array}{c}v\\ g{e}_3-{\displaystyle \frac{1}{m}}{K}_{t}v\\ W\left(\eta \right)\omega \\ -J^{-1}\{\omega \times [J\omega +h_r]+K_r\omega \}\end{array}\right],$$

(27)

$$G\left(x\right)=\left[\begin{array}{cc}{0}_{3\times 1}& 0_{3\times 3}\\ -{\displaystyle \frac{1}{m}}{R}_B^I{e}_3& 0_{3\times 3}\\ 0_{3\times 1}& 0_{3\times 3}\\ 0_{3\times 1}& J^{-1}diag(\ell ,\ell ,1)\end{array}\right],$$

(28)

and

$$D\left(x\right)=\left[\begin{array}{cc}{0}_{3\times 3}& 0_{3\times 3}\\ {\displaystyle \frac{1}{m}}{R}_B^I& 0_{3\times 3}\\ 0_{3\times 3}& 0_{3\times 3}\\ 0_{3\times 3}& J^{-1}\end{array}\right].$$

(29)

The model therefore captures nonlinear coupling, underactuation, disturbance inputs, and actuator constraints within a common plant.

2.1.7. Model and Actuator Parameters

The parameters in Table 2, adopted from the comparative model of Rinaldi et al. [6], are used for all controllers. The Lumenier QAV250 [32] is cited only as a representative small-quadrotor size class; the simulation parameters are not identified from that platform.

The hover rotor speed is

$${\Omega }_h=\sqrt{{\displaystyle \frac{mg}{4b}}}=639\mathrm{rad}·{\mathrm{s}}^{-1}\approx 6100\mathrm{rpm}.$$

(30)

The corresponding moment limits are

$${\overline{\tau }}_{\varphi }={\overline{\tau }}_{\theta }=\ell {\overline{U}}_{23}=1.766\mathrm{N}·\mathrm{m},$$

(31)

and

$${\overline{\tau }}_{\psi }={\overline{U}}_4=1.177\mathrm{N}·\mathrm{m}.$$

(32)

Thus, yaw authority is approximately $33.3\%$ lower than the maximum roll/pitch moment authority.

2.1.8. Actuator Saturation and Command-Shaping Constraints

Actuator dynamics are represented by instantaneous bounded virtual inputs, consistent with [6]. The applied control vector is

$$sat\left(u\right)={\left[{sat}_{0,{\overline{U}}_1}\left(U_1\right),{sat}_{-{\overline{U}}_{23},{\overline{U}}_{23}}\left(U_2\right),{sat}_{-{\overline{U}}_{23},{\overline{U}}_{23}}\left(U_3\right),{sat}_{-{\overline{U}}_4,{\overline{U}}_4}\left(U_4\right)\right]}^T,$$

(33)

where

$${sat}_{a,b}\left(s\right)=min\{b,max\{a,s\}\}.$$

(34)

The common actuator bounds are

$$\begin{array}{cc}\hfill 0\le U_1& \le 3mg=17.658\mathrm{N},\hfill \end{array}$$

(35a)

$$\begin{array}{cc}\hfill -1.5mg\le U_2,U_3& \le 1.5mg=8.829\mathrm{N},\hfill \end{array}$$

(35b)

$$\begin{array}{cc}\hfill -mg\ell \le U_4& \le mg\ell =1.177\mathrm{N}·\mathrm{m}.\hfill \end{array}$$

(35c)

Outer-loop commands are additionally constrained by

$$a_c\leftarrow {sat}_{-\overline{a},\overline{a}}\left(a_c\right),\overline{a}={[10,10,7]}^T\mathrm{m}·{\mathrm{s}}^{-2},$$

(36)

and

$$|{\varphi }_{\mathrm{ref}}|\le {25}^{\circ },\left|{\theta }_{\mathrm{ref}}\right|\le {25}^{\circ }.$$

(37)

For a steady tilt $\theta$, vertical force balance requires

$$U_1\ge {\displaystyle \frac{mg}{cos\theta }}.$$

(38)

With ${\overline{U}}_1=3mg$, the corresponding theoretical thrust-limited tilt is

$${\theta }_{max}={cos}^{-1}\left({\displaystyle \frac{mg}{{\overline{U}}_1}}\right)={cos}^{-1}\left({\displaystyle \frac{1}{3}}\right)\approx 70.5^{\circ }.$$

(39)

Similarly, steady disturbance rejection requires

$$|{\tau }_{w,x}^B|\le {\overline{\tau }}_{\varphi },|{\tau }_{w,y}^B|\le {\overline{\tau }}_{\theta },\left|{\tau }_{w,z}^B\right|\le {\overline{\tau }}_{\psi }.$$

(40)

These actuator constraints are applied identically to all controllers and are explicitly considered in the subsequent directional-robustness assessment.

2.2. Environmental Disturbance Model

The wind disturbance is represented as a body-frame force applied at an offset aerodynamic centre of pressure. An off-centre force generates both a translational force and a disturbance moment [7]:

$$F_w^I\left(t\right)=R_B^I\left(t\right)F_w^B\left(t\right),{\tau }_w^B\left(t\right)=r_{\mathrm{cp}}\times F_w^B\left(t\right).$$

(41)

Equivalently,

$${\tau }_w^B\left(t\right)=S\left(r_{\mathrm{cp}}\right)F_w^B\left(t\right),$$

(42)

where

$$S\left(r_{\mathrm{cp}}\right)=\left[\begin{array}{ccc}0& -l_z& l_y\\ l_z& 0& -l_x\\ -l_y& l_x& 0\end{array}\right],r_{\mathrm{cp}}={[l_x,l_y,l_z]}^T.$$

(43)

In the position-tracking scenarios, $F_w^I$ acts through the translational dynamics, whereas in the attitude-only directional study, the induced moment ${\tau }_w^B$ acts directly on the rotational dynamics. This distinction preserves the physical mechanism of each benchmark.

2.2.1. Fixed-Direction Disturbance-Window Study

A prescribed body-frame force is applied over a finite interval using the smooth activation function

$$g\left(t\right)={\displaystyle \frac{1}{2}}\left[1+tanh\left({\displaystyle \frac{t-t_{\mathrm{on}}}{{\tau }_g}}\right)\right]{\displaystyle \frac{1}{2}}\left[1+tanh\left({\displaystyle \frac{t_{\mathrm{off}}-t}{{\tau }_g}}\right)\right],$$

(44)

with

$$F_w^B\left(t\right)=g\left(t\right){\overline{F}}_w^B,$$

(45)

where ${\tau }_g=0.25$ s. Identical disturbance magnitude, direction, timing, references, initial conditions, and actuator limits are used for all controllers in the step-tracking and helical-trajectory scenarios.

2.2.2. Directional Attitude-Robustness Study

The directional study varies only the horizontal azimuth of a constant-magnitude body-frame disturbance:

$${\mathbf{F}}_w^B\left(\alpha \right)=\left[\begin{array}{c}{F}_{h}cos\alpha \\ F_{h}sin\alpha \\ F_z\end{array}\right],{∥{\mathbf{F}}_w^B\left(\alpha \right)∥}_2=\sqrt{F_h^2+F_z^2}.$$

(46)

For D3, $F_h=8$ N and $F_z=-1$ N, giving $\parallel {\mathbf{F}}_w^B{\parallel }_2=\sqrt{65}=8.06$ N. The azimuth convention is

$$\begin{array}{cccc}\hfill \alpha & =0^{\circ }:\hfill & \hfill & +x_B,\hfill \\ \hfill \alpha & ={90}^{\circ }:\hfill & \hfill & +y_B,\hfill \\ \hfill \alpha & ={180}^{\circ }:\hfill & \hfill & -x_B,\hfill \\ \hfill \alpha & ={270}^{\circ }:\hfill & \hfill & -y_B.\hfill \end{array}$$

(47)

The sweep is sampled every ${15}^{\circ }$:

$${\alpha }_k={15}^{\circ }k,k=0,1,\dots ,23.$$

(48)

The repeated ${360}^{\circ }$ point is used only to close directional plots.

The corresponding disturbance moment is

$${\tau }_w^B\left({\alpha }_k\right)={\mathbf{r}}_{cp}^B\times {\mathbf{F}}_w^B\left({\alpha }_k\right)=S\left({\mathbf{r}}_{cp}^B\right){\mathbf{F}}_w^B\left({\alpha }_k\right).$$

(49)

For each azimuth, attitude excursion, geometric tilt, angular rate, control effort, actuator saturation, and recovery time are recorded. Geometric tilt is

$${\Theta }_g\left(t\right)={cos}^{-1}\left[{\mathbf{e}}_3^T{\mathbf{R}}_{BI}\left(t\right){\mathbf{e}}_3\right].$$

(50)

A direction is classified as stable when all prescribed excursion, rate, saturation, and recovery criteria are satisfied. The stable-direction set is

$${\mathcal{A}}_c=\{{\alpha }_k\in [0^{\circ },{360}^{\circ }]:\mathrm{controller}c\mathrm{satisfies}\mathrm{all}\mathrm{criteria}\},$$

(51)

with coverage

$${\mathcal{C}}_c={\displaystyle \frac{\mu \left({\mathcal{A}}_c\right)}{{360}^{\circ }}}\times 100\%.$$

(52)

Figure 2 summarizes the directional disturbance framework used in this study.

2.3. Statement of the Control Problem

Consider the constrained nonlinear system of Equations (26)–(29), with the parameters of Table 2, the disturbance model above, and the common actuator constraints.

Let

$$\mathcal{K}=\{\mathrm{PID},\mathrm{LQR},{\mathcal{H}}_{\infty },\mathrm{SMC}\},\mathcal{M}=\{\mathrm{RCGA},\mathrm{PSO}\}.$$

For controller $K_i\in \mathcal{K}$ and optimizer $M_j\in \mathcal{M}$, the tuned parameter vector is

$$p_{ij}^{*}=arg\underset{p\in {\mathcal{P}}_i}{min}{J}_i\left(p\right),$$

(53)

where $J_i\left(p\right)$ is defined in Section 2.5.2.

Each optimized configuration is then assessed on the same plant, references, disturbance models, and constraints in terms of nominal tracking, three-dimensional trajectory tracking, disturbance recovery, and directional stability coverage. The closed loop must satisfy

$$x\left(t\right)\in {\mathcal{X}}_{\mathrm{adm}},sat\left(u\left(t\right)\right)\in {\mathcal{U}}_{\mathrm{adm}},t\in [0,T].$$

(54)

For the directional study, the principal quantity of interest is the stable-direction set ${\mathcal{A}}_c$ and its boundary:

$$\partial {\mathcal{A}}_c=\mathrm{boundary}\mathrm{of}\{\alpha \in [0^{\circ },{360}^{\circ }]:\mathrm{stability}\mathrm{and}\mathrm{recovery}\mathrm{criteria}\mathrm{are}\mathrm{satisfied}\}.$$

(55)

2.4. Controller Synthesis and Metaheuristic Tuning Framework

2.4.1. Cascaded Architecture and Outer Position Loop

All configurations use the cascaded architecture in Figure 3. The outer loop regulates position and generates acceleration, attitude, and collective-thrust commands; the inner loop is implemented using PID, LQR, ${\mathcal{H}}_{\infty }$, or SMC. The outer-loop structure is common to all configurations [1,33].

The position error is

$$e_p=r_d-r,$$

(56)

and the filtered outer-loop PID is

$$\begin{array}{cc}\hfill {\dot{I}}_p& =e_p,\hfill \end{array}$$

(57a)

$$\begin{array}{cc}\hfill {\dot{v}}_f& =N_f(-v-v_f),\hfill \end{array}$$

(57b)

$$\begin{array}{cc}\hfill a_c& =K_p^o{e}_p+K_i^o{I}_p+K_d^o{v}_f.\hfill \end{array}$$

(57c)

Using measured velocity avoids derivative kick for piecewise-constant references [2].

The acceleration command is limited by

$${\overline{a}}_c={sat}_{-\overline{a},\overline{a}}\left(a_c\right),\overline{a}={[10,10,7]}^T\mathrm{m}·{\mathrm{s}}^{-2},$$

(58)

and mapped to yaw-compensated attitude references:

$$\begin{array}{cc}\hfill {\varphi }_{\mathrm{ref}}& ={\displaystyle \frac{{\overline{a}}_{c,y}cos{\psi }_d-{\overline{a}}_{c,x}sin{\psi }_d}{g}},\hfill \end{array}$$

(59a)

$$\begin{array}{cc}\hfill {\theta }_{\mathrm{ref}}& =-{\displaystyle \frac{{\overline{a}}_{c,x}cos{\psi }_d+{\overline{a}}_{c,y}sin{\psi }_d}{g}}.\hfill \end{array}$$

(59b)

The collective thrust is

$$U_1=m(g-{\overline{a}}_{c,z}).$$

(60)

For ${\psi }_d=0$, this reduces to the mapping used by Rinaldi et al. [6]:

$${\varphi }_{\mathrm{ref}}={\displaystyle \frac{{\overline{a}}_{c,y}}{g}},{\theta }_{\mathrm{ref}}=-{\displaystyle \frac{{\overline{a}}_{c,x}}{g}}.$$

(61)

The common command limits are

$$|{\varphi }_{\mathrm{ref}}|\le {25}^{\circ },\left|{\theta }_{\mathrm{ref}}\right|\le {25}^{\circ },0\le U_1\le {\overline{U}}_1.$$

(62)

2.4.2. Inner-Loop Attitude Control Laws

Define

$${\eta }_{\mathrm{ref}}={[{\varphi }_{\mathrm{ref}},{\theta }_{\mathrm{ref}},{\psi }_{\mathrm{ref}}]}^T,e_{\eta }=\eta -{\eta }_{\mathrm{ref}},$$

(63)

with wrapped yaw error

$$e_{\psi }={wrap}_{(-\pi ,\pi ]}(\psi -{\psi }_{\mathrm{ref}}).$$

(64)

All controllers receive the same references, states, constraints, and disturbances.

PID Control  

The filtered PID controller is

$$\begin{array}{cc}\hfill {\dot{I}}_{\eta }& =e_{\eta },\hfill \end{array}$$

(65a)

$$\begin{array}{cc}\hfill {\dot{v}}_{\eta }& =N_f^a(\dot{\eta }-v_{\eta }),\hfill \end{array}$$

(65b)

with

$$u_{\eta }=-K_p^a{e}_{\eta }-K_i^a{I}_{\eta }-K_d^a{v}_{\eta },$$

(66)

where

$$u_{\eta }={[U_2,U_3,U_4]}^T.$$

(67)

The corresponding filtered PID form is

$$C_{\mathrm{PID}}\left(s\right)=K_p+{\displaystyle \frac{K_i}{s}}+{\displaystyle \frac{K_d{N}_f^{a}s}{s+N_f^a}}.$$

(68)

A common gain set is used across the three attitude channels [2,34]. No explicit anti-windup is included; therefore, prolonged saturation may cause integral accumulation and delayed recovery [2,34,35].

Linear Quadratic Regulator  

The LQR controller is designed from the hover-linearized attitude model

$${\dot{x}}_A=A{x}_A+B{u}_{\eta },$$

(69)

where

$$A=\left[\begin{array}{cc}{0}_{3\times 3}& I_3\\ 0_{3\times 3}& -J^{-1}{K}_r\end{array}\right],B=\left[\begin{array}{c}{0}_{3\times 3}\\ J^{-1}diag(\ell ,\ell ,1)\end{array}\right].$$

(70)

The quadratic objective is

$$J_{\mathrm{LQR}}={\int }_0^{\infty }\left(x_A^{T}Q{x}_A+u_{\eta }^{T}R{u}_{\eta }\right)dt,$$

(71)

with

$$Q=diag(q_{\varphi },q_{\theta },q_{\psi },q_p,q_q,q_r),R=diag(r_2,r_3,r_4).$$

(72)

The Riccati equation

$$A^{T}P+PA-PB{R}^{-1}{B}^{T}P+Q=0$$

(73)

gives

$$K_{\mathrm{LQR}}=R^{-1}{B}^{T}P,u_{\eta }=-K_{\mathrm{LQR}}{x}_A.$$

(74)

Under the standard stabilizability and detectability conditions, the linearized closed loop is asymptotically stable [36,37,38]. Metaheuristic tuning is therefore performed over the positive diagonal weights in Q and R. Since no integral action is included, constant disturbance moments may produce a nonzero steady-state error [36,37,39].

Robust ${\mathcal{H}}_{\infty }$ Control  

The ${\mathcal{H}}_{\infty }$ controller uses the same linearized attitude model and seeks to bound the worst-case gain from exogenous input w to weighted performance output z:

$$\begin{array}{cc}\hfill {\dot{x}}_A& =A{x}_A+B_{1}w+B_2{u}_{\eta },\hfill \end{array}$$

(75a)

$$\begin{array}{cc}\hfill z& =C_1{x}_A+D_{11}w+D_{12}{u}_{\eta },\hfill \end{array}$$

(75b)

$$\begin{array}{cc}\hfill y& =C_2{x}_A+D_{21}w+D_{22}{u}_{\eta }.\hfill \end{array}$$

(75c)

The generalized-plant matrices are

$$B_1=I_6,B_2=B,C_2=I_6,D_{21}=\epsilon I_6,D_{22}=0,\epsilon ={10}^{-6},$$

(76)

and

$$C_1=\left[\begin{array}{cc}{W}_a& 0\\ 0& W_r\\ 0& 0\end{array}\right],D_{11}=0,D_{12}=\left[\begin{array}{c}0\\ 0\\ W_u\end{array}\right].$$

(77)

Thus,

$$z={[{\left(W_a{e}_{\eta }\right)}^T,{\left(W_r\omega \right)}^T,{\left(W_u{u}_{\eta }\right)}^T]}^T.$$

(78)

The controller is synthesized such that

$$\parallel T_{zw}{\parallel }_{\infty }<\gamma ,$$

(79)

equivalently,

$${\parallel z\parallel }_2<\gamma {\parallel w\parallel }_2.$$

(80)

The corresponding DGKF Riccati conditions follow standard ${\mathcal{H}}_{\infty }$ theory [40], with LMI formulations available when required [41]. The resulting guarantee applies to the linearized generalized plant and excludes actuator saturation and unmodelled structured uncertainty [40,42].

Sliding-Mode Control  

For SMC, define

$$e_s={\eta }_{\mathrm{ref}}-\eta ,$$

(81)

and the sliding surface

$$s={\dot{e}}_s+B_{\beta }{e}_s,B_{\beta }=diag({\beta }_1,{\beta }_2,{\beta }_3).$$

(82)

On $s=0$,

$${\dot{e}}_s=-B_{\beta }{e}_s,$$

(83)

giving exponential error decay [43,44].

Using $\dot{\eta }\approx \omega$ near hover,

$$s\approx B_{\beta }{e}_s-\omega .$$

(84)

The rotational dynamics are written as

$$\dot{\omega }=f_{\omega }\left(\omega \right)+M_u{u}_{\eta }+J^{-1}{\tau }_w^B,$$

(85)

with

$$M_u=diag\left({\displaystyle \frac{\ell }{I_x}},{\displaystyle \frac{\ell }{I_y}},{\displaystyle \frac{1}{I_z}}\right).$$

(86)

The equivalent control is

$$u_{\eta ,\mathrm{eq}}=M_u^{-1}\left[-f_{\omega }\left(\omega \right)+B_{\beta }{\dot{e}}_s\right],$$

(87)

and the implemented SMC law is

$${\mathbf{u}}_{\eta }={\mathbf{u}}_{\eta ,\mathrm{eq}}+{\mathbf{M}}_u^{-1}{\mathbf{K}}_s{sat}_{\epsilon }\left(\mathbf{s}\right),{\mathbf{K}}_s\succ 0.$$

(88)

The continuous boundary-layer approximation is

$${sat}_{{\epsilon }_i}\left(s_i\right)={\displaystyle \frac{s_i}{\sqrt{s_i^2+{\epsilon }_i}}},{\epsilon }_i>0,$$

(89)

which reduces chattering at the expense of exact sliding [43,45]. The resulting dynamics are

$$\dot{\mathbf{s}}=-{\mathbf{K}}_s{sat}_{\epsilon }\left(\mathbf{s}\right)+{\mathbf{d}}_m,$$

(90)

where

$$d_m=-J^{-1}{\tau }_w^B+\delta (x,t).$$

(91)

Hence, direct disturbance moments are matched with respect to the attitude loop, whereas translational wind forces affect attitude indirectly through the cascaded outer loop [43,44,46,47].

2.5. Metaheuristic Tuning Framework

2.5.1. Decision Variables and Admissible Sets

Each inner-loop controller is tuned jointly with the common outer position controller, so the comparison concerns complete closed-loop configurations. Let $p_c\in {\mathcal{P}}_c$ denote the decision vector and admissible parameter box for controller configuration c. Shared gains are used across the three position channels and, where applicable, across the roll, pitch, and yaw channels. The adopted variables and bounds are listed in Table 3.

Positive LQR and ${\mathcal{H}}_{\infty }$ weights preserve their respective synthesis conditions. For SMC,

$$T_s={\displaystyle \frac{1}{\beta }}\in [0.0167,0.200]\mathrm{s},$$

(92)

and $\epsilon >0$ maintains the smooth boundary-layer approximation. Because large switching gains may exceed the actuator authority, all candidates are evaluated using the saturated nonlinear model.

2.5.2. Unified Objective Function

The same objective function, weights, thresholds, simulation horizon, and metric definitions are used for all controller families:

$$J\left(p\right)=\sum _{i\in \{X,Y,Z\}}{w}_i^P\left[{\mathrm{ITAE}}_i+{\lambda }_A{\mathrm{IAE}}_i\right]+\sum _{j\in \{\varphi ,\theta ,\psi \}}{w}_j^A\left[{\mathrm{ITAE}}_j+{\lambda }_A{\mathrm{IAE}}_j\right]+R\left(p\right).$$

(93)

For tracking error $e_k\left(t\right)$,

$${\mathrm{IAE}}_k={\int }_0^{T_f}|e_k\left(t\right)|dt,{\mathrm{ITAE}}_k={\int }_0^{T_f}t\left|e_k\left(t\right)\right|dt.$$

(94)

The penalty term is

$$R\left(p\right)={\lambda }_{\mathrm{os}}\sum _k{(M_{p,k}-{\overline{M}}_p)}_{+}+{\lambda }_{\mathrm{ts}}\sum _k{(t_{s,k}-{\overline{t}}_s)}_{+}+{\lambda }_{\omega }{\parallel \omega \parallel }_{\mathrm{RMS}},$$

(95)

where

$${\left(x\right)}_{+}=max(0,x),$$

(96)

and

$${\parallel \omega \parallel }_{\mathrm{RMS}}=\sqrt{{\displaystyle \frac{1}{T_f}}{\int }_0^{T_f}(p^2+q^2+r^2)dt}.$$

(97)

For non-settling responses, $t_s=2T_f$ is used as a finite penalty surrogate.

The common weights are

$$\begin{array}{cc}\hfill w^P& =[1,1,1.2],w^A=[0.6,0.6,1.0],\hfill \\ \hfill {\lambda }_A& =0.5,{\lambda }_{\mathrm{os}}=0.8,{\lambda }_{\mathrm{ts}}=0.6,{\lambda }_{\omega }=0.05,T_f=20\mathrm{s}.\hfill \end{array}$$

(98)

The optimal parameter vector is

$$p^{*}(c,o)=arg\underset{p\in {\mathcal{P}}_c}{min}{J}_c\left(p\right),c\in \{\mathrm{PID},\mathrm{LQR},{\mathcal{H}}_{\infty },\mathrm{SMC}\},o\in \{\mathrm{RCGA},\mathrm{PSO}\}.$$

(99)

Final assessment also uses physical metrics such as RMSE, IAE, geometric tilt, recovery time, saturation, and control effort.

2.5.3. Fairness of the Optimization Comparison

Within each controller family, RCGA and PSO use identical parameter bounds, objective functions, plant models, initial conditions, references, disturbance models, solvers, and random initialization seeds. Both use a population or swarm size of 30.

For PID, LQR, and SMC,

$$N_{\mathrm{eval}}=30\times 20=600,$$

(100)

whereas for the more computationally expensive ${\mathcal{H}}_{\infty }$ synthesis,

$$N_{\mathrm{eval},{\mathcal{H}}_{\infty }}=30\times 10=300.$$

(101)

Thus, RCGA and PSO are budget-matched within each controller family.

Each optimizer is executed once per configuration using a reproducible seed. The comparison therefore represents controlled single-run evidence rather than a statistical claim of optimizer superiority; multiple independent runs and nonparametric testing would be required for the latter [48]. Cross-controller conclusions are based jointly on J and the reported physical performance metrics.

2.5.4. Real-Coded Genetic Algorithm

RCGA operates directly on the continuous parameter vectors of Table 3 [49,50]. The population at generation g is

$${\mathcal{P}}^{\left(g\right)}=\{p_1^{\left(g\right)},p_2^{\left(g\right)},\dots ,p_N^{\left(g\right)}\},p_i^{\left(g\right)}\in {\mathcal{P}}_c,$$

(102)

with uniform initialization

$$p_i^{\left(0\right)}\sim \mathcal{U}[p^{-},p^{+}].$$

(103)

Fitness is evaluated through the nonlinear closed-loop simulation:

$$J_i^{\left(g\right)}=J\left(p_i^{\left(g\right)}\right).$$

(104)

Tournament selection and elitism are employed [51], with

$$J_{\mathrm{best}}^{(g+1)}\le J_{\mathrm{best}}^{\left(g\right)}.$$

(105)

BLX-$\alpha$ crossover generates offspring according to

$$p_j^c\sim \mathcal{U}[\underset{j}{min}-\alpha {\Delta }_j,\underset{j}{max}+\alpha {\Delta }_j],$$

(106)

where

$$\underset{j}{min}=min(p_j^{\left(1\right)},p_j^{\left(2\right)}),\underset{j}{max}=max(p_j^{\left(1\right)},p_j^{\left(2\right)}),{\Delta }_j=|p_j^{\left(1\right)}-p_j^{\left(2\right)}|,$$

(107)

with $\alpha =0.5$ [52,53].

Gaussian mutation is

$$p_j^m=p_j^c+{\sigma }_j{\xi }_j,{\xi }_j\sim \mathcal{N}(0,1),$$

(108)

where

$${\sigma }_j=s_m(p_j^{+}-p_j^{-}).$$

(109)

Candidates are projected onto the admissible box:

$$p\leftarrow {\Pi }_{{\mathcal{P}}_c}\left(p\right),{\left[{\Pi }_{{\mathcal{P}}_c}\left(p\right)\right]}_j=min\{p_j^{+},max\{p_j^{-},p_j\}\}.$$

(110)

The convergence criterion is

$${\displaystyle \frac{|J_{\mathrm{best}}^{\left(g\right)}-J_{\mathrm{best}}^{(g-1)}|}{max\{1,|J_{\mathrm{best}}^{(g-1)}\left|\right\}}}<{10}^{-3}.$$

(111)

The RCGA formulation follows [49,50,53,54].

2.5.5. Particle Swarm Optimization

PSO updates each particle using personal and global best information [55]. The velocity update is

$$v_i^{(k+1)}=w{v}_i^{\left(k\right)}+c_1{r}_1\odot (b_i-p_i^{\left(k\right)})+c_2{r}_2\odot (b^g-p_i^{\left(k\right)}),$$

(112)

followed by

$$p_i^{(k+1)}={\Pi }_{{\mathcal{P}}_c}\left[p_i^{\left(k\right)}+v_i^{(k+1)}\right].$$

(113)

When a component reaches a bound,

$$v_{i,j}^{(k+1)}=0\mathrm{if}{p}_{i,j}^{(k+1)}\in \{p_j^{-},p_j^{+}\}.$$

(114)

The personal and global best solutions are updated as

$$b_i\leftarrow p_i^{(k+1)}\mathrm{if}J\left(p_i^{(k+1)}\right)<J\left(b_i\right),$$

(115)

and

$$b^g=arg\underset{b_i}{min}J\left(b_i\right).$$

(116)

The adopted parameters are

$$w=0.7,c_1=c_2=1.5,$$

(117)

consistent with established PSO convergence guidance [56,57,58].

PSO uses the same objective, parameter bounds, population size, model, and stopping conditions as RCGA within each controller family. The two methods therefore differ primarily in their search mechanisms, and neither is assumed to dominate a priori [49].

2.6. Stability Analysis

This section summarizes the analytical stability properties of the closed-loop framework under the stated assumptions. Analytical results are distinguished from numerical synthesis conditions and simulation-based evidence.

2.6.1. Equilibrium and Hover Trim

For the undisturbed nonlinear model with inactive saturation, an equilibrium $(x^{*},u^{*})$ satisfies

$$f\left(x^{*}\right)+G\left(x^{*}\right)u^{*}=0.$$

(118)

At hover,

$${\omega }^{*}=0,U_2^{*}=U_3^{*}=U_4^{*}=0,$$

(119)

and translational force balance requires

$$mg{e}_3=U_1^{*}{R}_B^I\left({\eta }^{*}\right)e_3.$$

(120)

Hence, $U_1^{*}=mg$ and ${\varphi }^{*}={\theta }^{*}=0$, while position and yaw remain arbitrary. The hover-equilibrium set is therefore

$$\mathcal{E}=\left\{(x^{*},u^{*}):x^{*}={[{\left(r^{*}\right)}^T,0_3^T,0,0,{\psi }^{*},0_3^T]}^T,u^{*}={[mg,0,0,0]}^T\right\}.$$

(121)

The hover thrust occupies one third of the available collective-thrust range:

$${\displaystyle \frac{U_1^{*}}{{\overline{U}}_1}}={\displaystyle \frac{1}{3}}.$$

(122)

2.6.2. Local Linearization

Linearization about hover with ${\psi }^{*}=0$ gives

$$\begin{array}{cc}\hfill \delta \ddot{X}& =-g\delta \theta -{\displaystyle \frac{k_x}{m}}\delta \dot{X},\hfill \end{array}$$

(123a)

$$\begin{array}{cc}\hfill \delta \ddot{Y}& =g\delta \varphi -{\displaystyle \frac{k_y}{m}}\delta \dot{Y},\hfill \end{array}$$

(123b)

$$\begin{array}{cc}\hfill \delta \ddot{Z}& =-{\displaystyle \frac{1}{m}}\delta U_1-{\displaystyle \frac{k_z}{m}}\delta \dot{Z}.\hfill \end{array}$$

(123c)

Thus, vertical motion is directly controlled by $U_1$, whereas horizontal motion is generated through roll and pitch.

The attitude subsystem is

$$\delta {\dot{x}}_A=A\delta x_A+B\delta u_{\eta },$$

(124)

with A and B defined in Equation (70). The imposed ${25}^{\circ }$ attitude-command limit supports the local linear approximation during nominal operation, but does not provide a formal guarantee under large disturbances. Large-angle behaviour is therefore assessed using the nonlinear saturated model [39].

2.6.3. Controllability, Stabilizability, and Detectability

For the attitude subsystem, define

$$G_{\tau }=J^{-1}diag(\ell ,\ell ,1).$$

(125)

The controllability matrix contains

$$\left[\begin{array}{cc}B& AB\end{array}\right]=\left[\begin{array}{cc}{0}_{3\times 3}& G_{\tau }\\ G_{\tau }& -J^{-1}{K}_r{G}_{\tau }\end{array}\right].$$

(126)

Since $G_{\tau }$ is nonsingular,

$$rank\left[\begin{array}{cc}B& AB\end{array}\right]=6,$$

(127)

so $(A,B)$ is controllable and therefore stabilizable.

With $C_2=I_6$, $(C_2,A)$ is observable and detectable. These conditions ensure the existence of the stabilizing LQR solution for admissible positive weights [36,37].

The full hover-linearized model similarly satisfies

$$rank\left({\mathcal{C}}_{12}\right)=12\mathrm{for}{U}_1^{*}=mg>0.$$

(128)

For ${\mathcal{H}}_{\infty }$ synthesis,

$$D_{12}^T{D}_{12}=W_u^T{W}_u\succ 0,D_{21}{D}_{21}^T={\epsilon }^2{I}_6\succ 0,$$

(129)

and the required stabilizability, detectability, and regularity conditions are satisfied over the adopted positive weighting ranges [40,59]. These conditions apply to the linearized generalized plant and do not extend to saturation or large-angle motion.

2.6.4. Controller-Specific Stability Analysis

PID Control  

For attitude channel i,

$$G_i\left(s\right)={\displaystyle \frac{g_i}{s(s+d_i)}},d_i={\displaystyle \frac{k_r}{I_i}},$$

(130)

with

$$g_{\varphi }={\displaystyle \frac{\ell }{I_x}},g_{\theta }={\displaystyle \frac{\ell }{I_y}},g_{\psi }={\displaystyle \frac{1}{I_z}}.$$

(131)

Combining the plant with the filtered PID law gives

$${\Delta }_i\left(s\right)=s^4+a_3{s}^3+a_2{s}^2+a_{1}s+a_0,$$

(132)

where

$$\begin{array}{cc}\hfill a_3& =d_i+N_f^a,\hfill \end{array}$$

(133a)

$$\begin{array}{cc}\hfill a_2& =d_i{N}_f^a+g_i(K_p+K_d{N}_f^a),\hfill \end{array}$$

(133b)

$$\begin{array}{cc}\hfill a_1& =g_i(K_p{N}_f^a+K_i),\hfill \end{array}$$

(133c)

$$\begin{array}{cc}\hfill a_0& =g_i{K}_i{N}_f^a.\hfill \end{array}$$

(133d)

The quartic Routh–Hurwitz conditions are

$$a_3{a}_2>a_1,a_3{a}_2{a}_1>a_1^2+a_3^2{a}_0.$$

(134)

These conditions are checked for all three attitude channels. When satisfied, the linearized PID system is locally exponentially stable; the result excludes saturation and large-angle motion [2,34,35].

Linear Quadratic Regulator  

The LQR closed-loop matrix is

$$A_{\mathrm{cl}}=A-B{K}_{\mathrm{LQR}}.$$

(135)

For $V=x_A^{T}P{x}_A$,

$$\dot{V}=-x_A^T\left(Q+K_{\mathrm{LQR}}^{T}R{K}_{\mathrm{LQR}}\right)x_A<0.$$

(136)

Hence, the linearized attitude-error system is exponentially stable and the corresponding nonlinear hover equilibrium is locally exponentially stable by Lyapunov’s indirect method [38]. This does not establish global stability under nonlinearities or input constraints [60].

Robust ${\mathcal{H}}_{\infty }$ Control  

The synthesized controller internally stabilizes the generalized linear plant and satisfies

$$\parallel T_{zw}{\parallel }_{\infty }<\gamma .$$

(137)

Equivalently,

$${\parallel z\parallel }_2<\gamma {\parallel w\parallel }_2.$$

(138)

The guarantee applies to the linearized unsaturated generalized plant and represents weighted disturbance attenuation rather than robustness to arbitrary parametric uncertainty, saturation, or large-angle dynamics [40,41,42,59].

Sliding-Mode Control  

For sliding channel i,

$${\dot{s}}_i=-K_i{sat}_a\left(s_i\right)+d_i,\left|d_i\right|\le D_i.$$

(139)

With

$$V_i={\displaystyle \frac{1}{2}}{s}_i^2,$$

(140)

ideal switching gives

$${\dot{V}}_i\le -(K_i-D_i)\left|s_i\right|.$$

(141)

Thus, $K_i>D_i$ guarantees finite-time reaching, with

$$t_{r,i}\le {\displaystyle \frac{|s_i\left(0\right)|}{K_i-D_i}}.$$

(142)

For the implemented boundary layer,

$${\dot{V}}_i\le -{\displaystyle \frac{K_i{s}_i^2}{\sqrt{s_i^2+a_i}}}+D_i\left|s_i\right|,$$

(143)

yielding

$$|s_i|\le {\delta }_i,{\delta }_i={\displaystyle \frac{\sqrt{a_i}{D}_i}{\sqrt{K_i^2-D_i^2}}},$$

(144)

and

$$\underset{t\to \infty }{lim\; sup}|s_i\left(t\right)|\le {\delta }_i,\underset{t\to \infty }{lim\; sup}\left|e_i\left(t\right)\right|\le {\displaystyle \frac{{\delta }_i}{{\beta }_i}}.$$

(145)

These guarantees require sufficient actuator authority to realize the switching action [43,44,45,46].

2.6.5. Bounded Disturbances, Saturation, and Model Uncertainty

For locally exponentially stable unsaturated closed-loop dynamics, sufficiently small bounded disturbances yield a local input-to-state stability bound

$$\parallel x\left(t\right)\parallel \le \beta (\parallel x\left(0\right)\parallel ,t)+\chi \left(\underset{0\le \tau \le t}{sup}\parallel w\left(\tau \right)\parallel \right),$$

(146)

where $\beta$ and $\chi$ are class-$\mathcal{KL}$ and class-$\mathcal{K}$ functions, respectively [39].

When actuators saturate, the applied control differs from the unsaturated law and the preceding PID, LQR, and ${\mathcal{H}}_{\infty }$ guarantees no longer apply directly. Likewise, SMC reaching may fail if the required moment cannot be delivered.

A persistent disturbance moment cannot be statically balanced if

$$|{\tau }_{w,i}^B|>{\overline{\tau }}_i,i\in \{\varphi ,\theta ,\psi \}.$$

(147)

In that case, no constant-attitude equilibrium exists for that channel [61,62,63].

No controller is synthesized over an explicit parametric-uncertainty set. Accordingly, the analytical results support local hover stability and bounded local disturbance responses, while the large-disturbance and saturation behaviour is assessed numerically.

2.6.6. Directional Control Authority: Analytical Robustness Boundary

The directional robustness boundary is partly determined by actuator authority and therefore applies independently of the selected controller.

Proposition 1  

Consider a constant body-frame disturbance moment ${\tau }_w^B$ subject to

$$|U_2|,|U_3|\le {\overline{U}}_{23},\left|U_4\right|\le {\overline{U}}_4.$$

(148)

A necessary condition for a constant-attitude equilibrium with $\omega =0$ is

$$|{\tau }_{w,\varphi }^B|\le {\overline{\tau }}_{\varphi },|{\tau }_{w,\theta }^B|\le {\overline{\tau }}_{\theta },\left|{\tau }_{w,\psi }^B\right|\le {\overline{\tau }}_{\psi },$$

(149)

where

$${\overline{\tau }}_{\varphi }={\overline{\tau }}_{\theta }=\ell {\overline{U}}_{23}=1.766\mathrm{N}\mathrm{m},{\overline{\tau }}_{\psi }={\overline{U}}_4=1.177\mathrm{N}\mathrm{m}.$$

(150)

At equilibrium, the control moment must balance the disturbance moment. Violation of any inequality in Equation (149) therefore eliminates exact static trim, although finite-duration recovery may still be possible.

For D3,

$${\mathbf{r}}_{cp}^B={\left[\begin{array}{ccc}0.08& 0.08& 0.05\end{array}\right]}^T\mathrm{m},{\mathbf{F}}_w^B\left(\alpha \right)={\left[\begin{array}{ccc}8cos\alpha & 8sin\alpha & -1\end{array}\right]}^T\mathrm{N},$$

(151)

which gives

$${\tau }_w^B\left(\alpha \right)=\left[\begin{array}{c}-0.08-0.40sin\alpha \\ 0.08+0.40cos\alpha \\ 0.64(sin\alpha -cos\alpha )\end{array}\right]\mathrm{N}\mathrm{m}.$$

(152)

The maximum roll and pitch disturbance moments are

$$\underset{\alpha }{max}|{\tau }_{w,\varphi }^B|=\underset{\alpha }{max}\left|{\tau }_{w,\theta }^B\right|=0.48\mathrm{N}\mathrm{m}<1.766\mathrm{N}\mathrm{m}.$$

(153)

For yaw,

$${\tau }_{w,\psi }^B\left(\alpha \right)=0.64\sqrt{2}sin(\alpha -{45}^{\circ }),$$

(154)

so

$$\underset{\alpha }{max}\left|{\tau }_{w,\psi }^B\right|=0.905\mathrm{N}\mathrm{m}<1.177\mathrm{N}\mathrm{m}.$$

(155)

Hence, the necessary static-trim conditions are satisfied for all disturbance azimuths,

$$\alpha \in [0^{\circ },{360}^{\circ }),$$

(156)

including all sampled directions in ${\mathcal{D}}_{360}=\{0^{\circ },{15}^{\circ },\dots ,{345}^{\circ }\}$.

2.7. Simulation and Benchmarking Formulation

2.7.1. Numerical Setup

Three disturbance scenarios are considered: D1 evaluates disturbed step tracking, D2 evaluates disturbed three-dimensional helical tracking, and D3 evaluates directional attitude robustness. All controllers are simulated using the continuous-time nonlinear model and MATLAB’s adaptive Dormand–Prince solver (ode45). The numerical settings are summarized in Table 4.

D3 uses event-based termination at the prescribed near-flip or angular-rate threshold. All benchmarks are continuous-time simulations; discrete-time implementation effects are not modelled.

2.7.2. Reference Signals

Scenario S1: Nominal Step Response  

A simultaneous position and yaw step is applied at $t_0=0.1$ s:

$$r_d\left(t\right)={[2,1,1]}^T\mathbb{1}(t\ge t_0)\mathrm{m},$$

(157)

$${\psi }_d\left(t\right)={15}^{\circ }\mathbb{1}(t\ge t_0),{\varphi }_d={\theta }_d=0.$$

(158)

This scenario evaluates nominal transient and steady-state tracking performance.

Scenario S2: Three-Dimensional Helical Trajectory  

The reference trajectory is

$$r_d\left(t\right)={[Rcos\left({\omega }_{h}t\right)-R,Rsin\left({\omega }_{h}t\right),z_0+v_{z}t]}^T,$$

(159)

with

$$R=3\mathrm{m},{\omega }_h=0.6\mathrm{rad}·{\mathrm{s}}^{-1},v_z=0.25\mathrm{m}·{\mathrm{s}}^{-1},z_0=0.3\mathrm{m}.$$

(160)

The trajectory provides sustained coupled three-dimensional motion with a centripetal acceleration of

$$a_{\mathrm{cp}}=R{\omega }_h^2=1.08\mathrm{m}·{\mathrm{s}}^{-2}.$$

(161)

2.7.3. Disturbance Scenarios

D1 and D2 apply body-frame force disturbances to the full 12-state model, whereas D3 isolates attitude robustness under a direct aerodynamic disturbance moment.

Table 5. Disturbance scenarios.

Scenario	${\mathit{F}}_{\mathit{w}}^{\mathit{B}}$ [N]	Magnitude	Window [s]	Applied moment	Plant
D1: disturbed step	${[1.5,-1.0,0.6]}^T$	$1.90$ N	$[15,18]$	No; $r_{\mathrm{cp}}=0$	12-state
D2: disturbed helix	${[2.5,-1.0,-2.6]}^T$	$3.74$ N	$[18,21]$	No; $r_{\mathrm{cp}}=0$	12-state
D3: directional sweep	${[8cos\alpha ,8sin\alpha ,-1]}^T$	$8.06$ N	$[15,18]$	${\tau }_w^B=r_{\mathrm{cp}}^B\times F_w^B$	Attitude-only

Table 5. Disturbance scenarios.

Scenario	${\mathit{F}}_{\mathit{w}}^{\mathit{B}}$ [N]	Magnitude	Window [s]	Applied moment	Plant
D1: disturbed step	${[1.5,-1.0,0.6]}^T$	$1.90$ N	$[15,18]$	No; $r_{\mathrm{cp}}=0$	12-state
D2: disturbed helix	${[2.5,-1.0,-2.6]}^T$	$3.74$ N	$[18,21]$	No; $r_{\mathrm{cp}}=0$	12-state
D3: directional sweep	${[8cos\alpha ,8sin\alpha ,-1]}^T$	$8.06$ N	$[15,18]$	${\tau }_w^B=r_{\mathrm{cp}}^B\times F_w^B$	Attitude-only

For D1 and D2,

$$F_w^I=R_B^I{F}_w^B,$$

(162)

with $r_{\mathrm{cp}}=0$; hence, the disturbance acts directly on translation and affects attitude indirectly through the cascaded controller.

For D3,

$$F_w^B\left(\alpha \right)={[8cos\alpha ,8sin\alpha ,-1]}^T\mathrm{N}.$$

(163)

Using $F_D=\frac{1}{2}\rho C_D{A}_{\mathrm{ref}}{V}_{\mathrm{eq}}^2$ [64], with $\rho =1.225$ kg m⁻³, $A_{\mathrm{ref}}=0.040$ m², and $C_D=1$, the $8.06$ N resultant corresponds to

$$V_{\mathrm{eq}}=\sqrt{{\displaystyle \frac{2\parallel F_w^B{\parallel }_2}{\rho C_D{A}_{\mathrm{ref}}}}}=18.1\mathrm{m}{\mathrm{s}}^{-1}.$$

(164)

This value is an equivalent severity indicator rather than an identified airframe wind speed [65].

The force acts at

$$r_{\mathrm{cp}}^B={[0.08,0.08,0.05]}^T\mathrm{m},$$

(165)

producing

$${\tau }_w^B\left(\alpha \right)={\mathbf{r}}_{\mathrm{cp}}^B\times {\mathbf{F}}_w^B\left(\alpha \right).$$

(166)

D3 uses the attitude-only model with $U_1=mg$ and no propagated translational states.

2.7.4. Performance Metrics

For tracking error $e_k\left(t\right)=r_k\left(t\right)-y_k\left(t\right)$, the common integral metrics are

$${\mathrm{IAE}}_k={\int }_0^{T_f}\left|e_k\left(t\right)\right|dt,$$

(167)

$${\mathrm{ITAE}}_k={\int }_0^{T_f}t\left|e_k\left(t\right)\right|dt,$$

(168)

and

$${\mathrm{RMSE}}_k=\sqrt{{\displaystyle \frac{1}{T_f}}{\int }_0^{T_f}{e}_k^2\left(t\right)dt}.$$

(169)

For step responses, percentage overshoot is

$$M_{p,k}=100{\displaystyle \frac{max\{0,y_{\mathrm{peak},k}-y_{d,k}\}}{|y_{d,k}-y_k\left(0\right)|}},$$

(170)

and the 2% settling time is

$$t_{s,k}=inf\left\{t:|e_k\left(\tau \right)|\le 0.02|\Delta y_k|,\forall \tau \in [t,T_f]\right\}.$$

(171)

Non-settling responses are assigned $t_{s,k}=2T_f$.

The total control-effort index is

$$E_u={\int }_0^{T_f}{\displaystyle \sum _{i=1}^4}{U}_i^2\left(t\right)dt.$$

(172)

Because the virtual inputs have different physical units, $E_u$ is interpreted only as a relative control-activity index.

For D3, $U_1=mg$ is controller-independent; therefore, attitude-control activity is evaluated using

$$E_{\mathrm{diff}}\left({\alpha }_k\right)={\int }_0^{T_f}\left[U_2^2+U_3^2+U_4^2\right]dt,$$

(173)

with directional mean

$${\overline{E}}_{\mathrm{diff},c}={\displaystyle \frac{1}{M}}\sum _{{\alpha }_k\in {\mathcal{D}}_{360}}{E}_{\mathrm{diff},c}\left({\alpha }_k\right).$$

(174)

Actuator saturation is quantified by

$${\rho }_{\mathrm{sat}}={\displaystyle \frac{100}{T_f}}{\int }_0^{T_f}\mathbb{1}\left[\underset{i}{max}\left({\displaystyle \frac{|U_i\left(t\right)|}{{\overline{U}}_i}}\right)\ge 0.98\right]dt.$$

(175)

Attitude-excursion, recovery-time, and directional-robustness metrics are defined in Section 2.8.2.

2.8. Roll/Pitch Failure-Angle Sweep Protocol

The directional sweep evaluates how horizontal wind direction affects roll–pitch stability. For each azimuth angle, the body-frame force and centre-of-pressure moment are constructed, and the attitude-only closed-loop model is simulated. Attitude excursion, angular-rate magnitude, recovery time, control effort, and actuator saturation are then evaluated using common criteria. The overall directional roll/pitch failure-angle sweep protocol is illustrated in Figure 4. The protocol produces sampled stable and failed directional sectors rather than a single failure angle.

2.8.1. Direction Set and Disturbance Construction

The complete horizontal-azimuth grid is

$${\mathcal{D}}_{360}=\left\{0^{\circ },{15}^{\circ },{30}^{\circ },\dots ,{345}^{\circ }\right\},M=24,\Delta \alpha ={15}^{\circ }.$$

(176)

For polar visualization only, the response at $0^{\circ }$ is repeated at ${360}^{\circ }$ to close each envelope. This repeated endpoint is excluded from ${\mathcal{D}}_{360}$, the stable-direction set ${\mathcal{S}}_c$, and the coverage measure ${\Gamma }_c$.

Because the disturbance includes a fixed vertical component, the azimuth convention refers to its horizontal projection:

$$\begin{array}{cccccccc}\hfill \alpha =0^{\circ }& :\hfill & \hfill & +x_B,\hfill & \hfill \alpha ={90}^{\circ }& :\hfill & \hfill & +y_B,\hfill \\ \hfill \alpha ={180}^{\circ }& :\hfill & \hfill & -x_B,\hfill & \hfill \alpha ={270}^{\circ }& :\hfill & \hfill & -y_B.\hfill \end{array}$$

(177)

For each ${\alpha }_k\in {\mathcal{D}}_{360}$, the body-frame disturbance force is

$${\mathbf{F}}_w^B\left({\alpha }_k\right)=\left[\begin{array}{c}8cos{\alpha }_k\\ 8sin{\alpha }_k\\ -1\end{array}\right]\mathrm{N}.$$

(178)

Its magnitude is independent of direction:

$${∥{\mathbf{F}}_w^B\left({\alpha }_k\right)∥}_2=\sqrt{8^2+1^2}=\sqrt{65}=8.06\mathrm{N}.$$

(179)

The force acts at the body-frame centre-of-pressure offset

$${\mathbf{r}}_{\mathrm{cp}}^B={\left[\begin{array}{ccc}0.08& 0.08& 0.05\end{array}\right]}^T\mathrm{m},$$

(180)

and generates the disturbance moment

$${\tau }_w^B\left({\alpha }_k\right)={\mathbf{r}}_{\mathrm{cp}}^B\times {\mathbf{F}}_w^B\left({\alpha }_k\right).$$

(181)

Expanding Equation (181) gives

$${\tau }_w^B\left({\alpha }_k\right)=\left[\begin{array}{c}-0.08-0.40sin{\alpha }_k\\ \phantom{-}0.08+0.40cos{\alpha }_k\\ \phantom{-}0.64(sin{\alpha }_k-cos{\alpha }_k)\end{array}\right]\mathrm{N}\mathrm{m}.$$

(182)

The disturbance is applied smoothly over $15\le t\le 18$ s using the activation function defined in Equation (44).

The physical plant contains the six attitude states

$${\mathbf{x}}_a={\left[\begin{array}{cccccc}\varphi & \theta & \psi & p& q& r\end{array}\right]}^T,$$

(183)

with the necessary controller states appended. The collective thrust is fixed at $U_1=mg$, and the translational states are not propagated. The D3 sweep therefore evaluates directional attitude robustness under directly applied disturbance moments.

2.8.2. Attitude-Excursion, Angular-Rate, and Coverage Measures

For each direction ${\alpha }_k$, the maximum absolute roll and pitch excursions are

$${\varphi }_{max}\left({\alpha }_k\right)=\underset{t\in [0,T_f]}{max}|\varphi (t;{\alpha }_k)|,{\theta }_{max}\left({\alpha }_k\right)=\underset{t\in [0,T_f]}{max}\left|\theta (t;{\alpha }_k)\right|.$$

(184)

The maximum individual-axis roll–pitch excursion used by the stability classifier is

$${\Theta }_{\mathrm{exc},\max }\left({\alpha }_k\right)=max\left\{{\varphi }_{max}\left({\alpha }_k\right),{\theta }_{max}\left({\alpha }_k\right)\right\}.$$

(185)

Because the classification evaluates roll–pitch integrity, the roll–pitch angular-rate magnitude is defined as

$${\omega }_{rp}(t;{\alpha }_k)=\sqrt{p^2(t;{\alpha }_k)+q^2(t;{\alpha }_k)},{\omega }_{rp,max}\left({\alpha }_k\right)=\underset{t\in [0,T_f]}{max}{\omega }_{rp}(t;{\alpha }_k).$$

(186)

The maximum component rates are retained separately:

$$\begin{array}{cc}\hfill p_{max}\left({\alpha }_k\right)& =\underset{t\in [0,T_f]}{max}\left|p(t;{\alpha }_k)\right|,\hfill \\ \hfill q_{max}\left({\alpha }_k\right)& =\underset{t\in [0,T_f]}{max}\left|q(t;{\alpha }_k)\right|,\hfill \\ \hfill r_{max}\left({\alpha }_k\right)& =\underset{t\in [0,T_f]}{max}\left|r(t;{\alpha }_k)\right|.\hfill \end{array}$$

(187)

For diagnostic purposes, the full body-rate norm is also calculated:

$${∥\omega (t;{\alpha }_k)∥}_2=\sqrt{p^2(t;{\alpha }_k)+q^2(t;{\alpha }_k)+r^2(t;{\alpha }_k)},$$

(188)

$${∥\omega ∥}_{max}\left({\alpha }_k\right)=\underset{t\in [0,T_f]}{max}{∥\omega (t;{\alpha }_k)∥}_2.$$

(189)

The full body-rate norm is reported only as a diagnostic quantity because it can be dominated by yaw. Accordingly, ${\omega }_{rp,max}$ is used for roll–pitch classification, while $r_{max}$ is reported separately to identify yaw-dominated motion.

For complementary diagnostic reporting, the coordinate-independent geometric tilt is calculated as:

$${\Theta }_g(t;{\alpha }_k)={cos}^{-1}\left[cos\varphi (t;{\alpha }_k)cos\theta (t;{\alpha }_k)\right],$$

(190)

$${\Theta }_{g,max}\left({\alpha }_k\right)=\underset{t\in [0,T_f]}{max}{\Theta }_g(t;{\alpha }_k).$$

(191)

The geometric tilt ${\Theta }_{g,max}$ is reported as a coordinate-independent diagnostic measure and is not used in the stability classifier. Classification is based on the maximum individual-axis roll–pitch excursion ${\Theta }_{\mathrm{exc},\max }$ defined in Equation (185).

The recovery time is the earliest time after disturbance removal for which the roll–pitch excursion and roll–pitch rate remain within their recovery bounds:

$$t_{\mathrm{rec}}\left({\alpha }_k\right)=inf\left\{t-t_{\mathrm{off}}:\begin{array}{c}max\left\{\right|\varphi (s;{\alpha }_k)|,|\theta (s;{\alpha }_k)\left|\right\}\le 5^{\circ },\hfill \\ {\omega }_{rp}(s;{\alpha }_k)\le 1\mathrm{rad}{\mathrm{s}}^{-1},\hfill \\ \forall s\in [t,T_f]\hfill \end{array}\right\}.$$

(192)

For controller c, the sampled stable-direction set is

$${\mathcal{S}}_c=\left\{{\alpha }_k\in {\mathcal{D}}_{360}:\begin{array}{c}{\Theta }_{\mathrm{exc},max}\left({\alpha }_k\right)<{45}^{\circ },\hfill \\ {\omega }_{\mathrm{RP},max}\left({\alpha }_k\right)<30\mathrm{rad}{\mathrm{s}}^{-1},\hfill \\ t_{\mathrm{rec}}\left({\alpha }_k\right)\le 8\mathrm{s},\hfill \\ {\rho }_{\mathrm{sat}}\left({\alpha }_k\right)\le 35\%\hfill \end{array}\right\}.$$

(193)

The stable angular coverage and its percentage are

$${\Gamma }_c=\Delta \alpha \left|{\mathcal{S}}_c\right|,C_c={\displaystyle \frac{100{\Gamma }_c}{{360}^{\circ }}}.$$

(194)

Because $\Delta \alpha ={15}^{\circ }$, the identified sector boundaries have a directional resolution of ${15}^{\circ }$ and should be interpreted as sampled rather than exact continuous critical angles.

2.8.3. Failure Classification

The common classification thresholds are

$${\Theta }_{\mathrm{safe}}={45}^{\circ },{\Theta }_{\mathrm{flip}}={80}^{\circ },{\omega }_{rp,\lim }=30\mathrm{rad}{\mathrm{s}}^{-1},{\overline{\rho }}_{\mathrm{sat}}=35\%.$$

(195)

The ${45}^{\circ }$ safe-tilt boundary was selected to represent a practical multirotor operational limit rather than a universal instability boundary. This value is consistent with the default maximum in-air tilt specified by the PX4 multicopter position controller and the default maximum lean angle documented for ArduPilot Copter [66,67]. The remaining thresholds are explicit protocol choices introduced in this study to ensure identical classification across all controllers. Specifically, ${80}^{\circ }$ denotes a near-inverted roll–pitch condition, $30\mathrm{rad}{\mathrm{s}}^{-1}$ bounds excessive roll–pitch rate, and a saturation duty above $35\%$ identifies prolonged operation at the actuator limits. Similarly, the $8\mathrm{s}$ recovery window and the persistent recovery bands of $5^{\circ }$ and $1\mathrm{rad}{\mathrm{s}}^{-1}$ are study-defined evaluation criteria. They should therefore not be interpreted as universal physical limits of quadrotor flight.

The actuator-saturation duty ${\rho }_{\mathrm{sat}}$ is calculated using Equation (175). Each tested direction is assigned one mutually exclusive classification according to

$$S\left({\alpha }_k\right)=\left\{\begin{array}{cc}\mathrm{invalid},\hfill & \mathrm{if}\mathrm{nonfinite}\mathrm{states}\mathrm{or}\mathrm{control}\mathrm{inputs}\mathrm{occur},\hfill \\ \mathrm{flipped}\mathrm{or}\mathrm{near}-\mathrm{flipped},\hfill & \mathrm{if}{\Theta }_{\mathrm{exc},\max }\left({\alpha }_k\right)\ge {80}^{\circ },\hfill \\ \mathrm{roll}--\mathrm{pitch}\mathrm{excursion}\mathrm{exceeded},\hfill & \mathrm{if}{\Theta }_{\mathrm{exc},\max }\left({\alpha }_k\right)\ge {45}^{\circ },\hfill \\ \mathrm{roll}--\mathrm{pitch}\mathrm{rate}\mathrm{exceeded},\hfill & \mathrm{if}{\omega }_{\mathrm{RP},\max }\left({\alpha }_k\right)\ge 30\mathrm{rad}{\mathrm{s}}^{-1},\hfill \\ \mathrm{failed}\mathrm{to}\mathrm{recover},\hfill & \mathrm{if}{t}_{\mathrm{rec}}\left({\alpha }_k\right)>8\mathrm{s},\hfill \\ \mathrm{saturation}\mathrm{exceeded},\hfill & \mathrm{if}{\rho }_{\mathrm{sat}}\left({\alpha }_k\right)>35\%,\hfill \\ \mathrm{stable},\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(196)

The event-based solver terminates if the componentwise roll–pitch excursion ${\Theta }_{\mathrm{RP}}(t;{\alpha }_k)$ reaches ${80}^{\circ }$ or the roll–pitch angular-rate magnitude ${\omega }_{\mathrm{RP}}(t;{\alpha }_k)$ reaches $60\mathrm{rad}{\mathrm{s}}^{-1}$. The latter is a numerical safety threshold equal to twice the $30\mathrm{rad}{\mathrm{s}}^{-1}$ classification limit and is not used as the reported stability boundary.

Yaw error, $r_{max}$, and the full body-rate norm are reported separately but are not used as roll–pitch failure criteria. This distinction prevents yaw-authority-driven motion from being misclassified as roll–pitch instability.

2.8.4. Recovery Criterion

Disturbance resistance and post-disturbance recovery are evaluated separately. Let $t_{\mathrm{off}}=18\mathrm{s}$ and define the instantaneous componentwise roll–pitch excursion as

$${\Theta }_{\mathrm{RP}}(t;{\alpha }_k)=max\left\{|\varphi (t;{\alpha }_k)|,|\theta (t;{\alpha }_k)|\right\}.$$

(197)

Accordingly, ${\Theta }_{\mathrm{exc},\max }\left({\alpha }_k\right)={max}_{t\in [0,T_f]}{\Theta }_{\mathrm{RP}}(t;{\alpha }_k)$. The recovery time is

$$t_{\mathrm{rec}}\left({\alpha }_k\right)=inf\left\{t-t_{\mathrm{off}}:\begin{array}{cc}\hfill {\Theta }_{\mathrm{RP}}(\tau ;{\alpha }_k)& \le 5^{\circ },\hfill \\ \hfill {\omega }_{\mathrm{RP}}(\tau ;{\alpha }_k)& \le 1\mathrm{rad}{\mathrm{s}}^{-1},\hfill \\ \hfill & \forall \tau \in [t,T_f]\hfill \end{array}\right\},t\in [t_{\mathrm{off}},t_{\mathrm{off}}+8\mathrm{s}].$$

(198)

The persistence requirement prevents a temporary threshold crossing from being interpreted as recovery. If no admissible recovery time exists within the $8\mathrm{s}$ recovery window, $t_{\mathrm{rec}}=\infty$. Simulations terminated by near-flip or numerical rate-safety events are classified before recovery is evaluated. The yaw rate r and the full-body-rate norm ${\parallel \omega \parallel }_2$ are reported diagnostically but are not included in the roll–pitch recovery criterion.

2.8.5. Directional Stability Sectors and Coverage

For controller c, the stable, failed, and recovered direction sets are

$$\begin{array}{cc}\hfill {\mathcal{S}}_c& =\{{\alpha }_k\in {\mathcal{D}}_{360}:S\left({\alpha }_k\right)=\mathrm{stable}\},\hfill \end{array}$$

(199a)

$$\begin{array}{cc}\hfill {\mathcal{F}}_c& ={\mathcal{D}}_{360}\setminus {\mathcal{S}}_c,\hfill \end{array}$$

(199b)

$$\begin{array}{cc}\hfill {\mathcal{R}}_c& =\{{\alpha }_k\in {\mathcal{D}}_{360}:t_{\mathrm{rec}}\left({\alpha }_k\right)\le 8\mathrm{s}\}.\hfill \end{array}$$

(199c)

The recovered-direction set ${\mathcal{R}}_c$ is evaluated independently of the complete stable-direction set ${\mathcal{S}}_c$. Consequently, a direction may satisfy the recovery criterion while failing the complete stability protocol because of excessive actuator-saturation duty or another classification condition. Adjacent samples with the same classification are grouped into maximal circular sectors. A sector crossing ${360}^{\circ }/0^{\circ }$ is reported as a wrap-around sector.

The stable angular coverage is

$${\Gamma }_c=\Delta \alpha \left|{\mathcal{S}}_c\right|,C_c={\displaystyle \frac{100{\Gamma }_c}{{360}^{\circ }}},\Delta \alpha ={15}^{\circ }.$$

(200)

The same calculation is applied to the recovered-direction set. Because the sweep is discrete, the reported transition directions have a resolution of ${15}^{\circ }$ and represent boundaries between sampled sectors rather than exact continuous critical angles.

2.8.6. Scope of the Protocol

The protocol quantifies a sampled directional attitude-robustness envelope for the prescribed disturbance magnitude, centre-of-pressure offset, duration, actuator limits, and classification thresholds. It measures a different property from nominal tracking accuracy; therefore, both tracking and directional robustness metrics are required for controller comparison.

The results are conditional on the attitude-only D3 model with $U_1=mg$ and no propagated translational states. They characterize roll/pitch integrity, flip-like behaviour, and recovery under direct disturbance moments. Yaw error is excluded from failure classification but reported separately.

Accordingly, the sweep provides simulation evidence for the discrete set ${\mathcal{D}}_{360}$ with ${15}^{\circ }$ resolution. It does not constitute a continuous robustness margin, global region-of-attraction proof, or universal wind-withstand limit.

3. Results

This section reports the controller-optimization outcomes and the closed-loop results obtained under nominal step tracking, finite-duration wind disturbance, helical-trajectory tracking, and the directional roll–pitch robustness sweep.

3.1. Optimization Results and Controller Selection

The PID, LQR, ${\mathrm{H}}_{\infty }$, and SMC architectures were optimized using GA and PSO under the same objective function and simulation framework. For each architecture, the parameter set with the lowest objective-function value was retained for the subsequent evaluations. The resulting objective values and selected configurations are listed in Table 6, while the complete baseline, GA-optimized, and PSO-optimized parameter sets required for reproduction are reported in Appendix A.1.

PSO yielded the lowest objective value for all four architectures. The selected values were 11.70, 9.54, 5.23, and 4.28 for PID, LQR, ${\mathrm{H}}_{\infty }$, and SMC, respectively.

3.2. Nominal Step-Tracking Results

The selected PSO configurations and the baseline PID–PID controller were evaluated for $x_{\mathrm{ref}}=2\mathrm{m}$, $y_{\mathrm{ref}}=1\mathrm{m}$, $z_{\mathrm{ref}}=1\mathrm{m}$, ${\varphi }_{\mathrm{ref}}={\theta }_{\mathrm{ref}}=0^{\circ }$, and ${\psi }_{\mathrm{ref}}={15}^{\circ }$. Figure 5 presents the position and attitude responses, and Table 7 lists the aggregated response metrics.

The baseline position response did not satisfy the adopted settling criterion within the simulation interval. Among the optimized configurations, PID–SMC recorded the shortest aggregated position rise and settling times and the lowest position IAE and ITAE. PID–PID recorded the lowest position overshoot. For the aggregated attitude response, PID–SMC recorded the shortest rise and settling times and the lowest overshoot, while PID–${\mathrm{H}}_{\infty }$ recorded the lowest ITAE.

3.3. Step Tracking Under a Finite-Duration Wind Disturbance

In D1, the nominal step references were retained and the body-frame wind force ${\mathbf{F}}_w^B={[1.5,-1.0,0.6]}^T\mathrm{N}$, with magnitude $1.90\mathrm{N}$, was applied during $15\le t\le 18\mathrm{s}$. Figure 6 presents the resulting position and attitude responses.

Each row in Table 8 reports one controller under the same D1 conditions. The position metrics aggregate the x, y, and z responses, whereas the attitude metrics aggregate the $\varphi$, $\theta$, and $\psi$ responses. The reported $T_r$, $T_s$, and $M_p$ are the worst-case values across the corresponding three axes; IAE and ITAE are summed across those axes. Thus, the attitude entries are not yaw-specific quantities. The corresponding individual axis-wise position and attitude metrics are reported in Appendix A.2.

All controller configurations remained closed-loop stable during D1 and returned toward the commanded position after the disturbance interval. PID–SMC recorded the lowest position IAE and ITAE and the shortest aggregated position settling time. PID–${\mathrm{H}}_{\infty }$ and PID–SMC both recorded the lowest position overshoot. PID–LQR recorded the shortest aggregated attitude settling time, while PID–SMC recorded the shortest attitude rise time and the lowest attitude overshoot.

3.4. Helical-Trajectory Tracking Under Wind Disturbance

In D2, the four PSO-optimized controller pairs tracked the helical reference defined under Scenario S2 in Section 2.7.2. The body-frame wind force ${\mathbf{F}}_w^B={[2.5,-1.0,-2.6]}^T\mathrm{N}$, with magnitude $3.74\mathrm{N}$, was applied during $18\le t\le 21\mathrm{s}$.

Representative frames from the three-dimensional simulations are shown in Figure 7. Figure 8 compares the complete trajectories, and Figure 9 presents the corresponding position and attitude time responses. All four controllers completed the helical trajectory and remained closed-loop stable throughout the simulation.

For D2, let $\gamma \left(t\right)={\parallel {\mathbf{r}}_d\left(t\right)-\mathbf{r}\left(t\right)\parallel }_2$ denote the three-dimensional position-error magnitude. The D2 recovery threshold is calculated over the two-second pre-disturbance interval as ${\gamma }_{\mathrm{rec}}=mean\left(\gamma \right)+2std\left(\gamma \right)$. The D2 recovery time $t_{\mathrm{rec},\mathrm{D}2}$ is the earliest time after disturbance removal at which $\gamma \left(t\right)\le {\gamma }_{\mathrm{rec}}$ continuously for $1\mathrm{s}$. This moving-reference recovery metric differs from the attitude-based D3 recovery criterion. The normalized control effort is defined as

$$J_u={\int }_0^{T_f}{\displaystyle \sum _{i=1}^4}{\left({\displaystyle \frac{U_i\left(t\right)}{U_{i,lim}}}\right)}^{2}dt.$$

(201)

The D2 tracking, recovery, tilt, and normalized-control-effort metrics are listed in Table 9.

SMC recorded the lowest RMSE, MAE, IAE, ISE, wind-window RMSE, and maximum wind-window position error. LQR recorded the shortest recovery time and the lowest normalized control effort. The maximum geometric tilt ranged from $30.{417}^{\circ }$ to $30.{616}^{\circ }$, and, for every controller, the maximum value occurred within the D2 wind window.

Figure 10 reports the LQR control response, for which $J_u=3.6910\mathrm{s}$. The corresponding normalized effort values were $3.7113\mathrm{s}$ for PID, $3.7648\mathrm{s}$ for ${\mathrm{H}}_{\infty }$, and $3.7942\mathrm{s}$ for SMC.

3.5. Directional Roll–Pitch Robustness Results

The D3 disturbance direction was swept over $0^{\circ }\le {\alpha }_w<{360}^{\circ }$ using identical initial conditions, actuator constraints, disturbance duration, and classification criteria for all controllers. Table 10 summarizes the directional classification, attitude, angular-rate, effort, saturation, and recovery results.

LQR, ${\mathrm{H}}_{\infty }$, and SMC satisfied the complete stability protocol over the full directional sweep. PID satisfied the protocol over 58.33% of the sweep; its failed classifications occurred for ${105}^{\circ }\le {\alpha }_w<{180}^{\circ }$ and ${285}^{\circ }\le {\alpha }_w<{360}^{\circ }$ and corresponded to saturation duty above the 35% classification threshold. No controller reached the ${45}^{\circ }$ individual-axis excursion boundary or the ${80}^{\circ }$ near-flip boundary.

Figure 11 presents the maximum geometric tilt over wind azimuth. The recorded maxima were $39.{997}^{\circ }$, $18.{251}^{\circ }$, $2.{225}^{\circ }$, and $0.{046}^{\circ }$ for PID, LQR, ${\mathrm{H}}_{\infty }$, and SMC, respectively.

The maximum individual-axis roll–pitch excursions were $31.{141}^{\circ }$, $16.{454}^{\circ }$, $2.{214}^{\circ }$, and $0.{040}^{\circ }$ for PID, LQR, ${\mathrm{H}}_{\infty }$, and SMC, respectively, as shown in Figure 12.

The directional geometric-tilt envelopes are shown in Figure 13. The worst directions reported in Table 10 were ${15}^{\circ }$ for PID and LQR, ${90}^{\circ }$ for ${\mathrm{H}}_{\infty }$, and ${45}^{\circ }$ for SMC.

Figure 14 presents the differential attitude-control effort and saturation duty. The mean differential-effort values were 53.921, 59.003, 59.805, and 59.836 for PID, LQR, ${\mathrm{H}}_{\infty }$, and SMC, respectively. The maximum saturation duties were 38.452%, 15.964%, 11.971%, and 9.222%, respectively.

The stability and recovery classifications are shown in Figure 15. All controllers satisfied the roll–pitch recovery criterion over the complete sweep. The worst recovery times were $0.481\mathrm{s}$ for PID, $0.202\mathrm{s}$ for LQR, $1.95\times {10}^{-4}\mathrm{s}$ for ${\mathrm{H}}_{\infty }$, and $1.61\times {10}^{-4}\mathrm{s}$ for SMC. No near-flip event was recorded at the tested disturbance magnitude.

4. Discussion

4.1. Synthesis of the Principal Findings

The framework of Section 2 was designed to isolate controller structure as the principal comparison variable. All four inner-loop controllers act on the same nonlinear plant of Equations (26)–(29), use the same parameters in Table 2, share the outer-loop architecture of Figure 3, and are subject to the common actuator and command limits of Equations (Section 2.1.8), (36), and (37). They are also tuned with the common objective of Equations (93)–(98) over the admissible ranges in Table 3, using budget-matched RCGA and PSO settings as defined in Section 2.5.3. PSO returned the lower objective value for all controller families (Table 6); therefore, the subsequent comparisons use the PSO-selected configurations. Because the optimizer comparison is based on a controlled single run, it is interpreted as within-study evidence rather than a statistical claim of universal superiority [48].

The principal result is that controller separation depends strongly on the benchmark. Nominal step tracking, D1, and D2 show relatively small differences among the optimized controllers (Table 7, Table 8, and Table 9), whereas the D3 directional sweep separates them clearly. Stable directional coverage is $58.33\%$ for PID and $100\%$ for LQR, ${\mathcal{H}}_{\infty }$, and SMC (Table 10, Figure 15). The worst geometric tilt also ranges widely, from $39.{997}^{\circ }$ for PID to $0.{046}^{\circ }$ for SMC (Figure 11). Thus, the complete directional sweep of Equation (176) reveals differences that fixed-direction tests do not.

The failure mechanism is equally important. The D3 classifier of Equations (195) and (196) shows that no controller exceeded the ${45}^{\circ }$ individual-axis excursion limit or approached the ${80}^{\circ }$ near-flip threshold. PID’s ten failed directions resulted only from saturation duty exceeding the $35\%$ protocol limit, while its maximum individual-axis excursion remained $31.{141}^{\circ }$ (Figure 12). The study therefore characterizes directional robustness margins and saturation-based protocol failure rather than demonstrating an actual flip.

4.2. Nominal Response and the Finite-Duration Force Disturbance (D1)

The nominal results confirm that all optimized configurations achieve fast and well-damped tracking. PID–SMC produced the shortest position rise and settling times and the lowest aggregate position IAE and ITAE, whereas PID–PID and PID–${\mathcal{H}}_{\infty }$ produced the smallest overshoots (Table 7, Figure 5). These differences are consistent with the controller structures: the SMC law of Equation (88) provides strong corrective action near the sliding manifold, while the ${\mathcal{H}}_{\infty }$ formulation weights attitude, rate, and control activity through Equation (78) [40,42,43].

D1 applies a translational disturbance force to the full 12-state model with $r_{\mathrm{cp}}=0$ (Table 5); therefore, the attitude response arises indirectly through the outer-loop roll and pitch commands of Equations (59a) and (). The optimized controllers all recovered promptly after the disturbance, while the baseline response exhibited substantially larger accumulated error and longer settling (Table 8, Figure 6). Because the optimized responses converge to similar values once the disturbance window ends, D1 confirms effective transient-force rejection but only weakly discriminates among the optimized inner-loop designs. This motivates the direct-moment directional protocol of Section 2.8.

4.3. Helical Tracking Under Wind (D2)

D2 evaluates the same controllers during sustained coupled three-dimensional motion. The trajectory is defined in Section 2.7.2, and the resulting performance is summarized in Table 9. SMC achieved the lowest position RMSE, reducing it from $1.0790$ m for PID to $0.6768$ m, with the same ordering retained during the wind window. Since the D2 force enters the translational dynamics rather than the rotational input channels, the resulting SMC advantage reflects the complete cascaded configuration rather than only the matched-disturbance property of Equation (91) [43,44,47].

LQR achieved the lowest normalized D2 control effort defined by Equation (201), while SMC required slightly greater actuation for its improved tracking accuracy (Table 9, Figure 10). This is consistent with the LQR cost of Equation (71), which explicitly penalizes control through R, whereas the SMC switching contribution in Equation (88) maintains corrective action within the boundary-layer approximation of Equation (89) [37,38].

The D2 maximum geometric tilt is nearly identical for all controllers, spanning only $30.{417}^{\circ }$–$30.{616}^{\circ }$. This similarity should not be interpreted as equal attitude robustness because all configurations share the $\pm {25}^{\circ }$ roll/pitch command limits of Equations (37) and (62). The common command shaping therefore constrains the D2 tilt response and limits its ability to distinguish the inner-loop controllers. The direct-moment D3 benchmark was introduced specifically to expose these differences.

4.4. Directional Roll–Pitch Robustness (D3)

D3 differs from D1 and D2 in disturbance mechanism. The force defined by Equation (178) acts at the offset centre of pressure of Equation (180) and generates the direct moment of Equation (181). This moment is applied to the attitude-only model of Equation (183), with $U_1=mg$ and no propagated translational states. As stated in Section 2.8.6, D3 tilt, rate, effort, saturation, and recovery quantities are therefore complementary to, rather than directly comparable with, the D1 and D2 tracking metrics.

The D3 moment geometry explains the principal directional pattern. Proposition 1 in Section 2.6.6 gives the necessary static-trim condition of Equation (149). For the implemented geometry, Equation (182) gives maximum roll and pitch disturbance moments of $0.48$ N·m, below the corresponding $1.766$ N·m limits in Equation (150); the maximum yaw disturbance moment is $0.905$ N·m, below the $1.177$ N·m yaw limit. Hence, static trim remains feasible for every tested direction. However, the yaw channel reaches $76.9\%$ of its available authority, whereas roll and pitch reach only $27.2\%$, making yaw the most heavily loaded channel.

This asymmetry is reflected in Figure 14. From Equation (182), the yaw disturbance component is proportional to $sin(\alpha -{45}^{\circ })$, so it vanishes at ${45}^{\circ }$ and ${225}^{\circ }$ and reaches its largest magnitude at ${135}^{\circ }$ and ${315}^{\circ }$. The saturation-duty curves follow the same pattern, and PID’s failed sectors, ${105}^{\circ }\le {\alpha }_w<{180}^{\circ }$ and ${285}^{\circ }\le {\alpha }_w<{360}^{\circ }$, are centred around these high-yaw-loading directions. Thus, the observed saturation boundary is primarily a yaw-loading signature rather than a roll–pitch tilt boundary.

The attitude envelopes exhibit a different directional pattern. Figure 11, Figure 12, and Figure 13 show that the largest roll–pitch excursions occur near directions where the induced roll or pitch moments are largest, while the largest saturation duty occurs near the yaw-moment maxima. Consequently, the critical direction depends on the metric: the direction that maximizes attitude excursion is not the same as the direction that maximizes actuator saturation. This distinction is one of the main advantages of the full ${360}^{\circ }$ sweep.

The coverage results in Table 10 are consistent with this mechanism. LQR, ${\mathcal{H}}_{\infty }$, and SMC retained $100\%$ stable coverage, whereas PID exceeded the $35\%$ saturation-duty threshold in ten directions. Because Equation (149) remains satisfied at every azimuth, these failures cannot be attributed to loss of static control authority. A plausible contributor is the absence of explicit anti-windup in the implemented PID law: its integral state in Equation (65a) continues to evolve during actuator saturation, which can prolong limit contact after the disturbance has decayed [34,35]. This interpretation is consistent with the saturation-duty metric of Equation (175) and should be treated as a mechanism supported by the implemented model rather than as a universal property of PID control.

The angular-rate results provide further evidence that yaw dominates the full body-rate response. The roll–pitch rate metric of Equation (186) remains far below the 30 rad·s⁻¹ classification threshold for every controller, while the component values in Table 10(b) show much larger yaw-rate maxima. Because yaw is intentionally excluded from the roll–pitch failure classifier of Equation (196), these high yaw rates are diagnostic rather than failure-triggering quantities. This reinforces that the present D3 protocol evaluates roll–pitch integrity under a disturbance geometry that loads yaw most strongly.

Control effort and saturation duty likewise measure different phenomena. PID recorded the lowest mean differential effort, while SMC recorded the highest, yet PID also had the largest saturation duty and SMC the smallest (Table 10, Figure 14). The differential effort of Equation (173) measures integrated control activity, whereas Equation (175) measures time spent near an actuator bound. Therefore, effort alone cannot characterize actuator-limited robustness, and both quantities are required for meaningful interpretation.

Stability classification and post-disturbance recovery also represent distinct properties. The stable and recovered sets are defined independently by Equations (199a) and (). PID failed the stability protocol at ten of twenty-four directions but recovered in all twenty-four, and all four controllers achieved full recovery coverage (Figure 15). The recovery criterion of Equation (198) therefore should not be interpreted as equivalent to the complete stability classifier. The nearly immediate recovery values of ${\mathcal{H}}_{\infty }$ and SMC indicate that their roll–pitch states were already within the recovery band when the disturbance ended.

4.5. Comparison with Related Work, Practical Implications, and Scope

The contribution of the present work lies primarily in the evaluation protocol rather than in proposing a new controller. Previous comparative studies have assessed quadrotor controllers using nominal stabilization, trajectory tracking, and standard transient metrics [4,5,6], while wind-rejection studies have demonstrated effective disturbance attenuation using ${\mathcal{H}}_{\infty }$, observer-based, and sliding-mode approaches [14,18,27,28,68]. These studies generally consider a fixed direction, a small set of wind cases, or prescribed disturbance profiles. The present framework extends this literature by resolving the disturbance azimuth explicitly through the grid of Equation (176), reporting stable and failed sectors through the coverage definition of Equation (200), and linking the observed boundaries to attitude excursion, recovery, and actuator saturation. The resulting directional envelopes in Figure 13 and the stability maps in Figure 15 therefore provide information that cannot be obtained from a single-direction test.

Three practical implications follow within the stated simulation scope. First, directional coverage must always be reported together with the actuator constraints that produce it; the coverage values in Table 10 are conditional on the bounds of Equation (Section 2.1.8) and the parameters of Table 2. Second, the present platform has lower yaw moment authority than roll/pitch authority, as quantified by Equations (31) and (32); therefore, an off-axis force applied at an offset centre of pressure can become a yaw-saturation problem before it becomes a roll–pitch trim problem. Third, controller selection should be criterion-dependent: SMC provides the strongest attitude regulation in the present sweep, LQR gives the lowest D2 control effort, and ${\mathcal{H}}_{\infty }$ combines small directional excursions with rapid recovery, while PID requires sufficient saturation margin or additional anti-windup protection.

The conclusions remain conditional on the implemented model and protocol. The simulations assume continuous-time full-state feedback without sensor noise, estimation error, computation delay, or battery dynamics, and actuator dynamics are represented by the bounded virtual-input model of Equation (33). The disturbance is smoothly gated by Equation (44) and uses the scenario definitions of Table 5. D1 and D2 apply translational forces with $r_{\mathrm{cp}}=0$, whereas D3 applies the induced moment to the attitude-only plant; accordingly, their numerical metrics should not be compared directly. In addition, the ${15}^{\circ }$ grid of Equation (176) limits directional boundary resolution to one sampling interval, and the classifier thresholds of Equation (195) are protocol-dependent. The ${45}^{\circ }$ safe-excursion limit is consistent with established multirotor practice [66,67], whereas the saturation-duty, rate, and recovery thresholds are study-defined.

Within these limits, four conclusions are supported by the reported evidence: directional robustness distinguishes the optimized controllers more clearly than nominal tracking; the observed PID failures arise from sustained actuator saturation rather than unsafe roll–pitch excursion or loss of static trim authority; stability classification and post-disturbance recovery are distinct properties; and the critical direction depends on whether attitude excursion or actuator usage is used as the governing metric.

5. Conclusions and Future Work

This study presented a unified framework for comparing RCGA- and PSO-optimized PID, LQR, ${\mathcal{H}}_{\infty }$, and SMC quadrotor controllers under common plant dynamics, constraints, and disturbance scenarios. The controllers were evaluated under disturbed step tracking, wind-disturbed helical tracking, and a complete $0^{\circ }$–${360}^{\circ }$ directional attitude-robustness sweep at ${15}^{\circ }$ resolution.

The directional assessment provided the principal findings. PID achieved $58.33\%$ stable coverage, whereas LQR, ${\mathcal{H}}_{\infty }$, and SMC maintained $100\%$. The largest attitude excursions occurred near ${15}^{\circ }$, while the highest saturation demand occurred near ${135}^{\circ }$ and ${315}^{\circ }$, where the induced yaw moment was greatest. SMC provided the strongest directional attitude regulation, LQR offered the most balanced tracking–effort trade-off, and ${\mathcal{H}}_{\infty }$ exhibited low directional sensitivity and rapid recovery. PID was the least directionally robust, with failures caused exclusively by the $35\%$ saturation-duty criterion rather than unsafe tilt, excessive roll–pitch rate, or failed recovery. Thus, no controller was superior across all criteria.

The results further demonstrate that directional stability, recovery, control effort, and actuator saturation should be assessed separately. All controllers achieved full recovery coverage despite PID failing the stability protocol in ten of twenty-four directions, while lower control effort did not necessarily imply lower saturation susceptibility. No response approached the prescribed near-flip boundary; therefore, the study quantifies directional robustness margins and saturation-based failure mechanisms rather than actual flip or flip recovery. Overall, the proposed directional framework complements conventional tracking benchmarks by revealing controller limitations and actuator-critical disturbance directions that may remain hidden in fixed-direction evaluations.

Several extensions follow directly from the boundaries of the present study. Experimental and hardware-in-the-loop validation is the most important, since the evidence is simulation-based under idealized sensing and full-state feedback; the protocol of Section 2.8 transfers unchanged to a hardware-in-the-loop bench. The disturbance model should be broadened to multiple magnitudes, durations, and profiles, including stochastic turbulence spectra and aerodynamic uncertainty, and to a magnitude sweep that estimates the critical disturbance level at which flip-like behaviour begins, a quantity the present fixed-magnitude protocol cannot determine. Moreover, the plant can be extended with first-order rotor and motor dynamics, battery-voltage droop, sensor noise, estimation error, and computation delays, and evaluated under parametric uncertainty in mass, inertia, and aerodynamic coefficients, as well as under motor faults and actuator degradation.

Appendix A.1. Baseline and Optimized Controller Parameters

This appendix reports the complete baseline, GA-optimized, and PSO-optimized controller parameter sets used in the study. The outer-loop PID gains are denoted by $K_{p,o}$, $K_{i,o}$, and $K_{d,o}$, whereas the corresponding inner-loop gains are denoted by $K_{p,i}$, $K_{i,i}$, and $K_{d,i}$. The derivative-filter bandwidth and time constant satisfy $T_f=1/N_f$, subject to the numerical precision shown in the tables. The PSO parameter sets were retained for the subsequent evaluations because they produced the lowest objective-function value within every controller architecture. Table A1, Table A2, Table A3, and Table A4 report the complete parameter sets for the four controller architectures.

Table A1. Baseline and optimized parameters of the PID–PID controller.

Loop	Method	${\mathit{K}}_{\mathit{p}}$	${\mathit{K}}_{\mathit{i}}$	${\mathit{K}}_{\mathit{d}}$	${\mathit{N}}_{\mathit{f}}$(s−1)	${\mathit{T}}_{\mathit{f}}$ (s)
Outer position loop	Baseline (pidtune)	2.407	0.451	2.854	16.360	0.0611
	GA optimized	23.980	7.013	12.950	268.170	0.0037
	PSO optimized	25.998	0.100	16.371	300.000	0.0033
Inner attitude loop	Baseline (pidtune)	3.749	4.694	0.059	7.985	0.1250
	GA optimized	11.865	2.951	0.522	189.620	0.0053
	PSO optimized	9.765	0.123	1.399	10.701	0.0930

Note: The same outer-loop gain triplet is applied to the three position channels, and the same inner-loop gain triplet is applied to the three attitude channels, as defined in the controller formulation.

Table A1. Baseline and optimized parameters of the PID–PID controller.

Loop	Method	${\mathit{K}}_{\mathit{p}}$	${\mathit{K}}_{\mathit{i}}$	${\mathit{K}}_{\mathit{d}}$	${\mathit{N}}_{\mathit{f}}$(s−1)	${\mathit{T}}_{\mathit{f}}$ (s)
Outer position loop	Baseline (pidtune)	2.407	0.451	2.854	16.360	0.0611
	GA optimized	23.980	7.013	12.950	268.170	0.0037
	PSO optimized	25.998	0.100	16.371	300.000	0.0033
Inner attitude loop	Baseline (pidtune)	3.749	4.694	0.059	7.985	0.1250
	GA optimized	11.865	2.951	0.522	189.620	0.0053
	PSO optimized	9.765	0.123	1.399	10.701	0.0930

Note: The same outer-loop gain triplet is applied to the three position channels, and the same inner-loop gain triplet is applied to the three attitude channels, as defined in the controller formulation.

Table A2. Baseline and optimized parameters of the PID–LQR controller.

Parameter	Baseline	GA optimized	PSO optimized
$K_{p,o}$	20.000	20.000	30.010
$K_{i,o}$	1.200	0.240	0.100
$K_{d,o}$	11.000	10.950	12.850
$q_{\varphi }$	130.000	133.400	191.460
$q_{\theta }$	180.000	186.960	199.910
$q_{\psi }$	110.000	112.260	52.100
$q_p$	2.000	2.740	0.160
$q_q$	4.000	0.470	0.100
$q_r$	5.000	5.280	1.980
$r_{u_2}$	0.300	0.290	0.160
$r_{u_3}$	1.800	0.400	0.580
$r_{u_4}$	1.500	1.190	0.350

Table A2. Baseline and optimized parameters of the PID–LQR controller.

Parameter	Baseline	GA optimized	PSO optimized
$K_{p,o}$	20.000	20.000	30.010
$K_{i,o}$	1.200	0.240	0.100
$K_{d,o}$	11.000	10.950	12.850
$q_{\varphi }$	130.000	133.400	191.460
$q_{\theta }$	180.000	186.960	199.910
$q_{\psi }$	110.000	112.260	52.100
$q_p$	2.000	2.740	0.160
$q_q$	4.000	0.470	0.100
$q_r$	5.000	5.280	1.980
$r_{u_2}$	0.300	0.290	0.160
$r_{u_3}$	1.800	0.400	0.580
$r_{u_4}$	1.500	1.190	0.350

Table A3. Baseline and optimized parameters of the PID–${\mathrm{H}}_{\infty }$ controller.

Parameter	Baseline	GA optimized	PSO optimized
$K_{p,o}$	20.000	39.790	40.000
$K_{i,o}$	1.200	0.235	0.100
$K_{d,o}$	11.000	16.244	16.084
$w_{a,\varphi }$	130.000	276.920	400.000
$w_{a,\theta }$	180.000	155.760	319.360
$w_{a,\psi }$	110.000	29.940	211.040
$w_{r,p}$	2.000	1.380	0.100
$w_{r,q}$	4.000	18.810	20.000
$w_{r,r}$	5.000	19.480	10.370
$r_{u_2}$	0.300	0.818	3.000
$r_{u_3}$	1.800	1.779	1.935
$r_{u_4}$	1.500	1.601	1.825

Table A3. Baseline and optimized parameters of the PID–${\mathrm{H}}_{\infty }$ controller.

Parameter	Baseline	GA optimized	PSO optimized
$K_{p,o}$	20.000	39.790	40.000
$K_{i,o}$	1.200	0.235	0.100
$K_{d,o}$	11.000	16.244	16.084
$w_{a,\varphi }$	130.000	276.920	400.000
$w_{a,\theta }$	180.000	155.760	319.360
$w_{a,\psi }$	110.000	29.940	211.040
$w_{r,p}$	2.000	1.380	0.100
$w_{r,q}$	4.000	18.810	20.000
$w_{r,r}$	5.000	19.480	10.370
$r_{u_2}$	0.300	0.818	3.000
$r_{u_3}$	1.800	1.779	1.935
$r_{u_4}$	1.500	1.601	1.825

Table A4. Baseline and optimized parameters of the PID–SMC controller.

Parameter	Baseline	GA optimized	PSO optimized
$K_{p,o}$	12.500	35.3028	40.0000
$K_{i,o}$	7.500	0.12719	0.10000
$K_{d,o}$	8.500	14.3953	15.0651
$\beta$	27.000	56.9371	58.4464
$K_{d,\mathrm{smc}}$	530.000	442.473	730.898
a	0.100	0.194541	0.020000

Table A4. Baseline and optimized parameters of the PID–SMC controller.

Parameter	Baseline	GA optimized	PSO optimized
$K_{p,o}$	12.500	35.3028	40.0000
$K_{i,o}$	7.500	0.12719	0.10000
$K_{d,o}$	8.500	14.3953	15.0651
$\beta$	27.000	56.9371	58.4464
$K_{d,\mathrm{smc}}$	530.000	442.473	730.898
a	0.100	0.194541	0.020000

Appendix A.2. Axis-Wise Performance Metrics for the D1 Scenario

Table A5 and Table A6 report the individual axis-wise metrics used to construct the aggregated D1 results in Table 8. For each response group, the aggregated rise time, settling time, and overshoot are the maximum values across the corresponding axes, whereas the aggregated IAE and ITAE are obtained by summing the individual axis values. The steady-state error (SSE) is reported only in the present appendix because it is not included in the aggregated comparison.

Table A5. Axis-wise position-response metrics under the D1 finite-duration wind disturbance.

Controller	Axis	${\mathit{M}}_{\mathit{p}}$ (%)	${\mathit{T}}_{\mathit{r}}$ (s)	${\mathit{T}}_{\mathit{s}}$ (s)	SSE (m)	IAE	ITAE
PID–PID baseline	x	46.1	1.516	31.626	$-8.795\times {10}^{-2}$	9.3262	119.9020
	y	16.2	1.610	30.442	$\phantom{-}3.136\times {10}^{-2}$	3.8159	43.6879
	z	18.9	1.306	28.789	$-1.992\times {10}^{-2}$	2.9151	29.8242
PID–PID PSO	x	5.8	0.973	18.350	$\phantom{-}5.905\times {10}^{-3}$	2.6000	10.7412
	y	0.3	1.013	18.094	$\phantom{-}3.203\times {10}^{-3}$	0.9801	3.7690
	z	2.2	1.149	18.098	$\phantom{-}1.976\times {10}^{-3}$	0.7783	2.7292
PID–LQR PSO	x	5.0	0.705	18.309	$\phantom{-}4.258\times {10}^{-3}$	2.1657	8.4541
	y	0.3	0.745	18.076	$\phantom{-}2.415\times {10}^{-3}$	0.8519	3.0754
	z	1.9	0.674	1.165	$\phantom{-}1.191\times {10}^{-3}$	0.5792	1.9022
PID–${\mathrm{H}}_{\infty }$ PSO	x	3.8	0.732	18.227	$\phantom{-}3.612\times {10}^{-3}$	2.1288	6.8467
	y	0.2	0.730	17.888	$\phantom{-}1.741\times {10}^{-3}$	0.7891	2.3299
	z	1.4	0.676	1.185	$\phantom{-}9.607\times {10}^{-4}$	0.5439	1.4830
PID–SMC PSO	x	3.8	0.679	18.202	$\phantom{-}3.492\times {10}^{-3}$	2.1076	6.7369
	y	0.4	0.678	17.893	$\phantom{-}1.699\times {10}^{-3}$	0.7762	2.3005
	z	1.6	0.629	1.086	$\phantom{-}9.067\times {10}^{-4}$	0.5311	1.4449

Table A5. Axis-wise position-response metrics under the D1 finite-duration wind disturbance.

Controller	Axis	${\mathit{M}}_{\mathit{p}}$ (%)	${\mathit{T}}_{\mathit{r}}$ (s)	${\mathit{T}}_{\mathit{s}}$ (s)	SSE (m)	IAE	ITAE
PID–PID baseline	x	46.1	1.516	31.626	$-8.795\times {10}^{-2}$	9.3262	119.9020
	y	16.2	1.610	30.442	$\phantom{-}3.136\times {10}^{-2}$	3.8159	43.6879
	z	18.9	1.306	28.789	$-1.992\times {10}^{-2}$	2.9151	29.8242
PID–PID PSO	x	5.8	0.973	18.350	$\phantom{-}5.905\times {10}^{-3}$	2.6000	10.7412
	y	0.3	1.013	18.094	$\phantom{-}3.203\times {10}^{-3}$	0.9801	3.7690
	z	2.2	1.149	18.098	$\phantom{-}1.976\times {10}^{-3}$	0.7783	2.7292
PID–LQR PSO	x	5.0	0.705	18.309	$\phantom{-}4.258\times {10}^{-3}$	2.1657	8.4541
	y	0.3	0.745	18.076	$\phantom{-}2.415\times {10}^{-3}$	0.8519	3.0754
	z	1.9	0.674	1.165	$\phantom{-}1.191\times {10}^{-3}$	0.5792	1.9022
PID–${\mathrm{H}}_{\infty }$ PSO	x	3.8	0.732	18.227	$\phantom{-}3.612\times {10}^{-3}$	2.1288	6.8467
	y	0.2	0.730	17.888	$\phantom{-}1.741\times {10}^{-3}$	0.7891	2.3299
	z	1.4	0.676	1.185	$\phantom{-}9.607\times {10}^{-4}$	0.5439	1.4830
PID–SMC PSO	x	3.8	0.679	18.202	$\phantom{-}3.492\times {10}^{-3}$	2.1076	6.7369
	y	0.4	0.678	17.893	$\phantom{-}1.699\times {10}^{-3}$	0.7762	2.3005
	z	1.6	0.629	1.086	$\phantom{-}9.067\times {10}^{-4}$	0.5311	1.4449

Table A6. Axis-wise attitude-response metrics under the D1 finite-duration wind disturbance.

Controller	Angle	${\mathit{M}}_{\mathit{p}}$ (%)	${\mathit{T}}_{\mathit{r}}$ (s)	${\mathit{T}}_{\mathit{s}}$ (s)	SSE (rad)	IAE (rad)	ITAE (rad s)
PID–PID baseline	$\varphi$	—	0.000	20.972	$-3.010\times {10}^{-5}$	0.7337	10.9176
	$\theta$	—	0.000	21.082	$-4.478\times {10}^{-5}$	1.3937	16.3023
	$\psi$	22.240	0.074	0.916	$\phantom{-}1.106\times {10}^{-8}$	0.0377	0.0885
PID–PID PSO	$\varphi$	—	0.000	18.527	$-1.625\times {10}^{-7}$	0.8358	9.0812
	$\theta$	—	0.000	18.570	$\phantom{-}6.401\times {10}^{-7}$	1.2961	13.4645
	$\psi$	4.370	0.436	0.999	$\phantom{-}3.993\times {10}^{-4}$	0.0645	0.2907
PID–LQR PSO	$\varphi$	—	0.000	18.182	$-1.032\times {10}^{-7}$	0.8796	8.9917
	$\theta$	—	0.000	18.182	$\phantom{-}3.975\times {10}^{-7}$	1.4395	13.4589
	$\psi$	3.924	0.420	18.182	$\phantom{-}5.968\times {10}^{-14}$	0.0959	0.1864
PID–${\mathrm{H}}_{\infty }$ PSO	$\varphi$	—	0.000	18.206	$-4.714\times {10}^{-8}$	0.8787	8.9274
	$\theta$	—	0.000	18.228	$\phantom{-}2.487\times {10}^{-7}$	1.3987	13.3085
	$\psi$	6.213	0.083	0.941	$-1.221\times {10}^{-15}$	0.0256	0.0381
PID–SMC PSO	$\varphi$	—	0.000	18.206	$-4.654\times {10}^{-8}$	0.8943	8.9528
	$\theta$	—	0.000	18.228	$\phantom{-}2.407\times {10}^{-7}$	1.4259	13.3550
	$\psi$	1.923	0.067	0.905	$-8.610\times {10}^{-17}$	0.0198	0.0160

Note: Overshoot is not reported for roll and pitch because their reference commands are zero; hence, the conventional percentage overshoot definition is not applicable. Their reported rise times are zero for the same reason. The long roll–pitch settling times primarily reflect the finite-duration disturbance window and subsequent recovery rather than nominal command-following dynamics.

Table A6. Axis-wise attitude-response metrics under the D1 finite-duration wind disturbance.

Controller	Angle	${\mathit{M}}_{\mathit{p}}$ (%)	${\mathit{T}}_{\mathit{r}}$ (s)	${\mathit{T}}_{\mathit{s}}$ (s)	SSE (rad)	IAE (rad)	ITAE (rad s)
PID–PID baseline	$\varphi$	—	0.000	20.972	$-3.010\times {10}^{-5}$	0.7337	10.9176
	$\theta$	—	0.000	21.082	$-4.478\times {10}^{-5}$	1.3937	16.3023
	$\psi$	22.240	0.074	0.916	$\phantom{-}1.106\times {10}^{-8}$	0.0377	0.0885
PID–PID PSO	$\varphi$	—	0.000	18.527	$-1.625\times {10}^{-7}$	0.8358	9.0812
	$\theta$	—	0.000	18.570	$\phantom{-}6.401\times {10}^{-7}$	1.2961	13.4645
	$\psi$	4.370	0.436	0.999	$\phantom{-}3.993\times {10}^{-4}$	0.0645	0.2907
PID–LQR PSO	$\varphi$	—	0.000	18.182	$-1.032\times {10}^{-7}$	0.8796	8.9917
	$\theta$	—	0.000	18.182	$\phantom{-}3.975\times {10}^{-7}$	1.4395	13.4589
	$\psi$	3.924	0.420	18.182	$\phantom{-}5.968\times {10}^{-14}$	0.0959	0.1864
PID–${\mathrm{H}}_{\infty }$ PSO	$\varphi$	—	0.000	18.206	$-4.714\times {10}^{-8}$	0.8787	8.9274
	$\theta$	—	0.000	18.228	$\phantom{-}2.487\times {10}^{-7}$	1.3987	13.3085
	$\psi$	6.213	0.083	0.941	$-1.221\times {10}^{-15}$	0.0256	0.0381
PID–SMC PSO	$\varphi$	—	0.000	18.206	$-4.654\times {10}^{-8}$	0.8943	8.9528
	$\theta$	—	0.000	18.228	$\phantom{-}2.407\times {10}^{-7}$	1.4259	13.3550
	$\psi$	1.923	0.067	0.905	$-8.610\times {10}^{-17}$	0.0198	0.0160

Note: Overshoot is not reported for roll and pitch because their reference commands are zero; hence, the conventional percentage overshoot definition is not applicable. Their reported rise times are zero for the same reason. The long roll–pitch settling times primarily reflect the finite-duration disturbance window and subsequent recovery rather than nominal command-following dynamics.