11.1.1 Problem

We want to model the longitudinal dynamics of an aircraft.

11.1.2 Solution

The solution is to write the right-hand side function for the aircraft longitudinal dynamics differential equations.

11.1.3 How It Works

We summarized the symbols for the dynamical model in Table 11.1

Our aerodynamic model is very simple. The lift and drag are:

(11.1)

(11.2)

where S is the wetted area, the area that interacts with the airflow and is the area that is counted in computing the aerodynamic forces, and p is the dynamic pressure, the pressure on the aircraft caused by its velocity:

(11.3)

where ρ is the atmospheric density and v is the magnitude of the velocity. Atmospheric density is a function of altitude. For low-speed flight this is mostly the wings. Most books use q for dynamic pressure. We use q for pitch angular rate (also a convention), so we use p for pressure here to avoid confusion.

The lift coefficient, C_L is:

(11.4)

and the drag coefficient, C_D is:

(11.5)

The drag equation is called the drag polar. Increasing the angle of attack increases the aircraft lift, but also increases the aircraft drag. The coefficient k is:

(11.6)

Table 11.1

Aircraft Dynamics Symbols

Symbol	Description	Units
g	Acceleration of gravity at sea-level	9.806 m/s2
h	Altitude	m
k	Coefficient of lift induced drag
m	Mass	kg
p	Dynamic pressure	N/m2
q	Pitch angular rate	rad/s
u	x-velocity	m/s
w	z-velocity	m/s
C L	Lift coefficient
C D	Drag coefficient
D	Drag	N
I y	Pitch moment of inertia	kg-m2
L	Lift	N
M	Pitch moment (torque)	Nm
M e	Pitch moment due to elevator	Nm
r e	Elevator moment arm	m
S	Wetted area of wings (the area that contributes to lift and drag)	m2
S e	Wetted area of elevator	m2
T	Thrust	N
X	X force in the aircraft frame	N
Z	Z force in the aircraft frame	N
α	Angle of attack	rad
γ	Flight path angle	rad
ρ	Air density	kg/m3
θ	Pitch angle	rad

where 𝜖 ₀ is the Oswald efficiency factor, which is typically between 0.75 and 0.85. AR is the wing aspect ratio. The aspect ratio is the ratio of the span of the wing to its chord. For complex shapes, it is approximately given by the formula:

(11.7)

where b is the span and S is the wing area. Span is measured from wingtip to wingtip. Gliders have very high aspect ratios and delta-wing aircraft have low aspect ratios.

The aerodynamic coefficients are nondimensional coefficients that when multiplied by the wetted area of the aircraft, and the dynamic pressure, produce the aerodynamic forces.

The dynamical equations, the differential equations of motion, are [5]:

(11.8)

(11.9)

(11.10)

(11.11)

m is the mass, u is the x-velocity, w is the z-velocity, q is the pitch angular rate, θ is the pitch angle, T is the engine thrust, 𝜖 is the angle between the thrust vector and the x-axis, I_y is the pitch inertia, X is the x-force, Z is the z-force, and M is the torque about the pitch axis. The coupling between x and z velocities is due to writing the force equations in the rotating frame. The pitch equation is about the center of mass. These are a function of u, w, q, and altitude, h, which is found from:

(11.12)

The angle of attack, α is the angle between the u and w velocities and is:

(11.13)

The flight path angle γ is the angle between the vector velocity direction and the horizontal. It is related to θ and α by the relationship:

(11.14)

This does not appear in the equations, but it is useful to compute when studying aircraft motion. The forces are:

(11.15)

(11.16)

The moment, or torque, is assumed because of the offset of the center-of-pressure and center of mass, which is assumed to be along the x-axis:

(11.17)

where c_p is the location of the center of pressure. The moment due to the elevator is:

(11.18)

S_e is the wetted area of the elevator and r_E is the distance from the center of mass to the elevator. The dynamical model is in RHSAircraft. The atmospheric density model is an exponential model and is included as a subfunction in this function. RHSAircraft returns the default data structure if no inputs are given.

We will use a model of the F-16 aircraft for our simulation. The F-16 is a single engine supersonic multi-role combat aircraft used by many countries. The F-16 is shown in Figure 11.2.

