Modelling and robust control of an unmanned coaxial rotor helicopter with unstructured uncertainties

Abstract

A complete methodology for an unmanned coaxial rotor helicopter with unstructured uncertainties was proposed to achieve high-accuracy tracking performance from modelling to robust control. An integrative approach was introduced to systematically construct a whole dynamic model. The key parameters were selected carefully after iteratively being checked by empirical coefficients to decrease the budget and risk of programme. Moreover, a new control scheme is proposed to simultaneously incorporate six inputs to control six states based on the investment of singularity value responses and the general rule of relative gain array. Coprime factor uncertainty model is considered to represent a class of unstructured uncertainties, such as unmolded actuator dynamics and unpredicted interferences between two rotors. Furthermore, the $H_{\infty }$ loop-shaping control was proposed to apply the control design of the coaxial rotor helicopter to manage complicated uncertainties and multivariable coupling. Finally, simulation results show the effectiveness of the proposed controller design in the step response of the closed loop. The stable closed-loop plant is achieved and the tolerant size of unstructured uncertainty is up to 36.09%. Good step responses and satisfied decoupling were also investigated in detail.

Introduction

Interest in coaxial rotor helicopter has been increasing rapidly in recent years.^1–3 The coaxial rotor is a configuration that adopted two coaxially counter-rotating rotors. Compared with the traditional single rotor, coaxial rotor consumes less induced power while exhibits the same blade properties and solidity characteristics under different flight conditions. ⁴ This configuration has potential to achieve high speed, long range and heavy lift. This type of compact design also would increase efficiency and conserves space. However, a coaxial rotor helicopter is severely nonlinearly dynamic, has a poor model fidelity and strong coupling plant. In addition, there is a reciprocal and unpredicted aerodynamic effect between the upper and lower rotors. These characteristics considerably complicate the modelling and control of coaxial rotor helicopter to achieve high-performance design. Therefore, it will be a challenging work to model and robust control an unmanned coaxial rotor helicopter with uncertainties.

A comprehensive survey on coaxial rotor aerodynamics has been concluded before 1997 in most countries dealing with wake structure, rotor spacing and load sharing between the upper and lower rotors from the perspectives of computations and experiments. ⁵ A series of research works have analysed the coaxial rotor performance affected by space, blade twist and rotor rotation.^4,6 The individual blade geometries were investigated to seek the highest levels of efficiency utilizing an optimization method.^2,7 Moreover, aeromechanics optimum design for the coaxial compound helicopter was also taken into account.^1,3 In these previous works, however, researchers have mostly focused on performance analysis and optimum rotor. It lacks systematical modelling of the whole coaxial rotor helicopter. Therefore, modelling a coaxial rotor helicopter is needed to further promote.

Robust control law is generally considered an attractive approach to deal with inaccuracies and uncertainties. This method has been implemented to extensive successful actual flight applications, such as a multivariable controller in a vertical and short take-off landing aircraft, ⁸ an unmanned helicopter design ⁹ and the Bell 205. ¹⁰ In addition, Lin et al. ¹¹ introduced $H_{\infty }$ method to flight control of fixed-pitch coaxial rotor helicopter. A mixed sensitivity $H_{\infty }$ controllers have been developed from test benches to successful flights for a coaxial micro helicopter. ¹² Adaptive fuzzy tracking control scheme for marine vehicles with unknown parametric dynamics and uncertainties. ¹³ A control scheme combining the radial basis function neural networks and wave variable technique was proposed to simultaneously compensate for the effects caused by communication delays and dynamic uncertainties. ¹⁴ The neural networks approximation technique is employed to compensate for uncertainties in the robot arms and manipulated object. ¹⁵ Automatic motion control is investigated for wheeled inverted pendulum models with the unknown dynamics of subsystem. ¹⁶ However, these applications did not seriously consider the effect of uncertainties and nonlinear dynamic on performance. There always exist some complicated uncertainties in the coaxial rotor helicopter model, such as unmodelled actuator dynamics, sensor errors, high harmonic flapping dynamics, external disturbances and unpredicted interferences between the rotors. Uncertainties may lead to lose of tracking accuracy considerably, even leading to unstable condition and flight accidents. Therefore, it is starving for a control law to manage its complicated natural uncertainties and multivariable coupling of coaxial rotor helicopter, systematically.

Therefore, a complete methodology for an unmanned coaxial rotor helicopter with unstructured uncertainties was proposed to achieve high-accuracy tracking performance from modelling to robust control. An integrative approach was introduced to systematically construct a whole dynamic model. $H_{\infty }$ loop-shaping control was proposed to apply the control design of the coaxial rotor helicopter to manage complicated uncertainties and multivariable coupling. Furthermore, a novel control scheme was proposed including reasonable input–output pairing to simultaneously incorporate six inputs to control six states. Coprime factor uncertainty model is considered to represent a class of unstructured uncertainties. A mathematical prototype of coaxial rotor helicopter was constructed under hovering flight condition. Key parameters were identified by experiments or empirical formulations from the open literature. Finally, simulation results show the effectiveness of the proposed controller design in the step response of closed loop.

Modelling and trimming

Modelling

The dynamic model with a 6-degree-of-freedom (6-DOF) rigid body has appeared effective for controller design as well as the generation mechanisms of aerodynamic force and torque in the Bell 205 flight test. ¹⁰ This method was also considered as a reasonable representation of the key characteristics of a helicopter.^17,18 Therefore, this integrative approach is adopted to construct the dynamic rigid body combining with the sum of the forces and torques of fundamental component model. The basic 6-DOF rigid body is fixed at the centre of gravity $O_e$ . And body coordinates are shown in Figure 1.

Figure 1.

Body coordinates in the prototype.

