Data Science and Computing with Python for Pilots and Flight Test Engineers

Introduction

In this lesson we will explain how to create a linear time-invariant (LTI) control system described by an ordinary linear differential equation with constant coefficients in Python, such that we can do system analysis and design with it. This lesson will cover some of the theoretical background, and the next lesson, Control System Creation – Code, will show you how to actually implement all of this in Python code, which turns out to be pretty straightforward. Many links to various references are provided in the text, which we encourage you to follow for additional insights, as interest demands.

There are at least two convenient ways to perform such a system setup in Python:

• with the scipy.signal submodule of the SciPy module (which will use its own Python syntax), or
• with the Python Control Systems Library python-control (which has a MATLAB compatibility submodule control.matlab, which implements in Python MATLAB-style functionality and syntax for people who are familiar with the MATLAB Control System Toolbox).

We will show you how to work with both. The first solution is convenient, because we always have the SciPy library already installed anyway. The second solution comes handy, if you are taking a class or are working with people who use MATLAB, and you want to use the same functions and syntax as they do. Both of these solutions are available in Python for free, whereas you will notice at the link above that for the actual MATLAB Control Systems Toolbox you would have to pay, of course.

Whichever of the two libraries above you decide to use, there are three basic representations to choose from during the creation of such a control system. To create a new system, you will need to specify one of the following:

• 2 – tf: the polynomial coefficients of the numerator and denominator of the transfer function,
• 3 – zpk: the zeros, poles, and gain of the transfer function,
• 4 – ss: the state-space representation of the system’s differential equation (here we will need to specify the \(A, B, C, D\) matrices).

The scipy.signal submodule will be able to tell automatically, which of these three options we are choosing during system creation, because for each of them we will have to provide a different number of arguments to the system creation function of the submodule, as is indicated by the numbers in front of the three options above. The first option (tf) has two arguments (two lists of coefficients, one each for the numerator and denominator, respectively). The second option (zpk) has three arguments (list of zeros, list of poles, gain), while the third option has four arguments (the matrices \(A\), \(B\), \(C\), and \(D\)).

Let us first take a look at what such a control system looks like and then how to set it up in code in all three possible representations above. In subsequent lessons, we will perform system analysis and answer questions about the system, but in this lesson (and its code companion lesson) for now we will only learn how to create a system in your Python code and some of the theoretical background behind it.

System Description

Differential Equation describing the Control System

Let us take as a specific example to work with a linear time-invariant control system governed by the second-order linear ordinary differential equation with constant coefficients:
$$ \ddot x(t) + 2 \dot x(t) + 17 x(t) = f(t) $$

Transfer Function

The transfer function \(H(s)\) of a system is defined as the ratio of system output to system input in Laplace space \(H(s):=X(s)/F(s)\).

