A Predictive Explicit Guidance Scheme for Ballistic Missiles

A new approach to the design of ballistic missile guidance is presented in this paper. The proposed method uses the missile model to predict the likely impact point at every guidance cycle and apply course corrections based on the predicted impact point (PIP) deviations. The algorithm also estimates the in-flight thrust variation from nominal and accordingly updates the model to reduce the uncertainty in the prediction of the impact point. The performance of the algorithm is tested through 6-DOF simulation. The simulation results show excellent performance of the proposed guidance scheme in nominal & off nominal cases.


Keywords:  Ballistic , point mass, 6-DOF, down range, cross range,

rRadial position vector
VMissile velocity
γpFlight path angle in elevation plane
γyFlight path angle in azimuth plane
µGeocentric longitude
λGeocentric latitude
mVehicle mass
tburnoutPropulsion burnout time
TMissile thrust
DDrag
Cd, CAerodynamic coefficients
ρDensity
SEffective surface area
ωEarth angular speed
Ψ,Φ,θ Euler angles
α, β Angle of attack in elevation and azimuth  planes
∆DR Impact point down range error
∆CR Impact point cross range error

The role of a ballistic missile is to deliver one or more warheads to a predetermined stationary target.
The design and development of ballistic missile and its guidance has a long history. The problem of ballistic missile guidance has been discussed in numerous literatures such  by Siouris1 and Zarchan2.


A conventional guidance algorithm aiming for required height, velocity and flight path angles combination at the burnout to reach the desired impact point was presented by White3. A nominal trajectory following implicit guidance scheme popular in early days of ballistic missile guidance development was discussed by Schultz4, et al. In this scheme, the missile is guided to follow the nominal trajectory generated on ground and loaded to on-board computer (OBC) before missile lift-off. The guidance algorithm computes the deviation of the current missile position from the nominal position at each instant and commands the required course correction. The disadvantage in this method is that it requires high lateral accelerations pull to correct guidance errors in presence of disturbances such as wind, gust, etc.


In the present scenario, most of the ballistic missiles use explicit guidance scheme in which complete set of trajectory equations are solved on-board and the desired burnout conditions are obtained to hit the impact point. In5 an optimal explicit guidance scheme for a satellite launch vehicle was presented.


In this paper, a new approach has been presented which uses a missile prediction model residing in OBC to predict the likely impact point. The guidance works on the predicted impact point dispersions and brings the missile to optimum burnout states so that the subsequent ballistic phase ensures the missile impact at the desired impact point.


Most of the guidance schemes for ballistic missiles have been developed for the missiles having majority of flight duration out of sensible atmosphere where the effect of drag is negligible. The proposed guidance scheme is advantageous as it works well for short-range ballistic missiles as well where the missile spends significant duration within sensible atmosphere. This is possible since the missile prediction model provides flexibility to incorporate all the necessary data affecting vehicle dynamics. The high computational requirement of the proposed guidance scheme is met by the present day high speed processors.  The details of the algorithm, missile model, guidance command computations and simulation results are discussed in subsequent sections.


The guidance scheme assumes the input state vectors of the missile in the earth centered earth fixed frame (ECEF). The north-east-down (NED) frame acts as the guidance reference frame as the key guidance variables are defined in this frame. The diagram of ECEF co-ordinate frame along with NED frame is shown in Fig.1. The ECEF frame has its origin at the center of the earth and rotates with the earth. The X-axis lies in the equatorial plane going through the intersection of Greenwich meridian and the equator, the Z-axis goes through the North Pole and the Y-axis completes the right handed triad. The Local NED frame has its center at the missile’s center of gravity with its X-axis pointing towards local north, Y-axis pointing towards local east and the Z-axis pointing in the local vertical direction. The velocity vector and the flight path angles used as variables in missile model are shown in Fig. 2.



Figure 1. Reference frames- ECEF and local NED.



Figure 2. Flight path angles in elevation and azimuth plane.


The objective of the guidance algorithm is to guide the missile to reach the specified impact point on the ground. The guidance works on the predicted impact point (PIP) and gives the guidance corrections as a function of predicted terminal error.


Let