This body-fixed axis has three orthogonal axes, namely, x-axis forwards, y-axis rightwards and z-axis downwards. The applied aerodynamic forces and moments of various components were transformed to act along the body-fixed axis. The Newton–Euler equations are effectively applied in the form of nonlinear differential equations. The forces and moments are contributed by the coaxial rotor, fuselage, vertical fin and horizontal empennage.

Forces and moments generated by the each component are considered in detail. The coaxial rotor is the key component in constructing a simulation model. Because of the complexity and reciprocal effects between the rotors, calculating the thrusts and torques of the coaxial rotor is a complex and difficult work. However, in predicting the thrust and torque coefficients, blade element momentum theory (BEMT) has been validated to nearly consistent with the free-vortex prediction method ¹⁹ which offers a good approximation of the experimental data using full-scale and model-scale coaxial rotors. ²⁰ In addition, it is assumed that the slipstream wake of upper rotor interferes with the inflow of the lower rotor but the lower rotor does not affect the inflow of upper rotor. The advancing blade of the upper rotor was on the right side and the lower rotor rotated clock-wise. Therefore, an efficient BEMT is addressed to calculate thrust and torque coefficients of the rotors in the following way

$$C_{T_i}=\underset{r=0}{\overset{r=1}{\int }}d{C}_{T_i}$$

(1)

$$C_{Q_i}=\underset{r=0}{\overset{r=1}{\int }}{\lambda }_{i}d{C}_{T_i}+\underset{r=0}{\overset{r=1}{\int }}\left(\frac{1}{2}\right){\sigma }_i{C}_d{r}^3\mathrm{dr}$$

(2)

where $i=\{0,1\}$ stands for up and low rotor.

Therefore, the inflow velocity $\lambda$ and incremental thrust of the upper rotor are derived as follows

$${\lambda }_{\mathrm{up}}=\sqrt{{\left(\frac{{\sigma }_{\mathrm{up}}{C}_{l_{\alpha }}}{16F}\right)}^2+\frac{{\sigma }_{\mathrm{up}}{C}_{l_{\alpha }}}{8F}{\theta }_{\mathrm{up}}r}-\left(\frac{{\sigma }_{\mathrm{up}}{C}_{l_{\alpha }}}{16F}\right)$$

(3)

$$d{C}_{T_{\mathrm{up}}}=\left(\frac{{\sigma }_{\mathrm{up}}{C}_{l_{\alpha }}}{2}\right)({\theta }_{\mathrm{up}}{r}^2-{\lambda }_{\mathrm{up}}r)\mathrm{dr}$$

(4)

Blade is assumed to discretize into a series of small radial constant elements. ¹⁹ $\sigma$ is the solidity of rotor. $C_{l_{\alpha }}$ is lift curve slope. F is Prandtl’s tip loss function. $\theta$ is blade pitch. r is the non-dimensional radial span position. Thrust and torque coefficients of the upper rotor can be numerically integrated into each respective rotor disc

$${\lambda }_{\mathrm{up}}\left(r_n\right)=\frac{{\sigma }_{\mathrm{up}}{C}_{\mathrm{la}}}{16}\sqrt{1+\frac{32}{{\sigma }_{\mathrm{up}}{C}_{l_{\alpha }}F}{\theta }_{\mathrm{up}}\left(r_n\right)r_n}-\left(\frac{{\sigma }_{\mathrm{up}}{C}_{l_{\alpha }}}{16F}\right)$$

(5)

$$\mathrm{\Delta }{C}_{T_n}=\left(\frac{\sigma C_{l_{\alpha }}}{2}\right)(\theta \left(r_n\right)r_n^2-\lambda \left(r_n\right)r_n)\mathrm{\Delta }r$$

(6)

$$C_T=\sum _{n=1}^N\mathrm{\Delta }{C}_{T_n}$$

(7)

$$C_Q=\sum _{n=1}^N{\lambda }_{\mathrm{up}}\left(r_n\right)\mathrm{\Delta }{C}_{T_n}+\left(\frac{1}{2}\right){\sigma }_{\mathrm{up}}\sum _{n=1}^N{C}_d{r}_n^3\mathrm{\Delta }r$$

(8)

The mathematical principle of the upper rotor was also adopted to calculate thrust and torque coefficients of the lower rotor. In addition, the inflow region of the lower rotor is reasonably divided into two parts. The first one is called the inner part that is fully encountered by the slipstream wake of the upper rotor. The second one is the outer part that is affected freely by the slipstream wake. The boundary of two parts is defined by the ideal radial contraction of the wake of the upper rotor. Consequently, the inflow velocity $\lambda$ of the lower rotor is as follows.

• For $r<a$

$${\lambda }_{\mathrm{low}}=\sqrt{{(\frac{{\sigma }_{\mathrm{low}}{C}_{l_{\alpha }}}{16F}-\frac{{\lambda }_{\mathrm{up}}}{2a^2})}^2+\frac{{\sigma }_{\mathrm{low}}{C}_{l_{\alpha }}}{8F}{\theta }_{\mathrm{low}}r}-(\frac{{\sigma }_{\mathrm{up}}{C}_{l_{\alpha }}}{16F}-\frac{{\lambda }_{\mathrm{up}}}{2a^2})$$

(9)

• For $r\ge a$

$${\lambda }_{\mathrm{low}}=\sqrt{{\left(\frac{{\sigma }_{\mathrm{low}}{C}_{l_{\alpha }}}{16F}\right)}^2+\frac{{\sigma }_{\mathrm{low}}{C}_{l_{\alpha }}}{8F}{\theta }_{\mathrm{low}}r}-\left(\frac{{\sigma }_{\mathrm{low}}{C}_{l_{\alpha }}}{16F}\right)$$

