# Instantaneous Impact Point Guidance with Coast Arcs for Solid Rockets

Authors: Shengwei Jia; Fugui Li; Mingbo Tong
Journal: International Journal of Aerospace Engineering
 (2021)

---

## 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.


---

## Body

## 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 Figure1. 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.Figure 1
Typical launch process of the instantaneous impact point guidance for solid rockets.The three-dimensional point-mass equations of motion ascending through the atmosphere can be expressed in an inertial frame as(1)r.=v,v.=gr+T−A⋅xb−N⋅yb/mt.The thrust magnitudeT and the axial and normal forces A and N are given by
(2)T=Ispg0m.−Ph⋅smA=qSCAN=qSCN,where r and v are state vectors; xb is the unit vectors defining the rocket body longitudinal direction, and yb is the rocket body normal direction; Isp is the vacuum specific impulse, sm is the cross-sectional area of the engine nozzle,and S is the aerodynamic reference area; q is the dynamic pressure; CA and CN are the coefficient of axial force and the coefficient of normal force, respectively; g is the gravitational acceleration and g0=9.8m/s2.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 magnitudeT in Eq. (1) can be expressed as
(3)T=Ispg0⋅m.tig≤t<tig+Ts0t<tigt≥tig+Ts,where tig is the ignition time of the rocket motor, and Ts 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
(4)vf=vig+∫0Tsgr+T⋅xb/mtdt,rf=rig+vigTs+∫0Ts∫0tgr+T⋅xb/mτdτdt,where r0 and v0 are the current position and velocity vector, respectively; rig and vig are the position and velocity vector at ignition time, respectively; rf and vf are the terminal position and velocity vector at the time of engine depleted shutdown, respectively. According to the Ziorkovsky formula, the position increment RM and the velocity increment WM by consuming the whole fuel are determined by
(5)WM=Ispg0⋅lnm0/m0−ms,(6)RM=Ispg0ms−WMm0−ms/m..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 vector1rp can be expressed using the terminal position rf and velocity vf in Eq. (4) as
(7)1rp=cosβe+ϑlcosϑl1rf+sinβecosϑl1vf,where 1⋅ is the unit of vector; ϑl and βe are the flight path angle and range angle, respectively. The flight path angle ϑl is expressed as
(8)ϑl=arcsinrf⋅vfrf×vf.And the range angleβe is expressed as
(9)sinβe=A12A32−A12+A22A32−A22+A1A3/A12+A22,where the coefficients A1, A2, and A3 are given as [22]
(10)A1=−rfvf2μsinϑlcosϑl=−hμrfrf⋅vf,A2=hμrf−1;A3=hμrp−1,where μ is the geocentric gravitational constant; h is the angular momentum moment and h=rf×vf. The IIP longitude λp and latitude ϕp in the inertial frame can be obtained from the components of the IIP unit vector as follows:
(11)1rp∗=cosλpcosϕpλpcosϕpsinϕp.By considering the Earth’s rotation during the flight time, the IIP longitude in the Earth-centered Earth-fixed (ECEF) is given as(12)λp=λpE+ωet+tF,where the superscript E represents the values in the ECEF frame; ωe is the rotational rate of the Earth and t is the current time; tF 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 timetF. Alternatively, the problem can be defined as the determination of the velocity vector vf of an object at position rf that can transfer the object to position rp after time tF. Thus, considering the flight characteristics of solid rockets, the guidance goal is to find the control profile xbt that satisfies the dynamics (1), equality constraints (5) and (6), and the terminal condition (11).


