Now, assuming no error, you would eventually have in Scilab a continuous-time controller from synthesis, with the name Skr and Skr_tf for state-space and transfer function data format, respectively. Depending on your weight setup and gamma value, the resulting controller may be different. To be consistent, we stick to this third order controller from Part I

(1)

Step 6 : convert to discrete-time controller

To implement a controller, its discrete-time representation is needed. So (1) must be converted to a transfer function in Z-domain. Refer to Module 7 of Scilab Control Engineering Basics, bilinear transformation is normally a preferred choice to obtain . Conveniently, Scilab has a command cls2dls to do so. A linear system in state-state format and sampling period must be specified.

```
-->Ts = 0.01; // sampling period
-->Skd = cls2dls(Skr,Ts);
-->Skdtf=ss2tf(Skd)
Skdtf =
2 3
432.0685 - 1290.3395z + 1280.5196z - 422.2483z
----------------------------------------------
2 3
- 0.7323527 + 2.425178z - 2.6928211z + z
```

For the discrete-time controller to maintain the desired properties, it is imperative that the sampling period Ts passed to cls2dls, must match the timer interrupt period used on the target MCU. For our implementation, we select 0.01 sec as sampling period. In motion control it is sometime called “servo cycle.”

Well, if you have no experience with controller implementation, this is perhaps a good time to read Module 6.

Step 7 : (optional) rearrange to DF II structure

A technique that helps in the controller implementation is known as “direct-form conversion.” This discrete-time manipulation is commonly used in digital signal processing to reduce memory usage by half. Though that benefit is not significant in our low-order controller, it does help formulate a compact algorithm structure.

From now on, DF stands for Direct Form. We are interested in 2 variations: DF I and DF II. As a simple example, for a second-order transfer function

(2)

The conversion from DF I to an intermediate form, then to DF II is illustrated in Figure 6.

Notice that the number of delays ( blocks) are reduced from 4 to 2. This manipulation can be extended to a transfer function of any order.

Back to our discrete-time controller from synthesis

(3)

Multiply both numerator and denominator by

(4)

and compare to a general third order form

(5)

we obtain values for coefficients and . Having already the discrete-time transfer function Skdtf in Scilab, we can extract coefficient values using command coeff. Note, however, that Scilab orders the numerator and denominator polynomials in reverse order, and also that Scilab array index starts from 1, so care must be taken to extract the correct coefficient values.

To anticipate an even more confusing scenario, the algorithm will eventually be written in C. And array index in C starts from 0! Anyway, do not let this issue scare you. It’s just some sort of “bookkeeping” problem.

Let’s get back to the controller (5). Its DF II structure can be constructed using standard Xcos blocks as shown in Figure 7. Here we use Xcos only for diagram display purpose. It is left as an exercise to the reader to build an Xcos feedback model for simulation.

What we are more interested is to write a control algorithm from the DF II structure in Figure 7. Note the variables ] along the middle vertical path through the delays. These are defined as controller states.

Step 8 : coding the control algorithm

With the above development, we are ready to write a C algorithm for the controller (4). Somewhere at the top of C source file, declare global arrays for the coefficients, controller states, and input and output variables as follows

```
double coeffa[3] = {-2.692821, 2.425178, -0.732353};
double coeffb[4] = {-422.248307, 1280.519630, -1290.339489, 432.068496};
double cx[4]={0,0,0,0}; // controller states
double u_0; // output variable for positive sense
double u_0n; // output variable for negative sense
double e_0; // error input for the controller
```

and put the control algorithm in a timer interrupt routine. Here our example is for PIC24EP256MC202, a 16-bit MCU that uses Microchip’s XC16 compiler as development tool. The syntax for timer 1 interrupt is

```
void __attribute__((interrupt, auto_psv)) _T1Interrupt(void)
// Timer 1 interrupts every 0.01 second
{
// controller state updates
*(cx+4)=*(cx+3);
*(cx+3)=*(cx+2);
*(cx+2)=*(cx+1);
*(cx+1)=*cx;
e_0 = pcmda - dcmpos; // error input = command – actual pos
// ******* Control algorithm in DF II structure *******
*cx = e_0 -*(coeffa)*(*(cx+1))-*(coeffa+1)*(*(cx+2))-*(coeffa+2)*(*(cx+3));
u_0 = *(coeffb)*(*(cx))+*(coeffb+1)*(*(cx+1))+*(coeffb+2)*(*(cx+2))+* (coeffb+3)*(*(cx+3));
// ** The following code converts output to PWM and DIR signals **
if (u_0>=0) { // positive sense
if (u_0 < PWMMAX)
PWMVal = (unsigned int)u_0; // limit to PWM range
else PWMVal = PWMMAX;
DIR = 0; // rotate clockwise
}
else { // negative sense
u_0n = -u_0;
if (u_0n < PWMMAX)
PWMVal = (unsigned int)u_0n; // limit to PWM range
else PWMVal = PWMMAX;
DIR = 1; // rotate counter-clockwise
}
OC1R = PWMVal; // load to register of output compare 1 module
_T1IF = 0; // clear the interrupt flag
}
```

Some points in the C algorithm above need further explanation

- The following variables and macros also need to be declared at the top of source file
- pcmda: (double) angle command reference, from user input or generated from cubic spline function
- dcmpos: (double) actual motor angle from QEI module, updated every timer interrupt
- PWMVal: (unsigned int) PWM value from 0 to PWMMAX
- PWMMAX: (macro) maximum value of PWM, say, 65535 for 16-bit PWM
- DIR: (macro) output latch corresponding to the pin used as direction bit for DC motor drive. The prototype board uses pin 25 (RB14)

- Pointer arithmetic is used for faster execution. You may argue that the improvement is insignificant for a small order controller in this example. If you are more comfortable with reference by value programming style, feel free to modify it.
- The control algorithm consists of two lines of code around the middle of timer routine. The code after that is just a conversion process from a single output variable u_0, to a pair of output signals PWMVal and DIR for the H-bridge driver. PWMVal is always positive integer, while the value of DIR changes with the sign of controller output. Suppose 16-bit PWM is used, for example, PWMMAX = 65535. So when u_0 = -73214.34, PWMValue = 65535 and DIR = 1. The motor turns in a negative sense; in this case, counter-clockwise.
- Be careful about indexing. Acess anything outside an array results in system crash.
- Pay attention to motor sense and quadrature signals from your encoder unit. The QEI module in PIC MCU can be setup to take care of different phasing of A,B signals. Improper configuration can result in positive feedback and closed-loop unstable.
- Make sure that your timer interrupt period matches the value used at discretization in step 6. Here we use 0.01 sec. Suppose you are using a low performance MCU, it could happen that the control algorithm could not be finished within 0.01 sec duration. Then you must use longer interrupt period. Redo step 6 with Ts set to your chosen value.
- It is a good idea to have an LED blinking to indicate the algorithm is running. Then you’ll know when problems occur like program crashes, too short interrupt period, etc. What I do is setting a global variable msc, and put these lines at the end of timer routine
msc++; // msc defined as global variable of type volatile int if (msc>100) { LED1 = !LED1; msc = 0; }

For 0.01 sec timer interrupt period, this causes LED1 to turn on and off each second.

- Study the mechanisms of the timer module in your MCU. Apart from the initialization, it may need some other attention to work properly. For the PIC24, the interrupt flag _T1IF must be cleared at the end of timer 1 routine so that next interrupt could occur.

That’s it. Compile and download the algorithm to MCU. With proper initialization, the DC motor should lock at an angle determined by command reference value when the feedback loop is closed.

