Fuzzy Neural Network PID Control Used in Individual Blade Control

Abstract

In order to further reduce the vibration level of helicopters, the active vibration control technology of helicopters has been extensively studied. Among them, individual blade control (IBC) independently applies high-order harmonics to each blade with an actuator, which can improve the aerodynamic environment of the blade and effectively reduce the vibration load of the hub. The rotor structural dynamics model based on the Hamilton energy variation principle and the medium deformation beam theory were established firstly, and the aerodynamic model based on the dynamic inflow model and the Leishman–Beddoes unsteady aerodynamic model were also established. The structural finite element method and the direct numerical integration method were used to calculate the vibration response of the rotor to determine the vibration load of the hub. After these, the steepest descent-golden section combinatorial optimization algorithm was used to find the optimization parameters of IBC. Based on this, the input parameters of fuzzy neural network PID control were determined, and the rotor hub vibration load control simulation was conducted. Under the effect of IBC, the vibration loads of the hub could be reduced by about 60%. The article gives the best control laws of individual harmonic pitch control and their combinations. These results can theoretically be applied to the design of control law to reduce helicopter vibration loads.

1. Introduction

As a rotorcraft, the helicopter takes the rotor as the main power source for flight. It has many advantages that cannot be replaced by fixed wing aircraft, such as vertical take off and landing, hovering and flying in any direction. As a source of helicopter lift, the rotor also becomes a source of vibration. In the past, many passive methods, such as double-wire pendulums, were used to reduce helicopter rotor vibration loads [1]. Helicopter active control technology has been studied for helicopter vibration reduction since the 1970s, and this research has developed rapidly since the 1980s. Research into active control technology has also developed from the high-order harmonic control (HHC) to individual blade control (IBC), active flap control (AFC) [2], active torsion control (ATC) [3] and structural response active control (SRAC) [4]. AFC installs several flaps on the trailing edge of each blade. The flaps are driven by actuators embedded inside the blade, and generate flap angle changes according to the variation law provided by the control law. The ATC lays smart material on the blade, and the material will deform under the given electrical signal, so that during the rotation of the blade the angle of attack of different blade profiles will change according to the preset control law, so as to improve the aerodynamic environment of the blade and reduce the rotor vibration load. SRAC is achieved by using an actuator to excite the important force transmission components of a helicopter. By adjusting the amplitude, frequency, and phase of the excitation in real time through the controller, the vibration response generated by the active excitation force at key parts such as the cockpit is offset against the body vibration caused by the original excitation force of the helicopter.

IBC is a new method developed from HHC [5]. The basic principle of IBC is shown in Figure 1. HHC uses hydraulic actuators as active control actuators to apply excitation to the non-rotating ring of the swashplate. The frequency of the excitation force is determined by the number of blades. Compared to HHC technology, the actuator of the IBC is installed above the rotation ring of the swashplate, and each blade can independently achieve multi-harmonic pitch changes and arbitrary pitch motion.

Between 1977 and 1985, Ham conducted an early experimental study of IBC to improve blade bending stress, vibration loads, and stall-induced vibrations [6,7]. In the 1990s, Jacklin conducted full-scale tests of a BO-105 individual rotor blade pitch control system in the 40- by 80-Foot Wind Tunnel at the NASA Ames Research Center, with the main goal of reducing Blade Vortex Interaction (BVI) noise abatement [8]. The test results showed that the BVI noise could be reduced by more than 4 dB with 0.8 deg 2/rev IBC input at an advance ratio of 0.1. In 2015, Nitzsche F and D’Assuncao developed an independent blade control system centered on piezoelectric materials. In a systematic study, the control actuation of the actuator was achieved by overcoming the stroke limitation of piezoelectric materials. Experiments were conducted in the laboratory on rotors equipped with piezoelectric material independent blade actuators, verifying that piezoelectric material actuators can achieve vibration reduction in the rotor through independent blade control [9]. In 2017, Küfmann and Bartels designed an independent blade control system for a five bladed rotor composed of a multi-tilter rotating ring system [10]. In 2022, Yang used a genetic algorithm to optimize the vibration loads of a rotor with an individual blade control system [11]. The article considers the impact of individual harmonics of IBC on rotor vibration loads. Gao used a fuzzy neural network and particle swarm algorithm to search the optimal control parameters of individual blade control [12]. The 3/rev vertical vibration load of the helicopter in hover was reduced by 80% under the effect of IBC.