## 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 Figure1. 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.Figure 1
Typical launch process of the instantaneous impact point guidance for solid rockets.The three-dimensional point-mass equations of motion ascending through the atmosphere can be expressed in an inertial frame as(1)r.=v,v.=gr+T−A⋅xb−N⋅yb/mt.The thrust magnitudeT and the axial and normal forces A and N are given by
(2)T=Ispg0m.−Ph⋅smA=qSCAN=qSCN,where r and v are state vectors; xb is the unit vectors defining the rocket body longitudinal direction, and yb is the rocket body normal direction; Isp is the vacuum specific impulse, sm is the cross-sectional area of the engine nozzle,and S is the aerodynamic reference area; q is the dynamic pressure; CA and CN are the coefficient of axial force and the coefficient of normal force, respectively; g is the gravitational acceleration and g0=9.8m/s2.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 magnitudeT in Eq. (1) can be expressed as
(3)T=Ispg0⋅m.tig≤t<tig+Ts0t<tigt≥tig+Ts,where tig is the ignition time of the rocket motor, and Ts 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
(4)vf=vig+∫0Tsgr+T⋅xb/mtdt,rf=rig+vigTs+∫0Ts∫0tgr+T⋅xb/mτdτdt,where r0 and v0 are the current position and velocity vector, respectively; rig and vig are the position and velocity vector at ignition time, respectively; rf and vf are the terminal position and velocity vector at the time of engine depleted shutdown, respectively. According to the Ziorkovsky formula, the position increment RM and the velocity increment WM by consuming the whole fuel are determined by
(5)WM=Ispg0⋅lnm0/m0−ms,(6)RM=Ispg0ms−WMm0−ms/m..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 vector1rp can be expressed using the terminal position rf and velocity vf in Eq. (4) as
(7)1rp=cosβe+ϑlcosϑl1rf+sinβecosϑl1vf,where 1⋅ is the unit of vector; ϑl and βe are the flight path angle and range angle, respectively. The flight path angle ϑl is expressed as
(8)ϑl=arcsinrf⋅vfrf×vf.And the range angleβe is expressed as
(9)sinβe=A12A32−A12+A22A32−A22+A1A3/A12+A22,where the coefficients A1, A2, and A3 are given as [22]
(10)A1=−rfvf2μsinϑlcosϑl=−hμrfrf⋅vf,A2=hμrf−1;A3=hμrp−1,where μ is the geocentric gravitational constant; h is the angular momentum moment and h=rf×vf. The IIP longitude λp and latitude ϕp in the inertial frame can be obtained from the components of the IIP unit vector as follows:
(11)1rp∗=cosλpcosϕpλpcosϕpsinϕp.By considering the Earth’s rotation during the flight time, the IIP longitude in the Earth-centered Earth-fixed (ECEF) is given as(12)λp=λpE+ωet+tF,where the superscript E represents the values in the ECEF frame; ωe is the rotational rate of the Earth and t is the current time; tF 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 timetF. Alternatively, the problem can be defined as the determination of the velocity vector vf of an object at position rf that can transfer the object to position rp after time tF. Thus, considering the flight characteristics of solid rockets, the guidance goal is to find the control profile xbt 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.Figure 2
The principle of the pointing algorithm and the relationship between vectors.WherePimp is the intersection point between the suborbital orbit and target orbit; rimp is the geocentric distance at the intersection point Pimp; rsub.f, and vsub.f are the geocentric distance and inertial velocity vector in the suborbital Kepler orbit at rated operating time, respectively; vsub.imp and vorb.imp are the inertial velocity vectors at the point Pimp 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
(13)vf=vig+ΔvG+WMΓ,rf=rig+vigTs+ΔrG+RMΓ,where the variation on the velocity vector and position vector caused by gravity is formulated as
(14)ΔvG=∫0Tsgrdt;ΔrG=∫0TsΔvGdt.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.(15)ΔH=mfrf×vf−m0rig×vig.It is assumed that the body longitudinal directionxb 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
(16)rf=rsub.f+RM⋅Γ,vf=vsub.f+WM⋅Γ.Using Eq. (16), Eq. (17) can be rewritten as
(17)ΔH=mfrsub.f×vsub.f+rsub.f×WMΓ+RMΓ×vsub.f+RMΓ×WMΓ−m0rig×vig.Depending on the crossproduct of the vectors with the same directionRMΓ×WMΓ that is zero vector, Eq. (17) is reduced to the following form with some extra algebraic manipulations:
(18)ΔH=mf⋅rsub.f×vsub.f+mfrsub.f−RMWMvsub.f×WM⋅Γ−m0⋅rig×vig.In the vacuum phase, when the engine ignition time satisfies the time constraint, it is always possible to find an intersection lineRimp 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
(19)m0r0×v0=m0rimp×vsub.imp=m0rsub.f×vsub.f.Slimily, on the target orbit, the following equation exists(20)mfrf×vf=mfrimp×vorb.imp.Substituting Eq. (19) into Eq.(20), it can be obtained
(21)ΔH=mfrimp×vsub.imp+mfrsub.f−RMWMvsub.f×WMΓ−m0rimp×vsub.imp.The target orbit is constructed based on the requirement of the mission before launching; so, the intersection pointPimp is determined when the launch vehicle is on the suborbital orbit. Hence, at the Pimp point, the variation of momentum moment can be written in the concise form of
(22)ΔH=mfrimp×vsub.imp+WMΓ−m0rimp×vsub.imp,where the geocentric distance of the intersection point rimp is given by
(23)rimp=rsub.f−RM/WM⋅vsub.f.Note thatvorb.imp and rimp represent the same target orbit with vf and rf. 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 Pimp.


