Extended state observer–based nonlinear cascade proportional–integral–derivative control of the nano quadrotor

Abstract

The nano quadrotor is a nonlinear multi-input and multi-output system with strong coupling, which causes difficulties in control law design. In order to achieve a favorable performance, an extended state observer–based nonlinear cascade proportional–integral–derivative controller is proposed in this article. First, the nano quadrotor platform is built, and the dynamic model is established. Second, a novel and practical measuring method is given to obtain model parameters. Then, based on the active disturbance rejection control method, the design procedure of the extended state observer–based nonlinear cascade proportional–integral–derivative controller is presented. In the developed controller, a tracking differentiator is involved to extract the signals of gyroscope, and extended state observer is used to estimate the disturbance. To obtain a better performance of tracking differentiator and extended state observer, a systematic parameter-tuning method is studied. Finally, simulation results are given to demonstrate the efficiency of the proposed controller.

Introduction

The nano quadrotor is usually categorized as a type of light rotary wing unmanned aerial vehicle (UAV), which has quickly emerged in recent years. For its application in research, education and military, nano quadrotor has drawn increasing attention among scholars and geeks.^1,2 Except the advantages of common quadrotors, such as vertical takeoff and landing (VTOL) and hover above target, ³ nano quadrotor can operate in more confined spaces, do more difficult maneuver, and have a better safety. ² Therefore, nano quadrotor has many applications in military and civil use, such as reconnaissance, rescue, and so on.

The above unique features of the nano quadrotor are mainly related to its small scale (around $10\mathrm{cm}\times 10\mathrm{cm}$ ) and lightweight ( $<100\mathrm{g}$ ),^1,2 which will bring difficulties in hardware design, modeling, and control. ⁴ In order to scale down the quadrotor, the body is usually designed in a printed circuit board (PCB), and small brushed motors are used to avoid heavy electronic speed controller (ESC). The shorter control period is demanded because of the agility. Therefore, 32-bit microprogram control unit (MCU) is usually chosen instead of 16-bit MCU. The development of low-cost micro electro mechanical systems (MEMSs) has enabled the miniaturization of the quadrotor, and most six-axis or nine-axis MEMSs are packed just in around $4\mathrm{mm}\times 4\mathrm{mm}$ area, which is suitable for nano quadrotor. The mathematical model of nano quadrotor is similar to common quadrotor. However, the model parameters may be one or two order smaller, ² which cause difficulties in parameter measurement.^5,6 The nano quadrotor can also carry small sensors to obtain more information, such as first person view (FPV) cameras. ² Since the limits of size and weight of the NQP, the resolution of FPV camera is usually low and the image processing method becomes important. In order to analyse color image, the quaternion-type moments were proposed in Chen et al. ⁷ Motion estimation (ME) and disparity estimation (DE) were developed for low complexity multiview video coding, ⁸ and least significant bit (LSB) matching was tackled well by using differences between nonadjacent pixels.^9,10 Besides, the image segmentation suffers from the noise contained in low resolution images. Therefore, more robust algorithm for image segmentation was designed.^11,12

Except the difficulties in modeling, the quadrotor is a classical underactuated system with 6 degrees of freedom (DOFs) but only four inputs, and the dynamics of quadrotor are coupled for different channels, ¹³ which makes it difficult to control. Over the decades, many control strategies have been proposed to stabilize the quadrotor. Proportional–integral–derivative (PID) controller is a classical and widely used control strategy, and many PID control schemes have been designed and applied for quadrotor.^14–16 However, based on the theory that “using error to overcome error,” ¹⁷ conventional PID controller has the disadvantages in overshooting and fast response. An almost globally stabilizing nonlinear proportional–integral–derivative (NPID) was developed by Su et al., ¹⁸ which is more effective and robust than conventional PID. A decentralized proportional integral derivative neural network (PIDNN) control scheme is designed to solve the nonlinearity in quadrotor, and the connective weights of the PIDNN are trained online, which shows the controller has a good adaptability. ²¹ In order to response faster to the input instruction, angular velocity can be controlled in inner loop of PID controller, which is known as cascade proportional–integral–derivative (CPID). ²⁰ CPID has a good damping characteristic and a simple structure, but its setting time may be long. In degrade, the disturbance always exits, which will degrade the performance of the controller. In order to tackle the disturbance and uncertainty, a disturbance observer was designed to provide efficient estimation of the unknown compounded disturbance,^22,23 and adaptive control was applied to achieve a better robustness.^24–27 From the basic principle of classical PID, the active disturbance rejection controller (ADRC) has been proposed, ¹⁷ which consists of tracking differentiator (TD), extended state observer (ESO), and NPID. TD can deal with the conflict between the overshot and rapidity and has an excellent performance in filtering. ESO can estimate the unmodeled dynamics and the unknown disturbances, which can decouple the system. By utilizing the nonlinear feedback, the setting time is short. The ADRC is applied in quadrotor by Liu et al., ²⁸ and comparing with the CPID, the output with ADRC has no oscillation and has a small overshoot. The NPID in ADRC can be chosen in different forms for control needs and control precision, and the difficulty in parameter tuning is different among the forms. ²⁹ Since the ESO can decouple the quadrotor, the controller design can only consider the decoupling system. ³⁰