(10)

The a factor represents the maximum radial station of the wake of upper rotor. The thrust and power coefficient of the lower rotor can be calculated through the numerical integration as the same as the mathematical principles of upper rotor.

Therefore, the effective blade element theory is applied to calculate the thrust and torque of the upper and lower rotors over the blade span

$$\left[\begin{array}{c}{T}_i=\rho A{\left(\mathrm{\Omega }R\right)}^2{C}_{T_i}\hfill \\ \multicolumn{1}{c}{Q_i=\rho A{\left(\mathrm{\Omega }R\right)}^{2}R{C}_{Q_i}}\end{array}\right]$$

(11)

where T is the thrust, Q stands for torque, $\rho$ is air density, A is the disc area, $\mathrm{\Omega }$ is the rotational speed and R is radius of the rotor. $C_T$ and $C_Q$ are thrust and torque coefficients.

Rotor flapping dynamics is a primary factor in the process of control design. Thus, flapping equation for each rotor is described by

$${\stackrel{··}{\beta }}_i=\frac{1}{J_{{\beta }_i}}\underset{0}{\overset{R-e}{\int }}{\left(M_h\right)}_i\mathrm{dr}$$

(12)

where $\beta$ is out-of-plane flapping angle. $J_{\beta }$ is the moment of inertia of blade, R is the rotor radius, e is the blade hinge offset from the centre of the rotor, $M_h$ is the incremental moment coefficient with respect to the flapping hinge along the blade.

The solution of flapping equation was assumed to be periodic. Therefore, first harmonic terms formation of flapping dynamics was sufficiently derived under the quasi-static assumption in the following way

$${\beta }_i(t)=a_{i0}(t)-a_{i1}(t)\cos \Psi +b_{i1}(t)\sin \Psi$$

(13)

The three terms represent the coning, the pitching and the rolling motion of the tip-path plane, respectively. The first harmonic coefficients relationship between the flapping motion and blade pitching motion is represented as follows

$$\left\{\begin{array}{c}{a}_{i0}=\frac{\gamma }{8{\lambda }_{\beta }^2}({\theta }_{i0}-\frac{4}{3}{\lambda }_i)\hfill \\ \multicolumn{1}{c}{a_{i1}=\frac{1}{1+S_{\beta }^2}\{S_{\beta }{\theta }_{i1c}-{\theta }_{i1s}+(S_{\beta }\frac{16}{\gamma }-1)\overline{p}+(S_{\beta }+\frac{16}{\gamma })\overline{q}\}}\\ \multicolumn{1}{c}{b_{i1}=\frac{1}{1+S_{\beta }^2}\{S_{\beta }{\theta }_{i1s}+{\theta }_{i1c}-(S_{\beta }\frac{16}{\gamma }-1)\overline{q}+(S_{\beta }+\frac{16}{\gamma })\overline{p}\}}\end{array}\right\}$$

(14)

where $\gamma$ is Lock number, ${\lambda }_{\beta }=1+(K_{\beta }/I_{\beta }{\mathrm{\Omega }}^2)$ is the flapping frequency ratio, ${\theta }_0$ is collective pitch, ${\lambda }_i$ is the inflow velocity of rotor, $S_{\beta }=8({\lambda }_{\beta }^2-1)/\gamma$ is stiffness number, ${\theta }_{1c}$ and ${\theta }_{1s}$ are first harmonic coefficients of pitch motion, $\overline{p}=p/\mathrm{\Omega }$ and $\overline{q}=q/\mathrm{\Omega }$ .

Furthermore, the idea of modelling conventional helicopter’s fuselage was expanded to model fuselage of coaxial rotor helicopter. It was supposed that the force and moment of the fuselage only affected by downwash of lower rotor. Therefore, an empirical formulation was adopted to consider the effects of the lower rotor downwash in taking use of the isolated fuselage shape. ²¹ Thus, the forces of fuselage are calculated as follows

$$D_F=0.3\left(\frac{\mathrm{Projected\; area}}{\mathrm{Disc\; area}}\right)\times G.W$$

(15)

Thus, the forces and moments of the vertical fin and horizontal empennage were considered to be the same as in the literature. ²²

Consequently, the forces and moments of each component were summed up in body-axis system to construct the simulation model. The fidelity level of the constructed model is positioned between levels 1 and 2 that enables designing of the controller. This level of fidelity can interpret the cause and effect of model behaviour.

Parameter design and evaluation

Reliable design parameters of the coaxial rotor helicopter are needed to satisfy fundamental aerodynamic constraints. Empirically verified processes are comprehensively and systematically shown to decrease budget and risk in preliminary design stages in Figure 2.

Figure 2.

Process of the parameter design.

Gross weight requirement of the coaxial helicopter is 200 kg. Empirical value ranges of the parameter, such as thrust coefficients, tip speeds, solidity and figure of merit (FOM), are selected as the constrained conditions. These parameters reflect the key and fundamental aerodynamic constraints in structure design of the helicopter.

After iteratively being verified by empirical parameter, the reasonable parameters of the prototype are determined, which are shown in Table 1. The structural parameters are shown in Table 2. The configurations of the upper and lower rotors are identical in teetering two-blade including the NACA0012 blade sections, radius, platforms, twisted rates and cut-out.

Table 1.

Coaxial-helicopter properties.

Property	Value	Unit
Mass	200	kg
Diameter	3.6	m
Rotor space	0.18	m
Chord	0.15	m
Twist	−6	°
Off-set	0.1	m
Rotor speed	800	r/min
Lock number	4
Ixx	293.96	kg m2
Iyy	915.34	kg m2
Izz	881.4	kg m2
Ixz	34.37	kg m2

