Real-life Implementation of Temperature Control with PID | Embedded C Programming - Part 22

Hello again! This is the last part of the PID control system tutorial. The first two parts can be found here and here.

Finally, the time has come to do some practical work. So please assemble the system according to figure 4 of part 1. The difficulty may happen with attaching the DS18B20 sensor to the 5W resistor R5. I used a heat shrink tube in which I made a small hole to insert the sensor. After heating, I got the following construction (figure 9).

You can also add some thermal paste between the sensor and the resistor to decrease the thermal resistance.

When everything is assembled, you need to apply external power PWR to the board and optionally the VCC voltage. After that you can power the device on, and see the LCD view like in figure 8 of part 2.

Now our goal is to calculate the parameters of the regulator Kp, Ki, and Kd which depend on the plant parameters. We will consider the plant as a first order system with a delay.

The plant gain G(s) (Figure 1) is:

Where K_p is the plant gain, T_p is the plant time constant, 𝜏_p is the plant delay time. If the first two parameters may be familiar to you from this tutorial, then the last one is something new. Let’s consider the response of the first order system with delay to the unit step input (figure 10).

As you see, it’s a regular response to the step input of the first order system but it’s shifted along the time axis by the delay time 𝜏_p.

So the delay time is the time during which the system doesn’t react to the input because of its inertia. Temperature is a parameter with a significant amount of inertia, especially if the system is big enough that the delay time can reach several tens of seconds.

In the theoretical system presented in figure 10, 𝜏_p = 4s.

The time constant T_p, as explained here, is “the time required for the output to attain 63.2% of its steady state value” but I’ll show you another way to find its value.

So you need to draw the tangent to the graph curve at the point where the speed of the increasing of the output value is maximum (the derivative of the output is at its local maximum). In Figure 10 it’s point A. Then we prolong the tangent to the horizontal line to which the output is aimed. At the intersection of these two lines, there is point B. The time difference between points B and A is the time constant.

So in our case

And finally the plant gain K_p is the ratio of the difference between the end and initial values of the system output to the difference between the end and initial values of the system input. As the input is the unit step, its value is 1 - 0 = 1, and the output (according to figure 10) is 0.5 - 0 = 0.5.

So

Well, that’s all good, but these are theoretical values, and now let’s calculate these three parameters for our specific plant.

We will make two experiments and draw two graphs that will describe the plant parameters: static characteristics, and the response of the plant to the step input.

The static characteristic is the dependence between the output and input values of the system in a steady state. To draw it, we need to turn on our system and leave it in manual mode (the letter “M” in the right bottom corner of the LCD should be present). Then you should create a table somewhere with two columns: Y, % (input of the plant), and T_cur, degrees (output of the plant). Then write down the first row, where Y = 0% and T_cur is the environment temperature.

Then using the encoder, set Y = 20% and wait for several minutes before the system goes to the new steady state, and the temperature stops changing. Then write the second row.

Do the same things with the Y = 40%, 60%, 80%, and 100%. This may take some time but what can you do?

Finally, you should have the table like this (Table 2).

As you see, this characteristic is almost linear, which means that in all temperature ranges the response of the system to the input control impact Y will be the same. This is good for us because possible non-linearity of the plant may cause certain problems with the PID controller.

This characteristic is important to us because we can see from it, which temperature corresponds to which control impact Y. Also we can calculate the K_p value for any part of the characteristic. For instance, in the beginning, it’s:

And in the end it’s:

The difference is (0.4 - 0.35) / 0.4 x 100 % = 12.5%, so we can consider this nonlinearity as insignificant.

Now let’s do the second part of the experiment and measure the response of the plant to the step input. The first thing to do is to set the Y to 0% and wait till the resistor cools down to the initial temperature.

Now, when it’s done, we need to set the Y value as 100% as fast as possible to emulate the step input. We will not set it as 1% because the response will be too small. And as the static characteristic is almost linear we can set any step value, and just divide by it to calculate the K_p value.

