Distributed model prediction control method for unmanned aerial vehicle formation under discontinuous communication condition

Abstract

The invention discloses a distributed model prediction control method for an unmanned aerial vehicle formation under a non-continuous communication condition. The method comprises the following steps: establishing an unmanned aerial vehicle kinematics model and a trajectory tracking error model; simulating an independent random data loss condition of communication between unmanned aerial vehicles, and establishing a navigator and follower formation control target based on a navigation-following formation mode; the invention relates to a distributed model prediction controller design of an unmanned aerial vehicle formation under a non-continuous communication condition. A navigator unmanned aerial vehicle trajectory tracking controller, a follower unmanned aerial vehicle cooperative controller and a terminal feedback controller are included. Establishing a predicted value of a follower unmanned aerial vehicle cooperative controller for a reference trajectory of a navigator unmanned aerial vehicle and a predicted value of a reference trajectory of a neighbor node unmanned aerial vehicle in a prediction time domain; and solving the rolling optimization problem of each unmanned aerial vehicle at the current moment to obtain an optimal solution. According to the invention, under the condition of non-continuous communication, when data loss is caused, cooperative control can still be carried out among the unmanned aerial vehicles.

Description

Technical Field

[0001] The present application relates to the field of communication technologies, and particularly to a distributed model predictive control method for unmanned aerial vehicle (UAV) formation in the case of discontinuous communication. Background Art

[0002] When performing large-scale search and reconnaissance missions, UAVs often need to cooperate with multiple UAVs to complete the operation tasks. Compared with single-UAV operations, multi-UAV cooperative operations can improve the operation efficiency, expand the operation range, and enhance the search and reconnaissance accuracy, with stronger flexibility. Currently, there are mainly three structures for multi-UAV cooperative control: centralized, decentralized, and distributed. The centralized control structure requires the central node to communicate with all nodes, with a large amount of communication data, high requirements for communication bandwidth and processor computing performance, and poor scalability and anti-interference ability of the entire system. The decentralized control structure cannot handle the coupled cooperative objectives between multiple UAVs, and the system cannot achieve stable control performance. In the distributed control structure, each node only needs to exchange information with neighboring nodes, with a small amount of communication data and high timeliness. Through the independent calculations of each node, large-scale computing problems are dispersed into multiple small-scale computing problems, with a small computing burden on the system and good scalability and stability. Therefore, using the distributed control structure to achieve multi-UAV cooperative control is a more reliable method.

[0003] Currently, most of the methods for realizing multi-UAV distributed cooperative trajectory tracking control do not comprehensively consider problems such as nonlinear characteristics in the system, optimization of control effects, constraints in practical applications, and communication data loss between nodes, and there are many deficiencies compared with actual application scenarios. Summary of the Invention

[0004] In order to solve the problem of cooperative control of existing UAV formations in the case of discontinuous communication, the present application provides a distributed model predictive control method for UAV formations in the case of discontinuous communication, to achieve cooperative trajectory tracking control under the leader-follower formation of multiple UAVs when data loss occurs in the case of discontinuous communication.

[0005] The present application discloses a distributed model predictive control method for UAV formations in the case of discontinuous communication, which includes:

[0006] Step 1: Establish a UAV kinematic model and a trajectory tracking error model;

[0007] Step 2: Simulate the independent random data loss situation of communication between UAVs, and based on the leader-follower formation method, establish the formation control objectives of the leader and the followers.

[0008] Step 3: Design of the distributed model predictive controller for the UAV formation in the case of discontinuous communication, including the trajectory tracking controller of the leader UAV, the cooperative controller of the follower UAVs, and the terminal feedback controller;

[0009] Step 4: Establish the predicted values of the reference trajectory of the leader UAV and the predicted values of the reference trajectories of the neighbor node UAVs by the cooperative controller of the follower UAVs within the prediction time domain;

[0010] Step 5: Solve the rolling optimization problem of each UAV at the current moment to obtain the optimal solution.

[0011] Furthermore, in the said Step 1, the establishment process of the UAV motion model is as follows:

[0012] Ignoring the rolling motion of the UAV, establish the kinematic model of UAV i in the formation system as follows:

[0013]

[0014] In the above formula, x i , y i and z i respectively represent the three-dimensional coordinate information of UAV i in the inertial coordinate system, θ i and ψ i represent the pitch angle and heading angle of UAV i, v i represents the speed of UAV i, ω i and r i represent the pitch angular velocity and heading angular velocity of UAV i; x i = [x i y i z i θ i ψ i T represents the state quantity of UAV i, u i = [v i ω i r i T represents the control quantity of UAV i; the control quantity of UAV i during motion should be kept within the constraint range Ω i during motion, that is, u i ∈ Ω i , where Ω i is expressed as:

[0015] v min ≤ v i ≤ v max

[0016] ω min ≤ ω i ≤ ω max

[0017] r min ≤ r i ≤ r max (2)

[0018] In the above formula, v min and v max are respectively the minimum and maximum values of v i , ω min and ω max are respectively the minimum and maximum values of ω i , r min and r max are respectively the minimum and maximum values of r i ;

[0019] Discretize the mathematical model of formula (1). Take the sampling time as T, and the following discrete kinematic model is obtained:

[0020]

[0021] All variables in the above formula are the discrete forms of the variables in formula (1).

[0022] Furthermore, in the said step 1, the construction process of the trajectory tracking error model is as follows:

[0023] First, establish a reference trajectory model. The reference trajectory of the UAV is determined by a virtual reference UAV with the same model, that is:

[0024]

[0025] In the above formula, x d , y d and z d respectively represent the three-dimensional coordinate information of the reference trajectory in the inertial coordinate system, θ d and ψ d represent the pitch angle and heading angle of the reference trajectory, v d represents the reference trajectory speed, ω d and r d represent the pitch angular velocity and heading angular velocity of the reference trajectory, and their discrete forms are the same as those in formula (3);

[0026] Project the three-dimensional coordinate trajectory tracking error in the inertial coordinate system Oxyz onto the vehicle coordinate system O b x b y b z b . It mainly includes two steps, including:

[0027] 1) Rotate ψ around the z-axis of the ground coordinate system to the coordinate system {G1}. At this time, the x1 axis is located in the plane formed by the z-axis and the x b axis;

[0028] 2) Rotate by θ about the y-axis of the coordinate system G1 to the coordinate system {B}. At this time, it coincides with the vehicle coordinate system O b x b x b y b z b coincide;

[0029] In the case of not considering roll, the coordinate transformation matrix from the inertial coordinate system to the body coordinate system is:

[0030]

[0031] In the above formula, ψ is the heading angle, θ is the pitch angle, R z (ψ) and R y (θ) represent the coordinate transformation matrices for rotation about the z-axis and y b axis, respectively, and are

[0032] The three-dimensional coordinate trajectory tracking error in the inertial coordinate system Oxyz is projected onto the body coordinate system O b x b y b z b expressed as:

[0033]

[0034] Combining the pitch angle θ and the heading angle ψ trajectory tracking errors, the UAV trajectory tracking error is described as:

[0035]

[0036] In the above formula, x ie , y ie and z ie are the three-dimensional coordinate errors between the UAV i and the reference trajectory in the body coordinate system, θ ie and ψ ie are the pitch angle and heading angle errors between the UAV i and the reference trajectory;

[0037] After differentiating equation (7), we get:

[0038]

[0039] Linearize the non-linear trajectory tracking error model in equation (8) at the reference point to get:

[0040]

[0041] Define the control quantity as u ie = [u1 u2 u3] T, where \(u_1=(v i -v d )\cos2\theta d , and

[0042] Express equation (8) in the form of a general linear system and obtain after discretization:

[0043]

[0044] In the above equation T is the sampling time, I is the identity matrix with the same dimension as , The forms of A and B are as follows:

[0045]

[0046] Furthermore, in step 2, the leader-follower formation cooperative control objective, independent random data loss;

[0047] The leader-follower formation cooperative control objective includes the leader trajectory tracking control objective and the follower cooperative control objective;

[0048] b) Leader trajectory tracking control objective:

[0049] Define \(p l (t k ) = [x l (t k ) y l (t k ) z l (t k )] T to represent the position quantity of the leader UAV l, \(q l (t k ) = [\theta l (t k ) \psi l (t k )] T to represent the angular quantity of the UAV l, \(u l (t k ) = [v l (t k ) \omega l (t k ) r l (t k )] T to represent the control quantity of the UAV l, \(p d (t k ) = [x d (t k ) y d (tk )z d (t k )] T represents the position quantity of the reference trajectory, q d (t k )=[θ d (t k )ψ d (t k )] T represents the angular quantity of the reference trajectory, u d (t k )=[v d (t k )ω d (t k )r d (t k )] T represents the control quantity of the reference trajectory. Then, the control objective of the trajectory tracking problem of the leader UAV l is described as:

[0050] When t k →∞, p l (t k )-p d (t k )=0, q l (t k )-q d (t k )=0, u l (t k )-u d (t k )=0;

[0051] b) Follower cooperative control objective:

[0052] Define the set of neighbor node UAVs j of the follower UAV i as represents the desired position quantity of the leader UAV and the follower UAV i, q i (t k )=[θ i (t k )ψ i (t k )] T represents the angular quantity of the follower UAV i, p i (t k )=[x i (t k )y i (t k )z i (t k )] T represents the position quantity of the follower UAV i, pj (t k ) = [x j (t k ) y j (t k ) z j (t k )] T represents the position quantity of the neighbor node UAV j, represents the desired position quantity of the follower UAV i and the neighbor node j. Then, the cooperative control objective of the follower UAV i is described as:

[0053] When t k → ∞, p l (t k ) - p i (t k ) + d li = 0, q l (t k ) - q i (t k ) = 0, p i (t k ) - p j (t k ) + d ij = 0, u i (t k ) - u d (t k ) = 0;

[0054] In solving the distributed model predictive control optimization problem based on the leader-follower UAV formation method, an asynchronous solution order is set between the leader and the followers, and a synchronous solution order is adopted between the followers. That is, the leader first solves the optimization problem at time t k , and obtains the control sequence u k (τ; t k+1 ..., t k+N-1 ) and the corresponding state sequence x i (τ; t k ) when τ = t i (τ; t k ), and uses them as the reference values of the leader's position and control input in solving the followers' optimization problem at time t k .

[0055] Furthermore, the independent random data loss includes:

[0056] When the UAV has a data loss problem in the case of discontinuous communication, it is assumed that the communication data loss does not occur at two consecutive moments, and the probability of data loss at different moments is the same and independent of each other. Define the random variable δ to represent whether data is lost. Then, there is:

[0057]

[0058] Meanwhile, δ satisfies:

[0059] δ(t k )δ(t k+1 ) = 0 (12)

[0060] where δ(t k ) and δ(t k+1 ) respectively represent the values of the random variable δ at times t k and t k+1 . It is assumed that the data loss probability is the same among all UAVs.

[0061] Furthermore, step 3 includes:

[0062] The design of the trajectory tracking controller for the leader UAV includes:

[0063] The discrete model of equation (3) for the leader UAV l is abbreviated in the following form:

[0064] x l (t k+1 ) = f(x l (t k ), u l (t k )) (13)

[0065] where x l (t k ) = [x l (t k ) y l (t k ) z l (t k ) θ l (t k ) ψ l (t k )] T represents the state quantity of the leader UAV l at time t k , and u l (t k ) = [v l (t k ) ω l (t k ) r l (t k )] T represents the control quantity of the leader UAV l at time t k . The state constraint of the system is x l (t k ) ∈ χ land the input constraint u l (t k ) ∈ Ω l , where χ l and Ω l are closed sets containing the control objectives;

[0066] According to the control objective, define the variable p ld (t k ) = p l (t k ) - p d (t k ) represents the position tracking error of the leader UAV l, q ld (t k ) = q l (t k ) - q d (t k ) represents the angular measurement tracking error of the leader UAV l, u ld (t k ) = u l (t k ) - u d (t k ) represents the control quantity tracking error of the leader UAV l;

[0067] Design the cost function for the trajectory tracking control optimization problem of the leader UAV l at time t k as follows:

[0068]

[0069] where τ = t k , t k+1 ..., t k+N , N represents the prediction horizon of the distributed model predictive control algorithm, and the terminal cost function is related to the stability of the terminal feedback controller, and the rationality of the design will be proved in the stability analysis section. is the state variable, Q l1 and Q l2 are the weight matrices corresponding to the state variables, R l is the weight matrix of the control quantity, and Q' and R' are the weight matrices in the terminal cost with respect to at time t k and t k+1 . p ld (τ; t k ), q ld (τ; t k ) and u ld (τ; t k ) represent the predicted values of the state and control quantities at time τ at time t k ;

[0070] The optimization problem solved by the leader UAV l at time t is as follows: k The optimization problem solved at time t is as follows:

[0071]

[0072] where, Ω l is the control constraint domain, χ l is the state constraint domain, μ ∈ [0, 1), C l (ε l (τ; t k )) is the compatibility constraint, which includes the collision constraint of the UAV formation flight and the limiting constraint related to the algorithm stability, p ld (t k+N ; t k ), q ld (t k+N ; t k ) and u ld (t k+N ; t k ) represent the predicted values of the state and control variables at time t k for the state and control variables at time t k+N , p ld (t k ), q ld (t k ) and u ld (t k ) represent the actual values of the state and control variables at time t k ; at time t k , the optimization problem in equation (15) is solved to obtain the optimal control sequence as:

[0073]

[0074] In the above formula represents the optimal control sequence obtained after solving the optimization problem;

[0075] Take the optimized value of the control sequence at time t k as the control input of the system, that is:

[0076]

[0077] At each time t k Solve the optimization problem, and apply the optimized solution obtained from equation (17) to the system. At the same time, the prediction time domain is rolled from [t k , t k+N to [t k+1 , t k+1+N , and repeat the solution of the optimization problem in equation (15) at time t k+1 to achieve rolling optimization.

[0078] Furthermore, the design of the follower UAV cooperative controller includes:

[0079] The discrete model of the follower UAV i in Equation (3) is abbreviated as the following form

[0080] x i (t k+1 ) = f(x i (t k ), u i (t k )) (18)

[0081] where x i (t k ) = [x i (t k ) y i (t k ) z i (t k ) θ i (t k ) ψ i (t k )] T represents the state quantity of the follower UAV i at time t k , and u i (t k ) = [v i (t k ) ω i (t k ) r i (t k )] T represents the control quantity of the follower UAV i at time t k . The state constraint of the system x i (t k ) ∈ χ i and the input constraint u i (t k ) ∈ Ω i , where χ i and Ω i are closed sets containing the control objectives; according to the control objectives, the variable p li (t k ) = p l (t k ) - p i (t k ) + d li represents the position cooperation error between the leader UAV l and the follower UAV i, and q li (t k ) = q l (t k ) - qi (t k ) represents the angular coordination error between the leading UAV l and the follower UAV i, p ij (t k ) = p i (t k ) - p j (t k ) + d ij represents the position coordination error between the follower UAV i and the neighbor node UAV j, u id (t k ) = u i (t k ) - u d (t k ) represents the control quantity tracking error of the follower UAV i;

[0082] Design the k cost function for the cooperative control optimization problem of the follower UAV i at time t as:

[0083]

[0084] where τ = t k , t k+1 ..., t k+N , N represents the prediction time domain of the distributed model predictive control algorithm, represents the cooperative cost function of the follower UAV i, and the terminal cost function is related to the stability of the terminal feedback controller, is the state quantity in equation (10), Q i1 and Q i2 are the weight matrices corresponding to the state quantities, R i is the weight matrix of the control quantity, Q′ i and R′ i are for the at time t k and t k+1 in the terminal cost, Q ij is the cooperative weight matrix, p li (τ; t k ), q li (τ; t k ) and u li (τ; t k ) represent the predicted values of the state and control quantities at time τ at time t k ;

[0085] The design of the terminal feedback controller includes:

[0086] After solving the optimization problems of equations (15) and (27), The control sequence; when τ ∈ [t k ,..., t k+N , Ensure that the cost function decreases gradually. When τ > t k+N , it is necessary to design a terminal feedback controller so that when τ approaches ∞, the control objective is achieved;

[0087] At time t k , after solving the optimization problem, the value of u(t k+N ; t k+N ) and x(t k ; t k+N+1 ) at time t k+N+1 ; t k ) is already near the reference point. At this time, the model is linearized to replace the original system model. The linearized error model is shown in Equation (10);

[0088] Design the terminal feedback controller as:

[0089]

[0090] where δ is the random variable introduced above, satisfying the mathematical expectation E(δ) = n, n represents the data loss rate, and υ is obtained by solving Equation (46). For detailed reference, see Equation (10); Substituting Equation (37) into Equation (10) and arranging, we get:

[0091]

[0092] Furthermore, the step 4 includes:

[0093] If there is no data loss caused by communication interruption between the leader UAV and the follower UAV at time t k , since the asynchronous communication method is adopted between the leader UAV and the follower UAV, the follower UAV can obtain the true predicted state sequence and the true predicted control sequence Then the predicted trajectory at time t k+1 is expressed as:

[0094]

[0095] represents the true predicted control sequence obtained by the follower UAV i by solving the optimization problem at time t k . Take the first element of the sequence as the optimal control input to act on the system; represents the assumed predicted control input of the follower UAV i, which is used to generate the assumed predicted trajectory Solve the optimization problem; where, includes and

[0096] If there is no data loss caused by communication interruption between the follower UAV and the neighbor node UAV at time t k , since they use synchronous communication, the follower UAV cannot obtain the true predicted state sequence of the neighbor node UAV and the true predicted control sequence and can only design a hypothetical predicted control sequence, expressed as:

[0097]

[0098] where κ j (x j (·)) represents the terminal feedback controller, and its specific form is given in equation (37); the above equation means that the hypothetical predicted state sequence of neighbor node UAV j received at time t k is composed of the control sequences at times t k-1 calculated by solving the optimization problem at time t k ,..., t k-2+N and the control quantity at time t k-1+N calculated by the terminal feedback controller;

[0099] By assuming the predicted control sequence in equation (21) and the system model in equation (18), a hypothetical predicted state sequence is generated, expressed as:

[0100]

[0101] where is calculated from κ j (x j (·)) in equation (21) and the system model in equation (18);

[0102] If there is data loss caused by communication interruption between the leader UAV and the follower UAV at time t k , the follower UAV can only obtain the hypothetical predicted control sequence of the leader UAV, expressed as:

[0103]

[0104] By assuming the predicted control sequence in equation (23) and the system model in equation (18), a hypothetical predicted state sequence is generated, expressed as:

[0105]

[0106] If at time t kWhen data loss occurs due to communication loss between the follower UAV and neighbor node UAVs, the follower UAV can only obtain the hypothetical prediction control sequence of the neighbor node UAVs. According to the descriptions of the loss of independent random data δ in equations (11) and (12), the hypothetical prediction control sequence is expressed as:

[0107]

[0108] By assuming the prediction control sequence in equation (25) and the system model in equation (18), a hypothetical prediction state sequence is generated, expressed as:

[0109]

[0110] Furthermore, step 5 includes:

[0111] The optimization problem solved by the follower UAV i at time t k is as follows:

[0112]

[0113] where Ω i is the control constraint domain, χ i is the state constraint domain, μ ∈ [0, 1) is a constant, R is the safety distance during UAV formation flight, C i (ε i (τ; t k )) is the compatibility constraint, including the collision constraint of UAV formation flight and the limiting constraints related to algorithm stability, p li (t k+N ; t k ), q li (t k+N ; t k ) and u li (t k+N ; t k ) represent the predicted values of the state and control quantity at time t k for the state and control quantity at time t k+N , p li (t k ), q li (t k ) and u li (t k ) represent the actual values of the state and control quantity at time t k , Δ ij is added as a soft constraint to perform constraint contraction on ;

[0114] At time t k , the optimized control sequence obtained by solving the optimization problem in equation (27) is:

[0115]

[0116] In the above formula represents the optimal control sequence obtained after solving the optimization problem of formula (27);

[0117] Take the optimized value of the control sequence at time t k as the control input of the system, that is:

[0118]

[0119] At each time t k Solve the optimization problem of formula (27), and apply the optimized solution obtained from formula (29) to the system. At the same time, the prediction horizon of the algorithm is rolled from [t k , t k+N to [t k+1 , t k+1+N , and repeat the solution of the optimization problem of formula (27) at time t k+1 to achieve rolling optimization.

[0120] Furthermore, when solving the optimization problem of formula (27), it must be satisfied that:

[0121] ||p i (τ; t k ) - p j (τ; t k )|| > 2R, τ ∈ [t k ,..., t k+N (30)

[0122] where R is the safety distance during the formation flight of the UAVs;

[0123] There is an error ε j (τ; t k ) between p and j (τ; t k ), which is defined as:

[0124]

[0125] In the above formula is the estimated value of the position of neighbor UAV j at time τ ∈ [t k ,..., t k+N , and p j (τ; t k ) is the actual value of the position of neighbor UAV j at time τ ∈ [t k ,..., t k+N ;

[0126] Due to the error ε j (τ; t k ), the constraint in equation (30) needs to be shrunk, expressed as:

[0127]

[0128] where R is the safety distance during the formation flight of the UAVs, Δ ij is expressed as:

[0129]

[0130] In solving the optimization problem, the error between the assumed predicted state sequence and the actual predicted state sequence is restricted, expressed as:

[0131] ||ε j (τ; t k )|| ≤ min{Δ i (t k ), η i (t k )} (34)

[0132] The above equation is the expression form of the compatibility constraint C i (ε i (τ; t k )) in the optimization problem of equation (27), where Δ i (t k ) and η i (t k ) are respectively expressed as:

[0133]

[0134] where α ∈ [0, 1) is a constant related to the convergence speed, T is the step size between time t k and t k+1 , N is the prediction horizon length,

[0135] Q ij is the cooperation weight matrix, λ max (Q ij ) is the maximum eigenvalue of the matrix Q ij , η i (t k ) is an index related to the algorithm stability.

[0136] Due to the above technical solution, the present application has the following advantages: realizing the constrained distributed cooperative control of multiple unmanned aerial vehicles (UAVs) in a leader-follower formation mode when data loss occurs in non-continuous communication scenarios, i.e., non-persistent communication scenarios. BRIEF DESCRIPTION OF THE DRAWINGS

[0137] To more clearly illustrate the technical solutions in the embodiments of the present application, the following briefly introduces the drawings required for the description of the embodiments. Obviously, the drawings in the following description are only some embodiments recorded in the embodiments of the present application, and those of ordinary skill in the art can obtain other drawings based on these drawings.

[0138] Figure 1 Schematic diagram of the path of a distributed model predictive control method for UAV formation in non-continuous communication scenarios according to an embodiment of the present application;

[0139] Figure 2 Schematic diagram of the inertial coordinate system and body coordinate system according to an embodiment of the present application;

[0140] Figure 3 Schematic diagram of yaw-pitch coordinate transformation according to an embodiment of the present application;

[0141] Figure 4 Trajectory tracking diagram according to an embodiment of the present application;

[0142] Figure 5 Error diagram in the x-axis direction according to an embodiment of the present application;

[0143] Figure 6 Error diagram in the y-axis direction according to an embodiment of the present application;

[0144] Figure 7 Error diagram in the z-axis direction according to an embodiment of the present application;

[0145] Figure 8 Diagram of the θ state variable according to an embodiment of the present application;

[0146] Figure 9 Diagram of the ψ state variable according to an embodiment of the present application;

[0147] Figure 10 Diagram of the v control variable according to an embodiment of the present application;

[0148] Figure 11 Diagram of the ω control variable according to an embodiment of the present application;

[0149] Figure 12 Diagram of the r control variable according to an embodiment of the present application;

[0150] Figure 13 Schematic diagram of the UAV communication connection mode according to an embodiment of the present application. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0151] The present application will be further described in conjunction with the accompanying drawings and embodiments. The described embodiments are only a part of the embodiments of the present application, rather than all of the embodiments. All other embodiments obtained by those of ordinary skill in the art shall fall within the scope of protection of the embodiments of the present application.

[0152] Currently, most multi-UAV formation cooperative control methods generally cannot handle the cooperative control problems in the case of nonlinear systems, communication data loss, and state constraints. The unique control mechanism of the model predictive control method can handle the constrained control problems of nonlinear multi-input multi-output coupled systems, and has the characteristics of fast convergence of control objectives, good optimization, and strong robustness. Refer to Figure 1 , so when considering the discontinuous communication between nodes during UAV formation flight, the present application adopts a distributed control structure and provides an embodiment of a distributed model predictive control method for UAV formation in the case of discontinuous communication, which includes:

[0153] Step 1: Establish a UAV kinematic model and a trajectory tracking error model;

[0154] Step 2: Simulate the independent random data loss situation of communication between UAVs;

[0155] Step 3: Based on the leader-follower formation method, establish the formation control objectives of the leader and the followers;

[0156] Step 4: Design a distributed model predictive controller for UAV formation in the case of discontinuous communication, including a trajectory tracking controller for the leader UAV, a cooperative controller for the follower UAVs, and a terminal feedback controller;

[0157] Step 5: Establish the predicted values of the reference trajectory of the leader UAV and the predicted values of the reference trajectories of the neighbor node UAVs by the cooperative controller of the follower UAVs within the prediction time domain;

[0158] Step 6: Solve the rolling optimization problem of each UAV at the current moment to obtain the optimal solution of the problem. Repeat the process of Steps 5-6 at the next moment until the set upper limit of the number of iterations.

[0159] (1) UAV kinematic model

[0160] The motion of a UAV belongs to a six-degree-of-freedom spatial motion, including forward and backward, left and right, up and down, yaw, pitch, and roll motions, which is a complex nonlinear coupling system. Considering the simplification of the large UAV model in engineering practice and ignoring the roll motion of the UAV, the kinematic model of UAV i in the formation system is established as follows:

[0161]

[0162] In the above formula, x i , y i and z i respectively represent the three-dimensional coordinate information of the unmanned aerial vehicle i in the inertial coordinate system, θ i and ψ i represent the pitch angle and the heading angle of the unmanned aerial vehicle i, v i represents the speed of the unmanned aerial vehicle i, ω i and r i represent the pitch angular velocity and the heading angular velocity of the unmanned aerial vehicle i. x i =[x i y i z i θ i ψ i T represents the state quantity of the unmanned aerial vehicle i, u i =[v i ω i r i T represents the control quantity of the unmanned aerial vehicle i. During the movement of the unmanned aerial vehicle i, the control quantity should be kept within the constraint range Ω i , that is, u i ∈Ω i , where Ω i is expressed as:

[0163]

[0164] In the above formula, v min and v max are respectively the minimum and maximum values of v i , ω min and ω max are respectively the minimum and maximum values of ω i , r min and r max are respectively the minimum and maximum values of r i .

[0165] Discretize the mathematical model in formula (1). Take the sampling time as T, and the following discrete kinematic model is obtained:

[0166]

[0167] All variables in the above formula are the discrete forms of the variables in formula (1).

[0168] (2) Unmanned Aerial Vehicle Trajectory Tracking Error Model

[0169] First, establish a reference trajectory model. Here, the reference trajectory of the unmanned aerial vehicle is determined by a virtual reference unmanned aerial vehicle with the same model, that is:

[0170] ​​

[0171] In the above formula, x d , y d and z d respectively represent the three-dimensional coordinate information of the reference trajectory in the inertial coordinate system, θ d and ψ d represent the pitch angle and the heading angle of the reference trajectory, v d represents the reference trajectory speed, ω d and r d represent the pitch angular velocity and the heading angular velocity of the reference trajectory, and their discretized forms are the same as those in formula (3).

[0172] Project the three-dimensional coordinate trajectory tracking error in the inertial coordinate system Oxyz onto the vehicle coordinate system O b x b y b z b As shown in Figure 2 and Figure 3 , it is mainly divided into two steps, including:

[0173] 1) Rotate ψ around the z-axis of the ground coordinate system to the coordinate system {G1}, and at this time, the x1-axis is located in the plane formed by the z-axis and the x b axis;

[0174] 2) Rotate θ around the y b axis of the coordinate system G1 to the coordinate system {B}, and at this time, it coincides with the vehicle coordinate system O b x b y b z b .

[0175] Without considering roll, the coordinate transformation matrix from the inertial coordinate system to the body coordinate system is:

[0176]

[0177] In the above formula, ψ is the heading angle, θ is the pitch angle, R z (ψ) and R y (θ) represent the coordinate transformation matrices for rotation around the z-axis and the y b axis respectively, and are respectively

[0178] In summary, the projection of the three-dimensional coordinate trajectory tracking error in the inertial coordinate system Oxyz onto the body coordinate system O b x b y b z b can be expressed as:

[0179]

[0180] Combined with the pitch angle θ and the yaw angle ψ trajectory tracking errors, the UAV trajectory tracking error can be described as:

[0181]

[0182] In the above formula, x ie , y ie and z ie are the three-dimensional coordinate errors between the UAV i and the reference trajectory in the body coordinate system, θ ie and ψ ie are the pitch angle and yaw angle errors between the UAV i and the reference trajectory.

[0183] Taking the derivative of equation (7) gives:

[0184]

[0185] The relevant parameters in the above formula refer to equations (1), (4) and (7).

[0186] Considering the subsequent design of the terminal feedback controller, the nonlinear trajectory tracking error model in equation (8) is linearized at the reference point to obtain:

[0187]

[0188] Define the control quantity as u ie =[u1 u2 u3] T , where u1 = (v i - v d )cos2θ d , and

[0189] Express the above formula in the general linear system form and discretize it to obtain:

[0190]

[0191] In the above formula T is the sampling time, I is the identity matrix with the same dimension as , The forms of A and B are as follows:

[0192]

[0193] (III) Problem description of UAV formation cooperative control under discontinuous communication

[0194] 1) Problem description of leader-follower formation cooperative control:

[0195] Select the leader-follower UAV formation method. The formation control problem can be divided into two parts:

[0196] a) Leader trajectory tracking problem:

[0197] Define p l (t k ) = [x l (t k ) y l (t k ) z l (t k )] T represent the position quantity of the leader UAV l, q l (t k ) = [θ l (t k ) ψ l (t k )] T represent the angular quantity of the UAV l, u l (t k ) = [v l (t k ) ω l (t k ) r l (t k )] T represent the control quantity of the UAV l, p d (t k ) = [x d (t k ) y d (t k ) z d (t k )] T represent the position quantity of the reference trajectory, q d (t k ) = [θ d (t k ) ψ d (t k )] T represent the angular quantity of the reference trajectory, u d (t k ) = [v d (t k ) ω d (t k ) r d (t k )] T represent the control quantity of the reference trajectory. Then, the control objective description of the leader UAV l trajectory tracking problem is:

[0198] When t k → ∞, p l (t k ) - p d (tk ) = 0, q l (t k ) - q d (t k ) = 0, u l (t k ) - u d (t k ) = 0。

[0199] b) Follower collaborative control problem:

[0200] Define the set of neighbor node drones j of follower drone i as Denote the desired position quantity of the leader drone and follower drone i, q i (t k ) = [θ i (t k ) ψ i (t k )] T Denote the angular quantity of follower drone i, p i (t k ) = [x i (t k ) y i (t k ) z i (t k )] T Denote the position quantity of follower drone i, p j (t k ) = [x j (t k ) y j (t k ) z j (t k )] T Denote the position quantity of neighbor node drone j, Denote the desired position quantity of follower drone i and neighbor node j, then the collaborative control objective of follower drone i is described as:

[0201] When t k → ∞, p l (t k ) - p i (t k ) + d li = 0, q l (t k ) - q i (t k ) = 0, p i (t k ) - p j (t k) + d ij = 0, u i (t k ) - u d (t k ) = 0。

[0202] In solving the distributed model predictive control optimization problem based on the leader-follower UAV formation method, an asynchronous solution order is set between the leader and the followers, and a synchronous solution order is adopted between the followers. That is, the leader first solves the optimization problem at time t k to obtain the control sequence u k , t k+1 ..., t k+N-1 and the corresponding state sequence x i (τ; t k ) when τ = t i (τ; t k ), and use them as the reference values of the leader's position and control input in solving the followers' optimization problem at time t k . The design of the above leader-follower cooperative control optimization problem, the predicted values of the control input and the predicted values of the neighbor node state variables are described in detail in the subsequent controller design.

[0203] 2) Description of the independent random data loss problem:

[0204] UAVs may have data loss problems in the case of discontinuous communication. It is assumed that data loss in communication does not occur at two consecutive times, and the probability of data loss at different times is the same and independent of each other. Define the random variable δ to represent whether data loss occurs, then there is:

[0205]

[0206] At the same time, δ satisfies:

[0207] δ(t k )δ(t k+1 ) = 0 (12)

[0208] where δ(t k ) and δ(t k+1 ) represent the values of the random variable δ at times t k and t k+1 respectively. It is assumed that the data loss probability is the same among all UAVs.

[0209] (IV) Design of the distributed model predictive controller for UAV formation in the case of discontinuous communication

[0210] In the case of discontinuous communication, the problem of cooperative control of UAVs first needs to solve the design problem of the cooperative controller for the leader UAV. Then, according to the solution of the optimization problem by the leader UAV and the predicted values of the states of the neighbor node UAVs, solve the design problem of the cooperative controller for the follower UAVs.

[0211] 1) Design of the trajectory tracking controller for the leader UAV

[0212] The discrete model of the leader UAV l in equation (3) is abbreviated as the following form:

[0213] x l (t k+1 ) = f(x l (t k ), u l (t k )) (13)

[0214] where x l (t k ) = [x l (t k ) y l (t k ) z l (t k ) θ l (t k ) ψ l (t k )] T represents the state quantity of the leader UAV l at time t k , u l (t k ) = [v l (t k ) ω l (t k ) r l (t k )] T represents the control quantity of the leader UAV l at time t k . The state constraint of the system x l (t k ) ∈ χ l and the input constraint u l (t k ) ∈ Ω l , where χ l and Ω l are closed sets containing the control objectives.

[0215] According to the control objectives described in the previous section, define the variable p ld (t k ) = p l (t k ) - p d (tk ) represents the position tracking error of the leader UAV l, q ld (t k ) = q l (t k ) - q d (t k ) represents the angular measurement tracking error of the leader UAV l, u ld (t k ) = u l (t k ) - u d (t k ) represents the control input tracking error of the leader UAV l.

[0216] According to the above definition, design the k cost function of the trajectory tracking control optimization problem of the leader UAV l at time t as:

[0217]

[0218] where τ = t k , t k+1 ..., t k+N , N represents the prediction horizon of the distributed model predictive control algorithm, and the terminal cost function is related to the stability of the terminal feedback controller, and the rationality of the design will be proved in the stability analysis section. is the state variable in equation (10), Q l1 and Q l2 are the weight matrices corresponding to the state variables, R l is the weight matrix of the control input, and Q' and R' are the weight matrices in the terminal cost with respect to at time t k and t k+1 . p ld (τ; t k ), q ld (τ; t k ) and u ld (τ; t k ) represent the predicted values of the state and control input at time τ at time t k .

[0219] According to the above analysis, the optimization problem solved by the leader UAV l at time t k is as follows:

[0220]

[0221] where Ω l is the control constraint domain, χ l is the state constraint domain, μ ∈ [0, 1), C l (εl (τ; t k )) is defined as the compatibility constraint in Equation (34) below, including the collision constraint for the formation flight of UAVs and the constraint related to the algorithm stability. p ld (t k+N ; t k ), q ld (t k+N ; t k ) and u ld (t k+N ; t k ) represent the predicted values of the state and control variables at time t k for the state and control variables at time t k+N . p ld (t k ), q ld (t k ) and u ld (t k ) represent the actual values of the state and control variables at time t k . Then, solving the optimization problem of Equation (15) at time t k can obtain the optimal control sequence as follows:

[0222]

[0223] In the above formula represents the optimal control sequence obtained after solving the optimization problem of Equation (15).

[0224] According to the basic principle of the model predictive control algorithm, take the optimal value of the control sequence at time t k as the control input of the system, that is:

[0225]

[0226] At each time t k solve the optimization problem of Equation (15), and apply the optimal solution obtained from Equation (17) to the system. At the same time, the prediction horizon of the algorithm rolls from [t k , t k+N to [t k+1 , t k+1+N , and repeat solving the optimization problem of Equation (15) at time t k+1 to achieve rolling optimization.

[0227] 2) Design of the cooperative controller for the follower UAVs

[0228] Abbreviate the discrete model of the follower UAV i in Equation (3) as the following form

[0229] x i (t k+1 ) = f(xi (t k ),u i (t k )) (18)

[0230] where x i (t k ) = [x i (t k ) y i (t k ) z i (t k ) θ i (t k ) ψ i (t k )] T represents the state variables of the follower UAV i at time t k , and u i (t k ) = [v i (t k ) ω i (t k ) r i (t k )] T represents the control variables of the follower UAV i at time t k . The state constraint of the system x i (t k ) ∈ χ i and the input constraint u i (t k ) ∈ Ω i , where χ i and Ω i are closed sets containing the control objectives. According to the control objectives described in the previous section, define the variable p li (t k ) = p l (t k ) - p i (t k ) + d li to represent the position coordination error between the leader UAV l and the follower UAV i, and q li (t k ) = q l (t k ) - q i (t k ) to represent the angle coordination error between the leader UAV l and the follower UAV i, and p ij (t k ) = p i (t k ) - p j (t k ) + dij represents the position coordination error between follower UAV i and neighbor node UAV j, u id (t k )=u i (t k )-u d (t k ) represents the tracking error of the follower UAV i control quantity.

[0231] According to the above definition, design t k The cost function of the coordinated control optimization problem of the follower UAV i is:

[0232]

[0233] Where τ = t k ,t k+1 ...,t k+N , N represents the prediction time domain of the distributed model predictive control algorithm, represents the collaborative cost function of follower UAV i, and the terminal cost function It is related to the stability of the terminal feedback controller. The rationality of the design will be proved in the stability analysis section later. is the state quantity in formula (10), Q i1 and Q i2 is the weight matrix of the corresponding state quantity, R i is the control weight matrix, Q′ i and R′ i For the terminal cost In t k and t k+1 The weight matrix at the moment, Q ij is the collaborative weight matrix, p li (τ; t k ),q li (τ; t k ) and u li (τ; t k ) indicates that at t k The predicted value of the state and control quantity at time τ.

[0234] In the above cooperative control problem, since at each sampling time t k , each follower drone solves the optimization problem synchronously, so for follower drone i, it is impossible to obtain the true predicted trajectory of neighbor node follower drone j. In addition, if the communication data between the leader drone and the follower drone, and between the follower drone and the neighbor node is lost, the true predicted trajectory cannot be obtained. In order to solve this problem, neighbor node drone j needs to send a hypothetical predicted trajectory to follower drone i before solving the optimization problem. The form of the hypothetical predicted trajectory is explained below:

[0235] Indicates that the follower drone i is at t k The real predictive control sequence obtained by solving the optimization problem at all times, takes the first element of the sequence as the optimal control input to act on the system.

[0236] represents the hypothetical predicted control input of follower UAV i, which is used to generate the hypothetical predicted trajectory Solve optimization problems.

[0237] in, include and

[0238] If in t k There is no data loss caused by communication interruption between the leader UAV and the follower UAV at all times. Since they use asynchronous communication, the follower UAV can obtain the real predicted state sequence. and the true predictive control sequence Then t k+1 The moment prediction trajectory is expressed as:

[0239]

[0240] If in t k There is no data loss caused by communication interruption between the follower drone and the neighbor node drone. Since they use synchronous communication, the follower drone cannot obtain the real predicted state sequence of the neighbor node drone. and the true predictive control sequence We can only design a hypothetical predictive control sequence, which can be expressed as:

[0241]

[0242] Among them, κ j (x j (·)) represents the terminal feedback controller, and its specific form is given in equation (37). The above equation shows that at t k The hypothetical predicted state sequence of the neighboring node UAV j received at time t is given by k-1 Solve the optimization problem at the moment t k ,...,t k-2+N The control sequence at time t and the t calculated by the terminal feedback controller k-1+N The control quantity composition at each moment.

[0243] By assuming the predictive control sequence of formula (21) and the unified model of formula (18), the hypothetical predicted state sequence is generated, which is expressed as:

[0244]

[0245] Among them, is calculated from κ in Equation (21) j (x j (·)) and the system model of Equation (18).

[0246] If data loss occurs due to communication interruption between the leader UAV and the follower UAV at time t k , the follower UAV can only obtain the assumed prediction control sequence of the leader UAV, which is expressed as:

[0247]

[0248] By assuming the control sequence of Equation (23) and the control system model of Equation (18), an assumed prediction state sequence is generated, which is expressed as:

[0249]

[0250] If data loss occurs due to communication loss between the follower UAV and the neighbor node UAV at time t k , the follower UAV can only obtain the assumed prediction control sequence of the neighbor node UAV. According to the descriptions of the independent random data δ loss in Equations (11) and (12), the assumed prediction control sequence can be expressed as:

[0251]

[0252] By assuming the prediction control sequence of Equation (25) and the system model of Equation (18), an assumed prediction state sequence is generated, which is expressed as:

[0253]

[0254] According to the above analysis, the optimization problem solved by the follower UAV i at time t k is as follows:

[0255]

[0256] Among them, Ω i is the control constraint domain, χ i is the state constraint domain, μ ∈ [0, 1) is a constant, R is the safety distance during UAV formation flight, C i (ε i (τ; t k )) is the compatibility constraint, including the collision constraint of UAV formation flight and the defined constraints related to algorithm stability, which is defined in Equation (34) later, p li (t k+N ; tk ),q li (t k+N ;t k ) and u li (t k+N ;t k ) indicates that at t k The state and control quantity at time t k+N The predicted value at time, p li (t k ),q li (t k ) and u li (t k ) indicates that at t k The actual value of the state and control quantity at the moment, Δ ij Defined in formula (33), Here, considering that there is a difference between the assumed predicted state sequence and the actual predicted state sequence, the increase of Δ ij As a soft constraint, Perform constrained shrinkage.

[0257] Then in t k Solving the optimization problem of (27) at all times can obtain the optimal control sequence:

[0258]

[0259] In the above formula represents the optimal control sequence obtained after solving the optimization problem of formula (27)

[0260] According to the basic principle of model predictive control algorithm, the control sequence t k The optimal value at the moment is used as the control input of the system, which is:

[0261]

[0262] At every moment t k Solve the optimization problem of (27), and apply the optimization solution obtained from (29) to the system. At the same time, the prediction time domain of the algorithm is [t k ,t k+N ]Scroll to [t k+1 ,t k+1+N ], and at t k+1 Repeatedly solve the optimization problem (27) to achieve rolling optimization.

[0263] (V) Description of UAV formation flight constraints

[0264] When UAVs are flying in formation, collisions between them must be avoided. So when solving the optimization problem of equation (27), the following must be satisfied:

[0265] ||p i (τ; t k ) - p j (τ; t k )|| > 2R, τ ∈ [t k ,..., t k+N (30)

[0266] Among them, R is the safety distance during the formation flight of the UAVs.

[0267] Since the optimization algorithm adopts a synchronous solution mode and there is data loss caused by discontinuous communication, in the above design, an assumed predicted state sequence is used to replace the actual predicted state sequence. Then p j (τ; t k ) and There is an error ε j (τ; t k ), which is defined as:

[0268]

[0269] In the above formula Is the estimated value of the position quantity of neighbor node UAV j at the moment τ ∈ [t k ,..., t k+N , p j (τ; t k ) is the actual value of the position quantity of neighbor node UAV j at the moment τ ∈ [t k ,..., t k+N .

[0270] Due to the existence of the error ε j (τ; t k ), it is necessary to contract the constraint of formula (30), which is expressed as:

[0271]

[0272] Among them, R is the safety distance during the formation flight of the UAVs. Δ ij Is expressed as:

[0273]

[0274] Considering the formation constraint requirements and algorithm stability, the error between the assumed predicted state sequence and the actual predicted state sequence is restricted in the solution of the optimization problem, which is expressed as:

[0275] ||ε j (τ; t k )|| ≤ min{Δ i (t k ), η i (tk )} (34)

[0276] The above formula is the compatibility constraint C in the optimization problem of formula (27): i (ε i (τ; t k )) in the form of i (t k ) and η i (t k ) are respectively expressed as:

[0277]

[0278]

[0279] Among them, α∈[0,1) is a constant related to the convergence rate, T is t k to k+1 The step length between moments, N is the length of the prediction time domain, Q ij is the collaborative weight matrix, λ max (Q ij ) is the matrix Q ij The maximum eigenvalue of i (t k ) is an indicator related to the stability of the algorithm.

[0280] (VI) Terminal feedback controller design

[0281] After solving the optimization problems of equations (15) and (27), we can obtain The control sequence of k ,...,t k+N ]hour It can be ensured that the cost function is gradually decreasing, and when τ>t k+N It is necessary to design a terminal feedback controller so that when τ approaches ∞, the control objective can be achieved.

[0282] In t k After solving the optimization problem at all times, we get t k+N u(t k+N ;t k ) and t k+N+1 Time x(t k+N+1 ;t k ) value is already near the reference point. At this time, the model can be linearized to replace the original system model. The error model after linearization is shown in formula (10).

[0283] Design the terminal feedback controller as:

[0284]

[0285] Among them, δ is the random variable introduced above, satisfying the mathematical expectation E(δ)=n, where n represents the data loss rate, and υ is obtained by solving equation (46). Refer to equation (10) in detail. Substituting equation (37) into equation (10) and after arrangement, we get:

[0286]

[0287] Next, analyze the stability of the terminal feedback controller of the designed equation (37). Select the Lyapunov function as:

[0288]

[0289] Among them, P and Q are real symmetric positive definite weight matrices. Then we have ΔV(t k ) as:

[0290]

[0291] Substituting the error model of equation (10) into the above formula, we can get:

[0292]

[0293] Among them, In the above formula T is the sampling time, I is the identity matrix with the same dimension as , The forms of A and B are as follows:

[0294]

[0295] Take the mathematical expectation of ΔV(t k ) in equation (41), then:

[0296]

[0297] Among them,

[0298] To make equation (10) have t k →∞ under the control of the terminal feedback controller of equation (37), it is necessary to make E(ΔV(t k ))<0, then we have:

[0299]

[0300] That is:

[0301]

[0302] Converted to:

[0303]

[0304] The above formula can be transformed according to the Schur complement theorem into:

[0305]

[0306] where P and Q are real symmetric positive definite weight matrices, and E(δ) = n. The terminal feedback control gain υ can be obtained by solving the above linear matrix inequality.

[0307] (VII) Feasibility and Stability Analysis of Distributed Model Predictive Controller

[0308] 1) Feasibility Analysis

[0309] This part analyzes the feasibility of the distributed model predictive control algorithm for the UAV formation under discontinuous communication. Distributed model predictive control is a rolling horizon optimization method. Therefore, at each update moment, the algorithm must have a feasible solution that satisfies the constraint conditions of the optimization problem for the optimization problem to continue.

[0310] Here, it is assumed that at time t k , there exists an optimal solution for UAV i. Then the optimal solution must satisfy the constraint conditions of the optimization problems given by equations (15) and (27). Then, at the next sampling time t k+1 , the control sequence can be constructed as follows:

[0311]

[0312] where κ j (x j (·)) represents the terminal feedback controller, and the specific form is given in equation (37).

[0313] According to equation (47), when τ ∈ [t k+1 ,..., t k+N , Then the corresponding state trajectory is Because and are the optimal solutions obtained at time t k , obviously they satisfy the constraint conditions in the optimization problem. When τ = t k+N+1 , the system control law is At this time, the control quantity and the state quantity are already in a very small neighborhood near the reference point, and obviously they also satisfy the constraint conditions in the optimization problem. Then, through mathematical induction, it can be known that at each update moment t k , a feasible solution can be found for the optimization problem to satisfy the constraint conditions. Therefore, the optimization problem is feasible. Q.E.D.

[0314] 2) Stability analysis

[0315] For the UAV formation system, select the Lyapunov function at time t as follows: k

[0316]

[0317] where N a represents the set of follower UAVs i. In the above Lyapunov function, J l (·) represents the cost function of the leader UAV l, and J i (·) represents the cost function of the follower UAV i, and their specific forms are given by equations (14) and (19). According to the optimality of the cost function, then we have:

[0318]

[0319] For the terminal feedback controller designed according to equation (37), we know that:

[0320]

[0321] where

[0322]

[0323]

[0324]

[0325]

[0326] τ = t k , t k+1 ..., t k+N , N represents the prediction horizon of the distributed model predictive control algorithm, Q l1 and Q l2 are the weight matrices of the corresponding state variables of the leader UAV l, R l is the weight matrix of the control variable of the leader UAV l, p ld (t k+N+1 ; t k ) and q ld (t k+N+1 ; t k ) represent the predicted values of the state variables of the leader UAV l at time t k for the state variables at time t k+N+1 , and u ld (t k+N ; t k ​) represents the predicted value of the control quantity of the leading UAV l at time t k at time t k+N , p ld (t k ; t k ) and q ld (t k ; t k ) represent the actual value of the leading UAV l at time t k , u ld (t k ; t k ) represent the actual value of the leading UAV l at time t k , Q i1 and Q i2 are the weight matrices of the corresponding state quantities of the following UAV i, R i is the weight matrix of the control quantity of the following UAV i, Q ij is the collaborative weight matrix, p li (t k+N+1 ; t k ) and q li (t k+N+1 ; t k ) represent the predicted value of the state quantity of the following UAV i at time t k at time t k+N+1 , u li (t k+N ; t k ) represent the predicted value of the control quantity of the following UAV i at time t k at time t k+N , p li (t k ; t k ) and q li (t k ; t k ) represent the actual value of the following UAV i at time t k , u li (t k ; t k ) represent the actual value of the following UAV i at time t k .

[0327] Define variables

[0328]

[0329] Then we have:

[0330]

[0331] Arrange the above formula:

[0332]

[0333] Define Since x T Ax ≤ λ max (A)||x|| 2 , Equation (54) is rearranged as:

[0334]

[0335] According to the compatibility constraint C i (ε i (τ; t k )) where η i (t k ) is defined by Equation (36), combined with Equations (31) and (34), the above equation is transformed into:

[0336]

[0337] Also according to Equation (34), the above equation is transformed into:

[0338]

[0339] Returning to Equation (49), according to the constraint conditions of Equations (15) and (27), we have:

[0340]

[0341] After rearranging Equation (49), we get:

[0342]

[0343] According to the previous description, μ ∈ [0, 1), α ∈ [0, 1). Then when α + μ - 1 < 0,, V(t k+1 ) - V(t k ) ≤ 0, the Lyapunov function of the system decreases, so the system is closed-loop stable. Q.E.D.

[0344] To better explain the technical solution of the present invention, the present invention also provides a more specific implementation example:

[0345] The following is a detailed description in combination with a specific simulation example. Taking three unmanned aerial vehicles as an example, where unmanned aerial vehicle 1 is the leading unmanned aerial vehicle, and unmanned aerial vehicles 2 and 3 are follower unmanned aerial vehicles. The communication connection method between the unmanned aerial vehicles is as Figure 13 shown:

[0346] S1. Establish the kinematic model of the unmanned aerial vehicle and the trajectory tracking error model:

[0347] 1) Considering the simplification of large - scale UAV models in engineering practice and ignoring the roll motion of the UAV, the kinematic model of UAV i in the formation system is established as follows:

[0348]

[0349] where x i , y i and z i respectively represent the three - dimensional coordinate information of UAV i in the inertial coordinate system, θ i and ψ i represent the pitch angle and heading angle of UAV i, v i represents the speed of UAV i, ω i and r i represent the pitch angular velocity and heading angular velocity of UAV i. x i = [x i y i z i θ i ψ i T represents the state quantity of UAV i, and u i = [v i ω i r i T represents the control quantity of UAV i. The constraint range of UAV i during motion is:

[0350] - 5 ≤ v ≤ 5

[0351] - π / 4 ≤ ω ≤ π / 4

[0352] - π / 4 ≤ r ≤ π / 4

[0353] 2) UAV trajectory tracking error model:

[0354]

[0355] where the control quantity is u ie = [u1 u2 u3] T , u1 = (v i - v d )cos2θ d , and The reference trajectory is determined by a virtual reference UAV with the same parameters, that is:

[0356]

[0357] S2. Simulate the independent random data loss situation in communication between UAVs:

[0358] Define the δ random variable to represent whether data loss occurs, that is:

[0359]

[0360] Assume E(δ)=0.3, establish a sequence of δ random variables to simulate the communication state between UAVs.

[0361] S3. Based on the leader-follower formation method, establish the formation control objectives of the leader and followers:

[0362] 1) The trajectory tracking control objective of the leader:

[0363] When t k →∞, p l (t k )-p d (t k )=0, q l (t k )-q d (t k )=0, u l (t k )-u d (t k )=0.

[0364] Among them, p l (t k )=[x l (t k ) y l (t k ) z l (t k )] T represents the position quantity of the leader UAV l, q l (t k )=[θ l (t k ) ψ l (t k )] T represents the angular quantity of the UAV l, u l (t k )=[v l (t k ) ω l (t k ) r l (t k )] T represents the control quantity of the UAV l, p d (t k )=[x d (t k ) y d (t k ) zd (t k )] T represents the position quantity of the reference trajectory, q d (t k )=[θ d (t k ) ψ d (t k )] T represents the angular quantity of the reference trajectory, u d (t k )=[v d (t k ) ω d (t k ) r d (t k )] T represents the control quantity of the reference trajectory.

[0365] 2) Follower cooperative control objective:

[0366] When t k →∞, p l (t k ) - p i (t k ) + d li =0, q l (t k ) - q i (t k )=0, p i (t k ) - p j (t k ) + d ij =0, u i (t k ) - u d (t k )=0.

[0367] Among them, the set of neighbor node drones j of the follower drone i is represents the expected position quantity of the leader drone and the follower drone i, q i (t k )=[θ i (t k ) ψ i (t k )] T represents the angular quantity of the follower drone i, p i (t k )=[x i (t k ) y i (t k ) z i(t k )] T represents the position quantity of the follower UAV i, p j (t k ) = [x j (t k ) y j (t k ) z j (t k )] T represents the position quantity of the neighbor node UAV j, represents the desired position quantity between the follower UAV i and the neighbor node j, where are respectively set to -1 and 1, are respectively set to -2 and 2.

[0368] S4. Design of the distributed model predictive controller for UAV formation in the case of discontinuous communication:

[0369] 1) The trajectory tracking controller of the leader UAV is designed as follows:

[0370]

[0371] Among them, p l (t k ) = [x l (t k ) y l (t k ) z l (t k )] T is the position quantity of the leader UAV l, q l (t k ) = [θ l (t k ) ψ l (t k )] T is the angular quantity of the leader UAV l, p l (t k ) = [x l (t k ) y l (t k ) z l (t k )] T , p ld (t k ) = p l (t k ) - p d (t k ) represents the position tracking error of the leader UAV l, q ld (t k ) = ql (t k )-q d (t k ) represents the tracking error of the angle measurement of the leader UAV l, u ld (t k ) = u l (t k )-u d (t k ) represents the tracking error of the control quantity of the leader UAV l, τ = t k , t k+1 ..., t k+N , N represents the prediction horizon of 10, and the terminal cost function

[0372] where the initial value p l (t0) = [3 1 0.5] T , q l (t0) = [π*2 / 9 π / 8] T , u l (t0) = [0 0 0] T . Q l1 = diag([5,7,5]), Q l2 = diag([1,1]), R l = diag([0.5,0.5,0.5]), Q' = R' = diag([10,10,10,1,1]).

[0373] 2) The collaborative controller design of the follower UAV is as follows:

[0374]

[0375] where τ = t k , t k+1 ..., t k+N , N represents the prediction horizon of 10, represents the collaborative cost function of the follower UAV i, and the terminal cost function

[0376] where the initial value is Q i1 = diag([5,7,5]), Q l2 = diag([0.5,0.5]), R l = diag([0.01,0.5,0.5]), Q' = R' = diag([8,8,10,0.5,0.5]), Q ij = diag([0.5,0.5,0.5,0,0]).

[0377] 3) The terminal feedback controller is designed as follows:

[0378]

[0379] Among them, the terminal feedback control gain υ is obtained by solving the following matrix inequality.

[0380]

[0381] Among them, Q = diag([3, 3, 3]), P = diag([1.25, 1.25, 1.25]).

[0382] S5: Establish the predicted values of the reference trajectory of the leader UAV and the predicted values of the reference trajectory of the neighbor node UAV by the follower UAV collaborative controller within the prediction horizon N:

[0383] 1) Predicted value of the reference trajectory of the leader UAV:

[0384] If there is no data loss caused by communication interruption between the leader UAV and the follower UAV at time t k , since they use an asynchronous communication method, the follower UAV can obtain the true predicted state sequence and the true predicted control sequence Then the predicted trajectory at time t k+1 is expressed as:

[0385]

[0386] If there is no data loss caused by communication interruption between the follower UAV and the neighbor node UAV at time t k , since they use a synchronous communication method, the follower UAV cannot obtain the true predicted state sequence of the neighbor node UAV and the true predicted control sequence and can only design an assumed predicted control sequence. Then the predicted trajectory at time t k+1 is expressed as:

[0387]

[0388] Among them, κ j (x j (·)) represents the terminal feedback controller.

[0389] Through the above predicted control sequence and system model, an assumed predicted state sequence is generated, expressed as:

[0390]

[0391] 2) Predicted value of the reference trajectory of the neighbor node UAV:

[0392] If there is data loss caused by communication interruption between the leader UAV and the follower UAV at time t k , the follower UAV can only obtain the assumed prediction control sequence of the leader UAV, and the predicted trajectory at time t k+1 is expressed as:

[0393]

[0394] Through the above prediction control sequence and system model, an assumed prediction state sequence is generated, expressed as:

[0395]

[0396] If there is data loss caused by communication interruption between the follower UAV and the neighbor node UAV at time t k , the follower UAV can only obtain the assumed prediction control sequence of the neighbor node UAV, and the assumed prediction control sequence at time t k+1 can be expressed as:

[0397]

[0398] Through the above prediction control sequence and system model, an assumed prediction state sequence is generated, expressed as:

[0399]

[0400] S6: Solve the receding horizon optimization problem (15) of the leader UAV at time t k and the optimization problem (27) solved by the follower UAV i at time t k to obtain the optimal solution of the problem. Repeat the process of steps S5 - S6 at the next moment until the set upper limit of the number of iterations.

[0401] Figures 4 - 12 is the algorithm simulation effect diagram. Among them, Figure 5 the ordinate represents the error in the Ox - axis direction in the inertial coordinate system Oxyz, with the unit of meter (m), and the abscissa represents time, with the unit of second (s); Figure 6 the ordinate represents the error in the Oy - axis direction in the inertial coordinate system Oxyz, with the unit of meter (m), and the abscissa represents time, with the unit of second (s); Figure 7 the ordinate represents the error in the Oz - axis direction in the inertial coordinate system Oxyz, with the unit of meter (m), and the abscissa represents time, with the unit of second (s); Figure 8 the ordinate represents the pitch angle, with the unit of radian (rad), and the abscissa represents time, with the unit of second (s); Figure 9 the ordinate represents the heading angle, with the unit of radian (rad), and the abscissa represents time, with the unit of second (s); Figure 10The vertical axis represents speed, with the unit of meters per second (m / s), and the horizontal axis represents time, with the unit of seconds (s). Figure 11 The vertical axis represents pitch angular velocity, with the unit of radians per second (rad / s), and the horizontal axis represents time, with the unit of seconds (s). Figure 12 The vertical axis represents yaw angular velocity, with the unit of radians per second (rad / s), and the horizontal axis represents time, with the unit of seconds (s). The simulation results show that in the case of discontinuous communication, the distributed model predictive control UAV formation cooperative control algorithm has a good tracking effect on virtual UAVs with the same kinematic model. The three axes of the coordinate system converge, with a fast convergence speed, high accuracy, small overshoot, good stability, the control quantity meets the given constraint conditions, and the algorithm has a good control effect.

