State-Space Methods Part 3: Flight Controller System Identification

Introduction

Typically we think of system identification as the precursor to designing a control system, but for some systems in the early design phase there is no physical system yet. This article we will go through how to use system identification tools to identify controller dynamics rather than a plant, given a frequency response of the controller. We will then formulate the controller in a discrete state-space form that can be easily programmed into a digital system. For the purposes of this example we will be making the assumption that our controller is a single variable controller that controls a single longitudinal output. For the example aircraft this is not an accurate assumption but is necessary as explained in the analysis at the end of the article.

Controller

The controller that we will look at is the refueling pitch controller for the B-2 bomber. Refueling is one of the more difficult maneuvers to accomplish in the B-2 as it required hand-flying the aircraft during a tight control tracking task. We will look at a bode plot of the refueling controller and then reverse engineer a single-input single-output controller of similar performance.

This bode plot shows the feedback dynamics for the pitch controller during a refueling mode.

We can see that the bandwidth of the controller is 3.7 rad/s and the crossover frequency is 6 rad/s. This large bandwidth which is required for tightly-controlled tracking tasks like refueling. We can also see the rapid phase drop-off after the crossover frequency. For the rest of the article will cover how to make a digital controller with this desired frequency response characteristic.

Matlab System Identification

Now that we have frequency and magnitude data from a bode plot of the aircraft pitch controller, we can create a transfer function for it. Once we’ve interpolated the phase and magnitude data to the common set of frequencies we can view the new interpolated Bode plot to verify that we have the correct data points to feed to our identification tool.

We will be using the Matlab system identification toolbox as it has powerful features for model identification. Even though the tool has state-space identification modes We will be using the transfer function identification mode so that later on we can see how to implement the controller if we only have a desired transfer function.

Now we can open the Matlab system identification tool and input our data from the bode plot.

Now, back in the system identification window I select “estimate Transfer Function Models”

Our first model is a continuous system with one pole and no zeros. Lets name it p1z0.

When we check the identified model result we can see that the model fit is alright (70%), but I think that we can do a bit better.

I can go through and create different order models to fit with our data. Let’s do that for models up to the 5th-order.

Let’s compare our first-order models. We can deselect the others by clicking on them so that their lines are no longer bold. It looks like our best fit first-order model is the one with 1 pole and 1 zero with a fit of 71.12%.

Next lets look at our 2nd-order models. It looks like they are getting a much better fit with the best models at 75.65%.

Now let’s look at the 3rd-order models. Now we’re getting somewhere with model fits in the 80% range.

Let’s see how our 4th-order models are doing. Now at 88.42%, our model is getting pretty close.

Let’s see if we can do better with a 5th-order model. Not bad at all, the best model fits are above 96%.

Finally we can look at the 6th and 7th order models to see if we get any significant improvements. The improvements are there, but there are definitely diminishing returns.

6th-order models
7th-order models

By plotting the models out to the 16th-order we can see that we get drastically diminishing returns for the more complicated models. Past a 5th or 6th-order model, it seems like we may be overfitting to the data, especially because we don’t have a lot of other data to add to the model in the form of other test runs.

Now let’s take a closer look at the identified transfer function for the 6th-order model.

This looks a little rough so we can use Matlab’s zpk() function to factor this transfer function. When we factor it, we get a more intelligible transfer function
$$C(s)=\frac{0.00083601(s-120.9)(s+1.127)(s^2-0.07782s+0.007209)(s^2-26.33s+477.3)}{(s+29.4)(s-0.08225)(s+0.07858)(s-0.01802)(s^2+4.602s+14.15)}$$

We can see that there is one complex conjugate pair of poles in the denominator \((s^2+4.602s+14.15)\) this will be important later. Aside from that we can see that we have two sets of complex zeroes in the numerator.

When we look at the root-locus plot of this transfer function to gain insight into some of the terms that we can see in the denominator of the transfer function. We can see the motion of our two complex conjugate poles as shown by the green and red traces.

We can also see an unstable zero on the right hand side. In some controllers this zero is used to cancel out an unstable pole.

