Abstract
This paper proposed an analytical iterative guidance method with the desired instantaneous impact point constraint for solid rockets in “burn-coast-burn” trajectory mode. Solid rocket motors expect to remove the thrust termination mechanism to increase the structural strength and launch reliability, which induce new difficulties and challenges to the guidance problems. In terms of the “Hohmann transfer” principle, a pointing algorithm is deduced in depth to establish the theoretical relations among the ignition time, the required velocity vector, and the orbital element constraints and provides the analytical expression of the ignition time. Then, an analytic solution of the required velocity vector is derived based on orthogonal and nonorthogonal velocity vectors, and a complete guidance logic is used to solve the target orbit elements satisfying the desired instantaneous impact point. Finally, the application of the developed theoretical algorithm in this paper is conducted using a two-stage solid rocket. The proposed guidance method is verified by Monte Carlo simulations, and the testing results demonstrate the adaptability, strong robustness, and excellent performance for different desired impact point missions.
1. Introduction
Solid rockets without the thrust termination mechanism can be launched rapidly with high manoeuvrability and agility, fulfilling the longed-for requirement of responsive launching missions, which brings new difficulties and challenges to its guidance technology [1]. For impact missions on the earth’s surface, the flight trajectory of solid rockets is typical “burn-coast” mode, and the coast arc accounts for a large proportion [2, 3]. Once the engine is shut down, the flight states at the time of engine shutdown determine the impact point on the spherical surface, and the coast arc follows the Kepler orbit law in the process. An instantaneous impact point (IIP) is defined as the touchdown point of a rocket with an assumption of an immediate end of the propelled flight, which simplifies the constraints on the terminal target mission [4]. Thus, the guidance method for solid rockets changes the current IIP to the desired location on the surface of the Earth. In addition, to improve the load-to-structure mass ratio and increase the strength of structure, the thrust termination mechanism is removed from solid rocket motors [5], which will cause the solid motors to shut off only by the fuel exhaustion instead of the main engine cut-off controlled by the guidance commands.
Over the past few decades, numerous efforts have been devoted to improve the accuracy and robustness of ascent guidance for launch vehicles. For a liquid launch vehicle shutting off controlled by the guidance command, a strongly adaptive iterative guidance method is designed according to the optimal analytical solution under multiple constraints, which is derived by the optimal control principle with the optimal fuel as the performance criteria. Countering the main characteristics of liquid launch vehicles, reference [6] derived an iterative guidance method with the thrust vector direction as control variable and position constraint ability in three-dimensional space based on the simplified rocket dynamics equations and the analytic solution of the optimal control problem. References [7–9] are concerned with the theoretical and algorithmic development for the fast generation of optimal launch ascent trajectories through the atmosphere and guidance, in which the critical roles of the primer vector in the determination of optimal endoatmospheric ascent trajectories are explored. Reference [10] proposes computational guidance and control (CG&C) based on the rapid development of onboard computers in recent years, such as pseudospectral optimal control and convex optimization [11–15]. Compared with the liquid rocket engines, precise engine cut-off capability is necessary for accurate instantaneous impact point targeting. However, when it is not capable of engine cut-off, the rocket is subject to unnecessary and undesirable acceleration, which continues until the fuel is exhausted and results in the impact point deviation.
For solid rocket motors, removing the thrust termination mechanism can not only increase the strength of the structure but also reduce the cost, which is the inevitable trend of development. These have posed new difficulties and challenges on the ascent guidance in vacuum, which indicate that solid rockets must depend on the advanced guidance to satisfy the terminal constraints and shut off by fuel exhaustion. To adapt to the guidance problem of exhaustion shutdown, the closed-loop guidance method that is insensitive to the cut-off time connects the flight mission with the thrust vector direction of the rocket engine by introducing the concept of the required velocity vector [16]. The method in [17] applies the theoretical assumption that the velocity increment generated by the engine is regarded as the instantaneous pulse vector, and the influence on the rocket position during the continuous burning period is ignored due to the short nominal working time of the solid rocket; so, the thrust vector direction is just always controlled to keep consistent with the direction of the required velocity vector. Reference [18] generates the general energy management (GEM) method according to the relation between the arc and the central angle, in the process of which the velocity increment in the required velocity vector direction is controlled by the closed-loop commands. References [19, 20] and [21] use the linear programming method to design an open-loop alternating attitude energy management algorithm (AEM) and spline energy management (SEM) guidance method, which further neutralizes the additional velocity component to strengthen the terminal constraints. Although these classical methods have been reliably and widely applied to the current launch vehicles, there are some limitations of the above methods in the solid rocket with burn-coast-burn flight mode such as the improvement of mission requirements, especially under the characteristics of fast response mission, because the closed-loop guidance method does not have multiterminal constraint capability, and the iterative guidance method is only applicable to the propulsion process.
In order to satisfy the instantaneous impact point constraints of solid rocket with coast arcs, this paper presents an analytical instantaneous impact point guidance method that combines the solution of ignition time and the iterative calculation of the target orbit. The study is organized as follows: In Section II, the dynamic model for solid rockets and the instantaneous impact point guidance problem is established. In Section III, the pointing algorithm for burn arc is proposed, and the analytical solution of ignition time is obtained. In Section IV, the iterative calculation of the target orbit considering the deleted shutdown and the solution logic of the proposed guidance method are presented. In Section V, the numerical simulation results are discussed. In Section VI, conclusions are drawn.
2. Instantaneous Impact Point Guidance Problem
2.1. Dynamic Model for Solid Rockets
A multistage solid rocket is investigated in this paper, in which the last stage flies in the exoatmospheric, and the typical ascent launch process is as shown in Figure 1. After the launch vehicle clears the dense atmosphere, the aerodynamic forces and atmosphere pressure disappear, and the algorithmic development unconstrained by the aerodynamic angle and bending moment limit focuses on the autonomous ascent guidance. Such a vacuum flight environment facilitates the study on the guidance methods of depleted shutdown for solid rockets.
The three-dimensional point-mass equations of motion ascending through the atmosphere can be expressed in an inertial frame as
The thrust magnitude and the axial and normal forces and are given by where and are state vectors; is the unit vectors defining the rocket body longitudinal direction, and is the rocket body normal direction; is the vacuum specific impulse, is the cross-sectional area of the engine nozzle,and is the aerodynamic reference area; is the dynamic pressure; and are the coefficient of axial force and the coefficient of normal force, respectively; is the gravitational acceleration and .
Note that the centrifugal force, Coriolis force, and the J2 gravitational term are ignored or simplified in the design of the guidance model but treated as interference to verify the accuracy and robustness of the guidance algorithm in the simulation test. In the “burn-coast-burn” trajectory mode, the thrust magnitude in Eq. (1) can be expressed as where is the ignition time of the rocket motor, and is the rated operating time of the rocket motor. In the coast period, the launch vehicle passively flights along the Kepler orbit. After the ignition of rocket motors in vacuum, the equations of motion of solid rockets can be expressed as where and are the current position and velocity vector, respectively; and are the position and velocity vector at ignition time, respectively; and are the terminal position and velocity vector at the time of engine depleted shutdown, respectively. According to the Ziorkovsky formula, the position increment and the velocity increment by consuming the whole fuel are determined by
Due to the loss of throttle modulation and the thrust termination mechanism, solid rockets cannot control the thrust magnitude and the engine cut-off; so, Eq. (4) is constrained by Eqs.(5) and (6).
2.2. Instantaneous Impact Point Constraints
An instantaneous impact point of the solid rocket is the intersection point between the Earth’s surface and the Kepler trajectory of the rocket in the absence of specific acceleration. According to the properties of the Kepler orbit, the IIP unit vector can be expressed using the terminal position and velocity in Eq. (4) as where is the unit of vector; and are the flight path angle and range angle, respectively. The flight path angle is expressed as
And the range angle is expressed as where the coefficients , , and are given as [22] where is the geocentric gravitational constant; is the angular momentum moment and . The IIP longitude and latitude in the inertial frame can be obtained from the components of the IIP unit vector as follows:
By considering the Earth’s rotation during the flight time, the IIP longitude in the Earth-centered Earth-fixed (ECEF) is given as where the superscript E represents the values in the ECEF frame; is the rotational rate of the Earth and is the current time; is the flight time to arrive at the impact point.
Lambert’s problem seeks for the orbit that connects the two points in a prespecified transfer time . Alternatively, the problem can be defined as the determination of the velocity vector of an object at position that can transfer the object to position after time . Thus, considering the flight characteristics of solid rockets, the guidance goal is to find the control profile that satisfies the dynamics (1), equality constraints (5) and (6), and the terminal condition (11).
3. Algorithm for Burn Arc
Following the initially propelled ascending phase and already beyond the atmosphere, a pointing algorithm based on the concept of impulsive orbit transfer is established, which concentrates on how the given direction the continuous thrust vector is equivalent to the velocity vector impulse in the orbit transfer [23]. The schematic diagram of the pointing algorithm is shown in Figure 2.
Where is the intersection point between the suborbital orbit and target orbit; is the geocentric distance at the intersection point ; , and are the geocentric distance and inertial velocity vector in the suborbital Kepler orbit at rated operating time, respectively; and are the inertial velocity vectors at the point of the suborbital orbit and target orbit, respectively. The dash line represents Hohmann transfer orbit, and the solid line is the trajectory generated by employing the given direction continuous thrust vector.
3.1. Principle of Equivalent Pulse
In this process, solid rockets with uncontrolled thrust intensity and direction are injected to the target orbit. So, the equations of motion in Eq. (4) in an inertial coordinate system can be expressed as where the variation on the velocity vector and position vector caused by gravity is formulated as
In accordance with the character of Kepler orbits, the moment of momentum of the vehicle is conserved. Thus, the change of the moment of momentum caused by the continuous thrust is expressed as follows.
It is assumed that the body longitudinal direction is constant and coincides with the direction to determine the equivalence relationship. After substituting Eqs. (5) and (6) into Eq. (13), the simplified formula is expressed as
Using Eq. (16), Eq. (17) can be rewritten as
Depending on the crossproduct of the vectors with the same direction that is zero vector, Eq. (17) is reduced to the following form with some extra algebraic manipulations:
In the vacuum phase, when the engine ignition time satisfies the time constraint, it is always possible to find an intersection line that is on both the suborbital orbit plane and the target orbit plane. Due to the moment of momentum conservation of Kepler orbits, the following equations can be achieved. On the suborbital orbit, the following equation exists
Slimily, on the target orbit, the following equation exists
Substituting Eq. (19) into Eq.(20), it can be obtained
The target orbit is constructed based on the requirement of the mission before launching; so, the intersection point is determined when the launch vehicle is on the suborbital orbit. Hence, at the point, the variation of momentum moment can be written in the concise form of where the geocentric distance of the intersection point is given by
Note that and represent the same target orbit with and . The conclusion of Eq. (22) proves that there exists a specific ignition time [23, 24], which makes the effect of “directional continuous propulsion” on the trajectory equivalent to that of “instantaneous pulse velocity vector” applied at the intersection point .
3.2. Calculation of Velocity Vector
The relationship of the states (,) and the orbital elements can be obtained in the Earth centered inertial coordinate frame. In order to clearly reveal the vector relations in the PA guidance process, translate the vectors and to the guidance coordinate system , and Figure 3 illustrates the relevant vector relationships and definitions.
where , , is the normal direction of the target flight plane and .
According to the equivalent pulse theory of the PA method, the velocity vectors can be calculated by the orbital elements. It is known that in the Kepler orbit, the following equations hold where and are the velocity components along the and directions, respectively. The subscript represents the velocity vector in the suborbital orbit and the velocity vector in the target orbit, respectively; is the semilatus rectum and . Then, the true anomaly is given by
The velocities in the suborbital orbit and in the target orbit can be calculated by Eqs. (24) and (25). Similarly, the velocity and can be obtained, according to the true anomaly at the point . The velocity vectors and in the guidance coordinate system are expressed as where represents azimuth deviation caused by the deviation between the current orbital inclination and the target orbital inclination at the point and is given by
By substituting Eq. (26) into the expression of the angular momentum change Eq. (22) with some algebraic manipulation, the magnitude and the direction of the thrust vector are obtained by
3.3. Coas Time Calculation
For an elliptical orbit, the eccentric anomalies at the initial and the final points are calculated to determine the transfer time as follows.
Then, the transfer time between the two points and can be calculated by where mean motion . From the intersection point , the time that the vehicle moves from the current position to the orbit intersection point can be achieved by employing the theory of Kepler orbit: where and are the eccentric anomalies on the coasting arc at the time that the equivalent impulse occurs and at the current time, respectively. Thus, the actual ignition time of solid rockets is given by
In addition, it is noted that and represent the same target orbit with and Although the insertion point of the target orbit has changed, the terminal orbital element constraints are still satisfied. In order to reveal the position vector relations in Eq. (23) during the PA guidance process, the coast time from the shutdown position to the target point can be calculated by the pulse point , and the expression is as follows:
Then, the coast time from the current position to the target point can be calculated by
Because of the Earth’s rotation during the flight time, the IIP longitude in the Earth-centered Earth-fixed will deviate from the expected longitude , resulting in the deviation of the rocket terminal impact point. The flight time in the suborbital and target orbit can be calculated by orbital elements, and the vectors of the PA guidance method are also solved by orbital elements. The suborbital elements can be calculated from the current state vector of the rocket; however, solving the orbit elements satisfying the terminal instantaneous impact point is one of the typical Lambert problems after the ignition time of the suborbital (or the geocentric distance ) is determined. Therefore, the target orbit elements are the terminal constraints of the PA guidance method, which will be described in detail in the next section.
4. Calculation of Target Orbit considering Deleted Shutdown
The pointing algorithm establishes an equivalent theoretical relationship between the continuous propulsion process and the instantaneous velocity increment process by solving the ignition time. To achieve the desired impact point, the velocity vector at the equivalent point needs to satisfy the exhaustion shutdown constraint in Eq. (5) and terminal constraint in Eq. (11). The relationship of each vector in the three-dimensional space is shown in Figure 4.
4.1. Target Orbit Calculation
According to Eq. (7), an instantaneous impact point at the equivalent point with the velocity vector can be obtained, and instantaneous impact point at the equivalent point with the velocity vector can be expressed as follows:
Then, according to the equivalent relation in the principle of PA, an instantaneous impact point and instantaneous impact point are both at the equivalent point . Therefore, in order to satisfy the terminal constraint Eq. (7), the required velocity vector needs to be solved by
where the range angle can be obtained by
And the required velocity vector can be defined as follows:
In addition, the velocity increment generated by the solid rocket motor is uncontrollable during the course of exhausting shutdown, and Eq. (5) should be satisfied; so, the modulus of the required velocity vector can be calculated:
The velocity vector and the position vector are calculated by the navigation system, and the target impact point and the velocity increment are preinput parameters. Thus, the required velocity vector can be obtained by Eqs. (39) and (40), in which the flight path angle should be meet Eq. (37) constraint. After obtaining the required velocity vector , the target orbital elements of PA can be expressed in the geocentric inertial coordinate system as follows:
where function takes the third component of a vector; , , and are the terminal semimajor axis, eccentricity, and inclination angle, respectively.
4.2. Solution of an Iterator Variable
There are several possible iteration schemes to search for the required velocity vector satisfies the desired instantaneous impact point in Eq. (11). The geometry describing the position vectors and decomposition of the initial velocity is shown in Figure 5.
The required velocity vector can be expressed in the target flight plane as well without loss of generality as follows: where is a chordal unit vector defined from the difference of the two position vectors as follows:
Relationships among the velocity components in perpendicular axes (, ) and skewed axes (, ) are expressed as the following:
It is known that for the family of orbits that connect the two position vectors and , initial velocity components and are expressed as follows:
In Eq. (47), the product of the two velocity components in the skewed axes is expressed as follows:
Since and are precalculated, and the parameters and are related variables, the two velocity components are constant. By putting Eq. (46) into Eq. (48), one can obtain a relationship between the radial and circumferential velocity components as follows:
In order to obtain the analytic relationship between iteration variables and constrained variables, Eq. (49) can be rearranged to be expressed using the required velocity vector and flightpath angle as follows: where the flight path angle is used for parameterizing and updating. The main advantage of setting as an iteration parameter over other variables lies in the fact that it is bounded [25, 26]. In addition to the natural bound (), additional bounds for the minimum and maximum values for flight path angle can be obtained from the target missions and the engine energy. To determine the interval of iteration variables, the minimum flight path angle is determined by the direction of the vector. The minimum flight path angle is expressed by
The maximum flightpath angle is obtained from the escape condition and is expressed as follows:
In general, the desired instantaneous impact point of solid rockets should be within a certain range, especially the lower bound of range is one of the important mission indicators. Due to the influence of interference and uncertainty in the flight process, the rocket cannot reach the target even under the optimal flight command. Thus, to ensure the iteration variable is meaningful, the minimum value of the required velocity satisfying the terminal constraints is as follows:
According to Eqs. (51) and (52), the upper and lower bounds of iteration variable can be obtained from the geocentric distance of intersection point in Eq. (23) and the desired instantaneous impact point in Eq. (11). Then, there are many ways to find the root of the nonlinear equation (50) within the minimum and maximum value of the flightpath angle until the exhausting shutdown condition in Eq. (40) is satisfied. Finally, the required velocity vector is calculated by Eq. (27), so that the target orbital elements of the PA method are obtained according to Eqs. (41), (42), and (43).
4.3. Guidance Logic and Process
An instantaneous impact point guidance method is presented in the previous sections. The pointing algorithm establishes the relationship between the ignition time and the required velocity vector at the point by the principle of equivalent pulse. Then, an analytic iterative method is used to solve the target orbit elements satisfying the desired instantaneous impact point. The detailed procedures of the proposed guidance method are as follows.
Step I. Determining the equivalent pulse position , the ignition time is calculated by Eq. (33), and the suborbital elements of the pointing algorithm are obtained by the current state variables provided by the navigation system.
Step II. Calculate the range angle . The desired instantaneous impact point is given by Eq. (11), and the range angle can be obtained by Eq. (9), in which the desired longitude and latitude are predetermined by the mission.
Step III. Calculate the bounds for the flight path angle and . According to the results calculated in steps I and II, the flight path angle and is obtained by Eqs. (51) and (52), respectively.
Step IV. Judge whether the required velocity condition is satisfied or not. The minimum value of the required velocity satisfying the terminal constraints is given by Eq. (53). If , go to step V; otherwise, the step stops, and the output flight path angle is .
Step V. Update the iteration variable . Check the convergence of the iteration and exit if the exhaustion shutdown condition in Eq. (5) is satisfied.
Step VI. Calculate the required velocity vector . According to the flight path angle calculated in step IV and the velocity modulus in Eq. (40), the required velocity vector is obtained by Eq. (39).
Step VII. Calculate the target orbit elements. The terminal semimajor axis , eccentricity , and inclination angle are calculated from Eqs. (41), (42), and (43) according to the parameters and , respectively.
Step VIII. Calculate the rocket body longitudinal direction . The direction of the thrust vector can be calculated by Eq. (29) and .
Step IX. Go to step V.
5. Simulation Results
To evaluate the performance of the proposed guidance method, several cases are applied in this section. Monte Carlo simulations are widely used in science, engineering, and finance to assess risks and enable decision making. Reference [18] provides a very readable exposition of the subject matter pertinent to aerospace applications. The uncertainties in Monte Carlo simulations may be classified into two categories: epistemic and aleatory uncertainties. Epistemic parameters are those that are deterministic and have fixed values, but simply not precisely known (such as the trajectory state and aerodynamic coefficients). Aleatory variables are truly random in nature, such as the atmospheric density.
5.1. Performance of the IIP Method
The solid rocket consists of two booster stages with a total mass of 15,000 kg and a rated maximum range of 4600 km. The offline guidance commands are implemented in the first boost stage. Then, the proposed guidance method is applied in the 2nd stage, and the ignition time is calculated by Eq. (33). After the solid motor runs out of fuel, the rocket flights along the Kepler orbit to the inertial impact point. It is noted that the rocket will reenter the atmosphere in the process of the Kepler orbit. This paper focuses on the guidance of the inertial impact point in the ascent phase, thus ignoring the effect of aerodynamic force in the reentry process.
At present, the existing guidance methods are mainly GEM [18] and SEM [21] aiming at the guidance problem of solid rockets with instantaneous impact point constrained. In order to evaluate the performance of the proposed IIP method, the three methods are compared, and the state profiles of the trajectories in comparison are plotted in Figure 6.
(a)
(b)
The simulation results show that the IIP, GEM, and SEM methods for solid rockets can satisfy the terminal longitude and latitude constraints. However, the pitch angle of the proposed IIP method is smaller and smoother than that of GEM and SEM, which is beneficial to the stability of the control system of solid rocket. Then, the performance of the solid rocket engine is significantly affected by ambient temperature. The nominal operating time deviation reaches 10% between the highest and lowest temperature conditions, which brings difficulties and challenges to the adaptability and robustness of guidance methods. The deviation and uncertainty in the simulation are shown in Table 1.
In each Monte Carlo simulation run, random dispersion, and uncertainties, vehicle mass properties, aerodynamics, and atmosphere are unknown to the guidance algorithm. In the simulations, the guidance algorithm is called at a 10 Hz rate to provide the guidance commands. Based on the above discussion, the desired range was set at 4500 km, and 250 Monte Carlo simulations are conducted. The profiles of the states are shown in Figure 7.
(a)
(b)
(c)
(d)
Figure 7 demonstrates the state curves and command curves of the rocket flight trajectories under random deviation and uncertainty. The range between the desired impact point and the launch point of the solid rocket is 4500 km. Figure 7(a) shows that the terminal ranges (zero altitude condition) are near the desired target point, and the maximum range deviation is 45 km. Figures 7(b) and 7(c) show the change process of flight velocity and flight path angles, respectively, and the terminal values are near the desired values. Figure 7(d) illustrates the performance of guidance commands with deviation and uncertainty. In the second stage of the rocket, the guidance commands satisfying the terminal range constraint can be calculated according to the different trajectories onboard. The solid rocket has a long-time coast arc after depletion shutdown, and the state (altitude, velocity, and flight path angle) deviations will increase with the flight time of coast arc. In general, there is a reentry phase at the end of the rocket flight trajectory to correct the deviation.
5.2. Adaptability to Different Missions
To verify the accuracy and robustness of the proposed guidance approach under different desired impact point conditions, extensive Monte Carlo simulations [27] are presented in this section. Three different target ranges of 1500 km, 3000 km, and 4500 km are set in the simulations. In each mission, 250 Monte Carlo simulations are performed, and the total number of simulations is 750.
Figure 8 demonstrates the state curves of the rocket flight trajectories under three missions with random deviation and uncertainty. The ranges between the desired impact point and the launch point of the solid rocket are 1500 km, 3000 km, and 4500 km, respectively. Figures 8(a) and 8(b) show that the terminal ranges and the flight path angles are near the desired values. The performance of guidance commands with deviation and uncertainty under different missions is shown in Figure 9. Because the coast arc is much larger than the burn arc, the state deviation at the shutdown time will enlarge the range deviation along the Kepler orbit.
(a)
(b)
Monte Carlo simulation tests demonstrate the proposed guidance method has high accuracy and robustness, in which the mean deviation of the range reaches the magnitude of 15.84 (km), and the standard deviation is less than 42.31 (km); the mean deviation and standard deviation of longitude and latitude are magnitudes of 10-1 (°), respectively. The statistic results on the same 750 dispersed trajectories are listed in Table 2, and a scatter plot of the terminal deviations of longitude and latitude is depicted in Figure 10.
6. Conclusion
This paper presents an analytical guidance method with coast arcs for solving an instantaneous impact point constrained problem during exhausting shutdown. A pointing algorithm with the terminal orbital elements constraints is developed to establish the theoretical relationship between the ignition time and the required velocity vector in “burn-coast-burn” trajectory, which is applicable to solid motor shutting off by fuel exhaustion. Then, according to Lambert’s principle, the analytic solution of the required velocity vector is derived, and the iteration of orbit elements satisfying the desired impact point is constructed. The theoretical relationship between the burn arc and the coast arc of the rocket is established by the constraints of the orbital elements satisfying the desired impact point, which improves the adaptability of the guidance algorithm to the impact point missions.
Monte Carlo simulations using a two-stage solid rocket and different desired impact point (range) missions are performed. The guidance logic and process from the ignition time of the rocket to the desired impact point constraint are given in this paper, and the guidance commands are calculated onboard by the current state provided by the navigation. The simulation testing results demonstrate that the proposed instantaneous impact point guidance has high accuracy and strong robustness against to dispersion and uncertainty and strong adaptability to different desired impact points.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this article.