Next Article in Journal
Lane Level Positioning Method for Unmanned Driving Based on Inertial System and Vector Map Information Fusion Applicable to GNSS Denied Environments
Next Article in Special Issue
Fault-Tolerant Event-Triggrred Control for Multiple UAVs with Predefined Tracking Performance
Previous Article in Journal
Optimal Position and Target Rate for Covert Communication in UAV-Assisted Uplink RSMA Systems
Previous Article in Special Issue
Multi-Conflict-Based Optimal Algorithm for Multi-UAV Cooperative Path Planning
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Robust Control for UAV Close Formation Using LADRC via Sine-Powered Pigeon-Inspired Optimization

State Key Laboratory of Virtual Reality Technology and Systems, School of Automation Science and Electrical Engineering, Beihang University (BUAA), Beijing 100083, China
*
Author to whom correspondence should be addressed.
Drones 2023, 7(4), 238; https://doi.org/10.3390/drones7040238
Submission received: 25 February 2023 / Revised: 21 March 2023 / Accepted: 25 March 2023 / Published: 29 March 2023
(This article belongs to the Special Issue Swarm Intelligence in Multi-UAVs)

Abstract

:
This paper designs a robust close-formation control system with dynamic estimation and compensation to advance unmanned aerial vehicle (UAV) close-formation flights to an engineer-implementation level. To characterize the wake vortex effect and analyze the sweet spot, a continuous horseshoe vortex method with high estimation accuracy is employed to model the wake vortex. The close-formation control system will be implemented in the trailing UAV to steer it to the sweet spot and hold its position. Considering the dynamic characteristics of the trailing UAV, the designed control system is divided into three control subsystems for the longitudinal, altitude, and lateral channels. Using linear active-disturbance rejection control (LADRC), the control subsystem of each channel is composed of two cascaded first-order LADRC controllers. One is responsible for the outer-loop position control and the other is used to stabilize the inner-loop attitude. This control system scheme can significantly reduce the coupling effects between channels and effectively suppress the transmission of disturbances caused by the wake vortex effect. Due to the cascade structure of the control subsystem, the correlation among the control parameters is very high. Therefore, sine-powered pigeon-inspired optimization is proposed to optimize the control parameters for the control subsystem of each channel. The simulation results for two UAV close formations show that the designed control system can achieve stable and robust dynamic performance within the expected error range to maximize the aerodynamic benefits for a trailing UAV.

1. Introduction

In recent years, there has been growing interest in the close formation of unmanned aerial vehicles (UAVs) within the UAV formation community [1,2,3,4,5]. Its potential benefits include improving coordination, reducing radar cross-sections, enhancing obstacle avoidance, and saving energy. In close-formation flight, the trailing UAV can utilize the upwash of the wake vortex induced by the leading UAV to increase lift and reduce drag [6]. As a result, the trailing UAV can save energy, which has been demonstrated through theoretical analysis [7], observed in wind-tunnel experiments [8], and confirmed by flight tests [9]. These works also indicate that a robust formation flight-control system with high accuracy and excellent performance is crucial to both the implementation of UAV close-formation flight and maximization of the benefits of the formation aerodynamics. This is because the flight stability of the trailing UAV significantly deteriorates after being disturbed by the wake vortex effect, which can lead to flight safety issues. Furthermore, the region within 10 % of the wing span of the optimal relative position is defined as the position range of the trailing UAV allowed by close-formation flight. Once the trailing UAV is unable to remain in the region and hold its position, more than 30 % of the benefits of the formation aerodynamics will be lost [6]. It should be noted that the optimal relative position is also called the sweet spot of the trailing UAV in close-formation flight.
The design of robust close-formation control systems has been investigated using various methods, including adaptive control [10], backstepping control [11], sliding-mode control [12], robust control [13], etc. In adaptive control, the predictions for the wake vortex effect must be availabel to counteract it online. However, the adaptive law is only able to address the matched portion of the wake vortex effect. Backstepping control is generally employed to stabilize the nominal model of the trailing UAV since the method itself is not robust. To further realize the dynamic compensation for the wake vortex effect, backstepping control needs to be combined with a disturbance estimator, which will complicate the structure of the control system and increase the difficulty of tuning the control parameters. Compared with the above two nonlinear control methods, sliding-mode control and robust control are robust to the wake vortex effect. However, they require a priori knowledge of the gradients and boundaries of the wake vortex effect for proper control gains to guarantee stability. Furthermore, the robust close formation control systems developed using these two control methods cannot work in the optimal state most of the time since they are unable to dynamically compensate for the wake vortex effect. Therefore, these designs are relatively conservative for close formation flight.
In this paper, we aim to develop a robust close-formation control system with dynamic estimation and compensation to advance UAV close-formation flight to an engineer-implementation level. To characterize the wake vortex effect and analyze the sweet spot, a continuous-horseshoe vortex method with high estimation accuracy is employed to model the wake vortex [14]. Furthermore, the computational efficiency of this modeling method is sufficient for real-time implementation. The estimated wake vortex effect is integrated into a nonlinear high-fidelity UAV model to describe the dynamic characteristics of the trailing UAV under the disturbance of the wake vortex effect. Active-disturbance rejection control (ADRC) is a new feedback linearization method proposed by Han that is based on the classical PID control principle [15]. The core idea of ADRC is to consider the various internal and external uncertainties of the control object as the total disturbance. Thereafter, an extended-state observer (ESO) and state feedback controller are constructed to dynamically estimate and compensate for the total disturbance, respectively. This control method has promising engineering potential as it does not use any prior information about the disturbance and requires only minimal model knowledge of the control object [16,17]. However, the design of a nonlinear ADRC has some difficulties in parameter tuning and theory analysis due to the existence of strong nonlinear functions. Gao simplified the nonlinear ADRC to the linear ADRC (LADRC) [18] to make it more suitable for engineering applications. The convergence of the linear ESO (LESO) and the stability of the closed-loop LADRC have been proven under the assumption that the differentiation of the total disturbance is bounded [19]. Motivated by the facts stated above, we consider using the LADRC to design a robust close-formation control system. In close-formation flight, it is assumed that the leading UAV maintains a level and stable flight path while the trailing UAV is far away from the sweet spot at the beginning. The robust close-formation control system will be implemented in the trailing UAV to steer it to the sweet spot and hold its position. In light of the dynamic characteristics of the trailing UAV, the control system is divided into three control subsystems for the longitudinal, altitude, and lateral channels. Considering the strength of the wake vortex effect, the control subsystem of each channel is composed of two cascade controllers, namely an outer-loop position LADRC controller and an inner-loop attitude LADRC controller. According to the bandwidth theory proposed by Gao, the estimation and compensation of the LADRC controller depend on the bandwidth of the LESO and the gain of the state feedback controller, respectively. Furthermore, these two control parameters of the inner- and outer-loop LADRC controller must be considered together to balance their roles in the cascade structure. The manual tuning method of the control parameters is not only time-consuming but also challenging to achieve optimal control performance, especially when knowledge and experience are limited.
Unlike the manual tuning method, the evolutionary algorithm has proven to be feasible and effective for the parameter tuning of various research problems [20,21]. Pigeon-inspired optimization (PIO) was proposed by Duan [22] as a special evolutionary algorithm inspired by pigeon-homing behavior. Because of its easy implementation and fast convergence, PIO has been widely used to search for the optimal values of control parameters [23,24,25]. However, PIO often encounters the local optimum and prematurely converges since the appropriate adjustment of global exploration and local exploitation is not considered. For this reason, sine-powered controlled inertia weights are employed as an adjustment strategy to overcome the local optimum and enhance the global searching capability [26]. This improved pigeon-inspired optimization algorithm is known as SCPIO. Therefore, we utilize SCPIO to simultaneously tune the control parameters of the inner- and outer-loop LADRC controllers for the control subsystem of each channel. Using SCPIO, the developed robust close-formation control system can guarantee robust dynamic and error-bounded stable performance. Furthermore, the developed robust close-formation control system follows the wake vortex model with high estimation accuracy and the nonlinear high-fidelity UAV model, making it more effective and reliable for engineering applications.
The rest of this paper is organized as follows. Section 2 presents the wake vortex model and the six-degrees-of-freedom (6-DOF) nonlinear dynamic model of the trailing UAV under the disturbance of the wake vortex effect. Section 3 introduces the basic principle of the first-order LADRC. The robust close-formation control system with dynamic estimation and compensation is designed using the first-order LADRC in Section 4. The SCPIO algorithm-based control system parameter optimization is described in Section 5. The simulation results are presented in Section 5. Finally, the conclusions are presented in Section 6.

2. Close Formation Modeling

In close-formation flight, it is assumed that the leading UAV maintains a level and stable flight path and the trailing UAV is far away from the sweet spot at the beginning. The robust close-formation control system will be implemented in the trailing UAV to steer it to the sweet spot and hold its position, as shown in Figure 1. Under this condition, the model description of the whole formation can be broken down into two separate parts. The first is the wake vortex model, which should be able to accurately describe the velocity field of the wake vortex generated by the leading UAV. The second is the trailing UAV model, which can characterize the dynamics of a UAV affected by the wake vortex.

