Sunday, November 18, 2018

Accelerometer Position and Normal Acceleration Augmentation System


Imagine a very simple conventional aircraft, with very simple mechanical controls. From an initial straight flight condition of constant altitude, constant airspeed and wings leveled, in a given moment, for some reason, it is required that the aircraft climbs with a constant angle of trajectory. In order to perform this maneuver, the pilot will have to find the right combination of surfaces deflections (let’s say elevator and horizontal stabilizer) and possibly engine throttle setting.  Now imagine that there was a device interfacing the pilot and the surface deflection commands such that, somehow, it was possible to the pilot to simply “communicate” to the flying machine his/her intention to climb with a constant angle of trajectory. This is part of the idea behind the Fly-By-Wire (FBW) systems, available in most of the modern aircraft.

Bill Palmer, in his book Airbus A330 Normal Law, makes a very good summary of the philosophy of modern Fly-by-Wire systems: to make possible to the pilot to tell the aircraft safely  “what to do” instead of “how to do it”. Usually this will be related to parameters of the aircraft performance. Using the previous example to illustrate, we would command an allowed angle of trajectory instead of surface deflections.

A conventional input command would deflect the elevators and this would eventually cause a variation in the overall lift on the aircraft (as explained in this other article); in an augmented system (FBW), usually the longitudinal input commands a pitch rate or a rate of angle of trajectory (called gamma). A possible approach to implement a controller for this purpose would involve feeding back the normal load factor to the elevator command. But, in order to improve the handling qualities of the augmented system, we should find out how to get over sensor signal issues like, possibly, the Non-Minimum-Phase (NMP) zero of the natural response of load factor when measured at the aircraft C.G., as described by the aforementioned article.

Is it possible to get rid of Non-Minimum-Phase (NMP) zero in load factor response? Is it related to the position where the load factor is measured? Is there a better position to place the accelerometer to feedback purposes? These are some questions that this article intends to answer or at least to address a path to the answer. It can be seen as a second part of my previous article about NMP zero in normal load factor response.

Load factor response


Let us try to understand this phenomenon qualitatively first. As introduced in a previous article, we call load factor, the acceleration of the aircraft in the normal stability axis measured in g (the acceleration of gravity) or simply the ratio of the lift to the weight, such that, a level flight is said to be an 1-g flight, for example, once that the lift is balanced by the aircraft weight in such condition. If this is a steady condition, i.e. straight flight with wings leveled and no pitch rate, regardless of the position where we put our sensor (accelerometer), we will obtain the readout of 1 g. We are considering different positions here as the various stations throughout the longitudinal axis of the aircraft as shown in Figure 1.

Figure 1. Examples of different possible stations for accelerometer positioning. The white circle represents the CG position and the reference such that x=0ft at it.
What if, all of sudden, we apply a longitudinal command such that the body develops an angular acceleration, that is to say a pitch rate derivative (keeping no yawing, nor rolling moments)? Now the readout of an accelerometer placed at a distance of xa from the aircraft CG could be approximated by the following expression:

nz_xa = nz_CG + xa*qp/g0

nz_CG: load factor measured at aircraft CG using positive sign convention upwards (straight level flight is +1g, for example)
xa: distance of accelerometer position to aircraft CG using positive sign convention towards the aircraft nose
qp: pitch rate derivative using positive sign convention when the nose is lifting up related to the horizon
g0: acceleration of gravity

Equation 1. Load factor measurement at different stations

An interesting fact about the Nz response when we apply an elevator deflection to pitch up the aircraft, for example, is that, depending on the accelerometer position, we will see firstly a decrease of Nz before building up (the actual expected response). As explained in the other article, this happens because the elevator deflection produces a decrease in the overall lift on the aircraft before generating the moment to increase the angle-of-attack of the wings which will increase the lift and finally will overcome the initial lift decrease.

As usual, we need a model

As in the previous article, we will use here the F-16 model from Stevens & Lewis that I implemented as a GitHub project called McFlight. For our purposes, a state-space representation of the linearized aircraft model around a given flight condition will be enough. Find below the description of such model:

  • Steady flight condition: straight level flight at sea level, speed V = 502 ft/s, cg at 0.35 of MAC (Mean Aerodynamic Chord), and mass m = 20490 lb
  • State (X):
    • u_ftps: speed at longitudinal axis (x)
    • w_ftps: speed at vertical axis (z)
    • theta_rad: pitch angle in radians
    • q_rps: pitch rate in radians per second
  • Controls (U):
    • throttle: engine command (from 0 to 1)
    • elev_deg: elevator deflection in degrees
  • Outputs (Y):
    • nx_g: acceleration in longitudinal axis in g
    • ny_g: acceleration in lateral axis in g
    • nz_g: acceleration in vertical axis in g
    • alpha_deg: angle-of-attack in degrees
    • q_rps: pitch rate in radians per second
    • Q_lbfpft2: dynamic pressure in lb-force per square feet
    • mach: Mach number
    • nzs_g: load factor in stability axis, i.e., 1g at straight level flight condition

The data in matrices A, B, C and D can be seen directly at GitHub repository, at the trim sce file, specifically.

Prior to any further investigation, we have to make sure that our linear model  is as representative as the non-linear model in the surrounding area of the steady flight condition. For this purpose, let us compare theta and q outputs when the system is input by an elevator step.

(a)

(b)

(c)

Figure 2. Comparison between linear and non-linear models for an elevator step input (a) and corresponding theta (b) and pitch-rate (c) responses.

As shown in Figure 2, for an elevator deflection of one degree, our linear model presents a very good matching with the non-linear model for maneuvers within 8 deg/s of pitch rate and 16 deg of theta, which are good enough for our studies.


Measuring the load factor



From now on, we will use the linear model. As we described above, this model gives us the load factor measurement at CG position (nzs_g). However, we can apply the equation 1 to calculate the load factor response at different positions throughout the aircraft longitudinal axis. In the state-space representation we can perform some matrices operations to produce the new output.



