A rocket landing segment fast trajectory planning and guidance method based on flight airspace grid division

Abstract

This invention provides a rapid trajectory planning and guidance method for rocket landing phase based on flight airspace grid partitioning. It mainly includes establishing a rocket landing column coordinate system, flight airspace grid partitioning, a sequential convex optimization model for solving the absolute optimal trajectory, online interpolation to restore the absolute optimal trajectory, online sequential convex optimization trajectory optimization, and polynomial guidance. This invention aims to solve the problem of high-precision fixed-point landing of reusable launch vehicles' powered landing phase. By utilizing the landing column coordinate system, the three-dimensional motion of the rocket landing process is decomposed into two-dimensional motion in the longitudinal plane and tangential motion, reducing the complexity of the optimization model. An interpolation table is established using the flight airspace grid partitioning method, enabling rapid online acquisition of the absolute optimal trajectory. This transforms the global optimization problem into a local optimization problem, reducing the solution scale of the optimization model and improving the real-time performance of online optimization while satisfying various constraints and fuel optimization.

Description

Technical Field

[0001] This invention relates to a rapid trajectory planning and guidance method for rocket landing phase based on flight airspace grid division, applicable to trajectory planning and guidance tasks for the return landing phase of reusable launch vehicles, and belongs to the field of guidance, navigation and control technology for space transport vehicles. Background Technology

[0002] Online trajectory planning and guidance technology for rocket landing is one of the key technologies in the high-precision landing control system of reusable launch vehicles, directly affecting landing reliability and accuracy. Its main application is in the vertical recovery of the first stage of launch vehicles. The realization and engineering application of rocket recovery programs by Blue Origin and SpaceX in the United States signifies that vertical takeoff and landing vehicles have become a new approach to the development of reusable technology. Due to their own commercial interests, these companies have not disclosed many details about the online trajectory planning and guidance technology for recovery. Therefore, the design of my country's reusable launch vehicles must achieve breakthroughs in powered landing online trajectory planning and guidance technology. Space landing guidance technology has mainly gone through several stages of development. The first stage is reentry guidance based on nominal trajectory, which was used in early ballistic reentry landings of spacecraft. Since the spacecraft itself underwent uncontrolled reentry, the deviation was very large. The second stage is predictive guidance, used for landing guidance of currently operational spacecraft. It can adjust the landing point by adjusting the tilt angle, significantly reducing landing deviation. The third stage is polynomial guidance, mainly used for Mars-driven landings and lunar-driven soft landings. Polynomial guidance was used in the Chang'e-3 and Chang'e-4 missions. This method generates guidance commands by designing polynomial acceleration profiles, which is computationally simple and efficient. However, it cannot meet the various complex and strong constraints of the landing process, such as state constraints and fuel constraints, especially... For rocket reentry and soft landing missions within Earth's atmosphere, there are greater aerodynamic disturbances, and the constraints on engine thrust, thrust rate of change, and pitch angle are all more stringent. Furthermore, to increase rocket carrying capacity, it is essential to study how to land safely while minimizing fuel consumption. With the continuous improvement of onboard computer performance, a new stage of landing guidance technology has emerged, including computational guidance based on numerical methods, particularly landing guidance based on convex optimization theory. This has been extensively studied in both theory and engineering. Although convex optimization is renowned for its high efficiency, real-time performance remains the biggest obstacle to its practical engineering application, especially considering the comprehensive and stringent constraints of aerodynamic constraints, thrust constraints, and the requirement for high-precision soft landing.

[0003] Existing rocket landing trajectory planning and guidance methods often suffer from a trade-off between guidance accuracy, constraint satisfaction, and computational efficiency.

[0004] 1) Tracking guidance cannot adapt to large deviation conditions, making it difficult to achieve precise positioning and landing;

[0005] 2) Polynomial guidance considers fewer constraints, making it difficult to meet flight process constraints;

[0006] 3) Numerical guidance methods involve large amounts of computation, and under the conditions of onboard computers, real-time performance cannot be met; Summary of the Invention

[0007] The technical problem solved by this invention is to overcome the shortcomings of existing technologies and propose a rapid trajectory planning and guidance method for rocket landing phase based on flight airspace grid division. First, it can meet the dynamic constraints, state constraints, landing terminal constraints, engine thrust constraints, thrust change rate constraints, and yaw angle constraints during rocket landing, realizing trajectory planning under various comprehensive constraint conditions. Second, it can adapt to deviations in initial state, aerodynamic state, and other operating conditions, ensuring accurate landing under large deviation conditions while minimizing fuel consumption. Finally, this invention can reduce the complexity of numerical optimization and the solution scale of the optimization problem, and improve the real-time performance of online solutions.

[0008] The technical solution of this invention is:

[0009] A rapid trajectory planning and guidance method for rocket landing phase based on flight airspace grid division includes the following steps:

[0010] 1) Establish a landing column coordinate system and select the longitudinal section with a polar angle of zero as the reference section;

[0011] 2) Perform rectangular meshing on the reference profile selected in step 1) and store the position parameters of the mesh nodes;

[0012] 3) Establish a two-dimensional dynamic model of rocket landing in the reference profile selected in step 1);

[0013] 4) Based on the dynamic model established in step 3), construct a minimum fuel optimization model, and discretize and convexize it to obtain a convex optimization solution model for the standard landing trajectory planning problem sequence.

[0014] 5) Remove the initial velocity constraint from the standard optimization model established in step 4) to obtain the optimal ballistic solution model; substitute the grid node position parameters obtained in step 2) into the optimal ballistic solution model as initial position parameters for optimization and solution to obtain the absolute optimal ballistic parameters, and extract the absolute optimal ballistic initial velocity as the absolute optimal velocity parameter of the grid node, establish an interpolation table of grid node position parameters and absolute optimal velocity parameters, and bind it as parameters to the onboard computer;

[0015] 6) During flight, based on the rocket's current real-time position, the absolute optimal velocity is obtained by interpolating the rocket position parameters in step 5) using the absolute optimal velocity interpolation table. By performing trajectory integration, the rocket's position change is obtained. By performing cyclic interpolation and trajectory integration, the absolute optimal trajectory can be quickly reconstructed.

[0016] 7) Select a suitable point on the absolutely optimal trajectory obtained in step 6) as the target point; take the current position and current velocity of the rocket as the initial position and initial velocity, and the position and velocity of the target point as the final position and velocity, and substitute them into the sequential convex optimization solution model of the standard landing trajectory planning problem established in step 4). The local two-dimensional optimal trajectory of the rocket from the current point to the selected target point in the longitudinal flight plane can be obtained in the cylindrical coordinate system.

[0017] 8) Based on the two-dimensional trajectory obtained in step 7), the optimal trajectory is tracked using a fourth-order polynomial to obtain the magnitude and direction of the rocket's acceleration in the longitudinal plane of the cylindrical coordinate system; the remaining flight time of the rocket is estimated in the longitudinal plane; in the lateral direction of the cylindrical coordinate system, the tangential velocity component of the rocket is decelerated to zero in the remaining flight time using polynomial guidance to obtain the magnitude and direction of the tangential acceleration.

[0018] 9) Rotate the acceleration vector obtained in the cylindrical coordinate system to the fixed rectangular coordinate system at the landing point; calculate the rocket's three-axis attitude angle command based on the three-axis components of the acceleration command in the rectangular coordinate system; calculate the rocket's thrust magnitude command based on the magnitude of the acceleration command in the rectangular coordinate system.

[0019] 10) Determine if the current rocket flight altitude is lower than the preset altitude. If it is, end the online trajectory planning; if not, update the trajectory planning and guidance cycle, and return to step 5).

[0020] The absolute optimal ballistic solution model described in step 5) is in the following form:

[0021]

[0022] subject to r(n+1)=r(n)+0.5(v(n)+v(n+1))Δt+α r (n)

[0023] v(n+1)=v(n)+(0.5(T(n)-0.5ρS D C D ||v(n)||v(n)+T(n+1)-0.5ρS D C D ||v(n+1)||v(n+1)) / m(n))Δt+gΔt+α v (n)

[0024]

[0025] T min ≤Γ(n)≤T max

[0026]

[0027] m(1)=m0,r(1)=r0

[0028]

[0029] Where n is the number of the discrete point, n = 1, 2, 3, ..., N-1, N is the number of discrete points, r(n) is the position vector of the nth discrete point, v(n) is the velocity vector of the nth discrete point, T(n) is the rocket engine thrust vector at the nth discrete point, Γ(n) is the introduced scalar relaxation factor, m(n) is the rocket mass at the nth discrete point, ρ is the atmospheric density, and S... D For the rocket's aerodynamic reference area, C D Let α be the aerodynamic drag coefficient of the rocket, ||v(n)|| be the magnitude of the velocity vector at the nth discrete point, Δt be the time interval of the rocket optimization model, and α be the velocity vector at the nth discrete point. r (n) represents the position recursion virtual control term, α v (n) is the velocity recursion virtual control term, α m (n) is a quality recursive virtual control term.

[0030] In step 5), the position and velocity recursion of the discrete state equations are derived using the trapezoidal integral formula, as detailed below:

[0031] r(n+1)=r(n)+0.5(v(n)+v(n+1))Δt+α r (n)

[0032] v(n+1)=v(n)+(0.5(T(n)-0.5ρS D C D ||v(n)||v(n)+T(n+1)-0.5ρS D C D ||v(n+1)||v(n+1)) / m(n)+g)Δt+α v (n)

