Hypersonic aircraft trajectory planning method based on N-K iterative convex optimization algorithm

Abstract

The invention discloses a hypersonic aircraft trajectory planning method based on an N-K iterative convex optimization algorithm, and relates to the field of trajectory planning. The invention aims to solve the problem that an existing aircraft trajectory planning method is difficult to accurately plan an aircraft trajectory. The method comprises the following steps: obtaining a nondimensionalized hypersonic gliding section dynamic model; obtaining a linearized hypersonic gliding section dynamic model; performing discretization processing on parameters in the linearized hypersonic gliding section dynamic model to obtain a hypersonic gliding section dynamic model after parameter discretization; obtaining a control quantity constraint, a process constraint and a trust region constraint in hypersonic flight vehicle trajectory optimization; and constructing a trajectory planning problem based on iterative convex optimization by using the hypersonic glide dynamic model after parameter discretization, the control quantity constraint, the process constraint and the trust region constraint, and solving the trajectory planning problem based on the iterative convex optimization to obtain a hypersonic flight vehicle trajectory planning result. The method is used for planning the aircraft trajectory.

Description

Technical Field

[0001] The invention relates to the field of trajectory planning, and in particular to a hypersonic aircraft trajectory planning method based on a NK iterative convex optimization algorithm. Background Art

[0002] Hypersonic aircraft mainly refers to winged or wingless aircraft that fly in near-space, rely on aerodynamic forces to maintain long-term long-distance flight in near-space, and have a flight speed exceeding Mach 5. At present, the development direction of hypersonic aircraft is mainly divided into two: one is a cruise aircraft that uses a scramjet engine to achieve hypersonic flight, and the other is a glider aircraft that reenters from outside the atmosphere or near-space and relies on its own high lift-to-drag ratio to achieve long-term unpowered gliding in the atmosphere. Compared with traditional aircraft, boost-glide aircraft have the characteristics of high speed, long range, and strong maneuverability; at the same time, they also have the advantages of simple and diverse launch platforms and relatively low costs. These advantages make them have very high strategic value in the military field.

[0003] The near-space environment where the aircraft is located for a long time is complex and the atmosphere is dense. During its flight, the surface will produce intense friction with the surrounding air, and the aerodynamic heating is serious. At the same time, due to the limitations of structure and materials, and the fact that large maneuvers will have a significant impact on the terminal speed and range of the aircraft, the overload of the boost-glide aircraft during flight is greatly limited. In addition, due to the complex actual combat environment and the existence of anti-missile interception systems, the flight safety of hypersonic boost-glide aircraft is also greatly threatened. Therefore, in addition to considering the initial state constraints and terminal state constraints when planning the flight trajectory of the boost-glide aircraft, in order to ensure flight safety, it is also necessary to consider process constraints such as dynamic pressure constraints, overload constraints, heat flow constraints, and control quantity amplitude constraints, as well as the no-fly zone constraints of the interception system, which will further increase the difficulty of trajectory planning. Therefore, how to optimize the trajectory design of the aircraft under complex constraints will become a hot issue to be solved.

[0004] Aircraft trajectory optimization occupies an extremely important position in the overall design process of aircraft and is one of the key technologies of overall design. Its essence is a highly nonlinear optimal control problem. That is, to find a trajectory that can meet various constraints and can achieve the best performance indicators to describe the movement of the aircraft in space. Common performance indicators include: shortest arrival time, maximum terminal speed, optimal fuel, longest range, etc. The trajectory optimization problem is currently mainly solved by numerical methods. Numerical solutions are divided into direct methods and indirect methods. The principle of the indirect method is to use the Pondriagin maximum principle to transform the optimization problem into a two-point boundary value problem, and then obtain the optimal solution by deriving the relationship between the control quantity, state quantity and co-state variable. At present, the more commonly used indirect methods include shooting method, finite difference method and sequence gradient repair method. The direct solution method is to first transform the optimal control problem into a nonlinear programming problem using a parameterization method, and then use a numerical method to solve the transformed nonlinear programming problem. However, these methods have high computational complexity, difficult constraint processing, only local optimal solutions, and poor convergence. Convex optimization methods have become a hot topic in trajectory optimization research in recent years because they have advantages such as low requirements for the accuracy of the initial guess, the ability to solve the problem in a limited time, the local optimal solution of the convex problem is the global optimal solution, and a clear criterion for the infeasibility of the problem can be given within a limited number of iterations. However, the standard convex optimization problem only has a global optimal solution and is insensitive to the initial value. Therefore, the existing technology performs linear approximation and adds trust region constraints, so that the sequential convex optimization method can only search for the optimal solution near the expansion point, which will be very sensitive to the initial value and cannot guarantee the global optimal solution, resulting in low trajectory planning accuracy. Secondly, the subproblem after linearization is not equivalent to the original problem, and the approximation accuracy is closely related to the selection of expansion points. If the expansion point is not selected properly, the subproblem may deviate seriously from the original problem, and even fail to obtain a feasible solution, and it cannot be guaranteed that the solution obtained by solving the sequence of subproblems is convergent, resulting in the inability to achieve trajectory planning. Therefore, the aircraft trajectory planning method based on traditional convex optimization still has the problem of difficulty in accurately planning the aircraft trajectory. Summary of the invention

[0005] The purpose of the present invention is to solve the problem that the existing aircraft trajectory planning methods are still difficult to accurately plan the aircraft trajectory, and propose a hypersonic aircraft trajectory planning method based on the NK iterative convex optimization algorithm.

[0006] The trajectory planning method for hypersonic aircraft based on NK iterative convex optimization algorithm is as follows:

[0007] Step 1: Establish a hypersonic gliding stage dynamics model, and non-dimensionalize the hypersonic gliding stage dynamics model to obtain a non-dimensionalized hypersonic gliding stage dynamics model;

[0008] Step 2, linearizing the hypersonic gliding segment dynamics model after the dimensionless processing in step 1 to obtain a linearized hypersonic gliding segment dynamics model;

[0009] Step 3, discretizing the parameters in the linearized hypersonic gliding stage dynamics model to obtain the hypersonic gliding stage dynamics model after the discretization of the parameters;

[0010] Step 4: Obtain control quantity constraints, process constraints and trust region constraints in hypersonic vehicle trajectory optimization;

[0011] Step 5. Use the hypersonic gliding segment dynamics model obtained in step 3 after parameter discretization and the control quantity constraints, process constraints and trust region constraints obtained in step 4 to construct a trajectory planning problem based on iterative convex optimization, and solve the trajectory planning problem based on iterative convex optimization to obtain the hypersonic aircraft trajectory planning result.

