Aeroservoelastic Control Law Development for the Integrated Adaptive Wing Technology Maturation Wind-Tunnel Test

Aeroservoelastic Control Law Development for the Integrated

Adaptive Wing Technology Maturation Wind-Tunnel Test

∗ † ‡ § Josiah M. Waite, Jared Grauer, Robert E. Bartels, and Bret K. Stanford NASA Langley Research Center, Hampton, VA, 23681 I. Introduction he NASA Advanced Air Transport Technology (AATT) Project seeks to “explore and develop technologies and T concepts for improved energy efficiency and environmental compatibility for fixed wing subsonic transports” [ 1 ].

The IAWTM project is a subproject within AATT and is a joint Boeing/NASA effort to, in part, demonstrate the active control systems necessary to employ high-aspect ratio wings for improved aerodynamic efficiency. The aircraft configuration under study for the IAWTM project is a generic jet transport based on the NASA Common Research Model (CRM) with a wing aspect ratio of 13.5 (as opposed to nominal CRM aspect ratio of 9). The development of this modified CRM is given in Ref. [2].

A semi-span model of the high-aspect-ratio CRM with 10 trailing edge control surfaces is currently being fabricated and will be delivered to the NASA Langley Transonic Dynamics Tunnel (TDT). There, several tests will be conducted to demonstrate maneuver load alleviation (MLA), gust load alleviation (GLA), active flutter suppression (AFS), and drag reduction, all at transonic conditions. A layout of the wind tunnel model is shown in Fig. 1. The wind-tunnel Fig. 1 Layout of the wind-tunnel model showing the placement of all ten active control surfaces. Figure courtesy Boeing Research and Technology.

model is a spar/pod wing design with a relatively rigid fuselage attached to a mounting structure (seen in Fig. 1), which in turn is attached to a load balance behind the tunnel wall. Shims exist that can be attached to force the load to go through the fuselage rather than through the balance so as not to damage the balance in the event of flutter. The spar/pod design is similar to that used in the Boeing truss-braced wing aeroelastic model that was tested in the TDT in 2013 [ 3 ].

There is a central spar running the span of the wing to which the pods are attached that form the aerodynamic shape of the wing. The pods are constructed of ribs for attachment to the spar and skins to create the outer mold line. There are ten trailing-edge active control surfaces that will be used for the test objectives discussed above. Three of these (two outboard ailerons and one inboard aileron) are high-speed electrohydraulic actuators, and the remaining seven are low-speed mini-plain flaps, which have slower servo motors that will be used primarily for drag and maneuver load reductions.

∗ Research Aerospace Engineer, Aeroelasticity Branch, MS 340, Member AIAA, josiah.m.waite@nasa.gov.

† Research Engineer, Dynamic Systems and Control Branch, MS 308, Associate Fellow AIAA.

‡ Research Aerospace Engineer, Aeroelasticity Branch, MS 340, Senior Member AIAA.

§ Research Aerospace Engineer, Aeroelasticity Branch, MS 340, Associate Fellow AIAA.

Control laws must be developed for each test (MLA, GLA, AFS, and drag reduction) and must be able to work in tandem when two or more control strategies are being tested simultaneously (e.g., when control laws to mitigate maneuver loads and gust loads operate simultaneously). The development of these control laws requires accurate mathematical representations of the aerodynamics, structural dynamics, and control surface dynamics of the system.

Taken together, this mathematical representation—which is often in state-space form— is the aeroservoelastic (ASE) model.

Three low-order ASE models of differing fidelities were developed for the IAWTM test and will be presented here.

The differences in fidelity appeared in the aerodynamic modeling—one model assumed purely linear aerodynamics, while the other two models used different techniques to linearize the aerodynamics about a nonlinear mean flow.

Furthermore, simple observer-based full-state feedback controllers for MLA and AFS were designed from these low-order ASE models and implemented in a higher-order, nonlinear computational ASE simulation to assess the effectiveness of each controller. Observers were designed from each of the three lower-order ASE models, and each (linear) observer was then used as an estimate of the nonlinear dynamics of the system. Using the estimated states and the full-state feedback controller, control surface commands were fed back to the higher-order, nonlinear ASE system. The insight gained during this development will be applied to control law development for the physical IAWTM wind-tunnel test.

