System Identification and $H_\infty$-based Control of Quadrotor Attitude

TL;DR

This paper develops an $H_ abla$ control approach for quadrotor attitude stabilization, using system identification to handle uncertainties and demonstrating improved robustness and tracking performance over traditional controllers through real experiments.

Contribution

It introduces a novel $H_ abla$ control method based on system identification for quadrotors, enhancing robustness against uncertainties.

Findings

The $H_ abla$ controller outperforms PID and $mbda$ controllers in robustness.

Experimental results confirm improved tracking accuracy.

The approach effectively manages unmodelled dynamics and parameter uncertainties.

Abstract

The attitude control of a quadrotor is a fundamental problem, which has a pivotal role in a quadrotor stabilization and control. What makes this problem more challenging is the presence of uncertainty such as unmodelled dynamics and unknown parameters. In this paper, to cope with uncertainty, an $H_\infty$ control approach is adopted for a real quadrotor. To achieve $H_\infty$ controller, first a continuous-time system identification is performed on the experimental data to encapsulate a nominal model of the system as well as a multiplicative uncertainty. By this means, $H_\infty$ controllers for both roll and pitch angles are synthesized. To verify the effectiveness of the proposed controllers, some real experiments and simulations are carried out. Results verify that the designed controller does retain robust stability, and provide a better tracking performance in comparison with a well-tuned PID and a $\mu$ synthesis controller.

Figures

Tables

Table 1. Table 1 : List of symbols in Eq. 2

Symbol	Description
$I_{xx,yy,zz}$	Inertia moments
$J_r$	Rotor inertia
${\Omega }_i$	Propeller angular rate
$d$	Drag factor
$b$	Thrust factor
$l$	Horizontal distance: propeller center to CoG*
Center of Gravity

Equations

$$\mathcal{G}=\left\{G\left(s\right)\mid G=\left(1+\Delta W\right)G_0\right\},$$

$$\ddot{\varphi }=\frac{J_r\dot{\theta }\left({\Omega }_1+{\Omega }_3-{\Omega }_2-{\Omega }_4\right)}{I_{xx}}+\frac{I_{yy}-I_{zz}}{I_{xx}}\dot{\psi }\dot{\theta }+\frac{bl\left({\Omega }_2^2-{\Omega }_4^2\right)}{I_{xx}},$$

$$\ddot{\theta }=\frac{J_r\dot{\varphi }\left(-{\Omega }_1-{\Omega }_3+{\Omega }_2+{\Omega }_4\right)}{I_{yy}}+\frac{I_{zz}-I_{xx}}{I_{yy}}\dot{\psi }\dot{\varphi }+\frac{bl\left({\Omega }_3^2-{\Omega }_1^2\right)}{I_{yy}},$$

$$\ddot{\psi }=\frac{d\left({\Omega }_1^2+{\Omega }_3^2-{\Omega }_2^2-{\Omega }_4^2\right)}{I_{zz}}+\frac{I_{xx}-I_{yy}}{I_{zz}}\dot{\theta }\dot{\varphi },$$

$$\ddot{\varphi }=\frac{bl\left({\Omega }_2^2-{\Omega }_4^2\right)}{I_{xx}},$$

$$\ddot{\theta }=\frac{bl\left({\Omega }_3^2-{\Omega }_1^2\right)}{I_{yy}},$$

$$\ddot{\psi }=\frac{d\left({\Omega }_1^2+{\Omega }_3^2-{\Omega }_2^2-{\Omega }_4^2\right)}{I_{zz}}.$$

$$\varphi \left(s\right)=\frac{T_1^{2}bl}{s^2{\left(s+T_2\right)}^2{I}_{xx}}\left(u_4^2\left(s\right)-u_2^2\left(s\right)\right),$$

$$\theta \left(s\right)=\frac{T_1^{2}bl}{s^2{\left(s+T_2\right)}^2{I}_{yy}}\left(u_3^2\left(s\right)-u_1^2\left(s\right)\right),$$

$$\psi \left(s\right)=\frac{T_1^{2}d}{s^2{\left(s+T_2\right)}^2{I}_{zz}}\left(u_1^2\left(s\right)+u_3^2\left(s\right)-u_2^2\left(s\right)-u_4^2\left(s\right)\right),$$

$$\varphi \left(s\right)=\frac{A_1}{s^2}{U}_1;U_1=u_4^2\left(s\right)-u_2^2\left(s\right),$$

$$\theta \left(s\right)=\frac{A_2}{s^2}{U}_2;U_2=u_3^2\left(s\right)-u_1^2\left(s\right),$$

$$\psi \left(s\right)=\frac{A_3}{s^2}{U}_3;U_3=u_1^2\left(s\right)+u_3^2\left(s\right)-u_2^2\left(s\right)-u_4^2\left(s\right),$$

$$G_{roll}\left(s\right)=\frac{A_1\left(k_{d_r}{s}^2+k_{p_r}s+k_{i_r}\right)}{s^3+A_1\left(k_{d_r}{s}^2+k_{p_r}s+k_{i_r}\right)}{U}_1,$$

