Scenario-based Nonlinear Model Predictive Control for Building Heating Systems

State-of-the-art Model Predictive Control (MPC) applications for building heating adopt either a deterministic controller together with a nonlinear model or a linearized model with a stochastic MPC controller. However, deterministic MPC only considers one single realization of the disturbances and its performance strongly depends on the quality of the forecast of the disturbances, which can lead to low performance. In fact, inadequate building energy management can lead to high energy costs and CO$_2$ emissions. On the other hand, a linearized model can fail to capture some dynamics and behavior of the building under control. In this article, we combine a stochastic scenario-based MPC (SBMPC) controller together with a nonlinear Modelica model that is able to provide a richer building description and to capture the dynamics of the building more accurately than linear models. The adopted SBMPC controller considers multiple realizations of the external disturbances obtained through a statistically accurate model, so as to consider different possible disturbance evolutions and to robustify the control action. To this purpose, we present a scenario generation method for building temperature control that can be applied to several exogenous perturbations, e.g.\ solar irradiance, outside temperature, and that satisfies several important stastistical properties, in contrast with simpler and less accurate methods adopted in the literature. We show the benefits of our proposed approach through several simulations in which we compare our method against the standard ones from the literature, for several combinations of a trade-off parameter between comfort and energy cost. We show how our SBMPC controller approach outperforms the standard controllers available in the literature.


Introduction
Energy consumed in buildings for heating, ventilation, and air-conditioning (HVAC) purposes accounts for around half of total energy used in buildings [1][2][3][4][5]. For companies, especially if they are located in large buildings, it is therefore very important to limit the amount of energy wasted due to bad temperature control. Furthermore, it is also important to reduce as much as possible the energy waste in order to decrease emissions due to e.g. natural gas boilers [6]. Moreover, buildings have comfort temperature bounds that have to be respected during working hours, with few violations allowed [7]. The comfort violations should be limited while at the same time the energy cost has to be minimized. Simple solutions where heaters always run at maximum power are not acceptable due to the large costs and energy waste.
On top of that, buildings are subject to many exogenous disturbances, e.g. outside temperature, solar irradiance, and endogenous ones, e.g. occupancy. The profile of these disturbances, if not properly considered when determining the control actions, can lead to poor performance, both in terms of energy cost and discomfort. Moreover, the model of the building considered, due to approximations, may lead to further errors. However, a too complex model is also not useful for control purposes due to the high computational burden that it entails. Therefore, the problem of controlling the room temperature in large buildings is a complex task.

Literature review
In this section, we perform a brief literature review of the two main topics of this research: control for buildings and scenario generation.