2.1. Wake Vortex Model

Taking into account the estimation accuracy and computational efficiency for real-time implementation, a continuous-horseshoe vortex method is employed to model the wake vortex. According to its modeling approach, the wake vortex comprises an infinite number of continuously distributed semi-infinite horseshoe vortices. Considering the structure of the horseshoe vortex, the wake vortex is divided into a bound vortex and a free-trailing vortex. The bound vortex adheres to the wing surface and the filaments follow along the quarter-chord line. The free-trailing vortex falls off the wing surface and the filaments extend downstream to infinity parallel to the velocity vector of the leading UAV.
In the case of an approximately rectangular wing, the lift of a UAV presents an elliptic distribution along the quarter-chord line. Therefore, the circulation distribution of the wake vortex is assumed to be
Γ ( y c ) = Γ 0 1 2 y c b 2 , ( b / 2 y c b / 2 )
where y c is the lateral coordinate of the point on the quarter-chord line and Γ 0 = 2 V S C L / ( b π ) is the circulation at the wing root, with  V being the magnitude of the free-airflow velocity, S and b being the area and span of the wing, and  C L denoting the lift coefficient. For an arbitrary point P ( x , y , z ) inside the wake vortex, the velocity induced by a single straight vortex filament is obtained as follows according to the Kutta–Joukowski theorem and the Biot–Savart law,
v ( x , y , z ) = Γ ( y c ) 4 π h ( c o s θ 0 c o s θ ) · κ · n
where θ 0 and θ are the initial and final angles formed by point P and the vortex filament, h is the perpendicular distance from point P to the vortex filament, κ denotes the decrease in the strength of the wake vortex, and n is the unit vector of the induced velocity.
Following the continuous-horseshoe vortex method, the magnitude of the wake velocity is calculated by integrating the induced velocity of the single straight vortex filament along the quarter-chord line. Therefore, the body frame F = { O B , X B , Y B , Z B } of the leading UAV is established to determine the location and orientation of the vortex filament. As illustrated in Figure 2, the origin O B is defined at the center of gravity (CG) of the leading UAV. X B points toward the nose of the leading UAV, Y B points toward the right wing, and Z B points toward the fuselage belly. It should be noted that the body frame is fixed to the leading UAV. On this basis and taking into account the geometric relationship between the angle and position in (2), the velocity components at point P in the body frame of the leading UAV are given by
V f Y B = b / 2 b / 2 κ · μ · Γ ( y c ) z 4 π [ ( y y c ) 2 + z 2 + r c 2 ] · 1 x x 2 + ( y y c ) 2 + z 2 · d y c V f Z B = b / 2 b / 2 κ · μ · Γ ( y c ) ( y y c ) 4 π [ ( y y c ) 2 + z 2 + r c 2 ] · 1 x x 2 + ( y y c ) 2 + z 2 · d y c
V b X B = b / 2 0 κ · μ · Γ ( y c ) z 4 π ( x 2 + z 2 ) · y y c x 2 + ( y y c ) 2 + z 2 y + y c x 2 + ( y + y c ) 2 + z 2 · d y c V b Z B = b / 2 0 κ · μ · Γ ( y c ) x 4 π ( x 2 + z 2 ) · y y c x 2 + ( y y c ) 2 + z 2 y + y c x 2 + ( y + y c ) 2 + z 2 · d y c
where V f Y B and V f Z B are induced by all filaments of the free-trailing vortex on the Y B -axis and Z B -axis, respectively; V b X B and V b Z B are induced by all filaments of the bound vortex on the X B -axis and Z B -axis, respectively; r c is the core radius of the free-trailing vortex introduced to eliminate the possible singular problem; and μ is the interaction coefficient of filaments. The wake velocity at point P is defined as V B = [ V X B , V Y B , V Z B ] , where V X B , V Y B , and  V Z B are the velocity components on the X B , Y B , and  Z B axes, respectively. Therefore, V X B = V b X B , V Y B = V f Y B , and  V Z B = V b Z B + V f Z B . Note that V X B , V Y B , and  V Z B are also known as the backwash, sidewash, and upwash velocities, respectively.
The wake vortex effect refers to the forces and moments generated by the wake vortex of the leading UAV on the trailing UAV. Therefore, the wake vortex effect can be estimated according to the change in the attack angle for the trailing UAV. Let Δ α represent the increment in the attack angle at an arbitrary point on the quarter-chord line of the trailing wing, which can be obtained using the wake velocity calculated above [14]. Note that the aerodynamic and geometric parameters used in the following equation to estimate the wake vortex effect are for the trailing UAV.
Based on a statistical strategy, the induced lift and drag are written as
Δ L = q S C L α Δ α ¯ Δ D = q S C L + C L α Δ α ¯ 2 C L 2 / π A R
where Δ α ¯ = 1 N i = 1 N Δ α i is the average increment of the attack angle, with N being the number of statistical points; q represents the dynamic pressure; C L is the lift coefficient; C L α is the lift curve slope; S is the wing area; and A R is the aspect ratio of the wing that can be calculated by A R = b 2 / S .
Unlike the induced forces, the arm of force at the statistical point must be taken into account for the calculation of the induced moments. Therefore, one has
Δ L = q S 1 N i = 1 N C L α Δ α i y c i Δ M = q S 1 N i = 1 N C L α Δ α i x 0 | y c i | t a n Λ
where Δ L and Δ M are the induced rolling and pitching moments, respectively; y c i is the lateral arm of force at the i-th statistical point; Λ is the sweep angle of the wing; and x 0 is the lateral coordinate of the aerodynamic center.
The wake vortex effect is introduced into the dynamic and kinematic equations of the trailing UAV in the form of the induced forces and moments. It should be noted that this statistical strategy can improve the calculation efficiency while ensuring estimation accuracy compared to the continuous calculation method.
Remark 1. 
The generated side force and yawing moment are mainly caused by the sidewash velocity acting on the vertical tail of the trailing UAV. However, both the sidewash velocity and the area of the vertical tail are very small. Therefore, the induced side force and yawing moment can be ignored.

2.2. Trailing UAV Model

The highly nonlinear wake vortex effect encountered by the trailing UAV eliminates the possibility of using a simple linear model to achieve an accurate characterization of UAV dynamics while flying in close formation. Therefore, the 6-DOF nonlinear high-fidelity aircraft model, which was utilized in [27], is employed. To incorporate the wake vortex effect, we modified the dynamic and kinematic equations of the trailing UAV,
x ˙ E = ( u + Δ u ) ( cos θ cos ψ ) + ( v + Δ v ) ( sin ϕ cos ψ sin θ cos ϕ sin ψ ) + ( w + Δ w ) ( cos ϕ sin θ cos ψ + sin ϕ sin ψ ) y ˙ E = ( u + Δ u ) ( sin ψ cos θ ) + ( v + Δ v ) ( sin ψ sin θ sin ϕ + cos ψ cos ϕ ) + ( w + Δ w ) ( sin ψ sin θ cos ϕ cos ψ sin ϕ ) z ˙ E = ( u + Δ u ) sin θ ( v + Δ v ) ( cos θ sin ϕ ) ( w + Δ w ) ( cos θ cos ϕ )
V ˙ = 1 m ( ( D + Δ D ) + F T cos α cos β ) + g ( cos α cos β sin θ + sin β sin ϕ cos θ + sin α cos β cos ϕ cos θ ) α ˙ = q ( p cos α + r sin α ) tan β + 1 m V cos β ( ( L + Δ L ) F T sin α ) + 1 V cos β ( g ( sin α sin θ + cos α cos ϕ cos θ ) ) β ˙ = p sin α r cos α + 1 m V ( Y F T cos α sin β ) + 1 V ( g ( cos α sin β sin θ + cos β sin ϕ cos θ sin α sin β cos ϕ cos θ ) )
ϕ ˙ = p + tan θ ( q sin ϕ + r cos ϕ ) θ ˙ = q cos ϕ r sin ϕ ψ ˙ = sec θ ( q sin ϕ + r cos ϕ )
p ˙ = ( c 1 r + c 2 p ) q + c 3 ( L + Δ L ) + c 4 N q ˙ = c 5 p r c 6 ( p 2 r 2 ) + c 7 ( M + Δ M ) r ˙ = ( c 8 p c 2 r ) q + c 4 ( L + Δ L ) + c 9 N
where x E , y E , and z E are the position coordinates in the inertial frame; V represents the velocity in the wind frame; α and β are the aerodynamic angles; ϕ , θ , and ψ denote the attitude angles; p, q, and r are the angular rates; u, v, and w are the velocity components in the body frame, with  u = V cos α cos β , v = V sin β , and w = V sin α cos β ; Δ u , Δ v , and Δ w refer to the components of the wake velocity in the body frame of the trailing UAV, which can be obtained by transforming V B ; L, D, Y and L , M , N represent the aerodynamic forces and moments, whose calculation processes are presented in [28]; c 1 to c 9 are the inertia coefficients, whose expressions are also demonstrated in [28]; m and g are the mass and gravity acceleration, respectively; and Δ L , Δ D and Δ L , Δ M are the induced forces and moments. The control inputs of the model are defined as the thrust F T , the elevator δ e , the aileron δ a , and the rudder δ r .

