## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
  comment = "#>",
  fig.width = 9,
  fig.asp = 0.68,
  out.width = "90%",
  fig.align = "center"
)

## ----setup-lib----------------------------------------------------------------
library(airGRiwrm)
library(lubridate)

## ----load-data----------------------------------------------------------------
data(Severn)
Severn$BasinsInfo

## ----build-nodes--------------------------------------------------------------
# Extract base node information from Severn data
nodes_Severn <- Severn$BasinsInfo[, c(
  "gauge_id",
  "downstream_id",
  "distance_downstream",
  "area"
)]
names(nodes_Severn) <- c("id", "down", "length", "area")
nodes_Severn$model <- "RunModel_GR4J"

# Reduce to the two nodes of interest: 54095 (upstream) and 54001 (downstream)
nodes <- nodes_Severn[nodes_Severn$id %in% c("54095", "54001"), ]

# Node 54001 becomes the outlet -- no further downstream connection
nodes$down[nodes$id == "54001"] <- NA
nodes$length[nodes$id == "54001"] <- NA

# Add the Reservoir row: dam built directly on the river at 0 km downstream of 54095.
# The RunModel_GR4J row for 54095 now feeds into the Reservoir instead of 54001.
nodes <- rbind(
  nodes,
  data.frame(
    id = "Reservoir",
    down = "54001",
    length = 42,
    area = NA,
    model = "RunModel_Reservoir",
    stringsAsFactors = FALSE
  )
)

# Update 54095 to feed into the Reservoir (dam built at 0 km downstream)
nodes$down[nodes$id == "54095"] <- "Reservoir"
nodes$length[nodes$id == "54095"] <- 0

nodes

## ----create-griwrm------------------------------------------------------------
griwrm <- CreateGRiwrm(nodes)
griwrm

## ----plot-network, fig.asp = 0.7----------------------------------------------
plot(griwrm)

## ----prepare-meteo------------------------------------------------------------
BasinsObs <- Severn$BasinsObs
DatesR <- BasinsObs[[1]]$DatesR

PrecipTot <- cbind(sapply(BasinsObs, function(x) x$precipitation))
PotEvapTot <- cbind(sapply(BasinsObs, function(x) x$peti))
Qobs <- cbind(sapply(BasinsObs, function(x) x$discharge_spec))

# Convert to sub-catchment scale
Precip <- ConvertMeteoSD(griwrm, PrecipTot)
PotEvap <- ConvertMeteoSD(griwrm, PotEvapTot)

## ----prepare-quantities-------------------------------------------------------
Qrelease <- data.frame(Reservoir = rep(1E15, length(DatesR)))

## ----create-inputs------------------------------------------------------------
InputsModel <- suppressWarnings(
  CreateInputsModel(
    griwrm,
    DatesR,
    Precip,
    PotEvap,
    Qrelease = Qrelease
  )
)
str(InputsModel, max.level = 1)

## ----run-options--------------------------------------------------------------
date_start <- "2002-01-01"
date_end <- "2003-12-31"
IndPeriod_Run <- seq(
  which(format(InputsModel[[1]]$DatesR, format = "%Y-%m-%d") == date_start),
  which(format(InputsModel[[1]]$DatesR, format = "%Y-%m-%d") == date_end)
)
IndPeriod_WarmUp <- seq(1, IndPeriod_Run[1] - 1)

## ----params-------------------------------------------------------------------
Param <- list(
  `54095` = c(256.84415, 0.065046, 57.52368, 2.718095),
  `54001` = c(1.03, 24.77909, -1.904302, 21.75840, 1.378378),
  Reservoir = c(100E6, 1)
)

## ----filling-curve------------------------------------------------------------
curve <- approx(
  x = c(31 * 11 - 365, 30 * 6, 31 * 11, 366 + 30 * 6),
  y = c(20E6, 80E6, 20E6, 80E6),
  xout = 1:366
)$y

# Visualise the filling curve
plot(
  1:366,
  curve / 1E6,
  type = "l",
  lwd = 2,
  xlab = "Day of Year",
  ylab = expression("Target Volume (10"^6 * " m"^3 * ")"),
  main = "Annual Filling Curve",
  col = "steelblue"
)
abline(h = 100, lty = 2, col = "red") # Vmax reference
abline(h = 80, lty = 3, col = "darkgreen")
abline(h = 20, lty = 3, col = "darkgreen")
legend(
  "topright",
  legend = c("Target curve", "Vmax = 100 hm3", "Operating band"),
  col = c("steelblue", "red", "darkgreen"),
  lty = c(1, 2, 3),
  lwd = c(2, 1, 1)
)

