Neural Network-Based Trajectory Optimization Method for Aircraft Reentry

Abstract

The present invention discloses a method for optimizing the trajectory of an aircraft during reentry based on a neural network, which mainly solves the problems of poor real-time performance and adaptability in the prior art. The implementation scheme is as follows: establishing the continuous optimal control problem of aircraft reentry in the semi-speed coordinate system; transforming the continuous optimal control problem of aircraft reentry into a sequential convex optimal control problem; transforming the sequential convex optimal control problem into a sequential second-order cone programming problem; solving the sequential second-order cone programming problem; sampling from the solution results to obtain a state quantity data set and a control quantity data set; constructing a neural network and a loss function; using part of the state quantity data set as a training data set to perform offline training on the neural network until the loss function converges to a minimum value, obtaining a trained trajectory network; and using the trajectory network to online obtain the trajectory optimization result of the aircraft during reentry. The present invention has strong adaptability and good real-time performance, reduces the influence of parameter changes on the aircraft, and can be used for rocket recovery.

Description

Technical Field

[0001] The present invention belongs to the technical field of guidance and control, and particularly relates to a method for optimizing the trajectory of a vehicle during the reentry phase, which can be used for rocket recovery. Background Art

[0002] Hypersonic vehicles have important application backgrounds in fields such as aerospace return and long-range rapid strikes. However, due to the influence of the vehicle's own aerodynamic layout and the harsh flight environment, the dynamic model of hypersonic vehicles has the characteristics of high nonlinearity, strong coupling, and parameter uncertainty. At the same time, due to the complex flight constraints such as dynamic pressure, heat flux, and overload on the vehicle, the problem of optimizing the reentry trajectory is very complex. Therefore, the problem of optimizing the reentry flight trajectory has always been a key and difficult point in the research of the aerospace engineering field. With the development of aerospace technology, trajectory optimization with strong adaptability, fast response, and high reliability has become the pursuit goal.

[0003] Optimizing the trajectory of a vehicle during the reentry phase usually requires solving a continuous-time optimal control problem with one performance index and multiple constraints on state variables and control variables. Common performance indices include minimizing fuel consumption, minimizing the reentry time, minimizing the total heat load, or maximizing the range, speed, etc. The constraints include the equations of motion, necessary safety constraints such as dynamic pressure, heat flux, and overload, and terminal target state constraints.

[0004] There are many methods involved in the existing vehicle trajectory optimization problems, and the classification criteria are not unique. The indirect method and the direct method are the most commonly used classification criteria and also the two basic solution frameworks for trajectory optimization problems. Among them, the indirect method framework is based on the minimum principle, and the trajectory optimization problem of the entire reentry phase is transformed into a two-point boundary value problem to solve. The greatest advantage of the indirect method is the high accuracy of the solution. However, due to the high sensitivity of the two-point boundary value problem to the initial co-state and the strong nonlinear characteristics of the dynamic equations considering lift and drag in the atmospheric flight segment, the derivation process of the optimal control equation is complex. These defects limit the use of the indirect method. Although the direct method framework does not require the derivation of the optimal control equation and directly uses a nonlinear optimization algorithm to solve the optimal control problem, with relatively low sensitivity to the initial solution, its solution efficiency is low, and it is difficult to meet the real-time requirements for online applications. In addition, both of these frameworks are difficult to adapt to the influence of uncertain factors such as aerodynamic disturbances and thrust deviations during the flight process, which limits the further improvement of the overall performance of the vehicle. Summary of the Invention

[0005] The purpose of the present invention is to overcome the defects of the above-mentioned existing technologies, and propose a method for optimizing the trajectory of a vehicle during the reentry phase based on a neural network, so as to enhance the flight efficiency of the vehicle, improve the solution efficiency and accuracy, reduce the interference of various uncertain factors during the flight process, and improve the overall performance of the vehicle.

[0006] To achieve the above object, the technical solution of the present invention includes the following steps:

[0007] 1. A method for optimizing the trajectory of a vehicle during reentry based on a neural network, characterized by including the following steps:

[0008] (1) Describe the trajectory optimization of the vehicle during reentry as a continuous optimal control problem P0 composed of a mathematical model, boundary conditions, admissible control, performance index, and process constraints;

[0009] (2) Perform convexification processing on P0 by means of form transformation, slack variable, softening constraint, and successive linearization method to obtain a sequence of convex optimal control problems P1, and use the pseudospectral method to perform discrete parameterization processing on P1 to obtain a sequence of second-order cone programming problems P2;

[0010] (3) Solve the sequence of second-order cone programming problems P2 by the interior point method to obtain a nominal optimal reference trajectory including the vehicle state quantity curve and the vehicle control quantity curve;

[0011] (4) Establish a deviation model by offsetting the vehicle aerodynamic parameters, and perform offline solution on each set of aerodynamic parameters in the deviation model to obtain a non-nominal optimal reference trajectory including the vehicle state quantity curve and the vehicle control quantity curve under different parameter conditions;

[0012] (5) Sample the state quantity curves in the nominal optimal reference trajectory in step (3) and the non-nominal optimal reference trajectory in step (4) respectively to obtain a state quantity data set X including the four state variables of the vehicle's geocentric distance, longitude, latitude, and heading angle, and sample the control quantity curve in the optimal reference trajectory to obtain a control quantity data set Y1;

[0013] (6) Construct a neural network A composed of an input layer, two hidden layers, and an output layer cascaded in sequence, and set its loss function as: Loss(Y1,U1)=(Y1 - U1) 2 , where Y1 represents the control quantity data set, and U1 represents the output value of the neural network A, that is, the control quantity bank angle;

[0014] (7) Take out a part of the state quantity data set X in step (5) as a training trajectory set X1 and input it into the neural network A, and perform offline training on it using the BP algorithm. When the loss function of the neural network A converges to a minimum value, obtain a trained trajectory network A';

[0015] (8) Obtain the trajectory optimization result of the vehicle during reentry online:

[0016] During the flight of the aircraft, the on-board computer reads the real-time flight state variables measured by the aircraft navigation system, uses the real-time flight state variables as the input of the trained trajectory network A' for forward propagation, and obtains the real-time control variables;

[0017] Connect the state variables and control variables read each time into lines respectively to obtain the optimal trajectory of the aircraft during the reentry phase.

[0018] Compared with the prior art, the present invention has the following advantages:

[0019] 1) The present invention adopts a neural network model, transfers a large number of optimization calculations to the offline training process, is not affected by the complex non-linear characteristics of the aircraft reentry trajectory, not only reduces the online calculation amount, improves the solution efficiency and accuracy, but also has real-time performance.

[0020] 2) By establishing a deviation model, the present invention enables the trajectory planned by the neural network to adapt to parameter deviations within a large range, and reduces the interference of various uncertain factors during flight on the aircraft. Description of the Drawings

[0021] Figure 1 is the implementation flowchart of the present invention;

[0022] Figure 2 is the original non-convex set used in constructing the optimal trajectory of the aircraft in the present invention;

[0023] Figure 3 is the convex admissible control set obtained in convexifying the optimal trajectory of the aircraft in the present invention;

[0024] Figure 4 is the state variable envelope diagram of the present invention;

[0025] Figure 5 is the control variable envelope diagram of the present invention;

[0026] Figure 6 is the schematic diagram of optimal reference trajectory sampling of the present invention;

[0027] Figure 7 is the comparison diagram of the result fitted by the neural network in the present invention and the nominal optimal reference trajectory;

[0028] Figure 8 is the altitude curve diagram obtained by using the neural network in the present invention under nominal conditions;

[0029] Figure 9 is the longitude and latitude curve diagram obtained by using the neural network in the present invention under nominal conditions;

[0030] Figure 10 is the altitude curve obtained by using the neural network in the present invention under non-nominal parameter conditions;

[0031] Figure 11 This is the longitude and latitude curve obtained by using the neural network under non-nominal parameter conditions in the present invention. Detailed implementation manner

[0032] The following further elaborates on the embodiments and effects of the present invention with reference to the accompanying drawings.

[0033] Refer to Figure 1 , and the implementation steps of this example are as follows:

[0034] Step 1, establish the aircraft motion model.

[0035] The motion equation of the aircraft's center of mass is the basis for studying its motion characteristics. In this example, a reusable launch vehicle without power (RLV) is taken as the research object, and its re-entry center-of-mass motion equation is established in the semi-speed system.

[0036] The trajectory optimization problem usually does not consider the action of the control force, and since the additional Coriolis force level is small, its influence is ignored. The RLV re-entry process mainly relies on aerodynamic force and the earth's gravity to change the flight state. Assuming the earth is an ideal sphere, in the semi-speed coordinate system, the state differential equation f1 of the RLV is established as follows:

[0037]

[0038] Among them, r is the dimensionless geocentric distance of the RLV, is the dimensionless velocity of the RLV, Ω is the earth's rotation angular velocity, R0 is the earth's radius, γ is the flight path angle, ψ is the course angle, σ is the bank angle, L = 0.5R0ρV 2 SC L / m is the lift acceleration, D = 0.5R0ρV 2 SC D / m is the drag acceleration; e is the dimensionless energy of the RLV, C D and C L are the drag coefficient and the lift coefficient respectively. Both of these coefficients are functions of the angle of attack α and the Mach number Ma. S is the reference area of the RLV. ρ = ρ0exp(-(r - 1) / H) represents the atmospheric density, which is a function of the geocentric distance r. ρ0 is the atmospheric density at sea level, and H is the dimensionless scale height.

[0039] For the state differential equation, in this embodiment, the state variables are selected as s = (r, θ, φ, γ, ψ), and the control variable is the bank angle u = σ.

[0040] Step 2, establish the aircraft constraint conditions and performance indicators.

[0041] 2.1) Establish the constraint conditions of the aircraft:

[0042] 2.1.1) Establish the process constraint K:

[0043] When the reusable launch vehicle without power (RLV) flies in the atmosphere, a large amount of heat is generated due to the friction between the airframe and the atmosphere. Considering the structural strength of the RLV and the safe working conditions of the equipment, for flight safety, the dynamic pressure q, heat flux rate and normal overload n are constrained, denoted as K:

[0044]

[0045] where k Q is the heat flux rate calculation coefficient, g0 is the sea-level gravitational acceleration, q max , and n max are respectively the maximum dynamic pressure, maximum heat flux rate and maximum overload that the RLV can withstand;

[0046] 2.1.2) Establish the initial state constraint L and the terminal state constraint Φ1:

[0047] According to the flight mission of the RLV, the initial state information and terminal state information of the RLV can be determined, and then constraints are imposed on the initial state and terminal state of the RLV:

[0048] Initial state constraint L:

[0049]

[0050] where x represents the five state variables of the aircraft, represents the non-dimensional initial state of the aircraft, and e0 represents the initial non-dimensional energy of the RLV;

[0051] Terminal state constraint Φ1:

[0052]

[0053] where e f represents the terminal non-dimensional energy of the RLV, respectively represent the set values of the terminal non-dimensional geocentric distance, longitude, latitude, track angle and course angle of the RLV;

[0054] 2.1.3) Establish the admissible control C:

[0055] During the flight of the RLV, in order to ensure that the RLV will not perform large-angle flips, the amplitude of the control variable bank angle should also be restricted. Therefore, an admissible control constraint C is established to constrain the amplitude of the bank angle;

[0056] Admissible control C:

[0057] C: σ min ≤|σ|≤σ max

[0058] Among them, σ min and σ max respectively represent the lower and upper limits of the bank angle.

[0059] 2.2) Establish the performance index J

[0060] In the trajectory optimization of the vehicle reentry, the usual desired requirements for the vehicle include the shortest reentry time, the maximum terminal velocity, the longest flight range, and the minimum energy consumption. In this example, the performance index J is the shortest reentry time, that is:

[0061]

[0062] Step 3, construct the optimal control problem P0 of the RLV reentry trajectory optimization model.

[0063] According to the above state differential equation f1, process constraint K, initial constraint L, terminal constraint Φ1, admissible control C, and performance index J, the trajectory optimization model of the RLV reentry section is described as the optimal control problem P0 including these parameters:

[0064]

[0065] Step 4, transform the vehicle reentry continuous optimal control problem into a sequential convex optimal control problem P1.

[0066] 4.1) Form transformation of the process constraint K:

[0067] Redescribe the process constraint K in the form of a linear inequality about the state variable r to meet the requirements of a convex problem:

[0068]

[0069] The above formula is equivalent to

[0070]