[0033] The method for restoring the absolute optimal trajectory described in step 6) is as follows:

[0034]

[0035] ρ co (1)=ρ co_now,y(1)=y now

[0036] [v ρco (i),v y (i)]=interp(ρ co (i),y(i))

[0037]

[0038] Among them, [x now ,y now ,z now [ρ] represents the position of the rocket in a fixed rectangular coordinate system at the landing point. co_now ,y now ,z co_now These represent the radial, longitudinal, and tangential positions in the rocket's column coordinate system, where the tangential position z... co_now =0, interp() is a linear interpolation function for the rocket position parameter-absolute optimal velocity interpolation table, the values ​​in parentheses are the target interpolation points, θ co T is the polar angle in cylindrical coordinates; integral This is the ballistic integration step size.

[0039] The method for selecting the target point described in step 7) is as follows:

[0040] The closer the target point is to the rocket's current position, the fewer discrete points are needed for online trajectory optimization, resulting in faster computation. However, the cycle time for online trajectory optimization needs to be very small to ensure frequent trajectory updates. Conversely, the farther the target point is from the rocket's current position, the more discrete points are needed for online trajectory optimization, resulting in slower computation, but the cycle time can be appropriately increased. Therefore, the overall design of the target point height is as follows:

[0041] y target =y now -2T TP ·v y_now

[0042] Among them, y target Let y be the height of the target point. now v is the current height. y_now T represents the current falling speed. TP The trajectory planning period.

[0043] The local two-dimensional optimal trajectory optimization model from the current position to the target position mentioned in step 7) is as follows:

[0044]

[0045] subject to r(n+1)=r(n)+0.5(v(n)+v(n+1))Δt+αr (n)

[0046] v(n+1)=v(n)+(0.5(T(n)-0.5ρS D C D ||v(n)||v(n)+T(n+1)-0.5ρS D C D ||v(n+1)||v(n+1)) / m(n))Δt+gΔt+α v (n)

[0047]

[0048] m(1)=m0,r(1)=r co0 ,v(1)=v co0 ,r(N)=r cof +r err v(N)=v cof +v err Where s1, s2, and s3 are the coefficient weights of the fuel consumption index, position deviation index, and speed deviation index, respectively; k is the coefficient weight of the virtual control term; r err To optimize the deviation between the trajectory terminal position and the target position; v err To optimize the deviation between the trajectory terminal velocity and the target velocity; r co0 =[ρ co0 [y0] represents the current position vector of the rocket on the longitudinal plane of the cylindrical coordinate system; r is the rocket's current velocity vector on the longitudinal plane of the cylindrical coordinate system. cof =[ρ cof ,y f ] is the position vector of the rocket target point on the longitudinal plane of the cylindrical coordinate system; Let be the velocity vector of the rocket target point on the longitudinal plane of the cylindrical coordinate system.

[0049] In the local two-dimensional optimal trajectory optimization model described in step 7), the weight design of the performance index is as follows:

[0050]

[0051] s2=1000||r err || / ||r cof ||

[0052] s3 = 1000||v err || / ||v cof ||

[0053] k = 1000

[0054] s1, s2 and s3 are variable weight designs; The three-axis velocity vectors are fixed in a rectangular coordinate system at the landing point; These are the radial velocity, longitudinal velocity, and tangential velocity in the rocket's cylindrical coordinate system; V err_now This is the magnitude of the deviation between the velocity vector in the longitudinal plane of the rocket's current cylindrical coordinate system and the optimal velocity vector obtained by interpolation at this point; the meanings of the other variables are the same as above, and all have been explained.

[0055] In the optimization process, this invention decomposes the optimization model into dimensions. The original trajectory optimization problem in three-dimensional space is decomposed into two-dimensional motion in the longitudinal plane and tangential motion in the cylindrical coordinate system by establishing a landing cylindrical coordinate system.

[0056] This invention transforms the global trajectory optimization problem into a local optimization problem. In the vertical plane of the cylindrical coordinate system, the decomposed global two-dimensional trajectory optimization problem from the rocket's real-time position to the landing point is transformed into a local trajectory optimization problem from the rocket's real-time position to a point on the absolutely optimal trajectory through the transition of the absolute optimal trajectory.

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

[0058] 1) This invention constructs a landing cylindrical coordinate system. Based on this coordinate system, the trajectory optimization problem in the original three-dimensional space is decomposed into two-dimensional motion in the longitudinal plane of the cylindrical coordinate system and tangential motion in the tangential direction of the cylindrical coordinate system. The two-dimensional trajectory sequence convex optimization in the plane is used to replace the three-dimensional trajectory sequence convex optimization in the rectangular coordinate system, thereby reducing the model complexity of trajectory planning.

