Markov Chain Monte Carlo Based Energy Use Behaviors Prediction of Office Occupants

Prediction of energy use behaviors is a necessary prerequisite for designing personalized and scalable energy efficiency programs. The energy use behaviors of office occupants are different from those of residential occupants and have not yet been studied as intensively as residential occupants. This paper proposes a method based on Markov chain Monte Carlo (MCMC) to predict the energy use behaviors of office occupants. Firstly, an indoor electrical Internet of Things system (IEIoTS) for the office scenario is developed to collect the switching state time series data of selected user electrical equipment (desktop computer, water dispenser, light) and the historical environment parameters. Then, the Metropolis–Hastings (MH) algorithm is used to sample and obtain the optimal solution of the parameters for the office occupants’ behavior function, the model of which includes the energy action model, energy working hours model, and air-conditioner energy use behavior model. Finally, comparative experiments are carried out to evaluate the performance of the proposed method. The experimental results show that while the mean value performs similarly in estimating the energy use model, the proposed method outperforms the Maximum Likelihood Estimation (MLE) method on uncertainty quantification with relatively narrower confidence intervals.


Introduction
Building energy consumption has increased dramatically due to population growth, increased demand for building functions, and global climate change in recent decades [1]. The proportion of building energy consumption in China has increased to about 35 percent, in which office building energy consumption occupies an important part [2]. A large number of studies show that building energy consumption is not only affected by meteorological parameters, building shape, envelope structure, window-to-wall ratio, etc., it is also related to the behaviors of occupants [3]. By mining the hidden patterns from office occupants' equipment usage data and establishing a corresponding personnel energy use behavior model, the most preferred equipment control scheme according to the needs of users can be recommended and configured, and the experience of the office occupants can be improved and accurate demand response strategies for office buildings designed accordingly.
Many scholars have made considerable contributions in the field of personal behavior prediction and have obtained some excellent research results, some of which were reviewed in [4]. The traditional methods of personal behavior prediction mainly include naive Bayes (NB) [5], support vector machine (SVM) [6], k-nearest neighbor (KNN) [7], the hidden Markov model (HMM) [8], and convolutional neural networks (CNN) [9]. Withanage et al. [10] used the depth cuboid similarity feature (DCSF) algorithm to identify and predict the behavior characteristics of the elderly through video capture of personnel behavior information. Although this method has high accuracy, it is easy to invade experiments are carried out and the results are discussed in Section 4. Finally, conclusions are drawn in Section 5.

Methodology
The Markov chain Monte Carlo (MCMC) method has been widely used in machine learning, deep learning, natural language processing, and others [25]. The MCMC method consists of two MCs. The first MC is Monte Carlo sampling, and the second MC is the Markov chain. It is essentially a method of repeatedly drawing the random values of the assumed parameters according to the current value until it converges to the the real distribution. The samples of each value are randomly sampled, while the selection of values is limited only by the current state and the prior distribution of the assumed parameters [26][27][28][29].

Markov Chain
In the Markov chain, the probability of state transition from the current state to the next state depends only on its current state, independent of the previous states [30]. Suppose · · · , x i−3 , x i−2 , x i−1 , x i , · · · are the states in a Markov chain. The transition probability can be described as by: where P is the probability of the occurrence of the event; P(x i |x i−1 ) is the probability of transitioning to the state x i under the condition of x i−1 . A Markov chain that satisfies the detailed balance condition can converge to a stationary distribution, and it is given by: where π(x) is the stationary distribution of the system; π(x (i) ) represents the probability value corresponding to the current state; T(x) is the transfer matrix; T(x (i) |x (i−1) ) is the transition probability from state i − 1 to i. The existence of this property makes the Markov chain have a recurrence; therefore, the state of the Markov chain can return to any state with a probability greater than zero when the value of the variable iis large enough [31]. A Markov chain that satisfies the detailed balance condition can be converted from the i − 1 state to i for any two states and without any loss in the conversion. However, under normal circumstances, the two state transitions in the Markov chain are not equal and do not satisfy the detailed balance condition. How to construct the transfer matrix so that the stationary distribution is the target distribution is the key of Markov chain sampling.

Metropolis-Hastings Algorithm
The Metropolis-Hastings (MH) algorithm is the cornerstone of Markov chain Monte Carlo sampling [32]. In order to make the Markov chain meet the detailed balance condition, the acceptance rate r(x (i−1) , x (i) ) is introduced to make the original Markov chain jump from state i − 1 to i, which is received with the probability of r(x (i−1) , x (i) ) [33]. The original Markov chain is random. Assuming there is a stationary distribution of π(x), x is a sample obtained from the Markov chain, and the ratio of π(x (i) ) and π(x (i−1) ) can be calculated by Equation (2), while parameters x (i) and x (i−1) can be selected according to a certain probability. The parameters by multiple iteration sampling should be met π(x) [34,35]. The iterative steps of the MH algorithm are as follows: • The MH algorithm takes random values in the parameter space as the starting parameter x (0) . The random parameter x (1) is generated from the distribution function of parameters, and the probability density of the current parameter is calculated.