Note : you may experience a “runaway” problem when the loop is closed the first time. That’s because the error between command reference and actual encoder count may be too large. Some initialization trick might be needed. Since in this experiment the DC motor is turning free with no reference point, I simply load the encoder counter with the command reference before closing the loop. You can do the other way around; i.e., read the encoder count and preset the command to that value.

Let’s have a close look at the target board I use for this experiment. As claimed earlier, no luxury hardware is required. Figure 8 shows snapshots of the top and buttom of this prototype. All components are soldered and wired manually. I use this board for other embedded experiments so there are more peripherals than needed. Actually, in this control experiment only the MCU and FTDI USB-to-serial IC are used.

To obtain motor angle feedback, quadrature encoder signals A, B are fed to QEI1 module of PIC24EP256MC202. Only MC (Motor Control) family of PIC24EP has QEI module. The GP (General Purpose) family does not. Read our PIC24EP QEI Basics, in addition to Microchip datasheet and family reference manual, if you want to learn how to setup and use the module.

The experimental results that follow are captured from this board via UART. Simple C functions are written to aid data communication in ASCII format. I simply use Hyperterminal program to save data in text files and copy/paste to Scilab.

The board costs about USD 20. Despite its messy appearance, this prototype serves our purpose well. I am working on PCB design for this board. That would take a while.

Figure 9 shows step response and controller output data from the experiment. The motor angle was commanded to move from 10 degrees to 100 degrees. Rough assessment from the figure gives rise time = 0.5 second, 40% overshoot, and setting time = 2 seconds. The lower plot shows the controller output is well below 16-bit PWM limit (65535).

Notice small undershoot at the beginning of motion. This is common for a non-minimum phase plant.

Before ending this article, we would like to address a design question. How to improve this step response? Let’s say we want faster motion. That means the controller (and the closed-loop system) must have higher bandwidth.

In classical control, the best solution is to redesign the controller. Lead/lag factors may be added to alter the frequency response of loop transfer function. For scheme, the controller needs to be re-synthesized by adjusting weighting functions for and .

Referring to our H-infinty Problem Setup article, to achieve higher closd-loop bandwidth, we need to move parameter in to higher frequency. Then has to be adjusted accordingly to allow higher bandwidth for . In other words, the cutoff frequency of must be higher than . Otherwise, a stabilizing controller may not be admissible.

We demonstrate by synthesizing another 2 controllers, call them C_HG1 and C_HG2 where C_HG1 has higher gain than our previous controller, and C_HG2 has the highest gain of all. For HG1 synthesis, the weighting functions are selected as

(6)

(7)

Running the synthesis yields the controller

(8)

The weighing functions for HG2 are

(9)

(10)

resulting in the controller

(11)

All we have to do next is replace the coefficients declared in the C source file with new set of values from controller (8), compile and download the hex file to MCU, then run the experiment and capture data. Repeat this procedure for controller (11).

A laborious and error-prone task is to extract the coefficients from a controller and put them in the C source file in correct order. So why don’t we let Scilab do it for us?

With a discrete-time controller Skdtf obtained from dcm_hinf_st.sce, extract coefficients from the numerator and denominator.

```
-->an=coeff(Skdtf.den);
-->bn=coeff(Skdtf.num);
```

With desired values contained in bn and an, we have to print them on Scilab console in correct order.

```
-->msprintf('an = {%f, %f, %f}',an(3),an(2),an(1))
ans =
an = {-2.692821, 2.425178, -0.732353}
-->msprintf('bn = {%f, %f, %f, %f}',bn(4),bn(3),bn(2),bn(1))
ans =
bn = {-422.248299, 1280.519620, -1290.339491, 432.068499}
```

Our only job left is to copy the value in bracket and paste to the C declaration

```
double coeffa[3] = {-2.692821, 2.425178, -0.732353};
double coeffb[4] = {-422.248307, 1280.519630, -1290.339489, 432.068496};
```

Step response and controller output comparisons from the 3 controllers are shown in Figure 10. Clearly, high gain controllers can improve the response time, and also reduce overshoot and undershoot. The controller output always needs to be monitored to make sure that no saturation occurs.

Step response is a basic property commonly used to evaluate tracking performance of a feedback system. Beyond that, one may want to experiment with more advanced responses such as ramp or parabola. It is not difficult to execute and capture any type of input response on our prototype board, provided that you can write C function for such input.

To illustrate, suppose we want to observe tracking performance of the 3 controllers to reference points generated by cubic polynomial function. A command generator is coded in C to produce array of points, which are inputted to the controller in place of the pcmda variable in the algorithm above. Both input and response data are sent back to PC via UART. Figure 11 shows resulting tracking response comparison. The dash-dot curve is the command generated by cubic polynomial function.

In this 2-part article, the whole process of H-infinity control from synthesis to implementation is demonstrated using a DC motor plant and a low-cost MCU prototype board. Controller synthesis, order reduction, and discretization are performed using Scilab. Then the discrete-time controller is rearranged to a DF II structure ready to be implemented as a C algorithm. Since most embedded products can be developed using C language, the algorithm can be easily ported to another MCU. The important requirement is it has to run inside a timer routine, with interrupt period equals the value used in discretization step.

Scilab files used in Part II

- dcm_hinf_st_hg1.sce script for S/T mixed sensitivity synthesis HG1 version
- dcm_hinf_st_hg2.sce script for S/T mixed sensitivity synthesis HG2 version

While I would not argue that the theory side of this so-called post-modern control could be mathematically involved, once one put some effort to grab the essence, the practical aspect is quite systematic and hence suitable for computer-aided analysis and design. In fact, the word “design” may be replaced by “synthesis” to indicate problem solvability with less human intervention. Of course, it is still mandatory for us to specify the requirements to a synthesis algorithm, typically in some form of weighting functions. The process afterwards is automatic. At the end we get our controller, or software complaint that none could meet our specs. Perhaps we have to relax this and that a little, and rerun the algorithm. The chance of success depends on plant complexity as well as how stringent the desired stability and performance criteria.

This article and the sequel address how to apply control paradigm to a simple setup, angle control of a DC motor. The complete process is demonstrated, from constructing weighting functions from stability and performance requirements, forming generalized plant, synthesizing a controller, performing order reduction, and implementing on an embedded target board. Being in the model-based control category, control requires a math model of the plant. So the first step is to acquire one using some type of system identification such as least-square method.

All the details involved were discussed in our previous articles. So this one is more like a how-to demonstration. For more insight or background, the following posts might be helpful

- synthesis with Scilab Part I : Problem Setup
- synthesis with Scilab Part II : Single-loop Tracking
- synthesis with Scilab Part III : Cascade Design
- Offline Least-Square System Identification

The DC motor is the same setup used in the DC Motor Speed Control Part I: Open-Loop Command article, which is shown once again in Figure 1. The plant consists of those parts mounted on pine wood: DC motor, shaft encoder, and H-bridge driver. The controller is implemented on the target board on the right, using PIC24EP256MC202 from Microchip as MCU.

Step 1: estimate a plant model

