Digital Compensators using State Feedback Techniques

My last post demonstrated discrete root locus design. This post will demonstrate pole placement using state space methods. The end result will be a 4X faster transient response to a load disturbance. As in past posts, all the MATLAB, simulation, and dsPIC code will be shown.

Figure 14  Root Locus vs. State Space

Here is a preview of the results. 4X better settling, 1/2 the voltage drop.

The root locus design from the last post added an integrator, a pair of zeros, and a pair of poles. The difficulty of the technique lies in placing them to get an adequate response, but also dealing with the complexity. A simple second order system is easy to understand. Set the damping ratio and frequency and your done. But once you add the integrator, two zeros, and two poles to the two poles of the plant, you have 5 poles and 2 zeros! Unless the zeros cancel poles, which they did not, it is hard to predict the outcome other than trial and error. You can’t treat the response as a second order system, unless one pair of poles dominates the response.

So what is the alternative? It is to move the poles you have, rather than add more poles and zeros. Sound to good to be true? Well, let’s just do it and see.

Overview

Here is where we are going:

1. Writing the differential equations
2. Writing the state space equations
3. Simulating the plant
4. Structuring the feedback
5. Choosing the pole locations
6. Simulating the transient response
7. Simulating the finite math transient response
8. Implementing the dsPIC code
9. Comparing results with root locus
10. Wrap up

Writing the differential equations

The first step is to write out the differential equations for our buck converter. There are two configurations of the switch as shown here:

Figure 1  Buck Topology

In the upper position, the input is connected to the inductor. In the second position the inductor is grounded.

Figure 2  Differential Equations with Upper Switch Position

These are the equations when the input is connected to the inductor. V is the voltage on the capacitor, and I is the inductor current. When the switch is in the lower position, Vi will be zero, so the equations are exactly the same other than Vi is set to zero.

Writing the state space equations

Now the equations are put into state space form. I am using the implicit form, with the E on the left side next to dx/dt. I use this form because it is simpler to write out the state space from the differential equations.

Figure 3  State Space with Upper Switch Position

The basic idea is to treat the energy storage devices as states. So the voltage on the capacitor is one state, and the current in the inductor is the other state.

Simulating the plant

Once the equations are in the state space form, we need to create a model by averaging the state space equations for the two switch positions, using the duty cycle. All this can be done in MATLAB.

a1 = [(-1/(Rload+Resr)) (Rload/(Rload+Resr)); 
(-Rload/(Rload+Resr)) -(Resr*Rload/(Rload+Resr)+Rdsr)];
b1 = [0; 1];
c1 = [1 0; (Rload/(Rload+Resr)) (Rload*Resr/(Rload+Resr))];
d1 = [0; 0];
a2 = [(-1/(Rload+Resr)) (Rload/(Rload+Resr)); 
(-Rload/(Rload+Resr)) -(Resr*Rload/(Rload+Resr)+Rdsr)];
b2 = [0; 0];
c2 = [1 0; (Rload/(Rload+Resr)) (Rload*Resr/(Rload+Resr))];
d2 = [0; 0];
a = dc*a1 + dcb*a2;
b = dc*b1 + dcb*b2;
c = dc*c1 + dcb*c2;
d = dc*d1 + dcb*d2;
e = [C 0; 0 L];

plant = dss (a, b, c, d, e);
dc = Duty Cycle
dcb = 1 - Duty Cycle

We now have a model for the plant and we can do a bode plot to see if things look right. (Note, if you have never used MATLAB, you have to do a substitution of the variables. I am trying to prevent clutter to keep the post simple. Use comments if you get stuck and need help.)

Figure 4  Bode Plot of Plant

This shows the plant with two loads, a light load, and a heavy load. The heavy load has more peak. there are also two outputs, output 1 is voltage, and output 2 is current. The basic thing to notice is we are dealing with a classic two pole system. I am also neglecting the third output, the voltage at the load. The voltage here is the voltage on the capacitor.

We can also model this with a discrete state space.

plant = c2d(plant, T);

Figure 5  Bode Plot with Discrete State Space