### 3.2. Calculation of Velocity Vector

The relationship of the states (v,r) 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 vsub.imp and vorb.imp to the guidance coordinate system Op−1x1y1z, and Figure 3 illustrates the relevant vector relationships and definitions.Figure 3
The velocity vector relations of orbits at the pointPimp.where1y=1rimp, 1x=1z×1y, 1z is the normal direction of the target flight plane and 1z=rimp×rp∗/rimprp∗⋅sinβe..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(24)v⋅.x=μ/p⋅⋅1+e⋅cosf⋅,v⋅.y=μ/p⋅⋅e⋅sinf⋅,where v⋅.x and v⋅.y are the velocity components along the 1x and 1y directions, respectively. The subscript ⋅ represents the velocity vector vsub.imp in the suborbital orbit and the velocity vector vorb.imp in the target orbit, respectively; p is the semilatus rectum and p=a⋅1−e2. Then, the true anomaly f is given by
(25)f⋅=0,whene⋅=0,f⋅=arccosp⋅−r/e⋅⋅r,whene⋅≠0.The velocitiesvsub.x in the suborbital orbit and vorb.x in the target orbit can be calculated by Eqs. (24) and (25). Similarly, the velocity vsub.y and vorb.y can be obtained, according to the true anomaly f at the point Pimp. The velocity vectors vsub.imp and vorb.imp in the guidance coordinate system Op−1x1y1z are expressed as
(26)vsub.imp=vsub.imp.xcosδivsub.imp.yvsub.imp.xsinδiT,(27)vorb.imp=vorb.imp.xvorb.imp.y0T,where δi represents azimuth deviation caused by the deviation between the current orbital inclination and the target orbital inclination at the point Pimp and is given by
(28)δi=arcsincosiorbcosϕimp−arcsincosisubcosϕimp.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
(29)WPA=vorb.imp−vsub.imp2;Γ=vorb.imp−vsub.imp/WPA.


### 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.(30)E=arccosa⋅−R/a⋅⋅e⋅.Then, the transfer time between the two pointsa and b can be calculated by
(31)ta−b=ta+1n⋅Eb−Ea−esinEb−sinEa,where mean motion n⋅=μ/a⋅31/2. From the intersection point Pimp, the time t0−imp that the vehicle moves from the current position to the orbit intersection point can be achieved by employing the theory of Kepler orbit:
(32)t0−imp=t0+1nsubEimp−E0−esubsinEimp−sinE0,where E0 and Eimp 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
(33)tig=t0−imp+RM/WM−Ts.In addition, it is noted thatvorb.imp and rimp represent the same target orbit with vf and rf. 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 trf−p∗ from the shutdown position to the target point can be calculated by the pulse point Pimp, and the expression is as follows:
(34)trf−p∗=Ep∗−Eimp−eorb⋅sinEp∗−sinEimp/norb−RM/WM.Then, the coast timetF from the current position r0 to the target point rp∗ can be calculated by
(35)tF=tig+Ts+trf−p∗.Because of the Earth’s rotation during the flight time, the IIP longitudeλpE in the Earth-centered Earth-fixed will deviate from the expected longitude λp, 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 rimp) 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.


