Workflow Optimierung — Simultane Suche (alle Parameter gleichzeitig) • kwb.raindrop

Ziel

Die Vignette workflow_optimisation findet die günstigste Muldenkonfiguration je Überlaufziel per Bisektion: Parameter nacheinander, gestützt auf die Monotonie je Parameter (Monotonie-Analyse). Diese Vignette ist die unabhängige Gegenprobe: optimise_swale_design_simultaneous() optimiert alle Parameter gleichzeitig — Fläche, Muldentiefe und Speicherhöhe in einem Zug — und kommt dabei ohne die Monotonie-Annahme aus.

Warum das funktioniert: Unzulässige Designs (n_overflows > x) werden nicht ausgeschlossen, sondern bestraft — jedes unzulässige Design ist teurer als jedes zulässige, und überzählige Überlaufereignisse staffeln die Strafe. Das Optimum liegt (Kosten steigen monoton mit jedem Parameter) genau auf der Zulässigkeitsgrenze; die Straffunktion erlaubt der Suche, diese Grenze zu überqueren, auf beiden Seiten Information zu sammeln und die Parameter in einem einzigen Schritt gegeneinander zu tauschen (z. B. weniger Fläche gegen mehr Speicher). Ein Ausschluss unzulässiger Punkte würde der Suche jenseits der Grenze jede Richtungsinformation nehmen — sie würde blind an der Grenze abprallen, statt an ihr entlangzuwandern.

Drei Suchverfahren teilen sich dieselbe Infrastruktur (Straf-Ziel, Evaluations-Cache, Toleranz-Rasterung, Multi-Tal-Feinschliff) und unterscheiden sich nur darin, wie sie Kandidaten vorschlagen (method-Argument):

nelder_mead (Default, Empfehlung): Multistart-Simplex (stats::optim()) — Warmstart, Optimum des vorigen x-Ziels, je ein Anker-Start pro Speicherstufe, Raumfüller; jeder Start erhält einen gleichen Anteil am Laufbudget (max_evals).
diff_evolution: kompakte Differential Evolution (DE/rand/1/bin) — Vergleichsverfahren; deterministisch über einen internen Park-Miller-Generator (seed-Argument), Rs globaler Zufallsstrom (.Random.seed) bleibt unangetastet.
halton_search: quasi-zufällige, raumfüllende Halton-Stichprobe — bewusst naive Baseline, die zeigt, was die strukturierten Verfahren schlagen müssen.

Laufzeit: Die simultane Suche braucht je (Speichertyp, x)-Zelle mehr Engine-Läufe als die Bisektion (typisch 60–120 statt ~15; Suchphase plus Feinschliff). Über 6 Überlaufziele summiert sich das auf ~400–700 Engine-Läufe je Task; bei ~15 s je Lauf (Wien / Bad Aussee) sind das 2–3 h je Task — und wenn weniger freie Kerne als 6 Tasks verfügbar sind, entsprechend mehr Wandzeit (realistisch 2–5 h für den Nelder-Mead-Sweep). Eisenstadt (~2 s je Lauf) bleibt bei Minuten. Beide Rechen-Chunks zeigen deshalb einen Live-Fortschrittsbalken (ein Tick je Engine-Lauf, über progressr aus den Worker-Prozessen heraus) — ein über Minuten stehender Balken wäre ein echter Hänger, ein langsam wandernder ist Normalbetrieb. Der Methodenvergleich am Ende rechnet nur eine Zelle (x = 1) je Standort und Speichertyp. Wer zuerst einen schnellen Funktionstest will, rechnet nur Eisenstadt (Kommentar im Chunk site_config); der Laufzeit-Hebel für den vollen Sweep ist max_evals (Kommentar im selben Chunk).

