Adaptive Navigation, Guidance and Control Techniques Applied to Ballistic Projectiles and Rockets

Abstract

Accuracy and precision are the cornerstone for ballistic projectiles from the earliest days of this discipline. In the beginnings, impact point precision in artillery devices deteriorated when range were extended, particularly for non-propelled artillery rockets and shells. Later, inertial navigation and guidance systems are introduced and precision was unlinked from range increases. In the last 30 years, hybridization between inertial systems and GNSS devices has improved precision enormously. Unfortunately, during the last stages of flight, inertial and GNSS methods (hybridized or not) feature big errors on attitude and position determination. Low cost devices, which are precise on terminal guidance and do not feature accumulative error, such as quadrant photo-detector, seem to be appropriate to be included on the guidance systems. Hybrid algorithms, which combine GNSSs, IMUs and photodetectors, and a novel technic of attitude determination, which avoids the use of gyroscopes, are presented in this chapter. Hybridized measurements are implemented on modified proportional navigation law and a rotatory force control method. A realistic non-linear flight dynamics model has been developed to perform simulations to prove the accuracy of the presented algorithms.

1. Introduction

A precision-guided munition (PGM) is a guided munition intended to precisely hit a specific target, and to minimize collateral damage. Considering that the damage effects of explosive weapons decrease with distance, even modest improvements in accuracy enable a target to be attacked with fewer or smaller bombs. The precision of these weapons is dependent both on the precision of the measurement system used for location determination and the precision in setting the coordinates of the target. The latter critically depends on intelligence information, not all of which is accurate. If the targeting information is accurate, satellite-guided weapons (including inertial navigation in the event of signal loss) are significantly more likely to achieve a successful strike in any given weather conditions than any other type of precision-guided munition [1].

Development of low-cost navigation, guidance and control technologies for unguided rockets is a unique engineering challenge. Over the past several decades, numerous solutions have been proposed, primarily for large artillery projectiles or for slowly rolling airframes [1, 2].

GNSS/IMU hybridizing systems provide accurate solutions for PGMs but in some occasions these solutions might not be enough. For those systems, a circle error probable (CEP) is around 10–20 m in the best cases [3].

Development of algorithms for low-cost high-precision terminal guidance systems is a cornerstone in research on PGMs. [4] propose a robust guidance law which is mainly suitable for systems characterized by moderate glint levels such as electro-optical missiles [5] present a new precision guidance law for three-dimensional intercepts. In contrast to previously published guidance laws, it does not require knowledge of the range to the target. This makes it appropriate for use on platforms which have an imaging device, such as a video camera, as a primary sensor.

Whit the aim of reducing costs, many inertial navigation systems could be replaced with less accurate devices if it were guaranteed that GNSS signal would be continuously available to update the inertial system to limit its error growth.

However, given the uncertainty in most scenarios, an alternate way to reduce system cost and collateral damage is to lower the cost by developing lower-cost inertial sensors while improving their accuracy using other sources of navigation information such as Semi Active laser kits. The benefits of integrated data fusion have been demonstrated across the spectrum of antisubmarine, tactical air, and land warfare [6].

In the research described in this chapter, two measured quantities are used to obtain attitude information for high dynamic and spin rate vehicles: speed and gravity vectors. They are obtained in two different reference frames using a GNSS sensor and a strap-down accelerometer. After that, attitude determination algorithm is integrated in the global hybridized system together with SAL and inertial measurements. Non-linear flight simulations are performed to prove the applicability of the proposed approach for ballistic rocket navigation, guidance, and control.

1.1. Semi-active laser photodetectors (SAL)

Semi-active laser photodetectors (SAL), and particularly quadrant detector devices, have been developed to improve precision in guided weapons. Quadrant photo-detectors have been applied in many engineering ambits, such as measurement, control, laser collimation, target tracking, and particularly in PGM terminal guidance [7]. One of the greatest advantages of quadrant detector equipment is the high performance provided in terms of guidance, typically in the last stages of the trajectory, as compared to the low cost incurred. Coordination can be achieved without requiring lengthy transfer of coordinates which is susceptible to errors. But, constant lines of sight between the target, laser designator and the weapon must be maintained [3].

Quadrant photo-detector is a uniform disk made of silicon containing two gaps across the diameters. There are four independent and equal photo-detectors on the sensing surface, one on each quadrant. The centre of the detector is known very accurately since it is the mechanical intersection of the two gap lines and is not pretended to change with time or temperature. A symmetrical laser or other optical beam centered on the detector generates equal currents from each of the four sectors. If the beam moves from the centre of the detector, the currents from the four sectors change, and a processing method may give the coordinate displacements relative to the centre [8]. Precision on determining the coordinates of intersection of the beam with the photo-detector will determine the key points on the Navigation and Guidance algorithms for the terminal phase on a PGM. A wide dissertation on improving precision in this photodetector is presented in [3].

Modern laser guided ballistic rockets are integrating IMU, GPS and laser guidance capability, offering high precision, all-weather attack capability. For example, [9] design a missile target tracker using a filter/correlator based on forward-looking infrared sensor measurements. In this chapter, improvements on existing methods for terminal guidance are presented, which apply an effective hybridization algorithm in order to obtain an accurate vector between rocket and target from a combination of sensors previously mentioned, namely IMU, GPS and SAL.

1.2. Attitude determination technics and integration on the global guidance

Obtaining precise attitude information is essential for navigation and control. Its effectiveness is determined by the degree of precision of navigation and control systems, including inertial measurement units [10]. There is an extensive body of literature regarding attitude estimation using various sensor inputs [11].

