Dynamic Stiffness Enhancement of the Quadcopter Control System

Abstract

Vibrations often result in fatigue breakage in quadcopter flight operations. Reducing the vibration effect is the main important issue for quadcopter flight. Enhanced dynamic stiffness is required in the quadcopter control system for the vibration rejection ability. In this paper, the dynamic stiffness of the quadcopter control system is constructed and is used as an index for the performance of the quadcopter control system to resist an external oscillatory load. To rapidly reduce the vibration effect, a repetitive controller is introduced in the quadcopter control system, where the direct dynamic stiffness and the quadrature dynamic stiffness with variable frequencies are proposed. A theoretical model of the dynamic stiffness of the quadcopter control system is established by analyzing the definition of dynamic stiffness. Simulated and experimental results show that the magnitudes of the direct dynamic stiffness with the proposed repetitive controller are larger than those without the repetitive controller. The quadrature dynamic stiffness with the proposed repetitive controller is relatively smooth compared to that without the proposed repetitive controller, which can be used to verify the stability being improved by the designed repetitive controller. In addition, the magnitudes of the loss factor of the quadcopter control system with the repetitive controller are lower than those without the repetitive controller.

1. Introduction

Quadcopters have become a research topic in recent years because of their special functionality and industrial requirements. Quadcopters can not only increase the freedom of movement but also can be remotely controlled, so industrial applications have increased, such as for surveillance, search and rescue, building inspection, and some other applications. The quadcopter system has 6 degrees of freedom of movement in space (translation and rotation along three coordinate axes) but only 4 control degrees of freedom (the rotation speed of the four motors), which is called an under-actuated control system. It is well known that a complete actuated system is when the control degrees of freedom are equal to the movement degrees of freedom. However, for the attitude control (rotational motion along three axes), it is indeed fully driven. Compared with helicopters, quadcopters can achieve fewer flight attitudes, but basic forward, backward, translation and other states can be operated. However, the mechanical structure of the quadcopter is far simpler than that of the helicopter, and the maintenance and substitution costs are also very small, which gives the quadcopter a greater application advantage than the helicopter [1,2].

In order to achieve a stable flight operator, the quadcopter is supplied with three-directional gyroscopes and three-axis acceleration sensors to form an inertial navigation module, which can reckon the quadcopter’s attitude, acceleration, and angular velocity relative to the ground. The flight controller can be used to obtain the rotational force and lift required to maintain the motion state through an algorithm and uses the electronic controller to ensure that the motor can output the appropriate force.

With two pairs of motors in opposite directions, the counter-torque to the fuselage can be balanced. For the output power of the four motors to be increased, the rotor speed and the total pulling force can be increased. The total pulling force should be enough to supply the whole quadcopter’s weight. The rotorcraft will rise vertically from the ground; on the contrary, the output power of the four motors will be reduced. In addition, the quadrotor will fall vertically until it lands in a balanced manner, realizing vertical movement along the z-axis. Without the external disturbance, the lift generated by the rotor is equal to the quadcopter’s weight, and the quadcopter will remain in a hovering state. The key to vertical motion is to ensure that the rotation speed of the four rotors increases or decreases synchronously.

At present, the most common algorithm in the application of quadcopters is to use complementary filtering: that is, to calculate the angle change by combining the output of the acceleration sensor and the gyroscope. In addition, the environment of the quadcopter operation determines the MEMS sensor, which must be suitable for various harsh conditions while obtaining high-precision output. For example, the ideal output of a gyroscope is to only respond to changes in angular velocity, but in fact, due to the limitations of design and craftsmanship, gyroscopes are also sensitive to acceleration, which is the deg/sec/g indicator on gyroscope datasheets. For the application of a multi-rotor aircraft, this indicator is particularly important, because the motor in the aircraft generally brings relatively strong vibration. The change in the output of the gyroscope will cause the change in the angle, and the motor will malfunction. With the presence of external factors, vibration can be found in a quadcopters’ flight to cause a non-uniform flow field. In practice, the actual quadcopter lift force is made up of harmonics related to the angular velocity of the quadcopter propeller.

Many related works in the literature can be found to handle the problems of the disturbance, uncertainties, and motor faults for the quadcopter control system. In [3], the issue of the quadcopter with vibration mitigation is proposed. In [4], an adaptive sliding mode control is proposed to solve the problem of uncertainty and external disturbances for quadcopter flight. In [5], the detection and diagnosis of motor and propeller degradation is raised to reduce the vibration of the quadcopter operation. In [6], to solve the problem of motor faults, the fault-tolerant control is addressed. In [7], a model predictive controller and the disturbance observer are given to stable the quadcopter. In [8], a robust hybrid nonlinear control is applied to the quadcopter with the cross-coupling disturbance. In [9], to reduce the effect of load uncertainties, an observer-based attitude stabilization mechanism is proposed. In [10], an aerodynamic model is proposed to simulate the model of the external disturbance of the quadcopter. In [11], with parametric uncertainties and external disturbances, an adaptive fast finite-time control is adapted for a tilting quadcopter. In [12], a cooperative path following control is given to improve the stability of the quadcopter with unknown external disturbances. In [13], a bidirectional fuzzy brain emotional learning controller is proposed to solve the problem of the payload uncertainties and disturbances. In [14], a robust control with a disturbance observer is considered to improve the stability of the quadcopter. In [15], a finite-time control method for a multirotor UAV is proposed with parameter uncertainties and external disturbances. In [16], a fault-tolerant control is given to the damaged propellers. In [17], with the vibration of the multi-rotor arms, a fault detection using artificial intelligence is proposed. In [18], wind measurement and simulation methods for multi-rotor UAVs are proposed to achieve a stable operation. In [19], with motor fault and external disturbance, an adaptive attitude control for multi-rotor UAVs is given. In [20], a disturbance observer-based controller is proposed with parametric uncertainties and wind disturbance. Therefore, according to related works in the literature, how to reduce the vibration generated by external disturbance, wind disturbance, and motor/propeller faults is the key point for the quadcopter operation.