Plattenplatz: Jeder Engine-Lauf legt ein eigenes Szenario an (Kopie der base.h5 plus Output-HDF5s). make_swale_runner() löscht diese Dateien standardmäßig direkt nach dem Einlesen der dünnen Ergebniszeile (cleanup = TRUE) — ohne dieses Aufräumen füllen mehrere hundert Läufe je Task das Temp-Laufwerk und die Engine bricht mit No space left on device ab. Vor einem Neustart nach einem solchen Abbruch die raindrop_sim_*-Verzeichnisse unter tempdir() bzw. %LOCALAPPDATA%\Temp (Rtmp*) löschen.

Hinweis: Die Rechen-Chunks wurden übersprungen. Prüfungen: Windows: TRUE · CI/GitHub Actions: TRUE · base.h5 gefunden: TRUE (D:/a/_temp/Library/kwb.raindrop/extdata/models/eisenstadt-2005/base.h5). 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.

Standort-Konfiguration

Identisch mit der Vignette workflow_optimisation (gleiche Standorte, gleiche Suchräume, gleiche Kostensätze — nur so ist der Vergleich aussagekräftig):

# 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")
)
# Schneller Funktionstest (~15-30 min statt Stunden): nur Eisenstadt --
# 1-Jahres-Modell, ~2 s je Engine-Lauf. Zusaetzlich beschleunigt ein
# reduziertes Budget (z. B. max_evals <- 50) und, auf Windows deutlich,
# eine Virenscanner-Ausnahme fuer %LOCALAPPDATA%\Temp (jede Szenario-
# Datei und jeder Engine-Start wird sonst einzeln gescannt):
sites <- sites["Eisenstadt_2005"]

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

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()
cost_rates    <- default_cost_rates()
max_evals     <- 80           # Laufzeit-Hebel: Budget an frischen
                              # Engine-Laeufen je Zelle (Suchphase;
                              # der Feinschliff kommt obendrauf)

make_path_list <- function(modelname, model_dir) {
  list(
    modelname = modelname,
    root_path = file.path(tempdir(), paste0("raindrop_sim_", 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>"
  )
}

# Ein Task = eine komplette Optimierung (Standort x Speichertyp x
# Methode); gekapselt, damit Haupt-Sweep und Methodenvergleich denselben
# Code nutzen. tick/tick_cap: progressr-Fortschritt ueber die
# Worker-Grenze hinweg -- ein Tick je Engine-Lauf, am Task-Ende wird
# der Rest des Task-Kontingents aufgefuellt, damit der Balken exakt
# bei 100 % endet.
run_simultaneous_task <- function(site, type, method, x_targets,
                                  model_suffix,
                                  tick = NULL, tick_cap = Inf) {
  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(paste0(site, "_", model_suffix), cfg$dir),
    timeseries_rain = ts$rain,
    timeseries_et = ts$et
  )

  ticks_sent <- 0L
  if (!is.null(tick)) {
    inner_fn <- run_fn
    run_fn <- function(params) {
      if (ticks_sent < tick_cap) {
        ticks_sent <<- ticks_sent + 1L
        tick(sprintf("%s | %s | Lauf %d", site, type, ticks_sent))
      }
      inner_fn(params)
    }
  }

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

  t0 <- Sys.time()
  opt <- optimise_swale_design_simultaneous(
    run_fn, x_targets = x_targets,
    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,
    method = method,
    max_evals = max_evals,
    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)
  if (!is.null(tick) && is.finite(tick_cap) && tick_cap > ticks_sent) {
    tick(sprintf("%s | %s fertig (%d Laeufe, %.1f min)",
                 site, type, opt$n_runs_task[1], opt$minutes_task[1]),
         amount = tick_cap - ticks_sent)
  }
  opt
}

Simultane Optimierung aller Standorte (Nelder-Mead, parallel)

Wie im Bisektions-Workflow laufen Standort × Speichertyp = 6 unabhängige Tasks parallel; innerhalb eines Tasks teilen sich die x-Ziele den Evaluations-Cache.

t_nm_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)))

