Modeling and Flight Experiments for Swarms of High Dynamic UAVs: A Stochastic Configuration Control System with Multiplicative Noises

UAV Swarm with high dynamic configuration at a large scale requires a high-precision mathematical model to fully exploit its boundary performance. In order to instruct the engineering application with high confidence, uncertainties induced from either systematic measurement or the environment cannot be ignored. This paper investigates the Ito^ stochastic model of the UAV Swarm system with multiplicative noises. By combining the cooperative kinematic model with a simplified individual dynamic model of fixed-wing-aircraft for the first time, the configuration control model is derived. Considering the uncertainties in actual flight, multiplicative noises are introduced to complete the Ito^ stochastic model. Following that, the estimator and controller are designed to control the formation. The mean-square uniform boundedness condition of the proposed stochastic system is presented for the closed-loop system. In the simulation, the stochastic robustness analysis and design (SRAD) method is used to optimize the properties of the formation. More importantly, the effectiveness of the proposed model is also verified using real data of five unmanned aircrafts collected in outfield formation flight experiments.


Introduction
Swarms of UAVs, which can autonomously implement missions [1], have received significant attention in recent years. There are many application scenarios for UAV swarms, such as comprehensive combat [1], distributed reconfigurable sensor networks [2], surveillance [3], and reconnaissance systems [4]. The primary concern of the UAV swarm is the configuration control problem and related research mainly focuses on mathematical modeling [5], control strategies and methods [6,7] and collision and obstacle avoidance algorithms [8,9]. However, most of them are based on a deterministic system or a system with ideal Gaussian noises. High dynamic USCC model with multiplicative noise remains as one of the primary and practical issues for utilizing UAV swarms in engineering applications.
Generally, the stochastic differential equation refers to a stochastic process-driven system or an ordinary differential equation with a random coefficient [10][11][12][13][14]. Random factors are introduced into the system in the following three ways [10]: 1. the system's initial conditions or inputs are taken as random variables, 2. the system's random external disturbances, 3. the system's parameters and structures taken as random variables. The disturbances in the UAV swarm system are mainly induced Sensors 2019, 19, 3278 2 of 24 from the sensor measurement, internet transmission and the task environment, and they fall under point 2 mentioned above.
Existing work on the formation configuration control problem has been extensively investigated for low dynamic deterministic systems. In [15], the formation containment problems based on linear state equations of the multiagent systems were investigated. In [6,16], some theoretical and practical problems of multiple quadrotors were studied. In [17][18][19], the dynamical model based formation control problems of multirobot systems were studied. Most of the studies focus on low-speed vehicles, such as multiple quadrotors or multirobot systems, with an ideal environment; however, very few studies consider the formation control problems of high dynamic fixed-wing unmanned aircraft under complex environments. The mathematical modeling of stochastic high dynamic UAV swarms remains to be solved.
In this paper, we introduce the stochastic model for the following two reasons: (1) The UAV swarm, especially for high dynamic, dense configuration and large scale swarms, requires a high-precision mathematical model to describe the dynamic relationships among formation members and fully exploit its boundary performance; (2) When the UAV swarm is carrying out missions, the influence of various uncertainties (systematic measurement random interferences, network-induced random interferences and mission environment random interferences) cannot be ignored, and the relative movements of its members are usually random. Communication topology, caused by a network change, inevitably influences the process of information sending and receiving [20]. Therefore, it is necessary to combine the mathematical model of the USCC with the stochastic system to instruct the engineering practice such as cooperative detection under complex mission environments with higher confidence.
Undoubtedly, constructing a more adaptable stochastic model for multiple vehicle systems is an urgent task. The problem of Brownian motion-driven multiagent tracking was discussed in [21] and sufficient conditions for the tracking of multi-agents were obtained by using the auxiliary function of Brownian motion and random Itô integral technology. A time lag multiagent system model with measurement noise was set up in [22], and the stability theory of stochastic differential equations was used. In [23], the stochastic factor has been considered in the leader-following multiagent model based on the event-triggered control strategy, Itô formula and stability theory. Studies [20,[24][25][26][27] have considered the influence of stochastic disturbances on the multiagent system (MAS) and have established a stochastic model to control the MAS. However, existing works about the stochastic MAS are mainly based on assumed state matrices. Although the assumed formula can be applied across multiple levels, extra difficulties will occur in a certain practical system; for example, the process of combining the flight member's dynamics model and the formation's kinematics model is more complex; the modeling of process noises is more complex because the coefficients are presented in a very complex way in the practical system; the simulation and flight test are more difficult because of the complex system and high-risk environment. Therefore, numerous problems remain to be solved.
Typically, there are four main methods for formation coordination modeling reported in the literature: the leader-follower framework [28], virtual structure approach [29,30], behavior-based model [31,32] and graph theories [33]. Most of them focus on the consensus problems based on the kinematic model. In this paper, intending to improve both the formation properties and individual capabilities of UAV swarm in complex environments, we first use a simplified nonlinear dynamics model of the fixed-wing aircraft flight control system and then construct the group dynamic cooperative control system model together with the relative kinematic model. Based on Reynolds' three criteria [31], the model comprehensively considers the flight member's individual properties and the whole swarm's cooperativeness.
To the best of our knowledge, although few efforts have been devoted to the modeling of dynamic formation systems, there is almost no literature on the stochastic model of the UAV swarms. For example, a six-degree-of-freedom (DOF) unified nonlinear dynamic model of spacecraft formation was presented in [34]; In [35,36], the formation control laws for YF-22 aircraft models with six DOF dynamics plus kinematic equations were designed. Although these formation models are efficient for the cooperative of a group without stochastic disturbance, the control strategy may be invalid when control objects are moving in the noise environment.
To the best of the authors' knowledge, few papers discuss the quality of the configuration, and most of them have come up with an algorithm to control the formation, thereby achieving the desired configuration or avoiding collision [6][7][8]. However, they do not discuss the robustness of UAV swarms in much detail. In order to improve the robustness of the stochastic system, the parameters of the estimator and the controller are optimized by the stochastic robustness analysis and design (SRAD) method [37] in the simulation. Furthermore, few studies have carried out fixed-wing flight test experiments. A multi-UAV outfield flight experiment was carried out to verify the effectiveness of the formation collision forecast and coordination algorithm in [9]. In [36], a set of flights was performed to assess the performance of the formation control laws. To extend the previous outfield flight test results, the overall design in this paper is validated experimentally by flight testing using the leader-follower configuration.
Motivated by the discussions above, the stochastic USCC model with multiplicative noises is investigated in this paper. Compared with the existing literature, the main theoretical and experimental contributions of our work are summarized as follows: (a) Firstly, with the aim to instruct the engineering application of UAV swarms, we construct the nonlinear formation model by combining the nonholonomic individual dynamics model of a UAV swarm with the relative cooperative kinematics model. The model comprehensively considers the personality of the individual members and the cooperation of the whole formation. (b) Secondly, the stochastic state equation and output equation, together with the estimator and controller, finally constitute a state-feedback-based closed-loop Itô stochastic system, which makes full use of the platform's boundary performance and better matches the complex task environments to improve the cost-effectiveness of the UAV swarm system. (c) Thirdly, most studies of formation configuration control focus on theoretical analysis, while the technology proposed in this paper is for engineering applications and has been verified by outfield flight tests.
The rest of this paper is structured as follows: In Section 2, the problem's formulation and preliminary studies of the USCC stochastic system are presented. In Section 3, the mathematical model of formation control stochastic system is illustrated in detail. The estimator and controller are also designed to control the formation, and SRAD has been used to optimize the controller and estimator. The mean-square uniformly bounded condition of the proposed stochastic system is then presented. In Section 4, simulations and experiments are conducted to verify the effectiveness of the model. Finally, concluding remarks are given in Section 5.

Itô Stochastic System
Consider the linear Itô stochastic system as the model to be investigated as follows: Wiener processes, which are defined in the complete probability space (Ω, F, P) and are independent of each other. Define the following matrices: where N = n 2 , '⊕' denotes the Kronecker tensor for a matrix, and A ⊕ A = I n ⊗ A + A ⊗ I n , '⊗' is the Kronecker tensor product of a matrix.
where the element '1' appears in the first column, (n + 1)th column and [(n − 1)n + 1]th column of K N .

Mean-Square Uniform Boundedness of the Stochastic System
Since the stochastic system is complicated by external interferences, its stability condition is relatively strict and there is no trivial solution to the equation. To solve this problem, we take advantage of the mean-square uniform boundedness of the stochastic system. The condition for boundedness is slightly less strict than that of stability. Under the condition of boundedness, the states of the system are converging to bounded areas instead of certain stable points as time tends to infinity. Therefore, for the USCC stochastic system model, stability refers to the mean-square uniform boundedness.

Definition 1.
If there is a positive number c: Then the states X ij (t, t 0 , X ij0 ) are mean-square bounded. The subscript '0' represents the initial value, X ij0 denotes the initial states of the system which is composed of two members: i and j, t 0 denotes the initial time, t denotes the current time, · is the Euclidean norm, E{·} is the mathematical expectation, and sup is the minimum upper bound. Lemma 1. [12] The necessary and sufficient condition for the mean-square boundedness of the solution for the time-varying linear stochastic system (1) is that the following time-varying linear deterministic system is bounded: Lemma 3. [12] If B(t), G i (t), (i = 1, 2, · · · , m) are bounded and system (2) is mean-square uniformly asymptotically stable, then the solution of system (1) is mean-square uniformly bounded.
Based on Lemma 1, Lemma 2, Lemma 3, we present sufficient conditions for the stability of the Itô stochastic model and prove them. Theorem 1. The sufficient conditions for the stability of the Itô stochastic model (or the mean-square uniform boundedness of the stochastic system (1)) are: (1) B(t) and G i (t) are bounded.
(2) The real part of the eigenvalues of matrix M(t) are negative. (1), although we can use the conclusion of Lemma 1 to obtain the necessary and sufficient conditions for the mean-square boundedness directly, given that A and F k , k = 1 ∼ m are linear time-invariant matrices, we further simplify the mean-square boundedness conditions. According to Lemma 3, if it satisfies condition (1), the mean-square uniform boundedness of the system (1) is equivalent to system (2) and is mean-square uniformly asymptotically stable:

Proof. For Equation
According to Lemma 2, system (2) is mean-square uniformly asymptotically equivalent to condition (2). Proof completed.

Estimator of Itô Stochastic System
Lemma 4. [38] considering the Itô stochastic system in the form as follows: where x(t) ∈ R n , y(t) ∈ R m are system states and measured values, respectively.A, A 0 , H and H 0 are constant matrices (they can also be extended to time-varying matrices if needed). w(t) is a standard scalar Wiener process, as well as w 1 (t) and w 2 (t), where w 1 (t) ∈ R n and w 2 (t) ∈ R m . The initial state x(0) is a zero mean second-order stochastic process.
Assuming that x(0) is independent of w(t), w 1 (t), and w 2 (t), and it satisfies: Then the linear estimator with minimum mean square error is: x(0) = 0 The gain of the estimator is: where P(t) can be obtained as follow: Remark 1. In this study, according to the above lemmas, we can construct a more applicable and effective USCC stochastic model to improve the formation properties of the UAV swarm system. Moreover, the measurement equation with Gaussian noises, the optimal estimator and controller are ingeniously involved in the closed-loop model. The mean square uniform boundedness condition of USCC stochastic system can be obtained based on the lemmas and theorem proposed above.

USCC Stochastic System Modeling
The group dynamic cooperative control system model comprehensively considers the personality of the individual members and the cooperativeness of the whole formation based on Reynolds' three criteria [31]. Generally, the model is built with virtual forces: individuality and interaction forces. Individuality force describes the nodes' individual characteristics. Interaction force indicates the quality of group collaboration among nodes and describes the group dynamic cooperative characteristics, reflecting the ability to obey Reynolds' three criteria. We will use it as a general theory to guide the modeling of the USCC stochastic system.

Model of Individual Flight Control System
The individual flight control system of the UAV swarm adopts the north-up-east coordinate system. Assumption 1. The formation moves in a two-dimensional plane, thus the flight path inclination and pitch velocity are zero; the aircraft adopts side slip turning, thus the speed inclination angle, roll angle, roll angle velocity, angle of attack and side slip angle are all small values.

Assumption 2.
Thrust P is independent of velocity V.
The simplified nonlinear mathematical model of individual flight control system is: where V is the flight velocity; ϕ is the flight path declination; β is the lateral slip angle; ω is the rotational angle velocity of the body's coordinate system relative to the ground coordinate system. The subscript "y" denotes the y-component of ω. J y is the y-component of inertial moments of the body's coordinate system. M y is the y-component of moment caused by the external force (including thrust) on the mass center; P is the thrust; X is the resistance force; Z is the lateral force; δ is the rudder declination; K δ , T δ are gain and time constants of the control surface response, respectively (subscripts x,y, z are aileron, rudder, and elevator, respectively); K p , T p are gain and time constants for the thrust response, respectively; δ c , δ Pc are rudder angle command and thrust command, respectively. By performing a small-disturbance linearization on (16) [39], we can obtain: where X V = ∂X ∂V , the same as other elements which is in the same form with X V in (17).

Model of Formation Control System
For the convenience of outfield experiments and the safety of UAV swarm, we construct the model in a two-dimensional (2D) plane to make it more adaptive to complex task environments such as flat and dense formations, under which the aircraft would carry out missions at low altitude with almost no vertical maneuver space. Moreover, theoretical results can be covered fully and extended to the three-dimensional (3D) space. Therefore, we focus on the problem of USCC in the two-dimensional plane, as shown in Figure 1.
where V X X V ∂ = ∂ , the same as other elements which is in the same form with V X in (17).

Model of Formation Control System
For the convenience of outfield experiments and the safety of UAV swarm, we construct the model in a two-dimensional (2D) plane to make it more adaptive to complex task environments such as flat and dense formations, under which the aircraft would carry out missions at low altitude with almost no vertical maneuver space. Moreover, theoretical results can be covered fully and extended to the three-dimensional (3D) space. Therefore, we focus on the problem of USCC in the two-dimensional plane, as shown in Figure 1. In Figure 1, i v and j v are the two nodes adjacent to each other. The flight path coordinate of node i v is set as the relative coordinate system, in which the x-axis represents the direction of the velocity and the z-axis is perpendicular to the x-axis shown as in Figure 1.
The ground coordinate system is set as the fixed coordinate system. ij d is the distance between the two nodes; ij x and ij z are the relative distances in the forward and lateral directions of the flight path coordinate system, respectively. i V , j V , i ϕ and j ϕ represent their velocities and flight path declinations in the ground coordinate system, respectively. With the relationship from theoretical mechanics: absolute velocity = relative velocity + convected velocity, the following kinematics equation for node i v and node j v can be derived: The above equation can be decomposed in the flight path coordinate system of the node i v as: In Figure 1, v i and v j are the two nodes adjacent to each other. The flight path coordinate of node v i is set as the relative coordinate system, in which the x-axis represents the direction of the velocity and the z-axis is perpendicular to the x-axis shown as in Figure 1. The ground coordinate system is set as the fixed coordinate system. d ij is the distance between the two nodes;x ij and z ij are the relative distances in the forward and lateral directions of the flight path coordinate system, respectively. V i , V j , ϕ i and ϕ j represent their velocities and flight path declinations in the ground coordinate system, respectively.
With the relationship from theoretical mechanics: absolute velocity = relative velocity + convected velocity, the following kinematics equation for node v i and node v j can be derived: The above equation can be decomposed in the flight path coordinate system of the node v i as: Substituting the second formula in (17) into (19) and performing small perturbation linearization yields: Note that the aircraft flight momentum m i V i is relatively large. Moreover, β i has small values according to Assumption 1, and ϕ j − ϕ i ≈ 0, thus Ignoring these second-order small quantities and simplifying (20), we can get: Combining Equations (17) with (19) and (21) yields the formation control system model: Note that: the state coefficients V i , P i , x ij , z ij and Z i V i in (22) are obtained at the balanced point.

Process Noises
The formation could be easily affected by various forces in the atmosphere that cannot be accurately measured in advance. Therefore, the process noises cannot be ignored.
For the node ν i , assuming that the mass m i , velocity V j and flight path declination angle ϕ j of the adjacent node ν j , which are obtained from the supporting network, are given values (i.e., consider that V j and ϕ j are already estimated in ν j , and ignore the random transmission interference), but x ij , z ij , V i , iy are determined by the aircraft's instantaneous state (such as velocity, altitude, attack angle, yaw angle, etc.); These states and their influences are random in the real flight environment. Therefore, based on the central limit theorem, we assume that the above parameters approximately obey the Gaussian distribution, that is: The subscript "b" denotes that the values are determined and they are obtained at the balanced point. The formation states do not change much when they fly around the balanced point; thus, the variances of the random variables are approximately constant, and it can be assumed that the above parameters are independent of each other. In the following, we use n 1 , n 2 , · · · , n 12 to represent the random variables in (23). (2) (3) (6) in (23) into a 1 yields: Assume that w V i is relatively small compared to the aircraft speed V ib . Since V ib + w V i is in the denominator, the impact of w V i on a 1 is small and can be ignored; assuming that both w Z i V i and w z ij are small,w Z i V i w z ij is a second-order small quantity and can be ignored. Then we get: Assuming that all the random parts of (23) are small, and ignoring the second-order small quantity. For the same reason as a 1 , the expression of the coefficients a 2 ∼ a 15 could be derived. The results are given as follows: β i iy J iy m i V ib n 10 a 10b + a 10b4 n 4 + a 10b7 n 7 + a 10b9 n 9 + a 10b10 n 10 a 12b + a 12b8 n 8 + a 12b10 n 10 + a 12b12 n 12 The states of the system are X ij = [ ∆x ij ∆z ij ∆V i ∆ϕ i ∆β i ∆ω iy ∆P i ∆δ iy ] T ; the inputs are U ij = [ ∆δ iPc ∆δ iyc ∆V j ∆ϕ j ] T ; ∆V j and ∆ϕ j are determined random inputs. Substituting (25) to (36) into (22), then the open-loop state equation of the formation stochastic system could be get: . where For convenience of description, we define the square matrix T mn whose order is eight and has only one nonzero element; the value of the nonzero element is '1' and it lies in column m and row n.

Measurement Noises
The states of the stochastic system (37) are measured by the support network and relative navigation in the UAV swarm's information acquisition system. The measuring vector is defined as: T , assuming that ∆P i cannot be measured.
These measured values are mainly affected by sensor measurement, network transmission and random disturbances in the flight environment. Assuming that the measured noises of the system approximately obey the Gaussian distribution, whose mathematical expectation is 0 and variance is σ m 2 , the measurement equation is: The subscript "m" means the measured value of the system, n 13 , n 15 , · · · , n 19 are standard Gaussian white noises independent of each other. They are also independent of n 1 , n 2 , · · · , n 12 .
In summary, the measurement equation is: where,

Formation Estimator Design
For the state estimating problem of the formation control stochastic system (37) and (39), we use a novel Itô stochastic system estimator which has a fixed gain according to [38].
In Lemma 4, the gain of the estimator is time-varying, but in this paper, the USCC problem is investigated during cruising and the formation maintains a certain configuration with certain speed and height in that period, the formation does not change very much. Therefore, we use an estimator with fixed gain in (10) to estimate states, which can significantly simplify the computation and improve the real-time performance of the system.
The formation control stochastic system (37) and (39) can be described as the following Itô stochastic system: where W k (t) ( k = 1 ∼ 12) in (40)  The initial states X ij (0) satisfy: Then, substituting (40) into (41) and the fixed gain estimator of the formation control stochastic system can be obtained by Lemma 4: T is the state estimate vector; K f is fixed gain of the estimator. U ij is the control input.

Formation Controller Design
It can be seen from (22) that the system states have a high degree of coupling between each other (such as the forward distance ∆x ij and the sideslip angle ∆β i , the lateral distance ∆z ij and ∆V i , they all have high degree of coupling between each other), so the forward and lateral distance will be controlled with a couple in this paper. The PID-based formation control system structure is a commonly used design method in the engineering application at present [40]. According to the clustering algorithm proposed in [41], the PID formation controller we adopted is: K ωij = diag(ω ij , ω ij , 1, 1, 1, 1, 1, 1) (47) whereX ij is the output of the estimator; X * ij is the system command; superscript " * " indicates the command, the same as below; U jd is the determined interference input vector of the adjacent node; K c ∈ R 4×8 is the control law, in which the last two rows in K c are zero vectors because ∆V j and ∆ϕ j are the determined interference inputs in U ij ; K ωij is the adjacency adjustment matrix, 0 ≤ ω ij ≤ 1 are adjacency coefficients. The larger it is, the stronger the adjacency relationship between node ν i and ν j .
The key to designing a better PID controller is to contrive the proper PID gain parameters. In order to get a better parameter of the feedback coefficients, we use the SRAD method which combines the genetic algorithm and Monte Carlo simulation to improve the robustness of the stochastic system.
We can see from (45) that K c∆v i (∆V i − ∆V * i ) and K c∆ϕ i (∆φ i − ∆ϕ * i ) reflect the individuality forces of the individual aircraft, K c∆x ij ω ij (∆x ij − ∆x * ij ) and K c∆z ij ω ij (∆ẑ ij − ∆z * ij ) reflecting the interaction forces which represent the quality of the formation cooperation. Both of them can contribute to maintaining the configuration during formation maneuvering.

Closed-Loop USCC Stochastic System
In summary, the Itô stochastic system of the USCC is as follows: where the first equation is the state equation, the second is the measurement equation, the third is the control input, and the fourth is the state estimation equation. The above four equations together construct the expansion closed-loop equation for the stochastic system of USCC (as shown in Figure 2): ) is the bounded variance vector of measuring noises. Therefore, we can further deduce the following: Proposition 1. Under normal circumstances, a sufficient condition for the stability (or the mean square uniform boundedness of the stochastic system (50)) of the USCC system is that the real parts of the eigenvalues are negative, that is:

Simulation and Experiments
With the aim of cooperative detection under complex environments, we ran simulations with two aircraft to verify the effectiveness of the proposed stochastic model. Moreover, the equivalent outfield flight test was carried out to complete the mission of cooperative detection in a certain area. The results demonstrate that the formation could be achieved effectively and accurately.

Initial State
To make it convenient for analysis, we set the number of UAV swarm  The states of the expansion system are: where X ij ∈ R 16×1 ; X ij ∈ R 8×1 are original states;X ij ∈ R 8×1 are estimated states. The state transfer matrix of the expansion system is: where A ij ∈ R 16×16 ; A ij ∈ R 8×8 is the state transfer matrix of the original system; B ij ∈ R 8×4 is the input matrix of the original system; H ij ∈ R 7×8 is the estimate matrix of the original system; K ωij ∈ R 8×8 is the adjacent adjustment matrix in U ij ; K c ∈ R 4×8 is the control law; K f ∈ R 8×7 is the gain of the estimator. The input matrix of the expansion system is: where B ij (t) ∈ R 16×1 ; X * ij ∈ R 8×1 is the system command; U jd ∈ R 4×1 is the determined input of the adjacent node in U ij . Note that: B ij (t) is a time-varying matrix because X * ij is time-varying. The expansion stochastic state transfer matrix is: The stochastic input matrix of the expansion system is: The standard Wiener process is: where W 1 ∼ W 19 are the Wiener processes in (49); W is an independent 19-dimension standard Wiener process defined in the complete probability space: (Ω, F, P).