Neural network control is an intelligent control and identification method formed by simulating the physiological structure of the human brain. Fuzzy neural networks were proposed by Lee in 1975 [13] and were used for model parameter identification [14,15,16].

The Takagi Sugeno (TS) type fuzzy model was implemented by Gang [17] based on a genetic algorithm, and the control design was carried out by Mathiyalagan [18]. Lemos proposed a fast learning algorithm based on weighted non-norm neural networks in 2012 [19]. Afterwards, Ahn proposed a passive exponential filter with delay and external disturbances suitable for TS type neural networks [20]. In 2021, Fei designed a double hidden-layer recursive fuzzy neural network controller using terminal sliding mode control for nonlinear systems [21]. In 2021, Omar used the genetic algorithm optimization technique to optimize the distribution of the fuzzy membership functions in order to improve the performance of the suggested controller [22]. In 2023, Sabzalian proposed a new methodology based on an improved bidirectional recurrent neural network for lung cancer detection from images [23].

A mechanical actuator which can apply IBC to the rotor was designed, as shown in Figure 2. Previous studies only considered single harmonics in optimization analysis. This article considers the optimization results of the simultaneous action of second and third harmonics. The article adopts the steepest descent method based on the golden section to search for the optimal parameters of each harmonic when the helicopter is flying forward rather than just hovering. On this basis, a PID controller based on a fuzzy neural network was designed to control the rotor vibration load. Without selecting the theoretical optimal value, combined with the processing accuracy of the actuator, an approximate value of the amplitude, which can be achieved by machining, was used in the simulation.

2. Calculation Methods

2.1. Aeroelastic Model

2.1.1. Rotor Blade Airfoil Aerodynamics

The angle of attack of the blade is determined by the pitch angle$\theta$ and the flow angle $\varphi$ as shown in Figure 3 (where $\mathsf{\alpha }$ is angle of attack and $u_P$ and $u_t$ are vertical and lateral components of wind speed).

The pitch angle of a traditional helicopter is determined by the collective pitch angle,${ \theta }_0$, the cyclic pitch angles, ${\theta }_{1c}, {\theta }_{1s}$, and the blade linear pretwist, ${\theta }_w$. A helicopter, controled by IBC, adds highorder harmonic components of different orders, amplitudes and phases to the pitch angle of a traditional helicopter. The pitch angle $\theta$ can be written as in Equation (1).

$$\theta \left(t,\overline{r}\right)={\theta }_0+{\theta }_{1c}\cos (\Omega t)+{\theta }_{1s}\sin (\Omega t)+{\theta }_w(0.7-\overline{r})+A_2\cos (2\Omega t+{\phi }_2)+A_3\cos (3\Omega t+{\phi }_3)+\dots$$

(1)

where $\overline{r}$ is the relative radius of blade section position. $A_2, A_3$ and ${\phi }_2,{\phi }_3$ are the amplitudes and phases of higher order harmonics.

The profile airflow is more complex than traditional helicopters, requiring more accurate airfoil aerodynamic models. This paper calculates the unsteady aerodynamic forces of blade profiles based on the Leishman–Beddoes semi-empirical unsteady/dynamic stall aerodynamic model [24].

2.1.2. Elastic Model

The rotor blade can be analyzed as an elastic beam [25,26]. The equations of motion of the blade system were obtained through the method proposed by Hamilton, which can be expressed as:

$$\delta \Pi ={\int }_{t1}^{t2}\left(\delta U-\delta T-\delta W\right)dt=0$$

(2)

where $\delta U$ is the virtual variation of strain energy, $\delta T$ is kinetic energy and $\delta W$ is the virtual work performed by external forces. The variation of these can be expressed as:

$${\delta U}_i={\int }_0^R\iint \left(E{\epsilon }_{xx}\delta {\epsilon }_{xx} + G{\epsilon }_{x\eta }\delta {\epsilon }_{x\eta } + G{\epsilon }_{x\zeta }\delta {\epsilon }_{x\zeta }\right)d\eta d\zeta dx$$

(3)

$${\delta T}_i={\int }_0^R\iint {\rho }_s{\stackrel{⃑}{V}}_i·{\delta \stackrel{⃑}{V}}_{i}d\eta d\zeta dx$$

(4)

$$\delta W_i={\int }_0^R\left(L_u^A\delta u + L_v^A\delta v + L_w^A\delta w{ + L}_{\varphi }^A\delta \varphi \right)dx$$

