## Abstract

A linear quadratic Gaussian (LQG) controller for active vibratory loads reduction in helicopters is proposed based on a revisited higher harmonic control (HHC) input by active trailing-edge flaps (ATEFs). Conventional individual blade control (IBC) input is redefined using N − 1/rev interblade phase lead, N/rev collective, and N + 1/rev interblade phase lag signals where 1/rev frequency modulation originates from the multiblade coordinate (MBC) transform. A Mach-scaled flap blade is designed and analyzed by the multibody dynamics analysis DYMORE. A linear time-invariant representation is identified from N/rev envelopes of the input and output responses obtained by DYMORE analysis. A matlab/simulink closed-loop control simulation is designed using the identified state-space realization. The N/rev vibratory loads are reduced up to 52% with flap deflections and the linear control results match well with the nonlinear responses obtained from DYMORE. Furthermore, the multivariable closed-loop stability estimated by the loop transfer functions using disk margin analysis reveals sufficient gain and phase margins.

## Introduction

Higher harmonic control (HHC) for hub vibratory load reduction in helicopters has been studied and tested by many researchers during the past decades by means of hydraulic swashplate actuation, individual blade pitch actuation, active twisting rotor blades, or active trailing-edge flaps (ATEFs) [1–3]. The most popular HHC approach to date is the T-matrix assumption for the relation between control input and the hub load perturbations [4–7]. Another method is an HHC-based continuous-time controller using a notch filter targeting $N$/rev frequencies [8–10]. Although the previous T-matrix schemes showed sufficient reduction of vibratory loads by comprehensive rotor analysis and experiments, systematic design of controller simulations and comparisons with the closed-loop test results have been rarely presented. Further, as it is described by Hall, Bittanti et al., and Uddin et al., the conventional T-matrix approach lacks flexibility regarding the stability and robustness, since the relevant dynamics is formulated by a simple algebraic relationship [11–14].

ADS-33 requires linear stability analyses for aircraft handling qualities [15]. Furthermore, the linear system controller design is a convenient and well-established framework for the most control applications, and the evaluation of relative control margins provides control system robustness to model uncertainties. In an effort to incorporate linear control law to replace the T-matrix approach, Bittanti used collective control input and hub vertical force output to represent the rotor in a state-space form [12,13]. Uddin derived a state-space model of each blade deformation using only the first modes of flapping, lead lag, and torsion with simple two states for each mode [14]. However, the linear models used in those studies did not consider dynamics for $N$/rev frequencies and were not capable of reducing all the significant hub vibratory load components. Thus, it enforced direct comparisons with HHC controllers to be difficult.

As the helicopter rotor is rotating with $N$-blades, it may be characterized by a linear time-periodic (LTP) mathematical model with $N\Omega $-frequency at equilibrium. Then, a linear time-invariant (LTI) rotor representation is further extracted from the LTP rotor description. Siddiqi et al. used a frequency-domain estimation technique to find the relevant harmonic transfer functions [16–18]. Later, Cheng et al. and Padthe et al. explored a Fourier-expansion based modeling and extracted elaborate LTI representations [19,20]. However, those models provided only analytical formulations, which required specific aeroelastic equations of motion. To construct a realistic control simulation, LTI models need to be estimated from the input signals and output measurements in either the time or the frequency domain.

Eventually, it is essential to effectively define the input and output so that they satisfy a linear relationship. For example, the input and output frequencies should be matched. The conventional individual blade or on-blade control approaches, such as active twisting rotor or ATEF, are based on the individual blade pitch command, while the swashplate actuation is based on the fixed-frame manipulation. In general, control signal for each blade is defined in a rotating blade frame, but hub vibratory load is evaluated in the fixed hub frame due to the practical requirements of configuring the load sensor unit in the rotating frame. The individual blade control (IBC) assigns an interblade pitch phase of $+2\pi (i\u22121)/N,\u2009(i=1,\u20092,\u2026,\u2009N)$ to each blade relative to the reference blade [21,22]. The reference blade input consists of $(N\u22121)$, $N$, and $(N+1)$/rev harmonics, with cosine and sine basis functions. However, when the IBC inputs are transformed to the fixed frame by multiblade coordinate (MBC) transform, a shift of $\u22121$/rev frequency will occur. On the other hand, if the individual pitch command is assigned using a negative interblade phase of $\u22122\pi (i\u22121)/N,\u2009(i=1,\u20092,\u2026,\u2009N)$, a $+1$/rev frequency modulation will be obtained in the fixed frame input. There will be no frequency shift for the collective blade inputs. Consequently, an *N*/rev conventional IBC input generates an $(N\u22121)$/rev fixed frame control input, which results in an $(N\u22121)$/rev fixed frame hub load perturbation that is not a target vibratory load frequency. By the definition of LTI system, the frequency response is defined in terms of each corresponding output over input signal frequency. Since the frequency of the individual flap deflection command is different from that of the fixed frame hub load due to MBC transformation, a proper individual command to match the frequencies between the fixed frame command and output hub load is needed to define a LTI system. Based on this fact, this paper suggests a novel revised HHC-IBC input considering the MBC transform. To generate a fixed frame $N$/rev frequency input, an $N$/rev collective, an $(N+1)$/rev positive interblade phase, and an $(N\u22121)$/rev negative interblade phase input will be configured.

