Application of stochastic model predictive control for building energy systems using latent force models

2.1 State-space transformation of latent force models

Consider the continuous-time stochastic system

with state x ∈ R N x , control input u ∈ R N u , disturbance d ∈ R N d , system dynamics function f , a vector of stochastically independent Wiener processes w x ∈ R N w , and a diffusion matrix B x ∈ R N x × N w . The disturbance d cannot be measured directly, but it is assumed that each element can be modeled as a time-dependent GP

with a mean function m _i (t) and a kernel k _i (t, t′). In the following, m _i (t) = 0 is assumed without loss of generality [6].

GPs are not suitable for the control of systems with short sampling times due to the large computational effort. Therefore, an alternative approach is proposed in [7], where the GPs are reformulated as linear state-space models, for which the state can be estimated using the Kalman filter. In contrast to Gaussian Processes, the iterative procedure of the Kalman filter does not require saving large data sets and the computation time is time-independent. Provided that the spectral density of the kernel is of the form

the i-th element of the disturbance vector (2) can be modeled by the system [7]

with an internal state z i ∈ R N z , a Wiener process w _z,i , parameters of the system matrix a = [ a 1 , … , a N z ] , and the diffusion parameter σ. The number of internal states N _z as well as the parameters a and σ depend on the kernel k _i (t, t′).

A class of kernels that satisfies condition (3) is the Matérn kernel

with δ = |t − t′|, hyperparameters σ M 2 and l, Bessel function of second kind K _ν , gamma function Γ, and parameter ν = 0.5 + n with n ∈ N 0 [7]. Some kernels that do not fulfill condition (3) can be approximated by a state-space representation (4). For example, the periodic kernel

with hyperparameters σ P 2 , ω ₀, and l can be approximated using a series of stochastic resonators [15]. In the same way that kernels can be combined with each other, the state-space representations of the GP priors can be combined as well. The quasi-periodic kernel k _QP (t, t′) = k _M (t, t′)k _P (t, t′), for instance, can be approximated by the system matrix

with Kronecker product ⊗, identity matrix I , and number of internal states n _M and n _P , respectively [15]. If multiple disturbances affect system (1), each GP (2) can be modeled by a linear time invariant system (4), which can be combined to the extended system

where z = [ z 1 T , … , z N d T ] T denotes the augmented state vector, w z = [ w z , 1 T , … , w z , N d T ] T denotes the vector of independent Wiener processes, and blkdiag(⋅) denotes the block diagonal matrix.

2.2 Stochastic model predictive control

The system dynamics (1) involve Wiener processes and Gaussian Processes that cause uncertainty. Further uncertainty is introduced, if the states of the system cannot be measured directly but must be estimated. By combining the system dynamics (1) and the disturbance dynamics (8) to an extended system, the system state x and the disturbance state z can be estimated together, for example using the unscented Kalman filter [16]. Furthermore, a model predictive controller requires predicted trajectories of the disturbance. Since (8) is independent of x and u , the mean of the disturbance state μ z = E z and the covariance matrix of the disturbance state Σ _d = Cov[ d ] can be predicted by integrating

where the initial mean and covariance are given by the state estimator. These can be used to calculate the mean of the disturbance μ d = E d and the covariance matrix of the disturbance by

The estimated state, the predicted disturbance and the system model (1) can be utilized to construct a stochastic model predictive controller that solves an optimal control problem (OCP) in every time step. Since x is a random variable, constraints are formulated as chance constraints that must be fulfilled with a certain probability. Therefore, the stochastic OCP

is considered with terminal cost function V, integral cost function l, prediction horizon T, chance constraints h _i with probabilities α _i , and terminal chance constraints h _T,j with probabilities α _T,j . An exemplary control loop is shown in Figure 1. Here, the Observer estimates the mean and covariance of the states x and the internal states z . The disturbance prediction is computed by the Disturbance predictor via (9) and (10).

Figure 1:

Exemplary control loop for the LFM-SMPC. The observer provides state estimates for the mechanistic model and the latent states, which are used to compute the estimated disturbance trajectory d ̂ .

Solving the stochastic OCP (11) requires the propagation of the uncertain state x through the potentially nonlinear system dynamics function f . The temporal change of the probability density function of x can be represented by a partial differential equation known as the Fokker-Plank equation [17]. Since solving this equation is difficult for general functions f , several methods, such as the unscented transformation, Gaussian quadrature, and polynomial chaos expansion, have been proposed that approximate the solution of the Fokker-Plank equation. For example, the unscented transformation uses sigma-points to capture the uncertainty. They are computed with

where μ _x and Σ _x denote the mean and covariance of the state x . The additional parameters α _UT and κ specify the spread of the sigma points. To recompute the mean and covariance, the unscented transformation uses weights which are applied to the sigma points. These weights are obtained via

with the additional parameter β. Additional information about these methods can be found in [18].

