vignettes/How_to_autocalibrate_a_SWMM_model_with_swmmr.R

In dleutnant/swmmr: R Interface for US EPA's SWMM

## ----setup, include = FALSE---------------------------------------------------
calibration_res <- NULL
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>", 
  # https://community.rstudio.com/t/precompiling-vignette-with-devtools/1583/6
  eval = nzchar(Sys.getenv("BUILD_SWMM_VIGNETTE")) # to prevent eval on CRAN!
)

## ----load_libs----------------------------------------------------------------
#  library(swmmr)
#  library(DEoptim)

## ----model_setup--------------------------------------------------------------
#  # set path to inp
#  # If your operating system is Windows, the Example1.inp model is usually
#  # located at "C:\Users\your user name\Documents\EPA SWMM Projects\Examples".
#  # For convenience the Example1.inp model is also included in the swmmr package.
#  # Feel free to change this to your path of choice.
#  inp_file <- system.file("extdata", "Example1.inp", package = "swmmr", mustWork = TRUE)
#  
#  # both rpt and out files are temporary files
#  tmp_rpt_file <- tempfile()
#  tmp_out_file <- tempfile()
#  
#  # initiate the simulation
#  swmm_files <- run_swmm(
#    inp = inp_file,
#    rpt = tmp_rpt_file,
#    out = tmp_out_file
#  )

## ----obs----------------------------------------------------------------------
#  obs <- read_out(
#    file = swmm_files$out,
#    iType = 1,
#    object_name = "18",
#    vIndex = 4
#  )[["18"]]$total_inflow

## ----sim_and_read-------------------------------------------------------------
#  # read model structure
#  inp <- read_inp(swmm_files$inp)
#  
#  # show the original parameter values
#  inp$subcatchments[inp$subcatchments$Area > 10, ]
#  

## ----gof----------------------------------------------------------------------
#  # function calculates the goodness of fit value
#  # input x is a two column xts object, col1: obs, col2: sim
#  nse <- function(x) {
#    1 - sum((x[, 1] - x[, 2]) ^ 2) / sum((x[, 1] - mean(x[, 1])) ^ 2)
#  }

## ----obj_fun------------------------------------------------------------------
#  obj_fun <- function(x, inp, obs) {
#  
#    # set new parameters and update inp object
#    inp$subcatchments <- transform(
#      inp$subcatchments,
#      Perc_Imperv = ifelse(Area > 10, x, Perc_Imperv)
#    )
#  
#    # write new inp file to disk
#    tmp_inp <- tempfile()
#    write_inp(inp, tmp_inp)
#  
#    # run swmm with new parameter set
#    swmm_files <- suppressMessages(run_swmm(tmp_inp, stdout = NULL))
#  
#    # remove files when function exits to avoid heavy disk usage
#    on.exit(file.remove(unlist(swmm_files)))
#  
#    # read sim result
#    sim <- read_out(
#      file = swmm_files$out, # path to out file
#      iType = 1, # type: node
#      object_name = "18", # name of node
#      vIndex = 4 # parameter at node: total inflow
#    )[["18"]]$total_inflow # directly access to xts object
#  
#    # calculate goodness-of-fit
#    # note: multiply by minus one to have a real min problem (nse: +1 to -Inf)
#    nse(merge(obs, sim)) * -1
#  }
#  

## ----optim--------------------------------------------------------------------
#    set.seed(84) # to get reproducible results
#  
#    calibration_res <- DEoptim(
#      fn = obj_fun,
#      lower = c(0, 0),
#      upper = c(100, 100),
#      control = list(
#        itermax = 50, # maximum iterations
#        trace = 10, # print progress every 10th iteration
#        packages = c("swmmr"), # export packages to optimization environment
#        parVar = c("nse"), # export function to optimization environment
#        parallelType = 0 # set to 1 to use all available cores
#      ),
#      inp = inp, # 'inp' object
#      obs = obs # xts object containing observation data
#    )
#  
#    summary(calibration_res)