When we zoom in on these we can get a better idea of the structure closer to the origin. Notice how much smaller the scale is compared to the previous diagram. We can see that there are unstable poles in the right-hand plane. Unstable poles may be necessary for controlling an unstable plant, but are usually not a good idea to put in a controller.

We can also see that we have a set of complex zeros that correspond to poles that lie In the unstable plant on the real axis. This is the first \(s^2\) term in the numerator of the transfer function.

Now if this controller was acceptable and had acceptable closed-loop performance characteristics, we could implement the controller using the same state-space methods used to model any dynamic system.

Hypothetical Controller Synthesis

If we wanted to implement this controller, one way to do it is to use the modal state-space form. Modal form is a special form of the state-space representation where we diagonalize the A matrix. This is a good example of the output equation for state-space models as we can keep the outputs of the model in standard dimensions such as position and velocity, while having the actual states be in a non-dimensional eigenvalue or otherwise non-intuitive unit form.

The SISO case decomposes the A matrix of a state-space model into its eigenvectors
$$\begin{bmatrix}
x_1(n+1)\ x_1(n+1)\\
\vdots\\
x_N(n+1)
\end{bmatrix}=
\begin{bmatrix}
\lambda_1 & 0 & \dots& 0 \\
0 & \lambda_2 & \dots & 0 \\
\vdots & \vdots & \ddots & \vdots \\
0 & 0 & \dots &\lambda_N
\end{bmatrix}
\begin{bmatrix}
x_1(n)\ x_1(n)\\
\vdots\\
x_N(n)
\end{bmatrix}+
\begin{bmatrix}
b_1\\ b_2\\ \vdots \\ b_N
\end{bmatrix}u(n)$$
$$y=[C_1\ C_2\ \dots \ C_N]x(n)+Du(n)$$
The diagonal A matrix means that the new states are decoupled into N single-pole systems. Each of these systems represents a single mode. This only works if the eigenvectors are independent. The B matrix describes how the modes of the system are driven by the input. The C matrix then determines how the modes are mixed to form the output.

This formulation of a digital controller was used in the VAAC Harrier aircraft program. For that aircraft the controller order was reduced to an 8th-order controller and then implemented in modal form to minimize the number of additions and multiplications required per controller update.

The modalreal() function in Matlab can calculate this for us. Modal form is a special form of the state-space representation where we diagonalize the A matrix. This is a good example of the output equation for state-space models as we can keep the outputs of the model in standard dimensions such as position and velocity, while having the actual states be in a non-dimensional or otherwise non-intuitive unit form.

From the data that we extracted from the bode plot, the modal-state-space form of the controller is shown in the image to the right. We can see that the matrix is not entirely diagonal, and instead it has a block-diagonal structure. This means that we aren’t able to fully decouple the states. Complex conjugate pole pairs form block diagonal elements in the A matrix. With one complex pole pair in the denominator, we have one block-diagonal element in the A matrix.

We can also use the c2d() method to convert from a continuous controller to a discrete controller that has a refresh rate of 50 Hz.

This modification is important as digital controllers have a discrete sample rate that must be accounted for when designing the controller.

When we compare the frequency response of the discrete controller to the continuous transfer function shown above, we can see that it agrees well in magnitude, up to the higher frequencies, and it agrees pretty well in phase, with a slight drop off again at higher frequencies. The vertical line stops the discrete model as any signals above this frequency are aliased down to the lower frequency range due to the digital sampling.

Analysis

Without an accurate plant model, we can’t really examine the closed-loop performance of this controller. It is also not likely that for such a complicated airframe would have single-channel controllers. The actual implementation is guaranteed to be much more complicated and likely involves multi-dimensional control architectures.

The B-2 pitch controller uses a Load-factor(\(N_z\)) and Pitch-rate(\(q\)) PI controller design (NZQPPI). This allows it to have level 1 handling qualities throughout the flight envelope. A pure PI controller has a single zero to the left of the origin and a single pole at the zero. The 6th-order controller that we created is much more complicated than that so it is difficult to compare it.

For this example I wanted to show how we can easily generate a state-space model given a desired frequency response, and that this doesn’t have to be a model of the physical system, it can be the model of the controller instead. In future articles we will synthesize controllers in this form and examine their performance on actual systems.

Sources

Keep In Touch

Twitter Linkedin-in