In this plot I have added the discrete version in green. You can see how the phase delay from sampling rolls off the phase. The variable T used in c2d was 5uS. The discrete version will be used to design the compensator and simulate the system.

Structuring the feedback

The compensator structure uses Integral Control. You can read about this in “Digital Control of Dynamic Systems,” by Gene Franklin et al, Chapter 6, Section 6.

Figure 6  Compensation Structure

This structure comes right out of Franklin. The inner loop has K feedback. The outer loop has an integrator. The design I am presenting does not use N or H. The w input is the disturbance input. I call it windage so I can remember it. The Plant with the two greek letters are the A and B of the state space model. The C and D of the state space model are outputs generated from the states x.

Let’s first talk about K. The K box are constants multiplied by the states, in our case V and I. These constants move the poles of the plant. However, the resulting gain is small, so the output will not reach the r input, which is the regulator input. The integrator is what makes the output y achieve the desired value r.

You can see that the only thing added to the response of the plant is the integrator. No new poles or zeros. This will allow us to achieve a better transient response.

Choosing the pole locations

We have to decide where we want the poles.

Figure 7  Desired Poles

Let’s choose the poles shown here. You can see a pole at 0, that is the integrator. Then the two pole positions we want to move to. This would result in a frequency response like this:

Figure 8  Desired Frequency Response

DesiredResponseZ = zpk([], [0.2+0.15j 0.2-0.15j 0], 0.5, T);
[poles, zeros] = pzmap(DesiredResponseZ);

A nice two pole roll off with a 3dB point around 20KHz. This gives us a lot of bandwidth to play with, but we are not too close to the clock frequency of 200KHz. We have to have space to filter out the clock, etc.

To move the poles, we need the find the K values. However, we have to move the poles and deal with the integrator. So we have to manipulate the state equations by adding the integrator. Franklin gives the equation.

Figure 9  Modified State Equation

So let’s build this in MATLAB.

ab = [[1; 0; 0] [[plant.c(1) plant.c(4)]; plant.a]];
bb = [0; plant.b];
cb = [0 plant.c(1) plant.c(4)];
db = [plant.d(2)];
plant2 = ss(ab, bb, cb, db, T);

It might take a bit of work to figure out how that was done. Basically, it carves out the pieces from the plant and adds the pieces for the integrator. Once this far, we can build the compensator K values.

compensator = acker(plant2.a, plant2.b, poles);

This uses Ackermann’s Formula which is well documented by Franklin. The end result is:

294.8930  844.9357    8.3471

Three constants. The first one, K1, is the integrator constant. The second one, K2, is the voltage state, and the third one, K3, is the current state. We now have all the pieces in the Integral Control model and need to simulate it.

Simulating the transient response

integrator = ss(1,1,1,0,T);
smallCompensator = [compensator(2) compensator(3)];
innerLoop = feedback(plant, smallCompensator, [1], [1 2]);
outerLoop = feedback(compensator(1)*integrator*innerLoop, 1, [1], [1]);

What you should notice is that the loop is created with the original plant, not the one used for the acker function. The integrator is also introduced in the conventional way.

Figure 10  Step Response of Model

We end up with 30uS step response. Now we have to see if this is real by simulating the finite math.

Simulating the finite math transient response

The previous root locus post showed how to create the diff equations, so I am not going to do that here. However, I’ll cover some of the code. Since we need to simulate a load response, I created a plant model with two different loads. This means two sets of diff equations. I step the load by changing which set of equations is used.

First, the plant:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22

public void Calculate (double control)
        {
            u1 = control;
 
            double x1New = 0.0;
            double x2New = 0.0;
 
            if (load)
            {
                x1New = 0.953*x1 + 0.01124*x2 + 0.001103*u1;
                x2New = -2.135*x1 + 0.9409*x2 + 0.1878*u1;
            }
            else
            {
                 x1New = 0.9843 * x1 + 0.0116 * x2 + 0.001133 * u1;
                 x2New = -2.204 * x1 + 0.9402 * x2 + 0.1878 * u1;
            }
 
            x1 = x1New;
            x2 = x2New;
            y1 = x1;
        }