$$G_{pitch}\left(s\right)=\frac{A_2\left(k_{d_p}{s}^2+k_{p_p}s+k_{i_p}\right)}{s^3+A_2\left(k_{d_p}{s}^2+k_{p_p}s+k_{i_p}\right)}{U}_2,$$

$$G_{pitch}\left(s\right)=\frac{1547.4}{\left(s+5.373\right)\left(s^2+10.12s+390.4\right)},$$

$$G_{roll}\left(s\right)=\frac{2049.8}{\left(s+6.764\right)\left(s^2+19.03s+426.2\right)}.$$

$$|\frac{G\left(jw\right)}{G_0\left(jw\right)}-1|\le |W\left(jw\right)|\forall w.$$

$$W_{pitch}=\frac{1659.6\left(s^2+2.868s+60.44\right)}{\left(s+2.477e04\right)\left(s+9.678\right)},$$

$$W_{roll}=\frac{1.9017\left(s^2+3.813s+91.61\right)}{s^2+43.53s+545.3}.$$

$${\gamma }_{opt}=\underset{C}{\min }\Vert T_z{y}_d\Vert =\underset{C}{\min }{\Vert \begin{array}{c}{W}_{s}S\\ W_{u}U\\ WT\end{array}\Vert }_{\infty }.$$

$$T_{cl}=\frac{{\omega }_n^2}{s^2+2\zeta {\omega }_{n}s+{\omega }_n^2},$$

$$t_s\approx \frac{4}{\zeta {\omega }_n}=0.3,$$

$$M_p=e^{\frac{-\pi \zeta }{\sqrt{1-{\zeta }^2}}}=10^{-6},0<\zeta <1.$$

$$T_{id}=\frac{247.3}{s^2+30.67s+247.3},$$

$$W_s=a\frac{1}{S_{id}}=a\frac{1}{1-T_{id}}=a\frac{s^2+30.67s+247.3}{s\left(s+30.67\right)},$$

$$W_s=a\frac{\left(s^2+30.67s+247.3\right)}{\left(s+1000\right)\left(s+30.67\right)\left(s+0.001\right)}.$$

$$C_{pitch}\left(s\right)=\frac{3.3227e05\left(s+2.477e04\right)\left(s+461.2\right)\left(s+50.77\right)\left(s+9.678\right)\left(s+5.373\right)\left(s^2+10.12s+390.4\right)}{\left(s+2.477e04\right)\left(s+1.673e04\right)\left(s+1000\right)\left(s+30.67\right)\left(s+8.082\right)\left(s+0.001\right)\left(s^2+80.13s+3556\right)}.$$

$$C_{roll}\left(s\right)=\frac{4.3173e06\left(s+371.5\right)\left(s+59.89\right)\left(s+6.764\right)\left(s^2+43.53s+545.3\right)\left(s^2+19.03s+426.2\right)}{\left(s+2.175e05\right)\left(s+1000\right)\left(s+30.67\right)\left(s+27.95\right)\left(s+19.41\right)\left(s+0.001\right)\left(s^2+80.99s+3749\right)}.$$

$$C_{pitch}=\frac{3.3227e05(s+461.5)(s+45.87)(s+5.903)(s^2+9.836s+388.1)}{(s+1.673e04)(s+1000)(s+25.13)(s+0.001)(s^2+79.66s+3577)},$$

$$C_{roll}=\frac{4.3173e06(s+371.8)(s+53.33)(s+7.035)(s^2+18.34s+428.1)}{(s+2.175e05)(s+1000)(s+27.87)(s+0.001)(s^2+80.84s+3820)}.$$

$$K_p=2.6,K_i=0.2,K_d=0.65$$

$$Roll:$$

$$Pitch:$$

Full text

1 Introduction

Quadrotors are widely considered due to their remarkable features and applications in a variety of areas such as monitoring [1, 2], surveillance [3, 4], search and rescue [5, 6], agriculture [7, 8], and delivery of goods [9, 10]. Quadrotors are also of paramount interest owing to the advantages of low cost, maneuverable, small size, simple design structure. The quadrotor is a nonlinear system with coupled states. In general, there are six states including the robot position $\left(x,y,z\right)$ and the angles (yaw, pitch, roll) with only four inputs, resulting in an underactuated system. Despite the mentioned advantages of quadrotors, the open-loop instability, coupled states, and underactuation make their control challenging. A quadrotor control is accompanied with several other serious challenges, particularly when it turns from a simple model in the theory to a practical problem in which many other key factors such as parameters uncertainty, unmodelled dynamics, and input constraints need to be considered. The main focus of this paper is the stabilization and attitude control of a quadrotor considering such uncertainties and input constraints.