•
First, a random number µ t is extracted from the uniform distribution of [0, 1], and then, we determine whether to keep the current parameter according to whether the probability density ratio between the current parameter x (i) and the starting parameter x (i−1) is greater than µ t . If the reception probability of the current parameter is greater than µ t , µ t ≤ r(x (i−1) , x (i) ), the state is accepted, and it needs to be updated [36]. The acceptance probability r(x (i−1) , x (i) ) of MH can be described as: It is obvious that the acceptance probability r(x (i−1) , x (i) ) value is between [0, 1], and its purpose is accelerating the convergence of the Markov chain to the detailed balance condition.

Proposed Method
The energy use behavior of the office occupant has a certain regularity, which is related to the working time and thermal comfort of users. The daily on and off states and working hours of electrical equipment are related to personal energy use behavior habits, and the energy use behavior of air-conditioner is related to users' thermal comfort. The purpose of this paper is to predict the office occupant energy use behavior by mining the user's equipment usage data. This method regards the energy use behavior model as a time or temperature function and uses the MCMC method to fit the posterior probability distribution of the energy use action model, energy use working hours model, and air conditioning energy use behavior model. The prediction process of personnel energy use behavior can be roughly divided into three sections: data acquisition, data preprocessing, and model training.

Data Acquisition
In order to obtain the usage data of electrical equipment in the office, we build an indoor electrical Internet of Things system (IEIoTS), which can obtain the electrical equipment's state data and environmental parameters by installing smart sockets, smart switches, temperature and humidity sensors, and electrical parameter detection modules. The architecture of this system is shown in Figure 1.
As can be seen from Figure 1, this system includes electrical equipment, a gateway, and a server from the bottom up. The system was deployed in an office building of Shandong Jianzhu University, and the collected data were stored and managed in units of rooms. The electrical equipment of each room included smart sockets, smart switches, and other electrical equipment. The on-off state data of electrical equipment could be acquired by smart sockets and switches, and temperature data were acquired by temperature and humidity sensors. The gateway was the core equipment of the indoor electrical Internet of Things system, which could not only complete the two way transmission of control instructions and state data, but also complete the protocol conversion, and so on [37]. At the top of the system, there was a cloud server to manage the data in the building and store the effective information from the sensors in each room. The data of electrical equipment was uploaded to the database in the cloud server by the gateway for analysis and storage. The TCP/IP protocol was adopted, which could reliably transmit the data collected in the remote environment to the database in real time.
The dataset used in this paper was generated from the IEIoTS deployed in a laboratory building of Shandong Jianzhu University and an occupant in the office. The on/off state data of three pieces of equipment (water dispenser, desktop computer, and light) and the historical data of environmental parameters in the lab were collected and recorded as the original data. The sampling frequency was 1 sample per minute. The original data of half a year were used in this paper, and the total number of the data was 65,880, during which 98% were used as the training data and the remaining 2% as the test data.

Data Preprocessing
The preprocessing of the original data mainly included: time interception, removing the zero vector, excluding outliers, data interception, etc. The data preprocessing process is shown in Figure 2. (1) Time interception: The data should strictly require daily data starting from 00:00, with one minute as the sampling period and yyyy/m/d h:mm as the time format for date processing.
One-thousand-four-hundred-forty records could be obtained in one day. It was necessary to determine whether the data were complete. If there was a missing value, we assigned the previous state. (2) Removing the zero vector: Removing the state data of the powered device that was not working for a whole day. (3) Excluding outliers: For data outliers that occasionally occurred in the original data, outliers were detected and removed by the Laida algorithm [38].
(4) Data interception: By analyzing the time period of the equipment action daily, the state data during the device on and off time periods were respectively intercepted.