Modeling and control strategies for buildings
Many solutions have been proposed in the literature to control the room temperature in buildings using information available on the current temperature and forecasts of the disturbances. A simple solution involves a rule-based approach that is based on if-then-else rules and information about the current disturbances [8]. Although these controllers are simple to implement and may achieve a reasonable performance, they are not very efficient since they are based on user knowledge and rules-of-thumb and they do not actually perform an optimization. Model predictive control (MPC) [9][10][11][12] is a more advanced control strategy that is suitable for the room temperature control problem, since it can naturally include constraints in the control problem and since it can in general achieve a better performance [13][14][15]. Moreover, a building is subject to several disturbances, as mentioned before, and MPC can deal well with disturbances by using a robust or stochastic controller [11], which can achieve better performance than the deterministic counterpart. However, despite having a better constraint satisfaction, a robust MPC controller for buildings, e.g. [16], could be too conservative for the task of controlling the temperature of a room and it would result in a high amount of energy used as it would try to compensate every possible disturbance realization. Therefore, usually a stochastic MPC controller is preferred for building heating systems [12,17]. Indeed, by considering the stochastic properties of the disturbance or by considering several disturbance scenarios, stochastic MPC controllers can potentially achieve a better control performance compared to deterministic controllers, leading therefore to a reduced energy consumption while limiting the discomfort. In particular, scenario-based MPC (SBMPC) methods stand out as a useful tool in building heating systems, since they can use past data of the disturbances, which are available in the case of building heating systems, and they can successfully be applied to nonlinear models as well [12].
In this regard, several types of MPC algorithms have been applied to HVAC systems in the literature [13][14][15][17][18][19][20][21][22][23][24][25][26][27][28][29]; see also [30,31] and the references therein. In particular, [15] presents two stochastic MPC algorithms, i.e. a disturbance-feedback approach and a chance-constraint one. The results show that the stochastic controllers achieve a better performance than deterministic MPC and rule-based control. Authors of [17] develop an SBMPC controller that uses previous data of building occupancy and external ambient conditions forecast errors to solve a scenario-based optimal control problem. The scenarios are built through copulas that can be learned online and the method is applied to a room in a university building. The results show that the proposed controller is able to achieve a good performance in terms of energy cost, while having a larger computational complexity than standard deterministic methods. However, the method is applied to a single room and the authors of [17] suggest that a study needs to be carried out to asses whether the method can be applied to a larger building. A similar approach is presented in [18], where the main differences are that 1) slack variables are introduced in the cost function to improve the feasibility of the optimization problem obtained and 2) that different copula families are tested and compared (see also Section 1.1.2). In both papers, the authors do not rely on Gaussian assumptions for what concerns the properties of the disturbances. Given that the two previous methods might result in a large computational complexity, the authors of [19] extend the concept to an explicit SBMPC controller, such that the control inputs are computed offline and applied online. Experiments were performed once again for a single university room, showing better performance with respect to the standard methods. In order to deal with a multi-room setting, a distributed MPC controller is presented in [20], where a Lagrangian dual decomposition relaxation method is used to reduce the computational burden arising from the several rooms considered. Simulation results obtained considering a network of 10 rooms show an increased performance with respect to a baseline PID controller. In all these articles it is shown how stochastic MPC strategies can achieve a better performance than deterministic MPC. The article [21] presents an MPC algorithm in which a linear model is used to control a building including an active cold thermal storage in order to implement a demand response program. All these works, i.e. [15,[17][18][19][20][21], use a linearized model description and do not use a nonlinear model nor other more advanced modeling tools, e.g. Modelica [32,33], possibly leading to a decrease in the performance. Such tools and nonlinear models for building heating control usually include more features and components compared to a linear model and thus they can potentially provide a more accurate description of the building and of the influence of each disturbance, reducing therefore the modeling error and improving the overall performance. The article [22] adopts a nonlinear model arising from the heat pump and a battery inverter considered, but the considered MPC controller is a deterministic one. For what concerns nonlinear and more advanced models, while some works did consider their usage for HVAC systems, e.g. [13,14,23,24], all of them considered a deterministic setting instead of a stochastic one. To the best of our knowledge, no work has considered a nonlinear model description obtained through Modelica together with a stochastic controller, which would improve the performance by taking into account a more accurate model and the stochastic properties of the disturbances.

Scenario generation
Besides improving the state-of-the-art by proposing a control approach for more realistic models, i.e. nonlinear Modelica models, our work also contributes to the existing literature of scenario generation for buildings by improving upon the current state-of-the-art. In particular, scenarios of random variables that represent a time series, e.g. the ambient temperature for the next 24 hours with an hourly resolution, need to satisfy several important properties: 1. They should not be restricted to the standard assumption of Gaussian disturbances/forecasting errors as this assumption is quite restrictive when it comes to generating scenarios of heteroscedastic 1 processes, e.g. solar irradiance [19].
2. They need to consider the multivariate distribution of the random variables: if the scenarios represent a random variable at different time steps in the future, these scenarios should model the time correlation of the random variable [34]. 3. Besides modeling the time correlation, they should explicitly take into account the different time dependencies and avoid modeling a stationary distribution; i.e. the distribution of the random variable might be different at different hours of the day/times of the year or might change with the prediction horizon. 4. The methods to generate scenarios should be flexible enough to explicitly model the dependencies of the random variables on exogenous variables. 5. The computational burden of the scenario generation method should be small enough for online implementation.
In the context of building heating, while some scenario generation methods have been considered [15, 17-20, 25, 35-40], they have several problems. In particular, some of the existing methods [15,35,36] rely on the standard Gaussian assumptions [19]. In addition, although several works have addressed the Gaussian assumption [17-20, 25, 37, 39, 40], they still lack some of the required properties.
More specifically, in [17], a method based on empirical copulas is proposed. While the method satisfies some properties, e.g. time correlation, it fails to satisfy the following two: 1) it does not model time dependencies but it assumes that the marginal distributions are stationary, i.e. it assumes that the n-hours ahead distribution of a variable is the same at any hour of the day, any day of the year; 2) the scenarios are generated based on historical data without considering other possible exogenous inputs. The analytical copula method proposed in [20] overcomes some of these issues as it explicitly models the time dependency during a day. However, the distributions are still stationary, i.e. they vary within a day but they do not change along a year, and they are just based on historical data. In [18] and [19], a more general approach is proposed where different copula families are tested, and the best one is selected to generate scenarios. While the method is very general and flexible, it requires to compare different copula functions and can easily become computationally infeasible for online MPC. In addition, the method has two other problems: the best copula is selected by comparison with the empirical copula of [17], hence it has the same problems as [17]; moreover, the time dependencies considered by the copulas are not specified. In [37], scenarios from a weather meteorological model are employed. Even though the goal of weather models is to provide an ensemble of scenarios, to capture the uncertainty in the prediction, the method displays systematic errors, e.g. biases, and requires the application of advanced post-processing techniques based on copulas, e.g. ensemble copula coupling [41]. In [25], a method based on sampling historical forecasting errors is considered. Despite the method attempts to capture time correlation using an auto-regressive error model, that model is only used for error correction. In particular, to generate scenarios, the method samples from past historical errors and fails to satisfy properties 2-4. The recursive feasibility and stability of SBMPC is studied in [39]. To do so, it is assumed that disturbances can be represented by a scenario tree, and that the tree can be built using empirical samples from a discrete set of disturbances. This approach is clearly very limited as, it does not satisfy properties 2-4, and in addition it could have scalability and computational issues when the number of random variables increases. In [40], scenarios are used for modeling electricity prices and independent optimization problems are solved for each scenario; however, the method cannot be used to model uncertainty in other variables, e.g. weather variables, and fails to satisfy properties 2-4. In [38], two Poisson distributions are employed to model the occupancy in the building as a birth-death process. While such a parametric distribution might work well for occupancy, it has the same issue as the Gaussian assumption: the method cannot be generalized to other random variables.

Motivation and contributions of the paper
In this article, we focus on a scenario-based MPC (SBMPC) algorithm that includes both a nonlinear system description through Modelica, while using a scenario generation method based on probability distributions obtained empirically, making it a very suitable tool for a building heating control problem. The Modelica model description can lead to an improvement of the model accuracy; note that in the current literature of SBMPC for heating systems in buildings, a linearized model is always used. On top of that, by using a nonlinear Modelica model, we implement for the first time an SBMPC controller for HVAC systems in buildings that uses a nonlinear Modelica model description. In addition, we present a parametric Gaussian copula method to generate scenarios that it can satisfy the five required properties, whereas the methods used in the literature of HVAC systems for scenario generation suffer from some statistical problems, as mentioned earlier. We perform several simulations showing the benefits of the presented approach. Our contribution is therefore threefold: • We propose a control method for a building heating systems that considers a Modelica nonlinear model and an SBMPC controller.
• We generate scenarios using a new approach that, unlike the existing methods, satisfies all the important properties of scenario generation methods for time series data.
• We perform a comparison between several combinations of the couple model-controller: as models, we consider a Modelica and linear and as controller we consider a deterministic MPC and SBMPC.