## 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
(13)vf=vig+ΔvG+WMΓ,rf=rig+vigTs+ΔrG+RMΓ,where the variation on the velocity vector and position vector caused by gravity is formulated as
(14)ΔvG=∫0Tsgrdt;ΔrG=∫0TsΔvGdt.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.(15)ΔH=mfrf×vf−m0rig×vig.It is assumed that the body longitudinal directionxb 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
(16)rf=rsub.f+RM⋅Γ,vf=vsub.f+WM⋅Γ.Using Eq. (16), Eq. (17) can be rewritten as
(17)ΔH=mfrsub.f×vsub.f+rsub.f×WMΓ+RMΓ×vsub.f+RMΓ×WMΓ−m0rig×vig.Depending on the crossproduct of the vectors with the same directionRMΓ×WMΓ that is zero vector, Eq. (17) is reduced to the following form with some extra algebraic manipulations:
(18)ΔH=mf⋅rsub.f×vsub.f+mfrsub.f−RMWMvsub.f×WM⋅Γ−m0⋅rig×vig.In the vacuum phase, when the engine ignition time satisfies the time constraint, it is always possible to find an intersection lineRimp 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
(19)m0r0×v0=m0rimp×vsub.imp=m0rsub.f×vsub.f.Slimily, on the target orbit, the following equation exists(20)mfrf×vf=mfrimp×vorb.imp.Substituting Eq. (19) into Eq.(20), it can be obtained
(21)ΔH=mfrimp×vsub.imp+mfrsub.f−RMWMvsub.f×WMΓ−m0rimp×vsub.imp.The target orbit is constructed based on the requirement of the mission before launching; so, the intersection pointPimp is determined when the launch vehicle is on the suborbital orbit. Hence, at the Pimp point, the variation of momentum moment can be written in the concise form of
(22)ΔH=mfrimp×vsub.imp+WMΓ−m0rimp×vsub.imp,where the geocentric distance of the intersection point rimp is given by
(23)rimp=rsub.f−RM/WM⋅vsub.f.Note thatvorb.imp and rimp represent the same target orbit with vf and rf. 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 Pimp.


## 3.2. Calculation of Velocity Vector

The relationship of the states (v,r) 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 vsub.imp and vorb.imp to the guidance coordinate system Op−1x1y1z, and Figure 3 illustrates the relevant vector relationships and definitions.Figure 3
The velocity vector relations of orbits at the pointPimp.where1y=1rimp, 1x=1z×1y, 1z is the normal direction of the target flight plane and 1z=rimp×rp∗/rimprp∗⋅sinβe..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(24)v⋅.x=μ/p⋅⋅1+e⋅cosf⋅,v⋅.y=μ/p⋅⋅e⋅sinf⋅,where v⋅.x and v⋅.y are the velocity components along the 1x and 1y directions, respectively. The subscript ⋅ represents the velocity vector vsub.imp in the suborbital orbit and the velocity vector vorb.imp in the target orbit, respectively; p is the semilatus rectum and p=a⋅1−e2. Then, the true anomaly f is given by
(25)f⋅=0,whene⋅=0,f⋅=arccosp⋅−r/e⋅⋅r,whene⋅≠0.The velocitiesvsub.x in the suborbital orbit and vorb.x in the target orbit can be calculated by Eqs. (24) and (25). Similarly, the velocity vsub.y and vorb.y can be obtained, according to the true anomaly f at the point Pimp. The velocity vectors vsub.imp and vorb.imp in the guidance coordinate system Op−1x1y1z are expressed as
(26)vsub.imp=vsub.imp.xcosδivsub.imp.yvsub.imp.xsinδiT,(27)vorb.imp=vorb.imp.xvorb.imp.y0T,where δi represents azimuth deviation caused by the deviation between the current orbital inclination and the target orbital inclination at the point Pimp and is given by
(28)δi=arcsincosiorbcosϕimp−arcsincosisubcosϕimp.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
(29)WPA=vorb.imp−vsub.imp2;Γ=vorb.imp−vsub.imp/WPA.