Modeling and control of a quadrotor have been extensively studied in the literature. PID control is the most common controller widely used in commercial quadrotoros [11, 12, 13, 14]. In [12] three different PID structures are employed for the attitude control. In [13] a robust PID controller is designed using affine parameterization. In [14] a cascade control for the attitude control is employed and a robust compensator is added to reduce the wind disturbance effect. Linear Quadratic Regulator (LQR) is another linear-based controller used for quadrotors [15, 16]. In LQR controllers, setting proper weights is difficult and needs trial and error. Linear-based controllers have simple designs and have been tested successfully on experimental and simulated platforms. Linear controllers, however, are based on a linear approximation of the actual system model by neglecting the nonlinear behavior of the system which may increase the probability of the system failure. Moreover, they are limited to low velocities as well as small angles of deviation, and their stability is guaranteed only near the selected operating point.

To transcend the limitations of the linear controllers, many nonlinear control approaches have been proposed. In [17], a lyapunov-based controller is designed using singular perturbation theory for a quadrotor stabilization. Backstepping-based control has also been developed by many researchers [18, 19, 20, 21]. A backstepping controller using a lagrangian form of the quadrotor dynamics is designed in [19]. In [20] an optimal backstepping controller using $\text{H}_{\infty}$ is proposed. Sliding mode control is another nonlinear control method which is capable of dealing with some classes of uncertainties and external disturbances [22, 23, 24, 25, 26]. In [22, 23], a combination of the sliding mode control and backstepping control techniques are used for the control of quadrotors. Despite the mentioned capabilities, sliding mode control suffers from chattering, and it is necessary to know the upper bounds of uncertainties.

In order to reduce the chattering effect, different approaches like fuzzy gain scheduling [24] and higher order sliding mode have been proposed[25]. A fuzzy state observer is proposed in [27] to estimate the unknown nonlinear functions of the uncertain system model, then an integral terminal sliding mode controller is used. Many other robust nonlinear control techniques have also been proposed to cope with the quadrotor uncertainties. In [28], a robust compensating technique is used to control quadrotors with time-varying uncertainties and delays. The design of a nonlinear $\text{H}_{\infty}$ controller is studied in [29] to achieve the robustness. In [30], a nonlinear $\text{H}_{\infty}$ output feedback controller coupled with a high order sliding mode estimator is utilized to control the quadrotor in the presence of parameters uncertainty and external disturbances. The adaptive control is another approach for controlling systems with unknown model parameters and have been widely used for controlling various practical systems[31, 32, 33]. An adaptive controller is employed in [34] for trajectory tracking of a quadrotor with considering the input constraints and uncertain parameters in its nonlinear model.

In order to design a nonlinear-based controller, the quadrotor nonlinear model is required, which is hard to obtain. For instance, inertia matrix, center of gravity and some other parameters of a quadrotor are not easily accessible and it is difficult to obtain their accurate values. In most of the previous studies, a symmetric model for the quadrotor is assumed, however in real robots this assumption does not hold. Furthermore, some factors, like wind, cause environmental disturbances. Especially in an actual flight near obstacles or low attitude, undesirable wind effects on the stability and performance of a quadrotor is more obvious. Considering the wind model makes designing a nonlinear controller arduous and challenging. An approach to overcome the model uncertainty, such as unmodelled dynamics and parameter uncertainty, is using system identification techniques. In [35], a closed-loop multivariable extremum seeking algorithm (MESA) is suggested for parameter identification of a quadrotor. The black-box model identification using a continuous-time approach is adopted in [36] to obtain a linear model for the dynamics of a quadrotor.

Linear robust control is a popular approach in the control theory that can help to mitigate the effect of unmodeled dynamics in linear-based controllers. This approach has been used in many practical problems such as surgery robots[37], parallel robots [38], harmonic drive systems [39], and vibration control of elastic gantry crane and 3D bar structure [40, 41]; and the results confirm its remarkable ability to deal with uncertainty. In [42, 43], the linear robust $\text{H}_{\infty}$ controller is employed for attitude and altitude control of a helicopter. In order to utilize the $\text{H}_{\infty}$ controller, a comprehensive identification phase is required, which provides information about the linear model of the system and the associated uncertainty.

In this paper, there is no information available about the model and parameters (e.g. CAD model) of the quadrotor. Thus, in order to overcome the aforementioned limitations of linear and nonlinear controllers, embracing lack of system model and different sources of uncertainty, a linear robust controller has been derived and implemented on a real quadrotor. In this case, a continuous time system identification from the sampled data of the real robot is required to gather information about both the system model and uncertainty. The main contributions and phases of this paper are as follows:

• •

The first phase is performing a system identification based on the experimental frequency response estimates to obtain a nominal linear model aligned with a multiplicative uncertainty block. This phase is the principle and most challenging part of control design because of system identification difficulties and its importance in designing an effective controller. In this phase, the first step is designing appropriate experiments which provide us with informative and usable data for system identification. Then, by applying existing tools and methods such as MATLAB identification and CONTSID toolbox, the nominal models and uncertainty weighting functions are acquired.

• •

Having obtained information about the system in the previous phase, a linear $H_{\infty}$ control is synthesized. In this phase, it is important to choose proper sensitivity and input weighting functions to achieve desired tracking and regulating performance.

• •

Finally, the controller is implemented on the robot and its high performance and robustness are shown by simulation and experimental results. In addition, a PID and a $\mu$ analysis controllers are designed to compare their performance with $\text{H}_{\infty}$ controller.