II. Methods In total, four ASE models are discussed in this paper. One is a nonlinear computational ASE framework that was used as the testbed where control laws were implemented and assessed. The other three ASE models are of lower order and were used to generate state-space systems from which control laws were derived. In this section, the structural model used in all four ASE models will be presented. Next, a description of the nonlinear computational ASE framework in which the controllers were implemented will be given. Brief discussions of aerodynamic modeling for each lower-order ASE model will follow. Then, control law development and implementation will be discussed.

A. Structural Dynamics Modeling A schematic of the finite element model (FEM) is shown in Fig. 2. The structural model has a detailed fuselage and wind-tunnel mount modeled by solid plate elements. The FEM treats the spar as an equivalent beam with rigid elements to place concentrated masses and to attach control surfaces. The active control surfaces are modeled as flat plates.

Fig. 2 Finite element structural model. Figure courtesy Boeing Research and Technology.

Table 1 presents the frequencies and descriptions of the first ten structural modes. Because the wing structure is modeled with a sparse distribution of flexible and rigid beam elements, the method of modal displacement is based on zonal mathematical descriptions of the structural model displacements. A similar method was used for the Subsonic Ultra Green Aircraft Research (SUGAR) wind-tunnel model [ 4 ]. Because the CFD model has a symmetry plane at the fuselage centerline, y-displacements of the fuselage that cross the y=0 plane have been made zero in the CFD model.

Table 1 Frequencies and descriptions of the first ten (FEM) structural modes.

mode frequency description number (Hz) st 1 3.61 1 wing vertical bending st st 2 7.65 1 fuselage pitch / wing 1 vertical bending st 3 8.35 1 wing/balance fore/aft bending nd 4 11.08 2 wing vertical bending nd 5 14.98 2 wing / balance fore/aft bending nd 6 15.15 2 wing vertical bending / fuselage pitch st st 7 17.11 1 nacelle pitch / inboard wing 1 torsion 8 24.04 nacelle lateral rd 9 26.78 3 wing vertical bending 10 29.34 nacelle lateral / balance axial For coupling to the aerodynamic models, a linear modal representation of the structural dynamics was formed using 24 structural modes. Inputs to the structural model included generalized masses, modal frequencies, and modal damping.

For the cases presented in this paper, these parameters for the 24 structural modes came from MSC.NASTRAN [5].

B. Nonlinear Computational ASE Model The NASA Langley Research Center FUN3D (Fully Unstructured 3-Dimensional Navier-Stokes) [ 6 , 7 ] code was used to solve the steady and unsteady Reynolds-averaged Navier-Stokes (RANS) equations and couple them to modal structural dynamics equations. FUN3D is primarily a node-centered finite-volume unstructured CFD solver for compressible flows [ 6 , 7 ]. In this work, FUN3D solves the RANS equations using a second-order upwind scheme with the Spalart-Allmaras turbulence model [ 8 ]. FUN3D’s structural dynamics solver uses a second-order accurate predictor-corrector scheme. All aeroelastic FUN3D analyses presented in this paper, with the exception of the FUN3D linearized frequency-domain (LFD) [9] results (which will be discussed later), were solved using these schemes.

The volume grid used in the FUN3D analyses was a semi-infinite flow-field grid with a symmetry plane at the wind-tunnel mounting wall. It had 14,618,916 nodes, 3,598,892 tetrahedral cells, 184,622 pyramid cells and 27,598,875 prism cells. The overall surface mesh distribution is shown in Fig. 3a. It will be noted that there is clustering of nodes in the area of the control surfaces. Figure 3b shows the mesh layout over several of the inboard control surfaces.

A typical FUN3D aeroelastic analysis consists of sequentially obtaining three solutions: (1) the converged steady solution about a rigid wing, (2) the converged unsteady solution of an aeroelastically-deforming body with critical structural damping—this is the static aeroelastic solution, and (3) the converged unsteady solution of an aeroelastically- deforming body with zero, or low, structural damping. Typically, the static aeroelastic simulation is restarted from the converged steady solution about a rigid wing, and it serves as the starting point for the unsteady simulation with low (or zero) structural damping. When performing flutter analyses, a modal disturbance of some type—generally a modal velocity in one of the modes—is introduced in the third simulation. In the nonlinear computational ASE simulations st presented here, the 1 mode was given an initial velocity.

