State Feedback Simulation with Functional Language

Most power supply designers have never used a functional language for simulation, but then most programers don’t even know what a functional language is. So why would a hardware designer go where even the programmers won’t go? Quite simply because a functional language is more expressive and easier to use than conventional languages for simulation.

In this blog I will code the simulation that was used in my state feedback post

To keep things simple, the example uses floating point rather than fixed point, and has very little abstraction. A proficient functional programmer will build flexible and reusable simulation frameworks, but that will confuse the beginner, so I start with simple. My goal is to help you imagine the possible, not to provide a complete simulator. But like my other posts, you can duplicate what I have done and do real work.

The language I will use is Haskell. The reason for Haskell is that it is a pure functional language, unlike F#, and it supports lazy evaluation. Lazy evaluation allows simulation of infinite lists of states, which is not used in this example, but it makes it easier to code arbitrary time sequences. Therefore, if you bother to learn Haskell for your own simulations, your investment will pay off in the long run. Furthermore, Haskell is cross platform and can run on a Mac, Linux, Windows, etc, whereas F# mainly runs on Windows. For now you must trust that this is the right functional language.

Prelims

Note that the code in this post supports scrolling to see the end of long lines. Even code without line numbers has this behavior.

A Walk Through Simulation Code

If you want an introduction to Haskell, you can search the internet. Rather than regurtigate what others have already done, I will dive right into an example, explaining what the code does, and show the results. This example with run on the GHC Haskell Platform.

Here are a link to a (Haskell Introduction) if you want to dig deeper. You may want to dive into my code until you feel a little lost, then read the introduction, then come back and reread this post. This is not an easy post, but I think the value of functional languages is great enough that you will make some investment. You can also leave comments or email if you need some help. The main thing is that if you find this post a bit challenging, it is to be expected, so don’t feel bad. Few people have experience with functional languages and your brain is trained to think in terms of state, especially if you are a hardware/electronics engineer.

As in the C# example in my previous post, there is a state model of the plant and a compensator using integral control. The simulation produces an output file that can be plotted, in this case, using gnuplot.

Design Constants

vAmpGain =      1024.0 / 3.3            :: Float
iAmpGain =      0.3173 * 1024.0 / 3.3   :: Float
pwmGain =       1.0 / 400               :: Float
reference =     321.0                   :: Float
inputVoltage =  12.0                    :: Voltage
gain =          26.0                    :: Float
t =             4.0E-6                  :: Float
k1 =            380.1/vAmpGain          :: Float
k2 =            1089.2/vAmpGain         :: Float
k3 =            10.8/iAmpGain           :: Float

We start with some constants used for the simulation. The goal is to keep the simulation similar to the real design. For example, the voltage amplifier gain is 1024 codes divided by the voltage range. The PWM gain is in volts per code. The reference is the numerical value used in the dsPIC code. The constants k1-k3 are the state feedback constants. Don’t get too hung up on the math, the goal is to learn how to structure a simulator, not learn how my particular simulator works.

Defining the Data Model

type Voltage = Float
type Current = Float
type Length = Integer
type Time = Float

Now we need some basic types. We could code this with Float and Integer, but the code will be easier to read if we have types like Voltage and Current. In Haskell, they are equivalent and the compiler will compile with either.

data State = State Voltage Current
    deriving (Show)
data Integrator = Integrator Float
    deriving (Show)
data Input = Input Float
    deriving (Show)
data Mode = Normal | Load

These are called Algebraic Types. This syntax is not so obvious, so we need to pick it apart. Let’s look at the first one. On the left side we have “data” to indicate it is a type, then “State”. State is the name of the type. On the right side, we have “State” and two values.

You can create an instance of this type as follows:

(State 1.0 2.0)

The “State” on the right side can be any name, but in this usage it is clearer to use the same name. In other circumstances names may be different. For example:

data Shape = Circle Float | Rectangle Float Float

There are two possible instances:

Circle 1.0
Rectangle 2.0 3.0

This is similar to an enumerated type with data attached to each item. A more typical enumerated type would have no data”

