Analyzing Occupancy-Driven Thermal Dynamics in Smart Buildings

TL;DR

This paper introduces a probabilistic modeling framework for analyzing thermal dynamics and evaluating heating strategies in smart buildings, aiming to optimize energy efficiency while maintaining occupant comfort.

Contribution

It presents a novel Markov reward-based framework for modeling occupancy-driven thermal dynamics and compares multiple heating strategies using probabilistic model checking.

Findings

Selective heating strategy is 13.5 times more cost-effective.

Framework effectively evaluates and compares control strategies.

Occupancy patterns significantly influence thermal management efficiency.

Abstract

The fact that a proper HVAC control strategy can reduce the energy consumption of a building by up to 45% has driven significant research in demand-based HVAC control. This paper presents a novel framework for modeling and analysis of thermal dynamics in smart buildings that incorporates building's thermal properties, a stochastic occupancy model and heating strategies. Each zone of a building is modeled with the help of discrete time Markov rewards formalism where the states represent the occupancy of that zone (either occupied or empty), and the state rewards incorporate the thermal dynamics and heating strategy. To demonstrate the applicability of our proposed framework, we evaluate and compare six different heating strategies for the two zone scenario of a university building. The obtained quantitative results from the PRISM probabilistic model checker show that one of the evaluated control strategies (viz. selective strategy) satisfies our requirement in terms of maintaining the occupants' comfort while being up to 13.5 times more cost effective when compared to the other evaluated strategies. Such evaluations demonstrate the framework's ability to assist in selecting the control strategy tailored around the occupancy pattern and building's thermal property.

Tables

Table 1. Table 1 : Example of a time-dependent occupancy dataset

Day	9am				10am-4pm	5pm
	$v\text{ to }f$	$v\text{ to }v$	$f\text{ to }v$	$f\text{ to }f$	…..	$v\text{ to }f$	$v\text{ to }v$	$f\text{ to }v$	$f\text{ to }f$
1	✓		NA	NA	…..	✓
2		✓	NA	NA	…..	✓
3		✓	NA	NA	…..		✓
4		✓	NA	NA	…..		✓
5		✓	NA	NA	…..			✓
6	✓		NA	NA	…..			✓
7	✓		NA	NA		…..			✓
….	…..	…..	…..	…..	…..	…..	…..	…..	…..
….	…..	…..	…..	…..	…..	…..	…..	…..	…..
179	✓		NA	NA	…..				✓
180		✓	NA	NA	…..		✓

Table 2. Table 2 : Heating strategies

Strategies	Heating of zone 1	Heating of zone 2
S1	Off	Off
S2	On	On
S3	On	Off
S4	Off	On
S5	Alternating	Alternating
S6	Selective	Selective

Table 3. Table 4 : Cost analysis for different strategies

Strategies	Total electricity consumed			Total
	Economy	Off-peak	Peak	cost (£)
S1	NA	NA	NA	0
S2	2kW	6kW	8kW	2.70
S3	1kW	3kW	4kW	1.35
S4	1kW	3kW	4kW	1.35
S5	2kW	4kW	4kW	1.30
S6	2kW	NA	NA	0.20

Equations

$$C_i\frac{d{T}_i(t)}{dt}=\sum _{j\in {\mathcal{N}}_i}\frac{T_j(t)-T_i(t)}{R_i}+Q_i(t)$$

$$C_1\frac{d{T}_1(t)}{dt}$$

$$C_2\frac{d{T}_2(t)}{dt}$$

