Many-to-one collaborative interception rapid trajectory optimization method

Abstract

The invention relates to a many-to-one collaborative interception rapid trajectory optimization method, in particular to a many-to-one collaborative interception rapid trajectory optimization method. The objective of the invention is to solve the problems of unstable sight and overload saturation in the existing proportional guidance law method; and the pseudo-spectral method is relatively long in calculation time and cannot be applied online. The method comprises the following steps: establishing a three-dimensional dynamic model of a guidance section in the interception missile, and setting boundary constraints and process constraints of the interception missile; establishing a multi-interception missile time collaborative trajectory optimization original problem P0; representing the interception missile speed, the trajectory inclination angle, the trajectory deflection angle and the control quantity in the kinetic model by using a launching coordinate system; a Bezier function is used for representing the three-dimensional coordinates of the interception projectile under the launching coordinate system; converting the original problem P0 into a parameter optimization problem P1; and calling a nonlinear programming solver to solve all control points of the Bezier function, and solving the values of all state variables and control variables according to all control points of the Bezier function. The method is applied to the field of trajectory optimization.

Description

Technical Field

[0001] The invention relates to a fast trajectory optimization method for many-to-one collaborative interception. Background Art

[0002] The airspace between traditional aircraft and satellite orbits at 20 to 100 km from the ground is called near-space. In recent years, with the continuous improvement of science and technology, the performance of near-space aircraft has been continuously upgraded, especially many hypersonic aircraft with flight Mach numbers of 5 to 20, with fast flight speed, large maneuvering range, and variable maneuvering styles. This enables them to break through almost all current defense systems, which is a great challenge for the current interception system. Therefore, it is necessary to design an interception method for near-space hypersonic flight targets. When the flight speed and maneuverability of the interceptor missile have no obvious advantages over the hypersonic target, one interceptor missile cannot guarantee a successful interception. Multiple interceptors must be used to coordinate operations to complement their functions, implement multi-level and all-round strikes on the target, and improve the interception success rate.

[0003] Traditional methods for designing multi-missile coordinated interception trajectories mostly use proportional guidance laws, which have the advantages of simple structure and easy implementation. However, proportional guidance laws are mainly based on the design of guidance instructions based on relative kinematic models, without considering constraints such as angle of attack and overload. They are only suitable for intercepting low-speed targets, and may have problems such as unstable line of sight and overload saturation for highly maneuvering targets. Numerical optimization methods such as pseudo-spectral methods have become a means to solve the above problems due to their advantages of being able to use nonlinear models and handle complex constraints. However, pseudo-spectral methods are time-consuming to calculate and cannot be applied online. Summary of the invention

[0004] The purpose of the present invention is to solve the problems of unstable line of sight and overload saturation in the existing proportional guidance law method, as well as the long time-consuming pseudo-spectral calculation and inability to be applied online, and to propose a many-to-one collaborative interception rapid trajectory optimization method.

[0005] A method for fast trajectory optimization of many-to-one coordinated interception has the following specific process:

[0006] Step 1: Establish a three-dimensional dynamic model of the interceptor missile's mid-range guidance section and set boundary constraints and process constraints of the interceptor missile;

[0007] The maximum sum of terminal velocities of multiple interceptors is taken as the performance index, and the original problem P0 of multi-interceptor time coordinated trajectory optimization is established;

[0008] Step 2: Use x, y, z in the launch coordinate system to represent the interceptor missile velocity v, trajectory inclination angle θ, trajectory deviation angle ψ, and control quantities α, σ in the three-dimensional dynamic model of the interceptor missile's guidance section;

[0009] Step 3: Use Bezier function to characterize the three-dimensional coordinates x, y, and z of the interceptor missile in the launch coordinate system;

[0010] Find the first two and last two control points of the Bezier function according to the boundary conditions;

[0011] Transform the original problem P0 into a parameter optimization problem P1;

[0012] Step 4: Discretize and initialize the parameters of the Bezier functions of the first two and last two control points in the parameter optimization problem P1;

[0013] Step 5: Call the nonlinear programming solver to find all the control points of the Bezier function, and find the values ​​of all state quantities and control quantities based on all the control points of the Bezier function.

[0014] The beneficial effects of the present invention are:

[0015] The present invention designs a multi-interceptor time-coordinated trajectory optimization method based on Bezier curves to make up for the shortcomings of proportional guidance law and pseudo-spectral method; the method first uses a nonlinear dynamic model to model the multi-interceptor time-coordinated trajectory optimization problem, and then characterizes the flight trajectory of the interceptor as a high-order Bezier curve, that is, converts the infinite-dimensional trajectory optimization problem into a finite parameter optimization problem, so that it can be solved quickly and accurately. The present invention can handle the complex constraints in the flight process of the interceptor, solves the problem that the proportional guidance law easily causes line of sight instability and overload saturation, and because fewer Bezier curve parameters are used to replace infinite-dimensional trajectory information, the number of optimization variables is reduced, and the optimization efficiency is improved compared with the pseudo-spectral method. The present invention has the potential for online application.

[0016] The present invention proposes an efficient method for optimizing the trajectory of the mid-guidance section of an interceptor missile. Through the method, the three-dimensional coordinates of the interceptor missile are characterized by Bezier functions, so that the original infinite-dimensional trajectory optimization problem is transformed into a finite-dimensional parameter optimization problem, the number of optimization variables is greatly reduced, and the optimization efficiency is effectively improved. In addition, the present invention also considers allowing multiple interceptors to collaboratively intercept the same target, and sets the performance index to the sum of the terminal velocities of all interceptors. Compared with each interceptor missile intercepting the target separately, the present invention can increase the terminal velocity of the interceptor missile relatively, and effectively improve the strike capability. Compared with the existing trajectory optimization technology, the multi-interceptor missile time collaborative trajectory optimization method based on Bezier functions can not only handle complex constraints that traditional guidance laws cannot handle, but also greatly improve the computational efficiency compared to other numerical optimization methods, and has the potential for online application. BRIEF DESCRIPTION OF THE DRAWINGS

[0017] Figure 1 It is a flow chart of the present invention;

[0018] Figure 2 It is the Bezier collaborative interception three-dimensional trajectory diagram;

[0019] Figure 3 It is the interception missile speed curve;

[0020] Figure 4 It is the interceptor missile attack angle curve diagram;

[0021] Figure 5 This is the interceptor missile overload curve diagram;

[0022] Figure 6 It is the interception bomb dynamic pressure curve diagram;

[0023] Figure 7 This is the heat flux density curve of the interceptor missile;

[0024] Figure 8 This is the three-dimensional trajectory diagram of the collaborative interception of Bezier and pseudospectral methods. DETAILED DESCRIPTION

[0025] Specific implementation method 1: The specific process of a multi-to-one collaborative interception rapid trajectory optimization method in this implementation method is as follows:

[0026] Step 1: Establish a three-dimensional dynamic model of the interceptor missile's mid-range guidance section and set boundary constraints and process constraints of the interceptor missile;

[0027] The maximum sum of terminal velocities of multiple interceptors is taken as the performance index, and the original problem P0 of multi-interceptor time coordinated trajectory optimization is established;

[0028] Step 2: Use x, y, z in the launch coordinate system to represent the interceptor missile velocity v, trajectory inclination angle θ, trajectory deviation angle ψ, and control quantities α, σ in the three-dimensional dynamic model of the interceptor missile's guidance section;

[0029] Step 3: Use Bezier function to characterize the three-dimensional coordinates x, y, and z of the interceptor missile in the launch coordinate system;

