Workflow Optimierung (Eisenstadt · Wien · Bad Aussee)

# Im Quellbaum die Entwicklungsversion laden (immer aktuell, auch wenn
# das installierte Paket aelter ist); ausserhalb: installiertes Paket.
if (file.exists("../DESCRIPTION") &&
    requireNamespace("pkgload", quietly = TRUE)) {
  pkgload::load_all("..", quiet = TRUE)
} else {
  library(kwb.raindrop)
}

sites <- list(
 Eisenstadt_2005 = list(dir = "eisenstadt-2005", timeseries = FALSE,
   prior = "simulation_results_optimisation_Eisenstadt_2005.csv"),
 Wien = list(dir = "wien", timeseries = TRUE,
   prior = "simulation_results_optimisation_Wien.csv"),
 BadAussee = list(dir = "badaussee", timeseries = TRUE,
   prior = "simulation_results_optimisation_BadAussee.csv")
)

fixed <- list(connected_area = 1000,
              filter_height = 300,
              filter_hydraulicconductivity = 360,  # Rastermaximum, gratis
              bottom_hydraulicconductivity = 12)

make_path_list <- function(modelname, model_dir) {
  list(
    modelname = modelname,
    root_path = file.path(tempdir(), paste0("raindrop_opt_", model_dir)),
    dir_input  = "<root_path>/models/<modelname>/input",
    dir_output = "<root_path>/models/<modelname>/output",
    dir_target_output = "<dir_output>/<dir_target>",
    file_errors_hdf5 = "Fehlerprotokoll.h5",
    file_results_hdf5_element = "Mulde_Rigole.h5",
    file_results_hdf5_flaeche = "Dach.h5",
    file_results_hdf5_verschaltungen = "<dir_target>_Verschaltungen.h5",
    file_results_txt = "Mulde_Rigole_RAINDROP.txt",
    file_results_txt_multilayer = "Mulde_Rigole_RAINDROP_multi_layer.txt",
    file_target = "<dir_target>.h5",
    path_base = extdata_path("models", model_dir, "base.h5"),
    path_exe  = download_engine(),
    path_errors_hdf5 = "<dir_target_output>/<file_errors_hdf5>",
    path_results_hdf5_element = "<dir_target_output>/<file_results_hdf5_element>",
    path_results_hdf5_flaeche = "<dir_target_output>/<file_results_hdf5_flaeche>",
    path_results_hdf5_verschaltungen = "<dir_target_output>/<file_results_hdf5_verschaltungen>",
    path_results_txt = "<dir_target_output>/<file_results_txt>",
    path_results_txt_multilayer = "<dir_target_output>/<file_results_txt_multilayer>",
    path_target_input = "<dir_input>/<file_target>"
  )
}

area_bounds   <- c(25, 200)   # Muldenflaeche [m2], stufenlos
area_tol      <- 2            # Aufloesung der Flaechensuche [m2]
height_bounds <- c(100, 300)  # Muldentiefe [mm], stufenlos
height_tol    <- 10           # Aufloesung der Tiefensuche [mm]
storage_spec  <- default_storage_spec()
# Alternativ Hersteller-Stufen, z. B.:
# storage_spec <- default_storage_spec(
#   levels = sickerbox_level_presets()$graf_ecobloc_smart)

knitr::kable(tibble::tibble(
  Parameter = c("mulde_area [m2]", "mulde_height [mm]",
                "storage_height Sickerbox [mm]",
                "storage_height Schotterrigol [mm]",
                "filter_hydraulicconductivity [mm/h]",
                "filter_height [mm]"),
  Suchraum = c(
    sprintf("stufenlos %g bis %g (Toleranz %g)",
            area_bounds[1], area_bounds[2], area_tol),
    sprintf("stufenlos %g bis %g (Toleranz %g)",
            height_bounds[1], height_bounds[2], height_tol),
    paste("Stufen:",
          paste(storage_spec$infiltration_box$levels, collapse = " / ")),
    sprintf("stufenlos %g bis %g (Toleranz %g)",
            storage_spec$gravel_trench$bounds[1],
            storage_spec$gravel_trench$bounds[2],
            storage_spec$gravel_trench$tol),
    "fix = 360 (Rastermaximum; kostenfrei dominant)",
    "fix = 300"
  )
))

