Mars entry trajectory planning-oriented mapping pseudo-spectral model prediction convex optimization method

Abstract

The invention discloses a mapping pseudo-spectral model prediction convex optimization method for Mars entry trajectory planning, and the method comprises the steps: converting a non-uniform Chebyshev-Gauss-Labatto collocation point into an approximately uniform Mapped Chebyshev-Gauss-Labatto collocation point through employing Kosloff-Tal-Eaer mapping, improving the ill-conditioned property of a differential matrix, and reducing the numerical oscillation; a gravity center interpolation format is constructed on the basis of Mapped Chebyshev-Gauss-Labatto collocation points, the state increment and the control correction are accurately approached, and zigzag fluctuation of a constant value control profile is inhibited; according to a mixed trajectory updating strategy, a reference trajectory is updated in combination with Lagrange interpolation and numerical integration, linearization error accumulation is reduced, and the high-precision task requirement is met. According to the method, the smoothness, efficiency and reliability of mars entering trajectory optimization are remarkably improved.

Description

Technical Field

[0001] The present invention relates to the technical field of aerospace trajectory optimization, and specifically designs a Mars entry trajectory planning method based on convex optimization of a mapped Chebyshev pseudo-spectral model prediction. Background Art

[0002] Currently, the main convex optimization methods used for trajectory planning in the aerospace field include lossless convex programming, concave-convex programming (or convex difference programming), sequential convex optimization, and polynomial optimization. However, almost all of these methods linearize the differential equations of state variables and treat state variables and control variables as the variables to be optimized. This inevitably leads to the large number of variables to be processed, low optimization scale and computational efficiency, and poor numerical accuracy.

[0003] To reduce the optimization scale and computational complexity, the model-predictive static optimization (MPSP) method establishes a sensitivity equation between state increments and control corrections, approximating the state increments as linear functions of the control corrections, effectively improving computational efficiency. However, this method, which uses a uniform discretization strategy and pseudo-spectral polynomial approximation technique, still suffers from the drawback of expanding the optimization problem size and insufficient integration of path constraints. While existing model-predictive convex optimization (MPCP) methods can address 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, and some methods fail to fully consider path constraints. Furthermore, the generation of the optimal control profile during algorithm initialization relies on external methods, further increasing preprocessing costs.

[0004] In recent years, methods combining pseudospectral methods with convex optimization have significantly reduced the number of optimization variables through high-precision point placement and reconstruction of sensitivity equations. However, these methods are still limited by the inherent drawbacks of non-uniform point placement and may suffer from numerical instability or insufficient accuracy. Summary of the Invention

[0005] The purpose of the present invention is to overcome the shortcomings of the existing technology and propose a mapped Chebyshev pseudo-spectral model predictive convex optimization (MCGL-PMPCP) method for Mars entry trajectory planning. This method controls the uniformity of CGL pseudo-spectral point assignment by introducing conformal mapping and barycentric rational Lagrange interpolation, while improving the stability of Lagrange interpolation. It can improve the oscillation problem of the control quantity while ensuring the accuracy of the converged solution, thereby fundamentally improving the algorithm performance.

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

[0007] The present invention provides a convex optimization method for predicting a mapped pseudo-spectral model for Mars entry trajectory planning, comprising the following steps:

[0008] S1. Establishment of dimensionless kinetic model

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

[0010]

[0011] Wherein, all physical quantities 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:

[0012]

[0013] 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;

[0014] Introduce the equation:

[0015]

[0016] 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:

[0017]

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

[0019]

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

[0021] 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; the control matrix Reflects 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:

[0022]

[0023] The specific expression of the parameters in the A matrix is:

[0024] a 14 =sinγ,a 15 =Vcosγ, 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 Lagrange polynomial as interpolation basis, perform initialization processing and time domain transformation on the linearized error system;

[0027] Select the Lagrange polynomial L i (τ) is used as the interpolation basis. Since the collocation points of the CGL pseudospectral method are distributed on the interval [-1, 1], the time domain of the corresponding problem is transformed to the interval [-1, 1], and the time variable is transformed. The time domain transformation formula used is:

[0028]

[0029] Among them, t represents time, t0 represents the starting time of the trajectory, t f represents the termination time of the trajectory, τ represents the CGL collocation point;

[0030] S4, transform the second type of CGL collocation points into MCGL nodes through Kosloff-Tal-Eaer conformal mapping, and obtain the corresponding interpolation nodes and interpolation weights;