[0059] 2) This invention performs flight airspace meshing and utilizes the decomposition of motion in the landing column coordinate system to transform the three-dimensional spatial meshing into a two-dimensional planar meshing, reducing the computational load of mesh nodes and the amount of mesh parameter storage data. It also enables online restoration through planar rotation, achieving rapid generation of reference trajectories based on mesh node interpolation in three-dimensional space.

[0060] 3) This invention introduces the absolutely optimal trajectory and the absolutely optimal velocity as intermediate transition states. Based on the Bellman optimality principle, the global optimal problem is transformed into a local optimal problem, which reduces the solution scale of online trajectory planning.

[0061] 4) This invention organically combines sequential convex optimization method and polynomial guidance method to improve the real-time performance of online applications while satisfying process constraints, terminal constraints and minimizing fuel consumption. Attached Figure Description

[0062] Figure 1 This is a flowchart illustrating the implementation of the present invention.

[0063] Figure 2This is a schematic diagram of the cylindrical coordinate system of the present invention.

[0064] Figure 3 This is a schematic diagram of the spatial location of the cylindrical coordinate grid division section of the present invention.

[0065] Figure 4 This is a schematic diagram of the two-dimensional flight airspace grid division of the present invention. Detailed Implementation

[0066] This invention provides a rapid trajectory planning and guidance method for rocket landing phase based on flight airspace grid division, comprising the following steps:

[0067] 1) Establish a landing column coordinate system and select the longitudinal section with a polar angle of zero as the reference section;

[0068] 2) Perform rectangular meshing on the reference profile selected in step 1) and store the position parameters of the mesh nodes;

[0069] 3) Establish a two-dimensional dynamic model of rocket landing in the reference profile selected in step 1);

[0070] 4) Based on the dynamic model established in step 3), construct a minimum fuel optimization model, and discretize and convexize it to obtain a convex optimization solution model for the standard landing trajectory planning problem sequence.

[0071] 5) Remove the initial velocity constraint from the standard optimization model established in step 4) to obtain the optimal ballistic solution model; substitute the grid node position parameters obtained in step 2) into the optimal ballistic solution model as initial position parameters for optimization and solution to obtain the absolute optimal ballistic parameters, and extract the absolute optimal ballistic initial velocity as the absolute optimal velocity parameter of the grid node, establish an interpolation table of grid node position parameters and absolute optimal velocity parameters, and bind it as parameters to the onboard computer;

[0072] 6) During flight, based on the rocket's current real-time position, the absolute optimal velocity is obtained by interpolating the rocket position parameters in step 5) using the absolute optimal velocity interpolation table. By performing trajectory integration, the rocket's position change is obtained. By performing cyclic interpolation and trajectory integration, the absolute optimal trajectory can be quickly reconstructed.

[0073] 7) Select a suitable point on the absolutely optimal trajectory obtained in step 6) as the target point; take the current position and current velocity of the rocket as the initial position and initial velocity, and the position and velocity of the target point as the final position and velocity, and substitute them into the sequential convex optimization solution model of the standard landing trajectory planning problem established in step 4). The local two-dimensional optimal trajectory of the rocket from the current point to the selected target point in the longitudinal flight plane can be obtained in the cylindrical coordinate system.

[0074] 8) Based on the two-dimensional trajectory obtained in step 7), the optimal trajectory is tracked using a fourth-order polynomial to obtain the magnitude and direction of the rocket's acceleration in the longitudinal plane of the cylindrical coordinate system; the remaining flight time of the rocket is estimated in the longitudinal plane; in the lateral direction of the cylindrical coordinate system, the tangential velocity component of the rocket is decelerated to zero in the remaining flight time using polynomial guidance to obtain the magnitude and direction of the tangential acceleration.

[0075] 9) Rotate the acceleration vector obtained in the cylindrical coordinate system to the fixed rectangular coordinate system at the landing point; calculate the rocket's three-axis attitude angle command based on the three-axis components of the acceleration command in the rectangular coordinate system; calculate the rocket's thrust magnitude command based on the magnitude of the acceleration command in the rectangular coordinate system.

[0076] 10) Determine if the current rocket flight altitude is lower than the preset altitude. If it is, end the online trajectory planning; if not, update the trajectory planning and guidance cycle, and return to step 5).

[0077] The relationship between the landing column coordinate system and the fixed landing point coordinate system described in step 1) is as follows: Figure 2 As shown; where x, y, and z are the three coordinate systems of the fixed coordinate system at the landing point, designed to coincide with the North-Sky-East coordinate system, therefore x, y, and z also represent the north, sky, and east directions respectively; ρ co For the radial direction in cylindrical coordinates, y co θ represents the longitudinal direction of the cylindrical coordinate system. co θ is the polar angle in cylindrical coordinates.

