CN111196382A - Convergence-guaranteed real-time trajectory planning method for powered descent of rockets - Google Patents

Abstract

The invention discloses a real-time trajectory planning method for ensuring convergence of a rocket power descent segment, and belongs to the field of rocket recovery guidance. The implementation method of the invention comprises the following steps: considering the nonlinear aerodynamic drag and the constraint conditions of the magnitude and direction of the thrust, establishing a landing optimal control problem model of the rocket power descent segment with the magnitude and direction of the thrust as control quantities; the problem is simplified by problem dimension reduction, problem decomposition and advanced planning of certain state variables, and the original problem is converted into an unconstrained convex optimization problem and a second-order cone planning problem by combining a proper convex method; and finally, solving a convex optimization problem and a second-order cone planning problem to obtain a change strategy of the flight trajectory and the magnitude and direction of the thrust. The method fully utilizes a reliable and efficient convex optimization algorithm, the calculation time is about 15-30 milliseconds on a common desktop computer with the main frequency of 3.6 gigahertz, the obtained solution is close to the optimal fuel, the method is converged to a certain extent, and the reliable real-time planning of the landing track of the rocket power descent segment can be realized.

Description

Real-time trajectory planning method for rocket power descent segment capable of guaranteeing convergence

Technical Field

The invention belongs to the field of rocket recovery guidance, relates to a trajectory planning method for a rocket recovery power descent section, and particularly relates to a trajectory planning method for guaranteeing convergence and realizing real-time performance based on convex optimization.

Background

As early as the 60's of the last century, humans achieved lunar landing. For realizing the strong lifting, the soft landing technology of the manned lander on the moon is crucial. In the next decades, technological breakthroughs were made in landing landers on the rarefied stars in the atmosphere. However, landing a rocket on the earth remains challenging. In recent years, reusable rocket technology has attracted considerable attention worldwide. The technology can greatly reduce the rocket launching cost and greatly shorten the launching period. In order to achieve precise landing of the rocket, precise guidance of the power descent segment plays an extremely important role.

The landing trajectory planning of the dynamic descent segment needs to solve an optimal control problem. In order to improve the anti-interference capability and the maneuverability of the rocket in the power descent stage, the non-convex optimal control problem needs to be solved in real time so as to realize the accurate landing of the rocket. Solving this non-convex problem directly by existing methods (non-linear programming or targeting based on optimal control theory) is often difficult to apply in practical engineering, because such methods cannot guarantee their reliability and the solving efficiency is low. At present, the solution of the landing trajectory planning problem of rocket power descent mostly depends on direct method solution, such as nonlinear programming and sequence convex optimization. However, these algorithms cannot guarantee the reliability of the solution, and are currently difficult to apply in engineering.

Disclosure of Invention

Aiming at the defects of reliability and efficiency of the traditional rocket power descent trajectory planning method, the invention discloses a rocket power descent segment real-time trajectory planning method for ensuring convergence, which aims to solve the technical problems that: the problem is simplified through problem dimension reduction, problem decomposition and advanced planning of certain state variables, and then a proper convex method is combined, so that only one unconstrained convex optimization problem and at most two second-order cone planning problems need to be solved, namely, a reliable and efficient convex optimization algorithm is fully utilized to solve the optimal control problem which is originally difficult to solve, and the real-time trajectory planning efficiency and robustness of the rocket power descent segment are improved.

The purpose of the invention is realized by the following technical scheme.

The invention discloses a real-time trajectory planning method for ensuring convergence of a rocket power descent segment, which considers the nonlinear aerodynamic resistance and the constraint conditions of the magnitude and direction of thrust and establishes a landing optimal control problem model of the rocket power descent segment by taking the magnitude and direction of the thrust as control variables; the problem is simplified by problem dimension reduction, problem decomposition and advanced planning of certain state variables, and the original problem is converted into an unconstrained convex optimization problem and a second-order cone planning problem by combining a proper convex method; and finally, solving a convex optimization problem and a second-order cone planning problem to obtain a change strategy of the flight trajectory and the magnitude and direction of the thrust. The method fully utilizes a reliable and efficient convex optimization algorithm, the calculation time is about 15-30 milliseconds on a common desktop computer with the main frequency of 3.6 gigahertz, the obtained solution is close to the optimal fuel, the method is converged to a certain extent, and the reliable real-time planning of the landing track of the rocket power descent segment can be realized.

The invention discloses a rocket power descent segment real-time trajectory planning method capable of ensuring convergence, which comprises the following steps:

the method comprises the following steps: and performing dynamic modeling and dimensional normalization on the rocket power descending process to establish a three-dimensional dimensionless dynamic equation. And the constraint required by power descent flight is introduced to establish the optimal control problem of the power descent landing with optimal fuel.

The specific implementation method of the step one is as follows:

carrying out dynamic modeling on the rocket power descent flight, carrying out dimensional normalization, and expressing a dimensionless dynamic equation of the rocket power descent flight as

Wherein x, y and h are rocket positions, h is a height direction, the pointing direction is positive, x is a direction from a projection point of an initial position on a horizontal plane to a landing point, and y, h and x form a right-hand rule; v is the rocket speed; theta is a speed inclination angle, namely an included angle between the projection of the speed vector on an Oxh plane and the x axis, and the projection of the speed is positive above the x axis;

the track yaw angle is the included angle between the velocity vector and the Oxh plane, and the same side of the velocity vector and the y axis on the Oxh plane is positive; m is the rocket mass; e and sigma are used for expressing the direction of the thrust, wherein the e expresses the included angle between the thrust direction and the opposite direction of the speed; g is the acceleration of gravity at height h, g₀Corresponding to a gravitational acceleration of height 0; t represents the magnitude of the thrust; d represents the aerodynamic resistance; i is_spIs the specific impulse of the rocket engine. In the formula (1), variables other than the angle variable are subjected to dimensional normalization, and the position variables x, y, and h are normalized by the initial height h₀Initial velocity V for velocity V₀Initial mass m for mass m₀Time and contrast I_spBy using h₀/V₀Acceleration of gravity g

Thrust force T is used

To perform dimensional normalization respectively. Wherein the dimensionless resistance is expressed as

Where ρ is the dimensionless air density, as a function of altitudeTransformation, S_refIs a dimensionless reference area of the rocket, C_DIs the drag coefficient.

Introducing the constraints required for dynamic descent flight. First, the magnitude of thrust is constrained as follows

T_min≤T≤T_max(3)

Wherein T is_min>0 and T_maxIs the minimum and maximum allowable thrust magnitude. Further, the thrust direction is constrained as follows

0≤∈≤∈_max(4)

Wherein e_maxE [0, π/2) is the maximum allowable thrust direction and velocity counter direction angle. Finally, to achieve accurate landing, the following terminal constraints need to be satisfied

x(t_f)＝0 (5)

y(t_f)＝0 (6)

h(t_f)＝0 (7)

V(t_f)≤V_f(8)

θ(t_f)＝-π/2 (9)

Wherein V_fIs a small safe landing speed. Constraints (5) - (7) ensure that the rocket lands on the designated landing site, and constraints (8) - (10) ensure that the landing speed is less than V_fAnd is perpendicular to the landing site ground.

The optimization objective is to minimize fuel economy, thus creating an optimal control problem for power-down landing as follows

s.t. formula (1), (3) - (10)

In the optimal control problem P0, the kinetic equation is highly nonlinear, with time of flight being an optimization variable. Obviously, problem P0 is a non-convex problem.

Step two: since the time of flight is free, the independent variables of the kinetic equation are converted into altitude, and the terminal constraints and the objective function are correspondingly converted.

The second step is realized by the following concrete method:

the time of flight of the rocket during the power descent is unknown, but the initial and terminal altitudes are known, and during the power descent the altitude monotonically decreases over time, i.e., the rocket does not fly upwards. Thus converting the independent variable of formula (1) into a height