After that, we need to write down the T_cur value every 10 seconds until it becomes stable. Actually, the sign of ending the process is the equality of the last four values. Finally, you should get a table like this (Table 3).

Your T_cur values will be most likely different from what I have, so refer to Table 2 just as an example, not as the values to use.

Based on Table 2, let’s draw the static characteristic of the plant (figure 11). I used Microsoft Excel for this but you can use any software you’re comfortable with.

The blue line is the graph made based on the values from Table 3, and the orange graph is the approximated response based on the calculated plant parameters, which we will find very soon.

Let’s start with the gain. It’s better to calculate it based on Table 3 first and last values.

The time constant can be easily found in figure 12:

I wrote the "≈" sign here because it’s hard to measure it using this graph with high precision.

As for the delay time, it’s almost impossible to find it using Figure 12. So I zoomed in on the first part of the graph to find it (figure 13).

Here you can see clearer that the tangent is touching the green line at the moment of the highest temperature raising speed. Also even here you can see that the approximation orange curve is very close to the practical one.

From figure 3 we can find that 𝜏_p ≈ 4s.

Now, we have all the parameters of the plant. Let’s just out of curiosity make the equation of the approximated orange curve. Let’s convert equation (12) from the Laplacian form into the normal one.

The step input Y(s) in Laplace form is the following:

Then the plant output is

With the method of partial fractions (the same as in this tutorial) we can rewrite it:

Now let’s apply the inverse Laplace conversion to the equation (15):

This is the exact equation that approximates our curve but we need to add the initial value of the temperature T_i to it. Also, if the time is less that 𝜏_p the output should stay T_i. So the final equation is:

If we put the parameters of our plant into it, we can get the orange curve, in my case it’s:

Well, I hope your curiosity is now satisfied, and we can move on and find the PID controller parameters.

As you can see, for me the process took the 1020s, or 17 min!

Now, let’s draw this characteristic and determine the parameters of the plant (figure 12).

There are a lot of methods to calculate the parameters of the PID controller knowing the plant parameters, for instance, these ones: https://en.wikipedia.org/wiki/PID_controller#Overview_of_tuning_methods.

Some of them require quite complex calculations and are more precise, the others are approximate but quite simple. We will use the second approach because it’s enough calculations for today: I don’t want to take away Kushal Gowda’s bread and butter :)

So there is an approximate method that allows us to find the parameters of the PID controller based on the plant parameters and the type of the system response (Table 4).

1. Let’s start with the last column of Table 5 and set the following PID controller parameters (lines 23-25 of the “main.c” file):

const int32_t Kp = 223; //Proportional factor
const int32_t Ki = 39; //Integration factor
const int32_t Kd_inv = 186; //Inverted derivative factor (1 / Kd)

Then we recompile the firmware and download it into the MCU. After starting the system we press the encoder button to switch to the automated mode. We should see that the letter in the bottom right corner of the LCD has changed to “A”. As the initial value of the setpoint is the minimal possible temperature, the PID controller output Y will remain 0.

Now rotate the encoder handle fast to set T_s to about 15-25 degrees higher than it is now. For example, if my initial temperature is 35 degrees, I set the new value as 50. Then, the same as with the previous experiment we need to write down the values of the current temperature T_c, and additionally the values of the Y. We do this until the system doesn’t switch to the new stable state, or until we see that this is impossible. The results of this experiment are listed in Table 6.

As you can see, these parameters don't work very well because the system becomes unstable. Well, the situation is not that bad, the system doesn’t go to pieces, but it keeps oscillating forever which is unacceptable for a real control system. Let’s calculate the system response specifications described in this tutorial.

1. Delay time t_d = 80s (green cell in Table 6).
2. Rise time t_r = 180s (cyan cell in Table 6).
3. Peak time t_p = 220s (yellow cell in Table 6).
4. Peak overshoot
5. Setting time t_s = 200s (red cell in Table 6), if we consider the allowed tolerance