[0012] Furthermore, the hyperglide stage dynamics model is established in step 1, and the hyperglide stage dynamics model is dimensionless to obtain the dimensionless hyperglide stage dynamics model, specifically:

[0013]

[0014]

[0015] Among them, g 0 is the acceleration due to gravity, R 0 is the radius of the earth, r is the dimensionless distance from the center of the earth, v is the dimensionless speed of the aircraft, γ is the track angle, λ is the longitude, ψ is the heading angle, φ is the latitude, Ω v ,Ω γ ,Ω ψ is an intermediate variable, ω e is the dimensionless angular velocity of the Earth's rotation, is the angular velocity of the earth's rotation before becoming dimensionless, ρ is the atmospheric density, m is the mass of the aircraft, S is the cross-sectional area of ​​the aircraft, and L is the aerodynamic lift of the aircraft after becoming dimensionless. is the aerodynamic lift of the aircraft before non-dimensionalization, C L is the aerodynamic coefficient of aerodynamic lift, is the dimensionless aerodynamic drag of the aircraft, is the dimensionless aerodynamic drag of the forward vehicle, C D is the aerodynamic coefficient of aerodynamic drag, is the time before non-dimensionalization, t' is the time after non-dimensionalization, σ is the roll angle of the aircraft, is the aircraft speed before dimensionless conversion, is the distance from the center of the Earth before being dimensionless.

[0016] Furthermore, in the step 2, the hypersonic gliding section dynamics model obtained by the dimensionless processing in step 1 is linearized to obtain the linearized hypersonic gliding section dynamics model, which is specifically:

[0017] Step 21: Design the first control variable u:

[0018] u=[u 1 ,u 2 ] T

[0019]

[0020] Among them, u 1 is the climbing force coefficient, u 2 is the turning force coefficient;

[0021] Step 22: Introduce the first control variable u designed in step 21 to obtain the aircraft dynamics model, specifically:

[0022]

[0023] a 14 = sinγ,a 15 =vcosγ

[0024]

[0025]

[0026] Among them, x = [r, λ, φ, v, γ, ψ], x is the state of the aircraft, δx is the differential of x, yes The derivative of , f(x,u) is an intermediate function, is the symbol of partial derivative, δu is the differential of u, g(x), B(x), h(x) are intermediate variables, a 14 、a 15 、a 21 、a 23 、a 24 、a 25 、a 26 、a 31 、a 34 、a 35 、a 36 、a 41 、a 44 、a 45 、a 51 、a 54 、a 55 、a 61 、a 63 、a 64 、a 65 、a65 is an intermediate variable, L r , L v is an intermediate variable, h s is the height constant;

[0027] Step 23: Based on the aircraft dynamics model obtained in step 22, introduce the second independent variable τ' and the third control variable u t , thus obtaining the linearized dynamic model of the hypersonic gliding phase.

[0028] Furthermore, the intermediate variables g(x), B(x), and h(x) are specifically:

[0029]

[0030]

[0031] in, It is the aerodynamic drag of the aircraft after non-dimensionalization.

[0032] Furthermore, in step 23, the aircraft dynamics model obtained in step 22 introduces a second independent variable τ' and a third control variable u t , thus obtaining the linearized hypersonic gliding stage dynamic model, specifically:

[0033] First, define the second independent variable τ′∈[0,1] and the duration u of the trajectory t =t f -t 0 , change the trajectory optimization problem time from the interval [t 0 ,t f ] is mapped to the interval [0,1], and the time t is in the following form:

[0034] t=t 0 +(t f -t 0 )τ'

[0035] Among them, t 0 is the starting time of trajectory planning, t f is the end time of trajectory planning;

[0036] Then, use time t to obtain the expanded dimension state The aircraft dynamics model is rewritten as:

[0037]

[0038] The aircraft dynamics model is further rewritten as:

[0039]

[0040]

[0041] in, B1 and B2 are intermediate variables. yes The derivative with respect to time t, 0 1×2 is a zero matrix;

[0042] Finally, based on the further rewritten aircraft dynamics model, the linearized hypersonic gliding phase dynamics model is obtained:

[0043]

[0044] in, yes The differential of yes The differential of is an intermediate variable.

[0045] Furthermore, the discretization processing of the parameters in the linearized hypersonic gliding segment dynamics model in step 3 to obtain the hypersonic gliding segment dynamics model after parameter discretization is specifically as follows:

[0046] Step 31: Change the trajectory optimization problem time from the interval [t 0 ,t f ] is mapped to the interval [-1,1], and the time after trajectory mapping is obtained:

[0047]

[0048] Among them, τ is the time mapped to the interval [-1,1];

[0049] Step 32: Insert N time nodes in the time interval mapped to the [-1,1] interval, use the time nodes and the end boundary points as interpolation points, obtain N+1 interpolation points, use the Lagrange interpolation polynomial to approximate the state quantity corresponding to each interpolation point, and obtain the approximate value of the aircraft state quantity, specifically:

[0050]

[0051] Among them, x(τ) is the approximate value of the aircraft state at the current time τ, X(τ) is the Lagrange interpolation polynomial value at the current time τ, and X(τ n ) is the time τ corresponding to the nth interpolation point n The Lagrange interpolation polynomial value, L n (τ) is an intermediate variable, n=0,...,N, n is the number of the interpolation point, N+1 is the total number of interpolation points, τ m is the time corresponding to the mth interpolation point, τn is the moment corresponding to the nth interpolation point;

[0052] Step 33: Based on the approximate value of the aircraft state quantity obtained in step 32, obtain the discretized hypersonic gliding segment dynamic model.

[0053] Furthermore, the step 33 obtains the discretized hypersonic gliding segment dynamics model based on the approximate value of the aircraft state quantity obtained in step 32, specifically:

[0054] First, differentiate the formula for obtaining the approximate value of the aircraft state quantity and obtain the following formula:

[0055]

[0056] Where, j = 0, ..., N-1, j is the time node number, τ i is the time corresponding to the i-th interpolation point, τ j is the time corresponding to the jth time node, i and l are the numbers of the interpolation points, and D jn (τ j ) is an intermediate variable;

[0057] Then, the dynamic model of the hypersonic gliding phase after discretization is obtained:

[0058]

[0059] Among them, f(x(τ j ),u(τ j )) is an intermediate variable, u(τ j ) is the control variable At time node τ j The value at x(τ n ), x(τ j ) is an intermediate variable;

[0060] Finally, the dynamic model of the hypersonic gliding phase after discretization is simplified as follows:

[0061]

[0062] Where Δt is the length of the equidistant interval and D is D jn (τ j ), is the first state quantity at time τ, is an intermediate variable.

