Toward Smart Air Mobility: Control System Design and Experimental Validation for an Unmanned Light Helicopter

Abstract

Light helicopters are used for a variety of applications, attracting users from private and public market segments because of their agility and convenient storage capabilities. However, most light helicopters on the market today are designed and manufactured with technologies dating back to the 1980s, with safety issues to be addressed by advanced design methods, more powerful engines, and innovative solutions. In this regard, the DISRUPT (Development of an innovative and safe ultralight, two-seater turbine helicopter) project, led by Curti Aerospace Division (Italy) and co-funded by the EU H2020 program, is a state-of-the-art concept for a novel ultralight helicopter equipped with a ballistic parachute. In order to validate the first parachute ejection in a safe scenario, a dronization process was selected as a viable solution to be performed in collaboration with the University of Bologna. In the present paper, the steps followed to transform the helicopter into an unmanned vehicle are detailed according to the model-based design approach, with particular focus on mathematical modeling, control system design, and experimental validation. Obtained results demonstrate the feasibility of using a civil helicopter first as a remotely-piloted vehicle and then as a highly-automated personal transportation system in the framework of smart and sustainable air mobility.

1. Introduction

The interest in Urban Air Mobility (UAM) had a step increase over the last few years [1]. On the one hand, the slow growth rate of ground infrastructure led to critical traffic congestion in urban areas. On the other hand, the increasing demand for moving people and payloads further and faster drove the attention of the research community and stakeholders toward the exploitation of the vertical dimension [2]. For example, Amazon and Google pioneered the testing of urban parcel delivery by means of multirotor aircraft [3,4]. In such a way, they paved the way for a wide range of studies on highly-automated low-altitude vehicles as an alternative means of transportation, where “the regular Joe” is capable of performing a mission without having the skills of a licensed pilot [5,6,7]. In this respect, two early attempts that investigated concepts of operation and technologies for a new personal transportation system based on both an aerial platform and a ground infrastructure were, respectively, PPlane (2009–2013) and myCopter (2011–2014), projects funded by the European Commission under the 7th Framework Program (FP7) [8,9].

By taking advantage of consolidated experience in conventional aviation, high reliability of onboard systems, and rapid improvement of electrical propulsion performance, manufacturers and transport stakeholders (such as Airbus, Volocopter, and Uber) investigated concepts for personal air transportation systems. With the aim of playing a lead role in this new raising market, they considered electric platforms with Vertical Take-Off and Landing (VTOL) capabilities as key elements for the next generation of controlled airspace [10,11].

Among all the above-mentioned projects and applications, it is acknowledged that a cost-effective solution to sustainable Urban Air Mobility and Delivery (UAMD) is represented by the use of small/light aircraft, where onboard flight control systems, supported by Air Traffic Management (ATM) technology, will provide safe navigation in dynamic scenarios and weather conditions in the presence of other sky users [12]. Transforming a conventional aircraft (both fixed and rotary-wing) into a Remotely–Piloted Aerial System (RPAS) may represent a successful strategy for different reasons. First of all, available light/ultralight conventional aircraft have already passed through several design, test, and certification steps with the aim of fulfilling reliability, performance, and flying quality requirements [13]. Moreover, reversible control chains can be easily replaced by Electro-Mechanical Actuators (EMA), controlled by dedicated onboard avionics. Starting from this design bias, researchers can thus focus on the design and experimental validation of all other technologies allowing for UAMD (including Guidance, Navigation, and Control (GNC) systems, telemetry, communication, and ATM devices) in addition to ground handling facilities. In this respect, thanks to their compact size and peculiar VTOL configuration, civil ultralight helicopters represent suitable test-beds for performing the transition toward a highly-automated personal transportation system.

By focusing on the very recent past, examples of the transition of conventional helicopters into RPASs can be dated back to 2004, when the Unmanned Little Bird demonstrator, derived by Boeing from a civil MD 530F, made its first autonomous flight (with a safety pilot). In particular, a pre-programmed 20-min armed intelligence, surveillance, and reconnaissance mission was performed around the United States Army’s Yuma Proving Ground facility [14]. In 2006, Northrop Grumman introduced the MQ–8 Fire Scout unmanned helicopter family, obtained from Schweizer 333 and Bell 407, designed to provide reconnaissance, situational awareness, aerial fire, and precision targeting support for ground, air, and sea forces [15]. In 2008, an unmanned, highly-automated version of the Kaman K-MAX helicopter took its maiden flight, with the aim of operating in combat scenarios as well as in civilian situations involving chemical, biological, or radiological hazards [16]. Later on, Eurocopter launched a series of flights for a new rotary-wing solution designed to expand the mission capabilities of Eurocopter helicopters [17]. The Optionally–Piloted Vehicle (OPV) program, based on the EC145 helicopter platform (now Airbus Helicopters H145), was revealed during a demonstration flight: after an automatic takeoff, an EC145 flew a circuit via pre-programmed waypoints and performed a mid-route hover to deploy a load from the external sling. The EC145 continued on a return route segment representing a typical observation mission, followed by an automatic landing. Finally, Sikorsky demonstrated its OPV Matrix Technology on a modified S–76B helicopter called the Sikorsky Autonomy Research Aircraft (SARA). Since 2013, the program has made progress with more than 300 h of autonomous flight with the aim of improving decision-aiding for manned operations, while enabling both unmanned and reduced-crew operations [18].

This paper presents the results of a research work performed within DISRUPT (2016–2018), a collaborative project co-funded by the EU within the H2020 program and led by Curti Aerospace Division. Specifically, DISRUPT proposed a new light rotorcraft configuration, the two-seater Curti Zefhir helicopter, that features a turbine engine and an emergency ballistic parachute to respectively enhance flight performance and increase passenger safety (see Figure 1). PBS Velká Bíteš manufactures the turboshaft engine, derated from 160 to 105 kW of maximum continuous power. While ballistic parachutes have been certified on some fixed-wing aircraft, such as Cirrus light airplanes, their installation on helicopters is a challenging proposition due to the overhead presence of rotating blades. Contained in a non-rotating pod above the main rotor, the parachute solution proposed by Curti and Junkers ProFly thus becomes a backup for conditions where autorotation cannot be performed, such as (a) flight control failure or loss of maneuverability, (b) flying over an area where emergency landing cannot be safely performed, or (c) flight conditions that prevent restoring rotor rotation speed [19].

Although the main objectives of DISRUPT were not strictly related to the main topics of UAM, the need for a remotely-piloted configuration arose immediately; since the experimental validation of the parachute system with the full-scale helicopter was one of the main expected results, the transition toward an unmanned configuration became a mandatory activity to perform the ejection test without a human pilot on board. A crucial but challenging step of the process was the design of a stabilization system, intended as a flexible and reliable software/hardware solution allowing the pilot to manage the ejection task while reducing the workload required by control action. Helicopters generally show nonlinear, complex dynamics that might manifest some unstable flight characteristics in limited zones of the flight envelope. In the particular case of a radio-controlled rotorcraft, without the direct perception of linear accelerations and attitude motion, the remote piloting of a helicopter is indeed an extremely hazardous task [20,21]. Hence, an Automatic Flight Control System (AFCS) was designed, tested, and implemented, allowing the pilot to safely control the aircraft in terms of desired attitude.

The main goal of the paper is to present for the first time a detailed description of all the phases allowing the successful transition of a conventional light helicopter into a RPAS while investigating the validity of a rescue system in the framework of future UAM applications. According to the Model-Based Design (MBD) philosophy, (1) mission requirements are listed and (2) system architecture is defined. Furthermore, (3) an accurate 6DOF nonlinear model is implemented in the Matlab/Simulink environment, which includes helicopter subsystems, environmental effects, and sensor and actuator behavior. (4) The mathematical model is validated and refined by using flight data collected during an identification campaign. (5) After the analysis of open loop dynamic modes, (6) an attitude control system allowing the remote pilot to easily control the aircraft is designed, implemented, and validated by means of both (7) Hardware–In–the–Loop (HIL) techniques and (8) flight tests.

The paper is structured as follows. Section 2 addresses the outline of mission requirements and the selection of system components. The entire simulation model, the trim and stability analysis, and the model validation procedure are presented in Section 3 and Section 4, respectively. Control system design, implementation, and HIL validation are described in Section 5. Experimental results validating the AFCS performance and reporting the parachute recovery mission are finally summarized in Section 6. A section of concluding remarks ends this paper.