Traditionally, in order to obtain accurate values for determining attitude, expensive and/or weighty units, such as laser or fiber optic gyroscopes and accelerometers, or their MEMS equivalents, must be employed. Moreover, when high-demanding maneuvers are performed this equipment may become extremely expensive.

It is well-known that the attitude of an aero-vehicle may be determined, starting from an initial condition, integrating the angular rates (pitch, roll, and yaw rates) of the vehicle and propagating them forward in time. Nevertheless, accuracy requirements usually cannot be satisfied by using inexpensive sensors [10]. This problem becomes even more important when the vehicle cannot be reused: low-cost attitude determination systems are of key importance for these applications.

For example, [12] describe an attitude determination system that is based on two measurements of non-zero, non-co-linear vectors. Using the Earth’s magnetic field and gravity as the two measured quantities, a low-cost attitude determination system is proposed.

[13] develop an inexpensive Attitude Heading Reference System for general aviation applications by fusing low cost automotive grade inertial sensors with GPS. The inertial sensor suit consists of three orthogonally mounted solid-state rate gyros.

[14] describe an attitude estimation algorithm derived by post-processing data from a small low cost Inertial Navigation System recorded during the flight of a sub-scale commercial off the shelf UAV. Estimates of the UAV attitude are based on MEMS gyro, magnetometer, accelerometer, and pitot tube inputs.

[15] state that low-cost GNSS receivers and antennas can provide a precise attitude and drift-free position information, but accuracy is not continuous. Inertial sensors are robust to GNSS signal interruption and very precise over short time frames, which enables a reliable cycle slip correction. But low-cost inertial sensors suffer from a substantial drift. The authors propose a tightly coupled position and attitude determination method for two low-cost GNSS receivers, a gyroscope and an accelerometer and obtain a heading with an accuracy of 0.25° and an absolute position with an accuracy of 1 m.

Similar developments may be found within space vehicles, for example in [16]. In [17] the use of an inertial navigation system (INS) and a multiple GPS antenna system for attitude determination of an off-road vehicle is developed. And in [18], attitude determination using GPS carrier phase is successfully applied to aircraft in experiments.

Also, improved algorithms for estimating attitude in case of failures have been proposed in the literature. For example, [19] introduce algorithms with filter gain correction for the case of measurement malfunctions. Two different algorithms are proposed and applied for the attitude estimation process of a pico-satellite. The results of these algorithms are compared for different types of measurement faults in different estimation scenarios and recommendations about their applications are given.

However, as stated in [20], many of the presented methods, such as the ones employing local magnetic field vectors, are only valid for estimating the orientation of a slow-rotation body: for high spin rate bodies, electromagnetic interactions degrade magnetic measurements.

2. Non-linear flight dynamics model

This segment portrays the nonlinear flight dynamic model utilized as a part of this tests, including actuations, and navigation and control performances.

2.1. Rocket

The guidance and control detailing proposed in this investigation applies to a 140-mm axisymmetric turning rocket with wrap around balancing out blades. It highlights supersonic dispatch speed and a turn rate of roughly 150 Hz. The control system features a roll-decoupled fuse set at the nose of the rocket. This fuse is composed of four canard surfaces, decoupled 2 by 2. keeping in mind the end goal to produce control regulated in modulus and argument, situated in an orthogonal plane in respect to rocket, and its related moment as it is exposed in Figure 1.

Figure 1.

The 140-mm axisymmetric rocket with wrap around fins and a roll-decoupled fuse.

The non-controlled solid propellant thrust, mass, inertia moments ( $I_x$ and $I_y$ ) and centre of gravity coordinate from nose ( $X_{\mathit{CG}}$ ) versus time are shown in Table 1.

Time (s)	0.00	0.02	0.10	0.20	0.70	1.20	1.70	1.75	1.95	2.00	2.15	2.30	2.70	3.20	100.00
Thrust (kN)	0.00	25.00	22.50	23.00	24.00	25.50	28.50	29.16	15.00	10.00	5.00	2.50	0.00	0.00	0.00
Mass (kg)	62.40	62.30	61.52	60.58	55.76	50.67	45.13	44.54	42.72	42.47	42.00	41.77	41.57	41.57	41.57
$I_x\left(m·s^2\right)$	0.19	0.19	0.19	0.19	0.18	0.17	0.16	0.16	0.15	0.15	0.15	0.15	0.15	0.15	0.15
$I_y\left(m·s^2\right)$	18.85	18.83	18.71	18.57	17.80	16.96	16.02	15.91	15.59	15.54	15.46	15.42	15.38	15.38	15.38
$X_{\mathit{CG}}\left(m\right)$	1.13	1.13	1.13	1.13	1.11	1.10	1.07	1.07	1.06	1.06	1.06	1.06	1.06	1.06	1.06

Table 1.

140 mm axisymmetric rocket main parameters versus time.

Numerical simulations were employed to determine aerodynamic coefficients for the rocket under examination, which are showed in Figure 2.

Figure 2.

Aerodynamic coefficients vs. Mach number.

2.2. Coordinate systems definition

Two axis systems are defined along this paper: north east down axes (NED) and body axes (B). NED axes are defined by sub index NED. $x_{\mathit{NED}}$ pointing north, $y_{\mathit{NED}}$ perpendicular to $x_{\mathit{NED}}$ and pointing East, and $z_{\mathit{NED}}$ forming a clockwise trihedron. Body axes are defined by sub index B. $x_B$ pointing forward and contained in the plane of symmetry of the rocket, $z_B$ perpendicular to $x_B$ , pointing down and contained in the plane of symmetry of the rocket, and $y_B$ forming a clockwise trihedral. The origin of body axes is located at the centre of mass of the rocket and they are severely coupled to the roll-decoupled fuse. This concept is shown on Figure 3.

