Research on unmanned aerial vehicle modeling and control based on intelligent algorithms

Abstract

This article designs an automatic flight control system for an unmanned aerial vehicle helicopter. The differential evolution intelligent algorithm is used for a state-space model identification; the differential evolution method has an advantage of choosing initial point randomly. The accuracy of the identified model is verified by comparing the model-predicted responses with the responses collected during flight experiments. The reliability and efficiency of the differential evolution algorithm are demonstrated by the experimental results. A robust controller is designed based on the identified model for the unmanned aerial vehicle helicopter with two-loop control frame: the outer-loop is used to obtain the expected attitude angles through reference path and speed with guidance-based path-following control, and the inner-loop is used to control the attitude angles of helicopter tracking the expected ones with $H_{\infty }$ loop-shaping method. The greatest common right divisor method is used to choose the weighting matrix in loop shaping, in which the stability margin is larger and has a greater bandwidth of the unmanned aerial vehicle system. Finally, a space spiral curve trajectory tracking simulation is conducted to illustrate the efficiency of the proposed control systems, and the simulation results prove that the unmanned helicopter system achieves a top-level control performance.

Introduction

In recent years, unmanned vehicles have attracted a great deal of interest from the university, the industry, and the military world. The past decade has seen a golden age in the development of unmanned aerial vehicle (UAV), but there are only a few documented examples of small-scale unmanned helicopter applications in real-world scenarios; this is mainly due to the poor flight performance that can be achieved and guaranteed under automatic control.

The accurate model and proper control method are the important steps in flight control system design. The UAV helicopter model is a complicated system, and the system identification method is well suited to the rotorcraft problem. The UAV helicopter model is established by first principles. The nonlinear flight dynamic model of the unmanned helicopter is determined by choosing the body velocities, angular rates, and Euler angles as the state vector and the four helicopter control variables as the control input vector. A linear state-space model is derived from a nonlinear model for certain trim points using small disturbance principles. Conventional time-domain techniques are often not well suited to these difficult aspects of helicopter as discussed by Tischler and Ivler ¹ and Hensen and Steinbuch, ² while the frequency-domain methods are well suited to the rotorcraft problem.

Artificial intelligence (AI) is well suited to solving complex problems such as model identification, which is a logical reasoning process with many constrains. At the same time, AI can significantly enhance the intelligence level and autonomous capabilities of UAV, thus improving the level of cluster control. In this article, a nonlinear search based on differential evolution (DE) intelligent algorithm is conducted for a 6-degree-of-freedom (DoF) linear state-space model that matches the frequency-response dataset; the DE method has an advantage of choosing initial point randomly, and the reliability and efficiency of the DE algorithm are demonstrated by the experimental results.

A large set of control methods for UAV helicopter have been reported from classical control to neural-based adaptive control, and several typical methods are introduced as follows. Proportional–integral–derivative (PID) control ³ is the most popular method, but the high coupling of the UAV system cannot be solved. The optimal control method ⁴ which needs accurate model of UAV is not well suited. The neural-based adaptive control model ⁵ costs a long time to train, and the control frame is too complicated. H_∞ loop-shaping method can systematically handle the multivariable and uncertain nature of the aircraft, ⁶ which can solve the channel coupling by loop-shaping method, the accurate model is not needed due to the robustness of the control method, and do not need a long time to design the control system compared to the neural-based control method. The greatest common right divisor (GCRD) method is used to move the transfer function matrix from the real system to the target system, which is a very useful way to solve the difficulties in choosing a proper weighting matrix in loop shaping.

Frequency identification modeling method

Integrated window method

An integrated window method is used to combine the conditioned frequency response, achieving an integrated conditioned frequency response that has good coherence and low random error over the entire frequency range of interest, and the weighting cost function, J, to minimize at each discrete frequency is shown as follows

$$J=\sum _{i=1}^{n_w}{W}_i[{\left(\frac{G_{x{x}_c}-G_{x{x}_i}}{G_{x{x}_i}}\right)}^2+{\left(\frac{G_{y{y}_c}-G_{y{y}_i}}{G_{y{y}_i}}\right)}^2+{\left(\frac{G_{x{y}_c}-G_{x{y}_i}}{G_{x{y}_i}}\right)}^2+5.0{\left(\frac{{\gamma }_{x{y}_c}-{\gamma }_{x{y}_i}}{{\gamma }_{x{y}_i}}\right)}^2]$$