3. Structure of the First-Order LADRC

In this section, the first-order LADRC is investigated for the close-formation control system. The core idea of the LADRC is to estimate and compensate for various internal and external uncertainties of the control object by regarding them as a total disturbance. The LADRC is composed of two important parts, namely the linear extended-state observer (LESO) and the linear state-error feedback controller (LSEF).
Consider the following form of the first-order nonlinear system:
x ˙ = f ( x , t ) + d ( t ) + b u y = x
where f ( x , t ) is the internal uncertainty, d ( t ) represents the external disturbance, x and y are the state and output of the system, u is the control input, and b is the gain. For the above nonlinear system, the overall control structure of the first-order LADRC is illustrated in Figure 3. Let x 1 = x , x 2 = f ( x , t ) + d ( t ) + ( b b 0 ) u , where b 0 is an estimate of the real b. Regard x 2 as the total disturbance of the system and call it the extended state. By assuming x ˙ 2 = ξ , the original nonlinear system can be transformed into an integral chain system that includes the total disturbance, as follows:
x ˙ 1 = x 2 + b 0 u x ˙ 2 = ξ y = x 1
Then, the design process of the LESO and the LSEF are described for the transformed integral chain system.
LESO. The system state and total disturbance can be estimated using the LESO. A second-order extended-state observer is designed as
ε = z 1 y z ˙ 1 = z 2 + l 1 ε + b 0 u z ˙ 2 = l 2 ε
where z 1 and z 2 are the estimated values of the state x 1 and the total disturbance x 2 , respectively, and l 1 and l 2 are the observer gains to be adjusted. It should be noted that the transformed integral chain system can be observed by the LESO if l 1 and l 2 are selected appropriately. The bandwidth method proposed by Gao is employed to determine the observer gains,
l 1 = 2 ω l 2 = ω 2
where ω is the observer bandwidth.
LSEF. Based on the LESO, the total disturbance can be compensated for using the LSEF. The control signal after disturbance compensation is designed as
u = u 0 z 2 b 0
where u 0 is the control output of the LSEF. By ignoring the estimation error of z 2 x 2 and substituting (15) into (12), one has
x ˙ = u 0 y = x
The transfer relationship from u 0 to y becomes integral. So far, the dynamic linearization for the original nonlinear system is realized by real-time disturbance estimation and compensation. For the integral system, the LSEF is designed as
u 0 = k p ( r * z 1 )
where r * is the reference signal and k p is the feedback gain.
Remark 2. 
The convergence of the LESO and the stability of the closed-loop LADRC have been proven under the assumption that the differentiation of the total disturbance is bounded. Furthermore, the estimation and compensation capacity for the total disturbance depends on the bandwidth of the LESO and the gain of the LSEF, respectively.

4. Robust Control System Design

From the modified model, the flight dynamics of the trailing UAV would be seriously disturbed by the wake vortex effect. Therefore, the robust close-formation control system is essential for enhancing flight stability and acquiring the maximum aerodynamic benefits.

4.1. Control Objective

The robust close-formation control system is implemented on the trailing UAV to steer it to the sweet spot and hold its position. The position of the sweet spot relative to the leading UAV is defined as ( Δ x E d , Δ y E d , Δ z E d ) in the inertial frame. The tracking errors of the relative position components are calculated by
e x E = Δ x E d Δ x E e y E = Δ y E d Δ y E e z E = Δ z E d Δ z E
where Δ x E , Δ y E , and Δ z E are the actual relative position components of the trailing UAV. To acquire the maximum aerodynamic benefit, the control objectives for the three relative position components are specified as lim t e x E = 0 , lim t e y E = 0 , and lim t e z E = 0 . It should be noted that the control error of each relative position component is required to be less than 0.1 b for an engineering application.

4.2. Control System Design

Given the dynamic characteristics of the trailing UAV, the control system is divided into three control subsystems for the longitudinal, altitude, and lateral channels. Considering the strength of the wake vortex effect, the control subsystem of each channel is composed of two cascade controllers, namely an outer-loop position LADRC controller and an inner-loop attitude LADRC controller. The overall scheme of the control system is shown in Figure 4. Although these channels are coupled with each other, the coupling effect can be incorporated into the total disturbance of the channel. Therefore, the control subsystems for these channels are designed independently using the first-order LADRC. It should be noted that the action mechanism of the control system is to receive the relative position components of the sweet spot and translate them into the control inputs for the trailing UAV.
Longitudinal Channel. According to the modified model, the longitudinal position x E is mainly related to the velocity V, which mainly depends on the thrust F T . Therefore, the longitudinal channel can be divided into two cascade subsystems with strict-feedback form: the outer-loop longitudinal position subsystem and the inner-loop velocity subsystem. The outer-loop longitudinal position subsystem is defined by a first-order nonlinear system. The longitudinal position is considered the output and the velocity serves as the virtual input. The inner-loop velocity subsystem, on the other hand, is described by a first-order nonlinear system in which the velocity is the output variable and the input is provided by the thrust. Considering the principle of the first-order LADRC, the longitudinal channel is further transformed into a chain system including two first-order integral subsystems. They are
x ˙ E = f x E + b x E V V ˙ = f V + b V F T
where f x E and f V are the total disturbance of the outer-loop longitudinal position subsystem and the inner-loop velocity subsystem, respectively, and b x E = V c o s ( θ α ) and b V = c o s α m are the subsystem gains. The first-order LADRC controller is designed for each subsystem to dynamically estimate and compensate for the total disturbance, which can effectively suppress the transmission of the disturbance caused by the wake vortex effect. The specific forms are
ε x E = x E z 1 x E z ˙ 1 x E = z 2 x E + l 1 x E ε x E + b x E V * z ˙ 2 x E = l 2 x E ε x E l 1 x E = 2 ω x E l 2 x E = ω x E 2 u 0 x E = k p x E ( x E * z 1 x E ) V * = u 0 x E z 2 x E / b x E ε V = V z 1 V z ˙ 1 V = z 2 V + l 1 V ε V + b V F T z ˙ 2 V = l 2 V ε V l 1 V = 2 ω V l 2 V = ω V 2 u 0 V = k p V ( V * z 1 V ) F T = u 0 V z 2 V / b V
where x E * is the reference signal of the outer-loop longitudinal position LADRC controller, V * is not only the output of the outer-loop longitudinal position LADRC controller but also the reference signal of the inner-loop velocity LADRC controller, F T is the output of the inner-loop velocity LADRC controller, and the meanings of the other parameters are consistent with those in Section 3.
Altitude Channel. Based on the modified model, the altitude z E is mainly determined by the pitching angle θ , which is only controlled by the pitching angle rate q. Therefore, the altitude channel consists of two cascade subsystems with strict-feedback form: the outer-loop altitude subsystem and the inner-loop pitching-angle subsystem. Using the same transformation method as the longitudinal channel, the chain system of the altitude channel is
z ˙ E = f z E + b z E θ θ ˙ = f θ + b θ q
where f z E and f θ are the total disturbance of the outer-loop altitude subsystem and the inner-loop pitching-angle subsystem, respectively, and b z E = V ( π / 180 ) and b θ = 1 are the subsystem gains. Similarly, the first-order LADRC controller of each subsystem is designed as
ε z E = z E z 1 z E z ˙ 1 z E = z 2 z E + l 1 z E ε z E + b z E θ * z ˙ 2 z E = l 2 z E ε z E l 1 z E = 2 ω z E l 2 z E = ω z E 2 u 0 z E = k p z E ( z E * z 1 z E ) θ * = u 0 z E z 2 z E / b z E ε θ = θ z 1 θ z ˙ 1 θ = z 2 θ + l 1 θ ε θ + b θ q * z ˙ 2 θ = l 2 θ ε θ l 1 θ = 2 ω θ l 2 θ = ω θ 2 u 0 θ = k p θ ( θ * z 1 θ ) q * = u 0 θ z 2 θ / b θ
where z E * is the reference signal of the outer-loop altitude LADRC controller, θ * is not only the output of the outer-loop altitude LADRC controller but also the reference signal of the inner-loop pitching-angle LADRC controller, and q * is the output of the inner-loop pitching-angle LADRC controller.
Lateral Channel. Unlike the longitudinal and altitude channel, the lateral channel is further decomposed into two subchannels, namely the rolling subchannel and the yawing subchannel. Note that these two subchannels are closely coupled.
For the rolling subchannel, the lateral position y E is mainly related to the rolling angle ϕ , which is only affected by the rolling-angle rate p. Therefore, the rolling subchannel can be divided into two cascade subsystems with strict-feedback form: the outer-loop lateral-position subsystem and the inner-loop rolling-angle subsystem. Similarly, the chain system of the rolling subchannel is
y ˙ E = f y E + b y E ϕ ϕ ˙ = f ϕ + b ϕ p
where f y E and f ϕ are the total disturbance of the outer-loop lateral-position subsystem and the inner-loop rolling-angle subsystem, respectively, and b y E = V c o s α c o s β c o s θ ( π / 180 ) and b ϕ = 1 are the subsystem gains. The first-order LADRC controller of each subsystem is formulated as
ε y E = y E z 1 y E z ˙ 1 y E = z 2 y E + l 1 y E ε y E + b y E ϕ * z ˙ 2 y E = l 2 y E ε y E l 1 y E = 2 ω y E l 2 y E = ω y E 2 u 0 y E = k p y E ( y E * z 1 y E ) ϕ * = u 0 y E z 2 y E / b y E ε ϕ = ϕ z 1 ϕ z ˙ 1 ϕ = z 2 ϕ + l 1 ϕ ε ϕ + b ϕ p * z ˙ 2 ϕ = l 2 ϕ ε ϕ l 1 ϕ = 2 ω ϕ l 2 ϕ = ω ϕ 2 u 0 ϕ = k p ϕ ( ϕ * z 1 ϕ ) p * = u 0 ϕ z 2 ϕ / b ϕ
where y E * is the reference signal of the outer-loop lateral-position LADRC controller, ϕ * is not only the output of the outer-loop lateral-position LADRC controller but also the reference signal of the inner-loop rolling-angle LADRC controller, and p * is the output of the inner-loop rolling-angle LADRC controller.
The yawing subchannel only includes the yawing-angle system. The chain system of the yawing subchannel is given as
ψ ˙ = f ψ + b ψ r
where ψ is the yawing angle, r is the yawing-angle rate, and f ψ is the total disturbance of the yawing-angle system. Due to the coupling effect of the rolling and yawing subchannels, the yawing is disturbed as rolling occurs. Therefore, the first-order LADRC controller is designed to compensate for this disturbance. The specific form is
ε ψ = ψ z 1 ψ z ˙ 1 ψ = z 2 ψ + l 1 ψ ε ψ + b ψ r * z ˙ 2 ψ = l 2 ψ ε ψ l 1 ψ = 2 ω ψ l 2 ψ = ω ψ 2 u 0 ψ = k p ψ z 1 ψ r * = u 0 ψ z 2 ψ b ψ
where r * is the output of the yawing-angle LADRC controller.
Remark 3. 
According to the dynamic characteristics of the trailing UAV, the attitude rates of the altitude and lateral channels still need to be controlled. Therefore, a simple proportional (P) control is employed to add an attitude-rate control loop inside the attitude control loop. A trial-and-error tuning method is adopted to tune parameter P. The detailed design process of the attitude-rate control loop is omitted.

