joelnc
commented
3 years ago Hi Chao,
As mentioned on the SWMM list, there are many different ways to subset/swap values in a data frame. Below, I use within rather than transform. But assuming your upper and lower bounds correspond to slope, pct. imp., impervious n and previous n, the code below works for me with the rest of your example.
Also note the parameter names in subcatchments have underscores (Perc_Slope), while the ones in subareas have hyphens (N-Perv).
obj_fun <- function(x, inp, obs) {
## inp$subareas <- transform(
## inp$subareas,
## S_Imperv = x
## )
inp$subcatchments <- within(inp$subcatchments, {
Perc_Slope <- x[1]
Perc_Imperv <- x[2]
})
inp$subareas <- within(inp$subareas, {
"N-Imperv" <- x[3]
"N-Perv" <- x[4]
})I also made another slight modification in obj_fun to run on my system:
tmp_inp <- tempfile()
## out_dir <- "C:\Users\yeche\Documents\R_script\Tempfile"
## tmp_inp <- file.path(out_dir, "Example1_con.inp")
write_inp(inp, tmp_inp)Does this help?
CChaoE17
thwald
Hi, all,
I am working on my model calibration using the swmmr package. Firstly, I run the example script provided by Dominik on website and it works well, the Perc_Imperv parameter can be updated in new inp file during calibration. However, when I change the parameter to 'N_Imperv', 'N_Perv', 'S_Imperv', and 'S_Perv', from inp$subareas, the code can also work but the inp file did not get updated values of these parameters, rather, the value of 'Percentage Routed' changed even it is not the calibrated parameter. I also tried to calibrate the width parameter and it can get updated value for new inp file during process. I guess 'Perc_Imperv' and 'Width' are parameters from inp$subcatchments, so they can get updated value, however, the parameter from inp$subareas cannot get. I have no idea about this problem. The following are my script. If any of you happen to know what the problem is, please give me some ideas. Thank you. library(swmmr) library(DEoptim)
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 <- read_out( file = swmm_files$out, iType = 1, object_name = "18", vIndex = 4 )[["18"]]$total_inflow
read model structure
inp <- read_inp(swmm_files$inp)
show the original parameter values
inp$subcatchments inp$subareas inp$infiltration
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 <- function(x, inp, obs) {
set new parameters and update inp object
inp$subareas <- transform( inp$subareas,S_Imperv = x )
out_dir <- "C:\Users\yeche\Documents\R_script\Tempfile"
write new inp file to disk, I want to see the updated inp file, so I changed the path
tmp_inp <- file.path(out_dir, "Example1_con.inp") write_inp(inp, file.path(out_dir, "Example1_con.inp"))
run swmm with new parameter set
swmm_files <- suppressMessages(run_swmm(tmp_inp, stdout = NULL))
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 }
set.seed(1234) # to get reproducible results
calibration_res <- DEoptim( fn = obj_fun, lower = c(0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05), upper = c(0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5), control = list( itermax = 20, # maximum iterations trace = 5, # 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)