--- title: "Workflow Optimierung (Eisenstadt · Wien · Bad Aussee)" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Workflow Optimierung (Eisenstadt · Wien · Bad Aussee)} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE, eval = TRUE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 7.5, fig.height = 3.2 ) is_ghactions <- tolower(Sys.getenv("GITHUB_ACTIONS")) == "true" || tolower(Sys.getenv("CI")) %in% c("true", "1", "yes") # extdata robust aufloesen: zuerst installiertes Paket, sonst Quellbaum # (Knit direkt im Repo, auch ohne installiertes kwb.raindrop) extdata_path <- function(...) { p <- system.file("extdata", ..., package = "kwb.raindrop") if (nzchar(p)) return(p) src <- file.path("..", "inst", "extdata", ...) if (file.exists(src)) normalizePath(src) else "" } path_base <- extdata_path("models", "eisenstadt-2005", "base.h5") data_available <- nzchar(path_base) && file.exists(path_base) is_windows <- Sys.info()[["sysname"]] == "Windows" can_run <- data_available && is_windows && !is_ghactions t_vignette_start <- Sys.time() ``` ## Ziel Dieser Workflow findet für **alle drei Standorte** (Eisenstadt 2005, Wien, Bad Aussee) die **günstigste Muldenkonfiguration je Überlaufziel** x = 0…5 — pro Speichertyp (Sickerbox / Schotterrigol) — mit `optimise_swale_design()` statt eines Brute-Force-Rasters. Das Verfahren ist reine Bisektion ("Zahlenraten"): Fläche schrumpfen, dann Muldentiefe, Speicher nur erhöhen, wenn die Fläche am Anschlag klemmt. Voraussetzung ist die in der Vignette `monotonicity_analysis` belegte Monotonie ("größer = nie mehr Überläufe"); die dort abgeleiteten Absicherungen (Rand-Guard, Volumen-Schiedsrichter) sind in `find_min_feasible()` eingebaut. Die Filterdurchlässigkeit wird fest auf das Maximum gesetzt (kostenfrei dominant: gleiche Verdunstung, nie mehr Überläufe), Fläche und Muldentiefe werden stufenlos gesucht, die Sickerbox in den Default-Stufen 300/600/900/1200 mm (`sickerbox_level_presets()` bietet Hersteller-Alternativen), der Schotterrigol stufenlos im 3-fachen Box-Bereich. **Laufzeit:** Ein Engine-Lauf dauert ~2 s (Eisenstadt, 1 Jahr) bzw. ~15 s (Wien / Bad Aussee, 15-Jahres-Serien). Die 6 Tasks (Standort × Speichertyp) laufen parallel; die Gesamtdauer entspricht dem längsten Einzeltask — ca. 15–20 Minuten. ```{r availability_note, echo = FALSE, results = 'asis', eval = !can_run} cat(sprintf(paste0( "> **Hinweis:** Die Rechen-Chunks wurden übersprungen. Prüfungen: ", "Windows: %s · CI/GitHub Actions: %s · base.h5 gefunden: %s (%s). ", "Auf einem lokalen Windows-Rechner sollten alle drei Bedingungen ", "erfüllt sein — falls base.h5 fehlt: Vignette aus dem Paket-Repo ", "heraus rendern oder das Paket installieren.\n"), is_windows, is_ghactions, data_available, if (nzchar(path_base)) path_base else "weder installiert noch ../inst/extdata" )) ``` ## Standort-Konfiguration Die Standorte unterscheiden sich nur in zwei Punkten: dem Modell-Template (`base.h5`) und der Frage, ob eigene Regen-/ET0-Zeitreihen in die Engine geschrieben werden (Wien und Bad Aussee) oder die Kurven aus `base.h5` verwendet werden (Eisenstadt). Beides kapselt `make_swale_runner()`; die Zeitreihen-Aufbereitung (mm/h-Konvention, Serien-Angleichung) übernimmt `read_site_timeseries()`. ```{r site_config, eval = can_run} # 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 = "/models//input", dir_output = "/models//output", dir_target_output = "/", file_errors_hdf5 = "Fehlerprotokoll.h5", file_results_hdf5_element = "Mulde_Rigole.h5", file_results_hdf5_flaeche = "Dach.h5", file_results_hdf5_verschaltungen = "_Verschaltungen.h5", file_results_txt = "Mulde_Rigole_RAINDROP.txt", file_results_txt_multilayer = "Mulde_Rigole_RAINDROP_multi_layer.txt", file_target = ".h5", path_base = extdata_path("models", model_dir, "base.h5"), path_exe = download_engine(), path_errors_hdf5 = "/", path_results_hdf5_element = "/", path_results_hdf5_flaeche = "/", path_results_hdf5_verschaltungen = "/", path_results_txt = "/", path_results_txt_multilayer = "/", path_target_input = "/" ) } ``` ## Suchraum und Kostensätze Die Suchbereiche der variablen Parameter sind die Default-Argumente von `optimise_swale_design()` — bewusst identisch mit den min/max-Bereichen des Brute-Force-Rasters, denn nur dort ist die Monotonie geprüft. Hier stehen sie explizit, damit sie sichtbar und anpassbar sind: ```{r search_space, eval = can_run} 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" ) )) ``` Die Kostensätze sind die Defaults nach Leimgruber (2026-03-27, `default_cost_rates()`); einzelne Sätze lassen sich hier überschreiben — gerechnet wird zunächst mit den Defaults: ```{r cost_rates, eval = can_run} 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)" ) ``` ## Optimierung aller Standorte (parallel) Parallelisiert wird über **Standort × Speichertyp = 6 unabhängige Tasks**: Die Bisektion innerhalb einer Suche ist prinzipbedingt sequentiell, und die x-Ziele eines Speichertyps teilen sich den Evaluations-Cache — aber Box- und Rigol-Läufe teilen keinen einzigen Engine-Lauf (der Cache-Schlüssel enthält den Typ). Jeder Worker ist ein eigener R-Prozess mit eigenem `tempdir()`, Kollisionen sind damit ausgeschlossen. Die Wall-Time entspricht dem längsten Einzeltask (Wien/Rigol, ~15–20 min) statt der Summe aller Tasks (~65 min). Liegen die Rasterergebnisse der Workflow-Vignetten als CSV neben dieser Vignette, verengen sie als Warmstart die erste Flächensuche auf einen 25-m²-Rasterschritt. ```{r optimise_all, eval = can_run} 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() ``` ## Ergebnis: günstigstes Design je Standort und Überlaufziel ```{r results_table, eval = can_run} 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) ) ``` ```{r cost_curves, eval = can_run, fig.height = 3.4} 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() ``` ```{r export, eval = can_run} 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") ) ``` ```{r runtime_total, echo = FALSE, results = 'asis', eval = can_run} opt_wall_min <- as.numeric(difftime(t_opt_end, t_opt_start, units = "mins")) cat(sprintf(paste0( "**Gesamtlaufzeit Optimierung:** %d Engine-Läufe · Summe der ", "Task-Zeiten %.1f min · tatsächliche Laufzeit %.1f min ", "(paralleler Speedup %.1f×).\n"), sum(task_stats$n_runs_task), sum(task_stats$minutes_task), opt_wall_min, sum(task_stats$minutes_task) / max(opt_wall_min, 0.1) )) ``` ## Monte-Carlo-Analyse: Wie robust ist die Suche selbst? Die Bisektion ist deterministisch: gleiche Eingaben, gleicher Pfad, gleiches Ergebnis. Die Monte-Carlo-Frage lautet deshalb: **Hängt das gefundene Optimum vom Suchpfad ab?** Dazu wird der Teilungspunkt jeder Bisektion zufällig verschoben (`split_jitter = 0.3`: Teilung zufällig zwischen 20 % und 80 % des Intervalls statt exakt mittig) und die Optimierung `n_mc`-mal mit unterschiedlichen Seeds wiederholt — **Regen, Kostensätze und alle übrigen Eingaben bleiben unverändert** (x = 1, ohne Warmstart, damit jede Wiederholung den vollen Suchraum durchläuft). Gerechnet wird der volle Pool **3 Standorte × 2 Speichertypen × `n_mc` Wiederholungen = 60 unabhängige Tasks**, alle parallel (je Task eine komplette Neu-Optimierung mit ~15 Engine-Läufen; Wall-Time je nach Kernzahl ~20–35 min). Erwartung, wenn das Optimum eine Eigenschaft des Problems ist — und nicht des Wegs, den die Suche genommen hat: **identische Speicherstufe, Flächen-Spanne ≤ 2 × `area_tol`, Kostenspanne von wenigen Prozent**. Die Muldentiefe darf etwas weiter streuen als 2 × `height_tol`: Sie ist über die Hydraulik an die gefundene Fläche gekoppelt (eine um `area_tol` größere Fläche erlaubt eine entsprechend geringere Tiefe), sodass sich dort die Toleranzen beider Suchen addieren. ```{r mc_config, eval = can_run} n_mc <- 10 mc_seeds <- 1:n_mc ``` ```{r mc_search, eval = can_run} 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() ``` ```{r mc_search_summary, eval = can_run, fig.height = 3.4} 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") ) ``` ```{r mc_runtime, echo = FALSE, results = 'asis', eval = can_run} cat(sprintf(paste0( "**Laufzeit Such-Monte-Carlo:** %d Engine-Läufe in %.1f min ", "(parallel über %d Tasks: %d Standorte × %d Speichertypen × ", "%d Wiederholungen).\n"), sum(mc_search$n_runs_rep), as.numeric(difftime(t_mc_end, t_mc_start, units = "mins")), nrow(mc_tasks), length(sites), length(storage_spec), n_mc )) ``` ## Einordnung - **Plausibilität:** Die Optima müssen auf oder knapp unter den günstigsten zulässigen Rasterzellen liegen (Warmstart-Tabelle der Monotonie-Analyse); die Kostenkurven müssen monoton fallen. `monotonicity_warning = TRUE` in einer Zeile hieße: Zähler *und* Volumen sind gemeinsam gestiegen — dann Branch prüfen. - **Wien x = 0** ist der Stresstest: Im Raster war das Ziel mit minimalem Speicher teils unerreichbar — hier greift die Speicher-Eskalation automatisch; "infeasible_within_bounds" wäre ein reguläres Ergebnis, kein Fehler. - **Bad Aussee x = 1** trägt den bekannten +1-Zählwobble am Rasterrand; der Rand-Guard in `find_min_feasible()` deckt ihn ab. - **Nebenbedingung Tiefe:** `max_total_depth` (mm) begrenzt Muldentiefe + Filter + Speicher analytisch (DWA-A 138 / Überdeckung), ohne zusätzliche Läufe. ```{r vignette_runtime, echo = FALSE, results = 'asis', eval = can_run} cat(sprintf(paste0( "---\n\n**Gesamtlaufzeit dieser Vignette:** %.1f Minuten ", "(Optimierung %.1f min · Such-Monte-Carlo %.1f min · ", "Rest: Setup und Rendern).\n"), as.numeric(difftime(Sys.time(), t_vignette_start, units = "mins")), opt_wall_min, as.numeric(difftime(t_mc_end, t_mc_start, units = "mins")) )) ```