LG pseudo-spectral model prediction convex optimization method for aircraft trajectory planning

Abstract

An LG pseudo-spectrum model prediction convex optimization method for aircraft trajectory planning belongs to the technical field of aircraft trajectory optimization and guidance, and adopts a Legendre-Gauss pseudo-spectrum method discrete dynamics system, establishes a pseudo-spectrum sensitivity equation through non-uniform Legendre-Gauss pseudo-spectrum collocation and a Gauss integral rule, reduces the number of variables to be optimized, and guarantees the numerical integration precision at the same time; a reference trajectory updating strategy is designed, the reference trajectory is updated in a mode of combining Lagrange interpolation and numerical integration, and the convergence of the algorithm and the optimality of a numerical solution are improved. Simulation results show that compared with traditional model prediction convex optimization and sequence convex optimization methods, the calculation time is reduced by more than 45% under the same precision, the height error is reduced by 60%, and path constraints and control constraints can be efficiently processed. The method is suitable for the Mars entry trajectory planning problem of the aircraft, and has the advantages of high precision and light weight.

Description

Technical Field

[0001] The present invention belongs to the technical field of aircraft trajectory optimization and guidance, and particularly relates to an LG pseudo-spectral model prediction convex optimization method for aircraft trajectory planning. Background Art

[0002] In the field of aerospace trajectory planning, existing convex optimization methods (such as lossless convex programming, concave-convex programming, sequential convex optimization, and polynomial optimization) generally use the technical means of linearizing the differential equations of state variables and using both the state variables and the control variables as optimization variables. This technical approach requires processing a high-dimensional variable space during the solution process, significantly increasing the scale and computational complexity of the optimization problem. For example, when dealing with trajectory planning problems, traditional sequential convex optimization methods need to discretize the state variables and control variables at uniform discrete nodes. This causes the number of decision variables to grow exponentially with the number of discrete nodes, which in turn affects the ability to solve in real time.

[0003] To reduce variable dimensionality, the model-predictive static optimization (MPSP) method constructs sensitivity equations between state increments and control corrections, approximating the state increments as linear functions of the control corrections, effectively reducing the number of optimization variables. However, this method, which uses a uniform discretization strategy and pseudo-spectral polynomial approximation techniques, still suffers from the drawback of bloating the optimization problem size and insufficient integration of path constraints. Furthermore, while existing model-predictive convex optimization (MPCP) methods can handle constrained optimal control problems, their use of uniform discretization schemes and Euler's integration rule requires a large number of discrete nodes to ensure numerical accuracy. Furthermore, algorithm initialization relies on the generation of external optimal control profiles, further increasing preprocessing costs. Summary of the Invention

[0004] The purpose of the present invention is to overcome the shortcomings of the existing technology and propose an LG pseudo-spectral model predictive convex optimization method (PMPCP) for aircraft trajectory planning. This method first introduces the Legendre-Gauss pseudo-spectral method to improve the existing model predictive convex optimization (MPCP) method, reduces the number of variables to be optimized, gives a pseudo-spectral sensitivity equation, and then uses the optimized solution as the input of the original nonlinear dynamics. The reference trajectory of the next optimization is updated by numerical integration, which solves the problems of many variables to be optimized and large optimization scale in traditional methods.

[0005] The present invention is achieved through the following technical solutions.

[0006] The present invention provides a LG pseudo-spectral model prediction convex optimization method for aircraft trajectory planning, comprising the following steps:

[0007] S1. Establishment of dimensionless kinetic model

[0008] The dimensionless factors defining length, acceleration, time and velocity are the reference radius of Mars R0, the gravitational acceleration of Mars g0, and the time and speed Ignoring the rotation of Mars, we get the three-degree-of-freedom dimensionless particle dynamics model:

[0009]

[0010] All physical quantities in the formula are dimensionless physical quantities, r is the distance from the center of Mars, θ and φ are the longitude and latitude of Mars respectively, γ, ψ, and σ are the track angle, heading angle, and bank angle respectively, ω0 is the angular velocity of Mars rotation, V is the relative velocity vector, L and D are the aerodynamic lift acceleration and aerodynamic drag acceleration of the entry vehicle respectively. The expressions of the dimensionless lift and drag accelerations L and D are:

[0011]

[0012] Where R0 is the reference radius of Mars, ρ is the density of the Martian atmosphere, and S r is the reference area of ​​the inlet, C L 、C D are the lift coefficient and drag coefficient of the entry device respectively, and m is the mass of the entry device.

[0013] Introduce the equation:

[0014]

