# Helicopter Autorotation Trajectory Planning Method Using Functional Tensor-Train-Based Dynamic Programming Algorithms

Authors: Yunjie Wang; Chen Jiang; Shuai Deng; Haowen Wang
Journal: International Journal of Aerospace Engineering
 (2020)

---

## Abstract

Helicopter autorotation trajectory planning problems have been dealt within computationally expensive optimal control algorithms. This paper presents an efficient helicopter autorotation trajectory planning method, using functional tensor-train- (FT-) based dynamic programming (DP) algorithms. The autorotation trajectory planning method is shown real-time feasible, which involves general helicopter autorotation dynamics at the same time. To validate the dynamic feasibility of the trajectories, a trajectory-tracking controller using active disturbance rejection control (ADRC) is designed to ensure a helicopter model tracks the trajectories. Finally, a helicopter autorotation simulation with a six-degree-of-freedom high-fidelity multibody-based helicopter model is demonstrated for validation.


---

## Body

## 1. Introduction

Autorotation is a primary measure for helicopters, whether for manned or unmanned helicopters, to land safely after engine power failure. In autorotation, the main rotor is driven by the upward airflow through the rotor, making the flight similar to a gliding fixed-wing aircraft. To achieve a safe autorotation landing, unique control strategies are needed [1]. Generally speaking, the collective pitch should be carefully handled to maintain a sufficient and steady rotor speed before reaching the ground.Studies on means of achieving autorotation have been focused on solving optimal control problems [2–4]. Traditionally, simplified helicopter dynamic models are used, including 2-D point-mass models [4–6], three-degree-of-freedom rigid-body models [7, 8], and low-order six-degree-of-freedom rigid-body models [3]. With the dynamic model formulated, optimal autorotation problems are solved by numerical methods. Gradient algorithms, such as sequential gradient-restoration algorithms, are used in [5, 7, 9]. Other methods, such as direct methods, discretize the problem first and turn it into a nonlinear programming problem. Then, the nonlinear programming problem can be solved using algorithms such as sequential quadratic programming (SQP), as demonstrated in [6, 10, 11]. In general, these algorithms are off-line, computationally expensive for current computing power, which is for now unrealistic for on-line usages. Furthermore, there is a lack of validations either through high-fidelity simulations or experiments in the above studies.Recently, Taamallah [12] proposed the first real-time feasible, model-based trajectory planning method and designed a model-based trajectory-tracking controller to ensure the helicopter tracks the trajectories. In specific, they employed optimal planning based on differential flatness, assuming a helicopter as a rigid body. The trajectory planning problem is solved, regarding all the rotor forces and moments as plant inputs of the rigid body. Taamallah’s work provides novel directions for online trajectory planning of autorotation problems. However, as a consequence, forces and moments are simplified and decoupled with helicopter dynamics in such methods. It is necessary to mention that, in autorotation procedures, the main rotor’s ability to generate forces is highly restrained by the states of helicopters, especially by the main rotor’s rotating speed and the inflow state of the rotor [13]. The rotor speed is the crucial factor that determines whether a trajectory leads to a safe autorotation landing. Thus, the dynamics of autorotation is of great significance for trajectory planning and should be involved explicitly.In this paper, we present a real-time feasible autorotation trajectory planning method using functional tensor-train- (FT-) based dynamic programming (DP) algorithms, which ensures that general autorotation dynamics is satisfied along the trajectory. For validation of the dynamic feasibility of the trajectories, we also present a trajectory-tracking controller based on active disturbance rejection control (ADRC) to make a helicopter model track the trajectories. The validation of the trajectories using the controller is then implemented on a six-degree-of-freedom, high-fidelity, multibody-based helicopter simulation model.The trajectory planning method using functional tensor-train- (FT-) based dynamic programming (DP) algorithms will be presented in Section2. The ADRC-based trajectory-tracking controller is described in Section 3. In Section 4, we describe the six-degree-of-freedom, high-fidelity, multibody helicopter simulation model using the Tsinghua Rotorcraft Utility Simulation Tool (TRUST). In Section 5, autorotation trajectories are demonstrated with various initial conditions, and simulation results for validation are demonstrated and discussed. Finally, conclusions are presented in Section 6.


## 2. Autorotation Trajectory Planning Using FT-Based DP Algorithms

As mentioned in Section1, traditional trajectory planning methods for autorotation have been dealt within computationally expensive off-line algorithms. In 2017, the first real-time feasible trajectory planning method for autorotation was demonstrated in [12], which is based on differential flatness of the rigid-body dynamics. Such methods simplify forces as direct inputs, thus leaving the helicopter autorotation dynamics not included during trajectory planning procedures. In this section, we introduce a real-time feasible trajectory planning method, using functional tensor-train- (FT-) [14] based dynamic programming (DP) algorithms, which guarantees a strict satisfaction of helicopter dynamics along the trajectory.


### 2.1. FT-Based DP Algorithms

Functional tensor-train-based (FT-based) dynamic programming (DP) algorithms are newly proposed algorithms for solving high-dimensional stochastic optimal control (SOC) problems, which are mainly discounted-cost infinite-horizon Markov decision process (MDP) problems. Here, we give a brief review of the FT-based DP algorithms, and details can be found in [15, 16].Consider a system described by stochastic differential equations (SDE) as follows:(1)xt=bt,xt,utdt+σxtdwt,where xt∈Rn is the state vector, ut∈Rm is the control input, wt is a vector of independent unit Wiener processes, and b⋅ and σ⋅ are generally nonlinear functions. The cost functional is defined as
(2)Jxt0,u=E∫t0τe−βsgs,xs,usds+χτ<∞e−βτψJτ,xτ,where β is a discount factor, χτ<∞ is the indicator function of the state boundary, and g⋅ and ψJ⋅ are stage cost function and terminal cost function. The SOC problem is to find a control ut within a specified set on the time interval t0,τ, such that the cost Jxt0,u is minimized.Next, using the MCA method [17], continuous SOC problems are discretized into the discrete Markov decision processes (MDPs). The discretized problem is turned into searching for a value function that satisfies the following recursive equation:
(3)vz=minμgz,u+γ∑z′pz,z′uvz′,where vz is the optimal discretized value function and pz,z′u is the transition probability function.Then, a discounted-cost infinite-horizon MDP is formulated and can be solved by the FT-based DP algorithms, which are FT-based value iteration algorithm, FT-based policy iteration algorithm, and FT-based multilevel algorithm, respectively.For the traditional discrete-state Markov decision processes, computational requirements grow exponentially with dimensionality. For example, if a MDP has 8 dimensions, and each dimension has a discretization of 10 points, such a problem involves a search space of108 points. In order to mitigate such curse of dimensionality, FT-based DP algorithms use low-rank-functions, namely functional tensor-train, to represent value functions. The basic idea of function-train (FT) is to make a continuous analogue [14] of the tensor-train decomposition [18]. To be specific, it is a continuous version of tensor-train cross-approximation (TT-cross-approximation) [19], with the formulation as follows:
(4)fx1,⋯,xd=∏i=1dFixi,where fx1,⋯,xd is a d-dimensional multivariable function. Fixi is a set of univariate functions, which are also called cores:
(5)Fixi=f1,1ixi⋯f1,riixi⋮⋱⋮fri−1,riixi⋯fri−1,rikxi,where ri are the FT ranks evaluated by a continuous version of TT-rounding [14] and fi,jk are univariate hat functions.With methods described above, the exponentially growing computational complexity ofOnd for a typical dynamic programming problem is compressed to a polynomially growing complexity of
(6)Odnr2κnop+d2nr2+dnr3,where κ is the operations during value evaluation, nop is the step of operations within each Bellman equation step.To summarize, this method formulates a SOC problem as a dynamic programming problem and implements a continuous tensor decomposition method to compress such a problem, resulting in significant improvements of computing speed and storage savings. Such method has been shown to be able to work in real-time [15].


### 2.2. Formulating Autorotation Trajectory Planning Problems within FT-Based DP Algorithms

In this part, we describe a general method for solving autorotation trajectory planning problems in the framework of FT-based DP algorithms. First of all, we need to describe the dynamics of a helicopter in a SOC form, as equation (1) shows.