## 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.(30)E=arccosa⋅−R/a⋅⋅e⋅.Then, the transfer time between the two pointsa and b can be calculated by
(31)ta−b=ta+1n⋅Eb−Ea−esinEb−sinEa,where mean motion n⋅=μ/a⋅31/2. From the intersection point Pimp, the time t0−imp that the vehicle moves from the current position to the orbit intersection point can be achieved by employing the theory of Kepler orbit:
(32)t0−imp=t0+1nsubEimp−E0−esubsinEimp−sinE0,where E0 and Eimp 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
(33)tig=t0−imp+RM/WM−Ts.In addition, it is noted thatvorb.imp and rimp represent the same target orbit with vf and rf. 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 trf−p∗ from the shutdown position to the target point can be calculated by the pulse point Pimp, and the expression is as follows:
(34)trf−p∗=Ep∗−Eimp−eorb⋅sinEp∗−sinEimp/norb−RM/WM.Then, the coast timetF from the current position r0 to the target point rp∗ can be calculated by
(35)tF=tig+Ts+trf−p∗.Because of the Earth’s rotation during the flight time, the IIP longitudeλpE in the Earth-centered Earth-fixed will deviate from the expected longitude λp, 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 rimp) 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 vectorvorb.imp at the equivalent point Pimp 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.Figure 4
Calculation target orbit composed of the equivalent pulse pointPimp and the desired instantaneous impact point.


### 4.1. Target Orbit Calculation