[0015] Among them, the control quantity u represents the inclination acceleration, and the state quantity x is defined as [r,θ,φ,V,γ,ψ,σ] T , rewrite the dimensionless dynamic model into a general dynamic system form:

[0016]

[0017] S2. For nonlinear systems, establish a time domain [t0,t f ] on the linear error system:

[0018]

[0019] δx(t)=x(t)-x * (t),δu(t)=u(t)-u * (t)

[0020] Where δx(t) represents the state increment, that is, the difference between the current state x(t) and the reference trajectory state x * (t); δu(t) represents the control correction, that is, the difference between the current control input u(t) and the reference control input u * (t); the superscript * indicates the corresponding reference value. Control matrix To reflect the sensitivity of the control input to the change of the system state, the state matrix Reflecting the sensitivity of the system state to its own changes, we have:

[0021]

[0022] During the Mars entry process, since factors such as atmospheric density, drag and lift change over time, the elements of the matrix will also change over time, so it is not a constant matrix. The specific expression of the parameters in the A matrix is:

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

[0024] a 44 =-D V ,

[0025] Where: h s The Martian atmospheric density model elevation, ρ0 is the atmospheric density on the Martian surface.

[0026] S3. Select Legendre-Gauss pseudo-spectral collocation points according to Legendre polynomials, referred to as LG collocation points, discretize the state increment δx(t) and control correction δu(t) at the LG collocation points, and interpolate them using the Lagrange interpolation basis to obtain the corresponding state increment function δX(τ) and control correction function δU(τ);

[0027] S3.1 The Legendre polynomial The roots of LG are used as points, where the endpoints τ0 = -1 and τ N+1 =1 does not belong to the distribution point; N distribution points are unevenly distributed in the time interval, and are more concentrated at the two ends of the interval, which can better capture the dynamic changes of the aircraft;

[0028] S3.2 The linear error system is divided into N LG collocation points -1<τ1<…<τ N <1 and endpoint τ0=-1 and τ N+1 =1, and obtain N+2 discrete state increments and control corrections;

[0029] S3.3 The formula for interpolating discrete state increments and discrete control corrections is:

[0030]

[0031] Among them, δX(τ) represents the state increment function, δU(τ) represents the control correction function, τ represents the LG collocation point, τ i represents the i-th LG node, τ j represents the jth LG node; L i (τ) and is a polynomial built on these nodes, used to calculate the interpolation value of τ at any time point;

[0032] S4. Obtain the differential state expression based on the LG collocation points, perform time domain transformation on the linear error system model, and express the total state increment using the initial value of the state increment, the constant matrix, and the total control correction linear expression;

[0033] S4.1 Obtain the differential matrix based on LG collocation Element D ki , which is used to describe how the state increment δX changes with the control correction δU:

[0034]

[0035] where τ k represents the kth LG collocation point, which is the root of Legendre polynomial and is used to select the key time points in the optimization process; Derivative of the Legendre polynomial for the kth LG collocation point, D ki The value of is determined by the value of the reciprocal of the basis function at the collocation point;

[0036] S4.2 defines the time interval conversion through the transformation formula:

[0037]

[0038] The actual time interval [t0,t f ] is mapped to the standard interval [-1,1] to simplify subsequent mathematical operations and numerical calculations.

[0039] Among them, t represents time, t0 represents the starting time of the trajectory, t f Indicates the end time of the trajectory;

[0040] After S4.3 time domain transformation, the linear error system is transformed to obtain the discretized dynamic equation:

[0041]

[0042] Among them, δX(τ i ) is the state increment of the i-th collocation point, δX(τ k ) is the state increment of the k-th collocation point, δU(τ k ) is the control correction of the k-th collocation point, Jacobian matrices for state increments and control corrections;

[0043] By replacing and simplifying the equality constraint, we get a linear matrix equation as follows:

[0044]

[0045] in, Indicates the state increment at N LG points, represents the control correction at N LG points, is the identity matrix;

[0046] At the same time, the constant matrix at N LG collocation points and Expressed as:

[0047]

[0048] in, Indicates the point distribution with N LGs The constant matrix is ​​a block diagonal matrix with principal diagonal elements, Indicates the point distribution with N LGs The constant matrix is ​​a block diagonal matrix with principal diagonal elements;

[0049] S4.4 Combined with the corresponding constant matrix, the linear matrix equation is linearly transformed, and the total state increment is expressed by the initial value of the state increment, the constant matrix, and the total control correction, resulting in:

[0050]

