History: Bianca Maria Colosimo served as the senior editor for this article.
Funding: This work was partially supported by the National Science Foundation [Grants CMMI-1662553, CMMI-2226348, and CBET-1804321].
Data Ethics & Reproducibility Note: The code capsule is available on Code Ocean at https://codeocean.com/capsule/8623151/tree/v1 and in the e-Companion to this article (available at https://doi.org/10.1287/ijds.2023.0029 ).
]]>Recent advancements in computing have enabled the adoption of computer models as an essential tool to design and control systems. In general, a computer model consists of a set of mathematical functions that are based on complex physics-based first principals in order to closely resemble a real-world system. A computer model is usually utilized with a simulator. A simulator receives input variables that are typically observable or controllable, such as operational conditions, and it generates outputs through simulation. Besides the input variables, additional parameters need to be specified a priori to run simulation. These parameters are oftentimes not observable, thus their values should be estimated using either physical laws or data. When physical laws that identify 2 Article submitted to INFORMS Journal on Data Science; manuscript no. IJDS-2022-0008 appropriate parameter values are not available, one can use a data-driven approach using data collected from the physical system. Such a data-driven procedure is called parameter calibration in the literature (Kennedy and O'Hagan 2001).
Over the last couple of decades, parameter calibration has been studied in many fields such as biology, chemistry, climatology, and engineering (Santner et al. 2018). We provide a detailed literature review in Section 3. Existing studies, by and large, utilize limited datasets when physical trials are limited and computer models are expensive to run. To accommodate data scarcity, they build surrogate models for computer models and calibration is performed with surrogate models.
Typically, they generate data from a computer model a priori through a fixed set of simulation runs.
With the resulting data, they build a response surface (or emulator) for the computer model output using Gaussian processes (GPs) or other statistical models (Higdon et al. 2004, Tuo et al. 2021).
Consequently, calibration accuracy highly depends on the emulator accuracy. When the computer model data are generated at less informative design points, e.g., values far from (unknown) true parameter values, the emulator may not accurately characterize the response surface near the true parameter values, leading to inaccurate calibration results.
Recently, Liu et al. (2021) discussed a necessity to develop a new calibration procedure for big data settings. In many operational systems, thanks to the recent advances in sensing technology and data acquisition systems, massive amounts of observational data from an actual physical system have become available. Further, a sufficiently large number of simulation data can be generated through a medium-or low-fidelity simulator or using advanced computing facilities, e.g., supercomputers. When a large number of field observations and simulation data are available, Liu et al. (2021) presented a new approach to directly utilize data without building emulators.
The aforementioned existing literature calibrates all parameters in a similar manner despite of heterogeneous data requirement. That is, they calibrate all parameters at once with the same dataset. In some applications, however, some parameters are associated with certain specific operating conditions or a subset of datasets. Consider a building energy model (BEM). Among many parameters that need to be calibrated, some parameters are season-dependent, and thus they should be employed for simulating building energy operations only in specific seasons (Xu et al. 2021).
For example, "cooling supply air flow rate" is a cooling-season parameter, whereas "heating supply air flow rate" is associated with a heating season. On the other hand, other parameters, e.g., "fan total efficiency", are global parameters that need to be employed throughout the year. This BEM parameter calibration, as our motivating application, will be explained in more detail in Section 2.
From a data science point of view, these season-dependent (or block-dependent) parameters should be calibrated using data (both field operational and simulation-generated data) associated with the corresponding season (or block) only. That is, one needs to use a subset of data. When 3 the required dataset for each group of parameters is exclusive from one another, one can calibrate parameters separately for each group. However, when the subsets for different groups overlap, prior methods become inappropriate.
This study develops a new calibration approach where parameters can be divided into multiple groups and each group of parameters is associated with a subset (or block) of the entire dataset. Specifically, we estimate the parameters by minimizing multiple loss functions, each of which is associated with each group of parameters. Our proposed approach is useful when large computer experiments can be conducted, e.g., from a low-to medium-fidelity computer model.
We summarize and highlight the key contribution of this paper as follows.
• Our approach estimates unobservable parameters employed in a computer model where parameters are associated with different subsets of observational data that possibly overlap with one another. To the best of our knowledge, this is the first calibration study that takes this special problem structure into consideration in the calibration procedure.
• One salient feature of the proposed approach is that it adaptively generates data on the fly from the computer model during the calibration procedure so that the most appropriate data can be produced from the computer model. This is different from the traditional frameworks that generate simulation data at pre-selected data points a priori.
• To minimize multiple loss functions, we design a nonlinear optimization algorithm and derive its associated convergence guarantees. Under certain (reasonable) conditions, our analysis provides strong guarantees for our algorithm for different classes of functions.
• We conduct numerical experiments on a wide range of settings and a case study with real data to demonstrate the performance of the proposed approach. The results indicate its superiority with respect to multiple criteria, including estimation accuracy, uncertainty quantification, computational efficiency, and scalability, over alternative approaches.
A preliminary version of this work is presented in a short conference paper (Xu et al. 2021).
In this paper, we substantially extend the analysis to connect the calibration procedure with the statistical parameter estimation, derive and present convergence guarantees, quantify estimation uncertainties, and provide extensive empirical results via numerical examples and a case study.
The rest of this paper is organized as follows. Section 2 discusses the BEM application as our motivating example and other examples from various applications. Section 3 provides a literature review of related methods. Section 4 gives a formulation of the problem. Section 5 develops a new algorithm to solve the problem and studies the convergence properties. We evaluate the performance of our methodology using numerical examples and building energy simulation case study in Sections 6 and 7, respectively. Finally, Section 8 provides a brief conclusion.
Article submitted to INFORMS Journal on Data Science; manuscript no. IJDS-2022-0008
In this section, we introduce the motivating building energy application and discuss the wide applicability of our proposed method.
A physics-based BEM, which simulates the building energy operations, can play an important role in optimizing building design and operations. To effectively use the BEM, several parameters that specify building characteristics need to be estimated using operational data. Table 1 lists the parameters employed in the BEM simulator, called EnergyPlus, which is developed and maintained by the U.S. Department of Energy's National Renewable Energy Laboratories (U.S. Department of Energy 2019).
A unique aspect of the BEM is the seasonal dependency of the parameters, which is the focus of our study. Among the parameters, some parameters related to lighting, domestic hot water, window material, and ventilation can be considered as global parameters that are associated with the entire year-long operational data and simulation. On the contrary, cooling-and heating-related parameters are only associated with their seasonal portions of the year-long data. Interestingly, it is possible that different seasons may overlap in some areas. For example, in Texas, buildings operate heating devices from November to April, whereas cooling systems run from March to November.
Hence, during some periods of time, they operate both heating and cooling devices. Motivated by the BEM parameter calibration, this study develops a new calibration approach when parameters need to be calibrated with different portions of data.
There are various applications where our proposed method can be utilized. First, consider the hydrological computer model that simulates water flow to estimate natural system behaviors affected by climate change, land use, etc. This model can be used to investigate the rainfall-runoff relationship for the windrow compositing pad to estimate the amount of runoff (Duncan et al. 2013, Bhattacharjee et al. 2019). To accurately simulate rainfall-runoff dynamics, several parameters, including the depth of surface, depth of subsurface, saturated hydraulic conductivity of the gravel media, and saturated hydraulic conductivity of the supporting soil below the media, need to be calibrated. Among them, the depth of surface and depth of subsurface can be regarded as global parameters that are associated with the entire year-long data, whereas others are possibly season-dependent parameters associated only with their seasonal portion of data due to the fact that they are generally affected by weather.
Another application is the wind flow model that characterizes spatially heterogeneous wind patterns within a wind farm due to the interactions among turbines (You et al. 2017(You et al. , 2018)).
Recently, Howland et al. (2022) presented a flow model that consists of two sub-models: the poweryaw model and wake effect model. The power-yaw model predicts the power production of upwind turbines (yawed turbines) with a parameter λ p that adjusts the power generation amount in a a yaw-misalignment setting. The wake effect model predicts the wind speed deficit at downstream turbines (waked turbines) due to the shading effect of upwind turbines (Liu et al. 2021) so that it estimates the power production of the waked turbines. This wake effect model takes the wake spreading coefficient k w and the proportionality constant of the modeled Gaussian wake σ 0 as parameters. Among the parameters, λ p can be considered as a global parameter that is employed in all wind conditions. On the other hand, k w and σ 0 depend on atmospheric conditions (Howland et al. 2022), which are calibrated using a portion of data collected under specific conditions.
Lastly, consider a system with a hierarchical block structure. One example could be a pandemic simulation model for infectious diseases. When an outbreak, such as COVID-19, occurs, it is important to predict the spread of infectious diseases. In the simulation models, some parameters such as infection rate, reinfection rate, and positive test rate should be calibrated to achieve accurate predictions and they are usually associated with specific geographical areas. These areas overlap hierarchically from a country level, followed by state, county, city levels, etc. and each parameter should be calibrated using the data obtained from the associated area of interest. For instance, data to predict the infection rate at a city level overlap with that at a county level.
In this section, we provide a literature review of several relevant methodologies.
Article submitted to INFORMS Journal on Data Science; manuscript no. IJDS-2022-0008
The seminal work in modern parameter calibration is Bayesian calibration (Kennedy andO'Hagan 2001, Higdon et al. 2004). It uses a linear linkage model to connect the responses of physical processes with computer model's outputs, and the calibration parameters and hyperparameters in mean vectors and covariance matrices of the GP emulators are estimated from posterior distributions. The Bayesian calibration methodology has gained popularity due to its ability to quantify estimation uncertainties with limited field observations and/or scant simulation data (Kennedy and O'Hagan 2001, Higdon et al. 2004, Gramacy et al. 2015). Recently, Tuo and Wu (2016) suggested that the Bayesian approach could result in unreasonable estimates for imperfect computer models.
Here, an imperfect computer model implies that the outputs of the model, even with the optimally calibrated parameters, are different from the expected outcomes of a physical process. To address this limitation, Tuo and Wu (2015) presented a new calibration approach from a frequentist point of view.
Among these statistical methods, studies in the building energy literature have mainly adopted the Bayesian calibration approach (Chong andMenberg 2018, Coakley et al. 2014). However, the Bayesian approach has several limitations. First, it is computationally intensive, so studies mostly use aggregated data, e.g., monthly (Heo et al. 2012, 2013, 2015, Kim and Park 2016, Kristensen et al. 2017b, Li et al. 2016, Lim and Zhai 2017, Sokol et al. 2017, Tian et al. 2016) or annual (Booth et al. 2013) data. Some studies use hourly data, but they use a small subset of hourly data, instead of using the whole dataset collected over a sufficiently long period of time (Chong et al. 2017, Chong and Lam 2017, Manfren et al. 2013, Menberg et al. 2017, Kristensen et al. 2017a).
In Chong et al. (2017), only 80 hourly data samples are selected from the dataset collected for three months. To relieve the computational burden, a lightweight Bayesian calibration that uses a linear regression emulator is presented in Li et al. (2016). Menberg et al. (2017) employ the Hamiltonian Monte Carlo sampling to obtain the posterior distribution more efficiently. Despite these advances, the Bayesian calibration approach still remains computationally demanding, which limits the use of large-size datasets.
Second, the results from Bayesian calibration are sensitive to the prior specification (Liu et al. 2021). Our analysis, which will be discussed in more detail in Section 7, indicates that the posterior density is heavily affected by the prior specification. When a non-informative prior is employed, the resulting posterior tends to become non-informative, showing a relatively flat posterior distribution.
We believe this is because the use of limited data (either aggregated monthly or yearly data, or a small number of selected hourly data) provides insufficient information to generate a meaningful posterior. While a more informative prior leads to a sharper posterior, it requires domain knowledge for the appropriate prior setting.
Further, the Bayesian approach in the BEM calibration studies limits the number of calibration parameters to 2 to 5. This is because the hyperparameter estimation of the covariance function becomes computationally prohibitive as the problem size grows. Further, compounded with the aforementioned issues of limited data, when multiple parameters are calibrated, overparameterization may occur, resulting in a relatively flat posterior, which does not generate meaningful calibration results.
Most importantly, existing approaches cannot handle multiple blocks of overlapping data. In particular, when the block sizes are different, one might suggest that the problem is similar to the class imbalance problem where the data size of each class is substantially different (Byon et al. 2010). In such cases, one possible solution could be to impose different weights on each data block.
However, how to differentiate weights across blocks is not straightforward. In particular, when data blocks overlap, if the algorithm calibrates the parameters with their small portion of data with a higher weight, it would also affect other blocks with larger portions of data. Thus, the weight should consider the overlapping portion as well as the data size of each block; devising a proper weighting scheme with rigorous justification is challenging.
We also review two other streams of work that are relevant to this study, namely, adaptive sampling and multi-task learning. First, adaptive sampling, also known as sequential design in the statistical literature and active learning in the machine learning literature, has been actively studied. The idea is to collect informative data points sequentially in order to construct surrogates and/or minimize a black box function by accounting for the trade-off between exploration and exploitation (Liu et al. 2018). First, in surrogate modeling studies that often use GPs, adaptive sampling strategies choose the next design point that maximizes predictive variance or maximizes the average reduction in variance (Gramacy 2020, Liu et al. 2018). These studies aim to construct accurate surrogates.
Second, Bayesian optimization (BO) has received much attention recently for minimizing expensive black box function. BO proceeds by modeling the objective function via a GP, optimizing an acquisition function to find the next design point and then updating the posterior distribution of the GP. Acquisition functions, e.g., expected improvement (Jones et al. 1998) and lower confidence bound (Srinivas et al. 2010), attempt to strike a balance between the exploration of a new design point with high uncertainty and exploitation with a low objective function value. More details on BO can be found in Shahriari et al. (2015) and Frazier (2018). Typically, BO aims to minimize a single objective. Although the multi-task BO (Swersky et al. 2013) is proposed to handle multiple objective functions, it uses the same training data for all objectives to perform optimization.
Our proposed method solves a multi-objective problem with multiple blocks of data that possibly overlap with one another.
Article submitted to INFORMS Journal on Data Science; manuscript no. IJDS-2022-0008
On the other hand, multi-task learning aims to learn multiple tasks simultaneously by leveraging knowledge learned from one or some tasks while exploiting commonalities and differences (Zhang and Yang 2021). It helps to improve prediction accuracy and computational efficiency for learning multiple tasks, compared to individual task learning. In principle, the proposed multi-block calibration methodology can be regarded as a special case of multi-task learning, if we consider calibrating each block of parameters as each task; updated parameters associated with one block affect the loss function in other blocks so that loss function values benefit from or are hurt by others. However, multi-task learning is typically performed in a supervised learning setting with a fixed dataset, whereas our proposed method exploits adaptive sampling. Moreover, in multi-task learning, each task is learned with its own dataset which is disjoint from other tasks' data, however, our approach accounts for overlapping blocks of data.
We first present a calibration problem in a general single-block calibration setting in Section 4.1 and extend it to a multi-block parameter calibration problem in Section 4.2.
Let x ∈ Ω ⊆ R d denote the vector of physically observable input variables of dimension d ∈ Z + in a system, where Ω is a convex and compact region of input variables that is a subset of R d . Let θ ∈ Θ ⊆ R p denote a set of calibration parameters of dimension p ∈ Z + . Let y(x) ⊆ R denote the noisy field observational data of its true physical process ζ(x) at input x. Let η(x, θ) ⊆ R denote the response from the computer model, simulating the true process ζ(x) at input x. For instance, in the building energy simulation, x can be hourly environmental condition such as temperature, y(x) the actual hourly electricity consumption, and η(x, θ) the electricity consumption generated from BEM. The unknown calibration parameter vector θ needs to be estimated using y(x) and η(x, θ).
Assuming the simulator represents the physical process accurately, Higdon et al. (2004) present the regression model to connect physical process responses with computer model outputs as
where n denotes the number of observations and ϵ j an observation error that follows an independently and identically distributed (iid) normal distribution, ϵ j iid ∼ N (0, σ 2 ).
The likelihood function of θ is then
or more explicitly,
(3)
To estimate θ, let us consider the maximum likelihood estimation (MLE) method. The loglikelihood function is
We now maximize the log-likelihood to obtain θMLE given σ 2 as θMLE = arg max
Note that maximizing the likelihood function is equivalent to minimizing the empirical loss function that quantifies the discrepancy between the computer model output η(x, θ) and physical process data y(x) with L 2 norm. This calibration procedure that minimizes the empirical loss aligns with the idea of L 2 calibration studied in Tuo and Wu (2015). Tuo and Wu (2015) aimed to find the parameters that minimize the L 2 distance between the physical process responses and computer model outputs. Under the limited data setting, they built emulators for both computer model outputs and physical process responses. However, we consider a situation where a sufficiently large number of physical observational data is available and a computer model is relatively cheap to run.
With data availability, instead of generating data from a computer model at pre-defined design points, we want to generate computer model data more effectively so more useful data can be used for the calibration purpose and directly use it to estimate the parameters without building emulators. To this end, we minimize the loss functions in (5) with respect to θ. As the computer model output η(x, θ) is likely nonlinear, we consider nonlinear optimization techniques. Consider a first-order optimization method, i.e., gradient descent (GD). The GD method works by iteratively updating θ with θ ← θ -α∇F (θ) until a suitable termination criterion is satisfied, where F (θ) is the empirical loss function at θ, ∇F (θ) represents its gradient, and α is a step size. Here, because η(x, θ), as the output from a black box computer model, does not have a closed-form expression, ∇F (θ) can be approximated using the central finite difference
Article submitted to INFORMS Journal on Data Science; manuscript no. IJDS-2022-0008
where h > 0 is a small disturbance to the i th component of θ, e i denotes a p × 1 vector with i th element being one and others zero.
One key advantage of using nonlinear optimization is that it can adaptively generate computer model data η(x, θ). In other words, as θ is updated through the iterative GD process, we can run simulations to get η(x, θ) at the newly refined θ. It allows us to get more informative data toward minimizing the empirical loss (or maximizing the likelihood).
The single-block parameter calibration presented in Section 4.1 treats all parameters equally using the entire dataset. We propose a new block-wise calibration approach when some parameters are associated with a specific subset of the entire dataset.
Recall that in the BEM application, we have three distinct groups of parameters, denoted by θ g , θ c , and θ h in Table 1. The parameter set θ g consists of global parameters so it should be globally With this problem structure, we calibrate θ b using its corresponding portion of dataset, i.e., the b th data block. By using the right data block, it is expected to obtain more accurate estimates, compared to the single-block calibration that does not differentiate the block-dependent heterogeneity.
To account for the block dependency, we consider multiple likelihood functions, each associated with the b th block of datasets, as follows:
where
To maximize the log-likelihood, we consider the corresponding empirical loss function as
and obtain the estimates for the b th block of parameters θ b by
We refer to the proposed block-wise calibration approach as multi-block calibration (shortly, M-BC)
in the subsequent discussion.
It should be noted that the simulation output η(x j , θ b ; θ -b ) could possibly depend on θ -b , as well as θ b , when the b th block of data overlaps with other blocks associated with some parameters in θ -b .
For instance, let us consider the cooling-season parameter θ c . The simulation output during the cooling season (March-November) also depends on the setting of global parameters θ g . Also, during
March to April when the cooling and heating seasons overlap, heating-season parameters also affect the energy consumption output during a part of the cooling season. Therefore, in solving (10), the resulting θ * b may also depend on θ -b . Our ultimate goal is to find θ * b that minimizes F b , while at the same time, other block parameters minimize their own loss functions. In the next section, we provide an iterative algorithm and its convergence properties under certain conditions.
In this section, we propose a new algorithm to calibrate parameters in a block-wise manner. We also provide convergence guarantees for the algorithm under certain conditions. Finally, given convergence, we employ the MLE's asymptotic properties and construct confidence intervals to quantify estimation uncertainties.
Article submitted to INFORMS Journal on Data Science; manuscript no. IJDS-2022-0008
Recall that existing studies typically pre-generate data and build emulators for η(•) using the generated data. This approach is not directly applicable to the multi-block parameter calibration problem. One may attempt to build B emulators, say, η1 (x, θ 1 ; θ -1 ), . . . , ηB (x, θ B ; θ -B ), using each To this end, we design a new algorithm to solve the block-wise calibration problem. The fundamental idea is to devise a block-wise optimization algorithm so that the parameters in different blocks can be calibrated with their own objective functions over iterations. Specifically, we propose cyclically optimizing one block of parameters each time while keeping the parameters in other blocks fixed at their most up-to-date values. For instance, in the BEM calibration, we have θ 1 = θ g , θ 2 = θ c , and θ 3 = θ h with B = 3. We iteratively optimize θ g that minimizes F g (θ g ; θ -g ), while other parameters in θ c and θ h are fixed at their current iterates. Then we optimize θ c for F c (θ c ; θ -c )
with the previously optimized θ g and continue with a similar procedure for F h (θ h ; θ -h ) with the previously optimized θ g and θ c . This procedure is repeated until pre-specified stopping criteria are satisfied.
Let θ k b denote the iterate of θ b after finishing the k th iteration for b = 1, 2, . . . , B. At the (k + 1) th iteration, let F b (θ b ; θ k+1 -b ) in ( 9) denote the loss function for the b th block of parameters θ b , given the remaining parameter iterates θ k+1 -b . Without loss of generality, we optimize each block in an ascending order, starting from the first block, i.e., b = 1. Therefore, in optimizing θ b that minimizes
Here, note that the parameters before the b th block are previously updated, whereas those after the b th block are not refined yet. In this scheme, the loss function for the b th block in (9) at the (k + 1) th iteration can be written as
The computer model output η(•) is likely nonlinear, as most computer models consist of a set of complex mathematical functions. Hence, we consider nonlinear optimization methods to find θ b that minimizes F b (θ b ; θ k+1 -b ), given θ k+1 -b . Among various optimization methods, we use the most commonly used method, which is GD. Our approach can be easily extended to other second-order optimization methods such as Newton's method or quasi-Newton methods.
In minimizing F b (θ b ; θ k+1 -b ) over θ b , we further consider the inner iteration, given θ k+1 -b . Here, the inner iteration is for the GD-based parameter updates in each block, and the outer iteration implies the cycle for all blocks. Let θ k,m b denote the iterate of θ b at the m th inner iteration of the outer (k + 1) th iteration. We update θ k,m b by
where α k,m b > 0 is the step size and the
where θ b,i is the i th parameter in θ b for i = 1, 2, . . . , p b , with p b being the number of parameters in θ b for b = 1, 2, . . . , B. Since the loss functions have no mathematical closed-form expressions due to the black box nature of the simulator, similar to (6), we use a finite-differencing to approximate the gradient for each loss function as follows:
where h > 0 is a small perturbation to the i th parameter of each θ b and e b,i is a p b × 1 vector with its i th element being one and others zero.
We repeat the inner iteration for each block until some termination criteria are met. , we move to the next (b + 1) th block to minimize its loss function
After we finish the cycle with all blocks, we set k = k + 1 and repeat the procedure. We finish the iteration when either the norm of gradient of each loss function is less than a small tolerance τ 1 , the relative difference of loss function values is less than a small tolerance τ 2 , or the maximum number of iteration is larger than some value τ 3 . Algorithm 1 summarizes the proposed M-BC procedure.
We would like to highlight that the proposed cyclic calibration approach is different from the block coordinate descent (BCD) algorithm. The major difference is that our approach handles multiple objective functions, each of which is associated with its corresponding subset of the dataset.
On the other hand, BCD considers a single objective function.
As discussed earlier, an important feature of our approach is that it adaptively generates computer model outputs during the calibration procedure, unlike existing approaches that use precollected samples. As the optimization proceeds, newer parameters that make the computer model while convergence criterion not met do 7:
Run simulations for
) using ( 13)-( 14) and set a step size α k,m b 10:
end while As a remark, in multi-block calibration one can randomly choose the next block, but in this section we consider the fixed ordering that proceeds in a cyclical way. In Section 6, we empirically compare the cyclic ordering with stochastic ordering. Additionally, we note that, depending on the problem structure, the order of blocks may matter. However, in our numerical examples in a wide range of settings, we obtain comparable results in different orderings.
This section presents the convergence properties when each loss function F b , b = 1, 2, . . . , B, is a general nonconvex, strongly convex, or convex function. While we utilize the standard analysis procedure of GD, the block-wise treatment with multiple loss functions poses substantial challenges in analyzing the convergence properties of the proposed M-BC approach. In each of the three settings, we impose conditions on F b and possibly the iterates generated by Algorithm 1 to achieve certain convergence properties. In the special case where the number of blocks is equal to 1, our results recover the standard results for GD. The detailed proofs of all lemma and theorems are available in the online supplementary document.
Throughout the convergence analysis, we use the notations summarized in Table 2. We understand that the notations are complex, mainly due to the double (inner and outer) iteration nature in M-BC. For ease of understanding, Figure 2 In our analysis, we make the following assumption.
Article submitted to INFORMS Journal on Data Science; manuscript no. IJDS-2022-0008
Table 2 Nomenclature for convergence analysis
The iteration number for the outer iteration. m
The iteration number for the inner iteration, i.e.,
The number of inner iterations for the b th block at the k th outer iteration.
Fixed number of M b (k) regardless of the k th outer iteration.
Iterate
The iterate after all blocks are updated at the k th outer iteration, i.e.,
The unique global minimizer of F b in the strongly convex case; a point at which the minimum of F b is attained in the convex case, i.e.,
The iterate of the b th block parameters after finishing the inner iteration of the b th block at the k th outer iteration.
The unique global minimizer of F b given θ * -b in the strongly convex case; a point at which the minimum of
The iterate of the b th block parameter at the m th inner iteration of the (k + 1) th outer iteration, given
The unique global minimizer of F b given θ k+1 -b in the strongly convex case; a point at which the minimum of
The result above is analogous to that of GD. Namely, the progress made over a sequence of iterations (for M-BC, this is over the inner iterations) is proportional to the sum of the norms of the gradient squared. Using Lemma 1, we first analyze the convergence properties of M-BC on the general nonconvex case in Theorem 1.
Theorem 1 shows that under certain conditions, the iterate θ k converges to a stationary point θ * as k → ∞, i.e., the norm of the gradient of all objective functions vanishes in the limit.
Theorem 1 (General nonconvex case). Suppose Assumption 1 holds. Further, suppose that
We make a few remarks about the conditions in ( 16). From Lemma 1, it follows that
). Thus, the inequalities ( 16) imply that F b (θ k ) ≥ F b (θ k+1 ), i.e., the function value for the b th block decreases after all blocks have updated their associated parameters (across outer iterations). This seems to be a strong assumption in that the blocks do not compete over the course of optimization. However, this does not require that when every other block updates, the function value F b of the b th block is non-increasing. Instead, F b is non-increasing across outer iterations.
The result in Theorem 1 has an important implication. When we eventually find the stationary point θ * b for the parameters in the b th block, other blocks' parameters also converge to a stationary point θ * -b as well and
It implies that even though we calibrate parameters in a block-wise manner, our approach obtains a concurrent convergence of the whole parameter vector, that is, We make the possible extension of conditions in (16) at the end of this section. Additionally, we relax the conditions in ( 16) and prove the similar result in Corollary 1 by instead assuming the weaker technical condition in (18).
Corollary 1. Suppose Assumption 1 holds. Further, suppose that the infinitely cumulative sum of the decreases of the function values during the inner iterations is finite, i.e.,
Then, for all b,
Corollary 1 states that as long as the sum of loss function decreases is finite, M-BC converges to a stationary point. It allows other block calibrations to increase the b th block loss, so the noncompeting condition in ( 16) is not needed. Thus, this result is more general than that in Theorem 1.
Article submitted to INFORMS Journal on Data Science; manuscript no. IJDS-2022-0008
Next, we show the convergence rates for the strongly convex (Theorem 2) and convex (Theorem 3) cases.
Theorem 2 (Strongly convex case). Suppose Assumption 1 holds. Further, assume that the
for all k ≥ 0 and b, where
Then, for all k and b,
Both conditions in (20) bound the optimality gap between the progress made in inner iterations versus the progress made in outer iterations. Since the objective function of each block is changing across outer iterations, the absolute value safeguards against the possibility of a current iterate having lower function value than the optimal solution. We note that this is possible in our setting since the optimal solution of θ * = (θ * b ; θ * -b ) is defined as the point where all objectives are optimized and thus can potentially have a higher function value than the point (θ k, * b ; θ k -b ) has. Note that the conditions for Theorem 2 do not require that the function of the b th block is monotonically non-increasing across outer iterations. That being said, these assumptions are relatively weak and can be tightened so that the result still shows the linear convergence rate at possibly a slower rate. We provide the corollary about this.
Corollary 2. Suppose Assumption 1 holds. Further, assume that the function F b is µ b -strongly convex over θ b for any fixed θ -b for b = 1, 2, . . . , B. Moreover, suppose that
for all k ≥ 0 and b, where
for all k and b. Let α b ∈ (0, 1/L b ] for all b. Then, for all k and b,
Theorem 2 and Corollary 2 show linear convergence of the iterates generated by the M-BC algorithm, analogous to that of GD (in the special case where b = 1, we recover the result of GD). While convergence is derived across outer iterations, the rate of convergence depends on the number of inner iterations. Note that ρ = c 1 c 2 (1 -α b µ b ) M b < 1, thus, the larger the number of inner iterations, the faster the rate (constant), or alternatively, the larger the errors that are acceptable in ( 22) (larger constants c 1 and c 2 ). Note that compared to Theorem 2, Corollary 2 has somewhat stronger result that shows outer function value F b (θ k ) keep non-increasing with the linear rate
Additionally, Theorem 3 shows the convergence rate for the convex case with certain assumptions.
Theorem 3 (Convex case). Suppose Assumption 1 holds. Further, assume that the function
for all k ≥ 0 and b, where
Then, for all k ≥ 1 and b,
The first inequality in (25) is the same as the first condition in the strongly convex setting (see ( 20)). The second bounds the distance between the difference from the current iterate's function value to the optimal solution's function value across two outer iterations.
Note that the first assumption in ( 25) is relatively weak and can be tightened so that the result still shows the sublinear convergence rate at possibly a slower rate. We provide the result about this in Corollary 3.
Corollary 3. Suppose Assumption 1 holds. Further, assume that the function F b is convex for b = 1, 2, . . . , B over θ b for any fixed θ -b . Moreover, assume that
for all k ≥ 0 and b, where c 1 > 0. Additionally, assume that
for all k and b. Let α b ∈ (0, 1/L b ] for all b. Then, for all k and b,
Article submitted to INFORMS Journal on Data Science; manuscript no. IJDS-2022-0008
Theorem 3 and Corollary 3 show that in the convex setting, the iterates generated by the M-BC algorithm converge to a minimizer at a sublinear rate that depends on both the number of inner and outer iterations. Again, this is analogous to GD, and in the special case with b = 1, we recover the result of GD. Note that Corollary 3 has somewhat stronger result that shows outer iterates θ k keep non-increasing with the sublinear rate O 1/(kM b ) .
We make a comment about the number of iterations required for our algorithm to converge (i.e., number of inner iterations for a given block and number of outer iterations) and provide the computational complexity of M-BC. Of course, the number of iterations highly depends on the convergence criterion used in Algorithm 1. To give a concrete example, consider the case of strongly convex functions, where the convergence criterion in the inner loop is based on the optimality gap and suppose the inner and outer termination tolerances for relative differences are set as 0 < ϵ outer ≤ ϵ inner < 1. In this case, one can show that the number of inner iterations required for each block is O log (1/ϵ inner ) , and the total number of iterations is O log ϵ inner (1/ϵ outer ) . Thus, combining the time for data generation and for parameter estimation, the computational complexity becomes
O n log ϵ inner (1/ϵ outer ) , assuming that the simulation running time is proportional to the data size n. Similar results can be shown for the convex and nonconvex cases under appropriate termination conditions.
Finally, we would like to make additional remarks about the conditions in our analysis. While the conditions we impose in deriving convergence properties, for example, that the objective functions do not compete for all θ b , may seem to be strong, we make these assumptions for the analysis of the M-BC algorithm to be analogous to that of GD. Furthermore, the conditions in our analysis are less restrictive than that as they are over the course of the iterates generated by the algorithm, and not for the entire space. These conditions may not be verified except in relatively simple settings, as they require information about the optimal solution and/or problem specific parameters. In practice, after running the method on a specific problem, one can observe the trajectories to see if the conditions are satisfied. In our numerical examples in a wide range of settings and the BEM case study, we did not observe the situation that the non-competing conditions were severely violated.
That being said, we conclude this section by discussing possible extensions of the proposed approach and convergence results. Suppose that in the general nonconvex case, instead of assuming ( 16), we can relax the conditions as
for all k ≥ 0 and b = 1, 2, . . . , B with δ k b,1 > 0 and δ k b,2 > 0. That is, we allow the objective function of the block to increase by some levels (δ k b,1 and δ k b,2 ) after calibrating parameters in other blocks.
Under this assumption, if δ k b,1 = δ 1 and δ k b,1 = δ 2 for all b = 1, 2, . . . , B, then we conjecture that one can only show convergence to a neighborhood of the solution that depends on δ 1 and δ 2 . On the other hand, if δ k b,1 and δ k b,1 are summable, then a similar result to Theorem 1 can be proven. Similarly, the conditions for the results in the strongly convex and convex cases can be relaxed, and under appropriate conditions, convergence (and rates) can be proven. Our BEM case study is Section 7 demonstrates that the non-competing condition in ( 16) is satisfied. However, we plan to further investigate competing cases with (30) in our future study, as it may arise in other applications.
When θ b converges to a minimizer θ * b for each b based on the results in the previous section, we obtain the MLE for the calibration parameters. Thus, we can entertain its appealing theoretical properties such as consistency and normality. From the likelihood function in (8), the consistency property says θb,MLE converges in probability to a true parameter vector, denoted by θ b,0 , as n b → ∞, i.e., θb,MLE p -→ θ b,0 for each b. Further, the asymptotic normality indicates that the estimator
), where I(θ b,0 ) is an expected Fisher information matrix. We use the MLE's asymptotic properties to construct the confidence intervals of calibration parameters for uncertainty quantification.
Based on (8), the expected Fisher information matrix is represented by
where lb θ b |y(x), θ -b is a log-likelihood function from a single observation for the b th block as
The expected Fisher information matrix can be approximated by its empirical counterpart,
Using the chain rule, we obtain the (i, i ′ )-entry of ∂ 2 lb θ b |y(x j ), θ -b /∂θ 2 b as follows.
Article submitted to INFORMS Journal on Data Science; manuscript no. IJDS-2022-0008
Here, since the first-and second-order partial derivatives of η(x j , θ b ; θ -b ) are not analytically available, we obtain them numerically using the central finite difference (Abramowitz and Stegun 1972) as follows:
for i, i ′ = 1, 2, . . . , p b , where h, h ′ > 0 are small perturbation values. Further, assuming the physical process data y(x) have the same variance σ 2 across multiple blocks, we can find its MLE as
Then we can obtain the asymptotic 100(1 -α)% Wald confidence interval (CI) for each component of θb,MLE as follows:
for b = 1, 2, . . . , B, z 1-α/2 is a critical value of the normal distribution, and I -1
ii denotes the i th diagonal entry of the inverse of the Fisher information matrix I.
We provide guidance on defining the algorithmic parameters in M-BC (Algorithm 1).
• α k,m b > 0 (step size of the b th block at the (k + 1) th outer iteration and m th inner iteration):
The step size α k,m b can be obtained by a backtracking line search, as described in Algorithm 2 below. In the backtracking line search, the parameters can be set as ᾱ = 1, ρ = 1/2, c = 10 -4
(or 1/2) (Boyd andVandenberghe 2004, Nocedal andWright 2006). Note that, in principle, the step size α k,m b could also be set to a sufficiently small constant or a diminishing sequence (e.g., 1/m).
Algorithm 2 Backtracking Line Search
Update α ← ρα 4: end while
• h > 0 (bandwidth of finite-differencing when approximating the gradient): The bandwidth h is set to 10 -8 , which is the square root of machine precision. We note that in the setting where exact function values can be computed, the choice of bandwidth is optimal up to constants (Moré and Wild 2012). In practice, we note that a sufficiently small value h works well in most cases (e.g., 10 -4 ).
• τ 1 > 0 (termination tolerance: norm of the gradient): We suggest that τ 1 would be set as any sufficiently small value (e.g., 10 -4 ).
• τ 2 > 0 (termination tolerance: relative difference of loss function values): We suggest that τ 2 would be also set as any sufficiently small value (e.g., 10 -4 ).
• τ 3 > 0 (termination tolerance: maximum iterations): The maximum number of iterations is used to safeguard the algorithm from running forever. If the algorithm reaches this value, it exits the loop. We suggest that τ 3 would be set as sufficiently large value as long as computing resources are available. In our implementation, we set 2000 for the outer loop and 1000 for the inner loops of the algorithm.
In this section, we evaluate the calibration accuracy of our proposed method in comparison with other alternatives through numerical studies. We also report the half length of CI and coverage rate of each parameter to quantify estimation uncertainties. All experiments, except those for alternative methods in Section 6.4, are conducted with MATLAB 2014a. For the alternative methods in Section 6.4, R is used. All experiments use the following computing environment: 64-bit Windows OS with the Intel Xeon CPU E5-2697 @ 2.60 GHz processor and 128 GB RAM.
We consider three numerical examples with different settings of the following physical processes and computer models.
(a) Example I (perfect computer model, block-wise convex loss functions)
• Physical process
(b) Example II (perfect computer model, nonconvex loss functions)
• Physical process
(c) Example III (imperfect computer model, nonconvex loss functions)
• Physical process
In all examples, x ∼ U (0, 2π) and ε ∼ N (0, 0.01). In Examples I and II, the terms can be divided into three input domains x < 4π/5, 2π/5 ≤ x < 8π/5, and x ≥ 6π/5, respectively, to represent different blocks (e.g., seasonality in BEM). Hence, θ 1 , θ 2 , and θ 3 are the block parameters that need to be calibrated with dataset associated with the first, second, and third blocks of data collected when
x < 4π/5, 2π/5 ≤ x < 8π/5, and x ≥ 6π/5, respectively. Note that these data blocks are not disjoint.
For example, θ 1 and θ 2 together affect the computer model response for 2π/5 ≤ x < 4π/5. Thus, one should not calibrate each θ i separately. In Example III, the first term is globally employed regardless of the input value. The second and third terms are employed when x < π and x ≥ π, respectively, which mimic the block dependency. Thus, θ 1 should be calibrated with the entire dataset, whereas the block parameters θ 2 and θ 3 need to be calibrated with the second and third blocks of data collected when x < π and x ≥ π, respectively. Figure 3 In implementing the proposed approach, we choose the step size (α k,m b in ( 12)) using the backtracking line search in Algorithm 2. The termination condition is set as the norm of the gradient of the loss functions is less than 10 -4 in both inner and outer iterations.
We present the implementation results in comparison with the following two alternatives.
(a) Hybrid-block Calibration (H-BC): In this alternative, the computer model correctly simulates its response according to underlying block properties. Specifically, similar to M-BC, it considers the block dependency, using the corresponding subset of data for each block, rather than
Article submitted to INFORMS Journal on Data Science; manuscript no. IJDS-2022-0008 using the entire data. However, unlike M-BC that uses multiple loss functions, H-BC uses a single loss function in (5) and optimizes θ = (θ 1 , θ 2 , θ 3 ) altogether, namely,
In solving (40), uses the correct block of data to compute the gradient of the loss function, but instead of exploiting a block-wise minimization, it updates the parameters all at once, i.e., it updates θ 1 , θ 2 , θ 3 altogether at each GD iteration. only, but it also employ the other term with I(x ≥ π). Then it uses the single loss function, shown in ( 5), to calibrate parameters. We include S-BC because in the BEM simulation, seasonality can be ignored and one may blindly run the BEM simulator without correct scheduling.
When loss functions are convex, we can start from an arbitrarily selected starting point to attain a global minimum. In practice, however, we do not know the form of the loss functions, since the computer model is a black box. Thus, we use multiple initial points and select the best estimates that produce the smallest mean squared error (MSE) in a training set consisting of 1000 data samples. We then evaluate the performance using a test set with another 1000 data samples.
Table 3 summarizes the average of calibrated parameters (and standard deviation) in the third to fifth columns, obtained from 200 experiments, as well as MSEs in the test set in the last four columns. Each F i represents the MSE at the i th data block, and F all in the last column is overall MSE with all test data samples. Clearly, each calibrated value for θ 1 , θ 2 and θ 3 using M-BC in Examples I and II is closer to the true parameter values θ = (θ 1 , θ 2 , θ 3 ) = (-1, 10, 5).
We summarize some observations from Table 3.
• Since M-BC optimizes each block of parameters with the correct subset of data by minimizing the corresponding loss functions, it provides the most accurate calibration results with the lowest MSEs.
• H-BC also uses the correct portion of data, so its MSEs tend to be smaller than those from S-BC. But, its calibrated values deviate from the true values in some cases, e.g., θ 2 and θ 3 in Example I. It is because its single loss function for all three groups of parameters is quite complicated and highly nonconvex with interations between parameters (e.g., (θ 1 +1)(θ 2 -10)).
On the contrary, in the proposed M-BC, each loss function, given fixed parameters in other blocks, is simpler than the H-BC's loss function. For instance, in Example I, even though η(•) is nonconvex over (θ 1 , θ 2 , θ 3 ), the loss function for each parameter becomes block-wise convex.
This local convexification helps the M-BC procedure find the solution more effectively.
• The different performance of H-BC in Examples I and II is related to the problem structure where different weights are assigned to each term. In Example I, for the term including I(x < 2π/5) in the computer model, we assigned 1/2, while we assigned 1/10 and 1/20 as weights, respectively, for the terms including I(2π/5 ≤ x < 8π/5) and I(x ≥ 8π/5). Thus, θ 1 is the most influential parameter and the loss function changes rapidly as θ 1 varies. With this problem structure, H-BC weighs θ 1 more, leading to the correct calibration of θ 1 relative to θ 2 and θ 3 .
In Example II, the advantage of M-BC, compared to that of H-BC, in terms of calibration accuracy, is less clear. We believe it is because the weights for each term are not substantially different. These results indicate that the performance of H-BC is sensitive to the problem structure. All things considered, M-BC achieves the best accuracy in all examples.
• Obviously, some of the parameter estimates from S-BC largely deviate from their true values in many cases, leading to larger loss values. Recall that S-BC ignores the block property (by ignoring I(•) terms in the computer models), while observational data from the physical process are collected in different blocks. This mismatch between the physical process responses and simulation outputs leads to poor calibration results with large MSEs.
M-BC also has a computational advantage over the two benchmarks, H-BC and S-BC, in general. This is mainly because M-BC calibrates each block of parameters using its associated data only, unlike H-BC and S-BC that update all parameters simultaneously. computationally more efficient than H-BC and S-BC. Further, its running time is similar in all three examples, whereas the computation time of H-BC and S-BC varies significantly depending on the problem structure. We additionally note that comparing computation performance in terms of number of iterations between the three methods is far from trivial. This is because M-BC has the double-loop structure with inner and outer iterations, whereas H-BC and S-BC does not have inner iterations. Moreover, another difficulty is related to the way the data is utilized by the methods. Specifically, as discussed in Section 5.2, assuming that the simulation running time is proportional to the data size n, computational complexity for M-BC can be expressed as O n ϵ inner (1/ϵ outer ) in the strongly convex case. The analogue for H-BC and S-BC (assuming a tolerance of ϵ outer ) can be expressed as O n log(1/ϵ outer ) . While the terms in order notation look similar, hidden are condition number dependencies. For M-BC the dependence is on b κ b where κ b is the condition number for each block problem, whereas for H-BC and S-BC the dependence is on κ, the condition number of the full problem. Thus, their computational complexity cannot be readily and easily compared. do not exist, we do not report the coverage rates. Overall, the results suggest that M-BC provides its empirical coverage rates close to the nominal rate of 95%. On the contrary, the coverage rates of H-BC deviate from the nominal rates in some cases, e.g., those for θ 2 or θ 3 in Example I. The low coverage rates in S-BC are related to inaccurate estimations shown in Table 3.
We compare the performance of M-BC with two common calibration approaches in the literature, including L 2 calibration (Tuo and Wu 2015) and Bayesian calibration (Kennedy andO'Hagan 2001, Higdon et al. 2004). We briefly explain the two benchmark approaches. More details are available in Tuo and Wu (2015), Kennedy andO'Hagan (2001), andHigdon et al. (2004). Given physical observations {y(x j ), x j } n j=1 , L 2 calibration first obtains the estimated true physical process response ζ, where y(x) = ζ(x) + ϵ with an observation error ϵ, using kernel ridge regression (Wahba 1990) in the reproducing kernel Hilbert space. Similarly, using computer model responses {y(x j ′ ), (x j ′ , θ j ′ )} n ′ j ′ =1 generated at pre-specified design points, it builds an emulator η(•, •), often constructed by a GP (Santner et al. 2018), for the computer model. Then, θ is calibrated by solving the following optimization problem:
Next, the Bayesian calibration approach formulates the physical observation using the linear calibration and Bayesian calibration is excessive and they do not produce results in 5 days for each run. In summary, M-BC elicits the superiority of calibration accuracy as well as computational efficiency and scalability, compared to the two existing methods.
We make comments about the difficulty in theoretically comparing computational complexity between the M-BC and Bayesian approach. First, their data generation mechanisms are different.
The Bayesian approach generates a fixed set of data from pre-specified design points a priori, whereas M-BC collects computer response data on the fly. Next, in terms of computation after data is generated, in Bayesian calibration, the main computational issue is the numerical integration with respect to the posterior distribution of calibration parameters θ and matrix inversions for each θ value when using GPs in this numerical integration. If the number of observational data and computer experiments are small, the matrix inversion can be done with O (n + n ′ ) 3 in an iteration where n denotes the number of observational data and n ′ the number of computer experiments (Kennedy and O'Hagan 2001). However, when the data size gets larger, the matrix inversion becomes challenging. Further, it uses MCMC methods to explore posterior distributions. It is known that the MCMC computation generally depends on the number of parameters θ, proposal distribution, and number of simulation runs (or iterations). Chains' mixing time to its stationary distribution also affect the total computational complexity of the Bayesian approach. On the contrary, the computational complexity for M-BC is O n log ϵ inner (1/ϵ outer ) in the strongly convex case as discussed in Section 5.2. With these reasons, the direct comparison of the computational complexity between the M-BC and Bayesian approach does not seem to be adequate. Further, the termination conditions of the two approaches are different. However, their empirical performance shows the computational efficiency and scalability of M-BC over the Bayesian approach.
We consider some variants of the proposed calibration approach. The first variant removes the inner iteration in M-BC. We refer to this variant as M-BC-V. Even with no inner iterations, it still updates parameters in a manner, but unlike M-BC, it does not perform the block-wise minimization. It may seem to be similar to H-BC. However, the main difference is that H-BC updates all parameters simultaneously at once, whereas M-BC-V sequentially updates each parameter block based on the most up-to-date information of other blocks. We also consider another variant using a random block selection. This algorithm is the same as M-BC-V but the block is selected randomly, instead of sequentially (or cyclically). This algorithm is referred to as M-BC-VS with 'S' standing for 'stochastic'.
Table 7 summarizes calibration results from these variants. In Examples I and II, M-BC-V and M-BC-VS provide results comparable to the original M-BC. However, they yield inaccurate results in terms of MSE and unstable estimates with higher standard deviations in Example III. It appears that the block-wise minimization (i.e., the inner iteration for each block) in M-BC brings appealing improvements in some cases and provides consistently better results.
.00 (0.01) 9.99 (0.02) 5.00 (0.03) 0.01 (0.00) 0.01 (0.00) 0.01 (0.00) 0.01 (0.00) M-BC-V -1.00 (0.02) 9.99 (0.02) 5.00 (0.03) 0.01 (0.00) 0.01 (0.00) 0.01 (0.00) 0.01 (0.00) M-BC-VS -1.00 (0.02) 10.00 (0.03) 5.00 (0.03) 0.01 (0.00) 0.01 (0.00) 0.01 (0.00) 0.01 (0.00) II M-BC -1.00 (0.00) 10.00 (0.01) 5.00 (0.01) 0.01 (0.00) 0.01 (0.00) 0.01 (0.00) 0.01 (0.00) M-BC-V -1.00 (0.00) 10.00 (0.01) 5.00 (0.01) 0.01 (0.00) 0.01 (0.00) 0.01 (0.00) 0.01 (0.00) M-BC-VS -1.00 (0.00) 10.00 (0.02) 5.00 (0.02) 0.01 (0.00) 0.01 (0.00) 0.01 (0.00) 0.01 (0.00) III M-BC 0.38 (0.00) -5.75 (0.44) 8.18 (0.05) 6.89 (0.67) 11.35 (1.25) 2.47 (0.16) 6.89 (0.67) M-BC-V 0.36 (0.15) 0.12 (6.76) -0.99 (7.05) 15. 03 (12.43) 23.33 (16.15) 6.80 (18.80) 15.03 (12.43) M-BC-VS 0.36 (0.14) -0.92 (6.30) 1. 83 (6.89) 13.15 (10.15) 21.31 (13.98) 5.07 (11.96) 13.15 (10.15) Further, we investigate whether the block ordering affects the estimation performance in M-BC. With three blocks, there are six pathways in the cyclic block ordering. We do not notice Article submitted to INFORMS Journal on Data Science; manuscript no. IJDS-2022-0008 any significant differences in calibration results among the six pathways in all three examples. In the future, we plan to explore other variants and study their theoretical properties and empirical performance.
Additionally, we examine the scalability of the method in various settings such as the number of samples, the extent of block overlaps, and the number of blocks/parameters, using a newly designed problem structure. Due to the limited space, we omit the details of the problem structure and results but we discuss them in the online supplementary document. The results show that M-BC yields consistently good results in all settings in terms of calibration accuracy and uncertainty quantification. In particular, it demonstrates that the proposed method is suitable even for the problems with a significant data overlap. In the online supplementary document, we further examine the performance of M-BC when block domination occurs. Our approach provides superior results in terms of both calibration accuracy and uncertainty quantification across different levels of block dominance.
We evaluate the performance of our proposed approach in the BEM application. We use the hourly electrical energy consumption data and the simulation outputs obtained from the simulator, EnergyPlus 9.3.0 (U.S. Department of Energy 2019) on a municipal building of Mueller, Austin in Texas. EnergyPlus is a widely used BEM simulator to model both energy consumption-for heating, cooling, ventilation, lighting, plug and process loads-and water use in buildings. Each run for year-long simulation takes about 2 minutes with a standard desktop computer.
Among a large number of parameters employed in EnergyPlus, we choose important parameters in the building energy use, based on domain knowledge and previous studies (Manfren et al. 2013, Chong et al. 2017), including lighting, ventilation, domestic hot water, window material (optical properties), and heating and cooling systems, as summarized in Table 1. In this case study, we let θ g = (θ g,1 , . . . , θ g,6 ) T denote the vector of year-long global parameters and θ c = (θ c,1 , . . . , θ c,4 ) T and θ h = (θ h,1 , . . . , θ h,3 ) T , respectively, represent cooling-and heating-season parameter vectors.
We use the implementation setting similar to that in numerical examples. We terminate iterations when the relative difference between consecutive loss function values becomes smaller than a termination tolerance 10 -4 . We also use a small perturbation h = 10 -4 in computing the gradient.
The step size is decided using the backtracking line search. Figure 4 shows the example of loss function values for each seasonal block from the proposed M-BC approach. Each dashed vertical line indicates the iteration index where the inner iteration terminates, and each solid line marks the iteration index for the outer iteration. In the inner iteration, we first minimize the year-long loss F g , followed by F c and F h sequentially. Each plot in Figure 4 illustrates how these loss functions change over iterations. Interestingly, Figure 4(a) shows that F g decreases mostly when θ c is updated. This is because the cooling season during March to November (see Table 1) largely overlaps with the global block (January to December) in the studied Texas region. Similarly, the loss h during the heating season also decreases when θ c is calibrated. It is because heating and cooling seasons overlap during three months (see Table 1).
Overall, all of the three loss values converge within a few outer iterations, which indicates the important role of inner iterations that calibrates each block of parameters. LHS 1 through LHS 5. In all cases, M-BC achieves the smallest global (F g ), cooling-season (F c ), and heating-season (F h ) loss values (MSEs). In some cases (e.g., with LHS 3 and 4), the performance of H-BC is comparable to M-BC, but M-BC outperforms H-BC in other initial settings. S-BC performs worse than M-BC and H-BC in most cases. 3).
We further compare our approach with Bayesian calibration. To obtain posterior distributions, we exploit the No-U-Turn Sampler (NUTS), which is known to be useful for sampling from a high-dimensional posterior distribution as suggested in Chong and Menberg (2018). The NUTS is developed based on Hamiltonian Monte Carlo (HMC) which allows MCMC to converge faster in that it explores the high-dimensional posterior distribution by optimizing the tuning parameters such as a scaling factor and the number of leapfrog steps. More details are available in Chong and Menberg (2018). Even with NUTS, Bayesian calibration is computationally intensive, as discussed in Section 2. With the computational resources available to us, it takes an unrealistically long time if we use hourly year-long data. Thus, we use weekly data in the Bayesian analysis as in Kristensen et al. (2017b). Even with weekly data, one experiment requires about five days to calibrate the BEM parameters.
We consider two prior specifications in Bayesian approach: non-informative priors with uniform distributions and informative priors with Gaussian distributions, as suggested in Chong and Menberg (2018). In both prior settings, we use the settings in LHS 1 as the prior means for the parameters. We set the prior standard deviation as 0.2 in the informative prior settings, as in Chong and Menberg (2018).
Table 10 shows the MSEs in the test sets from 10 independent experiments where we use the posterior means as the calibrated estimates. Compared to M-BC, Bayesian calibration in both prior settings generates larger MSE values especially in the global and cooling-season blocks. We note that the MSE from M-BC in the heating-season block is slightly larger than that of the Bayesian approach with the Gaussian prior. One possible reason is that underlying heating-season loss is To further investigate, Figure 5 shows the posterior distributions of one of the parameterslightening level-when non-informative uniform and informative Gaussian priors are used. It turns out that the resulting posterior distribution when using the non-informative prior is still close to its prior distribution, showing large uncertainties. It implies that the calibrated parameter value from the posterior density does not deliver meaningful information. Meanwhile, when the informative prior is employed, the posterior distribution changes from its prior. Still, its MSEs are much larger in the global and cooling-season blocks (see Table 10).
In this study, we present a new multi-block parameter calibration methodology when computer model parameters need to be calibrated using different, possibly non-disjoint, subsets of data. We consider multiple loss functions, each for their associated block of parameters that use the corresponding dataset. In addition to the block property, our approach is different from the existing literature that generates the computer model responses at pre-selected parameter settings and builds emulators using the fixed dataset. Although being widely adopted, such traditional approach decouples the data generation stage and calibration stage. On the contrary, our approach adaptively runs the simulation at new parameter settings recommended by the nonlinear optimization procedure. As such, more instructive data is generated from the computer model, which makes the calibration procedure more effective.
Our implementation results in both numerical examples and BEM case study demonstrate that the discrepancy between the physical process and the computer model outputs can be significantly reduced, when the multi-block dependency of parameters is taken into consideration. Further, in our case study, the performance of the proposed M-BC algorithm is more robust, compared to the alternatives. We also construct CIs using asymptotic properties of ML estimators. Our approach yields more stable performance with narrower CIs while maintaining empirical coverage rates close to the nominal rates, compared to other methods. The comparison with the widely used Bayesian calibration elicits the advantage of our approach in terms of improved calibration accuracy.
In the future, we would like to extend the approach under more general settings, e.g., competing cases and time-variant processes (Byon et al. 2016). Developing multi-task learning and BO approaches that can handle overlapping data would also be a possible future research task. Additionally, we plan to use the well-calibrated BEM computer model to control building energy end use such as demand control (Jang et al. 2020, Li et al. 2020).