The successful outcome of the ejection test and the interest that has arisen in several journals and broadcast media prove the relevance of the research activity presented in this paper [22,23]. Zefhir is currently the only civil helicopter equipped with a ballistic parachute. Indeed, such a test has never been filmed or documented in the entire history of aerospace technology. However, due to the highly-classified nature of the data involved in the early stages of aircraft development, a detailed description of helicopter features and both numerical and experimental results is omitted in the present framework. The focus of the analysis is thus placed on the description of methodological aspects, with particular attention to both numerical and experimental validations, supported by results available in the literature. Furthermore, the comparison between experimental data and the results of simulations is possibly characterized in terms of relative errors, while the description of the technological setup is circumscribed to functional aspects. The uniqueness of the experiment and the absence of strict performance requirements finally vindicate the limits posed by the novelty of the proposed control approach. In this respect, the necessity to rapidly design a safe single-case ejection test necessarily restricts the degree of experimentation, driving the MBD workflow to focus on long-standing results in the field of PID control. Although the latter does not guarantee optimality, it takes advantage of (1) a reduced number of involved parameters; (2) simple implementation and low computational cost; (3) the possibility to perform dedicated flight tests aiming at characterizing the closed-loop dynamic behavior one axis at a time while evaluating the effects of single gain contribution; and (4) an intuitive sizing procedure, suitable for collaboration with the candidate pilot to pursue a set of prescribed handling qualities. Alternative control techniques, such as robust nonlinear and adaptive control that involve the stabilization of vehicle speed components, are currently under experimental validation by the authors, provided small-scale rotorcraft are adopted as test beds in the direction of safe, scalable, and high-performance air mobility and delivery scenarios [24].

2. Mission Requirements and System Architecture

2.1. Mission Requirements

Mission systems and subsystems are grouped into the ground segment and the flight segment:

• Ground segment or Ground Control Station (GCS): the complete set of ground-based systems used to control and monitor the flight segment. The main components include the human-machine interface, computer, telemetry, and aerials for the control, video, and data link to and from the unmanned vehicle.
• Flight segment: the helicopter is equipped with the necessary avionics to perform a remotely–piloted flight. The main components include sensors, actuators for rotor blade pitch angle control, an onboard computer, and aerials for the control, video, and data links to and from the ground segment.

The final mission is defined by the following phases (Figure 2):

• Pre-flight checks: the systems involved in the mission are prepared and visually checked. The helicopter is placed on flat terrain at a safety distance from the GCS. The airfield is required to be clear of obstacles while the mission airspace is circumscribed by a radius of 5 km and a height of 500 m with respect to the GCS.
• Avionics power-on: both the ground and the flight segment subsystems are activated. Telemetry data are received by the GCS, and software/hardware verification checks are performed. The pilot validates the correct actuation of control commands.
• Engine start: the ignition procedure is started by the pilot’s action and the turbine reaches the idle condition.
• Take-off and climb: the helicopter takes-off and climbs out of ground effect at a controlled rate until reaching 300 m above the airfield.
• Cruise: the helicopter is stabilized in steady level flight at about 30 kts.
• Engine shutdown and parachute ejection: the pilot performs the termination procedure, which includes engine shutdown and parachute ejection.
• Descent: the helicopter descends with a stabilized speed and lands within the prescribed area.

Figure 2. Mission phases definition.

Figure 2. Mission phases definition.

The mission is performed in the visual line of sight. However, telemetry information needs to be available to both the GCS crew and the pilot. Given the intrinsic dynamic instability of the helicopter, fuselage attitude stabilization algorithms are required to assist the pilot throughout the mission profile. Conversely, no closed-loop control is applied to the MR collective pitch. The use of an independent Flight Termination System (FTS) is mandatory to stop the engine in case of emergency.

2.2. System Architecture

The selection of components for both the ground and flight segments is performed on the basis of mission requirements. At the same time, the MBD approach is adopted to define systems and subsystems as the result of an iterative process, where the making of a simulation model represents the core of control system development (see Section 3 and Section 5). In what follows, the equipment list is presented, while the unmanned system layout is sketched in Figure 3.

The GCS is made of:

• the control module, where a modified commercial-off-the-shelf Radio Controller (RC) is used as a human-machine interface. Commands from the pilot, represented by stick deflections and switch activation inputs, are generated as Pulse-Width Modulated (PWM) signals and collected via the Pulse-Position Modulation (PPM) protocol. The PPM signal is finally provided to an integrated micro-controller board and output to the communication module according to serial protocol;
• the monitoring module, represented by a rugged laptop, where a graphical user interface is designed to display telemetry data, plan the mission, and send high-level commands via an Ethernet TCP/IP connection to a Real-Time Computer (RTC1) for data acquisition and processing;
• the communication module, which provides an RX/TX radio link to the flight segment. An ethernet switch is used to collect data from the monitoring module, while a ground-based radio modem is connected to a pair of 8 dBi $2.4$ GHz directional patch antennas (respectively characterized by right-hand circular and vertical polarization).

The helicopter is equipped with:

• a corresponding radio modem. Data are output via serial protocol and converted to a widespread standard industrial bus for communication with the Flight Management System (FMS);
• a Real-Time Computer (RTC2) performing FMS data acquisition and control tasks;
• a combined navigation and Attitude and Heading Reference System (AHRS) to estimate attitude information in a dynamic environment, along with position and velocity. Data are output via serial protocol and converted to a standard industrial bus for communication with the FMS;
• a set of 4 EMAs controls the collective, lateral, and longitudinal blade pitches of the MR and the collective pitch of the TR. An additional EMA is used for parachute deployment actuation. An FTS, based on a separate 868 MHz radio system, allows the Fuel Shut-Off Valve (FSOV) to close for emergency engine shutdown. The EMA and FTS selected for the experiment are devices available in the civil market.

3. System Modeling

Starting from the definition of reference frames, a 6 degrees-of-freedom model is adopted to represent the helicopter, with general expressions for the kinematics and dynamics of a rigid body with a center of gravity $CG$.

3.1. Reference Frames

Three right-handed orthogonal reference frames are introduced, according to the definitions in [25]:

• an Earth-fixed North-East-Down frame, ${\mathcal{F}}_E=\{O_E;{\mathit{x}}_E,{\mathit{y}}_E,{\mathit{z}}_E\}$: the origin, $O_E$, is arbitrarily fixed to a point on the Earth’s surface, ${\mathit{x}}_E$ aims in the direction of the geodetic North, ${\mathit{z}}_E$ points downwards along the Earth’s ellipsoid normal, and ${\mathit{y}}_E$ completes a right-handed triad. This frame is assumed to be inertial under the assumption of a flat and non-rotating Earth;
• a Local Vertical-Local Horizontal frame, ${\mathcal{F}}_H=\{CG;{\mathit{x}}_H,{\mathit{y}}_H,{\mathit{z}}_H\}$: the origin is located at the vehicle’s center of gravity, $CG$. Under the hypothesis of a flat Earth, ${\mathcal{F}}_H$ has axes parallel to ${\mathcal{F}}_E$;
• a body-fixed frame, ${\mathcal{F}}_B=\{CG;{\mathit{x}}_B,{\mathit{y}}_B,{\mathit{z}}_B\}$: the ${\mathit{x}}_B$–axis is positive out the nose of the rotorcraft in its plane of symmetry, ${\mathit{z}}_B$ is perpendicular to ${\mathit{x}}_B$ in the same plane of symmetry, pointing downwards, and ${\mathit{y}}_B$ completes a right-handed triad;
• an aircraft reference frame, ${\mathcal{F}}_A=\{O_A;{\mathit{x}}_A,{\mathit{y}}_A,{\mathit{z}}_A\}$, used to locate $CG$ and all helicopter components: axes are parallel to the body-fixed frame axes, such that ${\mathit{x}}_A=-{\mathit{x}}_B$, ${\mathit{y}}_A={\mathit{y}}_B$, and ${\mathit{z}}_A=-{\mathit{z}}_B$. The origin is located ahead and below the rotorcraft at some arbitrary point within the plane of symmetry. Stations ($ST$) are measured positive aft along the longitudinal axis. Buttlines ($BL$) are lateral distances, positive to the pilot’s right, and waterlines ($WL$) are measured vertically, positive upwards. A sketch of the rotorcraft, including the selected ${\mathcal{F}}_A$ frame, is reported in Figure 4. The positions of the main components, expressed in ${\mathcal{F}}_A$, are listed in

  Table 1. MR and TR relevant parameters.

  Parameter	Symbol	Computer Mnemonic	Value	Units
  Main Rotor
  MR radius	$R_{MR}$	ROTOR	3.8	m
  MR chord	$c_{MR}$	CHORD	0.195	m
  MR rotational speed	${\Omega }_{MR}$	OMEGA	528.5	rpm
  MR Lock number	${\gamma }_{MR}$	GAMMA	4.25	-
  MR hinge offset	$\epsilon$	EPSLN	0	percent/100
  MR flapping spring constant	$K_{\beta }$	AKBETA	0	N m/rad
  MR tangent of ${\delta }_3$	$K_1$	AKONE	0	-
  MR solidity	${\sigma }_{MR}$	SIGMA	0.0327	-
  MR hub stationline	${\mathit{STA}}_H$	STAH	2	m
  MR hub buttline	$B{L}_H$	BLH	0	m
  MR hub waterline	$W{L}_H$	WLH	2.4	m
  Tail Rotor
  TR radius	$R_{TR}$	RTR	0.57	m
  TR chord	$c_{TR}$	cTR	0.12	m
  TR rotational speed	${\Omega }_{TR}$	OMTR	3 061.8	rpm
  TR tangent of ${\delta }_3$	${K_1}_{TR}$	FKITR	1	-
  TR solidity	${\sigma }_{TR}$	STR	0.0382	-
  MR hub stationline	${\mathit{STA}}_{\mathit{TR}}$	STATR	6.4	m
  MR hub buttline	$B{L}_{TR}$	BLTR	−0.25	m
  MR hub waterline	$W{L}_{TR}$	WLTR	1.34	m

  and

  Table 2. Fuselage, empennages, and miscellaneous components location.

  Parameter	Symbol	Computer Mnemonic	Value	Units
  Fuselage (Fus.)
  Fus. aerodynamic ref. point stationline	${\mathit{STA}}_{RPF}$	STARPF	0	m
  Fus. aerodynamic ref. point buttline	$B{L}_{RPF}$	BLRPF	0	m
  Fus. aerodynamic ref. point waterline	$W{L}_{RPF}$	STARPF	0	m
  Horizontal stabilizer (HS)
  HS stationline	${\mathit{STA}}_{HS}$	STAHS	6.199	m
  HS buttline	$B{L}_{HS}$	BLHS	0.435	m
  HS waterline	$W{L}_{HS}$	WLHS	1.394	m
  Upper vertical fin (VF1)
  VF1 stationline	${\mathit{STA}}_{VF1}$	STAVF1	6.1	m
  VF1 buttline	$B{L}_{VF1}$	BLVF1	0.052	m
  VF1 waterline	$W{L}_{VF1}$	WLVF1	1.683	m
  Lower vertical fin (VF2)
  VF2 stationline	${\mathit{STA}}_{VF2}$	STAVF2	6.069	m
  VF2 buttline	$B{L}_{VF2}$	BLVF2	0.048	m
  VF2 waterline	$W{L}_{VF2}$	WLVF2	0.996	m
  Main rotor hub (MRH)
  MRH stationline	${\mathit{STA}}_{MRH}$	STAMRH	2	m
  MRH buttline	$B{L}_{MRH}$	BLMRH	0	m
  MRH waterline	$W{L}_{MRH}$	WLMRH	2.4	m
  Parachute canopy (PC)
  PC stationline	${\mathit{STA}}_{PC}$	STAPC	2	m
  PC buttline	$B{L}_{PC}$	BLPC	0	m
  PC waterline	$W{L}_{PC}$	WLPC	2.468	m

  , together with relevant helicopter data.

Let $s(·)=sin(·)$, $c(·)=cos(·)$. Vector transformation between ${\mathcal{F}}_H$ and ${\mathcal{F}}_B$ is provided by the rotation matrix [12]

$$\mathit{R}\left(\mathbf{\alpha }\right)=\left[\begin{array}{ccc}c\theta c\psi & c\theta s\psi & -s\theta \\ s\varphi s\theta c\psi -c\varphi s\psi & s\varphi s\theta s\psi +c\varphi c\psi & s\varphi c\theta \\ c\varphi s\theta c\psi +s\varphi s\psi & c\varphi s\theta s\psi -s\varphi c\psi & c\varphi c\theta \end{array}\right]$$

(1)

obtained by a 3-2-1 Euler rotation sequence where $\mathbf{\alpha }={[\varphi ,\theta ,\psi ]}^T$ describes the attitude of the rotorcraft in terms of the classical ‘roll’, ‘pitch’, and ‘yaw’ angles, respectively. The following notation is adopted: if w is an arbitrary vector, its components are transformed from ${\mathcal{F}}_H$ to ${\mathcal{F}}_B$ through ${\mathit{w}}_B=\mathit{R}{\mathit{w}}_H$. In what follows, the subscript B will be dropped for simplicity.

3.2. Rigid Body Dynamics

Vehicle dynamics is described by Newton–Euler equations of motion projected in ${\mathbb{F}}_B$, namely:

$$\dot{\mathbf{v}}=-\mathbf{\omega }\times \mathbf{v}+\mathbf{F}/m$$

(2)

$$\dot{\mathbf{\omega }}={\mathit{J}}^{-1}\left[-\mathbf{\omega }\times \left(\mathit{J}\mathbf{\omega }\right)+\mathbf{M}\right]$$

(3)

where $\mathit{v}={[u,v,w]}^T$ is linear velocity, $\mathbf{\omega }={[p,q,r]}^T$ is angular velocity,

$$\mathit{J}=\left[\begin{array}{ccc}{J}_{xx}& -J_{xy}& -J_{xz}\\ -J_{xy}& J_{yy}& -J_{yz}\\ -J_{xz}& -J_{yz}& J_{zz}\end{array}\right]$$

(4)

is the inertia tensor about $CG$ with respect to ${\mathcal{F}}_B$, and m is the total mass of the rotorcraft. $\mathbf{F}={[F_x,F_y,F_z]}^T$ and $\mathbf{M}={[M_x,M_y,M_z]}^T$ are the external force and moment vectors, respectively.

The external force acting on the rotorcraft is made of gravity, ${\mathit{F}}^{\left(g\right)}$, and aerodynamic, ${\mathit{F}}^{\left(a\right)}$, contributions. Taking into account Equation (1), gravity force vector expressed in the body frame is

$${\mathit{F}}^{\left(g\right)}=\mathit{R}\left(\mathbf{\alpha }\right)\left[\begin{array}{c}0\\ 0\\ mg\end{array}\right]=mg\left[\begin{array}{c}-sin\theta \\ sin\varphi cos\theta \\ cos\varphi cos\theta \end{array}\right]$$

(5)

where g is gravitational acceleration, described by means of WGS84 Taylor series model [26].

Rotorcraft attitude kinematics, that relates the generalized velocity $\dot{\mathbf{\alpha }}$ and the angular velocity $\mathbf{\omega }$ is given by [12]:

$$\dot{\mathbf{\alpha }}=\left[\begin{array}{ccc}1& sin\varphi tan\theta & cos\varphi tan\theta \\ 0& cos\varphi & -sin\varphi \\ 0& sin\varphi /cos\theta & cos\varphi /cos\theta \end{array}\right]\mathbf{\omega }$$

(6)

while the position of the helicopter ${\mathit{p}}_E={[x_E,y_E,z_E]}^T$, with components expressed in the inertial frame ${\mathcal{F}}_E$, is obtained from the equation:

$${\dot{\mathit{p}}}_E=\mathit{R}{\left(\mathbf{\alpha }\right)}^T\mathit{v}$$

(7)

3.3. Aerodynamic Forces and Moments

The characterization of aerodynamic force, ${\mathit{F}}^{\left(a\right)}={[X,Y,Z]}^T$, and moment, ${\mathit{M}}^{\left(a\right)}={[L,M,N]}^T$, is performed on the basis of the model detailed in [25], whose nomenclature is adopted in the present work. A conventional single MR helicopter with teetering configuration and counterclockwise rotation are considered. Contributions are provided by the main rotor (MR), tail rotor (TR), fuselage (F), horizontal stabilizer (HS), upper and lower vertical fins (VF1 and VF2), main rotor hub (MRH), and parachute canopy (PC). Air parameters are calculated from the International Standard Atmosphere (ISA) model as a function of rotorcraft altitude [27].