Overall, tracking the covariance for the stochastic OCP results at least in a computational demand of O n x 2 . The chance constraints (11c) and (11d) can be approximated as well by

where z(α _i ) and z(α _T,j ) are coefficients that can be chosen according to Chebyshev’s inequality [19].

3.1 Modeling solar radiation

Solar radiation has a major impact on the temperature inside buildings. Thus, considering the heat gain from radiation is a great opportunity to reduce the energy demand. Computing the aforementioned heat gain introduces several sources of uncertainty, starting at the radiation measurements. The required quantities are Direct Normal Irradiation (DNI), Diffuse Horizontal Irradiation (DHI) and Global Horizontal Irradiation (GHI). DNI is the radiation on a surface normal to the sun. It is measured with a solar tracker and a pyrheliometer which is labor-intensive. DHI is the diffuse, reflected radiation from the surrounding. Here, a pyranometer can be used with a shading device to block the direct sun beam. GHI is the global radiation and is measured with a pyranometer placed horizontally on a surface [23], [24]. These measured variables are related with

where z is the sun zenith angle [25]. Since GHI is significantly easier to obtain, the majority of on-site measurements are GHI and no information about DNI and DHI is available. Many correlations exist to obtain estimates for DNI and DHI from GHI. A comparison can be found in [25].

The tilted irradiation I _t in W m 2 on a 90° surface, e.g. exterior walls, can be computed with

with the angle θ as the sun altitude, R _d as the diffuse transposition factor of ground reflection, ρ for the foreground albedo and R _r as the transposition factor for ground reflection. Under the isotropic assumption, the transposition factors R _r and R _d can be computed with

where s is the tilt angle of the surface. Note that for R _d this is a rather strong assumption which is outlined in [25].

With the tilted irradiation I _t , one can compute the heat gain Q ̇ sol . The physically correct equations are more involved, which can be seen in [26]. Since the radiation is used as a disturbance in SMPC, a simple approximation with

is used, where A _win is the window surface area with the reflection factor R _win. Note that the explicit time dependence of (24) and (26) was dropped for brevity.

4.1 Use cases

Two BES are taken into account: An example building in EnergyPlus and a parametric model of a floor at the Bosch Research Campus in Renningen, Germany. Ideal HVAC equipment was used for both use cases, since only the energy demand of the zones were of interest. The differences are outlined in the following.

4.1.1 EnergyPlus model

EnergyPlus is a widely adopted BES simulation software used by engineers to e.g. estimate and configure the HVAC equipment [13]. It provides a complex calculation of the thermal properties of a building, including the various disturbances and interactions with the HVAC components. EnergyPlus supplies typical weather data for different regions on earth, which are used during simulation. Thus, it is a white-box modelling approach for BES.

For the simulations in EnergyPlus, the building in Figure 3 is considered. It consists of three adjacent zones, the south, hallway and north zone, where the south and north zones are offices. For the control, every zone is modeled by (20) with identified values for the heat capacities and resistances from simulated data. As said before, EnergyPlus provides typical weather data for different climate zones. Here, typical data from Munich was used. Through the typical weather data, values for DNI, DHI and GHI are available for computing the heat gain through solar radiation.

Figure 3:

3D view of the example building in EnergyPlus.

4.1.2 Renningen office floor

The Renningen model is based on the floor plan shown in Figure 4. It also uses the model in (20) without explicit zone couplings and identified parameters from real world measurements. These measurements include weather data, where only the ambient temperature and GHI are available. From all available zones, only 9, 10, 21 and 23 are used, since 11, 12 and 13 are not subject to solar radiation. Zone 22 and 24 are not used since the thermal characteristics match zone 23 and 21, respectively. From this real world example information about the used RB controller is available, which is used for comparison in the evaluation. It consists out of two PI-Controllers, one for the radiator and one for the ventilation. The setpoint changes between night and day are smoothed with a ramp, so no abrupt changes occur, like in [4].

Figure 4:

Zone layout of the floor at the Bosch research campus in Renningen.

4.2 Solar radiation prediction

From the building perspective, solar radiation, like ambient temperature, is a time dependent disturbance. As seen in Section 3.1, DNI, DHI and GHI are required for computing the heat gain Q ̇ sun . For Renningen only GHI is available, thus a correlation according to Orgill et al. [27] is used to compute DHI and (23) for computing DNI.

The computed heat gain for one zone in Renningen can be seen in Figure 5. Here, the signal properties of the disturbance are evident. A daily periodic trend with variations can be seen with a rough shape. Thus, a combination of a periodic and a Matérn kernel, a quasi-periodic kernel, is a suitable choice. For the Matérn kernel, ν = 1 2 is used because of the roughness. The periodic kernel is approximated with internal 22 states, which is the same for the used quasi-periodic kernel.

Figure 5:

Computed heat gain through solar radiation with α = 1.

4.3 SMPC settings