[0030] Find the first two and last two control points of the Bezier function according to the boundary conditions;

[0031] The original problem P0 is transformed into a parameter optimization problem P1, where the parameters to be optimized are the remaining control points of the Bezier function;

[0032] Step 4: Discretize and initialize the parameters of the Bezier functions of the first two and last two control points in the parameter optimization problem P1;

[0033] Step 5: Call the nonlinear programming solver to find all the control points of the Bezier function, and find the values ​​of all state quantities and control quantities based on all the control points of the Bezier function.

[0034] The many-to-one number is often greater than or equal to 2.

[0035] Specific implementation method 2: This implementation method is different from the specific implementation method 1 in that in step 1, a three-dimensional dynamic model of the guidance section of the interceptor missile is established, and boundary constraints and process constraints of the interceptor missile are set;

[0036] The maximum sum of terminal velocities of multiple interceptors is taken as the performance index, and the original problem P0 of multi-interceptor time coordinated trajectory optimization is established;

[0037] The specific process is:

[0038] In step 1, it is necessary to establish the dynamics model of the interceptor missile and analyze the constraints it is subject to during the interception process;

[0039] Step 1: Establish a three-dimensional dynamic model of the guidance section of the interceptor missile; the specific process is:

[0040] Treat the interception target as a fixed point;

[0041] Ignoring the influence of the earth's curvature and rotation, the interceptor missile adopts unpowered gliding in the mid-range guidance stage. The simplified three-dimensional dynamic model of the interceptor missile in the mid-range guidance stage in the launch coordinate system is described as:

[0042]

[0043] In the formula, x, y, z are the coordinate values ​​of the interceptor missile in the launch coordinate system;

[0044] v is the interceptor missile speed; θ is the trajectory inclination angle; ψ is the trajectory deviation angle;

[0045] g is the acceleration of gravity, σ is the roll angle, and m is the mass of the interceptor missile;

[0046] are the first-order derivatives of x, y, and z respectively;

[0047] is the first-order derivative of v; is the first-order derivative of θ; is the first-order derivative of ψ;

[0048] L and D are aerodynamic lift and aerodynamic drag, respectively;

[0049] The expressions of L and D are:

[0050]

[0051] Among them, ρ is the atmospheric density, S is the area of ​​the interceptor missile, and C L ,C Dare the aerodynamic lift coefficient and drag coefficient respectively. Their values ​​are only related to the flight Mach number and angle of attack. In practical applications, they can be obtained by interpolation of aerodynamic tables or function fitting. Here they are simplified into the expression of the angle of attack α;

[0052] C L ,C D The expression is:

[0053]

[0054] Where α is the angle of attack;

[0055] The expression of atmospheric density ρ is:

[0056]

[0057] Where, ρ0=1.226kg / m 3 is the atmospheric density at sea level, h ref =7254.24m is the reference altitude, y is the coordinate value of the interceptor missile on the y-axis of the launch coordinate system;

[0058] The expression of gravitational acceleration g is:

[0059] g=g0(R e / r) 2 (5)

[0060] Where g0 is the gravitational acceleration at sea level, g0 = 9.81 m / s 2 ;

[0061] R e is the radius of the Earth, R e =6378.14km;

[0062] r represents the distance from the interceptor to the center of the earth, r = R e +y;

[0063] Step 1 and 2: Set the boundary constraints and process constraints of the interceptor missile;

[0064] Step 13: Based on steps 11 and 12, the maximum sum of the terminal velocities of multiple interceptors is used as the performance indicator, and the original problem P0 of multi-interceptor time coordinated trajectory optimization is established.

[0065] The other steps and parameters are the same as those in the first embodiment.

[0066] Specific implementation method three: This implementation method is different from specific implementation methods one or two in that the boundary constraints and process constraints of the interceptor missile are set in steps one and two; the specific process is:

[0067] Step 121: Take the coordinates x, y, z of the interceptor missile in the launch coordinate system, the interceptor missile speed v, the ballistic inclination angle θ, and the ballistic deviation angle ψ as the state variable X = [xyzv θ ψ] T , the superscript T means to find the transpose;

[0068] The angle of attack α and the roll angle σ are used as the control variables U = [α σ] T ;

[0069] Step 122: Set the initial time and terminal time to t0 and t f , the position, trajectory inclination and trajectory deviation of the interceptor missile at the beginning of the mid-range guidance are determined, and the interception target position can generally be predicted in advance, that is, the position of the interceptor missile at the terminal moment is determined, then the boundary constraint of the trajectory optimization problem of the interceptor missile in the mid-range guidance stage is expressed as:

[0070]

[0071] Among them, x(t0), y(t0), z(t0) represent the coordinate values ​​of the interceptor missile in the launch coordinate system at the initial moment, x0, y0, z0 represent the initial coordinate values ​​of the given interceptor missile in the launch coordinate system, v(t0) represents the velocity of the interceptor missile at the initial moment, v0 represents the given initial velocity of the interceptor missile, θ(t0) represents the ballistic inclination of the interceptor missile at the initial moment, θ0 represents the given initial ballistic inclination of the interceptor missile, ψ(t0) represents the ballistic deviation angle of the interceptor missile at the initial moment, ψ0 represents the given initial ballistic deviation angle of the interceptor missile, x(t f ),y(t f ),z(t f ) represents the coordinate value of the interceptor missile in the launch coordinate system at the terminal moment, x f ,y f ,z f It represents the terminal coordinate value of a given interceptor missile in the launch coordinate system, that is, the predicted target position;

[0072] Step 1, 2, 3: During the flight of the interceptor missile in the mid-guidance phase, the process constraints are mainly: angle of attack constraint, overload constraint, dynamic pressure constraint and heat flux density constraint;

[0073] The maximum angle of attack, maximum overload, maximum dynamic pressure and maximum heat flux are represented by α max 、n max ,q max and Q max ,but:

[0074] -α max ≤α≤α max (7)

[0075]

[0076]

[0077] Q=Cρ 0.5 v 3.05 ≤Q max (10)

[0078] Among them, a y 、a z They represent the overload of the interceptor missile on the y-axis and z-axis of the launch coordinate system respectively;

[0079] q represents the dynamic pressure of the interceptor missile;

[0080] Q represents the heat flux density of the interceptor missile;

[0081] C is the coefficient of heat flux density, which is usually related to the material and size of the aircraft. Here, it is taken as C = 1.1813 × 10 -4 ;

[0082] Assume that there are m interceptors that cooperate to intercept the same target, m>1, that is, m interceptors are required to start from different initial positions and hit the target at the same time after flying along their respective interception trajectories. The cooperative interception time is T, and the superscript i represents the i-th interceptor, i=1,…,m. Then the flight time constraint becomes:

[0083]

[0084] in, Indicates the terminal time of the first interceptor missile, represents the terminal time of the i-th interceptor, represents the terminal moment of the mth interceptor;

[0085] In order to avoid the collision of two interceptors before hitting the target, it is necessary to add a collision avoidance constraint, that is, at a specific time t close to the terminal time d Previously, the distance between each two interceptors should always be greater than the minimum safety distance r:

[0086] Set the collision avoidance constraint. The collision avoidance constraint is:

[0087]

[0088] Among them, (x i ,y i ,z i ) represents the coordinate value of the i-th interceptor missile in the launch coordinate system;

[0089] (x β ,y β ,z β ) represents the coordinate value of the βth interceptor missile in the launch coordinate system, and β represents the βth interceptor missile;

