**Electro-Mechanical Actuators for Safety-Critical Aerospace Applications**

Editor

**Gianpietro Di Rito**

MDPI • Basel • Beijing • Wuhan • Barcelona • Belgrade • Manchester • Tokyo • Cluj • Tianjin

*Editor* Gianpietro Di Rito Dipartimento di Ingegneria Civile ed Industriale–Sede Aerospaziale Universita di Pisa ` Pisa Italy

*Editorial Office* MDPI St. Alban-Anlage 66 4052 Basel, Switzerland

This is a reprint of articles from the Special Issue published online in the open access journal *Aerospace* (ISSN 2226-4310) (available at: www.mdpi.com/journal/aerospace/special issues/ Electro mechanical actuators).

For citation purposes, cite each article independently as indicated on the article page online and as indicated below:

LastName, A.A.; LastName, B.B.; LastName, C.C. Article Title. *Journal Name* **Year**, *Volume Number*, Page Range.

**ISBN 978-3-0365-7933-7 (Hbk) ISBN 978-3-0365-7932-0 (PDF)**

© 2023 by the authors. Articles in this book are Open Access and distributed under the Creative Commons Attribution (CC BY) license, which allows users to download, copy and build upon published articles, as long as the author and publisher are properly credited, which ensures maximum dissemination and a wider impact of our publications.

The book as a whole is distributed by MDPI under the terms and conditions of the Creative Commons license CC BY-NC-ND.

## **Contents**


## **Chenfei She, Ming Zhang, Yibo Ge, Liming Tang, Haifeng Yin and Gang Peng**


## **Zhangbin Wu, Hongbai Bai, Guangming Xue and Zhiying Ren**


## **Jean-Charles Mar ´e**


## **About the Editor**

## **Gianpietro Di Rito**

Born in Massafra, Italy, October 2nd, 1974, he received his MSc Degree in Aerospace Engineering in 2001 and his PhD Degree in Aerospace Engineering in 2005, at the University of Pisa, Pisa, Italy. He was a Research Fellow at the Department of Aerospace Engineering of the University of Pisa from 2006 to 2011. From 2012 to 2021, he was an Assistant Professor of Aerospace Systems at the Department of Civil and Industrial Engineering of the University of Pisa, where in 2022 he became an Associate Professor of Aerospace Systems.

His main research interests are focused on the modelling, simulation, control, health-monitoring and testing of mechatronic systems for aerospace applications.

## **Preface to "Electro-Mechanical Actuators for Safety-Critical Aerospace Applications"**

Aircraft electrification is one of the most important and strategic initiatives currently supporting the innovation of the aviation industry. This manifests in the well-known more-electric aircraft concept (with the ultimate aim of achieving the all-electric long-term target), which aims to gradually replace onboard systems based on mechanical, hydraulic, or pneumatic power sources with electrically powered ones to reduce the weight and costs, optimize energy, and increase the eco-compatibility and reliability of future aircrafts.

A key technological enabler for pursuing these challenging objectives is electro-mechanical actuation. The applicability of electro-mechanical actuators (EMAs) in aerospace has been proved in terms of dynamic performances, but it still entails several concerns in terms of reliability/safety and operation in a harsh environment. In civil aircrafts, EMAs are often avoided for safety-critical functions (flight controls, brakes, landing gears, and nose wheel steering), essentially because the statistical database on the components' fault modes is poor.

This Special Issue is thus focused on advancements and innovations in the design, modelling/simulation, architectural definition, reliability/safety analysis, control, condition-monitoring, and experimental testing of EMAs developed for safety-critical aerospace applications.

The research papers included in this Special Issue will undoubtedly contribute to progress toward the objective of more electric flights.

> **Gianpietro Di Rito** *Editor*

## *Article* **Minimisation of Failure Transients in a Fail-Safe Electro-Mechanical Actuator Employed for the Flap Movables of a High-Speed Helicopter-Plane**

**Gianpietro Di Rito 1,\* , Romain Kovel <sup>2</sup> , Marco Nardeschi <sup>3</sup> , Nicola Borgarelli <sup>3</sup> and Benedetto Luciano <sup>4</sup>**

	- <sup>3</sup> Umbragroup Spa, Via Baldaccini 1, 06034 Foligno, Italy

**Abstract:** The work deals with the model-based characterization of the failure transients of a fail-safe rotary EMA developed by Umbragroup (Italy) for the flap movables of the RACER helicopterplane by Airbus Helicopters (France). Since the reference application requires quasi-static positiontracking with high disturbance-rejection capability, the attention is focused on control hardover faults which determine an actuator runaway from the commanded setpoint. To perform the study, a high-fidelity nonlinear model of the EMA is developed from physical first principles and the main features of health-monitoring and closed-loop control functions (integrating the conventional nested loops architecture with a deformation feedback loop enhancing the actuator stiffness) are presented. The EMA model is then validated with experiments by identifying its parameters by ad-hoc tests. Simulation results are finally proposed to characterize the failure transients in worst case scenarios by highlighting the importance of using a specifically designed back-electromotive damper circuitry into the EMA power electronics to limit the position deviation after the fault detection.

**Keywords:** health monitoring; electro-mechanical actuators; modelling; simulation; testing; flight control; reliability; fault-tolerant systems; failure transient analysis

## **1. Introduction**

The aircraft electrification is surely one of the most important and strategic initiatives currently supporting the innovation of the aviation industry [1,2]. In particular, the *more-electric aircraft* concept entails the gradual replacement of onboard systems based on mechanical, hydraulic, or pneumatic power sources with electrically powered ones, aiming to reduce weight and costs, to optimize energy and to increase the eco-compatibility of future aircrafts [3–6]. Electro-mechanical actuation clearly plays a key role for pursuing these challenging objectives. The applicability of Electro-Mechanical Actuators (EMAs) in aerospace is proven in terms of load and speed performances [7–12], but several reliability concerns still remain open [13–17]. The use of EMAs for safety-critical functions can be thus obtained only by fault-tolerant architectures, which apply hardware redundancies on electrical, electronic, or mechanical parts.

In general terms, depending on how the redundancy is applied, a fault-tolerant function can be maintained after a fault, or it can be lost while avoiding the extension of the fault effects to other functions, so that fail-operative or fail-safe functions are respectively obtained. With reference to flight control functions, this concept can be applied to both movables and actuators. To obtain a fail-operative flight control, different architectures can be used by applying load-level redundancy (splitting the movable into sub-movables and using a fail-safe EMA on each part), actuator-level redundancy (using multiple failsafe EMAs on a single movable), or subsystem-level redundancy (using a single failoperative EMA on a single movable). Any fault-tolerant system necessitates effective

**Citation:** Di Rito, G.; Kovel, R.; Nardeschi, M.; Borgarelli, N.; Luciano, B. Minimisation of Failure Transients in a Fail-Safe Electro-Mechanical Actuator Employed for the Flap Movables of a High-Speed Helicopter-Plane. *Aerospace* **2022**, *9*, 527. https:// doi.org/10.3390/aerospace9090527

Academic Editor: Mark Lowenberg

Received: 22 July 2022 Accepted: 16 September 2022 Published: 19 September 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

health-monitoring algorithms aiming to anticipate (i.e., avoid) or detect/isolate the fault, so that prognostic and diagnostic approaches are respectively defined. The prognostic solution, though potentially overwhelming [18–21], is nowadays far from being applicable to airworthy systems, and diagnostic approaches are typically preferred [22–24]. The diagnostic monitoring requires that the Fault Detection and Isolation (FDI) is implemented and executed in real time by onboard control electronics [25–28], so that, in case of fault, the redundant components or the isolation devices can be engaged [29,30].

The design and the validation of health-monitoring systems play a key role in this context. In a flight control EMA, the FDI output, including the consequent accommodation/compensation of the fault, must be provided with very small latency, and the failure transients must be adequately limited. The development of high-fidelity experimentally validated models of EMAs is of paramount importance for the validation of monitoring functions. Since nonlinearities, disturbances, environment, and loads can significantly affect the actuator response, an in-depth knowledge of both normal and faulty behaviours is required. The crucial problem entails the knowledge of faulty dynamics, especially in complex systems with a huge number of fault modes [31–33]. In the so-called data-driven techniques, this knowledge is achieved via experiments, by artificially injecting the major faults in the EMA and measuring its response [20,34–36]. This method provides accurate predictions, but rigging costs are often prohibitive. In addition, the FDI validation strongly depends on test conditions. As a relevant example, in [20], the mechanical degradation of the ball-screw elements of an aircraft EMA is investigated via a data-driven approach: the lifecycle of a rudder control actuator, including periodical maintenance checks, is simulated by testing a prototype EMA in laboratory environment with alternate endurance and monitoring trials. To accelerate the mechanical degradation, the prototype is intentionally modified with respect to the nominal design, by using a reduced number of recirculating paths in the ball-screw, by removing the anti-rotation device on the output shaft, by applying relevant radial loads, and by progressively removing the lubricant. Discrete-time and continuous-time fault symptoms are then computed by leveraging the EMA outputs via multivariate statistical methods (such as Hotelling's T<sup>2</sup> and Q techniques). The health monitoring demonstrates to be very effective, but the entire experimental activity required a specifically dedicated rig and took seven months. In addition, the experimental campaign did not take into account temperature effects.

In model-based techniques, the knowledge of the faulty dynamics is to a great extent obtained from mathematical models, capable of simulating the fault by physical first principles, and are experimentally validated with reference to normal and/or regime faulty conditions [26–28,37,38]. Oppositely to the data-driven case, this method generally provides less accurate predictions, but it is cost-effective, allows to verify the FDI functionalities in extreme conditions, and (above all) permits to generalize the validity of algorithms to similar equipment (i.e., governed by similar equations). As a relevant example, in [28], the major faults of a primary flight control movable driven by active-active EMAs are addressed via model-based approach: a set of monitoring algorithms are designed using a detailed nonlinear model of the system capable of fault simulation. Robust detection thresholds are determined taking into account parametric and input uncertainties, and the health-monitoring is verified through simulation, by injecting faults in an experimentally validated model of the system.

The basic objective of this work is to validate the monitoring algorithms of the failsafe EMA developed by UmbraGroup (Italy) for the flap movables of the RACER (Rapid And Cost-Efficient Rotorcraft) helicopter-plane by Airbus Helicopters with reference to the model-based analysis of the failure transients related to the hardover of the control electronics (major EMA fault mode). The paper is articulated as follows: the first part is dedicated to the system description and to the EMA modelling; successively, the main features of the closed-loop control and health-monitoring functions are presented. Finally, an excerpt of simulation results is proposed by characterizing the EMA failure transients in selected worst-case scenarios. The results are finally discussed by highlighting the

effectiveness and the most relevant criticalities of the proposed approach with suggestions of possible enhancements. fectiveness and the most relevant criticalities of the proposed approach with suggestions of possible enhancements.

the model-based analysis of the failure transients related to the hardover of the control electronics (major EMA fault mode). The paper is articulated as follows: the first part is dedicated to the system description and to the EMA modelling; successively, the main features of the closed-loop control and health-monitoring functions are presented. Finally, an excerpt of simulation results is proposed by characterizing the EMA failure transients in selected worst-case scenarios. The results are finally discussed by highlighting the ef-

*Aerospace* **2022**, *9*, x FOR PEER REVIEW 3 of 23

### **2. Materials and Methods 2. Materials and Methods**

### *2.1. System Description 2.1. System Description*

The reference EMA is used to control the six flap movables of the RACER helicopterplane, an innovative, high-speed, more-electric air vehicle developed by Airbus Helicopters, Figure 1. The RACER helicopter-plane is designed to reach maximum cruise speed 50% faster than a conventional helicopter (the Velocity Never Exceed, VNE, is 115 m/s) and to consume 15% less fuel per distance at reference cruise speed (90 m/s) [39]. The aerodynamic concept essentially merges a conventional helicopter with a low aspect-ratio box-wing airplane; at cruise speed, the two wing propellers generate thrust and the box-wing contributes to lift, generating low induced drag and minimized interactions with the main rotor flow [40,41], so that the rotor can be slowed by up to 15%, preventing the blades from working with transonic local flow (which reduces performances). The electrical system is based on high-voltage direct current power generation, assuring a consistent weight reduction [42]. The reference EMA is used to control the six flap movables of the RACER helicopterplane, an innovative, high-speed, more-electric air vehicle developed by Airbus Helicopters, Figure 1. The RACER helicopter-plane is designed to reach maximum cruise speed 50% faster than a conventional helicopter (the Velocity Never Exceed, VNE, is 115 m/s) and to consume 15% less fuel per distance at reference cruise speed (90 m/s) [39]. The aerodynamic concept essentially merges a conventional helicopter with a low aspect-ratio box-wing airplane; at cruise speed, the two wing propellers generate thrust and the boxwing contributes to lift, generating low induced drag and minimized interactions with the main rotor flow [40,41], so that the rotor can be slowed by up to 15%, preventing the blades from working with transonic local flow (which reduces performances). The electrical system is based on high-voltage direct current power generation, assuring a consistent weight reduction [42].

**Figure 1.** Flap movables on the RACER helicopter-plane by Airbus Helicopters. **Figure 1.** Flap movables on the RACER helicopter-plane by Airbus Helicopters.

The RACER flaps placed on both box-wing (four) and vertical stabilizers (two) are used to adapt the vehicle attitude, to enhance stability, to optimize the trim configuration, and to abate noise [43–45]. Depending on the vehicle weight, the airspeed, the altitude, and the rotors speed, the wing flaps are deflected to optimize the mean lift coefficient of the main rotor. On the other hand, the flaps on the vertical fins are used to eliminate the residual yawing torque generated by the propellers, by assuring that they only contribute to propulsion during cruise [43]. Given these basic flight control functions, the design of the closed-loop position control of the flap movables is mainly driven by disturbance rejection requirements (i.e., the capability to minimize the position deviation from the com-The RACER flaps placed on both box-wing (four) and vertical stabilizers (two) are used to adapt the vehicle attitude, to enhance stability, to optimize the trim configuration, and to abate noise [43–45]. Depending on the vehicle weight, the airspeed, the altitude, and the rotors speed, the wing flaps are deflected to optimize the mean lift coefficient of the main rotor. On the other hand, the flaps on the vertical fins are used to eliminate the residual yawing torque generated by the propellers, by assuring that they only contribute to propulsion during cruise [43]. Given these basic flight control functions, the design of the closed-loop position control of the flap movables is mainly driven by disturbance rejection requirements (i.e., the capability to minimize the position deviation from the commanded setpoint under external disturbances).

manded setpoint under external disturbances). It is worth noting that in the RACER helicopter-plane, the flaps are not used for manoeuvrability (trajectory control is managed through the cyclic stick, as for conventional It is worth noting that in the RACER helicopter-plane, the flaps are not used for manoeuvrability (trajectory control is managed through the cyclic stick, as for conventional helicopters), so that they are classified as secondary flight controls. The flap EMAs are thus designed to be fail-safe systems in such a way that, after a major fault, the actuator is still capable of maintaining the flap movable at a fixed deflection (last or neutral position, depending on the fault mode), providing an adequate torsional stiffness to avoid flutter concerns.

Each flap EMA is composed of two parts, Figure 2: an electromechanical rotary actuator (FLap Actuator, FLA) and a control electronic box (Actuator Control Electronics, concerns.

concerns.

ACE). The equipment locations on the RACER helicopter-plane are shown in Figure 3 (note that the reported layout also depicts flaps on the horizontal tail since they have been initially included in the flight control system [41,43], but then eliminated from the final design). (note that the reported layout also depicts flaps on the horizontal tail since they have been initially included in the flight control system [41,43], but then eliminated from the final design). ACE). The equipment locations on the RACER helicopter-plane are shown in Figure 3 (note that the reported layout also depicts flaps on the horizontal tail since they have been initially included in the flight control system [41,43], but then eliminated from the final design).

Each flap EMA is composed of two parts, Figure 2: an electromechanical rotary actuator (FLap Actuator, FLA) and a control electronic box (Actuator Control Electronics, ACE). The equipment locations on the RACER helicopter-plane are shown in Figure 3

Each flap EMA is composed of two parts, Figure 2: an electromechanical rotary actuator (FLap Actuator, FLA) and a control electronic box (Actuator Control Electronics,

helicopters), so that they are classified as secondary flight controls. The flap EMAs are thus designed to be fail-safe systems in such a way that, after a major fault, the actuator is still capable of maintaining the flap movable at a fixed deflection (last or neutral position, depending on the fault mode), providing an adequate torsional stiffness to avoid flutter

helicopters), so that they are classified as secondary flight controls. The flap EMAs are thus designed to be fail-safe systems in such a way that, after a major fault, the actuator is still capable of maintaining the flap movable at a fixed deflection (last or neutral position, depending on the fault mode), providing an adequate torsional stiffness to avoid flutter

*Aerospace* **2022**, *9*, x FOR PEER REVIEW 4 of 23

*Aerospace* **2022**, *9*, x FOR PEER REVIEW 4 of 23

**Figure 2.** Flap EMA layout: (**a**) FLA; (**b**) ACE. **Figure 2.** Flap EMA layout: (**a**) FLA; (**b**) ACE. **Figure 2.** Flap EMA layout: (**a**) FLA; (**b**) ACE.

**Figure 3.** Locations of the EMA parts on the RACER helicopter-plane (yellow: ACEs; blue: FLAs). **Figure 3.** Locations of the EMA parts on the RACER helicopter-plane (yellow: ACEs; blue: FLAs). **Figure 3.** Locations of the EMA parts on the RACER helicopter-plane (yellow: ACEs; blue: FLAs).

The ACE includes three electronic boards, Figures 4 and 5: The ACE includes three electronic boards, Figures 4 and 5: The ACE includes three electronic boards, Figures 4 and 5:


(Back Electro-Motive Force) damper circuitry, Figure 5. The activation of the BEMF damper circuitry in the PWR board is obtained by a logic signal named *system validity* (SV in Figures 4 and 5), which derives from an "AND" operator applied to the *local validity* signals provided by the two boards (LVCOM and LVMON in Figures 4 and 5). When SV is true, the power bridge thyristors are opened, and the damper thyristor is closed (Figure 5), so that the motor phases are shorted to the ground, and an electromagnetic damping torque is developed and transmitted to the EMA output shaft. This strategy permits significantly limiting the failure transients related to major faults (e.g., the control electronics hardover) since the unavoidable delays needed

to achieve the full engagement of the motor brakes can determine an excessive deviation from the commanded setpoint with potentially dangerous concerns due to the impact on mechanical stops. full engagement of the motor brakes can determine an excessive deviation from the commanded setpoint with potentially dangerous concerns due to the impact on mechanical stops. manded setpoint with potentially dangerous concerns due to the impact on mechanical stops.

**Legend**

RFj : j-th rotation feedback SV : System Validity Tk : k-th thyristor (k=1,…,6) TFj : j-th temperature feedback Vx : x-th phase voltage

The activation of the BEMF damper circuitry in the PWR board is obtained by a logic signal named *system validity* (SV in Figures 4 and 5), which derives from an "AND" operator applied to the *local validity* signals provided by the two boards (LVCOM and LVMON in Figures 4 and 5). When SV is true, the power bridge thyristors are opened, and the damper thyristor is closed (Figure 5), so that the motor phases are shorted to the ground, and an electromagnetic damping torque is developed and transmitted to the EMA output shaft. This strategy permits significantly limiting the failure transients related to major faults (e.g., the control electronics hardover) since the unavoidable delays needed to achieve the

The activation of the BEMF damper circuitry in the PWR board is obtained by a logic signal named *system validity* (SV in Figures 4 and 5), which derives from an "AND" operator applied to the *local validity* signals provided by the two boards (LVCOM and LVMON in Figures 4 and 5). When SV is true, the power bridge thyristors are opened, and the damper thyristor is closed (Figure 5), so that the motor phases are shorted to the ground, and an electromagnetic damping torque is developed and transmitted to the EMA output shaft. This strategy permits significantly limiting the failure transients related to major faults (e.g., the control electronics hardover) since the unavoidable delays needed to achieve the full engagement of the motor brakes can determine an excessive deviation from the com-

*Aerospace* **2022**, *9*, x FOR PEER REVIEW 5 of 23

*Aerospace* **2022**, *9*, x FOR PEER REVIEW 5 of 23

**Figure 4.** ACE: schematics of boards interfaces. **Figure 4.** ACE: schematics of boards interfaces. **Figure 4.** ACE: schematics of boards interfaces.

PDMON SV PDCOM

**Figure 5.** Schematics of the main sections of the PWR board. **Figure 5.** Schematics of the main sections of the PWR board. **Figure 5.** Schematics of the main sections of the PWR board.

Both COM and MON boards are controlled by Texas Instruments TMS570LC4357- EP ARM-based microcontrollers [46] using a 10 kHz sample rate for the digital signal processing. The FLA basically includes, Figure 6: Both COM and MON boards are controlled by Texas Instruments TMS570LC4357- EP ARM-based microcontrollers [46] using a 10 kHz sample rate for the digital signal processing. Both COM and MON boards are controlled by Texas Instruments TMS570LC4357-EP ARM-based microcontrollers [46] using a 10 kHz sample rate for the digital signal processing. The FLA basically includes, Figure 6:


• two power-off electromagnetic brakes used to block in position the EMA after a major fault detection; digital converter Analog Devices AD2S1210 ±π rad 2×10−4 rad Magnetic encoder

• a three-phase Permanent Magnet Synchronous Machine (PMSM) with surfacemounted magnets and sinusoidal back-electromotive forces, driven via Field-Ori-

• two motor rotation sensors: a resolver interfaced with the COM board and a magnetic

• a dual magnetic encoder for the output shaft rotation sensing, interfaced to both

• two power-off electromagnetic brakes used to block in position the EMA after a major

• an innovative Umbragroup-patented differential ball-screw mechanism implementing the mechanical power conversion from motor to output shaft, which, if compared with conventional gearboxes, assures a high gear ratio (more than 500) with mini-

Table 1 reports the main characteristics of the sensors used for the closed-loop control

**Component Model Range Accuracy**  Current sensor Allegro ACS723LLCTR-10AB-T ±10 A 0.1 A Resolver Tamagawa TS2610N171E64 ±π rad 4×10−4 rad

encoder interfaced with the MON board (RFCOM and RFMON in Figure 4);

mum backlash (less than 0.05 deg) and superior efficiency (about 95%).

*Aerospace* **2022**, *9*, x FOR PEER REVIEW 6 of 23

COM and MON boards (PFCOM and PFMON in Figure 4); • two temperature sensors (TFCOM and TFMON in Figure 4);

and health-monitoring functions of the actuator [47-51].

ented Control (FOC) technique;

fault detection;

**Table 1.** EMA control sensors data.

Resolver analog-to-

• an innovative Umbragroup-patented differential ball-screw mechanism implementing the mechanical power conversion from motor to output shaft, which, if compared with conventional gearboxes, assures a high gear ratio (more than 500) with minimum backlash (less than 0.05 deg) and superior efficiency (about 95%). (motor) Analog Devices ADA4571 ±π rad 4×10−4 rad Duplex magnetic encoder (output) RLS AksIM-2 ±0.157 rad 1.7×10−3 rad

**Figure 6.** Schematics of the FLA electro-mechanical section. **Figure 6.** Schematics of the FLA electro-mechanical section.

*2.2. Nonlinear Dynamic Modelling*  The EMA health-monitoring algorithms have been designed with the essential sup-Table 1 reports the main characteristics of the sensors used for the closed-loop control and health-monitoring functions of the actuator [47–51].



## *2.2. Nonlinear Dynamic Modelling*

The EMA health-monitoring algorithms have been designed with the essential support of the dynamic simulation by artificially injecting major system faults in a detailed nonlinear model of the actuator, developed via physical first principles and validated through experiments. This approach is of paramount importance, especially for the failure transient characterisation which is often unfeasible (or problematic) via testing.

The model of the RACER flaps EMA is essentially composed of:

	- # FOC current dynamics;
	- # multi-harmonic modelling of PMSM torque disturbances (due to cogging effects [52–55] and/or BEMF waveform distortions);
	- # 2-degree-of-freedom mechanical transmission with equations of motions related to motor and output rotations;
	- # sliding friction on motors and output shaft, described via combined "Coulomb– tanh" model [56,57];
	- # mechanical freeplay [21];
	- # internal stiffness dependence on output shaft position;
	- # Clarke-Park transforms for the FOC technique implementation;
	- # sensor errors and nonlinearities (bias, noise, resolution);
	- # command nonlinearities (saturation, rate limiting);
	- # digital signal processing at 10 kHz sampling rate;
	- # control hardover fault simulation, implying that the voltage demands on both quadrature and direct axes suddenly assume and maintain random values, so that the EMA motion is out of control (as a worst case scenario, the quadrature voltage is set to saturation value, while the direct voltage is set to zero).

## 2.2.1. Electro-Mechanical Section of the Model

The electro-mechanical section of the model, schematically represented in Figure 6, is governed by Equations (1)–(10),

$$\mathbf{V\_{abc}} = R\mathbf{I\_{abc}} + L\dot{\mathbf{I}\_{abc}} + \mathbf{e\_{abc}} \tag{1}$$

$$\mathbf{e\_{abc}} = \lambda\_m n\_d \dot{\theta}\_m \left[ \sin(n\_d \theta\_m) \quad \sin \left( n\_d \theta\_m - \frac{2}{3} \pi \right) \quad \sin \left( n\_d \theta\_m + \frac{2}{3} \pi \right) \right]^T,\tag{2}$$

$$J\_m \ddot{\theta}\_m = T\_m + T\_b - T\_{sfm} \tanh\left(\frac{\dot{\theta}\_m}{\omega\_{sfm}}\right) - d\_{vfm} \dot{\theta}\_m - \frac{d\_s \dot{\theta}\_s + k\_s \theta\_s}{\tau\_\mathcal{S}},\tag{3}$$

$$J\_o \ddot{\theta}\_o = T\_{aer} - T\_{sfo} \tanh\left(\frac{\dot{\theta}\_o}{\omega\_{sfo}}\right) - d\_{vfo} \dot{\theta}\_o + d\_s \dot{\theta}\_s + k\_s \theta\_{s\prime} \tag{4}$$

$$T\_{\mathfrak{M}} = \lambda\_{\mathfrak{M}} n\_d \left[ I\_d \sin(n\_d \theta\_{\mathfrak{m}}) + I\_b \sin\left(n\_d \theta\_{\mathfrak{m}} - \frac{2}{3}\pi\right) + I\_c \sin\left(n\_d \theta\_{\mathfrak{m}} + \frac{2}{3}\pi\right) \right] + \sum\_{j=1}^{M} T\_{\mathfrak{M}j} \cdot \sin\left(n\_{\mathrm{dif}} \theta\_{\mathfrak{m}}\right), \tag{5}$$

$$k\_s = k\_{\rm smin} + \gamma\_k (\theta\_o - \theta\_{o\rm max})^2,\tag{6}$$

$$T\_b = \begin{cases} 0 & t < t\_{FC} \\ -k\_b[\theta\_m - \theta\_m(t\_{FC})] - d\_b[\dot{\theta}\_m - \dot{\theta}\_m(t\_{FC})] & t \ge t\_{FC} \end{cases} \tag{7}$$

$$\dot{\theta}\_{\text{s}} = \begin{cases} -\frac{k\_s}{d\_s} \theta\_{\text{s}} & \left| \theta\_{\text{g}} - \theta\_o \right| < \varepsilon\_p \\ \dot{\theta}\_{\text{g}} - \dot{\theta}\_o & \left| \theta\_{\text{g}} - \theta\_o \right| \ge \varepsilon\_p \end{cases} \tag{8}$$

$$\theta\_s = \begin{cases} \begin{array}{c} \int \dot{\theta}\_s dt \\ \theta\_{\mathcal{S}} - \theta\_o - \varepsilon\_p \text{sgn}\left(\theta\_{\mathcal{S}} - \theta\_o\right) \\ \end{array} \begin{array}{c} \begin{array}{c} \left| \theta\_{\mathcal{S}} - \theta\_o \right| < \varepsilon\_p \\ \theta\_{\mathcal{S}} - \theta\_o \end{array} \right. \\ \end{array} \tag{9}$$

$$
\dot{\theta}\_{\mathcal{S}} = \frac{\theta\_m}{\tau\_{\mathcal{S}}} \, \, \, \, \tag{10}
$$

where **Vabc** = [*Va*, *V<sup>b</sup>* , *Vc*] *T* is the applied voltages vector, **Iabc** = [*Ia*, *I<sup>b</sup> , Ic*] *T* is the phase currents vector, **eabc** is the back-electromotive forces vector (sinusoidal BEMF waveforms are assumed), *λ<sup>m</sup>* is the magnet flux linkage, *R* and *L* are the resistance and the inductance of motor phases, respectively, *n<sup>d</sup>* is the number of rotor pole pairs, *θ<sup>m</sup>* is the motor rotation, *θ<sup>g</sup>* is the theoretical rotation imposed by a rigid ball-screw drivetrain, *θ<sup>s</sup>* is the torsional deformation referred to the first structural mode of the EMA, *θ<sup>o</sup>* is the output rotation, *J<sup>m</sup>* and *J<sup>o</sup>* are the motor and output inertias, respectively, *τ<sup>g</sup>* is the gear ratio of the differential ball-screw mechanism, *ε<sup>p</sup>* is the internal freeplay, *T<sup>m</sup>* is the motor torque, *Thd j* and *nhd j* are the amplitude and mechanical period indices, respectively, related to the *j*-th (*j* = 1, . . . , *M* where *M* is an integer number) harmonic torque disturbance contribution, *Taer* is the aerodynamic hinge moment, *T<sup>b</sup>* is the brakes torque, *k<sup>b</sup>* and *d<sup>b</sup>* are the torsional stiffness and damping of the brakes, respectively, *tFC* is the time at which the fault compensation occurs, *k<sup>s</sup>* and *d<sup>s</sup>* are the torsional stiffness and damping, respectively, referred to the first

.

structural mode of the EMA, *ks*min is the minimum internal stiffness, *γ<sup>k</sup>* is the parameter defining the stiffness variation with respect to output position, *dvfm* and *dvfo* are the viscous friction coefficients related to the motor and output shafts, respectively, while *Tsfm*, *ωsfm*, *Tsfo*, and *ωsfo* are the parameters of the "Coulomb–tanh" models simulating the sliding friction on motor and output shaft, respectively.

Concerning the aerodynamic hinge moment applied on the EMA (*Taer*, in Equation (4)), the requirements indicate that, apart from static loadings, two contributions of dynamic loads must be taken into account for the performance analysis: a deterministic one (related to the helicopter-plane motion, main rotor speed and angle, wing propeller speeds), in which harmonic loads of specific amplitudes and frequencies are superimposed, and a non-deterministic one, including gust loads and harmonic loads of constant amplitudes randomly applied along the the position-tracking frequency range. In this work, the study is focused on the vertical stabilizer flaps, since it represents the worst-case scenario for the EMAs employed in the RACER helicopter-plane, Table 2.

**Table 2.** Loads on vertical stabilizer FLA at VNE (derived from CFD analyses by Airbus Helicopters, worst-case scenario, all positions).


It is worth noting that the proposed model represents a balance between prediction accuracy, objectives of the study, and complexity of the model itself. More accurate simulations could include sophisticated friction models [56,57] and iron losses in the motor [58,59], but the inclusion of these features would entail minor effects for the examined application. In particular, the motor iron losses have been neglected because they depend on electrical frequency, which is relatively small in the position-tracking frequency range (<50 Hz, Table 2). On the other hand, more accurate friction models (including load and temperature dependence) could enhance the simulation, but a simplified approach has been preferred both for the lack of detailed information and to limit the number of model parameters.

The EMA model has been entirely developed in the Matlab–Simulink–Stateflow environment, and the numerical solution is obtained by a Runge–Kutta method with 10−<sup>6</sup> s integration step. The choice of a fixed-step solver is not strictly related to the objectives of this work in which the model (once experimentally validated) is used for "off-line" simulations characterising the EMA failure transients, but it has been selected for the next steps of the project, when the algorithms for the closed-loop control and the health-monitoring will be implemented in the ACE boards via the automatic Matlab compiler and executed in "real-time".

The parameters of the electro-mechanical section of the model are given in Table 3.


**Table 3.** Parameters of the electro-mechanical section of the model.


## **Table 3.** *Cont.*

#### 2.2.2. Electronic Section of the Model ( ) ( ) ( ) ( ) ( ) sgn *j j jj sat PI PI sat y y yy* <sup>⎪</sup> <sup>≥</sup> ⎪⎩

The closed-loop control of the RACER flap EMAs is schematically represented in Figure 7. The position-tracking architecture integrates the conventional three nested loops on motor currents, motor speed, and output position [8,9] with a deformation feedback loop ("Stiffness Enhancement System, SES" block in Figure 7) and a model-based correction of voltage commands, aiming to decouple the currents dynamics from the motor motion ("electro-mechanical decoupler" block in Figure 7). where is the discrete-time operator, () is the regulator input (tracking error), () is the regulator output, ூ () is the saturator block input (proportional–integral with respect to error, if no saturation is present), while () and ூ () are the proportional and integral gains, ௐ () is the back-calculation anti-windup gain, ௦௧ () is the saturation limit, and *Ts* is the sampling time.

Concerning the SES loop, its basic objective is to enhance the loads disturbance rejection in the frequency range where the first resonant pulsation of the ball-screw mechanism

**Figure 7.** Closed-loop control scheme of the actuator. **Figure 7.** Closed-loop control scheme of the actuator.

The FOC technique implemented in the COM board applies the direct and inverse Clark–Park transforms [60] via Equations (11)–(13),

$$\mathbf{x}\_{\mathbf{a}\otimes\mathbf{y}} = \mathbf{T}\_{\mathbf{C}} \mathbf{x}\_{\mathbf{a}\mathbf{b}\mathbf{c}} = \sqrt{\frac{2}{3}} \begin{bmatrix} 1 & -1/2 & -1/2 \\ 0 & \sqrt{3}/2 & -\sqrt{3}/2 \\ \sqrt{2}/2 & \sqrt{2}/2 & \sqrt{2}/2 \end{bmatrix} \mathbf{x}\_{\mathbf{a}\mathbf{b}\mathbf{c}} \tag{11}$$

$$\mathbf{x\_{dqz}} = \mathbf{T\_P} \mathbf{x\_{\alpha\beta}} \\ \mathbf{y} = \begin{bmatrix} \cos(n\_d \theta\_m) & \sin(n\_d \theta\_m) & 0 \\ -\sin(n\_d \theta\_m) & \cos(n\_d \theta\_m) & 0 \\ 0 & 0 & 1 \end{bmatrix} \mathbf{x\_{\alpha\beta}} \tag{12}$$

$$\mathbf{x\_{dqz}} = \mathbf{T\_P}\mathbf{T\_C}\mathbf{x\_{abc}} \Leftrightarrow \mathbf{x\_{abc}} = \left(\mathbf{T\_P}\mathbf{T\_C}\right)^T \mathbf{x\_{dqz}}\tag{13}$$

where **x**αβγ = [*xα*, *xβ*, *xγ*] *T* , **xabc** = [*xa*, *x<sup>b</sup>* , *xc*] *T* , and **xdqz** = [*x<sup>d</sup>* , *xq*, *xz*] *<sup>T</sup>* are generic three-phase vectors in the Clarke, Park, and stator reference frames, respectively, while **T<sup>C</sup>** and **T<sup>P</sup>** are the Clarke and Park transforms.

The digital regulators on position, speed, and currents implement proportional/integral actions on tracking error signals, plus anti-windup functions with back-calculation technique [61] to compensate for commands saturation. Each *j*-th (with *j* = *θ*, ω, and *I* indicating the position, speed, and currents loops, respectively) digital regulator is governed by Equations (14) and (15):

$$y\_{PI}^{(j)} = k\_P^{(j)} \varepsilon^{(j)} + \frac{k\_I^{(j)}}{z-1} \left[ \varepsilon^{(j)} + k\_{AW}^{(j)} \left( y^{(j)} - y\_{PI}^{(j)} \right) \right] \tag{14}$$

$$\mathbf{y}^{(j)} = \begin{cases} \mathbf{y}\_{PI}^{(j)} \\ \mathbf{y}\_{\text{sat}}^{(j)} \text{sgn}\left(\mathbf{y}\_{PI}^{(j)}\right) \end{cases} \begin{vmatrix} \mathbf{y}\_{PI}^{(j)} \\ \mathbf{y}\_{PI}^{(j)} \end{vmatrix} < \mathbf{y}\_{\text{sat}}^{(j)} \tag{15}$$

where *z* is the discrete-time operator, *ε* (*j*) is the regulator input (tracking error), *y* (*j*) is the regulator output, *y* (*j*) *PI* is the saturator block input (proportional–integral with respect to error, if no saturation is present), while *k* (*j*) *P* and *k* (*j*) *I* are the proportional and integral gains, *k* (*j*) *AW* is the back-calculation anti-windup gain, *y* (*j*) *sat* is the saturation limit, and *T<sup>s</sup>* is the sampling time.

Concerning the SES loop, its basic objective is to enhance the loads disturbance rejection in the frequency range where the first resonant pulsation of the ball-screw mechanism is located (according to FEM analyses performed by Umbragroup, from 70 to 90 Hz, depending on output shaft position, Equation (6)). The control task is achieved by superimposing to the current demand generated by the speed regulator (*Iqc*, in Figure 7) an additional one (*IqSES*, in Figure 7) that depends on the torsional deformation (*δ<sup>f</sup>* ) reconstructed by the motion feedbacks (*θmf* and *θof* in Figure 7), Equations (16)–(18).

$$I\_{qd} = I\_{qc} + I\_{q\text{SES}} \tag{16}$$

$$
\ddot{I}\_{qSES} = -a\_{SES}\dot{I}\_{qSES} - b\_{SES}I\_{qSES} - k\_{SES}\dot{\delta}\_{f\prime} \tag{17}
$$

$$
\delta\_f = \frac{\theta\_{mf}}{\tau\_\mathcal{R}} - \theta\_{of\ \prime} \tag{18}
$$

The structure of the current demand regulator (a second-order system responding to deformation rate input, Equations (17) and (18)) is defined by pursuing the following requirements:


quirements:

tracking, etc.);

to external loads.

tion (17) are set in such a way that

an increase of EMA stiffness, enhancing the disturbance rejection capabilities related to external loads. • by tuning *aSES* and *bSES*, the phase response of the SES current demand (*IqSES*) with respect to torsional deformation is about −180° from 70 to 90 Hz;

To fulfill these objectives, the positive-defined parameters, *kSES, aSES*, and *bSES* in Equation (17) are set in such a way that • by tuning *kSES*, the SES current demand (*IqSES*) implies an effective compensation without affecting the control stability.

To fulfill these objectives, the positive-defined parameters, *kSES, aSES*, and *bSES* in Equa-

*Aerospace* **2022**, *9*, x FOR PEER REVIEW 11 of 23

by the motion feedbacks (*θmf* and *θof* in Figure 7), Equations (16)–(18).

pending on output shaft position, Equation (6)). The control task is achieved by superimposing to the current demand generated by the speed regulator (*Iqc*, in Figure 7) an additional one (*IqSES*, in Figure 7) that depends on the torsional deformation (*δf*) reconstructed

*qSES SES qSES SES qSES SES f I aI bI k* =− − −

*mf f of g*

 θ

• the loop shall not affect the EMA low-frequency behaviour (maxima loads, position

• the loop shall generate demands only in the frequency range where the first resonant pulsation of the ball-screw mechanism is located, and the compensation shall imply an increase of EMA stiffness, enhancing the disturbance rejection capabilities related

The structure of the current demand regulator (a second-order system responding to deformation rate input, Equations (17) and (18)) is defined by pursuing the following re-

θ

τ

δ


Figure 8 shows the Bode diagram of the transfer function defined in Equation (19), which relates in the Laplace domain (i.e., *s* represents the complex variable) the SES current demand with the reconstructed deformation feedback, ( ) ( ) <sup>2</sup> *qSES SES f SES SES I s k s* δ *s s a sb* = − + + , (19) It can be noted that the design implies that from 70 to 90 Hz the phase response

$$\frac{I\_{q\text{SES}}(s)}{\delta\_f(s)} = -\frac{k\_{\text{SES}}s}{s^2 + a\_{\text{SES}}s + b\_{\text{SES}}},\tag{19}$$

*qd qc qSES I II* = + , (16)

= − , (18)

, (17)

δ

**Figure 8.** Bode diagram of the SES regulator (0 dB = 1 A/m): (**a**) behaviour from low to high frequencies; (**b**) detail in the frequency range where the first resonant pulsation of the ball-screw drivetrain is located.

It can be noted that the design implies that from 70 to 90 Hz the phase response ranges from −175◦ to −190◦ , and the regulator gain achieves its maximum, while it tends to be negligible at both low and high frequencies.

Regarding the currents-motion decoupling, it is obtained via Equations (20) and (21):

$$V\_d = V\_{d\varepsilon} - L\_q n\_d I\_{qf} \dot{\theta}\_{mf} \tag{20}$$

$$V\_q = V\_{q\varepsilon} + \left(\sqrt{\frac{3}{2}}\lambda\_m + L\_d I\_{df}\right) n\_d \dot{\theta}\_{mf\prime} \tag{21}$$

where *Vdc* and *Vqc* are the direct and quadrature voltage demands generated by the current regulators, and *L<sup>d</sup>* and *L<sup>q</sup>* are the inductances on the direct and quadrature axes (in the reference PMSM, having surface-mounted magnets, *L<sup>d</sup>* = *L<sup>q</sup>* = *L*). The currents dynamics imposed by the FOC technique implies that in the PMSM rotor frame (Equation (13)),

$$
\Delta \dot{I}\_d = V\_d - RI\_d + L n\_d I\_q \dot{\theta}\_{m\prime} \tag{22}
$$

$$L\dot{I}\_q = V\_q - RI\_q - \left(\sqrt{\frac{3}{2}}\lambda\_m + LI\_d\right)n\_d\dot{\theta}\_{m\nu} \tag{23}$$

Thus, by substituting Equations (20) and (21) into Equations (22) and (23), we have

$$L\dot{I}\_d = V\_{dc} - RI\_d + Ln\_d \left( I\_q \dot{\theta}\_m - I\_{qf} \dot{\theta}\_{mf} \right), \tag{24}$$

$$\mathrm{L}\dot{I}\_{q} = V\_{q\complement} - \mathrm{R}I\_{q} - \sqrt{\frac{3}{2}}\lambda\_{m}n\_{d}\left(\dot{\theta}\_{m} - \dot{\theta}\_{mf}\right) - \mathrm{L}n\_{d}\left(\mathrm{I}\_{d}\dot{\theta}\_{m} - \mathrm{I}\_{df}\dot{\theta}\_{mf}\right),\tag{25}$$

Now, if the sensor dynamics imply minor phase delays and/or attenuations (*θmf* ≈ *θm*, *Iqf* ≈ *Iq*, *Idf* ≈ *I<sup>d</sup>* ), the residuals terms at second hands in Equations (24) and (25) can be neglected, so that the currents dynamics on both direct and quadrature axes behave independently and are decoupled from the rotor motion.

To protect the system from major faults and to permit its reversion into a fail-safe configuration (EMA with engaged brakes, maintaining the flap at fixed deflection), the following set of health-monitoring algorithms are executed by the MON board:


For the examined application, the most feared EMA failure is the control hardover, i.e., an electronic fault for which the COM board applies and maintains random voltage demands on both quadrature and direct axes, so that the actuator motion is out of control. The coverage of this failure is here provided by the OSM, whose working flow chart is reported in Figure 9. The OSM fault flag (Fmon) is generated by elaborating as fault symptom the amplitude of the speed feedback signal (*ω*mon) at the *k*-th monitoring sample (processed at 10 kHz rate): if the fault symptom is greater than a pre-defined threshold (*ωth*), a fault counter (*cmon*) is increased by 2; if the threshold is not exceeded, the fault counter is decreased by 1 if it is positive at the previous step, otherwise it is held at 0. The fault is thus detected when the fault counter exceeds a pre-defined value (*cmon* max, which basically defines the OSM FDI latency).

The parameters of the electronic section of the model are given in Table 4.


**Table 4.** Parameters of the electronic section of the model.

ூ

ௐ

௦௧

ூ

ௐ

௦௧

ூ

ௐ

௦௧


**Table 4.** *Cont.*

**Figure 9.** Fault detection logics of the OSM. **Figure 9.** Fault detection logics of the OSM.

### The parameters of the electronic section of the model are given in Table 4. 2.2.3. Fault Simulation

2.2.3. Fault Simulation

**Table 4.** Parameters of the electronic section of the model. **Parameter Meaning Value Unit**  *Ts* Digital control sample time (all regulators) 10−4 s (ఏ) Proportional gain of the position regulator 1.58 × 104 1/s As previously mentioned, the basic objective of the work is to validate the EMA health-monitoring algorithms with reference to the control hardover fault (worst-case failure), by particularly focusing on the failure transient characterisation. The model has been developed as a finite-state machine by using Matlab-Simulink-Stateflow charts and logics so that the simulations of both hardover fault and the subsequent activation of the back-electromotive circuitry are integrated in the EMA simulator.

(ఏ) Integral gain of the position regulator 1.1 × 105 1/s2 (ఏ) Anti-windup gain of the position regulator 0.69 s The hardover fault is simulated by Equation (26), so that when the fault is injected (*t* = *tFI*), the direct and quadrature voltages are switched from the values demanded by the EMA control laws (Equations (20) and (21)) to zero and saturation values, respectively:

$$V\_d = \begin{cases} V\_{dc} - L n\_d I\_{d\bar{f}} \dot{\theta}\_{mf} & t < t\_{\text{FI}} \\ 0 & t \ge t\_{\text{FI}} \end{cases}; \; V\_q = \begin{cases} V\_{qc} + \left( \sqrt{\frac{2}{2}} \lambda\_{\text{in}} + L I\_{d\bar{f}} \right) n\_d \dot{\theta}\_{mf} & t < t\_{\text{FI}} \\ V\_{\text{max}} & t \ge t\_{\text{FI}} \end{cases},\tag{26}$$

(ఠ) Saturation limit of the speed regulator 4 A (ூ) Proportional gain of the current regulators 2.78 V/A (ூ) Integral gain of the current regulators 4.1 × 103 V/(A s) As described in Section 2.1, the PWR board of the EMA includes a BEMF damper circuitry, which, in case of a detected fault, imposes that the motor phases are shorted to the ground (so that an electromagnetic damping torque is developed and transmitted

> As previously mentioned, the basic objective of the work is to validate the EMA health-monitoring algorithms with reference to the control hardover fault (worst-case failure), by particularly focusing on the failure transient characterisation. The model has been developed as a finite-state machine by using Matlab-Simulink-Stateflow charts and logics so that the simulations of both hardover fault and the subsequent activation of the back-

(ூ) Anti-windup gain of the current regulators 150 A/V

(ூ) Saturation limit of the current regulators 28 V

*bSES* SES regulator parameter 2 2.37 × 105 rad2/s2 *kSES* SES regulator gain 103 A/(m s) *ωth* OSM fault symptom threshold 0.0175 rad/s *cmon* max OSM fault counter threshold 250 --

electromotive circuitry are integrated in the EMA simulator.

to the EMA output shaft). In the model, the BEMF damper activation is simulated via Equation (27), where *tFD* is the time at which the fault is detected by the OSM.

$$\mathbf{V\_{abc}} = \begin{cases} (\mathbf{T\_F T\_C})^T \mathbf{V\_{dqz}} & t < t\_{FD} \\ 0 & t \ge t\_{FD} \end{cases} \tag{27}$$

## *2.3. Experimental Test Campaign for the Model Validation*

To substantiate the failure transient analysis presented and discussed in Section 3, the EMA model has been experimentally validated through a specific test campaign carried out at the Umbragroup facilities. In particular, the following tests have been performed, aiming to identify the model parameters reported in Table 3:

	- # Test 1 (blocked motor with engaged brakes): chirp wave inputs are given as direct voltage demand, while the quadrature voltage is set to zero, aiming to identify motor phase resistance and inductance (*R* and *L*). The test is repeated at a different position of the PMSM rotor to verify that the phase inductance does not significantly depend on motor angle (assumption of the model);
	- # Test 2 (blocked motor with engaged brakes): step inputs of different amplitudes are given to the quadrature voltage demand, while the direct voltage is set to zero, aiming to confirm the values of motor phase resistance and inductance. The test is repeated at different position of the PMSM rotor;
	- # Test 3 (free-wheeling motor with disengaged brakes and open phases): the PMSM rotor is dragged by an external motor at different speed amplitudes and the phase-to-phase BEMF is measured, aiming to identify the motor flux linkage (*λm*) and to eventually highlight higher harmonic components in the BEMF waveform;
	- # Test 4 (blocked motor with engaged brakes): current loop tracking is tested by providing square-wave inputs of different amplitudes as quadrature current demand, while the direct current is set to zero, aiming to identify the damping and stiffness of the brakes (*d<sup>b</sup>* and *k<sup>b</sup>* );
	- # Test 5 (disengaged brakes): speed loop tracking is tested, by providing squarewave inputs of different amplitudes as speed demand, aiming to identify the torque disturbance parameters (*M, Thd*1, *nhd*1, *Thd*2, *nhd*2, *Thd*<sup>3</sup> and *nhd*3), the viscous damping coefficients (*dvfm* and *dvfo*), the parameters of the sliding friction models (*Tsfm*, *Tsfo ωsfm* and *ωsfo*), and the actuator inertias (*J<sup>m</sup>* and *Jo*).

At the current stage of the campaign, loaded tests and position-loop tests have not been carried out, but they have been planned for the future steps of the research, mainly to confirm the predictions of the resonant frequency of the ball-screw drivetrain, currently estimated through FEM analyses.

## **3. Results**

## *3.1. Experimental Validation of the Model*

An excerpt of the results obtained during the model validation campaign is reported from Figures 10–13.

Figures 10 and 11 are devoted to the identification of the electrical parameters of the motor phases (i.e., resistance and inductance), and it can be noted that the model succeeds in predicting the hardware response in both steady-state and dynamic operations. The repetition of tests at different positions of the PMSM rotor provided essentially identical results, thus confirming that the position-dependence of electrical parameters is negligible (basic assumption of the model).

**3. Results** 

**3. Results** 

voltage demand).

voltage demand).

(basic assumption of the model).

(basic assumption of the model).

*3.1. Experimental Validation of the Model* 

*3.1. Experimental Validation of the Model* 

*Aerospace* **2022**, *9*, x FOR PEER REVIEW 15 of 23

An excerpt of the results obtained during the model validation campaign is reported

**Figure 10.** Response to Test 1 (open-loop, engaged brakes, chirp wave input applied on direct voltage demand, ±2 V ranging from 10 Hz, at t = 0 s, to 10 kHz, at t = 10 s). **Figure 10.** Response to Test 1 (open-loop, engaged brakes, chirp wave input applied on direct voltage demand, ±2 V ranging from 10 Hz, at t = 0 s, to 10 kHz, at t = 10 s). **Figure 10.** Response to Test 1 (open-loop, engaged brakes, chirp wave input applied on direct voltage demand, ±2 V ranging from 10 Hz, at t = 0 s, to 10 kHz, at t = 10 s).

**Figure 11.** Response to Test 2 (open loop, engaged brakes, 2 V step input applied on quadrature **Figure 11.** Response to Test 2 (open loop, engaged brakes, 2 V step input applied on quadrature **Figure 11.** Response to Test 2 (open loop, engaged brakes, 2 V step input applied on quadrature voltage demand).

repetition of tests at different positions of the PMSM rotor provided essentially identical results, thus confirming that the position-dependence of electrical parameters is negligible

repetition of tests at different positions of the PMSM rotor provided essentially identical results, thus confirming that the position-dependence of electrical parameters is negligible

motor phases (i.e., resistance and inductance), and it can be noted that the model succeeds

Figures 10 and 11 are devoted to the identification of the electrical parameters of the

Figures 10 and 11 are devoted to the identification of the electrical parameters of the

Figure 12 is instead relevant for the identification of the magnet flux linkage as well

as for the characterisation of the BEMF waveform with respect to the motor angle. Again, the model succeeds in predicting the hardware response in terms of both magnet flux

**Figure 12.** Response to Test 3 (open loop, open phases, disengaged brakes, motor dragged at 80 rad/s, phase-to-phase BEMF measured). **Figure 12.** Response to Test 3 (open loop, open phases, disengaged brakes, motor dragged at 80 rad/s, phase-to-phase BEMF measured).

**Figure 13.** Response to Test 5 (closed loop on speed and currents, disengaged brakes, square-wave speed demand ±5 rad/sec). FFT analysis is performed on quadrature current signal (proportional to torque) and frequencies are normalized with the steady-state mechanical frequency (0.8 Hz). **Figure 13.** Response to Test 5 (closed loop on speed and currents, disengaged brakes, square-wave speed demand ±5 rad/sec). FFT analysis is performed on quadrature current signal (proportional to torque) and frequencies are normalized with the steady-state mechanical frequency (0.8 Hz).

formances have been verified using the experimentally validated model of the EMA by

reported in Table 2, the EMA position deviation shall be lower than the position sensor accuracy (0.1 deg). To demonstrate the effectiveness of the developed control system with special reference to the application of the SES deformation loop, an extensive simulation

**Figure 14.** EMA dynamic compliance response: simulation results obtained with the experimentally

validated model of the actuator (position setpoint corresponding to minimum stiffness).

campaign has been carried out, and a summary of the results is given in Figure 14.

The performance specification actually requires that under any of the load conditions

characterizing the disturbance rejection capability against external loads.

*3.2. Loads Disturbance Rejection Capability* 

Figure 12 is instead relevant for the identification of the magnet flux linkage as well as for the characterisation of the BEMF waveform with respect to the motor angle. Again, the model succeeds in predicting the hardware response in terms of both magnet flux linkage and BEMF waveform, which is essentially sinusoidal for the reference PMSM (basic assumption of the model).

*Aerospace* **2022**, *9*, x FOR PEER REVIEW 17 of 23

Figure 13 finally reports a closed-loop speed-tracking response, which is essentially relevant for the identification of system inertias (lower left plot in Figure 13), torque disturbance parameters (upper and lower right plots in Figure 13), and friction parameters. It is interesting to note that the FFT analysis performed on quadrature current signal clearly highlights the presence of three harmonic disturbances. The first two harmonic components are multiples of the fundamental electrical frequency, being 10 and 20 times the mechanical frequency (*n<sup>d</sup>* = 10, Table 3), and they derive from small deviations of the BEMF from the sinusoidal waveform. The resting harmonic component is instead located at 24 times the mechanical frequency, and the disturbance can be interpreted as an effect of cogging torque, due to assembly tolerances and/or magnet imperfections. As discussed in [52–55], these irregularities generate torque harmonics at frequencies that are multiple of the stator slots, which in the reference PMSM is 12. **Figure 13.** Response to Test 5 (closed loop on speed and currents, disengaged brakes, square-wave speed demand ±5 rad/sec). FFT analysis is performed on quadrature current signal (proportional to torque) and frequencies are normalized with the steady-state mechanical frequency (0.8 Hz). *3.2. Loads Disturbance Rejection Capability*  Before performing the failure transient analysis (Section 3.3), the control system per-

### *3.2. Loads Disturbance Rejection Capability* formances have been verified using the experimentally validated model of the EMA by

Before performing the failure transient analysis (Section 3.3), the control system performances have been verified using the experimentally validated model of the EMA by characterizing the disturbance rejection capability against external loads. characterizing the disturbance rejection capability against external loads. The performance specification actually requires that under any of the load conditions reported in Table 2, the EMA position deviation shall be lower than the position sensor

The performance specification actually requires that under any of the load conditions reported in Table 2, the EMA position deviation shall be lower than the position sensor accuracy (0.1 deg). To demonstrate the effectiveness of the developed control system with special reference to the application of the SES deformation loop, an extensive simulation campaign has been carried out, and a summary of the results is given in Figure 14. accuracy (0.1 deg). To demonstrate the effectiveness of the developed control system with special reference to the application of the SES deformation loop, an extensive simulation campaign has been carried out, and a summary of the results is given in Figure 14.

**Figure 14.** EMA dynamic compliance response: simulation results obtained with the experimentally validated model of the actuator (position setpoint corresponding to minimum stiffness). **Figure 14.** EMA dynamic compliance response: simulation results obtained with the experimentally validated model of the actuator (position setpoint corresponding to minimum stiffness).

The results clearly highlight that the use of the SES deformation loop implies a relevant enhancement of performances. The compliance response under deterministic loading profiles is marginal without using the SES loop, while it significantly diminishes if the ness);

*3.3. Failure Transient Analysis* 

SES loop is present. Further, the compliance under non-deterministic loads exceeds the requirement limit in the resonant pulsation region (70 Hz) if the SES loop is not used, while it becomes adequate if the SES loop is applied. The experimentally validated model of the EMA has been finally used to characterise the failure transients related to a control hardover fault, as simulated in Equation (26). The following worst-case scenario has been simulated:

The results clearly highlight that the use of the SES deformation loop implies a rele-

vant enhancement of performances. The compliance response under deterministic loading profiles is marginal without using the SES loop, while it significantly diminishes if the SES loop is present. Further, the compliance under non-deterministic loads exceeds the requirement limit in the resonant pulsation region (70 Hz) if the SES loop is not used,

### *3.3. Failure Transient Analysis* • the maximum static load plus the deterministic dynamic loads defined in Table 2 are

*Aerospace* **2022**, *9*, x FOR PEER REVIEW 18 of 23

while it becomes adequate if the SES loop is applied.

The experimentally validated model of the EMA has been finally used to characterise the failure transients related to a control hardover fault, as simulated in Equation (26). The following worst-case scenario has been simulated: applied to the output shaft; • the EMA is demanded to move to the maximum positive deflection (minimum stiff-

	- the EMA is demanded to move to the maximum positive deflection (minimum stiffness); setpoint (*t* = *tFI* = 0 s);
	- the brakes activation occurs with a predefined delay from the fault detection (*tFC* − *tFD* = 51 ms, Umbragroup information). = 51 ms, Umbragroup information). To evaluate the effectiveness of the BEMF damper circuitry, the simulation is per-

To evaluate the effectiveness of the BEMF damper circuitry, the simulation is performed by activating or not the system, obtaining the results in Figures 15 and 16. formed by activating or not the system, obtaining the results in Figures 15 and 16.

**Figure 15.** EMA failure transient with a control hardover fault without BEMF damper circuitry: simulation obtained with the experimentally validated model of the actuator (*tFI* = 0 s ; *tFD* = 13.4 ms; *tFC* = 64.4 ms). **Figure 15.** EMA failure transient with a control hardover fault without BEMF damper circuitry: simulation obtained with the experimentally validated model of the actuator (*tFI* = 0 s; *tFD* = 13.4 ms; *tFC* = 64.4 ms).

64.4 ms).

**4. Discussion** 

**5. Conclusions** 

**Figure 16.** EMA failure transient with a control hardover fault with BEMF damper circuitry: simulation obtained with the experimentally validated model of the actuator (*tFI* = 0 s ; *tFD* = 13.4 ms; *tFC* = **Figure 16.** EMA failure transient with a control hardover fault with BEMF damper circuitry: simulation obtained with the experimentally validated model of the actuator (*tFI* = 0 s; *tFD* = 13.4 ms; *tFC* = 64.4 ms).

It is interesting to note that, though the fault detection latency is extremely small (in both simulations, *tFD* = 13.4 ms), the position deviation without BEMF damper is excessive and the EMA reaches the mechanical endstroke at high speed, possibly implying permanent damages (Figure 15). On the other hand, the use of the BEMF damper permits to It is interesting to note that, though the fault detection latency is extremely small (in both simulations, *tFD* = 13.4 ms), the position deviation without BEMF damper is excessive and the EMA reaches the mechanical endstroke at high speed, possibly implying permanent damages (Figure 15). On the other hand, the use of the BEMF damper permits to strongly limit the position deviation during the failure transient, and the EMA can be blocked by the brakes in safety.

### strongly limit the position deviation during the failure transient, and the EMA can be **4. Discussion**

blocked by the brakes in safety. As highlighted by the results in Section 3.3, the control hardover fault can determine EMA damages and potentially unsafe operation for the FLA movables. The failure transient analysis conducted using the experimentally validated model of the EMA highlights that, in some special operating conditions, even if the fault is detected with extremely small latency (less than 15 ms), the actuator can reach the mechanical endstroke impacting As highlighted by the results in Section 3.3, the control hardover fault can determine EMA damages and potentially unsafe operation for the FLA movables. The failure transient analysis conducted using the experimentally validated model of the EMA highlights that, in some special operating conditions, even if the fault is detected with extremely small latency (less than 15 ms), the actuator can reach the mechanical endstroke impacting at high speed, thus causing permanent damages. To counteract this adverse situation, essentially caused by unavoidable delays in the activation of EMA brakes, the PWR board of the RACER flaps EMAs includes a specifically designed BEMF damper circuitry, which, immediately after the fault detection, opens the power bridge thyristors and connects all the motor phases to the ground, thus generating an electromagnetic damping torque.

The failure transients related to control hardover fault in the EMA employed for the

flap movables of the RACER helicopter-plane are characterised using an experimentally validated model of the system, which includes multi-harmonic torque disturbance simulation, "Coulomb–tanh" friction, mechanical freeplay, and position-dependant stiffness of the ball-screw drivetrain. The failure transient characterization performed in a worst-case scenario in terms of external loads and position setpoint demonstrates that, though the fault detection is executed with extremely small latency (less than 15 ms), a potentially

at high speed, thus causing permanent damages. To counteract this adverse situation, essentially caused by unavoidable delays in the activation of EMA brakes, the PWR board of the RACER flaps EMAs includes a specifically designed BEMF damper circuitry, which, immediately after the fault detection, opens the power bridge thyristors and connects all

## **5. Conclusions**

The failure transients related to control hardover fault in the EMA employed for the flap movables of the RACER helicopter-plane are characterised using an experimentally validated model of the system, which includes multi-harmonic torque disturbance simulation, "Coulomb–tanh" friction, mechanical freeplay, and position-dependant stiffness of the ball-screw drivetrain. The failure transient characterization performed in a worst-case scenario in terms of external loads and position setpoint demonstrates that, though the fault detection is executed with extremely small latency (less than 15 ms), a potentially dangerous actuator runaway can occur, causing high-speed impacts on the EMA mechanical endstroke, caused by the activation delay of the EMA brakes (about 50 ms). Simulation is thus used to point out that an effective solution can be obtained by including a BEMF damper circuitry in the EMA power electronics, which, immediately after the fault detection, opens the power bridge thyristors and connects all the motor phases to the ground, thus generating an electromagnetic damping torque.

The future developments of the research will be focused on:

	- # verify the actual location of the resonant pulsation of the ball-screw drivetrain (currently estimated via FEM analyses);
	- # characterise the actual disturbance rejection of external loads;

**Author Contributions:** Conceptualization, G.D.R.; methodology, G.D.R.; software, B.L.; validation, R.K., M.N. and N.B.; formal analysis, R.K., M.N. and N.B.; investigation, G.D.R. and B.L.; resources, R.K., M.N. and N.B.; data curation, B.L.; writing—original draft preparation, G.D.R.; writing—review and editing, G.D.R. and B.L.; visualization, R.K., M.N. and N.B.; supervision, R.K., M.N. and N.B.; project administration, M.N. and N.B.; funding acquisition, M.N. and N.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

## **References**


## *Article* **Preliminary Sizing of the Electrical Motor and Housing of Electromechanical Actuators Applied on the Primary Flight Control System of Unmanned Helicopters**

**Jeremy Roussel 1,2,\* , Marc Budinger <sup>1</sup> and Laurent Ruet <sup>2</sup>**


**Abstract:** Helicopter dronization is expanding, for example, the VSR700 project. This leads to the integration of electromechanical actuators (EMAs) into the primary flight control system (PFCS). The PFCS is in charge of controlling the helicopter flight over its four axes (roll, pitch, yaw, and vertical). It controls the blade pitch thanks to mechanical kinematics and actuators. For more than 60 years, the actuators have been conventionally using the hydraulic technology. The EMA technology introduction involves the reconsideration of the design practices. Indeed, an EMA is multidisciplinary. Each of its components introduces new design drivers and new inherent technological imperfections (friction, inertia, and losses). This paper presents a methodology to specify and pre-design critical EMAs. The description will be focused on two components: the electrical motor and the housing. This includes a data-driven specification, scaling laws for motor losses estimation, and surrogate modeling for the housing vibratory sizing. The tools are finally applied to two study cases. The first case considers two potential redundant topologies of actuation. The housing sizing shows that one prevails on the other. The second case considers the actuators of helicopter rotors. The electrical motor sizing highlights the importance of designing two separate actuators.

**Keywords:** specification; flight analysis; dimensional analysis; vibration; multidisciplinary optimization

## **1. Introduction**

*1.1. Context*

1.1.1. Helicopter Dronization

Today, we observe a fast increase in the number of projects about OPVs (Optional Pilot vehicles), UASs (Unmanned Aerial Systems), UAVs (Unmanned Aerial Vehicles), UAM (Urban Air Mobility), and VTOL (Vertical Take-Off Landing). For instance, the Ehang 184, the Vahana, the CityAirbus, and Boeing's self-piloted passenger drone can be seen on the civil range; the aerial fighter Northrop Grumman X-47B, the Airbus VSR700 (Figure 1), and the Leonardo AWHero can be seen on the military range. Furthermore, the market for aerial delivery already grows, with vehicles carrying parcels weighing less than 10 kg (DPD France drones in rural areas) up to 100 kg (Kawasaki K-Racer X1 drone prototype). It is clear that flying vehicles and future ones are already required to develop new functionalities to be more autonomous and to be safer. Reducing pilots' working load is part of the current helicopter development roadmaps for enhanced safety. Thus, today's market trend is globally facing a technological watershed toward more electrical solutions. Aircraft makers aim for the More Electric Aircraft (MEA) achievement [1]. The drone concept comes progressively by the implementation of new electrical solutions on already existing vehicles [2].

**Citation:** Roussel, J.; Budinger, M.; Ruet, L. Preliminary Sizing of the Electrical Motor and Housing of Electromechanical Actuators Applied on the Primary Flight Control System of Unmanned Helicopters. *Aerospace* **2022**, *9*, 473. https://doi.org/ 10.3390/aerospace9090473

Academic Editor: Gianpietro Di Rito

Received: 30 June 2022 Accepted: 19 August 2022 Published: 25 August 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

**Figure 1.** VSR700, the multi-mission naval UAS (©Airbus).

1.1.2. Primary Flight Control Systems: From Hydraulic to Electromechanical Technologies

The primary flight control system (PFCS) is in charge of controlling the helicopter flight over its four axes (roll, pitch, yaw, and vertical) by the control of the blade attack angles. The paper [3] describes the PFCS with pictures. More details can be found in the book [4] and the handbook [5]. The automatic pilot function is ensured through EMAs located in series and in parallel with the mechanical kinematic (Figure 2a).

**Figure 2.** PFCS architectures: sketch of principle. From (**a**–**c**), this figure shows basically the evolution of PFCS architectures as it applies to autonomous helicopters. This evolution clearly shows that the helicopter mass can be significantly reduced with the reduction in part numbers. As long as the actuation control loop is concerned, this part number reduction decreases the response delay sources. The quid pro quo for it is the increase in the actuator critical level because it gathers nearly all piloting functions.

Only one helicopter in the Airbus fleet (NH90) is fly-by-wire. The hydraulic actuators (Direct Drive Valve, DDV) are commanded directly by four electrical torque-motors connected to the FCC (Flight Control Computer), as shown in Figure 2b.

The hydraulic technology has been conventionally used in actuators for more than 60 years [1,6]. A new trend uses EMAs as substitutes for hydraulic actuators in the PFCS in actual helicopters or as part of fly-by-wire PFCSs in new autonomous helicopters (Figure 2c). This requires the reconsideration of the design practices right at the preliminary design phase. The VSR700 (Figure 1) is a use case. It is an already proved light helicopter (Cabri G2) turned into a drone by the integration of electrical components. These components include four EMAs in the PFCS.

## 1.1.3. Business Need

Today's business need is to perform EMA preliminary studies and learn more about this technology. Indeed, further to an application in the primary flight control system (main rotor and tail rotor), the EMA technology finds interest in many other helicopter applications, such as the automatic pilot actuators or the secondary flight control system (SFCS). This includes the actuation of landing gears and wing flaps and the propellers' pitch control (for high-speed helicopters, such as the RACER helicopter).

## *1.2. Helicopter Specificities*

The helicopter PFCS is a specific application for actuators. Linked to safety-critical surfaces, the actuator failure is qualified as "catastrophic" as it leads to the helicopter's crash. The actuation unit must comply with fail safe characteristics. Therefore, the actuation unit must have a redundant topology by force or position summing. The loading spectrum coming from the rotor spinning contains high-frequency content (fundamental at 20 hz minimum). Moreover, the operational pilot demand scenarios are difficult to predict. The required actuation performance is demanding as it combines a significant bandwidth response, reduced plays, high stiffness, and a low mass. The on-board environment is severe with temperature variations within [−40; 85] ◦C, humidity, a salty atmosphere, and vibrations (6 to 20 g, [7]). Finally, concerning the component lifespan, the time between overhauls is required to be from 3000 flight h onwards. This set of specificities highlights the critical function occupied by the actuator in terms of aviation safety.

## *1.3. Electromechanical Actuators*

An electromechanical actuator (EMA) includes components from multiple disciplines. There are mechanical parts (rod ends, bearings and rotary/linear conversion mechanisms, and a clutch), electrical parts (a motor, brake, and clutch), and electronic parts for power and control. The paper [3] presents an example of an EMA in detail and the different possible EMA architectures.

Except for low-power and/or less safety-critical applications (flaps, slats, spoilers, and a trim horizontal stabilizer), the EMA technology is not mature enough yet for primary flight controls [6]. This is essentially because of their lack of accumulated return of experience. The statistical database on components fault modes is poor [8]. EMAs entail some concerns in terms of reliability, risk of failures due to the jamming in the mechanical transmission components, health monitoring (HM) and assessment, and thermal management. The EMA applicability in aerospace has been proved in terms of dynamic performances [8]. In addition, EMAs offer interesting perspectives in terms of maintenance, integration, reconfiguration in case of failure, ease of operation, harsh running environment ([−50, 125] ◦C), and management of power [1,6].

## *1.4. Preliminary Sizing Method*

The preliminary design methodologies can be divided into two phases. The first phase is the system architecture choice, commonly guided by reliability studies, such as those presented in [9–11]. At this level, there are difficulties in taking into account the entire set of design criteria, important to evaluate an architecture. A study with a higher level of details is necessary, especially in the case of the EMA as it includes many constraints and interdisciplinary couplings between components. This is up to the second phase: the preliminary sizing. This phase is mainly based on multidisciplinary optimizations. The models used are usually analytical models or response surface models (RSM) to facilitate the design space exploration and the design optimization.

Some already existing preliminary sizing methodologies can be cited. The references [12,13] present a methodology to select the motor and gearhead of the actuators in the automotive field. The methodology includes a selection based on scaling laws. It outputs graphs showing all feasible motor/gear ratio combinations. In the aeronautic field and regarding the secondary flight control actuators, the paper [14] presents a methodology for the preliminary design of mechanical transmission systems. It is formalized as a constraint satisfaction problem (CSP) with an automated consistency checking and a pruning of the solution space. The mechanical components are modeled by scaling laws. In addition, the paper [10] presents a simulation and an optimization strategy to evaluate two concepts of actuation systems: the conventional hydraulic actuators and the electro hydrostatic actuators (EHAs). Moreover, the paper [15] describes a preliminary design method of EHAs based on multi-objective optimization (MOO) with Pareto dominance. Two objectives are set: the minimization of the mass and the maximization of the efficiency. The weight prediction is achieved using scaling laws, and the efficiency is calculated by a static energy loss model. The method outputs the design parameter, leading to a Pareto front in mass and efficiency.

The scaling laws are good candidates for preliminary design. They are an example of analytical models. They are interesting because they require few data to build them and validate them (existing industrial product ranges). In aerospace, the scaling laws are broadly used in the conceptual design of aircraft, especially regarding aerodynamics, propulsion, structure, and mass [16]. Moreover, the propellers are often described by scaling laws [17–19]. Furthermore, the scaling laws can be used in the field of robotic actuators where low speeds and high torques are usually required [20].

This paper suggests a preliminary sizing methodology applied in the aeronautic field, on helicopters, and on critical EMAs. It contains similar concepts as found in the previously cited literature. Indeed, it includes scaling laws. It includes an optimization where the component selection must satisfy a specification and design constraints while minimizing the actuator total mass. This paper adds its value by the introduction of two elements not considered in the literature yet. The first element is the motor heating based on dynamic criteria extracted from an equivalent representation applied on complex real mission profiles. The helicopter application is more dynamic than the aircraft application, as highlighted in the paper [3]. The second element is the vibratory environment through the actuator housing sizing.

## *1.5. Objective and Outline*

The objective is to set up a design methodology supported by tools estimating the actuator component characteristics for any helicopter PFCS application. The considered EMA architecture is a direct drive in-line EMA. Details regarding this architecture are provided in Section 7.

To address this topic, firstly, this paper briefly presents the proposed global methodology. Secondly, it briefly recalls the tools established in [3] to better understand the actuation need and develops the actuator specification from the given flight data records. Thirdly, the modeling is detailed, limited to the electrical motor and the housing. Finally, the toolbox is applied to two real use cases. One compares two potential redundant topologies of actuation in terms of housing. The other compares the electrical motor characteristics obtained for the specifications of a given helicopter (main rotor and tail rotor).

## **2. Global Methodology**

To answer the preliminary design need, this paper offers the global methodology presented in Figure 3. It consists of three main steps:


This methodology takes flight data and high-level project requirements as inputs. It outputs the characteristics of each component necessary to answer the actuation need with a minimized mass.

**Figure 3.** Proposed methodology for actuator sizing.

This paper addresses step 1 based on already published work in [3]. It details step 2 focusing only on the tools set up to model two components: the electrical motor and the housing. Step 3 is not explained. The mass minimization is a multidisciplinary design optimization (MDO) problem with numerous interdisciplinary couplings [21–23]. The MDO architecture is chosen to be with a single optimizer. Some variables and constraints are added into the sequencing. This optimization settings correspond to a hybrid individual disciplinary feasible (IDF) architecture and were inspired from [23–25].

## **3. Mission Profile Analysis**

To reach the "first time right" objective of a design office, the requirements of any application must be well understood. In that way, the specification is fully representative, and the design is well guided. The EMA application in a PFCS of a helicopter combines three types of difficulties which impede easily writing down its specification. These difficulties are: a set of multidisciplinary design drivers due to the different technologies of the components; an external loading spectrum, coming from the rotor spinning in air, difficult to model; and operational piloting scenarios difficult to predict.

At the time of the project development of unmanned vehicles, a reduced number of short flight-test records with an interesting sampling rate can be available. These data are called mission profiles. They are the most trustful inputs representative of the application requirements. The paper [3] offers to set up an EMA specification based on the analysis of these mission profiles. The offered methodology is briefly recalled in the following subsection.

## *3.1. Methodology*

The methodology offered by [3] has the objective to build the specification whilst simplifying the data analysis for the engineer and keeping them as the decision maker. The methodology is summed up in Figure 4. It is inspired from the work developed on railway trains and aircraft by [26,27]. It contains 4 main steps:


**Figure 4.** Proposed methodology for specification setup.

## *3.2. Estimation over Real Mission Profiles*

The paper [3] illustrates the specification methodology, applying it on real flight data, and emphasizes the specificity of a rotorcraft PFCS for EMA application. It compares the specification of a main rotor actuator (MRA) to the specification of a tail rotor actuator (TRA). The MRA and TRA integration is illustrated in Figure 5. The analysis of the data coming from a real helicopter flight (VSR700, Figure 1) over three flight phases (take-off, cruise, and landing) provides the specifications of the MRA and TRA. In the paper [3], both specifications are compared using ratios elaborated for comparisons. In this paper, both specifications are presented in Table 1. For confidentiality reasons, the values cannot be displayed. Therefore, Table 1 presents ratios such as *XMRA*/*XTRA*, where *X* is the concerned specification indicator. These specifications will be the inputs in Section 7.3.

**Figure 5.** Sketch of principle of the main rotor (**left**) and the tail rotor (**right**) on helicopters. The main rotor is responsible for lifting the helicopter weight. The tail rotor ensures the control of the helicopter yaw axis.

At first glance, Table 1 clearly shows that the MRA faces loads roughly 10 times higher than the TRA and moves with an equivalent continuous speed and acceleration roughly 40% lower than the TRA. In the paper [3], the comparison between both applications is emphasized regarding the continuous and maximum loading (from the rotor inertia and external load), the rolling fatigue of the mechanical components, the rotor inertia impact compared to the external load, and the motor loss identification in terms of the continuous motor torque contribution. This paper adds an additional indicator regarding the load-pitting analysis.

**Table 1.** VSR700 specification for MRA and TRA. *a* is the acceleration, *F* is the load; they come from measures on the VSR700 helicopter PFCS. *PR* = *F*(*t*) · *a*(*t*) is the power rate. The mean power rate *PR* = | < *F*(*t*) · *a*(*t*) > | takes an absolute value to be conservative. *viron* is a mean value of speed representative of motor iron losses. *FRMC* is an equivalent rolling fatigue load based on 10<sup>6</sup> cycles. (details in [3]).


The pitting load refers to the fluctuating load when the actuator is motionless. This paper suggests estimating the averaged contribution of the load upper frequencies. This can be performed compared to the entire load frequency content. Therefore, the load is separated from its low frequencies using a filter with a cut-off frequency set at the actuator bandwidth value. This resulting load is averaged by the RMS to estimate an equivalent continuous load. This result is normalized by the equivalent continuous load, including all load frequency content. The indicator is *Rp*:

$$R\_p = \frac{\left[hutter\_{HP}[F(t)]\right]\_{rms}}{[F(t)]\_{rms}}\tag{1}$$

The indicator *R<sup>p</sup>* is evaluated on a nontransient phase of the VSR700 mission profile (Table 2). The indicator shows a higher level of the pitting load in the MRA. This result is understood considering the helicopter rotor architectures (Figure 5). As a piece of information, the MRA selectively controls the blade attack angle during one rotor azimuth. Meanwhile, the TRA simultaneously controls the attack angle of all the blades whatever the azimuth.

**Table 2.** Comparison of the indicator *Rp* for MRA and TRA applications on VSR700 mission profiles. The indicator is evaluated over steady state flight phases (cruise).


## **4. Models for the Preliminary Sizing of EMA**

*4.1. Modeling Overview*

The EMA includes multidisciplinary components. Each of them has multiple and different key design drivers and operational scenarios. The selected components should comply with the actuator specification and ensure a minimized total mass. To answer the need of component selection, a knowledge-based process around component modeling must be set. It covers two levels [28]:

1. The component level: It deals with the determination of component characteristics from a reduced number of parameters to facilitate the optimization. The models involved in it are called the *estimation models* (Figure 6).

This paper proposes to focus only on the *estimation models*. To cover the tools setup, only two components are presented. They are those at stake within two onboard specificities: the thermal and vibratory environments. Therefore, the following sections deal with the modeling proposed for the electrical motor and the housing characteristics.

2. The system level: It deals with interactions between components, operational scenar-

## *4.2. Need of Estimation Models and Approach*

As shown in Figure 6, to start a first iteration in the sizing loop, the main characteristics of each component have to be identified from a reduced set of key parameters. The *estimation models* are introduced for this purpose. Per component, they directly link the primary characteristics, which define the component functionally, to the secondary characteristics, which can be seen as the dimensions and features of the imperfections. Thus, the *estimation models* provide the necessary parameters for the integration study, *simulation models*, and *evaluation models*.

Generally, at the component level, the models link the physical dimensions and characteristics of in-use materials to the primary and secondary characteristics. The design at the component level is an inverse problem which requires the primary characteristics as inputs.

In such a context of multidisciplinary modeling with optimization, a unified modeling approach is required. A dimensional analysis and the Vaschy–Buckingham theorem (Theorem 1) are good candidates for it [29,30]. Indeed, they are extensively used in aerodynamics and fluid mechanics because they provide a more physical and unified framework with a reduced number of parameters. This section shows how they can be extended to other domains, such as the electrical motor and the housing of an EMA.

**Theorem 1** (Vaschy–Buckingham theorem)**.** *Any physical equation dealing with n physical variables depending on k fundamental units (mass, length, time, temperature, charge) can be formulated as an equivalent equation with p* = *n* − *k dimensionless variables called "π-numbers" built from the initial variables.*

The development steps of an estimation model are presented in Figure 7. The starting point is the expression of one component characteristic *y* as an algebraic function *f* depending on the geometrical dimensions and material/physical properties *p<sup>i</sup>* . *L* is a length and *d<sup>i</sup>* the rest of dimensions.

$$y = f(\underbrace{\text{L}, d\_1, d\_2, \dots}\_{\text{dimensions}}, \underbrace{p\_1, p\_2, \dots}\_{\text{properties}}) \tag{2}$$

**Figure 7.** Estimation model development steps.

Applying the dimensional analysis and Theorem 1, the problem is rewritten into a reduced number of dimensionless parameters (Equation (3), [31]).

$$\pi\_{\mathcal{Y}} = \operatorname{g}(\underbrace{\pi\_{d\_1}, \pi\_{d\_2}, \dots, \pi\_{p\_1}, \pi\_{p\_2}, \dots}\_{p \text{ variables}}) \text{ with } \begin{cases} \pi\_{\mathcal{Y}} = \mathcal{Y}^{c\_{\mathcal{Y}}} \cdot L^{c\_{\mathcal{L}}} \cdot \prod\_{i} p\_i^{c\_{p\_i}} \\ \pi\_{d\_i} = \frac{d\_i}{L} \\ \pi\_{p\_i} = L^{c\_{p\_i 0}} \cdot \prod\_{j} p\_i^{c\_{p\_j j}} \end{cases} \tag{3}$$

The next step is to develop the estimation models based on the *π*-groups.

Scaling law formulations are undertaken when the dimensionless numbers *πd<sup>i</sup>* and *πp<sup>i</sup>* remain constant around a given component product range. This means the geometrical and/or material similarities are satisfied. This is applicable for the electrical motor of the actuator as detailed in the following Section 5.

When the dimensionless numbers *πd<sup>i</sup>* and *πp<sup>i</sup>* are not considered as constant, the approximation of the function *g* can be achieved by performing data regressions [31,32]. The data may come from manufacturer product data, test measurements, or finite element simulation results based on design of Experiment (DoE) as presented in Section 6 for the housing vibratory sizing.

## **5. Scaling Laws**

In this section, the scaling law theory is firstly developed and then applied to an electrical motor. The actuator sizing requires *estimation models* for its integration (motor dimensions) and its losses (copper and iron losses, inertia).

## *5.1. Fundamentals*

In the literature, scaling laws are also called *similarity laws* or *allometric models* [33]. They estimate the component main characteristics requested for their selection without requiring a detailed design.

Scaling laws are based on three hypotheses:


Once these assumptions are satisfied, *π<sup>y</sup>* is stated to be constant because it depends on constant variables:

$$
\pi\_y = \lg(\pi\_{d\_1}, \pi\_{d\_2}, \dots, \pi\_{p\_1}, \pi\_{p\_2}, \dots) = \text{constant} \tag{4}
$$

This gives the standard power-law shape of a scaling law:

$$
\pi\_y = y^{c\_y} \cdot L^{c\_L} \cdot \prod\_i p\_i^{c\_{p\_i}} = \text{constant} \implies y \approx L^c \tag{5}
$$

with *c* a constant. Then, as proposed by [34], the "star" notation is introduced. It indicates the scaling ratio of a desired component characteristic *x* by the same characteristic *xre f* of a component taken as a reference: *x* <sup>∗</sup> = *x*/*xre f* . This component of reference is picked up into the supplier range of the considered product.

Thus, Equation (5) becomes:

$$y^\* = L^{\*c} \iff \frac{y}{y\_{ref}} = \left(\frac{L}{L\_{ref}}\right)^c \tag{6}$$

From a single component of reference and a reduced number of parameters (no detailed design required), the scaling laws quickly extrapolate the main characteristics *yre f* of a known component toward the characteristic *y* of a possible component of the same technology:

$$y = y\_{ref} \cdot \left(\frac{L}{L\_{ref}}\right)^c \tag{7}$$

Consequently, the scaling laws level down the complexity of the inversion problem. All the useful relations are easily expressed as a function of a single *key design parameter* (also named *definition parameter*) that is associated with the component under design (Figure 6).

## *5.2. Electrical Motor Scaling Laws*

The previously mentioned approach is applied to the electrical brushless motor. The motor mass-law formulation is detailed hereafter. Beforehand, some hypotheses need to be stated: the main design driver is the maximum continuous winding temperature; the natural convection is the dominant thermal phenomenon; the mean induction in the airgap is constant for a given magnet technology; the number of magnets is constant over the considered product range; and the geometric similarities are verified, and the material and boundary limits similarities are satisfied.

In Figure 8, following the approach mentioned in Section 5.1, step by step, the torque evolution is obtained.

As diameters and lengths are supposed to evolve similarly (*d* ∗ = *L* ∗ , Figure 9), the mass *M* of the motor is basically approximated by:

$$M = \int \rho\_{\text{eq}} \, dV \implies M^\* = L^{\*3} \tag{8}$$

Using the torque expression (prerequisite of Figure 8), the motor mass *M* becomes:

$$M^\* = T^{\*3/3.5} \implies M = M\_{ref} \cdot \left(\frac{T}{T\_{ref}}\right)^{3/3.5} \tag{9}$$

Figure 10 compares the evolution of this law (Equation (9)) to real data from two manufactured motors: Parvex NK [35] and Kollmorgen RBE [36]. It is observed that a single reference of a motor allows to rebuild the evolution of the motor mass for a broad range of torque. This is possible with less than 10% of the mean relative error *e*.

**Figure 9.** Electrical motor: homothety hypothesis.

**Figure 10.** Motor mass: scaling law prediction compared to manufacturer catalogs (Parvex NX, Kollmorgen RBE).

In the same way, and following the same hypothesis, the other motor characteristic laws are formulated. Table 3 presents them with their prediction level compared to the Parvex NK catalog. The prediction levels compared to the RBE catalog are presented when data are missing in the Parvex catalog. Both catalogs validate the presented laws.


**Table 3.** Electrical motor scaling law sum up and prediction quality compared to manufacturer catalog data.

(1) Parvex specificity; (2) definition based on magnetic saturation criteria, missing data for validation; (3) definition based on thermal dissipation criteria; be careful that RBE catalog gathers different thermal integrations, *T* is dependent on them; (4) based on mechanical limitations; <sup>∗</sup> the 'star' notation: *x* <sup>∗</sup> = *x*/*xre f* .

In Table 3, the motor constant *K<sup>m</sup>* [(*N* · *m*) <sup>2</sup>/*W*] is introduced. It is the ratio of the squared torque provided per unit of heat generated. Moreover, it can be found in catalogs with another definition: [(*<sup>N</sup>* · *<sup>m</sup>*)/*W*0.5]. *<sup>K</sup><sup>m</sup>* is an interesting parameter because it directly links the application mechanical need (the required motor torque) with the motor losses without knowing the motor winding characteristics.

As far as the motor losses are concerned, [34,37] mention that the copper and iron losses bound the motor continuous operation domain. The operation domains of electrical motors can be found in [3,33]. Thus, at a steady state, the total heat generated by the electrical motor is the sum of Joules and iron losses, such as:

$$Q\_{th} = Q\_{Joules} + Q\_{ion} = \alpha \cdot T^2 + \beta \cdot f\_{elect}^b \tag{10}$$

where *α* and *β* are, respectively, copper and iron losses coefficients, and *felect* the electrical speed. *b* is a constant depending on the hysteresis and the Eddy's current contributions into iron losses, and [33,37] indicate a mean value of 1.5.

As far as the other components of the actuator are concerned, their scaling laws are found in the publications [20,33].

## **6. Regression Models**

In a helicopter context, the resonance frequencies and stress under a vibratory environment is an unavoidable check to perform. This section presents a preliminary vibratory study of a simplified actuator housing. The approach goes through the synthesis of a regression model (also called the surrogate model) to implement into the actuator preliminary sizing loop.

## *6.1. Introduction*

From a purely mechanical point of view, the design of the EMA housing has to focus on the elementary forces acting on the housing, which can be divided into two categories: the static stresses induced by the power transmission to the load, which have low frequencies, and the vibratory stresses induced by the vibratory environment, which have high frequencies (*f* ∈ [5; 2000] hz, [7]).

The path of the various static or dynamic loads is represented for a generic actuator in Figure 11.

The most significant static loadings are the tensile, compressive, buckling forces. They are transmitted through the rod to the nut, the screw, the bearings, and finally to the housing. The high number of cycles generally requires the fatigue limit of materials to be considered.

The dynamic stress is mainly generated by the transversal vibrations due to the vibratory environment which can generate important mechanical bending stresses. For a long actuator, as it can be for a direct-drive EMA, these stresses are prevailing.

## *6.2. Prior Considerations*

A high-fidelity model of vibratory stress in the housing would be difficult to develop. Indeed, ball bearings and roller screw stiffnesses and plays are unknown and not supplied in datasheets. The contact at the linear bushing level is unclear.

Therefore, we propose a simple model based on some simplifications. The first one concerns the potential mechanical backlash in the actuator assembly and their effect with respect to excitation frequencies. The vibratory amplitude relies on the acceleration level, as illustrated in Table 4.


**Table 4.** Vibratory amplitudes for different acceleration and frequency.


The vibratory amplitudes are estimated in Table 4 with Equation (11) regarding commonly used accelerations found in the DO160 standard [7] (6 g for in-cabin equipment, 20 g for an under-swashplate location). The amplitudes can be lower and close to the typical plays. In this case, the vibratory phenomena becomes even more complex. Typical plays at the linear bushing level are roughly 100 µm. The linear bushing or sleeve bearing guides the output rod inside the housing.

To model the effect of the play, we propose a lumped-parameter simulation. It associates one or two mass–spring system(s) excited through an elastogap where the play is modeled as a non-linear stiffness (very low value into the play, and high value far from the play).

In Figure 12, the accelerations stated by the DO–160 standard [7] are plotted in terms of amplitude and frequency.

**Figure 12.** DO–160 accelerations: amplitude evolutions with respect to frequencies.

Two cases are studied and illustrated thereafter: *case 1* at 250 hz with *x* ≤ 100 *µ*m and *case 2* at 50 hz with *x* ≥ 100 *µ*m. Both cases consider an acceleration of 6 g.

First of all, the model considers only one mass–spring system.

*Case 1* is simulated in Figure 13. The mass–spring resonance is set at 250 hz and the excitation is around this frequency. No resonance mode of the mass is observed, and the vibratory excitation is 'filtered' by the play.

**Figure 13.** Mass amplitudes of the simple mass–spring coupled with play (*f<sup>r</sup>* = 250 hz).

*Case 2* is simulated in Figure 14. The mass–spring resonance is set at 50 hz and the excitation is around this frequency. The mass vibrations are more important, and the resonance is observed.

Now, for the plays inside the actuator, the model considers two mass–spring systems linked by an elastogap. Using this model, the displacement of each mass is plotted in Figure 15. For a resonance frequency around 50 hz, the amplitudes are such that the vibrating parts interfaced by the play can be considered as one and a single part.

**Figure 14.** Mass amplitudes of the simple mass–spring coupled with play (*f<sup>r</sup>* = 50 hz).

**Figure 15.** Mass amplitudes of the double mass–spring coupled with play (*f<sup>r</sup>* = 50 hz).

For a resonance frequency around 250 hz, the vibratory amplitudes are much lower, and the two masses evolve within the play.

As a result, to keep the vibratory amplitudes smaller than the play between the parts, the resonance frequencies must be high (e.g., ≥200 hz; this limit depends on the amplitude of the acceleration). It is not the case for long and narrow actuators as direct-drive actuators can be. The resonance frequencies are low.

Thus, some simplifications can be introduced into the estimations of stresses of the actuator envelop under vibratory excitations. We propose to suppose, for housings with low resonance frequency, that the plays are negligible compared to the vibratory amplitudes. Consequently, the contact with linear bushing is modeled as an infinitely rigid contact.

## *6.3. Hypothesis*

Now, an FEM model of reduced parameters is developed in order to represent it by a surrogate model. Some more hypotheses are formulated.

A simplified geometry is considered (Figure 16): two hollowed cylinders, one in the other. The connection between them is assumed to be perfect. The set (motor, brake, connectors, and bearings) is supposed to be a cylinder of 1/3rd of *L<sup>a</sup>* with an equivalent density. This equivalent mass is modeled with Young's low modulus (1/10th of aluminium modulus). This choice is conservative; it is so as to not impact the stiffness of the housing.

**Figure 16.** Actuator simplified geometry for modal analysis.

The housing and output rod are set with the same material properties (aluminium). Their lengths are, respectively, supposed to be 2/3rd and 1/3rd of *La*.

The nut is modeled as a full cylinder with 90% steel density (7000 kg·m−<sup>3</sup> ) to consider the air content between the rollers. The geometry of the nut evolves with the geometrical similarity assumption (scaling law). The cylinder representing the nut is modeled using a low Young's modulus (1/10th of steel modulus) so as to not influence the stiffness of the structure.

The rod ends allow rotation with no friction.

The antirotation key and the sealing leap at the interface output rod with housing are not modeled.

The three following sections develop a surrogate model using the surrogate modeling technic suggested in the paper [31].

## *6.4. Problem Formulation*

Under vibration, the system can be associated to a basic damped mass–spring model. Figure 17 presents this model with *U* (m) the displacement of an equivalent mass *Meq* (kg) evolving according to an excitation force *<sup>F</sup>* (N), a stiffness *<sup>K</sup>eq* (N · <sup>m</sup>−<sup>1</sup> ), and a damping *<sup>C</sup>eq* (N · <sup>m</sup>−<sup>1</sup> · s).

The stress is linearly linked to the displacement:

$$
\sigma = k\_{\sigma} \cdot \mathcal{U} \tag{12}
$$

with *σ*<sup>0</sup> [*Pa*] and *U*<sup>0</sup> [m], the maximum stress and the maximum displacement at resonance frequency.

**Figure 17.** Mass–spring model.

Newton's second law applied to the moving body enables the Laplace function between the displacement *U*(*t*) and the excitation load *F*(*t*) to be estimated, such as:

$$\frac{\mathcal{U}(p)}{F(p)} = \frac{1}{M\_{\text{eq}} \cdot p^2 + \mathcal{C}\_{\text{eq}} \cdot p + \mathcal{K}\_{\text{eq}}} = \frac{1/\mathcal{K}\_{\text{eq}}}{p^2/\omega\_r^2 + 2 \cdot \tilde{\xi} \cdot p/\omega\_r + 1} \tag{13}$$

Considering an excitation of the mass with a sinusoidal force *F*(*t*) = *F*<sup>0</sup> · *sin*(*ω* · *t*), the maximum displacement at the first resonance mode is:

$$\mathrm{d}\mathcal{U}\_{0} = \frac{1}{\mathbf{2} \cdot \mathfrak{F}} \cdot \frac{\mathrm{F}\_{0}}{\mathrm{K}\_{\mathrm{eq}}} = \frac{\mathrm{Q}\_{\mathrm{m}} \cdot \mathrm{F}\_{0}}{\mathrm{K}\_{\mathrm{eq}}} \qquad\qquad\qquad\qquad\qquad\qquad\zeta = \frac{\mathrm{C}\_{\mathrm{eq}}}{\mathbf{2} \cdot \sqrt{\mathrm{K}\_{\mathrm{eq}} \cdot \mathrm{M}\_{\mathrm{eq}}}} \approx \frac{1}{\mathbf{2} \cdot \mathrm{Q}\_{\mathrm{m}}}\tag{14}$$

where *Q<sup>m</sup>* is the mechanical quality coefficient.

The article [38] reports that tests performed on industrial prototypes show a wide range of practical values for the equivalent mechanical quality coefficient *Qm*. In addition, it reports that experiments give typical values for *Q<sup>m</sup>* between 10 and 50, depending on the application and boundary conditions. For structural dynamic models, in the absence of better information, it is normally acceptable to assume a value of *Q<sup>m</sup>* = 30 (according to [39]).

The equivalent force *F*0of the acceleration effect can be evaluated thanks to an equivalent work W<sup>0</sup> [40,41]:

$$\mathcal{W}\_0 = \mathbf{F}\_0 \cdot \mathbf{U}\_0 = \iiint\_{\mathcal{V}} \boldsymbol{\mu}\_0 \cdot \boldsymbol{a} \cdot \boldsymbol{\rho} \cdot d\mathcal{V} \tag{15}$$

with *u*0(*x*) the deflection of the actuator envelop, *a* the amplitude of the vibratory sinusoidal acceleration. A mass *Macc* subjected to the acceleration can be defined:

$$F\_0 = M\_{\rm acc} \cdot a \tag{16} \\ \text{or} \\ \qquad \qquad M\_{\rm acc} = \frac{\iint \int\_{\mathcal{V}} u\_0 \cdot \rho \cdot d\mathcal{V}}{U\_0} \tag{16}$$

The mass subjected to the acceleration is not identical to the mass expressing the kinetic energy, *Meq*, defined such as [40,41]:

$$\frac{1}{2} \cdot M\_{\text{eq}} \cdot V\_0^2 = \iiint\_{\mathcal{V}} \frac{1}{2} \cdot \rho \cdot v\_0^2 \cdot d\mathcal{V} \tag{17}$$

with, at the first mode resonance, *V*<sup>0</sup> the speed of *Meq* and *v*0(*x*) the speed of each point of the actuator deflection. The speeds can be defined such as:

$$
v = w\_0 \cdot u \qquad\qquad\qquad V = w\_0 \cdot \mathcal{U}\tag{18}$$

Thus, we can easily define the equivalent mass such as:

$$M\_{eq} = \frac{\iint \int\_{\mathcal{V}} \rho \cdot u\_0^2 \cdot d\mathcal{V}}{\mathcal{U}\_0^2} \tag{19}$$

The following ratio is introduced:

$$k\_{\rm acc} = \frac{M\_{\rm acc}}{M\_{eq}} = \begin{cases} \mathcal{U}\_0 \cdot \frac{\iiint\_{\mathcal{V}} \, u\_0 \cdot d\mathcal{V}}{\iiint\_{\mathcal{V}} \, u\_0^2 \, d\mathcal{V}} & \text{if } \rho \text{ is constant} \\ \mathcal{U}\_0 \cdot \frac{\iiint\_{\mathcal{V}} \, u\_0 \cdot \rho \cdot d\mathcal{V}}{\iiint\_{\mathcal{V}} \, u\_0^2 \, \rho \cdot d\mathcal{V}} & \text{if } \rho \text{ is not constant} \end{cases} \tag{20}$$

and the maximum displacement at the resonance can be approximated as it follows:

$$\mathcal{U}\_0 = \frac{\mathcal{Q}\_m \cdot k\_{\rm acc} \cdot \mathcal{M}\_{\rm eq} \cdot a\_0}{\mathcal{K}\_{\rm eq}} = \frac{\mathcal{Q}\_m \cdot k\_{\rm acc} \cdot a\_0}{\omega\_0^2} \tag{21}$$

with *ω*<sup>2</sup> <sup>0</sup> = *Keq*/*Meq* the resonance angular frequency.

## *6.5. Dimensional Analysis*

As seen in Section 5, the use of a dimensional analysis and Buckingham's Theorem enable to reduce the number of variables expressing a physical problem. Here below, this approach is developed for the vibratory use case. By simplification and for a reduced number of parameters, a constant density *ρ* is assumed all along the actuator (Equation (20)).

The link between stress and displacement evolves according to the following variables:

$$\frac{\sigma}{\overline{\omega}I} = k\_{\sigma} = f(E, d\_{\rm rs}, L\_{a\prime}e\_1, e\_2, L\_{\rm rs}) \tag{22}$$

which can be rewritten with the following dimensionless numbers:

$$
\pi\_{k\_{\sigma}} = \frac{\sigma \cdot d\_{\rm rs}}{U \cdot E} = f(\frac{L\_a}{d\_{\rm rs}}, \frac{e\_1}{d\_{\rm rs}}, \frac{e\_2}{d\_{\rm rs}}, \frac{L\_{\rm rs}}{d\_{\rm rs}}) \tag{23}
$$

The resonance angular frequency evolves according to:

$$
\omega\_0 = \lg(E, \rho, d\_{rs}, L\_a, e\_1, e\_2, L\_{rs}) \tag{24}
$$

which can be rewritten with the following dimensionless numbers:

$$
\pi\_{\omega\_0} = \omega\_0 \cdot \left(\frac{\rho}{E}\right)^{1/2} \cdot d\_{\rm rs} = g(\frac{L\_a}{d\_{\rm rs}}, \frac{e\_1}{d\_{\rm rs}}, \frac{e\_2}{d\_{\rm rs}}, \frac{L\_{\rm rs}}{d\_{\rm rs}}) \tag{25}
$$

The stress under a vibratory acceleration can be expressed as:

$$
\sigma = k\_{\sigma} \cdot \mathcal{U} = k\_{\sigma} \cdot \frac{Q\_{m} \cdot k\_{\text{acc}} \cdot a}{\omega\_{r}^{2}} = \sigma\_{0} \cdot Q\_{m} \cdot a \cdot \frac{\iiint\_{\mathcal{V}} \,\mu\_{0} \cdot d\mathcal{V}}{\iiint\_{\mathcal{V}} \,\mu\_{0}^{2} \cdot d\mathcal{V}} \tag{26}
$$

The stress evolves according to:

$$
\sigma = h(k\_{\sigma}, \omega\_{0\prime}^2 a\_{\prime} Q\_{\mathfrak{m}}) \tag{27}
$$

which can be rewritten as:

$$
\pi\_0 = \frac{\sigma}{Q\_m \cdot a \cdot d\_{rs} \cdot \rho} = h(\frac{L\_a}{d\_{rs}}, \frac{e\_1}{d\_{rs}}, \frac{e\_2}{d\_{rs}}, \frac{L\_{rs}}{d\_{rs}}) \tag{28}
$$

The expression of the stress is thus only function of four aspect ratios. One of these aspect ratios, *Lrs*/*drs*, can be assumed to be constant because of the geometrical similarity assumption used for roller screw component sizing (scaling law).

The final expression of the stress is a function dependent of three dimensionless quantities:

$$\pi\_0 = \operatorname{g}(\pi\_1, \pi\_2, \pi\_3) \begin{cases} \pi\_0 = \sigma / (\operatorname{Q}\_{\text{m}} \cdot a \cdot d\_{\text{rs}} \cdot \rho) \\ \pi\_1 = \operatorname{L}\_a / d\_{\text{rs}} \\ \pi\_2 = e\_1 / d\_{\text{rs}} \\ \pi\_3 = e\_2 / d\_{\text{rs}} \end{cases} \tag{29}$$

It remains to determine this function *g*. The following section does the job.

## *6.6. FEM Software Model*

Using a software for computation by finite elements, a model is parametrized according to the previous considerations and hypotheses. It enables to obtain:


The intersection of both cylinders has been taken care of by smooth and arced geometries to avoid numerical stress constraints. The boundaries are pinned at each extremity of the actuator model. The deflection is allowed within the plane of the section presented in Figure 18.

**Figure 18.** Modal analysis of the simplified actuator geometry.

The tandem topology of actuation is largely used in the aeronautic field because it complies with safety aviation rules. This topology involves two actuators stuck to each other at their basement, and their extension is in opposite directions (more details provided in Section 7.2). This topology increases the total length of the actuation unit. To consider this use case, another surrogate model needs to be developed.

The FEM model for the single actuator is reused. The boundary conditions are changed. A symmetry constraint is applied to the geometry at the actuator basement. The stress is picked up onto two critical points: at the interface between both actuators and at the interface between the output rod and the housing. Two surrogate models are expressed to determine the value of *π*<sup>0</sup> for each of these points.

## *6.7. DOE and Surrogate Synthesis*

A design of experiments (DoE) is realized with *e*1, *e*2, *drs*, and *La*. The simulations (a modal analysis) generate the variable of interest *π*0.

The dependent variable of the problem is approximated thanks to a linear regression (response surface model (RSM)) where the development takes into account a mean value, a first-order member (which represents the main effects of the problem), a combined member (representing the interactions), and a second-order member to consider further effects. The development takes the following form:

$$
\pi\_0 = \underbrace{a\_0}\_{\text{mean value}} + \underbrace{\sum a\_l \pi\_l}\_{\text{main effect}} + \underbrace{\sum a\_{lj} \pi\_l \pi\_j}\_{\text{in iterations}} + \underbrace{\sum a\_{li} \pi\_l^2}\_{\text{high order effect}} \tag{30}
$$

The VPLM methodology (Variable Power-Law Metamodel) [31] is applied. A log transformation on variables is performed for the linearization which gives the form:

$$\log(\pi\_0) = a\_0 + \sum a\_i \log(\pi\_i) + \sum a\_{i\bar{j}} \log(\pi\_i) \log(\pi\_{\bar{j}}) + \sum a\_{i\bar{i}} \log(\pi\_i)^2 \tag{31}$$

and can be rewritten as:

$$\pi\_0 = 10^{a\_0} \prod\_{i=1}^{n} \pi\_i \tag{32}$$

This variable power-law form enables to deal with the large variation range of the dependent and independent variables.

The data set coming from the DoE is shared in two sets: one for the regression procedure so as to determine the coefficients *a<sup>i</sup>* and *aij* (Equation (31)) and the other for the test of the final surrogate.

The regression gives the following surrogate model which determines the value of *π*<sup>0</sup> for the housing of a single actuator:

$$\begin{split} \log\_{10}(\pi\_{0}) = & 68 \cdot \log\_{10}\left(\frac{L}{d}\right)^{2} \cdot \frac{1}{293} + 64 \cdot \log\_{10}\left(\frac{L}{d}\right) \cdot \log\_{10}\left(\frac{e\_{1}}{d}\right) \cdot \frac{1}{973} + \log\_{10}\left(\frac{L}{d}\right) \cdot \log\_{10}\left(\frac{e\_{2}}{d}\right) \cdot \frac{1}{1000} \\ &+ 1099 \cdot \log\_{10}\left(\frac{L}{d}\right) \cdot \frac{1}{984} + 86 \cdot \log\_{10}\left(\frac{e\_{1}}{d}\right)^{2} \cdot \frac{1}{657} + 229 \cdot \log\_{10}\left(\frac{e\_{1}}{d}\right) \cdot \log\_{10}\left(\frac{e\_{2}}{d}\right) \cdot \frac{1}{884} \\ &- 101 \cdot \log\_{10}\left(\frac{e\_{1}}{d}\right) \cdot \frac{1}{249} + 538 \cdot \log\_{10}\left(\frac{e\_{2}}{d}\right)^{2} \cdot \frac{1}{685} + 670 \cdot \log\_{10}\left(\frac{e\_{2}}{d}\right) \cdot \frac{1}{359} + \frac{527}{551} \end{split} \tag{33}$$

with *L* the length *L<sup>a</sup>* and *d* the diameter *drs*.

The two other surrogate models developed for the tandem actuator housing are shown in Appendix A.

## *6.8. Validation*

In Figure 19, the prediction of the surrogate model (Equation (33)) is compared to the FEM simulation result data set. The prediction level is satisfying with *R* <sup>2</sup> > 99%. The two other housing surrogate models show the same prediction level quality.

**Figure 19.** Surrogate model for single actuator housing: *π*<sup>0</sup> results (Equation (33)).

## **7. Results**

In this section, some real application cases are discussed to illustrate the previously presented methodologies.

The methodology mentioned in Section 2 is implemented as a graphical user interface (GUI) web application based on a Jupyter Notebook calling functions through Python scripts. The GUI is developed with a dashboard named Voila. The tools developed in Section 4 are included with an optimization algorithm (differential evolution). Different specifications from real application cases are executed into this sizing code and the output results are presented. These specifications are arbitrary or linked to redundant topologies of actuation and to two different helicopter use cases: the main rotor and the tail rotor.

The considered actuator architecture is provided in Figure 20. In a housing (H), a frameless PMSM electrical motor (EM) is guided by two bearings (BB1 and BB2). BB1 is an angular contact double-row ball bearing; it withstands the entire axial load. BB2 is a radial contact single-row ball bearing; it ensures the motor rotor alignment. The electrical motor (EM) is linked to a screw mechanism (SM). In this paper, a planetary roller-screw (PRS) technology is considered. The screw spins and the nut moves linearly. The output rod (OR) is fixed onto the nut of the SM. A rod end (RE) (also named a spherical bearing) ensures the connexion of the output rod with the loading device. In this paper, the loading device is the helicopter swashplate. On the other side, another rod end (RE) ensures the connexion of the actuator housing (H) with a frame. Finally, the electromagnetic brake (EMB) satisfies a safety function in the case of a electricity supply cut-off.

For confidentiality matters, the inputs and results are presented as ratios of quantities. The observation of the sizing evolutions is the focus point of this section.

**Figure 20.** Actuator architecture considered in results.

## *7.1. Arbitrary Sizings*

We suggest presenting some preliminary sizing results from an arbitrary specification. The arbitrary specification and the associated design hypothesis are given in Table 5. The actuator sizing results regarding the simplex topology (Figure 20) are available in Table 6, in the column named *sizing A*. To show the sensitivity of the values of the static load, stroke, and acceleration (RMS), the arbitrary specification is run again three times with, every time, a modification of one indicator value. *Sizing B* includes the arbitrary static load multiplied by 2. *Sizing C* includes the arbitrary stroke multiplied by 4. *Sizing D* includes the arbitrary RMS acceleration multiplied by 2.

The following lines comment the sizing results presented in Table 6. Although the mechanical components are not the focus point of this article, they must be mentioned to better understand the sizing choices performed by the optimization.

The arbitrary specification leads a design of the mechanical components mainly guided by the fatigue criteria (*sizing A*). This means that the static load criteria are satisfied with margins. Moreover, because of the low heat-transfer coefficient, the low emissivity, and the significant RMS acceleration level (see Table 5), the electrical motor design is driven by the heat dissipated through the actuator skin, by convection and radiation.


**Table 5.** Arbitrary specification and its design hypothesis.


### **Table 6.** Sizing results from the arbitrary specification.

In *sizing B*, the increase in the static load stands in the static load margin of the mechanical components; their sizing remains unchanged. However, the electromagnetic brake (EMB) must stop a higher torque. Indeed, the torque to stop in an emergency relies on two specified values: the static load converted by the indirect efficiency of the screw (SM) and the screw maximum speed to stop within the specified time. Doubling the static load, the EMB requires to develop a higher braking force. It is bigger then. This involves the increase in the EMB disk inertia that sums to the total rotating inertia. Thus, a higher electrical motor performance is required. Because of the increased length of the electrical components, the housing is heavier.

Increasing the specified stroke (*sizing C*) makes the actuator longer. The ratio diameter by length is decreased. The housing and the output rod must have their thickness increased to withstand the vibratory accelerations. The screw of the SM is longer, involving an additional mass. What is more, the increased actuator length offers a more extended outer surface for dissipating the heat generated by the motor. The electrical motor (EM) can have lower performances if the thread lead is increased. The EM has a reduced mass then.

The motor heat generation is based on a continuous torque. The RMS acceleration highly contributes to this torque. Doubling the specified RMS acceleration (*sizing D*) involves an important specified continuous torque. The motor size must be increased to satisfy this specification. Increasing the lead limits the motor size increase for a reduced mass. With a bigger motor, the rotating inertia is higher and the EMB is required to be bigger. Moreover, increasing the lead reduces the fatigue phenomenon applied to the mechanical components. The mechanical components are chosen with slightly smaller fatigue capabilities. They are slightly smaller then. Moreover, the sized electrical components make the actuator longer. Thus, the housing (H) is slightly heavier.

## *7.2. Sizing of Redundant Topologies of Actuation*

There are mainly two redundant topologies of actuation [6,42,43]. The first one is the force summing where the failed or passive actuator shall be free in motion. In this case, the actuator must be equipped with a clutch or any breaking fuse system. The second one is the position summing where the failed or passive actuator shall be locked from motion. In this case, the actuators must be equipped with a power-off brake. The position summing topology can be seen with actuators either installed in tandem (as presented in Figure 21a) or in parallel (Figure 21b). The tandem configuration is the one most commonly found in aeronautics. This equips the PFCS of fighter aircrafts and helicopters and the SFCS of commercial aircrafts (except the spoiler) [6]. Meanwhile, the parallel configuration is much less used [6].

**Figure 21.** Redundant topologies of actuator in position summing (passive/failed actuator shall be locked). (**a**) One tandem actuator, (**b**) two simplex actuators linked by a cross bar (*a* = *b* considered, *F*<sup>1</sup> = *F*<sup>2</sup> = *F*3/2).

This paper proposes to study the impact on the housing mass and output rod mass involved in the consideration of topology (b) compared to topology (a) (Figure 21). It is clear that topology (b) introduces potential additional plays and wear points compared to topology (a). However, it is important to estimate if this topology involves any mass gain that would compensate these drawbacks.

From tandem to parallel topologies, the force is halved; meanwhile, the stroke, speed, and acceleration are doubled. The work produced by the actuator remains the same.

The bar chart presented in Figure 22 presents the sizing results involved in both actuation topologies. The results are given as a ratio, with the sizing results obtained for a simplex actuator not equipped with any electromagnetic brake.

**Figure 22.** Mass evolution of tandem or parallel topologies compared to simplex topology.

As the parallel topology requires a doubled stroke and half load, the actuator is longer and smaller in diameter than in the tandem topology. The thicknesses of the housing and the output rod are nearly doubled. The mass of the set (housing and rod) results to be half heavier than the one of the tandem topology. The parallel topology shows electrical components and mechanical ones with reduced characteristics and reduced masses compared to the tandem topology. Therefore, the contribution of the set (housing and rod) mass on the total actuator mass is much higher for the parallel topology. Figure 23 confirms it with a contribution of a third of the actuator mass concerning the tandem topology against more than half the actuator mass concerning the parallel topology.

**Figure 23.** Component mass distribution for one actuator.

The parallel configuration seemed to be a relevant choice because of the load reduction. However, this finally penalizes the actuator because the load reduction involves a small actuator diameter. With the increase in the stroke, the diameter-to-length ratio is not interesting anymore. There is not any mass gain on the housing and output rod. The cross bar has to be considered in the mass statement too.

Finally, the total actuation mass gain involved by the parallel topology is not significant enough relatively to the potential drawbacks it introduces.

## *7.3. Sizing of Main and Tail Rotor Actuators*

Section 3 showed that the TRA application was distinguished from the MRA in terms of dynamism, load, and especially in terms of induced motor losses. To illustrate this difference between the MRA and TRA applications (Table 1), their specifications are executed into the sizing code. The sized motor characteristics and an actuator mass distribution are displayed in Figures 24 and 25.

Figure 25 shows the important mass decrease among components induced by the TRA application. The TRA total mass shows to be 80% smaller than the actuator, satisfying both applications.

For confidentiality reasons, the results of the MRA and TRA sizing are normalized. The reference is taken as the actuator sizing which satisfies both applications at the same time.

Figure 24 globally shows that the selected motor for the TRA has reduced characteristics compared to the one selected for the MRA. The rotor inertia, which induces motor losses, is drastically lowered compared to the MRA-selected motor.

**Figure 24.** Motor characteristic evolution of specific actuators for the MRA and TRA applications compared to a unique actuator for both applications (reference).

**Figure 25.** Mass distribution of specific actuators for the MRA and TRA applications compared to a unique actuator for both applications.

Section 3 suggests the mean power rate *PR* and the RMS acceleration *arms* as a way to consider the induced motor losses. In Figure 24, the orange line shows the sizing resulting from the reference actuator specification regardless of *PR* and *arms*. The losses result to be estimated at roughly 40% of what they are for the reference actuator, and the continuous torque is at roughly 70%. The performance difference is also seen through the actuator mass statement in Figure 25. It is clear that not considering the values *PR* and *arms* in the actuator specification involves a significant risk of undersizing the electrical motor of the actuator.

The MRA and TRA of the considered helicopter PFCS lead to two significantly different sizes. The MRA and TRA are two applications to be distinguished. Designing a specific actuator for each application will benefit the helicopter mass and its electrical network.

## **8. Discussion**

This paper presents a design methodology supported by tools for the preliminary sizing of critical actuators. This methodology can be applied to any actuator architectures. This methodology finds an extension to the design of multirotor drones as presented in [44]. Moreover, it completes, at the component level, the methodology proposed by [45] at the vehicle level. The author of [45] developed a power system architecture sizing process for the preliminary design phase of civil aircraft.

The first tool presented is the one from [3]. It is a data-driven specification methodology which draws a parallel between measurement data on flight and EMA technologies using indicators to estimate over mission profiles.

Nowadays, the industrial context comes to a digital twin. Data content is globally exponentially increasing. As a result, there is a necessity to develop such data analysis methodologies based on synthetic values. What is more, the industrial trend is to reduce the number of helicopter flight tests to shrink development costs. Consequently, the methodology must extract as much added value as possible from any available data. The paper [3] presents a statistical approach getting interest in this area. From a reduced data set, the statistical laws are formulated. The laws are used then to express the load, speed, and acceleration domain limits of the application.

Furthermore, as seen in this paper, the specification methodology applied on main and tail rotor applications taught about the helicopter specificities. It showed a clear dynamism and load difference between the MRA and the TRA. The indicators quantified the importance of the motor losses on the TRA compared to the MRA. These losses are induced by the mean power rate *PR* and especially the equivalent continuous acceleration *arms*. The last part of the paper showed the risk of undersizing if *arms* and *PR* are not taken into account. In addition, it concludes that the helicopter mass could benefit from a specific sizing for the MRA and the TRA. The TRA shall be designed with a rotor inertia as low as possible.

The second tool presented is the scaling laws based on a dimensional analysis. The electrical motor is chosen as an applicative example. The scaling laws capture the main physical phenomena driving the component design and easily estimate the main component

characteristics from a reduced number of parameters. As illustrated by the last part of this paper, the scaling laws are useful for the exploration study of a design domain: sizing with optimization, a scenario analysis, an integration study (mass and dimensions), and a technology or architecture comparison. In addition, scaling laws can be useful to replace component catalogs when they are unavailable or to complete missing contents. In the context of negotiation, scaling laws are easy enough to exchange with suppliers and to challenge them.

The third tool presented is a surrogate model or response surface model (RSM). It is set up using dimensionless numbers (or *π*-numbers). This has several advantages. First, it decreases the number of variables to be manipulated and therefore it drastically decreases the number of physical or numerical experiments to be carried out. Secondly, it increases the regression robustness [32], in particular if the RSM is built within the logarithmic space. Paper [17] mentions that the logarithmic shows good results in interpolation and in extrapolation because of the power-law form. This is the case for the VPLM methodology (Variable Power-Law Metamodel) [31] used in the surrogate model setup in this paper.

The complex vibratory problem of the actuator housing is suggested to be addressed by a surrogate model because it is adapted to preliminary studies. Indeed, it is a simple way of modeling with a reduced number of parameters. The suggested model easily provides design trends for architecture decision making. This is shown in the last part of this paper where the actuator topologies are compared. Obviously, the final actuator housing design will require laboratory experiments on a vibratory test bench.

Finally, as perspectives, the presented design methodology can be applied to any other actuator architectures. Thus, for a given application, it allows to study and select the best actuator architecture. For critical actuators, an analysis among different redundant topologies of actuation and a safety analysis with failure trees are important features to develop. In addition, a price assessment is a relevant feature to implement.

**Author Contributions:** Conceptualization, J.R., M.B. and L.R.; methodology, J.R., M.B. and L.R.; software, J.R. and M.B.; validation, J.R., M.B. and L.R.; formal analysis, J.R. and M.B.; investigation, J.R. and M.B.; resources, J.R. and M.B.; data curation, J.R. and M.B.; writing—original draft preparation, J.R.; writing—review and editing, J.R. and M.B.; visualization, J.R.; supervision, M.B. and L.R.; project administration, J.R.; funding acquisition, J.R. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Airbus Helicopters and the French Minister of Higher Education, Research and Innovation (MESRI) through the French National Association of Research and Technology (ANRT).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The authors would like to acknowledge Airbus and the ANRT who funded the PhD thesis that produced the work presented in this paper.

**Conflicts of Interest:** The authors declare no conflict of interest.

## **Nomenclature**



## **Abbreviations**

The following abbreviations are used in this manuscript:


## **Appendix A. Surrogate Model: Housing of the Tandem Topology of Actuation**

For the tandem topology of actuation and at the interface between output rod and housing, the regression gives the following surrogate model:

$$\begin{split} \log\_{10}(\pi\_{0}) &= 180 \cdot \log\_{10}\left(\frac{L}{d}\right)^{2} \cdot \frac{1}{983} + 23 \cdot \log\_{10}\left(\frac{L}{d}\right) \cdot \log\_{10}\left(\frac{e\_{1}}{d}\right) \cdot \frac{1}{243} + 43 \cdot \log\_{10}\left(\frac{L}{d}\right) \cdot \log\_{10}\left(\frac{e\_{2}}{d}\right) \cdot \frac{1}{880} \\ &+ 1187 \cdot \log\_{10}\left(\frac{L}{d}\right) \cdot \frac{1}{766} + 149 \cdot \log\_{10}\left(\frac{e\_{1}}{d}\right)^{2} \cdot \frac{1}{993} + 179 \cdot \log\_{10}\left(\frac{e\_{1}}{d}\right) \cdot \log\_{10}\left(\frac{e\_{2}}{d}\right) \cdot \frac{1}{712} \\ &- 175 \cdot \log\_{10}\left(\frac{e\_{1}}{d}\right) \cdot \frac{1}{429} + 625 \cdot \log\_{10}\left(\frac{e\_{2}}{d}\right)^{2} \cdot \frac{1}{791} + 1229 \cdot \log\_{10}\left(\frac{e\_{2}}{d}\right) \cdot \frac{1}{694} + \frac{922}{999} \end{split} \tag{A1}$$

with *L* the length *L<sup>a</sup>* and *d* the diameter *drs*.

For the tandem topology of actuation and at the interface between both actuators, the regression gives the following surrogate model:

$$\begin{split} \log\_{10}(\pi\_{0}) &= 96 \cdot \log\_{10}\left(\frac{L}{d}\right)^{2} \cdot \frac{1}{899} + 25 \cdot \log\_{10}\left(\frac{L}{d}\right) \cdot \log\_{10}\left(\frac{e\_{1}}{d}\right) \cdot \frac{1}{974} + 8 \cdot \log\_{10}\left(\frac{L}{d}\right) \cdot \log\_{10}\left(\frac{e\_{2}}{d}\right) \cdot \frac{1}{223} \\ &+ 1012 \cdot \log\_{10}\left(\frac{L}{d}\right) \cdot \frac{1}{585} - 8 \cdot \log\_{10}\left(\frac{e\_{1}}{d}\right)^{2} \cdot \frac{1}{819} - 43 \cdot \log\_{10}\left(\frac{e\_{1}}{d}\right) \cdot \log\_{10}\left(\frac{e\_{2}}{d}\right) \cdot \frac{1}{974} \\ &- 84 \cdot \log\_{10}\left(\frac{e\_{1}}{d}\right) \cdot \frac{1}{811} + 2 \cdot \log\_{10}\left(\frac{e\_{2}}{d}\right)^{2} \cdot \frac{1}{15} - 259 \cdot \log\_{10}\left(\frac{e\_{2}}{d}\right) \cdot \frac{1}{372} - \frac{175}{731} \end{split} \tag{A2}$$

with *L* the length *L<sup>a</sup>* and *d* the diameter *drs*.

## **References**


## *Article* **Incremental Nonlinear Dynamic Inversion Attitude Control for Helicopter with Actuator Delay and Saturation**

**Shaojie Zhang \*, Han Zhang and Kun Ji**

College of Automation Engineering, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China **\*** Correspondence: zhangsj@nuaa.edu.cn

**Abstract:** In this paper, an incremental nonlinear dynamic inversion (INDI) control scheme is proposed for the attitude tracking of a helicopter with model uncertainties, and actuator delay and saturation constraints. A finite integral compensation based on model reduction is used to compensate the actuator delay, and the proposed scheme can guarantee the semi-globally uniformly ultimately bounded tracking. The overall attitude controller is separated into a rate, an attitude, and a collective pitch controller. The rate and collective pitch controllers combine the proposed method and INDI to enhance the robustness to actuator delay and model uncertainties. Considering the dynamic of physical actuators, pseudo-control hedging (PCH) is introduced both in the rate and attitude controller to improve tracking performance. By using the proposed controller, the helicopter shows good dynamics under the multiple restrictions of the actuators.

**Keywords:** incremental nonlinear dynamic inversion (INDI); actuator compensation; model reduction; pseudo-control hedging (PCH); helicopter attitude control

## **1. Introduction**

The helicopter, aircraft with one or more power-driven horizontal propellers or rotors that enable it to take off and land vertically, to move in any direction, or to remain stationary in the air, has become very popular for a wide range of services, including air–sea rescue, firefighting, traffic control, oil platform resupply, and business transportation [1]. However, these tasks often bring heavy workloads to pilots, especially in situations of high crosswinds or low light. Furthermore, subject to a complicated dynamic response, multiple flight modes, system uncertainties, and rapidly varying flight conditions, the helicopter is a highly complex system. For the reasons above, a highly reliable and effective flight control system which allows the helicopter to execute multiple tasks in adverse flying conditions becomes more in demand [2].

In the past few years, various fight control methodologies are developed for the helicopter [3–6]. Feedback linearization, as the most widely used nonlinear control method in the aircraft control systems, is often combined with adaptive control to deal with the model uncertainties [7–10]. However, it is hard to guarantee that the control system can recover from a failure in adaptation [11]. Therefore, whether it can be applied in the systems with high security requirements is worthy of consideration. In order to overcome the shortcomings above, the INDI technique is adopted in this paper for helicopter flight control.

By producing the incremental form of the control command by calculating the error between the virtual control law and the acceleration of the system state, INDI is robust to model uncertainties [12–19]. In [17], the stability and robustness of the INDI technique has been proven. The INDI scheme was first used in the design of a six-degree-of-freedom helicopter controller in [18], and its robustness to model uncertainties was verified by simulation. Ref. [19] uses the INDI technique to redesign the existing Apache flight control and improve the handling quality.

However, the weakness of the INDI controller is the accuracy of the onboard measurement and actuator delay. The current measurement technology combined with data

**Citation:** Zhang, S.; Zhang, H.; Ji, K. Incremental Nonlinear Dynamic Inversion Attitude Control for Helicopter with Actuator Delay and Saturation. *Aerospace* **2023**, *10*, 521. https://doi.org/10.3390/ aerospace10060521

Academic Editor: Gianpietro Di Rito

Received: 24 October 2022 Revised: 16 May 2023 Accepted: 24 May 2023 Published: 1 June 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

processing algorithms (such as Kalman filtering) has been able to reduce the uncertainties brought by these sensors. However, there is no effective solution to the poor performance caused by the existence of actuator delay. According to [1], the delay time is approximately 100 ms in the helicopter. In a control system operating at 100 Hz, there will be a difference of about ten samples between the command signal and the actuator, which seriously deteriorates the tracking performance of the system and even puts the system stability in risk. Therefore, various control approaches are proposed for the systems subject to actuator delay [20–27]. Ref. [25] extends the Artstein model reduction method in [26] to nonlinear systems, and designs a compensation control law for known and unknown systems, respectively. Then, the Lyapunov–Krasovskii functional is adopted to prove the stability of them. Ref. [27] designs a feedback robust tracking controller with delay compensation for a class of systems with actuator delay and external disturbance, and achieved the desired effect. In [28,29], Rohollah Moghadam proposed an ADP-based solution to the optimal adaptive adjustment problem of systems with state delay and input delay, which can be applied to the optimal control problem of a class of nonlinear time-delay systems.

Besides the constraint of time delay, actuator saturation is very common and a more general problem in the helicopter due to the limitation of actuators in terms of the position and rate saturation. Therefore, many recent works have been carried out on actuator saturation nonlinearity as it causes the windup phenomenon [30–35]. Based on nonlinear partial differential inequalities, an optimal saturation compensator was developed in [32]. In [11], the pseudo-control hedging method is proposed to offset the virtual control input when the actuators are saturated. Ref. [33] further expanded the PCH theory and [34] applied the method to the development of the Boeing 747 flight control system.

In this paper, combining with model reduction and the PCH technique, a novel INDI-based controller is proposed for helicopter attitude control with actuator delay and saturation. The main contributions of this paper are the following: (1) a novel INDIbased actuator delay compensation control scheme with guaranteed stability is proposed, which can be applied to nonlinear systems with actuator delay in a certain range; and (2) aiming at helicopter characteristics, the proposed method and PCH are adopted to design the controller for the helicopter which is subjected to model uncertainty, actuator delay, and saturation.

The overall structure of this paper is as follows: The problem is formulated in Section 2. Section 3 presents an actuator delay compensation scheme for the INDI controller. The introducing of the helicopter model and the design of the anti-windup helicopter attitude controller by the proposed method is given in Section 4. Section 5 focuses on the display of the simulation results and related explanations. The conclusions are presented in Section 6.

## **2. Problem Formulation**

Consider a class of the nonlinear system with input delay in the following form:

$$\begin{cases} \dot{\mathbf{x}} = f(\mathbf{x}) + \mathbf{g}(\mathbf{x}, u\_{\tau}) \\ \mathbf{y} = \mathbf{x} \end{cases} \tag{1}$$

where *x* is the state vector in R*<sup>n</sup>* . *f*(·) and *g*(·) are nonsingular, bounded continuous smooth nonlinear functions in R*<sup>n</sup>* . The delayed input vector *u<sup>τ</sup>* = *u*(*t* − *τ*) and output vector *y* are both functions in R*<sup>n</sup>* , where *<sup>τ</sup>* <sup>∈</sup> <sup>R</sup><sup>+</sup> is the time-delay constant. The following assumptions are used for further development.

**Assumption 1.** *The desired trajectory vector x<sup>d</sup>* (*t*) *is bounded and has bounded time derivatives up to i-th for i* = 1, 2, 3*; that is, there exists η<sup>i</sup> such that x* (*i*) *d* (*t*) <sup>≤</sup> *<sup>η</sup><sup>i</sup> for all t, where* k·k *represents the standard Euclidean norm.*

**Assumption 2.** *Both state x and the delayed input u<sup>τ</sup> are measurable. Furthermore, the first derivative of <sup>x</sup> in the last sampling time, denoted by* . *x*0*, can also be measured by acceleration sensors.* **Assumption 3.** *If x*(*t*), . *x*(*t*)*, and u*(*t*) *are bounded, then f*(*x*) *and g*(*x*, *u*) *are bounded. Moreover, the first partial derivatives of f*(*x*) *and g*(*x*, *u*) *with respect to their arguments x and u exist and are bounded. The infinity norm of the partial derivative of g*(*x*, *u*) *to u around* (*x*0, *u*0) *and its inverse can be upper-bounded as*

$$\left\| \left. \frac{\partial \mathbf{g}(\mathbf{x}, \boldsymbol{\mu})}{\partial \boldsymbol{u}} \right|\_{\mathbf{x}\_{0}, \boldsymbol{\mu}\_{0}} \right\|\_{\infty} \leq \mathfrak{f}\_{1}, \left\| \left. \left( \frac{\partial \mathbf{g}(\mathbf{x}, \boldsymbol{\mu})}{\partial \boldsymbol{u}} \right) \right\|\_{\mathbf{x}\_{0}, \boldsymbol{\mu}\_{0}}^{-1} \right\|\_{\infty} \leq \mathfrak{f}\_{2} \tag{2}$$

*where <sup>ξ</sup>*1, *<sup>ξ</sup>*<sup>2</sup> <sup>∈</sup> <sup>R</sup><sup>+</sup> *are known constants.*

## **3. Control Law Design**

Before the development of the control law, we define some variables for subsequent analysis. We denote the tracking error as

$$e = \mathfrak{x}\_d(t) - \mathfrak{x}(t) \tag{3}$$

Then, we give the virtual control from the tracking error, denoted by

$$\boldsymbol{\nu} = \dot{\boldsymbol{x}}\_d + k\_\nu \boldsymbol{e} \tag{4}$$

where *<sup>k</sup><sup>ν</sup>* <sup>∈</sup> <sup>R</sup>*n*×*<sup>n</sup>* denotes a positive gain matrix.

To facilitate the subsequent analysis, an auxiliary tracking error which is inspired from the model reduction method is defined as

$$
\dot{\mathbf{r}} = \dot{e} + k\_r e + \mathbf{g}(\mathbf{x}, \boldsymbol{\mu}(t-\tau)) - \mathbf{g}(\mathbf{x}, \boldsymbol{\mu}(t)) \tag{5}
$$

where *<sup>k</sup><sup>r</sup>* <sup>∈</sup> <sup>R</sup>*n*×*<sup>n</sup>* is a known positive gain matrix.

Submitting the expressions in (1) and (3), the transformed open-loop tracking error system can be represented in an input delay free form as

$$
\dot{\sigma} = \dot{\mathbf{x}}\_d + k\_r e - f(\mathbf{x}) - \mathbf{g}(\mathbf{x}, \boldsymbol{\mu}) \tag{6}
$$

We rewrite (6) as its partial first-order Taylor series expansion around the current solution of the system, denoted by (*x*0, *u*0):

$$\begin{aligned} \dot{\tau} &= \dot{\mathbf{x}}\_d + k\_r e - f(\mathbf{x}\_0) - \mathbf{g}(\mathbf{x}\_0, \boldsymbol{\mu}\_0) - \left. \frac{\partial}{\partial \mathbf{x}} [f(\mathbf{x}) + \mathbf{g}(\mathbf{x}, \boldsymbol{\mu})] \right|\_{\mathbf{x}\_0, \boldsymbol{\mu}\_0} (\mathbf{x} - \mathbf{x}\_0) \\ &- \frac{\partial}{\partial \mathbf{u}} [f(\mathbf{x}) + \mathbf{g}(\mathbf{x}, \boldsymbol{\mu})] \Big|\_{\mathbf{x}\_0, \boldsymbol{\mu}\_0} (\boldsymbol{\mu} - \mathbf{u}\_0) \end{aligned} \tag{7}$$

where *x*<sup>0</sup> and *u*<sup>0</sup> denote the values in the last control step of *x* and *u*, respectively. Based on the assumption on INDI, the variation *x* − *x*<sup>0</sup> can be neglected for each small time increment. Then, (7) can be simplified by

$$
\dot{\boldsymbol{\sigma}} = \dot{\boldsymbol{\pi}}\_d + k\_r \boldsymbol{e} - f(\mathbf{x}\_0) - \boldsymbol{g}(\mathbf{x}\_0, \boldsymbol{\mu}\_0) - G(\mathbf{x}\_0, \boldsymbol{\mu}\_0)(\boldsymbol{u} - \boldsymbol{u}\_0) \tag{8}
$$

where *<sup>G</sup>*(*x*0, *<sup>u</sup>*0) = *<sup>∂</sup>g*(*x*,*u*) *∂u x*0,*u*<sup>0</sup> . Based on (8) and the INDI control law, the control law *u*(*t*) can be given by

$$u = u\_0 + G(\mathbf{x}\_0, u\_0)^{-1} \left[ (\nu - \dot{\mathbf{x}}\_0) + k\_\mu \int\_{t\_0}^t \left[ \dot{e} + k\_\tau e + \mathbf{g}(\mathbf{x}, u(\theta - \tau)) - \mathbf{g}(\mathbf{x}, u(\theta)) \right] d\theta \right] \tag{9}$$

where *<sup>k</sup><sup>u</sup>* <sup>∈</sup> <sup>R</sup>*n*×*<sup>n</sup>* is a positive control gain matrix. The control law in (9) can be thought of as a combination of an INDI-based term through an online state measurement and a predictor term which can stabilize the system and compensate the input delay. Note that the control law in (9) does not directly depend on *f*(*x*) anymore, which means the controller is robust to the uncertainty of the model that only depends on the states of the system. However, the uncertainty in the control effectiveness matrix ∆*G* should meet k∆*G*k ≤ 0.5k*G*k [17]. Compared with the model-based feedforward control method in [27], the proposed controller achieves a good control effect even when the model is uncertain. After applying the Expression (4) and (9), (8) can be expressed as

$$\sigma = (\mathbf{k}\_{\varGamma} - \mathbf{k}\_{\varGamma})\mathbf{e} + \mathbf{g}(\mathbf{x}\_{0\varGamma}\mathbf{u}\_{\varPi}) - \mathbf{g}(\mathbf{x}\_{0\varGamma}\mathbf{u}\_{0}) - \mathbf{k}\_{\varPi} \int\_{t\_{0}}^{t} [\dot{\mathbf{e}} + \mathbf{k}\_{\varGamma}\mathbf{e} + \mathbf{g}(\mathbf{x}, \mathbf{u}(\theta - \mathsf{T})) - \mathbf{g}(\mathbf{x}, \mathbf{u}(\theta))] d\theta \tag{10}$$

where *g*(*x*0, *uτ*<sup>0</sup> ) is the control effectiveness function under the last state vector *x*<sup>0</sup> and the last delayed input vector *uτ*<sup>0</sup> . Then, the time derivative of (10) can be obtained by

$$
\dot{r} = (k\_r - k\_\nu)\dot{e} - k\_\nu r \tag{11}
$$

In addition, we can also get the time derivative of (9) by using (8):

$$
\dot{\mu} = G(\mathbf{x}\_0 \,\mu\_0)^{-1} (\dot{\nu} + k\_\mu r) \tag{12}
$$

Recall the Assumption 3: one thing that can be determined is that the discrepancy between the *g*(*x*, *u*) and *g*(*x*, *uτ*) upper bound is

$$\log(\mathbf{x}, \boldsymbol{\mu}) - \mathbf{g}(\mathbf{x}, \boldsymbol{\mu}\_{\tau}) \le \rho(\|\boldsymbol{e}\_{\boldsymbol{\mu}}\|) \|\boldsymbol{e}\_{\boldsymbol{\mu}}\|\tag{13}$$

where *<sup>e</sup><sup>u</sup>* <sup>∈</sup> <sup>R</sup>*<sup>n</sup>* is defined as

$$e\_{\mathfrak{u}} = \mathfrak{u} - \mathfrak{u}\_{\mathsf{T}} = \int\_{t-\mathsf{T}}^{t} \dot{\mathfrak{u}}(\theta)d\theta \tag{14}$$

In addition, the bounding function *ρ*(*eu*) ∈ R is a known positive globally invertible nondecreasing function.

**Theorem 1.** *The control law given in (9) can ensure the semi-globally uniformly ultimately bounded (SUUB) tracking in the case that*

$$0 \le \tau \le \frac{\chi}{2\xi\varrho^2} \tag{15}$$

*where χ is subsequently defined control gains, provided the control gains kν*, *k<sup>r</sup>* , *ku*, *ω are selected with the following sufficient conditions:*

$$\begin{cases} \begin{array}{c} k\_r > \frac{\chi^2 + 2}{4} \\ k\_V = k\_r \end{array} \\ \begin{array}{c} \frac{1 - \delta}{2\omega \tau \xi\_2 2} < k\_{\mu} < \frac{1 + \delta}{2\omega \tau \xi\_2 2} \\ \tau \le \omega \le \frac{\chi^2}{4\tau \xi\_2 2} \end{array} \end{cases} \tag{16}$$

*where δ* = q <sup>1</sup> <sup>−</sup> <sup>4</sup> *<sup>χ</sup>*<sup>2</sup> *ωτξ*<sup>2</sup> <sup>2</sup> *and ω are subsequently defined constants.*

**Proof.** Define *<sup>y</sup><sup>a</sup>* <sup>∈</sup> <sup>R</sup>2*n*+<sup>1</sup> as

$$y\_a = \begin{bmatrix} e^T & r^T & \sqrt{Q} \end{bmatrix}^T \tag{17}$$

where the positive definite LK functional *Q* is defined by

$$Q = \omega \int\_{t-\tau}^{t} \left( \int\_{s}^{t} ||\dot{\mu}(\gamma)||^{2} d\gamma \right) ds \tag{18}$$

where *ω* ∈ R is a positive constant. Then, the positive definite Lyapunov functional *V* is defined as

$$V(y\_a) = \frac{1}{2}e^T e + \frac{1}{2}r^T r + Q \tag{19}$$

satisfying the following inequalities

$$c\_1 \|y\_a\|^2 \le V \le c\_2 \|y\_a\|^2 \tag{20}$$

where *<sup>c</sup>*1, *<sup>c</sup>*<sup>2</sup> <sup>∈</sup> <sup>R</sup><sup>+</sup> are known constants.

After submitting the Equations (5) and (11), we have the time derivative of (19) as

$$\begin{array}{l}\dot{V} = e^{T}r - k\_{r}e^{T}e + e^{T}(\operatorname{g}(\mathbf{x},\boldsymbol{u}) - \operatorname{g}(\mathbf{x},\boldsymbol{u}\_{\tau})) \\ + (k\_{r}k\_{\nu} - k\_{r}^{\prime})r^{T}e + (k\_{r} - k\_{\nu})r^{T}(\operatorname{g}(\mathbf{x},\boldsymbol{u}) - \operatorname{g}(\mathbf{x},\boldsymbol{u}\_{\tau})) \\ - (k\_{\nu} - k\_{r})r^{T}r - k\_{\nu}r^{T}r + \omega\tau||\dot{\boldsymbol{u}}||^{2} - \omega\int\_{t-\tau}^{t}||\dot{\boldsymbol{u}}(\boldsymbol{\gamma})||d\boldsymbol{\gamma} \end{array} \tag{21}$$

Combining (13) and canceling common terms, we can upper-bound (21) in the case of *k<sup>r</sup>* = *k<sup>v</sup>* as

$$\dot{V} \le \left\| e \right\| \left\| r \right\| - k\_{\tau} \left\| e \right\|^{2} - k\_{\nu} \left\| r \right\|^{2} + \omega \tau \left\| \dot{u} \right\|^{2} + \rho(\left\| e\_{\bar{u}} \right\|) \left\| e \right\| \left\| e\_{\bar{u}} \right\| - \omega \int\_{t-\tau}^{t} \left\| \dot{u}(\gamma) \right\|^{2} d\gamma \tag{22}$$

According to Young's inequality, the following relation can be obtained

$$||e||\|r|| \le \frac{\chi^2}{4}||e||^2 + \frac{1}{\chi^2}||r||^2\tag{23}$$

$$\rho(||e\_{\boldsymbol{u}}||)||\boldsymbol{e}||\,||\boldsymbol{e}\_{\boldsymbol{u}}|| \leq \frac{1}{2}||\boldsymbol{e}||^{2} + \frac{\rho^{2}(||\boldsymbol{e}\_{\boldsymbol{u}}||)}{2}||\boldsymbol{e}\_{\boldsymbol{u}}||^{2} \tag{24}$$

where *χ* ∈ *R* <sup>+</sup> is a known constant. Moreover, using the Cauchy–Schwarz inequality, the terms in (24) can be further upper-bounded as

$$||e\_{\mu}|| \le \tau \int\_{t-\tau}^{t} ||\dot{u}(\gamma)||^{2} d\gamma \tag{25}$$

By adding and subtracting *τ* R *t t*−*τ* . *u*(*γ*) 2 *dγ* , inequality (22) can be upper-bounded as

$$\begin{array}{l} \dot{V} \le -(k\_{\tau} - \frac{\chi^2}{4} - \frac{1}{2})||e||^2 + \omega \tau ||\dot{u}||^2 - \left(k\_{\mu} - \frac{1}{\chi^2}\right)||r||^2\\ -\tau \int\_{t-\tau}^{t} ||\dot{u}(\gamma)||^2 d\gamma - \left(\frac{\omega}{\tau} - \frac{1}{2}\rho^2(||e\_{\boldsymbol{u}}||) - 1\right)||e\_{\boldsymbol{u}}||^2 \end{array} \tag{26}$$

Recall the Equation (12) and Assumption 3: there exists a positive constant *λ* ∈ R such that

$$\left\|\left\|\dot{u}\right\|\right\|^2 \le \lambda + \mathfrak{F}\_2^{\ 2} k\_u^{\ 2} \left\|r\right\|^2 \tag{27}$$

Then, Equation (26) can be rewritten as

$$\dot{V} \le -\kappa \|y\_a\|^2 - \left(\frac{\omega}{\tau} - \frac{1}{2}\rho^2(\|e\_{\mu}\|) - 1\right) \|e\_{\mu}\|^2 + \omega\tau\lambda \tag{28}$$

where the function *<sup>κ</sup>* <sup>∈</sup> <sup>R</sup><sup>+</sup> in (28) is defined as

$$\kappa = \min\left[ (k\_r - \frac{\chi^2}{4} - \frac{1}{2}), (k\_u - \frac{1}{\chi^2} - \omega \tau \xi\_2^\* 2^2 k\_u^2), (\frac{1}{\omega}) \right] \tag{29}$$

where the inequation

$$\begin{split} \boldsymbol{Q} &\leq \omega \tau \sup\_{\boldsymbol{s}\in[t, t-\tau]} \left[ \int\_{\boldsymbol{s}}^{t} \left|| \dot{\boldsymbol{u}}(\boldsymbol{\gamma}) \right||^{2} d\boldsymbol{\gamma} \right] \\ &= \omega \tau \int\_{t-\tau}^{t} \left|| \dot{\boldsymbol{u}}(\boldsymbol{\gamma}) \right||^{2} d\boldsymbol{\gamma} \end{split} \tag{30}$$

is used. Substituting the bound given in (20), the inequality in (28) can be further upperbounded as .

$$
\dot{V} \le -\frac{\kappa}{c\_2} V + \omega \tau \lambda \tag{31}
$$

in the set *S* defined by

$$\mathcal{S} = \left\{ e\_{\mathcal{U}}(t) \, \vert \, \vert \, \vert e\_{\mathcal{U}} \vert \, < \rho^{-1} \left( \sqrt{\frac{2\omega}{\tau} - 2} \right) \right\} \tag{32}$$

Hence, Expression (31) can lead to the solution as

$$V \le V(0)e^{-\frac{\mathcal{E}}{c\_2}t} + \omega \tau \lambda \frac{c\_2}{\mathcal{\varphi}} (1 - e^{-\frac{\mathcal{E}}{c\_2}t}) \tag{33}$$

In the case of *eu*(0)= 0 ∈ S and according to the definition of *V* in (19) and the solution in (19), it can be concluded that *e*,*r* are bounded. From Assumption 1 and the bounded *r*, we can infer that the variable *e<sup>u</sup>* = R *t t*−*τ* . *u*(*θ*)*dθ* = *G*(*x*0, *u*0) −1 R *t t*−*τ* ( . *x<sup>d</sup>* + *kνe* − .. *x*<sup>0</sup> + *kur*)*dθ* is bounded. Then, using the definition of *r* in (5) and combining the bounded tracking error *e* and Assumption 1, we can infer that both *<sup>x</sup>*(*t*) and . *x*(*t*) are also bounded. Finally, we can obtain from Equation (9) that *u* is bounded with the initial condition *u*(0) = 0 . Therefore, all of variables in the closed loop is guaranteed to be bounded by the proposed control law.

## **4. Attitude Controller Design for Helicopter**

## *4.1. Helicopter Model*

The attitude control of the Messerschmitt–Bölkow–Blohm (MBB) Bö-105 helicopter is considered in this paper. Here, we give the attitude model in the equations of motion from

$$
\dot{\omega}\_{\hbar} = \begin{bmatrix} \dot{p} \\ \dot{q} \\ \dot{r} \end{bmatrix} = J^{-1}(m - \omega\_{\hbar} \times J \omega\_{\hbar})\tag{34}
$$

$$\dot{\boldsymbol{\theta}} = \begin{bmatrix} \dot{f} \\ \dot{\theta} \\ \dot{\psi} \end{bmatrix} = \boldsymbol{\Omega}\boldsymbol{\omega}\_{\boldsymbol{h}} = \begin{bmatrix} 1 & \sin\phi\tan\theta & \cos\phi\tan\theta \\ 0 & \cos\phi & -\sin\phi \\ 0 & \sin\phi/\cos\theta & \cos\phi/\cos\theta \end{bmatrix} \boldsymbol{\omega}\_{\boldsymbol{h}} \tag{35}$$

where *m* = *mf us*(*x*) + *mht*(*x*) + *mvt*(*x*) + *mmr*(*x*, *u*) + *mtr*(*x*, *u*) represents the total moments with respect to the gravity center of the helicopter. It consists of the moments produced by the fuselage, horizontal tail, vertical tail, main rotor, and tail rotor, which are represented by the subscript *f us*, *ht*, *vt*, *mr*, and *tr*, respectively. *ω<sup>h</sup>* = - *p q r<sup>T</sup>* and *θ* = - *φ θ ψ<sup>T</sup>* indicate the roll, pitch and yaw angular rate and attitude angle of helicopter, respectively. *J* is the inertia matrix of the helicopter which is given by

$$J = \begin{bmatrix} 1433 & 0 & -660 \\ 0 & 4973 & 0 \\ -660 & 0 & 4099 \end{bmatrix} \tag{36}$$

The controller proposed in this paper is based on the principle of time-scale separation, which assumes that the state variables in the inner loop are preforming fast while the related parameters in the outer loop change more slowly under the same actuator inputs. Table 1 can verify that this hypothesis is reasonable.

**Table 1.** Time-scale separation between attitude angles and rates.


From Table 1 we can see that there exists a time-scale separation between the time derivative of attitude angles and rates. Therefore, we divide these six state variables into two loops for the controller design, namely, the rate loop and the attitude loop. This type of assumption is often carried out for flight dynamics and control applications. Between two loops, the parameters associated with the slow dynamics are treated as constants by the fast dynamics and its dynamic inversion is performed assuming that the states controlled in the inner loop achieve their commanded values instantaneously. The fast variables are thus used as control inputs to the slow dynamics. are not considered, the tracking performance of the control system will be severely degraded and even face stability problems. To overcome these problems, the proposed method and PCH are adopted in the next sections. **Table 2.** Actuators' rate limits of Bö-105.

However, what needs to be considered is the dynamic of the actuators. In fact, there exists a time delay between the controller delivering the signal to the actuators. Moreover, actuators are often limited by their moving rate, which is shown in Table 2. If these issues are not considered, the tracking performance of the control system will be severely degraded and even face stability problems. To overcome these problems, the proposed method and PCH are adopted in the next sections. **Actuator Name Variable Name Maximum Rate Limit**  Main rotor θ<sup>0</sup> 16.0deg/*s* Longitudinal cyclic θ<sup>1</sup>*<sup>s</sup>* 28.8deg/*s* Lateral cyclic θ<sup>1</sup>*<sup>c</sup>* 16.0deg/*s*

<sup>0</sup>*tr* 32.0deg/*s*

ν

in (4)

(37)

**Table 2.** Actuators' rate limits of Bö-105. Tail rotor


θ

### *4.2. Anti-Windup Design* error signal of the control system is too radical for actuators, it allows the production of a

To overcome the effect caused by the actuator saturation, the pseudo-control hedging (PCH) scheme is adopted. The PCH solves the effect of actuator saturation by modifying the virtual control input instead of directly influencing the actuator input. When the input error signal of the control system is too radical for actuators, it allows the production of a signal *ν<sup>s</sup>* which is opposite to the virtual control law to the first-order reference model, so that it can prevent the system from still trying to track the commanded references when actuator saturation occurs. signal ν *<sup>s</sup>* which is opposite to the virtual control law to the first-order reference model, so that it can prevent the system from still trying to track the commanded references when actuator saturation occurs. In order to achieve the control hedging, a first-order reference model is introduced as Figure 1. *com I* represents the command signal and *sat I* is the filtered signal to ensure ν

In order to achieve the control hedging, a first-order reference model is introduced as Figure 1. *Icom* represents the command signal and *Isat* is the filtered signal to ensure the input is within the acceptable range of the system. *ν<sup>s</sup>* is defined, corresponding to the input under the actuator dynamics. The maximum rate change allowed for this helicopter follows the ADS-33E-PRF standard, that is, the limits of 40 degrees per second on the roll and pitch rates and 80 degrees per second on the yaw rates. By using the reference model, we can obtain the time derivative of *Icom* such that the virtual control *ν* in (4) can be made easily when no saturation occurs. the input is within the acceptable range of the system. *<sup>s</sup>* is defined, corresponding to the input under the actuator dynamics. The maximum rate change allowed for this helicopter follows the ADS-33E-PRF standard, that is, the limits of 40 degrees per second on the roll and pitch rates and 80 degrees per second on the yaw rates. By using the reference model, we can obtain the time derivative of *com I* such that the virtual control can be made easily when no saturation occurs.

**Figure 1.** First-order reference model with saturation filter. **Figure 1.** First-order reference model with saturation filter.

The PCH signal *ν<sup>s</sup>* in the affine nonlinear systems can be calculated by

*s des*

$$\begin{array}{l} \nu\_{\rm s} = \left[ \dot{\mathbf{x}}\_{0} + \mathbf{G} (\mathbf{x}\_{0}, \boldsymbol{\mu}\_{0}) (\boldsymbol{u}\_{\rm des} - \boldsymbol{u}\_{0}) \right] - \left[ \dot{\mathbf{x}}\_{0} + \mathbf{G} (\mathbf{x}\_{0}, \boldsymbol{\mu}\_{0}) (\boldsymbol{u} - \boldsymbol{u}\_{0}) \right] \\ = \mathbf{G} (\mathbf{x}\_{0}, \boldsymbol{u}\_{0}) (\boldsymbol{u}\_{\rm des} - \boldsymbol{u}) \end{array} \tag{37}$$

where *des u* represents the desired actuator input, which can be produced by the rate con-

signal in real time, which requires an error *rat e* to be defined between the reference sig-

*rat des rat e y* = − ω

(38)

For the rate loop, it is expected that the helicopter can track the input angular rate

0 0 ( , )( ) *des Gx u u u* = − where *udes* represents the desired actuator input, which can be produced by the rate controller.

nal and the system output, yielding

troller.

*4.3. Rate Loop* 

## *4.3. Rate Loop*

For the rate loop, it is expected that the helicopter can track the input angular rate signal in real time, which requires an error *erat* to be defined between the reference signal and the system output, yielding

$$
\omega\_{\rm rat} = \omega\_{\rm des} - \mathcal{Y}\_{\rm rat} \tag{38}
$$

where *ωdes* is the reference command which can be produced by the attitude controller. The rate of angular change between the body frame to the North–East–Down (NED) co-ordinate system *yrat* = *x* = *ω<sup>h</sup>* = - *p q r<sup>T</sup>* represents the angular rate output of helicopter.

Parallel to the INDI design procedure, we differentiate the output expression to obtain its explicit dependence on the actuator inputs *u*. This yields

$$
\dot{y}\_{\rm rat} = \dot{\omega}\_{\rm l} = \int^{-1} (\mathfrak{m} - \omega \mathfrak{e}\_{\rm l} \times \mathfrak{J} \omega\_{\rm l}) \tag{39}
$$

(39) can be rewritten as

.

$$
\dot{y}\_{\text{rat}} = \dot{\omega}\_{\text{h}} = f\_{\text{rat}}(\mathbf{x}) + \mathbf{g}\_{\text{rat}}(\mathbf{x}, \boldsymbol{\mu}) \tag{40}
$$

where

$$f\_{rat}(\mathbf{x}) = f^{-1}[m\_{fus}(\mathbf{x}) + m\_{lt}(\mathbf{x}) + m\_{\upsilon l}(\mathbf{x}) - \omega \times I\omega] \tag{41}$$

$$\log\_{rat}(\mathbf{x}, \boldsymbol{\mu}) = f^{-1}[m\_{mr}(\mathbf{x}, \boldsymbol{\mu}) + m\_{tr}(\mathbf{x}, \boldsymbol{\mu})] \tag{42}$$

Note that the number of actuators of the helicopter is four, which is not equal to the output state number in the rate loop. This means a control allocation scheme must be used to deal with this overdetermined system. However, because the change of the collective pitch of the main rotor, denoted by *θ*0, is always accompanied by the alteration of total lift, the value of *θ*<sup>0</sup> can be determined by the velocity on the *z*-axis under the NED reference frame. Therefore, we separately give the control law of the main rotor in the subsequent design. Now, we define the other three actuators as *u* 0 = - *<sup>θ</sup>*1*<sup>s</sup> <sup>θ</sup>*1*<sup>c</sup> <sup>θ</sup>*0*tr<sup>T</sup>* .

Combining the analysis of the helicopter model before, we can obtain the rotational dynamics under actuator delay:

$$
\dot{y}\_{\text{rat}} = \dot{\omega}\_{\text{h}} = f\_{\text{rat}}(\mathbf{x}) + \mathbf{g}\_{\text{rat}}(\mathbf{x}, u\_{\text{\tau}}) \tag{43}
$$

Choosing the virtual control *νrat* = . *ωdes* + *kνreret* with the control gain matrix *kν<sup>r</sup>* ∈ *R* 3×3 , the controller *u* 0 can be given by

$$
\mu' = \mu'\_0 + \Delta\mu'\_n + \Delta\mu'\_c \tag{44}
$$

$$
\Delta u'\_{\,\,n} = G\_{\rm rot}(\chi\_0, \mu\_0)^{-1} (\nu - \dot{\omega}\_{\hbar\_0}) \tag{45}
$$

$$
\Delta u'\_{\varepsilon} = k\_{\text{tr}} G\_{\text{rat}}(\mathbf{x}\_{0}, u\_{0})^{-1} \int\_{t\_{0}}^{t} [\dot{e}\_{\text{rat}} + k\_{\text{rr}} e\_{\text{rat}} + g\_{\text{rat}}(\mathbf{x}, u\_{\text{\pi}}) - g\_{\text{rat}}(\mathbf{x}, u)] d\theta \tag{46}
$$

where *Grat*(*x*0, *u*0) = *J* −1 *∂ ∂u* 0 [*mmr*(*x*, *u*) + *mtr*(*x*, *u*)] *x*0,*u*<sup>0</sup> and *kur*, *krr* ∈ *R* <sup>3</sup>×<sup>3</sup> are also diagonal matrices. Note that there is a very complicated relationship between control input *u* and the moment generated by the main rotor and tail rotor because of the aerodynamics of the rotors. Hence, we adopt the central finite differences to calculate the control effectiveness matrix denoted by *D*, yielding

*D* = *∂ ∂u* 0 [*mmr*(*x*, *u*) + *mtr*(*x*, *u*)] *x*0,*u*<sup>0</sup> = *m<sup>T</sup> mr x*0,*u* 0 <sup>0</sup>+ h *δθ*1*<sup>s</sup>* 0 0i*<sup>T</sup>* <sup>−</sup>*m<sup>T</sup> mr x*0,*u* 0 0− h *δθ*1*<sup>s</sup>* 0 0i*<sup>T</sup>* <sup>2</sup>*δθ*1*<sup>s</sup> m<sup>T</sup> mr x*0,*u* 0 <sup>0</sup>+ h 0 *δθ*1*<sup>c</sup>* 0 i*T* <sup>−</sup>*m<sup>T</sup> mr x*0,*u* 0 0− h 0 *δθ*1*<sup>c</sup>* 0 i*T* <sup>2</sup>*δθ*1*<sup>c</sup> m<sup>T</sup> tr x*0,*u* 0 <sup>0</sup>+ h 0 0 *<sup>δ</sup>θ*0*tr* <sup>i</sup>*<sup>T</sup>* <sup>−</sup>*m<sup>T</sup> tr x*0,*u* 0 0− h 0 0 *<sup>δ</sup>θ*0*tr* <sup>i</sup>*<sup>T</sup>* <sup>2</sup>*δθ*0*tr* (47)

where *δθ*1*<sup>s</sup>* , *δθ*1*<sup>c</sup>* , and *δθ*1*<sup>c</sup>* are a small percent of each actuator input value.

For actuator dynamics, a pseudo-control hedge command is generated to provide the control system from trying to track the reference command when saturation occurs. According to (37), the pseudo-control hedge command *ν<sup>r</sup>* for the rate loop can be obtained by

$$\begin{array}{l} \nu\_r = (\dot{\mathbf{x}}\_0 + \mathbf{J}^{-1} D(u'\_{\rm des} - u'\_0)) - (\dot{\mathbf{x}}\_0 + \mathbf{J}^{-1} D(u' - u'\_0)) \\ = \mathbf{J}^{-1} D(u'\_{\rm des} - u') \end{array} \tag{48}$$

where *u* 0 *des* represents the desired input vector produced by the rate controller.

After completing this, the whole rate-loop control system is accomplished. However, the helicopter still faces stability issues, for the reason that the Euler angle in the attitude loop is not closed-loop stable.

## *4.4. Attitude Loop*

Then, for the attitude loop, we can use the NDI control on account of no model uncertainty existing here. In this loop, the system can be given by

$$
\dot{\theta} = \begin{bmatrix}
1 & \sin\phi \tan\theta & \cos\phi \tan\theta \\
0 & \cos\phi & -\sin\phi \\
0 & \sin\phi / \cos\theta & \cos\phi / \cos\theta
\end{bmatrix} \omega\_h \tag{49}
$$

$$y\_{att} = \mathbf{x} = \boldsymbol{\theta} = \begin{bmatrix} \boldsymbol{\phi} & \boldsymbol{\theta} & \boldsymbol{\psi} \end{bmatrix}^{T} \tag{50}$$

where *yatt* represents the attitude angle output of the helicopter.

Unlike the rate loop, there is no model uncertainty or time delay in the attitude loop. Therefore, the attitude controller only needs to convert the attitude angle tracking error into the desired angular rate command as the inner loop input signal through the NDI method. Considering the state Equation (49), the reference input signal of rate loop *ωre f* can be obtained by

$$
\omega\_{ref} = \begin{bmatrix}
1 & 0 & -\sin\theta \\
0 & \cos\phi & \sin\phi\cos\theta \\
0 & -\sin\phi & \cos\phi\cos\theta
\end{bmatrix} \nu\_{att} \tag{51}
$$

where *νatt* is the virtual control for the attitude loop and it can be given by the attitude tracking error *eatt* as .

$$\nu\_{att} = \theta\_{des} + k\_{\nu a} e\_{att} \tag{52}$$

in which *<sup>k</sup>ν<sup>a</sup>* <sup>∈</sup> <sup>R</sup>3×<sup>3</sup> is the control gain matrix and *θ des* is the attitude reference command for the helicopter. Note that the inverse of the transformation matrix always exists for detΩ = 1/ cos *θ* 6= 0.

In the attitude loop, the pseudo-control hedge command *ν<sup>a</sup>* is

$$\boldsymbol{\nu}\_{a} = \boldsymbol{\Omega} (\boldsymbol{\omega}\_{ref} - \boldsymbol{\omega}\_{h}) \tag{53}$$

## *4.5. Control Law for Collective Pitch of Main Rotor*

As mentioned in the rate loop, the operation of the collective pitch of the main rotor *θ*<sup>0</sup> will change the lift of the helicopter directly. Hence, the actuator *θ*<sup>0</sup> should be related to the vertical velocity of the helicopter. In the *z*-axis direction, the following equation is made

$$\dot{V}\_z - \mathbf{g} = \begin{bmatrix} -\sin\theta & \cos\theta\sin\phi & \cos\theta\cos\phi \end{bmatrix} \frac{F}{m} \tag{54}$$

where g is the gravitational acceleration and *F* = - *F<sup>x</sup> F<sup>y</sup> F<sup>z</sup> T* represents the total force in the three axes. It contains the contributions of all the main parts of the helicopter, yielding

$$F = F\_{fus} + F\_{ht} + F\_{vt} + F\_{mr} + F\_{tr} \tag{55}$$

Once again, *F* consists of the force produced by the fuselage, horizontal tail, vertical tail, main rotor, and tail rotor, respectively.

Note that, although the total force *F* contains many parts, the main part is the force generated by the main rotor since it resists gravity while allowing the helicopter to move flexibly. Based on the assumption above, we can obtain the direct dependence about . *V<sup>z</sup>* to the delayed input *θ*0*<sup>τ</sup>* as

$$\dot{V}\_z = g + \frac{1}{m} \begin{bmatrix} -\sin\theta\\ \cos\theta\sin\phi\\ \cos\theta\cos\phi \end{bmatrix}^T F\_{mr}(\mathbf{x}, u\_\tau) \tag{56}$$

Choosing the virtual control *νvz* = . *Vzdes* + *kνzevz* with the control gain constant *kν<sup>z</sup>* ∈ *R* <sup>+</sup>, the controller for *θ*<sup>0</sup> can be given by

$$
\theta\_0 = \theta\_{0,0} + \Delta\theta\_{0,\text{Il}} + \Delta\theta\_{0,\text{c}} \tag{57}
$$

$$
\Delta\theta\_{0,n} = m \left( \begin{bmatrix} -\sin\theta \\ \cos\theta \sin\phi \\ \cos\theta \cos\phi \end{bmatrix}^T H \right)^{-1} (\nu\_z - \dot{V}\_{z\_0}) \tag{58}
$$

$$\Delta\theta\_{0,\mathcal{E}} = k\_{\text{lz}}m \left( \begin{bmatrix} -\sin\theta\\ \cos\theta\sin\phi\\ \cos\theta\cos\phi \end{bmatrix}^T H \right)^{-1} \int\_{t\_0}^{t} \left[ \dot{e}\_z + k\_{\text{rz}}e\_z + F\_{\text{mr}}(\mathbf{x}, \boldsymbol{\mu}\_{\tau}) - F\_{\text{mr}}(\mathbf{x}, \boldsymbol{\mu}) \right] d\theta \tag{59}$$

where *θ*0,0 also represents the previous sampling value of *θ*0; *kuz*, *krz* ∈ *R* are positive control gains; and *H* is the control effectiveness matrix, which can be expressed as

$$H = \left. \frac{\partial F\_{mr}(\mathbf{x}, \boldsymbol{\mu})}{\partial \theta\_0} \right|\_{\mathbf{x}\_0, \boldsymbol{\mu}\_0} \tag{60}$$

It can also be calculated by using the central finite difference as

$$H = \left[\frac{F\_{mr}^{\,\,T}(\mathbf{x}\_{0\prime}\theta\_{0,0} + \delta\_{\theta\_0}) - F\_{mr}^{\,\,T}(\mathbf{x}\_{0\prime}\theta\_{0,0} - \delta\_{\theta\_0})}{2\delta\_{\theta\_0}}\right]^T \tag{61}$$

where *δθ*<sup>0</sup> is a small percent of *θ*0.

Now, the helicopter attitude control system under the multiple constraints of the actuators has been designed, which is shown in Figure 2.

tuators has been designed, which is shown in Figure 2.

is a small percent of

**Figure 2.** Schematic of the attitude control system based on the proposed method. **Figure 2.** Schematic of the attitude control system based on the proposed method.

### **5. Simulation Results 5. Simulation Results**

where

θ 0 δ

In this section, several simulations will be given to verify the advantages of the proposed control law. We will simulate the attitude control of the helicopter in a hovering state. The model uncertainties are given as , , 0.1 , 0.1 , *L LL L tr tr C CC C* α αα α Δ =− Δ =− , , 0.5 , *L L ht ht C C* α α Δ = and , , 0.5 *L L vt vt C C* α α Δ = , where *CL*α , *<sup>L</sup>* ,*tr C* α , *<sup>L</sup>* ,*ht C* α , and *<sup>L</sup>* ,*vt C* α are the lift curve slope of the blade of the main rotor, the blade of the tail rotor, the horizontal tail, and the vertical tail, respectively. The delay time is τ =100*ms* , initial input In this section, several simulations will be given to verify the advantages of the proposed control law. We will simulate the attitude control of the helicopter in a hovering state. The model uncertainties are given as ∆*CL<sup>α</sup>* = −0.1*CL<sup>α</sup>* , ∆*CLα*,*tr* = −0.1*CLα*,*tr*, ∆*CLα*,*ht* = 0.5*CLα*,*ht* , and ∆*CLα*,*vt* = 0.5*CLα*,*vt* , where *CL<sup>α</sup>* , *CLα*,*tr*, *CLα*,*ht* , and *CLα*,*vt* are the lift curve slope of the blade of the main rotor, the blade of the tail rotor, the horizontal tail, and the vertical tail, respectively. The delay time is *τ* = 100*ms*, initial input *u*<sup>0</sup> = - 0.2484 0.0275 <sup>−</sup>0.0135 0.1289*<sup>T</sup>* , diagonal element of control gain matrix *kur* is - 1 1 1.5*<sup>T</sup>* , and *krr* is - 2 2 3*<sup>T</sup>* .

<sup>0</sup> <sup>0</sup> 0 ,

0

θ

δ

2 *<sup>T</sup> T T Fx Fx mr mr <sup>H</sup>* θ

Now, the helicopter attitude control system under the multiple constraints of the ac-

 +− − <sup>=</sup> 

0 0,0 0 0,0 (, ) (, )

*x u*

0 0

 θδ

<sup>∂</sup> <sup>=</sup> (60)

 θ

(61)

(,) *mr*

*F xu <sup>H</sup>* θ

∂

θδ

It can also be calculated by using the central finite difference as

θ0 .

<sup>0</sup> [ ] 0.2484 0.0275 0.0135 0.1289 *<sup>T</sup> <sup>u</sup>* = − , diagonal element of control gain matrix *ur k* is [ ] 1 1 1.5 *<sup>T</sup>* , and *rr k* is [ ] <sup>223</sup> *<sup>T</sup>* . The delay of the actuators will degrade the tracking performance and increase the control effort provided by the actuators, which always leads to state overshoot and actu-The delay of the actuators will degrade the tracking performance and increase the control effort provided by the actuators, which always leads to state overshoot and actuator saturation, and even causes the system to become unstable when the delay time gradually increases. This phenomenon can be observed in Figure 3. *Aerospace* **2023**, *10*, x FOR PEER REVIEW 13 of 18

**Figure 3.** The response of the INDI-based control system under different actuators delay**. Figure 3.** The response of the INDI-based control system under different actuators delay.

1s(deg) 1c(deg)In this simulation, the operating frequency of the control system is set at 100 Hz. When the actuator delay is 50 ms, the INDI-based control system can barely maintain its tracking performance. A small range of oscillation appears in response when the delay time is 100 ms. However, actuators need to change frequency to maintain this steady state in the INDI scheme, which is shown in Figure 4. This can be understood as, when time delay exists, the control input does not correspond to the current input error, and it needs to be adjusted constantly within the whole time delay. In the case of the delay of 150 ms, the

In the next simulation, the delay time between the controller and actuators is set at

100 ms. In Figure 5, it can be seen that, when the delay compensation is applied, the state overshoot is significantly reduced and the system**'**s rapidity is also improved compared to the original response. Figure 6 shows that, in addition to the need for a larger control

effect, the original phenomenon of rapid changes has disappeared.

**Figure 4.** Actuator input of helicopter.

0tr(deg)

(deg)

(deg)

(deg)

system response has already experienced a relatively large oscillation, and its dynamic characteristics are seriously degraded. **Figure 3.** The response of the INDI-based control system under different actuators delay**.** 

**Figure 3.** The response of the INDI-based control system under different actuators delay**.** 

*Aerospace* **2023**, *10*, x FOR PEER REVIEW 13 of 18

*Aerospace* **2023**, *10*, x FOR PEER REVIEW 13 of 18

**Figure 4.** Actuator input of helicopter. **Figure 4.** Actuator input of helicopter. In the next simulation, the delay time between the controller and actuators is set at

In the next simulation, the delay time between the controller and actuators is set at 100 ms. In Figure 5, it can be seen that, when the delay compensation is applied, the state overshoot is significantly reduced and the system**'**s rapidity is also improved compared to the original response. Figure 6 shows that, in addition to the need for a larger control In the next simulation, the delay time between the controller and actuators is set at 100 ms. In Figure 5, it can be seen that, when the delay compensation is applied, the state overshoot is significantly reduced and the system's rapidity is also improved compared to the original response. Figure 6 shows that, in addition to the need for a larger control effect, the original phenomenon of rapid changes has disappeared. 100 ms. In Figure 5, it can be seen that, when the delay compensation is applied, the state overshoot is significantly reduced and the system**'**s rapidity is also improved compared to the original response. Figure 6 shows that, in addition to the need for a larger control effect, the original phenomenon of rapid changes has disappeared.

**Figure 5.** The response after delay compensation.

On the basis of delay compensation, we carry out the PCH design for the system. The advantages of PCH can be seen in Figure 7, which not only eliminates the overshoot, but also reduces the 0.7 s settling time of the system. Figure 8 is also given here to compare the changes of the actuators between the three cases.

also reduces the 0.7 s settling time of the system. Figure 8 is also given here to compare

**Figure 6.** Actuator input after delay compensation**. Figure 6.** Actuator input after delay compensation. the changes of the actuators between the three cases.

**Figure 5.** The response after delay compensation**.** 

**Figure 5.** The response after delay compensation**.** 

*Aerospace* **2023**, *10*, x FOR PEER REVIEW 14 of 18

**Figure 7.** The response after the introduction of PCH**. Figure 7.** The response after the introduction of PCH.

Finally, the response of the three channels of the attitude angle under the proposed

change regularly corresponding to the tracking of the three Euler angular rates. Furthermore, the system can also track the command signal well when the actuator input is satu-

**Figure 8.** Actuator input after PCH**. Figure 8.** Actuator input after PCH.

**Figure 9.** Three Euler angle responses based on the proposed method**.** 

rated and delayed, which can be observed in Figure 10.

Finally, the response of the three channels of the attitude angle under the proposed control scheme is given in Figure 9. It can be seen that the three Euler angles can be decoupled and show good dynamic characteristics. In the Figure 10, the angular rates also change regularly corresponding to the tracking of the three Euler angular rates. Furthermore, the system can also track the command signal well when the actuator input is saturated and delayed, which can be observed in Figure 10. control scheme is given in Figure 9. It can be seen that the three Euler angles can be decoupled and show good dynamic characteristics. In the Figure 10, the angular rates also change regularly corresponding to the tracking of the three Euler angular rates. Furthermore, the system can also track the command signal well when the actuator input is saturated and delayed, which can be observed in Figure 10.

0 5 10 15 20 25

0 5 10 15 20 25

0 5 10 15 20 25 Time(deg)

Delay compensation PCH-Delay compensation

Finally, the response of the three channels of the attitude angle under the proposed

**Figure 8.** Actuator input after PCH**.** 


0

10

0tr(deg)

20


1c(deg)

0

1s(deg)

*Aerospace* **2023**, *10*, x FOR PEER REVIEW 15 of 18

**Figure 9.** Three Euler angle responses based on the proposed method**. Figure 9.** Three Euler angle responses based on the proposed method.

**Figure 10.** Actuator input corresponding to the responses in Figure 8**. Figure 10.** Actuator input corresponding to the responses in Figure 8.

### **6. Conclusions 6. Conclusions**

script.

**References** 

https://doi.org/10.2514/2.5039.

https://doi.org/10.1155/2017/5309403.

In this paper, an INDI-based actuator compensation attitude controller is developed for the helicopter subject to the time delay, position, and rate saturations in actuators. The controller is composed of a rate controller which ensures the rate performance of the helicopter, an attitude controller which guarantees the attitude tracking performance, and a collective pitch controller which meets the needs of the vertical changes in the *z* -axis direction. The model reduction method is used to design an INDI-based controller for the rate loop and collective pitch of the main rotor, which improves the robustness to the timevarying actuator delay. Meanwhile, the PCH technique is introduced for both the rate loop In this paper, an INDI-based actuator compensation attitude controller is developed for the helicopter subject to the time delay, position, and rate saturations in actuators. The controller is composed of a rate controller which ensures the rate performance of the helicopter, an attitude controller which guarantees the attitude tracking performance, and a collective pitch controller which meets the needs of the vertical changes in the *z*-axis direction. The model reduction method is used to design an INDI-based controller for the rate loop and collective pitch of the main rotor, which improves the robustness to the time-varying actuator delay. Meanwhile, the PCH technique is introduced for both the

and attitude loop to provide a filter such that the following commands hold within the capability of the controllers. Finally, simulations demonstrate the effectiveness and ro-

copter will be designed and a maneuver flight simulation analysis that meets the require-

**Author Contributions:** Conceptualization, S.Z., H.Z. and K.J.; methodology, S.Z. and H.Z.; formal analysis, H.Z.; investigation, H.Z. and K.J.; resources, H.Z. and K.J.; data curation, H.Z.; writing original draft preparation, H.Z.; writing—review and editing, H.Z. and K.J.; supervision, S.Z.; project administration, S.Z. All authors have read and agreed to the published version of the manu-

**Funding:** This research has received funding from the Natural Science Foundation of Jiangsu Province under Grant BK20201291, Aeronautical Science Foundation of China under Grant 201957052002, and Research and Practice Innovation Program of Nanjing University of Aeronautics

**Conflicts of Interest:** The authors declare no conflict of interest.

1. Padfield, G.D. *Helicopter Flight Dynamics: The Theory and Application of Flying Qualities and Simulation Modelling*, 2nd ed.; Ameri-

2. Prouty, R.W.; Curtiss, H.C. Helicopter Control Systems: A History. *J. Guid. Control. Dyn.* **2003**, *26*, 12–18.

3. Hu, J.; Gu, H. Survey on Flight Control Technology for Large-Scale Helicopter. *Int. J. Aerosp. Eng.* **2017**, *2017*, 5309403.

4. Kim, N.; Calise, A.J.; Hovakimyan, N.; Prasad, J.V.R.; Corban, E. Adaptive Output Feedback for High-Bandwidth Flight Control.

ments of ADS-33E-PRF will be performed.

and Astronautics under Grant xcxjh20210317. **Data Availability Statement:** Not applicable.

can Institute of Aeronautics and Astronautics: Reston, VA, USA, 2007; Volume 113.

*J. Guid. Control. Dyn.* **2002**, *25*, 993–1002. https://doi.org/10.2514/2.4985.

rate loop and attitude loop to provide a filter such that the following commands hold within the capability of the controllers. Finally, simulations demonstrate the effectiveness and robustness of the proposed controller. In the future, the outer loop controller for the helicopter will be designed and a maneuver flight simulation analysis that meets the requirements of ADS-33E-PRF will be performed.

**Author Contributions:** Conceptualization, S.Z., H.Z. and K.J.; methodology, S.Z. and H.Z.; formal analysis, H.Z.; investigation, H.Z. and K.J.; resources, H.Z. and K.J.; data curation, H.Z.; writing original draft preparation, H.Z.; writing—review and editing, H.Z. and K.J.; supervision, S.Z.; project administration, S.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research has received funding from the Natural Science Foundation of Jiangsu Province under Grant BK20201291, Aeronautical Science Foundation of China under Grant 201957052002, and Research and Practice Innovation Program of Nanjing University of Aeronautics and Astronautics under Grant xcxjh20210317.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

## **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

## *Article* **An Improved Fault Identification Method for Electromechanical Actuators**

**Gaetano Quattrocchi , Pier C. Berri , Matteo D. L. Dalla Vedova \* and Paolo Maggiore**

Department of Aerospace and Mechanical Engineering (DIMEAS), Politecnico di Torino, 10129 Turin, Italy; gaetano.quattrocchi@polito.it (G.Q.); pier.berri@polito.it (P.C.B.); paolo.maggiore@polito.it (P.M.) **\*** Correspondence: matteo.dallavedova@polito.it

**Abstract:** Adoption of electromechanical actuation systems in aerospace is increasing, and so reliable diagnostic and prognostics schemes are required to ensure safe operations, especially in key, safetycritical systems such as primary flight controls. Furthermore, the use of prognostics methods can increase the system availability during the life cycle and thus reduce costs if implemented in a predictive maintenance framework. In this work, an improvement of an already presented algorithm will be introduced, whose scope is to predict the actual degradation state of a motor in an electromechanical actuator, also providing a temperature estimation. This objective is achieved by using a properly processed back-electromotive force signal and a simple feed-forward neural network. Good prediction of the motor health status is achieved with a small degree of inaccuracy.

**Keywords:** prognostics; electromechanical actuators; neural network; temperature

## **1. Introduction**

Electromechanical actuators (EMA) are widely used in industry, and in recent decades are seeing increasing adoption in the aviation sector, especially when the more electric or all-electric [1] design philosophies are adopted, which seeks to reduce secondary power type usage, using mostly or exclusively electrical power for systems actuation [2].

An interesting application of the more electric philosophy is in the Airbus A380 where two electrical lines are used as backup sources of power in case of loss of the main hydraulic lines. However, electrohydrostatic actuators (EHAs) are used rather than EMAs, since one of EMAs' most common failure modes is jamming [3]. Jamming is characterized by a sudden actuator stop in a certain position, thus locking the flight control surface position, creating a dangerous moments imbalance and thus possibly uncontrollable yaw, pitch or roll. Furthermore, the estimation of jamming probability is harder than for EHAs, where the knowledge regarding current hydraulic actuators can be used [4].

For this reason, as of today the use of EMAs is still mainly limited to non-safety-critical systems, such as secondary flight systems, high-lift devices or airbrakes. However, EMAs have some merits that make them preferable over other architecture, such as complexity, weight and maintenance requirements, especially in low-power applications [5]. In fact, the simplest method for providing mechanical power using an electrical supply is an electromechanical actuator, which in its most basic form is an electrical motor with a mechanism converting rotary motion to linear (using for example a ball screw). Usually, in aviation, given the high actuation torques needed and the volume constraints, some form of a reducer is incorporated between the motor and the rotation-linear converter to multiply the motor torque. A schematic view of an EMA is shown in Figure 1.

Another important aspect to consider is the lack of extensive failure datasets in operating conditions, which further discourages the adoption of EMAs given the severity of some failure modes [6]. Even though EMAs, as previously stated, are already used in secondary flight control actuation, data obtained for this use case are not directly transferable for primary flight control due to the very different operating nature of the two applications.

**Citation:** Quattrocchi, G.; Berri, P.C.; Dalla Vedova, M.D.L.; Maggiore, P. An Improved Fault Identification Method for Electromechanical Actuators. *Aerospace* **2022**, *9*, 341. https://doi.org/10.3390/ aerospace9070341

Academic Editor: Gianpietro Di Rito

Received: 29 April 2022 Accepted: 23 June 2022 Published: 25 June 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

**Figure 1.** Block schematics of an EMA.

To increase safety and reliability, Prognostics and Health Management (PHM) could prove beneficial. The traditional definition of PHM can be found in [7], where it is defined as the ability to evaluate the current health status and predict the future behavior of a system, using knowledge of the current state; additionally, maintenance can be planned accordingly to properly maintain the system.

Several approaches to PHM exist and can be classified into three different categories: model-based, hybrid or data-driven, depending on the degree of knowledge of the system studied and the availability of data.

Model-based PHM uses a representation of the system (or component) of interest, using a set of differential Equations [8–10]. If properly modeled, the system status can be assessed with good accuracy and fault propagation is generally embedded in the fault model [11,12]. The main drawback of model-based approaches is the difficulty in creating accurate and representative models for complex systems; in particular in aeronautics, it is hard to model, or present in a simplified way, complex interactions such as aerodynamics or mechanical backlash.

The opposite approach is data-driven, where the system is treated as a black box, i.e., without detailed knowledge. In this case, machine learning methods are widely used to detect variations in the system's health status. One of the most used tools is Neural Networks of various kinds, such as in [13–16]. Some other tools used are combinations of filters and neural networks, as in [17,18], Bayesian inference [19] and Markov chains [20]. The main weakness of data-driven approaches is the need for a large dataset, which is generally not readily available, especially if run-to-failure of complex and expensive systems are needed.

Hybrid approaches use both model-based and data-driven techniques. The idea is to create a physical model of the system, extract relevant health indicators using physical quantities, and then use machine learning techniques to obtain an estimation of the actual system health status. Examples can be found in [21–24].

In this work, a hybrid prognostic method to evaluate the damage of a BLDC is presented, using the back-EMF signal as a prognostic indicator. In particular, the paper will describe an evolution of the algorithm presented in [25], which has been improved and can now also give a temperature estimate of the motor phases, besides the damage to individual stator coils and the static eccentricity of the rotor.

## **2. Materials and Methods**

## *2.1. Scope of the Work and Novelties Introduced*

The scope of this work is to present an evolution of the prognostic algorithm presented in [25]. The use of the back-EMF coefficient is supported by the fact that it is sensitive to the faults of interest (partial inter-turn shorts, static eccentricity), while not being affected by the command or load imposed on the actuator, as demonstrated in [26].

As previously stated, in this work the effect of temperature on phase resistance is also considered. The inclusion of temperature dependency requires a new approach to the evaluation of system status.

The algorithm can be described as follows:


The first two steps are self-explanatory; the faults are modeled as in [25], using a fault vector *f* = [*Na*, *N<sup>b</sup>* , *Nc*, *ξx*, *ξy*, ∆*T*], with 6 different components: *Na*, *N<sup>b</sup>* , *N<sup>c</sup>* which are the fraction of turns shorted for each phase; *ξx*, *ξ<sup>y</sup>* which are the components of static eccentricity in cartesian coordinates (in [25] polar coordinates are used) and finally ∆*T* which is the temperature deviation from reference conditions *T*<sup>0</sup> = 20 ◦C. The variation ranges for the variables are the same as those in [25]. For stator windings temperature, the range ∆*T* = [−50; + 70] ◦C is chosen, since it is representative of conditions encountered in aeronautical applications. However, the range chosen might be too restrictive to include sudden transitory temperature spikes and might need adjustments in future developments. The faults are injected into the model prior to each simulation.

In the third step, the estimation of the back-EMF coefficient is obtained using the following equation:

$$V - k'\omega = V - (k + k\_\varepsilon)\omega = R\_0 i \tag{1}$$

where *k* 0 is the estimated back-EMF coefficient and *k<sup>e</sup>* is the estimation error. In this case, values of *V*, *ω*, *i* are those that can be measured from the simulation (and by extension, from a real system), while *R*<sup>0</sup> is the nominal resistance of a phase. There is an error in the estimation of the back-EMF coefficient since the nominal resistance is used, thus not considering the effect of both a temperature variation and partial phase to phase short, which changes this value.

The actual (or true) system condition is described by:

$$V - k\omega = R \cdot i = (R\_0 + \Delta R)i \tag{2}$$

where *k* is the true back-EMF coefficient, *R* = (*R*<sup>0</sup> + ∆*R*) is the true resistance, i.e., the effective resistance of the coil in the instantaneous conditions of temperature, fault, etc. and ∆*R* is the deviation of actual resistance from nominal value.

Now, subtracting Equation (2) from Equation (1) and rearranging, one can obtain:

$$k\_{\ell} = \Delta R \frac{i}{\omega} \tag{3}$$

Assuming that ∆*R* is constant (i.e., *R* varies slowly), which is reasonable in the framework presented, since each measurement is very short (in the order of one second), derivating Equation (3) we have:

$$\frac{\partial k\_{\ell}}{\partial (i/\omega)} = \Delta R \tag{4}$$

Furthermore, assuming that *k* is constant, which implies that the fault does not change during the simulation, we can obtain:

$$\frac{\partial(k+k\_{\ell})}{\partial(i/\omega)} = \frac{\partial k'}{\partial(i/\omega)} = \Delta R \tag{5}$$

So, we have demonstrated that *k* 0 = (*k* + *ke*) is linearly dependent to *i*/*ω* with a slope equal to ∆*R*. Now, to obtain the real values of back-EMF coefficient and resistance, we have to iteratively reconstruct the value of *k* (Equation (2)), using a temporary *R* ∗ variable to optimize to make Equation (5) equal to zero. This stems from the definition of back-EMF:

$$|BEMF| = \frac{\partial \Phi}{\partial t} = \frac{\partial \Phi}{\partial \theta}\frac{\partial \theta}{\partial t} = \frac{\partial \Phi}{\partial \theta}\omega \quad \rightarrow \quad k\_j(n\_j, \theta, T) = \frac{\partial \Phi(n\_j, \theta, T)}{\partial \theta} \tag{6}$$

where Φ is the concatenated magnetic flux, *n* is the number of shorted turns and *T* is the temperature. In other words, the concatenated magnetic flux is a function of angle, number of turns and temperature and so is the back-EMF coefficient *k<sup>j</sup>* . It does depend only on motor geometry and on the magnetic properties of the magnets, and thus the temperature dependency. In this preliminary work, the temperature-induced variation of magnetic properties is not considered and will be added in further developments.

The following problem must then be solved:

$$\begin{cases} \quad V - k^\* \omega = R^\* \cdot i \\ \quad \frac{\partial (k^\*)}{\partial (i/\omega)} = 0 \end{cases} \tag{7}$$

At convergence, we obtain *R* ∗ = *R* which implies *k<sup>e</sup>* = 0 (from Equation (2)): we are calculating true resistance and true back-EMF coefficient. In this work, a simple bisection method is used to perform the optimization.

Up to this point, we have considered *kj*(*θ*, *i*/*ω*), as visible in Figure 2, where the subscript *j* indicates one of the three phases; in this case, the number of variables to optimize will be 3 · *n*, where *n* is the number of subdivisions along the *θ* axis. The total number of variables will thus be 6 · *n*, i.e., 3*n* values of resistance and 3*n* values of back-EMF coefficient.

However, in order to simplify the computation, the dependence on *θ* has been dropped, thus collapsing the 3D graph into a 2D plot of *kj*(*i*/*ω*), i.e., Figure 3. Now, a 'global' (or *generalized*) *k<sup>j</sup>* approximation can be calculated, using least square fit. Values close to zero have been discarded, since they provide no additional information, and an absolute value on *k<sup>j</sup>* has been applied. The final result is a reduction of the number of variables from 6 · *n* to 6, i.e., 3 generalized resistances (one for each phase) and 3 generalized back-EMF coefficients.

These 6 values are used in a simple feed-forward neural network to perform an estimation of the fault vector *f* .

Figures 2 and 3 have been obtained using a parabolic position command (i.e., a speed ramp) with a constant angular acceleration of 0.3 rad/s<sup>2</sup> ; initial conditions are zero angular position and zero angular velocity (*θ* = ˙*θ* = 0). The following fault vector was seeded: *<sup>f</sup>* = [0.0375, 0.0504, 0.0507, 0.0059, 2.8 <sup>×</sup> <sup>10</sup>−<sup>5</sup> , 4].

**Figure 2.** Reconstructed *ka* as function of both *θ* and *i*/*ω*.

**Figure 3.** Reconstructed *ka* as function of *i*/*ω* (non optimized).

## *2.2. Brief System Overview*

As previously stated, the system used to test the algorithm is the same as that used in [25,26], with minor modifications (Figure 4); it represents an electromechanical actuator acting on a flight surface.

The model is a high fidelity representation of a trapezoidal EMA acting on secondary flight control. It is a detailed, component-level model with very limited assumptions, i.e., lumped parameters. Many non-linear phenomena are modeled, including but not limited to electronic noise, dry friction and current and speed saturations. The model has been validated on literature data as reported in [27,28]. For further details on the model, please refer to [26].

With respect to previous iterations, the main enhancement is the implementation of a temperature dependency on phase resistance (each of the three RL branches in Figure 5), using the classical equation:

$$R = R\_{ref}(1 + \alpha \Delta T) \tag{8}$$

where *α* = 0.00404 ◦C −1 is the resistance temperature coefficient of copper (the materials of which the coils are made), *Rre f* is the reference temperature resistance (in this case *Tre f* = 20 ◦C) and ∆*T* is the temperature difference from *Tre f* .

**Figure 5.** Phases model detail (a subsystem of 'BLDC electromagnetic model').

Furthermore, the external load is set to zero, supposing an actuation during ground operations, while the actuation command is now a parabolic position command (i.e., a velocity ramp, or constant acceleration), since using this type of command better covers the (*i*/*ω*) space and thus embed more information in the same simulation time as opposed to more classical commands such as position ramps or steps.

Some examples of the data obtained after simulations can be found in Figures 6 and 7. It can be noted how the first ∼40 ms present strong fluctuations in the current and voltages

values (but also angular velocity): this is due to the strong non-linear effects modeled in the system, e.g., dry friction. For this reason, the first ∼50 ms of each simulation are discarded before applying the algorithm presented.

**Figure 6.** Graphs of phase current for different fault conditions.

**Figure 7.** Graphs of phase voltage for different fault conditions.

## *2.3. Dataset Used*

The number of different fault conditions simulated is 600; the dataset has been randomly divided into 70%, 15%, 15% subsets for training, validation and testing, respectively.

Dataset size has been empirically chosen to be a good representation of the 6-dimensional fault space; the dataset, regarding the first 5 variables of the fault vector (*Na*, *N<sup>b</sup>* , *Nc*, *ξx*, *ξy*) is the same presented in [25], with the two eccentricity components now transformed into cartesian coordinates from polar. For each fault vector, a temperature difference value (∆*T*), randomly sampled from the interval considered, has been appended.

Using an AMD 5600X, each simulation takes about 50 s to run; parallel pooling has been used to reduce the total dataset generation time, combined with Simulink Accelerator mode. For this study, a full software workflow has been used. In other words, the full system is simulated and relevant data logged. The algorithm presented is then applied and neural network outputs are compared to the fault vectors injected in the model before simulation, which are considered ground truth.

In the future, we hope to implement the faults considered in this study on a real test bench, even though it is a complex task, especially for the partial phase shorts.

## **3. Results**

In this section, the network used for fault regression will be described. Several tests on different hyperparameter combinations have been carried out, and the best performing set has been used for error distribution calculation. The network architecture is very similar to that presented in [25], with minor modifications to the values of the hyperparameters.

The architecture chosen is a feed-forward neural network with a single hidden layer of size 12 (Figure 8); the training function uses the Levenberg–Marquardt algorithm; the activation function for each neuron is the hyperbolic tangent sigmoid. The maximum size of failed validation checks is set to 10.

As expected, the network inputs vector is of size 6—including 3 generalized back-EMF coefficients and 3 generalized resistances, while the output is again of size 6 and is the fault vector used to generate the simulation.

**Figure 8.** Neural network topology.

In Figure 9, the mean absolute error box plot for each variable is shown. For visualization purposes, a new subset of 10 rows is used and called an external validation set. This dataset has never been used during training of the neural network so it is a good representation of the predictive capabilities of the network with new data during operations. As visible, the mean absolute error is very small, in the order of 0.02 on normalized data. The error distribution is however uneven between variables and this might be caused by the relatively small dataset used in training.

In absolute terms, the mean error for the ∆*T* estimation, on the external verification dataset, is about 1.8 ◦C, which is a good result.

**Figure 9.** Mean absolute error for external validation set.

## **4. Discussion**

As described in the previous section, the prediction accuracy of the system is promising. Even though the problem complexity has increased with respect to the previous algorithm [25], the neural network is capable of predicting the system status with adequate accuracy.

The network regression task is probably aided by the extensive pre-processing applied to raw data before being used as an input dataset, including filtering. In a real system, properly calibrated filtering will be needed to smooth out high-frequency noise in the signals of interest, before the application of the algorithm is presented.

Furthermore, the data compression techniques described in Section 2.1 are helpful in reducing the dimensions of the regression problem, even though the estimation obtained is a generalization of the health status rather than an estimation for each possible angular position.

## **5. Conclusions**

In this work, an improved prognostic algorithm for EMA has been presented. The strong points of this paper include the ability to estimate the current health status of the motor in terms of fault variables including partial phase shorts, static eccentricity and temperature deviation from ambient conditions. Furthermore, no additional sensors besides those needed for normal operations are required. Motor damage estimation is carried out using a feed-forward neural network after the application of the algorithm to raw data.

As with any other work that includes neural networks in the pipeline, a parametric optimization of the network could yield additional benefits in the form of higher accuracy or a simpler network if the algorithm is to be implemented on embedded hardware.

Furthermore, several assumptions have been made to simplify the algorithm, mainly the temperature independence of magnetic properties. Even though it is a reasonable assumption for small variations of temperature (magnetic flux variation of ca. −4% for 100 ◦C for SaCo magnets), a generalized algorithm should include and simulate such variations.

An increase in the simulation dataset size, after considering what the best combination of actuation command and load is, could prove beneficial in further increasing the network accuracy and robustness.

Finally, an empirical validation on properly calibrated equipment is mandatory to test that the assumptions made are reasonable and can effectively represent the system status.

**Author Contributions:** Conceptualization, P.C.B.; methodology, G.Q. and P.C.B.; software, G.Q.; validation, G.Q.; formal analysis, G.Q. and P.C.B.; investigation, G.Q.; resources, G.Q.; data curation, G.Q.; writing—original draft preparation, G.Q.; writing—review and editing, G.Q. and M.D.L.D.V. visualization, G.Q.; supervision, M.D.L.D.V. and P.M.; project administration, M.D.L.D.V. and P.M.; funding acquisition, M.D.L.D.V. and P.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The datasets needed to replicate what was presented in this study are openly available in FigShare at doi:10.6084/m9.figshare.19635378, accessed on 24 June 2022.

**Conflicts of Interest:** The authors declare no conflict of interest.

## **Abbreviations**

The following abbreviations are used in this manuscript:

BEMF Back Electro-Motive Force;


## **References**


## *Article* **Reliability-Oriented Configuration Optimization of More Electrical Control Systems**

**Zirui Liao 1,2,3, Shaoping Wang 1,2, Jian Shi 1,2,\*, Dong Liu 1,2 and Rentong Chen 1,2**


**Abstract:** More electrical vehicles adopt dissimilar redundant control systems with dissimilar power supplies and dissimilar actuators to achieve high reliability and safety, but this introduces more intricacy into the configuration design. Currently, it is difficult to identify the optimum configuration via the conventional trial-and-error approach within an acceptable timeframe. Hence, it is imperative to discover novel methods for the configuration design of more electrical vehicles. This paper introduced the design specification of more electric vehicles and investigated the contribution of different kinds of actuators, presenting a new multi-objective configuration optimization approach on the foundation of system reliability, weight, power, and cost. By adopting the non-dominated sorting genetic algorithm-II (NSGA-II), the Pareto optimization design set was obtained. Then, the analytic hierarchy process (AHP) was introduced to make a comprehensive decision on the schemes in the Pareto set and determine the optimal system configuration. Eventually, numerical results indicated that the reliability of our designed configuration increased by 5.89% and 55.34%, respectively, compared with dual redundancies and single redundancy configurations, which verified the effectiveness and practicability of the proposed method.

**Keywords:** more electric vehicles; dissimilar redundant actuation system; NSGA-II algorithm; optimization design

## **1. Introduction**

Safety critical systems in aircraft [1], carrier rockets [2], ships [3] and large machines [4] usually adopt redundant actuation systems to guarantee high reliability and safety [5]. In such designs, the rest systems can complete the work when one or more actuators fail. To reduce the number of common-cause faults and common-mode faults, dissimilar redundancy systems are adopted consequentially [6], in which different actuators with the same performances are used. With the rapid development of more electrical technologies, and more electric vehicles based on different power supply systems, actuation systems and computers have become increasingly common. In terms of actuation systems, hydraulic actuators (HAs), electro-hydrostatic actuators (EHAs) and electromechanical actuators (EMAs) [7] are leveraged in more electric vehicles. For instance, more electrical aircrafts are adopting heterogeneous actuation systems to move control surfaces based on HAs and EHAs [8]. Since HAs utilize hydraulic power supply systems while EHAs and EMAs adopt electrical power supply systems, the configuration design of more electric vehicles (MEVs) faces substantial challenges. The first one is how to design more electric control systems to maintain high mission reliability without permissive cost and weight. The second one is how to strike a balance between power supply systems and actuator distribution. The last one is how to determine the types of actuators for more electric control systems,

**Citation:** Liao, Z.; Wang, S.; Shi, J.; Liu, D.; Chen, R. Reliability-Oriented Configuration Optimization of More Electrical Control Systems. *Aerospace* **2022**, *9*, 85. https://doi.org/ 10.3390/aerospace9020085

Academic Editor: Gianpietro Di Rito

Received: 24 November 2021 Accepted: 28 January 2022 Published: 6 February 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

which include HAs, EHAs, and EMAs. Given the aforementioned factors, the design of safety-critical system configurations becomes a daunting task. Therefore, an appropriate and effective approach should be adopted to solve the combination-explosion problem of system configurations, thus achieving the optimization design of a more electric control system (MECS) [9,10].

The optimization design of an electric control system needs to realize multi-objective optimization and comprehensively consider various contradictory indicators. Traditional methods include the weighting method [11,12], the constraint method [13,14], and goal programming [15–17]. However, these methods require the designer to master corresponding background knowledge to determine weights, energy dissipation, and other indicators, and the explosive configuration is another challenge that arises from those methods. In addition, the program needs to be run independently many times, which may lead to inconsistent results each time and render it more difficult for designers to make eventual decisions. In recent decades, various multi-objective intelligent optimization algorithms have emerged and vigorously developed, for example, the genetic algorithm (GA) [18], differential evolution algorithm (DE) [19], particle swarm optimization algorithm (PSO) [20,21], and non-dominated sorting genetic algorithm (NSGA) [22]. Among them, NSGA is famous for its good performance due to the adoption of non-dominated sorting and fitness-sharing strategy. The former increases the possibility of superior characteristics inherited by the next generation, while the latter sustains population diversity, overcoming the overmuch reproduction of super individuals and preventing prematurity convergence [23]. Nevertheless, there are also disadvantages such as high computational complexity in NSGA. Thus, in 2002, Deb et al. proposed NSGA-II [24]. In contrast to NSGA, the rapid non-dominated sorting method, elitism preservation strategy and congestion comparison were introduced in NSGA-II, which greatly reduced the computational complexity of the algorithm. Moreover, NSGA-II also expanded the distribution space of the solution set in the Pareto frontier, thus maintaining population diversity [25]. At present, NSGA-II has been widely applied to tackle multi-objective optimization problems in various fields. Xia et al. [26] investigated the multi-objective optimization problem for AUV conceptual design, where NSGA-II was applied to find the optimal Pareto frontier, with a comprehensive consideration of cost, effectiveness, and risks. The result verified the effectiveness of the algorithm. Alam et al. [27] studied the problem of AUV design and construction and employed NSGA-II to determine the optimum design of a torpedo-like AUV with a total length of 1.3 m. In addition, a kind of heavier-than-water underwater vehicle (HUV) was regarded as a research object by Liu et al. [28], and NSGA-II was adopted to establish a global approximation model, thus assisting the eventual optimal design of HUV. Nevertheless, these methods ignored the optimization designs of the global architectures of dissimilar control systems. As mentioned above, the optimization design of MECS is to solve the multi-objective optimization issue, where some contradictory indicators require all-round considerations. Therefore, the NSGA-II, combined with the AHP [29–32] method, was adopted in this paper to solve this problem, in which the binary encoding mode was employed and the decision variables were confined in discrete space. Moreover, the mutation process was modified to better match the characteristics of the problem. Eventually, the case study verified the validity of our approach.

The other sections of the paper are arranged as follows: Section 2 demonstrates the system descriptions. Section 3 presents the optimization approach on the foundation of NSGA-II and AHP method. Section 4 displays the case study, while the conclusions are presented in Section 5.

## **2. Mathematical Modeling of More Electric Control System (MECS)**

According to the design specification of safety critical systems, more electric control systems require high reliability and safety. For example, the possibility of mission failure per flight owing to certain material damages in the flight control system should not surpass the upper limit. Normally, the failure rate of a flight control system is 10−9/*h* to 10−10/*h* for a commercial aircraft. Hence, it is necessary to utilize a dissimilar flight control system to maintain high reliability.

## *2.1. Single Control System Structure*

A schematic diagram of the basic structure of a more electric control system, in which the controller, actuator, and sensor are main components, is illustrated in Figure 1. Since the electric power supply and hydraulic power supply are provided simultaneously, the actuator can be selected from HAs, EHAs, or EMAs, as shown in Figure 1. A HA is powered by a hydraulic supply system while an EHA or EMA is powered by an electric supply system. Since a centralized hydraulic power supply system has the characteristics of highpower density and fast response, HAs are widely used in control systems. However, the control system is heavy and easily exposed to oil contamination when HAs are adopted as they need pipeline transmission from the centralized power supply to the actuator. An EHA consists of the motor, pump, and cylinder that replaces the pipeline transmission by wire transmission. Hence, the weight of an EHA is light, whereas its heat dissipation is poor. An EMA consists of the motor and the ball screws. The command is transmitted from the wire controlling the motor and ball screws to drive the load. Therefore, an EMA is simple and light weight, but it gets stuck easy. Hence, to determine which type of actuator should be used, one has to comprehensively consider the weight, performance, and reliability.

**Figure 1.** Structure of typical actuation system.

## *2.2. Redundant Configuration of Actuation System*

The more electric control of safety critical systems not only requires the high-precision performance but also requires extremely high reliability and safety. Thus, redundant designs are often adopted, such as dual-redundant actuators and triple-redundant actuators. In order to integrate the advantages of hydraulic power and electric power, the typical isomorphic actuation system and heterogeneous actuation system shown in Figures 2–4 are used in real-world applications. In an aircraft, the HA/HA, HA/EHA, and HA/EMA are classic dual-redundant actuation systems, as shown in Figures 2–4.

**Figure 2.** Configuration of actuation system with two HAs.

**Figure 3.** Configuration of actuation system with one HA and one EMA.

**Figure 4.** Configuration of actuation system with one HA and one EHA.

Generally, the regulation of a dissimilar redundant actuator is to design the different actuator powered by a different power supply for one surface. In such design, when any actuator or power supply system fails, the other can fulfill the task through fault switching. Table 1 presents the typical configuration of a dissimilar redundant actuation system.


**Table 1.** Typical configuration of dissimilar redundant actuation system.

**Remark 1.** *Triple redundancy was considered in this research because when the number of redundancies increases, the performance indicators of an MECS will increase accordingly, which is not beneficial for system operation. Besides, an actuation system with more than three redundancies is rarely used in practice*.

## *2.3. Redundant Configuration of More Electric Control System*

As mentioned above, dissimilar redundancy technology is widely used in safety critical system design in order to improve system reliability. In a commercial aircraft, various controllers, actuators, and power supplies are applied in the aircraft actuation system. A typical redundant actuation system configuration based on high reliability for a commercial aircraft is illustrated in Figure 5. Although multiple control computers, actuators, and power supply systems can achieve high reliability, weights and costs will increase correspondingly, with maintainability decreased as well. Therefore, it is imperative to optimize the quantity of redundant control computers, power supply systems, and actuators. At same time, designers have to solve the common faults of dissimilar power supplies and dissimilar actuators as shown in Figures 2–4.

**Figure 5.** Configuration of commercial aircraft redundant actuation system.

In Figure 5, *I* expresses the information nodes of the control computer; *P* describes the power supply nodes; *A* is the actuator nodes; and *S* shows the surface of the aircraft. Herein, the actuator includes an HA, EHA, and EMA, as shown in Figures 2–4. Both the fonts of actuation nodes and rudder surfaces are bold since they are our main considerations.

**Definition 1.** *Power supply module set is* **P** = {*Pi*}(*i* = 1, 2, . . . , *n*). *For aircraft shown in Figure 5, P*<sup>1</sup> = *H means the power supply; P*<sup>1</sup> *is the hydraulic power; and P*<sup>2</sup> = *E means that the power supply P*<sup>2</sup> *is electric power*.

**Definition 2.** *Actuator module set is* **A** = *Aj* (*j* = 1, 2, . . . , *m*). *For aircraft shown in Figure 5, A*<sup>1</sup> = *HA means that the actuator A*<sup>1</sup> *is a hydraulic actuator. A*<sup>2</sup> = *EHA means that the actuator A*<sup>2</sup> *is an electro-hydrostatic actuator, and A*<sup>3</sup> = *EMA means that the actuatorA*<sup>3</sup> *is a mechatronic actuator*.

**Remark 2.** *Here, the commercial aircraft redundant actuation system is introduced as a typical MESC merely. In fact, our method is extensible and can be applied to more kinds of systems. In Section 4, we abstract a typical MECS with 3 hydraulic power supplies and 2 electric power supplies as the research object, which can be also extended to more combinations.*

## **3. Multi-Objective Optimization of MECS Based on NSGA-II and AHP**

Though an MECS can improve system reliability and ensure safety, when the system structure is more sophisticated, the corresponding structural indicators such as weight and cost will increase enormously. The design of an MECS may also increase the system volume and faults derived from the heavy weight. Therefore, we need to optimize factors such as the weight, power dissipation, and reliability of an MECS to ensure a compromise between these indicators.

## *3.1. Multi-Objective Optimization Modeling of MECS*

The primary task for optimization design is to transfer the optimization problem into mathematical descriptions. Generally, an optimization problem is composed of three elements: objective functions, decision variables, and constraint conditions. The optimization problem of an MECS can be expressed as

$$\begin{array}{c} \mathcal{I} = \max \mathcal{R}\_s(t) \\ \mathcal{R}\_s(t) = f(H\_i, \mathcal{E}\_i, HA\_j, EHA\_j, EMA\_j) \\ \text{s.t.} 0 \le m(H\_i, \mathcal{E}\_i, HA\_j, EHA\_j, EMA\_j) \le \mathcal{M} \\ 0 \le \psi(H\_i, \mathcal{E}\_i, HA\_j, EHA\_j, EMA\_j) \le \mathcal{V} \\ 0 \le c(H\_i, \mathcal{E}\_i, HA\_j, EHA\_j, EMA\_j) \le \mathcal{C} \\ \vdots \\ 0 \le g(H\_i, \mathcal{E}\_i, HA\_j, EHA\_j, EMA\_j) \le \mathcal{G} \\ i \in = \mathcal{L} = \{1, \dots, m\} \subset \mathbb{N} \\ j \in = \mathcal{F} = \{1, \dots, n\} \subset \mathbb{N} \end{array} \tag{1}$$

where function *J* is the maximum reliability; *m*(*H<sup>i</sup>* , *E<sup>i</sup>* , *HA<sup>j</sup>* , *EHA<sup>j</sup>* , *EMAj*) is the actual evaluation mass of the MECS; *M* is the superior limit of mass; *ψ*(*H<sup>i</sup>* , *E<sup>i</sup>* , *HA<sup>j</sup>* , *EHA<sup>j</sup>* , *EMAj*) expresses the actual evaluation power dissipation of the MECS; Ψ is the superior limit of power dissipation; *c*(*H<sup>i</sup>* , *E<sup>i</sup>* , *HA<sup>j</sup>* , *EHA<sup>j</sup>* , *EMAj*) describes the actual evaluation cost of the MECS; *c* represents the superior limit of cost; *g*(*H<sup>i</sup>* , *E<sup>i</sup>* , *HA<sup>j</sup>* , *EHA<sup>j</sup>* , *EMAj*) represents the other actual evaluation indicators of MECS; and *G* represents the superior limit of other indicators. The specific evaluation methods of objective functions and constraints are given below.

## 3.1.1. Objective Function

According to Equation (1), reliability is the most essential indicator and the objective function. System reliability refers to the ability of the system to meet the specified functions within the specified time and under the specified conditions. Only when specific functions of an MECS are controlled effectively at the same time can the overall function of the system be guaranteed. Therefore, the overall functional reliability of an MECS is defined as the ability to effectively control the specific functions simultaneously within the specified time and under the specified conditions.

Based on the definition above and by considering an MECS with eight actuation functions, the functional reliability of an MECS can be expressed as

$$R\_S = \Pr\{F\_1 \cap F\_2 \cap \dots \cap F\_8\} = \Pr\{S\} \tag{2}$$

where *F*1, . . . , *F*<sup>8</sup> represent the pivotal functions that MECS should accomplish.

By considering the multiple control surface combinations involved in the realization of system functions, Equation (2) can be rewritten as

$$R\_S = \Pr\left\{ \left( \stackrel{p}{\underset{x=1}{\cup}} F\_1 \right) \cap \left( \stackrel{q}{\underset{y=1}{\cup}} F\_2 \right) \dots \cap \left( \stackrel{m}{\underset{z=1}{\cup}} F\_8 \right) \right\} \tag{3}$$

where *<sup>p</sup>* ∪ *x*=1 *F*1 ∩ *q* ∪ *y*=1 *F*2 · · · ∩ *m* ∪ *z*=1 *F*8 represents a certain combination that realizes the pivotal functions.

To simplify the calculation, the program first calculated the overall reliability of electric power supplies *E*1, *E*<sup>2</sup> and then calculated the overall reliability of hydraulic power supplies *H*1, *H*2. Next, we calculated the minimum path and the program performing disjoint operation. Through decoupling analysis, nine minimal paths were obtained. Consequently, the functional reliability of an MECS can be calculated by

$$R\_S = \Pr\left\{ \bigcup\_{i=1}^{9} \mathcal{S}\_i \right\} = q\_{E1} q\_{E2} q\_{E5} q\_{E6} p\_{E9} p\_{ED} p\_{E1} p\_{E1} p\_{EN} p\_{E0} p\_{EP} p\_{E9} p\_{E8} p\_{ES} p\_{ES} p\_{E1}$$

*pEV pEW pEX pH*<sup>1</sup> *pH*<sup>3</sup> *pP*<sup>1</sup> *pP*<sup>2</sup> *pK*<sup>1</sup> *pK*<sup>2</sup> *pDD pYY pC*<sup>1</sup> *pC*<sup>2</sup> *pC*<sup>3</sup> + *qE*<sup>1</sup> *pE*2*qE*5*qE*<sup>6</sup> *pEB pED pEH pEJ pEN pEO pEP pEQ pER pES pET pEU pEV pEW pEX pH*<sup>3</sup> *pP*<sup>1</sup> *pP*<sup>2</sup> *pK*<sup>1</sup> *pK*<sup>2</sup> *pDD pYY pC*<sup>1</sup> *pC*<sup>2</sup> *pC*<sup>3</sup> + · · · + *qE*<sup>1</sup> *pE*<sup>2</sup> *pE*<sup>5</sup> *pEB pED pEH pEJ pEN pEO pEP pEQ pER pES pET pEU pEV pEW pEX pH*<sup>1</sup> *pH*<sup>3</sup> *pP*<sup>1</sup> *pP*<sup>2</sup> *pK*<sup>1</sup> *pK*<sup>2</sup> *pDD pC*<sup>1</sup> *pC*<sup>2</sup> *pC*<sup>3</sup> + *pE*<sup>1</sup> *pE*<sup>5</sup> *pEB pED pEH pEJ pEN pEO pEP pEQ pER pES pET pEU pEV pEW pEX pH*<sup>1</sup> *pH*<sup>3</sup> *pP*<sup>1</sup> *pP*<sup>2</sup> *pK*<sup>1</sup> *pK*<sup>2</sup> *pDD pC*<sup>1</sup> *pC*<sup>2</sup> *pC*<sup>3</sup> (4)

**Remark 3.** *Equation (4) gives the relationship between system reliability and component system. In fact, system reliability is determined by each component reliability, that is, the actuator reliability, while actuator reliability is associated with weight, power consumption, and cost. These indicators will eventually have an influence on system reliability.*

## 3.1.2. Constraint Conditions

The aim of optimization design is to identify the optimal solution from the practical solutions. The optimal solution can meet the goal of design as far as possible. The elevation of weight will restrict the MECS installation power; hence, weight is a vital index for system property. Power dissipation is another vital property index. In the optimization design of an MECS, particularly for a highly-efficient MECS, how to decrease the power dissipation with the existent technology and apparatus is a challenge which need to be tackled by designers. Meanwhile, as an MECS has to harbor high performance and reliability, so the reliability of an MECS should be considered. For practical application, the cost is an indispensable indicator that must be evaluated. Eventually, weight, power dissipation, cost, and reliability are determined as constraint conditions.

## • Weight

Weight is a vital evaluation index and decisive factor of an MECS. When we evaluate weight, the difficulties usually include various components, many of which exhibit a nonlinear growth relationship with design requirements, and serious coupling with other systems with different types of actuators. For an EHA, since it consists of the integration block, cylinder, pump, and motor, its weight evaluation can be presented as

$$M\_{EHA} = M\_{block} + M\_{cylinder} + M\_{pump} + M\_{motor} \tag{5}$$

where *Mblock*, *Mcylinder*, *Mpump* refer to the weight of the integration block, cylinder, and pump, respectively. *Mmotor* denotes the weight of the brushless direct-current motor (BLDC). All of them have the same unit, kilograms.

A permanent magnet brushless direct-current motor (BLDC) is utilized in the EHA due to its satisfactory control property and great power-to-mass ratio. The weight forecast of the BLDC is expressed by

$$M\_{\text{motor}} = 0.628T^{3/3.5} + 0.783\tag{6}$$

in which *T* represents the torque of BLDC.

The pump weight is proportionate to pump output. Thus, *Mpump* can be expressed as [33]

$$M\_{pump} = 0.339D + 2.038\tag{7}$$

The integration block is the EHA frame and involves the indispensable parts such as the nonreturn valve, filter, and accumulator. The overall weight is speculated by the EHA power, which is expressed as [33]

$$M\_{block} = 0.105P\_{EHA} + 2\tag{8}$$

where *PEHA* denotes the maximal power of EHA.

The fluid cylinder weight includes the four parts stated below:

$$M\_{\text{cylinder}} = M\_{\text{cover}} + M\_{\text{shell}} + M\_{\text{piston}} + M\_{\text{rod}} \tag{9}$$

where *Mcover*, *Mshell*, *Mpiston*, *Mrod* are the weights of the cylinder cover, shell, piston, and plunger rod, respectively. All of them can be calculated by

$$\begin{cases} M\_{\text{rod}} = \frac{\pi}{4} \times d\_{\text{rod}}^2 \times L\_{\text{rod}} \times \rho\_{\text{steel}}\\ M\_{\text{piston}} = A \times t\_{\text{piston}} \times \rho\_{\text{copper}}\\ M\_{\text{shell}} = \frac{\pi}{4} \times \left(d\_{\text{shell}}^2 - \frac{4A}{\pi}\right) \times L\_{\text{shell}} \times \rho\_{\text{steel}}\\ M\_{\text{cover}} = 2 \times \frac{\pi}{4} \times d\_{\text{shell}}^2 \times t\_{\text{cover}} \times \rho\_{\text{steel}} \end{cases} \tag{10}$$

Thus, the total weight can be evaluated as

$$M = \sum\_{i=1}^{n\_{HA}} M\_{HA} + \sum\_{i=1}^{n\_{EHA}} M\_{EHA} + \sum\_{i=1}^{n\_{EMA}} M\_{EMA} + M\_{elec} + M\_{pipe} + M\_{wire} \tag{11}$$

where *MHA*, *MEHA*, *MEMA* represent the weights of HA, EHA, and EMA, respectively, and *Melec*, *Mpipe*, *Mwire* represent the weights of electric equipment, pipe and wire, respectively.

• Power efficiency

For an EHA, the motor converts electric energy to hydraulic power, pushing the cylinder to realize the movement of the control surface. For an EHA system, when its control surface load is certain, it saves more energy, resulting in lower power dissipation. The control surface torque *TS*, control surface velocity *ωS*, and pressure *P<sup>S</sup>* are known

design specifications. The aim of decreasing the power dissipation of an EHA is to decrease the power output of the motor *P*, which is written as

$$P = T \times \mathcal{W} \tag{12}$$

The motor is directly connected to the pump, so the torque and speed of the motor and pump are the same. The torque of the pump is expressed as

$$T = f\_{pm} \times \frac{d\omega}{dt} + \mathcal{K}\_{fric} \times \omega + (p\_f + 2p\_{pipe}) \times \frac{D}{2\pi} \tag{13}$$

where *Jpm* is the pump rotational inertia; *Kf ric* denotes the viscosity factor, *Jpm* × *dω*/*dt* = 0; the pump rotation velocity is constant; *p<sup>f</sup>* denotes the differential pressure between the 2 cylinder cavities; and *ppipe* denotes the pressure consumption within the pipe. Hence, Equation (13) is written as

$$T = K\_{fric} \times \omega + (p\_f + 2p\_{pipe}) \times \frac{D}{2\pi} \tag{14}$$

in which *p<sup>f</sup>* = *F*/*A*, *F* is the loading force of an EHA and *A* denotes the piston effective area, which can be calculated by

$$A = \frac{kF}{P\_S} \tag{15}$$

in which *k* > 1 denotes a logical excess margin, and *F* can be described by

$$F = \frac{T\_S \times \sin(\theta + 30^\circ)}{L} \tag{16}$$

For an EMA, the power dissipation is determined by the ball screw. Therefore, Equation (17) gives the calculation method for the drive torque of the ball screw:

$$T\_a = \frac{F\_a \times I}{2\pi \times n\_1} \tag{17}$$

where *T<sup>a</sup>* represents the drive torque; *I* is the lead screw; *n*<sup>1</sup> is the positive efficiency of the feed screw; and *F<sup>a</sup>* refers to the axial load, which can be expressed by

$$F\_a = F + \mu mg\tag{18}$$

where *F* means the axial cutting force of the lead screw, *µ* is the comprehensive friction coefficient, and *m* refers the weight of the worktable and workpiece.

After the calculation of the drive torque *Ta*, the motor power can be determined accordingly; thus, the corresponding power dissipation of the motor can be obtained.

For the power dissipation evaluation of a HA, the volumetric efficiency of the pump *η*<sup>0</sup> is the main affecting factor, which can be described by

$$
\eta\_0 = \frac{Q}{Q\_0} \times 100\% \tag{19}
$$

where *Q*<sup>0</sup> means the theoretical flow of the pump, and *Q* is the actual pump flow, which is written as

$$Q = Q\_{pump} + Q\_{pipe} + Q\_{cylinder} \tag{20}$$

where *Qpump* denotes the pump leak, *Qpipe* denotes the loss of flow within the hydraulic tube, and *Qcylinder* denotes the cylinder flow.

$$\begin{cases} Q\_{pump} = \mathfrak{f} \times (p\_1 - p\_2) = \mathfrak{f} \times \mathcal{F} / A \\ Q\_{cylinder} = A \times v \\ Q\_{pipe} = 2 \times \mathfrak{f} \times p\_{pipe} \end{cases} \tag{21}$$

where *ξ* is the pump leakance coefficient, *ppipe* denotes the pipe pressure drop, and *v* denotes the cylinder speed.

• Cost

For an MECS, cost is another index worthy of consideration. In this paper, the total costs are divided into two parts:

$$\mathbf{C} = \mathbf{C}\_{manu} + \mathbf{C}\_{oper} \tag{22}$$

where *C* represents the total cost of an MECS, *Cmanu* refers to the manufacturing cost, and *Coper* represents the operation cost.

However, when we evaluated the total cost, the component cost exhibited a nonlinear growth relationship with the requirement. To simplify this process, this paper utilized different configurations to obtain the costs of all components. In the evaluation of manufacturing costs, the main components considered included the HA, EHA, EMA, oil boxes, pipelines, cables, engine-driven pump (EDP), and electric motor driven pump (EMP). The evaluation was conducted mainly based on the similarity principle. Thus, the manufacturing costs is expressed as

$$\mathbf{C}\_{\text{manu}} = \sum\_{i=1}^{n\_{\text{edp}}} \mathbf{C}\_{\text{edp}} + \sum\_{i=1}^{n\_{\text{emp}}} \mathbf{C}\_{\text{emp}} + \sum\_{i=1}^{n\_{\text{tank}}} \mathbf{C}\_{\text{tank}} + \sum\_{i=1}^{n\_{\text{motor}}} \mathbf{C}\_{\text{water}} + \sum\_{i=1}^{n\_{\text{att}}} \mathbf{C}\_{\text{atf}} + \sum\_{i=1}^{n\_{\text{pipe}}} \mathbf{C}\_{\text{pipe}} + \sum\_{i=1}^{n\_{\text{wire}}} \mathbf{C}\_{\text{wire}} \tag{23}$$

where *Cedp* represents the cost of the EDP; *Cemp* refers to the cost of the EMP; *Ctank* is the cost of the tank; *Cmotor* is the cost of the motor; *Cact* describes the cost of the actuator; *Cpipe* expresses the cost of the pipe; and *Cwire* denotes the wire cost.

The operation expense primarily denotes the Direct Operating Costs (DOCs), which are changeable costs directly derived from operating the aircraft [34]. The fuel expense and maintenance expenditure of DOCs are changeable. Hence, the DOCs are predominantly related to fuel expense and maintenance expenditure herein. The operation expense speculation is written as

$$\mathbb{C}\_{oper} = \mathbb{C}\_{oil} + \mathbb{C}\_{main} \tag{24}$$

where *Coil* and *Cmain* represent the fuel cost and maintenance cost, respectively.

Based on the statistical proportion of main flight control maintenance cost in the total maintenance expenditure, the fuel expense speculation primarily considers the fuel consumption during the life of the commercial aircraft. The function of fuel cost with respect to time is expressed as

$$\mathbb{C}\_{oil}(t) = \mathbf{50} \times \mathbb{C}\_0 \times \text{Weight} \times (1.02^t - 1) \tag{25}$$

where *C*<sup>0</sup> means the average annual fuel cost of the aircraft and *Weight* represents the fuel weight.

The maintenance cost of the aircraft in the whole life cycle mainly has two parts [35]. One is called the cyclic cost *CC*, which is related to the take-off and landing of the aircraft, such as the maintenance of its braking device, flap system, and landing gear. The other cost is related to the flight time of the aircraft, such as the regular loss and replacement of parts caused by each hour of flight. This part of the cost is called hourly cost *CH*. The evaluation of the maintenance cost is usually accomplished through the statistic of man-hour cost. Equation (26) provides the method for man-hour cost calculation

$$\mathbb{C}\_{\text{main}} = \mathbb{C}\_{H} \times t + \mathbb{C}\_{\mathbb{C}} \times t + M \times (1/\mathfrak{a}) \times \mathbb{R} \times t \tag{26}$$

## 3.1.3. Design Variables

The expression of a goal suggests that there are substantial that which ought to be identified in configuration designs. Nevertheless, considering all those variables during optimization makes convergence difficult and is quite time-consuming. Hence, merely vital variables having a remarkable influence on MECS property ought to be optimized. Overall, the selective principles of the design variates are stated below:


According to the aforesaid principles, the optimization variates are chosen via:

$$\mathbf{x} = \left\{ \begin{array}{c} \mathbf{x}\_1 \\ \mathbf{x}\_2 \\ \mathbf{x}\_3 \end{array} \right\} = \left\{ \begin{array}{c} n\_{HA} \\ n\_{EHA} \\ n\_{EMA} \end{array} \right\} \tag{27}$$

## *3.2. Multi-Objective Optimization Based on NSGA-II*

The optimization issue herein aimed at the configuration design of an MECS. NSGA-II is proper for optimizing multi-objective issues as it exhibits strong distributed capability and rapid convergence. The flow chart of NSGA-II is shown in Figure 6. Subsequently, the encoding/decoding mode and fast non-dominated sorting process are elaborated.

**Figure 6.** Flow chart of NSGA-II.

## 3.2.1. Encoding and Decoding

In our research, the binary encoding mode was adopted in NSGA-II, which was simple, efficient, and easy to design and use for implementing crossover and mutation operations. Two types of energy information were encoded through this mode, while the configuration of each actuation system was also expressed in binary. Figure 7 illustrates this coding mode.

**Figure 7.** Binary encoding mode adopted in NSGA-I.

## 3.2.2. Fast Non-Dominated Sorting

Compared with non-dominated sorting approach utilized in NSGA arithmetic, the rapid non-dominated sorting utilized in NSGA-II algorithm requires lower computational complexity, which is only about *O*(*MN*<sup>2</sup> ). The overall algorithm flow was presented as follow:

**Algorithm 1:** Fast non-dominated sorting

1: Fast-non-dominated-sort (*P*) 2: *for each p* ∈ *P* 3: *for each q* ∈ *P* 4: **if** *p* ≤*<sup>n</sup> q*, **then** # if *p* is dominated by *q*, then add *q* to *S<sup>p</sup>* 5: *S<sup>p</sup>* = *S<sup>p</sup>* ∪ *q* 6: **else if** *q* ≤*<sup>n</sup> p*, **then** 7: *n<sup>p</sup>* = *n<sup>p</sup>* + 1 8: **if** *n<sup>p</sup>* = 0, **then** 9: *prank* = 1, *F*<sup>1</sup> = *F*<sup>1</sup> ∪ *P* **#** when *n<sup>p</sup>* of the individual is 0, then this individual is the first level of Pareto 10: The comparison of dominating relationships between individuals, *Sp* and *np*, are introduced for storage and records, respectively; ≤*<sup>n</sup>* represents the comparison of dominating relationships. The solution of *n<sup>p</sup>* = 0 is stored in the records of level 1, and the solution of level 1 has higher priority than that of level 2. 11: *i* = 1. 12: **while** *F<sup>i</sup>* = ∅ **do** 13: H=∅ 14: *for each p* ∈ *F<sup>i</sup>* 15: *for each q* ∈ *S<sup>p</sup> #* Sort all the individuals in *S<sup>p</sup>* 16: *n<sup>q</sup>* = *n<sup>q</sup>* − 1 17: **if** *n<sup>q</sup>* = 0, **then # when** *n<sup>q</sup>* of the individual is 0, it is a non-dominated individual 18: *qrank* = *i* + 2, *H* = *H* ∪ *q* # The Pareto level of this individual is the current highest level plus 1. At this moment, the initial value of i was 0, so we added 2. 19: **end while** 20: *i* = *i* + 1, *F<sup>i</sup>* = *H* 21: Loop the program to obtain level 2, level 3 . . . The computational complexity is *O*(*MN*<sup>2</sup> )

## 3.2.3. Crowding Degree Calculation

Crowing degree *I<sup>d</sup>* refers to the local crowding distance between any two adjacent points whose level are the same in the target space. The purpose of introducing crowding degree is to improve the distribution uniformity of the Pareto solution set. Furthermore, it could not only boost the diversity of population but also enhance system robustness. The crowding degree *I<sup>d</sup>* could be expressed as the length of the maximum rectangle of individual *i* on both sides, where the rectangle only includes *i* itself.

The main procedure to determine the degree of individual congestion involves three steps:


• Define the crowding degree of marginal individuals as a larger number *L*(1)*<sup>d</sup>* = *L*(*l*)*<sup>d</sup>* = *M* to prioritize individuals on the sorting edge; thus, the crowding degree of any other individual *I<sup>d</sup>* can be expressed as

$$I\_d = \sum\_{j=1}^{m} \left( \left| f\_j^{i+1} - f\_j^{i-1} \right| \right) \tag{28}$$

where *j* refers to each evaluation indicator; *f i*+1 *j* represents the *jth* evaluation indicator value of individual *i* + 1; and *f i*+1 *j* represents the *jth* evaluation indicator value of individual *i* − 1.

Through the aforesaid steps and calculation, every individual was endowed with two attributes, i.e., the crowding distance and the rank. This laid the foundation for the follow-up processes.

## 3.2.4. Optimal Selection

The optimal selection avoids the loss of effective factors and ensures the survival rate of high-performance individuals, so it can ensure that the Pareto optimal solution is continuously optimized. The optimal selection can not only improve the efficiency and convergence of the optimization but also ensure the uniformity of the optimization process. The selection process is completed by comparing the results of fast non-dominated sorting and crowding calculation, and better individuals are selected after comparison. Two steps are involved during this process.


## 3.2.5. Crossover

The main function of crossover is to simulate gene recombination in the process of heredity and evolution. There are various approaches to implement it, such as uniform crossover, multi-point crossover, and binary crossover. In our research, the simulated binary crossover is adopted for crossover operation, which is expressed as

$$\begin{array}{l} y\_{1j} = 0.5 \times \left[ \left( 1 + \gamma\_{j} \right) \mathbf{x}\_{1j} + \left( 1 - \gamma\_{j} \right) \mathbf{x}\_{2j} \right] \\ y\_{2j} = 0.5 \times \left[ \left( 1 - \gamma\_{j} \right) \mathbf{x}\_{1j} + \left( 1 + \gamma\_{j} \right) \mathbf{x}\_{2j} \right] \end{array} \\ \text{(29)}$$
 
$$\begin{cases} \left( 2u\_{j} \right)^{\frac{1}{\eta + 1}} \text{if } u\_{j} \le 0.5 \\ \left( \frac{1}{2\left( 1 - u\_{j} \right)} \right)^{\frac{1}{\eta + 1}} \text{else} \end{cases} \tag{29}$$

where *u<sup>j</sup>* ∈ (0, 1); *x*1*<sup>j</sup>* and *x*2*<sup>j</sup>* are parent individuals; *y*1*<sup>j</sup>* and *y*2*<sup>j</sup>* are offspring individuals; and *η* > 0 refers to the cross-distribution index. Generally, *η* = 20 is the best value of the cross-distribution index by default.

## 3.2.6. Mutation

Mutation is widely applied to simulate variation links in the process of biological heredity and evolution, which refers to the replacement between genes, that is, substitute the other gene values with the alleles on the locus to create new individuals. Through mutation, not only could the local searchability be improved but also the population diversity is guaranteed.

This paper introduced the polynomial mutation method into the research on multiobjective optimization problems. The form of mutation operator is expressed as

$$V\_k' = V\_k + \delta(u\_k - L\_k)$$

$$\delta = \begin{cases} \left[ 2u + (1 - 2u)(1 \mid -\delta\_1)^{\eta\_m + 1} \right]^{\frac{1}{\eta\_m + 1}} - 1 \text{ if } u \le 0.5\\ 1 - \left[ 2(1 - u) + 2(u - 0.5)(1 - \delta\_2)^{\eta\_m + 1} \right]^{\frac{1}{\eta\_m + 1}} \text{else if } u > 0.5 \end{cases} \tag{30}$$

where *δ*<sup>1</sup> = (*Vk*−*L<sup>k</sup>* ) (*uk*−*L<sup>k</sup>* ) , *δ*<sup>2</sup> = (*uk*−*V<sup>k</sup>* ) (*uk*−*L<sup>k</sup>* ) , *u* denotes a stochastic number in interval [0,1], *η<sup>m</sup>* denotes the distribution index, and *V<sup>k</sup>* is a parent individual.

## *3.3. Comprehensive Evaluation of System Configuration Based on AHP*

The analytic hierarchy process (AHP) is a multi-objective decision analysis approach that combines qualitation and quantitation analyses, where elements associated with task decisions are divided into the object level, criterion level and scheme level. The AHP mathematizes the decision-making via a few quantitation data based on deep analyses, such as the influential factors and inner association of the decision-making problems. AHP can offer an easy decision approach for intricate decision issues with several standards and no evident structure features. Questions to be evaluated in MECS optimization design ought to be methodized and layered, and the structure model of hierarchy analysis ought to be constructed. A three-level hierarchy structure of AHP was presented by Figure 8. Level 1 is the general objective of decision-making, which denotes the optimum design of MECS. Level 2 denotes the criterion layer, where the weight, power dissipation, expense, and dependability are utilized as evaluation standards of system multi-objective optimization design. Level 3 denotes the scheme layer, where the solutions are determined in the Pareto frontier acquired by NSGA-II. Besides, the process of AHP mainly includes the judgment matrix structure, weight coefficient calculation, consistency of judgment matrix verification, and weight coefficients of goal level calculation. The flow chart of AHP was shown in Figure 9.

**Figure 8.** Three-level hierarchical framework of AHP.

**Figure 9.** Flowchart of AHP.


$$J = \begin{bmatrix} 1 & 3 & 9 & 1 \\ \frac{1}{3} & 1 & \frac{1}{3} & \frac{1}{3} \\ \frac{1}{9} & \frac{1}{3} & 1 & \frac{1}{9} \\ 1 & 3 & 9 & 1 \end{bmatrix} \tag{31}$$

3. Validate the judgement matrix coherence. Coherence index is computed via:

$$\begin{cases} \ C I = \frac{\lambda\_{\text{max}} - n}{n - 1} \\ \ C R = \frac{C I}{\mathbb{Z} \Gamma} \end{cases} \tag{32}$$

where *CI* denotes the coherence index of the judgment matrix, *RI* denotes the average stochastic coherence index of the matrix (specific values were presented by Table 2), *CR* denotes the stochastic coherence ratio of the matrix, *λ*max denotes the maximal value of the matrix characteristic value, and *n* denotes the matrix order.

**Table 2.** Values of the average stochastic coherence index.


By computation, if *λ*max = 3.7638 and *CI* = (3.7638 − 4)/(4 − 1) = −0.0787, then *CR* = −0.0787/0.90 = −0.0875 < 0.1, which suggests that the weight matrix is coherent. When the consistence test is not met, return to step 2 and reconstruct the judgement matrix.

4. Compute the weighted coefficient between the contrasted elements with the relevant standards. Compute the continued product *M<sup>i</sup>* of each row element in *A*, the product of every row element, and its *n*-th root *w*1.

$$\begin{cases} \begin{array}{c} M\_i = \prod\_{i=1}^n \mathbf{x}\_{ij} = \begin{bmatrix} .27 & 0.0370 & 0.0041 & .27 \end{bmatrix}^T, i = 1, \cdots, n \\\ \varpi\_1 = \sqrt[n]{M\_i} = \begin{bmatrix} .2.2795 & 0.4387 & 0.2533 & .2.2795 \end{bmatrix}^T, i = 1, \cdots, n \end{cases} \end{cases} \tag{33}$$

Normalize *w*<sup>1</sup> as *W<sup>i</sup>* = *w*1/∑ *n <sup>i</sup>*=<sup>1</sup> *w*1, where *W<sup>i</sup>* denotes the weighted coefficient of every factor.

$$\mathcal{W} = \begin{bmatrix} \ 0.4341 & 0.0835 & 0.0482 & 0.4341 \end{bmatrix} \tag{34}$$

5. Speculate the design in the Pareto frontier as per the weighted coefficients of every standard, and afterwards get the optimum design of MECS.

## **4. Case Study and Discussion**

The diagram of an MECS with three fluid power supplies and two electric power supplies (3H2E) was presented by Figure 10. The tendency toward more electricity causes the transformation from hydraulic power into standby energy, or the replacement of hydraulic power with electrical energy, which produces more options for all control surfaces. For every actuator, it can be an HA linked to hydraulic source or an EHA linked to electrical source. In addition, each actuator is controlled by at least one controller. Therefore, the number of alternative configurations of the MECS is exceedingly high.

**Figure 10.** A typical 3H2E more electric control system (MECS).

In this section, we introduced the proposed optimization method into the case study of the MECS illustrated in Figure 10. By setting the maximum optimization iterations of NSGA-II as 200, and population size as 500, we obtained the Pareto frontier of the optimization result illustrated in Figure 11. *W* refers to the weight of the MECS, *P* is the power consumption of the MECS, and *C* represents the expense of the MECS. The redround-dot-curved surface is a 3D Pareto optimum front curved surface. Every dot denotes an actuation design. The green line, marked by hollow pentagram; blue line, marked by cross symbol; and yellow line, marked by triangle were projected on three planes. The Pareto front was utilized to realize multi-objective optimization and to achieve a certain decrease in configuration quantity.

**Figure 11.** 3D-Pareto frontier of multi-objective optimization of MECS.

Moreover, the graph in Figure 11 presents the relationship among the three objectives. The line in x, y plane presents that weight and power have the relation of promotion. The line in x, z plane presents that weight and expense are contradictory. The line in y, z plane displays contradictory power and expense of an MECS. When the power decreases, the expense is elevated. The design objective is to minimize power and expense; hence, those two functions conflict with one another. Due to the conflicts between design indexes, it is hard to make each index optimum simultaneously.

In addition, to better illustrate the relationship between these three objectives and system reliability, we projected Figure 11 onto three 2-D planes, respectively, and obtained the 2D-pareto frontiers in Figures 12–14. Reliability was used as the scale parameter to color different system configurations, i.e., the brighter the color, the higher the reliability.

**Figure 12.** 2D Pareto frontier between weight, power consumption, and reliability.

**Figure 13.** 2D Pareto frontier among weight, cost, and reliability.

**Figure 14.** 2D Pareto frontier among power consumption, cost, and reliability.

Through multi-objective optimization, the Pareto optimization frontier was acquired. Nevertheless, finding the best way to identify the optimal solution in the Pareto frontier remains a daunting challenge. Thus, AHP was introduced to solve this problem. According to Section 3.3, diverse optimum solutions can be realized via structuring diverse judgment matrices when we use AHP to perform decision analyses. Thus, to better observe the association between the goals more, the judgment matrix was presented by Table 3, and the eventual result was illustrated in Figures 15 and 16.


**Table 3.** The association between the objectives.

**Figure 15.** AHP assessment outcomes of Pareto frontier.

**Figure 16.** Reliability comparison among different configurations.

Through a comprehensive analysis of AHP, the score of each solution on the Pareto frontier was acquired. The solution with the greatest score is the best one. As presented by Figures 15 and 16, the optimal design strategy is number 96 of the Pareto frontiers, with triple redundancy configuration used in the MCES. The reliability of our designed configuration increased by 5.89% and 55.34% respectively, compared with the dual redundancy and single redundancy configurations. The results indicated that our approach combining multi-objective optimization and decision-making could realize the multi-objective optimization design of an MECS. The solutions on the Pareto frontier acquired by NSGA-II arithmetic were uniformly distributed, which suggested that we could offer designers more diversified choices. Subsequently, AHP was employed for the eventual decision analyses, which enabled designers to introduce predilections and experiences for different goals.

## **5. Conclusions and Future Work**

The present research demonstrated the multi-objective optimization for an MECS. By designing the level length, the objectives are optimized. However, it is hard to minimize every goal simultaneously due to the conflicts between these objectives. Thus, an appropriate multi-objective optimization algorithm, which can compromise multiple objectives, is required to acquire the Pareto optimum solution set. Besides, since the Pareto frontier solution is not a solution but a solution set, an approach of multi-objective decision-making analysis is required to identify the optimum solution in the Pareto frontier. Herein, a combination of multi-objective optimization and multi-objective decision-making was utilized. The NSGA-II was employed to determine the optimum solution set, i.e., the Pareto frontier. Subsequently, the optimum design was acquired via an AHP. The outcome revealed the practicability of our optimization approach. In addition, the present research discovered that the multi-objective optimization and multi-objective decision-making approaches could be utilized in the optimization design of an MECS, and that this approach suited other components and systems as well.

Future work should encompass the improvement of real-time performance for an optimization algorithm, with more evaluation indicators considered.

**Author Contributions:** Conceptualization, Z.L. and R.C.; methodology, Z.L. and R.C.; software, D.L.; validation, Z.L. and D.L.; writing—original draft preparation, Z.L.; writing—review and editing, S.W.; supervision, S.W. and J.S.; project administration, J.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China (No. 51620105010) and the program of China Scholarship Council (No. 202106020106).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** This study does not report any data.

**Conflicts of Interest:** The authors declare no conflict of interest.

## **References**


## *Article* **Novel Approach to Fault-Tolerant Control of Inter-Turn Short Circuits in Permanent Magnet Synchronous Motors for UAV Propellers**

**Aleksander Suti \* , Gianpietro Di Rito and Roberto Galatolo**

Department of Civil and Industrial Engineering, University of Pisa, Largo Lucio Lazzarino 2, 56122 Pisa, Italy; gianpietro.di.rito@unipi.it (G.D.R.); roberto.galatolo@unipi.it (R.G.)

**\*** Correspondence: aleksander.suti@dici.unipi.it; Tel.: +39-05-0221-7211

**Abstract:** This paper deals with the development of a novel fault-tolerant control technique aiming at the diagnosis and accommodation of inter-turn short circuit faults in permanent magnet synchronous motors for lightweight UAV propulsion. The reference motor is driven by a four-leg converter, which can be reconfigured in case of a phase fault by enabling the control of the central point of the motor Y-connection. A crucial design point entails the development of fault detection and isolation (FDI) algorithms capable of minimizing the failure transients and avoiding the short circuit extension. The proposed fault-tolerant control is composed of two sections: the first one applies a novel FDI algorithm for short circuit faults based on the trajectory tracking of the motor current phasor in the Clarke plane; the second one implements the fault accommodation, by applying a reference frame transformation technique to the post-fault commands. The control effectiveness is assessed via nonlinear simulations by characterizing the FDI latency and the post-fault performances. The proposed technique demonstrates excellent potentialities: the FDI algorithm simultaneously detects and isolates the considered faults, even with very limited extensions, during both stationary and unsteady operating conditions. In addition, the proposed accommodation technique is very effective in minimizing the post-fault torque ripples.

**Keywords:** all-electric propulsion; electric machines; fault diagnosis; fault-tolerant control; inter-turn short circuit; modelling; simulation

## **1. Introduction**

The design of next-generation long-endurance UAVs is undoubtedly moving towards the full-electric propulsion. Though immature today, full-electric propulsion systems are expected to obtain large investments in the forthcoming years, to replace conventional internal combustion engines, as well as hybrid or hydrogen-based solutions. Full-electric propulsion would guarantee smaller CO2-emissions, lower noise, a higher efficiency, a reduced thermal signature (crucial for military applications), and simplified maintenance. However, several reliability and safety issues are still open, especially for long-endurance UAVs flying in unsegregated airspaces.

In this context, the Italian Government and the Tuscany Regional Government cofunded the project TERSA (*Tecnologie Elettriche e Radar per Sistemi aeromobili a pilotaggio remoto Autonomi*) [1], led by Sky Eye Systems (Italy) in collaboration with the University of Pisa and other Italian industries. The TERSA project aims to develop an Unmanned Aerial System (UAS) with a fixed-wing UAV, Figure 1, having the following main characteristics:


**Citation:** Suti, A.; Di Rito, G.; Galatolo, R. Novel Approach to Fault-Tolerant Control of Inter-Turn Short Circuits in Permanent Magnet Synchronous Motors for UAV Propellers. *Aerospace* **2022**, *9*, 401. https://doi.org/10.3390/ aerospace9080401

Academic Editor: Cengiz Camci

Received: 1 July 2022 Accepted: 25 July 2022 Published: 26 July 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

• Landing system: parachute and airbags; Propulsion system: Permanent Magnet Synchronous Motor (PMSM) powering a twin-blade fixed-pitch propeller;

 Take-off system: pneumatic launcher; Landing system: parachute and airbags;

• Propulsion system: Permanent Magnet Synchronous Motor (PMSM) powering a twin-blade fixed-pitch propeller; Innovative sensing systems:

*Aerospace* **2022**, *8*, x FOR PEER REVIEW 2 of 25

• Innovative sensing systems: o Synthetic aperture radar, to support surveillance missions in adverse environ-

Range: >3 km;


**Figure 1.** Rendering of the TERSA UAV. **Figure 1.** Rendering of the TERSA UAV.

With reference to the activities related to the TERSA UAV propulsion system development (which this work refers to), a special attention has been devoted to the architecture definition for obtaining fault-tolerant capabilities. In particular, propulsion systems based on three-phase PMSMs driven by conventional three-leg converters typically have failure rates from 100 to 200 per million flight hours [2], which are not compatible with the reliability and safety levels required for the airworthiness certification [3]. The failure rate of PMSMs is mostly driven by motor phase faults and converter faults (open-switch in a converter leg, open-phase, inter-turn, phase-to-leg, phase-to-ground, or capacitor shortcircuits [4]), which cover from 60% to 70% of the total fault modes [5]. Provided that the weight and envelopes required by UAV applications impede the extensive use of hardware redundancy (e.g., redundant motors), the reliability enhancement of full-electric pro-With reference to the activities related to the TERSA UAV propulsion system development (which this work refers to), a special attention has been devoted to the architecture definition for obtaining fault-tolerant capabilities. In particular, propulsion systems based on three-phase PMSMs driven by conventional three-leg converters typically have failure rates from 100 to 200 per million flight hours [2], which are not compatible with the reliability and safety levels required for the airworthiness certification [3]. The failure rate of PMSMs is mostly driven by motor phase faults and converter faults (open-switch in a converter leg, open-phase, inter-turn, phase-to-leg, phase-to-ground, or capacitor short-circuits [4]), which cover from 60% to 70% of the total fault modes [5]. Provided that the weight and envelopes required by UAV applications impede the extensive use of hardware redundancy (e.g., redundant motors), the reliability enhancement of full-electric propulsion systems can be achieved only through the redundancy of phases or stator modules [6] or by unconventional converters. By reducing design complexity, weight and envelopes, the use of four-leg converters [7–9], Figure 2, could represent a suitable solution for UAV applications.

pulsion systems can be achieved only through the redundancy of phases or stator modules [6] or by unconventional converters. By reducing design complexity, weight and envelopes, the use of four-leg converters, [7–9], Figure 2, could represent a suitable solution for UAV applications. In this converter topology, a couple of power switches are added, as stand-by devices, to the conventional three-leg bridge, enabling the control of the central point of the "star" connection (often Y-connection) of the motor. Nevertheless, to benefit from their faulttolerance capability, these devices must be promptly reconfigured when a fault occurs, and extremely fast Fault-Detection and Isolation (FDI) algorithms must be developed. Especially for PMSMs operating at high speeds (such as the ones used for aircraft propellers), phase faults can generate abnormal torque ripples and the related failure transients can potentially cause unsafe conditions [10,11].

**Figure 2.** Four-leg converter driving a three-phase PMSM with access to the central point. **Figure 2.** Four-leg converter driving a three-phase PMSM with access to the central point.

In this converter topology, a couple of power switches are added, as stand-by devices, to the conventional three-leg bridge, enabling the control of the central point of the "star" connection (often Y-connection) of the motor. Nevertheless, to benefit from their Many pieces of research on the development of FDI algorithms for PMSMs have been carried out in the last decades, and a special attention has been recently addressed to the Inter-Turn Short Circuit (ITSC) fault [12–15], which typically contributes to 10% of the PMSM failure rate [5,9,16].

fault-tolerance capability, these devices must be promptly reconfigured when a fault occurs, and extremely fast Fault-Detection and Isolation (FDI) algorithms must be developed. Especially for PMSMs operating at high speeds (such as the ones used for aircraft Faiz [12] carried out a comprehensive literature review on FDI methods for ITSC faults in PMSMs. Depending on the measurements processed by the algorithms, they can be categorized as torque-based [17,18], flux-based [19,20], electromagnetic parametersbased [21], voltage-based, and current-based.

propellers), phase faults can generate abnormal torque ripples and the related failure transients can potentially cause unsafe conditions [10,11]. Many pieces of research on the development of FDI algorithms for PMSMs have been carried out in the last decades, and a special attention has been recently addressed to the Inter-Turn Short Circuit (ITSC) fault [12–15], which typically contributes to 10% of the PMSM failure rate [5,9,16]. Faiz [12] carried out a comprehensive literature review on FDI methods for ITSC faults in PMSMs. Depending on the measurements processed by the algorithms, they can be categorized as torque-based [17,18], flux-based [19,20], electromagnetic parametersbased [21], voltage-based, and current-based. The most common strategy is given by the voltage-based methods. Immovilli and Bianchini [22] propose reconstructing harmonic patterns of the voltage components in the Clarke plane. The technique works very well in transient conditions, but it is not able to identify the shorted phase. The FDI algorithm proposed in [23], implements an estimator of back electromotive force and detects the ITSC by comparing the estimates to a reference The most common strategy is given by the voltage-based methods. Immovilli and Bianchini [22] propose reconstructing harmonic patterns of the voltage components in the Clarke plane. The technique works very well in transient conditions, but it is not able to identify the shorted phase. The FDI algorithm proposed in [23], implements an estimator of back electromotive force and detects the ITSC by comparing the estimates to a reference model. The method is very robust against load disturbances, but it requires a very accurate model of the Back Electromotive Force (BEMF). Other approaches use as fault symptoms the harmonic components of the Zero Sequence Voltage Component (ZSVC) signal. Urresty et al. [24] proposes a Vold–Kalman Filtering Order Tracking (VKF-OT) algorithm to track the ZSVCs first harmonic component. The method offers reliable indicators, which are applicable in transient phases too, but it needs to be extended in terms of fault isolation. A robust algorithm is proposed by Hang et al., in [25], where the fault indicators are defined to remove the influence of rotor speed variation and a frequency tracking algorithm based on reference frame transformation is used to extract the fault symptoms. A possible drawback is the dependency of the fault indicator on the PMSM parameters, which may vary with environmental temperature. Boileau [26] considered as a fault indicator the amplitude of the positive voltage sequence while Meinguet [27] used as an indicator the ratio between the amplitudes of negative and positive voltage sequences. These methods are not able to locate the shorted phase, neither.

model. The method is very robust against load disturbances, but it requires a very accurate Concerning the current-based methods, the most relevant one is based on the socalled Park Vector Approach (PVA) [28–31]. When the PMSM is operating in normal

model of the Back Electromotive Force (BEMF). Other approaches use as fault symptoms the harmonic components of the Zero Sequence Voltage Component (ZSVC) signal. Ur-

are applicable in transient phases too, but it needs to be extended in terms of fault isolation. A robust algorithm is proposed by Hang et al., in [25], where the fault indicators are defined to remove the influence of rotor speed variation and a frequency tracking algorithm based on reference frame transformation is used to extract the fault symptoms. A possible drawback is the dependency of the fault indicator on the PMSM parameters, which may vary with environmental temperature. Boileau [26] considered as a fault indicator the amplitude of the positive voltage sequence while Meinguet [27] used as an conditions, the current phasor in the Clarke plane draws a circular trajectory, while it follows an elliptic one if a stator fault occurs. In PVAs, the fault index is thus defined by a geometrical distortion factor, obtained by currents measurements. Though demonstrating to be excellent for online ITSC FDI, the PVA is not able to locate the shorted phase.

All the current methods have in common the observation that in case of ITSCs, the stator currents are not balanced, and higher harmonic components in the signals can be used as fault symptoms. Different signal processing methods are employed to extract these symptoms (or patterns). Frequency-domain transforms such as Fast-Fourier Transforms (FFT) are applied in [13,32–34], while mixed frequency/time-domain transforms such as Short-Time Fourier Transform (STFT) are used in [14], Wavelet Transform (WT) in [15,35,36], and Hilbert–Huang Transform (HHT) in [37,38]. Frequency-domain transforms entail the loss of transient events such as speed or loads variations, so FFT analyses, though effective in stationary conditions, must be avoided if the FDI is required during transients. By requiring a higher computational effort, STFT analyses partially compensate the FFT limits, but FDI capabilities are still limited in unsteady conditions. On the other hand, WT methods permit to decompose time domain signals into stationary and non-stationary contributions, but the appropriate choice of wavelets is a drawback. Many limits of WT techniques are removed by HHT techniques, even if they are effective when applied to signals characterised by a narrow frequency content.

Other methods, based on artificial intelligence applications (Convolutional Neural Networks are used in [39]) have been also explored, but, though the results are remarkable, these techniques are not preferred for airworthiness certification [40].

Within this research context, this paper aims to develop an innovative Fault-Tolerant Control (FTC), capable of prompt FDI and accommodation of ITSC faults in three-phase PMSMs for lightweight UAVs propulsion, in both stationary and transient conditions. The proposed FTC is composed of two sections: the first one addresses the FDI problem via an original current-based method (here referred as Advanced-PVA, APVA), the second is dedicated to the fault accommodation, obtained by activating the stand-by leg of a four-leg converter and controlling the central point of the motor Y-connection. It is worth noting that the ITSC accommodation is developed by applying the Rotor Current Frame Transformation Control (RCFTC) technique, already applied by the authors in [9] for the FTC of open-circuit faults.

The paper is articulated as follows: in the first part, the nonlinear model of the propulsion system is presented; successively, the FDI algorithms and the fault accommodation technique are described by highlighting their basic design criteria; finally, a summary of the simulation results is proposed, by characterizing the failure transients and the post-fault behaviours after the injection of ITSC faults.

## **2. Materials and Methods**

## *2.1. System Description*

The reference propulsion system, designed for the full-electric propulsion of a lightweight UAV, is composed of (Figure 3):

	- # Three-phase surface-mounted PMSM with phase windings in Y-connection;
	- # Twin-blade fixed-pitch propeller [41];
	- # Mechanical coupling joint.
	- # CONtrol/MONitoring (CON/MON) module, for the implementation of the closed-loop control and health-monitoring functions;
	- # Four-leg converter;
	- # Three current sensors (CSa, CSb, CSc), one per each motor phase;
	- # One Angular Position Sensor (APS), measuring the motor angle;

o Three current sensors (CSa, CSb, CSc), one per each motor phase; o One Angular Position Sensor (APS), measuring the motor angle;

*Aerospace* **2022**, *8*, x FOR PEER REVIEW 5 of 25

o Four-leg converter;

The CON module operates the closed-loop control of the system by implementing two nested loops on propeller speed and motor currents (via Field-Oriented Control, FOC). All the regulators implement proportional/integral actions on the tracking error signals, plus anti-windup functions with back-calculation algorithms to compensate for commands saturation. The MON module executes the health-monitoring algorithms, including the FTC proposed in this work. The CON module operates the closed-loop control of the system by implementing two nested loops on propeller speed and motor currents (via Field-Oriented Control, FOC). All the regulators implement proportional/integral actions on the tracking error signals, plus anti-windup functions with back-calculation algorithms to compensate for commands saturation. The MON module executes the health-monitoring algorithms, including the FTC proposed in this work.

### *2.2. Model of the Aero-Mechanical Section 2.2. Model of the Aero-Mechanical Section*

The dynamics of the aero-mechanical section of the propulsion system, providing the UAV with the thrust is schematically depicted in Figure 4a and is modelled by [7,9]: The dynamics of the aero-mechanical section of the propulsion system, providing the UAV with the thrust is schematically depicted in Figure 4a and is modelled by [7,9]:

$$\begin{cases} I\_p \ddot{\theta}\_p = -Q\_p - \mathbb{C}\_{gb} \left( \dot{\theta}\_p - \dot{\theta}\_m \right) - \mathbb{K}\_{gb} \left( \theta\_p - \theta\_m \right) + Q\_d\\ I\_m \ddot{\theta}\_m = Q\_m + \mathbb{C}\_{gb} \left( \dot{\theta}\_p - \dot{\theta}\_m \right) + \mathbb{K}\_{gb} \left( \theta\_p - \theta\_m \right) + Q\_c\\ Q\_p = \mathbb{C}\_{Q\_p} \left( \dot{\theta}\_p, AR \right) \rho D\_p^5 \dot{\theta}\_p^2\\ AR = \mathcal{V}\_d / D\_p \theta\_p\\ Q\_c = Q\_{c\max} \sin(n\_h n\_d \theta\_m) \end{cases} \tag{1}$$

the motor, respectively, is the propeller-resistant torque, is a gust-induced disturbance torque, is the motor torque, , is the cogging torque, and is the maximum cogging torque, is the pole pairs number, is the harmonic index of the cogging disturbance, is the nondimensional torque coefficient of the propeller, is the propeller advance ratio, is the propeller diameter, is the air density, is the UAV forward speed, while and are the stiffness and the damping of the mechanical coupling joint. where *J<sup>p</sup>* and *θp*, *Jm*, and *θ<sup>m</sup>* are the inertia and the angular position of the propeller and the motor, respectively, *Q<sup>p</sup>* is the propeller-resistant torque, *Q<sup>d</sup>* is a gust-induced disturbance torque, *Q<sup>m</sup>* is the motor torque, *Qc*, is the cogging torque, and *Qcmax* is the maximum cogging torque, *n<sup>d</sup>* is the pole pairs number, *n<sup>h</sup>* is the harmonic index of the cogging disturbance, *CQ<sup>p</sup>* is the nondimensional torque coefficient of the propeller, *AR* is the propeller advance ratio, *D<sup>p</sup>* is the propeller diameter, *ρ* is the air density, V*<sup>a</sup>* is the UAV forward speed, while *Kgb* and *Cgb* are the stiffness and the damping of the mechanical coupling joint.

where and , , and are the inertia and the angular position of the propeller and

**Figure 4.** Full electric propulsion system: (**a**) mechanical scheme; (**b**) three-phase PMSM schematics with shorted turns in phase a (one pole pair and accessible neutral point). **Figure 4.** Full electric propulsion system: (**a**) mechanical scheme; (**b**) three-phase PMSM schematics with shorted turns in phase a (one pole pair and accessible neutral point).

### *2.3. Model of the PMSM with ITSC fault 2.3. Model of the PMSM with ITSC Fault*

The mathematical modelling of PMSM with ITSC faults that will be used for this work has been previously presented and applied in literature, demonstrating it to be very accurate [13,42–44]. The mathematical modelling of PMSM with ITSC faults that will be used for this work has been previously presented and applied in literature, demonstrating it to be very accurate [13,42–44].

The occurrence of a short circuit in a stator phase causes an asymmetry in the motor magnetic flux, due to an additional circuit with insulation resistance , in which the The occurrence of a short circuit in a stator phase causes an asymmetry in the motor magnetic flux, due to an additional circuit with insulation resistance *R<sup>f</sup>* , in which the shorted current *i <sup>f</sup>* flows, Figure 4b:

shorted current flows, Figure 4b: Once indicated, the fault extension along the winding with the parameter = / (defined as the ratio of the number of shorted turns to the total one), the winding affected by ITSC can be split into an undamaged and a damaged part, Figure 4b. Hence, the electrical equations can be written as: Once indicated, the fault extension along the winding with the parameter *µ* = *N<sup>f</sup>* /*N* (defined as the ratio of the number of shorted turns to the total one), the winding affected by ITSC can be split into an undamaged and a damaged part, Figure 4b. Hence, the electrical equations can be written as:

$$\mathbf{V}\_{\rm xn} = \mathbb{R}\mathbf{I} + \mathbb{L}\frac{d}{dt}\mathbf{I} + \mathbf{e}.\tag{2}$$

 = ℝ + <sup>+</sup> . (2) In Equation (2), = [ − , − , − , 0 ] is the applied voltages vector, = , , , is the currents vector, e is the BEMF vector, which can be expressed as in Equation (3) (in Equation (3) and in the following equations, "s" and "c" briefly indicates "sine" and "cosine" functions): In Equation (2), *Vxn* = [*V<sup>a</sup>* − *Vn*, *V<sup>b</sup>* − *Vn*, *V<sup>c</sup>* − *Vn*, 0 ] T is the applied voltages vector, *I* = h *ia*, *i<sup>b</sup>* , *ic*, *i f* iT is the currents vector, e is the BEMF vector, which can be expressed as in Equation (3) (in Equation (3) and in the following equations, "s" and "c" briefly indicates "sine" and "cosine" functions):

$$\mathbf{e} = \begin{bmatrix} (1 - \mu)e\_a \\ e\_b \\ e\_c \\ \mu e\_a \end{bmatrix} = -\lambda\_m \dot{\theta}\_e \begin{bmatrix} (1 - \mu)\mathbf{s}(\theta\_\varepsilon) \\ \mathbf{s}(\theta\_\varepsilon - 2/3\pi) \\ \mathbf{s}(\theta\_\varepsilon + 2/3\pi) \\ \mu \mathbf{s}(\theta\_\varepsilon) \end{bmatrix} \tag{3}$$

where = is the electrical angle of the motor and is the rotor magnet flux linkage. The resistance and inductance matrixes are: where *θ<sup>e</sup>* = *ndθ<sup>m</sup>* is the electrical angle of the motor and *λ<sup>m</sup>* is the rotor magnet flux linkage. The resistance and inductance matrixes are:

$$\mathbb{R} = \begin{bmatrix} (1-\mu)\mathbb{R} & 0 & 0 & \mathbb{R}\_f \\ 0 & \mathbb{R} & 0 & 0 \\ 0 & 0 & \mathbb{R} & 0 \\ \mu\mathbb{R} & 0 & 0 & -\left(\mu\mathbb{R} + \mathbb{R}\_f\right) \end{bmatrix}, \quad \mathbb{L} = \begin{bmatrix} (1-\mu)^2 L & 0 & 0 & 0 \\ 0 & L & 0 & 0 \\ 0 & 0 & L & 0 \\ \mu^2 L & 0 & 0 & -\mu^2 L \end{bmatrix}, \tag{4}$$

tion resistance, given by:

in which *R* is the phase resistance, *L* is the phase self-inductance, and *R<sup>f</sup>* is the insulation resistance, given by:

$$R\_f = k\_{R\_f} R (1 - \mu) \,\tag{5}$$

with *kR<sup>f</sup>* a factor depending on the insulation material.

By considering an inter-turn fault, the PMSM torque can be calculated as [44]:

$$Q\_m = \left( (1 - \mu)e\_a i\_a + e\_b i\_b + e\_c i\_c - \mu e\_a i\_f \right) / \dot{\theta}\_{m\nu} \tag{6}$$

Since the neutral point voltage *V<sup>n</sup>* is not null when a short circuit occurs, it is convenient to reformulate the electrical equations by extrapolating *Vn*. Hence, by applying Kirchhoff laws to the circuit in Figure 4b and by substituting in Equation (2), we have:

$$\mathbf{V} = \mathcal{R}\mathbf{I} + \mathcal{L}\frac{d}{dt}\mathbf{I} + \mathbf{E},\tag{7}$$

where the reformulated BEMF vector *E* is:

$$\mathbf{E} = \frac{1}{3} \begin{bmatrix} 2e\_{a\_f} - (e\_b + e\_c) \\ 2e\_b - \left(e\_{a\_f} + e\_c\right) \\ 2e\_c - \left(e\_{a\_f} + e\_b\right) \\ 3e\_f \end{bmatrix}'\tag{8}$$

The new resistance and inductance matrixes are:

$$
\mathcal{R} = \begin{bmatrix}
(1 - 2/3\mu)\mathcal{R} & 0 & 0 & 2/3\mathcal{R}\_f \\
\mu\mathcal{R}/3 & \mathcal{R} & 0 & -\mathcal{R}\_f/3 \\
\mu\mathcal{R}/3 & 0 & \mathcal{R} & -\mathcal{R}\_f/3 \\
\mu\mathcal{R} & 0 & 0 & -\left(\mu\mathcal{R} + \mathcal{R}\_f\right)
\end{bmatrix}, \quad \mathcal{L} = \begin{bmatrix}
(1 - 4/3\mu + 2/3\mu^2)L & 0 & 0 & 0 \\
(2 - \mu)\mu L/3 & L & 0 & 0 \\
(2 - \mu)\mu L/3 & 0 & L & 0 \\
\mu^2 L & 0 & 0 & -\mu^2 L
\end{bmatrix}, \tag{9}
$$

while *V* = [*Va*, *V<sup>b</sup>* , *Vc*, 0 ] T is a vector defining the voltage commands sent by the converter to the motor terminals (first three components) and the voltage on the shorted turns. Finally, since the PMSM is controlled via FOC technique, it is convenient to express the voltage vector via the Clarke–Park transformation, as a function of direct and quadrature voltages *V<sup>d</sup>* and *Vq*, generated by the closed-loop control algorithms:

$$
\begin{bmatrix} V\_a \\ V\_b \\ V\_c \\ V\_n \end{bmatrix} = \sqrt{\frac{2}{3}} \begin{bmatrix} \mathbf{c}(\theta\_\ell) & -\mathbf{s}(\theta\_\ell) \\ \mathbf{c}(\theta\_\ell - 2/3\pi) & -\mathbf{s}(\theta\_\ell - 2/3\pi) \\ \mathbf{c}(\theta\_\ell + 2/3\pi) & -\mathbf{s}(\theta\_\ell + 2/3\pi) \\ 0 & 0 \end{bmatrix} \begin{bmatrix} V\_d \\ V\_q \end{bmatrix} \tag{10}
$$

## **3. Fault-Tolerant Control System**

*3.1. FDI Algorithm Conceptualization*

When a PMSM operates at constant speed without faults, the current phasor in the Clarke plane draws a circular trajectory, while it follows an elliptical trajectory when an ITSC fault occurs [28–31]. Torque oscillations appear, at the motor level, with amplitudes that depend on the short-circuit extension.

This section aims to demonstrate that the geometrical parameters of this elliptical trajectory (major, minor axes length, and axes inclination) are univocally related to the location and extension of the ITSC. In particular, we will demonstrate that:


It is worth noting that, during constant speed operations, any elliptical trajectory of the motor current phasor in the Clarke plane (including the circular one, as a special case) can be reconstructed via the Fortescue decomposition [45], as the sum of three symmetrical and balanced rotating systems: *positive*-, *negative*-, and *zero-sequence* systems. *Aerospace* **2022**, *8*, x FOR PEER REVIEW 8 of 25

> When an ITSC fault occurs, the central point of the Y-connection is isolated, so the *zerosequence* current is zero, and the ellipse tracked by the current phasor (ˆ*I*) can be considered as the sum of two counter-rotating phasors: *positive*- (ˆ*I* <sup>+</sup>) and *negative-sequence* phasors ( ˆ*I* −), as given by: When an ITSC fault occurs, the central point of the Y-connection is isolated, so the *zero-sequence* current is zero, and the ellipse tracked by the current phasor () can be considered as the sum of two counter-rotating phasors: *positive-* () and *negative-sequence* phasors (), as given by:

$$\hat{I} = |\hat{I}(t)|e^{j\varphi(t)} = |\hat{I}^+|e^{j(\dot{\theta}\_\ell t + \varphi^+)} + |\hat{I}^-|e^{-j(\dot{\theta}\_\ell t - \varphi^-)}\,\,\,\,\,\,\tag{11}$$

where *ϕ*, *ϕ* <sup>+</sup>, and *ϕ* − are the phase angles of the resultant, *positive*-, and *negative*-*sequence* current phasors, respectively. As illustrated in Figure 5, the semi major and semi minor axis lengths are obtained by: Where , , and are the phase angles of the resultant, *positive-,* and *negative*-*sequence* current phasors, respectively. As illustrated in Figure 5, the semi major and semi minor axis lengths are obtained by:

$$\mathbf{s}\_{M,m} = |\mathbf{f}^+| \pm |\mathbf{f}^-|. \tag{12}$$

**Figure 5.** Decomposition of the current phasor in the Clarke plane into *positive-* and *negative-sequence* components. **Figure 5.** Decomposition of the current phasor in the Clarke plane into *positive*- and *negativesequence* components.

While the inclination of the ellipse major axis corresponds to the angle that maximizes the current phasor amplitude, up to: While the inclination of the ellipse major axis corresponds to the angle that maximizes the current phasor amplitude, up to:

$$\max\left(|\hat{I}|\right) = |\hat{I}^+| + |\hat{I}^-|.\tag{13}$$

The condition in Equation (13) is reached when the *positive-* and *negative-sequence* phasors have the same angular orientation: = ̇ |(||) + = −̇ |(||) + , (14) The condition in Equation (13) is reached when the *positive*- and *negative-sequence* phasors have the same angular orientation:

$$\left.\dot{\varphi}\_{\varepsilon l} = \dot{\theta}\_{\varepsilon} t\right|\_{\max(\left|\dot{I}\right|)} + \left.\varphi^{+} = -\dot{\theta}\_{\varepsilon} t\right|\_{\max(\left|\dot{I}\right|)} + \left.\varphi^{-}\right|\_{\max(\left|\dot{I}\right|)}\tag{14}$$

which leads to: which leads to:

m leads to:

$$
\varphi\_{el} = \left(\boldsymbol{\varrho}^{-} + \boldsymbol{\varrho}^{+}\right)/2.\tag{15}
$$

*positive-* and *negative-sequence* phasors, expressed by: <sup>∗</sup> = atan <sup>∗</sup> <sup>∗</sup> ⁄ , (16) To identify the ellipse inclination, it is thus crucial to evaluate the phase angles of the *positive*- and *negative-sequence* phasors, expressed by:

By introducing a complex analysis notation, the phasor ∗ can be also represented in

$$\boldsymbol{\sigma}^\* = \operatorname{atan} \left( i\_{\boldsymbol{\beta}}^\* / i\_{\boldsymbol{\alpha}}^\* \right) . \tag{16}$$

terms of a balanced three-phase space vector system, as:

the and axes in the Clarke plane.

where "∗" stands for "+" or "−", and *i* ∗ *<sup>α</sup>* and *i* ∗ *β* are the projections of the phasor ˆ*I* ∗ on the *α* and *β* axes in the Clarke plane.

By introducing a complex analysis notation, the phasor ˆ*I* ∗ can be also represented in terms of a balanced three-phase space vector system, as:

$$\hat{I}^\* = \sqrt{2/3} \left( \text{Re} \left( \hat{I}\_a^\* \right) + \text{Re} \left( \hat{I}\_b^\* \right) e^{j2\pi/3} + \text{Re} \left( \hat{I}\_c^\* \right) e^{j4\pi/3} \right) \tag{17}$$

Then, by expressing the space vector ˆ*I* ∗ into its real (*i* ∗ *α* ) and imaginary (*i* ∗ *β* ) components, we have:

$$\begin{split} \mathbf{i}\_a^\* &= \sqrt{2/3} \left[ \mathrm{Re} \left( \mathbf{\hat{I}}\_a^\* \right) - 1/2 \left( \mathrm{Re} \left( \mathbf{\hat{I}}\_b^\* \right) + \mathrm{Re} \left( \mathbf{\hat{I}}\_c^\* \right) \right) \right] = \sqrt{3/2} \mathrm{Re} \left( \mathbf{\hat{I}}\_a^\* \right), \\ \mathbf{i}\_\beta^\* &= \sqrt{2/3} \left[ \sqrt{3}/2 \left( \mathrm{Re} \left( \mathbf{\hat{I}}\_b^\* \right) - \mathrm{Re} \left( \mathbf{\hat{I}}\_c^\* \right) \right) \right] = \ast - \sqrt{3/2} \mathrm{Im} \left( \mathbf{\hat{I}}\_a^\* \right), \end{split} \tag{18}$$

The phase angles of the *positive*- and *negative-sequence* phasors with ITSC faults can be thus calculated via the *positive*- and *negative-sequence* phasors related to phase *a* only ( ˆ*I* + *<sup>a</sup>* and ˆ*I* − *a* ), given by the Fortescue transform:

$$
\begin{bmatrix}
\hat{I}\_a^+ \\
\hat{I}\_a^- \\
0 \\
\hat{I}\_f
\end{bmatrix} = \frac{1}{3} \begin{bmatrix}
1 & e^{j2\pi/3} & e^{j4\pi/3} & 0 \\
1 & e^{j4\pi/3} & e^{j2\pi/3} & 0 \\
1 & 1 & 1 & 0 \\
0 & 0 & 0 & 3
\end{bmatrix} \begin{bmatrix}
\hat{I}\_a \\
\hat{I}\_b \\
\hat{I}\_c \\
\hat{I}\_f
\end{bmatrix} \tag{19}
$$

The vector of phasors ˆ*I* = h ˆ*Ia*, ˆ*Ib* , ˆ*Ic*, ˆ*If* i can be obtained from:

$$\mathbf{J} = \left(\mathcal{R} + j\dot{\theta}\_{\varepsilon}\mathcal{L}\right)^{-1}(\mathcal{V} - \mathbf{E}),\tag{20}$$

where the voltage and BEMF phasors are:

$$
\hat{\mathbf{V}} = \sqrt{\frac{2}{3}} (V\_d + jV\_q) \begin{bmatrix} 1 \\ \begin{pmatrix} -1/2 - j\sqrt{3}/2 \\ -1/2 + j\sqrt{3}/2 \\ 0 \end{pmatrix} \\\tag{21}
$$

$$
\hat{\mathbf{E}} = j\lambda\_m \dot{\theta}\_e \begin{bmatrix} 1 - 2/3\mu \\ \mu - 1/2 - j\sqrt{3}/2 \\ \mu - 1/2 + j\sqrt{3}/2 \\ \mu \end{bmatrix} \tag{22}
$$

The solution of Equation (20) in terms of current phasors are presented by the polarcoordinate plot in Figure 6a, in which, given along the radius the ITSC fault extension (*µ* = *N<sup>f</sup>* /*N*), the inclination of the ellipse major axis of the current phasor trajectory is reported along the phase angle for each motor phase. In addition, Figure 6b shows the ratio of ellipse semi-axes lengths as a function of the ITSC fault extension parameter *µ*. It is worth noting that, if an ITSC fault occurs on phase *a* (*b*, *c*), the ellipse inclination will be 0 ◦ (240◦ , 120◦ ) ±15◦ .

**Figure 6.** Inclination of the ellipse major axis (**a**) and ellipses axis ratio (**b**) under ITSC faults of dif‐ ferent extensions and locations. **Figure 6.** Inclination of the ellipse major axis (**a**) and ellipses axis ratio (**b**) under ITSC faults of different extensions and locations.

The above‐mentioned discussion thus supports the development of an FDI algorithm based on the estimation of the geometrical characteristics of (generic) elliptical trajectories of current phasors in the Clarke plane, so that: The above-mentioned discussion thus supports the development of an FDI algorithm based on the estimation of the geometrical characteristics of (generic) elliptical trajectoriesof current phasors in the Clarke plane, so that:


The FDI concept essentially entails the identification of geometrical properties of an ellipse. Clearly, at a single monitoring time, a single point of the ellipse is known (ఈ, ఉ Figure 5); therefore, the fitting problem is under‐determined. On the other hand, the prob‐ lem becomes over‐determined if a sufficient number of measurements are used. The over‐ determined problem of fitting an ellipse to a set of data points arises in many applications such as computer graphics [46], hydraulic engineering [47], and statistics [48]. Many dif‐ ferent methods are proposed in literature for fitting the geometric primitives of an ellipse: voting/clustering methods (e.g., the Hough transform [49] and the RANSAC technique [50]) are robust but they have poor accuracy and require large computing resources, opti‐ mization methods and direct least‐square methods [51] are accurate, but, in case of non‐ linear applications such as an ellipse fitting, they apply iterative techniques. The FDI concept essentially entails the identification of geometrical properties of an ellipse. Clearly, at a single monitoring time, a single point of the ellipse is known (*iα*, *i<sup>β</sup>* Figure 5); therefore, the fitting problem is under-determined. On the other hand, the problem becomes over-determined if a sufficient number of measurements are used. The over-determined problem of fitting an ellipse to a set of data points arises in many applications such as computer graphics [46], hydraulic engineering [47], and statistics [48]. Many different methods are proposed in literature for fitting the geometric primitives of an ellipse: voting/clustering methods (e.g., the Hough transform [49] and the RANSAC technique [50]) are robust but they have poor accuracy and require large computing resources, optimization methods and direct least-square methods [51] are accurate, but, in case of nonlinear applications such as an ellipse fitting, they apply iterative techniques.

Here we instead use a direct least‐square method, as proposed by Halìř and Flusser [52], who, starting from the quadratically‐constrained least‐squares minimization postu‐ lated by Fitzgibbon [53], reformulated the problem via partitioning techniques, by obtain‐ ing a very robust and computationally‐effective algorithm. Here we instead use a direct least-square method, as proposed by Halìˇr and Flusser [52], who, starting from the quadratically-constrained least-squares minimization postulated by Fitzgibbon [53], reformulated the problem via partitioning techniques, by obtaining a very robust and computationally-effective algorithm.

Once given the vectorial conic definition of an ellipse, Once given the vectorial conic definition of an ellipse,

$$
\Gamma \cdot \gamma = 0,\tag{23}
$$

in which ൌ ሾଶ, , ଶ, , , 1ሿ and ൌ ሾ, , , , , ሿ are vectors defining the Carte‐ sian coordinates of the ellipse points and the ellipse coefficients, respectively, the over‐ determined fitting problem to a set of coordinate points and (where ൌ 1, … , , in which **Γ** = - *α* 2 , *αβ*, *β* 2 , *α*, *β*, 1 and *γ* = [*A*, *B*, *C*, *D*, *E*, *F*] T are vectors defining the Cartesian coordinates of the ellipse points and the ellipse coefficients, respectively, the overdetermined fitting problem to a set of coordinate points *α<sup>i</sup>* and *β<sup>i</sup>* (where *i* = 1, . . . , *n*, and

and is greater than the number of ellipse coefficients, i.e., >6) can be solved by the

following eigenvalues problem (for more details, see Appendix A):

*n* is greater than the number of ellipse coefficients, i.e., *n* > 6) can be solved by the following eigenvalues problem (for more details, see Appendix A):

$$\begin{cases} \mathbb{M}\gamma\_1 = \lambda\gamma\_1\\ \begin{array}{l} \gamma\_1^\mathrm{T} \mathbb{C}\_1 \gamma\_1 = 1\\ \gamma\_2 = -\mathbb{S}\_3^{-1} \mathbb{S}\_2^\mathrm{T} \gamma\_1\\ \gamma = \left(\gamma\_1 \ \gamma\_2\right)^\mathrm{T} \end{array} \tag{24}$$

where M is the reduced scatter matrix,

$$\mathbb{M} = \mathbb{C}\_1^{-1} \left( \mathbb{S}\_1 - \mathbb{S}\_2 \mathbb{S}\_3^{-1} \mathbb{S}\_2^{\mathrm{T}} \right), \tag{25}$$

$$\mathbb{S}\_1 = \mathbb{D}\_1^\mathrm{T} \mathbb{D}\_1, \quad \mathbb{S}\_2 = \mathbb{D}\_1^\mathrm{T} \mathbb{D}\_2, \quad \mathbb{S}\_3 = \mathbb{D}\_2^\mathrm{T} \mathbb{D}\_2,\tag{26}$$

$$\mathbb{D}\_1 = \begin{bmatrix} \alpha\_1^2 & \alpha\_1 \beta\_1 & \beta\_1^2 \\ \vdots & \vdots & \vdots \\ \alpha\_n^2 & \alpha\_n \beta\_n & \beta\_n^2 \end{bmatrix}, \begin{bmatrix} \mathbb{D}\_2 = \begin{bmatrix} \alpha\_1 & \beta\_1 & 1 \\ \vdots & \vdots & \vdots \\ \alpha\_n & \beta\_n & 1 \end{bmatrix} \end{bmatrix} \tag{27}$$

$$\mathbb{C}\_1 = \begin{bmatrix} 0 & 0 & 2 \\ 0 & -1 & 0 \\ 2 & 0 & 0 \end{bmatrix} \, \tag{28}$$

$$\boldsymbol{\gamma} = \begin{bmatrix} \boldsymbol{\gamma} \ \boldsymbol{\gamma} \end{bmatrix}^{\mathrm{T}} \tag{29}$$

with the ellipse coefficient vector segmented in *γ*<sup>1</sup> = - *A B C*<sup>T</sup> and *γ*<sup>2</sup> = - *D E F*<sup>T</sup> . The solution of Equation (24) corresponds to the eigenvector *γ* yielding a minimal nonnegative eigenvalue *λ*. Once obtained *γ*, the lengths of major and minor semi-axes *s<sup>M</sup>* and *s<sup>m</sup>* are given by [54]:

$$s\_{M,m} = \frac{\sqrt{\frac{2(AE^2 + \mathbb{C}D^2 - BDE + (B^2 - 4AC)F)}{\left(A + \mathbb{C} \pm \sqrt{(A - \mathbb{C})^2 + B^2}\right)^{-1}}}}{4AC - B^2}}{4AC - B^2} \tag{30}$$

while the major axis inclination *ϕel* is [54]:

$$\varphi\_{el} = \begin{cases} 0 & \text{for } \mathcal{B} = 0, \ A < \mathcal{C} \\ \pi/2 & \text{for } \mathcal{B} = 0, \ A > \mathcal{C} \\ 1/2 \cot^{-1}((A - \mathcal{C}) \text{ / } \mathcal{B}) & \text{for } \mathcal{B} \neq 0, \ A < \mathcal{C} \\ \pi/2 + 1/2 \cot^{-1}((A - \mathcal{C}) \text{ / } \mathcal{B}) & \text{for } \mathcal{B} \neq 0, \ A > \mathcal{C} \end{cases} \tag{31}$$

## *3.2. FDI Algorithm Design and Implementation*

One of the basic limitations of the proposed FDI concept is that the technique is defined by assuming constant speed motor operations. During unsteady conditions, the current phasor changes its amplitude while rotating in the Clarke plane. Consequently, the ellipse identification does not work appropriately, and it could cause false alarms.

To overcome the problem, it is worth noting that, differently from the constant speed operations with ITSC faults, the axes of the ellipse reconstructed in unsteady conditions rotate. For this reason, the developed FDI algorithm detects an ITSC fault only if the value of the ellipse major axis is constant and stationary. This FDI logic has been implemented as represented by the flow chart in Figure 7: at each *k*-th monitoring sample, if the Boolean variable *mon*(*k*) is true (Equation (32)), a fault counter *n* (*k*) *c* is increased by two, otherwise it is reduced by one. When the fault counter reaches a predefined maximum value (*nth*), the algorithm outputs a true Boolean fault flag *fx*, which indicates that a fault on the phase *x* (= *a*, *b*, *c*) is detected and isolated.

$$
gamma^{(k)} = \left(\Delta\_d^{(k)} \ge \varepsilon\_{\rm{dtlh}} \land \ \Delta\_i^{(k)} \le \varepsilon\_{\rm{itlh}}\right) \tag{32}
$$

where:

$$\Delta\_d^{(k)} = \left| s\_M^{(k)} - s\_m^{(k)} \right|,\tag{33}$$

$$
\Delta\_i^{(k)} = \min \left( \left. \Delta\_i^{(k)} \right|\_{a'} \left. \Delta\_i^{(k)} \right|\_{c'} \left. \Delta\_i^{(k)} \right|\_{b} \right), \tag{34}
$$

$$
\left| \Delta\_i^{(k)} \right|\_a = \left| \varphi\_{el}^{(k)} \right|\_\prime \ \Delta\_i^{(k)} \Big|\_b = \left| \varphi\_{el}^{(k)} \right| - \frac{2\pi}{3} \ \ \Delta\_i^{(k)} \Big|\_c = \left| \varphi\_{el}^{(k)} \right| + \frac{2\pi}{3} \ \tag{35}
$$

**Figure 7.** Flow chart for the FDI logic. **Figure 7.** Flow chart for the FDI logic.

where: Δ () <sup>=</sup> () <sup>−</sup> () , (33) The details and a dedicated discussion about the FDI parameters design are proposed in Section 4, while tables reporting the FDI parameters are reported in Appendix B.

## *3.3. Fault Accommodation Algorithm*

Δ () <sup>=</sup> min Δ () , Δ () , Δ () , (34) Δ () = () , Δ () = () <sup>−</sup> 2 <sup>3</sup> , Δ () = () + 2 <sup>3</sup> , (35) The details and a dedicated discussion about the FDI parameters design are proposed in Section 4, while tables reporting the FDI parameters are reported in Appendix B. After the FDI algorithm detects and isolates the ITSC fault, the motor still has two active phases. To maintain the performances after the fault, the FTC should accommodate the system by restoring the operation of the current phasor in the Clarke plane as it was in healthy conditions. This is here achieved by enabling the control of the central point of the motor Y-connection, via the stand-by leg of the four-leg converter, and by applying an RCFTC technique [9,55].

*3.3. Fault Accommodation Algorithm* After the FDI algorithm detects and isolates the ITSC fault, the motor still has two active phases. To maintain the performances after the fault, the FTC should accommodate the system by restoring the operation of the current phasor in the Clarke plane as it was To maintain the torque performances with only two phases, the currents amplitudes in the healthy phases in the rotor reference frame (*i* # *x f* , *i* # *y f*) must increase by <sup>√</sup> 3 w.r.t. those before the fault, and they must shift by 60◦ along the electrical angle. In addition, the amplitude of the current flowing into the central point (*i* # *n* ) must be <sup>√</sup> 3 times those in the healthy phases:

$$\begin{cases} i\_{xf}^{\#} = \sqrt{2} \left( i\_d^{\#} \text{c} \left[ \theta\_\varepsilon + \frac{2\pi}{3} \left( m + \frac{7}{4} \right) \right] - i\_q^{\#} \text{s} \left[ \theta\_\varepsilon + \frac{2\pi}{3} \left( m + \frac{7}{4} \right) \right] \right) \\ i\_{yf}^{\#} = \sqrt{2} \left( i\_d^{\#} \text{c} \left[ \theta\_\varepsilon + \frac{2\pi}{3} \left( m + \frac{5}{4} \right) \right] - i\_q^{\#} \text{s} \left[ \theta\_\varepsilon + \frac{2\pi}{3} \left( m + \frac{5}{4} \right) \right] \right) \\\ i\_{wf}^{\#} = 0 \\\ i\_n^{\#} = \sqrt{6} \left( i\_d^{\#} \text{c} \left( \theta\_\varepsilon + \frac{2\pi}{3} m \right) - i\_q^{\#} \text{s} \left( \theta\_\varepsilon + \frac{2\pi}{3} \right) \right) \end{cases} \tag{36}$$
