Development of a robust attitude control for nonidentical rotor quadrotors using sliding mode control

Abstract

This work proposes a novel development of sliding mode control (SMC) for quadrotor helicopters utilizing offset cancellation technique. Most of the works on sliding mode control for quadrotors assumed that all rotors were identical, but in fact, the rotors even from the same factory are not exactly the same in aerodynamics properties. In addition, rotor deformation due to shock or heavy usage during operation may cause changes in parameters such as lift and drag coefficients, and so on, which, lead to offset error in the responses of attitude control. By applying the proposed offset cancellation technique to the sliding mode control, the attitude error is forced to have a zero average. The proposed controller is used to control attitude and position of a quadrotor helicopter which is composed of four direct current (DC) brushless motors to generate lift forces independently. Performance comparison with the other control algorithms such as the proportional–integral–derivative controller and the conventional sliding mode control is simulated and experimented on two real quadrotor platforms and results are discussed.

Introduction

Nowadays, quadrotor helicopters are easily available and gain attention from many users and researchers. With its simple platform, a quadrotor helicopter can provide more maneuvering ability than the conventional single or double rotor helicopter. Instead of controlling the blade pitch angle, which requires complicated mechanism, speeds of the four rotors are controlled independently to achieve the desired lift force and moment. The researches on the quadrotor helicopter can be decomposed into four main activities: construction of the helicopter, analysis of its aerodynamics, identification of the dynamic model, and design of the control algorithm. The basic structure of quadrotor helicopter consists of cross configuration of four DC brushless motors. Two pairs of propellers are installed and rotated in opposite directions in order to cancel vertical angular momentum of each pair of propellers during hovering.

Even though the quadrotor helicopter has six degrees of freedom in its motion, only four rotors are sufficient to control the quadrotor helicopter to reach the desired operation point. ¹ Aerodynamics of the propellers and the frame are very important. Various aerodynamics effects of the propellers mounting on the frame always interfere and produce a coupled dynamic model. ²

Regarding the control aspect, classical proportional–integral–derivative controller (PID) control is one of the most popular algorithms applied in the quadrotor helicopters due to its simplicity. ^{1,3 –6} However, the PID controller may not provide the best performance due to sensor noises, disturbance, and parameter variations. Other algorithms, such as accelerometer feedback, ⁷ disturbance compensation, ⁴ sliding mode control (SMC), ⁸ and time optimum control, ⁹ were applied to control the quadrotor helicopters. SMC is widely used in several applications, such as navigation of a mobile robot ^10,11 and stability control of an underwater vehicle. ¹² To achieve good performance under uncertain conditions, SMC has been used to control the quadrotor helicopter. However, chattering of control signal in the SMC degrades control performance of the conventional SMC. It produces fluctuation of control signal. To reduce this effect, a chattering-free technique was proposed ¹³ to obtain better performance. Synthesis control, ¹⁴ active disturbance rejection technique, ¹⁵ and terminal SMC ¹⁶ were applied and able to demonstrate fast tracking with small steady state error for unmanned aerial vehicles.

Most researches which used SMC to control quadrotor helicopters assumed that all rotors were identical. This article discusses the effect of nonidentical rotors and proposes a novel SMC using an offset cancellation technique to control attitude and position of a quadrotor helicopter when the properties of rotor are not identical.

In the following sections, dynamics model of quadrotor will be introduced, followed by derivation of the control algorithm (conventional and the proposed SMC). System requirements will then be shown. Simulation and experimental results are finally listed.

Quadrotor helicopter dynamics model

All forces acting on a helicopter frame are considered in the body reference frame ( $\widehat{x}$ , $\widehat{y}$ , and $\widehat{z}$ ), which moves and rotates with the helicopter motion. The body reference frame is transformed to the fixed earth reference frame by the Euler’s transformation matrix. Figure 1 shows the body reference frame of the model. The unit vectors $\widehat{x}$ , $\widehat{y}$ , and $\widehat{z}$ are taken on the body reference frame, and $\widehat{X}$ , $\widehat{Y}$ , and $\widehat{Z}$ represent the unit vectors taken on the fixed earth reference frame. Roll (ϕ), pitch (θ), and yaw (φ) rotation axes lie on the $\widehat{x}$ , $\widehat{y}$ , and $\widehat{z}$ direction, respectively.

Figure 1.

Quadrotor helicopter and its body reference coordination.

The forces acting on a quadrotor helicopter consist of lifting forces from the rotor F_L , H-forces of the rotor F_H , drag force on the frame body F_D , gravitation force F_g , and external force F_E . The total force F can be expressed by

$F=F_L+F_H+F_E+F_D+F_g$

1

The rotor lift force F_L points to $-\widehat{z}$ direction in the body reference frame as shown in equation (2). The lift force from rotor i is determined by F_Li

${\mathbf{\text{F}}}_{\mathbf{\text{L}}}={\mathbf{\text{F}}}_{\mathbf{\text{L1}}}+{\mathbf{\text{F}}}_{\mathbf{\text{L2}}}+{\mathbf{\text{F}}}_{\mathbf{\text{L3}}}+{\mathbf{\text{F}}}_{\mathbf{\text{L4}}}$

2

The H-force F_H or rotor drag is a function of the rotor speed and linear velocity motion of the helicopter. The F_H is neglected in this work due to the small linear velocity motion of the helicopter compared to the rotor speed. When the helicopter moves with both the linear and angular motions, the external force due to centripetal force denoted by F_E exerts on the helicopter body. A relation between centripetal force and velocity of a helicopter can be expressed as follows

$F_E=m\left(\Omega \times v\right)$

3

The vector v in equation (3) is the linear velocity, while Ω is the angular velocity of the helicopter around the body axis. To simplify the dynamic model, the external force F_E is neglected since when a helicopter moves in linear motion, the helicopter angular velocity Ω is controlled to be zero. The frame dragging force F_D is small and negligible. The gravitation force F_g points to $\widehat{Z}$ direction of the fixed earth reference frame. A relation between body reference frame and fixed earth reference frame can be written following the Euler’s angle definition