3.3.1. MR and TR Modeling

The following assumptions and simplifications are made about the MR model: (a) rotor blades are rigid in bending and torsion; (b) flapping angles are small, and the analysis follows the simple strip theory [28]; (c) the effects of aircraft motion on blade flapping are limited to those related to the angular accelerations $\dot{p}$ and $\dot{q}$, the angular rates p and q, and the normal acceleration component $\dot{w}$; (d) blade flow stall is disregarded; (e) rotor inflow is uniform, and no inflow dynamics is modeled; (f) main rotor blade flapping is approximated by the first harmonic terms with time-varying coefficients, that is

$$\beta \left(t\right)=a_0-a_{1}cos\xi -b_{1}sin\xi$$

(8)

where $a_0$ is treated as a preset constant (coning angle) and $\xi$ is blade azimuth. Coefficients $a_1\left(t\right)$ and $b_1\left(t\right)$ respectively represent the longitudinal and lateral tilt of the rotor tip–path plane, obtained as solutions to the equations in Appendix C of [25] with null hinge offset ratio, $ϵ=0$, flapping spring constant, $K_{\beta }=0$, and pitch-flap coupling ratio, $tan{\delta }_3=0$. Finally, the MR shaft is aligned with ${\mathit{z}}_B$.

The tail rotor is modeled according to a teetering configuration without cyclic pitch. Provided that the flapping frequency is typically much higher than that of the MR system, TR tip-path plane dynamics are neglected, no flapping spring constant is considered, and the pitch-flap coupling ratio, ${{\delta }_3}_{TR}$, is characterized by a non-null value (see Appendix D in [25]).

Contrary to some of the assumptions provided in [25], the blades of both MR and TR are characterized by cambered airfoils with a lift-curve slope $a<2\pi$ 1/rad and a zero-lift angle of attack ${\alpha }_0\ne 0$. The rotor blade profile drag coefficient, $C_d$, is calculated as

$$C_d=0.008+0.3{\left({\displaystyle \frac{6C_T}{\sigma a}}\right)}^2+\Delta C_d$$

(9)

where $C_T$ is the rotor thrust coefficient, $\sigma$ is rotor solidity, and $\Delta C_d$ is the extra drag coefficient determined by flow compressibility effects. Let $M_{90}$ be the Mach number evaluated at the tip of the advancing blade, where $\xi =90$ deg. In order to estimate the extra drag, the approximate model proposed by Prouty and described in [28] is adopted, where