This is equivalent to saying that the transfer function \(H(s)\) is the Laplace transform of the impulse response function of the system, i.e. the Laplace transform \(X(s)\) of the response \(x(t)\) of the system to an impulse input \(f(t)=\delta(t)\), where \(\delta(t)\) is the Dirac delta function. (This equivalence, i.e. the fact that \(H(s)=X(s)\) if \(F(s)=\delta(s)\), with \(\delta(s)\) being the Laplace transform of impulse input \(f(t)=\delta(t)\) is readily seen, if one remembers that \(\delta(s)=1\) (the Laplace transform of the Dirac delta function (\delta(t)\) is equal to 1).

The transfer function \(H(s)\) for the above example control system is obtained readily from the Laplace transform of its defining differential equation and turns out to be:

$$ H(s) = \frac{1}{s^2 + 2s + 17}. $$

Coefficients of Transfer Function Numerator and Denominator

It is very easy to read of the lists of coefficients of the numerator and denominator of the transfer function above. The numerator, N, is simply \(N=1\). Therefore the list of coefficients is [1]. The denominator \(D\) (not to be confused with the state-space matrix \(D\) used below), is \(D = s^2 + 2s + 17\). The list of coefficients of this polynomial in \(s\) is [1, 2, 17]. Notice that we order them in the list starting from the highest order term to the lowest order term. If a coefficient of the polynomial is zero, we must still indicate is as zero in the list (unless all higher order coefficients are also zero – the list starts with the highest non-zero coefficient).

Zeros, Poles, and Gain

Sometimes the numerator and denominator of the transfer function might not be written as an expanded polynomial with terms multiplied out, but rather as a product. (What we are referring to here is that you can, for instance, write the polynomial \(s^2+2s+2\) also as a product in the form \((s+2)(s+1)\), which makes it very easy to see that its roots are \(s=-2\) and \(s=-1\), whereas to get the expanded form, you would need to multiply out the brackets first.)

In this case, it is easier to specify the transfer function by its zeros and poles, rather than by its polynomial coefficients. The zeros are the roots of the numerator (i.e. the values of \(s\) for which the numerator becomes zero). The poles of the transfer function are the roots of the denominator (i.e. the values of \(s\) for which the denominator becomes zero). The gain is just the constant prefactor of the transfer function (when we write the transfer function in a way that the highest order coefficient of the denominator is equal to 1).

In the above example system, we see that the transfer function has no zeros (not even complex ones), because the numerator is the constant polynomial 1.

But the transfer function does have poles. They are the roots of the polynomial of the denominator (the so called characteristic polynomial of the transfer function or differential equation). These roots are located at \(-1+4i\) and \(-1-4i\), where \(i\) denotes the imaginary number defined by equation \(i^2=-1\) (see complex numbers, if you are unfamiliar). In order to see this in retrospect, note that \((s+1-4i)(s+1+4i) = s^2 + 2s + 17\), which is exactly the denominator of our transfer function above.

Every non-constant polynomial with complex coefficients (which includes real coefficients) will always have at least one complex root (see the fundamental theorem of algebra), and \(s\) in our case is a complex variable, because it comes from the definition of the Laplace transform. For a second-order polynomial, a complex root (with non-zero imaginary part) will always appear as a complex conjugate pair, just as we see above.

In summary, specifying the zeros, poles, and the gain of the transfer function is also a way to define the system, if one does not want to bother to compute the transfer function polynomial coefficients.

State-Space Representation

Yet another way to represent a linear time-invariant control system is the state-space representation. It is founded on the fact that one can turn every higher-order linear ordinary differential equation into a system of first-order linear differential equations.

We can then write our single higher-order linear ordinary differential equation as a system of equations, and this system of equations can be represented as one equation in matrix form, with suitable definition of the matrices and vectors:

\begin{eqnarray}
\dot{\mathbf{x}}(t) &=& A\mathbf{x}(t) + B\mathbf{u}(t)\\
\mathbf{y}(t) &=& C\mathbf{x}(t) + D\mathbf{u}(t)
\end{eqnarray}

\(A\) is the state matrix (plant matrix), governing the free response of the system. \(B\) is the input matrix (control matrix), indicating how control variables affect the time evolution of the system (being passed into the system).

\(C\) is the output matrix, which constructs output variables from the system state (the output variables of the system do not have to be the state variables of the system in general. \(D\) is feedthrough (feedforward) matrix, if any control input influences outputs directly, without passing the controls through the system first.

In our cases, we will have in general \(C\) being the identity matrix, because we want our output variables simply to be the state variables of the system, and \(D\) will be the zero matrix, because we will not have any control variable changing the output directly, without passing through the system first.

The state space representation of the above system is

\begin{equation} \underbrace{\begin{pmatrix}\ddot x \\ \dot x\end{pmatrix}}_{=:\dot{\mathbf{x}}(t)} = \underbrace{\begin{pmatrix}-2 & -17 \\ 1 & 0\end{pmatrix}}_{=:A}\underbrace{\begin{pmatrix}\dot x \\ x\end{pmatrix}}_{=:\mathbf{x}(t)} + \underbrace{\begin{pmatrix}1 \\ 0\end{pmatrix}}_{=:B} \underbrace{f(t)}_{=:\mathbf{u}(t)} \end{equation}

\begin{equation} \underbrace{x(t)}_{=:\mathbf{y}(t)} = \underbrace{\pmatrix{0 & 1}}_{=:C} \underbrace{\pmatrix{\dot x \\ x}}_{=:\mathbf{x}(t)} + \underbrace{\pmatrix{0}}_{=:D} \underbrace{f(t)}_{=:\mathbf{u}(t)} \end{equation}

The first row is just the above differential equation. The second row just states \(\dot x = \dot x\). Variable \(x\) has no influence on this, and neither does \(f(t)\), so those two coefficients in the lower row are zero.

For the output, we are only interested in \(x(t)\) here. If we output also \(\dot x(t)\), we would get another transfer function.