Let us use the suffix p to indicate a derivative, then we can write the state-space representation as:

Xp = A.X + B.u
Y = C.X + D.u

Algebraically we have to calculate some new lines at matrices C and D such that we obtain the new required output, which we will call nz_acc_g (load factor at accelerometer). Considering the order of states, controls and output variables as specified previously, we know that the fourth element of Xp will give us the derivative of pitch rate, i.e., Xp[4] = qp, that can also be written as:

Xp[4] = A[4,:] * X + B[4,:] * u

where A[4,:] represents the full fourth line of matrix A and B[4,:] the full fourth line of matrix B. We also know that the eighth element of Y will give us nzs_g, i.e.:

Y[8] = C[8,:] * X + D[8,:] * u

In order to calculate the load factor at the accelerometer position, using equation 1, we will have:

nz_acc_g = nzs_g + (xa/g)*qp
Y[9] = Y[8] + (xa/g)*Xp[4]
Y[9] = (C[8,:]*X + D[8,:]*u) + (xa/g) * (A[4,:]*X + B[4,:]*u)
Y[9] = (C[8,:] + (xa/g)*A[4,:])*X + (D[8,:] + (xa/g)*B[4,:])*u

In other words, we will augment the matrices C and D with the new lines from the calculations above. In Scilab code, for instance, as can be seen at GitHub, we will have:

C_body(9,:) = C_body(8,:) + l_arm_ft/params.g0_ftps2*A_body(4,:);
D_body(9,:) = D_body(8,:) + l_arm_ft/params.g0_ftps2*B_body(4,:);

where l_arm_ft = xa and params.g0_ftps2 is the acceleration of gravity.

The Figure 3 shows the load factor response for three different values of xa positions. Note that as xa increases, towards the aircraft nose, the initial nz response in the opposite direction (negative) decreases such that when xa is 6.1ft, the nz response is always greater than zero throughout the time. In the frequency domain, we will observe that the zero at the right-half of s-plane for measurements at CG will move out towards the infinity as xa increases. Then, at xa=6.1ft, the NMP zero effect disappears or at least is not noticeable (red curve at Figure 3). The point where this happens is considered an “instantaneous center of rotation”.

Figure 3. Load factor measured at different stations of the aircraft.
For the sake of accuracy, we didn’t produce the xa=6.1ft by trial-and-error. McRuer tells us that the “instantaneous center of rotation” can be found by the quotient of two control stability derivatives, i.e.,

xa = Z_delev/M_delev

Z_delev is the dimensional derivative of the force in Z axis of body reference due to elevator deflection
M_delev is the dimensional derivative of pitching moment due to elevator deflection

These derivatives are calculated by the script stability_deriv_body at McFlight.

Accelerometers position


As you may expect, I paved all the way down here to say that the best place to install the accelerometers is in the “instantaneous center of rotation”. Indeed, but as noticed by McRuer, this location varies as CG position and effective control arm change depending on the aircraft configuration. One possible alternative would be choosing somewhere in the neighbourhood and/or estimate its position taking into account the accuracy of the other parameters involved in the calculation (as pitch rate derivative, for example).

In case of fighters (as our F-16), in order to provide the pilot with a very precise control of g-load in high-g maneuvers, the nz control law can use the signal provided by an accelerometer placed at the pilot station (Stevens & Lewis, 2003).

In summary, the best choice to place the accelerometers will take into account several requirements, including the aeroservoelasticity ones to avoid unintended interactions between the flight control law and the structural modes of the aircraft. This leads, for example, to place the sensor close to a node of the fuselage bending moment. Anyway, once again, aircraft development is not achievable by the decisions of one single discipline.

Wednesday, May 30, 2018

Six DoF Nonlinear Aircraft Model: Challenges in Software Development


Recently I wrote an article where I needed an aircraft model to simulate the phenomenon I was trying to describe. I’ve used a six-degrees-of-freedom (6 DoF) nonlinear model of an F-16 provided by the book Aircraft Control and Simulation (Stevens & Lewis, 2003). The codes were in Fortran and I translated them to Scilab scripts doing a little bit of software engineering and turn all of them into a small project at GitHub named McFlight. Through the next lines I will detail what is supposed to be a six-degrees-of-freedom nonlinear aircraft model, using the structure at McFlight as an example. And also I will try to assess how far (or close) such model is from the ones used by aerospace industry.

What is a 6 DoF Nonlinear Aircraft Model?


A very good approach to make science is to find a mathematical model to represent the piece of world you are interested in. Once your model is validated, you are able to make good predictions on that specific field. However, the beauty of a model is that you focus your spyglass on the relevant details of the observed phenomenon. For example, if I am studying the aerodynamics of wings, I would certainly be interested to model the flow detachment developing throughout the wing surface. But in a flight mechanics model I just need to know how aerodynamic forces change due to variables like angle-of-attack, Mach number, angle-of-sideslip, control surfaces deflections, regardless of how the flow detachment or any other specific aerodynamic phenomenon happens.

So we have an aircraft model to study the flight mechanics. This model is a set of ordinary differential equations and we may represent it as shown in Figure 1.

Figure 1. Differential equation for dynamic system model.

Basically we have a variable called state (X), a vector with the variables which describe the dynamics of our system, in this case, for example, we would have angle-of-attack, airspeed, pitch angle, rolling angle, sideslip angle, etc. And a variable called control (U) which represents the inputs in our system, in this case, the deflection of the aerodynamic surfaces (elevator, rudder, aileron, flaps, etc), engine setting or any other control input. Our model is the solution of the differential equations system shown in Figure 1. It is important to note that the input variables must be independent, otherwise it will be a state or output variable. The identification of the inputs is not always straightforward such that we have to be careful in order to avoid ending up with an incorrectly determined system of equations.

