An example workflow of building energy model calibration
Hongyuan Jia
Adrian Chong
2021-07-16
Building energy model description
A short description has been given in the Synthetic meter data creation workflow.
Observed inputs, outputs and calibration parameters
We choose AMY (Annual Meteorological Year) as the observed input. AMY together with other separate local weather variables are commonly used as the observed inputs for calibration.
For calibration parameters, we select 7 commonly used ones for demonstration, including:
- Material properties. Here materials include the exterior wall insulation and exterior window
- Infiltration rate
- Occupancy density
- Lighting power density
- Electric equipment power density
- Cooling setpoint
- Heating setpoint
The mapping between those parameters and corresponding EnergyPlus inputs for the baseline model are listed in the table blow:
| # | Calibration parameter | Class | Object | Field | Value |
|---|---|---|---|---|---|
| 1 | Material properties | Material |
Steel Frame NonRes Wall Insulation |
Conductivity {W/m-K} |
0.049 |
| 2 | Density {kg/m3} |
265.0 | |||
| 3 | Specific Heat {J/kg-K} |
836.8 | |||
| 4 | WindowMaterial:SimpleGlazingSystem |
NonRes Fixed Assembly Window |
U-Factor {W/m2-K} |
3.237 | |
| 5 | Solar Heat Gain Coefficient |
0.39 | |||
| 6 | Infiltration rate | ZoneInfiltration:DesignFlowRate |
All *_Infiltration objects |
Flow per Exterior Surface Area {m3/s-m2} |
0.000302 |
| 7 | Occupancy density | People |
All objects for offices | Zone Floor Area per Person {m2/persion} |
18.58 |
| 8 | Lighting power density | Lights |
All *_Lights objects for offices |
Watts per Zone Floor Area {W/m2} |
10.76 |
| 9 | Electric equipment power sensity | ElectricEquipment |
All *_Equip objects for offices |
Watts per Zone Floor Area {W/m2} |
10.76 |
| 10 | Cooling setpoint | Schedule:Compact |
CLGSETP_SCH |
Field 6 and Field 13 |
24 |
| 11 | Heating setpoint | Schedule:Compact |
HTGSETP_SCH |
Field 6, Field 16 and Field 21 |
21 |
For observed output, we choose hourly building electricity consumption. The corresponding output variable in EnergyPlus is Electricity:Facility with Hourly reporting frequency.
Create the initial test model to be calibrated
Read the model
# path of the model
path <- here("data-raw/idf/RefBldgMediumOfficeNew2004_Chicago.idf")
# read the model as the starting point for the test model
init <- read_idf(path)Create the initial test model
In order to create a test model for calibration, we modified the values of 11 calibration parameters in the model.
For demonstration, we set the initial parameter values before calibration to a random value in predefined ranges.
# actual parameter values in the model
real_conductivity <- 0.049
real_density <- 265.0
real_specific_heat <- 836.8
real_u_value <- 3.237
real_SHGC <- 0.39
real_infiltration <- 0.000302
real_people <- 18.58
real_lpd <- 10.76
real_epd <- 10.76
real_clg_sp <- 24
real_htg_sp <- 21
# define the upper and lower offset
offset <- 0.50
lower <- 1 - offset
upper <- 1 + offset
# define the lower and upper limits of parameters
param <- dplyr::tribble(
~parameter , ~real , ~lower , ~upper ,
"conductivity" , real_conductivity , real_conductivity * lower , real_conductivity * upper ,
"density" , real_density , real_density * lower , real_density * upper ,
"specific_heat" , real_specific_heat , real_specific_heat * lower , real_specific_heat * upper ,
"u_value" , real_u_value , real_u_value * lower , real_u_value * upper ,
"SHGC" , real_SHGC , 0.1 , 0.6 ,
"infiltration" , real_infiltration , 0 , 0.0005 ,
"people" , real_people , 10 , 25 ,
"lpd" , real_lpd , 2 , 20 ,
"epd" , real_epd , 2 , 20 ,
"cooling" , real_clg_sp , 22 , 26 ,
"heating" , real_htg_sp , 15 , 21.5
)
# for reproducibility
set.seed(4)
# initial parameter values for calibration
param <- param %>%
mutate(init = map2_dbl(lower, upper, runif, n = 1L)) %>%
select(parameter, real, init, everything())The ranges and initial values of parameters are listed in the table below.
tbl_param <- param %>%
mutate(name = c(
"Insulation conductivity (W/(m K))",
"Insulation density (kg/m<sup>3</sup>)",
"Insulation specific heat (J/(kg K))",
"Window U value [W/(m<sup>2</sup> K)]",
"Window SHGC",
"Infiltration per floor area (m<sup>3</sup>/(s m<sup>2</sup>))",
"People density (m<sup>2</sup>/person)",
"Lighting power density (W/m<sup>2</sup>)",
"Equipment power density (W/m<sup>2</sup>)",
"Cooling setpoint (°C)",
"Heating setpoint (°C)")) %>%
select(name, parameter, everything())
tbl_param %>%
gt::gt() %>%
gt::cols_hide(vars(real)) %>%
gt::fmt_number(vars(init, lower, upper)) %>%
gt::tab_header(title = "List of calibration parameters and their ranges") %>%
gt::cols_label(
name = "Calibration parameter",
parameter = "Symbol",
init = "Initial value",
lower = "Min",
upper = "Max") %>%
gt::fmt_markdown(vars(name)) %>%
gt::tab_style(
style = gt::cell_text(weight = "bold"),
locations = gt::cells_column_labels(everything()))## Warning: The `.dots` argument of `group_by()` is deprecated as of dplyr 1.0.0.| List of calibration parameters and their ranges | ||||
|---|---|---|---|---|
| Calibration parameter | Symbol | Initial value | Min | Max |
Insulation conductivity (W/(m K)) |
conductivity | 0.05 | 0.02 | 0.07 |
Insulation density (kg/m3) |
density | 134.87 | 132.50 | 397.50 |
Insulation specific heat (J/(kg K)) |
specific_heat | 664.20 | 418.40 | 1,255.20 |
Window U value [W/(m2 K)] |
u_value | 2.52 | 1.62 | 4.86 |
Window SHGC |
SHGC | 0.51 | 0.10 | 0.60 |
Infiltration per floor area (m3/(s m2)) |
infiltration | 0.00 | 0.00 | 0.00 |
People density (m2/person) |
people | 20.87 | 10.00 | 25.00 |
Lighting power density (W/m2) |
lpd | 18.31 | 2.00 | 20.00 |
Equipment power density (W/m2) |
epd | 19.08 | 2.00 | 20.00 |
Cooling setpoint (°C) |
cooling | 22.29 | 22.00 | 26.00 |
Heating setpoint (°C) |
heating | 19.91 | 15.00 | 21.50 |
Below we made some necessary to the model, including adding the output meter and setting the begin year of RunPeriod objects based on our synthetic meter data.
After this, we update the model based on the initial parameter values generated above and saved the model as Init.idf in the data/idf folder.
# remove all existing outputs
init$`Output:Variable` <- NULL
init$`Output:Meter` <- NULL
# add hourly building electricity output meter
init$add(Output_Meter = list("Electricity:Facility", "Hourly"))
# make sure weater file is used
init$SimulationControl$Run_Simulation_for_Weather_File_Run_Periods <- "Yes"
# update RunPeriod to correctly indicate an AMY EPW file is used
init$RunPeriod$annual$set(
name = "Philadelphia 2014",
begin_year = 2014,
day_of_week_for_start_day = "Wednesday"
)
init_val <- function(var) param %>% filter(parameter == !!var) %>% pull()
update_insulation(init,
conductivity = init_val("conductivity"),
density = init_val("density"),
specific_heat = init_val("specific_heat")
)
update_window(init, u_value = init_val("u_value"), SHGC = init_val("SHGC"))
update_infiltration(init, flow_per_area = init_val("infiltration"))
update_people(init, people = init_val("people"))
update_lights(init, lpd = init_val("lpd"))
update_equip(init, epd = init_val("epd"))
update_setpoint(init, cooling = init_val("cooling"), heating = init_val("heating"))
# save
init$save(here("data/idf/Init.idf"), overwrite = TRUE)Evaluate the performance of the initial init model
Run initial simulation with TMY3 weather
# file path of TMY3 EPW file
path_epw_tmy3 <- here("data-raw/epw/TMY3/USA_PA_Philadelphia.Intl.AP.724080_TMY3.epw")
# run annual simulation
init$run(path_epw_tmy3, here("data/sim/Init"), echo = FALSE)Extract the simulation results
meter_init <- init %>% extract_electricity()Visualize the trend
# select one week in July to examine the discrepancy
plot_init <- meter_init %>% weekly_compare(7, 1)
ggsave(here("figures/Init.png"), plot_init, height = 6, width = 10, dpi = 300)
knitr::include_graphics(here("figures/Init.png"))
One week electricity comparison between synthetic data and simulation
Calculate the statistical indicators
stats_init <- meter_init %>% cal_stats()Current calibration creteria in ASHRAE Guideline 14 for NMBE and CVRMSE is 10% and 30% respectively.
The initial NMBE and CV(RMSE) is 51.86% and 67.67%, respectively.
Use AMY
# read the initial test model
amy <- read_idf(here("data/idf/Init.idf"))
# path of AMY EPW
path_epw_amy <- here("data-raw/epw/AMY/PA_PHILADELPHIA_720304_14-13.epw")
# run annual simulation
amy$run(path_epw_amy, here("data/sim/AMY"), echo = FALSE)Extract the simulation results
meter_amy <- amy %>% extract_electricity()Visualize the trend
# select one week in July to examine the discrepancy
plot_init <- meter_amy %>% weekly_compare(7, 1)
ggsave(here("figures/AMY.png"), plot_init, height = 6, width = 10, dpi = 300)
knitr::include_graphics(here("figures/AMY.png"))
One week electricity comparison between synthetic data and simulation after using AMY weather
Calculate the statistical indicators
stats_amy <- meter_amy %>% cal_stats()Ater using the AMY weather file, the NMBE and CV(RMSE) is 50.80% and 62.12%, respectively.
Use optimization to calibrate the model
# define setup values
MU <- 4L # number of individuals per generation
LAMBDA <- 4L # number of offspring
MAX_GEN <- 20L # maximum of generation
MAX_CVRMSE <- 0.3 # creteria for CVRMSE
MAX_NMBE <- 0.1 # creteria for NMBE
P_MUTATE <- 0.1 # mutation probability
P_RECOMB <- 0.7 # crossover probability
# define an ecr control object to store information on the objective function
control <- ecr::initECRControl(calib_fitness, n.objectives = 2, minimize = TRUE)
# initial a logger to store population and fitness
log <- ecr::initLogger(control, log.pop = TRUE)
# initial archive of Pareto front
pareto <- ecr::initParetoArchive(control)
# set evolutionary operators
# use no priors and set all parameters to follow uniform distribution
control <- ecr::registerECROperator(control, "mutate", ecr::mutUniform,
lower = param$lower, upper = param$upper)
control <- ecr::registerECROperator(control, "recombine", ecr::recCrossover)
control <- ecr::registerECROperator(control, "selectForMating", ecr::selSimple)
control <- ecr::registerECROperator(control, "selectForSurvival", ecr::selNondom)
# initialize population of MU random values
pop <- ecr::initPopulation(
MU, ecr::genReal, n.dim = nrow(param),
lower = param$lower, upper = param$upper)
# evaluate fitness of initial population
fit <- evaluate_fitness(control, pop, 1, workers = 4)
# update log
ecr::updateLogger(log, pop, fit, MU)
# run evolutionary loop
repeat {
# generate offspring
offspring <- ecr::generateOffspring(control, pop, fit, LAMBDA, p.recomb = P_RECOMB, p.mut = P_MUTATE)
# evaluate fitness of the offspring
fit_offspring <- evaluate_fitness(control, offspring, log$env$n.gens + 1L, workers = 4)
# do logging
ecr::updateLogger(log, offspring, fit_offspring, MU)
# prepare next generation
res <- ecr::replaceMuPlusLambda(control, pop, offspring, fit, fit_offspring)
# extract population and fitness
pop <- res$population
fit <- res$fitness
# update pareto archieve
ecr::updateParetoArchive(pareto, pop, fit)
# check whether terminator conditions are met
stop_obj <- ecr:::doTerminate(log, list(
# stop if CVRMSE <= 30% and NMBE <= 10%
stopOnMeetCreteria(cvrmse = MAX_CVRMSE, nmbe = MAX_NMBE),
# stop if max generation matched
ecr::stopOnIters(MAX_GEN)
))
# stop if creteria are met
if (length(stop_obj) > 0L) break
}
# extract optimization results
result <- ecr:::makeECRResult(control, log, pop, fit, stop_obj)Calibration results
Pareto front
After the optimization completed, we extracted all the population together with the fitness values, and saved them into a table in data/population.csv. The is_pareto column indicates whether current individual is a Pareto front or not.
# get unique Pareto front
front <- get_front(pareto)
# get all population together with fitness
population <- get_population(log)
# rename columns using parameter names
names(population)[seq_len(nrow(param)) + 2L] <- param$parameter
# save results
write_csv(population, here("data/population.csv"))The figure blow shows the distribution of the Pareto front. Since we use synthetic meta data and models that only vary several parameters, it only takes 3 generations for the NSGA-II algorithm to converge.
# plot Pareto front
plot_pareto <- population %>%
ggplot(aes(nmbe, cvrmse)) +
geom_rect(aes(xmin = 0, xmax = MAX_NMBE, ymin = 0, ymax = MAX_CVRMSE),
color = "grey70", size = 1, linetype = 2, fill = "grey90", alpha = 0.1) +
geom_point() +
geom_point(data = front, color = "red") +
scale_x_continuous("NMBE / %", labels = scales::label_percent(),
breaks = seq(0, 1, 0.1), expand = c(0, 0), limits = c(0, max(population$nmbe) * 1.05)) +
scale_y_continuous("CV(RMSE) / %", labels = scales::label_percent(),
breaks = seq(0, 1, 0.1), expand = c(0, 0), limits = c(0, max(population$cvrmse) * 1.05)) +
theme_bw()
ggsave(here::here("figures/pareto.png"), plot_pareto, height = 6, width = 6, dpi = 300)
knitr::include_graphics(here("figures/pareto.png"))
Pareto front of optimization
Distribution of calibrated parameters
population %>%
filter(is_pareto) %>%
distinct(!!!syms(c(param$parameter, "nmbe", "cvrmse"))) %>%
mutate(index = seq_len(n())) %>%
select(index, everything()) %>%
gt::gt() %>%
gt::fmt_number(all_of(param$parameter)) %>%
gt::fmt_percent(vars(nmbe, cvrmse)) %>%
gt::cols_label(
index = "#",
conductivity = gt::html("Insulation conductivity (W/(m K))"),
density = gt::html("Insulation density (kg/m<sup>3</sup>)"),
specific_heat = gt::html("Insulation specific heat (J/(kg K))"),
u_value = gt::html("Window U value [W/(m<sup>2</sup> K)]"),
SHGC = gt::html("Window SHGC"),
infiltration = gt::html("Infiltration per floor area (m<sup>3</sup>/(s m<sup>2</sup>))"),
people = gt::html("People density (m<sup>2</sup>/person)"),
lpd = gt::html("Lighting power density (W/m<sup>2</sup>)"),
epd = gt::html("Equipment power density (W/m<sup>2</sup>)"),
cooling = gt::html("Cooling setpoint (°C)"),
heating = gt::html("Heating setpoint (°C)"),
nmbe = "NMBE",
cvrmse = "CV(RMSE)") %>%
gt::tab_style(
style = gt::cell_text(weight = "bold"),
locations = gt::cells_column_labels(everything()))| # | Insulation conductivity (W/(m K)) | Insulation density (kg/m3) | Insulation specific heat (J/(kg K)) | Window U value [W/(m2 K)] | Window SHGC | Infiltration per floor area (m3/(s m2)) | People density (m2/person) | Lighting power density (W/m2) | Equipment power density (W/m2) | Cooling setpoint (°C) | Heating setpoint (°C) | NMBE | CV(RMSE) |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | 0.07 | 204.81 | 904.41 | 4.72 | 0.23 | 0.00 | 14.09 | 10.05 | 8.53 | 22.63 | 19.03 | 8.12% | 11.84% |
| 2 | 0.03 | 204.81 | 904.41 | 4.72 | 0.54 | 0.00 | 11.66 | 15.46 | 10.41 | 25.50 | 15.30 | 0.39% | 29.71% |
population %>%
filter(is_pareto) %>%
distinct(!!!syms(c(param$parameter, "nmbe", "cvrmse")), .keep_all = TRUE) %>%
select(index_gen, index_ind, nmbe, cvrmse) %>%
mutate(path_plot = here(sprintf("figures/Gen%i/Gen%i_Ind%i.png", index_gen, index_gen, index_ind))) %>%
pull(path_plot) %>%
knitr::include_graphics()
Session Info
sessionInfo()## R version 4.1.0 (2021-05-18)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19042)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=English_United States.1252
## [2] LC_CTYPE=English_United States.1252
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.1252
## system code page: 936
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] here_1.0.1 eplusr_0.14.2 lubridate_1.7.10 forcats_0.5.1
## [5] stringr_1.4.0 dplyr_1.0.5 purrr_0.3.4 readr_1.4.0
## [9] tidyr_1.1.3 tibble_3.1.0 ggplot2_3.3.3 tidyverse_1.3.0
##
## loaded via a namespace (and not attached):
## [1] fs_1.5.0 bit64_4.0.5 RColorBrewer_1.1-2 progress_1.2.2
## [5] httr_1.4.2 rprojroot_2.0.2 ecr_2.1.0 tools_4.1.0
## [9] backports_1.2.1 utf8_1.2.1 R6_2.5.0 lazyeval_0.2.2
## [13] DBI_1.1.1 colorspace_2.0-0 withr_2.4.1 tidyselect_1.1.0
## [17] prettyunits_1.1.1 processx_3.5.0 mco_1.15.6 bit_4.0.4
## [21] compiler_4.1.0 parallelMap_1.5.0 cli_2.3.1 rvest_1.0.0
## [25] gt_0.2.2 xml2_1.3.2 plotly_4.9.4.1 labeling_0.4.2
## [29] sass_0.3.1 scales_1.1.1 checkmate_2.0.0 plot3D_1.3
## [33] callr_3.5.1 commonmark_1.7 digest_0.6.27 rmarkdown_2.7
## [37] smoof_1.6.0.2 pkgconfig_2.0.3 htmltools_0.5.1.1 parallelly_1.24.0
## [41] dbplyr_2.1.0 fastmap_1.1.0 highr_0.8 htmlwidgets_1.5.3
## [45] rlang_0.4.10 readxl_1.3.1 rstudioapi_0.13 RSQLite_2.2.4
## [49] BBmisc_1.11 generics_0.1.0 farver_2.1.0 jsonlite_1.7.2
## [53] magrittr_2.0.1 Rcpp_1.0.6 munsell_0.5.0 fansi_0.4.2
## [57] lifecycle_1.0.0 stringi_1.5.3 yaml_2.2.1 RJSONIO_1.3-1.4
## [61] plyr_1.8.6 misc3d_0.9-0 grid_4.1.0 blob_1.2.1
## [65] listenv_0.8.0 parallel_4.1.0 crayon_1.4.1 haven_2.3.1
## [69] hms_1.0.0 knitr_1.31 ps_1.6.0 pillar_1.5.1
## [73] tcltk_4.1.0 future.apply_1.7.0 codetools_0.2-18 reshape2_1.4.4
## [77] fastmatch_1.1-0 reprex_1.0.0 glue_1.4.2 evaluate_0.14
## [81] ParamHelpers_1.14 data.table_1.14.0 modelr_0.1.8 vctrs_0.3.6
## [85] cellranger_1.1.0 gtable_0.3.0 future_1.21.0 assertthat_0.2.1
## [89] cachem_1.0.4 xfun_0.21 formattable_0.2.1 broom_0.7.5
## [93] viridisLite_0.3.0 memoise_2.0.0 units_0.7-1 globals_0.14.0
## [97] ellipsis_0.3.1