Wherein k is_D＝0.5ρS_refC_DThe superscript (') indicates the derivation of the height h, and h is the normalized height from 1 to 0. The new arguments enable the terminal constraints (5) - (10) to be translated into

x(h_f)＝0 (13)

y(h_f)＝0 (14)

V(h_f)≤V_f(15)

θ(h_f)＝-π/2 (16)

Conversion of the objective function (11) into

The time-independent kinetic equation (1) is converted into the height-independent kinetic equation (12), and the corresponding terminal constraints and objective functions are also converted into the height-independent forms (13) - (18).

Step three: reducing the non-linearity of the kinetic equation. Part of nonlinearity in a kinetic equation is transferred into a constraint, then the nonlinearity is further reduced by reducing the dimension of the kinetic equation, two problems with smaller dimension are further established, namely a longitudinal problem and a lateral problem, and finally the lateral problem is simplified into a polynomial coefficient solving problem.

The third concrete implementation method comprises the following steps:

the primary nonlinearity of the original problem P0 is in the kinetic equation. To reduce the nonlinearity of the kinetic equation, three new control variables are defined as follows

u₁:＝T cos∈/m (19)

u₂:＝T sin∈cosσ/m (20)

u₃:＝T sin∈sinσ/m (21)

And a new state variable

ω:＝ln m (22)

The control variable nonlinearity in the third to fifth equations in the kinetic equation (12) is eliminated according to the defined variables, and the mass-cost equation is converted into

The thrust magnitude and direction constraints (3) - (4) are converted into

In addition, the optimization objective (18) is translated into

A reduction in the non-linearity of the kinetic equation results in an increase in the non-linearity of the constraint and objective function.

The non-linearity of the kinetic equation is further reduced by reducing the kinetic equation to dimensionality. The kinetic equation is divided into three parts, namely the longitudinal kinetic equation with respect to the state quantities x, V and theta

With respect to states y and

equation of lateral dynamics

And mass consumption equation (23).

The purpose of the kinetic equation decomposition is to transform the original problem P0 into a less-dimensional longitudinal optimal control problem and a lateral planning problem. The longitudinal dynamics equation (27) is included in the longitudinal problem. Since only u₁And u₂The constraints (24) - (25) on the magnitude and direction of thrust existing in the longitudinal dynamical equations, and therefore, in the longitudinal problem, need to be removed₃Converted into the following form

-u₁tan(∈_max-△_∈)≤u₂≤u₁tan(∈_max-△_∈) (30)

△ therein_TAnd △_∈Is the magnitude and direction of thrust reserved for lateral movement. When u is obtained by a lateral problem₃After u, u₁，u₂And u₃The constraints (24) - (25). △ that need to satisfy the original thrust magnitude and direction_TAnd △_∈Is calculated as follows

△ therein_∈,maxIs an avoidance △_∈An excessive threshold value, parameter κ used to reflect the relative relationship between longitudinal and lateral maneuvers, is

Where { x_f,soft,y_f,softAnd the terminal position of a soft landing point is obtained from the soft landing track obtained in the fourth step.

The longitudinal optimal control problem is written as

s.t. x′＝cotθ (35)

-u₁tan(∈_max-△_∈)≤u₂≤u₁tan(∈_max-△_∈) (39)

x(h_f)＝0,V(h_f)≤V_f,θ(h_f)＝-π/2 (40)

Problem(s)

Dependent on the unknown state quantities m and

the m acquisition method is realized in step four.

The planning problem of establishing a track yaw profile, namely the lateral problem, is as follows

Problem P_latThe objective is to plan a feasible track yaw profile to meet the terminal constraints of position y. Clearly, there are many possible track yaw profiles. The track yaw profile must satisfy certain properties. First, a feasible track yaw profile must meet

And

second, the corresponding state y should satisfy y (h)₀)＝y(h_f) 0. According to the Rollo theorem, by y (h)₀)＝y(h_f) When the value is 0, a point h exists_i∈[h_f,h₀]So that y' (h)_i) 0, according to formula (41),

by

It is known that there is a point h_e∈[h_f,h_i]So that

Thus, a feasible flight pathYaw angle profile satisfies

And

dividing the flight path yaw angle into two sections for planning, namely h belongs to [ h ∈ [_e,h₀]And h e [ h ∈_f,h_e]. For the convenience of calculation, will

Consider a two-stage Bernstein polynomial, the first stage being a fourth order polynomial and the second stage being a second order polynomial. Thus, it is possible to provide

Is written as

Wherein the Bernstein coefficient { ζ_1,i}_i＝0,…,4And { ζ_2,i}_i＝0,…,2Needs to be solved. The first stage polynomial has five coefficients, but only two conditions

And

therefore, in order to solve all the coefficients, three additional conditions need to be included,

and y (h)_f) 0. The Bernstein coefficient of the first-segment polynomial is expressed as

The Bernstein coefficient of the second-segment polynomial may pass through the condition

And

to obtain

ζ_2,0＝0 (49)

The way these coefficients are solved is based on the properties of Bernstein polynomials. And (6) substituting all coefficients into an equation (43) to obtain a track yaw angle profile.

Part of the non-linearity of the kinetic equation is transferred to constraints (24) - (25) and then the non-linearity is further reduced by reducing the dimensionality of the kinetic equation, thereby creating two smaller dimensional problems, namely longitudinal problems

And lateral problems P_lat(θ), last lateral problem P_lat(theta) is converted into a Bernstein polynomial coefficient solving problem, the coefficient passing equation(44) - (51) obtaining.

Step four: further simplifying the longitudinal optimal control problem. The speed inclination angle is separated from the longitudinal kinetic equation, and the speed inclination angle is optimized independently, so that the nonlinearity of the longitudinal kinetic equation is reduced.

The fourth concrete implementation method comprises the following steps:

longitudinal optimal control problem

The kinetic equations in (1) are still highly non-linear and need further simplification. The non-linearity of the longitudinal dynamical equation is mainly due to the trigonometric function term of the velocity tilt angle θ. Thus solving the speed tilt angle from the longitudinal direction

And (4) separating, and independently optimizing the speed inclination angle. Thus, longitudinal problems

Is simplified into

-u₁tan(∈_max-△_∈)≤u₂≤u₁tan(∈_max-△_∈) (56)

V(h_f)≤V_f(57)

With respect to states x and thetaThe kinetic equation is used to plan the velocity-tilt profile and is therefore not a problem

In (1). Compared with the problem

Problem(s)

The non-linearity is less but still a non-convex problem. Step five will be right to problem

And (4) carrying out convex formation.

It is easy to obtain a feasible velocity tilt profile such that the state x satisfies the terminal constraints. However, a feasible but unreasonable velocity tilt angle profile may cause problems

It is not feasible. Problem(s)

Is mainly due to the thrust direction constraint (56) not being satisfied. | u₂Too large | easily results in the constraint (56) not being satisfied. Therefore, | u that needs to be generated when planning the velocity tilt angle₂| is as small as possible. Therefore, the following optimal control problem is obtained

s.t.x′＝cotθ (59)

x(h_f)＝0,θ(h_f)＝-π/2 (60)

Problem(s)

Dependent on V and

due to the problems

Has not been solved, so the velocity V is unknown, and furthermore due to the problem P_lat(theta) depends on theta, so track yaw angle

Is also unknown. But V and

present only in the objective function, hence V and

is sufficient to solve the problem

First, an approximate velocity profile is planned

Considering a soft landing problem, the thrust magnitude profile of the soft landing is T_min-T_maxThe bang-bang structure of (1) has a conversion time of t_sAnd the thrust direction is always the opposite direction of speed, i.e., ∈ ═ 0. Transition time t_sIs the variable to be sought. Using the thrust profile to integrate the kinetic equation (1) until V (t) is less than or equal to V_fAnd (5) stopping. Note that this time is t_fAnd the height at that moment is h_f(t_s). If h is_f(t_s) When t is equal to 0, then t is_sIs a solution to the soft landing problem. The soft landing problem is written as

P_soft:find t_s

s.t. formula (1), T is T_min-T_maxForm e is 0

h(t_f)＝0 (61)

Problem P_softIs a one-dimensional root finding problem. Using the velocity profile of the soft landing trajectory as an approximate velocity profile

Further, the profile described by the following formula is used as an approximate track yaw profile

Is called as in the above formula "

Initialization ".

When a problem arises

V in (1) and

quilt

And

after the substitution, problems

It is still a non-convex problem. In order to improve the efficiency of solving the problem, it must be transformed into a convex optimization problem. First, the term cot θ is expressed as a fifth-order Bernstein polynomial as follows

Therefore, equation (64) holds

Substituting the above formula into the problem

In the objective function of

The Bernstein polynomial B in the objective function (65) has six coefficients { ζ }_i}_i＝0,…,5Where four coefficients are represented analytically. According to the problems

The following equations (66) to (68) are clearly true

B(h₀)＝cotθ₀(66)

B(h_f)＝cotθ_f(67)

B′(h₀)＝θ₀′csc²θ₀(68)

Wherein theta is₀:＝θ(h₀) And theta_f:＝θ(h_f) Is a known quantity, θ₀′:＝θ′(h₀) Is the initial derivative of the velocity tilt angle. According to the properties of the polynomials (66) to (68) and Bernstein

ζ₀＝0 (69)

ζ₅＝cotθ₀(71)

Furthermore, the kinetic equation for state x is rewritten as

According to the integral property of Bernstein polynomial

ζ₃＝6x₀-[(ζ₀+ζ₄+ζ₅)+(ζ₁+ζ₂)](73)

