*Article* **Applying Neural Networks in Aerial Vehicle Guidance to Simplify Navigation Systems**

#### **Raúl de Celis \*, Pablo Solano and Luis Cadarso**

Aerospace Systems and Transport Research Group, European Institute for Aviation Training and Accreditation (EIATA), Rey Juan Carlos University, Fuenlabrada, 28943 Madrid, Spain; pablo.solano@urjc.es (P.S.); luis.cadarso@urjc.es (L.C.) **\*** Correspondence: raul.decelis@urjc.es; Tel.: +34-914888775

Received: 12 November 2020; Accepted: 10 December 2020; Published: 11 December 2020 -

**Abstract:** The Guidance, Navigation and Control (GNC) of air and space vehicles has been one of the spearheads of research in the aerospace field in recent times. Using Global Navigation Satellite Systems (GNSS) and inertial navigation systems, accuracy may be detached from range. However, these sensor-based GNC systems may cause significant errors in determining attitude and position. These effects can be ameliorated using additional sensors, independent of cumulative errors. The quadrant photodetector semiactive laser is a good candidate for such a purpose. However, GNC systems' development and construction costs are high. Reducing costs, while maintaining safety and accuracy standards, is key for development in aerospace engineering. Advanced algorithms for getting such standards while eliminating sensors are cornerstone. The development and application of machine learning techniques to GNC poses an innovative path for reducing complexity and costs. Here, a new nonlinear hybridization algorithm, which is based on neural networks, to estimate the gravity vector is presented. Using a neural network means that once it is trained, the physical-mathematical foundations of flight are not relevant; it is the network that returns dynamics to be fed to the GNC algorithm. The gravity vector, which can be accurately predicted, is used to determine vehicle attitude without calling for gyroscopes. Nonlinear simulations based on real flight dynamics are used to train the neural networks. Then, the approach is tested and simulated together with a GNC system. Monte Carlo analysis is conducted to determine performance when uncertainty arises. Simulation results prove that the performance of the presented approach is robust and precise in a six-degree-of-freedom simulation environment.

**Keywords:** nonlinear-flight-mechanics; neural networks; guidance, navigation, and control; machine learning; model; matlab-simulink

#### **1. Introduction**

Global Navigation Satellite Systems (GNSS) signals are widely utilized for aerospace applications. However, reliability decreases as the requirement of the application for which it is designed increases. The main cause producing this effect is the reduced signal/noise relationship caused by the attenuation and loss of the GNSS signal. This basically means that independent sources of data for navigation are needed to ameliorate these negative effects and reduce interference. Inertial Navigation Systems (INS) are a good example of devices which are independent of external perturbations. Particularly, an inertial estimation unit (IMU) is an electronic gadget that measures and reports a body's specific force, angular rate, and orientation, utilizing a blend of accelerometers, gyroscopes, and sometimes magnetometers. IMUs are normally used in airplanes, including unmanned aeronautical vehicles (UAVs), and spacecraft. However, these systems also feature important lacks, such as frequent incorrect initialization, accelerometer and gyroscope imperfections, which are trigger for cumulative errors

and imperfections in implemented dynamics model. Despite of this fact, Inertial Navigation Systems, when hybridized with GNSS receivers to minimize drift, are excellent for GNC data acquisition [1,2].

However, precision and cost are counterposed objectives. Reducing costs, while maintaining safety and accuracy standards, is key for development in aerospace engineering. Advanced algorithms for getting such standards while cutting down costs are are cornerstone. For example, to maintain an acceptable precision level while reducing costs, less precise devices may substitute expensive systems as long as GNSS signal is reachable and persistent to update the inertial system. However, many scenarios feature high uncertainty and alternatives are needed. An option to satisfy accuracy needs and budget limitations is to merge data of a few low cost sensors, which makes possible increases in accuracy levels.

The advantages of coordinated combination of information have appeared in numerous air applications [3]. For example, information combination strategies for six degrees of freedom rockets are depicted in [4]. The main issues in using various sorts of INS augmented with GNSS updates have been considered by [5]. Notwithstanding INS/GNSS hybridization, a set of nonlinear observers are presented by [6]. Note that, in case there are various sensors available, they may be additional contributions to a filter, e.g., the Kalman filter [1,2].

As it is shown in [1,2], the need to develop new Guidance, Navigation and Control (GNC) frameworks has fostered research on stability and controllability of aerospace vehicles. A novel guidance law is presented in [7], where only observations of line-of-sight angle and its rate of change coming from a seeker are employed. Ref. [8,9] present GNC cooperative techniques based on the conventional Proportional Navigation (PN). In [10] a target-follower engagement is considered, in which the target is followed while it tries to prevent interception. An attitude control-framework device for a spinning sounding rocket, which depends on a proportional, integral, and derivative (PID) controller, is created in [11]. Proportional-derivative GNC laws for the terminal phases of flight are proposed in [12,13]. In [14], a limited time concurrent sliding-mode GNC law is introduced. An overall scheme concerning the guidance and autopilot modules for a class of spin-stabilized balance controlled devices is introduced in [15].

Yet, even in GNSS/IMU hybrid devices, there exist negative influences, such as irregular estimations, which might be predominant during terminal guidance. Other methods, which are based on image recognition using multispectral cameras and other sensors, may be used in navigation for aerospace applications [16]. However, they usually feature high costs. Hence, advancement on new algorithms which may easily fulfill the required precision levels and budget limitations is a foundation in research. For instance, there are recent advances which consist of incorporating IMU, GPS, and laser guidance capacity, offering high accuracy and all-weather capacity [17,18].

Laser guidance may be provided by means of Semi Active Laser Kits (SAL). These devices are applied in many designing areas, such as calculating rotational speed of objects and estimating dynamics of laser spots [19,20]. The bonus of these kits is their favorable position during the last periods of the guidance, when they can provide high precision for GNC systems.

Therefore, it can be stated that sensor hybridization techniques [16,21] for viable and robust estimations are a current need when autonomy, accuracy, and minimal cost are to be achieved. However, also note that as the number of sensors to employed increases, the cost of system also increases. In this sense, Machine Learning (ML) techniques come onto the scene. They offer multitudinous options and innovative solutions of particular interest for GNC applications, where their foray is still latter and shallow, yet with no doubt promising. The utilization of ML strategies for the estimation of parameters dependent on the dynamics of aerial vehicles presents the bit of leeway that once the algorithm is calibrated or trained, it is not important to know the physical-mathematical establishments that rule the flight mechanics. Given the input signals, ML algorithms may restore the data that can later be utilized within the GNC system, such that the subsequent solutions will fit the genuine output [22,23]. Taking benefit of these facts, a reduced set of sensors may be selected to work together

with ML algorithms, all while safety and accuracy standards are matched, and complexity and costs are decreased.

However, the application of these strategies to a wide set of scenarios, which may also include uncertain conditions, depends largely on the representativity and amount of input and output data employed for training ML algorithms. This fact implies that desired performance stability and convergence is to be restricted to the trained mission envelope. Other approaches, which could ensure convergence and stability parameters under the proposed uncertain conditions, might also be employed for this type of application. For instance, adaptive control that uses adaptation laws to online estimate unknown system parameter variations for various mission envelopes [24–26].

Altogether, the objective of this paper is to improve current guidance strategies applying a powerful hybridization approach, which also introduces ML to enable attitude determination with a reduced availability of sensors, namely GNSS, accelerometer and semiactive laser quadrant photo-detectors. In particular, neural networks (NN) are implemented to precisely estimate the gravity vector to be combined with velocity and line of sight vectors in order to determine the attitude or rotation of the vehicle without needing gyroscopes. Note that the mentioned vectors need to be obtained in two different reference frames because otherwise the attitude determination problem cannot be solved.

#### *Contributions*