Outline
The outline of the article is as follows. In Section 2, we describe the problem under consideration. We present the control algorithms used throughout the article in Section 3. In Section 4, the adopted scenario generation method is presented. We present the simulation results on a case study in Section 5 and lastly we present some conclusions and suggestions for future work in Section 6.

Buildings
In this paper, we focus our attention on buildings with local heat production units. The type of building considered can be controlled via two control inputs, i.e. u = q heat q cool ⊤ , where q heat is the amount of heating power transferred to the building and q cool is the cooling power provided to the building. We assume that the building can be modeled using an RC-model with two states [13]: T zone as the average temperature of the rooms and T wall as the average temperature of the walls. In addition, it is assumed that the building is affected by three disturbances: solar irradiance I, outdoor temperature T amb , and building occupancy θ occ . While past measurements of external disturbances, e.g. solar irradiance and outdoor temperature, are available, we do not have any measurement of the occupancy of the building. Note that, although this is an important disturbance to consider, it is also difficult to measure in practice [42,43]. Therefore, to estimate the models and to perform simulations, we assume that the occupancy profile is fixed for every day of the week, i.e. we assume that the building is fully occupied during working hours and empty outside of these hours.

Modelica
We have modeled the buildings, comprising also the heating, cooling, and ventilation units, with Modelica [32,33], which is an object-oriented and equation-oriented language that is designed to model the behavior of physical systems. In particular, the building is modeled based on an RC-model, which has been identified through the Grey-Box Buildings toolbox [44]. The building has also been extensively validated using data collected from the building as in [13,44]. Such data has been gathered between 2016 and 2018 and includes e.g. internal temperature, domestic hot water usage, external temperature, and solar irradiance. The adoption of Modelica in our work provides the benefit that we can improve the amount of detail and accuracy of the model w.r.t. linear models. Note that e.g. some of the HVAC components modeled in Modelica result in nonlinear model components. Therefore, although many other works, e.g. [15,[17][18][19][20], use indeed a linearization of a nonlinear model, in this work we directly use a nonlinear model and we obtain therefore a more meaningful representation of the real building. Readers interested in the modeling procedure of buildings in a Modelica environment are referred to [13,44].
Note that other high-fidelity simulation tools exist for buildings, e.g. TRNSYS, EnergyPlus, ESP-r, IDA ICE; see [45,46] for a complete review. However, compared to Modelica, these software tools lack in modularity and flexibility for prototyping and simulating new technologies [47]. Moreover, in our particular setting, we have data available from the buildings to be controlled that allows us to estimate the parameters of the building. However, the data is not comprehensive enough for satisfactory whitebox modeling. Therefore, as highlighted in [44], we adopt a grey-box model due to necessity of estimating certain parameters of the model of the building that are not known a priori. Whereas in the aforementioned white-box modeling tools it could in theory be possible to calibrate the parameters of the model, the lack of a more detailed physical knowledge of the building makes the choice of a whitebox model, for our specific setting, less desirable. Lastly, note that Modelica is an open-source tool, making it particularly appealing for commercial applications. For more advantages of using Modelica for HVAC systems we refer to [48].
Remark 1. Note that, as mentioned in [44], the actual difference between white-box and grey-box models does not depend on the complexity of the model. For instance, even a single-state model can be a white-box model if all its parameters can be determined based solely on physical knowledge. However, a white-box model becomes grey when one or more of its parameters are estimated based on a fitting of the model to measurement data, irrespectively of the complexity of the white-box model.
Remark 2. The model of the building used in this paper has been extensively validated following the exact same procedure reported in [44]. Therefore, the validation process is omitted here as it is already explained in the aforementioned paper. Multiple-years data has been used to identify and validate the model. The interested readers might find the whole identification and validation procedure in [13,44].

Model predictive control
MPC is a control tool that has been extensively studied in the last forty years [9][10][11]. The main strength of MPC is to use a model of the system under control to find optimal inputs for that system with respect to a certain objective, which can be either stabilizing the system or minimizing an economical goal. Since in the MPC framework the problem of finding the optimal inputs for the system is an optimization problem, constraints can easily be included, as well as several performance criteria. Compared to simpler strategies e.g. rule-based control, MPC requires a larger computational effort and requires to solve an optimization problem online but it can provide an increase in performance. Due to its versatility and good performance achieved, MPC has been applied to various systems and fields, e.g. power systems [49,50], traffic networks [51], water networks [52], aerospace [53], among others. For a survey that includes also MPC applications for large-scale and industrial systems, we refer the reader to [54].
In recent years, several MPC strategies have been developed to cope with external disturbances, in particular, robust and stochastic MPC strategies. In this paper we focus on an SBMPC algorithm, which is presented in Section 3.2.

Control loop and practical implementation
Many operations have to be carried out on the real building by the building energy management system; the overall control scheme is presented in Figure 1. The operations are [13]: 1. Monitoring: some measurements, e.g. water temperature, heat flux, are performed by the building energy control and management system. 2. Forecasting: weather forecasts are obtained as will be explained in Section 4. 3. State estimation: some states, e.g. internal wall temperatures, cannot be measured and therefore they have to be estimated. 4. Optimization of the control inputs: an optimization problem, explained in Section 3, is solved at every time step with a sampling time of 1h. Only the first inputs of the optimal sequence are applied to the system.
The real building is therefore controlled by all these steps, carried out in sequence. The optimization problem discussed in step 4 above is solved through JModelica.org [55]. The direct collocation method is used to discretize time so that the optimization problem is reduced to a nonlinear programming problem [56]. CasADi [57] is used to obtain the first-order and second-order derivatives of the expressions in the nonlinear programming problem with respect to the decision variables, which are required by the solvers used by JModelica.org. We use IPOPT [58] to solve the nonlinear programming problem, together with the sparse linear solver MA57 [59].
Remark 3. Note that in this section we consider simulations instead of experiments in the real building. This implies a small difference w.r.t. Figure 1: instead of applying the inputs to the real building, we simulate its behavior through Modelica for one time step. In other words, the loop is "closed" by applying the optimal inputs to a model of the building rather than to the building itself, using the actual values of the disturbances instead of the forecasts. The model used for simulating the building behavior is the same Modelica model used for the optimal control problem in the MPC framework.  Figure 1: Scheme of the MPC framework (adapted from [13]).

Linear model estimation
To compare the nonlinear MPC controller with the standard linear counterpart, a linear model of the building is needed. To obtain such a model, data from the building is considered and a linear model is estimated using linear least squares. In detail, considering the same inputs, state space, and disturbances as for the nonlinear model (see Section 2.1), we can assume that the building dynamics are of the form: Then, using the same data as those used for estimating the Modelica nonlinear model, we solve a linear least square problem and estimate the values of the matrices A, B 1 , and B 2 , using the mean absolute error as key performance indicator.
Remark 4. In this article, as highlighted in Remark 3, we use the nonlinear Modelica model to compute the next state for the closed-loop simulations. This can introduce some bias when comparing the performance of the controllers that use the linear model with respect to the ones that use the Modelica model, since there would be no model mismatch in the latter case. Nevertheless, we present here also a linearized model as a reference, to show how such model would behave in the considered scenario.
Remark 5. Related to the previous remark, note also that this kind of nonlinear Modelica model has been already applied to a building heating system in [13]. Indeed, in that reference the same modeling procedure was applied to a similar building and experimental results were carried out. No large mismatch between the Modelica nonlinear model and the actual physical building nor other fundamental flaws were noticed. Therefore, we can argue that the nonlinear Modelica model is a very good approximation of the real building and can therefore be used in the closedloop simulations to simulate the evolution of the system.

Control algorithms and models
In this section we present the two considered control algorithms, i.e. MPC and SBMPC.

Deterministic MPC
In Deterministic MPC, the external disturbances, e.g. temperature or solar irradiance, are predicted with a point forecasting technique and in which the predictions are then assumed to represent the expected value. In this context, at each time step, the MPC optimization problem is solved, yielding an optimal control input sequence u * . Then the first element of the sequence is applied, the horizon is moved one time step forward, the system is sampled, and the optimization problem is solved again.
Given the task of controlling the room temperature in a building while minimizing both the energy costs and the discomfort, the optimization problem solved at each time step by a deterministic MPC controller is given by: subject to where: • N is the prediction horizon.
• The system state is defined by , with T zon k and T wall k as the room and wall temperatures.
• T 1 is the current temperature.
• The input control is defined by q k = [q heat k , q cool k ], with q heat k and q cool k as the input heating/cooling power.
• The cost function represents the weighted average between the energy cost J e k and the discomfort cost J d k : and α is the weighting parameter that defines the relative importance of each cost.
• The building dynamics are defined by (2c), where f (·) represents the Modelica model of the building.
• The building is disturbed by some uncontrollable inputs d k = [T amb k , I k , θ occ k ], with T amb k the ambient temperature, I k the solar irradiance, and θ occ k the building occupancy.
• The upper and lower comfort temperature bounds are respectively defined by T max k and T min k , and they vary in time depending on the hour of the day and day of the week.
• Q heat max , Q cool max , η cool , η gas , c gas , and c ele are constant parameters and are, respectively, the maximum heating power, the maximum cooling power, the cooling efficiency, the heating efficiency, the gas cost, and the electricity cost.
It is important to note that the role of the discomfort cost J d is to act as a soft constraint so that it penalizes the deviations of the temperature outside the comfort bounds, but remains 0 if the temperature is inside the bounds. The controller can therefore choose to implement a control action that leads to a violation of the comfort bounds if this can lead to a lower total cost. Lastly, note that in (2a) the final cost is related only to the states. This is standard in the MPC framework (see e.g. [11], Eq. 2.3), since the inputs are applied from k + 1 until k + N , thus the evolution of the system is considered from k + 2 until k + N + 1. Therefore, there is no input cost considered for time-step k + N + 1, but there is a state cost that is the one we obtain by applying the inputs at k + N .

Scenario-based MPC
It is possible to improve the performance of the deterministic MPC of the previous subsection by considering several scenarios of the disturbances acting into the system. This approach, known as scenario-based MPC (SBMPC), considers multiple realizations/scenarios for the disturbances, different system states for each scenario, and a cost function that consists of the average of the original cost functions across all scenarios. For the control inputs, two possibilities exist: different control inputs for each scenario (as with the system state) and shared control inputs across all scenarios. While the former has the advantage of being less conservative, the latter is more computational friendly. For the case of building control, we consider shared control inputs across all scenarios as this reduces the computational complexity.
Defining M different scenarios for the disturbances. i.e.
, the SBMPC optimization problem solved at each time step can be defined as: subject to where: • T k = [T zon k,1 , T wall k,1 , . . . , T zon k,M , T wall k,M ] represents the state at time step k for each of the M scenarios.
, represents the i th disturbance scenario at time step k.
• The cost function is the average across all scenarios of the weighted average between the energy cost J e k and the discomfort cost J d k of each specific scenario.
• The input control q k = [q heat k , q cool k ] remains equal across all scenarios.
• The building dynamics are represented independently for each scenario by (5c).
• The constant parameters are the same as for deterministic MPC.
Remark 6. Note that other stochastic formulations exist that can provide guarantees on the feasibility of the obtained solution, e.g. [19]. In this paper, we adopt the scenario-based formulation (also referred to as multiscenario formulation) as in e.g. [60][61][62][63]. The implementation of other stochastic methods will be investigated as future work.

Linear MPC
We will compare the SBMPC approach against two linear MPC approaches: deterministic linear MPC and linear SBMPC. In both cases, the optimization problems solved at each time step are the same as the ones defined by (2) and (5) but with a minor modification. Instead of using the nonlinear dynamics (2c) and (5c), the dynamics are given by the linear model defined in (1). In particular, for linear deterministic MPC, constraint (2c) is replaced by: Similarly, for linear SBMPC, constraint (5c) is replaced by: for k = 1, . . . , N.

Scenario generation method
In this section, we describe the scenario generation method for modeling the uncertainty in the system disturbances.

Introduction
Let us define a random variable X representing some time series process, e.g. external temperature, and the related multidimensional random variable X representing the distribution of X in a time grid of N time steps, i.e. X = [X 1 , . . . , X N ] ⊤ . To generate scenarios, we will build the multivariate distribution of X, i.e. F (X), so that by sampling from F (X) we can obtain M scenarios of X, i.e. x 1 , . . . , x M .
When building F (X), in order to satisfy the desired properties of scenario generation methods (see Section 1.1.2), several requirements need to be satisfied: • F (X) should not be substituted by the N marginal distributions F (X 1 ), . . . , F (X N ). In particular, F (X i ) only represents the distribution of X at time step i but does not consider the correlation between X 1 , . . . , X N .
• F (X) should not be built as a stationary distribution. Instead, the distribution should consider the properties of the underlying random variable X. For example, in the case of temperature or solar irradiance, it is clear that F (X) should vary with the day of the year d as well as the hour of the day h, i.e. F (X) := g(X, d, h).
• The distribution F (X) should include any external dependency of X. For instance, if X represents the ambient temperature, F (X) needs to explicitly include the dependency w.r.t. factors like the solar irradiance I, i.e. F (X) := g(X, I).

Scenario generation method
The proposed method consists of four steps: 1. generation of a deterministic forecastx of the random variable X.

Generation of the marginal probability distribution
F (X 1 ), . . . , F (X N ) along the horizon N . 3. Generation of the distribution F (X) using a parametric copula and the marginal distributions F (X 1 ), . . . , F (X N ). 4. Sampling of scenarios using F (X).
In this section, we explain the four steps in detail.

Deterministic forecast
To build a deterministic/point forecast of the variable of interest we employ state-of-the-art methods for each variable of interest: • For the solar irradiance, we consider two forecasting models: the deep neural network proposed in [64] for the short-term predictions (anything below 6 hours), and the ECMWF weather forecast [65] for long-term predictions (anything beyond 6 hours). This distinction is made because, in the context of solar irradiance forecasting, machine learning techniques perform better for the short-term horizons, while numerical weather forecasts are more accurate for long-term horizons [64] (see also Remark 7). For what concerns the long-term prediction forecasts, note that weatherbased models are highly complex models based on weather patterns; for practical applications, these forecasts are not produced by researcher but rather purchased from three main providers of weather forecasts worldwide and for our study we purchased them from ECMWF [65].
• For the ambient temperature, considering the recent success of deep learning methods for forecasting energy-related variables [66][67][68][69][70][71][72], we develop a deep neural network that uses as inputs the past values of the ambient temperature (hourly values over the last three days), the ECMWF weather forecast of the solar irradiance on hourly resolution over the forecasting horizon, and the hour of the day and day of the year when the prediction is made. For the deep neural network, we consider a two-hidden layer architecture whose parameters are optimized using hyperopt [66], a Bayesian optimization algorithm. We optimize the number of neurons per layer and the activation function. As a result of the optimization we considered a deep neural network with 240 (first hidden layer) and 135 (second hidden layer) neurons. The network uses the rectifier linear unit as activation function and is optimized using the Adam [67] optimizer with early stopping.
It is important to note that the other steps to generate scenarios are independent of the method employed to generate the deterministic forecast. As such, while we advocate for the use of state-of-the-art methods to obtain the most accurate scenarios, the proposed methodology would work as well with any deterministic forecast.
Remark 7. The splits between horizon regarding the forecast of the irradiance are a well studied problem in the literature (see e.g. [64,68]). Note that the exact split (4, 5, 6 hours) might be specific to the location.

Marginal distributions
To generate the marginal distributions, considering its simplicity yet high accuracy, we employ the method of empirical quantiles [69]. In detail, to generate the marginal distribution F (X) of a variable X, the simplest version of this method consists of four steps: 1. Consider deterministic forecasts of the variable in the past, e.g.x 1 , . . . ,x n . 2. Compute the associated historical forecasting errors of the deterministic forecast, e.g.ǭ 1 , . . . ,ǭ n . 3. Compute the empirical quantiles of the errors and its associated empirical distribution F (ǫ). 4. Model the marginal distributions as the point forecast plus the marginal distribution of the errors: For the proposed approach, the method is modified in order to model non-stationary marginal distributions. In particular, defining X k,h,d as the random variable representing the value of X at time h of day d that is predicted k time steps ahead, the proposed approach estimates each distribution F (X k,h,d ) independently. To do so, the distributions of the errors ǫ k,h,d are independently considered for each time step k, time h, and day d, and the distribution F (X k,h,d ) is estimated accordingly: In addition, the following two considerations are made: • To obtain non-stationary marginal distributions that explicitly model the variability of the distribution along a year, F (ǫ k,h,d ) is estimated using the historical errors of the last 60 days.
• To explicitly model the variability of the distribution with the time step and time of the day, F (ǫ k,h,d ) is estimated using past errors of the deterministic forecasts made for the same time step k and time of the day h. Particularly, F (ǫ k,h,d ) is estimated using the historical errorsǭ k,h,d−1 ,ǭ k,h,d−2 , . . . ,ǭ k,h,d−n .
In other words, assuming that the copula function is known, the multivariate cumulative distribution can be easily obtained if the marginal distributions are known. Using this theorem, to generate scenarios, we employ one of the copulas functions that requires fewer computational time: the Gaussian copula. This selection is done for three reasons: i) the method to generate scenarios should be fast for real time implementation; ii) empirically, we observed the Gaussian copula to be a good fit for the disturbances considered, i.e. ambient temperature and irradiance; iii) the Gaussian copula is a well established method that has been used to generate scenarios for different energy-based applications [70,73].
For the sake of simplicity, we refer to [70] for details on the estimation of the Gaussian copula. Here, we simply outline the main idea of the method, which relies on two random variables transformations: Then, to generate M scenarios of X at time step k, i.e. In particular, they follow the original marginal distributions F 1 (X 1 ), . . . , F N (X N ) and they model the inter-correlation in X = [X 1 , . . . , X N ] ⊤ .
As it was done for the marginal distributions, the Gaussian copula method is modified in order to model nonstationary distributions. In particular, defining X k,h,d as the random variable representing the value of X at time h of day d and predicted with a time step k, the proposed approach estimates the copula function with the two following modifications: • To have non-stationary distributions that explicitly model the variability of the distribution along a year, the copula function is estimated using the historical data of the last 60 days. We show 10 temperature scenarios and 10 solar irradiance scenarios in Figures 2 and 3, respectively.
To evaluate the marginal distributions, we compute their 90% and 80% coverage. That is, we compute the percentage of historical elements that do indeed fall in the interval [5%, 95%] and [10%, 90%]. As can be seen from Table 1, the coverage of the distributions is acceptable. As could be expected, the distribution of the irradiance is more far off, but overall the coverage is within expected errors.
Remark 8. Let us explain here how we capture the correlation between the different variables. Since we generate scenarios using copula functions and given a set of individual variables with their marginal distribution, we can use a copula function to represent the full distribution of all the variables. For that, we use the individual marginal distri-butions and some map/inference process. The individual variables can be either the same variable at different time points, e.g. the temperature at different hours, or different variables altogether. From the perspective of the method it does not make any difference. In both cases we consider the marginal distributions of the individual variables (something that we can compute with historical data), and we map them to a copula function that captures the full distribution of the variables (including correlation).