[0078] The reference profile mentioned in step 1) is a longitudinal plane with a polar angle of zero, and its specific spatial location is as follows: Figure 3 As shown, where The longitudinal plane formed is the reference section, and the cylindrical coordinate system can be regarded as being formed by rotating the reference section, with the rotation angle equal to θ. co If so, Then it can be assumed that the planar positions of the rockets are the same, if Therefore, we can assume that the planar motion of the rockets is the same. Furthermore, if... (or Therefore, the motion of the rocket relative to the landing point can be considered equivalent. Among them... It is the velocity component of the rocket's vertical longitudinal plane, also known as tangential velocity, and its effect is to rotate the longitudinal plane.

[0079] Step 2) describes the rectangular grid division as follows: Figure 4 As shown. The ramp angle is a fixed, inverted cone-shaped airspace at the landing point, designed to prevent the rocket from landing too low and colliding with ground obstacles. The area above the ramp angle is the flight airspace, and the area below it is the no-fly zone.

[0080] The two-dimensional dynamic model of the rocket established in step 3) is as follows:

[0081]

[0082] Where r(t)=[ρ co (t),y co (t)] T is the position vector on the longitudinal plane of the cylindrical coordinate system; It is the velocity vector on the longitudinal plane of the cylindrical coordinate system; It is the thrust vector in the longitudinal plane of the cylindrical coordinate system; g = [0, -g0] is the constant gravitational acceleration vector in the longitudinal plane of the cylindrical coordinate system; g0 is the magnitude of the gravitational acceleration at sea level; I sp ρ is the engine specific impulse; ρ is the atmospheric density; S D C is the aerodynamic reference area of ​​the rocket. D This represents the rocket drag coefficient.

[0083] The fuel-efficient optimization model for step 4) is:

[0084]

[0085] m(0)=m0, r(0)=r0

[0086]

[0087] Where, θ gs It is the slope angle, θ max This is the maximum engine tilt angle. The meanings of the other variables are the same as those described above, and all have been explained.

[0088] Next, we need to construct a convex optimization model to perform lossless convexization of the thrust. By adding the intermediate quantity Γ(t), we can construct the convex function constraint of the thrust as: ||T(t)||≤Γ(t), where T min ≤Γ(t)≤T max

[0089] Discretization yields the following two-dimensional sequential convex optimization model in the longitudinal plane of the cylindrical coordinate system:

[0090]

[0091] subject to r(n+1)=r(n)+v(n)Δt

[0092] v(n+1)=v(n)+((T(n)-0.5ρS D ·

[0093] C D ||v(n)||v(n)) / m(n)+g)Δt

[0094]

[0095] ||Γ(n+1)-Γ(n)||≤dT limit Δt n=1,2,3,...,N-1

[0096]

[0097] ||T(n)||≤Γ(n),T min ≤Γ(n)≤T max

[0098]

[0099] m(1)=m0, r(1)=r0, v(1)=v0

[0100]

[0101] Where n is the number of the discrete point, n = 1, 2, 3, ..., N-1, N is the number of discrete points, r(n) is the position vector of the nth discrete point, v(n) is the velocity vector of the nth discrete point, T(n) is the rocket engine thrust vector at the nth discrete point, Γ(n) is the introduced scalar relaxation factor, m(n) is the rocket mass at the nth discrete point, ρ is the atmospheric density, and S... D For the rocket's aerodynamic reference area, C D Let ||v(n)|| be the aerodynamic drag coefficient of the rocket, ||v(n)|| be the magnitude of the velocity vector at the nth discrete point, and Δt be the time interval of the rocket optimization model.

[0102] In the practical application of the aforementioned optimization model, in order to improve the accuracy of discrete optimization, the position recursion and velocity recursion of the discrete state equations adopt the trapezoidal integral formula, as follows:

[0103] r(n+1)=r(n)+0.5(v(n)+v(n+1))Δt+α r (n)

[0104] v(n+1)=v(n)+(0.5(T(n)-0.5ρS D C D ||v(n)||v(n)+T(n+1)-0.5ρS D C D ||v(n+1)||v(n+1)) / m(n)+g)Δt+α v (n),

[0105] In step 5), removing the initial velocity constraint yields the design form of the absolutely optimal ballistic solution model as follows:

[0106]

[0107] subject to r(n+1)=r(n)+0.5(v(n)+v(n+1))Δt+α r (n)

[0108] v(n+1)=v(n)+(0.5(T(n)-0.5ρS D C D ||v(n)||v(n)+T(n+1)-0.5ρS D C D ||v(n+1)||v(n+1)) / m(n))Δt+gΔt+α v (n)

[0109]

[0110] ||Γ(n+1)-Γ(n)||≤dT limit Δt n=1,2,3,...,N-1

[0111]

[0112] ||T(n)||≤Γ(n),T min ≤Γ(n)≤T max

[0113]

[0114] m(1)=m0,r(1)=r0

[0115]

