Optimal PI control parameters for accurate underfloor heating temperature control

In low energy buildings, the effect of internal and solar gains on heat balance of rooms is large. As a result, the heating systems, designed assuming steady-state conditions with no heat gains, are over-dimensioned for most of the heating period. This poses a challenge for room-based control systems, especially for thermostatic valves, but also for PI controllers. Using over-dimensioned room units might result in room temperature fluctuations. For finding solutions to this problem by using simulations, correct modelling of the control system together with the room is crucial. Therefore, the aim of this research was to determine the challenges that occur while matching measured and simulated temperature profiles and test the effect of PI control parameters on the calibrated model control accuracy. The experiments were carried out for the underfloor heating system of a test building. The building was simulated in IDA-ICE software and calibration via minimising root mean square error of energy consumption in GenOpt was carried out. The PI parameters were fit by optimisation with objective to simulate the measured temperatures accurately. The effect of the optimized PI parameters was determined by comparison to IDA-ICE default parameters and parameters from Cohen-Coon method.


Introduction
It is common practise to over-dimension the heating systems in buildings for the sake of safety margins. Moreover, the current load calculation standard does not take into account the internal gains making the estimated power demand even higher [1]. Most of the year, these oversized heating systems result in reduction of efficiency [2] and one of its causes is control inaccuracy due to large unnecessary control region while using typical underfloor heating controllers such as thermostats and PI controllers. Although, recently literature discusses self-learning, adaptive and predictive controllers, still many of the "smart" controllers make use of the classical controllers to reach and maintain the calculated setpoint they suggest [3].
Tuning the simple controllers has been a relevant topic for decades [3] and current self-learning PI controllers make use of it by running several tests and calculations to auto-tune the parameters online. However, in building energy simulations, varying the parameters is not common as ideal systems are used or parameters are set to default values and the effect is often under-estimated [4]. Finding better parameters than default ones, requires expert knowledge in control theory and calculations for identifying suitable parameters. Therefore, we aim to optimize the parameters via optimization software GenOpt [5] to enable tuning the parameters for simulations without any prior knowledge of the field. Optimising the parameters in simulation does not enable to omit the experts setting the parameters in the actual buildings if the model of the building is faulty. In this paper, we first calibrate a small test house model based on measurements. That being a complex, often trialand-error problem, we introduce optimization algorithms and parametric runs here as well. In this research, we used the nearly zero energy test house located at Tallinn University of Technology campus as object for measurements and simulations. It is a 110-m 2 building with five test rooms, all equipped with underfloor heating; the design heating power at -22℃ for each room is shown in Fig. 1 and the floor construction is shown in Fig. 2. The house in equipped with a heat recovery ventilation with the design air exchange of 1.8 h -1 . Both walls and roof include several partial elements with test constructions that possible introduce additional leakage. Maivel et al. have described the building in detail in [6]. We collected the important characteristics into Table 1.

Fig. 2. Floor Construction layers with wet installation
underfloor heating location specified in red and rebar in black thick line

Experimental setup
Measurements were carried out to test two different kinds of control systems. Self-learning PI control with unknown parameters and simple thermostats with on-off control:  Test 1: PI control (Feb 23 -Mar 11)  Test 2: On-off control (Mar 11-21) During the experiments, the doors in between the rooms and the windows were closed; the venetian blinds were drawn for the conference room (R9) windows except the northern window but not in other rooms. Ventilation was set at level 2 that should correspond to the design air exchange; however, the airflow rates were not measured. Internal gains were generated with dummies in R5, R6, and R9. Their input power was measured. Underfloor heating system was used during the tests and radiators were turned off. All the circuits were maximally open due to the initial experiment characteristics. For more detailed information, see [7].
Local weather station data are available from site. Solar data measurements were taken from a high, unshaded roof 500 meters from the test house. Diffuse and direct radiation on horizontal surface was measured. Python library pysolar was used to translate the direct radiation from horizontal to normal [8].

Simulation setup
The model of the test building was created in IDA-ICE 4.8 [9] and the calibration process was carried out as described in section 2.4. Based on the calibrated model and the initial model, new models were built with internal gains of the Estonian Standard [10] and climate as in Estonian Test Reference Year (TRY) [11]. On these models, PI parameters were optimized and tested. This process is described in detail in section 2.5.

