PI Parameter Inﬂuence on Underﬂoor Heating Energy Consumption and Setpoint Tracking in nZEBs

: In rooms with underﬂoor heating (UFH), local on–o ﬀ controllers most often regulate the air temperature with poor accuracy and energy penalties. It is known that proportional–integral (PI) controllers can regulate most processes more precisely. However, hydronic UFH systems have long time constants, especially in low-energy buildings, and PI parameters are not easy to set manually. In this work, several potential PI parameter estimation methods were applied, including optimizing the parameters in GenOpt, calculating the parameters based on simpliﬁed models, and tuning the parameters automatically in Matlab. For all found parameter combinations, the energy consumption and control precision were evaluated. Simpler methods were compared to the optimal solutions to ﬁnd similar parameters. Compared with an on–o ﬀ controller with a 0.5 K dead-band, the best PI parameter combination found was with a proportional gain of 18 and an integration time of 2300 s, which could decrease the energy consumption for heating by 9% and by 5% compared with default PI parameters. Moreover, while GenOpt was the best method to ﬁnd the optimal parameters, it was also possible with a simple automatic test and calculation within a weekend. The test can be, for example, 6-h setbacks applied during the nights or weekend-long pseudo-random changes in the setpoint signal. The parameters can be calculated based on the simpliﬁed model from these tests using any well-known simple method. Results revealed that the UFH PI controller with the correct parameters started to work in a predictive fashion and the resulting room temperature curves were practically ideal.


Introduction
The change towards nearly zero-energy buildings (nZEBs) and renewable energy sources influences the technologies used for heating and its control [1,2]. The intermittent production of renewable electricity calls for flexibility in all consumers, including buildings [3]. Space heating is responsible for up to 70% of the final energy demand in residential buildings [4]. Therefore, it has a high potential for flexibility. In modern buildings, the use of heat pumps has intensified [5]. Only electricity-based heating is relevant to the power grid, therefore, heat pumps are a clear target.
To be exploited when the grid needs it, heat pumps should use an electricity price or other signal for optimizing their performance. Some of the heat pumps already optimize their behaviour according to the price. As one solution to improve the flexibility, model predictive control (MPC) can be used [6,7]. It enables the use of historic and forecasted data to predict the most optimal course of action. At the occurrence of renewed data, the optimization can be corrected. For an MPC for a single-family house, parameters for UFH in nZEB and analyzes their effect on the energy performance and indoor air temperature of the building. PI performance is compared with a traditional thermostat's performance in the same situation. Both an accurate temperature tracking performance and a considerable energy saving compared with conventional control are expected. The results may be utilized in the design of UFH systems with accurate temperature control and energy savings compared with conventional UFH systems.

The Building
The work is based on a test building at TalTech University campus, which is described in detail in several previous publications [25][26][27]. Two almost identical rooms with a floor area of 10.4 m 2 were analyzed, except that one of them (Room 6 or R6) has two 4 m 2 windows facing south and west, while the windows of the other (Room 5 or R5) face north and west. The floor plan of the building is shown in Figure 1 with the two test rooms highlighted with red rectangles. Previously, the test house model in IDA ICE 4.8 software [28] was calibrated against measured air temperatures in the test room R5 during temperature setback cycles with varying durations [27]. As a result, the heat losses and thermal mass of the room structures are adequately defined in the model. This model was used for the simulations in the current work. In the simulations, all of the other rooms were heated constantly with ideal heaters to the setpoint of 21 • C.
Energies 2020, 13, x FOR PEER REVIEW  3 of 21 thermostat's performance in the same situation. Both an accurate temperature tracking performance and a considerable energy saving compared with conventional control are expected. The results may be utilized in the design of UFH systems with accurate temperature control and energy savings compared with conventional UFH systems.

The Building
The work is based on a test building at TalTech University campus, which is described in detail in several previous publications [25][26][27]. Two almost identical rooms with a floor area of 10.4 m 2 were analyzed, except that one of them (Room 6 or R6) has two 4 m 2 windows facing south and west, while the windows of the other (Room 5 or R5) face north and west. The floor plan of the building is shown in Figure 1 with the two test rooms highlighted with red rectangles. Previously, the test house model in IDA ICE 4.8 software [28] was calibrated against measured air temperatures in the test room R5 during temperature setback cycles with varying durations [27]. As a result, the heat losses and thermal mass of the room structures are adequately defined in the model. This model was used for the simulations in the current work. In the simulations, all of the other rooms were heated constantly with ideal heaters to the setpoint of 21 °C. The building has wooden-frame walls, a wooden-frame roof, and concrete floors with a crawl space below. The total heat-up time constant for the rooms is around 85 hours and the effective time constant for temporary setbacks is around 12 h [27]. The absolute cool-down time constant of one test room is around 24 h when the other rooms are heated constantly. The time constant for the whole building cool-down is ca. 100 h. The time constants are long mainly due to the concrete floor and highly insulated building envelope. The values were confirmed by the experimental data presented in [27].