[0031] S4.1 uses the second type of collocation to construct the CGL pseudo-spectral method, that is, the N-order Chebyshev polynomial P N (τ) = cos(Ncos -1 τ) zero point τ0<τ1<…<τ N, these zero points correspond to the key state points in the Mars entry process and are used to accurately describe and optimize the entry trajectory:

[0032] τ k =cos(πk / N), k=0,1,…,N

[0033] Where, t0=-1,t N =1,τ k represents the kth CGL collocation point for the description of the normalized time point;

[0034] S4.2 uses the Kosloff-Tal-Eaer conformal mapping to transform the coordinates of the second type of CGL collocation points to obtain the interpolation nodes:

[0035]

[0036] Where: α represents the conformal mapping parameter, which characterizes the uniformity of the mapped CGL points. When the value of α is closer to 1, the distribution of the mapped CGL points is more uniform. α is selected as 0.99.

[0037] In order to ensure the numerical accuracy of the mapped CGL pseudospectral method, Kosloff and Tal-Ezer proposed the following scheme for selecting the conformal mapping parameter α:

[0038]

[0039] Where: ε is the expected numerical calculation accuracy, N is the number of collocation points;

[0040] S4.3 uses the Clenshaw–Curtis integration rule and takes different integration weights according to the parity of the N order;

[0041] When N is an odd number, the corresponding integral weight is:

[0042]

[0043] When N is an even number, the corresponding integral weight is:

[0044]

[0045] The superscript " indicates that both the first and last terms in the summation sequence need to be divided by 2 to avoid repeated calculation of endpoints during discrete summation and improve integration accuracy. s =w N-s It reflects the symmetry of the nodes and realizes the accelerated convergence of the even function integral;

[0046] S5. Use the center of gravity Lagrange interpolation method to discretize the linear error system, approximate the state increment and control correction, and obtain the discrete state increment function δX(τ) and control correction function δU(τ) in the linearized error system;

[0047] S5.1 discretizes the linear error system at N+1 CGL collocation points to obtain N+1 discrete state increments and control corrections;

[0048] According to the interpolation node -1=τ0<···<τ N =1 and the corresponding function values ​​f(τ0),···,f(τ N+1 ), the function f(τ) is interpolated with the N-order Lagrange interpolation polynomial f N (τ) to approximate:

[0049]

[0050] S5.2 Based on the idea of ​​Lagrange interpolation, the definition of centroid weight is added to derive the improved centroid Lagrange interpolation formula;

[0051] Let l(τ)=(τ-τ0)(τ-τ1)…(τ-τ N ), and define the center of gravity weight w i for:

[0052]

[0053] Where τ is the CGL collocation point, τ i represents the i-th CGL collocation point, τ j represents the jth CGL point. During the rapid altitude descent phase of Mars entry, the CGL nodes are densely distributed, and the center of gravity weight will automatically amplify the contribution of the nodes in this area to improve the dynamic response accuracy. At the same time, the interpolation oscillation at the densely distributed nodes is suppressed to ensure a smooth transition of the roll angle.

[0054] The improved centroid Lagrange interpolation formula is obtained as follows:

[0055]

[0056] S5.3 considers reducing the oscillation of the interpolation solution at both ends and ensuring the approximate accuracy of the variables. The obtained center-of-gravity Lagrange interpolation formula is used, combined with the mapping CGL node after coordinate transformation, to place the state increment δx and the control correction δu at the interpolation nodes λ0,…,λ N The N-order interpolation polynomial is used to express the oscillation problem of the controlled quantity using weighted summation:

[0057]

[0058] Among them, λ represents the mapping CGL node, λ i represents the i-th mapped CGL node, which optimizes the uniformity of CGL point distribution. i Represents the centroid weight of the mapping CGL collocation point, and satisfies:

[0059] w0=0.5,w N =(-0.5) N ,w k =(-1) k ,k=1,2,…,N-1

[0060] Among them, w0 strengthens the contribution of the left endpoint to adapt to the high dynamic changes in the early stage of Mars entry; and w N Then the signs of the weights of the right endpoints are alternated to balance the accumulated error at the end;

[0061] S6. Discretize the linear error system, transform the differential equation constraints into a series of equality constraints, derive the mapping CGL pseudo-spectral sensitivity equation, and then define the objective function and constraints of the convex optimization problem;