Main Results
Define the following matrix: According to Theorem 1, the sufficient conditions of the stability of USCC stochastic model are: (1) B ij (t), G ijk , k = 1 ∼ 19 are bounded.
(2) The real part of the eigenvalues of matrix M ij are negative.
It can be seen from (53) and (55) that all the elements in the matrices B ij (t), G ijk are bounded according to their definition, because in B ij , K ip , T ip are constant parameters of the thrust response and K iδ y , T iδ y are constant parameters of the elevator response; K c ∈ R 4×8 is the control law which can be obtained after simulation; K ωij is an adjacency adjustment matrix whose elements are in the range of (0, 1); X * ij are bounded command values; U jd is the determined interference input vector of the adjacent node; K f ∈ R 8×7 is the fixed gain of the estimator which can be determined from simulation; the nonzero elements in the matrix H ij ∈ R 7×8 are '1'; and in G ijk , E ijk ( k = 13 ∼ 19) is the bounded variance vector of measuring noises. Therefore, we can further deduce the following: Proposition 1. Under normal circumstances, a sufficient condition for the stability (or the mean square uniform boundedness of the stochastic system (50)) of the USCC system is that the real parts of the eigenvalues are negative, that is: where the max Reλ(M ij ) is the maximum real part of the eigenvalue of M ij .

