Compensation control for the quadrotor based on the active model in underground environment

Abstract

With its flexibility in operation, rapid inspection capability, omnidirectional detection, and strong terrain adaptability, the quadrotor has become an effective tool for underground inspection. The underground environment is characterized by narrow passages and turbulent airflow. When encountering strong gusts in underground spaces, the quadrotor maintains a large attitude angle to ensure stability. There is an significant model deviation between the actual dynamic model of the quadrotor and the hovering dynamic model. In general, the controller of the quadrotor is designed based on the hovering dynamic model in still air. The deviation between the model in the controller and the actual model can lead to the loss of control of the quadrotor. To overcome the problem above, at the environmental perception level, a real-time wind field prediction model is constructed based on the information obtained by wind speed sensors. By integrating this predicted information with an improved adaptive extended set membership filtering algorithm, a multi-step incremental prediction mechanism is introduced. At the control architecture level, a dynamic model compensation mechanism is proposed for nonlinear system characteristics. By identifying model bias online, an active model correction strategy is established. The incremental predictive control value is used as the feedforward compensation input, and an optimal control problem with constraints is established. Finally, a simulation example is presented to verify the effectiveness of the quadrotor control method.

Introduction

The working environment in coal mines is extremely harsh, ¹ and there are great risks for inspection personnel using traditional manual inspection methods. The use of robots for inspections is becoming a mainstream trend in underground environments. At present, coal miners need to carry professional equipment to conduct inspections in underground tunnels every day, in order to identify potential hazards promptly. This labor-intensive job exposes workers to constant high risks.

The inspection of coal mine tunnels mainly relies on ground robots, which include tracked and wheeled configurations. However, ground robots have limited mobility on complex terrains, and often have blind spots in the field of view. It is difficult for these configurations to meet requirements in complex and diverse environments. With the rapid development of unmanned aerial vehicle(UAV) technology, applications of quadrotors have become a research focus. Quadrotors are widely used in industrial inspection, aerial photography, agricultural operations and other fields.^2,3 Benefiting from flexible operation, rapid inspection, omnidirectional detection, and strong terrain adaptability, quadrotors enable inspections in complex underground environments. Thus, quadrotors are poised to become a significant robot configuration for underground inspections.

The quadrotor is classified as an underactuated system since it possesses only four control inputs yet generates six motion outputs. Input signals consist of the desired rotation rates of the four motors, while the outputs encompass three translational coordinates and three rotational coordinates relative to the quadrotor’s body axes. The quadrotor also exhibits nonlinear coupling dynamics. Three factors contribute to this. First, attitude representation is nonlinear. It evolves on SO(3)—the mathematical group of all three-dimensional rotations—which means Euler angles suffer from trigonometric singularities. Second, translation and attitude are coupled: thrust must be rotated from the body frame to the inertial frame using a rotation matrix, and horizontal motion requires tilting the vehicle to direct thrust. Third, gyroscopic effects arise from the cross product of angular velocity and angular momentum, introducing quadratic nonlinearities. Finally, the quadratic mapping from rotor speeds to thrust and moments compounds these effects.

For quadrotors, existing control methods can mainly be divided into three categories: linear control, nonlinear control, and intelligent control methods. Li and Li ⁴ designed a Proportional-Integral-Derivative (PID) controller for three-dimensional position and attitude control of quadrotors. Simulation results show that the quadrotor can stably fly under wind disturbance, but with slight overshoot. Pan et al. ⁵ similarly developed an optimal PID controller with a Kalman filter for target trajectory estimation. Their results indicate that the system can quickly reach the desired height without steady-state error, yet exhibits significant overshoot.

To overcome the above overshoot problem, Okyere et al. ⁶ propose a step-by-step LQR controller to control the altitude motion of quadrotors. Their numerical experiments show that if the R value is high and the K value is low, the response speed will be slower. When the K value is high, the state response reaches zero faster, but the controller response has large steady-state error. Shah et al. ⁷ study the trajectory tracking control of quadrotors using gain scheduling integral LQR, where an integral term is added to improve the tracking performance of the quadrotor. Their results indicate that the continuous gain control law can solve the problems of jitter and discontinuity.

In order to control the attitude response of the quadrotor, Han et al. ⁸ develop an H_∞ controller based on the gray-box method, and conduct related experiments. ⁹ Rich et al. ¹⁰ study a robust H_∞controller for quasi-hover conditions using the Glover-McFarlane loop shaping method. Bonna and Camino^11–13 study linearization of the feedback for achieving attitude control and stabilizing the quadrotor. Lu et al. ¹⁴ propose a backstepping controller, which uses online optimization techniques to control the position and attitude of quadrotors. Their strategy can simultaneously solve path planning and tracking control problems, as traditional control designs only address these problems separately. First, accessible trajectories are generated by the planning generator, and then the tracking controller is adjusted and optimized to achieve the new planned path. The experiment shows that the accuracy of the quadrotor trajectory is within an acceptable range. Dydek et al. ¹⁵ investigate an adaptive controller for the quadrotor with the baseline trajectory tracking controller to enhance system performance. The effectiveness of this approach is validated through both simulations and indoor experiments. Their adaptive controller demonstrates strong robustness against parameter uncertainties. The controller is particularly effective in mitigating the impact of abnormal loss of thrust caused by component failures or physical damage.

Shao et al. ¹⁶ investigate a collision-free elliptical target encirclement cooperative control method with improved convergence characteristics for UAVs under wind disturbances. First, an enhanced fixed-time encirclement guidance law is designed, enabling UAVs to reach user-preset elliptical trajectories within explicit time constraints independent of initial state conditions. This guidance law achieves rapid error decay through a carefully designed radial error feedback mechanism (with power terms dynamically switching between fractional exponents and polynomial terms). Specifically, they ¹⁷ employ a simple filtering-based unknown system dynamics estimator to precisely compensate for uncertainties such as wind disturbances in a computationally efficient manner.

Cömert et al. ¹⁸ and Castillo Zamora et al. ¹⁹ study the sliding mode controller (SMC) for the position and altitude control of quadrotors. The results indicate that the settling time of the SMC controller is shorter, but it can generate large pitch and roll angles exceeding 20°, which is unfavorable in practice. As an alternative approach, the model predictive controller (MPC) is used to control processes while satisfying a series of constraints. This controller relies on a linearized dynamic model obtained through the system identification process. Ganga and Dharmana ²⁰ conduct a comparative study between traditional PID controllers and MPC for the altitude motion control of the quadrotor. Simulation results show that compared with PID, MPC has better control characteristics, which has shorter settling time and no overshoot.

Intelligent controllers are mainly designed based on bionic principles, such as fuzzy theory or neural network theory. Raza and Gueaieb ²¹ propose a fuzzy logic controller to control the position and attitude angle of the quadrotor under disturbance conditions. The results show that the fuzzy logic controller can provide good performance even under sensor noise and wind disturbance. Zare et al. ²² study the trajectory tracking control of the quadrotor using a fuzzy logic controller. The simulation results show that the fuzzy logic controller performs well in tracking the desired position and attitude angle of the quadrotor.

The sizes of mine tunnels are only 2–5 m in width and 3–6 m in height. Due to narrow spaces, abrupt airflow changes in turning areas, and ground effect interference, the quadrotor inspection in underground environment still remains a challenge. For instance, in tunnels, drones must maintain a minimum safe distance of 0.5 m from coal walls to prevent collisions and stay at least 0.5 m above the ground to mitigate ground effect disturbances. There is no natural wind in complex underground environment, and fans are used for air exchange usually. The fans usually have a high power and produce a large local airflow. The underground mining area has a complex network containing various T-shaped intersections and crossroads. Due to the narrowness of the channels, strong wind disturbances consisting of steady winds and gusts readily form. Given the quadrotor’s underactuated nature, multivariable coupling, and nonlinear dynamic characteristics, its control stability may significantly decrease under strong external disturbances during underground inspection.

The existing control schemes (such as traditional PID algorithm and sliding mode control algorithm) have obvious shortcomings in dealing with sudden underground disturbances. On the one hand, insufficient state prediction accuracy leads to control command oscillations and trajectory deviations; On the other hand, the quadrotor maintains a large attitude angle to ensure stability when encountering strong gusts in underground spaces. There is a significant model deviation between the actual dynamic model of the quadrotor and the hovering dynamic model. The controller of the quadrotor is usually designed based on the hovering dynamic model in still air. The deviation between the model in the controller and the actual model can lead to loss of control. In particular, the trajectory deviation caused by wind disturbances becomes the main threat to quadrotor flight safety.

To overcome the problem above, this study proposes corresponding improvement methods. At the environmental perception level, a real-time wind field prediction model is constructed through wind speed sensor data and proper orthogonal decomposition (POD) method.^23–25 Combined with an improved adaptive extended set membership filtering (ASMF) algorithm, a multi-step incremental prediction mechanism is introduced. This mechanism effectively compensates for model bias caused by flight mode switching and operating point migration by predicting the optimal control effect over multiple future control cycles, while an online update strategy is adopted to reduce computational complexity. At the control architecture level, considering online identification of model bias, an active model correction strategy is established, and an optimal control problem with constraints is constructed. By solving the above optimal control problem for the compensation amount in real-time, effective control of the quadrotor in underground environments is achieved.