[0090] r represents the minimum safety distance set, r>0;

[0091] t represents time, t0 represents the initial time, t d Indicates the specific time of setting, t0<t d <t f .

[0092] The other steps and parameters are the same as those in the first or second embodiment.

[0093] Specific implementation method 4: This implementation method is different from the specific implementation methods 1 to 3 in that, in step 13, based on step 11 and step 12, the maximum sum of the terminal velocities of multiple interceptors is used as the performance index to establish the original problem P0 of multi-interceptor time coordinated trajectory optimization; the specific process is as follows:

[0094] In order to ensure the high survivability and kill probability of the subsequent terminal guidance phase of the interceptor missile, the maximum sum of the terminal velocities of m interceptor missiles is taken as the performance index (unknown);

[0095] The original problem P0 of multi-interceptor time coordinated trajectory optimization is expressed as:

[0096]

[0097] Among them, v i (t f ) represents the terminal velocity of the i-th interceptor, J represents the performance index, α i represents the angle of attack of the ith interceptor, σ i represents the roll angle of the i-th interceptor missile.

[0098] The other steps and parameters are the same as those in Specific Embodiments 1 to 3.

[0099] Specific implementation method 5: This implementation method is different from any one of specific implementation methods 1 to 4 in that in step 2, the interceptor missile velocity v, ballistic inclination angle θ, ballistic deviation angle ψ, and control quantities α, σ in the three-dimensional dynamic model of the guidance section of the interceptor missile are represented by x, y, z in the launch coordinate system;

[0100] The specific process is:

[0101] According to the first three equations in formula (1), we can get:

[0102]

[0103] Derivation of formula (14) yields:

[0104]

[0105] in, are the second-order derivatives of x, y, and z respectively;

[0106] The control quantity can be obtained by inversely solving the last two equations in equation (1) and equations (2) and (3):

[0107]

[0108] Note that the interceptor missile in the mid-stage guidance phase is a gliding flight, and there is no thrust effect. The fourth equation in the dynamic equation (1) is not used in the process of calculating the control quantity, so it needs to be added as a constraint condition in the solution of the problem.

[0109] Formulas (14), (15), and (16) must satisfy the constraint of formula (17);

[0110]

[0111] Among them, δ is a constant, a very small number.

[0112] The other steps and parameters are the same as those in Specific Embodiments 1 to 4.

[0113] Specific implementation method 6: This implementation method is different from any one of specific implementation methods 1 to 5 in that in step 3, the three-dimensional coordinates x, y, z of the interceptor missile in the launch coordinate system are characterized by using Bezier functions;

[0114] Find the first two and last two control points of the Bezier function according to the boundary conditions;

[0115] The original problem P0 is transformed into a parameter optimization problem P1, where the parameters to be optimized are the remaining control points of the Bezier function;

[0116] The specific process is:

[0117] Step 31: Use Bezier function to characterize the three-dimensional coordinates x, y, z of the interceptor missile in the launch coordinate system; the specific process is:

[0118] First, consider a single interceptor missile and express the three-dimensional coordinates x, y, z of the interceptor missile in the launch coordinate system as Bezier function form:

[0119]

[0120] Among them, x(τ), y(τ), z(τ) are the Bezier functions corresponding to the three-dimensional coordinates x, y, z of the interceptor missile in the launch coordinate system; τ is the time interval;

[0121] n x ,n y ,n z is the order of the Bezier function;

[0122] B x,j ,B y,j ,B z,j are the Bernstein basis functions under the x, y, and z axes of the emission coordinate system respectively;

[0123] B x,j (τ) is the Bernstein basis function in the x direction of the emission coordinate system, B y,j (τ) is the Bernstein basis function in the y direction of the emission coordinate system, B z,j (τ) is the Bernstein basis function in the z direction of the emission coordinate system;

[0124] P x,j ,P y,j ,P z,j They are the control points of the Bezier function under the x, y, and z axes of the launch coordinate system respectively;

[0125] When P x,j ,P y,j ,P z,j After that, the Bezier function is determined, so in the subsequent solution process, P x,j ,P y,j ,P z,j is the variable to be optimized;

[0126] Taking the x direction as an example, for the x direction of the emission coordinate system, the Bernstein basis function is:

[0127]

[0128] Where j is the jth order of the Bezier function;

[0129] Step 32: Find the first two and last two control points of the Bezier function according to the boundary conditions; the specific process is:

[0130] To use the boundary conditions to find the four control points before and after the Bezier function, we first need to transform the flight time t∈[0,T] into the interval τ∈[0,1], where

[0131] The boundary conditions are:

[0132]

[0133] Wherein, the superscripts ' and · denote the derivative of the dimensionless time parameter τ and the flight time t, respectively;

[0134] x(τ=0), y(τ=0), z(τ=0) respectively represent the Bezier function values ​​corresponding to the three-dimensional coordinates x, y, z of the interceptor missile in the launch coordinate system when τ=0;

[0135] x(τ=1), y(τ=1), z(τ=1) respectively represent the Bezier function values ​​corresponding to the three-dimensional coordinates x, y, z of the interceptor missile in the launch coordinate system when τ=1;

[0136] x′(τ=0), y′(τ=0), z′(τ=0) represent the first-order derivative values ​​of the Bezier function corresponding to the three-dimensional coordinates x, y, z of the interceptor missile in the launch coordinate system when τ=0 respectively;

[0137] x′(τ=1), y′(τ=1), z′(τ=1) represent the first-order derivative values ​​of the Bezier function corresponding to the three-dimensional coordinates x, y, and z of the interceptor missile when τ=1 respectively;

[0138] represents the first-order derivative of x0, represents the first-order derivative of y0, represents the first-order derivative of z0;

[0139] Represents x f The first derivative of Represents y f The first derivative of Indicates z f The first derivative of ;

[0140] Taking the x direction as an example, the derivative of equation (18) is obtained in the x direction, and the derivative of x(τ) is:

[0141]

[0142] Among them, B′ x,j (τ) is derived from formula (19) to obtain:

[0143]

[0144] Among them, B′ x,j (τ) represents B x,j The first derivative of (τ);

[0145] Substituting τ = 0 and τ = 1 into equations (19) and (22), the parameter form of the boundary conditions is obtained as follows:

[0146]

[0147] Formula (23) combined with formula (20) and (18) yields:

[0148]

[0149] Among them, P x,0 Represents the first control point of the Bezier function in the x direction, P x,1Represents the second control point of the Bezier function in the x direction;

[0150] represents the penultimate control point of the Bezier function in the x direction, Indicates the last control point of the Bezier function in the x direction;

[0151] The first two and last two control points of the Bezier function are determined by formula (24):

[0152]

[0153] Step 3: Based on the first two and last two control points of the Bezier function (Equation (25)), the remaining (n x -3) control points are the variables to be optimized (a total of (n x +1) control point); the same applies to y and z directions;

[0154] For the single interceptor trajectory optimization problem, there are (n x +n y +n z -9) variables to be optimized, expressed as

[0155] in,

[0156] Indicates the x direction to be optimized (n x -3) control points, P x,2 represents the third control point of the Bezier function in the x direction, Represents the third to last control point of the Bezier function in the x direction;

[0157] Indicates the y direction to be optimized (n x -3) control points, P y,2 Represents the third control point of the Bezier function in the y direction, Represents the third to last control point of the Bezier function in the y direction;

[0158] Indicates the z direction to be optimized (n x -3) control points, P z,2 represents the third control point of the Bezier function in the z direction, Represents the third to last control point of the Bezier function in the z direction;