A function to generate PRBS input is implemented on the target board using (a C code example is provided on page 12 our System Modeling 101 article, then the input and output data is collected as shown in Figure 2. Note that the motor angle fluctuates around its initial angle at 180 degree by the PRBS input with strength +/- 5000. With 16-bit PWM and one direction bit, this translates to about 7% duty cycle, just enough to overcome shaft friction. Too large PRBS input strength might cause nonlinear motion or saturation.

Passing this input and output data to LSID algorithm in Scilab yields the following discrete-time plant

```
-->Pzsys_id
Pzsys_id =
2 3 4
- 0.0000014 - 0.0000006z + 0.0000054z + 0.0000158z + 0.0000102z
-----------------------------------------------------------------
2 3 4 5
- 0.0415585 + 0.1184330z + 0.2822806z + 0.1315280z - 1.4906741z + z
```

A standard controller synthesis is performed in continuous-time domain. Using a customized function as described in our Convert Transfer Function from Discrete to Continuous Time article results in a continuous-time plant model

```
-->P_id
P_id =
2 3 4 5
- 0.0009402 - 0.0000079s + 2.941D-08s + 1.146D-10s + 2.901D-13s - 9.479D-17s
--------------------------------------------------------------------------
2 3 4 5
- 0.0002877 - 0.0367550s - 0.0048099s - 0.0000557s - 0.0000002s - 2.500D-10s
```

with zeros at

```
-->roots(P_id.num)
ans =
3436.017
- 240.67648 + 307.73657i
- 240.67648 - 307.73657i
200.
- 94.567972
```

and poles at

```
-->roots(P_id.den)
ans =
- 306.919 + 216.09937i
- 306.919 - 216.09937i
- 123.56295
- 8.4363
- 0.0078361
```

Note that the plant model is non-minimum phase, and has a low frequency pole in place of an integrator that it should have in theory. Its Bode plot is shown in Figure 3.

Note: All the process above is contained in Scilab script dcm_lsid.sce

Before we move on, it is worth mentioning that in a real application, one should verify that the model achieved is a good candidate of the physical plant. A way to do so is to feed another set of rich input data (different from the one used in the identification process) to both the plant and model, and then compare the output data to see how well they match. Since system identification is not the main focus of our discussion, I choose to omit this verification step to save space and time.

Step 2: create weighting functions and form generalized plant

Standard S/T mixed sensitivity scheme is used for this problem. For initial design, let us make a crude stability and performance requirements, say

- normalized tracking error less than 0.1 degree at low frequency region below 0.1 Hz
- closed-loop bandwidth about 10 Hz
- sufficient stability margin (no high peak on )

With some trial and error on weighting function adjustments, we are eventually satisfied with

(1)

(2)

as weighting functions on and , respectively. In Scilab, these can be created as follows.

```
s=poly(0,'s');
a=0.002;
m=2;
wb=5;
w1_n=((1/m)*s+wb);
w1_d=(s+wb*a);
w1 = syslin('c',w1_n,w1_d);
w2_d=2*(s/500 + 1);
w2_n=s/100+1;
w2 = syslin('c',w2_n,w2_d);
```

See H-infinity Synthesis with Scilab Part I : Problem Setup article on the roles of parameters a, m, and wb and how to select them.

To observe the criteria on and , The inverses of weighting functions (1) and (2) are plotted

```
gainplot([1/w1; 1/w2],0.001,100,['';'']);
```

as shown in Figure 4. You can verify that they are fair candidates to the requirements listed above.

I cannot emphasize enough that weight selection is the crucial step of H_\infty control design. It is the last moment that we could tweak things to our desire, before the machine takes control. Some control theory background is helpful. If you are already good at classical control design, that same knowledge is inherited in the synthesis process by means of weighting functions. Otherwise, you may want to have a peek on Module 2 and Module 3 of our Scilab control engineering basics study modules.

Now that we have the plant model and weighting functions, a generalized plant can be formed. From my experience, converting all transfer functions to state-space systems and performing matrix augmentation gives best numerical result.

```
[Ap,Bp,Cp,Dp]=abcd(Pqv_a); // plant model (called P_id previously)
[Aw1,Bw1,Cw1,Dw1]=abcd(w1); // weighting function for S
[Aw2,Bw2,Cw2,Dw2]=abcd(w2); // weighting function for T
Agp=[Ap zeros(size(Ap,1),size(Aw1,2)) zeros(size(Ap,1),size(Aw2,2));
-Bw1*Cp Aw1 zeros(size(Aw1,1),size(Aw2,2));
Bw2*Cp zeros(size(Aw2,1),size(Aw1,2)) Aw2];
Bgp = [zeros(size(Bp,1),size(Bw1,2)) Bp;
Bw1 zeros(size(Bw1,1),size(Bp,2));
zeros(size(Bw2,1),size(Bw1,2)) zeros(size(Bw2,1),size(Bp,2))];
Cgp = [-Dw1*Cp Cw1 zeros(size(Cw1,1),size(Cw2,2));
Dw2*Cp zeros(size(Cw2,1),size(Cw1,2)) Cw2;
-Cp zeros(size(Cp,1),size(Cw1,2)) zeros(size(Cp,1),size(Cw2,2))];
Dgp = [Dw1 0.001; 0 0;1 0]; // makes D12 full rank
Pgen=syslin('c',Agp,Bgp,Cgp,Dgp);
```

Step 3: run the synthesis

A couple of Scilab commands could be used for synthesis. hinf and h_inf, for example, are two separate commands with different syntax. To use the former,

```
[Ak,Bk,Ck,Dk]=hinf(Agp,Bgp,Cgp,Dgp,1,1,5);
```

where 1,1 are the number of control input and measured output. In this case our controller is SISO. The last argument is gamma variable you can adjust. Roughly speaking, smaller gamma means more stringent criteria. If you get a controller that yields unstable closed-loop, try relaxing gamma to larger value.

Note there is no iteration in hinf. The command just runs once and returns a controller. In contrast,

```
[Sk,ro]=h_inf(Pgen, [1,1],0, 200000, 100);
```

iterates to yield a sub-optimal controller. Type help h_inf to learn what the arguments are.

Regardless of which command used, at the end you get a controller in state-space format (packed or unpacked).

Step 4: reduce controller order

synthesis generally returns a controller with higher order than necessary. From this DC motor example we get a controller of order 7. In an embedded application the resources are limited. High order means computation complexity and hence longer sampling period. Worse, the coefficients of higher order controller are subjected to numerical error, especially when using fixed-point, low performance MCUs.

So, we want to trim our controller to a lowest order one possible, without losing its properties significantly. How can we achieve such task?

A common approach is known as balanced truncation. Google it to learn the details. Here I just provide a Scilab function kreduced.sci ready for use. Load this function to Scilab workspace

```
exec('kreduced.sci',-1);
```

Pack the controller if necessary

```
Sk=syslin('c',Ak,Bk,Ck,Dk);
```

and pass it to kreduced. You will see on the console some message such as below

```
-->Skr = kreduced(Sk);
Warning :
matrix is close to singular or badly scaled. rcond = 1.9157D-10
150422.14
7713.887
115.87572
0.0339754
0.0003412
0.0000053
2.023D-09
Enter number of states to keep:
```

The computed Hankel norms might be different depending on your controller. The idea is to find the first large gap in the values and truncate there. In the numbers above, we observe that between 115.87572 and 0.0339754, the controller states have noticable energy drop. This says the balanced states above 3 are quite insignificant. Hence we can reduce the controller order to 3. Type the number and press [Enter] to yield the reduced controller contained in Skr. To observe this controller in transfer function form, use Scilab ss2tf command

```
-->Skr_tf = ss2tf(Skr)
Skr_tf =
2 3
384.79556 + 46179.923s + 1146.8162s - 500s
------------------------------------------
2 3
4.9087826 + 461.63448s + 31.256351s + s
```

Before getting into implementation on Part II, we end this article with the last step that should never be omitted.

Step 5: checking closed-loop stability

An unstable feedback system is not only useless, it could also be hazardous to the environment and humans in the workplace. Despite what the theory says, synthesis algorithms may give a controller that is not stabilizing, mostly when the weighting functions are not chosen properly, or too stringent criteria requested. Or, a full-order controller from synthesis could be degraded by reduction scheme to the point that closed-loop stability is lost. In any case, we always have to check this property right after synthesis.

The easiest method for a continuous-time SISO feedback system is to check whether the closed-loop poles are Hurwitz; i.e., have negative real parts.

In Scilab, construct the , , and transfer functions from the plant model and the controller *after model reduction*

```
-->L=Skr_tf*Pqv_a;
-->S=1/(1+L);
-->T = 1-S;
```

and check the poles of (or )

```
-->roots(T.den)
ans =
- 306.74046 + 215.68462i
- 306.74046 - 215.68462i
- 124.3588
- 13.073325 + 5.6429085i
- 13.073325 - 5.6429085i
- 6.9173679 + 1.1537611i
- 6.9173679 - 1.1537611i
- 0.0083338
```

The closed-loop frequency responses can be observed

```
-->gainplot([S; T],0.001,100,['';'']);
```

as shown in Figure 5. Verify whether they meet our specs given earlier. If not, adjust the weighting functions and rerun the synthesis.

In Control of DC Motor Part II : Controller Implementation we will discuss conversion of the controller to discrete-time and implementation on the target board.

Scilab files used in this article

- dcm_lsid.sce least-square system identification of DC motor
- dcm_hinf_st.sce script for H-infinity S/T mixed sensitivity synthesis
- kreduced.sci controller model reduction script
- download all as dcm_hinf.zip

- To lessen the effect of measurement noise, derivative part is implemented as a filter with parameter
- Back calculation anti-windup scheme is implemented with tracking gain
- Setpoint weightings for proportional and derivative paths can be adjusted via and , respectively

A feedback diagram with this advanced PID controller is constructed using Xcos palettes as in Figure 1.

In equation form, this controller can be described as

(1)

with

(2)

where , , and , are reference command, plant output, controller output, and saturated controller output, respectively. As described in our Discrete-time PID Controller Implementation article, using backward difference relationship

(3)

Equation (1) can be converted to z-domain as

(4)

Rewrite (4) in terms of

(5)

To implement this PID scheme as a computer algorithm, we have to convert (5) to a difference equation. It is straightforward to show that (5) can be rewritten as

(6)

with coefficients

(7)

So the corresponding difference equation is

(8)

Equation (8) is ready for implementation on a target processor. Before that phase, we want to make sure that our equation and coefficients are without error. One easy way is to perform simulation on Xcos and compare the response to the original continuous-time PID controller. For this purpose, we construct a model advpid_imp.zcos as shown in Figure 2, consisting of 2 feedback loops. The upper loop is controlled by discrete-time PID in the form (6), and the lower loop contains the continuous-time PID. The simulation results from the two closed-loop systems are then compared to verify how well they match.

Note that the discrete-time PID in the upper loop is contained in a superblock. The internal details are shown in Figure 3, which corresponds to the discrete-time controller (6).

Also, at the output of discrete-time PID controller, a LPF transfer function is inserted to prevent an algebraic loop error, normally occurred with hybrid simulation. The LPF pole is chosen well above the closed-loop bandwidth so the filter does not have noticeable effect on the responses.

For easy editing, all the parameters in advpid_imp.zcos are initialized using a script file advpid_imp.sce. The plant is chosen as a third-order lag transfer function

(9)

which can be created in Scilab by

```
s = poly(0,'s');
P = syslin('c',1/(s+1)^3);
```

Then, controller parameters are assigned values. These can be chosen as you like since the purpose of this simulation is to compare the responses. Here the PID gains are obtained from Zieglers Nichols frequency domain tuning method, and others are assigned some practical values.

```
kp = 4.8; // PID gains
ki = 2.7;
kd = 2.1;
N = 10; // filter coefficient
kt = 1.2; // back calculation gain for anti-windup
wp = 0.7; // setpoint weight for proportional term
wd = 0.1; // setpoint weight for derivative term
Ts = 0.01; // sampling peroid
```

For sampling period Ts, the value should match the simulation sampling period in the clock.

The parameters left to be assigned are the limits in saturation block. Put in some reasonable values such that some saturation effect happens during transient, since we prefer response comparison with the back calculation term activated. Too small the limit range would cause severe performance degradation. By some trial and error, we are finally satisfied with these values for saturation block

```
ulim = 2000;
llim = -2000;
```

Finally, the coefficients in (7) need to be computed. We introduce additional variables x1 and x2 for terms that appear in several places.

```
x1 = (1+N*Ts);
x2 = (2+N*Ts);
a1 = x2/x1;
a2 = -1/x1;
b1 = kp;
b2 = -kp*x2/x1;
b3 = kp/x1;
c1 = ki*Ts;
c2 = -ki*Ts/x1;
c3 = kt*Ts;
c4 = -kt*Ts/x1;
d1 = kd*N/x1;
d2 = -2*kd*N/x1;
d3 = kd*N/x1;
```

After all parameters are assigned, interactively or by executing advpid_imp.sce, we proceed by clicking on the simulation start button. The simulation results in Figure 4 show that the plant and controller outputs from continuous and discrete PID are almost identical. This makes us confident that the discrete-time PID and its coefficients are derived correctly.

After the verification process by Xcos simulation, we are now ready to implement our advanced PID controller on a target processor, in this case a PIC24EP256MC202 by Microchip. The plant is a DC motor with H-bridge drive, as described in the DC Motor Control Part I article. Figure 5 shows the experimental setup used. The rightmost board is our controller prototype, where the PID algorithm will be downloaded and executed.

The source code is written entirely in C. The whole code, consisting of a couple of source files, is rather long and messy due to supporting functions such as UART communication, command handling, etc. Below we discuss only the parts related to our PID controller implementation.

The sampling period, controller parameters and resulting coefficients are defined as global variables, with some initial values assigned

```
// sampling period
double Ts = 0.01; // sampling time
// -- these parameters are user-adjustable
double Kp = 1272; // proportional gain
double Ki = 8777; // integral gain
double Kd = 46; // derivative gain
double Kt = 10; // tracking gain
double Wp = 0.5; // proportional weight
double Wd = 0; // derivative weight
int N = 20; // filter coefficient
// ----- coefficients of PID algorithm --------------
double a1, a2, b1, b2, b3, c1, c2, c3, c4;
double d1, d2, d3;
```

and also variables to keep previous values of controller inputs and outputs

```
double ep_2, ep_1, ep_0, e_1, e_0, eus_1, eus_0, ed_2, ed_1, ed_0 ;
double u_2, u_1, u_0, u_0n; // variables used in PID computation
```

Now, the coefficients have to be computed before the algorithm starts, and every time the user changes any parameter involved. So it is convenient to put the computation in a function

```
void PIDSetup(void) // PID coefficient setup
// -- this function must be invoked anytime
// -- any parameter involved is changed by user
{
double x1, x2;
_T1IE = 0; // disable timer 1
x1 = 1 + N*Ts;
x2 = 2 + N*Ts;
a1 = x2/x1;
a2 = -1/x1;
b1 = Kp;
b2 = -Kp*x2/x1;
b3 = Kp/x1;
c1 = Ki*Ts;
c2 = -Ki*Ts/x1;
c3 = Kt*Ts;
c4 = -Kt*Ts/x1;
d1 = Kd*N/x1;
d2 = -2*Kd*N/x1;
d3 = Kd*N/x1;
_T1IE = 1; // enable timer 1
_T1IF = 0; // reset timer 1 interrupt flag
}
```

As usual, the actual PID algorithm is placed in a timer interrupt, in this case timer 1.

```
void __attribute__((interrupt, auto_psv)) _T1Interrupt(void)
// Timer 1 interrupt every Ts second
{
// perform position read from QEI module of PIC24EP256MC202
QEIpVal.half[0] = POS1CNTL; // read lsw
QEIpVal.half[1] = POS1HLD; // read msw from hold register
dcmpos = QEIpVal.full*360/ENCPPMx4; // position in degree
if (SysFlag.PosLoop == CLOSED) // closed loop PID control
{
u_2 = u_1;
u_1 = u_0;
ep_2 = ep_1;
ep_1 = ep_0;
ep_0 = Wp*pcmd - dcmpos; // weighted proportional error
e_1 = e_0;
e_0 = pcmd - dcmpos; // true error
eus_1 = eus_0; // back calculation error
if (abs(u_0) <= PWMMAX) eus_0 = 0;
else if (u_0>PWMMAX) eus_0 = PWMMAX - u_0;
else eus_0 = -u_0 - PWMMAX;
ed_2 = ed_1;
ed_1 = ed_0;
ed_0 = Wd*pcmd - dcmpos;
u_0 = a1*u_1+a2*u_2+b1*ep_0+b2*ep_1+b3*ep_2+c1*e_0+c2*e_1+c3*eus_0+c4*eus_1+d1*ed_0+d2*ed_1+d3*ed_2;
if (u_0>=0) { // positive sense
if (u_0 < PWMMAX) PWMVal = (unsigned int)u_0; // limit to PWM range
else PWMVal = PWMMAX;
DIR = 0;
}
else { // negative sense
u_0n = -u_0;
if (u_0n < PWMMAX) PWMVal = (unsigned int)u_0n; // limit to PWM range
else PWMVal = PWMMAX;
DIR = 1;
}
OC1R = PWMVal;
} // if (SysFlag.PosLoop == CLOSED)
_T1IF = 0; // reset interrupt flag
}
```

Note that our H-brige driver is commanded by a pair of signals: PWM and DIR, as explained in the DC Motor Control Part I article. The motor turns clockwise and counter-clockwise when DIR = 0 and 1, respectively, and the motor speed corresponds to the duty cycle of PWM signal, dictated by PWMVal variable.

An initial set of PID gains is achieved by performing some naive automatic tuning based on the Ziegler-Nichols frequency domain method. The C code is not shown in this article, though it is quite similar to that given in our Digital PID Controllers document. The automatic tuning yields . Other parameters are set to . As shown in Figure 6, This set of PID gains with rather high Ki value results in oscillary response (dashed red). So we begin fine tuning by reducing to , and , resulting in the dotted blue, and black, respectively. The overshoot and oscillation lessen with decreasing .

Next, we try incresing from , to , and , resulting in the responses in Figure 7. Oscillation is reduced with increasing .

Our last experiment for this article is varying the proportional setpoint weight . The PID gains are fixed at , , . Figure 8 shows the responses with , and . Interestingly, the overshoot can be reduced by decreasing , though the rise time becomes slower.

We leave it to anyone who implements this advanced PID controller on his/her system to experiment with the anti-windup back calculation gain and derivative setpoint weighting . It is suggested in some literature that be set to 0 so that abrupt change in command would not affect the derivative action of the controller. We set , and for all the experimental responses in this article.

]]>