[0116] Where, α r (n) represents the position recursion virtual control term; α v (n) is the velocity recursion virtual control term; α m (n) is the quality recursive virtual control term; k is the weight coefficient of the virtual control term; the meanings of the other variables are the same as those above, and have all been explained.

[0117] The specific steps for optimization are as follows:

[0118] The first step is to generate an initial solution using polynomial guidance, as follows:

[0119]

[0120]

[0121] The above equations are iteratively applied until the position vector reaches the landing point within the allowable deviation range, at which point the loop stops. This yields the trajectory versus time curve, velocity versus time curve, mass versus time curve, and total flight time t. total_init ;

[0122] Based on the selected number of discrete points N, the initial time interval Δt = t for the optimization model is obtained. total_init / (N-1);

[0123] Based on this time, time interpolation can be performed on the obtained trajectory-time curve, velocity-time curve, and mass-time curve to generate a set of initial solutions that satisfy thrust constraints and dynamic constraints.

[0124] The second step is to substitute the initial solution into the sequential convex optimization model and perform optimization to obtain a new set of optimized solutions [r,v,T,Γ,m,α]. r ,α v ,α m ].

[0125] The third step is to determine the deviation of virtual control item parameters;

[0126] if, That is, when all the values ​​of the virtual control terms solved in the (q+1)th iteration are less than the set virtual control term deviation, the iterative solution process ends immediately, and the solution of the (q+1)th iteration is the final solution;

[0127] Otherwise, proceed to step four;

[0128] The fourth step is to determine the deviation between the state parameters and control parameters;

[0129] if, If the changes in the optimized state parameters and control parameters in the (q+1)th iteration relative to the optimized state parameters and control parameters in the qth iteration are less than the set deviation range, then the optimization solution is considered to have converged, and the iterative solution process ends.

[0130] Otherwise, take the solution of the (q+1)th optimization iteration as the initial solution and return to the second step.

[0131] In step 5), the interpolation table is constructed by taking the feasible spatial grid nodes obtained in step 2). Substituting r0 as the initial position constraint into the absolute optimal ballistics optimization model in turn, the following results are obtained sequentially. and establish about The interpolation table.

[0132] In step 6), the method for restoring the absolutely optimal trajectory is as follows:

[0133]

[0134] ρ co (1)=ρ co_now ,y(1)=y now

[0135]

[0136] Among them, [x now ,y now ,z now [ρ] represents the position of the rocket in a fixed rectangular coordinate system at the landing point. co_now ,y now ,z co_now These represent the radial, longitudinal, and tangential positions in the rocket's column coordinate system, where the tangential position z... co_now =0, interp() is a linear interpolation function for the rocket position parameter-absolute optimal velocity interpolation table, the values ​​in parentheses are the target interpolation points, θ co T is the polar angle in cylindrical coordinates; integral This is the ballistic integration step size.

[0137] In step 7), the target point is selected as follows:

[0138] The closer the target point is to the rocket's current position, the fewer discrete points are needed for online trajectory optimization, resulting in faster computation. However, the cycle time for online trajectory optimization needs to be very small to ensure frequent trajectory updates. Conversely, the farther the target point is from the rocket's current position, the more discrete points are needed for online trajectory optimization, resulting in slower computation, but the cycle time can be appropriately increased. Therefore, the overall design of the target point height is as follows:

[0139] y target =y now -2T TP ·v y_now

[0140] Among them, y target Let y be the height of the target point. now v is the current height. y_now T represents the current falling speed. TP The trajectory planning period.

[0141] In step 7), the local two-dimensional optimal trajectory optimization model from the current position to the target position is as follows:

[0142]

[0143] subject to r(n+1)=r(n)+0.5(v(n)+v(n+1))Δt+α r (n)

[0144] v(n+1)=v(n)+(0.5(T(n)-0.5ρS D C D ||v(n)||v(n)+T(n+1)-0.5ρS D C D ||v(n+1)||v(n+1)) / m(n))Δt+gΔt+α v (n)

[0145]

[0146] T min ≤Γ(n)≤T max

[0147]

[0148] m(1)=m0,r(1)=r co0 ,v(1)=v co0 ,r(N)=r cof +r err v(N)=v cof +v err Where s1, s2, and s3 are the coefficient weights of the fuel consumption index, position deviation index, and speed deviation index, respectively; k is the coefficient weight of the virtual control term; r err To optimize the deviation between the trajectory terminal position and the target position; v err To optimize the deviation between the trajectory terminal velocity and the target velocity; r co0 =[ρ co0 [y0] represents the current position vector of the rocket on the longitudinal plane of the cylindrical coordinate system; r is the rocket's current velocity vector on the longitudinal plane of the cylindrical coordinate system. cof =[ρ cof ,y f ] is the position vector of the rocket target point on the longitudinal plane of the cylindrical coordinate system; Let be the velocity vector of the rocket target point on the longitudinal plane of the cylindrical coordinate system.

[0149] Its specific iterative solution method is the same as the iterative process of optimization solution in step 5).