To address challenges posed by complex underground environments to quadrotor inspection, this study proposes a series of methods. The contributions are summarized as follows. (1) Conventional control approaches typically rely on fixed mathematical models, rendering them prone to failure when dynamic models severely mismatch actual quadrotor states due to turbulent airflow and confined spaces in underground environments. In contrast, the proposed method dynamically identifies model deviations by real-time monitoring of discrepancies between system states and predicted values, subsequently adjusting control parameters online for rapid adaptation. This eliminates dependence on offline recalibration, significantly enhancing control robustness. (2) Traditional methods lack active compensation mechanisms for disturbance rejection. Upon encountering sudden disturbances such as strong airflow, they experience feedback delays and rely on fixed parameters, making timely attitude adjustment impossible. The proposed method employs multi-step incremental prediction results as feedforward signals to preemptively compensate for future disturbance effects, overcoming inherent latency limitations of conventional feedback control. Such an approach substantially improves dynamic response and reduces control deviations. As quadrotors are strongly coupled nonlinear systems, traditional linear methods struggle to handle complex coupling between attitude and position, particularly under model mismatch. In this study, a constrained optimal control problem is formulated that comprehensively considers actuator limits and safety boundaries, ensuring effective and safe control. This avoids control commands exceeding actuator limits due to neglected constraints in conventional methods. (3) Furthermore, whereas traditional prediction methods exhibit weak environmental perception, inability to perceive dynamic changes in real-time, and low multi-step prediction accuracy with error accumulation, we propose a method to compute compensation values efficiently in real-time. This guarantees that the nominal controller maintains high-performance control under model mismatch conditions, sustaining its precise control throughout operation. Consequently, it satisfies stringent demands of underground inspection for long-term trajectory planning and precise control, thereby enabling precise inspection in complex underground environments.

The remainder of this paper is organized as follows. In Section 2, a real-time wind field prediction model is constructed through wind sensor information and the POD method. In Section 3, the dynamic modeling of the quadrotor is introduced. In Section 4, an optimal control state prediction algorithm is proposed, integrating predicted wind data with a multi-step incremental prediction mechanism. In Section 5, an active model correction strategy is introduced to adaptively refine the dynamic model based on real-time deviations, followed by the formulation of the optimal control problem. In Section 6, numerical simulations are given to validate the effectiveness and robustness of the proposed methods. Finally, the conclusions of this study are presented in Section 7.

Dynamic modeling of the quadrotor

We begin by defining the coordinate frames. The local inertial frame $O_g{X}_g{Y}_g{Z}_g$ is established following the right-hand rule. The origin $O_g$ is fixed at the takeoff point of the quadrotor. The axes are defined as follows: $O_g{X}_g$ points toward geographic north, $O_g{Y}_g$ points east, and $O_g{Z}_g$ points downward toward the Earth’s center. The body frame is usually represented as $O_b{X}_b{Y}_b{Z}_b$ , the origin $O_b$ is located at the center of mass of the quadrotor body. $O_b{X}_b$ points forward along the direction of flight, $O_b{Y}_b$ points to the right when viewed from the cockpit looking forward, and the $O_b{Z}_b$ points downward. For convenience, parameter descriptions are provided in Appendix Table 1.

$\varphi ,\theta ,\psi$ represent the roll angle, pitch angle, and yaw angle of the quadrotor, respectively. Roll angle $\varphi$ refers to rotation about $O_b{X}_b$ ; the pitch angle $\theta$ represents rotation about $O_b{Y}_b$ ; the pitch angle $\psi$ represents rotation about $O_b{Z}_b$ . Following the right-hand rule, positive rotations are defined as clockwise when looking along the positive axis direction.

When the axes $O_b{X}_b,O_b{Y}_b,$ and $O_b{Z}_b$ independently rotate, the corresponding transformation matrices from the local inertial frame to the body frame are given by:

$${\mathit{R}}_g^b(x,\varphi )=\left[\begin{array}{ccc}1\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& \cos \varphi \hfill & -\sin \varphi \hfill \\ \multicolumn{1}{c}{0}& \sin \varphi \hfill & \cos \varphi \hfill \end{array}\right]$$

(1)

$${\mathit{R}}_g^b(y,\theta )=\left[\begin{array}{ccc}\cos \theta \hfill & 0\hfill & \sin \theta \hfill \\ \multicolumn{1}{c}{0}& 1\hfill & 0\hfill \\ \multicolumn{1}{c}{-\sin \theta }& 0\hfill & \cos \theta \hfill \end{array}\right]$$

(2)

$${\mathit{R}}_g^b(z,\psi )=\left[\begin{array}{ccc}\cos \psi \hfill & -\sin \psi \hfill & 0\hfill \\ \multicolumn{1}{c}{\sin \psi }& \cos \psi \hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 1\hfill \end{array}\right]$$

(3)

The transformation matrix from the local inertial frame to the body frame for three -axis rotation can be expressed as

$$\begin{array}{c}\multicolumn{1}{c}{{\mathbf{R}}_g^b={\mathbf{R}}_g^b(z,\psi ){\mathbf{R}}_g^b(y,\theta ){\mathbf{R}}_g^b(x,\varphi )=\left[\begin{array}{ccc}\cos \theta \cos \psi \hfill & \cos \theta \sin \psi \hfill & -\sin \theta \hfill \\ \multicolumn{1}{c}{\sin \varphi \sin \theta \sin \psi -\cos \varphi \sin \psi }& \sin \varphi \sin \theta \sin \psi +\cos \varphi \cos \psi \hfill & \sin \varphi \cos \theta \hfill \\ \multicolumn{1}{c}{\cos \varphi \sin \theta \cos \psi +\sin \varphi \sin \psi }& \cos \varphi \sin \theta \sin \psi -\sin \varphi \cos \psi \hfill & \cos \varphi \cos \theta \hfill \end{array}\right]}\end{array}$$

(4)

The above matrix represents the relationship that quadrotor first rotates $\varphi$ about $O_b{X}_b$ , then rotates $\theta$ about $O_b{Y}_b$ , finally rotates $\psi$ about $O_g{Z}_g$ . The transformation matrix from the body frame to the local inertial frame may be given by

$$\begin{array}{c}\multicolumn{1}{c}{{\mathit{R}}_b^g=({\mathit{R}}_g^b{)}^{-1}=\left[\begin{array}{ccc}\cos \theta \cos \psi \hfill & \sin \varphi \sin \theta \sin \psi -\cos \varphi \sin \psi \hfill & \cos \varphi \sin \theta \cos \psi +\sin \varphi \sin \psi \hfill \\ \multicolumn{1}{c}{\cos \theta \sin \psi }& \sin \varphi \sin \theta \sin \psi +\cos \varphi \cos \psi \hfill & \cos \varphi \sin \theta \sin \psi -\sin \varphi \cos \psi \hfill \\ \multicolumn{1}{c}{-\sin \theta }& \sin \varphi \cos \theta \hfill & \cos \varphi \cos \theta \hfill \end{array}\right]}\end{array}$$

(5)

During the flight, the quadrotor is mainly subjected to the following forces and moments: gravity, lifts generated by rotor rotation, moments caused by these lifts, the torques generated by rotor rotation, and additional force brought by the air. As shown in Figure 1, the wheelbase between adjacent rotors is $l_a$ . And $T_i,i=1,2,3,4$ represents the lift generated by each rotor, $M_{\mathrm{ri}},i=1,2,3,4$ represents the counter torque generated by each of them. In addition, $m$ represents the overall mass of the quadrotor, $\mathit{f}(t)$ represents the local gravitational acceleration. The lift forces generated by the rotors in the body frame of the quadrotor are

