Loading required pacakges.

suppressPackageStartupMessages(library(tidyverse)) # data processing
suppressPackageStartupMessages(library(forcats)) #data processing
suppressPackageStartupMessages(library(rlang))  #data processing
suppressPackageStartupMessages(library(readr)) #data processing
suppressPackageStartupMessages(library(lubridate)) #datetime date handling
suppressPackageStartupMessages(library(ggforce)) #visualization
suppressPackageStartupMessages(library(patchwork)) #visualization
suppressPackageStartupMessages(library(expm)) #fast matrix exponential 
suppressPackageStartupMessages(library(lhs))  #latin-hyper cube sampling 
suppressPackageStartupMessages(library(Rfast)) #fast matrix algebra
suppressPackageStartupMessages(library(parallel)) #parallel particle filter
suppressPackageStartupMessages(library(DEoptim)) #parallel particle filter
suppressPackageStartupMessages(library(kdevine)) #kernel density approximation for 4.4

1. Introduction

This is accompanies toy example code of section 3.2 in paper “Ham, S. woo, & Karava, P. (2020). Online model for unit-level heating and cooling energy use prediction in multi-family residential buildings for eco-feedback design. Submitted to Journal of Energy and Buildings.”. This is a simple demonstaration of Liu-West filter for a simple building gray-box model. We generated a synthetic data and applied Liu-West filter to see if the filter can be applicable for this problem. All code is written in R language. The main purpose of this document is to provide reproducible example.

The pacakge depedency of this code is managed by renv package. You can look at renv.lock file to see the required package. However, for the simplicity, just run following script on this rproject.

After clone the repository and run building_lw_filter.Rproj file (you must have Rstudio) with R>3.5.3.

git clone https://github.com/ecosang/building_lw_filter.git

In R console, run following script.

# Run this code in Rstudio
install.packages('renv',repos="https://cran.rstudio.com")
renv::equip() #install required software
renv::restore()

Technically, this installs all required packages, and you can reproduce all codes below. However, if this doesn’t work, please report issues.

All functions used in this code is in code/utility.R. Also, all generated data and trained model are stored in data folder.

2. Model description

2.1 A simple gray-box model with two states.

The gray-box model is composed of two states. Below figure shows R-C circuit diagram of this model.

All variables are listed below.

Variables

  • \(x_{\text{i}}\): indoor temperature state [\(^\circ\text{C}\)]
  • \(x_{\text{e}}\): envelope temperature state [\(^\circ\text{C}\)]
  • \(T_{\text{a}}\): outdoor air temperature [\(^\circ\text{C}\)]
  • \(y_{x_\text{i},t_k}\): measured temperature of state \(x_{\text{i}}\) [\(^\circ\text{C}\)]
  • \(R_{\text{ie}}\): thermal resistance between \(x_{\text{i}}\) and \(x_{\text{e}}\) [\(\text{K/W}\)]
  • \(R_{\text{ea}}\): thermal resistance between \(x_{\text{e}}\) and \(T_{\text{a}}\) [\(\text{K/W}\)]
  • \(C_{\text{e}}\): thermal capacitance of \(x_{\text{e}}\) [\(\text{J/K}\)]
  • \(C_{\text{i}}\): thermal capacitance of \(x_{\text{i}}\) [\(\text{J/K}\)]
  • \(\dot{Q}_{\text{hc}}\): heating or cooling flow rate [\(\text{W}\)]
  • \(t\) and \(t_k\) are time in seconds and its discretized time, respectively.
  • \(d\omega/dt\) is standard weiner process.
  • \(\sigma_\text{i}^2\) and \(\sigma_\text{e}^2\) are process noise variance of states \(x_{\text{i}}\) and \(x_{\text{e}}\), repsectively.
  • \(\sigma_\text{d,y}^2\) is observational noise variance of measurement \(y\).

2.2 Governing differential equations

