A no-fly zone avoidance guidance method for hypersonic aircraft based on pseudospectral method

Abstract

The invention discloses a hypersonic aircraft no-fly zone avoidance guidance method based on pseudo-spectral method, including: establishing a hypersonic aircraft unpowered reentry guidance dynamics equation; establishing an ambient atmospheric density and gravity acceleration model; analyzing a hypersonic aircraft aerodynamic model; establishing a hypersonic aircraft heat flow, dynamic pressure, overload, and no-fly zone and other process constraint models; according to the above models and constraints, discretizing the reentry process state variables and control variables based on the Radau pseudo-spectral principle; for a constructed optimization problem, using the fourth-order Ronge-Kutta method to integrate and solve the initial state variables and control variable estimates; based on discrete states and control variables, taking minimum heat flow and trajectory oscillation as optimization targets, establishing a nonlinear programming problem, and the solution process adaptively corrects the discrete state according to the accuracy requirements. The present invention can enhance the success rate of avoiding multiple no-fly zones in a complex environment during the reentry guidance process.

Description

Technical Field

[0001] The present invention relates to the field of aerospace technology, and in particular to a hypersonic aircraft no-fly zone avoidance guidance method based on pseudo-spectral method. Background Art

[0002] As an important guarantee of my country's national defense industry defense system, hypersonic aircraft is a national weapon with an unshakable position. Generally, hypersonic aircraft refers to aircraft with a flight speed greater than 5Ma, which adopts a rotating body or wave-riding body with a relatively large lift-drag ratio. It is dropped from a space combat aircraft at high altitude or boosted to a predetermined altitude by a rocket booster. After that, the gliding warhead relies on aerodynamic lift to achieve long-distance non-ballistic reentry flight, and has the ability to reach long-range quickly, thus realizing long-range precision strikes. The flight span of this aircraft realizes long-range high-speed flight in the atmosphere and across the atmosphere, that is, below 20km and at an altitude of 20 to 100km. It has many advantages such as fast flight speed, good concealment, and strong penetration capability. It can greatly expand the battlefield space and has an important deterrent effect in the military. Future wars reflect the complexity, informatization and intelligence of the battlefield environment, and air strikes will rely on speed and altitude to win. Hypersonic aircraft with ultra-high altitude and ultra-high speed give it the following advantages:

[0003] (1) Wide combat range, able to reach the entire world in a short time. Hypersonic vehicles can reach any point on the earth within 2 hours, and can therefore quickly strike various military targets thousands or tens of thousands of kilometers away, greatly expanding the battlefield space.

[0004] (2) Good concealment and difficult to detect. With its extremely fast flight speed and super maneuverability, hypersonic aircraft can quickly pass through the combat airspace of the opponent's defense system, and the number of accumulated echoes is small, making it much more difficult for ground air defense weapon systems to detect it, which is convenient for intelligence reconnaissance and battlefield surveillance.

[0005] (3) Strong penetration and survivability. On the basis of good concealment, hypersonic flight also makes the opponent's interception weapons have insufficient reaction time and flight speed to implement interception, so the penetration capability is extremely strong and can achieve unexpected attack effects. In addition, its high flight altitude, extremely fast flight speed and strong lateral maneuverability make it impossible for any current air defense weapon to hit it. If stealth technology is applied to hypersonic aircraft, its concealment will be better, and its penetration and survivability will be stronger.

[0006] (4) The range is relatively long and the combat effectiveness is high. At present, the range of hypersonic missiles being studied abroad is hundreds of kilometers or even thousands of kilometers. They use high-speed flight to obtain huge kinetic energy and kill targets by direct collision. Their combat effectiveness is extremely high. If the warhead is designed to be of the same mass as that of a subsonic aircraft, its warhead will be more powerful.

[0007] (5) Significant strategic deterrence. The significant military advantages of hypersonic vehicles will surely make them a new generation of aerospace "trump card" with strong strategic deterrence.

[0008] However, the flight environment of hypersonic vehicles is extremely complex and subject to multiple constraints, especially the no-fly zones formed by geopolitical or military monitoring factors during the mission, which further increases the complexity of the reentry trajectory optimization process. The reentry process of hypersonic vehicles is already subject to multiple process constraints, and due to the characteristics of unpowered gliding, its energy needs to be carefully managed to ensure the subsequent strike and destructive capabilities. The emergence of no-fly zones will further restrict the reentry trajectory planning, reduce the feasible flight area of ​​the aircraft, and pose further challenges to the reentry process of the aircraft.

[0009] Therefore, how to design a reentry state and control quantity sequence that can generate a path under multiple constraints that can satisfy both process constraints and terminal constraints, reach the designated area, pass through a complex battlefield environment, and satisfy smaller heat flow indicators and trajectory oscillation conditions has become a problem that needs to be solved urgently. Summary of the invention