This is very simple. Just use a conditional to turn the load on and off by changing diff equations.

Second, the compensator:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

public short Calculate(short x1, short x2, short u1)
        {
            qx += u1;
 
            int accA = 0;
            accA = DSP.MAC(accA, qk1, qx);
            accA = accA << 1;
            accA = DSP.MAC(accA, qk2, x1);
            accA = accA << 2;
            accA = DSP.MAC(accA, qk3, x2);
            accA = accA << 2;
 
            qy = DSP.Get(accA);
 
            X1.Add(qx);
            X2.Add(x1);
            X3.Add(x2);
 
            return qy;
        }

Now things are a bit tricky. There is some scaling going on. I’ll cover that in detail when we get to the dsPIC code, but the left shifts are multiplying up the results to match the scaling down of the constants.

Lastly, the simulation code:

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

public void SimulateFixedPoint(bool useNoise)
        {
            double vAmpGain = 1024.0/3.3;
            double iAmpGain = 0.3173 * 1024.0/3.3;
            double pwmGain = 1.0/400;
            Q15 k1 = new Q15((294.8930 / (vAmpGain)) / 8.0);
            Q15 k2 = new Q15(((844.9357 / vAmpGain) / 4.0));
            Q15 k3 = new Q15(8 * 8.3471 / iAmpGain);
            double reference = 321;
 
            Plant plant = new Plant();
            Compensator compensator = new Compensator(k1, k2, k3);
            
            for (int i = 0; i < simSize; i++)
            {
 
                if (i == 300)
                /   reference = 0x0200;
                /f (i == 600)
                /   reference = 0x0080;
                if (i == 900)
                    plant.Load = true;
                if (i == 1200)
                    plant.Load = false;
                short error = (short)(-((plant.Y1*vAmpGain-reference)*1.0));
                short control = compensator.Calculate(
                         (short)(plant.X1 * vAmpGain),
                        (short)(plant.X2 * iAmpGain), 
                        error
                        );
                double PDC = control*0.8;
                plant.Calculate(12.0 * pwmGain * PDC);
             }
        }

From here you can see the scaling of the constants. The math is done in order so that after each MAC instruction you shift left. When all done, the scaling is in the ballpark, and the final control value is multiplied by a constant, in this case 0.8, to tune the gain.

So let’s see how this turned out:

Figure 11  Fixed Point Simulation

Not quite what MATLAB said, but much better than the root locus design. I did not try to figure out why it is different. That seems like chasing your tail.

Implementing the dsPIC code

So now we have to code the real thing and see how it performs. Let’s start with the state table and constants:

1
2
3
4
5
6
7
8
9
10
11
12

StateTable[0] = 0; // u1
    StateTable[1] = 0; // x1
    StateTable[2] = 0; // x2
    StateTable[3] = 0; // x3
    StateTable[4] = 0; // y1
 
    double vAmpGain = 1024.0/3.3;
    double iAmpGain = 0.3173 * 1024.0/3.3;
    
    ConstTable[0] = Q15((294.8930 / (1.0*vAmpGain)) / 8.0);             // k1
    ConstTable[1] = Q15(-(844.9357 / vAmpGain) / 4.0);            // k2
    ConstTable[2] = Q15(-8.0 * (8.3471 / iAmpGain));            // k3

The constants are scaled down to less than one. We will operate on them on the given order. After we multiply X1 by K1, we can shift left, then multiply X2 by K2, shift, etc.

Here is the calculation of the control value:

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

