Performance comparison of quadratic, nonlinear, and mixed integer nonlinear MPC formulations and solvers on an air source heat pump hydronic floor heating system

There is a gap in literature on comparisons between different MPC optimal control formulations and solver choices for the same building HVAC system. Mixed Integer Nonlinear (MINL) formulations are rarely considered, despite being the most physically accurate way to represent HVAC systems. This work compares several MPC formulations, including Quadratic, Nonlinear, and MINL, applied to a case study building and investigates benefits and challenges of MINL MPCs from practical perspectives. Ten different MPC formulations were developed and implemented using Pyomo. Then, a detailed emulator model was developed using open-source Modelica libraries and used with BOPTEST to assess the performance of each MPC. Results show that convergence and control switching behaviours of MINL MPCs are sensitive to formulations, initialization approaches, solver selections, and solver parameters. Thus, they require significant effort for tuning. However, a very well-tuned MINL MPC performed similarly to successful Nonlinear MPC formulations.


Introduction
HVAC systems account for 20% of the total primary energy consumption in Major Economies Forum (MEF) countries (Metrics 2015). Therefore, advanced controls for those systems can help reduce environmental impact as well as help renewable penetration by unlocking load flexibility potential in buildings (Roth et al. 2002;del Mar Castilla et al. 2014). For the last decades there has been a lot of research on advanced control, including Model Predictive Control (MPC), for the optimal operation of building HVAC systems (Drgoňa et al. 2020;Kathirgamanathan et al. 2021;Rockett and Hathway 2017).
For MPC, the optimal control problem can be formulated mathematically in a variety of different ways, even for identical HVAC systems. The primary reasons for the variety include: • There are various approximation techniques, e.g. the McCormick convex relaxation of a bilinear function (McCormick 1976), piecewise linear approximation for a nonlinear function, and linear programming (LP) relaxation of a mixed integer linear program (MILP), which relaxes the integer constraints. • MPC itself has different theoretical approaches, such as centralized versus decentralized/distributed or stochastic versus deterministic. • The objective and constraints of a MPC can be formulated in different ways regarding the quantification, relative importance, and constraints on energy use, carbon emissions, energy cost, load flexibility, and thermal comfort. • Each component of an HVAC system can be modelled in several ways, which affects the performance of the prediction accuracy, robustness and computing time. For example, modelling the Coefficient of Performance (COP) of an air to water heat pump ranges from the constant COP approach to the DOE-2 (Research Group LBL, S 1991) like performance mapping as a function of part load ratio, outdoor air temperature and supply water temperature.
• The selection of MPC optimization variables is a control design factor. For example, one may select either controllable inputs that are the same as those of the physical system, such as percent valve position, or that are abstracted from the system model, such as heat flow rates. In the latter case, a strategy is needed to convert optimal solutions to control inputs that are available in the physical system.
The significance of these variations of MPC formulations is that they could change not only the accuracy/physical reliability of a model and computation time, but also the class of optimization problems (e.g. from convex to nonconvex, from NLP to LP or MILP, and vice versa), affecting mathematical properties of global optimality, feasibility, uniqueness of a solution, convergence of numerical optimization algorithms, and MPC closedloop stability. Consequently, it could considerably impact the overall MPC performance (e.g. energy consumption, comfort, computational time, and the rate of change of control inputs). As a consequence, a significant amount of time could be spent iterating to identify the most suitable optimization formulation for each HVAC system. Despite this importance, there are very few well-documented papers that investigate different MPC formulations. This is especially true for non convex MINL MPC, which is one of the most natural and straightforward ways to formulate optimal control problems for many HVAC applications. This paper tries to partially fill this gap in the literature by providing comparisons of multiple MPC formulations for an HVAC system that has both binary and continuous control variables. In addition, this paper investigates the practical applicability of MINL MPC approaches, but not MILP and more in general Mixed Integer Convex Programming (MICP) problems. The reasons for excluding MICP are twofold. First Mixed Integer Linear (MIL) and Mixed Integer Convex (MIC) problems can be solved reliably with commercial and open source tools. Second, there are several ways to approximate an MINLP into an MICP or MILP, such as piecewise linearization or McComik methods. Therefore, to avoid dispersing the manuscript focus, a separate work should be carried out to test the different approximation methods on building HVAC control related problems.
Section 2 reports the case study details, MPC formulations, optimization solvers and the co-simulation setup. Section 3 presents the results for the MPC formulations comparison. Finally, Section 4 reports the main conclusions of this study.