From our previous article DC Motor Speed Control Part I: Open-loop Command, an experimental setup for DC motor is constructed and tested. The result shows that the command and actual motor speed differ significantly. Not a surprise. An open-loop control scheme cannot regulate the speed because no feedback is applied. So in this article, a PID control is deployed.

Figure 1 shows a simplified block diagram of closed-loop DC motor speed control. This diagram is drawn as a big picture for general audience to understand. Speed feedback is actually read from the QEI module of PIC24EP, so the leftmost summing junction should be inside the microcontroller. Anyway, we draw it explicitly to show the feedback structure more clearly.

The PID algorithm is a discrete-time version with derivative term replaced by filter transfer function

(1)

How to convert (1) to computer software algorithm is explained in Discrete-time PID Controller Implementation. The algorithm is put into ISR of timer 1. Note that the controller coefficients are functions of PID gains , filter coefficient , and sampling time . So whenever any of these parameters changes, the controller coefficients need to be recomputed. This can be written as a function

```
// ------------ PID control function ----------------------
void vPIDSetup(void)
// -- this function must be invoked anytime
// -- any parameter involved is changed by user
{
_T1IE = 0; // disable timer 1
a0v = (1+Nv*Ts);
a1v = -(2 + Nv*Ts);
// a2 = 1;
b0v = Kpv*(1+Nv*Ts) + Kiv*Ts*(1+Nv*Ts) + Kdv*Nv;
b1v = -(Kpv*(2+Nv*Ts) + Kiv*Ts + 2*Kdv*Nv);
b2v = Kpv + Kdv*Nv;
ku1v = a1v/a0v;
ku2v = a2v/a0v;
ke0v = b0v/a0v;
ke1v = b1v/a0v;
ke2v = b2v/a0v;
_T1IE = 1; // enable timer 1
_T1IF = 0; // reset timer 1 interrupt flag
SysFlag.VPIDchg = FALSE; // reset PID change flag
}
```