According to Eq. (7), an instantaneous impact point 1rp.sub at the equivalent point Pimp with the velocity vector vsub.imp can be obtained, and instantaneous impact point 1rp.orb at the equivalent point Pimp with the velocity vector vorb.imp can be expressed as follows:
(36)1rp.orb=cosβe+ϑlcosϑl1rimp+sinβecosϑl1vorb.imp.Then, according to the equivalent relation in the principle of PA, an instantaneous impact point1rp.sub and instantaneous impact point 1rp.orb are both at the equivalent point Pimp. Therefore, in order to satisfy the terminal constraint Eq. (7), the required velocity vector vorb.imp needs to be solved by
(37)1rp∗=cosβe+ϑlcosϑl1rimp+sinβecosϑl1vorb.imp,where theβe range angle can be obtained by
(38)βe=arcsinrimp×rp∗/rimprp∗.And the required velocity vectorvorb.imp can be defined as follows:
(39)vorb.imp=vx⋅1x+vy⋅1y,=vorb.imp⋅cosϑl⋅1x+sinϑl⋅1y.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 vorb.imp can be calculated:
(40)vorb.imp=vsub.imp2+WM2.The velocity vectorvsub.imp and the position vector rimp are calculated by the navigation system, and the target impact point 1rp∗ and the velocity increment WM are preinput parameters. Thus, the required velocity vector vorb.imp can be obtained by Eqs. (39) and (40), in which the flight path angle ϑl should be meet Eq. (37) constraint. After obtaining the required velocity vector vorb.imp, the target orbital elements of PA can be expressed in the geocentric inertial coordinate system as follows:
(41)aorb=−μrimp/2μ−rimp⋅vorb.imp2,(42)eorb=vorb.imp2−μ/rimprimp−rimp⋅vorb.impvorb.imp/μ,(43)iorb=arccos1z3,where function⋅3 takes the third component of a vector; aorb, eorb, and iorb 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 vectorvorb.imp 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.Figure 5
Calculation required velocity vector according to different velocity components.The required velocity vectorvorb.imp can be expressed in the target flight plane as well without loss of generality as follows:
(44)vorb.imp=vx⋅1x+vy⋅1y=vc⋅1c+vρ⋅1y,where 1c is a chordal unit vector defined from the difference of the two position vectors as follows:
(45)rc=rp∗−rimp=rc⋅1c.Relationships among the velocity components in perpendicular axes (vx, vy) and skewed axes (vc, vρ) are expressed as the following:
(46)vc=rc⋅vxrp∗⋅sinβe;vρ=vy−rp∗⋅cosβe−rimprp∗sinβe⋅vx.It is known that for the family of orbits that connect the two position vectorsrimp and rp∗, initial velocity components vc and vρ are expressed as follows:
(47)vc=cμprp∗rimpsinβe;vρ=μp⋅1−cosβesinβe.In Eq. (47), the product of the two velocity components in the skewed axes is expressed as follows:
(48)vc⋅vρ=μc2rp∗rimpcos2βe/2.Sincerimp and rp∗ are precalculated, and the parameters βe and c are related variables, the two velocity components vc⋅vρ are constant. By putting Eq. (46) into Eq. (48), one can obtain a relationship between the radial and circumferential velocity components as follows:
(49)rimp/rp∗−cosβe⋅vx2+sinβe⋅vxvy−μ1−cosβe/rimp=0.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 vorb.imp and flightpath angle ϑl as follows:
(50)rimp/rp∗−cosβe⋅cos2ϑl+sinβe⋅sinϑlcosϑl−μ⋅1−cosβe/vorb.imp2⋅rimp=0,where the flight path angle ϑl is used for parameterizing and updating. The main advantage of setting ϑl as an iteration parameter over other variables lies in the fact that it is bounded [25, 26]. In addition to the natural bound (−π∕2<ϑl<π∕2), additional bounds for the minimum and maximum values for flight path angle ϑl 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 c vector. The minimum flight path angle ϑl.min is expressed by
(51)tanϑl.min=1c⋅1rimp1c⋅1x=rp∗cosβe−rimp/rcrp∗cosβe−π/2/rc=cosβe−rimp/rp∗/sinβe.The maximum flightpath angleϑl.max is obtained from the escape condition and is expressed as follows:
(52)tanϑl.max=sinβe+21−cosβe⋅rimp/rp∗1−cosβe.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ϑl is meaningful, the minimum value of the required velocity vorb.imp satisfying the terminal constraints is as follows:
(53)vorb.imp.min=2μ1−cosβerimpsinβe⋅tanϑl.min.According to Eqs. (51) and (52), the upper and lower bounds of iteration variable ϑl can be obtained from the geocentric distance of intersection point rimp in Eq. (23) and the desired instantaneous impact point rp∗ 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 vorb.imp 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 pointPimp 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 rimp, the ignition time tig 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 βe. The desired instantaneous impact point is given by Eq. (11), and the range angle βe can be obtained by Eq. (9), in which the desired longitude λp and latitude ϕp are predetermined by the mission.Step III. Calculate the bounds for the flight path angle ϑl.min and ϑl.max. According to the results calculated in steps I and II, the flight path angle ϑl.min and ϑl.max 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 vorb.imp.min satisfying the terminal constraints is given by Eq. (53). If vsub.imp2+WM2≥vorb.imp.min, go to step V; otherwise, the step stops, and the output flight path angle is ϑl.min.Step V. Update the iteration variable ϑl. 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 vorb.imp. According to the flight path angle ϑl calculated in step IV and the velocity modulus in Eq. (40), the required velocity vector vorb.imp is obtained by Eq. (39).Step VII. Calculate the target orbit elements. The terminal semimajor axis aorb, eccentricity eorb, and inclination angle iorb are calculated from Eqs. (41), (42), and (43) according to the parameters rimp and vorb.imp, respectively.Step VIII. Calculate the rocket body longitudinal direction xb. The direction of the thrust vector can be calculated by Eq. (29) and xb=Γ.Step IX. Go to step V.


## 4.1. Target Orbit Calculation

