Robustness Analysis of the Model-Predictive Position Control of an Electro-Mechanical Actuator for Primary Flight Surfaces

1. Introduction

In recent years, the aerospace sector has undergone a progressive shift toward the electrification of onboard systems, driven by the need to improve energy efficiency, reduce environmental impact, and optimize system integration [1]. The More Electric Aircraft (MEA) concept seeks to replace traditional hydraulic and pneumatic systems with electrically powered alternatives, offering benefits in terms of weight reduction, maintenance costs, and sustainability [2,3,4]. Nevertheless, given the present power density storage technology [5], the feasibility of achieving the All Electric Aircraft (AEA) is currently limited to Urban Air Mobility (UAM) applications, such as the emerging market of Urban Air Taxi (UAT) and electrical Vertical Take-Off and Landing (eVTOL) vehicles [6]. A key enabling technology within this transformation is the Electro-Mechanical Actuator [7], which converts electrical energy into mechanical motion using an electric motor and a mechanical transmission, hence without relying on hydraulic fluids, eliminating leakage risks and simplifying system architecture. Significant progress in enabling components—such as electric machines [8], power electronics [9], and materials of mechanical elements [10]—has notably increased the Technology Readiness Level (TRL) of these systems. Although the feasibility of EMAs in aerospace applications has been demonstrated in terms of load capacity, dynamic performance, and health monitoring capabilities [11,12,13,14,15], the deployment on safety-critical flight controls or landing gears still raises concerns related to reliability, fault tolerance, and compliance with airworthiness certification standards [16,17]. To address these functions, several architectural solutions have been proposed, including simplex configurations with passive fail-safe mechanisms and more complex redundant fault-tolerant architectures involving dual motors, independent control paths, or differential mechanisms. Nevertheless, actuator reliability depends not only on hardware redundancy but also on the control and monitoring strategies implemented [18]. Robust and high-performance controllers are essential to ensure system stability and disturbance rejection under variable conditions and partial faults [19]. Classical control approaches, such as Proportional-Integral (PI) or Linear Quadratic Regulator (LQR) schemes, are widely employed in aerospace applications due to their simplicity, reliability and deterministic behaviour. However, these methods operate on current or past states and measurements, hence limited in managing future disturbances or constraints. In contrast, Model Predictive Control has emerged as a powerful alternative, enabling the anticipation of future system behaviour through internal prediction models and the explicit handling of input and state constraints within the control action optimization window [20]. Eren et al. [22] provide a comprehensive review of MPC techniques for aircraft and spacecraft systems applications, focusing on safety-critical applications and fault-tolerant control, observing that the complexity of nonlinear formulations is a challenge for practical implementation in safety-critical functions. Although RTCA DO-178C does not restrict the choice of control techniques from a certification standpoint [23], the software implementation of PI controllers is generally easier to validate due to their interpretability and limited branching logic, making Modified Condition/Decision Coverage (MC/DC) relatively straightforward. In contrast, for MPC involving model nonlinearities or constrained real-time optimization, the highly complex logic, numerical challenges, and recursive solver behaviour may significantly increase the verification effort—particularly in applications subject to higher Development Assurance Level (DAL) classification.

MPC has been studied in the context of aircraft flight mechanics for load alleviation functions [24] and dynamic system modelling [25], spacecraft enhancement dynamics [26], and UAVs manoeuvres optimization [27,28]. MPC has also demonstrated promising results in governing Permanent Magnet Synchronous Motors (PMSM), functioning as a speed controller [29], MIMO Direct-Quadrature currents regulator [30], and direct torque controller (DTC) [31]. Additionally, Hasan et al. [32] employed MPC as a speed regulator for electromagnetic actuators, integrated with a PI position controller. Nevertheless, the use of MPC for position control in safety-critical EMAs remains a promising yet open research topic. The main barriers are associated with the nonlinear and coupled dynamics of electro-mechanical systems, the need for real-time implementation on embedded hardware with limited computational resources, and the rigorous certification requirements in aviation environments.

The main objective of this work is to assess the feasibility and advantages of applying MPC to a rotary EMA developed by UMBRAGROUP for primary flight control applications, by comparing its performance and robustness with those of conventional control laws.

The paper is articulated as follows: the first part is dedicated to the description of the EMA architecture and to the development of a detailed model integrating both electrical and mechanical subsystems from first physical principles. Subsequently, the classical PI-based cascade of position, speed and current regulators is designed through a multi-objectives optimization. A Linear Time Invariant (LTI) MPC then replaces the PI position regulator, while keeping other loops unchanged. Comprehensive simulations with a nonlinear model in Simulink assess the performance of PI- and MPC-controlled EMAs in time and frequency domains, also evaluating robustness against individual and coupled parameters deviations and model uncertainties, as well as the power absorption and computational resource requirements. Finally, key findings are discussed, highlighting differences between the two approaches, along with suggestions for future refinements.

2. Materials and Methods

2.1. System Description

The reference architecture for this study is a rotary fail-safe EMA (Figure 1), based on a differential ball-screw patented by UMBRAGROUP (UG) [33], within the CleanSky 2 program [34]. The closed-loop control of EMAs designed for secondary flight surfaces primarily focuses on disturbance rejection requirements, minimizing position deviations from the commanded setpoint when external loads act on the surface. Conversely, the position tracking behaviour under load is the key performance requirement that an actuator specifically designed for primary flight controls must meet. In this last case, the position control loop should operate at higher bandwidths (e.g., 2–10 Hz), ensuring precise tracking, while the output shaft speed typically ranges between 10 and 60 degrees per second.

The specific case study investigates the rotary EMA architecture integrated into the rudder actuation system of a medium-sized eVTOL vehicle featuring a V-tail configuration. In such a layout, the rudder functions as a primary flight control surface, contributing to both yaw and pitch control during cruise and transition phases. Consequently, the actuator must ensure high dynamic performance, particularly in terms of frequency response, position tracking accuracy, and robustness against aerodynamic disturbances.

2.2. System Dynamics Modelling

High-fidelity nonlinear models of EMAs are widely used to accurately reproduce the behaviour of the real system. Such model-based approaches have become increasingly relevant in the context of digital twin development, enabling simulation-based analysis and verification for monitoring and control algorithms. The developed nonlinear model of the EMA, illustrated in Figure 2, includes two main sections: the electro-mechanical section and the digital electronics section, each comprising multiple functional blocks. The electro-mechanical section consists of the following components:

• Three-phase Surface mounted Permanent Magnet Synchronous Motor (S-PMSM) current dynamics in the direct-quadrature (d-q) reference frame;
• UMBRAGROUP patented mechanical transmission based on a high transmission ratio differential ball-screw mechanism, including its internal stiffness and resonance;
• Coulomb and viscous friction;
• Rotational dynamics of the output shaft, including the reference aircraft surface inertia and external load effects.

The digital electronics section includes:

• Three nested control loops acting on motor direct and quadrature currents, motor speed, and output shaft position;
• Control nonlinearities (saturations, rate-limiter on position command);
• Sensors dynamics, latencies and resolution;
• Resolver sensors for motor and output shaft position feedback;
• Motor speed reconstruction from angular position;
• Hall-effect sensors for current measurement;
• Digital signal processing performed at 10 kHz.

The model is implemented in Simulink using the ODE4 Runge-Kutta solver with a fixed integration step of 10^-5 s.

2.2.1. Electro-Mechanical Section

The electro-mechanical section includes three main dynamic contributions: the electric motor with d-q axes current dynamics (1), the rotor-side one degree of freedom (DOF) rotational motion (2) and the output shaft with 1-DOF rotational dynamics (3). These equations are coupled through the compliant mechanical transmission, modelled with lumped spring-damper parameters, assuming low structural damping.

$$\begin{gathered} V_d=R_{ph}{I}_d+L_{ph}{\dot{I}}_d-n_d{L}_{ph}{\dot{\theta }}_m{I}_q \\ {V}_q=R_{ph}{I}_q+L_{ph}{\dot{I}}_q+n_d{L}_{ph}{\dot{\theta }}_m{I}_d+K_e{\dot{\theta }}_m \end{gathered}$$

(1)

$$\begin{gathered} J_m{\ddot{\theta }}_m=-K_g/{\tau }_g\left({\theta }_m/{\tau }_g-{\theta }_o\right)-C_g/{\tau }_g\left({\dot{\theta }}_m/{\tau }_g-{\dot{\theta }}_o\right)-T_f+k_t{I}_q \\ {T}_f=B_m{\dot{\theta }}_m+T_c\tanh \left({\dot{\theta }}_m/{\omega }_c\right) \end{gathered}$$

(2)

$$J_o{\ddot{\theta }}_o=-K_g\left({\theta }_o-{\theta }_m/{\tau }_g\right)-C_g\left({\dot{\theta }}_o-{\dot{\theta }}_m/{\tau }_g\right)+T_{ext}$$

(3)

V_d, V_q and I_d, I_q are respectively the direct and quadrature voltages and currents. L_ph and R_ph are the motor phase inductance and resistance, n_d is the number of pole pairs, and K_e is the back-electromotive force constant of the PMSM. θ_m and θ_o are respectively the angular positions of the motor and the output shaft, and J_m and J_o are the associated inertias. K_g and C_g represent the torsional stiffness and damping of the differential ball-screw drivetrain with gear ratio τ_g. The motor torque constant k_t is equal to K_e for S-PMSM. B_m is the viscous friction coefficient, while T_c and ω_c are the parameters of the Coulomb-tanh friction model illustrated in Figure 3. Finally, T_ext is the external torque acting on the output shaft, accounting for load disturbances and aerodynamic actions.

The proposed model represents a balance between accuracy and complexity, consistent with the current objective of comparing control strategies rather than developing a high-fidelity simulation framework. More advanced EMA models usually include three-phase current dynamics, with Space Vector Pulse Width Modulation (SVPWM) and field-oriented control (FOC) [35], to represent real-time switching behaviour accurately. PMSM iron losses, detailed friction models and sensor noise were not included, since experimentally identified data were not available; however, these effects would impact both controllers equally. The parameters used in the electro-mechanical section of the simulation model are listed in Table 1.

2.2.2. Digital Electronics Section

The digital control section of the model includes two main main parts.

• Control laws: the EMA output position command is first rate-limited at the maximum output shaft speed ω_o^max, then compared to the resolver feedback. The position error is fed to the position regulator which generates a motor speed request. This is compared to the value derived from the motor resolver position feedback and fed to the motor speed regulator which outputs a quadrature current demand. The direct current instead has a fixed null reference to achieve Maximum Torque Per Ampere (MTPA) control for S-PMSM [36], assuming L_d = L_q = L_ph. On the baseline design, position, speed and current regulators have a Proportional-Integral structure with saturation and anti-windup scheme. In the MPC-based control, only the position regulator is replaced with an MPC, maintaining the other control laws unchanged. Additionally, current regulator integrates the decoupling technique. By reformulating eq. (1), the relationships in (4) can be established, linking the voltages V_d^reg, V_q^reg, calculated by the current regulator, to their post-decoupling counterparts V_d^dec, V_q^dec. This approach ensures that the current dynamics along the d-q axes remain unaffected by motor speed variations, leading to independent behaviour of the d and q currents. The equation (4) assumes that discrepancies between actual and measured values of currents and speed remain minimal, alongside with low variations in electrical constants, thus ensuring an effective decoupling. Finally, the voltages are dynamically saturated according to the constraint (5).

$$\begin{gathered} V_d^{dec}=V_d^{reg}+n_d{L}_{ph}{\dot{\theta }}_m{I}_q=R_{ph}{I}_d+L_{ph}{\dot{I}}_d, \\ {V}_q^{dec}=V_q^{reg}-n_d{L}_{ph}{\dot{\theta }}_m{I}_d-K_e{\dot{\theta }}_m=R_{ph}{I}_q+L_{ph}{\dot{I}}_q \end{gathered}$$

(4)

$$\sqrt{{\left(V_q^{dec}\right)}^2+{\left(V_d^{dec}\right)}^2}=V_{DC}^{max}/\sqrt{3}$$

(5)

• Sensors processing: the rotor and output shaft positions are measured through resolver angular sensors, while d-q currents are acquired with Hall-effect current sensors. On the real EMA, each motor phase current is measured, then d-q currents are calculated in the digital control section via Clarke-Park transformation (6).

$$\begin{gathered} \left[ \begin{array}{c}{I}_d\\ I_q\\ I_z\end{array} \right]={\mathbb{T}}_{PC}\left({\theta }_e\right)\left[ \begin{array}{c}{I}_a\\ I_b\\ I_c\end{array} \right], {\theta }_e\left(t\right)=n_d{\theta }_m\left(t\right) \\ {\mathbb{T}}_{PC}\left({\theta }_e\right)=\sqrt{{\displaystyle \frac{2}{3}}}·\left[ \begin{array}{ccc}\cos \left({\theta }_e\right)& \cos \left({\theta }_e-2\pi /3\right)& \cos \left({\theta }_e+2\pi /3\right)\\ -\sin \left({\theta }_e\right)& -\sin \left({\theta }_e-2\pi /3\right)& -\sin \left({\theta }_e+2\pi /3\right)\\ \sqrt{2}/2& \sqrt{2}/2& \sqrt{2}/2\end{array} \right] \end{gathered}$$

(6)

Sensor dynamics are modelled via continuous-time first order low-pass filters, whose characteristics are reported in Table 2. Sensor feedback are then digitized with sample and hold functions at 10 kHz rate.

2.3. Conventional control Laws Optimization

Conventional Proportional-Integral controllers are commonly used for their simplicity and capability to achieve zero steady-state tracking error. The general continuous-time transfer function is given by equation (7).