#### 2.2.1. Helicopter Dynamic Model Formulation

The helicopter dynamic function in the SOC equations is given by a nonlinear three-degree-of-freedom rigid-body helicopter autorotation dynamic model. Note that in this paper, two different helicopter dynamic models are described. The first one is the nonlinear three-degree-of-freedom rigid-body helicopter dynamic model, which is used in this section as the system dynamics in equation (1). The second one is the six-degree-of-freedom, high-fidelity, multibody-based helicopter simulation model, which will be used as a validation model in simulation and will be described in Section 5. The three-degree-of-freedom dynamic model is chosen here for computation considerations, and such model is capable of predicting steady collective pitch manipulations and rotor power consumption [13].States variables for the nonlinear three-degree-of-freedom rigid-body helicopter dynamic model are chosen as(7)s=z,vx,vz,θ,q,Ω,which include height z, horizontal velocity vx, rate of descent vz, pitch angle θ, pitch rate q, and rotating speed of the main rotor Ω. The control variables are collective pitch θ0 and longitudinal cyclic pitch θ1s. The dynamic model equations are formulated as follows [20]:
(8)ż=vzv̇x=−1mTBsinθ+HBcosθ+DFcosαFv̇z=−1mTBcosθ−HBsinθ+DFsinαF+gθ̇=qq̇=1IFYMY−TB⋅lR+HB⋅hR−LHcosθ−γ⋅lHΩ̇=−1IRQ,where TB and HB are rotor forces resolved in the body axes. Note that all the axes are defined in accordance with [13]. MY and Q are pitch moment and rotor torque, respectively. DF and LH are the drag of the fuselage and the lift of the horizontal tail resolved in the body axes. αF and γ are the angle of attack of the fuselage and the flight-path angle. hR and lR are vertical and horizontal distance of rotor hub to the center of gravity. lH is the horizontal distance of the horizon tail from the center of gravity. Details of forces generated by the fuselage and the horizontal tail can be seen in [20].Next, the rotor performance during autorotation is given by the following procedures. The induced velocityν is calculated by [9]:
(9)ν=κνhfIfG,where κ is the empirical correction factor of nonuniform inflow and νh is the induced velocity at hover. fI is the induced velocity parameter which has included the influence of the vortex-ring state, given by [9]:
(10)fI=1.0/v¯xb2+−v¯zb+fI2,if−2v¯zb+32+v¯xb2≥1.0−v¯zb0.373v¯zb2+0.598v¯xb2−1.991,otherwise,where v¯xb and v¯zb are normalized velocities resolved in body axes:
(11)v¯xbv¯zb=1vhcosθ−sinθsinθcosθvxvz.The ground effect factorfG is expressed as [13]
(12)fG=1−11−R/4h2,where h=z+z0 and z0 is the rotor height when the aircraft is on the ground. Note that we ignore the influence of horizontal speed on the ground effect, and such assumption is reasonable for a typical autorotation procedure.The rotor thrust coefficientCT is expressed as [13]
(13)CT=12aσ13+12v¯xb2θ0−34θtwB3+141+v¯xb2B4−12λ1+12v¯xb2B2,where a is the blade lift-curve slope, σ is the rotor solidity ratio, θtw is the blade twist, and B is the blade tip-loss factor.The rotor drag force coefficientCH is made up of the rotor profile drag CH0 and the induced drag CHi. The rotor pitch moment coefficient is CMY. CH and CMY are expressed as follows:
(14)CH=CH0+CHiCH0=cdσ83v¯xb+1.98v¯xb2.7CHi=aσ2−13βc+12v¯xbλθ0−18θtwv¯xbλ+−14v¯xbβ1c+14λθ1s+34λβ1c+16β0β1c+14β02+β1c2CMY=−12εnbΩ2β1cSb,where cd is the blade’s mean profile drag coefficient, ε is the flap offset, nb is the number of blades, and Sb is the blade’s first moment about the flap hinge. β0 is the coning angle of rotor disk, and β1c is the first harmonic cyclic flap, given by
(15)β0=γ18θ0−34θtw1+v¯xb2B4+θtw101+56v¯xb2B5−16λ−v¯xbθ1sB3−32gΩ2Rβ1c=−8/3θ0−3/4λ−v¯xbθ1s1−1/2v¯xb2−θ1s,where λ is the normalized inflow velocity, expressed as
(16)λ=ν/ΩR−v¯zb.The power coefficient required by the main rotorCP is equal to the rotor torque coefficient CQ, which is made up of the rotor induced power, the rotor profile power, the rotor parasite power, and the rotor climb power. A general expression is given by [13]
(17)CP=κCTνΩR+18σcd1+4.6v¯xb2+12v¯xb3fπR2−v¯zbCT,where f is the fuselage equivalent parasite drag area.


#### 2.2.2. Formulation of SOC Problem for Autorotation Trajectory Planning

Following the helicopter autorotation dynamic model described above, and for numerical considerations [5], the stochastic differential equations are obtained by normalizing equations of (8). The normalized states and inputs are as follows:
(18)x1=110Rz,x2=100Ω0Rvx,x3=100Ω0Rvz,x4=10θ,x5=1000Ω0q,x6=1Ω0Ωu1=10θ0,u2=10θ1s.In addition, the Jacobian matrix of control inputs∂ẋi/∂uj,i=1,⋯,6,j=1,2 is calculated by numerical difference method for any state x and input u. The diffusion function σ⋅ is set to 10−9 for each term.Next, the cost functions are designed. The stage costgx,u is expressed by
(19)gx,u=W1x12+W2x22+W3x32+W4x42+W61−x62+W7u12+W8u22.The terminal costψJx is given by setting an absorbing region [15], which encourages a safe landing, as shown below:
(20)ψJx=0ifsafelylandedWT1x12+WT2x2−x2glide2+WT3x32+WT4x42+WT5x52otherwise,where W· and WT· denote weighting factors of stage cost and terminal cost, respectively. We add cost functions to control inputs in order to avoid fierce manipulations. x2glide is the average horizontal speed, estimated by empirical flight test results and numerical simulation results. A safe autorotation landing is defined by
(21)vzτ≤vzmax,vxτ≤vxmax,θτ≤θmax,qτ≤qmax,z=0.The design of cost functions is shown to be important for such autorotation problems, and parameters are needed to be adjusted by numerical experiments.


#### 2.2.3. Trajectory Generation Using FT-Based Algorithms

After the problem is formulated in the form of a SOC problem in Section2.2.2, solutions are then obtained using FT-based algorithms.Although FT-based algorithms are proposed for stochastic optimal control problems, we slightly alter such algorithms to generate a trajectory, rather than to obtain control inputs directly. The reason is that even though the nonlinear dynamics of helicopter described in Section2.2 gives nice results in terms of predicting steady collective pitch manipulations and rotor power consumption, such model is not able to make good predictions of the helicopter dynamic response [21]. Therefore, the control inputs generated by the SOC controller are not adopted, and we make use of the trajectory instead. As mentioned before, the results of rotor power and collective pitch along the trajectory can be regarded reasonable.We employ the FT-based one-way multigrid algorithm proposed in [15, 22]. For each discretization level, Niter=100 steps of FT-based policy iterations are applied. The trajectory of normalized states x is obtained by an integration procedure of given initial conditions to the controller:
(22)xti=∫ti−1tibxti−1,uxi−1,ti∈t0,τ,xt0=x0.In specific, a fourth-order order Runge-Kutta method is implemented. Thus, preliminary trajectories are obtained by makingx unnormalized by equation (18). Because the trajectories are generated by integrations based on the dynamic equations of the helicopter, they satisfy the specified helicopter dynamics by nature.However, although such trajectories are generated by solving dynamic equations, the trajectories still need to be smoothed. There are two main reasons. The first reason is that for real-time realizability, the integration time step (for example, 0.05 s in our study) is not small enough for a controller to generate sufficient differential information. The second reason is that although we have cost functions on control inputs, the control values still result in discontinuous variations, making the trajectory not smooth enough for our controller. Thus, we apply a fast interpolation method, i.e., cubic Hermite interpolation, to obtain smooth trajectories.


