**2. Methodology**

In this study, the goal of the MPC was to minimize the total energy cost for the building thermal demand satisfaction. The typical structure of an MPC is shown in Figure 1. It is mainly composed of two parts: the building predictive model and the optimizer. The building predictive model should be able to dynamically forecast the building's energy response in a certain period (prediction horizon, *ph*), while its inputs can vary both in a controlled (manipulated variables) and in an uncontrolled (disturbances) way. To solve the optimization problem, it is important to define a proper objective function and to respect the system constraints; in this way, the optimizer has the possibility to select the best control actions to maximize the performance. Following a "receding horizon" logic, the MPC updates the best control action at each timestep, moving the prediction horizon forward and repeating the optimization [5].

**Figure 1.** Architecture of a typical building's model predictive control (MPC).

As anticipated, two different approaches for the building prediction model are used in this work: a data-driven model based on an ANN and a physical-based model with a thermal–electrical analogy structure. In order to obtain training data for the ANN and to check the MPC effectiveness, a detailed building model designed in TRNSYS [14] was adopted as a control building.

The MPC routine was written in MATLAB [32], and for each time step the controlled building started the MATLAB engine to run the controller. The uncontrolled inputs of the models were weather conditions and heat gains, while the manipulated variable was the hourly building energy demand. Since the objective of the work was to focus on the comparison between the two approaches for building modeling in an operative MPC, an ideal HVAC system was considered in the building, i.e., the thermal energy demand was treated as a control action (Figure 1).

The optimizer solved the optimization problem in the prediction horizon, *ph*, in order to minimize the total energy cost with constraints on the internal comfort conditions. As an incentive for the exploitation of the energy flexibility, a dynamic energy cost tariff was considered [33]. In addition, to amplify the cost variations, a penalty signal was used in the MPC optimizer. It was obtained with a statistical method based on mean and standard deviation:

$$p(t) = \frac{c(t) - \mu\_{\mathcal{C}}}{\sigma\_{\mathcal{C}}} \tag{1}$$

where *p* is the penalty signal at each time *t*, *c* is the energy cost, while μc and σc are the cost signal mean value and standard deviation, respectively. Additional details about the penalty signal will be discussed in Section 3.

The performance of each model was assessed by evaluating the root mean square error, *RMSE*, which is defined as:

$$RMSE = \sqrt{\frac{1}{n} \sum\_{j=1}^{n} \left( y\_{\text{model},j} - y\_{\text{data},j} \right)^2} \tag{2}$$

where *y* is the variable being evaluated and *n* is the number of points considered. Another index that will be used in the results is the root square error, *RSE*, defined as:

$$RSE\_i = \left| y\_{\text{model},j} - y\_{\text{data},j} \right| \tag{3}$$

The first analysis consisted of evaluating the deviation of the models with respect to the building reference data. Then, in order to assess the accuracy of the models in the operative MPC, the *RMSE* was calculated between the prediction of the MPC building model and the actual thermal behavior of the building. Due to the different mathematical formulations of the two models, the optimization problem was different for the two cases. The following Sections 2.1 and 2.2 discuss the mathematical formulations of the two approaches.

#### *2.1. MPC with Physical-Based System Model*

A lumped-parameter model based on the thermal–electrical analogy was chosen as the building prediction model in the physical-based MPC. The building thermal dynamics is represented by an equivalent circuit of thermal resistances and capacitances [34]. A third order model was selected since it represented a good compromise between network complexity and capability of predicting the short-term dynamics of the building [35].

As shown in Figure 2, three thermal nodes were used. Each node is described by a temperature (*T*) and a thermal capacitance (C). Then, four thermal resistances (R) were used to model the heat transfer between the nodes.

**Figure 2.** Thermal network for building model.

All the numerical values of the parameters (R and C) were deducted by the knowledge of the thermal and geometrical building features. Specifically, the first node (*T*e, Ce) represents the external building thermal mass, the second node (*T*air, Cair) is the indoor air node, while the last node (*T*i, Ci) is the internal building thermal mass. As suggested by EN ISO 13790 [36], Ce and Ci are calculated by summing the heat capacitances of all the building elements (up to the thermal insulation) in direct thermal contact with the internal air zone.