4.3. Stability Analysis of the Control System

To ensure the stability of the multiloop control system, this subsection mainly analyzes the stability of the entire error-feedback closed-loop control subsystem designed in Section 4.2. Before proceeding to the formal analysis, a theorem about the LESO of the LADRC is introduced.
Theorem 1. 
Assume that the total disturbance f of the system is differentiable and let h = f ˙ . If h is bounded, there exists a constant σ i > 0 and a finite T > 0 such that | x ˜ i ( t ) | σ i , i = 1 , 2 , , n + 1 , t > T > 0 and ω > 0 .
Proof. 
The proof of the theorem was previously established by Zhiqiang Gao [19] and is omitted here for brevity. □
According to Theorem 1, the observers designed in (20), (22), (24), and (26) are all stable. Building on this foundation, the stability analysis of the entire error-feedback closed-loop subsystem is described below. Construct a Lyapunov function of the following form for the multiloop control system of the UAV close-formation model:
V ( x ) = 1 2 e 1 T e 1 + 1 2 e 2 T e 2 + 1 2 e 3 T e 3
where x represents the state variables of the multiloop control system and e 1 , e 2 , and e 3 are the tracking errors of the position loop, attitude loop, and angular-rate loop, respectively. Therefore, the entire multiloop control system is asymptotically stable as long as V ˙ ( x ) < 0 . Next, this is proven by adding each loop incrementally from the inner loop to the outer loop.
Angular-Rate Loop Stability Analysis. The tracking error of the angular-rate loop is e 3 = x 3 c x 3 , where x 3 c = [ P c , Q c , R c ] T is the command signal. For this error, the Lyapunov function V 3 ( x 3 ) is constructed as follows:
V ( x 3 ) = 1 2 e 3 T e 3
To accelerate the response speed, the attitude angular-rate loop adopts the P control. By tuning the appropriate values k 3 = [ k p , k Q , k R ] for the P parameters of each control channel, it can be guaranteed that V ˙ ( x 3 ) < 0 . This implies that the tracking error e 3 converges asymptotically to zero and x 3 converges asymptotically to x 3 c .
Attitude Loop Stability Analysis. The tracking error of the attitude loop is defined as e 2 = x 2 c x ^ 2 , where x 2 = [ V , ϕ , θ , ψ ] is the state variable of the attitude loop and x ^ 2 denotes the estimated values of the state variables from the designed LESO. For the attitude loop, construct the following Lyapunov function
V ( x 2 , x 3 ) = 1 2 e 2 T e 2 + 1 2 e 3 T e 3 = 1 2 e 2 T e 2 + V ( x 3 )
It should be pointed out that because the attitude loop includes the angular-rate loop, its tracking error must be taken into account when designing the Lyapunov function for the attitude loop. Assuming that the observation errors of the states and total disturbances are σ 2 1 and σ 2 2 , then x ^ 2 = x 2 + σ 2 1 and f ^ 2 = f 2 + σ 2 2 . For the error e 3 , the state x 3 of the angular-rate loop converges asymptotically to the command signal x 3 c . Assuming that the convergence error is ε 3 , then u 2 = x 3 c = x 3 + ε 3 , i.e.,  x 3 = u 2 ε 3 . Combining the control variable u 2 of the attitude loop in (20), (22), (24), and (26), the derivative of the Lyapunov function is calculated by
V ˙ ( x 2 , x 3 ) = e 2 T e ˙ 2 + V ˙ ( x 3 ) = e 2 T ( x ˙ 2 c x ˙ 2 σ ˙ 2 1 ) + V ˙ ( x 3 ) = e 2 T ( x ˙ 2 c f 2 b 2 x 2 σ ˙ 2 1 ) + V ˙ ( x 3 ) = e 2 T ( x ˙ 2 c f 2 b 2 ( u 2 ε 3 ) σ ˙ 2 1 ) + V ˙ ( x 3 ) = e 2 T ( x ˙ 2 c b 2 b 2 1 ( k 2 e 2 + x ˙ 2 c f 2 σ 2 2 ) + b 2 ε 3 σ ˙ 2 1 ) + V ˙ ( x 3 ) = e 2 T ( k 2 e 2 + b 2 ε 3 + σ 2 2 σ ˙ 2 1 ) + V ˙ ( x 3 ) e 2 T ( k 2 e 2 b 2 ε 3 + σ 2 2 σ ˙ 2 1 ) + V ˙ ( x 3 )
where b 2 is the control-gain vector of the attitude loop and k 2 is the error-proportional coefficient vector of the attitude angles. It should be noted that due to the existence of V ˙ ( x 3 ) , the error-proportional coefficient of the attitude loop also includes k 3 . Due to the convergence of the LESO and the angular-rate loop, ε 3 , σ 2 2 , and σ ˙ 2 1 are all bounded. Therefore, by tuning the appropriate values for k 2 and k 3 , it can be ensured that V ˙ ( x 2 , x 3 ) < 0 . This implies that the tracking error e 2 asymptotically converges to zero and x ^ 2 asymptotically converges to x 2 c . Furthermore, as  x ^ 2 asymptotically converges to x 2 , the attitude-angle state x 2 asymptotically converges to x 2 c . It is worth noting that the velocity loop of the attitude loop does not include the angular-rate loop but rather takes the thrust of the trailing UAV as input. However, the error analysis process for the velocity loop is similar to that of other attitude-angle variables. Therefore, there is no separate analysis for the velocity loop.
Position Loop Stability Analysis. Similarly, the tracking error of the position loop is e 1 = x 1 c x ^ 1 , where x 1 = [ x E , y E , z E ] . For this error, the Lyapunov function for the position loop is constructed using the same structure as previously mentioned
V ( x 1 , x 2 , x 3 ) = 1 2 e 1 T e 1 + 1 2 e 2 T e 2 + 1 2 e 3 T e 3 = 1 2 e 1 T e 1 + V ( x 2 , x 3 )
It is essential to point out that the Lyapunov function of the position loop is exactly the Lyapunov function of the multiloop control system (27). Assuming that the observation errors of the states and the total disturbances are σ 1 1 and σ 1 2 , then x ^ 1 = x 1 + σ 1 1 and f ^ 1 = f 1 + σ 1 2 . For the error e 2 , the state x 2 of the angular-rate loop converges asymptotically to the command signal x 2 c . Assuming that the convergence error is ε 2 , then u 1 = x 2 c = x 2 + ε 2 , i.e.,  x 2 = u 1 ε 2 . Similar to the analysis process for the attitude loop, the derivative of the Lyapunov function is given by
V ˙ ( x 1 , x 2 , x 3 ) = e 1 T ( k 1 e 1 + b 1 ε 2 + σ 1 2 σ ˙ 1 1 ) + V ˙ ( x 2 , x 3 ) e 1 T ( k 1 e 1 b 1 ε 2 + σ 1 2 σ ˙ 1 1 ) + V ˙ ( x 2 , x 3 )
Therefore, one can draw a similar conclusion that by tuning the appropriate values for k 1 , k 2 , and k 3 , it can be ensured that V ˙ ( x 1 , x 2 , x 3 ) < 0 , which means that the derivative of the Lyapunov function of the multiloop control system for the UAV close-formation model is V ˙ ( x ) < 0 .
In summary, the entire multiloop control system for the UAV close-formation model is asymptotically stable.