cost_rates <- default_cost_rates()
# Beispiel fuer eine Anpassung (auskommentiert):
# cost_rates$excavation_eur_per_m3       <- 85
# cost_rates$infiltration_box_eur_per_m3 <- 400

knitr::kable(
  tibble::tibble(Kostensatz = names(cost_rates),
                 `EUR je m2/m3` = unlist(cost_rates)),
  caption = "Kostensaetze (inkl. Einbau)"
)

t_opt_start <- Sys.time()

tasks <- expand.grid(site = names(sites), type = names(storage_spec),
                     stringsAsFactors = FALSE)

future::plan(future::multisession,
             workers = min(nrow(tasks),
                           max(1, parallel::detectCores() - 1)))

opt_list <- future.apply::future_lapply(seq_len(nrow(tasks)), function(i) {
  site <- tasks$site[i]
  type <- tasks$type[i]

  # Worker = eigener Prozess: Paket dort genauso laden wie im Hauptprozess
  if (file.exists("../DESCRIPTION") &&
      requireNamespace("pkgload", quietly = TRUE)) {
    pkgload::load_all("..", quiet = TRUE)
  } else {
    library(kwb.raindrop)
  }

  cfg <- sites[[site]]
  ts <- if (cfg$timeseries) {
    read_site_timeseries(
      extdata_path("models", cfg$dir, "rain.csv.gz"),
      extdata_path("models", cfg$dir, "et.csv"),
      verbose = FALSE
    )
  } else {
    NULL
  }

  run_fn <- make_swale_runner(make_path_list(site, cfg$dir),
                              timeseries_rain = ts$rain,
                              timeseries_et = ts$et)

  prior <- if (file.exists(cfg$prior)) {
    readr::read_csv(cfg$prior, show_col_types = FALSE)
  } else {
    NULL
  }

  t0 <- Sys.time()
  opt <- optimise_swale_design(run_fn, x_targets = 0:5,
                               area_bounds = area_bounds,
                               area_tol = area_tol,
                               height_bounds = height_bounds,
                               height_tol = height_tol,
                               storage_spec = storage_spec[type],
                               fixed = fixed,
                               prior_results = prior,
                               cost_rates = cost_rates,
                               verbose = FALSE)
  opt$site <- site
  opt$n_runs_task <- attr(opt, "n_runs_total")
  opt$minutes_task <- round(as.numeric(
    difftime(Sys.time(), t0, units = "mins")), 1)
  opt
}, future.seed = TRUE)

future::plan(future::sequential)
opt_all <- dplyr::bind_rows(opt_list)

t_opt_end <- Sys.time()

knitr::kable(
  opt_all[, c("site", "x", "storage_type", "status", "mulde_area",
              "mulde_height", "storage_height", "n_overflows",
              "overflow_volume_m3", "et_pct", "cost_total")],
  digits = c(NA, 0, NA, NA, 1, 0, 0, 0, 1, 1, 0)
)

library(ggplot2)

ok <- opt_all[opt_all$status == "ok", ]
ggplot(ok, aes(x, cost_total / 1000, colour = storage_type)) +
  geom_line() +
  geom_point(size = 2) +
  facet_wrap(~ site, scales = "free_y") +
  scale_x_continuous(breaks = 0:5) +
  labs(title = "Kosten-Wirksamkeits-Kurven: Was kostet ein Ueberlauf weniger?",
       x = "Ueberlaufziel x (zulaessige Ereignisse)",
       y = "Kosten Optimum [Tsd. EUR]",
       colour = "Speichertyp",
       caption = cost_rates_caption("de", cost_rates)) +
  theme_bw()