X m = f 1 ( X m ,u(t),t) MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipq0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaWaaCbiaeaaca WGybWaaSbaaSqaaiaad2gaaeqaaaqabeaacqGHIaYTaaGccqGH9aqp caWGMbWaaSbaaSqaaiaaigdaaeqaaOGaaiikaiaadIfadaWgaaWcba GaamyBaaqabaGccaGGSaGaamyDaiaacIcacaWG0bGaaiykaiaacYca caWG0bGaaiykaaaa@4580@ (1)

X m = f 2 ( X m ,u(t),h) MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipq0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaWaaCbiaeaaca WGybWaaSbaaSqaaiaad2gaaeqaaaqabeaacqGHIaYTaaGccqGH9aqp caWGMbWaaSbaaSqaaiaaikdaaeqaaOGaaiikaiaadIfadaWgaaWcba GaamyBaaqabaGccaGGSaGaamyDaiaacIcacaWG0bGaaiykaiaacYca caWGObGaaiykaaaa@4575@ (2)

be the equations of motion of the missile with time and altitude as the independent variables respectively. If the initial conditions, the independent guidance vector ‘u’ and the missile state vector Xm (from Onboard INS) are known at every guidance cycle, the equations of motion can be integrated forward to find the probable impact point. In other words,


X m f = t c t f 1 ( X m , u ( t ) , t ) d t + h min 0 f 2 ( X m , u ( t ) , h ) d h MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipq0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaamiwamaaBa aaleaacaWGTbGaamOzaaqabaGccqGH9aqpdaWdXbqaaiaadAgadaWg aaWcbaGaaGymaaqabaaabaGaamiDamaaBaaameaacaWGJbaabeaaaS qaaiaadshaa0Gaey4kIipakiaacIcacaWGybWaaSbaaSqaaiaad2ga aeqaaOGaaiilaiaadwhacaGGOaGaamiDaiaacMcacaGGSaGaamiDai aacMcacaWGKbGaamiDaiabgUcaRmaapehabaGaamOzamaaBaaaleaa caaIYaaabeaaaeaacaWGObWaaSbaaWqaaiGac2gacaGGPbGaaiOBaa qabaaaleaacaaIWaaaniabgUIiYdGccaGGOaGaamiwamaaBaaaleaa caWGTbaabeaakiaacYcacaWG1bGaaiikaiaadshacaGGPaGaaiilai aadIgacaGGPaGaamizaiaadIgaaaa@6095@ (3)                                        

The time-based integration is done from current time until the predicted height becomes less than prefixed altitude in the descent phase of the trajectory. The integration of trajectory equations continue with h as the independent variable from the prefixed altitude in the descent phase until the predicted height becomes zero. This is to ensure that the prediction continues until impact on the ground. The propagation using both ‘time’ and ‘height’ as independent variables is advantageous in terms of computation time as propagation using ‘time’ alone requires very small integration step size at the end.


The desired impact point Xt is known before the launch of the missile. Therefore, the predicted impact point deviation can be calculated by taking the difference between predicted missile impact point and the desired impact point i.e., Xt- Xmf. The predicted impact point error is resolved to get the down range and cross range errors (∆DR, ∆CR). Finally, the guidance variables are updated based on predicted errors to minimize the impact point deviation. The detailed block diagram of the guidance algorithm is shown below in Fig. 3 for better understanding of the flow of the algorithm.
The success of the algorithm mainly depends on the missile model, choice of the guidance variables and the law for updating the guidance variables in minimizing the impact point deviations.



Figure 3. Block diagram of guidance algorithm.


3.1 Guidance Variables for Minimising Impact Point Deviations



The guidance variables chosen are  ( γ p , γ p ) MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipq0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaaiikamaaxa cabaGaeq4SdC2aaSbaaSqaaiaadchaaeqaaaqabeaacqGHIaYTaaGc caGGSaWaaCbiaeaacqaHZoWzdaWgaaWcbaGaamiCaaqabaaabeqaai abgkci3caakiaacMcaaaa@4106@ , where γ p MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipq0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaWaaCbiaeaacq aHZoWzdaWgaaWcbaGaamiCaaqabaaabeqaaiabgkci3caaaaa@3A5E@ γ y MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipq0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaWaaCbiaeaacq aHZoWzdaWgaaWcbaGaamyEaaqabaaabeqaaiabgkci3caaaaa@3A67@ ,are the flight path angle rates of the missile in elevation and azimuth plane respectively. The flight path rates are the basic trajectory parameters, which can be altered for lateral correction by the guidance system to achieve the desired impact point.