Figure 3.

Reference systems.

2.3. Mathematical equations

Total forces and moments on the rocket are given (expressed in body axes) by (1) and (2), respectively:

$$\overrightarrow{F_B}=\overrightarrow{D_B}+\overrightarrow{L_B}+\overrightarrow{M_B}+\overrightarrow{P_B}+\overrightarrow{T_B}+\overrightarrow{W_B}+\overrightarrow{C_B}+\overrightarrow{C{F}_B},$$ E1

$$\overrightarrow{M_B}=\overrightarrow{O_B}+\overrightarrow{P_{M_B}}+\overrightarrow{M_{M_B}}+\overrightarrow{S_B}+\overrightarrow{C{M}_B}$$ E2

where $\overrightarrow{{\mathrm{D}}_{\mathrm{B}}}$ is the drag force, $\overrightarrow{{\mathrm{L}}_{\mathrm{B}}}$ is the lift force, $\overrightarrow{{\mathrm{M}}_{\mathrm{B}}}$ is the Magnus force, $\overrightarrow{{\mathrm{P}}_{\mathrm{B}}}$ is the pitch damping force, $\overrightarrow{{\mathrm{T}}_{\mathrm{B}}}$ is the thrust force, $\overrightarrow{{\mathrm{W}}_{\mathrm{B}}}$ is the weight force and $\overrightarrow{{\mathrm{C}}_{\mathrm{B}}}$ is the Coriolis force, $\overrightarrow{C{F}_B}$ is the control force executed by the airfoils, $\overrightarrow{{\mathrm{O}}_{\mathrm{B}}}$ is the overturn moment, $\overrightarrow{{\mathrm{P}}_{{\mathrm{M}}_{\mathrm{B}}}}$ is the pitch damping moment, $\overrightarrow{{\mathrm{M}}_{{\mathrm{M}}_{\mathrm{B}}}}$ is the Magnus moment and $\overrightarrow{{\mathrm{S}}_{\mathrm{B}}}$ is the spin damping moment and $\overrightarrow{C{M}_B}$ is the control moment executed by the airfoils. Rocket forces in body axes include contributions from drag, lift, Magnus, pitch damping, thrust, weight and Coriolis forces, which are described by the following expressions:

$$\overrightarrow{D_B}=-\frac{\pi }{8}{d}^2\rho \left(C_{D_0}+C_{D_{{\alpha }^2}}{\alpha }^2\right)‖\overrightarrow{v_B}‖\overrightarrow{v_B},$$ E3

$$\overrightarrow{L_B}=-\frac{\pi }{8}{d}^2\rho \left(C_{L_{\alpha }}·\alpha +C_{L_{{\alpha }^3}}{\alpha }^2\right)\left({‖\overrightarrow{v_B}‖}^2\overrightarrow{x_B}-\left(\overrightarrow{x_B}·\overrightarrow{v_B}\right)\overrightarrow{v_B}\right),$$ E4

$$\overrightarrow{M_B}=-\frac{\pi }{8}{d}^3\rho \frac{C_{\mathit{mf}}}{I_x}\left(\overrightarrow{L_B}·\overrightarrow{x_B}\right)\left(\overrightarrow{x_B}\times \overrightarrow{v_B}\right),$$ E5

$$\overrightarrow{P_B}=\frac{\pi }{8}{d}^3\rho \frac{C_{\mathit{Nq}}}{I_y}{‖\overrightarrow{v_B}‖}^2\left(\overrightarrow{L_B}\times \overrightarrow{x_B}\right),$$ E6

$$\overrightarrow{T_B}=T\left(t\right)\overrightarrow{x_B},$$ E7

$$\overrightarrow{W_B}=m\overrightarrow{g_B},$$ E8

$$\overrightarrow{C_B}=-2m\overrightarrow{\mathrm{\Omega }}\times \overrightarrow{v_B},$$ E9

where $\mathrm{d}$ is the rocket caliber, $\mathrm{\rho }$ is the air density, ${\mathrm{C}}_{{\mathrm{D}}_0}$ is the drag force linear coefficient, ${\mathrm{C}}_{{\mathrm{D}}_{{\mathrm{\alpha }}^2}}$ is the drag force square coefficient, $\mathrm{\alpha }$ is the total angle of attack, ${\mathrm{C}}_{{\mathrm{L}}_{\mathrm{\alpha }}}$ is the lift force linear coefficient, ${\mathrm{C}}_{{\mathrm{L}}_{{\mathrm{\alpha }}^3}}$ is the lift force cubic coefficient, ${\mathrm{C}}_{\mathrm{mf}}$ is the Magnus force coefficient, $\overrightarrow{L_B}$ is the rocket angular momentum expressed in body axes, ${\mathrm{I}}_{\mathrm{x}}\mathit{\text{and}}{\mathrm{I}}_{\mathrm{y}}$ are the rocket inertia moments in body axes, ${\mathrm{C}}_{\mathrm{Nq}}$ is the pitch damping force coefficient, $\overrightarrow{{\mathrm{x}}_{\mathrm{B}}}$ is the rocket nose pointing vector expressed in body axes, $\overrightarrow{{\mathrm{g}}_{\mathrm{B}}}$ is the gravity vector in body axes, $\overrightarrow{\mathrm{\Omega }}$ is the earth angular speed vector, and $\overrightarrow{v_B}$ is the rocket velocity expressed in body axes.