As concerns the thermal resistance, Rw is the thermal resistance from the indoor air node temperature to the ambient air temperature (*T*amb), due to air changes and windows. Ree and Rie are the thermal resistances between the external building thermal mass node and *T*amb and *T*air, respectively. They are calculated as equivalent thermal resistances due to the conductive heat transfer of all the building envelope layers, from outdoor to the thermal insulation for Ree, and from thermal insulation to indoor for Rie. These thermal resistances also take into account the convective heat transfer phenomena between the external surface and ambient temperature (Ree) and between the internal building envelope surface and indoor air temperature (Rie). In the same fashion, Ri considers the thermal resistance between the indoor air node and the internal thermal mass *T*i. The heat fluxes, directly applied to the internal thermal nodes *T*air and *T*i, are the cooling power derived by an ideal HVAC system (*Q*) and the total heat gains (*G*). The latter includes both solar and internal contributions, which are provided with a scalar factor (f) for both the internal air and the internal mass node.

The dynamics of the resistance–capacitance (RC) model can be represented by the following equations:

$$\mathbf{C\_{air}}\frac{\mathbf{d}T\_{\rm air}}{\mathbf{d}t} = \frac{\left(T\_{\rm e} - T\_{\rm air}\right)}{\mathbf{R\_{ice}}} + \frac{\left(T\_{\rm amb} - T\_{\rm air}\right)}{\mathbf{R\_{W}}} + \frac{\left(T\_{\rm i} - T\_{\rm air}\right)}{\mathbf{R\_{i}}} + \mathbf{f\_{air}}\mathbf{G} + \mathbf{Q} \tag{4}$$

$$\text{C}\_{\text{e}} \frac{\text{d}T\_{\text{e}}}{\text{d}t} = \frac{(T\_{\text{amb}} - T\_{\text{e}})}{R\_{\text{ee}}} + \frac{(T\_{\text{air}} - T\_{\text{e}})}{R\_{\text{ice}}} \tag{5}$$

$$\mathbf{C}\_{\rm i} \frac{\mathbf{d}T\_{\rm i}}{\mathbf{d}t} = \frac{(T\_{\rm air} - T\_{\rm i})}{\mathbf{R}\_{\rm i}} + \mathbf{f}\_{\rm i}\mathbf{G} \tag{6}$$

Using these relations, a discrete time invariant state–space formulation can be set up:

$$\mathbb{E}\left[X(k+1)\right] = \left[\mathbf{A}\right]\left[X(k)\right] + \left[\mathbf{B}\right]\left[\mathbf{U}(k)\right] \tag{7}$$

where [*X*] = [*T*air *T*e *T*i]<sup>T</sup> represents the system state at each discrete time *k*, [*U*] = [*T*amb *G Q*]<sup>T</sup> is the vector of the inputs, and [A] and [B] are coefficient matrices. A discrete time (*k*) of 1 hour is adopted, as well as for the simulation time step (*t*).

Along with the building prediction model, the MPC must include an optimizer. In the present work, the objective of the MPC was to minimize the total energy cost for the building thermal demand satisfaction while the indoor temperature remains in an allowed comfort range. For the physical-based MPC, a simple linear programming (LP) problem was formulated. In this optimization problem, both the objective function and constraints must be linear functions of the decision variables. The objective function (OF) is calculated as the sum of the building thermal energy consumption (*Qk*) multiplied by the penalty signal *pk* at each time step *k* in the whole prediction horizon (*ph*):

$$\text{OF} = \sum\_{k=1}^{pl} p\_k Q\_k \tag{8}$$

Therefore, the minimization problem can be written as:

$$\text{minOF} \tag{9}$$

subject to the following comfort and ideal HVAC constraints:

$$\forall \; k = 1, \dots, \; ph \qquad T\_{\min} \le T\_{\text{air},k} \le T\_{\text{max}} \tag{10}$$

$$\forall \; k = 1, \dots, \text{ } \text{ } \text{ } \text{ } \text{ } \quad \text{ } 0 \le Q\_k \le Q\_{\text{max}} \tag{11}$$

The LP optimization problem is solved at each time step within a MATLAB script, according to a "dual-simplex" algorithm. The actual air zone temperature, *T*air(*t*), is passed to the MPC as the starting condition for the optimization. Based on the receding horizon principle, the control action at the controlled building time, *t*, is the first value of *Qk* in the optimal decision variable sequence. At the following time step, *t*+1, the optimization problem is re-solved, moving forward the dataset of the disturbances by a time step. Figure 3 shows the scheme of the operation.

**Figure 3.** MPC receding horizon scheme.

#### *2.2. MPC with Data-Driven System Model*