void __attribute__((__interrupt__, no_auto_psv)) _ADCP3Interrupt()  
{
    _ADCP3IF = 0;                                // Clear the ADC pair3 
 
    register int regA asm("A");
        
       StateTable[0] = (RefOutputVoltage - VOUT_F);                    // u1 = deltaV
 
    regA = __builtin_lac(StateTable[1],0);
    regA = __builtin_add(regA, StateTable[0], 0);
    StateTable[1] = __builtin_sac(regA,0);                            // x1 += x1
 
    StateTable[2] = VOUT_F;                                        // v = vout
 
    StateTable[3] = IVSENS_P - IVSENS_N - 0x30;                    // i = iout+ - iout-
 
    regA = __builtin_movsac(&StateTablePtr, &xPrefetch, 2,
                            &ConstTablePtr, &yPrefetch, 2, 0, regA);
    regA = __builtin_lac(0,0);                                        // acc = 0
 
    regA = __builtin_mac(regA, xPrefetch, yPrefetch,
                         &StateTablePtr, &xPrefetch, 2,
                         &ConstTablePtr, &yPrefetch, 2, 0, regA);    // acc += k1 * x1
    regA = __builtin_sftac(regA, -1);                                // acc = acc << 1
 
    regA = __builtin_mac(regA, xPrefetch, yPrefetch,
                         &StateTablePtr, &xPrefetch, 2,
                         &ConstTablePtr, &yPrefetch, 2, 0, regA);    // acc += k2 * x2
    regA = __builtin_sftac(regA, -2);                                // acc = acc << 2
 
    regA = __builtin_mac(regA, xPrefetch, yPrefetch,
                         &StateTablePtr, &xPrefetch, -6,
                         &ConstTablePtr, &yPrefetch, -6, 0, regA);    // acc += k3 * x3
    regA = __builtin_sftac(regA, -2);                                // acc = acc << 2
 
    StateTable[4] = __builtin_sac(regA,0);                            // y1 = acc
 
    regA = __builtin_mpy(StateTable[4], Gain,
                         &StateTablePtr, &xPrefetch, 0,
                         &ConstTablePtr, &yPrefetch, 0);
    PDC1 = __builtin_sac(regA,0);
}

We can pick this apart. First, we calculate the error, then we integrate it. We use __builtin_add instruction and pull the value from the accumulator. When then get the voltage, which ultimately comes from the C, and we get the current, which ultimately comes from the L. These are the states that will be multiplied by the K’s.

Then we setup the prefetch and load 0 into the accumulator. Finally we MAC and left shift, accounting for the scaling of the constants. When all this is done, we store the value and multiply by a Gain. This is same as tuning the control in the finite math simulation.

Then we update the PWM.

You should be able to analyze the finite math simulation and PIC code to see how it all works. The result is pretty much a match with the simulation.

Figure 12  State Space Design Load Response

This matches the finite math simulation quite well. A 100mV dip that gets back to the output voltage in 60uS.

Comparing results with root locus

For comparison, let’s look at the root locus results again.

Figure 13  Root Locus Design Load Response

Seems we are about 4X better on settling time. We also have a much smaller dip. Not only that, but the root locus design clocked at 250KHz, and the state space at 200KHz. Also, the root locus had a much smaller delay from measurement to PWM change. In the state space design, all the math was done before the PWM change. Yet it is still much better.

Let’s put them side by side:

Figure 14  Root Locus vs. State Space

Notice the time scales and voltage scales are different. The state space design is on top.

Figure 15  Kelly and Rinne Results

This scope shot came from the paper “Control of DC-DC Converters by Direct Pole Placement and Adaptive Feedforward Gain Adjustment,” by kelly and Rinne 2005.

Recovery is about 60-70uS. Dip around 80mV. Their design used a 25% to 70% load change, for a 1A max output, meaning 0.5A change. My design stepped from 200mA to A.

The design in the paper did not use an integrator, but instead used feed forward. There is basically a control loop that runs slower than the PWM loop and adjusts the gain of the feed forward. I would guess that this is what prevents the overshoot in my design, but it also accounts for the 20mV drop in the output. The gain loop can’t react fast enough to correct the load change in the scope shot.

This gives a point of reference for the performance achieved by the design presented..

Wrap up

So what’s the bottom line here? Basically, moving poles gives a better response than adding poles and zeros. If your an analog designer, you have change your mind set from adding to shifting. However, all your skills about choosing where you want the poles to be applies, as long as you can transfer those skills from the s-domain to the z-domain.