In general, the ability of lower frequency disturbance rejection can be regarded as static stiffness, which can be achieved by a high mathematical design [21]. Therefore, to achieve better performance of the quadcopter control system, the dynamic stiffness [21,22,23] should be enhanced. In this paper, a methodology to analyze the dynamic stiffness of the quadcopter control system is proposed. The dynamic stiffness is used to be an index for the performance of the quadcopter control system to resist an external oscillatory load. The proposed dynamic stiffness is varying under excitations at different frequencies. The outline of this paper is developed as follows. Section 2 contains the mathematical modeling of the quadcopter. In Section 3, the dynamic stiffness of the quadcopter control system with a repetitive controller is proposed. Simulations and experimental results are performed in Section 4. Finally, concluding remarks are stated in Section 5.

2. System Description and Modeling

The quadcopter’s modeling is shown in Figure 1, where the four motors are used to force the machine. The inertial frame $\mathsf{\xi }=\left(x,y,z\right)$ with three axes of x, y, and z is defined, while the angular position can be obtained as $\mathsf{\eta }=\left(\varphi ,\theta ,\phi \right)$. The three Euler angles of $\varphi$ (defined as roll angle), $\theta$ (defined as pitch angle), and $\psi$ (defined as yaw angle) are represented as the rotations around the x, y, and z-axis, respectively. A vector $q=\left(\mathsf{\xi },\mathsf{\eta }\right)$ is obtained by the linear vector of $\mathsf{\xi }=\left(x,y,z\right)$ and the angular position vector of $\mathsf{\eta }=\left(\varphi ,\theta ,\phi \right)$. In the frame, the linear velocities of $V=\left(v_x,v_y,v_z\right)$ and the angular velocities of $\mathsf{\nu }=\left(p,q,r\right)$ are defined. The rotation matrix of $R$ from the body frame to the inertial frame can be determined by

$$R=\left[\begin{array}{ccc}{\mathrm{c}}_{\phi }{\mathrm{c}}_{\theta }& {\mathrm{c}}_{\phi }{\mathrm{s}}_{\theta }{\mathrm{s}}_{\varphi }-{\mathrm{s}}_{\phi }{\mathrm{c}}_{\varphi }& {\mathrm{c}}_{\phi }{\mathrm{s}}_{\theta }{\mathrm{c}}_{\varphi }+{\mathrm{s}}_{\phi }{\mathrm{s}}_{\varphi }\\ {\mathrm{s}}_{\phi }{\mathrm{c}}_{\theta }& {\mathrm{s}}_{\phi }{\mathrm{s}}_{\theta }{\mathrm{s}}_{\varphi }+{\mathrm{c}}_{\phi }{\mathrm{c}}_{\varphi }& {\mathrm{s}}_{\phi }{\mathrm{s}}_{\theta }{\mathrm{c}}_{\varphi }-{\mathrm{c}}_{\phi }{\mathrm{s}}_{\varphi }\\ -{\mathrm{s}}_{\theta }& {\mathrm{c}}_{\theta }{\mathrm{s}}_{\varphi }& {\mathrm{c}}_{\theta }{\mathrm{c}}_{\varphi }\end{array}\right]$$

(1)

where $s_i=\sin i$ and $c_i=\cos i$ for $i=\varphi ,\theta ,\phi$ are defined.

$W_{\eta }$ is defined as the transformation matrix for angular velocities from the inertial frame to the body frame. Then, we have

$$W_{\eta }=\left[\begin{array}{ccc}1& 0& -{\mathrm{s}}_{\theta }\\ 0& {\mathrm{c}}_{\varphi }& {\mathrm{c}}_{\theta }{\mathrm{s}}_{\varphi }\\ 0& -{\mathrm{s}}_{\varphi }& {\mathrm{c}}_{\theta }{\mathrm{c}}_{\varphi }\end{array}\right]$$

(2)

and

$$\mathsf{\nu }=\left[\begin{array}{c}p\\ q\\ r\end{array}\right]=W_{\eta }\left[\begin{array}{c}\frac{d\varphi }{dt}\\ \frac{d\theta }{dt}\\ \frac{d\phi }{dt}\end{array}\right]$$

(3)

Assume that the four arms of the quadcopter are aligned to the x- and y-axes of the device body. Thus, the inertia matrix of $I=diag\left(I_{xx},I_{yy},I_{zz}\right)$ is a diagonal matrix with $I_{xx}=I_{yy}$. The angular velocity of the rotor axis is given as ${\omega }_i$, and $f_i$ is the generated force in the direction of the rotor axis. ${\tau }_{M_i}$ is the torque of the rotor axis, where $i=1,2,3,4$ are for the four-rotor axis. Then, we have $f_i=k{\omega }_i{}^2$ and ${\tau }_{M_i}=b{\omega }_i{}^2+I_M\frac{d{\omega }_i}{dt}$, where $k$ is the lift constant, $b$ is the drag constant, and the rotor inertia moment is defined as $I_M$. In the z-axis of the body frame, the combined thrust of rotors is given as $T$. The torque of ${\mathsf{\tau }}_B=\left(\begin{array}{ccc}{\tau }_{\varphi }& {\tau }_{\theta }& {\tau }_{\phi }\end{array}\right)$ is defined in the direction of $\left(\varphi ,\theta ,\phi \right)$. Therefore, we have