$$H_{PI}\left(\mathrm{s}\right)=K_P+{\displaystyle \frac{K_I}{s}}={\displaystyle \frac{K_P\left(s+K_I/K_P\right)}{s}}=K_P{\displaystyle \frac{\left(s+z_I\right)}{s}}:=K_P·{\widehat{H}}_{PI}\left(s\right)$$

(7)

The two parameters (for each control loop) to be selected are:

• the proportional gain: K_P;
• the zero of the integrator: z_I.

These parameters are tuned to meet the stability and closed-loop performance requirements reported in Table 3. It is worth noting that the bandwidths of the three control loops are designed to be separated by one decade, in order to obtain disjointed dynamics.

PI controllers are designed through an optimization procedure based on a reduced-order Linear Time-Invariant model of the EMA, derived from equations (1)–(4) by imposing zero Coulomb friction, since its linearization would generate unrealistically high values for ω > ω_c, and by neglecting the d-current dynamics (i.e. using a mono-phase equivalent model of the electrical machine (8)). Figure 4 illustrates the key degrees of freedom of this reduced electro-mechanical model.

$$V_q=R_{ph}{I}_q+L_{ph}{\dot{I}}_q+K_e{\dot{\theta }}_m$$

(8)

This approach enables the systematic definition and minimization of performance-oriented cost functions, with the goal of identifying the optimal controller parameters that guarantee accurate position tracking and effective load disturbance rejection with minimal energy consumption. This also provides a consistent basis for comparison with the MPC-based control strategy. The controller design is framed as a multi-objective optimization problem, where each cost function J_i​ quantifies a specific performance metric of the EMA under predefined test conditions. A total of six cost functions are defined based on the EMA response in the following two scenarios:

• Position tracking test: the commanded output position is a step of θ_os = 1 deg, with no external load applied. The cost function measures the squared tracking error between the output position and the reference command (9). Similarly, $J_{\mathsf{\Omega }}^{{\theta }_{os}}$ is linked to the deviation of the motor speed from the regulator reference (10). The cost function $J_P^{{\theta }_{os}}$ quantifies the normalized electric power absorption: since it includes the open-loop input voltage V_q, this metric reflects the required effort to achieve the desired tracking performance.

  $$J_{\mathsf{\Theta }}^{{\theta }_{os}}={{\displaystyle \int }}_0^t{\left(1-{\theta }_{o_{fb}}^{{\theta }_{os}}\left(\tau \right)\right)}^{2}d\tau$$

  (9)

  $$J_{\mathsf{\Omega }}^{{\theta }_{os}}={{\displaystyle \int }}_0^t{\left[\left({\omega }_{m_{dem}}^{{\theta }_{os}}\left(\tau \right)-{\omega }_{m_{fb}}^{{\theta }_{os}}\left(\tau \right)\right)/{\omega }_m^{max}\right]}^{2}d\tau$$

  (10)

  $$J_P^{{\theta }_{o_s}}={{\displaystyle \int }}_0^t\left|I_q^{{\theta }_{o_s}}\left(\tau \right)\cdot V_q^{{\theta }_{o_s}}\left(\tau \right)\right|/P_{el}^{max}d\tau$$

  (11)
• Load disturbance rejection test: the position command is set to the equilibrium condition (i.e. θ_o^cmd = 0 deg) while an external torque step of T_ext = 1 Nm is applied. The three corresponding cost functions are analogously defined (12)-(14).

  $$J_{\mathsf{\Theta }}^{T_{ext}}={{\displaystyle \int }}_0^t{\left({\theta }_o^{T_{ext}}\left(\tau \right)\right)}^{2}d\tau$$

  (12)

  $$J_{\mathsf{\Omega }}^{T_{ext}}={{\displaystyle \int }}_0^t{\left[\left({\omega }_{m_{dem}^{T_{ext}}}\left(\tau \right)-{\omega }_{m_{fb}^{T_{ext}}}\left(\tau \right)\right)/{\omega }_m^{max}\right]}^{2}d\tau$$

  (13)

  $$J_P^{T_{ext}}={{\displaystyle \int }}_0^t\left|I_q^{T_{ext}}\left(\tau \right)\cdot V_q^{T_{ext}}\left(\tau \right)\right|/P_{el}^{max}d\tau$$

  (14)

Finally, the global cost function is defined as a weighted linear combination of the individual terms J_i​, normalized by their respective minimum values to ensure comparability among different metrics (15). The minimum of equation (15) yields the optimal placement of the PI controller zeros.

$$J_w={\displaystyle \sum }_{i=1}^6{w}_i\cdot {\displaystyle \frac{J_i}{\min J_i}}-1 , 0<w_i<1, {\displaystyle \sum }_i{w}_i=1$$

(15)

To facilitate further results interpretation and to clarify the adopted weight assignment convention, it is useful to group each pair of cost functions associated to a specific performance metric: position tracking (16), speed tracking (17) and power consumption (18).

$$J_{\mathsf{\Theta }}=w_{\mathsf{\Theta }}^{{\theta }_{os}}\cdot {\displaystyle \frac{J_{\mathsf{\Theta }}^{{\theta }_{os}}}{\min \left(J_{\mathsf{\Theta }}^{{\theta }_{os}}\right)}}+w_{\mathsf{\Theta }}^{T_{ext}}\cdot {\displaystyle \frac{J_{\mathsf{\Theta }}^{T_{ext}}}{\min \left(J_{\mathsf{\Theta }}^{T_{ext}}\right)}}$$

(16)

$$J_{\mathsf{\Omega }}=w_{\mathsf{\Omega }}^{{\theta }_{os}}\cdot {\displaystyle \frac{J_{\mathsf{\Omega }}^{{\theta }_{os}}}{\min \left(J_{\mathsf{\Omega }}^{{\theta }_{os}}\right)}}+w_{\mathsf{\Omega }}^{T_{ext}}\cdot {\displaystyle \frac{J_{\mathsf{\Omega }}^{T_{ext}}}{\min \left(J_{\mathsf{\Omega }}^{T_{ext}}\right)}}$$

(17)

$$J_P=w_P^{{\theta }_{os}}\cdot {\displaystyle \frac{J_P^{{\theta }_{os}}}{\min \left(J_P^{{\theta }_{os}}\right)}}+w_P^{T_{ext}}\cdot {\displaystyle \frac{J_P^{T_{ext}}}{\min \left(J_P^{T_{ext}}\right)}}$$

(18)

The controller optimization is carried out in the MATLAB environment, based on a state-space representation of the EMA LTI model. The procedure consists of the following main steps:

• Zeros placement variation: a common variation vector for zero placement is defined and applied equally across all control loops. The zeros location is varied as fractions of the respective closed-loop design bandwidths, within the range [0.1–0.5], to ensure a minimal impact of the zero on the bandwidth itself. Adopting a Monte Carlo-like optimization algorithm, a large number of zeros placement combinations are generated for the three nested control loops.
• Feed-forward design and stability verification: for each configuration, the feed-forward transfer function is constructed by cascading ${\widehat{H}}_{PI}\left(s\right)$ with the open-loop transfer function from its output to the EMA system output in the corresponding control loop. The proportional gain is then automatically calculated by enforcing that the cut-off frequency of the feed-forward transfer function matches the predefined closed-loop bandwidth. Stability margins are checked, and if they meet the required specifications, the state-space LTI closed-loop model of the EMA system is constructed.
• Time-domain simulation and cost functions evaluation: the position closed-loop EMA is simulated under position tracking and disturbance rejection tests; consequently, the values of the six cost functions are evaluated for each zeros configuration.
• Normalization of cost functions: for each cost function, the minimum value across all combinations is identified to allow its normalization. This ensures consistent comparison across metrics with different physical units and scales.
• Weight definition for multi-objectives cost function: given the focus on primary flight control applications, the optimization strategy prioritizes position tracking performance over secondary objectives. Accordingly, the weights associated with motor speed tracking and power consumption are fixed to small values (19); conversely, the weights associated with position tracking cost functions are allowed to vary within a discrete set of values (20) to be treated as tuneable parameters within the optimization process.

  $$w_{\mathsf{\Omega }}^{{\theta }_{os}}=w_{\mathsf{\Omega }}^{T_{ext}}=w_P^{{\theta }_{os}}=w_P^{T_{ext}}=0.1$$

  (19)

  $$w_{\mathsf{\Theta }}^{{\theta }_{o_s}}\in \left\{0.3, 0.4, 0.5\right\}, w_{\mathsf{\Theta }}^{T_{d_s}}\in \left\{0.1, 0.2, 0.3\right\}$$

  (20)

Only combinations satisfying the constraint $\sum w_i=1$ are valid. This approach effectively incorporates weights tuning into the overall optimization framework, thereby expanding the dimensionality of the parameters space and enabling a balanced trade-off among the multiple performance objectives.

6.

Computation of the global cost function: for each valid combination of zero placements and cost function weights, the global cost function J_w is evaluated according to (15).

7.

Optimal design identification: the optimal controllers configuration is identified as the one that minimizes the global cost function​, ensuring the best trade-off among the considered performance metrics.

2.4. Model Predictive Controller development

An MPC for EMA position regulation was developed based on the standard theory of discrete-time, LTI, state-space formulation for single-input-single-output model-based predictive control [37]. This control approach represents a novel strategy in aerospace EMA applications. Its objective is to optimize the system performance by predicting future behaviour and computing control actions accordingly. This technique relies on a dynamic model of the plant and solves an optimization problem at each time step to minimize a cost function while satisfying system constraints. The implementation of an MPC algorithm requires the following elements:

• Information on the current and previous states of the system;
• A mathematical model of the plant to predict its future evolution;
• An optimization criterion (typically a quadratic cost function).

The MPC strategy is based on three core concepts, displayed in Figure 5:

• Prediction Horizon (N_p) defines the finite time span (in terms of discrete time steps k) over which the future system output $\widehat{y}\left(k\right)$ is predicted, using the plant model and current state information. A longer prediction horizon allows the controller to capture a more comprehensive understanding of the system evolution, though increasing the complexity of the algorithm.
• Control Horizon (N_c) denotes the finite time window over which the future control inputs u(k) are optimized, i.e. for k=k_i,…,k_i+N_c-1, where k_i is the initial time. N_c is constrained by computational costs and the need to balance predictive accuracy with real-time implementation. Thus, N_c is usually smaller than N_p ​, and for k>k_i+N_c-1 the control input is assumed to remain constant. The selection of an appropriate control horizon is crucial as it directly influences the controller ability to anticipate and regulate the system behaviour effectively.
• Receding Horizon control: a dynamic optimization framework enables the controller to adapt its control strategy in real-time by recalculating control inputs at each time step, based on the most recent state information, to enhance robustness against model inaccuracies and external perturbations and to maintain system stability. Mathematically, at time k, the MPC computes an optimal sequence of control actions over the full control horizon, however, only the first input u(k) is applied to the plant. At the next time step, the system state is updated, and a new optimization is performed.

2.4.1. Discrete-time LTI MPC formulation

The design of a Model Predictive Controller is based on a mathematical model of the plant, with its algorithm complexity increasing as the model becomes more detailed or nonlinear. A basic MPC formulation relies on a discrete-time, LTI, SISO, state-space model (the extension to multiple-input multiple-output (MIMO) systems is straightforward), which embedded dynamic equations are formulated as (21), (22), where k is the actual time sample, x_m(k) is the vector of states at time k with length n₁, u(k) and y(k) are respectively the input and output at the current time (both scalar for SISO), A_m(n₁ x n₁), B_m(n₁ x 1) and C_m(1 x n₁) are the plant state-space representation matrices with constant parameters. Due to the receding horizon technique, the D_m matrix is always null, since the input cannot affect directly the output at the same time step. If the dynamic model is nonlinear, the matrices can be obtained via linearization around an equilibrium point.

$$x_m\left(k+1\right)=A_m{x}_m\left(k\right)+B_{m}u\left(k\right)$$

(21)

$$y\left(k\right)=C_m{x}_m\left(k\right)$$

(22)

Expressing the variation of the input and states within two consecutive time samples (23), it is possible to rearrange the dynamic equations (21),(22) in incremental terms (24),(25).

$$\begin{gathered} \mathsf{\Delta }{x}_m\left(k+1\right):=x_m\left(k+1\right)-x_m\left(k\right) \\ \mathsf{\Delta }u\left(k\right):=u\left(k\right)-u\left(k-1\right) \end{gathered}$$

(23)

$$\mathsf{\Delta }{x}_m\left(k+1\right)=A_m\mathsf{\Delta }{x}_m\left(k\right)+B_m\mathsf{\Delta }u\left(k\right)$$

(24)

$$y\left(k+1\right)=y\left(k\right)+C_m{A}_m\mathsf{\Delta }{x}_m\left(k\right)+C_m{B}_m\mathsf{\Delta }u\left(k\right)$$

(25)

Defining the so-called augmented states vector (26), the augmented state-space LTI model of the plant is obtained (27)-(29) used to design the MPC controller and for the calculation of future states and outputs.

$$x\left(k\right):={\left[\mathsf{\Delta }{x}_m{\left(k\right)}^T y\left(k\right)\right]}^T$$

(26)

$$\begin{gathered} x\left(k+1\right)=A_{a}x\left(k\right)+B_a\Delta u\left(k\right) \\ y\left(k\right)=C_{a}x\left(k\right) \end{gathered}$$

(27)