Outline of the Work
The PI parameters were estimated for the two test rooms in several different ways. Firstly, they were optimized in GenOpt with the aim of minimal setpoint tracking errors both for the constant and variable setpoints (Section 2.5). Secondly, they were calculated and estimated using simplified models. The data used for the model fitting are described in Section 2.3 and the model fitting process is described in Section 2.4. The models were used to either autotune the parameters in Matlab or to calculate the parameters using well-known methods such as AMIGO, SIMC, and Cohen-Coon. Both of these approaches are also clarified in Section 2.5. The performance of all the gained parameters was cross-checked in both rooms over the whole heating period. The analysis is described in detail in Section 2.6. The building has wooden-frame walls, a wooden-frame roof, and concrete floors with a crawl space below. The total heat-up time constant for the rooms is around 85 h and the effective time constant for temporary setbacks is around 12 h [27]. The absolute cool-down time constant of one test room is around 24 h when the other rooms are heated constantly. The time constant for the whole building cool-down is ca. 100 h. The time constants are long mainly due to the concrete floor and highly insulated building envelope. The values were confirmed by the experimental data presented in [27].

Outline of the Work
The PI parameters were estimated for the two test rooms in several different ways. Firstly, they were optimized in GenOpt with the aim of minimal setpoint tracking errors both for the constant and variable setpoints (Section 2.5). Secondly, they were calculated and estimated using simplified models. The data used for the model fitting are described in Section 2.3 and the model fitting process is described in Section 2.4. The models were used to either autotune the parameters in Matlab or to calculate the parameters using well-known methods such as AMIGO, SIMC, and Cohen-Coon. Both of these approaches are also clarified in Section 2.5. The performance of all the gained parameters was cross-checked in both rooms over the whole heating period. The analysis is described in detail in Section 2.6.

Input Data
All the data used for the PI parameter calculations are summarized in Table 1. In this section, only the grey area is described, the rest is tackled in the following sections. Here, the data from [27] were used, where the authors performed temperature setbacks with different lengths in the test building. The air temperature during setbacks with durations of 2 days and 3 days was measured in room 5, where the temperature setpoint was normally kept at 21 • C and during the setbacks was lowered to 18 • C. In the calibrated IDA ICE model, shorter setbacks of 1, 3, 6, 12, and 24 h were simulated using a constant outdoor temperature of 0 • C, with no solar and internal gains. Between the setbacks, the initial temperature of 21 • C was stabilized. Without solar gains, the two test rooms are equivalent and therefore, the PI parameters estimation is based on only one of them. In addition, an ideal-like step test was simulated with the same constant outdoor conditions. A step from no heating to full power heating was performed. The simulation period was prolonged for so long that the stability of the indoor air temperature was achieved both before and after the step. This meant two months in simulation to stabilize at the balance temperature, and one month after the step for reaching a steady state.
Additionally, simulations with Estonian test reference year (TRY) [29] and pseudo-random binary signal (PRBS) as setpoints were used. For the PRBS temperature setpoint, the zero level was set on 18 • C and the maximum level on 24 • C. The simulations were done for two separate weeks, one in March and one in February:
The model fitting was done both on the entire weeks and only on the weekends of these weeks (12 p.m. Friday to 12 p.m. Sunday).
For the optimization (the last two rows in Table 1), the same two weeks of Estonian TRY were used as well as the whole heating period from 1 October to 30 April. The setpoints for the optimization cases are the same as used for the evaluation and are described in Section 2.6.

Model Fitting
A simplified process model of the system is needed to use most of the PI parameter calculation methods. Based on the generated input data, a first order process model with a time delay was fitted. Therefore, the temperature response of an input step change is where θ(t) is room air temperature in • C at time t seconds after the step, θ(0) is the initial temperature before the step, K p is the process gain (unitless), T is the time constant, and L is the time delay, both in seconds. The model fitting was performed in Matlab using System Identification Toolbox [30].

