- Email: [email protected]
Parameter identifiability in Bayesian inference for building energy models
Energy & Buildings 198 (2019) 318–328
Contents lists available at ScienceDirect
Energy & Buildings journal homepage: www.elsevier.com/locate/enbuild
Parameter identifiability in Bayesian inference for building energy models Dong Hyuk Yi a, Deuk Woo Kim b, Cheol Soo Park c,∗ a
Department of Architecture and Architectural Engineering, College of Engineering, Seoul National University, 1, Gwanak-ro, Gwanak-gu, Seoul, 08826, South Korea b Department of Living and Built Environment Research, Korea Institute of Civil Engineering and Building Technology, 283, Goyang-daero, Ilsanseo-gu, Goyang, Gyeonggi, 10223, South Korea c Department of Architecture and Architectural Engineering, Institute of Construction and Environmental Engineering, Institute of Engineering Research, College of Engineering, Seoul National University, 1, Gwanak-ro, Gwanak-gu, Seoul, 08826, South Korea
a r t i c l e
i n f o
Article history: Received 15 January 2019 Revised 13 May 2019 Accepted 5 June 2019 Available online 6 June 2019 Keywords: Bayesian inference Parameter identifiability Likelihood confidence interval Likelihood confidence region Biplot
a b s t r a c t Parameter identifiability is the concept of whether uncertain parameters can be correctly estimated from the observed data. The main cause of parameter unidentifiability in Bayesian inference is known as ‘overparameterization’. In this study, the likelihood confidence interval (CI) and the likelihood confidence region (CR) were introduced for quantifying the parameter identifiability. The likelihood CI and CR can be regarded as the parameter range (one-dimensional) and parameter space (two-dimensional or higher) that can identify parameter values, respectively. For this purpose, an EnergyPlus reference office building provided by the US DOE was used in this study. Four estimation parameters in the EnergyPlus model were analyzed using the likelihood CI and CR. It was found that the closer the likelihood CI of a parameter is to the prior’s parameter range, the more unidentifiable the parameter. In addition, a biplot analysis was conducted to examine a correlation between two parameters. The more correlated a parameter is with others, the more unidentifiable the parameter. It is suggested that the visual assessment of likelihood CIs and CRs can help in investigating whether Bayesian inference results can be accurately obtained. © 2019 Elsevier B.V. All rights reserved.
1. Introduction Even though building simulation disciplines have been developed over the past several decades, the performance gap between simulation prediction and measured data is still significant [1–3]. One of the reasons for the gap is the “stochastic nature” in the dynamic behavior of building systems. Therefore, it is important to reflect such stochastic characteristics of simulation inputs in a simulation model [4,5]. Hence, Bayesian inference has been highlighted owing to its capability to estimate probability distribution of uncertain inputs through the combination of domain knowledge (i.e., prior information) and observed data (i.e., likelihood) [6]. Bayesian inference has been adopted for energy retrofit analysis [7–13] and urban-scaled building energy modeling [8,9,14–17]. However, the following issues remain unsolved: (1) observed data, (2) prior choice, (3) simulation running time and (4) parameter identifiability.
∗
Corresponding author. E-mail addresses: [email protected] (D.H. Yi), [email protected] (D.W. Kim), [email protected] (C.S. Park). https://doi.org/10.1016/j.enbuild.2019.06.012 0378-7788/© 2019 Elsevier B.V. All rights reserved.
• Observed data: This is to check whether the observed data are suitable for parameter estimation of a simulation model. It has been reported that Bayesian inference can be influenced by insufficient data [18], low precision of observed data [19], low temporal resolution of observed data [20], the lack of correlation between measured data and simulation output [21], etc. • Prior choice: The posterior inference, or the outcome of Bayesian inference, is influenced by the prior set that is usually guessed at by an analyst. In other words, the posterior can be biased by the prior distribution (e.g., Gaussian prior, uniform prior, triangular prior). In addition, the result of Bayesian inference is also affected by the prior informativeness (e.g., the prior with large variance vs. the prior with small variance) [22,23]. • Simulation running time: Since Bayesian inference is based on Monte Carlo simulation (e.g., random walk Markov chain Monte Carlo [MCMC] [24–26], importance sampling [27–29], Hamiltonian Monte Carlo [30,31], and No-U-Turn sampler [32,33]), the simulation running time is demanding. For this reason, many Bayesian inference studies have used meta models as a surrogate model (e.g., linear regression model [13,15,17,34,35], artificial neural network model [36], random forest model
D.H. Yi, D.W. Kim and C.S. Park / Energy & Buildings 198 (2019) 318–328
319
Fig. 1. Scope of this study (shown in the dotted box) vs. those of previous studies.
[36], Gaussian process model [10,12,37,38]) or normative quasisteady-state model [7–9,11]. In addition, because the simulation running time is also influenced by the amount of observed data, a recent study [18] has introduced a Kullback-Leibler divergence method [39], which can select the subset from the entire observed dataset with low information loss. • Parameter identifiability: Uncertain parameters could be wrongly estimated when the parameters to be estimated are overparameterized (i.e., too many parameters to be estimated) [20,23,31]. In such cases, some parameters of the model become unidentifiable. In addition, when parameters are correlated with each other, the degree of the identifiability is also reduced [40]. The previous studies [20,23,31] have not described (1) how the overparameterization can change the parameter’s likelihood or (2) a causality between the parameter’s likelihood and Bayesian inference results. Thus, this study can be regarded as a follow-up study of [20,23,31] and aims to investigate the parameter identifiability problem. As shown in Fig. 1, several studies have addressed simulation models [7–13,15,17,34–38], observed data [18–21], prior [22,23], MCMC [23,31] and posterior [20,23,31] in Bayesian inference. With regard to the parameter identifiability, an ‘indirect’ method has been applied, e.g., monitoring MCMC convergence [23,31] or comparing prior to posterior [20,23,31]. Rather than using the ‘indirect’ method, this study suggests a ‘direct’ method using the likelihood confidence interval (hereinafter referred to as likelihood CI) and the likelihood confidence region (hereinafter referred to as likelihood CR) as means to measure whether parameters can be estimated from the observed data [41,42]. Previous studies [22,31,40] found that the number of parameters (i.e., overparameterization) can cause an unidentifiability problem. In this regard, this paper shows how the likelihood CI and CR can vary according to the number of parameters in the model as well as how the likelihood CI and CR can be used as a measure for parameter identifiability. In this study, the authors selected a reference office building presented by the US DOE [43] and generated a meta-model using an artificial neural network (ANN). Then, we analyzed the parameter identifiability based on the likelihood CIs and CRs. The parameters’ correlation, one of the reasons for parameter identifiability, was discussed using a biplot analysis, which is a visual representation of the correlation between parameters. Finally, the impact of parameter identifiability on the Bayesian inference results was addressed.
2. Methodology 2.1. Bayesian inference and MCMC Bayesian inference is based on a conditional probability theorem, as shown in Eq. (1). The posterior (P(θ |D)) is the conditional probability of uncertain parameter θ given data D. P(θ |D) can be derived from the product of the prior (P(θ )) and the likelihood (l(D|θ )) [44]. The prior, a set of uncertain parameters, is usually subjectively assumed by an analyst [45]. The level of the analyst’s certainty for the parameter value can vary. Depending on the level of the certainty, the prior can be classified as informative prior (specifying the expected parameter value in the prior) or uninformative prior (expressing the ambiguity about the parameter value in the prior) [46]. Typical examples are Gaussian prior (informative prior) and uniform prior (uninformative prior). When the uninformative prior is applied, an equal weight is applied throughout the parameter range, so that the posterior is proportional to the likelihood [47,48]. More details on the calculation of the posterior from the prior can be found in [44,49].
P (θ |D ) ≈ P (θ )l (D|θ ),
(1)
where
θ : Uncertain parameter D: Observed data P(θ |D): Posterior P(θ ): Prior l(D|θ ): Likelihood The likelihood is expressed as a function of the observed data under each parameter value as shown in Eq. (1). Unlike the prior, the analyst’s assumption is not involved in defining the likelihood. This can be regarded as objective information in the Bayesian framework as opposed to the prior [45]. The likelihood is equivalent to the probabilistic expression of the error term (e in Eq. (2)), including both the observation error and systematic gap between the observed data and the model [44]. Assuming that the error follows the Gaussian form, the likelihood can be defined as shown in Eq. (3).
D = f ( θ ) + e,
(2)
1 T l (D|θ ) ∝ exp − (D − f (θ ) ) C −1 (D − f (θ ) ) , 2 where
(3)
320
D.H. Yi, D.W. Kim and C.S. Park / Energy & Buildings 198 (2019) 318–328
f(θ ): Model output (i.e., physical model) e: Error between the data and the model C: Covariance matrix of the likelihood For sampling, an adaptive metropolis algorithm [50–52], an MCMC algorithm, was used in this study (Eq. (4)). This algorithm is a variation of the random walk metropolis algorithm [24–26] and improves convergence by updating the covariance matrix ( n in Eq. (4)) of the subsamples generated at regular update intervals in the whole sampling process and updating it (λ n in Eq. (4)). In this study, 50,0 0 0 MCMC samples (burn-in samples: 25,0 0 0, posterior samples: 25,0 0 0, regular update intervals: 50 0 0) were generated for each Bayesian inference case, and the posterior was derived by the Gaussian kernel density estimation with the MCMC samples. Burn-in refers to discarding an initial portion of a Markov chain sample so that the effect of initial values on the posterior inference is minimized [53]. For this MCMC calculation, the Python MCMC module (PyMC 2 [54]) was used.
Yn+1 ∼ N (Xn ,
λn ),
(4)
where Yn+1 : Proposal distribution used at next period (n + 1) N: Multivariate Gaussian distribution Xn : Posterior samples λ: Scaling factor (= 2.38 / parameter dimensions) n : Covariance matrix of subsamples n (subscript): Current posterior sampling period
of the likelihood CIs and CR for two parameters (X1, X2). The likelihood CR projected on the X1 and X2 axes represents the likelihood CIs for X1 and X2, respectively. The exact method to calculate Eqs. (7) and ((8) can be found in [42].
θˆ = argmax(l (D|θ ) ),
(6)
χ 2 (θ ) = (D − f (θ ))T C −1 (D − f (θ )) ≈ −2 ln (l (D|θ ) ),
(7)
⎧ ⎫ ⎨ ⎬ l D | θ ( ) θCR −2ln = χ 2 (θ ) − χ 2 θˆ < α , ⎩ ⎭ l D|θˆ
2.2. Likelihood-based parameter identifiability analysis: CI, CR and biplot The parameter identifiability analysis investigates whether the model parameters can be estimated based on the observed data [40]. In other words, this analyzes whether the parameter values can be determined within a finite parameter space. To explain the concept of the parameter identifiability, l(D|θ ) in Eq. (1) is taken again. Let a model be completely specified by the parameter set θ , and let x be a vector of observed random variables, where x ∈ ( means all possible sets of the observations). The parameters of the model are unidentifiable if there exists some θ 1 = θ 2 yielding the same distributions (Eq. (5)) for all x ∈ [40,55–57]. In other words, a combination of two different parameters can produce the same likelihood. These unidentifiable parameters would make it difficult to estimate the other parameters’ values [40].
l (x|θ1 ) = l (x|θ2 ),
Fig. 2. Likelihood CI and CR for the parameter pair (X1, X2).
(5)
The parameter identifiability is influenced by the following factors: the number of parameters, parameter redundancy, the range (min, max) of parameters, correlation between parameters, and the quantity and quality of the observed data. Please note that these factors can cause the parameter identifiability problem individually or interconnectedly [40,41,57–63]. In this study, the likelihood CI and the likelihood CR were introduced to analyze the parameter identifiability. The likelihood CI and CR are the parameter range (one-dimensional) and the parameter space (two-dimensional or higher) that can identify the parameter likelihood values, respectively [41,42]. The parameter value (θˆ ) that maximizes the likelihood can be defined as shown in Eq. (6). According to the asymptotic sampling distribution theory [41,42], the normal distributed likelihood (l(D|θ )) can be converted to a chi-squared distribution (χ 2 (θ )) Eq. (7)). Then, a parameter subspace that satisfies the condition that the difference (χ 2 (θ ) − χ 2 (θˆ )) is less than the likelihood threshold (α ) is the likelihood CR (θ CR in Eq. (8)). Fig. 2 shows a visual representation
(8)
where
θˆ : Parameter value that maximizes the likelihood θ CR : Parameter value within the likelihood CR χ 2 (θ ): Chi-squared distribution α : Likelihood threshold (inverse cumulative distribution function of the Chi-squared distribution for α confidence level) In addition, biplot [64–66] analysis was used to examine the correlation between parameters. The biplot can visualize a dataset with more than two parameters in a low-dimensional graph (in general, two dimensions are chosen) [66]. The dimensional reduction of the dataset can be performed through a multivariate analysis based on a singular value decomposition (e.g., principal component analysis [67,68], correspondence analysis [69], etc.). In the biplot analysis, each parameter can be represented as a vector on a 2-D graph. The parameter vectors are obtained with respect to the first two principal components (hereinafter referred to as PC1 and PC2, respectively, selected in the order of their eigenvalues), each weighted by the standard deviation of that component or the square root of the corresponding eigenvalue [66]. Fig. 3 shows an example from a three-dimensional likelihood CR to a two-dimensional biplot. For the biplot analysis, Pandas [70], Matplotlib [71], and Scikit-Learn [72] modules were used. A normalization process (e.g., the min-max method and Z-score method) is required because the units/values of each parameter can be different from one another. For the dimensional reduction, it is recommended that the sum of the explanatory power (i.e., the ratio of the eigenvalue of each principal component to the eigenvalue sum of two principal components) for two principal components must be greater than 70% [73] (in Fig. 3, the sum is 58.8% + 23.6% = 82.4%). Please note that the correlation between two parameters is equal to the cosine of the angle between the vectors (Fig. 3(c)),
D.H. Yi, D.W. Kim and C.S. Park / Energy & Buildings 198 (2019) 318–328
321
Fig. 3. Biplot Example.
ing power density, and infiltration rate. The minimum, median and maximum values of the parameters were determined from [75,76]. The original model values are defined as ‘true values’. Please note that the authors did not use an on-site measured data because the on-site data are usually influenced by unknown disturbances/noises. This study was intended to suggest the need for the parameter identifiability analysis as a preliminary process. That’s why the authors used a simple ‘toy experiment’ clean model (EnergyPlus) as shown in Fig. 4. In addition, what matters for parameter identifiability is not “the number of output variables”, but “the informativeness of the measured data”, which can be checked by the likelihood CI, CRs and correlations between parameters. Since the MCMC in Bayesian inference demands significant simulation running time, the artificial neural network (ANN) model was employed as a surrogate model to the EnergyPlus model. The process to generate the surrogate model was as follows:
Fig. 4. Target building.
as shown in Eq. (9). The absolute value ranges of the correlation coefficient (i.e., cosine value between two vectors of the biplot) are classified as 0 to 0.3 (very weak), 0.3 to 0.5 (weak), 0.5 to 0.7 (moderate), and 0.7 to 1.0 (strong) according to [74]. For example, the cosine value between X1 and X2 is 0.55 (Fig. 3(c)), indicating a moderate positive correlation between X1 and X2. The cosine value between X1 and X3 is −0.28, indicating a weak negative correlation. The cosine value between X2 and X3 is −0.96, meaning a strong negative correlation. If the directions between two vectors are opposite to each other, it means a strong negative correlation.
( v1 · v2 ) Cos(v1, v2) = , (|v1||v2| )
(9)
3. Target building The EnergyPlus model for a mid-size reference office building developed by the US DOE [43] was chosen as the target building (Fig. 4). The building is located in Atlanta, Georgia, USA, and has three stories above ground, a 4982 m2 floor area, and a 33% window-wall ratio. Each floor consists of four perimeter zones and one interior zone. The target building’s plant is a furnace for heating and packaged air-conditioning units for cooling. The air handing unit of the building is a VAV system with reheat coils. The building uses gas for heating and electric energy for cooling, lighting, fan, pump, and plug loads. In this study, the authors selected both annual gas energy use (kWh) and annual electric energy use (kWh) as the model outputs, or the observed data. As shown in Table 1, the following four parameters were chosen as Bayesian inference objects: fenestration U value, fenestration SHGC, light-
• Step 1: perform Latin hypercube sampling [77] based on Table 1. The number of sampled EnergyPlus simulation models is 10 0 0. • Step 2: Based on the generated EnergyPlus simulation models (Step 1), conduct a series of EnergyPlus simulation runs. • Step 3: Collect simulation results such as annual electric energy and gas uses. • Step 4: Construct ANN models using inputs and outputs of the EnergyPlus models; 750 and 250 EnergyPlus simulation runs were used for training and test, respectively. For steps 1 to 4, the following Python modules were used: pyDOE [78] for Latin hypercube sampling, Eppy [79] for generating EnergyPlus models, and scikit-learn [72] for training the ANN model (optimization solver: Adam, activation function: ReLU function, the number of hidden layers: 1, the number of neurons: 50). The generated ANN models are sufficiently accurate when compared to EnergyPlus model outputs (MBE: −0.01%, CVRMSE: 0.13%, R2 : 0.97 for gas; MBE: 0.02%, and CVRMSE: 0.35%, R2 : 0.98 for electric energy) (Fig. 5). In this study, the ANN model was saved as a class variable for Python. The authors used the Pickle module (object serialization module for Python) to export the ANN class variable to a binary file that can be read in Python for the purpose of using the ANN model in MCMC iteration. 4. Results Three comparison cases were made, as shown in Table 2. The baseline case estimates all four parameters, presented in Table 1. 30 0,0 0 0 Monte Carlo simulation runs using Latin hypercube sampling were completed for each case in Table 2. The computational
322
D.H. Yi, D.W. Kim and C.S. Park / Energy & Buildings 198 (2019) 318–328 Table 1 List of estimation parameters. Parameters
True value
Min
Median
Max
References
Fenestration U value (W/m2 K) Fenestration SHGC (0.0–1.0) Lighting power density (W/m2 ) Infiltration rate (air changes per hour, 1/h)
4.09 0.25 16.89 0.35
1.03 0.13 7.0 0.15
3.705 0.495 14.0 0.45
6.38 0.86 21.0 0.75
[75] [75] [75] [76]
Fig. 5. Comparison between EnergyPlus and ANN models (250 test data used).
Fig. 6. Biplot analysis.
Table 2 Comparison cases. Number of parameters Baseline Case I Case II
4 3 2
costs of using the method proposed in this study were as follows: 308.3, 304.7, and 316.9 s for Baseline, Case I, and Case II, respectively (CPU: AMD Ryzen 5 1600X 3.60 GHz, RAM: 8 GB). When one parameter is removed, its true value in Table 1 was used for Cases I and II. In this study, 95% likelihood CIs and 95% likelihood CRs were used. 4.1. Likelihood CI and biplot analysis Table 3 shows the ratios of likelihood CIs for four parameters to their parameter ranges, Table 4 shows a visual comparison of the likelihood CIs for four parameters to their ranges, Table 5 shows
the degrees of correlation between two parameters, and Fig. 6 shows the biplot analyses for the three cases. The results are summarized as follows: • Baseline case: Except lighting power density, the ratios of the likelihood CIs of the other three parameters to their parameter ranges are greater than 0.90 (Table 3). In other words, the three parameters are difficult to estimate because their likelihood CIs are close to their parameter ranges. Particularly, the likelihood CI of the infiltration rate is the same as its parameter range, and its likelihood shape is flat, indicating that its likelihood over its parameter range is equally distributed (Table 4). As shown in Table 5 and Fig. 6, strong correlations exist between ‘Fenestration U value and Fenestration SHGC (0.83)’ and ‘Fenestration SHGC and Lighting power density (−0.93)’. Because Fenestration SHGC is a shared parameter included in the aforementioned parameter pairs, it can be inferred that Fenestration SHGC may cause an identifiability problem due to its correlation with other parameters [40]. In contrast, Fenestration U value, which is difficult to identify (B/A = 0.93, Table 3), has a moderate or weak correlation (−0.56 with lighting power
D.H. Yi, D.W. Kim and C.S. Park / Energy & Buildings 198 (2019) 318–328
323
Table 3 Ratios of parameters’ likelihood CIs to parameter ranges. Parameter
Parameter range (A) (Table 1)
Cases
Likelihood CI (B) (refer to Table 4)
B/A
Fenestration U value (W/m2 K)
6.38–1.03 = 5.35
Baseline Case I Case II
6.38–1.39 = 4.99 6.38–1.83 = 4.55 6.38–2.23 = 4.15
0.93 0.85 0.78
Fenestration SHGC (0–1)
0.86–0.13 = 0.73
Baseline Case I Case II
0.84–0.13 = 0.71 0.83–0.13 = 0.70 n/a
0.97 0.96 n/a
Lighting power density (W/m2 )
21.0–7.0 = 14.00
Baseline Case I Case II
20.88–11.58 = 9.30 20.64–12.26 = 8.38 20.10–13.49 = 6.61
0.66 0.60 0.47
Infiltration rate (ACH)
0.75–0.15 = 0.60
Baseline Case I Case II
0.75–0.15 = 0.60 n/a n/a
1.00 n/a n/a
The bold numbers indicate that the corresponding parameter is excluded in the next case study. Table 4 Comparison of likelihood CIs to parameter ranges.
1.0
0.5
3.70
6.38
1.0
0.5
0.0 1.03
Likelihood CI (0.13, 0.84)
1.0
0.5
0.0 0.13
0.50
0.86
0.5
0.0 21.0
Log likelihood (Normalized)
Log likelihood (Normalized)
Likelihood CI (11.58, 20.88)
Lighting power density Log likelihood (Normalized)
0.5
0.0 1.03
3.70
6.38
Fenestration U value
Likelihood CI (0.13, 0.83)
excluded in Case II
0.5
0.0 0.13
0.50
0.86
Fenestration SHGC
1.0
14.0
1.0
1.0
Fenestration SHGC
7.0
6.38
Likelihood CI (2.23, 6.38)
Fenestration U value Log likelihood (Normalized)
Log likelihood (Normalized)
Fenestration U value
3.70
Likelihood CI (12.26, 20.64)
1.0
0.5
0.0
14.0
7.0
21.0
Lighting power density
Log likelihood (Normalized)
0.0 1.03
Case II
Likelihood CI (1.83, 6.38)
Log likelihood (Normalized)
Case I Log likelihood (Normalized)
Log likelihood (Normalized)
Baseline Likelihood CI (1.39, 6.38)
Likelihood CI (13.49, 20.10)
1.0
0.5
0.0 7.0
14.0
21.0
Lighting power density
Likelihood CI (0.15, 0.75) 1.0
excluded in Case I
0.5
0.0 0.15
0.45
Infiltration rate
excluded in Case II
0.75
density, −0.34 with infiltration rate) except for the Fenestration SHGC. The infiltration rate has a moderate, weak, or very weak correlation (0.24, −0.34, −0.59) with the other three parameters. Hence, it can be said that Fenestration U value and Infiltration rate cause the identifiability problem independently when compared with the Fenestration SHGC. Based on the above
analysis, it was decided that infiltration rate should be excluded from the list of estimation parameters for Case I. • Case I: As shown in Table 3, the Fenestration U value (0.93 → 0.85) and Fenestration SHGC value (0.97 → 0.96) are found to have high ratios of the likelihood CIs to their parameter ranges, indicating that the two parameters are difficult
324
D.H. Yi, D.W. Kim and C.S. Park / Energy & Buildings 198 (2019) 318–328 Table 5 The degrees of correlation between two parameters (cosine values between two parameters). Two vectors
Baseline
Case I
Case II
Fenestration U value vs. Fenestration SHGC Fenestration U value vs. Lighting power density Fenestration U value vs. Infiltration rate Fenestration SHGC vs. Lighting power density Fenestration SHGC vs. Infiltration rate Lighting power density vs. Infiltration rate
0.83 −0.56 −0.34 −0.93 0.24 −0.59
0.55 −0.28 excluded −0.96 excluded excluded
excluded −0.21 excluded excluded excluded excluded
Table 6 Six priors (P1–P6). mean Parameters Fenestration U value (W/m2 K) Fenestration SHGC (0–1) Lighting power density (W/m2 ) Infiltration rate (ACH) ∗
P1
P2
P3
P4
3.705 0.495 14.0 0.45
5.80 0.69 10.77 0.72
4.12 0.82 15.80 0.33
2.68 0.36 9.31 0.51
P5
P6
4.36 0.53 14.15 0.59
1.61 0.21 20.08 0.16
Std. Dev.
min∗
max∗
1.37 0.19 3.57 0.15
1.03 0.13 7.00 0.15
6.38 0.86 21.0 0.75
The minimum and maximum values are the same as those in Table 1.
Table 7 Posterior modes. True value
Parameter range (A) (from Table 1)
Fenestration U value (W/m2 K)
4.09
6.38–1.03 = 5.35
Fenestration SHGC (0–1)
0.25
Lighting power density (W/m2 ) Infiltration Rate (ACH)
Parameter
Posterior modes Case
(Max – min) of six posterior modes (C)
C/A
B/A (from Table 3)
3.51 3.38 3.38
2.68 2.78 1.71
0.50 0.52 0.32
0.93 0.85 0.78
0.44 0.48 n/a
0.16 0.15 n/a
0.57 0.47 n/a
0.78 0.64 n/a
0.97 0.96 Excluded
14.80 15.22 15.99
14.94 15.06 16.53
17.56 17.48 17.57
3.80 2.91 2.00
0.27 0.21 0.14
0.66 0.60 0.47
0.57 n/a n/a
0.63 n/a n/a
0.24 n/a n/a
0.48 n/a n/a
0.80 n/a n/a
1.00 Excluded Excluded
P1
P2
P3
P4
P5
P6
Baseline Case I Case II
4.26 3.98 3.98
6.19 6.16 5.09
4.97 4.77 4.15
3.65 3.89 3.60
4.50 5.02 4.10
0.86–0.13 = 0.73
Baseline Case I Case II
0.44 0.37 n/a
0.64 0.62 n/a
0.61 0.61 n/a /
0.39 0.36 n/a
16.89
21.0–7.0 = 14.00
Baseline Case I Case II
14.96 15.35 16.37
13.76 14.57 15.57
14.99 15.23 16.60
0.35
0.75–0.15 = 0.60
Baseline Case I Case II
0.48 n/a n/a
0.72 n/a n/a
0.38 n/a n/a
to estimate. The likelihood shapes of the two parameters are similar to those of Case I (Table 4). According to Table 5 and Fig. 6, there is a strong correlation (−0.96) between ‘Fenestration SHGC and Lighting power density’. It can be inferred that the identifiability problem exists due to the correlation between Fenestration SHGC and Lighting power density. Therefore, it was determined that Fenestration SHGC should be excluded for Case II. • Case II: After two parameters (Infiltration rate and Fenestration SHGC) were excluded from the list of estimation parameters, Case II was conducted. The ratios of the likelihood CIs of Fenestration U value and lighting power density to their parameter ranges are reduced (0.93 → 0.78, 0.66 → 0.47). In addition, distinct forms of peaks are observed in the likelihood shapes of the two parameters (Fenestration U value and lighting power density), compared to those of the baseline case and Case I (Table 4). The cosine value between the two parameters (Fenestration U value and Lighting power density) is −0.21, implying that the correlation between the two parameters is very weak [74] (Table 5, Fig. 6). It is expected that the estimation for the two parameters would be correct partly because the parameter subspace is reduced from the entire parameter space. This result will be substantiated in Section 4.2. 4.2. Likelihood CR analysis To investigate how posterior can be biased depending on prior, the following six priors (P1-P6) were made as shown in Table 6.
The six priors were assumed to follow a Gaussian distribution. P1 represents the median values in Table 1, while P2–P6 were sampled using Latin hypercube sampling [77]. The standard deviations of P1-P6 for each parameter were assumed to be equal to (mean of P1 – min of P1) /1.96. Table 7 tabulates the six posterior modes, the difference between their maximum and minimum modes (denoted as ‘C’), and the ratio of C to A (parameter range). As explained in Section 4.1, the ratio (B/A) can be regarded as a measure for the parameter ‘unidentifiability’. The greater the ratio (B/A) is, the more unidentifiable the parameter. Similar findings can be drawn by comparing the ratio (C/A) to (B/A). It can be concluded that as the ratio (C/A) increases, it is more likely that the Bayesian results are farther away from the true value. Fig. 7 shows the true values (a blue star), their posterior modes (green circles), and their prior means (white squares) on the likelihood CRs (red dots) in a 2-D space. Because the correlations between ‘Fenestration U value and Fenestration SHGC’ are positively linear (the cosine values are 0.83 in Baseline case and 0.55 in Case I, Table 5), the likelihood CRs and posterior modes of the two are also positively related, as shown in Fig. 7(a). In contrast, ‘Fenestration SHGC and Lighting power density’ are negatively correlated (the cosine values between the two are −0.93 in Baseline case and −0.96 in Case I, Table 5), and the likelihood CRs and posterior modes of the two are found to be negatively linear (Fig. 7(b)). In other words, the posterior modes of two parameters tend to follow the correlation between the two.
D.H. Yi, D.W. Kim and C.S. Park / Energy & Buildings 198 (2019) 318–328
Fig. 7. Likelihood CRs for estimation parameters. (For interpretation of the references to color in the text, the reader is referred to the web version of this article.)
325
326
D.H. Yi, D.W. Kim and C.S. Park / Energy & Buildings 198 (2019) 318–328
Fig. 8. The parameter unidentifiability increases as the likelihood CI increases. (For interpretation of the references to color in the text, the reader is referred to the web version of this article.)
It is noteworthy that the greater the ratios of the parameters’ likelihood CIs to the parameter ranges (Table 3), the wider the parameters’ likelihood CRs and posterior modes were estimated. For example, the likelihood CR and posterior modes for the Infiltration rate, of which the ratio is 1.00 (Table 3), are widely spread over its parameter range (Fig. 7(c), (e) and (f)). Similarly, the likelihood CR and posterior modes of Fenestration SHGC, of which the ratios are 0.97 in the Baseline case and 0.96 in Case I (Table 3), are widely placed across its parameter range (Fig. 7(a), (d) and (e)). From the Baseline case via Case I to Case II, the likelihood CRs of the parameters are gradually reduced owing to excluded parameters as well as removed parameter correlations. For this reason, the posterior modes of lighting power density and Fenestration U value tend to converge to their true values (Case II in Fig. 7(b), Table 7). In particular, owing to the exclusion of Fenestration SHGC (the parameter that caused the strong correlation with other parameters, Table 5), a significant change in the shape of the likelihood CR and the posterior modes occurred from Case I to Case II (Fig. 7(b)). The black solid line in Fig. 8 shows the relationship between the two ratios (B/A and C/A, Table 7). As B/A increases, the difference between the minimum and maximum posterior modes increases (corresponding to the vertical height of the dotted red circles in Fig. 8), indicating that the Bayesian inference results are widespread. In other words, the parameters might be incorrectly estimated, not sufficiently close to the true value (Table 7, Fig. 8). As presented in Table 7 and Fig. 8, it can be inferred that the greater the ratio of the parameter’s likelihood CI to the parameter range is, the greater the parameter unidentifiability. The degree of the uncertainty in posterior modes of the parameters increases in proportion to the ratio of the parameter’s likelihood CI to the parameter range. Conclusively, this study shows that the likelihood CIs and CRs can be used to investigate the parameter identifiability in the Bayesian inference problems, as illustrated in Fig. 8.
5. Conclusion In this paper, the authors suggested a methodology to investigate the parameter identifiability whether uncertain parameters can be accurately estimated. The parameter identifiability was investigated in terms of the overparameterization as well as the correlation between estimation parameters. The reference office building represented in an EnergyPlus simulation model was chosen for this case study. Four parameters in the EnergyPlus simulation model were chosen to be estimated, and a meta model, or ANN model, was also developed. For investigating the parameter identifiability, the likelihood CI and CR values were introduced. In addition, biplot analysis was introduced to check the parameter correlation. The findings can be summarized as follows: • The likelihood CIs can be used to investigate the parameter identifiability. If the likelihood CI of a parameter is close to its parameter range, the parameter would not be identifiable (Tables 3 and 4). As we remove unidentifiable parameters from the parameter list, the other parameters become more identifiable. • It was found that the correlation between two parameters also influences the identifiability (Table 5 and Fig. 6). The more correlated a parameter with others, the more unidentifiable the parameter would be. As we remove the correlated parameter, other parameters become more identifiable. • The CR tends to be widespread when a parameter is not identifiable. When two parameters are linearly correlated, their posterior modes tend to follow a linear relationship. As correlated parameters are excluded, the posterior modes of remained parameters would be close to their true values (Case II in Fig. 7, Table 7). • The parameter unidentifiability increases as the parameter’s likelihood CI increases (Table 7).
D.H. Yi, D.W. Kim and C.S. Park / Energy & Buildings 198 (2019) 318–328
• Based on the aforementioned findings of this study, what follows can be suggested for Bayesian inference of building energy models: • It is a prerequisite to check the parameter identifiability using the likelihood CI and CRs as well as the biplot analysis, as presented in this paper. If a certain parameter is found to be unidentifiable, it is suggested to remove the unidentifiable parameter in the order of the highest CI. • If a parameter is found to be unidentifiable, it is necessary to examine its correlation with other parameters. • The visual assessment of likelihood CIs and CRs can help in the investigation of whether the Bayesian inference results are accurately obtained. Declaration of Competing Interest None. Acknowledgments This research was supported by Institute of Construction and Environmental Engineering at Seoul National University. The authors wish to express their gratitude for the support. The Institute of Engineering Research at Seoul National University provided research facilities for this work. References [1] G. Augenbroe, Trends in building simulation, Build. Environ. 37 (8/9) (2002) 891–902. [2] IBPSA (1987–2017), Proceedings of the IBPSA conferences (1987–2017). [3] P. de Wilde, The gap between predicted and measured energy performance of buildings: a framework for investigation, Autom. Constr 41 (2014) 40–49. [4] Y. Heo, R. Choudhary, G.A. Augenbroe, Calibration of building energy models for retrofit analysis under uncertainty, Energy Build. 47 (2012) 550–560. [5] C. Hopfe, J. Hensen, Uncertainty analysis in building performance simulation for design support, Energy Build. 43 (2011) 2798–2805. [6] W. Tian, Y. Heo, P. de Wilde, Z. Li, D. Yan, C.S. Park, X. Feng, G. Augenbroe, A review of uncertainty analysis in building energy assessment, Renew. Sustain. Energy Rev. 93 (2018) 285–301. [7] Y. Heo, G. Augenbroe, R. Choudhary, Risk analysis of energy-efficiency projects based on Bayesian calibration of building energy models, in: Building Simulation 2011, Proceedings, 2011, pp. 2579–2586. [8] A. Booth, R. Choudhary, D. Spiegelhalter, Handling uncertainty in housing stock models, Build. Environ. 48 (2012) 35–47. [9] A. Booth, R. Choudhary, Decision making under uncertainty in the retrofit analysis of the UK housing stock: implications for the Green Deal, Energy Build. 64 (2013) 292–308. [10] Y. Heo, V.M. Zavala, Gaussian process modeling for measurement and verification of building energy savings, Energy Build. 53 (2012) 7–18. [11] Y. Heo, G. Augenbroe, R. Choudhary, Quantitative risk management for energy retrofit projects, J. Build. Perform. Simul. 6 (4) (2013) 257–268. [12] M. Manfren, N. Aste, R. Moshksar, Calibration and uncertainty analysis for computer models–a meta-model based approach for integrated building energy simulation, Appl. Energy 103 (2013) 627–641. [13] W. Tian, Q. Wang, J. Song, S. Wei, Calibrating dynamic building energy models using regression model and Bayesian analysis in building retrofit projects, eSim, May 7 to 10, 2014. [14] R. Choudhary, W. Tian, Influence of district features on energy consumption in non-domestic buildings, Build. Res. Inf. 42 (1) (2014) 32–46. [15] J. Sokol, C.C. Davila, C.F. Reinhart, Validation of a Bayesian-based method for defining residential archetypes in urban building energy models, Energy Build. 134 (2017) 11–24. [16] W. Tian, R. Choudhary, A probabilistic energy model for non-domestic building sectors applied to analysis of school buildings in greater London, Energy Build. 54 (2012) 1–11. [17] F. Zhao, S.H. Lee, G. Augenbroe, Reconstructing building stock to replicate energy consumption data, Energy Build. 117 (2016) 301–312. [18] A. Chong, K.P. Lam, M. Pozzi, J. Yang, Bayesian calibration of building energy models with large datasets, Energy Build. 154 (2017) 343–355. [19] Y. Kang, M. Krarti, Bayesian-Emulator based parameter identification for calibrating energy models for existing buildings, Build. Simul. (2016) 411–428 Springer. [20] M. Kristensen, R. Choudhary, S. Petersen, Bayesian calibration of building energy models: comparison of predictive accuracy using metered utility data of different temporal resolution, Energy Procedia 122 (2017) 277–282. [21] W. Tian, S. Yang, Z. Li, S. Wei, W. Pan, Y. Liu, Identifying informative energy data in Bayesian calibration of building energy models, Energy Build. 119 (2016) 363–376.
327
[22] S. Yoon, C. Park, Impact of prior distributions of energy model inputs on prediction of building energy retrofits, in: Proceedings of the CIB World Building Congress 2016, Volume V (Advancing products and services), 2016, pp. 728–736. [23] A. Chong, K. Menberg, Guidelines for the Bayesian calibration of building energy models, Energy Build. 174 (2018) 527–547. [24] W. Hastings, Monte Carlo sampling methods using Markov chains and their applications, Biometrika 57 (1) (1970) 97–109. [25] N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller, E. Teller, Equation of state calculations by fast computing machines, J. Chem. Phys. 21 (6) (1953) 1087–1092. [26] C. Sherlock, P. Fearnhead, G.O. Roberts, The random walk Metropolis: linking theory and practice through a case study, Stat. Sci. 25 (2) (2010) 172–190. [27] J. Geweke, Bayesian inference in econometric models using Monte Carlo integration, Econometrica 57 (6) (1989) 1317–1339. [28] C.J. Geyer, Importance sampling, simulated tempering, and umbrella sampling, in: Handbook of Markov Chain Monte Carlo, Chapman & Hall/CRC, Boca Raton, 2011, pp. 295–311. [29] T.P. Minka, Expectation propagation for approximate Bayesian inference, in: Proceedings of the Seventeenth conference on Uncertainty in artificial intelligence, Morgan Kaufmann Publishers Inc., 2001, pp. 362–369. [30] A. Chong, K.P. Lam, A comparison of MCMC algorithms for the Bayesian calibration of building energy models, in: Proceedings of the 15th IBPSA Building Simulation Conference, 2017. [31] K. Menberg, Y. Heo, R. Choudhary, Efficiency and reliability of Bayesian calibration of energy supply system models, Build. Simul. (2017) (2017) 1212–1221. [32] M.D. Hoffman, A. Gelman, The No-U-turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo, J. Mach. Learn. Res. 15 (1) (2014) 1593–1623. [33] Q. Liu, D. Wang, Stein variational gradient descent: a general purpose Bayesian inference algorithm, in: Advances in Neural Information Processing Systems, 2016, pp. 2378–2386. [34] Q. Li, G. Augenbroe, J. Brown, Assessment of linear emulators in lightweight Bayesian calibration of dynamic building energy models for parameter estimation and performance prediction, Energy Build. 124 (2016) 194–202. [35] Q. Li, L. Gu, G. Augenbroe, C.J. Wu, J. Brown, Calibration of dynamic building energy models with multiple responses using Bayesian inference and linear regression models, Energy Procedia 78 (2015) 979–984. [36] S. Nagpal, C. Mueller, A. Aijazi, C.F. Reinhart, A methodology for auto-calibrating urban building energy models using surrogate modeling techniques, J. Build. Perform. Simul. 12 (1) (2018) 1–16. [37] Y. Kim, Comparative study of surrogate models for uncertainty quantification of building energy model: Gaussian process emulator vs. polynomial chaos expansion, Energy Build. 133 (2016) 46–58. [38] M. Riddle, R.T. Muehleisen, A guide to Bayesian calibration of building energy models, ASHRAE/IBPSA-USA Building Simulation Conference, 2014, 2014. [39] S. Kullback, R.A. Leibler, On information and sufficiency, Ann. Math. Stat. 22 (1) (1951) 79–86. [40] B. Rannala, Identifiability of parameters in MCMC Bayesian inference of phylogeny, Syst. Biol. 51 (5) (2002) 754–760. [41] A. Raue, C. Kreutz, T. Maiwald, J. Bachmann, M. Schilling, U. Klingmüller, J. Timmer, Structural and practical identifiability analysis of partially observed dynamical models by exploiting the profile likelihood, Bioinformatics 25 (15) (2009) 1923–1929. [42] W. Meeker, L. Escobar, Teaching about approximate confidence regions based on maximum likelihood estimation, Am. Stat. 49 (1) (1995) 48–53. [43] M. Deru, K. Field, D. Studer, K. Benne, B. Griffith, P. Torcellini, B. Liu, M. Halverson, D. Winiarski, M. Rosenberg, US Department of Energy commercial reference building models of the national building stock, 2011. [44] D. Higdon, C.S. Reese, J.D. Moulton, J.A. Vrugt, C. Fox, Posterior exploration for computationally intensive forward models, in: Handbook of Markov Chain Monte Carlo, Chapman & Hall/CRC, Boca Raton, 2011, pp. 401–418. [45] S. Press, Subjective and Objective Bayesian Statistics: Principles, Models, and Applications, John Wiley & Sons, 2009. [46] A. Gelman, A. Jakulin, M.G. Pittau, Y.-S. Su, A weakly informative default prior distribution for logistic and other regression models, Ann. Appl. Stat. 2 (4) (2008) 1360–1383. [47] D. Zwickl, M.T. Holder, Model parameterization, prior distributions, and the general time-reversible model in Bayesian phylogenetics, Syst. Biol. 53 (6) (2004) 877–888. [48] A. Rossman, T. Short, M. Parks, Bayes estimators for the continuous uniform distribution, J. Stat. Educ. 6 (3) (1998). https://doi.org/10.1080/10691898.1998. 11910622. [49] C.J. Geyer, Introduction to Markov Chain Monte Carlo, in: Handbook of Markov Chain Monte Carlo, Chapman & Hall/CRC, Boca Raton, 2011, pp. 3–48. [50] H. Haario, E. Saksman, J. Tamminen, An adaptive Metropolis algorithm, Bernoulli 7 (2) (2001) 223–242. [51] G. Roberts, J.S. Rosenthal, Examples of adaptive MCMC, J. Comput. Graph. Stat. 18 (2) (2009) 349–367. [52] J. Rosenthal, Optimal proposal distributions and adaptive MCMC, in: Handbook of Markov Chain Monte Carlo, Chapman & Hall/CRC, Boca Raton, 2011, pp. 93–111. [53] SAS institution, “Burn-in, thinning, and Markov chain samples”, 2009, Available online from: https://support.sas.com/documentation/cdl/en/statug/ 63033/HTML/default/viewer.htm#statug_introbayes_sect007.htm, accessed on 2019.04.04.
328
D.H. Yi, D.W. Kim and C.S. Park / Energy & Buildings 198 (2019) 318–328
[54] A. Patil, D. Huard, C.J. Fonnesbeck, PyMC: Bayesian stochastic modelling in Python, J. Stat. Softw. 35 (4) (2010) 1–81. [55] M. Little, W.F. Heidenreich, G. Li, Parameter identifiability and redundancy: theoretical considerations, PLoS One 5 (1) (2010) e8915. [56] T. Rothenberg, Identification in parametric models, Econometrica 39 (3) (1971) 577–591. [57] S. Silvey, Statistical Inference, Routledge, 2017. [58] R. Brun, P. Reichert, H.R. Künsch, Practical identifiability analysis of large environmental simulation models, Water Resour. Res. 37 (4) (2001) 1015–1030. [59] E.A. Catchpole, B.J. Morgan, Detecting parameter redundancy, Biometrika 84 (1) (1997) 187–196. [60] D. Cole, B. Morgan, D. Titterington, Determining the parametric structure of models, Math. Biosci. 228 (1) (2010) 16–30. [61] Y. Luo, E. Weng, X. Wu, C. Gao, X. Zhou, L. Zhang, Parameter identifiability, constraint, and equifinality in data assimilation with ecosystem models, Ecol. Appl. 19 (3) (2009) 571–574. [62] S. Marsili-Libelli, M.B. Beck, P. Brunner, B. Croke, J. Guillaume, A. Jakeman, J. Jakeman, K. Keesman, H. Stigter, Practical identifiability analysis of environmental models, in: 7th International Congress on Environmental Modelling and Software, June 15 to 19, San Diego, California, USA, 2014, pp. 681–693. [63] E. Rosero, Z.L. Yang, T. Wagener, L.E. Gulden, S. Yatheendradas, G.Y. Niu, Quantifying parameter sensitivity, interaction, and transferability in hydrologically enhanced versions of the Noah land surface model over transition zones during the warm season, J. Geophys. Res. 115 (D3) (2010). https://doi.org/10.1029/ 2009JD012035. [64] K. Gabriel, The biplot graphic display of matrices with application to principal component analysis, Biometrika 58 (3) (1971) 453–467. [65] S. Holland, Principal Components Analysis (PCA), University of Georgia, 2008 Available online from: http://strata.uga.edu/software/pdf/pcaTutorial.pdf Last accessed 01/02/2019. [66] D. Torres-Salinas, N. Robinson-Garcia, E. Jiménez-Contreras, F. Herrera, E. López-Cózar, On the use of biplot analysis for multivariate bibliometric and scientific indicators, J. Am. Soc. Inf. Sci. Technol. 64 (7) (2013) 1468–1479.
[67] J. Shlens, A tutorial on principal component analysis, arXiv preprint arXiv:1404. 1100, 2014. https://arxiv.org/abs/1404.1100, accessed on 2018.01.02. [68] L. Smith, A tutorial on principal components analysis, 2002, Available online from: https://ourarchive.otago.ac.nz/bitstream/handle/10523/7534/ OUCS- 2002- 12.pdf. Last accessed 01/02/2019. [69] M. Greenacre, Correspondence analysis, in: The Oxford Handbook of Quantitative Methods. Statistical Analyses, Oxford University press, 2013, pp. 142–153. [70] W. McKinney, pandas: a Python data analysis library, 2015, Available online from: http://pandas.sourceforge.net. [71] J. Hunter, Matplotlib: a 2D graphics environment, Comput. Sci. Eng. 9 (3) (2007) 90–95. [72] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, Scikit-learn: machine learning in Python, J. Mach. Learn. Res. 12 (Oct) (2011) 2825–2830. [73] I.T. Jolliffe, J. Cadima, Principal component analysis: a review and recent developments, Phil. Trans. R. Soc. A 374 (2065) (2016) 20150202. [74] D. Mindrila, P. Balentyne, Scatterplots and correlation, 2013, Available online from: https://www.westga.edu/academics/research/vrc/assets/docs/ scatterplots_and_correlation_notes.pdf. [75] ASHRAE, ASHRAE Handbook: Fundamentals 2013, ASHRAE, 2013. [76] CIBSE, Environmental design, The Chartered Institution of Building Services Engineers, London, 2006. [77] M. McKay, R. Beckman, W. Conover, Comparison of three methods for selecting values of input variables in the analysis of output from a computer code, Technometrics 21 (2) (1979) 239–245. [78] M. Baudin, pyDOE, 2015, Available online from: http://pythonhosted.org/ pyDOE. [79] S. Philip, T. Tran, L. Tanjuatco, Eppy: scripting language for E+, EnergyPlus (version 0.5.48) [Software-GNU Affero General Public License], 2011, Available online from: https://pypi.org/project/eppy/.
Contents lists available at ScienceDirect
Energy & Buildings journal homepage: www.elsevier.com/locate/enbuild
Parameter identifiability in Bayesian inference for building energy models Dong Hyuk Yi a, Deuk Woo Kim b, Cheol Soo Park c,∗ a
Department of Architecture and Architectural Engineering, College of Engineering, Seoul National University, 1, Gwanak-ro, Gwanak-gu, Seoul, 08826, South Korea b Department of Living and Built Environment Research, Korea Institute of Civil Engineering and Building Technology, 283, Goyang-daero, Ilsanseo-gu, Goyang, Gyeonggi, 10223, South Korea c Department of Architecture and Architectural Engineering, Institute of Construction and Environmental Engineering, Institute of Engineering Research, College of Engineering, Seoul National University, 1, Gwanak-ro, Gwanak-gu, Seoul, 08826, South Korea
a r t i c l e
i n f o
Article history: Received 15 January 2019 Revised 13 May 2019 Accepted 5 June 2019 Available online 6 June 2019 Keywords: Bayesian inference Parameter identifiability Likelihood confidence interval Likelihood confidence region Biplot
a b s t r a c t Parameter identifiability is the concept of whether uncertain parameters can be correctly estimated from the observed data. The main cause of parameter unidentifiability in Bayesian inference is known as ‘overparameterization’. In this study, the likelihood confidence interval (CI) and the likelihood confidence region (CR) were introduced for quantifying the parameter identifiability. The likelihood CI and CR can be regarded as the parameter range (one-dimensional) and parameter space (two-dimensional or higher) that can identify parameter values, respectively. For this purpose, an EnergyPlus reference office building provided by the US DOE was used in this study. Four estimation parameters in the EnergyPlus model were analyzed using the likelihood CI and CR. It was found that the closer the likelihood CI of a parameter is to the prior’s parameter range, the more unidentifiable the parameter. In addition, a biplot analysis was conducted to examine a correlation between two parameters. The more correlated a parameter is with others, the more unidentifiable the parameter. It is suggested that the visual assessment of likelihood CIs and CRs can help in investigating whether Bayesian inference results can be accurately obtained. © 2019 Elsevier B.V. All rights reserved.
1. Introduction Even though building simulation disciplines have been developed over the past several decades, the performance gap between simulation prediction and measured data is still significant [1–3]. One of the reasons for the gap is the “stochastic nature” in the dynamic behavior of building systems. Therefore, it is important to reflect such stochastic characteristics of simulation inputs in a simulation model [4,5]. Hence, Bayesian inference has been highlighted owing to its capability to estimate probability distribution of uncertain inputs through the combination of domain knowledge (i.e., prior information) and observed data (i.e., likelihood) [6]. Bayesian inference has been adopted for energy retrofit analysis [7–13] and urban-scaled building energy modeling [8,9,14–17]. However, the following issues remain unsolved: (1) observed data, (2) prior choice, (3) simulation running time and (4) parameter identifiability.
∗
Corresponding author. E-mail addresses: [email protected] (D.H. Yi), [email protected] (D.W. Kim), [email protected] (C.S. Park). https://doi.org/10.1016/j.enbuild.2019.06.012 0378-7788/© 2019 Elsevier B.V. All rights reserved.
• Observed data: This is to check whether the observed data are suitable for parameter estimation of a simulation model. It has been reported that Bayesian inference can be influenced by insufficient data [18], low precision of observed data [19], low temporal resolution of observed data [20], the lack of correlation between measured data and simulation output [21], etc. • Prior choice: The posterior inference, or the outcome of Bayesian inference, is influenced by the prior set that is usually guessed at by an analyst. In other words, the posterior can be biased by the prior distribution (e.g., Gaussian prior, uniform prior, triangular prior). In addition, the result of Bayesian inference is also affected by the prior informativeness (e.g., the prior with large variance vs. the prior with small variance) [22,23]. • Simulation running time: Since Bayesian inference is based on Monte Carlo simulation (e.g., random walk Markov chain Monte Carlo [MCMC] [24–26], importance sampling [27–29], Hamiltonian Monte Carlo [30,31], and No-U-Turn sampler [32,33]), the simulation running time is demanding. For this reason, many Bayesian inference studies have used meta models as a surrogate model (e.g., linear regression model [13,15,17,34,35], artificial neural network model [36], random forest model
D.H. Yi, D.W. Kim and C.S. Park / Energy & Buildings 198 (2019) 318–328
319
Fig. 1. Scope of this study (shown in the dotted box) vs. those of previous studies.
[36], Gaussian process model [10,12,37,38]) or normative quasisteady-state model [7–9,11]. In addition, because the simulation running time is also influenced by the amount of observed data, a recent study [18] has introduced a Kullback-Leibler divergence method [39], which can select the subset from the entire observed dataset with low information loss. • Parameter identifiability: Uncertain parameters could be wrongly estimated when the parameters to be estimated are overparameterized (i.e., too many parameters to be estimated) [20,23,31]. In such cases, some parameters of the model become unidentifiable. In addition, when parameters are correlated with each other, the degree of the identifiability is also reduced [40]. The previous studies [20,23,31] have not described (1) how the overparameterization can change the parameter’s likelihood or (2) a causality between the parameter’s likelihood and Bayesian inference results. Thus, this study can be regarded as a follow-up study of [20,23,31] and aims to investigate the parameter identifiability problem. As shown in Fig. 1, several studies have addressed simulation models [7–13,15,17,34–38], observed data [18–21], prior [22,23], MCMC [23,31] and posterior [20,23,31] in Bayesian inference. With regard to the parameter identifiability, an ‘indirect’ method has been applied, e.g., monitoring MCMC convergence [23,31] or comparing prior to posterior [20,23,31]. Rather than using the ‘indirect’ method, this study suggests a ‘direct’ method using the likelihood confidence interval (hereinafter referred to as likelihood CI) and the likelihood confidence region (hereinafter referred to as likelihood CR) as means to measure whether parameters can be estimated from the observed data [41,42]. Previous studies [22,31,40] found that the number of parameters (i.e., overparameterization) can cause an unidentifiability problem. In this regard, this paper shows how the likelihood CI and CR can vary according to the number of parameters in the model as well as how the likelihood CI and CR can be used as a measure for parameter identifiability. In this study, the authors selected a reference office building presented by the US DOE [43] and generated a meta-model using an artificial neural network (ANN). Then, we analyzed the parameter identifiability based on the likelihood CIs and CRs. The parameters’ correlation, one of the reasons for parameter identifiability, was discussed using a biplot analysis, which is a visual representation of the correlation between parameters. Finally, the impact of parameter identifiability on the Bayesian inference results was addressed.
2. Methodology 2.1. Bayesian inference and MCMC Bayesian inference is based on a conditional probability theorem, as shown in Eq. (1). The posterior (P(θ |D)) is the conditional probability of uncertain parameter θ given data D. P(θ |D) can be derived from the product of the prior (P(θ )) and the likelihood (l(D|θ )) [44]. The prior, a set of uncertain parameters, is usually subjectively assumed by an analyst [45]. The level of the analyst’s certainty for the parameter value can vary. Depending on the level of the certainty, the prior can be classified as informative prior (specifying the expected parameter value in the prior) or uninformative prior (expressing the ambiguity about the parameter value in the prior) [46]. Typical examples are Gaussian prior (informative prior) and uniform prior (uninformative prior). When the uninformative prior is applied, an equal weight is applied throughout the parameter range, so that the posterior is proportional to the likelihood [47,48]. More details on the calculation of the posterior from the prior can be found in [44,49].
P (θ |D ) ≈ P (θ )l (D|θ ),
(1)
where
θ : Uncertain parameter D: Observed data P(θ |D): Posterior P(θ ): Prior l(D|θ ): Likelihood The likelihood is expressed as a function of the observed data under each parameter value as shown in Eq. (1). Unlike the prior, the analyst’s assumption is not involved in defining the likelihood. This can be regarded as objective information in the Bayesian framework as opposed to the prior [45]. The likelihood is equivalent to the probabilistic expression of the error term (e in Eq. (2)), including both the observation error and systematic gap between the observed data and the model [44]. Assuming that the error follows the Gaussian form, the likelihood can be defined as shown in Eq. (3).
D = f ( θ ) + e,
(2)
1 T l (D|θ ) ∝ exp − (D − f (θ ) ) C −1 (D − f (θ ) ) , 2 where
(3)
320
D.H. Yi, D.W. Kim and C.S. Park / Energy & Buildings 198 (2019) 318–328
f(θ ): Model output (i.e., physical model) e: Error between the data and the model C: Covariance matrix of the likelihood For sampling, an adaptive metropolis algorithm [50–52], an MCMC algorithm, was used in this study (Eq. (4)). This algorithm is a variation of the random walk metropolis algorithm [24–26] and improves convergence by updating the covariance matrix ( n in Eq. (4)) of the subsamples generated at regular update intervals in the whole sampling process and updating it (λ n in Eq. (4)). In this study, 50,0 0 0 MCMC samples (burn-in samples: 25,0 0 0, posterior samples: 25,0 0 0, regular update intervals: 50 0 0) were generated for each Bayesian inference case, and the posterior was derived by the Gaussian kernel density estimation with the MCMC samples. Burn-in refers to discarding an initial portion of a Markov chain sample so that the effect of initial values on the posterior inference is minimized [53]. For this MCMC calculation, the Python MCMC module (PyMC 2 [54]) was used.
Yn+1 ∼ N (Xn ,
λn ),
(4)
where Yn+1 : Proposal distribution used at next period (n + 1) N: Multivariate Gaussian distribution Xn : Posterior samples λ: Scaling factor (= 2.38 / parameter dimensions) n : Covariance matrix of subsamples n (subscript): Current posterior sampling period
of the likelihood CIs and CR for two parameters (X1, X2). The likelihood CR projected on the X1 and X2 axes represents the likelihood CIs for X1 and X2, respectively. The exact method to calculate Eqs. (7) and ((8) can be found in [42].
θˆ = argmax(l (D|θ ) ),
(6)
χ 2 (θ ) = (D − f (θ ))T C −1 (D − f (θ )) ≈ −2 ln (l (D|θ ) ),
(7)
⎧ ⎫ ⎨ ⎬ l D | θ ( ) θCR −2ln = χ 2 (θ ) − χ 2 θˆ < α , ⎩ ⎭ l D|θˆ
2.2. Likelihood-based parameter identifiability analysis: CI, CR and biplot The parameter identifiability analysis investigates whether the model parameters can be estimated based on the observed data [40]. In other words, this analyzes whether the parameter values can be determined within a finite parameter space. To explain the concept of the parameter identifiability, l(D|θ ) in Eq. (1) is taken again. Let a model be completely specified by the parameter set θ , and let x be a vector of observed random variables, where x ∈ ( means all possible sets of the observations). The parameters of the model are unidentifiable if there exists some θ 1 = θ 2 yielding the same distributions (Eq. (5)) for all x ∈ [40,55–57]. In other words, a combination of two different parameters can produce the same likelihood. These unidentifiable parameters would make it difficult to estimate the other parameters’ values [40].
l (x|θ1 ) = l (x|θ2 ),
Fig. 2. Likelihood CI and CR for the parameter pair (X1, X2).
(5)
The parameter identifiability is influenced by the following factors: the number of parameters, parameter redundancy, the range (min, max) of parameters, correlation between parameters, and the quantity and quality of the observed data. Please note that these factors can cause the parameter identifiability problem individually or interconnectedly [40,41,57–63]. In this study, the likelihood CI and the likelihood CR were introduced to analyze the parameter identifiability. The likelihood CI and CR are the parameter range (one-dimensional) and the parameter space (two-dimensional or higher) that can identify the parameter likelihood values, respectively [41,42]. The parameter value (θˆ ) that maximizes the likelihood can be defined as shown in Eq. (6). According to the asymptotic sampling distribution theory [41,42], the normal distributed likelihood (l(D|θ )) can be converted to a chi-squared distribution (χ 2 (θ )) Eq. (7)). Then, a parameter subspace that satisfies the condition that the difference (χ 2 (θ ) − χ 2 (θˆ )) is less than the likelihood threshold (α ) is the likelihood CR (θ CR in Eq. (8)). Fig. 2 shows a visual representation
(8)
where
θˆ : Parameter value that maximizes the likelihood θ CR : Parameter value within the likelihood CR χ 2 (θ ): Chi-squared distribution α : Likelihood threshold (inverse cumulative distribution function of the Chi-squared distribution for α confidence level) In addition, biplot [64–66] analysis was used to examine the correlation between parameters. The biplot can visualize a dataset with more than two parameters in a low-dimensional graph (in general, two dimensions are chosen) [66]. The dimensional reduction of the dataset can be performed through a multivariate analysis based on a singular value decomposition (e.g., principal component analysis [67,68], correspondence analysis [69], etc.). In the biplot analysis, each parameter can be represented as a vector on a 2-D graph. The parameter vectors are obtained with respect to the first two principal components (hereinafter referred to as PC1 and PC2, respectively, selected in the order of their eigenvalues), each weighted by the standard deviation of that component or the square root of the corresponding eigenvalue [66]. Fig. 3 shows an example from a three-dimensional likelihood CR to a two-dimensional biplot. For the biplot analysis, Pandas [70], Matplotlib [71], and Scikit-Learn [72] modules were used. A normalization process (e.g., the min-max method and Z-score method) is required because the units/values of each parameter can be different from one another. For the dimensional reduction, it is recommended that the sum of the explanatory power (i.e., the ratio of the eigenvalue of each principal component to the eigenvalue sum of two principal components) for two principal components must be greater than 70% [73] (in Fig. 3, the sum is 58.8% + 23.6% = 82.4%). Please note that the correlation between two parameters is equal to the cosine of the angle between the vectors (Fig. 3(c)),
D.H. Yi, D.W. Kim and C.S. Park / Energy & Buildings 198 (2019) 318–328
321
Fig. 3. Biplot Example.
ing power density, and infiltration rate. The minimum, median and maximum values of the parameters were determined from [75,76]. The original model values are defined as ‘true values’. Please note that the authors did not use an on-site measured data because the on-site data are usually influenced by unknown disturbances/noises. This study was intended to suggest the need for the parameter identifiability analysis as a preliminary process. That’s why the authors used a simple ‘toy experiment’ clean model (EnergyPlus) as shown in Fig. 4. In addition, what matters for parameter identifiability is not “the number of output variables”, but “the informativeness of the measured data”, which can be checked by the likelihood CI, CRs and correlations between parameters. Since the MCMC in Bayesian inference demands significant simulation running time, the artificial neural network (ANN) model was employed as a surrogate model to the EnergyPlus model. The process to generate the surrogate model was as follows:
Fig. 4. Target building.
as shown in Eq. (9). The absolute value ranges of the correlation coefficient (i.e., cosine value between two vectors of the biplot) are classified as 0 to 0.3 (very weak), 0.3 to 0.5 (weak), 0.5 to 0.7 (moderate), and 0.7 to 1.0 (strong) according to [74]. For example, the cosine value between X1 and X2 is 0.55 (Fig. 3(c)), indicating a moderate positive correlation between X1 and X2. The cosine value between X1 and X3 is −0.28, indicating a weak negative correlation. The cosine value between X2 and X3 is −0.96, meaning a strong negative correlation. If the directions between two vectors are opposite to each other, it means a strong negative correlation.
( v1 · v2 ) Cos(v1, v2) = , (|v1||v2| )
(9)
3. Target building The EnergyPlus model for a mid-size reference office building developed by the US DOE [43] was chosen as the target building (Fig. 4). The building is located in Atlanta, Georgia, USA, and has three stories above ground, a 4982 m2 floor area, and a 33% window-wall ratio. Each floor consists of four perimeter zones and one interior zone. The target building’s plant is a furnace for heating and packaged air-conditioning units for cooling. The air handing unit of the building is a VAV system with reheat coils. The building uses gas for heating and electric energy for cooling, lighting, fan, pump, and plug loads. In this study, the authors selected both annual gas energy use (kWh) and annual electric energy use (kWh) as the model outputs, or the observed data. As shown in Table 1, the following four parameters were chosen as Bayesian inference objects: fenestration U value, fenestration SHGC, light-
• Step 1: perform Latin hypercube sampling [77] based on Table 1. The number of sampled EnergyPlus simulation models is 10 0 0. • Step 2: Based on the generated EnergyPlus simulation models (Step 1), conduct a series of EnergyPlus simulation runs. • Step 3: Collect simulation results such as annual electric energy and gas uses. • Step 4: Construct ANN models using inputs and outputs of the EnergyPlus models; 750 and 250 EnergyPlus simulation runs were used for training and test, respectively. For steps 1 to 4, the following Python modules were used: pyDOE [78] for Latin hypercube sampling, Eppy [79] for generating EnergyPlus models, and scikit-learn [72] for training the ANN model (optimization solver: Adam, activation function: ReLU function, the number of hidden layers: 1, the number of neurons: 50). The generated ANN models are sufficiently accurate when compared to EnergyPlus model outputs (MBE: −0.01%, CVRMSE: 0.13%, R2 : 0.97 for gas; MBE: 0.02%, and CVRMSE: 0.35%, R2 : 0.98 for electric energy) (Fig. 5). In this study, the ANN model was saved as a class variable for Python. The authors used the Pickle module (object serialization module for Python) to export the ANN class variable to a binary file that can be read in Python for the purpose of using the ANN model in MCMC iteration. 4. Results Three comparison cases were made, as shown in Table 2. The baseline case estimates all four parameters, presented in Table 1. 30 0,0 0 0 Monte Carlo simulation runs using Latin hypercube sampling were completed for each case in Table 2. The computational
322
D.H. Yi, D.W. Kim and C.S. Park / Energy & Buildings 198 (2019) 318–328 Table 1 List of estimation parameters. Parameters
True value
Min
Median
Max
References
Fenestration U value (W/m2 K) Fenestration SHGC (0.0–1.0) Lighting power density (W/m2 ) Infiltration rate (air changes per hour, 1/h)
4.09 0.25 16.89 0.35
1.03 0.13 7.0 0.15
3.705 0.495 14.0 0.45
6.38 0.86 21.0 0.75
[75] [75] [75] [76]
Fig. 5. Comparison between EnergyPlus and ANN models (250 test data used).
Fig. 6. Biplot analysis.
Table 2 Comparison cases. Number of parameters Baseline Case I Case II
4 3 2
costs of using the method proposed in this study were as follows: 308.3, 304.7, and 316.9 s for Baseline, Case I, and Case II, respectively (CPU: AMD Ryzen 5 1600X 3.60 GHz, RAM: 8 GB). When one parameter is removed, its true value in Table 1 was used for Cases I and II. In this study, 95% likelihood CIs and 95% likelihood CRs were used. 4.1. Likelihood CI and biplot analysis Table 3 shows the ratios of likelihood CIs for four parameters to their parameter ranges, Table 4 shows a visual comparison of the likelihood CIs for four parameters to their ranges, Table 5 shows
the degrees of correlation between two parameters, and Fig. 6 shows the biplot analyses for the three cases. The results are summarized as follows: • Baseline case: Except lighting power density, the ratios of the likelihood CIs of the other three parameters to their parameter ranges are greater than 0.90 (Table 3). In other words, the three parameters are difficult to estimate because their likelihood CIs are close to their parameter ranges. Particularly, the likelihood CI of the infiltration rate is the same as its parameter range, and its likelihood shape is flat, indicating that its likelihood over its parameter range is equally distributed (Table 4). As shown in Table 5 and Fig. 6, strong correlations exist between ‘Fenestration U value and Fenestration SHGC (0.83)’ and ‘Fenestration SHGC and Lighting power density (−0.93)’. Because Fenestration SHGC is a shared parameter included in the aforementioned parameter pairs, it can be inferred that Fenestration SHGC may cause an identifiability problem due to its correlation with other parameters [40]. In contrast, Fenestration U value, which is difficult to identify (B/A = 0.93, Table 3), has a moderate or weak correlation (−0.56 with lighting power
D.H. Yi, D.W. Kim and C.S. Park / Energy & Buildings 198 (2019) 318–328
323
Table 3 Ratios of parameters’ likelihood CIs to parameter ranges. Parameter
Parameter range (A) (Table 1)
Cases
Likelihood CI (B) (refer to Table 4)
B/A
Fenestration U value (W/m2 K)
6.38–1.03 = 5.35
Baseline Case I Case II
6.38–1.39 = 4.99 6.38–1.83 = 4.55 6.38–2.23 = 4.15
0.93 0.85 0.78
Fenestration SHGC (0–1)
0.86–0.13 = 0.73
Baseline Case I Case II
0.84–0.13 = 0.71 0.83–0.13 = 0.70 n/a
0.97 0.96 n/a
Lighting power density (W/m2 )
21.0–7.0 = 14.00
Baseline Case I Case II
20.88–11.58 = 9.30 20.64–12.26 = 8.38 20.10–13.49 = 6.61
0.66 0.60 0.47
Infiltration rate (ACH)
0.75–0.15 = 0.60
Baseline Case I Case II
0.75–0.15 = 0.60 n/a n/a
1.00 n/a n/a
The bold numbers indicate that the corresponding parameter is excluded in the next case study. Table 4 Comparison of likelihood CIs to parameter ranges.
1.0
0.5
3.70
6.38
1.0
0.5
0.0 1.03
Likelihood CI (0.13, 0.84)
1.0
0.5
0.0 0.13
0.50
0.86
0.5
0.0 21.0
Log likelihood (Normalized)
Log likelihood (Normalized)
Likelihood CI (11.58, 20.88)
Lighting power density Log likelihood (Normalized)
0.5
0.0 1.03
3.70
6.38
Fenestration U value
Likelihood CI (0.13, 0.83)
excluded in Case II
0.5
0.0 0.13
0.50
0.86
Fenestration SHGC
1.0
14.0
1.0
1.0
Fenestration SHGC
7.0
6.38
Likelihood CI (2.23, 6.38)
Fenestration U value Log likelihood (Normalized)
Log likelihood (Normalized)
Fenestration U value
3.70
Likelihood CI (12.26, 20.64)
1.0
0.5
0.0
14.0
7.0
21.0
Lighting power density
Log likelihood (Normalized)
0.0 1.03
Case II
Likelihood CI (1.83, 6.38)
Log likelihood (Normalized)
Case I Log likelihood (Normalized)
Log likelihood (Normalized)
Baseline Likelihood CI (1.39, 6.38)
Likelihood CI (13.49, 20.10)
1.0
0.5
0.0 7.0
14.0
21.0
Lighting power density
Likelihood CI (0.15, 0.75) 1.0
excluded in Case I
0.5
0.0 0.15
0.45
Infiltration rate
excluded in Case II
0.75
density, −0.34 with infiltration rate) except for the Fenestration SHGC. The infiltration rate has a moderate, weak, or very weak correlation (0.24, −0.34, −0.59) with the other three parameters. Hence, it can be said that Fenestration U value and Infiltration rate cause the identifiability problem independently when compared with the Fenestration SHGC. Based on the above
analysis, it was decided that infiltration rate should be excluded from the list of estimation parameters for Case I. • Case I: As shown in Table 3, the Fenestration U value (0.93 → 0.85) and Fenestration SHGC value (0.97 → 0.96) are found to have high ratios of the likelihood CIs to their parameter ranges, indicating that the two parameters are difficult
324
D.H. Yi, D.W. Kim and C.S. Park / Energy & Buildings 198 (2019) 318–328 Table 5 The degrees of correlation between two parameters (cosine values between two parameters). Two vectors
Baseline
Case I
Case II
Fenestration U value vs. Fenestration SHGC Fenestration U value vs. Lighting power density Fenestration U value vs. Infiltration rate Fenestration SHGC vs. Lighting power density Fenestration SHGC vs. Infiltration rate Lighting power density vs. Infiltration rate
0.83 −0.56 −0.34 −0.93 0.24 −0.59
0.55 −0.28 excluded −0.96 excluded excluded
excluded −0.21 excluded excluded excluded excluded
Table 6 Six priors (P1–P6). mean Parameters Fenestration U value (W/m2 K) Fenestration SHGC (0–1) Lighting power density (W/m2 ) Infiltration rate (ACH) ∗
P1
P2
P3
P4
3.705 0.495 14.0 0.45
5.80 0.69 10.77 0.72
4.12 0.82 15.80 0.33
2.68 0.36 9.31 0.51
P5
P6
4.36 0.53 14.15 0.59
1.61 0.21 20.08 0.16
Std. Dev.
min∗
max∗
1.37 0.19 3.57 0.15
1.03 0.13 7.00 0.15
6.38 0.86 21.0 0.75
The minimum and maximum values are the same as those in Table 1.
Table 7 Posterior modes. True value
Parameter range (A) (from Table 1)
Fenestration U value (W/m2 K)
4.09
6.38–1.03 = 5.35
Fenestration SHGC (0–1)
0.25
Lighting power density (W/m2 ) Infiltration Rate (ACH)
Parameter
Posterior modes Case
(Max – min) of six posterior modes (C)
C/A
B/A (from Table 3)
3.51 3.38 3.38
2.68 2.78 1.71
0.50 0.52 0.32
0.93 0.85 0.78
0.44 0.48 n/a
0.16 0.15 n/a
0.57 0.47 n/a
0.78 0.64 n/a
0.97 0.96 Excluded
14.80 15.22 15.99
14.94 15.06 16.53
17.56 17.48 17.57
3.80 2.91 2.00
0.27 0.21 0.14
0.66 0.60 0.47
0.57 n/a n/a
0.63 n/a n/a
0.24 n/a n/a
0.48 n/a n/a
0.80 n/a n/a
1.00 Excluded Excluded
P1
P2
P3
P4
P5
P6
Baseline Case I Case II
4.26 3.98 3.98
6.19 6.16 5.09
4.97 4.77 4.15
3.65 3.89 3.60
4.50 5.02 4.10
0.86–0.13 = 0.73
Baseline Case I Case II
0.44 0.37 n/a
0.64 0.62 n/a
0.61 0.61 n/a /
0.39 0.36 n/a
16.89
21.0–7.0 = 14.00
Baseline Case I Case II
14.96 15.35 16.37
13.76 14.57 15.57
14.99 15.23 16.60
0.35
0.75–0.15 = 0.60
Baseline Case I Case II
0.48 n/a n/a
0.72 n/a n/a
0.38 n/a n/a
to estimate. The likelihood shapes of the two parameters are similar to those of Case I (Table 4). According to Table 5 and Fig. 6, there is a strong correlation (−0.96) between ‘Fenestration SHGC and Lighting power density’. It can be inferred that the identifiability problem exists due to the correlation between Fenestration SHGC and Lighting power density. Therefore, it was determined that Fenestration SHGC should be excluded for Case II. • Case II: After two parameters (Infiltration rate and Fenestration SHGC) were excluded from the list of estimation parameters, Case II was conducted. The ratios of the likelihood CIs of Fenestration U value and lighting power density to their parameter ranges are reduced (0.93 → 0.78, 0.66 → 0.47). In addition, distinct forms of peaks are observed in the likelihood shapes of the two parameters (Fenestration U value and lighting power density), compared to those of the baseline case and Case I (Table 4). The cosine value between the two parameters (Fenestration U value and Lighting power density) is −0.21, implying that the correlation between the two parameters is very weak [74] (Table 5, Fig. 6). It is expected that the estimation for the two parameters would be correct partly because the parameter subspace is reduced from the entire parameter space. This result will be substantiated in Section 4.2. 4.2. Likelihood CR analysis To investigate how posterior can be biased depending on prior, the following six priors (P1-P6) were made as shown in Table 6.
The six priors were assumed to follow a Gaussian distribution. P1 represents the median values in Table 1, while P2–P6 were sampled using Latin hypercube sampling [77]. The standard deviations of P1-P6 for each parameter were assumed to be equal to (mean of P1 – min of P1) /1.96. Table 7 tabulates the six posterior modes, the difference between their maximum and minimum modes (denoted as ‘C’), and the ratio of C to A (parameter range). As explained in Section 4.1, the ratio (B/A) can be regarded as a measure for the parameter ‘unidentifiability’. The greater the ratio (B/A) is, the more unidentifiable the parameter. Similar findings can be drawn by comparing the ratio (C/A) to (B/A). It can be concluded that as the ratio (C/A) increases, it is more likely that the Bayesian results are farther away from the true value. Fig. 7 shows the true values (a blue star), their posterior modes (green circles), and their prior means (white squares) on the likelihood CRs (red dots) in a 2-D space. Because the correlations between ‘Fenestration U value and Fenestration SHGC’ are positively linear (the cosine values are 0.83 in Baseline case and 0.55 in Case I, Table 5), the likelihood CRs and posterior modes of the two are also positively related, as shown in Fig. 7(a). In contrast, ‘Fenestration SHGC and Lighting power density’ are negatively correlated (the cosine values between the two are −0.93 in Baseline case and −0.96 in Case I, Table 5), and the likelihood CRs and posterior modes of the two are found to be negatively linear (Fig. 7(b)). In other words, the posterior modes of two parameters tend to follow the correlation between the two.
D.H. Yi, D.W. Kim and C.S. Park / Energy & Buildings 198 (2019) 318–328
Fig. 7. Likelihood CRs for estimation parameters. (For interpretation of the references to color in the text, the reader is referred to the web version of this article.)
325
326
D.H. Yi, D.W. Kim and C.S. Park / Energy & Buildings 198 (2019) 318–328
Fig. 8. The parameter unidentifiability increases as the likelihood CI increases. (For interpretation of the references to color in the text, the reader is referred to the web version of this article.)
It is noteworthy that the greater the ratios of the parameters’ likelihood CIs to the parameter ranges (Table 3), the wider the parameters’ likelihood CRs and posterior modes were estimated. For example, the likelihood CR and posterior modes for the Infiltration rate, of which the ratio is 1.00 (Table 3), are widely spread over its parameter range (Fig. 7(c), (e) and (f)). Similarly, the likelihood CR and posterior modes of Fenestration SHGC, of which the ratios are 0.97 in the Baseline case and 0.96 in Case I (Table 3), are widely placed across its parameter range (Fig. 7(a), (d) and (e)). From the Baseline case via Case I to Case II, the likelihood CRs of the parameters are gradually reduced owing to excluded parameters as well as removed parameter correlations. For this reason, the posterior modes of lighting power density and Fenestration U value tend to converge to their true values (Case II in Fig. 7(b), Table 7). In particular, owing to the exclusion of Fenestration SHGC (the parameter that caused the strong correlation with other parameters, Table 5), a significant change in the shape of the likelihood CR and the posterior modes occurred from Case I to Case II (Fig. 7(b)). The black solid line in Fig. 8 shows the relationship between the two ratios (B/A and C/A, Table 7). As B/A increases, the difference between the minimum and maximum posterior modes increases (corresponding to the vertical height of the dotted red circles in Fig. 8), indicating that the Bayesian inference results are widespread. In other words, the parameters might be incorrectly estimated, not sufficiently close to the true value (Table 7, Fig. 8). As presented in Table 7 and Fig. 8, it can be inferred that the greater the ratio of the parameter’s likelihood CI to the parameter range is, the greater the parameter unidentifiability. The degree of the uncertainty in posterior modes of the parameters increases in proportion to the ratio of the parameter’s likelihood CI to the parameter range. Conclusively, this study shows that the likelihood CIs and CRs can be used to investigate the parameter identifiability in the Bayesian inference problems, as illustrated in Fig. 8.
5. Conclusion In this paper, the authors suggested a methodology to investigate the parameter identifiability whether uncertain parameters can be accurately estimated. The parameter identifiability was investigated in terms of the overparameterization as well as the correlation between estimation parameters. The reference office building represented in an EnergyPlus simulation model was chosen for this case study. Four parameters in the EnergyPlus simulation model were chosen to be estimated, and a meta model, or ANN model, was also developed. For investigating the parameter identifiability, the likelihood CI and CR values were introduced. In addition, biplot analysis was introduced to check the parameter correlation. The findings can be summarized as follows: • The likelihood CIs can be used to investigate the parameter identifiability. If the likelihood CI of a parameter is close to its parameter range, the parameter would not be identifiable (Tables 3 and 4). As we remove unidentifiable parameters from the parameter list, the other parameters become more identifiable. • It was found that the correlation between two parameters also influences the identifiability (Table 5 and Fig. 6). The more correlated a parameter with others, the more unidentifiable the parameter would be. As we remove the correlated parameter, other parameters become more identifiable. • The CR tends to be widespread when a parameter is not identifiable. When two parameters are linearly correlated, their posterior modes tend to follow a linear relationship. As correlated parameters are excluded, the posterior modes of remained parameters would be close to their true values (Case II in Fig. 7, Table 7). • The parameter unidentifiability increases as the parameter’s likelihood CI increases (Table 7).
D.H. Yi, D.W. Kim and C.S. Park / Energy & Buildings 198 (2019) 318–328
• Based on the aforementioned findings of this study, what follows can be suggested for Bayesian inference of building energy models: • It is a prerequisite to check the parameter identifiability using the likelihood CI and CRs as well as the biplot analysis, as presented in this paper. If a certain parameter is found to be unidentifiable, it is suggested to remove the unidentifiable parameter in the order of the highest CI. • If a parameter is found to be unidentifiable, it is necessary to examine its correlation with other parameters. • The visual assessment of likelihood CIs and CRs can help in the investigation of whether the Bayesian inference results are accurately obtained. Declaration of Competing Interest None. Acknowledgments This research was supported by Institute of Construction and Environmental Engineering at Seoul National University. The authors wish to express their gratitude for the support. The Institute of Engineering Research at Seoul National University provided research facilities for this work. References [1] G. Augenbroe, Trends in building simulation, Build. Environ. 37 (8/9) (2002) 891–902. [2] IBPSA (1987–2017), Proceedings of the IBPSA conferences (1987–2017). [3] P. de Wilde, The gap between predicted and measured energy performance of buildings: a framework for investigation, Autom. Constr 41 (2014) 40–49. [4] Y. Heo, R. Choudhary, G.A. Augenbroe, Calibration of building energy models for retrofit analysis under uncertainty, Energy Build. 47 (2012) 550–560. [5] C. Hopfe, J. Hensen, Uncertainty analysis in building performance simulation for design support, Energy Build. 43 (2011) 2798–2805. [6] W. Tian, Y. Heo, P. de Wilde, Z. Li, D. Yan, C.S. Park, X. Feng, G. Augenbroe, A review of uncertainty analysis in building energy assessment, Renew. Sustain. Energy Rev. 93 (2018) 285–301. [7] Y. Heo, G. Augenbroe, R. Choudhary, Risk analysis of energy-efficiency projects based on Bayesian calibration of building energy models, in: Building Simulation 2011, Proceedings, 2011, pp. 2579–2586. [8] A. Booth, R. Choudhary, D. Spiegelhalter, Handling uncertainty in housing stock models, Build. Environ. 48 (2012) 35–47. [9] A. Booth, R. Choudhary, Decision making under uncertainty in the retrofit analysis of the UK housing stock: implications for the Green Deal, Energy Build. 64 (2013) 292–308. [10] Y. Heo, V.M. Zavala, Gaussian process modeling for measurement and verification of building energy savings, Energy Build. 53 (2012) 7–18. [11] Y. Heo, G. Augenbroe, R. Choudhary, Quantitative risk management for energy retrofit projects, J. Build. Perform. Simul. 6 (4) (2013) 257–268. [12] M. Manfren, N. Aste, R. Moshksar, Calibration and uncertainty analysis for computer models–a meta-model based approach for integrated building energy simulation, Appl. Energy 103 (2013) 627–641. [13] W. Tian, Q. Wang, J. Song, S. Wei, Calibrating dynamic building energy models using regression model and Bayesian analysis in building retrofit projects, eSim, May 7 to 10, 2014. [14] R. Choudhary, W. Tian, Influence of district features on energy consumption in non-domestic buildings, Build. Res. Inf. 42 (1) (2014) 32–46. [15] J. Sokol, C.C. Davila, C.F. Reinhart, Validation of a Bayesian-based method for defining residential archetypes in urban building energy models, Energy Build. 134 (2017) 11–24. [16] W. Tian, R. Choudhary, A probabilistic energy model for non-domestic building sectors applied to analysis of school buildings in greater London, Energy Build. 54 (2012) 1–11. [17] F. Zhao, S.H. Lee, G. Augenbroe, Reconstructing building stock to replicate energy consumption data, Energy Build. 117 (2016) 301–312. [18] A. Chong, K.P. Lam, M. Pozzi, J. Yang, Bayesian calibration of building energy models with large datasets, Energy Build. 154 (2017) 343–355. [19] Y. Kang, M. Krarti, Bayesian-Emulator based parameter identification for calibrating energy models for existing buildings, Build. Simul. (2016) 411–428 Springer. [20] M. Kristensen, R. Choudhary, S. Petersen, Bayesian calibration of building energy models: comparison of predictive accuracy using metered utility data of different temporal resolution, Energy Procedia 122 (2017) 277–282. [21] W. Tian, S. Yang, Z. Li, S. Wei, W. Pan, Y. Liu, Identifying informative energy data in Bayesian calibration of building energy models, Energy Build. 119 (2016) 363–376.
327
[22] S. Yoon, C. Park, Impact of prior distributions of energy model inputs on prediction of building energy retrofits, in: Proceedings of the CIB World Building Congress 2016, Volume V (Advancing products and services), 2016, pp. 728–736. [23] A. Chong, K. Menberg, Guidelines for the Bayesian calibration of building energy models, Energy Build. 174 (2018) 527–547. [24] W. Hastings, Monte Carlo sampling methods using Markov chains and their applications, Biometrika 57 (1) (1970) 97–109. [25] N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller, E. Teller, Equation of state calculations by fast computing machines, J. Chem. Phys. 21 (6) (1953) 1087–1092. [26] C. Sherlock, P. Fearnhead, G.O. Roberts, The random walk Metropolis: linking theory and practice through a case study, Stat. Sci. 25 (2) (2010) 172–190. [27] J. Geweke, Bayesian inference in econometric models using Monte Carlo integration, Econometrica 57 (6) (1989) 1317–1339. [28] C.J. Geyer, Importance sampling, simulated tempering, and umbrella sampling, in: Handbook of Markov Chain Monte Carlo, Chapman & Hall/CRC, Boca Raton, 2011, pp. 295–311. [29] T.P. Minka, Expectation propagation for approximate Bayesian inference, in: Proceedings of the Seventeenth conference on Uncertainty in artificial intelligence, Morgan Kaufmann Publishers Inc., 2001, pp. 362–369. [30] A. Chong, K.P. Lam, A comparison of MCMC algorithms for the Bayesian calibration of building energy models, in: Proceedings of the 15th IBPSA Building Simulation Conference, 2017. [31] K. Menberg, Y. Heo, R. Choudhary, Efficiency and reliability of Bayesian calibration of energy supply system models, Build. Simul. (2017) (2017) 1212–1221. [32] M.D. Hoffman, A. Gelman, The No-U-turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo, J. Mach. Learn. Res. 15 (1) (2014) 1593–1623. [33] Q. Liu, D. Wang, Stein variational gradient descent: a general purpose Bayesian inference algorithm, in: Advances in Neural Information Processing Systems, 2016, pp. 2378–2386. [34] Q. Li, G. Augenbroe, J. Brown, Assessment of linear emulators in lightweight Bayesian calibration of dynamic building energy models for parameter estimation and performance prediction, Energy Build. 124 (2016) 194–202. [35] Q. Li, L. Gu, G. Augenbroe, C.J. Wu, J. Brown, Calibration of dynamic building energy models with multiple responses using Bayesian inference and linear regression models, Energy Procedia 78 (2015) 979–984. [36] S. Nagpal, C. Mueller, A. Aijazi, C.F. Reinhart, A methodology for auto-calibrating urban building energy models using surrogate modeling techniques, J. Build. Perform. Simul. 12 (1) (2018) 1–16. [37] Y. Kim, Comparative study of surrogate models for uncertainty quantification of building energy model: Gaussian process emulator vs. polynomial chaos expansion, Energy Build. 133 (2016) 46–58. [38] M. Riddle, R.T. Muehleisen, A guide to Bayesian calibration of building energy models, ASHRAE/IBPSA-USA Building Simulation Conference, 2014, 2014. [39] S. Kullback, R.A. Leibler, On information and sufficiency, Ann. Math. Stat. 22 (1) (1951) 79–86. [40] B. Rannala, Identifiability of parameters in MCMC Bayesian inference of phylogeny, Syst. Biol. 51 (5) (2002) 754–760. [41] A. Raue, C. Kreutz, T. Maiwald, J. Bachmann, M. Schilling, U. Klingmüller, J. Timmer, Structural and practical identifiability analysis of partially observed dynamical models by exploiting the profile likelihood, Bioinformatics 25 (15) (2009) 1923–1929. [42] W. Meeker, L. Escobar, Teaching about approximate confidence regions based on maximum likelihood estimation, Am. Stat. 49 (1) (1995) 48–53. [43] M. Deru, K. Field, D. Studer, K. Benne, B. Griffith, P. Torcellini, B. Liu, M. Halverson, D. Winiarski, M. Rosenberg, US Department of Energy commercial reference building models of the national building stock, 2011. [44] D. Higdon, C.S. Reese, J.D. Moulton, J.A. Vrugt, C. Fox, Posterior exploration for computationally intensive forward models, in: Handbook of Markov Chain Monte Carlo, Chapman & Hall/CRC, Boca Raton, 2011, pp. 401–418. [45] S. Press, Subjective and Objective Bayesian Statistics: Principles, Models, and Applications, John Wiley & Sons, 2009. [46] A. Gelman, A. Jakulin, M.G. Pittau, Y.-S. Su, A weakly informative default prior distribution for logistic and other regression models, Ann. Appl. Stat. 2 (4) (2008) 1360–1383. [47] D. Zwickl, M.T. Holder, Model parameterization, prior distributions, and the general time-reversible model in Bayesian phylogenetics, Syst. Biol. 53 (6) (2004) 877–888. [48] A. Rossman, T. Short, M. Parks, Bayes estimators for the continuous uniform distribution, J. Stat. Educ. 6 (3) (1998). https://doi.org/10.1080/10691898.1998. 11910622. [49] C.J. Geyer, Introduction to Markov Chain Monte Carlo, in: Handbook of Markov Chain Monte Carlo, Chapman & Hall/CRC, Boca Raton, 2011, pp. 3–48. [50] H. Haario, E. Saksman, J. Tamminen, An adaptive Metropolis algorithm, Bernoulli 7 (2) (2001) 223–242. [51] G. Roberts, J.S. Rosenthal, Examples of adaptive MCMC, J. Comput. Graph. Stat. 18 (2) (2009) 349–367. [52] J. Rosenthal, Optimal proposal distributions and adaptive MCMC, in: Handbook of Markov Chain Monte Carlo, Chapman & Hall/CRC, Boca Raton, 2011, pp. 93–111. [53] SAS institution, “Burn-in, thinning, and Markov chain samples”, 2009, Available online from: https://support.sas.com/documentation/cdl/en/statug/ 63033/HTML/default/viewer.htm#statug_introbayes_sect007.htm, accessed on 2019.04.04.
328
D.H. Yi, D.W. Kim and C.S. Park / Energy & Buildings 198 (2019) 318–328
[54] A. Patil, D. Huard, C.J. Fonnesbeck, PyMC: Bayesian stochastic modelling in Python, J. Stat. Softw. 35 (4) (2010) 1–81. [55] M. Little, W.F. Heidenreich, G. Li, Parameter identifiability and redundancy: theoretical considerations, PLoS One 5 (1) (2010) e8915. [56] T. Rothenberg, Identification in parametric models, Econometrica 39 (3) (1971) 577–591. [57] S. Silvey, Statistical Inference, Routledge, 2017. [58] R. Brun, P. Reichert, H.R. Künsch, Practical identifiability analysis of large environmental simulation models, Water Resour. Res. 37 (4) (2001) 1015–1030. [59] E.A. Catchpole, B.J. Morgan, Detecting parameter redundancy, Biometrika 84 (1) (1997) 187–196. [60] D. Cole, B. Morgan, D. Titterington, Determining the parametric structure of models, Math. Biosci. 228 (1) (2010) 16–30. [61] Y. Luo, E. Weng, X. Wu, C. Gao, X. Zhou, L. Zhang, Parameter identifiability, constraint, and equifinality in data assimilation with ecosystem models, Ecol. Appl. 19 (3) (2009) 571–574. [62] S. Marsili-Libelli, M.B. Beck, P. Brunner, B. Croke, J. Guillaume, A. Jakeman, J. Jakeman, K. Keesman, H. Stigter, Practical identifiability analysis of environmental models, in: 7th International Congress on Environmental Modelling and Software, June 15 to 19, San Diego, California, USA, 2014, pp. 681–693. [63] E. Rosero, Z.L. Yang, T. Wagener, L.E. Gulden, S. Yatheendradas, G.Y. Niu, Quantifying parameter sensitivity, interaction, and transferability in hydrologically enhanced versions of the Noah land surface model over transition zones during the warm season, J. Geophys. Res. 115 (D3) (2010). https://doi.org/10.1029/ 2009JD012035. [64] K. Gabriel, The biplot graphic display of matrices with application to principal component analysis, Biometrika 58 (3) (1971) 453–467. [65] S. Holland, Principal Components Analysis (PCA), University of Georgia, 2008 Available online from: http://strata.uga.edu/software/pdf/pcaTutorial.pdf Last accessed 01/02/2019. [66] D. Torres-Salinas, N. Robinson-Garcia, E. Jiménez-Contreras, F. Herrera, E. López-Cózar, On the use of biplot analysis for multivariate bibliometric and scientific indicators, J. Am. Soc. Inf. Sci. Technol. 64 (7) (2013) 1468–1479.
[67] J. Shlens, A tutorial on principal component analysis, arXiv preprint arXiv:1404. 1100, 2014. https://arxiv.org/abs/1404.1100, accessed on 2018.01.02. [68] L. Smith, A tutorial on principal components analysis, 2002, Available online from: https://ourarchive.otago.ac.nz/bitstream/handle/10523/7534/ OUCS- 2002- 12.pdf. Last accessed 01/02/2019. [69] M. Greenacre, Correspondence analysis, in: The Oxford Handbook of Quantitative Methods. Statistical Analyses, Oxford University press, 2013, pp. 142–153. [70] W. McKinney, pandas: a Python data analysis library, 2015, Available online from: http://pandas.sourceforge.net. [71] J. Hunter, Matplotlib: a 2D graphics environment, Comput. Sci. Eng. 9 (3) (2007) 90–95. [72] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, Scikit-learn: machine learning in Python, J. Mach. Learn. Res. 12 (Oct) (2011) 2825–2830. [73] I.T. Jolliffe, J. Cadima, Principal component analysis: a review and recent developments, Phil. Trans. R. Soc. A 374 (2065) (2016) 20150202. [74] D. Mindrila, P. Balentyne, Scatterplots and correlation, 2013, Available online from: https://www.westga.edu/academics/research/vrc/assets/docs/ scatterplots_and_correlation_notes.pdf. [75] ASHRAE, ASHRAE Handbook: Fundamentals 2013, ASHRAE, 2013. [76] CIBSE, Environmental design, The Chartered Institution of Building Services Engineers, London, 2006. [77] M. McKay, R. Beckman, W. Conover, Comparison of three methods for selecting values of input variables in the analysis of output from a computer code, Technometrics 21 (2) (1979) 239–245. [78] M. Baudin, pyDOE, 2015, Available online from: http://pythonhosted.org/ pyDOE. [79] S. Philip, T. Tran, L. Tanjuatco, Eppy: scripting language for E+, EnergyPlus (version 0.5.48) [Software-GNU Affero General Public License], 2011, Available online from: https://pypi.org/project/eppy/.