2 System and uncertainty encapsulation

In many systems, uncertainty emerges due to the dynamical perturbation, such as unmodelled and high frequency dynamics, in various parts of a system. To capture uncertainty in the system, we utilize the classical multiplicative perturbation model. This model helps to encapsulate the various sources of uncertainty, e.g. parametric and unmodelled dynamics, in a full block of multiplicative uncertainty. Considering the nominal system transfer function as $G_{0}$, the family of uncertain systems may be shown as follows:

[TABLE]

where $W$ is the uncertainty weighting function, and $\Delta$ is a stable perturbation satisfying $\lVert\Delta\rVert_{\infty}<1$.

Hence, for designing a robust control for a nonlinear system, an approximate linear model along with a weighting function providing information about the uncertainty profile, are required. In order to obtain a controller having proper performance in practice, it is necessary to perform a system identification on the experimental data generated from the real system. There are some challenges, however, for a quadrotor system identification:

• •

The quadrotor robots are inherent unstable, thus a closed loop identification method is recommended.

• •

For the robust control design a continuous-time model is needed.

• •

The quadrotor system is a multi-input multi-output (MIMO) system, but due to its open-loop instability and piloting difficulty, the experimental data is better to be obtained for each channel separately, then perform the system identification.

In what follows, we dwell on the closed loop system identification. For this, first we derive a theoretical model for the system, then identify the system based on the gathered experimental data.

2.1 Theoretical model

Modeling the dynamic of quadrotors has been studied extensively in the literature. Many models attempt to consider parameters and factors that are difficult to calculate or their values may vary in different situations such as the presence of wind. The goal of this paper is to design a controller that is largely independent of modeling parameters. To this end, first we model the dynamic of the robot. As shown in Fig. 1 for a drone with four rotors, which rotate in opposite directions (two rotors rotate in a clockwise direction and two other rotate counterclockwise), the angular accelerations pertain to the pitch, roll, and yaw angles are modeled by the following nonlinear models:

[TABLE]

where $\phi$, $\theta$, and $\psi$ stand for the roll, pitch and yaw angles, respectively. The remaining parameters are listed in table 1. $\left(I_{yy}-I_{zz}\right)\dot{\psi}\dot{\theta}$, $\left(I_{zz}-I_{xx}\right)\dot{\psi}\dot{\phi}$ and $\left(I_{xx}-I_{yy}\right)\dot{\theta}\dot{\phi}$ are body gyro effects. Considering $\Omega_{r}=\left(\Omega_{1}+\Omega_{3}-\Omega_{2}-\Omega_{4}\right)$, $J_{r}\dot{\theta}\Omega_{r}$ and $J_{r}\dot{\phi}\Omega_{r}$ represent propeller gyro effects.

In order to analyze the quadrotor model in the presence of a PID controller, the gyroscopic effects can be ignored compared to the motors’ action. Thus, Eq. (2) is rewritten as:

[TABLE]

The dynamics of the rotor is considered as $\frac{T_{1}}{s+T_{2}}$, in which $T_{1},T_{2}\in\mathbb{R}^{+}$ and $T_{2}$ is the pole of the system and $\nicefrac{{T_{1}}}{{T_{2}}}$ is the DC gain. Thus, the model of the robot is obtained as follows:

[TABLE]

where $u_{1}$, $u_{2}$, $u_{3}$, and $u_{4}$ are motors inputs. Assuming the motors have fast response, the Eq. (4) can be rewritten as follows, which is a double-integrator unstable system:

[TABLE]

where $A_{1},A_{2}$, and $A_{3}\in\mathbb{R}^{+}$. The system model in the presence of a PID controller, as shown in Fig. 2, is obtained as:

[TABLE]

where $K_{p_{r,p}}$, $K_{d_{r,p}}$, and $K_{i_{r,p}}$ are the nonnegative proportional, derivative, and integral gains, respectively. These equations give a useful insight into the number of poles and zeros of the system in the identification phase.

2.2 Control problem

Attitude control is the integral part of the stabilization problem and helps to reach the desired 3d orientation. In this paper, we focus on the roll and pitch angles because they play the main role in the stability of the quadrotor. Due to the coupling of the system’s states, the performance of roll and pitch has considerable effect on other states. For example, when the quadrotor moves forward, the pitch angle changes and then returns to zero. Thus, a fast and reliable control on roll and pitch angles is crucial. In addition, the disturbance exerts significant influence on the performance and stability of roll and pitch angles and consequently the overall performance of the robot.

By considering $X=\left(\phi,\dot{\phi},\theta,\dot{\theta}\right)$, the control problem is to reach $X_{d}=\left(\phi_{d},0,\theta_{d},0\right)$ in which $\phi_{d}$ and $\theta_{d}$ are desired roll and pitch angles. In our case, the $\phi_{d}$ and $\theta_{d}$ are set to zero to stabilize the quadrotor. However, as mentioned before, overcoming model uncertainty is the main goal of this paper which should be achieved with considering two following conditions.

Having a high tracking performance and disturbance rejection