readr::write_csv(opt_all, "optimisation_results_all_sites.csv")
for (site in unique(opt_all$site)) {
  readr::write_csv(opt_all[opt_all$site == site, ],
                   sprintf("optimisation_results_%s.csv", site))
}

task_stats <- unique(opt_all[, c("site", "storage_type", "n_runs_task",
                                 "minutes_task")])
knitr::kable(
  task_stats,
  col.names = c("Standort", "Speichertyp", "Engine-Laeufe", "Minuten")
)

n_mc <- 10
mc_seeds <- 1:n_mc

t_mc_start <- Sys.time()

mc_tasks <- expand.grid(site = names(sites), type = names(storage_spec),
                        rep = seq_len(n_mc), stringsAsFactors = FALSE)

future::plan(future::multisession,
             workers = min(nrow(mc_tasks),
                           max(1, parallel::detectCores() - 1)))

mc_search <- future.apply::future_lapply(seq_len(nrow(mc_tasks)), function(i) {
  site <- mc_tasks$site[i]
  type <- mc_tasks$type[i]
  rep  <- mc_tasks$rep[i]

  if (file.exists("../DESCRIPTION") &&
      requireNamespace("pkgload", quietly = TRUE)) {
    pkgload::load_all("..", quiet = TRUE)
  } else {
    library(kwb.raindrop)
  }
  set.seed(mc_seeds[rep])

  cfg <- sites[[site]]
  ts <- if (cfg$timeseries) {
    read_site_timeseries(
      extdata_path("models", cfg$dir, "rain.csv.gz"),
      extdata_path("models", cfg$dir, "et.csv"),
      verbose = FALSE
    )
  } else {
    NULL
  }
  run_fn <- make_swale_runner(make_path_list(paste0(site, "_MC"), cfg$dir),
                              timeseries_rain = ts$rain,
                              timeseries_et = ts$et)

  opt <- optimise_swale_design(
    run_fn, x_targets = 1,
    area_bounds = area_bounds, area_tol = area_tol,
    height_bounds = height_bounds, height_tol = height_tol,
    storage_spec = storage_spec[type],
    fixed = fixed,
    split_jitter = 0.3,
    cost_rates = cost_rates, verbose = FALSE
  )
  opt$site <- site
  opt$rep <- rep
  opt$n_runs_rep <- attr(opt, "n_runs_total")
  opt
}, future.seed = TRUE)

future::plan(future::sequential)
mc_search <- dplyr::bind_rows(mc_search)

t_mc_end <- Sys.time()

ok_mc <- mc_search[mc_search$status == "ok", ]

ggplot(ok_mc, aes(rep, cost_total / 1000, colour = storage_type)) +
  geom_point(size = 2) +
  facet_wrap(~ site, scales = "free_y") +
  scale_x_continuous(breaks = seq_len(n_mc)) +
  labs(title = sprintf(
         "Such-Monte-Carlo (x = 1, %d zufaellige Suchpfade je Standort und Typ)",
         n_mc),
       x = "Wiederholung (Seed)",
       y = "Kosten Optimum [Tsd. EUR]",
       colour = "Speichertyp") +
  theme_bw()

knitr::kable(
  ok_mc %>%
    dplyr::group_by(site, storage_type) %>%
    dplyr::summarise(
      flaeche_spanne_m2 = max(mulde_area) - min(mulde_area),
      tiefe_spanne_mm = max(mulde_height) - min(mulde_height),
      speicher_identisch = dplyr::n_distinct(storage_height) == 1,
      kosten_min = min(cost_total),
      kosten_max = max(cost_total),
      kosten_spanne_pct = round(100 * (max(cost_total) - min(cost_total)) /
                                  min(cost_total), 2),
      .groups = "drop"
    ),
  digits = 1,
  caption = paste("Streuung ueber die Suchpfade -- erwartet: identische",
                  "Speicherstufe, Flaechen-Spanne <= 2 x area_tol,",
                  "Kostenspanne wenige Prozent; die Tiefe streut wegen",
                  "der Flaechen-Kopplung weiter")
)