Computational Intelligence Technologies for Occupancy Estimation and Comfort Control in Buildings

Abstract

This paper focuses on the development of a multi agent control system (MACS), combined with a stochastic based approach for occupancy estimation. The control framework aims to maintain the comfort levels of a building in high levels and reduce the overall energy consumption. Three independent agents, each dedicated to the thermal comfort, the visual comfort, and the indoor air quality, are deployed. A stochastic model describing the CO₂ concentration has been studied, focused on the occupancy estimation problem. A probabilistic approach, as well as an evolutionary algorithm, are used to provide insights on the stochastic model. Moreover, in order to induce uncertainty, parameters are treated in a fuzzy modelling framework and the results on the occupancy estimation are investigated. In the control framework, to cope with the continuous state-action space, the three agents utilise Fuzzy $\mathcal{Q}$-learning. Simulation results highlight the precision of parameter and occupancy estimation, as well as the high capabilities of the control framework, when taking into account the occupancy state, as energy consumption is reduced by $55.9\%$, while the overall comfort index is kept in high levels, with values close to one.

1. Introduction

The problem of occupancy estimation and prediction, in the context of comfort control in buildings, is of crucial importance and has been studied by many authors [1]. Our motivation for this work is to combine an algorithm that estimates the number of occupants inside a room, given the CO₂ state concentration, with a multi agent control system that utilises the available information and controls the comfort indexes: thermal comfort, visual comfort, and indoor air quality. Most heating, ventilation, and air conditioning systems (HVAC) function on the assumption of maximum occupancy, not taking into account the variations of the occupancy level. There are various factors to be studied, in order for an estimation algorithm to be developed, that computes the actual number of occupants that entered a room during a time period.

In [2], Dong et al. study the control of HVAC systems using three approaches for the occupancy prediction: the expectation maximisation algorithm, a finite state automata method and introduces a method based on uncertain basis functions. A method based on a hybrid neural network is proposed in [3], where an extreme learning machine tuned with differential evolution is proposed. Models of occupancy in a single person office are studied in [4] using a non-homogeneous Poisson process to simulate the occupancy sequence. In [5] Feng et al. study the simulation of occupancy in buildings in four levels: the number of occupants in a building, the occupancy status in a space, the number of occupants in a space and the space location of an occupant. The correlation between occupancy and energy performance of a building is presented in [6], together with the impact on the quality of the indoor environment. A study that takes into account the stochastic nature of occupants behaviour can be found in [7], where a Markov chain is proposed. The Markov process is also studied in the context of occupancy models in [8]. Occupancy prediction based on the carbon dioxide concentration is studied in [9], where Jiang et al. study a Bayesian filtering approach. In [10], authors study models for CO₂ concentration prediction: A deterministic mass balance differential equation is firstly used and, moreover, a linear model is constructed using system identification techniques. The importance of occupancy for building control is highlighted in [11], where a model predictive controller utilises the occupancy information.

Occupancy estimation models can be classified into three categories. (a) statistical (black-box) models, (b) physical/deterministic (white-box) models, and (c) grey-box models. The first category includes hidden Markov models, machine learning methods, and artificial neural networks [2,3,4]. The second category, the physical models, are most often based on a CO₂ mass balance equation [9,12,13]. The benefits of the first and second category, i.e., black and white box models, are incorporated in the category of grey-box models. The latter models synthesise physical knowledge with statistical and data-driven parameter estimation. A grey-box model can describe variations in a room’s carbon dioxide level and estimate room occupancy [3,14,15,16,17].

Multi agent control systems are a group of agents which, in order to solve a problem, interact with the environment and with each other. MACS have been proposed for achieving high levels of thermal comfort inside a building and reducing the consumption of the heating/cooling system by efficiently controlling the indoor temperature. In [18], a MACS for a building with two rooms is proposed in order to control the thermal comfort of the two rooms through a central HVAC system. Intelligent control of HVAC systems [19] is also studied in [20] using a fuzzy based approach for personal comfort satisfaction. In [21] a MACS with multiple layers and a central agent for coordination is proposed. The system aims to reduce the overall efficiency of the building and simultaneously to satisfy the thermal needs of the occupants. Additionally, MACS have been proposed for achieving air quality, and thermal and visual comfort in the interior of the buildings with minimisation of energy consumption [22]. The MACS controls the actuators regarding the CO₂ concentration, the illumination, and the temperature in conjunction with reducing the energy consumption. These approaches are based on multi-layer MACS with local agents and supervisors which can lead to failure if the central agents fail. Moreover, they are based on the expert knowledge or on offline optimisation algorithms for adjusting parameters and they do not use any learning mechanism [23]. Only recently, few methods using reinforcement learning have been proposed but these methods are only limited to achieve thermal comfort [24]. The main drawback of this method is that it cannot be applied for continuous state-action space. Fuzzy logic systems can be used in order to overcome this drawback. Moreover, none of the above mentioned papers, combine occupancy estimation models with a multi agent control system, in order to reduce energy waste.

In the proposed system, a grey-box model consists of stochastic differential equation (SDE) and fuzzy sets is able to quantify uncertainties, which improves the occupancy estimation. The contribution of this study is to:

• Provide a fully decentralised control architecture;
• Provide a control architecture based on reinforcement learning;
• Incorporate a fuzzy logic system with reinforcement learning;
• Combine a multi agent control system with a grey-box modelling;
• Combine a multi agent control system with occupancy estimation;
• Combine fuzzy modelling with occupancy estimation.

The paper is organised as follows: The problem under study is stated in Section 2, in addition to a discussion on general notions of stochastic processes and the analysis of the CO₂ adopted stochastic model. Stochastic differential equations’ parameter estimation and differential evolution are also analysed. The use of the occupancy estimation algorithm, proposed in [14] is being studied for the estimated and the fuzzy parameters of the stochastic model, using two different occupancy profiles. Furthermore, a thorough analysis of Fuzzy $\mathcal{Q}$-learning and the MACS is presented, together with the mathematical modelling of the building under study. Simulation results are demonstrated in Section 3 for each of the methods used and conclusions are discussed.

2. Materials and Methods

2.1. Problem Statement

The problem studied in the present work focuses on the development of a multi agent control system for the control of the artificial lighting, ventilation, and heating/cooling subsystems, in order for the comfort indexes to be maintained in the desired levels. The MACS should utilise the occupancy information, estimated using carbon dioxide concentration data. The carbon dioxide generation is considered to be given by the solution of a stochastic differential equation. By numerically solving the stochastic differential equation the training data are generated, representing the CO₂ levels in the room under study. Following that, an occupancy estimation algorithm is implemented to give information on the number of occupants. Two main assumptions have been made:

• The training set consists of 1440 state measurements, corresponding to one working day, 24 h, times 60 min per hour;
• The state, i.e., the CO₂ concentration level, is fully observable.

The work flow of the proposed method, graphically illustrated in Figure 1, is the following: Extending the CO₂ mass balance equation to a stochastic differential equation. The stochastic mass balance equation is numerically solved, to generate data of CO₂ concentration. Assuming that the solution time-series represents a stochastic process with unknown parameters, the parameters are estimated by maximising the likelihood function of the data with a differential evolution algorithm. Given the data and the estimated parameters, the occupancy estimation algorithm is implemented to compute the unknown number of occupants, that generated the data. Since the method of maximum likelihood provides only point estimates of the unknown parameters, fuzzy modelling is employed to study the impact of the parameters’ uncertainty on the occupancy estimation. The estimated profile of the building occupancy is going to be used as input in the MACS system. The MACS system consists of three decentralised local agents and controls the actuators, each responsible for the indoor temperature, the indoor horizontal illuminance, and the carbon dioxide concentration. Occupancy estimation defines whether the building is unoccupied and whether the actuators need not to operate, in order to reduce energy waste in time periods where there are no occupants.

2.2. Stochastic Processes and CO₂

The generation and decay of CO₂ is governed by the mass balance ordinary differential equation [25,26] given by:

$$V\frac{d{x}^{C{O}_2}\left(t\right)}{dt}=q(x_{out}^{C{O}_2}-x^{C{O}_2}\left(t\right))+G_{C{O}_2}$$

(1)

where V denotes the volume of the room, $x^{C{O}_2}\left(t\right)$ the carbon dioxide concentration, $x_{out}^{C{O}_2}$ is the outdoor carbon dioxide concentration, q is a term associated with the ventilation and infiltration and window air exchange. $G_{C{O}_2}$ is a term associated with the generation of carbon dioxide in the room. The latter contains two terms: the carbon dioxide generation per occupant $c_{occ}^{C{O}_2}$ multiplied by a positive integer valued number $n_{occ}\left(t\right)$, denoting the number of persons, coming from the fact that the more the persons the higher the concentration in the room. If a closed window and mechanical ventilation free room are considered, the q term in Equation (1) is characterised only by the infiltration rate $q_{inf}$.