Calibration process
For evaluation, we chose a period of three consecutive full days (midnight to midnight) of data when the solar data from local measurements existed and showed a similar pattern over the two test periods. For on-off control, it is March 12-14 and for PI control February 27 to March 1.
For the initial model, set temperature during on-off tests was kept constant for the test rooms whereas for the PI tests the night periods from 9 p.m. to 3 a.m. were with slightly lower setpoints than daytime periods as the selflearning PI carried out some tests for parameter improvement. During the night periods, test rooms' setpoints were set to the measured temperatures and for other rooms the temperatures were optimised. For corridor (R2), and rooms R7 and R8, the measured temperatures were used as setpoints at all times.
In our simulation, the internal gains from dummies were forced to the measured values. However, in the simulation the internal gains were delayed by a first order transfer function to model the imperfect mixing of the air in the room. In rooms R2 and technical room (R3) the generated internal gains from computers and other equipment were not known but calibrated with optimisation methods.
The supply temperature of underfloor heating was set to the measured profile. Infiltration was measured but due to several known leakages that were taped during the blower door test, the actual infiltration should be calibrated.
The procedure for the calibration of the model was comprised of the following steps: 1. define and simulate model with all known design parameters, 2. identify possible causes of errors in temperature profiles and energy consumption, remodel with best estimations, 3. vary set temperatures to fit temperature profiles, 4. optimise varying envelope and ventilation parameters to fit total energy consumption and power profile.
The first two steps contain classical methods of reading the building plans and project and building the house in IDA-ICE.
Step 2 involves initial simulation tests and studying the results for reasons of mismatch. A parametric study is carried out in step 3, testing a range of temperatures automatically. For step 4, we applied GenOpt®'s hybrid method. Although, comparison of several different methods is suggested in [12], for the initial tests for the parameter analysis, we consider the accuracy of one method sufficient.
The energy consumption was optimised for the chosen three-day periods. The varied parameters are the infiltration rate, ventilation air exchange rate, thermal bridges and extra infiltration added to R9, because of some holes in the wall.

Control analysis methods
Estonian TRY and the Estonian standard profiles for internal gains were used for a full heating period of simulation to analyse the control parameters. Set temperature was 21 ℃ in all rooms; controls in all rooms had the same parameters at the first stage. The ventilation air exchange rate was 0.42 l/s/m 2 everywhere and the minimum exhaust air temperature was set to 5 ℃. The external blinds were omitted.
The heating period was estimated using other parameters from the improved base model of the PI case. The beginning of the next period was counted when the heating had been on/off for more than a week consecutively. The test simulations showed that the heating period is from October 1 to April 30.
To analyse the effect of a calibrated and noncalibrated models on control parameters, several envelope characteristics were deliberately falsified. Two wrong models were generated by exchanging parameter values in the calibrated case: 1. design case: design (not optimised) infiltration, air exchange rate, and thermal bridges 2. over-dimensioned case: over-dimensioned heating power (to test parameters optimised for significantly lower power) For each wrong as well as the calibrated model, control parameters were optimised to achieve the best fit of room air temperatures. The parameters obtained from the design case were tested on the calibrated and overdimensioned models to achieve the effect of calibration. The objective function is the sum of the root mean square errors (RMSEs) of the room air temperatures.
PI parameters were also calculated for all models. For rooms R5 and R9, step tests were carried out at outdoor the temperature of -22 ℃ (design temperature) and no solar radiation. No internal gains were applied. From the step response, parameters were calculated using Cohen-Coon method [13]. The gain (K) and time constant (Ti) were calculated as follows: where, gp is the process gain (percent change in the temperature divided by percent change in the control parameter), τ is the time constant (time from the end of dead time to the time when temperature reaches 63% of its total change), and td is the dead time (time from changing the control parameter value to the time that temperature rise tangential meets the initial temperature level). All of these variables are determined as explained in [13]. Steps turning heating signal from 0 to 1 as well as 0.5 to 1 were carried out to achieve higher quality average values because of the flaw of the Cohen-Coon method to result in large changes in the calculated parameters when small changes in the input data are made. The calculated parameters for the design case were tested on the over-dimensioned case to obtain the maximum error except for using default parameters. The control parameters and performance indicators obtained for the calibrated model were compared with those from design calculations. The aim was to determine if the optimization would give similar or better results. This would enable use of this method for parameter obtaining without step tests.
Comparing energy consumption is not reasonable for the case of all setpoints set to 21 ℃ as the time of temperature below the setpoint depends on control. For a realistic case where the lower limit of the temperature is kept, setpoint correction is performed. The correction allows the temperature to stay below the setpoint 3 % of the time. Therefore, if the temperature value at 3 % time of cumulative temperature graph is θ ℃ below 21 ℃ the corrected setpoint is 21+θ ℃.