C. Aerodynamics Modeling in Lower-order ASE Models 1. Linearized Aerodynamics using Eigensystem Realization Algorithm from CFD Data ASE state-space models can be generated from experimental data using system identification techniques such as the eigensystem realization algorithm (ERA). To use this technique, prescribed inputs are applied to the system (which, in the current case, may be control surface rotations applied to the physical wind-tunnel model) and outputs of the system are measured. The system identification technique then operates on the set of input-output data to generate a state-space (a) CFD surface mesh with refined control surfaces. (b) CFD refined inboard control surface mesh.

Fig. 3 model of the dynamics of the system, whose inputs and outputs are of the same units as the applied inputs and measured outputs of the physical model during the system identification process. Since, for the IAWTM tests, the physical model will be outfitted with accelerometers along the wing (see Fig. 4), the ASE model derived from the system identification process will be in units of acceleration.

Fig. 4 Visualization of sensor and actuator layout for the wind tunnel model. Figure courtesy Boeing Research and Technology.

ASE state-space models can also be generated from computational simulation data using the same system identification techniques. One advantage to using simulation data is the ability to monitor any variable within the simulation as output. Access to such a wide variety of variables is often infeasible, and sometimes impossible, when using experimental data. For example, in a computational environment, it is relatively easy to apply a control surface rotation to the aerodynamic grid and monitor forces over the surface of the wing. Exercising a system identification algorithm on that input-output data set leads to a reduced-order aerodynamic state-space model without the interaction of the structural dynamics. This is useful since this purely aerodynamic model could then be coupled to a structural dynamics state-space model through a gain representing dynamic pressure, and the result would be an aeroservoelastic state-space model that is valid over a wide range of dynamic pressures. To derive an ASE model that is valid over a similar range of dynamic pressures experimentally, one would either have to instrument the wind tunnel model with additional sensors to measure the purely aerodynamic response, or apply inputs at several dynamic pressures and interpolate the results. This is a major advantage of deriving ASE models, and their associated control laws, from simulation data. Yet, in practice, it is often safer to derive ASE models from experimental data since small errors in the computational model can propagate through to the control laws, causing them to be ineffective. This is especially true in transonic flight where aerodynamic nonlinearities have been notoriously difficult to predict.

The process of generating ASE state-space models from computational simulation data using ERA has been previously demonstrated [ 10 , 11 ]. A similar process was followed in this work. A converged static aeroelastic solution was obtained at a dynamic pressure of 275 psf and at Mach 0.8. Then, several input-output data sets were generated from a series of aerodynamic training simulations. These training simulations were done using the finite-volume, time-marching version of FUN3D. Control modes were generated and included in these training simulations that represented control surface rotations in a manner similar to the structural mode representations. In these training simulations, the displacements and velocities of each structural and control mode were prescribed and the resulting aerodynamic forces induced on every structural mode were calculated.

Once the aerodynamic training simulations were completed, the input-output data, which in this case were modal displacement and modal aerodynamic force data, were passed to the AEROM (AeroElastic Reduced-Order Modeling) software [ 12 ]. AEROM is a suite of tools developed at the NASA Langley Research Center particularly for generating aeroelastic reduced-order models. AEROM is largely based on the System/Observer/Controller Identification Toolbox [ 13 ], or SOCIT, another code developed at the NASA Langley Research Center. Two tools in particular from SOCIT were used in this work: its “pulse” algorithm and its eigensystem realization algorithm (ERA). The “pulse” algorithm [ 13 ] was used to generate individual impulse responses of all inputs on all outputs from the input-output data produced by the training simulations. A brief overview of how this algorithm works in aeroelastic applications is given in Ref. [ 12 ].

A discrete-time, state-space model of the system was then generated via ERA [ 13 , 14 ]. For the present case, SOCIT’s “pulse” and ERA were used to generate aerodynamic state-space models whose inputs were modal displacements (of structural and control modes) and whose outputs were modal aerodynamic forces on the 24 structural modes.

Lastly, these aerodynamic state-space models were combined with a linear representation of the structural dynamics, which was described in Section II.A, resulting in a reduced-order ASE model where the aerodynamic models pass generalized aerodynamic forces (GAFs) to the structural model, and the structural model returns generalized displacements back to the aerodynamic model. This aeroelastic ROM had two sets of inputs: (1) the dynamic pressure, which simply acts as a gain on the GAFs from the aerodynamic models, and (2) the control surface displacements. Note that no control system model was included here. Only the aerodynamic influence of control surfaces was captured, which implies that the control surfaces themselves were assumed to be perfectly rigid with infinite stiffness.