[0150] In step 7), the weight design of performance metrics in the local two-dimensional optimal trajectory optimization model is as follows:

[0151]

[0152] s2=1000||r err || / ||r cof ||

[0153] s3 = 1000||verr || / ||v cof

[0154] k = 1000

[0155] s1, s2 and s3 are variable weight designs; The three-axis velocity vectors are fixed in a rectangular coordinate system at the landing point; These are the radial velocity, longitudinal velocity, and tangential velocity in the rocket's cylindrical coordinate system; V err_now This is the magnitude of the deviation between the velocity vector in the longitudinal plane of the rocket's current cylindrical coordinate system and the optimal velocity vector obtained by interpolation at this point; the meanings of the other variables are the same as above, and all have been explained.

[0156] In step 8), the specific steps of the polynomial guidance tracking method are as follows:

[0157] The first step is to select the guidance target point. The height of the guidance target point is set to y. co_next =y co_now -2T guid ·v y_now Based on the optimized trajectory from step 7), the height is interpolated to obtain ρ. co_next , These are the radial position, longitudinal velocity, radial velocity, longitudinal acceleration, and radial acceleration of the guided target point, respectively.

[0158] The second step is to obtain the acceleration command.

[0159]

[0160]

[0161] in, These are longitudinal acceleration requirements, radial acceleration requirements, and tangential acceleration requirements.

[0162] The third step is to convert the acceleration requirements in the cylindrical coordinate system to the landing and fixed rectangular coordinate system;

[0163]

[0164] The third step is to generate guidance commands;

[0165] γ CX =0

[0166]

[0167] Where, γ CX , ψ CX , These are the roll program angle, yaw program angle, and pitch program angle, respectively; T cmdThis is the engine thrust command.

[0168] Fourth step: Update the guidance cycle and return to step one;

[0169] The contents not described in detail in this specification are common knowledge to those skilled in the art.

Claims

1. A rapid trajectory planning and guidance method for rocket landing phase based on flight airspace grid division. Its features include the following steps: 1) Establish a landing column coordinate system and select the longitudinal section with a polar angle of zero as the reference section; 2) Perform rectangular meshing on the reference profile selected in step 1) and store the position parameters of the mesh nodes; 3) Establish a two-dimensional dynamic model of rocket landing in the reference profile selected in step 1); 4) Based on the dynamic model established in step 3), construct a minimum fuel optimization model, and discretize and convexize it to obtain a convex optimization solution model for the standard landing trajectory planning problem sequence. 5) Remove the initial velocity constraint from the standard optimization model established in step 4) to obtain the optimal ballistic solution model; substitute the grid node position parameters obtained in step 2) into the optimal ballistic solution model as initial position parameters for optimization and solution to obtain the absolute optimal ballistic parameters, and extract the absolute optimal ballistic initial velocity as the absolute optimal velocity parameter of the grid node, establish an interpolation table of grid node position parameters and absolute optimal velocity parameters, and bind it as parameters to the onboard computer; 6) During flight, based on the rocket's current real-time position, the absolute optimal velocity is obtained by interpolating the rocket position parameters in step 5) using the absolute optimal velocity interpolation table. By performing trajectory integration, the rocket's position change is obtained. By performing cyclic interpolation and trajectory integration, the absolute optimal trajectory can be quickly reconstructed. 7) Select a suitable point on the absolutely optimal trajectory obtained in step 6) as the target point; take the current position and current velocity of the rocket as the initial position and initial velocity, and the position and velocity of the target point as the final position and velocity, and substitute them into the sequential convex optimization solution model of the standard landing trajectory planning problem established in step 4). The local two-dimensional optimal trajectory of the rocket from the current point to the selected target point in the longitudinal flight plane can be obtained in the cylindrical coordinate system. 8) Based on the two-dimensional trajectory obtained in step 7), the optimal trajectory is tracked using a fourth-order polynomial to obtain the magnitude and direction of the rocket's acceleration in the longitudinal plane of the cylindrical coordinate system; the remaining flight time of the rocket is estimated in the longitudinal plane; in the lateral direction of the cylindrical coordinate system, the tangential velocity component of the rocket is decelerated to zero in the remaining flight time using polynomial guidance to obtain the magnitude and direction of the tangential acceleration. 9) Rotate the acceleration vector obtained in the cylindrical coordinate system to the fixed rectangular coordinate system at the landing point; calculate the rocket's three-axis attitude angle command based on the three-axis components of the acceleration command in the rectangular coordinate system; calculate the rocket's thrust magnitude command based on the magnitude of the acceleration command in the rectangular coordinate system. 10) Determine if the current rocket flight altitude is lower than the preset altitude. If it is, end the online trajectory planning; if not, update the trajectory planning and guidance cycle, and return to step 5).