$$Q_{in{t}_i}(t)=\{\begin{array}{ll}{Q}_{in{t}_i},& \text{if the room is occupied}\\ 0,& \text{otherwise}\end{array}$$

$$Q_{ra{d}_i}(t)=\{\begin{array}{ll}{Q}_{ra{d}_i},& \text{if the radiator is on}\\ 0,& \text{otherwise}\end{array}$$

$$\dot{T}=\hat{A}T+\hat{B}U$$

$${\hat{A}}_{ij}=\{\begin{array}{ll}0,& \text{if }j\ne i,(i,j)\notin \mathcal{E}\\ \frac{1}{(C_i{R}_i)},& \text{if }j\ne i,(i,j)\in \mathcal{E}\\ -{\sum }_{k\ne i}{\hat{A}}_i,& \text{if }j=i,(i,j)\in \mathcal{E}\end{array}$$

$${\hat{B}}_{ij}=\{\begin{array}{ll}\frac{1}{(C_i)},& \text{if }i=j,(i,j)\in \mathcal{E}\\ 0,& \text{otherwise}\end{array}$$

$$\begin{split}\left[\begin{array}[]{cc}\dot{T}_{1}\\ \dot{T}_{2}\end{array}\right]&=\left[\begin{array}[]{cc}\frac{-1}{R_{1}C_{1}}&\frac{1}{R_{1}C_{1}}\\ \frac{1}{R_{2}C_{2}}&\frac{-1}{R_{2}C_{2}}\end{array}\right]\left[\begin{array}[]{cc}T_{1}\\ T_{2}\end{array}\right]+\left[\begin{array}[]{cc}\frac{1}{C_{1}}&0\\ 0&\frac{1}{C_{2}}\\ \end{array}\right]\left[\begin{array}[]{cc}{Q_{1}}\\ {Q_{2}}\end{array}\right]\\ &=\left[\begin{array}[]{cc}\frac{-1}{R_{1}C_{1}}&\frac{1}{R_{1}C_{1}}\\ \frac{1}{R_{2}C_{2}}&\frac{-1}{R_{2}C_{2}}\end{array}\right]\left[\begin{array}[]{cc}T_{1}\\ T_{2}\end{array}\right]+\left[\begin{array}[]{cc}\frac{Q_{1}}{C_{1}}\\ \frac{Q_{2}}{C_{2}}\\ \end{array}\right]\\ \end{split}$$

$$A=\left[\begin{array}{cc}0.7001& 0.2999\\ 0.3007& 0.6993\end{array}\right],B=\left[\begin{array}{cc}0.7299& 0\\ 0& 1\end{array}\right]$$

$$T[1]=AT[0]+Q[0]$$

$$T[k]=\underset{\text{step 0}}{\underbrace{A^{k}T[0]}}+\underset{\text{steps from 1 to k-1}}{\underbrace{\sum _{m=1,n=2}^{m=k-1,n=k}{A}^{m}Q[k-n]}}+\underset{\text{step k}}{\underbrace{Q[k-1]}}$$

$$p_{ij}^k=\{\begin{array}{ll}{p}_{vf}^k,& \text{if }i=even,j=odd,j-i=1\\ p_{vv}^k,& \text{if }i=even,j=even,j-i=2\\ p_{ff}^k,& \text{if }i=odd,j=odd,j-i=2\\ p_{fv}^k,& \text{if }i=odd,j=even,j-i=3\\ 1,& \text{if }i\ge 2n-1\text{ and }j=2n+1\\ 0,& \text{otherwise}\end{array}$$

$$\rho (s_i^k)=\{\begin{array}{ll}{\text{A}}^{\theta }T[0],& \text{if }{s}_i^k\in s_{init}\\ {\text{A}}^{\theta -k}{\mathcal{Q}}_i^k,& \text{if }\theta ,K>k>0\\ J_i^k,& \text{if }\theta =K\\ 0,& \text{otherwise}\text{par}\end{array}$$

$${\xi }_m(s_i^k)=\{\begin{array}{ll}\frac{Q_{int}}{C_m},& \text{if }occupied\in L(s_i^{k-1})\\ 0,& \text{otherwise}\end{array}$$

$${\zeta }_m(s_i^k)=\{\begin{array}{ll}\frac{Q_{rad}}{C_m},& \text{if }heatin{g}_{on}\in L(s_i^{k-1})\\ 0,& \text{otherwise}\end{array}$$

$$P({\mathcal{Z}}^{k+1}=\mathcal{X}|{\mathcal{Z}}^k=f)=\{\begin{array}{ll}\frac{{\sum }_{i=0}^{\mathcal{W}}{f}_i^k\hookrightarrow f_i^{k+1}}{{\sum }_{i=0}^{\mathcal{W}}{f}_i^k\hookrightarrow f_i^{k+1}+{\sum }_{i=0}^{\mathcal{W}}{f}_i^k\hookrightarrow v_i^{k+1}},& \text{if }\mathcal{X}=f\\ 1-P({\mathcal{Z}}^{k+1}=f|{\mathcal{Z}}^k=f),& \text{if }\mathcal{X}=v\end{array}$$

$$P({\mathcal{Z}}^{k+1}=\mathcal{X}|{\mathcal{Z}}^k=v)=\{\begin{array}{ll}\frac{{\sum }_{i=0}^{\mathcal{W}}{v}_i^k\hookrightarrow f_i^{k+1}}{{\sum }_{i=0}^{\mathcal{W}}{v}_i^k\hookrightarrow f_i^{k+1}+{\sum }_{i=0}^{\mathcal{W}}{v}_i^k\hookrightarrow v_i^{k+1}},& \text{if }\mathcal{X}=f\\ 1-P({\mathcal{Z}}^{k+1}=f|{\mathcal{Z}}^k=v),& \text{if }\mathcal{X}=v\end{array}$$

Full text

1 Introduction

The Internet-of-things (IoT) has enabled the realization of Smart Buildings, which aim to deliver useful building services that are efficient, cost effective, reliable, an ubiquitous ensuring occupants’ comfort. Smart building automation systems aim to reduce energy consumption of buildings by operating the Heating, Ventilation and Air Conditioning (HVAC) system in a more efficient way, thereby to addressing energy and environmental concerns [16]. The first step to achieve optimal performance is to construct a quantitative modeling framework of the building itself from energy consumption perspective. This is a difficult task due to the inherent complexity of all the different thermal interactions in the rooms, together with all the unmeasured disturbances as a result of changes in weather or occupancy [25]. Furthermore, the cost of energy adds more complexity to the problem. For instance, in the winter season, given a building topology, occupancy pattern and other factors (e.g., solar heat gain), what should be an effective heating strategy while considering the heating cost and occupants’ comfort? An intelligent heating strategy can be to predict the expected temperature of a zone ahead of time based on the expected occupancy and underlying thermal model of the zone, and turn the heating on/off accordingly. This also allows to schedule heating when energy prices are lower, hence reducing the overall building operating cost. Heating a room when it is empty is of course a counter-intuitive choice in many cases, however, if preheating a room is cheaper (based on the electricity cost) and capable of keeping the room warm for next few hours, then suddenly such a strategy can become a preferred option. To assess such strategies, a novel building thermal model and analysis framework is thus required, which motivates our work in this paper.

Related works. Different modelling techniques for building energy management have been proposed in the literature  [29, 30, 36], the most popular of which is to use a Resistance Capacitance (RC) network to describe the underlying system process dynamics of buildings and zones  [6, 38]. In [19], a seventh-order RC network is used to model two zones followed by a model order reduction technique based on the underlying physical quantities. Using this order reduction methodology, a library of thermal models with varying orders of complexity has further been presented. In [17], the authors proposed a second-order RC-network for modeling the heat conduction through walls, while [12] adopted the RC-network framework and then proposed a technique for aggregation-based model reduction for large buildings. The heat gained within a zone in a building is highly dependent on the occupancy of the zone. Consequently, several advanced models have been proposed in the literature to model building occupancy [37, 7]. In [15], proposed a two-tier deconvolutional approach is proposed to model occupancy profile of a room by using pilot data to identify models for the CO2 and temperature from which occupancy is estimated. The occupancy estimation problem is modeled in [22] as a regression framework. Alternative techniques based on Markov principles have been presented in [31, 28]. The authors in [7] make use of inhomogenous Markov chains to detect the presence of occupants where the states represented the increment of occupancy in a zone. In the work [14], Gaussian mixture models were first used to identify the sensor which yielded the highest information gain for detecting the number of occupants in the zones. The selected features are then fed into a hidden Markov model to estimate the number of occupants. Since occupancy also stimulates the demand of other heat sources, such as lighting or electronic devices, occupancy has become an important factor in estimating building energy consumption [35, 27]. To summarize, researchers from academia and industry has made substantial efforts to in model building thermal dynamics and occupancy, but they have mostly dealt with them as two separate problems instead of an unified model that captures both the thermal dynamics and stochastic occupancy pattern.

Contributions of this paper. In this paper, we aim to bridge this gap and propose a novel framework, and develop an unconventional model for analyzing the temperature evolution of zones in a building. The proposed model incorporates the thermal dynamics, stochastic occupancy pattern, and heating sources (e.g., radiators) in a zone. More specifically, the thermal dynamics of each zone is modeled using the discrete time-inhomogeneous Markov reward model (DTIMRM) formalism from the given building topology as RC equations 111An RC network can be generated from the building information models (BIM), or by using the state-of-the-art tools such as EnergyPlus [10], occupancy pattern of a zone, and the radiator status from the given heating strategy. The DTIMRM states represent the occupancy pattern (calculated from a given dataset) of the zone (occupied or empty) and the thermal dynamics of the zone including the given heating strategy is captured as state rewards. The time-inhomogeneity of the model is then resolved by unfolding it over a specified number of time steps, thus converting an DTIMRM in a discrete time-homogeneous Markov reward model (DTMRM). Afterwards, we perform parallel composition of these DTMRMs, each of which represent a zone in a building.

The composed model is then evaluated to quantify the expected temperature of each zone with respect to different heating control strategies. The heating strategies are preset based on the predefined thermal conditions. This falls under the framework of rule based control which is commonly employed in the heating of buildings [1]. Performing a heating strategy analysis allows us to highlight the performance of the developed framework under different settings, while demonstrating the ability of the framework to guide us in selecting the correct control strategy depending on the specific needs. Using our proposed framework, computation of the probabilistic reachability does not require a large state-space, as witnessed when using probabilistic reachability analysis tools such as FAUST2 [32] which requires gridding over the individual zone temperature range. Such quantitative analysis can help to adopt the most effective and cost efficient heating strategy that maintains a specified range of temperature (for occupants’ comfort) during the intended duration of operation as shown in our obtained results. A key advantage of our proposed model is its capability to deal with the time-inhomogeneity, that none of the above mentioned existing approaches can handle. The framework is also generic enough and can be easily extended to model additional number of zones. As mentioned earlier, the building topology is modeled as RC equations (treated as inputs) which are converted to Markov models later. Hence, different building typologies can be evaluated using our approach with minimal effort. This is also true for evaluating different occupancy patterns and heating strategies.

For automated analysis of the proposed models, we use the Probabilistic Model Checking (PMC) technique [24], and more specifically, the PRISM model checker [23]. PMC is mainly based on the construction and analysis of a probabilistic model, typically a Markov chain or a Markov process. These models are constructed in an exhaustive fashion. In discrete-event simulations, approximate results are generated by averaging results from a large number of random samples. In contrast, PMC applies numerical computations to provide accurate results which makes our approach more reliable. Although the approach in [12] also used a Markov chain for modeling thermal dynamics of a building, where temperature of zones is represented as states, such a model cannot incorporate thermal gains from internal or external sources. In other words, it can only model the conduction between the zones via the surface between them (e.g., walls). Our proposed framework overcomes this limitation by using the Markov reward model formalism. In addition, due to the use of formal verification tools, the results obtained using our framework can aid in certifying safety-critical building facilities, such as hospitals, medical labs, military and defense-related infrastructures, nuclear power plants, food processing plants and so on.

Paper organization. The outline of the paper is as follows. Section 2 presents some background about probabilistic model checking. In Section 3, we describe the framework in detail while Section 4 outlines the obtained results. In Section 5, we discuss some important aspects of our proposed framework approach, and finally, Section 6 concludes the work with directions of future research.

2 Preliminaries

2.1 Probabilistic model checking

Model checking [9] is a well-established technique used in both industry and academia to verify the correctness of finite-state concurrent models. Probabilistic Model Checking (PMC) [24] deals with systems that exhibit stochastic behavior, and is based on the construction and analysis of a probabilistic model of the system, typically a Markov chain. We make use of the DTMRM [33] formalism based on a discrete time Markov chain (DTMC) for modelling the occupancy and thermal dynamics.

Definition 1 (Discrete time Markov chain (DTMC))

Let AP be a fixed, finite set of atomic propositions. A (labelled) DTMC $\mathcal{D}$ is a tuple $(S,s_{init},TL,P,L)$ where S is a finite set of states, $s_{init}\in S$ is the initial state, $TL$ is a set of transition labels denoting the transitions between states, $P:S\times S\rightarrow[0,1]$ is a probability matrix such that $\sum_{s^{\prime}\in S}p(s,s^{\prime})=1$ for all $s\in S$, and $L:S\rightarrow 2^{AP}$ is a labelling function that assigns to each state $s\in S$ the set $L(s)$ of atomic propositions that are valid in $s$. The time-inhomogeneous version of DTMC $\mathcal{D}$ (denoted as $\mathcal{D}_{ih}$) has the probability matrix as a function of time, e.g., $P^{k}:S\times S\rightarrow[0,1]$ given $k\in\mathbb{N}$.

Definition 2 (Discrete time Markov reward model (DTMRM))

A discrete-time reward model (DTMRM) $\mathcal{M}$ is a pair $\mathcal{(D,\rho)}$ where DTMC $\mathcal{D}=(S,s_{init},TL,P,L)$ and $\rho:S\rightarrow\mathbb{R}_{\geq 0}$ is a reward assignment function that associates a real reward to any state in $S$. A Real number $\rho(s)$ denotes the state reward $s$. If DTMC $\mathcal{D}$ in the DTMRM pair is a time-inhomogeneous DTMC ($\mathcal{D}_{ih}$), then such a model is considered as a discrete time-inhomogeneous reward model (DTIMRM).

Path and cumulative rewards. Let $\mathcal{M}$ be a DTMRM with underlying DTMC $\mathcal{D}=(S,s_{init},TL,P,L)$ and reward function $\rho$. A finite path $\sigma$ with length $l$ is a sequence $s_{0}\rightarrow\dots\rightarrow s_{l-1}$ with $s_{i}\in S$ and $p(s_{i},s_{i+1})>0$ for $0\leq i<l-1$. Let $Pr_{s}$ denote the unique probability measure as defined in [3] on the set of paths that start in state $s$. The cumulative reward along a finite path $\sigma$ with length $l$ is defined as $\rho(\sigma)=\sum_{i\geq 0}^{l-1}\rho(s_{i})$. In the Probabilistic model checking approach using DTMCs or DTMRMs, the properties are usually expressed in some form of temporal logic such as Probabilistic Computation Tree Logic (PCTL) [20], a stochastic variant of the well-known Computation Tree Logic (CTL) [8].

2.2 PRISM modeling and property specification

PRISM [23] is a well known probabilistic model checking tool for the formal modeling and verification of stochastic systems. The current version of the tool supports four types of probabilistic models: discrete-time Markov chains (DTMCs), continuous-time Markov chains (CTMCs), discrete-time Markov decision processes (MDPs) and probabilistic timed automata (PTA). Models in PRISM are specified using a high-level modelling language based on the Reactive Modules formalism [2].

To analyze a PRISM model, it is required to specify one or more properties. Since PRISM is a probabilistic model checker, the property specification language in PRISM is based on the temporal logic and includes operators from PCTL, CSL [5] and its extensions[21]. An example of a PCTL property in PRISM with its natural language translation is: $\mathcal{\textbf{P}}_{\geq 0.95}[\textbf{F}~{}complete]$ - “The probability of the system eventually completing its execution successfully is at least 0.95”. Note, the F operator in PRISM property specification language refers to the eventually ($\Diamond$) operator in PCTL [24]. Additional properties can be specified by adding the notion of rewards. Each state (and/or transition) of the model is assigned a real-valued reward, allowing queries such as: $\mathcal{\textbf{R}}_{=?}[\textbf{F}~{}success]$ - “What is the expected reward accumulated before the system successfully terminates?” The R operator in PRISM property specification language extends the PCTL with the notion of rewards [24]. Rewards can be used to specify a wide range of measures of interest, for example, the total operational costs and the total percentage of time the system is available. In this paper, we utilize the notion of cumulative rewards (using the C operator) to evaluate the temperature of a zone in a building. A detail description of the PRISM language, property specification and probabilistic model checking using the PRISM tool can be found in [24].

3 Thermal model of a Building

A high-level description of the proposed modeling framework is depicted in Figure 1. Given a building and its corresponding topology, we first construct the equivalent RC-network model made up of $N$ zones. Note, the zones should be connected to each other through walls. This RC-network model is converted into a state-space models which allows for ease of conversion into a Markov model. We also perform an occupancy pattern calculation which computes the probability of occupancy in each zone during each time step. For each zone, we construct an equivalent DTIMRM. Afterwards, we unfold them to a specified number of time steps which converts them into equivalent DTMRMs. These models incorporate the calculated occupancy pattern and the heating strategy to evaluate. We further perform the *parallel composition *of these DTMRMs to obtain the whole building Markov model. Quantitative analysis is then performed on this model through the use of PMC technique (using a probabilistic model checker). For this paper, we consider two teaching rooms (zones) within the Department of Computer Science, University of Oxford and restrict ourselves to having two heat sources in each zone: (i) radiators 222Use of radiators for heating rooms are still common in many European universities and offices, including the University of Oxford. However, this does not restrict the proposed framework to be applied to the forced-air heating systems, which is very common in the US. This will be further clarified in Remark 2. and (ii) occupants. The teaching rooms are connected back to back through walls. The rooms are used between 8 am and 5 pm during term time. The stochastic occupancy pattern of these zones are calculated from our developed dataset containing the past occupancy record of a zone. A layout of these rooms is shown in Figure 2(a) and Figure 2(b).

3.1 RC-network representation

A building thermal model is constructed by combining parameters of thermal interaction between two zones separated by a solid surface (e.g., walls, windows, ceilings, and floors). For the sake of simplicity, here we ignore the inter-zone convective heat transfer that occurs through the open doors and hallways. A lumped parameter model of combined heat flow across a surface is modeled as a simple RC-network, with current and voltage being analogous to heat flow and temperature, respectively. In this modeling framework, the capacitance are used to model the total thermal capacity of the wall, and the resistances are used to represent the total resistance that the wall offers to the flow of heat from one side to the other.

A building topology can be represented as a directed graph $\mathcal{G}=(\mathcal{V},\mathcal{E})$, where $\mathcal{V}:=\{1,...,N\}$ denotes the set of $N$ nodes of the graph. A node represents a physical zone (e.g., a room or a hallway). Each node $i\in\mathcal{V}$ corresponds to a temperature $T_{i}$ and has an associated capacitance $C_{i}$. The set $\mathcal{E}\subset\mathcal{V}\times\mathcal{V}$ denotes the set of all edges which represent pathways for conductive heat transports. More specifically, each edge $(i,j)\in\mathcal{E}$ represents the conductive thermal interaction between the nodes $i$ and $j$ represented by a resistance $R_{i}\in\mathbb{R}_{+}$. This is a different approach compared to [13], where the building topology was represented using an undirected graph and the authors considered a point on the surface (such as the wall) as a physical zone, whereas we ignore that. The thermal dynamics of a multi-zone building, described by a directed graph $\mathcal{G}$, is modeled by a system of first order coupled linear differential equations as follows, for $i=1,...,N$:

[TABLE]

where $\mathcal{N}_{i}:=\{j\in\mathcal{V}:j\neq i,(i,j)\in\mathcal{E}\}$ denotes the set of neighbors connecting to the node $i$. The heat gain term ${Q}_{i}(t)$ is the rate of thermal energy entering the node $i$ from all sources, such as occupants, lighting etc., and radiators other than the conduction from neighboring nodes. For instance, let us consider two rooms connected back to back to one another and are depicted using Figure 2(c) following [26], where $V_{ref}$ is a reference temperature. The thermal model is composed of two zones having an average zone temperature $T_{1},\;T_{2}$, which exchange heat with one another assuming that the heat is only gained in a zone via its occupants and radiators. To build the equivalent RC network the following assumptions have been applied:

Similar thermal properties between heat zones. 2. 2.

Uniform temperature across the walls. 3. 3.

Negligible solar gain and heat exchange with outside through walls. (Note that this is not a limitation of the proposed modeling approach. How to incorporate these factors in the proposed model has been addressed in Section 5.) 4. 4.

The air in each zone having one temperature across its volume.

Using Equation (1) we obtain the equations for each zone:

[TABLE]

where, ${i\in\{1,2\}}$ represents the respective zone and ${Q}_{i}(t)={Q}_{rad_{i}}(t)+{Q}_{int_{i}}(t)$, in which stochastic terms ${Q}_{int_{i}}(t)$ and ${Q}_{rad_{i}}(t)$ represent the internal heat gain (i.e., the rate of heat generated by the occupants, equipment, lights, etc.) and heat gains from the radiators respectively, and can be specified as follows:

[TABLE]

[TABLE]

Equation (2) can be rewrittten in the form of a linear time invariant state space model as follows:

[TABLE]

where the state vector ${T}:=[{T}_{1},...,{T}_{N}]^{T}$ denotes the temperature of each zone and $U$ denotes the heat gain vector ${Q}:=[{Q}_{1},...,{Q}_{N}]^{T}$. Both the transition rate matrices $\hat{A}$ and $\hat{B}$ have dimensions of $N\times N$. The entries of $\hat{A}$ matrix are given by:

[TABLE]

and the entries of matrix $\hat{B}$ can be expressed as:

[TABLE]

For instance, the example shown in Figure 2 (c) with two zones, can be expressed as:

[TABLE]

Based on the data obtained from [19] ($C_{1}=1.3700,C_{2}=1.000,R_{1}=1.7429,R_{2}=5.5897$), we discretize the model using the forward Euler method [4] with a sampling time of 1 hour ($\Delta=1~{}hr$), and obtain the A and B matrices (discretized version of $\hat{A}$  and  $\hat{B}$ matrix, respectively)

[TABLE]

Remark 1

We are interested in modeling a university building with classrooms (zones), where the duration of each class is around 1 hour. Moreover, as mentioned earlier, the thermal gain of a zone is assumed only from the occupants or from the radiator. Since thermal dynamics of zones are relatively slow with time constants in tens of minutes [17] to hours [34], the dynamics of the radiators are replaced by static gains in this paper without significant loss of accuracy. This also justifies the time interval of one hour while discretizing from continuous-time to discrete-time.

3.2 Markov modeling of occupancy and thermal dynamics

In this section, we show how the thermal dynamics of zones can be encoded using Markov chains. For such modeling, we use the discretized A and B matrices as in the previous section. The equation for the expected temperature of zones at time step 1 can be written as :

[TABLE]

where, ${Q}[0]$ represents the combined gain from the occupants and radiator from the previous time step (time step 0, e.g. initial values in this case). Similarly, the expected temperature for the zones in $k$-th time step can be derived as:

[TABLE]

Since occupancy is a stochastic term, e.g., a zone can be either occupied or empty at a given time step, the Markov chain analogy is now clear. Before proceeding further, it is worth mentioning the assumptions of the proposed model as follows:

• •

A zone can be either occupied or empty. If occupied, it is assumed to be fully occupied, otherwise empty.

• •

Since we are modeling classrooms, each time step (class duration) is assumed to be equal to 1 hour.

• •

The internal gain (per hour) of a zone from the occupants and zone’s radiator is assumed to be constant.

• •

Each zone has one radiator for heating purpose.

Remark 2

In contrast to the use of radiators for heating, the forced-air heating system is more popular in the US. In a forced-air system, the heating rate is modeled as a continuous variable, however, it can be abstracted to a single number. In order to do that, one needs to discretize the equation $\Delta H=C_{pa}\dot{m}_{i}^{in}(t)(T^{s}-T_{i}(t))$, where $C_{pa}$ is the specific heat capacitance of the supplied air at constant pressure, $T^{s}$ is the supply air temperature (usually constant for a variable air volume system [11]), and $\dot{m}_{i}^{in}$ denotes the mass flow rate of the supply air entering the $i$-th zone.

Given that a building $\mathcal{B}$ has $N$ zones, each zone can be modeled as a separate discrete-time inhomogeneous Markov reward models (DTIMRM) as a compact representation. In the DTIMRM model, the $states$ represent the occupancy of a zone (either occupied or empty) and we use the notion of $state~{}rewards$ to encode the thermal dynamics in the proposed Markov model based on Equation (5). The rewards are collected at each time step based on the current state, and hence the expected temperature at $k$-th time step can be calculated as cumulative rewards up to time step $k$. The DTIMRM of a zone is shown in Figure 4. For clarity, the figure is color coded - the dark green node represents an empty zones and light yellow colored node represents an occupied zone. We assume that the zone was initially empty encoded by state 0.

Definition 3

(zone-DTIMRM) The occupancy and temperature evaluation of a zone can be defined as discrete-time inhomogeneous Markov reward model (DTIMRM) $\mathcal{R}_{ih}$ which is a pair $(\mathcal{D}_{ih},\rho^{k})$ with time-inhomogeneous DTMC $\mathcal{D}_{ih}=(S,s_{init},L,P^{k})$ and a reward function $\rho^{k}$, where:

• •

$S=\{s_{1},s_{2}\}$, represents if a zone is either occupied or empty;

• •

$s_{init}=\{s_{1}\}$ is the initial state;

• •

$L:S\rightarrow 2^{AP}$, a labeling function that assigns each state $s_{i}\in S$ a label from the set of atomic propositions $AP=\{occupied,empty,heating_{on},heating_{off}\}$;

• •

$P^{k}:S\times S\rightarrow[0,1]$ is the time dependent stochastic (occupancy) matrix that associates each pair $(s,s^{\prime})\in S\times S$ the transition probability where $\sum_{s^{\prime}\in S}p^{k}(s,s^{\prime})=1$ for all $s\in S$. More specifically: $P^{k}=\begin{bmatrix}p_{ff}^{k}&p_{fv}^{k}\\ p_{vf}^{k}&p_{vv}^{k}\\ \end{bmatrix}$ where $f$ and $v$ stands for occupied and empty states respectively, for instance, $p_{vf}^{k}=p^{k}(v,f)$.

• •

$\rho^{k}:S\rightarrow\mathbb{R}_{\geq 0}$, a time dependent reward function that assigns a reward to each state to encode the thermal dynamics and gains (from heating and occupant) to the states at a corresponding time step $k$ (explained in detail in Definition 4).

To resolve the time-inhomogeneity, we expand the model by unrolling the DTIMRM up to $K+1$ steps ($k\in K$), where $K$ denotes the time window (operating time) of the zones. Thus the DTIMRM is converted into a discrete-time Markov reward model (DTMRM). Figure 3 shows the Markov model of a zone for $K+1$ steps. From each empty state, there are two possible transitions labelled as $t_{k}$ with probability $P_{vv}^{k}$ and $P_{vf}^{k}$, which respectively represent the probability of being empty or occupied in the next time step (how we derive these transition probability is explained in the next subsection). Similarly, from each occupied state, there are two possible transitions labelled as $t_{k}$ with probability $P_{fv}^{k}$ and $P_{vv}^{k}$ that respectively denote the probability of being empty or occupied in the next time step. The grey colored node ($K+1$-th state) denotes the sink state. Note that, we use the same labels for labelling transitions at a specific time step, however only one of them will be triggered at a given time step. Formally, a zone can be redefined as follows:

Definition 4

(zone-DTMRM) The occupancy and temperature evolution of a zone can be defined as discrete-time Markov reward model (DTMRM) $\mathcal{R}$ which is a pair $(\mathcal{D},\rho)$ with time-homogeneous DTMC $\mathcal{D}=(S,s_{init},TL,L,P)$ and a reward function $\rho:S\rightarrow\mathbb{R}_{\geq 0}$ that assigns a reward to each state $s_{i}^{k}\in S$, where:

• •

$S=\{s_{0}^{0},s_{1}^{1},s_{2}^{1},s_{3}^{2},s_{4}^{2},....,s_{2n}^{K},s_{2n+1}^{K+1},\}$, where $s_{i}^{k}\in S$ represents $i^{th}$ state in $k^{th}$ time step and $k={1,2,...,K+1}$;

• •

$s_{init}=\{s_{0}^{0}\}$ is the initial state;

• •

$TL=\{t_{1},t_{2},...,t_{K+1}\}$ is a set of transition labels denoting the time steps and used for synchronization (explained in Definition 5);

• •

$L:S\rightarrow 2^{AP}$ is a labeling function that assigns each state $s_{i}^{k}\in S$ a label from the set of atomic propositions $AP=\{occupied,empty,heating_{on},heating_{off}\}$;

• •

$P:S\times S\rightarrow[0,1]$ is the stochastic (occupancy) matrix that associates each pair $(s,s^{\prime})\in S\times S$ the transition probability where $\sum_{s^{\prime}\in S}p(s,s^{\prime})=1$ for all $s\in S$. Since, in our model each odd numbered states denote an occupied zone and each even numbered states denote an empty state, for each time step $k$ the entries of the stochastic matrix $P=(p_{ij}^{k})$ are given by:

[TABLE]

• •

$\rho$ is a reward function defined as follows:

[TABLE]

where $\textbf{A}^{\theta}=[X_{mj},X_{m(j+1)},...,X_{mN}]$ is $1\times N$ matrix for the $m$-th zone given that $X=A^{\theta}$. The term $\theta$ defines the specific time of interest to evaluate the model (explained in Example 1) and the discretized $A$ matrix is derived from thermal dynamics as described in the previous subsection. For instance, if a building has two zones, then $\textbf{A}^{\theta}$ can be defined as $\textbf{A}^{\theta}=[X_{11},X_{12}]$ for zone 1 ($m=1$) and $\textbf{A}^{\theta}=[X_{21},X_{22}]$ for zone 2 ($m=2$). $T[0]$ is a $N\times 1$ matrix specified as $T[0]=[T_{1}[0],T_{2}[0],...,T_{N}[0]]^{T}$, where $T_{m}[0]$ represents the initial temperature of $m$-th zones. $\mathcal{Q}^{k}_{i}$ is also a $N\times 1$ matrix specified as $\mathcal{Q}^{k}_{i}=[\mathcal{C}^{k}_{1i},\mathcal{C}^{k}_{2i},...,\mathcal{C}^{k}_{Ni}]^{T}$, where $\mathcal{C}^{k}_{mi}$ represents the combined (coupled) gain of the $m$-th zone at time step $k$ for the $i$-th state. $\mathcal{C}^{k}_{mi}$ can be expressed as $\mathcal{C}^{k}_{mi}=\{\xi_{m}(s_{i}^{k})+\zeta_{m}(s_{i}^{k})\}$, and ${J}^{k}_{i}=\mathcal{C}^{k=K}_{m(i=2n)}$ for the $m$-th zone. $\xi:S\rightarrow\mathbb{R}_{\geq 0}$ and $\zeta:S\rightarrow\mathbb{R}_{\geq 0}$ are functions that map each state $s_{i}^{k}\in S$ (from the DTMC of $m$-th zone ) as follows:

[TABLE]

and,

[TABLE]

where $Q_{int}$ and $Q_{rad}$ are the per time step temperature gain of an occupied zone and from radiator from previous time step of a respective zone, and $C_{m}$ is the capacitance of the $m$-th zone.

Remark 3

Note that, we use heating as a deterministic input to the model unlike the stochastic one as occupancy. In Section 4 we evaluate different heating strategies by providing the information if the heating is on in that time step, which sets the value of $\zeta_{m}(s_{i}^{k})$ for different time steps $k$ in the model for the $m$-th zone.

Example 1. Let us assume that we build the model for 9 time steps. So, using this model we can reason about the expected temperature at different time steps starting from time step 1 to time step 9. Consider a scenario where we are interested in reasoning about the expected temperature of the zone at time step $\theta=3$ and $\theta=8$ using the model that we built for 9 time steps. Since the thermal dynamics depends on the elapsed time, the state reward model needs to be dynamic as well. Hence, the parameter $\theta$ is used to assign the dynamic reward structure based on a specific time of interest. For the case $\theta=3$, the reward of state 0 is set to $\textbf{A}^{3}T[0]$, for states 1 and 2 the rewards are set to $\textbf{A}^{2}\mathcal{Q}^{1}_{1}$ and $\textbf{A}^{2}\mathcal{Q}^{1}_{2}$, for states 3 and 4 the rewards are set to $\textbf{A}\mathcal{Q}^{2}_{3}$ and $\textbf{A}\mathcal{Q}^{2}_{4}$, and finally, for state 5 and 6, the state rewards are set to $\mathcal{Q}^{3}_{5}$ and $\mathcal{Q}^{3}_{6}$. The reward for the rest of the states should be set to 0 since we are reasoning about the temperature at $\theta=3$, and hence we are interested in the cumulative reward from $k=0$ to $k=3$ while ignoring the rewards for other states. However, if we want to reason about the temperature of the same zone at time step $\theta=8$, the rewards for the same states will be set to different values. For this specific case the rewards from state 0 to state 4 need to be set as $\textbf{A}^{8}T[0]$, $\textbf{A}^{7}\mathcal{Q}^{1}_{1}$, $\textbf{A}^{7}\mathcal{Q}^{1}_{2}$, $\textbf{A}^{6}\mathcal{Q}^{2}_{3}$, $\textbf{A}^{6}\mathcal{Q}^{2}_{4}$ respectively, the rewards for state 5 to state 18 will also be assigned according to the Equation (6). Based on this example we can see that parameter $\theta$ can to be adjusted based on the time step of interest.

Given the DTMRM of each zone $\mathcal{R}_{i}\in\mathcal{B}$, the building model then can be constructed by using the parallel composition. Hence, if the building has $N$ rooms, the building model can be defined as $\mathcal{B}=\mathcal{R}_{1}||\mathcal{R}_{2}||...||\mathcal{R}_{N}$.

Definition 5

(Parallel Composition) Given two DTMRM pairs $\mathcal{R}_{1}=(\mathcal{D}_{1},\rho_{1})$ and $\mathcal{R}_{2}=(\mathcal{D}_{2},\rho_{2})$ with DTMC denoted as $\mathcal{D}_{1}=(S_{1},s_{init2},TL_{1},L_{1},P_{1})$ and $\mathcal{D}_{2}=(S_{2},s_{init2},$ $TL_{2},L_{2},P_{2})$ with reward functions $\rho_{1}$ and $\rho_{2}$, their parallel composition is the system $\mathcal{R}_{1}||\mathcal{R}_{2}=(\tilde{S},\tilde{s}_{0},\tilde{T}L,\tilde{L},\tilde{P})$ with reward function $\tilde{\rho}$ where,

• •

$\tilde{S}=(S_{1}\times S_{2});$

• •

$\tilde{s}_{0}=(s_{init1},s_{init2});$

• •

$\tilde{TL}=TL_{1}\cup TL_{2};$

• •

$\tilde{L}=L_{1}\cup L_{2};$

• •

$\tilde{P}=\begin{cases}\frac{s_{1}\xrightarrow{\alpha,r}s_{1}^{\prime}~{}\text{and}~{}s_{2}\xrightarrow{\alpha,q}s_{2}^{\prime}}{(s_{1},s_{2})\xrightarrow{\alpha,rq}(s_{1}^{\prime},s_{2}^{\prime})},&\mbox{if }TL_{1}=TL_{2}=\alpha\\ \\ \frac{s_{1}\xrightarrow{\alpha_{1},r}s_{1}^{\prime}}{(s_{1},s_{2})\xrightarrow{\alpha_{1},r}(s_{1}^{\prime},s_{2})},\frac{s_{2}\xrightarrow{\alpha_{2},q}s_{2}^{\prime}}{(s_{1},s_{2})\xrightarrow{\alpha_{2},q}(s_{1},s_{2}^{\prime})},&\mbox{otherwise}\\ \end{cases}$

• •

$\tilde{\rho}(s_{i})$ = $\rho_{1}(s_{i})\cup\rho_{2}(s_{i})$ (refer to Remark 4)

where $s_{1},s_{1}^{\prime}\in S_{1}$, $\alpha_{1}\in TL_{1}$, $\alpha_{2}\in TL_{2}$, $\alpha\in TL_{1}\cap TL_{2}$, $r=p_{1}(s_{1},s_{1}^{\prime})$, $s_{2}$,$s_{2}^{\prime}\in S_{2}$, $q=p_{2}(s_{2},s_{2}^{\prime})$.

Remark 4

Note that, we are interested in reasoning about the expected temperature of zones given that a reward function captures the temperature dynamics of a zone at a specific time step. Hence, each state in the composed model of a building $\mathcal{B}$ has $n$ separate reward functions given that the building has $N$ zones in it. This is unlike the approach introduced in the Markov reward automata [18] where the rewards are added in the composed model.

As mentioned in section 3.1, zones in a building are connected via walls between them, and hence their thermal dynamics are also coupled. An alternative approach for modeling building dynamics with occupancy pattern will be to model both the zone explicitly into a single Markov chain (for the two zone case), instead of building two separate Markov chains and composing them together. However, the complexity of developing such models will increase dramatically with the increasing number of zones. Our proposed modeling approach is modular since each zone is modeled separately at first, and composed later. We use PRISM modeling language (the detail syntax of PRISM modeling language, and the detail of PRISM tool can be found in [23]) for modeling, and encode each zone of a building as separate modules. This makes our model easily extendable to $N$ zones using the module renaming feature in PRISM with minimal effort. The corresponding probabilistic model of the building with $N$ zones in it (composed DTMC in our case) is constructed as the parallel composition of modules specified in PRISM language which is useful for the automated analysis of such models.

3.3 Transition matrix for zone occupancy

In this section, we show how the transition matrix $P$ of a zone can be derived from the occupancy dataset. The occupancy dataset of a zone should contain per-hour occupancy record of each zone in the building for a period of time (say, for the last 6 months). This dataset can be updated periodically with recent data to reflect any changes in the occupancy pattern. Before calculating the probability, it is required to process the dataset in a specific format as shown in Table 1, where we enlist (for each hour) how many times a zone was observed to be empty or occupied from either occupied or empty state for a period of time. In our case, the classes in the university rooms start from 9 am, so the rooms are unoccupied at 8 am which clarifies the NA (not applicable) entries in the table. Given that a zone is occupied at time step $k$, the probability that the zone will be either occupied or empty in time step $k+1$ can be calculated from the dataset as follows:

[TABLE]

where $\mathcal{W}$ represents the total number of days in the dataset, $\hookrightarrow$ represents a transition between occupied (denoted as ‘$f$’) and empty zone (denoted as ‘$v$’) on $i$-th day from time step $k$ to $k+1$. Thus, using Equation (7), we obtain the value of $P_{{fv}}^{k}$ and $P_{{ff}}^{k}$ at time step $k$. Similarly, if the room is currently empty ($\mathcal{Z}_{k}=v$), then the probability that the room will be either occupied or empty in the next time step, can be calculated as follows:

[TABLE]

We obtain the value of $P_{{vv}}^{k}$ and $P_{{vf}}^{k}$ at time step $k$ with the help of Equation (8). For instance, Table 1 shows an example of such a table with hourly occupancy data from 9 am to 5 pm for 6 months ($\mathcal{W}=$180 days). Let us consider the first 7 days record from Table 1. Given that a zone is empty at 8 am and using Equations (7) and (8), we can derive the probability that the zone will be occupied at 9 am is 0.43 and the probability that the zone will be empty is 0.57.

4 Quantitative Analysis

In this section we present the quantitative results obtained from the two zone scenario as shown in Section 2. We use the PRISM probabilistic model checker for specification and analysis of our model. Since the developed PRISM model (using the PRISM modelling language) is parametric, it is possible to evaluate the model for different parameter values without any modification of the core model. After the parallel composition, our approach resulted into 132 states (for the two zone scenario). From our occupancy dataset, we evaluate the occupancy pattern as shown in Figure 5 for zone 1 and zone 2 from 8 am to 5 pm. The initial temperature (at 8 am) of zone 1 and zone 2 are initialized to 18°C and 16°C respectively, the internal gain of a zone from occupants is assumed to be 0.7°C per hour and the gain from the radiator (e.g., heating is on) is assumed to be 1.5°C per hour. Both the zones are empty at 8 am and we evaluate the expected temperature of the zones from 9 am to 5 pm. In PRISM, such analysis can be performed by specifying a PCTL property for cumulative rewards as $\mathcal{\textbf{R}}_{=?}[\mathcal{\textbf{C}}\leq\Psi]$, where ($9\leq\theta\leq 17$) and $\Psi=\theta+1$. Given that we intend to maintain the expected temperature of both zones between 20°C to 22°C, we evaluate six potential heating strategies as shown in Table 3. Strategies from S1 to S4 are self explanatory. In strategy S5, we turn on the radiator for heating in every alternating hours, e.g., the heating is on at 9 am, 11 am, 1 pm, and 3 pm. In contrast, in strategy S6 we turn on the heating only for two time steps, at 8 am (pre-heating the zones) and 9 am, and keep it off for the rest of time. As mentioned earlier, such heating strategies are preset based on the predefined thermal conditions. This falls under the framework of rule based control which is commonly employed in the heating of buildings [1].

We plot the obtained results in Figure 6. In Figure 6(a) and Figure 6(b) we observe that if strategy S1 is followed, e.g., heating is off for both the zones, then the expected temperature of these zones gradually increases (due to the internal gain from occupancy), but stays below 20°C till 5 pm. In contrast, if S2 is followed, e.g., heating is turned on for both zones, then the temperature is expected to cross 22°C by 1 pm in zone 1, and by 12 pm in zone 2. Meanwhile, by the end of the day (5 pm), both zones are expected to reach approximately 30°C. While applying strategies S3 and S4 to our model as shown in Figure 6(c) and in Figure 6(d), we also observe that the expected temperature of both zones eventually crosses the maximum value of our intended threshold of 22°C, even though it never reaches 30°C when compared to the heating strategy S2. With strategy S5 in Figure 6(e), the expected temperature of zone 1 and zone 2 crosses 22°C after 2 pm, whereas both the zones stay below 20°C before 12 pm. Therefore, all of these five strategies can cause discomfort to the occupants. Finally, we apply strategy S6 (we turn on heating only for two hours from 8 am) to our model the obtained result is shown in Figure 6(f). Interestingly, we observe that the temperature of both zones reaches almost 20°C by 9 am, and stay below 22°C till the end of the day. Hence, based on the comfort of the occupants, heating strategy 6 best fits our requirement. Similarly, any given heating strategy can be evaluated using the proposed approach without any major changes to the proposed model.

Cost analysis: We also evaluate these six heating strategies with respect to their cost. Let us consider the following electricity tariff as shown in Table 3, in which the hours of a day have been divided in economy, peak and off-peak hours based on the electricity demand. We only show the tariffs for 8 am to 5 pm based on our operating window. Economy tariff is the cheapest since the electricity demand is the lowest in that part of the day. In contrast, during peak hours the electricity demand is higher, hence the price is most expensive as well. Assuming that each zone has a 1 kW radiator installed for heating, the cost of the evaluated strategies are shown in Table 4, where column 1 enlists the strategies, columns 2, 3, and 4 shows how much energy was consumed in total by both zones in those corresponding hours (economy, peak, off-pick) and column 5 shows the total cost incurred by applying the associated heating strategies. The cost calculation process is straightforward: Total cost = kWs consumed in economy hours * price per kW in economy tariff + kWs consumed in off-peak hours * price per kW in off-peak tariff + kWs consumed in peak hours * price per kW in peak tariff. According to the cost analysis, among the evaluated strategies S6 is not only the most appropriate one in terms of maintaining the comfort level, but also the most cost effective one. Compared to S2, S3, S4 and S5, strategy S6 is 13.5, 6.75, 6.75 and 6.5 times more cost efficient respectively.

5 Discussion

As mentioned in the quantitative analysis section, our model for the two zone scenario requires 132 states. In contrast to our approach, the work in [32] requires the state-space of the two zones to be gridded over the possible temperature values the zone temperatures can be in. If we consider the temperature of the zones to range between 15 °C and 30 °C (reasonable temperature range the zones will be in) and assume a uniform gridding procedure with a maximal error of 0.5 °C, it will require $900$ states. Hence, we argue that our approach is more scalable in such cases. The proposed framework (as shown in Figure 1) is generic, which means $N$ number of zones can be modeled using the demonstrated approach, given that the total number of states does not cause a state space explosion of the PRISM tool. PRISM model checker includes multiple model checking engines, many of which are based on symbolic implementations (using binary decision diagrams and their extensions) which enables the probabilistic verification of models of up to $10^{10}$ states.

It is worth mentioning that for the sake of simplicity, in this paper we ignored the local weather prediction outside the building (heat exchange with outside environment through walls). This can be addressed by including the outside environment as a separate zone in our model, and thus modeling the heat exchange with the adjacent zones. Note that, our modeling can express both heating and cooling of zones, where the heating requires to be encoded as a positive reward, and cooling requires to be encoded as a negative reward in the DTMRM of a zone. The PRISM tool allows modeling and analysis of non-negative rewards/costs only, and hence there is no straightforward way to analyze the model with cooling using the PRISM tool. However, one way to overcome this limitation is to encode the cooling also as a positive reward, and calculate the cumulative reward until the target state is reached. Finally, the cumulative reward needs to be deducted from the initial temperature of a zone.

6 Conclusions

In this paper we presented a framework for modeling and automated analysis of thermal dynamics of a building using the discrete-time Markov reward model formalism. The model is built using the stochastic occupancy pattern, building’s thermal properties and the heating strategies to be evaluated. We evaluated six different heating strategies using the PRISM model checker and concluded that the “selective strategy” is the best option for maintaining the occupants’ comfort. Also, we performed a cost analysis that shows that the “selective strategy” is the most cost-efficient option among the evaluated heating strategies. As a part of our future works, we plan to incorporate solar gain and heat exchange with outside environment in our model. We also aim to extend the framework for automated synthesis of heating strategies while minimizing the cost and maximizing the occupants’ comfort.