## 2.1. FT-Based DP Algorithms

Functional tensor-train-based (FT-based) dynamic programming (DP) algorithms are newly proposed algorithms for solving high-dimensional stochastic optimal control (SOC) problems, which are mainly discounted-cost infinite-horizon Markov decision process (MDP) problems. Here, we give a brief review of the FT-based DP algorithms, and details can be found in [15, 16].Consider a system described by stochastic differential equations (SDE) as follows:(1)xt=bt,xt,utdt+σxtdwt,where xt∈Rn is the state vector, ut∈Rm is the control input, wt is a vector of independent unit Wiener processes, and b⋅ and σ⋅ are generally nonlinear functions. The cost functional is defined as
(2)Jxt0,u=E∫t0τe−βsgs,xs,usds+χτ<∞e−βτψJτ,xτ,where β is a discount factor, χτ<∞ is the indicator function of the state boundary, and g⋅ and ψJ⋅ are stage cost function and terminal cost function. The SOC problem is to find a control ut within a specified set on the time interval t0,τ, such that the cost Jxt0,u is minimized.Next, using the MCA method [17], continuous SOC problems are discretized into the discrete Markov decision processes (MDPs). The discretized problem is turned into searching for a value function that satisfies the following recursive equation:
(3)vz=minμgz,u+γ∑z′pz,z′uvz′,where vz is the optimal discretized value function and pz,z′u is the transition probability function.Then, a discounted-cost infinite-horizon MDP is formulated and can be solved by the FT-based DP algorithms, which are FT-based value iteration algorithm, FT-based policy iteration algorithm, and FT-based multilevel algorithm, respectively.For the traditional discrete-state Markov decision processes, computational requirements grow exponentially with dimensionality. For example, if a MDP has 8 dimensions, and each dimension has a discretization of 10 points, such a problem involves a search space of108 points. In order to mitigate such curse of dimensionality, FT-based DP algorithms use low-rank-functions, namely functional tensor-train, to represent value functions. The basic idea of function-train (FT) is to make a continuous analogue [14] of the tensor-train decomposition [18]. To be specific, it is a continuous version of tensor-train cross-approximation (TT-cross-approximation) [19], with the formulation as follows:
(4)fx1,⋯,xd=∏i=1dFixi,where fx1,⋯,xd is a d-dimensional multivariable function. Fixi is a set of univariate functions, which are also called cores:
(5)Fixi=f1,1ixi⋯f1,riixi⋮⋱⋮fri−1,riixi⋯fri−1,rikxi,where ri are the FT ranks evaluated by a continuous version of TT-rounding [14] and fi,jk are univariate hat functions.With methods described above, the exponentially growing computational complexity ofOnd for a typical dynamic programming problem is compressed to a polynomially growing complexity of
(6)Odnr2κnop+d2nr2+dnr3,where κ is the operations during value evaluation, nop is the step of operations within each Bellman equation step.To summarize, this method formulates a SOC problem as a dynamic programming problem and implements a continuous tensor decomposition method to compress such a problem, resulting in significant improvements of computing speed and storage savings. Such method has been shown to be able to work in real-time [15].


## 2.2. Formulating Autorotation Trajectory Planning Problems within FT-Based DP Algorithms

In this part, we describe a general method for solving autorotation trajectory planning problems in the framework of FT-based DP algorithms. First of all, we need to describe the dynamics of a helicopter in a SOC form, as equation (1) shows.


### 2.2.1. Helicopter Dynamic Model Formulation

The helicopter dynamic function in the SOC equations is given by a nonlinear three-degree-of-freedom rigid-body helicopter autorotation dynamic model. Note that in this paper, two different helicopter dynamic models are described. The first one is the nonlinear three-degree-of-freedom rigid-body helicopter dynamic model, which is used in this section as the system dynamics in equation (1). The second one is the six-degree-of-freedom, high-fidelity, multibody-based helicopter simulation model, which will be used as a validation model in simulation and will be described in Section 5. The three-degree-of-freedom dynamic model is chosen here for computation considerations, and such model is capable of predicting steady collective pitch manipulations and rotor power consumption [13].States variables for the nonlinear three-degree-of-freedom rigid-body helicopter dynamic model are chosen as(7)s=z,vx,vz,θ,q,Ω,which include height z, horizontal velocity vx, rate of descent vz, pitch angle θ, pitch rate q, and rotating speed of the main rotor Ω. The control variables are collective pitch θ0 and longitudinal cyclic pitch θ1s. The dynamic model equations are formulated as follows [20]:
(8)ż=vzv̇x=−1mTBsinθ+HBcosθ+DFcosαFv̇z=−1mTBcosθ−HBsinθ+DFsinαF+gθ̇=qq̇=1IFYMY−TB⋅lR+HB⋅hR−LHcosθ−γ⋅lHΩ̇=−1IRQ,where TB and HB are rotor forces resolved in the body axes. Note that all the axes are defined in accordance with [13]. MY and Q are pitch moment and rotor torque, respectively. DF and LH are the drag of the fuselage and the lift of the horizontal tail resolved in the body axes. αF and γ are the angle of attack of the fuselage and the flight-path angle. hR and lR are vertical and horizontal distance of rotor hub to the center of gravity. lH is the horizontal distance of the horizon tail from the center of gravity. Details of forces generated by the fuselage and the horizontal tail can be seen in [20].Next, the rotor performance during autorotation is given by the following procedures. The induced velocityν is calculated by [9]:
(9)ν=κνhfIfG,where κ is the empirical correction factor of nonuniform inflow and νh is the induced velocity at hover. fI is the induced velocity parameter which has included the influence of the vortex-ring state, given by [9]:
(10)fI=1.0/v¯xb2+−v¯zb+fI2,if−2v¯zb+32+v¯xb2≥1.0−v¯zb0.373v¯zb2+0.598v¯xb2−1.991,otherwise,where v¯xb and v¯zb are normalized velocities resolved in body axes:
(11)v¯xbv¯zb=1vhcosθ−sinθsinθcosθvxvz.The ground effect factorfG is expressed as [13]
(12)fG=1−11−R/4h2,where h=z+z0 and z0 is the rotor height when the aircraft is on the ground. Note that we ignore the influence of horizontal speed on the ground effect, and such assumption is reasonable for a typical autorotation procedure.The rotor thrust coefficientCT is expressed as [13]
(13)CT=12aσ13+12v¯xb2θ0−34θtwB3+141+v¯xb2B4−12λ1+12v¯xb2B2,where a is the blade lift-curve slope, σ is the rotor solidity ratio, θtw is the blade twist, and B is the blade tip-loss factor.The rotor drag force coefficientCH is made up of the rotor profile drag CH0 and the induced drag CHi. The rotor pitch moment coefficient is CMY. CH and CMY are expressed as follows:
(14)CH=CH0+CHiCH0=cdσ83v¯xb+1.98v¯xb2.7CHi=aσ2−13βc+12v¯xbλθ0−18θtwv¯xbλ+−14v¯xbβ1c+14λθ1s+34λβ1c+16β0β1c+14β02+β1c2CMY=−12εnbΩ2β1cSb,where cd is the blade’s mean profile drag coefficient, ε is the flap offset, nb is the number of blades, and Sb is the blade’s first moment about the flap hinge. β0 is the coning angle of rotor disk, and β1c is the first harmonic cyclic flap, given by
(15)β0=γ18θ0−34θtw1+v¯xb2B4+θtw101+56v¯xb2B5−16λ−v¯xbθ1sB3−32gΩ2Rβ1c=−8/3θ0−3/4λ−v¯xbθ1s1−1/2v¯xb2−θ1s,where λ is the normalized inflow velocity, expressed as
(16)λ=ν/ΩR−v¯zb.The power coefficient required by the main rotorCP is equal to the rotor torque coefficient CQ, which is made up of the rotor induced power, the rotor profile power, the rotor parasite power, and the rotor climb power. A general expression is given by [13]
(17)CP=κCTνΩR+18σcd1+4.6v¯xb2+12v¯xb3fπR2−v¯zbCT,where f is the fuselage equivalent parasite drag area.


