CN111196382B - Real-time trajectory planning method for rocket power descent segment capable of guaranteeing convergence - 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 weight of 0A force acceleration; 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, varying with altitude, 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 make the terminal constraints (5) - (10) changeInto

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)

Wherein Δ_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 original thrust magnitude and direction constraints (24) - (25) need to be satisfied. Delta_TAnd Δ_∈Is calculated as follows

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

Where { x_f,soft,y_f,softIs the terminal position of a soft landing site,and acquiring 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]. To countIs calculated conveniently, 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(θ) is converted into a Bernstein polynomial coefficient solving problem, the coefficients being obtained by equations (44) - (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.

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)

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)

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 first one will notThe constraint determined by equal sign is approximately 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

Is restricted in the thrust direction of (1), 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 quality profile is calculated (93) again, and the trajectory is determined 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.

To verify the feasibility of the method, rocket powers were chosenThe descent landing task is verified. 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 of_∈,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; theta is the velocity tilt angle, as shown in FIG. 2, i.e., the projection of the velocity vector onto the Oxh plane and the x-axisThe 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_max15 ° is the maximum allowable thrust direction and speed 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_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)

Wherein Δ_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 original thrust magnitude and direction constraints (24) - (25) need to be satisfied. Delta_TAnd Δ_∈Is calculated as follows

Wherein Δ_∈,maxAt 5 deg., find Delta_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

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. Wherein order u₃(h₀) Is equal to 0, to obtain

By the five conditions, theBernstein coefficients for a segment of a polynomial may 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.

Partial non-linearity of kinetic equationsIs transferred to constraints (24) - (25) and then further reduces the non-linearity by reducing the dimensionality of the kinetic equations, 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 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)

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)

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 time of dayHas a height of 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. And simultaneously, the delta of the task is obtained according to 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)

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

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

Is an unconstrained convex optimization problem, and can quickly and reliably solve the problem by a quasi-Newton methodAnd (5) solving. 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

Is restricted in the thrust direction of (1), and_∈cannot fully protectSyndrome 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

As follows two Bernstein polypeptidesPolynomial

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, in equation (93)

And u₃Are respectively as

And

based onTwo situations can occur in the two steps. 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. A real-time trajectory planning method for ensuring convergence of a rocket power descent segment is characterized by comprising the following steps: comprises the following steps of (a) carrying out,

the method comprises the following steps: carrying out dynamic modeling and dimensional normalization on the rocket power descending process, and establishing a three-dimensional dimensionless dynamic equation; introducing the constraint required by power descent flight, and establishing the optimal control problem of the power descent landing with optimal fuel;

step two: because the flight time is free, the independent variable of the dynamic equation is converted into the height, and the terminal constraint and the objective function are correspondingly converted;

step three: reducing the non-linearity of the kinetic equation; transferring part of nonlinearity in a kinetic equation into a constraint, then further reducing the nonlinearity by reducing the dimension of the kinetic equation, further establishing two problems with smaller dimension, namely a longitudinal problem and a lateral problem, and finally simplifying the lateral problem into a polynomial coefficient solving problem;

step four: the longitudinal optimal control problem is further simplified; separating the speed inclination angle from the longitudinal kinetic equation, and independently optimizing the speed inclination angle so as to reduce the nonlinearity of the longitudinal kinetic equation;

step five: the simplified longitudinal optimal control problem is emphasized; firstly, processing a nonlinear kinetic equation and a convex non-convex constraint condition, secondly, performing equivalent transformation on a target function, and finally obtaining a convex longitudinal optimal control problem;

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

2. A rocket power descent segment real-time trajectory planning method to ensure convergence as recited in claim 1, further comprising: and step seven, the problem in the step three is reduced in dimension, the problem in the step four is decomposed, the problem in the step five is highlighted, and the trajectory of the rocket power descent section can be planned in the step six, namely, the rocket power descent section guidance is carried out through the rocket power descent section trajectory planned in real time in the step six, so that the real-time trajectory planning efficiency and robustness of the rocket power descent section are improved.

3. A method for guaranteed convergence rocket power descent segment real-time trajectory planning as defined in claim 1 or 2, further comprising: the specific implementation method of the step one is that,

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 respectively perform dimension normalization; wherein the dimensionless resistance is expressed as

Where ρ is the dimensionless air density, varying with altitude, S_refIs a dimensionless reference area of the rocket, C_DIs the drag coefficient;

introducing the constraint required by the 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, pi/2) is the maximum allowable included angle between the thrust direction and the speed 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_fIs a small safe landing speed; constraints (5) - (7) ensure that the rocket lands on the fingerLanding site is fixed, and constraints (8) - (10) ensure that 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.

4. A rocket power descent segment real-time trajectory planning method to ensure convergence as recited in claim 3, further comprising: the concrete implementation method of the second step is that,

the time of flight of the rocket during the power descent phase is unknown, but the initial and terminal altitudes are known, and during the power descent phase, 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).