Simulation and Experiments
With the aim of cooperative detection under complex environments, we ran simulations with two aircraft to verify the effectiveness of the proposed stochastic model. Moreover, the equivalent outfield flight test was carried out to complete the mission of cooperative detection in a certain area. The results demonstrate that the formation could be achieved effectively and accurately.

Initial State
To make it convenient for analysis, we set the number of UAV swarm n = 2. The configuration of the formation during cruising is x 12b = 100 m, z 12b = −173.2 m (that is, ν 2 located 100 meters forward and 173.2 meters left of ν 1 ), and the cruising speeds are V 1b = V 2b = 100m/s. The cruising trajectory declination is ϕ 1b = ϕ 2b = 0 rad. The current configuration of the formation is the cruising formation, the current speed is V 1 = V 2 = 100 m/s, and the current flight path declination is ϕ 1 = ϕ 2 = 0 rad. The original mass is m 1 = m 2 = 1400 Kg, the y-components of inertial moments are I y1 = I y2 = 3980 Kg · m 2 .

Formation Parameters
Assume that the supporting network is strongly connected and that the adjacency coefficient ω ij in (47) is 1.

Formation Commands
The expected configuration of the system are x * ij = 50 m, z * ij = −73.2 m, the expected formation speed is V f = 100 m/s, and the expected formation declination is ϕ f = 0 rad.