### 2.2.2. Formulation of SOC Problem for Autorotation Trajectory Planning

Following the helicopter autorotation dynamic model described above, and for numerical considerations [5], the stochastic differential equations are obtained by normalizing equations of (8). The normalized states and inputs are as follows:
(18)x1=110Rz,x2=100Ω0Rvx,x3=100Ω0Rvz,x4=10θ,x5=1000Ω0q,x6=1Ω0Ωu1=10θ0,u2=10θ1s.In addition, the Jacobian matrix of control inputs∂ẋi/∂uj,i=1,⋯,6,j=1,2 is calculated by numerical difference method for any state x and input u. The diffusion function σ⋅ is set to 10−9 for each term.Next, the cost functions are designed. The stage costgx,u is expressed by
(19)gx,u=W1x12+W2x22+W3x32+W4x42+W61−x62+W7u12+W8u22.The terminal costψJx is given by setting an absorbing region [15], which encourages a safe landing, as shown below:
(20)ψJx=0ifsafelylandedWT1x12+WT2x2−x2glide2+WT3x32+WT4x42+WT5x52otherwise,where W· and WT· denote weighting factors of stage cost and terminal cost, respectively. We add cost functions to control inputs in order to avoid fierce manipulations. x2glide is the average horizontal speed, estimated by empirical flight test results and numerical simulation results. A safe autorotation landing is defined by
(21)vzτ≤vzmax,vxτ≤vxmax,θτ≤θmax,qτ≤qmax,z=0.The design of cost functions is shown to be important for such autorotation problems, and parameters are needed to be adjusted by numerical experiments.


### 2.2.3. Trajectory Generation Using FT-Based Algorithms

After the problem is formulated in the form of a SOC problem in Section2.2.2, solutions are then obtained using FT-based algorithms.Although FT-based algorithms are proposed for stochastic optimal control problems, we slightly alter such algorithms to generate a trajectory, rather than to obtain control inputs directly. The reason is that even though the nonlinear dynamics of helicopter described in Section2.2 gives nice results in terms of predicting steady collective pitch manipulations and rotor power consumption, such model is not able to make good predictions of the helicopter dynamic response [21]. Therefore, the control inputs generated by the SOC controller are not adopted, and we make use of the trajectory instead. As mentioned before, the results of rotor power and collective pitch along the trajectory can be regarded reasonable.We employ the FT-based one-way multigrid algorithm proposed in [15, 22]. For each discretization level, Niter=100 steps of FT-based policy iterations are applied. The trajectory of normalized states x is obtained by an integration procedure of given initial conditions to the controller:
(22)xti=∫ti−1tibxti−1,uxi−1,ti∈t0,τ,xt0=x0.In specific, a fourth-order order Runge-Kutta method is implemented. Thus, preliminary trajectories are obtained by makingx unnormalized by equation (18). Because the trajectories are generated by integrations based on the dynamic equations of the helicopter, they satisfy the specified helicopter dynamics by nature.However, although such trajectories are generated by solving dynamic equations, the trajectories still need to be smoothed. There are two main reasons. The first reason is that for real-time realizability, the integration time step (for example, 0.05 s in our study) is not small enough for a controller to generate sufficient differential information. The second reason is that although we have cost functions on control inputs, the control values still result in discontinuous variations, making the trajectory not smooth enough for our controller. Thus, we apply a fast interpolation method, i.e., cubic Hermite interpolation, to obtain smooth trajectories.


## 2.2.1. Helicopter Dynamic Model Formulation

The helicopter dynamic function in the SOC equations is given by a nonlinear three-degree-of-freedom rigid-body helicopter autorotation dynamic model. Note that in this paper, two different helicopter dynamic models are described. The first one is the nonlinear three-degree-of-freedom rigid-body helicopter dynamic model, which is used in this section as the system dynamics in equation (1). The second one is the six-degree-of-freedom, high-fidelity, multibody-based helicopter simulation model, which will be used as a validation model in simulation and will be described in Section 5. The three-degree-of-freedom dynamic model is chosen here for computation considerations, and such model is capable of predicting steady collective pitch manipulations and rotor power consumption [13].States variables for the nonlinear three-degree-of-freedom rigid-body helicopter dynamic model are chosen as(7)s=z,vx,vz,θ,q,Ω,which include height z, horizontal velocity vx, rate of descent vz, pitch angle θ, pitch rate q, and rotating speed of the main rotor Ω. The control variables are collective pitch θ0 and longitudinal cyclic pitch θ1s. The dynamic model equations are formulated as follows [20]:
(8)ż=vzv̇x=−1mTBsinθ+HBcosθ+DFcosαFv̇z=−1mTBcosθ−HBsinθ+DFsinαF+gθ̇=qq̇=1IFYMY−TB⋅lR+HB⋅hR−LHcosθ−γ⋅lHΩ̇=−1IRQ,where TB and HB are rotor forces resolved in the body axes. Note that all the axes are defined in accordance with [13]. MY and Q are pitch moment and rotor torque, respectively. DF and LH are the drag of the fuselage and the lift of the horizontal tail resolved in the body axes. αF and γ are the angle of attack of the fuselage and the flight-path angle. hR and lR are vertical and horizontal distance of rotor hub to the center of gravity. lH is the horizontal distance of the horizon tail from the center of gravity. Details of forces generated by the fuselage and the horizontal tail can be seen in [20].Next, the rotor performance during autorotation is given by the following procedures. The induced velocityν is calculated by [9]:
(9)ν=κνhfIfG,where κ is the empirical correction factor of nonuniform inflow and νh is the induced velocity at hover. fI is the induced velocity parameter which has included the influence of the vortex-ring state, given by [9]:
(10)fI=1.0/v¯xb2+−v¯zb+fI2,if−2v¯zb+32+v¯xb2≥1.0−v¯zb0.373v¯zb2+0.598v¯xb2−1.991,otherwise,where v¯xb and v¯zb are normalized velocities resolved in body axes:
(11)v¯xbv¯zb=1vhcosθ−sinθsinθcosθvxvz.The ground effect factorfG is expressed as [13]
(12)fG=1−11−R/4h2,where h=z+z0 and z0 is the rotor height when the aircraft is on the ground. Note that we ignore the influence of horizontal speed on the ground effect, and such assumption is reasonable for a typical autorotation procedure.The rotor thrust coefficientCT is expressed as [13]
(13)CT=12aσ13+12v¯xb2θ0−34θtwB3+141+v¯xb2B4−12λ1+12v¯xb2B2,where a is the blade lift-curve slope, σ is the rotor solidity ratio, θtw is the blade twist, and B is the blade tip-loss factor.The rotor drag force coefficientCH is made up of the rotor profile drag CH0 and the induced drag CHi. The rotor pitch moment coefficient is CMY. CH and CMY are expressed as follows:
(14)CH=CH0+CHiCH0=cdσ83v¯xb+1.98v¯xb2.7CHi=aσ2−13βc+12v¯xbλθ0−18θtwv¯xbλ+−14v¯xbβ1c+14λθ1s+34λβ1c+16β0β1c+14β02+β1c2CMY=−12εnbΩ2β1cSb,where cd is the blade’s mean profile drag coefficient, ε is the flap offset, nb is the number of blades, and Sb is the blade’s first moment about the flap hinge. β0 is the coning angle of rotor disk, and β1c is the first harmonic cyclic flap, given by
(15)β0=γ18θ0−34θtw1+v¯xb2B4+θtw101+56v¯xb2B5−16λ−v¯xbθ1sB3−32gΩ2Rβ1c=−8/3θ0−3/4λ−v¯xbθ1s1−1/2v¯xb2−θ1s,where λ is the normalized inflow velocity, expressed as
(16)λ=ν/ΩR−v¯zb.The power coefficient required by the main rotorCP is equal to the rotor torque coefficient CQ, which is made up of the rotor induced power, the rotor profile power, the rotor parasite power, and the rotor climb power. A general expression is given by [13]
(17)CP=κCTνΩR+18σcd1+4.6v¯xb2+12v¯xb3fπR2−v¯zbCT,where f is the fuselage equivalent parasite drag area.