The main contribution of this scientific research is the application of Machine Learning techniques, i.e., neural network (NN) algorithms, to hybridize GNSS, accelerometer and semiactive laser quadrant photo-detectors signals. The role of the neural networks is to predict the gravity vector to estimate the attitude of the vehicle. Consequently, the advantage of such a hybrid system over the traditional ones, which are usually based on GNSS and IMUs, is the capability of eliminating gyros, which may be expensive and too sensitive for high demanding maneuvers and not reliable at all during some stages of flight.

The presented approach relies on neural networks and training algorithms to predict the gravity vector in body fixed axes while the vehicle is flying. The three components of the acceleration of the vehicle in body fixed axes are the inputs for the NN. After that, the predicted gravity vector is processed together, by means of a hybridisation algorithm, with velocity and line of sight vector to determine body rotation or attitude.

To reproduce the flight dynamics of an aerial vehicle, a nonlinear mathematical model is proposed, which considers nonlinear aerodynamic forces and moments and that has been validated to build up realistic conditions for simulation experiments [1,2]. On top of that, a robust double-input double-output control algorithm is employed to manage coupling among the normal and lateral nonlinear dynamics.

Note that the presented approach depends on the amount of available data for training, which means stability and convergence may be restricted to the trained mission envelope. However, note that training has been performed for a wide variety of launching, flight, and destination point conditions to resemble realistic settings, i.e., for a comprehensive set of missions. Overall, the methodology results in good enough quality results, even including good response to uncertainty in several conditions and characteristics, i.e., showing good GNC performance. Therefore, the presented research poses a path for a generalized and systematic application of NN/Machine Learning in GNC systems.

The rest of this paper is organized as follows. In Section 2, the system modeling is described in detail. Section 3 describes the navigation, guidance and control algorithms. Section 4 exposes simulations results. Finally, discussion and conclusions are presented.

#### **2. Vehicle Modeling**

This section is dedicated to the vehicle description, the flight dynamics model, and sensor and actuation models.

#### *2.1. Definition of the Vehicle*

The proposed GNC approach is applied to an aerial vehicle which features a maneuvering system composed of four canard surfaces, which is roll-decoupled from the main body of the vehicle. The motivation for this aerodynamic configuration is deeply explained in [1]. Note a canard is a small winglike surface attached to an aircraft forward of the main wing to provide extra stability or control, usually replacing the tail. Here, canards are decoupled 2 by 2, to produce control force and its related torque (see [1] for more details on this).

Table 1 shows some characteristic parameters of the vehicle, including thrust typical parameter values, vehicle and fuel mass, inertia, and aerodynamic parameters. These parameters are obtained from fluid dynamics numerical simulations, experimental measurements, and wind tunnel verification (see [1] for clarifications). Note that, to keep continuity and derivability on aerodynamic coefficients and thrust, a cubic splines based interpolation method has been employed at intermediate points. According to the shown moments of inertia, the vehicle features planes of symmetry.


**Table 1.** Aerial vehicle parameters.

#### *2.2. Equations of Flight Mechanics*

To construct the equations of flight, three reference frames are defined to project forces and moments: NED axes, working axes and body axes. NED axes, which are ground axes, are depicted by sub index *NED*. They are defined by *xNED* pointing north, *zNED* orthogonal to *xNED* and pointing nadir, and *yNED* yielding a clockwise trihedron. Working axes are represented by sub index *w*. They are given by *xw* pointing to the destination point, *yw* orthogonal to *xw* and pointing apex, and *zw* forming a clockwise trihedron. *AZ*<sup>0</sup> is the initial azimuth, i.e., the azimuth between *xe* and *xw*. Body axes are depicted by sub index *b*. *xb* pointing forward and contained in the plane of symmetry of the vehicle, *zb* orthogonal to *xb* pointing down and contained in the plane of symmetry of the vehicle, and *yb* shaping a clockwise trihedron. The origin of body axes is located at the gravity center of the vehicle

and they are rigid coupled to the roll-decoupled control device. Figure 1 shows the previously defined axes systems.

**Figure 1.** The three employed reference frames.

Next, flight dynamics equations are described. Note that these equations are compliant with the standard convention in [27]. Because the vehicle is assumed to be rigid, classical mechanics theory is employed. The Newton–Euler equations describe the combined translational and rotational dynamics of a rigid body. These laws, which are given by six equations, relate the motion of the center of gravity of a rigid body with the sum of forces and moments acting on the rigid body.

$$
\begin{bmatrix}
\overrightarrow{F\_{ext}} \\
\overrightarrow{M\_{ext}}
\end{bmatrix} = \begin{bmatrix}
\overrightarrow{L}^{\flat} + \overrightarrow{D} + \overrightarrow{P} + \overrightarrow{M} + \overrightarrow{T} + \overrightarrow{W} + \overrightarrow{C} \\
\overrightarrow{P\_M} + \overrightarrow{O} + \overrightarrow{M\_M} + \overrightarrow{S}
\end{bmatrix},
\tag{1}
$$

Equation (1) shows total external forces and moments acting on the vehicle: −→*<sup>L</sup>* is the lift force, −→*<sup>D</sup>* is the drag force, −→*<sup>P</sup>* is the pitch damping force, −→*<sup>M</sup>* is the Magnus force, −→*<sup>T</sup>* is the thrust force, −→*<sup>W</sup>* is the weight force, −→*<sup>C</sup>* is the Coriolis force, −→*PM* the is pitch damping moment, −→*<sup>O</sup>* is the overturn moment, −−→*MM* is the Magnus moment, and −→*<sup>S</sup>* is the spin damping moment.