The controller performance depends on not only the structure of controller but also the parameters of controller. In classical PID, the parameters are mostly tuned by experience. By combining the PID with intelligent control,^31,32 the difficulty in parameter tuning can be overcome. However, the intelligent control is hard to be realized in embedded system ³³ . By classifying the control objects in timescale, ADRC can settle this problem well. ¹⁷ Control objects in similar timescales can share one set of control parameters, which enables the advent of unified parameter-tuning theory. Each parameter in ADRC was analyzed in MATLAB to seek rule of confirming parameters. ³⁴ A systematic parameter-tuning scheme was given in Liu and Shao ³⁵ and applied in asynchronous motor. Considering the similarity of a class of systems, a computation method of timescale in a second-order system was given, and the feasibility of a class of systems using the same ADRC was verified. ³⁶ The controller scaling, parameterizations and optimizations were presented in Gao, ³⁷ and a controller tuning method in frequency domains was summarized. When the parameters are well tuned, the controller can be efficient and has a good antidisturbance performance.

Motivated by the above studies, an extended state observer–based nonlinear cascade proportional–integral–derivative (ENCPID) controller is proposed in this article, which has several modifications listed as follows. The first modification is combining the nonlinear cascade proportional–integral–derivative (NCPID) with the ESO and TD. By changing the NPID in ADRC with NCPID, the dynamics of response will be faster. The second modification is using TD to filter the signal of gyroscope, and prediction is introduced to compensate for the phrase loss. The third is a systematic parameter-tuning procedure that is given based on the previous tuning methods.

This article is organized as follows: First, nano quadrotor platform (NQP) is introduced, and the dynamic model of the nano quadrotor is given. Second, the proposed ENCPID is designed, and a procedure of parameters turning is presented. Finally, a series of simulation results and conclusions are given in the last section.

Platform design and modeling

The self-designed NQP is shown in Figure 1, which consists of the nano quadrotor, the remote control, the ultra wide band (UWB) location system, and the ground station. The nano quadrotor can operate on remote control mode or autonomous mode. In remote control mode, the nano quadrotor receives command from a six-channel remote control and adjusts its flight attitude in real time. In autonomous mode, a preset flight path is designed on ground station and then the quadrotor will track the path autonomously. In both modes, the real-time attitude and position information are displayed on ground station and recorded synchronously for further analysis.

Figure 1.

The self-designed NQP.

The nano quadrotor

The nano quadrotor includes a flight control board, four coreless motors, a lithium battery, and a protective cage. The flight control board is shown in Figure 2, which comprises a nanoprocessor, an inertial measurement unit (IMU), a wireless module interface, and a reserved interface for first-person view (FPV) camera, UWB location system or global positioning system (GPS) module. The nanoprocessor onboard is STM32F405RGT6, with a frequency of 168 MHz, a build-in float point unit (FPU), and a special digital signal processing (DSP) library. An InvenSense IMU MPU9150 is integrated in the center of the board, which contains a tri-axis gyroscope, a tri-axis accelerometer, and a tri-axis magnetometer. Comparing with the separation scheme such as MPU6050 + HMC5883, the MPU9150 has a better coaxality. The wireless interface is connected to a 2.4-GHz Zigbee, and the antenna is printed on the circus board to reduce the weight. The effective communication distance is 150 m, which is enough for the nano quadrotor working in sight. In order to determinate the position, an reserved interface can be connected to one of the three modules: UWB tag, FPV camera, and GPS. The UWB tag and FPV camera are used in GPS-denied area, such as indoor environment, and GPS is equipped in open field.

Figure 2.

The flight control board.

For the balance of torque and mass, four coreless motors are fixed on the flight control board symmetrically, and plastic sleeves are involved to link them. The motors are driven by SI2302 MOSFET to provide lift, and 55-mm rotor is chosen to match the motor. Instead of step-less speed adjustment, more simple and practical pulse width modulation (PWM) control is adopted. To ensure an enough smoothness and accuracy of control, the motor speed range is divided into 500 levels. A 350-mA h lithium battery is fixed below the flight control board to supply electric power. The voltage monitor is designed to avoid over discharging, and max flight time is up to about 7 min. The nano quadrotor is surrounded by a resilient protection cage, which can prevent the damage of collision. The total weight of the nano quadrotor is 43.60 g (without position sensor), and the mass distribution is shown in Figure 3.

Figure 3.

The mass distribution of the nano quadrotor.

System dynamic modeling

To simplify the modeling, some assumptions are made: ³⁸

1. The nano quadrotor is a rigid body where vibration and deformation are ignored.

2. The shape and the distribution of mass are symmetrical.

3. The form and the parameter of the dynamic model are time invariant.

4. The gyroscopic moment of the rotor is ignored.

The nano quadrotor is a 6-degree-of-freedom (DOF) rigid body, including 3-DOF translation and 3-DOF rotation. In order to describe the movement, the body frame (denoted as $R^b$ ) and the navigation frame (denoted as $R^n$ ) are defined in Figure 4.

Figure 4.

The frame of nano quadrotor.

The rotation matrix and transformation matrix from $R^b$ to $R^n$ are given by ³⁹