In mathematics, we say that a system is nonlinear when the output is not necessarily proportional to the input. For example, we cannot say that when we deflect the elevators, the lift force will increase by the same amount, differing only by a multiplication factor or that if we deflect the elevators in 5°, the increase in the lift force will correspond to the ones equivalent to 3° plus 2° deflections. In other words, the system will not necessarily respect the mathematical properties of linearity, which are:

  • Additivity: f(a+b) = f(a) + f(b)
  • Homogeneity: f(a*x) = a*f(x)

Nonlinear systems are much more representative of physical phenomena than linear systems. However, we have strong mathematical tools to work out on linear systems, including some of the most used design techniques for automatic control. Fortunately, we can link these two modeling approaches finding a linear approximation of a nonlinear model in the surroundings of an equilibrium point. Then we can make use of the linear mathematical tools, always paying attention to the assumptions taken on linearization process.

Finally, when we say six-degrees-of-freedom (or 6 DoF) we mean that our model is considering 3 translational and 3 rotational directions of movement as shown in the Figure 2 borrowed from this other article in which I discuss the “Aircraft diving before climbing or NMP zero in load factor response”.

Figure 2. Six degrees of freedom. The straight lines mean translational movements of back and forth following the arrows directions; the arches mean rotational movements clockwise and counter-clockwise around the center of gravity of the aircraft (assumed to be in the circle with the cross in the pictures). (Embraer super-tucano three views were found here).

Software requirements for the implementation


The main objective of the model we are describing here is to assess the flight mechanics and eventually to support the design of automatic flight controls. Then let us follow briefly a typical cycle for flight controls system (FCS) development so that we can tackle the main needs from a software engineering perspective.

Considering a classical approach, as stated by Rauw, the very first stage would be running (guess what) a 6 DoF nonlinear aircraft model starting from different trim conditions (equilibrium points, as explained below), mainly to get the feeling about the plant dynamics and aircraft mechanics as a whole. Also this nonlinear model will be the input for the linearization process which, in turn, will provide the linearized models for the controller design. Notice that both linearization and simulation operations will require a trim condition (provided by a trimming operation), but they are not related to each other. In other words, the linearization of the nonlinear model does not depend on the simulation of this model.

Breaking down the operations in further details:

  • We need to trim the aircraft model. This means we have to find the inputs and states values to “balance” the plant, such that the states derivatives are equal to zero. Finding these equilibrium conditions can be done by building a cost function whose minimal value gives us a trim condition as a solution (this is a possible approach and the one we are using here);
  • We need to run the model. This means solving numerically the differential equation representing the model (Figure 1), using the trim conditions found above to initialize the model;
  • We need to linearize the model for different trim conditions. This means finding a first order approximation for the model function in the surrounding of a given equilibrium point. This will support the controller design (assuming we are using a linear design approach).

Trimming and running the model enable us to simulate (off-line, which means here that is not real time simulation as Rauw called) many conditions and maneuvers throughout the flight envelope, helping even to clarify the requirements for the controller itself. Also we will need to assess how to discretize the flight envelope for the controller design, which means choosing the points to linearize the model. But what we are interested to highlight here is that the activities listed above (mainly trimming and running), which use the nonlinear aircraft model, will be executed very often, so it is desirable that they can perform computationally as fast as possible. We don’t want to wait minutes to trim the model or that 1 second of off-line simulations run in 1 or more real second.

Once we finished the linear design of the controller (using linear analysis with robustness, comfort and some performance criteria) we have to keep going with the development adding the controller in the full nonlinear aircraft model “to ensure that the design works in this domain” (as shown in this GARTEUR book). Notice that we will have to take our designed linear controller to the nonlinear (more representative) full model, doing the necessary adaptations. And as taught by master McRuer, the next step “usually involves a series of increasingly realistic tests that attempt to thoroughly analyze the system in both normal and abnormal operating conditions”. This may imply the use of different test environments: from the real FCS software plugged into the nonlinear aircraft model to the “Iron bird” tests. We have to assure the equivalence between each test/development environment such that, in a final extent, we can be confident that the FCS eventually flying is the same designed on ground.

The development of an aircraft model is the result of a team work. We have propulsion data, aerodynamic data, mass distribution data, eventually we have to be aware about the sensors, the actuators and the dynamic modes of the structure, and any further data from other disciplines involved in the process. During the development these data will be iteratively updated. We have to be able to isolate these categories of data in our aircraft model such that we can quickly update each part separately and possibly add new parts.

Now hopefully we have a more clear picture of the main concerns when implementing our model:

  • Performance. A software that runs quickly is always very welcomed. Throughout our design cycle we will need to run and re-run the model many times both in off-line or real time simulations, possibly in many different environments. But we have to be careful in how much performance we really need. We may waste too much time and efforts trying to improve code performance beyond our needs, risking to end up spoiling some other important code properties, such as readability. Although it is very difficult to figure out how much performance we will need in early stages of development, we can modularize the model so that we can improve the performance bottlenecks in later stages if it is necessary. Also, we can use different models in different environments (not all of them requires the same level of performance and mutability), but in this case we must have strong procedures to assure the equivalence between the models or, for example, trusted qualified tools to translate from one model to the others;
  • Modularization. In software engineering history, we have an extensive list of reasons why you should modularize your code. In this specific case we can point out few but very strong justifications for doing so. We will have many people working in the same model at the same time, merging their work can be easier if they are changing separated modules. We will have parts of our model coming from other teams (possibly using different software tools) and we will also have to provide parts of our model to other teams; again, if we have well defined modules, this will be more efficiently done. And finally, as we stated in the previous item, if, for example, we find necessary to improve performance, this will be easier to change if we tackle specific modules of our model. Notice that most of the advantages of modularization comes from a very well defined interface between the modules, this means an interface that is very unlikely to change. For example, no matter how the aerodynamicists change their data and calculations internally, we will have our integrated model still working as long as their module keeps receiving alpha, beta, Mach, etc, as inputs, and providing the aerodynamic forces and coefficients as outputs;
  • Readability. Once again, an aircraft model is a work of many hands. When caring about the readability, we reduce the amount of time required for the understanding of a given piece of code, even for the author of that code when looking at it again in the future. This also reduces the possibility of coding errors, making changes faster and safer. We talked so far about coding, but it is important to say that the widespread approach in industry for dynamic systems modeling is the model-based design which relies on its visual method to improve readability.