The thermal dynamics are expressed with below system of differntial equations. Here, \(R\) is thermal resistance and \(C\) is thermal capacitance.

\[ dx_{\text{i},t} =\left( \frac{1}{R_{\text{ie}} C_{\text{i}}}(x_{\text{e},t}-x_{\text{i},t}) +\frac{1}{C_{\text{i}}} \dot{Q}_{\text{hc},t} \right)dt+\sigma_i d\omega_{\text{i},t} \]

\[dx_{\text{e},t}=\left( \frac{1}{R_{\text{ie}} C_{\text{e}}}(x_{\text{i},t}-x_{\text{e},t})+ \frac{1}{R_{\text{ea}} C_{\text{e}}}+(T_{\text{out},t}-x_{\text{e},t}) \right)dt+\sigma_ed\omega_{\text{e}}\]

\[y_{x_\text{i},t_k}=x_{\text{i},t_k}+\varepsilon_{y,t_k} \text{ where } \,\varepsilon_{y,t_k}\sim\mathcal{N}(0,\sigma_{\text{d},y})\]

2.3 System equations in a probablistic format

The system of matrix can be expressed in matrix form.

\[\begin{bmatrix} \dot{x}_{\text{i},t} \\ \dot{x}_{\text{e},t} \end{bmatrix}= \underbrace{ \begin{bmatrix} -\frac{1}{R_\text{ie}C_\text{i}} & \frac{1}{R_\text{ie}C_\text{i}} \\ \frac{1}{R_\text{ie}C_\text{e}} & -\frac{1}{R_\text{ie}C_\text{e}}-\frac{1}{R_\text{ea}C_\text{e}} \end{bmatrix} }_{\textbf{A}} \begin{bmatrix} x_{\text{i},t} \\ x_{\text{e},t} \end{bmatrix} + \underbrace{\begin{bmatrix} 0 & \frac{1}{C_\text{i}} \\ \frac{1}{R_\text{ea}C_\text{e}} & 0 \end{bmatrix} }_{\textbf{B}} \begin{bmatrix} T_{\text{out},t} \\ \dot{Q}_{\text{hc},t} \end{bmatrix}+ \begin{bmatrix} \sigma_{i}\dot{\boldsymbol{\omega}}_{\text{i},t} \\ \sigma_{e}\dot{\boldsymbol{\omega}}_{\text{e},t} \end{bmatrix}\]

\[ y_{x_\text{i},t_k}=\underbrace{\begin{bmatrix}1&0\end{bmatrix}}_{\textbf{C}}\begin{bmatrix}x_{\text{i},t_k}\\x_{\text{e},t_k}\end{bmatrix}+\varepsilon_{y,t_k} \]

Since this system is a linear gaussian model, it can be discretized without integration (see section 3.1 of paper) (or this link). Here, subscript \(_d\) indicates discretization. \(\textbf{x}_{t_k}=\left[x_{\text{i},t_k},\,x_{\text{e},t_k}\right]^\intercal\), \(\textbf{u}_{t_k}=\left[T_{\text{a}},\dot{Q}_{\text{hc}}\right]^\intercal\), and \(\boldsymbol{\theta}=\{R_{\text{ie}}, R_{\text{ea}}, C_{\text{e}}, C_{\text{i}}, \sigma_\text{i}, \sigma_\text{e}\}\)

\[P(\textbf{x}_{t_k+1}|\textbf{x}_{t_k})=\mathcal{N}(\textbf{x}_{t_k+1}|f_d\left(\textbf{x}_{t_k},\textbf{u}_{t_k},\boldsymbol{\theta}\right),\boldsymbol{\sigma}_{\text{d},x})\]

\[P(y_{x_\text{i},t_k}|\textbf{x}_{t_k})=\mathcal{N}(y_{x_\text{i},t_k}|g_d\left(\textbf{x}_{t_k},\textbf{u}_{t_k},\boldsymbol{\theta}\right),\sigma_{\text{d},y})\]