Table 2.

Structural parameters.

symbol	Physical meaning	Value (m)
$h_M$	Distance with respect to the origin of the body-fixed axes in z-axis	0.3
$y_M$	Distance with respect to the origin of the body-fixed axes in y-axis	0.05
$l_M$	Distance with respect to the origin of the body-fixed axes in x-axis	−0.05
$l_F$	Horizontal distance between the centres of fuselage and body	−0.1
$l_H$ , $l_V$	$X_B$ distance between the origin of fin and empennage	−0.2
$F_L,F_W,F_H$	Length, width and height of fuselage	3, 0.35, 0.5
$b_i$	Blade number	2, 2

In order to show the reasonable and satisfied requirements of the selected parameters, the verified calculation is implemented using the parameters in Tables 1 and 2 as follows:

Thrust coefficient

It is a non-dimensional coefficient to conveniently calculate the rotor thrust. Empirical value range of $C_T$ is 0.0036–0.0144. ²³ Taking equations (5)–(7), thrust coefficients of upper and lower rotors are, respectively, 0.0046 and 0.0035, which satisfy the empirical value range.

Blade solidity

It stands for the ratio of blade area to disc area as follows

$$\sigma =\frac{\mathrm{bc}}{\pi R}$$

(16)

Empirical value range of $\sigma$ is 0.05–0.12, ²¹ coaxial rotor was 0.106, and that of upper and lower rotors was 0.053. They all belong to empirical value range.

Disc loading

It is defined as thrust per unit area of disc

$$\mathrm{DL}=\frac{T}{\pi R^2}$$

(17)

Empirical value range of DL is 191.394–574.182 N/m². ²¹ Disc area of one of the rotors was considered. The calculated value is 198.8 N/m² which belongs to the range.

Tip speed

It is used to avoid tip stall and compressibility. The formulation is represented as follows

$$V_{\mathrm{tip}}=\mathrm{\Omega }·R$$

(18)

where $\mathrm{\Omega }=2\pi n/60$ . Empirical value range of tip speed is 150.4–243.84 m/s. ²³ The value of angular velocity $\mathrm{\Omega }$ of the rotor is selected as 83.7 rad/s. Then, the value of tip speed is 150.72 m/s and Mach is 0.49. It means that this lower tip speed has the lower noise and better hovering performance.

Lock number

It represents the ratio of aerodynamic to centrifugal forces. The heavier the blade, the lower the Lock number ²¹

$$\gamma =\frac{\rho \mathrm{ac}{R}^4}{I_b}$$

(19)

$$I_b=\left(\frac{m{R}^3}{3}\right){(1-\left(\frac{e}{R}\right))}^3$$

(20)

Empirical value range of Lock number is 2–10. Carbon fibre was selected as the material of blade. And its density is 1.75 g/cm³. Then, the mass of the single blade is 1.785 kg when the value of thickness of the blade is designed as 0.4 cm. Therefore, the value of Lock number is 4 using equations (19) and (20).

FOM

It is defined as the ratio of induced power to actual power as follows ²⁴

$$\mathrm{FOM}=\frac{C_T^{3/2}}{\sqrt{2}{C}_Q}$$

(21)

$$\mathrm{FOM}\_\mathrm{coaxial}=\frac{{\kappa }_{\mathrm{int}}(C_{T_{\mathrm{up}}}^{3/2}+C_{T_{\mathrm{up}}}^{3/2})}{\sqrt{2}{C}_Q}$$

(22)

Empirical value range of FOM is 0.75–0.80. ²¹ Using equation (21), the values of upper and lower rotors are, respectively, 0.7316 and 0.4961. The value of FOM is 0.7770 through equations (22). The other two parameters belong to empirical value range.

Power calculation and engine selection

Power calculation is used to select a suitable engine. The total consumed power of the coaxial rotor is calculate by

$$P_{\mathrm{up}}=C_p\rho A{\left(\mathrm{\Omega }R\right)}^3$$

(23)

$$P_{\mathrm{rotor}}=P_{\mathrm{up}}+P_{\mathrm{low}}$$

(24)

Then, the total power of the coaxial rotor is 26.4 kW. After survey of the aircraft engine, Rotax 582 is a perfect choice to satisfy the power requirement.

Gross weight requirement

This requirement is used to guarantee the coaxial rotor helicopter to get off the ground. The trim equation in z-axis is as follows

$$-(T_{\mathrm{up}}+T_{\mathrm{low}})+G.W+D_F=0$$

(25)

To calculate equations (7) and (11), thrusts of upper and low rotors are 1141.5 and 881 N, respectively. The drag of fuselage is 60.5 N. Therefore, the gross weight requirement of the coaxial helicopter with 200 kg is satisfied.

Consequently, the e-mathematical calculations show that the parameters are reasonable and sufficient. In addition, the reliable design parameters will effectively decrease the budget and risk of programme (Table 3).

Table 3.

Values of fundamental aerodynamic parameters.

Parameter	Value
	Upper rotor	Lower rotor
Thrust coefficient	0.0046	0.0035
Blade solidity	0.053	0.053
Disc loading	198.8 N/m2	198.8 N/m2
Tip speed	150.72 m/s	150.72 m/s
Lock number	4	4
Figure of merit	0.7316	0.4961
Power calculation	13.2 kW	13.2 kW
Shared gross weight	1141.5 N	881 N

Trimming

Trimming is a significant procedure to equilibrate the nonlinear dynamic equation under the specified flight condition. To establish a square equilibrium equation, the same cyclic pitch inputs for the upper and lower rotors were implemented to trim the nonlinear equation of the coaxial helicopter. ²⁵