$$\begin{gathered} \stackrel{\underset{⏞}{x\left(k+1\right)}}{\left[\begin{array}{c}\mathsf{\Delta }{x}_m\left(k+1\right)\\ y\left(k+1\right)\end{array}\right]}=\stackrel{\underset{⏞}{A_a}}{\left[\begin{array}{cc}{A}_m& 0_m^T\\ C_m{A}_m& 1\end{array}\right]}\stackrel{\underset{⏞}{x\left(k\right)}}{\left[\begin{array}{c}\mathsf{\Delta }{x}_m\left(k\right)\\ y\left(k\right)\end{array}\right]}+\stackrel{\underset{⏞}{B_a}}{\left[\begin{array}{c}{B}_m\\ C_m{B}_m\end{array}\right]}\mathsf{\Delta }u\left(k\right) \\ y\left(k\right)=\underset{C_a}{\underbrace{\left[0_m 1\right]}}\left[\begin{array}{c}\mathsf{\Delta }{x}_m\left(k\right)\\ y\left(k\right)\end{array}\right] \end{gathered}$$

(28)

$$\begin{array}{l}{A}_a=\left[\begin{array}{cc}{A}_m& 0_m^T\\ C_m{A}_m& 1\end{array}\right], \dim \left(A_a\right)=\left(n_1+1\right)\times \left(n_1+1\right)\\ B_a=\left[\begin{array}{c}{B}_m\\ C_m{B}_m\end{array}\right], \dim \left(B_a\right)=\left(n_1+1\right)\times \left(n_1+1\right)\\ C_a=\left[0_m 1\right], \dim \left(C_a\right)=1\times \left(n_1+1\right)\\ 0_m:=\stackrel{n_1}{\overbrace{\left[0 0 \dots 0\right]}}\end{array}$$

(29)

2.4.2. Prediction Within One Optimization Window

At each discrete time step k_i > 0, the MPC algorithm solves an optimization problem to compute the N_c​ future increments of the control variable Δu(k_i),…,Δu(k_i + N_c - 1). Based on this control sequence, using the augmented model, the algorithm then predicts sequentially the future states (30), and the output (31) of the system over the prediction horizon N_p (i.e. the optimization window). The notation ​x(k | k_i) denotes the predicted state for future time k = k_i + 1, …, k_i + N_p carried out at time k_i. By defining two vectors containing all N_p values of input (32) and output (33), the dynamic equation used for output prediction can be rewritten in a compact, state-space-like form (34). Remarkably, the F and Φ matrices depend only on the augmented model matrices, defined as constant, enabling a significant simplification of the MPC code generation.

$$\begin{array}{ll}x\left(k_i+1 | k_i\right)& =A_{a}x\left(k_i\right)+B_a\mathsf{\Delta }u\left(k_i\right)\\ x\left(k_i+2 | k_i\right)& =A_{a}x\left(k_i+1 | k_i\right)+B_a\mathsf{\Delta }u\left(k_i+1\right)=\\ \dots & ={\left(A_a\right)}^{2}x\left(k_i\right)+A_a{B}_a\mathsf{\Delta }u\left(k_i\right)+B_a\mathsf{\Delta }u\left(k_i+1\right)\\ & ⋮\\ x\left(k_i+N_p | k_i\right)& ={\left(A_a\right)}^{N_p}x\left(k_i\right)+A_a^{N_p-1}B\Delta u\left(k_i\right)+A_a^{N_p-2}B\Delta u\left(k_i+1\right)+\dots \\ & +A_a^{N_p-N_c}B\Delta u\left(k_i+N_c-1\right)\end{array}$$

(30)

$$\begin{array}{ll}y\left(k_i+1 | k_i\right)& =C_a{A}_{a}x\left(k_i\right)+C_a{B}_a\mathsf{\Delta }u\left(k_i\right)\\ y\left(k_i+2 | k_i\right)& =C_a{\left(A_a\right)}^{2}x\left(k_i\right)+C_a{A}_a{B}_a\Delta u\left(k_i\right)+C_a{B}_a\mathsf{\Delta }u\left(k_i+1\right)\\ & ⋮\\ y\left(k_i+N_p | k_i\right)& =C_a{\left(A_a\right)}^{N_p}x\left(k_i\right)+C_a{A}_a^{N_p-1}B\mathsf{\Delta }u\left(k_i\right)+C_a{A}_a^{N_p-2}{B}_a\mathsf{\Delta }u\left(k_i+1\right)+\dots \\ & +C_a{A}_a^{N_p-N_c}{B}_a\mathsf{\Delta }u\left(k_i+N_c-1\right)\end{array}$$

(31)

$$\mathsf{\Delta }U={\left[\mathsf{\Delta }u\left(k_i\right) \mathsf{\Delta }u\left(k_i+1\right)\dots \mathsf{\Delta }u\left(k_i+N_c-1\right)\right]}^T$$

(32)

$$Y={\left[y\left(k_i+1|k_i\right)y\left(k_i+2|k_i\right)\dots y\left(k_i+N_p|k_i\right)\right]}^T$$

(33)

$$\begin{gathered} Y=F x\left(k_i\right)+\mathsf{\Phi }\mathsf{\Delta }U \\ F=\left[\begin{array}{c}{C}_a{A}_a\\ C_a{A}_a^2\\ \begin{array}{c}:\\ C_a{A}_a^{N_p}\end{array}\end{array}\right] \\ \mathsf{\Phi }=\left[\begin{array}{ccccc}{C}_a{B}_a& 0& 0& \dots & 0\\ C_a{A}_a{B}_a& C_a{B}_a& 0& \dots & 0\\ C_a{A}_a^2{B}_a& C_a{A}_a{B}_a& C_a{B}_a& \dots & 0\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ C_a{A}_a^{N_p-1}& C_a{A}_a^{N_p-2}{B}_a& C_a{A}_a^{N_p-3}{B}_a& \dots & C_a{A}_a^{N_p-N_c}{B}_a\end{array}\right] \end{gathered}$$

(34)

2.4.3. Optimization Criterion

The control objective of the MPC is to minimize the deviation between the predicted output and a reference signal r(k_i), which is assumed to remain constant over the prediction horizon (35). This goal is achieved by computing the optimal sequence of control increments ΔU that minimizes a convex quadratic cost function (36). The first part aims to minimize the tracking error e(k_i) = r(k_i) – y(k_i), while the second term penalizes excessive control effort by weighting the magnitude of ΔU. This trade-off is regulated by a tuneable parameter r_w ≥ 0, embedded in the diagonal input-weighting matrix ${\overline{R}}_w=I_{N_c\times N_c}\cdot r_w$, where $I_{N_c\times N_c}$ is the identity matrix. The input-weight parameter r_w governs the balance between tracking performance and control effort. Setting r_w = 0 results in a purely performance-oriented controller, neglecting the cost of control variation. Conversely, larger r_w​ values yield smoother control inputs by limiting abrupt variations, thus producing a more conservative response.

$$R_r^T:={\overline{R}}_r\cdot r\left(k_i\right), {\overline{R}}_r=\stackrel{N_p}{\overbrace{\left[1, 1,\dots ,1\right]}}$$

