Nomenclature
-
$\tau_{m}$
-
response time constant of the primary rotor
-
$a_{1}$
-
primary rotor longitudinal flapping coefficient
-
$b_{1}$
-
primary rotor lateral flapping coefficient
-
$A_{1}$
-
primary rotor lateral periodic pitch input
-
$B_{1}$
-
primary rotor longitudinal periodic pitch input
-
$\tau_{s}$
-
stable aileron response time constant
-
$c_{1}$
-
lateral flap coefficient of the aileron
-
$d_{1}$
-
stabilizing aileron longitudinal flapping coefficient
-
$L_{b1}$
-
flapping moment derivative of the lateral primary rotor
-
$M_{a1}$
-
flapping moment derivative of the longitudinal primary rotor
-
$\tau_{\text{sp}}$
-
lateral stability aileron response time constant
-
$\tau_{sq}$
-
longitudinal stability aileron response time constant
-
$D_{\alpha}$
-
best positions of
$\alpha$
wolf -
$D_{\beta}$
-
best positions of
$\beta$
wolf -
$D_{\delta}$
-
best positions of
$\delta$
wolf -
$\alpha_{l}$
-
initial value of
$\alpha$
-
$R^{2}$
-
Coefficient of determination
-
${\hat{y}}_{i}$
-
predicted value
-
$y_{i}$
-
actual value
-
${\overset{\overline{}}{y}}_{i}$
-
average value
- WLS
-
Weighted Least Squares
- IGWO
-
Improved Grey Wolf Optimization
- WLS-IGWO
-
Weighted Least Squares-Improved Grey Wolf Optimization
1.0 Introduction
A small unmanned helicopter, due to its abilities to hover, perform vertical take-offs and landings, and cruise at low speeds, as well as its advantages of flexibility, ease of control and low cost, is widely used in both military and civilian fields. However, due to the nonlinear dynamic characteristics of unmanned helicopters and the influence of various uncertainties in flight, system identification and modeling are important research topics in unmanned helicopter control. An accurate dynamics model of small unmanned helicopters is key to flight control and autonomous flight. However, due to the complex mechanical structure of small unmanned helicopters, such as the rotor system, there are many difficulties in modeling small unmanned helicopters. Therefore, it is difficult to establish a dynamic model of small unmanned helicopters quickly, accurately and at low cost.
The dynamic characteristics of a small unmanned helicopter are complex, featuring underactuation, strong coupling and multivariable properties, making it crucial to establish an accurate dynamic model. The dynamic modeling methods of small unmanned helicopters can be mainly divided into two types: mechanism-based modeling and system identification. Mechanism-based modeling is based on disciplines such as aerodynamics and flight mechanics, obtaining a more accurate mathematical model through rigorous theoretical calculations, which requires extensive professional knowledge and experience. System identification can effectively deal with the complex structural characteristics of the system by analysing the linear relationship between the helicopter input and output data and establishing an equivalent input and output model [Reference Leshikar, Valasek and McQuinn1, Reference Hosseini, Steinert, Hofmann, Fang, Steffensen, Holzapfel and Göttlicher2]. For example, the study on the physical modeling and simulation verification of small fixed-wing unmanned aircraft provides a reliable control foundation by accurately describing the flight state [Reference Venkataraman and Seiler3, Reference Kapeel, Safwat, Kamel, Khalil, Elhalwagy and Hendy4]. Experimental results indicate that, whether it is a small helicopter, a multirotor unmanned aerial vehicle, or a fixed-wing unmanned aircraft, system identification and modeling play a critical role in formulating and optimising flight control strategies [Reference Dayhoum, Zakaria and Abdelhamid5, Reference Tischler, Scepanovic and Gubbels6]. Therefore, system identification is an important means of understanding the dynamic characteristics of an aircraft, providing an accurate model foundation for the design and optimisation of controllers. However, traditional system identification methods often suffer from low accuracy, slow convergence and sensitivity to noise when dealing with the complexity of unmanned helicopters.
In recent years, researchers have adopted various methods for system identification and modeling of multirotor UAVs and single-rotor helicopters. For instance, Hoshu et al. [Reference Hoshu, Wang, Ansari, Sattar and Bilal7] and Niki et al. [Reference Maleki, Karimi, Mohammadi and Ashenayi8] studied system identification methods for multirotor UAVs, which achieved certain improvements in identification accuracy, but their models struggled to handle complex nonlinear dynamic characteristics. In the study of small helicopters, an improved internal model control (IMC) strategy has also been used for yaw control, enhancing helicopter manoeuverability through precise system identification; however, its performance under data noise interference still requires improvement [Reference Huang, Pei and Cheng9, Reference Xia10]. Steen [Reference Steen11] and Tischler et al. [Reference Tischler12] based their work on frequency domain system identification methods, which demonstrated good stability, though their convergence speed and robustness were limited in scenarios with rapid dynamic changes. System identification based on the prediction error method has been widely used in multirotor UAV modeling, effectively improving flight control stability and accuracy [Reference Hoshu, Wang, Ansari, Sattar and Bilal7, Reference Maleki, Karimi, Mohammadi and Ashenayi8, Reference Steen11]. Geluardi et al. [Reference Geluardi, Nieuwenhuizen, Venrooij, Pollini and Bülthoff13] performed system identification for the hovering state of the Robinson R44 helicopter, using frequency domain methods to model and analyse its dynamic characteristics. Khadeeja Nusrath et al. [Reference TK, VP and Singh14] studied the system identification of swashplateless rotor UAVs; although the method successfully modeled the dynamic characteristics of such rotors, its accuracy and stability were insufficient in situations with significant data noise and high model complexity. Bauer and Nagy [Reference Bauer and Nagy15] researched high-fidelity system identification methods based on flight data, demonstrating good identification capabilities. Jianhong and Ramirez-Mendoza [Reference Jianhong and Ramirez-Mendoza16] proposed a cascaded estimation method, providing a new theoretical perspective for system identification; however, in practical applications, they still face challenges brought by complex environments. Further research indicates that existing methods still need improvement in their identification effectiveness when dealing with the strong nonlinearity and parameter uncertainties of unmanned helicopters [Reference Juhasz, Celi and Tischler17, Reference Babcock18]. Moreover, data-driven methods such as deep learning and neural networks have gradually been applied to UAV system identification. For example, research on rotation modeling and control strategies for small, unmanned helicopters based on deep learning has demonstrated the advantages of deep learning in handling complex nonlinear systems [Reference Xia10]. Methods combining the prediction error approach and deep learning, when dealing with multivariable and strongly coupled nonlinear systems, have high model complexity and computational costs, and are prone to getting stuck in local optima. Although domain adversarial neural networks have alleviated this issue to some extent, their robustness in the face of uncertainties and noise interference still needs improvement [Reference Chen, Zhang, Wang, Yu, Wang and Chen19].
The application of optimisation algorithms in system identification has also received widespread attention. Among them, the grey wolf optimisation (GWO) algorithm, with its strong global search capability and convergence speed, has demonstrated outstanding performance in nonlinear system optimisation. GWO simulates the hunting behaviour of grey wolves, achieving a global optimal solution through the interaction between individuals. However, the traditional GWO algorithm, when solving multi-objective optimisation problems, tends to fall into premature convergence, leading to inaccurate identification results [Reference Liu, Wei and Huang20–Reference Nadimi-Shahraki, Taghian and Mirjalili24]. Therefore, improved versions of the GWO algorithm have been proposed. For instance, Altay et al. [Reference Luo, He and Gao25] and Luo et al. [Reference Ou, Yin and Mo26] combined fusion strategies and multilayer perceptron neural networks to further enhance the accuracy and speed of identification. The improved GWO algorithm proposed by Qiu et al. [Reference Qiu, Yang and Chen27] demonstrated superior performance in functional optimisation and engineering design problems, with good global search capability and convergence speed. Altay et al. [Reference Altay and Altay28] also combined the improved GWO with multilayer perceptron neural networks, further improving the optimisation effect, especially in nonlinear system modeling applications. Although existing methods have achieved some success in the system identification of unmanned helicopters, there are still areas that need improvement.
Therefore, this paper proposes a new identification algorithm combining the weighted least squares method and the improved grey wolf optimisation algorithm (WLS-IGWO) for system identification of small unmanned helicopters. The weighted least squares (WLS) method, as a linear least squares algorithm, has the advantages of being simple and highly flexible [Reference Chen, Zhang, Wang, Yu, Wang and Chen19], while the GWO algorithm is a heuristic optimisation algorithm inspired by the collective behaviour of grey wolves in nature. The improved grey wolf optimisation (IGWO) algorithm features faster convergence speed and better global search capability, avoiding getting trapped in local optima and enhancing the robustness of the algorithm [Reference Nadimi-Shahraki, Taghian and Mirjalili24–Reference Ou, Yin and Mo26]. In the WLS-IGWO algorithm, the WLS method is used to fit the frequency domain response and obtain the initial model structure, followed by parameter optimisation of the initial model using the IGWO algorithm. By introducing nonlinear convergence factors, improving grey wolf grouping rules, and dynamic collaboration rules based on the basic GWO algorithm, the global search capability of the algorithm is enhanced, overcoming the tendency of the basic GWO algorithm to fall into local optima. The WLS-IGWO algorithm can ensure identification accuracy while reducing computational complexity, thereby providing an effective new approach for precise modeling and control of small unmanned helicopters. Compared with the traditional identification methods, the proposed method shows significant advantages in model fitting accuracy, response speed and anti-interference ability, which not only provides more reliable model support for the control and navigation of unmanned helicopters, but also has wide applicability, providing a new idea for the modeling and control of similar aircraft, robots or other nonlinear systems.
2.0 Small, unmanned helicopter mathematical model
This article uses the blade element method to conduct aerodynamic analysis on the primary rotor and stability aileron of a small, unmanned helicopter [Reference Tischler, Scepanovic and Gubbels6]. The dynamic characteristics of a small, unmanned helicopter are affected by the flapping motion of the primary rotor and stabiliser ailerons. The primary rotor flapping equation of a small, unmanned helicopter is:
\begin{align}{\tau _m}{{\dot {a}}_1} = - {a_1} - {\tau _m}q + {A_b}{b_1} + {B_1}\\[-26pt] \nonumber \end{align}
\begin{align}{\tau _m}{\dot b_1} = - {b_1} - {\tau _m}p - {B_a}{a_1} - {A_1}\end{align}
Where:
${\tau _s}$
is the response time constant of the primary rotor,
${a_1}$
is the primary rotor longitudinal flapping coefficient,
${b_1}$
is the primary rotor lateral flapping coefficient;
${A_1}$
is the primary rotor lateral periodic pitch input,
${B_1}$
is the primary rotor longitudinal periodic pitch input.
The stable ailerons of the small, unmanned helicopter are connected to the rotor system through the main shaft, and the stable ailerons will produce a gyroscopic effect to maintain a stable state when rotating. Therefore, the stabilised aileron can be thought of as a lift-less rotor aircraft. The flapping equation of the small, unmanned helicopter’s stabilised aileron is:
\begin{align}{\tau _s}{\dot {c}_1} = - {c_1} - {\tau _s}q + {C_1} \\[-26pt] \nonumber\end{align}
\begin{align}{\tau _s}{\dot {d}_1} = - {d_1} - {\tau _s}p + {D_1}\end{align}
Where:
${\tau _s} = \frac{{16}}{{{\gamma _s}{\rm{\Omega }}}}$
,
${C_1} = {C_{lon}}{\delta _{lon}}$
,
${D_1} = {D_{lat}}{\delta _{lat}}$
.
${\tau _s}$
is the stable aileron response time constant,
${c_1}$
is the lateral flap coefficient of the aileron,
${d_1}$
is the stabilising aileron longitudinal flapping coefficient.
The lateral channel flapping motion equation for the small, unmanned helicopter is as follows when the stable input for aileron control is substituted into the primary rotor flapping equation:
\begin{align}{\tau _m}{\dot {b}_1} = - {b_1} - {\tau _m}p + {B_d}{d_1} + {B_{lat}}{\delta _{lat}}\end{align}
In the steady state
${\dot {b}_1} = 0$
, we get:
\begin{align}{b_1} = - {\tau _m}p + {B_d}{d_1} + {B_{lat}}{\delta _{lat}}\end{align}
After computing the Laplace transform of the stabilised aileron lateral channel’s flapping equation, the transfer function is as follows:
\begin{align}{d_1}\!\left( s \right) = \frac{{ - {\tau _s}p\!\left( s \right) + {D_{lat}}{\delta _{lat}}\!\left( s \right)}}{{{\tau _s}s + 1}}\end{align}
Substituting Equation (7) into Equation (6) we get:
\begin{align}{b_1}\!\left( s \right) = - {\tau _m}p\!\left( s \right) + {B_{{d_1}}}\frac{{ - \tau p\!\left( s \right) + {D_{lat}}{\delta _{lat}}}}{{{\tau _s}s + 1}} + {B_{lat}}{\delta _{lat}}\!\left( s \right)\end{align}
To obtain the lateral angular rate equation, compute and arrange the inverse Laplace transform of using Equation (7):
\begin{align}{b_1}\!\left( s \right) = \frac{{ - \!\left( {{\tau _m} + {B_{{d_1}}}{\tau _s}} \right)p\!\left( s \right)}}{{s + 1/{\tau _s}}} + \frac{{\!\left( {{B_{lat}} + {B_{lat}}{D_{lat}}} \right){\delta _{lat}}\!\left( s \right)}}{{s + 1/{\tau _s}}}\end{align}
The lateral rotational motion equation of the small, unmanned helicopter is
$\dot {p} = {L_{{b_1}}}{b_1}$
, perform Laplace transformation and substitute Equation (8) into it, and the small, unmanned helicopter’s lateral angular speed transfer function model is gained:
\begin{align}\frac{{p\!\left( s \right)}}{{{\delta _{lat}}\!\left( s \right)}} = \frac{{{L_{{b_1}}}\!\left( {{B_{lat}} + {B_{{d_1}}}{D_{lat}}} \right)/{\tau _s}}}{{{s^2} + \!\left( {1/{\tau _s}} \right)s + {L_{{b_1}}}\!\left( {{\tau _m} + {B_{{d_1}}}{\tau _s}} \right)/{\tau _s}}}\end{align}
In the same way, adding the stabilising aileron control input, the small, unmanned helicopter’s longitudinal channel flapping equation is:
\begin{align}{\tau _m}{\dot {a}_1} = - {a_1} - {\tau _m}q + {A_c}{c_1} + {A_{lon}}{\delta _{lon}}\end{align}
In the steady state
${\dot {a}_1} = 0$
,we get:
\begin{align}{\tau _s}{\dot c_1} = - {c_1} - {\tau _s}q + {C_{lon}}{\delta _{lat}}\end{align}
The longitudinal rotational motion equation of a small, unmanned helicopter is:
$\dot q = {M_a}{a_1}$
. The small, unmanned helicopter’s longitudinal channel angular rate transfer function model is as follows, similar to the lateral channel method:
\begin{align}\frac{{q\!\left( s \right)}}{{{\delta _{lon}}\!\left( s \right)}} = \frac{{{M_{{a_1}}}\!\left( {{A_{lon}} + {A_{{c_1}}}{C_{lon}}} \right)/{\tau _s}}}{{{s^2} + \!\left( {1/{\tau _s}} \right)s + {M_{{a_1}}}\!\left( {{\tau _m} + {A_{{c_1}}}{\tau _s}} \right)/{\tau _s}}}\end{align}
Since small, unmanned helicopters performing real flight data acquisition have rudder dynamics, the system identification structure is analysed by substituting the executive rudder dynamics into the transfer function model. The dynamic characteristics of the unmanned helicopter rudder studied in this paper can be represented as a second order system:
\begin{align}{H_{servo}} = \frac{{\omega _s^2}}{{{s^2} + 2{\zeta _s}{\omega _s}s + \omega _s^2}}\end{align}
Therefore, ignoring the cross-coupling effects of flapping motion and considering the characteristics of the steering gear, the transfer functions of the small, unmanned helicopter’s lateral and longitudinal channels are:
\begin{align}\frac{{p\!\left( s \right)}}{{{\delta _{lat}}\!\left( s \right)}} = \frac{{{L_{{b_1}}}\!\left( {{B_{lat}} + {B_{{d_1}}}{D_{lat}}} \right)/{\tau _{sp}}}}{{{s^2} + \!\left( {1/{\tau _{sp}}} \right)s + {L_{{b_1}}}\!\left( {{\tau _{mp}} + {B_{{d_1}}}{\tau _{sp}}} \right)/{\tau _{sp}}}} \times \frac{{\omega _s^2}}{{{s^2} + 2{\zeta _s}{\omega _s}s + \omega _s^2}}\end{align}
\begin{align}\frac{{q\!\left( s \right)}}{{{\delta _{lon}}\!\left( s \right)}} = \frac{{{M_{{a_1}}}\!\left( {{A_{lon}} + {A_{{c_1}}}{C_{lon}}} \right)/{\tau _{sp}}}}{{{s^2} + \!\left( {1/{\tau _{sp}}} \right)s + {M_{{a_1}}}\!\left( {{\tau _{mp}} + {A_{{c_1}}}{\tau _{sp}}} \right)/{\tau _{sp}}}} \times \frac{{\omega _s^2}}{{{s^2} + 2{\zeta _s}{\omega _s}s + \omega _s^2}}\end{align}
Where:
$\;{L_{{b_1}}}$
is the flapping moment derivative of the lateral primary rotor,
${M_{{a_1}}}$
is the flapping moment derivative of the longitudinal primary rotor;
${\tau _{sp}}{\rm{\;}}$
is the lateral stability aileron response time constant,
${\tau _{sq}}{\rm{\;}}$
is the longitudinal stability aileron response time constant. The aforementioned four parameters must be examined and determined via the system identification approach since they cannot be determined theoretically.
3.0 Weighted least squares and improved grey wolf optimisation identification algorithm
3.1 Weighted least squares initial model
Through the frequency response curves of the small, unmanned helicopter’s lateral and longitudinal channels, this paper uses the WLS to fit its frequency response curves to obtain the initial model structure.
Assume that a small, unmanned helicopter’s the transfer function has the following form:
\begin{align}G\!\left( s \right) = \frac{{{b_m}{s^m} + {b_{m - 1}}{{(s)}^{m - 1}} + \cdots + {b_1}s + {b_0}}}{{{a_n}{s^n} + {a_{n - 1}}{s^{n - 1}} + \cdots + {a_1}s + 1}}\end{align}
Its frequency characteristics are:
\begin{align}G\!\left( {j\omega } \right) = \frac{{\!\left( {{b_0} - {b_2}{\omega ^2} + \cdots } \right) + j\omega \!\left( {{b_1} - {b_3}{\omega ^2} + \cdots } \right)}}{{\!\left( {1 - {a_2}{\omega ^2} + \cdots } \right) + j\omega \!\left( {{a_1} - {a_3}{\omega ^2} + \cdots } \right)}} = \frac{{\alpha \!\left( \omega \right) + j\omega \beta \!\left( \omega \right)}}{{\sigma \!\left( \omega \right) + j\omega \eta \!\left( \omega \right)}} = \frac{{N\!\left( {j\omega } \right)}}{{D\!\left( {j\omega } \right)}}\end{align}
Where:
\begin{align}\alpha \!\left( \omega \right) = {b_0} - {b_2}{\omega ^2} + \ldots, \beta \!\left( \omega \right) = {b_1} - {b_3}{\omega ^2} + \ldots, \sigma \!\left( \omega \right) = 1 - {a_2}{\omega ^2} + \ldots, \eta \!\left( \omega \right) = {a_1} - {a_3}{\omega ^2} + \ldots \end{align}
Decompose the system frequency characteristics into real frequency characteristics and imaginary frequency characteristics, then:
\begin{align}{G^{\rm{*}}}\!\left( {j{\omega _i}} \right) = A\!\left( {{\omega _i}} \right){e^{j\varphi \!\left( {{\omega _i}} \right)}} = Re\left( {{\omega _i}} \right) + jIm\!\left( {{\omega _i}} \right){\rm{\;\;\;\;\;\;\;\;}}i = 1,2, \ldots, L\end{align}
Where:
$\;R\!\left( {{\omega _i}} \right) = A\!\left( {{\omega _i}} \right){\rm{cos}}\varphi \!\left( {{\omega _i}} \right)$
;
$Im\!\left( {{\omega _i}} \right) = A\!\left( {{\omega _i}} \right){\rm{sin}}\phi \!\left( {{\omega _i}} \right)$
Among them,
$\;A\!\left( {{\omega _i}} \right)$
and
$\varphi \!\left( {{\omega _i}} \right)$
are the amplitude ratio and phase difference of the measured output curve and input curve, respectively, and L is the number of frequency points taken in the experiment.
Assuming at the frequency point
${\omega _i}$
, the error between the measured frequency characteristics and the estimated frequency characteristics is:
\begin{align}\varepsilon \!\left( {{\omega _i}} \right) = \left[ {Re\!\left( {{\omega _i}} \right) + jIm\!\left( {{\omega _i}} \right)} \right] - \frac{{N\!\left( {j\omega } \right)}}{{D\!\left( {j\omega } \right)}}\end{align}
Since it is a nonlinear minimisation problem, it is difficult to minimise its error. For this reason, the following error criterion is introduced:
\begin{align}J = \mathop \sum \limits_{i = 1}^L {\| { D\!\left( {j{\omega _i}} \right)\varepsilon \!\left( {j{\omega _i}} \right)} \| }^{2} = \mathop \sum \limits_{i = 1}^L \| {D\!\left( {j{\omega _i}} \right)\left[ {Re\!\left( {{\omega _i}} \right) + jIm\!\left( {{\omega _i}} \right)} \right] - N\!\left( {j{\omega _i}} \right)}\| ^{2}\end{align}
Construct a linear minimum problem by transforming the initial nonlinear minimum problem. Substitute Equation (18) into Equation (22) to get:
\begin{align} J = \mathop \sum \limits_{i = 1}^L &\| {\sigma \!\left( {{\omega _i}} \right) + j{\omega _i}\eta \!\left( {{\omega _i}} \right)\left[ {Re\!\left( {{\omega _i}} \right) + jIm\!\left( {{\omega _i}} \right)\left] - \right[\alpha \!\left( {{\omega _i}} \right) + j\omega \beta \!\left( {{\omega _i}} \right)} \right]} \| {^2} \nonumber \\& {} = \mathop \sum \limits_{i = 1}^L \left[ {\sigma \!\left( {{\omega _i}} \right)Re\!\left( {{\omega _i}} \right) + jIm\!\left( {{\omega _i}} \right)RE\!\left( {{\omega _i}} \right) - {\omega _i}\eta \!\left( {{\omega _i}} \right)Im\!\left( {{\omega _i}} \right) - \alpha \!\left( {{\omega _i}} \right)} \right] \end{align}
In order to minimise J, Let J take the respective derivatives of parameters
${a_k}$
and
${b_k}$
, and let them be 0, get:
\begin{align} {{\partial J} \over {{\partial _1}}}{|_{{a_1} = {{\hat \alpha }_1}}} = & \mathop \sum \limits_{i = 1}^L 2\left[ {\sigma \!\left( {{\omega _i}} \right)\!Re\!\left( {{\omega _i}} \right) - {\omega _i}\eta \!\left( {{\omega _i}} \right)\!Im\!\left( {{\omega _i}} \right) - \alpha \!\left( {{\omega _i}} \right)} \right]\left[ { - {\omega _i}Im\!\left( {{\omega _i}} \right)} \right]|{a_1} = \hat \alpha {}\nonumber\\& + \mathop \sum \limits_{i = 1}^L 2\left[ {{\omega _i}\eta \!\left( {{\omega _i}} \right)\!RE\!\left( {{\omega _i}} \right) + \sigma \!\left( {{\omega _i}} \right)\!Im\!\left( {{\omega _i}} \right) - {\omega _i}\beta \!\left( {{\omega _i}} \right)} \right]{\omega _i}Re\!\left( {{\omega _i}} \right){|_{{a_1} = {{\hat \alpha }_1}}} = 0 \end{align}
\begin{align} {{\partial J} \over {\partial {a_2}}}{|_{{a_2} = {{\hat \alpha }_2}}} = {}&\mathop \sum \limits_{i = 1}^L 2\left[ {\sigma \!\left( {{\omega _i}} \right)\!Re\!\left( {{\omega _i}} \right) - {\omega _i}\eta \!\left( {{\omega _i}} \right)\!Im\!\left( {{\omega _i}} \right) - \alpha \!\left( {{\omega _i}} \right)} \right]\!Re\!\left( {{\omega _i}} \right)\!\left( { - \omega _i^2} \right)|{a_2} = {{\hat \alpha }_2} \nonumber\\& {} + \mathop \sum \limits_{i = 1}^L 2\left[ {{\omega _i}\eta \!\left( {{\omega _i}} \right)\!RE\!\left( {{\omega _i}} \right) + \sigma \!\left( {{\omega _i}} \right)\!Im\!\left( {{\omega _i}} \right) - {\omega _i}\beta \!\left( {{\omega _i}} \right)} \right]\!Im{\omega _i}\!\left( { - \omega _i^2} \right){|_{{a_1} = \hat \alpha }} = 0 \end{align}
$\vdots $
\begin{align}{{\partial J} \over {\partial {b_1}}}{|_{{b_1} = {{\hat b}_1}}} = \mathop \sum \limits_{i = 1}^L 2\left[ {\sigma \!\left( {{\omega _i}} \right)\!Re\!\left( {{\omega _i}} \right) - {\omega _i}\eta \!\left( {{\omega _i}} \right)\!Im\!\left( {{\omega _i}} \right) - \alpha \!\left( {{\omega _i}} \right)} \right]\!\left( { - 1} \right){|_{{b_1} = {{\hat b}_1}}} = 0\end{align}
\begin{align}{{\partial J} \over {\partial {b_2}}}{|_{{b_2} = {{\hat b}_2}}} = \mathop \sum \limits_{i = 1}^L 2\left[ {{\omega _i}\eta \!\left( {{\omega _i}} \right)\!Re\!\left( {{\omega _i}} \right) + \sigma \!\left( {{\omega _i}} \right)\!Im\!\left( {{\omega _i}} \right) - {\omega _i}\beta \!\left( {{\omega _i}} \right)} \right]\!\left( {{ - _i}} \right){|_{{b_2} = {{\hat b}_2}}} = 0\end{align}
$\vdots $
Substituting Equation (19) into the above formula, one may consider the small, unmanned aircraft utilised in this article’s transfer function model to be a fourth-order system. Assume that the function of transfer is:
\begin{align}G\!\left( s \right) = \frac{{{b_0}}}{{{a_4}{s^4} + {a_3}{s^3} + {a_2}s_ + ^2{a_1}s + 1}}\end{align}
To determine the WLS method’s initial model structure, solve the following:
\begin{align}\left[ \begin{array}{c@{\quad}c@{\quad}c@{\quad}c@{\quad}c}{{W_0}} {}&{{T_1}} {}&{{S_2}} {}&{ - {T_3}} {}&{{S_4}}\\{{T_1}} {}&{{R_2}} {}&0 {}&{ - {R_4}} {}&0\\{{S_2}} {}&0 {}&{{R_4}} {}&0 {}&{ - {R_6}}\\{{T_3}} {}&{{R_4}} {}&0 {}&{ - {R_6}} {}&0\\{{S_4}} {}&0 {}&{{R_6}} {}&0 {}&{ - {R_8}}\end{array} \right]\left[ {\begin{array}{*{20}{l}}{{b_0}}\\{{a_1}}\\{{a_2}}\\{{a_3}}\\{{a_4}}\end{array}} \right] = \left[ {\begin{array}{*{20}{l}}{{S_0}}\\0\\{{R_2}}\\0\\{{R_4}}\end{array}} \right]\end{align}
3.2 Grey wolf optimisation algorithm
3.2.1 Basic grey wolf optimisation algorithm
The motivation for using the improved GWO algorithm in the system identification of small, unmanned helicopters lies in its strong global search capability, adaptability to complex nonlinear systems, avoidance of local optimality, fast convergence speed and strong algorithm robustness. The GWO algorithm can effectively solve the problem of dynamic parameter identification of small, unmanned helicopter systems, especially in noisy and uncertain environments, and can improve the accuracy and robustness of model identification, meeting the requirements of real-time performance and computational efficiency.
GWO algorithm is a population-based algorithm. Different from other population algorithms, the GWO algorithm has a hierarchical structure and is able to establish the hierarchical advantages of the objective function [Reference Xiaolin, Gang and Guohua21].
In the GWO algorithm,
$\alpha $
,
$beta$
and
$\delta $
are used to represent the three wolves in the initial search process, which have the highest priority in decision-making. In addition, the other dominant one is
$\omega $
wolf, which does not participate in any decision-making. In the grey wolf algorithm, only the alpha wolf can identify the prey location and guide the path planning.The principle of the grey wolf algorithm is shown in Fig. 1.
Figure 1. Grey wolf optimisation algorithm principle.
From a behavioural point of view, the GWO algorithm simulates the hunting process of wolves collectively searching for prey.Therefore, it is assumed that the places of
$\alpha $
wolf,
$\beta $
wolf and
$\delta $
wolf are the best positions, and their mathematical expression is as shown in Equation (30):
\begin{align}{D_\alpha } = \left| {{C_1}{X_\alpha } - X} \right|\end{align}
\begin{align}{D_\beta } = \left| {{C_1}{X_\beta } - X} \right|\end{align}
\begin{align}{D_\delta } = \left| {{C_1}{X_\delta } - X} \right|\end{align}
Where:
${D_\alpha }$
,
${D_\beta }$
and
${D_\delta }$
the best positions of
$\alpha $
wolf,
$\beta $
wolf, and
$\delta $
wolf, respectively. The average distance between them can be obtained by
$\omega $
wolf, and its mathematical equation is as follows:
\begin{align}{X_1} = {X_\alpha } - {A_1}{D_\alpha }\end{align}
\begin{align}{X_2} = {X_\beta } - {A_2}{D_\beta }\end{align}
\begin{align}{X_3} = {X_\delta } - {A_3}{D_\delta }\end{align}
\begin{align}X = \frac{{{X_1} + {X_2} + {X_3}}}{3}\end{align}
Where:
${X_\alpha }$
represents the position of
$\alpha $
,
${X_\beta }$
represents the position of
$\beta $
,
${X_\delta }$
represents the position of
$\delta $
,
${A_1}$
,
${A_2}$
and
${A_3}$
represent the position adjustment coefficients, forming the vector A,
${C_1}$
,
${C_2}$
and
${C_3}$
form vector C, X represents the places of the current solution.
Vector A and vector C make the GWO algorithm a random algorithm. Since the vector has a specified range, the vector fluctuation within this range can avoid local minima. The expressions of vector A and vector C are as follows:
\begin{align}\left\{ {\begin{array}{*{20}{l}}{A = 2{\alpha _1} - \alpha } {}{}\\{\alpha = 2\!\left( {1 - \frac{t}{T}} \right)} \end{array}} \right.\\[-24pt] \nonumber\end{align}
\begin{align}C = 2{r_2}\end{align}
Where: t is the number of iterations, T is the maximum number of iterations,
${r_1}$
and
${r_2}$
are random vectors between [0,1],
${\rm{\;}}\alpha $
is linearly reduced from 2 to 0 during the iteration process, so the value range of
$\alpha $
is set to [−2, 2].
3.2.2 Improved GWO algorithm
In contrast to the basic GWO algorithm, the IGWO has stronger local search and global optimisation capabilities. The local search capability makes the algorithm converge faster, and the global optimisation capability makes the population richer, its performance and efficiency are higher, and it improves the problem-solving ability of the algorithm [Reference Li, Lin and Liu22, Reference Ou, Yin and Mo26].
To improve the accuracy of the algorithm, the weighted least squares is combined with the improved grey wolf optimisation algorithm (WLS-IGWO), and to improve the basic GWO algorithm as follows:
-
1. Nonlinear convergence factor
In the GWO algorithm, under the action of the convergence factor, its change interval is [−a, a]. By controlling the change of A, the search range of the grey wolf algorithm can be adjusted. If
$\left| A \right| \gt 1$
, then the grey wolves have a larger search range and strong global search capabilities. If
$\left| A \right| \lt 1$
, the search range of the wolf pack is concentrated locally, which is suitable for local refined optimisation. However, in the changing trend of
$\alpha $
linear decrease, to avoid the GWO algorithm slipping into global optimisation and local optimisation,
$\alpha $
nonlinear improvement is made:(39)
\begin{align}\alpha = {\alpha _l} \times \left[ {[1 - {{\!\left( {\frac{{{2^{t/T}} - 1}}{{e - 1}}} \right)}^k}} \right]\end{align}
Where:
${\alpha _l}$
is the initial value of
$\alpha $
, k is the control factor that determines the attenuation amplitude of
$\alpha $
, the value range is [Reference Leshikar, Valasek and McQuinn1, Reference Xia10]. The size of k is negatively related to the decay speed of
$\alpha $
. -
2. Improve grey wolf grouping rules
In the standard GWO algorithm, all individuals participate in the optimisation process during the calculation process, and the entire population will move in the direction of pack
$\alpha $
. But when the position of a is
$\alpha $
local optimum, the algorithm will fall into a local optimum. Improved as shown in Fig. 2, the wolf pack is split into chasing group and bounding group. The chasing group continuously participates in iterations in the algorithm, while the bounding group does not participate in iterations. After each iteration has been performed, the least adapted individual in the chasing group withdraws from it and participates in the enclosing group, and the most adapted individual in the enclosing swarm participates in the catching up swarm to start the iteration of the algorithm. -
3. Dynamic collaboration rules
In the GWO algorithm, the location update of individuals in the wolf group is carried out under the joint action of the positions of
$\alpha $
wolves,
$\beta $
wolves and
$\delta $
wolves. However, it lacks information exchange between individuals in the population, which may also lead to the GWO algorithm local optimum. Therefore, this paper uses dynamic weighting rules to increase the influence of adjacent individuals in the wolf pack, while maintaining the guiding role of
$\alpha $
wolves on grey wolves in the population.(40)
\begin{align}{C_1} = \frac{{ran\!\left( {0,1} \right)}}{2} \\[-24pt] \nonumber\end{align}
(41)
\begin{align}{C_2} = \frac{{\!\left( {1 - ran\!\left( {0,1} \right)} \right)}}{2} \\[-24pt] \nonumber \end{align}
(42)
\begin{align}{C_3} = 0.5 \\[-24pt] \nonumber \end{align}
(43)
\begin{align}X\!\left( {t + 1} \right) = {C_3} \times {X_1} + {C_2} \times {X_2} + {C_1} \times {X_3} \end{align}
Where:
${C_1}$
,
${C_2}$
and
${C_3}$
are distance adjustment parameters, respectively, where
${C_3}$
is a constant, ran(0,1) represents a random number between 0 and 1. In the dynamic weighting rules, the alpha wolf, as the alpha wolf, should maintain its guiding role in the population, so its weight is the largest. For
$\beta $
wolves and
$\delta $
wolves, in order to increase the information exchange and population diversity of different individuals in the grey wolf population, random inertia weights are added to avoid local optimality and enhance the global optimisation capability of the algorithm.
3.3 Algorithm flow diagram
The IGWO algorithm is used to optimise the parameters in the initial model Equation (29) obtained by the (WLS) method. Calculate the global optimal parameters and substitute the optimal parameters into Equation (29) to obtain the transfer function model, and use the mean square error of the transfer function prediction result as the fitness value.
The primarily steps to improve the GWO algorithm to optimise WLS model parameters are:
Step 1: The initial population number is m = 50, the maximum number of iterations is T = 500. The upper bound of the search range is
${u_b}$
, and the lower bound is
${l_b}$
. The mean square error of the predicted results of the transfer function is taken as the fitness value, and the MSE calculation formula is:
$\;{e_{{\rm{MS}}}} = {1 \over N}\mathop \sum \nolimits_{i = 1}^N {\!\left( {y\!\left( {{x_i}} \right) - f\!\left( {{x_i}} \right)} \right)^2}$
.
Where: N is the number of data;
$f\!\left( {{x_i}} \right)$
is the model prediction data;
$y\!\left( {{x_i}} \right)$
is the actual output data.
Step 2: Build an initial model and import parameters into the initial model.
Step 3: Generate a random wolf group according to the parameters, and calculate the fitness value of the chasing group.
Step 4: Use the improved grey wolf grouping rule to get the wolves
$\alpha $
,
$\beta $
,
$\delta $
with the best fitness value and the individual with the worst fitness value. Initialise the wolves surrounding the group to get the individual
$\mu $
with the best fitness value. Use individual
$\mu $
instead of chasing the wolf with the worst fitness value in the group.
Step 5: Calculate the update parameter
$\alpha $
, add random inertia weight, update parameters A and C; update
$\omega $
wolves.
Step 6: If the termination condition is met, it ends and the optimal solution of the
$\alpha $
wolf position and corresponding parameters is output; if it is not met, it returns to continue the optimisation iteration.
The WLS-IGWO identification process is shown in Fig. 3:
4.0 System identification and model verification
4.1 Flight test data
This article uses the independently developed Raptor-50 small, unmanned helicopter system as the flight test verification platform, as shown in Fig. 4. In order to obtain identification data, the required flight verification data is obtained by scanning the lateral channel and longitudinal channel of the small single-rotor unmanned helicopter to input excitation signals.
4.2 Data preprocessing
Raw flight data acquisition will receive external environmental interference and flight environment and other unfavorable factors interference, so that the raw flight verification data generated in the systematic error and random error. Therefore, data preprocessing is performed on the raw flight data collected prior to validation. In this paper, the detrending term and the rejection and correction of wild values are used to eliminate the bias generated by the sensor when acquiring data and to improve the confidence of the data, and in order to reduce the influence of interference signals, processing techniques include data smoothing and low-pass filtering. Taking the transverse channel as an example, the flight data acquisition input, output and pre-processing results of the test data applied to the system identification are shown in the Fig. 5 and Fig. 6.
Figure 2. Principle of improved grey wolf optimisation algorithm grouping rules.
Figure 3. WLS-IGWO algorithm flow chart.
Figure 4. Small, unmanned helicopter.
Figure 5. Input data preprocessing.
As can be seen from Fig. 5 and Fig. 6 above, the input and output data of small, unmanned helicopter can be preprocessed effectively to reduce noise interference and random errors in flight test data, make the data smoother, and enhance the robustness of the system.
4.3 Result verification and analysis
Select a set of frequency sweep flight verification data for the lateral and longitudinal channels, respectively. After calculating the initial model structure through the WLS method, and optimisation of WLS initial model parameters using IGWO algorithm. The results are shown in Table 1:
Substituting the parameter results in Table 1 into Equation (29), the lateral channel transfer function model is obtained:
\begin{align}G\!\left( s \right) = \frac{{ - 6.9715 \times {{10}^7}}}{{{s^4} - 54.1734{s^3} + 4669{s^2} + 60249s + 193891}}\end{align}
The longitudinal channel transfer function model is:
\begin{align}G\!\left( s \right) = \frac{{1.7476 \times {{10}^7}}}{{{s^4} - 6.567{s^3} + 6541{s^2} - 15568s + 284218}}\end{align}
The transfer function is solved to obtain a set of more similar characteristic roots for the small, unmanned helicopter actuated servos. For constant values of the parameters related to the small, unmanned helicopter’s structure in Equations (44) and (45), Table 2 displays the results of the identification.
To verify the validity of the model, the traditional method, WLS-GWO and WLS-IGWO algorithms are used to compare the model verification results of small, unmanned helicopters in lateral and longitudinal channels as well as the model identification errors through the actual flight data, and the results of the comparative verification are shown in Figs 7–12.
Where Ipwm is the PWM control input, p and q are the roll angle rate and pitch angle rate, respectively. It can be seen from the Figs 7–12, that the WLS-IGWO method better reflects the dynamic characteristics of the system among the model predictions obtained by the three methods.
By reducing the amplitude and phase difference between the model structure and the associated frequency domain response, the conventional approach finds the unknown parameters in the model structure [Reference Geluardi, Nieuwenhuizen, Venrooij, Pollini and Bülthoff13, Reference TK, VP and Singh14]. The traditional method is an algorithm widely used and effective in identifying small, unmanned helicopter systems.
In order to verify the adaptability and accuracy of the model, ten sets of verification samples for all of the horizontal channel and the longitudinal channel were chosen using the traditional method, WLS-GWO and WLS-IGWO methods to verify and compare the identification accuracy of the algorithm, then determine the lateral channel and the longitudinal channel’s coefficient of determination
${R^2}$
. The coefficient of determination is used to evaluate the quality of the model. The closer the value is to 1, the higher the identification accuracy of the model; conversely, the closer the value is to 0, the worse the identification accuracy is
\begin{align}{R^2} = 1 - {{\mathop \sum \limits_{i = 1}^n {{({y_i} - y_i^\prime)}^2}} \over {\mathop \sum \limits_{i = 1}^n {{({y_i} - \bar y)}^2}}}\end{align}
where
${\hat y_i}$
is the predicted value;
${y_i}$
is the actual value;
${\bar y_i}$
is the average value.
Table 1. Results of parameter optimisation
Table 2. Values of model parameters
Figure 6. Output data preprocessing.
It can be seen from Fig. 13 that under ten different sets of verification samples for the lal channel, the R value predicted by the WLS-IGWO method is 10% higher on average than the WLS-GWO method, and 11% higher than the traditional method.
Figure 7. Lateral channel control input.
Figure 8. Lateral channel angular rate model verification.
Figure 9. Lateral channel identification error.
Figure 10. Longitudinal channel control input.
Figure 11. Longitudinal channel angular rate model verification.
Figure 12. Longitudinal channel identification error.
Figure 13. Lateral channel prediction accuracy analysis.
It can be seen from Fig. 14 that under ten different sets of verification samples for the longitudinal channel, the predicted output value of the WLS-IGWO method is improved by an average of 10% compared with the WLS-GWO method and by an average of 12% compared with the traditional method.
Figure 14. Longitudinal channel prediction accuracy analysis.
In response to changes in complex data, the weighted least squares method improves the accuracy of parameter estimation by adding appropriate weights during the process of fitting the frequency domain response curve. Combined with the GWO algorithm and inspired by natural organisms, the IGWO algorithm has strong local search and global optimisation capabilities, the algorithm converges quickly, and the algorithm has certain adaptability advantages. Therefore, applying the WLS-IGWO algorithm to the system identification problem of establishing a flight dynamics model for small, unmanned helicopters can better adapt to complex practical problems.
From the above analysis, it can be seen that the transfer function model obtained from the identification of small, unmanned helicopter systems based on the WLS-IGWO algorithm can better predict the relationship between input and output, has higher accuracy, better model performance and more realistically reflects the small, unmanned helicopter system dynamic characteristics of human helicopter.
5.0 Conclusion
Considering the issue of accurately identifying the frequency domain prediction response curve fitting model for small, unmanned helicopters, a new algorithm using the WLS-IGWO approach to identify model parameters is proposed. Initially, the fitting model’s starting parameters are obtained using the WLS method, and then the IGWO algorithm is used to optimise the initial parameters. Combined with IGWO’s strong local search and global optimisation capabilities and fast convergence speed to improve model accuracy and reduce model identification errors, the lateral and longitudinal channels of the small, unmanned helicopter’s dynamic model was acquired by means of computational optimisation. Through flight test data sample verification, the WLS-IGWO method, WLS-GWO method and traditional method were used to compare the model prediction results. The comparison results show that the WLS-IGWO algorithm identification model results proposed in this article are more accurate in predicting the output results of small, unmanned helicopters. It has certain accuracy and adaptability in establishing dynamic models of small, unmanned helicopters, and provides an application reference for system identification problems in dynamic modeling of small, unmanned helicopters.
However, system identification methods are usually computationally heavy, especially when applied to complex multidimensional nonlinear systems, and there are still some shortcomings in real-time identification and control performance. Therefore, the next step is to study the online system identification algorithm for adaptive control, so that the unmanned helicopter can update the model parameters in real time during flight, and enhance the response and control ability of the system in a dynamic environment. With the development of unmanned aerial vehicle group cooperative flight task, future research can extend the cooperative identification of multiple unmanned aerial systems to support system identification in swarm intelligence and cooperative control.
Acknowledgments
Thanks to all the co-authors for their contributions to the manuscript writing, review and revision of this paper and the collation of experimental data.
Funding
The author declares that this manuscript was supported by the China Scholarship Council.