(1)

where $W_i$ is the weighting function, $n_w$ is the number of windows, subscript is the individual window data, and subscript c is the “integrated” window results.

A function that varies inversely with the random error, $({\epsilon }_r{)}_i$ , is

$$W_i={\left[\frac{{({\epsilon }_r)}_i}{{({\epsilon }_r)}_{min}}\right]}^{-4}$$

so that the data with the minimum error $({\epsilon }_r{)}_{min}$ are given a weight of $W_i=1$ and the windows with higher errors are de-weighted accordingly.

The frequency-response function $H(f)$ can be calculated by the following formula

$$H(f)=\frac{G_{x{y}_c}(f)}{G_{x{x}_c}(f)}$$

(2)

The integrated result shown in Figure 1 indicates excellent identification ( ${\gamma }^2>0.8$ ) at a wide frequency range (0.5–30 rad/s), and the frequency responses are shown in Table 1.

Figure 1.

Coherence of lateral channel.

Table 1.

Frequency response and ranges selected for identification.

Out/in	${\delta }_{\mathrm{lat}}$ (rad/s)	${\delta }_{\mathrm{lon}}$ (rad/s)	${\delta }_{\mathrm{ped}}$ (rad/s)	${\delta }_{\mathrm{col}}$ (rad/s)
u	–	0.6–24	–	–
v	0.7–28	–	–	–
w	–	–	–	0.2–4
p	0.5–30	1.0–10	1.2–8.0	–
q	–	0.6–25	–	–
r	2.0–10	–	0.7–18	0.8–3
$a_x$	–	0.6–24	–	–
$a_y$	0.7–28	0.9–12	1.0–10	–
$a_z$	–	–	–	0.2–4

Simplified model based on frequency-response table

The state-space model of UAV helicopter is established, and the stability and control derivatives are identified by matching the conditioned frequency responses

$$\begin{array}{c}\stackrel{·}{x}=\mathrm{Ax}+\mathrm{Bu}\hfill \\ y=\mathrm{Cx}+\mathrm{Du}\hfill \end{array}$$

(3)

where

$$\begin{array}{c}x={\left[\begin{array}{ccccccccc}\Delta u\hfill & \Delta v\hfill & \Delta w\hfill & \Delta p\hfill & \Delta q\hfill & \Delta r\hfill & \Delta \varphi \hfill & \Delta \theta \hfill & \Delta \psi \hfill \end{array}\right]}^T\hfill \\ u={\left[\begin{array}{cccc}\Delta {\delta }_{\mathrm{col}}\hfill & \Delta {\delta }_{\mathrm{ped}}\hfill & \Delta {\delta }_{\mathrm{lat}}\hfill & \Delta {\delta }_{\mathrm{lon}}\hfill \end{array}\right]}^T\hfill \end{array}$$

where u, v, and w are the speed weight of the UAV; p, q, and r are the angular rate weight of the UAV; $\varphi$ , θ, and ψ are the Euler angle of the UAV; and ${\delta }_{\mathrm{col}}$ , ${\delta }_{\mathrm{ped}}$ , ${\delta }_{\mathrm{lat}}$ and ${\delta }_{\mathrm{lon}}$ are the four control inputs of the UAV: collective control input, pedal control input, latitude control input, and longitude control input; the matrices of $\mathbf{A}$ and $\mathbf{B}$ are stability derivatives and control derivatives which are shown as follows