[0402] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present application and are not intended to limit them. Although the present application has been described in detail with reference to the above embodiments, those of ordinary skill in the art should understand that: the specific implementation manners of the present application can still be modified or equivalently replaced, and any modification or equivalent replacement that does not depart from the spirit and scope of the present application shall be covered by the protection scope of the claims of the present application.

Claims

1. A distributed model predictive control method for UAV formation in the case of discontinuous communication, characterized in that Including: Step 1: Establish the kinematic model and trajectory tracking error model of the UAV. Step 2: Simulate the independent random data loss situation in the communication between UAVs, and based on the leader-follower formation method, establish the formation control objectives of the leader and followers. Step 3: Design the distributed model predictive controller for the UAV formation in the case of discontinuous communication, including the trajectory tracking controller of the leader UAV, the cooperative controller of the follower UAVs, and the terminal feedback controller. Step 4: Establish the predicted values of the reference trajectory of the leader UAV and the predicted values of the reference trajectories of the neighboring node UAVs by the cooperative controller of the follower UAVs within the prediction time domain. Step 5: Solve the rolling optimization problem of each UAV at the current moment to obtain the optimal solution.

2. The method according to claim 1, wherein In the above-mentioned Step 1, the establishment process of the UAV kinematic model is as follows: Ignoring the roll motion of the UAV, the kinematic model of the UAV i in the formation system is established as follows: In the above formula, x i , y i and z i respectively represent the three-dimensional coordinate information of the UAV i in the inertial coordinate system, θ i and ψ i represent the pitch angle and heading angle of the UAV i, v i represents the speed of the UAV i, ω i and r i represent the pitch angular velocity and heading angular velocity of the UAV i; x i =[x i y i z i θ i ψ i T represents the state quantity of the UAV i, u i =[v i ω i r i T represents the control quantity of the UAV i;​​ During the movement of the unmanned aerial vehicle i, the control quantity should be kept within the constraint range Ω due to hardware factors i That is, u i ∈Ω i , where Ω i is expressed as: v min ≤v i ≤v max ω min ≤ ω i ≤ ω max r min ≤r i ≤r max (2) In the above formula, v min and v max are the minimum and maximum values of v i respectively, ω min and ω max are the minimum and maximum values of ω i respectively, r min and r max are the minimum and maximum values of r i respectively; Discretize the mathematical model in Equation (1), taking the sampling time as T, to obtain the following discretized kinematic model: All variables in the above equation are the discrete forms of the variables in Equation (1).

3. The method according to claim 2, wherein In the above-mentioned Step 1, the construction process of the trajectory tracking error model is as follows: First, establish the reference trajectory model. The reference trajectory of the UAV is determined by a virtual reference UAV with the same model, that is: In the above formula, x d , y d and z d respectively represent the three-dimensional coordinate information of the reference trajectory in the inertial coordinate system, θ d and ψ d represent the pitch angle and heading angle of the reference trajectory, v d represents the speed of the reference trajectory, ω d and r d represent the pitch angular velocity and heading angular velocity of the reference trajectory, and their discretized forms are the same as those in formula (3); Project the three-dimensional coordinate trajectory tracking error in the inertial coordinate system Oxyz onto the vehicle coordinate system O b x b y b z b There are mainly two steps, including: 1) Rotate by ψ about the z-axis of the ground coordinate system to the coordinate system {G1}, where the x1-axis lies in the plane formed by the z-axis and the x b -axis; 2) Rotate by θ about the y-axis of the coordinate system G1 to the coordinate system {B}, at this time it coincides with the vehicle coordinate system O b x b x b y b z b ; Coordinate transformation matrix from the inertial coordinate system to the body coordinate system without considering roll is as follows: In the above formula, ψ is the heading angle, θ is the pitch angle, and R z (ψ) and R y (θ) represent the coordinate transformation matrices for rotation about the z-axis and the y b axis, respectively, and are The projection of the three-dimensional coordinate tracking error in the inertial coordinate system Oxyz onto the body coordinate system O b x b y b z b is expressed as: Combining the pitch angle θ and the heading angle ψ trajectory tracking errors, the UAV trajectory tracking error is described as: In the above formula, x ie , y ie and z ie are the three-dimensional coordinate errors between the UAV i and the reference trajectory in the body coordinate system, and θ ie and ψ ie are the pitch angle and heading angle errors between the UAV i and the reference trajectory; After differentiating Equation (7), we get: Linearize the non-linear trajectory tracking error model in Equation (8) at the reference point to obtain: Define the control quantity as u ie = [u1 u2 u3] T , where u1 = (v i - v d ) cos2θ d , and Express Equation (8) in the general linear system form and discretize it to obtain: In the above formula where T is the sampling time, I is the identity matrix with the same dimension as and A and B are in the following forms:

4. The method according to claim 1, wherein In the above-mentioned Step 2, the leader-follower formation cooperative control objective and independent random data loss; The leader-follower formation cooperative control objective includes the leader trajectory tracking control objective and the follower cooperative control objective; a) Leader trajectory tracking control objective: Define p l (t k ) = [x l (t k ) y l (t k ) z l (t k )] T represents the position quantity of the leader UAV l, q l (t k ) = [θ l (t k ) ψ l (t k )] T represents the angular quantity of the UAV l, u l (t k ) = [v l (t k ) ω l (t k ) r l (t k )] T represents the control quantity of the UAV l, p d (t k ) = [x d (t k ) y d (t k ) z d (t k )] T represents the position quantity of the reference trajectory, q d (t k ) = [θ d (t k ) ψ d (t k )] T represents the angular quantity of the reference trajectory, u d (t k ) = [v d (t k ) ω d (t k ) r d (t k )] T represents the control quantity of the reference trajectory. Then, the control objective description of the trajectory tracking problem of the leader UAV l is: When t k → ∞, p l (t k ) - p d (t k ) = 0, q l (t k ) - q d (t k ) = 0, u l (t k ) - u d (t k ) = 0; b) Follower cooperative control objective: Define the set of neighbor node drones j of the follower drone i as Denote the desired position quantities of the leader drone and the follower drone i as q i (t k ) = [θ i (t k ) ψ i (t k )] T Denote the angular quantity of the follower drone i as p i (t k ) = [x i (t k ) y i (t k ) z i (t k )] T Denote the position quantity of the follower drone i as p j (t k ) = [x j (t k ) y j (t k ) z j (t k )] T Denote the position quantity of the neighbor node drone j Denote the desired position quantity between the follower drone i and the neighbor node j. Then, the cooperative control objective of the follower drone i is described as: When t k → ∞, p l (t k ) - p i (t k ) + d li = 0, q l (t k ) - q i (t k ) = 0, p i (t k ) - p j (t k ) + d ij = 0, u i (t k ) - u d (t k ) = 0; In solving the distributed model predictive control optimization problem based on the leader-follower UAV formation method, an asynchronous solution order is set between the leader and the followers, and a synchronous solution order is adopted between the followers, that is, the leader first solves the optimization problem at time t k to obtain the control sequence u k , t k+1 ..., t k+N-1 and the corresponding state sequence x i (τ; t k ) when τ = t i (τ; t k ), and use them as the reference values of the leader's position and control input in solving the followers' optimization problem at time t k .