(35)

$$J=\underset{\begin{array}{c}Minimization of\\ tracking error\end{array}}{\underbrace{{\left(R_r-Y\right)}^T\left(R_r-Y\right)}}+\underset{\begin{array}{c}Minimization of\\ input "effort"\end{array}}{\underbrace{ \mathsf{\Delta }{U}^T {\overline{R}}_w \mathsf{\Delta }U }}$$

(36)

By substituting equation (34) into (36), the expanded expression shown in (37) is obtained.

$$J={\left(R_r-F x\left(k_i\right)\right)}^T\left(R_r-F x\left(k_i\right)\right)-2\mathsf{\Delta }{U}^T{\mathsf{\Phi }}^T\left(R_r-F x\left(k_i\right)\right)+\mathsf{\Delta }{U}^T\left({\mathsf{\Phi }}^T\mathsf{\Phi }+{\overline{R}}_w\right)\mathsf{\Delta }U$$

(37)

By setting the gradient of the cost function with respect to ΔU to zero, and considering the problem’s convex and quadratic structure, a unique global minimum is guaranteed, which corresponds to the optimal control trajectory (38). The term ${\left({\mathsf{\Phi }}^T\mathsf{\Phi }+{\overline{R}}_w\right)}^{-1}$ in (38) is commonly referred to as the Hessian matrix, whose existence is guaranteed by the problem construction. Since both Φ and ${\overline{R}}_w$ are constant matrices in this formulation, the resulting Hessian matrix is also constant.

$$\Delta U_{opt}={\left[{\mathsf{\Phi }}^T\mathsf{\Phi }+\underset{\begin{array}{c}\Delta U\\ weight\end{array}}{\underbrace{ {\overline{R}}_w}}\right]}^{-1}{\mathsf{\Phi }}^T\left[\underset{Set-point}{\underbrace{ {\overline{R}}_s r\left(k_i\right)}}-\underset{\begin{array}{c} Current\\ state\end{array}}{F\underbrace{x\left(k_i\right)}}\right]$$

(38)

2.4.4. Model Predictive Position Controller Design

The first design decision involves choosing how to integrate the Model Predictive Controller: since classical EMA architectures consist of three nested loops, one option would be to replace all of them with a single MPC regulator. However, this strategy introduces significant challenges: internal saturation limits (e.g., maximum motor speed and current), typically enforced by external saturation blocks, must be explicitly addresses as hard constraints on system states, inputs and outputs. This approach would require solving a constrained optimization problem with linear inequalities constraints. Although MPC algorithms can address this task using Quadratic Programming methods such as Hildreth’s algorithm [38], the associated computational complexity and potential constraints conflicts make this solution unsuitable for safety-critical applications. Preliminary control laws design and further optimization of PI controllers highlighted a strong dependence of the EMA closed-loop performance on the tuning of the position regulator, while motor speed and current regulators have a secondary influence. Based on these evidence, the MPC was employed as the position regulator, in cascade with the optimized PI controllers for speed and current loops. Aerospace applications impose a trade-off between embedded model accuracy and necessary computational resources: certification requirements in safety-critical systems limit the use of complex or unverified algorithms. Therefore, the embedded model must be kept as simple as possible while still capturing the essential EMA dynamics. To fulfil these requirements, the MPC was developed without internal constraints and based on a simplified LTI model. Specifically, since the MPC operates as a position regulator, its embedded model corresponds to the speed closed-loop EMA shown in Figure 6, represented by the transfer function (39) and discretized at a 10 kHz sampling rate. Referring to Figure 6, the open-loop model of the EMA corresponds to the five-states representation used also for PIs optimization (Figure 4). The inclusion of current and speed loops integrators adds two additional states. An additional state results from the discrete-time derivative and filtering used in speed reconstruction, yielding a total of eight states (40). These additional states are not directly measurable but can be computed from the PI regulators' internal variables and sensor processing at the previous time step. The output speed ω_o is estimated using the same approach applied to calculate the motor speed from the position feedback.

$$EM{A}_{MPC}\left(s\right)={\displaystyle \frac{{\theta }_o}{{{\omega }_m}_{dem}}}\left(s\right)$$

(39)

$$x_{EM{A}_{MPC}}={\left[I_q, {\theta }_m, {\omega }_m, {\theta }_o, {\omega }_{\mathrm{o}}, \int I_q, \int {\omega }_m, {\omega }_m^{d/dt}\right]}^T$$

(40)

In the absence of a formal procedure to determine an optimal configuration for MPC parameters (N_c, N_p and r_w) a-priori, the design was carried out through an iterative refinement. A series of closed-loop simulations were conducted, progressively adjusting the tuning parameters based on the observed performance, until satisfactory tracking accuracy and closed-loop stability were achieved with minimal control effort. This iterative process led to some key considerations:

• a short control horizon (N_c​ = 5) proved sufficient, as longer values did not significantly improve performance;
• a long prediction horizon (N_p​ ​= 500) was necessary to properly capture the system dynamics and generate effective control actions;
• Given the incremental nature of the MPC law, the input-weighting factor should remain small (r_w = 10^-3) to prevent over-smooth inputs, thereby preserving actuator responsiveness and ensuring adequate control authority.

In the nonlinear Simulink model, the MPC is implemented via a custom MATLAB function that computes the optimal control trajectory as expressed in (38). The necessary matrices (41), being constant in this formulation, are precomputed and loaded into the code. By opting for a significantly reduced control horizon, these matrices are of small dimension. To compute ΔU_opt, the MPC block requires the reference input r(k) = θ_o^cmd(k) and the augmented states vector (26); therefore, state values at both the current and previous time steps are needed. Finally, the MPC applies only the first element of the control increment vector to compute the optimal speed reference, as defined in equation (42). Finally, the MPC output is externally saturated, to avoid a more complex constrained optimization procedure, and it is sent to the motor speed regulator in cascade.

$$\begin{gathered} Inv={\left({\mathsf{\Phi }}^T\mathsf{\Phi }+{\overline{R}}_w\right)}^{-1}, \dim \left(Inv\right)=N_c\times N_c \\ Phi{T}_{Rr}={\mathsf{\Phi }}^T\cdot {\overline{R}}_r, \dim \left(Phi{T}_{Rr}\right)=N_c\times 1 \\ PhiT\_F={\mathsf{\Phi }}^T\cdot F, \dim \left(PhiT\_F\right)=N_c\times \left(n_1+1\right) \end{gathered}$$

(41)

$${{\omega }_m^{opt}|}_k={{\omega }_m^{opt}|}_{k-1}+{\mathsf{\Delta }{U}_{opt}|}_k\left(1\right)$$

(42)

2.5. Performance Comparison and Robustness Analysis

To substantiate the differences between the two control strategies, the following simulations in time and frequency domain have been performed.