Keeping in mind the end goal to demonstrate the control forces and moments in body reference frame for each of the four fins, it must be viewed as first the effective incidence aerodynamic speed on each of the four control surfaces. The expressions for control force on each of the four control surfaces are characterized in the accompanying equations:

$$\overrightarrow{C{F}_B}=\sum _{i=1}^4\left[\frac{1}{8}\frac{{\alpha }_{E{f}_i}}{‖{\alpha }_{E{f}_i}‖}{d}^2\mathit{\rho \pi }{‖\overrightarrow{v_{\mathit{xE}{f}_i}}‖}^2\left[{C_N}_{\mathit{\alpha w}}\text{cos}{\alpha }_{E{f}_i}+\frac{2}{d^2\pi }{S}_{\mathit{exp}}{\text{sin}}^2{\alpha }_{E{f}_i}\right]\left(\overrightarrow{u_{F{N}_i}}\mathit{cos}{\delta }_i-\overrightarrow{x_B}\mathit{sin}{\delta }_i\right)\right],$$ E10

where ${C_N}_{\mathit{\alpha w}}$ is the aerodynamic coefficient of the normal force for a fin, $S_{\mathit{exp}}$ is the reference surface of the fin, ${\delta }_i$ responds to fin deflection angle, $\overrightarrow{u_{F{N}_i}}$ depends on fin orientation, concretely in body axes, $\overrightarrow{u_{F{N}_1}}=\left[010\right]$ , $\overrightarrow{u_{F{N}_3}}=\left[0-10\right]$ , $\overrightarrow{u_{F{N}_2}}=\left[001\right]$ , $\overrightarrow{u_{F{N}_4}}=\left[00-1\right]$ , and, ${\alpha }_{E{f}_i}$ and $\overrightarrow{v_{\mathit{xE}{f}_i}}$ are the effective angle of attack and the effective aerodynamic speed on each of the four fins, respectively. These last two magnitudes are modeled as it is expressed in the following equations:

$$\overrightarrow{v_{\mathit{xE}{f}_i}}=\overrightarrow{v_B}-\left[\overrightarrow{v_B}·\overrightarrow{u_{b_i}}\right]\overrightarrow{u_{b_i}},$$ E11

$${\alpha }_{E{f}_i}=\mathit{\text{acos}}\left[\frac{\overrightarrow{v_B}-\left[\overrightarrow{v_B}·\overrightarrow{u_{b_i}}\right]\overrightarrow{u_{b_i}}}{‖\overrightarrow{v_B}-\left[\overrightarrow{v_B}·\overrightarrow{u_{b_i}}\right]\overrightarrow{u_{b_i}}‖}·\overrightarrow{x_B}\right]+{\delta }_i,$$ E12

where $\overrightarrow{u_{b_i}}$ depends again on fin orientation, concretely in body axes, $\overrightarrow{u_{b_2}}=\left[010\right]$ , $\overrightarrow{u_{b_4}}=\left[0-10\right]$ , $\overrightarrow{u_{b_3}}=\left[001\right]$ , $\overrightarrow{u_{b_1}}=\left[00-1\right]$ .

Likewise, rocket moments in body axes include contributions from overturning, pitch damping, Magnus, and spin damping moments, which are described by the following:

$$\overrightarrow{O_B}=\frac{\pi }{8}{d}^3\rho \left(C_{M_{\alpha }}+C_{M_{{\alpha }^3}}{\alpha }^2\right){‖\overrightarrow{v_B}‖}^2\left(\overrightarrow{v_B}\times \overrightarrow{x_B}\right)$$ E13

$$\overrightarrow{P_{M_B}}=\frac{\pi }{8}\frac{d^3\rho }{I_y}{C}_{M_q}‖\overrightarrow{v_B}‖\left(\overrightarrow{L_B}-\left(\overrightarrow{L_B}·\overrightarrow{x_B}\right)\overrightarrow{x_B}\right)$$ E14

$$\overrightarrow{M_{M_B}}=-\frac{\pi }{8}\frac{d^4\rho }{I_x}{C}_{\mathit{mm}}\left(\left(\overrightarrow{L_B}·\overrightarrow{x_B}\right)\left(\left(\overrightarrow{v_B}·\overrightarrow{x_B}\right)\overrightarrow{x_B}\right)-\overrightarrow{v_B}\right)$$ E15

$$\overrightarrow{S_B}=\frac{\pi }{8}\frac{d^4\rho }{I_x}{C}_{\mathit{\text{spin}}}‖\overrightarrow{v_B}‖\left(\overrightarrow{L_B}·\overrightarrow{x_B}\right)\overrightarrow{x_B}$$ E16

where ${\mathrm{C}}_{{\mathrm{M}}_{\mathrm{\alpha }}}$ is the overturning moment linear coefficient, ${\mathrm{C}}_{{\mathrm{M}}_{{\mathrm{\alpha }}^3}}$ is the overturning moment cubic coefficient, ${\mathrm{C}}_{{\mathrm{M}}_{\mathrm{q}}}$ is the pitch damping moment coefficient, ${\mathrm{C}}_{\mathrm{mm}}$ is the Magnus moment coefficient and ${\mathrm{C}}_{\text{spin}}$ is the spin damping moment coefficient.

The control moment provided by the control surfaces may be expressed as follows:

$$\overrightarrow{C{M}_B}=\sum _{i=1}^4\left[\left[d_{\mathit{ax}}\overrightarrow{x_B}+d_{\mathit{lat}}\overrightarrow{u_{b_i}}\right]\times \overrightarrow{C{F}_{B_i}}\right],$$ E17

where $d_{\mathit{ax}}$ is the longitudinal distance, parallel to $\overrightarrow{x_B}$ , of fin centre of pressure (CP) to rocket centre of mass (CG), which depends on Mach number; $d_{\mathit{lat}}$ is the lateral distance, which is orthogonal to $\overrightarrow{x_B}$ and parallel to $\overrightarrow{u_{b_i}}$ for each fin, from fin centre of pressure to rocket centre of mass, which is supposed to be constant in this model.

To solve the motion of the rocket a body reference frame, which is coupled to the fuse, is used. Note that, because the fuse is uncoupled from the back part, which turns at high rates, Magnus force and moment and gyroscopic effects coming from the rear part must be modeled and included in the equations of motion. The turn rate of the back piece of the rocket is modeled as follows:

$$p_r=-\int \left({\mathrm{K}}_{\mathrm{s}}\delta \left(t_0\right)-\frac{\frac{\pi }{8}{d}^4\rho }{I_x}{C}_{\mathit{\text{spin}}}‖\overrightarrow{v_B}‖‖\left(\overrightarrow{L_B}·\overrightarrow{x_B}\right)\overrightarrow{x_B}‖\right)\mathit{dt},$$ E18

where $\delta \left(t_0\right)$ is a Dirac’s delta and K_s an experimental constant. Note that initial spin speed is modeled as an impulse which correlates to experimental data. It is accepted that the fuse mass is unimportant, which infers that non-apparent responses are included amongst fuse and aft part. Then, considering the aft impact is communicated as additional forces and moments to the Newton-Euler equations expressed in the B reference framework, the equations of motion may be expressed as follows: $\overrightarrow{F_B}=\frac{d\overrightarrow{{\mathit{mv}}_B}}{\mathit{dt}}+\overrightarrow{{\omega }_B^{\prime }}\times \overrightarrow{m{v}_B}$ and $\overrightarrow{M_B}=\frac{d\overrightarrow{L_B^{\prime }}}{\mathit{dt}}+\overrightarrow{{\omega }_B^{\prime }}\times \overrightarrow{L_B^{\prime }}$ , where $\overrightarrow{{\omega }_B^{\prime }}=\left(\begin{array}{c}{p}_r+p\\ q\\ r\end{array}\right)$ and $\overrightarrow{L_B^{\prime }}=\left[\begin{array}{ccc}{I}_x& 0& 0\\ 0& I_y& 0\\ 0& 0& I_y\end{array}\right]\overrightarrow{{\omega }_B^{\prime }}=\overline{\overline{I}}·\overrightarrow{{\omega }_B^{\prime }}$ are the angular speed and momentum, respectively, of the joint body, namely rocket and fuse. The aerodynamic and gyroscopic contributions of the aft part are computed separately and moved to the left part of Newton-Euler equations as follows:

$${\overrightarrow{M}}_r=-\frac{\pi }{8}{d}^3\rho \frac{C_{\mathit{mf}}}{I_x}\left[\overrightarrow{{\omega }_B^{\prime }}\overline{\overline{I}}·\overrightarrow{x_b}\right]\left(\overrightarrow{x_b}\times \overrightarrow{v_b}\right),$$ E19

$${\overrightarrow{M}}_{\mathit{Mr}}=-\frac{\frac{\pi }{8}{d}^4\rho }{I_x}{C}_{\mathit{mm}}\left(\left[\overrightarrow{{\omega }_B^{\prime }}\overline{\overline{I}}·\overrightarrow{x_b}\right]\left(\left(\overrightarrow{v_B}·\overrightarrow{x_B}\right)\overrightarrow{x_B}\right)-\overrightarrow{v_B}\right),$$ E20

$$\overrightarrow{G_r}=-\overline{\overline{I}}\frac{d}{\mathit{dt}}\overrightarrow{{\omega }_B^{\prime }}+\overrightarrow{{\omega }_B}\times \overrightarrow{L_B^{\prime }},$$ E21

$$\overrightarrow{F_{\mathit{ext}}}+{\overrightarrow{M}}_r=\frac{d\overrightarrow{{\mathit{mv}}_B}}{\mathit{dt}}+\overrightarrow{{\omega }_B}\times \overrightarrow{m{v}_B,}$$ E22

$$\overrightarrow{M_{\mathit{ext}}}+\overrightarrow{M_{\mathit{Mr}}}+\overrightarrow{G_r}=\frac{d\overrightarrow{L_B}}{\mathit{dt}}+\overrightarrow{{\omega }_B}\times \overrightarrow{L_B},$$ E23

where ${\overrightarrow{\mathrm{M}}}_{\mathrm{r}}$ is the Magnus force of the rotating part of the rocket, $\overrightarrow{{\mathrm{M}}_{\mathrm{Mr}}}$ is the Magnus moment of the rotating part of the rocket, $\overrightarrow{{\mathrm{G}}_{\mathrm{r}}}$ is the gyroscopic moment of the rotating part of the rocket, $\mathrm{p},\mathrm{q}\text{and}\mathrm{r}$ are the angular speed components of the fuse, $\overrightarrow{{\omega }_B}$ .

The conditions of movement given by Eq. (22) and Eq. (23) are integrated forward in time employing a fixed time step Runge-Kutta of fourth order to acquire a single flight trajectory.

3. Semi-active laser quadrant photodetector model