5. The method according to claim 4, characterized in that The above-mentioned independent random data loss includes: When the UAV has a data loss problem in the case of discontinuous communication, it is assumed that the communication data loss situation does not occur at two consecutive moments, and the probability of data loss at different moments is the same and independent of each other. Define the random variable δ to represent whether data is lost, then we have: At the same time, δ satisfies: δ(t k )δ(t k+1 ) = 0 (12) where, δ(t k ) and δ(t k+1 ) respectively represent the values of the random variable δ at times t k and t k+1 , assuming that the data loss probabilities are the same among all UAVs.

6. The method according to claim 5, characterized in that The above-mentioned Step 3 includes: The design of the trajectory tracking controller of the leader UAV includes: Abbreviate the discrete model in Equation (3) for the leader UAV l as the following form: x l (t k+1 ) = f(x l (t k ), u l (t k )) (13) where x l (t k ) = [x l (t k ) y l (t k ) z l (t k ) θ l (t k ) ψ l (t k )] T represents the state variables of the leader UAV l at time t k , u l (t k ) = [v l (t k ) ω l (t k ) r l (t k )] T represents the control variables of the leader UAV l at time t k . The state constraints of the system x l (t k ) ∈ χ l and the input constraints u l (t k ) ∈ Ω l , where χ l and Ω l are closed sets containing the control objectives; Define the variable p according to the control objective ld (t k ) = p l (t k ) - p d (t k ) represents the position tracking error of the leader UAV l, q ld (t k ) = q l (t k ) - q d (t k ) represents the angle measurement tracking error of the leader UAV l, u ld (t k ) = u l (t k ) - u d (t k ) represents the control quantity tracking error of the leader UAV l; Design t k The cost function for the optimization problem of the trajectory tracking control of the Moment Navigator UAV l is as follows: where, τ = t k , t k+1 ..., t k+N , N represents the prediction horizon of the distributed model predictive control algorithm, and the terminal cost function is related to the stability of the terminal feedback controller, and the rationality of the design will be proved in the stability analysis section is the state variable, Q l1 and Q l2 are the weight matrices corresponding to the state variables, R l is the weight matrix of the control variable, and Q' and R' are the weight matrices in the terminal cost with respect to at t k and t k+1 moments, and p ld (τ; t k )、q ld (τ; t k ) and u ld (τ; t k ) represent the predicted values of the state and control variables at τ moment at t k moment; Pilot drone l at t k The optimization problem solved at the moment is as follows: where, Ω l is the control constraint domain, χ l is the state constraint domain, μ ∈ [0, 1), C l (ε l (τ; t k )) is the compatibility constraint, which includes the collision constraint of the UAV formation flight and the defining constraints related to the algorithm stability. p ld (t k+N ; t k ), q ld (t k+N ; t k ) and u ld (t k+N ; t k ) represent the predicted values of the state and control variables at time t k for time t k+N . p ld (t k ), q ld (t k ) and u ld (t k ) represent the actual values of the state and control variables at time t k . At time t k , solve the optimization problem in Equation (15) to obtain the optimal control sequence as follows: In the above formula represents the optimal control sequence obtained after solving the optimization problem; Take the optimized value of the control sequence t k at the moment as the control input of the system, that is: At each moment t k Solve the optimization problem, and apply the optimal solution obtained from equation (17) to the system. Meanwhile, the prediction time domain rolls from [t k , t k+N to [t k+1 , t k+1+N , and repeat to solve the optimization problem of equation (15) at moment t k+1 to achieve rolling optimization.

7. The method according to claim 6, wherein The design of the cooperative controller of the follower UAVs includes: Abbreviate the discrete model in Equation (3) for the follower UAV i as the following form x i (t k+1 ) = f(x i (t k ), u i (t k )) (18) where x i (t k ) = [x i (t k ) y i (t k ) z i (t k ) θ i (t k ) ψ i (t k )] T represents the state variables of the follower UAV i at time t k , u i (t k ) = [v i (t k ) ω i (t k ) r i (t k )] T represents the control variables of the follower UAV i at time t k . The state constraint of the system x i (t k ) ∈ χ i and the input constraint u i (t k ) ∈ Ω i , where χ i and Ω i are closed sets containing the control objectives. According to the control objectives, the variable p li (t k ) = p l (t k ) - p i (t k ) + d li represents the position coordination error between the leader UAV l and the follower UAV i, q li (t k ) = q l (t k ) - q i (t k ) represents the angle coordination error between the leader UAV l and the follower UAV i, p ij (t k ) = p i (t k ) - p j (t k ) + d ij represents the position coordination error between the follower UAV i and the neighbor UAV j, u id (t k ) = u i (t k ) - u d (t k ) represents the tracking error of the control quantity of the follower UAV i; Design t k The cost function for the cooperative control optimization problem of the follower UAV i at all times is as follows: where τ = t k , t k+1 ..., t k+N , N represents the prediction horizon of the distributed model predictive control algorithm represents the collaborative cost function of the follower UAV i, and the terminal cost function is related to the stability of the terminal feedback controller is the state variable in equation (10), Q i1 and Q i2 are the weight matrices corresponding to the state variables, R i is the weight matrix of the control variable, Q′ i and R′ i are the weight matrices in the terminal cost with respect to at time t k and t k+1 , Q ij is the collaborative weight matrix, p li (τ; t k ), q li (τ; t k ) and u li (τ; t k ) represent the predicted values of the state and control variables at time τ at time t k ; The design of the terminal feedback controller includes: After solving the optimization problems in equations (15) and (27), we obtain the control sequence; when τ ∈ [t k ,..., t k+N , it is ensured that the cost function gradually decreases. When τ > t k+N , it is necessary to design a terminal feedback controller to achieve the control objective as τ approaches ∞. In t k After solving the optimization problem at all times, we get t k+N u(t k+N ;t k ) and t k+N+1 Time x(t k+N+1 ;t k ) value is already near the reference point. At this time, the model is linearized to replace the original system model. The error model after linearization is shown in formula (10); Design the terminal feedback controller as: Among them, δ is the random variable introduced above, satisfying the mathematical expectation E(δ) = n, where n represents the data loss rate, and υ is obtained by solving equation (46). For details, refer to equation (10); substituting equation (37) into equation (10) and after rearrangement, we get:

8. The method according to claim 7, wherein The above-mentioned Step 4 includes: If there is no data loss caused by communication interruption between the leader UAV and the follower UAV at time t k Since the asynchronous communication method is adopted between the leader UAV and the follower UAV, the follower UAV can obtain the true predicted state sequence and the true predicted control sequence Then at time t k+1 the predicted trajectory is expressed as: Denote the true predictive control sequence obtained by the follower UAV \(i\) solving the optimization problem at time \(t\). Take the first element of the sequence as the optimal control input to act on the system; k Denote the assumed predictive control input of the follower UAV \(i\), which is used to generate the assumed predictive trajectory Solve the optimization problem; where includes and ​ If at time t k there is no data loss caused by communication interruption between the follower UAV and the neighbor node UAV, since they use synchronous communication, the follower UAV cannot obtain the true predicted state sequence of the neighbor node UAV and the true predicted control sequence Only a hypothetical predicted control sequence can be designed, which is expressed as: where, κ j (x j (·)) represents the terminal feedback controller, and its specific form is given in Equation (37); the above equation means that the assumed predicted state sequence of neighbor node UAV j received at time t k is composed of the control sequences at times t k-1 obtained by solving the optimization problem at time t k ,..., t k-2+N and the control quantity calculated by the terminal feedback controller at time t k-1+N . By assuming the predictive control sequence in Equation (21) and the system model in Equation (18), generate the assumed predictive state sequence, denoted as: Among them, is calculated from formula (21) κ j (x j (·)) and system model formula (18); If at time t k data loss occurs due to a communication interruption between the leader UAV and the follower UAV, the follower UAV can only obtain the assumed prediction control sequence of the leader UAV, which is expressed as: By assuming the predictive control sequence in Equation (23) and the system model in Equation (18), generate the assumed predictive state sequence, denoted as: If at time t k data loss occurs due to the loss of communication between the follower UAV and the neighbor node UAV, the follower UAV can only obtain the assumed prediction control sequence of the neighbor node UAV. According to the descriptions of the loss of independent random data δ in equations (11) and (12), the assumed prediction control sequence is expressed as: By assuming the predictive control sequence in Equation (25) and the system model in Equation (18), generate the assumed predictive state sequence, denoted as:

9. The method according to claim 8, wherein The above-mentioned Step 5 includes: The follower UAV i solves the following optimization problem at time t k as follows: Among them, Ω i is the control constraint domain, χ i is the state constraint domain, μ ∈ [0, 1) is a constant, R is the safety distance during the formation flight of the UAVs, C i (ε i (τ; t k )) is the compatibility constraint, including the collision constraint of the UAV formation flight and the limiting constraint related to the algorithm stability, p li (t k+N ; t k ), q li (t k+N ; t k ) and u li (t k+N ; t k ) represent the predicted values of the state and control variables at time t k for the state and control variables at time t k+N . p li (t k ), q li (t k ) and u li (t k ) represent the actual values of the state and control variables at time t k . Adding Δ ij as a soft constraint to perform constraint contraction on ; Solve the optimization problem of equation (27) at time t k to obtain the optimized control sequence as follows: In the above formula represents the optimal control sequence obtained after solving the optimization problem of formula (27); Take the optimized value of the control sequence t k at the moment as the control input of the system, that is: At each moment \(t\) k Solve the optimization problem in equation (27), and apply the optimal solution obtained from equation (29) to the system. Meanwhile, the prediction horizon of the algorithm rolls from \([t k ,t k+N \) to \([t k+1 ,t k+1+N \), and repeat the solution of the optimization problem in equation (27) at time \(t k+1 \) to achieve rolling optimization.

10. The method according to claim 9, characterized in that, When solving the optimization problem in Equation (27), it must be satisfied that: ||p i (τ;t k ) - p j (τ;t k )|| > 2R, τ ∈ [t k ,..., t k+N (30) Wherein, R is the safety distance during the formation flight of the UAVs; p j (τ; t k ) and there is an error ε j (τ; t k ), defined as: In the above formula is the estimated value of the position of neighbor node UAV j at time τ ∈ [t k ,..., t k+N , and p j (τ; t k ) is the actual value of the position of neighbor node UAV j at time τ ∈ [t k ,..., t k+N ; Due to the error ε j (τ; t k ), the constraints of equation (30) need to be contracted, expressed as: where R is the safety distance during the formation flight of the drones, Δ ij is expressed as: In solving the optimization problem, the error between the assumed predicted state sequence and the actual predicted state sequence is restricted, expressed as: ||ε j (τ;t k )||≤min{Δ i (t k ),η i (t k )} (34) The above equation is the compatibility constraint C in the optimization problem of equation (27). i (ε i (τ; t k )) is the expression form, where Δ i (t k ) and η i (t k ) are respectively expressed as: Among them, α ∈ [0, 1) is a constant related to the convergence speed, T is the step size between t k and t k+1 moments, N is the prediction horizon length, Q ij is the collaborative weight matrix, λ max (Q ij ) is the maximum eigenvalue of matrix Q ij , η i (t k ) is an index related to the algorithm stability.

Images