$$\begin{array}{c}\mathbf{A}=\left[\begin{array}{ccccccccc}{X}_u\hfill & X_v\hfill & X_w\hfill & X_p\hfill & X_q-W_0\hfill & X_r+V_0\hfill & 0\hfill & -g\cos {\Theta }_0\hfill & 0\hfill \\ \multicolumn{1}{c}{Y_u}& Y_v\hfill & Y_w\hfill & Y_p+W_0\hfill & Y_q\hfill & Y_r-U_0\hfill & g\cos {\Theta }_0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{Z_u}& Z_v\hfill & Z_w\hfill & Z_p-V_0\hfill & Z_q+U_0\hfill & Z_r\hfill & 0\hfill & -g\sin {\Theta }_0\hfill & 0\hfill \\ \multicolumn{1}{c}{L_u}& L_v\hfill & L_w\hfill & L_p\hfill & L_q\hfill & L_r\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{M_u}& M_v\hfill & M_w\hfill & M_p\hfill & M_q\hfill & M_r\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{N_u}& N_v\hfill & N_w\hfill & N_p\hfill & N_q\hfill & N_r\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 1\hfill & \sin {\Phi }_0\tan {\Theta }_0\hfill & \cos {\Phi }_0\tan {\Theta }_0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & \cos {\Phi }_0\hfill & -\sin {\Phi }_0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & \sin {\Phi }_0\sec {\Theta }_0\hfill & \cos {\Phi }_0\sec {\Theta }_0\hfill & 0\hfill & 0\hfill & 0\hfill \end{array}\right]\hfill \\ \mathbf{B}=\left[\begin{array}{cccc}{X}_{{\delta }_{\mathrm{col}}}\hfill & X_{{\delta }_{\mathrm{lat}}}\hfill & X_{{\delta }_{\mathrm{lon}}}\hfill & X_{{\delta }_{\mathrm{ped}}}\hfill \\ \multicolumn{1}{c}{Y_{{\delta }_{\mathrm{col}}}}& Y_{{\delta }_{\mathrm{lat}}}\hfill & Y_{{\delta }_{\mathrm{lon}}}\hfill & Y_{{\delta }_{\mathrm{ped}}}\hfill \\ \multicolumn{1}{c}{Z_{{\delta }_{\mathrm{col}}}}& Z_{{\delta }_{\mathrm{lat}}}\hfill & Z_{{\delta }_{\mathrm{lon}}}\hfill & Z_{{\delta }_{\mathrm{ped}}}\hfill \\ \multicolumn{1}{c}{L_{{\delta }_{\mathrm{col}}}}& L_{{\delta }_{\mathrm{lat}}}\hfill & L_{{\delta }_{\mathrm{lon}}}\hfill & L_{{\delta }_{\mathrm{ped}}}\hfill \\ \multicolumn{1}{c}{M_{{\delta }_{\mathrm{col}}}}& M_{{\delta }_{\mathrm{lat}}}\hfill & M_{{\delta }_{\mathrm{lon}}}\hfill & M_{{\delta }_{\mathrm{ped}}}\hfill \\ \multicolumn{1}{c}{N_{{\delta }_{\mathrm{col}}}}& N_{{\delta }_{\mathrm{lat}}}\hfill & N_{{\delta }_{\mathrm{lon}}}\hfill & N_{{\delta }_{\mathrm{ped}}}\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill \end{array}\right]\hfill \end{array}$$

The derivatives in the state-space model (3) are very small to be neglectable when the helicopter is in hover condition; if they are not dropped from the model, the identified results will be redundant. We will show how to reduce the model structure based on frequency-response table.

If the off-axis response is poor which means that the energy from input to output is weak, as r is the heading angle rate, ${\delta }_{\mathrm{lon}}$ is the longitude control input, and the energy from the control input ${\delta }_{\mathrm{lon}}$ to output r is weak, then we can see from Table 1 that

$$\frac{r}{{\delta }_{\mathrm{lon}}}(s)\approx 0$$

(4)

From equation (3), we can get equation (5) as follows

$$\begin{array}{c}\multicolumn{1}{c}{\frac{r}{{\delta }_{\mathrm{lon}}}(s)=\left(\frac{1}{s-N_r}\right)\{N_{{\delta }_{\mathrm{lon}}}+[N_p\frac{p}{{\delta }_{\mathrm{lon}}}(s)]+[N_q\frac{q}{{\delta }_{\mathrm{lon}}}(s)]}\\ \multicolumn{1}{c}{+[N_u\frac{u}{{\delta }_{\mathrm{lon}}}(s)]+[N_v\frac{v}{{\delta }_{\mathrm{lon}}}(s)]+[N_w\frac{w}{{\delta }_{\mathrm{lon}}}(s)]\}}\end{array}$$

(5)

Combining equations (4) and (5), if equation (4) is satisfied, we obtain $N_{{\delta }_{\mathrm{lon}}}\approx 0$ , $N_q\approx 0$ and $N_u\approx 0$ .