5. A rocket power descent segment real-time trajectory planning method to ensure convergence as recited in claim 4, further comprising: the third step is realized by the concrete method that,

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 decrease in the non-linearity of the kinetic equation brings about a nonlinear increase in the constraint and objective function;

further reducing the non-linearity of the kinetic equation 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 convert the original problem P0 into a longitudinal optimal control problem and a lateral planning problem with smaller dimensions; the longitudinal dynamics equation (27) is contained 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)

Wherein Δ_TAnd Δ_∈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) of the original thrust magnitude and direction need to be satisfied; delta_TAnd Δ_∈Is calculated as follows

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

Where { x_f,soft,y_f,softAcquiring the terminal position of a soft landing point 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 the fourth step;

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

P_lat.(θ):find

Problem P_lat.(θ) the objective is to plan a feasible track yaw profile to meet the terminal constraints for position y; obviously, there are manyA feasible track yaw profile; 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 Rolle's 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](ii) a For the convenience of calculation, will

Regarding the polynomial as two sections of Bernstein polynomials, the first section is a fourth-order polynomial, and the second section 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; 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; substituting all coefficients into a formula (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).

6. A method for guaranteed convergence rocket powerdown leg real-time trajectory planning as recited in claim 5 further comprising: the concrete implementation method of the step four is that,

longitudinal optimal control problem

The kinetic equations in (1) are still highly non-linear and need to be further simplified; the non-linearity of the longitudinal dynamical equation is mainly caused by a trigonometric function term of the velocity inclination angle theta; thus solving the speed tilt angle from the longitudinal direction

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

Performing the following steps; phase (C)Is more problematic than

Problem(s)

Less non-linearity, but still a non-convex problem; step five will be right to problem

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

Is not feasible; problem(s)

Is mainly due to the thrust direction constraint (56) not being satisfied; | u₂Too large tends to cause the constraint (56) to be unsatisfied; 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 the speed, namely the epsilon is 0; transition time t_sIs a variable to be solved; using the thrust profile to integrate the kinetic equation (1) until V (t) is less than or equal to V_fStopping; note that this time is t_fAnd the height at that moment is h_f(t_s) (ii) a 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; will be softVelocity profile of land trajectory as 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 solving efficiency of the problem, the problem needs to be converted 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,...,5Wherein the 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, theta'₀:＝θ′(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 ζ₅(ii) a When its value is substituted into the objective function (65), the objective function depends only on ζ₁And ζ₂(ii) a Therefore, problems arise

Unconstrained optimization problem expressed as

The function F is actually a function of the variable ζ₁And ζ₂A convex function of (d); 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,...,5Obtaining a velocity tilt angle profile in place of formula (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

7. A rocket power descent segment real-time trajectory planning method to ensure convergence as recited in claim 6, further comprising: the concrete implementation method of the step five is that,

in the simplified longitudinal direction

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

firstly, converting a nonlinear kinetic equation (53) and an equation constraint (54) into linearity; definition V²As new state variables as follows

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

Formulas (76) - (77) are all linear;

then, processing the non-convex constraints (55); there are two approaches to the convex processing of the constraint (55), first introduced; 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₁(ii) a Thus, the constraint (55) is approximated as

The second processing method is further simplified based on the formula (78); simplifying the second-order cone constraint into a linear constraint; thus, the non-convex constraint (55) is approximated as a linear constraint

Carrying out equivalent transformation on the non-convex target function; the original objective function is replaced by a linear objective function

The objective functions (80) and (52) have approximate optimization effects, even under certain conditions, the optimization effects are the same;

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 dynamical equation (53) and the non-convex constraints (54) and (55) are embossed into a linear dynamical equation (76), a linear equality constraint (77), and either a convex constraint (78) or a linear constraint (79); original longitudinal optimal control problem

Is embossed into

8. A rocket power descent segment real-time trajectory planning method to ensure convergence as recited in claim 7, further comprising: the concrete realization method of the sixth step is that,

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

And lateral problems P_lat.(θ); 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 a final track, the problem needs to be solved according to a certain sequence; first, solve the soft landing problem P_softAnd execute "

Initialization "get

And

then based on

And

solving a velocity dip angle planning problem

Obtaining theta and x; 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₂(ii) a Through the solving process, the state quantities x, y, V, theta,

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

Due to u₃Is not included in the problem

In the thrust direction of (2), and_∈cannot completely guarantee u₁，u₂And u₃The original thrust direction constraint (25) is satisfied; due to the fact thatCalculating the epsilon according to the formula (88), and further detecting whether the constraint (25) is met;

if e exceeds the boundary e_maxNeed to be aligned with

The section 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(ii) a 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; 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 the following two sections Bernstein polynomial

Wherein the Bernstein coefficient is obtained by the method in the third step; 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

Replacing; 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

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

a critical value of the difference; if neither of the two conditions occurs, the trajectory is found to be

Wherein m is m^dReplacing; if one of the two situations occurs, the track yaw angle profile and the quality profile need to be updated; then solve the problem again

And calculating new u₃(ii) a Then, whether the element belongs to the boundary element is judged_max(ii) a 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

Images