## ----controller-factory-------------------------------------------------------
fn_filling_curve_factory <- function(sv, curve, QrelMax, QrelMin) {
  # Shared state across controller calls
  if (is.null(sv$state)) {
    sv$state <- list(prevError = 0, release_fraction = 1.0)
  }

  # PD gains
  Kp <- 1
  Kd <- 0.1

  function(Y) {
    Vsim <- Y[1] + Y[2] # Current simulated volume with inflows
    j <- as.numeric(format(sv$ts.date, "%j")) # Day of year
    Vobj_ts <- approx(x = 1:366, y = curve, xout = j)$y # Target for today

    error <- Vsim - Vobj_ts
    dError <- error - sv$state$prevError

    # Step 1: Compute Qrelease from PD controller
    ctrl_action <- Kp * error + Kd * dError
    sv$state$prevError <- error
    Qrelease <- ctrl_action

    # Step 2: Apply hard constraints -- clamp between scaled minimum and maximum
    frac <- sv$state$release_fraction
    Qrelease <- max(frac * QrelMin, min(QrelMax, Qrelease))

    # Return a single-column matrix with the release command
    return(cbind(Qrelease))
  }
}

## ----supervisor---------------------------------------------------------------
sv <- CreateSupervisor(InputsModel)

## ----controller-fun-----------------------------------------------------------
fn_filling_curve <- fn_filling_curve_factory(
  sv = sv,
  curve = curve,
  QrelMax = 200 * 86400, # Maximum release rate (m3/day)
  QrelMin = 30 * 86400 # Mandatory minimum release rate (m3/day)
)

## ----controller---------------------------------------------------------------
CreateController(
  sv,
  ctrl.id = "FillingCurve",
  Y = c("Reservoir$Vsim", "54095"), # Read current reservoir volume with inflows
  U = "Reservoir$Qrelease", # Release command only
  FUN = fn_filling_curve # PD control logic
)

## ----time-loop----------------------------------------------------------------
# Define simulation period: from 01/01/2000 up to 31/12/2009
dfTS <- data.frame(
  DatesR = DatesR,
  yearmonth = format(DatesR, "%Y-%m")
)
tactical_months <- unique(dfTS$yearmonth[
  dfTS$yearmonth >= format(as.Date(date_start), "%Y-%m") &
    dfTS$yearmonth <= format(as.Date(date_end), "%Y-%m")
])

first_run <- TRUE

# Store the tactical release fraction; default to 1 (unrestricted)
release_fraction <- 1.0

for (ym in tactical_months) {
  message(sprintf("Processing period: %s", ym), appendLF = FALSE)

  # Extract the current month's indices
  ym_IndPeriod_Run <- which(dfTS$yearmonth == ym)

  if (first_run) {
    message(" -> Release fraction = 1.0 (initial run)")
    RunOptions <- CreateRunOptions(
      InputsModel,
      IndPeriod_WarmUp = IndPeriod_WarmUp,
      IndPeriod_Run = ym_IndPeriod_Run,
      IniStates = list(Reservoir = c("Reservoir.V" = 0)),
      warnings = FALSE
    )
    OM <- suppressMessages(
      RunModel(
        sv,
        RunOptions = RunOptions,
        Param = Param
      )
    )
    first_run <- FALSE
  } else {
    # --- Tactical rule: compute release fraction from volume ratio ---
    V_current <- OM$Reservoir$StateEnd$Reservoir$V

    # Get the day-of-year for the first day of this month
    first_day_ym <- as.Date(paste0(ym, "-01"))
    j_first <- as.numeric(format(first_day_ym, "%j"))
    V_objective <- approx(x = 1:366, y = curve, xout = j_first)$y

    # Compute the release fraction: ratio of current to objective volume
    # When V_current < V_objective, fraction < 1, reducing the minimum release
    # When V_current >= V_objective, fraction = 1, full minimum applies
    release_fraction <- min(1, V_current / V_objective)

    # Update the supervisor state so the controller reads the new fraction
    sv$state$release_fraction <- release_fraction

    message(sprintf(
      " -> Release fraction: %.2f (V_cur=%.1f, V_obj=%.1f)",
      release_fraction,
      V_current / 1E6,
      V_objective / 1E6
    ))

    # Run the model for the current month with supervisor control
    OM <- suppressWarnings(
      RunModel(
        OM,
        InputsModel = sv,
        RunOptions = RunOptions,
        IndPeriod_Run = ym_IndPeriod_Run
      )
    )
  }
}

## ----plot-results, fig.asp = 0.7----------------------------------------------
# Expand the annual filling curve to cover the full simulation period
# by repeating it for each year
ydoys <- as.numeric(format(as.Date(OM[[1]]$DatesR), "%j"))
curve_full <- approx(x = 1:366, y = curve, xout = ydoys)$y
# Plot Reservoir volume, inflows, and release
plot(OM$Reservoir, Vobs = curve_full)

## ----plot-downstream, fig.asp = 0.6-------------------------------------------
plot(attr(OM, "Qm3s"))