data Boolean = True | False

In our case we have four types: State, Integrator, Input, and Mode. Using these types will make the code more readable. We also have “deriving (Show).” This just allows the type to be printed. You can learn about that later.

Modeling the Plant

State is used for the state of the plant. The integrator is part of the controller. Input is the reference value. Mode is used to control a load transient. Mode works like an enum and has no value for each item.

plant :: Mode -> Input -> State -> State
plant Normal (Input u) (State v i) = 
    (State (0.9843*v + 0.01160*i + 0.001133*u) (-2.204*v + 0.9402*i + 0.1878*u))
plant Load   (Input u) (State v i) = 
    (State (0.9530*v + 0.01124*i + 0.001103*u) (-2.135*v + 0.9409*i + 0.1878*u))

Now we define a function that models the plant. The first line defines the type. The name of the function is “plant.Text to the right of" is the type. This function has three input parameters: Mode, Input, and State. And there is one output parameter State. In other languages it might look like:

State plant (Mode m, Input i, State s)

Why the funny syntax where the arrow is used for input and output parameters? In a functional language you can partially apply a function. So we could conceptualize the function like this:

plant::Mode->(function::Input->State->State)

The function plant, with one input of type Mode returns another function with two inputs. You can then define function as:

function::Input->(function2::State->State)

The function2 takes a State and returns a State.

So evaluation proceeds by applying parameters one by one, where each application returns a function until the last parameter is reached where a value is returned. This function could be written like this:

plant::(Mode, Input, State)->State

In this case, the function behaves more like a traditional language in that all three input parameters must be provided at the same time. There is no way to do a partial application.

When a partial application is used, the returned value is a function. Therefore in Haskell, a function can be a value.

After the type there are two lines that implement the function. This is using pattern matching. Notice that one uses Normal, and one Load:

plant Normal...
plant Load...

At run time when the function is evaluated, it will match on the Mode parameter and evaluate one of the two definitions. It will also match on the Input and State types in such a way that the values are assigned to the variables u, v, and i. These are used to calculate the next State. Basically the stuff to the left of the equals is matching, and the stuff to the right of the equals is generating the return value.

Now, let’s step back a moment. We are modeling state. If we were using a procedural language like C, we would hold state in a variable. A variable will map to a register in a microprocessor or to memory. As the state changes, the value in the register or memory changes. Haskell is a language without state. When a function is evaluated, it applies a value in a calculation, and returns a new value. It is often said that functional languages have no side effects. That refers to this aspect of values, not states. There are benefits to a language without side effects and you can research those later. Most of the benefits are tied to robustness of code, ability to prove correctness of code, etc. But for the purpose of the simulation example, just keep in mind that state is managed by consuming states and producing new states, and since these states are added to a list, they can be plotted. In a procedural language, you would have one variable that maintained state, and an array of copies of the value in the variable. Once you have a fuller understanding of functional languages, you can come to your own conclusions about this property of Haskell.

Modeling the Compensator

compensator :: Input -> Integrator -> State -> (Integrator, Float)
compensator (Input u)  (Integrator integ) (State v i) = 
    (Integrator integrator, integrator*k1 - v*k2 - i*k3)
    where
        integrator = integ+u

Now we need a compensator. This takes an Input, an Integrator, the State from the plant, and returns an Integrator and Float. The return values are in parens like (Integrator,Float). This is called a tuple. Think of it as a value with two items in it, or a struct without a name.

The compensator here is using state feedback. It first integrates as seen in the where clause, and then multiplies the integrator and state variables (v, i) by the constants k1-k3.

Modeling the System

system :: Integrator -> State -> [(Time,Integrator,State)]
system integrator state = doSystem 0 input mode integrator state
    where
        doSystem :: Integer -> [Float] -> [Mode] -> Integrator -> State -> [(Time,Integrator,State)]
        doSystem step (ref:refs) (md:mds) integrator state@(State v i) =
            ((fromInteger step * t, integ, result):(doSystem (step+1) refs mds integ result))
            where
                error = ref - plantOutput state * vAmpGain
                (integ, control) =
                    compensator (Input error) integrator (State (v*vAmpGain) (i*iAmpGain))
                pdc = control * gain
                result = plant md (Input (inputVoltage*pwmGain*pdc)) state