Properties
As a final remark, we outline how the proposed method satisfies each of the required properties mentioned in Section 1.1.2: 1. As the distribution of the disturbances are modelled with non-parametric quantile functions, the generated scenarios are not restricted to the standard assumption of Gaussian forecasting errors. 2. Since the scenarios are generated using a copula function, the multivariate distribution is explicitly considered and the scenarios include the time correlations. 3. As the marginal distributions are estimated for each hour of the day, time of the year, and time-step, the resulting multivariate distribution is non-stationary and captures all time dependencies. 4. As the point forecast considers external factors, the method is not limited to historical data of the variable of interest. 5. Since the Gaussian copula and the empirical quantile methods have low computational costs, the method is especially suitable for online optimization.
Remark 9. In this section, we presented the properties that scenario generation methods should have and we adopt the method presented in Section 4.2.3, qualitatively comparing it against other scenario generation methods used in the literature for SBMPC. However, a comparative and quantitative study comprehending several scenario generation methods aimed at investigating which one provides the best control performance is out of scope for this work and is left as a suggestion for future work in Section 6.

Case study
We present in this section the simulation results in which we compare 5 different controllers:  • PIMPC: perfect-information MPC, obtained using the values of the measurements of the disturbances as if they were known in advance. It is of course not possible to have the real values of the actual measurements beforehand in practice, but this controller can be used as a benchmark for the ideal theoretical achievable performance.
• DetMPC-Mod: deterministic MPC controller presented in Section 3.1 together with the nonlinear Modelica model.
• SBMPC-Mod: SBMPC controller presented in Section 3.2 together with the nonlinear Modelica model.
• DetMPC-Lin: deterministic MPC controller together with the linearized model presented in Section 2.5.
First, the simulation setup is discussed, then the results and the discussions are presented.