In order to adopt a more realistic model for the carbon dioxide concentration in the room we must treat the CO₂ time series as a stochastic process. By adopting a stochastic model, e.g., considering the CO₂ concentration data as a stochastic process, uncertain factors and other random unmodelled phenomena are taken into account [26]. The stochastic nature emerges not only from known factors, e.g., varying exhalation rates of the occupants in the room, the air exchange through doors, the non-homogeneous nature of the room’s air, as well as unknown factors that cannot be described in mathematical terms. As referred to in [26], the CO₂ differential equation contains difficult to exact modelling factors, even in the deterministic case, such as surface area on which indoor pollutants are deposited and the deposition velocity of the pollutant, among others. Such terms can be considered as unmodelled, since their exact mathematical expression in the equation may lack accuracy. The Brownian motion, which will be introduced to the mass balance differential equation, can model both known but uncertain, and unknown or unmodelled phenomena, that govern the indoor carbon dioxide concentration.

A formal definition of a stochastic process is given below:

Definition 1

([27]). Let the probability space (Ω,$\mathcal{F}$,$\mathcal{P}$), the ordered set T and the measurable space (E,$\mathcal{G}$). A stochastic process is a collection of random variables $X=\{X_t;t\in T\}$, where for each fixed $t\in T$, $X_t$ is a random variable in (Ω,$\mathcal{F}$,$\mathcal{P}$) to (E,$\mathcal{G}$). Ω denotes the sample space, and E denotes the state space of the stochastic process $X_t$.

Stochastic processes can be viewed as a function of both $t\in T$ and $\omega \in \Omega$, e.g., for fixed t we have a random variable, $X_t$, and for fixed $\omega$ in the sample space $\Omega$ we find the function $X_t\left(\omega \right):T\to E$, which is the realisation of the process. One of the most important stochastic processes is the Brownian motion [28].

Therefore, considering the CO₂ generation as a stochastic process the concentration is associated to a random variable, due to the unknown influencing random effects. Since a deterministic differential equation model has been used, Equation (1), the randomness must be encapsulated in the model. By adding a noise term, the ordinary differential equation can be extended to a stochastic differential equation. Stochastic differential equations are used to model dynamics which exhibit random behaviour. It is now important to state the connection of stochastic differential equations and stochastic processes, in order to justify why adopting such an approach makes perfect sense.

Our purpose is to solve the following stochastic differential equation of the general form

$$dX\left(t\right)=\underset{\mathit{drift}}{\underbrace{\mu (t,X(t);\theta )dt}}+\underset{\mathit{diffusion}}{\underbrace{g(t,X(t);\theta )d{B}_t}}$$

(2)

where $X\left(t\right)\in \mathbb{R}$, $\mu :[0,\infty )\times \mathbb{R}\times \Theta \to \mathbb{R}$ is the drift term and $g:[0,\infty )\times \mathbb{R}\times \Theta \to \mathbb{R}$ is the diffusion term. $B_t$ is a 1-dimensional standard Brownian motion. An initial condition is given $X\left(t\right)=x$, where $t=0$.

Let the filtration ${\mathcal{F}}_t$ generated by the Brownian motion $B_t$, i.e., ${\mathcal{F}}_t:=\sigma \left(\{B_s:s\in [0,t]\}\right)$, for every $t\ge 0$. Under mild conditions (Similarly to Lipschitz continuity and linear growth) on the functions $\mu (t,X(t);\theta )$ and $g(t,X(t);\theta )$ the solution of the stochastic differential equation exists and it is a unique process, $X\left(t\right),t\ge 0$ satisfying:

• It has continuous paths;
• It is ${\mathcal{F}}_t$-adapted;
• Almost surely, it holds that

  $$X\left(t\right)=x_0+{\int }_0^t\mu (s,X\left(s\right);\theta )ds+{\int }_0^{t}g(s,X\left(s\right);\theta )d{B}_s$$

  (3)

The solutions of stochastic differential equations are stochastic processes called diffusion processes. The mathematical treatment of stochastic differential equations is quite technical, however in the case of linear one-dimensional equations, closed form solutions could be derived. Moreover, as discussed in the next section, the solutions of stochastic differential equations are Markov processes [29], a property that turns out to be crucial in practical implementations.

2.3. The Hull–White Model as Mass Balance Equation

Since, we wish to treat the unknown process of carbon dioxide indoor generation and decay as a stochastic process, it is useful to model the mass balance equation with a stochastic differential equation. A model which originates from the financial mathematics literature, called the Hull–White Model (HWM) [29,30], has the suitable mathematical form in order for Equation (1) to be extended into a stochastic framework. This model has been introduced by John Hull and Alan White in 1994, as an extension to the Vasicek model for interest rates modelling. The Hull–White model is given by the following equation:

$$dR\left(t\right)=\left[\mathsf{\Phi }\left(t\right)-\alpha R\left(t\right)\right]dt+\sigma d{B}_t$$

(4)

where $R\left(t\right)$ denotes the interest rate, $\mathsf{\Phi }\left(t\right)$ denotes a time varying function, and $\alpha$ is the rate at which the rate is pushed towards the long term level.

Borrowing the mathematical structure of the HWM, a stochastic differential equation that is a generalisation of the Ornstein–Uhlenbeck model, can be written:

$$d{X}^{C{O}_2}\left(t\right)=[\underset{\mathsf{\Phi }\left(t\right)}{\underbrace{q_{inf}{X}_{out}^{C{O}_2}+c_{occ}^{C{O}_2}{n}_{occ}\left(t\right)}}-q_{inf}{X}^{C{O}_2}\left(t\right)]dt+\sigma d{B}_t$$

(5)

where $X^{C{O}_2}\left(t\right)$ is the carbon dioxide concentration in ppm (Appendix A), $q_{inf}$ denotes the infiltration rate in the room, i.e., the air change per hour, $\sigma$ denotes the constant volatility of changes in CO₂, due to the Gaussian white noise (Brownian motion’s increments are $\mathcal{N}$ distributed.) $d{B}_t$, $X_{out}^{C{O}_2}$ is the outdoor CO₂ concentration in ppm and $c_{occ}^{C{O}_2}$ is the CO₂ generation per person in ppm/h. Equation (5) is the stochastic mass balance differential equation.

We are dealing with a first order linear non autonomous stochastic differential equation, thus an analytical solution can be derived. The solution of Equation (5) can be found by applying Itô’s lemma on the function $e^{-q_{inf}}{X}^{C{O}_2}\left(t\right)$ and then integrating to get:

$$X^{C{O}_2}\left(t\right)=X^{C{O}_2}\left(s\right)e^{-q_{inf}(t-s)}+{\int }_s^t{e}^{-q_{inf}(t-\tau )}\mathsf{\Phi }\left(\tau \right)d\tau +\sigma {\int }_s^t{e}^{-q_{inf}(t-\tau )}d{B}_{\tau }$$

(6)

therefore, the expectation, mean solution, conditioned on the $\sigma$-algebra ${\mathcal{F}}_s$ is given by:

$$\mathbb{E}\left[X^{C{O}_2}\left(t\right)|{\mathcal{F}}_s\right]=X^{C{O}_2}\left(s\right)e^{-q_{inf}(t-s)}+{\int }_s^t{e}^{-q_{inf}(t-\tau )}\mathsf{\Phi }\left(\tau \right)d\tau$$

(7)

since, $\mathbb{E}\left[\sigma {\int }_s^t{e}^{-q_{inf}(t-\tau )}d{B}_{\tau }|{\mathcal{F}}_s\right]=0$. ${\mathcal{F}}_s$ is the filtration generated by the process up to time s. Moreover, the variance of the process is

$$\begin{array}{cc}\hfill \mathrm{Var}\left[X^{C{O}_2}\left(t\right)|{\mathcal{F}}_s\right]& =\mathbb{E}\left[{\left(X^{C{O}_2}\left(t\right)-\mathbb{E}\left[X^{C{O}_2}\left(t\right)|{\mathcal{F}}_s\right]\right)}^2|{\mathcal{F}}_s\right]\hfill \\ \hfill & =\mathbb{E}\left[{\left(\sigma {\int }_s^t{e}^{-q_{inf}(t-\tau )}d{B}_{\tau }\right)}^2\right]\hfill \\ \hfill & =\mathbb{E}\left[{\sigma }^2{\int }_s^t{e}^{-2q_{inf}(t-\tau )}{d}_{\tau }\right],\mathrm{due}\mathrm{to}\mathrm{Itô}’\mathrm{s}\mathrm{Isometry}\hfill \\ \hfill & =\frac{{\sigma }^2}{2q_{inf}}\left(1-e^{-2q_{inf}(t-s)}\right)\hfill \end{array}$$

(8)

The mean solution of the stochastic differential equation, describing $X^{C{O}_2}\left(t\right)$, is composed by an exponential term which forces the state to decay and an integral term which is due to the non autonomous nature of the differential equation: this term reflects the presence of occupants in the room at each time t, i.e., a piecewise constant function, such that $n_{occ}:[0,\infty )\to {\mathbb{Z}}^{+}$. The integral is fully tractable, thus a closed form solution is derived. The presence of occupants inside the room should, from a physics perspective, increase the levels of the carbon dioxide concentration, a fact that is reflected naturally by the above mathematical modelling.

The computation of the mean and variance describing the behaviour of the stochastic process $X^{C{O}_2}=\{X^{C{O}_2}\left(t\right);t\ge 0\}$ allows to write analytical expressions of the Likelihood function as the product of the transition densities from a state at time s to the state at time t, with $t\ge s$.