Note that variable names are ended with v to indicate they are used in velocity feedback loop, to distinquish from position feedback in case it is also implemented in the same system. To guard against read-modify-write problem, timer 1 is disabled during coefficient update.

For the PID algorithm, what is different from the code in our previous article is the output part. Now we have to convert a single control variable to a pair of signals PWMVal and DIR (see Part I). The conversion is straightforward. When the output is positive, it is sent out directly (limited by maximum possible value) and DIR is set to 0. On the other hand, when the output is negative, it is converted to positive value and DIR is set to 1. The code below shows how to implement the PID algorithm for our DC motor setup

```
e2v=e1v; // update variables
e1v=e0v;
u2v=u1v;
u1v=u0v;
e0v = scmd - srpm; // compute new error
u0v = -ku1v*u1v - ku2v*u2v + ke0v*e0v + ke1v*e1v + ke2v*e2v;
if (u0v>=0) { // positive sense
if (u0v>PWMMAX) u0v = PWMMAX; // limit to PWM range
PWMVal = (unsigned int)u0v;
DIR = 0;
}
else { // negative sense
u0vn = -u0v;
if (u0vn>PWMMAX) u0vn = PWMMAX;
PWMVal = (unsigned int)u0vn;
DIR = 1;
}
OC1R = PWMVal;
```

A global variable is defined so that the feedback loop could be open or closed by user command. Figure 2 shows step response of closed-loop system when the motor is commanded to rotate at speed 10 RPM. The controller gains are set at . We see that at steady state, the motor speed is regulated around 10 RPM, the commanded value. The response is not very smooth due to the speed readout is rounded to integer value. Comparing this to the open-loop response in part I, we see clearly the benefit of feedback.

Another important property of feedback control is disturbance rejection performance. Disturbance could be the result of, say, load change. An easy way to demonstrate output disturbance for this small DC motor is by grabbing the motor shaft with your fingers ( for your own safety, never attempt this on a high power motor). With no feedback, you can momentarily slow down or stop the rotation. In contrast, a closed-loop speed control commands the motor to generate more torque to combat the speed loss. You could feel the torque increase on your fingers.

Figure 3 shows disturbance rejection property of the closed-loop motor speed control. The motor is commanded to rotate at 100 RPM. Then at t = 5 seconds, I apply more load by grabbing the motor shaft. The speed regulation is disturbed, but then compensated by the controller so that the response remains regulated at 100 RPM.

In this article, we discuss speed control of a simple brushed DC motor with incremental encoder feedback. A standard PID algorithm is implemented on PIC24EP256MC202, a 70 MIPS microcontroller from Microchip. We focus on constructing an experimental platform using a small DC motor commercially available and off-the-shelf H-bridge driver. This motion control setup will be used in our subsequent articles