[0062] S6.1 solves the derivatives of the state increment δX and the control correction δU at the mapping CGL collocation point λ and obtains:

[0063]

[0064] Among them, D ki It is the element of the kth row and ith column of the first-order differential matrix of the CGL pseudo-spectral method, reflecting the contribution of the state deviation and control correction at different collocation points to the change rate of the state deviation and control correction at the current collocation point, and satisfies:

[0065]

[0066] Among them, the weight ratio Reflects the node λ i and λ k The larger the weight, the stronger the dynamic coupling and the higher the local dependence of the derivative calculation; while the diagonal term D KK It serves as a compensation term contributed by other nodes to maintain the stability of the discretized system;

[0067] S6.2 transforms the linear error system into a series of algebraic equality constraints at the mapped CGL nodes:

[0068]

[0069] Among them, δX(λ i ) is the state increment of the i-th node, δ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;

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

[0071]

[0072] in, Represents the state increment at N+1 mapped CGL nodes, Represents the control correction at N+1 mapped CGL nodes:

[0073]

[0074] is the identity matrix;

[0075] At the same time, the constant matrix at the N+1 mapping CGL nodes and Expressed as:

[0076]

[0077] in, Indicates that N+1 mappings are made to the CGL nodes. The constant matrix is ​​a block diagonal matrix with principal diagonal elements, Indicates that N+1 mappings are made to the CGL nodes. The constant matrix is ​​a block diagonal matrix with principal diagonal elements;

[0078] S6.3 discretizes the linear error system at the mapped CGL collocation points and constructs the mapped CGL pseudo-spectral sensitivity equation:

[0079]

[0080] in: and are the discrete sequences of state increments and control corrections, and is the coefficient matrix of the linearized system;

[0081] S6.4 Define the objective function and constraints of a convex optimization problem;

[0082] The quadratic objective function is defined in the framework of the mapped CGL pseudo-spectral model predictive convex optimization method (MCGL-PMPCP), and the corresponding optimal control problem can be defined as:

[0083]

[0084] g(X k ,U k )≤0,h(X k ,U k )=0

[0085] Among them: Q, R, R u are respectively the semi-positive definite weight matrices, and represent the desired initial and final states, respectively, and represent the initial and terminal states of the reference trajectory, g(X k ,U k ) and h(X k ,U k ) is expressed as a convexified path constraint, which can be convexified using the first-order Taylor expansion method. 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 amount, and the terminal control correction has no effect on the optimal solution.

[0086] S7. Calculate the defect constant matrix, and then use the primal-dual interior point method to solve the corresponding optimal control problem to obtain the control correction δU k and state increment δX k ;

[0087] And update the state sequence Determine the state increment δX k Whether the set tolerance is met, that is, the following convergence conditions are met, and the solution is iterated;

[0088] sup‖δX k ‖ ∞ ≤ε,k=0,1,…,N

[0089] Where, ε>0, is the user-given value;

[0090] S8, if the state increment δX k If the convergence condition is not met, then Will U k As the input of the nonlinear system, the reference trajectory X is updated by combining Lagrange interpolation and numerical integration. k , the steps are as follows:

[0091] S8.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 a uniform time node, and the control sequence at this series of nodes is estimated;

[0092] S8.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 ;

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

[0094] S9, the discretized solution is continuously approximated. When the state increments meet the convergence conditions, the solution ends. Obtain the solution to the optimal control problem and optimize the trajectory.

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

[0096] (1) The present invention transforms the non-uniform CGL collocation into an approximately uniform MCGL collocation by introducing the Kosloff-Tal-Eaer conformal mapping and the barycentric Lagrange interpolation technique, thereby ensuring the uniformity of the discrete collocations and reducing the number of discrete nodes, thereby reducing the scale of the optimization problem and improving the computational efficiency.

[0097] (2) The present invention uses numerical integration to update the reference trajectory, 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 has a significant improvement effect on the numerical accuracy of the converged solution.

[0098] (3) The present invention addresses the problem of oscillation at both ends of the numerical solution caused by the pathological state of the original CGL pseudo-spectral differential matrix. Through improved pseudo-spectral point matching and interpolation techniques, the oscillation phenomenon of the constant control quantity when solving the Mars entry trajectory optimization problem is effectively suppressed. The algorithm also improves the interpolation stability of the stationary function and ensures the accuracy of the converged solution. BRIEF DESCRIPTION OF THE DRAWINGS

[0099] Figure 1 Principle flow chart of the method of the present invention.