$$\{\begin{array}{c}\multicolumn{1}{c}{F_{\mathrm{Tbxi}}=0}\\ \multicolumn{1}{c}{F_{\mathrm{Tbyi}}=0}\\ \multicolumn{1}{c}{F_{\mathrm{Tbzi}}=-T_i}\end{array},i=1,2,3,4$$

(6)

Where $F_{\mathrm{Tbxi}},F_{\mathrm{Tbyi}},F_{\mathrm{Tbzi}},i=1,2,3,4$ is the component of lift of each rotor. $T_i,i=1,2,3,4$ is lift generated by each rotor, and the magnitude of the lift is proportional to the rotor lift coefficient and the square of the angular velocity of each rotor:

$$T_i=C_T{\omega }_i^2,i=1,2,3,4$$

(7)

where $C_T$ is the lift coefficient, ${\omega }_i,i=1,2,3,4$ denotes the rotational speed of each rotor. Thus, the moment in different directions $M_{\mathrm{Tbxi}},M_{\mathrm{Tbyi}},M_{\mathrm{Tbzi}}$ , $i=1,2,3,4$ generated by the lift can be determined by

$$\{\begin{array}{c}\multicolumn{1}{c}{M_{\mathrm{Tbxi}}=l_a\frac{\sqrt{2}}{2}{T}_i\cos (\frac{i\pi }{2}+\frac{\pi }{4})}\\ \multicolumn{1}{c}{M_{\mathrm{Tbyi}}=l_a\frac{\sqrt{2}}{2}{T}_i\sin (\frac{i\pi }{2}+\frac{\pi }{4})}\\ \multicolumn{1}{c}{M_{\mathrm{Tbzi}}=0}\end{array},i=1,2,3,4$$

(8)

Figure 1.

Force analysis of quadrotor.

The counter torque of each rotor can be determined by the following equation

$$M_{\mathrm{ri}}=C_r{\omega }_i^2,i=1,2,3,4$$

(9)

where $C_r$ is the counter torque coefficient. The component of each rotor alone different axes can be obtained by

$$\{\begin{array}{c}\multicolumn{1}{c}{M_{\mathrm{rbxi}}=0}\\ \multicolumn{1}{c}{M_{\mathrm{rbyi}}=0}\\ \multicolumn{1}{c}{M_{\mathrm{rbzi}}=M_{\mathrm{ri}}\tan (\frac{i\pi }{2}-\frac{\pi }{4})}\end{array},i=1,2,3,4$$

(10)

The three-axis components of the body’s gravity in the body frame can be transformed through coordinate transformation relationships as follows

$$G_b=R_g^b\left[\begin{array}{c}0\hfill \\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{\mathrm{mg}}\end{array}\right]=\{\begin{array}{c}{G}_{\mathrm{bx}}=-\mathrm{mg}\sin \theta \hfill \\ \multicolumn{1}{c}{G_{\mathrm{by}}=\mathrm{mg}\cos \theta \sin \varphi }\\ \multicolumn{1}{c}{G_{\mathrm{bz}}=\mathrm{mg}\cos \theta \cos \varphi }\end{array}$$

(11)

Thus, $T_{\mathrm{bx}},T_{\mathrm{by}},T_{\mathrm{bz}}$ , the forces in various directions of the quadrotor, in the body frame can be obtained by

$$\begin{array}{c}{T}_{\mathrm{bx}}=\sum _{i=1}^4{F}_{\mathrm{Tbxi}}+G_{\mathrm{bx}}+D_{\mathrm{bx}}\hfill \\ T_{\mathrm{by}}=\sum _{i=1}^4{F}_{\mathrm{Tbyi}}+G_{\mathrm{by}}+D_{\mathrm{by}}\hfill \\ T_{\mathrm{bz}}=\sum _{i=1}^4{F}_{\mathrm{Tbzi}}+G_{\mathrm{bz}}+D_{\mathrm{bz}}\hfill \end{array}$$

(12)

where $D_{\mathrm{bx}},D_{\mathrm{by}},D_{\mathrm{bz}}$ is the additional force generated by the external airflow,which can be calculated using the following equation:

$$\begin{array}{c}{D}_{\mathrm{bx}}=C_{\mathrm{bx}}(\alpha )\left|v_{\mathrm{rex}}\right|v_{\mathrm{rex}}{S}_{\mathrm{refx}}\hfill \\ D_{\mathrm{by}}=C_{\mathrm{by}}(\beta )\left|v_{\mathrm{rey}}\right|v_{\mathrm{rey}}{S}_{\mathrm{refx}}\hfill \\ D_{\mathrm{bz}}=C_{\mathrm{bz}}(\alpha ,\beta )\left|v_{\mathrm{rey}}\right|v_{\mathrm{rey}}{S}_{\mathrm{refz}}\hfill \end{array}$$

(13)

In the above equation, $C_{\mathrm{bx}}(\alpha )$ is the axial aerodynamic coefficient alone the $x$ axis, which is a function of the axial angle of attack $\alpha$ . $C_{\mathrm{by}}(\beta )$ is the axial aerodynamic coefficient alone the $y$ axis, which is a function of the axial angle of attack $\beta$ . The aerodynamic coefficient $C_{\mathrm{bz}}(\alpha ,\beta )$ is a function of both $\alpha$ and $\mathit{f}(t)$ . The parameters $S_{\mathrm{refx}}$ , $S_{\mathrm{refy}}$ , and $S_{\mathrm{refz}}$ correspond to the reference areas of the quadrotor in the x, y, and z axes, respectively. ${\mathit{v}}_{\mathrm{re}}=\mathit{V}-\mathit{v}=[v_{\mathrm{rex}}{v}_{\mathrm{rey}}{v}_{\mathrm{rez}}]$ is the relative velocity between the external flow and the quadrotor body. $\mathit{V}=[V_x{V}_y{V}_z]$ the airflow velocity. $\mathit{v}={\left[v_x{v}_y{v}_z\right]}^T$ is the velocity vector of the quadrotor. Accordingly, the torque of the quadrotor in various directions can be obtained in the body frame as follows

$$\begin{array}{c}{M}_{\mathrm{bx}}=\sum _{i=1}^4{M}_{\mathrm{Tbxi}}+\sum _{i=1}^4{M}_{\mathrm{rbxi}}\hfill \\ M_{\mathrm{by}}=\sum _{i=1}^4{M}_{\mathrm{Tbyi}}+\sum _{i=1}^4{M}_{\mathrm{rbyi}}\hfill \\ M_{\mathrm{bz}}=\sum _{i=1}^4{M}_{\mathrm{Tbzi}}+\sum _{i=1}^4{M}_{\mathrm{rbzi}}\hfill \end{array}$$

(14)

The overall external force vector, the velocity vector and angular velocity vector of the body frame are defined by ${\mathit{T}}_b=[T_{\mathrm{bx}},T_{\mathrm{by}},T_{\mathrm{bz}}{]}^T$ , $\mathit{v}={\left[\begin{array}{ccc}{v}_x\hfill & v_y\hfill & v_z\hfill \end{array}\right]}^T$ , and $\mathbf{\omega }={\left[\begin{array}{ccc}p\hfill & q\hfill & r\hfill \end{array}\right]}^T$ , respectively. According to Newton’s second law, we have

$$m\frac{d\mathit{v}}{\mathrm{dt}}={\mathit{T}}_b$$

(15)

The derivative $\frac{d\mathit{v}}{\mathrm{dt}}$ can be expressed as

$$\frac{d\mathit{v}}{\mathrm{dt}}=\stackrel{·}{\mathit{v}}+\mathbf{\omega }\times \mathit{v}$$

(16)

Substituting equation (16) into equation (15) yields

$$\stackrel{·}{\mathit{v}}=\frac{{\mathit{T}}_b}{m}-\mathbf{\omega }\times \mathit{v}=\{\begin{array}{c}\multicolumn{1}{c}{\stackrel{·}{v_x}=\frac{T_{\mathrm{bx}}}{m}+r{v}_y-q{v}_z}\\ \multicolumn{1}{c}{\stackrel{·}{v_y}=\frac{T_{\mathrm{by}}}{m}+p{v}_z-r{v}_x}\\ \multicolumn{1}{c}{\stackrel{·}{v_z}=\frac{T_{\mathrm{bz}}}{m}+q{v}_x-p{v}_y}\end{array},$$

(17)

where $p,q$ and $r$ represent roll, pitch, and yaw velocities, respectively. According to the conversion relationship between the local inertial frame and the body frame, the following conversion relationship can be obtained for the velocity vector and acceleration vector:

$$\left[\begin{array}{c}\stackrel{··}{x}\hfill \\ \multicolumn{1}{c}{\stackrel{··}{y}}\\ \multicolumn{1}{c}{\stackrel{··}{z}}\end{array}\right]={\mathit{R}}_b^g\stackrel{·}{\mathit{v}}={\mathit{R}}_b^g\left[\begin{array}{c}\frac{T_{\mathrm{bx}}}{m}+r{v}_y-q{v}_z\hfill \\ \multicolumn{1}{c}{\frac{T_{\mathrm{by}}}{m}+p{v}_z-r{v}_x}\\ \multicolumn{1}{c}{\frac{T_{\mathrm{bz}}}{m}+p{v}_z-r{v}_x}\end{array}\right]$$

(18)

According to the law of angular momentum, it can be concluded that

$$\frac{d(\mathit{J}\omega )}{\mathrm{dt}}={\mathit{M}}_b,$$

(19)

where $M_b=[M_{\mathrm{bx}},M_{\mathrm{by}},M_{\mathrm{bz}}]$ is the overall external force vector, $J=\left[\begin{array}{ccc}{J}_x\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& J_y\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & J_z\hfill \end{array}\right]$ , $J_x,J_y,J_z$ are the inertias about the axes. According to the derivative relationship, we have:

$$\frac{d(\mathit{J}\omega )}{\mathrm{dt}}=\mathit{J}\stackrel{·}{\mathbf{\omega }}+\mathit{\omega }J\omega$$

(20)

By substituting equation (20) into equation (19), we can obtain

$$\mathit{J}\stackrel{·}{\mathbf{\omega }}={\mathit{M}}_b-\mathit{\omega }J\omega$$

(21)

Expanding the above equation yields the rotational dynamics in component form:

$$\begin{array}{c}\stackrel{·}{p}=\frac{M_{\mathrm{bx}}}{J_x}+\frac{(J_y-J_z)}{J_x}\mathrm{qr}\hfill \\ \stackrel{·}{q}=\frac{M_{\mathrm{by}}}{J_y}+\frac{(J_z-J_x)}{J_y}\mathrm{pr}\hfill \\ \stackrel{·}{r}=\frac{M_{\mathrm{bz}}}{J_z}+\frac{(J_x-J_y)}{J_z}\mathrm{pq}\hfill \end{array}$$

(22)

The relationship between the angular velocity of the body and the Euler angular velocity can be expressed as:

$$\begin{array}{c}\multicolumn{1}{c}{\left[\begin{array}{c}p\hfill \\ \multicolumn{1}{c}{q}\\ \multicolumn{1}{c}{r}\end{array}\right]={[\mathit{R}(y,\theta )\mathit{R}(x,\varphi )]}^{-1}\stackrel{·}{\psi }{\mathit{e}}_3+{[\mathit{R}(x,\theta )]}^{-1}\stackrel{·}{\theta }{\mathit{e}}_2+\stackrel{·}{\varphi }{\mathit{e}}_1}\\ =\left[\begin{array}{ccc}1\hfill & 0\hfill & -\sin \theta \hfill \\ \multicolumn{1}{c}{0}& \cos \varphi \hfill & \cos \theta \sin \varphi \hfill \\ \multicolumn{1}{c}{0}& -\sin \varphi \hfill & \cos \theta \cos \varphi \hfill \end{array}\right]\left[\begin{array}{c}\stackrel{·}{\varphi }\hfill \\ \multicolumn{1}{c}{\stackrel{·}{\theta }}\\ \multicolumn{1}{c}{\stackrel{·}{\psi }}\end{array}\right],\hfill \end{array}$$

(23)

where $e_1={\left[\begin{array}{ccc}1\hfill & 0\hfill & 0\hfill \end{array}\right]}^T$ , $e_2={\left[\begin{array}{ccc}0\hfill & 1\hfill & 0\hfill \end{array}\right]}^T$ , and $e_3={\left[\begin{array}{ccc}0\hfill & 0\hfill & 1\hfill \end{array}\right]}^T$ . Assuming that

$${\mathit{S}}_g=\left[\begin{array}{ccc}1\hfill & 0\hfill & -\sin \theta \hfill \\ \multicolumn{1}{c}{0}& \cos \varphi \hfill & \cos \theta \sin \varphi \hfill \\ \multicolumn{1}{c}{0}& -\sin \varphi \hfill & \cos \theta \cos \varphi \hfill \end{array}\right],$$

(24)

we have

$$\left[\begin{array}{c}p\hfill \\ \multicolumn{1}{c}{q}\\ \multicolumn{1}{c}{r}\end{array}\right]={\mathit{S}}_g^b\left[\begin{array}{c}\stackrel{·}{\varphi }\hfill \\ \multicolumn{1}{c}{\stackrel{·}{\theta }}\\ \multicolumn{1}{c}{\stackrel{·}{\psi }}\end{array}\right]$$

(25)

The angular velocity conversion matrix from the body frame to the local inertial frame is as follows:

$${\mathit{S}}_b^g=({\mathit{S}}_g^b{)}^{-1}=\left[\begin{array}{ccc}1\hfill & \sin \varphi \tan \theta \hfill & \cos \varphi \tan \theta \hfill \\ \multicolumn{1}{c}{0}& \cos \varphi \hfill & -\sin \varphi \hfill \\ \multicolumn{1}{c}{0}& \sin \varphi /\cos \theta \hfill & \cos \varphi /\cos \theta \hfill \end{array}\right]$$

(26)

Combining equation (22) and equation (25), the dynamic equation in the local inertial frame may be given by

$$\left[\begin{array}{c}\stackrel{··}{\varphi }\hfill \\ \multicolumn{1}{c}{\stackrel{··}{\theta }}\\ \multicolumn{1}{c}{\stackrel{··}{\psi }}\end{array}\right]={\mathit{S}}_b^g\left[\begin{array}{c}\multicolumn{1}{c}{\stackrel{·}{p}}\\ \multicolumn{1}{c}{\stackrel{·}{q}}\\ \multicolumn{1}{c}{\stackrel{·}{r}}\end{array}\right]={\mathit{S}}_b^g\left[\begin{array}{c}\multicolumn{1}{c}{\frac{M_{\mathrm{bx}}}{J_x}+\frac{(J_y-J_z)}{J_x}\mathrm{qr}}\\ \multicolumn{1}{c}{\frac{M_{\mathrm{by}}}{J_y}+\frac{(J_z-J_x)}{J_y}\mathrm{pr}}\\ \multicolumn{1}{c}{\frac{M_{\mathrm{bz}}}{J_z}+\frac{(J_x-J_y)}{J_z}\mathrm{pq}}\end{array}\right]$$

(27)

Gust load prediction based on observation data and POD method

The underground environment is complex and time-changing, and it is easy to generate random local gusts of wind. In underground mines, gusts of wind are usually caused by the acceleration of airflow in some complex tunnel terrains. The biggest uncertainty of the quadrotor in underground inspections is the loss of control caused by gusts. Therefore, the real-time capability, robustness, and accuracy of gusts prediction are extremely high.

To meet these requirements, we introduce a gust loads prediction method based on the POD method.^23–25 First, wind sensors at different locations in underground tunnels and on quadrotors are installed to collect information on environmental wind, and then the intensity and direction of gusts based on the collected data are predicted. $\mathit{A}$ can be defined based on the snapshots $V_j^i(\mathit{x}),1\le j\le M,1\le i\le L$ in the following matrix form

$$\mathit{A}=\left(\begin{array}{cccc}{V}_1^1(\mathit{x})\hfill & V_1^2(\mathit{x})\hfill & \cdots \hfill & V_1^L(\mathit{x})\hfill \\ \multicolumn{1}{c}{V_2^1(\mathit{x})}& V_2^2(\mathit{x})\hfill & \cdots \hfill & V_2^L(\mathit{x})\hfill \\ \multicolumn{1}{c}{\cdots }& \cdots \hfill & \cdots \hfill & \cdots \hfill \\ \multicolumn{1}{c}{V_M^1(\mathit{x})}& V_M^2(\mathit{x})\hfill & \cdots \hfill & V_M^L(\mathit{x})\hfill \end{array}\right)$$

(28)

where M denotes the number of grid nodes, and L represents the number of snapshots.

The mean value of these instantaneous images ${\overline{V}}_j(\mathit{x})$ can be given by the following formulation

$${\overline{V}}_j(\mathit{x})=\sum _{i=1}^L{V}_j^i(\mathit{x})/L,1\le j\le M$$

(29)

Using the difference between each sample and the average value of the sample, we obtain a new set of instantaneous images as follows

$${\hat{V}}_j^i(\mathit{x})={\overline{V}}_j(\mathit{x})-{\overline{V}}_j^i(\mathit{x}),1\le j\le M,1\le i\le L$$

(30)

Let $\mathit{B}=({\hat{V}}_j^i(\mathit{x}){)}_{i=1}^L$ be a matrix. The inner product $<\mathit{a},\mathit{b}>={\mathit{a}}^T\mathit{b}$ is defined in $L^2$ space; $\mathit{a},\mathit{b}$ is the column vector defined in measure space $\Omega$ . Assume that the velocity $V(\mathit{x},t)$ is approximated as a linear combination of known snapshots.

$$V(\mathit{x},t)\approx {\overline{V}}_j(\mathit{x})+\sum _{j=1}^L{a}_j(t){\varphi }_j(\mathit{x}),$$

(31)

where ${\phi }_j(\mathit{x}),1,2,3,L$ forms an orthonormal basis. Following the reference,^23–25 the eigenvectors ${\phi }_j(\mathit{x}),1,2,3,L$ are computed via the generalized eigenvalue problem of $\mathit{C}=\mathit{B}{\mathit{B}}^T$ . In order to solve the above eigenvalue problem, we use the singular value decomposition (SVD) method to decompose the matrix into the following forms:

$$\mathit{B}=\mathit{U}\left(\begin{array}{cc}S\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill \end{array}\right){\mathit{W}}^T,$$

(32)

where, $\mathit{U}=({\varphi }_1,{\varphi }_2,\cdots .{\varphi }_M)\in R^{M\times M}$ , $\mathit{W}=({\phi }_1,{\phi }_2,\cdots .{\phi }_L)\in R^{L\times L}$ are orthogonal eigenvectors on ${\mathit{B}}^T\mathit{B}$ and $\mathit{B}{\mathit{B}}^T$ , respectively. $S=\mathrm{diag}({\sigma }_1,{\sigma }_2,..{\sigma }_l)\in R^{l\times l}$ is a diagonal matrix. ${\sigma }_i(i=1,2..l)$ is a positive value. As the number of discrete grid points is generally much larger than the number of instantaneous images ( $M>>L$ ), yet their nonzero eigenvalues are identical. Thus, we can first solve the eigenvalue problem of the matrix ${\mathit{B}}^T\mathit{B}$ to obtain the eigenvector with non-zero eigenvalue as follows

$${\mathit{B}}^T\mathit{B}{\phi }_j(\mathit{x})={\lambda }_j{\phi }_j(\mathit{x}),1\le j\le l$$

(33)

Since the eigenvalues of the matrix $\mathit{B}{\mathit{B}}^T$ and ${\mathit{B}}^T\mathit{B}$ have the above relationship, we can see that ${\lambda }_j={\sigma }_j^2,1\le j\le l$ . Therefore, we can get the corresponding eigenvector about the non-zero eigenvalue of the matrix $\mathit{B}{\mathit{B}}^T$ , where $\mathit{U}=({\varphi }_1,{\varphi }_2,\cdots .{\varphi }_M)$ are the orthogonal basis of the required solution, in which

$${\varphi }_j(x)=\frac{1}{{\sigma }_j}\mathit{B}{\varphi }_j,1\le j\le l.$$

(34)

After the POD method is used to obtain the orthogonal basis, it is necessary to select sufficiently large truncation order of basis vectors for approximating $V(\mathit{x},t)$ , ${\lambda }_i$ represents the energy contained in the $\mathit{f}(t)$ th mode. Generally, the following truncation is used to select the orthogonal basis.

$$E(N)=\sum _{i=1}^N{\lambda }_i/\sum _{i=1}^l{\lambda }_i\ge E$$

(35)

where $E$ is usually set to 0.95. When $E(N)\ge 0.95$ , the error between the approximate physical field distribution and the real distribution is within acceptable bounds, that is,

$$V_j(\mathit{x},t)\approx {\overline{V}}_j(\mathit{x})+\sum _{j=1}^N{a}_j(t_i){\varphi }_j(\mathit{x})$$

(36)

By substituting equation (36) into equation (13), the impact of external wind on the quadrotor can be estimated through the POD method. In order to weaken the coupling relationship between the force and torque controlled by each motor, it is assumed that:

$$\left[\begin{array}{c}{U}_1\hfill \\ \multicolumn{1}{c}{U_2}\\ \multicolumn{1}{c}{U_3}\\ \multicolumn{1}{c}{U_4}\end{array}\right]=\left[\begin{array}{c}\multicolumn{1}{c}{-C_T\sum _{i=1}^4{\omega }_i^2}\\ \multicolumn{1}{c}{0.5C_T{l}_a(-{\omega }_3^2+{\omega }_1^2+{\omega }_2^2-{\omega }_4^2)}\\ \multicolumn{1}{c}{0.5C_T{l}_a({\omega }_3^2-{\omega }_1^2-{\omega }_4^2+{\omega }_2^2)}\\ \multicolumn{1}{c}{C_r({\omega }_1^2-{\omega }_2^2+{\omega }_3^2-{\omega }_4^2)}\end{array}\right]$$

(37)

Thus, the control equation for the quadrotor can be written as

$$\begin{array}{c}\stackrel{··}{x}=\frac{U_1}{m}{U}_x+\frac{D_{\mathrm{gx}}}{m}\hfill \\ \stackrel{··}{y}=\frac{U_1}{m}{U}_y+\frac{D_{\mathrm{gy}}}{m}\hfill \\ \stackrel{··}{z}=\frac{U_1}{m}\cos \theta \cos \varphi +g+\frac{D_{\mathrm{gz}}}{m}\hfill \\ \stackrel{··}{\varphi }=U_2/J_x+\stackrel{·}{\theta }\stackrel{·}{\psi }(J_y-J_z)/J_x\hfill \\ \stackrel{··}{\theta }=U_3/J_y+\stackrel{·}{\varphi }\stackrel{·}{\psi }(J_z-J_x)/J_y\hfill \\ \stackrel{··}{\psi }=U_4/J_z+\stackrel{·}{\varphi }\stackrel{·}{\theta }(J_x-J_y)/J_z\hfill \end{array}$$

(38)

where $U_x$ , $U_y$ represent virtual control inputs of quadrotors in the horizontal direction.

$$\begin{array}{c}{U}_x=\cos \varphi \sin \theta \cos \psi +\sin \varphi \sin \psi \hfill \\ U_y=\cos \varphi \sin \theta \sin \psi -\sin \varphi \cos \psi \hfill \end{array}$$

(39)

in which $D_{\mathrm{gx}},D_{\mathrm{gy}},D_{\mathrm{gz}}$ is the external disturbances, including air drag, which can be expressed as

$$\begin{array}{c}{U}_s=(\sin \varphi \sin \theta \cos \psi -\cos \varphi \sin \psi )\hfill \\ U_t=(\sin \varphi \sin \theta \sin \psi +\cos \varphi \cos \psi )\hfill \\ D_{\mathrm{gx}}=D_{\mathrm{bx}}\cos \theta \cos \psi +D_{\mathrm{by}}{U}_s+D_{\mathrm{bz}}{U}_x\hfill \\ D_{\mathrm{gy}}=D_{\mathrm{bx}}\cos \theta \sin \psi +D_{\mathrm{by}}{U}_t+D_{\mathrm{bz}}{U}_y\hfill \\ D_{\mathrm{gz}}=-D_{\mathrm{bx}}\sin \theta +D_{\mathrm{by}}\sin \varphi \cos \theta +D_{\mathrm{bz}}\cos \varphi \cos \theta \hfill \end{array}$$

(40)

where $D_{\mathrm{bx}},D_{\mathrm{by}},D_{\mathrm{bz}}$ are the aerodynamic drags in the body frame. From equation (13), $D_{\mathrm{bx}},D_{\mathrm{by}},D_{\mathrm{bz}}$ are directly proportional to the square of velocity. Moreover, this velocity $\mathit{V}$ can be rapidly estimated in real time using the POD method at time $t$ and location $\mathit{x}$ .

The three-axis position control of quadrotors relies on a single control input $U_1$ . If the quadrotor does not perform the rolling, pitching, or yawing actions, the quadrotor is constrained to vertical motion only. On the contrary, if a quadrotor performs any of the rolling, pitching, or yawing actions, it will disturb the equilibrium and cause displacement in the horizontal direction. As the control system of the quadrotor is an underactuated system, the dual closed-loop control strategy ²⁶ is adopted for the quadrotor. The inner loop is used to ensure that the quadrotor maintains a stable attitude during flight. The outer loop is used to ensure that the quadrotor accurately flies along the predetermined trajectory. The structural diagram of the entire system is as shown in Figure 2. $x_d,y_d,z_d,{\varphi }_d,{\theta }_d,{\psi }_d$ represent the desired position and attitude angle, and $x,y,z,\varphi ,\theta ,\psi$ represent the actual position and attitude angle. In the position controller, the control input $U_1$ is first determined by the z-axis acceleration, and then $U_x,U_y$ can be obtained by substituting $U_1$ into (38). Then, the desired control input is obtained through the attitude solver based on equation (39). Assuming that $\mathit{X}=[x_1,x_2,x_3,x_4,x_5,x_6,x_7,x_8,x_9,x_{10},x_{11},x_{12}]\in R^{12\times 1}$ where $x_1=\varphi ,x_2=\stackrel{·}{\varphi },x_3=\theta ,x_4=\stackrel{·}{\theta },x_5=\psi ,$

Figure 2.

Control diagram of the quadrotor.

$x_6=\stackrel{·}{\psi },x_7=x,x_8=\stackrel{·}{x},$ $x_9=y,x_{10}=\stackrel{·}{y},x_{11}=z,x_{12}=\stackrel{·}{z}$ . Then equation (38) can be rewritten as

$$\stackrel{·}{\mathit{X}}=\left[\begin{array}{c}\stackrel{·}{x_1}\hfill \\ \multicolumn{1}{c}{\stackrel{·}{x_2}}\\ \multicolumn{1}{c}{\stackrel{·}{x_3}}\\ \multicolumn{1}{c}{\stackrel{·}{x_4}}\\ \multicolumn{1}{c}{\stackrel{·}{x_5}}\\ \multicolumn{1}{c}{\stackrel{·}{x_6}}\\ \multicolumn{1}{c}{\stackrel{·}{x_7}}\\ \multicolumn{1}{c}{\stackrel{·}{x_8}}\\ \multicolumn{1}{c}{\stackrel{·}{x_9}}\\ \multicolumn{1}{c}{{\stackrel{·}{x}}_{10}}\\ \multicolumn{1}{c}{\begin{array}{c}{\stackrel{·}{x}}_{11}\hfill \\ \multicolumn{1}{c}{{\stackrel{·}{x}}_{12}}\end{array}}\end{array}\right]=\left[\begin{array}{c}{x}_2\hfill \\ \multicolumn{1}{c}{x_4{x}_6{a}_1+U_2/J_x}\\ \multicolumn{1}{c}{x_4}\\ \multicolumn{1}{c}{x_2{x}_6{a}_2+U_3/J_y}\\ \multicolumn{1}{c}{x_6}\\ \multicolumn{1}{c}{x_2{x}_4{a}_3+U_4/J_z}\\ \multicolumn{1}{c}{x_8}\\ \multicolumn{1}{c}{\frac{U_1}{m}{U}_x+\frac{D_{\mathrm{gx}}}{m}}\\ \multicolumn{1}{c}{x_{10}}\\ \multicolumn{1}{c}{\frac{U_1}{m}{U}_y+\frac{D_{\mathrm{gy}}}{m}}\\ \multicolumn{1}{c}{x_{12}}\\ \multicolumn{1}{c}{\frac{U_1}{m}\cos \theta \cos \varphi +g+\frac{D_{\mathrm{gz}}}{m}}\end{array}\right],$$

(41)

where $a_1=\frac{J_y-J_z}{J_x},a_2=\frac{J_z-J_x}{J_y},a_3=\frac{J_x-J_y}{J_z}$ .

The state equation of the quadrotor can be expressed as follows:

$$\{\begin{array}{c}{\mathit{X}}_{t+1}={\mathit{X}}_t+\mathit{\Delta }t\stackrel{·}{{\mathit{X}}_t}+{\mathit{W}}_t={\mathit{A}}_d{\mathit{X}}_t+{\mathit{B}}_d{\mathit{u}}_t+{\mathit{W}}_t\hfill \\ \multicolumn{1}{c}{{\mathit{y}}_t={\mathit{C}}_d{\mathit{X}}_t}\end{array},$$

(42)

in which, ${\mathit{X}}_t\in R^{12}$ is the state of the system at time $t$ ; ${\mathit{u}}_t\in R^4$ is control input; ${{\mathit{C}}_d}^{l\times 12}$ is the system measurement matrix; ${\mathit{y}}_t\in R^l$ is the system output. $W_t\in R^{12}$ is the system process noise, l and m are numbers related to system inputs and outputs, ${\mathrm{X}}_d$ is the identity matrix. ${\mathit{u}}_t\in R^4$ is control input designed based on the state equation.

${\mathit{A}}_d\in {R^{12\times }}^{12},{\mathit{B}}_d\in {R^{12\times }}^4$ can be obtained by taking the partial derivatives as follows

$${\mathit{A}}_d=\frac{\partial ({\mathit{X}}_t+\mathit{\Delta }t\stackrel{·}{{\mathit{X}}_t})}{\partial {\mathit{X}}_t},{\mathit{B}}_d=\frac{\partial ({\mathit{X}}_t+\mathit{\Delta }t\stackrel{·}{{\mathit{X}}_t})}{\partial u_t}$$

(43)

Thus, the prediction method is as follows. Firstly, open-loop prediction is performed based on the following equation

$$\begin{array}{c}{\hat{\mathit{X}}}_{t+i|t}={\mathit{A}}_d{\hat{\mathit{X}}}_{t+i-1|t}+{\mathit{B}}_d{\mathit{u}}_{t+i}\hfill \\ ={\mathit{A}}_d{\hat{\mathit{X}}}_{t+i-1|t}^1+\sum _{n=0}^i{\mathit{A}}_d^n{\mathit{B}}_d{\mathit{u}}_{t+n}\hfill \\ ={\mathit{A}}_d{\hat{\mathit{X}}}_{t+i|t}^1+\sum _{n=0}^i{\mathit{A}}_d^n{\mathit{B}}_d{\mathit{u}}_{t+n},1\le i\le p\hfill \end{array}$$

(44)

where $\mathit{f}(t)$ represents the predicted time range and ${\hat{\mathit{X}}}_{t+i-1|t}^1$ represents the part of $\mathit{f}(t)$ ,which is independent of current and future inputs. Based on the above prediction expression (44), the control input vector $\mathit{f}(t)$ can be obtained by minimizing the following indicator function for prediction,

$$minJ=({\mathit{R}}_t^x-{\mathit{X}}_t^v{)}^T({\mathit{R}}_t^x-{\mathit{X}}_t^v)+{\mathit{U}}_t^T\mathbf{\gamma }{\mathit{U}}_t$$

(45)

The optimal control input solution above is

$${\mathit{U}}_t^{*}=({\mathit{G}}_0^T{\mathit{G}}_0+\mathbf{\gamma }{)}^{-1}{\mathit{G}}_0^T({\mathit{R}}_t^x-{\mathit{X}}_t^1)$$

(46)

where $\mathit{f}(t)$ is the system prediction matrix, $\mathit{f}(t)$ is the predicted state vector, $\mathit{f}(t)$ is the known vector part of $\mathit{f}(t)$ , $\mathit{f}(t)$ is the control input weight matrix, and ${\mathit{R}}_t^x$ is the reference value of the system state. The vector composed of the first-row elements of the optimal control matrix ${\mathit{U}}_t^{*}$ is used as the control input vector. Then (45) and (46) are repeatedly used to estimate the system state during the continuous sampling period.

Predictive control method based on active model

The flight modes of the quadrotors can be divided into the hovering mode (speed below 1 m/s), the cruising mode (speed between 1 and 15 m/s), and the ground approach mode (takeoff and landing or flying close to the ground). In addition, the transitions between modes are not clearly bounded and involve complex transient processes. The dynamic characteristics of the aircraft in hovering mode are relatively stable. In this mode, the effect of the slight changes on the dynamic model is not significant. Because model identification and verification are effective in this mode, the hovering dynamic model is usually selected for the control of the quadrotor. When the quadrotor flies at high speeds or near the ground, the dynamic models of the quadrotor may have significant differences from those in hovering mode.

To describe model uncertainty, we consider linearization errors, internal and external disturbances, and unmodeled dynamics as additive process noise. This noise term not only depends on the system state vector and external disturbances, but also has a coupling effect with the inherent noise characteristics of the system, which needs to be compensated in the controller design. A nonlinear function $\mathit{f}(t)$ is introduced to describe the model difference during flight. Let $\mathit{X}\in R^{12\times 1}$ be the system state and let $\mathit{W}\in R^{12\times 1}$ be the process noise of the hovering model. The model difference description of this system can be expressed as:

$$\begin{array}{c}{B}_f\mathit{f}(t)=(\stackrel{·}{\stackrel{~}{\mathit{X}}})(t)-\stackrel{·}{\mathit{X}}(t)\hfill \\ \stackrel{·}{\mathit{f}}(t)=\mathit{h}(t)\hfill \end{array}$$

(47)

in which $\mathit{X}(t)$ is the state of the system reference model, $\tilde{\mathit{X}}(t)$ is the actual dynamic state of the system, ${\mathit{B}}_f$ is the coefficient adjustment matrix, which can be manually selected, and $\mathit{h}(t)$ is the rate of change of the model mismatch. The uncertain dynamics of actual systems can be expressed in the following form:

$$\begin{array}{c}\tilde{\mathit{X}}(t+1)={\mathit{A}}_t\mathit{X}(t)+{\mathit{B}}_t\mathit{U}(t)+{\mathit{B}}_{\mathit{f}}\mathit{f}(t)\hfill \\ \mathit{X}(t+1)={\mathit{A}}_t\mathit{X}(t)+{\mathit{B}}_t\mathit{u}(t)\hfill \\ \mathit{Y}(t+1)=\mathit{CX}(t)+V(t).\hfill \end{array}$$

(48)

The nominal system matrix ${\mathit{A}}_t$ and control matrix ${\mathit{B}}_t$ at time can be obtained by substituting time $t$ into (43). Define $\mathit{U}(t)$ as the final control input of a real system, define $\mathit{u}(t)$ as the input designed based on state equations. Define $\mathit{V}(t)$ , $\mathit{C}$ as the measurement noise of the system, the measurement matrix of the system, respectively. Define $\mathit{Y}(t)$ as the true output of the system, which can be obtained through sensors. In order to eliminate the influence of model differences described in equation (48) on control, we have

$$\Vert \tilde{\mathit{X}}(t)-\mathit{X}(t)\Vert =0.$$

(49)

The above equation can be rewritten as the following discretization form

$$\Vert {\mathit{B}}_t\mathit{U}(t)+{\mathit{B}}_{\mathit{f}}\mathit{f}(t)-{\mathit{B}}_t\mathit{u}(t)\Vert =0.$$

(50)

This study designs a controller based on the active model (as shown in Figure 3), which obtains compensation $\mathit{U}(t)$ through the optimal prediction method to maximize the performance of the controller. As the fact that the quadrotors are underactuated systems, the dimensionality of the system input $\mathit{U}(t)$ is much smaller than that of $\mathit{f}(t)$ . $\mathit{U}(t)$ may be obtained by using optimization methods. Equation (50) can be transformed into the following form:

$$\begin{array}{c}\multicolumn{1}{c}{{\mathit{U}}_t^{*}=\underset{{\mathit{U}}_t}{\arg min}{J}_t({\mathit{U}}_t)}\\ \multicolumn{1}{c}{J_t({\mathit{U}}_t)={({\mathit{B}}_t\mathit{U}(t)+{\mathit{B}}_{\mathit{f}}\mathit{f}(t)-{\mathit{B}}_t\mathit{u}(t))}^T\mathit{H}({\mathit{B}}_t\mathit{U}(t)}\\ \multicolumn{1}{c}{+{\mathit{B}}_{\mathit{f}}\mathit{f}(t)-{\mathit{B}}_t\mathit{u}(t)),}\end{array}$$

(51)

in which $\mathit{H}$ is a positive definite weight matrix that can be customized and adjusted. Based on the principle of adaptive group membership estimation, ²⁷ the defined indicators ${\delta }_{t+1}$ may be given by

$$\begin{array}{c}\multicolumn{1}{c}{{\delta }_{t+1}=1-{({\mathit{Y}}_{t+1}-{\mathit{C}}_d({\mathit{A}}_d{\hat{\mathit{X}}}_{\mathit{t}|t}+{\mathit{B}}_d{\mathit{U}}_{\mathit{t}}))}^T}\\ \multicolumn{1}{c}{{\mathit{W}}_{\mathit{t}}^{-1}({\mathit{Y}}_{t+1}-{\mathit{C}}_d({\mathit{A}}_d{\hat{\mathit{X}}}_{\mathit{t}|t}+{\mathit{B}}_d{\mathit{U}}_{\mathit{t}}))}\end{array}$$

(52)

where ${\delta }_{t+1}$ denotes the effectiveness indicator of model prediction. When ${\delta }_{t+1}>0$ , the prediction is effective. Then we have

$$\begin{array}{c}\multicolumn{1}{c}{J_t^{\delta }(U_t,Y_{t+1})={({\mathit{Y}}_{t+1}-{\mathit{C}}_d({\mathit{A}}_d{\hat{\mathit{X}}}_{\mathit{t}|t}+{\mathit{B}}_d{\mathit{U}}_{\mathit{t}}))}^T}\\ \multicolumn{1}{c}{{\mathit{W}}_t^{-1}({\mathit{Y}}_{t+1}-{\mathit{C}}_d({\mathit{A}}_d{\hat{\mathit{X}}}_{\mathit{t}|t}+{\mathit{B}}_d{\mathit{U}}_{\mathit{t}})).}\end{array}$$

(53)

Then we have an optimization problem as follows

$$\begin{array}{c}\multicolumn{1}{c}{J_t^{\delta *}(U_t,Y_{t+1})=min{({\mathit{Y}}_{t+1}-{\mathit{C}}_d({\mathit{A}}_d{\hat{\mathit{X}}}_{\mathit{t}|t}+{\mathit{B}}_d{\mathit{U}}_{\mathit{t}}))}^T}\\ \multicolumn{1}{c}{{\mathit{W}}_t^{-1}({\mathit{Y}}_{t+1}-{\mathit{C}}_d({\mathit{A}}_d{\hat{\mathit{X}}}_{\mathit{t}|t}+{\mathit{B}}_d{\mathit{U}}_{\mathit{t}})).}\end{array}$$

(54)

By taking into account the index of equation system correction (51) and the stability index of system prediction (52), the following weighted optimization index may be given by

$$\begin{array}{c}\underset{U_t}{min}{\overline{J}}_t({\mathit{U}}_t,{\mathit{Y}}_{t+1})\hfill \\ {\overline{J}}_t({\mathit{U}}_t)=J_t({\mathit{U}}_t)+k_{\alpha }{J}_t^{\delta *}({\mathit{U}}_t,{\mathit{Y}}_{t+1}).\hfill \end{array}$$

(55)

where $k_{\alpha }$ is the weight coefficient selected independently. If satisfied by the optimal solution $\frac{\partial {\overline{J}}_t(U_t)}{\partial U_t}=0$ , the optimal solution of the above equation can be obtained as

$$\begin{array}{c}{U}_{\mathit{t}}^{*}(Y_{t+1})={({\mathit{B}}_d^T\mathit{H}{\mathit{B}}_d+k_{\alpha }{\mathit{B}}_d^T{\mathit{C}}_d^T{\mathit{W}}_t^{-1}{\mathit{C}}_{\mathit{d}}{\mathit{B}}_d)}^{-1}\hfill \\ \multicolumn{1}{c}{[k_{\alpha }{\mathit{B}}_d^T{\mathit{C}}_d^T{\mathit{W}}_t^{-1}({\mathit{Y}}_{t+1}-{\mathit{C}}_d{\mathit{A}}_d{\hat{\mathit{X}}}_{t|t}{)-B}_d^T\mathit{H}({\mathit{B}}_f{\hat{\mathit{f}}}_t-{\mathit{B}}_d{\mathit{u}}_t)].}\end{array}$$

(56)

where ${\hat{\mathit{f}}}_{\mathit{t}}$ is the prediction vector for the model difference, which can be obtained based on the dynamic model (41) and the relative wind speed prediction equation (36). This solution ignores the effect of ${\mathit{Y}}_{t+1}$ , here ${\mathit{Y}}_{t+1}\in [\mathit{Y}(t),\mathit{C}{\hat{\mathit{X}}}_{\mathit{t}|t}]$ belongs to a certain numerical value between the current state and the predicted state. This can be solved by substituting $U_{\mathit{t}}^{*}(Y_{t+1})$ into (48) to obtain the following linear optimization problem, and then obtain ${\mathit{Y}}_{t+1}^{*}$ of the optimal solution

$$\begin{array}{c}\underset{U_t}{min}{\overline{J}}_t({\mathit{U}}_t,{\mathit{Y}}_{t+1})=\underset{U_t}{min}{\overline{J}}_t(U_{\mathit{t}}^{*}(Y_{t+1}),{\mathit{Y}}_{t+1}).\hfill \\ s.t.\mathit{Y}(t+1)\in [\mathit{Y}(t),\mathit{C}{\hat{\mathit{X}}}_{\mathit{t}|t}]\hfill \end{array}$$

(57)

Set $k_{\alpha }$ , $B_f=\mathit{I}$ , the control correction amount can be calculated by the following formula

$$U_t^{*}={\left({\mathit{B}}_d^T\mathit{H}{\mathit{B}}_d\right)}^{-1}{\mathit{B}}_d^T\mathit{H}({\mathit{B}}_f{\hat{\mathit{f}}}_t-{\mathit{B}}_d{\mathit{u}}_t).$$

(58)

Assuming that

$${\mathit{V}}_t={\tilde{\mathit{X}}}_t^T{\tilde{\mathit{X}}}_t.$$

(59)

By substituting equation (62) into equation (48), we have

$${\tilde{\mathit{X}}}_{t+1}={\mathit{X}}_{t+1}+(\mathit{I}-Z){\mathit{f}}_t{-Z(K}_t^e{\mathit{V}}_t-{\mathit{h}}_t).$$

(60)

where ${\mathit{X}}_{\mathit{t}+1}$ is the abbreviation for $\mathit{X}(t+1)$ , ${\mathit{h}}_t$ is the abbreviation for $\mathit{h}(t)$ . ${\mathit{K}}_{\mathit{t}}^{\mathit{e}}$ and $\mathit{Z}={\mathit{B}}_d{\left({\mathit{B}}_d^T\mathit{H}{\mathit{B}}_d\right)}^{-1}{\mathit{B}}_d^T\mathit{H}$ is the gain for the model difference in the adaptive set membership estimator, ${\mathit{K}}_{\mathit{t}}^{\mathit{e}}$ is generally positive. Then we have

$$\begin{array}{c}{\mathit{V}}_{t+1}={\tilde{\mathit{X}}}_t^T{\tilde{\mathit{X}}}_t=(\overline{\mathit{A}}{\mathit{X}}_t+(\mathit{I}-Z){\mathit{f}}_{\mathit{t}}-\mathit{Z}({\mathit{K}}_t^e{\mathit{V}}_t-{\mathit{h}}_t){)}^T\hfill \\ \multicolumn{1}{c}{(\overline{\mathit{A}}{\mathit{X}}_t+(\mathit{I}-Z){\mathit{f}}_{\mathit{t}}-\mathit{Z}({\mathit{K}}_t^e{\mathit{V}}_t-{\mathit{h}}_t))\le {\mathit{X}}_{\mathit{t}}^{\mathit{T}}{\overline{\mathit{A}}}^{\mathit{T}}\overline{\mathit{A}}{\mathit{X}}_{\mathit{t}}}\\ \multicolumn{1}{c}{{+2X}_{\mathit{t}}^{\mathit{T}}{\overline{\mathit{A}}}^{\mathit{T}}(I-Z){\mathit{f}}_{\mathit{t}}{+f}_{\mathit{t}}^{\mathit{T}}{(I-Z)}^{\mathit{T}}(I-Z){\mathit{f}}_{\mathit{t}}{-2X}_{\mathit{t}}^{\mathit{T}}{\overline{\mathit{A}}}^{\mathit{T}}\mathit{Z}{\mathit{f}}_{\mathit{t}}}\\ \multicolumn{1}{c}{{-2f}_{\mathit{t}}^{\mathit{T}}{(I-Z)}^{\mathit{T}}\mathit{Z}{\mathit{f}}_{\mathit{t}},{+f}_t^T{\mathit{Z}}^T\mathit{Z}{\mathit{f}}_{\mathit{t}}{+V}_t^T{\mathit{K}}_t^{\mathrm{eT}}{\mathit{Z}}^T{\mathit{ZK}}_t^e{\mathit{V}}_t{+h}_t^T{\mathit{h}}_t}\end{array}$$

(61)

in which $\overline{\mathit{A}}$ is the system matrix of the equivalent nominal model, satisfied

$${\mathit{A}}_d{\mathit{X}}_t+{\mathit{B}}_d{\mathit{u}}_t=\overline{\mathit{A}}{\mathit{X}}_t$$

(62)

Choose $\mu >0$ , we obtain

$${\mathit{V}}_{\mathit{t}}^{\mathit{T}}{\mathit{K}}_{\mathit{t}}^{\mathit{eT}}{\mathit{Z}}^{\mathit{T}}{\mathit{ZK}}_{\mathit{t}}^{\mathit{e}}{\mathit{V}}_{\mathit{t}}{+h}_{\mathit{t}}^{\mathit{T}}{\mathit{h}}_{\mathit{t}}\le \mu {\mathit{X}}_{\mathit{t}}^{\mathit{T}}{\overline{\mathit{A}}}^{\mathit{T}}\overline{\mathit{A}}{\mathit{X}}_{\mathit{t}}$$

(63)

Figure 3.

The control process based on active modelling.

Then

$$\begin{array}{c}{\mathit{V}}_{t+1}={\tilde{\mathit{X}}}_t^T{\tilde{\mathit{X}}}_t=\le {\mathit{X}}_{\mathit{t}}^{\mathit{T}}{\overline{\mathit{A}}}^{\mathit{T}}\overline{\mathit{A}}{\mathit{X}}_{\mathit{t}}{+X}_{\mathit{t}}^{\mathit{T}}{\overline{\mathit{A}}}^{\mathit{T}}(2I-4Z){\mathit{f}}_{\mathit{t}}\hfill \\ {+f}_{\mathit{t}}^{\mathit{T}}{(\mathit{I}-Z)}^{\mathit{T}}(\mathit{I}-3Z){\mathit{f}}_{\mathit{t}}+\mu {\mathit{X}}_{\mathit{t}}^{\mathit{T}}{\overline{\mathit{A}}}^{\mathit{T}}\overline{\mathit{A}}{\mathit{X}}_{\mathit{t}}\hfill \end{array}$$

(64)

Choose $\mathit{H}$ to satisfy

$$\Vert \mathit{Z}\Vert =\Vert {\mathit{B}}_d{{(\mathit{B}}_d^T\mathit{H}{\mathit{B}}_d)}^{-1}{\mathit{B}}_d^T\mathit{H}\Vert \le 1/\lambda$$

(65)

where $\lambda >0$ , $\Vert {\mathit{f}}_t\Vert$ satisfies

$$\Vert {\mathit{f}}_t\Vert \le \lambda \Vert \overline{\mathit{A}}{\mathit{X}}_{\mathit{t}}\Vert$$

(66)

Then we have

$$\begin{array}{c}{\tilde{\mathit{V}}}_{t+1}-{\tilde{\mathit{V}}}_t={\mathit{X}}_{\mathit{t}}^{\mathit{T}}({\overline{\mathit{A}}}^{\mathit{T}}\overline{\mathit{A}}-I){\mathit{X}}_{\mathit{t}}+\lambda \Vert {\mathit{X}}_{\mathit{t}}^{\mathit{T}}{\overline{\mathit{A}}}^{\mathit{T}}(2I-4Z)\overline{\mathit{A}}{\mathit{X}}_{\mathit{t}}\Vert \hfill \\ +{\mathbf{\lambda }}^{\mathit{2}}\Vert {\mathit{X}}_{\mathit{t}}^{\mathit{T}}{\overline{\mathit{A}}}^{\mathit{T}}{(I-Z)}^{\mathit{T}}(I-3Z)\overline{\mathit{A}}{\mathit{X}}_{\mathit{t}}\Vert {+X}_{\mathit{t}}^{\mathit{T}}{\overline{\mathit{A}}}^{\mathit{T}}\mathit{\mu }\overline{\mathit{A}}{\mathit{X}}_{\mathit{t}}\hfill \\ {\mathit{X}}_{\mathit{t}}^{\mathit{T}}({\overline{\mathit{A}}}^{\mathit{T}}\overline{\mathit{A}}+\lambda {\overline{\mathit{A}}}^{\mathit{T}}\Vert (2I-4Z)\Vert \overline{\mathit{A})}{\mathit{X}}_{\mathit{t}}-{\mathbf{\lambda }}^{\mathit{2}}{\overline{\mathit{A}}}^{\mathit{T}}\Vert {(I-Z)}^{\mathit{T}}(I-3Z)\Vert \overline{\mathit{A})}{\mathit{X}}_{\mathit{t}}\hfill \\ {+X}_{\mathit{t}}^{\mathit{T}}{\overline{\mathit{A}}}^{\mathit{T}}\mathit{\mu }\overline{A}{\mathit{X}}_{\mathit{t}}=(-3-{\mathbf{\lambda }}^{\mathit{2}}+3/{\mathbf{\lambda }}^{\mathit{2}}{+\mu )X}_{\mathit{t}}^{\mathit{T}}({\overline{\mathit{A}}}^{\mathit{T}}\overline{\mathit{A})}{\mathit{X}}_{\mathit{t}}\hfill \end{array}$$

(67)

If $\lambda >0$ is large enough, we obtain

$${\tilde{\mathit{V}}}_{t+1}-{\tilde{\mathit{V}}}_t\le (-3-{\lambda }^2+3/{\lambda }^2+\mu ){\mathit{X}}_{\mathit{t}}^{\mathit{T}}({\overline{\mathit{A}}}^{\mathit{T}}\overline{\mathit{A})}{\mathit{X}}_{\mathit{t}}>0$$

(68)

This proves that the system is asymptotically stable.

Numerical examples and discussion

Based on field-collected gust data, the POD method is employed for gust prediction. The prediction accuracy is evaluated by comparing the predicted results with actual measurements, as illustrated in Figures 4 to 7. The results demonstrate that the proposed approach achieves accurate gust prediction, validating the effectiveness of the POD-based methodology.

Figure 4.

Comparison of wind field data and prediction data along X axis of roadway.

Figure 5.

Comparison of wind field data and prediction data along Y axis of roadway.

Figure 6.

Comparison of wind field data and prediction data along Z axis of roadway.

Figure 7.

Comparison of wind field data and prediction data of total velocity along roadway.

The above gust data were adopted as inputs to simulate hover control performance. Figures 8 to 10 demonstrate that the proposed active model control method effectively attenuates gust disturbances. Quantitative comparison reveals position deviations of merely 0.1 m in both X and Y directions for the proposed method—substantially lower than PID (0.35 m) and adaptive control (0.17 m). In the Z direction, the proposed method maintains altitude without significant descent, whereas PID control may cause altitude loss, posing significant safety risks.

Figure 8.

X-axis drift under the control of different controllers affected by gust during hovering.

Figure 9.

Y-axis drift of different controllers affected by gust during hovering.

Figure 10.

Z-axis drift of different controllers affected by gust during hovering.

The above gust data were employed as inputs to simulate hover control performance. Figures 11 to 13 demonstrate that the proposed active model control method effectively attenuates gust-induced disturbances. Quantitative comparison reveals position deviations of merely 0.1 m in X and Y directions—substantially lower than PID (0.35 m) and adaptive control (0.17 m). In the Z direction, the proposed method maintains altitude without significant descent, whereas PID control may cause altitude loss, posing significant safety risks.

Figure 11.

X-axis drift of constant speed flight by different controllers affected by gust.

Figure 12.

Y-axis drift of constant speed flight by different controllers affected by gust.

Figure 13.

Z-axis drift of constant speed flight by different controllers affected by gust.

Under gust-free conditions, the proposed method, PID and adaptive controllers all achieve accurate elliptical trajectory tracking (Figures 14 and 15). Figures 16 and 17 illustrate the tracking performance under gust disturbances. The proposed method maintains precise path following with minimal deviation, whereas the conventional PID controller exhibits significant trajectory fluctuations. The adaptive control method demonstrates intermediate performance.

Figure 14.

The flight trajectories without the gust disturbance.

Figure 15.

The RMS values of flight trajectory errors without the gust disturbance.

Figure 16.

The flight trajectories with the gust disturbance.

Figure 17.

The RMS values of flight trajectory errors with the gust disturbance.

The PID controller operates via a passive feedback mechanism that reacts only after errors occur, suffering from integral windup and phase lag that induce substantial tracking error and significant overshoot. The noise sensitivity of the derivative term necessitates filtering, which exacerbates phase lag. By contrast, the active model control method leverages precise system models to predict future state trajectories, enabling anticipatory feedforward compensation and proactive disturbance rejection for zero-phase-error tracking with minimal overshoot. Explicit constraint handling within the optimization framework facilitates smooth control planning and steady-state error mitigation. Adaptive control methods exhibit performance intermediate between the reactive PID and the predictive active model control. Temporal delays in parameter convergence, coupled dual-loop dynamics between estimation and state regulation, and transient model inaccuracies during adaptation collectively degrade transient response, yielding moderate overshoot superior to PID yet inferior to the proposed method.

The proposed methodology comprises three core modules: POD, ASMF, and online optimization. Real-time feasibility on resource-constrained underground platforms demands rigorous complexity analysis. The POD module serves as a dimensionality reduction tool. Computational costs concentrate in the offline phase, where singular value decomposition of the snapshot matrix incurs O( $W{M}^2$ or $W^{2}M$ ) complexity,where W represents the number of time steps. The online phase requires only matrix-vector multiplication between the state vector and pre-stored basis functions, thereby reducing complexity to O( $\mathrm{WM}$ ). When $W<<M$ this module incurs negligible runtime overhead, with memory occupancy maintained within several 100 kb—well-suited for storage-constrained embedded environments. The ASMF module encompasses sliding surface computation, state estimation, and adaptive gain updating, with per-step complexity of O(M²)–O(M³). Costs arise mainly from observer updates and integrating high-frequency switching dynamics. Substituting full adaptive laws with fixed gains or introducing boundary layers to suppress chattering substantially alleviates computational load, thus enabling millisecond-scale control periods. However, the online optimization module constitutes the principal bottleneck. Conventional Model Predictive Control (MPC) with rolling optimization necessitates solving quadratic programming problems at each step, entailing O( $W^3{M}^3$ ) complexity—prohibitive for ARM Cortex-M processors (<500 MHz). Explicit MPC, wherein polytopic partition control laws are computed offline, reduces online operations to table lookup or linear interpolation. Through algorithmic simplification coupled with hardware/software co-optimization—such as enabling hardware floating-point units and leveraging CMSIS-DSP libraries—the proposed method achieves real-time execution on resource-constrained underground platforms.

Conclusions

This study addresses the quadrotor control challenge in complex underground environments. A six-degree-of-freedom nonlinear dynamic model is developed. To mitigate gust effects, a real-time wind field prediction model is constructed based on wind speed sensor data and the POD method. In this framework, the prediction model and ASMF algorithm are combined, through which a multi-step incremental prediction mechanism is introduced. An active model correction framework is formulated as a constrained optimal control problem with online bias identification. By solving the above optimal control problem for the compensatory control inputs in real-time, effective control of the quadrotor in underground environments is achieved. Furthermore, convergence of the compensation algorithm is rigorously proved. Finally, the effectiveness of the algorithm is verified through numerical examples.