[0051] Except for the end point τ N+1 The state increment δX N+1 Except for the initial state increment, the rest of the state increments can be modified by the initial state increment δX0 and the control Linear table output;

[0052] S5. Based on Gauss integral, the linear relationship between the terminal state increment, the initial state increment and the control correction is expressed, and the LG pseudo-spectral sensitivity equation is obtained; a convex optimization model is established to rewrite the optimal control problem;

[0053] S5.1 According to Gauss integral, the state increment function is in the actual time interval [t0,t f After the integral approximation on ], the state increment at the terminal N+1th matching point can be obtained:

[0054]

[0055] Among them, the weight matrix Indicates the weight of each collocation point in the integral; is the integral weight at the kth Gauss point, which is used to calculate the weighted sum of state increment and control correction; is the Nth-order Legendre polynomial P N (τ) at node τ k The first derivative at ;

[0056] S5.2 Combined with the total state increment expression, the terminal state increment δX is obtained. N+1 With the initial state increment δX0 and control correction The linear relationship between , that is, the LG pseudo-spectral sensitivity equation is:

[0057]

[0058] By using the LG pseudo-spectral sensitivity equation, the changes in state and control can be approximately calculated without directly solving the complex dynamic equations in the entire time interval, thereby simplifying the calculation process and improving the solution efficiency.

[0059] S5.3. Define a quadratic objective function under the LG-PMPCP framework of the LG pseudo-spectral model prediction convex optimization method. The objective function is to minimize the control correction δU k With the state increment δX k , the constraints include pseudo-spectral sensitivity equation and path constraint, then the corresponding optimal control problem can be defined as:

[0060]

[0061] in, represents the desired initial state, represents the desired end state, represents the initial state of the reference trajectory, Represents the terminal state of the reference trajectory, and the three weight matrices Q, R, R d is a semi-positive definite matrix; g(X k ,U k )≤0,h(X k ,U k )=0 represent the inequality and equality path constraints on the state and control variables, respectively, and can be convexified by a first-order Taylor expansion. In addition, the first term in the objective function represents the minimization of the state increment and control correction, while the second term is used to smooth the control variable, and the terminal control correction has no effect on the optimal solution.

[0062] S6. Calculate the defect constant matrix and use the primal-dual interior point method to solve the optimal control problem under the LG pseudo-spectral model predictive convex optimization control framework to obtain the control correction δU k and state increment δX k , and update the state sequence Determine the state increment δX k Whether the tolerance sup‖δX is met k ‖ ∞ ≤ε,k=1,…,N,ε>0, which is the user-set value. Iterative solution;

[0063] S7, if the state increment δX k If the set tolerance is not met, Will U k As the input of the nonlinear system, the reference trajectory is updated by combining Lagrange interpolation and numerical integration;

[0064] The reference trajectory update method combining Lagrange interpolation and numerical integration has the following specific steps:

[0065] S7.1 After the control sequence at the non-uniform pseudo-spectral distribution point is obtained by current optimization, the control sequence is first interpolated using Lagrange interpolation to map the control quantity at the non-uniform pseudo-spectral distribution point to uniform time nodes, and the control sequence at this series of nodes is estimated;

[0066] S7.2 uses the fixed-step Runge-Kutta integration rule to forward integrate the original nonlinear system based on the control sequence on the uniform node to generate a high-precision state sequence X k ;

[0067] S7.3 uses Lagrange interpolation to estimate the state sequence at the LG collocation point within the next optimization time interval based on the uniformly distributed state sequence, which is the reference trajectory for the next optimization;

[0068] S8. Continuously approximate the discretized solution. When the state increments meet the convergence conditions, the solution ends. Obtain the solution to the optimal control problem and optimize the trajectory.

[0069] The beneficial effects of the present invention compared with the prior art are:

[0070] (1) The present invention reduces the number of discrete nodes by introducing the LG pseudospectral method, adopting non-uniform point assignment and high-precision integration rules, thereby reducing the scale of the optimization problem and improving computational efficiency.

[0071] (2) The present invention updates the reference trajectory by numerical integration, avoiding the linearization error accumulation problem caused by directly using the optimized state sequence and control sequence as the reference trajectory for the next iteration, and improving the numerical accuracy of the converged solution. BRIEF DESCRIPTION OF THE DRAWINGS

[0072] Figure 1 1 is a flow chart of the principle of the method of the present invention.

[0073] Figure 2 is the rotation relationship between the coordinate systems.

[0074] Figure 3 Longitude-latitude profiles and velocity-altitude profiles for Mars entry for the four methods.