[0010] The purpose of the present invention is to provide a hypersonic aircraft no-fly zone avoidance guidance method based on pseudospectral method to solve the above problems.

[0011] The present invention solves the technical problem by adopting the following technical solution:

[0012] A hypersonic aircraft no-fly zone avoidance guidance method based on pseudospectral method comprises the following steps:

[0013] S1. Establish the unpowered reentry guidance dynamic equations of hypersonic vehicles;

[0014] S2. Establish the ambient atmospheric density and gravity acceleration model;

[0015] S3. Analyze the aerodynamic model of hypersonic vehicle based on S1 and S2;

[0016] S4. Establish process constraint models for hypersonic vehicles such as heat flow, dynamic pressure, overload, and no-fly zones;

[0017] S5. According to the above model and constraints, the reentry process state variables and control variables are discretized based on the Radau pseudo-spectral method principle;

[0018] S6. For the constructed optimization problem, the fourth-order Ronge-Kutta method is used to integrate and solve the initial state variables and the estimated values ​​of the control variables;

[0019] S7. Based on discrete states and control variables, a nonlinear programming problem is established with minimum heat flow and trajectory oscillation as optimization objectives. The discrete states are adaptively corrected according to the accuracy requirements during the solution process.

[0020] Furthermore, the equation established in step S1 is:

[0021]

[0022]

[0023]

[0024]

[0025]

[0026]

[0027]

[0028] Among them, r is the distance from the center of the earth, λ is the longitude, φ is the latitude, S e represents the range angle, V is the flight speed, σ is the roll angle, θ is the ballistic inclination angle, and ψ is the heading angle, which is the angle between the local horizontal projection of the velocity vector and the true north direction, and clockwise rotation is positive; L, D represent the aerodynamic lift and drag respectively, and m0, g represent the mass of the aircraft and the earth's gravity acceleration at the height of the aircraft:

[0029]

[0030]

[0031] Where R0 represents the radius of the earth, S ref ,m,ρ represents the reference area, mass and atmospheric density of the aircraft, C L ,C D represents the lift coefficient and the drag coefficient.

[0032] Furthermore, the equation established in step S2 is:

[0033] ρ=ρ0e -h / 7200 (10)

[0034]

[0035] Where, ρ0 = 1.225 kg / m 3 is the standard atmospheric density, g0=9.81m / s 2 is the standard earth gravity acceleration, h is the aircraft's flight altitude, R0, r is the earth's radius and the distance from the aircraft to the earth's center.

[0036] Furthermore, in step S3, the method for analyzing the aerodynamic model of the hypersonic vehicle based on S1 and S2 is:

[0037] Assuming that the reentry glide missile torque is in a balanced state and there is no side force throughout the process, the aerodynamic force and aerodynamic torque of the aircraft during the reentry process are mainly determined by the following parameters:

[0038]

[0039] Among them, F and M represent the aerodynamic force and aerodynamic torque of the aircraft respectively, Re represents the Reynolds number to measure the magnitude of friction; Ma is the Mach number, α, β, Indicates the airflow angle and airflow angular velocity of the current gliding missile; Characterize the angular velocity of the gliding missile body; Characterizes the current rudder deflection angle of the gliding missile; wing reference area S ref , current dynamic pressure q, reference characteristic length l, after normalization processing, the function expression of aerodynamic parameters is obtained; shape, scale, and power represent the shape characteristics, size, and engine characteristics of the aircraft respectively;

[0040] Current gliding missile lift coefficient:

[0041]

[0042] Current glide missile drag coefficient:

[0043]

[0044] The above two physical quantities are determined by the following parameters:

[0045] C L |C D =f{M,Re,α,β,shape,poweron / off} (15)

[0046] After reasonable simplification, the final form is:

[0047] C L |C D =f{M,α,β} (16)

[0048] Under hypersonic flight conditions, the following aerodynamic coefficient fitting results are obtained by combining the aerodynamic data table:

[0049] C L =C L00 +C L10 α+C L01 M+C L20 α2 +C L11 αM+C L02 M 2 +C L21 α 2 M+C L12 αM 2 +C L03 M 3 +C L22 α 2 M 2 +C L13 αM 3 +C L04 M 4

[0050] C D =C D00 +C D10 α+C D01 M+C D20 α 2 +C D11 αM+C D02 M 2 +C D21 α 2 M+C D12 αM 2 +C D03 M 3 +C D22 α 2 M 2 +C D13 αM 3 +C D04 M 4 (17)

[0051] Where C Lij represents the lift polynomial fitting coefficient, C Dij Represents the resistance polynomial fitting coefficient.

[0052] Furthermore, the method for establishing the process constraint models of the hypersonic vehicle heat flow, dynamic pressure, overload, and no-fly zone in step S4 is:

[0053] The reentry process of a hypersonic vehicle is mainly limited by heat flow, dynamic pressure and overload constraints.

