FB pixel

Embedded C Programming with the PIC18F14K50 - 20. Automated Temperature Control System with PID Control. Part 3 - Practical WorkNew

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

Assembled plant
Figure 9 - Assembled Plant

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

First order system response delay unit step input
Figure 10 - Response of the first order system with delay to the unit step input

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

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.

Plant static characteristic
Figure 11 - Static Characteristic of the Plant

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

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

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

Plant response 100 percent step input
Figure 12 - Response of the plant to the 100% step input

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

Plant zoomed response 100 percent step input
Figure 13 - Zoomed response of the plant to the 100% step input

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.

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

Table of Equations
Table 4 - Calculation of the PID Controller Parameters

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.

Table of Equations
Table 5 - PID Controller 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.

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.

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.

System response minimal integral error
System response minimal integral error
Figure 14 - System response with the minimal integral error

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.

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.

System response 20 percent overshoot
System response 20 percent overshoot
Figure 15 - System response with the 20% overshoot

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.

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.

System response optimal parameters
System response optimal parameters
Figure 16 - System response with the optimal parameters

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.

Make Bread with our CircuitBread Toaster!

Get the latest tools and tutorials, fresh from the toaster.

What are you looking for?