Setup
The closed-loop control is applied as explained in Section 2.4, i.e. the MPC problem is solved and the first input is applied to the system. Then, for all the controllers, the evolution of the real building between sampling times is simulated through Modelica.
We perform simulations for one month in the winter season with 1h sampling time, i.e. we solve 24 · 30 = 720 optimization problems for each controller. The prediction horizon is N p = 24, i.e. corresponding to one day. We consider an office building in Brussels, Belgium, with 7 floors and a total surface of 10000 m 2 . A nonlinear model of the building is estimated using Modelica based on the considerations of Section 2.1 and using data from the real building. In addition, a linear counterpart is also estimated using regular linear least squares as explained in Section 2.5. The heating system consists of 2 gas boilers of 500 kW each and one chiller of 500 kW. We consider thermal comfort bounds that change throughout time, i.e. the lower and upper comfort bounds are set respectively to 21.5 • C and 24 • C during occupation hours and 18 • C and   Table 4: Contribution of the subcosts Je and αJ d to the total closed-loop for all the controllers considered in the case study. The first number represents the percentage of the energy cost Je in the total cost, while the second number represents the percentage of the discomfort cost multiplied by α, i.e. αJ d . This table is also depicted as a stacked bar plot in Figure 4. 26 • C during the non-occupation hours, as shown in Figure  5. Furthermore, the building occupancy profile follows the temperature comfort bounds, i.e. the occupancy is set to 1 when the comfort bounds are tight and 0 when they are loose. The solver and software tools used to solve the optimization problem are as explained in Section 2.4. Lastly, the parameters of the building presented in Section 3 are shown in Table 2.
For what concerns the SBMPC controllers, we choose 4 different number of scenarios: 10, 20, 30, and 40. Moreover, we perform the same simulations varying the parameter α in Section 3, choosing the values in the set {50, 100, 200, 500}; recall that a higher α means a higher focus on the comfort of the occupants rather than on the economical cost. We indicate respectively by SBMPC-Mod-n s and by SBMPC-Lin-n s , n s ∈ {10, 20, 30, 40}, the SBMPC controller with the Modelica model and the SBMPC controller with the linear model, considering n s scenarios.