Model Training
The energy use action model refers to the probability that the electrical equipment enters the power-on and power-off state within a certain temperature; the energy use working hours model refers to the probability of the electrical equipment duration of energy use; and the energy use behavior model of air-conditioner refers to the probability that the air-conditioner enters into the cooling and heating state at a given temperature. Since time and temperature are continuous variables, which are difficult to determine for the entire posterior distribution, the MCMC was used to study the approximate distribution and obtain the fitting curve of the data distribution.
The energy use action model and the energy use working hours model reflect the time rule of user's energy use behavior. The air-conditioner energy use behavior model reflects user's thermal comfort, which depends on user's comfortable temperature range. The process of model training was as follows, taking the energy use action model as an example to elaborate in detail. At first, we needed to observe the electrical equipment's state distribution based on the time series to determine the fitting function with the parameters that needed to be trained, set the number of samples, and give the initial value to the parameters for the sampling. Secondly, we assumed the a priori distribution of the parameters, and MH random sampling was carried out by searching the distribution space of the parameters. Based on Markov characteristics, new values were randomly assigned to the parameters according to the current state. Thirdly, we checked the new random values; if they satisfied the MH sampling theorem, then we accepted them as a new current state; if not, we rejected and returned them to the previous state. Lastly, we repeated the previous two steps for the specified number of iterations, until it reached the detailed condition of a Markov chain. This method could construct the energy use behavior model based on different time ranges (daily, weekly, annually, etc.), and the calibration frequency of the model needed to be adjusted accordingly. We took the construction of the day based energy use behavior model as an example to evaluate the performance of the proposed method in this paper. The flowchart of energy use behavior prediction is shown in Figure 3. The algorithm returned all the values generated by the parameters, and the mean of these values was used as the most likely final value for the parameters in the model. We used the Bernoulli function to calculate the probability of the equipment action time. The MCMC returned an approximation of the distribution instead of the accurate results.

Action Model
(1) Function model selection: This study took the desktop computer as an example to introduce the algorithm in detail. First, a fitting function was required to model the posterior probability distribution of the energy use action behaviors (on and off) for the desktop computer. The distribution of the historical state data of the desktop computer was observed. The power-on time period was selected as 6:00-9:00, and the power-off time period was selected as 14:00-17:00. The action time probability distribution is shown in Figure 4.
As shown in the Figure 4, it can be concluded that users usually powered on and powered off after 7:30 and 15:00, respectively, and the model was created to simulate the transition process of the user between power-on and power-off in the form of a probability. After many experiments, the logistic function was chosen as the posterior probability distribution model of the on-off action. It can be described as: where t represents the action time. The parameter β affects the direction and steepness of the curve; the parameter α affects the position. β and α are the two model parameters that need to be trained in the MCMC process jointly. The logistic function with different parameters is shown in Figure 5. The logistic function satisfied the experimental data, because the probability of the power-on and power-off of the desktop computer would gradually change, which could reflect the change of the on-off action model.   (2) MH sampling: The procedure of the MH algorithm was randomly sampling from the probability distribution model of the parameters. Each time, the sampling values of the parameters (β and α) of the function were changed, and the response was observed. The mean value of all sample values was taken to obtain the optimal approximation.
First, in order to plot the random values of the parameters β and α, it was necessary to assume the prior distribution of the parameters. In this paper, a 2D Gaussian distribution was selected as the prior distribution of the parameters. The two edge distributions of the 2D Gaussian distribution were all in the form of a univariate Gaussian distribution, and this can be expressed as: where µ is the mean of the normal distribution; τ is the accuracy; f (x) is Gaussian distribution function. The accuracy (τ) determines the concentration of the distribution; the mean (µ) determines the position of the distribution. The higher the accuracy, the more centralized the presentation data and the smaller the change. The mean can be positive or negative, and the precision must be positive. In order to broaden the search area of the acquired β and α values, choose µ = 0 and τ = 0.01. The parameter space of the prior distribution is shown in Figure 6. It can be seen from Figure 6 that not every point in the parameter distribution space could be attempted, but the most fitting model could be created by random sampling from a higher probability region (the red part in Figure 6). In order to make the test results more accurate, a step size as long as possible should be chosen. In this paper, the step size was 5000. The MH algorithm evaluated the random values of β and α each time. If the MH sampling condition was satisfied, the current values of β and α would be accepted, and the state was transferred. If not, the values of β and α were random values of the previous time. The sampling trajectory of the parameters β and α of the desktop computer power-on model are shown in Figure 7.  The Bernoulli variable was adopted to model the on or off action of equipment, where the power-on was set to 1 and the power-off was set to 0 [39]. The Bernoulli variable of the desktop computer state data depended on the time, which could be defined by the logistic function. It can be described as: where p(t i ) is the logistic function with time as the argument, representing the probability of action at a certain time. Therefore, the probability of power-on and power-off was converted by Equation (4), and it can be written in the form: The goal of MCMC was to find the optimal values of parameters β and α based on the data from the assumption of normal prior distribution. In order to find the most optimal posterior probability distribution of the turn-on and turn-off action model, the mean values of β and α samples were taken as the optimal parameters of desktop computer, water dispenser, and light. The probability distribution of the energy use action model is shown in Figure 8.