To achieve this goal, we need that the controller provides a short settling time (about $0.3s$) and a small overshoot (about $10^{-6}rad$). 2. 2.

Avoiding actuator saturation

Due to the limitations of the motor, the control effort should be in a feasible range. In addition, actuator saturation may cause instability of the system.

Briefly, the problem is finding a controller providing the stability for the quadrotor by striking a balance between robustness and performance of the system, in the presence of uncertainty and the actuator saturation.

2.3 Experimental setup

The quadrotor used in the experiments (Fig 3(a)) is based on a typical quadrotor design with some structural modifications. The power management circuit board and the carbon fiber tubes are used as the frame hub and motors arm, respectively. The main goal in the design of this quadrotor is achieving a compact integration of the mechanical structure and electronic devices in order to reduce the weight while maintaining a symmetric mass distribution with a center of gravity close to the center of the cross. The vehicle propeller-tip to propeller-tip distance is $38\,cm$. This robot, embracing all the modules, sensors, battery, motors, and props has a mass of $570\,g$, and provides a total thrust of 850 g. The flight time is approximately 9 minutes with a 3-cell $1300\,mAh$ Li-Po battery.

Three main units which play significant roles in the control of the robot are explained below, and the block diagram of their structure and interconnections are illustrated in Fig. 3(b).

Embedded electronics: The robot configuration consists of a flight control, Inertial Measurement Unit (IMU), radio control unit, telemetry module, and external IO. The flight control board, which is the main board of the robot, runs on an ARM STM32F407 micro-controller, which can operate up to $168\,MHz$. Other units interface with this controller via I2C and UART protocols. Each brushless motor has its own ESC (Electronic Speed Control) units, which is connected to the main board. PWM (Pulse-Width Modulation) is used to communicate with ESC at the speed of $400\,Hz$, and the main board sends control signals to ESC to adjust the rotational speed of the motor.

Inertial Measurement Unit (IMU): The IMU utilized in this quadrotor is MPU-6050, which contains a 3-axis gyroscope, 3-axis accelerometer, 3-axis magnetometer. The main board communicates with the MPU-6050 via I2C bus and runs a sensor fusion algorithm with the frequency of $500\,Hz$. The quadrotor attitude is computed by using Mahony filter algorithm [44], which estimates the absolute roll and pitch angles via fusing the acceleration vector and angular speeds measured by the gyroscopes. The angular speed which is measured by the gyroscope associated with yaw axis is used for estimating the yaw angle.

Wireless communication systems: We established two types of system to support real time communications between the flying quadrotor and our personal computer as a ground control system (GCS): a digital radio telemetry unit ($915\,MHz$) and an analog radio link ($2.4\,GHz$). We use NRF modules as a digital radio telemetry for data acquisition. The telemetry wireless interfaces with the main board using a UART serial protocol. Necessary flight data can be either saved on the internal memory of the micro-controller or sent to the GCS at rates of up to $115200\,bps$ using telemetry communication. Futaba remote controller is used as the analog radio link. The radio control (RC) receiver is connected to the main board through the SBUS protocol that allows to receive flight commands at the speed of $100\,Kbps$.

2.4 System identification

The problem of system identification of an aircraft or rotorcraft is a topic of interest and have been widely studied. In most of the studies, the system identification methods are used to obtain parameters value in the model of the system. However, in this paper, a black-box system identification, merely based on input and output experimental data, is performed.

The first and basic step in the identification phase is selecting a proper input for applying to the system. In the frequency-domain identification, which is desired in this paper, a periodic excitation input is much preferred because it helps to cover and sweep the desired frequency spectrum effectively. The excitation inputs are better to be applied automatically (e.g programmed on the microcontroller) rather than manually (e.g radio controller). Thus, in order to obtain a continuous time model, a frequency-domain identification method is performed in which periodic chirp functions with different frequencies and amplitudes are implemented on the on-board micro-controller, and are used as exciting inputs. During the experiments, the identification for the pitch and roll angles are executed separately, and the quadrotor is supposed to track the reference chirp input. For example, if the chirp signal is selected as a reference input for the roll angle, the quadrotor uses the implemented simple PID controller for the roll angle to track the reference. Setting proper frequencies and amplitudes for chirp functions is important, because an improper frequency or amplitude, which are beyond the quadrotor ability to track, may cause inappropriate tracking and produces misleading information for identification. In this identification, the frequency and amplitude range vary in $\left[0.05Hz,5Hz\right]$ and $\left[5^{\circ},11^{\circ}\right]$, respectively, which are feasible and reasonable pitch and roll angles.

Prior to the identification phase, the probable existence of delay is studied using correlation analysis method, and according to its results, the delay of the system is negligible. Secondly, the effect of exciting the roll angle on the pitch angle, and vice versa, is analyzed. In Fig. 4(a), the pitch angle, the desired reference pitch angle, and the roll angle are plotted, for one of the experiences. As shown in this figure, the coupling is negligible and the roll angle is almost independent of the pitch angle, thus their coupling can be considered as an unmodelled system uncertainty. The same condition holds for the pitch angle as well (Fig. 4(b)). As a result, off-diagonal elements of the system model can be considered zero, which implies that the system can be represented by two SISO subsystems. Although this assumption may make the controller task more difficult at high frequencies, the controller design becomes much easier.

