diff --git a/.Rbuildignore b/.Rbuildignore index aaee986..b2336fa 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -9,5 +9,6 @@ ^index\.md$ ^README\.md$ ^vignettes/index\.Rmd$ +^vignettes/monotonicity_analysis$ ^\.positai$ ^\.claude$ diff --git a/.gitignore b/.gitignore index 5c9c9ee..4024089 100644 --- a/.gitignore +++ b/.gitignore @@ -5,3 +5,17 @@ docs inst/doc .positai + +# Von den Vignetten erzeugte Ergebnisdateien (Workflows: simulation_results_*; +# monotonicity_analysis.Rmd: mono_*). Reproduzierbar durch erneutes Rendern +# -> nicht einchecken. HTML- und R-Dateien unter vignettes/ ignoriert +# bereits vignettes/.gitignore (*.html, *.R). +vignettes/*.csv +vignettes/*.pdf +vignettes/figure/ +Rplots.pdf + +# Deploy-Ordner der Monotonie-Analyse: komplett unversioniert. Achtung: +# Das betrifft auch den handgeschriebenen Ergebnisbericht index.html +# (liegt lokal, als claude.ai-Artefakt und deployt auf dem Server). +vignettes/monotonicity_analysis/ diff --git a/DESCRIPTION b/DESCRIPTION index 56be511..3f6e263 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -41,9 +41,11 @@ Suggests: plotly, readr, rmarkdown, + testthat (>= 3.0.0), writexl -VignetteBuilder: +VignetteBuilder: knitr +Config/testthat/edition: 3 Remotes: github::kwb-r/kwb.event, github::kwb-r/kwb.utils diff --git a/NAMESPACE b/NAMESPACE index 6e8323e..08ebe14 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -7,7 +7,10 @@ export(cost_rates_caption) export(default_canonical_wb_variables) export(default_cost_rates) export(default_param_labels) +export(default_storage_spec) +export(default_storage_types) export(download_engine) +export(find_min_feasible) export(find_single_param_variations) export(get_simulation_results_all) export(get_simulation_results_optim) @@ -18,6 +21,8 @@ export(h5_read_values) export(h5_validate_write) export(h5_write_values) export(list_h5_datasets) +export(make_swale_runner) +export(optimise_swale_design) export(plot_cost_overflow_boxplot) export(plot_cost_vs_evaporation) export(plot_cost_vs_overflow_volume) @@ -31,8 +36,11 @@ export(read_hdf5_connections) export(read_hdf5_scalars) export(read_hdf5_timeseries) export(read_raindrop_errors) +export(read_site_timeseries) export(run_model) export(run_scenarios) +export(sickerbox_level_presets) +export(stack_levels) importFrom(dplyr,"%>%") importFrom(dplyr,across) importFrom(dplyr,all_of) diff --git a/NEWS.md b/NEWS.md index 3f6bb77..fd71f28 100644 --- a/NEWS.md +++ b/NEWS.md @@ -2,6 +2,70 @@ ## New features +* New **swale-design optimiser** — finds the cost-minimal design per + overflow target (`n_overflows <= x`) with roughly ten engine runs per + (storage type, target) cell instead of a full factorial sweep, at finer + resolution (2 m² / 10 mm instead of 25 m² / 100 mm grid steps): + - `find_min_feasible()` — the single primitive everything is built + from: bisection for the smallest feasible value of one design + parameter, over continuous bounds (`lower`/`upper`/`tol`) or discrete + stack `levels` (binary search). Evaluations are memoised and two + safety rules from the monotonicity analysis are built in: an **edge + guard** for the +1 event-counting wobble of the 4-h event separation + (a descending ladder below an "infeasible by exactly one event" upper + edge, so the Bad Aussee 175-m²-feasible/200-m²-infeasible pattern + cannot eat a solution) and a **volume referee** that warns — and + flags the result — iff the overflow count *and* the overflow volume + increase together (real non-monotonicity; never observed in the + 5 112 validation comparisons). + - `optimise_swale_design()` — coordinate descent in cost order: shrink + `mulde_area` (the expensive lever) first, then `mulde_height` (the + cheap one); the storage layer starts at its smallest level and is + escalated only when the area is stuck at its upper bound. One shared + evaluation cache spans all `x_targets` and both storage types (a run + classifies itself for every target at once), warm-start brackets are + derived from prior brute-force results (CSV schema of the + workflows), `max_total_depth` adds an analytic depth constraint + (DWA-A 138 groundwater clearance / cover requirements) that costs no + simulation runs, and "infeasible within bounds" is a regular result + status, not an error. Costs are attached via `compute_costs()`; all + evaluated designs ship as attribute `"evaluations"`. + - `make_swale_runner()` — package-level refactoring of the `run_one()` + function previously duplicated across the three case-study + vignettes: one closure factory covering both variants (Eisenstadt: + `base.h5` rain curve scaled by `rain_factor`; Wien / Bad Aussee: own + rain + ET0 series in mm/h incl. the Growth/Shading end-time fix). + Returns the thinned one-row optimisation result augmented with + `overflow_volume_m3` (= `sum_overflows` [mm] × `mulde_area` / 1000). + - `stack_levels()`, `sickerbox_level_presets()`, + `default_storage_spec()`, `default_storage_types()` — storage-layer + search spaces: achievable stack heights from module heights (incl. + mixed combinations such as Rigofill full + half block), manufacturer + presets (GRAF, Fränkische, ACO, Wavin; verify against data sheets + before productive runs) with the brute-force grid levels + 300/600/900/1200 mm as the default, and the gravel-trench range + coupled at 3 × the box range (usable-porosity ratio 0.95 / 0.3). + +* New conditional vignette `monotonicity_analysis` — validates the + optimiser's core assumption on the three brute-force result sets + (5 112 neighbour comparisons): `n_overflows` is quasi-monotone in + every design parameter (13 violations, all +1 counting artefacts of + the 4-h event separation), the overflow volume (in m³) is monotone + without exception, ET depends on `mulde_area` only, and the filter + conductivity is a cost-free dominant lever (fix at maximum). Renders + after the three workflow vignettes into + `vignettes/monotonicity_analysis/` (deploy unit with the plain-language + report `index.html` and the exported `mono_*` detail tables as CSV + + interactive HTML). + +* **testthat suite added** (edition 3; `tests/testthat/`): unit tests for + the bisection primitive (threshold accuracy, run counts, wobble guard, + volume referee, discrete levels) and end-to-end optimiser tests + against a synthetic monotone hydraulic model, verified against a fine + brute-force reference (cost within 5 %, monotone cost-effectiveness + curve, storage escalation, infeasibility handling, warm-start + savings, `max_total_depth`). + * New exported plot `plot_cost_vs_evaporation()` — third cost view: scatters `cost_total` (EUR, x) against the element evapotranspiration share (`element.WB_Evapotranspiration_`, %, y). Points share the diff --git a/R/cost_tooltip.R b/R/cost_tooltip.R index cb492ec..00de3da 100644 --- a/R/cost_tooltip.R +++ b/R/cost_tooltip.R @@ -250,7 +250,7 @@ storage_type_shapes <- function(storage_type, lang = c("de", "en")) { ) } -#' Usable storage volume of the storage layer [m3] per row +#' Usable storage volume of the storage layer (m3) per row #' #' Area x height x usable porosity (thetaS - thetaFC) of the storage type. #' Taken from a precomputed `storage_volume_m3` column when available, @@ -282,7 +282,7 @@ storage_volume_from_df <- function(df) { #' @param df Data frame with the columns listed above. #' @param tt Label list from `cost_tooltip_labels()`. #' @param digits Integer. Rounding for the numeric tooltip values. -#' @param evap_min Numeric. Minimum element evapotranspiration share [%] of +#' @param evap_min Numeric. Minimum element evapotranspiration share (%) of #' the complete model run — the reference for the cost-per-percent line. #' `NULL` falls back to the minimum within `df` (identical as long as `df` #' is unfiltered). diff --git a/R/find_min_feasible.R b/R/find_min_feasible.R new file mode 100644 index 0000000..ef6c228 --- /dev/null +++ b/R/find_min_feasible.R @@ -0,0 +1,196 @@ +#' Smallest feasible parameter value via bisection (monotone threshold search) +#' +#' Core building block of the swale-design optimiser: finds the smallest +#' value of one design parameter for which the overflow target is met +#' (`n_overflows <= x_max`), assuming quasi-monotone feasibility (larger +#' value = never more overflows; verified for the RAINDROP model in the +#' `monotonicity_analysis` vignette). Each evaluation halves the search +#' interval, so `ceiling(log2(range / tol))` evaluations suffice. +#' +#' Two safety rules from the monotonicity analysis are built in: +#' \itemize{ +#' \item \strong{Edge guard}: if the upper bound is infeasible by no more +#' than `wobble` events (the +1 counting artefact of the 4-h event +#' separation), a descending ladder below the edge searches for a +#' feasible anchor before the branch is declared infeasible. +#' \item \strong{Volume referee}: whenever `n_overflows` increases with +#' the parameter (a counting flip), the overflow volume must have +#' decreased; if the volume increased as well, a warning is emitted and +#' `monotonicity_violation` is set (real non-monotonicity -- never +#' observed at the three validation sites). +#' } +#' +#' @param evaluate `function(value)` returning a list / one-row data.frame +#' with at least `n_overflows`; if it also contains `volume_column`, the +#' volume referee is active. Evaluations are memoised per value. +#' @param x_max Feasibility target: feasible iff `n_overflows <= x_max`. +#' @param lower,upper Numeric search bounds (continuous mode). +#' @param tol Resolution of the continuous search (same unit as the value). +#' @param levels Sorted numeric vector of discrete candidate values +#' (discrete mode, e.g. Sickerbox stack heights). If given, `lower`, +#' `upper` and `tol` are ignored and a binary search over the levels is +#' performed. +#' @param wobble Maximum counting-artefact size tolerated by the edge guard +#' (default 1, matching the observed +1 flips). +#' @param volume_column Name of the volume element in the `evaluate` result +#' used by the volume referee (default `"overflow_volume_m3"`). +#' @param verbose Print one line per evaluation. +#' +#' @return List with +#' \describe{ +#' \item{value}{smallest feasible value, or `NA` if infeasible} +#' \item{n_overflows}{overflow count at `value`} +#' \item{status}{`"ok"`, `"at_lower_bound"` (already feasible at the lower +#' end -- caller may widen the bracket) or `"infeasible"`} +#' \item{evaluations}{tibble of all evaluated values (value, n_overflows, +#' volume), sorted by value} +#' \item{n_evaluations}{number of distinct evaluations} +#' \item{monotonicity_violation}{`TRUE` if the volume referee fired} +#' } +#' +#' @examples +#' # synthetic monotone step function: feasible from 137.4 m2 on +#' f <- function(v) list(n_overflows = if (v >= 137.4) 0L else 10L) +#' find_min_feasible(f, x_max = 0, lower = 25, upper = 200, tol = 2)$value +#' +#' @export +find_min_feasible <- function(evaluate, + x_max, + lower = NULL, + upper = NULL, + tol = 1, + levels = NULL, + wobble = 1L, + volume_column = "overflow_volume_m3", + verbose = FALSE) { + + discrete <- !is.null(levels) + if (discrete) { + grid <- sort(unique(levels)) + axis_lo <- 0 # virtual: below the smallest level + axis_hi <- length(grid) + axis_tol <- 1 + to_value <- function(i) grid[[i]] + } else { + stopifnot(is.numeric(lower), is.numeric(upper), upper > lower, tol > 0) + axis_lo <- lower + axis_hi <- upper + axis_tol <- tol + to_value <- identity + } + + evals <- new.env(parent = emptyenv()) + eval_at <- function(axis_pos) { + v <- to_value(axis_pos) + key <- format(v, digits = 15) + if (!is.null(evals[[key]])) return(evals[[key]]) + res <- as.list(evaluate(v)) + if (!"n_overflows" %in% names(res)) { + stop("find_min_feasible(): evaluate() must return an element 'n_overflows'") + } + row <- list( + value = v, + n_overflows = as.numeric(res$n_overflows), + volume = if (volume_column %in% names(res)) { + as.numeric(res[[volume_column]]) + } else { + NA_real_ + } + ) + assign(key, row, envir = evals) + if (isTRUE(verbose)) { + message(sprintf(" eval %s -> n_overflows = %s", + format(v), format(row$n_overflows))) + } + row + } + feasible <- function(axis_pos) { + n <- eval_at(axis_pos)$n_overflows + !is.na(n) && n <= x_max + } + + status <- "ok" + + # 1) Upper edge: with monotone feasibility the whole range is infeasible + # if the upper edge is -- unless the edge is a +1 counting wobble, in + # which case a descending ladder looks for a feasible anchor below. + if (!feasible(axis_hi)) { + n_hi <- eval_at(axis_hi)$n_overflows + anchor <- NA_real_ + if (!is.na(n_hi) && n_hi <= x_max + wobble) { + offset <- axis_tol + repeat { + p <- axis_hi - offset + if (discrete) p <- ceiling(p) + if (p <= axis_lo) break + if (feasible(p)) { + anchor <- p + break + } + offset <- offset * 2 + } + } + if (is.na(anchor)) { + status <- "infeasible" + } else { + axis_hi <- anchor + } + } + + # 2) Bisection: invariant lo infeasible (or virtual/edge), hi feasible. + best <- NA_real_ + if (!identical(status, "infeasible")) { + lo <- axis_lo + hi <- axis_hi + if (!discrete && feasible(lo)) { + hi <- lo # optimum at (or below) the lower end + } + while ((hi - lo) > axis_tol) { + mid <- (lo + hi) / 2 + if (discrete) mid <- floor(mid) + if (mid <= lo || mid >= hi) break + if (feasible(mid)) hi <- mid else lo <- mid + } + best <- hi + at_lower <- if (discrete) best <= 1 else best <= axis_lo + if (at_lower) status <- "at_lower_bound" + } + + # 3) Volume referee over all evaluations of this search + ev <- do.call(rbind, lapply(ls(evals), function(k) { + e <- get(k, envir = evals) + data.frame(value = e$value, n_overflows = e$n_overflows, + volume = e$volume) + })) + ev <- ev[order(ev$value), , drop = FALSE] + rownames(ev) <- NULL + monotonicity_violation <- FALSE + if (nrow(ev) >= 2) { + dn <- diff(ev$n_overflows) + dv <- diff(ev$volume) + flips <- which(!is.na(dn) & dn > 0) + real <- flips[!is.na(dv[flips]) & dv[flips] > 1e-9] + if (length(real) > 0) { + monotonicity_violation <- TRUE + warning(sprintf( + paste0("find_min_feasible(): n_overflows AND overflow volume ", + "increase between value %s and %s -- real non-monotonicity, ", + "bisection result unreliable for this branch."), + format(ev$value[real[1]]), format(ev$value[real[1] + 1]) + ), call. = FALSE) + } + } + + list( + value = if (identical(status, "infeasible")) NA_real_ else to_value(best), + n_overflows = if (identical(status, "infeasible")) { + NA_real_ + } else { + eval_at(best)$n_overflows + }, + status = status, + evaluations = tibble::as_tibble(ev), + n_evaluations = nrow(ev), + monotonicity_violation = monotonicity_violation + ) +} diff --git a/R/make_swale_runner.R b/R/make_swale_runner.R new file mode 100644 index 0000000..ab766aa --- /dev/null +++ b/R/make_swale_runner.R @@ -0,0 +1,188 @@ +#' Capillary suction from hydraulic conductivity (Rawls fit) +#' +#' Suction head Psi_s in mm as a function of the saturated hydraulic +#' conductivity in mm/h, as used by the workflow vignettes. +#' +#' @param kf_mmh Saturated hydraulic conductivity in mm/h. +#' @return Suction head in mm. +#' @keywords internal +#' @noRd +psi_s_mm <- function(kf_mmh) { + (3.237 * (kf_mmh / 25.4)^(-0.328)) * 25.4 +} + +#' Create a site-specific single-scenario runner for the optimiser +#' +#' Factors the `run_one()` function that was duplicated across the three +#' workflow vignettes (Eisenstadt 2005, Wien, Bad Aussee) into one +#' package-level closure factory. The returned function runs the RAINDROP +#' engine for one parameter set and returns the thinned one-row +#' optimisation result (overflow events + water balance), augmented with +#' the input parameters and the overflow volume in m3. +#' +#' Site differences are covered by the arguments: Eisenstadt scales the +#' rain curve shipped in `base.h5` by `rain_factor` (leave +#' `timeseries_rain` = `NULL`), Wien and Bad Aussee replace the rain and +#' ET0 curves entirely (`timeseries_rain` / `timeseries_et`, values in +#' mm/h as written by the vignettes). +#' +#' @param path_list Path definition list as used by the workflow vignettes +#' (resolvable with `kwb.utils::resolve()`, must contain `path_base`, +#' `path_exe`, `dir_input`, `dir_output`, `dir_target_output`, +#' `path_target_input`, `path_results_hdf5_element`, +#' `path_results_hdf5_flaeche`, `file_target`). +#' @param timestep_hours Engine time step in hours (default 0.1). +#' @param timeseries_rain Optional data.frame `time`/`value` (mm/h) written +#' to `//Kurven/Regen`; when given, `//Kurven/Growth_1` and +#' `//Kurven/Shading_1` end times are extended to the rain series end and +#' `rain_factor` is ignored. +#' @param timeseries_et Optional data.frame `time`/`value` (mm/h) written +#' to `//Kurven/ET0`. +#' @param storage_types Soil presets of the storage layer per storage type, +#' see [default_storage_types()]. +#' @param event_separation_hours Event separation for overflow counting +#' (default 4, as in the vignettes and the monotonicity analysis). +#' @param scenario_prefix Prefix for generated scenario names (default +#' `"o"` -> `o00001`, `o00002`, ... -- distinct from the grid runs +#' `s00001` ...). +#' @param debug Passed on to the engine/reader helpers. +#' +#' @return `function(params)` where `params` is a named list (or one-row +#' data.frame) with `mulde_area`, `mulde_height` (mm), `storage_type`, +#' `storage_height` (mm), `connected_area` (m2), `filter_height` (mm), +#' `filter_hydraulicconductivity` (mm/h), `bottom_hydraulicconductivity` +#' (mm/h) and optionally `rain_factor` (default 1) and `lai` +#' (default 3.9). It returns a one-row tibble with the parameters, the +#' scenario name and the optimisation metrics (`n_overflows`, +#' `sum_overflows` in mm, `overflow_volume_m3`, water-balance shares). +#' +#' @seealso [optimise_swale_design()], [find_min_feasible()] +#' @export +make_swale_runner <- function(path_list, + timestep_hours = 0.1, + timeseries_rain = NULL, + timeseries_et = NULL, + storage_types = default_storage_types(), + event_separation_hours = 4, + scenario_prefix = "o", + debug = FALSE) { + + counter <- 0L + + function(params) { + params <- as.list(params) + required <- c("mulde_area", "mulde_height", "storage_type", + "storage_height", "connected_area", "filter_height", + "filter_hydraulicconductivity", + "bottom_hydraulicconductivity") + missing <- setdiff(required, names(params)) + if (length(missing) > 0) { + stop("make_swale_runner(): params is missing: ", + paste(missing, collapse = ", ")) + } + st <- storage_types[[params$storage_type]] + if (is.null(st)) { + stop("make_swale_runner(): unknown storage_type '", + params$storage_type, "'") + } + rain_factor <- if (is.null(params$rain_factor)) 1 else params$rain_factor + lai <- if (is.null(params$lai)) 3.9 else params$lai + + counter <<- counter + 1L + s_name <- sprintf("%s%05d", scenario_prefix, counter) + paths <- kwb.utils::resolve(path_list, dir_target = s_name) + + fs::dir_create(paths$dir_input, recurse = TRUE) + fs::dir_create(paths$dir_output, recurse = TRUE) + fs::dir_create(paths$dir_target_output, recurse = TRUE) + + fs::file_copy(path = paths$path_base, + new_path = paths$path_target_input, + overwrite = TRUE) + + h5 <- hdf5r::H5File$new(paths$path_target_input, mode = "a") + on.exit(try(h5$close_all(), silent = TRUE), add = TRUE) + + new_path <- stringr::str_c( + normalizePath(fs::path_abs(paths$dir_target_output)), "\\" + ) + + vals <- h5_read_values(h5) + + vals$`//Berechnungsparameter/Ergebnispfad` <- new_path + vals$`//Berechnungsparameter/Zeitschritt_Infiltration` <- timestep_hours + vals$`//Berechnungsparameter/Zeitschritt_ET` <- timestep_hours + vals$`//Berechnungsparameter/Zeitschritt_Verschaltungen` <- timestep_hours + vals$`//Berechnungsparameter/R-Plots` <- 0 + vals$`//Berechnungsparameter/Ausgabemodus` <- "Optimierung" + vals$`//Berechnungsparameter/Evapotranspiration_aktiv` <- 1 + + vals$`//Massnahmenelemente/Dach/Berechnungsparameter/Evapotranspiration_aktiv` <- 1 + vals$`//Massnahmenelemente/Dach/Allgemein/Flaeche` <- params$connected_area + + vals$`//Massnahmenelemente/Mulde_Rigole/Berechnungsparameter/Evapotranspiration_aktiv` <- 1 + vals$`//Massnahmenelemente/Mulde_Rigole/Allgemein/Regen-Skalierungsfaktor` <- 1 + vals$`//Massnahmenelemente/Mulde_Rigole/Allgemein/Flaeche` <- params$mulde_area + vals$`//Massnahmenelemente/Mulde_Rigole/Eigenschaften_Oberflaeche/Ueberlaufhoehe` <- params$mulde_height + vals$`//Massnahmenelemente/Mulde_Rigole/Bodenschichtung/Startwerte_theta_ActualSoilMoisture` <- + c(0.3, st$Startwerte_theta_ActualSoilMoisture) + vals$`//Massnahmenelemente/Mulde_Rigole/Bodenschichtung/Schichtdicken` <- + c(params$filter_height, params$storage_height) + vals$`//Bodenarten/Speicher/thetaWP_MoistureAtWiltingPoint` <- st$thetaWP_MoistureAtWiltingPoint + vals$`//Bodenarten/Speicher/thetaFC_MoistureAtFieldCapacity` <- st$thetaFC_MoistureAtFieldCapacity + vals$`//Bodenarten/Speicher/thetaS_MoistureAtSaturation` <- st$thetaS_MoistureAtSaturation + vals$`//Massnahmenelemente/Mulde_Rigole/Allgemein/Endversickerungsrate` <- + params$bottom_hydraulicconductivity + vals$`//Massnahmenelemente/Mulde_Rigole/Parameter_Evapotranspiration/LAI_LeafAreaIndex` <- lai + + vals$`//Bodenarten/Bodenfilter/Ks_HydraulicConductivity` <- + params$filter_hydraulicconductivity + vals$`//Bodenarten/Bodenfilter/Psi_Saugspannung_CapillarySuction` <- + psi_s_mm(params$filter_hydraulicconductivity) + + if (!is.null(timeseries_et)) { + vals$`//Kurven/ET0` <- timeseries_et + } + if (!is.null(timeseries_rain)) { + vals$`//Kurven/Regen` <- timeseries_rain + vals$`//Kurven/Growth_1`$time[2] <- max(timeseries_rain$time) + vals$`//Kurven/Shading_1`$time[2] <- max(timeseries_rain$time) + } else if (is.data.frame(vals[["//Kurven/Regen"]])) { + vals[["//Kurven/Regen"]]$value <- + vals[["//Kurven/Regen"]]$value * rain_factor + } + + h5_write_values(h5, vals, resize = TRUE, + scalar_strategy = "error", verbose = FALSE) + h5$close_all() + + run_model(path_exe = paths$path_exe, + path_input = paths$path_target_input, + debug = debug) + + # Thin immediately (lean read), exactly like the vignettes' run_one() + sim_one <- get_simulation_results_optim( + paths = paths, + path_list = path_list, + simulation_names = s_name, + debug = debug, + lean = TRUE + ) + + row <- add_overflow_events_and_waterbalance( + simulation_results = sim_one, + event_separation_hours = event_separation_hours, + canonical_variables = default_canonical_wb_variables() + ) + + dplyr::bind_cols( + tibble::as_tibble(params[required]), + tibble::tibble(rain_factor = rain_factor, lai = lai), + row + ) %>% + dplyr::mutate( + # sum_overflows is mm water column over the swale area + overflow_volume_m3 = .data$sum_overflows * .data$mulde_area / 1000 + ) + } +} diff --git a/R/optimise_swale_design.R b/R/optimise_swale_design.R new file mode 100644 index 0000000..a5ae954 --- /dev/null +++ b/R/optimise_swale_design.R @@ -0,0 +1,282 @@ +#' Warm-start area bracket from prior (brute-force) results +#' +#' Narrows the mulde_area search interval to one grid step around the +#' cheapest feasible grid cell of the matching branch, if prior results +#' contain it. Falls back to the full bounds otherwise. +#' +#' @keywords internal +#' @noRd +area_bracket_from_prior <- function(prior, type, h_s, h_m, x, bounds) { + needed <- c("storage_type", "storage_height", "mulde_height", "mulde_area", + "n_overflows", "filter_hydraulicconductivity") + if (is.null(prior) || !all(needed %in% names(prior))) return(bounds) + kf_max <- suppressWarnings( + max(prior$filter_hydraulicconductivity, na.rm = TRUE) + ) + d <- prior[prior$storage_type == type & + prior$filter_hydraulicconductivity == kf_max & + prior$storage_height == h_s & + prior$mulde_height == h_m, , drop = FALSE] + if (nrow(d) < 2) return(bounds) + ok <- d$mulde_area[!is.na(d$n_overflows) & d$n_overflows <= x] + if (length(ok) == 0) return(bounds) + first_ok <- min(ok) + areas <- sort(unique(d$mulde_area)) + step <- if (length(areas) > 1) min(diff(areas)) else diff(bounds) + lo <- max(bounds[1], first_ok - step) + hi <- min(bounds[2], first_ok) + if (hi <= lo) bounds else c(lo, hi) +} + +#' Find the cost-optimal swale design per overflow target +#' +#' Coordinate-descent optimiser built from a single primitive +#' ([find_min_feasible()], bisection over one parameter): shrink the +#' expensive lever first (`mulde_area`), then the cheap one +#' (`mulde_height`); the storage layer starts at its smallest level and is +#' only escalated when the area is stuck at its upper bound. The filter +#' conductivity is expected to be fixed at the maximum via `fixed` (it is +#' cost-free and dominant, see the `monotonicity_analysis` vignette). Every +#' engine run is cached, so the sweep over all `x_targets` and both storage +#' types shares evaluations. +#' +#' @param run_fn `function(params)` running one scenario and returning at +#' least `n_overflows` plus `sum_overflows` (mm) or `overflow_volume_m3`; +#' typically created with [make_swale_runner()]. `params` is a named list +#' of `mulde_area`, `mulde_height`, `storage_type`, `storage_height` plus +#' everything in `fixed`. +#' @param x_targets Integer vector of overflow targets (feasible :<=> +#' `n_overflows <= x`), default `0:5`. +#' @param area_bounds,area_tol Search range (m2) and resolution for +#' `mulde_area`. +#' @param height_bounds,height_tol Search range (mm) and resolution for +#' `mulde_height`. +#' @param storage_spec Storage search space per type, see +#' [default_storage_spec()]: discrete `levels` (infiltration box) or +#' continuous `bounds` + `tol` (gravel trench). +#' @param fixed Named list of parameters passed unchanged to `run_fn` +#' (connected area, filter geometry, kf at maximum, ...). Must contain +#' `filter_height` for the cost model. +#' @param prior_results Optional data.frame with prior (grid) results in +#' the workflow CSV schema, used as warm start (narrows the first area +#' bracket to one grid step). +#' @param max_total_depth Optional analytic depth constraint in mm: +#' `mulde_height + filter_height + storage_height <= max_total_depth` +#' (e.g. from DWA-A 138 groundwater clearance or cover requirements). +#' Enforced without any simulation runs. +#' @param cost_rates Unit costs, see [default_cost_rates()]. +#' @param verbose Print one progress line per solved cell. +#' +#' @return Tibble with one row per (storage type, x): the optimal design +#' (`mulde_area`, `mulde_height`, `storage_height`), its metrics +#' (`n_overflows`, `overflow_volume_m3`, `et_pct`), cost columns from +#' [compute_costs()], a `status` (`"ok"` or `"infeasible_within_bounds"`), +#' `monotonicity_warning` (volume referee) and `n_runs_new` (fresh engine +#' runs spent on this cell). All evaluated designs are attached as +#' attribute `"evaluations"`. +#' +#' @seealso [find_min_feasible()], [make_swale_runner()], +#' [default_storage_spec()] +#' @export +optimise_swale_design <- function(run_fn, + x_targets = 0:5, + area_bounds = c(25, 200), + area_tol = 2, + height_bounds = c(100, 300), + height_tol = 10, + storage_spec = default_storage_spec(), + fixed = list( + connected_area = 1000, + filter_height = 300, + filter_hydraulicconductivity = 360, + bottom_hydraulicconductivity = 12 + ), + prior_results = NULL, + max_total_depth = NULL, + cost_rates = default_cost_rates(), + verbose = TRUE) { + + stopifnot(is.function(run_fn), !is.null(fixed$filter_height)) + filter_height <- fixed$filter_height + + # --- shared evaluation cache (one engine run per distinct design) ------- + cache <- new.env(parent = emptyenv()) + runs_executed <- 0L + + eval_design <- function(type, area, h_m, h_s) { + key <- paste(type, format(area, digits = 10), format(h_m, digits = 10), + format(h_s, digits = 10), sep = "|") + hit <- cache[[key]] + if (!is.null(hit)) return(hit) + params <- c(list(mulde_area = area, mulde_height = h_m, + storage_type = type, storage_height = h_s), fixed) + res <- as.list(run_fn(params)) + if (!"n_overflows" %in% names(res)) { + stop("optimise_swale_design(): run_fn() must return 'n_overflows'") + } + vol <- res[["overflow_volume_m3"]] + if (is.null(vol) && !is.null(res[["sum_overflows"]])) { + vol <- res[["sum_overflows"]] * area / 1000 # mm x m2 / 1000 = m3 + } + et <- res[["element.WB_Evapotranspiration_"]] + out <- list(storage_type = type, mulde_area = area, mulde_height = h_m, + storage_height = h_s, + n_overflows = as.numeric(res$n_overflows), + overflow_volume_m3 = if (is.null(vol)) NA_real_ else as.numeric(vol), + et_pct = if (is.null(et)) NA_real_ else as.numeric(et)) + runs_executed <<- runs_executed + 1L + assign(key, out, envir = cache) + out + } + + # --- analytic depth constraint ------------------------------------------ + hm_upper <- function(h_s) { + up <- rep(height_bounds[2], length(h_s)) + if (!is.null(max_total_depth)) { + up <- pmin(up, max_total_depth - filter_height - h_s) + } + up + } + + # area search with warm-start bracket; widens the bracket when the + # optimum turns out to lie below it + search_area <- function(eval_a, x, bracket) { + res <- find_min_feasible(eval_a, x_max = x, + lower = bracket[1], upper = bracket[2], + tol = area_tol) + if (identical(res$status, "at_lower_bound") && + bracket[1] > area_bounds[1]) { + res <- find_min_feasible(eval_a, x_max = x, + lower = area_bounds[1], upper = bracket[1], + tol = area_tol) + } + if (identical(res$status, "infeasible") && + bracket[2] < area_bounds[2]) { + # warm start was too optimistic -> retry up to the full upper bound + res <- find_min_feasible(eval_a, x_max = x, + lower = bracket[1], upper = area_bounds[2], + tol = area_tol) + } + res + } + + # --- solve one (storage type, x) cell ------------------------------------ + solve_cell <- function(type, x) { + runs_before <- runs_executed + spec <- storage_spec[[type]] + if (is.null(spec)) { + stop("optimise_swale_design(): storage_spec has no entry '", type, "'") + } + discrete <- !is.null(spec$levels) + mono_warn <- FALSE + + infeasible_row <- function() tibble::tibble( + x = x, storage_type = type, status = "infeasible_within_bounds", + mulde_area = NA_real_, mulde_height = NA_real_, + storage_height = NA_real_, n_overflows = NA_real_, + overflow_volume_m3 = NA_real_, et_pct = NA_real_, + monotonicity_warning = mono_warn, + n_runs_new = runs_executed - runs_before + ) + + if (discrete) { + levels_all <- sort(spec$levels) + levels_all <- levels_all[hm_upper(levels_all) >= height_bounds[1]] + if (length(levels_all) == 0) return(infeasible_row()) + h_s <- levels_all[1] + } else { + gb <- spec$bounds + if (!is.null(max_total_depth)) { + gb[2] <- min(gb[2], max_total_depth - filter_height - height_bounds[1]) + } + if (gb[2] <= gb[1]) return(infeasible_row()) + gravel_tol <- if (is.null(spec$tol)) 25 else spec$tol + h_s <- gb[1] + } + + a_star <- NA_real_ + repeat { + h_m_up <- hm_upper(h_s) + eval_a <- function(a) eval_design(type, a, h_m_up, h_s) + bracket <- area_bracket_from_prior(prior_results, type, h_s, h_m_up, + x, area_bounds) + res_a <- search_area(eval_a, x, bracket) + mono_warn <- mono_warn || res_a$monotonicity_violation + if (!identical(res_a$status, "infeasible")) { + a_star <- res_a$value + break + } + # area stuck at the upper bound -> escalate the storage layer + eval_s <- function(h) eval_design(type, area_bounds[2], hm_upper(h), h) + if (discrete) { + rest <- levels_all[levels_all > h_s] + if (length(rest) == 0) return(infeasible_row()) + res_s <- find_min_feasible(eval_s, x_max = x, levels = rest) + } else { + if (h_s >= gb[2]) return(infeasible_row()) + res_s <- find_min_feasible(eval_s, x_max = x, + lower = h_s, upper = gb[2], + tol = gravel_tol) + } + mono_warn <- mono_warn || res_s$monotonicity_violation + if (identical(res_s$status, "infeasible")) return(infeasible_row()) + h_s <- res_s$value + } + + # shrink the cheap lever last: mulde_height at fixed (a*, h_s). + # A second area pass is provably redundant: a smaller mulde_height + # only weakens the hydraulics, so the minimal feasible area cannot + # decrease any further. + h_m_up <- hm_upper(h_s) + h_m_star <- h_m_up + if (h_m_up > height_bounds[1]) { + res_h <- find_min_feasible( + function(h) eval_design(type, a_star, h, h_s), + x_max = x, lower = height_bounds[1], upper = h_m_up, + tol = height_tol + ) + mono_warn <- mono_warn || res_h$monotonicity_violation + if (!identical(res_h$status, "infeasible")) h_m_star <- res_h$value + } + + final <- eval_design(type, a_star, h_m_star, h_s) + if (isTRUE(verbose)) { + message(sprintf( + "[%s | x = %d] area %s m2, height %s mm, storage %s mm (%d neue Laeufe)", + type, x, format(a_star), format(h_m_star), format(h_s), + runs_executed - runs_before + )) + } + tibble::tibble( + x = x, storage_type = type, status = "ok", + mulde_area = a_star, mulde_height = h_m_star, storage_height = h_s, + n_overflows = final$n_overflows, + overflow_volume_m3 = final$overflow_volume_m3, + et_pct = final$et_pct, + monotonicity_warning = mono_warn, + n_runs_new = runs_executed - runs_before + ) + } + + # --- sweep all cells (shared cache makes repeats cheap) ------------------ + cells <- expand.grid(type = names(storage_spec), + x = sort(unique(as.integer(x_targets))), + stringsAsFactors = FALSE) + out <- dplyr::bind_rows( + lapply(seq_len(nrow(cells)), + function(i) solve_cell(cells$type[i], cells$x[i])) + ) + + out$filter_height <- filter_height + out <- compute_costs(out, cost_rates = cost_rates) + out <- dplyr::arrange(out, .data$storage_type, .data$x) + + evaluations <- dplyr::bind_rows( + lapply(ls(cache), function(k) tibble::as_tibble(get(k, envir = cache))) + ) + attr(out, "evaluations") <- dplyr::arrange( + evaluations, .data$storage_type, .data$mulde_area + ) + attr(out, "n_runs_total") <- runs_executed + out +} diff --git a/R/read_site_timeseries.R b/R/read_site_timeseries.R new file mode 100644 index 0000000..266abf2 --- /dev/null +++ b/R/read_site_timeseries.R @@ -0,0 +1,83 @@ +#' Read and prepare site rain/ET0 time series for the engine (mm/h) +#' +#' Factors the time-series preparation duplicated in the Wien and Bad +#' Aussee workflow vignettes into one helper: reads the shipped GeoSphere +#' rain series (`rain.csv.gz`: columns `time` (datetime), `rr` (mm per +#' interval), `station`, further columns tolerated) and reference ET0 +#' series (`et.csv`: +#' `date;value` with `dd.mm.yyyy`, mm per day), converts both to hours +#' since series start, aligns the series ends (the shorter series is +#' extended to the longer one's end, repeating its last value) and +#' converts the values to the engine's **mm/h** rate convention (rain: +#' mm per interval / interval hours; ET0: mm per day / 24). +#' +#' @param path_rain Path to the rain CSV (may be gzipped). +#' @param path_et Path to the ET0 CSV (semicolon separated). +#' @param verbose Print alignment messages (default TRUE). +#' +#' @return List with data.frames `rain` and `et` (columns `time` = hours +#' since start, `value` = mm/h) ready for +#' `make_swale_runner(timeseries_rain = , timeseries_et = )`. +#' +#' @seealso [make_swale_runner()] +#' @export +read_site_timeseries <- function(path_rain, path_et, verbose = TRUE) { + + if (!requireNamespace("readr", quietly = TRUE)) { + stop("read_site_timeseries() requires the 'readr' package") + } + + timeseries_et <- readr::read_delim(path_et, delim = ";", + col_types = "cd") %>% + dplyr::mutate( + date = lubridate::dmy(.data$date), + time = as.integer(difftime(.data$date, min(.data$date), + units = "hours")) + ) %>% + dplyr::select(-"date") %>% + dplyr::filter(!is.na(.data$value)) %>% + dplyr::relocate("time", .before = "value") + timeseries_et$time[nrow(timeseries_et)] <- + ceiling(timeseries_et$time[nrow(timeseries_et)]) + + timeseries_rain <- readr::read_csv(path_rain, show_col_types = FALSE) %>% + dplyr::rename(datetime = "time", value = "rr") %>% + dplyr::mutate( + time = as.double(difftime(.data$datetime, min(.data$datetime), + units = "secs")) / 3600 + ) %>% + dplyr::filter(!is.na(.data$value)) %>% + # robust gegen zusaetzliche Spalten (Bad Aussee hat z.B. "substation"): + # die Engine erwartet exakt time + value + dplyr::select("time", "value") + timeseries_rain$time[nrow(timeseries_rain)] <- + ceiling(timeseries_rain$time[nrow(timeseries_rain)]) + + # Serien-Enden angleichen: die kuerzere Serie wird mit ihrem letzten + # Wert bis zum Ende der laengeren verlaengert (wie in den Vignetten). + extend_to <- function(df, t_end, label) { + if (t_end <= max(df$time)) return(df) + if (isTRUE(verbose)) { + message(sprintf( + "%s series extended by %.1f hours to %.0f h (last value %.4f)", + label, t_end - max(df$time), t_end, df$value[nrow(df)] + )) + } + dplyr::bind_rows(df, tibble::tibble(time = t_end, + value = df$value[nrow(df)])) + } + t_end <- max(max(timeseries_rain$time), max(timeseries_et$time)) + timeseries_et <- extend_to(timeseries_et, t_end, "ET0") + timeseries_rain <- extend_to(timeseries_rain, t_end, "Rain") + + # mm je Intervall -> mm/h (Engine liest beide Kurven als mm/h-Rate; + # ET0-Tageswerte ohne /24 wuerden 24x zu hoch integriert) + period_rain <- c(diff(timeseries_rain$time), + mean(diff(timeseries_rain$time))) + timeseries_rain$value <- timeseries_rain$value / period_rain + period_et <- c(diff(timeseries_et$time), mean(diff(timeseries_et$time))) + timeseries_et$value <- timeseries_et$value / period_et + + list(rain = as.data.frame(timeseries_rain), + et = as.data.frame(timeseries_et)) +} diff --git a/R/stack_levels.R b/R/stack_levels.R new file mode 100644 index 0000000..b13562a --- /dev/null +++ b/R/stack_levels.R @@ -0,0 +1,115 @@ +#' Achievable storage-layer stack heights from module heights +#' +#' Enumerates all storage-layer heights that can be built by stacking +#' (and mixing) the given module heights, e.g. full blocks combined with +#' at most one half block. +#' +#' @param modules Numeric vector of module heights in mm (e.g. `c(660, 350)` +#' for a full block plus a half block). +#' @param max_count Integer vector (recycled to `length(modules)`): maximum +#' number of modules of each type in one stack. Defaults to 7 for every +#' module (cf. GRAF EcoBloc smart, stackable up to 7 layers). +#' @param max_height Maximum total stack height in mm (default 2600). +#' +#' @return Sorted numeric vector of achievable stack heights in mm. +#' +#' @examples +#' stack_levels(360) # 360, 720, ..., 2520 +#' stack_levels(c(660, 350), max_count = c(7, 1)) # Rigofill full + half block +#' +#' @export +stack_levels <- function(modules, + max_count = rep(7L, length(modules)), + max_height = 2600) { + stopifnot(is.numeric(modules), all(modules > 0)) + max_count <- rep_len(as.integer(max_count), length(modules)) + counts <- Map(function(m, k) 0:min(k, floor(max_height / m)), + modules, max_count) + grid <- do.call(expand.grid, counts) + h <- as.vector(as.matrix(grid) %*% modules) + sort(unique(h[h > 0 & h <= max_height])) +} + +#' Sickerbox storage-height presets (brute force default + manufacturers) +#' +#' Named list of storage-height level vectors (mm) for the infiltration-box +#' storage layer. `brute_force` is the default used by the workflow +#' vignettes (300/600/900/1200 mm -- itself a combination of several box +#' types). The manufacturer presets are generated with [stack_levels()] +#' from typical module heights of commercial block systems; verify against +#' the current data sheets before productive optimisation runs. +#' +#' @param max_height Maximum total stack height in mm passed to +#' [stack_levels()] (default 2600). +#' +#' @return Named list of sorted numeric vectors (mm). +#' +#' @export +sickerbox_level_presets <- function(max_height = 2600) { + list( + brute_force = c(300, 600, 900, 1200), + graf_ecobloc_smart = stack_levels(360, max_height = max_height), + graf_ecobloc_420 = stack_levels(660, max_height = max_height), + fraenkische_rigofill = stack_levels(c(660, 350), max_count = c(7L, 1L), + max_height = max_height), + aco_stormbrixx_hd = stack_levels(614, max_height = max_height), + aco_stormbrixx_sd = stack_levels(342, max_height = max_height), + wavin_aquacell = stack_levels(400, max_height = max_height) + ) +} + +#' Default storage specification for the swale-design optimiser +#' +#' Storage-layer search space per storage type: the infiltration box uses +#' discrete stack levels (default: the brute-force grid levels), the gravel +#' trench is continuous with bounds coupled to the box level range by +#' `coupling_factor` (default 3, approximating the usable-porosity ratio +#' 0.95 / 0.3). +#' +#' @param levels Numeric vector of infiltration-box stack heights in mm. +#' @param coupling_factor Factor between gravel-trench bounds and the box +#' level range. +#' @param gravel_tol Bisection tolerance for the continuous gravel-trench +#' height in mm. +#' +#' @return Named list with entries `infiltration_box` (with `levels`) and +#' `gravel_trench` (with `bounds` and `tol`). +#' +#' @export +default_storage_spec <- function(levels = sickerbox_level_presets()$brute_force, + coupling_factor = 3, + gravel_tol = 25) { + list( + infiltration_box = list(levels = sort(unique(levels))), + gravel_trench = list(bounds = coupling_factor * range(levels), + tol = gravel_tol) + ) +} + +#' Default storage-type soil presets (Speicher layer) +#' +#' Soil parameters of the storage (2nd Bodenschichtung) layer per storage +#' type, as used by the workflow vignettes: infiltration box ("Sickerbox", +#' thetaS 0.95) and gravel trench ("Schotterrigol", thetaS 0.3). +#' +#' @return Named list (per storage type) of lists with +#' `Startwerte_theta_ActualSoilMoisture`, `thetaWP_MoistureAtWiltingPoint`, +#' `thetaFC_MoistureAtFieldCapacity`, `thetaS_MoistureAtSaturation`. +#' +#' @export +default_storage_types <- function() { + list( + infiltration_box = list( + Startwerte_theta_ActualSoilMoisture = 0, + thetaWP_MoistureAtWiltingPoint = 0, + thetaFC_MoistureAtFieldCapacity = 0, + thetaS_MoistureAtSaturation = 0.95 + ), + gravel_trench = list( + Startwerte_theta_ActualSoilMoisture = 0, + thetaWP_MoistureAtWiltingPoint = 0, + thetaFC_MoistureAtFieldCapacity = 0, + thetaS_MoistureAtSaturation = 0.3 + ) + ) +} diff --git a/man/default_storage_spec.Rd b/man/default_storage_spec.Rd new file mode 100644 index 0000000..6e7347a --- /dev/null +++ b/man/default_storage_spec.Rd @@ -0,0 +1,32 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/stack_levels.R +\name{default_storage_spec} +\alias{default_storage_spec} +\title{Default storage specification for the swale-design optimiser} +\usage{ +default_storage_spec( + levels = sickerbox_level_presets()$brute_force, + coupling_factor = 3, + gravel_tol = 25 +) +} +\arguments{ +\item{levels}{Numeric vector of infiltration-box stack heights in mm.} + +\item{coupling_factor}{Factor between gravel-trench bounds and the box +level range.} + +\item{gravel_tol}{Bisection tolerance for the continuous gravel-trench +height in mm.} +} +\value{ +Named list with entries \code{infiltration_box} (with \code{levels}) and +\code{gravel_trench} (with \code{bounds} and \code{tol}). +} +\description{ +Storage-layer search space per storage type: the infiltration box uses +discrete stack levels (default: the brute-force grid levels), the gravel +trench is continuous with bounds coupled to the box level range by +\code{coupling_factor} (default 3, approximating the usable-porosity ratio +0.95 / 0.3). +} diff --git a/man/default_storage_types.Rd b/man/default_storage_types.Rd new file mode 100644 index 0000000..7e951c0 --- /dev/null +++ b/man/default_storage_types.Rd @@ -0,0 +1,18 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/stack_levels.R +\name{default_storage_types} +\alias{default_storage_types} +\title{Default storage-type soil presets (Speicher layer)} +\usage{ +default_storage_types() +} +\value{ +Named list (per storage type) of lists with +\code{Startwerte_theta_ActualSoilMoisture}, \code{thetaWP_MoistureAtWiltingPoint}, +\code{thetaFC_MoistureAtFieldCapacity}, \code{thetaS_MoistureAtSaturation}. +} +\description{ +Soil parameters of the storage (2nd Bodenschichtung) layer per storage +type, as used by the workflow vignettes: infiltration box ("Sickerbox", +thetaS 0.95) and gravel trench ("Schotterrigol", thetaS 0.3). +} diff --git a/man/find_min_feasible.Rd b/man/find_min_feasible.Rd new file mode 100644 index 0000000..b903494 --- /dev/null +++ b/man/find_min_feasible.Rd @@ -0,0 +1,83 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/find_min_feasible.R +\name{find_min_feasible} +\alias{find_min_feasible} +\title{Smallest feasible parameter value via bisection (monotone threshold search)} +\usage{ +find_min_feasible( + evaluate, + x_max, + lower = NULL, + upper = NULL, + tol = 1, + levels = NULL, + wobble = 1L, + volume_column = "overflow_volume_m3", + verbose = FALSE +) +} +\arguments{ +\item{evaluate}{\verb{function(value)} returning a list / one-row data.frame +with at least \code{n_overflows}; if it also contains \code{volume_column}, the +volume referee is active. Evaluations are memoised per value.} + +\item{x_max}{Feasibility target: feasible iff \code{n_overflows <= x_max}.} + +\item{lower, upper}{Numeric search bounds (continuous mode).} + +\item{tol}{Resolution of the continuous search (same unit as the value).} + +\item{levels}{Sorted numeric vector of discrete candidate values +(discrete mode, e.g. Sickerbox stack heights). If given, \code{lower}, +\code{upper} and \code{tol} are ignored and a binary search over the levels is +performed.} + +\item{wobble}{Maximum counting-artefact size tolerated by the edge guard +(default 1, matching the observed +1 flips).} + +\item{volume_column}{Name of the volume element in the \code{evaluate} result +used by the volume referee (default \code{"overflow_volume_m3"}).} + +\item{verbose}{Print one line per evaluation.} +} +\value{ +List with +\describe{ +\item{value}{smallest feasible value, or \code{NA} if infeasible} +\item{n_overflows}{overflow count at \code{value}} +\item{status}{\code{"ok"}, \code{"at_lower_bound"} (already feasible at the lower +end -- caller may widen the bracket) or \code{"infeasible"}} +\item{evaluations}{tibble of all evaluated values (value, n_overflows, +volume), sorted by value} +\item{n_evaluations}{number of distinct evaluations} +\item{monotonicity_violation}{\code{TRUE} if the volume referee fired} +} +} +\description{ +Core building block of the swale-design optimiser: finds the smallest +value of one design parameter for which the overflow target is met +(\code{n_overflows <= x_max}), assuming quasi-monotone feasibility (larger +value = never more overflows; verified for the RAINDROP model in the +\code{monotonicity_analysis} vignette). Each evaluation halves the search +interval, so \code{ceiling(log2(range / tol))} evaluations suffice. +} +\details{ +Two safety rules from the monotonicity analysis are built in: +\itemize{ +\item \strong{Edge guard}: if the upper bound is infeasible by no more +than \code{wobble} events (the +1 counting artefact of the 4-h event +separation), a descending ladder below the edge searches for a +feasible anchor before the branch is declared infeasible. +\item \strong{Volume referee}: whenever \code{n_overflows} increases with +the parameter (a counting flip), the overflow volume must have +decreased; if the volume increased as well, a warning is emitted and +\code{monotonicity_violation} is set (real non-monotonicity -- never +observed at the three validation sites). +} +} +\examples{ +# synthetic monotone step function: feasible from 137.4 m2 on +f <- function(v) list(n_overflows = if (v >= 137.4) 0L else 10L) +find_min_feasible(f, x_max = 0, lower = 25, upper = 200, tol = 2)$value + +} diff --git a/man/make_swale_runner.Rd b/man/make_swale_runner.Rd new file mode 100644 index 0000000..b9e346d --- /dev/null +++ b/man/make_swale_runner.Rd @@ -0,0 +1,74 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/make_swale_runner.R +\name{make_swale_runner} +\alias{make_swale_runner} +\title{Create a site-specific single-scenario runner for the optimiser} +\usage{ +make_swale_runner( + path_list, + timestep_hours = 0.1, + timeseries_rain = NULL, + timeseries_et = NULL, + storage_types = default_storage_types(), + event_separation_hours = 4, + scenario_prefix = "o", + debug = FALSE +) +} +\arguments{ +\item{path_list}{Path definition list as used by the workflow vignettes +(resolvable with \code{kwb.utils::resolve()}, must contain \code{path_base}, +\code{path_exe}, \code{dir_input}, \code{dir_output}, \code{dir_target_output}, +\code{path_target_input}, \code{path_results_hdf5_element}, +\code{path_results_hdf5_flaeche}, \code{file_target}).} + +\item{timestep_hours}{Engine time step in hours (default 0.1).} + +\item{timeseries_rain}{Optional data.frame \code{time}/\code{value} (mm/h) written +to \verb{//Kurven/Regen}; when given, \verb{//Kurven/Growth_1} and +\verb{//Kurven/Shading_1} end times are extended to the rain series end and +\code{rain_factor} is ignored.} + +\item{timeseries_et}{Optional data.frame \code{time}/\code{value} (mm/h) written +to \verb{//Kurven/ET0}.} + +\item{storage_types}{Soil presets of the storage layer per storage type, +see \code{\link[=default_storage_types]{default_storage_types()}}.} + +\item{event_separation_hours}{Event separation for overflow counting +(default 4, as in the vignettes and the monotonicity analysis).} + +\item{scenario_prefix}{Prefix for generated scenario names (default +\code{"o"} -> \code{o00001}, \code{o00002}, ... -- distinct from the grid runs +\code{s00001} ...).} + +\item{debug}{Passed on to the engine/reader helpers.} +} +\value{ +\verb{function(params)} where \code{params} is a named list (or one-row +data.frame) with \code{mulde_area}, \code{mulde_height} (mm), \code{storage_type}, +\code{storage_height} (mm), \code{connected_area} (m2), \code{filter_height} (mm), +\code{filter_hydraulicconductivity} (mm/h), \code{bottom_hydraulicconductivity} +(mm/h) and optionally \code{rain_factor} (default 1) and \code{lai} +(default 3.9). It returns a one-row tibble with the parameters, the +scenario name and the optimisation metrics (\code{n_overflows}, +\code{sum_overflows} in mm, \code{overflow_volume_m3}, water-balance shares). +} +\description{ +Factors the \code{run_one()} function that was duplicated across the three +workflow vignettes (Eisenstadt 2005, Wien, Bad Aussee) into one +package-level closure factory. The returned function runs the RAINDROP +engine for one parameter set and returns the thinned one-row +optimisation result (overflow events + water balance), augmented with +the input parameters and the overflow volume in m3. +} +\details{ +Site differences are covered by the arguments: Eisenstadt scales the +rain curve shipped in \code{base.h5} by \code{rain_factor} (leave +\code{timeseries_rain} = \code{NULL}), Wien and Bad Aussee replace the rain and +ET0 curves entirely (\code{timeseries_rain} / \code{timeseries_et}, values in +mm/h as written by the vignettes). +} +\seealso{ +\code{\link[=optimise_swale_design]{optimise_swale_design()}}, \code{\link[=find_min_feasible]{find_min_feasible()}} +} diff --git a/man/optimise_swale_design.Rd b/man/optimise_swale_design.Rd new file mode 100644 index 0000000..690f6f4 --- /dev/null +++ b/man/optimise_swale_design.Rd @@ -0,0 +1,83 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/optimise_swale_design.R +\name{optimise_swale_design} +\alias{optimise_swale_design} +\title{Find the cost-optimal swale design per overflow target} +\usage{ +optimise_swale_design( + run_fn, + x_targets = 0:5, + area_bounds = c(25, 200), + area_tol = 2, + height_bounds = c(100, 300), + height_tol = 10, + storage_spec = default_storage_spec(), + fixed = list(connected_area = 1000, filter_height = 300, filter_hydraulicconductivity = + 360, bottom_hydraulicconductivity = 12), + prior_results = NULL, + max_total_depth = NULL, + cost_rates = default_cost_rates(), + verbose = TRUE +) +} +\arguments{ +\item{run_fn}{\verb{function(params)} running one scenario and returning at +least \code{n_overflows} plus \code{sum_overflows} (mm) or \code{overflow_volume_m3}; +typically created with \code{\link[=make_swale_runner]{make_swale_runner()}}. \code{params} is a named list +of \code{mulde_area}, \code{mulde_height}, \code{storage_type}, \code{storage_height} plus +everything in \code{fixed}.} + +\item{x_targets}{Integer vector of overflow targets (feasible :<=> +\code{n_overflows <= x}), default \code{0:5}.} + +\item{area_bounds, area_tol}{Search range (m2) and resolution for +\code{mulde_area}.} + +\item{height_bounds, height_tol}{Search range (mm) and resolution for +\code{mulde_height}.} + +\item{storage_spec}{Storage search space per type, see +\code{\link[=default_storage_spec]{default_storage_spec()}}: discrete \code{levels} (infiltration box) or +continuous \code{bounds} + \code{tol} (gravel trench).} + +\item{fixed}{Named list of parameters passed unchanged to \code{run_fn} +(connected area, filter geometry, kf at maximum, ...). Must contain +\code{filter_height} for the cost model.} + +\item{prior_results}{Optional data.frame with prior (grid) results in +the workflow CSV schema, used as warm start (narrows the first area +bracket to one grid step).} + +\item{max_total_depth}{Optional analytic depth constraint in mm: +\code{mulde_height + filter_height + storage_height <= max_total_depth} +(e.g. from DWA-A 138 groundwater clearance or cover requirements). +Enforced without any simulation runs.} + +\item{cost_rates}{Unit costs, see \code{\link[=default_cost_rates]{default_cost_rates()}}.} + +\item{verbose}{Print one progress line per solved cell.} +} +\value{ +Tibble with one row per (storage type, x): the optimal design +(\code{mulde_area}, \code{mulde_height}, \code{storage_height}), its metrics +(\code{n_overflows}, \code{overflow_volume_m3}, \code{et_pct}), cost columns from +\code{\link[=compute_costs]{compute_costs()}}, a \code{status} (\code{"ok"} or \code{"infeasible_within_bounds"}), +\code{monotonicity_warning} (volume referee) and \code{n_runs_new} (fresh engine +runs spent on this cell). All evaluated designs are attached as +attribute \code{"evaluations"}. +} +\description{ +Coordinate-descent optimiser built from a single primitive +(\code{\link[=find_min_feasible]{find_min_feasible()}}, bisection over one parameter): shrink the +expensive lever first (\code{mulde_area}), then the cheap one +(\code{mulde_height}); the storage layer starts at its smallest level and is +only escalated when the area is stuck at its upper bound. The filter +conductivity is expected to be fixed at the maximum via \code{fixed} (it is +cost-free and dominant, see the \code{monotonicity_analysis} vignette). Every +engine run is cached, so the sweep over all \code{x_targets} and both storage +types shares evaluations. +} +\seealso{ +\code{\link[=find_min_feasible]{find_min_feasible()}}, \code{\link[=make_swale_runner]{make_swale_runner()}}, +\code{\link[=default_storage_spec]{default_storage_spec()}} +} diff --git a/man/read_site_timeseries.Rd b/man/read_site_timeseries.Rd new file mode 100644 index 0000000..d80aebb --- /dev/null +++ b/man/read_site_timeseries.Rd @@ -0,0 +1,35 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/read_site_timeseries.R +\name{read_site_timeseries} +\alias{read_site_timeseries} +\title{Read and prepare site rain/ET0 time series for the engine (mm/h)} +\usage{ +read_site_timeseries(path_rain, path_et, verbose = TRUE) +} +\arguments{ +\item{path_rain}{Path to the rain CSV (may be gzipped).} + +\item{path_et}{Path to the ET0 CSV (semicolon separated).} + +\item{verbose}{Print alignment messages (default TRUE).} +} +\value{ +List with data.frames \code{rain} and \code{et} (columns \code{time} = hours +since start, \code{value} = mm/h) ready for +\code{make_swale_runner(timeseries_rain = , timeseries_et = )}. +} +\description{ +Factors the time-series preparation duplicated in the Wien and Bad +Aussee workflow vignettes into one helper: reads the shipped GeoSphere +rain series (\code{rain.csv.gz}: columns \code{time} (datetime), \code{rr} (mm per +interval), \code{station}, further columns tolerated) and reference ET0 +series (\code{et.csv}: +\verb{date;value} with \code{dd.mm.yyyy}, mm per day), converts both to hours +since series start, aligns the series ends (the shorter series is +extended to the longer one's end, repeating its last value) and +converts the values to the engine's \strong{mm/h} rate convention (rain: +mm per interval / interval hours; ET0: mm per day / 24). +} +\seealso{ +\code{\link[=make_swale_runner]{make_swale_runner()}} +} diff --git a/man/sickerbox_level_presets.Rd b/man/sickerbox_level_presets.Rd new file mode 100644 index 0000000..43eafdb --- /dev/null +++ b/man/sickerbox_level_presets.Rd @@ -0,0 +1,23 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/stack_levels.R +\name{sickerbox_level_presets} +\alias{sickerbox_level_presets} +\title{Sickerbox storage-height presets (brute force default + manufacturers)} +\usage{ +sickerbox_level_presets(max_height = 2600) +} +\arguments{ +\item{max_height}{Maximum total stack height in mm passed to +\code{\link[=stack_levels]{stack_levels()}} (default 2600).} +} +\value{ +Named list of sorted numeric vectors (mm). +} +\description{ +Named list of storage-height level vectors (mm) for the infiltration-box +storage layer. \code{brute_force} is the default used by the workflow +vignettes (300/600/900/1200 mm -- itself a combination of several box +types). The manufacturer presets are generated with \code{\link[=stack_levels]{stack_levels()}} +from typical module heights of commercial block systems; verify against +the current data sheets before productive optimisation runs. +} diff --git a/man/stack_levels.Rd b/man/stack_levels.Rd new file mode 100644 index 0000000..82c00fa --- /dev/null +++ b/man/stack_levels.Rd @@ -0,0 +1,31 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/stack_levels.R +\name{stack_levels} +\alias{stack_levels} +\title{Achievable storage-layer stack heights from module heights} +\usage{ +stack_levels(modules, max_count = rep(7L, length(modules)), max_height = 2600) +} +\arguments{ +\item{modules}{Numeric vector of module heights in mm (e.g. \code{c(660, 350)} +for a full block plus a half block).} + +\item{max_count}{Integer vector (recycled to \code{length(modules)}): maximum +number of modules of each type in one stack. Defaults to 7 for every +module (cf. GRAF EcoBloc smart, stackable up to 7 layers).} + +\item{max_height}{Maximum total stack height in mm (default 2600).} +} +\value{ +Sorted numeric vector of achievable stack heights in mm. +} +\description{ +Enumerates all storage-layer heights that can be built by stacking +(and mixing) the given module heights, e.g. full blocks combined with +at most one half block. +} +\examples{ +stack_levels(360) # 360, 720, ..., 2520 +stack_levels(c(660, 350), max_count = c(7, 1)) # Rigofill full + half block + +} diff --git a/tests/testthat.R b/tests/testthat.R new file mode 100644 index 0000000..3dc1a40 --- /dev/null +++ b/tests/testthat.R @@ -0,0 +1,4 @@ +library(testthat) +library(kwb.raindrop) + +test_check("kwb.raindrop") diff --git a/tests/testthat/test-find_min_feasible.R b/tests/testthat/test-find_min_feasible.R new file mode 100644 index 0000000..1dde70e --- /dev/null +++ b/tests/testthat/test-find_min_feasible.R @@ -0,0 +1,87 @@ +test_that("kontinuierliche Bisektion findet die Schwelle innerhalb der Toleranz", { + calls <- 0L + f <- function(v) { + calls <<- calls + 1L + list(n_overflows = if (v >= 137.4) 0L else 10L) + } + res <- find_min_feasible(f, x_max = 0, lower = 25, upper = 200, tol = 2) + expect_identical(res$status, "ok") + expect_gte(res$value, 137.4) + expect_lte(res$value, 137.4 + 2) + # log2(175 / 2) ~ 6.5 -> hoechstens ~9 Laeufe inkl. Randtests + expect_lte(res$n_evaluations, 9) + # Memoisierung: jede Stelle nur einmal evaluiert + expect_equal(calls, res$n_evaluations) + expect_false(res$monotonicity_violation) +}) + +test_that("unterer Rand zulaessig -> at_lower_bound", { + f <- function(v) list(n_overflows = if (v >= 10) 0L else 10L) + res <- find_min_feasible(f, x_max = 0, lower = 25, upper = 200, tol = 2) + expect_identical(res$status, "at_lower_bound") + expect_equal(res$value, 25) +}) + +test_that("tief unzulaessiger Rand -> sofort infeasible (1 Lauf)", { + f <- function(v) list(n_overflows = 99L) + res <- find_min_feasible(f, x_max = 0, lower = 25, upper = 200, tol = 2) + expect_identical(res$status, "infeasible") + expect_true(is.na(res$value)) + expect_equal(res$n_evaluations, 1L) +}) + +test_that("Rand-Guard: +1-Zaehl-Wobble am oberen Rand frisst keine Loesung", { + # zulaessig in [100, 180), am oberen Rand springt der Zaehler auf x+1 + n_fun <- function(v) if (v >= 180) 2L else if (v >= 100) 1L else 50L + f <- function(v) list(n_overflows = n_fun(v)) + res <- find_min_feasible(f, x_max = 1, lower = 25, upper = 200, tol = 2) + expect_identical(res$status, "ok") + expect_gte(res$value, 100) + expect_lte(res$value, 102) +}) + +test_that("Volumen-Schiedsrichter: Warnung nur bei echter Nicht-Monotonie", { + # Zaehler springt bei 150 von 1 auf 2 -- Volumen steigt MIT: echte Verletzung + f_bad <- function(v) list( + n_overflows = if (v >= 150) 2L else if (v >= 100) 1L else 9L, + overflow_volume_m3 = if (v >= 150) 200 else if (v >= 100) 100 else 5000 + ) + expect_warning( + res_bad <- find_min_feasible(f_bad, x_max = 1, lower = 25, upper = 200, + tol = 2), + "non-monotonicity" + ) + expect_true(res_bad$monotonicity_violation) + + # gleicher Zaehler-Sprung, aber Volumen faellt weiter: harmloser Wobble + f_ok <- function(v) list( + n_overflows = if (v >= 150) 2L else if (v >= 100) 1L else 9L, + overflow_volume_m3 = 5000 - 20 * v + ) + expect_no_warning( + res_ok <- find_min_feasible(f_ok, x_max = 1, lower = 25, upper = 200, + tol = 2) + ) + expect_false(res_ok$monotonicity_violation) +}) + +test_that("diskrete Levels: Binaersuche ueber Stufen", { + lv <- c(300, 600, 900, 1200) + f <- function(v) list(n_overflows = if (v >= 900) 0L else 7L) + res <- find_min_feasible(f, x_max = 0, levels = lv) + expect_identical(res$status, "ok") + expect_equal(res$value, 900) + expect_lte(res$n_evaluations, 3) + + # schon die kleinste Stufe reicht + f2 <- function(v) list(n_overflows = 0L) + res2 <- find_min_feasible(f2, x_max = 0, levels = lv) + expect_identical(res2$status, "at_lower_bound") + expect_equal(res2$value, 300) + + # keine Stufe reicht (tief unzulaessig) + f3 <- function(v) list(n_overflows = 9L) + res3 <- find_min_feasible(f3, x_max = 0, levels = lv) + expect_identical(res3$status, "infeasible") + expect_equal(res3$n_evaluations, 1L) +}) diff --git a/tests/testthat/test-optimise_swale_design.R b/tests/testthat/test-optimise_swale_design.R new file mode 100644 index 0000000..f40cc23 --- /dev/null +++ b/tests/testthat/test-optimise_swale_design.R @@ -0,0 +1,129 @@ +# Synthetisches, monotones Hydraulikmodell: Ueberlaeufe fallen mit der +# Rueckhaltekapazitaet cap = Flaeche x (Muldentiefe + Porositaet x +# Speicherhoehe). Kein Engine-Aufruf noetig -> der Optimierer laesst sich +# gegen eine Brute-Force-Referenz auf feinem Raster verifizieren. +synthetic_run_factory <- function(demand) { + porosity <- c(infiltration_box = 0.95, gravel_trench = 0.3) + function(params) { + cap <- params$mulde_area * + (params$mulde_height + + porosity[[params$storage_type]] * params$storage_height) + ratio <- demand / cap + list( + n_overflows = max(0, floor(ratio) - 3), + sum_overflows = 800 * max(0, ratio - 3), + element.WB_Evapotranspiration_ = 0.1 * params$mulde_area + ) + } +} + +test_fixed <- list(connected_area = 1000, filter_height = 300, + filter_hydraulicconductivity = 360, + bottom_hydraulicconductivity = 12) + +# Brute-Force-Referenz: guenstigstes zulaessiges Design auf feinem Raster +reference_optimum <- function(run, type, x, storage_values) { + grid <- expand.grid( + mulde_area = seq(25, 200, by = 0.5), + mulde_height = seq(100, 300, by = 5), + storage_height = storage_values, + stringsAsFactors = FALSE + ) + grid$storage_type <- type + porosity <- c(infiltration_box = 0.95, gravel_trench = 0.3) + cap <- grid$mulde_area * + (grid$mulde_height + porosity[[type]] * grid$storage_height) + grid$n <- pmax(0, floor(environment(run)$demand / cap) - 3) + grid <- grid[grid$n <= x, , drop = FALSE] + if (nrow(grid) == 0) return(NULL) + grid$filter_height <- 300 + costs <- compute_costs(grid) + costs[which.min(costs$cost_total), , drop = FALSE] +} + +test_that("Optimierer findet das Kostenminimum (Vergleich mit Brute-Force)", { + run <- synthetic_run_factory(demand = 3.6e5) + out <- optimise_swale_design(run, x_targets = 0:3, + fixed = test_fixed, verbose = FALSE) + + expect_true(all(out$status == "ok")) + expect_false(any(out$monotonicity_warning)) + # Zulaessigkeit: Ueberlaufziel eingehalten + expect_true(all(out$n_overflows <= out$x)) + + for (i in seq_len(nrow(out))) { + type <- out$storage_type[i] + stor <- if (type == "infiltration_box") c(300, 600, 900, 1200) + else seq(900, 3600, by = 25) + ref <- reference_optimum(run, type, out$x[i], stor) + expect_false(is.null(ref)) + # innerhalb 5 % des (quasi-kontinuierlichen) Brute-Force-Optimums + expect_lte(out$cost_total[i], ref$cost_total * 1.05) + } + + # Kosten-Wirksamkeits-Kurve: lockereres Ziel ist nie teurer + for (type in unique(out$storage_type)) { + cc <- out$cost_total[out$storage_type == type][order(out$x[out$storage_type == type])] + expect_true(all(diff(cc) <= 1e-9)) + } + + # Laufbudget: alle 8 Zellen zusammen deutlich unter Brute-Force-Groesse + expect_lte(attr(out, "n_runs_total"), 200) +}) + +test_that("Speicher-Eskalation greift, wenn die Flaeche am Anschlag klemmt", { + run <- synthetic_run_factory(demand = 8e5) + out <- optimise_swale_design(run, x_targets = 0, + fixed = test_fixed, verbose = FALSE) + box <- out[out$storage_type == "infiltration_box", ] + expect_identical(box$status, "ok") + # bei Minimal-Speicher 300 ist selbst 200 m2 unzulaessig -> Eskalation + expect_gt(box$storage_height, 300) + expect_lte(box$n_overflows, 0) +}) + +test_that("unloesbar innerhalb der Bounds ist ein regulaeres Ergebnis", { + run <- synthetic_run_factory(demand = 5e6) + out <- optimise_swale_design(run, x_targets = 0, + fixed = test_fixed, verbose = FALSE) + expect_true(all(out$status == "infeasible_within_bounds")) + expect_true(all(is.na(out$mulde_area))) + expect_true(all(is.na(out$cost_total))) +}) + +test_that("Warmstart aus Rasterergebnissen spart Laeufe", { + run <- synthetic_run_factory(demand = 3.6e5) + + # Prior im CSV-Schema der Workflows (kf = 360, h_m = 300, Rasterschritt 25) + prior <- expand.grid( + mulde_area = seq(25, 200, by = 25), + mulde_height = 300, + storage_type = c("infiltration_box", "gravel_trench"), + stringsAsFactors = FALSE + ) + prior$storage_height <- ifelse(prior$storage_type == "infiltration_box", + 300, 900) + prior$filter_hydraulicconductivity <- 360 + prior$n_overflows <- vapply(seq_len(nrow(prior)), function(i) { + as.numeric(run(c(as.list(prior[i, ]), test_fixed))$n_overflows) + }, numeric(1)) + + cold <- optimise_swale_design(run, x_targets = 0:3, + fixed = test_fixed, verbose = FALSE) + warm <- optimise_swale_design(run, x_targets = 0:3, + fixed = test_fixed, prior_results = prior, + verbose = FALSE) + + # identisches Ergebnis, weniger Laeufe + expect_equal(warm$cost_total, cold$cost_total, tolerance = 0.02) + expect_lt(attr(warm, "n_runs_total"), attr(cold, "n_runs_total")) +}) + +test_that("max_total_depth wirkt als analytische Nebenbedingung", { + run <- synthetic_run_factory(demand = 3.6e5) + out <- optimise_swale_design(run, x_targets = 0, fixed = test_fixed, + max_total_depth = 1200, verbose = FALSE) + ok <- out[out$status == "ok", ] + expect_true(all(ok$mulde_height + ok$filter_height + ok$storage_height + <= 1200 + 1e-9)) +}) diff --git a/tests/testthat/test-stack_levels.R b/tests/testthat/test-stack_levels.R new file mode 100644 index 0000000..8ebdd63 --- /dev/null +++ b/tests/testthat/test-stack_levels.R @@ -0,0 +1,22 @@ +test_that("stack_levels: Vielfache einer Modulhoehe bis zur Obergrenze", { + expect_equal(stack_levels(360), seq(360, 2520, by = 360)) + expect_equal(stack_levels(614), c(614, 1228, 1842, 2456)) + expect_true(all(stack_levels(400) <= 2600)) +}) + +test_that("stack_levels: Mischkombination Vollblock + max. 1 Halbblock", { + expect_equal( + stack_levels(c(660, 350), max_count = c(7L, 1L)), + c(350, 660, 1010, 1320, 1670, 1980, 2330) + ) +}) + +test_that("Presets und Default-Spezifikation", { + presets <- sickerbox_level_presets() + expect_equal(presets$brute_force, c(300, 600, 900, 1200)) + expect_true(all(vapply(presets, function(p) all(diff(p) > 0), logical(1)))) + + spec <- default_storage_spec() + expect_equal(spec$infiltration_box$levels, c(300, 600, 900, 1200)) + expect_equal(spec$gravel_trench$bounds, c(900, 3600)) +}) diff --git a/vignettes/monotonicity_analysis.Rmd b/vignettes/monotonicity_analysis.Rmd new file mode 100644 index 0000000..544f0de --- /dev/null +++ b/vignettes/monotonicity_analysis.Rmd @@ -0,0 +1,455 @@ +--- +title: "Monotonie-Check der Optimierungsergebnisse" +output: rmarkdown::html_vignette +vignette: > + %\VignetteIndexEntry{Monotonie-Check der Optimierungsergebnisse} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r setup, include = FALSE, eval = TRUE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>", + fig.width = 7, + fig.height = 3.2 +) + +# Diese Vignette wertet die Ergebnis-CSVs der drei Workflow-Vignetten aus +# (workflow_badaussee, workflow_eisenstadt-2005, workflow_wien). Sie ist +# bedingt: alle Analyse-Chunks laufen nur, wenn die CSVs neben dieser +# Datei liegen -- also NACH einem lokalen Lauf der Workflows (Windows + +# Engine). Auf CI/GitHub Actions wird nur der Text gerendert. +result_files <- c( + BadAussee = "simulation_results_optimisation_BadAussee.csv", + Eisenstadt = "simulation_results_optimisation_Eisenstadt_2005.csv", + Wien = "simulation_results_optimisation_Wien.csv" +) +results_available <- all(file.exists(result_files)) + +# Auf GH Actions (pkgdown-Deploy) liegt der Ergebnisbericht (index.html) +# nicht neben dem gerenderten Artikel -- der Link darauf wird dann nicht +# gesetzt (nur beim lokalen Rendern in den Ordner monotonicity_analysis/). +is_ghactions <- tolower(Sys.getenv("GITHUB_ACTIONS")) == "true" || + tolower(Sys.getenv("CI")) %in% c("true", "1", "yes") +``` + +## Worum geht es? + +Der geplante Optimierer sucht die günstigste Muldenkonfiguration per +**Bisektion** (Intervallhalbierung): "Zu klein" → größer probieren, +"reicht" → kleiner probieren. Das funktioniert nur, wenn eine Regel gilt: + +> **Monotonie:** Wird das Bauwerk größer (mehr Fläche, mehr Tiefe, mehr +> Speicher, durchlässigerer Filter), darf es niemals *mehr* +> Überlaufereignisse geben als vorher. + +Diese Vignette prüft die Regel an allen Nachbarpaaren des Brute-Force-Rasters +(je Standort 576 Läufe, 1 704 Vergleichspaare — insgesamt 5 112). +Verglichen werden immer zwei Läufe, die sich **nur in einem Parameter um eine +Stufe** unterscheiden. Zusätzlich geprüft: das Überlaufvolumen in m³ +(`sum_overflows` ist mm Wassersäule über der Muldenfläche, daher +m³ = mm × Muldenfläche / 1000), die Verdunstung sowie zwei +Struktur-Eigenschaften, die der Optimierer nutzt (Ambiguitätsband, +Schwellentreppe, kf-Dominanz). + +## Voraussetzungen + +Diese Vignette **nach den drei Workflow-Vignetten ausführen** — sie liest +deren Ergebnisdateien (`simulation_results_optimisation_*.csv`) aus dem +`vignettes/`-Verzeichnis. + +```{r availability_note, echo = FALSE, results = 'asis', eval = !results_available} +cat(paste0( + "> **Hinweis:** Es wurden keine Ergebnisdateien gefunden — die ", + "Analyse-Chunks wurden übersprungen. Bitte zuerst die Workflow-Vignetten ", + "(Bad Aussee, Eisenstadt 2005, Wien) lokal ausführen (Windows + Engine) ", + "und diese Vignette anschließend erneut rendern.\n" +)) +``` + +### Verwendete Dateien und Ablage-Struktur + +Die Links folgen der Deploy-Struktur auf dem Server: Diese Analyse, der +Ergebnisbericht und die Detailtabellen liegen unter +`…/optimisation/monotonicity_analysis/`, die Ergebnisse der +Workflow-Vignetten unter `…/optimisation/brute-force/`. Damit die Links auch +lokal stimmen, wird diese Vignette in den gleichnamigen Unterordner +gerendert: + +```r +rmarkdown::render("monotonicity_analysis.Rmd", + output_dir = "monotonicity_analysis") +``` + +(Eingelesen werden die CSVs unverändert aus `vignettes/`, wo die Workflows +sie ablegen — nur die Links zeigen auf die Deploy-Orte.) + +| Standort | Eingangsdaten (CSV) | Interaktive Ergebnistabelle | +|----------|---------------------|-----------------------------| +| Bad Aussee | [CSV](../brute-force/simulation_results_optimisation_BadAussee.csv) | [HTML-Tabelle](../brute-force/simulation_results_optimisation_BadAussee.html) | +| Eisenstadt 2005 | [CSV](../brute-force/simulation_results_optimisation_Eisenstadt_2005.csv) | [HTML-Tabelle](../brute-force/simulation_results_optimisation_Eisenstadt_2005.html) | +| Wien | [CSV](../brute-force/simulation_results_optimisation_Wien.csv) | [HTML-Tabelle](../brute-force/simulation_results_optimisation_Wien.html) | + +```{r report_link, echo = FALSE, results = 'asis', eval = !is_ghactions} +cat(paste0( + "Der zusammenfassende, allgemeinverständliche Ergebnisbericht (Artefakt, ", + "Titel: „Monotonie-Analyse RAINDROP“) liegt als [index.html](index.html) ", + "direkt neben dieser Seite — als Verzeichnis-Startseite von ", + "`…/optimisation/monotonicity_analysis/`. Nicht zu verwechseln mit dem ", + "„Brute-Force“-Linkhub: Dessen Quelle `index.Rmd` liegt in `vignettes/`, ", + "sein gerendertes `index.html` eine Ebene über dieser Seite.\n" +)) +``` + +```{r libraries, eval = results_available, message = FALSE} +library(dplyr) +library(ggplot2) + +geom_params <- c("mulde_area", "mulde_height", "filter_hydraulicconductivity", + "storage_type", "storage_height") +check_params <- setdiff(geom_params, "storage_type") + +results <- lapply(result_files, function(f) { + readr::read_csv(f, show_col_types = FALSE) %>% + # Die Engine liefert die Ueberlaufrate in mm/h bezogen auf die + # Muldenflaeche; sum_overflows ist damit mm Wassersaeule. + # Volumen: m3 = mm x Muldenflaeche / 1000. + mutate(overflow_volume_m3 = sum_overflows * mulde_area / 1000) +}) + +# Alle Nachbarschritte entlang eines Parameters: alle uebrigen Parameter +# fixieren ("Gruppe"), nach dem Parameter sortieren, jeden Uebergang zum +# Nachbarwert als Schritt ausgeben. Verletzung :<=> dn > 0. +step_table <- function(d, v) { + grp <- setdiff(geom_params, v) + d %>% + group_by(across(all_of(grp))) %>% + arrange(.data[[v]], .by_group = TRUE) %>% + mutate(val_from = lag(.data[[v]]), + val_to = .data[[v]], + n_from = lag(n_overflows), + vol_from = lag(overflow_volume_m3), + et_from = lag(element.WB_Evapotranspiration_), + is_last = row_number() == n()) %>% + ungroup() %>% + filter(!is.na(val_from)) %>% + mutate(dn = n_overflows - n_from, + dvol = overflow_volume_m3 - vol_from, + det = element.WB_Evapotranspiration_ - et_from) +} + +steps_all <- purrr::map_dfr(names(results), function(site) { + purrr::map_dfr(check_params, function(v) { + step_table(results[[site]], v) %>% mutate(site = site, param = v) + }) +}) +``` + +## 1 Schrittverteilung: fällt, Plateau oder steigt? + +Plateaus (keine Änderung) sind unkritisch — sie sind die flachen Stufen der +Treppenfunktion und entstehen vor allem dort, wo bereits n = 0 +erreicht ist. Kritisch sind nur Anstiege ("mehr Überläufe trotz größer"). + +```{r step_distribution, eval = results_available} +step_distribution <- steps_all %>% + group_by(site, param) %>% + summarise(steps = n(), + faellt = sum(dn < 0), + plateau = sum(dn == 0), + steigt = sum(dn > 0), + steigt_anteil_pct = round(100 * mean(dn > 0), 2), + max_sprung = max(dn), + .groups = "drop") + +knitr::kable(step_distribution) +``` + +## 2 Die Verletzungen im Detail + +Jede Verletzung wird mit ihrer Volumen-Gegenprobe (in m³) gezeigt: Fällt das +Überlaufvolumen am selben Schritt weiter, ist der Anstieg des Zählers ein +Artefakt der Ereignistrennung (Pausen > 4 h teilen ein langes +Überlaufereignis in zwei gezählte), keine echte Verschlechterung. + +```{r violations, eval = results_available} +violations <- steps_all %>% + filter(dn > 0) %>% + transmute(site, param, + kontext = paste0("h_m ", mulde_height, + " | kf ", filter_hydraulicconductivity, + " | ", storage_type, + ifelse(param == "storage_height", "", + paste0(" ", storage_height))), + von = val_from, nach = val_to, + n = paste0(n_from, " → ", n_overflows), + volumen_m3 = paste0(round(vol_from, 1), " → ", + round(overflow_volume_m3, 1)), + volumen_delta_pct = round(100 * dvol / vol_from, 1), + am_rasterrand = is_last) + +DT::datatable(violations, filter = "top", + options = list(pageLength = 25, autoWidth = TRUE)) +``` + +Zur Anschauung dieselbe Bauweise an allen drei Standorten (Muldentiefe +300 mm, kf 360 mm/h, Sickerbox 300 mm): Der Zähler fällt +überall steil — nur in Bad Aussee springt er am letzten Rasterschritt von +1 auf 2 (roter Punkt), während das Volumen auch dort weiter fällt. + +```{r example_series, eval = results_available, warning = FALSE} +example_series <- purrr::map_dfr(names(results), function(site) { + results[[site]] %>% + filter(mulde_height == 300, filter_hydraulicconductivity == 360, + storage_type == "infiltration_box", storage_height == 300) %>% + arrange(mulde_area) %>% + mutate(site = site, + verletzung = n_overflows > lag(n_overflows, default = Inf)) +}) + +ggplot(example_series, aes(mulde_area, n_overflows)) + + geom_line(colour = "grey50") + + geom_point(aes(colour = verletzung), size = 2, show.legend = FALSE) + + scale_colour_manual(values = c(`FALSE` = "steelblue4", `TRUE` = "firebrick")) + + scale_y_sqrt() + + facet_wrap(~ site, scales = "free_y") + + labs(title = "Ueberlaufereignisse je Muldenflaeche (Wurzel-Skala)", + x = "Muldenflaeche [m2]", y = "Anzahl Ueberlaufereignisse") + + theme_bw() + +ggplot(example_series, aes(mulde_area, overflow_volume_m3)) + + geom_line(colour = "grey50") + + geom_point(colour = "springgreen4", size = 2) + + scale_y_sqrt() + + facet_wrap(~ site, scales = "free_y") + + labs(title = "Ueberlaufvolumen je Muldenflaeche (Wurzel-Skala): faellt ausnahmslos", + x = "Muldenflaeche [m2]", y = "Ueberlaufvolumen [m3]") + + theme_bw() +``` + +## 3 Gegenprobe: das Überlaufvolumen (m³) + +Auch in m³ gerechnet (mm Wassersäule × Muldenfläche / 1000 — entlang der +Fläche wächst der Umrechnungsfaktor mit, die Monotonie ist also nicht +automatisch übertragbar) gilt das Ergebnis: + +```{r volume, eval = results_available} +volume_check <- steps_all %>% + group_by(site) %>% + summarise(steps = n(), + volumen_steigt = sum(dvol > 1e-9), + plateau = sum(abs(dvol) <= 1e-9), + volumen_faellt = sum(dvol < -1e-9), + .groups = "drop") + +knitr::kable(volume_check) +``` + +## 4 Verdunstung: hängt nur an der Fläche + +```{r et, eval = results_available} +et_check <- steps_all %>% + group_by(site, param) %>% + summarise(steigt_pct = round(100 * mean(det > 1e-9), 1), + flach_pct = round(100 * mean(abs(det) <= 1e-9), 1), + faellt_pct = round(100 * mean(det < -1e-9), 1), + .groups = "drop") + +knitr::kable(et_check) +``` + +Konsequenz für das Sekundärziel "Verdunstung maximieren": Der Trade-off +Kosten ↔ Verdunstung verläuft eindimensional entlang der +Muldenfläche — mehr Verdunstung gibt es nur über mehr Fläche. + +## 5 Ambiguitätsband und Schwellentreppe + +Für die Bisektion relevant: Gibt es Konfigurationsgruppen, in denen oberhalb +der ersten zulässigen Fläche wieder eine unzulässige liegt (Ambiguitätsband)? +Und ist die Schwellentreppe a*(x) — die kleinste zulässige Fläche je +Überlaufziel x — monoton (ein lockereres Ziel verlangt nie mehr Fläche)? + +```{r ambiguity, eval = results_available} +grp_area <- setdiff(geom_params, "mulde_area") + +ambiguity <- purrr::map_dfr(names(results), function(site) { + purrr::map_dfr(0:5, function(x) { + results[[site]] %>% + group_by(across(all_of(grp_area))) %>% + summarise(first_ok = ifelse(any(n_overflows <= x), + min(mulde_area[n_overflows <= x]), NA), + last_bad = ifelse(any(n_overflows > x), + max(mulde_area[n_overflows > x]), NA), + .groups = "drop") %>% + summarise(site = site, x = x, + gruppen = n(), + mit_loesung = sum(!is.na(first_ok)), + ambig = sum(!is.na(first_ok) & !is.na(last_bad) & + last_bad > first_ok), + max_band_m2 = max(c(0, (last_bad - first_ok)[ + !is.na(first_ok) & !is.na(last_bad)]), na.rm = TRUE)) + }) +}) + +knitr::kable(ambiguity) + +staircase <- purrr::map_dfr(names(results), function(site) { + results[[site]] %>% + group_by(across(all_of(grp_area))) %>% + reframe(x = 0:5, + a_star = sapply(0:5, function(x) + ifelse(any(n_overflows <= x), + min(mulde_area[n_overflows <= x]), NA))) %>% + group_by(across(all_of(grp_area))) %>% + summarise(treppen_verletzungen = { + a <- a_star[order(x)] + a <- a[!is.na(a)] + if (length(a) > 1) sum(diff(a) > 0) else 0L + }, .groups = "drop") %>% + summarise(site = site, branches = n(), + treppen_verletzungen = sum(treppen_verletzungen)) +}) + +knitr::kable(staircase) +``` + +## 6 kf-Dominanz: der Filter ist ein Gratis-Hebel + +Die Filterdurchlässigkeit kostet nichts (sie taucht in `compute_costs()` +nicht auf). Wenn sie gleichzeitig die Verdunstung nicht verändert und die +Überläufe nie erhöht, kann der Optimierer sie fest auf das Maximum setzen. + +```{r kf_dominance, eval = results_available} +kf_dominance <- purrr::map_dfr(names(results), function(site) { + results[[site]] %>% + group_by(filter_hydraulicconductivity) %>% + summarise(site = site, + mittlere_ET_pct = round(mean(element.WB_Evapotranspiration_), 2), + mittlere_n_overflows = round(mean(n_overflows), 1), + mittleres_volumen_m3 = round(mean(overflow_volume_m3), 1), + .groups = "drop") +}) + +knitr::kable(kf_dominance) +``` + +## 7 Warmstart: billigstes zulässiges Raster-Design je Ziel x + +Diese Tabelle ist der Startpunkt der Bisektion: das jeweils günstigste +zulässige Design aus dem vorhandenen Raster (bei kf = Maximum), je +Überlaufziel x und Speichertyp. Der Optimierer verfeinert nur noch im +25-m²-Bracket darunter. + +```{r warmstart, eval = results_available} +warmstart <- purrr::map_dfr(names(results), function(site) { + d <- results[[site]] %>% + filter(filter_hydraulicconductivity == max(filter_hydraulicconductivity)) + purrr::map_dfr(0:5, function(x) { + d %>% + filter(n_overflows <= x) %>% + group_by(storage_type) %>% + slice_min(cost_total, n = 1, with_ties = FALSE) %>% + ungroup() %>% + transmute(site = site, x = x, storage_type, mulde_area, mulde_height, + storage_height, cost_total, + overflow_volume_m3 = round(overflow_volume_m3, 1), + ET_pct = round(element.WB_Evapotranspiration_, 1)) + }) +}) + +DT::datatable(warmstart, filter = "top", + options = list(pageLength = 12, autoWidth = TRUE)) +``` + +## Fazit: Regeln für den Optimierer + +Stand der letzten vollständigen Auswertung (2026-07-10): 5 112 +Vergleichspaare, 13 Verletzungen (0,25 %), alle mit Sprunghöhe genau +1 +und fallendem Volumen — Eisenstadt 0, Wien 1 (bei n ≈ 287, +irrelevant für x ≤ 5), Bad Aussee 12 (alle am Rasterrand +175 → 200 m², Niveau n = 1). Das +Überlaufvolumen stieg in keinem einzigen Vergleich. **Die +Monotonie-Voraussetzung der Bisektion ist damit erfüllt**, abgesichert durch +drei Regeln: + +1. **Grid-Warmstart:** Bisektion nur im 25-m²-Bracket um das bekannte + `first_ok` aus dem Raster verfeinern — dort ist die Grenze an allen drei + Standorten eindeutig. +2. **Rand-Guard statt Abbruch:** Fällt der obere Intervallrand nur um +1 + über das Ziel (n = x + 1), erst Cache-/Rasterpunkte darunter + prüfen, bevor "keine Lösung" gemeldet wird. +3. **Volumen als Schiedsrichter:** Bei jedem nicht-monotonen Flip prüfen, ob + `sum_overflows` weiter gefallen ist. Ja → bekannter Zähl-Wobble, + weiterrechnen. Nein → Warnung (bisher nie aufgetreten) — der eingebaute + Rauchmelder für künftige Standorte. + +Beim Hinzufügen eines neuen Standorts: Workflows laufen lassen und diese +Vignette erneut rendern. Zeigt Abschnitt 1 Sprünge > +1 oder Abschnitt 3 +steigendes Volumen, ist die Bisektion für diesen Standort nicht abgesichert. + +```{r export, eval = results_available, message = FALSE} +save_table_html <- function(df, file, title) { + htmlwidgets::saveWidget( + DT::datatable(df, filter = "top", + options = list(pageLength = 25, autoWidth = TRUE)), + file = file, selfcontained = TRUE, title = title + ) + # saveWidget() laesst bei Zielpfaden ausserhalb des Arbeitsverzeichnisses + # den "_files"-lib-Ordner stehen, obwohl die Datei selfcontained + # ist -> redundante Duplikate aufraeumen. + unlink(sub("\\.html$", "_files", file), recursive = TRUE) +} + +exports <- list( + mono_step_distribution = list(df = step_distribution, + titel = "Monotonie: Schrittverteilung"), + mono_violations_detail = list(df = violations, + titel = "Monotonie: Verletzungen im Detail"), + mono_volume = list(df = volume_check, + titel = "Monotonie: Ueberlaufvolumen (m3)"), + mono_et = list(df = et_check, + titel = "Monotonie: Verdunstung"), + mono_ambiguity = list(df = ambiguity, + titel = "Ambiguitaetsband je Ziel x"), + mono_staircase = list(df = staircase, + titel = "Schwellentreppe a*(x)"), + mono_kf_dominance = list(df = kf_dominance, + titel = "kf-Dominanz"), + mono_warmstart_designs = list(df = warmstart, + titel = "Warmstart-Designs") +) + +# Exporte in den Unterordner des Ergebnisberichts (index.html), damit der +# gesamte Ordner monotonicity_analysis/ als eine Einheit deploybar ist. +out_dir <- "monotonicity_analysis" +dir.create(out_dir, showWarnings = FALSE) + +for (name in names(exports)) { + readr::write_csv(exports[[name]]$df, file.path(out_dir, paste0(name, ".csv"))) + save_table_html(exports[[name]]$df, file.path(out_dir, paste0(name, ".html")), + title = exports[[name]]$titel) +} + +# Den Rmd-Quellcode mit in den Deploy-Ordner kopieren, damit der Link +# "monotonicity_analysis.Rmd" des Ergebnisberichts (index.html) dort +# funktioniert. +invisible(file.copy("monotonicity_analysis.Rmd", + file.path(out_dir, "monotonicity_analysis.Rmd"), + overwrite = TRUE)) +``` + +Die exportierten Detailtabellen liegen neben dieser Seite (im Unterordner +`monotonicity_analysis/`), jeweils als CSV und als interaktive +HTML-Tabelle: + +| Tabelle | Inhalt | CSV | HTML | +|---------|--------|-----|------| +| Schrittverteilung | fällt / Plateau / steigt je Standort × Parameter | [CSV](mono_step_distribution.csv) | [HTML](mono_step_distribution.html) | +| Verletzungen | alle 13 Fälle mit Volumen-Gegenprobe (m³) | [CSV](mono_violations_detail.csv) | [HTML](mono_violations_detail.html) | +| Überlaufvolumen | Monotonie des Volumens in m³ | [CSV](mono_volume.csv) | [HTML](mono_volume.html) | +| Verdunstung | ET-Richtung je Parameter | [CSV](mono_et.csv) | [HTML](mono_et.html) | +| Ambiguitätsband | Eindeutigkeit der Zulässigkeitsgrenze je Ziel x | [CSV](mono_ambiguity.csv) | [HTML](mono_ambiguity.html) | +| Schwellentreppe | Monotonie von a*(x) je Branch | [CSV](mono_staircase.csv) | [HTML](mono_staircase.html) | +| kf-Dominanz | ET / Überläufe / Volumen je kf-Stufe | [CSV](mono_kf_dominance.csv) | [HTML](mono_kf_dominance.html) | +| Warmstart-Designs | billigstes zulässiges Raster-Design je x und Speichertyp | [CSV](mono_warmstart_designs.csv) | [HTML](mono_warmstart_designs.html) | diff --git a/vignettes/workflow_optimisation.Rmd b/vignettes/workflow_optimisation.Rmd new file mode 100644 index 0000000..6c184d3 --- /dev/null +++ b/vignettes/workflow_optimisation.Rmd @@ -0,0 +1,293 @@ +--- +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 +``` + +## 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); der komplette Workflow +liegt damit bei ca. 30–40 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 + +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} +opt_all <- purrr::map_dfr(names(sites), function(site) { + cfg <- sites[[site]] + message("=== ", site, " ===") + + ts <- if (cfg$timeseries) { + read_site_timeseries( + extdata_path("models", cfg$dir, "rain.csv.gz"), + extdata_path("models", cfg$dir, "et.csv") + ) + } 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, + fixed = fixed, + prior_results = prior, + cost_rates = cost_rates, + verbose = TRUE) + opt$site <- site + opt$n_runs_site <- attr(opt, "n_runs_total") + opt$minutes_site <- round(as.numeric( + difftime(Sys.time(), t0, units = "mins")), 1) + message(sprintf("%s fertig: %d Engine-Laeufe in %.1f min", + site, opt$n_runs_site[1], opt$minutes_site[1])) + opt +}) +``` + +## 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)) +} + +knitr::kable(unique(opt_all[, c("site", "n_runs_site", "minutes_site")]), + col.names = c("Standort", "Engine-Laeufe", "Minuten")) +``` + +## 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.