For missiles with liquid propulsion, thrust shut off time can be used as another guidance variable whereas the same is not possible for missiles with solid propulsion and hence the thrust shut off time is not used as the guidance variable. However, the burnout time is essential for the proposed predictive guidance as the missile model requires it for propagation. The burnout time is initialized with nominal burnout time and is updated in-flight based on the propulsion performance variation and frozen before the start of closed loop guidance.


The in-flight variation of thrust is accounted in the prediction model by re-scaling the thrust ~ time curve, computed based on the sensed acceleration measured over the nominal acceleration assuming total impulse of the solid motor constant. The problem of ballistic missile guidance assuming constant total impulse of the solid rocket motor is solved in6. In the realistic scenario, the constant impulse of the solid propulsion system is not ensured therefore a velocity trimming guidance scheme is required as presented7
.


3.2 Guidance Command Calculation

The guidance assumes the proper initialization of the variables ( γ p , γ p ) MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipq0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaaiikamaaxa cabaGaeq4SdC2aaSbaaSqaaiaadchaaeqaaaqabeaacqGHIaYTaaGc caGGSaWaaCbiaeaacqaHZoWzdaWgaaWcbaGaamiCaaqabaaabeqaai abgkci3caakiaacMcaaaa@4106@ . This initialization is done on ground before the launch of the missile. The missile model for the first time uses these initialized values and subsequently uses the previous guidance cycle commanded values to calculate PIP. The following procedure is carried out at every guidance cycle for minimizing the PIP deviations.


From the Eqn. (3) and the missile prediction model, it is clear that the final state is in polar form (r, λ, μ). The final state is computed in ECEF Cartesian form, which is the PIP in ECEF frame. This PIP then transformed into launch point east-north-vertical (ENV) frame and subsequently into down range-cross range-altitude (DCH) frame. The down range axis is the line joining launch point (LP) and target point (TP). The cross range axis is the line perpendicular to down range axis. Similar sequence of transformations is applied on the desired ECEF state of the target point to get the desired impact point in DCH frame i.e.,


[ X mf ] DRCR = [T] lpenv drcr . [T] ecef lpenv . [ X mf ] ecef MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipq0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaWaamWaaeaaca WGybWaaSbaaSqaaiaad2gacaWGMbaabeaaaOGaay5waiaaw2faamaa BaaaleaacaWGebGaamOuaiaadoeacaWGsbaabeaakiabg2da9iaacU facaWGubGaaiyxamaaDaaaleaacaWGSbGaamiCaiaadwgacaWGUbGa amODaaqaaiaadsgacaWGYbGaam4yaiaadkhaaaGccaGGUaGaai4wai aadsfacaGGDbWaa0baaSqaaiaadwgacaWGJbGaamyzaiaadAgaaeaa caWGSbGaamiCaiaadwgacaWGUbGaamODaaaakiaac6cacaGGBbGaam iwamaaBaaaleaacaWGTbGaamOzaaqabaGccaGGDbWaaSbaaSqaaiaa dwgacaWGJbGaamyzaiaadAgaaeqaaaaa@5F99@

[ X t ] DRCR = [T] lpenv drcr . [T] ecef lpenv . [ X t ] ecef MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipq0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaWaamWaaeaaca WGybWaaSbaaSqaaiaadshaaeqaaaGccaGLBbGaayzxaaWaaSbaaSqa aiaadseacaWGsbGaam4qaiaadkfaaeqaaOGaeyypa0Jaai4waiaads facaGGDbWaa0baaSqaaiaadYgacaWGWbGaamyzaiaad6gacaWG2baa baGaamizaiaadkhacaWGJbGaamOCaaaakiaac6cacaGGBbGaamivai aac2fadaqhaaWcbaGaamyzaiaadogacaWGLbGaamOzaaqaaiaadYga caWGWbGaamyzaiaad6gacaWG2baaaOGaaiOlaiaacUfacaWGybWaaS baaSqaaiaadshaaeqaaOGaaiyxamaaBaaaleaacaWGLbGaam4yaiaa dwgacaWGMbaabeaaaaa@5DD1@