[0071] Among them, l Q (e), l q (e), l n (e) respectively represent the functions of the heat flux rate, dynamic pressure, and overload of the vehicle with respect to the dimensionless geocentric distance r, denotes the maximum value among l Q (e), l q (e), l n (e);.

[0072] 4.2) Relaxation of the control constraint:

[0073] 4.2.1) Transform the control variable

[0074] When the bank angle σ is selected as the control variable, undesirable oscillations are likely to occur during the successive linearization of the differential equations. Therefore, it is necessary to change the control variable to eliminate such undesirable oscillations, as detailed below:

[0075] First, based on the non-convex nature of the centroid motion equation f1 with respect to both the state variable x and the control variable u = σ when the bank angle σ is selected as the control variable, the successive linearization method is used to convexify the equation f1. The solution at the k-th step in the successive solution process includes the state variable solution x (k) (e) and the control variable solution u (k) (e), denoted as {x (k) (e); u (k) (e)}.

[0076] Second, based on the solution obtained in the k-th iteration, the equation for the (k + 1)-th iteration is linearized as:

[0077] x′ = F x (x (k) , u (k) , e)x + F u (x (k) , u (k) , e)u + b(x (k) , u (k) , e)

[0078] where x′ represents the derivative of the five state variables with respect to the dimensionless energy e, represents the Jacobian matrix of the differential equation system f with respect to x, represents the Jacobian matrix of the differential equation system f with respect to u, and b(x (k) , u (k) , e) = f(x (k) , u (k) , e) - F x (x (k) , u (k) , e)x (k) - F u (x (k) , u (k) , e)u (k) represents the co-vector.

[0079] Then, based on the fact that the above coefficient matrices F x , F u and the co-vector b are related to u (k) , oscillations are likely to occur in the control variable u (k) during the iteration process. That is, during the iteration process, the oscillation of the control variable u (k) is used to control the control variable u (k +1), as the iteration progresses, the oscillation may be further amplified. These oscillations are not inherent in the solution of the problem and are manifestations of the instability of the numerical solution process. To eliminate the relationship between the linearized dynamic equation and u (k) Introduce the following control variable to replace the original control quantity: where u1 represents the cosine value of the tilt angle and u2 represents the sine value of the tilt angle. Separate the replaced control variable from the centroid motion equation, linearize the differential equation system based on the new control variable, and its Jacobian coefficient matrix is approximately a zero matrix, which can eliminate the adverse oscillation transmission.

[0080] 4.2.2) Establish new control variable constraints

[0081] After the control variable transformation, the control constraint is expressed as: cosσ max ≤u1≤cosσ min , since the control quantities u1 and u2 are not independent of each other, the following equation also needs to be satisfied:

[0082]

[0083] This equality constraint is a quadratic constraint, and the original admissible control set determined by the constraint is as Figure 2 shown. This original admissible control set is two arcs on the unit circle between u1 = ω l = cosσ max and u1 = ω h = cosσ min . The set composed of these two arcs is non-convex and needs to be convexified, that is, the constraint is relaxed, and the original constraint is relaxed to the following second-order cone inequality constraint:

[0084]

[0085] The set determined by relaxing the constraint expands the original admissible control set, making the original non-convex set become a convex set, as Figure 3 shown in the shaded area.

[0086] 4.3) Successive linearization of the performance index:

[0087] 4.3.1) Substitute the drag acceleration expression into the performance index and express the performance index as a function of the following state variables:

[0088]

[0089] where,, D = 0.5R0ρV 2 SC D / m is the drag acceleration, and the coefficient does not contain any state variables, where \(V\) is the dimensionless velocity of the RLV, \(m\) is the mass of the vehicle, \(\Omega\) is the angular velocity of the Earth's rotation, \(R_0\) is the radius of the Earth, \(e\) is the dimensionless energy of the RLV, and \(C\) D is the drag coefficient, which is a function of the angle of attack \(\alpha\) and the Mach number \(Ma\), \(S\) is the reference area of the RLV, \(\rho=\rho_0\exp(-(r - 1) / H)\) represents the atmospheric density, which is a function of the geocentric distance \(r\), \(\rho_0\) is the atmospheric density at sea level, and \(H\) is the dimensionless scale height.

[0090] 4.3.2) Linearization of performance index:

[0091] The atmospheric density \(\rho\) in the performance index formula is the only quantity that can be optimized, and the relevant state component is the geocentric distance \(r\). The atmospheric density approximately adopts an exponential model. The integrand of this performance index is a convex function of \(r\). However, to finally obtain a second-order cone programming problem, it needs to be transformed into a linear function of \(r\). Therefore, near the solution \(r\) (k) (e) of the previous iteration, a first-order Taylor expansion is performed on the integrand of this performance index:

[0092]

[0093] where \(\rho\) (k) represents the atmospheric density calculated using \(r\) (k) (e), represents the value of the derivative of the atmospheric density \(\rho\) with respect to \(r\) at \(r\) (k) (e), respectively represent the first-order term and the remainder term of the performance index with respect to the dimensionless geocentric distance \(r\) near the solution \(r\) (k) (e) of the previous iteration;

[0094] 4.3.3) Introducing a regularization term to ensure constraint invariance:

[0095] After linearizing the original performance index, a regularization term also needs to be introduced to ensure that the positive effect is produced after relaxing the constraints in step 4.2), that is, introducing

[0096] where \(\varepsilon\) ψ ≠0 is a constant. To ensure that the obtained solution is approximately the solution of the original problem, the magnitude of \(\varepsilon\) ψ needs to be set small enough.

[0097] 4.4) Relaxation of terminal constraints:

[0098] Although the terminal constraints are all linear equality constraints, in the initial optimization process, it may be difficult to satisfy the two hard constraints of longitude and latitude. Therefore, the following penalty terms are introduced in the performance index to replace these two terminal constraints:

[0099]

[0100] Among them, d(θ(e f ), φ(e f )) represents the penalty term, c θ > 0 and c φ > 0 are given constants,

[0101] Since this penalty term is non-linear, two slack variables and need to be introduced to transform it into a linear function and incorporate it into the performance index: And it is subject to linear inequality constraints: The performance index is updated to:

[0102]

[0103] Among them, c1 and c0 respectively represent the first-order term and the remainder term of the dimensionless geocentric distance r in the successive linearization process of the performance index, respectively represent the set values of the dimensionless geocentric distance, longitude, latitude, track angle, and heading angle of the RLV terminal.