Estimating PI Parameters
The PI parameters K and T i were estimated, where K is the proportional factor and T i is the integration time of the integral part of the PI in its ideal form: where u is the control signal (unitless) and E is the difference between the setpoint and measured air temperature in • C that is feedback to the control. For all the cases in Table 1, the PI parameters were estimated by one or more of the following methods: 1. Optimized using GenOpt; 2.
Calculated from an applicable simple method.
In the optimization method, the PI parameters were optimized in GenOpt using a hybrid GPS algorithm [31]. The optimization was carried out for the three different periods described previously and two different setpoint profiles, which are also used for the evaluation and are described below in Section 2.6. The objective of the optimization was to minimize the average absolute difference between the setpoint temperature and the simulated temperature.
In the second method, the PI parameters were auto-tuned in Matlab ® /Simulink for the previously fitted simplified models (described in Section 2.4). The tuning was performed aiming for a short rise time (speed) and overshoot of no more than 5% of the desired temperature increase.
In the third method, all the models that had been fitted based on the different input data were used to calculate the PI parameters. Three widely known methods-Cohen-Coon, Skogestad IMC (SIMC), and AMIGO-were used for that. The PI parameters K and Ti are calculated according to these methods as follows [32]: Cohen-Coon (CC): Skogestad IMC (SIMC): AMIGO: (8) where Kp, L and T are the parameters from the fitted models with the general representation in Equation (1). The parameters a and τ are unitless parameters: Energies 2020, 13, 2068 6 of 20

The Evaluation Tests
All the estimated PI parameter combinations were tested in simulations in both test rooms. The accuracy of the setpoint tracking was assessed on both the constant and variable setpoints. The constant setpoint was chosen to be 21 • C and the variable setpoint was calculated from price data 2017-2018 [33], based on the simple algorithm given in [34] that does not perform the best for their purpose of load shifting but gives us an hourly changing setpoint profile. In the price-based control, the air temperature setpoint is changed hourly between 20, 21, and 24°C. The lower two levels are meant for comfort and have to be met at all times, the highest level is implemented for load shifting and does not need to be tracked. All evaluations were done for the whole heating period (01 October-30 April). All combinations of PI parameters, both rooms, and both setpoint profiles were evaluated based on: • The average absolute error (AAE) of the air temperature from the setpoint; • The heating energy consumption per square meter of the floor area.
For the energy consumption comparison, it is important that no parameter combinations would result in temperatures lower than the given comfort setpoints. In most cases, this was not achieved and, therefore, the setpoints had to be shifted. The goal was to achieve temperatures equal or above the setpoint for at least 97% of the time, as suggested in the thermal comfort standard EN 16798-2 [35]. Based on the initial simulations, cumulative temperature graphs were generated. In the constant setpoint case, the setpoint was shifted exactly as much as the cumulative graph was, below the setpoint at 3% of the time. For the variable setpoints, shifts for both the two 20 • C and 21 • C setpoints were calculated. The 3% of the 20 • C was at 1.3% of the total time and for the 21 • C setpoint at 45.2% of the total heating period length. The maximum of the shifts calculated for these two points was applied to the whole profile.

Benchmarks
The simulation software IDA ICE's default PI parameter values K = 0.3 and T i = 300 s were used for the benchmark simulations. Furthermore, on-off controls with four different dead-band widths were evaluated for the comparison. A modern one with a dead-band of 0.5 • C was used, but also close to ideal versions, with dead-bands of 0.16 • C and 0.05 • C and a conservative one with a 1 • C dead-band, were used as well.

PI Implementation in IDA ICE and PI Mechanics
As the PI controller can be implemented in various formats, the implementation in IDA ICE is shown in Figure 2. The example code in Figure 2A is modified for the case where error filtering is turned off, the mode is heating, and the conversion unit equals 1. The parameter tt, the tracking time, is set to 30 s.
The hilimit and lolimit are the limits for the PI output signal. In this work, the PI output signal is the fraction of the nominal mass flow to the UFH and is, therefore, limited from 0 to 1. In Figure 2B, an increase in the sample air temperature over the setpoint, i.e., due to solar gains, can be observed. In Figure 2C,D, the calculation of the script can be followed. The lines are colored according to the variable text colors in the script.
In Figure 2C, it can be observed that, even though the temperature is over the setpoint between 3 a.m. and 5 a.m. (Figure 2B), the PI signal is not zero. It only gets to zero when the integral part also decreases so much so that the sum of the integral and error parts is less or equal to zero. Although the OutSignal is limited, the negative values of OutSignalTemp are still used for the calculation. This enables the effect, which looks like prediction in some cases. This effect is further discussed in Section 3.3.