[0063] Furthermore, the step 4 of obtaining the control quantity constraints, process constraints and trust region constraints in the trajectory optimization of the hypersonic vehicle is specifically as follows:

[0064] Step 41: Introduce auxiliary variable u 3 , obtain the control quantity constraints in the trajectory optimization of the hypersonic vehicle as follows:

[0065] (1) Introducing auxiliary variable u 3 , get u 3 The magnitude constraint is:

[0066] (u 3 ) min ≤u 3 ≤(u 3 ) max

[0067]

[0068] Among them, (u 3 ) min is u 3 The lower bound of (u 3 ) max is u 3 The upper bound of u 1 is the climbing force coefficient, u 2 is the turning force coefficient;

[0069] Climbing force coefficient u 1 The constraints are as follows:

[0070] u 3 cosσ max ≤u 1 ≤u 3 cosσ min

[0071] Turning force coefficient u 2 The constraints are as follows:

[0072]

[0073] (2) The constraint on the rate of change of the roll angle is as follows:

[0074]

[0075] in, is the rate of change of the roll angle, is the minimum rate of change of the roll angle, is the minimum rate of change of the roll angle, u 1 (τ n ) is the interpolation point τ n The climbing force coefficient, u 2 (τ n ) is the interpolation point τ n The turning force coefficient;

[0076] Step 42: Obtain process constraints in hypersonic vehicle trajectory optimization:

[0077] (1) Geocentric distance constraint:

[0078]

[0079] In the formula, For the The lower bound of the distance from the center of the earth corresponding to the state of the aircraft obtained in the iteration; According to The upper bound of the distance from the center of the earth corresponding to the state of the aircraft obtained in the iteration;

[0080] (2) No-fly zone restrictions:

[0081]

[0082] in, It is The state quantity obtained by the iteration corresponds to the The coordinates of the no-fly zone, p w is the coordinate of the center point of the no-fly zone, R w is the radius of the no-fly zone, N p is the number of no-fly zones, It is the no-fly zone sign;

[0083] Step 43: Set trust region constraints:

[0084]

[0085] Where η is the trust region radius, which is the same as the x dimension, It is The aircraft state quantity obtained in the iteration is x, and x is the selected aircraft state quantity.

[0086] Furthermore, in step 5, the hypersonic gliding segment dynamics model obtained in step 3 after parameter discretization and the control quantity constraints, process constraints and trust region constraints obtained in step 4 are used to construct a trajectory planning problem based on iterative convex optimization, and the trajectory planning problem based on iterative convex optimization is solved to obtain the trajectory planning result of the hypersonic aircraft, which is specifically:

[0087] First, let the optimization index function be the linear deviation index function J 1 And solve the trajectory planning problem based on iterative convex optimization to obtain the aircraft state quantity and the first control quantity

[0088] The trajectory planning problem based on iterative convex optimization is as follows:

[0089]

[0090] Then, determine whether the iteration stop condition J is met 1 <ε1, if the iteration stop condition is met, the iteration output will be stopped and If the stop iteration condition is not met, it is determined whether the preset iteration round is reached. If the preset number of iterations is reached, the iteration output is stopped. and If the preset number of iterations is not reached, update and Continue to iterate;

[0091] renew and Use the following formula:

[0092]

[0093] in, is the iteration number, is the total number of iterations, ε is the iterative convergence evaluation index value, Indicates The state quantity of the hypersonic aircraft in the iteration, It is The state quantity of the hypersonic aircraft in the iteration, δx is the differential of x, Δt is the time difference, Indicates and The first control quantity of the iteration, δu is the differential of u, is the initial state of the hypersonic vehicle, is the state of the aircraft at the end, x 0 is the initial state of the hypersonic vehicle, x f is the end state of the hypersonic vehicle, x min is the lower bound of the aircraft state, x max is the upper bound of the aircraft state, u 2 (τ n+1 ) is τ n+1 The climbing force coefficient value at 2 (τ n ) is τ n The climbing force coefficient value at 1 (τ n ) is τ n u 1 value, is the preset slack variable value, ε1 is the first iteration convergence evaluation index value;

[0094] Then, the aircraft state quantity and the first control quantity Initial value, and then let the optimization index function be the state deviation index function J 2 , and then solve the trajectory planning problem based on iterative convex optimization to obtain the trajectory planning result of the hypersonic vehicle.

[0095] Furthermore, the optimization index function is set to be the state deviation index function J 2 , and then solve the trajectory planning problem based on iterative convex optimization, specifically:

[0096] The aircraft status and the first control quantity As the optimization J 2 In and The initial value of , and then solve the following trajectory planning problem based on iterative convex optimization:

[0097]

[0098] Then, determine whether the iteration stop condition J is met 2 <ε2, if the iteration stop condition is met, the iteration is stopped and the solution is output; if the iteration stop condition is not met, it is determined whether the preset iteration round is reached. If the preset number of iterations is reached, the iteration is stopped; if the preset number of iterations is not reached, the solution is updated. and Continue iterating:

[0099]

[0100] In the formula, is the preset slack variable value, c 1 、c 2 and q are constant coefficients, For SP2 The state quantity obtained by the iterative optimization is the same as that obtained by SP1. The difference in state quantity obtained by iterative optimization, u t is the trajectory duration, and ε2 is the first iteration convergence evaluation index value.

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

[0102] The present invention firstly uses the Newton-Kantorovich method to linearize the dynamic equation of the hyperglide stage; secondly, uses Gaussian interpolation to discretize the dynamic equation of the hyperglide stage, and converts an infinite-dimensional nonlinear programming problem into a finite-dimensional linear programming problem; thirdly, convexifies the process constraints; finally, proposes a convex optimization method based on Newton-Kantorovich iteration, and uses it to solve the trajectory optimization problem. The present invention can quickly converge to a feasible solution to the original problem while eliminating the initial value sensitivity, and can ensure that the solution obtained by solving the sequence sub-problem is convergent, so the present invention can quickly and accurately plan the trajectory of the aircraft. BRIEF DESCRIPTION OF THE DRAWINGS

[0103] Figure 1 is the ground track curve of the aircraft;

[0104] Figure 2 is the speed change curve;

[0105] Figure 3 is the height change curve;

[0106] Figure 4 is the angle of attack variation curve;

[0107] Figure 5 is the roll angle variation curve;

[0108] Figure 6 is the track angle variation curve;

[0109] Figure 7 is the dynamic pressure variation curve;

[0110] Figure 8 is the heat flow variation curve;

[0111] Fig. 9 is the overload variation curve;