$$T={\displaystyle \sum _{i=1}^4{f}_i}=k{\displaystyle \sum _{i=1}^4{\omega }_i^2}$$

(4)

and

$${\mathsf{\tau }}_B=\left[\begin{array}{c}{\tau }_{\varphi }\\ {\tau }_{\theta }\\ {\tau }_{\phi }\end{array}\right]=\left[\begin{array}{c}lk\left(-{\omega }_2{}^2+{\omega }_4{}^2\right)\\ lk\left(-{\omega }_1{}^2+{\omega }_3{}^2\right)\\ {\displaystyle \sum _{i=1}^4{\tau }_{M_i}}\end{array}\right]$$

(5)

where $l$ is the distance between the mass center and the rotor.

Therefore, the dynamics characteristic is determined by Newton–Euler equations. In the body frame, $m\frac{d{V}_B}{dt}$ is defined as the force for the mass acceleration, $\mathsf{\nu }\times m{V}_B$ is defined by the centrifugal force being equal to the gravity of $R^{T}G$. Then, the whole thrust $T_B$ generated by the rotors motor can be obtained by

$$m\frac{d{V}_B}{dt}+\mathsf{\nu }\times m{V}_B=R^{T}G+T_B$$

(6)

Assume that the centrifugal force is neglected. Thus, (7) can be simplified as

$$m\frac{d^2\mathsf{\xi }}{d{t}^2}=G+R{T}_B$$

(7)

In the body frame, the angular acceleration of the inertia can be obtained by $I\frac{d\mathsf{\nu }}{dt}$, the centripetal forces and the gyroscopic forces are determined by $\mathsf{\nu }\times \left(I\mathsf{\nu }\right)$ and $\Gamma$, respectively, and $\mathsf{\tau }$ is represented by the external torque. Therefore, we have

$$I\frac{d\mathsf{\nu }}{dt}+\mathsf{\nu }\times \left(I\mathsf{\nu }\right)+\Gamma =\mathsf{\tau }$$

(8)

The angular accelerations in the inertial frame can be given as

$$\frac{d^2\mathsf{\eta }}{d{t}^2}=\frac{d\left(W_{\eta }{}^{-1}\mathsf{\nu }\right)}{dt}=\frac{d{W}_{\eta }{}^{-1}}{dt}\mathsf{\nu }+W_{\eta }{}^{-1}\frac{d\mathsf{\nu }}{dt}$$

(9)

The Lagrangian $\Lambda$ can be given as

$$\Lambda \left(q,\frac{dq}{dt}\right)=\frac{m}{2}{\left(\frac{d\mathsf{\xi }}{dt}\right)}^T\frac{d\mathsf{\xi }}{dt}+\frac{{\mathsf{\nu }}^{T}I\mathsf{\nu }}{2}-mgz$$

(10)

The external forces and torques of the device can be determined by Euler–Lagrange equations, i.e.,

$$\left[\begin{array}{c}f\\ \mathsf{\tau }\end{array}\right]=\frac{d}{dt}\left(\frac{\partial \Lambda }{\partial \left(\frac{dq}{dt}\right)}\right)-\frac{\partial \Lambda }{\partial q}$$

(11)

The Jacobian matrix of $J\left(\mathsf{\eta }\right)$ from $\mathsf{\nu }$ to $\frac{d\mathsf{\eta }}{dt}$ is

$$J\left(\mathsf{\eta }\right)=W_{\eta }{}^{T}I{W}_{\eta }=\left[\begin{array}{ccc}{I}_{xx}& 0& -I_{xx}{\mathrm{s}}_{\theta }\\ 0& I_{yy}{\mathrm{c}}_{\phi }^2+I_{zz}{\mathrm{s}}_{\varphi }^2& \left(I_{yy}-I_{zz}\right){\mathrm{c}}_{\varphi }{\mathrm{s}}_{\varphi }{\mathrm{c}}_{\theta }\\ -I_{xx}{\mathrm{s}}_{\theta }& \left(I_{yy}-I_{zz}\right){\mathrm{c}}_{\varphi }{\mathrm{s}}_{\varphi }{\mathrm{c}}_{\theta }& I_{xx}{\mathrm{s}}_{\theta }^2+I_{yy}{\mathrm{s}}_{\varphi }^2{\mathrm{c}}_{\theta }^2+I_{zz}{\mathrm{c}}_{\varphi }^2{\mathrm{c}}_{\theta }^2\end{array}\right]$$

(12)

Thus, in the inertial frame, the rotational energy $E_{rot}$ is obtained by

$$E_{rot}=\frac{1}{2}{\mathsf{\nu }}^{T}I\mathsf{\nu }=\frac{1}{2}{\left(\frac{d^2\mathsf{\eta }}{d{t}^2}\right)}^{T}J\left(\frac{d^2\mathsf{\eta }}{d{t}^2}\right)$$

(13)

The external angular force is the torque outputs of the rotor motors, i.e.,

$$\mathsf{\tau }={\mathsf{\tau }}_B=J\frac{d^2\mathsf{\eta }}{d{t}^2}+\frac{dJ}{dt}\frac{d\mathsf{\eta }}{dt}-\frac{1}{2}\frac{\partial }{\partial \mathsf{\eta }}\left[{\left(\frac{d\mathsf{\eta }}{dt}\right)}^{T}J\left(\frac{d\mathsf{\eta }}{dt}\right)\right]=J\frac{d^2\mathsf{\eta }}{d{t}^2}+C\left(\mathsf{\eta },\frac{d\mathsf{\eta }}{dt}\right)\frac{d\mathsf{\eta }}{dt}$$