PI Implementation in IDA ICE and PI Mechanics
As the PI controller can be implemented in various formats, the implementation in IDA ICE is shown in Figure 2. The example code in Figure 2A is modified for the case where error filtering is turned off, the mode is heating, and the conversion unit equals 1. The parameter tt, the tracking time, is set to 30 s.
D. The hilimit and lolimit are the limits for the PI output signal. In this work, the PI output signal is the fraction of the nominal mass flow to the UFH and is, therefore, limited from 0 to 1. In Figure 2B, an increase in the sample air temperature over the setpoint, i.e., due to solar gains, can be observed. In Figure 2C, and 2D, the calculation of the script can be followed. The lines are colored according to the variable text colors in the script.
In Figure 2C, it can be observed that, even though the temperature is over the setpoint between 3 a.m. and 5 a.m. (Figure 2B), the PI signal is not zero. It only gets to zero when the integral part also decreases so much so that the sum of the integral and error parts is less or equal to zero. Although the OutSignal is limited, the negative values of OutSignalTemp are still used for the calculation. This enables the effect, which looks like prediction in some cases. This effect is further discussed in section 3.3.

Found Simplified Models
The simplified model of the system that is needed for the parameter calculation was estimated for 16 different cases. All three parameters of the gained models varied between all cases. The used cases and exact parameter values are included in Table 2 with parameter values also visualized in Figure 3. The process gain (K p ) has two clearly different orders and altogether three different levels. The values were around 1 for all cases where the PRBS signal was used as the setpoint and were much larger for  The T values varied least of the parameters, i.e., between 10,000 and 100,000 s (between around 4 and 15 h). Only in the same model 16 case, where an extra-large L value occurred, the T value was a lot lower at a bit less than 5000. So exceptionally, for this model, L is larger than T.
Based on mostly the Kp and L values, the models are divided into four groups, shown in Table  2. The setbacks and longer step groups are self-evident from above. The PRBS models are divided into models with a short L (PRBS sL) and a long L (PRBS lL). These groups will be used below for visualization.  Table 2.   Table 2.
The time delay (L) values for the PRBS cases had around a 100 times difference between the March week and February week values in R5, and the same difference was larger than 1000 times in R6, the southern room with more solar gains. L was smaller than 10 s for the two shortest setbacks, between 10 and 30 s for the February PRBS tests in R6, and larger than 100 in all other cases. There was Energies 2020, 13, 2068 9 of 20 ranging from 140 to 4000 (around 2 m to 1 h) and in one case (1-week PRBS in March for R6, model 16) it was over 50,000 s (around 14 h).
The T values varied least of the parameters, i.e., between 10,000 and 100,000 s (between around 4 and 15 h). Only in the same model 16 case, where an extra-large L value occurred, the T value was a lot lower at a bit less than 5000. So exceptionally, for this model, L is larger than T.
Based on mostly the K p and L values, the models are divided into four groups, shown in Table 2. The setbacks and longer step groups are self-evident from above. The PRBS models are divided into models with a short L (PRBS sL) and a long L (PRBS lL). These groups will be used below for visualization.

Identified PI Parameters
In total, 68 PI parameter value pairs were obtained. All the parameter values are included in Appendix A, Table A1. However, all the parameter combinations are also visualized in Figure 4a, where each point on the graph is a parameter combination. The scales are the logarithms of the parameter values with base 10. The graphs in Figures 4b and 5 follow the same logic. In Figure 4a, the parameter estimation method is shown by the marker shape and the model group is shown by the marker color. In the logarithmic scales, the tendency in the parameter estimation results seems to be roughly linear, so the lower the integration time the higher the proportional gain.
by the marker color. In the logarithmic scales, the tendency in the parameter estimation results seems to be roughly linear, so the lower the integration time the higher the proportional gain.
For the very small proportional gain, the integration time varies significantly from this otherwise linear behavior in the log10-log10 scale. The reason for this is depicted partly in Figure 5a. As can be seen, this covers the four cases calculated or optimized for March. Actually, these were all achieved for Room 6. This means that the solar peaks have been severe and almost no heating was needed. Therefore, these cases resulted in obscure parameters.
The clear separation between parameters is evident. The two sets of parameters with both blue and red (optimal) results made up one group and both green ones the other. This is also the difference in outdoor conditions, as can be seen in Figure 5a. The first group was generated at dynamic outdoor temperatures and realistic solar irradiation, while the second group bordered constant outdoor temperatures and no solar radiation. Here, also the separation between the March and Jan/Feb periods is clear, so it can be assumed that more solar gains causes the K parameter to be smaller and Ti to be longer. For the optimal cases, the combinations closer to the blue ones are optimized for the variable setpoint, the lower values for the constant setpoint.
In Figure 4b, the parameter combinations, which do not achieve the needed setpoints in Room 6 for at least 97% of the time (with a slack of 0.05 °C), are colored black. Both the one constant and two variable setpoint levels are checked and the coloring shows if any of the three are violated. If the graph would be for R5, all of the points, except the one with a dashed circle around it, would be black. This means that only one parameter combination would achieve the required temperatures in R5, if the setpoints were not shifted, as it was described in Section 2.6.
In Figure 5b, all the K-Ti pairs are colored by the log10 (K/Ti) value. This logarithm is further used for describing the pairs, as this is a clear indicator whether the pair is in the lower right or upper left corner of the log10-log10 graph.