Initial model
First, all the parameters defined in design were used as inputs. Some parameters were not known, such as the internal gains in rooms R2 and R3 as well as time constants for internal gains and on-off controller. The quality of input parameters for the simulation is shown in Table 2. The measured and simulated data did not correspond from the beginning. Initially the total energy consumption for the three test days was 29 % higher than measured in the PI case and 23 % higher for the on-off case.

Improved model
Observing the forced room temperatures in R7, R8 and R2, we can see that none of the rooms keeps it well. One example is shown in Fig. 3 from R7. To ensure the rooms R7, R8 and R2 as realistic boundary conditions for the test rooms, the temperature profiles of these rooms have to fit. For that, an initial search for better control parameters was carried out. For PI case, a couple of parameters close to the default ones were tested to find better-performing set. The RMSE in R7 decreased while increasing gain and decreasing integration time. It is logical as a fast-reacting PI is needed. Gain of 1 and minimum integration time of 10 seconds were chosen to keep the simulation times in reasonable limits. Until further notice, these parameters were used for all rooms. In the on-off case, thermostat dead band was varied and logically the lower it is the more exact the profile. A minimum value of 0.05 ℃ was used. From the results, we can see that not all the rooms achieved their set temperatures at all times (see 26.02 in Fig. 3). We assumed that the design values for underfloor heating reflected the demand and the losses downwards were not included. To take that into account the resistances of the floor construction below and above the underfloor heating pipes were calculated and correspondingly the power was raised by 5 %. It improved the situation only partly but from design, we have no reason to use other values. Different boundary conditions could evoke this effect, for example different temperatures for some rooms, but we will see this later.
The temperatures in R2 are on the contrary too high at night times both in PI and on-off cases (see Fig. 4). This probably occurred due to the massiveness of the floorheating while in the experiment there were computers in that room generating free heat. We fitted the available constant power by a parametric run with objective function being RMSE of the temperature profile. Parametric run resulted in 120 W constant heat gain for PI case and 250 W for on-off case. As this is mostly constant in practice, a power value was chosen that had the minimum sum of the errors at these two cases, which was at 120 W.  Initially, we set the temperatures for rooms R1, R3-R6 and R9 to 21 ℃ but most probably, the technical room (R3) and bathroom were having higher temperatures, although not measured. Moreover, for the PI cases, we can see from mass flow data that the controller was running different tests and learning during the nights. Therefore, in the test rooms we force the set temperature to follow the measured temperature from 9 p.m. to 3 a.m. and daily keep 21 ℃. As the actual circuit in R3 and R4 is the same, the control output from R3 is directed also to R4. The same applies to R1 and R2, so the control output from R2 is assigned to R1. For room R3, set temperatures were 23 ℃ for the day and 19 ℃ for the night periods to simulate the higher temperatures and the learning PI. After implementing all these improvements, the total energy use difference still existed but it had decreased to 19 % for PI and 13 % for on-off case.