Implementation: McFlight architecture


As we explained previously, our McFlight is a possible implementation of the model provided by Stevens and Lewis. For further details about their model, please refer to the section 3.5 (Aircraft models for simulation) of the second edition (2003). Details about the equations will be found in the chapter 2 of this book. The implementation of the 6 DoF nonlinear F16 model is given by the authors using structured programming in Fortran, according to the following set of functions:
  • Subroutine F: is the actual dynamic model containing the equations of motion, equivalent to the mathematical function f at Figure 1. This subroutine makes use of all the other functions below;
  • Subroutine ADC: meaning airdata computer and responsible for the atmospheric calculations (mach, dynamic pressure, etc);
  • Functions TGEAR, PDOT, RTAU and THRUST: responsible for the engine calculations;
  • Functions for aerodynamic data: a set of functions representing the aerodynamic coefficients and forces.
  • A cost function for trimming.
It is also provided some numerical algorithms for trimming, differential equation solving and linearizing.

Our (re)implementation of this model will try to pursuit the requirements described in the previous section. The high-level description is illustrated by the UML (Unified Modeling Language) diagram of components shown in Figure 3. For those not familiar with UML, this diagram aims to focus on the interfaces between the different parts of the software, each of these parts with a very well-defined concern (a component). The components are represented by a square with two small other squares on the left. The dotted-line arrow indicates a relation of dependency between the components, meaning that the origin depends on the destination. We also enclosed some components with similar concerns in what UML calls package (the square with an upper small square with a label). Notice that we are able to replace a component for another one with a different internal implementation but providing the same services to their dependent components, and keeping the whole system still working. In other words, as long as we keep the contract assumed between two components (the interface), they are transparently interchangeable by other ones implementing and using the same interfaces.

Figure 3. UML component diagram of McFlight.
We used Scilab scripts for the implementation and despite of not coding a canonical description of component (how it would be in Java, for example) our modules try to follow the same concept. In the following items we link the diagram at Figure 3 to the code at GitHub.
  • Atmosphere (atmosphere/atmosphere.sci): model of the International Standard Atmosphere, providing temperature, pressure and density data based on altitude;
  • Aerodata (eqm/aerodata_f16.sci): functions providing aerodynamic forces and coefficients data based on aircraft angles and control surfaces deflections;
  • Engine (eqm/engine_f16.sci): functions providing engine thrust and dynamics based on throttle setting;
  • Params (eqm/params_f16.sci): function returning geometric and mass properties of the aircraft;
  • EQM (eqm/eqm.sci): actual implementation of the equations of motion of the aircraft, i.e., the function f represented in Figure 1. It requires the implementation of all components enclosed in the aircraft and nature packages shown in Figure 3.
  • Trimmer (trim/trim_f16.sci): algorithms and functions for trimming the dynamic system represented by EQM;
  • Simulator (sim/sim_f16.sci): algorithms and functions to solve the dynamic system represented by EQM, i.e., to run the flight simulation;
  • Linearizer (sim/lin_f16.sci): algorithms and functions to find a linear model in a given trim condition.
I am not an advanced user of Scilab (yet) but I used two approaches to decouple the dependency between the modules (note: by “decoupling” here, we mean that the dependency relies on the well-defined interfaces as much as possible): (1) calling different .sci files which implement the functions with the same signature expected by the script (for example, in order to change the engine, wherever a file is executing engine_f16.sci on top, I could replace by, let’s say, engine_delorean.sci); or (2) passing a struct as a function parameter, assuring that it follows the expected data structure (for example, the parameter params in eqm.sci). If we succeeded in decoupling the modules implementations, we would be able to do things like, for example:
  • Replace the module Aerodata for a newer version of aerodynamic data;
  • Change the internal algorithms used by Trimmer, Simulator or Linearizer;
  • Few changes enhancing the model to include sensors, actuators and a feedback controller;
  • Replace the whole package Nature for other planet atmospheric model (and also gravity model which is currently tangled inside eqm and then not enough decoupled) enabling us to simulate our aircraft somewhere else in the universe;
As we said before, an alternative approach to implement dynamic models (and actually the de facto standard in aerospace industry) is the model-based design. We just want to make clear that it is possible to use this approach still keeping the decoupling philosophy behind the diagram at Figure 3. Indeed we are starting the equivalent implementation using Scilab Xcos.

To finish this section I would like to share a strategy in the methodology I followed when writing the scripts. For each module I developed, I tried to build a piece of automatic test. For example, for the atmospheric model, I wrote a CSV file with some data directly extracted from the NASA document of the standard atmosphere, and made a script to automatically compare each line of this CSV with the values returned by the implemented function (check the folder called tests). This reminds me what I heard from a senior engineer once: prior to any assigned task, figure out how you are going to validate it.

How far are we from the real world of the industry?