[ E ] DRCR = [ X t ] drcr [ X mf ] drcr MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipq0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaWaamWaaeaaca WGfbaacaGLBbGaayzxaaWaaSbaaSqaaiaadseacaWGsbGaam4qaiaa dkfaaeqaaOGaeyypa0Jaai4waiaadIfadaWgaaWcbaGaamiDaaqaba GccaGGDbWaaSbaaSqaaiaadsgacaWGYbGaam4yaiaadkhaaeqaaOGa eyOeI0Iaai4waiaadIfadaWgaaWcbaGaamyBaiaadAgaaeqaaOGaai yxamaaBaaaleaacaWGKbGaamOCaiaadogacaWGYbaabeaaaaa@4E53@

∆DR = E(1)  ∆CR = E(2)


The down range error (∆DR) and the cross range error (∆CR) are the functions of the missile burnout parameters in terms of its position, velocity and the flight path angles in elevation and azimuth planes. The missile prediction model estimates the burnout states such as position, velocity and flight path angles. The flight path angle change required at burnout to hit the desired impact point taking estimated burnout position and velocity are as follows


∆γp = (∂γp/∂DR)* ∆DR (4)

   ∆γy = (∂ γy/∂CR)*∆CR (5)

 

where (∂γp/∂DR) and (∂γy/∂CR) are the sensitivity coefficients determining the per unit change requirement in flight path angle for unit change in down range and cross range respectively. ∆DR and ∆CR are the impact point down range and cross range errors as computed by the missile prediction model.


The sensitivity coefficients are computed at every guidance cycle by perturbing flight path angles at burnout and evaluating its effect on impact down range and cross range. The angular flight path angle change requirements ∆γp and ∆γy are distributed from current time to estimated propulsion burnout time (tburnout). The guidance demanded flight path rates are computed as,


γ p . (i+1)= γ . p (i) + Δ γ p ( t burnout t ) MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipq0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGceaqabeaadaWfGa qaamaaxacabaGaeq4SdC2aaSbaaeaacaWGWbaabeaaaSqabeaacaGG UaaaaaGcbeqaaaaacaGGOaGaamyAaiabgUcaRiaaigdacaGGPaGaey ypa0ZaaCbiaeaadaWfGaqaaiabeo7aNbWcbeqaaiaac6caaaGcdaWg aaqaaiaadchaaeqaaiaacIcacaWGPbGaaiykaaqabeaaaaGaey4kaS YaaSaaaeaacqqHuoarcqaHZoWzdaWgaaqaaiaadchaaeqaaaqaamaa bmaabaGaamiDamaaBaaaleaacaWGIbGaamyDaiaadkhacaWGUbGaam 4BaiaadwhacaWG0baabeaakiabgkHiTiaadshaaiaawIcacaGLPaaa aaaabaaaaaa@5570@ (6)

γ y . (i+1)= γ . y (i) + Δ γ y ( t burnout t ) MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipq0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaWaaCbiaeaada WfGaqaaiabeo7aNnaaBaaabaGaamyEaaqabaaaleqabaGaaiOlaaaa aOqabeaaaaGaaiikaiaadMgacqGHRaWkcaaIXaGaaiykaiabg2da9m aaxacabaWaaCbiaeaacqaHZoWzaSqabeaacaGGUaaaaOWaaSbaaeaa caWG5baabeaacaGGOaGaamyAaiaacMcaaeqabaaaaiabgUcaRmaala aabaGaeuiLdqKaeq4SdC2aaSbaaeaacaWG5baabeaaaeaadaqadaqa aiaadshadaWgaaWcbaGaamOyaiaadwhacaWGYbGaamOBaiaad+gaca WG1bGaamiDaaqabaGccqGHsislcaWG0baacaGLOaGaayzkaaaaaaaa @5584@ (7)

                                       

The guidance demanded flight path rate is to be achieved by axial thrust and aerodynamics. The flight path angle rates are translated into an equivalent attitude command and communicated to the attitude autopilot for execution. It can be seen that lateral acceleration, flight path angle rate and the angle of attack are all closely related. The relation among them is given below.