Optimal model
The test rooms' set points had been a little off from the aimed 21 ℃, so the setpoints of rooms R5, R6, R9 were optimised by minimising the RMSE of the temperature profiles. For PI case, the nightly setpoints were kept as measured and the optimal daily setpoints were 21.4 ℃, 21.7 ℃ and 21.1 ℃, respectively. For on-off, these temperatures were constant all the time with values 21.3 ℃, 21.4 ℃, and 21.1 ℃. The different solar peaks of the periods could cause the gap between the setpoints of room R6 that is a small room with large south and west windows. The parameters to vary for fitting the total energy consumption to the measured one were:  setpoints for room R3,  internal gains in R3,  ventilation air exchange rate,  infiltration,  thermal bridges,  leakage in R9 wall (there are some holes in the wall).
For fitting the power profile, the objective function was the RMSE of the energy meter values (the cumulative power profile). To ensure that the internal gains in R3 were not too high, the RMSE of R3 temperature profile with respect to the set temperatures was added to the objective. The values were in similar order of magnitude but mostly, energy dominated in the objective function value.
Optimising power profile in this way did not work for PI, as the resulting energy consumption was different from the measured case (see Fig. 5). The simulated energy demand would have been 11 % lower from the measured one (almost as different as not optimizing at all). The power profile itself did not fit well either as there were still many unknowns such as bypass and pump control. In the on-off case, the total energy consumption difference was however only 0.2 %. Nevertheless, here the RMSE of the power profile was 2146 W. We can deduce that with the available knowledge of the building, it is extremely complicated to fit the power profile in an exact manner. Therefore, further we will focus on the total energy consumption and temperature profiles.
Next, the sum of absolute difference between the measured and simulated total energy consumption and the RMSE of R3 temperature profile was used as the objective function. The difference from total energy consumption for PI optimal case is 0.1 %, for on-off 1.6 %. These results were obviously better than for the power optimisation case.
As most of the varied parameters were characteristic to the building, a mutual optimum was looked for. For that, parametric run around the optimal values was carried out and the case with the lowest sum of the two objective function values was chosen where both errors separately were below 2. The optimal temperatures in R3 were 26 ℃ in on-off case and 25 ℃ for PI. In the mutual optimum PI case, the optimal temperatures changed to 23 ℃ at night and 25 ℃ in the daytime.  The mutual optimum parameter values are shown in Table 3. The resulting temperature profiles (from R9 as an example) and energy meter data are shown in Fig. 6 and Fig. 7. For mutual optimum, the difference from total energy consumption for PI case increased to 0.2 % and for on-off case, it decreased to 1.3 %. As the differences were very small, these could be counted as almost equal. The reason for increase could be the fact that using only one optimization method, for simulation based optimization the optimum was not correctly found at all times. The GenOpt manual [12] suggests to use several different methods and to improve the optimums with parametric runs around them. This will remain as further work in following papers.

On-off control
For all the on-off cases, optimal dead band value was the possible minimal (we used 0.05 ℃) and the same applied for the time constant (minimal was 0 s). As the parameters were the same, only the error value and energy outputs could be compared after setpoint correction. The corrected setpoints varied from room to room, ranging from 21.00 ℃ to 21.04 ℃ for all except R1 where corrected setpoints were 21.05 ℃ to 21.07 ℃. After setpoint correction, sum of all RMSE values of temperature profiles for design case was 3.34 ℃, Table 4 shows. For the calibrated version, 3.23 ℃ and for overdimensioned case, 3.27 ℃ RMSE values apply. Total energy consumptions for underfloor heating were 48.0, 53.2, and 53.4 kWh/m 2 /y for design, calibration and overdimensioning, respectively. The difference of energy usage was therefore 10 % for the design case and 0.4 % for the over-dimensioned case.