$\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]=R\left(\varphi ,\theta ,\phi \right)\left[\begin{array}{c}\widehat{x}\\ \widehat{y}\\ \widehat{z}\end{array}\right]$

4

Let c and s denote cosine and sine functions, respectively. The rotation matrix can be expressed as

$R\left(\varphi ,\theta ,\phi \right)=\left[\begin{array}{ccc}\text{c}\phi \text{c}\theta & -\text{s}\phi \text{c}\varphi +\text{c}\phi \text{s}\theta \text{s}\varphi & \text{s}\phi \text{s}\varphi +\text{c}\phi \text{c}\varphi \text{s}\theta \\ \text{s}\phi \text{c}\theta & \text{c}\phi \text{c}\varphi +\text{s}\phi \text{s}\theta \text{s}\varphi & -\text{c}\phi \text{s}\varphi +\text{s}\theta \text{s}\phi \text{c}\varphi \\ -\text{s}\theta & \text{c}\theta \text{s}\varphi & \text{c}\theta \text{c}\varphi \end{array}\right]$

Because the lift forces from all rotors are along the $\widehat{z}$ direction of the body reference frame, they can be transformed into the fixed earth reference frame using

$F=\left[\begin{array}{c}{F}_X\\ F_Y\\ F_Z\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ mg\end{array}\right]-F_L\left[\begin{array}{c}\text{s}\phi \text{s}\varphi +\text{c}\phi \text{c}\varphi \text{s}\theta \\ -\text{c}\phi \text{s}\varphi +\text{s}\theta \text{s}\phi \text{c}\varphi \\ \text{c}\theta \text{c}\varphi \end{array}\right]$

5

The gravitational force, mg, is along the $\widehat{Z}$ direction of the fixed earth reference frame. The rotor lift force, F_L , is a function of rotor speed ω. Define ω_i as the rotor speed of motor i, where i = 1, 2, 3, 4, and rotor lift coefficients are defined as k_i

$\left|F_L\right|=k_1{\omega }_1^2+k_2{\omega }_2^2+k_3{\omega }_3^2+k_4{\omega }_4^2$

6

where k_i is always positive, equation (5) is rewritten as

$F=\left[\begin{array}{c}{F}_X\\ F_Y\\ F_Z\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ mg\end{array}\right]-\left[\begin{array}{c}\text{s}\phi \text{s}\varphi +\text{c}\phi \text{c}\varphi \text{s}\theta \\ -\text{c}\phi \text{s}\varphi +\text{s}\theta \text{s}\phi \text{c}\varphi \\ \text{c}\theta \text{c}\varphi \end{array}\right]k_1{\omega }_1^2+k_2{\omega }_2^2+k_3{\omega }_3^2+k_4{\omega }_4^2$

7

Total torque acting on the helicopter is considered in the body reference frame and can be written as

$\tau ={\tau }_{\mathit{L}}+{\tau }_{\mathit{r}}+{\tau }_{\mathit{g}}+{\tau }_{\mathit{d}}+{\tau }_{\mathit{m}}+{\tau }_{\mathit{f}}$

8

τ_L from equation (8) is the torque due to lift force from each motor

${\tau }_{\mathit{L}}=\left[\begin{array}{c}{\tau }_{L\varphi }\\ {\tau }_{L\theta }\\ {\tau }_{L\phi }\end{array}\right]=\left[\begin{array}{c}l\left(k_4{\omega }_4^2-k_2{\omega }_2^2\right)\\ l\left(k_1{\omega }_1^2-k_3{\omega }_3^2\right)\\ 0\end{array}\right]$

9

where l is the distance from each rotor to the center of gravity of a helicopter, which is considered to be symmetrical. The torque τ_r from equation (8) is the gyroscopic torque due to rotor rotation

${\tau }_{\mathit{r}}=\Omega \times J\omega$

10

where the variable J is the moment of inertia of the rotor along rotor rotation axis and is the angular velocity of the frame along the roll (ϕ), pitch (θ), and yaw (φ) axes

$\Omega =\left[\begin{array}{c}\dot{\varphi }\\ \dot{\theta }\\ \dot{\phi }\end{array}\right]$

11

From equation (11), the gyroscopic torque due to rotor rotation can be expressed as

${\tau }_{\mathit{r}}=\left[\begin{array}{c}{\tau }_{r\varphi }\\ {\tau }_{r\theta }\\ {\tau }_{r\phi }\end{array}\right]=J\left[\begin{array}{c}\dot{\theta }\\ -\dot{\varphi }\\ 0\end{array}\right]\left(e_1{\omega }_1+e_2{\omega }_2+e_3{\omega }_3+e_4{\omega }_4\right)$

12

where e_i is the direction of rotor rotation corresponding to the ith rotor along $\widehat{z}$ axis, e_i = 1 for clockwise, and e_i = −1 for counterclockwise direction.

The variable τ_g in equation (8) is the gyroscopic torque due to rotation of frame

${\tau }_{\mathit{g}}=\Omega \times I\Omega$

13

The variable I is the moment of inertia of the frame along $\widehat{x}$ , $\widehat{y}$ , and $\widehat{z}$ body axes, and equation (13) can be rewritten as

${\tau }_{\mathit{g}}=\left[\begin{array}{c}{\tau }_{g\varphi }\\ {\tau }_{g\theta }\\ {\tau }_{g\phi }\end{array}\right]=\left[\begin{array}{c}\left(I_z-I_y\right)\dot{\theta }\dot{\phi }\\ \left(I_x-I_z\right)\dot{\varphi }\dot{\phi }\\ \left(I_y-I_x\right)\dot{\varphi }\dot{\theta }\end{array}\right]$

14

The variable τ_d from equation (8) is the torque due to rotor drag force about $\widehat{z}$ axis. It is produced from the rotor drag force F_di, which is a function of drag force coefficient, β_i, and rotor speed. Define F_di as the drag force from the ith rotor which can be written as