[0075] Figure 4 Track angle and heading angle profiles for the four Mars entry methods.

[0076] Figure 5 Mars entry tilt angle and tilt angular velocity profiles for the four methods.

[0077] Figure 6 Mars entry path constraint profiles for the four approaches.

[0078] Figure 7 Convergence histories of Mars entry control energy and virtual control for four methods.

[0079] Figure 8 Mars entry altitude profiles and altitude error profiles for the four methods. DETAILED DESCRIPTION

[0080] In order to verify the advantages of the LG-PMPCP method in terms of computational efficiency and numerical accuracy, the present invention will be further described in detail below with reference to the accompanying drawings and implementation examples.

[0081] The present invention provides an efficient convex optimization predictive control method integrating LG pseudo-spectrum improvement strategy, the flow chart of which is as follows: Figure 1 As shown, it includes the following steps:

[0082] S1. Establishment of dimensionless kinetic model

[0083] The dimensionless factors defining length, acceleration, time and velocity are the reference radius of Mars R0, the gravitational acceleration of Mars g0, and the time and speed Ignoring the rotation of Mars, we get the three-degree-of-freedom dimensionless particle dynamics model:

[0084]

[0085] All physical quantities in the formula are dimensionless physical quantities, r is the distance from the center of Mars, θ and φ are the longitude and latitude of Mars respectively, γ, ψ, and σ are the track angle, heading angle, and bank angle respectively, ω0 is the angular velocity of Mars rotation, V is the relative velocity vector, L and D are the aerodynamic lift acceleration and aerodynamic drag acceleration of the entry vehicle respectively. The expressions of the dimensionless lift and drag accelerations L and D are:

[0086]

[0087] Where R0 is the reference radius of Mars, ρ is the density of the Martian atmosphere, S r is the reference area of ​​the inlet, C L 、C D are the lift coefficient and drag coefficient of the entry device respectively, and m is the mass of the entry device.

[0088] Introduce the equation:

[0089]

[0090] Among them, the control quantity u represents the inclination acceleration, and the state quantity x is defined as [r,θ,φ,V,γ,ψ,σ] T , rewrite the dimensionless dynamic model into a general dynamic system form:

[0091]

[0092] S2. For nonlinear systems, establish a time domain [t0,t f ] on the linear error system:

[0093]

[0094] δx(t)=x(t)-x * (t),δu(t)=u(t)-u * (t)

[0095] Where δx(t) represents the state increment, that is, the difference between the current state x(t) and the reference trajectory state x * (t); δu(t) represents the control correction, that is, the difference between the current control input u(t) and the reference control input u * (t); the superscript * indicates the corresponding reference value. Control matrix To reflect the sensitivity of the control input to the change of the system state, the state matrix Reflecting the sensitivity of the system state to its own changes, we have:

[0096]

[0097] During the Mars entry process, since factors such as atmospheric density, drag and lift change over time, the elements of the matrix will also change over time, so it is not a constant matrix. The specific expression of the parameters in the A matrix is:

[0098] a 14 =sinγ,a 15 =Vcosγ,

[0099] a44 =-D V ,

[0100] Where: h s The Martian atmospheric density model elevation, ρ0 is the atmospheric density on the Martian surface.

[0101] S3. Select Legendre-Gauss pseudo-spectral collocation points according to Legendre polynomials, referred to as LG collocation points, discretize the state increment δx(t) and control correction δu(t) at the LG collocation points, and interpolate them using the Lagrange interpolation basis to obtain the corresponding state increment function δX(τ) and control correction function δU(τ);

[0102] S3.1 The Legendre polynomial The roots of LG are used as points, where the endpoints τ0 = -1 and τ N+1 =1 does not belong to the distribution point; N distribution points are unevenly distributed in the time interval, and are more concentrated at the two ends of the interval, which can better capture the dynamic changes of the aircraft;

[0103] S3.2 The linear error system is divided into N LG collocation points -1<τ1<…<τ N <1 and endpoint τ0=-1 and τ N+1 =1, and obtain N+2 discrete state increments and control corrections;

[0104] S3.3 The formula for interpolating discrete state increments and discrete control corrections is:

[0105]

[0106] Among them, δX(τ) represents the state increment function, δU(τ) represents the control correction function, τ represents the LG collocation point, τ i represents the i-th LG node, τ j represents the jth LG node; L i (τ) and is a polynomial built on these nodes, used to calculate the interpolation value of τ at any time point;

[0107] S4. Obtain the differential state expression based on the LG collocation points, perform time domain transformation on the linear error system model, and express the total state increment using the initial value of the state increment, the constant matrix, and the total control correction linear expression;

