Calculation method of UAV flight safety envelope based on nonlinear programming and reachability analysis

Abstract

The present invention discloses a method for calculating the flight safety envelope of a UAV based on nonlinear programming and reachability analysis, comprising the following steps: 1. Establishing a corresponding nonlinear model according to the dynamic characteristics of the UAV, mainly including nonlinear state equations, state and input constraints, etc.; 2. Based on the reachable equilibrium set theory, using the alternating direction multiplier method to calculate the reachable equilibrium set of the UAV; 3. Taking the calculated reachable equilibrium set as the basic set, respectively calculating the forward reachable set and the backward reachable set of the basic set, and taking the intersection of the forward reachable set and the backward reachable set as the flight safety envelope of the UAV. Aiming at the problem of UAV flight safety, the present invention proposes a method for determining the flight safety boundary, determines the state restrictions during the flight of the UAV, and provides a necessary basis for the boundary protection control of the UAV.

Description

Technical Field

[0001] The invention relates to a method for calculating a flight safety envelope of an unmanned aerial vehicle, and in particular to a method for calculating a flight safety envelope of an unmanned aerial vehicle based on nonlinear programming and reachability analysis. Background Art

[0002] As the performance of drones improves, their maneuverability becomes stronger and stronger, and the resulting flight safety issues have also attracted people's attention. Flight safety refers to the occurrence of casualties or damage to the aircraft due to improper operation or external reasons during the operation of the aircraft. Unlike other traffic accidents, the probability of flight accidents is low but the destructiveness is huge. Looking at the various flight accidents that have occurred in history, it can be seen that different flight accidents occur at each stage of the aircraft's taxiing, takeoff, climb, cruising to the final landing, and the main inducing factors are also different, but they can be classified into driver factors, environmental factors, and system failures.

[0003] There is still no unified concept of loss of control and the cause is unknown. However, it is internationally recognized that a flight state that exceeds a certain flight safety boundary (i.e., the flight safety envelope) can be judged as loss of control, and loss of control does not necessarily lead to a flight accident, but if it is not detected or left unchecked, an accident will definitely occur.

[0004] In order to give full play to the maneuverability of UAVs, avoid flight accidents, and reveal the mechanism of loss of control, it is of great theoretical significance and engineering value to study the calculation of the flight safety boundary of UAVs. The traditional flight safety envelope refers to the speed and altitude level flight envelope, but even if the aircraft is within this envelope, it will still induce loss of control due to external disturbances. Therefore, scholars have given a new concept to the flight envelope. They believe that the flight envelope consists of an environmental envelope, a structural envelope, and a dynamic envelope. The dynamic envelope has attracted the attention of many scholars because it takes into account the dynamic characteristics of the aircraft and is complex to calculate. This is also the flight envelope that the present invention focuses on. Summary of the invention

[0005] Purpose of the invention: The purpose of the present invention is to provide a method for calculating the flight safety envelope of a UAV based on nonlinear programming and reachability analysis, so as to provide a necessary basis for flight safety boundary protection control.

[0006] Technical solution: The method for calculating the flight safety envelope of a UAV based on nonlinear programming and reachability analysis of the present invention comprises the following steps:

[0007] S1. Establish a nonlinear state equation model of the UAV based on its dynamic characteristics, and determine the value range of its state quantity and input quantity by simulation, wind tunnel experiment, flight experiment or theoretical calculation;

[0008] S2. Based on the reachable equilibrium set theory, the calculation of the reachable equilibrium set of the nonlinear state equation of the UAV is transformed into an equivalent optimization problem, and then the corresponding optimization problem is solved using the alternating multiplier method to obtain the reachable equilibrium set of the nonlinear state equation of the UAV;

[0009] S3. Taking the reachable equilibrium set calculated in step S2 as the basic set, based on the reachability analysis theory, respectively calculate the forward reachable set and the backward reachable set of the basic set, and then take their intersection as the flight safety envelope of the UAV.

[0010] Furthermore, the nonlinear state space model of the UAV in step S1 is:

[0011]

[0012] in, is the system state quantity, is the system input, is a nonlinear function, and They are constraints on state quantity and input quantity respectively, and their specific forms are as follows:

[0013]

[0014]

[0015] Among them, x i,min is the i-th state quantity x i The minimum value of x i,max is the i-th state quantity x i The maximum value of u, i=1,2,...,n j,min is the jth input quantity u j The minimum value of u j,max is the jth input quantity u j The maximum value of , j = 1, 2, ..., m.

[0016] Furthermore, in step S2, the calculation of the reachable equilibrium set of the nonlinear state space model of the UAV is transformed into an equivalent optimization problem as follows:

[0017] The reachable equilibrium set is expressed as:

[0018]

[0019] Among them, S ε is the reachable equilibrium set of the UAV; x ε ,u ε are the state quantity and input quantity at the equilibrium point respectively; f(·,·) is a nonlinear function, and are the constraints on state quantity and input quantity respectively.

[0020] Finding the reachable equilibrium set of the drone is the root-finding problem of the nonlinear equation under constraints, as shown in the following formula:

[0021]

[0022] In order to facilitate the handling of constraints, the above nonlinear equation root-finding problem is transformed into an equivalent optimization problem.

[0023] First, construct the objective function as follows:

[0024]

[0025] in, is the set of all n-order symmetric positive definite matrices, J(x,u) is the objective function, and f(x,u) is a nonlinear function. In this case, f(x,u)=0 if and only if J(x,u)=0.

[0026] Secondly, for constraint processing, the constraint types of state quantities and input quantities are linear inequality constraints as shown below:

[0027] x i -x i,max ≤0,x i,min -x i ≤0i=1,2,…,n

[0028] u j -u j,max ≤0,u j,min -u j ≤0j=1,2,…,m

[0029] Among them, x i,min is the i-th state quantity x i The minimum value of x i,max is the i-th state quantity x i The maximum value of u j,min is the jth input quantity u j The minimum value of u j,max is the jth input quantity u j The maximum value of .

[0030] Introduce X=(x T ,u T ) T represents the optimization variable, h(X)=J(x,u) represents the objective function, c i (X)≤0,i∈Ω, where Ω={1,2,…,2(n+m)} is the index set, representing the 2(n+m) linear inequality constraints in the above formula, c i(X)≤0,i∈Ω can be expressed in the following matrix form:

[0031] AX≤b

[0032] in, A=diag(A1,A2,…,A n+m ),A i =(1,-1) T ,i=1,2,…,(n+m), b j =(X j,max ,-X j,min ) T , j=1,2,…,(n+m).

[0033] For each linear inequality constraint c i (X)≤0,i∈Ω, introduce slack variable s i ≥0, so that the linear inequality constraint becomes a linear equality constraint:

[0034] c i (X)+s i =0,i∈Ω

[0035] Let S = (s1, s2, ..., s 2(n+m) ) T , I is a 2(n+m)-dimensional identity matrix, and the linear equality constraint above is expressed in matrix form:

[0036] AX+IS=b

[0037] At this point, the root-finding problem of the nonlinear equation is transformed into the following optimization problem:

[0038] min X,S h(X)

[0039] stAX+IS=b

[0040] s i ≥0,i∈Ω

[0041] Furthermore, in step S2, the corresponding optimization problem solved by using the alternating multiplier method is specifically:

[0042] definition is a first-order continuously differentiable function and satisfies: when S ≥ 0, φ(S) = 0, otherwise φ(S) > 0, definition for:

[0043] g(S)=Mφ(S)

[0044] Among them, M>0 is a large positive number. Adding g(S) to the objective function, we get the following optimization problem with equality constraint structure separation:

[0045] min X,S h(X)+g(S)

[0046] stAX+IS=b

[0047] make As the Lagrange multiplier of the equality constraint, construct the following Lagrange function:

[0048] L(x,y,λ)=f(X)+g(S)-λ T (AX+IS-b)

[0049] Solving the optimization problem of equality constraint structure separation is to solve the saddle point (X * ,S * ,λ * ), satisfying the following expression:

[0050]

[0051] The alternating direction multiplier method is used to solve the problem shown in the above formula. The solution process includes two stages: prediction and correction.

[0052] In the prediction stage, the k-th iteration generates X, S, and the Lagrange multiplier λ as X k , S k , and λ k , by the iterative value sequence v of the kth step k =(X k ,S k ,λ k ) Generate the predicted value of the k-th iteration sequence The specific implementation is as follows.

[0053] Given (S k ,λ k ), is the solution to the following optimization problem, we can get:

[0054]

[0055] Where β>0 is a constant.

[0056] Lagrange multiplier λ k The predicted value of is:

[0057]

[0058] Where β>0 is a constant.

[0059] According to the obtained To solve S k The predicted value of The specific expression is as follows:

[0060]

[0061] Where β>0 is a constant.

[0062] In the correction phase, according to the iteration value v generated by the k-th iteration k =(X k ,S k ,λ k ) and predicted values To generate the iterative value X of the k+1th step X, S, Lagrange multiplier λ k+1 , S k+1 , and λ k+1 And let v k+1 =(X k +1 ,S k+1 ,λ k+1 ), there is the following relationship:

[0063]

[0064] Where γ∈(0,2) is a constant.

[0065] Furthermore, step S3 is specifically as follows:

[0066] Use the polyhedron to calculate the reachable equilibrium set S in step S2 ε As an approximation, let K be the reachable equilibrium set S ε The polyhedron approximation of K can be expressed as follows:

[0067]

[0068] in, is t linear functions of the state x.

[0069] Given T f >0 is a time constant, and the forward reachable set and invariant set of the nonlinear state equation of the UAV are expressed as:

[0070]

[0071]

[0072] Among them, R(T f ,K) is the forward reachable set of the UAV, I(T f ,K) is the invariant set of the UAV, x(0) is the initial condition of the UAV, and are the constraints of state quantity and input quantity respectively, For the permissible input, τ is [0,Tf ] is a time point within the time period, x(τ) is the state of the drone at time τ.

[0073] The calculation of the invariant set and the forward reachable set is expressed by the following dual relation:

[0074] R(T f ,K)=(I(T f ,K c )) c

[0075] This formula shows that the forward reachable set can be determined by calculating the invariant set.

[0076] Nonlinear state equations for UAVs Another system that is reversed in time is:

[0077]

[0078] The reachable set of the nonlinear state equation of the UAV is called the forward reachable set R f (T f ,K), and the forward reachable set of the other system in the reverse direction in time is called the backward reachable set R b (T f ,K); The dynamic flight envelope of the UAV is defined as the intersection of the forward reachable set and the backward reachable set, as shown in the following formula:

[0079] E(T f ,K)=R f (T f ,K)∩R b (T f ,K)

[0080] Among them, E(T f , K) is the dynamic flight envelope.

[0081] From the dual relationship between the invariant set and the forward reachable set, we can determine the forward reachable set by calculating the invariant set of the drone. The following is an explanation of the calculation of the invariant set by the level set method. As mentioned above, the polyhedron approximation K of the reachable equilibrium set is:

[0082]

[0083] Then the invariant set is expressed in the form of level set as:

[0084]

[0085] where Φ(τ; x(0), u(·)) is the state trajectory of the system under given initial conditions x(0) and admissible control, and the level set function V(x, T f) can be expressed as the viscosity solution of the following Hamilton-Jacobi-Bellman partial differential equation:

[0086]

[0087] in, It is solved by numerical calculation method.

[0088] The present invention provides a calculation system for a UAV flight safety envelope based on nonlinear programming and reachability analysis, comprising:

[0089] The model building module is used to establish the nonlinear state equation of the UAV according to its dynamic characteristics, and determine the value range of its state quantity and input quantity by simulation, wind tunnel experiment, flight experiment or theoretical calculation;

[0090] The reachable equilibrium set solving module converts the calculation of the reachable equilibrium set into an equivalent optimization problem based on the reachable equilibrium set theory, and then uses the alternating multiplier method to solve the corresponding optimization problem;

[0091] The flight safety envelope determination module uses the calculated reachable equilibrium set as the basic set, calculates the forward reachable set and backward reachable set of the basic set respectively, and then takes their intersection as the flight safety envelope of the UAV.

[0092] A device of the present invention includes a memory and a processor, wherein:

[0093] A memory for storing computer programs that can be run on the processor;

[0094] The processor is used to execute the steps of the above-mentioned method for calculating the flight safety envelope of a drone based on nonlinear programming and reachability analysis when running the computer program.

[0095] A storage medium of the present invention stores a computer program, and when the computer program is executed by at least one processor, the steps of the method for calculating the flight safety envelope of a drone based on nonlinear programming and reachability analysis are implemented.

[0096] Beneficial effect: Compared with the traditional method of calculating the UAV flight safety envelope, the present invention transforms the solution of the reachable equilibrium set from the root-finding problem of nonlinear equations under constraints into an optimization problem under constraints based on the optimization method. With the help of the alternating direction multiplier method, the optimization variables can be grouped and iterated. Only part of the optimization variables are calculated in each iteration, which can effectively reduce the amount of calculation in the iterative process, thereby improving the calculation efficiency and reducing the time for calculating the UAV flight safety envelope. BRIEF DESCRIPTION OF THE DRAWINGS

[0097] Figure 1It is a flow chart of the method of the present invention. DETAILED DESCRIPTION

[0098] The following will describe the implementation methods of the present invention in detail in combination with the embodiments, so as to fully understand and implement the implementation process of how the present invention applies technical means to solve technical problems and achieve technical effects. It should be noted that as long as there is no conflict, the various embodiments of the present invention and the various features in the embodiments can be combined with each other, and the technical solutions formed are all within the protection scope of the present invention.