[0100] Figure 2 Mars entry longitude-latitude profiles and velocity-altitude profiles for the four methods.

[0101] Figure 3 Mars entry track angle and heading angle profiles for the four methods.

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

[0103] Figure 5 Mars entry path constraint profiles for the four approaches.

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

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

[0106] In order to verify the advantages of the MCGL-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.

[0107] The present invention provides a Mars entry trajectory planning method based on the convex optimization of the mapped Chebyshev pseudo-spectral model prediction, which transforms the non-uniform CGL nodes into relatively uniform mapped CGL nodes and improves the pathological characteristics of the CGL differential matrix by expanding the analytical region of the function to be approximated. The flow chart is as follows Figure 1 As shown, it includes the following steps:

[0108] S1. Establishment of dimensionless kinetic model

[0109] The dimensionless factors that define length, acceleration, time, and velocity are the Mars reference radius R0, the Martian surface gravity acceleration g0, and the dimensionless time and dimensionless velocity Ignoring the rotation of Mars, we get the three-degree-of-freedom dimensionless particle dynamics model:

[0110]

[0111] Wherein, all physical quantities 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:

[0112]

[0113] 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;

[0114] Introduce the equation:

[0115]

[0116] 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:

[0117]

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

[0119]

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

[0121] 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; the control matrix Reflects 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:

[0122]

[0123] The specific expression of the parameters in the A matrix is:

[0124] a 14 =sinγ,a 15 =Vcosγ, a 44 =-D V ,

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

[0126] S3, select Lagrange polynomial as interpolation basis, perform initialization processing and time domain transformation on the linearized error system;

[0127] Select the Lagrange polynomial L i(τ) is used as the interpolation basis. Since the collocation points of the CGL pseudospectral method are distributed on the interval [-1, 1], the time domain of the corresponding problem is transformed to the interval [-1, 1], and the time variable is transformed. The time domain transformation formula used is:

[0128]

[0129] Among them, t represents time, t0 represents the starting time of the trajectory, t f represents the termination time of the trajectory, τ represents the CGL collocation point;

[0130] S4, transform the second type of CGL collocation points into MCGL nodes through Kosloff-Tal-Eaer conformal mapping, and obtain the corresponding interpolation nodes and interpolation weights;

[0131] S4.1 uses the second type of collocation to construct the CGL pseudo-spectral method, that is, the N-order Chebyshev polynomial P N (τ) = cos(Ncos -1 τ) zero point τ0<τ1<…<τ N , these zero points correspond to the key state points in the Mars entry process and are used to accurately describe and optimize the entry trajectory:

[0132] τ k =cos(πk / N), k=0,1,…,N

[0133] Where, t0=-1,t N =1,τ k represents the kth CGL collocation point for the description of the normalized time point;

[0134] S4.2 uses the Kosloff-Tal-Eaer conformal mapping to transform the coordinates of the second type of CGL collocation points to obtain the interpolation nodes:

[0135]

[0136] Where: α represents the conformal mapping parameter, which characterizes the uniformity of the mapped CGL points. When the value of α is closer to 1, the distribution of the mapped CGL points is more uniform. α is selected as 0.99.

[0137] In order to ensure the numerical accuracy of the mapped CGL pseudospectral method, Kosloff and Tal-Ezer proposed the following scheme for selecting the conformal mapping parameter α:

[0138]

[0139] Where: ε is the expected numerical calculation accuracy, N is the number of collocation points;

[0140] S4.3 uses the Clenshaw–Curtis integration rule and takes different integration weights according to the parity of the N order;

[0141] When N is an odd number, the corresponding integral weight is:

[0142]

[0143] When N is an even number, the corresponding integral weight is:

[0144]

[0145] The superscript " indicates that both the first and last terms in the summation sequence need to be divided by 2 to avoid repeated calculation of endpoints during discrete summation and improve integration accuracy. s =w N-s It reflects the symmetry of the nodes and realizes the accelerated convergence of the even function integral;

[0146] S5. Use the center of gravity Lagrange interpolation method to discretize the linear error system, approximate the state increment and control correction, and obtain the discrete state increment function δX(τ) and control correction function δU(τ) in the linearized error system;

[0147] S5.1 discretizes the linear error system at N+1 CGL collocation points to obtain N+1 discrete state increments and control corrections;

