PID controller in Haskell

Or, when will the screaming stop? When?

Lately I’ve been interested in learning how PID controllers work. To that end:

Learning is the product of the activity of learners.
John Holt

This post is actually written in literate Haskell, so you can download the source here and compile it as-is. You’ll need to ensure you have installed the latest tubes package (cabal install tubes should suffice).

1 PID controller theory

We are trying to control something - say, a rotor on a quadcoptor. We have a measurement of how level we are along some axis, along with a desired value.

The error is the difference between my measured value and my desired value. PID controllers try to minimize this error by emitting a correction value.

The error function is very simple: e(t) = desired_value - measured_value(t), where measured_value is a function giving the measurement at any given time, and t is time. This is a little pedantic but we’ll use this in a moment.

The control algorithm makes three considerations:

  • The correction should be in proportion to the size of the error (small errors should beget small corrections, etc);

  • The integral of the error function from start to the current time which, informally, tracks the amount of error we have accumulated over the runtime of the controller; and

  • The correction should consider the derivative of the error function at the current time which gives a forecast about the general direction the function is heading.

Hence, PID. Each of these three computed values is multiplied by a different constant, called a gain, allowing PID controllers to be tuned to correct behavior. In our example we are measuring how level a rotor is but the thing we are controlling is power sent to the motor. So we must have these configurable gain values.

The function that governs a PID controller is this:

u(t) = K[p] e(t) + K[i] (Integral 0 -> t of e(t)dt) + K[d] de(t)/dt

where
    K[p] = Proportional gain
    K[i] = Integral gain
    K[d] = Derivative gain
    e(t) = Error function: desired value - measured value at time t

I can feel the excitement welling within you!

2 Setup and a little background

First, the modules and language extensions I’ll be using:

{-# LANGUAGE Arrows #-}
import Tubes
import Prelude hiding (map)
import Control.Arrow

I hear you ask, “What’s the arrow crap?” There are a lot of excellent resources for arrows and I won’t bother trying to retread. The simple answer is that an arrow models a process that transforms some value a into a value b.

And again I hear you: “Isn’t that what a function does?” Indeed. And actually functions are arrows. But there are other kinds. Arrows can, for instance, perform multiple computations simultaneously, built by combining other Arrows. You can use the functions in Control.Arrow to do this or you can be lazy and have the Arrows language extension use them for you.

The tubes library defines a base Tube type, which is a computation that can suspend itself to await upstream values or yield values downstream. A Channel is a restricted, specialized variety of Tube that can do both. It also happens to be an Arrow. The tubes documentation explains this in greater detail.

3 Let’s write some actual fucking code

3.1 Integrals

First we will define a function computes a running integral. The integral term allows us to take into account how much error we have accumulated over time. If significant error persists for a while, the I term will change proportionally to reflect that.

integral :: (Fractional a, Monad m) => Channel m (a,a) a
integral = Channel $ loop 0 where
    loop sumErr = do
        (dt, err) <- await
        let result = sumErr + err*dt
        yield result
        loop result

The input is a pair of values. The first item in the pair is the number of timesteps since our last sample. The second item is the actual measured value.

We will take into account the time difference because there are a number of reasons why we might miss a few samples in a real time, complex system. Hence we scale the error by the number of timesteps since the last reading to approximate the integral.

3.2 Derivatives

And now for the derivative:

deriv :: (Fractional a, Monad m) => Channel m (a,a) a
deriv = Channel $ loop 0 where
    loop lastErr = do
        (dt, err) <- await
        yield $ (err - lastErr) / dt
        loop err

Again the input is a pair of the time since we last saw a value, and a new actual value. The derivative of a function is basically a measurement of the slope of the curve at a given point. Our algorithm is crude but will get the job done: we subtract the last value from the current value and divide it by the amount of time between the two.

4 And now the point of our arrow

Our fake readings will be a sequence of pairs: the first item in the pair will be a number indicating a timestamp, and the second will be the actual value read. Intuitively this is more or less the shape of the data I could expect from a real system.

However, we needed the time differences for our intrepid integral and derivative functions. So let’s write a Channel that turns the timestamps into time differentials:

timeDiff :: (Fractional a, Monad m) => a -> Channel m (a, a) (a, a)
timeDiff startTime = Channel $ loop startTime where
    loop lastTime = do
    (t, v) <- await
    let dt = t - lastTime
    yield (dt, v)
    loop t

Finally, we can write out PID controller using the functions we’ve written and the Arrows language extension:

pid :: (Fractional a, Monad m)
    => a -- ^ proportional gain
    -> a -- ^ integral gain
    -> a -- ^ derivative gain
    -> a -- ^ desired value
    -> Channel m (a, a) a

pid kp ki kd desired = timeDiff 0 >>> proc pv -> do
    let pv' = fmap (desired -) pv
    i <- integral -< pv'
    d <- deriv -< pv'
    returnA -< kp*(snd pv') + ki*i + kd*d

The pv term means “process variable” and is a pair containing the time differential and the measured value.

We use the (>>>) function from Control.Arrow to define pid as accepting the output of timeDiff. The Arrows syntax further makes it very simple and intuitive to feed a value through two concurrent Channels and combine their results at the end.

Let’s test it out.

main :: IO ()
main = do
    let readings = [(1, 5) , (3 , 7), (4 , 9), (7 , 14)
                   ,(8, 11), (10, 9), (11, 8), (13,12), (15, 9)]
    let target_value = 10 -- why not?
    runTube $ each readings
           >< tune (pid 0.6 0.1 0.15 target_value)
           >< map show
           >< pour display

The output from running this programThe functions tune and pour correspond to the types Channel and Sink, respectively. Both are wrappers around a more fundamental Tube type, and after they have been constructed safely they can be unwrapped using those functions. display is a Sink that comes with the tubes package.

:

4.25
2.75
1.5000000000000002
-2.65
-0.25
0.85
1.65
-1.6
0.9249999999999999

Not too shabby, honestly. While this could use some fine-tuning the correction values are more or less solid attempts to keep the output near 10.

My goal would be to write control software in Haskell using tubes, with sensors emitting streams of values and control routines operating on them in real time and constant memory usage.

Anyway I had fun and I hope you did too!