# Live-Fortschritt ueber die Worker hinweg: 1 Tick = 1 Engine-Lauf
# (tick_cap = grosszuegiges Kontingent je Task; der Rest wird am
# Task-Ende aufgefuellt, der Balken endet also exakt bei 100 %)
progressr::handlers(progressr::handler_txtprogressbar())
tick_cap_nm <- 1000

nm_all <- progressr::with_progress({
  p <- progressr::progressor(steps = nrow(tasks) * tick_cap_nm)
  nm_list <- future.apply::future_lapply(seq_len(nrow(tasks)), function(i) {
    run_simultaneous_task(tasks$site[i], tasks$type[i],
                          method = "nelder_mead", x_targets = 0:5,
                          model_suffix = "NM",
                          tick = p, tick_cap = tick_cap_nm)
  }, future.seed = TRUE)
  dplyr::bind_rows(nm_list)
})

future::plan(future::sequential)

t_nm_end <- Sys.time()

Ergebnis: günstigstes Design je Standort und Überlaufziel

knitr::kable(
  nm_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 <- nm_all[nm_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 (simultane Suche, Nelder-Mead)",
       x = "Ueberlaufziel x (zulaessige Ereignisse)",
       y = "Kosten Optimum [Tsd. EUR]",
       colour = "Speichertyp",
       caption = cost_rates_caption("de", cost_rates)) +
  theme_bw()

readr::write_csv(nm_all, "optimisation_results_simultaneous_all_sites.csv")

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

Gegenprobe: Vergleich mit der Bisektion

Liegt der Export der Bisektions-Vignette (optimisation_results_all_sites.csv) neben dieser Vignette, werden beide Optima Zelle für Zelle verglichen. Beide Verfahren müssen — bis auf die Suchtoleranzen, also wenige Prozent — auf dieselben Kosten kommen. Fände die simultane Suche systematisch günstigere Designs, wäre das ein Hinweis auf Parameter-Wechselwirkungen, die die Koordinatensuche nicht sieht (und ein Fall für die Monotonie-Analyse); fände sie nur teurere, hat der Simplex sein Laufbudget nicht ausgeschöpft oder klemmt in einem lokalen Tal (n_starts / max_evals erhöhen).

bisect_all <- readr::read_csv("optimisation_results_all_sites.csv",
                              show_col_types = FALSE)

vergleich <- dplyr::full_join(
  dplyr::select(bisect_all, site, storage_type, x,
                status_bisektion = status, kosten_bisektion = cost_total),
  dplyr::select(nm_all, site, storage_type, x,
                status_simultan = status, kosten_simultan = cost_total),
  by = c("site", "storage_type", "x")
) %>%
  dplyr::mutate(
    delta_pct = round(100 * (kosten_simultan - kosten_bisektion) /
                        kosten_bisektion, 1)
  ) %>%
  dplyr::arrange(site, storage_type, x)

knitr::kable(
  vergleich, digits = 0,
  caption = paste("Gegenprobe Bisektion vs. simultane Suche:",
                  "delta_pct < 0 heisst, die simultane Suche hat ein",
                  "guenstigeres Design gefunden")
)

Alternative Optimierer im Vergleich

Dieselbe Zelle (x = 1, beide Speichertypen, alle Standorte), drei Suchverfahren: Nelder-Mead (aus dem Haupt-Sweep oben), Differential Evolution und die Halton-Baseline. Erwartung: Nelder-Mead und DE liegen innerhalb weniger Prozent beieinander; die naive Halton-Stichprobe bleibt trotz des gemeinsamen Feinschliffs messbar dahinter — der Abstand zeigt, wie viel die strukturierte Suche beiträgt. Parallelisiert über Standort × Speichertyp × Methode = 12 Tasks.

t_cmp_start <- Sys.time()

method_tasks <- expand.grid(site = names(sites),
                            type = names(storage_spec),
                            method = c("diff_evolution", "halton_search"),
                            stringsAsFactors = FALSE)

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

tick_cap_cmp <- 300   # eine Zelle je Task

cmp_all <- progressr::with_progress({
  p <- progressr::progressor(steps = nrow(method_tasks) * tick_cap_cmp)
  cmp_list <- future.apply::future_lapply(seq_len(nrow(method_tasks)),
                                          function(i) {
    run_simultaneous_task(method_tasks$site[i], method_tasks$type[i],
                          method = method_tasks$method[i], x_targets = 1,
                          model_suffix = toupper(substr(
                            method_tasks$method[i], 1, 2)),
                          tick = p, tick_cap = tick_cap_cmp)
  }, future.seed = TRUE)
  dplyr::bind_rows(cmp_list)
})

future::plan(future::sequential)

methoden <- dplyr::bind_rows(
  dplyr::filter(nm_all, x == 1),
  cmp_all
)

t_cmp_end <- Sys.time()

methoden_breit <- methoden %>%
  dplyr::select(site, storage_type, method, cost_total, n_runs_new) %>%
  tidyr::pivot_wider(names_from = method,
                     values_from = c(cost_total, n_runs_new))

knitr::kable(
  methoden_breit, digits = 0,
  caption = paste("Methodenvergleich bei x = 1: Kosten des Optimums und",
                  "frische Engine-Laeufe der Zelle je Suchverfahren.",
                  "Der Nelder-Mead-Wert stammt aus dem Haupt-Sweep",
                  "(profitiert dort vom Cache der uebrigen x-Ziele).")
)

ggplot(methoden[methoden$status == "ok", ],
       aes(method, cost_total / 1000, fill = storage_type)) +
  geom_col(position = "dodge") +
  facet_wrap(~ site, scales = "free_y") +
  labs(title = "Kosten des gefundenen Optimums je Suchverfahren (x = 1)",
       x = "Suchverfahren",
       y = "Kosten Optimum [Tsd. EUR]",
       fill = "Speichertyp",
       caption = cost_rates_caption("de", cost_rates)) +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 20, hjust = 1))