Based on the small angle assumption, the linear equilibrium equations are obtained. Then, these equations are solved using analytical method. The trim values are listed in Table 4. Table 5 outlines the thrust and state variables of the coaxial rotor helicopter in hovering flight.

Table 4.

Trim values.

$(B_1{)}_{\mathrm{up}}$	$(A_1{)}_{\mathrm{up}}$	$(B_1{)}_{\mathrm{low}}$	$(A_1{)}_{\mathrm{low}}$	$({\theta }_{\mathrm{col}}{)}_{\mathrm{up}}$	$({\theta }_{\mathrm{col}}{)}_{\mathrm{low}}$
−0.1138	0.1211	1.1138	0.1211	0.242	0.2731

Table 5.

State variables in hovering flight.

$T_{\mathrm{up}}$	$T_{\mathrm{low}}$	$\theta$	$\varphi$	p	q	r
1141.5	881	0.1147	−0.1249	0	0	0

In Table 4, the collective pitch of the lower rotor is slightly higher than the upper rotor. This is because the part of lower rotor’s collective pitch compensates for the larger downwash wake interference from the upper rotor. It is found that lateral cyclic trim values are larger than longitudinal cyclic pitch values to balance the equilibrium equation in the constructed model.

Another interesting consequence is that the upper rotor shared 56.4% of the total thrust under torque-balanced condition (Table 5). This is the reason that some of induced power of lower rotor is lost. In addition, the thrusts ratio exhibits a good agreement with other unrestricted experimental results. ²⁰ Consequently, the trim and state values show that the constructed model is an effective representation of the main characteristics of the coaxial rotor in hovering condition.

Dynamic stability and singularity

Dynamic stability characteristics of a coaxial rotor helicopter absolutely restrict the goal of control design. Small perturbation theory was adopted to calculate the stable and aerodynamic derivatives of the coaxial rotor helicopter at hovering condition. Furthermore, the superposition of individual effects of the upper and lower rotors is supposed to depict changes in the dynamic system of the constructed model. Effects of the coaxial rotors are linearly added by assuming the helicopter with two main rotors. Therefore, the linear model is derived as follows

$$\stackrel{·}{\mathbf{x}}=\mathbf{Ax}+\mathbf{Bu}$$

(26)

where the state and control matrices were written in terms of $\mathbf{A}$ and $\mathbf{B}$ . The elements of the state and control matrices are the sum of the effects of the upper and lower rotors, respectively. The input vector is

$$\mathbf{u}=\left[\begin{array}{c}{(A_1,B_1,{\theta }_{\mathrm{col}})}_{\mathrm{up}},{(A_1,B_1,{\theta }_{\mathrm{col}})}_{\mathrm{low}}\hfill \end{array}\right]$$

(27)

The state vector is

$$\mathbf{x}=\left[uw{\theta }_{q}v{\varphi }_{p}r\right]$$

(28)

The natural modes of motion are indicated by the eigenvalues of the state matrix as shown in Table 6. Pitch subsidence, roll subsidence and heave subsidence modes of the constructed model are stable. However, the model is unstable in the yaw mode, which most likely lacks a tail rotor. Two complex conjugate pairs on the right half plane exhibit unstable oscillations of longitudinal and lateral dynamics. The first conjugate pair attributes to the coupling of pitch and longitudinal velocity by speed stability. The second is caused by the coupling of roll and lateral velocity. These two results are also similar to the phenomenon in a conventional helicopter. The distribution trend of natural modes of the constructed model also coincides well with the results of other research. ¹ Consequently, the constructed model is unstable without an automatic augmentation system.

Table 6.

Eigenvalues of the model.

Eigenvalues	Natural modes
−0.563	Roll
−0.4806	Pitch
−0.4368	Heave
0.2125 ± 0.3908i	Unstable
0.6308 ± 0.4509i	Unstable
0.6925	Unstable

After calculating the state and control matrices, it is necessary to consider control scheme as diagonally dominant as possible with input–output pairing. Time-domain method is not efficient to analyse the relationships between inputs and outputs because of divergences of the unstable model response. Therefore, singularity value response was implemented to contribute to the controller structure design in frequency domain.

Figure 3 indicates that considerable couplings are observed between longitudinal velocity and mixed-rate pitch. The differences of responses that longitudinal velocity and mixed-rate pitch with respect to the longitudinal cyclic pitch of upper and lower rotor are minimal. Moreover, the longitudinal cyclic $(B_1{)}_{\mathrm{up}}$ is more effective on ${\theta }_q$ than other state variables except u and v, as shown in Figure 3(a). Then, it indicates that $(B_1{)}_{\mathrm{low}}$ can be used to control u, as shown in Figure 3(b).

Figure 3.

The longitudinal cyclic to states: (a) upper rotor and (b) lower rotor.

It shows that there are considerable couplings among the controlled variables with respect to lateral cyclic input, as shown in Figure 4. The difference between lateral velocity and mixed-rate roll is small. The longitudinal velocity u and lateral velocity v are clearly affected by the longitudinal $(B_1{)}_{\mathrm{low}}$ and $(B_1{)}_{\mathrm{up}}$ as well as the lateral $(B_1{)}_{\mathrm{up}}$ and $(A_1{)}_{\mathrm{up}}$ . In addition, Figure 4(a) shows that $(A_1{)}_{\mathrm{up}}$ is more effective on ${\varphi }_p$ than ${\theta }_q$ . The lateral cyclic pitch is a little bit more effective on v than u, as shown in Figure 4(b).

Figure 4.

The lateral cyclic to states: (a) upper rotor and (b) lower rotor.

As shown in Figures 3 and 4, inherent cross couplings are observed between longitudinal and lateral motions. These findings indicate that the constructed dynamic model has an ability to depict the basic axis symmetry characteristic of a coaxial rotor helicopter.