## 2.2.2. Formulation of SOC Problem for Autorotation Trajectory Planning

Following the helicopter autorotation dynamic model described above, and for numerical considerations [5], the stochastic differential equations are obtained by normalizing equations of (8). The normalized states and inputs are as follows:
(18)x1=110Rz,x2=100Ω0Rvx,x3=100Ω0Rvz,x4=10θ,x5=1000Ω0q,x6=1Ω0Ωu1=10θ0,u2=10θ1s.In addition, the Jacobian matrix of control inputs∂ẋi/∂uj,i=1,⋯,6,j=1,2 is calculated by numerical difference method for any state x and input u. The diffusion function σ⋅ is set to 10−9 for each term.Next, the cost functions are designed. The stage costgx,u is expressed by
(19)gx,u=W1x12+W2x22+W3x32+W4x42+W61−x62+W7u12+W8u22.The terminal costψJx is given by setting an absorbing region [15], which encourages a safe landing, as shown below:
(20)ψJx=0ifsafelylandedWT1x12+WT2x2−x2glide2+WT3x32+WT4x42+WT5x52otherwise,where W· and WT· denote weighting factors of stage cost and terminal cost, respectively. We add cost functions to control inputs in order to avoid fierce manipulations. x2glide is the average horizontal speed, estimated by empirical flight test results and numerical simulation results. A safe autorotation landing is defined by
(21)vzτ≤vzmax,vxτ≤vxmax,θτ≤θmax,qτ≤qmax,z=0.The design of cost functions is shown to be important for such autorotation problems, and parameters are needed to be adjusted by numerical experiments.


## 2.2.3. Trajectory Generation Using FT-Based Algorithms

After the problem is formulated in the form of a SOC problem in Section2.2.2, solutions are then obtained using FT-based algorithms.Although FT-based algorithms are proposed for stochastic optimal control problems, we slightly alter such algorithms to generate a trajectory, rather than to obtain control inputs directly. The reason is that even though the nonlinear dynamics of helicopter described in Section2.2 gives nice results in terms of predicting steady collective pitch manipulations and rotor power consumption, such model is not able to make good predictions of the helicopter dynamic response [21]. Therefore, the control inputs generated by the SOC controller are not adopted, and we make use of the trajectory instead. As mentioned before, the results of rotor power and collective pitch along the trajectory can be regarded reasonable.We employ the FT-based one-way multigrid algorithm proposed in [15, 22]. For each discretization level, Niter=100 steps of FT-based policy iterations are applied. The trajectory of normalized states x is obtained by an integration procedure of given initial conditions to the controller:
(22)xti=∫ti−1tibxti−1,uxi−1,ti∈t0,τ,xt0=x0.In specific, a fourth-order order Runge-Kutta method is implemented. Thus, preliminary trajectories are obtained by makingx unnormalized by equation (18). Because the trajectories are generated by integrations based on the dynamic equations of the helicopter, they satisfy the specified helicopter dynamics by nature.However, although such trajectories are generated by solving dynamic equations, the trajectories still need to be smoothed. There are two main reasons. The first reason is that for real-time realizability, the integration time step (for example, 0.05 s in our study) is not small enough for a controller to generate sufficient differential information. The second reason is that although we have cost functions on control inputs, the control values still result in discontinuous variations, making the trajectory not smooth enough for our controller. Thus, we apply a fast interpolation method, i.e., cubic Hermite interpolation, to obtain smooth trajectories.


## 3. Autorotation Trajectory-Tracking Controller Based on ADRC

In order to show the practicability of trajectories generated by FT-based DP methods, we demonstrate a trajectory-tracking controller based on active disturbance rejection control (ADRC) and use the controller to make a high-fidelity helicopter model track the trajectories. Tracking of position and velocity is vital for a successful autorotation landing, and the dynamics during autorotation may be more complicated than the formulations described in Section2.2. Besides, there also exist unexpected disturbances during real flights. Thus, we use active disturbance rejection control (ADRC) methods to design the tracking controller. Active disturbance rejection control (ADRC) is proposed by Han [23], and the ADRC controller is capable of estimating inner modelling errors or outer disturbances, thus making compensation accordingly.In specific, we implement the ADRC-based trajectory-tracking controller described in [24], and such controller has been successfully validated through flight tests [24, 25]. The structure of the controller is shown in Figure 1: where xtraj and vtraj are position vectors and velocity vectors of the trajectory and ur represents the reference velocities that needed to be tracked by the inner loop controller.Figure 1
Trajectory-tracking controller structure.The order of the ADRC controller differs with different channels. For the autorotation application, we use a 3rd-order ADRC controller for forward, vertical, lateral channels, and a PI controller for yaw channel. In general, an ADRC controller is made up of a tracking differentiator (TD), an extended state observer (ESO), and a nonlinear state error feedback (NLSEF). Disturbances or modelling errors of the system are observed by the ESO, and then NLSEF is utilized to restrain them. A typical architecture of an ADRC controller is shown in Figure2.Figure 2
ADRC architecture.Wherev is the input signal, b0 is the input gain factor, z· is the output of the ESO. For brevity, only the 3rd-order ADRC controller is presented here. The TD of a 3rd-order is formulated as [25]
(23)fh1=fhanv1−v0,v2,r,h0v1=v1+h0⋅v2v2=v2+h0⋅fh1fh2=fhanv2,v3,r,h0v3=v3+h0⋅fh2,where fhan⋅ is an optimal control function defined in [23] and r and h0 are controller parameters. The 3rd-order ESO is formulated as
(24)e=z1−yz1=z1+hz2−β01ez2=z2+hz3−β02⋅fale,0.5,δz3=z3+hz4−β03⋅fale,0.25,δ+b0uz4=z4+h−β04⋅fale,0.125,δ,where h is the integration step of the system, β· is the observer gain, and function fal⋅ is defined in [23]. The 3rd-order NLSEF is formulated as
(25)e1=v1−z1e2=v2−z2e3=v3−z3u0=β1fale1,α1,δ+β2fale2,α2,δ+β3fale3,α3,δu=u0−z4/b0.Note that because in this paper we focus on the longitudinal performance of the controller, the yaw controller can be regarded as a yaw stabilizer.The input signal of the ADRC is processed by the following procedure:(26)vxb_refvyb_refvzb_ref=cosθ0−sinθ0cosϕ0sinθ0cosθvx_traj−vxvy_traj−vyKzztraj−z+vz_traj−vz,where xtraj, ztraj, vx_traj, and vz_traj are positions and velocities along the trajectory, Kz is the vertical factor, and vxb_ref,vxb_ref,vxb_ref is the input signal of the ADRC controller.The inner loop controller is a PI controller with damping feedbacks, and such controller is proven to be working well with ADRC in trajectory-tracking applications [25].Following the above process, a trajectory controller is obtained. It is needed to mention that when the helicopter is near the ground, the controller is reset to make the helicopter descent slowly. Such techniques to ensure a stable landing is also reported by [12].


## 4. The Tsinghua Rotorcraft Utility Simulation Tools (TRUST)

This section briefly introduces the Tsinghua Rotorcraft Utility Simulation Tools (TRUST), with which we model the S-58 helicopter and validate the trajectories using the ADRC controller. Tsinghua Rotorcraft Utility Simulation Tools (TRUST) stems from [26], which is based on the framework of multibody dynamics. The TRUST simulator is able to simulate many kinds of rotorcrafts, including regular helicopters, tandem helicopters, and compound helicopters, with articulated or flexible rotor hubs. Regular helicopter models consist of main rotor dynamics, tail rotor dynamics, rigid-body fuselage dynamics, horizontal tails, and vertical tails. Such models have proven to be of good agreement with various test data [26].


### 4.1. Model Description

Considering both fidelity and simulation feasibility, we adopt the following modelling methods.


#### 4.1.1. Main Rotor