The objective of this paper is to design a linear quadratic Gaussian (LQG) regulator based on the identified LTI representation to evaluate the closed-loop stability and sensitivity. Since the vibratory load transfer structure is a multi-input and multi-output (MIMO) system, the LQG control is preferred for initial control and effective reduction of multimode vibrations. The remaining sections of this paper introduce: (1) baseline multibody aeroelastic rotor analysis containing trailing-edge flap components and relevant analyses; (2) LTI identification of the flap rotor; (3) LQG-based design of the closed-loop vibratory load controller; and (4) vibration reduction control simulation and robust stability analysis.

## Mach-Scale Seoul National University Flap Rotor Analysis

A brief introduction of a comprehensive rotor analysis for a Mach-scale flap rotor, named Seoul National University Flap (SNUF) rotor, will be given in this section. Figure 1 shows the fabricated prototype blade. Table 1 summarizes the specifications of the SNUF rotor blade.

Description | Property | Value |
---|---|---|

Solidity ratio | $\sigma $ | 0.1146 |

Number of blades | $N$ | 4 |

Airfoil | NACA | 0015 |

Tip Mach number | $Ma$ | 0.6 |

Radius (m) | $R$ | 1.5 |

Chord (m) | $c$ | 0.135 |

Lock number | $\gamma $ | 4.8 |

Thrust (N) | $T$ | 2500 |

Rotating speed (rpm) | $\Omega $ | 1300 |

Flap center location | — | 75%$R$ |

Flap length | — | 20%$R$ |

Flap chord length | — | 15%$R$ |

Description | Property | Value |
---|---|---|

Solidity ratio | $\sigma $ | 0.1146 |

Number of blades | $N$ | 4 |

Airfoil | NACA | 0015 |

Tip Mach number | $Ma$ | 0.6 |

Radius (m) | $R$ | 1.5 |

Chord (m) | $c$ | 0.135 |

Lock number | $\gamma $ | 4.8 |

Thrust (N) | $T$ | 2500 |

Rotating speed (rpm) | $\Omega $ | 1300 |

Flap center location | — | 75%$R$ |

Flap length | — | 20%$R$ |

Flap chord length | — | 15%$R$ |

### Multibody Aeroelastic Modeling of Seoul National University Flap Rotor.

The one-dimensional geometrically exact beam elements and two-dimensional cross-sectional properties of the hingeless SNUF blade are analyzed using DYMORE and VABS [23,24]. Figure 2 shows the components of the DYMORE modeling employed in this study. The hingeless configuration is represented by a pitch bearing and equivalent bearing stiffness calculated by the pitch link stiffness of the rotor test stand. The blade is divided into three beam sections, and the middle section is attached with the bracket and flap assembly. The flap mechanism is simplified by revolute, universal, prismatic, and spherical joints. Twenty-eight and five quadratic elements are used for the main blade and the flap sections, respectively. Aerodynamic loads are estimated by the internal aerodynamic interfaces of DYMORE that include the three-dimensional finite state dynamic inflow model by Peters and He [25,26]. In this paper, 15 modes of the dynamic inflow are selected, which are equivalent to 136 states. The number of aerodynamic interface nodes is set to 30 with Gaussian distribution.

Figure 3 shows the cross-sectional design of the representative section. An asymptotically exact cross-sectional analysis is performed using VABS and PreVABS. The cross-sectional design optimization framework used in this study was initially developed in previous studies [27–29]. The spar structure is modified to readily incorporate the flap drive shaft bearing in the prototype blade. Figure 4 shows the planform of the blade and X12 to X21 indicate cross-sectional analysis stations. The sectional properties are linearly interpolated between the discrete spanwise stations. Figure 5 depicts the sectional properties distribution.