[0148] According to the interpolation node -1=τ0<···<τ N =1 and the corresponding function values ​​f(τ0),···,f(τ N+1 ), the function f(τ) is interpolated with the N-order Lagrange polynomial f N (τ) to approximate:

[0149]

[0150] S5.2 Based on the idea of ​​Lagrange interpolation, the definition of centroid weight is added to derive the improved centroid Lagrange interpolation formula;

[0151] Let l(τ)=(τ-τ0)(τ-τ1)…(τ-τ N ), and define the center of gravity weight w i for:

[0152]

[0153] Where τ is the CGL collocation point, τ i represents the i-th CGL collocation point, τ jrepresents the jth CGL point. During the rapid altitude descent phase of Mars entry, the CGL nodes are densely distributed, and the center of gravity weight will automatically amplify the contribution of the nodes in this area to improve the dynamic response accuracy. At the same time, the interpolation oscillation at the densely distributed nodes is suppressed to ensure a smooth transition of the roll angle.

[0154] The improved centroid Lagrange interpolation formula is obtained as follows:

[0155]

[0156] S5.3 considers reducing the oscillation of the interpolation solution at both ends and ensuring the approximate accuracy of the variables. The obtained center-of-gravity Lagrange interpolation formula is used, combined with the mapping CGL node after coordinate transformation, to place the state increment δx and the control correction δu at the interpolation nodes λ0,…,λ N The N-order interpolation polynomial is used to express the oscillation problem of the controlled quantity using weighted summation:

[0157]

[0158] Among them, λ represents the mapping CGL node, λ i represents the i-th mapped CGL node, which optimizes the uniformity of CGL point distribution. i Represents the centroid weight of the mapping CGL collocation point, and satisfies:

[0159] w0=0.5,w N =(-0.5) N ,w k =(-1) k ,k=1,2,…,N-1

[0160] Among them, w0 strengthens the contribution of the left endpoint to adapt to the high dynamic changes in the early stage of Mars entry; and w N Then the signs of the weights of the right endpoints are alternated to balance the accumulated error at the end;

[0161] S6. Discretize the linear error system, transform the differential equation constraints into a series of equality constraints, derive the mapping CGL pseudo-spectral sensitivity equation, and then define the objective function and constraints of the convex optimization problem;

[0162] S6.1 solves the derivatives of the state increment δX and the control correction δU at the mapping CGL collocation point λ and obtains:

[0163]

[0164] Among them, D kiIt is the element of the kth row and ith column of the first-order differential matrix of the CGL pseudo-spectral method, reflecting the contribution of the state deviation and control correction at different collocation points to the change rate of the state deviation and control correction at the current collocation point, and satisfies:

[0165]

[0166] Among them, the weight ratio Reflects the node λ i and λ k The larger the weight, the stronger the dynamic coupling and the higher the local dependence of the derivative calculation; while the diagonal term D KK It serves as a compensation term contributed by other nodes to maintain the stability of the discretized system;

[0167] S6.2 transforms the linear error system into a series of algebraic equality constraints at the mapped CGL nodes:

[0168]

[0169] Among them, δX(λ i ) is the state increment of the i-th node, δ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;

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

[0171]

[0172] in, Represents the state increment at N+1 mapped CGL nodes, Represents the control correction at N+1 mapped CGL nodes:

[0173]

[0174] is the identity matrix;

[0175] At the same time, the constant matrix at the N+1 mapping CGL nodes and Expressed as:

[0176]

[0177] in, Indicates that N+1 mappings are made to the CGL nodes. The constant matrix is ​​a block diagonal matrix with principal diagonal elements, Indicates that N+1 mappings are made to the CGL nodes. The constant matrix is ​​a block diagonal matrix with principal diagonal elements;

[0178] S6.3 discretizes the linear error system at the mapped CGL collocation points and constructs the mapped CGL pseudo-spectral sensitivity equation:

[0179]

[0180] in: and are the discrete sequences of state increments and control corrections, and is the coefficient matrix of the linearized system;

[0181] S6.4 Define the objective function and constraints of a convex optimization problem;

[0182] The quadratic objective function is defined in the framework of the mapped CGL pseudo-spectral model predictive convex optimization method (MCGL-PMPCP), and the corresponding optimal control problem can be defined as:

[0183]

[0184] g(X k ,U k )≤0,h(X k ,U k )=0