Imagine that you have decided to spend your fortune of billions of dollars developing an aircraft by your own or at least the FCS of an existing one. Would it be possible of doing so with the very simple application described above? Well, certainly our model would have to be enhanced. Let us infer what else we need when it comes to the real world of the industry.
  • The data available would be massive. For example, the aerodynamic data would include the contribution of each surface in each different configuration of the aircraft, possibly in ice condition, throughout the whole flight envelope. One possible way to implement the nonlinear effects is using lookup tables and this certainly would give us thousand of tables to be interpolated. This means that we would have at least to include many other inputs to our aerodata module and be able to handle quickly the lookups in the tables;
  • Care about the architecture trade-offs wisely. It is very beautiful to have a decoupled and componentized architecture but do not overengineer it. For example why would we isolate the codes related to atmosphere and gravity to the extreme if our aircraft surely will never fly in outer space? This can even compromise the software performance.
  • Care about the software performance and ways to improve it when it is necessary. Being very specific, we shall always have means to compile our model or at least parts of it to make it faster, for example. As I said before, remember that we have to linearize around several points of the flight envelope and simulate before that, to check in what regions we should increase the number of points (transonic region, for instance) to improve our design.
  • Eventually we will have to mix discrete and continuous time simulations. This means we would have to face all sort of problems with sampling, for example. The model-based design approach makes it easier.
  • Have a robust version and configuration control of each module, or artifact. Eventually we will have to reproduce some simulation environments with different versions of engine, aerodata, mass properties, etc, for example. Having a rigorous tracking framework for our artifacts changes, helps us avoiding and fixing problems. Additionally this would make our lives easier when we need to talk to the aeronautics certification authorities.
Bad news (or good news) are that you won’t be able to make it by yourself. The pieces of information in each of the blocks in Figure 3 would be the result of many hours of work of many engineers. In the real world, the cherry on top of our challenge here would be the management effort to keep all the teams working synergically to achieve what was planned, constrained by the predicted time and budget. Piece of cake!

References


Magni, J.F., Bennani, S., Terlouw, J. (Eds.): Robust Flight Control: A Design
Challenge. Lecture Notes in Control and Information Services 224, Springer Ver-
lag, 1997.

McRuer, D., Ashkenas, I., Graham, D.: Aircraft Dynamics and Automatic Con-
trol. Princeton University Press, Princeton, New Jersey, USA, 1973.

Rauw, M. FDC 1.4 - A Simulink Toolbox for Flight Dynamics and Control Analysis.

Stevens, B.L., Lewis, F.L.: Aircraft Control and Simulation. John Wiley & Sons
Inc., 2003.

Wednesday, April 18, 2018

Aircraft diving before climbing or NMP zero in load factor response


Source: pixabay

Imagine a conventional aircraft, like the ones which are flying across the world carrying passengers, for example. There is a very subtle and intriguing fact (at least for flight mechanics geeks) when the aircraft is commanded “to climb up”: the very first response of the aircraft body, for a short time period, is in the other way round, i.e., it dives a very little bit before climbing up. This fact is well explained by the aeronautical engineering experts and I’ll try to build up this explanation throughout the following lines. I’m assuming the reader has only some Mechanics background in the beginning, but the last section will require some Control Theory background, once that we will try to illustrate this puzzling concept called Non-Minimum-Phase (NMP) zero.


What actually happens


Let us start picking up some facts very well explained by the giants (Stevens And Lewis, John Anderson, Robert Nelson, etc). A conventional aircraft, flying straight with no changes in its altitude, has wings with this magical shape such that, roughly speaking, the greater the angle between a centerline through the wing section and the wind direction (called angle-of-attack), the greater the magnitude of the supernatural force upwards, called lift (Figure 1). The basics of the aircraft longitudinal commands (upwards and downwards) will be provided by the balance of the forces generated by the main wing and the horizontal stabilizer (little wing at the rear of the fuselage). More specifically, the longitudinal dynamics of the aircraft will be performed by the elevators, the small moving surfaces at the rear of the horizontal stabilizer. Moving up or down (both left and right simultaneously), the elevators are able to vary the vertical forces on the horizontal stabilizer making the aircraft climb or dive in a short term (assuming constant the other aircraft commands inputs) (Figure 2). 

Figure 1. Angle-of-attack and lift force, assuming straight level flight.


Figure 2. Elevator commands and longitudinal dynamics. Left sketch: elevators up, force in horizontal stabilizer downwards, counter-clockwise momentum, nose-up movement; right sketch: elevators down, force in horizontal stabilizer upwards, clockwise momentum, nose-down movement.
Considering our aircraft, flying straight and steadily with constant speed and altitude, if the elevators are moved upwards, we have the following sequence of events:

  1. The final elevators position causes a decrease in the lift force on the horizontal stabilizer. The force vector can even possibly turn to downward sense;
  2. The overall lift force acting in the aircraft decreases and the weight now is greater than the lift force, making the whole body accelerate downwards;
  3. The variation of the force in the tail generates an angular momentum, resulting in a pitching up movement (aircraft nose up);
  4. The nose up movement increases the angle-of-attack of the wing and therefore the lift on the wing;
  5. The angular momentum variation due to the new angle-of-attack of the wing counteracts the momentum variation due to the vertical force on the stabilizer, eventually breaking the pitching up movement;
  6. The increase in the lift on the wing is greater than the decrease in the lift on the horizontal stabilizer such that the overall lift force results greater than the weight, making the whole body finally accelerate upwards.

This sequence of events is assuming a time frame of few seconds after the elevators command, because in a long term, the aircraft will decelerate and require other input commands to restore the equilibrium condition or to settle in a new one.

But notice that the very first dynamic response of the aircraft (described by item 2) is in the opposite sense of the intended response, which is increasing the overall lift force of the aircraft. And so far we have a hypothesis. We need a model to confirm our suspicion and to dig into the details a bit deeper.

We need a model


Making science is not an easy job, neither is it cheap. So let’s trust the giants Stevens and Lewis and borrow their F-16 model. Just considering the hours of wind tunnel tests, the hours of engineering work compiling the aerodynamic data and the hours of real flight to validate the data make the costs go beyond the hundreds of thousands of US dollars (the hourly operating costs of the F-16 is around US$ 8,000). Luckily to me it only costs few days translating the Fortran code to Scilab and typing the aerodynamic data to the computer. I made the code available at GitHub, named mcflight.

Our model is a dynamic representation of an aircraft able to move in so-called 6 degrees of freedom (shortly 6 DOF): roughly speaking, 3 translational and 3 rotational movement possibilities (Figure 3). Although, for the simulation we are interested in, we only need 3 degrees of freedom, once that we are dealing only with longitudinal axis (i.e., translational movements upward/downward and backward/forward, and rotational movement of aircraft nose up/down).