Working Hours Model
According to the characteristics of the data distribution, we adopted the skewed distribution as a posterior distribution of the working hours model. The skewed distribution is the frequency distribution of deviating from the symmetrical value of variables, which can be expressed as: where φ(x; µ, σ) and Φ(x; µ, σ) are the density function and distribution function of the standard normal distribution, respectively [40,41]. The MCMC method was used to obtain the optimal values of the mean values µ, variance σ 2 , and skewness λ in the skewed normal distribution, as well as fitting the maximum probability value of the user electrical equipment's used time. The energy use working hours models of the desktop computer, light, and water dispenser were established, which are shown in Figure 9.

Air-Conditioner Energy Use Behavior Model
The air-conditioner can effectively improve the indoor thermal environment and meet the desired thermal comfort of office occupants. Through the collection of the environmental parameters (temperature, humidity) of the office air outlet, the on-off state and the set temperature of the air-conditioner could be analyzed, then the air-conditioner energy use behaviors of office occupants could be predicted.
The MCMC algorithm was used to build the energy use behavior model of the manual on and off air-conditioner and predict the user's air-conditioner use behavior based on temperature data. The energy use behavior model of the air-conditioner refers to the probability that the air-conditioner enters into the cooling or heating state in a given temperature. The air-conditioner energy use behavior model included air-conditioner heating-on, heating-off, cooling-on, and cooling-off action model. The training process of the air-conditioner energy use behavior model was similar to the example of the energy use action model in Section 3.3.1. The difference was that the air-conditioner energy use behavior model was based on temperature, while the energy use action model was based on time series of the equipment state. Specifically, the steps of MCMC algorithm for the air-conditioner energy use behavior model are as follows. At first, we needed to observe the air-conditioner state distribution based on the temperature to determine the fitting function with the parameters that needed to be trained; the logistic function was selected similarly, and we set the number of samples and gave the initial value to the parameters for the sampling. Secondly, we assumed the a priori distribution of the parameters was 2D Gaussian, and MH random sampling was carried out by searching the distribution space of the parameters. Based on Markov characteristics, new values were randomly assigned to the parameters according to the current state. Thirdly, we checked the new random values; if they satisfied the MH sampling theorem, then we accepted them as a new current state; if not, we rejected and returned them to the previous state. Lastly, we repeated the previous two steps for the specified number of iterations, until it reached the detailed condition of a Markov chain, which was the probability distributions of the air-conditioner energy use behavior.
The experimental period was from 22 January 2019 to 4 June 2019. The highest, lowest, and average outdoor temperatures of Jinan in each month during the experiment are shown in Figure 10. The temperatures in May and June were higher, with the average temperatures of 24 • C and 27.5 • C, respectively. The highest temperature in May was 38 • C; the lowest temperature in January was −7 • C. Figure 11 shows the temperature change of a user's office with the air-conditioner from 22 to 25 May 2019. It can be seen from Figure 11 that the temperature of the air outlet was about 10-15 • C, which was lower than the indoor temperature when the air-conditioner was running. Based on this, the working states of the air-conditioner could be judged.   Figure 11 shows several typical air-conditioner operation processes. Processes 2 and 4 show that after the user turned on the air-conditioner, the air outlet temperature dropped to about 17 • C, and with little fluctuation, so the air-conditioner compressor was in constant frequency operation. After the operation of Processes 1, 3, and 5 for a period of time, due to the frequency conversion of the compressor, the temperature was fluctuating continuously. Before the air-conditioner was started manually, the temperature was basically stable and high. Therefore, the above judgment method could distinguish the manual on and off condition of the air-conditioner from the automatic on and off process, so as to obtain the record of the manual on and off condition of the air-conditioner. Through the above MCMC algorithm, we could get the probability model of the on and off condition of the air-conditioner at a certain temperature, as shown in Figure 12c.

Evaluation Indices
The root mean squared error (RMSE) and mean absolute error (MAE) are usually used as the loss function of linear regression in machine learning. RMSE is the square root of the squared sum of the deviation between the observed value and the true value with the ratio of the observed times and is used to measure the deviation between the observed value and the true value. MAE is the average of absolute errors, which can better reflect the actual situation of the predicted value error. Therefore, RMSE and MAE were selected to verify the accuracy of the experiment, which can be described as: where m is the number of sample data; y i is the true value; y i represents the predicted value. y i and y i are both probability values, and the range is [0,1]. Using the above evaluation indices, the loss between the energy use action time model and the real value of each electrical equipment could be obtained, and the accuracy of the experiment could be evaluated. The MCMC method was used to obtain 5000 sample values, and the average of all sample values was taken as the final parameter of the model. The results of RMSE and MAE are shown in Table 1; the value of RMSE and MAE of the light power-on model was smaller than the desktop computer power-on model, so the light power-on model was more accurate than the desktop computer power-on model.