[0185] Among them: Q, R, R u are respectively the semi-positive definite weight matrices, and represent the desired initial and final states, respectively, and represent the initial and terminal states of the reference trajectory, g(X k ,U k ) and h(X k ,U k ) is expressed as a convexified path constraint, which can be convexified using the first-order Taylor expansion method. 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 amount, and the terminal control correction has no effect on the optimal solution.

[0186] S7. Calculate the defect constant matrix, and then use the primal-dual interior point method to solve the corresponding optimal control problem and obtain the control correction δU k and state increment δX k ;

[0187] And update the state sequence Determine the state increment δX kWhether the set tolerance is met, that is, the following convergence conditions are met, and the solution is iterated;

[0188] sup‖δX k ‖ ∞ ≤ε,k=0,1,…,N

[0189] Where, ε>0, is the user-given value;

[0190] S8, if the state increment δX k If the convergence condition is not met, then Will U k As the input of the nonlinear system, the reference trajectory X is updated by combining Lagrange interpolation and numerical integration. k , the steps are as follows:

[0191] S8.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 a uniform time node, and the control sequence at this series of nodes is estimated;

[0192] S8.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 ;

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

[0194] S9, the discretized solution is continuously approximated. When the state increments meet the convergence conditions, the solution ends. Obtain the solution to the optimal control problem and optimize the trajectory.

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

[0196] In order to verify the performance of the proposed MCGL-PMPCP method in avoiding the oscillation of the constant control profile, combined with Figures 2 to 7 The example verification of the present invention is described. In the simulation, the YALMIP toolbox is used to describe the convex optimization problem, and the solver MOSEK is used to solve the convex optimization problem.

[0197] 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.9m2 ,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.

[0198] Parameter description: R0 is the reference radius of Mars, g0 is the gravitational acceleration on the surface of Mars, ρ0 is the atmospheric density on the surface of Mars, h s is the altitude of the Martian atmospheric density model, m is the mass of the entry vehicle, S r is the reference area of ​​the inlet, k Q is the aerodynamic thermal coefficient, R n is the nose cone radius of the entry device, C L is the lift coefficient of the entry device, C D is the entry resistance, is the maximum heat flux, q max is the maximum dynamic pressure, a max is the maximum overload.

[0199] The MCGL-PMPCP method imposes virtual control and adds a penalty term for the virtual control amount in the objective function. The penalty weight w v Set to 10 6 The conformal mapping parameter α is 0.99. The simulation results of LG-PMPCP, LGR-PMPCP, CGL-PMPCP and MCGL-PMPCP are compared when the reference trajectory is updated using the numerical integration method.

[0200] It is not difficult to see from the simulation results that both the CGL-PMPCP and MCGL-PMPCP methods are very close to the simulation results of LG-PMPCP and LGR-PMPCP. This shows that the four methods can converge to the optimal solution. The main difference lies in the computational efficiency of each method and the smoothness of the control volume profile. Figure 4It can be seen that only the roll angular velocity profile of MCGL-PMPCP is relatively smooth, while the roll angular velocity profiles of the other three methods all show a sawtooth shape. The reason is that the MCGL-PMPCP method uses the Kosloff-Tal-Eaer conformal mapping to improve the non-uniformly distributed second-type CGL points into approximately uniform points, and then uses the centroid Lagrange interpolation to replace the Lagrange interpolation method, thereby alleviating the oscillation at both ends of the numerical solution caused by the ill-conditioned original CGL pseudo-spectral differential matrix, and also making the algorithm more stable for the interpolation of stationary functions. This improvement is reflected in the simulation results that the initial (0~80s) and middle (130~170s) parts of the roll angular velocity profile are smooth curves instead of sawtooth shapes. As for the roll angular velocity profile after 200s, since the roll angle actually shows an upward trend, the roll angular velocity is not a constant value, so the roll angular velocity oscillation degree of the four methods is relatively small. But from the perspective of computational efficiency, Figure 6 It shows that CGL-PMPCP requires the longest number of iterations, reaching 21 times, while MCGL-PMPCP only needs 6 iterations to converge. The data in Table 2 also show that the calculation time of the CGL-PMPCP method exceeds 10s, while the calculation time of MCGL-PMPCP and other Legendre pseudo-spectral model prediction convex optimization methods are basically the same. However, since the MCGL-PMPCP method uses conformal mapping to improve the uniformity of pseudo-spectral points, it will reduce the calculation accuracy ε, which in turn leads to a decrease in the accuracy of the algorithm convergence solution. This is in Figure 7 This is also verified in , where the MCGL-PMPCP method has the largest height error. Table 2 also shows that its corresponding height error is 0.1277 km, the maximum among the four methods. In summary, the MCGL-PMPCP method can effectively suppress the oscillation of the constant control variable, but the accuracy of the numerical solution is somewhat reduced.

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