[0108] S4.1 Obtain the differential matrix based on LG collocation Element D ki , which is used to describe how the state increment δX changes with the control correction δU:

[0109]

[0110] where τ k represents the kth LG collocation point, which is the root of Legendre polynomial and is used to select the key time points in the optimization process; Derivative of the Legendre polynomial for the kth LG collocation point, D ki The value of is determined by the value of the reciprocal of the basis function at the collocation point;

[0111] S4.2 defines the time interval conversion through the transformation formula:

[0112]

[0113] The actual time interval [t0,t f ] is mapped to the standard interval [-1,1] to simplify subsequent mathematical operations and numerical calculations.

[0114] Among them, t represents time, t0 represents the starting time of the trajectory, t f Indicates the end time of the trajectory;

[0115] After S4.3 time domain transformation, the linear error system is transformed to obtain the discretized dynamic equation:

[0116]

[0117] Among them, δX(τ i ) is the state increment of the i-th collocation point, δX(τ k ) is the state increment of the k-th collocation point, δU(τ k ) is the control correction of the k-th collocation point, Jacobian matrices for state increments and control corrections;

[0118] By replacing and simplifying the equality constraint, we get a linear matrix equation as follows:

[0119]

[0120] in, Indicates the state increment at N LG points, represents the control correction at N LG points, is the identity matrix;

[0121] At the same time, the constant matrix at N LG collocation points and Expressed as:

[0122]

[0123] in, Indicates the point distribution with N LGs The constant matrix is ​​a block diagonal matrix with principal diagonal elements, Indicates the point distribution with N LGs The constant matrix is ​​a block diagonal matrix with principal diagonal elements;

[0124] S4.4 Combined with the corresponding constant matrix, the linear matrix equation is linearly transformed, and the total state increment is expressed by the initial value of the state increment, the constant matrix, and the total control correction, resulting in:

[0125]

[0126] Except for the end point τ N+1 The state increment δX N+1 Except for the initial state increment, the rest of the state increments can be modified by the initial state increment δX0 and the control Linear table output;

[0127] S5. Based on Gauss integral, the linear relationship between the terminal state increment, the initial state increment and the control correction is expressed, and the LG pseudo-spectral sensitivity equation is obtained; a convex optimization model is established to rewrite the optimal control problem;

[0128] S5.1 According to Gauss integral, the state increment function is in the actual time interval [t0,t f After the integral approximation on ], the state increment at the terminal N+1th matching point can be obtained:

[0129]

[0130] Among them, the weight matrix Indicates the weight of each collocation point in the integral; is the integral weight at the kth Gauss point, which is used to calculate the weighted sum of state increment and control correction; is the Nth-order Legendre polynomial P N (τ) at node τ k The first derivative at ;

[0131] S5.2 Combined with the total state increment expression, the terminal state increment δX is obtained. N+1 With the initial state increment δX0 and control correction The linear relationship between , that is, the LG pseudo-spectral sensitivity equation is:

[0132]

[0133] By using the LG pseudo-spectral sensitivity equation, the changes in state and control can be approximately calculated without directly solving the complex dynamic equations in the entire time interval, thereby simplifying the calculation process and improving the solution efficiency.

[0134] S5.3. Define a quadratic objective function under the LG-PMPCP framework of the LG pseudo-spectral model prediction convex optimization method. The objective function is to minimize the control correction δU k With the state increment δX k , the constraints include pseudo-spectral sensitivity equation and path constraint, then the corresponding optimal control problem can be defined as:

[0135]

[0136] in, represents the desired initial state, represents the desired end state, represents the initial state of the reference trajectory, Represents the terminal state of the reference trajectory, and the three weight matrices Q, R, R d is a semi-positive definite matrix; g(X k ,U k )≤0,h(X k ,U k )=0 represent the inequality and equality path constraints on the state and control variables, respectively, and can be convexified by a first-order Taylor expansion. In addition, the first term in the objective function represents the minimization of the state increment and control correction, while the second term is used to smooth the control variable, and the terminal control correction has no effect on the optimal solution.

[0137] S6. Calculate the defect constant matrix and use the primal-dual interior point method to solve the optimal control problem under the LG pseudo-spectral model predictive convex optimization control framework to obtain the control correction δU k and state increment δX k , and update the state sequence Determine the state increment δX k Whether the tolerance sup‖δX is met k ‖ ∞ ≤ε,k=1,…,N,ε>0, which is the user-set value. Iterative solution;