Now we glue it all together into a system. The system function takes an integrator, state, and returns a list of time, integrator, and state. These three return values will be plotted. The function has a supporting function named doSystem in the where clause. The doSystem function is a recursive function that evaluates the system at each step. The value named step is the time step that is increased by one one each recursive call. Let’s break down the function line by line.

doSystem :: Integer -> [Float] -> [Mode] -> Integrator -> State -> [(Time,Integrator,State)]
system integrator state = doSystem 0 input mode integrator state

The first line defines the signature of the types using the non-tuple form of a function. The second line defines the function. To the left of equals is the function name followed by two named input values, integrator and state. To the right is the function definition, which says system is a function doSystem that takes a starting step of 0, an input, mode, integrator, and a state. These are defined in the where clause.

doSystem :: Integer -> [Float] -> [Mode] -> Integrator -> State -> [(Time,Integrator,State)]
        doSystem step (ref:refs) (md:mds) integrator state@(State v i) =

The function doSystem has its own signature, but the second line is structured for pattern matching. The first parameter step will receive the initial value of 0, and then will be be evaluated with increments, like 0,1,2,3… The second parameter is written in two pieces:

(ref:refs)

The colon separates the list, defined as [Float], into to pieces. ref is the head of the list, and single float value. refs is the tail of the list, as list of the remaining values. This allows the function to use the first element in the list and then pass the rest in a recursive call to doSystem to evaluate the next step in the simulation.

The pattern match is the state. It is written:

state@(State v i)

The tag to the left of the refers to the whole value. The piece to the right of the breaks the state into its constituent pieces, v and i. This allows the function to operate on the pieces, or on the whole.

error = ref - plantOutput state * vAmpGain (integ, control) =
compensator (Input error) integrator (State (v*vAmpGain) (i*iAmpGain))
pdc = control * gain
result = plant md (Input (inputVoltage*pwmGain*pdc)) state

Before we look at the main doSystem definition, we need to look at the where clause. The error line calculates the error that is fed into the plant. The compensator is fed the error, the integrator, the state, and returns an integrator and control value. The parens in (inter, control) take the return value as a tuple. The state is constructed from the v and i of the previous state using gains, just like the input is constructed by the error. Input and State are types. We could pass floats, but then the code will not be as descriptive.

Eventually the result value is calculated from the plant, which applies the state feedback to the state. This leads to the main definition.

doSystem step (ref:refs) (md:mds) integrator state@(State v i) =
            ((fromInteger step * t, integ, result):(doSystem (step+1) refs mds integ result))

This builds the return tuple followed by the recursive call. Let’s break this down. The tuple is:

((fromInteger step * t, integ, result):

The step times the t value is the time, then the integrator value, then the result in the form of the new state. The : divides the head, on the left in the form of the tuple, from the tail which comes from the recursive call:

(doSystem (step+1) refs mds integ result))

The tail calls doSystem with an incremented step, and passes the tail of its input values, the integrator, and the current state.

Most functional languages handle list processing with this head : tail structure, which is also found in Lisp/Scheme. It is used in matching patterns to dissect a list, and in constructing a list. One thing that makes Haskell unique is that it is a lazy language. This means values are not generated until needed. This allows you to write functions that generate infinite lists of values. This is useful for building simulators, because you can control the simulation time independently from the definition of the plant, compensator, etc.

Simulation

Now that we have a system, we need to simulate:

simulate :: Integer -> [(Time,Integrator,State)]
simulate len = doSimulate 0 len (system compensatorInit plantInit)
    where
        doSimulate :: Integer -> Integer -> [(Time,Integrator,State)] -> [(Time,Integrator,State)]
        doSimulate start finish (out:outs)
            | start == finish = []
            | otherwise = out:doSimulate (start+1) finish outs