(5)

where ${\epsilon }_{xx}$ is axial strain, and ${\epsilon }_{x\eta }$and ${\epsilon }_{x\zeta }$are engineering shear strains. $x\eta \zeta$ is the rotating deformed blade coordinate system. $L_u^A$, $L_v^A$, and $L_w^A$ are the distributed airloads in the x, y, z directions, respectively, and $L_{\varphi }^A$ is the aerodynamic pitching moment about the undeformed elastic axis.

For the i-th blade, $\delta \Pi$ is expressed as:

$$\delta {\Pi }_i={\int }_{\psi I}^{\psi F}[\sum _{j=1}^N(\delta U_j-\delta T_j-\delta W_j)]d\psi =0$$

(6)

Using the notation

$${\mathsf{\Delta }}_j=\delta U_j-\delta T_j-\delta W_j$$

(7)

Each section of the blade is simplified as a small beam. Each beam has twelve degrees of freedom. The elemental nodal displacement vector is defined as:

$$q_j^T=[u_1,{v_1,v}_1^{\prime }{,w}_1{,w}_1^{\prime },{\varphi }_1,u_2,v_2,v_2^{\prime }{,w}_2{,w}_2^{\prime },{\varphi }_2]$$

(8)

Using appropriate shape functions, ${\Delta }_j$ can be written as:

$${\Delta }_j=\delta q_j^T\left({\left[M\right]}_j{\ddot{q}}_j+{\left[C\right]}_j{\dot{q}}_j+{\left[K\right]}_j{q}_j-{\left\{F\right\}}_j\right)$$

(9)

where ${\left[M\right]}_j$, ${\left[C\right]}_j$ ${\left[K\right]}_j$and ${\left\{F\right\}}_j$ are mass matrices, damping matrices, stiffness matrices and load matrices of each element.

By assembling elemental matrices, the total energy can be expressed as:

$$\delta {\Pi }_i={\int }_{\psi I}^{\psi F}\delta q^T(\left[M\right]\ddot{q}+\left[C\right]\dot{q}+\left[K\right]q-\{F\})d\psi =0$$

(10)

The above formula can be transformed into the following finite element equation of the motion of the blade.

$$\left[M\right]\ddot{q}+\left[C\right]\dot{q}+\left[K\right]q=\left\{F\right\}$$

(11)

2.1.3. Response Solution

The response of the blade is calculated based on the previous aerodynamic and dynamic calculations. The modal superposition method, which obtains the response of the system by overlaying the contributions of various modes [27,28], was adopted to solve the problem. The method can be expressed as in Equation (12).

$$Y\left(t\right)=\sum _{i=1}^N{y}_i{\gamma }_i\left(t\right)=\varphi \gamma$$

(12)

where $Y\left(t\right)$ is the displacement of the blade; $y_i$ is the i-th natural vibration mode of the blade; ${\gamma }_i$ is the response coefficient corresponding to the i-th mode; and $\varphi$ is the modal matrix of the blade. $\varphi$ was derived from the $\left[M\right], \left[C\right], \left[K\right]$ in Equation (11) by using the EIG function in matlab.

The modal superposition method requires the truncation of the modal, which will reduce the computational accuracy. In order to balance the calculation accuracy and time cost, this paper selected the first 8 modes of the blade for calculation.

2.2. Hub Loads

In the analysis, the blade was divided into several segments and the blade load was calculated from the superimposition of each segment. With the assumption that the helicopter is flying steady, the hub loads in the flap plane can be seen in Figure 4. The loads mainly consist of gravity, lift and centrifugal forces. The loads in the lag plane can be seen in Figure 5. The loads mainly include aerodynamics, centrifugal forces, Coriolis forces and inertial forces.

The aerodynamic force of the airfoil profile can be obtained from Section 2.1.1, and the aerodynamics of the blade section are shown in Equation (13).

$$\begin{array}{c}dY=\frac{1}{2}{C}_L\rho V^{2}cdr\\ dX=\frac{1}{2}{C}_D\rho V^{2}cdr\\ dL=dYcos\varnothing +dXsin\varnothing \\ dQ=dXcos\varnothing -dYsin\varnothing \end{array}$$

(13)

The forces in each segment can be expressed as follows:

$$\begin{array}{l}{\mathrm{d}\mathrm{T}}_{\mathrm{s}}=\mathrm{d}\mathrm{L}\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{w}\mathrm{\prime }\\ {\mathrm{d}\mathrm{H}}_{\mathrm{H}}=\mathrm{d}\mathrm{Q}\mathrm{s}\mathrm{i}\mathrm{n}\mathsf{\psi }-\mathrm{d}\mathrm{L}\mathrm{c}\mathrm{o}\mathrm{s}\mathsf{\psi }\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{w}\mathrm{\prime }\\ {\mathrm{d}\mathrm{S}}_{\mathrm{s}}=\mathrm{d}\mathrm{Q}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{c}\mathrm{o}\mathrm{s}\mathsf{\psi }-\mathrm{d}\mathrm{L}\mathrm{s}\mathrm{i}\mathrm{n}\mathsf{\psi }\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{w}\mathrm{\prime }\\ {\mathrm{d}\mathrm{N}}_{\mathrm{s}}=\mathrm{m}{\mathsf{\Omega }}^2\mathrm{x}\mathrm{d}\mathrm{r}·\mathrm{s}\mathrm{i}\mathrm{n}\mathsf{\psi }+\mathrm{m}{\mathsf{\Omega }}^2\mathrm{y}\mathrm{d}\mathrm{r}·\mathrm{c}\mathrm{o}\mathrm{s}\mathsf{\psi }\\ {\mathrm{d}\mathrm{N}}_{\mathrm{H}}=\mathrm{m}{\mathsf{\Omega }}^2\mathrm{x}\mathrm{d}\mathrm{r}·\mathrm{c}\mathrm{o}\mathrm{s}\mathsf{\psi }-\mathrm{m}{\mathsf{\Omega }}^2\mathrm{y}\mathrm{d}\mathrm{r}·\mathrm{s}\mathrm{i}\mathrm{n}\mathsf{\psi }\\ {\mathrm{d}\mathrm{g}}_{\mathrm{s}}=2\mathrm{m}\mathsf{\Omega }\dot{\mathrm{x}}\mathrm{d}\mathrm{r}·\mathrm{c}\mathrm{o}\mathrm{s}\mathsf{\psi }\\ {\mathrm{d}\mathrm{g}}_{\mathrm{H}}=2\mathrm{m}\mathsf{\Omega }\dot{\mathrm{x}}\mathrm{d}\mathrm{r}·\mathrm{s}\mathrm{i}\mathrm{n}\mathsf{\psi }\\ {\mathrm{d}\mathrm{C}}_{\mathrm{t}}=\mathrm{m}\ddot{\mathrm{z}}\mathrm{d}\mathrm{r}\\ {\mathrm{d}\mathrm{C}}_{\mathrm{s}}=\mathrm{m}\ddot{\mathrm{y}}\mathrm{d}\mathrm{r}·\mathrm{c}\mathrm{o}\mathrm{s}\mathsf{\psi }\\ {\mathrm{d}\mathrm{C}}_{\mathrm{H}}=\mathrm{m}\ddot{\mathrm{y}}\mathrm{d}\mathrm{r}·\mathrm{s}\mathrm{i}\mathrm{n}\mathsf{\psi }\\ {\mathrm{d}\mathrm{T}=\mathrm{d}\mathrm{T}}_{\mathrm{s}}+\mathrm{d}{\mathrm{C}}_{\mathrm{t}}\\ {\mathrm{d}\mathrm{S}=\mathrm{d}\mathrm{S}}_{\mathrm{s}}+\mathrm{d}{\mathrm{N}}_{\mathrm{s}}+\mathrm{d}{\mathrm{g}}_{\mathrm{s}}+\mathrm{d}{\mathrm{C}}_{\mathrm{s}}\\ {\mathrm{d}\mathrm{H}=\mathrm{d}\mathrm{H}}_{\mathrm{H}}+\mathrm{d}{\mathrm{N}}_{\mathrm{H}}+\mathrm{d}{\mathrm{g}}_{\mathrm{H}}+\mathrm{d}{\mathrm{C}}_{\mathrm{H}}\end{array}$$

In the analysis, the blade was divided into several segments and each segment was very short, so it was assumed that the forces were distributed evenly in each segment. Therefore, the hub forces can be expressed as shown in Equation (14). $T,S \mathrm{a}\mathrm{n}\mathrm{d} H$ are the thrust, lateral force and backward force of the hub.