Optimal Design of Estimator Gain and Control Law
In order to improve the robustness of the stochastic system, the estimator gain in (43) and the control law in (45) are optimized by the SRAD [37].
The SRAD design flow is shown as in Figure 3, which is composed of two parts: a modern optimization algorithm and a control structure design. MCE denotes Monte Carlo evaluation, SRA denotes stochastic robustness analysis.
The cost function we designed is J ij = 12 i=1 w i I i 2 (q i ), where q i represents the 12 indicators which are shown in Table 1, w i is the weight of each indicator, and I i (·) is the membership function of each indicator which obeys the rising-ridge distribution (59) or 0-1 distribution (60).   For indicators 1-9 and 12, the membership function obeys the rising-ridge distribution, i.e., where a, b are the best and allowable value of the indicators respectively. x is the simulation result. For indicators 10 and 11, the membership function obeys the 0-1 distribution, i.e., where (a, b) is the allowable range of the indicators, x is the simulation result. p i (K c , K f ) is the probability of indicator i that cannot satisfy the stability and requirements whose probability distribution function is I i (·).
The minimum cost function J reflects the minimum probability that any indicator cannot satisfy the stability and requirements. It also indicates the minimum errors of the selected properties. Therefore, the obtained controller has high-quality robustness, and the probability that the control system does not meet the requirements is significantly reduced after multiple simulations.
The design steps in Figure 3 are as follows: (1) Design the controller G c (K c ) and estimator G f (K f ) for the controlled object H(n i ); (2) Define the indicators for SRAD: I(H(n i ), G c (K c ), G f (K f )); (3) Carry out Monte Carlo simulation on the closed-loop system to obtain the probabilityp i (K c , K f ) that cannot satisfy the stability and performance; (4) Constitute a random cost functionĴ(K c , K f ) to satisfy both of robust stability and performance; (5) Apply a modern optimization algorithm to get an optimal value. After we get the minimum value ofĴ(K c , K f ), then we obtain a stochastic robust optimal controller and an optimal estimator.
The optimization process of the designed parameters K f and K c is shown in Figure 4. After iterating 15 times with the genetic algorithm and running the Monte Carlo simulation 100 times per iteration, we got the optimal cost value: J = 4.56 and the optimal parameters:     The framework of the simulation model is shown in Figure 5.