One way to construct a cheap experimental setup for motion control is to find a used DC servomotor with encoder equipped. Quality control of second-hand stores is a major problem; not everyone would be lucky to find a set in good-condition. There is no simple way to verify that the encoder is working okay. Fortunately, brand new DC motor with encoder unit can now be bought online within budget. An example is DC motor product line from www.pololu.com . Here we select their 19:1 Metal Gearmotor 37Dx52L mm with 64 CPR Encoder product as shown in Figure 1. This motor is equipped with a bare incremental encoder at the bottom as in Figure 2. A 12V power source is needed to drive the motor, and 5V for the encoder circuit. For other specs, please refer to the data on Pololu website.

There are many H-bridge motor driver products from various vendors, including Pololu, that could be used in this project. It happens that I have an old board model HMD-10A from Anycontrol Ltd, Thailand handy. It works fine in the desired setup. What makes this driver different from other products is the inverse polarity of PWM signal; i.e., the motor speed is proportional to the duty cycle of logic 0. So to make the controller algorithm less confusing, I use 1/6 of 74LS04 hex inverter at the drive input to make motor speed proportional to duty cycle of logic 1. The inverter also functions as a 3.3 to 5 volt level shifter and an buffer to the PWM output.

Since our simple experimental setup does not need any load at the motor output, it is convenient to mount DC motor, H-bridge drive and interface protoboard on a piece of wood as shown in Figure 3.

The board on the right (not screwed on wood) is our PIC32EP256MC202 module. The prototype is soldered and wired as in Figure 4. Such board could be assembled easily with budget around USD 20. Actually, not all components shown in Figure 4 are needed for the experiment in this article.

To perform an initial test for our DC motor setup, we attempt to drive the motor in an open-loop fashion. This involves sending a pair of signals: PWM for motor speed, and DIR for direction (0 = clockwise and 1 = counter-clockwise). The PWM module is setup for 16-bit resolution as explained in “A Note on Output Compare (PWM) Module of PIC24E” article, with one exception. In that article the PWM module obtains the clock source from the peripheral clock, whose frequency is too high for our DC motor driver. This makes it difficult to command low-speed turning of the motor. So we instead assign a timer output as clock source to the PWM. For example, to use timer 2, setup OCTSEL bits in the OC1CON1 register as follows

```
OC1CON1bits.OCTSEL = 0b000; // select Timer 2 as clock
```

Then initialize timer 2 to the desired frequency. Say, to set timer 2 interrupt at 8.7 KHz rate,

```
void initT2() // Initialize Timer 2 to use as PWM clock
{
_T2IE = 0; // disable first
T2CON = 0x8010 ; // Internal Clock (FOSC/2)=70MHZ,
PR2 = 1000; //ts/(tcy*64) where tcy = 1/70M
TMR2 = 0x0000 ;
_T2IF = 0;
_T2IE = 1;
}
```

Do not forget that, even though timer 2 has no particular task in its service routine, we need to reset _T2IF so that the next timer interrupt could occur

```
void __attribute__((interrupt, auto_psv)) _T2Interrupt(void)
{
_T2IF = 0;
}
```

The open-loop test needs no feedback. Even though the actual RPM is read from the QEI module of PIC24EP, the data is used for plotting a comparison only.

So, to create a reasonable speed command, we need to somehow find the relationship between the PWM value and motor speed. From the DC motor datasheet, at 12 Vdc the maximum speed of geared output shaft is 500 RPM with no load. This should be achieved when PWM value is at 100%. The PWM resolution is 16-bit. So let us define the following parameters

```
#define PWMMAX 65535 // maximum PWM for 16-bit resolution
#define MAXRPM 500 // maximum motor speed is 500 RPM
```

Construct a speed command variable scmd, in the range of +/- 500 RPM. Whenever this variable is updated by the user in real-time, the DIR value is determined by the sense, and PWM value is computed accordingly and written to the OC1R register. The C-code for such task can be written as follows (PWMVal is defined globally and DIR is the latch register of digital output pin used to control direction)

```
if (scmd < 0) { // negative sense
PWMVal = (unsigned int)-scmd*PWMMAX/MAXRPM;
if (PWMVal>PWMMAX) PWMVal = PWMMAX;
DIR = 1; // CCW direction
}
else { // positive sense
PWMVal = (unsigned int)scmd*PWMMAX/MAXRPM;
if (PWMVal>PWMMAX) PWMVal = PWMMAX;
DIR = 0; // CW direction
} // if (scmd < 0)
OC1R = PWMVal;
```

Now, to obtain the actual motor speed, we simply read from the VEL1CNT register of QEI1 module, and convert to RPM

```
QEIvVal = VEL1CNT; // read velocity each vsamp
srpm = 60*QEIvVal/(Ts*4*ENCPPM); // actual motor speed in RPM
```

where Ts is the sampling period, and ENCPPR is defined as total number of encoder pulses per round, 1200 for this motor. Note that the QEI module is set up to detect each rising and falling edges of quadrature signals A and B. This increases the resolution 4 times.

Set the speed command variable scmd to 10 RPM. The plot in Figure 5 compares the actual motor speed to the command. Well, from this comparison, the open-loop speed control performs poorly. The speed error is more than 150 %. This might come from certain factors. Our linear estimation of the PWM versus motor speed relationship could be inaccurate, or data from the manufacture is obtained from a different setup. etc.

In any case, without feedback one could not hope for a reliable speed regulation. In the next article, we will exploit PID feedback to cope with this speed error and combat disturbances such as when load changes.

]]>

In our article Scilab Recipe 4: GUI Basics, we discuss the process of designing a simple GUI for our PID feedback simulation. The prime function of this PIDGUI window is to allow easy adjustment of the controller parameters and saturation limits, then update the corresponding variables in Scilab workspace. The actual simulation and response plots are performed in the Xcos model.

Figure 1 shows a more complicated Xcos PID diagram advpid.zcos used for simulation. The diagram has 2 feedback loops for comparison. The upper loop contains a simple PID block provided by Xcos. The controller used in the lower loop, on the other hand, is an advanced PID with the following additional features

- To lessen the effect of measurement noise, derivative part is implemented as a filter with parameter
- Back calculation anti-windup scheme is implemented with tracking gain
- Setpoint weightings for proportional and derivative paths can be adjusted via and , respectively

Together with the 3 original PID gains, there are a total of 7 adjustable parameters in the advanced PID controller. Moreover, it would be nice if we could easily adjust the limits of saturation blocks to see the integrator windup effects. It is therefore convenient to handle all these 9 parameters by a GUI window.

Our PIDGUI application is shown in Figure 2. This window is created by executing the script pidgui.sce. The design is discussed in Scilab Recipe 4: GUI Basics on our Scilab Ninja website.

```
-->exec('pidgui.sce',-1);
```

The user interface should appear, with some default values for the parameters and their ranges. You can edit the script file to put in your choice of values. Note also that running the script file closes all other graphic windows. If that is undesirable, remove the command xdel(winsid()); at the top.

A nicer way to override the parameter values is to create another setup file with your plant and parameters tailored to that plant. The script advpid.sce contains the following commands that create a plant and parameters tuned by ZNFD method.

```
s = poly(0,'s');
// plant used in Astrom and Hagglund book
P = syslin('c',1/(s+1)^3);
// pid gains for above plant
// tuned by ziegler nichols frequency response method
ku = 8; // ultimate gain
tu = 3.5; // ultimate period
kp = 0.6*ku;
ki = kp/(0.5*tu);
kd = 0.125*kp*tu;
N = 10;
//
kt = 1.2; // back calculation gain for anti-windup
```