[0112] Fig.10 Comparison curve between integral and planned trajectory;

[0113] Fig.11 It is the comparison curve between integral and planned height. DETAILED DESCRIPTION

[0114] Specific implementation method 1: The specific process of the hypersonic vehicle trajectory planning method based on the NK iterative convex optimization algorithm in this implementation method is as follows:

[0115] Step 1: Establish a hypersonic gliding stage dynamics model, and make the hypersonic gliding stage dynamics model dimensionless to obtain a dimensionless hypersonic gliding stage dynamics model, specifically:

[0116] First, when using the convex optimization method for trajectory optimization, the influence of the J2 perturbation term is ignored. At this time, the dynamic model of the aircraft in the half-speed system can be written as follows:

[0117]

[0118] in, All about time The derivative of is the distance from the center of the earth, λ is the longitude, φ is the latitude, is the unit speed, γ is the track angle, ψ is the heading angle, is the angular velocity of the Earth's rotation, R 0 =6378km is the radius of the earth, ρ is the atmospheric density, ρ 0 is the standard sea level atmospheric density, h is the current altitude of the aircraft, and h s =7110m is the height constant, is the aerodynamic drag of the aircraft, is the aerodynamic lift of the aircraft, σ is the roll angle of the aircraft, S is the cross-sectional area of ​​the aircraft, C L is the aerodynamic coefficient of aerodynamic lift, C D is the aerodynamic coefficient of aerodynamic drag, C L , C D It can be expressed as a function of angle of attack and Mach number.

[0119] Then, in order to ensure the calculation accuracy, the dynamic equation of the unpowered section is dimensionless, and the distance from the center of the earth before the variable is dimensionless is used. Velocity before dimensionless conversion The angular velocity of the Earth's rotation before non-dimensionalization Time before dimensionless After dimensionless processing, the dynamic model formula (1) of the aircraft is rewritten as follows:

[0120]

[0121] Among them, g 0 is the acceleration due to gravity, R 0 =6378km is the radius of the earth, r is the dimensionless distance from the center of the earth, v is the dimensionless speed of the aircraft, γ is the track angle, λ is the longitude, ψ is the heading angle, φ is the latitude, Ω v ,Ω γ ,Ω ψ is an intermediate variable, ω e is the dimensionless angular velocity of the Earth's rotation, is the angular velocity of the earth's rotation before becoming dimensionless, ρ is the atmospheric density, m is the mass of the aircraft, S is the cross-sectional area of ​​the aircraft, and L is the aerodynamic lift of the aircraft after becoming dimensionless. is the aerodynamic lift of the aircraft before non-dimensionalization, C L is the aerodynamic coefficient of aerodynamic lift, is the dimensionless aerodynamic drag of the aircraft, is the dimensionless aerodynamic drag of the forward vehicle, C D is the aerodynamic coefficient of aerodynamic drag, is the time before dimensionless conversion, t' is the time after dimensionless conversion, is the aircraft speed before dimensionless conversion, is the distance from the center of the Earth before being dimensionless.

[0122] Step 2: linearizing the hypersonic gliding segment dynamics model after the dimensionless processing in step 1 to obtain a linearized hypersonic gliding segment dynamics model, including the following steps:

[0123] At present, most algorithms for solving nonlinear differential equations are Newton's method or variants based on Newton's method. When solving a set of nonlinear differential equations, a set of equivalent algebraic equations are usually first obtained by using discretization and polynomial approximation. Then, Newton's method or variants based on Newton's method are used to iterate from a given initial value until a numerical solution that satisfies a certain accuracy is obtained. Using Newton's method requires calculating the Jacobian matrix of the equation group, and the calculation of the Jacobian matrix is ​​cumbersome and time-consuming. In order to reduce the complexity of calculating the Jacobian matrix, the present invention directly applies Newton's method to the original differential equation group, rather than to the algebraic equation group obtained using discretization and polynomial approximation methods. This method of directly applying Newton's method to differential equations and iteratively linearizing nonlinear differential equations is called the "Newton-Kantorovich (NK)" method.

[0124] Without loss of generality, assume that there exists a first-order differential equation:

[0125]

[0126] Assume there is an iteration variable s (k) , which is an approximation of the true solution s of the differential equation. Then, Taylor expansion of equation (15) at s yields:

[0127]

[0128] Among them, f s represents the derivative of the function f(s) with respect to s;

[0129] The recursive formula for the next iteration is:

[0130]

[0131] Among them, f and f s Both are s=s(k) The value at time;

[0132] Define the iteration variable s (k) The difference from the true solution s of the equation is δs:

[0133] δs=ss (k) (18)

[0134] The Frechet derivative is defined as:

[0135]

[0136] in, is a scaling factor, F(s) is a nonlinear operator including differentials or partial differentials, and F s δs is the Frechet derivative of the nonlinear operator F(s) with respect to s in the direction of δs, F s is a linear operator;

[0137] By the above definition, the scaling factor can be Transition from complex functional analysis to the simpler realm of solving calculus. The importance of the Frechet derivative is that it gives the linear terms in the generalized Taylor expansion for any nonlinear operator F(s), as follows:

[0138] F(s+δs)=F(s)+F s (s)δs+o[(δs) 2 ](20)

[0139]

[0140] Combining equation (19), equation (20) and equation (21), we get:

[0141]

[0142] In the unpowered flight phase, the angle of attack and the bank angle are the control variables in the trajectory optimization problem. It can be observed that the angle of attack is implicit in the aerodynamic force in the unpowered phase dynamic equation, and the bank angle exists in the form of a trigonometric function. Therefore, if the angle of attack and the bank angle are directly used as control variables, the corresponding control matrix must be strongly nonlinear, which will greatly increase the difficulty of linearizing the dynamic equations, and after linear post-processing, the solved control variables are also prone to "jittering". In order to simplify the problem, consider introducing other variables as new control variables. The newly selected control variable should meet the following two conditions: 1) The control matrix of the new control variable can be written in a simple linear form about the state quantity; 2) There is a unique mapping relationship between the new control variable and the original control variable angle of attack and bank angle, which can well reflect the changes in the angle of attack and bank angle.

[0143] Step 21: Introduce the first control quantity u and the component climbing force coefficient u 1 , turning force coefficient u 2 :

[0144] u=[u 1 ,u 2 ] T (twenty three)

[0145]

[0146] Step 22: Obtain the aircraft dynamics model based on the first control variable u introduced in step 21:

[0147] The aircraft dynamics model is as follows:

[0148]

[0149]

[0150] Among them, g(x), B(x), h(x) are intermediate variables, f(x,u) is an intermediate function, and x = [r, λ, φ, v, γ, ψ] is the state of the aircraft;