Einordnung

Konsistenz: Bisektion und simultane Suche bestätigen sich gegenseitig, wenn ihre Kosten je Zelle nur um wenige Prozent differieren — dann ist das Optimum eine Eigenschaft des Problems, nicht des Suchwegs.
Wenn die simultane Suche systematisch günstiger ist, gibt es zwei mögliche Ursachen: (a) die Monotonie-Annahme der Bisektion ist verletzt (echter Modell-Alarm → die Monotonie-Analyse prüfen), oder
1. der Spezifikkosten-Proxy, aus dem die Bisektion ihre Suchreihenfolge herleitet (€ je mm Speicherkapazität, Kapazitätsmodell V ≈ Fläche × (Tiefe + Porosität × Speicher)), greift zu kurz — etwa weil ein Hebel nicht kapazitäts-additiv wirkt (z. B. eine variable Filterdurchlässigkeit, die die Hydraulik nichtlinear verändert) oder das Kapazitätsmodell die Standort-Hydraulik schlecht beschreibt. Die simultane Suche trägt die cost_rates direkt in ihrer Zielfunktion und braucht weder Proxy noch Kapazitätsmodell: Sie ist die annahmefreie Instanz und das primäre Verfahren, sobald Parameter ohne saubere Grenzkosten-je-Kapazität ins Spiel kommen.
monotonicity_warning = TRUE heißt hier: Unter den Evaluationen der Zelle liegt ein strikt größeres Design mit mehr Überläufen und mehr Überlaufvolumen — echte Nicht-Monotonie; dann verdient die Zelle einen Blick in das "evaluations"-Attribut.
Methodenwahl: nelder_mead bleibt die Empfehlung (beste Präzision je Engine-Lauf). diff_evolution ist die Absicherung gegen Simplex-Artefakte, halton_search die Messlatte von unten.
Budget: max_evals begrenzt die Suchphase je Zelle; der Feinschliff (Multi-Tal-Musterabstieg) kommt obendrauf. Wer Laufzeit sparen muss, reduziert zuerst x_targets, dann max_evals.