5. Sine-Powered Pigeon-Inspired Optimization

In the designed control system, each first-order LADRC controller has two control parameters, namely the feedback gain k p * and the observer bandwidth ω * ( * = x E , V , z E , θ , y E , ϕ , ψ ). Therefore, the control subsystem of the longitudinal and altitude channels contains four control parameters and the control subsystem of the lateral channel includes six control parameters. Due to the cascade structure of the control subsystem, the correlation among the control parameters is very high. To consider all control parameters together and find their optimal values, the tuning of the control parameters is restructured as an optimization problem. Sine-powered controlled PIO is utilized to optimize the control parameters of the control subsystem of each channel.

5.1. Standard PIO Algorithm

The standard PIO algorithm is divided into the map and compass operator and the landmark operator to perform a search based on the homing behavior of pigeons. Using the two operators can achieve better performance in searching for the global optimal position. It is assumed that there are N p i g pigeons in the D p i g -dimension search space. The position and velocity of the i-th pigeon are X i p i g = [ x i 1 p i g , x i 2 p i g , , x i D p i g p i g ] and V i p i g = [ v i 1 p i g , v i 2 p i g , , v i D p i g p i g ] , respectively. They are updated in each iteration.
Map and Compass Operator. In the map and compass operator, the new position and velocity of the i-th pigeon at the N C -th iteration are calculated by
V i p i g ( N C ) = V i p i g ( N C 1 ) e R N C + r a n d · X g b e s t p i g X i p i g ( N C 1 ) X i p i g ( N C ) = X i p i g ( N C 1 ) + V i p i g ( N C )
where R is the map and compass factor; r a n d represents a random number with r a n d [ 0 , 1 ] ; and X g b e s t p i g is the current global optimal position, which can be obtained by comparing the fitness value with all the positions of the pigeons. When the number of iterations reaches N C 1 m a x p i g , stop the current operator and proceed to the landmark operator.
Landmark Operator. In the landmark operator, the pigeon with the lowest fitness value will be abandoned. This decreases the number of pigeons by half in each iteration and the remaining pigeons update their positions according to the center of the pigeons, which can be described as
N p i g ( N C ) = N p i g ( N C 1 ) 2 X c p i g ( N C ) = i = 1 N p i g ( N C ) X i p i g ( N C 1 ) F i t ( X i p i g ( N C 1 ) ) i = 1 N p i g ( N C ) F i t ( X i p i g ( N C 1 ) ) X i p i g ( N C ) = X i p i g ( N C 1 ) + r a n d · ( X c p i g ( N C ) X i p i g ( N C 1 ) )
where X c p i g denotes the center position of all the pigeons. Note that F i t ( X i p i g ( N C 1 ) ) = 1 / f i t ( X i p i g ( N C 1 ) ) , where f i t ( · ) is the fitness function. The whole algorithm ends after another N C 2 m a x p i g iterations.

5.2. SCPIO Algorithm

In the map and compass operator of the standard PIO, the map and compass factor R can balance the global exploration and local exploitation of the whole pigeon group. Global exploration is strengthened if R is relatively large; conversely, local exploitation can be enhanced. However, a fixed R is always adopted in the map and compass operator of the standard PIO. As a result, R has not made a contribution to the balance of global exploration and local exploitation. Therefore, a sine-powered controlled improvement strategy is employed to adjust R dynamically with the iterations. The adjustment formula is written as
R ( N C ) = r ( N C ) · sin ( π R ( N C 1 ) ) r ( N C ) = r m a x N C · ( r m a x r m i n ) / N C m a x
where r is the control parameter in the range of ( 0 , 1 ) ; r m a x and r m i n are the specified maximum and minimum values, respectively; and N C m a x p i g represents the maximum number of iterations. The random nature of R is gradually weakened with the iterations. This indicates that the dynamic adjustment strategy of R retains the advantages of traversing and randomizing the weight in the early iterations, which urges the pigeons to perform extensive global exploration. In the later iterations, the pigeons gather around the global optimal position and the low flight velocity enables the pigeons to perform a fine search. Furthermore, the new position of each pigeon is obtained by simply adding the previous position X i p i g ( N C 1 ) and the current velocity V i p i g ( N C ) in the map and compass operator of the standard PIO. Although the global exploration and local exploitation of the algorithm have been balanced by adjusting R dynamically, this balance may be conservative since the relationship between the position and velocity is not considered. According to [26], the position-update formula with the dynamic weight is given as follows
X i p i g ( N C ) = ω i p i g ( N C ) · X i p i g ( N C 1 ) + ω i p i g ( N C ) · V i p i g ( N C ) + r a n d · σ ( N C ) · X g b e s t p i g
where ω i p i g and ω i p i g are the dynamic weights of the position and velocity of the i-th pigeon, respectively, and σ is the acceleration coefficient. They are calculated by
ω i p i g ( N C ) = σ ( N C ) = e x p f i t X i p i g ( N C ) f i t ( N C ) ¯ 1 + e x p f i t X i p i g ( N C ) f i t ( N C ) ¯ N C ω i p i g ( N C ) = 1 ω i p i g ( N C )
where f i t ( N C ) ¯ is the average fitness value of all pigeons at the N C -th iteration.
Note that the standard PIO that has been improved by using the sine-powered controlled strategy is abbreviated as SCPIO.
Remark 4. 
In the landmark operator of the standard PIO, all pigeons will move toward the center position, which implies that the movement direction of each pigeon is definite. Therefore, there is little room for improvement according to the action mechanism of the landmark operator. Considering this fact, SCPIO still adopts the landmark operator of the standard PIO.

5.3. Construction of the Fitness Function

According to the control objective, the integrated-time absolute error (ITAE) criterion is employed to construct the performance index, namely the fitness function, for the optimal design of the control system. On the one hand, the ITAE criterion includes two factors: time and error, which can take into account both the dynamic performance and stable performance of the control system and simultaneously ensure the response speed and stable accuracy. On the other hand, the ITAE criterion emphasizes the effect of the recent response, which can minimize the influence of the large initial error on the dynamic process. Furthermore, the attitude ITAE criterion of the trailing UAV is also incorporated into the fitness function to reduce the attitude oscillation in the position responses. Since the control system is composed of three control subsystems for different channels, the parameter optimization of the control system is transformed into the parameter optimization of three control subsystems.
Therefore, the fitness functions of three control subsystems for different channels are given as follows
J L o n = k e x E 0 t T t | e x E | d t + k Δ V 0 t T t | Δ V | d t J A l t = k e z E 0 t T t | e z E | d t + k Δ θ 0 t T t | Δ θ | d t J L a t = k e y E 0 t T t | e y E | d t + k Δ ϕ 0 t T t | Δ ϕ | d t + k Δ ψ 0 t T t | Δ ψ | d t
where f i t L o n , f i t A l t , and f i t L a t are the fitness functions of the longitudinal, altitude, and lateral channel, respectively; e x E , e z E , and e y E are the tracking errors between the sweet spot and the actual relative position; Δ V , Δ θ , Δ ϕ , and Δ ψ are the velocity and attitude disturbances; k * ( = e x E , e z E , e y E , Δ V , Δ θ , Δ ϕ , Δ ψ ) represents the weight coefficient; and t T denotes the simulation termination time of the control system. By minimizing the fitness function, the optimal control parameters of the control subsystem can be obtained.

5.4. Optimization Procedure