2. The rapid trajectory planning and guidance method for rocket landing phase based on flight airspace grid division according to claim 1, characterized in that... In step 5), the form of the absolutely optimal ballistic solution model is as follows: Where n is the number of the discrete point, n = 1, 2, 3, ..., N-1, N is the number of discrete points, r(n) is the position vector of the nth discrete point, v(n) is the velocity vector of the nth discrete point, T(n) is the rocket engine thrust vector at the nth discrete point, Γ(n) is the introduced scalar relaxation factor, m(n) is the rocket mass at the nth discrete point, ρ is the atmospheric density, and S... D For the rocket's aerodynamic reference area, C D Let α be the aerodynamic drag coefficient of the rocket, ||v(n)|| be the magnitude of the velocity vector at the nth discrete point, Δt be the time interval of the rocket optimization model, and α be the velocity vector at the nth discrete point. r (n) represents the position recursion virtual control term, α v (n) is the velocity recursion virtual control term, α m (n) is a quality recursive virtual control term.

3. The rapid trajectory planning and guidance method for rocket landing phase based on flight airspace grid division according to claim 2, characterized in that... In the optimization model, the position and velocity recursion of the discrete state equations adopt the trapezoidal integral formula, as follows:

4. The rapid trajectory planning and guidance method for rocket landing phase based on flight airspace grid division according to claim 1, characterized in that... In step 6), the method for restoring the absolutely optimal ballistic trajectory is as follows: r co (1)=r co_now ,y(1)=y now [v ρco (i),v y (i)]=interp(ρ co (i),y(i)) Among them, [x now ,y now ,z now [ρ] represents the position of the rocket in a fixed rectangular coordinate system at the landing point. co_now ,y now ,z co_now These represent the radial, longitudinal, and tangential positions in the rocket's column coordinate system, where the tangential position z... co_now =0, interp() is a linear interpolation function for the rocket position parameter-absolute optimal velocity interpolation table, the values ​​in parentheses are the target interpolation points, θ co T is the polar angle in cylindrical coordinates; integral This is the ballistic integration step size.

5. The rapid trajectory planning and guidance method for rocket landing phase based on flight airspace grid division according to claim 1, characterized in that... In step 7), the target point is selected in the following way: The closer the target point is to the rocket's current position, the fewer discrete points are needed for online trajectory optimization, resulting in faster computation. However, the cycle time for online trajectory optimization needs to be very small to ensure frequent trajectory updates. Conversely, the farther the target point is from the rocket's current position, the more discrete points are needed for online trajectory optimization, resulting in slower computation, but the cycle time can be appropriately increased. Therefore, the overall design of the target point height is as follows: y target =y now -2T TP ·v y_now Among them, y target Let y be the height of the target point. now v is the current height. y_now T represents the current falling speed. TP The trajectory planning period.

6. The rapid trajectory planning and guidance method for rocket landing phase based on flight airspace grid division according to claim 1, characterized in that... In step 7), the local two-dimensional optimal trajectory optimization model from the current position to the target position is as follows: Where s1, s2, and s3 are the coefficient weights of the fuel consumption index, position deviation index, and speed deviation index, respectively; k is the coefficient weight of the virtual control term; r err To optimize the deviation between the trajectory terminal position and the target position; v err To optimize the deviation between the trajectory terminal velocity and the target velocity; r co0 =[ρ co0 [y0] represents the current position vector of the rocket on the longitudinal plane of the cylindrical coordinate system; r is the rocket's current velocity vector on the longitudinal plane of the cylindrical coordinate system. cof =[ρ cof ,y f ] is the position vector of the rocket target point on the longitudinal plane of the cylindrical coordinate system; Let be the velocity vector of the rocket target point on the longitudinal plane of the cylindrical coordinate system.

7. The rapid trajectory planning and guidance method for rocket landing phase based on flight airspace grid division according to claim 6, characterized in that... In the local two-dimensional optimal trajectory optimization model in step 7), the weight design of the performance index is as follows: s2＝1000||r err || / ||r cof || s3=1000||v err || / ||v cof || k＝1000 s1, s2, and s3 are variable weight designs; [v xnow ,v ynow ,v znow [v] represents the three-axis velocity vector of the landing point in a fixed Cartesian coordinate system; ρco_now ,v ynow ,v zco_now These represent the radial velocity, longitudinal velocity, and tangential velocity in the rocket's cylindrical coordinate system, respectively; V err_now This is the magnitude of the deviation between the velocity vector in the longitudinal plane of the rocket's current cylindrical coordinate system and the optimal velocity vector obtained by interpolation at this point; the meanings of the other variables are the same as above, and all have been explained.

Images