Because the structure and control surfaces are “trained” about the nonlinear static aeroelastic CFD solution, the GAFs computed by this reduced-order model capture the influence of nonlinearities in the mean flow that linear methods will not, which may lead to better accuracy. This is especially important in transonic flow, where shocks on the vehicle can significantly influence the GAFs.

2. Linearized Aerodynamics using Linearized Frequency-Domain CFD In addition to being a node-centered, finite-volume fluid dynamics solver, FUN3D also contains a streamline upwind Petrov-Galerkin (SUPG) stabilized finite-element (SFE) solver [ 15 , 16 ]. Recently, an implementation of the linearized frequency-domain (LFD) approach was added within the SUPG SFE solver in FUN3D [ 9 ]. In this approach, the structural mode shapes are transferred to the aerodynamic surface and a static aeroelastic solution is derived.

For its application to the IAWTM configuration, this static aeroelastic solution was computed with the finite-volume time-marching FUN3D solver at a dynamic pressure of 275 psf and Mach 0.8 (the same static solution used for the ROM). Then, in the SFE solver within FUN3D, the governing equations of the flow are linearized about the (nonlinear) static aeroelastic solution and the responses of a unit displacement perturbation applied to each structural mode are computed. These responses, which are oscillatory surface forces applied on the aerodynamic surface, can be projected to the structural modes to calculate the generalized aerodynamic forces (GAF). This implementation of the LFD approach in the SFE solver within FUN3D will subsequently be referred to as FUN3D/SFE LFD. The advantages of this method over linear methods are similar to those for the ROM: because FUN3D/SFE LFD linearizes the governing equations of the flow about the nonlinear static aeroelastic CFD solution, the GAFs computed by FUN3D/SFE LFD capture the influence of nonlinearities in the mean flow that linear methods will not. This is done through an exact linearization of the flow equations about the static aeroelastic solution, whereas the ROM’s linearization is influenced, to some degree, by the nonlinearity arising from physically perturbing the modes through finite motions in space. Additionally, like the ROM, the FUN3D/SFE LFD approach allows for a single model to be used over a range of dynamic pressures.

After computing GAFs using FUN3D/SFE LFD, a state-space model was generated. A discussion of how a state-space system was generated from LFD-generated GAFs will be included here.

3. Linear Aerodynamics The last ASE model developed here was based on linear aerodynamics and was constructed in ZAERO [ 17 ]; the aerodynamic paneling (shown in Fig. 5) consisted of flat plate elements for the wing and pylon, and body elements for the fuselage and flow-through nacelle. Specifically, the ZONA6 oscillatory subsonic lifting method was used to model the wing, and the BODY7 non-lifting method to model the fuselage and nacelle. The wing and nacelle paneling were splined to the structure with an infinite plate spline, and the motion of the fuselage body was ignored.

A discussion of how control surfaces were modeled in ZAERO will be included here, as well as how the state-space ASE model was generated.

Fig. 5 Aerodynamic paneling used in the ZAERO software.

D. Observer and Control Law Development Full-state feedback controllers were developed for MLA and AFS. Full-state feedback was chosen due to its simple mathematical properties over other control schemes. Full-state feedback involves multiplying the states (or estimates of the states) by some feedback gain matrix, K , to produce some control input vector, u . To do this using the aerodynamic states from FUN3D was not directly possible as the number of states was extremely large (modal structural data, but also velocity, pressure, and density information at each node), and the requisite matrices were not easily accessible. Instead, a small amount of data was taken from FUN3D at each time step (24 modal accelerations, ‹ η , converted to a number ‹ of vertical acceleration measurements along the wing, δ ), and used to estimate the full state vector of a stand-alone linear aeroelastic model. This linear aeroelastic state was then used to produce the commanded control surface rotations, which were fed to FUN3D. The development of these observers and controllers followed the methods presented in Ref.

