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

Published

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

## System Assembly Recommendations

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.

## Plant Parameters Definition

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 Kp is the plant gain, Tp 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 Tp, 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 Kp 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.

## Measurement of The Plant Parameters

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).

### Table 2 - Static Characteristics of the Plant

 Y,% T_cur, C 0 26 20 34 40 42 60 50 80 57 100 64

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 Kp 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 Kp 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).

### Table 3 - Response of the Plant to the 100% Step Input

 t, s T_cur t, s T_cur t, s T_cur t, s T_cur t, s T_cur 0 26.3 210 47.9 420 57.3 630 61.5 840 63.3 10 27.3 220 48.6 430 57.6 640 61.6 850 63.4 20 28.9 230 49.1 440 57.9 650 61.7 860 63.5 30 30.5 240 49.8 450 58.1 660 61.8 870 63.5 40 31.9 250 50.3 460 58.5 670 61.9 880 63.5 50 33.2 260 50.9 470 58.6 680 62 890 63.6 60 34.5 270 51.5 480 58.9 690 62.1 900 63.6 70 35.6 280 52 490 59.1 700 62.2 910 63.6 80 36.8 290 52.5 500 59.3 710 62.3 920 63.7 90 37.8 300 52.9 510 59.5 720 62.4 930 63.7 100 38.9 310 53.4 520 59.8 730 62.5 940 63.8 110 39.8 320 53.8 530 60 740 62.6 950 63.8 120 40.8 330 54.3 540 60.1 750 62.7 960 63.8 130 41.7 340 54.6 550 60.3 760 62.8 970 63.9 140 42.6 350 55 560 60.5 770 62.8 980 63.9 150 43.5 360 55.4 570 60.6 780 62.9 990 64 160 44.3 370 55.7 580 60.8 790 63 1000 64 170 45 380 56.1 590 61 800 63 1010 64 180 45.8 390 56.4 600 61.1 810 63.1 1020 64 190 46.5 400 56.7 610 61.2 820 63.2 200 47.3 410 57 620 61.3 830 63.3

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 Ti to it. Also, if the time is less that 𝜏p the output should stay Ti. 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).

## Calculation of the PID Controller Parameters

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).

## Experiments with the PID Controller

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 td = 80s (green cell in Table 6).
2. Rise time tr = 180s (cyan cell in Table 6).
3. Peak time tp = 220s (yellow cell in Table 6).
4. Peak overshoot
5. Setting time ts = 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 Kp.

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 td = 90s (green cell in Table 7).
2. Rise time tr = 230s (cyan cell in Table 7).
3. Peak time ts = 430s (yellow cell in Table 7).
4. Peak overshoot
5. Setting time ts = 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 td = 120s (green cell in Table 8).
2. Rise time tr = 260s (cyan cell in Table 8).
3. Peak time ts = 330s (yellow cell in Table 8).
4. Peak overshoot
5. Setting time ts = 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.

### Table 6 - System response with the minimal integral error

 t, s T_c Y t, s T_c Y t, s T_c Y 0 35 20 130 45.8 100 260 49.9 83.3 10 35.9 100 140 46.4 100 270 50.1 16.2 20 37 100 150 47 100 280 49.8 100 30 38.1 100 160 47.5 100 290 50.1 38.2 40 39.1 100 170 48 100 300 49.8 91.2 50 40 100 180 48.5 100 310 49.8 83.8 60 40.9 100 190 49 100 320 50.1 16.2 70 41.7 100 200 49.5 100 330 49.7 100 80 42.5 100 210 49.9 83 340 50 67.7 90 43.1 100 220 50.1 25 350 49.9 57.4 100 43.9 100 230 49.8 100 360 49.8 100 110 44.5 100 240 50.1 38.2 370 50.1 29.4 120 45.3 100 250 49.9 91.2 380 49.8 100

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

### Table 7 - System response with the 20% overshoot

 t, s T_c Y t, s T_c Y t, s T_c Y 0 35 20 180 47 85.5 360 49.8 52.2 10 35.5 100 190 47.5 72.6 370 49.7 82 20 36.1 100 200 47.8 67.7 380 49.8 66 30 37.1 100 210 48 90 390 49.9 72 40 38.1 100 220 48.3 64.3 400 49.8 76.2 50 39.1 100 230 48.5 66.5 410 49.9 59 60 40 100 240 48.6 89.7 420 49.9 62 70 40.7 100 250 48.8 57.9 430 50 43 80 41.5 100 260 49 63.5 440 49.9 90 90 42.3 100 270 49 93.9 450 50 48.8 100 42.9 100 280 49.3 52.7 460 49.9 72.2 110 43.6 100 290 49.3 73.3 470 50 52.3 120 44.3 93.9 300 49.3 94.7 480 49.9 75.5 130 44.8 100 310 49.5 62.7 490 50 54.7 140 45.2 82.1 320 49.5 79.2 500 49.9 78.5 150 46 70.8 330 49.6 69.2 510 50 57.1 160 46.3 87.8 340 49.7 57.4 520 50 57.1 170 46.7 83 350 49.6 89.9 530 50 59.5

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

### Table 8 - System response with the optimal parameters

 t, s T_c Y t, s T_c Y t, s T_c Y 0 27 0 130 39.5 100 260 47.8 100 10 27.6 100 140 40.3 100 270 48.3 97 20 28.9 100 150 41 100 280 48.7 100 30 30.1 100 160 41.8 100 290 49.1 99 40 31.3 100 170 42.5 100 300 49.5 75 50 32.1 100 180 43.2 100 310 49.7 81 60 33 100 190 43.8 100 320 49.9 83 70 34 100 200 44.5 100 330 50 55 80 35 100 210 45.1 100 340 49.9 80 90 36 100 220 45.6 100 350 50 66 100 36.9 100 230 46.2 100 360 50 66 110 37.8 100 240 46.7 100 370 50 66 120 38.6 100 250 47.3 100 380 50 66

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