[0099] Aiming at the problem of UAV flight safety, the present invention proposes a method for determining the flight safety boundary to solve the problems that the existing technology is single and cannot effectively handle constraints, and determines the state limit of the UAV during flight, which provides the necessary basis for the boundary protection control of the UAV. Figure 1 As shown, a method for calculating the flight safety envelope of a UAV based on nonlinear programming and reachability analysis of the present invention comprises the following steps:

[0100] Step 1: Establish a nonlinear state equation model of the UAV according to its dynamic characteristics, and determine the value range of its state quantity and input quantity by simulation, wind tunnel experiment, flight experiment or theoretical calculation;

[0101] In order to facilitate modeling and subsequent calculations, the following assumptions are introduced:

[0102] 1) The drone is a rigid body;

[0103] 2) The earth is a flat surface, ignoring the curvature and rotation of the earth;

[0104] 3) The atmosphere is a static standard atmosphere, without considering wind interference, etc.;

[0105] 4) Ignore the fuel consumption during the UAV maneuvering process, that is, the mass of the UAV is constant;

[0106] 5) Each state quantity and input quantity of the drone is constrained between two constant values;

[0107] 6) The maneuver constraint type of the UAV is linear or can be linearized.

[0108] Based on the above assumptions, the general form of the UAV nonlinear state equation model can be established as:

[0109]

[0110] in, is the system state quantity, is the system input, is a nonlinear function, and They are constraints on state quantity and input quantity respectively, and their specific forms are as follows:

[0111]

[0112] Among them, x i,min is the i-th state quantity x i The minimum value of x i,max is the i-th state quantity x i The maximum value of u, i=1,2,...,n j,min is the jth input quantity u j The minimum value of u j,max is the jth input quantity u j The maximum value of j = 1, 2, ..., m. According to the constraint of formula (2), the flight safety envelope of the nonlinear state equation of the UAV of formula (1) is calculated. There are many different understandings and definitions of the flight safety envelope of the UAV. The present invention mainly focuses on the dynamic or dynamic envelope.

[0113] Step 2: Based on the reachable equilibrium set theory, the calculation of the reachable equilibrium set is transformed into an equivalent optimization problem, and then the alternating multiplier method is used to solve the corresponding optimization problem;

[0114] The concept of reachable equilibrium set refers to the set of all equilibrium points that the drone can reach under the available bounded input, that is, under the control of the input Under the restriction, the state quantity The set of all balancing points in can be expressed as:

[0115]

[0116] Among them, S ε is the reachable equilibrium set of the UAV; x ε ,u ε are the state quantity and input quantity at the equilibrium point respectively; f(·,·) is a nonlinear function, and are the constraints of state quantity and input quantity respectively. Finding the reachable equilibrium set of the UAV is the root-finding problem of the nonlinear equation under the constraints, as shown in the following formula:

[0117]

[0118] Among them, x ε ,u ε are the state quantity and input quantity at the equilibrium point respectively; f(·,·) is a nonlinear function, and are the constraints on state quantity and input quantity respectively.

[0119] In order to facilitate the handling of constraints, the root-finding problem of the nonlinear equation shown in (4) is transformed into an equivalent optimization problem.

[0120] First, construct the objective function as follows:

[0121]

[0122] in, is the set of all n-order symmetric positive definite matrices, J(x,u) is the objective function, and f(x,u) is a nonlinear function. In this case, f(x,u)=0 if and only if J(x,u)=0.

[0123] Next, the constraint processing is performed. The constraint type of formula (2) is a linear inequality constraint as shown below:

[0124]

[0125] Among them, x i,min is the i-th state quantity x i The minimum value of x i,max is the i-th state quantity x i The maximum value of u j,min is the jth input quantity u j The minimum value of u j,max is the jth input quantity u j The maximum value of .

[0126] Introduce X=(x T ,u T ) T represents the optimization variable, h(X)=J(x,u) represents the objective function, c i (X)≤0,i∈Ω, where Ω={1,2,…,2(n+m)} is the index set, representing the 2(n+m) linear inequality constraints of formula (6). Formula (6) can be expressed in the following matrix form:

[0127] AX≤b (7)

[0128] in, A=diag(A1,A2,…,A n+m ),A i =(1,-1) T ,i=1,2,…,(n+m), b j =(X j,max ,-X j,min ) T , j=1,2,…,(n+m).

[0129] For each linear inequality constraint c i (X)≤0,i∈Ω, introduce slack variable s i ≥0, so that the linear inequality constraint becomes a linear equality constraint:

[0130] c i (X)+s i =0,i∈Ω (8)

[0131] Let S = (s1, s2, ..., s 2(n+m) ) T , I is a 2(n+m)-dimensional identity matrix, and equation (8) can be expressed in the following matrix form:

[0132] AX+IS=b (9)