Comparison between MCMC Estimation and MLE
The Maximum Likelihood Estimate (MLE) is a statistical method based on the principle of maximum likelihood, which is the application of probability theory in statistics. It provides a way to estimate a model for a given observation. The steps of solving the maximum likelihood estimation value mainly include: writing out the likelihood function (θ) at first, taking the logarithm of the likelihood function, and finally, solving the maximum likelihood function estimation value θ. (θ) and θ can be described as: where x 1 , x 2 , ..., x N are samples; is the parameter vector. We compared the MCMC and MLE algorithm for the analysis of the energy use behavior models corresponding to the optimal parameters β and α.  Table 2, for the desktop computer power-on model; the model interval difference (d) fitted by the MLE algorithm and the MCMC algorithm was 1.6074 and 0.9780, respectively. The smaller the d was, the more the data in the confidence interval, and the better the fitting result. Therefore, the MCMC algorithm was better than the MLE algorithm. By comparing the mean (v), the 90% credible interval (L), and the difference of the interval (d) for the model, it could be obtained that while the mean value performed similarly in estimating the energy use model, the proposed method outperformed the Maximum Likelihood Estimation (MLE) algorithm on uncertainty quantification with relatively narrower confidence intervals.

Experimental Results and Analysis
In order to analyze energy use behavior intuitively, the probability distribution model of energy use behavior in Figure 8 was transformed into the probability density distribution model, as shown in Figure 12a,b.
As can be seen from the Figure 12a,b, the energy use behaviors had a certain regularity. In the morning, the user tended to turn on the light between 7:40 and 7:50 firstly, then turn on the desktop between 7:50 and 8:00, and finally, turn on the water dispenser between 8:00 and 8:10. In the afternoon, the user turned off the light between 15:10 and 15:20 at first, then turned off the desktop and the water dispenser between 15:20 and 15:30 PM. However, the analysis of the office occupant energy use behavior habits had a high time dependence and uncertainty. The energy use behavior of the office occupants can be deduced from Table 3.
It can be seen from Figure 12c that the average temperature of air-conditioner cooling-on was 27.5 • C, and the average off temperature was 23.3 • C; the average temperature of air-conditioner heating-on was 18.2 • C, and the average off temperature was 25.7 • C. The on and off behavior of the air-conditioner also reflected the thermal comfort of the users. When the indoor temperature in summer exceeded 27.5 • C, the thermal endurance of users exceeded the threshold value, and the air-conditioner was turned on for cooling; when the indoor temperature dropped to 23.3 • C, the users felt cold, and the air-conditioner was turned off at this time. Therefore, 23.3-27.5 • C could be regarded as the comfortable temperature range for the experimenters under the working condition of the air-conditioner in summer.

Conclusions
In this paper, a new method for predicting personal energy use behavior was proposed. Based on the data of electrical equipment used by an office occupant for half a year, the prediction of energy use behavior was realized by using the MCMC method. Firstly, the on-off state data of electrical equipment and temperature history data were obtained by building the IEIoTS for office application, and the experimental data were obtained by data preprocessing. Secondly, the sampling parameters were determined by selecting an a priori function model from observation data. Then, the optimal parameters were obtained through continuous sampling with the MH algorithm in MCMC, and the energy use action model, the energy use working hours model, and the air-conditioner energy use behavior model were fitted. Finally, comparative analysis of the experiments verified that the proposed method had high accuracy and stability. The method proposed in this paper was a generic statistical modeling framework. As long as the user behavior had certain statistical patterns, this method could be used for prediction. Therefore, the model could be applied to other types of buildings, such as residential buildings, industrial buildings, and so on.
In some cases, equipment states and performance data may not be available. For example, considering the user's privacy, we may not be able to access equipment states of residential buildings, but it is possible to infer the working state information of the equipment by using the aggregated data from the household electricity meters, and then, the method in this paper could be used to predict the user behaviors. Although the proposed method could obtain the rules of personal energy use behaviors, it could not predict the energy use behavior rules with a high resolution. In the future study, we will further improve the accuracy of the algorithm, for example by acquiring longer electrical equipment state data. In addition, we will design energy saving strategies on the basis of office occupants' energy use behaviors.

Conflicts of Interest:
The authors declare no conflict of interest.