I have decided to create a library for solving common control engineering problems in Python, available here on GitHub.

The goal is to be able to rapidly protype controllers in Python, in a way comparable to how one can do it in Matlab (but, without the big price tag). There is a similar project (by Richard Murray, no less) called the Python Control Systems Library.

Why re-invent the wheel? The project by Richard Murray requires Slycot, which was (for me, at least) a major pain to install. Controlpy should rely only on standard Python libraries such as NumPy as SciPy, making it much easier to install. Finally, it’s a nice opportunity for me to improve my Python skills, and to perhaps refine my knowledge on controls. Perhaps it may be useful for other people too.

What can it do? At the moment it can solve for LQR controllers and analyse the controllability of an LTI system. I hope to add H2 and H-infinity synthesis soon too. I will be adding features as I need them.

UPDATE: Please see here for an update: Python control library: controlpy

I have been using Python exclusively for my scientific computing for about half a year (having been frustrated by Matlab’s awkward syntax for complex programs, and annoying licensing). One area where I struggled was to compute steady state LQR controllers.

I’ve learnt that SciPy offers a Ricatti equation solver (`scipy.linalg.solve_continuous_are` and `scipy.linalg.solve_discrete_are`). Since solving the Ricatti equation is the hard part of solving for an LQR gain, this implies that one can compute infinite horizon LQR controllers straight-forwardly using only SciPy.linalg. Similarly, one can compute steady state Kalman filters.

Below are my wrapper functions for continuous and discrete time LQR controllers. The equations come from Bertsekas “Dynamic Programming and Optimal Control”. I’ve done some basic sanity checks, and it seems to work.

from __future__ import division, print_function

import numpy as np
import scipy.linalg

def lqr(A,B,Q,R):
"""Solve the continuous time lqr controller.

dx/dt = A x + B u

cost = integral x.T*Q*x + u.T*R*u
"""
#ref Bertsekas, p.151

#first, try to solve the ricatti equation
X = np.matrix(scipy.linalg.solve_continuous_are(A, B, Q, R))

#compute the LQR gain
K = np.matrix(scipy.linalg.inv(R)*(B.T*X))

eigVals, eigVecs = scipy.linalg.eig(A-B*K)

return K, X, eigVals

def dlqr(A,B,Q,R):
"""Solve the discrete time lqr controller.

x[k+1] = A x[k] + B u[k]

cost = sum x[k].T*Q*x[k] + u[k].T*R*u[k]
"""
#ref Bertsekas, p.151

#first, try to solve the ricatti equation
X = np.matrix(scipy.linalg.solve_discrete_are(A, B, Q, R))

#compute the LQR gain
K = np.matrix(scipy.linalg.inv(B.T*X*B+R)*(B.T*X*A))

eigVals, eigVecs = scipy.linalg.eig(A-B*K)

return K, X, eigVals

An alternative method is to use the Python Control toolbox. This offers a Matlab-like syntax for a variety of control methods. The downside is that it requires Slycot library, which is annoying to get running on Windows. It also does not appear to support discrete time systems.