Serious couplings between heave and yaw motions are observed (Figure 5). The couplings between longitudinal-lateral dynamic and heave-direction dynamic also exist. However, the collective pitches of two rotors are more effective on vertical velocity w and yaw rate r than the other controlled variables. Moreover, the collective pitch of the two rotors significantly affects heave motion than yaw motion, as shown in Figure 5. Therefore, one of the two collective pitch can control vertical velocity. The remaining input could be used to control yaw rate.

Figure 5.

Collective pitch to states: (a) upper rotor collective and (b) lower rotor collective.

In summary, the singularity value response simulations show that the constructed model experiences complicated couplings among channels, as shown in Figures 3–5. The relationship between input and output is unclear. Some situations are even confusing. In order to contribute to make the control scheme as diagonally dominant as possible, it is needed to investigate the characteristics of the constructed model further.

Robust control design

Control theory

$H_{\infty }$ loop shaping exhibits following advantages: (1) no iteration to achieve optimal $\gamma$ ; (2) no pole-zero cancellations and (3) robust to instabilities, uncertainties and perturbations. ²⁶ Therefore, $H_{\infty }$ loop shaping was proposed to apply to the robust control design of coaxial rotor helicopter. The main idea of the method is that it minimizes the infinity norm of the transfer function matrix between exogenous signal and objective signal vectors

$$\underset{\mathrm{stab}{K}_{\infty }}{max}\Vert \left[\begin{array}{c}{w}_1\hfill \\ \multicolumn{1}{c}{w_2}\end{array}\right]\to \left[\begin{array}{c}{z}_1\hfill \\ \multicolumn{1}{c}{z_2}\end{array}\right]\Vert =\epsilon$$

(29)

The nominal plant G is shown in state-space form

$$G:=\left[\begin{array}{cc}A\hfill & B\hfill \\ \multicolumn{1}{c}{C}& D\hfill \end{array}\right]$$

(30)

The shaped plant $G_s$ is achieved through multiplying two weighting matrices

$$G_s=W_{2}G{W}_1$$

(31)

Coprime factor uncertainty model is adopted to consider a class of unstructured uncertainties. Then, the perturbation model of normalized factorization formula of plant is described as follows

$$G_p={(M+{\mathrm{\Delta }}_M)}^{-1}(N+{\mathrm{\Delta }}_N)$$

(32)

where ${\mathrm{\Delta }}_M$ and ${\mathrm{\Delta }}_N$ are a class of uncertainties transfer function. They stand for a family of unpredicted or unconsidered uncertainties in additive and multiplicative uncertainty formulation bounded by $H_{\infty }$ norm. Then, stability margin is

$${\epsilon }_{max}^{-1}={\{1-\Vert \left[NM\right]{\Vert }_{\mathrm{H}}^2\}}^{-1/2}={(1+\rho (\mathrm{XZ}))}^{1/2}$$

(33)

where $\epsilon$ is the stability margin, $||•|{|}_{\mathrm{H}}$ is the Hankel norm, $\rho$ is the maximum eigenvalue, X and Z are unique positive solutions to algebraic Riccati equation

$$\begin{array}{c}(A-B{F}^{-1}{D}^{T}C)Z+Z{(A-B{F}^{-1}{D}^{T}C)}^T\hfill \\ -Z{C}^T{E}^{-1}\mathrm{CZ}+B{F}^{-1}{B}^T=0\hfill \end{array}$$

(34)

$$\begin{array}{c}{(A-B{F}^{-1}{D}^{T}C)}^{T}X+X(A-B{F}^{-1}{D}^{T}C)\hfill \\ -\mathrm{XB}{F}^{-1}{B}^T+C^T{E}^{-1}C=0\hfill \end{array}$$

(35)

To simplify plan G, assumed that $D=0$ with a strictly proper plant. A minimal state-space representation of the normalized left coprime factorization of the shaped plant $G_s={\tilde{M}}^{-1}\tilde{N}$ is derived as follows

$$[\tilde{N},\tilde{M}]=\left[\begin{array}{ccc}A+\mathrm{LC}\hfill & B+\mathrm{LD}\hfill & L\hfill \\ \multicolumn{1}{c}{E^{-1/2}C}& E^{-1/2}D\hfill & E^{-1/2}\hfill \end{array}\right]$$

(36)

where $L=-(Z{C}^T)E^{-1}$ , $E=I$ and $I_n$ is the $n\times n$ identity matrix.

The performance index is given by minimizing the infinity norm of the transfer function matrix between exogenous signal vectors and objective signal vectors

$$min{\Vert \left[\begin{array}{c}K\hfill \\ \multicolumn{1}{c}{I}\end{array}\right]{(I+\mathrm{GK})}^{-1}{\tilde{M}}^{-1}\Vert }_{\infty }<\gamma$$

(37)

For a given $\gamma >{\gamma }_{min}$ , $H_{\infty }$ loop-shaping controller is obtained as

$$K_{\infty }=\left[\begin{array}{cc}A+\mathrm{BP}+{\gamma }^2{\left(Q^T\right)}^{-1}Z{C}^T(C+\mathrm{DP})\hfill & {\gamma }^2{\left(Q^T\right)}^{-1}Z{C}^T\hfill \\ \multicolumn{1}{c}{B^{T}X}& -D^T\hfill \end{array}\right]$$

(38)

Consequently, the controller law for a coaxial rotor helicopter is $K=W_1{K}_{\infty }{W}_2$ .

Control structure design