Results and discussion
We focus our attention on four different aspects: 1. an analysis for the performance of the controller with different values of α. 2. A comparison between the nonlinear Modelica model and the linear model.

A comparison between SBMPC strategies and
DetMPC strategies. 4. A comparison between the SBMPC strategies with different numbers of scenarios.
Lastly, we pick a single representative optimization result and discuss it in more detail. We show the results of the simulations in Figures 4, 6-9 and Tables 3-5. In Table 3, we show the total closed-loop costs for each strategy and each different α. In Table 4, we show the percentage of each subcost with respect to the total closed-loop cost for each strategy and each different α, i.e. we show Je Je+αJ d and αJ d Je+αJ d . In Table 5 we show the total amount of discomfort using the unit measure K · h, as standard in the literature [13,15,17,25], i.e. we show the integral of the comfort bounds violation. Figure 4 reports the results presented in Table 4 in a graphic way; the same holds for Figure 6 and Table 5  9 show respectively the temperature evolution, the heating power, and the temperature and irradiance profiles for the representative simulation.

Performance with different values of α
In this section, we analyze how the different values of α alters the performance of the controllers in terms of energy cost and discomfort. This analysis should always be carried out when considering a case study as the one presented here, in order to understand which is the range of values for α that provide the best trade-off between comfort and energy cost reduction.
From Tables 3-5 we can notice that the larger the α, the larger, in general, the total costs and the lower the discomfort. This is as expected, since the role of α is to penalize the discomfort and a larger value means that we aim for a lower discomfort. We have noticed through simulations that, for this specific case study, a value of α lower than 50 yields a very high and unacceptable discomfort cost, while for values larger than 500, the energy costs increase highly without yielding a high reduction in the discomfort costs. Therefore, we focus our analysis on α in the range [50,500].
We can observe that for α = 50 the discomfort cost is high and that the comfort could be improved by increasing α. For α ≥ 100, the discomfort reaches more acceptable levels and this happens by consuming a larger quantity of energy in heating and thus increasing the total costs. However, the small decrease in the discomfort between the case α = 500 and α ∈ {100, 200} does not seem to justify the large increase in the total cost observed for α = 500. Therefore, we can claim that, for this case study, the optimal values for α are in the range [100,200]. Note also that from Figure 4 and Table 4 we can notice how the two costs, i.e. energy and comfort, compose the total cost. As α increases, we notice an increase in the αJ d cost with respect to the J e , as a higher α penalizes more even the small deviations from the comfort bounds. Therefore, we should not wrongly conclude that less energy is consumed for higher values of α, but rather that small deviations are penalized more.