Since the control subsystem is designed independently according to the channel, the parameter optimization of the control subsystem for the longitudinal channel is taken as an example. Using the SCPIO algorithm, the detailed optimization procedures for the optimal control parameters are described as follows:
Step 1. 
Initialize the SCPIO parameters, including the number of pigeons N p i g , the dimension of thesearch space D p i g , the maximum and minimum values of the map and compass factor r m a x and r m i n , the iteration numbers of two operators N C 1 m a x p i g and N C 2 m a x p i g , and the position X i p i g and velocity V i p i g of all pigeons.
Step 2. 
Drive the close-formation simulation system using the pigeons in Step 1 to calculate the fitness function. Compare the fitness value and find the current optimal position.
Step 3. 
Conduct the iteration. If  N C N C 1 m a x p i g , perform the improved map and compass operator to update the pigeons. Then, drive the close-formation simulation system using the updated pigeons to calculate the fitness function. Update the optimal position by comparing the new fitness values with the current optimal one. When N C 1 m a x p i g < N C , perform the landmark operator to continue the similar optimization process.
Step 4. 
Once the iteration time reaches N C 1 m a x p i g + N C 2 m a x p i g , terminate the algorithm and output the optimal position X g b e s t p i g .
The pseudocode of the above steps is given in Algorithm 1.
Algorithm 1: SCPIO.
Drones 07 00238 i001

6. Simulation Results and Analysis

In this section, the simulation verification is conducted for a leading-trailing UAV close formation based on the F-16 aircraft model that serves two main purposes. The first is to determine the location of the sweet spot for the trailing UAV in close-formation flight. This is achieved by analyzing the variations in the wake vortex effect with the relative position to the leading UAV based on the established wake vortex model. The second is to validate the tracking performance of the designed robust close-formation control system via the SCPIO under the wake vortex effect. Furthermore, the robust close-formation control system based on a PI structure is used for comparison to more intuitively display the advantages of the designed control system. Both the inner and outer loops of the longitudinal and lateral channels adopt the P controller. The PI controller is developed for the outer loop of the altitude channel and the inner loop still adopts the P controller. Similarly, the parameter of each controller for the three channels also needs to be tuned.
In the simulation results, it is assumed that the leading UAV maintains a level and stable flight path at a speed of 152 m/s and an altitude of 4605 m, with attack and pitch angles of 4.5 deg. The initial flight states of the trailing UAV are x E ( 0 ) = 0 m, y E ( 0 ) = 0 m, z E ( 0 ) = 4572 m, V ( 0 ) = 152 m/s, α ( 0 ) = 4.5 deg, β ( 0 ) = 0 deg, ϕ ( 0 ) = 0 deg, θ ( 0 ) = 4.5 deg, ψ ( 0 ) = 0 deg, p ( 0 ) = 0 deg/s, q ( 0 ) = 0 deg/s, and r ( 0 ) = 0 deg/s. The initial relative position between the leading and trailing UAV is set to Δ P ( Δ x E , Δ y E , Δ z E ) = [ 14 b , 2.9 b , 3.6 b ] . The modeling parameters of the wake vortex are specified as V = 152 m/s, μ = 1.6 , C L = 0.2134 , and C L α = 0.0422 , N = 100 . The aerodynamic derivatives, mass properties, and structural parameters of the F-16 aircraft model are given in [27]. The robust close-formation control system is applied to the trailing UAV to steer it to the sweet spot.

6.1. Analysis of the Sweet Spot

According to the continuous-horseshoe vortex method, the intensity of the wake vortex can be described as a function of the position in the body frame of the leading UAV. The variation of the dimensionless velocity field induced by the wake vortex with the longitudinal position is illustrated in Figure 5. We can see that the intensity of the induced velocity field decreases significantly with the negative increase in the longitudinal position. Furthermore, the intensity decays almost to zero when x < 10 b , which implies that the wake vortex effect disappears completely. Therefore, the downstream longitudinal position should not exceed 10 b if one wants to utilize the wake vortex effect of the leading UAV. The position of the sweet spot relative to the leading UAV is denoted as ( Δ x B d , Δ y B d , Δ z B d ) in the body frame of the leading UAV. Considering the fuselage length of the UAV and the minimum safe spacing of close formation, the longitudinal relative position component Δ x B d of the sweet spot is set at 3 b to maximize the wake vortex effect. The variation of the induced velocity with the lateral and vertical positions calculated at x = 3 b is demonstrated in Figure 6. Two sweet spots A s 1 and A s 2 can be seen and their positions are consistent with earlier findings [6]. At the sweet spot, the maximum induced upwash velocity occurs and the induced sidewash velocity is relatively small. Therefore, the lateral and vertical relative position components are Δ y B d = ± 0.75 b and Δ z B d = 0 , respectively. The trailing UAV should be driven toward the sweet spot and maintained in that position during close-formation flight. Note that the region designated by B is the area with the maximum downwash velocity, where the trailing UAV should be forbidden to appear. This means that the robust close-formation control system needs to have high control accuracy for the trailing UAV. Note that the position of the sweet spot is analyzed by the separation distance relative to the leading UAV, which makes the position independent of the coordinate frame. Therefore, the relative position ( Δ x E d , Δ y E d , Δ z E d ) of the sweet spot in the inertial frame is still ( 3 b , ± 0.75 b , 0 ) .

6.2. Implementation of Control System Optimization

Since the tracking and anti-disturbance performance of the designed robust control system depends on the control parameters of the control subsystem of each channel, the proposed SCPIO is used to search for their optimal values rather than a manual search. Then, the optimization results are compared to those of the standard PIO [22] and particle swarm optimization (PSO) [29] to validate the proposed SCPIO, which ensures the optimization of the control parameters. The parameters of each optimization algorithm are presented in Table 1. Considering the amplitude range of the wake vortex effect and the response characteristics of the outer-loop position and the inner-loop attitude, the search sets of the observer bandwidths of the position and attitude loops are given as { ω * ( = x E , y E , z E ) | 0 < ω * < 1 } and { ω * ( = V , θ , ϕ , ψ ) | 5 < ω * < 10 } , respectively. According to the expected dynamic characteristics of the control system, the controller gains are set to k p * ( = x E , z E ) [ 0.055 , 0.085 ] , k p y E [ 0.01 , 0.04 ] , and k p * ( = θ , ϕ , ψ ) [ 0.7 , 1.4 ] , k p V [ 0.1 , 0.4 ] , respectively. Note that the optimization of the control parameters is performed according to the divided channel. Furthermore, only the wake vortex effect of the corresponding channel is added to the UAV dynamic model in the optimization. The position and attitude weight coefficients of each channel are given as k * ( = e x E , Δ V ) = 1 , k * ( = e z E , Δ θ ) = 1 and k e y E = 1 , k Δ ϕ = 10 , and k Δ ψ = 100 . The simulation termination time of the control system is set to t T = 400 s. Each optimization algorithm is run 10 times and the optimization result with the smallest fitness value was selected as the result of the optimization algorithm. The optimal values are presented in Table 2 and the evolution curves are illustrated in Figure 7. The proposed SCPIO has the strongest optimization ability, indicating that control parameters optimized with SCPIO can lead to the best performance of the control system in terms of tracking accuracy and rejection ability. By observing the optimal values of each optimization algorithm, a common phenomenon is apparent, that is, the bandwidths of the inner- and outer-loop extended-state observers of the altitude and lateral channels are larger than the bandwidth of the longitudinal channel. This indicates that the wake vortex effect is quite strong in these two channels.
In the robust close-formation control system based on a PI structure, the proposed SCPIO is also utilized to search for the optimal values of the control parameters to ensure that the comparison results are more fair and convincing. The fitness function of each channel and the parameters of the algorithm are the same as those of the designed control system. Similar to the above analysis and optimization process, the obtained optimal values are presented in Table 3. Note that the fitness value of the lateral channel at the optimal position is less than those of the other two channels. This is because the optimal values of the control parameters can make the position and attitude accurately converge to the given command and initial state, respectively.

6.3. Tracking-Performance Validation