\[f_d\left(\textbf{x}_{t_k},\textbf{u}_{t_k},\boldsymbol{\theta}\right)=\textbf{A}_\text{d}\textbf{x}_{t_k}+\textbf{B}_\text{d}\textbf{u}_{t_k}\,\text{ and }\,g_d\left(\textbf{x}_{t_k},\boldsymbol{\theta}\right)=\textbf{C}_\text{d}\textbf{x}_{t_k}\] where \(\textbf{A}_\text{d}\), \(\textbf{B}_\text{d}\), and \(\textbf{C}_\text{d}\) are discretized parameters of \(\textbf{A}\), \(\textbf{B}\), and \(\textbf{C}\) matrix given above.

3. Synthetic Data

3.1 True parameter

Based on the above model, we assigned true parameter values to generate synthetica data. From the synthetic data, we apply Liu-West filter, and the final posterior distribution of parameters will be compared to those true parameter values.

\(x_{\text{i},0}\) \(x_{\text{e},0}\) \(C_{\text{i}}\) \(C_{\text{e}}\) \(R_{\text{ie}}\) \(R_{\text{ea}}\)
true 21.0 15.0 500000 6000000 1/55 1/55

To assign, true values, we assume following properties.

  • Total surface are of exterior wall and window are \(200 \text{m}^2\) and \(50 \text{m}^2\), respectively.
  • U-value of exterior wall and window are \(0.2 \text{W/m}^2\text{-K}\) and \(1.4 \text{W/m}^2\text{-K}\), respectively.
  • The total R-value is \(\frac{1}{UA_{\text{ext-wall}}+UA_{\text{window}}}=\frac{1}{110}\). Thus, we split \(1/110\) by \(1/55\) and \(1/55\) because only the sum of \(UA_{\text{ext-wall}}+UA_{\text{window}}\) is important.
  • Assuming all envelope mass is concerete. Specific heat (\(C_p\): 1000 \(\text{J/kg-K}\)) and density (\(\rho\): \(2400\text{ kg/m}^3\)). Volume with \(0.1\text{m}\) thickness: \(250\text{m}^2\times0.1\text{m}=25\text{m}^3\).
  • The exterior wall heat capacitance (\(C_{\text{e}}\)): \(C_p\times\rho\times V\approx6000000\text{ J/K}\)
  • Assuming air specific heat is \(1\text{J/K-m}^3\) and volume \(1000\text{ m}^3\), the indoor heat capacitance (\(C_{\text{i}}\)) is \(500000\text{ J/K}\).

Noise parameters

\(\sigma_{\text{d},y}\) \(\sigma_{x_i}\) \(\sigma_{x_e}\)
\(0.25\) \(0.01/\sqrt{900}\) \(0.01/\sqrt{900}\)

Here, \(\sigma_{\text{d},y}\) is set to \(0.25\) becuase the sensor measurement accuracy is \(\pm{0.5} ^\circ\text{C}\). Specifically, \(\mathcal{N}(0,0.25)\) generates data about \(\left[-0.5,0.5\right]\) thinking 95% quantiles.

The continuous process noise \(\sigma_i\) and \(\sigma_e\) are set to 0.25/30, respectively. When we discretize the continuous system, the discretized standard deviation is approximately an order of \(\sqrt{\Delta t}\). Therefore, assuming our process noise is \(0.25\) in a discrete time-scale, it is approximately \(0.01/\sqrt{\Delta t}\). Our dataset has \(900\text{s}\) time-scale. Therefore, they are set to \(0.01/\sqrt{900}\)

In the model, we will use fixed value of \(\sigma_{\text{d},y}\) because the measurement noise value can be obtained from the sensor information. Actually, this is helpful to stabilize Liu-West filter operation.

3.2 Read input data and generate synthetic data

Based on input data (\(\textbf{u}_{1:t_K}\)), we will generate synthetic observation data by using a stochastic simulation.