Figure 3. Six degrees of freedom. The straight lines mean translational movements of back and forth following the arrows directions; the arches mean rotational movements clockwise and counter-clockwise around the center of gravity of the aircraft (assumed to be in the circle with the cross in the pictures). (Embraer super-tucano three views were found here).
By dynamic representation we mean a set of differential equations whose solution will provide us with the description of the aircraft movement throughout the time (for a very detailed explanation about the equations please check chapter 2 of Stevens and Lewis). In order to validate our hypothesis, we will need to “make” the aircraft fly steadily (wings level, straight flight) and then apply a small deflection on elevators to make the nose up. If we are right, we will observe a small decrease in lift force before building up.

To observe the lift force, in flight mechanics we use a physical quantity called load factor which is simply the ratio between the Lift force and the Weight (Figure 4). Although it is dimensionless, we express the load factor in g (referring to acceleration of gravity) and it can be positive, when Lift is in upward sense, or negative, otherwise. A very good g-meter is our stomach and you probably have experienced some g’s (positive) when riding a roller coaster. When an aircraft is flying straight and you are indulging yourself with a hot beverage served by the flight attendant, the load factor is 1 g (i.e. weight equals to lift) at that moment. During a normal airline flight you will experience slight variations around 1g, but the load factor will always be greater than 0. The negative g’s are experienced in aerobatic flights, for example, and notice that it is not only a simply free fall (which would be 0 g) but an acceleration downwards.

Figure 4. Load factor. Ratio between Lift Force and Weight.

Finally simulating


We are going to run a set of scilab scripts, analyze the time history of some variables (such as load factor) and use a lot of imagination to see an F-16 flying around.

The first step is to put our aircraft in a steady flight condition, which in our case means finding the elevator deflection and throttle setting (i.e., thrust provided by the engines) to allow a flight at constant speed, constant altitude, constant angle-of-attack. In other words, this is an equilibrium condition. The translation of this problem to our physical model is finding the solution of the differential equation using the proper set of constraints. In Aeronautics this is called trimming the aircraft. So let’s trim the aircraft in straight flight with wings level at sea level (altitude 0 ft) and with a speed of 200 ft/s. Asking scilab to solve it for us, we find the condition shown in Figure 5.

Altitude (h): 0 ft (sea level)
Speed (V): 200 ft/s
Throttle setting (engine): 0.287 (from 0 to 1)
Elevator deflection: 0.72° (-25.00° to 25.00°)
Angle-of-attack (alpha): 19.70º

Figure 5. Trim condition.

Let’s then apply the elevator deflection. In a time history simulation, we can say that it is a good practice to leave the aircraft fly steadily for a fraction of second (at least) just to check if there is an equilibrium point, before applying the input change. The following curves will show a time frame of 3 seconds of flight simulation with an elevator deflection of -1°, applied from 0.50s to 0.53s.

Figure 6. Load factor (Nz) and elevator deflection curves.
As we expected, the first dynamic response of the elevator deflection is a decrease in Nz. Then (after 0.53s) there is a monotonic increase in Nz until the end of the simulation. The Figure 7 shows the correspondence of the explained items in the previous section and the moments in the simulation zoomed in the first second of simulation.

Figure 7. Nz response right after elevator deflection. The numbers in parentheses indicate the corresponding item of explanation in the first section.
Notice that the amount of variation of Nz is around 0.004g, decreasing from 1.000g to 0.996g and increasing back in less 0.2s. This would be probable unperceivable by someone in the aircraft. A fighter like the F-16 is designed to be agile, to have a very fast dynamic response. If we had a model of a bigger guy (a cargo aircraft, for instance), the variation of Nz certainly would be slower and with a bigger magnitude (hopefully we will have another aircraft model available at mcflight soon). Figure 8 shows an illustrative animation of the F-16 flight.


Figure 8. Simulation F-16 model. Notice that the animation illustrates the drop in altitude, which is not coincident with the drop in Nz in the events timeline. Also the drop in altitude was (very) exaggerated for matters of illustration.

We could play with some parameters of our aircraft model to artificially make it have a behaviour similar to the bigger aircrafts. For example, we could change its ability to spin around the axis of the wings, i.e., the ability to nose up or down. Physically this means changing the moment of inertia around this axis. Try this out with mcflight code, varying the value of AYY.

From a Control Theory Point of View


From now on, it will be necessary some Control Theory background to fully understand it. A brave reader can keep going without this.

Considering the trim condition described in the previous section, we performed a linearization in order to find the state-space representation (Figure 9). In our Scilab scripts, this means using the command lin to find the matrices A, B, C and D and syslin to build a Scilab state-space representation. Mathematically we are reducing our flight dynamics to a first-order linear system representation of the equations. This will allow us to use the very well established theory of linear systems but the drawback is that this simplification is only compatible with the original (non-linear) model in the surroundings of the trim condition. So, before doing anything, we have to check the matching between our linear and non-linear models.

Figure 9. State-space representation, where x is the state variables vector, u is the control variables vector and y is the output variables vector.

Our linear representation is assuming the following variables:

  • Control variables (u): throttle setting (engines) and elevator deflection;
  • State variables (x): air-speed (assume that this equals to aircraft speed), angle-of-attack (alpha), attitude angle (theta, angle related to the horizon, initially equals to alpha), altitude, engine power setting.
  • Output variables (y): load factor (Nz). 
We can now simulate our linear model (in Scilab we can use csim command) with the same input we used in the previous section, using as input only elevator deflection and watching only the Nz as output. Figure 10 shows an overplot of linear and non-linear responses of Nz. Notice that after 1.4s our linear model starts drifting away from the non-linear model. However, the time frame from 0.0s to 1.4s is enough for our purposes.

Figure 10. Comparison between linear and non-linear simulations.