Main rotor is modelled as an articulated rotor with rigid blades of certain twist. Each blade is attached to the hub through pitch, lag, and flap hinges. The inflow is modelled as the three-state Pitt-Peters dynamic inflow [27]. Blade lift and drag forces are calculated using the blade element theory, and lift-curve slope is obtained from a nonlinear tabulate of angle of attack and relative airflow speed.In addition, because the main rotor may traverse from powered lift state into vortex-ring state (VRS) during autorotation, a modification [28] of vortex-ring state is applied to the three-state Pitt-Peters model.Due to the multibody dynamic framework, each blade’s motion is solved individually by the Newton/Euler equations. Thus, flap angles of blades are not assumed to be small, nor the blades motions are treated as a consequence of a periodical disk, which involves inevitable accuracy loss for blade motion responses [29].


#### 4.1.2. Tail Rotor

Tail rotor is modelled similar to the main rotor, except that the rotor hub is rigid and no VRS correction is applied.


#### 4.1.3. Fuselage

Fuselage is regarded as a rigid body, with linear aerodynamic lift and drag coefficients.


#### 4.1.4. Horizontal and Vertical Tails

The aerodynamic forces of horizontal and vertical tails are calculated by flat plate models.For the purpose of validating autorotation trajectory planning and tracking, the nonlinear helicopter model described above is capable of predicting dynamic responses of the helicopter [28, 30].


### 4.2. Dynamics of TRUST Model

As mentioned in Section2, there are two different helicopter dynamic models used in our study. The nonlinear three-degree-of-freedom rigid-body helicopter dynamics is capable of calculating steady collective pitch manipulations, rotor forces, and power, but not sufficient for dynamic response predictions. Due to the fact that the rotor forces in such models are essentially zero-order processes, dynamic responses of the rotor cannot be well described.In the TRUST helicopter model, a more sophisticated dynamic inflow model is involved, and rotor blade motions are accurately described in the multibody framework. Blade lift and drag forces are from nonlinear tabulates of experiment data. Furthermore, helicopter dynamics is well known for coupling effects between different channels [31]. Hence, for validation purposes, a high-order, six-degree-of-freedom, helicopter dynamic model should be considered.


## 4.1. Model Description

Considering both fidelity and simulation feasibility, we adopt the following modelling methods.


### 4.1.1. Main Rotor

Main rotor is modelled as an articulated rotor with rigid blades of certain twist. Each blade is attached to the hub through pitch, lag, and flap hinges. The inflow is modelled as the three-state Pitt-Peters dynamic inflow [27]. Blade lift and drag forces are calculated using the blade element theory, and lift-curve slope is obtained from a nonlinear tabulate of angle of attack and relative airflow speed.In addition, because the main rotor may traverse from powered lift state into vortex-ring state (VRS) during autorotation, a modification [28] of vortex-ring state is applied to the three-state Pitt-Peters model.Due to the multibody dynamic framework, each blade’s motion is solved individually by the Newton/Euler equations. Thus, flap angles of blades are not assumed to be small, nor the blades motions are treated as a consequence of a periodical disk, which involves inevitable accuracy loss for blade motion responses [29].


### 4.1.2. Tail Rotor

Tail rotor is modelled similar to the main rotor, except that the rotor hub is rigid and no VRS correction is applied.


### 4.1.3. Fuselage

Fuselage is regarded as a rigid body, with linear aerodynamic lift and drag coefficients.


### 4.1.4. Horizontal and Vertical Tails

The aerodynamic forces of horizontal and vertical tails are calculated by flat plate models.For the purpose of validating autorotation trajectory planning and tracking, the nonlinear helicopter model described above is capable of predicting dynamic responses of the helicopter [28, 30].


## 4.1.1. Main Rotor

Main rotor is modelled as an articulated rotor with rigid blades of certain twist. Each blade is attached to the hub through pitch, lag, and flap hinges. The inflow is modelled as the three-state Pitt-Peters dynamic inflow [27]. Blade lift and drag forces are calculated using the blade element theory, and lift-curve slope is obtained from a nonlinear tabulate of angle of attack and relative airflow speed.In addition, because the main rotor may traverse from powered lift state into vortex-ring state (VRS) during autorotation, a modification [28] of vortex-ring state is applied to the three-state Pitt-Peters model.Due to the multibody dynamic framework, each blade’s motion is solved individually by the Newton/Euler equations. Thus, flap angles of blades are not assumed to be small, nor the blades motions are treated as a consequence of a periodical disk, which involves inevitable accuracy loss for blade motion responses [29].


## 4.1.2. Tail Rotor

Tail rotor is modelled similar to the main rotor, except that the rotor hub is rigid and no VRS correction is applied.


## 4.1.3. Fuselage

Fuselage is regarded as a rigid body, with linear aerodynamic lift and drag coefficients.


## 4.1.4. Horizontal and Vertical Tails

The aerodynamic forces of horizontal and vertical tails are calculated by flat plate models.For the purpose of validating autorotation trajectory planning and tracking, the nonlinear helicopter model described above is capable of predicting dynamic responses of the helicopter [28, 30].


## 4.2. Dynamics of TRUST Model

As mentioned in Section2, there are two different helicopter dynamic models used in our study. The nonlinear three-degree-of-freedom rigid-body helicopter dynamics is capable of calculating steady collective pitch manipulations, rotor forces, and power, but not sufficient for dynamic response predictions. Due to the fact that the rotor forces in such models are essentially zero-order processes, dynamic responses of the rotor cannot be well described.In the TRUST helicopter model, a more sophisticated dynamic inflow model is involved, and rotor blade motions are accurately described in the multibody framework. Blade lift and drag forces are from nonlinear tabulates of experiment data. Furthermore, helicopter dynamics is well known for coupling effects between different channels [31]. Hence, for validation purposes, a high-order, six-degree-of-freedom, helicopter dynamic model should be considered.


## 5. Numerical Experiments

In the first part, we show the trajectory planning results of various autorotation initial conditions based on the Sikorsky S-58. In the second part, to validate the dynamic feasibility of the trajectory planning results, we implement the ADRC controller and validate the trajectory planner by making a 6-DOF nonlinear high-fidelity S-58 model track the trajectory.The Sikorsky S-58 is a single-engined helicopter equipped with a 4-bladed articulated rotor with a radius of 8.535 m. The take-off weight used in our research is 4500 kg. Detailed parameters can be found in [32, 33], and main parameters among which can be found in Table 1.Table 1
S-58 specifications.
ParameterValue usedParameterValue usedRotor radius, m8.53Fuselage inertiaIyy, kg m248796Blade chord, m0.417Fuselage inertiaIxx, kg m216265Blade twist, deg-8Fuselage inertiaIzz, kg m248796Number of blades4HT area, m21.15Flap hinge offset, m0.3048Lift-curve slope of HT3.73Rotating inertia of the MR, kg m27101.2Rotor hub height above ground, m4.36Tilting angle, deg0Fuselage flat plate area, m23.55


### 5.1. Autorotation Trajectory Planning Results