Simulation Framework
The framework of the simulation model is shown in Figure 5.
In the framework shown in Figure 5, the flight members in the formation are referred to nodes in the supporting network. The node obtains information through the formation support network and the sensor system, including neighbor nodes' information and environment information. The decision management system allocates missions to flight members and plans the flight route for the formation. Finally, formation control system and member flight control system carry out missions using the Itô stochastic system model (49) and the parameters presented above. The framework of the simulation model is shown in Figure 5.

Simulation Results
In order to better reflect the performance of the USCC stochastic system, we define the following indicator.

Definition 2.
The weighted variance of the estimate error X ij = X ij −X ij is: where the diagonal matrix W v ∈ R 8×8 is the weighted matrix of weighted variances, and the diagonal elements correspond to the statesX ij of the estimator. The larger the weight is, the more accurate the estimation of the corresponding state becomes. The trace of W v is tr(W v ) = 1. Assuming that the estimation accuracies of ∆x ij , ∆ẑ ij , ∆V i and ∆φ i are required to be higher, the weighting matrix we designed is: W v = diag(0.2, 0.2, 0.2, 0.2, 0.1, 0.05, 0, 0.05).
The simulation results correspond with the optimal cost value and the optimal parameters are shown as in Table 1 and Figure 6. Note that the fluctuation in the table refers to the standard deviation of the difference between the real-time configuration and the expected configuration.
It can be observed from Table 1 that all the indicators are in the appropriate range. The real part of the maximum eigenvalue is negative and the weighted variance of the estimation error calculated from simulation results is 3.86. From Figure 6 we can observe that the two aircraft achieved the desired configuration after 54.7 s, and the forward distance steady error is 0.0332 m, the lateral distance steady error is 0.0944 m. The formation speed and the formation declination meet the designed requirements well. Thus the system is mean-square uniform bounded according to Proposition 1. Assuming that the estimation accuracies of and ˆi ϕ Δ are required to be higher, the weighting matrix we designed is: The simulation results correspond with the optimal cost value and the optimal parameters are shown as in Table 1 and Figure 6. Note that the fluctuation in the table refers to the standard deviation of the difference between the real-time configuration and the expected configuration.