[0104] 4.5) Successive linearization of the centroid motion equation:

[0105] According to the replaced control quantity, the centroid motion equation can be expressed as:

[0106]

[0107] Among them, f1(x, e) is the coefficient matrix of the terms without the control quantity in the original differential equation system, and B(x, e) is the coefficient matrix of the terms with the control quantity in the original differential equation system;

[0108] Using the successive linearization method to convexify the centroid motion equation, the solution at the k-th time in the successive solution process includes the state variable solution x (k) (e) and the control quantity solution u (k) (e), denoted as {x (k) (e); u (k) (e)};

[0109] According to the solution obtained in the k-th iteration, the equation is linearized in the (k + 1)-th iteration as:

[0110] x′ = F x (x (k) , u (k) , e)x + B(x (k) , e)u + b(x (k) , u (k) , e)

[0111] Among them, represents the Jacobian matrix;

[0112] b(x (k) , u (k) , e) = f1(x (k) , e) - F x (x (k) , u (k) , e)x (k) represents the cotangent vector;

[0113]

[0114] represents the partial derivative of B(x, e)u with respect to x, represents the partial derivative of f1(x, e) with respect to x;

[0115] L / D = C L / C D , both of which can be regarded as single-variable functions of energy e. Therefore, the only element in matrix B(x, e) related to the state variable x is cosγ, resulting in having a non-zero term. Since the flight path angle |γ| is very small during the reentry of the aircraft, that is, this non-zero term is very small, the linearized equation can be simplified as:

[0116] x' = F x (x (k) , e)x + B(x (k) , e)u + b(x (k) , e);

[0117] Furthermore, the terms related to the Earth's rotation, which are very small in magnitude, are separated from the equations of motion of the center of mass, and the equations of motion of the center of mass are updated to the following concise form:

[0118] x' = f1(x, e) + B(x, e)u = f0(x, e) + f Ω (x, e) + B(x, e)u

[0119] where, represents the terms in the differential equation that do not contain the control variable, and f Ω (x, e) is the term related to the Earth's angular rotation rate;

[0120] Still further, the updated equations of motion of the center of mass are linearly approximated successively, and the dimensionless energy e in the symbols is omitted. According to the solution obtained in the k-th iteration, the equation for the (k + 1)-th iteration is linearized as:

[0121] x' = A(x (k) )x + B(x (k) )u + b(x (k) )

[0122] where, b(x (k) ) = f0(x(k) ) - A(x (k) )x (k) +f Ω (x (k) ) represents the co - vector, A(x (k) ) and B(x (k) ) represent the Jacobian matrices of f0(x (k) ) with respect to x and u respectively, f0(x (k) ) represents the state differential equation obtained from the k - th iteration of f0(x, e), f Ω (x (k) ) represents the state differential equation obtained from the k - th iteration of f Ω (x, e);

[0123] To ensure the local validity of successive linearization, a trust - region constraint is imposed: |x - x (k) | ≤ δ, where δ is a given five - dimensional constant vector.

[0124] 4.6) According to the convexification treatment of the constraints and performance indices in steps 4.1) to 4.5), the trajectory optimization model of the RLV re - entry section is re - described as a sequential convex optimal control problem P1:

[0125]

[0126]

[0127] s.t. x′ = A(x (k) )x + B(x (k) )u + b(x (k) )

[0128] |x(e) - x (k) (e)| ≤ δ

[0129]

[0130]

[0131]

[0132]

[0133]

[0134] Step 5, transform the sequential convex optimal control problem P1 into a sequential second - order cone programming problem P2.

[0135] According to the principle of the Gauss pseudospectral method, the sequential convex optimal control problem P1 is discretized into a sequential convex programming problem P2, and its implementation steps are as follows:

[0136] 5.1) Set collocation points and discrete points:

[0137] Convert the energy independent variable e to the pseudo - energy independent variable E ∈ [-1, 1]:

[0138]

[0139] where e0, e f represent the initial non - dimensional energy and the terminal non - dimensional energy of the aircraft respectively;

[0140] Set N LG points (E1, E2, …, E m , …, E N ) as collocation points within E ∈ (-1, 1). The coordinates of each collocation point are the roots of the N - th order Legendre polynomial , where m = 1, 2, …, N represents each LG point;

[0141] Set discrete points, including the two end - points E0, E f of the pseudo - energy independent variable and all LG points inside the intervals, that is, a total of N + 2 discrete points E0, E1, …, E n , …, E N+1 , where n = 0, 1, …, N + 1 represents each discrete point, and E N+1 = E f represents the right - hand end - point of the pseudo - energy independent variable.

[0142] 5.2) Discretization of the differential equation of the center - of - mass motion:

[0143] 5.2.1) Use the Lagrange interpolation polynomial to represent the constraint of the differential equation of the center - of - mass motion for the pseudo - energy independent variable E ∈ [-1, 1):

[0144] Use N + 1 Lagrange interpolation polynomials as basis functions to approximate the state variables at each LG point , where i = 0, 1, …, N represents each interpolation polynomial;

[0145] Take the derivative of the approximate state variable x(E m ) with respect to the pseudo - energy independent variable E, that is where, represents the derivative of the basis function L i (E m ) at the LG point;

[0146] Use to replace x′ in the original differential equation of the center - of - mass motion, and transform the constraint of the differential equation of the center - of - mass motion into an algebraic equation constraint at the LG collocation points:

[0147]

[0148] Among them, represents the solution of the state variable at the k-th iteration at the m-th LG point, x m represents the solution of the state variable at the m-th LG point, u m represents the solution of the control quantity at the m-th LG point;

[0149] 5.2.2) Use Gauss integration to estimate the constraint of the centroid motion differential equation at the right endpoint of the pseudo-energy independent variable:

[0150] Since the pseudo-energy interval corresponding to the approximate expression of the state variable is [-1, 1), the terminal state variable x N+1 is not included. Therefore, Gauss integration is used to estimate x N+1 , that is, the constraint of the centroid motion differential equation that the terminal state x N+1 needs to satisfy:

[0151]

[0152] Among them, w m represents the LG integration weight.

[0153] 5.3) Set other discrete constraint conditions:

[0154] 5.3.1) Trust region constraint:

[0155]

[0156] Among them, represents the solution of the state variable at the k-th iteration at the n-th discrete point;

[0157] 5.3.2) Control constraint:

[0158]

