Configurations of model predictive control to exploit energy flexibility in building thermal loads

Abstract— A model predictive control (MPC) framework is developed in the present study, with the final objective to improve the energy flexibility of building thermal loads through demand-side management. Three different configurations are tested and tuned, with the following objective functions: minimizing the delivered energy to the building, the electrical energy used by the HVAC system (heat pump) or the cost of this electricity use. To validate these MPC configurations, a MatlabTrnsys co-simulator is also created, in order to run the MPC on a virtual plant composed of a detailed building model. The MPC strategy manages to run effectively on the chosen study case (a residential building with heat pump in Spain), and the differences between configurations are discussed. I. INTRODUCTION Demand-side management (DSM) is becoming a promising solution to solve the balancing issues created by the increasing penetration of non-dispatchable sources of energies (such as wind and solar) in the national grids. In particular, the thermal loads of buildings (heating, cooling and production of domestic hot water or DHW) are partly shiftable in time and thus represent good candidates for DSM schemes [1]. To unlock this flexibility potential of buildings, smart control strategies are however necessary. Indirect control strategies are generally preferred over direct control of the heating/cooling devices, for the satisfaction of the users. Among indirect controls, rule-based strategies reacting on penalty signals such as the time-variant electricity price have shown a satisfactory performance in terms of load-shifting while having a low complexity level. Model predictive control (MPC) outperforms such simpler strategies, but presents more difficulties in its implementation, in great part due to the need for accurate linear models [2]. MPC has been extensively studied for building climate control [3], [4], however not necessarily for unlocking flexibility. The majority of the reviewed literature on the topic resorted to economic MPC, e.g. a control strategy intending to minimize the energy costs. Few articles have considered other objective functions, for instance rewarding the completion of a flexibility request coming from an aggregator at an upper level [5]. No direct comparison between different MPC configurations for energy flexibility were carried out, to the knowledge of the authors, and the studies considering both heating and cooling on the same study case are rare [6]. The present work intends to develop an integrated MPC approach to harness the flexibility potential of the thermal loads provided to the building by a heat pump system. Different configurations of the MPC are presented, in T. Péan and J. Salom are with the Catalonia Institute for Energy Research (IREC), 08930 Sant Adriá de Besós (Barcelona), Spain (corresponding author phone: +34 933562615, email: [email protected], [email protected]). particular with different objective functions aiming to reduce the delivered thermal energy, the used electrical energy, or the associated economical costs. These configurations are tested through co-simulation with a dynamic building simulation tool, and evaluated with the same Key Performance Indicators (KPIs). II. MODELLING OF THE SUBSYSTEMS A. Presentation of the study case The MPC configuration is developed for a test building, typical of the Spanish and Mediterranean climate areas. It corresponds to a flat for a family of four members, within a multi-family building block. A refurbished version of the building is considered, with an additional 12 cm layer of insulation. A reversible air-to-water heat pump provides the heating or cooling needed to condition the space or the water tank used to store domestic hot water (DHW). The heated or cooled water is circulated through Fan Coil Units (FCU) installed in each room to ensure thermal comfort of the occupants. A detailed model was realized in TRNSYS and validated through metered data [7]. This model was utilized to create training and validation datasets, to obtain the reduced order model described in the following sections. To this end, the building model was virtually excited with a heat input following a Pseudo-Random Binary Signal (PRBS) [8]. B. Building envelope Figure 1. Scheme of the RC simplified building model. The building envelope is represented through a simple RC (resistance-capacity) grey-box model, analog to electrical networks [8]. The scheme is presented in Figure 1. In the present case, two states are defined: - 𝑇"#$ ∈ ℝ'[℃] is the indoor operative temperature in the zone (average of the temperatures in the different rooms of the apartment). 𝑇"#$ represents the main output of the model and should be contained within T. Péan and R. Costa Castelló are with the Department of Automatic Control (ESAII), Universitat Politècnica de Catalunya, 08028 Barcelona, Spain (e-mail: [email protected]). Configurations of model predictive control to exploit energy flexibility in building thermal loads Thibault Péan, Jaume Salom and Ramon Costa-Castelló, Member, IEEE the comfort boundaries defined by the occupants 𝑇"#$;'𝑇"#$ which differ in heating and cooling mode and according to the time of the day. - 𝑇,∈ ℝ'[℃] is an intermediate temperature state within the wall envelope (non-observed state). Two states enable to represent both the dynamics of the internal air node (fast response) and of the building fabric itself (walls, furniture, indoor thermal mass with higher inertia). More complex models can be developed, but the second-order model has shown sufficient performance for the type of application described in this work [9]. The evolution of the temperature states within the building are governed by the following set of equations: 𝐶"#$ ∙ 𝑇/#$ =1 𝑅"#$ 𝑇,− 𝑇"#$ +𝑔𝐴 ∙ 𝐼8+ 𝑄:+ 𝑄;<< (1) 𝐶,∙ 𝑇,=1 𝑅"#$ 𝑇"#$ − 𝑇,+1 𝑅, 𝑇=− 𝑇, (2) The exogenous input variables (or disturbances) to the building model are: - 𝑄:∈ ℝ''[𝑘𝑊]: the thermal power delivered to the space. In space heating mode, 𝑄:> 0, and in space cooling mode, 𝑄:< 0. 𝑄: is a controllable input, which value will be chosen by the MPC controller. - 𝑇=∈ ℝ [°C]: the ambient outdoor temperature; - 𝐼8∈ ℝC'[𝑘𝑊/𝑚F]: the ground horizontal solar irradiation. - 𝑄;<< ∈ ℝC'[𝑘𝑊]: the uncontrollable heat inputs within the building envelope, mostly due to the occupants and equipment. These heat inputs are represented through deterministic profiles. The parameters of the model are the following, their values obtained from the model identification process are presented in Table 1 : - 𝐶"#$, 𝐶,∈ ℝC F'[kWh/K]: the thermal capacities associated with the two states, representing the thermal mass of the building; - 𝑅"#$, 𝑅,∈ ℝC F'[𝐾/𝑘𝑊]: are respectively the thermal resistance between the two states, and the wall resistance between the wall node and the outside temperature; - 𝑔𝐴 ∈ ℝC'[𝑚F]: the coefficient (aperture area) defining the proportion of solar irradiation 𝐼8 entering the building through the windows (solar gains). C. Thermal Energy Storage (TES) To store the DHW, a water tank of 260 liters is integrated into the heat pump internal unit. Let us define the following state variable: - 𝑇MN: ∈ ℝ''[℃] is the water temperature in the TES and is the other output of the model. By law and to avoid spread of legionella disease [10], the TES water must be kept above a certain temperature 𝑇MN: > 𝑇MN:; however there is no upper boundary for this state. The evolution of the state temperature 𝑇MN: is governed by the following equation: 𝐶MN: ∙ 𝑇MN: = 𝑄MN: +1 𝑅MN: 𝑇"#$ − 𝑇MN: − 𝑄O8P (3) Where: - 𝑄MN: ∈ ℝC''[𝑘𝑊]: the thermal heating power delivered by the heat pump to the TES tank. 𝑄MN: is a controllable input (continuous variable). Similarly to 𝑄:, 𝑄MN: is constrained within the operation boundaries 𝑄8Q;'𝑄8Q of the heat pump system. - 𝑄O8P ∈ ℝC''[𝑘𝑊]: the DHW tapping from the occupants (deterministic standard tapping profile). - 𝑅MN: ∈ ℝC''[𝐾/𝑘𝑊]: the thermal resistance between the water in the tank and its surrounding environment (due to the tank insulation). - 𝐶MN: ∈ ℝC''[𝑘𝑊ℎ/𝐾]: the thermal capacity of the water tank, defined by 𝐶MN: = 𝑀,∙ 𝑐U,,, with 𝑀, the mass of the water contained in the tank 𝑐U,, the specific heat of water. TABLE 1. PARAMETERS AND THEIR VALUES. Model parameters with their values and units 𝐶"#$ = 0.26'𝑘𝑊ℎ/𝐾 𝑅"#$ = 0.42'𝐾/𝑘𝑊 𝐶,=19.1'𝑘𝑊ℎ/𝐾 𝑅,= 8.86'𝐾/𝑘𝑊 𝐶MN: = 0.290'𝑘𝑊ℎ/𝐾 𝑅MN: =601'𝐾/𝑘𝑊 𝑔𝐴 = 1.92'𝑚F 𝛼]= 0.0259 𝛾 = 0.222'𝑘𝑊/𝐾 𝛼`= 0.00838'℃b` 𝛽 = 0.557'𝑘𝑊/𝐾 𝛼F= −0.00704℃b` D. Emitters 𝑄: and 𝑄MN: (i.e. the thermal powers transmitted respectively to the indoor space and to the TES tank) have been chosen as the controllable variables to simplify the optimal control problem. However, a real heat pump system cannot be fed directly with a command of delivered thermal heating (or cooling) thermal power. In general, it is the supply temperature which is sent as an input to the system. In order to transform 𝑄: and 𝑄MN: into supply temperature signals, the following linear equations are used: 𝑄:= 𝛾 𝑇efU,: − 𝑇"#$ (4) 𝑄MN: = 𝛽 𝑇efU,MN: − 𝑇MN: (5) With: - 𝑇efU,: ∈ ℝ'[℃]: the water supply temperature for space heating/cooling (i.e. the temperature of the water flow leaving the heat pump and entering the fan coil units); - 𝑇efU,MN: ∈ ℝ'[℃]: the water supply temperature for TES charging (i.e. temperature of the water leaving the heat pump and entering the TES heat exchanger). Given that the time step used 𝑡e=12'𝑚𝑖𝑛 is higher than the dynamics of the heat pump, it is considered that the system can reach this set-point value within the allowed time step. The dynamics of the heat pump are thus not considered, the aim of the MPC consists mainly in determining the operation plan of the system over the time horizon. Equation (4) represents the final emitter to the building space. In the present case, FCU are used so that both heating and cooling operation are possible (radiators would have a slightly different behavior for instance). A linear regression from data generated with a detailed FCU model in TRNSYS, sized for the present building study case, gave the value of 𝛾 presented in TABLE 1. 𝛾 depends on the water and air flow rates on both sides of the FCU, and the efficiency of the heat exchange between the two heat carriers. Equation (5) represents the heat exchange from the heat pump circuit to the TES tank. 𝛽 = 𝑚,∙ 𝑐U,,, with 𝑚, the water mass flow rate on the heat pump side. E. Heat pump performance According to [11], the coefficient of performance (COP) of a heat pump can be approximated with a quadratic function. Since only the quantity 1𝐶𝑂𝑃 is used in the objective function (to obtain the electrical power used by the heat pump system), a linear model of this quantity is derived: Figure 2. Comparison between the quadratic model of the COP (full lines and points), and the linear model of 1/COP (dashed lines). The COP is defined as the ratio between the thermal power 𝑄: or 𝑄MN: and the electricity consumption of the heat pump 𝑃8Q. Its value depends on both the supply temperature 𝑇efU and the outside temperature 𝑇=. A comparison between the quadratic model of the COP and the linear model of 1/COP is presented in Figure 2. A satisfactory fit is observed, except for few high values of 𝑇= and low values of 𝑇efU. III. OPTIMAL CONTROL PROBLEM A. Overall state-space model Equations (1) to (5) can be summarized into a generic statespace model of the form: 𝒙 = 𝑨𝒙 + 𝑩𝑪𝒖𝑪+ 𝑩𝑿𝒖𝑿 𝒚 = 𝑪𝒙''''''''''''''''''''''''''''''''' (7) It should be noted that the traditional formulation 𝒙 = 𝑨𝒙 +𝑩𝒖 has been modified, so that the vector of the controllable input variables 𝒖𝑪 is separated from the vector of the exogenous input variables 𝒖𝑿. These vectors are as follows: 𝒙 = 𝑇"#$ 𝑇, 𝑇MN: , 𝒖𝑪=𝑄: 𝑄MN: , 𝒖𝑿= 𝑇= 𝐼8 𝑄;<< 𝑄O8P and 𝒚 = 𝑇"#$ 𝑇MN: And the model matrices: 𝑨 = −1 𝑅"#$𝐶"#$ 1 𝑅"#$𝐶"#$ 0 1 𝑅"#$𝐶, −1 𝑅,𝐶, −1 𝑅"#$𝐶, 0 1 𝑅MN:𝐶MN: 0 − 1 𝑅MN:𝐶MN: 𝑩𝑪= ` stuv 0 0 0 0` swxy , 𝑩𝑿= 0z{ stuv ` stuv 0 ` |}s}0 0 0 0 0 0 − ` swxy , 𝑪 = 100 001 The model is discretized with a discretization time step 𝑡e= 12'𝑚𝑖𝑛. B. Specifications Regarding the input constraints, the thermal power delivered by the heat pump (𝑄: or 𝑄MN:) is constrained within the operation boundaries 𝑄8Q;'𝑄8Q = [2; 4'𝑘𝑊] when the heat pump is turned on [12]. This modulation is made possible by the presence of an inverter controlling the compressor of the heat pump, however for operational reasons, the system cannot work at very low part-load ratios. TABLE 2. COMFORT BOUNDARIES ACCORDING TO OCCUPANCY AND SEASON. 𝑻𝒊𝒏𝒕;'𝑻𝒊𝒏𝒕 Winter Summer Occupied periods 20℃;24℃ 22℃;25℃ Unoccupied periods 18℃;24℃ 22℃;30℃ As for the output constraints, the indoor temperature 𝑇"#$ must stay within the comfort range 𝑇"#$;'𝑇"#$ . The values of the comfort boundaries are presented in Table 2 and are derived from the European norm [13]; a setback of the temperature is generally introduced when the occupants are not present. The lower boundary for the TES water temperature is fixed to 𝑇MN: =50℃ [10]. C. Optimal control problem formulation In the aforementioned model, the space heating/cooling operation (indexed 𝑆) and the TES charging operation (indexed 𝑇𝐸𝑆) have been separated. However, only one heat pump system covers these two different needs, and only one mode can be operated at a time (in practice, the output flow of the heat pump is directed through different circuits by means 1 𝐶𝑂𝑃 =𝑃8Q 𝑄:+ 𝑄MN: = 𝛼]+ 𝛼`𝑇=+ 𝛼F𝑇efU (6) of a divergent valve). To account for these constraints, two binary variables 𝛿: and 𝛿MN: ' ∈ {0,1} are introduced, describing the operation of the heat pump in both modes. The optimal control problem (OCP) is described as follows: Problem 1: min f‹,Œ 𝐽 = 𝛼;Ž• ∙ 𝐽∆f + (1 − 𝛼;Ž•) ∙ 𝐽;Ž• s.t. ∀'𝑘' ∈ 1, 𝑁 : Model: 𝒙(𝑘 + 1) = 𝑨 ∙ 𝒙(𝑘) + 𝑩𝑪∙ 𝒖𝑪(𝑘) + 𝑩𝑿∙ 𝒖𝑿(𝑘) 𝒚 𝑘 + 1 = 𝑪 ∙ 𝒙 𝑘 ''''''''''''''''''''''''''''''''''''''''''''''''''' Input constraints: 𝛿:(𝑘) ∙ 𝑄8Q ≤ 𝑄:(𝑘) ≤ ' 𝛿:(𝑘) ∙ 𝑄8Q 𝛿MN:(𝑘) ∙ 𝑄8Q ≤ 𝑄MN:(𝑘) ≤ ' 𝛿MN:(𝑘) ∙ 𝑄8Q 𝛿:𝑘 + 𝛿MN: 𝑘 ≤ 1 Output constraints: 𝑇"#$(𝑘) ≤ 𝑇"#$(𝑘) ≤ ' 𝑇"#$(𝑘) 𝑇MN: ≤ 𝑇MN:(𝑘) D. Objective function The general objective function has the form 𝐽 = 𝛼;Ž• ∙ 𝐽∆f + (1 − 𝛼;Ž•) ∙ 𝐽;Ž•, as shown in the OCP formulation. 𝐽∆f is a smoothing term that is always included, in order to penalize the changes in control actions (and thus avoid power spikes), and to facilitate the OCP solving: 𝐽;Ž• represents the real objective that the MPC controller intends to minimize. In the present study, three different objectives have been studied: minimization of the thermal energy delivered by the heat pump with 𝐽;Ž• = 𝐽•#, minimization of the electricity used by the heat pump with 𝐽;Ž• = 𝐽•–, and minimization of the cost of the electricity use with 𝐽;Ž• = 𝐽<;e$. These terms are explicitly described in (9) to (11). 𝐽•# considers directly the minimization of the control actions 𝑄: and 𝑄MN:, in a similar fashion than for instance [5]: 𝐽•# = 𝒖𝒄(𝑘) ˜ ™š` = 𝑄:(𝑘) + 𝑄MN:(𝑘) ˜ ™š` (9) However, this formulation does not take into account the performance of the heat pump, which changes significantly over time (see (6)). For this reason, the term 𝐽•– was introduced: it calculates the electricity needed by the heat pump to deliver its required thermal output: 𝐽•– =𝒖𝒄𝑘 𝐶𝑂𝑃(𝑘) ˜ ™š` =𝑄:(𝑘) + 𝑄MN:(𝑘) 𝐶𝑂𝑃(𝑘) ˜ ™š` (10) Further than the electricity used by the heat pump, the cost of this energy use is taken into account with 𝐽<;e$ (11). If a constant tariff for electricity is considered, 𝐽•– and 𝐽<;e$ are equivalent. However, the large rollout of smart meters in Europe enables the application of dynamic hourly tariffs for electricity. In fact, such tariffs already exist in Spain (denominated as “voluntary price for small consumers”) [14]. The variations of the electricity price create opportunities for economical savings if the heat pump is operated during low price periods instead of high price periods (load-shifting). 𝐽<;e$ =𝒖𝒄𝑘 𝐶𝑂𝑃(𝑘) ˜ ™š` ∙𝑃𝑟𝑖𝑐𝑒•–(𝑘) (11) The normalized electricity price 𝑃𝑟𝑖𝑐𝑒•–(𝑘) is here considered as the penalty signal. In fact, other penalty signals such as the CO2 emissions of the electricity production at national scale could also be used with the same structure, hence the generalization potential of the present work. E. Pareto fronts Figure 3. Pareto fronts for the different objectives. The objective term 𝐽;Ž• is balanced with the smoothing term 𝐽∆f by means of the weighting factor 𝛼;Ž• (which obviously depends on the chosen objective). To determine appropriate values of the weighting factors, the Pareto fronts are represented in the cases of the three different objectives aforementioned (see Figure 3). In all three cases, values between 0.2 and 0.5 seem a satisfactory compromise between the two objectives. Higher values of 𝛼;Ž• (i.e. giving more weight to the smoothing term) lead to a faster solution, but where the heat pump is turned on all the time to avoid changes in the control action, which is not desirable. The cases with lower values of 𝛼;Ž• require more computation time to encounter a solution and lead to more frequent switching of the heat pump (sometimes during as low as one time step), which is not desirable either. For the thermal energy minimization case, the value 𝛼•# = 0.01 is chosen since computation time is not an issue in this case, and 𝛼•– = 𝛼<;e$ = 0.25 for the other two cases. F. Solving of the OCP The MILP OCP problem is formulated in MATLAB, with the help of the Yalmip tool [15], and solved with the Gurobi solver [16]. The discretization time step is chosen as 𝑡e= 12'𝑚𝑖𝑛, and the prediction horizon 𝑁 = 120'time steps (24 hours), which enable to cover the daily patterns observed both in occupancy and weather. The control horizon is set to 1 hour (5 time steps of 12 minutes), meaning that the controller receives feedback and calculates a new optimal trajectory at this frequency. A shorter control horizon would ensure a 0.01 0.1 0.3 0.5 0.7 0.9 0.99 0.01 0.1 0.3 0.5 0.7 0.9 0.99 0.01 0.1 0.3 0.5 0.7 0.9 0.99 0 30 60 90 120 150 180 210 0 50 100 150 200 250 300 350 010 20 30 40 50 Objective Jel (electricity) Objective Jen (thermal energy) Objective Jcost (cost of electricity) Objective JΔu(smoothing) Thermal power Cost 𝐽∆f = 𝒖𝒄𝑘 − 𝒖𝒄𝑘 − 1 ˜ ™šF (8) faster reaction in case of unexpected event within the 1 hour open loop, but would also increase the computational effort. IV. RESULTS OF THE MPC APPLIED TO THE TEST CASE A. Comparison on a one-day prediction The MPC controller previously described was applied by simulations to the building study case described in section II.A. The three different objectives are tested separately and compared in Figure 4 and TABLE 3. Figure 4. Time series of the inputs and outputs for the 3 studied cases. On the top graph, the inputs QS (in blue shades) and QTES (in orange shades) are shown. In the middle graph, the input Tint is shown with its lower constraint, and on the bottom graph, the output TTES with its lower constraint. TABLE 3. RESULTS OF THE OCP CALCULATION OVER THE TIME HORIZON FOR A TEST CASE OF ONE WINTER DAY (JANUARY 2016), FOR THE THREE DIFFERENT STUDIED OBJECTIVES. Case Jen Jel Jcost Thermal energy Q [kWh] 40.32 41.65 43.38 [%] - +3.3% +7.6% Electrical energy P [kWh] 13.22 12.62 12.70 [%] - -4.5% -3.9% Cost [EUR] 1.59 1.53 1.53 [%] - -3.9% -3.6% Computation time [s] 9.10 290.92 388.42 [%] - +3096.9% +4168.4% All the cases perform satisfactorily and achieve their objective while complying with all the imposed constraints. The computation time is greatly increased by the introduction of the COP expression in the objective function, which makes it non-linear. Case 𝐽•# presents the lower thermal energy delivered, and case 𝐽•– presents the lower electricity use (which are their respective objectives). Unexpectedly, case 𝐽•– performs slightly better in terms of cost than 𝐽<;e$, probably due to the different balance between the smoothing term and the objective term. However these two cases present a very similar behavior. Given these observations, the 𝐽•– case was chosen for further evaluation with dynamic simulations. The load shifting can be observed in Figure 4: for instance, space heating operation is anticipated during the night in the 𝐽•– and 𝐽<;e$ cases compared to the 𝐽•# case, so as to benefit from a better COP or a lower price at this moment. B. Dynamic co-simulation with a building simulation tool To evaluate the functioning of the MPC in a more realistic setup, a dynamic co-simulation is carried out. The MPC controller in MATLAB is coupled with the building simulation tool TRNSYS which serves as the controlled plant in this case, in a similar manner than [17]. This setup enables to test the control strategy without having to implement it in a real building. For the co-simulation, the MPC controller in MATLAB is called every hour and determines the optimal operation for the next day (𝑁 = 120 time steps of 12 minutes). Only the first 5 control actions (supply temperature set-points and on-off signals) are then sent to TRNSYS, which simulates the dynamic behavior with a more accurate model and a higher time resolution (3 minutes time steps). One typical winter week is simulated, with the weather data of January 2015. To avoid infeasibility, additional conditions are set to the constraints. At every MPC computation, it can happen that the initial conditions of the outputs 𝒚 do not verify the imposed constraints. In such case, the constraints of the first 3 time steps are relaxed (i.e. set below the values of the outputs, so that the first time steps verify the constraints). Even with these additional settings, the hard constraints are met for 99.7% of the time for the 𝑇MN: constraint, and between 93.4% (𝐽•# case) and 99.1% of the time (𝐽•– case) for the 𝑇"#$ constraint, during the considered week. To deal with the unfeasibility issue, a more elaborated manner would consist in softening the constraints, i.e. to permit constraint violations but associated to a certain penalty, reflected in the objective function (this is the subject of further work). A sample of time series is presented in Figure 5 (two days). Overall, it is observed than in the 𝐽•– minimization case, the heat pump tends to be operated during longer periods, but at lower power (i.e. at lower supply temperature and hence higher COP). The TES charging presents more flexibility in its time of activation, since its constraint is constant, while the space heating has to operate most often at the time where the 𝑇"#$ constraint is stricter. Compared to the 𝐽•# case, the 𝐽•– case increased the delivered thermal energy by +3.5% during the considered week, decreased the electricity use by -0.7% and increased the cost by +0.7% (because the constraints were met less often in the 𝐽•# case). 0 5 Thermal power [kW] Qs (Jen) Qtes (Jen) Qs (Jel) Qtes (Jel) Qs (Jcost) Qtes (Jcost) 17 18 19 20 21 Temperature [°C] Tint (Jen) Tint (Jel) Tint (Jcost) Constraint 45 50 55 60 65 0 3 6 9 12 15 18 21 24 Temperature [°C] Time[h] Ttes (Jen) Ttes (Jel) Ttes (Jcost) Constraint Figure 5. Time series of the co-simulation for the Jen and Jel cases. V. CONCLUSION In the present work, an MPC framework was developed and tested with the aim to enhance the energy flexibility of thermal loads in buildings equipped with heat pumps. In the proposed setup and among the three tested configurations, the minimization of the thermal energy delivered to the building presented the lowest computation effort. However, it does not take into account the efficiency of the heat pump, therefore the minimization of the electricity use and the associated costs were also tested. These two cases lead to similar results, and cost savings of up to 4% compared to the first MPC configuration. The high computation time of the two latter configurations could present some obstacles to their largescale implementation. However, with regards to the slow dynamics of building thermal loads, the solving time of the OCP is considered acceptable. The testing setup constitutes another valuable output of the present research. The connection between an optimization tool (controller in MATLAB) and a virtual plant (building detailed model in TRNSYS) enables to test a variety of predictive control strategies. In further research, the variation of the COP with the load could be included in the heat pump model, since its performance varies at part-load conditions [18]. Furthermore, the MPC configurations will be tested with a real heat pump placed in a semi-virtual environment setup (i.e. the heat pump will be connected to thermal benches which emulate the loads calculated by a dynamic building simulation tool). ACKNOWLEDGMENT This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 675318 (INCITE). Part of this work stems from the activities carried out in the framework of the IEA-EBC Annex 67 (International Energy Agency – Energy in Buildings and Communities program) about Energy-Flexible Buildings. REFERENCES [1] S. Ø. Jensen, A. Marszal-Pomianowska, R. Lollini, W. Pasut, A. Knotzer, P. Engelmann, A. Stafford, and G. Reynders, “IEA EBC Annex 67 Energy Flexible Buildings,” Energy Build., vol. 155, pp. 25–34, Aug. 2017. [2] D. Fischer and H. Madani, “On heat pumps in smart grids: A review,” Renew. Sustain. Energy Rev., vol. 70, no. October 2016, pp. 342–357, 2017. [3] Y. Ma, A. Kelman, A. Daly, and F. Borrelli, “Predictive Control for Energy Efficient Buildings with Thermal Storage: Modeling, Stimulation, and Experiments,” IEEE Control Syst., vol. 32, no. 1, pp. 44–64, 2012. [4] F. Oldewurtel, A. Parisio, C. N. Jones, D. Gyalistras, M. Gwerder, V. Stauch, B. Lehmann, and M. Morari, “Use of model predictive control and weather forecasts for energy efficient building climate control,” Energy Build., vol. 45, pp. 15–27, 2012. [5] G. Bianchini, M. Casini, D. Pepe, A. Vicino, and G. G. Zanvettor, “An Integrated MPC Approach for Demand-Response Heating and Energy Storage Operation in Smart Buildings,” IEEE Conf. Decis. Control 2017, no. Cdc, pp. 3865–3870, 2017. [6] T. Q. Péan, J. Salom, and R. Costa-Castelló, “Review of control strategies for improving the energy flexibility provided by heat pump systems in buildings,” J. Process Control, no. Special Issue on Efficient Energy Management, Apr. 2018. [7] J. Ortiz, F. Guarino, J. Salom, C. Corchero, and M. Cellura, “Stochastic model for electrical loads in Mediterranean residential buildings: Validation and applications,” Energy Build., vol. 80, pp. 23–36, 2014. [8] P. Bacher and H. Madsen, “Identifying suitable models for the heat dynamics of buildings,” Energy Build., vol. 43, no. 7, pp. 1511–1522, 2011. [9] G. Reynders, J. Diriken, and D. Saelens, “Quality of grey-box models and identified parameters as function of the accuracy of input and observation signals,” Energy Build., vol. 82, 2014. [10] Ministerio de Industria, “Guia tecnica. Agua caliente sanitaria central.,” 2010. [11] C. Verhelst, D. Degrauwe, F. Logist, J. Van Impe, and L. Helsen, “Multi-objective optimal control of an air-to-water heat pump for residential heating,” Build. Simul., vol. 5, no. 3, pp. 281–291, 2012. [12] Hitachi Air Conditioning Products Europe, “Technical Catalogue - Yutaki Series,” 2016. [13] CEN, “EN 15251: Indoor environmental input parameters for design and assessment of energy performance of buildings addressing indoor quality, thermal environment, lighting and acoustic.” European Committee for Standardization, Brussels, Belgium, 2007. [14] Red Electrica de España, “ESIOS – Sistema de información del operador del sistema.,” 2018. [Online]. Available: https://www.esios.ree.es/en (accessed 07/03/2018). [Accessed: 04Apr-2017]. [15] J. Lofberg, “YALMIP : a toolbox for modeling and optimization in MATLAB,” 2004 IEEE Int. Conf. Robot. Autom. (IEEE Cat. No.04CH37508), pp. 284–289, 2004. [16] Gurobi Optimization, “Gurobi.” http://www.gurobi.com, 2018. [17] N. Alibabaei, A. S. Fung, and K. Raahemifar, “Development of Matlab-TRNSYS co-simulator for applying predictive strategy planning models on residential house HVAC system,” Energy Build., vol. 128, pp. 81–98, 2016. [18] Y. J. Kim, E. Fuentes, and L. K. Norford, “Experimental Study of Grid Frequency Regulation Ancillary Service of a Variable Speed Heat Pump,” IEEE Trans. Power Syst., vol. 31, no. 4, pp. 3090–3099, 2016. 0 1 2 3 4 5 Power [kW] Q_S (Jen) Q_TES (Jen) Q_S (Jel) Q_TES (Jel) 17 18 19 20 21 22 Temperature [°C] T_int (Jen) Lower bound T_int T_int (Jel) 45 50 55 60 65 70 24 30 36 42 48 54 60 66 72 Temperature [°C] Time [h] T_TES (Jen) Lower bound T_TES T_TES (Jel)