After performing several experiments for both pitch and roll angles, a linear transfer function should be fitted for each of them to obtain the stable closed-loop system model. We use MATLAB identification toolbox for obtaining continuous-time transfer function using time-domain data. In this step, it is important to select a proper number of poles and zeros. For this purpose, in each identification, we use different combinations of zeros and poles. However, to reduce the number of these combinations, the number of poles and zeros can be assumed close to the approximate analytical transfer function of the system, obtained in Eq. (6). Finally, for the model identification, the number of poles and zeros are considered in the ranges of $\left[2,5\right]$ and $\left[0,3\right]$, respectively. Consequently, for each experiment, fifteen proper transfer functions are obtained, and one of them is selected as the best-fitted model. In order to select a transfer function, fitness percentage of different fitted models has a decisive role. However, there are two other factors which must be considered. First, it is preferred that the selected transfer functions have similar frequency responses, meaning that the zero/pole numbers should be close to each other. Second, the transfer functions with right hand-side zeros are not selected because the quadrotor system does not show the characteristics of a system with right half-plane zeros.

Based on the obtained transfer functions, for the pitch angle, the real part of poles mainly lie in the range of $\left[-6,0\right]$, and for the roll angle, real part of poles are in the ranges of $\left[-7,0\right]$ and $\left[-10,-8\right]$. In both cases, the zeros lie in the range of $\left[-1,1\right]$. These zeros, which are located near the imaginary axis, mostly appear in pairs with poles very close to them, and therefore, they may be ignored due to the zero-pole cancellation.

By analyzing the obtained transfer functions according to the three aforementioned criteria, it is observed that in the most cases, transfer functions with three poles and no zero can suitably describe the quadrotor system behavior. Figure 5(a) and 5(b) show the histogram of the real part of poles in the transfer functions of the pitch and roll angles, with three poles and no zero, respectively. According to these figures, it is expected that poles lie in the range of $\left[-7,-3\right]$ and $\left[-10,-8\right]$. In Fig. 6, the Bode diagram of the transfer functions obtained for the pitch and roll angle are depicted. It is clear from Fig. 6 that the behavior of the transfer functions is similar, especially in frequencies less than $15Hz$ where the robot mostly operates in. Now, by having the system transfer functions pertain to both axes of pitch and roll, we select one of them as the nominal transfer function. Hence, the transfer function which lies near to the median of all responses, is considered as the nominal model. The selected nominal model transfer function for the pitch and roll angle are as follows:

[TABLE]

In order to verify the nominal models, the real experimental data (dashed) and the time-domain response (solid) are depicted in Fig. 7(a) and 7(b). For the sake of brevity, only two of the experiment results are shown. In these figures, the applied inputs to the nominal model are the same as the real experiments.

2.5 Uncertainty weighting function

As mentioned before, based on Eq. (1), $W$ is a fixed transfer function for normalizing $\Delta$. In other words, $W\Delta$ expresses the normalized system variation away from unity at each frequency, and based on Eq. (1), may be obtained by $G\left(jw\right)/G_{0}\left(jw\right)-1=\Delta\left(jw\right)W\left(jw\right)$. Since $\lVert\Delta\rVert_{\infty}<1$, the following can be concluded, which provide a practical approach for obtaining the uncertainty weighting function:

[TABLE]

After selecting a nominal model between all best-fitted transfer functions, other systems can be considered as a perturbed model of the nominal system, and $\lvert G\left(jw\right)/G_{0}\left(jw\right)-1\rvert$ is plotted for each experiment (as illustrated in Fig. 8). As a result, by estimating a transfer function, which is an upper bound for those variations, the uncertainty weighting function (the dashed red line in Fig. 8(a) and 8(b)) is acquired. The uncertainty weighting functions for pitch and roll angles are obtained as follows:

[TABLE]

The above-mentioned steps for obtaining the nominal models of the systems and the uncertainty weighting functions are summarized as below:

1. Step 1.

Designing appropriate experiments and obtaining 8 and 9 data packages for pitch and roll angles, respectively. 2. Step 2.

Fitting 15 transfer functions for each data package. 3. Step 3.

Selecting the best-fitted transfer function for each data package. 4. Step 4.

Selecting the nominal models among 8 and 9 transfer functions obtained for pitch and roll angles. 5. Step 5.

Considering the other transfer functions as perturbed models, obtaining the $\lvert G\left(jw\right)/G_{0}\left(jw\right)-1\rvert$s, and plotting them for roll and pitch angles. 6. Step 6.

Obtaining the uncertainty weighting functions, which are an upper bound for the variations obtained in Step 5, for each roll and pitch angle.

3 Robust ${H}_{\infty}$ Control ̵ُSynthesis

The goal of ${H}_{\infty}$ controller is providing a robust stable closed-loop system with high performance tracking and disturbance rejection, in the presence of uncertainty and the actuator saturation. Based on Fig. 9, in order to achieve this goal the following objectives shall be met simultaneously:

$\rVert T_{z_{1}y_{d}}\lVert_{\infty}=\rVert W_{s}S\lVert_{\infty}\leq 1$, where $T_{z_{1}y_{d}}=\frac{W_{s}}{1+CP}=W_{s}S$ is the transfer function from $y_{d}$ to $z_{1}$, which is the nominal tracking performance in an asymptotic sense. $W_{s}$ is a frequency dependent sensitivity weighting function for normalizing and shaping the closed loop performance specification. 2. 2.

$\rVert T_{z_{2}y_{d}}\lVert_{\infty}=\rVert W_{u}U\lVert_{\infty}\leq 1$, in which $T_{z_{2}y_{d}}=\frac{W_{u}C}{1+CP}=W_{u}CS=W_{u}U$ is the transfer function from $y_{d}$ to $z_{2}$, and expresses the nominal performance of the control effort. $W_{u}$ is a frequency dependent weighting function for normalizing and shaping the control input. 3. 3.

$\rVert T_{z_{3}y_{d}}\lVert_{\infty}=\rVert WT\lVert_{\infty}\leq 1$, which is the result of the small gain theorem. $T_{z_{3}y_{d}}=\frac{WCP}{1+CP}=WT$ is the transfer function from $y_{d}$ to $z_{3}$, and $W$ represents the uncertainty weighting function.

Considering this system with input, $y_{d}$, and output vector $\mathbf{z}=[z_{1},z_{2},z_{3}]^{T}$, the aforementioned conditions can be merged in an induced norm of the transfer function from $y_{d}$ to $\mathbf{z}$, in which the goal is to find a controller to minimize this norm. This problem is called a mixed-sensitivity problem and is formulated as follows:

[TABLE]

In order to define the sensitivity weighting function, first an ideal closed-loop transfer function is designed. In this case, based on the standard second order systems, $\omega_{n}$ and $\zeta$ are chosen such that the following overshoot and the settling time are satisfied

[TABLE]

To achieve $W_{s}$, we use the inverse of the desired sensitivity function:

[TABLE]

where $a\leq 1$ is a tuning parameter. Eq. (17) is modified by inserting a non-dominant pole and by slightly transferring its purely imaginary pole to the left half plane in order to convert $W_{s}$ to a strictly proper and stable weighting function, respectively. Hence, $W_{s}$ is corrected as:

[TABLE]

In the obtained nominal system, Eq. (7), the input of the system is the desired angle, which should be less than $20^{\circ}$. Therefore, $W_{u}$ is set to $0.05$. Having set the weighting functions, the mixed sensitivity problem can be solved using the robust toolbox of MATLAB. After a few correction steps, the controller for the pitch angle by tuning $a=0.88$ with $\gamma_{opt}=0.9929$ is synthesized as follows:

[TABLE]

The singular values of the closed-loop system, and the frequency response of each I/O pair are depicted in Fig. 10. The solid blue line shows the maximum singular value of the closed loop system, which is flat within a large frequency range, and is less than one in all range of the frequency domain. The dotted orange line shows the Bode diagram of $WT$ in which the maximum value is 0.42 indicating the robustness of the closed loop system with a margin of greater than two. The performance transfer function and the control effort transfer function are shown by the red and green lines, respectively. Based on Fig. 10, the magnitude of $WT$ starts to reduce significantly in frequencies higher than about 30 Hz due to the limited bandwidth of the system. About this frequency, a slight decrease also appears in the magnitude of $W_{s}S$ until roughly 900 Hz where a dramatic drop occurs. Conversely, the magnitude of $W_{u}U$ increases in these frequencies to compensate for the lack of the system stability and performance. However, since the robot working frequency is less than approximately 20 Hz, high frequencies performance of the system can be ignored. In low frequencies, as expected, the magnitude of $W_{s}S$ increases when that of $WT$ decreases.

As mentioned in section 2, a PD controller has been employed to identify the system. Therefore, to implement the final controller a cascade architecture must be used in which the PD controller is in the inner loop, and the $H_{\infty}$ controller is functioning in the outer loop as illustrated in Fig. 11. The closed-loop step response of the nominal (red) and some uncertain samples of the system are shown in Fig. 12(a). It can be seen that the controller can robustly stabilize the system and provide a fast response with a reasonable control effort, which is depicted in Fig. 12(b).

The same procedure is applied to obtain the $H_{\infty}$ controller for the roll angle. $W_{u}$ is set to 0.05, and $W_{s}$ is same as that of the pitch angle, in Eq. (18). After solving the mixed-sensitivity problem, $\gamma_{opt}$ value is obtained 0.9925, for $a=0.92$ and the following controller:

[TABLE]

The singular values of the closed-loop system and the frequency response of each I/O pair are illustrated in Fig. 13. The analysis of this plot is the same as that of the pitch angle, which is omitted for the sake of brevity. The step responses of the nominal and perturbed systems in the presence of $H_{\infty}$ controller are depicted in Fig. 14(a). The control efforts are also shown in Fig. 14(b).