[0054]

[0055]

[0056]

[0057] In the formula, is the heat flux density at the stagnation point; is the maximum heat flux density allowed; ρ, ρ0 represent the atmospheric density corresponding to the actual height and the standard atmospheric density at sea level, respectively; C1, R d are the coefficients related to the aircraft respectively; Respectively represent the current flight speed and normalized speed; q max is the maximum allowable dynamic pressure, n is the total overload constraint; max is the maximum total overload allowed;

[0058] The no-fly zone is regarded as an infinitely high cylinder. The aircraft trajectory can only avoid it from the left and right sides, without considering avoiding it from the top or bottom. Assuming λ m ,φ m ,λ z ,φ z ,R z ,R e Respectively represent the current longitude and latitude of the aircraft, the longitude and latitude of the center of the no-fly zone, the radius of the no-fly zone and the radius of the earth, then the path constraints shown in the formula are satisfied:

[0059]

[0060] Furthermore, the method of step S5 is:

[0061] Through the Radau pseudo-spectral method, the control variables and state variables are discretized at a series of Legendre-Gauss-Radau points, the control variables and state variables are approximated by interpolation polynomials, and the dynamic differential equations are approximated by the derivation of the state variables, thereby transforming the continuous optimal control problem into a series of NLP problems subject to algebraic constraints, and then using the corresponding numerical method to solve the problem; the main steps are as follows:

[0062] (1) Time domain transformation;

[0063] (2) Based on time domain transformation, polynomial approximation of state variables and control variables;

[0064] (3) Transform the constraints of differential equations based on time domain transformation;

[0065] (4) Transform terminal constraints based on time domain transformation;

[0066] (5) Based on the time domain transformation, the performance index function is approximated.

[0067] Furthermore, the method of step S6 is:

[0068] The Ronge-Kutta integral formula is as follows:

[0069]

[0070] Where yn ,y n+1 They represent the current and next moment values ​​of the dependent variable of the function relationship f, respectively, n ,h d Respectively represent the value of the independent variable at the current moment and the update step size;

[0071] By presetting the roll angle control amount and the predetermined angle of attack profile, as well as the heading angle deviation corridor guidance, the estimated values ​​of the state variables and control variables that can roughly complete the reentry mission are obtained, providing initial value guesses for solving the subsequent nonlinear programming problem.

[0072] Furthermore, the method of step S7 is:

[0073] After discretizing the state variables, control variables and constraints, the NLP problem is obtained. The objective function is set to the minimum heat flow and trajectory oscillation, and the optimization problem is solved. The adaptive correction of the pseudo-spectral orthogonal nodes is considered during the solution process. The main steps are as follows:

[0074] (1) Approximation accuracy calculation;

[0075] (2) Based on the approximation accuracy calculation, design a subdivided grid or add a point matching criterion;

[0076] (3) Determine the number of points according to the evaluation results of the discrimination criteria;

[0077] (4) Based on the evaluation results of the discrimination criteria, the number and location of the divided time zones are determined.

[0078] Beneficial effects of the present invention:

[0079] The present invention discloses a hypersonic aircraft re-entry no-fly zone avoidance guidance method based on an adaptive pseudo-spectral method. Aiming at the high uncertainty, high nonlinearity, high complexity, high constraint environment of the hypersonic aircraft re-entry process, as well as the minimum heat flux requirement of re-entry and the requirement of avoiding the no-fly zone, the adaptive pseudo-spectral method is used to discrete state variables and control variables, and on this basis, the original complex trajectory planning problem is transformed into a nonlinear planning problem that can be solved, so as to realize the multi-constraint requirements of the re-entry trajectory planning and the requirement of avoiding the no-fly zone, and finally obtain a smaller aerodynamic heating effect. BRIEF DESCRIPTION OF THE DRAWINGS

[0080] Figure 1 The figure is a flow chart of the method of the present invention.

[0081] Figure 2 Flowchart for trajectory planning.

[0082] Figure 3 This is the coordinate system transformation diagram.

[0083] Figure 4This is the lift coefficient fitting effect diagram.

[0084] Figure 5 This is the drag coefficient fitting effect diagram.

[0085] Figure 6 This is the lift-to-drag ratio fitting effect diagram.

[0086] Figure 7 It is the height and speed change curve.

[0087] Figure 8 It is the curve of ballistic inclination and ballistic deviation.

[0088] Fig. 9 It is the curve of attack angle and roll angle change.

[0089] Fig.10 Schematic diagram of constraint satisfaction.

[0090] Fig.11 The attached figure is a schematic diagram of no-fly zone avoidance 1.

[0091] Fig.12 The attached figure is a schematic diagram of no-fly zone avoidance 2. DETAILED DESCRIPTION