[0138] S7, if the state increment δX k If the set tolerance is not met, Will U k As the input of the nonlinear system, the reference trajectory is updated by combining Lagrange interpolation and numerical integration;

[0139] The reference trajectory update method combining Lagrange interpolation and numerical integration has the following specific steps:

[0140] S7.1 After the control sequence at the non-uniform pseudo-spectral distribution point is obtained by current optimization, the control sequence is first interpolated using Lagrange interpolation to map the control quantity at the non-uniform pseudo-spectral distribution point to uniform time nodes, and the control sequence at this series of nodes is estimated;

[0141] S7.2 uses the fixed-step Runge-Kutta integration rule to forward integrate the original nonlinear system based on the control sequence on the uniform node to generate a high-precision state sequence X k ;

[0142] S7.3 uses Lagrange interpolation to estimate the state sequence at the LG collocation point within the next optimization time interval based on the uniformly distributed state sequence, which is the reference trajectory for the next optimization;

[0143] S8. Continuously approximate the discretized solution. When the state increments meet the convergence conditions, the solution ends. Obtain the solution to the optimal control problem and optimize the trajectory.

[0144] Examples of the method of the present invention:

[0145] Combine Figures 3 to 8 To demonstrate the practical application of this invention, the YALMIP toolbox was used in the simulation to formulate the convex optimization problem, and the MOSEK solver was used to solve it. Because the PMPCP method employs a state error model, it can only optimize state increments and control corrections, not directly optimizing state variables or control variables. Therefore, it cannot solve the problem of maximizing the terminal altitude or minimizing the terminal velocity at Mars entry. Therefore, this simulation considers optimizing the control correction energy, with a fixed terminal time of 360 seconds.

[0146] The parameters in the simulation are set as R0 = 3397.2 km, g0 = 3.711 m / s 2 ,ρ0=0.0158kg / m 3 , h s =9354.5m, m=2804kg, S r =15.9m 2 ,k Q =1.9027×10 -4 ,R n =6.476,C L =0.36,C D =1.45, q max =8.5kPa,a max =18.5g0, and the initial values ​​and allowable variation ranges of other state quantities can be found in Table 1.

[0147] from Figure 3 and Figure 4 It can be seen that the optimal flight trajectories obtained by the four algorithms are almost the same, and the profiles including altitude, speed, track angle and heading angle are almost overlapping.

[0148] from Figure 5 It can be seen that the optimal roll angle and roll angular velocity profiles obtained by the MPCP and PMPCP methods are very close, and exhibit a profile similar to Bang-Bang control. However, since the objective function defined by the SCP method is different from that of the other two methods, its control profile has a consistent trend with the results of these two methods, but is not the same. In addition, the roll angular velocity of the PMPCP method is sawtooth-shaped, while the corresponding curves of the MPCP and SCP methods are relatively smooth. This phenomenon is that the pseudo-spectral method applies Lagrange interpolation to approximate the function at non-uniform points. It has high approximation accuracy for smooth functions with obvious fluctuations, but is prone to oscillation when approximating relatively stable functions (such as constant functions).

[0149] Figure 7 The control energy profile is given, where PMPCP and MPCP adopt Calculation of the corresponding control energy shows that the total amount of control energy required by MPCP and SCP is larger. This is actually because these two methods use a large number of uniform discrete nodes, resulting in a larger sum of the objective function.

[0150] Figure 8 The differences between the optimal trajectory obtained by the four methods and the optimal integral trajectory in the height channel are given. It can be seen that the MPCP and PMPCP methods have the highest maximum height error Δ|h| max The difference in this aspect is not large, which cannot reflect the advantages of the pseudo-spectral discrete format; however, SCP has the smallest height error. However, SCP needs to be optimized many times to obtain a result less than the convergence threshold, and the calculation time is more than twice that of other methods.

[0151] To more intuitively demonstrate the performance of the four algorithms, Table 2 lists the corresponding key parameters. For the SCP method, using numerical integration to update the reference trajectory did not significantly improve the numerical accuracy of the converged solution. The maximum altitude error was only reduced from 0.08309 km to 0.07954 km, the number of iterations was only reduced from 9 to 7, and the computation time was 7.6805 seconds. Overall, the different reference trajectory updating methods had little effect on the converged solution of the SCP method. However, for the MPCP and PMPCP methods, using numerical integration to update the reference trajectory significantly improved the accuracy of the converged solution. The maximum altitude error of the MPCP method was reduced to 0.05327 km, while that of the LG-PMPCP and LGR-PMPCP methods was reduced to 0.03162 km and 0.009845 km, respectively, significantly lower than that of the direct reference trajectory update method. However, using numerical integration to update the reference trajectory also increased the number of sequential optimizations for the MPCP and PMPCP methods, slightly increasing the computational complexity. However, the CPU time consumed by the PMPCP method was still lower than that of the MPCP and SCP methods.