[0159] Among them, (u1) m and (u2) m respectively represent the values of the control quantities u1 and u2 at the m-th LG point, ω l and ω h respectively represent the minimum and maximum values of the cosine of the bank angle;

[0160] 5.3.3) Process constraint:

[0161]

[0162] Among them, r m represents the dimensionless geocentric distance of the aircraft at the m-th LG point, represents the maximum values of the heat flux rate, dynamic pressure and overload of the aircraft at the m-th LG point;

[0163] 5.3.4) Boundary Conditions:

[0164]

[0165]

[0166]

[0167]

[0168]

[0169]

[0170] Among them, represents the initial value set for the aircraft state variables, r N+1 represents the terminal state of the dimensionless geocentric distance of the aircraft, θ N+1 represents the terminal state of the longitude of the aircraft, φ N+1 represents the terminal state of the latitude of the aircraft, γ N+1 represents the terminal state of the flight path angle of the aircraft, ψ N+1 represents the terminal state of the heading angle of the aircraft, respectively represent the terminal dimensionless geocentric distance, longitude, latitude, flight path angle and heading angle of the aircraft, and respectively represent the relaxation variables of longitude and latitude;

[0171] 5.3.5) Performance Index:

[0172] The performance index function is approximated by Gauss integration as:

[0173]

[0174] Among them, c θ , respectively represent the penalty coefficients of longitude and latitude, respectively represent the first-order term and the remainder term with respect to the dimensionless geocentric distance r in the successive linearization process of the performance index, ε ψ represents the heading angle regularization term coefficient, ψ m represents the value of the heading angle at the m-th LG point.

[0175] 5.4) According to the discretization processing of each constraint and performance index in Steps 5.1) to 5.3), the trajectory optimization model of the RLV reentry section is re-described as a sequential convex optimal control problem P2:

[0176]

[0177]

[0178]

[0179]

[0180]

[0181]

[0182]

[0183]

[0184]

[0185]

[0186] Among them, n = 0, 1, 2... N + 1 represents each discrete point, and m = 1, 2, 3... N represents each LG point.

[0187] Step 6: Solve the off-line second-order cone programming problem P2 of the re-entry section sequence of the non-powered reusable launch vehicle RLV.

[0188] The specific implementation of this step is as follows:

[0189] 6.1) Let the iteration count k = 0, and give the initial state quantity at the discrete point represents the solution of the state variable at the k-th iteration at the n-th discrete point, n = 0, 1, 2... N + 1;

[0190] 6.2) Increment the iteration count k by one, and use the interior point method to solve problem P2 to obtain the solution

[0191] where z (k+1) represents the solution vector of the (k + 1)-th iteration,

[0192] represents the solution of the state variable of the (k + 1)-th iteration, represents the solution of the state variable at the (k + 1)-th iteration at the n-th discrete point, n = 0, 1, 2... N + 1,

[0193] represents the solution of the control quantity of the (k + 1)-th iteration, represents the solution of the control quantity at the (k + 1)-th iteration at the m-th discrete point, m = 1, 2, 3... N,

[0194] and represents the solution of the slack variable of the (k + 1)-th iteration;

[0195] 6.3) Let a five - dimensional column vector ε represent the convergence judgment condition, and judge whether the state variable solutions at each discrete point satisfy the convergence condition:

[0196] If not, let x (k) = x (k+1) , and return to step 6.2;

[0197] Otherwise, output the optimal solution The calculation is completed.

[0198] Finally, according to u in the optimal solution (k+1) obtain the nominal optimal reference trajectory of the control quantity, and according to x in the optimal solution (k +1) obtain the nominal optimal reference trajectory of the state quantity, which includes 4 state quantity curves of the four state variables of the vehicle's geocentric distance, longitude, latitude, and heading angle and 1 control quantity curve.

[0199] Step 7, establish a deviation model and solve the non - nominal optimal reference trajectory.

[0200] 7.1) Assume that the true aerodynamic parameters follow a normal distribution relative to the nominal aerodynamic parameters, and take the nominal aerodynamic parameters as the mean of the normal distribution. Set the standard deviation σ of the normal distribution to satisfy that the range corresponding to 3σ is [-20%, 20%], and obtain the normal distribution function of the aerodynamic parameters as the deviation model according to the mean and the standard deviation;

[0201] 7.2) Randomly generate 500 groups of aerodynamic parameters in the deviation model, and for the second - order cone programming problems of the sequences corresponding to these aerodynamic parameters, repeat step 6 for offline solution to obtain the non - nominal optimal reference trajectories under different parameter conditions, which include 4 state quantity curves of the four state variables of the vehicle's geocentric distance, longitude, latitude, and heading angle and 1 control quantity curve. Combine with the nominal optimal reference trajectory obtained in step 6 to obtain the control quantity envelope of these 501 groups of optimal reference trajectories, as Figure 5 shown, and obtain the state quantity envelope of these 501 groups of optimal reference trajectories, as Figure 4 shown, where Figure 4 (a) is the altitude envelope, Figure 4 (b) is the longitude - latitude envelope, Figure 4 (c) is the heading - angle envelope.

[0202] Step 8, obtain the state quantity data set X and the control quantity data set Y1.

[0203] 8.1) Sample the 4 state quantity curves and 1 control quantity curve in the nominal optimal reference trajectory obtained in step 6 respectively to obtain the nominal state quantity data set X containing the four state variables of the vehicle's geocentric distance, longitude, latitude, and heading angleb and the nominal control quantity dataset Y b , for convenience of representation, only the schematic diagrams of one state component and the control quantity are drawn, as shown in Figure 6 shown;

[0204] 8.2) Sample the 4 state quantity curves and 1 control quantity curve in the non-nominal optimal reference trajectory obtained in step 7 respectively to obtain the non-nominal state quantity dataset X containing the four state variables of the geocentric distance, longitude, latitude, and heading angle of the aircraft c and the non-nominal control quantity dataset Y c ;

[0205] 8.3) Combine the nominal state quantity dataset X b and the non-nominal state quantity dataset X c to obtain the state quantity dataset X, and combine the nominal control quantity dataset Y b and the non-nominal control quantity dataset Y c to obtain the control quantity dataset Y1.

[0206] Step 9, construct a neural network.