In the author's previous study, bench experiments on the flap mechanism, tensile testing of the SNUF blade, and modal tests results are presented [30]. Figure 6 shows that natural frequencies and the first mode shapes of the prototype blade match well with the design values. Further, Fig. 7 depicts rotating natural frequencies of SNUF blade, HART-II blade, and BK-117 flapped blade [10,31]. The first torsional frequency is 3.8/rev, which is designed to be lowered to obtain sufficient control authority as other active blades are [10,31,32].

### Baseline Trim Analysis.

To select a typical analysis condition representing significant vibratory loads, 4/rev vibratory loads are identified for a steady forward flight state at advance ratios of 0.1 to 0.38, where $CT/\sigma $ is 0.06 and the shaft tilt angle $\alpha s$ is 6 deg. A wind-tunnel trim with target thrust and zero hub moment is implemented using AUTOPILOT function in DYMORE. The initial collective pitch angle is set to obtain an average of zero vertical hub loads in the perturbation identification period, resulting in a fast convergence and reliable trim results.

The finite state dynamic wake model in the present DYMORE analysis is compared with the prescribed and free wake level of CAMRAD-II analysis to verify the trim analysis aerodynamic environment [33]. As described by Peters, the finite state dynamic wake model is essentially a prescribed-wake analysis and, as shown in Figs. 8 and 9, it correlates well with the CAMRAD-II prescribed-wake analysis. Although it cannot account for the wake roll-up or the detailed blade-vortex interactions, it is useful for higher harmonic controller design.

As shown in Fig. 10, the maximum VI is obtained at advance ratios between 0.15 and 0.2. There, an advance ratio of 0.15 is selected as a representative trim condition in this study.

## Development of a Linear Quadratic Gaussian Control Law

This section provides a systematic procedure to identify the present LTI flap rotor and design technique of an LQG-based closed-loop HHC control law.

### Multiblade Coordinate Transform and Individual Flap Control Input Allocation.

*N*-bladed rotor, the transformation matrix $M$ for the

*n*th harmonic is defined as Eqs. (4)–(6). Subscripts $R$ and $NR$ indicate rotating and nonrotating frames, respectively. The $n$th harmonic is specified by the number of blades, where $n=1,2,\u2026,(N\u22122)/2$ for even $N$ and $n=1,2,\u2026,(N\u22121)/2$ for odd $N$. By applying an MBC transformation to individual blade frame flap deflections, one can interpret the relation between fixed frame flap perturbation and fixed frame hub loads

*N*/rev hub load perturbations

Equation (7) describes the individual blade frame input signal for a $N$-bladed rotor, where the superscripts in parentheses indicate the blade number, and superscripts without parentheses are the amplitudes of the $(N\u22121)$, $N$, and $(N+1)$/rev components. The input signal is divided into cosine and sine basis functions to represent all the phases over the period. The fixed frame perturbed output signal is also decomposed into cosine and sine basis functions, as described in Eq. (8), by finding the phase lag from the least-square analysis.

### Linear Time-Invariant Identification Algorithm of an on-Blade Control System.

Off-line identification is performed by first setting the state-space representation of individual flap inputs and hub $N$/rev outputs, as shown in Eqs. (7) and (8). Since each signal component of the input vector $[\theta cN\theta sN\theta cN+1\theta sN+1\theta cN\u22121\theta sN\u22121]T$ is orthogonal in its cosine and sine directions, input signal for system identification is assigned independently. For example, $[100000]T$ gives the output perturbed step response vector $[\delta FxcN\u2009\delta FxsN\u2009\delta FycN\u2009\delta FysN\u2009\delta FzcN\u2009\delta FzsN\u2009\delta MxcN\u2009\delta MxsN\u2009\delta MycN\u2009\delta MysN\u2009\delta MzcN\u2009\delta MzsN]T$ to the collective $N$/rev cosine input. The number of output $No$ is the number of output vector components to control, and the number of input $Nu$ is 6. As a result, a $No\xd7Nu$ transfer matrix will be constructed, from which a state-space model will be obtained.

If the target rotor has $N$ blades, $N$/rev perturbed output load will be found by splining the maximum or minimum points of the oscillatory component corresponding to $N$/rev frequency of the output signal. To search accurately for the maximum or minimum points in $N$/rev frequency, the minimum distance constraint of 1/$N$ rotor period length can be used. In this paper, matlab command “envelope” in the signal processing toolbox^{TM} is used. The linear step response shape is then generated by applying cubic splines to the $N$/rev peak locations. Finally, an LTI transfer function is found using the instrumental variable method, which is implemented by matlab “tfest” command [38].

