Bode Plots in Python

This notebook will go through a pratical example of how to do bode plots in Python and how to find several other control related charateristics. We’ll look at an actual transfer function and derive all the information from that.

Let’s start off by importing the Python libraries that we are going to need. If you are not able to import the library, please install it using pip in your terminal or command prompt.

import matplotlib.pyplot as plt
import matplotlib
import control
from scipy import signal
from control import tf

%matplotlib inline
matplotlib.rcParams['figure.facecolor'] = 'w'

Now let’s define the transfer function that we are going to look at. We’ll look at an open loop transfer function from a mechanical system that is regulated. The method will work for transfer function in mechanics and electronics though

$$OL(s)=20\cdot \frac{300}{(0,01s+1)(0,1s+1)}\cdot 0,0033$$

Now to actually use the python libraries, we want the transfer function to be on a different form. The form should be $s^2+s+k$. So we’ll make sure to dissolve the parantheses in the denominator

$$OL(s)=20\cdot \frac{300}{0,001s^2+0,011s+1}\cdot 0,0033$$

Now that we have the expression in the correct form, we can put it into our first Python function. The function comes from the classsignal from the scipy library. We imported it just before by saying from scipy import signal. Now we are going to use a method from that class called lti. lti takes in several parameters, but we are only going to use two of them in this example. The first parameter is the numerator of our transfer function. The input should be in brackets [] and we are going to multiply everything together. So the first input would be [20*300*0.0033].

The second input will be the denominator of our transfer function. So we have $0,001s^2+0,11+1$. If we remove the $s$’s we get the second parameter which would then be [0.001, 0.11, 1]. As you can see we seperate the values by a comma.

Now let’s put all of this into our class method and assign it to the variable system

system = signal.lti([20*300*0.0033], [0.001, 0.11, 1])

Now the signal.lti doesn’t do much for us. We want to have a bode plot, so we’ll have to use yet another method from the signal class to create the bode plot. The method signal.bode will do the trick. This returns three things: a frequency array in rad/s, a magnitude array in dB, and a phase array in degrees.

So when we call this method we want to assign in to three different variables. The first return value will be the frequency array and we’ll assign that to the variable w. The next return value will be the magnitude array and we’ll assign that to the variable mag. The last return value will be the phase array and we’ll assign that to the variable phase.

We’ll also need to tell the program the range of frequencies we want to look at (this is optional though). To do this we assign the variable r to the built in python function range(0, 10000) which will give us a bode plot of frequencies between 0 and 10.000. So when we call signal.bode we’ll first input our variable system. Then we’ll set the optional parameter w equal to r.

Let’s try it out

r = range(0, 10000)
w, mag, phase = signal.bode(system, w=r)

Okay so now we have our variables w, mag and phase assigned. But they just contain a really long list of a lot of numbers. What we want to do then is to plot all those numbers. To do that, we’ll use the pyplot class from the matplotlib library which we imported earlier by saying import matplotlib.pyplot as plt which means we can now just reference it by calling plt.someMethod().

So we’ll create a new figure by calling plt.figure(). Then we’ll make sure that we can see the logarithmic lines in the plot by calling plt.grid(True, which="both"). Then the last thing to do is the actual plot. We’ll start off by plotting the magnitude. So we want the angular speed on the x-axis, which was w and the magnitude on the y-axis, which was mag. So to plot it we’ll call plt.semilogx(w, mag) since it’s only the x-axis that’s logarithmic.

To plot the phase we’ll do exactly the same. Only when we call plt.semilogx() we’ll input w and phase instead.

To show the plot’s call plt.show()

plt.figure()
plt.grid(True, which="both")
plt.semilogx(w, mag)    # Bode magnitude plot
plt.ylabel("[dB]")
plt.xlabel("[rad/s]")
plt.figure()
plt.grid(True, which="both")
plt.semilogx(w, phase)  # Bode phase plot
plt.ylabel("[degrees]")
plt.xlabel("[rad/s]")
plt.show()

All of this is really great. But it can be kind of messy to actually read values of a semi logarithmic plot. So to get around that, we’ll use another library that really comes in handy in both mechanics and electronics. It’s call control and we imported it at the beginning of the notebook. It gives us a lot of capabilities to find characteristics of our transfer function.

In the beginning of the notebook we called the following from control import tf. Which means that we can now use all the functionality of the tf class. So we’ll create a tf object and assign it to the variable sys. The tf object basicly behaves like the signal.lti object that we used before, so we’ll give it exactly the same input

sys = tf([20*300*0.0033], [0.001, 0.11, 1])

Now comes the fun part. We can use the function control.margin to find out all the margins of our system. Which is all of the things that we would otherwise have had to decifer from our bode plots. So we’ll call control.margin(sys) and assign it to four variables gm, pm, wg, wp.

gm, pm, wg, wp = control.margin(sys)

Now gm will output our Gain Margin which in this case is $\infty$

gm

    inf

pm will be our Phase Margin which will have a value in degrees

pm

    43.51277688046244

In this case we do not have a Frequency For Gain margin which would be wg

wg

    nan

Now the last variable is wp which gives us the Frequency For Phase Margin which we do have in this case. The value is in dB

wp

    123.93293689252836

The control library can also be used to find the poles of the transfer function. Simply call control.pole() and input the sys variable into the parantheses.

control.pole(sys)

    array([-100.+0.j,  -10.+0.j])

In this case we have a pole at $-10$ and $-100$. To find zeroes we can just call control.zero() and input the sys variable. But there are no zeroes in this transfer function.

The last thing we want to see is the DC Gain of the system. This is found by calling control.dcgain() and inputting the sys variable in the parantheses.

control.dcgain(sys)

    19.8

And we are done