α= mV γ p +gcos γ p T+0.5ρ V 2 S C Nα MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipq0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaeqySdeMaey ypa0ZaaSaaaeaacaWGTbGaamOvamaaxacabaGaeq4SdCgaleqabaGa eyOiGClaaOWaaSbaaSqaaiaadchaaeqaaOGaey4kaSIaam4zaiGaco gacaGGVbGaai4Caiabeo7aNnaaBaaaleaacaWGWbaabeaaaOqaaiaa dsfacqGHRaWkcaaIWaGaaiOlaiaaiwdacqaHbpGCcaWGwbWaaWbaaS qabeaacaaIYaaaaOGaam4uaiaadoeadaWgaaWcbaGaamOtaiabeg7a Hbqabaaaaaaa@5224@ (8)

  β= mVcos γ p γ y T+0.5ρ V 2 S C Nα MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipq0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaeqOSdiMaey ypa0ZaaSaaaeaacaWGTbGaamOvaiGacogacaGGVbGaai4Caiabeo7a NnaaBaaaleaacaWGWbaabeaakmaaxacabaGaeq4SdC2aaSbaaSqaai aadMhaaeqaaaqabeaacqGHIaYTaaaakeaacaWGubGaey4kaSIaaGim aiaac6cacaaI1aGaeqyWdiNaamOvamaaCaaaleqabaGaaGOmaaaaki aadofacaWGdbWaaSbaaSqaaiaad6eacqaHXoqyaeqaaaaaaaa@504C@ (9)

Now, the attitude command to the autopilot to achieve (or equivalently) can be communicated.  Attitude command in terms of euler angles from NED to body is,


ψ= π 2 β γ y ;φ=π;θ=( γ p +α ) MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipq0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaeqiYdKNaey ypa0ZaaSaaaeaacqaHapaCaeaacaaIYaaaaiabgkHiTiabek7aIjab gkHiTiabeo7aNnaaBaaaleaacaWG5baabeaakiaacUdacqaHgpGAcq GH9aqpcqaHapaCcaGG7aGaeqiUdeNaeyypa0JaeyOeI0YaaeWaaeaa cqaHZoWzdaWgaaWcbaGaamiCaaqabaGccqGHRaWkcqaHXoqyaiaawI cacaGLPaaaaaa@5209@ (10)  

  

The missile is modeled as a point mass traveling over a spherical rotating earth. The propulsive, aerodynamic and gravity forces are modeled. The detailed derivation of trajectory equations is discussed by Vinh8. The trajectory equations are presented by Song & Tahk9 as a function of time is given below


    v ˙ =k. ( TcosαcosβD ) m gsin γ p MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipq0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGabmODayaaca Gaeyypa0Jaam4Aaiaac6cadaWcaaqaamaabmaabaGaamivaiGacoga caGGVbGaai4Caiabeg7aHjGacogacaGGVbGaai4Caiabek7aIjabgk HiTiaadseaaiaawIcacaGLPaaaaeaacaWGTbaaaiabgkHiTiaadEga ciGGZbGaaiyAaiaac6gacqaHZoWzdaWgaaWcbaGaamiCaaqabaaaaa@4EF8@     (11)


γ . p = ( γ . p ) demand MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipq0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaWaaCbiaeaacq aHZoWzaSqabeaacaGGUaaaaOWaaSbaaSqaaiaadchaaeqaaOGaeyyp a0JaaiikamaaxacabaGaeq4SdCgaleqabaGaaiOlaaaakmaaBaaale aacaWGWbaabeaakiaacMcadaWgaaWcbaGaamizaiaadwgacaWGTbGa amyyaiaad6gacaWGKbaabeaaaaa@4593@          if  (t < tburnout) (12)


γ ˙ p = gcos γ p v + v r cos γ p + ω 2 r v cosλ( cos γ p cosλ+sin γ p sin γ y sinλ ) +2ωcos γ y cosλ MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipq0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGceaqabeaacuaHZo WzgaGaamaaBaaaleaacaWGWbaabeaakiabg2da9maalaaabaGaeyOe I0Iaam4zaiGacogacaGGVbGaai4Caiabeo7aNnaaBaaaleaacaWGWb aabeaaaOqaaiaadAhaaaGaey4kaSYaaSaaaeaacaWG2baabaGaamOC aaaaciGGJbGaai4BaiaacohacqaHZoWzdaWgaaWcbaGaamiCaaqaba aakeaacqGHRaWkdaWcaaqaaiabeM8a3naaCaaaleqabaGaaGOmaaaa kiaadkhaaeaacaWG2baaaiGacogacaGGVbGaai4CaiabeU7aSnaabm aabaGaci4yaiaac+gacaGGZbGaeq4SdC2aaSbaaSqaaiaadchaaeqa aOGaci4yaiaac+gacaGGZbGaeq4UdWMaey4kaSIaci4CaiaacMgaca GGUbGaeq4SdC2aaSbaaSqaaiaadchaaeqaaOGaci4CaiaacMgacaGG UbGaeq4SdC2aaSbaaSqaaiaadMhaaeqaaOGaci4CaiaacMgacaGGUb Gaeq4UdWgacaGLOaGaayzkaaaabaGaey4kaSIaaGOmaiabeM8a3jGa cogacaGGVbGaai4Caiabeo7aNnaaBaaaleaacaWG5baabeaakiGaco gacaGGVbGaai4CaiabeU7aSbaaaa@7F19@          (t>= tburnout) (13)