[0092] In order to make the purpose, technical solution and advantages of the embodiments of the present invention clearer, the technical solution in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0093] refer to Figure 1 The present invention discloses a hypersonic vehicle no-fly zone avoidance guidance method based on pseudo-spectral method, comprising the following steps:

[0094] S1. Establish the unpowered reentry guidance dynamic equations of hypersonic vehicles;

[0095] S2. Establish the ambient atmospheric density and gravity acceleration model;

[0096] S3. Analyze the aerodynamic model of hypersonic vehicle based on S1 and S2;

[0097] S4. Establish process constraint models for hypersonic vehicles such as heat flow, dynamic pressure, overload, and no-fly zones;

[0098] S5. According to the above model and constraints, the reentry process state variables and control variables are discretized based on the Radau pseudo-spectral method principle;

[0099] S6. For the constructed optimization problem, the fourth-order Ronge-Kutta method is used to integrate and solve the initial state variables and the estimated values ​​of the control variables;

[0100] S7. Based on discrete states and control variables, a nonlinear programming problem is established with minimum heat flow and trajectory oscillation as optimization objectives. The discrete states are adaptively corrected according to the accuracy requirements during the solution process.

[0101] Figure 2 First, the entire reentry trajectory planning process is demonstrated, including the initial value generation, pseudo-spectral method discrete variables and solution steps.

[0102] The method of step S1 is as follows:

[0103] according to Figure 3 The relevant angles defined by the coordinate system transformation relationship are used to establish the motion equation of the hypersonic vehicle reentry process as follows:

[0104]

[0105]

[0106]

[0107]

[0108]

[0109]

[0110]

[0111] Among them, r is the distance from the center of the earth, λ is the longitude, φ is the latitude, S e represents the range angle, V is the flight speed, σ is the roll angle, θ is the ballistic inclination angle, and ψ is the heading angle, which is the angle between the local horizontal projection of the velocity vector and the true north direction, and clockwise rotation is positive; L, D represent the aerodynamic lift and drag respectively, and m0, g represent the mass of the aircraft and the earth's gravity acceleration at the height of the aircraft:

[0112]

[0113]

[0114] Where R0 represents the radius of the earth, S ref ,m,ρ represents the reference area, mass and atmospheric density of the aircraft, C L ,C D represents the lift coefficient and the drag coefficient.

[0115] The equation established in step S2 is:

[0116] ρ=ρ0e -h / 7200 (10)

[0117]

[0118] Where, ρ0 = 1.225 kg / m 3 is the standard atmospheric density, g0=9.81m / s 2 is the standard earth gravity acceleration, h is the aircraft's flight altitude, R0, r is the earth's radius and the distance from the aircraft to the earth's center.

[0119] The method for analyzing the aerodynamic model of the hypersonic vehicle based on S1 and S2 in step S3 is:

[0120] Assuming that the reentry glide missile torque is in a balanced state and there is no side force throughout the process, the aerodynamic force and aerodynamic torque of the aircraft during the reentry process are mainly determined by the following parameters:

[0121]

[0122] Among them, F and M represent the aerodynamic force and aerodynamic torque of the aircraft respectively, Re represents the Reynolds number to measure the magnitude of friction; Ma is the Mach number, α, β, Indicates the airflow angle and airflow angular velocity of the current gliding missile; Characterize the angular velocity of the gliding missile body; Characterizes the current rudder deflection angle of the gliding missile; wing reference area S ref , current dynamic pressure q, reference characteristic length l, after normalization processing, the function expression of aerodynamic parameters is obtained; shape, scale, and power represent the shape characteristics, size, and engine characteristics of the aircraft respectively;

[0123] Current gliding missile lift coefficient:

[0124]

[0125] Current glide missile drag coefficient:

[0126]

[0127] The above two physical quantities are determined by the following parameters:

[0128] C L |C D =f{M,Re,α,β,shape,poweron / off} (15)

[0129] After reasonable simplification, the final form is:

[0130] CL |C D =f{M,α,β} (16)

[0131] Under hypersonic flight conditions, the following aerodynamic coefficient fitting results are obtained by combining the aerodynamic data table:

[0132] C L =C L00 +C L10 α+C L01 M+C L20 α 2 +C L11 αM+C L02 M 2 +C L21 α 2 M+C L12 αM 2 +C L03 M 3 +C L22 α 2 M 2 +C L13 αM 3 +C L04 M 4

[0133] C D =C D00 +C D10 α+C D01 M+C D20 α 2 +C D11 αM+C D02 M 2 +C D21 α 2 M+C D12 αM 2 +C D03 M 3 +C D22 α 2 M 2 +C D13 αM 3 +C D04 M 4 (17)

[0134] Where C Lij represents the lift polynomial fitting coefficient, C Dij Represents the resistance polynomial fitting coefficient.

[0135] like Figure 4 , Figure 5 , Figure 6 As shown, by substituting the relevant fitting coefficients of the aircraft, the fitting results of lift, drag and lift-to-drag ratio corresponding to the aircraft are obtained through simulation.