$$\begin{array}{c}T\left(\psi \right)={\displaystyle \sum _{i=1}^N}\underset{e}{\overset{R}{\int }}d{T}_i\left({\psi }_i\right)dr\\ S\left(\psi \right)={\displaystyle \sum _{i=1}^N}\underset{e}{\overset{R}{\int }}d{S}_i\left({\psi }_i\right)dr\\ H\left(\psi \right)={\displaystyle \sum _{i=1}^N}\underset{e}{\overset{R}{\int }}d{H}_i\left({\psi }_i\right)dr\end{array}$$

(14)

$d{T}_i\left({\psi }_i\right)$ is the the microsection thrust of the ith blade at azimuth ${\psi }_i$.${\psi }_i$can be expressed as in Equation (15).

$${\psi }_i=\psi +\frac{2\pi }{N}(i-1)$$

(15)

Reference [15] studied the effect of IBC on the vertical force of the hub. In order to analyze the comprehensive impact of IBC on the thrust, lateral force and backward force of the hub, the vibration load $Fvb$ was defined as in Equation (16).

$$\begin{array}{c}Fvb=Max\left(T\left(\psi \right)\right)-Min\left(T\left(\psi \right)\right)+Max\left(H\left(\psi \right)\right)-Min\left(H\left(\psi \right)\right)\\ +Max\left(S\left(\psi \right)\right)-Min\left(S\left(\psi \right)\right)\end{array}$$

(16)

Fvb0 is the baseline of $Fvb$. It represents the vibration load without applying IBC. Fv, which is the ratio of Fvb to Fvb0, can represent the impact of IBC on the vibration loads of the hub. Fv is given as in Equation (17).

$$Fv=\frac{Fvb}{Fvb0}$$

(17)

This paper determined whether the applied IBC reduced the rotor vibration load by calculating whether Fv was less than 1. The smaller the Fv, the better the vibration reduction efficiencies.

2.3. Optimization Algorithm

Due to the complexity of the calculation of the helicopter hub loads, the steepest descent-golden section combinatorial optimization algorithm was used to find the optimal parameters to reduce rotor vibration. The calculation process of this method is shown in Figure 6 (where $J$ represents the vibration load Fvb). The amplitudes and phases of second and third harmonics were used as the optimization variable, which was defined as in Equation (18).

$${\delta }_f=\left[A_2,{\phi }_2,A_3,{\phi }_3\right]$$

(18)

At the beginning of the calculation, give n = 0 and calculate the gradient $G_{grad}^n$.

$$G_{grad}^n=\nabla J^k\left({\delta }_f\right)=(\frac{\partial F_{vb}}{\partial A_2},\frac{\partial F_{vb}}{\partial {A\phi }_2},\frac{\partial F_{vb}}{\partial A_3},\frac{\partial F_{vb}}{\partial {\phi }_3})$$

(19)

[a, b] are the range of values for variable ${\delta }_f$. The entire calculation process stops when $G_{grad}^n$ meets the accuracy.

2.4. Fuzzy Neural Network Combined with PID Controller

As the actuator to be designed is an electromechanical type, the amplitude of higher order harmonics that can be applied is determined by the height of the actuator’s structural surface, and the phase can be controlled by the motor. Therefore, in the control simulation, the amplitude of the higher harmonic was selected based on the optimization algorithm. In this system, the input ${ r}_{in}\left(k\right)$was the increment of phases of higher order harmonics, and the vibration load $F_v$was used as the system output ${ r}_{out}\left(k\right)$. The incremental PID control algorithm was adapted in the article [29], and can be described as in Equation (20)

$$u\left(k\right)=u\left(k-1\right)+∆u\left(k\right)$$

(20)

where $∆u\left(k\right)=k_p\left(e\left(k\right)-e\left(k-1\right)\right)+k_{i}e\left(k\right)+k_d(e\left(k\right)-2e\left(k-1\right)+e\left(k-2\right))$.

The optimization process is shown in Figure 7.

The paper uses a fuzzy neural network to automatically adjust the three parameters of the PID controller, as shown in Figure 8 (where $k_p$ is proportional coefficient, $k_i$ is integral time constant and $k_d$ is differential time constant).

The second layer in the figure represents the fuzzification layer, which is used to receive numerical values passed down from the input layer and to perform fuzzification processing. The membership of each input variable is defined by the Gaussian function, which is as follows:

$${\mu }_i^j=e^{-\frac{{(x_i-c_{ij})}^2}{{\sigma }_{ij}{}^2}}$$

(21)