[0202]

[0203] Table 2 Comparison of key parameters of four methods

[0204]

[0205]

Claims

1. A convex optimization method for mapping pseudo-spectral model prediction for Mars entry trajectory planning, characterized by The following steps are involved: S1. Establishment of dimensionless kinetic model The dimensionless factors that define length, acceleration, time, and velocity are the Mars reference radius R0, the Martian surface gravity acceleration g0, and the dimensionless time and dimensionless velocity Ignoring the rotation of Mars, we get the three-degree-of-freedom dimensionless particle dynamics model: Wherein, all physical quantities 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, 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; the control matrix Reflects 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 elevation of the Martian atmospheric density model, ρ0 is the atmospheric density on the surface of Mars; S3, select the Lagrange polynomial as the interpolation basis, perform initialization processing and time domain transformation on the linearized error system; Select the Lagrange polynomial L i (τ) is used as the interpolation basis. Since the collocation points of the CGL pseudospectral method are distributed on the interval [-1, 1], the time domain of the corresponding problem is transformed to the interval [-1, 1], and the time variable is transformed. The time domain transformation formula used is: Among them, t represents time, t0 represents the starting time of the trajectory, t f represents the termination time of the trajectory, τ represents the CGL collocation point; S4, transform the second type of CGL collocation points into MCGL nodes through Kosloff-Tal-Eaer conformal mapping, and obtain the corresponding interpolation nodes and interpolation weights; S4.1 uses the second type of collocation to construct the CGL pseudo-spectral method, that is, the N-order Chebyshev polynomial P N (τ) = cos(Ncos -1 τ) zero point τ0<τ1<…<τ N , these zero points correspond to the key state points in the Mars entry process and are used to accurately describe and optimize the entry trajectory: t k =cos(πk / N),k=0,1,…,N Where, t0=-1,t N =1,τ k represents the kth CGL collocation point for the description of the normalized time point; S4.2 uses the Kosloff-Tal-Eaer conformal mapping to transform the coordinates of the second type of CGL collocation points to obtain the interpolation nodes: Where: α represents the conformal mapping parameter, which characterizes the uniformity of the mapped CGL points. When the value of α is closer to 1, the distribution of the mapped CGL points is more uniform. α is selected as 0.