# load functions
source("../code/utility.r")
# load input u data.
dat=readr::read_csv("../data/syn_data.csv")
dt=dat$time[2]-dat$time[1] #discrete time
# define all true parameters
x0_true=c(21.0, 15.0) # initial states

#           Ci     Ce        Rie       Rea  #sigma_x_i #sigma_x_e
par_true=c(500000, 6000000,  1/55,     1/55, 0.01/sqrt(dt),      0.01/sqrt(dt) )
dims=c(2,2,1) # dimension of x, u, y, which is [nx, nu, ny]
nx=dims[1] # x dim
nu=dims[2] # u dim
ny=dims[3] # y dim
t_K=dim(dat)[1] #total time time is (1:t_K)

# create u_matrix [nu x t_K], (T_a, Qhc)
u_mat=rbind(t_a=dat$t_a,q_hc=dat$Q_hc) #input u matrix [nu x NT].
# weplit data into two parts to separate prior generation part, filter part.
index_split=ceiling(t_K/2)

Create data synthetic dataset. The data is created once and stored into data/synthetic_data.rds fiel. Then, it is loaded for next time because we generate the data with random noise from \(\sigma_{x_\text{i}}\), \(\sigma_{x_\text{e}}\), and \(\sigma_{\text{d},y}\). Here the dataset are splitted into two parts. The first part (prior_data) is used to create priors. The second part (filter_data) is used to used particle filtering process.

## not run this code. Data is generated once, and we reuse it.
# create discretized system matrix
dsys=discretize_system(par=par_true,dims=dims,dt=dt)
# create data set with stochastic process
simulated_data=create_data_set(dsys=dsys,u_mat=u_mat,x0=x0_true,dims=dims,dt=dt,seed_num=1234,stochastic=T)
y_t_i=simulated_data$y_mat[1,] # observation data

# split data by train_data, test_data, all_data
prior_data=list(u_mat=u_mat[,1:index_split],
                y_mat=y_t_i[1:index_split])
filter_data=list(u_mat=u_mat[,(index_split+1):t_K],
                y_mat=y_t_i[(index_split+1):t_K])
all_data=list(u_mat=u_mat,
                y_mat=y_t_i)

# store synthetic data.
write_rds(list(prior_data=prior_data,filter_data=filter_data,all_data=all_data,index_split=index_split),
          paste0("../data/synthetic_data.rds"))

To see how it looks like, visualize the data with non-stochastic simulation.

## [1] "non-stochastic"

4. Create priors

As described in section 3.4 Model initialization (Prior generation) in the paper, we created priors to initialize Liu-West filter.

Fisrt, we will use optimization to generate various set of parameter solutions on zero-process noise model.

Define parameter lower and upper bounds that create priors.

4.1 Upper and lower bounds of states and parameter

\(T_{\text{i},0}\) \(T_{\text{e},0}\) \(C_{\text{i}}\) \(C_{\text{e}}\) \(R_{\text{ie}}\) \(R_{\text{ea}}\)
true 21.0 15.0 500000 6000000 1/55 1/55
min 0 0 100000 1000000 1e-9 1e-9
max 30.0 30.0 1000000 10000000 0.1 0.1 
# define upper and lower bounds of x0 and theta for prior_data.
max_x_par=c(c(30,30,c(1000000,10000000,.1,0.1)))
min_x_par=c(10,0,100000,1000000,1e-9,1e-9)

4.2 Optimization to get set of parameters

Within the range, n-step ahead prediction (Eq. 14 in the paper) is done to get several local optimal parameter sets.

# Not run this code because it takes time. 
# We will reuse the store results.