• The EMA unloaded position tracking response is assessed in the frequency domain with harmonic commands of 1 deg amplitude.
• In the loaded position tracking tests, aerodynamic actions on the control surface are modelled as an equivalent torsional spring with constant stiffness K_aero, acting on the EMA output shaft (43).

$$T_{aero}\left(t\right)=-K_{aero}\cdot {\theta }_o\left(t\right)$$

(43)

The tests are conducted using sinusoidal commands with amplitudes of 0.1 deg, 1 deg and 10 deg.

3.

The EMA dynamic compliance is derived imposing harmonic and step-like loads of 1 Nm amplitude.

To quantify the energetic benefit of MPC, the percentage energy reduction is defined as in (44), where the energy terms (E_MPC, E_PI) are evaluated according to the respective power integral (45).

$$E_{cons.}\%={\displaystyle \frac{E_{MPC}-E_{PI}}{E_{PI}}}\cdot 100$$

(44)

$$E=\underset{o}{\overset{t}{{\displaystyle \int }}}{P}_{el} d\tau =\underset{o}{\overset{t}{{\displaystyle \int }}}{\displaystyle \frac{3}{2}}\left(V_q^{ref}{I}_q^{fb}+V_d^{ref}{I}_d^{fb}\right) d\tau$$

(45)

To assess the robustness of the control strategies, a reference test sequence is designed. The actuator is first commanded to track and hold a 1 deg position step command under aerodynamic spring-load. Subsequently, an external load step of 1 Nm is applied on the output shaft to evaluate the ability to reject gust disturbances. The robustness analysis is conducted in two phases. First, a comprehensive single-parameter sensitivity analysis is carried out by individually varying several parameters of the EMA nonlinear model. Due to the high inherent robustness of both controllers, only the most relevant results are reported, related to the variation of the following parameters:

• PMSM torque constant k_t ;
• friction torque multiplier K_f ;
• equivalent aerodynamic spring load stiffness K_aero.

Figure 7 displays the ten selected points for each parameter along with their mean values (used as baseline) and the standard deviation.

As a second step, focusing on MPC, a multi-parameter robustness analysis is performed to account for potential nonlinear interactions among uncertainties. A Monte Carlo-like approach is employed to generate 10³ combinations of the parameters above. The MPC robustness is then evaluated in comparison to the baseline performance. Key statistics metrics include the deviation of the position tracking error (46) and the deviation in the electrical power consumption (47), providing a quantitative assessment of tracking accuracy and energy efficiency under uncertainty. For each Monte Carlo run and for both performance metrics, the four statistical moments (the mean, the standard deviation, the skewness and the kurtosis) are computed over the simulation timespan N_sim. This analysis intended to support a statistical robustness evaluation by estimating the likely distribution of performance metrics across a representative range of uncertain conditions.

$${\delta }_{{e_{\theta }}_o}={{\overline{e}}_{\theta }}_o\left(k_t,K_{aero},K_f\right)-{\overline{e}}_{{\theta }_o^{baseline}} , {{\overline{e}}_{\theta }}_o={\displaystyle \frac{1}{N_{sim}}}{\displaystyle \sum }_{k=1}^{N_{sim}}\left|{\theta }_o^{cmd}\left(k\right)-{\theta }_o^{fb}\left(k\right)\right|$$

(46)

$${\delta }_P={\overline{P}}_{el}\left(k_t,K_{aero},K_f\right)-{\overline{P}}_{el}^{baseline}, {\overline{P}}_{el}={\displaystyle \frac{1}{N_{sim}}}{\displaystyle \sum }_{k=1}^{N_{sim}}{P}_{el}\left(k\right)$$

(47)

3. Results

3.1. Proportional-Integral Controllers Optimization

Figure 8, Figure 9 and Figure 10 provide a concise overview of the main trends observed in the PI cost functions. Figure 8 shows the individual cost terms plotted against the normalised locations of the zeros placement: the x axis reports $z_{I_q}/{\mathsf{\Omega }}_{I_q}^{BW}$, the y axis $z_{{\mathsf{\Omega }}_m}/{\mathsf{\Omega }}_{{\mathsf{\Omega }}_m}^{BW}$, and the different surfaces correspond to different values of $z_{{\mathsf{\Theta }}_o}/{\mathsf{\Omega }}_{{\mathsf{\Theta }}_o}^{BW}$. Results from both the position-tracking and disturbance-rejection tests suggested that, for each performance metric, the corresponding pair of cost functions exhibit similar trend, hence the definitions (16)-(18) illustrated in Figure 9. Conversely, position- and speed-related cost functions displayed opposite trends: combinations that minimizes $J_{\mathsf{\Theta }}$ inevitably increased $J_{\mathsf{\Omega }}$, emphasising the necessity of a multi-objectives optimization. The current loop zero had only a marginal influence on the power cost index of disturbance rejection tests. Finally, Figure 10 locates every controller design in the three-dimensional space ($J_{\mathsf{\Theta }}$, $J_{\mathsf{\Omega }}$, $J_P$). The plane (48) is tangent to the Pareto front of the multi-objective optimization exactly at a single point (highlighted), which therefore represents the global optimum and the recommended zero placement configuration for PIs design, reported numerically in Table 4.

$$J_{\mathsf{\Theta }}+J_{\mathsf{\Omega }}+J_P=\min \left(J_{\mathsf{\Theta }}+J_{\mathsf{\Omega }}+J_P\right)$$

(48)

3.2. PI and MPC-Based Control Performance Comparison

3.2.1. Unloaded Position Tracking Tests

The Bode diagram in Figure 11 compares the MPC- and PI-based systems frequency responses during unloaded position tracking tests. The acceptance mask for the frequency response is defined as follows.

• At frequencies lower than 0.2 Hz the magnitude shall remain within ±1 dB, while the phase shall be higher than -10 deg.
• For frequencies between 0.2 Hz and 2 Hz, a linear increase in the magnitude is allowed till ±3 dB at the design bandwidth, whereas the phase shall remain above -60 deg.
• Beyond 2 Hz the magnitude is expected to decay with a slope between –20 dB/dec and –40 dB/dec, while the phase must be greater than -180 deg.

While both architectures achieved the nominal position-loop bandwidth, the MPC-controlled EMA exhibited a flatter gain profile and reduced phase lag up to approximately 6 Hz, suggesting enhanced dynamic performance.

The computational effort required by both control strategies was evaluated through a one-second simulation with step-like position command, performed under identical hardware conditions. Execution times were comparable (1.77 s for the PI and 1.83 s for the MPC), indicating the feasibility of MPC real-time implementation without extra processing resources.

3.2.2. Loaded Position Tracking Tests

Figure 12 shows the loaded frequency response at increasing levels of command amplitude. At 0.1 deg, typical of aircraft stability augmentation functions, the PI controller suffered early gain attenuation and phase lag, while MPC maintained full sensitivity at small commands. At 1 deg, representative of typical autopilot commands, both controllers met the performance requirements. Under extreme manoeuvres (10 deg), saturation effects caused both to exceed the acceptance limits.