$$R_b^n=\left[\begin{array}{ccc}{c}_{\theta }{c}_{\psi }& s_{\varphi }{s}_{\theta }{c}_{\psi }-c_{\varphi }{s}_{\psi }& c_{\varphi }{s}_{\theta }{c}_{\psi }+s_{\varphi }{s}_{\psi }\\ c_{\theta }{s}_{\psi }& s_{\varphi }{s}_{\theta }{s}_{\psi }+c_{\varphi }{c}_{\psi }& c_{\varphi }{s}_{\theta }{c}_{\psi }-s_{\varphi }{s}_{\psi }\\ -s_{\theta }& s_{\varphi }{c}_{\theta }& c_{\varphi }{c}_{\theta }\end{array}\right]$$

(1)

$$T=\left[\begin{array}{ccc}1& s_{\varphi }{t}_{\theta }& c_{\varphi }{t}_{\theta }\\ 0& c_{\varphi }& -s_{\varphi }\\ 0& s_{\varphi }/c_{\theta }& c_{\varphi }/c_{\theta }\end{array}\right]$$

(2)

where $s_{*}=\sin {(}^{*})$ , $c_{*}=\cos {(}^{*})$ , and $t_{*}=\tan {(}^{*})$ .

According to the Newton–Euler equation, we have ³⁹

$$\begin{array}{c}{\overrightarrow{F}}^b=m{\stackrel{·}{V}}^b\hfill \\ \multicolumn{1}{c}{{\overrightarrow{\tau }}^b=J{\stackrel{·}{\omega }}^b}\end{array}$$

(3)

where ${\overrightarrow{F}}^b$ is the amount of external force; m is the total mass of the nano quadrotor; $V^b$ are the linear velocities in body frame; ${\overrightarrow{\tau }}^b$ is the resultant external moment; $J=\mathrm{diag}\{J_x,J_y,J_z\}$ is the moment of inertia matrix with $J_x,J_y,J_z$ are the moment of inertia about $x,y,z$ axis, respectively; and ${\omega }^b$ is the angular velocity.

The forces on the nano quadrotor are the gravity, the lift of rotor, and the air drag. The gravity is

$$G=-\mathrm{mg}{\overrightarrow{z}}^b$$

(4)

where g is the gravitational acceleration and ${\overrightarrow{z}}^b$ is the unit vector of the $Z-\mathrm{axis}$ of body frame.

The lift relates to the air density, the aerodynamic shape, and the speed of the rotor. ⁴⁰ Define $T^b=[\begin{array}{cccc}{T}_1^b\hfill & T_2^b\hfill & T_3^b\hfill & T_4^b\hfill \end{array}{]}^T$ as the lift, then we have

$$T_i^b=C_l{\omega }_{\mathrm{ri}}^2,i=1,2,3,4$$

(5)

where $C_l$ is the lift coefficient and ${\omega }_{\mathrm{ri}}$ is the speed of the $i\mathrm{th}$ rotor.

The air drag is given by

$$D^n=-C_d\left|V^n\right|V^n$$

(6)

where $D^n=[\begin{array}{ccc}{D}_x^n\hfill & D_y^n\hfill & D_z^n\hfill \end{array}{]}^T$ is the air drag and $C_d$ is the drag coefficient. Then, the amount of external force can be written as

$$\stackrel{{\rightharpoonup }^n}{F}=-mg\stackrel{{\rightharpoonup }^n}{z}-D^n+T^b{R}_n^b\stackrel{{\rightharpoonup }^n}{z}$$

(7)

The pitching moment is caused by the lift difference of rotors 1 and 2 and rotors 3 and 4, which can be calculated as follows

$${\tau }_{\theta }=l(T_1^b+T_2^b-T_3^b-T_4^b)$$

(8)

where l is the half distance of the adjacent motors.

Similarly, the rolling moment is

$${\tau }_{\varphi }=l(T_2^b+T_4^b-T_1^b-T_3^b)$$

(9)

The yawing moment is determined by the twisting moment difference of rotors 1 and 3 and rotors 2 and 4. The twisting moment of each rotor is written as

$${\tau }_{\mathrm{ti}}=C_t{\omega }_{\mathrm{ri}}^{2}i=1,2,3,4$$

(10)

where $C_t$ is the twisting moment coefficient.

Thus, the yawing moment is

$${\tau }_{\gamma }={\tau }_{t1}+{\tau }_{t3}-{\tau }_{t2}-{\tau }_{t4}$$

(11)

Let U be the control input, and combining equations (5), (8)–(11), we have

$$\begin{array}{c}\multicolumn{1}{c}{U={\left[\begin{array}{cccc}{U}_1\hfill & U_2\hfill & U_3\hfill & U_4\hfill \end{array}\right]}^T}\\ \multicolumn{1}{c}{={\left[\begin{array}{cccc}{T}^b\hfill & {\tau }_{\theta }\hfill & {\tau }_{\varphi }\hfill & {\tau }_{\gamma }\hfill \end{array}\right]}^T}\\ \multicolumn{1}{c}{=\left[\begin{array}{c}{C}_l({\omega }_{r1}^2+{\omega }_{r2}^2+{\omega }_{r3}^2+{\omega }_{r4}^2)\hfill \\ \multicolumn{1}{c}{C_{l}l({\omega }_{r1}^2+{\omega }_{r2}^2-{\omega }_{r3}^2-{\omega }_{r4}^2)}\\ \multicolumn{1}{c}{C_{l}l(-{\omega }_{r1}^2+{\omega }_{r2}^2+{\omega }_{r3}^2-{\omega }_{r4}^2)}\\ \multicolumn{1}{c}{C_t({\omega }_{r1}^2-{\omega }_{r2}^2+{\omega }_{r3}^2-{\omega }_{r4}^2)}\end{array}\right]}\end{array}$$