According to the same analysis, when the helicopter is in hover condition, the frequency-response pairs $q/{\delta }_{\mathrm{col}}$ , $q/{\delta }_{\mathrm{lat}}$ , $u/{\delta }_{\mathrm{lat}}$ , $u/{\delta }_{\mathrm{ped}}$ , $v/{\delta }_{\mathrm{col}}$ , $w/{\delta }_{\mathrm{lon}}$ , $w/{\delta }_{\mathrm{lat}}$ and $w/{\delta }_{\mathrm{ped}}$ have low coherence, so we can conclude that these transfer function pairs are nearly zero. As a result, we can derive the simplified state-space model based on frequency-response table as follows

$$\begin{array}{c}\mathbf{A}=\left[\begin{array}{ccccccccc}{X}_u\hfill & 0\hfill & 0\hfill & X_p\hfill & X_q\hfill & 0\hfill & 0\hfill & -g\cos {\Theta }_0\hfill & 0\hfill \\ \multicolumn{1}{c}{Y_u}& Y_v\hfill & 0\hfill & Y_p\hfill & 0\hfill & 0\hfill & g\cos {\Theta }_0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & Z_w\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill & -g\sin {\Theta }_0\hfill & 0\hfill \\ \multicolumn{1}{c}{L_u}& L_v\hfill & L_w\hfill & L_p\hfill & L_q\hfill & L_r\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{M_u}& 0\hfill & 0\hfill & 0\hfill & M_q\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& N_v\hfill & N_w\hfill & N_p\hfill & 0\hfill & N_r\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 1\hfill & \sin {\Phi }_0\tan {\Theta }_0\hfill & \cos {\Phi }_0\tan {\Theta }_0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & \cos {\Phi }_0\hfill & -\sin {\Phi }_0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill & \sin {\Phi }_0\sec {\Theta }_0\hfill & \cos {\Phi }_0\sec {\Theta }_0\hfill & 0\hfill & 0\hfill & 0\hfill \end{array}\right]\hfill \\ \mathbf{B}=\left[\begin{array}{cccc}{X}_{{\delta }_{\mathrm{col}}}\hfill & 0\hfill & X_{{\delta }_{\mathrm{lon}}}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& Y_{{\delta }_{\mathrm{lat}}}\hfill & Y_{{\delta }_{\mathrm{lon}}}\hfill & Y_{{\delta }_{\mathrm{ped}}}\hfill \\ \multicolumn{1}{c}{Z_{{\delta }_{\mathrm{col}}}}& 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{L_{{\delta }_{\mathrm{col}}}}& L_{{\delta }_{\mathrm{lat}}}\hfill & L_{{\delta }_{\mathrm{lon}}}\hfill & L_{{\delta }_{\mathrm{ped}}}\hfill \\ \multicolumn{1}{c}{M_{{\delta }_{\mathrm{col}}}}& M_{{\delta }_{\mathrm{lat}}}\hfill & M_{{\delta }_{\mathrm{lon}}}\hfill & M_{{\delta }_{\mathrm{ped}}}\hfill \\ \multicolumn{1}{c}{N_{{\delta }_{\mathrm{col}}}}& N_{{\delta }_{\mathrm{lat}}}\hfill & 0\hfill & N_{{\delta }_{\mathrm{ped}}}\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill \end{array}\right]\hfill \end{array}$$

Model identification using intelligent algorithm

The following transfer function is obtained by taking the Laplace transfer of equation (3) ⁷

$$\mathbf{T}(s)=[\mathbf{C}{(s\mathbf{I}-\mathbf{A})}^{-1}\mathbf{B}+\mathbf{D}]$$

(6)

A weighted function J is determined by taking the error between the identification frequency response and the model response, and the unknown parameters in state-space model are obtained through minimizing J over a selected frequency range

$$J=\sum _{l=1}^{n_{\mathrm{TF}}}\left\{\frac{20}{n_{\omega }}\sum _{{\omega }_1}^{{\omega }_{n_w}}{W}_{\gamma }[W_g{(\left|{\hat{T}}_c^l\right|-\left|T^l\right|)}^2+W_p{(\angle {\hat{T}}_c^l-\angle T^l)}^2]\right\}$$

(7)

where $||$ is the magnitude (dB); $\angle$ is the phase (degree); $n_{\omega }$ is the number of frequency points; and ${\omega }_1$ and ${\omega }_{n_{\omega }}$ are the starting and ending frequencies, respectively