[0207] 9.1) Construct a neural network A with an input layer, two hidden layers, and an output layer connected in sequence, and set its loss function as: Loss(Y1,U1)=(Y1-U1) 2 , where Y1 represents the control quantity dataset, and U1 represents the output value of the neural network A, that is, the control quantity bank angle;

[0208] 9.2) Set the number of neurons in the input layer of the neural network A to 4, the number of neurons in the output layer of the neural network A to 1, and the activation function to a linear function. The two hidden layers each contain 30 neurons and 15 neurons, and the activation function is the tanh function.

[0209] Step 10, use the BP algorithm to perform offline training on the neural network A to obtain the trained trajectory network A'.

[0210] 10.1) Initialize the weight parameter β of the neural network A;

[0211] 10.2) Take out a part of the state quantity dataset X as the training trajectory set X1 and input it into the neural network A for control quantity prediction to obtain the output result U1 = A(X1), and calculate the loss value Loss of the neural network A according to the neural network output result U1 and the control quantity dataset Y1;

[0212] 10.3) Adopt the backpropagation method to calculate the network parameter gradient of the neural network A through the loss value Loss of the neural network A, and use the gradient descent algorithm to update the network parameters of the neural network A according to the network parameter gradient;

[0213] 10.4) Repeat steps 10.2) and 10.3) until the loss function of neural network A converges to the minimum value, and the trained trajectory network A' is obtained.

[0214] Step 11: Online obtain the trajectory optimization results of the unpowered reusable launch vehicle (RLV) during the reentry phase.

[0215] During the flight phase of the RLV, the on-board computer reads the real-time flight state variables measured by the RLV navigation system, and uses the real-time flight state variables as the input of the trained trajectory network A' for forward propagation to obtain the real-time control variables.

[0216] Connect the state variables and control variables read each time into lines respectively to obtain the optimal trajectory of the RLV during the reentry phase.

[0217] The technical effects of the present invention are further described below in combination with simulation experiments:

[0218] I. Simulation conditions

[0219] The initial conditions of the aircraft are as follows:

[0220] Initial geocentric distance Initial longitude Initial latitude Initial track angle Initial heading angle Initial speed

[0221] The terminal conditions of the aircraft are as follows:

[0222] Terminal geocentric distance Terminal longitude Terminal latitude Terminal track angle Terminal heading angle Terminal speed

[0223] All numerical simulations are carried out on a PC with an I7-8700F (3.20 GHz) CPU, and the software environment is MATLAB.

[0224] II. Simulation content

[0225] Simulation experiment 1: Under the above simulation conditions, input the nominal state variable dataset X b into the trajectory network A', and test the output of the neural network. The results are as Figure 7 shown. As can be seen from Figure 7 , the bank angle output by the network is highly consistent with the bank angle obtained by convex programming solution, indicating that the neural network has high accuracy under nominal parameter conditions.

[0226] Simulation Experiment 2: Under the above simulation conditions, input the nominal state variable dataset X b into the trajectory network A′, and according to the four error functions in the following formula, count the error of the output of the trajectory network A′ relative to the nominal control variable dataset Y b :

[0227]

[0228]

[0229] where n represents the number of sampling points of this group of trajectories, the superscript i corresponds to the sampling point ordinal number, and σ Optim represents the optimal bank angle obtained by the convex programming method, and σ NN represents the bank angle obtained by neural network fitting. Δσ AAE , Δσ MAE , Δσ ARE , Δσ MRE are respectively the mean absolute error, maximum absolute error, mean relative error, and maximum relative error of the bank angle of this group of trajectories. The results are shown in Table 1:

[0230] Table 1 Four errors of the bank angle fitted by neural network A′ under nominal conditions

[0231]

[0232] As can be seen from Table 1, under nominal parameter conditions, the error of the neural network is small and the accuracy is high.

[0233] Simulation Experiment 3: Under the above simulation conditions, use the remaining part of the state variable dataset X that has not been used as the training set as the test dataset. Input the training dataset, test dataset, and all state variable datasets into the trajectory network A′ respectively, and use the error function of Simulation Experiment 2 to count the error of the output of the trajectory network A′ relative to the control variable dataset Y1. The test results are shown in Table 2.

[0234] Table 2 Respective average values and maximum values of the four errors under non - nominal parameter conditions

[0235]

[0236] It can be seen from the statistical data in the above table that under non - nominal parameter conditions, the overall fitting effect of the neural network is good.

[0237] Simulation Experiment 4: Under the above simulation conditions, simulate the generation of the altitude trajectory of the trajectory network under nominal parameters, and the results are as Figure 8 shown.

[0238] From Figure 8It can be seen that the altitude trajectory fitted by the neural network basically coincides with the optimal altitude trajectory obtained by the convex programming method under nominal parameters.

[0239] Simulation Experiment 5: Under the above simulation conditions, simulate the generation of the longitude and latitude trajectories of the trajectory network under nominal parameters, and the results are as Figure 9 shown.

[0240] From Figure 9 it can be seen that the longitude and latitude trajectories fitted by the neural network basically coincide with the optimal longitude and latitude trajectories obtained by the convex programming method under nominal parameters.

[0241] Simulation Experiment 6: Under the above simulation conditions, bias the lift parameter and drag parameter by +20% and +10% respectively, and simulate the generation of the altitude trajectory of the trajectory network under non-nominal parameters. The results are as Figure 10 shown.

[0242] From Figure 10 it can be seen that the altitude trajectory quickly generated by the neural network has a good degree of coincidence with the optimal altitude trajectory obtained by the convex programming method under non-nominal parameters.

[0243] Simulation Experiment 7: Under the above simulation conditions, bias the lift parameter and drag parameter by +20% and +10% respectively, and simulate the generation of the longitude and latitude trajectories of the trajectory network under non-nominal parameters. The results are as Figure 11 shown.

[0244] From Figure 11 it can be seen that the longitude and latitude trajectories quickly generated by the neural network have a good degree of coincidence with the optimal longitude and latitude trajectories obtained by the convex programming method under non-nominal parameters.

[0245] Simulation Experiment 8: Under the above simulation conditions, test the terminal error of the trajectories generated online by the trajectory network for several groups of different lift parameters and drag parameters with respect to the optimal trajectories obtained by the convex programming method in the same situation, and the time required to generate the trajectories. The results are shown in Table 3.

[0246] Table 3 Terminal Error of Neural Network Trajectories Relative to Optimal Trajectories and Time Required to Generate Trajectories

[0247]