In order to validate the tracking performance of the designed robust close-formation control system under the wake vortex effect, two qualitative criteria are introduced, dynamic performance and stable performance. Furthermore, the tracking performance of the robust close-formation control system based on a PI structure under the wake vortex effect is also discussed for comparative analysis. On the other hand, the reference-tracking responses of the two control systems are generated without considering the wake vortex effect. Note that in the following demonstration of the simulation results, SCPIO-LADRC represents the proposed design and SCPIO-PI denotes the PI structure. SCPIO-LADRC Ref and SCPIO-PI Ref are the respective reference tracking responses.
The relative position tracking responses of the trailing UAV with respect to the leading UAV are demonstrated in Figure 8. The proposed design and the PI structure have almost the same reference-tracking responses. However, the dynamic performance of the PI structure deteriorates considerably when the wake vortex effect is added, indicating that the PI structure has almost no anti-disturbance ability. Even worse, the relative position not only significantly deviates from the sweet spot but also shows no sign of convergence at the end of the simulation. Compared to the PI structure, there is no significant difference between the tracking responses of the proposed design and its reference responses. This demonstrates the strong ability of the proposed design to suppress the disturbances caused by the wake vortex effect. Furthermore, it should also be noted that the dynamic performance of the proposed design deteriorates in terms of the setting time. For the expected control error of ± 0.1 b , however, the setting time increases slightly compared to that of the reference response. From the perspective of stable performance, the proposed design can guarantee that the relative position accurately converges to the sweet spot. On the other hand, the relative position components of the disturbed altitude and lateral channels under the PI structure have a large amplitude oscillation around the sweet spot. Therefore, these two channels need a high frequency for the extended-state observer, which further validates the optimization results of SCPIO.
The inner-loop state responses of the trailing UAV, including the velocity V, aerodynamic angles α , β , attitude angles ϕ , θ , ψ , and angular rates p , q , r , are illustrated in Figure 9. Similar to the results of the tracking responses, dynamic performance under the PI structure cannot be guaranteed. This is because the disturbances caused by the wake vortex effect can be divided into the induced forces acting on the outer-loop position subsystem and the induced moments acting on the inner-loop attitude subsystem. Owing to the induced moments, there are serious high-frequency oscillations and spikes in the inner-loop attitude subsystem. Furthermore, it should be noted that the poor inner-loop dynamic performance is a leading cause of the significant deterioration of the outer-loop dynamic performance. In contrast, the inner-loop state responses of the trailing UAV based on the proposed design exhibit smooth and stable convergence behavior. This means that the influences of the inner-loop disturbances on the dynamic performances are not transmitted to the outer loop. In conclusion, the LADRC controller is also used in the inner loop to form a cascade control system with the outer-loop LADRC controller, which can realize hierarchical suppression for the strong external disturbances caused by the wake vortex effect and significantly improve the tracking performance of the outer loop.
According to the established wake vortex model, the wake vortex effect can be considered a function of the relative position and attitude between the leading and trailing UAVs. If the trailing UAV has different tracking trajectories to converge to the sweet spot, the time histories of the wake vortex effect are also different, as shown in Figure 10. Under the control of the PI structure, the trailing UAV experienced more aggressive variations in the wake vortex effect. This is due to the fact the PI structure control has poor dynamic performance after being disturbed by the wake vortex effect, which makes the trailing UAV oscillate considerably around the sweet spot. The oscillations make the trailing UAV periodically fly through the downwash and upwash wake vortex regions, thereby inducing the more aggressive wake vortex effect. In turn, the more aggressive wake vortex effect can cause further deterioration in the control performance. The time responses of the control inputs under the two different control systems are shown in Figure 11. For the PI structure, severe control inputs are generated to counteract the disturbance of the more aggressive wake vortex effect as much as possible, which can cause actuator faults. However, through the proposed design, the disturbance of the wake vortex effect can be suppressed with less control effort. From the perspective of the control inputs, therefore, the proposed design is much more efficient than the PI structure. Furthermore, the proposed design leads to a 21.59 % decrease in the thrust input at the sweet spot, indicating that the trailing UAV could potentially save about 21.59 % of energy during close-formation flight.

7. Conclusions

In this paper, a robust close-formation control system with dynamic estimation and compensation is designed for UAV close-formation flight. The designed control system is divided into three control subsystems for the longitudinal, altitude, and lateral channels. The control subsystem of each channel is composed of two cascaded first-order LADRC controllers. One is responsible for the outer-loop position control and the other is used to stabilize the inner-loop attitude. This control system scheme can significantly reduce the coupling effect between channels and effectively suppress the transmission of the disturbance caused by the wake-vortex effect. A continuous-horseshoe vortex method with high estimation accuracy is employed to estimate the wake-vortex effect. The wake-vortex effect decreases significantly with the negative increase in the longitudinal position. Therefore, the longitudinal separation distance between the leading and trailing UAV should not exceed 10 b if one wants to utilize the wake-vortex effect. The estimated wake-vortex effect is integrated into a nonlinear high-fidelity UAV model to describe the dynamic characteristics of the trailing UAV under the disturbance of the wake-vortex effect. SCPIO is utilized to simultaneously optimize the control parameters for the control subsystem of each channel, which can help the subsystem to achieve optimal performance. Compared to the conventional PI structure, the designed robust close-formation control system achieves error-bounded stable and robust dynamic performance. Furthermore, the control system follows the wake-vortex model with high estimation accuracy and the nonlinear high-fidelity UAV model, making it more suitable for engineering applications. In the future, engineering verification will be conducted to assess its practicality.

Author Contributions

Author Contributions: Conceptualization, G.Y. and H.D.; methodology, G.Y. and H.D.; software, G.Y.; validation, G.Y. and H.D.; formal analysis, G.Y.; investigation, G.Y.; resources, G.Y.; data curation, G.Y.; writing—original draft preparation, G.Y.; writing—review and editing, H.D.; visualization, G.Y.; supervision, H.D.; project administration, H.D.; funding acquisition, H.D. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Natural Science Foundation of China grant numbers T2121003, 91948204, and U20B2071.

Institutional Review Board Statement

Not applicable for studies not involving humans or animals.

Informed Consent Statement

This study did not involve humans.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Zhan, G.; Gong, Z.; Lv, Q.; Zhou, Z.; Wang, Z.; Yang, Z.; Zhou, D. Flight Test of Autonomous Formation Management for Multiple Fixed-Wing UAVs Based on Missile Parallel Method. Drones 2022, 6, 99. [Google Scholar] [CrossRef]
  2. Li, S.; Li, Y.; Zhu, J.; Liu, B. Predefined Location Formation: Keeping Control for UAV Clusters Based on Monte Carlo Strategy. Drones 2022, 7, 29. [Google Scholar] [CrossRef]
  3. Xu, D.; Guo, Y.; Yu, Z.; Wang, Z.; Lan, R.; Zhao, R.; Xie, X.; Long, H. PPO-Exp: Keeping Fixed-Wing UAV Formation with Deep Reinforcement Learning. Drones 2022, 7, 28. [Google Scholar] [CrossRef]
  4. Jiang, Y.; Bai, T.; Wang, Y. Formation Control Algorithm of Multi-UAVs Based on Alliance. Drones 2022, 6, 431. [Google Scholar] [CrossRef]
  5. Luo, Y.; Bai, A.; Zhang, H. Distributed formation control of UAVs for circumnavigating a moving target in three-dimensional space. Guid. Navig. Control 2021, 1, 2150014. [Google Scholar] [CrossRef]
  6. Zhang, Q.; Liu, H.H. Aerodynamics Modeling and Analysis of Close Formation Flight. J. Aircr. 2017, 54, 2192–2204. [Google Scholar] [CrossRef]
  7. Kent, T.E.; Richards, A.G. Analytic Approach to Optimal Routing for Commercial Formation Flight. J. Guid. Control Dyn. 2015, 38, 1872–1884. [Google Scholar] [CrossRef]
  8. Bangash, Z.A.; Sanchez, R.P.; Ahmed, A.; Khan, M.J. Aerodynamics of Formation Flight. J. Aircr. 2006, 43, 907–912. [Google Scholar] [CrossRef]
  9. Hanson, C.E.; Pahle, J.; Reynolds, J.R.; Andrade, S.; Nelson, B. Experimental Measurements of Fuel Savings During Aircraft Wake Surfing. In Proceedings of the 2018 Atmospheric Flight Mechanics Conference, Atlanta, GA, USA, 25–29 June 2018. [Google Scholar]
  10. Zhang, Q.; Liu, H.H. Aerodynamic Model-Based Robust Adaptive Control for Close Formation Flight. Aerosp. Sci. Technol. 2018, 79, 5–16. [Google Scholar] [CrossRef]
  11. Zhang, Q.; Liu, H.H. UDE-Based Robust Command Filtered Backstepping Control for Close Formation Flight. IEEE Trans. Ind. Electron. 2018, 65, 8818–8827. [Google Scholar] [CrossRef]
  12. Galzi, D.; Shtessel, Y. Closed-Coupled Formation Flight Control Using Quasi-Continuous High-Order Sliding-Mode. In Proceedings of the 2007 American Control Conference, New York, NY, USA, 9–13 July 2007. [Google Scholar]
  13. Liu, C.; Jiang, B.; Zhang, K. Adaptive Fault-Tolerant H-Infinity Output Feedback Control for Lead–Wing Close Formation Flight. IEEE Trans. Syst. Man Cybern. Syst. 2020, 50, 2804–2814. [Google Scholar] [CrossRef]
  14. Yuan, G.; Xia, J.; Duan, H. A Continuous Modeling Method via Improved Pigeon-Inspired Optimization for Wake vVortices in UAVs Close Formation Flight. Aerosp. Sci. Technol. 2022, 120, 107259. [Google Scholar] [CrossRef]
  15. Han, J. From PID to Active Disturbance Rejection Control. IEEE Trans. Ind. Electron. 2009, 56, 900–906. [Google Scholar] [CrossRef]
  16. Zhang, Y.; Chen, Z.; Zhang, X.; Sun, Q.; Sun, M. A Novel Control Scheme for Quadrotor UAV Based upon Active Disturbance Rejection Control. Aerosp. Sci. Technol. 2018, 79, 601–609. [Google Scholar] [CrossRef]
  17. Ahi, B.; Nobakhti, A. Hardware Implementation of an ADRC Controller on a Gimbal Mechanism. IEEE Trans. Control Syst. Technol. 2018, 26, 2268–2275. [Google Scholar] [CrossRef]
  18. Gao, Z. Scaling and Bandwidth-Parameterization Based Controller Tuning. In Proceedings of the American Control Conference, Minneapolis, MN, USA, July 2003; Available online: https://www.semanticscholar.org/paper/Scaling-and-bandwidth-parameterization-based-tuning-Gao/4f4d7b650767e01138a6969b97e5e2601779fb4c (accessed on 20 March 2023).
  19. Zheng, Q.; Gaol, L.Q.; Gao, Z. On Stability Analysis of Active Disturbance Rejection Control for Nonlinear Time-Varying Plants with Unknown Dynamics. In Proceedings of the 46th IEEE Conference on Decision and Control, New Orleans, LA, USA, 12–14 December 2007. [Google Scholar]
  20. Qi, Y.; Liu, J.; Yu, J. Dynamic Modeling and Hybrid Fireworks Algorithm-Based Path Planning of an Amphibious Robot. Guid. Navig. Control 2022, 2, 2250002. [Google Scholar] [CrossRef]
  21. Zhu, K.; Han, B.; Zhang, T. Multi-UAV Distributed Collaborative Coverage for TargetSearch Using Heuristic Strategy. Guid. Navig. Control 2021, 1, 2150002. [Google Scholar] [CrossRef]
  22. Duan, H.; Qiao, P. Pigeon-Inspired Optimization: A New Swarm Intelligence Optimizer for Air Robot Path Planning. Int. J. Intell. Comput. 2014, 7, 24–37. [Google Scholar] [CrossRef]
  23. Tong, B.; Wei, C.; Shi, Y. Fractional Order Darwinian Pigeon-Inspired Optimization for Multi-UAV Swarm Controller. Guid. Navig. Control 2022, 2, 2250010. [Google Scholar] [CrossRef]
  24. Duan, H.; Zhao, J.; Deng, Y.; Shi, Y.; Ding, X. Dynamic Discrete Pigeon-Inspired Optimization for Multi-UAV Cooperative Search-Attack Mission Planning. IEEE Trans. Aerosp. Electron. Syst. 2021, 57, 706–720. [Google Scholar] [CrossRef]
  25. Huo, M.; Duan, H.; Fan, Y. Pigeon-Inspired Circular Formation Control for Multi-UAV System with Limited Target Information. Guid. Navig. Control 2021, 1, 2150004. [Google Scholar] [CrossRef]
  26. Chen, K.; Zhou, F.; Liu, A. Chaotic Dynamic Weight Particle Swarm Optimization for Numerical Function Optimization. Knowl. Based Syst. 2018, 139, 23–40. [Google Scholar] [CrossRef]
  27. Sonneveldt, L. Nonlinear F-16 Model Description; Delft University of Technology: Delft, The Netherlands, 2006. [Google Scholar]
  28. Richard, S.R. Nonlinear F-16 Simulation Using Simulink and Matlab; University of Minnesota: Minneapolis, MN, USA, 2003. [Google Scholar]
  29. Poli, R.; Kennedy, J.; Blackwell, T. Particle Swarm Optimization: An Overview. Swarm Intell. 2007, 1, 33–57. [Google Scholar] [CrossRef]