To test the LTI representation, only the hub vertical shear, pitching, and rolling moments are selected as outputs, which gives output vector $[\delta FzcN\u2009\delta FzsN\u2009\delta MxcN\u2009\delta MxsN\u2009\delta MycN\u2009\delta MysN\u2009]T$. Flap actuation is applied once the dynamic rotor response converges to the steady-state. Figures 13–15 show the perturbed hub vertical load responses to several inputs and the estimated responses from the LTI transfer function multiplied by the 4/rev sinusoid found by least-square fitting. The initial transient behaviors are accurately estimated by the LTI transfer functions. However, due to the LTP characteristics of the 4/rev response, the 0 and 8/rev frequency components may remain at equilibrium. Thus, there are small deviations between the DYMORE prediction and the LTI response.

Figure 16 shows the frequency response of the state-space representation and the reduced-order system obtained with the balanced singular perturbation approximation algorithm, which is implemented by matlab “balred” command [39]. The initial system estimate has 216 states and the reduced-order model has 44 states. The resulting state-space representation of the flap rotor features has fully populated the A, B, C, and D matrices.

### Linear Quadratic Gaussian Synthesis for Vibratory Load Reduction.

In general, linear quadratic regulator (LQR) controllers exhibit benign stability margins. The gain margin will be infinite for a stable plant, or at least $\u22126\u2009$dB of the gain reduction margin for an unstable plant. Furthermore, the minimum phase margin will become $\xb160\u2009deg$ if the system modeling is correct. However, if the set of the measurable state variables is not selected appropriately, the guaranteed margins will not be satisfied [40]. In the present LTI represented IBC system, the dynamics of the internal states are unknown because of the complicated transfer function set and the order reduction scheme used. To avoid such problem, a dynamic state feedback using the Kalman Filter as the state observer can be used.

Controllability and observability are evaluated by constructing the controllability and observability matrices using the identified flap rotor state space representation. Because the system matrices are fully populated, the corresponding controllability and observability matrices will be in a full rank. Furthermore, the Ricatti differential equation is solved with suboptimal fashion, as the flap rotor is reachable and controllable. The control input is estimated by the algebraic Ricatti equations in Eqs. (11) and (12).

## Closed-Loop Vibration Control Simulation

In this section, the results of the LQG-based HHC control response from LTI representation in matlab/simulink and the nonlinear DYMORE response to the identical flap deflection input are presented, as well as the closed-loop stability analysis.

### Controller Evaluation in Steady Flight Conditions.

The wind-tunnel trim isolated rotor configuration in the present DYMORE analysis has invariant primary pitch control settings. By following the system identification and controller design procedure mentioned in the Development of a Linear Quadratic Gaussian Control Law section, a linear control law is established for the trim condition.

At first, it is important to check whether any steady hub load values are influenced by the flap controls so that the performance of the rotor do not change. Figures 18(a)–18(c) show steady-state one-revolution closed-loop responses from the matlab/simulink simulation, which uses identified LTI representation of the rotor and DYMORE nonlinear analysis results using the same flap deflection in Fig. 18(d). It is verified that the present control scheme does not degrade the original rotor performance, for the steady values of the thrust and hub moments remain the same level in the baseline rotor.

As it is explained in the Development of a Linear Quadratic Gaussian Control Law section, there may exist discrepancy from DYMORE results due to the periodic components, for the present LTI model only reproduces $N$/rev components. Nevertheless, the LTI representation of the flap rotor is valid as both outputs coincide in terms of magnitude and phase response of dominant $N$/rev signals. Figure 19 shows the 4/rev vibratory loads component of responses in Fig. 18. All the components in the hub vibratory loads are reduced by an average of 52%, even when the hub in-plane shear loads are not controlled.

### Closed-Loop Stability Analysis.