where $c_{ij}$ and ${\sigma }_{ij}$ represent the center and width of the Gaussian function.

The third layer is the fuzzy inference layer, which can be written as in Equation (22):

$${\alpha }_j={\mu }_1^{i1}{\mu }_2^{i2}\cdots {\mu }_n^{in}$$

(22)

The function of the fourth layer is to transform dimensional expressions into dimensionless expressions in the inference process of fuzzy neural networks. The expression is as follows:

$${\beta }_j={\alpha }_j/{\displaystyle \sum _j^m{\alpha }_j},j=1,2,\cdots m$$

(23)

The fifth layer is the defuzzification layer, which is used to convert the blurring quantity into an accurate quantity; that is:

$$y_i={\displaystyle \sum _{j=1}^r{w}_{ij}{\beta }_j},j=1,2\cdots r$$

(24)

where $w_{ij}$ represents the connection weight between the third and fourth layer.

In the structure of fuzzy neural networks, the parameters that can be self-adjusted are $c_{ij}$, ${\sigma }_{ij}$ and $w_{ij}$. Assuming that the expected output of the controlled object is $y_k\prime$ and the actual output of the controlled object is $y_k$, the expression for the error performance index is:

$$E=\frac{1}{2}{\displaystyle \sum _{k=1}^{r}e{(k)}^2}=\frac{1}{2}{\displaystyle \sum _{k=1}^r{\left[y(k)-y\prime (k)\right]}^2}$$

(25)

By defining $x_i^{(j)}$ and $y_i^{(j)}$ as the input and output of the ith node in jth layer, the following can be obtained:

$$\begin{array}{l}{\delta }_i^{(5)}=-\frac{\partial E}{\partial f_i^{(5)}}=-\frac{\partial E}{\partial y_i}=y_i^{\prime }-y_i\\ \frac{\partial E}{\partial w_{ij}}=-\frac{\partial E}{\partial f_i^{(5)}}\frac{\partial f_i^{(5)}}{\partial w_{ij}\delta }=-{\delta }_i^{(5)}{y}_j^{(4)}=-(y_i^{\prime }-y_i){\beta }_j\\ {\delta }_j^{(4)}=-\frac{\partial E}{\partial f_i^{(4)}}=-{\displaystyle \sum _{i=1}^r\frac{\partial E}{\partial f_i^{(5)}}}\cdot \frac{\partial f_i^{(5)}}{\partial g_j^{(4)}}\cdot \frac{\partial g_j^{(4)}}{\partial f_i^{(4)}}={\displaystyle \sum _{i=1}^r{\delta }_j^{(5)}{w}_{ij}}\\ {\delta }_j^{(3)}=-\frac{\partial E}{\partial f_i^{(3)}}=-\frac{\partial E}{\partial f_i^{(4)}}\cdot \frac{\partial f_i^{(4)}}{\partial g_j^{(3)}}\cdot \frac{\partial g_j^{(3)}}{\partial f_i^{(3)}}={\delta }_j^{(4)}{\displaystyle \sum _{\begin{array}{l}i=1\\ i\ne j\end{array}}^m{x}_i^{(3)}}/({\displaystyle \sum _{i=1}^m{x}_i^{(3)}}{)}^2\\ ={\delta }_j^{(4)}{\displaystyle \sum _{\begin{array}{l}i=1\\ i\ne j\end{array}}^m{\alpha }_i}/({\displaystyle \sum _{i=1}^m{\alpha }_i}{)}^2\\ {\delta }_j^{(2)}=-\frac{\partial E}{\partial f_j^{(2)}}=-{\displaystyle \sum _{k=1}^m\frac{\partial E}{\partial f_k^{(3)}}\cdot \frac{\partial f_k^{(3)}}{\partial g_{ij}^{(2)}}\cdot \frac{\partial g_{ij}^{(2)}}{\partial f_{ij}^{(2)}}}\\ ={\displaystyle \sum _{k=1}^m{\delta }_j^{(3)}{S}_{ij}{e}^{f_{ij}^{(2)}}}={\displaystyle \sum _{k=1}^m{\delta }_j^{(3)}{S}_{ij}{e}^{-\frac{(x_i-c_{ij})}{{\sigma }_{ij}^2}}}\\ \frac{\partial E}{\partial c_{ij}}=\frac{\partial E}{\partial f_{ij}^{(2)}}\cdot \frac{\partial f_{ij}^{(2)}}{\partial c_{ij}}=-{\delta }_{ij}^{(2)}\cdot \frac{2(x_i-c_{ij})}{{\sigma }_{ij}^2}\\ \frac{\partial E}{\partial {\delta }_{ij}}=\frac{\partial E}{\partial f_{ij}^{(2)}}\cdot \frac{\partial f_{ij}^{(2)}}{\partial {\delta }_{ij}}=-{\delta }_{ij}^{(2)}\cdot \frac{2(x_i-c_{ij})}{{\sigma }_{ij}^3}\end{array}$$