99. In order to ensure the numerical accuracy of the mapped CGL pseudospectral method, Kosloff and Tal-Ezer proposed the following scheme for selecting the conformal mapping parameter α: Where: ε is the expected numerical calculation accuracy, N is the number of collocation points; S4.3 uses the Clenshaw–Curtis integration rule and takes different integration weights according to the parity of the N order; When N is an odd number, the corresponding integral weight is: When N is an even number, the corresponding integral weight is: The superscript " indicates that both the first and last terms in the summation sequence need to be divided by 2 to avoid repeated calculation of endpoints during discrete summation and improve integration accuracy. s =w N-s It reflects the symmetry of the nodes and realizes the accelerated convergence of the even function integral; S5. Use the center of gravity Lagrange interpolation method to discretize the linear error system, approximate the state increment and control correction, and obtain the discrete state increment δx and control correction δu in the linear error system; S5.1 discretizes the linear error system at N+1 CGL collocation points to obtain N+1 discrete state increments and control corrections; According to the interpolation node -1=τ0<···<τ N =1 and the corresponding function values ​​f(τ0),···,f(τ N+1 ), the function f(τ) is interpolated with the N-order Lagrange interpolation polynomial f N (τ) to approximate: S5.2 Based on the idea of ​​Lagrange interpolation, the definition of centroid weight is added to derive the improved centroid Lagrange interpolation formula; Let l(τ)=(τ-τ0)(τ-τ1)…(τ-τ N ), and define the center of gravity weight w i for: Where τ is the CGL collocation point, τ i represents the i-th CGL collocation point, τ j represents the jth CGL point. During the rapid altitude descent phase of Mars entry, the CGL nodes are densely distributed, and the center of gravity weight will automatically amplify the contribution of the nodes in this area to improve the dynamic response accuracy. At the same time, the interpolation oscillation at the densely distributed nodes is suppressed to ensure a smooth transition of the roll angle. The improved centroid Lagrange interpolation formula is obtained as follows: S5.3 considers reducing the oscillation of the interpolation solution at both ends and ensuring the approximate accuracy of the variables. The obtained center-of-gravity Lagrange interpolation formula is used, combined with the mapping CGL node after coordinate transformation, to place the state increment δx and the control correction δu at the interpolation nodes λ0,…,λ N The N-order interpolation polynomial is used to express the oscillation problem of the controlled quantity using weighted summation: Among them, λ represents the mapping CGL node, λ i represents the i-th mapped CGL node, which optimizes the uniformity of CGL point distribution. i Represents the centroid weight of the mapping CGL collocation point, and satisfies: w0=0.5,w N =(-0.5) N ,In k =(-1) k ,k=1,2,…,N-1 Among them, w0 strengthens the contribution of the left endpoint to adapt to the high dynamic changes in the early stage of Mars entry; and w N Then the signs of the weights of the right endpoints are alternated to balance the accumulated error at the end; S6. Discretize the linear error system, transform the differential equation constraints into a series of equality constraints, derive the mapping CGL pseudo-spectral sensitivity equation, and then define the objective function and constraints of the convex optimization problem; S6.1 solves the derivatives of the state increment δX and the control correction δU at the mapping CGL collocation point λ and obtains: Among them, D ki It is the element of the kth row and ith column of the first-order differential matrix of the CGL pseudo-spectral method, reflecting the contribution of the state deviation and control correction at different collocation points to the change rate of the state deviation and control correction at the current collocation point, and satisfies: Among them, the weight ratio Reflects the node λ i and λ k The larger the weight, the stronger the dynamic coupling and the higher the local dependence of the derivative calculation; while the diagonal term D KK It serves as a compensation term contributed by other nodes to maintain the stability of the discretized system; S6.2 transforms the linear error system into a series of algebraic equality constraints at the mapped CGL nodes: Among them, δX(λ i ) is the state increment of the i-th node, δ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, Represents the state increment at N+1 mapped CGL nodes, Represents the control correction at N+1 mapped CGL nodes: is the identity matrix; At the same time, the constant matrix at the N+1 mapping CGL nodes and Expressed as: in, Indicates that N+1 mappings are made to the CGL nodes. The constant matrix is ​​a block diagonal matrix with principal diagonal elements, Indicates that N+1 mappings are made to the CGL nodes. The constant matrix is ​​a block diagonal matrix with principal diagonal elements; S6.3 discretizes the linear error system at the mapped CGL collocation points and constructs the mapped CGL pseudo-spectral sensitivity equation: in: and are the discrete sequences of state increments and control corrections, and is the coefficient matrix of the linearized system; S6.4 Define the objective function and constraints of a convex optimization problem; The quadratic objective function is defined in the framework of the mapped CGL pseudo-spectral model predictive convex optimization method (MCGL-PMPCP), and the corresponding optimal control problem can be defined as: g(X k ,U k )≤0,h(X k ,U k )=0 Among them: Q, R, R u are respectively the semi-positive definite weight matrices, and represent the desired initial and final states, respectively, and represent the initial and terminal states of the reference trajectory, g(X k ,U k ) and h(X k ,U k ) is expressed as a convexified path constraint, which can be convexified using the first-order Taylor expansion method. 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 amount, and the terminal control correction has no effect on the optimal solution. S7. Calculate the defect constant matrix, and then use the primal-dual interior point method to solve the corresponding optimal control problem and obtain the control correction δU k and state increment δX k ; And update the state sequence Determine the state increment δX k Whether the set tolerance is met, that is, the following convergence conditions are met, and the solution is iterated; sup‖δX k ‖ ∞ ≤ε,k=0,1,…,N Where, ε>0, is the user-given value; S8, if the state increment δX k If the convergence condition is not met, then Will U k As the input of the nonlinear system, the reference trajectory X is updated by combining Lagrange interpolation and numerical integration. k , the steps are as follows: S8.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 a uniform time node, and the control sequence at this series of nodes is estimated; S8.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 ; S8.3 uses Lagrange interpolation to estimate the state sequence at the CGL collocation point within the next optimization time interval based on the uniformly distributed state sequence, which is the reference trajectory for the next optimization; S9, the discretized solution is continuously approximated. When the state increments meet the convergence conditions, the solution ends. Obtain the solution to the optimal control problem and optimize the trajectory.

Images