We are now ready to take a look on the poles and zeros of the transfer function of the system elevator-to-nz (Figure 10). Notice that we have a zero in the right-half-plane, also known as non-minimum-phase (NMP) zero (actually, we have three positive zeros, but the ones nearby the origin are due to the fact we inserted altitude and engine power as states, you can check this, eliminating them as state variables). And this guy is the guilty of the phenomenon we are describing here as stated by Stevens and Lewis (again): “these types of zeros occur when there are two or more different paths to the system output, or two or more different physical mechanisms, producing competing output components”. Yes, we have horizontal stabilizer and wings producing lift and affecting the load factor.

Figure 11. Poles and zeros of linear system elevator-to-nz.

Well, next time you are in a straight level flight of an airliner and all of a sudden you feel a quick and slight “floating” out of your seat (Nz less than 1g) before “sticking” to it (Nz greater than 1g) when the aircraft starts climbing, you can tell the passenger at your side that such sequence of events was caused by an NMP zero in the plant. Actually this is very unlikely, unless you have a g-meter with you. Anyway, at least the explanation you know ;-)

Why should I learn programming?


During every aircraft development, in order to assure that the flying machine behaves exactly as it was designed for, the flight test is performed by the manufacturer, also to prove to the competent authorities that the aircraft is able to fly safely. For the flight test purposes, the aircraft is equipped with all sorted of sensors (let’s say around 100,000 parameters), whose measured data we could assume recorded each 10 milisseconds. Doing a rough calculation we conclude that a single hour of flight test produces around 270 Gigabytes of data (in real life is actually much more than that).

The main point here is that for a single flight, there is a huge amount of data available to the aeronautical engineers validate and improve their design. But this only can happen in a phase of data analysis, which comes after few steps of processing the raw data, obviously unthinkable without some programming work. In a worst case scenario, the aeronautical engineers with their skills provided by a semester-long course of structured programming will struggle to write a little script from scratch to process the data and to allow them to check (for the thousandth time) if the aircraft behaved as expected.

This is not an isolated case of aeronautical engineering. Nowadays, in any field of science or engineering, our capacity to generate data has increased dramatically. For instance, the European Organization of Nuclear Research (CERN) announced (already in 2013) that their data centre had 100 petabytes of stored data, “equivalent 700 years of full HD-quality movies”. And the amount of general data produced daily is even more impressive. In every minute YouTube users share 400 hours of new videos, Facebook users share 216,302 photos, Instagram users like 2,430,555 posts.

This tells us that a problem-solving in this world overwhelmed by data is more and more walking towards formulate the problem for the computer helping us to solve it, no matter the field of knowledge we are talking about. This requires knowing how to ask the machine to solve our problem. As stated by an important Brazilian researcher, “every science is a computer science”.


Computer Science gap in other Sciences/Engineerings


If we compare a programming course in a standard engineering graduation today and 30 years ago, we would be impressed how topics are pretty much the same, sometimes only replacing Fortran by C, Python or even Matlab. Of course the basics of algorithms and problem-solving is very important to be kept, but the new engineers should at least be aware of the newest capabilities of programming and even the tools and techniques supporting the programming activity.

Resuming our example in the beginning of this text, aeronautical engineers 20 years ago could handle well (or, let’s say, using cutting edge technology by that time) the data available with their piece of Fortran code. Today, after getting (with a couple of clicks) the hundreds of gigabytes of data to analyze, they write a piece of code not that different than the one 20 years ago (in a different language, hopefully), but that runs in a hardware brutally faster and eventually achieve the expected results. Certainly, their job could be done faster by improving the code or by some different problem-solving technique; or perhaps many small repetitive tasks in the middle could be automated.

We could summarize this scenario as follows: in one hand we have engineers that eventually can make their job done, but is not fully aware that it could be improved by the latest programming capabilities available; on the other hand, we have computer engineers able to apply the best techniques but not required to help the other (non-computer) engineers.

Here we plotted the issue in the engineering field that is certainly far from being the worst case scenario, once that, at least, they have a programming course as part of their graduation. If we think about Human and Health sciences, the computer science gap may be much bigger.


Coding at School


Some people believe that in order to get over the aforementioned gap, we should go further and teach programming and related computer science topics earlier, at basic school. Hadi Partovi, the CEO of Code.org (a non-profit organization that aims expanding the access to computer science) said that “we don’t teach biology or chemistry to kids because they’re going to become surgeons or chemists, we teach them about photosynthesis and that water is H2O, or how lightbulbs work, just to understand the world around us. You don’t use any of it, but you do on a day-to-day basis use public-key encryption, and the average American has absolutely no idea what that is” (full text).

Besides initiatives as Code.org, programming contests worldwide has been becoming popular among youngsters before grad schools. The Brazilian Computer Society, for example, since 2002, offers modalities even for students under high school. The students shall solve computing and logical reasoning problems using just pencil and paper and, according to the Society, the goal is “to awaken the interest in computing problems and discover potential talents to programming”. The students with the best results in Brazilian Olympiad in Informatics will represent the country in the International Olympiad in Informatics which started at 1989 with 46 contestants from 13 countries and in its latest edition (2016) had 308 contestants from 80 countries.

Many countries have already included programming or “computational thinking” in their elementary school curricula. At European Union, 15 countries made it. And, at least in Brazil, where programming is not yet at the young students curriculum, private schools for “programming and robotics for teens” are rising in great cities like São Paulo, picturing a similar phenomenon which happened many years ago with the english private schools created all over the country. Sounds like the world is requiring the kids to learn a new language.


What about the future?


Studies published at the World Economic Forum Annual Meeting last year (2016) pointed some interesting facts about the jobs of the future: "The study shows that workers who successfully combine mathematical and interpersonal skills in the knowledge-based economies of the future should find many rewarding and lucrative opportunities."

If you realise that you can describe the routine of your day-to-day job in repetitive and well defined steps to execute some task, trust me, a computer program can do it way better than you can. This may sound alarming, but I can assure you that a significant percentage of a standard day at work of even very specialized science/engineering fields are covered by very repetitive tasks.