$F_{di}={\beta }_i{\omega }_i^2$

15

The maximum torque from a rotor, ${\tau }_{di\_\text{max}}$ , is obtained when the rotor rotates to the position illustrated in Figure 2, where ${\tau }_{di\_\text{max}}$ is obtained as

${\tau }_{di\_\text{max}}=-e_i{F}_{di}\left(L+r\prime \right)+e_i{F}_{di}\left(L-r\prime \right)$

16

Figure 2.

Torque due to the drag force acting on a rotor.

We can visualize the drag force on a rotor as shown in Figure 3, where rotor drag force is a function of rotor coefficients and radius r′. ² However, since the rotor rotates from 0 to 2 π rad as shown in Figure 3, the torque due to drag force varies with the angle of rotation. The radius of the left and right rotor blades is equal; therefore, the total torque at the rotor position between 0 − π and π − 2π is identical. The drag torque at any rotor angles can be written as

${\tau }_{di\_\theta }=-e_i{F}_{di}\left(L+r\prime \text{sin}\theta \right)+e_i{F}_{di}\left(L-r\prime \text{sin}\theta \right)$

17

${\tau }_{di\_\theta }=-e_{i}2r\prime F_{di}\left|\text{sin}\theta \right|$

18

Figure 3.

The relation between rotor torque due to drag force and rotor angle.

The torque in equation (18) is the absolute value of sine wave, as shown in Figure 3. Therefore, the average torque from equation (18) can be calculated

${\tau }_{di}=\frac{1}{\pi }\underset{0}{\overset{\pi }{\int }}-e_{i}2r\prime F_{di}\text{ sin}\theta \text{d}\theta$

19

${\tau }_{di}=\frac{-4e_{i}r\prime F_{di}}{\pi }=\frac{-4e_{i}r\prime }{\pi }{\beta }_i{\omega }_i^2$

20

${\tau }_{di}=-{\alpha }_i{e}_i{\omega }_i^2$

21

The parameter α_i in equation (21) is drag torque coefficient of rotor i which is a positive constant. From equation (21), total drag torque acting on the body frame can be written as

${\tau }_{\mathit{d}}=\left[\begin{array}{c}{\tau }_{d\varphi }\\ {\tau }_{d\theta }\\ {\tau }_{d\phi }\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right]\left(-e_1{\alpha }_1{\omega }_1^2-e_2{\alpha }_2{\omega }_2^2-e_3{\alpha }_3{\omega }_3^2-e_4{\alpha }_4{\omega }_4^2\right)$

22

The torque acting on the rotor axis resulting from change of angular momentum of the rotors, τ_m , in equation (8) can be expressed as

${\tau }_{\mathit{m}}=\left[\begin{array}{c}{\tau }_{m\varphi }\\ {\tau }_{m\theta }\\ {\tau }_{m\phi }\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right]J\left(e_1{\dot{\omega }}_1+e_2{\dot{\omega }}_2+e_3{\dot{\omega }}_3+e_4{\dot{\omega }}_4\right)$

23

The variable τ_f in equation (8) is the torque due to the difference of air velocity at the advancing and receding of the rotor blades when the helicopter is moving. However, at low speed linear motion, τ_f , is negligible. The total torque τ acting on the body frame is expressed by

$\tau =\left[\begin{array}{c}l\left(k_4{\omega }_4^2-k_2{\omega }_2^2\right)\\ l\left(k_1{\omega }_1^2-k_3{\omega }_3^2\right)\\ \left(-e_1{\alpha }_1{\omega }_1^2-e_2{\alpha }_2{\omega }_2^2-e_3{\alpha }_3{\omega }_3^2-e_4{\alpha }_4{\omega }_4^2\right)\end{array}\right]+\left[\begin{array}{c}\dot{\theta }J\left(e_1{\omega }_1+e_2{\omega }_2+e_3{\omega }_3+e_4{\omega }_4\right)+\left(I_z-I_y\right)\dot{\theta }\dot{\phi }\\ -\dot{\varphi }J\left(e_1{\omega }_1+e_2{\omega }_2+e_3{\omega }_3+e_4{\omega }_4\right)+\left(I_x-I_z\right)\dot{\varphi }\dot{\phi }\\ J\left(e_1{\dot{\omega }}_1+e_2{\dot{\omega }}_2+e_3{\dot{\omega }}_3+e_4{\dot{\omega }}_4\right)+\left(I_y-I_x\right)\dot{\varphi }\dot{\theta }\end{array}\right]$

24

After substitution of all torques in equations (8) to (24), the dynamic model of the quadrotor helicopter can be expressed as

$\left[\begin{array}{c}\ddot{z}\\ \ddot{\varphi }\\ \begin{array}{c}\ddot{\theta }\\ \ddot{\phi }\end{array}\end{array}\right]=\left[\begin{array}{c}g-\frac{1}{m}\left(k_1{\omega }_1^2+k_2{\omega }_2^2+k_3{\omega }_3^2+k_4{\omega }_4^2\right)\text{cos}\varphi \text{cos}\theta \\ \frac{l}{I_x}\left(k_4{\omega }_4^2-k_2{\omega }_2^2\right)\\ \begin{array}{c}\frac{l}{I_y}\left(k_1{\omega }_1^2-k_3{\omega }_3^2\right)\\ \frac{1}{I_z}\left(-e_1{\alpha }_1{\omega }_1^2-e_2{\alpha }_2{\omega }_2^2-e_3{\alpha }_3{\omega }_3^2-e_4{\alpha }_4{\omega }_4^2\right)\end{array}\end{array}\right]+\left[\begin{array}{c}\frac{\dot{\theta }J}{I_x}\left(e_1{\omega }_1+e_2{\omega }_2+e_3{\omega }_3+e_4{\omega }_4\right)+\frac{\left(I_z-I_y\right)}{I_x}\dot{\theta }\dot{\phi }\\ -\frac{\dot{\varphi }J}{I_y}\left(e_1{\omega }_1+e_2{\omega }_2+e_3{\omega }_3+e_4{\omega }_4\right)+\frac{\left(I_x-I_z\right)}{I_y}\dot{\varphi }\dot{\phi }\\ \frac{J}{I_z}\left(e_1{\dot{\omega }}_1+e_2{\dot{\omega }}_2+e_3{\dot{\omega }}_3+e_4{\dot{\omega }}_4\right)+\frac{\left(I_y-I_x\right)}{I_z}\dot{\varphi }\dot{\theta }\end{array}\right]$

