Building energy model description

Observed inputs, outputs and calibration parameters

Create the initial test model to be calibrated

Read the model
Create the initial test model

Evaluate the performance of the initial init model

Run initial simulation with TMY3 weather
Extract the simulation results
Visualize the trend
Calculate the statistical indicators

Use AMY

Extract the simulation results
Visualize the trend
Calculate the statistical indicators

Use optimization to calibrate the model

Calibration results

Pareto front
Distribution of calibrated parameters

Session Info

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 (&deg;C)",
        "Heating setpoint (&deg;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/m³)	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/(m² K)]	u_value	2.52	1.62	4.86
Window SHGC	SHGC	0.51	0.10	0.60
Infiltration per floor area (m³/(s m²))	infiltration	0.00	0.00	0.00
People density (m²/person)	people	20.87	10.00	25.00
Lighting power density (W/m²)	lpd	18.31	2.00	20.00
Equipment power density (W/m²)	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 (&deg;C)"),
        heating = gt::html("Heating setpoint (&deg;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/m³)	Insulation specific heat (J/(kg K))	Window U value [W/(m² K)]	Window SHGC	Infiltration per floor area (m³/(s m²))	People density (m²/person)	Lighting power density (W/m²)	Equipment power density (W/m²)	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

An example workflow of building energy model calibration

Hongyuan Jia

Adrian Chong

2021-07-16

Building energy model description

Observed inputs, outputs and calibration parameters

Create the initial test model to be calibrated

Read the model

Create the initial test model

Evaluate the performance of the initial init model

Run initial simulation with TMY3 weather

Extract the simulation results

Visualize the trend

Calculate the statistical indicators

Use AMY

Extract the simulation results

Visualize the trend

Calculate the statistical indicators

Use optimization to calibrate the model

Calibration results

Pareto front

Distribution of calibrated parameters

Session Info