Some researchers have been investigating control structure of a helicopter. Heave–yaw and roll–pitch controllers were designed for a coaxial micro helicopter in the literature. ²⁷ Cai et al. ²⁸ proposed a hierarchically divided three-layer flight control structure law for a small-scale unmanned helicopter. Inner- and outer-loop architectures were introduced into the design and flight testing of the controller of a robotic helicopter. ²⁹ However, the coaxial rotor helicopter has additional input DOFs to control structure design. It is normally required to adopt a different control strategy and scheme. Therefore, it is implemented to select the reasonable input–output pairings with additional inputs. Relative gain array (RGA) is adopted to provide a suggestive rule to arrange dominant elements as diagonal as possible following the analysis of the mentioned singularity value responses. Its formulation is as follows

$$\mathrm{RGA}(G)=G\times {\left(G^{-1}\right)}^T$$

(39)

where the transfer function form of the plant is represented as follows

$$G(s)=C{(\mathrm{sI}-\mathrm{A})}^{-1}\mathrm{B}+\mathrm{D}$$

(40)

where $\mathrm{A}$ and $\mathrm{B}$ are state and control matrices, $\mathrm{D}=0_{6\times 6}$ and output matrix

$$\mathrm{C}=\left[\begin{array}{cccccccc}1\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 1\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0.3\hfill & 1\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 1\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0.15\hfill & 1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & 1\hfill \end{array}\right]$$

(41)

A mixed-rate technique is adopted. ${\theta }_q$ and ${\varphi }_p$ are the combination of the pitch-pitch rate and the roll-roll rate, respectively, in order to construct square plant.^30,31

Taking equations (39) and (40), the steady-state RGA matrix of the constructed model is rearranged using the rule of RGA

$$\mathrm{RGA}=\left[\begin{array}{cccccc}1.7958\hfill & 0.0001\hfill & -0.8384\hfill & 0.0752\hfill & -0.0326\hfill & -0.0000\hfill \\ \multicolumn{1}{c}{0.0004}& 3.1960\hfill & -0.0013\hfill & -0.0005\hfill & 0.0014\hfill & -2.1961\hfill \\ \multicolumn{1}{c}{-0.8348}& -0.0013\hfill & 1.7980\hfill & -0.0366\hfill & 0.0741\hfill & 0.0005\hfill \\ \multicolumn{1}{c}{0.0755}& -0.0000\hfill & -0.0329\hfill & 1.7960\hfill & -0.8386\hfill & 0.0000\hfill \\ \multicolumn{1}{c}{-0.0368}& 0.0029\hfill & 0.0742\hfill & -0.8340\hfill & 1.7954\hfill & -0.0017\hfill \\ \multicolumn{1}{c}{-0.0002}& -2.1977\hfill & 0.0004\hfill & -0.0002\hfill & 0.0002\hfill & 3.1974\hfill \end{array}\right]$$

(42)

Based on the singularity value response simulation results in section ‘Dynamic stability and singularity’ and rearranged RGA matrix in equation (42), the control architecture of the coaxial rotor helicopter is constructed, as shown in Figure 6.

Figure 6.

Block diagram of the controller.

The rearranged controlled variables order is $u,w,{\theta }_q,v,{\varphi }_p,r$ , and the input vector is $\mathbf{u}=[\begin{array}{c}{(A_1,B_1,{\theta }_{\mathrm{col}})}_{\mathrm{up}},{(A_1,B_1,{\theta }_{\mathrm{col}})}_{\mathrm{low}}\hfill \end{array}]$ . This new control scheme allows the simultaneous control of six states with six inputs. It is obvious that the structure distinguishes from the previous four-input control strategy. ³²

Control synthesis

Based on the theory in section ‘Control theory’, the implementation process of the control law is presented in Table 7.

Table 7.

Implementation process of theory.

Input: a minimal state-space realization G

Output: the $H_{\infty }$ loop-shaping controller

Process:

1. Scaling inputs and outputs, and shaping the plant.

2. Update the $W_1$ and $W_2$ to satisfy desired loop shape.

3. Robust stabilization. Solve the stability margin

4. Check if $\epsilon >0.3$ . If it does, calculate the controller K, otherwise re-back to the step 1.

5. Form the closed loop, and check the time responses for rise-time, overshoot, settle-time, and interaction level by performing the unit step responses. If the performance is not satisfied, re-run the step 1).

The ranks of $(A,B)$ and $(A,C)$ are as follows

$$\mathrm{rank}(A,B)=6$$

(43)

$$\mathrm{rank}(A,C)=6$$

(44)

It shows that the state-space realization of the model is truly minimal realization. Therefore, the control synthesis can be implemented in the following procedure.

Input and output scaling

Each trim value was approximately chosen as the maximum allowable input to scale the responding input. ⁸ The output scaling process was based on the effectiveness of the open-loop responses and performance requirements. In this study, the selected output maximum scale value was high when the output singular response value was high. ³³

Selecting pre- and post-compensators

The diagonal elements of pre-compensators are selected to boost a low-frequency loop gain and to reduce gain at high frequency. ²⁹ The integral terms are chosen to improve the performance of the closed loop, such as the roll-off rate at cross-over frequency and the low-frequency gain, while the proportional parts adjust to the lead or lag phase. Altering pre- and post-compensators can seriously affect trade-off between robust stability and closed-loop performance. After numerous iterative trials, the pre-compensator is finally selected as follows

$$\begin{array}{c}{W}_1=\mathrm{diag}\{\frac{1.15s+4.25}{s+0.01},\frac{1.15s+5}{s},\frac{1.65s+6.17}{s},\hfill \\ \frac{2.25s+4.75}{s+0.01},\frac{2.54s+5.55}{s},\frac{1.05s+4.75}{s}\}\hfill \end{array}$$

(45)

The diagonal elements of post-compensator are chosen as a constant in the following way