$$W_r={[1.58(1-e^{-{\gamma }_{\mathrm{xy}}^2})]}^2,W_g=1.0,W_p=0.01745$$

where W_r is the weighting coefficient of the frequency response according to the coherence value; W_g is the weighting coefficient of the magnitude of the frequency response; and W_p is the weighting coefficient of the phase (degree) of the frequency response, where W_p = pi/180 = 0.01745, changing the degree to radian.

Each input–output pair in good coherence is selected for the identification, and an acceptable level of accuracy requires the cost function $J_{\mathrm{ave}}=J/n_{\mathrm{TF}}\le 100$ . ⁸ An identification vector $\Theta =[\begin{array}{cccc}{\theta }_1\hfill & {\theta }_2\hfill & \cdots \hfill & {\theta }_{\mathrm{np}}\hfill \end{array}]$ is determined by collecting the parameters to be identified in the model matrices. A nonlinear search based on DE algorithm is conducted for the linear state-space model, which has an advantage of choosing initial point randomly. Suppose that $x_i$ is a candidate solution in solution space, the DE’s population of the nth iteration is represented by $\mathbf{X}(n)=[\begin{array}{cccc}{\mathbf{x}}_1\hfill & {\mathbf{x}}_2\hfill & \cdots \hfill & {\mathbf{x}}_{\mathrm{NP}}\hfill \end{array}]$ . Mutation, recombination, and selection are the three evolutionary operations of DE algorithm, which are described as follows. ⁹

Mutation

Mutation is to generate totally different individual randomly by different individuals as given in equation (8), and the mutation process is described in Figure 2

$${\mathbf{v}}_i={\mathbf{x}}_{r_1}+F·({\mathbf{x}}_{r_2}-{\mathbf{x}}_{r_3})$$

(8)

Figure 2.

Mutation process of DE.

Recombination

Recombination operation is to copy individual as given in equation (9), and the recombination process is shown in Figure 3