Simulate takes an integer that defines the length of the simulation. It is defined as doSimulate so that you can hide the value that tracks time.

doSimulate 0 len (system compensatorInit plantInit)

The doSimulate takes a 0 and length. The 0 value is incremented on the recursive call. This hides the counting variable from the user of the function simulate. This call uses a language feature called guards (|).

doSimulate start finish (out:outs)
            | start == finish = []
            | otherwise = out:doSimulate (start+1) finish outs

Let’s unpack this a bit. We give doSimulate a start and finish value. We give a list which comes from the function system, and we peel off the head into the value named out, and the tail into the value named outs. When start is equal to finish, we give an empty list as the result. For all other values we evaluate doSimulate with an incremented start value, passing the new start value and the tail outs. This recursive call recurses until start and finish are equal. When the recursion unwinds after the last value, the output list is constructed by:

out:doSimulate

where the colon divides the new value on the left, and the list it is hooked to it on the right. This may seem inefficient compared to the loop. However, good compilers can recognize this pattern and eliminate the overhead from the evaluation such that when start equals finish, it can build the list directly and there is no need to unwind the recursive definition. This allows you to define the function recursively, but the actual implementation is still efficient.

Modeling System Inputs

We have to feed the system with an input, so let’s look at how we can define the input as in infinite list with a step.

input :: [Float]
input = doInput 0
    where
        doInput :: Integer -> [Float]
        doInput step
            | step > 100 = reference:doInput (step+1)
            | step > 50 = 1.1*reference:doInput (step+1)
            | otherwise = reference:doInput (step+1)

This definition uses a recursive doInput with a value that tracks the number of steps. The guards are used to control the step. Values past step 100 use the value reference, which is one of the constants defined at the beginning of the code. Between 50 and 100 the reference is 10% higher. All other values are the reference. Notice there is no guard like this one:

| start == finish = []

This means there is no natural end to the list. The list is infinite. However, the simulate function limits the number of steps, so the actual list generated by this infinite list definition is limited to the simulation size. Infinite lists allow us to define the input without any model of the simulation length. The model of the input is decoupled from the simulation length.

This can be done in non-functional languages with other mechanisms, but other functional languages such as ML don’t allow infinite lists. ML is an eager language, so if you write code like this, evaluation of the input function would generate an infinite list and the program would hang. This is the main reason I recommend Haskell as a simulation language over ML.

The Code

Now, I have left out some of the details of modules, writing to files, etc. However, here is the complete code so you can run this yourself with GHC.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155

-----------------------------------------------------------------------------
--
-- Module      :  Main
-- Copyright   :  Proclivis Product Development 2012
-- License     :  GNU
--
-- Maintainer  :  Mike Jones
-- Stability   :  Stable
-- Portability :
--
-- A simple simulation example for a buck
--
-----------------------------------------------------------------------------
 
module Main (
    main
) where
 
import System.IO
 
-- Design constants
 
vAmpGain =      1024.0 / 3.3            :: Float
iAmpGain =      0.3173 * 1024.0 / 3.3   :: Float
pwmGain =       1.0 / 400               :: Float
reference =     321.0                   :: Float
inputVoltage =  12.0                    :: Voltage
gain =          26.0                    :: Float
t =             4.0E-6                  :: Float
k1 =            380.1/vAmpGain          :: Float
k2 =            1089.2/vAmpGain         :: Float
k3 =            10.8/iAmpGain           :: Float
 
-- Basic types
 
type Voltage = Float
type Current = Float
type Length = Integer
type Time = Float
 
-- Algebraic types for simulation
 
data State = State Voltage Current
    deriving (Show)
data Integrator = Integrator Float
    deriving (Show)
data Input = Input Float
    deriving (Show)
 
-- Algebraic types for controlling the load
 
data Mode = Normal | Load
 
-- The plant is a state equation for a simple Buck
 
-- Create a plant
plantInit :: State
plantInit = (State 0.0 0.0)
 