Setpoint Temperature Tracking and PI Output Signal Behaviour
Each parameter combination results in different air temperature profiles and PI output signal profiles. There are four examples of the temperature and PI output profile combinations shown in Figure 6 for the constant setpoint cases and in Figure 7 for the variable setpoint cases in Room 6. In both figures, the Jan/Feb week is depicted on the left and the March week on the right. The parameter (b) (a) Figure 5. Graph (a) shows the underlying climate data and graph (b) shows the log-ratio values of all the PI parameter pairs. In (a), the constant climate is at 0 • C with no solar radiation, HP stands for heating period and all the dates are covered in Section 2.
For the very small proportional gain, the integration time varies significantly from this otherwise linear behavior in the log10-log10 scale. The reason for this is depicted partly in Figure 5a. As can be seen, this covers the four cases calculated or optimized for March. Actually, these were all achieved for Room 6. This means that the solar peaks have been severe and almost no heating was needed. Therefore, these cases resulted in obscure parameters.
The clear separation between parameters is evident. The two sets of parameters with both blue and red (optimal) results made up one group and both green ones the other. This is also the difference in outdoor conditions, as can be seen in Figure 5a. The first group was generated at dynamic outdoor temperatures and realistic solar irradiation, while the second group bordered constant outdoor temperatures and no solar radiation. Here, also the separation between the March and Jan/Feb periods is clear, so it can be assumed that more solar gains causes the K parameter to be smaller and T i to be longer. For the optimal cases, the combinations closer to the blue ones are optimized for the variable setpoint, the lower values for the constant setpoint.
In Figure 4b, the parameter combinations, which do not achieve the needed setpoints in Room 6 for at least 97% of the time (with a slack of 0.05 • C), are colored black. Both the one constant and two variable setpoint levels are checked and the coloring shows if any of the three are violated. If the graph would be for R5, all of the points, except the one with a dashed circle around it, would be black. This means that only one parameter combination would achieve the required temperatures in R5, if the setpoints were not shifted, as it was described in Section 2.6.
In Figure 5b, all the K-T i pairs are colored by the log10 (K/T i ) value. This logarithm is further used for describing the pairs, as this is a clear indicator whether the pair is in the lower right or upper left corner of the log10-log10 graph.

