Study of a neural network-based system for stability augmentation of an airplane

Study of a neural network-based system for stability augmentation of an airplane Author: Roger Isanta Navarro Report Supervisors: Oriol Lizandra Dalmases Fatiha Nejjari Akhi-Elarab Aeronautical Engineering September 2013 This page intentionally left blank i Contents 1 Aim of the study .................................................................................................... 1 2 Scope .................................................................................................................... 2 3 Justification ........................................................................................................... 3 4 State of the art ...................................................................................................... 4 5 Background ........................................................................................................... 7 5.1 Airplane dynamics .............................................................................................. 7 5.2 Fuzzy Logic, Neural networks and ANFIS ............................................................ 9 5.3 Optimal control and dynamic programming .................................................... 14 6 Neural network structure .................................................................................... 17 6.1 Structure .......................................................................................................... 17 6.2 Behavior ........................................................................................................... 18 6.3 Implementation ................................................................................................ 20 6.4 Code validation ................................................................................................ 20 7 Numerical simulation .......................................................................................... 27 7.1 Fourth order Runge-Kutta method ................................................................... 27 7.2 Implementation ................................................................................................ 27 7.3 Code validation ................................................................................................ 28 8 Stability augmentation system proposal ............................................................. 35 8.1 Optimal control algorithm for offline training .................................................. 35 8.2 Simulation results for offline control ................................................................ 36 8.3 Simulation results for offline control with unstable aircraft............................. 44 8.4 Integration in phugoid PID controller ............................................................... 50 9 Viability of online critic-adaptive stability augmentation system ........................ 55 9.1 Online stability augmentation system proposal ............................................... 55 9.2 Viability of the proposal ................................................................................... 56 10 Conclusions and further work ............................................................................. 58 11 Environmental effects ......................................................................................... 61 12 Budget ................................................................................................................. 62 12.1 Work hours ....................................................................................................... 62 12.2 Used software .................................................................................................. 62 12.3 Used hardware ................................................................................................. 62 ii 13 Source codes ....................................................................................................... 63 14 References .......................................................................................................... 64 iii List of figures Figure 4.1 Example of pitch damper stability augmentation system .............................. 4 Figure 5.1 Artificial neuron model [2] ............................................................................. 9 Figure 5.2 Membership evaluation of traditional logic ................................................. 10 Figure 5.3 Membership evaluation of fuzzy logic using bell functions ......................... 11 Figure 5.4 Membership evaluation of fuzzy logic using linear functions ...................... 11 Figure 5.5 Diagram representation of a 3-inputs 3-rules ANFIS network ..................... 12 Figure 5.6 Example of optimal path from to .......................................................... 15 Figure 5.7 Backwards dynamic programming example to find optimal path from to ......................................................................................................................................... 15 Figure 6.1 Diagram representation of a 2-inputs 2-rules network................................ 17 Figure 6.2 Example of partial activation of network paths ........................................... 19 Figure 6.3 Test approximation to single-variable function using linear membership functions .......................................................................................................................... 21 Figure 6.4 Mean Square Error vs. number of membership functions in test approximation to single variable function using linear membership functions ............... 21 Figure 6.5 Training time vs. number of membership functions in test approximation to single variable function using linear membership functions ............................................ 22 Figure 6.6 Two-variable test approximation function................................................... 23 Figure 6.7 Test approximation to two-variable function using 5 linear membership functions per input ........................................................................................................... 23 Figure 6.8 Test approximation to two-variable function using 9 linear membership functions per input ........................................................................................................... 24 Figure 6.9 Test approximation to two-variable function using 25 linear membership functions per input ........................................................................................................... 24 Figure 6.10 Two-variable test approximation function for online training ................... 25 Figure 6.11 Online training test approximation to two-variable function using 25 linear membership functions per input ..................................................................................... 26 Figure 7.1 Phugoid oscillation for evaluation of numerical integration method .......... 29 Figure 7.2 Short period oscillation for evaluation of numerical integration method ... 30 Figure 7.3 response to elevator ( ) ............................................................. 31 Figure 7.4 Theoretical response to elevator ( ) .......................................... 31 Figure 7.5 response to elevator ( ) ............................................................. 32 Figure 7.6 Theoretical response to elevator ( ) .......................................... 32 Figure 7.7 response to elevator ( ) ............................................................... 33 Figure 7.8 Theoretical response to elevator ( ) ............................................. 33 Figure 7.9 response to elevator ( ) ............................................................... 34 iv Figure 7.10 Theoretical response to elevator ( ) ........................................... 34 Figure 8.1 Discretization example and optimal control algorithm behavior ................. 36 Figure 8.2 Block diagram of offline stability augmentation system (training phase) .... 37 Figure 8.3 Block diagram of offline stability augmentation system (application phase) ......................................................................................................................................... 37 Figure 8.4 Simulation results of short period offline stabilization augmentation system for Boeing 747. Response to initial vertical speed ( ) ( ). ........... 38 Figure 8.5 Simulation results of short period offline stabilization augmentation system for Boeing 747. Response to initial vertical speed ( ) Variation of R. ( ). .. 39 Figure 8.6 Simulation results of short period offline stabilization augmentation system for Boeing 747. Response to commanded ( ) ( ). ............. 40 Figure 8.7 Full system simulation results of short period offline stabilization augmentation system for Boeing 747. Response to commanded ( ) ( ). Short time simulation. ................................................................... 42 Figure 8.8 Full system simulation results of short period offline stabilization augmentation system for Boeing 747. Response to commanded ( ) ( ). Large time simulation. ................................................................... 44 Figure 8.9 Unstable behavior of modified Boeing 747 with no control action. ............ 45 Figure 8.10 Simulation results of short period offline stabilization augmentation system for Boeing 747. Response to initial vertical speed ( ) ( ). ......................................................................................................................................... 46 Figure 8.11 Simulation results of short period offline stabilization augmentation system for unstable Boeing 747. Response to commanded ( ) ( ). ................................................................................................................. 47 Figure 8.12 System simulation results of short period offline stabilization augmentation system for unstable Boeing 747. Response to commanded ( ) ( ). ....................................................................................................... 49 Figure 8.13 Block diagram of offline stability augmentation system integrated with a PID controller. (Application phase) .................................................................................. 50 Figure 8.14 Full system simulation results of short period offline stabilization augmentation system for unstable Boeing 747 with phugoid PID control. Response to commanded ( ) ( ). Short time simulation ....................... 52 Figure 8.15 Full system simulation results of short period offline stabilization augmentation system for unstable Boeing 747 with phugoid PID control. Response to commanded ( ) ( ). Large time simulation. ...................... 54 Figure 9.1 Block diagram representation of a critic-adaptive control system .............. 55 Figure 9.2 Block diagram of offline training phase in online stability augmentation system .............................................................................................................................. 56 v Figure 9.3 Block diagram of online stability augmentation system application ............ 56 vi List of tables Table 7.1 Validation data for numerical integration method [1] .................................. 29 Table 12.1 Budget breakdown ...................................................................................... 62 vii List of abbreviations Airplane dynamics: State matrix Input matrix Wing span  Mean aerodynamic chord Drag coefficient Lift coefficient Elevator deflection Throttle position Pitch angle Longitudinal speed variation Moments of inertia Products of inertia Moments in body axes Aircraft mass Balance rate Pitch rate Yaw rate Wing area Subscript Refers to lineal speed Subscript Refers to lineal speed Subscript Refers to linear speed Subscript Refers to angular speed Subscript Refers to angular speed Subscript Refers to angular speed Subscript Refers to control action Subscript Refers to elevator deflection Subscript Refers to throttle position Longitudinal speed Control vector Lateral speed 5 Usually, the stability augmentation systems are integrated in the Flight Management System and its operation varies according to the manufacturer and the model. Thus, the European manufacturer Airbus, whose majority of airplane models use the Fly-by-wire technology, incorporate in its signal the action of the stability augmentation systems, among others. The American manufacturer Boeing, in some of its models, opts for a different approach: not allowing a direct actuation of these systems on the aircraft control surfaces, but designing the system to restrict control authority beyond certain range by increasing the back pressure once the desired limit has been reached. This is achieved by electronically controlled actuators according to certain protections and augmentation systems among which the stability augmentation system are included. The stability augmentation systems base their operation on PID controllers; however, because of their adaptive capability and the experience that has been collected over the last decades in other control systems, there is a growing interest in developing stability augmentation systems using neural networks. These systems are, theoretically, able to control the evolution of all the variables that govern a certain behavior of the airplane while permitting real-time adaptability, allowing greater control and stability optimization. The term “theoretically” has been added because this idea, and subsequently this project, deals with a problem whose ideal solution, to this day, has not yet been implemented. By ideal solution, the author refers to the one that can simultaneously control four state variables governing the longitudinal motion of an aircraft, by using neural networks. Neural networks grow exponentially with the number of inputs of the problem and with the accuracy or discretization performed on each of them. A neural network with four inputs ( and ) and good discretization, results in a structure too large to be trained in a reasonable time. Offline training (prior to the inclusion of the network in the system) can take months, and obviously, given these times of training, a real-time online training is not even worth to be considered. If to these inputs one would want to add other data like the airplane altitude, flight speed or other aspects of the flight that might be taken into consideration, the problem becomes untreatable. Given this current inability to develop this ideal system, there are many research groups working on the development of alternative methodologies. On the contrary, control systems of one or two variables are commonly used in a variety of applications such as, for example, control of the temperature variation in a furnace, the control of flow through a conduit or just the traditional inverted pendulum problem. 6 Among the research groups that have developed or are developing multivariable control systems using neural networks in the aeronautical sector, one article [5] stands out, presented by S. Ferrari and R. F. Stengel in which an adaptive critical system for the control of a six degrees of freedom simulation of a jet business is developed. As is the case of this study, the authors of the article use the techniques of adaptive critic systems as base for their proposed solution. This article, even though its complexity is greater than that encompassed by this study, has been a source of inspiration when developing it. [5] presents a nonlinear control system consisting on a network of networks, in which the neural network controller is algebraically trained offline by the observation that its gradients must equal corresponding gain matrices at some chosen operating points. The online training is performed by a dual heuristic adaptive critic architecture capable of optimizing the performance incrementally over time by accounting for plant dynamics, unmodeled dynamics, parameter variations and control failures. 7 5 Background The following section aims to provide a basic understanding of the main tools used throughout this study. Only a general description is provided here in the report section. More developed explanations can be found in the corresponding annexes and references. For the development of this project, knowledge of three main areas has been necessary. Firstly, the aim being the study of a stability augmentation system for aircraft, it is essential to develop the equations governing the airplane dynamics. Secondly, one of the main basic requirements of this project is the use of neural networks for the functioning of stability augmentation system. This has required studying important disciplines of artificial intelligence such as fuzzy logic, neural networks, and more specifically, the ANFIS neural networks. Finally, the proposed solution is an optimal control system, therefore it is necessary to understand the operation of these systems as well as the dynamic programming and their primary development tool. 5.1 Airplane dynamics To develop the training and simulation programs it is necessary to obtain a system of equations of the form: 󰇗 (5.1) where and are the state variables and control vectors, respectively, and and are the state and input matrices, respectively. Generally, the state variables from which the system is developed are, arranged as shown in the Equation (5.2): longitudinal speed variation , vertical speed , aircraft pitch rate and pitch angle . The control variables are the elevator deflection and the throttle position . The system of equations is obtained after linearization of the kinematic and dynamic equations of motion of the airplane. Detail of the obtaining procedure and linearization of these equations is given in Annex 2 - Airplane Dynamics, as well as the definition of all the variables appearing. The linearized system of equations for the airplane’s longitudinal motion yields to: 8 { 󰇗 󰇗 󰇗 󰇗} [ 󰇗 󰇗 󰇗 󰇗 [ 󰇗 󰇗] [ 󰇗 󰇗] [ 󰇗( ) 󰇗] 󰇗 ( 󰇗) ] { } { 󰇗 󰇗 󰇗 } (5.2) The control terms include within their definition the elevator deflection and throttle position . In the case concerning this study, however, the throttle position will not be considered for airplane control for the high inertia of the plane against the thrust results in a response with a period too large to control its longitudinal stability. The control vector can be written, therefore, { 󰇗 󰇗 󰇗 } (5.3) Since a short period stability augmentation system is to be designed, the simplified equations for this type of motion are also presented. These are obtained by neglecting the X-force equation and setting , for the speed is substantially constant while the airplane pitches relatively rapid, and taking into consideration that 󰇗 is small compared to and is also small when compared to the product . After a final consideration of , Equation (5.2) yields to: [ 󰇗 󰇗] [ [ 󰇗 ] [ 󰇗 ] ] [ ] (5.4) 9 5.2 Fuzzy Logic, Neural networks and ANFIS The approximator tasks functions, interpolation, online learning and adaptability that appear in this study fall largely on mathematical-logical structures called neural networks. Neural networks mimic in a mathematical form the neuronal function at a level of synaptic relationship between neurons or nodes, creating kinship relationship (parent - child) through which values are propagated. This achieves converting input signals to output after a few transformations either predefined or acquired through learning or training. The basic structural model of an artificial neuron is presented below: Figure 5.1 Artificial neuron model [2] The operation of an artificial neuron is divided into two main phases: the algebraic sum of the inputs from the synapse and evaluation of a linear function ( ), which results in the output of the neuron. Within the sum of inputs a threshold value may be included, whose synaptic weight is equal to one, which allows the modulation of the activation potential of the neuron. This model can be summarized by the following formula: ( ) (∑ ( ) ( ) ( )) (5.5) The activation function ( ) can be defined in very different ways depending on the interest and the specific role of the neuron in the network. Commonly, linear functions, linear within a range and constant outside, with Gaussian bell shape or sigmoidal functions are widely used. 10 For the development of this study ANFIS type networks [2] (Artificial Neuro-Fuzzy Inference System) have been chosen, which have a much more rigid structure than other networks such as a multilayer perception (MLP), but, given the type of problem and considering a same allowable error in the output, require less neural complexity and achieve better results in terms of execution and training time. ANFIS networks mix the bases of neural networks with fuzzy logic taking advantage of both technologies. Fuzzy logic is a probabilistic form of logic which rather than using binary subdivisions (true - false), evaluates the membership of an input to continuous subsets with values ranging from to . Traditionally, computer algorithms define strict boundaries between subsets. Thus, to determine if a certain temperature is considered cold, warm or hot boundary values are set that separate these subgroups. Figure 5.2 Membership evaluation of traditional logic In the previous image two clear boundaries are established for temperatures of and . Algorithmically, the implementation of such logic is very simple, using conditional statements such as if-else the belonging of an event to one of the three subsets can be quickly set and acts accordingly. Instead of establishing strict boundaries and define membership in a subset by a binary assertion, fuzzy logic gives membership degrees or "truth values" of the event with respect to each membership function. These values correspond to bell-type functions, linear, or the ones the user considers appropriate. In Figures 5.3 and 5.4 the 0 0.2 0.4 0.6 0.8 1 1.2 -20 0 20 40 60 80 100 𝑇 [º𝐶] Cold Warm Hot 11 membership functions for the above example are shown, using, first, bell-type functions, and linear, in the second. Figure 5.3 Membership evaluation of fuzzy logic using bell functions Figure 5.4 Membership evaluation of fuzzy logic using linear functions Therefore, fuzzy logic does not establish a belonging to a single subset, but defines partial belongings to all or various subsets within the range of the variable. This leads to the fact that the output of the system is not defined for a particular situation, but it is a 0 0.2 0.4 0.6 0.8 1 1.2 -20 0 20 40 60 80 100 𝑇 [º𝐶] Cold Warm Hot 0 0.2 0.4 0.6 0.8 1 1.2 -20 0 20 40 60 80 100 𝑇 [º𝐶] Cold Warm Hot 12 mixture of all solutions in terms of the different degrees of membership of the input. This reasoning is much more similar to the one taking place in human mind: the temperature of the sample may not be cold or hot, but between cold and warm. The parameters that define the membership functions may also be variable, thereby changing over time the degree of membership of an input to some subsets, and providing the system with learning capacity and adaptability. The basic structure of an ANFIS network is as follows: Figure 5.5 Diagram representation of a 3-inputs 3-rules ANFIS network The network consists of five layers. The first evaluates the membership of each input to the corresponding membership functions, i.e. performs the fuzzy logic part of the process. For each path is assigned a membership factor resulting from the combination of membership functions and evaluates, in the fourth layer, the input values according to this membership coefficient. Second and third layers are the ones charged of determining the factor. The fifth layer only performs an addition of the approximations calculated by each path and returns the result as output value of the network. Obviously, the number of entries is variable and so is the number of paths. The greater the number of paths is, the larger the accuracy capability of the network, i.e. the smaller the error. However, increasing the complexity or the number of inputs incurs in an exponential growth in runtime. 13 An extremely important aspect within ANFIS neural networks is the training method of these, for this defines the convergence speed of the network. Before reviewing this aspect of ANFIS networks, it is worth mentioning which are the variable parameters of the structure, that is, which are trainable. Neurons in layers and contain trainable parameters. The output of the neurons of the first layer results from the evaluation of the degree of membership of the input value to a membership function. These membership functions can vary their parameters over time to allow for greater adaptability of the network, and to better approximate the function. Typically, the membership functions are bell-shaped and are defined according to the expression: ( ) | | (5.6) where parameters , and define the width of the bell function, its maximum slope, and its center, respectively. These parameters are called premise parameters and are trainable. The output of the network is not linear with respect to them, which will determine its training. Regarding the fourth layer, in it the value of the inputs is calculated on a hyperplane defined by the so-called consequence parameters. These hyperplanes modeled after the belonging to different membership functions are precisely the ones who allow the network to approximate nonlinear functions. The expression that is evaluated in a fourth-layer neuron is: ( )   (∑ ) (5.7) where  is the relative weight of the path and is the number if input variables. The consequence parameters are and . They are also trainable and in this case, the output value of the network is linear with respect to them. This last fact is of great significance in the ANFIS networks when training. The training of the network seeks to minimize the Mean Squared Error of it against a training set, i.e., versus a group of input-output pairs. In a general network, where all layers of all neurons have trainable synaptic weights, the training method usually performed is called backpropagation training; consisting in propagating the error upstream by means of partial derivatives of the output with respect to the different network parameters. The ANFIS network, however, because it presents a different and 14 much more defined structure has certain training alternatives that can accelerate its training process. Given the non-linearity of the output with respect to the premise parameters, to minimize the mean square error one can use optimization algorithms such as the Steepest Descent Method [2]. To do so, one must only obtain expressions for the partial derivatives of the output of the network directly with respect to these premise parameters, thus obtaining the gradient that will allow, through a series of iterations, the optimization of these parameters. To train the consequence parameters the fact that the output is linear with respect to them can be exploited. This allows the use of algorithms such as Least-Squares Method or its recursive version, much more rapid and which only requires one iteration to obtain optimized values. These two training methods can be used together in what is called hybrid training, consisting of alternative executions of both until an acceptable Mean Square Error is obtained. It is also possible, as it is the case in this study, that based on the actual use of the network within the system it might be more interesting to use only one method (the steepest-descent method is also applicable to the optimization of the parameters accordingly ), or to perhaps train only one type of the parameters. Regarding this, the training of the consequence parameters it is much more essential than the training of the premise ones; the former largely determine the resulting value of the approach, while the latter only allows for greater accuracy. 5.3 Optimal control and dynamic programming The chosen approach to meet the requirements of this project is based on the design of an optimal control system [3]. Such systems seek to minimize a cost function (or maximize a reward function) defined after the temporal evolution of the state and control variables: ( ) ( ) ∫ ( ) (5.8) Equation which, when discretized, yields to: ∑( ) (5.9) where , and are matrices whose values determine the weight of the state and control variables over the total calculated cost. Thus, modifying , and one can 21 Figure 6.3 Test approximation to single-variable function using linear membership functions It can be seen that both the approach using and membership functions are not acceptable. Instead, the membership functions approximation is very good. Figures 6.4 and 6.5 how, first, the variation of Mean Square Error, and second, the training time, compared to the number of membership functions. Figure 6.4 Mean Square Error vs. number of membership functions in test approximation to single variable function using linear membership functions -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 -1 -0.5 0 0.5 1 𝑓(𝑥) 𝑥 f(x) Lin 3 Lin 4 Lin 5 1.E-04 1.E-03 1.E-02 1.E-01 1.E+00 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 𝑀𝑆𝐸 #𝑀𝐹𝑠 22 Figure 6.5 Training time vs. number of membership functions in test approximation to single variable function using linear membership functions Note, first, that the ordinate axis of Figure 6.4 is presented in logarithmic scale. The Mean Square Error decreases considerably until membership functions are reached. In general, increasing them continues to decrease the error, although not in a so remarkable way. In terms of training time, this increases until the number of membership functions is high enough to quickly adapt to the function to approximate. It should be noted that the behavior observed in these two graphs is closely related to the function to approximate, and will vary over other cases. The structure presented does however, in general, have a better behavior when the number of membership functions is increased (within a logical limit), both regarding execution times and in the error obtained. Next, the network is validated against the following two variable function: ( ) [ ] (6.4) It is important to check not only that the network is able to perform an approximation to a function, but that it accepts multiple inputs (variables). The function to approximate is presented in Figure 6.6. 0 2 4 6 8 10 12 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 𝑡 [𝑠] #𝑀𝐹𝑠 23 Figure 6.6 Two-variable test approximation function The representations of the approximations obtained with , and membership functions per input follow: Figure 6.7 Test approximation to two-variable function using 5 linear membership functions per input -8 -6 -4 -2 0 2 4 6 8 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 24 Figure 6.8 Test approximation to two-variable function using 9 linear membership functions per input Figure 6.9 Test approximation to two-variable function using 25 linear membership functions per input Just as it happens with the single-variable function, the training time for the case of membership functions ( ) is lower than in the cases of and functions ( and respectively). In the case of functions a Mean Square Error of is acchieved. -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 25 Finally, the online adaptability of the network is tested. A network trained with the previous inputs will be presented with a slightly modified input consisting in the addition of a three dimensional bell centered in point ( ) ( ): ( ) ( ) ( ) (6.5) The ANFIS Network will be forced to load the saved data from the previous test and an online training will be performed of the new expression: ( ) ( ) ( ) (6.6) The new modified function and the online training result follow. Notice how the central farthest peak is increased due to the addition of the bell function. Figure 6.10 Two-variable test approximation function for online training The training result again matches the validation function. The online training phase took to be completed. This online training time will always be similar to the offline time if the modification over the training data is important, for the dimensionality of it does not change. However, if the changes are smaller, the online training time also becomes smaller. -8 -6 -4 -2 0 2 4 6 8 10 26 Figure 6.11 Online training test approximation to two-variable function using 25 linear membership functions per input -8 -6 -4 -2 0 2 4 6 8 10 27 7 Numerical simulation To perform the numerical simulation of the differential equation governing the behavior of the system (Equation (6.2)) or any other system of the same type, a Runge-Kutta method of fourth order will be used. Its implementation, as well as the one of the neural network, has been done as a C++ class, and is provided as a Dynamic Link Library (DLL). The class, called RK4, is capable of simulating a system of variables given the matrices and , an initial state vector and a control vector. Although in this study only one control parameter will be considered, the class has been created looking forward to a possible extension or different uses, and allows considering a control parameter vector and hence is a matrix. 7.1 Fourth order Runge-Kutta method As a reminder the algorithm used is briefly described. Given an equation of the type: 󰇗 ( ) (7.1) The Runge-Kutta method allows to determine the state vector at time , known the state at . Using a fourth-order numerical scheme: ( ) ( ) ( ) (7.2) where { ( ( ) ) ( ( ) ) ( ( ) ) ( ( ) ) (7.3) 7.2 Implementation Similarly as occurs with the neural network, the numerical integration algorithm is presented as a C++ class, so instead of a continuous execution, the code responds to calls from the code that uses the said class. 28 The class is named RK4, since it is an implementation of the fourth-order Runge Kutta method. The constructor of the class, i.e. the command used to create the type RK4 variable, supports two variables and which are the number of state and control variables of the system, respectively. Consequently, they also determine the dimensions of the matrices and , of Equation (5.1). To work correctly this code requires the matrices , and time differential over which the integration is to be performed: . Thus there are three functions whose call is mandatory: void defineA(double **A), void defineB(double **B) and void defineDeltaT(double deltaT). Through the function called NextStep(double *X, double *U, double *Xnext), providing a state vector and a control vector , the vector corresponding to ( ) is obtained. Function operation is simple, for it basically implements the Equation (7.2). To do so it makes use of two private functions which determine the values of the four ’s defined in Equation (7.3). Besides determining the next state vector ( ), two more public functions have been defined, which in addition to compute vector ( ) calculate the gradient of ( ) and the gradient of ( ). These functions are called, respectively, void dxdt(double *X, double *U, double *Xnext, double *grad) and void dxdtNext(double *X, double *U, double *Xnext, double *grad). These functions are necessary for the optimal control algorithm that will be discussed later. 7.3 Code validation To validate the simulation code, the behavior of the system will be compared to the solution proposed in [1], when provided with the stability derivatives of a Boeing 747, which can also be found in the same reference. The same stability derivatives and other characteristics of the B747 can be found in Annex 4 – Data for Boeing 747. First, the behavior of the open loop system without control actions will be studied. Given the stability derivatives and other characteristics of the airplane, the code returns a matrix identical to the proposed in [1] but using units of the International System, which allows you to check that Equation (5.2) is correct. It should therefore be expected that if the system is the same, and if the numerical simulation code is correct, the results are the same. The matrix obtained in International System units is: 29 [ ] (7.4) Next, the following initial perturbations introduced: firstly, in first, later in , to check if the frequency and attenuation of the two modes correspond to the ones given in the literature. The simulations were performed with . Reference values which will serve to validate the code follow: Period: [ ] [ ] Phugoid Short-Period Table 7.1 Validation data for numerical integration method [1] For the case of the unitary initial disturbance in , one obtains a large oscillation period and a low attenuation characteristic of a phugoid mode (Figure 7.1). Figure 7.1 Phugoid oscillation for evaluation of numerical integration method The maximums identified in the graph correspond to the points ( ) and ( ). The period between the two peaks is thus , a value which practically coincides with that proposed in the literature. The logarithmic decay between the two peaks is: -1 -0.5 0 0.5 1 1.5 0 100 200 300 400 500 600 Δ𝑢 [𝑚∕𝑠] 𝑡 [𝑠] 30 (7.5) Therefore, (7.6) Value that also resembles the one used to validate. The results corresponding to short period mode is presented in Figure 7.2: Figure 7.2 Short period oscillation for evaluation of numerical integration method The minimums identified in the graph correspond to the points ( ) and ( ). The period between the two peaks is thus , a value which also practically coincides with that proposed in the literature [1]. The logarithmic decay between the two peaks is, in this case: (7.7) Hence, (7.8) The attenuation obtained is higher than the theoretical one. This is a consequence of the influence of the phugoid mode on the vertical direction, which is smaller than the one of short period mode, yet large enough for its effects to be observable. A smaller -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 0 2 4 6 8 10 12 14 16 18 𝑤 [𝑚∕𝑠] 𝑡 [𝑠] 37 The corresponding system structure is shown in the two following pictures. Figure 8.2 shows the structure corresponding to the offline training phase. Within it, the optimal control uses the plant (or a simulation in this case) to obtain the system response for every pair ( ) and every value of and uses a neural network as the discretized space shown in Figure 8.1. In parallel, the policies leading to the optimized cost are stored in a second neural network, which will be the one used to control the system. Figure 8.2 Block diagram of offline stability augmentation system (training phase) Once the policies neural network has been trained, in the second phase it can be used for offline stabilization of the airplane. Figure 8.3 shows a block diagram of the offline stabilization augmentation system proposed. The system receives an input representing the desired elevator deflection , which is processed by the policies network and returns the real deflection to be applied in order to, in time, achieve the desired deflection, while following the optimal path calculated in the previous phase. The neural network requires also the values of and , which are obtained from the airplane response. Figure 8.3 Block diagram of offline stability augmentation system (application phase) Some results obtained with this system are shown in the following pages. First of all, a simple simulation of the short period response to an initial value of ⁄ is shown. The simulation corresponds to a choice of [ ] and . This choice is one of the best ones at optimizing the behavior, that is the angle of attack of the airplane (see Figure 8.4). Next, four simulations of the same example are presented simultaneously in Figure 8.5, but for different values of , i.e. giving more or less cost-weight to the control action. As can be seen, a value of , provides one of the best convergences of the system. 38 Figure 8.4 Simulation results of short period offline stabilization augmentation system for Boeing 747. Response to initial vertical speed ( ⁄) ( [ ] ). -2.0E+00 0.0E+00 2.0E+00 4.0E+00 6.0E+00 8.0E+00 1.0E+01 1.2E+01 0 5 10 15 20 𝑤 [𝑚∕𝑠] 𝑡[𝑠] w -1.6E-02 -1.4E-02 -1.2E-02 -1.0E-02 -8.0E-03 -6.0E-03 -4.0E-03 -2.0E-03 0.0E+00 2.0E-03 0 5 10 15 20 𝑞 [𝑟𝑎𝑑∕𝑠] 𝑡[𝑠] q -1.4E-02 -1.2E-02 -1.0E-02 -8.0E-03 -6.0E-03 -4.0E-03 -2.0E-03 0.0E+00 2.0E-03 4.0E-03 0 5 10 15 20 δe [𝑟𝑎𝑑] 𝑡[𝑠] δe 39 Figure 8.5 Simulation results of short period offline stabilization augmentation system for Boeing 747. Response to initial vertical speed ( ⁄) Variation of R. ( [ ]). -2.0E+00 0.0E+00 2.0E+00 4.0E+00 6.0E+00 8.0E+00 1.0E+01 1.2E+01 0 5 10 15 20 𝑤 [𝑚∕𝑠] 𝑡[𝑠] w (R = 100) w (R = 10) w (R = 1) w (R = 0.1) -3.5E-02 -3.0E-02 -2.5E-02 -2.0E-02 -1.5E-02 -1.0E-02 -5.0E-03 0.0E+00 5.0E-03 0 5 10 15 20 𝑞 [𝑟𝑎𝑑∕𝑠] 𝑡[𝑠] q (R = 100) q (R = 10) q (R = 1) q (R = 0.1) -6.0E-02 -4.0E-02 -2.0E-02 0.0E+00 2.0E-02 4.0E-02 6.0E-02 8.0E-02 1.0E-01 0 5 10 15 20 δe [𝑟𝑎𝑑] 𝑡[𝑠] δe (R = 100) δe (R = 10) δe (R = 1) δe (R = 0.1) 40 Figure 8.6 Simulation results of short period offline stabilization augmentation system for Boeing 747. Response to commanded ( ) ( [ ] ). -2.0E+01 -1.8E+01 -1.6E+01 -1.4E+01 -1.2E+01 -1.0E+01 -8.0E+00 -6.0E+00 -4.0E+00 -2.0E+00 0.0E+00 0 5 10 15 20 𝑤 [𝑚∕𝑠] 𝑡[𝑠] w (with SAS) w (without SAS) -6.0E-02 -5.0E-02 -4.0E-02 -3.0E-02 -2.0E-02 -1.0E-02 0.0E+00 0 5 10 15 20 𝑞 [𝑟𝑎𝑑∕𝑠] 𝑡[𝑠] q (with SAS) q (without SAS) 0.0E+00 1.0E-02 2.0E-02 3.0E-02 4.0E-02 5.0E-02 6.0E-02 7.0E-02 0 5 10 15 20 δe [𝑟𝑎𝑑] 𝑡[𝑠] δe (with SAS) δe (without SAS) 41 The third set of graphs (Figure 8.6), presents the response to a commanded . As can be seen the system adjusts in order to stabilize the aircraft, while reaching the commanded deflection value after the stabilization. Next, a set of graphs show the full system response to the same commanded . This example simulates the whole four state variables system while controlling and . The first five graphs and the last five are equal, the later showing a larger time scale for better appraisal of the phugoid mode. The results shows how the short period mode is stabilized as were the previous examples, while the characteristic response of phugoid mode is maintained. -1.0E+01 0.0E+00 1.0E+01 2.0E+01 3.0E+01 4.0E+01 5.0E+01 6.0E+01 7.0E+01 8.0E+01 0 5 10 15 20 Δ𝑢 [𝑚∕𝑠] 𝑡[𝑠] Δu (with SAS) Δu (without SAS) -2.0E+01 -1.8E+01 -1.6E+01 -1.4E+01 -1.2E+01 -1.0E+01 -8.0E+00 -6.0E+00 -4.0E+00 -2.0E+00 0.0E+00 0 5 10 15 20 𝑤 [𝑚∕𝑠] 𝑡[𝑠] w (with SAS) w (without SAS) 42 Figure 8.7 Full system simulation results of short period offline stabilization augmentation system for Boeing 747. Response to commanded ( ) ( [ ] ). Short time simulation. -6.0E-02 -5.0E-02 -4.0E-02 -3.0E-02 -2.0E-02 -1.0E-02 0.0E+00 1.0E-02 2.0E-02 0 5 10 15 20 𝑞 [𝑟𝑎𝑑∕𝑠] 𝑡[𝑠] q (with SAS) q (without SAS) -4.0E-01 -3.0E-01 -2.0E-01 -1.0E-01 0.0E+00 1.0E-01 2.0E-01 3.0E-01 0 5 10 15 20 Δ𝜃 [𝑟𝑎𝑑] 𝑡[𝑠] Δθ (with SAS) Δθ (without SAS) 0.0E+00 1.0E-02 2.0E-02 3.0E-02 4.0E-02 5.0E-02 6.0E-02 7.0E-02 0 5 10 15 20 δe [𝑟𝑎𝑑] 𝑡[𝑠] δe (with SAS) δe (without SAS) 43 -1.0E+01 0.0E+00 1.0E+01 2.0E+01 3.0E+01 4.0E+01 5.0E+01 6.0E+01 7.0E+01 8.0E+01 0 100 200 300 400 500 600 Δ𝑢 [𝑚∕𝑠] 𝑡[𝑠] Δu (with SAS) Δu (without SAS) -2.0E+01 -1.8E+01 -1.6E+01 -1.4E+01 -1.2E+01 -1.0E+01 -8.0E+00 -6.0E+00 -4.0E+00 -2.0E+00 0.0E+00 0 100 200 300 400 500 600 𝑤 [𝑚∕𝑠] 𝑡[𝑠] w (with SAS) w (without SAS) -6.0E-02 -5.0E-02 -4.0E-02 -3.0E-02 -2.0E-02 -1.0E-02 0.0E+00 1.0E-02 2.0E-02 3.0E-02 4.0E-02 0 100 200 300 400 500 600 𝑞 [𝑟𝑎𝑑∕𝑠] 𝑡[𝑠] q (with SAS) q (without SAS) 44 Figure 8.8 Full system simulation results of short period offline stabilization augmentation system for Boeing 747. Response to commanded ( ) ( [ ] ). Large time simulation. 8.3 Simulation results for offline control with unstable aircraft The results obtained by simulating an aircraft model Boeing 747, in which the stability derivative has been changed to render its behavior unstable, are represented in this section. Specifically, has been increased from an original value of to a positive one: . Figure 8.8 shows two graphs with the response of the open loop system without control action. In the first graph the divergent behavior of short period mode. The second graph, which shows a larger time simulation clearly shows the divergence of the system. The structure of the stability augmentation system applied in this case is unchanged from the previous case, detailed in Figures 8.9 and 8.10. -4.0E-01 -3.0E-01 -2.0E-01 -1.0E-01 0.0E+00 1.0E-01 2.0E-01 3.0E-01 0 100 200 300 400 500 600 Δ𝜃 [𝑟𝑎𝑑] 𝑡[𝑠] Δθ (with SAS) Δθ (without SAS) 0.0E+00 1.0E-02 2.0E-02 3.0E-02 4.0E-02 5.0E-02 6.0E-02 7.0E-02 0 100 200 300 400 500 600 δe [𝑟𝑎𝑑] 𝑡[𝑠] δe (with SAS) δe (without SAS) 45 Figure 8.9 Unstable behavior of modified Boeing 747 with no control action. Within the following pages the results obtained with the offline system, when presented with the same situation as in previous sections with the stable model, are shown. First, the response to an initial value of ⁄ and values of [ ] are presented in Figure 8.10. The results for the same condition with a commanded deflection are shown in Figure 8.11. The results show how the system is able to stabilize the airplane. As it may be appreciated the required action to control the unstable airplane is greater, and the stabilization is a bit slower; nevertheless, this is a result to be expected. It is worth mentioning that the training of the unstable aircraft model has required a higher resolution. Up to membership functions have been assigned to both the and variables and divisions to the values within the optimal control algorithm. The training time for this set up has resulted in over . -1.0E+02 -8.0E+01 -6.0E+01 -4.0E+01 -2.0E+01 0.0E+00 2.0E+01 4.0E+01 6.0E+01 8.0E+01 1.0E+02 0 5 10 15 20 𝑤 [𝑚∕𝑠] 𝑡[𝑠] w -6.0E+06 -4.0E+06 -2.0E+06 0.0E+00 2.0E+06 4.0E+06 6.0E+06 8.0E+06 1.0E+07 020 40 60 80 100 𝑤 [𝑚∕𝑠] 𝑡[𝑠] w 46 Figure 8.10 Simulation results of short period offline stabilization augmentation system for Boeing 747. Response to initial vertical speed ( ⁄) ( [ ] ). -2.0E+01 -1.5E+01 -1.0E+01 -5.0E+00 0.0E+00 5.0E+00 1.0E+01 1.5E+01 2.0E+01 0 5 10 15 20 𝑤 [𝑚∕𝑠] 𝑡[𝑠] w (with SAS) w (without SAS) -2.0E-02 -1.5E-02 -1.0E-02 -5.0E-03 0.0E+00 5.0E-03 1.0E-02 0 5 10 15 20 𝑞 [𝑟𝑎𝑑∕𝑠] 𝑡[𝑠] q (with SAS) q (without SAS) -2.5E-02 -2.0E-02 -1.5E-02 -1.0E-02 -5.0E-03 0.0E+00 0 5 10 15 20 δe [𝑟𝑎𝑑] 𝑡[𝑠] δe (with SAS) 53 -2.0E+00 0.0E+00 2.0E+00 4.0E+00 6.0E+00 8.0E+00 1.0E+01 1.2E+01 1.4E+01 1.6E+01 0 100 200 300 400 500 600 𝑤 [𝑚∕𝑠] 𝑡[𝑠] w (with SAS & PID) -1.0E-02 -5.0E-03 0.0E+00 5.0E-03 1.0E-02 1.5E-02 2.0E-02 2.5E-02 0 100 200 300 400 500 600 𝑞 [𝑟𝑎𝑑∕𝑠] 𝑡[𝑠] q (with SAS & PID) 0.0E+00 1.0E-02 2.0E-02 3.0E-02 4.0E-02 5.0E-02 6.0E-02 7.0E-02 8.0E-02 0 100 200 300 400 500 600 Δ𝜃 [𝑟𝑎𝑑] 𝑡[𝑠] Δθ (with SAS & PID) 54 Figure 8.15 Full system simulation results of short period offline stabilization augmentation system for unstable Boeing 747 with phugoid PID control. Response to commanded ( ) ( [ ] ). Large time simulation. The obtained results show how the stability augmentation system reaches and stabilization despite their natural unstable behavior, even when subject to the variations of the desired deflection requested by the PID controller. Within the first seconds of the simulation, if compared to previous results, the variations of and take longer to stabilize. This is a result of the PID controller continuously requesting a different deflection in order to stabilize , and the stability augmentation system trying to reach this commanded at its every change. The values of and are, however, contained within very short ranges which, considering the unstable behavior of the aircraft, is a success. Considering the long time behavior, the system reaches a stable state resulting from both the action of the PID controller and the stabilization augmentation system. The order of precision of the integration and differentiation methods used for the PID is much simpler than the ones used to develop the stability augmentation system, for this example only serves as a theoretical test run of a possible integration of two systems. The constants used for the PID control are for the integral part, for the proportional and for the derivative. -6.0E-02 -5.0E-02 -4.0E-02 -3.0E-02 -2.0E-02 -1.0E-02 0.0E+00 1.0E-02 0 100 200 300 400 500 600 δe [𝑟𝑎𝑑] 𝑡[𝑠] δe (with SAS & PID) 55 9 Viability of online critic-adaptive stability augmentation system This section aims to assess the viability of using the proposed solution in a criticaladaptive system, with online functionality. We approach here the part of the problem of most current research and which, as of today, poses a major problem. The objective tried to solve is the one of stabilizing the system as has been done in the previously presented sections, with the addition that the aircraft can change its behavior, ie, the dynamic equations can vary. 9.1 Online stability augmentation system proposal The suggested online system is based on the well-known critical-adaptive system structure, which features a neural network defining the policy to be followed, the plant, and a block, or set of blocks, which evaluate the performance of the first neural network and update it accordingly. Figure 9.1 Block diagram representation of a critic-adaptive control system Even though a system based on the previous structure will eventually reach good training results, it is always encouraged to perform a previous offline training using a simulation of the plant. The structure of such offline training method, which is an extension of the one shown in Figure 8.2, follows. A new neural network has been added to the training, for it will be needed later when running the online version. This network simulates the airplane response to a given and and must be trained using a constant time step. Since the output of this network may consist of two or four values, one traditional ANFIS network will not suffice. Instead two or four networks will be required, 56 or as an alternative a M-ANFIS network, which is a modification of the normal version which can output multiple results. Figure 9.2 Block diagram of offline training phase in online stability augmentation system Once the offline training has finished presented structure can be used as the performance evaluation block shown in Figure 9.1. The resulting system is presented in the following diagram. Figure 9.3 Block diagram of online stability augmentation system application It is important to use a neural network to simulate the plant, because the data to perform the optimal control algorithm need to cover the whole range of the state variables, while the plant only offers data of the instantaneous state variables. 9.2 Viability of the proposal The proposed system is unable to perform the online training within logical execution times. While training times during the offline phase are on the order of days, mainly because of the inclusion of new neural networks in the system; the online implementation of the system cannot perform fast enough training to fit a real-time variant system. We find ourselves with the problem that is described at beginning of the study, according to which the stabilization of a system of this kind, with the use of neural networks, is beyond the reach of the tools (computers and methodology ) that are 57 actually available;, and which is the reason of existence of several research groups in this area. This result therefore expected. In order to accelerate the process, a multicore implementation of the optimal control algorithm and the neural networks might be developed to run on multiple computers. However, the number of computing units would be constricted to the capacity of the airplane. Nevertheless, although this would divide the run time, it is usually not the solution to such problem; instead, the consideration of variations or different approaches might be considered, as well as, although already optimized, further optimizations should be conducted over the codes and algorithms in order to decrease its execution time 58 10 Conclusions and further work Within this project, an airplane stability augmentation system for short period mode has been developed, using neural networks, fuzzy logic and optimal control. Moreover, its functionality has been contrasted with initial perturbations of the state variables and desired input deflection of the elevator, both constant and variable within a PID control system. The results were initially obtained for a model of Boeing 747. Later the approach has been tested on an unstable system obtained from the previous model. The proposed system offers a good response in its offline version, but due to the high dimensionality of the problem and technical limitations the solution developed cannot train a system in online mode. The online system, whose structure is derived from the offline version, requires the use of at least three neural networks running in parallel, plus a dynamic programming algorithm also running to continuously update the neural network values as the system behavior changes. While the use of multiple computers could render the system more agile (each neural network uses two strands in its online operation and, although not developed in this project, the dynamic programming algorithm is parallelizable), this would reduce the execution time, but does not constitute a real solution to the exponential growth of the dimensionality of the problem. Although the network structure presents several variations with respect to the classic neural networks, which make it faster, it would be necessary a further optimization or the consideration of a different approach. During the development of this study the linearized equations of the dynamics of the aircraft in its longitudinal mode have been used. However, the code is developed in a way so that the equations governing the system can easily be replaced by other data sources, whether they are non-linearized equations or by data obtained directly from the airplane in flight. However, although it is possible to train the system exclusively from data extracted directly from the plant, it is always advisable to make a prior training through simulation, so as to ensure a clean sweep of all possible states. The tools used by the system have considerable freedom and portability. They are easily includable in other projects, and their characteristics are slightly dimensionally linked to the problem itself: both the structure of the programmed neural network and the numerical integration code allow the use of a variable number of inputs. The dynamic programming code that performs optimal control calculations is in turn easily modifiable 59 and with a few changes can, for example, work with problems where more than a control variable is taken into account. This means that, despite the proposed system is developed within the framework of aeronautical engineering and as a stability augmentation system for an aircraft, using very few modifications, the solution developed is also applicable to many different control systems for very diverse problems. This project, because of its characteristics of study as well as for representing an introductory work to some of the main tools of artificial intelligence, presents many possibilities for future actions for it to be expanded and complemented. From a generic point of view, the study of alternative methods or optimizations for the proposed solution looks very attractive. These optimizations can vary from improving the training algorithms of the networks, the study of a possible restructuration or redefinition of them, or the improvement of the optimal control algorithm. Although all these aspects have been designed to operate in a fast and quite optimal way, as observed in the results of offline simulation of the system, the high complexity representing the online problem makes the research in these areas is very interesting. The study of different alternatives to the Performance evaluation block (see Figure 9.1) is also a source of many possibilities for future efforts; it is not coincidence much of the actual research of the field is centered on this point. From a more specific point of view, the number of additional factors to consider which may be of interest in this study is quite large: the response to gusts, the introduction of a delay in the system that represents the time between the neural network decision and the action on the elevator, or simply the sudden fail of components which alter the behavior of the aircraft; all them, amongst others, are factors whose consideration within this project may yield interesting results and present future lines of work. The proposed solution has not been developed from a unique source of information, but is the result of the consideration of many possible solutions and ideas, the refutation of some and acceptance of others. It results from the search for alternatives to problems faced during its development and the lecture of various sources that have served as inspiration and from which many details have been extracted to configure the final proposed solution. This process has enabled the author to gain valuable knowledge in all the fields in which the study is based; particularly in areas of a special personal interest but never studied 60 by him such as fuzzy logic, neural networks and optimal control. It has also allowed improving the knowledge of certain concepts of airplane flight mechanics as is, for example, the nature of the response of an aircraft and its importance in the design of a control system. This project also serves for both the reader and the author as link or introduction to more complex studies such as those who are nowadays being conducted by leading research groups in this area. These studies, whether or related to optimal control systems or adaptive critical control systems, are essentially based on the principles outlined in this work. 61 11 Environmental effects This study has been completely developed using a computer, therefore its environmental impact is virtually null. The use of a stability augmentation system, as discussed in the justification for the project, enables the creation of aircraft whose center of gravity is closer to the pressure center or exceeds it, thus representing a decrease on the negative lift in the horizontal stabilizer, and consequently a reduction in fuel consumption and emissions. Although this project represents a study, its approach is intended for possible use in real aircraft. This would require lengthy periods of training and testing which would have an associated cost in the form of emissions. However, given the nature of the proposed system, learning could take place during commercial flights so its effect on the environment would not be accounted as an extra, but would fall within the one incurred in the scheduled flight. Finally, depending on the accuracy with which one wants to provide the system and the input variables to be monitored by neural networks, the computing power required could be quite large and, as a consequence, the energy required for its execution would be considerable. 62 12 Budget The budget associated with this project is broken down into the following: work hours, required software and used hardware. 12.1 Work hours The working time spent on this project has been about hours per week during approximately months. Considering an hourly rate of ⁄, 12.2 Used software All the software has been used under academic licensing, that is, at no cost. Nevertheless the licenses costs are presented  Microsoft Visual Studio license:  Microsoft Office Hogar y Empresas: 12.3 Used hardware This project has been developed on a years old personal computer with a current value of . The total project budget follows in the table below: Item Cost (academic licenseng) Cost (normal licensing) Work hours whork hours Software Microsoft Visual Studio Microsoft Office Hardware Personal computer Total Table 12.1 Budget breakdown