[0136] The method for establishing the process constraint models of the hypersonic vehicle heat flow, dynamic pressure, overload, and no-fly zone in step S4 is:

[0137] The reentry process of a hypersonic vehicle is mainly limited by heat flow, dynamic pressure and overload constraints.

[0138]

[0139]

[0140]

[0141] In the formula, is the heat flux density at the stagnation point; is the maximum heat flux density allowed; ρ, ρ0 represent the atmospheric density corresponding to the actual height and the standard atmospheric density at sea level, respectively; C1, R d are the coefficients related to the aircraft respectively; Respectively represent the current flight speed and normalized speed; q max is the maximum allowable dynamic pressure, n is the total overload constraint; max is the maximum total overload allowed;

[0142] A no-fly zone is an area where reentry vehicles are affected by radar detection, electromagnetic interference, interception, terrain, geopolitical factors, etc. The reentry trajectory should try to avoid the no-fly zone. In order to facilitate design analysis, the no-fly zone is regarded as an infinitely high cylinder. The aircraft trajectory can only avoid it from the left and right sides, without considering avoiding it from the top or bottom; assuming λ m ,φ m ,λ z ,φ z ,R z ,R e Respectively represent the current longitude and latitude of the aircraft, the longitude and latitude of the center of the no-fly zone, the radius of the no-fly zone and the radius of the earth, then the path constraints shown in the formula are satisfied:

[0143]

[0144] The method of step S5 is:

[0145] Through the Radau pseudo-spectral method, the control variables and state variables are discretized at a series of Legendre-Gauss-Radau points, the control variables and state variables are approximated by interpolation polynomials, and the dynamic differential equations are approximated by the derivation of the state variables, thereby transforming the continuous optimal control problem into a series of NLP problems subject to algebraic constraints, and then using the corresponding numerical method to solve the problem; the main steps are as follows:

[0146] (1) Time domain transformation;

[0147] (2) Based on time domain transformation, polynomial approximation of state variables and control variables;

[0148] (3) Transform the constraints of differential equations based on time domain transformation;

[0149] (4) Transform terminal constraints based on time domain transformation;

[0150] (5) Based on the time domain transformation, the performance index function is approximated.

[0151] The method of step S5 specifically includes:

[0152] For the traditional optimal control problem, define the state variable x(·)∈R n , control variable u(·)∈R m , static parameter p∈R q Initial test time t0, terminal time t f (free or fixed), minimizing a general Bolza-type performance criterion:

[0153]

[0154] Satisfy the differential equation constraints:

[0155]

[0156] Path constraints:

[0157] C(x(t),u(t),t;p)≤0 (25)

[0158] Boundary constraints:

[0159]

[0160] Where p=[p1,p2,...,p q ] T ,x=[x1,x2,...,x n ] T ,u=[u1,u2,...,u m ] T .

[0161] The idea of ​​RPM to solve the optimal control problem is to discretize the control variables and state variables at a series of Legendre-Gauss-Radau (LGR) points, approximate the control variables and state variables by interpolation polynomials, and approximate the dynamic differential equations by taking the derivative of the state variables, thereby transforming the continuous optimal control problem into a series of NLP problems subject to algebraic constraints, and then using the corresponding numerical method to solve the problem. The steps are as follows:

[0162] (1) Time domain transformation

[0163] The Radau pseudospectral method requires the time interval of the trajectory optimization problem to be changed from the interval [t0,t f ] to [-1,1], the transformation is as follows:

[0164]

[0165] Where τ is the current time after transformation, t,t0,t f Indicates the current time, initial time and terminal time;

[0166] (2) Polynomial approximation of state variables and control variables

[0167] The collocation points of the Radau pseudospectral method are LGR points, which are values ​​in the interval (-1,1]. The k-th order LGR point is the polynomial P k (ζ)+P k-1 The root of (ζ), where P k (ζ) is a Legendre polynomial of order k, defined as follows:

[0168]

[0169] In order to make the collocation points cover [-1,1], the nodes of Radau pseudo-spectral method are selected as collocation points and initial time nodes ζ0 = -1. That is, if the number of discrete nodes of the optimization problem is N, the LGR collocation with K = N-1 is used, and the state variables are discretized in combination with the initial nodes as follows:

[0170]

[0171] τ i is the discretized point, N is the number of discrete points, x and X represent the actual state quantity and the discrete approximate state quantity respectively.

[0172] Among them are Lagrange interpolation basis functions

[0173]

[0174] And make the approximate state at the node equal to the actual state, there is X i =X(τi )=x(τ i ).

[0175] The derivative of the control variable is not used in the Radau pseudospectral method, so u(τ k )=U k , (k=1,...N-1), LGR nodes are still used, and Lagrange interpolation basis functions are used to approximate the control variables as u, U represents the actual control quantity and the approximate control quantity:

[0176]

[0177] (3) Differential equation constraint transformation

[0178] Differentiate the equation, and we have

[0179]

[0180] The differential matrix satisfies:

[0181]

[0182] Thus, the dynamic differential equation constraints are transformed into algebraic constraints:

[0183]

[0184] Where k=1,...N-1.

[0185] (4) Terminal constraints

[0186] The terminal state should satisfy the kinetic equation constraints as follows:

[0187]

[0188] The terminal constraints can be discretized and approximated by using Gauss integral, and we can get:

[0189]

[0190] Where:

[0191]

[0192] τ f Indicates the terminal time after time domain conversion;

[0193] (5) Performance index function approximation

[0194] By using Gauss integral to approximate the integral term of the Bolza-type performance index function, the approximate performance index function in the Radau pseudospectral method is obtained as follows:

[0195]

[0196] Where Φ() and g() represent the terminal index function and the integral index function.

[0197] Through the above mathematical transformation, the original continuous optimal control problem is transformed into a general nonlinear programming problem.

[0198] Solve for the state X at the discrete node i (i=0,1...,N-1) and control variable U k (k=1,2,...,N-1), and the initial terminal time t0,t f , minimize the performance index and satisfy the boundary conditions of the original optimal control problem.

[0199] The method of step S6 is:

[0200] The Ronge-Kutta integral formula is as follows:

[0201]

[0202] Where y n ,y n+1 They represent the current and next moment values ​​of the dependent variable of the function relationship f, respectively, n ,h d Respectively represent the value of the independent variable at the current moment and the update step size;

[0203] By presetting the roll angle control amount and the predetermined angle of attack profile, as well as the heading angle deviation corridor guidance, the estimated values ​​of the state variables and control variables that can roughly complete the reentry mission are obtained, providing initial value guesses for solving the subsequent nonlinear programming problem.

[0204] Step S7 specifically includes:

[0205] After discretizing the state variables, control variables and constraints using the Radau pseudospectral method, the NLP problem is obtained. The objective function is set to the minimum heat flow and trajectory oscillation, and the optimization problem is solved. The longitudinal trajectory of the aircraft reentry process is often a "Qian Xuesen" trajectory that jumps forward, with strong maneuverability and penetration capabilities, which also leads to a relatively oscillating reentry trajectory. The Radau pseudospectral method described above can be divided into a global pseudospectral method (p-RPM) and a local pseudospectral method (h-RPM). The global RPM has a faster convergence speed for smooth trajectory optimization, but if the state is more oscillating, it is difficult to ensure the accuracy of the optimization process; the local RPM effectively improves the accuracy by dividing the time domain and fitting each sub-time domain with a low-order polynomial, but the convergence speed is slow. The hp adaptive pseudospectral method developed in recent years can effectively combine the advantages of the two, reasonably divide the entire time domain into smooth sub-time domains, and is suitable for the reentry problem of aircraft.

[0206] The main steps of hp adaptive pseudo-spectrum are as follows:

[0207] (1) Approximation accuracy calculation:

[0208] Compute the first-order differential residuals of the state variables and the process constraint residuals at each collocation point in each time subinterval:

[0209]

[0210]

[0211] Where i, j represent the points in the interval and the labels from the interval, respectively, f, C represent the dynamic equations and process constraints;

[0212] From the above, we can see that the maximum residual in this interval is:

[0213]

[0214] Set the maximum allowable deviation ε, if e j <ε, indicating that the state and control quantity in interval j satisfy the differential equation constraints and process constraints, and the iterative calculation of this interval is stopped. Otherwise, the optimization accuracy should be further improved by subdividing the time region or increasing the number of points in this region. When all the self-intervals meet the accuracy requirements, a feasible solution is obtained.

[0215] (2) Subdividing the grid or adding points to the judgment criteria

[0216] The subdivision time (h) or the increase in the number of distribution points (p) is determined by the curvature of the sub-time domain interval. If the curvature is greater than the set value, it means that the interval is more oscillating and needs to be further divided into time areas. Otherwise, the number of distribution points is increased.

[0217] The relative curvature of the j-th sub-time region is defined as:

[0218]

[0219] in They represent the maximum curvature and average curvature of the matching points in the interval respectively. The curvature of the matching points is calculated as follows:

[0220]

[0221] (3) Determination of the number of points

[0222] If the interval is relatively smooth, increase the number of matching points, which is determined by the following formula:

[0223]

[0224] In the formula, ceil means rounding up, and A is an adjustable parameter greater than 0. It can be seen from the above formula that the larger the residual of the subinterval is, the more points need to be added, and vice versa, the fewer points need to be added.

[0225] (4) Determination of the number and location of time zones

[0226] The number of divided time zones is determined by the following formula:

[0227]