Wherein ζ₀，ζ₄And ζ₅Obtained from formulae (69) to (71). Problem(s)

Is used to obtain ζ₀，ζ₃，ζ₄And ζ₅. When its value is substituted into the objective function (65), the objective function depends only on ζ₁And ζ₂. Therefore, problems arise

Unconstrained optimization problem expressed as

The function F is actually a function of the variable ζ₁And ζ₂A convex function of (a). Thus, problems arise

The method is an unconstrained convex optimization problem and can be quickly and reliably solved by a quasi-Newton method. When ζ is₁And ζ₂After obtaining, the Bernstein coefficient { ζ_i}_i＝0,…,5A velocity gradient profile is obtained in place of equation (63).

In this step, the longitudinal optimal control problem is further simplified

Separating the velocity tilt angle from the longitudinal dynamics equation, reducing the non-linearity of the longitudinal optimal control problem to form a problem

For the separated velocity tilt angle, a problem is created

And then the problem is converted into an unconstrained convex optimization problem

Step five: and the simplified longitudinal optimal control problem is emphasized. Firstly, processing a nonlinear kinetic equation and a convex non-convex constraint condition, secondly, carrying out equivalent transformation on an objective function, and finally obtaining a convex longitudinal optimal control problem.

The concrete implementation method of the step five is as follows:

in the simplified longitudinal direction

In (2), the kinetic equation (53) is non-linear, the constraints (54) - (55) are non-convex, and the objective function (52) is also non-convex.

First, nonlinear kinetic equations (53) and equation constraints (54) are converted to linear. Definition V²As new state variables as follows

Thus, the kinetic equation (53) and the equality constraint (54) are converted to

Formulae (76) to (77) are all linear.

Non-convex constraints are then processed (55). There are two approaches to the convex representation of the constraint (55), the first of which is described first. In the constraint (55), the inequality determined by the second inequality sign is a second order cone constraint and is a convex constraint. Due to | u₂Is far less than | u₁I, so the constraint determined by the first inequality sign is approximated by T_min/m≤u₁. Thus, the constraint (55) is approximated as

The second processing method is a further simplification based on equation (78). The second order cone constraint is also reduced to a linear constraint. Thus, the non-convex constraint (55) is approximated as a linear constraint

And carrying out equivalent transformation on the non-convex objective function. The original objective function is replaced by a linear objective function

The objective functions (80) and (52) have an approximate optimization effect, which is the same even under certain conditions.

When the dynamic equation and the constraint are convexly processed and the objective function is converted, the original longitudinal problem

Is converted into a convex optimal control problem

-u₁tan(∈_max-△_∈)≤u₂≤u₁tan(∈_max-△_∈) (85)

In this step, the nonlinear kinetic equation (53) and the non-convex constraints (54) and (55) are embossed into a linear kinetic equation (76), a linear equation constraint (77), and either a convex constraint (78) or a linear constraint (79). Original longitudinal optimal control problem

Is embossed into

Step six: and sequentially solving the lateral problem, the speed inclination angle optimization problem and the longitudinal problem obtained in the third step to the fifth step. And then judging whether the track yaw angle profile needs to be improved and the quality profile needs to be updated according to the obtained solution, and further judging whether the longitudinal problem needs to be solved again. If not, outputting the obtained track; and if necessary, solving the longitudinal problem again, and updating the corresponding variable to obtain the track. The track is the track of the rocket power descending section planned in real time.

The sixth specific implementation method comprises the following steps:

in steps one through five, the problem P0 is converted into a longitudinal problem

And lateral problems P_lat(theta). The longitudinal problem is then decomposed into a problem for planning the velocity-tilt angle profile

And problems with

To improve the solution efficiency, question

Convex optimization problem being convex to be unconstrained

Problem(s)

Is embossed into a problem

In order to obtain the final trajectory, the problem needs to be solved in a certain order. First, solve the soft landing problem P_softAnd execute "

Initialization "get

And

then based on

And

solving a velocity dip angle planning problem

θ and x are obtained. Then substituting theta into the flight path yaw angle planning problem P_lat(theta) obtaining

And y. Finally, the sum of the values of theta,

and

substitution problem

To obtain V, u₁And u₂. Through the solving process, the state quantities x, y, V, theta,

And a control quantity u₁、u₂. U is calculated by the following formula (87)₃

Due to u₃Is not included in the problem

And △, and_∈cannot completely guarantee u₁，u₂And u₃The original thrust direction constraint (25) is satisfied. Therefore, the e is calculated by the formula (88), and whether the constraint (25) is satisfied is detected.

If e exceeds the boundary e_maxNeed to be aligned with