# check optimizer runs. Single run.
# ss=DEoptim::DEoptim(fn=get_rmse, lower=rep((-0.5+1e-9),length(L_val)),
#                      upper=rep(0.5,length(L_val)),control=list(NP=500, itermax=10000,trace=TRUE,reltol=1e-18,steptol=40,parallelType=1,
#                                                   parVar=list("discretize_system","nstep_cpp","expm")),
#                      dt=dt,u_mat=train_data$u_mat,y_mat=train_data$y_mat,dims=dims,output="cost",normalize=TRUE,L_val=L_val)

sol_list=list() #store optimization results in a list.
iii=1  # to create a random seed
# sometimes optimization fails because cost is NaN. Thus, have this trycatch loop with while so that the optimizer learns again until we get satisfactory number of outcomes.

while(length(sol_list)<=100){
  set.seed(iii+30000)
  # sometimes optimizer fails since it gives NaN during the nstep-ahead prediction. Thus, we use try-catch.
  try_result=tryCatch(DEoptim::DEoptim(fn=get_rmse, lower=min_x_par/max_x_par,
                     upper=max_x_par/max_x_par,
                     control=list(NP=100, itermax=10000,trace=TRUE,reltol=1e-16,steptol=15,parallelType=1,
                                                  parVar=list("discretize_system","nstep","expm")),
                     dt=dt,u_mat=prior_data$u_mat,y_mat=prior_data$y_mat,dims=dims,output="cost",normalize=TRUE,L_val=max_x_par),
                       error=function(e){"error"} )
  
  if((try_result=="error")){
    
  }else{
    sol_list[[(length(sol_list)+1)]]=try_result
  }
  iii=iii+1
}
# store solution.
write_rds(sol_list,paste0("data/sol_list.rds"))

4.3 Visualize obtained set of parameters

To run Liu-West filter on the filter_data, we need both prior for initial states and parameters. In 3.3.2 section of the research paper, initial states and parameters of prior_data are obtained from optimizer. The initial states of filter_data is actually the last states of prior_data. Therefore, we do n-step ahead prediction on prior_data with set of parameters from the obtimizer to have initial states of filter_data.

# load optimization results
sol_list=readr::read_rds(paste0("../data/sol_list.rds"))

# sol_list contains intial states and parameters of prior_data.
# we will have initial state and parameters for filter_data in sol_list_temp
sol_list_temp=list()

for (i in 1:length(sol_list)){
  sol_list_temp[[i]]=(sol_list[[i]]$optim$bestmem)*max_x_par #extract optimizer solution and unnormalize it.
  # store n-step ahead prediction result
  # 1:index_split: prior data, index_split+1:index_split*2: filter_data
  tempp=get_rmse(par=sol_list_temp[[i]],dt=dt,
                 u_mat=all_data$u_mat[,1:(index_split+1)],
                 y_mat=all_data$y_mat[1:(index_split+1)],
                 dims=dims,output="data",normalize=FALSE)
  # put indepx_split+1 states into initial states into sol_list. 
  sol_list_temp[[i]][1:2]=tempp$x_pred[,index_split+1]
}
# sol_list_temp to dataframe
all_df=do.call(rbind,sol_list_temp)%>%as_tibble()

state_names<-c("Ti0","Te0")
par_name<-c("Ci","Ce","Rie","Rea")
# split data_frame into state/parameter
par_df=all_df[,-c(1:2)]%>%set_names(par_name)
state_df=all_df[,1:2]%>%set_names(state_names)

After getting initial states for filter_data, visualize the results. As you can see the red dot (true parameter) are in a range of various local optimal solutions.

One distrinct advantage of this approach is it can capture the correlation betwen parameters in a prior distribution. For example, \(R_{\text{ie}}\) is highly correlated with \(R_{\text{ea}}\).

4.4 Kernel density approximation of set of pameters

From the optimization results, we will generate initial prior distribution. To do this, we will approximate the optimization results by using multivariate Kernel density approximation.

## Not run this code because we will use generated particles.
# multivariate kernel density approximation
kde_par <- kdevine::kdevine(as.matrix(par_df)) #kernel model K(x)
kde_state <- kdevine::kdevine(as.matrix(state_df)) #kernel model K(theta)