[0159] The superscript T means to find the transpose;

[0160] The present invention does not fix the terminal velocity, terminal trajectory inclination and terminal trajectory deviation of each interceptor missile, so the terminal velocity v of each interceptor missile is set to f 、Terminal trajectory inclination angle θ f and the terminal trajectory angle ψ f Also used as the variable to be optimized, expressed as [v f θ f ψ f ] T , a total of (n x +n y +n z -6) variables to be optimized;

[0161] If there are m interceptors for time-coordinated trajectory optimization, and m>1, then the variable to be optimized becomes Total (n x +n y +n z -6)×m variables to be optimized;

[0162] in, represents the x-direction of the i-th interceptor to be optimized (n x -3) array of control points, Indicates the y direction of the i-th interceptor to be optimized (n x -3) array of control points, represents the z direction of the i-th interceptor to be optimized (n x -3) an array of control points;

[0163] represents the terminal velocity of the ith interceptor, represents the terminal ballistic inclination angle of the i-th interceptor missile, represents the terminal ballistic deviation angle of the i-th interceptor missile;

[0164] For each interceptor missile, according to equations (14), (15), (16), use Bezier functions to describe x, y, z to represent other state variables v, θ and ψ as well as control variables α and σ, so that they satisfy constraints (7), (8), (9), (10), (17). In addition, constraints (11) and (12) need to be added.

[0165] Step 34: Based on step 32 and step 33, the multi-interceptor time coordinated trajectory optimization problem P0 is transformed into a parameter optimization problem P1.

[0166] The other steps and parameters are the same as those in Specific Implementation Methods 1 to 5.

[0167] Specific implementation method 7: This implementation method is different from the specific implementation methods 1 to 6 in that, in step 34, based on step 32 and step 33, the multi-interceptor time coordinated trajectory optimization problem P0 is transformed into a parameter optimization problem P1; the specific process is:

[0168]

[0169] Among them, x i (t0),y i (t0),z i (t0) respectively represent the coordinate values ​​of the i-th interceptor missile at the initial moment in the launch coordinate system;

[0170] They respectively represent the initial coordinate values ​​of the i-th interceptor missile given in the launch coordinate system;

[0171] v i (t0) represents the velocity of the i-th interceptor at the initial moment, θ i (t0) represents the ballistic inclination angle of the i-th interceptor at the initial moment, ψ i (t0) represents the ballistic deviation angle of the i-th interceptor missile at the initial moment;

[0172] represents the initial velocity given by the i-th interceptor, represents the initial ballistic inclination angle given by the i-th interceptor missile, represents the initial trajectory deviation angle given by the i-th interceptor missile;

[0173] x i (t f ),y i (t f ),z i (t f ) represent the coordinate values ​​of the i interceptor missiles at the terminal moment in the launch coordinate system;

[0174] They respectively represent the terminal coordinate values ​​of the i-th interceptor missile given in the launch coordinate system;

[0175] α i represents the angle of attack of the ith interceptor missile, represents the overload of the ith interceptor missile on the y-axis of the launch coordinate system, represents the overload of the i-th interceptor missile on the z-axis of the launch coordinate system;

[0176] L i represents the aerodynamic lift of the ith interceptor, m i represents the mass of the ith interceptor, g i represents the gravitational acceleration of the i-th interceptor;

[0177] ρ i represents the atmospheric density of the ith interceptor, v i represents the speed of the i-th interceptor, q i represents the dynamic pressure of the i-th interceptor;

[0178] Q i represents the heat flux density of the i-th interceptor, represents the first-order derivative of the velocity of the i-th interceptor missile, D i represents the aerodynamic drag of the i-th interceptor, θ i represents the ballistic inclination angle of the i-th interceptor missile, and δ is a given constant.

[0179] The other steps and parameters are the same as those in Specific Embodiments 1 to 6.

[0180] Specific implementation eight: This implementation differs from any one of specific implementations one to seven in that in step four, the Bezier functions of the first two and last two control points in the parameter optimization problem P1 are discretized and the parameters are initialized; the specific process is:

[0181] Step 41: Discretize the Bezier functions of the first two and last two control points in the parameter optimization problem P1; the specific process is:

[0182] Problem P1 is still an infinite-dimensional optimization problem, which needs to be solved through discretization and parameter initialization in step 3. The specific operation of discretization is: by selecting K discrete points on the interceptor trajectory, all control quantities and state quantities of the interceptor are expressed in discrete form, and the constraints in P1 are established at these K discrete points, so that it can be converted into a finite-dimensional optimization problem;

[0183] The present invention adopts the root discretization interval τ of the K-order Legendre polynomial:

[0184] 0=τ1<τ2<…<τ K-1 <τ K =1 (27)

[0185] Among them, in order from small to large, τ1 represents the first root of the Legendre polynomial, τ2 represents the second root of the Legendre polynomial, τ K represents the Kth root of the Legendre polynomial, i.e. the last root, τ K-1 represents the K-1th root of the Legendre polynomial; K ≥ 1;

[0186] For a single interceptor missile, in the x direction, use represents all the control points of the Bezier function describing x(τ), that is, Then the interceptor missile coordinate x in the launch coordinate system is expressed at K discrete points as:

[0187]

[0188] in, is the weight matrix of all control points, for The first derivative of for The second derivative of

[0189] [x] K×1 is the value of the interceptor missile coordinate x at K discrete points in the launch coordinate system, [x′] K×1 for [x] K×1 The first derivative of, [x″] K×1 for [x] K×1 The second derivative of

[0190] To k , 1≤k≤K is substituted into formula (19) to obtain:

[0191]

[0192] Among them, B x,0 (τ1) represents the value of the Bernstein basis function in the x direction of the emission coordinate system when j = 0, τ = τ1, The Bernstein basis function in the x direction of the emission coordinate system is j = n x , the value when τ=τ1, B x,0 (τ K ) represents the Bernstein basis function in the x direction of the emission coordinate system at j = 0, τ = τ K The value of The Bernstein basis function in the x direction of the emission coordinate system is j = n x ,τ=τ K The value at time;

[0193] B x ' ,0 (τ1) represents the value of the first derivative of the Bernstein basis function in the x direction of the emission coordinate system when j = 0, τ = τ1, The first derivative of the Bernstein basis function in the x direction of the emission coordinate system is j = n x , the value when τ=τ1, B x ' ,0 (τ K ) represents the first-order derivative of the Bernstein basis function in the x-direction of the emission coordinate system at j = 0, τ = τ K The value of The first derivative of the Bernstein basis function in the x direction of the emission coordinate system is j = n x ,τ=τ K The value at time;

[0194] B″ x,0 (τ1) represents the value of the second-order derivative of the Bernstein basis function in the x-direction of the emission coordinate system when j = 0, τ = τ1, The second derivative of the Bernstein basis function in the x direction of the emission coordinate system is j = n x , the value when τ=τ1, B″ x,0 (τ K ) represents the second-order derivative of the Bernstein basis function in the x-direction of the emission coordinate system at j = 0, τ = τ K The value of The second derivative of the Bernstein basis function in the x direction of the emission coordinate system is j = n x ,τ=τ K The value at time.

[0195] When the Bezier function order n x And when the number of discrete points K is determined, and The first and second order derivatives of are constant matrices; the same is true for the y and z directions;

[0196] Substituting equation (28) into equations (14) and (15) yields the discretized form of all state variables and the first-order derivatives of the state variables [v] m×1 , [ψ] m×1 , [θ] m×1 ,