[ 18 ], and is essentially a linear quadratic Gaussian (LQG) controller. The stand-alone linear aeroelastic observer was integrated in time alongside FUN3D, and was assumed to be of the form: ˆ ¤ ‹ ˆ x = A ˆ x + Bu + L ( ‹ δ − δ ) (1) where u is a vector of control surface rotations and ˆ x is a vector of linear aeroelastic states, generated from the linearized ‹ aeroelastic ROM (or another linear aeroelastic model); δ is the vector of vertical accelerations along the wing, which ˆ ‹ comes from FUN3D, and δ is its estimate generated from the linear aeroelastic model. As previously mentioned, L is an observer gain matrix whose derivation is described below. When L is properly designed, the error between the FUN3D ˆ ‹ ‹ accelerations and the linear accelerations ( δ − δ ) will quickly decay to zero during the unsteady simulation.

Having obtained a compact aeroelastic state, feedback was assumed of the form: u = − K ˆ x

The feedback gain matrix (i.e., the controller), K , was found such that a performance index, J , was minimized: π ∞ ( ) T T J = ˆ x Q ˆ x + u Ru dt (2) ∫ ∫ ∞ ∞ 1 1 T T where ˆ x Q ˆ x dt represents the energy of some combination of the states, and u Ru dt represents the energy 2 0 2 0 necessary for control; Q is a state weighting matrix; R is a control weighting matrix. Larger values of R relative to Q penalize the controls more. The goal was to minimize the energy in the states without using too much energy in the control effort. The L matrix from Equation 1 was generated using a similar equation to the J performance index with R and Q being white noise inputs.

When minimizing Eq. 2 to find a controller, K , the choices for the Q and R matrices have a significant influence on the positions of the closed-loop plant poles. Similarly, when minimizing Eq. 2 to find an observer gain matrix, L , the choices for the Q and R matrices influence the closed-loop observer poles. Thus, there were two sets of Q and R matrices associated with K and L , respectively. The method used here to find Q and R follows that found in Ref. [ 18 ], which was to assume some performance output — which was to be kept small — of the form z = H ˆ x . To find the Q and R matrices associated with the controller, K , z was chosen to be the modal velocities of the structural modes; therefore, H was simply a matrix that isolated the modal velocity terms from the estimated state vector ˆ x and zeroed out all other terms. Then, the performance index was redefined as: π ∞ ( ) T T ¯ J = z Qz + u Ru dt ¯ where Q was a diagonal matrix containing q = and z was some maximum allowable deviation in each i max 2 i z ma x i ¯ element of z for i equal to 1 to the number of elements in z . In Ref. [ 18 ], R is defined in a similar fashion as Q , but for control inputs. However, in order to yield satisfactory behavior from the combined system, R was tuned to ensure that the closed-loop observer poles were faster than the closed-loop plant poles. This was done simply by defining R = ρ I and changing ρ until the desired time response was obtained. In this work, the performance of several controllers ¯ derived from different linear ASE models was to be compared. In order to show a fair comparison, the values for Q and R associated with both K and L were identical across all controllers.

III. Results A. Comparison of Lower-Order ASE Models Figure 6 shows the root loci of the aeroelastic system for each model. A magnified view of this same plot is shown in Fig. 7, which shows that the flutter mechanism predicted by the ROM is mode 4, by FUN3D/SFE LFD is mode 5, and by ZAERO is mode 6. These differences may be caused by variances in the root-tracking methods used in the codes, and it may simply indicate that the flutter mechanism includes a coupling of modes 4, 5, and 6. Figure 7 also includes frequency and damping values extracted from FUN3D (time-domain) time histories for modes 4, 5, and 6 (again, with the matrix pencil method mentioned previously) at four dynamic pressures. These extracted values indicate that the predictions of the linear and linearized models are capturing the frequencies of the critical flutter mechanism well.

B. Controller Effectiveness The final paper will include additional discussion and conclusions regarding the effectiveness of the con- trollers for: • Maneuver load alleviation • Active flutter suppression Fig. 6 Comparison of aeroelastic root loci as predicted by several codes. Note that real values greater than zero indicate instability in the system.

Fig. 7 A magnified comparison of aeroelastic root loci as predicted by several codes. Note that real values greater than zero indicate instability in the system. Also note that the symbols represent an estimate of frequency and damping from time-domain data obtained by FUN3D for modes 4 and 5 at four dynamic pressures ( Q ).