25

where I_x, I_y, and I_z are the moments of inertia of the helicopter along $\widehat{x}$ , $\widehat{y}$ , and $\widehat{z}$ axes of the body reference frame, respectively.

Control algorithm

Since moment of inertia of rotor is small compared to that of the frame, the ratio of J/I is very small. Thus, the gyroscopic torques along the roll and pitch axes and angular momentum along the yaw axis, according to equation (25), are small and negligible. When all of the rotor coefficients are identical, the dynamic model is simplified to

$\left[\begin{array}{c}\ddot{z}\\ \ddot{\varphi }\\ \begin{array}{c}\ddot{\theta }\\ \ddot{\phi }\end{array}\end{array}\right]=\left[\begin{array}{c}g-\frac{1}{m}\left(k_1{\omega }_1^2+k_2{\omega }_2^2+k_3{\omega }_3^2+k_4{\omega }_4^2\right)\text{cos}\varphi \text{cos}\theta \\ \frac{l}{I_x}\left(k_4{\omega }_4^2-k_2{\omega }_2^2\right)+\frac{\left(I_z-I_y\right)}{I_x}\dot{\theta }\dot{\phi }\\ \begin{array}{c}\frac{l}{I_y}\left(k_1{\omega }_1^2-k_3{\omega }_3^2\right)+\frac{\left(I_x-I_z\right)}{I_y}\dot{\varphi }\dot{\phi }\\ \frac{1}{I_z}\left(-e_1{\alpha }_1{\omega }_1^2-e_2{\alpha }_2{\omega }_2^2-e_3{\alpha }_3{\omega }_3^2-e_4{\alpha }_4{\omega }_4^2\right)\end{array}\end{array}\right]$

26

Sliding mode control

SMC for the roll angle when all rotors have the same lift and drag coefficients is expressed by the sliding surface and the dynamics equation

$S={\lambda }_r\left(\varphi -{\varphi }_d\right)+\dot{\varphi }$

27

$\ddot{\varphi }=\frac{l}{I_x}{u}_2+\frac{\left(I_z-I_y\right)}{I_x}\dot{\theta }\dot{\phi }$

28

Whereas ϕ_d is the designed roll angle. The derivative of the sliding surface becomes

$\dot{S}={\lambda }_r\dot{\varphi }+\ddot{\varphi }$

29

$\dot{S}={\lambda }_r\dot{\varphi }+\frac{l}{I_x}{u}_2+\frac{\left(I_z-I_y\right)}{I_x}\dot{\theta }\dot{\phi }$

30

The equivalent control input $\widehat{u}$ for the roll angle is determined as

${\widehat{u}}_2=-\frac{{\lambda }_r\dot{\varphi I_x}}{l}-\frac{\left(I_z-I_y\right)}{l}\dot{\theta }\dot{\phi }$

31

The control input u for the roll angle becomes

$u_2=-\frac{{\lambda }_r\dot{\varphi }{I}_x}{l}-\frac{\left(I_z-I_y\right)}{l}\dot{\theta }\dot{\phi }-M_2\text{sign}\left({\lambda }_r\left(\varphi -{\varphi }_d\right)+\dot{\varphi }\right)$

32

Similarly, the dynamic models of the altitude, pitch, and yaw angles can be expressed as

$\ddot{z}= g+\frac{1}{m}{u}_1\text{cos}\varphi \text{cos}\theta$

33

$\ddot{\theta }=\frac{l}{I_y} u_3+\frac{\left(I_x-I_z\right)}{I_y}\dot{\varphi }\dot{\phi }$

34

$\ddot{\psi }=\frac{u_4}{I_z}$

35

where the control laws for the altitude, pitch, and yaw angles can be written as

$u_1=\frac{-m\left(g+{\lambda }_z\dot{z}\right)}{\text{cos}\varphi \text{cos}\theta }-M_1\text{sign}\left({\lambda }_z\left(z-z_d\right)+\dot{z}\right)$

36

$u_3=-\frac{{\lambda }_p\dot{\theta }{I}_y}{l}-\frac{\left(I_x-I_z\right)}{l}\dot{\varphi }\dot{\phi }-M_3\text{sign}\left({\lambda }_p\left(\theta -{\theta }_d\right)+\dot{\theta }\right)$

37

$u_4=-{\lambda }_y\dot{\phi }{I}_z-M_4\text{sign}\left({\lambda }_y\left(\phi -{\phi }_d\right)+\dot{\phi }\right)$

38

where $\dot{\varphi }$ , $\dot{\theta }$ , and $\dot{\phi }$ are the angular velocities around the body rotation axis. The angular velocities are obtained from a three-axis rate gyro sensor. The data from the rate gyro are fused with the data from an accelerometer using Kalman filter. After that, the information of the Euler’s body angles (ϕ, θ, and φ) are obtained and used in the control laws as expressed in equations (32) and (36) to (38).

The proposed SMC

The control law derived in the previous section is based on the assumption that all rotor coefficients are identical. Preliminary experiments have shown that the state outputs (i.e. angle and altitude) were stable but they did not converge to the desired point. To investigate this phenomenon, the Lyapunov candidate function is introduced

$V=\frac{1}{2}{S}^2$

39

The derivative of equation (39) is

$\dot{V}=S\dot{S}$

40

From equations (27) and (30)

$\dot{V}=\left({\lambda }_r\left(\varphi -{\varphi }_d\right)+\dot{\varphi }\right)\left({\lambda }_r\dot{\varphi }+\frac{l}{I_x}{u}_2+\frac{\left(I_z-I_y\right)}{I_x}\dot{\theta }\dot{\phi }\right)$