[0152] Table 1 Numerical integration parameters of the dimensionless energy kinetic model

[0153]

[0154] Table 2 Comparison of key parameters of four methods when updating reference trajectory through numerical integration

[0155]

Claims

1. A LG pseudo-spectral model prediction convex optimization method for aircraft trajectory planning, characterized by The following steps are involved: S1. Establishment of dimensionless kinetic model The dimensionless factors defining length, acceleration, time and velocity are the reference radius of Mars R0, the gravitational acceleration of Mars g0, and the time and speed Ignoring the rotation of Mars, we get the three-degree-of-freedom dimensionless particle dynamics model: All physical quantities in the formula are dimensionless physical quantities, r is the distance from the center of Mars, θ and φ are the longitude and latitude of Mars respectively, γ, ψ, and σ are the track angle, heading angle, and bank angle respectively, ω0 is the angular velocity of Mars rotation, V is the relative velocity vector, L and D are the aerodynamic lift acceleration and aerodynamic drag acceleration of the entry vehicle respectively. The expressions of the dimensionless lift and drag accelerations L and D are: Where R0 is the reference radius of Mars, ρ is the density of the Martian atmosphere, and S r is the reference area of ​​the inlet, C L 、C D are the lift coefficient and drag coefficient of the entry device, respectively, and m is the mass of the entry device; Introduction equation: Among them, the control quantity u represents the inclination acceleration, and the state quantity x is defined as [r,θ,φ,V,γ,ψ,σ] T , rewrite the dimensionless dynamic model into a general dynamic system form: S2. For nonlinear systems, establish a time domain [t0,t f ] on the linear error system: δx(t)=x(t)-x * (t),δu(t)=u(t)-u * (t) Where δx(t) represents the state increment: the difference between the current state x(t) and the reference trajectory state x * (t); δu(t) represents the control correction: the current control input u(t) is different from the reference control input u * (t); the superscript * indicates the corresponding reference value; the control matrix To reflect the sensitivity of the control input to the change of the system state, the state matrix Reflecting the sensitivity of the system state to its own changes, we have: The specific expression of the parameters in the A matrix is: Where: h s The Martian atmospheric density model elevation, ρ0 is the atmospheric density on the Martian surface; S3. Select Legendre-Gauss pseudo-spectral collocation points according to Legendre polynomials, referred to as LG collocation points, discretize the state increment δx(t) and control correction δu(t) at the LG collocation points, and interpolate them using the Lagrange interpolation basis to obtain the corresponding state increment function δX(τ) and control correction function δU(τ); S3.1 The Legendre polynomial The roots of LG are used as points, where the endpoints τ0 = -1 and τ N+1 =1 does not belong to the distribution point; N distribution points are unevenly distributed in the time interval, and are more concentrated at the two ends of the interval, which can better capture the dynamic changes of the aircraft; S3.2 The linear error system is divided into N LG collocation points -1<τ1<…<τ N <1 and endpoint τ0=-1 and τ N+1 =1, and obtain N+2 discrete state increments and control corrections; S3.3 The formula for interpolating discrete state increments and discrete control corrections is: Among them, δX(τ) represents the state increment function, δU(τ) represents the control correction function, τ represents the LG collocation point, τ i represents the i-th LG node, τ j represents the jth LG node; L i (τ) and is a polynomial built on these nodes, used to calculate the interpolation value of τ at any time point; S4. Obtain the differential state expression based on the LG collocation points, perform time domain transformation on the linear error system model, and express the total state increment using the initial value of the state increment, the constant matrix, and the total control correction linear expression; S4.1 Obtain the differential matrix based on LG collocation Element D ki , which is used to describe how the state increment δX changes with the control correction δU: where τ k represents the kth LG collocation point, which is the root of Legendre polynomial and is used to select the key time points in the optimization process; Derivative of the Legendre polynomial for the kth LG collocation point, D ki The value of is determined by the value of the reciprocal of the basis function at the collocation point; S4.2 defines the time interval conversion through the transformation formula: The actual time interval [t0,t f ] is mapped to the standard interval [-1,1] to simplify subsequent mathematical operations and numerical calculations; Among them, t represents time, t0 represents the starting time of the trajectory, t f Indicates the end time of the trajectory; After S4.3 time domain transformation, the linear error system is transformed to obtain the discretized dynamic equation: Among them, δX(τ i ) is the state increment of the i-th collocation point, δX(τ k ) is the state increment of the k-th collocation point, δU(τ k ) is the control correction of the k-th collocation point, Jacobian matrices for state increments and control corrections; By replacing and simplifying the equality constraint, we get a linear matrix equation as follows: in, Indicates the state increment at N LG points, represents the control correction at N LG points, is the identity matrix; At the same time, the constant matrix at N LG collocation points and Expressed as: in, Indicates the point distribution with N LGs The constant matrix is ​​a block diagonal matrix with principal diagonal elements, Indicates the point distribution with N LGs The constant matrix is ​​a block diagonal matrix with principal diagonal elements; S4.4 Combined with the corresponding constant matrix, the linear matrix equation is linearly transformed, and the total state increment is expressed by the initial value of the state increment, the constant matrix, and the total control correction, resulting in: Except for the end point τ N+1 The state increment δX N+1 Except for the initial state increment, the rest of the state increments can be modified by the initial state increment δX0 and the control Linear expression; S5. Based on Gauss integral, the linear relationship between the terminal state increment, the initial state increment and the control correction is expressed, and the LG pseudo-spectral sensitivity equation is obtained; a convex optimization model is established to rewrite the optimal control problem; S5.1 According to Gauss integral, the state increment function is in the actual time interval [t0,t f After the integral approximation on ], the state increment at the terminal N+1th matching point can be obtained: Among them, the weight matrix Indicates the weight of each collocation point in the integral; k = 1, 2, ..., N is the integral weight at the kth Gauss collocation point, which is used to calculate the weighted sum of the state increment and control correction; is the Nth-order Legendre polynomial P N (τ) at node τ k The first derivative at ; S5.2 Combined with the total state increment expression, the terminal state increment δX is obtained. N+1 With the initial state increment δX0 and control correction The linear relationship between , that is, the LG pseudo-spectral sensitivity equation is: By using the LG pseudo-spectral sensitivity equation, the changes in state and control can be approximately calculated without directly solving the complex dynamic equations in the entire time interval, thereby simplifying the calculation process and improving the solution efficiency. S5.