References [1] “Advanced Air Transport Technology (AATT) Project Technical Challenges,” https://www.nasa.gov/aeroresearch/ programs/aavp/aatt/technical-challenges , 2017. Accessed: 2019-10-20.

[2] Brooks, T. R., Kenway, G. K., and Martins, J. R., “Benchmark aerostructural models for the study of transonic aircraft wings,” AIAA Journal , Vol. 56, No. 7, 2018, pp. 2840–2855.

[3] Allen, T. J., Sexton, B. W., and Scott, M. J., “SUGAR Truss Braced Wing Full Scale Aeroelastic Analysis and Dynamically Scaled Wind Tunnel Model Development,” AIAA SciTech Forum , 2015.

[4] Bartels, R. E., Scott, R. C., Allen, T. J., and Sexton, B. W., “Aeroelastic Analysis of SUGAR Truss-Braced Wing Wind-Tunnel Model Using FUN3D and a Nonlinear Structural Model,” AIAA SciTech Conference , 2015.

[5] Rodden, W. P., and Johnson, E. H., MSC/NASTRAN Aeroelastic Analysis: User’s Guide; Version 68 , MacNeal-Schwendler Corporation, 1994.

[6] Biedron, R., Carlson, J.-R., Derlaga, J., Gnoffo, P., Hammond, D., Jones, W., Kleb, B., Lee-Rausch, E., Nielsen, E., Park, M., Rumsey, C., Thomas, J., and Wood, W., “FUN3D Manual: 13.4,” Tech. Rep. TM-2018-220096, NASA, 2018.

[7] Anderson, W. K., and Bonhaus, D. L., “An Implicit Upwind Algorithm for Computing Turbulent Flows on Unstructured Grids,” Computers and Fluids , Vol. 23, No. 1, 1994, pp. 1–22. URL https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.

gov/20040112057.pdf .

[8] Spalart, P. R., and Allmaras, S. R., “One-Equation Turbulence Model for Aerodynamic Flows,” 30th AIAA Aerospace Sciences Meeting and Exhibit , Reno, NV, 1992. https://doi.org/https://arc.aiaa.org/doi/pdf/10.2514/6.1992-439, URL https: //arc.aiaa.org/doi/pdf/10.2514/6.1992-439.

[9] Jacobson, K., Stanford, B., Wood, S., and Anderson, W. K., “Flutter Analysis with Stabilized Finite Elements based on the Linearized Frequency-domain Approach,” AIAA Scitech 2020 Forum , 2020, p. 0403.

[10] Waite, J., Stanford, B., Bartels, R. E., and Silva, W. A., “Reduced Order Modeling for Transonic Aeroservoelastic Control Law Development,” AIAA Scitech 2019 Forum , 2019, p. 1022.

[11] Waite, J., Stanford, B., Bartels, R. E., and Silva, W. A., “Active Flutter Suppression Using Reduced Order Modeling for Transonic Aeroservoelastic Control Law Development,” AIAA Aviation 2019 Forum , 2019, p. 3025.

[12] Silva, W. A., “AEROM: NASA’s Unsteady Aerodynamic and Aeroelastic Reduced-Order Modeling Software,” Aerospace , Vol. 5, No. 2, 2018.

[13] Juang, J.-N., Horta, L. G., and Phan, M., “System/Observer/Controller Identification Toolbox,” Tech. Rep. TM-107566, NASA, 1992.

[14] Juang, J.-N., and Pappa, R. S., “An Eigensystem Realization Algorithm for Modal Parameter Identification and Model Reduction,” Journal of Guidance, Control, and Dynamics , Vol. 8, No. 5, 1985, pp. 620–627.

[15] Anderson, W. K., and Newman, J., “High-Order Stabilized Finite Elements on Dynamic Meshes,” 2018 AIAA Aerospace Sciences Meeting , 2018, p. 1307.

[16] Anderson, W. K., Newman, J. C., and Karman, S. L., “Stabilized Finite Elements in FUN3D,” Journal of Aircraft , Vol. 55, No. 2, 2018, pp. 696–714.

[17] Zona Technology, “ZAERO Version 9.2 User’s Manual,” 2017.

[18] Stevens, B. L., Lewis, F. L., and Johnson, E. N., Aircraft Control and Simulation: Dynamics, Controls Design, and Autonomous Systems, 3rd Edition , John Wiley and Sons, 2015.