PI control
The PI model in IDA-ICE enables to change in addition to gain and integration time also the tracking time and the time constant. The time constant characterizes the control difference filtering. We set it to 0 seconds representing no filtering. Tracking time (Tt) represents the integration time of the values over the maximum or below the minimum PI output limits. We set it to be equal to the integration time as is usual in classical PI and as the parameter' calculation algorithm assumes as well. Tested on calibrated model with otherwise default parameters, setting tracking time to default 30 seconds versus to integration time (300 seconds in the default case) changed the RMSE from 3.237 K to 3.263 K. The energy consumption per heating period increased by 8 kWh (from 5272 to 5280 kWh). As we expected these differences to be marginal, we stayed with the integration time although the default gave a bit better results for RMSE. On the two remaining parameters, the gain and integration time, we ran the optimisation by minimising the sum of temperature RMSEs via GenOpt, the results are shown in Table 5 together with the default parameters as well as parameters calculated from Cohen-Coon method. It can be seen that optimal and calculated gain (K) and time constant (Ti) values are much higher than the default ones. The default values are aimed for any kind of heat emitter and mostly for fast-reacting ones. The higher gain characterises the need for more aggressive control for the underfloor heating compared to other heaters. In longer integration time, the inertia of the underfloor heating shows. The reason for the difference between calculated and optimal cases is not obvious. However, we can assume that the step test at low constant ambient temperature and dynamic operation could result in fundamental control differences. Table 6 shows all the simulated results for the heating period for three model types and parameter obtaining methods. In addition, the three crossed cases for effect determination are shown in last rows. For each, both initial and after setpoint correction results are shown. Only for the last three rows, it does not apply, as here only the RMSE is crucial. The corrected setpoints were similar to on-off cases ranging from 21.00 ℃ to 21.07 ℃. After the correction for over-dimensioned case with PI default parameters, the 3 % limit was exactly 21 ℃ for all rooms.
The PI parameter effect on energy demand was marginal. The energy consumption results were very similar both for initial and corrected cases comparing between different parameters on same model. The difference stayed below 0.5 % from the optimal value.
The RMSE values in both Table 4 and Table 6 are higher than 3 ℃ as it is the sum of RMSEs of nine rooms. Therefore, a RMSE of 3.6 ℃ would mean on average 0.9 ℃ RMSE per room. The values were that high due to excessive solar gains in autumn and spring season as there was no cooling nor external shading applied. The temperature profile across the heating period for R9 is depicted in Fig. 8.  Nevertheless, a tendency in RMSE can be highlighted. The RMSE for default parameters was 3-5 % higher than for optimal parameters, and 1-3 % higher for the calculated parameter cases. If we compare the RMSE values from 21 ℃ also for the corrected case (red values in tables), these percentages were 9-13 % and -3 % to 1 % due to different temperature levels. The calibration, which effect is illustrated by the cases where design case optimal parameters were tested on calibrated and overdimensioned models, showed a logical increase of the RMSEs. However, the values were extremely low, 0.3 %. This could be because only a small number of variables had been changed in the falsified models as well as high RMSE from sunny days hiding the results from other periods.
The realistic maximum error that could be made in design stage is described by the last case where parameters calculated for the design case were used for the over-dimensioned case. This showed the RMSE increase of 2.2 %. As the energy consumption was similar for both, the optimisation method could be used instead of the calculation method. All the described RMSE results are depicted in Fig. 9. However, to compare the performance on cloudy cold days, the analysis was done for the calibrated case in January. Corrected set temperatures were used. The RMSE error values were here around 1 ℃ compared to more than 3 ℃ for the total heating period. The temperature profile differences in January can be seen in Fig. 10. For the shown week in R9, the RMSE values for the default and of-off cases were 0.2 ℃ and 0.19 ℃, whereas the PI control with calculated and optimal parameters showed RMSE values of 0.14 ℃ and 0.13 ℃.

Comparison
The difference in total energy consumption for on-off and PI cases was small. However, the temperatures fluctuated significantly less for calculated and optimized PI parameter control than for default PI parameter or on-off control. The RMSE of on-off cases was 6-9 % higher than of optimal PI parameter cases over the heating period. This is illustrated also by the temperature profiles in January in Fig. 10. One of the reasons for small difference in energy consumption could be the very small value of the on-off dead band. This presents an idealistic on-off controller. It ensured that the temperature fluctuations were very small and it could result in lower energy consumption than realistically possible.

Conclusions
Obtaining the PI control parameters through optimisation proves successful as the RMSE results show a small improvement even compared to the parameters calculated from the Cohen-Coon method. Testing the different control parameters has shown that the effect of the parameter changes on the temperature profile is up to 5 %. Optimal and calculated PI parameters had significantly higher gain and longer integration time compared to default ones, and resulted in improved control accuracy characterised by lower RMSE from the set temperature. Despite of differences in control accuracy, in energy use the on-off control with the optimal dead band and time constant is comparable to the PI with default parameters for the observed underfloor heating case. This research also shows that the calibration of a model of the building together with its heating system is possible but challenging. The energy consumption and temperature profiles result in a good fit while a significant difference remains in the power profile.