41

Based on control law in equation (32), when the controlled states converge to the set points, mean of force u₂ from equation (32) becomes zero to produce zero torque, the error term $\left(\varphi -{\varphi }_d\right)$ converges to zero, and all angular velocities become zero. This yields asymptotic stable

${\dot{V}}_{\left(u_2=0\right)}={\dot{V}}_{\left(\varphi -{\varphi }_d=0\right)}=0$

42

Normally, u₂ is a direct result from the lift coefficient of the rotors, k, and the rotor speed, ω

$u_2=k\left({\omega }_4^2-{\omega }_2^2\right)$

43

In practice, if the coefficient k has been changed according to the actual usage, it can cause the nonzero average of u ₂. According to our experiment, rotor speed ω can be precisely controlled; this leads to my conclusion that the coefficient k of each rotor must be different. Therefore, equation (43) should be rewritten as

${u\prime }_2=k_4{\omega }_4^2-k_2{\omega }_2^2=u_2+u_{\xi }$

44

Equation (44) indicates that with nonidentical lift coefficients, it does not provide zero torque even though u ₂ from equation (32) is zero, and thus equation (42) is invalid in real system. Since the coefficients k in the dynamic model (25) are not the same, the dynamic model cannot be simplified as expressed in equation (26). For the roll axis, the dynamic model is expressed by

$\ddot{\varphi }=\frac{l}{I_x}\left(k_4{\omega }_4^2-k_2{\omega }_2^2\right)+\frac{\left(I_z-I_y\right)}{I_x}\dot{\theta }\dot{\phi }$

45

The derivative of the sliding surface is expressed by

$\dot{S}={\lambda }_r\dot{\varphi }+\frac{l}{I_x}\left(k_4{\omega }_4^2-k_2{\omega }_2^2\right)+\frac{\left(I_z-I_y\right)}{I_x}\dot{\theta }\dot{\phi }$

46

In the actual quadrotor helicopter, the lift coefficients k ₂ and k ₄ are not identical. Define the lift coefficients as the summation of an approximate value k and the correction value ε

$k_2=k+{\epsilon }_2,k_4=k+{\epsilon }_4$

47

The derivative of the sliding surface now becomes

$\dot{S}={\lambda }_r\dot{\varphi }+\frac{l}{I_x}\left(k{\omega }_4^2-k{\omega }_2^2\right)+\frac{l}{I_x}\left({\epsilon }_{4}k{\omega }_4^2-{∊}_{2}k{\omega }_2^2\right)+\frac{\left(I_z-I_y\right)}{I_x}\dot{\theta }\dot{\phi }$

48

Separate part of the correction value into a new function

$f_{\epsilon }=\frac{l}{I_x}\left({\epsilon }_{4}k{\omega }_4^2-{∊}_{2}k{\omega }_2^2\right)$

49

The derivative of the sliding surface becomes

$\dot{S}={\lambda }_r\dot{\varphi }+\frac{l}{I_x}\left(k{\omega }_4^2-k{\omega }_2^2\right)+f_{\epsilon }+\frac{\left(I_z-I_y\right)}{I_x}\dot{\theta }\dot{\phi }$

50

Denote the steady state position error as e_ϕ

$e_{\varphi }=\varphi -{\varphi }_d$

51

$\underset{\varphi \to {\varphi }_d}{\text{lim}}{e}_{\varphi }=0$

52

Nonidentical parameters cause offset error at steady state. When the mean of e_ϕ is close to zero, the correction function, f_ε, converges to the constant value. This constant value is the correction force that is lost due to the approximation k. Let f_ε denote the mean of correction function

$e_{\varphi }\to 0\text{ then }{f}_{\epsilon }\to \text{constant}$

53

When f_ε is constant, its derivative is zero

$\frac{\text{d}}{\text{d}t}{f}_{\epsilon }\to 0$

54

The position error is proportion to the derivative of f_ε

$e_{\varphi }\infty \frac{\text{d}}{\text{d}t}{f}_{\epsilon }$

55

$C_{\varphi }{e}_{\varphi }=\frac{\text{d}}{\text{d}t}{f}_{\epsilon }$

56

The variable C_ϕ in equation (56) represents a positive constant value. From equation (56), the mean of the parameter error is the integration of position error

$f_{\epsilon }=\int C_{\varphi }{e}_{\varphi }\text{ d}t$

57

When f_ε and e_ϕ are uniform-distributed functions, the mean function in equation (57) can be rewritten as

$\frac{1}{T}\underset{0}{\overset{T}{\int }}{f}_{\epsilon }\text{ d}t=\frac{1}{T}\underset{0}{\overset{T}{\int }}\left(\int C_{\varphi }{e}_{\varphi }\text{ d}t\right)\text{d}t$

58

$\frac{\text{d}}{\text{d}t}\left(\frac{1}{T}\underset{0}{\overset{T}{\int }}{f}_{\epsilon }\text{ d}t\right)=\frac{\text{d}}{\text{d}t}\left(\frac{1}{T}\underset{0}{\overset{T}{\int }}\left(\int C_{\varphi }{e}_{\varphi }\text{ d}t\right)\text{d}t\right)$

59

From equation (59), the relationship between parameter error and position error becomes

$f_{\epsilon }=\int C_{\varphi }{e}_{\varphi }\text{ d}t$

60

Now substitute equation (60) into equation (48)

$\dot{S}={\lambda }_r\dot{\varphi }+\frac{l}{I_x}\left(k{\omega }_4^2-k{\omega }_2^2\right)+\frac{\left(I_z-I_y\right)}{I_x}\dot{\theta }\dot{\phi }+\int C_{\varphi }{e}_{\varphi }\text{ d}t$

61

The equivalent control input for the roll angle is determined as

${\widehat{u}}_2=-\frac{{\lambda }_r\dot{\varphi }{I}_x}{l}-\frac{\left(I_z-I_y\right)}{l}\dot{\theta }\dot{\phi }-\int C_{\varphi }{e}_{\varphi }\text{ d}t$