[0133] At this point, the root-finding problem of the nonlinear equation (4) is transformed into the following optimization problem:

[0134]

[0135] The variables and their definitions refer to the above.

[0136] definition is a first-order continuously differentiable function and satisfies: when S ≥ 0, φ(S) = 0, otherwise φ(S) > 0, definition for:

[0137] g(S)=Mφ(S) (11)

[0138] Among them, M>0 is a large positive number. Adding g(S) to the objective function, we get the following optimization problem with equality constraint structure separation:

[0139]

[0140] This form of optimization problem is a typical form that can be solved by the alternating direction multiplier method. As the Lagrange multiplier of the equality constraint, construct the following Lagrange function:

[0141] L(x,y,λ)=f(X)+g(S)-λ T (AX+IS-b) (13)

[0142] Solving problem (12) is to find the saddle point (X * ,S * ,λ * ), satisfying the following expression:

[0143]

[0144] The alternating direction multiplier method based on customized PPA can solve the problem shown in (14). The solution process includes two stages: prediction and correction. In the prediction stage, X, S, and Lagrange multiplier λ generated by the k-th iteration are represented as X k, S k , and λ k , by the iterative value sequence v of the kth step k =(X k ,S k ,λ k ) Generate the predicted value of the k-th iteration sequence The specific implementation is as follows. Given (S k ,λ k ), is the solution to the following optimization problem, we can get:

[0145]

[0146] Where β>0 is a constant.

[0147] The predicted values ​​of the Lagrange multipliers are:

[0148]

[0149] Where β>0 is a constant.

[0150] According to the obtained To obtain The specific expression is as follows:

[0151]

[0152] Where β>0 is a constant.

[0153] In the correction phase, according to the iteration value v generated by the k-th iteration k =(X k ,S k ,λ k ) and predicted values To generate the iterative value X of the k+1th step X, S, Lagrange multiplier λ k+1 , S k+1 , and λ k+1 And let v k+1 =(X k+1 ,S k+1 ,λ k+1 ), there is the following relationship:

[0154]

[0155] Where γ∈(0,2) is a constant.

[0156] Step 3: Take the reachable equilibrium set calculated in step 2 as the basic set, calculate the forward reachable set and backward reachable set of the basic set respectively, and then take their intersection as the flight safety envelope of the UAV;

[0157] Reachability analysis refers to the process of determining whether an effective control strategy can be used to make the state within the initial set reach the specified target set within a specific time range for a differential dynamic system. This method has outstanding performance in safety analysis in the aviation field and is widely used. In general, reachability analysis studies the reachability problem of differential dynamic systems through positive / reverse reachable sets to solve aviation safety problems in different situations.

[0158] Consider the UAV system, whose dynamic equation is: The variables and their meanings are the same as above. S calculated in step 2 ε For the convenience of reachability analysis, we use polyhedrons to represent S ε As an approximation, let K be the corresponding polyhedron, T f >0 is a given constant value, and the reachable set and invariant set of the nonlinear state equation of the UAV can be expressed as:

[0159]

[0160]

[0161] Among them, R(T f ,K) is the forward reachable set of the UAV, I(T f ,K) is the invariant set of the UAV, x(0) is the initial condition of the UAV, and are the constraints of state quantity and input quantity respectively, For the permissible input, τ is [0,T f ] is a time point within the time period, x(τ) is the state of the drone at time τ.

[0162] The calculation of invariant sets and reachable sets can be related to each other through the following relationship, that is, the duality principle, which can be expressed as:

[0163] R(T f ,K)=(I(T f ,K c )) c (twenty one)

[0164] A similar definition can be used for another system that is reversed in time as follows:

[0165]

[0166] The reachable set of the nonlinear state equation of the UAV is called the forward reachable set R f (T f ,K), and the forward reachable set of the other system in the reverse direction in time is called the backward reachable set R b (Tf ,K); The dynamic flight envelope of the UAV is defined as the intersection of the forward reachable set and the backward reachable set, as shown in the following formula:

[0167] E(T f ,K)=R f (T f ,K)∩R b (T f ,K) (23)

[0168] From the dual relationship between the invariant set and the forward reachable set, we can determine the forward reachable set by calculating the invariant set of the drone. The following is an explanation of the calculation of the invariant set by the level set method. As mentioned above, the polyhedron approximation K of the reachable equilibrium set is:

[0169]

[0170] Then the invariant set can be expressed in level set form as:

[0171]

[0172] where Φ(τ; x(0), u(·)) is the state trajectory of the system under given initial conditions x(0) and admissible control, and the level set function V(x, T f ) can be expressed as the viscosity solution of the following Hamilton-Jacobi-Bellman partial differential equation:

[0173]

[0174] in, It can be solved with the help of Matlab's level set toolbox. After calculating the invariant set, the forward reachable set and the backward reachable set can be calculated respectively by the duality principle, and their intersection is taken as the flight safety envelope of the drone.

Claims

1. A method for calculating the flight safety envelope of a UAV based on nonlinear programming and reachability analysis, characterized in that: The following steps are involved: S1. According to the dynamic characteristics of the UAV, a nonlinear state equation model is established, and the value range of its state quantity and input quantity is determined by simulation, wind tunnel experiment, flight experiment or theoretical calculation; the nonlinear state equation model of the UAV is: in, is the system state quantity, is the system input, is a nonlinear function, and They are constraints on state quantity and input quantity respectively, and their specific forms are as follows: Among them, x i,min is the i-th state quantity x i The minimum value of x i,max is the i-th state quantity x i The maximum value of, i=1,2,...,n,u j,min is the jth input quantity u j The minimum value of u j,max is the jth input quantity u j The maximum value of , j = 1, 2, ..., m; S2. Based on the reachable equilibrium set theory, the calculation of the reachable equilibrium set of the nonlinear state equation of the UAV is transformed into an equivalent optimization problem, and then the corresponding optimization problem is solved by the alternating multiplier method to obtain the reachable equilibrium set of the nonlinear state equation of the UAV; the calculation of the reachable equilibrium set of the nonlinear state equation of the UAV is transformed into an equivalent optimization problem, specifically: The reachable equilibrium set is expressed as: Among them, S ε is the reachable equilibrium set of the UAV, x ε ,u ε are the state quantity and input quantity at the equilibrium point respectively; Finding the reachable equilibrium set of the drone is the root-finding problem of the nonlinear equation under constraints, as shown in the following formula: In order to facilitate the handling of constraints, the above nonlinear equation root-finding problem is transformed into an equivalent optimization problem; First, construct the objective function as follows: in, is the set of all n-order symmetric positive definite matrices, J(x,u) is the objective function, f(x,u) is a nonlinear function, and f(x,u)=0 if and only if J(x,u)=0; Secondly, for constraint processing, the constraint types of state quantities and input quantities are linear inequality constraints as shown below: x i -x i,max ≤0,x i,min -x i ≤0i=1,2,…,n u j -u j,max ≤0,u j,min -u j ≤0j=1,2,…,m Introduce X=(x T ,u T ) T represents the optimization variable, h(X)=J(x,u) represents the objective function, c i (X)≤0,i∈Ω, where Ω={1,2,…,2(n+m)} is the index set, representing the 2(n+m) linear inequality constraints in the above formula, c i (X)≤0,i∈Ω is expressed in the following matrix form: AX≤b in, A=diag(A1,A2,…,A n+m ),A i =(1,-1) T ,i=1,2,…,(n+m), b j =(X j,max ,-X j,min ) T , j = 1, 2, ..., (n + m); For each linear inequality constraint c i (X)≤0,i∈Ω, introduce slack variable s i ≥0, so that the linear inequality constraint becomes a linear equality constraint: c i (X)+s i =0,i∈Ω Let S = (s1, s2, ..., s 2(n+m) ) T , I is a 2(n+m)-dimensional identity matrix, and the linear equality constraint above is expressed in matrix form: AX+IS=b At this point, the root-finding problem of the nonlinear equation is transformed into the following optimization problem: my X,S h(X) stAX+IS=b s i ≥0,i∈Ω The corresponding optimization problem solved by the alternating multiplier method is as follows: definition is a first-order continuously differentiable function and satisfies: when S ≥ 0, φ(S) = 0, otherwise φ(S) > 0, definition for: g(S)=Mφ(S) Among them, M>0 is a large positive number. Adding g(S) to the objective function, we get the following optimization problem with equality constraint structure separation: min X,S h(X)+g(S) stAX+IS=b make As the Lagrange multiplier of the equality constraint, construct the following Lagrange function: L(x,y,λ)=f(X)+g(S)-λ T (AX+IS-b) Solving the optimization problem of equality constraint structure separation is to solve the saddle point (X * ,S * ,λ * ), satisfying the following expression: The alternating direction multiplier method is used to solve the problem shown in the above formula. The solution process includes two stages: prediction and correction. In the prediction stage, the k-th iteration generates X, S, and the Lagrange multiplier λ as X k , S k , and λ k , by the iterative value sequence v of the kth step k =(X k ,S k ,λ k ) Generate the predicted value of the k-th iteration sequence The specific implementation is as follows: Given (S k ,λ k ), is the solution to the following optimization problem, we can get: Where β>0 is a constant; Lagrange multiplier λ k The predicted value of is: According to the obtained To solve S k The predicted value of The specific expression is as follows: In the correction phase, according to the iteration value v generated by the k-th iteration k =(X k ,S k ,λ k ) and predicted values To generate the iterative value X of the k+1th step X, S, Lagrange multiplier λ k+1 , S k+1 , and λ k+1 And let v k+1 =(X k+1 ,S k+1 ,λ k+1 ), there is the following relationship: Where γ∈(0,2) is a constant; S3, taking the reachable equilibrium set calculated in step S2 as the basic set, and based on the reachability analysis theory, respectively calculating the forward reachable set and the backward reachable set of the basic set, and taking their intersection as the flight safety envelope of the UAV; Specifically: Use the polyhedron to calculate the reachable equilibrium set S in step S2 ε As an approximation, let K be the reachable equilibrium set S ε The polyhedron approximation of K is expressed as follows: in, is t linear functions about the state x; Given T f >0 is a time constant, and the forward reachable set and invariant set of the nonlinear state equation of the UAV are expressed as: Among them, R(T f ,K) is the forward reachable set of the UAV, I(T f ,K) is the invariant set of the UAV, x(0) is the initial condition of the UAV, For the permissible input, τ is [0,T f ] is a time point in the time period, x(τ) is the state of the drone at time τ; The calculation of the invariant set and the forward reachable set is expressed by the following dual relation: R(T f ,K)=(I(T f ,K c )) c Nonlinear state equations for UAVs Another system that is reversed in time is: The reachable set of the nonlinear state equation of the UAV is called the forward reachable set R f (T f ,K), and the forward reachable set of the other system in the reverse direction in time is called the backward reachable set R b (T f ,K); The dynamic flight envelope of the UAV is defined as the intersection of the forward reachable set and the backward reachable set, as shown in the following formula: E(T f ,K)=R f (T f ,K)∩R b (T f ,K) Among them, E(T f , K) is the dynamic flight envelope.