Launch this script

```
-->exec('advpid.sce',-1);
```

and import parameter values into PIDGUI by clicking the [Restore] button.

Now launch the Xcos model advpid.zcos in Figure 1. Click Simulation-Start. We would see a response comparison between the upper loop (standard PID) and lower loop (advanced PID) . They do not differ much at this time because the saturation limits are large. Try reducing them to, say, a limit for 12-bit DAC. In the Saturation Limits section of PIDGUI, put -2047 and 2047 to the lower and upper edit boxes, respectively. Run the simulation again. You ‘d see the result shown in Figure 3.

Here we list some experiments that could be done more conveniently with PIDGUI.

Change the variable names in Scilab To Workspace blocks for each simulation run, and plot them after finished.

At first launch, the saturation limits are set to +/- 1e6, which are so high like there is virtually no limit. Reduce them to +/- 10000, 5000, 2000, 1500, while maintaining other parameters. Figure 4 shows that saturation effects in this case does not cause larger overshoot, but introduce lags in the response.

The plot shown in Figure 6 of our article PID Anti Windup Schemes can be observed by adjusting the back calculation gain , click [update] and run the simulation.

In a sense, setpoint weighting can be thought of as a modest feedforward control. In the Xcos diagram, the error input to the proportional section of the controller is

(1)

i.e., the command is weighted by , with range 0 – 1. Figure 5 shows the responses with 5 values of : 1, 0.75, 0.5, 0.25, 0. We see that decreasing could improve overshoot, though too small a value worsens the rise time.

Similar to the proportional case above, the error input entering the derivative part is given by

(2)

with ranges between 0 – 1 as well. Figure 6 shows the responses with = 1, 0.75, 0.5, 0.25, 0. From this result, the advantage of using derivative weight on this plant and PID gain setup is not obvious.

It is left as an exercise for the reader to simulate responses from combination of these parameters.

We call this controller a decorated PID because it contains a couple more features than standard PID controller. The anti-windup function makes this controller nonlinear. This fact makes feedback analysis more challenging, but at least simulation results could help with the learning process.

PIDGUI is a simple user interface designed to aid parameter adjustments before simulation, so the reader could concentrate more on feedback system behavior. All files are provided below for you to download and experiment by yourself. The Xcos file accepts any continuous-time plant transfer function P created by syslin command.

Scilab and Xcos files used in this article

- pidgui.sce the user interface in Figure 2
- advpid.sce setup file for Xcos model
- advpid.zcos Xcos simulation model in Figure 1

Download all as zip file: pidgui.zip

]]>

Integrator windup in PID feedback control is well-known to control engineers. It happens when the integral term is used with some nonlinear saturation in the loop, such as when a physical variable reach its limit. When that happens, the feedback loop breaks, causing error accumulation in the integrator. This results in worse step and disturbance responses than the case of pure linear system.

To combat with the detrimental windup effects, a commercial PID controller often has some additional function called anti-windup. There are variations of implementation out there. Here we only investigate two basic schemes: * conditional integration* and *back calculation*.

Before further discussion, we have to review the two commonly-used forms of PID. First,

(1)

is often referred to as the “textbook” form. The second one is the parallel form used in Scilab Xcos block

(2)

Some tuning methods such as Ziegler Nichols use the form (1). In this article we prefer (2). This variation should not pose a problem since the parameters between the two can be related by

(3)

The plant used is the same one from our Digital PID Controllers article.

(4)

There we already performed PID parameter tuning. Those parameters will be used in the simulation.

To see how integrator windup deteriorates the step response, we use an Xcos diagram mypid.zcos shown in Figure 1. An output limit block is connected to the PID controller output, which can be switched on/off by a selector. The saturation is set up to + 511, a limit for 10-bit DAC or PWM. Further discussion of this setup can be read from our Scilab Recipe 3: Xcos Block Seasoning.

Setup the PID block by parameters achieved from the ZNFD tuning in Digital PID Controllers, and converted by (3).

```
ku = 8; // ultimate gain
tu = 3.5; // oscillation period
kp = 0.6*ku; // PID controller gains
ki = kp/(0.5*tu);
kd = 0.125*kp*tu;
```

The command input is a switching step command that starts from 0 to 1000 units at t=0.1 sec and switch to – 1000 units at t=20 secs. Run the simulation twice, first time as pure linear system (switch connected to input 1) and second time with output limit in effect (switch connected to input 2). Save the data to different variable names and plot a comparison. The result is shown in Figure 2. Notice that the windup effect for this system does not cause higher overshoot, but degrades significantly the rise time of step response.

To implement anti-windup, we will have to add extra mechanisms into our PID controller. So instead of using the PID block from Scilab palette, we have to construct our own, which is not difficult to do. Furthermore, the modified structure from our article Discrete-time PID Controller Implementation is used; i.e., the derivative part is replaced by Low-pass Filter for better noise immunity.

The idea of conditional integration scheme is simple. Whenever there is some indication that saturation is causing error accumulation, the integrator in PID controller is turned off. There is however some variation in deciding when to switch off. The easiest one to implement is perhaps monitoring the controller output and comparing with the limits. Whenever saturation occurs, the integrator is turned off.

The Xcos model awupid_m1.zcos in Figure 3 is built to experiment with this conditional integration scheme, implemented in the lower feedback diagram.

The plant blocks now accept data from a linear system P in Scilab workspace. So you can easily use the model with some other plant by constructing a new P. At this point we still use (4)

```
-->s = poly(0,'s');
-->P = syslin('c',1/(s+1)^3);
```

Notice that the derivative is implemented as an LPF with additional parameter N. We choose this as some large value, say

```
-->N = 1000;
```

Other PID gains remain the same as above.

Note: For your convenience, all setups are contained in a script file awupid.sce. Run the script once before simulation.

```
-->exec('awupid.sce',-1);
```

To see a comparison, connect the selector switch in the upper diagram to input 2. That routes the PID controller output through the output DAC limit block. Then click run. Figure 4 is the resulting scope view of output responses. Tracking performance is significantly improved with the conditional integration scheme, in the sense of overshoot reduction. However, it could not improve the rise-time of step response that degrades from the pure linear system.

There is no parameter to adjust in the conditional integration scheme. The integrator is simply switched on or off depending on whether a signal saturates.

Back calculation anti-windup scheme can be viewed as supplying a supplementary feedback path around the integrator. This feedback becomes active and helps stabilize the integrator only when the main feedback loop is open due to saturation. To experiment with this scheme, use an Xcos model awupid_m2.zcos in Figure 5. The back calculation feedback is implemented in the lower diagram, with additional compensation gain kt.

Figure 6 shows five output responses with kt = 0.1, 0.5, 1, 2, 5. The overshoot lessens with increasing kt, though too large a value causes slower rise-time, and even induces undershoot for the case kt = 5.

Similar to the conditional integration scheme, rise-time cannot be improved to match the case of pure linear system.

In this article, we discuss two anti-windup schemes in a PID controller and simulate some closed-loop step responses with a chosen plant transfer function. The audience can download Xcos files and perform further investigation with other choices of parameters and plants. We would also be interested to learn how these schemes behave under more realistic situations, such as tracking performance with other command inputs, with added input/output disturbances, or operating under noisy condition, etc.

- K.J. Astrom and T.Hagglund. Advanced PID Control Instrumentation, Systems, and Automation Society, 2006.
- V. Toochinda. Digital PID Controllers, www.controlsystemslab.com, 2011.
- V. Toochinda. Module 4: PID Control, in Scilab Control Engineering Basics, Scilab.Ninja, 2014.