(14)

where the matrix $C\left(\mathsf{\eta },\frac{d\mathsf{\eta }}{dt}\right)$ being the Coriolis term contains the gyroscopic and centripetal terms. Then, we have

$$C\left(\mathsf{\eta },\frac{d\mathsf{\eta }}{dt}\right)=\left[\begin{array}{ccc}0& \left[\begin{array}{l}\left(I_{yy}-I_{zz}\right)\left(\frac{d\theta }{dt}{\mathrm{c}}_{\varphi }{\mathrm{s}}_{\varphi }+\frac{d\phi }{dt}{\mathrm{s}}_{\varphi }^2{\mathrm{c}}_{\theta }\right)\\ +\left(I_{zz}-I_{yy}\right)\left(\frac{d\phi }{dt}{\mathrm{c}}_{\varphi }^2{\mathrm{c}}_{\theta }-\frac{d\phi }{dt}{I}_{xx}{\mathrm{c}}_{\theta }\right)\end{array}\right]& \left(I_{zz}-I_{yy}\right)\frac{d\phi }{dt}{\mathrm{c}}_{\varphi }{\mathrm{s}}_{\varphi }{c}_{\theta }^2\\ \left[\begin{array}{l}\left(I_{zz}-I_{yy}\right)\left(\frac{d\theta }{dt}{\mathrm{c}}_{\varphi }{\mathrm{s}}_{\varphi }+\frac{d\phi }{dt}{\mathrm{s}}_{\varphi }{\mathrm{c}}_{\theta }\right)\\ +\left(I_{yy}-I_{zz}\right)\left(\frac{d\phi }{dt}{\mathrm{c}}_{\varphi }^2{\mathrm{c}}_{\theta }+\frac{d\phi }{dt}{I}_{xx}{\mathrm{c}}_{\theta }\right)\end{array}\right]& \left(I_{zz}-I_{yy}\right)\frac{d\varphi }{dt}{\mathrm{c}}_{\varphi }{\mathrm{s}}_{\varphi }& \left[\begin{array}{l}-I_{xx}\frac{d\phi }{dt}{\mathrm{s}}_{\theta }{\mathrm{c}}_{\theta }+I_{yy}\frac{d\phi }{dt}{\mathrm{s}}_{\phi }^2{\mathrm{s}}_{\theta }{\mathrm{c}}_{\theta }\\ +I_{zz}\frac{d\phi }{dt}{\mathrm{c}}_{\phi }^2{\mathrm{s}}_{\theta }{\mathrm{c}}_{\theta }\end{array}\right]\\ \left(I_{yy}-I_{zz}\right)\frac{d\phi }{dt}{\mathrm{c}}_{\theta }^2{\mathrm{s}}_{\varphi }{\mathrm{c}}_{\phi }-I_{xx}\frac{d\theta }{dt}{\mathrm{c}}_{\theta }& \left[\begin{array}{l}\left(I_{zz}-I_{yy}\right)\left(\frac{d\theta }{dt}{\mathrm{c}}_{\varphi }{\mathrm{s}}_{\varphi }{\mathrm{s}}_{\theta }+\frac{d\varphi }{dt}{\mathrm{s}}_{\varphi }^2{\mathrm{c}}_{\theta }\right)\\ +\left(I_{yy}-I_{zz}\right)\left(\frac{d\varphi }{dt}{\mathrm{c}}_{\varphi }^2{\mathrm{c}}_{\theta }\right)+I_{xx}\frac{d\phi }{dt}{\mathrm{s}}_{\theta }{\mathrm{c}}_{\theta }\\ -I_{yy}\frac{d\phi }{dt}{\mathrm{s}}_{\varphi }^2{\mathrm{s}}_{\theta }{\mathrm{c}}_{\theta }-I_{zz}\frac{d\phi }{dt}{\mathrm{c}}_{\varphi }^2{\mathrm{s}}_{\theta }{\mathrm{c}}_{\theta }\end{array}\right]& \left[\begin{array}{l}\left(I_{yy}-I_{zz}\right)\frac{d\varphi }{dt}{\mathrm{c}}_{\varphi }{\mathrm{s}}_{\varphi }{c}_{\theta }^2+I_{xx}\frac{d\theta }{dt}{\mathrm{s}}_{\theta }{\mathrm{c}}_{\theta }\\ -I_{yy}\frac{d\theta }{dt}{\mathrm{s}}_{\varphi }^2{\mathrm{s}}_{\theta }{\mathrm{c}}_{\theta }-I_{zz}\frac{d\theta }{dt}{\mathrm{c}}_{\varphi }^2{\mathrm{s}}_{\theta }{\mathrm{c}}_{\theta }\end{array}\right]\end{array}\right]$$

(15)

From (17), we have

$$\frac{d^2\mathsf{\eta }}{d{t}^2}=J^{-1}\left({\mathsf{\tau }}_B-C\left(\mathsf{\eta },\frac{d\mathsf{\eta }}{dt}\right)\frac{d\mathsf{\eta }}{dt}\right)$$

(16)

The simplified model of the quadcopter can be obtained as

$$m\frac{d^{2}X}{d{t}^2}+D\frac{dX}{dt}+mG=T$$

(17)