3.2.3. Load Disturbance Rejection Tests

The response of the EMA to a step-like gust disturbance is shown Figure 13. Despite a comparable peak in output position, the MPC-controlled actuator returned to equilibrium significantly faster, demonstrating superior damping. This improved rejection capability was further confirmed by the motor speed and current profiles: while PI induced strong oscillations, the MPC ensured smoother transients. In addition, energy consumption was reduced by 20%, highlighting MPC as a promising solution also for aircraft secondary control surfaces.

The Bode diagram of the EMA dynamic compliance, defined as the transfer function from external torque input to the angular displacement, is illustrated in Figure 14. Around the design bandwidth, the MPC achieved slightly greater load attenuation and a higher phase lead, enabling faster compensation of the disturbance. At low frequencies, PI outperformed MPC due to better phase anticipation. In the high frequency range, the responses of the two controllers converged, as the actuator resonance mode dominated the overall dynamics.

3.3. Robustness Analysis

The time-domain responses of MPC- and PI-controlled EMAs under single-parameter variation are illustrated in Figure 15, Figure 16 and Figure 17. Considering the variation of the motor torque constant (Figure 15), both control strategies exhibited strong robustness in the output position response. However, significant sensitivity is observed in the electrical quantities. At lower values of the torque constant k_t​, a higher quadrature current is required to generate the same motor torque, resulting in increased power peaks in both systems. Nevertheless, in this extreme scenario, the relative energy saving achieved by the MPC strategy became even more evident, highlighting its superior efficiency.

Uncertainty in the friction coefficients (Figure 16) introduced noticeable deviations not only in the electrical but also in the mechanical responses. Deviations in the output position became significant, particularly for the PI controller, highlighting its higher sensitivity to friction-induced drag. The MPC controller, instead, exhibited improved consistency in position tracking and faster recovery following the disturbance, while showing reduced energy consumption, especially under severe parameter deviations.

Unlike the previous cases, variations in K_aero​ influenced the external load model rather than the intrinsic characteristics of the actuator. The MPC controller was able to leverage this variation more effectively, progressively enhancing its energy efficiency from the unloaded condition up to the maximum aerodynamic load.

The second phase of the analysis focused on a detailed robustness characterization of the MPC with respect to combined variations of the parameters k_t​, K_f​ and K_aero​. Figure 18 shows the statistical distribution of the position error deviation across all Monte Carlo simulations, in terms of mean, standard deviation, skewness, and kurtosis. The distribution of the mean value, centered near zero with a slight positive bias, indicated that tracking performance is generally preserved. The standard deviation revealed a multi-modal profile, suggesting the presence of distinct sensitivity clusters due to the discrete influence of the K_aero​​. Skewness values were broadly distributed around zero, showing no dominant direction asymmetry. The high kurtosis observed across the dataset confirmed a leptokurtic distribution, characterized by a non-negligible occurrence of extreme outliers. These cases represented critical parameter combinations that lead to significant deviations in tracking performance.

To better understand the interaction between uncertain parameters and their impact on control performance, the statistical metrics of the position error deviation were mapped over the k_t​ - K_f​ plane and layered with K_aero​ values (Figure 19). The mean and standard deviation plots revealed a reduction in control effectiveness under conditions of low electro-mechanical coupling combined with high friction and external loads.

The power consumption deviation metrics are summarized in Figure 20. The distribution of the mean value was predominantly non-negative and slightly right-skewed, indicating that most uncertain configurations tend to increase average power consumption compared to the nominal scenario. Notably, a secondary mode with negative values emerged, suggesting that a subset of parameter combinations can also lead to net energy savings. This behaviour aligned with the bimodal trend observed in the skewness plot. The heavy-tailed standard deviation suggested nonlinear amplification of power fluctuations under uncertainty. The kurtosis appeared strongly leptokurtic, with frequent extreme outliers. This reflected the system susceptibility to large energetic deviations under specific uncertainty conditions, in contrast to the more constrained behaviour observed in the position error analysis.

Figure 21 is useful for identifying specific regions where strong power performance degradation occurred. Similarly to the position trend, the mean and standard deviation​ increased significantly for low k_t​ and high K_f​​, highlighting a critical region for EMA robustness. Skewness confirmed highly asymmetric and heavy-tailed behaviour in the same region. Interestingly, the combined influence of extreme values of K_aero​​​ and k_t​ appeared secondary in terms of mean and standard deviation but had a significant impact on kurtosis, suggesting that certain parameter combinations might have amplified rare yet extreme power deviations.

4. Discussion

This study provides a comprehensive assessment of the performance and robustness of conventional PI-based and advanced MPC strategies applied to an electromechanical actuator for primary flight control surfaces. The findings highlight distinct response patterns between conventional PI-based and model predictive control strategies, driven by their respective structural characteristics and interaction with system dynamics. To ensure a fair assessment, the baseline PI controller was first tuned through a multi-objective optimization procedure, balancing dynamic response and power consumption. Nevertheless, the MPC controller consistently outperformed the PI architecture in terms of both dynamic response and energy efficiency, across test scenarios representative of primary and secondary flight control applications. Notably, this was obtained without a significant increase in computational cost, as simulation times remained comparable.

The comparison of PI and MPC robustness under single-parameter variations highlighted a major degradation in PI performance under high friction, low torque constant, and stiff aerodynamic loads. In contrast, the MPC demonstrated greater overall stability, enabled by its predictive formulation, even relying on a non-adaptive embedded model.

Robustness characterization of the MPC design was further assessed via Monte Carlo analysis over a multidimensional uncertainty space. The results confirmed that control effectiveness is not uniformly distributed across the uncertainty space, but rather influenced by nonlinear interactions between electromechanical and external load conditions. It also enabled the identification of critical parameter regions and rare but impactful deviations on position tracking accuracy and power consumption, which are especially relevant for control system design targeting certification and fault resilience in aerospace applications.

Although the study covered a broad range of single and combined parameter variations, some limitations remain. Sensor noise, mechanical degradation and failure modes were not included in the current model and may further influence the robustness envelope. Simulation model refinements are planned to account for power electronics switching effects and provide a more accurate description of friction dynamics. While the current MPC formulation does not incorporate real-time linearization or uncertainty estimation, such features could enhance robustness, although their computational impact must be carefully evaluated. Experimental testing on physical EMA hardware is expected to complement the simulation-based findings and provide insight into real-world performance under operational constraints.

5. Conclusions

This study proposed a simplified and real-time implementable MPC strategy for electro-mechanical actuators in aerospace applications. The control architecture was validated through an extensive simulation campaign, demonstrating improved tracking performance, energy efficiency, and load disturbance rejection compared to an optimized PI benchmark. The robustness-oriented methodology confirmed the MPC reliability under both individual and combined parameters variation, supporting its suitability for both primary and potentially secondary flight control surfaces.