The first step degree of each parameter can be derived from the above formula. Therefore, for the three parameters that need to be adjusted, the expression for their variation is as shown in Equation (16):

$$\begin{array}{c}{w}_{ij}(k+1)=w_{ij}(k)-\eta \frac{\partial E}{\partial w_{ij}}+\lambda \Delta w_{ij}(k) i=1,2,\cdots r j=1,2,\cdots m\\ c_{ij}(k+1)=c_{ij}(k)-\eta \frac{\partial E}{\partial c_{ij}}+\lambda \Delta c_{ij}(k) i=1,2,\cdots n j=1,2,\cdots m_i\\ {\sigma }_{ij}(k+1)={\sigma }_{ij}(k)-\eta \frac{\partial E}{\partial {\sigma }_{ij}}+\lambda \Delta {\sigma }_{ij}(k) i=1,2,\cdots n j=1,2,\cdots m_i\end{array}$$

(26)

where $\Delta {\sigma }_{ij}(k)={\sigma }_{ij}(k)-{\sigma }_{ij}(k-1)$, $\Delta w_{ij}(k)=w_{ij}(k)-w_{ij}(k-1)$, $\Delta c_{ij}(k)=c_{ij}(k)-c_{ij}(k-1)$ where $\eta$ represents the learning rate; $\lambda$ is the inertia coefficient. In this paper, $\eta$ and $\lambda$ are set to 0.1.

3. Results

3.1. Model Validation

In order to verify the correctness of the dynamic model, modal measurement experiments were conducted, as shown in Figure 9. Two rubber ropes with appropriate length and stiffness were chosen to hang the blades at both ends, so that the blades were approximately in a free-free state. Ten excitation points and one measurement point were determined along the spanwise and chordal directions on the blade, a hammer was used to excite each excitation point, and the acceleration sensor was placed at the measurement point to obtain the vibration response of the blade. The flap and lag plane at each excitation point was struck three times with the hammer. The collection and processing of experimental data were based on the platform of LMS Test.Lab. Measurement data were used to calculate the frequency response function and correlation function. The frequency response function is shown in Figure 10, from which the modal frequency of the blade can be obtained. By determining the vibration mode at the corresponding frequency, it can be determined whether the modified modal is the blade elastic vibration modal.. The blade was made of composite, and its shape is shown in Figure 11. The typical profile of the blade is shown in Figure 12. According to the structure of the blade profile, the blade properties of each section were calculated with variational asymptotic beam sectional analysis (VABS). The blade properties are shown in

Table 1. Properties of the blade section.

r/m	EIf/N·kg	EIl/N·kg	GJ/N·m2	M/kg/m	YG/m
0–0.042	164.1	890.1	70.3	0.643	0
0.042–0.082	159.2	744.2	68.4	0.624	0
0.082–0.892	12.8	562.6	7.47	0.283	0
0.892–0.912	13.4	596.9	7.49	0.255	0

, in which r means radial station, EIf is flap stiffness, EIl is lag stiffness, GJ is torsional stiffness, M is blade sectional mass and YG is chordwise blade c.g. location.

The parameters were substituted into the model established in Section 2 to obtain the mode frequencies of the blade. The calculated and experimental values are shown in

Table 2. Mode frequencies.

Mode	Calculation/Hz	Experiment/Hz	Error/%
1st flap	26.21	26.32	−0.42
2nd flap	74.89	75.93	−1.37
3rd flap	150.00	144.07	4.12
1st lag	172.78	180.76	−4.41
2nd lag	490.51	517.54	−5.22
1st torsion	199.52	203.32	−1.87

. The calculating values were in fairly good agreement with the experimental data.

3.2. Optimal Parameters

This paper analyzed the impact of the amplitudes and phases of different order harmonics on the rotor vibration load. The helicopter parameters in reference [11] were used as example parameters. Firstly, the paper investigated the influence of single harmonics on vibration loads. When the helicopter was flying at an advance ratio of 0.26, the influence of phases of second and third harmonics on vibration loads is shown in Figure 13 and Figure 14.