[0248] As can be seen from Table 3, the trajectories generated by the neural network have high terminal accuracy and can well adapt to the influence of changes in aerodynamic parameters. Moreover, the time required to generate trajectories under different lift parameter and drag parameter deviations does not exceed 1.06 s, which is significantly shorter than the calculation time of more than ten seconds for traditional offline calculation methods.

[0249] The results of the above simulation experiments show that the solution proposed by the present invention has a fast solving speed and better real-time performance, and can also adapt to parameter deviations such as lift and drag during the flight of the aircraft.

[0250] The specific embodiments described above are only examples of the present invention and do not constitute a limitation to the present invention. Without departing from the principle of the present invention, various deductions and changes in form or details are still within the scope of the claims of the present invention.

Claims

1. A method for optimizing the trajectory of an aircraft during reentry based on a neural network, characterized in that, It includes the following steps: (1) Describe the trajectory optimization of the reentry section of the aircraft as a continuous optimal control problem P0 composed of a mathematical model, boundary conditions, admissible controls, performance indicators, and process constraints; (2) Conduct convexification processing on P0 by means of form replacement, slack variable, softened constraint, and successive linearization method to obtain a sequential convex optimal control problem P1, and use the pseudospectral method to perform discrete parameterization processing on this P1 to obtain a sequential second-order cone programming problem P2; (3) Solve the sequential second-order cone programming problem P2 by the interior point method to obtain a nominal optimal reference trajectory including 4 aircraft state quantity curves and 1 aircraft control quantity curve; (4) Establish a deviation model by perturbing the aircraft aerodynamic parameters, and perform offline solution for each set of aerodynamic parameters in the deviation model to obtain a non-nominal optimal reference trajectory including 4 aircraft state quantity curves and 1 aircraft control quantity curve under different parameter conditions; (5) Sample the state quantity curves in the nominal optimal reference trajectory in step (3) and the non-nominal optimal reference trajectory in step (4) respectively to obtain a state quantity data set X containing four state variables of the aircraft's geocentric distance, longitude, latitude, and heading angle, and sample the control quantity curve in the optimal reference trajectory to obtain a control quantity data set Y1; (6) Construct neural network A composed of cascaded input layer, two hidden layers, and output layer in sequence, and set its loss function as: Loss(Y1, U1) = (Y1 - U1) 2 , where Y1 represents the control quantity data set, and U1 represents the output value of neural network A, that is, the control quantity tilt angle; (7) Take out a part of the state quantity data set X in step (5) as a training trajectory set X1 and input it into neural network A, and perform offline training on it using the BP algorithm. When the loss function of neural network A converges to a minimum value, obtain a trained trajectory network A'; (8) Obtain the trajectory optimization result of the reentry section of the aircraft online: During the flight of the aircraft, the on-board computer reads the real-time flight state quantity measured by the aircraft navigation system, uses the real-time flight state quantity as the input of the trained trajectory network A' for forward propagation, and obtains the real-time control quantity; Connect the state quantity and control quantity read each time into lines respectively to obtain the optimal trajectory of the reentry section of the aircraft.

2. The method according to claim 1, wherein: The continuous optimal control problem P0 in step (1) is expressed as: Among them, f1 is the state differential equation, which is expressed as: r, V, Ω, L, and D are dimensionless geocentric distance, velocity, angular velocity of the Earth's rotation, and lift and drag accelerations, respectively. J represents the performance index of the continuous optimal control problem, and x represents the five state variables of the aircraft. and represent the initial and final dimensionless states of the aircraft, respectively. e represents the dimensionless energy of the aircraft, and e0 and e f represent the initial and final dimensionless energies of the aircraft, respectively. θ and φ are longitude and latitude, respectively. R0 is the radius of the Earth, γ is the flight path angle, ψ is the heading angle, σ is the bank angle, and σ min and σ max are the lower and upper limits of the bank angle, respectively. q and n are the heat flux rate, dynamic pressure, and overload of the aircraft, respectively. q max , n max are the maximum values of the heat flux rate, dynamic pressure, and overload that the aircraft can withstand, respectively. k Q and ρ represent the heat flux rate calculation coefficient and atmospheric density, respectively.

3. The method according to claim 1, characterized in that: The sequential convex optimal control problem P1 obtained in step (2) is expressed as: s.t. x′ = A(x (k) )x + B(x (k) )u + b(x (k) ) |x(e)-x (k) (e)|≤δ where \(u = [u_1, u_2]=[\cos\sigma,\sin\sigma]\) represents the control variable, \(r\) and \(V\) are the dimensionless geocentric distance and velocity respectively, \(J\) represents the performance index of the continuous optimal control problem, \(x\) represents the five state variables of the aircraft, represents the dimensionless initial state of the aircraft, \(e\) represents the dimensionless energy of the aircraft, \(e_0\) and \(e\) f represent the initial and final dimensionless energies of the aircraft respectively, \(\theta\) and \(\varphi\) are the longitude and latitude respectively, \(\gamma\) is the flight path angle, \(\psi\) is the heading angle, \(\sigma\) is the bank angle, \(\sigma\) min and \(\sigma\) max are the lower and upper limits of the bank angle respectively, \(\lambda_{\theta}\) and \(\lambda_{\varphi}\) are the slack variables of the longitude and latitude respectively, the superscript \((k)\) indicates that the variable is in the \(k\)-th iteration, \(k\) represents the iteration count, \(A(x\) (k) ) and \(B(x\) (k) ) represent the coefficient matrices of the state variable and the control variable respectively, \(b(x\) (k) ) represents the co-vector matrix, \(c\) θ 、 \(c_{\theta}\) and \(c_{\varphi}\) are the penalty coefficients of the longitude and latitude respectively, \(c_1\) and \(c_0\) are the first-order term and the remainder term of the dimensionless geocentric distance \(r\) in the successive linearization process of the performance index, \(\varepsilon\) ψ represents the coefficient of the heading angle regularization term, \(\delta\) represents the trust region constraint vector, \(g\) represents the convexified process constraint, \(\omega\) l 、\(\omega\) h are the cosine values corresponding to the upper and lower limits of the bank angle respectively, \(r_f\), \(\theta_f\), \(\varphi_f\), \(\gamma_f\) and \(\psi_f\) represent the dimensionless geocentric distance, longitude, latitude, flight path angle and heading angle of the final state of the aircraft respectively.