[0197] Substituting equation (28) into equation (16) gives the discretized form [α] that represents all control quantities: m×1 , [σ] m×1 ;

[0198] will [v] m×1 , [ψ] m×1 , [θ] m×1 , [α] m×1 , [σ] m×1 Substituting into P1 will transform P1 into a finite-dimensional discretized form;

[0199] Step 42: Initialize the parameters of the Bezier functions of the first two and last two known control points in P1.

[0200] The other steps and parameters are the same as those in Specific Embodiments 1 to 7.

[0201] Specific implementation method 9: This implementation method is different from any one of specific implementation methods 1 to 8 in that in step 42, the parameters of the Bezier function of the first two and last two known control points in P1 are initialized;

[0202] The specific process is:

[0203] Set the variable to be optimized [v f θ f ψ f ] T The initial value of the variable to be optimized [v f θ f ψ f ] T The initial value of is given based on experience;

[0204] The calculation process of the initial value of is:

[0205] Generally speaking, the Bezier curve order needs to be greater than 3 to more accurately characterize the trajectory of the interceptor missile;

[0206] Here it is assumed that the three-dimensional coordinates of the interceptor missile are completely determined by the third-order Bezier curve, that is:

[0207]

[0208] In the formula, x APP (τ),y APP (τ),z APP (τ) respectively represent the initial values ​​of the three-dimensional coordinates of the interceptor missile in the launch coordinate system represented by the third-order Bezier curve;

[0209] P x,0 Represents the first control point of the Bezier function in the x direction in the launch coordinate system, P y,0 Indicates the first control point of the Bezier function in the y direction in the launch coordinate system, P z,0 Represents the first control point of the Bezier function in the z direction in the launch coordinate system;

[0210] P x,1 Represents the second control point of the Bezier function in the x direction in the launch coordinate system, P y,1 Represents the second control point of the Bezier function in the y direction in the launch coordinate system, P z,1 Represents the second control point of the Bezier function in the z direction in the launch coordinate system;

[0211] P x,2 represents the third (second to last) control point of the Bezier function in the x direction in the launch coordinate system, P y,2represents the third (second to last) control point of the Bezier function in the y direction in the launch coordinate system, P z,2 Represents the third (second to last) control point of the Bezier function in the z direction in the launch coordinate system;

[0212] P x,3 represents the fourth (last) control point of the Bezier function in the x direction in the launch coordinate system, P y,3 Represents the fourth (last) control point of the Bezier function in the y direction in the launch coordinate system, P z,3 Represents the fourth (last) control point of the Bezier function in the z direction in the launch coordinate system;

[0213] Similar to the solution method of formula (25), the values ​​of all control points in formula (30) are obtained:

[0214]

[0215] In actual situations, the order of the Bezier curve is generally greater than 3. Substituting equation (30) into equation (28) yields the initial values ​​of the Bezier curve parameters:

[0216]

[0217] In the formula, the superscript - 1 represents matrix inversion;

[0218] Indicates the x direction to be optimized in the launch coordinate system (n x -3) The initial value of the array consisting of control points, Indicates the x direction to be optimized in the launch coordinate system (n x -3) weight matrix of control points, [x APP ] K×1 It represents the initial value of the x-axis coordinate of the interceptor missile represented by the third-order Bezier curve in the launch coordinate system;

[0219] Indicates the y direction to be optimized in the launch coordinate system (n y -3) The initial value of the array consisting of control points, Indicates the y direction to be optimized in the launch coordinate system (n y -3) weight matrix of control points, [y APP ] K×1 It represents the initial value of the y-axis coordinate of the interceptor missile represented by the third-order Bezier curve in the launch coordinate system;

[0220] Indicates the z direction to be optimized in the launch coordinate system (n z-3) The initial value of the array consisting of control points, Indicates the z direction to be optimized in the launch coordinate system (n z -3) weight matrix of control points, [z APP ] K×1 It represents the initial value of the z-axis coordinate of the interceptor missile represented by the third-order Bezier curve in the launch coordinate system.

[0221] The other steps and parameters are the same as those in Specific Embodiments 1 to 8.

[0222] Specific implementation method ten: This implementation method is different from any one of specific implementation methods one to nine in that in step five, a nonlinear programming solver is called to find all control points of the Bezier function, and the values ​​of all state quantities and control quantities are found according to all control points of the Bezier function;

[0223] The specific process is:

[0224] Step 51: Call a nonlinear programming solver, such as the sqp algorithm, to solve P1 and obtain all control points of the Bezier function;

[0225] Step 52: Based on all the control points of the Bezier function, the Bezier function form (18) represented by the three-dimensional coordinates x, y, z of the interceptor missile in the launch coordinate system, as well as equations (14) and (16) are solved to obtain the values ​​of all state quantities and control quantities.

[0226] Complete multiple interceptor missiles to strike the same target, coordinating the terminal speed with the optimal trajectory in time.

[0227] The other steps and parameters are the same as those in Specific Embodiments 1 to 9.

[0228] Example

[0229] The initial positions, velocities, ballistic inclination angles and ballistic deviation angles of the interceptor missiles 1, 2 and 3 used in this embodiment are shown in Table 1. The interceptor missiles 1, 2 and 3 are represented by superscript 1 , 2 , 3 The target position, angle of attack, overload, dynamic pressure, heat flux density constraint and coordinated interception time in the embodiment are shown in Table 1.

[0230] Table 1 Parameter values

[0231]

[0232] Consider the case where interceptors 1, 2 and 3 cooperate to intercept the same target. The time cooperative trajectory optimization method based on Bezier curve proposed by the present invention is adopted and solved by fmincon function in MATLAB platform. It takes 2.62 seconds to obtain the three-dimensional trajectory of cooperative interception. Figure 1 As shown:

[0233] The interceptor missile speed curve is as follows Figure 2 As shown:

[0234] The interceptor missile attack angle curve is as follows Figure 3 As shown:

[0235] The interceptor missile overload curve is as follows Figure 4 As shown:

[0236] The interception dynamic pressure curve is as follows Figure 5 As shown:

[0237] The heat flux curve of the interceptor missile is as follows: Figure 6 As shown:

[0238] Depend on Figure 3 , 4 As can be seen from Figures 5 and 6, in this embodiment, the angle of attack, overload, dynamic pressure and heat flux density solved by the trajectory optimization method based on Bezier curve proposed in the present invention all meet the constraints.

[0239] The proposed method and the pseudo-spectral method were used for comparative simulation experiments, in which the proposed method provided initial values ​​for the pseudo-spectral method, and the results of three interceptors coordinating to intercept the same target were shown in Table 2. The Bezier curve-based method only took 2.62 seconds, which only took up 45.02% of the time taken to start the pseudo-spectral method with the Bezier initial value, but was able to achieve 91.46% of the effect of the pseudo-spectral method, which illustrates the high efficiency and feasibility of the method.

[0240] Table 2 Comparison of Bezier and pseudo-spectral results (three interceptors)

[0241]

[0242] The comparison of the three-dimensional trajectory of cooperative interception obtained by the two methods is as follows: Figure 8 As shown:

[0243] In addition, the performance indicators of the three interceptor missiles in the two cases of coordinated interception and separate interception are calculated respectively. As shown in Table 3, the sum of the terminal velocities of the three interceptor missiles in the case of coordinated interception is 6080.6 m / s, and the sum of the terminal velocities of the three interceptor missiles in the case of separate interception is 5787.2 m / s. The use of coordinated interception improves the performance indicators by 5.07%, indicating that coordinated interception can effectively increase the terminal velocities of the three interceptor missiles compared with separate interception.