62

Thus, the control law of the roll angle can be written as

$u_2=-\frac{{\lambda }_r\dot{\varphi }{I}_x}{l}-\frac{\left(I_z-I_y\right)}{l}\dot{\theta }\dot{\phi }-\int C_{\varphi }{e}_{\varphi }\text{ d}t-M_2\text{sign}\left({\lambda }_r\left(\varphi -{\varphi }_d\right)+\dot{\varphi }\right)$

63

With similar derivations, the control laws for the remaining attitude control can be determined

$u_1=-\frac{m\left(g+{\lambda }_z\dot{z}\right)}{\text{cos}\varphi \text{cos}\theta }-\int C_z{e}_z\text{ d}t-M_1\text{sign}\left({\lambda }_z\left(z-z_d\right)+\dot{z}\right)$

64

$u_3=-\frac{{\lambda }_p\dot{\theta }{I}_y}{l}-\frac{\left(I_x-I_z\right)}{l}\dot{\varphi }\dot{\phi }-\int C_{\theta }{e}_{\theta }\text{ d}t-M_3\text{sign}\left({\lambda }_p\left(\theta -{\theta }_d\right)+\dot{\theta }\right)$

65

$u_4=-{\lambda }_y\dot{\phi }{I}_z-\int C_{\phi }{e}_{\phi }\text{ d}t-M_4\text{sign}\left({\lambda }_y\left(\phi -{\phi }_d\right)+\dot{\phi }\right)$

66

Parameter identification

Several control algorithms are designed for quadrotors to control dynamics behavior of the quadrotors. In general, the output from control law is the required force of each motor which is determined from the rotor speed using a simple relation between rotor speed and thrust force. However, there is no force control loop to ensure that each rotor produces the force as required by the control law. This problem causes the ideal performance quadrotor differs from the real performance and requires improvement. This behavior was confirmed from the experiments, which showed that the state outputs (angle and altitude) were stable but did not converge to the desired values.

To confirm that the quadrotor rotors have different behaviors, the lift and drag coefficients k and α from all rotors are investigated. From the experimental results, the relations between rotor speed and lift force are obtained as shown in Figure 4. According to the rotor aerodynamic, lift force, F_Li , is the multiplication of square of the rotor speed, ω_i, and the lift coefficient, k_i , of the ith rotor

$F_{Li}=k_i{\omega }_i^2$

67

Figure 4.

Experimental result of lift force and rotor speed.

The lift coefficient of each rotor is identified using the experimental result in Figure 4 and equation (67). The results from statistical regression are obtained as follows

$F_{L1}=1.289\times {10}^{-5}{\omega }_1^2$

68

$F_{L2}=0.972\times {10}^{-5}{\omega }_2^2$

69

$F_{L3}=1.239\times {10}^{-5}{\omega }_3^2$

70

$F_{L4}=1.014\times {10}^{-5}{\omega }_4^2$

71

The drag coefficient is identified from the power and rotor speed relationship. The parameter α_i in the dynamic model (21) depends on the rotor characteristic and can be determined by

$P_i={\alpha }_i{\omega }_i^3$

72

where P_i and ω_i are the power and the speed of the ith rotor, respectively. From the experimental results, the relation between the motor power and the rotor speed is shown in Figure 5.

Figure 5.

Experimental result of motor power and rotor speed.

Similar to the lift coefficient, the drag torque coefficient of each rotor is determined. Using the experimental results in Figure 5, equation (72), the drag torque acting on yaw axis in each rotor can be computed as follows

${\tau }_{d\phi \_1}=0.192\times {10}^{-6}{\omega }_1^2$

73

${\tau }_{d\phi \_2}=0.40\times {10}^{-6}{\omega }_2^2$

74

${\tau }_{d\phi \_3}=0.180\times {10}^{-6}{\omega }_3^2$

75

${\tau }_{d\phi \_4}=0.149\times {10}^{-6}{\omega }_4^2$

76

From the parameter identification, four different lift and drag coefficients are obtained. However, most of the researches of quadrotor control always assumed identical lift and drag coefficients for all rotors using their average values.

System requirements and design

System requirements

In order to verify that the proposed offset cancellation technique performs efficiently regardless of the size of quadrotor, two quadrotor helicopters are used in this research. They have different sizes, weights, and physical parameters.

Brushless DC motors are used in this work to drive the rotors. The motor power must be sufficient to generate lift force under full load. The selection of motor’s power has been described in Appendix 1. In this work, the motor’s output power is 100 W for the large quadrotor and 45 W for the small quadrotor. A lithium polymer (LiPo) battery with 12-V 3 Ah capacity is used in the small quadrotor and 12-V 5 Ah capacity is used in the large quadrotor. The LiPo battery is suitable as the power source in quadrotor because it provides a high energy but light weight. ¹⁷

A 32-bit microcontroller, PIC32MX460F512L-80I/PT, is selected for the processing unit of the quadrotor. It is a powerful microcontroller available in the market and operates at 80 MHz which is sufficiently fast to control the attitude and heading of the quadrotor. The flight data are periodically monitored and sent to the ground computer via Xbee communication module.

Mechanical system design

The gyro sensor and accelerometer modules are installed at the center of both helicopters. The controller board interfaces with the accelerometer via RS-232. Angular rate measurement in three axes of the gyro sensor is obtained through ADC module of the microcontroller. The platforms of the two quadrotor helicopters are built as shown in Figure 6. Their parameters are shown in Tables 1 and 2. Moment of inertia along roll, pitch, and yaw axes of the body frame (including components such as the batteries, motors, drivers, and controller board) are obtained from 3-D drawing in SolidWorks [Version 2010]. The moment of inertia of the rotors and motors are obtained by calculation from the motors and rotors size and weight.

Figure 6.

Two quadrotor helicopter platforms.

Table 1.

Large quadrotor helicopter parameters.

Parameter	Value	Unit
Ix	0.087	kg m2
Iy	0.087	kg m2
Iz	0.172	kg m2
J	0.001	kg m2
m	1.82	kg
l	0.4	m