Setpoint Temperature Tracking and PI Output Signal Behaviour
Each parameter combination results in different air temperature profiles and PI output signal profiles. There are four examples of the temperature and PI output profile combinations shown in Figure 6 for the constant setpoint cases and in Figure 7 for the variable setpoint cases in Room 6. In both figures, the Jan/Feb week is depicted on the left and the March week on the right. The parameter combinations are chosen as the ones with minimum and maximum log10 ratios of the parameters, the IDA ICE default combination, and the one which resulted in optimal energy consumption (see Section 3.5). The combinations are ordered by the log10 ratio of the parameters with the minimum ratio at the top and the maximum ratio at the bottom. The IDA ICE default combination is the second (0.3/300) and the optimal is the third from the top (18/2300). Here, the parameter values were rounded to two significant numbers.
In the first column of Figure 6, most of the controllers show results that suggest maintaining a constant setpoint in the situation with no solar gains is an easy task. The small fluctuations are largest when a very small proportional gain (K = 0.012 in Figure 6a) with a large integration time is applied. This controller changes the signal too slowly, as its PI output signal in black shows. The signal stays almost constant throughout the day and even throughout the week. Due to the same effect, temperatures drop below the constant setpoint in March in Figure 6b and the setpoint tracking is poor in the variable cases. The level at which the signal is constant depends on the season, as there is a clear difference between February and March. abrupt the changes, as the cumulative graph indicates behaviors close to on-off signals. As shown in Figure 9, a zoom-in on R6's constant setpoint graph, the higher temperatures at the high-temperature end are clearly dependent on the Ti value. The low-temperature end seems to be more dependent on the K value. Therefore, the energy consumption of the parameter combination is mostly dependent on the K value and avoiding over-heating at the disturbances is more dependent on the Ti value. This effect was also observed in the analysis.   The constant setpoint cases in Figure 6 show that 2400/42 manages to maintain the constant setpoint the best. However, there is no significant difference for the variable setpoint cases. However, the PI output signal in the same case changes most rapidly. Both a large proportional gain and relatively small integration time contribute to this. Such switching reduces the life span of most of the devices, so this would not be acceptable in practice. For the case with also a large proportional gain but with a large integration time as well (18/2300), the signal is a bit smoother. In the long integration time cases, the heating starts earlier and stops sooner than for the shorter integration time. It can be observed that the PI signal turns on before the temperature lowers below the setpoint generating a prediction effect. This is especially clear for 18/2300 during the March week.
The variable setpoint cases in February in Figure 7's first column show that in cold weather with no solar peaks, the 24 • C setpoint peaks were not reached due to the short duration of the setpoint increase. Therefore, setpoint tracking during high setpoints is clearly not good but is also not required. However, the PI signal is 1 during these times, which means the heater is fully on as is the aim for load shifting. In this figure, again controllers 18/2300 and 2400/42 both maintain the lower setpoint well. However, the latter is switching on and off often and has almost no other state. In March, the solar peaks govern the temperatures. However, the second column of Figure 7 shows that the heating is turned on as well.
All the cumulative profiles over the heating period are shown in Figure 8. For the PI signal, only R6 is shown as the profiles look very similar for the two rooms. The switching behavior indicated before is clearly dependent on the log10 ratio of the PI parameters. The higher the ratio, the more abrupt the changes, as the cumulative graph indicates behaviors close to on-off signals. As shown in Figure 9, a zoom-in on R6's constant setpoint graph, the higher temperatures at the high-temperature end are clearly dependent on the Ti value. The low-temperature end seems to be more dependent on the K value. Therefore, the energy consumption of the parameter combination is mostly dependent on the K value and avoiding over-heating at the disturbances is more dependent on the Ti value. This effect was also observed in the analysis.

Setpoint Shifting
It is clear, that some of the parameter combinations did not achieve the required temperature setpoint and some resulted in higher temperatures above the setpoint. Especially at the high temperature end, there was also a clear difference between rooms R5 and R6, as can be seen from Figure 8. This was caused by the room orientations as the R6 faces south-west and gets more solar gains than the north-west orientated R5. As declared in Section 2.6, the setpoints were shifted for all cases in the way that temperatures would reach the required setpoint for at least 97% of the time. The shift values were different for R5 and R6 as well as for the constant and variable setpoint cases. As a result, all temperatures reached the given setpoints at around 95-97% of the heating period. This accuracy was considered satisfactory. The shifts are shown together with the energy consumption evaluation in Figure 10.

Setpoint Shifting
It is clear, that some of the parameter combinations did not achieve the required temperature setpoint and some resulted in higher temperatures above the setpoint. Especially at the high temperature end, there was also a clear difference between rooms R5 and R6, as can be seen from Figure 8. This was caused by the room orientations as the R6 faces south-west and gets more solar gains than the north-west orientated R5. As declared in Section 2.6, the setpoints were shifted for all cases in the way that temperatures would reach the required setpoint for at least 97% of the time. The shift values were different for R5 and R6 as well as for the constant and variable setpoint cases. As a result, all temperatures reached the given setpoints at around 95-97% of the heating period. This accuracy was considered satisfactory. The shifts are shown together with the energy consumption evaluation in Figure 10.