[0244] Table 3 Comparison of performance indicators in cooperative interception and separate interception

[0245]

[0246] The present invention may also have many other embodiments. Without departing from the spirit and essence of the present invention, those skilled in the art may make various corresponding changes and modifications based on the present invention, but these corresponding changes and modifications should all fall within the scope of protection of the claims attached to the present invention.

Claims

1. A fast trajectory optimization method for many-to-one coordinated interception, characterized by: The specific process of the method is: Step 1: Establish a three-dimensional dynamic model of the interceptor missile's mid-range guidance section and set boundary constraints and process constraints of the interceptor missile; The maximum sum of terminal velocities of multiple interceptors is taken as the performance index, and the original problem P0 of multi-interceptor time coordinated trajectory optimization is established; Step 2: Use x, y, z in the launch coordinate system to represent the interceptor missile velocity v, trajectory inclination angle θ, trajectory deviation angle ψ, and control quantities α, σ in the three-dimensional dynamic model of the interceptor missile's guidance section; Step 3: Use Bezier function to characterize the three-dimensional coordinates x, y, and z of the interceptor missile in the launch coordinate system; Find the first two and last two control points of the Bezier function according to the boundary conditions; Transform the original problem P0 into a parameter optimization problem P1; Step 4: Discretize and initialize the parameters of the Bezier functions of the first two and last two control points in the parameter optimization problem P1; Step 5: Call the nonlinear programming solver to find all the control points of the Bezier function, and find the values ​​of all state quantities and control quantities based on all the control points of the Bezier function.

2. The method for fast trajectory optimization of many-to-one coordinated interception according to claim 1, characterized in that: In the step 1, a three-dimensional dynamic model of the guidance section of the interceptor missile is established, and boundary constraints and process constraints of the interceptor missile are set; The maximum sum of terminal velocities of multiple interceptors is taken as the performance index, and the original problem P0 of multi-interceptor time coordinated trajectory optimization is established; The specific process is: Step 1: Establish a three-dimensional dynamic model of the guidance section of the interceptor missile; the specific process is: Treat the interception target as a fixed point; Ignoring the influence of the earth's curvature and rotation, the interceptor missile adopts unpowered gliding in the mid-range guidance stage. The three-dimensional dynamic model of the interceptor missile in the mid-range guidance stage in the launch coordinate system is described as: In the formula, x, y, z are the coordinate values ​​of the interceptor missile in the launch coordinate system; v is the interceptor missile speed; θ is the trajectory inclination angle; ψ is the trajectory deviation angle; g is the acceleration of gravity, σ is the roll angle, and m is the mass of the interceptor missile; are the first-order derivatives of x, y, and z respectively; is the first-order derivative of v; is the first-order derivative of θ; is the first-order derivative of ψ; L and D are aerodynamic lift and aerodynamic drag, respectively; The expressions of L and D are: Among them, ρ is the atmospheric density, S is the area of ​​the interceptor missile, and C L ,C D are the aerodynamic lift coefficient and drag coefficient respectively; C L ,C D The expression is: Where α is the angle of attack; The expression of atmospheric density ρ is: p=ρ0e -y / href (4) Where, ρ0=1.226kg / m 3 is the atmospheric density at sea level, h ref =7254.24m is the reference altitude, y is the coordinate value of the interceptor missile on the y-axis of the launch coordinate system; The expression of gravitational acceleration g is: g=g0(R e / r) 2 (5) Where g0 is the gravitational acceleration at sea level; R e is the radius of the Earth; r represents the distance from the interceptor to the center of the earth, r = R e +y; Step 1 and 2: Set the boundary constraints and process constraints of the interceptor missile; Step 13: Based on steps 11 and 12, the maximum sum of the terminal velocities of multiple interceptors is used as the performance indicator, and the original problem P0 of multi-interceptor time coordinated trajectory optimization is established.

3. The method for fast trajectory optimization of many-to-one coordinated interception according to claim 2, characterized in that: In the steps 1 and 2, the boundary constraints and process constraints of the interceptor missile are set; the specific process is: Step 121: Take the coordinates x, y, z of the interceptor missile in the launch coordinate system, the interceptor missile speed v, the ballistic inclination angle θ, and the ballistic deviation angle ψ as the state variable X = [xyzv θ ψ] T , the superscript T means to find the transpose; The angle of attack α and the roll angle σ are used as the control variables U = [α σ] T ; Step 122: Set the initial time and terminal time to t0 and t f , then the boundary constraint of the trajectory optimization problem of the interceptor missile in the mid-stage guidance phase is expressed as: Among them, x(t0), y(t0), z(t0) represent the coordinate values ​​of the interceptor missile in the launch coordinate system at the initial moment, x0, y0, z0 represent the initial coordinate values ​​of the given interceptor missile in the launch coordinate system, v(t0) represents the velocity of the interceptor missile at the initial moment, v0 represents the given initial velocity of the interceptor missile, θ(t0) represents the ballistic inclination of the interceptor missile at the initial moment, θ0 represents the given initial ballistic inclination of the interceptor missile, ψ(t0) represents the ballistic deviation angle of the interceptor missile at the initial moment, ψ0 represents the given initial ballistic deviation angle of the interceptor missile, x(t f ),y(t f ),z(t f ) represents the coordinate value of the interceptor missile in the launch coordinate system at the terminal moment, x f ,y f ,z f It represents the terminal coordinate value of a given interceptor missile in the launch coordinate system, that is, the predicted target position; Step 1, 2, 3: During the flight of the interceptor missile in the mid-guidance phase, the process constraints are: angle of attack constraint, overload constraint, dynamic pressure constraint and heat flux density constraint; The maximum angle of attack, maximum overload, maximum dynamic pressure and maximum heat flux are expressed as α max 、n max ,q max and Q max ,but: -α max ≤α≤α max (7) Q=Cρ 0.5 v 3.05 ≤Q max (10) Among them, a y 、a z They represent the overload of the interceptor missile on the y-axis and z-axis of the launch coordinate system respectively; q represents the dynamic pressure of the interceptor missile; Q represents the heat flux density of the interceptor missile; C is the coefficient of heat flux; Assume that there are m interceptors cooperating to intercept the same target, m>1, the cooperative interception time is T, and the superscript i represents the i-th interceptor, i=1,…,m, then the flight time constraint becomes: in, Indicates the terminal time of the first interceptor missile, represents the terminal time of the i-th interceptor, represents the terminal moment of the mth interceptor; Set the collision avoidance constraint. The collision avoidance constraint is: Among them, (x i ,y i ,z i ) represents the coordinate value of the i-th interceptor missile in the launch coordinate system; (x β ,y β ,z β ) represents the coordinate value of the βth interceptor missile in the launch coordinate system, and β represents the βth interceptor missile; r represents the minimum safety distance set, r>0; t represents time, t0 represents the initial time, t d Indicates the specific time of setting, t0<t d <t f .

4. The method for fast trajectory optimization of many-to-one coordinated interception according to claim 3 is characterized in that: In step 13, based on step 11 and step 12, the maximum sum of the terminal velocities of multiple interceptors is used as the performance index to establish the original problem P0 of multi-interceptor time coordinated trajectory optimization; the specific process is as follows: The maximum sum of the terminal velocities of m interceptor missiles is taken as the performance index; The original problem P0 of multi-interceptor time coordinated trajectory optimization is expressed as: Among them, v i (t f ) represents the terminal velocity of the i-th interceptor, J represents the performance index, α i represents the angle of attack of the ith interceptor, σ i represents the roll angle of the i-th interceptor missile.