If at this point you are not convinced that you should improve your computational thinking skills, I quote the disturbing statement of this Brazilian researcher: “at this very moment, someone in the world is programming your job”.

Pragmatism doesn't mean reckless trial-and-error


When I finished my studies at university, I was eager to apply my technical knowledge in a real life context. And I confess that the very first contact with the crazy pace of job market was a bit shocking. At the time, as a purist, I incorrectly thought that using a pragmatic approach to solve a problem was akin to being reckless or irresponsible, which could only result in providing a sub-optimal solution. Later with experience it was eventually discovered that by wisely discarding parts of my theoretical knowledge it was possible to arrive more quickly at simpler solutions, without necessarily losing robustness.

According to Oxford dictionary, a purist is “a person who insists on absolute adherence to traditional rules or structures”. Also, it defines pragmatist as “a person who is guided more by practical considerations than by ideals”. This text intends to show that a purist approach can be relaxed and “more practically guided” still delivering successful solutions. This will be illustrated by a few examples in software and engineering contexts, but the technical details will be presented as lightly as possible so the reader from other fields can apply the examples to their own situations.

Would Physics equations from 1687 to 1905 be wrong?


In 1905, Albert Einstein proposed the ideas of his Special Theory of Relativity by the paper “On the Electrodynamics of Moving Bodies”, which pointed out some inconsistencies in Newtonian mechanics (1687) when speed approaches a significant fraction of the speed of light. Again, when the speed of the moving body is a significant fraction of the speed of light. Aircraft, for example, were (and still are) built using pure Newtonian physics with a successful history that speaks by itself. And how come, if these equations are not with the best precision? Well, Special Theory of Relativity adds a precision that is not required whatsoever by the models used for flying machines design. To give us an idea, the fastest speed ever attained by humans (11,1 km/s by Apollo 10 astronauts) will make Newtonian equations differ from Einstein equations by a factor of 0.9999999993 (being 1.0 corresponding to an identical result from both theories).

This bit of history helps us to illustrate how science works or more specifically how the scientific method works. When observing the problem we are interested to solve, we develop mathematical models, which are approximations of the natural phenomenon we are studying. The quality of a model in guiding us towards a solution will be determined by the degree of precision required in our solution. In other words, if our moving body speed is an insignificant fraction of the speed of light, we don’t need the degree of correctness provided by the Special Relativity Theory.

On the other hand our well-known Global Positioning System (also known as GPS) makes an elegant and crucial use of relativity theories to guarantee its very accurate precision, showing that each context dictates the required level of detail in the theoretical model.

What if aircraft were designed by the best purists in each discipline?


The concept and development of these amazing machines, weighing sometimes hundreds of tons, flying away from the ground and crossing oceans in a few hours, is the result of the hard work of thousands of people in several different disciplines. To enumerate just a part of them: aerodynamics, structures, flight controls, propulsion, electrical systems, hydraulic systems, etc.
If each of these disciplines worked strictly to attain their own requirements, we would have aircraft like the ones shown in the picture below, and none of them would fulfil all the requirements of the product at the same time.

As you can imagine, aircraft design is a huge effort of tuning the needs of all the disciplines together in one single product. And this has been working successfully for more than a century now.

Source: The Aerospace Design Process. A note from a reviewer friend: "I would prefer to see the glider associated with the aerodynamics and a parachute associated with noise" (Cap. Kev)

A bright tool called Object Oriented Programming (OOP)


The Object Oriented Programming (or shortly OOP), roughly speaking, is a way of programming such that the world is modelled by the concept of “objects”. An object encapsulates data (through the “attributes” concept) and the operations allowed over this data (through the “methods” concept). This paradigm, among many other features, allows the program to evolve and to be maintained more smoothly, using concepts like data encapsulation and code re-use, for example.

As far as I know, there is no hard, formal, mathematical definition for OOP (or any other programming paradigm), such that we say a programming language has more or less OOP features. Once I heard from a co-worker “if there is no explicit private/public keywords in Python, then it is not an object-oriented language”. By then, I started to think how an extreme purist view can make one blind. What are the questions we are really interested in answering? (1) Do we need an object-oriented language? or (2) do we need a software easy to be maintained and to evolve? I can assure you that is perfectly acceptable to use a multi-paradigm language that is not fully object-oriented but fulfil nicely your requirements. On the contrary, doing so, means studying carefully all the extent of your problem and then proposing a solution and not the other way round, i.e., picking up one solution and trying to make it fit to your problem.

For instance, if you have assessed all the main risks of your project, all the main test case scenarios (including the very critical ones) and mapped them in your development process, it really doesn’t matter what language or paradigm you are supposed to use. Unless your customer has a weird requirement that says you have to use a language with an explicit private/public keyword, you can use whatever your team is comfortable with. I would not be impressed with a very well written pure C code, implementing objects (yes, try to use function pointers into your structs and have fun).

The emerging agile philosophy all around


In this very well-written comparison between purists and pragmatist approaches, the BreakingSmart makes an enlightening statement about the agile philosophy: “today, the then-minority tradition of pragmatic hacking has matured into agile development, the dominant modern approach for making software. But the significance of this bit of history goes beyond the Internet. Increasingly, the pragmatic, agile approach to building things has spread to other kinds of engineering and beyond, to business and politics” (learn more).

Once the trial-and-error approach got the cost of a click on a computer, followed by results a few seconds later on the screen, hands-on creativity could walk along with productivity. But this triggers the uncle Ben’s warning: “with great power there must also come great responsibility”. It hasn’t been rare to hear that the agile methodologies are these magical tools which will make any development fit on any schedule. This is tricky and may possibly drive to an erratic trial-and-error approach without a minimal of planning or modeling to support it. As the examples above could make clear (hopefully) a successful pragmatic approach comes with sufficient background knowledge to assure which pieces of information can be discarded without compromising the final solution.