-- Increment one step in the plant simulation
plant :: Mode -> Input -> State -> State
plant Normal (Input u) (State v i) =
    (State (0.9843*v + 0.01160*i + 0.001133*u) (-2.204*v + 0.9402*i + 0.1878*u))
plant Load   (Input u) (State v i) =
    (State (0.9530*v + 0.01124*i + 0.001103*u) (-2.135*v + 0.9409*i + 0.1878*u))
 
-- Get the output of the plant
plantOutput :: State -> Voltage
plantOutput (State v i) = v
 
-- The compensator is a state feedback with integral control
 
-- Create a compensator
compensatorInit :: Integrator
compensatorInit = Integrator 0.0
 
-- Increment one step in the compensation simulation
compensator :: Input -> Integrator -> State -> (Integrator, Float)
compensator (Input u)  (Integrator integ) (State v i) =
    (Integrator integrator, integrator*k1 - v*k2 - i*k3)
    where
        integrator = integ+u
 
-- Input
input :: [Float]
input = doInput 0
    where
        doInput :: Integer -> [Float]
        doInput step
            | step > 100 = reference:doInput (step+1)
            | step > 50 = 1.1*reference:doInput (step+1)
            | otherwise = reference:doInput (step+1)
            
-- Mode
mode :: [Mode]
mode = doMode 0
    where
        doMode :: Integer -> [Mode]
        doMode step
            | step > 200 = Normal:doMode (step+1)
            | step > 150 = Load:doMode (step+1)
            | otherwise = Normal:doMode (step+1)
 
-- The system             
system :: Integrator -> State -> [(Time,Integrator,State)]
system integrator state = doSystem 0 input mode integrator state
    where
        doSystem :: Integer -> [Float] -> [Mode] -> Integrator -> State -> [(Time,Integrator,State)]
        doSystem step (ref:refs) (md:mds) integrator state@(State v i) =
            ((fromInteger step * t, integ, result):(doSystem (step+1) refs mds integ result))
            where
                error = ref - plantOutput state * vAmpGain
                (integ, control) =
                    compensator (Input error) integrator (State (v*vAmpGain) (i*iAmpGain))
                pdc = control * gain
                result = plant md (Input (inputVoltage*pwmGain*pdc)) state      
 
-- Simulate 
simulate :: Integer -> [(Time,Integrator,State)]
simulate len = doSimulate 0 len (system compensatorInit plantInit)
    where
        doSimulate :: Integer -> Integer -> [(Time,Integrator,State)] -> [(Time,Integrator,State)]
        doSimulate start finish (out:outs)
            | start == finish = []
            | otherwise = out:doSimulate (start+1) finish outs
 
-- Formatting functions
 
-- Take a list of integration and state values and returns a list of strings
formatResults :: [(Time, Integrator,State)] -> [String]
formatResults [] = []
formatResults ((time, (Integrator integ), (State v i)):iss) =
    (show time ++ "\t" ++ show integ ++ "\t" ++ show v ++ "\t" ++ show i) : (formatResults iss)
 
-- Print a list of strings to stdio
printResults :: [String] -> IO ()
printResults [] = return ()
printResults (s:ss) =
    do
        putStrLn s
        printResults ss
 
-- Write a list of strings to a file handle
fileResults :: Handle -> [String] -> IO ()
fileResults h [] = return ()
fileResults h (s:ss) =
    do
        hPutStrLn h s
        fileResults h ss
        
-- The module entry point
main = do
    h <- openFile "/Users/michaeljones/Documents/data.csv" WriteMode
    fileResults h (formatResults (simulate 250))
    hClose h

Graphing

The above simulation creates a CSV file. You can plot this with gnuplot or your favorite plotting tool. For gnuplot:

set xlabel 'Time'
plot '/Users/michaeljones/Documents/data.csv' using 1:3 with lines title 'Buck Output'

p.

Figure 1  Simulation Results

This shows the first state variable v. The first bump is a reference change. The second two transients are a load change.

Summing Up Simulation Results