[0228] B is an adjustable parameter; the above formula shows that the larger the residual, the more oscillating the trajectory is, and more sub-time regions need to be divided. The position of the divided time region is determined according to the average curvature. The larger the curvature, the more oscillating the trajectory is, and the time period should be narrower to ensure the optimization accuracy. Define the curvature density coefficient constant c of the jth time region j satisfy:

[0229]

[0230] Then the i-th subdivision point in the j-th time zone should be arranged at the distribution point The place makes:

[0231]

[0232] Based on the above adjustments, the trajectory planning problem after discretization by adaptive pseudo-spectral method is started, such as Figure 7 As shown, the changes in altitude and speed during the reentry process are shown, indicating that the altitude and speed constraints of the terminal can be met; Figure 8 As shown in the figure, the changes of ballistic inclination and ballistic deviation during the reentry process are shown; Fig. 9 As shown in Figure 1, the changes in the angle of attack and the roll angle during the reentry process are shown; Fig.10As shown in the figure, the process constraints during the reentry process are satisfied. The curves in the figure show that the process constraints can be satisfied. Fig.11 and Fig.12 As shown, it indicates that the reentry process can avoid one or more no-fly zones and meet the reentry trajectory requirements.

[0233] The present invention aims at the high uncertainty, high nonlinearity, high complexity, high constraint environment of the reentry process of a hypersonic aircraft, as well as the minimum heat flux requirement for reentry and the requirement for avoiding no-fly zones. The present invention discretizes state variables and control variables by an adaptive pseudo-spectral method, and on this basis transforms the original complex trajectory planning problem into a solvable nonlinear planning problem, so as to realize the multi-constraint requirements of reentry trajectory planning and the requirement for avoiding no-fly zones, and ultimately obtain a smaller aerodynamic heating effect.

[0234] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit it. Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art should understand that they can still modify the technical solutions described in the aforementioned embodiments, or make equivalent replacements for some of the technical features therein. However, these modifications or replacements do not deviate the essence of the corresponding technical solutions from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A hypersonic vehicle no-fly zone avoidance guidance method based on pseudospectral method, characterized in that: The steps include: S1. Establish the unpowered reentry guidance dynamic equations of hypersonic vehicles; S2. Establish the ambient atmospheric density and gravity acceleration model; S3. Analyze the aerodynamic model of hypersonic vehicle based on S1 and S2; S4. Establish the hypersonic vehicle thermal flow, dynamic pressure, overload, and no-fly zone process constraint model; S5. According to the process constraint model and constraints, the reentry process state variables and control variables are discretized based on the Radau pseudo-spectral method principle; S6. For the constructed optimization problem, the fourth-order Ronge-Kutta method is used to integrate and solve the initial state variables and the estimated values ​​of the control variables; S7. Based on discrete states and control variables, a nonlinear programming problem is established with minimum heat flow and trajectory oscillation as optimization objectives. The discrete states are adaptively corrected according to the accuracy requirements during the solution process.

2. The method for avoiding a no-fly zone for a hypersonic vehicle based on pseudospectral method according to claim 1, characterized in that: The equation established in step S1 is: Among them, r is the distance from the aircraft to the center of the earth, λ is the longitude, φ is the latitude, S e represents the range angle, V is the current flight speed, σ is the roll angle, and θ is the ballistic inclination angle; ψ is the heading angle, which is the angle between the local horizontal projection of the velocity vector and the true north direction, and clockwise rotation is positive; L, D represent the aerodynamic lift and drag respectively, and m0, g represent the mass of the aircraft and the earth's gravity acceleration at the height of the aircraft: Where R0 represents the radius of the earth, S ref ,m,ρ represents the reference area, mass and actual altitude of the aircraft corresponding to the atmospheric density, C L ,C D Indicates the current lift coefficient and drag coefficient of the gliding missile.

3. The method for avoiding a no-fly zone for a hypersonic vehicle based on pseudospectral method according to claim 2, characterized in that: The equation established in step S2 is: p=ρ0e -h / 7200 (10) Where, ρ0 = 1.225 kg / m 3 is the standard atmospheric density at sea level, g0 = 9.81 m / s 2 is the standard earth gravity acceleration, h is the aircraft's flight altitude, R0, r is the earth's radius and the distance from the aircraft to the earth's center.