The parameters of the cost functions in equations (19) and (20) are set as follows:
(27)W1=50.0,W2=12.5,W3=5.0,W4=1.0,W5=0.5,W6=10.0W7=1.0,W8=1.0WT1=400.0,WT2=0.25,WT3=0.25,WT4=0.2,WT5=0.03.Note that costs for control inputs are added to avoid extreme manipulations. The empirical gliding speed is set asx2glide=5.0. The definition of a successful autorotation landing is defined as
(28)vz≤2m/s,vx≤3m/s,θ≤12deg.In addition, we also implement a start cost function before the value iterations. A start cost function is used to initialize the global value function before the iterations begin. The start cost function can be expressed as follows:(29)Jstart=ψJifz>hpψJ+x1hp−1215x22+15x32+20x42otherwise,where hp is set to 4.0 m in this case. We use the one-way method of N=20,40 points for each dimension and set the max rank approximation of the core functions to rmax=15. Using the one-way method, the initial value function of N=40 is obtained from the solutions of N=20. Convergence is plotted in Figure 3.Figure 3
FT-based algorithms value iteration plots. The blue line denotes values ofN=20, and the black line denotes values of N=40.The left panel shows whenN=40, the value function converges to a certain value around 2000 quickly. The relative error of the value function between different iterations is plotted in the right panel. It is shown that for the discretization of N=40, the relative error is small for a value of about 3×10−5, which is near the converging threshold of 10−5.Various autorotation cases using the trajectory planner are tested, with different initial height and horizontal speed:Case 1.
s0=−55m,20m/s,0m/s,−3.7deg,0deg/s,235RPM.Case 2.
s0=−70m,10m/s,0m/s,−2.4deg,0deg/s,235RPM.Case 3.
s0=−70m,20m/s,0m/s,−3.7deg,0deg/s,235RPM.Case 4.
s0=−80m,10m/s,0m/s,−2.4deg,0deg/s,235RPM.Trajectories are obtained following procedures described in Section2. We apply a Hermite interpolation of 20 points sampled from the preliminary trajectories. Note that such number of points can be regarded as sufficient, considering there are 10 points in [11] and 16 points in [12]. Results are shown in Figure 4.Figure 4
Autorotation trajectories. Black lines for Case1, blue lines for Case 2, green lines for Case 3, and red lines for Case 4.As Figure4 shows, the helicopter enters a steady descent quickly, and begins to decrease at about 20 meters above the ground. Such behaviours agree generally with [9]. For high initial speed cases, the forward speed is decreased immediately after the engine fails, in order to get the helicopter prepared for the final landing. Such behaviours are in accordance with [5]. All the trajectories terminate within the landing constraints defined in (28). Thus, the trajectory planning method using FT-based DP algorithms is able to generate successful autorotation trajectories that satisfy the helicopter autorotation dynamics described in Section 2. Furthermore, comparisons are made with results obtained from Sequential Quadratic Programming (SQP) using the SNOPT software package [34]. In specific, besides the landing specifications assumed in (28), path constraints are required, mainly as 0m/s≤vz≤9m/s,−30deg≤θ≤30deg,Ω¯≥70%. A typical result is demonstrated using initial conditions of Case 1 and is shown in Figure 5.Figure 5
Comparisons of trajectories with results from SQP methods. Black lines are trajectories using FT-based DP method, and blue lines are trajectories using SQP method.Both trajectories in Figure5 lead to safe autorotation landing, considering landing constraints are all satisfied. Besides, the two trajectories show similarities in terms of the rate of descent. The main difference is that the trajectory using the FT-based DP method tends to maintain the rotating speed by entering the steady autorotation right after the engine failure, while the trajectory obtained by SQP method appears to be more agile since the rotor speed does not reach the lower limit of the path constraint. In this respect, the autorotation trajectory using the FT-based DP method is more conservative in terms of keeping a higher rotor speed.Both the above trajectory planning results are obtained using one core of a 3.20 GHz Inter i7-8700 CPU. Time costs for each trajectory generation using the FT-based DP method (shown as Algorithm 1) and SQP method (shown as Algorithm 2) are listed in Table2. The average time cost of the FT-based DP method is 0.60 s, which indicates that the method is real-time feasible.Table 2
Time cost for each trajectory generation.
Case number1234AverageAlgorithm 1 time cost, s0.620.590.610.590.60Algorithm 2 time cost, s9.1211.5411.2113.7411.40Before entering autorotation, the time delay of disengaging the clutch is a key factor that determines whether the subsequent maneuvers lead to a safe landing. A long time delay makes autorotation a tough task because of a low rotor speed. Here, the influence of time delay is considered by setting lowered rotating speeds. Initial conditions are given bys0=−90m,20m/s,0m/s,−3.7deg,0deg/s,Ω, where Ω=100%,95%,90%⋅235RPM for each trajectory.From Figure6, we can see that for lowered rotor speeds caused by different time delays, the trajectory planner is still capable of generating trajectories without modifying any cost function parameter. The rotating speed is recovered along the trajectory by certain maneuvers.Figure 6
Entering autorotation with different rotor speed. Black lines for100%Ω, red lines for 95%Ω, and blue lines for 90%Ω.As for real-time applications, because the trajectory planner is able to generate a global trajectory by any specified initial condition, the computing speed shown in Table2 can be considered real-time applicable with such hardware. The trajectory computing time cost in our study is about 0.6 seconds, which is similar to the time of a fast human pilot reaction. However, as mentioned above, the time delay is a key factor of the autorotation procedure, which should be shortened with best efforts. Therefore, although the trajectory planner is shown to be able to deal with situations of different time delays, an update frequency of at least 1.5 Hz is recommended for real-time applications, considering the time cost of other parts (such as detection of engine failure). Thus, we recommend running such algorithms with CPU of at least 3.1 GHz or higher.


### 5.2. Trajectory-Tracking Simulation with the TRUST 6-DOF Helicopter Dynamic Model

In this part, for validation of the dynamic feasibility of the trajectory, we demonstrate a six-degree-of-freedom flight simulation, with the trajectory planner described in Section2 and the trajectory-tracking controller described in Section 3. For brevity, control parameters can be found in the Appendix. As mentioned before, the simulation is based on a six-degree-of-freedom nonlinear multibody-based helicopter model of S-58, with a VRS correction of the Pitt-Peters dynamic inflow and an accurate dynamic description of the blade’s flap motions.We demonstrate the simulation results of initial conditions from Case1. The initial conditions of states s are given by
(30)s0=−55m,20m/s,0m/s,−0.7deg,0deg/s,235RPM.We assume the engine failure starts att=0s in the forward flight, and the simulations are terminated once the height descends to zero.Simulation results are shown in Figure7. Remember that vzb is positive downward in the body axes and θ is positive in the nose-up direction. The red lines denote the reference values that are needed to be tracked. The black lines denote the 6-DOF nonlinear helicopter model responses. The blue line indicates that the helicopter is close to the ground, and the controller is reset to make a slow landing. The value of the reset height is chosen as 1.0 m by simulation experiments.Figure 7
Height, rotor speed, attitude of Euler angles, and velocities in the body axes during simulation. Black lines denote simulation results, red lines denote reference values of the trajectory-tracking controller, and blue lines denote controller reset for landing.As the results shown in Figure7, the velocities and attitude angles of the helicopter are vxb=−0.14m/s, vyb=−0.42m/s, vz=1.15m/s, ϕ=3.21deg, and θ=−2.33deg when landing, which satisfy the safe autorotation landing specifications defined above. The height history and rotating speed of the main rotor are shown in lower right and the upper right panel of Figure 7, which indicate that the rotating speed is still beyond 75%, although the ground effect is neglected in the six-degree-of-freedom nonlinear helicopter model. Notice that the performance of horizontal speed is not as good as the vertical channel. However, we do not set horizontal distance as any kind of target variable in this study. We also find that increasing horizontal controlling gains affects the vertical channel instead. Apparently, a more sophisticated model-based controller should improve the tracking performance, but our purpose of validating the trajectories is achieved. Simulation results show that following the trajectories generated by the trajectory planner, the helicopter is successfully guided to a safe landing. Thus the trajectory is both real-time feasible and dynamic feasible.


## 5.1. Autorotation Trajectory Planning Results