Semi-active laser kit consists of a quadrant photo detector that may be modeled as it is shown in Figure 4, where the external circle models the locator and the inward one the laser spot.

Figure 4.

Quadrant photo-detector configuration used.

In order to estimate laser footprint spot centre coordinates, electric intensities given by each of the photo-diodes (I1; I2; I3 and I4), which depend on area lit up by the laser spot, might be utilized. The following conditions characterize the most reasonable calculation, where $x_{\mathit{\text{quad}}},y_{\mathit{\text{quad}}}$ are the calculated laser footprint spot centre coordinates.

$$\left[\begin{array}{c}{x}_{\mathit{\text{quad}}}\\ y_{\mathit{\text{quad}}}\end{array}\right]=\left[\begin{array}{c}\mathit{ln}\frac{I_2}{I_4}\\ \mathit{ln}\frac{I_1}{I_3}\end{array}\right]$$ E24

The following mathematical relationship is always kept: $\frac{y_c}{x_c}=\frac{y_{\mathit{\text{quad}}}}{x_{\mathit{\text{quad}}}}$ , i.e., the transformation is conformal as showed in (25), where $x_c$ and $y_c$ are the genuine spot focus positions, not positions gotten by (24). Radial measurements may be interpolated introducing desired radius using the equivalences showed by the accompanying Table 2, where genuine and measured radial distances, $r_c$ and $r_{\mathit{\text{quad}}}$ , respectively, are given by Eqs. (26) and (27):

$${\theta }_c={\theta }_{\mathit{\text{quad}}}=\text{atan}\frac{y_{\mathit{\text{quad}}}}{x_{\mathit{\text{quad}}}}$$ E25

$$r_{\mathit{\text{quad}}}=\sqrt{x_{\mathit{\text{quad}}}^2+y_{\mathit{\text{quad}}}^2}$$ E26

$$r_c=f\left(r_{\mathit{\text{quad}}}\right)$$ E27

$r_{\mathit{\text{quad}}}$	0.48	0.99	1.50	2.01	2.67	3.68	5.88
$r_c$	0.1	0.2	0.3	0.4	0.5	0.6	0.7

Table 2.

Interpolation between measured radial distance, $r_{\mathit{\text{quad}}}$ and real radial distance, $r_c$ .

Then, the measurement output of the quadrant detector sensor may be expressed as it is indicated in (28), where $R_{\mathit{\text{quad}}}$ is the physical radius of the quadrant detector:

$$\left[\begin{array}{c}{x}_c\\ y_c\end{array}\right]=R_{\mathit{\text{quad}}}\left[\begin{array}{c}{r}_c\mathit{cos}{\theta }_c\\ r_c\mathit{sin}{\theta }_c\end{array}\right]$$ E28

4. Trajectory control algorithms

This section describes in detail the proposed navigation, guidance, and control algorithms.

Navigation for this vehicle alludes to the assurance, amid the totality of flight, of the rocket position and attitude, and target position. So as to decide rocket position integration of accelerometers’ data and GNSS sensor measurements might be hybridized and utilized as a part of the routing procedure. The elements of these estimations are not objective of this section.

As expressed before, one way for calculating attitude might be the estimation of various vectors in two reference frameworks.

On the off chance that a GNSS sensor gadget is prepared on the flying machine, velocity vector can be specifically calculated from sensor data in the NED axis. Another estimation of the velocity vector in body triad can be acquired from an arrangement of accelerometers prepared on the ship, one on each of the axis. These gadgets can quantify variations in speed. After integrating along time, from an initial condition given, velocity vector can be obtained.

A vector which might be utilized to characterize the rotation of the vehicle is the gravity vector. It is extremely easy to be resolved in NED ternion as it is constantly parallel to $z_{\mathit{NED}}$ . Note that accuracy might be expanded utilizing more entangled models, i.e., it can be demonstrated relying upon latitude and longitude.

The cornerstone of the displayed attitude calculation technique is estimating gravity vector in body axis. For instance, by deciding the constant component of the measured acceleration employing a low pass filter, where Jerk in body axes is calculated by derivation of acceleration; at that point, it is integrated so as to get the non-steady part of increasing speed, and, by subtracting this non-steady segment from the measured acceleration, gravity vector is evaluated. Be that as it may, this technique is not legitimate when the air ship pivots. Another technique to acquire gravity vector is integrating the mechanization equations; at that point, control thusly the subsequent conditions. Once more, gyros are required to implement this method. Basically, it consists on subtracting the contribution of the specific aerodynamic and inertial forces, which can be obtained from expressions (3) to (10) and dividing them by the vehicle mass, from the acceleration measurements of the accelerometers. The Eq. (29) aims to express this fact:

$$\overrightarrow{g_B}=\overrightarrow{A_B}-\frac{1}{m}\left(\overrightarrow{D_B}+\overrightarrow{L_B}+\overrightarrow{M_B}+\overrightarrow{P_B}+\overrightarrow{T_B}+\overrightarrow{C_B}+\overrightarrow{C{F}_B}\right),$$ E29

where $\overrightarrow{g_B}$ is the gravity vector expressed in body axes and $\overrightarrow{A_B}$ is the acceleration measured by accelerometers.

Another vector that may be expressed in both reference systems is the line of sight vector, which is especially useful during terminal phase, is the line of sight vector. It can be expressed in NED axes by subtracting target position from vehicle position obtained by GNSS measurements, and in body reference frame from SAL measurements as it is expressed in (30), where $d_p$ is the distance from the quadrant detector to the centre of mass of the rocket.