Setpoint Shifting
It is clear, that some of the parameter combinations did not achieve the required temperature setpoint and some resulted in higher temperatures above the setpoint. Especially at the high temperature end, there was also a clear difference between rooms R5 and R6, as can be seen from Figure 8. This was caused by the room orientations as the R6 faces south-west and gets more solar gains than the north-west orientated R5. As declared in Section 2.6, the setpoints were shifted for all cases in the way that temperatures would reach the required setpoint for at least 97% of the time. The shift values were different for R5 and R6 as well as for the constant and variable setpoint cases. As a result, all temperatures reached the given setpoints at around 95-97% of the heating period. This accuracy was considered satisfactory. The shifts are shown together with the energy consumption evaluation in Figure 10.

Energy Performance and Total Setpoint Tracking Accuracy
The energy consumption results after setpoint shifting are shown in Figures 10 and 11. It is clear that the variable setpoint cases consumed less energy. This is because the average room temperatures were lower. The setpoints were also higher than the constant cases in some periods but coincidentally the higher setpoint temperatures often occurred during the day and the lower setpoints occurred during the night, so this does not influence heating energy use much. Also, the high setpoints were not actually reached. In the constant temperature cases, a clear optimum emerged between the log10 ratio of −3 and −1. This means that in optimal cases, the K value was 10 to 1000 times smaller than Ti.
the on-off cases with different dead-bands. From top to bottom (yellow to blue) the corresponding dead-bands are 1 K, 0.5 K, 0.16 K, and 0.05 K. The optimal PI parameter combinations result in a lower energy consumption than even the lowest of the lines with an unrealistically small dead-band. The commonly used dead-band of 0.5 K consumes 2-3 kWh/m 2 /year more energy than the PI cases for the variable setpoint. For the constant setpoint, the lowest PI results are up to 7 kWh/m 2 /year or 9% lower than for the on-off with a 0.5 K dead-band, which, for example, in R6 is at 81 kWh/m 2 /year. Omitting the extreme poorly performing cases, the total variation in energy consumption is more than 10 kWh/m 2 /year or 12% in the constant setpoint case. Figure 11 shows the same data colored by the model group. The IDA ICE default parameter is at one edge of the optimum range with exactly 1000 times difference. The energy consumption is already around 5 kWh/m 2 /year or 5% higher on that edge compared with the optimal case. The parameters optimized for setpoint tracking are also close to an optimal energy consumption. The PRBS sL group performs well almost in all cases but not optimally, while in all other groups some combinations perform poorly. The optimal range of parameters is shown in detail in Table 3. Most of the optimal values were calculated using TRY climate data but the methods varied. Figure 11. Influence of the log10 of the PI parameters ratio K/Ti on energy consumption; colors visualize the underlying model group. The horizontal lines in Figure 10 represent the shifted energy performance at the benchmark for the on-off cases with different dead-bands. From top to bottom (yellow to blue) the corresponding dead-bands are 1 K, 0.5 K, 0.16 K, and 0.05 K. The optimal PI parameter combinations result in a lower energy consumption than even the lowest of the lines with an unrealistically small dead-band. The commonly used dead-band of 0.5 K consumes 2-3 kWh/m 2 /year more energy than the PI cases for the variable setpoint. For the constant setpoint, the lowest PI results are up to 7 kWh/m 2 /year or 9% lower than for the on-off with a 0.5 K dead-band, which, for example, in R6 is at 81 kWh/m 2 /year. Omitting the extreme poorly performing cases, the total variation in energy consumption is more than 10 kWh/m 2 /year or 12% in the constant setpoint case. Figure 11 shows the same data colored by the model group. The IDA ICE default parameter is at one edge of the optimum range with exactly 1000 times difference. The energy consumption is already around 5 kWh/m 2 /year or 5% higher on that edge compared with the optimal case. The parameters optimized for setpoint tracking are also close to an optimal energy consumption. The PRBS sL group performs well almost in all cases but not optimally, while in all other groups some combinations perform poorly. The optimal range of parameters is shown in detail in Table 3. Most of the optimal values were calculated using TRY climate data but the methods varied. The AAE of the temperatures for rooms R5 and R6 are shown in Figure 12. The AAE is clearly dependent on the room and setpoint but not on the parameter combination. The AAE is constantly at 0.5 K for R5 and around 0.7 for R6 in the variable setpoint cases. The accuracy here depends mostly on the solar gains. For the constant setpoint case, the optimal region is everything, with a Ti lower than 10 4 and a K higher than 10 0.5 . The error is around 0.2 K for all the simulations in R5, for R6 the error ranges from 0.25 to 0.6 K, and in extreme cases to 1 K.

Discussion
Different PI parameter estimation methods were applied on various periods and control profiles. An optimal region of the parameter ratio was determined where the energy consumption was the lowest. Half of the parameter combinations in the optimal region for energy consumption were found