The parameters of the cost functions in equations (19) and (20) are set as follows:
(27)W1=50.0,W2=12.5,W3=5.0,W4=1.0,W5=0.5,W6=10.0W7=1.0,W8=1.0WT1=400.0,WT2=0.25,WT3=0.25,WT4=0.2,WT5=0.03.Note that costs for control inputs are added to avoid extreme manipulations. The empirical gliding speed is set asx2glide=5.0. The definition of a successful autorotation landing is defined as
(28)vz≤2m/s,vx≤3m/s,θ≤12deg.In addition, we also implement a start cost function before the value iterations. A start cost function is used to initialize the global value function before the iterations begin. The start cost function can be expressed as follows:(29)Jstart=ψJifz>hpψJ+x1hp−1215x22+15x32+20x42otherwise,where hp is set to 4.0 m in this case. We use the one-way method of N=20,40 points for each dimension and set the max rank approximation of the core functions to rmax=15. Using the one-way method, the initial value function of N=40 is obtained from the solutions of N=20. Convergence is plotted in Figure 3.Figure 3
FT-based algorithms value iteration plots. The blue line denotes values ofN=20, and the black line denotes values of N=40.The left panel shows whenN=40, the value function converges to a certain value around 2000 quickly. The relative error of the value function between different iterations is plotted in the right panel. It is shown that for the discretization of N=40, the relative error is small for a value of about 3×10−5, which is near the converging threshold of 10−5.Various autorotation cases using the trajectory planner are tested, with different initial height and horizontal speed:Case 1.
s0=−55m,20m/s,0m/s,−3.7deg,0deg/s,235RPM.Case 2.
s0=−70m,10m/s,0m/s,−2.4deg,0deg/s,235RPM.Case 3.
s0=−70m,20m/s,0m/s,−3.7deg,0deg/s,235RPM.Case 4.
s0=−80m,10m/s,0m/s,−2.4deg,0deg/s,235RPM.Trajectories are obtained following procedures described in Section2. We apply a Hermite interpolation of 20 points sampled from the preliminary trajectories. Note that such number of points can be regarded as sufficient, considering there are 10 points in [11] and 16 points in [12]. Results are shown in Figure 4.Figure 4
Autorotation trajectories. Black lines for Case1, blue lines for Case 2, green lines for Case 3, and red lines for Case 4.As Figure4 shows, the helicopter enters a steady descent quickly, and begins to decrease at about 20 meters above the ground. Such behaviours agree generally with [9]. For high initial speed cases, the forward speed is decreased immediately after the engine fails, in order to get the helicopter prepared for the final landing. Such behaviours are in accordance with [5]. All the trajectories terminate within the landing constraints defined in (28). Thus, the trajectory planning method using FT-based DP algorithms is able to generate successful autorotation trajectories that satisfy the helicopter autorotation dynamics described in Section 2. Furthermore, comparisons are made with results obtained from Sequential Quadratic Programming (SQP) using the SNOPT software package [34]. In specific, besides the landing specifications assumed in (28), path constraints are required, mainly as 0m/s≤vz≤9m/s,−30deg≤θ≤30deg,Ω¯≥70%. A typical result is demonstrated using initial conditions of Case 1 and is shown in Figure 5.Figure 5
Comparisons of trajectories with results from SQP methods. Black lines are trajectories using FT-based DP method, and blue lines are trajectories using SQP method.Both trajectories in Figure5 lead to safe autorotation landing, considering landing constraints are all satisfied. Besides, the two trajectories show similarities in terms of the rate of descent. The main difference is that the trajectory using the FT-based DP method tends to maintain the rotating speed by entering the steady autorotation right after the engine failure, while the trajectory obtained by SQP method appears to be more agile since the rotor speed does not reach the lower limit of the path constraint. In this respect, the autorotation trajectory using the FT-based DP method is more conservative in terms of keeping a higher rotor speed.Both the above trajectory planning results are obtained using one core of a 3.20 GHz Inter i7-8700 CPU. Time costs for each trajectory generation using the FT-based DP method (shown as Algorithm 1) and SQP method (shown as Algorithm 2) are listed in Table2. The average time cost of the FT-based DP method is 0.60 s, which indicates that the method is real-time feasible.Table 2
Time cost for each trajectory generation.
Case number1234AverageAlgorithm 1 time cost, s0.620.590.610.590.60Algorithm 2 time cost, s9.1211.5411.2113.7411.40Before entering autorotation, the time delay of disengaging the clutch is a key factor that determines whether the subsequent maneuvers lead to a safe landing. A long time delay makes autorotation a tough task because of a low rotor speed. Here, the influence of time delay is considered by setting lowered rotating speeds. Initial conditions are given bys0=−90m,20m/s,0m/s,−3.7deg,0deg/s,Ω, where Ω=100%,95%,90%⋅235RPM for each trajectory.From Figure6, we can see that for lowered rotor speeds caused by different time delays, the trajectory planner is still capable of generating trajectories without modifying any cost function parameter. The rotating speed is recovered along the trajectory by certain maneuvers.Figure 6
Entering autorotation with different rotor speed. Black lines for100%Ω, red lines for 95%Ω, and blue lines for 90%Ω.As for real-time applications, because the trajectory planner is able to generate a global trajectory by any specified initial condition, the computing speed shown in Table2 can be considered real-time applicable with such hardware. The trajectory computing time cost in our study is about 0.6 seconds, which is similar to the time of a fast human pilot reaction. However, as mentioned above, the time delay is a key factor of the autorotation procedure, which should be shortened with best efforts. Therefore, although the trajectory planner is shown to be able to deal with situations of different time delays, an update frequency of at least 1.5 Hz is recommended for real-time applications, considering the time cost of other parts (such as detection of engine failure). Thus, we recommend running such algorithms with CPU of at least 3.1 GHz or higher.


## 5.2. Trajectory-Tracking Simulation with the TRUST 6-DOF Helicopter Dynamic Model

In this part, for validation of the dynamic feasibility of the trajectory, we demonstrate a six-degree-of-freedom flight simulation, with the trajectory planner described in Section2 and the trajectory-tracking controller described in Section 3. For brevity, control parameters can be found in the Appendix. As mentioned before, the simulation is based on a six-degree-of-freedom nonlinear multibody-based helicopter model of S-58, with a VRS correction of the Pitt-Peters dynamic inflow and an accurate dynamic description of the blade’s flap motions.We demonstrate the simulation results of initial conditions from Case1. The initial conditions of states s are given by
(30)s0=−55m,20m/s,0m/s,−0.7deg,0deg/s,235RPM.We assume the engine failure starts att=0s in the forward flight, and the simulations are terminated once the height descends to zero.Simulation results are shown in Figure7. Remember that vzb is positive downward in the body axes and θ is positive in the nose-up direction. The red lines denote the reference values that are needed to be tracked. The black lines denote the 6-DOF nonlinear helicopter model responses. The blue line indicates that the helicopter is close to the ground, and the controller is reset to make a slow landing. The value of the reset height is chosen as 1.0 m by simulation experiments.Figure 7
Height, rotor speed, attitude of Euler angles, and velocities in the body axes during simulation. Black lines denote simulation results, red lines denote reference values of the trajectory-tracking controller, and blue lines denote controller reset for landing.As the results shown in Figure7, the velocities and attitude angles of the helicopter are vxb=−0.14m/s, vyb=−0.42m/s, vz=1.15m/s, ϕ=3.21deg, and θ=−2.33deg when landing, which satisfy the safe autorotation landing specifications defined above. The height history and rotating speed of the main rotor are shown in lower right and the upper right panel of Figure 7, which indicate that the rotating speed is still beyond 75%, although the ground effect is neglected in the six-degree-of-freedom nonlinear helicopter model. Notice that the performance of horizontal speed is not as good as the vertical channel. However, we do not set horizontal distance as any kind of target variable in this study. We also find that increasing horizontal controlling gains affects the vertical channel instead. Apparently, a more sophisticated model-based controller should improve the tracking performance, but our purpose of validating the trajectories is achieved. Simulation results show that following the trajectories generated by the trajectory planner, the helicopter is successfully guided to a safe landing. Thus the trajectory is both real-time feasible and dynamic feasible.


## 6. Conclusions

This paper demonstrates a trajectory planning method for autorotation, which is real-time feasible and guarantees a strict satisfaction of specified helicopter dynamics along the trajectory. A trajectory-tracking controller is demonstrated to ensure the helicopter fly along the trajectories. Successful autorotation trajectories with various initial conditions are shown and discussed, and the time cost of trajectory generation procedures is around 0.60 s. A comprehensive validation is made based on a six-degree-of-freedom high-fidelity nonlinear S-58 model using the TRUST simulator. Simulation results show that although the horizontal velocities are not tracked as well as the vertical channel, the trajectories generated by the trajectory planner and the trajectory tracker are capable of guiding a helicopter into a successful autorotation. Thus, the trajectory planner is real-time feasible, and the generated trajectories are dynamic feasible. Further potential improvements for the trajectory planning method could be selecting the derivatives of the control as control inputs of the three-degree-of-freedom helicopter dynamics, extending the two-dimensional helicopter dynamics to a three-dimensional helicopter dynamic model, and automatic tuning of the weight functions.



---

*Source: 1372052-2020-07-05.xml*