(12)

Define $P=[x,y,z{]}^T$ as the position in navigation frame and $\vartheta ={[\theta ,\varphi ,\gamma ]}^T$ as the attitude angle in body frame. Then, the dynamic equation of micro quadrotor can be described as

$$\{\begin{array}{c}\multicolumn{1}{c}{\stackrel{··}{x}=\frac{1}{m}{D}_x^n+\frac{U_1}{m}(c_{\varphi }{s}_{\theta }{c}_{\psi }+s_{\varphi }{s}_{\psi })}\\ \multicolumn{1}{c}{\stackrel{··}{y}=\frac{1}{m}{D}_y^n+\frac{U_1}{m}(c_{\varphi }{s}_{\theta }{c}_{\psi }-s_{\varphi }{c}_{\psi })}\\ \multicolumn{1}{c}{\stackrel{··}{z}=-g+\frac{1}{m}{D}_z^n+\frac{U_1}{m}\left(c_{\varphi }{c}_{\theta }\right)}\\ \multicolumn{1}{c}{\stackrel{··}{\theta }=\left(\frac{J_y-J_z}{J_x}\right)\stackrel{·}{\varphi }\stackrel{·}{\gamma }+\frac{U_2}{J_x}}\\ \multicolumn{1}{c}{\stackrel{··}{\varphi }=\left(\frac{J_z-J_x}{J_y}\right)\stackrel{·}{\theta }\stackrel{·}{\gamma }+\frac{U_3}{J_y}}\\ \multicolumn{1}{c}{\stackrel{··}{\gamma }=\left(\frac{J_x-J_y}{J_z}\right)\stackrel{·}{\theta }\stackrel{·}{\varphi }+\frac{U_4}{J_z}}\end{array}$$

(13)

Parameter measurement

The mathematic model of the NQP has been derived in equation (13). However, some of the parameters are still unknown. In this section, some experiments are designed to obtain the model parameter.

The nano quadrotor is an irregular-shaped body. Hence, the moment of inertia cannot be directly calculated. In order to measure the parameter, a trifilar torsional pendulum (TTP) is designed, which is shown in Figure 5.

Figure 5.

Experiment setup for TTP.

In Figure 5, the level meter is used to ensure that the upper and lower disks are in horizontal level, and a smooth link is beneath the lower disk to avoid the side rotation. The structure of the three-line pendulum is shown in Figure 6, where R is the distance between suspension point and disk center, L is the length of suspension line, and $\mu$ is the pendulum angle.

Figure 6.

The structure of TTP.

The moment of inertia around the pendulum axis is given by ⁴¹

$$J_t={\left(\frac{T_p}{2\pi }\right)}^2\frac{(m_0+m_1)g{R}^2}{L}$$

(14)

where $T_p$ is the pendulum period, $m_0$ is the mass of lower disk, $m_1$ is the mass of load, and $J_t$ is the total moment of inertia of the lower disk and load around the pendulum axis.

The mass and the length can be measured precisely by electronic balance and vernier caliper, respectively. However, the precision of $T_p$ is hard to guarantee. In Ringegni et al., ⁴¹ the period was obtained using a Hewlett-Packard 5308a period counter, which was triggered when the small protuberance in pendulum shadows the light, and the precision was limited by the size of the protuberance. In Tang and Shangguan, ⁴² high-precision photoelectronic sensor was involved, but the high-speed camera is much expensive. In this article, the inner timer in the nano quadrotor is chosen to measure the period, which is a 16-bit counter and can operate at 84 MHz to give the time signal. Here, an example is given as follows: The periodic signal along with the time signal is shown in Figure 7.

Figure 7.

The periodic signal with time.

In order to get rid of the noise of magnetometer in Figure 7, a median filter is used, and the result is shown in Figure 8.

Figure 8.

The periodic signal with time after filter.

A binary conversion is made according to the mean of the signal, and the periods are the difference of the time of the rising edges as shown in Figure 9. The mean period in Figure 10 is the T_p.

Figure 9.

The binarization of periodic signal.

Figure 10.

The mean period of the signal.

The center of gravity (CG) can be obtained by suspending the nano quadrotor. ⁴² When the quadrotor is fixed on the lower disk, the CG should be right above the center of the disk.

The measurement of $J_z$ is given as an example:

1. Fix one nano quadrotor on the center of lower disk horizontal.

2. Start measuring the periodic signal (here the magnetic of x-axis is chosen).

3. Give an initial pendulum angle and record the signal for enough time.

4. Intercept reposeful record to obtain $T_p$ and use equation (14) to calculate $J_{z1}$ .

5. Fix two nano quadrotors, and repeat steps 2 to 4 to obtain $J_{z2}$ .

6. Calculate the moment of inertia of the nano quadrotor and disk around z-axis as

$$J_z=J_{z2}-J_{z1}$$

(15)

$$J_{\mathrm{disk}}=J_{z1}-J_z$$

(16)

where the $J_{\mathrm{disk}}$ is the moment of inertia of the disk in TTP. The lift coefficient can be obtained by the experiment in Wang et al. ⁵

In Figure 11, the $\mathrm{lift}\mathrm{-}\mathrm{spee}{\mathrm{d}}^2$ curve is given. The regression coefficient is lift coefficient of the rotor. The measurement of twisting moment is difficult. In Wang et al., ⁵ a cantilever is fixed on motor to measure the twisting and calculates the twisting moment through the pressure on the end of cantilever. When the motor is running at high speed, the pressure on balance is hard to read because of the shaking. In Wang et al., ¹³ the twisting moment is measured by resistance strain gage, which transforms the moment to electrical signal. After the filter and amplifying the signal, the results are satisfactory. However, the strain gage needs calibration before using, and its output is sensitive to temperature. In this article, a simple and useful method is designed to measure the twisting moment. The force analysis of one pendulum point is demonstrated as shown in Figure 12, where H is the vertical distance between upper and lower disks, C is the length that pendulum point moves, $F_L$ is the tension of the line, $G_{\mathrm{total}}$ is the total gravity of the lower disk and load, and $F_T$ is the twisting force.

Figure 11.

The $\mathrm{lift}\mathrm{-}\mathrm{spee}{\mathrm{d}}^2$ curve.

Figure 12.

Force analysis of TTP.

According to the similarity of force triangular and geometric triangular, and considering of the lift of rotor, the proportion relationship exists

$$\frac{H}{C}=\frac{G_{\mathrm{total}}-F_{\mathrm{Lift}}}{F_T}$$

(17)

where $F_{\mathrm{Lift}}$ is the lift of rotors.

The pendulum angle $\mu$ can be measured by IMU, and C can be calculated by

$$C=\sqrt{2R^2(1-\cos \mu )}$$

(18)

And the twisting torque can be calculated as follows

$${\tau }_T=\frac{F_T}{l_1}$$

(19)

In Figure 13, the $\mathrm{twisting}\mathrm{torque}\mathrm{-}\mathrm{spee}{\mathrm{d}}^2$ curve is given. The regression coefficient is twisting moment coefficient of the rotor.

Figure 13.

The $\mathrm{twisting}\mathrm{torque}\mathrm{-}\mathrm{spee}{\mathrm{d}}^2$ curve.

So far, all the model parameters of the nano quadrotor are obtained and listed in Table 1.

Table 1.

The model parameters of the nano quadrotor.

Parameter	Value	Units
$J_x,J_y$	1.916e−5	$\mathrm{kg}{\mathrm{m}}^2$
$J_z$	2.7810e−5	$\mathrm{kg}{\mathrm{m}}^2$
$C_l$	2.478e−08	$\mathrm{N}\mathrm{mi}{\mathrm{n}}^2/\mathrm{ra}{\mathrm{d}}^2$
$C_t$	1.193e−09	$\mathrm{N}\mathrm{mi}{\mathrm{n}}^2/\mathrm{ra}{\mathrm{d}}^2$

Controller design

The principle of ENCPID

Considering the parameter uncertainties and unmodeled disturbances in equation (13), the principle of ENCPID controller is presented in this section.

The framework of the ENCPID flight control system is shown in Figure 14, where ${(}^{*}{)}_d$ is the desired value of ${(}^{*})$ . An ENCPID controller consists of three parts: TD, ESO, and NCPID. The structure is presented in Figure 15. ⁴³

Figure 14.

The framework of ENCPID flight control system.

Figure 15.

The structure of ENCPID.

Take pitch channel as an example and the three parts in Figure 15 are designed as follows:

• 1. The TD can extract the signal and its differential efficiently. Considering that the gyro measured value contains noise and its differential will amplify the noise, the TD is designed as

$$\{\begin{array}{cc}\multicolumn{1}{c}{{\omega }_{\theta 1}(k+1)=}& {\omega }_{\theta 1}\left(k\right)+T{\omega }_{\theta 2}\left(k\right)\hfill \\ \multicolumn{1}{c}{{\omega }_{\theta 2}(k+1)=}& {\omega }_{\theta 2}\left(k\right)+\mathrm{Tfst}({\omega }_{\theta 1}\left(k\right)\hfill \\ \multicolumn{1}{c}{}& -{\omega }_{\theta }\left(k\right),{\omega }_{\theta 2}\left(k\right),r,h)\hfill \end{array}$$

(20)

where ${\omega }_{\theta }$ is the gyro measured value; ${{\omega }_{\theta }}_1$ and ${\omega }_{\theta 2}$ track ${\omega }_{\theta }$ and its differential, respectively; T is the control period; r is the rapid factor; h is the filter factor; and $\mathrm{fst}{(}^{*})$ is defined as ⁴³

$$\begin{array}{c}\mathrm{fst}(x_1,x_2,r,h)=-\{\begin{array}{c}\frac{\mathrm{ra}}{d},a\le d_0\hfill \\ \multicolumn{1}{c}{\mathrm{rsgn}\left(a\right),a>d_0}\end{array}\hfill \\ \multicolumn{1}{c}{\begin{array}{c}d=\mathrm{rh}\hfill \\ \multicolumn{1}{c}{d_0=\mathrm{dh}}\end{array}}\\ \multicolumn{1}{c}{y_0=x_1+h{x}_2,a_0={(d^2+8r\left|y_0\right|)}^{1/2}}\\ \multicolumn{1}{c}{a=\{\begin{array}{c}{x}_2+\frac{a_0-d}{2},\left|y_0\right|>d_0\hfill \\ \multicolumn{1}{c}{x_2+\frac{y_0}{h},\left|y_0\right|\le d_0}\end{array}}\end{array}$$

(21)

• 2. ESO can estimate the states and the disturbance. In equation (13), take pitch channel as example, the relation between $\theta$ and $U_2$ can be written as

$$\begin{array}{c}\stackrel{··}{\theta }=f_{\theta }+\frac{U_2}{J_x}\hfill \end{array}$$

(22)

where $f_{\theta }=((J_y-J_z)/J_x)\stackrel{·}{\varphi }\stackrel{·}{\gamma }+d_{\theta }$ is the total disturbance of pitch channel and $d_{\theta }$ is the external disturbance. The internal disturbance is mainly caused by the coupling and the uncertainty of the parameters, and the external disturbance can also be generated from the wind, the collisions, and so on.

In order to compensate for the disturbance $f_{\theta }$ , the ESO is designed as follows

$${\epsilon }_1=z_1\left(k\right)-\theta \left(k\right)\{\begin{array}{c}{z}_1\left(k+1\right)=z_1\left(k\right)+T\left(z_2\left(k\right)-{\beta }_1{\epsilon }_1\right)\\ z_2\left(k+1\right)=z_2\left(k\right)+T\left(z_3\left(k\right)-{\beta }_{2}fal\left({\epsilon }_1,{\alpha }_{z1},{\delta }_z\right)+b{U}_2\left(k\right)\right)\\ z_3\left(k+1\right)=z_3\left(k\right)-T{\beta }_{3}fal\left({\epsilon }_1,{\alpha }_{z2},{\delta }_z\right)\end{array}$$

(23)

where $\theta (k)$ is the pitch angle; $z_1$ and $z_2$ are the estimation of $\theta (k)$ and its differential, respectively; $z_3$ is the estimation of disturbance $f_{\theta }$ ; ${\alpha }_{z1}$ , ${\alpha }_{z2}$ , ${\delta }_z$ , ${\beta }_1$ , ${\beta }_2$ , and ${\beta }_3$ are parameters to be designed; b is the compensate coefficient; and $\mathrm{fal}{(}^{*})$ is given by ⁴³

$$\mathrm{fal}(\epsilon ,\alpha ,\delta )=\{\begin{array}{c}\left|\epsilon \right|\mathrm{sgn}\left(\epsilon \right),\left|\epsilon \right|>\delta \hfill \\ \multicolumn{1}{c}{\frac{\epsilon }{{\delta }^{1-\alpha }},\left|\epsilon \right|\le \delta }\end{array}\delta >0$$

(24)

Using $z_3$ to compensate equation (22) yields

$$\stackrel{··}{\theta }=\frac{U_{20}}{J_x}$$

(25)

where $U_{20}=U_2+J_x{z}_3/b$ is the middle control variable.

• 3. The NCPID of pitch channel contains two loops, $\theta -\mathrm{loop}$ and ${\omega }_x-\mathrm{loop}$ , which are shown in Figure 16.

Figure 16.

The diagram of NCPID.

Considering that the differential of Euler angles can be obtained directly from the IMU rather than making forward difference, the controller is designed as

$$\{\begin{array}{cc}{U}_{21}=\hfill & k_p\mathrm{fal}(({\theta }_d-\theta ),{\alpha }_{u1},{\delta }_{u1})\hfill \\ \multicolumn{1}{c}{}& +k_d\mathrm{fal}((-{\omega }_1),{\alpha }_{u1},{\delta }_u)\hfill \\ \multicolumn{1}{c}{U_{20}=}& k_{\mathrm{pp}}\mathrm{fal}((U_{21}-{\omega }_1),{\alpha }_{u0},{\delta }_{u0})\hfill \\ \multicolumn{1}{c}{}& +k_{\mathrm{dd}}\mathrm{fal}((-{\omega }_2),{\alpha }_{u0},{\delta }_u)\hfill \end{array}$$

(26)

where $k_p$ , $k_d$ , $k_{\mathrm{pp}}$ , $k_{\mathrm{dd}}$ , ${\alpha }_{u1}$ , ${\delta }_{u1}$ , ${\alpha }_{u0}$ , and ${\delta }_{u0}$ are parameters to be tuned. The controller can be simplified without the nonlinear characteristic of $\mathrm{fal}{(}^{*})$ as

$$\{\begin{array}{c}{U}_{21}=k_p({\theta }_d-\theta )-k_d{\omega }_1\hfill \\ \multicolumn{1}{c}{U_{20}=k_{\mathrm{pp}}(U_{21}-{\omega }_1)-k_{\mathrm{dd}}{\omega }_2}\end{array}$$

(27)

It is noted that ${\omega }_1\approx \stackrel{·}{\theta }$ and ${\omega }_2\approx \stackrel{··}{\theta }$ , and substituting equation (27) into equation (25) yields the closed-loop transfer function of pitch channel

$$\frac{\theta }{{\theta }_d}=\frac{1}{\frac{J_x+k_{\mathrm{dd}}}{k_{\mathrm{pp}}{k}_p}{s}^2+\frac{k_{\mathrm{pp}}{k}_p+k_{\mathrm{pp}}}{k_{\mathrm{pp}}{k}_p}s+1}$$

(28)

Parameters tuning

In this section, the influence of each controller part related to performance is given. Then, the tuning procedure is presented in detail. Finally, a general turning method is summarized.

The parameters of each controller part are of crucial importance to the performance. According to separation principle, ⁴³ the parameters are independent among TD, ESO, and NCPID, which means that the parameters in the three modules can be turned individually.

TD gives the estimate deferential signals ${\omega }_1$ and ${\omega }_2$ to provide damping, and if the rapidity of TD does not satisfy the system, the response will be in weak damping with overshoot. In addition, the noise in deferential signal can cause an unsmooth control input, which causes an adverse impact on the performance. Hence, TD should have a good filter property. ESO estimates the total disturbance including the coupling of the channels to compensate, and the output of ESO determines the ability of resisting noise interference. Both under-compensation and over-compensation should be avoided. NCPID has a decisive influence to the control characteristics. By tuning the NCPID, the pole points of the system can be placed to required area, and the closed-loop performance is satisfied. The parameters to be determined are r; h in TD; ${\alpha }_{z1}$ , ${\alpha }_{z2}$ , ${\delta }_z$ , ${\beta }_1$ , ${\beta }_2$ , and ${\beta }_3$ in ESO; and $k_p$ , $k_d$ , $k_{\mathrm{pp}}$ , and $k_{\mathrm{dd}}$ in NCPID.

The relation between r and transient dynamics $T_d$ is ²⁸

$$r=\frac{4|x_e-x_s|}{T_d^2}$$

(29)

where $x_s$ is the initial value and $x_e$ is the final value.

In equation (29), if we need $T_d$ to be small enough, r should be increased as much as possible to meet the rapidity. h is the filter factor, and when increasing h, the effect of filter is better. However, if h is too large, the phase loss of $v_1$ and $v_2$ cannot be neglected. In practical, h is chosen as ²⁸

$$h=n_{h}T3\le n_h\le 7$$

(30)

where $n_h$ is the filter multiples.

The parameters in TD are highly correlated with the input signal. The fast change dynamic signal requires increasing r and a small phase loss, and much noisy measurement needs a good filter property. Relation equations of parameters of input signals and TD are given by ⁴⁴

$$\{\begin{array}{ccc}{A}_0\approx \frac{r}{0.81{{\omega }_h}^2}<A,\hfill & \frac{r}{{{\omega }_h}^2{\beta }_{h}A}\hfill & <0.81\hfill \\ \multicolumn{1}{c}{A_0\approx {\beta }_{h}A\le A,}& \frac{r}{{{\omega }_h}^2{\beta }_{h}A}\hfill & \ge 0.81\hfill \end{array}$$

(31)

where ${\beta }_h=1/(1+n_h(n_h-1)({\omega }_{h}T{)}^2)\le 1$ ; $A_0$ and A are the amplitudes of the input and the output of TD, respectively; and ${\omega }_h$ is the frequency of input signal.

In equation (31), when $r/({{\omega }_h}^2{\beta }_{h}A)<0.81$ , signal will have an attenuation in amplitude. Therefore, $r/({{\omega }_h}^2{\beta }_{h}A)\ge 0.81$ should be guaranteed to have a good tracking performance.

Increasing $n_h$ will cause a less amplitude loss; however, the phase loss cannot be neglected when $n_h$ is large. Considering that the TD gives the differentiated signal, the phase loss of the filtered signal can be compensated by the following predicting with its differential ¹⁹

$${\omega }_0={\omega }_1+n_{h}h{\omega }_2$$

(32)

where $n_h$ is the steps of prediction; ${\omega }_0$ is the filtered signal after compensation; and ${\omega }_1$ and ${\omega }_2$ are the filtered signal and its differential, respectively.

${\alpha }_{z1}$ , ${\alpha }_{z2}$ , ${\alpha }_1$ , and ${\alpha }_2$ determine the shape of $\mathrm{fal}{(}^{*})$ ; in application, generally, the values are chosen as: ${\alpha }_{z1}=0.5$ , ${\alpha }_{z2}=0.25$ , ${\alpha }_1=0.75$ , and ${\alpha }_2=1.25$ .^28,43

${\delta }_z$ and ${\delta }_{\tau }$ are the width of linear interval of $\mathrm{fal}{(}^{*})$ . If they are too large, the effect of nonsmooth feedback is lost, and when they are too small, chattering will be caused around original point. The reasonable values are ${\delta }_z=0.01$ and ${\delta }_u=0.05$ .

b is the compensation coefficient, and from equation (13), $b=1/J_x$ ideally. However, in practice, b can be tuned to change the compensation proportion. ²⁸

The convergence and stability of ESO are mainly influenced by ${\beta }_1$ , ${\beta }_2$ , and ${\beta }_3$ . In general, when the disturbance is large, the parameters should be increased. When the ESO is well turned, $z_3$ can estimate the total disturbance, and the system will be transformed into a certain series-connection integral object. By turning NCPID, the pole of closed-loop system can be placed to meet the requirement.

Some rules are given as follows: ³⁵

1. In linear ESO, when the control period T is settled, ${\beta }_1$ , ${\beta }_2$ , and ${\beta }_3$ of ESO can be chosen as $1/h$ , $1/(3h^2)$ , and $1/(32h^3)$ .

2. In nonlinear ESO, ${\beta }_1$ can be chosen as $1/h$ first, then tuning ${\beta }_2$ and ${\beta }_3$ alternately. When $z_3$ lags behind the disturbance, try to increase ${\beta }_3$ , and when ESO becomes unstable, increase ${\beta }_2$ properly.

In equation (28), the closed-loop transform function is a second-order system. The adjust time is

$$\{\begin{array}{c}\frac{(k_{\mathrm{pp}}{k}_d+k_{\mathrm{pp}})}{(J_x+k_{\mathrm{dd}})}=2\xi {\omega }_c\hfill \\ \multicolumn{1}{c}{\frac{k_{\mathrm{pp}}{k}_p}{(J_x+k_{\mathrm{dd}})}={\omega }_c^2}\\ \multicolumn{1}{c}{t_s=\frac{4}{\xi {\omega }_c}}\end{array}$$

(33)

where $\xi$ is damping ratio, ${\omega }_c$ is natural frequency, and $t_s$ is the adjust time.

First, the expected performance index $t_s$ should be given, and ${\omega }_c$ can be determined with $\xi$ chosen as the optimal damped ratio $0.6\le \xi \le 0.8$ . Then, proportion between $k_p$ , $k_d$ , $k_{\mathrm{pp}}$ , and $k_{\mathrm{dd}}$ is given by

$$\{\begin{array}{c}\frac{k_p}{1+k_d}=\frac{{\omega }_c}{2\xi }\hfill \\ \multicolumn{1}{c}{\frac{k_{\mathrm{pp}}}{(J_x+k_{\mathrm{dd}})}=\frac{{\omega }_c^2}{k_p}}\end{array}$$

(34)

If $k_{\mathrm{dd}}$ is tuned as $k_{\mathrm{dd}}/J_x\ge 10$ , $J_x$ can be ignored; otherwise, $J_x$ should use the measurement, and fine-tuning $k_{\mathrm{dd}}$ can compensate the measurement error. After pole placement, the system is stable, and the nonlinear characteristic can be added. Then, fine-tune the $k_d$ to overcome the possibly overshoot caused by $\mathrm{fal}{(}^{*})$ .

Simulations

In this section, extensive simulations are presented to evaluate the proposed ENCPID. By applying the tuning method in the previous section, the parameters of ENCPID are given in Table 2.

Table 2.

The parameters of improved ADRC.

Controller part	Parameter	Pitch channel	Roll channel	Yaw channel
TD	r	5,000	5,000	5,000
	h	0.03	0.03	0.03
	$n_h$	2	2	2
ESO	${\beta }_1$	200	200	200
	${\beta }_2$	2,000	2,000	2,000
	${\beta }_3$	13,420	13,420	13,420
	${\alpha }_{z1}$	0.5	0.5	0.5
	${\alpha }_{z2}$	0.25	0.25	0.25
	${\delta }_z$	0.01	0.01	0.01
	b	12	12	12
NCPID	$k_p$	5	5	5
	$k_d$	2.8	2.8	1.8
	$k_{\mathrm{pp}}$	0.0002	0.0002	0.0004
	$k_{\mathrm{dd}}$	0.00002	0.00002	0.00005
	${\alpha }_{u0}$	0.75	0.75	0.75
	${\alpha }_{u1}$	1.25	1.25	1.25
	${\delta }_u$	0.05	0.05	0.05

TD: tracking differentiator; ESO: extended state observer; NCPID: nonlinear cascade proportional–integral–derivative.

Taking pitch channel as example to demonstrate the performance of TD and ESO, the simulation results are shown in Figures 17–23.

Figure 17.

The estimation of angular velocity.

Figure 18.

The estimation of angular acceleration.

Figure 19.

The output of TD in different $n_h$ .

Figure 20.

The output of controller.

Figure 21.

The estimation of pitch.

Figure 22.

The estimation of pitch derivative.

Figure 23.

The estimation of disturbance.

As shown in Figures 17 and 18, the outputs of TD can track the angular velocity and angular acceleration well, and the noise in IMU data is suppressed effectively. In Figure 19, the relation between prediction steps $n_h$ and phrase loss are given, and the larger the $n_h$ , the lesser the phrase loss. When control input is mixed with white noise of high frequency, the output may be smooth because of the inertia of the system; however, it may give an impact on the motor. Figure 20 demonstrates that the noise in control input can be handled by introducing TD. Figures 21–23 show that the three outputs of ESO can track pitch, the derivative of pitch, and the total disturbance well.

In order to demonstrate the performance of the proposed ENCPID, the simulation results of ADRC, CPID, and ENCPID are shown in Figures 24–26. The adjusting time is ENCPID ≈ ADRC > CPID, which is because of the nonlinear function $\mathrm{fal}{(}^{*})$ in ENCPID and ADRC. However, in 1–1.5 s, the ENCPID and CPID response faster than ADRC because the transition brought by TD exists in ADRC, and the $\omega -\mathrm{loop}$ in ENCPID and CPID can guarantee the rapid response of the input instruction.

Figure 24.

The step response of pitch channel.

Figure 25.

The step response of roll channel.

Figure 26.

The step response of yaw channel.

Conclusion

In this article, the attitude control of the NQP using the ENCPID is studied. A nano quadrotor platform, including the nano quadrotor, the remote, and the ground station, is designed first. For further research, the mathematical model is deduced, and its parameters are measured precisely by several experiments. Then, the ENCPID controller is proposed, and a systematic parameter-tuning method is summarized. Finally, the simulation results show that the ENCPID has a rapid response over the adjusting time and a good noise rejection capability. The nano quadrotor platform built in this article has a great potential in deeper research, and some application, such as trajectory control and autonomous cooperative control, will be studied in the future.