where $X={\left[\begin{array}{ccc}x& y& z\end{array}\right]}^T$, $D=diag\left(D_x,D_y,D_z\right)$, $G={\left[\begin{array}{ccc}0& 0& g\end{array}\right]}^T$, and $T=\left[\begin{array}{c}\left({\mathrm{c}}_{\phi }{\mathrm{s}}_{\theta }{\mathrm{c}}_{\varphi }+{\mathrm{s}}_{\phi }{\mathrm{s}}_{\varphi }\right)T\\ \left({\mathrm{s}}_{\phi }{\mathrm{s}}_{\theta }{\mathrm{c}}_{\varphi }-{\mathrm{c}}_{\phi }{\mathrm{s}}_{\varphi }\right)T\\ \left({\mathrm{c}}_{\theta }{\mathrm{c}}_{\varphi }\right)T\end{array}\right]$ are given. $D_x$, $D_y$ and $D_z$ are the drag force coefficients for velocities in $\left(x,y,z\right)$.

3. Dynamic Stiffness Analysis of the Quadcopter Control System

A closed loop control system is shown in Figure 2, where $G\left(s\right)$ and $K\left(s\right)$ are the controlled plant and loop controller, respectively. In general, $d$ is the system vibration disturbance, $y$ is the system output, $u$ is the controller output, and $e$ is the error. Note that $P\left(s\right)=K\left(s\right)G\left(s\right)$. Let the frequency responses of $P\left(j\omega \right)=P_1\left(\omega \right)+j{P}_2\left(\omega \right)$ and $K\left(j\omega \right)=C_1\left(\omega \right)+j{C}_2\left(\omega \right)$ be given, where we have