$$\Delta C_d\left(M_{90}\right)=\left\{\begin{array}{cc}12.5{\left(M_{90}-M_{dr}\right)}^3\hfill & \mathrm{for}{M}_{90}\ge M_{dr}\hfill \\ 0\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(10)

and $M_{dr}=0.74$ is the drag-rise Mach number. With respect to the characterization of rotor inflow, a number of non-ideal effects are considered, based on the approach in [28], for the characterization of forces and moments. A constant tip-loss factor $B<1$ is adopted to account for blade tip losses. Other non-ideal effects, including nonuniform inflow, wake swirl and contraction, and blade interference, are accounted for by an induced power factor $k_i$, assumed to be a constant. The MR in-ground effect is provided by the model in ref. [29], and the inflow iterative scheme is solved according to Halley’s method with a damping coefficient equal to $0.01$ [30].

Cockpit/RC control of MR is provided by pilot commands in terms of lateral cyclic ${\delta }_a$, longitudinal cyclic ${\delta }_e$, and collective ${\delta }_c$. All commands are expressed in terms of non-dimensional variables, such that ${\delta }_a\in [-1,+1]$ (positive direction: right to generate $L>0$), ${\delta }_e\in [-1,+1]$ (positive direction: aft to generate $M>0$), and ${\delta }_c\in [-1,+1]$ (positive direction: up to generate $Z<0$). Onboard control of the tail rotor is performed by pedal commands, expressed as ${\delta }_p\in [-1,+1]$ (positive direction: right pedal forward to generate $N>0$). The transformation of pilot commands into blade pitch angles is provided by a set of low-order polynomial functions, ${A_1}_s=C_1\left({\delta }_a\right)$, ${B_1}_s=C_2\left({\delta }_e\right)$, ${\theta }_0=C_3\left({\delta }_c\right)$, and ${{\theta }_0}_{TR}=C_4\left({\delta }_p\right)$, provided by the manufacturer. ${A_1}_s$ and ${B_1}_s$, respectively, represent the lateral and longitudinal cyclic pitch angles measured from the MR hub plane in ${\mathcal{F}}_B$. Rotor blades are modeled with a linear twist, such that ${\theta }_0$ is the blade collective pitch ideally extrapolated to the rotor center and ${\theta }_{tw}$ is the total blade twist angle (tip minus root pitch angle). No twist characterizes TR blades, where collective pitch is identified by ${{\theta }_0}_{TR}$.

An additional degree of freedom is related to the power plant made of free turbines, MR, and TR transmissions. In particular, MR and TR rotational speeds vary according to the current torque requirements and the engine power available. Changes in speed cause the free turbine governor to vary fuel flow to change the available power and maintain the desired angular rate. The engine dynamic model is found in [25]. For the sake of brevity, details are not provided in the present paper. Modeling parameters in terms of maximum available power, engine dynamics, specific fuel consumption, and mechanical transmission efficiency are provided by the manufacturer.

3.3.2. Fuselage, Empennages, and Miscellaneous Components

With respect to fuselage aerodynamics, it is assumed that longitudinal forces and moments are dependent on fuselage angle of attack and lateral forces and moments are dependent on angle of sideslip. The exception is the drag force, which is assumed to have a contribution from both angles of attack and sideslip. The modeling is based on a low- and a high-angle representation of forces and moments, according to Appendix F in [25], with data obtained through a detailed computational fluid dynamics characterization. Phasing between the two approximations is performed by means of cubic spline interpolation, with improved performance with respect to the proposed linear transition.

The modeling of the two vertical empennages and of the horizontal stabilizer also follows the approach in [25]. The aerodynamics of the MR hub and parachute pod are assessed by the equivalent flat plate area model. As an example, the force vector generated by MRH is expressed as:

$${\mathit{F}}_{MRH}=-{\displaystyle \frac{1}{2}}\rho \left(\left[{A_x}_{MRH},{A_y}_{MRH},{A_z}_{MRH}\right]{\mathit{V}}_{MRH}\right){\mathit{V}}_{MRH}$$

(11)

where ${\mathit{V}}_{MRH}={[u_{MRH},v_{MRH},w_{MRH}]}^T$ is the velocity, relative to the air mass, of the main rotor hub and includes the contribution of MR downwash, according to [31]. ${A_x}_{MRH}$, ${A_y}_{MRH}$, and ${A_z}_{MRH}$ are the equivalent flat plate drag areas, respectively orthogonal to ${\mathit{x}}_B$, ${\mathit{y}}_B$, and ${\mathit{z}}_B$. The moment generated by ${\mathit{F}}_{MRH}$ about $CG$ is given by ${\mathit{M}}_{MRH}={\mathit{d}}_{MRH}\times {\mathit{F}}_{MRH}$, where

$${\mathit{d}}_{MRH}=\left[\begin{array}{c}ST{A}_{CG}-ST{A}_{MRH}\\ B{L}_{MRH}-B{L}_{CG}\\ W{L}_{CG}-W{L}_{MRH}\end{array}\right]$$

(12)

is the vector directed from $CG$ to MRH position, assumed to be coincident with its center of pressure, with constant components expressed in ${\mathcal{F}}_B$.

4. Trim and Stability Analysis

The nonlinear model described in Section 3 is implemented in the Matlab/Simulink environment, where differential equations are solved by the Dormand-Prince ode8 method with a frequency of 1000 Hz [32]. In what follows, (1) the trim conditions are determined for different cruise speeds, (2) a linearization procedure is applied to the complete model about such equilibria, and (3) an open-loop dynamic analysis is performed to investigate the helicopter control and stability properties.

4.1. Trim Analysis

The helicopter model is numerically trimmed for straight-and-level flight at $h=50$ m in standard atmospheric conditions. Different values of forward speed are considered, ranging from 0 km/h (hover) to 180 km/h (approximately the never-exceed speed), with steps of 5 km/h. For the sake of brevity, the results of both the static and the following dynamic analysis are summarized only for the hovering condition, for which dedicated flight tests were performed for validation purposes.

The main results of trim analysis for the hovering condition are given in

Table 3. Trim analysis for the hovering flight.

Parameter	Symbol	Value	Units	Est. Error $\left|·\right|$
Main Rotor
Long. first–harmonic flapping coeff.	$a_1$	2.92	deg	N/A
Lat. first–harmonic flapping coeff.	$b_1$	−1.10	deg	N/A
Induced speed	$v_i$	7.79	m/s	N/A
Aerodynamic torque	Q	1 579.5	Nm	3.1%
Tail Rotor
Long. first–harmonic flapping coeff.	${a_1}_{TR}$	0.24	deg	N/A
Lat. first–harmonic flapping coeff.	${b_1}_{TR}$	−0.24	deg	N/A
Induced speed	${v_i}_{TR}$	11.64	m/s	N/A
Aerodynamic torque	$Q_{TR}$	26.1	Nm	4.4%
Fuselage
Roll angle	$\varphi$	−2	deg	17.6%
Pitch angle	$\theta$	−1.96	deg	6.7%
Control Pitch Angles
MR lat. cyclic pitch	${A_1}_s$	−1.10	deg	4.8%
MR lon. cyclic pitch	${B_1}_s$	−2.91	deg	42.0%
MR collective pitch	${\theta }_0$	12.62	deg	2.1%
TR collective pitch	${{\theta }_0}_{TR}$	8.28	deg	2.1%

and compared with the data available from flight tests performed with the same vehicle configuration (deviations with respect to measured data are reported in terms of absolute values of percentage errors). To this end, the helicopter was equipped with a set of sensors, including: (a) potentiometers for cockpit command acquisition and blade pitch measurement, (b) torque–meters for MR and TR torque analysis, and (c) a AHRS providing rigid body attitude, angular rate, acceleration, speed, and position information.

According to Table 3, good agreement is found between predicted and measured values, showing the validity of the modeling approach. A major difference characterizes the longitudinal cyclic pitch, with a $42\%$ error. It must be noted that a degree of uncertainty characterizes the knowledge of $CG$ position (especially the $ST{A}_{CG}$ parameter) in the actual flight configuration, which is estimated by means of CAD analysis and suspension techniques. Uncertainty also characterizes the aerodynamics of the fuselage, especially in the case of hovering and low-speed forward flight, where MR wake envelops a large portion of the fuselage. For the aim of the present analysis, the model adopted for both MR inflow and fuselage aerodynamics necessarily represents a compromise solution, which allows for satisfactory accuracy in terms of the dynamic characterization of rotorcraft without the cost of excessively-complex aerodynamic models.

The match expected at hover between MR cyclic pitch angles and flapping coefficients, namely ${A_1}_s=b_1$ and ${B_1}_s=-a_1$, holds almost exactly in Table 3. Slight differences occur for the simulated hover condition, which is actually obtained by flying the helicopter at a residual forward speed of $0.1$ m/s. With respect to the experimental campaign, effective environmental conditions were also monitored, provided the helicopter was maintained in upwind hover while estimating a maximum wind speed of 20 km/h.

4.2. Dynamic Analysis

Consider the equations of motion introduced in Section 3 and detailed in [25]. In nonlinear form, it is

$$\dot{\mathit{x}}=\mathit{f}\left(\mathit{x},\mathit{u},t\right)$$

(13)

provided x is rigid–body state vector, namely

$$\mathit{x}={[u,w,q,\theta ,v,p,\varphi ,r]}^T$$

(14)

while time evolution of $x_E$, $y_E$, $z_E$, and $\psi$ is not accounted in the framework of system linearization. Control vector u has four components, expressed in terms of pilot commands as:

$$\mathit{u}={[{\delta }_c,{\delta }_e,{\delta }_a,{\delta }_p]}^T$$

(15)

Using small perturbation theory [29], helicopter motion is described in terms of perturbation from the equilibrium condition, ${\mathit{x}}_e={[U_e,W_e,Q_e,{\Theta }_e,V_e,P_e,{\Phi }_e,R_e]}^T$ and ${\mathit{u}}_e={[{U_1}_e,{U_2}_e,{U_3}_e,{U_4}_e,]}^T$, written in the form $\mathit{x}={\mathit{x}}_e+\delta \mathit{x}$ and $\mathit{u}={\mathit{u}}_e+\delta \mathit{u}$. By following the approach and the nomenclature of [29], given the trim conditions in Section 4.1, the model in Equation (13) is linearized at all considered speeds to respectively obtain system and input matrices

$$\mathit{A}={\left({\displaystyle \frac{\partial \mathit{f}}{\partial \mathit{x}}}\right)}_{{\mathit{x}}_e,{\mathit{u}}_e},\mathit{B}={\left({\displaystyle \frac{\partial \mathit{f}}{\partial \mathit{u}}}\right)}_{{\mathit{x}}_e,{\mathit{u}}_e}$$

(16)

as a function of aerodynamic derivatives. The latter are estimated by numerical differencing in the Matlab/Simulink environment [29]. To this end, aerodynamic forces and moments are positively perturbed by each of the state and input vector components in turn, with amplitude equal to $0.02$ (respectively intended in terms of m/s for u, v, w, rad/s for p, q, r, rad for $\varphi$, $\theta$, and non-dimensional units for control inputs). State and control derivatives are written in the form:

$$X_u={\displaystyle \frac{1}{m}}{\displaystyle \frac{\partial X}{\partial u}}$$

(17)

and

$$L_p^{\prime }={\displaystyle \frac{J_{zz}}{J_{xx}{J}_{zz}-J_{xz}^2}}{\displaystyle \frac{\partial L}{\partial p}}+{\displaystyle \frac{J_{xz}}{J_{xx}{J}_{zz}-J_{xz}^2}}{\displaystyle \frac{\partial N}{\partial p}}$$

(18)

$$N_r^{\prime }={\displaystyle \frac{J_{xz}}{J_{xx}{J}_{zz}-J_{xz}^2}}{\displaystyle \frac{\partial L}{\partial r}}+{\displaystyle \frac{J_{xx}}{J_{xx}{J}_{zz}-J_{xz}^2}}{\displaystyle \frac{\partial N}{\partial r}}$$

(19)

A total of 36 stability derivatives and 24 control derivatives are determined in the standard 6DOF representation for each flight condition. Due to the highly-classified nature of the data involved in the project, only one sample derivative is analyzed in the present paper at hover. A qualitative discussion about the behavior of the most significant derivatives is provided in what follows.

The effect of linear velocity on aerodynamic forces is principally taken into account by $X_u$, $Y_v$, and $Z_w$. The force damping derivatives $X_u<0$ and $Y_v<$ 0, which respectively reflect the drag and side force on rotor–fuselage combination, steadily increase in absolute value and are practically linear with speed beyond 50 km/h. At low speed, the effect of disc tilt following perturbations in u and v becomes predominant. Similar considerations hold for the heave damping derivative $Z_w$, which is mostly influenced by the fuselage and horizontal empennage in high-speed flight. At low speed, the MR tends to dominate $Z_w$ through a reduction in $C_T$ determined by a vertical speed perturbation. In order to validate the numerical linearization routine, a comparison is performed with the analytical results obtained for the stability and control derivatives according to formulas available in the literature. As an example, the MR contribution only to $Z_w$ can be analytically estimated as [29]

$$Z_w=-{\displaystyle \frac{\rho \left(\Omega R\right)\pi R^2}{m}}{\displaystyle \frac{\partial C_T}{\partial {\mu }_z}}$$

(20)

where ${\mu }_z=w/\left(\Omega R\right)$ is MR climb ratio and

$${\displaystyle \frac{\partial C_T}{\partial {\mu }_z}}\approx {\displaystyle \frac{2a\sigma \left|\lambda \right|}{16\left|\lambda \right|+a\sigma }}$$

(21)

Based on the data in Table 1 and Table 3, it is $\lambda =-0.0371$ at hover, such that $\partial C_T/\partial {\mu }_z\approx 0.018$. It follows $Z_w\approx -0.317$ 1/s, which is close to the value numerically obtained in the same condition for the full helicopter, namely $-0.345$ 1/s. In such a case, the estimation error obtained according to literature results is $-8.1\%$, provided that fuselage and appendages contributions are disregarded.

The speed stability effect is observed in $M_u>0$ and $L_v^{\prime }<0$, the latter showing a practically linear behavior with speed. $M_w>0$ is representative of the incidence static stability effect, which increases non-monotonically with speed and approximately tracks $M_u$, being influenced by MR inflow on helicopter components. Finally, $N_v^{\prime }>0$ accounts for the weathercock effect by means of TR and vertical fins (stabilizing with speed) and the fuselage (destabilizing).

The damping derivatives $L_p^{\prime }<0$, $M_q<0$, and $N_r^{\prime }<0$ reflect short-term, small, and moderate-amplitude handling characteristics. If, on the one hand, $L_p^{\prime }$ and $M_q$ principally account for MR flapping motion in the presence of roll and pitch rate perturbations, $N_r^{\prime }$ is dominated by loads on TR and vertical fins, with a stronger yaw-damping effect at high forward speeds.

Given the stability and control derivatives obtained above, the complete system and input matrices A and B are generated according to the structure provided on page 277 in [29]. Note that, with the idea of designing closed-loop control systems, the input matrix B is configured for application to the non-dimensional pilot commands. The formulation in terms of blade pitch control angles is however possible by means of the mapping functions $C_1$, $C_2$, $C_3$, and $C_4$ introduced in Section 3.3.1.

For the aim of the present work, however, the results of a decoupled analysis are first discussed. Based on the approximate separation between the longitudinal and the lateral-directional dynamics, the decoupled representation is available in ref. [29], where input matrices are applied to blade pitch control angles. The longitudinal dynamics are described by the forced system:

$$\begin{array}{ll}{\displaystyle \frac{d}{dt}}\left[\begin{array}{c}u\\ w\\ q\\ \theta \end{array}\right]=& \underset{{\mathit{A}}_{lon}}{\underbrace{\left[\begin{array}{cccc}{X}_u& X_w& X_q-W_e& -gcos{\Theta }_e\\ Z_u& Z_w& Z_q+U_e& -gsin{\Theta }_e\\ M_u& M_w& M_q& 0\\ 0& 0& 1& 0\end{array}\right]}}\left[\begin{array}{c}u\\ w\\ q\\ \theta \end{array}\right]\\ & +\underset{{\mathit{B}}_{lon}}{\underbrace{\left[\begin{array}{cc}{X}_{{\theta }_0}& X_{{B_1}_s}\\ Z_{{\theta }_0}& Z_{{B_1}_s}\\ M_{{\theta }_0}& M_{{B_1}_s}\\ 0& 0\end{array}\right]}}\left[\begin{array}{c}{\theta }_0\\ {B_1}_s\end{array}\right]\end{array}$$

(22)

A pair of complex-conjugate poles is determined from ${\mathit{A}}_{lon}$, which is related to an unstable phugoid mode with natural frequency ${\omega }_{ph}$ and time constant ${\tau }_{ph}$ (calculated as the reciprocal of the real part of the poles in its absolute value). Two real stable modes are also evaluated, namely heave and pitch subsidence effects. The first pole, identified by $p_{hv}<0$, is practically determined by the vertical damping derivative $Z_w$. The second pole, $p_{ps}<0$, accounts for the fundamental contribution of both $Z_w$ and $M_q$ and is characterized by a time constant approximately estimated as ${\tau }_{ps}\approx -1/\left(Z_w+M_q\right)$ [29]. The decoupled lateral-directional dynamics are defined by the system:

$$\begin{array}{ll}{\displaystyle \frac{d}{dt}}\left[\begin{array}{c}v\\ p\\ r\\ \varphi \end{array}\right]=& \underset{{\mathit{A}}_{lat}}{\underbrace{\left[\begin{array}{cccc}{Y}_v& Y_p+W_e& Y_r-U_e& gcos{\Phi }_{e}cos{\Theta }_e\\ L_v^{\prime }& L_p^{\prime }& L_r^{\prime }& 0\\ N_v^{\prime }& N_p^{\prime }& N_r^{\prime }& 0\\ 0& 1& cos{\Phi }_{e}tan{\Theta }_e& 0\end{array}\right]}}\left[\begin{array}{c}v\\ p\\ r\\ \varphi \end{array}\right]\\ & +\underset{{\mathit{B}}_{lat}}{\underbrace{\left[\begin{array}{cc}{Y}_{{A_1}_s}& Y_{{{\theta }_0}_{TR}}\\ L_{{A_1}_s}^{\prime }& L_{{{\theta }_0}_{TR}}^{\prime }\\ N_{{A_1}_s}^{\prime }& N_{{{\theta }_0}_{TR}}^{\prime }\\ 0& 0\end{array}\right]}}\left[\begin{array}{c}{A_1}_s\\ {{\theta }_0}_{TR}\end{array}\right]\end{array}$$

(23)

A pair of complex-conjugate poles is derived from ${\mathit{A}}_{lat}$ with the natural frequency ${\omega }_{dr}$. Such poles characterize the dutch-roll mode, which is unstable but slowly develops with a time constant ${\tau }_{dr}$. The roll subsidence mode, mostly determined by the damping derivative $L_p^{\prime }$, is related to the real pole $p_{roll}<0$. The spiral subsidence mode at hover is stable, $p_{spiral}<0$, and dampens with a time constant ${\tau }_{spiral}$.

The analysis of coupled representation behind state matrix A is also considered, and the obtained poles are marked by a superscript ‘c’. A comparison with the corresponding values derived through the decoupled analysis is provided where possible. Two real poles are first extracted. The roll subsidence effect is recognized in the first pole, $p_1=0.85·p_{roll}$, provided $p_1\approx L_p^{\prime }$. The same consideration holds for vertical damping mode, identified by $p_8^{\left(c\right)}\approx p_{hv}\approx Z_w$. Three pairs of complex-conjugate poles complete the analysis. The first pair, $p_{2,3}^{\left(c\right)}$, is stable with a real part proportional to $M_q+Z_w$. It is representative of a damped oscillation with a natural frequency ${\omega }_{py}^{\left(c\right)}$ and time constant ${\tau }_{py}^{\left(c\right)}$, determined by the coupling of pitch and yaw subsidence modes. The second pair, $p_{4,5}^{\left(c\right)}$, characterizes the unstable phugoid mode, which develops with a time constant ${\tau }_{ph}^{\left(c\right)}=0.75·{\tau }_{ph}$ and shows natural frequency ${\omega }_{ph}^{\left(c\right)}=1.38·{\omega }_{ph}$. The last pair of complex poles, $p_{6,7}$ characterizes the dutch roll motion, which is unstable and develops with natural frequency ${\omega }_{dr}^{\left(c\right)}=0.92·{\omega }_{dr}$ and time constant ${\tau }_{dr}^{\left(c\right)}=1.07·{\tau }_{dr}$.

4.3. Model Validation

In Section 4, a comparison is provided between simulated and measured variables regarding the static characterization of hovering conditions. In what follows, predicted dynamic properties about the same equilibrium are validated through identification methods [33]. To this end, flight data are collected and eventually filtered after performing frequency sweep maneuvers about the hover, according to the approach described in [34]. The frequency response for each selected input–output pair is then identified during an optimization process driven by the difference between the computed and the predicted frequency responses. The fidelity of the model is finally established using time domain verification, according to which time response predicted by the identified model is compared with the response recorded during flight tests.

Different maneuvers and data pairs are considered for the identification of transfer functions, such as $p\left(s\right)/{A_1}_s\left(s\right)$ and $q\left(s\right)/{B_1}_s\left(s\right)$, with the aim of validating the predicted dynamic information. For the sake of brevity, the adopted procedure is detailed for the characterization of the heave subsidence mode, whose dominant derivative $Z_w$ is discussed above. In particular, the first input-output data pair describes the effect of MR collective pitch angle ${\theta }_0$ on vertical acceleration, $a_z=\dot{w}$, expressed in a body-fixed frame. A detail of the data taken into account for such an identification procedure is reported in Figure 5. The predicted transfer function as obtained from the state-space representation in Equation (22) is:

$${\left.{\displaystyle \frac{a_z\left(s\right)}{{\theta }_0\left(s\right)}}\right|}_{mdl}={\displaystyle \frac{Z_{{\theta }_0}s\left(s-z_1\right)\left(s-z_2\right)\left(s-z_3\right)}{\left(s-p_{hv}\right)\left(s-p_{ps}\right)\left(s^2-2/{\tau }_{ph}s+{\omega }_{ph}^2\right)}}$$

(24)

where a set of 4 zeros is determined. The first one is located at the origin, $z_1=0.9983p_{ps}$ is real negative, and $z_2$, $z_3$ are a complex–conjugate pair such that $z_2{z}_3=1.008{\omega }_{ph}^2$ and $z_2+z_3=1.059·2/{\tau }_{ph}$. Heave subsidence mode evidently dominates the motion along ${\mathit{z}}_B$, provided that almost perfect pole-zero cancellation characterizes the terms depicted in gray color. It follows:

$${\left.{\displaystyle \frac{a_z\left(s\right)}{{\theta }_0\left(s\right)}}\right|}_{mdl}\approx {\displaystyle \frac{Z_{{\theta }_0}s}{s-p_{hv}}}$$

(25)

Numerical identification is performed by using a Prediction Error Minimization (PEM) method focused on simulation [35], provided the transfer function in Equation (24) is assumed as the initial guess model. The identified transfer function is:

$$\begin{array}{ll}{\left.{\displaystyle \frac{a_z\left(s\right)}{{\theta }_0\left(s\right)}}\right|}_{id}& ={\displaystyle \frac{0.9964Z_{{\theta }_0}\left(s-{\pi }_1\right)\left(s-{\pi }_2\right)}{\left(s-1.0421p_{hv}\right)\left(s-0.9912{\pi }_2\right)}}\\ & ·{\displaystyle \frac{\left(s^2+{\pi }_{3}s+{\pi }_4\right)}{\left(s^2+1.0090{\pi }_{3}s+0.9964{\pi }_4\right)}}\\ & \approx {\displaystyle \frac{0.9964Z_{{\theta }_0}s}{\left(s-1.0421p_{hv}\right)}}\end{array}$$

(26)

where pole–zero cancellation can evidently be performed for the gray terms. It must be noted that ${\pi }_1\approx 0$, such that the zero at the origin is also recovered. The estimation error between model-predicted and identified parameters is provided in the second line of Equation (26), where the updated values of $\left|Z_{{\theta }_0}\right|$ and $\left|p_{hv}\right|$, respectively, result in being $0.36\%$ smaller and $4.21\%$ bigger than the model-predicted ones in Equation (25). A sample comparison between measured and refined-simulation data after heave subsidence mode characterization is finally provided in Figure 6 for the acceleration.

Encouraging results are indeed obtained for other input-output pairs, thus validating the modeling approach. In all cases, in fact, very good agreement is found between the dynamic properties obtained through numerical simulation and identification techniques.

5. Control System Design and Test

In what follows, the control system design phase is described based on the mathematical model in Section 3 and the analysis performed in Section 4. The closed-loop system is first analyzed by Model-In-the-Loop (MIL) simulations, where linear controllers are directly designed and validated in the nonlinear framework by means of an extensive campaign of simulations performed in collaboration with the candidate pilot. HIL tests are then performed to refine the control gains and validate the software/hardware setup [36].

5.1. Model–in–the–Loop Validation

Pilot commands, here named ${\delta }_a^{\left(pilot\right)}$, ${\delta }_e^{\left(pilot\right)}$, ${\delta }_c^{\left(pilot\right)}$, and ${\delta }_p^{\left(pilot\right)}$ are an input to the control system and follow the same convention described in Section 3.3.1. According to the given requirements, no closed-loop control is designed for MR collective pitch, such that ${\delta }_c\equiv {\delta }_c^{\left(pilot\right)}$.

In the framework of control system design, simulation models for selected AHRS and actuators are also developed. Modeling parameters in terms of accuracy and performance are obtained from both datasheets and dedicate experiments performed in laboratory facilities.

The first controller is designed to stabilize yaw motion by the actuation of TR collective pitch angle ${{\theta }_0}_{TR}$ through the closed-loop feedback of yaw rate r. Let $e_r={\xi }_r{\delta }_p^{\left(pilot\right)}-r$ be the error between the desired and the measured angular rate, provided that ${\xi }_r>0$ is a prescribed constant that transforms the non-dimensional command provided by the remote pilot into the desired yaw rate. The control scheme is described by the equation:

$${\delta }_p=k_p^{\left(r\right)}{e}_r+k_i^{\left(r\right)}{\int }_0^t{e}_r\left(s\right)ds$$

(27)

where $k_p^{\left(r\right)}>0$ and $k_i^{\left(r\right)}>0$ are control gains, respectively, providing proportional and integral contributions related to the error signal $e_r\left(t\right)$.

The second controller is used to stabilize the fuselage’s attitude in terms of roll and pitch angles by the actuation of MR lateral and longitudinal cyclic control angles, respectively. With respect to roll angle stabilization, it is:

$${\delta }_a=k_p^{\left(\varphi \right)}{e}_{\varphi }+k_i^{\left(\varphi \right)}{\int }_0^t{e}_{\varphi }\left(s\right)ds+k_d^{\left(\varphi \right)}p$$

(28)

where $k_p^{\left(\varphi \right)}>0$ and $k_i^{\left(\varphi \right)}>0$. A derivative–like contribution is also provided by the direct feedback of roll rate p through the gain $k_d^{\left(\varphi \right)}<0$. The error between desired and measured roll angle is calculated as $e_{\varphi }={\xi }_{\varphi }{\delta }_a^{\left(pilot\right)}-\varphi$, where ${\xi }_{\varphi }>0$ is a prescribed constant. Controller structure for the stabilization of pitch angle follows the same approach, namely:

$${\delta }_e=k_p^{\left(\theta \right)}{e}_{\theta }+k_i^{\left(\theta \right)}{\int }_0^t{e}_{\theta }\left(s\right)ds+k_d^{\left(\theta \right)}q$$

(29)

where $k_p^{\left(\theta \right)}>0$, $k_i^{\left(\theta \right)}>0$, and $k_d^{\left(\theta \right)}<0$. The error between desired and measured pitch angle is $e_{\theta }={\xi }_{\theta }{\delta }_e^{\left(pilot\right)}-\theta$, where ${\xi }_{\theta }>0$.

In Figure 7 and Figure 8 the results of a sample maneuver are reported. Simulation is started at $h=50$ m with null attitude of the helicopter (${\varphi }_0={\theta }_0={\psi }_0=0$ deg) and an initial angular rate about the yaw axis, such that $p_0=q_0=0$ deg/s and $r_0=-5$ deg/s. MR collective pitch angle is kept constant and equal to the value obtained in Table 3 for the hovering condition, namely ${\theta }_0=12.62$ deg, corresponding to ${\delta }_c^{\left(pilot\right)}=0.495$. Let ${\xi }_r=40·\pi /180$ rad/s, ${\xi }_{\varphi }=25·\pi /180$ rad, and ${\xi }_{\theta }=12·\pi /180$ rad. Input values to the controllers are ${\delta }_p^{\left(pilot\right)}=0$, ${\delta }_a^{\left(pilot\right)}=-0.08$, and ${\delta }_e^{\left(pilot\right)}=-0.163$, which respectively provide the desired values $\varphi =-2$ deg, $\theta =-1.96$, and $r=0$ deg/s necessary to hover. In Figure 7, state variables describing fuselage attitude are plotted as a function of time, showing the stabilizing effect of implemented controllers. The corresponding control pitch angles are depicted in Figure 8, where the hover trim variables reported in Table 3 are retrieved. Given the highly-classified nature of the data involved during the dronization process, the adopted first-guess controller gains are omitted.

5.2. Hardware–in–the–Loop Validation

The simulation setup described above is deployed to a HIL laboratory facility according to the scheme outlined in Figure 9. System components are set up as follows:

• The software developed in Matlab/Simulink for the mathematical modeling of helicopter dynamics and AHRS devices is automatically coded and deployed to a high-performance Real-Time Target Machine (RTTM) by Simulink Real-Time^TM tools. Solver frequency is set at 20 kHz, while AHRS model data are generated at 100 Hz. Software coding and deployment are performed through a host desktop PC, where the FlightGear open-source application is used to represent simulation data through a 3D graphical interface.
• The output of RTTM is provided via a dedicated standard industrial bus I/O module with two isolated ports. The first port is used to output the emulated AHRS data. The second port is used to generate repeatable control commands for HIL validation only, as if they were provided by the pilot on the ground.
• AHRS data and pilot commands from the RC device are the inputs to the onboard computer. At the time of the HIL experiments the RTC2, was already mounted on the helicopter as a pure acquisition device for an extensive campaign of manned flight tests. Hence, laboratory HIL tests were performed by arranging the RTC1 as an onboard computer. To emulate the presence of the radio modem, signals from the RC are converted from serial to standard industrial protocol by a micro-controller board equipped with a dedicated conversion shield. The control laws designed in Matlab/Simulink are coded and deployed to the RTC1 through the rugged laptop. The code developed for the onboard computer makes use of proprietary libraries for PID control implementation, acquisition and processing of input signals (including the application of Butterworth filters with order 1 and a cut-off frequency of 5 Hz to measured data), and real-time monitoring of selected variables.
• Control signals are acquired through a terminal board by a dedicated I/O module, a 16 bit analog input device selected to close the control loop. An ad hoc test bench is also provided where 1 EMA is controlled, in turn, by a voltage signal. Information about the linear motion are acquired and made available to evaluate the actuation performance.

Different maneuvers are performed during HIL simulations to validate the control strategy in Section 5.1 and the hardware implementation. A sample case is reported in what follows. Starting from a hovering condition, a step input ${\delta }_a^{\left(pilot\right)}=0.08$ is generated via the RTTM in order to reach a desired roll angle of 2 deg while keeping the other inputs unaltered. In Figure 10 and Figure 11, the results of the maneuver are reported in terms of variation with respect to the hover trim variables.

It can be noted that, for the same maneuver, the error between HIL and MIL simulations always remains bounded and smaller than $0.001$ deg (roll angle) and $0.005$ deg/s (roll rate). Furthermore, discretization and quantization effects of signals are investigated, which, however, do not affect controller efficacy. This and many other simulation tests validate the quality of the simulation software and the correct implementation of acquisition, actuation, and control system protocols in the presence of real flight hardware in a controlled environment prior to flight.

6. Flight Tests with the Unmanned Helicopter

After an extensive campaign of HIL simulations aiming at the fine-tuning of controller gains and the correct setup of hardware implementation, the helicopter is finally configured for unmanned flight tests and equipped with the ballistic parachute canopy. In order to simulate the presence of on board passengers, sandbags are put on the two seats, thus replicating the inertial configuration analyzed in Section 4, with the exception of the canopy. The campaign, performed in June 2018 at the airport of Oristano–Fenosu (Sardinia, Italy) in 4 days, is organized according to the following steps:

• Step 1. Direct control of onboard actuators by remote pilot commands, such that ${\delta }_a={\delta }_a^{\left(pilot\right)}$, ${\delta }_e={\delta }_e^{\left(pilot\right)}$, ${\delta }_c={\delta }_c^{\left(pilot\right)}$, and ${\delta }_p={\delta }_p^{\left(pilot\right)}$. This piloting configuration allows for validation of the overall actuation setup and represents a reversion mode in case of AHRS failure (manual mode).
• Step 2. The controller in Equation (27) is activated in order to stabilize the yaw rate. Different flight tests are performed and control gains are refined according to remote pilot recommendations, such that $k_p^{\left(r\right)}$ and $k_i^{\left(r\right)}$ are respectively increased by about $15\%$ and $14\%$ with respect to the first-guess values in Section 5.1.
• Step 3. Before activating the controllers in Equations (28) and (29), an intermediate test is performed in order to evaluate the damping contribution only provided by gains $k_d^{\left(\varphi \right)}$ and $k_d^{\left(\theta \right)}$ to the flying qualities about the roll and the pitch axis, respectively. To this end, the yaw rate is stabilized as in Step 2, while the direct control action of the pilot on lateral and longitudinal cyclic commands is supported by roll and pitch damper controllers, configured as follows:

  $${\delta }_a={\delta }_a^{\left(pilot\right)}+k_d^{\left(\varphi \right)}p$$

  (30)

  $${\delta }_e={\delta }_e^{\left(pilot\right)}+k_d^{\left(\theta \right)}q$$

  (31)
  At the end of Step 3, control gains are fine-tuned such that $k_d^{\left(\varphi \right)}$ and $k_d^{\left(\theta \right)}$ are respectively increased by about $2\%$ and $13\%$ with respect to the first-guess values.
• Step 4. The attitude controllers in Equations (28) and (29) are investigated, leaving the pilot with direct control of MR collective pitch only. Control gains are corrected such that $k_p^{\left(\theta \right)}$ and $k_i^{\left(\theta \right)}$ are respectively increased by $25\%$ and $60\%$ with respect to the precautionary small values proposed in Section 5.1. Finally, $k_p^{\left(\varphi \right)}$ and $k_i^{\left(\varphi \right)}$ are left unaltered.

Some flight data is reported, which describes the tests performed after Step 4 with the unmanned system in its definitive mission configuration.

In Figure 12, the commanded value of yaw rate, calculated as ${\xi }_r{\delta }_p^{\left(pilot\right)}$ (black line), is compared with the corresponding value measured by the AHRS (gray line). The data are expressed in deg/s and show the correlation between the desired and achieved attitude motion while the pilot performs oscillatory yawing maneuvers.

In Figure 13a,b the stabilization of roll and pitch angles is also analyzed over the same time period (80 s). In particular, roll angle oscillates with a standard deviation of $0.78$ deg about the mean value of $-1.91$ deg. Similar considerations hold for the pitch angle, characterized by a standard deviation of $0.74$ deg and a mean value of $0.35$ deg. If, on the one hand, the roll angle is consistent with the simulation results obtained in Table 3, the pitch angle shows major difference. This is caused by the presence of light tail wind and the fact that the inertial and aerodynamic configuration of the unmanned helicopter differs because of the presence of the parachute canopy over MRH. Collective command, characterized by a standard deviation of $0.01$, remains almost constant and equal to $0.66$ (corresponding to $13.67$ deg pitch angle).

The final experiment, performed on 22 June, is described in Figure 14, where helicopter trajectory is plotted in a $3D$ environment. Position data are obtained from GPS measurements provided by the AHRS and recorded by the RTC2. After the initial phase required for pre–flight checks and turbine engine warm up, the take–off occurs at time $t_0$. The climb phase to the height $h_1=330$ m is performed in $t_1-t_0=97$ s in the presence of South–West (SW) wind, with an average climb rate of about $3.4$ m/s. In particular, during the first 40 s the climb rate is stabilized at 2 m/s by pilot’s action, and then pushed to $4.5$ m/s until reaching the maximum height. At time $t_1$, the prescribed flight termination procedure is activated by switching-off the engine and commanding parachute ejection at time $t_2=t_1+4$ s. Complete parachute deployment is performed in about 5 s, at time $t_3=t_2+5$ s (see Figure 15). During the helicopter accelerated free fall the total height loss is $h_3-h_2=-146$ m, with an average vertical speed of $-16.2$ m/s. After $t=t_3$ the rate of descent stabilizes to a practically constant value of $7.5$ m/s until the helicopter safely lands at $t_4=t_3+27$ s. The effect of wind is visible in Figure 14, where helicopter trajectory deviates in the North–East direction and stops near the runway at about 285 m from the take–off point. Upon impact with the ground, acceleration peaks are recorded that fall within the parameters of crash tests in both the aeronautical and automotive sectors. Test data show that the system is likely to achieve its goal of saving lives, even at a lower altitude of just 150 m.

7. Conclusions

In the present paper, the complete procedure adopted to transform a light helicopter into an unmanned rotorcraft is described. By adopting the MBD approach, mission requirements were first outlined, and the design of the control system was addressed in terms of system architecture definition. Particular attention was devoted to the mathematical model of the helicopter and its subsystems, made on the basis of geometric, inertial, and aerodynamic data provided by the manufacturer and refined by identification techniques.

With the purpose of validating an innovative ballistic parachute rescue system, a closed-loop controller was developed to allow stable maneuvering in the field of view of a remote pilot. To this end, attitude stabilization algorithms were first tested in a Model-In-the-Loop environment. Furthermore, laboratory experiments allowed for (1) Hardware-In-the-Loop validation of involved equipment and (2) control software deployment on real-time target machines. Dedicated flight tests were performed to prove the effectiveness of the approach and the achievement of the desired closed–loop flying qualities. The final mission successfully showed the feasibility of the proposed termination procedure by securing a safe helicopter landing in the event of engine failure. The experiment allowed researchers to focus on the design and experimental validation of technologies at the core of future UAM, envisaging a more efficient, safe, and possibly sustainable exploitation of the vertical dimension.