Scilab and Xcos files used

- mypid.zcos Xcos diagram in Figure 1
- awupid_m1.zcos Xcos diagram in Figure 3
- awupid_m2.zcos Xcos diagram in Figure 5
- awupid.sce Scilab setup script

Download all as zipped file: awupid.zip

]]>For the continuous-time PID, we start with the so-called parallel form

(1)

with its Laplace transform

(2)

For practical implementation, it is quite common to modify the derivative term to an LPF filter, to make it less noisy

(3)

A straightforward way to discretize this controller is to convert the integral and derivative terms to their discrete-time counterpart. There are commonly 3 variations to do so, by means of forward Euler, backward Euler, and trapezoidal methods. Below we give representation for each.

Given a sampling period , the integral term can be represented in discrete-from by

* Forward Euler: *

(4)

* Backward Euler: *

(5)

* Trapezoidal: *

(6)

Similarly, the derivative term in (3) can be discretized as

* Forward Euler: *

(7)

* Backward Euler: *

(8)

* Trapezoidal: *

(9)

Obviously for all the terms above, the sampling period affects the gains of integral and derivative terms.

As an example, suppose we use backward Euler methods for both the integral and derivative terms, the resulting discrete-time PID controller is represented by

(10)

We want to simulate how this controller performs compared to its continuous-time version. We will use the setup in Figure 10 from our Module 4: PID Control. The plant consists of a robot joint driven by DC motor and a LPF at its input. Add the second feedback loop with discrete-time PID controller as shown in the Xcos diagram in Figure 1, or download dpidsim.zcos. Note that the controller is assembled from Xcos gain blocks, and the DLR blocks from discrete-time systems palette.

First we initialize some variables used by the model. You can do so using a script file or type in Scilab console directly.

```
Kp = 200; // these are values obtained from ZNFD tuning
Ki = 400;
Kd = 42;
Ts = 0.01; // must match sampling period used in simulation
N = 20;
```

Then the second and third terms from the controller expression in (10) is inputted to the DLR blocks for integral and derivative as shown in Figure 2 and 3, respectively.

Running the simulation yields the result in Figure 4. Notice that the step responses from the continuous (green) and discrete (red) PID controllers do not match, though they bear quite similar behavior.

Of course, one would not expect a good match since the standard continuous-time PID block provided by Xcos does not have a filter in the derivative term. So to make a fair comparison, it takes some more work to construct a custom continuous-time PID controller according to (3) as shown in Figure 5.

Download dpidsim2.zcos and run the simulation. The responses in Figure 6 depicts almost exact match between continuous and discrete PID control.

In order to implement on a computer, a discrete-time controller in the Z-domain must be transformed to its difference equation from as explained in our previous article Digital PID Controllers. Only now the process becomes more involved because the sampling time is embedded in the gains. Also, our controller now has the filter in its derivative term. It is straightforward for the reader to verify that the discrete-time PID controller (10) can be manipulated into the form

(11)

where and are controller output and input, respectively, and the coefficients are described by

(12)

From (11), we rearrange

(13)

to yield the difference equation for controller output

(14)

To make sure that our coefficient computation is correct, we construct another Xcos model dpidsim3.zcosas shown in Figure 7, where the discrete-time PID is formed using (12). It is convenient to do the computation in a script file dpidsim3_setup.sce and execute it before simulation.

The response comparison in Figure 8 shows quite a good match between the algorithmic PID and its continuous-time representation, though the discrepancy is somewhat more noticeable than in Figure 6.

Equation (12) can be programmed easily using your choice of computer language. For example, in C it would appear as follow

```
// global variables
double e2, e1, e0, u2, u1, u0; // variables used in PID computation
double r; // command
double y; // plant output
// -- these parameters should be user-adjustable
double Kp = 10; // proportional gain
double Ki = 1; // integral gain
double Kd = 1; // derivative gain
N = 100; // filter coefficients
Ts = 0.01; // This must match actual sampling time PID
// controller is running
// -- the following coefficients must be recomputed if
// --- any parameter above is changed by user
a0 = (1+N*Ts);
a1 = -(2 + N*Ts);
a2 = 1;
b0 = Kp*(1+N*Ts) + Ki*Ts*(1+N*Ts) + Kd*N;
b1 = -(Kp*(2+N*Ts) + Ki*Ts + 2*Kd*N);
b2 = Kp + Kd*N;
ku1 = a1/a0; ku2 = a2/a0; ke0 = b0/a0; ke1 = b1/a0; ke2 = b2/a0;
void pid(void) // implement this as timer interrupt service routine
{
e2=e1; e1=e0; u2=u1; u1=u0; // update variables
y = read(); // read plant output
e0 = r – y; // compute new error
u0 = -ku1*u1 – ku2*u2 + ke0*e0 + ke1*e1 + ke2*e2; // eq (12)
if (u0 > UMAX) u0 = UMAX; // limit to DAC or PWM range
if (U0 < UMIN) u0 = UMIN;
write(u0); // sent to output
}
```

In this article we discuss a practical discrete-time PID implementation, where the PID parameters are also functions of sampling time. The derivative term is commonly changed to an LPF to make it less noisy. Three conversion methods from continuous-time to discrete-time are generally used: forward Euler, backward Euler, and trapezoidal. Here we experiment only the backward Euler method. The reader is encouraged to simulate other methods and see which one gives the best match to continuous-time PID control.

We show in the last part how to implement a discrete-time PID controller as an algorithm on a computer or embedded system. In fact, the process can be applied to any type of controller described as a discrete-time transfer function in Z domain.

- Show that a continuous-time PID controller in parallel form (2); i.e., without a filter in derivative term, can be discretized using trapezoidal (also called bilinear transform or Tustin) method as
(15)

Hint: substitute in (2) and collect terms.

- Implement the controller (15) using your choice of computer language.

- dpid.zip : all Scilab and Xcos files used in this article

March 2015 : Correct the controller coefficient derivation and add dpidsim3.zcos to verify the result by response comparivon. Thanks to *Manuel Brüderlin, Chair for Computational Analysis of Technical Systems (CATS) Aachen*, who noticed that was missing from the coefficient equations of (11).

*Scilab-related articles are moved to our companion site: Scilab Ninja*

Scilab is a numerical computation software developed by Scilab Enterprises and has become my favorite choice for control engineering. Recently, so many posts on this site depend on this open-source software in some way. I now decide to construct another website Scilab Ninja specially for Scilab-related materials, including theRTSX (Robotic Tools for Scilab and Xcos). Some hightlights are Scilab Control Engineering Basics study modules that I use as supplement in my control design courses, and the use of scripts and Xcos models in certain control applications such as hardware-in-the-loop model development.

]]>Our articles on PID control are scattered around related websites. Here we provide links to them, which will be updated periodically.

- Digital PID Controllers: original article since June 2011
- eMotor Phones Home: implement the PID algorithm from above article on an electronic motor simulation platform.
- Discrete-time PID Controller Implementation: derive relationships between continuous and discrete PID gains as functions of sampling time, a discussion that is not covered in the original paper.
- PID Anti-Windup Schemes: two schemes for integrator-windup compensation are discussed and simulated.
- A Decorated PID Controller: Simulation of a full-featured PID controller made convenient with PIDGUI, a simple Scilab graphical user interface .

*Scilab Control Engineering Basics*- Module 4: PID Control: Learn how to simulate a PID feedback system with Xcos, effect of saturation, and parameter tuning with Ziegler-Nichols frequency domain method.

For reader unfamiliar with Scilab and Xcos, the software we use in most articles, check out our Scilab Recipe Series.