According to Eq. (7), an instantaneous impact point 1rp.sub at the equivalent point Pimp with the velocity vector vsub.imp can be obtained, and instantaneous impact point 1rp.orb at the equivalent point Pimp with the velocity vector vorb.imp can be expressed as follows:
(36)1rp.orb=cosβe+ϑlcosϑl1rimp+sinβecosϑl1vorb.imp.Then, according to the equivalent relation in the principle of PA, an instantaneous impact point1rp.sub and instantaneous impact point 1rp.orb are both at the equivalent point Pimp. Therefore, in order to satisfy the terminal constraint Eq. (7), the required velocity vector vorb.imp needs to be solved by
(37)1rp∗=cosβe+ϑlcosϑl1rimp+sinβecosϑl1vorb.imp,where theβe range angle can be obtained by
(38)βe=arcsinrimp×rp∗/rimprp∗.And the required velocity vectorvorb.imp can be defined as follows:
(39)vorb.imp=vx⋅1x+vy⋅1y,=vorb.imp⋅cosϑl⋅1x+sinϑl⋅1y.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 vorb.imp can be calculated:
(40)vorb.imp=vsub.imp2+WM2.The velocity vectorvsub.imp and the position vector rimp are calculated by the navigation system, and the target impact point 1rp∗ and the velocity increment WM are preinput parameters. Thus, the required velocity vector vorb.imp can be obtained by Eqs. (39) and (40), in which the flight path angle ϑl should be meet Eq. (37) constraint. After obtaining the required velocity vector vorb.imp, the target orbital elements of PA can be expressed in the geocentric inertial coordinate system as follows:
(41)aorb=−μrimp/2μ−rimp⋅vorb.imp2,(42)eorb=vorb.imp2−μ/rimprimp−rimp⋅vorb.impvorb.imp/μ,(43)iorb=arccos1z3,where function⋅3 takes the third component of a vector; aorb, eorb, and iorb 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 vectorvorb.imp 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.Figure 5
Calculation required velocity vector according to different velocity components.The required velocity vectorvorb.imp can be expressed in the target flight plane as well without loss of generality as follows:
(44)vorb.imp=vx⋅1x+vy⋅1y=vc⋅1c+vρ⋅1y,where 1c is a chordal unit vector defined from the difference of the two position vectors as follows:
(45)rc=rp∗−rimp=rc⋅1c.Relationships among the velocity components in perpendicular axes (vx, vy) and skewed axes (vc, vρ) are expressed as the following:
(46)vc=rc⋅vxrp∗⋅sinβe;vρ=vy−rp∗⋅cosβe−rimprp∗sinβe⋅vx.It is known that for the family of orbits that connect the two position vectorsrimp and rp∗, initial velocity components vc and vρ are expressed as follows:
(47)vc=cμprp∗rimpsinβe;vρ=μp⋅1−cosβesinβe.In Eq. (47), the product of the two velocity components in the skewed axes is expressed as follows:
(48)vc⋅vρ=μc2rp∗rimpcos2βe/2.Sincerimp and rp∗ are precalculated, and the parameters βe and c are related variables, the two velocity components vc⋅vρ are constant. By putting Eq. (46) into Eq. (48), one can obtain a relationship between the radial and circumferential velocity components as follows:
(49)rimp/rp∗−cosβe⋅vx2+sinβe⋅vxvy−μ1−cosβe/rimp=0.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 vorb.imp and flightpath angle ϑl as follows:
(50)rimp/rp∗−cosβe⋅cos2ϑl+sinβe⋅sinϑlcosϑl−μ⋅1−cosβe/vorb.imp2⋅rimp=0,where the flight path angle ϑl is used for parameterizing and updating. The main advantage of setting ϑl as an iteration parameter over other variables lies in the fact that it is bounded [25, 26]. In addition to the natural bound (−π∕2<ϑl<π∕2), additional bounds for the minimum and maximum values for flight path angle ϑl 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 c vector. The minimum flight path angle ϑl.min is expressed by
(51)tanϑl.min=1c⋅1rimp1c⋅1x=rp∗cosβe−rimp/rcrp∗cosβe−π/2/rc=cosβe−rimp/rp∗/sinβe.The maximum flightpath angleϑl.max is obtained from the escape condition and is expressed as follows:
(52)tanϑl.max=sinβe+21−cosβe⋅rimp/rp∗1−cosβe.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ϑl is meaningful, the minimum value of the required velocity vorb.imp satisfying the terminal constraints is as follows:
(53)vorb.imp.min=2μ1−cosβerimpsinβe⋅tanϑl.min.According to Eqs. (51) and (52), the upper and lower bounds of iteration variable ϑl can be obtained from the geocentric distance of intersection point rimp in Eq. (23) and the desired instantaneous impact point rp∗ 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 vorb.imp 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 pointPimp 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 rimp, the ignition time tig 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 βe. The desired instantaneous impact point is given by Eq. (11), and the range angle βe can be obtained by Eq. (9), in which the desired longitude λp and latitude ϕp are predetermined by the mission.Step III. Calculate the bounds for the flight path angle ϑl.min and ϑl.max. According to the results calculated in steps I and II, the flight path angle ϑl.min and ϑl.max 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 vorb.imp.min satisfying the terminal constraints is given by Eq. (53). If vsub.imp2+WM2≥vorb.imp.min, go to step V; otherwise, the step stops, and the output flight path angle is ϑl.min.Step V. Update the iteration variable ϑl. 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 vorb.imp. According to the flight path angle ϑl calculated in step IV and the velocity modulus in Eq. (40), the required velocity vector vorb.imp is obtained by Eq. (39).Step VII. Calculate the target orbit elements. The terminal semimajor axis aorb, eccentricity eorb, and inclination angle iorb are calculated from Eqs. (41), (42), and (43) according to the parameters rimp and vorb.imp, respectively.Step VIII. Calculate the rocket body longitudinal direction xb. The direction of the thrust vector can be calculated by Eq. (29) and xb=Γ.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.Figure 6
Profiles of the state parameters compared IIP and GEM with the SEM method. (a) The pitch angle vs. time. (b) The altitude vs. range.
(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 Table1.Table 1
Parameter deviation and uncertainty conditions.
Parameter errorsUnit1st stage2nd stageTotal massKg±100±20Engine working times±10±3.0Average thrust of engine%±12±5Aerodynamic coefficients%±15—Atmospheric density%±15—WindStationary wind and random directionIn 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 Figure7.Figure 7
Profiles of the state parameters of 250 Monte Carlo runs for solid rockets. (a) The altitude. (b) The velocity. (c) The flight path angle. (d) The guidance commands.
(a)(b)(c)(d)Figure7 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.Figure8 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.Figure 8
Profiles of the state parameters of 750 Monte Carlo runs under different missions. (a) Altitude. (b) Flight path angle.
(a)(b)Figure 9
The guidance commands the history with deviation and uncertainty under different missions.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.Table 2
Statistics results of 750 Monte Carlo runs for solid rockets under different missions.
Parameter errorsUnitMeanStdv1500 kmRange errorKm15.8425.08Longitude error°/10-11.592.58Latitude error°/10-1-0.660.943000 kmRange errorKm8.7642.22Longitude error°/10-10.493.96Latitude error°/10-1-0.851.734500 kmRange errorKm1.1742.31Longitude error°/10-16.063.51Latitude error°/10-1-3.822.04Figure 10
Drogue deployment accuracy under the proposed guidance method.


## 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.Figure 6
Profiles of the state parameters compared IIP and GEM with the SEM method. (a) The pitch angle vs. time. (b) The altitude vs. range.
(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 Table1.Table 1
Parameter deviation and uncertainty conditions.
Parameter errorsUnit1st stage2nd stageTotal massKg±100±20Engine working times±10±3.0Average thrust of engine%±12±5Aerodynamic coefficients%±15—Atmospheric density%±15—WindStationary wind and random directionIn 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 Figure7.Figure 7
Profiles of the state parameters of 250 Monte Carlo runs for solid rockets. (a) The altitude. (b) The velocity. (c) The flight path angle. (d) The guidance commands.
(a)(b)(c)(d)Figure7 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.Figure8 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.Figure 8
Profiles of the state parameters of 750 Monte Carlo runs under different missions. (a) Altitude. (b) Flight path angle.
(a)(b)Figure 9
The guidance commands the history with deviation and uncertainty under different missions.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.Table 2
Statistics results of 750 Monte Carlo runs for solid rockets under different missions.
Parameter errorsUnitMeanStdv1500 kmRange errorKm15.8425.08Longitude error°/10-11.592.58Latitude error°/10-1-0.660.943000 kmRange errorKm8.7642.22Longitude error°/10-10.493.96Latitude error°/10-1-0.851.734500 kmRange errorKm1.1742.31Longitude error°/10-16.063.51Latitude error°/10-1-3.822.04Figure 10
Drogue deployment accuracy under the proposed guidance method.


## 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.



---

*Source: 5571877-2021-05-31.xml*