Table 2.

Small quadrotor helicopter parameters.

Parameter	Value	Unit
Ix	0.0175	kg m2
Iy	0.0175	kg m2
Iz	0.0226	kg·m2
J	0.001	kg·m2
m	1.08	kg
l	0.22	m

Electronic boards including microcontroller, communication, and sensor modules are installed as shown in Figures 7 and 8 for the large and small quadrotor helicopters, respectively. Sonar sensor is used for altitude measurement when the helicopter flies at low altitude.

Figure 7.

Electronic boards and sensors of the large quadrotor helicopter.

Figure 8.

Electronic boards and sensors of the small quadrotor helicopter.

Simulation results

To evaluate the performance of the conventional SMC algorithm, as expressed by the dynamic model in equation (26) and control law in equations (32) and (36) to (38) in comparison with the proposed SMC with offset cancellation technique as expressed in equation (63) to (66), simulations are conducted. Since the sensors always contain noises in the output readings, noises are considered in the simulation as well. According to the result from experiment, the measurement noise is zero-mean Gaussian distribution with variance of 0.34. To simulate the real word environment, the noises are therefore added to all measurement readings.

The value of the control parameters is roughly obtained using Lyapunov criteria and fine tuned manually to achieve the maximum control performance. In this work, the same sliding surface and chattering gain parameter are used for the SMC.

Analytic approximation of the control parameter is based on the Lyapunov candidate function in equation (40), and the quadrotor is stable when

$S\dot{S}<0$

77

This relation is valid by introducing a positive constant M where

$\dot{S}=-M\text{sign}\left(S\right)$

78

This yields to

$S\dot{S}=-\left|S\right|M$

79

Apply the sliding surface in equations (27) and (29) into the Lyapunov candidate function in equation (40)

$S\dot{S}=\left({\lambda }_r\left(\varphi -{\varphi }_d\right)+\dot{\varphi }\right)\left({\lambda }_r\dot{\varphi }+\frac{l}{I_x}{u}_2+\frac{\left(I_z-I_y\right)}{I_x}\dot{\theta }\dot{\phi }\right)$

80

Define E as control distance, which is the distance between the required and the actual positions

$E=\left(\varphi -{\varphi }_d\right)$

81

Equation (80) can be rewritten as

$S\dot{S}=E\dot{\varphi }{\lambda }_r^2+\left(E\frac{l}{I_x}{u}_2+e\frac{\left(I_z-I_y\right)}{I_x}\dot{\theta }\dot{\phi }+{\dot{\varphi }}^2\right){\lambda }_r+\left(\frac{l}{I_x}\dot{\varphi }{u}_2+\frac{\left(I_z-I_y\right)}{I_x}\dot{\varphi }\dot{\theta }\dot{\phi }\right)$

82

The control force, u₂, and the control distance, E, must be in opposite direction

$S\dot{S}=E\dot{\varphi }{\lambda }_r^2+\left(-E\frac{l}{I_x}{u}_2+e\frac{\left(I_z-I_y\right)}{I_x}\dot{\theta }\dot{\phi }+{\dot{\varphi }}^2\right){\lambda }_r+\left(-\frac{l}{I_x}\dot{\varphi }{u}_2+\frac{\left(I_z-I_y\right)}{I_x}\dot{\varphi }\dot{\theta }\dot{\phi }\right)$

83

Equation (83) is a quadratic function which has a minimum point at

$\frac{\text{d}}{\text{d}\lambda }\left(S\dot{S}\right)=0$

84

From equation (84), the sliding variable, λ_r, can be determined by

${\lambda }_r=-\frac{\left(-E\frac{l}{I_x}{u}_2+E\frac{\left(I_z-I_y\right)}{I_x}\dot{\theta }\dot{\phi }+{\dot{\varphi }}^2\right)}{2e\dot{\varphi }}$

85

The positive value, M, in equation (79) is the chattering gain which can be obtained by

$0<M\le -\frac{S\dot{S}}{\left|S\right|}$

86

The procedure can be extended to the other control parameters which are determined by the same analytic technique. The maximum angular velocity ( $\dot{\varphi }$ , $\dot{\theta }$ , and $\dot{\phi }$ ), the linear velocity ( $\dot{z}$ ), the control force (u₁, u₂, and u₃), and the control torque (u₄) of the quadrotor are limited by the capability of the rotor and the motor power. They are also limited by the capability of the sensors. The normal operation of the quadrotor must be designed to operate within the sensor and actuator limitation as shown in Table 3 and the control parameters as shown in Table 4.

Table 3.

The maximum operating condition of the quadrotor.

Operating parameter	Quadrotor model		Unit
	Small	Large
$\dot{\varphi }$	1.4	1.4	rad/s
$\dot{\theta }$	1.4	1.4	rad/s
$\dot{\phi }$	1.0	0.25	rad/s
$\dot{z}$	0.25	0.25	m/s
Control distance (E) on roll and pitch	1.57	1.57	rad
Control distance (E) on yaw	6.28	6.28	rad
Control distance (E) on altitude	2.0	2.0	m
Maximum force for roll and pitch control (u2 and u3)	1.0	2.0	N
Maximum torque for yaw control (u4)	0.2	0.3	N·m
Maximum force for altitude control (u1)	12	20	N

Table 4.

The control parameters.

Control parameter	Small model	Large model
λz	2.4	2.6
λr	3.5	2.7
λp	3.5	2.7
λy	4.0	3.0
M 1	0.6	0.6
M 2	2.0	2.0
M 3	2.0	2.0
M 4	2.0	0.8
Cz	1.0	2.0
Cϕ, Cθ	2.7	5.4
Cψ	0.25	0.5