Literature review
Equivalent formulations for building applications can be found in many papers. One example is when peak demand is considered in a control objective. This approach replaces the maximum power or demand cost over a prediction horizon with linear constraints and a slack variable, namely the target peak or demand cost (ASHRAE 2019, Chapter 43). Depending on applications, this approach could convert NLP to LP or MINLP to MILP .
Approximation techniques have also been widely applied to HVAC systems. Risbeck et al. (2017) applied a piecewise linearization technique for optimal scheduling of operations of chillers, pumps, cooling towers, boilers and thermal energy storages. The introduced approach approximates a nonlinear chiller performance map with a set of piecewise linear models, converting MINLP to MILP. Kim et al. (2015) applied a linear programming relaxation approach that relaxes integer constraints for coordinating operations of multiple rooftop units. This approach converts an IP to LP for better computation efficiency. Atam and Helsen (2015) applied a convex relaxation technique to handle the bilinearity that naturally appears in thermal energy systems. The proposed method converts a nonconvex optimization problem to a convex problem to ensure global optimality.
Considering MPC architecture and the theoretical approach, Scherer et al. (2014) and Walker et al. (2017) compared centralized and distributed MPC architectures, highlighting that distributed approaches have slightly worse KPI performance but better computational time. Oldewurtel et al. (2012), Drgoňa et al. (2013), Ma, Matuško, and Borrelli (2014) and Maasoumy et al. (2014) compared deterministic versus robust or stochastic MPC, showing that a robust or stochastic MPC performs better in scenarios of high uncertainty and is comparable in other cases. Rather than focussing on architecture and theoretical approach, Cigler et al. (2013) and Drgona and Kvasnica (2013) instead analysed the formulation of the MPC problem, focusing on different cost functions and constraints, assessing which formulations are more robust and computationally efficient, but limiting their analysis to LP, QP and MILP.
Considerable effort was also put into analysing different building envelope thermal modelling approaches. Prívara et al. (2013) compared several black-box and greybox model structures to model building envelope systems and concluded that black box models are more computationally efficient for larger case studies but become less reliable for longer prediction horizons. Sourbron, Verhelst, and Helsen (2013) analysed the effects of grey box model order on the performance of MPC for concrete core activated buildings. Blum et al. (2019) also shows that model order has a strong influence on the model quality. Furthermore, Blum et al. (2019) identified seven factors that play an important role in the accuracy of the building envelope model. Picard et al. (2017) and Picard et al. (2016) show that a purely physical driven white box approach can be viable in certain building types. Kim et al. (2016 pointed out that a typical identification algorithm with any model structure would likely result in a biased model when significant unmeasured disturbances (e.g. unmeasured internal heat gain, in/exfiltration, door/window openings, zonal plug load) presented in a training dataset, and proposed a new identification approach for a typical grey box model structure to mitigate this negative effect. Drgoňa et al. (2020) and Serale et al. (2018) provided comprehensive reviews on building MPC literature. For the papers that they reviewed, many works address the benefit and applicability of their own MPC formulations compared with rule-based controls. Less comprehensive work is available on comparisons between different HVAC modelling approaches, optimization variable choices, and their impacts on the resulting MPC performance. Verhelst et al. (2012) performed an extensive analysis of different COP formulations in the MPC problem leading to LP and NLP problems, highlighting the potential benefit of a nonlinear formulation. Pčolka et al. (2016) compares a linear time invariant MPC, a linear time variant and a nonlinear MPC in a case study for a heat pump and domestic hot water system. It reports that the nonlinear solution is the best, but the linear time variant gets close and remains more robust. In both studies of Pčolka et al. (2016) and Verhelst et al. (2012), binary variables were not taken into account to avoid MILNP formulations, although MINLP arises fairly naturally when dealing with HVAC systems. Indeed, very few studies can be found comparing MINLP with other formulations. To the authors' knowledge, only Burger et al. (2018) introduces a custom MINLP solver compared with Bonami et al. (2008) for a solar thermal system. Furthermore, it is hard to crosscompare to different works due to the unique case study system and lack of common metrics.

Objectives and contributions
This work aims to partially fill this literature gap by presenting comparisons of ten MPC formulations with a greater focus on MINL MPCs, for a relatively common HVAC system that requires control decisions on valve on/off and supply water temperature setpoint. The diversity of formulations is due to two issues that appear in a broad range of HVAC systems: (1) the nonlinearity arising from modelling heat pump COP and (2) the binary onoff physical control inputs for distribution circuit valves. Depending on modelling and approximation approaches to handle them, the resulting optimal control problem formulations become QP, NLP or MINLP. Each formulation encompasses a trade-off between accuracy in the prediction, robustness to find an optimal solution, computational requirements, the rate of change of control inputs, energy consumption, and comfort violation.
The contributions of this paper are: • Present a well-documented work on how HVAC performance can vary with MPC formulations • Understand the benefits of increased prediction accuracy from increased model complexity and the corresponding trade-offs • Introduce a new Key Performance Indicator (KPI) to quantify the rate of change of a control input • Survey, introduce, and test available optimization solvers of each problem formulation, especially novel MINLP-specific solvers, that were not comprehensively investigated, but could potentially be useful for typical building HVAC optimal control problems • Share lessons-learned for designing MINL MPC.

Case study description
The chosen case study is a newly built two-room apartment in Milan, Italy. The HVAC system is a two-circuit radiant floor heating system connected to an air source heat pump. A diagram of the HVAC system is presented in Figure 1.
Each thermal zone is independently controlled via its own on/off valve. The pump works with a constant head tuned to provide the nominal water mass flow rate to each floor heating circuit when a valve is open. Control decisions are the two valves' statuses and the heat pump supply water set point. Despite its simplicity, the case study includes the HVAC hydronic system components that allow for several optimal control models, and lead to three different optimization problem classes to be solved: QP, NLP, and MINLP.
The emulator model for the apartment and HVAC system was developed in Modelica using the IBPSA 3.0 (master branch commit 8a0d237) (Wetter, Blum, and Hu 2019), Buildings 8.0 (master branch commit 69bb7cf) (Wetter et al. 2014) and IDEAS 2.2.1 (master branch commit 32860ea) (Jorissen et al. 2018) libraries. The MPC algorithms were implemented using the Python-based optimization modelling language Pyomo (Bynum et al. 2021). Finally the various MPC formulations were co-simulated with the emulator model thanks to the Application Performance Interfaces (APIs) and run-time environment provided by the BOPTEST software framework (Blum et al. 2021), which also provides as output a standard subset of KPIs, including thermal discomfort, energy consumption, cost of the energy and computational time ratio. In this way, it will be possible to consistently compare all MPC formulations on the same emulator, highlight the pros and cons of each approach, and make the emulator publicly available for continued usage and further comparison of control approaches.

Introduction of optimization solvers
Different solvers and not different algorithms were compared on the same platform because comparing the solvers can be considered the state of the art of available algorithmic implementations, while a custom implementation could be criticized. Furthermore, great care was taken in choosing the optimization and co-simulation toolchain to make the comparison as fair as possible. Pyomo allows to easily couple different solvers through the AMPL interface (AMPL 2003) which is supported by a wide variety of solvers. Table 1 summarizes all the solver used for the study.
The QP and NLP solver, IPOPT (Wachter 2002) was chosen due to the popularity and widespread usage. IPOPT uses MA57 as linear subsolver. The MINLP solvers were chosen by looking at the results from Kronqvist et al. (2019), which analysed solver performance on a set of 335 convex MINLP problems and included both open source and commercial solvers. The subset of solvers chosen for this study includes some of the best performing and most popular choices for open source and commercial alternatives. BONMIN (Bonami et al. 2008) is an open source project belonging to the project COIN-OR foundation (COIN-OR Foundation 2006), as does IPOPT, and it is a MINLP solver mainly used for convex MINLP problems. MINLP solvers divide the optimization problem into a MILP or MIQP problem and NLP sub problems. In our setup, BONMIN uses CBC (Forrest and Lougee-Heimer 2005) as the MIQP subsolver and IPOPT as the NLP solver. In particular, two algorithms are tested from BON-MIN. The first is BONMIN-BB, which uses a variation of the Branch and Bound algorithm to convert the problem. The second one is BONMIN-Hyb, which is a hybrid approach between Branch and Cut and Outer Approximation algorithms, which is faster than BONMIN-BB, but suffers more from the issue of falling into a local minima when the objective function is not convex. As a commercial alternative, Baron (Kılınç and Sahinidis 2018) was used since, differently from BONMIN, it should be able to guarantee a close to global optimum even when the objective function is not convex. Baron uses a variation of the Branch and Bound algorithm to reduce the MINLP problem into a subset of NLP problems and a MIQP problem. The MIQP solver is CPLEX and the NLP solver is IPOPT. Lastly an open source alternative to Baron, SCIP, is also tested, which is a global solver that uses a variation of the Branch and Bound algorithm. In our setup, SCIP uses CBC as the MIQP subsolver and IPOPT as the NLP solver.

Emulator model
In Figure 2, a schematic of the apartment is presented, while in Figures 3-5, the yearly frequency plot of the dry bulb temperature, humidity ratio and global horizontal radiation for the location (Milan, Italy) are shown. Milan can be considered a continental temperate humid climate. The maximum T DryBulb is 32 [ • C], the minimum is −7.4 [ • C] and the average 11.7 [ • C]. In Figure 6, a brief validation of the emulator model is shown, which used experimental data coming from a globe thermometer positioned in the centre of the living room, while the boundary conditions were determined from local weather stations and localized forecast services. For an in depth description of the envelope, the reader can access the   test case (two zone hydronic apartment) documentation at the BOPTEST repository Two zone hydronic apartment (IBPSA 2019). In Table 2, a summary of the main features of the case study is reported. The thermal zones and floor heating system are modelled using the Buildings library. The remainder of the hydronic system is modelled using the IBPSA library apart from the heat pump, where a dynamic performance map model from the IDEAS library is used. The baseline controller in the emulator is an on-off controller with a 1[ • C] hysteresis on the room set point temperature. The zone valve fully opens when the hysteresis controller gives an on signal and remains closed otherwise. Each thermal zone has its own thermostat and is controlled independently. The pump works to provide a constant head, tuned to provide each floor heating circuit the respective nominal water mass flow rate. The heat pump supply water set point temperature is calculated via a climatic curve that depends on the external temperature. The  baseline results shown in the Results Section 3, refer to the emulator running the simulation with this baseline controller. The external Python-Pyomo based MPC is able to override both the zone valve on-off signal, u i , and the heat pump supply water set point T in,set .

Reduced order model
A grey-box model based on the resistance-capacitance (RC) analogy was identified using the Matlab identification toolbox (Ljung and Singh 2012) for use within the MPC controller. Looking at Table 2, the apartment can be considered to be well insulated and with heavy construction. Different combinations of resistors and capacitors were tried, leading to a 3C7R scheme that was adopted for each thermal zone. The scheme of the resulting RC circuit for each thermal zone is shown in Figure 7. The three capacities are related to the room temperature T r , wall temperature T w and floor temperature T f . Resistances  connect the capacities nodes to each other and furthermore, two resistances connect C r and C w to the external temperature T ext . The wall has a resistance that connects also with the sky temperature T sky . The sky temperature allows the low-order model to better treat the radiative heat exchange with the external environment, especially in the presence or absence of clouds. Lastly, the capacitors of the rooms are connected to each other through a resistor as a proxy for air exchange between the two thermal zones. The solar heat source s [W/m 2 ] is the hemispherical global radiation hitting the external wall and window. It is divided between the wall and the floor and multiplied respectively by the opaque area A wall and the windows area A win . a and c are tuning parameters that can be assumed as proxy of absorptance and trasmittance.
int [W] are the internal gains divided between sensible and radiative by the parameter b. The sensible part goes to the room C r and radiative goes to the wall C w . Finally, the heat flow rate to the floor heating system, injected in the floor capacity C f , is shown in Figure 7 as Ḣ = H in −Ḣ out . It is modelled in four different ways resulting in different classes of optimization problems: (1) Linear formulation where Ḣ itself is treated as an optimization variable (in this case, we let Ḣ :=Q).
(2) Linear formulation where Ḣ is modelled with the supply water temperature T in , return water temperature T out and the nominal value of mass flowṁ fnom (in this case Ḣ :=ṁ fnom c w (T in − T out )) Figure 7. Thermal zone reduced order model, the red dots are the temperatures, the blue parallel lines are the capacitors associated with the temperature states, the resistances are the thermal resistances between the temperatures and the red lines indicate a heat flow into the node.
(3) Nonlinear formulation where Ḣ is modelled with the energy balance in 2., but multiplied by the valve position u i , and with the continuous relaxation of the integer constraint, i.e. u i ∈ [0, 1] (in this case Ḣ := u iṁfnom c w (T in − T out )) (4) Mixed Integer Nonlinear formulation where Ḣ is modelled as 3. but without the continuous relaxation.
For formulation 1 the optimal control variable will be the thermal power provided to the floor heatingQ. For formulation 2 the optimal control variables will be the floor heating inlet temperature T in . For formulations 3-4 the optimal control variables are T in and the zone valve status u i . For the formulations 2-4 using the supply temperature T in as an optimization variable, the return/outlet temperature T out was modelled with the following linear equation to correlate T out with T in and floor temperature T f as shown in Equation (1).
w f is the weighting factor for the identification process. The rationale behind this linear relationship is that the water mass flow rateṁ floor is constant and so if we consider the floor heating system as an heat exchanger w f would be equivalent to a constant effectiveness. This equation is valid for this case study since the zones valves can only be open or closed and do not provide variable flow control. The authors would like to point out that this approximation might be more problematic for formulation 3., because with the continuous relaxation the valve will be able to modulate the flow and could lead to a larger error between the reduced order model and the emulator model and instability in the MPC.
Since the training of RC models is not the focus of the present work, the identification procedure is briefly reported. All the reduced order model parameters shown in Figure 7, so not including w f , were trained using two weeks of free floating data, simulated using the detailed emulator model, where the boundary conditions were derived from a synthetic profile obtained through a Fourier analysis of the typical year data (Smart and Ballinger 1984). In this way all major frequency components are present. After identifying all the model parameters, an additional week of data where the heating system is excited by turning it on and off is used to find the weighting parameter w f . The result of the overall identification process leads to a Normalized Root Mean Square Error (NRMSE) goodness of fit of 82%, defined in Equation (2), in open loop simulation for the whole heating season which for Milan is from the 15th of October to the 15th of April, where the heating system was functional.

MPC formulations
Besides the modelling approaches, other constraints and objective functions were specified to complete MPC problem formulations. In Table 3, all formulation elements and MPC implementation parameters that were commonly used for the ten different MPCs are shown. Furthermore, Table 4 lists all other optimization variables, constraints and objectives. Finally, in Table 5, complete MPC formulations are succinctly summarized with the notations of Table 4. Table 3 contains all common elements across the formulations, including dynamic states, disturbances and implementation choices. In Table 3, only the final choice for control horizon and time step is shown. However, four different control horizons were tried for the optimal control problems from 6 up to 72 [h]. The final choice was 24 [h] since a longer prediction horizon did not show significant improvement on KPIs. This time scale aligns with the fact that we are dealing with a heavy construction and a floor heating with a high thermal inertia. Another parametric study was carried out to identify a suitable time step where a MPC solution updates. Several different time steps were tested ranging 5,10,15 and 30 [min], the results show that bringing it below 15 [min] did not give any significant benefits.The reason is that in this case study the authors considered a floor heating system where the time constant of the floor heating is in the order of magnitude of 3 [h], in case of air systems or radiators the time step should be lower to avoid numerical integration problems. A direct collocation method was implemented Energy cost NLP using the Pyomo problem statement to convert the continuous optimal control problems into discrete programming problems. The six states include all temperatures in both thermal zones and the disturbances as reported in Section 2.4, plus the two room set points and the energy price p e [e] set to 0.20 [e/kWh]. Table 4 lists all other optimization variables, constraints, and objective functions used in different formulations (each MPC formulation is expressed with a combination of those in Table 5). δ liv/bed is an auxiliary variable representing a temperature deviation from a setpoint, and is coupled with the constraint C comf and the objective j comf . By looking at the constraint C comf , the value of δ will be higher than zero if the room temperature T r is lower than the setpoint temperature T set . In this case, it will be penalized by including δ 2 in the objective j comf . This will push the MPC to keep T r higher than the setpoint temperature.
The other optimization variables, as well as the constraints, are related to the floor heat flow rate models explained in Section 2.4. T in is the supply temperature and can go up to the maximum temperature of 45 [ • C] to avoid high temperatures in the floor, down to a minimum temperature, defined as the adiabatic mixing temperature in constraint C Tmin . The formulation of constraint C Tmin comes from a local energy and mass balance at the return outlet of the floor heating system under the assumption that the nominal flow rate is the same for all circuits. T f ,liv/bed is the floor temperature and u liv/bed is the floor heating circuit valve control. u R represents a continuous relaxation on the valve control, so that the valve can continuously modulate the flow from totally closed to totally opened. u B indicates the valve controls without the relaxation and are consistent with the actual system and emulator model. u open means the valve controls under the assumption that the valve is always open and the modulation is carried out only at the supply temperature T in . Finally,Q liv/bed is used directly as an optimization variable when the heat flow rate is not modelled explicitly. It can go from zero to a maximum value determined by constraint C Q max . The constraint on the heat flow rate C Q max imposes the maximum heat flow rate linearly with the floor temperature as a function of the maximum supply temperature T in max and the weighting parameter w f . This constraint helps to model the behaviour of the floor slab, where, at constant supply temperature, the higher the floor temperature, the lower the heat transfer rate to the floor. Figure 8 is given as a visual representation of the idea.
The last section of Table 5 presents all the components that can make up a complete control objective function. The complete objective function J tot to be minimized is the sum of these different objective components with some weights, denoted as k i where i corresponds to a specific objective component. The weighting parameters k i need to be tuned to balance the impact of each objective  on the total objective function J tot . To find the best values of k i , several iterative studies were performed for each MPC formulation and priority was given to the comfort constraint. The objectives j COP,L and j COP,NL are the energy cost and are calculated as the energy price p e multiplied by the total heat flow rate provided by the heat pump, Ḣ liv (t) + Ḣ bed (t), divided by the heat pump COP. The COP is a function of the external temperature COP(T ext ) for j COP,L , so the underlying problem remains quadratic.
In j COP,NL , the COP is a function of the external and supply temperatures COP(T ext , T in ). Here, the optimization problem becomes nonlinear because a control variable is present in the denominator of a fraction. j comf is the temperature mismatch between room temperature T r and setpoint T set and works as a comfort proxy. The switching frequency objective j switch is the sum of the squared valve control derivatives. This objectives serves the purpose of penalizing undesirable sudden changes in the control variables. Since the derivative is discretized with forward Euler in the direct collocation method, there is no difference between the formulations that use valve control state as integer u B and as continuous u R . Finally, the binary constraint, j B , forces u liv and u bed to be close to either 0 or 1 to avoid having an objective greater than 0. The reasoning behind this constraint is to approximate a MINLP as a NLP. However, as was found in this study, care should be taken when initializing this optimization problem because the introduction of the binary constraint causes a significant discontinuity in the solution space. Starting from Tables 3 and 4, several MPC formulations can be defined, ranging from QP to MINLP. Table 5 reports the MPC formulations coupled with the solvers and relative options used for the study. Table 5 shows the ten MPC formulations that were tested. The Tag column reports the formulation names. The Formulation column shows the corresponding optimization variables, constraints, and objectives described in Table 4. Note that the auxiliary variable δ, the temperature constraint C comf , the temperature mismatch j comf , and the switching objective j switch are present in all formulations, apart from MPC2 where j switch is not needed since there is no valve control. The Problem Type column reports the optimization problem type for each MPC formulation, which can be either QP, NLP or MINLP. The Solver column shows the solver chosen, including the MINLP handling algorithm option if present. The Tolerance column is the solver tolerance which was determined through a parametric study as a compromise between quality of the solution and computational time for each solver. The same tolerance must be used for the same solver class to guarantee fairness in the solution and computational time. Furthermore, it is increased for the MINLP solvers since it considers both the reduction in the cost function steps for the NLP sub problems and the integer approximation for the MINLP problem. The Initialization column defines how the optimization problem variables were initialized. Free-floating initialization means that a simulation runs using the reduced order model subject to the same boundary conditions, as in the forecasts used for the optimal control, with the floor heat flow rate set to zero. In other formulations, a slightly randomized solution of a different MPC formulation is used as initialization as indicated. The Post Process column shows the steps needed to convert the optimal control trajectory into the physical control inputs used in the emulator model. The Subsolvers column refers to the solvers used by the solver indicated in the Solver column. Another parameter not included in the table is the Timeout. It corresponds to the maximum time between each control horizon optimization and was fixed at two minutes for all formulations as a compromise between giving the solvers enough time to converge to a solution and the overall computational time.
Below, a summary for each MPC formulation presented in Table 5 is provided: • MPC1: This formulation uses heat flow rates directly as optimization variablesQ liv/bed , and j COP,L as the energy objective. In this way, the final constraints are linear in optimization variables, and the objective function is quadratic, making a QP problem. The solver of choice for QP was IPOPT with the MA57 linear subsolver. The tolerance was set to 10 −6 from the default value of 10 −8 and the timeout time to 120 [s], and the initialization is free-floating. Some post processing is required to convert the optimal control trajectory into the physical control variables used in the emulator. IfQ liv/bed is higher than a threshold value, equivalent to the minimum cutoff power of the heat pump set as 20% of the nominal value 800 [W], the valves u liv/bed will be opened otherwise they remain closed. The supply temperature setpoint is calculated using the previous step return temperature plus the delta given byQ liv/bed .  (3) in Section 2.4 and in Table 4. The final problem formulation is nonlinear in terms of constraints and optimization variables due to the multiplication between supply temperature and valve control. The objective is also nonlinear due to the presence of COP as a function of the supply temperature T in and the multiplication of two optimization variables in the calculation of the heat flow rate Ḣ . The added nonlinearity of MPC3 compared to MPC2 makes the problem nonconvex because the solver can change T in or u R to modulate the heat flow rate, making the process of finding a global optimum harder. Solver settings were identical to MPC2. MPC also needs to convert the optimal control trajectory into the physical control variables used in the emulator. If the circuit valve u liv/bed value is higher than a threshold value, the most common choice would be 0.4. However, the authors found out that playing around with this parameter and reducing it to 0.3 lead to a more stable solution. So in case the valve control is higher than the threshold the valves u liv/bed will be opened, else they remain closed.
Then T in is used as setpoint for the supply temperature in the emulator.
• MPC4: This formulation is identical to MPC3, apart from the addition of the binary objective j B . The idea behind j B is to force u R to be either closed u R = 0, or open u R = 1, by penalizing all solutions that modulate the flow rate in a continuous manner. This allows for a smaller control space by reducing the optimal operating range of the zone valves. This will serve as an NLP approximation of a MINLP formulation. In this case u liv/bed will be used directly in the emulator, because u liv/bed in the raw solution when using j B are very close to 0 or very close to 1, making the rounding process trivial and T in is used as setpoint for the supply temperature. • MPC5: This formulation is the same as MPC3 but replacing the continuous relaxation (u R ) with the binary constraint (u B ). In this way the final problem formulation is mixed integer nonlinear due to the multiplication between contnious supply temperature and on/off valve control variables, making a MINLP problem. The additional complexity of MPC5 compared to the previous QP and NLP formulations requires a dedicated MINLP solver. The solver of choice for MPC5 was BONMIN-BB. CBC was used as MIQP subsolver and IPOPT as NLP solver. The tolerance was set to 10 −4 from the default value of 10 −6 and the timeout time to 120 [s]. The initialization is done by taking the solution of MPC3 after rounding the values of u liv/bed . The MINLP solution can be directly applied to the emulator with no need for post processing of the solution.In this case u liv/bed will be directly used in the emulator and T in is used as setpoint for the supply temperature. • MPC6: This formulation is identical to MPC5 apart from the addition of the binary objective j B and the initialization done with MPC4 solution. The rationale is similar as in the transition from MPC3 to MPC4. However, instead of a single NLP problem, it is extended to all subsets of NLP problems generated by the MINLP solver.In this case u liv/bed will be directly used in the emulator and T in is used as setpoint for the supply temperature. • MPC7: This formulation is identical to MPC6, where a different MINLP algorithm option was used for the Bonmin solver BONMIN-Hyb. In this case u liv/bed will be directly used in the emulator and T in is used as setpoint for the supply temperature. • MPC8: This formulation is identical to MPC5, where a different MINLP solver was used named Baron. The MIQP solver is CPLEX and the NLP solver is IPOPT.In this case u liv/bed will be directly used in the emulator and T in is used as setpoint for the supply temperature. • MPC9: This formulation is identical to MPC8 with the addition of the binary constraint j B .In this case u liv/bed will be directly used in the emulator and T in is used as setpoint for the supply temperature.
• MPC10: This formulation is identical to MPC9 with the difference that the MINLP solver of choice was SCIP.In this case u liv/bed will be directly used in the emulator and T in is used as setpoint for the supply temperature.

Co-simulation setup
All elements shown in previous sections are coupled in a co-simulation, a graphical representation is given in Figure 9. The optimal control routine runs on the Pyomo Python toolbox Pyomo (Bynum et al. 2021) and the detailed emulator model is wrapped in a Docker container using the BOPTEST software (Blum et al. 2021). all simulations were carried out on a Linux Ubuntu 18 laptop with 16GB of RAM and an Intel(R) Core(TM) i7-8650U CPU @ 1.90GHz. all solvers have multi thread capability so up to 8 threads were used for the simulations. This does not mean that all solvers use multi-thread with equal effectiveness. Different solvers will have different core-allocation mechanisms which may result in very different multi-thread efficiency. However, this is also part of the evaluation in the performance of a solver, since recent core architectures make parallel computing very important. All cases mentioned in Table 5 were directly implemented in Python using a concrete instance modelling feature of Pyomo. The solvers were compiled externally and coupled with Pyomo using the AMPL interface. Finally the Kalman filter from Labbe (2018) was used to update the states of the reduced order model at each time initialization.
BOPTEST provides an easy to use API interface that allows the optimization scripts to manipulate the control variables of the detailed model, access sensor data and access forecasts, and the KPIs calculated by BOPTEST. The control variables are the supply temperature setpoint and the zone valve open or closed signal. The forecasts are the disturbance variables reported in Table 3 and are considered deterministic, which means that the same disturbances will also be used in the emulator model. The measurements are the room temperature, the water supply temperature and the return temperature from each zone, and are used by the Kalman filter to estimate the initial value of room, wall, and floor temperatures for each zone for each prediction horizon. To estimate the performance of each MPC formulation, BOPTEST can calculate many KPIs, though this paper considers the thermal discomfort, computational time ratio, and energy cost. Furthermore, four additional KPIs were used to evaluate the performance of these MPC formulations, namely the total computational time [s] of the MPC solver, the thermal energy used [kWh], the control arc length, and the number of MPC solver time-out or error events [%] that occurred case study specific throughout the evaluation period. The equations and descriptions of the KPIs are reported in Table 6, and the descriptions of their variables are reported below the table. For a fair KPI estimation, the emulator and reduced order model thermal envelope states were initialized to 20[ • C] temperature, and a warm up period of 5 [days] was used. Only valid solutions are used to update the MPC control trajectory. Time-out solutions are valid when variables values are within the bounds, however, they may be not fully converged, meaning that constraints may not be completely satisfied. Time-out solutions are considered invalid if outside variables bounds. Solutions with an error in the solver status are discarded. When a solution is discarded, the previous solution current time step is used until a new valid solution is found. In the discomfort KPI K dis calculation, T r [K] is the room operative temperature and T r,set [K] the heating setpoint, N is the total number of zones, and t o and t f are the initial and final time of the evaluation period. In the computational time ratio KPI K timr δT k [s] is the MPC computational time at step k, M is the number of control steps and δt k [s] is the time interval of control step k. In the cost KPI K cost , p el [EUR/kWh] is the electricity price considered as constant, Q tot [kWh] is the total energy supplied by the heat pump and A tot [m 2 ] is the total floor area of the apartment. Q tot [kWh] is also used to calculate the thermal energy KPI K en . K err is the ratio between the number of MPC iterations that had either timeout N timeouts or errors N errors and the total number of iterations N total , 2976, for the period considered, the month of January with a time step of 15 [min]. Lastly, a new KPI is introduced in this manuscript called control arc length K conlen . The idea of this KPI is to quantify the frequency of switching for the control system throughout the evaluation period. K conlen is the ratio between the length of the actual control trajectory. It is calculated on the total heat flow rate provided by the heat pumpQ hp versus a fictional reference trajectory u ref that consider the control variable u to be constant for the evaluation period. A visual representation is given in Figure 10.

Results
Simulation evaluation results for the month of January are shown because similar conclusions were drawn from the rest of the heating season. As mentioned, all MPCs were updated for every 15 [min] time step with the control horizon of 24 [h]. Each MPC was calculated 2976 times for the simulation period and has around 1000 constraints and 600 optimization variables for the linear problems and 1200 optimization variables for the nonlinear problems.

KPIs comparison
Figures 11 and 12 report the KPI results calculated by BOPTEST. From the Discomfort KPI K dis in Figure 11, most of the formulations outperform the rule based controller with more than 90% decrease in discomfort. The main reason is the ability of MPC to predict the step change in zone heating setpoint temperature and compensate for the delayed response of the floor heating system. However, the MPC7 solution using Bonmin with the Hybrid method led to discomfort similar to the baseline controller. The authors managed to make MPC7 work properly, not shown in the chart. However, it required significant manual tuning effort for the weight on comfort constraints and solver internal options (MINLP approximation relaxation, integer tolerance) to guide the solution. This highlights that Bonmin-Hyb is probably not robust enough for this type of problem.
Looking at the computational time ratio K timr in Figure  11 and total computational time K ttot in Figure 12, QP and NLPs (MPC1 -MPC4) required less computing time than MINLPs (MPC5-MPC10) as expected. However, all MPC formulations have a K timr value much lower than one, meaning that MINLPs could be used for a real time application from the computation perspective. Looking at a Figure 11. BOPTEST KPIS for the month of January as thermal discomfort K dis (top), computational time ratio K timr in logarithmic scale (middle), and energy cost K cost (bottom). The results are shown for all MPC combinations explained in Table 5.  Table 6. The results are shown for all MPC combinations explained in Table 5 This big variation of computational times shows that MINLP MPCs are sensitive to a formulation, initialization, and solver selection. In fact, the differences between MPC5 and MPC6 are only the initialization approach and the inclusion of j B , and the difference between MPC9 and MPC10 is only the solver. However, MPC5 and MPC10 result in much longer computational time (a factor of 10) caused by a lot of errors or timeouts for around 30% of the iterations (see the last subfigure in Figure 12), meaning that the MINLP solvers are struggling to converge. Furthermore, SCIP was not able to find a solution without the j B binary objective. It is very interesting that all MINLP solvers benefit from the introduction of j B in terms of solver reliability and time, because the objective seems to be redundant for integer formulations and solvers. This might be explained by the fact that by introducing j B , each NLP approximation of the MINLP is itself an approximation of a MINLP problem. Furthermore, the addition of the binary constraint j B reduces the available control space, making the MINLP objective function more convex and bringing it closer to a Mixed Integer Convex Programming (MICP) problem.
The comparisons of energy cost K cost and thermal energy K en are shown in Figure 11 and in Figure 12, respectively. For K en , all formulations and the baseline are within 6 % of each other, with MPC2, MPC4, MPC6, MPC7, and MPC9 being marginally better than the baseline, and MPC1, MPC3, MPC8 and MPC10 being marginally worse. The MPCs perform in a similar fashion with respect to the RBC because the case study considered is a newly built low transmittance building with low consumption. Therefore, not much energy can be wasted by a properly tuned RBC. Furthermore, the RBC is under heating the building as it can be seen by looking at the discomfort Kdis. Instead, the MPCs do a better than the RBC in shifting the thermal energy supplied to anticipate the demand of the setpoint to drastically reduce discomfort, while keeping the thermal input lower thanks to better management of the partial loads When looking at cost K cost however, MPC1 (QP) is the worst performer with 5% increase with respect to the baseline RBC 2 , while MPC2 (NLP), MPC4 (NLP), MPC6 (MINLP) and MPC9 (MINLP) show the best performance with a 13% decrease in cost with respect to the baseline, which is aligned with other literature results. This is due to the fact that MPC1 adopts a simpler COP model,that is only a function of external temperature since the supply temperature is not a control variable, while the other formulations use a more precise COP definition including supply temperature allowing a more efficient use of the heat pump by keeping on average a lower supply temperature. It is not clear from this analysis alone why the MPC2, MPC4, MPC6 and MPC9 seem to outperform the other nonlinear formulations for these KPIs. For that, we consider a more detailed analysis of time series data in Section 3.2. Although there are cost savings and comfort improvements for those MPCs (see Figure 11), they may not be applicable in practice due to short cycling which increases the wear of components. It is also interesting to see the   Table 6. In this case the total power supplied by the heat pump was used to calculate K conlen .

Typical day analysis
influence of relaxation schemes for integer constraints: The only difference of MPC3 (NLP) compared with MPC2 (LNP) and MPC4 (NLP) that are in the smooth group (middle figures) is in rounding the MPC decision to enforce it to either 0 or 1. This rounding clearly introduces prediction errors and hence pushes MPC3 to solve an unexpected problem at the next time step. The aggressive result of MPC3 indicates that nonlinear MPC could be sensitive to prediction errors at least for this case study. Lastly, the frequent on-off of the system causes a higher discrepancy between the prediction of the low order model versus the emulator model. The main reason is that the Buildings envelope model is made up of hundreds of states while the reduced order model uses only a handful. The consequence is that the high frequency components for the aggressive case impact the temperature nodes in different ways between the emulator and reduced order models, increasing the error of the prediction. From the comprehensive analysis using KPIs of discomfort K dis , cost K cost and on-off switching K conlen , MPC2 (NLP), MPC4 (NLP), MPC6 (MINLP) and MPC9 (MINLP) that are in the smooth group are superior than other MPCs and the baseline, and MPC9 performs the best. However, as pointed out in Section 3.1, MINLP MPC formulations are sensitive to formulations, initialization approaches, solver selections and solver parameters, and the incremental savings of MPC9 are insignificant. Overall, we conclude that increasing the complexity of MINLP-MPC approaches does not bring substantial benefits at least for this case study but only adds computational burden, makes the MPC less robust, and requires more manual tuning efforts. The open-source and commercially available, general purpose MINLP solvers are not yet sufficiently reliable for real time MPC applications for the HVAC system. The newly introduced KPI K conlen described in Table 6 can compactly indicate the level of short cycling. Indeed, from Figure 14, which compares K conlen , the same conclusion from the analysis of Figure 13 can be drawn: the smooth group (middle figures: MPC2, MPC4, MPC6 and MPC9) has half the control arc length of the aggressive group (right figures: MPC3, MPC5, MPC8 and MPC10). However, the difference between baseline and MPC1 which use an on-off approach and MPC2,MPC4,MPC6 and MPC9 which modulate the load more is not highlighted. The reason is that the control variable considered is the total heat pump powerQ, which partially masks the behaviour of separate control variables, water supply setpoint T in,set and valves on-off u i . The authors are aware of this simplification and will explore a more comprehensive approach in future work.

Conclusions
This manuscript compared ten combinations of MPC formulations and solvers ranging between a QP formulation (MPC1), NLP formulations (MPC2-MPC4), and MINLPs (MPC5-MPC10) with a high focus on the assessment of MINLP MPCs for a residential radiant floor heating system. The performances were evaluated with a detailed emulator model where all MPCs exhibit unmodeled dynamics. The conclusions are summarized as below.
• MPC1 and MPC2 are both linear approximations, for the constraints, of the original MINLP problem. However, there is a large performance gap between the two. MPC1 is the most intuitive way to linearize the heat flow rate by combining valve and temperature into thermal powerQ. However, in this way, the MPC is not able to accurately estimate the heat pump COP as a function of supply temperature, leading to an on-off operation of the system. MPC2, instead, used a less intuitive approach, keeping the valves always open and modulating only the supply temperature in order to maximize COP. Here, the more detailed COP formulation led to an overall better solution. • MINLP MPC performance is sensitive to formulations, MINLP solvers, solver options and tolerances, and the corresponding tuning process is not trivial. However, a MINLP formulation is natural and more straightforward to implement, since the optimization variables could be consistent with the physical controllable variables and modelling approaches for simulations (e.g. the detailed COP model or use of physical valve positions as optimization variables). • To make MINLP MPC properly work, it often required solving another, approximated MPC problem. In our case, an NLP MPC was solved for initializing an MINLP MPC. • For this case study, no substantial benefits of any MINLP MPC formulation were found over an NLP MPC. Meanwhile, MINLP formulations dramatically increased computational time and were less stable, with time outs or errors occurring for over 30% of MPC control steps. In terms of generalization of the results, the authors think that for residential newly built buildings, the insulation between apartments is high enough, so that they can be considered almost adiabatic. Therefore, controlling the whole building would be almost equivalent to having a specific MPC for each apartment, making the problem very similar to the one dealt with in the manuscript. However, in the case of older buildings where the heat transfer between apartments is important, the coupling between the different MPCs must be done more carefully if using a distributed approach. Instead, a centralized approach would lead to a similar conclusion since more complex problem formulations are needed and MINLP solvers suffer exponentially more, in terms of computational time and robustness, the larger the optimal control problem becomes. • This paper introduced and compared the state-ofthe-art, open source MINLP solvers (SCIP, Bonmin) and commercial solver (Baron). They could be viable options for building HVAC optimal control problems with appropriate tuning, but they were not the best solution for this case study. • From the overall analysis, it is shown that the initial effort of finding a suitable linear constraint formulation or reduce the unwanted control space (jB), (in modelling the radiant floor heat transfer) leads to identical energy and discomfort performance compared with a very well tuned MINLP MPC, while being faster and more robust.In fact, a lot of effort went into the optimal control community, other than the building control field, to try and locally approximate the MINLP as MILP or MICP with various approximation methods. Future work will investigate the performance of different approximation methods applied to building HVAC control problems. Therefore, the authors conclude that using a linear formulation for the constraints with a nonlinear objective function and proper initialization method, such as slightly randomized free floating solutions, gives a good balance between the accuracy of the prediction, computational requirements, and robustness. Furthermore, the outcome of this manuscript shows that the introduction of additional sensors, control inputs and network infrastructure to allow for a more complex predictive control system may not necessarily bring significant improvements with respect to a simpler, more robust formulation approach. Therefore, it is important for implementations to consider the reliability of the whole infrastructure.