# Generate 
set.seed=104
NP=10000 #number of particles
# Generate particles from the approximated kernel density
prior_par=Rfast::transpose(abs(kdevine::rkdevine(NP, kde_par)))
prior_state=Rfast::transpose(kdevine::rkdevine(NP, kde_state))
# store generate particles 
write_rds(kde_par,paste0("../data/kde_par.rds"))
write_rds(kde_state,paste0("../data/kde_state.rds"))
write_rds(prior_par,paste0("../data/prior_par.rds"))
write_rds(prior_state,paste0("../data/prior_state.rds"))

Visualization of generated particles. The generated prior particles (gray dots) cover wider ranges than optimization results while it captures the correlation between parameters.

5. Liu-West filter

Since we have priors, we can run Liu-West filter. For complete algorithm of the filter, refer Table 2 in the paper. The filter code is given lw_ss_filter in utility.R file.

5.1 Parameter normalization and generate noise parameters.

# L_values of parameter to normalize parameters within a range of (0,1].
# After dividing by L_values, we subtract 0.5 to make it a range of (-0.5,0.5]
# This is obtained maximum values of generated iniital priors
L_val_par=apply(prior_par,MARGIN=c(1),max) 
prior_par_n=(prior_par)/L_val_par-0.5 #normalized parameter 
num_state=dim(prior_state)[1] # number of state

# define L_values of noise parameter (sigma_x_i,sigma_x_e)
L_val_sd=rep(0.1/sqrt(dt),num_state) 
L_val=c(L_val_par,L_val_sd) #L_val for parameter +noise parameter

# normalized priors for noise parameter by using Latin-Hyper cube sampling
set.seed(1234)
prior_sd_n=(Rfast::transpose(lhs::randomLHS(NP,length(L_val)-(dim(prior_par_n)[1] ))))-0.5

# all normalized parameters.
prior_all_n=rbind(prior_par_n,prior_sd_n)
num_par=dim(prior_all_n)[1]

# initial x0
x=prior_state #x0
pii=w=rep(1/NP,NP) #initilal weight

# inputs for the model
inputs=list()
inputs$L_val=L_val #L_val for normalization/unnormalization.
inputs$dt=dt #time interval

inputs$y_mat=filter_data$y_mat # y_data
inputs$u_mat=filter_data$u_mat # input u_data

# parameter names
par_names=c("Ci","Ce","Rie","Rea","sdte","sdti")
inputs$par_names=par_names

# initial parameter prior particles
theta=prior_all_n
delta=0.9 # Filter tuning parameter 
seed_num=13244 #seed number

5.2 Liu-West filter calculation

Run the filter and store the results.

# Not execute code. Takes time.
# filter run
res=lw_ss_filter(inputs=inputs,NP=NP,pii=pii,x=x,theta=theta,delta=delta,seed=seed_num,training=TRUE)

# split data to reduce file size
for(i in 1:dim(res$r_thetas)[1]){
  write_rds(res$r_thetas[i,,],paste0("../data/res_thetas_",i,".rds"))
}
for(i in 1:dim(res$xs)[1]){
  write_rds(res$xs[i,,],paste0("../data/res_xs_",i,".rds"))
}
res$thetas<-NULL
res$r_thetas<-NULL
res$xs<-NULL
write_rds(res,"../data/res.rds")

5.3 Visualize filter results

Visualize the filter results and compare with true paraemter values. When the filter moves, the posterior distribution updates and is narrower with data. As can be seen the graph, the true parameter value (red horizontal line) is within the posterior distribution, and we can confirm the proposed method can be applicable to the simple gray-box model of a building.

6. Conclusion

In this document, we create a synthetic dataset from a simple building gray-box model and apply Liu-West filter to find posterior distribution of parameters. By making a comparison with the true parameter, we demonstrate the filter can be used in this applicaiton.