**4. Practical implementation of the Bayesian predictor**

According to what reported in sub-section 4.2, the conditional probability tables of Bayesian Networks (BNs) might be defined from a combination of a priori knowledge about the involved physical processes and a dataset of collected measures. Part of the dataset is usually used for probability estimation, while the remaining part is usually used for validating the BN. Such validation may be done according to the following steps:


Once each instantaneous index is computed, a global performance index for the whole validation dataset might be computed for each node.

## **4.1. The metrics used for validation**

The performance of a prediction model may be evaluated by using different kinds of indices. As well known, referring to prediction of a variable at a certain time *i*, the difference between the predicted value *X* ^ *i* and the actual value *Xi* is defined as error *Ei* ≜ *X* ^ *<sup>i</sup>* - *Xi* . Since the validation of a prediction model is interested just in the magnitude of the error, either its absolute value called absolute error *AEi* ≜|*Ei* | or its squared error *SEi* ≜*Ei* 2 can be used. Percentage error is here defined as the ratio (%) between the error and the actual value of the variable: *PEi* ≜100⋅*Ei* / *Xi* . In order to have a global performance index to be evaluated over the whole validation dataset made up of *K* samples, these instantaneous indices must be combined into global indices. Therefore, the mean absolute error is defined as:

$$\text{MAE} \triangleq \frac{1}{K} \sum\_{1}^{K} \mid \stackrel{\wedge}{X}\_{i} \text{ - } X\_{i} \,\vert \,\tag{12}$$

and the root mean square error is:

Thanks to this algorithm, which is easily implemented by computer programs, the distribution re-adjust its shape according to data provided by available sample data. It exploits the notion of experience, that is quantitative memory which can be based both on quantitative expert judgment and past cases. The prior parameters in the first version of the network are set with a particular equivalent sample size based on expert knowledge (by tuning the values *N'*ijk). Then subsequent data are added, starting from the given equivalent sample size, so that the two Dirichlet parameters reshape the probability density function according to observations included in the dataset, which determine the values of the set of parameters *N*ijk. Furthermore,

*=k* and *Π*<sup>i</sup>

*i*=1 *n* ∏ *j*=1 *qi Nijk* ' + *Nijk Nij* '+ *Nij*

According to what reported in sub-section 4.2, the conditional probability tables of Bayesian Networks (BNs) might be defined from a combination of a priori knowledge about the involved physical processes and a dataset of collected measures. Part of the dataset is usually used for probability estimation, while the remaining part is usually used for validating the BN.

**3.** evaluation of estimation/prediction performance by computing a proper performance

Once each instantaneous index is computed, a global performance index for the whole

The performance of a prediction model may be evaluated by using different kinds of indices. As well known, referring to prediction of a variable at a certain time *i*, the difference between

validation of a prediction model is interested just in the magnitude of the error, either its

Percentage error is here defined as the ratio (%) between the error and the actual value of the

is defined as error *Ei* ≜ *X*


. In order to have a global performance index to be evaluated over

^ *<sup>i</sup>* - *Xi*

2

. Since the

can be used.

*=j* in the next case *C*m+1 to be seen in the

(11)

it can be shown that the probability that *X*<sup>i</sup>

14 Dynamic Programming and Bayesian Inference, Concepts and Applications

, and *Nij* = ∑

*k*=1 *ri Nijk* .

**4. Practical implementation of the Bayesian predictor**

Such validation may be done according to the following steps:

**2.** evidence propagation through the Bayesian reasoning engine;

and the actual value *Xi*

**1.** instantiation of the BN with collected evidences;

validation dataset might be computed for each node.

**4.1. The metrics used for validation**

^ *i*

absolute value called absolute error *AEi* ≜|*Ei*

where *Nij*

'

index.

the predicted value *X*

variable: *PEi* ≜100⋅*Ei* / *Xi*

= ∑ *k*=1 *ri Nijk* '

database (i.e. expected value) is computed by means of the equation:

*P*(*Cm*+1 | *D*, *Bs*, *ξ*) =∏

$$\text{RMSE} \triangleq \sqrt{\frac{1}{K} \sum\_{i}^{K} \left(\widehat{\hat{\mathbf{X}}\_{i}} - \mathbf{X}\_{i}\right)^{2}}.\tag{13}$$

The indices of the predicted variables are related to different physical quantities with different units. For this reason, these indices must be normalized with respect to their typical range of variation. Two more indices are defined and practically used for evaluating prediction models:

$$\text{NMMAE} \triangleq \frac{\frac{1}{K} \Sigma\_1^K \parallel \stackrel{\wedge}{X}\_i \cdot X\_i \parallel}{\lfloor X\_{\text{max}} \cdot X\_{\text{min}} \rceil} \tag{14}$$

$$\text{INRMSE} \triangleq \frac{\sqrt{\frac{1}{K} \sum\_{i}^{K} (\hat{X}\_{i} \cdot X\_{i})^{2}}}{\lfloor X\_{\text{max}} \cdot X\_{\text{min}} \rfloor} . \tag{15}$$

ASHRAE Guideline 14-2002 [24] establishes that for calibrated simulations, the *CVRMSE* and *NMBE* of energy models shall be determined for each calibration parameter by comparing simulation-predicted data to the utility data used for calibration. The proposed indices are the coefficients of variation of the root mean square error (*CVRMSE*) and normalized mean bias error (*NMBE*). All these indices are normalized with respect to the arithmetic mean of the variable. Following this guideline, the *RMSE* has been selected as the main performance index for evaluating the accuracy of a BN. However, the mean value of a variable is a good normal‐ ization factor only if the variable is always positive (or always negative). For this reason, the range of the considered variable has been taken as a normalization factor and the *NRMSE* has been selected as final index for the design process of the BN: it includes information about both bias and variance of the error.

### **4.2. Instantiation of BNs for MPC**

The control block in the framework depicted in Fig. 1 needs to evaluate a cost function for each candidate control policy. However, since discrete nodes are often used in BNs for describing nonlinearities, the resulting cost function may also present significant nonlinearities. One of the main problems in MPC consists in guaranteeing closed-loop stability for the controlled system. As stated in [25], stability about MPC requires some continuity assumptions. The usual approach to ensure stability is to consider the value function of the MPC cost as a candidate Lyapunov function. Even if closed-loop stability is guaranteed in some way without the continuity assumption, the absence of a continuous Lyapunov function may result in a closedloop system that has no robustness. This implies that, whenever possible, a continuous cost function is highly recommended in MPC: it implies better stability and robustness and faster convergence of the optimization algorithm. Moreover, from a modelling point of view, the real buildings are usually continuous systems, or at least hybrid systems (i.e. mixed discretecontinuous) but never pure discrete: a discrete approximation would produce a big variance of the prediction error.

The cumulative distribution function (cdf) is then given by *Φ*(*x*)=1 / 2*π∫*

<sup>1</sup> <sup>+</sup> *<sup>e</sup>* <sup>2</sup>*<sup>z</sup>* ' , where *<sup>z</sup>* '

<sup>1</sup> <sup>+</sup> *<sup>e</sup>* <sup>2</sup>*<sup>z</sup>* , where *<sup>z</sup>* =0.7988( *<sup>x</sup>* - *<sup>μ</sup>*

\_

*σ* =max (*σ*

\_ , *σmin*

For each interval of a standard discrete node, probability is given by *qi* <sup>=</sup> *qi*

=Φ(*xi*+1) - Φ(*xi*

For a discrete periodic variable (like a direction) defined in the interval *x*0; *xL* (that for a direction is 0;2*π* ), the pdf can be split into periodic intervals shifted each other by integer

If standard deviation is sufficiently smaller than the period *σ*≪(*xL* - *x*0), this infinite series can be approximated by the three central intervals, that is for *k* = - 1,0, 1. Moreover, when *σ* ≅(*xL* - *x*0), probability is almost equally distributed all over the interval and the approxi‐ mation introduced by neglecting the stubs for *k* < - 1, *k* >1 can be simply compensated by renormalizing the resulting distribution. Thus, in this case, probability of each interval is given as in the previous algorithm but the right and left stubs of the Gaussian are overlapped in order to consider periodicity of the discrete variable. For each interval of a periodic discrete node,

*qi* '

'

' where:

∑*i*=0 *<sup>L</sup>* -1 *pi*

*Φ* '

<sup>Φ</sup>(*x*)<sup>≈</sup> *<sup>e</sup>* <sup>2</sup>*<sup>z</sup>*

*<sup>x</sup>* must be replaced with *<sup>x</sup>* - *<sup>μ</sup>*

Given any standard deviation *σ*

\_ falls:

; *xl*+1 in which *μ*

multiples of period (*xL* - *x*0):

probability is given by *pi* <sup>=</sup> *pi*

*xl*

(*x*)<sup>≈</sup> *<sup>e</sup>* <sup>2</sup>*<sup>z</sup>* '

*<sup>σ</sup>* :

close form. An effective and sufficiently accurate way for approximating it is provided by [26]:

In order to consider a general distribution with mean *μ* and standard deviation *σ*, the variable

>0, and mean value *μ*

standard deviation must be at least greater than the standard deviation *σmin* of the interval

), *<sup>σ</sup>min* <sup>=</sup> *xl* +1 - *xl*

*<sup>σ</sup>* )(1 <sup>+</sup> 0.04417( *<sup>x</sup>* - *<sup>μ</sup>*

Bayesian Networks for Supporting Model Based Predictive Control of Smart Buildings

*x*<sup>0</sup> + *k*\*(*xL* - *x*0); *xL* + *k*\*(*xL* - *x*0) , *k* ∈Z (20)

\_


http://dx.doi.org/10.5772/58470

=0.7988*x*(1 + 0.04417*x* 2) (16)

*<sup>σ</sup>* )<sup>2</sup>

*dt* that has no

17

) . (17)

>0 to be approximated, the value of

<sup>12</sup> . (18)

∑*i*=0 *<sup>N</sup>* -1 *qi*

) (19)

'

' where:

Due to these reasons, the BNs must be continuous with respect to continuous variables. When BNs contain discrete interval nodes (this is often the case in practice, since it allows to model more complex relationships between nodes), an instantiation as "*hard evidence*" (i.e. one state only of a node is assigned 100% likelihood while the other states are given 0%) produces an abrupt variation of the outputs only when passing from one interval to the next one, thus showing a discontinuous behaviour. Therefore, an instantiation must be always done as "*soft evidence*" (i.e. more than one state of a node is instantiated, in fact shaping a probability density function). The idea, with reference to Fig. 5, is to approximate each evidence with nominal value *μ* and estimated standard deviation *σ* as a Gaussian distribution *N* (*μ*;*σ*). This distribu‐ tion is then approximated by an interval discrete distribution that is actually instantiated into the BN as "soft evidence". When the evidence has been propagated, the beliefs retrieved by the network are used to iterate the prediction for an arbitrary number *P* of steps (see Fig. 7). The final prediction sequence is then summarized as a sequence of normal continuous distributions *N* (*μi* ;*σi* ), *i* =1, …, *P*.

Technically, the key challenge is approximating a continuous probability density function by means of a discrete one, as reported in the following. Given mean value and standard deviation of a normal distribution *N*( *μ* \_ *i*;*σ* \_ *i*), and an interval discrete random variable de‐ fined by the values *x*0, *x*1, …, *xL* (that corresponds to the *L* + 1 bounds of the *L* inter‐ vals), the probabilities of each interval must be found such that mean value and standard deviation of the discrete distribution are as close as possible to *μ* \_ and *σ* \_ respectively. As well known, a normal Gaussian probability density function (pdf) takes the form *f* (*x*)=1 / 2*πe*-*<sup>x</sup>* <sup>2</sup> /2 , - *∞* < *x* <*∞*.

**Figure 5.** Instantiation of a measure as "soft evidence" and corresponding information flow. Instantiation and propa‐ gation is repeated in the prediction loop (dashed), as described in Fig. 7.

The cumulative distribution function (cdf) is then given by *Φ*(*x*)=1 / 2*π∫* -*∞ <sup>x</sup> e*-*<sup>t</sup>* <sup>2</sup> /2 *dt* that has no close form. An effective and sufficiently accurate way for approximating it is provided by [26]:

$$\mathbf{O}^{'}(\mathbf{x}) \approx \frac{e^{\frac{\mathbf{x}^{\cdot 2^{\cdot 2^{\cdot 2^{\cdot 2^{\cdot 2^{\cdot 2^{\cdot 2^{\cdot 2^{\cdot 2^{\cdot 2^{\cdot 2^{\cdot 2^{\cdot 2}}}}}}}}}}}}}}}}}}} \ \ \ \ \text{where } \text{z} \text{ } = 0.7988 \text{x} \{1 + 0.04417 \text{x} } ^ { \text{2} } \text{} \}} \tag{16}$$

In order to consider a general distribution with mean *μ* and standard deviation *σ*, the variable *<sup>x</sup>* must be replaced with *<sup>x</sup>* - *<sup>μ</sup> <sup>σ</sup>* :

$$\Phi(\mathbf{x}) \approx \frac{\epsilon^{2z}}{1 + \epsilon^{2z}}, \quad \text{where } z = 0.7988 \left(\frac{\mathbf{x} \cdot \boldsymbol{\mu}}{\sigma}\right) \left(1 + 0.04417 \left(\frac{\mathbf{x} \cdot \boldsymbol{\mu}}{\sigma}\right)^2\right). \tag{17}$$

Given any standard deviation *σ* \_ >0, and mean value *μ* \_ >0 to be approximated, the value of standard deviation must be at least greater than the standard deviation *σmin* of the interval

*xl* ; *xl*+1 in which *μ* \_ falls:

function is highly recommended in MPC: it implies better stability and robustness and faster convergence of the optimization algorithm. Moreover, from a modelling point of view, the real buildings are usually continuous systems, or at least hybrid systems (i.e. mixed discretecontinuous) but never pure discrete: a discrete approximation would produce a big variance

Due to these reasons, the BNs must be continuous with respect to continuous variables. When BNs contain discrete interval nodes (this is often the case in practice, since it allows to model more complex relationships between nodes), an instantiation as "*hard evidence*" (i.e. one state only of a node is assigned 100% likelihood while the other states are given 0%) produces an abrupt variation of the outputs only when passing from one interval to the next one, thus showing a discontinuous behaviour. Therefore, an instantiation must be always done as "*soft evidence*" (i.e. more than one state of a node is instantiated, in fact shaping a probability density function). The idea, with reference to Fig. 5, is to approximate each evidence with nominal value *μ* and estimated standard deviation *σ* as a Gaussian distribution *N* (*μ*;*σ*). This distribu‐ tion is then approximated by an interval discrete distribution that is actually instantiated into the BN as "soft evidence". When the evidence has been propagated, the beliefs retrieved by the network are used to iterate the prediction for an arbitrary number *P* of steps (see Fig. 7). The final prediction sequence is then summarized as a sequence of normal continuous

Technically, the key challenge is approximating a continuous probability density function by means of a discrete one, as reported in the following. Given mean value and standard

fined by the values *x*0, *x*1, …, *xL* (that corresponds to the *L* + 1 bounds of the *L* inter‐ vals), the probabilities of each interval must be found such that mean value and standard

known, a normal Gaussian probability density function (pdf) takes the form

*i*), and an interval discrete random variable de‐

respectively. As well

Gaussian approx.

> Continuous distribution

\_ and *σ* \_

Interval discrete distribution

Instantiation & propagation

Prediction Loop

Interval discrete distribution

**Figure 5.** Instantiation of a measure as "soft evidence" and corresponding information flow. Instantiation and propa‐

\_ *i*;*σ* \_

Discrete approx.

Continuous distribution

of the prediction error.

distributions *N* (*μi*

*f* (*x*)=1 / 2*πe*-*<sup>x</sup>* <sup>2</sup>

;*σi*

deviation of a normal distribution *N*( *μ*

, - *∞* < *x* <*∞*.

/2

Measure Gaussian

Singleton value

approx.

gation is repeated in the prediction loop (dashed), as described in Fig. 7.

), *i* =1, …, *P*.

16 Dynamic Programming and Bayesian Inference, Concepts and Applications

deviation of the discrete distribution are as close as possible to *μ*

$$\sigma = \max \left( \stackrel{\smile}{\sigma}\_{\prime} \quad \sigma\_{\min} \right)\_{\prime} \quad \sigma\_{\min} = \frac{\underset{x\_{l+1} \quad x\_{l}}{\sqrt{12}}}{\sqrt{12}} \,. \tag{18}$$

For each interval of a standard discrete node, probability is given by *qi* <sup>=</sup> *qi* ' ∑*i*=0 *<sup>N</sup>* -1 *qi* ' where:

$$\dot{q}\_{i} = \mathbb{O}(\boldsymbol{x}\_{i+1}) \text{ - } \mathbb{O}(\boldsymbol{x}\_{i}) \tag{19}$$

For a discrete periodic variable (like a direction) defined in the interval *x*0; *xL* (that for a direction is 0;2*π* ), the pdf can be split into periodic intervals shifted each other by integer multiples of period (*xL* - *x*0):

$$\left\|\mathbf{x}\_{0} + k^{\*}(\mathbf{x}\_{L} \cdot \mathbf{z}\_{0})\right\| \mathbf{x}\_{L} + k^{\*}(\mathbf{x}\_{L} \cdot \mathbf{z}\_{0}) \right\| \cdot k \in \mathbb{Z} \tag{20}$$

If standard deviation is sufficiently smaller than the period *σ*≪(*xL* - *x*0), this infinite series can be approximated by the three central intervals, that is for *k* = - 1,0, 1. Moreover, when *σ* ≅(*xL* - *x*0), probability is almost equally distributed all over the interval and the approxi‐ mation introduced by neglecting the stubs for *k* < - 1, *k* >1 can be simply compensated by renormalizing the resulting distribution. Thus, in this case, probability of each interval is given as in the previous algorithm but the right and left stubs of the Gaussian are overlapped in order to consider periodicity of the discrete variable. For each interval of a periodic discrete node, probability is given by *pi* <sup>=</sup> *pi* ' ∑*i*=0 *<sup>L</sup>* -1 *pi* ' where:

$$\begin{aligned} p\_i \stackrel{\circ}{=} \{ \mathsf{op}(\mathbf{x}\_{i+1}) \cdot \mathsf{op}(\mathbf{x}\_i) \} + \{ \mathsf{op}(\mathbf{x}\_{i+1} + \begin{pmatrix} \mathbf{x}\_L & \mathbf{x}\_0 \end{pmatrix}) \cdot \mathsf{op}(\mathbf{x}\_i + \begin{pmatrix} \mathbf{x}\_L & \mathbf{x}\_0 \end{pmatrix}) \} + \\ \{ \mathsf{op}(\mathbf{x}\_{i+1} \cdot \begin{pmatrix} \mathbf{x}\_L & \mathbf{x}\_0 \end{pmatrix}) \cdot \mathsf{op}(\mathbf{x}\_i \cdot \begin{pmatrix} \mathbf{x}\_L & \mathbf{x}\_0 \end{pmatrix}) \} \end{aligned} \tag{21}$$

because for such complex domains eliciting expert knowledge to learn conditional probability tables is not feasible. Hence several sets of data were generated through simulations prior to the application of the EM learning process whose main features were presented in sub-section

Bayesian Networks for Supporting Model Based Predictive Control of Smart Buildings

http://dx.doi.org/10.5772/58470

19

**•** the first one was made up of randomly generated data, which means that the inputs (e.g. weather, heat gains, occupancy figures etc..) were allowed to vary without additional

**•** the second sample, called "likely" dataset, was generated through simulations whose inputs were allowed to vary within their same ranges cited above, but their differential variations being constrained: it means that the difference between the value of each variable at the present time step and the value of the same variable at the previous time step was limited

**•** the third "typical" sample was built through simulations, whose inputs were taken from real measurements, such as real weather conditions in Barcelona, number of people passing through the station and trains recorded by occupancy surveys, internal heat gains estimated

As the learning process got information from all these samples, the resulting networks not only were expected to be able to respond to sudden variations of external actions and to consider quite unlikely events (which is the case of the first dataset) but to be very sensitive to

The whole building model used for running simulations was developed within the Seam4us research as a lumped parameter model in the DymolaTM simulation environment, that is based on the Modelica language [27]. Starting from a validated library for building simulation developed by the Lawrence Berkeley National Laboratory [28], a specific library for under‐

However, such a model cannot be run in real-time when the controller needs to determine the best candidate control strategies, so it was reduced into the less computationally demanding form of Bayesian Networks. In order for the PdG-L3 model to predict the environmental status of the station, it was split into two Bayesian Networks, each relative to a different physical

**•** temperature prediction dynamic Network (TP-DBN), which is in the form of a DBN, because it forecasts expected temperature in the station given inputs about current and past time

**•** air flow prediction Bayesian network (AF-BN), which is in the form of a regular BN, because it estimates variables relative to air flow in the station and energy consumption of the fans,

Once available, the two networks were run according to the scheme outlined in Fig. 7. Here a portion of the whole process in Fig. 1 is depicted. At every iteration the controller will opportunely query the two of the networks to get future estimations about the variables

the more likely scenarios included in the second and third datasets.

ground station was developed, which is so far in its validation phase.

4.2. Three datasets were generated:

constraints within their range;

on the basis of real measurements etc..

by a threshold;

phenomenon:

steps;

given its current status.

The instantiation strategy described above was used according to the scheme in Fig. 5 for implementing the discrete approximation and smoothing the behaviour of the BNs when inputting evidences. The first Gaussian approximation of Fig. 5 has mean value equal to the measure (eventually calibrated) and standard deviation represented by the uncertainty of the measurement process, that can be determined by the data-sheet of the installed sensor. The second Gaussian approximation is trivially implemented by computing mean and standard deviation of the input distribution. Fig. 6 shows examples of Gaussian distributions approxi‐ mated so as to be assigned as evidences of a Bayesian Network.

Fig. 6: Comparison between the desired Gaussian distribution and the approximated one to **Figure 6.** Comparison between the desired Gaussian distribution and the approximated one to be instantiated in an interval discrete chance node of a BN: example of a non-periodic variable (a) and of a periodic variable (b).

be instantiated in an interval discrete chance node of a BN: example of a non-periodic variable (a)

#### 6.1 The procedure Basically, the development of both Dynamic Bayesian Networks (DBNs) and regular BNs is not an **5. Development of the Bayesian models and cost function**

#### easy task and usually consists of three main phases: 1. definition of the network topology (structural learning); **5.1. The procedure**

2. preparation of the training set and learning of the conditional probability tables; 3. final assessment of the network. In the case object of this chapter, the behaviour of the metro station Passeig de Gracia was first simulated through whole building analyses, which provided datasets encompassing all the possible Basically, the development of both Dynamic Bayesian Networks (DBNs) and regular BNs is not an easy task and usually consists of three main phases:

environmental conditions, including those which are considered as not very likely, then that **1.** definition of the network topology (structural learning);

on the basis of real measurements etc..

6. Development of the Bayesian models and cost function


a threshold;

included in the second and third datasets.

and of a periodic variable (b).

generated: - the first one was made up of randomly generated data, which means that the inputs (e.g. weather, heat gains, occupancy figures etc..) were allowed to vary without additional constraints within their range; - the second sample, called "likely" dataset, was generated through simulations whose inputs In the case object of this chapter, the behaviour of the metro station Passeig de Gracia was first simulated through whole building analyses, which provided datasets encompassing all the possible environmental conditions, including those which are considered as not very likely, then that knowledge was transferred into Bayesian Networks. This was deemed necessary

> were allowed to vary within their same ranges cited above, but their differential variations being constrained: it means that the difference between the value of each variable at the present time step and the value of the same variable at the previous time step was limited by


As the learning process got information from all these samples, the resulting networks not only were expected to be able to respond to sudden variations of external actions and to consider quite unlikely events (which is the case of the first dataset) but to be very sensitive to the more likely scenarios

The whole building model used for running simulations was developed within the Seam4us research as a lumped parameter model in the DymolaTM simulation environment, that is based on the Modelica because for such complex domains eliciting expert knowledge to learn conditional probability tables is not feasible. Hence several sets of data were generated through simulations prior to the application of the EM learning process whose main features were presented in sub-section 4.2. Three datasets were generated:

*pi* '

and of a periodic variable (b).

6.1 The procedure

**5.1. The procedure**

0 0.01 0.02 0.03 0.04

> 0 0.2 0.4 0.6 0.8 1

generated:

=(Φ(*xi*+1) - Φ(*xi*

18 Dynamic Programming and Bayesian Inference, Concepts and Applications

mated so as to be assigned as evidences of a Bayesian Network.

6. Development of the Bayesian models and cost function

1. definition of the network topology (structural learning);

**5. Development of the Bayesian models and cost function**

easy task and usually consists of three main phases:

not an easy task and usually consists of three main phases: **1.** definition of the network topology (structural learning);

3. final assessment of the network.


Discrete approximated ditribution (µ=2.5/σ=10.8)


Desired normal distribution N(2.5;10)

constraints within their range;

on the basis of real measurements etc..

a threshold;

**3.** final assessment of the network.

included in the second and third datasets.

)) + (Φ(*xi*+1 + (*xL* - *x*0)) - Φ(*xi* + (*xL* - *x*0))) +

The instantiation strategy described above was used according to the scheme in Fig. 5 for implementing the discrete approximation and smoothing the behaviour of the BNs when inputting evidences. The first Gaussian approximation of Fig. 5 has mean value equal to the measure (eventually calibrated) and standard deviation represented by the uncertainty of the measurement process, that can be determined by the data-sheet of the installed sensor. The second Gaussian approximation is trivially implemented by computing mean and standard deviation of the input distribution. Fig. 6 shows examples of Gaussian distributions approxi‐

(a) (b)

Fig. 6: Comparison between the desired Gaussian distribution and the approximated one to

0 0.2 0.4 0.6 0.8 1

0 0.01 0.02 0.03 0.04


Discrete approximated ditribution (µ=27.7/σ=14.9)


Desired normal distribution N(30;10)

be instantiated in an interval discrete chance node of a BN: example of a non-periodic variable (a)

**Figure 6.** Comparison between the desired Gaussian distribution and the approximated one to be instantiated in an

interval discrete chance node of a BN: example of a non-periodic variable (a) and of a periodic variable (b).

Basically, the development of both Dynamic Bayesian Networks (DBNs) and regular BNs is not an

In the case object of this chapter, the behaviour of the metro station Passeig de Gracia was first simulated through whole building analyses, which provided datasets encompassing all the possible environmental conditions, including those which are considered as not very likely, then that knowledge was transferred into Bayesian Networks. This was deemed necessary because for such complex domains eliciting expert knowledge to learn conditional probability tables is not feasible. Hence several sets of data were generated through simulations prior to the application of the EM learning process whose main features were presented in sub-section 4.2. Three datasets were




As the learning process got information from all these samples, the resulting networks not only were expected to be able to respond to sudden variations of external actions and to consider quite unlikely events (which is the case of the first dataset) but to be very sensitive to the more likely scenarios

The whole building model used for running simulations was developed within the Seam4us research as a lumped parameter model in the DymolaTM simulation environment, that is based on the Modelica

2. preparation of the training set and learning of the conditional probability tables;

**2.** preparation of the training set and learning of the conditional probability tables;

Basically, the development of both Dynamic Bayesian Networks (DBNs) and regular BNs is

In the case object of this chapter, the behaviour of the metro station Passeig de Gracia was first simulated through whole building analyses, which provided datasets encompassing all the possible environmental conditions, including those which are considered as not very likely, then that knowledge was transferred into Bayesian Networks. This was deemed necessary

(Φ(*xi*+1 - (*xL* - *<sup>x</sup>*0)) - <sup>Φ</sup>(*xi* - (*xL* - *<sup>x</sup>*0))) (21)


As the learning process got information from all these samples, the resulting networks not only were expected to be able to respond to sudden variations of external actions and to consider quite unlikely events (which is the case of the first dataset) but to be very sensitive to the more likely scenarios included in the second and third datasets.

The whole building model used for running simulations was developed within the Seam4us research as a lumped parameter model in the DymolaTM simulation environment, that is based on the Modelica language [27]. Starting from a validated library for building simulation developed by the Lawrence Berkeley National Laboratory [28], a specific library for under‐ ground station was developed, which is so far in its validation phase.

However, such a model cannot be run in real-time when the controller needs to determine the best candidate control strategies, so it was reduced into the less computationally demanding form of Bayesian Networks. In order for the PdG-L3 model to predict the environmental status of the station, it was split into two Bayesian Networks, each relative to a different physical phenomenon:


Once available, the two networks were run according to the scheme outlined in Fig. 7. Here a portion of the whole process in Fig. 1 is depicted. At every iteration the controller will opportunely query the two of the networks to get future estimations about the variables relevant to select the most opportune control policy to be adopted at each running step. To this aim, the networks need to be instantiated first: the current temperature in the station's platform (PL3) and weather conditions will be provided by the permanent monitoring network installed in the station, along with candidate fan frequencies. Given these inputs, the controller is allowed to query the AF-BN in order to estimate fans consumption and air changes in the station at each time step. Such a prediction step takes a few seconds and is performed by the software HuginTM through algorithms for belief propagation, of the kind already reported in sub-section 4.1. Then, the TP-DBN will take these variables as inputs, along with other state variables (e.g. current PL3 temperature, temperature difference between inside and outside and forecasted weather, people, train arrival etc..) in order to predict PL3 temperature at the next time step (which is typically every hour, unless a different control step is required). Again belief propagation is performed with this second network. Then, the same loop will be repeated at each iteration. This loop is based on the assumption that temperature variations between two consecutive hours is so small to be considered constant by the AF-BN without big prejudice. In other words, to the AF-BN purposes temperature at the next time step (T1) was considered equal to temperature at the current time step (T0). This assumption did not compromise the reliability of future estimations, as will be shown by the validation at the end of this sub-section. Both the networks were built following the same methodology:

of nodes required to be included in the Bayesian graphs. This final set was naturally grouped into two sub-clusters of variables: one of them including those related to air flow processes, and the other one including those related to the environmental temperature dynamics. Finally, the qualitative relationships among variables (represented by arcs in the networks) were worked out partly by expert knowledge, and partly using the "structural learning" tool

Bayesian Networks for Supporting Model Based Predictive Control of Smart Buildings

Number of people in station at t+k Number of trains in station at t+k Weather data at t+k

Air change at t+k Temperature at t+k

PREDICTION MODEL

DISTURBANCES MODEL

CONTROLLER STATION

http://dx.doi.org/10.5772/58470

21

supplied by HuginTM, which used the information included in the "random" dataset.

AF-BN TP-DBN

New Cycle: k=k+1

output variables given by the indices outlined in sub-section 5.1;

**Figure 7.** The basic loop of the predictive cycle adopted for PdG-L3 which involves both Bayesian Networks. It imple‐ ments the prediction model in Fig.1, which starts by setting a counter at 1 and loops until the counter reaches any

The second step was the most critical in the development process, and helped in several tasks:

**•** meaning and dependencies between nodes have been reviewed according to the relation‐ ships suggested by physical laws (e.g. PL3 temperature is affected by the amount of outside air supplied by ventilation and by the temperature difference between indoors and out‐

**•** the number of intervals for discretizing the state space pertaining to every node, which was refined through an iterative process, with the final purpose of minimizing the errors of the

**•** a few links have been rearranged, in order to improve the performances of the two networks.

Fig. 8-a depicts the final structure of the dynamic Bayesian network (i.e. TP-DBN), which was used to predict PL3's temperature in PdG station in the next step (node TemPL3\_p01), starting from inputs such as: forecasted number of people in the station at the next step (NPeSta\_p01), forecasted internal gains supplied by trains at the next step (GaiTr\_p01), current PL3's temperature (TemPL3), forecasted outdoor temperature (TOuMet\_p01), forecasted air changes per hour (ACOPL3\_p01) and deviation of temperature from the past time step (DTePL3). The

Fan absorbed power at t+k

Fan driving frequencies at t+k

Weather data at t+k

desired prediction horizon.

doors, which became its parents);

Variation of temperature at t+k-1 Temperature at t+k-1


The first step started from the analysis of 81 variables included in the DymolaTM dataset, properly filtered and re-sampled at a 60 sec rate. Iterative cluster analyses [29] were useful to group those variables into clusters and determine those which were redundant [30], because providing the same information given by others (e.g. surface temperatures in PL3 were highly correlated and grouped in the same level, hence just one of them was kept and the remaining ones were dropped out).

This first step allowed us to cut down the number of involved variables to 39. Then, a further reduction was done, where only those variables showing the strongest dependence relation‐ ships were kept, cutting the whole number down to 25, hence finding the minimum number of nodes required to be included in the Bayesian graphs. This final set was naturally grouped into two sub-clusters of variables: one of them including those related to air flow processes, and the other one including those related to the environmental temperature dynamics. Finally, the qualitative relationships among variables (represented by arcs in the networks) were worked out partly by expert knowledge, and partly using the "structural learning" tool supplied by HuginTM, which used the information included in the "random" dataset.

relevant to select the most opportune control policy to be adopted at each running step. To this aim, the networks need to be instantiated first: the current temperature in the station's platform (PL3) and weather conditions will be provided by the permanent monitoring network installed in the station, along with candidate fan frequencies. Given these inputs, the controller is allowed to query the AF-BN in order to estimate fans consumption and air changes in the station at each time step. Such a prediction step takes a few seconds and is performed by the software HuginTM through algorithms for belief propagation, of the kind already reported in sub-section 4.1. Then, the TP-DBN will take these variables as inputs, along with other state variables (e.g. current PL3 temperature, temperature difference between inside and outside and forecasted weather, people, train arrival etc..) in order to predict PL3 temperature at the next time step (which is typically every hour, unless a different control step is required). Again belief propagation is performed with this second network. Then, the same loop will be repeated at each iteration. This loop is based on the assumption that temperature variations between two consecutive hours is so small to be considered constant by the AF-BN without big prejudice. In other words, to the AF-BN purposes temperature at the next time step (T1) was considered equal to temperature at the current time step (T0). This assumption did not compromise the reliability of future estimations, as will be shown by the validation at the end

20 Dynamic Programming and Bayesian Inference, Concepts and Applications

of this sub-section. Both the networks were built following the same methodology:

and discern the strongest relationships among the variables themselves;

constraints on all the possible weather conditions and disturbance actions;

the final shape of the networks;

networks.

ones were dropped out).

**1.** first structural learning: it was determined first by the a-priori knowledge from the researchers, who ordered and connected the variables according to the physical relation‐ ships among them; then statistical analyses were used to improve the assumed structure

**2.** improvement of the network's structure: this was carried out through analysis of its performance indices, after learning conditional probability tables from the random dataset, which was the first set of data generated by running the DymolaTM model without

**3.** final refinement using the two additional datasets: adding more datasets allowed the developers to quantify even probabilistic relationships among the variables and to define

**4.** final evaluation of the networks: as for the three steps above, also in this final step the indices defined in paragraph 5.1 were used to evaluate the prediction capabilities of the

The first step started from the analysis of 81 variables included in the DymolaTM dataset, properly filtered and re-sampled at a 60 sec rate. Iterative cluster analyses [29] were useful to group those variables into clusters and determine those which were redundant [30], because providing the same information given by others (e.g. surface temperatures in PL3 were highly correlated and grouped in the same level, hence just one of them was kept and the remaining

This first step allowed us to cut down the number of involved variables to 39. Then, a further reduction was done, where only those variables showing the strongest dependence relation‐ ships were kept, cutting the whole number down to 25, hence finding the minimum number

**Figure 7.** The basic loop of the predictive cycle adopted for PdG-L3 which involves both Bayesian Networks. It imple‐ ments the prediction model in Fig.1, which starts by setting a counter at 1 and loops until the counter reaches any desired prediction horizon.

The second step was the most critical in the development process, and helped in several tasks:


Fig. 8-a depicts the final structure of the dynamic Bayesian network (i.e. TP-DBN), which was used to predict PL3's temperature in PdG station in the next step (node TemPL3\_p01), starting from inputs such as: forecasted number of people in the station at the next step (NPeSta\_p01), forecasted internal gains supplied by trains at the next step (GaiTr\_p01), current PL3's temperature (TemPL3), forecasted outdoor temperature (TOuMet\_p01), forecasted air changes per hour (ACOPL3\_p01) and deviation of temperature from the past time step (DTePL3). The

The third step was aimed at performing further refinement using the "typical" and "likely" datasets, which were generated through DymolaTM. Technically, that means that the EM learning algorithm was performed by adding the information included in these two datasets to the information already derived from the "random" dataset. This process allowed to include more information in the networks, mainly about those scenarios which are likely to occur more often. So it helped make the model more accurate to predict those states which are more likely

Bayesian Networks for Supporting Model Based Predictive Control of Smart Buildings

http://dx.doi.org/10.5772/58470

23

The refinement was mainly performed in terms of tuning the subdivision into intervals of the nodes and in terms of converting discrete variable into continuous variables, should they work better. The steps no. 2 and 3 required many iterations of learning, refining and validating.

Each loop was characterized by either a modification of the structure or of any network's node and conditional dependencies. Then the modified network's version was evaluated to define whether the modification had to be kept (in case it decreased the errors) or had to be rejected

On the whole and just to give an example, 140 cycles were made with the TP-DBN (Tab. 1), which was useful to reduce the error from 4.98 ° to 0.72 °C (in the *RMSE* case) and from 17% to 4% (in the *NRMSE* case). The trend during the refinement process led to a continuous increment of performances, as shown in Tab. 1. In addition, 82 cycles were needed to optimize the AF-BN: for the control variable, Station Fan Power (PElSF1) *RMSE* fell down from 1858 W

In the TP-DBN all the nodes were represented by discrete variables. In the AF-BN all the variables were continuous except the following ones: frequencies of fans (DFreTF1\_p01,

1 4.98 17

54 3.32 12

98 1.81 6

114 1.00 4

140 0.72 4

**Table 1.** Gradual improvement of performances during continuous refinement of the TP-DBN network.

**RMSE (°C) NRMSE (%)**

(in case it increased the errors), being the errors estimated as explained in section 5.1.

to 377 W, whereas *NRMSE* fell from 10.3% down to 2.3%.

…

…

…

…

DFreTF2\_p01, DFreSF1\_p01) and wind direction (WiDMet\_p01).

**Cycle no. TemPL3\_p01**

to occur.

**Figure 8.** Predictive and dynamic Bayesian Network relative to temperature in PL3 (a) and predictive Bayesian Net‐ work relative to air flow changes in the station (b).

network's intermediate variables are useful to perform computations and simplify conditional probabilistic relationships among variables.

Similarly holds with the AF-BN network (Fig. 8-b), whose inputs are: forecasted frequencies of fans in the station and tunnels at the next time step (DFreTF1\_p01, DFreTF2\_p01, DFreSF1\_p01), forecasted internal gains by trains (GaiTr1\_p01), forecasted wind direction and speed (WiDMet\_p01, WiSMet\_p01), outdoor temperature (TOuMet\_p01) and current temper‐ ature (TemPL3). The main outputs are the power consumption of fans – in the station (PElSF1) and in the tunnels (PElTF1, PElTF2)-and air flow rates expected across the corridors leading to PL3: AFlCNl\_p01 (corridor CNl), AFlCNop\_p01 (sum of corridors CNo and CNp), AFlCNq\_p01 (corridoio CNq) and AFlSlb\_p01 (station link). These estimated airflows are then summed up coherently to the Air Mass Balance for computing the overall air change in PL3 (ACOPL3), needed as input from the TP-DBN (Fig. 7).

During this process, whenever expert knowledge deemed more than one structure as reason‐ able for a certain network or fragment of it, the best one was chosen by evaluating which of them returned the highest performances according to the indices described in sub-section 5.1, after learning from the random datasets, which is the one where the widest variations and fastest dynamics of variables were considered. Even the other aforementioned amendments (e.g. rearrangement of links, optimal discretization of the networks' nodes etc..) were tested through optimization of the performances indices defined in sub-section 5.1.

The third step was aimed at performing further refinement using the "typical" and "likely" datasets, which were generated through DymolaTM. Technically, that means that the EM learning algorithm was performed by adding the information included in these two datasets to the information already derived from the "random" dataset. This process allowed to include more information in the networks, mainly about those scenarios which are likely to occur more often. So it helped make the model more accurate to predict those states which are more likely to occur.

The refinement was mainly performed in terms of tuning the subdivision into intervals of the nodes and in terms of converting discrete variable into continuous variables, should they work better. The steps no. 2 and 3 required many iterations of learning, refining and validating.

Each loop was characterized by either a modification of the structure or of any network's node and conditional dependencies. Then the modified network's version was evaluated to define whether the modification had to be kept (in case it decreased the errors) or had to be rejected (in case it increased the errors), being the errors estimated as explained in section 5.1.

On the whole and just to give an example, 140 cycles were made with the TP-DBN (Tab. 1), which was useful to reduce the error from 4.98 ° to 0.72 °C (in the *RMSE* case) and from 17% to 4% (in the *NRMSE* case). The trend during the refinement process led to a continuous increment of performances, as shown in Tab. 1. In addition, 82 cycles were needed to optimize the AF-BN: for the control variable, Station Fan Power (PElSF1) *RMSE* fell down from 1858 W to 377 W, whereas *NRMSE* fell from 10.3% down to 2.3%.

In the TP-DBN all the nodes were represented by discrete variables. In the AF-BN all the variables were continuous except the following ones: frequencies of fans (DFreTF1\_p01, DFreTF2\_p01, DFreSF1\_p01) and wind direction (WiDMet\_p01).

network's intermediate variables are useful to perform computations and simplify conditional

**Figure 8.** Predictive and dynamic Bayesian Network relative to temperature in PL3 (a) and predictive Bayesian Net‐

Similarly holds with the AF-BN network (Fig. 8-b), whose inputs are: forecasted frequencies of fans in the station and tunnels at the next time step (DFreTF1\_p01, DFreTF2\_p01, DFreSF1\_p01), forecasted internal gains by trains (GaiTr1\_p01), forecasted wind direction and speed (WiDMet\_p01, WiSMet\_p01), outdoor temperature (TOuMet\_p01) and current temper‐ ature (TemPL3). The main outputs are the power consumption of fans – in the station (PElSF1) and in the tunnels (PElTF1, PElTF2)-and air flow rates expected across the corridors leading to PL3: AFlCNl\_p01 (corridor CNl), AFlCNop\_p01 (sum of corridors CNo and CNp), AFlCNq\_p01 (corridoio CNq) and AFlSlb\_p01 (station link). These estimated airflows are then summed up coherently to the Air Mass Balance for computing the overall air change in PL3

During this process, whenever expert knowledge deemed more than one structure as reason‐ able for a certain network or fragment of it, the best one was chosen by evaluating which of them returned the highest performances according to the indices described in sub-section 5.1, after learning from the random datasets, which is the one where the widest variations and fastest dynamics of variables were considered. Even the other aforementioned amendments (e.g. rearrangement of links, optimal discretization of the networks' nodes etc..) were tested

through optimization of the performances indices defined in sub-section 5.1.

probabilistic relationships among variables.

22 Dynamic Programming and Bayesian Inference, Concepts and Applications

work relative to air flow changes in the station (b).

(ACOPL3), needed as input from the TP-DBN (Fig. 7).


**Table 1.** Gradual improvement of performances during continuous refinement of the TP-DBN network.

*J* =∑*<sup>k</sup>* =1

Subject to constraints:

*ACOPl*3(*k*) > *ACOPl*3*Min*

*TemPl*3(*k*) <*TemPl*3*Max*

*<sup>H</sup> <sup>α</sup>PT* ( <sup>|</sup>*PElTF* 1(*<sup>k</sup>* ) <sup>+</sup> *PElTF* 2(*<sup>k</sup>* )|

<sup>+</sup>*αDT* ( *TOuWS* (*<sup>k</sup>* ) - *TemPl*3(*<sup>k</sup>* ) *DT*

*AC*

<sup>+</sup>*αAC*( *ACOPl*<sup>3</sup>

*FreSFMin* < *FreSF* 1(*k*)= *FreSF* 2(*k*) < *FreSFMax*

*TemPl*3*Max*, *FreSFMin*, *FreSFMax*.

**5.3. The networks and their estimates**

the desired values denoted with bar notation *TemPl*3

\_

2*PT*


<sup>~</sup> )

<sup>~</sup> )

2

2

<sup>~</sup> ) <sup>+</sup> *<sup>α</sup>PS* ( <sup>|</sup>*PElSF* 1(*<sup>k</sup>* ) <sup>+</sup> *PElSF* 2(*<sup>k</sup>* )|

<sup>+</sup> *<sup>α</sup><sup>T</sup>* ( *TemPl*<sup>3</sup> \_

The variables marked with "tilde" (~) are the normalisation coefficients that corresponds to the typical values of the corresponding variable. At sampling time *t* ∈N, *PElTF* 1(*k*) represents the estimation of the value of variable *PElTF* 1 at time *t* + *k* (*k* ∈N) evaluated at time *t* (i.e. *PElTF* 1(*t* + *k* |*t*)). The design parameters of the MPC controller are the prediction horizon *H* ,

objective in the cost function *α*…, and the bounds of the given constraints: *ACOPl*3*Min*,

The "*soft evidence*" instantiation strategy presented in section 5.2 has been implemented in a java library that wraps HuginTM reasoning engine and allows also for other high-level func‐ tionalities, like multiple network iterations and interconnection between different networks that shares the same variables. This library is able to get the set of variables describing the current state of the station and to initialize with them the first network to be queried. Then, it is able to perform probabilistic inference by running the HuginTM application and to extract the outputs, which will be transferred to the second network to be queried according to the procedure already shown in Fig. 7. This is done *H* times, if *H* is the desired prediction horizon. The same functions has been also integrated in an excel spreadsheet for validating BNs. The two networks were combined according to the scheme depicted in Fig. 7 and used to simulate the predictive control. Fig. 10 shows the prediction results for one typical day (in blue), i.e. prediction horizon *H* =24 hours, compared with the simulation results achievable from the Dymola model (shown in red), regarding the main outputs: expected energy consumption by one of the station' fans (6a), overall airflows coming from outdoor air (6b) and future temper‐ ature (6c) plots in the platform of Line 3. In order to assess the level of accuracy of such predictions, Table 2 shows the corresponding *RMSE* and *NRMSE* relative to the variables

\_

, *ACOPl*3 \_

, the weights of each single

2*PS*

Bayesian Networks for Supporting Model Based Predictive Control of Smart Buildings


<sup>~</sup> )

<sup>~</sup> )

*T*

<sup>+</sup> *<sup>α</sup>DF* ( *FreSF* 1(*<sup>k</sup>* ) - *FreSF* 1(*<sup>k</sup>* - 1) *DF*

<sup>~</sup> ) <sup>+</sup>

http://dx.doi.org/10.5772/58470

2

(23)

25

2 +

**Figure 9.** Qualitative comparison between the real temperature plot computed by DymolaTM and forecasts by the Bayesian Network TP-DBN.

Finally, and as fourth step, the performances of the two networks were verified also through simulations. Fig. 9 shows the good agreement between the real temperature simulated by DymolaTM in PL3 and the forecasted plot of PL3 as predicted by TP-DBN. The simulations performed by the Bayesian Networks in this case were carried out according to what already described. The input values at the first time step were instantiated as evidences taken by the Dymola model. Then, the outputs from the networks were used as inputs for the next time step in the networks and the simulations were iterated in the same way all over the period shown in the diagram. It's clear that the predictive and dynamic Bayesian network are able to accurately model the temperature plot in PL3 and to give the right inputs to the controller, in order to evaluate the best control policy.