$$\begin{bmatrix} \overrightarrow{\overline{\boldsymbol{D}}}\\ \overrightarrow{\overline{\boldsymbol{D}}}\\ \overrightarrow{\overline{\boldsymbol{D}}}\\ \overrightarrow{\overline{\boldsymbol{M}}}\\ \overrightarrow{\overline{\boldsymbol{W}}}\\ \overrightarrow{\overline{\boldsymbol{W}}} \end{bmatrix} = \begin{bmatrix} \left(\mathsf{C}\_{L\_{u}}(\boldsymbol{M})\cdot\boldsymbol{\alpha}+\mathsf{C}\_{L\_{x^{2}}}(\boldsymbol{M})\boldsymbol{\alpha}^{2}\right)\left(\|\overrightarrow{\boldsymbol{\overline{w}}}\_{w}^{\top}\|^{2}\overrightarrow{\boldsymbol{x}}\_{w}^{\top}-\left(\overrightarrow{\boldsymbol{x}}\_{w}^{\top}\cdot\overrightarrow{\boldsymbol{\overline{w}}}\_{w}^{\top}\right)\overrightarrow{\boldsymbol{v}}\_{w}^{\top}\right) \\\ \left(\mathsf{C}\_{D\_{0}}(\boldsymbol{M})+\mathsf{C}\_{D\_{x^{2}}}(\boldsymbol{M})\boldsymbol{\alpha}^{2}\right)\|\overrightarrow{\boldsymbol{v}}\_{w}^{\top}\|\overrightarrow{\boldsymbol{v}}\_{w}^{\top}\| \\\ \left.-\mathsf{d}\frac{\mathsf{C}\_{N\_{0}}(\boldsymbol{M})\|\overrightarrow{\boldsymbol{w}}\_{w}^{\top}\|^{2}\left(\overrightarrow{\boldsymbol{L}}\_{w}^{\top}\times\overrightarrow{\boldsymbol{x}}\_{w}^{\top}\right)}{\mathsf{d}\frac{\mathsf{C}\_{M\_{x^{2}}}(\boldsymbol{M})}{\mathsf{d}\boldsymbol{x}}\left(\overrightarrow{\boldsymbol{L}}\_{w}^{\top}\times\overrightarrow{\boldsymbol{w}}\_{w}^{\top}\right)}\right)\right] \\\ \left.\begin{array}{c} \mathsf{C}\_{w}(\boldsymbol{M})\overset{\textstyle \$$

For the computational experiments in this paper, the external forces in working axes are shown in expression (2), where *CL<sup>α</sup>* is the lift force linear coefficient, *CLα*<sup>3</sup> is the lift force cubic coefficient, *α* is the total angle of attack, *CD*<sup>0</sup> is the drag force linear coefficient, *CDα*<sup>2</sup> is the drag force square coefficient, −→*Lw* is the vehicle angular momentum in working axes, *Ix and Iy* are the vehicle inertia

moments in body axes, *CNq* is the pitch damping force coefficient, *Cm f* is the Magnus force coefficient, −*x* →*<sup>w</sup>* is vehicle pointing vector in working axes, −*g* <sup>→</sup>*<sup>w</sup>* is the gravity vector in working axes, −→<sup>Ω</sup> is earth's angular speed vector, −*v* →*<sup>w</sup>* is vehicle velocity in working axes, *d* is the reference surface of the vehicle, *ρ* is the air density, and *m* is the mass of the vehicle. Note that, −*x* →*<sup>w</sup>* is the unitary vector of the *xw* axis, and −*g* →*<sup>w</sup>* is the gravity vector in working axes. Be aware they are nonlinear functions of the variables describing the movement of the vehicle, such as aerodynamic speed, total angle of attack, Mach number, and aerodynamic parameters.

$$\begin{bmatrix} \overrightarrow{P\_{M}} \\ \overrightarrow{O} \\ \overrightarrow{O} \\ \overrightarrow{M}\_{M} \\ \overrightarrow{S} \end{bmatrix} = \frac{\pi}{8} d^{3} \rho \begin{bmatrix} \frac{1}{L\_{y}} \mathbb{C}\_{M\_{\text{q}}}(M) \|\overrightarrow{v\_{w}}\| \| \left( \overrightarrow{L\_{w}} - \left( \overrightarrow{L\_{w}} \cdot \overrightarrow{\mathbf{x}\_{w}^{\text{s}}} \right) \overrightarrow{\mathbf{x}\_{w}^{\text{s}}} \right) \\ \left( \mathbb{C}\_{M\_{\text{k}}}(M) + \mathbb{C}\_{M\_{\text{k}}}(M) \mathbf{a}^{2} \right) \|\overrightarrow{v\_{w}}\| \|^{2} \left( \overrightarrow{v\_{w}} \times \overrightarrow{\mathbf{x}\_{w}^{\text{s}}} \right) \\ - \frac{d}{\mathbb{L}\_{x}} \mathbb{C}\_{\text{num}}(M) \left( \left( \overset{\leftharpoonup}{L\_{w}} \cdot \overrightarrow{\mathbf{x}\_{w}^{\text{s}}} \right) \left( \left( \overrightarrow{v\_{w}} \cdot \overrightarrow{\mathbf{x}\_{w}^{\text{s}}} \right) \overrightarrow{\mathbf{x}\_{w}^{\text{s}}} \right) - \overrightarrow{v\_{w}} \right) \\ \frac{d}{\mathbb{L}\_{x}} \mathbb{C}\_{\text{s} \text{pin}}(M) \|\overrightarrow{v\_{w}}\| \left( \overrightarrow{L\_{w}} \cdot \overrightarrow{\mathbf{x}\_{w}^{\text{s}}} \right) \overrightarrow{\mathbf{x}\_{w}} \end{bmatrix}, \tag{3}$$

Similarly, the equations in (3) show the mathematical expressions for the moments, including overturning, pitch damping, Magnus, spin damping, and variables and parameters. Here, *CMq* is the pitch damping moment coefficient, *CM<sup>α</sup>* is the overturning moment linear coefficient, *CMα*<sup>3</sup> is the overturning moment cubic coefficient, *Cmm* is the Magnus moment coefficient and *Cspin* is the spin damping moment coefficient.

Control forces (−→*CF*) and moments (−→*CM*) are obtained from the maneuvering system, which is composed of four canard surfaces. Therefore, the contribution of each of them is summed to obtain the total control forces and moments.

$$\begin{bmatrix} \overrightarrow{\text{CF}} \\\\ \overrightarrow{\text{C}\dot{\text{M}}} \end{bmatrix} = \sum\_{i=1}^{i=4} \begin{bmatrix} \frac{1}{8} d^2 \rho \pi \|\overrightarrow{\text{v}}\_b^\flat\|^2 (\mathcal{C}\_{Na\_w}(\mathcal{M}) \delta\_i) \overrightarrow{n\_{ci}} \\\\ \frac{1}{8} d^3 \rho \pi \|\overrightarrow{\text{v}}\_b^\flat\|^2 (\mathcal{C}\_{Ma\_w}(\mathcal{M}) \delta\_i) (\overrightarrow{\text{x}}\_b^\flat \times \overrightarrow{n\_{ci}}) \end{bmatrix} \tag{4}$$

The expressions in (4) show the mathematical functions for control forces and moments, where *CN<sup>α</sup><sup>w</sup>* and *CM<sup>α</sup><sup>w</sup>* are the force and moments coefficients of the canard surface respectively, <sup>−</sup>*<sup>n</sup>* →*ci* is the normal vector of each canard, and *δ<sup>i</sup>* is the deflection angle of the canard surface. Here, −→*vb* is vehicle velocity in body axes.

$$
\begin{bmatrix}
\overrightarrow{C}\overrightarrow{F} + \overrightarrow{F\_{ext}} \\
\overrightarrow{CM} + \overrightarrow{M\_{ext}}
\end{bmatrix} = \begin{bmatrix}
\frac{dm\overrightarrow{\upsilon\_b}}{dt} + \overrightarrow{\omega\_b} \times m\overrightarrow{\upsilon\_b} \\
\frac{d\overrightarrow{L\_b}}{dt} + \overrightarrow{\omega\_b} \times \overrightarrow{L\_b}
\end{bmatrix} \tag{5}
$$

As stated before, a Newton-Euler approach is used to formulate the equations of motion of the aerial vehicle. These equations are in (5). Note that the body-fixed coordinate system (denoted by frame *b*) and the flat-Earth coordinate system (denoted by frame *e*) are related by Euler yaw (*ψ*), pitch (*θ*), and roll (*φ*) angles.

In Equations (5), −→*vb* stands for vehicle speed expressed in body axes, −*ω* →*<sup>b</sup>* for angular speed of the vehicle in body axes, and −→*Lb* for the angular momentum also in body axes. Recall that control and external forces and control and external moments must be expressed in body axes also to be employed in Equations (5).

#### *2.3. Sensors Models*

As exposed in the introduction, this research aims at simplifying navigation systems. Here, this means reducing the need for complex and/or expensive sensors. A gyroscope is a device used for measuring orientation and angular velocity, and it is widely employed in navigation systems. However, their precision downgrades for high-dynamics aerial platforms meanwhile costs increase if performance is to be maintained. Therefore, the objective is to avoid them by fusing information from

GNSS sensors, accelerometers, and photo-detector signals to improve vehicle navigation performance in terms of accuracy. This section aims at describing the employed models for these sensors.

#### 2.3.1. Global Navigation Satellite System (GNSS)

The GNSS sensor is modeled as a random noise and a bias added to the model calculated position. Note that these systems have typical accuracy of 3 m; therefore, the random noise and bias parameters have been adjusted to satisfy that performance. Because it is not the objective of this paper to model such a sensor, the reader is referred to [2,28] for more details on this.

These kind of sensors provide good performance during intermediate phases of flight and are employed to determine the line of sight vector expressed in NED axes. Note that GNSS sensors also provide velocity vector in NED axes. This vector is also modeled as a random noise and a bias of 0.1 ms<sup>−</sup>1, which resembles real performance of these sensors.

#### 2.3.2. Accelerometers

An accelerometer is a device that measures acceleration, i.e., the rate of change of velocity of the vehicle in its own instantaneous reference frame. They are modeled as a random noise and bias of 0.001 ms−<sup>2</sup> which resembles real performance of these sensors. Because they provide the acceleration vector expressed in body axes, velocity vector expressed in body triad can be obtained after integration of each of its components along time.

In addition, the magnitudes obtained from accelerometers will be used in the gravity vector estimation approach. As it is explained in the following sections, the velocity vector module is required to estimate it.

#### 2.3.3. Semiactive Laser Kit

Laser guidance may be used to guide a vehicle to a target by means of a laser beam. With this method, a laser is kept pointed at the objective and the laser radiation bobs off the objective and is dispersed every which way. At the point when the vehicle is close enough for a portion of the reflected laser energy from the objective to arrive at it, a laser seeker detects which direction this energy is coming from and provides a signal to correct the trajectory towards the source. The device seeking the laser and providing the signal is a semiactive laser kit.

The signal provided by this sensor features the centroid of the laser footprint in the photo-detector of the kit, which is composed of four photodiodes that convert light into an electrical current. To estimate its coordinates, the produced electrical intensities in the photodiodes (*I*1, *I*2, *I*<sup>3</sup> and *I*4) are employed, which rely upon the illuminated area. These coordinates can be determined as [*ln <sup>I</sup>*<sup>4</sup> *I*2 , *ln <sup>I</sup>*<sup>1</sup> *I*3 ] [17], and from them, it is possible to obtain the measured radial distance, *rquad*. Notwithstanding, real coordinates differ from those calculated, although the transformation is conformal [17]. To obtain the real radial distance, *rc*, the following mathematical functional relationship may be employed: *rc* = *f*(*rquad*). Then, cubic splines based interpolation is applied to obtain a continuous relationship. Equations (6) are utilized to estimate *xc* and *yc* (see [17] fore more details on this), which are the real spot center coordinates:

$$
\begin{bmatrix} x\_c \\ y\_c \end{bmatrix} = R\_{quad} \cdot \frac{r\_c}{r\_{quad}} \begin{bmatrix} \ln \frac{I\_4}{I\_2} \\ \ln \frac{I\_1}{I\_3} \end{bmatrix} \tag{6}
$$

where the physical radius of the photo-detector of the kit is given by *Rquad*. Consequently, the line of sight vector projected in body axes may be calculated from *xc* and *yc* and also from the distance of the photo-detector of the kit to the center of mass of the vehicle.

Note that the signal of this sensor is only available during the terminal phase of the flight. However, it is during final stages of flight when errors of 3 m in positioning target and vehicle induces enormous errors. In this way, an accurate terminal guidance sensor, for example, a semiactive laser kit,

is suggested for these last flight stages. This semiactive laser sensor, in combination with GNSS and accelerometers, will provide an accurate determination of the line of sight, especially in the terminal phase. For that purpose, the signals of these sensors' must be hybridized.

Next, GNC algorithms are presented. At their kernel, neural networks are implemented to determine gravity vector in two reference frames in order to determine vehicle attitude. In addition, hybridization algorithms are applied to sensors' signals to improve precision.

#### **3. Guidance, Navigation and Control (GNC) Algorithm Definition**

This section details the proposed navigation, guidance and control algorithms. A scheme of the overall process is depicted in Figure 2. The navigation function determines the position and attitude of the vehicle by means of the information sensed by the sensors. The position is determined through the integration of the signals provided by the accelorometers and the hybridization of the signal from a GNSS device. The determination of attitude is the core of the research in this paper. From the information provided by the accelerometers and the GNSS, the neural network determines the gravity vector in two different axes systems. Then Euler angles are devised from a triad algorithm. The guidance function compares the information from the navigation function with a reference and computes a desired action to the control function. The control function processes this desired action and transforms it into a command to the actuators of the plant (i.e., the vehicle), which execute the action. The action taken is again measured by the sensors, which closes the loop.

**Figure 2.** Scheme of the GNC process.

#### *3.1. Navigation*

Navigation refers to the determination during the flight, of the position and attitude of the vehicle, and target position.

Determining the position of the vehicle may consist of integrating accelerometers' measurements to be hybridized with GNSS sensor information. The details of these calculations are out of the scope of this paper, see [1,2,28,29] for details on this approach.

As it was mentioned before, determining attitude involves knowing two or more different vectors in two different reference systems. The velocity vector and the line of sight vectors can be the two needed vectors. If a GNSS sensor device is equipped on the aircraft, velocity vector can be directly obtained from sensor measurements in the NED triad, which can be expressed as shown in (7), where *vxNED* , *vyNED* and *vzNED* are the components of this velocity in NED axes.

$$
\overrightarrow{\boldsymbol{\upsilon}\_{NED}} = [\boldsymbol{\upsilon}\_{\boldsymbol{\upsilon}\_{NED}}, \boldsymbol{\upsilon}\_{\boldsymbol{\upsilon}\_{NED}}, \boldsymbol{\upsilon}\_{\boldsymbol{\varepsilon}\_{NED}}]^T \tag{7}
$$

In parallel, the same velocity vector can also be calculated in body triad from a set of accelerometers, one on each of the axes. Integrating each of their measures along time, the velocity vector is obtained as shown in (8). Here, *axB* , *ayB* and *azB* are the components of the acceleration in body axes as measured by the accelerometers and −→*ω<sup>b</sup>* is the estimated angular speed expressed in body axes. Note that, at this point, −→*ω<sup>b</sup>* is unknown, and the algorithm for estimating it will be shown in the following sections.

$$
\overrightarrow{\boldsymbol{w}}\_{B}^{\flat} = \int \left\{ [\boldsymbol{a}\_{\boldsymbol{x}\_{B}}, \boldsymbol{a}\_{\boldsymbol{y}\_{B'}}, \boldsymbol{a}\_{\boldsymbol{z}\_{B}}]^{T} + \overrightarrow{\boldsymbol{\omega}}\_{b}^{\sharp} \times \overrightarrow{\boldsymbol{\varpi}}\_{B}^{\flat} \right\} dt \tag{8}
$$

Similarly, the line of sight vector must be obtained in NED and body reference systems, −−−−−→ *LOSNED* and −−−→ *LOSB*, respectively. −−−−−→ *LOSNED* can be easily obtained from GNSS sensor information. However, the semiactive laser kit is needed to derive −−−→ *LOSB*, and this sensor signal is not available until the vehicle is close enough to the target. This means another vector is needed to successfully estimate the attitude of the vehicle during all the phases of flight.

The gravity vector poses as a natural candidate for such a challenge. Notice that determining the gravity vector in NED triad is straightforward. It is always parallel to −−→ *zNED*. Its expression is shown in (9), where *g* is the gravity acceleration, which is a fixed constant in this model (9.81 m/s2). Note that precision may be increased using more sophisticated models, i.e., it can be made variable with longitude, latitude, and altitude.

$$\overrightarrow{\text{g}\_{N\to D}} = \lg[0,0,1]^T \tag{9}$$

However, the gravity vector expressed in another reference system, i.e, body axes, is also needed. However, although there are multiple available approaches to obtain it, none of them is simple and/or require additional sensors. For example, it can be estimated determining the constant component of the measured acceleration employing a low pass filter, where Jerk in body axes is calculated by derivation of the acceleration; then, it is integrated to obtain the nonconstant component of acceleration and, by subtracting this nonconstant component from the measured acceleration, gravity vector may be estimated. However, this method is not valid when the aircraft rotates. Another method to obtain the gravity vector is to integrate the mechanization equations [30] to obtain it from the resulting expressions. However, gyros are needed to implement this method. Therefore, the keystone of the presented attitude determination method is determining gravity vector in body axes.

An estimation method for the gravity vector, which is valid for nonrotating and rotating aircraft and which is only based on accelerometers, is presented in the following subsection.

#### 3.1.1. Neural Network Based Gravity Vector Estimation

Among the numerous applications that machine learning offers to exemplary and current GNC issues (see [23,31–34]), its potential to precisely estimate the gravity vector from sensor information is one of the main unexplored settings. The utilization of neural networks (NN) to understand the evolution of nonlinear equations has been demonstrated before [35], regardless of uncertainty. Scientifically, this infers NN will learn flight mechanics equations [36] and produce an equal outcome.

The gravity vector estimation method presented here depends on the aerial platform on which it will be employed. This means that the method must be adjusted for the aircraft of interest. Without loss of generality, the estimation method detailed in this section is particularized for a four canard controlled aerial vehicle. However, using the appropriate neural network training, it may be applied to other aircraft.

The estimated gravity vector may be expressed as shown in (10), where its components in body axes are displayed. The point is to recoup a high precision gravity vector by consolidating the estimations from the accelerometers and the potential offered by machine learning.

$$\overline{\widehat{\mathbf{g}}}\_B^t = \begin{bmatrix} \widehat{\mathbf{g}}\_{\overline{x}\_{B'}}^- \widehat{\mathbf{g}}\_{\overline{y}\_{B'}}^- \widehat{\mathbf{g}}\_{\overline{z}\_B}^- \end{bmatrix}^T \tag{10}$$

In order to prove the suitability of the proposed approach two different methods or strategies are proposed, as they can be visualized in Table 2:


**Method 1** *axB* −→*g<sup>B</sup> ayB azB* **Method 2** *axB* −→*g<sup>B</sup> ayB azB*

**Table 2.** Neural network schemes for the two different methods or strategies.

The choice of the number and shape of neurons as well as the amount of training and validation data selected for the presented two strategies is a result of the literature review (specially from [35]) and a performed hyperparametric study. This study provided two points of interest: around 100 neurons and 50 neurons. Increasing the number or neurons or layers translated only into an increase of computation time for a limited improvement in terms of error of approximation. A further and detailed hyperparametric study will be performed in future work to precisely determine the optimal working point but is not the objective of this research to formalize this statement. The preliminary results of this hyperparametric study suggest that there is a limit number of neurons in the intermediate layer (estimated at about 100 neurons), and over this limit, results do not get improved.

Then, neural networks are trained replicating the flight dynamics problem. Two examples of the available 3 · <sup>10</sup><sup>8</sup> rows of data, which are obtained from 12,000 simulations, are showed in Table 3.


**Table 3.** Neural network input and target values.

Regarding the training process, three different backpropagation algorithms are employed: Scaled Conjugate Gradient (SCG) [37], Bayesian regularization (BR) [38,39], and Levenberg-Marquardt backpropagation (LM) [40,41]. The choice of these algorithms is a result of literature study. The percentage of data employed in this training is 70%. As it is common practice, a representative amount of sensor data and its corresponding gravity vector are left aside for validation purposes. In this case, a 15% of the available data corresponds to validation data. Note that the total amount of methods and training algorithms provide six different combinations which are analyzed next.

The performance of each of the six approaches can be quantified by means of the Mean Squared Error (MSE) and the Regression (R) parameter values. The MSE is the average squared difference between outputs and targets. Lower value means lower error. Zero means no error. R values measure the correlation between outputs and targets. An R value of 1 means a close relationship and 0 a random relationship. Other kind of error indicators (such as Mean Average Error, MAE) may also be used to monitor and validate the training to avoid overfitting.

For each of the the training processes, a maximum number of 1000 iterations has been established. As it is common practice in the field, classified by epochs. For the LM and SCG algorithms, training automatically stops when generalization stops improving, as indicated by an increase in the MSE of the validation samples. In the case of the BR algorithm, training stops according to adaptive weight minimization (regularization). In both cases the MAE is controlled as usual to avoid overfitting.

In addition, the trained NN is tested with the independent data (15% of the collected data), and MSE and R values are also calculated to validate the presented strategies.

Table 4 summarizes the obtained results for the training, validation and tests. It shows the values for the MSE and the R parameters. In the first column, the "Set" of data is defined, i.e., train (70%), validation (15%), or test (15%) data. The second column displays the employed training algorithm. The third and forth columns present the MSE and R values for the methods 1 and 2 showed in Table 2.


**Table 4.** MSE and R values for neural network based gravity vector estimator.

Analyzing the results in Table 4, it may be concluded that the best results are obtained for the combination of Method 1 and LM algorithm, which yields a MSE value of 5.63 · <sup>10</sup>−<sup>5</sup> and a Regression value of 0.998. Additionally, the combinations of Method 2 and LM and BR algorithms also result in acceptable values for MSE and R values, they are of the same order of magnitude. Consequently, we may conclude that Methods 1 and 2 provide good results when the LM and BR training algorithms are used. Nevertheless, the SCG training algorithm is not appropriate for this application, as the best results for this algorithm are 2 orders of magnitude worse as compared to the the rest of the algorithms.

Next, the attitude determination algorithm is presented. It is based on the estimated gravity vector by the neural networks (NN). In addition, note that, because there is information regarding

two additional vectors during terminal flight, i.e., the speed vector and the line of sight vector, a hybridization approach is also presented to improve performance.

#### 3.1.2. Attitude Determination Algorithm

Attitude determination can be determined by solving a classical Wahba's problem [42]. An orthonormal base must be defined in both axes systems, B and NED. This orthonormal base is defined for intermediate phases of flight, when signal of the photo-detector is not available and for the terminal phase of flight, when it is available, by unitary vectors*i*,*j* and*k* expressed in both bases. For the intermediate phases, these unitary vectors are calculated using the speed vector and the gravity vector. Furthermore, for the terminal flight, the line of sight vector and the gravity vector are to be employed. For the intermediate phases, *f l*, the unitary vectors are defined by expressions (11) and (12).

$$\overrightarrow{i\_{NED\_{fl}}} = \frac{\overrightarrow{v\_{NED}}}{\|\overrightarrow{i\_{NED}}\|}, \quad \overrightarrow{j\_{NED\_{fl}}} = \frac{\overrightarrow{v\_{NED}} \times \overrightarrow{g\_{NED}}}{\|\overrightarrow{v\_{NED}} \times \overrightarrow{g\_{NED}}\|}, \quad \overrightarrow{k\_{NED\_{fl}}} = \frac{\overrightarrow{i\_{NED\_{fl}}} \times \overrightarrow{j\_{NED\_{fl}}}}{\|\overrightarrow{i\_{NED\_{fl}}} \times \overrightarrow{j\_{NED\_{fl}}}\|}\tag{11}$$

$$\overrightarrow{i\_{B\_{fl}}} = \begin{array}{c} \overrightarrow{v\_{B}} \\ \frac{\overrightarrow{\langle v\_{B} \rangle}}{\|\overrightarrow{v\_{B}}\|} \end{array}, \quad \overrightarrow{j\_{B\_{fl}}} = \begin{array}{c} \overrightarrow{v\_{B}} \times \overrightarrow{\xi\_{B}} \\ \frac{\langle \overrightarrow{v\_{B}} \times \overrightarrow{\xi\_{B}} \rangle}{\|\overrightarrow{v\_{B}} \times \overrightarrow{\xi\_{B}}\|} \end{array}, \quad \overrightarrow{k\_{B\_{fl}}} = \begin{array}{c} \overrightarrow{i\_{B\_{fl}} \times \overrightarrow{j\_{B\_{fl}}}} \\ \frac{\langle \overrightarrow{i\_{B\_{fl}} \times \overrightarrow{j\_{B\_{fl}}}} \rangle}{\|\overrightarrow{i\_{B\_{fl}} \times \overrightarrow{j\_{B\_{fl}}}}\|} \end{array} \tag{12}$$

During the terminal guidance phase, *t f* , when the photo-detector is receiving information, a new set of unitary vectors is obtained by Equations (13) and (14).

$$\overrightarrow{j\_{NED\_{tf}}} = \frac{\overrightarrow{IDS\_{NED}}}{\|\overrightarrow{IDS\_{NED}}\|} \quad \overrightarrow{j\_{NED\_{tf}}} = \frac{\overrightarrow{IDS\_{NED}} \times \overrightarrow{j\_{NED}}}{\|\overrightarrow{IDS\_{NED}} \times \overrightarrow{j\_{NED}}\|} \quad \overrightarrow{k\_{NED\_{tf}}} = \frac{\overrightarrow{i\_{NED\_{tf}}} \times \overrightarrow{j\_{NED\_{tf}}}}{\|\overrightarrow{i\_{NED\_{tf}}} \times \overrightarrow{j\_{NED\_{tf}}}\|} \tag{13}$$

$$\overrightarrow{i\_{B\_{tf}}} = \begin{array}{c} \overrightarrow{\overline{LOS}\_{B}}^{\flat} \\ \left| \overrightarrow{LOS}\_{B}^{\flat} \right| \end{array}, \quad \overrightarrow{j\_{B\_{tf}}} = \begin{array}{c} \overrightarrow{\overline{LOS}\_{B}^{\flat} \times \overrightarrow{g\_{B}}^{\flat}} \\ \left| \overrightarrow{LOS}\_{B}^{\flat} \times \overrightarrow{g\_{B}}^{\flat} \right| \end{array}, \quad \overrightarrow{k\_{B\_{tf}}} = \begin{array}{c} \overrightarrow{i\_{B\_{tf}}} \times \overrightarrow{j\_{B\_{tf}}} \\ \left| \overrightarrow{i\_{B\_{tf}}} \times \overrightarrow{j\_{B\_{tf}}} \right| \end{array} \tag{14}$$

Note that to determine the attitude of a vehicle with respect to a reference frame, the direction cosine matrix (DCM) must be determined. It represents the attitude of the body frame (B) relative to the reference frame (NED). It is specified by a 3 × 3 rotation matrix, where the columns represent unit vectors in the body axes projected along the reference axes. Therefore, the expression to determine the DCM is as shown in Equation (15), where −→*iBi* , −→*jBi* , −→*kBi* is a 3 × 3 square matrix composed of orthonormal vectors in body triad, −−→ *iNEDi* , −−−→ *jNEDi* , −−−→ *kNEDi* expresses the same concept in NED triad, and *DCMB*,*NEDi* is the director cosine matrix that transforms NED triad into body triad. Notice that depending on the phase of flight, i.e., intermediate (*f l*) or terminal (*t f*), the matrix may be calculated with different inputs.

$$\left[ \overleftarrow{i\_{B\_{i'}}} \overrightarrow{j\_{B\_{i'}}} \overrightarrow{k\_{B\_{i}}} \right] = D \text{CM}\_{\text{B\_{i}NED\_{i}}} \left[ \overrightarrow{i\_{NED\_{i'}}} \overrightarrow{j\_{NED\_{i'}}} \overrightarrow{k\_{NED\_{i}}} \right] \forall i \in \{tf, fl\} \tag{15}$$

The DCM matrix can be solved from Equation (15) as it is shown in Equation (16). Employing an orthonormal base simplifies the calculation of the inverse matrix as it is the transposed matrix.

$$D\mathbb{C}M\_{\mathbb{B},\mathbb{N}\boldsymbol{ED}\_{i}} = \left[\overrightarrow{i\_{\mathbb{B}\_{i'}}}\overrightarrow{j\_{\mathbb{B}\_{i'}}}\overrightarrow{k\_{\mathbb{B}\_{i}}}\right] \left[\overrightarrow{i\_{\mathbb{N}\boldsymbol{ED}\_{i'}}}\overrightarrow{j\_{\mathbb{N}\boldsymbol{ED}\_{i'}}}\overrightarrow{k\_{\mathbb{N}\boldsymbol{ED}\_{i}}}\right]^{T} \forall i \in \{tf,fl\} \tag{16}$$

After obtaining the two different director cosine matrices, which will be essentially similar matrices, the rotation is characterized. The most suitable method to express this rotation is through quaternions, as they avoid any possible singularities on the poles of rotation. It is widely known that quaternions themselves are enough to express rotations without singularities, but it is also known that conceptually they are difficult to be visualized. An easier manner to define these rotations is by means of Euler angles. Concretely, the most common aeronautical rotation is defined by roll (*φi*), pitch (*θi*), and yaw (*ψi*) angles for *i* ∈ {*t f* , *f l*}. A method to obtain a quaternion solution is explained

in [2]. Note that two different values for each quaternion are obtained, *i* ∈ {*t f* , *f l*}. This fact requires an hybridization between them in order to only obtain one value for each quaternion.

#### Hybridization Algorithm

The Euler angles (or their corresponding quaternions) values obtained for *i* ∈ {*t f* , *f l*} are hybridized applying the recursive algorithm described in (17) and (18):

$$\left\{ \overrightarrow{\operatorname{Eul}} \right\}\Big|\_{n} = \begin{cases} \left. \left\{ \overrightarrow{\operatorname{Eul}} \right\} \right|\_{n-1} + \kappa|\_{n} \left[ \left. \left\{ \overrightarrow{\operatorname{Eul}} \right\} \right|\_{n} \right|\_{n} - \left. \left\{ \overrightarrow{\operatorname{Eul}} \right\} \right|\_{n-1} \right] \text{ if } \exists \left\{ \overrightarrow{\operatorname{Eul}} \right\} \Big|\_{n} \Big|\_{n} \\\left. \left\{ \overrightarrow{\operatorname{Eul}} \right\} \right|\_{n-1} + \kappa|\_{n} \left[ \left. \left\{ \overrightarrow{\operatorname{Eul}} \right\} \right|\_{n} - \left. \left\{ \overrightarrow{\operatorname{Eul}} \right\} \right|\_{n-1} \right] \text{ if } \exists \left\{ \overrightarrow{\operatorname{Eul}} \right\} \Big|\_{n} \end{cases} \tag{17}$$
 
$$\kappa|\_{n} = \Gamma \cdot [\Gamma + \Lambda]^{-1} \tag{18}$$

where −→*Eul* are the Euler angles (*φ*, *<sup>θ</sup>*, *<sup>ψ</sup>*), and <sup>Γ</sup> and <sup>Λ</sup> are the error covariance matrices for *<sup>i</sup>* = *t f* and for *<sup>i</sup>* = *f l* measurements, which are set to 1.3 · <sup>10</sup>−<sup>6</sup> and 0.95 · <sup>10</sup>−3, respectively. Those values were determined empirically.

The Euler angles obtained in (17) may be used to characterize rotations and angular speeds in navigation, guidance, and control algorithms. This basically means that −→*ω<sup>b</sup>* is now known. Furthermore, from these Euler angles, the hybridized director cosine matrix (*DCMB*,*NED*) may be calculated.

#### *3.2. Guidance Law*

Guidance is given in two stages. The first comprises of a constant angle glide trajectory, while the second one is based on a modified proportional law.

#### 3.2.1. Constant Angle Glide Trajectory

Equation (19) proposes a law which is chosen to increase range. It adjusts the longitudinal axis of the vehicle (*xb*) with a vertical flight plane, orthogonal to ground, parallel to the line joining the gravity center of the vehicle and the destination target and containing the gravity center of the vehicle. The line of sight is expressed in working axes is given by vector −−−→ *LOSw* = [*LOS*1*<sup>w</sup>* , *LOS*2*<sup>w</sup>* , *LOS*3*<sup>w</sup>* ]. *xb* in working axes is represented by vector −*x* →*bw* = [*xb*<sup>1</sup>*<sup>w</sup>* , *xb*<sup>2</sup>*<sup>w</sup>* , *xb*<sup>3</sup>*<sup>w</sup>* ]. Consequently, the lateral correction to be applied (*ψdem*) is calculated by the first component of Equation (19), while the correction in the vertical plane (*θdem*) with respect to a constant glide angle trajectory given by *C*<sup>1</sup> [1] is given by the second component. Note that guidance effectively starts after apogee, which is determined by the pitch angle (*θ*) and after fuel burn time.

$$\begin{aligned} \left[ \begin{array}{c} \Psi\_{dem} \\ \theta\_{dem} \end{array} \right] = \left\{ \begin{array}{c} \left( \arctan \frac{\text{LOS}\_{\text{W}}}{\text{LOS}\_{\text{I}\_{\text{W}}}} - \operatorname\*{atanx} \frac{\text{x}\_{\text{I}\_{\text{W}}}}{\text{x}\_{\text{I}\_{\text{I}\_{\text{W}}}}} \right) \\\\ \mathbf{C}\_{1} \end{array} \right] \text{ if } t > 5 \text{ and } \theta \le 0 \end{aligned} \tag{19}$$

#### 3.2.2. Modified Proportional Law

The guidance for the terminal phase of flight is formulated as a modified proportional law ruled by expression (20). Equation (21) calculates time to target, *tgo*. Guidance is activated only when the vertical coordinate of the line of sight vector is greater than a given constant (*C*2) [1].

$$
\begin{bmatrix}
\Psi\_{dem} \\
\theta\_{dem}
\end{bmatrix} = \begin{cases}
\begin{bmatrix}
\frac{\overrightarrow{LOS}\_w^\flat - \overrightarrow{v\_w}^\flat t\_{\mathcal{B}^\flat}}{t\_{\mathcal{B}^\flat}^2} \cdot \begin{bmatrix}
\frac{\overrightarrow{\Phi}\_{\mathcal{W}}^\flat}{\theta}
\end{bmatrix} \text{if } atan \frac{LOS\_{\text{NED}}}{\sqrt{LOS\_{\text{NED}}^2 + LOS\_{\text{NED}}^2}} \le \mathcal{C}\_2 \\
\\ \begin{bmatrix}
0 \\ 0
\end{bmatrix} \text{ else}
\end{cases}
\tag{20}
$$

$$t\_{\mathcal{S}^0} = \max\left[\begin{array}{c} \frac{1}{\mathcal{S}} \left( \overrightarrow{v\_w} \cdot \overrightarrow{j\_w} + \sqrt{(\overrightarrow{v\_w} \cdot \overrightarrow{j\_w})^2 + 2\text{g}LOS\_{2w}} \right) \\ \frac{1}{\mathcal{S}} \left( \overrightarrow{v\_w} \cdot \overrightarrow{j\_w} - \sqrt{(\overrightarrow{v\_w} \cdot \overrightarrow{j\_w})^2 + 2\text{g}LOS\_{2w}} \right) \end{array} \right] \tag{21}$$

Here, −→*iw* , −→*jw* , and −→*kw* represent the orthonormal basis of the working axes, and −−−→ *LOSw* <sup>=</sup> [*LOS*1*<sup>w</sup>* , *LOS*2*<sup>w</sup>* , *LOS*3*<sup>w</sup>* ] *<sup>T</sup>* and −−−−−→ *LOSNED* = [*LOS*1*NED* , *LOS*2*NED* , *LOS*3*NED* ] *<sup>T</sup>* are the vectors of the line of sight and its components in working and NED axes, respectively.

#### *3.3. Control System*

The control law presented in [29] is employed in the current research. Two control conditions are presented in the actuation device: the modulus and the angle for the control force. Control is computed by a double loop feedback system. The inner loop aims at augmenting the stability of the vehicle. Equation (22) characterizes modulus (*τc*) and angle (*φc*) of the control force. Its inputs are pitch (*θdem*) and yaw (*ψdem*) errors. *Ki*, *Kd* and *Kp* are the integral, derivative and proportional constants of the controller, *Kmod* is a constant to adjust force module and *L*1 and *L*2 are experimental gains. The procedure to decide these constant values, which appear in Table 5, is clarified in [1]. Note that the acceleration, without accounting for gravity, of the vehicle in body axes is defined by [*accxb*, *accyb*, *acczb*]. Euler angles as introduced before are represented by [*φ*, *θ*, *ψ*].

$$
\begin{bmatrix}
\Phi\_{\rm c} \\
\texttt{\tau}\_{\rm c}
\end{bmatrix} = \begin{bmatrix}
K\_{\rm p} \left(E1 - E2\right) + K\_{\rm i} \int \left(E1 - E2\right) dt + K\_{\rm d} \frac{d}{dt} \left(E1 - E2\right) + E1 \\
K\_{\rm md} \sqrt{(\theta\_{\rm dem} - L\_{1}\theta)^{2} + (L\_{2}(\Psi\_{\rm dem} - L\_{1}\Psi))^{2}}
\end{bmatrix}
$$

$$
\text{where}
\quad
\begin{cases}
E1 = \operatorname{atan} \frac{\theta\_{\rm dem} - L\_{1}\theta}{L\_{2}(\Psi\_{\rm dem} - L\_{1}\Psi)} \\
E2 = \operatorname{atan} \frac{\theta \mathcal{E}\_{\rm cik}}{\theta \mathcal{E}\_{\rm yb}}
\end{cases}
\tag{22}
$$

**Table 5.** Values for the constants of the control systems for each flight phase.


Summarizing, the control law works as follows. The controller determines the required pointing angle of the aerodynamic force. This is calculated obtaining the arc-tangent of the quotient of the pitch and yaw error, which provides an angle, in the *yb* − *zb* plane, at which the aerodynamic force should point to reach the objective. However, due to gyroscopic effects, the response of the vehicle is hard to govern, i.e., pointing the control force upwards may not make the vehicle to react upwards. Consequently, knowing the acceleration of the vehicle, without accounting for gravity, is a must. Similarly, the difference between *φ<sup>c</sup>* and the angle the projection of the aerodynamic force in the *yb* − *zb* plane forms with *yb* needs to be determined [1].

The aforementioned parameters of control are transformed into canard surface deflections, i.e., *δ*1, *δ*2, *δ*<sup>3</sup> and *δ*4, which are ruled by two different actuators, as it is shown in Equation (23).

$$\begin{bmatrix} \delta\_1\\ \delta\_2\\ \delta\_3\\ \delta\_4 \end{bmatrix} = \mathbf{r}\_c \begin{bmatrix} \sin \phi\_c\\ \cos \phi\_c\\ \sin \phi\_c\\ \cos \phi\_c \end{bmatrix} \tag{23}$$

#### **4. Numerical Simulations**

The described nonlinear dynamics are integrated forward in time utilizing a fixed time step Runge–Kutta method of fourth grade to get a single flight path. [1] shows the validation of this modeling and solving approach for aerial platforms. To demonstrate the precision of the novel methodology introduced in this research, which is based on neural networks, the obtained results are compared to the obtained outcomes in [29]. The methodology in [29] features a Kalman based hybridization [43,44] of GNSS, IMU and semiactive laser quadrant photo-detector sensors. To integrate the equations of motion and their interactions with GNC system and environment, MATLAB/Simulink R2020a on a personal computer with a processor of 2.8 Ghz and 32 GB RAM is used.

The remainder of this section is separated in two subsections. The first one presents the noncontrolled trajectories to which the developed navigation, guidance, and control algorithms will be applied. The second one performs Monte Carlo simulations of ballistic flights, Kalman hybridization based controlled flights, and neural networks based controlled flights. Moreover, an ideal controller without induced errors in the line of sight is also developed to compare results with ideal results.

#### *4.1. Noncontrolled Trajectories*

Three nominal trajectories are established to test the developed approach. Each one differs in its initial pitch angle: 20◦, 30◦ and 45◦. Destination points are at 18,790 m, 23,007 m, and 26,979 m, respectively. In order to compensate Coriolis and gyroscopic forces, initial lateral correction is set to 0.15◦, 0.19◦ and 0.31◦, respectively. Figure 3 shows many of these trajectories in a 3D setting for different settings.

**Figure 3.** Noncontrolled flights for 20◦ (black), 30◦ (red) and 45◦ (blue) initial pitch angles.

#### *4.2. Monte Carlo Simulations*

Noncontrolled flights have been validated against real data provided by the Spanish air force. Monte Carlo simulations are performed to calculate closed-loop performance over a full range of uncertainty settings, which have been defined with the support of the Spanish air force. These settings model the potential uncertainty that can arise in aspects such as initial conditions, sensor information procurement, weather conditions, and thrust properties. Note that, details on uncertainty models for sensors are given in the previous sections. Table 6 shows mean and standard deviation for the rest of the considered uncertain parameters.

A set of 2000 flights is performed for each of the following approaches: noncontrolled flights, Kalman hybridization based controlled flights, neural network based controlled flights, and ideally controlled flights. Simulations are run for each of the initial pitch angles. Note that the previously used data for neural network training is different from the data employed in the simulations in this section.


**Table 6.** Monte Carlo simulation parameters.

#### *4.3. Discussion*

Results for noncontrolled trajectories are shown in Figure 3, which depicts destination point dispersion patterns. This figure shows the ballistic trajectory that the vehicle follows for three different launch angles without using the control system at all, that is, it shows the flight of the system before implementing the improvements provided by the GNC system. The circular error probable (CEP), which is defined as the radius of a circle, centered on the mean, whose boundary is expected to include the landing points of 50% of the flights. The CEP is employed as a quality check parameter at the final step of the simulation as it is a valid reference for any utilized method. Indeed, the lower the CEP is the better the global GNC device is.

Table 7 displays the CEP for noncontrolled flights, Kalman hybridization based controlled flights and ideally controlled flights for each of the initial pitch angles.

**Table 7.** The CEP for noncontrolled flights, Kalman hybridization based controlled flights and ideally controlled flights for 20◦, 30◦ and 45◦ initial pitch angles.


Analyzing Table 7, it may be concluded that the CEP increases with target distance for noncontrolled flights, as expected. However, it almost remains constant for either Kalman based or ideally controlled flights, which means the use of an appropriate GNC device eliminates the correlation between the CEP and the distance to the objective. Recall, the main purpose is to show that these results are reproducible when employing machine learning and when reducing the availability of sensors.

Table 8 shows again the CEP for different trajectories, now as obtained by the novel presented approach. Each row, excluding the first one, which shows the headings of the columns, displays the resulting CEPs for every combination of trajectory and NN training algorithm. The first column in the table identifies the trajectory, the second one the training algorithm, the third one the CEP for the NN architecture in method 1, and the last one the CEP for the NN architecture in method 2. One of the main conclusions drawn from Table 8 is that the results for the SCG training algorithm present an unacceptable big CEP, which means low accuracy when reaching destination. Consequently, this training algorithm should be discarded in this case. Indeed, we recover results equivalent to a defective GNC system, even showing worse performance than ballistic flights. Several tests and a hyperparametric analysis were conducted, and it is observed that these kinds of errors were systematic, concluding that this training algorithm does not match well with the fed data to the NN. However, the results from BR and LM training algorithms show a good behavior throughout the trajectory, which means level of accuracy at destination is high. Diving into the numerical results in Table 8, it can be stated that the presented novel approach, which relaxes sensor requirements, is even able of outperforming the Kalman hybridization based approach. It can also be observed that the results are coherent with the training results depicted in Table 4. For example, the poor MSE and R values for the SCG training algorithm are reflected in the unacceptable GNC device results. Comparing the

obtained CEPs to what it was obtained in [1,22,29], it can be concluded that the results here are of the same order of magnitude and that the NN algorithms are viable for these type of applications.

As a summary, results for ballistic trajectories and comparisons between different approaches are shown in Figure 4. It is composed of four columns and three rows of subfigures. Each row features a different initial pitch angle. The first column of subfigures compares ballistic flights against Kalman hybridization assisted flights, the second one compares Kalman hybridization against neural network hybridization, the third one neural network hybridization against ideal controller, and the last one ballistic flights against neural network hybridization assisted flights. Note that, even with an ideal controller, which features perfect information on the attitude angles, there are still errors associated to the aerodynamic response of the vehicle.

**Figure 4.** Detailed shots for different algorithms.

**Table 8.** The CEP in neural network based controlled flights for the different methods and training algorithms for 20◦, 30◦ and 45◦ initial pitch angles.


#### **5. Conclusions**

A novel methodology, which depends on an innovative hybridization among several sensor signals, has been proposed. At the core of the approach, neural networks are employed to get estimations of the gravity vector, which allows avoiding the use of gyroscopes. Traditional GNSS/IMU frameworks feature little errors of up to one meter, which may imply huge mistakes in line of sight vector computation when separation to the objective is small. With the proposed approach the exactness of line of sight calculation can be improved during the terminal GNC, enhancing the accuracy at the destination point, while sensor needs are lowered.

The proposed approach employs information gathered from GNSS, acceloremeters, and a semiactive laser kit. With that information, two different neural network architectures are applied to estimate the gravity vector in order to determine the attitude of the vehicle. Three training algorithms have been addressed to tune the parameters in the neural networks. In total, six different strategies are developed for estimating the gravity vector along the trajectory. In addition, because the methodology allows for determining attitude in two ways, the information is hybridized with the aim of augmenting precision.

This innovative methodology is integrated into a two-phase guidance algorithm for aerial vehicles, which provides the required input data for the GNC system. The guidance law is founded on a constant glide angle and on a modified proportional law. The control algorithm is based on a robust and effective but simple double-input double-output device. Overall, the resulting GNC system presents excellent values regarding dispersion at the destination objective, significantly increasing the precision for noncontrolled flights, as expected, but also matching accuracy as provided by other GNC systems requiring more sensors on-board. Note that results also show good behavior of the system under uncertainty conditions. Summarizing, the developed approach, which is based on neural networks, shows that precision levels can be matched or improved as compared to other methodologies.

Future research will address the increase of the use of neural networks in other modules of GNC algorithms to further simplify the overall architecture.

**Author Contributions:** Conceptualization, R.d.C.; methodology, R.d.C. and L.C.; software, L.C.; validation, R.d.C.; formal analysis, R.d.C. and P.S.; investigation, R.d.C. and P.S.; resources, L.C.; data curation, R.d.C. and P.S.; writing–original draft preparation, L.C.; writing–review and editing, R.d.C.; visualization, R.d.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Project Grant F663—AAGNCS by the "Dirección General de Investigación e Innovación Tecnológica, Consejería de Ciencia, Universidades e Innovación, Comunidad de Madrid" and "Universidad Rey Juan Carlos".

**Acknowledgments:** The authors would like to thank Lieutenant Colonel Jesús Sánchez (NMT) of the National Institute for Aerospace Technology (INTA) for the solid modeling of the concept.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