Figure 13 and Figure 14 show that although the amplitudes were different, the trend of the influence of phases on vibration load was similar, but the vibration reduction effect achieved was different. They all decreased, firstly, and then increased with the increases of phases. Only a certain range of phases of IBC could reduce the rotor vibration load. It can be seen from Figure 14 that, under larger amplitudes of third harmonics, IBC with any phases of harmonics could not achieve vibration reduction.

Then, the paper analyzed the impact of amplitudes of different harmonics on vibration loads. The paper calculated the variation trend of the vibration load with the amplitude of each harmonic under a specific phase of different harmonic. The results are shown in Figure 15 and Figure 16.

Figure 15 and Figure 16 show that although the phases were different, the trend of the influence of amplitudes on vibration load was similar. Under different phases, the amplitude of the harmonic varied slightly when the vibration load reached its minimum value.

The previous analysis results showed that it was hard to find the optimal parameters directly with IBC. Therefore, the steepest descent-golden section combinatorial optimization algorithm was used to find the optimal parameters for IBC. It firstly analyzed the single harmonic control law for vibration load reduction. The results of the second harmonic control are shown in Figure 17.

When the helicopter was flying at a low speed, the optimal harmonic amplitude firstly decreased and then increased with the increase in the forward flying speed. When the forward speed was high, the optimal amplitude increased rapidly with the increase in the forward speed. According to the results of References [30,31], when the SA349/2 helicopter was flying at a low speed, the higher the speed, the smaller the in-plane vibration hub forces. Therefore, as the forward flight speed increased, the higher-order harmonic amplitude required for vibration reduction decreased. Vibration hub forces increased significantly when helicopters fly at high speeds, and the required higher-order harmonic amplitude also increased.

The results of the third harmonic control are shown in Figure 18.

The optimal phase of the third harmonic also increased slowly, The difference is that the optimal amplitude of the third harmonic increased with the increase in the forward flying speed both at low and high speeds. According to the results of Reference [16], out-plane hub forces increase as the forward flight speed increases, and the required higher-order harmonic amplitude also increases.

Then, this article analyzed the results of combined second and third harmonic control. They are shown in Figure 19 and Figure 20.

It can be seen from the results that the optimal states parameters of multiple harmonics were not the combination of those obtained directly from single harmonic control. The optimal amplitudes were lower than those of single states, which were easier to put into practice.

3.3. Results under Fuzzy Neural Network PID Control

This article took the helicopter’s forward flight speed of $\mu =0.26$ as the control object. Based on the results of Section 2.2, the amplitude and phase of the second harmonic were set to $0.2°$ and $180°$. The expected output was set to 40% of the initial vibration value; that is, the model helicopter vibration load was set to reduce the vibration reduction effect by 60%, and obtain the hub vibration load through the fuzzy neural network PID control. The change trend of the vibration load ratio is shown in Figure 21.

The changes of PID parameter values are shown in Figure 22. The parameters of PID were $k_p=0.594,k_i=0.709,k_d=0.820$.

The fuzzy neural network PID controller could quickly calculate the manipulation quantity that meets the expected vibration effect based on the minimum performance function. The phase of the second order harmonic is shown in Figure 23.

4. Conclusions

The paper analyzed the impact of different orders, amplitudes, and phases of IBC on vibration loads. The constrained optimization problem was established, with the high-order harmonic-related parameters as design variables, the rotor vibration load coefficient as the objective function, and the harmonic amplitude and phase variation range as constraints. The optimization problem was solved using the steepest descent-golden section combinatorial optimization algorithm. If combined harmonics are used, smaller amplitudes of second and third order harmonics are needed. On this basis, fuzzy neural network PID control was carried out to achieve the rapid control of the vibration load of the hub. According to the results obtained by the algorithms, the vibration level was 60% less than that without the applied IBC.

The paper only theoretically analyzed the impact of IBC on rotor vibration load, and did not conduct relevant experimental research. The paper studied the influence of IBC on the vibration load of the hub, and in the future studying the influence of IBC on the blade profile load, which has a significant impact on blade fatigue, could be considered. The mechanical actuator designed in the article has immutable amplitude, so only the phase of high order harmonics was controlled during fuzzy neural network PID control. With the development of new actuators, it may be necessary to simultaneously control the amplitude and phase of each harmonic.