γ . y = ( γ ˙ y ) demand    MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipq0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaWaaCbiaeaacq aHZoWzaSqabeaacaGGUaaaaOWaaSbaaSqaaiaadMhaaeqaaOGaeyyp a0ZaaeWaaeaacuaHZoWzgaGaamaaBaaaleaacaWG5baabeaaaOGaay jkaiaawMcaamaaBaaaleaacaqGKbGaaeyzaiaab2gacaqGHbGaaeOB aiaabsgaaeqaaOGaaeiiaiaabccaaaa@461D@          if (t < tburnout) (14)


γ ˙ y = vcos γ p cos γ y tanλ r +2ω( tan γ p sin γ y cosλsinλ ) ω 2 r Vcos γ p cos γ y sinλcosλ MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipq0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGceaqabeaacuaHZo WzgaGaamaaBaaaleaacaWG5baabeaakiabg2da9iabgkHiTiaadcca caWGGaWaaSaaaeaacaWG2bGaci4yaiaac+gacaGGZbGaeq4SdC2aaS baaSqaaiaadchaaeqaaOGaci4yaiaac+gacaGGZbGaeq4SdC2aaSba aSqaaiaadMhaaeqaaOGaciiDaiaacggacaGGUbGaeq4UdWgabaGaam OCaaaacqGHRaWkcaaIYaGaeqyYdC3aaeWaaeaaciGG0bGaaiyyaiaa c6gacqaHZoWzdaWgaaWcbaGaamiCaaqabaGcciGGZbGaaiyAaiaac6 gacqaHZoWzdaWgaaWcbaGaamyEaaqabaGcciGGJbGaai4Baiaacoha cqaH7oaBcqGHsislciGGZbGaaiyAaiaac6gacqaH7oaBaiaawIcaca GLPaaaaeaacqGHsisldaWcaaqaaiabeM8a3naaCaaaleqabaGaaGOm aaaakiaadkhaaeaacaWGwbGaci4yaiaac+gacaGGZbGaeq4SdC2aaS baaSqaaiaadchaaeqaaaaakiGacogacaGGVbGaai4Caiabeo7aNnaa BaaaleaacaWG5baabeaakiGacohacaGGPbGaaiOBaiabeU7aSjGaco gacaGGVbGaai4CaiabeU7aSbaaaa@8214@          if(t>= tburnout) (15)


r . =v s in γ p MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipq0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaWaaCbiaeaaca WGYbWaaSbaaSqaaaqabaaabeqaaiaac6caaaGccqGH9aqpcaWG2bWa aSraaSqaaaqabaGcciGGZbGaaiyAaiaac6gacqaHZoWzdaWgaaWcba GaamiCaaqabaaaaa@3FC8@ (16)


  μ . = vcos γ p cos γ y rcosλ MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipq0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaWaaCbiaeaacq aH8oqBaSqabeaacaGGUaaaaOGaeyypa0ZaaSaaaeaacaWG2bGaci4y aiaac+gacaGGZbGaeq4SdC2aaSbaaSqaaiaadchaaeqaaOGaci4yai aac+gacaGGZbGaeq4SdC2aaSbaaSqaaiaadMhaaeqaaaGcbaGaamOC aiGacogacaGGVbGaai4CaiabeU7aSbaaaaa@4B70@ (17)


   λ . = vcos γ p sin γ y r MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipq0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaWaaCbiaeaacq aH7oaBaSqabeaacaGGUaaaaOGaeyypa0ZaaSaaaeaacaWG2bGaci4y aiaac+gacaGGZbGaeq4SdC2aaSbaaSqaaiaadchaaeqaaOGaci4Cai aacMgacaGGUbGaeq4SdC2aaSbaaSqaaiaadMhaaeqaaaGcbaGaamOC aaaaaaa@46EC@ (18)

