To teach myself how to use Python instead of MATLAB, I could think of no better way than to dig up one of my old MATLAB homework assignments and do the whole thing in Python. My way of coding things may not be the best, this is my first attempt at using Python, but I was able to do the entire assignment and this tutorial will share how I did it. Let’s start by going over the assignment.
For a lab in class, we had a simple, physical second order system set up and LabView recorded the time and position of our simple mass. We checked both the free and step responses of this system and LabView recorded the data.
The goal for this assignment is to find which parameters for a model of a second order system will match the ones physically measured in this lab. There are four main parts to the assignment:
- Use a least-squares fit to identify parameters
- Estimate the parameters using logarithmic decrement
- Add dry friction to the model
- See if the model matches when a step input is applied
Getting Started in Python
In this tutorial, I’m going to assume you’re using Spyder as installed by the pythonxy package, although I’m sure everything works just peachy using any other appropriate Python set up. (Spyder is my favorite IDE, my next favorite that I’m forced to use when I boot into Linux is the PythonTookit (PTK), Eclipse with PyDev also works well, but I prefer the others).
When using Spyder, it’s very helpful to test the code you’re going to write in the live console before you write it into your script. There are a couple of advantages to doing things this way:
- You always know that you’ve got the syntax correct
- You have live access to your variables in the variables explorer. For some reason the variables don’t stay in the explorer when you execute your script, they get cleared out as soon as your script finishes running (if you run in debug mode, you might be able to pause your code and look at the variables, but I think it’s just easier this way…)
Spyder with some important areas highlighted in red. The bottom right is the live console, the top right can be changed between the Object Inspector and the Variable Explorer. The Object Inspector is really helpful for on-the-fly documentation and help. The Variable Explorer is always handy to see what variables you've used and what kind of values they have.
Now, when you create a file in Spyder for the first time it will automatically add a few lines of code at the beginning, this is mainly for documentation. In Python the pound sign (#) is used for comments (as opposed to % in MATLAB). This isn’t really too big of a change, except that it’s very useful in Unix/Linux systems (possibly OS X as well, I’m not sure) because you can add the line:
#!/usr/bin/python to the top of your script and then it can be run as its own executable! The other thing to note is that the the triple quotes are just a special kind of comment that can be used as documentation for your script.
The first thing we have to do to make our Python script work right is to import the proper libraries. There are a few ways you can do this, I’ll go over a few of the options, since it makes a slight difference as to how you write your code. First let’s import the numpy library two different ways (numpy is a library that contains a lot of the functions we’ll be needing, others we’ll be using will be scipy and pylab).
from numpy import *
This will blindly import everything in the numpy library. This is nice because you will have direct access to all the numpy commands. It’s not so nice because if you import another library the same way, you may have duplicate commands. An alternative is to import numpy like this:
import numpy as np
This will require a little more typing, because every time you use a command from the numpy library you will have to type
np.command() instead of just
command(). While this may be a pain, it’s also very nice because now you can specify which library your commands come from. I find it particularly useful when I’m in linux using PTK because PTK’s Namespace Browser has a live variable display, but it also includes all of your available commands in that display. If you import as np, then instead of listing every command available in numpy, it only lists np, and once you click on np it shows you all available commands from within that library.
There are a few other tricks to importing libraries. Often, commands are nested within a library, and how you import them can affect how you may call the commands later. For example:
from scipy.integrate import odeint
now allows me to use the command
odeint whenever I want, without having to call
integrate.odeint. Finally, you can use an import command on its own. If I load pylab like this:
then every call to a command from within pylab must be accessed like this:
I’m sure there’s a lot more that could be said about importing libraries in Python, but I’ve already said too much about them for the purpose of this tutorial, so let’s go on.
Processing External Data
The first thing we need to do in this program is to load our external data. Both numpy and scipy have a few commands available for this, depending on what format your data is in. Before we can load any data, though, we need to set the proper working directory. Your working directory is shown in the top right, you can browse to a new folder by clicking the folder button, then set the working directory by clicking the button directly to the right of it.
The boxes highlighted in red will help you set your working directory in Spyder
My data is a simple tab-delimited text file, so the numpy function loadtxt works wonderfully.
free_response = loadtxt("free_response.lvm")
Note: type “loadtxt” into the Spyder Object Inspector to look at the documentation for this function. If you scroll down a ways there is a nice “See Also” box that provides some helpful tips for a few other functions that might be useful, particularly scipy.io.loadmat
It’s helpful to have a quick look at this data to see how we’ll want to process it. We can plot the data from the interactive terminal just to see what it looks like.
free_response = loadtxt("free_response.lvm")
Looking at my code to plot this data, you can see a few simple differences between Python and MATLAB.First of all, we access arrays slightly differently than in MATLAB. In MATLAB, parentheses are used (), but in Python we use brackets .
It’s also noteworthy that in Python our arrays start at 0, while in MATLAB they start at 1.
The plot command, taken from Pylab, was modeled after MATLAB’s, so there shouldn’t be much difference there.
Now, when you look at this graph, you can see that the first few points are just stagnant, and that we don’t really want them when we’re dealing with how to fit the model. Let’s delete the first 20 points. Before we can delete a series of points, it’s helpful to go over a couple of useful commands,
arange will create an array based on a start value, an end value, and a given step size, this can be really useful for handling indeces in arrays.
linspace is very similar, except that instead of defining a step size you define the desired number of points, which can be really useful for graphing.
arange to help us delete the first 20 points:
a = arange(0,20,1)
Now we can call the delete command to delete those points from the array:
free_response = delete(free_response, a, 0)
As you type the delete command in the terminal, you can see the documentation above that it takes an array, an object, and an axis. We give it the name of our array, and as an “object” we give it the array of indeces we want to delete, since we want to delete along the first axis, we pass 0 as the axis.
The next step is to re-zero our time vector, which is easy enough to do by subtracting the minimum time from the time vector:
Step 1 – Approximating the model with least-squares
The theory of this method is not really what we’re trying to explain in this tutorial, but the following equations give you an idea of what we’re trying to do with a least-squares fit:
Equations for the least squares method we will use.
To set up our least-squares approximation, we’ll need vectors for position, velocity, and acceleration. We can find these using the
time = free_response[:,0]
pos = free_response[:,1]
vel = diff(pos)/diff(time)
time = delete(time,-1)
accel = diff(vel)/diff(time)
diff command is easy enough to understand, you can learn more about it in the Object Inspector documentation. One interesting thing in the above commands, however, deals with indexing arrays in Python. Passing an index of -1 is simlar to passing
end in MATLAB.
Now we need to adjust the vectors so that they’re all the same size:
time = delete(time,-1)
vel = delete(vel,-1)
pos = delete(pos, [pos.size-1, pos.size-2], None)
Now we get to have a little bit of fun with Python’s way to smash arrays into matrices. We want to stack our velocity and position vectors on top of each other to make a matrix. We do this with the
vstack command. This vertically stacks two arrays on top of each other (yes, there is a very similar
hstack command to stack arrays or matrices next to each other horizontally). We then transpose this just so it will be in the correct orientation for the least squares method we’ll use. You can transpose any array just by adding
.T to the end of it.
A = vstack((vel,pos))
A = A.T
And now we’re ready to run our least squares fit to get some parameters out to model. The least squares function lives inside the numpy.linalg namespace, so we could call it like this:
phi = linalg.lstsq(A,accel)
This won’t give you quite everything in the way you would expect, however. If we look at the lstsq documentation, we see that this command actually returns quite a bit more than just the least squares parameters we’re interested in.
The values returned by the least squares function
To put things in a little bit nicer form for us, let’s use this instead:
phi,residues,rank,s = linalg.lstsq(A,accel)
Now we can extract the second-order parameters:
omegan = sqrt(-phi)
zeta = phi/(-2*omegan)
And now for the fun stuff, we get to define a function to run through an ode solver. Python is a full-fledged programming language, so defining functions is really quite simple, we can define a function anywhere we like, no need for a separate file or anything silly like that. The key things to defining a function are first, we use the
def command to let Python know we’re defining a function, then we name it, include the values it can take in parenthesis, then we indent all the following lines that are part of our function (indentation is important in Python). We should also
return some value before we finish our function. Here’s what we’d do to have a function that acts like our second order system:
return[ x, (-2*omegan*zeta*x - omegan**2*x)]
And that’s how easy functions are in Python! Now on to solving ODE’s. The first thing our ODE will need is the initial conditions, which are easy enough:
x0 = [pos,vel]
Now, technically to get the same algorithm that MATLAB’s ode45 uses (or in other words, to use the fourth-order Runge-Kutta method), we need to use the function
scipy.integrate.ode. This is a very advanced ode solving function that will allow us to choose one of several methods to solve an ode, unfortunately I still haven’t figured out how to use it, so in the mean time we can use the much simpler
scipy.integrate.ode. To make typing easier, we can import the library like this:
from scipy import eye
from scipy import integrate.odeint
This allows us to call the command
odeint without needing to type the
scipy.integrate at the beginning. Now, using the command is as simple as:
pos_calc = odeint(stage1, x0, time)
If you are family with plotting in MATLAB, using pylab plotting functions in Python will look very familiar to you. I imported pylab without any *, so I need to prefix all my pylab commands with
pylab.plot(time,pos, label='recorded data')
pylab.plot(time,pos_calc[:,0], label='simulated ode')
And voila, we’ve completed stage one and we have a nice graph to prove it! Just a few notes on the side, if you typed all this into a script and see no graph after you hit “run”, don’t be afraid, there is one more pylab command you need before the graph will appear:
The graphical result of our efforts thus far.
This let’s you control when in your script your graphs will appear. One thing to note, if you have a long program but you only want to execute part of it (or if you’ve written your code in “block code”) you can hit F9 and only the selected code or block will run. At least on my system, the default iPython console does this very slowly, I found it much more speedy to open up a normal Python console (Run -> Open Interpreter) and then hit F9.
Files used in this tutorial: files_used