Figure 1. UAV close-formation flight.
Figure 1. UAV close-formation flight.
Drones 07 00238 g001
Figure 2. Composition of the wake vortex.
Figure 2. Composition of the wake vortex.
Drones 07 00238 g002
Figure 3. General structure block of the first-order LADRC.
Figure 3. General structure block of the first-order LADRC.
Drones 07 00238 g003
Figure 4. Overall scheme of the control system.
Figure 4. Overall scheme of the control system.
Drones 07 00238 g004
Figure 5. Sectional views of dimensionless velocity field induced by the wake vortex at different longitudinal positions.
Figure 5. Sectional views of dimensionless velocity field induced by the wake vortex at different longitudinal positions.
Drones 07 00238 g005
Figure 6. Variation of the induced velocity with lateral and vertical position (x = −3b).
Figure 6. Variation of the induced velocity with lateral and vertical position (x = −3b).
Drones 07 00238 g006
Figure 7. Evolution curves of each optimization algorithm.
Figure 7. Evolution curves of each optimization algorithm.
Drones 07 00238 g007
Figure 8. Relative position tracking responses of the trailing UAV with respect to the leading UAV.
Figure 8. Relative position tracking responses of the trailing UAV with respect to the leading UAV.
Drones 07 00238 g008
Figure 9. Velocity, aerodynamic angles, attitude angles, and attitude rate responses of the trailing UAV.
Figure 9. Velocity, aerodynamic angles, attitude angles, and attitude rate responses of the trailing UAV.
Drones 07 00238 g009
Figure 10. Time histories of the wake vortex effect experienced by the trailing UAV.
Figure 10. Time histories of the wake vortex effect experienced by the trailing UAV.
Drones 07 00238 g010
Figure 11. Time responses of the control inputs.
Figure 11. Time responses of the control inputs.
Drones 07 00238 g011
Table 1. Parameters of each optimization algorithm.
Table 1. Parameters of each optimization algorithm.
AlgorithmParameterDescriptionValue
PSO N C m a x Maximum iterative number50
N P S O Number of particles100
ω P S O Inertia weight0.4
c 1 Self-learning factor2
c 2 Group-learning factor2
PIO, SCPIO N C 1 m a x p i g Iteration number of the map and compass operator30
N C 2 m a x p i g Iteration number of the landmark operator20
N p i g Number of pigeons100
RMap and compass factor0.4
[ r m i n , r m a x ] Range of the control parameter of the sine map[0.1, 0.9]
Table 2. Optimal values of each optimization algorithm (The SCPIO algorithm and its optimized results are highlighted in bold).
Table 2. Optimal values of each optimization algorithm (The SCPIO algorithm and its optimized results are highlighted in bold).
ChannelControl ParametersAlgorithmOptimal ValuesFitness Value
Longitudinal [ K p x E , ω x E , K p V , ω V ] SCPIO[0.0712, 0.11, 0.26, 6.12]60,126
PIO[0.0634, 0.16, 0.23, 6.56]72,151
PSO[0.0607, 0.20, 0.21, 6.81]78,136
Altitude [ K p z E , ω z E , K p θ , ω θ ] SCPIO[0.06854, 0.71, 1.05, 8.21]139,842
PIO[0.07345, 0.64, 1.16, 7.78]199,774
PSO[0.07562, 0.61, 1.24, 7.46]239,729
Lateral [ K p y E , ω y E , K p ϕ , ω ϕ , K p ψ , ω ψ ] SCPIO[0.0254, 0.27, 0.98, 7.04, 1.10, 7.58]93,251
PIO[0.0207, 0.30, 0.87, 7.43, 1.21, 7.42]134,696
PSO[0.0318, 0.25, 1.25, 6.78, 1.00, 7.63]113,974
Table 3. Optimal values for the robust close-formation control system based on a PI structure.
Table 3. Optimal values for the robust close-formation control system based on a PI structure.
ChannelControl ParametersAlgorithmOptimal ValuesFitness Value
Longitudinal [ K p x E P I , K p V P I ] SCPIO[0.1024, 250.7534]28,146
Altitude [ K p z E P I , K i z E P I , K p θ P I ] SCPIO[0.0021, 0.0005, 1.1568]39,084
Lateral [ K p y E P I , K p ϕ P I , K p ψ P I ] SCPIO[0.08791, 1.1326, 2.9736]10,211
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Yuan, G.; Duan, H. Robust Control for UAV Close Formation Using LADRC via Sine-Powered Pigeon-Inspired Optimization. Drones 2023, 7, 238. https://doi.org/10.3390/drones7040238

AMA Style

Yuan G, Duan H. Robust Control for UAV Close Formation Using LADRC via Sine-Powered Pigeon-Inspired Optimization. Drones. 2023; 7(4):238. https://doi.org/10.3390/drones7040238

Chicago/Turabian Style

Yuan, Guangsong, and Haibin Duan. 2023. "Robust Control for UAV Close Formation Using LADRC via Sine-Powered Pigeon-Inspired Optimization" Drones 7, no. 4: 238. https://doi.org/10.3390/drones7040238

APA Style

Yuan, G., & Duan, H. (2023). Robust Control for UAV Close Formation Using LADRC via Sine-Powered Pigeon-Inspired Optimization. Drones, 7(4), 238. https://doi.org/10.3390/drones7040238

Article Metrics

Back to TopTop