Comparison between the nonlinear Modelica model and the linear model
Recalling Remarks 4-5, we present here a comparison between the controllers that use the Modelica model with respect to the ones that use the linear model. While using the Modelica model in the closed-loop for computing the evolution of the system results in a bias towards the nonlinear controllers, it is nevertheless of interest analyzing how good the linear controllers are compared to the ones using the Modelica model.
From Table 3, we can see that all the controllers that use Modelica perform better than their linear counterparts for all the values of α ≥ 100. For α = 50, instead, the linear model yields a lower total cost than the one of Modelica in 3 out of 5 cases. Nevertheless, the linear model might seem to work better, i.e. to have a lower total cost, for a lower value of α because it always allows a large discomfort cost and it cannot manage well to keep the temperature within or close to the comfort bounds. The total cost can therefore be lower than for the Modelica-based controllers, but this occurs because the energy cost is low and the discomfort cost, although high, does not have a large impact on the total cost for a small α. This also explains why for a large α, i.e. for α ≥ 100, the total cost of the linear model controllers can become much higher than the one of the controllers that use Modelica. Indeed, the discomfort cost is always high, but the larger penalization, i.e. the larger α, makes the total cost much higher. This fact can also be observed from Table 5, where we observe a much higher discomfort for all the controllers with the linear models. The same conclusions can be drawn by analyzing Figure  6, which displays the results of Table 5 Table 4.
tions. This concept will be discussed also in Section 5.2.5. We can therefore conclude that, independently of the bias mentioned in Remark 4, the linear model in this case fails to provide satisfactory performance in terms of comfort for the occupants of the building.

Comparison between SBMPC strategies and
DetMPC strategies By checking again Table 3 we can compare the SBMPC strategies to the DetMPC ones. It can be noted from the table that, for all the values of α, the SBMPC con-trollers perform almost always better than their deterministic counterpart, both for the linear and the nonlinear Modelica model. Note that the reduction, although not very large, is still consistent and it ranges from a minimum of 1.98% to a maximum of 14.03% for the controllers that use the Modelica model. Furthermore, by checking Table  5 and Figure 6, we can notice that the SBMPC strategies perform better than the DetMPC ones also in terms of comfort. Therefore, the SBMPC controllers can improve the overall performance, both in terms of total costs and discomfort, with respect to their DetMPC counterparts.   By analyzing the results of Table 3, there does not seem to be a value for the number of scenarios that outperforms the other values, i.e. the performance does not seem to increase by increasing the number of scenarios. In 2 out of 4 columns of Table 3, the SBMPC-Mod-20 achieves the best performance among the SBMPC-Mod controllers and in the other two cases, a number of scenarios equal to respectively 10 and 40 appear to be better than the other values. Therefore, increasing the number of scenarios does not seem to directly lead to a decrease in the total cost. This could be related to the fact that, while increasing the number of scenarios makes the system more robust to disturbances, it also makes the problem more complex to solve. Therefore, it might happen that, the larger the number of scenarios, more local minima exists, and the more likely it is that the solver converges to a suboptimal local minimum. Note that this issue does not affect the controllers with the linear model, as the problem solved in that case is a quadratic programming one; thus it is convex, and does not suffer from local minima issues.

