Model predictive control of distributed energy resources in residential buildings considering forecast uncertainties

Forecast uncertainties pose a considerable challenge to the success of model predictive control (MPC) in buildings. Numerous possibilities for considering forecast uncertainties in MPCs are available, but an in-depth comparison is lacking. This paper compares two main approaches to consider uncertainties: robust and stochastic MPC. They are benchmarked against a deterministic MPC and an MPC with perfect forecast. The MPCs utilize a holistic building model to reﬂect modern smart homes that include photovoltaic power generation and storage, thermally controlled loads, and smart appliances. Real-world data are used to identify the thermal building model. The performance of the various controllers is investigated under three levels of uncertainty for two building models with diﬀerent envelope performance. For the highly insulated building, the deterministic MPC achieves satisfactory thermal comfort when the forecast error is medium or low, but the thermal comfort is compromised for high forecast errors. In the poorly insulated building, thermal comfort is compromised at medium and high forecast errors. Compared to the deterministic MPC, the robust MPC increases the electricity cost by up to 4.5% and provides complete temperature constraint satisfaction while the stochastic MPC increases the electricity cost by less than 1% and fulﬁlls the thermal comfort requirements.


Introduction
Climate change makes it urgent for people to improve energy efficiency and increase power generation from renewable sources.In the European Union (EU), approximately 40% of the final energy is consumed by buildings [1], over two-thirds of which is used for heating, ventilation, and air conditioning (HVAC) [1].To increase the energy efficiency of buildings, the EU has established the Energy Performance of Buildings Directive, which requires all new buildings achieve zeroemission as of 2030 and the worst-performing 15% of the building stock be upgraded [2].In addition to energy efficiency, the EU has set a new target of 40% renewables in the energy mix and a benchmark of 49% of renewables in buildings, as specified in its revised Renewable Energy Directive [3].Because of the inherent intermittency of renewable resources, increasing penetration of renewable power generation in the distribution brings a big challenge to resilient and reliable operation of the electric power grid.On the other hand, distributed energy resources, such as solar photovoltaic (PV), energy storage, electric vehicle chargers, together with the slow thermal dynamics of buildings, create the opportunity of managing building electric loads caused by HVAC systems, lighting and plug loads, leading to the so-called demand flexibility in literature [4,5].
A promising method to exploit the considerable potential of demand flexibility is MPC.As an advanced control algorithm that uses future predictions to derive optimal decisions, MPC has been demonstrated superior to conventional building controls via simulations and experiments (e.g.[6][7][8][9]).Many publications (e.g.[7][8][9][10][11]) applied the MPC to the control of HVAC systems only.Some studies (e.g., [12][13][14][15][16]) extended the MPC to integrate additional resources such as PV modules and battery energy storage.There are relatively few MPC studies ( [17,18]) that have considered a comprehensive set of distributed energy resources.
Since MPC relies on a process model to predict the future behavior of the controlled system, the uncertainties associated with the process model and the forecast of disturbances can deteriorate the control performance if they are not appropriately accounted for in the predictive optimization.Different MPC formulations are available to consider uncertainties and address their impact on the MPC performance [19].The https://doi.org/10.1016/j.enbuild.2023.113753Received 26 May 2023; Received in revised form 12 October 2023; Accepted 13 November 2023 three most commonly used approaches are deterministic MPC (DMPC), robust MPC (RMPC), and stochastic MPC (SMPC) [20]: • DMPC: DMPC is an MPC without explicit consideration of uncertainties [7].However, it is still able to cope with a certain degree of uncertainties because of its rolling horizon approach.This approach solves an open-loop problem over a finite time horizon but only implementing the derived inputs for the first time step.
Subsequently, another open-loop problem is formulated and the procedure is repeated.This repeated solving of the open-loop problems essentially introduces feedback into the system, which is the mechanism that the DMPC uses to compensate for forecast errors without explicitly considering them.• RMPC: RMPC explicitly accounts for uncertainties, expressed as the error bounds of model accuracy and exogenous disturbances, and addresses the uncertainties by assuming the worst-case scenario in the optimization [21].The worst-case scenario is the maximum possible uncertainty in the disturbances.Therefore, RMPC is also frequently called a minimax controller since it minimizes the cost function based on the assumption of the maximum possible uncertainty in the disturbances.Although the RMPC approach guarantees the satisfaction of constraints if the realized disturbances are within the predefined bounds, it may be overly conservative [22] and degrade the control performance.• SMPC: SMPC exploits knowledge about the statistical distribution of the forecast errors to satisfy the constraints with a user-defined probability [7].SMPC leads to a less conservative control strategy than RMPC because it allows infrequent constraint violations in minimizing the cost function.Although SMPC requires more implementation effort than both DMPC and RMPC, it is regarded as a very suitable control strategy for buildings [23], to deal with the constraint of thermal comfort, which can be violated for a small fraction of time during the year.
A variety of studies (e.g.[7,17,20,[24][25][26][27]) have applied robust or stochastic model predictive controls for buildings.Typically, these studies investigate whether explicit consideration of uncertainty in MPC improves the control performance (e.g., thermal comfort and energy consumption, renewabe energy use and cost) compared to the DMPC or conventional control.A selection of representative studies are discussed as state of the art in the following.
For the RMPC, Yang et al. [24] and Maasoumy et al. [25] used the same research methodology to handle uncertainties in model predictive control for energy efficiency in buildings.They incorporated adaptive modeling technique and robust control into MPC, leading to the so-called adaptive robust MPC (ARMPC).The adaptive feature refers to the regular update of the parameters of the building thermal model based on the online measurements of building operation data.The robust feature refers to the consideration of the uncertainties of disturbances in the predictive model and the use of the worst case during optimization to ensure that constraints are satisfied despite uncertainties.Yang et al. [24] compared the performance of ARMPC with AMPC (adaptive model and basic MPC formulation), basic MPCs with a static building model, and a conventional thermostat-based control under five uncertainty levels of internal loads.Based on the simulation in April using an energy-saving-biased objective function, they found that basic MPCs and AMPC consistently achieved about 20% energy savings relative to the thermostat-based control while the ARMPC decreased energy saving from 20% to 15% as the uncertainty level was increased from 0 to 60%.However, the ARMPC achieved the thermal comfort compliance ratio more than 90% for all different uncertainty levels but the other cases mostly achieved the thermal comfort compliance ratio between 54% and 70%.Similarly, Maasoumy et al. [25] compared the performance of ARMPC with a DMPC and a rule-based controller (RBC) under different levels of uncertainty due to weather forecasts and internal heat gains.They found that the DMPC performed best when the modeled uncertainties were lower than 30%, the RMPC performed best when the uncertainties were in between 30% and 67%, and the RBC performed best when the uncertainties were higher than 67%.Kim [26] investigated the use of RMPC for controlling the passive building thermal mass and the mechanical thermal energy storage and compared its performance with a DMPC.They found that both RMPC and DMPC had their strengths, and no superior one could be determined.When the forecast error was small, the conservative nature of the RMPC led to inferior performance compared to the DMPC.On the other hand, the conservative nature of RMPC was beneficial when the forecast error was large.
For the SMPC, Oldewurtel et al. [7] applied it for building climate control and compared its performance against a rule-base control and a DMPC for a number of cases with different buildings, HVAC systems, and weather conditions.Their results showed that SMPC was superior to RBC and DMPC with respect to energy consumption and thermal comfort.They attributed the excellent performance of SMPC to its capability to directly account for the uncertainty of weather forecasts in making control decisions.Nagpal et al. [17] investigated the performance of an SMPC in comparison to a DMPC for the operation of smart sustainable buildings that included a comprehensive set of distributed energy resources, such as PV, battery storage, thermostatically controlled loads, deferrable smart appliances and an electric vehicle.Based on a one-day simulation for a numerical case study that considered the uncertainties of ambient temperature and solar irradiance, they found that the electricity cost increased from 1.27e for the DMPC to 1.41e for the SMPC but meanwhile the thermal comfort constraint violation decreased from 1.27 • C to 0.02 • C. Amadeh et al. [27] proposed an SMPC-based framework to quantify demand flexibility accounting for the uncertainties of weather forecasts (ambient temperature and solar irradiance) and the building model.Based on a residential building with radiant floor heating, they compared the SMPC-and DMPC-based approaches for demand flexibility quantification over long periods (8-12 hours) under different winter conditions.They found that the DMPC-based approach tends to overestimate the potential of demand flexibility and the thermal comfort may be significantly compromised.For example, for the case of downward flexibility over 12 hours, the flexibility potential estimated by the SMPC was 9.5%-12.5% lower than the results from the DMPC as the ambient temperature uncertainty was increased from one to three times of standard deviation under the average winter weather conditions.However, the thermal comfort violation was ∼90 • C h for the DMPC while it was no more than ∼30 • C h for the SMPC.
Our literature review reveals two main research gaps: First, the two prevalent MPC formulations to explicitly consider uncertainties, SMPC and RMPC, are typically benchmarked against a DMPC and not against each other in literature.Thus, an in-depth comparison between RMPC and SMPC is still missing [23].The second research gap lies in the scarce consideration of holistic models to reflect the load-shifting capabilities of modern and future smart homes.We aim to address the above research gaps by considering three MPC variants (i.e., DMPC, RMPC, and SMPC) for demand-side management in residential buildings that include a rich set of controllable loads and distributed energy resources such as the HVAC system, schedulable and interruptible appliances, a hot water heater, photovoltaic power generation, and a battery storage system.The RMPC and SMPC are benchmarked against a DMPC and a performance bound (PB), which is a conceptual case where the uncertainty in the disturbance forecast is entirely known in advance.Hence, the performance of PB is not affected by the uncertainties.Moreover, RMPC and SMPC are compared based on the performance metrics of energy cost and thermal comfort.
The remainder of this paper is organized as follows.The modeling of individual components is presented first in Section 2. Based on the component models, three different MPCs are formulated in Section 3. The quantification of uncertainties and the test cases used to compare the performance of MPCs are discussed sequentially in Sections 4 and 5.The results are presented and analyzed in Section 6.The paper ends with conclusions and future work in section 7.

Component modeling
The building considered in this work consists of several components (see Fig. 1): the rooftop PV system for power generation, a battery for electrical energy storage, a heat pump for space heating, a domestic hot water tank (DHW), and deferrable electrical appliances (i.e., dishwasher, washing machine, and clothes dryer).The MPC optimizes the operation of the components while being able to purchase power from the grid and sell generated PV power to the grid.We define hereinafter the set  = {1, 2, ..., } as the prediction horizon,  ∈  as the time step, and Δ as the sample time.

Building thermal model
The MPC requires a thermal building model to predict its behavior in order to determine optimal heating sequences.A third-order resistancecapacitance model is applied to capture the thermal dynamics of the building, based on [28,29].This single-zone gray-box model consists of three state variables: the space air temperature  air , and the internal and external temperatures of the building envelope  w,in and  w,out .The model regards the heating or cooling output of the heat pump Φ h as the input variable and the global solar radiation Φ global , the ambient air temperature  amb and the internal heat gains Φ internal as the three where  conv is the convective fraction of the heat pump output and the internal heat gains and the value of 0.8 is used,  sol is an identifiable parameter that captures the influence of the global solar radiation Φ global on the air temperature,  air and  w are the respective heat capacities of the air and building envelope,  is the thermal resistance,  amb,eq represents the equivalent ambient temperature (i.e., sol-air temperature) after accounting for the effect of solar radiation on the exterior surface of the building envelope.The sol-air temperature is calculated as [29]: where ℎ f = 0.5 is the short-wave absorption coefficient of the exterior surface, and ℎ A = 25 W∕(m 2 K) the exterior heat transfer coefficient [29].The influence of the solar radiation on the interior side of the building envelope is neglected to simplify the model.
The heating or cooling output Φ h of the heat pump is calculated as: where  HP is the power consumption and  is the coefficient of performance (COP) of the heat pump.The continuous differential equations (Eq.( 1)-( 3)) are discretized with the sample time Δ, which yields the following time-discrete model: Where , ,  and  respectively denote the state variables (i.e.,  air ,  w,in and  w,out ), the model output (i.e.,  air ), the controlled input (i.e., Φ h ), the disturbances (i.e. amb , Φ global and Φ internal ).
The heating and cooling system is constrained by a specified temperature range for thermal comfort (Eq.( 8)) and the maximum power consumption of the heat pump (Eq.( 9)): air,t,min −  t ≤  air,t ≤  air,t,max +  t (8)  HP,min ≤  HP,t ≤  HP,max (9) where  t is a nonnegative slack variable to prevent numerical infeasibility in case the constraint without relaxation cannot be satisfied.
The occupancy schedules from [30], including one for weekdays and another for weekends, are used to calculate the internal heat gains from occupants.The schedules consider three occupant statuses: active at home (06:00 to 08:00, 17:00 to 23:00), sleep at home (23:00 to 06:00), and away (08:00 to 17:00).According to the German standard DIN V 18599 [31], the total heat emitted from occupancy, lighting, and electrical devices for a detached house over the course of one day is approximately 45 Wh∕m 2 .In this work, the daily 45 Wh∕m 2 is distributed as following: 10.5 Wh∕m 2 for the status of sleep, 6.75 Wh∕m 2 for the status of being away from home to account for electrical devices that are still running and 27.75 Wh∕m 2 for the status of being active at home.

Photovoltaic power generation
The photovoltaic power generation  PV at time  is described according to [32,33]: where  module is the number of PV modules,  module is the area of a single PV module (m 2 ),  STC is the electrical efficiency of the PV module at standard testing conditions,  L&T denotes the combined efficiency of the inverter and the maximum power point tracking controller,  t represents the solar irradiation incident on the PV module at time  (W∕m 2 ) and it is calculated according to [34],  T refers to the temperature coefficient of PV power (1∕ • C),  amb,t is the ambient air temperature ( • C) at time , and  represents the nominal operating cell temperature ( • C).

Battery storage
The energy stored in the battery is modeled according to [12,17]: where  B,t is the energy stored in the battery at time ,  ch,t and  d,t represent, respectively, the charging and discharging power at the time  and  ch and  d are respectively the charging and discharging efficiency.
Furthermore, the state of charge of the battery is calculated according to: where  B,max is the maximum energy that can be stored in the battery.
The battery storage system is constrained by: where  ch,t and  d,t are two binary variables used to indicate the charging or discharging status of the battery.For example, if  ch,t = 0 (no charging) then Eq. ( 14) yields  ch t = 0. Similarly, if  ch,t = 1 (charging occurs) then Eq. ( 14) yields  ch,min ≤  ch,t ≤  ch,max .Eq. ( 16) prevents simultaneous charging and discharging of the battery.

Domestic hot water heater
The model of the domestic hot water heater is based on [35]: where  tank ,  in and  air represent the temperature of the tank water ( • C), inlet water ( • C), the ambient air ( • C), respectively,   t is the hot water demand (m 3 ∕s),  tank is the tank volume (m 3 ),  tank is the heat capacity of the tank water (J∕ • C).  el represents the power consumption of the water heater (W),  el is the efficiency of converting electrical power into a heat flux,  tank denotes the skin heat loss coefficient of the tank (W∕m 2• C), and  tank is the tank's surface area (m 2 ).The hot water heater system is constrained by: F. Langner, W. Wang, M. Frahm et al.
The room temperature  air is uncertain due to forecast inaccuracies which could cause the violation of the hard temperature constraints for  tank .Therefore, the temperature constraint is softened with the slack variable .The hot water demand profile is retrieved from [36] and shifted by one hour to meet the occupancy schedule better.

Deferrable electric appliances
Model predictive controllers can schedule appliances to low power price periods to reduce the electricity bill.The occupant provides the controller with a time window for appliance operation, and the controller determines the optimal times to run the appliance.For example, the residents leave for work in the morning and want the dishes cleaned upon their return.The MPC would schedule the dishwasher to the most suitable times under the constraint that the dishwasher must have finished when the residents return.Electrical appliances can be categorized into interruptible and non-interruptible appliances.The operation of interruptible appliances such as dishwashers can be stopped and continued at a later point in time.This can contribute to increasing the flexibility of the power demand in buildings.
This work considers three appliances: a dishwasher, a washing machine, and a clothes dryer.The clothes dryer and the dishwasher are selected because they offer the largest potential of demand response amongst all appliances [37].Though the washing machine usually has a lower potential of demand response [37], it is included in the present study.
Two major approaches are available to model electric appliances: The first approach is to divide the appliance operation into different phases (e.g., the phases of pre-washing, washing, rinsing, and drying for dishwashers) with the power consumption being constant during each phase.The appliance operation can be suspended between phases but cannot be interrupted in the middle of a phase [38].The second approach is to refrain from subdividing the appliance's run time into uninterruptible phases and to allow interruptions at any time [39,40].The first approach is adopted in our study.
The present work uses the power consumption profiles from [41], where the power consumption of the dishwasher and the washing machine is structured in phases.For the MPC to schedule the appliances in the right time window and to enable interruptibility, constraints need to be formulated.The mathematical description of the constraints is based on [38].Let  be the set of appliances where each appliance  has a set of power phases  i .Each appliance  ∈  can be scheduled to operate in a user-defined time window the earliest starting time and  i F the latest finishing time.The power consumption of the -th appliance during the power phase  ∈  i is denoted as  i j .
For the mathematical formulation of the problem, two binary decision variables are necessary:  i t,j and  i t,j . i t,j indicates the on/off status of the -th appliance in its -th phase at time  and  i t,j indicates whether the -th power phase of the -th appliance is finished at time .

Constraint for power consumption
The following constraint ensures that the power consumption of the appliance  during the power phase  equals the power consumption profile  i j , which is expressed as: where  i t,j is the actual power consumption of the appliance .

Constraints for operational time
The decision variable  i t,j is set to 1 after the power phase  of appliance  finishes.Let  0 denote the last time step when phase  is processed.Then, for all time steps after  0 , the following condition holds: Eq. ( 22) protects the user-defined appliance operation time window from being jeopardized: Each phase of an appliance needs a specific length of time to complete during the time set  .Therefore, the following constraint is defined: where   j is the total runtime needed by phase  of appliance .

Constraints for continuous phase operation
Once phase  of appliance  starts, it will operate continuously until it finishes.The continuous operation of each individual phase is satisfied with the following: In the above constraints, Eq. ( 24) indicates that the status of a power phase can not be both running and finished.Eq. ( 25) forces the appliance to continue its operation in phase  at time step  if the phase has not finished yet, and Eq. ( 26) specifies that  i j remains to be 1 after phase  is finished.

Constraint for sequential operation
The operation of different phases of an electric appliance must follow a particular order.For example, for dishwashers, rinsing must be performed after washing and before drying.Such sequential operation of different phases is satisfied with the following equation: After minor changes (e.g., the superscripts for two dependent appliances), Eq. ( 27) can also be used to specify the sequential operation between different appliances.For example, the clothes dryer can run only after the washing machine finished its last phase.

Constraints for the time delay between phases
Interruptible appliances allow a phase to be delayed for operation after its precedent finishes.However, there may be limits on the time delay between consecutive power phases.For example, the time interval between the power phases of a washing machine cannot be too long to avoid ruining the clothes.For convenience, an auxiliary binary variable  is defined first in Eq. ( 28) for the time delay between two consecutive phases.The auxiliary variable  is then used in Eq. ( 29) to establish the time delay constraint: Eq. (28) shows that  t,j = 1 if phase  − 1 is finished but phase  has not started yet and  t,j = 0 otherwise.In Eq. ( 29),  i max,j and  i min,j are respectively the maximum and the minimum allowed number of time steps between phase  − 1 and phase  of appliance .Similar equations can be established to enforce the time delay constraint between two different appliances (e.g., the washing machine and the clothes dryer) as needed.

MPC problem formulation
To ensure the first law of thermodynamics, the power consumption in the building must equal the sum of the provided power by the grid, battery, and PV system.
In Eq. ( 30),  buy,t and  sell,t are the bought and sold power from and to the grid.Since the future solar radiation and the ambient temperature are uncertain, the future generated power by the photovoltaic module is uncertain.For  > 1, the forecast value of the PV power generation is used as prediction which is calculated with Eq. ( 10) and the weather forecast.The power exchange with the grid is constrained by:  buy,t ⋅  buy,min ≤  buy,t ≤  buy,t ⋅  buy,max (31)  sell,t ⋅  sell,min ≤  buy,t ≤  sell,t ⋅  sell,max (32) 0 ≤  buy,t +  sell,t ≤ 1 (33) where Eq. ( 33) prevents simultaneous purchasing and selling of power from and to the grid.All of the presented constraints have to be considered in the optimization.The optimization is conducted to minimize the financial expenses which is expressed in the following objective function.
The objective function is to minimize the cost of energy consumption over the entire time horizon.The objective function is penalized if the soft constraints are violated.It is calculated as follows min

⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞ ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞ ⏟
Cost of power transactions subject to equations ( 6) - (33) where Ω buy,t is the tariff for purchasing power and Ω sell,t the tariff for selling power to the grid. is a large positive number used to penalize the soft constraint violations as indicated by the slack variables  t and  t .Since the optimization includes binary variables, the resulting optimization problem is a mixed-integer linear problem.

Deterministic MPC
The DMPC is the standard MPC formulation that does not model uncertainties and thereby not explicitly consider them in the optimization.Even though the DMPC copes with the occurring forecast errors implicitly through the rolling horizon approach, it cannot address the impact of inaccurate forecasts and will inevitably cause the problem of constraint violations.

Stochastic MPC
Stochastic model predictive control allows constraint violations with a certain probability.In a stochastic MPC with additive noise on the disturbance forecast and input, the controller model becomes t =  t (36) where  =  HP ⋅ Δ is the uncertainty of the heat pump output due to the uncertainty of COP Δ,  are the stochastic forecast errors associated with the disturbance forecasts  (i.e. amb , Φ global , Φ internal ).The noises  t ,  t ∼  (0,  2 t ) are assumed to be independent and identically distributed and to follow a normal distribution with zero mean and variances  2 t .For the purpose of allowing constraint violations with a certain probability, the constraints are formulated in a probabilistic way as so-called chance constraints.Two different versions of chance constraints can be distinguished [20]: joint and individual chance constraints.Joint constraints enforce the fulfillment of multiple constraints jointly with the probability of 1 − , whereas individual constraints enforce each constraint independently with the probability of 1 − , where the symbol  indicates the acceptable constraint violation level.Because a tractable representation of joint constraints does not necessarily exist [42], individual chance constraints are applied in this work.Thus, the thermal comfort constraint (Eq.( 8)) is rewritten as ℙ( air,t = ( t-1 + ( t-1 +  t-1 ) + ( t-1 +  t-1 )) ≥  air,t,min ) = 1 −  (37) ℙ( air,t = ( t-1 + ( t-1 +  t-1 ) + ( t-1 +  t-1 )) ≤  air,t,max ) = 1 −  (38) For the compact notation of the reformulated constraints, we define the following matrices: describes the uncertainty propagation throughout the prediction horizon  , since the forecast uncertainties accumulate over  .This means that the forecast uncertainty increases as the forecast extends further into the future.To address this, the SMPC increasingly tightens the chance constraints across  (see Eq. ( 41) and Eq. ( 42)).According to [42], the chance constraints Eq. ( 37) and Eq. ( 38) can be analytically reformulated to their equivalent deterministic constraint as follows so that the SMPC problem is tractable: t ≤  air,t,max − Φ −1 (1 − ) where Φ −1 (⋅) is the inverse cumulative probability function of the standard normal distribution,  is the covariance matrix of the disturbances  and , and  t represents the -th row of .The slack variable  t is added after the reformulation to avoid infeasibility, analogously to Eq. ( 8).
The probability distribution of the disturbances leads to probabilistic states, the expected values of which are used to calculate the expected value of the objective function:

Robust MPC
The core idea of RMPC is to assume the worst-case scenario in the optimization [43].RMPC uses the same controller model as SMPC (Eq.(35) and Eq. ( 36)) but the uncertain variables ,  are bounded, and the optimization is conducted considering the maximum possible disturbances.Therefore, RMPC prevents all constraint violations if uncertain variables are realized within the provided bounds.This hedging against the worst-case scenario can lead to overly conservative results if the worst-case scenario is unlikely to enter.In this work, the bounds of the uncertain variables ,  are selected to be three times their respective standard deviation: The optimization in RMPC is conducted while assuming the maximum disturbances , , leading to the following objective function:

Uncertainty quantification
The major sources of uncertainty come from weather forecasts, occupant behavior, and the mismatch between the thermal model and the real building.In the present work, we focus on the uncertainties of weather (including ambient air temperature  amb and solar radiation Φ global ) and occupant behavior (i.e., the internal heat gains Φ internal ).
For an air-source heat pump, the uncertainty of ambient air temperature also results in an uncertainty of heat pump COP .
The considered uncertainties affect both indoor air temperature and photovoltaic power generation.A deviation of the observed ambient air temperature, solar radiation, or internal heat gains from their forecasts can result in mismatch between the predicted indoor air temperature and PV power generation and their actually observed values, leading to inferior MPC performance for power management and compromised thermal comfort.Therefore, accurate uncertainty quantification is an important part of RMPC and SMPC.
The weather uncertainties are modeled by comparing the forecasted data obtained from ECMWF [44] for Karlsruhe, Germany to the measured data from the German Weather Service [45] for the weather station in Rheinstetten.This weather station is the closest one (approximately 8 kilometers) to Karlsruhe, where the experimental building is located.Both data sources have an hourly time step.Based on the 12-hour forecast window, the records from 2015 to 2020 are used to calculate the forecast error.
Fig. 2a shows the frequency of forecast errors of the ambient air temperature for all data, without distinguishing seasons and times of day.The fitted normal distribution has a skewness of -0.67 and a kurtosis of 4.46, which supports the use of a normal distribution to represent the forecast errors [46].The fitted normal distribution has a mean value  amb = −0.066• C and a standard deviation of  amb = 1.61 • C. To simplify the derivation of chance-constraints to be discussed later, the mean value is approximated with 0 • C.
It has been observed that the solar radiation forecast error varies significantly with season and time [47].The solar radiation forecast error is much lower in winter than in summer and also varies during the day.Therefore, the solar radiation data set is split into four seasons and each season is further divided into data sets containing the individual hours of the day, resulting in 96 data sets (four seasons and 24 hours per day).The winter season (December-February) is used in this work.Fig. 2b presents the standard deviation of the forecast errors of the solar radiation throughout the day.The forecast error's standard deviation is zero at night and follows a sinusoidal pattern throughout the day, peaking at noon.Except for 07:00, 08:00, and 17:00, the skewness and kurtosis of the forecast error distribution indicate that a normal distribution is appropriate to represent the forecast errors.During sunrise (07:00, 08:00) and sunset (17:00), the kurtosis exceeds its threshold because the forecast error is mostly zero.This is because the time at which the sun rises and sets is constantly changing.Nevertheless, a normal distribution is still used to simplify the derivation of the chance constraints.The mean values of the solar radiation forecast error are in the range between −46 W∕m 2 and 1.3 W∕m 2 and they are approximated with zero, similarly to what has been done for the ambient air temperature.
The uncertainty of air temperature causes the uncertainty of heat pump COP ().For the air source heat pump considered in this work, its COP is modeled as a third-order polynomial function of the ambient air temperature based on the data points provided by the manufacturer (Fig. 3).Fig. 3 shows that  depends nonlinearly on the ambient air temperature which causes the uncertainty of  to increase with the ambient air temperature  amb .To account for this changing uncertainty, we derive the statistical distributions of the uncertainty of  for various ambient air temperatures from −15 • C to 15 • C with a step of 0.1 • C.This enables the MPC to estimate the uncertainty of  based on the temperature forecast.For all fitted normal distributions, the kurtosis is below 3.3 and the skewness below 0.46, supporting the use of normal distributions [46].The mean values are close to zero (| COP | ≤ 0.0052) and thus approximated with  COP = 0. Fig. 4 presents the derived probability density distributions for the uncertainty of  for ambient air temperatures of −15 • C and 15 • C.
Due to the lack of data, the uncertainty of internal heat gains is simply modeled as a bounded normal distribution with the mean value  internal = 0, the standard deviation of  internal,t = 1 6 Φ internal,t , the bounds in between − 1 3 Φ internal,t and  within a day.After sunrise, the uncertainty in the solar radiation forecast increases until 12:00 and subsequently decreases until sunset (see Fig. 2b).In contrast, the uncertainty in ambient air temperature remains constant throughout the day.The COP, , for the heat pump behaves similarly to the ambient air temperature as they are correlated.However, the uncertainty of  varies with the ambient air temperature and is slightly higher in high-temperature conditions.For instance, at 14:00 when the ambient air temperature reaches its peak of 8 • C, the standard deviation range of  is slightly broader in conditions with lower ambient air temperatures.Similarly, the uncertainty of the internal heat gains forecast is assumed to depend on the forecast values.A high value of Φ internal has a high standard deviation interval.

System identification of the building thermal model
Real-world data are utilized to identify the building thermal model presented in Section 2.1.The data are gathered in a small two-story low-energy research facility (100 m 2 ), which is located in Karlsruhe in Germany [48].The building has five rooms, a kitchen and a bathroom.Air temperature measurements were collected in each room in the building.For the identification of a single-zone model, the measurements were combined using a weighted average.The measurements from each room are weighted by the room volume and the weighted average is used for the identification of a single-zone model.
The training data consists of three weeks of data from December 20, 2021, to January 09, 2022, when occupants were absent for the Christmas holiday and active system excitations could be conveniently performed.The identified model has been validated on one week of operation from January 10 to January 16 in 2022, when occupants returned to the office.Fig. 6 compares the measured space air temperature and 24-hour-ahead simulation results from the identified thermal building model.The validation yields satisfying results with a root mean square error of 0.2 • C.

Performance metrics
Two key performance metrics (KPIs) are used to evaluate the performance of the different controllers.These two metrics reflect conflicting goals that need to be reconciled.The first metric evaluates the degree of violations to the predefined thermal comfort for occupants and it is defined as follows: The above equation shows that the magnitude of thermal comfort violation is defined as the sum product of the deviation from the temperature constraint (Eq.( 8)) and the time duration of the constraint violation.
The second metric evaluates the energy cost:

Scenarios for the comparison
Four different controllers (i.e., PB, DMPC, SMPC, RMPC) are compared on three different days.The three days are selected to have different levels of forecast errors.Considering that the ambient temperature forecast error is the main driver of temperature constraint violations, it has been selected as the decision criterion to select days with high, medium, and low forecast errors.Table 1 presents the three days selected for simulation and their daily average forecast error of the ambient air temperature.The negative sign of the average error indicates that the forecast overestimates the actual ambient air temperatures.Only days with an overestimation are selected because they represent a critical scenario to challenge the controller in the heating season.
For all three days, the same electricity price signal and space air temperature setpoint schedule are used across all simulations.We assume  an electricity tariff Ω buy where the prices for the next 24 h are known without uncertainty.For Ω buy , we use the price signal from [49], representing the wholesale electricity price in Germany on January 13, 2022.Since this price does not reflect what a typical customer pays, it is rescaled so that the mean price of the signal matches the actual electricity retail price in Germany.The tariff Ω sell is assumed to be constantly 0.086 e∕kWh which is the current financial compensation for selling PV power to the grid in Germany [39].

Parameter settings
The MPC problem is formulated and implemented in a generic manner.Therefore, the MPC problem can be instantiated with different parameter values.For this work, the time horizon for prediction and control  is 24 hours and the sample time Δ is 15 minutes.Table 2 lists the parameter values for the heat pump, PV, battery storage, and domestic water heater.

Table 1
The three days for simulation with different levels of forecast errors.For deferrable electric appliances, the present work uses the power consumption profiles from [41], where the power consumption of the dishwasher and the laundry machine is structured in phases.The maximum delay times in between phases were derived from [50].All necessary numerical values are presented in Table 3 to enable reproducibility.

Implementation
Fig. 7 outlines the MPC implementation for the present simulation study.Uncertain forecasts are denoted as ( ⋅).At every time step , the MPCs obtain 24-hour-ahead weather forecasts for ambient air temperature ( Tamb ) and global solar radiation ( Φglobal ).The heat pump coefficient of performance ε and the internal heat gains Φinternal are predicted based on Tamb and the assumed occupancy schedule, respec- tively.
Solving these optimization problems yields an optimal input sequence (i.e.,  HP ,  ch ,  d ,  el ,  i j ) for times , ...,  +  .However, only the first input (for time ) is applied to the building model.The building model considers real-time values of the uncertain variables. amb and Φ global are obtained from weather measurements,  is calculated from  amb , and for Φ internal , a random forecast error is drawn due to a lack of measurements for Φ internal .The whole optimization process is repeated rollingly with a sample time of Δ, ensuring continuous optimization of the building's performance based on real-time and forecast data.The forecast and real-time data for the uncertain variables are presented in Appendix A.1.
The mixed-integer linear program is implemented in MATLAB R2021b with the help of the toolbox YALMIP [51] and solved with CPLEX 12.10 [52].The optimization is conducted on a computer equipped with an Intel Core i7-10700 (2.9 GHz) and 32 GB RAM.The computation time for a single time step of 15 min is around 30 sec, most of which attributes to the appliance model.

Results of equipment operation scheduling
The cost-saving strategy of the MPCs is illustrated using the results of DMPC for the day of December 19, 2017 with a low level of forecast errors.The strategy is similar for the other MPCs and days.Thus, it suffices to consider this exemplary day and controller in detail.
First, we investigate the heat pump operation.Fig. 8 presents the profiles of the heat pump power consumption, the indoor air temperature, the lowest comfortable air temperature and the electricity tariff under the DMPC for the day of December 19, 2017, when a low level of forecast errors occurred.Apparently, the controller has managed to heat up the building preemptively during time steps with low electricity prices to reduce necessary heating in subsequent time steps with high electricity prices.For example, at 03:45, the indoor air temperature is 18.2 • C.Even though the indoor air temperature is higher that the heating setpoint, the heat pump is controlled to run until 05:00 to heat up the indoor air temperature to 20 • C. The preheating is needed because both the heating setpoint and the electricity price are increased in the forthcoming time steps.Thus, the MPC activates the heat pump at appropriate time steps when the electricity price is low, enables the use of building thermal mass to store energy, and releases the stored energy to the space when the electricity price is high to respect the temperature constraints.This phenomenon can be demonstrated by another time period from 12:00 to 15:00 when the controller preheats the building and the electricity prices are low.
The operation of the battery is presented in Fig. 9. Fig. 9a depicts the charging and discharging power of the battery and the electricity price over the simulation period, and Fig. 9b the resulting state of charge of the battery.The MPC charges the battery from the grid when the electricity price is low and discharges the battery to operate the heat pump when the electricity price is high.For example, from 12:30 to 14:00, the battery is charged when the electricity price is low because the MPC anticipates the electricity price to rise in the future.From 17:00 to 20:00, the MPC discharges the battery to operate the heat pump for space heating and domestic hot water.This strategy prevents purchasing electricity when the electricity price is increased.
The scheduling of deferrable electric appliances follows a similar logic as that for the heat pump.All three appliances are scheduled to run between 07:00 and 22:00.The MPC controller activates them during the time steps with the lowest electricity prices (see Fig. 10).Specifically, the dishwasher starts its first phase at 12:00 and finishes all four phases at 14:00.There is no break between the phases even though time delays are allowed (see Table 3).Similarly, the laundry machine starts its first phase at 11:15 and continues its operation without interruptions since no financial benefits can be achieved by suspending a phase.At 12:45 when the laundry machine completes, the enforced delay between the laundry machine and the dryer can be observed.This delay allows the occupant to move the clothes from the laundry machine to the dryer.The clothes dryer starts at 13:00 and it continues without interruption because it has only one phase.The scheduling results of the heat pump, battery, and electric appliances are derived from the optimization after considering multiple factors such as the electricity prices, the PV power generation, and the state of charge of the battery storage.Therefore, it is worthwhile to explore the interactions between the whole building and the power grid.Fig. 11 shows the purchased power, the sold power, and the electricity tariffs.This figure further demonstrates that power is usually purchased when the price is low.The purchased power peaks around 13:00 due to the operation of the appliances and the domestic hot water.The controller chooses to prioritize self-consumption and does not sell any power to the grid, as the financial compensation is much lower than the cost of purchasing electricity.

Performance comparison of MPC variants
All three MPCs (i.e., SMPC, RMPC, and DMPC) together with the theoretical reference case of PB are simulated for three days with different levels of forecast errors.The results are presented in Table 4.For the SMPC, a value of  = 0.1 is selected, corresponding to a 90% probability of temperature constraint satisfaction.The simulation results are compared using thermal comfort and electricity cost, the two metrics described in Section 5.2.Based on reference [7], the maximum acceptable comfort violation is 70 Kh∕year.If the comfort violation is assumed to be evenly distributed over 365 days of a year, the comfort violation should satisfy Δ comfort ≤ 0.19 Kh∕day.low forecast errors, but results in a slightly higher electricity cost than DMPC and SMPC for the day with high forecast error.DMPC and SMPC underestimate the heating load due to the high forecast error, resulting in insufficient heating to meet the temperature constraints.Thus, they lead to a reduced electricity costs compared to PB, but at the expense of temperature constraint violations.• Except for PB, the DMPC achieves the lowest electricity cost.The comfort violation is acceptable (i.e., Δ comfort ≤ 0.19 Kh) for the days with low and medium forecast error.For the day with high forecast error, the comfort violation is precisely on the threshold in between acceptable and unacceptable.
• The SMPC made a tradeoff between the two performance metrics.
Relative to the DMPC, the SMPC has a negligible increase (0.2%-0.3%) of electricity cost but reduces the comfort violation by 50% for the day with high forecast errors and to almost zero for the days with medium and low forecast error.For all days, the comfort violation is well below the limit of 0.19 Kh. • Similarly to PB, RMPC achieves optimal thermal comfort.However, SMPC also provides satisfactory thermal comfort and results in 1.1%-1.3%lower electricity costs.Therefore, between the controllers that explicitly consider uncertainty, RMPC and SMPC, SMPC appears to be more cost-effective.
Table 4 shows that the DMPC has the highest thermal comfort violations among all studied controls.However, the degree of comfort violations is low, which may make it difficult to select between SMPC  and DMPC because the former entails significant efforts to implement.The low comfort violations of DMPC might be caused by the following factors: • The considered building is highly insulated and energy efficient.Thus, its thermal resistance is considerably high, which that small forecast in the ambient temperature has on the room temperature in any given time step.• In winter, the forecast error solar radiation is significantly low, as observed in our data analysis for uncertainty quantification.• The relaxed temperature constraints during daytime also contribute to low comfort violations.It has been assumed that the occupants are not at home during working hours and the indoor air temperature constraints are set accordingly.When the occupants are assumed to return from work at 17:00, the sun has set, leading to a perfect solar radiation forecast of 0 W∕m 2 .Thus, the relaxed temperature constraints, together with an early sunset reduce the likelihood of compromising the thermal comfort constraints.
To investigate the impact of building envelope performance on the comparison of MPCs, a poorly insulated building is emulated by reducing the resistances of the thermal building model by 50% (see Table A.6 for the parameters of the thermal building models).The simulation is rerun for the poorly insulated building model and Table 5 provides the results.
By comparing Table 5 with Table 4, we have the following observations: • The electricity cost increases significantly by 50%-70% for all cases because more heating is required for the poorly insulated building.• The ranking of the three MPCs with respect to their performance on the two metrics stays the same between the original model and the poorly insulated building model.The DMPC achieves the lowest cost but with the highest thermal comfort violation while the RMPC achieves the lowest thermal comfort violation but with the highest cost.The performance of SMPC lies in between the other two MPCs for both metrics.• For the DMPC, the comfort violation increases by 142%, 162%, and 31% respectively for the days with high, medium, and low forecast error.This results in compromised thermal comfort for the days with medium and high forecast error since the comfort violations exceed the acceptable threshold of 0.19 Kh.Thus, thermal comfort in poorly insulated buildings appears vulnerable to forecast uncertainties, and it seems necessary to explicitly consider these uncertainties.• RMPC entirely avoids comfort violations for all days.The cost difference between the RMPC and the DMPC increases with the declining level of forecast errors.Compared with the DMPC, the RMPC increases the electricity cost by 2% for the case with the highest forecast error to 4.5% for the case with the lowest forecast error.The above difference can be explained by the hedging of the RMPC against the worst-case scenario (i.e., the case with the most unfavorable deviations from the forecasts).As the forecast error decreases, the possibility of incurring the worst-case scenario becomes very low.• Compared to the DMPC, the SMPC reduces the comfort violation by 43% and 62% for the high and medium forecast error days, while the electricity cost increases negligibly by 0.3%.Consequently, the SMPC achieves satisfactory thermal comfort for the low and medium forecast error days, but the thermal comfort is compromised for the high forecast error day.
One particular advantage of the SMPC not shown in Table 4 and Table 5 is its flexibility.The "weight" of thermal comfort is tunable by the probability of constraint violations 1 −  (see Eq. ( 37) and (38)).
For example, if the comfort violation of the SMPC is regarded too high, the  value can be decreased.Similarly, the  value can be increased to prioritize cost optimization.To illustrate the use of  to balance the cost and thermal comfort, simulations with different values of  are performed for the poorly insulated building and the day with high forecast error.Fig. 12 shows the results together with the comparison with PB, DMPC, and RMPC.Fig. 12 shows the SMPC's trade-off between electricity cost and thermal comfort violation can be adjusted by modifying  to achieve satisfactory thermal comfort for the day with high forecast error.As  decreases, the SMPC reduces comfort violation at the ex- pense of increased cost.For  = 0.001, the comfort violations are close to zero, while the electricity cost is still e0.16lower than the RMPC's.Thus, the RMPC results based on the worst-case scenario seem overly conservative to satisfy the comfort constraints while leading to high electricity cost.If  is set at 0.5, the reformulated chance constraints simplify into their deterministic equivalent [20] since Φ −1 (0.5) = 0 (see Eq. ( 41) and ( 42)).Thus, the SMPC with  = 0.5 will coincide with the DMPC.The appropriate value of  depends on the level of forecast accuracy, where a high forecast error requires a lower  to satisfy the comfort requirements.addition, building degraded envelope performance likely lower value of  than a high-performance building.If  is set too low, SMPC will be unnecessarily conservative, while a too high  will compromise thermal comfort.A dynamic adjustment of  based on currently observed forecast errors, known building insulation characteristics, and occupant preferences could further improve SMPC and is the subject of future work.

Conclusion
This paper presents a comparison study of deterministic, robust and stochastic model predictive controls for residential buildings with distributed energy resources, including PV, battery storage, and controllable loads (e.g., heat pump, domestic hot water and smart electric appliances).All considered MPCs are formulated and implemented in a generic manner, including the thermal resistance and capacitance method for building thermal modeling and the physical models for different distributed energy resources.Two thermal building models are considered, one identified with real data from a research building and a modified thermal building model to emulate a poorly insulated building.The uncertainties associated with weather and internal heat gains are quantified and considered in the SMPC and RMPC.The SMPC is reformulated into a tractable representation via the use of chanceconstraints.
The comprehensive set of distributed energy resources is successfully implemented in all MPCs.It contributes to reducing the cost function by shifting the power consumption to times with low electricity prices, demonstrating the potential for demand flexibility in buildings.
The simulations are conducted for three days with different levels of forecast error.For the highly insulated building, the DMPC performs well under medium and low forecast errors because it achieves satisfactory thermal comfort and low electricity costs.However, when the forecast error is high or the building is poorly insulated, the DMPC results in compromised thermal comfort.
The SMPC has a slightly (less than 1%) higher electricity cost than the DMPC and meanwhile it achieves satisfactory thermal comfort for all cases if the probability of constraint satisfaction is selected appropriately.It provides a flexible trade-off between energy cost and thermal comfort that can be adjusted by modifying the probability of constraint satisfaction.In the case of the poorly insulated building with a high forecast error, a higher constraint satisfaction probability is required to achieve satisfactory thermal comfort.
The RMPC has up to 4.5% higher energy cost than the DMPC, but never violates the temperature constraints.It tends to be overly conservative and less cost-effective than the SMPC in improving thermal comfort, especially on days with low forecast error.
In conclusion, the DMPC is well-suited for an application in highly insulated buildings when the forecast uncertainty is low or moderate.However, explicit consideration of forecast uncertainty is necessary for poorly insulated buildings or high forecast error.The strengths of RMPC are the ease of implementation and complete constraint satisfaction for a cost increase of up to 4.5%.The SMPC provides satisfactory thermal comfort for a negligible cost increase of less than 1% and is the preferred option for achieving optimal cost-effectiveness.
For SMPC, the choice of the appropriate constraint violation probability  is an open question that could be addressed by dynamically adapting  during operation based on the monitored comfort violations and forecast errors.Such an investigation is foreseen for future work.Furthermore, real-world experiments could be conducted to explore the differences between multiple controllers.In the present work, a predetermined electricity tariff was assumed.Future research could investigate the effects of uncertain electricity prices on MPC performance and explore strategies for managing price uncertainty.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.Please note that the ambient air temperature is used as decision criterion to classify the days as high, medium, or low forecast error.Thus, the forecast error in the solar radiation does not necessarily correlate with the classification.

Fig. 2 .
Fig. 2. Forecast error of the solar radiation and ambient temperature with derived error distributions.

Fig. 4 .
Fig. 4. Probability density of the forecast error of the heat pump coefficient of performance.

Fig. 5 .
Fig. 5. Forecast of the uncertain variables and the interval for one standard deviation.

Fig. 8 .
Fig. 8. Indoor air dynamics and power consumption of the heat pump under the DMPC for the day of December 19, 2017.

Fig. 9 .
Fig. 9. Battery schedule and state of charge under the DMPC for the day of December 19, 2017.

Fig. 10 .
Fig. 10.Appliance scheduling of the DMPC on the day of December 19, 2017.

Fig. A. 14 .
Fig. A.14. Forecast and measured ambient air temperature for the three scenarios.

Figs
Figs. A.13-A.15  present the inputs of the internal heat gains, ambient air temperature, and solar radiation, respectively.Please note that the ambient air temperature is used as decision criterion to classify the days as high, medium, or low forecast error.Thus, the forecast error in the solar radiation does not necessarily correlate with the classification.

Table 2
Parameter values for the building thermal model, PV modules, battery storage and domestic hot water heater.

Table 3
Parameter values for the dishwasher, the washing machine and the clothes dryer.

Table 4
Comparison of performance metrics from different MPCs.
Table 4 shows the following: • Across all days, PB achieves perfect comfort satisfaction because it assumes perfect forecasts of weather and internal heat gains.PB incurs the lowest electricity cost for the two days with medium and F. Langner, W. Wang, M. Frahm et al.

Table 5
Comparison of performance metrics from different MPCs for the poorly insulated building.

Table A . 6
Parameters of the thermal building model for the highly insulated and poorly insulated building.