4. The method for avoiding a no-fly zone for a hypersonic vehicle based on pseudospectral method according to claim 3, characterized in that: The method for analyzing the aerodynamic model of the hypersonic vehicle based on S1 and S2 in step S3 is: Assuming that the reentry glide missile torque is in a balanced state and there is no side force throughout the process, the aerodynamic force and aerodynamic torque of the aircraft during the reentry process are mainly determined by the following parameters: Among them, F and M represent the aerodynamic force and aerodynamic torque of the aircraft respectively, Re represents the Reynolds number to measure the magnitude of friction; Ma is the Mach number, α, β, Indicates the airflow angle and airflow angular velocity of the current gliding missile; Characterize the angular velocity of the gliding missile body; Characterizes the current rudder deflection angle of the gliding missile; wing reference area S ref , current dynamic pressure q, reference characteristic length l, after normalization, the function expression of aerodynamic parameters is obtained; shape, scale, and power represent the shape characteristics, size, and engine characteristics of the aircraft respectively; Current gliding missile lift coefficient: Current glide missile drag coefficient: The above two physical quantities are determined by the following parameters: C L |C D =f{M,Re,α,β,shape,poweron / off} (15) After reasonable simplification, the final form is: C L |C D =f{M,α,β} (16) Under hypersonic flight conditions, the following aerodynamic coefficient fitting results are obtained by combining the aerodynamic data table: C L =C L00 +C L10 α+C L01 M+C L20 α 2 +C L11 αM+C L02 M 2 +C L21 α 2 M+C L12 αM 2 +C L03 M 3 +C L22 α 2 M 2 +C L13 αM 3 +C L04 M 4 C D =C D00 +C D10 α+C D01 M+C D20 α 2 +C D11 αM+C D02 M 2 +C D21 α 2 M+C D12 αM 2 +C D03 M 3 +C D22 α 2 M 2 +C D13 αM 3 +C D04 M 4 (17) Where C Lij represents the lift polynomial fitting coefficient, C Dij Represents the resistance polynomial fitting coefficient.

5. The method for avoiding a no-fly zone for a hypersonic vehicle based on pseudospectral method according to claim 4, characterized in that: The method for establishing the hypersonic vehicle heat flow, dynamic pressure, overload, and no-fly zone process constraint model in step S4 is: The reentry process of a hypersonic vehicle is mainly limited by heat flow, dynamic pressure and overload constraints. In the formula, is the heat flux density at the stagnation point; is the maximum allowable heat flux density; ρ, ρ0 represent the atmospheric density corresponding to the actual height and the standard atmospheric density at sea level, respectively; C1, R d are the coefficients related to the aircraft respectively; Respectively represent the current flight speed and normalized speed; q max is the maximum allowable dynamic pressure, n is the total overload constraint; max is the maximum total overload allowed; The no-fly zone is regarded as an infinitely high cylinder. The aircraft trajectory can only avoid it from the left and right sides, without considering avoiding it from the top or bottom. Assuming λ m ,φ m ,λ z ,φ z ,R z ,R e Respectively represent the current longitude and latitude of the aircraft, the longitude and latitude of the center of the no-fly zone, the radius of the no-fly zone and the radius of the earth, then the path constraints shown in the formula are satisfied:

6. The method for avoiding a no-fly zone for a hypersonic vehicle based on pseudospectral method according to claim 5, characterized in that: The method of step S5 is: Through the Radau pseudo-spectral method, the control variables and state variables are discretized at a series of Legendre-Gauss-Radau points, the control variables and state variables are approximated by interpolation polynomials, and the dynamic differential equations are approximated by the derivation of the state variables, thereby transforming the continuous optimal control problem into a series of NLP problems subject to algebraic constraints, and then using the corresponding numerical method to solve the problem; the main steps are as follows: (1) Time domain transformation; (2) Based on time domain transformation, polynomial approximation of state variables and control variables; (3) Transform the constraints of differential equations based on time domain transformation; (4) Transform the terminal constraints based on time domain transformation; (5) Based on the time domain transformation, the performance index function is approximated.

7. The method for avoiding a no-fly zone for a hypersonic vehicle based on pseudospectral method according to claim 6, characterized in that: The method of step S6 is: The Ronge-Kutta integral formula is as follows: Where y n ,y n+1 They represent the current and next moment values ​​of the dependent variable of the function relationship f, respectively, n ,h d Respectively represent the value of the independent variable at the current moment and the update step size; By presetting the roll angle control amount and the predetermined angle of attack profile, as well as the heading angle deviation corridor guidance, the estimated values ​​of the state variables and control variables that can roughly complete the reentry mission are obtained, providing initial value guesses for solving the subsequent nonlinear programming problem.

8. The method for avoiding a no-fly zone for a hypersonic vehicle based on pseudospectral method according to claim 7, characterized in that: The method of step S7 is: After discretizing the state variables, control variables and constraints, we get the NLP problem. We set the objective function to the minimum heat flow and trajectory oscillation, and start solving the optimization problem. During the solution process, we consider the adaptive correction of the pseudo-spectral orthogonal nodes. The main steps are as follows: (1) Approximation accuracy calculation; (2) Based on the approximation accuracy calculation, design a subdivided grid or add a point matching criterion; (3) Determine the number of points according to the evaluation results of the discrimination criteria; (4) Based on the evaluation results of the discrimination criteria, the number and location of the divided time zones are determined.

Images