4. The method according to claim 1, wherein: The sequential second-order cone programming problem P2 obtained by performing discrete parameterization processing on this P1 using the pseudospectral method in step (2) is expressed as: Among them, \(u = [u_1, u_2]=[\cos\sigma,\sin\sigma]\) represents the control variable, \(r\) and \(V\) are the dimensionless geocentric distance and velocity respectively, \(J\) represents the performance index of the continuous optimal control problem, \(x\) represents the five state variables of the aircraft, represents the dimensionless initial state of the aircraft, \(e\) represents the dimensionless energy of the aircraft, \(e_0\) and \(e\) f represent the initial and final dimensionless energies of the aircraft respectively, \(\theta\) and \(\varphi\) are the longitude and latitude respectively, \(\gamma\) is the track angle, \(\psi\) is the course angle, \(\sigma\) is the bank angle, \(\sigma\) min and \(\sigma\) max are the lower and upper limits of the bank angle respectively, represent the slack variables of longitude and latitude respectively, the superscript \((k)\) indicates that the variable is in the \(k\) - th iteration, \(k\) represents the iteration count, \(A(x\) (k) ) and \(B(x\) (k) ) represent the coefficient matrices of the state variable and the control variable respectively, \(b(x\) (k) ) represents the co - vector matrix, \(c\) θ 、 represent the penalty coefficients of longitude and latitude respectively, \(c_1\) and \(c_0\) represent the first - order term and the remainder term of the dimensionless geocentric distance \(r\) in the successive linearization process of the performance index, \(\varepsilon\) ψ represents the coefficient of the course - angle regularization term, \(\delta\) represents the trust - region constraint vector, represents the convexified process constraint, \(\omega\) l 、\(\omega\) h represent the cosine values corresponding to the upper and lower limits of the bank angle respectively, represent the dimensionless geocentric distance, longitude, latitude, track angle and course angle of the final state of the aircraft respectively, \(n = 0,1,2,\cdots,N + 1\) represents each discrete point, \(m = 1,2,3,\cdots,N\) represents each LG point, \(N\) represents the number of selected LG points, \(w\) represents the integration weight corresponding to each LG point, \(D\) i represents the derivative of the basis function at the LG point.

5. The method according to claim 1, wherein: Step (3) uses the interior point method to solve the sequential second-order cone programming problem P2, and the implementation is as follows: (3a) Let the iteration count k = 0, and the initial state quantity at the given discrete points where represents the solution of the state variable at the k-th iteration at the n-th discrete point, n = 0, 1, 2... N + 1; (3b)Increment the iteration count k by one, and solve problem P2 using the interior point method to obtain the solution vector z for the (k + 1)-th iteration (k+1) : Among them, represents the state variable solution at the (k + 1)-th iteration, represents the state variable solution at the (k + 1)-th iteration at the n-th discrete point, where n = 0, 1, 2... N + 1. Denotes the control quantity solution at the (k + 1)-th iteration, Denotes the control quantity solution at the (k + 1)-th iteration at the m-th discrete point, where m = 1, 2, 3... N, and represents the solution of the slack variable at the (k + 1)-th iteration; (3c) Let a five-dimensional column vector ε represent the convergence judgment condition, and judge whether the state variable solutions at each discrete point satisfy the convergence condition: If not satisfied, let x (k) = x (k+1) , and return to (3b); Otherwise, output the optimal solution The calculation is completed.

6. The method according to claim 1, characterized in that: In step (4), establish a deviation model and perform offline solution to obtain a non-nominal optimal reference trajectory including the aircraft state quantity curve and the aircraft control quantity curve under different parameter conditions. The implementation is as follows: Assume that the true aerodynamic parameters follow a normal distribution relative to the nominal aerodynamic parameters, and the nominal aerodynamic parameters are the mean of the normal distribution. Set the standard deviation σ of the normal distribution to satisfy the range corresponding to 3σ as [-20%, 20%], and obtain the normal distribution function of the aerodynamic parameters according to the mean and standard deviation as the deviation model; Using this deviation model, 500 sets of aerodynamic parameters are randomly generated, and for the sequence second-order cone programming problems corresponding to these aerodynamic parameters, step (4) is repeated for offline solution to obtain the non-nominal optimal reference trajectories including 4 aircraft state quantity curves and 1 aircraft control quantity curve under different parameter conditions.

7. The method according to claim 1, characterized in that: In step (5), the state quantity data set X and the control quantity data set Y1 are obtained, and the implementation is as follows: (5a) Sample the four state quantity curves and one control quantity curve in the nominal optimal reference trajectory obtained in step (3) respectively to obtain a nominal state quantity dataset X containing four state variables: the geocentric distance, longitude, latitude, and heading angle of the aircraft b and a nominal control quantity dataset Y b ; (5b) Sample the four state variable curves and one control variable curve in the non-nominal optimal reference trajectory obtained in step (4) respectively, to obtain a non-nominal state variable dataset X containing the four state variables of the geocentric distance, longitude, latitude, and heading angle of the aircraft c and a non-nominal control variable dataset Y c ; (5c) Combine the nominal state quantity data set X b and the non-nominal state quantity data set X c to obtain the state quantity data set X. Combine the nominal control quantity data set Y b and the non-nominal control quantity data set Y c to obtain the control quantity data set Y1.

8. The method according to claim 1, wherein: The parameters of each layer in the neural network A constructed in step (6) are as follows: The number of neurons in the input layer is 4; The number of neurons in the output layer is 1 and the activation function is a linear function; The two hidden layers respectively contain 30 neurons and 15 neurons, and the activation functions of the two hidden layers are both tanh functions.

9. The method according to claim 1, characterized in that: In step (7), the BP algorithm is used to perform offline training on the neural network A, and the implementation is as follows: (9a) Initialize the weight parameters β of the neural network A; (9b) Use the training trajectory set X1 as the input of the neural network A for control quantity prediction, obtain the output result U1 = A(X1), and calculate the loss value Loss of the neural network A according to the neural network output result U1 and the control quantity data set Y1; (9c) Adopt the backpropagation method to calculate the network parameter gradient of the neural network A through the loss value Loss of the neural network A, and use the gradient descent algorithm to update the network parameters of the neural network A according to the network parameter gradient; (9d) Repeat steps 9b) and 9c) until the loss function of the neural network A converges to the minimum value to obtain the trained trajectory network A'.

Images