$$\begin{array}{c}{u}_{i,j}=\hfill \\ \{\begin{array}{c}{v}_{i,j},\hfill \\ \multicolumn{1}{c}{x_{i,j},}\end{array}\begin{array}{c}\mathrm{if}\mathrm{randb}\le C_R\mathrm{or}j=r_j\hfill \\ \multicolumn{1}{c}{\mathrm{if}\mathrm{randb}>C_R\mathrm{or}j\ne r_j}\end{array}(i=1,\dots ,N_p;j=1,\dots ,N_p)\hfill \end{array}$$

(9)

where C_R is a constant in the interval $[0,1]$ and $r_j$ is an integer by $[1,D]$ .

Figure 3.

Recombination process of DE.

Selection

Selection operator is to choose better individuals using equation (10), which makes to converge to the best individuals

$${{\mathbf{x}}^{t+1}}_i=\{\begin{array}{c}{\mathbf{u}}_i^t,\hfill \\ \multicolumn{1}{c}{{\mathbf{x}}_i^t,}\end{array}\begin{array}{c}J({\mathbf{u}}_i^t)<J({\mathbf{x}}_i^t)\hfill \\ \multicolumn{1}{c}{\mathrm{otherwise}}\end{array}$$

(10)

Design of the frequency sweep test

Raptor-50 UAV helicopter with a complete auto-pilot system is used to do frequency-sweep test as shown in Figure 4. Low frequencies are important for the identification of the speed derivatives, and high frequency excitations are important for the identification of the coupled dynamics. To guarantee that the flight data capture the dynamics, a frequency-sweep technique is used, where the pilot gradually increases the frequency of the input from 0.1 to 5 Hz. Figure 5 shows the flight data of input and output by frequency-sweep test. The input is the control input, outputs are the helicopter states, and the sampling rate of the test is 100 Hz.

Figure 4.

Instrumented UAV helicopter in experiment: (a) ground test and (b) frequency-sweep test.

Figure 5.

Input and output by frequency-sweep test.

Identification results and model verification

The frequency responses computed from the flight data (solid) and the identified model are compared in Figure 6, which indicates that identified model and flight data are in good agreement. The model identified should be verified to ensure it can predict the real dynamic response accurately. Figure 7 shows the verification results in the time domain, which indicates that the model identified is well matched.

Figure 6.

Frequency response from the flight data and identified model.

Figure 7.

Identified model verification results.

Robust controller design

$H_{\infty }$ loop-shaping design process

The $H_{\infty }$ loop-shaping standard block diagram is shown in Figure 8, where G is the linear model plant and W₁ and W₂ are the weight matrices.

Figure 8.

$H_{\infty }$ loop-shaping standard block diagram.

$G(s)$ is a left coprime factorization of the linear plant expressed as follows

$$G(s)={\tilde{M}}^{-1}(s)\tilde{N}(s)$$

where ${\tilde{M}}^{-1}(s)$ and $\tilde{N}(s)$ are the stable coprime transfer functions.

${\tilde{M}}^{-1}(s)$ and $\tilde{N}(s)$ are normalized left coprime if they satisfy the following normalization condition

$$\tilde{M}{\tilde{M}}^{*}+\tilde{N}{\tilde{N}}^{*}=I,\forall s\in \mathrm{jR}$$

where ${\tilde{M}}^{*}(s)={\tilde{M}}^T(-s)$ and ${\tilde{N}}^{*}(s)={\tilde{N}}^T(-s)$ .

Coprime factor uncertainty is obtained by introducing perturbations to a left coprime factorization of G

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

Stabilizing controllers, $K_{\infty }$ , are obtained by maximizing the following equation

$$\begin{array}{c}{max}_{\mathrm{stab}{K}_{\infty }}\Vert \left[\begin{array}{c}{w}_1\hfill \\ \multicolumn{1}{c}{w_2}\end{array}\right]\to \left[\begin{array}{c}{z}_1\hfill \\ \multicolumn{1}{c}{z_2}\end{array}\right]{\Vert }^{-1}\hfill \\ =\underset{\mathrm{stab}{K}_{\infty }}{max}{\Vert \left[\begin{array}{c}I\hfill \\ \multicolumn{1}{c}{K_{\infty }}\end{array}\right]{(I-G_s{K}_{\infty })}^{-1}\left[\begin{array}{cc}I\hfill & G_s\hfill \end{array}\right]{\Vert }_{\infty }}^{-1}=\epsilon \hfill \end{array}$$

(11)

where $\epsilon$ is the stability margin in the interval [0, 1], which is the measure of robustness and performance; a margin greater than 0.3 is considered good, based on theory and practical experience, and the stability is guaranteed for

$$\Vert \left[{\tilde{\Delta }}_N{\tilde{\Delta }}_M\right]\Vert <\epsilon$$

New method for weight matrix selection

Design engineers usually obtain good loop-shaping weights and controllers by choosing proper diagonal weights $W_1$ and $W_2$ to achieve a desired performance. However, the design experience shows that the diagonal weights cannot satisfy the requirements very well, especially for plants with strong cross-coupling, such as helicopter. Then, the non-diagonal weights are required to obtain good performance. ¹⁰

This article uses GCRD method to transfer the transfer function matrix from the real system G to the target system $G_s$ and is described below.

Definition 1

Let $G_1,G_2\in \mathit{M}(\wp )$ and have the same number of columns. Then

1. $R\in \mathit{M}(\wp )$ is said to be a common right divisor (CRD) of $G_1$ and $G_2$ if there exist $N_1,N_2\in \mathit{M}(\wp )$ such that $G_1=N_{1}R$ and $G_2=N_{2}R$ .

2. $R\in \mathit{M}(\wp )$ is said to be a GCRD of $G_1$ and $G_2$ if (1) it is a CRD of $G_1$ and $G_2$ and (2) any other CRD of $G_1$ and $G_2$ , say $\overline{R}$ is a right divisor of R, that is, there exists $\overline{N}\in \mathit{M}(\wp )$ such that $R=\overline{N}\overline{R}$ .

Theorem 1

Let $N_{\Omega }\in \mathit{M}(\wp )$ have the following arbitrary partitioning

$$N_{\Omega }=\left(\begin{array}{c}{N}_{1\Omega }\hfill \\ \multicolumn{1}{c}{N_{2\Omega }}\end{array}\right)$$

where $N_{1\Omega }$ is $r_1\times m$ , $N_{2\Omega }$ is $r_2\times m$ and $r_1+r_2=r$ . Then, $N_{\Omega }(s)$ has a UR-decomposition $N_{\Omega }(s)=U_{\Omega }\left(\begin{array}{c}{R}_{\Omega }\hfill \\ \multicolumn{1}{c}{0}\end{array}\right)$ , and furthermore, $R_{\Omega }$ is a GCRD of $N_{1\Omega }$ and $N_{2\Omega }$ ; for the expression of $R_{\Omega }$ and $U_{\Omega }$ , see theorem 2.

Theorem 2

Let $N_{\Omega }(s)$ has full-column rank and the following state-space realization

$$N_{\Omega }(s)=\left[\begin{array}{cc}A\hfill & B\hfill \\ \multicolumn{1}{c}{C}& D\hfill \end{array}\right]$$

where $A\in {\mathbb{C}}^{n\times n}$ , $B\in {\mathbb{C}}^{n\times m}$ and $C\in {\mathbb{C}}^{r\times n}$ .

In the following, we define

$$\begin{array}{c}\multicolumn{1}{c}{\left[\begin{array}{cc}\tilde{A}\hfill & \tilde{B}\hfill \\ \multicolumn{1}{c}{\tilde{C}}& \tilde{D}\hfill \end{array}\right]=\left[\begin{array}{cc}\tilde{T}(A-B{D}^{*}C){\tilde{T}}^{-1}\hfill & \tilde{T}B\hfill \\ \multicolumn{1}{c}{D_{\perp }{D}_{\perp }^{*}C{\tilde{T}}^{-1}}& D\hfill \end{array}\right]}\\ \multicolumn{1}{c}{=\left[\begin{array}{ccccc}{\tilde{A}}_{00}\hfill & 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{{\tilde{A}}_{10}}& {\tilde{A}}_{11}\hfill & 0\hfill & 0\hfill & {\tilde{B}}_1\hfill \\ \multicolumn{1}{c}{{\tilde{A}}_{20}}& 0\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{{\tilde{A}}_{30}}& {\tilde{A}}_{31}\hfill & {\tilde{A}}_{32}\hfill & {\tilde{A}}_{33}\hfill & {\tilde{B}}_3\hfill \\ \multicolumn{1}{c}{{\tilde{C}}_0}& {\tilde{C}}_1\hfill & 0\hfill & 0\hfill & D\hfill \end{array}\right]\begin{array}{c}{\tilde{n}}_0\hfill \\ \multicolumn{1}{c}{{\tilde{n}}_1}\\ \multicolumn{1}{c}{{\tilde{n}}_2}\\ \multicolumn{1}{c}{{\tilde{n}}_3}\end{array}}\end{array}$$

(12)

where $D_{\perp }\in {\mathbb{C}}^{r\times (r-m)}$ such that $[\begin{array}{cc}D\hfill & D_{\perp }\hfill \end{array}]$ is unitary and $\tilde{T}$ is an observable controllable canonical transformation matrix chosen so that $(\begin{array}{cccc}{\tilde{A}}_{11},\hfill & {\tilde{B}}_1,\hfill & {\tilde{C}}_1,\hfill & D\hfill \end{array})$ is a minimal realization of $(\begin{array}{cccc}{\tilde{A}}_{11},\hfill & {\tilde{B}}_1,\hfill & {\tilde{C}}_1,\hfill & D\hfill \end{array})$ .

Then, we can get the minimal realization of $R_{\Omega }$ and $U_{\Omega }$ as follows

$$U_{\Omega }=\left[\begin{array}{ccc}{\overline{A}}_{11}-{\tilde{B}}_1\tilde{F}\hfill & {\tilde{B}}_1\hfill & {\tilde{H}}_1{D}_{\perp }\hfill \\ \multicolumn{1}{c}{{\tilde{C}}_1-D{\tilde{F}}_1}& D\hfill & D_{\perp }\hfill \end{array}\right]$$

(13)

$${3\mathrm{ptR}}_{\Omega }(s)=\left[\begin{array}{cc}A\hfill & B\hfill \\ \multicolumn{1}{c}{{\tilde{F}}_1{\tilde{T}}_1+D^{*}C}& I\hfill \end{array}\right]$$

(14)

where $\tilde{F}$ and ${\tilde{H}}_1$ are, respectively, any $m\times {\tilde{n}}_1$ state feedback and ${\tilde{n}}_1\times r$ output injection which make $(-\mathrm{Is}+{\tilde{A}}_{11}-{\tilde{B}}_1{\tilde{F}}_1{)}^{-1}$ and $(-\mathrm{Is}+{\tilde{A}}_{11}-{\tilde{H}}_1{\tilde{C}}_1{)}^{-1}$ $\Omega$ -stable.

Next, we will use the GCRD method to transfer the transfer function matrix from the real system G to the target system $G_s$ . First, we get the target system $G_s(s)=\mathrm{diag}[\begin{array}{cccc}7/s\hfill & 7.5/s\hfill & 9/s\hfill & 8/s\hfill \end{array}]$ based on the requirement that the crossover frequencies of the open-loop system should not be less than 7 rad/s.

Let $N_{1\Omega }$ and $N_{2\Omega }$ of theorem 1 be the real system G and the target system $G_s$ in the $H_{\infty }$ loop shaping; then, we can get the transfer matrix W based on equation (3).

Simulation of control performance

Simulation of inner-loop control performance

The traditional $H_{\infty }$ loop-shaping method and the GCRD method are used to compare the control performance of the UAV system, and the Bode plots for UAV single loop are shown in Figures 9 and 10.

Figure 9.

Bode plots for UAV single loop (traditional method): (a) w channel, (b) r channel, (c) φ channel and (d) θ channel.

Figure 10.

Bode plots for UAV single loop (GCRD method): (a) w channel, (b) r channel, (c) φ channel and (d) θ channel.

When the loop-shaping design was finished, equation (13) was maximized to achieve the 24-order controller $K_{\infty }$ , and the stability margin $\epsilon =0.6957$ . The Bode plots with gain and phase margins are shown in Figure 9, and the robustness specifications of all the loops are satisfied; moreover, the loop has the greater stability margin compared with the traditional method shown in Figure 10.

The simulation results proved that the stability performance of the unmanned helicopter system achieves a top-level performance in accordance with the ADS-33E-PRF. ¹¹ Table 2 compares bandwidth ${\omega }_{\mathrm{BW}}$ which is obtained from Bode plots of Figures 9 and 10; it shows that the new method is better than the traditional method in decoupling performance. Figure 11 shows the responses for each of the four axes when sharp pulses are applied directly at correspondent plant’s input; the response time of the system with the new method is faster, and the oscillatory is smaller.

Table 2.

Bandwidth evaluations with ADS-33E.

Channel	${\omega }_{\mathrm{BW}}$ (rad/s) ADS-33E	${\omega }_{\mathrm{BW}}$ (traditional method)	${\omega }_{\mathrm{BW}}$ (new method)
$\theta$	≥2.5	18.21	>30
$\phi$	≥2.5	13.24	>30
r	≥3.5	13.98	>30

Figure 11.

Normalized responses to pulse disturbance: (a) traditional method and (b) GCRD method.

Simulation of outer-loop control performance

The inner-loop closure makes available a high-bandwidth, robust, and decoupled system with four inputs ( $w_r$ , $r_r$ , ${\varphi }_r$ , ${\theta }_r$ ) and four outputs (w, r, $\varphi$ , $\theta$ ). The decoupling justifies the use of cascaded single-input, single-output (SISO) outer loops to allow the helicopter to fly a specified trajectory. This article uses guidance-based path-following method ¹² to design the outer loop. Figure 12 shows the simulation results of the robotic helicopter system that tracks the climbing helix trajectory: the center of the circle is (100, 0), the radius is 100 m, the expected climbing height is 150 m, the expected climbing speed is 1 m/s, and the expected forward speed is 10 m/s; adding the random wind disturbance in the lateral direction with the amplitude ±2 m/s, the simulation results show that the control result is very effective.

Figure 12.

Simulation results of the climbing helix trajectory tracking.

Conclusion

This article designed an automatic control system for a UAV helicopter based on intelligent algorithms; the frequency identification method is used to extract the UAV helicopter model, and DE intelligent algorithm is used to get the unknown matrix derivatives of the UAV model, using the $H_{\infty }$ loop-shaping method to design the inner loop as to satisfy the requirements of control performance in the ADS33E. The GCRD method is used to choose the weighting matrix of the UAV control system which exhibits a larger robust stability margin. The simulation results of a space spiral curve trajectory tracking proved the efficiency of the proposed control systems.