$$\overrightarrow{\mathit{LO}{S}_B}=\frac{\left[d_p{x}_c{y}_c\right]}{\sqrt{d_p^2+x_c^2+y_c^2}}$$ E30

Attitude can be determined operating with these pairs of vectors with matrixial algorithms as it is stated on [21].

The guidance consists of a modified proportional law, governed by the following equations: Eq. (31) gives the yaw error; Eq. (32) determines the pitch error; and Eq. (33) estimates time to impact ( $t_{\mathit{go}}$ ).

$${\psi }_{\mathit{err}}=\frac{\overrightarrow{\mathit{LO}{S}_{\mathit{NED}}}-\overrightarrow{v_{\mathit{NED}}}{t}_{\mathit{go}}}{t_{\mathit{go}}^2}·\left[001\right]$$ E31

$${\theta }_{\mathit{err}}=\frac{\overrightarrow{\mathit{LO}{S}_{\mathit{NED}}}-\overrightarrow{v_{\mathit{NED}}}{t}_{\mathit{go}}}{t_{\mathit{go}}^2}·\left[-100\right]$$ E32

$$t_{\mathit{go}}=\frac{1}{g}\overrightarrow{v_{\mathit{NED}}}·\left[010\right]+\frac{1}{g}\sqrt{{\left(\overrightarrow{v_{\mathit{NED}}}·\left[010\right]\right)}^2+2g\overrightarrow{\mathit{LO}{S}_{\mathit{NED}}}·\left[010\right]}$$ E33

Next, the utilized control law is presented, which gets as result two control parameters to be employed by the actuation system. Control is handled by a double loop feedback system, which uses accelerations and angular speed in body axes. The inner loop is only used as a system of stability augmentation. The two control parameters are the control angle for the rotating force ( ${\varphi }_c$ ) and the module of the control force ( ${\tau }_c)$ . The control angle for the rotating force is defined in Eq. (34), taking pitch ( ${\theta }_{\mathit{err}}$ ) and yaw ( ${\psi }_{\mathit{err}}$ ) errors as inputs. The module of the control force produced is also controlled. It is calculated in Eq. (35); note that this is done by processing the quadratic average of pitch and yaw errors. In these expressions $L1$ and $L2$ are experimental gains, $K_i$ , $K_d$ and $K_p$ are the integral, derivative, and proportional constants of the controller, and $K_{\mathit{mod}}$ is a constant to adjust the control force module.

$$\begin{array}{l}{\varphi }_c=K_p\text{[}atan\frac{{\theta }_{err}-L1\theta }{L2\left({\psi }_{err}-L1\psi \right)}-atan\frac{ac{c}_{zb}}{ac{c}_{yb}}\text{]}\hfill \\ +K_i\int \text{[}atan\frac{{\theta }_{err}-L1\theta }{L2\left({\psi }_{err}-L1\psi \right)}-atan\frac{ac{c}_{zb}}{ac{c}_{yb}}\text{]}dt\hfill \\ +K_d\frac{d}{dt}\text{[}atan\frac{{\theta }_{err}-L1\theta }{L2\left({\psi }_{err}-L1\psi \right)}-atan\frac{ac{c}_{zb}}{ac{c}_{yb}}\text{]}+atan\frac{{\theta }_{err}-L1\theta }{L2\left({\psi }_{err}-L1\psi \right)}\hfill \end{array}$$ E34

$${\tau }_c=K_{\mathit{mod}}\sqrt{{\left({\theta }_{\mathit{err}}-L1\theta \right)}^2+{\left(L2\left({\psi }_{\mathit{err}}-L1\psi \right)\right)}^2}$$ E35

Figure 5 shows the logic of the controller. It has three fundamental sources of info: the acceleration of the rocket in body axis, communicated by its three components [ $\mathit{ac}{c}_{\mathit{xb}}$ ; $\mathit{ac}{c}_{\mathit{yb}}$ ; $\mathit{ac}{c}_{\mathit{zb}}$ ], the pitch and yaw errors and the measurements from gyros on each axis, [ $\varphi ,\theta ,\psi$ ].

Figure 5.

Control system scheme.

Generally, the controller ascertains the required pointing angle of the aerodynamic force figuring the arc-tangent of the quotient of the pitch and yaw error. This gives an angle at which the aerodynamic force, in the $y_b-z_b$ plane, must point to reach the target. However, the gyroscopic effect due to the spinning part of the rocket makes the response difficult to govern, i.e., imposing a ${\varphi }_c$ of 90° will not make the rocket to respond upwards. Subsequently, the acceleration of the rocket must be likewise measured, without representing gravity, so as to have the effect between the difference between the angle that forms the projection of the aerodynamic force in the $y_b-z_b$ plane with $y_b$ and ${\varphi }_c$ .

In order to translate these control parameters into fin deflections, i.e., ${\delta }_1,{\delta }_2,{\delta }_3\mathit{\text{and}}{\delta }_4$ managed by two actuators, the relationships in (36) are applied.

$${\delta }_1={\delta }_3={\tau }_c\mathit{sin}{\varphi }_c;{\delta }_2={\delta }_4={\tau }_c\mathit{cos}{\varphi }_c$$ E36

5. Simulation results

MATLAB/Simulink R2016a on a desktop computer with a processor of 2.8 Ghz and 8 GB RAM was used. The rest of this section is divided in three different subsections. The first one presents the ballistic flights of the nominal trajectories to which the navigation, guidance and control algorithms developed will be applied. The second one describes the Monte Carlo simulations to be performed. And the last one compares the results of ballistic flights, controlled flights with GNSS/Accelerometer guided trajectory and controlled flights with GNSS/Accelerometer/Photo-detector guidance.