5. The method for fast trajectory optimization of many-to-one coordinated interception according to claim 4 is characterized in that: In the step 2, the interceptor missile velocity v, trajectory inclination angle θ, trajectory deviation angle ψ, and control quantities α and σ in the three-dimensional dynamic model of the guidance section of the interceptor missile are represented by x, y, and z in the launch coordinate system; The specific process is: According to the first three equations in formula (1), we can get: Derivation of formula (14) yields: in, are the second-order derivatives of x, y, and z respectively; The control quantity can be obtained by inversely solving the last two equations in equation (1) and equations (2) and (3): Formulas (14), (15), and (16) must satisfy the constraint of formula (17); Here, δ is a constant.

6. The method for fast trajectory optimization of many-to-one coordinated interception according to claim 5, characterized in that: In the step 3, the three-dimensional coordinates x, y, and z of the interceptor missile in the launch coordinate system are characterized by using the Bezier function; Find the first two and last two control points of the Bezier function according to the boundary conditions; Transform the original problem P0 into a parameter optimization problem P1; The specific process is: Step 31: Use Bezier function to characterize the three-dimensional coordinates x, y, z of the interceptor missile in the launch coordinate system; the specific process is: First, consider a single interceptor missile and express the three-dimensional coordinates x, y, z of the interceptor missile in the launch coordinate system as Bezier function form: Among them, x(τ), y(τ), z(τ) are the Bezier functions corresponding to the three-dimensional coordinates x, y, z of the interceptor missile in the launch coordinate system; τ is the time interval; n x ,n y ,n z is the order of the Bezier function; B x,j ,B y,j ,B z,j are the Bernstein basis functions under the x, y, and z axes of the emission coordinate system respectively; B x,j (τ) is the Bernstein basis function in the x direction of the emission coordinate system, B y,j (τ) is the Bernstein basis function in the y direction of the emission coordinate system, B z,j (τ) is the Bernstein basis function in the z direction of the emission coordinate system; P x,j ,P y,j ,P z,j They are the control points of the Bezier function under the x, y, and z axes of the launch coordinate system respectively; For the x direction of the emission coordinate system, the Bernstein basis function is: Where j is the jth order of the Bezier function; Step 32: Find the first two and last two control points of the Bezier function according to the boundary conditions; the specific process is: Convert the flight time t∈[0,T] to the interval τ∈[0,1], where The boundary conditions are: Wherein, the superscripts ' and · denote the derivative of the dimensionless time parameter τ and the flight time t, respectively; x(τ=0), y(τ=0), z(τ=0) respectively represent the Bezier function values ​​corresponding to the three-dimensional coordinates x, y, z of the interceptor missile in the launch coordinate system when τ=0; x(τ=1), y(τ=1), z(τ=1) respectively represent the Bezier function values ​​corresponding to the three-dimensional coordinates x, y, z of the interceptor missile in the launch coordinate system when τ=1; x′(τ=0), y′(τ=0), z′(τ=0) represent the first-order derivative values ​​of the Bezier function corresponding to the three-dimensional coordinates x, y, z of the interceptor missile in the launch coordinate system when τ=0 respectively; x′(τ=1), y′(τ=1), z′(τ=1) represent the first-order derivative values ​​of the Bezier function corresponding to the three-dimensional coordinates x, y, and z of the interceptor missile when τ=1 respectively; represents the first-order derivative of x0, represents the first-order derivative of y0, represents the first-order derivative of z0; Represents x f The first derivative of Represents y f The first derivative of Indicates z f The first derivative of ; In the x direction, the derivative of equation (18) is: Among them, B′ x,j (τ) is derived from formula (19) to obtain: Among them, B′ x,j (τ) represents B x,j The first derivative of (τ); Substituting τ = 0 and τ = 1 into equations (19) and (22), the parameter form of the boundary conditions is obtained as follows: Formula (23) combined with formula (20) and (18) yields: Among them, P x,0 Represents the first control point of the Bezier function in the x direction, P x,1 Represents the second control point of the Bezier function in the x direction; represents the penultimate control point of the Bezier function in the x direction, Indicates the last control point of the Bezier function in the x direction; The first two and last two control points of the Bezier function are determined by formula (24): Step 3. Based on the first two and last two control points of the Bezier function, the remaining (n x -3) control points are variables to be optimized; For the single interceptor trajectory optimization problem, there are (n x +n y +n z -9) variables to be optimized, expressed as in, Indicates the x direction to be optimized (n x -3) control points, P x,2 represents the third control point of the Bezier function in the x direction, Represents the third to last control point of the Bezier function in the x direction; Indicates the y direction to be optimized (n x -3) control points, P y,2 Represents the third control point of the Bezier function in the y direction, Represents the third to last control point of the Bezier function in the y direction; Indicates the z direction to be optimized (n x -3) control points, P z,2 represents the third control point of the Bezier function in the z direction, Represents the third to last control point of the directional Bezier function; The superscript T means to find the transpose; The terminal velocity v of each interceptor missile is f 、Terminal trajectory inclination angle θ f and the terminal trajectory angle ψ f Also used as the variable to be optimized, expressed as [v f θ f ψ f ] T , a total of ( n x +n y +n z -6 ) variables to be optimized; If there are m interceptors for time-coordinated trajectory optimization, and m>1, then the variable to be optimized becomes Total (n x +n y +n z -6)×m variables to be optimized; in, represents the x-direction of the i-th interceptor to be optimized (n x -3) array of control points, Indicates the y direction of the i-th interceptor to be optimized (n x -3) array of control points, represents the z direction of the i-th interceptor to be optimized (n x -3) an array of control points; represents the terminal velocity of the ith interceptor, represents the terminal ballistic inclination angle of the i-th interceptor missile, represents the terminal ballistic deviation angle of the i-th interceptor missile; Step 34: Based on step 32 and step 33, the multi-interceptor time coordinated trajectory optimization problem P0 is transformed into a parameter optimization problem P1.

7. The method for fast trajectory optimization of many-to-one coordinated interception according to claim 6, characterized in that: In step 34, based on step 32 and step 33, the multi-interceptor time coordinated trajectory optimization problem P0 is transformed into a parameter optimization problem P1; the specific process is: Among them, x i (t0),y i (t0),z i (t0) respectively represent the coordinate values ​​of the i-th interceptor missile at the initial moment in the launch coordinate system; They respectively represent the initial coordinate values ​​of the i-th interceptor missile given in the launch coordinate system; v i (t0) represents the velocity of the i-th interceptor at the initial moment, θ i (t0) represents the ballistic inclination angle of the i-th interceptor at the initial moment, ψ i (t0) represents the ballistic deviation angle of the i-th interceptor missile at the initial moment; represents the initial velocity given by the i-th interceptor, represents the initial ballistic inclination angle given by the i-th interceptor missile, represents the initial trajectory deviation angle given by the i-th interceptor missile; x i (t f ),y i (t f ),z i (t f ) represent the coordinate values ​​of the i interceptor missiles at the terminal moment in the launch coordinate system; They respectively represent the terminal coordinate values ​​of the i-th interceptor missile given in the launch coordinate system; α i represents the angle of attack of the ith interceptor missile, represents the overload of the ith interceptor missile on the y-axis of the launch coordinate system, represents the overload of the i-th interceptor missile on the z-axis of the launch coordinate system; L i represents the aerodynamic lift of the ith interceptor, m i represents the mass of the ith interceptor, g i represents the gravitational acceleration of the i-th interceptor; ρ i represents the atmospheric density of the ith interceptor, v i represents the speed of the i-th interceptor, q i represents the dynamic pressure of the i-th interceptor; Q i represents the heat flux density of the i-th interceptor, represents the first-order derivative of the velocity of the i-th interceptor missile, D i represents the aerodynamic drag of the i-th interceptor, θ i represents the ballistic inclination angle of the i-th interceptor missile, and δ is a given constant.