Autonomous Flight Experiments
In order to observe the performance of the model and verify the effectiveness of the USCC stochastic system, we conducted an equivalent outfield autonomous formation flight test by using multiple UAVs. As shown in Figure 7, we carried out experiments with seven nonholonomic UAVs. The UAV swarm can cooperatively search the certain area with different configurations. It can be observed from Table 1 that all the indicators are in the appropriate range. The real part of the maximum eigenvalue is negative and the weighted variance of the estimation error calculated from simulation results is 3.86. From Figure 6 we can observe that the two aircraft achieved the desired configuration after 54.7 s, and the forward distance steady error is 0.0332 m, the lateral distance steady error is 0.0944 m. The formation speed and the formation declination meet the designed requirements well. Thus the system is mean-square uniform bounded according to Proposition 1.

Autonomous Flight Experiments
In order to observe the performance of the model and verify the effectiveness of the USCC stochastic system, we conducted an equivalent outfield autonomous formation flight test by using multiple UAVs. As shown in Figure 7, we carried out experiments with seven nonholonomic UAVs. The UAV swarm can cooperatively search the certain area with different configurations. The loads in each cabin of the UAV are shown in Figure 8, including the power module, formation cooperative guidance module, autopilot module, detection module, formation communication module and flight data transmission module. In the experiment, we adopted the proposed USCC stochastic model into the formation cooperative guidance module to instruct the flight members to reach the desired position and maintain a steady configuration. The framework of the whole system is shown in Figure 9. The loads in each cabin of the UAV are shown in Figure 8, including the power module, formation cooperative guidance module, autopilot module, detection module, formation communication module and flight data transmission module. In the experiment, we adopted the proposed USCC stochastic model into the formation cooperative guidance module to instruct the flight members to reach the desired position and maintain a steady configuration. The framework of the whole system is shown in Figure 9.
The loads in each cabin of the UAV are shown in Figure 8, including the power module, formation cooperative guidance module, autopilot module, detection module, formation communication module and flight data transmission module. In the experiment, we adopted the proposed USCC stochastic model into the formation cooperative guidance module to instruct the flight members to reach the desired position and maintain a steady configuration. The framework of the whole system is shown in Figure 9.    The formation ground station is set to observe the real-time formation flight process, and upload commands to instruct the formation to change its configuration or adjust relative distances. It can also control the pod to capture the configuration. The flight member's digital monitor station is set to observe the real-time flight member's flight statuses and to ensure the safety of the whole flight process.
Two configurations of five drones we designed in the experiments are shown in Figure 10a,b. The lateral and forward distance between neighbor UAVs was set to 50 m. Moreover, the configurations in Figure 10c,d were captured by the UAV that was flying higher with a pod.
(a) (b) Figure 9. The framework of the USCC outfield flight system.
The formation ground station is set to observe the real-time formation flight process, and upload commands to instruct the formation to change its configuration or adjust relative distances. It can also control the pod to capture the configuration. The flight member's digital monitor station is set to observe the real-time flight member's flight statuses and to ensure the safety of the whole flight process.
Two configurations of five drones we designed in the experiments are shown in Figure 10a,b. The lateral and forward distance between neighbor UAVs was set to 50 m. Moreover, the configurations in Figure 10c,d were captured by the UAV that was flying higher with a pod.
It can be observed from Figure 10 that the UAVs achieved the desired configuration smoothly and maintained the formation effectively under the proposed model.
For the convenience of analysis, the actual flight data could be saved through the formation monitoring station. In this paper, we take the flight data of the wedge configuration with and without the USCC stochastic system model in the experiment to invert the flight process and evaluate the effectiveness of the USCC stochastic system. The results are shown in Figure 11.
The curves shown as in Figure 11 demonstrate that the five UAVs achieved the desired configuration and steadily maintained the formation. The flight paths in the rectangle that the arrows "1" and "2" point to in Figure 11a are amplified in Figure 11c,d. The data are summarized in Table 2, from which we can observe that in the flight test not using the model, the relative height of the UAVs should be maintained at the very least at 50 m to avoid the risk of collision, while the heights of the five UAVs with USCC stochastic system model converged to 200 m and the relative distance in height was zero. The average 3D distance between flight members in the flight test with the proposed model was shortened by 32.14% compared to the test without the model.
The formation ground station is set to observe the real-time formation flight process, and upload commands to instruct the formation to change its configuration or adjust relative distances. It can also control the pod to capture the configuration. The flight member's digital monitor station is set to observe the real-time flight member's flight statuses and to ensure the safety of the whole flight process.
Two configurations of five drones we designed in the experiments are shown in Figure 10a,b. The lateral and forward distance between neighbor UAVs was set to 50 m. Moreover, the configurations in Figure 10c,d were captured by the UAV that was flying higher with a pod.  It can be observed from Figure 10 that the UAVs achieved the desired configuration smoothly and maintained the formation effectively under the proposed model.
For the convenience of analysis, the actual flight data could be saved through the formation monitoring station. In this paper, we take the flight data of the wedge configuration with and without the USCC stochastic system model in the experiment to invert the flight process and evaluate the effectiveness of the USCC stochastic system. The results are shown in Figure 11.   The curves shown as in Figure 11 demonstrate that the five UAVs achieved the desired configuration and steadily maintained the formation. The flight paths in the rectangle that the arrows "1" and "2" point to in Figure 11a are amplified in Figure 11d and Figure 11c. The data are summarized in Table 2, from which we can observe that in the flight test not using the model, the relative height of the UAVs should be maintained at the very least at 50 m to avoid the risk of collision, while the heights of the five UAVs with USCC stochastic system model converged to 200 m and the relative distance in height was zero. The average 3D distance between flight members in the flight test with the proposed model was shortened by 32.14% compared to the test without the model.  The results show that the introduction of multiplicative noises improves the formation maintenance performance and extends the boundary properties of the formation flight, such as the safety distance is shorter and the relative height can be eliminated. Therefore, the proposed stochastic model could provide the UAV swarm with a larger maneuvering space, and then improve the efficiency and quality of mission execution, enhance the operational capability in high-risk environments and improve the adaptation of the system to the complex environment.

Conclusions
In this paper, the problem of the state estimation and control of the UAV swarm system with the consideration of multiplicative noises is studied. The closed-loop Itô stochastic system we constructed is the combination of a state equation introduced from group kinematic model and individual dynamic model considering the multiplicative noises, an observation equation considers the measurement noises, an estimator and a controller. Following that, the proof to verify the mean-square uniform boundedness of the system is presented. The optimal estimator and controller are obtained with the use of SRAD in the simulation. Finally, simulation results show that the system is stable and the selected indicators meet the requirements. The outfield experiment results demonstrate that the configuration with the proposed model could be significantly condensed by 32.14% compared to the test with the traditional model. Therefore, the stochastic system of USCC with multiplicative noises proposed in this paper could contribute to effectively exploiting the boundary performances of the system and constructing a high dynamic formation in practical application, thus better matching the actual environment.
However, there is much to be researched further in this area. For example, in the practical application, the time delay cannot be ignored, especially for large scale formations. The modeling of the USCC stochastic system considering time delay is currently under investigation.