The stochastic model predictive controller is implemented using the framework GRAMPC-S [28]. GRAMPC-S provides several approximation methods to propagate the uncertainties of the predicted states. In our evaluation, the unscented transformation with α _UT = 0.1, β = 2 and κ = 1 is used to track the variance of the system dynamics. A diffusion term of 1E−5 K is added to both systems for the uncertainty propagation which resembles uncertain system dynamics. The parameters for the MPC models were identified from simulation data.

The temperature of each zone is subject to the constraints

where T _u is the upper and T _l the lower bound. The bounds are found in Table 1 which specify the comfort range during the day and energy efficiency for the night. Lower temperatures can result in condensation and thus yield mold, which is undesired. For the SMPC, Gaussian chance constraint approximations with the probabilities α ∈ 95 % , 80 % , 50 % are used. Note that for a Gaussian chance constraint the 50 % case is equivalent to a deterministic MPC [18] since the quantity z(α _i ) is equal to zero. Specifying the comfort range as chance constraints is sensible, since small constraint violations are tolerable. The integral cost function is given by

for every control input u _i to ensure energy efficiency.

The solar radiation is considered in two different ways: As a constant disturbance and with an LFM. An SMPC without disturbance knowledge is used as a baseline. Combined with the different probabilities of constraint fulfillment, 9 different controllers are considered and simulated for a whole year. They are compared based on the zones’ energy demand and thermal comfort levels. The discomfort is assessed with the integrated constraint violation given in K h. The different settings of the SMPC such as sample time and prediction horizon are shown in Table 2.

4.4 Renningen results

The results for Renningen are given in Figure 6. Additionally, a time series plot is shown in Figure 7 which shows the SMPC with α = 95 % with a constant disturbance prediction. Overall for α = 95 % result in the highest energy consumption followed by α = 80 % and α = 50 %. When looking at the discomfort, α = 95 % achieve the lowest followed again by α = 80 % and α = 50 %. Including a constant disturbance prediction reduces both energy demand and discomfort, except for zone 21. Here a small increase in Energy consumption occurs. With the constant disturbance prediction, α = 50 % results in the lowest energy consumption for each zone. For α = 95 % in the zones 9 and 10, this effectively results in no discomfort. The effect for zone 21 and 23 is much smaller than expected. Further investigation shows, that the ventilation for these zones is at its upper control limit. This can also be seen in Figure 7. The influence of the solar radiation is greater, since the windows of zone 21 and 23 are facing south. Nevertheless, the discomfort is at least halved with the constant disturbance prediction.

Figure 6:

Renningen results for a prediction horizon of 1 h. The RB controller is depicted by .

Figure 7:

Time series result for zone 21 with α = 95 % and constant disturbance prediction. The lower plot shows the controls u _a () and u _w (). The shaded area depicts the admissible temperature range.

Using an LFM to stochastically predict the solar heat gain results in a higher or comparable energy consumption for the zones. For zone 9 and 10 this results in a marginally higher energy consumption than the constant prediction. In zone 21 and 23, the discomfort for α ∈ 95 % , 80 % , 50 % is about 30 K h lower than the constant disturbance prediction. On the other hand, the energy consumption is 200–300 kW h higher.

Comparing with the RB controller, the SMPC at least reduces the energy demand for every zone. Including disturbance information leads to an additional reduction in discomfort.

Overall, the LFM results in a lower discomfort at the price of marginal higher energy consumption compared to the constant disturbance prediction. When a low discomfort is desired, the LFM is clearly the best choice. If energy efficiency is the biggest concern, a constant disturbance prediction with α = 50 % performs best. Note again that for α = 50 % and a Gaussian constraint approximation, the SMPC result is equivalent to a deterministic MPC, see Section 4.3.

4.5 EnergyPlus results

The EnergyPlus use case was simulated with a time horizon of 2 h with the results given in Figure 8. A similar trend like in the Renningen case can be seen. The constant disturbance prediction achieves a lower discomfort for α = 50 % compared to the LFM. In contrast, the energy consumption is slightly larger. Using an LFM results in a considerably lower discomfort for α = 95 % and α = 80 % with a higher energy consumption of about 4–8 %.

Figure 8:

EnergyPlus results for a prediction horizon of 2 h.

To further assess the performance of the SMPC, the complementary cumulative constraint violation is shown in Figure 9. For the south, the SMPC without disturbance information has a significantly greater proportion of high constraint violations. Including a constant disturbance prediction reduces this probability significantly. This is different for the northern zone. Here, the SMPC is grouped in clusters of their chance constraint approximation percentage, e.g. for α = 95 % the performance between no disturbance prediction, constant prediction and LFM is similar. Thus, the advantage of adding disturbance information is smaller for the north zone than for the south zone.This meets the expectation that solar radiation has a lower impact on the northern zone. Moreover, model uncertainties are more dominant in the north, resulting in a higher discomfort for α = 50 %.

Figure 9:

Complementary cumulative distribution of the constraint violation for 2 h.