5.1. Ballistic trajectories

To test the developed algorithms, three nominal trajectories will be employed, which differ in their launch or initial pitch angle: 20°, 30° and 45°. Table 3 shows the characteristic parameters for these shots: initial pitch angle in the first column, initial lateral correction in the second one, and impact point in the last one. Initial lateral correction is performed in order to compensate Coriolis force and gyroscopic effects.

Initial pitch angle (^°)	Initial lateral correction (^°)	Impact point (m)
20	0.1524	18790.38
30	0.1989	23007.26
45	0.3082	26979.00

Table 3.

Nominal trajectories’ parameters.

The results for the ballistic trajectories for the three proposed initial pitch angles are shown in Figure 6. It shows impact point dispersion patterns for each of the ballistic cases. Also, the circular error probable (CEP) may be observed for each of the initial shot pitch angle.

Figure 6.

Ballistic shots for 20°, 30° and 45° initial pitch angles.

5.2. Monte Carlo simulations

Monte Carlo analysis is conducted to determine closed-loop performance across a full spectrum of uncertainty in initial conditions, sensor data acquisition, atmospheric conditions, and thrust properties. For atmospheric conditions variations in turbulence are considered using the specification MIL-F-8785C and the Dryden Wind turbulence model. Monte Carlo simulation distribution parameters are listed in the next Table 4. A set of 2000 shots is performed for each of the following combinations: ballistic shots, GNSS/Accelerometer assisted shots and GNSS/Accelerometer/Photo-Detector assisted shots. Initial shot angles of 20°, 30° and 45° are performed. Note that a total of 18,000 simulation shots are performed at the end of simulation campaign.

Parameter	Mean	Standard deviation
Initial φ	0°	20^°
Initial pitch	Nominal (20°, 30°, 45°)	0.01^°
Wind speed	10 m/s	5 m/s
Wind direction	0°	20°
Thrust at each time instant	T(t)	10 N
Initial azimuth deviation	Nominal lat. correction	0.01°

Table 4.

Monte Carlo simulation parameters.

5.3. Discussion

Values for navigation, guidance and control parameters defined on previous sections ( $C1$ , $C2$ , $K_i$ , $K_p$ , $K_d$ , $K_{\mathit{mod}}$ , $L1$ and $L2$ ) are expressed in Table 5. These parameters where selected experimentally in the model in order to obtain stable flight conditions.

Parameter	C1	C2	Ki	Kp	Kd	Kmod	L1	L2
Value	−7.5^°	−19 deg	1	0.3	0.025	0.08	0.01	1

Table 5.

Values for the constants.

Figure 7 shows detailed information about comparisons between different approaches. On the top, middle and bottom rows, shots with launch angles of 20, 30 and 45° are presented, respectively. Furthermore, on the left column ballistic flights and GNSS/Accelerometer assisted flights are compared, on the middle column GNSS/Accelerometer and GNSS/Accelerometer/Photo-detector assisted flights, and finally on the right column ballistic flights and GNSS/Accelerometer/Photo-detector assisted flights are compared for each of the three-initial pitch or launch angles. Controlled flights exhibit tighter impact groupings, getting tighter for the GNSS/Accelerometer/Photo-detector controller. Spread in the impact distribution does remain in the guided flights with GNSS/Accelerometer controller due to the difficulties discussed before, especially on sensors subsection, where it is explained the typical error of GNSS sensors and its associated accuracy problems during terminal guidance phase.

Figure 7.

Detailed shots for different algorithms.

The circular error probable (CEP) for each of the targets and for ballistic and controlled flights is shown in Table 6. The first column shows the initial pitch angle, the second one the CEP for the ballistic flight, the third one the CEP for the GNSS/Accelerometer Controlled Flight, and the last column the CEP for the GNSS/Accelerometer/Photo-detector Controlled Flight. The CEP for ballistic shots increases as initial pitch angle increases, while for controlled flights it remains stable, obtaining much better results for GNSS/Accelerometer/Photo-detector controller. Note that improvements or reductions on the CEP are above the 95%.

Initial pitch angle (^°)	Ballistic flight (m)	GNSS/accelerometer controlled flight (m)	GNSS/accelerometer/photo-detector controlled flight (m)
20	169.34	78.27	1.75
30	239.37	73.80	1.39
45	281.59	78.84	0.97

Table 6.

Circle error probable for the different cases.

6. Conclusions

A novel approach for navigation, guidance and control of high-rate spinning ballistic rockets, which is based on an innovative hybridization between GNSS/Accelerometer and semi-active laser quadrant photo-detector, has been developed.

Because ballistic rockets cannot be reused and high precision is of key importance, expensive equipment, e.g., sensors, has been traditionally employed. In this paper, however, it has been demonstrated that high precision may be obtained while using commercial off the shelf equipment, which is not usually highly precise.

Attitude determination is based on an algorithm which hybridizes data coming from multiple sensors and on a gravity vector estimator, avoiding the use of gyros. This approach is embedded in a two-phase guidance algorithm and a novel control technique for high-rate spinning rockets. The guidance algorithm is based on a modified proportional law while the control algorithm is based on a simple but effective and robust double-input double-output controller.

The proposed algorithms improve enormously accuracy by mixing those inaccurate signals in the terminal trajectory, with the signals of a precise semi-active laser quadrant detector, which is able to determine line of sight with high fidelity in body axes. Using the proposed hybridized algorithm during the last phases of flight, improves accuracy nearly to the ideal case as it was proved in simulations.