where


D= 1 2 ρ v 2 s C d MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipq0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaamiraiabg2 da9maalaaabaGaaGymaaqaaiaaikdaaaGaeqyWdiNaamODamaaCaaa leqabaGaaGOmaaaakiaadohacaWGdbWaaSbaaSqaaiaadsgaaeqaaa aa@3FAC@

The aerodynamic drag coefficient Cd is a function of mach no. and height. The air density (ρ) and gravitational acceleration (g) are functions of height. 


The equations of motion need to be integrated numerically to obtain the final missile states. A fourth-order Runge-kutta (RK4) method has been chosen for solving the trajectory equations. The RK4 method is discussed in details by Sastry10.


4.1 Estimation of Net Acceleration Variation


The proposed predictive guidance is sensitive to modeling and data inaccuracies as it will lead to higher predicted impact errors hence higher control effort. Therefore, it is important to identify all the necessary parameters affecting missile in-flight dynamics. The main parameters that can vary in real time are the vehicle propulsion and aerodynamics. The variation in these two parameters is accounted by introducing a scale factor (k) in the Eqn. (11). The scale factor is estimated during missile flight at very guidance cycle by taking the ratio of sensed acceleration obtained from onboard INS and the nominal acceleration. The nominal acceleration is calculated using stored nominal thrust and drag profiles, resolved in the current velocity direction. For instance, a scale factor of 1.07 would indicate a net acceleration (i.e., (Thrust – Drag)/Mass) variation of 7 % over that of the nominal.


Therefore, the equation for the net acceleration along the flight direction will be


V =k.( TD M ) MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipq0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaWaaCbiaeaaca WGwbaaleqabaGaeyOiGClaaOGaeyypa0Jaam4Aaiaac6cadaqadaqa amaalaaabaGaamivaiabgkHiTiaadseaaeaacaWGnbaaaaGaayjkai aawMcaaaaa@4028@

here k is the scale factor.


The proposed guidance algorithm is thoroughly tested in 6-DOF simulation under nominal and off nominal cases. The 6-DOF simulation test bed is written in ‘FORTRAN’ programming language. The nominal missile data used for evaluating the guidance algorithm is shown in Table 1. Off-nominal cases involve perturbation on the vehicle thrust, weight and aerodynamics. The algorithm is also tested under wind conditions. The perturbation cases used for testing the algorithm is mentioned in Table 2. The simulation results show an excellent performance of the algorithm and the impact error is observed well below 50 m in all the cases.


The simulation results of nominal and perturbation cases are discussed below. The guidance predicted errors (down range and cross range errors) in nominal, thrust-up and thrust down cases are shown in Fig.4 and fig.5 and are observed well within a dead zone of 20 m. It takes few guidance cycles (8-10 guidance cycles) to converge into guidance tunnel of specified value. The simulation study is carried out taking guidance cycle of 100 ms. The trajectories in azimuth and elevation plane are shown in Fig.6 and fig.7.


Table 1. Missile data



Table 2. Vehicle data perturbation cases



Figure 4. Predicted down range error at impact.



Figure 5. Predicted cross range error at impact.



5.1 Hardware-in-loop Simulation Results


Hardware-in-loop simulation (HILS) test bed consists of 6-dof model, On board computer (OBC), Missile interface unit (MIU) and Inertial navigation system (INS). The test bed was based on 1553 bus communication protocol. The mission software consisting of guidance algorithm is validated in HILS environment using OBC in loop configuration. It is also validated through using sensor in loop, actuator in loop and sensor-actuator in loop. A typical actuator in loop simulation results are shown below in Fig.8 and fig.9.



Figure 6. Trajectory in azimuth plane.



Figure 7. Trajectory in elevation plane.


A novel explicit guidance scheme has been developed for the ballistic missiles. A single stage flight vehicle with solid propulsion system was used to evaluate the performance of the guidance algorithm in this paper. The guidance algorithm can be used for multiple stage flight vehicles by incorporating necessary data of different stages in the prediction model. The guidance algorithm developed is advantageous especially for short-range ballistic missiles having significant flight within atmosphere. The guidance algorithm has been tested and validated through HILS test bed and all the simulation results were satisfactory.