My goal for this post was to stimulate ideas. A lazy functional language has attributes that aid writing simulation code, such as recursion and lazy evaluation. You can take this code and play with it using the GHC compiler.

However, an industrial strength simulator would not be coded so directly. For my own simulations I have a set of reusable functions so that I can define a particular simulation with just a few lines of code for the constants. I have avoided that here and presented a direct and concrete simulation because it is simpler to understand how the language works, yet you can copy it and do real work and decide for yourself if you can get value out of using a functional language.

If you prefer a language with commercial support and are willing to be tied to one vendor/architecture, you might take a look at F#. F# is not a pure functional language, so you will have to learn how to mix two programing paradigms, however there are benefits to this and you may also be familiar with Visual Studio.

Now for Some Subtleties

It turns out that the simulation does not produce optimal results. I discovered this quite accidentally when coding a simulation framework to generalize the simulator. The same issue applies to the simulation in C# and the dsPIC code used in my other posts. So if you made it this far, you get to see the fix. Let’s begin with the better results compared to the original results.

Figure 2  Simulation Results

The black line is the original simulation, and the blue line is the corrected simulation.

The change is shown in this code:

compensator :: Input -> Integrator -> State -> (Integrator, Float)
compensator (Input u)  (Integrator integ) (State v i) =
    (Integrator integrator, integ*k1 - v*k2 - i*k3)
    where
        integrator = integ+u

Notice the difference in these two lines:

(Integrator integrator, integ*k1 - v*k2 - i*k3)      NEW
 (Integrator integrator, integrator*k1 - v*k2 - i*k3)     OLD

All that is different is that the k1 constant is multiplied with the old integrator value instead of the updated one. Practically speaking, this means that during a reference change or startup, as the output approaches the reference, the integrator value will be a bit larger because the previous value has more error. If you go back and review the equations in Figure 9 of the state feedback post, you will see that this correction matches the equations. My bad :-)

A Simulation Framework

It is always easier to use a framework rather than recode concrete versions of simulators. So in this section I’ll show you what my code looks like when using my floating point simulator. We a bit of learning, you can build your own simulator in Haskell.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43

module Main (
    main
) where
 
import System.IO
import Input
import State
import Plant
import Mode
import Integrator
import Reference
import Compensator
import System
import Support
 
-- Design constants
 
vAmpGain    = 1024.0 / 3.3            :: Float
iAmpGain    = 0.3173 * 1024.0 / 3.3   :: Float
pwmGain     = 1.0 / 400               :: Float
ref         = 321.0                   :: Float
inputVoltage= 12.0                    :: Float
gain        = 26.0                    :: Float
dt          = 4.0E-6                  :: Float
k1          = 380.1/vAmpGain          :: Float
k2          = 1089.2/vAmpGain         :: Float
k3          = 10.8/iAmpGain           :: Float
 
input       = initialInput 50 100 inputVoltage (1.1*inputVoltage)
reference   = initialReference 150 200 ref (1.1*ref)
mode        = initialMode 250 300
plant       = initialPlant 0.9843 0.01160 0.001133 (-2.204) 0.9402 0.1878
                0.9530 0.01124 0.001103 (-2.135) 0.9409 0.1878
integrator  = initialIntegrator
state       = initialState
compensator = initialCompensator vAmpGain iAmpGain k1 k2 k3
system      = initialSystem input reference mode plant integrator compensator vAmpGain iAmpGain pwmGain gain dt
 
-- The module entry point
main = do
    h <- openFile "/Users/michaeljones/Documents/data.csv" WriteMode
    fileResults h (formatResults (simulate system 350))
    hClose h

Figure 3  Simulation Results

This simulator supports changes in the input, reference, and load. You can see the input change at 50 * dt, the reference change at 150 * dt, and the load change at 250 * dt. As in the original example, the simulation time is decoupled from the models by modeling with infinite lists.

Finishing Off

I have shown how to build a simple Haskell simulator, and what a reusable one looks like. Hopefully this stimulates some interest in functional programming. Let me know what you think.