4 Experimental results

4.1 Real Implementation

In this section, we verify the performance and the applicability of the controllers which were designed in the preceding section by implementing on the real robot. The control structure is the same as Fig. 11 for both pitch and roll angles. By setting the time-step of the system to 0.001, the continuous controllers are discretized with the frequency of $1\,kHz$ for the practical implementation. To make this task easier, an order reduction is performed that decreases the order of the controllers from 8 to 6. Finally, the controllers for both axes after applying the order reduction are given by:

[TABLE]

Based on the frequency response of each I/O pair, this reduction has a negligible effect.

Having implemented the controllers, some tests are performed in the outdoor and indoor environments to evaluate the $H_{\infty}$ controller ability in coping with uncertainties. In the first experiment, the operator drives the quadrotor using a radio-controller to assess the qualitative performance of the robot in comparison with a well-tuned PID controller. The tuning parameters of the PID controller for both roll and pitch angles are selected as follows:

[TABLE]

According to this experience, the maneuverability has been enhanced, and the quadrotor can quickly regulate the roll and the pitch angle in the presence of the external disturbances such as the unexpected wind effects (i.e. caused by the drone itself while flying near the walls or the weather) or some manual disturbances applied by the operator. Since the states of the quadrotor are coupled, a refinement on the altitude control of the robot is also observed. The video of this experiment is also available online111https://youtu.be/vAi4x2_XSQQ.

To evaluate the robustness, stability and tracking performance of the proposed controllers, an experiment is done in which the quadrotor should track the reference input in the outdoor environment, in the presence of light wind. Fig. 15(a) shows the tracking performance of the pitch controller and the ability of the roll controller in regulating the roll angle close to zero. In addition, as it is shown by data-tip in Fig. 15(a), the performance of the system (e.g. settling time and overshoot) is close to that of the desired transfer function obtained in Eq. (16). The performance of the well-tuned PID controller is also depicted in Fig 15(b). According to Fig. 15, the $H_{\infty}$ controller has significantly enhanced the performance of the system.

4.2 Simulation Results

After validating the practicability and efficiency of the proposed controllers through real experiments, in what follows, the performance of the controllers are evaluated in the simulation and is compared with a well-tuned PID and another robust controller obtained by $\mu$-synthesis method. The robust stability and robust performance of the controllers are also analyzed using the structural singular values. The gains of PID controller are as below:

[TABLE]

Fig. 16(a) illustrates the step-response of an uncertain system in the presence of the disturbances and uniform random sensor noise. The amplitude of the disturbances is $0.1$, and the maximum value of the random sensor noise is $0.02$. Fig. 16(b), depicts the roll angle performance in a similar experiment. Based on these figures, the $H_{\infty}$ controller’s overshoots for the pitch and roll systems are about $1\%$ and $0.1\%$, respectively, which are less than that of PID and $\mu$ synthesis controllers which are approximately $10\%$. In addition, the $H_{\infty}$ controller regulates the system faster, and its settling time is roughly $0.35\,sec$. The settling times of both PID and $\mu$ synthesis controllers are about $1\,sec$. $H_{\infty}$ controller also has better performance when disturbances occur.

In Fig. 17, the robust stability of the $H_{\infty}$, $\mu$ synthesis, and PID controllers are evaluated by structural singular values. As expected, the $H_{\infty}$ and $\mu$ synthesis controllers fulfill the robust stability condition and their $\mu$ values are less than $1$. The robust stability condition is also met for the PID controller, however, in frequencies about $20\,rad/sec$ gets close to $1$. Moreover, when the PID parameters are set in such a way that the system has a faster response, the robust stability will be lost. The robust performance analysis of the roll and pitch systems using structural singular values are depicted in Fig. 18. Predictably, the $\mu$ synthesis controller offers the robust performance in contrast to $H_{\infty}$ and PID controllers. The $\mu$ synthesis controller is a conservative approach to have the robust performance, which results in poor time-domain performance in comparison with the $H_{\infty}$ controller. Moreover, $\mu$ synthesis method produces a high-order controller which makes its implementation on a real robot challenging [45].

5 Conclusion

This research has studied the robust attitude control of a quadrotor with taking into account the uncertainty and inputs constraints. Knowing that the model and parameters of the system do not exist, a continuous-time black-box model identification based on the real sampled data is adopted to acquire both nominated linear model of the system and the uncertainty model encapsulating the deviation of the nonlinear system from the nominal model. Now, the essential prerequisites for applying the robust $H_{\infty}$ controller is provided. By solving the mixed-sensitivity problem, a robust stabilizing $H_{\infty}$ controller is obtained which satisfies tracking, disturbance attenuation and input saturation objectives. Having accomplished that, the pitch and roll controllers are implemented on the robot, and qualitative and quantitative experiments are performed to verify the performance of $H_{\infty}$ controllers. According to both experimental and simulation results, the controllers successfully fulfill the aforementioned objectives. Finally, by comparing the performance of the proposed controllers with PID and $\mu$ synthesis controllers, it is demonstrated that the $H_{\infty}$ controller has a better performance.