The inertia matrix is found by taking this model, distributing the mass amongst all the vertices, and computing the inertia from the formulas

(11.19)

(11.20)

(11.21)

where N is the number of nodes and r_k is the vector from the origin (which is arbitrary) to node k.

There are many limitations to this model. First of all, the thrust is applied immediately with 100% accuracy. The thrust is also not a function of airspeed or altitude. Real engines take some time to achieve the commanded thrust and the thrust levels change with airspeed and altitude. In the model, the elevator also responds instantaneously. Elevators are driven by motors, usually hydraulic, but sometimes pure electric, and they take time to reach a commanded angle. In our model, the aerodynamics are very simple. In reality, lift and drag are complex functions of airspeed and angle of attack and are usually modeled with large tables of coefficients. We also model the pitching moment by a moment arm. Usually, the torque is modeled by a table. No aerodynamic damping is modeled although this appears in most complete aerodynamic models for aircraft. You can easily add these features by creating functions:

11.2.1 Problem

We want to determine the equilibrium state for the aircraft. This is the orientation at which all forces and torques balance.

11.2.2 Solution

The solution is to compute the Jacobian for the dynamics. The Jacobian is a matrix of all first-order partial derivatives of a vector valued function, in this case the dynamics of the aircraft.

11.2.3 How It Works

We want to start every simulation from an equilibrium state. This is done using the function EquilibriumState. It uses fminsearch to minimize:

(11.22)

given the flight speed, altitude, and flight path angle. It then computes the elevator angle needed to zero the pitch angular acceleration. It has a built-in demo for equilibrium level flight at 10 km.

CostFun is the cost functional given below.

The vector of values is the first input. Our first guess is that thrust equals drag. The vertical velocity and thrust are solved for by fminsearch. fminsearch searches over thrust and vertical velocity to find an equilibrium state.

The results of the demo are:

The initial and final costs show how successful fminsearch was in achieving the objective of minimizing the w and u accelerations.

11.3.1 Problem

We want to simulate the aircraft.

11.3.2 Solution

The solution is to create a script that calls the right-hand side of the dynamical equations, RHSAircraft in a loop, and plot the results.

11.3.3 How It Works

The simulation script is shown below. It computes the equilibrium state, then simulates the dynamics in a loop by calling RungeKutta. It applies a disturbance to the aircraft. It then uses PlotSet to plot the results.

The applied external acceleration puts the aircraft into a slight climb with some noticeable oscillations.

The simulation results are shown in Figure 11.3. The aircraft climbs steadily. Two oscillations are seen. A high frequency one primarily associated with pitch and a low frequency one with the velocity of the aircraft.

11.4.1 Problem

We are going to implement a neural net so that our aircraft control system can learn. We need an activation function to scale and limit measurements.

11.4.2 Solution

Use a sigmoid function as our activation function.

11.4.3 How It Works

The neural net uses the following sigmoid function

(11.23)

The sigmoid function with k = 1 is plotted in the following script.

Results are shown in Figure 11.4.

11.5.1 Problem

We want to use a neural net to add learning to the aircraft control system.

11.5.2 Solution

Use a sigma-pi neural net function. A sigma-pi neural net sums the inputs and products of the inputs to produce a model.

11.5.3 How It Works

The adaptive neural network for the pitch axis has seven inputs. The output of the neural network is a pitch angular acceleration that augments the control signal coming from the dynamic inversion controller. The control system is shown in Figure 11.5. The left-most box produces the reference model given the pilot input. The output of the reference model is a vector of the desired states that are differenced with the true states and fed to the PID controller and the neural network. The output of the PID is differenced with the output of the neural network. This is fed into the model inversion block that drives the aircraft dynamics.

The sigma-pi neural net is shown in Figure 11.6 for a two-input system.

The output is:

(11.24)

The weights are selected to represent a nonlinear function. For example, suppose we want to represent the dynamic pressure:

(11.25)

We let x ₁ = ρ and x ₂ = v ². Set and all other weights to zero. Suppose we didn’t know the constant . We would like our neural net to determine the weight through measurements. Learning for a neural net means determining the weights so that our net replicates the function it is modeling. Define the vector z, which is the result of the product operations. In our two-input case this would be:

(11.26)

c is a constant. The output is:

(11.27)

We could assemble multiple inputs and outputs:

(11.28)

where z_k is a column array. We can solve for the weights w using least squares given the outputs, y, and inputs, x. Define the vector of y to be Y  and the matrix of z to be Z. The solution for w is:

(11.29)

The least squares solution is:

(11.30)

This gives the best fit to w for the measurements Y  and inputs Z. Suppose we take another measurement. We would then repeat this with bigger matrices. As a side note, you would really compute this using an inverse. There are better numerical methods for doing least squares. MATLAB has the pinv function. For example:

As you can see, they all agree! This is a good way to initially train your neural net. Collect as many measurements as you have values of z and compute the weights. Your net is then ready to go.

The recursive approach is to initialize the recursive trainer with n values of z and y.

(11.31)

(11.32)

The recursive learning algorithm is:

(11.33)

(11.34)

(11.35)

RecursiveLearning demonstrates recursive learning or training. It starts with an initial estimate based on a four-element training set. It then recursively learns based on new data.

Figure 11.7 shows the results. After an initial transient, the learning converges. Every time you run this you will get different answers because we initialize with random values.

You will notice that the recursive learning algorithm is identical in form to the Kalman Filter given in Section 4.​1.​3. Our learning algorithm was derived from batch least squares, which is an alternative derivation for the Kalman Filter.

11.6.1 Problem

One issue with a sigma-pi neural network is the number of possible nodes. For design purposes, we need a function to enumerate all possible sets of combinations of inputs. This is to determine the limitation of the complexity of a sigma-pi neural network.

11.6.2 Solution

Write a combination function that computes the number of sets.

11.6.3 How It Works

In our sigma-pi network, we hand-coded the products of the inputs. For more general code, we want to enumerate all combinations of inputs. If we have n inputs and want to take them k at a time, the number of sets is:

(11.36)

The code to enumerate all sets is in the function Combinations.

This handles two special cases on input and then calls itself recursively for all other cases. Here are some examples:

You can see that if we have four inputs and want all possible combinations, we end up with 14 in total! This indicates a practical limit to a sigma-pi neural network, as the number of weights will grow fast as the number of inputs increases.

11.7.1 Problem

We need a sigma-pi net function for general problems.

11.7.2 Solution

Use a sigma-pi function.

11.7.3 How It Works

The functionality is distributed among sub-functions called from the switch statement.

The demo shows an example of using the function to model dynamic pressure. Our inputs are the altitude and the square of the velocity. The neural net will try to fit:

(11.37)

to

(11.38)

We first get a default data structure. Then we initialize the filter with an empty x. We then get the initial weights by using batch learning. The number of columns of x should be at least twice the number of inputs. This gives a starting p matrix and initial estimate of weights. We then perform recursive learning. It is important that the field kSigmoid is small enough so that valid inputs are in the linear region of the sigmoid function. Note that this can be an array so that you can use different scalings on different inputs.

The batch results are as follows for five examples of dynamic pressures at low altitude. As you can see the truth model and neural net outputs are quite close:

The recursive learning results are shown in Figure 11.8. The results are pretty good over a wide range of altitudes. You could then just use the “update” action during aircraft operation.

11.8.1 Problem

We want a PID controller to control the aircraft.

11.8.2 Solution

Write a function to implement PID control. The inputs will be the pitch angle error.

11.8.3 How It Works

Assume we have a double integrator driven by a constant input:

(11.39)

where u = u_d + u_c.

The result is:

(11.40)

The simplest control is to add a feedback controller

(11.41)

where K is the forward gain and τ is the damping time constant. Our dynamical equation is now:

(11.42)

The damping term will cause the transients to die out. When that happens, the second derivative and first derivatives of x are zero and we end up with an offset:

(11.43)

This is generally not desirable. You could increase K until the offset were small, but that would mean your actuator would need to produce higher forces or torques. What we have at the moment is a proportional derivative (PD) controller. Let’s add another term to the controller:

(11.44)

This is now a proportional integral derivative (PID) controller. There is now a gain proportional to the integral of x. We add the new controller, and then take another derivative to get:

(11.45)

Now in steady state:

(11.46)

If u is constant, the offset is zero. Define s as the derivative operator.

(11.47)

Then:

(11.48)