Figure 8. Missile axial (Ax) and accelerations (Ay and Az).



Figure 9. Missile body rates- roll, pitch and Yaw.


The authors acknowledge Shri N.V. Kadam for his motivation and support during the course of the work. The authors also acknowledge the immense involvement and guidance of Dr S.K. Choudhury during the testing and validation of the algorithms through HILS test bed.


1.     Siouris, George M. Missile guidance and control systems. Springer-Verlag New York, 2004, pp.365-489.

2.     Zarchan, Paul. Tactical and strategic missile guidance. In Progress in astronautics and aeronautics, Edn 2. 1994, 157, pp.9-23, 238-263. [Full text via CrossRef]

3.     White, John E. Guidance and targeting for the strategic target system. J. Guidance, Control Dyn., 1992, 15(6), 1313-1319.
[Full text via CrossRef]

4.     Schultz, P.R.; soufl, Robert V. & Grubin, Carl. A boost guidance scheme for following a given trajectory and satisfying injection constraints. J. Spacecraft, 1966, 3(8),  1209-1215.
[Full text via CrossRef]

5.     Sinha, S.K.; Shrivastava, S.K.; M.S. & Prabhu, K.S.  Optimal explicit guidance for three-dimensional launch trajectory. Acta Astronautica, 1989, 19(2), 115-123.
[Full text via CrossRef]

6.     Padhi, Radhakant. An optimal explicit guidance scheme for ballistic missiles with solid motors. 1999, AIAA-99-4140. [Full text via CrossRef]

7.     Thomas, Tessy. Guidance scheme for solid propelled vehicle during atmospheric phase. Def. Sci. J., 2005,  55(3), 253-264.
[Full text PDF]

8.     Vinh, Nguyen X. Optimal trajectories in atmospheric flight, ESPC. Amsterdam-Oxford-New York, 1981, pp. 47-62.  

9.     Song, Eun-Jung & Tahk, Min-Jea. Suboptimal midcourse guidance for interception of free-fall targets. 1999, AIAA-99-4067.  [Full text via CrossRef]

10.   Sastry, S.S. Introductory methods of numerical analysis. Ed 3. 1998,  pp. 254-256.

, *,

Mr N. Prabhakar obtained his BE (Electrical and Electronics) from Annamalai University and ME (Aeronautics) from Indian Institute of Science, Bangalore, in 1975 and 1978, respectively. Currently working as Scientist ‘H’/Outstanding Scientist at DRDL and Project Director AD (Missions), Programme AD. He received, DRDO Scientist of the Year Award in 2001, the Path Breaking Technology Award in the year 2007 and DRDO Performance Excellence Award in 2009. He has been working in areas of trajectory optimisation, missile system design, modeling and simulation. His main areas of interest are optimisation, distributed simulation, and system design.

Mr Indra Deo Kumar obtained his BTech (Electronics and Communication) from NIT Kurukshetra and MTech (Digital Systems) from MNNIT Allahabad, in 2000 and 2003 respectively. He is working as Scientist ‘D’ at Programme AD, Research Centre Imarat. He received DRDO AGNI Award for Excellence in Self-reliance in 2007 and DRDO Performance Excellence Award in 2009 for his contribution to ‘AD’ Programme and PRITHVI/ DHANUSH respectively. His significant contributions are in the area of missile system, guidance design, trajectory optimisation, modeling and simulation.

Mr Sunil Kumar Tata obtained Msc (Mathematics) from Jawaharlal Nehru Technological University (JNTU) in 1999. He is presently working as Scientist- ‘D’ in DRDL. He received DRDO AGNI Award for Excellence in Self-reliance in 2007 for his contribution to ‘AD’ Programme. His significant contributions are in the areas of missile system, guidance design, 6DOF modeling and simulation.

Dr V. Vaithiyanathan obtained MSc from Bharathidasan University, in 1989; MPhil from Madras University, in 1991; MCA from Bharathidasan University, in 2000 and PhD from Alagappa University, in 2004. Presently he is Associate Dean – Research at School of Computing, SASTRA University. His research areas of interest include: Cryptography, image processing and soft computing techniques.