2. Now let’s switch to the third column of Table 5 and change the PID controller parameters (lines 23-25 of the “main.c” file):

const int32_t Kp = 207; //Proportional factor
const int32_t Ki = 3; //Integration factor
const int32_t Kd_inv = 172; //Inverted derivative factor (1 / Kd)

In this set of parameters, the Ki value differs the most from the previous one. Let’s see how this will affect the system response.

The method of the experiment is the same as the previous one. Its results are listed in Table 7.

So here are the formulas of the main PID controller parameters: proportional Kpr, integral Ki, and derivative Kd factors. Here I changed the designation of the proportional factor to Kpr not to mix up it with the plant’s gain K_p.

There are several types of system responses that mainly differ with the overshoot. You can read more about the system response characteristics here.

The aperiodic process doesn’t have any overshoot at all but it has the longest settling time. So the system is very stable but quite slack. The process with the minimal integral error has a 40% overshoot but it has the shortest settling time, and as follows from its name, the minimal integral error. The other two process types are transitional between the first two.

Let’s now calculate the parameters of the PID regulator for our specific case. Please note that your values will be different because you will have other plant parameters.

As the method is quite approximate, the experimental tuning is a normal part of it. So we first set some parameters, see how the controller works, then change them and see how that affects the system response.

Now, the system response finishes without the oscillation despite the chaotic control impact Y. Also there is no overshoot. Moreover, you may notice that the response time becomes much longer. Let’s calculate the system response specifications for this case.

1. Delay time t_d = 90s (green cell in Table 7).
2. Rise time t_r = 230s (cyan cell in Table 7).
3. Peak time t_s = 430s (yellow cell in Table 7).
4. Peak overshoot

5. Setting time t_s = 430s (yellow cell in Table 7, merged with the peak time).

So, as you can see, the system performance got worse by almost double but the oscillations have gone.

3. As I said initially, we need to adjust the parameters in an empirical way to get the best performance. Let’s analyze the two experiments that we’ve made. The Kpr and Kd values changed insignificantly, but the Ki reduced more than 10 times: from 39 to 3. As follows from the experiments’ results, it affects the system response crucially. So to get the best system performance without the oscillation we need to set some interim value of Ki. I made several experiments and got the best results with the parameters that are presented in the initial program code:

const int32_t Kp = 200; //Proportional factor
const int32_t Ki = 10; //Integration factor
const int32_t Kd_inv = 200; //Inverted derivative factor (1 / Kd)

Let’s see the results of this experiment.

As you see, now the system has even better performance than in the first experiment. Even though the time seems to be longer, this time the initial temperature value is 27 degrees instead of 35, so some time is required to reach this value. Also, you may notice that the Y value became less random and got to the required value very quickly, so no overshoot happened as well.

Let’s calculate the system response specifications for this case.

1. Delay time t_d = 120s (green cell in Table 8).
2. Rise time t_r = 260s (cyan cell in Table 8).
3. Peak time t_s = 330s (yellow cell in Table 8).
4. Peak overshoot

5. Setting time t_s = 330s (yellow cell in Table 8, merged with the peak time).

Well, that’s almost all I wanted to tell you in this tutorial. What I want to mention is that to increase the accuracy of the PID controller, it’s better to reduce step h. In our case, we can’t do that significantly, as the DS18B20 temperature measurement time is 750 ms. You can decrease the sensor resolution (please read in its datasheet how to do this) and thus decrease the measurement time to 87ms but lower sensor resolution leads to more chaotic changes of the Y value. So for best performance, it’s better to use a faster sensor, even one with an analog output.

And now that’s really all. This tutorial seems really gigantic in comparison to the previous ones (but not to the next ones).

As homework, I suggest you first repeat all the steps that I made but with your own parts. And second, make more experiments trying to change the Kpr, Kd, and h values and see how they affect the system response.

Let’s now draw the graphs that correspond to Table 6.

Let’s now draw the graphs that correspond to Table 7.

Let’s now draw the graphs that correspond to Table 8.