2. The method for calculating the flight safety envelope of an unmanned aerial vehicle based on nonlinear programming and reachability analysis according to claim 1, characterized in that: The forward reachable set is determined by calculating the invariant set of the drone; the invariant set is calculated by the level set method, and the polyhedron approximation K of the reachable equilibrium set is: Then the invariant set is expressed in the form of level set as: where Φ(τ; x(0), u(·)) is the state trajectory of the system under given initial conditions x(0) and admissible control, and the level set function V(x, T f ) is expressed as the viscosity solution of the following Hamilton-Jacobi-Bellman partial differential equation: in, It is solved by numerical calculation method.

3. A calculation system for the UAV flight safety envelope based on nonlinear programming and reachability analysis, characterized in that: include: The model building module is used to establish the nonlinear state equation of the UAV according to its dynamic characteristics, and determine the value range of its state quantity and input quantity by simulation, wind tunnel experiment, flight experiment or theoretical calculation; the nonlinear state equation model of the UAV is: in, is the system state quantity, is the system input, is a nonlinear function, and They are constraints on state quantity and input quantity respectively, and their specific forms are as follows: Among them, x i,min is the i-th state quantity x i The minimum value of x i,max is the i-th state quantity x i The maximum value of, i=1,2,...,n,u j,min is the jth input quantity u j The minimum value of u j,max is the jth input quantity u j The maximum value of , j = 1, 2, ..., m; The reachable equilibrium set solving module, based on the reachable equilibrium set theory, transforms the calculation of the reachable equilibrium set into an equivalent optimization problem, and then uses the alternating multiplier method to solve the corresponding optimization problem to obtain the reachable equilibrium set; transforms the calculation of the reachable equilibrium set of the nonlinear state equation of the UAV into an equivalent optimization problem, specifically: The reachable equilibrium set is expressed as: Among them, S ε is the reachable equilibrium set of the UAV, x ε ,u ε are the state quantity and input quantity at the equilibrium point respectively; Finding the reachable equilibrium set of the drone is the root-finding problem of the nonlinear equation under constraints, as shown in the following formula: In order to facilitate the handling of constraints, the above nonlinear equation root-finding problem is transformed into an equivalent optimization problem; First, construct the objective function as follows: in, is the set of all n-order symmetric positive definite matrices, J(x,u) is the objective function, f(x,u) is a nonlinear function, and f(x,u)=0 if and only if J(x,u)=0; Secondly, for constraint processing, the constraint types of state quantities and input quantities are linear inequality constraints as shown below: x i -x i,max ≤0,x i,min -x i ≤0i=1,2,…,n u j -u j,max ≤0,u j,min -u j ≤0j=1,2,…,m Introduce X=(x T ,u T ) T represents the optimization variable, h(X)=J(x,u) represents the objective function, c i (X)≤0,i∈Ω, where Ω={1,2,…,2(n+m)} is the index set, representing the 2(n+m) linear inequality constraints in the above formula, c i (X)≤0,i∈Ω is expressed in the following matrix form: AX≤b in, A=diag(A1,A2,…,A n+m ),A i =(1,-1) T ,i=1,2,…,(n+m), b j =(X j,max ,-X j,min ) T , j = 1, 2, ..., (n + m); For each linear inequality constraint c i (X)≤0,i∈Ω, introduce slack variable s i ≥0, so that the linear inequality constraint becomes a linear equality constraint: c i (X)+s i =0,i∈Ω Let S = (s1, s2, ..., s 2(n+m) ) T , I is a 2(n+m)-dimensional identity matrix, and the linear equality constraint above is expressed in matrix form: AX+IS=b At this point, the root-finding problem of the nonlinear equation is transformed into the following optimization problem: my X,S h(X) stAX+IS=b s i ≥0,i∈Ω The corresponding optimization problem solved by the alternating multiplier method is as follows: definition is a first-order continuously differentiable function and satisfies: when S ≥ 0, φ(S) = 0, otherwise φ(S) > 0, definition for: g(S)=Mφ(S) Among them, M>0 is a large positive number. Adding g(S) to the objective function, we get the following optimization problem with equality constraint structure separation: min X,S h(X)+g(S) stAX+IS=b make As the Lagrange multiplier of the equality constraint, construct the following Lagrange function: L(x,y,λ)=f(X)+g(S)-λ T (AX+IS-b) Solving the optimization problem of equality constraint structure separation is to solve the saddle point (X * ,S * ,λ * ), satisfying the following expression: The alternating direction multiplier method is used to solve the problem shown in the above formula. The solution process includes two stages: prediction and correction. In the prediction stage, the k-th iteration generates X, S, and the Lagrange multiplier λ as X k , S k , and λ k , by the iterative value sequence v of the kth step k =(X k ,S k ,λ k ) Generate the predicted value of the k-th iteration sequence The specific implementation is as follows: Given (S k ,λ k ), is the solution to the following optimization problem, we can get: Where β>0 is a constant; Lagrange multiplier λ k The predicted value of is: According to the obtained To solve S k The predicted value of The specific expression is as follows: In the correction phase, according to the iteration value v generated by the k-th iteration k =(X k ,S k ,λ k ) and predicted values To generate the iterative value X of the k+1th step X, S, Lagrange multiplier λ k+1 , S k+1 , and λ k+1 And let v k+1 =(X k+1 ,S k+1 ,λ k+1 ), there is the following relationship: Where γ∈(0,2) is a constant; The flight safety envelope determination module uses the calculated reachable equilibrium set as the basic set, calculates the forward reachable set and the backward reachable set of the basic set respectively, and then takes their intersection as the flight safety envelope of the UAV; specifically: Use the polyhedron to calculate the reachable equilibrium set S in step S2 ε As an approximation, let K be the reachable equilibrium set S ε The polyhedron approximation of K is expressed as follows: in, is t linear functions about the state x; Given T f >0 is a time constant, and the forward reachable set and invariant set of the nonlinear state equation of the UAV are expressed as: Among them, R(T f ,K) is the forward reachable set of the UAV, I(T f ,K) is the invariant set of the UAV, x(0) is the initial condition of the UAV, For the permissible input, τ is [0,T f ] is a time point in the time period, x(τ) is the state of the drone at time τ; The calculation of the invariant set and the forward reachable set is expressed by the following dual relation: R(T f ,K)=(I(T f ,K c )) c Nonlinear state equations for UAVs Another system that is reversed in time is: The reachable set of the nonlinear state equation of the UAV is called the forward reachable set R f (T f ,K), and the forward reachable set of the other system in the reverse direction in time is called the backward reachable set R b (T f ,K); The dynamic flight envelope of the UAV is defined as the intersection of the forward reachable set and the backward reachable set, as shown in the following formula: E(T f ,K)=R f (T f ,K)∩R b (T f ,K) Among them, E(T f , K) is the dynamic flight envelope.

4. A device, characterized in that: comprising a memory and a processor, wherein: A memory for storing computer programs that can be run on the processor; A processor is used to execute the steps of the method for calculating the flight safety envelope of a drone based on nonlinear programming and reachability analysis as described in any one of claims 1-2 when running the computer program.

5. A storage medium, characterized in that: The storage medium stores a computer program, which, when executed by at least one processor, implements the steps of the method for calculating the flight safety envelope of a drone based on nonlinear programming and reachability analysis as described in any one of claims 1-2.

Images