Representative simulation
In this section, we present a representative simulation, i.e. one week of simulation of the building, by showing the temperature evolution, heating power, and external disturbances for a specific value of α and of the number of scenarios. While the analysis of the results shown in this section are related to a specific case, the results can be generalized and the analysis of a case with a different α would be similar to what is here presented.
We compare three different control strategies here, namely SBMPC-Mod, SBMPC-Lin and PIMPC. We show in Figure 7 the temperature evolution inside the room, in Figure 8 the heating power, and in Figure 9 the temperature and irradiance profile, for one week of simulation with 20 scenarios and α = 100. For Figure 7, we also show the lower comfort bounds.
By analyzing the Figures 7-8, we can note that the PIMPC manages to keep the temperature within the comfort bounds by using properly the heating power, thanks to the knowledge of the actual values of the future disturbances. For what concerns SBMPC-Mod and SBMPC-Lin, we can notice in Figure 7 what we have already under- lined in Section 5.2.2, i.e. the fact that a controller that uses a linear model is not able to keep the temperature within the comfort bounds. We see indeed that, for most of the time, SBMPC-Lin yields the temperature profile that has the lowest value. This can also be observed from Figure 8, where we can notice that SBMPC-Lin heats less than SBMPC-Mod and it also starts heating later. On the other hand, SBMPC-Mod is able to maintain a larger temperature in the room and to be closer to the temperature comfort bounds. It uses more heating power, as can be seen from Figure 8, which leads to a higher energy cost compared to SBMPC-Lin, but also leads to an overall lower total cost and higher comfort from the user, as can be observed from Tables 3-5.
In the first two days of simulation, which correspond to the weekend days, the heating power is turned off due to the low comfort bounds and to the prediction horizon of 24 h. However, it is possible to observe an increase in the room temperature, both for the day 16/12/2017 and 17/12/2017; for the latter, the rise in the temperature is even steeper. This can be explained by looking at Figure  9, where we can observe that the day 17/02/2017 was a particularly sunny day. Hence, the steep increase in the temperature is due to large value of the solar irradiance in that particular day. This shows how important the influence of the external disturbances on the building can be, which once again underlines the importance of having good forecasts and it corroborates our choice of a scenariobased approach.
Note also that, while a longer prediction horizon could have been more beneficial and could have led to a higher performance of the controllers, we chose to focus our anal-ysis on a horizon of 24 h as it provides a good trade-off between performance and computation time of all the simulations. The analysis performed here is still valid even when considering other prediction horizons. Moreover, carrying out a study on which prediction horizon is more beneficial for the considered building is beyond the scope of this paper.

Summary
We can summarize the observations obtained from the results of the simulations: • too low values of α, i.e. α < 100, yield a high discomfort and too high values of α, i.e. α ≥ 500 yield a very large total cost without leading to a large improvement of the comfort. A trade-off between the two costs seems to be well achieved by a value of α between these two extrema, i.e. α ∈ [100, 200]. Moreover, tuning the parameter α is of great importance in order not to obtain a very large discomfort for the occupants of the building.
• For all the values of α, the linear model shows a very large value of discomfort. In general, the linear model fails to capture many model dynamics, by heating much less and thus leading to temperatures that are well outside the comfort bounds.
• SBMPC can improve the performance with respect to DetMPC strategies, for both the linear and the Modelica models. This is due to the fact that SBMPC strategies consider different external disturbances scenarios and they have a large impact on the temperature evolution of the room.
• Increasing the number of scenarios does not seem to lead to a large decrease in the cost. This can be related to the fact that increasing the number of scenarios also increases the optimization complexity.
Remark 10. Note that all the simulations performed in this work refer to a single, specific building. However, the controller designed can be applied to any building, as long as the model is adapted to the specific building under control. For instance, a similar controller, using a Modelica model and a deterministic MPC algorithm, has been applied to a different building with successful results in [13]. Moreover, the scenario generation method presented here can also be applied to several disturbances that affect buildings. Therefore, while the results presented in this section refer to a single building, the method can be applied to any building. The results presented in Tables 3-5 will change if another building is considered, but we expect in any case similar results when applying the presented method to other buildings.

Conclusions
We have presented a stochastic SBMPC controller using a Modelica nonlinear model that can be applied to building heating in buildings and that overcomes the limitations of both deterministic and linear MPC approaches. The building under control is affected by several external disturbances, e.g. outside temperature, solar irradiance, and we have proposed a new approach for generating disturbance scenarios that, unlike the existing methods from the literature, satisfies all the important properties of scenario generation methods for time series data. This proposed scenario generation method can be used in the SBMPC controllers.
To analyze and study the control approach, we considered a real building and performed several simulations to compare the controller that uses the linearized model against the controller that uses the Modelica model, different cost weights in the MPC cost function, deterministic MPC against SBMPC, and lastly different number of scenarios. Based on the results, we showed that SBMPC controllers outperform the deterministic MPC controllers both with and without the Modelica model. At the same time, the linear model has shown not to be able to capture many model dynamics and this leads to poor performance, justifying the usage of more advanced models, e.g. Modelica ones.
As future work, we will develop a model with e.g. En-ergyPlus to be used in the closed-loop simulations to improve the comparison performed between the linear and nonlinear model in Section 5. After such step, we will also perform experiments in the real building. Moreover, a quantitative study on different scenario generation methods can be considered, in order to assess how the performance of the controller varies with the different scenario generation methods, as well as other stochastic MPC algorithms. Furthermore, a thorough comparison could be performed between our proposed method and other ones present in the literature of building heating systems, e.g. H ∞ , fuzzy control, rule-based control. On top of that, a distributed or decentralized MPC controller can be developed to control independently each room, which might be more beneficial for very large buildings compared to controlling all the rooms with a single centralized controller. Lastly, it is known that occupancy can largely affect the energy performance of buildings. Occupancy data was not available in this study, but we suggest to carry out a characterization study of such disturbance, so that occupancy scenarios can be generated and the overall performance of the proposed SBMPC controllers can be further assessed.