The data-driven system used in the present study is based on an artificial neural network (ANN), a mathematical model that reflects the functioning of a biological brain [37]. In an ANN, the inputs (*x*) have the same role of biological dendrites, while the outputs (*y*) can be regarded as biological axons. The processing of the data takes place in the neurons, which, in an ANN, apply a nonlinear activation function, g, on the input data. Being a pure mathematical model without physical meaning, an ANN needs to be trained with existing data. The purpose of the training, which can be carried out with di fferent error minimization techniques, is to determine the coe fficient weights and biases of the network, i.e., the parameters that fully describe an ANN. Taking into account a feedforward ANN consisting of only one layer of neurons (also called a hidden layer because its activation values are not directly accessible from outside the network) and a linear activation function in the output layer (with just one output), the mapping carried out by the ANN on the input can be expressed as follows:

$$y = \sum\_{j=1}^{m} \left( \left\| \left( \sum\_{i=1}^{d} \mathbf{w}\_{ji} \mathbf{x}\_{i} + \mathbf{b} \right) + \hat{\mathbf{b}} \right\| \right) \tag{12}$$

where d is the number of inputs, <sup>w</sup>*ji* is the weights matrix of the inputs, b is the bias vector of the inputs, m is the number of neurons in the hidden layer, *w*ˆ *j* is the weights matrix of the hidden layer, and b is the bias vector of the hidden layer. ˆ

Generally, training data are divided into inputs and targets, and a well-trained ANN is expected to determine its outputs with a low deviation in respect to the provided targets. The available data of the system under study need to be examined carefully in order to train the ANN only with the inputs that most influence the objective targets. If the physics of the system under study is complex, a method to individuate the most e ffective inputs involves the use of statistical approaches such as factor analysis. Instead, if the physics of the system is not entirely unknown, the operator can try to select the input variables that mostly influence the desired target. In the present work, since the thermal behavior of the building is known, the four input variables (outdoor temperature, solar gains, internal gains, and indoor temperature) that most significantly influence the target variable (the thermal power required by the building) were selected.

The ANN was trained with 168 hourly-based input/target data referred to as a typical week. These data, provided by the building simulation environment designed in TRNSYS, were not obtained with a fixed set-point condition. In this case, in fact, the ANN's performance would have been insignificant, as there would have been no correlation between the output and the controlled variable (the indoor temperature). Indeed, the control unlocks the energy flexibility of the building and allows indoor temperature variations in a given comfort range, as explained previously. Therefore, to improve the ANN's prediction capability, the building simulation environment was allowed to work with multiple indoor set-point temperatures, varying within a reasonable comfort range. To avoid the overfitting of the data, i.e., an exaggerated interpolating behavior of the ANN, only a fraction of the dataset was actually used to train the network. Specifically, the initial dataset of 168 points was randomly divided into three subsets: a training set (60%), a validation set (20%), and a test set (20%). While the training data were used by the ANN to complete its training, the validation set was used as an internal interrupt criterium to end the training if overfitting occurred. The test set, instead, was used to evaluate the ANN's performance after training, in order to check its prediction capability with new data.

ANNs are available in di fferent architectures, based on the physical–mathematical problem that is being studied. In the present work, since the goal was to estimate the thermal power required by a building, a fitting ANN was chosen. In MATLAB, fitting ANNs have a feedforward architecture and are trained according to a Levenberg–Marquardt backpropagation algorithm, which uses regression analysis and *RMSE* to evaluate the performance. As it is well-known that even one layer of neurons is su fficient to represent complex, nonlinear problems [37], in this study, one hidden layer with five neurons was used. The neurons used a hyperbolic tangent sigmoid as activation function. The ANN therefore had the architecture as represented in Figure 4.

**Figure 4.** Artificial neural network (ANN) architecture of building prediction model.

After training, the ANN model was used in the MPC to predict the building's thermal demand for a given prediction horizon, *ph*. To this purpose, the ANN inputs were divided into uncontrolled (outdoor temperature, solar gains, and internal gains) and controlled (indoor temperature) variables (Figure 1). In this way, it was possible to manipulate the ANN as a function of the indoor temperature and to carry out an optimization process in order to minimize the overall energy cost in the time period *ph*. The objective function of the optimization algorithm can be therefore written as in Equation (8), subject to the constraints defined in Equations (10) and (11). Since the ANN function for the thermal demand was nonlinear, the optimization algorithm chosen in MATLAB was a programming solver based on the gradient method that uses an initial value for the indoor temperature as first attempt of solution. In the same fashion as the LP optimization problem defined for the physical-based model, the control action of the ANN-based MPC works according to the receding horizon principle (Figure 3).