Discussion
Different PI parameter estimation methods were applied on various periods and control profiles. An optimal region of the parameter ratio was determined where the energy consumption was the lowest. Half of the parameter combinations in the optimal region for energy consumption were found via GenOpt, although they were optimized for the minimal temperature setpoint tracking error. Although most reliably well-performing, this approach is not always suitable in practice as it requires an advanced model of the building. Therefore, it is practical that the other half of the parameter combinations in the optimal region were found using only short tests and simple calculations.
For all these other methods, simplified models were identified. In the optimal region, all the tested simplified methods were represented: Cohen-Coon, AMIGO, and SIMC. The results tuned in Matlab were not represented, probably due to the chosen goal being speed for that methodology. The models underlying these calculations were obtained from the week or weekend pseudo-random temperature setpoint (PRBS) data or setbacks of 6, 12, or 24 h. It is clear that the longer the setback, the easier it is to identify a simple model on it. This is probably the reason why the 1-and 3-h setbacks resulted in less desirable parameters. Still, conducting 24-h setbacks would probably not be comfortable for the occupants. Therefore, it is beneficial that 6-h setbacks could suffice. For example, these could be conducted during the night when the outdoor conditions are less variable with no solar gains. The suitable PRBS cases included both the January and March data, indicating that it is possible to get quality parameters in various weather conditions.
The optimal parameter combinations resulted in an annual heating energy reduction of up to 9% or 7 kWh/m 2 /year. The comparison of heat emitters and controllers in the European standard room shows similar results with 5% to 10% savings for the PI controlled UFH compared with the on-off control [20]. This does not compare to the 32% achieved for radiators in [23], however, the actual difference is difficult to compare as the baselines are different. The reduction of 7 kWh/m 2 /year here can be seen as highly significant as this can be achieved with only parameter correction, which does not require intensive computation when the simple tests are applied. Accounting for the more expensive thermostat head with variable parameters option, the payback time of this change is around 5 years. This saving can be achieved without setpoint reductions, which means no penalty on comfort. On the contrary, due to less fluctuation, comfort could even improve.
The methodology used here could be applied in any UFH system. In public and office buildings, a detailed model often exists and optimization of the parameters could be possible. Due to the large floor areas in these buildings, the absolute savings could be significant compared with the on-off control. Even more evident would be the saving in outdoor UFH systems installed under garage runways or stadiums to keep them clear from ice and snow.
Evidently, the parameter value results apply to the studied building, and future research can determine possible variation of the parameters in buildings with a smaller or higher thermal mass, insulation level, and maximum heating power. However, the wide range of well-performing parameter combinations and the fact that the suitable region is the same for both the north and south facing rooms provides an indication that the parameters from this region could be suited to different buildings as well. This should be confirmed by future studies on the subject.

Conclusions
Several combinations of the input data and PI parameter estimation methods were applied with the aim to improve UHF temperature control, resulting in 68 different PI parameter combinations. Based on the results and discussion above, most importantly concluded is that: • For the first time in the scientific literature, it is shown that UFH can operate with determined PI parameters similar to ideal control; • A performance close to optimal could also be achieved by parameters achieved from shorter tests, e.g., weekend pseudo-random setpoints, and 6-to 24-h setbacks which were shown to be suitable; • The optimal PI parameters improved the room temperature control accuracy considerably, and that the results show that the UFH PI control with the correct parameters started to work in a predictive fashion and the resulting room temperature curves were practically ideal; • The optimal PI parameters reduced the energy consumption for heating by up to 9% (7 kWh/m 2 /year) in comparison with the on-off control (at around 80 kWh/m 2 /year) and by 5% in comparison with the default PI parameters; • The variation amplitude of the heating energy needed using different estimated (not random) parameters was more than 15 kWh/m 2 /year for the constant setpoint, which stresses the importance of having the correct PI parameters; • The optimal PI parameters included combinations with log10 (K/Ti) between −3 and −1, in these combinations, the proportional gain K ranged from 2 to 100 and the integration time Ti from 500 to 6700 s, and thus higher gain and longer integration time values than are conventionally used are recommended; • For the variable setpoint, using the PI control had a similar effect to decreasing the dead-band and the variation in the PI parameters did not have a significant further effect on the energy consumption, except for when they were extremely poorly tuned; • The average absolute error for the air temperatures from the setpoint was well below 0.5 K for the constant setpoints, but above for the variable setpoints.