The original linear stochastic differential equation is equivalent [31] to the following discrete time system:

$$X_{t_{k+1}}^{C{O}_2}=f_k{X}_{t_k}^{C{O}_2}+u_k+q_k,\mathrm{where}{q}_k\sim \mathcal{N}(0,{\Sigma }_k)$$

(9)

where the coefficients $f_k$, $u_k$ and ${\Sigma }_k$ reads as:

$$\begin{array}{}\mathrm{(10)}& \hfill f_k& =\psi (t_{k+1},t_k)\hfill \mathrm{(11)}& \hfill u_k& ={\int }_{t_k}^{t_k+1}\psi (t_{k+1},\tau )\mathsf{\Phi }\left(\tau \right)d\tau \hfill \hfill \mathrm{(12)}& \hfill {\Sigma }_k& ={\sigma }^2{\int }_{t_k}^{t_k+1}{\psi }^2(t_{k+1},\tau )d\tau \hfill \end{array}$$

where ${\psi }_k(t_{k+1},t_k)=e^{-q_{inf}\Delta t}$, and $\Delta t=t_{k+1}-t_k$ These expressions can be analytically evaluated, to give the mean and variance of the original stochastic differential equation.

As discussed earlier, our purpose is to model the carbon dioxide concentration by a stochastic process. This is performed through the assumption that the generation and decay of carbon dioxide is governed by a linear non autonomous stochastic differential equation. The solution of the latter is a stochastic process and can be obtained either by using the discrete representation Equation (9), or by solving the original stochastic differential Equation (5) using a numerical integration scheme. To generate the training data, on which the parameters and occupancy estimation algorithm is based, the differential equation is solved using the method of Euler–Maruyama (EM), see Algorithm 1 [31].

Algorithm 1 Euler Maruyama

1: Given a stochastic differential equation of the form $dX\left(t\right)=\mu (t,X\left(t\right);\theta )dt+g(t,X\left(t\right);\theta )d{B}_t$, with known parameters $\theta$, an initial condition $X\left(t_0\right)$ and a time interval $[t_0,t]$. 2: Divide the time interval $[t_0,t]$ into N steps of $\Delta t$ length. 3: for Every time step k do 4:       Draw a random variable such as $\Delta B_k\sim (0,\Delta t)$ 5:       Compute the approximation: $\widehat{X}\left(t_{k+1}\right)=\widehat{X}\left(t_k\right)+\mu (t_k,X\left(t_k\right);\theta )\Delta t+g(t_k,X\left(t_k\right);\theta )\Delta B_k$ 6: end for 7: return The approximation $\widehat{X}\left(t_k\right)$ of the solution.

In this paper, CO₂ data are generated by simulating Equation (5), using the EM method. In order to be able to perform this, the parameters $\theta$, are considered to be known. In Figure 2, the evolution of the stochastic differential equation’s solution, with a single occupant, is depicted: Ten solution paths, each associated with a realisation of the Brownian motion, the mean solution and the $95\%$ quantiles.

In the following section, the method for the estimation of true parameters will be analysed, if an assumption of unknown parameters lies.

2.4. Estimation of the Stochastic Process Parameters

Since in real world applications the data of the carbon dioxide concentration may come in a time series, i.e., a set of available data, one might ask what the parameters of the process would be. Thus, an algorithm must be implemented in order to be able to estimate the true parameters of the problem. Henceforth, it is assumed that the solution time series of the carbon dioxide concentration is given, and the unknown parameters $\theta$ of the stochastic process are to be estimated.

We, therefore, have a problem of parameter estimation on stochastic differential equations: Given the stochastic differential Equation (5), the task is to estimate the parameters vector $\theta$ from a finite sample of available observations of the stochastic process $X^{C{O}_2}\left(t\right)=\{X^{C{O}_2}\left(t\right),\mathrm{where}t=t_0,\dots ,t_N\}$. The model is parametrised by the following quantities: the infiltration rate $q_{inf}$, the outdoor carbon dioxide concentration $X_{out}^{C{O}_2}$, the generation rate of carbon dioxide $c_{occ}^{C{O}_2}$ by a single occupant and the volatility of the noise $\sigma$.

Parameter estimation problem in probabilistic models has been investigated by many authors and several methods have been proposed. For our task, the well established method of maximum likelihood in the context of stochastic differential equations, is adopted. The likelihood function must be maximised, with respect to the unknown parameter vector $\theta$.

In light of the Markov property of the stochastic differential equation’s solution, an analytical expression of the likelihood of the observed data given the parameters is written by:

$$p(X^{C{O}_2}\left(t_0\right)\dots X^{C{O}_2}\left(t_N\right)|\theta )=\prod _{k=0}^{N-1}p\left(X^{C{O}_2}\left(t_{k+1}\right)\right|X^{C{O}_2}\left(t_k\right),\theta )$$

(13)

where, in general, $p\left(X\right(t\left)\right|X\left(s\right),\theta )$ is the transition density, describing the probability of the stochastic process to be from $X\left(s\right)$, to the state $X\left(t\right)$ in time $t-s$. Since the solution of the stochastic differential equation is a linear transformation of the Brownian motion, which is a Gaussian Process, it follows that the solution will also be Gaussian, thus, the transition densities are given by:

$$p\left(X^{C{O}_2}\left(t\right)\right|X^{C{O}_2}\left(s\right),\theta )=\mathcal{N}(\mathbb{E}\left(X^{C{O}_2}\left(t\right)\right|X^{C{O}_2}\left(s\right)),\mathrm{Var}\left(X^{C{O}_2}\left(t\right)\right|X^{C{O}_2}\left(s\right)))$$

(14)

where the Markov property is explicitly used, and the expressions for the distribution’s mean and variance are given by Equations (7) and (8), respectively.

The parameter vector $\widehat{\theta }=\left({\widehat{\theta }}_1,{\widehat{\theta }}_2,{\widehat{\theta }}_3,{\widehat{\theta }}_4\right)=\left({\widehat{q}}_{inf},{\widehat{X}}_{out}^{C{O}_2},{\widehat{c}}_{occ}^{C{O}_2},\widehat{\sigma }\right)$ must be computed, such that the likelihood to be maximised. In simple terms, the parameter space $\Theta$, is searched, for those parameters that generate the actual data, which in our case are the data from the numerical solution of the stochastic mass balance differential equation. For practical reasons, the parameters are determined in order for the negative log-likelihood function to be minimised, i.e.,

$$\widehat{\theta }=\underset{\begin{array}{c}{\theta }_i\in [l_{{\theta }_i},u_{{\theta }_i}]\\ i=1\dots 4\end{array}}{arg\; min}L\left(X^{C{O}_2}\left(t_0\right)\dots X^{C{O}_2}\left(t_N\right),\theta \right)$$

(15)

where $l_{{\theta }_i},u_{{\theta }_i}$ denote the lower and upper bounds of the feasible domain for each of the four parameters. $L(.)$ denotes the negative log-likelihood function given by:

$$L(X^{C{O}_2}\left(t_0\right)\dots X^{C{O}_2}\left(t_N\right),\theta )=-\sum _{k=0}^{N-1}logp\left(X^{C{O}_2}\left(t_{k+1}\right)\right|X^{C{O}_2}\left(t_k\right),\theta )$$

(16)

In this paper a differential evolution algorithm is proposed, in order to numerically minimise the negative log-likelihood, without the need of derivative information.

Differential Evolution (DE) [32], belongs to the class of evolutionary algorithms (EAs), that stochastically explore the space of the optimisation problem’s decision variables. DE uses a structure of a population of candidate solutions and guides the search procedure through a number of generations, until termination criteria have been met. The algorithm is parametrised by a minimum number of hyper-parameters: the mutation scaling F, the crossover probability $p_{cr}$ and the population size $n_p$.

The mutation of each individual comes in terms of perturbation of a target vector from the population with two weighted differential, described by:

$${\theta }_{i,\mathrm{Mutant}}^{g+1}={\theta }_{r_1}^g+F({\theta }_{r_2}^g-{\theta }_{r_3}^g)+F({\theta }_{r_4}^g-{\theta }_{r_5}^g)$$

(17)

for $i=1,\dots ,n_p$, where $r_1$ to $r_5$$\in \{1,2,\dots ,n_p\}$ and $i\ne r_1\ne r_2\ne r_3\ne r_4\ne r_5$. The mutation operator produces a population of perturbed model parameter vectors, called mutant vectors. Crossover is described by:

$$\begin{array}{c}\hfill {\theta }_{ji,\mathrm{Trial}}^{g+1}=\left\{\begin{array}{c}{\theta }_{ji,\mathrm{Mutant}}^{g+1}\mathrm{if}\mathbb{U}(0,1)\le p_{cr}\mathrm{or}j={\mathbb{U}}_i\hfill \\ {\theta }_{ji}^g\mathrm{otherwise}\hfill \end{array}\right.\end{array}$$

(18)

where ${\mathbb{U}}_i$ denotes the uniform distribution in integers $\{1,\dots ,4\}$ and $\mathbb{U}$ the uniform distribution in $(0,1)$. The candidate solutions will be represented by a population of individuals. For every generation g the population of the parameters is described:

$$\begin{array}{}\\ {\mathcal{P}}^g={\left\{{\mathcal{P}}_i\right\}}_{i=1\dots M}=\end{array}\begin{array}{ccc}& \begin{array}{c}{\mathcal{P}}_1^g& {\mathcal{P}}_2^g& \dots & {\mathcal{P}}_{n_p}^g\end{array}& \\ (& \begin{array}{cccc}{\theta }_1^1& {\theta }_1^2& \dots & {\theta }_1^{n_p}\\ ⋮& ⋮& ⋮& ⋮\\ {\theta }_4^1& {\theta }_4^2& \dots & {\theta }_4^{n_p}\end{array}& )\end{array}$$

(19)

The proposed method is summarised in Algorithm 2.

Algorithm 2 Parameter estimation using likelihood function and differential evolution

1: Generate Data by numerically integrating the stochastic mass balance differential Equation (6) using Euler–Maruyama method. 2: Compute The mean and variance solutions of the SDE analytically 3: Compute The transition densities and the negative log-likelihood function. 4: Given Data from EM solution of the Stochastic Differential Equation: $\mathcal{D}={\left\{(t_k,X^{C{O}_2}\left(t_k\right))\right\}}_{k=1,\dots ,N}$. 5: Initialise Differential evolution variant DE/Rand/2/Bin, where Rand describes the mutation operator, 2 indicates the number of differentials, and Bin the crossover scheme. 6: Set Mutation parameter F, crossover probability $p_{cr}$, population size $n_p$ and a maximum number of generations: MaxGen. 7: Generate Initial population of candidate solutions 8: ${\mathcal{P}}^0\leftarrow Uniform(l_{{\theta }_i},u_{{\theta }_i},n_p)$ 9: forEvery Generation gdo 10:     for all individuals in ${\mathcal{P}}^g$ do 11:           Mutate 12:           ${\mathcal{P}}_{mutation}^{g+1}\leftarrow$ Mutation(${\mathcal{P}}^g,F$) 13:           if Bound constraints $[l_{{\theta }_i},u_{{\theta }_i}]$ are violated then 14:                 ${\mathcal{P}}_{mutation}^{g+1}\leftarrow$ Resampling(${\mathcal{P}}_{mutation}^{g+1},l_{{\theta }_i},u_{{\theta }_i})$ 15:           end if 16:           Crossover 17:           ${\mathcal{P}}_{trial}^{g+1}\leftarrow$ Crossover(${\mathcal{P}}^g,{\mathcal{P}}_{mutation}^{g+1},p_{cr}$) 18:           Selection 19:           if $L\left({\mathcal{P}}^g\right)<L\left({\mathcal{P}}_{\mathrm{trial}}^{g+1}\right)$ then ${\mathcal{P}}^{g+1}\leftarrow {\mathcal{P}}^g$ 20:           else${\mathcal{P}}^{g+1}\leftarrow {\mathcal{P}}_{\mathrm{trial}}^{g+1}$ 21:           end if 22:       end for 23:       $g\leftarrow g+1$ 24: end for 25: Return The Estimated Parameter vector $\widehat{\theta }=\left({\widehat{q}}_{inf},{\widehat{X}}_{out}^{C{O}_2},{\widehat{C}}_{occ}^{C{O}_2},\widehat{\sigma }\right)$

In the latter algorithm, it is assumed that the data are available to us. To achieve this, we synthetically generate them, by solving the stochastic differential equation using the Euler–Maruyama method. The function Resampling implements the resampling of the DE solutions into the feasible domain, since the mutation operation may produce solutions that are outside of the parameter’s bound constraints. The maximum number of generations, MaxGen is set to 100, the population size $n_p=80$, the crossover probability $p_{cr}=0.2$ and a mutation scaling F sampled uniformly in $[0.2,2]$, in every generation.

The justification for the choice of the hyper parameters is based on [32]: Empirical studies on the simulations performed, show that the convergence of the algorithm occurs at around 70th generation, hence a maximum number of 100 generations would suffice. The probability of crossover, $p_{cr}$ is chosen to be small in order to avoid premature convergence and moreover, to increase the search robustness. The mutation scaling is chosen to be given by a uniform distribution in the range $[0.2,2]$, in order to achieve a trade-off between exploring tight valleys and search diversity. As far as the population size is concerned, the choice of 80 achieves a reasonably fast convergence, and it is chosen on a trial and error basis.

2.5. Fuzzy Occupancy Estimation

It should be clearly stated, that up to this point the unknown (to be estimated) parameters were considered to be: the infiltration rate $q_{inf}$, the outdoor carbon dioxide concentration $X_{out}^{C{O}_2}$, the carbon dioxide generation per occupant $c_{occ}^{C{O}_2}$ and the volatility of the noise term $\sigma$. Since the number of occupants in Equation (5), can be considered as a parameter too, the problem is addressed by likewise asking the same question as in the parameter estimation context: which occupancy vector generates the observed data.

Hence, if the model is re-parametrised with a single unknown parameter $n_{occ}\left(t\right)$ and the estimated parameters $\widehat{\theta }$ are considered to generate the true data in high accuracy, then a likelihood approach can be used to determine the occupancy vector. We followed and implemented the occupancy estimation algorithm proposed in [14] and estimated the occupancy vector at every time instance. The occupancy vector $n_{occ}$ is very high dimensional, thus an optimisation based approach even in terms of evolutionary algorithms has, to our knowledge, not been studied. Wolf et al. propose a greedy algorithm, based on maximum likelihood’s principle, which initially considers the occupancy vector as a zero vector, representing an empty room, and iteratively adds integers values to the vector points, i.e., adds an occupant in each time step, where the likelihood function is maximised. For completeness, an implementation of the algorithm is given in Algorithm 3.

The reason that the estimated parameters $\widehat{\theta }$ serve as inputs in the algorithm, is that in a real world setting the data observation will be given as a time-series, where the model’s parameters will be unknown. In the case when the estimated parameters are close to the ground truth, then the estimated model will generate very similar data and the only parameter that will be to be determined is the occupancy vector. The performance of the algorithm has been tested under two different occupancy profiles, with very accurate results. It is crucial to note that the method of parameter estimation, in our case the differential evolution algorithm, must result in close values of the true parameters $\theta$ otherwise, as explored next, the occupancy estimation algorithm will not output the correct estimation vector.

Algorithm 3 [14] Occupancy estimation algorithm

1: Inputs: Observed Data$X^{C{O}_2}\left(t\right)$, Estimated Parameters$\widehat{\theta }$, Likelihood Function$L\left(n_{occ}\left(t\right)\right)$ 2: Initialise ${\widehat{n}}_{occ}\leftarrow {\mathbf{0}}_{1\times \mathrm{Number}\mathrm{of}\mathrm{Observations}}$ $\mathcal{A}\leftarrow {\mathcal{I}}_{\mathrm{Number}\mathrm{of}\mathrm{Observations}\times \mathrm{Number}\mathrm{of}\mathrm{Observations}}$ Compute the value of the negative log-likelihood function: $L_{{\widehat{n}}_{occ}}\leftarrow L(X^{C{O}_2}\left(t_0\right)\dots X^{C{O}_2}\left(t_N\right),\widehat{\theta },{\widehat{n}}_{occ})$ 3: for a number of iterations do 4:       for every time step, $i=1\dots N$ do 5:             $n_{occ}^i\leftarrow {\widehat{n}}_{occ}+\mathcal{A}(i,:)$ 6:             $L_{n_{occ}^i}\leftarrow L(X^{C{O}_2}\left(t_0\right)\dots X^{C{O}_2}\left(t\right),\widehat{\theta },n_{occ}^i)$ 7:       end for 8:       $n_{occ}^{min}\leftarrow {arg\; min}_{\begin{array}{c}{n}_{occ}^i\end{array}}{L}_{n_{occ}^i}$ 9:       Greedy Selection 10:      if $L_{n_{occ}^{min}}<L_{{\widehat{n}}_{occ}}$ then ${\widehat{n}}_{occ}\leftarrow n_{occ}^{min}$ 11:      else ${\widehat{n}}_{occ}\leftarrow {\widehat{n}}_{occ}$ 12:      end if 13: end for 14: Output: Estimation of the occupancy vector ${\widehat{n}}_{occ}$

The methodology of parameter $\theta$ estimation, via the maximum likelihood approach, results in point estimates. Since the true parameters that generated the data of carbon dioxide concentration are unknown, a question of how certain about these values we could be, emerges. As earlier discussed there are two kinds of uncertainty: the uncertainty that stems from random effects and unmodelled dynamics and the structural uncertainty, which corresponds to possible fluctuations of the model parameters. Stochastic differential equations can be considered to model the first case. Pertaining to the second type of uncertainty, a solution can be the use of fuzzy modelling [33]. The subject of stochastic differential equations with fuzzy parameters, leading to Fuzzy stochastic differential equations, have been studied, from a solution point of view, in [34,35].

In this paper, the effect of this structural parametric uncertainty onto the occupancy estimation algorithm, is investigated. The estimated parameters are mapped onto the fuzzy space using triangular fuzzy numbers, thus an uncertainty is induced on the model’s parametric space. Fuzzy numbers of the following form have been used:

$$\mu (x;x_L,x_C,x_R)=max\left(0,min\left(\frac{x-x_L}{x_C-x_L},\frac{x_R-x}{x_R-x_C}\right)\right)$$

(20)

where x denotes a point in the domain $\mathcal{X}$ of the parameter to be fuzzified, $[x_L,x_R]$ denotes the support of the fuzzy number and $x_C$ denotes the core of the fuzzy number. Note that $\mu$ here symbolises the membership function of the fuzzy set, not to be confused with the drift term in Equation (2).

The assumption of structural uncertainty has been made for two of the estimated parameters: the infiltration rate $q_{inf}$ and the carbon dioxide generation per occupants $c_{occ}^{C{O}_2}$. This reflects the fact that in a real scenario the infiltration may have small deviations from the true value [36]. Similarly, the CO₂ generation per occupant in the room, $c_{occ}^{C{O}_2}$, in any case cannot be a crisp value, since it strongly depends on factors, such as the type of work that an occupant performs during its presence in the room, the age of the occupant, the gender of the occupant, and various physical characteristics [37,38]: a young child’s carbon dioxide generation differs from an adult heavy worker’s.

Smaller uncertainty bounds are induced on the infiltration parameter, and larger on the CO₂ generation per occupant. The fuzzification of the latter parameters is implemented by a symmetric triangular fuzzy number, given by Equation (18), where the supports have been calculated by ${\widehat{\theta }}_1\pm 10\%$ and ${\widehat{\theta }}_3\pm 30\%$ for the infiltration parameter and the carbon dioxide generation, respectively. The estimated parameters ${\widehat{\theta }}_1$ and ${\widehat{\theta }}_3$ are the cores of the fuzzy sets that describe the uncertainty of the infiltration and CO₂ generation. Figure 3 depicts the fuzzy sets, their support and the crisp uncertainty set, corresponding to the $0.5$-level.

To study the result of the fuzzy parameters in the occupancy estimation algorithm, the a-cuts $A_a$ have been computed using the definition: $A_a=\{x\in \mathcal{X},m\left(x\right)\ge a\}$. A $0.5$-level is considered, meaning that a trust region is given, in the values of the fuzzified parameters with membership $\mu \ge 0.5$.

The limit points of the $a-cuts$ are given by:

$$\begin{array}{}\hfill \mathrm{(21)}& \hfill {\underline{\widehat{\theta }}}_1& =a\left({\widehat{\theta }}_1-{\tilde{l}}_{{\widehat{\theta }}_1}\right)+{\tilde{l}}_{{\widehat{\theta }}_1}\hfill \hfill \mathrm{(22)}& \hfill {\overline{\widehat{\theta }}}_1& ={\tilde{r}}_{{\widehat{\theta }}_1}-a({\tilde{r}}_{{\widehat{\theta }}_1}-{\widehat{\theta }}_1)\hfill \hfill \mathrm{(23)}& \hfill {\underline{\widehat{\theta }}}_3& =a\left({\widehat{\theta }}_3-{\tilde{l}}_{{\widehat{\theta }}_3}\right)+{\tilde{l}}_{{\widehat{\theta }}_3}\hfill \hfill \mathrm{(24)}& \hfill {\overline{\widehat{\theta }}}_3& ={\tilde{r}}_{{\widehat{\theta }}_3}-a({\tilde{r}}_{{\widehat{\theta }}_3}-{\widehat{\theta }}_3)\hfill \end{array}$$

where ${\tilde{l}}_{{\widehat{\theta }}_3},{\tilde{r}}_{{\widehat{\theta }}_3}$ are the left and right limits of the support for the parameter ${\widehat{c}}_{occ}^{C{O}_2}$ and ${\tilde{l}}_{{\widehat{\theta }}_1},{\tilde{r}}_{{\widehat{\theta }}_1}$ are the left and right limits of the support for the parameter ${\widehat{q}}_{inf}$. It should be clear that the sets $A_a$ are uncountable and we are focusing on investigating the boundary case, i.e., the model with parameters chosen firstly as ${\underline{\widehat{\theta }}}_1,{\underline{\widehat{\theta }}}_3$ and secondly ${\overline{\widehat{\theta }}}_1,{\overline{\widehat{\theta }}}_3$.

The implementation of the occupancy estimation algorithm, results in different occupancy vectors $n_{occ}\left(t\right)$, reflecting the fact that when fuzzy parameters are considered, the estimation gives an upper and lower estimation of the number of occupants in the room at each time.

2.6. Fuzzy $\mathcal{Q}$-Learning and Multi Agent Control System

Reinforcement learning refers to a group of algorithms that have been inspired from the learning mechanism of humans or other leaving beings. Reinforcement learning algorithms aim to extract a policy via exploring the possible state-action pairs. The term policy refers to the mapping between states and actions. All reinforcement learning algorithms need to utilise a signal called reward, which defines a measure of how well an agent performs at a given state.

A tabular form of reinforcement learning, is $\mathcal{Q}$-learning [39]. In the $\mathcal{Q}$-learning framework, a $\mathcal{Q}$-function, that tries to estimate the future discounted rewards for actions in given states, is gradually built by the agent. The $\mathcal{Q}$-function output for a given state x and an action a is defined as $\mathcal{Q}(x,a)$ and the values of the $\mathcal{Q}$-function are stored in a $\mathcal{Q}$-table. $\mathcal{Q}$-learning’s main drawback is that it cannot be applied for continuous state-action space. Fuzzy $\mathcal{Q}$-learning [40] can be used to overcome this difficulty and it is summarised in the following steps:

• Statexobservation;
• Output selection, for each fired fuzzy rule, based on the exploration/exploitation algorithm;
• Global output’s $\alpha \left(x\right)$ and corresponding value’s of $\mathcal{Q}(x,\alpha )$ computation by:

  $$\alpha \left(x\right)=\frac{{\sum }_{i=1}^M{\alpha }_i\left(x\right){\alpha }_i}{{\sum }_{i=1}^M{\alpha }_i\left(x\right)}$$

  (25)

  $$\mathcal{Q}\left(x,\alpha \right)=\frac{{\sum }_{i=1}^M{\alpha }_i\left(x\right){\alpha }_{i}q[i,i^{†}]}{{\sum }_{i=1}^M{\alpha }_i\left(x\right)}$$

  (26)
  where M is the number of the fired fuzzy rules,${\alpha }_i\left(x\right)$is the fired degree for each fuzzy rule,${\alpha }_i$is the action that is selected for each fuzzy rule and$q[i,i^{†}]$is the corresponding$\mathcal{Q}$value for the pair of the fuzzy rule i and the action$i^{†}$. The action corresponds to the action that the exploration/exploitation algorithm selects.
• The global action $\alpha \left(x\right)$ is applied and the next state $x^{\prime }$ is observed;
• Reward R computation;
• $\mathcal{Q}$-values are updated according to the formula:

  $$\Delta q[i,i^{*}]=\eta \Delta \mathcal{Q}\frac{{\alpha }_i\left(x\right)}{{\sum }_{i=1}^M{\alpha }_i\left(x\right)}$$

  (27)
  where $\Delta \mathcal{Q}=R+\gamma V\left(x^{\prime }\right)-\mathcal{Q}(x,\alpha )$ and

  $$V\left(x^{\prime }\right)=\frac{{\sum }_{i=1}^M{\alpha }_i\left(x^{\prime }\right){\alpha }_{i}q[i,i^{*}]}{{\sum }_{i=1}^M{\alpha }_i\left(x^{\prime }\right)}$$

  (28)
  and $q[i,i^{*}]$ is the corresponding $\mathcal{Q}$ value for the pair of the fuzzy rule i and the action $i^{*}$. γ denotes the discount factor and η the learning rate. The action $i^{*}$ corresponds to the action that has the maximum corresponding $\mathcal{Q}$ value for the fuzzy rule i.

For solving complex or physically distributed problems, a multi agent control system can be used. $\mathcal{Q}$-learning and consequently fuzzy $\mathcal{Q}$-learning can be applied to MACS with the most common approach being independent learners [41], where each agent ignores the existence of other agents and interacts independently with the environment.

2.7. System Modelling and Description

2.7.1. Heating and Cooling System

The space is a building office with a surface of 50 m² and volume 150 m³ (5 m × 10 m × 3 m) and a 3 m² (3 m × 1 m) window which is located in the north side of the office. The heating/cooling system of the building introduce either hot or cold air into the space [42,43]. The heat flow into the space can be computed by the next equation:

$$\frac{dQ\left(t\right)}{dt}=\pm \left(T_{\frac{heat}{cool}}-T_{room}\right)\dot{M}c+H{G}_o{n}_{occ}\left(t\right)+H{G}_E{n}_{occ}\left(t\right)+H{G}_L+H{G}_S$$

(29)

where the left hand side denotes the heat flow rate into the space. $T_{\frac{heat}{cool}}$ denotes the temperature of the hot and the cold air 40 °C or 16 °C, respectively. $T_{room}$ is the temperature of the room’s air, $\dot{M}$ is the mass air flow of the air from the heating/cooling system (2600 kg/h) and c is the air heat capacity at constant pressure 1005.4 J/kg °C. $H{G}_o$ denotes the heating gain from the occupants (115 W/occupant) [44], $n_{occ}\left(t\right)$ is the number of occupants and $H{G}_E$ is the heating gain from the electrical loads (5.38 W/m${}^2\times 50$m ${}^2=269$ W) [36], $H{G}_L$ is the heating gain from the artificial lighting $(0.5\times 800=400\mathrm{W})$ and $H{G}_S$ is the heating gain from the solar irradiance through the window [45] which is computed as:

$$H{G}_S=A_w\times r_w\times D\times I_w-g_w\times A_w\times (T_{room}-T_{out})$$

(30)

where $A_w$ denotes the area of window, $r_w$ denotes the transparency of the window, $D=1$ denotes the window shadowing factor, $I_w$ is the diffused solar irradiance that reaches the window, $g_w$ is the glass conductance (2 W/m² °C) and $T_{out}$ is the temperature of the outdoor air. The thermal model of the building does not take into consideration the effect of the solar irradiance through the building opaque surfaces. The derivative of the space temperature over time, is expressed as:

$$\frac{d{T}_{room}}{dt}=\pm \frac{1}{M_a·c}\left(\frac{dQ\left(t\right)}{dt}-\frac{d{Q}_l\left(t\right)}{dt}\right)$$

(31)

where $M_a$ is the air mass inside the building (183.75 kg) and:

$$\frac{d{Q}_l\left(t\right)}{dt}=\frac{T_{room}-T_{out}}{R_{eq}}$$

(32)

where $R_{eq}$ is the equivalent thermal resistance of the building envelope ($1.0306\times {10}^{-6}$°C/W). The equivalent thermal resistance is computed by the following set of equations [42].

$$\begin{array}{c}\hfill R_{eq}=\frac{R_{w\alpha }{R}_{wi}}{R_{w\alpha }+R_{wi}},R_{w\alpha }=\frac{T{h}_{w\alpha }}{k_{w\alpha }{A}_{w\alpha }},R_{wi}=\frac{T{h}_{wi}}{k_{wi}{A}_w}\end{array}$$

(33)

where $A_w$ is the area of the windows as given in Equation (25), $T{h}_{wi}$ is the thickness of the windows (0.01 m), $k_{wi}$ is the thermal conductivity of the windows (0.78 W/m °C), $A_{w\alpha }$ is the area of the walls (187 m₂), $T{h}_{w\alpha }$ is the thickness of the walls (0.2 m), $k_{w\alpha }$ is the thermal conductivity of the walls (0.038 W/m °C).

The heat loss/gain through the building envelope (walls area, windows area, and type of insulation) is embedded inside the equivalent thermal resistance [42]. The positive or negative sign on Equations (24) and (26), refers to the process of heating or cooling, respectively.

2.7.2. Lighting System

In order to compute the diffused illuminance provided by the sky inside the building the following equation is used. The horizontal illuminance in a point p [46] which is located in the interior of the building is computed as:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill E_{p,d}=\frac{r_w{L}_s}{2}\frac{z}{\sqrt{h_p^2+z^2}}\left(ta{n}^{-1}\frac{x_w+w_w-x_p}{\sqrt{h_p^2+z^2}}+ta{n}^{-1}\frac{x_w+x_p}{\sqrt{h_p^2+z^2}}\right)\\ \hfill -\frac{z}{\sqrt{{(h_p+h_w)}^2+z^2}}(ta{n}^{-1}\frac{x_w+w_w-x_p}{\sqrt{{(h_p+h_w)}^2+z^2}}\\ \hfill +ta{n}^{-1}\frac{x_p+x_w}{\sqrt{{(h_p+h_w)}^2+z^2}})\end{array}\end{array}$$

(34)

where $E_{p,d}$ is the horizontal diffused illuminance of the point p, $L_s$ is the illuminance from the sky, z is the horizontal distance between the point p and the window (10 m close to the south wall of the building which is amongst the most shading points inside the building). $h_p$ denotes the height between the lower edge of the window and the point p (0.5 m), $h_w$ is the height of the window (1 m), $w_w$ is the width of the window (3 m), $x_w$ is the distance between the left wall and the left edge of the window (3.5 m) and $x_p$ is the distance between the point p and the left wall (2.5 m).

It is assumed that the horizontal diffused illuminance, provided by the sky, is distributed uniformly and equals to the value of the point p, and consequently the control decisions are taken according to this value. Fluorescence lamps are used for the artificial lighting and the ratio between the horizontal luminance and the consumed power is 100 lm/W [47]. The electro-chromic window [48], can change its transparency by providing a small amount of voltage. The nominal power of the artificial lighting is 400 W which can provide a maximum illuminance of 800 $lux$ in a surface of 50 m${}^2$. The relation between the $lux$ and the consumption is assumed to be linear. Moreover, it is assumed that the horizontal illuminance provided by the lamps is distributed uniformly all over the surface of the office.

2.7.3. Ventilation System

The CO₂ mass balance with no noise consideration and with a given ventilation rate is given by:

$$\frac{d{x}^{C{O}_2}\left(t\right)}{dt}=\frac{(q_{vent}+q_{inf})(x_{out}^{C{O}_2}-x^{C{O}_2}\left(t\right))+c_{occ}^{C{O}_2}{n}_{occ}\left(t\right)}{V}$$

(35)

where V is the volume of the space (150 × 10${}^3$ L), and $q_{vent}$ is the ventilation rate, which is now introduced by the MACS. According to [36,44] for a building office, the air change rate must be between 2 and 8. For a space of $150\times {10}^3$ L the ventilation rate of the ventilation system can be $600\times {10}^3$ L/h. The installed power for such a system is approximately 165 W by assuming a relationship of 0.5 W per flow rate (L/s). As we will discuss in the simulation section, the generation rate per person has been estimated to $26.9427$ L/h, the outdoor carbon dioxide concentration has been estimated to 359.4813 ppm and the infiltration rate has been estimated to 0.1167 air change rate which equals to $17.535$ L/h. At this point, it should be noted that mechanical ventilation affects the indoor temperature, however nowadays mechanical ventilation systems incorporate heat recovery that can reach efficiency of almost $90\%$. Based on that fact, the effect of the ventilation flow rate to the indoor temperature is neglected from the building’s thermal model.

2.8. Multi Agent Control System

The MACS consists of a group of three agents: $A=\{A{G}_T,A{G}_L,A{G}_{CD}\}$, where $A{G}_T$ is the temperature agent, $A{G}_L$ is the illuminance agent, and $A{G}_{GD}$ the CO₂ agent. The states of the agents are defined by a total number of six fuzzy state variables $X_i$, i.e., two variables for each agent. For each positive values input, five triangular membership functions (MFs) have been used. In addition to that, for each both positive and negative values input, seven triangular MFs have been used (Figure 4). Membership functions are denoted as linguistic variables: PVB, PB, PM, PS, Z, NS, NM, and NB, e.g., positive very big, positive big, positive medium, positive small, zero, negative small, and negative medium, respectively.

$A{G}_T$ has two variable inputs: the outdoor temperature which is in the range $[0,1]$, and the error $e_T$ between the set point and the indoor temperature, normalised in the range $[-1,1]$. The total states are 35, which are represented by an equal number of fuzzy rules. The output vector of the temperature agent has five fuzzy singleton sets $A_T=\{\frac{1}{-0.05}+\frac{1}{-0.1}+\frac{1}{0}+\frac{1}{0.01}+\frac{1}{0.05}\}$.

The global action is the control signal, that defines the percentage of the air flow rate by the heating/cooling system, according to its maximum capacity. Positive signal actuate the heating system while negative signal actuate the cooling system. The reward $R_{A_T}$ of this agent is defined as follows:

$$R_{A_T}(x,\alpha ,x^{\prime })=-\left|e_T\right|-0.1·P_{HC}$$

(36)

where $P_{HC}$ is the power consumption of the heating/cooling system.

$A{G}_L$ has two variable inputs, the indoor horizontal illuminance, normalised in the range $[0,1]$ and the error $e_L$ between the set point and the indoor illuminance normalised in the range $[-1,1]$. The total states are 35, which are represented by an equal number of rules. The output vector of the $A{G}_L$ has five fuzzy singleton sets, $A_L=\{\frac{1}{-0.3}+\frac{1}{-0.02}+\frac{1}{0}+\frac{1}{0.02}+\frac{1}{0.3}\}$. The global action defines the percentage of the power change to be consumed by the artificial lighting system according to its nominal operating power for positive signal, while negative signal change the transparency of the electro-chromic window. The reward $R_{A_L}$ of this agent is defined as follows:

$$R_{A_L}(x,\alpha ,x^{\prime })=-\left|e_L\right|-0.1·P_{AL}$$

(37)

where $P_{AL}$ is the power consumption of the artificial lighting system.

$A{G}_{CD}$ has two variable inputs, the number of the occupants normalised in the range $[0,1]$ and the error $e_{CO}$ between the set point and the indoor CO₂ concentration, normalised in the range $[0,1]$. The total states (Five membership functions for each input. See plot (a) on Figure 4) are 25, which are represented by an equal number of rules. The output vector of the agent has five fuzzy singleton sets, $A_{CD}=\{\frac{1}{-0.4}+\frac{1}{-0.02}+\frac{1}{0}+\frac{1}{0.02}+\frac{1}{0.04}\}$. The global action is the control signal which represents the percentage of the power, change to be consumed by the ventilation system according to its maximum power. The reward $R_{A_{CD}}$ of this agent is defined as follows:

$$R_{A_{CD}}(x,\alpha ,x^{\prime })=-\left|e_{co}\right|-0.1·P_V$$

(38)

where $P_V$ is the power consumption of the ventilation system.

All three agents use the same exploration/exploitation algorithm. When a new state is visited, the agent performs exploration for 500 rounds. Following the 500 rounds, the agent checks and performs the actions that has not been performed at all. By this way it is assured that all possible actions are explored. Consequently, the agent performs exploitation for $99\%$ and exploration for $1\%$ for a given state.

Additionally, $A{G}_T$ and $A{G}_L$ operate only when occupancy is detected inside the building space. For periods, where the occupancy vector is zero, i.e., ${\widehat{n}}_{occ}\left(t\right)=0$, the agent’s control signal is set to zero and the learning mechanism stops functioning. Thus, when the building is empty a waste of energy is avoided.

3. Results

3.1. Parameters and Occupancy Estimation

The proposed methodology has been implemented to estimate the parameters $\theta$, as well as the occupancy vector for a set of given data, that stem from the EM solution of the stochastic differential Equation (5). All mentioned algorithms have been programmed in MATLAB. Figure 5 illustrates ten solution paths of the stochastic differential equation, where the parameters are chosen as: $q_{inf}=0.12{\mathrm{h}}^{-1}$ [36], $X_{out}^{C{O}_2}=400$ (ppm) [10], $c_{occ}^{C{O}_2}=180$ (ppm/h), [14], and $\sigma =2$ [14]. A profile $n_{occ}\left(t\right)$ materialises a simple scenario of presence in a room.

In order to replicate a realistic problem, a single realisation of the solution is considered, and the time series is used to estimate the parameters using the method of evolutionary maximisation of the likelihood function. In the last figure, it is clear that the solution of the stochastic mass balance equation reflects the fact that when occupants are present in the room, e.g., approximately at 13:40 to almost 21:00, the CO₂ increases as expected. When $n_{occ}\left(t\right)$ vanishes, then the solution has a mean reverting behaviour (Ornstein–Uhlenbeck Process.), analogous to the infiltration rate.

The estimation algorithm has been implemented using a maximum number of generations $MaxGen=100$, a population of $n_p=80$, a crossover probability $p_{cr}=0.2$ and a mutation scaling F sampled uniformly from $[0.2,2]$. The estimated parameters of the proposed method are given in

Table 1. True and Estimated model parameters.

Parameters	True Values $\theta$	Estimated Values $\widehat{\theta }$
Infiltration rate	0.12 h${}^{-1}$	0.1166 h${}^{-1}$
Outdoor CO2 concentration	400 ppm	391.4339 ppm
CO2 generation per occupant	180 ppm/h	179.1345 ppm/h
$\sigma$	2	3.8800

.

Figure 6, illustrates the solution of the stochastic model, i.e., the training data, and the model’s solution corresponding to the estimated parameters. The solution of the differential equation using the estimated parameters $\widehat{\theta }$ is in very close vicinity to the true solution.

The occupancy estimation algorithm, has been implemented in the latter case and the estimation ${\widehat{n}}_{occ}\left(t\right)$ is identical to the ground truth, $n_{occ}\left(t\right)$. The result is illustrated in Figure 7.

The performance of both the parameter estimation method, as well as the occupancy estimation algorithm is excellent. For the implementation of the occupancy algorithm, 1440 points of CO₂ concentration are used, after having generated them from the numerical solution of the mass balance equation; using the method of Euler–Maruyama. Furthermore, the estimated parameter $\widehat{\theta }$ has been used, as described in the section of occupancy estimation algorithm. The number of iterations for the algorithm has been set to 1000.

To study the result of parametric uncertainty on the occupancy estimation algorithm, the estimated infiltration rate and the estimated CO₂ generation per occupant, are modelled using fuzzy numbers. The parameters of the fuzzy case are computed using the a-cut method, as previously described, and are illustrated in Figure 3. The stochastic differential equation solutions using the estimated parameters, as well as the fuzzy estimated parameters (

Table 2. Estimated and Fuzzy parameters.

Estimated Parameters $\widehat{\theta }$	Fuzzy Parameters $\underline{\widehat{\theta }}$	Fuzzy Parameters $\overline{\widehat{\theta }}$
${\widehat{\theta }}_1=0.1166$ h${}^{-1}$	$0.1107$ h${}^{-1}$	$0.1224$ h${}^{-1}$
${\widehat{\theta }}_2=391.4339$ ppm	$391.4339$ ppm	$391.4339$ ppm
${\widehat{\theta }}_3=179.1345$ ppm/h	$152.2644$ ppm/h	$206.0047$ ppm/h
${\widehat{\theta }}_4=3.8800$	$3.8800$	$3.8800$

) are illustrated in Figure 8. It is clear that a parametric uncertainty bound is introduced in the model, since fuzzy parameters are considered.

The occupancy estimation algorithm is implemented in the case of both $\underline{\widehat{\theta }}$ and $\overline{\widehat{\theta }}$ parameters. The results are depicted in Figure 9. The estimation result produced by the algorithm differs from the real number of occupants, that generate the CO₂ data. By considering the infiltration and the carbon dioxide generation by the occupants as fuzzy, a further uncertainty in the model other than the noise, is induced. The occupancy vector in the case of $\underline{\widehat{\theta }}$ is, in some time points, overestimated. This reflects the fact that since CO₂ generation is much lower than the estimated value, then the algorithm, based on the maximum likelihood principle, tries to determine what $n_{occ}\left(t\right)$ should be in order for the training data to be generated. Therefore, it computes an estimation vector ${\widehat{n}}_{occ}\left(t\right)$, so that the solution, if considered parameters $\underline{\widehat{\theta }}$, would match the true data. In the case of $\overline{\widehat{\theta }}$ the estimated occupancy vector is, in some time points, underestimated reflecting the fact that in some time points the number of occupants should be smaller than 4, in order for the true data to be produced.

For the MACS development, the proposed methodology is used, for the estimation of an occupancy profile that resembles a working office: an increasing number of occupants are present in the room’s interior until late afternoon. Following to that, a single occupant is present for two hours, until midnight. The occupancy profile and the associated CO₂ concentration, are illustrated in Figure 10.

The true parameters used in the simulation, as well as their estimation, using differential evolution, are given in

Table 3. True and Estimated Model parameters.

Parameters	True Values $\theta$	Estimated Values $\widehat{\theta }$
Infiltration rate	0.12 h${}^{-1}$	0.1167 h${}^{-1}$
Outdoor CO2 concentration	420 ppm	359.4813 ppm
CO2 generation per occupant	180 ppm/h	179.6181 ppm/h
$\sigma$	2	2.8917

.

The algorithm of occupancy estimation, performs very accurately and the true occupant vector $n_{occ}\left(t\right)$ is identical to the estimated, ${\widehat{n}}_{occ}\left(t\right)$. The results are given in Figure 11.

Similarly to the previous case, the parametric uncertainty on the parameters is addressed, by assuming that the infiltration and the CO₂ generation per occupant are modelled by triangular fuzzy numbers, with core given by ${\widehat{\theta }}_1$ for $q_{inf}$ and ${\widehat{\theta }}_3$ for $c_{occ}^{C{O}_2}$. Applying the a-cut methodology the left and right parameters were computed to be (

Table 4. Estimated and Fuzzy parameters.

Estimated Parameters $\widehat{\theta }$	Fuzzy Parameters $\underline{\widehat{\theta }}$	Fuzzy Parameters $\overline{\widehat{\theta }}$
${\widehat{\theta }}_1=0.1167$ h${}^{-1}$	$0.1108$ h${}^{-1}$	$0.1225$ h${}^{-1}$
${\widehat{\theta }}_2=359.4813$ ppm	$359.4813$ ppm	$359.4813$ ppm
${\widehat{\theta }}_3=179.6181$ ppm/h	$152.6754$ ppm/h	$206.5608$ ppm/h
${\widehat{\theta }}_4=2.8917$	$2.8917$	$2.8917$

).

The application of the occupancy algorithm in the case of fuzzy parameters is shown in Figure 12. The analysis follows the same philosophy as previously described. It is quite interesting, that fuzzy parameters increase the uncertainty in the model, and the occupancy estimation algorithm results in overestimated or underestimated vectors. However, it should be clearly stated that this is not a weak point of the algorithm, but a direct result from the maximum likelihood nature of it.

3.2. Multi Agent Control Framework

As far as the MACS is concerned, the simulation time has been set to one year with a simulation step time of 0.001 h. The data, regarding the ambient temperature, the horizontal diffused illuminance, and the solar irradiance are acquired from the database of the Institute for Environmental Research and Sustainable Development (IERSD), National Observatory of Athens, and they are associated with the year 2019, for Athens, Greece with a sample time of 1 h. For the outdoor CO₂ concentration, the constant value of 359.4813 ppm that has been estimated, is used. The set point for the heating/cooling system has been set to 23 °C, when the outdoor temperature falls under 21 °C, and 26 °C when the outdoor temperature exceeds 29 °C. The set point for the indoor illuminance has been set to 700 $lux$ and the set point for the CO₂ concentration has been set to 380 ppm.

All the sub-models lighting, thermal, CO₂ concentration models and the algorithms of the MACS have been implemented in MATLAB/Simulink, each sub-model and each agent is an embedded MATLAB function with multiple inputs and outputs. The relation between inputs and outputs are described by the corresponding algorithms and equations. For example, the CO₂ concentration sub-model is a two input (number of occupants and ventilation rate) and one output system (CO₂ concentration).

The following figures illustrate the outdoor temperature, the indoor temperature, the control signal of the heating/cooling system and the occupancy of the building. After the extensive exploration, the indoor temperature remains almost constant at 23 °C for the winter and at 26 °C for the summer, during the hours when occupants are present inside the building. Figure 13 and Figure 14 are associated with one random day of winter and one random day of summer, respectively. When the occupancy level of the building is zero, the control signal of the agent is zero and consequently the indoor temperature is lower from the set point for the winter and above the set point for the summer. When occupants are entering the building the agent starts to operate again and the indoor temperature is reaching the set point with small oscillations.

A similar philosophy follows $A{G}_L$. At first, the corresponding agent visits new states and mainly explores the state-action space. Following the exploration, the indoor illuminance remains almost constant at 700 $lux$, at the hours of existing occupancy. Figure 15 depicts the outdoor horizontal diffused illuminance, the total indoor illuminance, the control signal, the transparency of the electro-chromic window, the provided illuminance by the artificial lighting system and the number of people inside the building, for a random day. When the occupancy of the building is zero, the control signal of the agent is zero, and consequently the indoor illuminance is proportional to the outdoor horizontal illuminance. When occupants enter the building, the agent starts to operate again and the indoor horizontal illuminance is reaching the set point with small oscillations, by either changing the transparency of the window or by changing the illuminance of the artificial lighting.

Figure 16 illustrate the occupancy of the building, the control signal of the agent and the indoor CO₂ concentration, for one random day. The agent explores in the beginning and after the exploration, the indoor concentration of the CO₂ remains lower than the value of 800 ppm. The ventilation system operates when there are people inside the building, and either reduces or stop its operation, when the occupancy level vanishes, and consequently the concentration of the CO₂ is reduced.

The overall energy consumption of the three systems: heating/cooling (H/C) system, artificial lighting (A/L) system and ventilation (V) system for one year period is computed as 1663 kWh in relation to the much larger value 3770 kWh, which is the overall consumption of these systems for the same year without taking into account the occupancy of the building: the heating/cooling system, as well as the lighting system, operate only for hours where ${\widehat{n}}_{occ}\ne 0$. This concludes an overall reduction of $55.9\%$.

The consumptions of each system, with and without regarding the occupancy of the building, are presented in

Table 5. Energy consumption of the three sub systems in numbers.

	MACS Taking into Account Occupancy	MACS without Taking into Account Occupancy
Heating / Cooling System	747 kWh	1478 kWh
Artificial Lighting System	414 kWh	1780 kWh
Ventilation System	502 kWh	502 kWh
Total	1663 kWh	3770 kWh

and graphically illustrated in Figure 17.

The overall comfort index $C_i$ arises from three indexes; the thermal comfort index $T_i$, the visual comfort index $L_i$, and the indoor air quality index $A_i$. The trapezoidal membership functions, that calculate the respecting indexes according to the indoor temperature, the indoor horizontal illuminance and the CO₂ concentration, are depicted in Figure 18.

The overall comfort index is computed by the following equation:

$$\begin{array}{c}\hfill C_i\left(x\right)=\left\{\begin{array}{cc}0.4T_i+0.4L_i+0.2A_i,\hfill & {\widehat{n}}_{occ}\left(t\right)>0\hfill \\ 0,\hfill & \mathrm{otherwise}\hfill \end{array}\right.\end{array}$$

(39)

Figure 19 illustrates $C_i$ and the number of occupants for a random day. $C_i$ is zero when no occupants are inside the building and equals, for most of the time, to one when there are occupants inside the building.

The MACS performance, pertaining to the case of fuzzy parameters, is presented in the following figures. Figure 20 and Figure 21, illustrate the outdoor temperature, the indoor temperature, the control signal, and the occupancy level estimated when using $\underline{\widehat{\theta }}$ and $\overline{\widehat{\theta }}$, respectively.

Figure 22 depicts the MACS results, associated with the occupancy estimation, pertaining to $\underline{\widehat{\theta }}$. The outdoor illuminance, the total indoor illuminance, the control signal, the transparency of the window, the artificial lighting and the occupancy level, are illustrated. The results when using the $\overline{\widehat{\theta }}$-based estimated occupancy vector, are given in Figure 23.

Figure 24 and Figure 25, illustrate the occupancy of the building, the control signal of the agent and the indoor CO₂ concentration, for one random day for the different occupancy vectors.

The total consumption, as well as the consumptions of the heating/cooling system, the artificial lighting system and the ventilation system for two different occupancy profiles, are presented in

Table 6. Energy consumption of the three sub systems in numbers, for the MACS using the occupancy profiles estimated using the fuzzy parameter.

	MACS Taking into Account Occupancy, Estimated Using $\underline{\widehat{\theta }}$	MACS Taking into Account Occupancy, Estimated Using $\overline{\widehat{\theta }}$
Heating / Cooling System	740 kWh	614 kWh
Artificial Lighting System	414 kWh	414 kWh
Ventilation System	379 kWh	500 kWh
Total	1533 kWh	1528 kWh

. In the first profile, the overall consumption is 1533 kWh, meaning a reduction of 130 kWh (Comparison with the first column of Table 5). This reduction is a result of the ventilation system. The agent of the ventilation system learns a policy that reduces the consumption of the actuator but, as it is illustrated in Figure 24, the indoor air of the building, is of lower quality, as the CO₂ concentration reaches values up to 900 ppm. Small reduction is also observed to the heating/cooling system (Figure 20), but it is justified by the increased number of the occupants: More occupants are heating the thermal space, due to the occupant’s heating gains.

In the second profile, the overall consumptions is 1528 kWh, meaning a reduction of 135 kWh (Comparison with the first column of Table 5). In this case, the reduction is a result from the heating/cooling system (Figure 21). The agent of this system learns a better policy that reduces the consumption of the system, without particularly reducing the thermal comfort. A small reduction of 2 kWh is observed to the ventilation system (Figure 25), which is justified by the lower number of the occupants. The consumption of the lighting system, as it is illustrated in Figure 22 and Figure 23, remains the same, since the lighting system is not affected by the occupant number.

4. Conclusions

A MACS with a parameter and occupancy estimation algorithm has been studied, with respect to the control performance and maintenance of the building’s comfort levels. A stochastic process models the CO₂ generation and decay, by including the unmodelled dynamics. Given CO₂ concentration data, the parameters of the model, as well as the occupancy level, are estimated. Maximum likelihood methodology follows the estimation. In order to induce further uncertainty in the assumed model, the infiltration and carbon dioxide concentration are treated using fuzzy numbers, and the effect of the occupancy estimation algorithm in the new model is studied. By using two boundary cases of the $\alpha$-cut method, for the fuzzy parameters, two occupancy estimations are produced. These occupancy estimated profiles can be interpreted to describe an uncertainty bound on the occupancy level which stems from the parametric uncertainty on the infiltration rate and the carbon dioxide generation rate per occupant.

The simulation results highlight the successful application of the independent learners approach of the MACS. This control framework succeeds stability over the desired set points for all the control parameters, in a fully decentralised framework. Each agent acts independently, thus a possible failure in one agent will not lead to a failure of the total control system. The enhanced approach, which takes into account the occupancy of the building, leads to even better performance because not only the set points of the parameters are achieved, but also an energy reduction of about $55.9\%$ is observed. This reduction is a result of the integration of the estimated occupancy vector, which dictates which time periods the thermal and lighting agents will function. Therefore, a waste of energy is avoided for the time periods where there are no present occupants in the building. This approach is also tested for different occupancy profiles arising by the fuzzy parameters with similar successful results.

In conventional building control techniques the control algorithm receives immediate feedback and generates a control action. In the proposed MACS for building control, the control algorithm receives a delayed feedback, which results each agent to incorporate previous knowledge. Thus, an agent applies the optimum policy depending on the long term reward.

The over/under estimated occupancy level, when considering fuzzy parameters, indicate that further research should focus on the field of parametric uncertainty and on fuzzy stochastic differential equations. The incorporation of uncertainty bounds on occupancy estimation, could be included in the MACS control design to further reduce the energy consumption. Future research directions will also focus on the occupancy prediction problem, as well as the development of evolutionary based algorithms for occupancy estimation based on CO₂ concentration data.

A quiet interesting topic that we surely plan to investigate involves more complex approaches of cooperative MAS, which combine the indoor environmental quality parameters which, in combination with artificial general intelligence (Deep AI), may lead to further improvement of the proposed system. Furthermore, the impact of the ventilation system to the indoor temperature will be introduced into the thermal model of the building for providing more accurate results. Finally, the effect of ventilation systems, with and without heat recovery, on thermal comfort and energy consumption, will be studied and compared.