$$[\begin{array}{l}{P}_1\left(\omega \right)=G_1{C}_1\left(\omega \right)-G_2{C}_2\left(\omega \right)\\ P_2\left(\omega \right)=G_1{C}_2\left(\omega \right)+G_1{C}_2\left(\omega \right)\end{array}$$

(18)

Therefore, the dynamic stiffness of Figure 2 can be denoted by $D\left(s\right)=1+P\left(s\right)$. We have $D\left(j\omega \right)=D_1\left(\omega \right)+j{D}_2\left(\omega \right)$, and ${\left[D\left(j\omega \right)\right]}^{-1}$ is defined as the dynamic compliance. Note that $D_1\left(\omega \right)=1+P_1\left(\omega \right)$ and $D_2\left(\omega \right)=P_2\left(\omega \right)$ are defined as direct dynamic stiffness and the quadrature dynamic stiffness, respectively. Note that $D_1\left(\omega \right)$ can be represented as the conservative properties of the system, and $D_2\left(\omega \right)$ can be represented as the dissipative properties. Then, the magnitude of the dynamic stiffness $D\left(j\omega \right)$ can be given as $\left|D\left(j\omega \right)\right|=\sqrt{{\left[1+P_1\left(\omega \right)\right]}^2+{\left[P_2\left(\omega \right)\right]}^2}$, which can be rewritten as $D\left(j\omega \right)=D_1\left(\omega \right)\left[1+jH\left(\omega \right)\right]$, and $H\left(\omega \right)$ can be defined as the loss factor of the control system, i.e.,

$$H\left(\omega \right)=\frac{D_2\left(\omega \right)}{D_1\left(\omega \right)}=\frac{P_2\left(\omega \right)}{1+P_1\left(\omega \right)}$$

(19)

For the required performance, it can be expected that $D_1\left(\omega \right)$ is increasing with $\omega$ and $H\left(\omega \right)$ is decreasing with $\omega$. From the above analysis, the dynamic stiffness of the closed-loop control system in Figure 3 is indeed more than that without a loop controller.

For the vibration rejection, a repetitive control system is induced, as shown in Figure 3. Assume that the vibration disturbance $d$ is a single period of $T_d$. The dynamic stiffness of Figure 3 can be denoted by

$$D_r\left(s\right)=1+R\left(s\right)P\left(s\right)=1+\frac{P\left(s\right)}{1-e^{-s{T}_d}}.$$

(20)

Then, the frequency response of $D_r\left(j\omega \right)$ can be given as $D_r\left(j\omega \right)=D_{r1}\left(\omega \right)+j{D}_{r2}\left(\omega \right)$. Note that $D_{r1}\left(\omega \right)=1+P_1\left(\omega \right){\displaystyle \sum _{k=0}^n{\mathrm{c}}_k}+P_2\left(\omega \right){\displaystyle \sum _{k=0}^n{\mathrm{s}}_k}$ and $D_{r2}\left(\omega \right)=P_2\left(\omega \right){\displaystyle \sum _{k=0}^n{\mathrm{c}}_k}-P_1\left(\omega \right){\displaystyle \sum _{k=0}^n{\mathrm{s}}_k}$ are represented as the direct dynamic stiffness and the quadrature dynamic stiffness, respectively, where $c_k=\cos \left(k\omega T_d\right)$ and $s_k=\sin \left(k\omega T_d\right)$ are defined.

The loss factor $H_r\left(\omega \right)$ of Figure 3 can be given as

$$H_r=\frac{D_{r2}\left(\omega \right)}{D_{r1}\left(\omega \right)}=\frac{P_2\left(\omega \right){\displaystyle \sum _{k=0}^n{\mathrm{c}}_k}-P_1\left(\omega \right){\displaystyle \sum _{k=0}^n{\mathrm{s}}_k}}{1+P_1\left(\omega \right){\displaystyle \sum _{k=0}^n{\mathrm{c}}_k}+P_2\left(\omega \right){\displaystyle \sum _{k=0}^n{\mathrm{s}}_k}}$$

(21)

As $\omega =\frac{2k\pi }{T_d}=k{\omega }_d$, $k\in {\rm N}$, the loss factor $H_r\left(k{\omega }_d\right)$ can be given as

$$H_r\left(k{\omega }_d\right)=\frac{\left(1+n\right)P_2\left(k{\omega }_d\right)}{1+\left(1+n\right)P_1\left(k{\omega }_d\right)}$$

(22)

Due to the poles of $R\left(s\right)=\frac{1}{1-e^{-s{T}_d}}$ being located at an imaginary axis, to enhance the stability of the repetitive control system in Figure 3, a low-pass filter $Q\left(s\right)$ is induced as shown in Figure 4, where $\Vert Q\left(j\omega \right)\Vert \approx 1$ for $\omega \le {\omega }_q$ is given.

The dynamic stiffness of Figure 4 can be denoted by

$$D_{qr}\left(s\right)=1+R_q\left(s\right)P\left(s\right)=1+\frac{P\left(s\right)}{1-e^{-s{T}_d}Q\left(s\right)}$$

(23)

Its frequency response can be given as

$$D_{qr}\left(j\omega \right)=1+P\left(j\omega \right)\left[1+{}_{q}e^{-j\omega T_d}+{}_{q}e^{-2j\omega T_d}+\dots \right]$$

(24)

where ${}_{q}e^{-jn\omega T_d}=e^{-jn\omega T_d}{\left[Q\left(j\omega \right)\right]}^n$, $n=1,2,3,\dots$ is defined.

Note that $Q\left(j\omega \right)=Q_1\left(\omega \right)+j{Q}_2\left(\omega \right)$ is given.

In addition, we have

$$[\begin{array}{l}{D}_{qr}\left(j\omega \right)=1+\left(P_1\left(\omega \right)+j{P}_2\left(\omega \right)\right)\left[1+{\displaystyle \sum _{k=1}^n{\left[e^{-j\omega T_d}Q\left(j\omega \right)\right]}^k}\right]\\ =1+\left(P_1\left(\omega \right)+j{P}_2\left(\omega \right)\right)\left[{\displaystyle \sum _{k=0}^n{\left[\left(\cos \left(\omega T_d\right)-j\sin \left(\omega T_d\right)\right)\left(Q_1\left(\omega \right)+j{Q}_2\left(\omega \right)\right)\right]}^k}\right]\\ =1+{\displaystyle \sum _{k=0}^n\left(P_1\left(\omega \right)+j{P}_2\left(\omega \right)\right){\left(\tilde{c}-j\tilde{s}\right)}^k}\end{array}$$

(25)

where $\tilde{c}=c_1{Q}_1\left(\omega \right)+s_1{Q}_2\left(\omega \right)$ and $\tilde{s}=s_1{Q}_1\left(\omega \right)-c_1{Q}_2\left(\omega \right)$ are given.

We also have $\left|Q\left(j\omega \right)\right|=\sqrt{{\tilde{c}}^2+{\tilde{s}}^2}$. Note that $D_{qr1}\left(\omega \right)=1+P_1\left(\omega \right){\displaystyle \sum _{k=0}^n{\tilde{c}}_k}+P_2\left(\omega \right){\displaystyle \sum _{k=0}^n{\tilde{s}}_k}$ and $D_{qr2}\left(\omega \right)=P_2\left(\omega \right){\displaystyle \sum _{k=0}^n{\tilde{c}}_k}-P_1\left(\omega \right){\displaystyle \sum _{k=0}^n{\tilde{s}}_k}$ are defined as the direct dynamic stiffness and the quadrature dynamic stiffness, respectively.

The loss factor $H_{qr}\left(\omega \right)$ of Figure 4 can be given as

$$H_{qr}=\frac{D_{qr2}\left(\omega \right)}{D_{qr1}\left(\omega \right)}=\frac{P_2\left(\omega \right){\displaystyle \sum _{k=0}^n{\tilde{c}}_k}-P_1\left(\omega \right){\displaystyle \sum _{k=0}^n{\tilde{s}}_k}}{1+P_1\left(\omega \right){\displaystyle \sum _{k=0}^n{\tilde{c}}_k}+P_2\left(\omega \right){\displaystyle \sum _{k=0}^n{\tilde{s}}_k}}$$

(26)

As $\omega =\frac{2k\pi }{T_r}=k{\omega }_r$, $k\in {\rm N}$, the loss factor $H_{qr}\left(k{\omega }_r\right)$ can be given as

$$H_{qr}\left(k{\omega }_r\right)=\frac{P_2\left(k{\omega }_r\right){\displaystyle \sum _{k=0}^n{\left[Q_1\left(k{\omega }_r\right)\right]}^k}-P_1\left(k{\omega }_r\right){\displaystyle \sum _{k=0}^n{\left[-Q_2\left(k{\omega }_r\right)\right]}^k}}{1+P_1\left(k{\omega }_r\right){\displaystyle \sum _{k=0}^n{\left[Q_1\left(k{\omega }_r\right)\right]}^k}+P_2\left(\omega \right){\displaystyle \sum _{k=0}^n{\left[-Q_2\left(k{\omega }_r\right)\right]}^k}}$$

(27)

A general multi-periodic disturbance of $d\left(t\right)={\displaystyle \sum _{i=1}^m{d}_i\left(t\right)}$ is given, where $d_i\left(t\right)$ is a periodic signal with period $T_{di}$, $\forall i=1,2,\dots ,m$. Therefore, we have a repetitive control system with $d\left(t\right)={\displaystyle \sum _{i=1}^m{d}_i\left(t\right)}$, as shown in Figure 5.

The repetitive controller is given as $R_{qm}\left(s\right)={\displaystyle \sum _{i=1}^m\frac{1}{1-Q{e}^{-s{T}_{di}}}}$, and its dynamic stiffness can be given as ${\widehat{D}}_{qr}\left(s\right)=1+P\left(s\right){\displaystyle \sum _{i=1}^m\frac{1}{1-Q\left(s\right)e^{-s{T}_{di}}}}$, which can be rewritten in the frequency domain, i.e.,

$${\widehat{D}}_{qr}\left(j\omega \right)=1+{\displaystyle \sum _{k=0}^n\left(P_1\left(\omega \right)+j{P}_2\left(\omega \right)\right){\left(\widehat{c}-j\widehat{s}\right)}^k}$$

(28)

where $\widehat{c}=Q_1\left(\omega \right){\displaystyle \sum _{i=1}^m\cos \left(\omega T_{di}\right)}+Q_2\left(\omega \right){\displaystyle \sum _{i=1}^m\sin \left(\omega T_{di}\right)}$ and $\widehat{s}=Q_1\left(\omega \right){\displaystyle \sum _{i=1}^m\sin \left(\omega T_{di}\right)}-Q_2\left(\omega \right){\displaystyle \sum _{i=1}^m\cos \left(\omega T_{di}\right)}$ are given.

Note that ${\widehat{D}}_{qr1}\left(\omega \right)=1+P_1\left(\omega \right){\displaystyle \sum _{k=0}^n{\widehat{c}}^k}+P_2\left(\omega \right){\displaystyle \sum _{k=0}^n{\widehat{s}}^k}$ and ${\widehat{D}}_{qr2}\left(\omega \right)=P_2\left(\omega \right){\displaystyle \sum _{k=0}^m{\widehat{c}}^k}-P_1\left(\omega \right){\displaystyle \sum _{k=0}^m{\widehat{s}}^k}$ are defined as the direct dynamic stiffness and the quadrature dynamic stiffness, respectively.

The loss factor ${\widehat{H}}_{qr}\left(\omega \right)$ of Figure 5 can be given as

$${\widehat{H}}_{qr}=\frac{{\widehat{D}}_{qr2}\left(\omega \right)}{{\widehat{D}}_{qr1}\left(\omega \right)}=\frac{P_2\left(\omega \right){\displaystyle \sum _{k=0}^n{\widehat{c}}^k}-P_1\left(\omega \right){\displaystyle \sum _{k=0}^n{\widehat{s}}^k}}{1+P_1\left(\omega \right){\displaystyle \sum _{k=0}^n{\widehat{c}}^k}+P_2\left(\omega \right){\displaystyle \sum _{k=0}^n{\widehat{s}}^k}}$$

(29)

4. Simulated Results Analysis

To analyze the dynamic stiffness of the quadcopter control system, an illustrated example of a quadcopter is shown in Figure 6, where its specifications are listed in

Table 1. Specification of the quadcopter.

Model weight	4100 g
Number of rotors	4
Frame size	690 mm
Power	22.8~26.0 V
Average current	15 A
Propeller	Diameter 14 inch 5.5 pitch
Max thrust	2100 g
Average thrust	1030 g
Hover thrust	1025 g
Idle thrust	50 g

. The controlled plant for the given quadcopter is measured by an analyzer, as shown in Figure 7, where the control plant of $G\left(s\right)=\frac{3.329s^2+340.19s+42051}{s^3+71.05s^2+5248.8s+56.31}$ can be calculated by the curve fitting. A gain of $K=0.95$ is pre-determined for the closed loop of $P\left(s\right)/\left(1+P\left(s\right)\right)$, which is stable, where $P\left(s\right)=KG\left(s\right)$ is given. In the designed case, the low-pass filter of $Q\left(s\right)=\frac{1}{\frac{s^2}{{\omega }_q{}^2}+\frac{2{\zeta }_{q}s}{{\omega }_q}+1}$, ${\omega }_q=20\mathrm{Hz}$, and ${\zeta }_q=\sqrt{2}/2$ can be designed to achieve the performance of the vibration rejection [24]. Under the proposed quadcopter flight operation, the vibration is measured as shown in Figure 8, where the harmonics of the measured results can be given as 3–11 Hz.

In this study, the flight controller mainly assists the controlled body, while the gimbal image stabilization controller provides the source of vibration interference during flight. The vibration of the gimbal provides a source of interference for the controlled object of the drone to perform experiments. As shown in Figure 6, the red arrow is the gimbal image stabilization controller, which simulates the environment when the system on the drone is abnormal by actively creating vibration in the roll axis. Figure 8 shows the measured FFT signal (x–y–z) of the IMU angle (Figure 8a) and angular acceleration (Figure 8b) to obtain the vibration signal of the proposed quadcopter flight operation. From the measured vibration harmonics in Figure 8, the harmonic frequencies of the vibration are determined as 1–20 Hz.

From the above analysis results, a disturbance with a period $T_d=2\mathrm{s}$ is given, and the dynamic stiffness of the proposed quadcopter control system is shown in Figure 8. From Figure 9a, it can be found that the magnitudes of the direct dynamic stiffness with the repetitive controller are larger than those without the repetitive controller, and the maximum of the direct dynamic stiffness can be found at harmonics of the periodic signal, i.e., $\frac{k}{T_d}=\frac{k}{2}=0.5k$, $k=1,2,3,\cdots$. From Figure 9b, the quadrature dynamic stiffness with $Q\left(s\right)$ is relatively smooth compared to that without $Q\left(s\right)$, which can be used to verify the stability being improved with $Q\left(s\right)$. Figure 10 shows the comparisons of loss factor, where the magnitudes of loss factor of the quadcopter control system with the repetitive controller are lower than those without the repetitive controller.

Figure 11 and Figure 12 shows the comparisons of the dynamic stiffness and loss factor with a periodic disturbance of a period $T_d=0.5\mathrm{s}$. Similar to the results of Figure 9 and Figure 10, it can be found that the dynamic stiffness with the repetitive controller is larger than that without the repetitive controller. The maximum of the direct dynamic stiffness can be found at harmonics of the periodic signal, i.e.,$\frac{k}{T_d}=\frac{k}{0.5}=2k$, $k=1,2,3,\dots$. In Figure 12, the loss factor of the quadcopter control system with the repetitive controller is lower than that without the repetitive controller.

Based on Figure 4 with the given parameters, a periodic disturbance is given as

$$[\begin{array}{l}d\left(t\right)=\sin \left(2\pi t\right)+\sin \left(4\pi t\right)\left[u\left(t\right)-u\left(t-5\right)\right]+\\ \hspace{1em}\hspace{1em}\hspace{1em}\sin \left(8\pi t\right)\left[u\left(t-3\right)-u\left(t-6\right)\right]+\sin \left(10\pi t\right)\left[u\left(t-6\right)-u\left(t-20\right)\right]\end{array}$$

(30)

Based on the disturbance input in (30), as shown Figure 13a, the error responses of the control systems in Figure 2 and Figure 4 can be obtained in Figure 13b. It can be found that the error of Figure 4 is smaller than that of Figure 2, and a repetitive error marked blue dashed line can be found in the control system of Figure 2. The rapid decay rate of Figure 4 can be found, where the larger error can be found at the commands’ switched time point, i.e., 0 s, 3 s, 5 s, and 6 s.

For the dynamic stiffness with the general periodic disturbance, considering the measured vibration harmonics in Figure 8, the harmonic frequencies are determined as 1–20 Hz, i.e., $T_d=0.05\mathrm{s}~1\mathrm{s}$. From Figure 14a, it can be found that the magnitudes of the direct dynamic stiffness with the repetitive controller are larger than those without the repetitive controller. From Figure 14b, the quadrature dynamic stiffness with $Q\left(s\right)$ is relatively smooth to compare that without $Q\left(s\right)$, which can be used to verify the stability being improved with $Q\left(s\right)$. Figure 15 shows the comparisons of loss factor, where the magnitudes of the loss factor of the quadcopter control system with the repetitive controller are lower than those without the repetitive controller. Based on Figure 5 with the given parameters, a multi-periodic disturbance is given in Figure 16a, i.e.,

$$d\left(t\right)=\sin \left(4\pi t\right)+2\sin \left(10\pi t\right)+0.5\sin \left(14\pi t\right)$$

(31)

Based on the disturbance input in (31), the error responses of the control systems in Figure 2 and Figure 5 can be obtained in Figure 16b. It can be found that the error of Figure 5 is smaller than that of Figure 2, and a repetitive error marked by the blue dashed line can be found in the control system of Figure 2. The rapid decay rate of Figure 5 can be found in Figure 16b.

5. Experimental Results Analysis

Figure 17 is the experimental setup of the study with the illustrated example of Figure 6. To provide a vibration situation under the quadcopter flight operation, a camera stabilizer constructed in the quadcopter is used to provide the vibration source. The control algorithm is implemented by the 32-bit ARM with 168 MHz installed on the quadcopter, and the inertial measurement unit is used to feedback the angular velocity and acceleration of X–Y–Z. Three kinds of circumstances are given to verify the proposed method. Note that the experiment environment wind speed is 0 to 1 knot in calm weather.

For the first case, the quadcopter control system only with feedforward + PID controller is given. The camera stabilizer is still, and any vibration is not provided by the stabilizer. Figure 18, Figure 19, Figure 20 and Figure 21 show the experimental results under the stabilizer still. It can be found that the acceleration oscillations of the x-axis and y-axis are smaller than $\pm 1\mathrm{m}/{\mathrm{s}}^2$. The angle of oscillation is less than $\pm 4$ degrees. A larger acceleration can be found at the time of taking off or landing of the quadcopter.

For the second case, the vibration is produced by the stabilizer, and the repetitive controller is given in the quadcopter control system. Figure 22, Figure 23, Figure 24 and Figure 25 show the experimental results. The acceleration oscillations of the x-axis and y-axis are smaller than $\pm 1\mathrm{m}/{\mathrm{s}}^2$. The angle of oscillation is less than $\pm 5$ degrees.

For the third case, only with feedforward control + PID controller, the vibration is provided. Figure 26, Figure 27, Figure 28 and Figure 29 show the experimental results. The acceleration oscillations of the x-axis and y-axis are larger than $\pm 3\mathrm{m}/{\mathrm{s}}^2$ and the angle of oscillation is less than $\pm 10$ degrees. At 14:51:00:000, however, the quadcopter system almost becomes out of control, even though the controller seems to stabilize the quadcopter system. At 14:55:00:000, the accelerations are larger than $\pm 40\mathrm{m}/{\mathrm{s}}^2$ and the angles are larger than $\pm 80$ degrees. This eventually leads to the crash of the experimental machine.

6. Conclusions

This paper proposed the dynamic stiffness analysis method of the quadcopter control system with the repetitive controller. A theoretical model of the dynamic stiffness of the quadcopter control system is established through analyzing the definition of the dynamic stiffness. The direct dynamic stiffness and the quadrature dynamic stiffness with variable frequencies are used to verify the effect of the quadcopter control system. Simulated and experimental results show that the magnitudes of the direct dynamic stiffness with the proposed repetitive controller are larger than those without the repetitive controller. The quadrature dynamic stiffness with the proposed repetitive controller is relatively smooth compared to that without the proposed repetitive controller, which can be used to verify the stability being improved by the designed repetitive controller. In addition, the magnitudes of the loss factor of the quadcopter control system with the repetitive controller are lower than those without the repetitive controller.