Note that:

(11.49)

where τ_d is the rate time constant, which is how long the system will take to damp and τ_i is how fast the system will integrate out a steady disturbance.

where s = jω and . The closed loop transfer function is:

(11.50)

The desired closed loop transfer function is:

(11.51)

or

(11.52)

The parameters are:

(11.53)

(11.54)

(11.55)

This is a design for a PID. However, it is not possible to write this in the desired state space form:

(11.56)

(11.57)

because it has a pure differentiator. We need to add a filter to the rate term so that it looks like:

(11.58)

instead of s. We aren’t going to derive the constants and will leave it as an exercise for the reader. The code for the PID is in PID.

It is interesting to evaluate the effect of the integrator. This is shown in Figure 11.9. The code is the demo in PID. Instead of numerically integrating the differential equations we convert them into sampled time and propagate them. This is handy for linear equations. The double integrator equations are in the form:

(11.59)

(11.60)

This is the same form as the PID controller.

It takes about 2 minutes to drive x to zero, which is close to the 100 seconds specified for the integrator.

11.9.1 Problem

We want to control the pitch angle of an aircraft with a PID controller.

11.9.2 Solution

Write a script to implement the controller with the PID controller and pitch dynamic inversion compensation.

11.9.3 How It Works

The PID controller changes the elevator angle to produce a pitch acceleration to rotate the aircraft. The elevator is the moveable horizontal surface that is usually on the tail wing of an aircraft. In addition, additional elevator movement is needed to compensate for changes in the accelerations due to lift and drag as the aircraft changes its pitch orientation. This is done using the pitch dynamic inversion function. This returns the pitch acceleration, which must be compensated for when applying the pitch control.

The simulation incorporating the controls is AircraftSim below. There is a flag to turn on control and another to turn on the learning control. We command a 0.2-radian pitch angle using the PID control. The results are shown in Figure 11.10, Figure 11.11 and Figure 11.12.

The maneuver increases the drag and we don’t adjust the throttle to compensate. This will cause the airspeed to drop. In implementing the controller we neglected to consider coupling between states, but this can be added easily.

11.10.1 Problem

We want a nonlinear inversion controller with a PID controller and the sigma-pi neural net.

11.10.2 Solution

Train the neural net with a script that takes the angle and velocity squared input and computes the pitch acceleration error.

11.10.3 How It Works

The PitchNeuralNetTraining script computes the pitch acceleration for a slightly different set of parameters. It then processes the delta-acceleration. The script passes a range of pitch angles to the function and learns the acceleration. We use the velocity squared as an input because the dynamic pressure is proportional to the velocity squared. The base acceleration (in dRHSL) is for our “a-priori” model. dRHS is the measured values. We assume that these are obtained during flight testing.

The script first finds the equilibrium state using EquilibriumState. It then sets up the sigma-pi neural net using SigmaPiNeuralNet. PitchDynamicInversion is called twice, once to get the model aircraft acceleration aM (the way we want the aircraft to behave) and once to get the true acceleration aT. The delta acceleration, dA is used to train the neural net. The neural net produces aNN. The resulting weights are saved in a .mat file for use in AircraftSim. The simulation uses dRHS, but our pitch acceleration model uses dRHSL. The latter is saved in another .mat file.

As can be seen, the neural net reproduces the model very well. The script also outputs DNN. mat, which contains the trained neural net data.

11.11.1 Problem

We want to demonstrate our learning control system for controlling the longitudinal dynamics of an aircraft.

11.11.2 Solution

Enable the control functions to the simulation script described in AircraftSimOpenLoop.

11.11.3 How It Works

After training, the neural net in the previous recipe we set addLearning to true. The weights are read in. We command a 0.2-radian pitch angle using the PID learning control. The results are shown in Figure 11.14, Figure 11.15, and Figure 11.16. The figures show without learning control on the left and with learning control on the right.

Learning control helps the performance of the controller. However, the weights are fixed throughout the simulation. Learning occurs prior to the controller becoming active. The control system is still sensitive to parameter changes since the learning part of the control was computed for a pre-determined trajectory. Our weights were determined only as a function of pitch angle and velocity squared. Additional inputs would improve the performance. There are many opportunities for you to try to expand and improve the learning system.