3. Define a quadratic objective function under the LG-PMPCP framework of the LG pseudo-spectral model prediction convex optimization method. The objective function is to minimize the control correction δU k With the state increment δX k , the constraints include pseudo-spectral sensitivity equation and path constraint, then the corresponding optimal control problem can be defined as: in, represents the desired initial state, represents the desired end state, represents the initial state of the reference trajectory, Represents the terminal state of the reference trajectory, and the three weight matrices Q, R, R d is a semi-positive definite matrix; g(X k ,U k )≤0,h(X k ,U k )=0 represent the inequality and equality path constraints on the state and control variables, respectively, and can be convexified by a first-order Taylor expansion. In addition, the first term in the objective function represents the minimization of the state increment and control correction, while the second term is used to smooth the control variable, and the terminal control correction has no effect on the optimal solution. S6. Calculate the defect constant matrix and use the primal-dual interior point method to solve the optimal control problem under the LG pseudo-spectral model predictive convex optimization control framework to obtain the control correction δU k and state increment δX k , and update the state sequence Determine the state increment δX k Whether the tolerance sup‖δX is met k ‖ ∞ ≤ε, k=1,…,N, ε>0, which is the user-set value; iterative solution; S7, if the state increment δX k If the set tolerance is not met, Will U k As the input of the nonlinear system, the reference trajectory is updated by combining Lagrange interpolation and numerical integration; The reference trajectory update method combining Lagrange interpolation and numerical integration has the following specific steps: S7.1 After the control sequence at the non-uniform pseudo-spectral distribution point is obtained by current optimization, the control sequence is first interpolated using Lagrange interpolation to map the control quantity at the non-uniform pseudo-spectral distribution point to uniform time nodes, and the control sequence at this series of nodes is estimated; S7.2 uses the fixed-step Runge-Kutta integration rule to forward integrate the original nonlinear system based on the control sequence on the uniform node to generate a high-precision state sequence X k ; S7.3 uses Lagrange interpolation to estimate the state sequence at the LG collocation point within the next optimization time interval based on the uniformly distributed state sequence, which is the reference trajectory for the next optimization; S8. Continuously approximate the discretized solution. When the state increments meet the convergence conditions, the solution ends. Obtain the solution to the optimal control problem and optimize the trajectory.

Images