[0151] make According to the deduction:

[0152] F(x+δx,u+δu)=F(x,u)+F x (x,u)δx+F u (x,u)δu(26)

[0153]

[0154] Among them, δu is the differential of u, δx is the differential of x, yes The derivative of x (x,u)δx、F u (x,u)δu is an intermediate variable, f(x+λδx,u) and f(x,u+λδu) are obtained by formula (25), is the scaling factor;

[0155] Substituting formula (27) and (28) into formula (26), the aircraft dynamics model is rewritten as follows:

[0156]

[0157]

[0158] in, L r , L vis an intermediate variable.

[0159] The trajectory optimization problem studied in this invention is considered to be a two-point boundary value problem with a fixed initial state and terminal state and free time. If the terminal time needs to be controlled, the method of expanding the dimension of the dynamic equation can be used to treat time as a new state quantity and the terminal time as a new control quantity.

[0160] Step 2: Introduce a new second independent variable τ' and a third control variable u t , thus obtaining the linearized dynamic model of the hypersonic gliding phase:

[0161] First, define a new second independent variable τ′∈[0,1] and the duration u of the trajectory t =t f -t 0 , change the trajectory optimization problem time from the interval [t 0 ,t f ] is mapped to the interval [0,1], and the time t is in the following form:

[0162] t=t 0 +(t f -t 0 )τ'(30)

[0163] Among them, t 0 is the starting time of trajectory planning, t f is the end time of trajectory planning;

[0164] Then, the converted time t is used as the new state quantity to expand the system equation and obtain the expanded state quantity At this time, the aircraft dynamics model can be written as:

[0165]

[0166] Similarly, the aircraft dynamics model can also be written as follows:

[0167]

[0168] in, B1 and B2 are intermediate variables. yes The derivative with respect to time t, 0 1×2 is a zero matrix;

[0169] Thus, the linearized hypersonic gliding stage dynamic model is obtained as follows:

[0170]

[0171]

[0172] in, yes The differential of yes The differential of is an intermediate variable.

[0173] Step 3: discretizing the parameters in the linearized hypersonic gliding segment dynamics model to obtain the hypersonic gliding segment dynamics model after parameter discretization, including the following steps:

[0174] Step 31: Change the trajectory optimization problem time from the interval [t 0 ,t f ] is mapped to the interval [-1,1], and the time after trajectory mapping is obtained:

[0175]

[0176] Among them, τ is the time mapped to the interval [-1,1]; τ = -1 corresponds to the trajectory start time t 0 , when τ = 1, the corresponding trajectory end time is t f ;

[0177] Step 32: Insert N time nodes (LG nodes) in the time interval mapped to the [-1,1] interval to obtain N+1 interpolation points (time nodes and end boundary points). Use the Lagrange interpolation polynomial to approximate the state quantity corresponding to each interpolation point to obtain the approximate value of the aircraft state quantity, specifically:

[0178]

[0179] Among them, x(τ) is the approximate value of the aircraft state at the current time τ, X(τ) is the Lagrange interpolation polynomial value at the current time τ, and X(τ n ) is the time τ corresponding to the nth interpolation point n The Lagrange interpolation polynomial value, L n (τ) is the intermediate variable, n is the number of the interpolation point, N+1 is the total number of interpolation points, τ m is the time corresponding to the mth interpolation point, τ n is the time corresponding to the nth interpolation point, n = 0, ..., N;

[0180] Step 33: Based on the approximate value of the aircraft state quantity obtained in step 32, a discretized hypersonic gliding segment dynamics model is obtained, specifically:

[0181] First, differentiate formula (35) to obtain the following formula:

[0182]

[0183] Where, j = 0, ..., N-1, j is the time node number, τ i is the time corresponding to the i-th interpolation point, τ j is the time corresponding to the jth time node, i and l are the numbers of the interpolation points, and D jn (τ j ) is an intermediate variable;

[0184] Then, the dynamic model of the hypersonic gliding phase after discretization is obtained:

[0185]

[0186] Formula (37) can be simplified as follows:

[0187]