8. The method for fast trajectory optimization of many-to-one coordinated interception according to claim 7, characterized in that: In the step 4, the Bezier functions of the first two and last two control points in the parameter optimization problem P1 are discretized and the parameters are initialized; The specific process is: Step 41: Discretize the Bezier functions of the first two and last two known control points in the parameter optimization problem P1; The specific process is: The interval τ is discretized using the roots of the K-order Legendre polynomial: 0=τ1<τ2<…<τ K-1 <t K =1(27) Among them, in order from small to large, τ1 represents the first root of the Legendre polynomial, τ2 represents the second root of the Legendre polynomial, τ K represents the Kth root, or the last root, of the Legendre polynomial, τ K-1 represents the K-1th root of the Legendre polynomial; K ≥ 1; For a single interceptor missile, in the x direction, use represents all the control points of the Bezier function describing x(τ), that is, Then the interceptor missile coordinate x in the launch coordinate system is expressed at K discrete points as: in, is the weight matrix of all control points, for The first derivative of for The second derivative of [x] K×1 is the value of the interceptor missile coordinate x at K discrete points in the launch coordinate system, [x′] K×1 for [x] K×1 The first derivative of, [x″] K×1 for [x] K×1 The second derivative of To k , 1≤k≤K is substituted into formula (19) to obtain: Among them, B x,0 (τ1) represents the value of the Bernstein basis function in the x direction of the emission coordinate system when j = 0, τ = τ1, The Bernstein basis function in the x direction of the emission coordinate system is j = n x , the value when τ=τ1, B x,0 (τ K ) represents the Bernstein basis function in the x direction of the emission coordinate system at j = 0, τ = τ K The value of The Bernstein basis function in the x direction of the emission coordinate system is j = n x ,τ=τ K The value at time; B′ x,0 (τ1) represents the value of the first derivative of the Bernstein basis function in the x direction of the emission coordinate system when j = 0, τ = τ1, The first derivative of the Bernstein basis function in the x direction of the emission coordinate system is j = n x , the value when τ=τ1, B′ x,0 (τ K ) represents the first-order derivative of the Bernstein basis function in the x-direction of the emission coordinate system at j = 0, τ = τ K The value of The first derivative of the Bernstein basis function in the x direction of the emission coordinate system is j = n x ,τ=τ K The value at time; B x ″ ,0 (τ1) represents the value of the second-order derivative of the Bernstein basis function in the x-direction of the emission coordinate system when j = 0, τ = τ1, The second derivative of the Bernstein basis function in the x direction of the emission coordinate system is j = n x , the value when τ=τ1, B x ″ ,0 (τ K ) represents the second-order derivative of the Bernstein basis function in the x-direction of the emission coordinate system at j = 0, τ = τ K The value of The second derivative of the Bernstein basis function in the x direction of the emission coordinate system is j = n x ,τ=τ K The value at time. Substituting equation (28) into equations (14) and (15) yields the discretized form of all state variables and the first-order derivatives of the state variables [v] m×1 , [ψ] m×1 , [θ] m×1 , Substituting equation (28) into equation (16) gives the discretized form [α] that represents all control quantities: m×1 , [σ] m×1 ; will [v] m×1 , [ψ] m×1 , [θ] m×1 , [α] m×1 , [σ] m×1 Substituting into P1 will transform P1 into a finite-dimensional discretized form; Step 42: Initialize the parameters of the Bezier functions of the first two and last two known control points in P1.

9. The method for fast trajectory optimization of many-to-one coordinated interception according to claim 8, characterized in that: In the step 42, the parameters of the Bezier functions of the first two and last two control points in P1 are initialized; The specific process is: Set the variable to be optimized [v f θ f ψ f ] T The initial value of The calculation process of the initial value of is: Assume that the three-dimensional coordinates of the interceptor missile are completely determined by the third-order Bezier curve, that is: In the formula, x APP (τ),y APP (τ),z APP (τ) respectively represent the initial values ​​of the three-dimensional coordinates of the interceptor missile in the launch coordinate system represented by the third-order Bezier curve; P x,0 Indicates the first control point of the Bezier function in the x direction in the launch coordinate system, P y,0 Represents the first control point of the Bezier function in the y direction in the launch coordinate system, P z,0 Represents the first control point of the Bezier function in the z direction in the launch coordinate system; P x,1 Represents the second control point of the Bezier function in the x direction in the launch coordinate system, P y,1 Represents the second control point of the Bezier function in the y direction in the launch coordinate system, P z,1 Represents the second control point of the Bezier function in the z direction in the launch coordinate system; P x,2 Represents the third control point of the Bezier function in the x direction in the launch coordinate system, P y,2 Represents the third control point of the Bezier function in the y direction in the launch coordinate system, P z,2 Represents the third control point of the Bezier function in the z direction in the launch coordinate system; P x,3 Represents the fourth control point of the Bezier function in the x direction in the launch coordinate system, P y,3 Represents the fourth control point of the Bezier function in the y direction in the launch coordinate system, P z,3 Represents the fourth control point of the Bezier function in the z direction in the launch coordinate system; Get the values ​​of all control points in formula (30): Substituting formula (30) into formula (28) yields the initial values ​​of the Bezier curve parameters: In the formula, the superscript -1 represents matrix inversion; Indicates the x direction to be optimized in the launch coordinate system (n x -3) The initial value of the array consisting of control points, Indicates the x direction to be optimized in the launch coordinate system (n x -3) weight matrix of control points, [x APP ] K×1 It represents the initial value of the x-axis coordinate of the interceptor missile represented by the third-order Bezier curve in the launch coordinate system; Indicates the y direction to be optimized in the launch coordinate system (n y -3) The initial value of the array consisting of control points, Indicates the y direction to be optimized in the launch coordinate system (n y -3) weight matrix of control points, [y APP ] K×1 It represents the initial value of the y-axis coordinate of the interceptor missile represented by the third-order Bezier curve in the launch coordinate system; Indicates the z direction to be optimized in the launch coordinate system (n z -3) The initial value of the array consisting of control points, Indicates the z direction to be optimized in the launch coordinate system (n z -3) weight matrix of control points, [z APP ] K×1 It represents the initial value of the z-axis coordinate of the interceptor missile represented by the third-order Bezier curve in the launch coordinate system.

10. The method for fast trajectory optimization of many-to-one coordinated interception according to claim 9, characterized in that: In the step 5, the nonlinear programming solver is called to find all the control points of the Bezier function, and the values ​​of all state quantities and control quantities are found according to all the control points of the Bezier function; the specific process is: Step 51: Call the nonlinear programming solver to solve P1 and obtain all control points of the Bezier function; Step 52: Based on all the control points of the Bezier function, the Bezier function form (18) represented by the three-dimensional coordinates x, y, z of the interceptor missile in the launch coordinate system, as well as equations (14) and (16) are solved to obtain the values ​​of all state quantities and control quantities.

Images