The profile is improved, so that the original thrust direction constraint is met. First, it is necessary to find the height point h at which ∈ reaches the maximum value_sep. Then for h e [ h ∈ [_sep,h₀]Segment, make ∈ equal to ∈_maxAnd calculating u inversely by the following formula₃

Wherein

Symbol sum u₃The same is true. Improved track yaw profile

By integration

The kinetic equation is obtained, i.e.

For h e [ h [ [ h ]_f,h_sep]Section, according to the method for planning track yaw angle in step three

Is a two-stage Bernstein polynomial

Wherein the Bernstein coefficient can be obtained by the method in step three. While improving track yaw profile

After being obtained, the corresponding

Calculated by equation (92).

In addition, when solving the problem

Mass is an approximate mass profile with a soft landing

Instead. When u is obtained₁，u₂And u₃Then, a more accurate quality profile is obtained

If the track yaw needs to be improved, in equation (93)

And u₃Are respectively as

And

based on the above two steps, two situations can occur. Firstly, the track yaw angle needs to be improved, and secondly, the updated quality profile m^dAnd approximate mass profile

With a large difference, i.e.

Wherein delta_mIs judgment m^dAnd

threshold value of the gap size. If neither of the two conditions occurs, the trajectory is found to be

Wherein m is m^dInstead. If one of the two situations occurs, the track yaw profile and the quality profile need to be updated. Then solve the problem again

And calculating new u₃. Then, whether the element belongs to the boundary element is judged_max. If the boundary is exceeded, the task is regarded as unreachable; if the boundary is not exceeded, the trajectory planning is successful,a new mass profile is calculated again using (93), the trajectory being found as

The method also comprises the seventh step: through the problem dimension reduction in the third step, the problem decomposition in the fourth step and the problem camouflaging in the fifth step, the trajectory of the rocket power descent section can be planned in the sixth step, namely, the rocket power descent section guidance is carried out through the rocket power descent section trajectory planned in real time in the sixth step, and the real-time trajectory planning efficiency and robustness of the rocket power descent section are improved.

Has the advantages that:

1. the invention discloses a real-time trajectory planning method for a rocket power descent segment, which simplifies the problem by problem dimension reduction, problem decomposition and advanced planning of certain variables, then uses a convexity method to convexity the simplified problem, namely converts the optimal control problem of the original non-convex rocket power descent segment into a series of problems which can be efficiently and reliably solved, and solves an unconstrained convex optimization and at most two second-order cone planning problems to ensure that the rocket power descent landing trajectory meeting the constraint can be obtained in a convergent manner.

2. The real-time trajectory planning method for the rocket power descent segment for ensuring convergence disclosed by the invention can obtain a rocket power descent landing trajectory with fuel close to the optimal value by optimizing the fuel consumption, so that the carrying capacity of a rocket can be improved.

Drawings

FIG. 1 is a diagram of the processing steps and relationships thereof for the problems involved in the method for real-time trajectory planning for a rocket powerdown leg to ensure convergence of the present invention;

FIG. 2 is a schematic diagram of a coordinate system and the definition of the velocity tilt angle and the track yaw angle in step one of the present invention;

FIG. 3 is a schematic diagram of a possible track yaw angle in step three of the present invention;

FIG. 4 is a graph of an objective function and a contour plot of the velocity-tilt optimization problem in step four of the present invention;

FIG. 5 is a schematic diagram of the feasible regions determined by the original thrust magnitude and direction constraint and the approximate thrust magnitude and direction constraint in step five of the present invention;

FIG. 6 is a pseudo code diagram of the rocket powerdown leg real-time trajectory planning method of the present invention ensuring convergence;

FIG. 7 is a variation of the velocity bank angle and track yaw angle profiles in an example of the invention;

FIG. 8 is a graph of the magnitude of the components of thrust in the direction of velocity, the direction of the velocity pitch angle, and the direction of the track yaw angle for an example of the present invention;

FIG. 9 is a graph showing the specific thrust magnitude and the change in the angle between the thrust direction and the reverse direction of the velocity in an example of the present invention;

FIG. 10 is a graph of velocity change, and changes in thrust acceleration and aerodynamic drag acceleration for an example of the present invention;

FIG. 11 is a graph of the initial planned track yaw angle and the improved track yaw angle, the change in the included angle between the thrust direction and the opposite direction of the speed, and the change in the mass profile difference between the first and second steps, according to an embodiment of the present invention;

FIG. 12 is a three-dimensional trajectory and thrust vector diagram of an example of the present invention.

Detailed Description

For a better understanding of the objects and advantages of the present invention, reference should be made to the following detailed description taken in conjunction with the accompanying drawings and examples.

In order to verify the feasibility of the method, a rocket power descent landing task is selected for verification. Initial mass of rocket is m₀26000kg, and the engine specific impulse is I_sp270s, coefficient of aerodynamic drag C_D1.8, reference area S_ref＝9m²Maximum thrust of the engine is T_max310 kN. The adjustable range of the thrust is 0.5T_maxTo T_max. The maximum allowable included angle between the thrust direction and the reverse direction of the speed is epsilon_max＝15°,△_∈Critical value △_∈,max5 deg. is equal to. In addition, a terminal speed limit is set to V_f2 m/s. Parameter h for planning track yaw angle_e300 m. Setting delta_m＝0.1％m₀. The initial state of the tasks in the example is shown in table 1.

TABLE 1 initial states of tasks in the example

As shown in FIG. 1, the real-time trajectory planning method for ensuring the convergence of the rocket power descent segment disclosed by the invention comprises the following steps:

the method comprises the following steps: and performing dynamic modeling and dimensional normalization on the rocket power descending process to establish a three-dimensional dimensionless dynamic equation. And the constraint required by power descent flight is introduced to establish the optimal control problem of the power descent landing with optimal fuel.

Carrying out dynamic modeling on the rocket power descent flight, carrying out dimensional normalization, and expressing a dimensionless dynamic equation of the rocket power descent flight as

Wherein x, y and h are rocket positions, h is a height direction, the pointing direction is positive, x is a direction from a projection point of an initial position on a horizontal plane to a landing point, y, h and x form a right-hand rule, and a coordinate system is defined as shown in fig. 2; v is the rocket speed; θ is the velocity tilt angle, as shown in fig. 2, i.e. the angle between the projection of the velocity vector on the Oxh plane and the x-axis, and the projection of the velocity is positive above the x-axis;

is the track yaw angle, as shown in fig. 2, i.e. the included angle between the velocity vector and the Oxh plane, and the same side of the velocity vector and the y-axis on the Oxh plane is positive; m is the rocket mass; e and sigma are used for expressing the direction of the thrust, wherein the e expresses the included angle between the thrust direction and the opposite direction of the speed; g is the acceleration of gravity at height h, g₀Corresponding to a gravitational acceleration of height 0, i.e.

T represents the magnitude of the thrust; d represents the aerodynamic resistance;

I_sp＝(270s)/(h₀/V₀) 16.784 is the specific impulse of a rocket engine. In the formula (1), variables other than the angle variable are subjected to dimensional normalization, and the position variables x, y, and h are normalized by the initial height h₀3700m, the velocity V is the initial velocity V₀230m/s, mass m is the initial mass m₀26000kg, time and beta_spBy using h₀/V₀For gravitational acceleration g, 16.087s

Thrust force T is used

To perform dimensional normalization respectively. Wherein the dimensionless resistance is expressed as

Where ρ is the dimensionless air density, varying with altitude,

is a dimensionless reference area of the rocket, C_D1.8 is the drag coefficient.

Introducing the constraints required for dynamic descent flight. Firstly, the magnitude of the thrust is restricted as follows,

T_min≤T≤T_max(3)

wherein T is_min＝0.5T_max0.417 and T_max0.834 is the minimum and maximum allowable thrust magnitude. Further, the thrust direction is constrained as follows

0≤∈≤∈_max(4)

Wherein e_maxMaximum allowed push at 15 °The force direction and the speed are at an angle in the opposite direction. Finally, to achieve accurate landing, the following terminal constraints need to be satisfied

x(t_f)＝0 (5)

y(t_f)＝0 (6)

h(t_f)＝0 (7)

V(t_f)≤V_f(8)

θ(t_f)＝-π/2 (9)

Wherein V_f＝8.696×10^-3. Constraints (5) - (7) ensure that the rocket lands on the designated landing site, and constraints (8) - (10) ensure that the landing speed is less than V_fAnd is perpendicular to the landing site ground.

The optimization objective is to minimize fuel economy, thus creating an optimal control problem for power-down landing as follows

s.t. formula (1), (3) - (10)

In the optimal control problem P0, the kinetic equation is highly nonlinear, with time of flight being an optimization variable. Obviously, problem P0 is a non-convex problem.

Step two: since the time of flight is free, the independent variables of the kinetic equation are converted into altitude, and the terminal constraints and the objective function are correspondingly converted.

The time of flight of the rocket during the power descent is unknown, but the initial and terminal altitudes are known, and during the power descent the altitude monotonically decreases over time, i.e., the rocket does not fly upwards. Thus converting the independent variable of formula (1) into a height

Wherein k is_D＝0.5ρS_refC_D＝5.917×10^-7ρ, superscript (') denotes the derivation of height h, and h is the normalized height from 1 to 0. The new arguments enable the terminal constraints (5) - (10) to be translated into

x(h_f)＝0 (13)

y(h_f)＝0 (14)

V(h_f)≤V_f(15)

θ(h_f)＝-π/2 (16)

Conversion of the objective function (11) into

The time-independent kinetic equation (1) is converted into the height-independent kinetic equation (12), and the corresponding terminal constraints and objective functions are also converted into the height-independent forms (13) - (18).

Step three: reducing the non-linearity of the kinetic equation. Part of nonlinearity in a kinetic equation is transferred into a constraint, then the nonlinearity is further reduced by reducing the dimension of the kinetic equation, two problems with smaller dimension are further established, namely a longitudinal problem and a lateral problem, and finally the lateral problem is simplified into a polynomial coefficient solving problem.

The primary nonlinearity of the original problem P0 is in the kinetic equation. In order to reduce the nonlinearity of the kinetic equation, three new control quantities are defined,

u₁:＝T cos∈/m (19)

u₂:＝T sin∈cosσ/m (20)

u₃:＝T sin∈sinσ/m (21)

and a new state variable,

ω:＝ln m (22)

the control variable nonlinearity in the third to fifth equations in the kinetic equation (12) is eliminated according to the defined variables, and the mass-cost equation is converted into

The thrust magnitude and direction constraints (3) - (4) are converted into

In addition, the optimization objective (18) is translated into

A reduction in the non-linearity of the kinetic equation results in an increase in the non-linearity of the constraint and objective function.

The non-linearity of the kinetic equation is further reduced by reducing the kinetic equation to dimensionality. The kinetic equation is divided into three parts, namely the longitudinal kinetic equation with respect to the state quantities x, V and theta

With respect to states y and

equation of lateral dynamics

And mass consumption equation (23).

The purpose of the kinetic equation decomposition is to transform the original problem P0 into a less-dimensional longitudinal optimal control problem and a lateral planning problem. The longitudinal dynamics equation (27) is included in the longitudinal problem. Since only u₁And u₂The constraints (24) - (25) on the magnitude and direction of thrust existing in the longitudinal dynamical equations, and therefore, in the longitudinal problem, need to be removed₃Converted into the following form

-u₁tan(∈_max-△_∈)≤u₂≤u₁tan(∈_max-△_∈) (30)

△ therein_TAnd △_∈Is the magnitude and direction of thrust reserved for lateral movement. When u is obtained by a lateral problem₃After u, u₁，u₂And u₃The constraints (24) - (25). △ that need to satisfy the original thrust magnitude and direction_TAnd △_∈Is calculated as follows

△ therein_∈,maxDetermine △ at 5 °_T＝0.0341T_max. The parameter κ is used to reflect the relative relationship between longitudinal and lateral maneuvers as

Where { x_f,soft,y_f,softAnd the terminal position of a soft landing point is obtained from the soft landing track obtained in the fourth step.

The longitudinal optimal control problem is written as

s.t. x′＝cotθ (35)

-u₁tan(∈_max-△_∈)≤u₂≤u₁tan(∈_max-△_∈) (39)

x(h_f)＝0,V(h_f)≤V_f,θ(h_f)＝-π/2 (40)

Problem(s)

Dependent on the unknown state quantities m and

the m acquisition method is realized in step four.

The planning problem of establishing a track yaw profile, namely the lateral problem, is as follows

Problem P_latThe objective is to plan a feasible track yaw profile to meet the terminal constraints of position y. Clearly, there are many possible track yaw profiles. The track yaw profile must satisfy certain properties. First, a feasible track yaw profile must meet

And

second, the corresponding state y should satisfy y (h)₀)＝y(h_f) 0. According to the Rollo theorem, by y (h)₀)＝y(h_f) When the value is 0, a point h exists_i∈[h_f,h₀]So that y' (h)_i) 0, according to formula (41),

by

It is known that there is a point h_e∈[h_f,h_i]So that

Thus, a feasible track yaw profile satisfies

And

dividing the flight path yaw angle into two sections for planning, namely h belongs to [ h ∈ [_e,h₀]And h e [ h ∈_f,h_e]. For the convenience of calculation, will

Viewed as a two-stage Bernstein polynomial, the first stage being a fourthThe second stage is a second order polynomial. Thus, it is possible to provide

Is written as

Wherein the Bernstein coefficient { ζ_1,i}_i＝0,…,4And { ζ_2,i}_i＝0,…,2Needs to be solved. The first stage polynomial has five coefficients, but only two conditions

And

therefore, in order to solve all the coefficients, three additional conditions need to be included,

and y (h)_f) 0. Wherein order u₃(h₀) Is equal to 0, to obtain

With the five conditions, Bernstein coefficients of the first-segment polynomial can be expressed as

The Bernstein coefficient of the second-segment polynomial may pass through the condition

And

to obtain

ζ_2,0＝0 (49)

The way these coefficients are solved is based on the properties of Bernstein polynomials. And (6) substituting all coefficients into an equation (43) to obtain a track yaw angle profile.

Part of the non-linearity of the kinetic equation is transferred to constraints (24) - (25) and then the non-linearity is further reduced by reducing the dimensionality of the kinetic equation, thereby creating two smaller dimensional problems, namely longitudinal problems

And lateral problems P_lat.(θ), last lateral problem P_lat.(θ) is converted into a Bernstein polynomial coefficient solving problem, and the coefficients are obtained by the equations (44) to (51).

Step four: further simplifying the longitudinal optimal control problem. The speed inclination angle is separated from the longitudinal kinetic equation, and the speed inclination angle is optimized independently, so that the nonlinearity of the longitudinal kinetic equation is reduced.

Longitudinal optimal control problem

The kinetic method ofThe range is still highly non-linear and needs further simplification. The non-linearity of the longitudinal dynamical equation is mainly due to the trigonometric function term of the velocity tilt angle θ. Thus solving the speed tilt angle from the longitudinal direction

And (4) separating, and independently optimizing the speed inclination angle. Thus, longitudinal problems

Is simplified into

-u₁tan(∈_max-△_∈)≤u₂≤u₁tan(∈_max-△_∈) (56)

V(h_f)≤V_f(57)

The kinetic equations for states x and θ are used to plan the velocity tilt angle profile and are therefore not a problem

In (1). Compared with the problem

Problem(s)

Less non-linearity, but still oneNon-convex problems. Step five will be right to problem

And (4) carrying out convex formation.

It is easy to obtain a feasible velocity tilt profile such that the state x satisfies the terminal constraints. However, a feasible but unreasonable velocity tilt angle profile may cause problems

It is not feasible. Problem(s)

Is mainly due to the thrust direction constraint (56) not being satisfied. | u₂Too large | easily results in the constraint (56) not being satisfied. Therefore, | u that needs to be generated when planning the velocity tilt angle₂| is as small as possible. Therefore, the following optimal control problem is obtained

s.t.x′＝cotθ (59)

x(h_f)＝0,θ(h_f)＝-π/2 (60)

Problem(s)

Dependent on V and

due to the problems

Has not been solved, so the velocity V is unknown, and furthermore due to the problem P_lat(theta) depends on theta, so track yaw angle

Is also unknown. But V and

present only in the objective function, hence V and

is sufficient to solve the problem

First, an approximate velocity profile is planned

Considering a soft landing problem, the thrust magnitude profile of the soft landing is T_min-T_maxThe bang-bang structure of (1) has a conversion time of t_sAnd the thrust direction is always the opposite direction of speed, i.e., ∈ ═ 0. Transition time t_sIs the variable to be sought. Using the thrust profile to integrate the kinetic equation (1) until V (t) is less than or equal to V_fAnd (5) stopping. Note that this time is t_fAnd the height at that moment is h_f(t_s). If h is_f(t_s) When t is equal to 0, then t is_sIs a solution to the soft landing problem. The soft landing problem is written as

P_soft:find t_s

s.t. formula (1), T is T_min-T_maxForm e is 0

h(t_f)＝0 (61)

Problem P_softIs a one-dimensional root finding problem. The solution to the problem can be found as t_s0.502. △ for the task is determined from the soft landing point _∈5 deg. is equal to. Using the velocity profile of the soft landing trajectory as an approximate velocity profile

Further, the profile described by the following formula is used as an approximate track yaw profile

Is called as in the above formula "

Initialization ".

When a problem arises

V in (1) and

quilt

And

after the substitution, problems

It is still a non-convex problem. In order to improve the efficiency of solving the problem, it must be transformed into a convex optimization problem. First, the term cot θ is expressed as a fifth-order Bernstein polynomial as follows

Therefore, equation (64) holds

Substituting the above formula into the problem

In the objective function of

The Bernstein polynomial B in the objective function (65) has six coefficients { ζ }_i}_i＝0,…,5Where four coefficients are represented analytically. According to the problems

The following equations (66) to (68) are clearly true

B(h₀)＝cotθ₀(66)

B(h_f)＝cotθ_f(67)

B′(h₀)＝θ₀′csc²θ₀(68)

Wherein theta is₀:＝θ(h₀) And theta_f:＝θ(h_f) Is a known quantity, θ₀′:＝θ′(h₀) Is the initial derivative of the velocity tilt angle. According to the properties of the polynomials (66) to (68) and Bernstein

ζ₀＝0 (69)

ζ₅＝cotθ₀(71)

Furthermore, the kinetic equation for state x is rewritten as

From the integral nature of the Bernstein polynomial, one can derive

ζ₃＝6x₀-[(ζ₀+ζ₄+ζ₅)+(ζ₁+ζ₂)](73)

Wherein ζ₀，ζ₄And ζ₅Obtained from formulae (69) to (71). Problem(s)

All ofConstraining for obtaining ζ₀，ζ₃，ζ₄And ζ₅. When its value is substituted into the objective function (65), the objective function depends only on ζ₁And ζ₂. Therefore, problems arise

Expressed as the following unconstrained optimization problem.

Fig. 4 is a graph and contour plot of the function F. The function F being a function of the variable ζ₁And ζ₂A convex function of (a). Thus, problems arise

The method is an unconstrained convex optimization problem and can be quickly and reliably solved by a quasi-Newton method. When ζ is₁And ζ₂After obtaining, the Bernstein coefficient { ζ_i}_i＝0,…,5A velocity gradient profile is obtained in place of equation (63).

In this step, the longitudinal optimal control problem is further simplified

Separating the velocity tilt angle from the longitudinal dynamics equation, reducing the non-linearity of the longitudinal optimal control problem to form a problem

For the separated velocity tilt angle, a problem is created

And then the problem is converted into an unconstrained convex optimization problem

Step five: and the simplified longitudinal optimal control problem is emphasized. Firstly, processing a nonlinear kinetic equation and a convex non-convex constraint condition, secondly, carrying out equivalent transformation on an objective function, and finally obtaining a convex longitudinal optimal control problem.

In the simplified longitudinal direction

In (2), the kinetic equation (53) is non-linear, the constraints (54) - (55) are non-convex, and the objective function (52) is also non-convex.

First, nonlinear kinetic equations (53) and equation constraints (54) are converted to linear. Definition V²As new state variables as follows

Then, the kinetic equation (53) and the equality constraint (54) are converted into

Formulae (76) to (77) are all linear.

Non-convex constraints are then processed (55). There are two approaches to the convex representation of the constraint (55), the first of which is described first. In the constraint (55), the inequality determined by the second inequality sign is a second order cone constraint and is a convex constraint. Due to | u₂Is far less than | u₁I, so the constraint determined by the first inequality sign is approximated by T_min/m≤u₁. Thus, the constraint (55) is approximated as

The regions determined by constraints (55) - (56) and the regions determined by constraints (78) and (56) are shown in fig. 5. The second processing method is a further simplification based on equation (78). The second order cone constraint is also reduced to a linear constraint. Then the non-convex constraint (55) is approximated as a linear constraint

In this example, a first method is chosen to handle the constraint (55).

And carrying out equivalent transformation on the non-convex objective function. The original objective function is replaced by a linear objective function

The objective functions (80) and (52) have an approximate optimization effect, which is the same even under certain conditions.

When the dynamic equation and the constraint are convexly processed and the objective function is converted, the original longitudinal problem

Is converted into a convex optimal control problem

-u₁tan(∈_max-△_∈)≤u₂≤u₁tan(∈_max-△_∈) (85)

In this step, the nonlinear kinetic equation (53) and the non-convex constraints (54) and (55) are embossed into a linear kinetic equation (76), a linear equation constraint (77), and either a convex constraint (78) or a linear constraint (79). Original longitudinal optimal control problem

Is embossed into

Step six: and sequentially solving the lateral problem, the speed inclination angle optimization problem and the longitudinal problem obtained in the third step to the fifth step. And then judging whether the track yaw angle profile needs to be improved and the quality profile needs to be updated according to the obtained solution, and further judging whether the longitudinal problem needs to be solved again. If not, outputting the obtained track; and if necessary, solving the longitudinal problem again, and updating the corresponding variable to obtain the track. The track is the track of the rocket power descending section planned in real time.

In steps one through five, the problem P0 is converted into a longitudinal problem

And lateral problems P_lat(theta). The longitudinal problem is then decomposed into a problem for planning the velocity-tilt angle profile

And problems with

To improve the solution efficiency, question

Convex optimization problem being convex to be unconstrained

Problem(s)

Is embossed into a problem

In order to obtain the final trajectory, the problem needs to be solved in a certain order. First, solve the soft landing problem P_softAnd execute "

Initialization "get

And

then based on

And

solving a velocity dip angle planning problem

θ and x are obtained. Then substituting theta into the flight path yaw angle planning problem P_lat(theta) obtaining

And y. Finally, the sum of the values of theta,

and

substitution problem

To obtain V, u₁And u₂. Through the solving process, the state quantities x, y, V, theta,

And a control quantity u₁、u₂. U is calculated by the following formula (87)₃

Due to u₃Is not included in the problem

And △, and_∈cannot completely guarantee u₁，u₂And u₃The original thrust direction constraint (25) is satisfied. Therefore, the e is calculated by the formula (88), and whether the constraint (25) is satisfied is detected.

If e exceeds the boundary e_maxThen need to be paired

The profile is improved, so that the original thrust direction constraint is met. First, it is necessary to find the height point h at which ∈ reaches the maximum value_sep. Then for h e [ h ∈ [_sep,h₀]Segment, make ∈ equal to ∈_maxAnd calculating u inversely by the following formula₃

Wherein

Symbol sum u₃The same is true. Improved track yaw profile

Can be integrated by

The kinetic equation is obtained, i.e.

For h e [ h [ [ h ]_f,h_sep]Section, according to the method for planning track yaw angle in step three

Is a two-stage Bernstein polynomial

Wherein the Bernstein coefficient can be obtained by the method in step three. While improving track yaw profile

After being obtained, the corresponding

Calculated by the following equation.

In addition, when solving the problem

Mass is an approximate mass profile with a soft landing

Instead. When u is obtained₁，u₂And u₃Then, a more accurate quality profile is obtained

If the track yaw needs to be improved, equation (93)In (1)

And u₃Are respectively as

And

based on the above two steps, two situations can occur. Firstly, the track yaw angle needs to be improved, and secondly, the updated quality profile m^dAnd approximate mass profile

With a large difference, i.e.

If neither of the two conditions occurs, the trajectory is found to be

Wherein m is m^dInstead. If one of the two situations occurs, the track yaw profile and the quality profile need to be updated. Then solve the problem again

And calculating new u₃. Then, whether the element belongs to the boundary element is judged_max. If the boundary is exceeded, the task is regarded as unreachable; if the boundary is not exceeded, the trajectory planning is successful, a new quality profile is calculated (93) again, and the trajectory is determined as

The solving step is summarized as pseudo code see fig. 6.

The method also comprises the seventh step: through the problem dimension reduction in the third step, the problem decomposition in the fourth step and the problem camouflaging in the fifth step, the trajectory of the rocket power descent section can be planned in the sixth step, namely, the rocket power descent section guidance is carried out through the rocket power descent section trajectory planned in real time in the sixth step, and the real-time trajectory planning efficiency and robustness of the rocket power descent section are improved.

By calculation, this example solves the longitudinal problem for the first time

Then e exceeds the boundary e_maxThere is a need for an improvement in track yaw. After the track yaw angle is improved, the longitudinal problem is solved again

And the trajectory planning of the rocket power descending process is realized.

The optimum fuel consumption was 4031.52kg and the flight time was 39.81 s. Fig. 7 illustrates the change in the speed bank angle and track yaw angle profiles. FIG. 8 shows thrust magnitude, T, in the velocity direction, velocity bank angle direction, and track yaw angle direction_i＝u_im, i is 1,2, 3. From the figure, it can be seen that | T₂I and I T₃Far less than | T₁L. Figure 9 shows the specific thrust magnitude and the change in the thrust direction angle against the speed. Fig. 10 shows the change in velocity, and the change in thrust acceleration and aerodynamic drag acceleration. As can be seen from the figure, aerodynamic drag plays a major role in the deceleration of the rocket in the early stages of landing. FIG. 11 shows the first planned and improved track yaw, first solving the problem

The rear thrust direction epsilon and the thrust direction epsilon after the track yaw angle is improved are changed, and the mass section difference of the front and the rear two times is changed. It can be seen from the figure that when the thrust direction constraint is not satisfied, it is satisfied by improving the track yaw angle. Fig. 12 shows a three-dimensional trajectory diagram of rocket flight and the direction of thrust. It can be seen from the figure that the rocket is landed in a vertical attitude.

The above detailed description is intended to illustrate the objects, aspects and advantages of the present invention, and it should be understood that the above detailed description is only exemplary of the present invention and is not intended to limit the scope of the present invention, and any modifications, equivalents, improvements and the like within the spirit and principle of the present invention should be included in the scope of the present invention.

Claims (8)

1. the real-time trajectory planning method of the rocket power descending section of guaranteeing convergence, is characterized in that: comprise the steps, Step 1: Carry out dynamic modeling and dimension normalization of the rocket power descent process, and establish a three-dimensional dimensionless dynamic equation; and introduce the constraints required for the power descent flight to establish the optimal control of the fuel-optimized power descent landing. question; Step 2: Since the flight time is free, convert the independent variables of the dynamic equation into heights, and convert the terminal constraints and objective functions accordingly; Step 3: Reduce the nonlinearity of the dynamic equation; first transfer part of the nonlinearity in the dynamic equation to the constraints, and then further reduce the nonlinearity by reducing the dimension of the dynamic equation, and then establish two Problems with smaller dimensions, namely longitudinal and lateral problems, and finally the lateral problem is reduced to a polynomial coefficient solution problem; Step 4: Further simplify the longitudinal optimal control problem; separate the velocity inclination angle from the longitudinal dynamics equation, and optimize the velocity inclination angle separately, thereby reducing the nonlinearity of the longitudinal dynamic equation; Step 5: Convex the simplified longitudinal optimal control problem; first deal with nonlinear dynamic equations and convex non-convex constraints, secondly perform the equivalent transformation of the objective function, and finally obtain a convex longitudinal optimal control problem ; Step 6: Solve the lateral problem, the velocity tilt angle optimization problem and the longitudinal problem obtained in the above steps 3 to 5 in sequence; Then it is judged whether it is necessary to re-solve the longitudinal problem; if not, output the obtained trajectory; if necessary, re-solve the longitudinal problem and update the corresponding variables to obtain the trajectory; the trajectory is the trajectory of the real-time planned rocket-powered descent segment. 2. the rocket power descending section real-time trajectory planning method of guaranteeing convergence as claimed in claim 1, is characterized in that: also comprises step 7, by the problem dimensionality reduction of step 3, the problem decomposition of step 4, the problem convex in step 5 In step 6, the trajectory of the rocket-powered descending segment can be planned, that is, the rocket-powered descending segment is guided through the rocket-powered descending segment trajectory planned in real time in step 6, so as to improve the efficiency and robustness of the real-time trajectory planning of the rocket-powered descending segment. 3. the rocket power descending section real-time trajectory planning method of guaranteeing convergence as claimed in claim 1 or 2, is characterized in that: step 1 concrete realization method is, The dynamics of the rocket-powered descent flight is modeled and the dimension is normalized. The dimensionless dynamic equation of the rocket-powered descent flight is expressed as

Among them, x, y and h are the position of the rocket, h is the height direction, pointing upward is positive, x is the direction from the projection point of the initial position on the horizontal plane to the landing point, y and h, x constitute the right-hand rule; V is the rocket's Velocity; θ is the velocity inclination angle, that is, the angle between the projection of the velocity vector on the Oxh plane and the x-axis, and the projection of the velocity is positive above the x-axis;

is the track yaw angle, that is, the angle between the velocity vector and the Oxh plane, the velocity vector and the y-axis are positive on the same side of the Oxh plane; m is the mass of the rocket; ∈ and σ are used to indicate the direction of thrust, where ∈ means The angle between the thrust direction and the opposite direction of the speed; g is the gravitational acceleration at the height h, g ₀ corresponds to the gravitational acceleration at a height of 0; T represents the magnitude of the thrust; D represents the aerodynamic resistance; I _sp is the specific impulse of the rocket engine ; In formula (1), except for the angle variable, all other variables are dimensionally normalized, the position variables x, y and h use the initial height h ₀ , the speed V uses the initial speed V ₀ , and the mass m uses the initial mass m ₀ , h ₀ /V ₀ for time and specific impulse I _sp , and gravitational acceleration g for

For thrust T

to perform dimension normalization separately; where the dimensionless resistance is expressed as

where ρ is the dimensionless air density, which varies with height, _Sref is the dimensionless reference area of the rocket, and C _D is the drag coefficient; Introduce the constraints required for powered descent flight; first, the thrust is constrained as follows T _min ≤T≤T _max (3) where T _min >0 and T _max are the minimum and maximum allowable thrust magnitudes; in addition, the thrust direction is constrained as follows 0≤∈≤∈ _max (4) where ∈ _max ∈[0,π/2) is the maximum allowable angle between the thrust direction and the reverse direction of the speed; finally, in order to achieve accurate landing, the following terminal constraints need to be satisfied x(t _f )=0 (5) y(t _f )=0 (6) h(t _f )=0 (7) V(t _f )≤V _f (8) θ(t _f )=-π/2 (9)

where V _f is a small safe landing speed; constraints (5)-(7) ensure that the rocket lands on the designated landing site, while constraints (8)-(10) ensure that the landing speed is less than V _f and is perpendicular to the ground at the landing site ; The optimization goal is to make the most fuel-efficient, so the optimal control problem of power descent and landing is established as follows

s.t. Formulas (1), (3)-(10) In the optimal control problem P0, the dynamic equations are highly nonlinear, and the flight time is an optimization variable; obviously, the problem P0 is a non-convex problem. 4. the rocket power descending section real-time trajectory planning method of guaranteeing convergence as claimed in claim 3 is characterized in that: the concrete realization method of step 2 is, The flight time of the rocket in the powered descent stage is unknown, but the initial and terminal altitudes are known, and in the powered descent stage, the altitude decreases monotonically with time, that is, the rocket will not fly upward; therefore, formula (1) The independent variable of is converted into height, we get

where k _D = 0.5ρS _ref C _D , the superscript (') represents the derivation of the height h, and h is the normalized height from 1 to 0; the new independent variable makes the terminal constraints (5)-( 10) Convert to x(h _f )=0 (13) y(h _f )=0 (14) V(h _f )≤V _f (15) θ(h _f )=-π/2 (16)

The objective function (11) is transformed into

The kinetic equation (1) with time as the independent variable is transformed into the dynamic equation with height as the independent variable (12), and the corresponding terminal constraints and objective functions are also transformed into the form with height as the independent variable (13)-(18 ). 5. the rocket power descending section real-time trajectory planning method of guaranteeing convergence as claimed in claim 4 is characterized in that: the concrete realization method of step 3 is, The main nonlinearity of the original problem P0 lies in the dynamic equation; in order to reduce the nonlinearity of the dynamic equation, three new control variables are defined as follows u ₁ :=T cos∈/m (19) u ₂ :=T sin∈cosσ/m (20) u ₃ :=T sin∈sinσ/m (21) and a new state variable ω:=ln m (22) According to the defined variables, the non-linearities of the control variables in equations 3 to 5 in the kinetic equation (12) are eliminated, and the mass consumption equation is transformed into

The thrust magnitude and direction constraints (3)-(4) are transformed into

Furthermore, the optimization objective (18) is transformed into

The decrease of nonlinearity of the dynamic equation brings about the increase of the nonlinearity of constraints and objective function; The nonlinearity of the dynamic equation is further reduced by reducing the number of dimensions of the dynamic equation; the dynamic equation is divided into three parts, namely the longitudinal dynamic equation with respect to the state quantities x, V and θ

about state y and

The lateral dynamics equation of

and mass consumption equation (23); The purpose of the dynamic equation decomposition is to transform the original problem P0 into a longitudinal optimal control problem and a lateral planning problem with smaller dimensions; the longitudinal dynamic equation (27) is included in the longitudinal problem; since only u ₁ and u ₂ exist Because of the longitudinal dynamics equation, the constraints (24)-(25) on the magnitude and direction of the thrust in the longitudinal problem need to remove u ₃ and transform into the following form

-u ₁ tan(∈ _max -△ _∈ )≤u ₂ ≤u ₁ tan(∈ _max -△ _∈ ) (30) where _ΔT and _Δ∈ are the thrust magnitude and direction reserved for lateral motion; when u ₃ is obtained through the lateral problem, u ₁ , u ₂ and u ₃ need to satisfy the original thrust magnitude and direction constraints (24 )-(25); △ _T and △ _∈ are calculated as follows

where △ _∈,max is a critical value to avoid △ _∈ being too large, and the parameter κ is used to reflect the relative relationship between longitudinal and lateral maneuvering, which is

where {x _f,soft ,y _f,soft } is the terminal position of a soft landing point, obtained from the soft landing trajectory obtained in step 4; The vertical optimal control problem is written as

s.t.x′=cotθ(35)

-u ₁ tan(∈ _max -△ _∈ )≤u ₂ ≤u ₁ tan(∈ _max -△ _∈ ) (39) x(h _f )=0, V(h _f )≤V _f , θ(h _f )=-π/2 (40) question

depends on the unknown state quantities m and

The acquisition method of m is realized in step 4; The planning problem of establishing the yaw angle profile of the track, that is, the lateral problem is as follows P _lat. (θ):

The purpose of the problem P _lat. (θ) is to plan a feasible track yaw angle profile that satisfies the terminal constraints at position y; obviously, there are many feasible track yaw angle profiles; the track yaw angle profile Certain properties must be satisfied; first, the feasible track yaw angle profile must satisfy

and

Secondly, the corresponding state y should satisfy y(h ₀ )=y(h _f )=0; according to Rollo’s theorem, from y(h ₀ )=y(h _f )=0, there exists a point h _i ∈[h _f , h ₀ ], so that y'( _hi )=0, according to equation (41),

Depend on

Knowing that there exists a point he _e ∈[h _f , _hi ] such that

Therefore, a feasible track yaw angle profile satisfies

and

The track yaw angle is divided into two sections for planning, namely h∈[h _e , h ₀ ] and h∈[h _f , h _e ]; for the convenience of calculation, the

Treated as a two-stage Bernstein polynomial, the first is a fourth-order polynomial and the second is a second-order polynomial; therefore

written as

where the Bernstein coefficients {ζ _1,i } _i=0,...,4 and {ζ _2,i } _i=0,...,2 need to be solved; the first polynomial has five coefficients, but only two conditions

and

Therefore, in order to solve all the coefficients, three additional conditions need to be included,

and y(h _f )=0; through the five conditions, the Bernstein coefficient of the first stage polynomial is expressed as

The Bernstein coefficient of the second segment polynomial can be obtained by the condition

and

get ζ _2,0 = 0 (49)

The solution method of these coefficients is based on the properties of Bernstein polynomials; Substitute all coefficients into equation (43) to obtain the track yaw angle profile; Part of the nonlinearity of the dynamics equations is transferred to constraints (24)-(25), and then the nonlinearity is further reduced by reducing the dimensionality of the dynamics equations, which in turn establishes two smaller dimensional problems, the longitudinal question

and the lateral problem _Plat. (θ), and finally the lateral problem _Plat. (θ) is transformed into a Bernstein polynomial coefficient solution problem, and the coefficients are obtained by equations (44)-(51). 6. the rocket power descending section real-time trajectory planning method of guaranteeing convergence as claimed in claim 5 is characterized in that: the concrete realization method of step 4 is, Longitudinal optimal control problem

The dynamic equations in

separated from the velocity pitch and optimized separately for the velocity inclination angle; therefore, the longitudinal problem

is simplified to

-u ₁ tan(∈ _max -△ _∈ )≤u ₂ ≤u ₁ tan(∈ _max -△ _∈ ) (56) V(h _f )≤V _f (57) The kinetic equations for states x and θ are used to plan the velocity inclination profile and are therefore not a problem

Medium; compared to the problem

question

Less nonlinear, but still a non-convex problem; step 5 will correct the problem

Convex; It is easy to obtain a feasible velocity inclination profile such that state x satisfies the terminal constraints; however, a feasible but unreasonable velocity inclination profile may cause problems

infeasible; problem

is infeasible mainly because the thrust direction constraint (56) is not satisfied; if |u ₂ | is too large, it is easy to cause the constraint (56) not to be satisfied; therefore, when planning the velocity inclination angle, it is necessary to make the |u ₂ | is as small as possible; therefore, the following optimal control problem is obtained

s.t.x′=cotθ(59) x(h _f )=0, θ(h _f )=-π/2 (60) question

depends on V and

due to problems

has not been solved yet, so the velocity V is unknown, and since the problem _Plat .(θ) depends on θ, the track yaw angle

is also unknown; but V and

exists only in the objective function, so with V and

An approximation of instead of suffices to solve the problem

First plan an approximate velocity profile

Consider a soft landing problem, the thrust magnitude profile of the soft landing is a bang-bang structure of T _min -T _max , its transition time is t _s , and the thrust direction is always the opposite direction of speed, that is, ∈=0; transition time t _s is the variable to be determined; use the thrust profile to integrate the dynamic equation (1) until V(t)≤V _f stops; denote this moment as t _f , and the height at this moment as h _f (t _s ); if h _f (t _s )=0, then the _ts is the solution of the soft landing problem; the soft landing problem is written as P _soft : find t _s st formula (1), T is in the form of T _min -T _max , ∈=0 h(t _f )=0 (61) The problem _Psoft is a one-dimensional root-finding problem; consider the velocity profile of the soft landing trajectory as an approximate velocity profile

In addition, the profile described by the following equation is used as the approximate track yaw angle profile

The above formula is called "

initialization"; when the problem

in V and

quilt

and

After replacing, the problem

is still a non-convex problem; in order to improve the solution efficiency of this problem, it must be transformed into a convex optimization problem; first, the term cotθ is expressed as a fifth-order Bernstein polynomial as follows

Therefore, Equation (64) holds

Substitute the above formula into the problem

In the objective function of , we get

The Bernstein polynomial B in the objective function (65) has six coefficients {ζ _i } _i=0,...,5 , four of which are expressed analytically; according to the problem

constraints, the following equations (66)-(68) obviously hold B(h ₀ )=cotθ ₀ (66) B(h _f )=cotθ _f (67) B'(h ₀ )=θ' ₀ csc ² θ ₀ (68) where θ ₀ :=θ(h ₀ ) and θ _f :=θ(h _f ) are known quantities, and θ′ ₀ :=θ′(h ₀ ) is the initial derivative of the velocity inclination angle; according to equation (66)- (68) and the properties of Bernstein polynomials, we get ζ ₀ = 0 (69)

ζ ₅ =cotθ ₀ (71) Furthermore, rewrite the kinetic equation for state x as

According to the integral properties of Bernstein polynomials, we get ζ ₃ =6x ₀ -[(ζ ₀ +ζ ₄ +ζ ₅ )+(ζ ₁ +ζ ₂ )] (73) where ζ ₀ , ζ ₄ and ζ ₅ are obtained from equations (69)-(71); the question

The full constraints of ζ 0 , ζ 3 , ζ 4 and ζ 5 are used to obtain ζ ₀ , ζ ₃ , ζ ₄ and ζ ₅ ; when its values are substituted into the objective function (65), the objective function depends only on ζ ₁ and ζ ₂ ; therefore, the problem

It can be expressed as the following unconstrained optimization problem

The function F is actually _a convex function with respect to the variables ζ1 and ζ2 _; hence the question

is an unconstrained convex optimization problem, which can be solved quickly and reliably by the quasi-Newton method; when ζ ₁ and ζ ₂ are obtained, the Bernstein coefficient {ζ _i } _i=0,...,5 is substituted into equation (63) to obtain the velocity slope angular profile; In this step, the longitudinal optimal control problem is further simplified

Separating the velocity inclination angle from the longitudinal dynamics equation reduces the nonlinearity of the longitudinal optimal control problem and forms the problem

For the separated velocity inclination angle, establish the problem

It is then transformed into an unconstrained convex optimization problem

7. the rocket power descending section real-time trajectory planning method of guaranteeing convergence as claimed in claim 6 is characterized in that: the concrete realization method of step 5 is, in the simplified longitudinal problem

, the dynamic equation (53) is nonlinear, the constraints (54)-(55) are non-convex, and the objective function (52) is also non-convex; First, transform the nonlinear dynamic equation (53) and equation constraints (54) into linear ones ^; define V2 as a new state variable as follows

Therefore, the kinetic equation (53) and the equation constraint (54) are transformed into

Equations (76)-(77) are all linear; Then, the non-convex constraint (55) is processed; there are two processing methods to convexize the constraint (55), the first one is introduced first; in the constraint (55), the inequality determined by the second inequality sign is a second order Cone constraint, which is a convex constraint; since |u ₂ | is much smaller than |u ₁ |, the constraint determined by the first inequality sign is approximated as T _min /m≤u ₁ ; thus, constraint (55) is approximated as

The second processing method is to further simplify on the basis of equation (78); the second-order cone constraint is also simplified to a linear constraint; therefore, the non-convex constraint (55) is approximated as the following linear constraint

Equivalently transform the non-convex objective function; replace the original objective function with the following linear objective function

The objective functions (80) and (52) have approximate optimization effects, and even under certain conditions, the optimization effects are the same; When the dynamic equations and constraints are convexed and the objective function is transformed, the original longitudinal problem

is transformed into a convex optimal control problem

-u ₁ tan(∈ _max -△ _∈ )≤u ₂ ≤u ₁ tan(∈ _max -△ _∈ ) (85)

In this step, the nonlinear dynamic equation (53) and the non-convex constraints (54) and (55) are convexized into linear dynamic equations (76), linear equality constraints (77), and convex constraints (78) or linear constraints (79); the original longitudinal optimal control problem

be convexed as

8. the rocket power descending section real-time trajectory planning method of guaranteeing convergence as claimed in claim 7 is characterized in that: the concrete realization method of step 6 is, In steps 1 to 5, problem P0 is transformed into a vertical problem

and the lateral problem P _lat. (θ); then, the longitudinal problem is decomposed into a problem for planning the velocity inclination profile

and question

In order to improve the solution efficiency, the problem

Convex into an unconstrained convex optimization problem

question

convexity as a problem

In order to get the final trajectory, the problem needs to be solved in a certain order; first solve the soft landing problem _Psoft and execute "

"initialize" to get

and

then based on

and

Solving the Velocity Inclination Planning Problem

Obtain θ and x; then, substitute θ into the track yaw angle planning problem _Plat. (θ) to obtain

and y; finally, set θ,

and

Substitution problem

Obtain V, u ₁ and u ₂ ; through the above solution process, obtain the state quantities x, y, V, θ,

and control quantities u ₁ , u ₂ ; u ₃ is calculated by the following formula (87)

Since u ₃ is not included in the question

In the thrust direction constraint of , and △ _∈ cannot fully guarantee that u ₁ , u ₂ and u ₃ satisfy the original thrust direction constraint (25); therefore, calculate ∈ by formula (88), and then check whether the constraint (25) is satisfied;

If ∈ exceeds the bound ∈ _max , we need to

The profile is improved so that the original thrust direction constraint is satisfied; first, it is necessary to find the height point h _sep where ∈ reaches the maximum value; then for the h∈[h _sep , h ₀ ] segment, make ∈ equal to ∈ _max , and inverse the following formula find u ₃

in

Same sign as u ₃ ; improved track yaw profile

by points

The kinetic equation is obtained, namely

For the h∈[h _f ,h _sep ] segment, according to the method of planning the track yaw angle in step 3, let

is the following two-segment Bernstein polynomial

The Bernstein coefficient is obtained by the method in step 3; when the improved track yaw angle profile

After getting, the corresponding

Calculated by formula (92);

Furthermore, when solving the problem

, the mass is replaced by the approximate mass profile m of the soft landing; when u ₁ , u ₂ and u ₃ are obtained, a more accurate mass profile is obtained

If the track yaw angle needs to be improved, in Eq. (93)

and u ₃ are respectively

and

Two situations will arise; one is that the track yaw angle needs to be improved, and the other is that the updated mass profile m ^d and the approximate mass profile

big difference, that is

where δ _m is the judgment m ^d and

critical value for the size of the gap; if neither of the above occurs, the trajectory is obtained as

where m is replaced by m ^d ; if one of the above two conditions occurs, the track yaw profile and mass profile need to be updated; then the problem is re-solved

and calculate the new u ₃ ; then judge whether ∈ exceeds the boundary ∈ _max ; if it exceeds the boundary, the task is considered unreachable; if it does not exceed the boundary, the trajectory planning is successful, and (93) is used to calculate the new mass profile again, and the trajectory is obtained as

Images