To investigate the influences of the control parameters on the attitude control of the quadrotor, several values of the parameters are applied. The simulation results are shown in Figure 9. From the results, the positive constant M has an influence on the vibration of the quadrotor. With large value of M, the quadrotor vibrates with large amplitude. The sliding variable, λ, has an influence on the response time. With low value of λ, the quadrotor has slow response. In the proposed SMC, the compensation parameter, C, is introduced. With its appropriate value, the attitude response is compensated to have zero (SolidWorks, Version 2010) average error. However, if the value of C is too high, the quadrotor might be unstable.

Figure 9.

The influences of the control parameters.

Figure 10 shows the simulation results of the conventional SMC control when all rotor coefficients of the quadrotor helicopter are identical and nonidentical.

Figure 10.

Simulation results of roll angle when the rotors are identical and nonidentical.

When identical rotors are used, the average value of the roll angle error is only −0.0584°. In contrary, the nonidentical rotors cause the average error rise up to −8.74°.

The comparison of control performance of different controllers is shown in Figure 11. With nonidentical rotors along with noise disturbance to the system, both PID and SMCs can stabilize the quadrotor. However, the system still has an offset error. The average value of the roll angle error with the offset cancellation technique is suppressed only to 0.0616°, while the results of the conventional SMC, the adaptive SMC, and the PID are 1.21°, −1.16°, and −1.45°, respectively.

Figure 11.

Comparison of different controllers.

From both results, fluctuations of the output are observed. These fluctuations are caused from the noises. To eliminate the fluctuation of the output response, a low-pass filter is applied.

The cutoff frequency of the low-pass filter is designed corresponding to the sampling rate. The response of the output when the low-pass filter is applied to the control system is shown in Figure 12. The output with the low-pass filter shows less fluctuation on the output angle.

Figure 12.

Simulation results of roll angle with and without a low-pass filter.

The performance of the proposed SMC with offset cancellation technique is seen from the simulation results as shown in Figure 13. The simulation is based on nonidentical rotor parameters. The conventional SMC results in offset error in the output, while the proposed algorithm successfully controls the output to the desired angle of 0°.

Figure 13.

Performance comparison of the conventional SMC and the proposed SMC with offset cancellation technique. SMC: sliding mode control.

Experimental results

The experiments are conducted on two quadrotor helicopters with nonidentical rotor parameters. Both helicopters are controlled using the same algorithm.

The experimental results in Figures 14 and 15 show the roll and pitch angles from the large quadrotor helicopter when controlled by the conventional SMC compared with the proposed SMC with offset cancellation technique. Due to nonidentical characteristic of the rotors, the results of conventional SMC shows 1.35° offset error in roll angle and 7.47° offset in pitch angle, while the proposed SMC with offset cancellation technique reduces the offset error to −0.25° and 0.08° for roll and pitch angles, respectively.

Figure 14.

Experimental results of roll angle from the large quadrotor helicopter when applying the conventional SMC and the proposed SMC with offset cancellation technique. SMC: sliding mode control.

Figure 15.

Experimental results of pitch angle from the large quadrotor helicopter when applying the conventional SMC and the proposed SMC with offset cancellation technique. SMC: sliding mode control.

From Figures 14 and 15, the experimental results of the conventional SMC stop after 25 s because the offset on row and pitch angles causes the quadrotor becomes unstable and tends to fall; therefore, emergency stop was activated to protect the quadrotor from crash.

Similarly, Figures 16 and 17 show performance comparison for small quadrotor helicopter. The results of conventional SMC show 3.87° and −0.15° for roll and pitch angles, respectively, while the proposed SMC with offset cancellation technique reduces the offset error to −0.09° and −0.14° for roll and pitch angles, respectively.

Figure 16.

Experimental results of roll angle from the small quadrotor helicopter when applying the conventional SMC and the proposed SMC with offset cancellation technique. SMC: sliding mode control.

Figure 17.

Experimental results of pitch angle from the small quadrotor helicopter when applying the conventional SMC and the proposed SMC with offset cancellation technique. SMC: sliding mode control.

The experimental results from the large and the small quadrotors show that the proposed SMC with offset cancellation technique can force the average tracking error to zero regardless of the quadrotor size.

For the altitude and yaw control, the results are shown in Figures 18 and 19, respectively. The altitude is set to 1.0 m above the ground, and the yaw angle is set at 0° from north.

Figure 18.

Experimental results of altitude of SMC with offset cancellation technique from the large quadrotor helicopter. SMC: sliding mode control.

Figure 19.

Experimental results of yaw angle of SMC with offset cancellation technique from the large quadrotor helicopter. SMC: sliding mode control.

To evaluate position control performance of the proposed SMC with offset cancellation technique, two experiments are conducted. The first experiment is on hovering control, while the second experiment is on tracking control. The hovering control result is shown in Figure 20. A fixed target position is provided to the controller. The experimental results show that the quadrotor helicopter successfully hovers around the target position.

Figure 20.

Experimental results of position control of the proposed SMC with offset cancellation technique. SMC: sliding mode control.

The tracking control result is shown in Figure 21. Four target positions are set as the tracking positions as illustrated by points A, B, C, and D, respectively. With the proposed SMC with offset cancellation technique, the quadrotor helicopter successfully tracks along the predefined path. In contrary with the conventional SMC, the quadrotor helicopter flies aimlessly and does not track the predefined path due to the nonidentical rotor parameters.

Figure 21.

Experimental results of tracking control of the proposed SMC with offset cancellation technique. SMC: sliding mode control.

Conclusion

This research proposed SMC with offset cancellation technique to control attitude and altitude of the quadrotor helicopter. Four DC brushless motors are used as the main actuators to generate lift forces to stabilize and control the helicopter motion. The nonlinear dynamic model including the internal and external forces acting on the quadrotor helicopter was derived. The dynamic model was simplified by eliminating the terms with small values such as H-force due to rotor drag and drag force of the frame body. The lift and drag coefficients of the rotor characteristics were identified from the experimental data. The performance of the proposed SMC with offset cancellation technique was then compared with the conventional SMC.

From both simulation and experimental results, the conventional SMC could not remove the offset error, while the proposed SMC with offset cancellation technique could successfully remove the offset error and control the attitude and position of the quadrotor helicopter.