Conduct sensitivity analysis and Bayesian calibration of EnergyPlus models.
Currently, epluspar is not on CRAN. You can install the development version from GitHub.
install.packages("epluspar", repos = "https://hongyuanjia.r-universe.dev")The basic workflow is basically:
- Adding parameters for sensitivity analysis using
$param()or$apply_measure(). - Check parameter sampled values and generated parametric models using
$samples()and$models()respectively. - Run EnergyPlus simulations in parallel using
$run(). - Gather EnergyPlus simulated data using
$report_data()or$tabular_data(). - Evaluate parameter sensitivity using
$evaluate().
Create a SensitivityJob object:
# use an example file from EnergyPlus v8.8 for demonstration here
path_idf <- file.path(eplusr::eplus_config(8.8)$dir, "ExampleFiles", "5Zone_Transformer.idf")
path_epw <- file.path(eplusr::eplus_config(8.8)$dir, "WeatherData", "USA_CA_San.Francisco.Intl.AP.724940_TMY3.epw")
# create a `SensitivityJob` class which inheris from eplusr::ParametricJob class
sen <- sensi_job(path_idf, path_epw)
#> Adding an object in class 'Output:SQLite' and setting its 'Option Type' to 'SimpleAndTabular' in order to create SQLite output file.
#> Moving all objects in class 'Output:Meter:MeterFileOnly' to 'Output:Meter' in order to use csv for data extraction.Set sensitivity parameters using $param() or $apply_measure().
- Using
$param()
# set parameter using similar syntax to `Idf$set()` in eplusr
sen$param(
# For adding a single object field as parameter
# Syntax: Object = list(Field = c(Min, Max, Levels))
`Supply Fan 1` = list(Fan_Total_Efficiency = c(0.1, 1.0, 5)),
# For adding a class field as parameter
Material := list(
Thickness = c(min = 0.01, max = 0.08, levels = 5),
Conductivity = c(min = 0.01, max = 0.6, levels = 6)
),
# use `.names` to give names to each parameter
.names = c("efficiency", "thickness", "conductivity"),
# See `r` and `grid_jump` in `sensitivity::morris`
.r = 8, .grid_jump = 1
)- Using
$apply_measure()
# first define a "measure"
my_actions <- function (idf, efficiency, thickness, conducitivy) {
idf$set(
`Supply Fan 1` = list(Fan_Total_Efficiency = efficiency),
Material := list(Thickness = thickness, Conductivity = conducitivy)
)
idf
}
# then apply that measure with parameter space definitions as function arguments
sen$apply_measure(my_actions,
efficiency = c(0.1, 1.0, 5),
thickness = c(0.01, 0.08, 5),
conducitivy = c(0.1, 0.6, 6),
.r = 8, .grid_jump = 1
)
#> -- EnergPlus Sensitivity Analysis Job ------------------------------------
#> * Seed: 'C:\Users\hongy\AppData\Local\EnergyPlusV8-8-0\ExampleFiles\5...
#> * Weather: 'C:\Users\hongy\AppData\Local\EnergyPlusV8-8-0\WeatherData...
#> * EnergyPlus Version: '8.8.0'
#> * EnergyPlus Path: 'C:\Users\hongy\AppData\Local\EnergyPlusV8-8-0'
#> Applied Measure: 'my_actions'
#> Parameters [3]
#> [1]: 'efficiency' [0.1, 1] (lvl: 5)
#> [2]: 'thickness' [0.01, 0.08] (lvl: 5)
#> [3]: 'conducitivy' [0.1, 0.6] (lvl: 6)
#> Parametric Models [32]:
#> [01]: 'Case01.idf'
#> [02]: 'Case02.idf'
#> [03]: 'Case03.idf'
#> [04]: 'Case04.idf'
#> [05]: 'Case05.idf'
#> [06]: 'Case06.idf'
#> [07]: 'Case07.idf'
#> [08]: 'Case08.idf'
#> [09]: 'Case09.idf'
....Get samples
sen$samples()
#> case efficiency thickness conducitivy
#> <int> <num> <num> <num>
#> 1: 1 0.325 0.0275 0.4
#> 2: 2 0.100 0.0275 0.4
#> 3: 3 0.100 0.0100 0.4
#> 4: 4 0.100 0.0100 0.3
#> 5: 5 0.100 0.0625 0.3
#> 6: 6 0.325 0.0625 0.3
#> 7: 7 0.325 0.0800 0.3
#> 8: 8 0.325 0.0800 0.4
#> 9: 9 0.100 0.0800 0.3
#> 10: 10 0.325 0.0800 0.3
#> 11: 11 0.325 0.0800 0.2
#> 12: 12 0.325 0.0625 0.2
#> 13: 13 0.100 0.0450 0.2
#> 14: 14 0.325 0.0450 0.2
#> 15: 15 0.325 0.0450 0.1
#> 16: 16 0.325 0.0625 0.1
#> 17: 17 0.550 0.0100 0.5
#> 18: 18 0.775 0.0100 0.5
....Run simulations and calculate statistic indicators
# run simulations in temporary directory
sen$run(dir = tempdir(), echo = FALSE)
#> ── EnergPlus Parametric Simulation Job ─────────────────────────────────────────
#> * Seed: '/home/hongyuanjia/.local/EnergyPlus-8-8-0/ExampleFiles/5Zone_Trans...
#> * Weather: '/home/hongyuanjia/.local/EnergyPlus-8-8-0/WeatherData/USA_CA_Sa...
#> -- EnergPlus Sensitivity Analysis Job ------------------------------------
#> * Seed: 'C:\Users\hongy\AppData\Local\EnergyPlusV8-8-0\ExampleFiles\5...
#> * Weather: 'C:\Users\hongy\AppData\Local\EnergyPlusV8-8-0\WeatherData...
#> * EnergyPlus Version: '8.8.0'
#> * EnergyPlus Path: 'C:\Users\hongy\AppData\Local\EnergyPlusV8-8-0'
#> Applied Measure: 'my_actions'
#> Parameters [3]
#> [1]: 'efficiency' [0.1, 1] (lvl: 5)
#> [2]: 'thickness' [0.01, 0.08] (lvl: 5)
#> [3]: 'conducitivy' [0.1, 0.6] (lvl: 6)
#> Parametric Models [32]:
#> [01]: 'Case01.idf' <-- SUCCEEDED
#> [02]: 'Case02.idf' <-- SUCCEEDED
#> [03]: 'Case03.idf' <-- SUCCEEDED
#> [04]: 'Case04.idf' <-- SUCCEEDED
#> [05]: 'Case05.idf' <-- SUCCEEDED
#> [06]: 'Case06.idf' <-- SUCCEEDED
#> [07]: 'Case07.idf' <-- SUCCEEDED
#> [08]: 'Case08.idf' <-- SUCCEEDED
#> [09]: 'Case09.idf' <-- SUCCEEDED
....
# extract output
# here is just am example
eng <- sen$tabular_data(table_name = "site and source energy",
column_name = "energy per total building area",
row_name = "total site energy")[, as.numeric(value)]
# calculate sensitivity
(result <- sen$evaluate(eng))
#>
#> Call:
#> sensitivity::morris(model = NULL, factors = l$param$param_name, r = .r, design = list(type = "oat", levels = l$param$levels, grid.jump = .grid_jump), binf = l$param$min, bsup = l$param$max, scale = .scale)
#>
#> Model runs: 32
#> mu mu.star sigma
#> efficiency -0.88000 0.88000 0.80796040
#> thickness -0.51000 0.51000 0.36893670
#> conducitivy 0.16875 0.16875 0.09613049
# extract data
attr(result, "data")
#> index name mu mu.star sigma
#> <int> <char> <num> <num> <num>
#> 1: 1 efficiency -0.88000 0.88000 0.80796040
#> 2: 2 thickness -0.51000 0.51000 0.36893670
#> 3: 3 conducitivy 0.16875 0.16875 0.09613049Plot
# plot
plot(result)
- Setting input and output variables using
$input()and$output()respectively. Input variables should be variables listed in RDD while output variables should be variables listed in RDD and MDD. - Adding parameters to calibrate using
$param()or$apply_measure(). - Check parameter sampled values and generated parametric models using
$samples()and$models()respectively. - Run EnergyPlus simulations in parallel using
$eplus_run(). - Gather simulated data of input and output parameters using
$data_sim(). - Specify field measured data of input and output parameters using
$data_field(). - Specify field measured data of input and output parameters using
$data_field(). - Specify input data for Stan for Bayesian calibration using
$data_bc(). - Run bayesian calibration and get predictions using stan using
$stan_run().
Create a BayesCalibJob object:
# use an example file from EnergyPlus v8.8 for demonstration here
path_idf <- file.path(eplusr::eplus_config(8.8)$dir, "ExampleFiles", "RefBldgLargeOfficeNew2004_Chicago.idf")
path_epw <- file.path(eplusr::eplus_config(8.8)$dir, "WeatherData", "USA_CA_San.Francisco.Intl.AP.724940_TMY3.epw")
# create a `SensitivityJob` class which inherits from eplusr::ParametricJob class
bc <- bayes_job(path_idf, path_epw)$read_rdd() and $read_mdd() can be used to get RDD and MDD for
current seed model.
(rdd <- bc$read_rdd())
#> Initializing RDD...
#> Initializing RDD... [SUCCESSFUL]
#> == EnergyPlus Report Data Dictionary File ================================
#> * EnergyPlus version: 8.8.0 (7c3bbe4830)
#> * Simulation started: 2020-09-10 18:55:00
#>
#> -- Details ---------------------------------------------------------------
#> index reported_time_step report_type
#> <int> <char> <char>
#> 1: 1 Zone Average
#> 2: 2 Zone Average
#> 3: 3 Zone Average
#> 4: 4 Zone Average
#> 5: 5 Zone Average
#> ---
#> 616: 616 HVAC Average
#> 617: 617 HVAC Sum
#> 618: 618 HVAC Average
#> 619: 619 HVAC Sum
#> 620: 620 HVAC Average
#> variable
#> <char>
....
(mdd <- bc$read_mdd())
#> == EnergyPlus Meter Data Dictionary File =================================
#> * EnergyPlus version: 8.8.0 (7c3bbe4830)
#> * Simulation started: 2020-09-10 18:55:00
#>
#> -- Details ---------------------------------------------------------------
#> index reported_time_step report_type
#> <int> <char> <char>
#> 1: 1 Zone Meter
#> 2: 2 Zone Meter
#> 3: 3 Zone Meter
#> 4: 4 Zone Meter
#> 5: 5 Zone Meter
#> ---
#> 173: 173 Zone Meter
#> 174: 174 Zone Meter
#> 175: 175 Zone Meter
#> 176: 176 Zone Meter
#> 177: 177 Zone Meter
#> variable units
#> <char> <char>
....Input variables and output variables can be set by using $input() and
$output(), respectively. For $input(), only variables listed in RDD
are supported. For $output(), variables listed in RDD and MDD are both
supported.
By default, they are all empty and $input(), $output() will return
NULL.
bc$input()
#> NULL
bc$output()
#> NULL
bc$models()
#> Unable to create parametric models because input variables are not set. Please use '$input()' to set input variables.
#> NULLYou can specify input and output parameters using RddFile, MddFile
and data.frames.
# using RDD and MDD
bc$input(rdd[1:3])
#> index class key_value variable_name
#> <int> <char> <char> <char>
#> 1: 1 Output:Variable * Site Outdoor Air Drybulb Temperature
#> 2: 2 Output:Variable * Site Outdoor Air Dewpoint Temperature
#> 3: 3 Output:Variable * Site Outdoor Air Wetbulb Temperature
#> reporting_frequency
#> <char>
#> 1: Timestep
#> 2: Timestep
#> 3: Timestep
bc$output(mdd[1])
#> index class key_value variable_name reporting_frequency
#> <int> <char> <char> <char> <char>
#> 1: 1 Output:Meter <NA> Electricity:Facility Timestep
# using data.frame
bc$input(eplusr::rdd_to_load(rdd[1:3]))
#> index class key_value variable_name
#> <int> <char> <char> <char>
#> 1: 1 Output:Variable * Site Outdoor Air Drybulb Temperature
#> 2: 2 Output:Variable * Site Outdoor Air Dewpoint Temperature
#> 3: 3 Output:Variable * Site Outdoor Air Wetbulb Temperature
#> reporting_frequency
#> <char>
#> 1: Timestep
#> 2: Timestep
#> 3: Timestep
bc$output(eplusr::mdd_to_load(mdd[1]))
#> index class key_value variable_name reporting_frequency
#> <int> <char> <char> <char> <char>
#> 1: 1 Output:Meter <NA> Electricity:Facility TimestepYou can set append to NULL to remove all existing input and output
parameters.
bc$input(append = NULL)
#> NULL
bc$output(append = NULL)
#> NULLYou can also directly specify variable names:
bc$input("CoolSys1 Chiller 1", paste("chiller evaporator", c("inlet temperature", "outlet temperature", "mass flow rate")), "hourly")
#> index class key_value
#> <int> <char> <char>
#> 1: 1 Output:Variable CoolSys1 Chiller 1
#> 2: 2 Output:Variable CoolSys1 Chiller 1
#> 3: 3 Output:Variable CoolSys1 Chiller 1
#> variable_name reporting_frequency
#> <char> <char>
#> 1: Chiller Evaporator Inlet Temperature Hourly
#> 2: Chiller Evaporator Outlet Temperature Hourly
#> 3: Chiller Evaporator Mass Flow Rate Hourly
bc$output("CoolSys1 Chiller 1", "chiller electric power", "hourly")
#> index class key_value variable_name
#> <int> <char> <char> <char>
#> 1: 1 Output:Variable CoolSys1 Chiller 1 Chiller Electric Power
#> reporting_frequency
#> <char>
#> 1: HourlyNote that variable cannot be set as both an input and output variable.
bc$output("CoolSys1 Chiller 1", name = "chiller evaporator inlet temperature", reporting_frequency = "hourly")
#> Error: Variables specified have already been set as input: 'CoolSys1 Chiller 1:Chiller Evaporator Inlet Temperature'Also, note that input and output variables should have the same reporting frequency.
bc$output(mdd[1], reporting_frequency = "daily")
#> Error: Object specified does not have the same reporting frequency as existing input ('Hourly'). Invalid input reporting frequency: 'daily'For $output(), both variables in RDD and MDD are supported. However,
for $input(), only variables in RDD are allowed.
Similarly like SensitivityJob, parameters can be added using either
$param() or $apply_measure().
Here use $param() for demonstration. Basically there are 3 format of
defining a parameter:
-
object = list(field1 = c(min, max), field2 = c(min, max), ...)This is the basic format.
field1andfield2inobjectwill be added as two different parameters, with minimum and maximum value specified asminandmax. -
class := list(field1 = c(min, max), field2 = c(min, max), ...)This is useful when you want to treat
field1andfield2in all objects inclassas two different parameters. Please note the use of special notion of:=instead of=. -
.(objects) := list(field1 = c(min, max), field2 = c(min, max), ...)Sometimes you may not want to treat a field in all objects in a class but only a subset of objects. You can use a special notation on the left hand side
.(). In the parentheses can be object names or IDs.
bc$param(
`CoolSys1 Chiller 1` = list(reference_cop = c(4, 6), reference_capacity = c(2.5e6, 3.0e6)),
.names = c("cop1", "cap1"), .num_sim = 5
)
#> -- EnergPlus Parametric Simulation Job -----------------------------------
#> * Seed: 'C:\Users\hongy\AppData\Local\EnergyPlusV8-8-0\ExampleFiles\R...
#> * Weather: 'C:\Users\hongy\AppData\Local\EnergyPlusV8-8-0\WeatherData...
#> * EnergyPlus Version: '8.8.0'
#> * EnergyPlus Path: 'C:\Users\hongy\AppData\Local\EnergyPlusV8-8-0'
#> Applied Measure: ''
#> Parametric Models [5]:
#> [1]: 'Case1.idf'
#> [2]: 'Case2.idf'
#> [3]: 'Case3.idf'
#> [4]: 'Case4.idf'
#> [5]: 'Case5.idf'
#> << Job has not been run before >>Parameter values can be retrieved using $samples().
bc$samples()
#> case cop1 cap1
#> <int> <num> <num>
#> 1: 1 5.568475 2581922
#> 2: 2 4.806140 2875407
#> 3: 3 4.794320 2753627
#> 4: 4 4.145201 2628566
#> 5: 5 5.923980 2946315Generated Idfs can be retrieved using $models().
names(bc$models())
#> [1] "Case1" "Case2" "Case3" "Case4" "Case5"$eplus_run() runs all parametric models in parallel. Parameter
run_period can be given to insert a new RunPeriod object. In this
case, all existing RunPeriod objects in the seed model will be
commented out.
bc$eplus_run(dir = tempdir(), run_period = list("example", 7, 1, 7, 3), echo = FALSE)
#> Reset 'Run Simulation for Weather File Run Periods' in 'SimulationControl' from 'No' to 'Yes' to make sure input run period can take effect.
#> Reset 'Run Simulation for Weather File Run Periods' in 'SimulationControl' from 'No' to 'Yes' to make sure input run period can take effect.
#> Reset 'Run Simulation for Weather File Run Periods' in 'SimulationControl' from 'No' to 'Yes' to make sure input run period can take effect.
#> Reset 'Run Simulation for Weather File Run Periods' in 'SimulationControl' from 'No' to 'Yes' to make sure input run period can take effect.
#> Reset 'Run Simulation for Weather File Run Periods' in 'SimulationControl' from 'No' to 'Yes' to make sure input run period can take effect.
#> ── EnergPlus Parametric Simulation Job ─────────────────────────────────────────
#> * Seed: '/home/hongyuanjia/.local/EnergyPlus-8-8-0/ExampleFiles/RefBldgLarg...
#> * Weather: '/home/hongyuanjia/.local/EnergyPlus-8-8-0/WeatherData/USA_CA_Sa...
#> -- EnergPlus Parametric Simulation Job -----------------------------------
#> * Seed: 'C:\Users\hongy\AppData\Local\EnergyPlusV8-8-0\ExampleFiles\R...
#> * Weather: 'C:\Users\hongy\AppData\Local\EnergyPlusV8-8-0\WeatherData...
#> * EnergyPlus Version: '8.8.0'
#> * EnergyPlus Path: 'C:\Users\hongy\AppData\Local\EnergyPlusV8-8-0'
#> Applied Measure: ''
#> Parametric Models [5]:
#> [1]: 'Case1.idf' <-- SUCCEEDED
#> [2]: 'Case2.idf' <-- SUCCEEDED
#> [3]: 'Case3.idf' <-- SUCCEEDED
#> [4]: 'Case4.idf' <-- SUCCEEDED
#> [5]: 'Case5.idf' <-- SUCCEEDED
#> Simulation started at '2020-09-10 18:55:43' and completed successfully after 13.28 secs.$data_sim() returns a data.table (when merge is TRUE) or a list
of 2 data.table (when merge is FALSE) which contains the simulated
data of input and output parameters. These data will be stored
internally and used during Bayesian calibration using Stan.
The resolution parameter can be used to specify the time resolution of
returned data. Note that input time resolution cannot be smaller than
the reporting frequency, otherwise an error will be issued.
bc$data_sim("1 min")
#> Error: Input resolution should not be smaller than reporting frequency (1 hour). Invalid resolution found: '1 min'bc$data_sim("6 hour")
#> $input
#> case Date/Time
#> <int> <char>
#> 1: 1 07/01 06:00
#> 2: 1 07/01 12:00
#> 3: 1 07/01 18:00
#> 4: 1 07/01 24:00
#> 5: 1 07/02 06:00
#> 6: 1 07/02 12:00
#> 7: 1 07/02 18:00
#> 8: 1 07/02 24:00
#> 9: 1 07/03 06:00
#> 10: 1 07/03 12:00
#> 11: 1 07/03 18:00
#> 12: 1 07/03 24:00
#> 13: 2 07/01 06:00
#> 14: 2 07/01 12:00
#> 15: 2 07/01 18:00
#> 16: 2 07/01 24:00
#> 17: 2 07/02 06:00
....$data_field() takes a data.frame of measured value of output
parameters and returns a list of data.tables which contains the
measured value of input and output parameters, and newly measured value
of input if applicable.
For input parameters, the values of simulation data for the first case are directly used as the measured values.
For demonstration, we use the seed model to generate fake measured output data.
# clone the seed model
seed <- bc$seed()$clone()
# remove existing RunPeriod objects
seed$RunPeriod <- NULL
# set run period as the same as in `$eplus_run()`
seed$add(RunPeriod = list("test", 7, 1, 7, 3))
#> $test
#> <IdfObject: 'RunPeriod'> [ID:846] `test`
#> Class: <RunPeriod>
#> +- 01 : "test", !- Name
#> |- 02*: 7, !- Begin Month
#> |- 03*: 1, !- Begin Day of Month
#> |- 04*: 7, !- End Month
#> |- 05*: 3, !- End Day of Month
#> |- 06 : "UseWeatherFile", !- Day of Week for Start Day
#> |- 07 : "Yes", !- Use Weather File Holidays and Special Da...
#> |- 08 : "Yes", !- Use Weather File Daylight Saving Period
#> |- 09 : "No", !- Apply Weekend Holiday Rule
#> |- 10 : "Yes", !- Use Weather File Rain Indicators
#> \- 11 : "Yes"; !- Use Weather File Snow Indicators
seed$SimulationControl$set(
`Run Simulation for Sizing Periods` = "No",
`Run Simulation for Weather File Run Periods` = "Yes"
)
#> <IdfObject: 'SimulationControl'> [ID:2]
#> Class: <SimulationControl>
#> +- 1: "YES", !- Do Zone Sizing Calculation
#> |- 2: "YES", !- Do System Sizing Calculation
#> |- 3: "YES", !- Do Plant Sizing Calculation
#> |- 4: "No", !- Run Simulation for Sizing Periods
#> |- 5: "Yes", !- Run Simulation for Weather File Run Periods
#> |- 6: "YES", !- Do HVAC Sizing Simulation for Sizing Periods
#> \- 7: 2; !- Maximum Number of HVAC Sizing Simulation Passes
# save the model to tempdir
seed$save(tempfile(fileext = ".idf"))
# run
job <- seed$run(bc$weather(), echo = FALSE)
# get output data
fan_power <- epluspar:::report_dt_aggregate(job$report_data(name = bc$output()$variable_name, all = TRUE, day_type = "normalday"), "6 hour")
fan_power <- epluspar:::report_dt_to_wide(fan_power)
# add Gaussian noice
fan_power <- fan_power[, -"Date/Time"][
, lapply(.SD, function (x) x + rnorm(length(x), sd = 0.05 * sd(x)))][
, lapply(.SD, function (x) {x[x < 0] <- 0; x})
]
# set field data
bc$data_field(fan_power)
#> $input
#> case Date/Time
#> <int> <char>
#> 1: NA 07/01 06:00
#> 2: NA 07/01 12:00
#> 3: NA 07/01 18:00
#> 4: NA 07/01 24:00
#> 5: NA 07/02 06:00
#> 6: NA 07/02 12:00
#> 7: NA 07/02 18:00
#> 8: NA 07/02 24:00
#> 9: NA 07/03 06:00
#> 10: NA 07/03 12:00
#> 11: NA 07/03 18:00
#> 12: NA 07/03 24:00
#> COOLSYS1 CHILLER 1:Chiller Evaporator Inlet Temperature [C](6 Hour)
#> <num>
#> 1: 8.186813
#> 2: 8.065758
#> 3: 8.195386
....$data_bc() takes a list of field data and simulated data, and returns
a list that contains data input for Bayesian calibration using the Stan
model
n: Number of measured parameter observations.m: Number of simulated observations.p: Number of input parameters.q: Number of calibration parameters.yf: Data of measured output after z-score standardization using data of simulated output.yc: Data of simulated output after z-score standardization.xf: Data of measured input after min-max normalization.xc: Data of simulated input after min-max normalization.tc: Data of calibration parameters after min-max normalization.
str(bc$data_bc())
#> List of 11
#> $ n : int 12
#> $ n_pred: int 12
#> $ m : int 60
#> $ p : int 3
#> $ q : int 2
#> $ yf : num [1:12] -0.879 1.111 1.498 -0.848 -0.879 ...
#> $ yc : num [1:60] -0.876 1.149 1.648 -0.879 -0.879 ...
#> $ xf :Classes 'data.table' and 'data.frame': 12 obs. of 3 variables:
#> ..$ COOLSYS1 CHILLER 1:Chiller Evaporator Inlet Temperature [C](6 Hour) : num [1:12] 0.251854 0.000112 0.269682 0.25723 0.25723 ...
#> ..$ COOLSYS1 CHILLER 1:Chiller Evaporator Outlet Temperature [C](6 Hour): num [1:12] 0.816 0 0.0627 0.8792 0.8792 ...
#> ..$ COOLSYS1 CHILLER 1:Chiller Evaporator Mass Flow Rate [kg/s](6 Hour) : num [1:12] 0.00247 0.87075 1 0 0 ...
#> ..- attr(*, ".internal.selfref")=<externalptr>
#> $ xc :Classes 'data.table' and 'data.frame': 60 obs. of 3 variables:
#> ..$ COOLSYS1 CHILLER 1:Chiller Evaporator Inlet Temperature [C](6 Hour) : num [1:60] 0.251854 0.000112 0.269682 0.25723 0.25723 ...
#> ..$ COOLSYS1 CHILLER 1:Chiller Evaporator Outlet Temperature [C](6 Hour): num [1:60] 0.816 0 0.0627 0.8792 0.8792 ...
#> ..$ COOLSYS1 CHILLER 1:Chiller Evaporator Mass Flow Rate [kg/s](6 Hour) : num [1:60] 0.00247 0.87075 1 0 0 ...
#> ..- attr(*, ".internal.selfref")=<externalptr>
#> $ x_pred:Classes 'data.table' and 'data.frame': 12 obs. of 3 variables:
#> ..$ COOLSYS1 CHILLER 1:Chiller Evaporator Inlet Temperature [C](6 Hour) : num [1:12] 0.251854 0.000112 0.269682 0.25723 0.25723 ...
....You can also supply your own field data and simulated data and using
$data_bc() to construct the input for the Stan model. Input
data_field and data_sim should have the same structure as the output
from $data_field() and $data_sim(). If data_field and data_sim
is not specified, the output from $data_field() and $data_sim() will
be used.
You can save the internal Stan code using $stan_file(). If no path is
specified, a character vector that contains the stan code will be
returned.
bc$stan_file()
#> [1] "/*"
#> [2] " The MIT License (MIT)"
#> [3] ""
#> [4] " Copyright (c) 2019-2020 Hongyuan Jia and Adrian Chong"
#> [5] ""
#> [6] " Permission is hereby granted, free of charge, to any person obtaining a copy"
#> [7] " of this software and associated documentation files (the \"Software\"), to deal"
#> [8] " in the Software without restriction, including without limitation the rights"
#> [9] " to use, copy, modify, merge, publish, distribute, sublicense, and/or sell"
#> [10] " copies of the Software, and to permit persons to whom the Software is"
#> [11] " furnished to do so, subject to the following conditions:"
#> [12] ""
#> [13] " The above copyright notice and this permission notice shall be included in"
#> [14] " all copies or substantial portions of the Software."
#> [15] ""
#> [16] " THE SOFTWARE IS PROVIDED \"AS IS\", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR"
#> [17] " IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,"
#> [18] " FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE"
#> [19] " AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER"
#> [20] " LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,"
....You can run Bayesian calibration using Stan using $stan_run().
If data argument is not specified, the output of $data_bc() is
directly used.
options(mc.cores = parallel::detectCores())
bc$stan_run(iter = 300, chains = 3)Instead of using builtin Stan model, you can also provide your own Stan
code using file argument.
bc$stan_run(file = bc$stan_file(tempfile(fileext = ".stan")), iter = 300, chains = 3)You can also use custom data set
res <- bc$stan_run(data = bc$data_bc(), iter = 300, chains = 3)
#> Warning: The largest R-hat is NA, indicating chains have not mixed.
#> Running the chains for more iterations may help. See
#> http://mc-stan.org/misc/warnings.html#r-hat
#> Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
#> Running the chains for more iterations may help. See
#> http://mc-stan.org/misc/warnings.html#bulk-ess
#> Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
#> Running the chains for more iterations may help. See
#> http://mc-stan.org/misc/warnings.html#tail-essThe stan model is store in fit
rstan::stan_trace(res$fit)
#> 'pars' not specified. Showing first 10 parameters by default.
rstan::stan_hist(res$fit)
#> 'pars' not specified. Showing first 10 parameters by default.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
The predicted values is stored in y_pred.
str(res$y_pred)
#> Classes 'data.table' and 'data.frame': 5400 obs. of 8 variables:
#> $ index : int 1 1 1 1 1 1 1 1 1 1 ...
#> $ sample : int 1 2 3 4 5 6 7 8 9 10 ...
#> $ Date/Time : chr " 07/01 06:00" " 07/01 06:00" " 07/01 06:00" " 07/01 06:00" ...
#> $ COOLSYS1 CHILLER 1:Chiller Evaporator Inlet Temperature [C](6 Hour) : num 8.19 8.19 8.19 8.19 8.19 ...
#> $ COOLSYS1 CHILLER 1:Chiller Evaporator Outlet Temperature [C](6 Hour): num 8.08 8.08 8.08 8.08 8.08 ...
#> $ COOLSYS1 CHILLER 1:Chiller Evaporator Mass Flow Rate [kg/s](6 Hour) : num 0.0683 0.0683 0.0683 0.0683 0.0683 ...
#> $ COOLSYS1 CHILLER 1:Chiller Electric Power [W](6 Hour) : num 0 0 0 0 0 0 0 0 0 0 ...
#> $ COOLSYS1 CHILLER 1:Chiller Electric Power [W](6 Hour) [Prediction] : num -65.3 56.4 -67.2 -23.1 83.5 ...
#> - attr(*, ".internal.selfref")=<externalptr>