$$W_2=\mathrm{diag}\{0.1,1,1,0.1,1,1\}$$

(46)

Controller synthesis and check time responses

$H_{\infty }$ loop-shaping controller is synthesized with stable margin of 0.3691 (please refer to the ‘Support materials’ section). It ensures up to 36.91% uncertainty of the normalized coprime factors. Comparing the results in the literature,^29,33 this stable margin shows a competitive performance. Therefore, the feedback controller of the plant G is constructed. The eigenvalues of the closed loop are presented in Table 8.

Table 8.

Eigenvalues of the closed loop.

Eigenvalues
−20.4999 ± 0.8045i	−4.7777
−18.1799 ± 3.4587i	−4.0394
−10.7468 ± 0.7599i	−3.4645
−4.6146 ± 4.4976i	−3.5821
−2.6380 ± 4.2527i	−2.5885
−3.0333 ± 3.9859i	−2.3625
−4.0722 ± 3.9449i	−2.1015
−3.7291 ± 0.3931i
−2.5816 ± 1.0776i
−2.0040 ± 0.4213i

In Table 8, the real parts of eigenvalues are all negative. It shows that unstable model becomes stable after applying $H_{\infty }$ loop-shaping control. Time response analysis will be analysed in detail in next section.

Numerical simulations

Specifications

Guidelines for performance specifications were obtained from successful assessments of quality management of unmanned helicopters.^9,29 Economic specifications were explored and developed to satisfy the requirements of controller design as follows:

1. Stability margin is larger than 0.3.

2. The overshoot of a step-time response is less than 20% and does not exhibit residual oscillation.

3. The decoupling and signal tracking are acceptable to a certain extent.

Therefore, these specifications are set as constraints to achieve sufficient control for our preliminary stage.

Performance analysis and discussion

According to the specifications, numerous iterative trials have been conducted to alter pre- and post-compensators until sufficient performances were achieved. Since couplings among six channels are complex, it must be considered carefully to achieve trade-off among the specifications for robust stability margin and performance.

A class of good nominal performances is illustrated in Figure 7, including small-overshoot, appropriate rise-time, satisfied settle-time and reasonable decoupling. Rise-times of the step responses of two variables are approximately 0.7 and 0.5 s, respectively. These achieved rise-times also agreed well with the bandwidth of the mode of mixed-rate pitch and roll. A slight residual couplings remain in velocities but being suppressed with the prolonging of time, as shown in Figure 7.

Figure 7.

Step response on attitude: (a) pitch and pitch rate and (b) roll-roll rate.

In Figure 8, the step responses of the longitudinal and lateral velocities exhibit less overshoot, satisfied settle-time and reasonable decoupling. Rise-times of the step responses of longitudinal and lateral velocities are approximately 0.25 and 0.2 s, respectively. Rise-time tends to be faster than that the dynamic model bandwidth can handle. This phenomenon is mainly caused by high stability margin. It exhibits slightly higher overshoot but no oscillation. Satisfied decoupling is achieved.

Figure 8.

Step response of velocity: (a) lateral velocity and (b) longitudinal velocity.

In Figure 9, step responses of the heave and yaw channels have a little overshoot but less than 20%. Rise-times of step responses of the heave and yaw variables are approximately 0.7 and 0.5 s, respectively. The perfect decoupling performance in the two channels is also similar to that of the attitude channel. In Figure 9, the couplings of step responses among the controlled variables are suppressed with the prolonging of time. The satisfied settle-times are also achieved, as shown in Figure 9. Based on experiences, altering pre- and post-compensators is a complicated task, particularly when settling six channels simultaneously to achieve a perfect trade-off.

Figure 9.

Step responses of vertical and yaw velocities: (a) vertical velocity and (b) yaw rate.

In summary, the obtained stability margin is a perfect result to satisfy the performance specification (1). The overshoot of step response of each channel tends to calm. And they exhibit no residual oscillation in simulation results. Step responses of state variables also exhibit a perfect decoupling. And rise-times of step response are also reasonable. Moreover, the performance of achieved closed loop is relatively good although not perfect. Based on our experience, altering the pre- and post-compensators of a controller is a complicated task, particularly when settling six channels simultaneously to achieve a perfect trade-off. Adopting a hierarchical control structure provides a good idea on which possible solutions should be used to improve performance.

Conclusion

This article proposed a complete methodology for an unmanned coaxial rotor helicopter with unstructured uncertainties to achieve high-accuracy tracking performance from modelling to robust control. Conclusions are as follows:

1. A reasonable coaxial rotor helicopter prototype is constructed. This prototype has an ability to depict main characteristics of this type helicopter. Parameters agree well with experiments or empirical results from the open literature. In the parameter design, an empirical and verified process is implemented to decrease the budget and risk of programme. This design process can also contributed to other new rotorcraft design. Trim condition is calculated and analysed. Stability and control matrices are calculated through the quasi-steady assumption.

2. The $H_{\infty }$ loop-shaping control is proposed to apply to a coaxial rotor helicopter with unstructured uncertainties. Based on the investment of singularity value responses and RGA, a new control scheme is proposed to simultaneously incorporate six inputs to control six states. Coprime factor uncertainty model is considered. The implementation process of the control law is presented in detail.

3. Simulation results show the effectiveness of the proposed controller for a coaxial rotor helicopter with unstructured uncertainties to achieve high-accuracy tracking performance. The resulting controller provided a stability margin of 0.3619. It means that tolerated sizes of a class of uncertainties are up to 36.19. Step responses exhibit no residual oscillation, an overshoot response of less than 20%, reasonable decoupling and good rise-time response. This proposed methodology can also be applied to address the high tracking performance of controller in the presence of unstructured uncertainties in other new type rotorcraft.