As the original nonlinear rotor dynamics are linearized at an equilibrium point to obtain the LTI representation, there will always be some uncertainty, which will become worse in practical situations under sensor calibration errors or uncertain noise. Therefore, the control law used should satisfy the minimum gain and phase margin requirements. However, according to Doyle, there were no guaranteed stability margins for a combined Kalman filter and LQR control law [41]. In this paper, to investigate the present MIMO control system margins, the disk margin analysis based on the structured singular value decomposition is directly applied to LQG simulink block using matlab Robust Control Toolbox [42]. Loop transfer function is constructed by setting a break point at the plant input. There are two functions in the disk margin analysis. The loop-at-a-time disk margin analysis perturbs gain and phase of a single channel while others are hold. On the other hand, multiple disk margin analysis introduces simultaneous, independent variations in each channel that generally estimates smaller gain and phase margin of the closed-loop than the classic disk margin analysis does. Detailed formulation of the multiloop margin analysis is included in Ref. [43]. Figure 20 shows the corresponding matlab/simulink block diagram for the loop transfer analysis. Figure 21 depicts the allowable gain and phase variations of the present control system. Disk margin analysis for input channel 1 ($\theta cN$) results in gain margin of 12.7 dB and phase margin of $\xb1$63.8$\u2009deg$, whereas the multiple disk margin analysis gives the minimum gain margin of 10.7 dB and the minimum phase margin of $\xb1$57.5$\u2009deg$, which is analogous to the minimum phase margin of $\xb1$60$\u2009deg$ in a separate LQR-controlled system.

## Conclusions

An improved LQG-based higher harmonic control law is investigated by applying the principle of MBC transformations. Individual blade input signals are redefined according to the phase shift characteristics of the frame transformation. The time-invariant system representations are obtained from every step response of the $N$/rev input and output measurements. A reduced order state-space realization is conducted and a LQG controller is designed for the present plant. The major conclusions drawn from this study are as follows:

Since the LTI representation is based on the one-to-one correspondence of the input–output frequency, it will be required to match the frequencies of the fixed-frame input and output. Accordingly, the presented individual higher-harmonic control input will be simple to implement in real situations.

A detailed LTI state-space representation procedure of the flap rotor from the improved individual HHC flap input is presented and validated by comparison with the nonlinear aeroelastic analysis. Then, using the well-established linear control design rules, relative margins required by the regulations are obtained in a straightforward fashion. Furthermore, the LQG controller design is preferred to the original LQR design because a set of internal states cannot be selected by measurement. The LQG controller is expected to provide more margins thanks to its dynamic state feedback capability.

The closed-loop responses obtained from the LTI simulation and the DYMORE nonlinear analysis match well. This indicates sufficient reliability to implement the controller within the control margins.

Robust control stability analysis theory is readily applicable to the LQG-based HHC controller. The closed-loop multiple disk margin stability analysis reveals the minimum values of 10.7 dB for the gain margin and $\xb1$57.5 deg for the phase margin, and appropriate selection on the control gain based on the vibration reduction performance or the control stability will be studied in future works.

In the future, a hover test stand will be manufactured. Prior to wind tunnel tests, endurance and rotational testing of the SNUF blades at several collective pitch angles will be conducted. Then, effectiveness of the present HHC controller will be evaluated accurately by conducting the dynamic characteristics measurement of SNUF. Further, in the upcoming wind-tunnel tests, LTI characteristics of SNUF rotor will be identified and the present HHC control algorithm will be compared against the experimental result.

## Acknowledgment

This work was conducted at the High-Speed Compound Unmanned Rotorcraft (HCUR) research laboratory with the support of the Agency for Defense Development (ADD). It was supported by the Advanced Research Center Program (NRF-2013R1A5A1073861) through a grant from the National Research Foundation of Korea (NRF), funded by the Korean government (MSIP), and contracted through the Advanced Space Propulsion Research Center at Seoul National University.

## Funding Data

National Research Foundation of Korea (Funder ID: 10.13039/501100003725).

## Nomenclature

- $CT$ =
thrust coefficient

- $Fx,NP$ =
*N*/rev hub forward vibratory shear load - $Fy,NP$ =
*N*/rev hub sideward vibratory shear load - $Fz,NP$ =
*N*/rev hub vertical vibratory shear load - $MNR\u2192R$ =
nonrotating to rotating frame multiblade coordinate transform matrix

- $MR\u2192NR$ =
rotating to nonrotating frame multiblade coordinate transform matrix

- $Mx,NP$ =
*N*/rev hub vibratory rolling moment - $My,NP$ =
*N*/rev hub vibratory pitching moment - $N$ =
number of blades

- $Ni$ =
number of inputs

- $No$ =
number of outputs

- $Q$ =
vibratory loads weighting matrix

- $R$ =
flap control input weighting matrix

- $R$ =
blade radius

- $VI$ =
vibration index

- $\alpha s$ =
hub shaft tilt angle

- $\alpha i$ =
vibratory loads weighting coefficient,

*i*th load - $\beta i$ =
flap control input weighting coefficient,

*i*th input - $\sigma $ =
solidity ratio

- $\psi $ =
blade azimuth angle

- $\Omega $ =
rotational speed