[0188] Among them, u(τ j ) is the control variable At time node τ j , Δt is the length of the equidistant interval (the time difference of each small segment), and D is D jn (τ j ), x(τ n ), x(τ j ) is the intermediate variable obtained according to formula (35), is the first state quantity at time τ, f(x(τ j ),u(τ j ))、f(x(τ), According to formula (25),

[0189] The control variables exist directly in the kinematic equations and do not exist in differential form, so there is no need to use interpolation polynomials for approximation.

[0190] Step 4: Convexify the control quantity constraints and process constraints in the trajectory optimization of the hypersonic vehicle, and introduce the trust region constraints:

[0191] Step 41: Introduce auxiliary variable u 3 , obtain the control quantity constraints in the trajectory optimization of hypersonic aircraft:

[0192] Introducing auxiliary variables Auxiliary variable u 3 It is only related to the angle of attack and the state of the aircraft, and has nothing to do with the roll angle. Determine u according to the state quantity corresponding to the trajectory at the interpolation node and the value range of the angle of attack. 3 The value range of u 3 The magnitude constraint can be expressed as:

[0193] (u 3 ) min ≤u 3 ≤(u 3 ) max (38)

[0194] Among them, (u 3 ) min is u 3 The lower bound of (u 3 ) max is u 3 The upper bound of u 1 is the climbing force coefficient, u 2 is the turning force coefficient;

[0195] u 1 The value range of can be determined according to the auxiliary variable u 3 The value range of the roll angle is determined by:

[0196] u 3 cosσ max ≤u 1 ≤u 3 cosσ min (39)

[0197] Among them, σ max is the maximum tilt angle in the constraint, σ min is the minimum roll angle in the constraint;

[0198] Using relaxation techniques, we can obtain 2 Constraints:

[0199]

[0200] In the process of trajectory optimization, in addition to constraining the amplitude of the angle of attack and the roll angle, it is usually necessary to constrain the rate of change of the roll angle:

[0201] u 2 Taking the derivative with respect to time, we get:

[0202] u 2 ′=(u 3 sinσ)′=u 3 cosσσ′=u 1 σ′ (41)

[0203] Where σ′ is the derivative of the roll angle σ;

[0204] At this time there are:

[0205]

[0206] in, is the rate of change of the roll angle, is the minimum rate of change of the roll angle, is the minimum rate of change of the roll angle, u 1 (τ n ) is the interpolation point τ n The climbing force coefficient, u 2 (τ n ) is the interpolation point τ n The turning force coefficient;

[0207] Through the above processing, the initial constraints on the angle of attack and the roll angle are transformed into constraints on the newly introduced control quantities. And the constraints obtained through the above process are all convex constraints.

[0208] Step 42: Convexify the process constraints in the trajectory optimization of the hypersonic vehicle:

[0209] Dynamic pressure constraint, overload constraint, heat flow constraint, balanced glide constraint and no-fly zone constraint are all nonlinear functions of state quantities, and they are not convex functions. Therefore, when using convex optimization algorithm to solve trajectory optimization problem, it is necessary to convexify it. Now based on the height and speed expression of dynamic pressure constraint, overload constraint and heat flow constraint and the balanced glide constraint expression, it is converted into a constraint on the distance from the center of the earth, which can be abbreviated as:

[0210]

[0211] In the formula, For the The lower bound of the distance from the center of the earth corresponding to the state of the aircraft obtained in the iteration; According to The upper bound of the distance from the center of the earth corresponding to the state of the aircraft obtained in the iteration;

[0212] The no-fly zone constraint is also non-convex. It is linearized and transformed into an affine function about the state quantity. Taking a no-fly zone as an example, the specific form is:

[0213]

[0214] in, It is The state quantity obtained by the iteration corresponds to the The coordinates of the no-fly zone, p w is the coordinate of the center point of the no-fly zone, R w is the radius of the no-fly zone, N p is the number of no-fly zones, It is the no-fly zone sign;

[0215] Step 4.3: In order to reduce the error after linearization of the dynamic equation, the trust region constraint is introduced, which is expressed as:

[0216]

[0217] Where η is the trust region radius, which is the same as the x dimension, It is The aircraft state quantity obtained by the iteration, x is the aircraft state quantity selected manually;

[0218] Step 5: Using the hypersonic glide segment dynamics model obtained in step 3 after parameter discretization and the constraints obtained in step 4, a trajectory planning problem based on iterative convex optimization is constructed, and the trajectory planning problem based on iterative convex optimization is solved to obtain the trajectory planning result of the hypersonic aircraft, specifically:

[0219] Step 51: Taylor expansion is performed on the hypersonic gliding segment dynamics model obtained in step 3 after the parameters are discretized, and based on the hypersonic gliding segment dynamics model obtained in step 4 after the parameters are discretized after Taylor expansion, a trajectory planning problem based on iterative convex optimization is constructed, specifically:

[0220] When the conditions are met When , solve the following problem to minimize J:

[0221]

[0222] Then, determine whether the iteration stop condition is met If the iteration stop condition is met, the iteration is stopped; if the iteration stop condition is not met, it is determined whether the preset iteration round is reached. If the preset number of iterations is reached, the iteration is stopped; if the preset number of iterations is not reached, the iteration is continued according to the following content:

[0223]

[0224] in, is the iteration number, is the total number of iterations, It is The change in the state of the hypersonic aircraft in the iteration, ε is a small quantity, is the iterative convergence evaluation index value, J is the optimization index function, Indicates The state quantity of the hypersonic aircraft in the iteration, It is The state quantity of the hypersonic aircraft in the iteration, It is The state quantity of the hypersonic aircraft in the iteration, δx is the differential of x, Δt is the time difference, Indicates and The first control quantity of the iteration, u is the first control quantity, δu is the differential of u, is the initial state of the hypersonic vehicle, is the state of the aircraft at the end, x 0 is the initial state of the hypersonic vehicle, x f is the end state of the hypersonic vehicle, x min is the lower bound of the aircraft state, x max is the upper bound of the aircraft state, u 2 (τ n+1 ) is τ n+1 The climbing force coefficient value at 2 (τ n ) is τ n The climbing force coefficient value at 1 (τ n ) is τ n u 1 value.

[0225] Step 52: The constraints that need to be considered in the gliding phase of the aircraft mainly include: initial state constraint, terminal state constraint, dynamic pressure constraint, overload constraint, thermal flow constraint, no-fly zone constraint, control quantity constraint, etc. According to the trajectory planning problem formula based on iterative convex optimization obtained in step 51, the time optimality is used as the optimization index to solve the trajectory optimization problem and obtain the trajectory planning result of the hypersonic aircraft:

[0226] First, let the optimization index function be the linear deviation index function J 1 , and solve the trajectory planning problem based on iterative convex optimization to obtain and The details are as follows:

[0227]

[0228] Then, determine whether the iteration stop condition J is met 1 <ε1, if the iteration stop condition is met, then stop the iteration and output x (k) and If the stop iteration condition is not met, determine whether the preset iteration round has been reached. If the preset number of iterations has been reached, stop the iteration. If the preset number of iterations has not been reached, update and Continue iterating:

[0229]

[0230] Then, the aircraft state quantity and the first control quantity Initial value, and then let the optimization index function be the state deviation index function J 2, and then solve the trajectory planning problem based on iterative convex optimization to obtain the trajectory planning result of the hypersonic vehicle, which is specifically:

[0231]

[0232] Then, determine whether the iteration stop condition J is met 2 <ε2, if the iteration stop condition is met, the iteration is stopped and the solution is output; if the iteration stop condition is not met, it is determined whether the preset iteration round is reached. If the preset number of iterations is reached, the iteration is stopped; if the preset number of iterations is not reached, the solution is updated. and Continue iterating:

[0233]

[0234] In the formula, is the preset slack variable value, c 1 、c 2 and q are constant coefficients, For SP2 The state quantity obtained by the iterative optimization is the same as that obtained by SP1. The difference in state quantity obtained by iterative optimization, u t is the trajectory duration, and ε2 is the first iteration convergence evaluation index value.

[0235] Embodiment: In order to verify the beneficial effects of the present invention, the present invention conducted the following simulation experiments:

[0236] The purpose of the first step SP1 is to ensure that the linearization error of the dynamic equation is relatively small and to obtain a feasible solution that preliminarily satisfies the constraints under a set of poor initial values.

[0237] The second step SP2 introduces optimization indicators and ensures the convergence speed. To ensure continuity, the terminal state of the planned pull-up segment is used as the initial value of the gliding segment optimization. The simulation parameters are shown in Table 1:

[0238] Table 1 Simulation parameters

[0239]

[0240] The trust region constraints are set as:

[0241]

[0242] The convergence condition is set as:

[0243]

[0244] The simulation results are as follows Figure 1-Figure 11As shown in the figure, it can be seen from the simulation results that the trajectory optimization based on the convex optimization method of Newton-Kantorovich iteration can obtain a trajectory that satisfies complex constraints such as dynamic pressure constraints, overload constraints, thermal flow constraints, no-fly zone constraints, initial state and terminal state. The initial guess of the gliding time is 1156.8s, and the optimized time is 1053s, indicating that the optimal solution based on time optimization is obtained. In addition, the angle of attack and the roll angle keep changing continuously throughout the process without sudden changes, which can reduce the control difficulty of the control system.

[0245] It can be seen that after interpolating the control quantity generated by the convex optimization method, the result obtained by integration can be well consistent with the result produced by convex optimization, so the accuracy of the result solved by this method can be proved.

Claims

1. A hypersonic vehicle trajectory planning method based on NK iterative convex optimization algorithm, characterized in that The specific process of the method is: Step 1: Establish a hypersonic gliding stage dynamics model, and non-dimensionalize the hypersonic gliding stage dynamics model to obtain a non-dimensionalized hypersonic gliding stage dynamics model; Step 2, linearizing the hypersonic gliding segment dynamics model after the dimensionless processing in step 1 to obtain a linearized hypersonic gliding segment dynamics model; Step 3, discretizing the parameters in the linearized hypersonic gliding stage dynamics model to obtain the hypersonic gliding stage dynamics model after parameter discretization; Step 4: Obtain control quantity constraints, process constraints and trust region constraints in hypersonic vehicle trajectory optimization; Step 5. Use the hypersonic gliding segment dynamics model obtained in step 3 after parameter discretization and the control quantity constraints, process constraints and trust region constraints obtained in step 4 to construct a trajectory planning problem based on iterative convex optimization, and solve the trajectory planning problem based on iterative convex optimization to obtain the hypersonic aircraft trajectory planning result.

2. The hypersonic vehicle trajectory planning method based on the NK iterative convex optimization algorithm according to claim 1 is characterized in that: The step 1 of establishing a hyperglide stage dynamics model and non-dimensionalizing the hyperglide stage dynamics model to obtain a non-dimensionalized hyperglide stage dynamics model is specifically as follows: Among them, g0 is the acceleration of gravity, R0 is the radius of the earth, r is the dimensionless distance from the center of the earth, v is the dimensionless speed of the aircraft, γ is the track angle, λ is the longitude, ψ is the heading angle, φ is the latitude, Ω v ,Ω γ ,Ω ψ is an intermediate variable, ω e is the dimensionless angular velocity of the Earth's rotation, is the angular velocity of the earth's rotation before becoming dimensionless, ρ is the atmospheric density, m is the mass of the aircraft, S is the cross-sectional area of ​​the aircraft, and L is the aerodynamic lift of the aircraft after becoming dimensionless. is the aerodynamic lift of the aircraft before non-dimensionalization, C L is the aerodynamic coefficient of aerodynamic lift, is the dimensionless aerodynamic drag of the aircraft, is the dimensionless aerodynamic drag of the forward vehicle, C D is the aerodynamic coefficient of aerodynamic drag, is the time before non-dimensionalization, t' is the time after non-dimensionalization, σ is the roll angle of the aircraft, is the aircraft speed before non-dimensionalization, is the distance from the center of the Earth before dimensioning.

3. The hypersonic vehicle trajectory planning method based on the NK iterative convex optimization algorithm according to claim 2 is characterized in that: The step 2 is to linearize the hypersonic gliding segment dynamics model after the dimensionless processing in step 1 to obtain the linearized hypersonic gliding segment dynamics model, specifically: Step 21: Design the first control variable u: u=[u1,u2] T Among them, u1 is the climbing force coefficient, u2 is the turning force coefficient; Step 22: Introduce the first control variable u designed in step 21 to obtain the aircraft dynamics model, specifically: a 14 =sinγ,a 15 =vcosγ Among them, x = [r, λ, φ, v, γ, ψ], x is the state of the aircraft, δx is the differential of x, yes The derivative of , f(x,u) is an intermediate function, is the symbol of partial derivative, δu is the differential of u, g(x), B(x), h(x) are intermediate variables, a 14 、a 15 、a 21 、a 23 、a 24 、a 25 、a 26 、a 31 、a 34 、a 35 、a 36 、a 41 、a 44 、a 45 、a 51 、a 54 、a 55 、a 61 、a 63 、a 64 、a 65 、a 65 is an intermediate variable, L r , L v is an intermediate variable, h s is the height constant; Step 23: Based on the aircraft dynamics model obtained in step 22, introduce the second independent variable τ' and the third control variable u t , thus obtaining the linearized dynamic model of the hypersonic gliding phase.

4. The hypersonic vehicle trajectory planning method based on the NK iterative convex optimization algorithm according to claim 3 is characterized in that: The intermediate variables g(x), B(x), and h(x) are specifically: in, It is the aerodynamic drag of the aircraft after non-dimensionalization.

5. The hypersonic vehicle trajectory planning method based on the NK iterative convex optimization algorithm according to claim 4 is characterized in that: In step 23, the second independent variable τ' and the third control variable u are introduced based on the aircraft dynamics model obtained in step 22. t , thus obtaining the linearized hypersonic gliding stage dynamic model, specifically: First, define the second independent variable τ′∈[0,1] and the duration u of the trajectory t =t f -t0, change the trajectory optimization problem time from the interval [t0,t f ] is mapped to the interval [0,1], and the time t is in the following form: t=t0+(t f -t0)τ' Among them, t0 is the starting time of trajectory planning, t f is the end time of trajectory planning; Then, use time t to obtain the expanded dimension state The aircraft dynamics model is rewritten as: The aircraft dynamics model is further rewritten as: in, B1 and B2 are intermediate variables. yes The derivative with respect to time t, 0 1×2 is a zero matrix; Finally, based on the further rewritten aircraft dynamics model, the linearized hypersonic gliding phase dynamics model is obtained: in, yes The differential of yes The differential of is an intermediate variable.

6. The hypersonic vehicle trajectory planning method based on the NK iterative convex optimization algorithm according to claim 5 is characterized in that: The step 3 is to discretize the parameters in the linearized hypersonic gliding segment dynamics model to obtain the hypersonic gliding segment dynamics model after the discretization of the parameters, specifically: Step 31: Change the trajectory optimization problem time from the interval [t0,t f ] is mapped to the interval [-1,1], and the time after trajectory mapping is obtained: Among them, τ is the time mapped to the interval [-1,1]; Step 32: Insert N time nodes in the time interval mapped to the [-1,1] interval, use the time nodes and the end boundary points as interpolation points, obtain N+1 interpolation points, use the Lagrange interpolation polynomial to approximate the state quantity corresponding to each interpolation point, and obtain the approximate value of the aircraft state quantity, specifically: Among them, x(τ) is the approximate value of the aircraft state at the current time τ, X(τ) is the Lagrange interpolation polynomial value at the current time τ, and X(τ n ) is the time τ corresponding to the nth interpolation point n The Lagrange interpolation polynomial value, L n (τ) is an intermediate variable, n=0,...,N, n is the number of the interpolation point, N+1 is the total number of interpolation points, τ m is the time corresponding to the mth interpolation point, τ n is the moment corresponding to the nth interpolation point; Step 33: Based on the approximate value of the aircraft state quantity obtained in step 32, obtain the discretized hypersonic gliding segment dynamic model.

7. The hypersonic vehicle trajectory planning method based on the NK iterative convex optimization algorithm according to claim 6 is characterized in that: The step 33 obtains the discretized hypersonic gliding segment dynamics model based on the approximate value of the aircraft state quantity obtained in step 32, specifically: First, the formula for obtaining the approximate value of the aircraft state quantity is differentiated to obtain the following formula: Where, j = 0, ..., N-1, j is the time node number, τ i is the time corresponding to the i-th interpolation point, τ j is the time corresponding to the jth time node, i and l are the numbers of the interpolation points, and D jn (τ j ) is an intermediate variable; Then, the dynamic model of the hypersonic gliding phase after discretization is obtained: Among them, f(x(τ j ),u(τ j )) is an intermediate variable, u(τ j ) is the control variable At time node τ j The value at x(τ n ), x(τ j ) is an intermediate variable; Finally, the dynamic model of the hypersonic gliding phase after discretization is simplified as follows: Where Δt is the length of the equidistant interval and D is D jn (τ j ), is the first state quantity at time τ, is an intermediate variable.

8. The hypersonic vehicle trajectory planning method based on the NK iterative convex optimization algorithm according to claim 7 is characterized in that: The step 4 of obtaining the control quantity constraints, process constraints and trust region constraints in the trajectory optimization of the hypersonic vehicle is specifically as follows: Step 41: Introduce auxiliary variable u3 to obtain the control quantity constraint in the trajectory optimization of the hypersonic vehicle, as follows: (1) Introduce the auxiliary variable u3 and obtain the amplitude constraint of u3 as follows: (u3) min ≤u3≤(u3) max Among them, (u3) min is the lower bound of u3, (u3) max is the upper bound of u3, u1 is the climbing force coefficient, and u2 is the turning force coefficient; The constraints on the climbing force coefficient u1 are as follows: u3 cosσ max ≤u1≤u3 cosσ min The constraints on the turning force coefficient u2 are as follows: (2) The constraint on the rate of change of the roll angle is as follows: in, is the rate of change of the roll angle, is the minimum rate of change of the roll angle, is the minimum rate of change of the roll angle, u1(τ n ) is the interpolation point τ n The climbing force coefficient, u2(τ n ) is the interpolation point τ n The turning force coefficient; Step 42: Obtain process constraints in hypersonic vehicle trajectory optimization: (1) Geocentric distance constraint: In the formula, For the The lower bound of the distance from the center of the earth corresponding to the state of the aircraft obtained in the iteration; According to The upper bound of the distance from the center of the earth corresponding to the state of the aircraft obtained in the iteration; (2) No-fly zone restrictions: in, It is The state quantity obtained by the iteration corresponds to the The coordinates of the no-fly zone, p w is the coordinate of the center point of the no-fly zone, R w is the radius of the no-fly zone, N p is the number of no-fly zones, It is the no-fly zone sign; Step 43: Set trust region constraints: Where η is the trust region radius, which is the same as the x dimension, It is The aircraft state quantity obtained in the iteration is x, and x is the selected aircraft state quantity.

9. The hypersonic vehicle trajectory planning method based on the NK iterative convex optimization algorithm according to claim 8, characterized in that: In the step 5, the hypersonic gliding segment dynamics model obtained in the step 3 after the parameter discretization and the control quantity constraint, process constraint and trust region constraint obtained in the step 4 are used to construct a trajectory planning problem based on iterative convex optimization, and the trajectory planning problem based on iterative convex optimization is solved to obtain the trajectory planning result of the hypersonic aircraft, which is specifically: First, let the optimization index function be the linear deviation index function J1 and solve the trajectory planning problem based on iterative convex optimization to obtain the aircraft state quantity and the first control quantity The trajectory planning problem based on iterative convex optimization is as follows: Then, determine whether the iteration stop condition J1 < ε1 is met. If the iteration stop condition is met, stop the iteration output. and If the stop iteration condition is not met, it is determined whether the preset iteration round is reached. If the preset number of iterations is reached, the iteration output is stopped. and If the preset number of iterations is not reached, update and Continue to iterate; renew and Use the following formula: in, is the iteration number, is the total number of iterations, ε is the iterative convergence evaluation index value, Indicates The state quantity of the hypersonic aircraft in the iteration, It is The state quantity of the hypersonic aircraft in the iteration, δx is the differential of x, Δt is the time difference, Indicates and The first control quantity of the iteration, δu is the differential of u, is the initial state of the hypersonic vehicle, is the state of the aircraft at the end, x0 is the initial state of the hypersonic aircraft, and x f is the end state of the hypersonic vehicle, x min is the lower bound of the aircraft state, x max is the upper bound of the aircraft state, u2(τ n+1 ) is τ n+1 The climbing force coefficient value at n ) is τ n The climbing force coefficient value at n ) is τ n The u1 value at is the preset slack variable value, ε1 is the first iteration convergence evaluation index value; Then, the aircraft state quantity and the first control quantity The initial value is taken, and then the optimization index function is set to be the state deviation index function J2. Then, the trajectory planning problem based on iterative convex optimization is solved to obtain the trajectory planning result of the hypersonic aircraft.

10. The hypersonic vehicle trajectory planning method based on the NK iterative convex optimization algorithm according to claim 9, characterized in that: The optimization index function is set as the state deviation index function J2, and then the trajectory planning problem based on iterative convex optimization is solved, specifically: The aircraft status and the first control quantity As an optimization in J2 and The initial value of , and then solve the following trajectory planning problem based on iterative convex optimization: Then, determine whether the iteration stop condition J2<ε2 is met. If the iteration stop condition is met, stop the iteration and output the solution result; If the stop iteration condition is not met, determine whether the preset iteration round is reached. If the preset number of iterations is reached, stop the iteration. If the preset number of iterations is not reached, update and Continue iterating: In the formula, is the preset slack variable value, c1, c2 and q are constant coefficients, For SP2 The state quantity obtained by the iterative optimization is the same as that obtained by SP1. The difference in state quantity obtained by iterative optimization, u t is the trajectory duration, and ε2 is the first iteration convergence evaluation index value.

Images