Resolve "Feature request: feedback control"

Closes #19 (closed)

R/CreateSupervisor.R

@@ -24,6 +24,7 @@
 #' InputsModel <- CreateInputsModel(griwrm, DatesR, Precip, PotEvap, Qobs)
 #' sv <- CreateSupervisor(InputsModel)
 CreateSupervisor <- function(InputsModel) {
+  if(!inherits(InputsModel, "GRiwrmInputsModel")) stop("`InputsModel` parameter must of class 'GRiwrmInputsModel' (See ?CreateInputsModel.GRiwrm)")
   # Create Supervisor environment in the parent of GlobalEnv
   e <- new.env(parent = parent.env(globalenv()))
   class(e) <- c("Supervisor", class(e))

R/RunModel.GR.R

@@ -7,7 +7,6 @@
 #' @export
 #'
 RunModel.GR <- function(x, RunOptions, Param, ...) {
-  message("RunModel.GR")
 
   if (inherits(x, "SD")) {
     # Lag model take one parameter at the beginning of the vector

R/RunModel.GRiwrmInputsModel.R

@@ -8,32 +8,25 @@
 #' @return \emph{GRiwrmOutputsModel} object which is a list of \emph{OutputsModel} objects (See \code{\link[airGR]{RunModel}}) for each node of the semi-distributed model.
 #' @export
 RunModel.GRiwrmInputsModel <- function(x, RunOptions, Param, ...) {
-  message("RunModel.GRiwrmInputsModel")
 
-  # Run runoff model for each sub-basin
-  OutputsModel <- lapply(X = x, FUN = function(IM) {
-    RunModel.GR(IM,
-                RunOptions = RunOptions[[IM$id]],
-                Param = Param[[IM$id]])
-    })
-  class(OutputsModel) <- append(class(OutputsModel), "GRiwrmOutputsModel")
+  checkRunModelParameters(x, RunOptions, Param)
 
-  # Loop over sub-basin using SD model
-  for(id in getSD_Ids(x)) {
-    IM <- x[[id]]
-    message("RunModel.GRiwrmInputsModel: Treating sub-basin ", id, "...")
+  OutputsModel <- list()
+  class(OutputsModel) <- c("GRiwrmOutputsModel", class(OutputsModel))
+
+  for(IM in x) {
+    message("RunModel.GRiwrmInputsModel: Treating sub-basin ", IM$id, "...")
 
     # Update x$Qupstream with simulated upstream flows
     if(any(IM$UpstreamIsRunoff)) {
-      IM <- UpdateQsimUpstream(IM, RunOptions[[id]]$IndPeriod_Run, OutputsModel)
+      IM <- UpdateQsimUpstream(IM, RunOptions[[IM$id]]$IndPeriod_Run, OutputsModel)
     }
 
-    # Run the SD model for the sub-basin
-    OutputsModel[[id]] <- RunModel.SD(
+    # Run the model for the sub-basin
+    OutputsModel[[IM$id]] <- RunModel.InputsModel(
       IM,
-      RunOptions = RunOptions[[id]],
-      Param = Param[[id]],
-      OutputsModel[[id]]
+      RunOptions = RunOptions[[IM$id]],
+      Param = Param[[IM$id]]
     )
 
   }

R/RunModel.R

@@ -6,6 +6,5 @@
 #' @return Either a [list] of OutputsModel object (for GRiwrmInputsModel) or an OutputsModel object (for InputsModel)
 #' @export
 RunModel <- function(x, ...) {
-  message("RunModel")
   UseMethod("RunModel", x)
 }

R/RunModel.SD.R

 #' Run SD Model from run-off model outputs
 #'
 #' @inheritParams airGR::RunModel_Lag
-#' @param x `InputsModel` used as `InputsModel` parameter for [airGR::RunModel]
-#' @param OutputsModel `OutputsModel` object returned by a GR model by [airGR::RunModel]
+#' @param x `InputsModel` used as `InputsModel` parameter for [airGR::RunModel_Lag]
+#' @param QSimDown a [numeric] corresponding to the runoff of the sub-basin (Typically the `Qsim` outputs of the GR model)
 #' @param ... further arguments passed to or from other methods
 #'
 #' @return `OutputsModel` object. See [airGR::RunModel_Lag]
 #' @export
 #'
-RunModel.SD <- function(x, RunOptions, Param, OutputsModel, ...) {
-  message("RunModel.SD")
-  x$OutputsModel <- OutputsModel
+RunModel.SD <- function(x, RunOptions, Param, QsimDown, ...) {
+  x$OutputsModel <- list(Qsim = QsimDown)
   RunModel_Lag(x, RunOptions = RunOptions, Param = Param[1])
 }

R/RunModel.Supervisor.R

@@ -14,32 +14,53 @@ RunModel.Supervisor <- function(x, RunOptions, Param, ...) {
   })
 
   # Run runoff model for each sub-basin
-  OutputsModel <- lapply(X = x$InputsModel, FUN = function(IM) {
+  x$OutputsModel <- lapply(X = x$InputsModel, FUN = function(IM) {
     RunModel.GR(IM,
                 RunOptions = RunOptions[[IM$id]],
                 Param = Param[[IM$id]])
     })
-  class(OutputsModel) <- append(class(OutputsModel), "GRiwrmOutputsModel")
+  class(x$OutputsModel) <- append(class(x$OutputsModel), "GRiwrmOutputsModel")
+  # Save Qsim for step by step simulation
+  Qsim <- lapply(x$OutputsModel, function(OM) {
+    OM$Qsim
+  })
 
-  # Loop over time steps
-  # Loop over sub-basin using SD model
+  # Adapt RunOptions to step by step simulation
   for(id in getSD_Ids(x$InputsModel)) {
-    IM <- x$InputsModel[[id]]
-    message("RunModel.GRiwrmInputsModel: Treating sub-basin ", id, "...")
+    RunOptions[[id]]$IndPeriod_WarmUp <- 0L
+    RunOptions[[id]]$Outputs_Sim <- "StateEnd"
+  }
 
-    # Update InputsModel$Qupstream with simulated upstream flows
-    if(any(IM$UpstreamIsRunoff)) {
-      IM <- UpdateQsimUpstream(IM, RunOptions[[id]]$IndPeriod_Run, OutputsModel)
-    }
+  # Loop over time steps
+  for(iTS in RunOptions[[1]]$IndPeriod_Run) {
+    # Run regulation on the whole basin for the current time step
+    x$ts.index <- iTS
+    x$ts.date <- x$InputsModel[[1]]$DatesR[iTS]
+    doSupervision(x)
 
-    # Run the SD model for the sub-basin
-    OutputsModel[[id]] <- RunModel.SD(
-      IM,
-      RunOptions = RunOptions[[id]],
-      Param = Param[[id]],
-      OutputsModel[[id]]
-    )
+    # Loop over sub-basin using SD model
+    for(id in getSD_Ids(x$InputsModel)) {
 
+      # Update InputsModel$Qupstream with simulated upstream flows
+      for(i in which(x$InputsModel[[id]]$UpstreamIsRunoff)) {
+        x$InputsModel[[id]]$Qupstream[iTS, i] <-
+          x$OutputsModel[[x$InputsModel[[id]]$UpstreamNodes[i]]]$Qsim[iTS - x$ts.index0]
+      }
+
+      # Run the SD model for the sub-basin and one time step
+      RunOptions[[id]]$IndPeriod_Run <- iTS
+      RunOptions[[id]]$IniStates <- unlist(x$OutputsModel[[id]]$StateEnd)
+      x$OutputsModel[[id]] <- RunModel.SD(
+        x$InputsModel[[id]],
+        RunOptions = RunOptions[[id]],
+        Param = Param[[id]],
+        QsimDown = Qsim[[id]][iTS - x$ts.index0]
+      )
+      Qsim[[id]][iTS - x$ts.index0] <- x$OutputsModel[[id]]$Qsim
+    }
+  }
+  for(id in getSD_Ids(x$InputsModel)) {
+    x$OutputsModel[[id]]$Qsim <- Qsim[[id]]
   }
-  return(OutputsModel)
+  return(x$OutputsModel)
 }

R/RunModel_Lag.R

+RunModel_Lag <- function(InputsModel, RunOptions, Param) {
+  NParam <- 1
+
+  ##Arguments_check
+  if (!inherits(InputsModel, "InputsModel")) {
+    stop("'InputsModel' must be of class 'InputsModel'")
+  }
+  if (!inherits(InputsModel, "SD")) {
+    stop("'InputsModel' must be of class 'SD'")
+  }
+  if (!inherits(RunOptions, "RunOptions")) {
+    stop("'RunOptions' must be of class 'RunOptions'")
+  }
+  if (!is.vector(Param) | !is.numeric(Param)) {
+    stop("'Param' must be a numeric vector")
+  }
+  if (sum(!is.na(Param)) != NParam) {
+    stop(paste("'Param' must be a vector of length", NParam, "and contain no NA"))
+  }
+  if (is.null(InputsModel$OutputsModel)) {
+    stop(
+      "'InputsModel' should contain an 'OutputsModel' key containing the output of the runoff of the downstream subcatchment"
+    )
+  }
+  if (is.null(InputsModel$OutputsModel$Qsim)) {
+    stop(
+      "'InputsModel$OutputsModel' should contain a key 'Qsim' containing the output of the runoff of the downstream subcatchment"
+    )
+  }
+  if (sum(!is.na(InputsModel$OutputsModel$Qsim)) != length(RunOptions$IndPeriod_Run)) {
+    stop(
+      "'InputsModel$OutputsModel$Qim' should have the same lenght as 'RunOptions$IndPeriod_Run' and contain no NA"
+    )
+  }
+
+  OutputsModel <- InputsModel$OutputsModel
+  OutputsModel$QsimDown <- OutputsModel$Qsim
+
+  if (inherits(InputsModel, "hourly")) {
+    TimeStep <- 60 * 60
+  } else if (inherits(InputsModel, "daily")) {
+    TimeStep <- 60 * 60 * 24
+  } else {
+    stop("'InputsModel' should be of class \"daily\" or \"hourly\"")
+  }
+
+  # propagation time from upstream meshes to outlet
+  PT <- InputsModel$LengthHydro / Param[1L] / TimeStep
+  HUTRANS <- rbind(1 - (PT - floor(PT)), PT - floor(PT))
+
+  NbUpBasins <- length(InputsModel$LengthHydro)
+  LengthTs <- length(OutputsModel$QsimDown)
+  OutputsModel$Qsim <- OutputsModel$QsimDown * InputsModel$BasinAreas[length(InputsModel$BasinAreas)] * 1e3
+
+  IniSD <- RunOptions$IniStates[grep("SD", names(RunOptions$IniStates))]
+  if (length(IniSD) > 0) {
+    if (sum(floor(PT)) + NbUpBasins != length(IniSD)) {
+      stop(
+        sprintf(
+          "SD initial states has a length of %i and a length of %i is required",
+          length(IniSD),
+          sum(floor(PT)) + NbUpBasins
+        )
+      )
+    }
+    IniStates <- lapply(seq_len(NbUpBasins), function(x) {
+      iStart <- 1
+      if (x > 1) {
+        iStart <- iStart + sum(floor(PT[1:x - 1]) + 1)
+      }
+      IniSD[iStart:(iStart + PT[x])]
+    })
+  } else {
+    IniStates <- lapply(seq_len(NbUpBasins), function(x) {
+      rep(0, floor(PT[x] + 1))
+    })
+  }
+
+  for (upstream_basin in seq_len(NbUpBasins)) {
+    Qupstream <- InputsModel$Qupstream[RunOptions$IndPeriod_Run, upstream_basin]
+    if (!is.na(InputsModel$BasinAreas[upstream_basin])) {
+      # Upstream flow with area needs to be converted to m3 by time step
+      Qupstream <- Qupstream * InputsModel$BasinAreas[upstream_basin] * 1e3
+    }
+    Qupstream1 <- c(IniStates[[upstream_basin]][-length(IniStates[[upstream_basin]])], Qupstream[1:(LengthTs - floor(PT[upstream_basin]))])
+    Qupstream2 <- IniStates[[upstream_basin]]
+    if(LengthTs - floor(PT[upstream_basin]) - 1 > 0) Qupstream2 <- c(Qupstream2, Qupstream[1:(LengthTs - floor(PT[upstream_basin]) - 1)])
+    OutputsModel$Qsim <- OutputsModel$Qsim +
+                         Qupstream1 * HUTRANS[1, upstream_basin] +
+                         Qupstream2 * HUTRANS[2, upstream_basin]
+  }
+  # Warning for negative flows
+  if (any(OutputsModel$Qsim < 0)) {
+    warning(length(which(OutputsModel$Qsim < 0)), " time steps with negative flow, set to zero.")
+    OutputsModel$Qsim[OutputsModel$Qsim < 0] <- 0
+  }
+  # Convert back Qsim to mm
+  OutputsModel$Qsim <- OutputsModel$Qsim / sum(InputsModel$BasinAreas, na.rm = TRUE) / 1e3
+
+  if ("StateEnd" %in% RunOptions$Outputs_Sim) {
+    OutputsModel$StateEnd$SD <- lapply(seq(NbUpBasins), function(x) {
+      LengthTs <- tail(RunOptions$IndPeriod_Run,1)
+      InputsModel$Qupstream[(LengthTs - floor(PT[x])):LengthTs, x]
+    })
+  }
+
+  return(OutputsModel)
+}

R/utils.R

@@ -53,7 +53,6 @@ setDataToLocation <- function(control, supervisor) {
 #' @return [NULL]
 doSupervision <- function(supervisor) {
   for (id in names(supervisor$controllers)) {
-    ctrlr <- supervisor$controllers[[id]]
     # Read Y from locations in the model
     supervisor$controllers[[id]]$Y$v <-
       sapply(supervisor$controllers[[id]]$Y$loc, getDataFromLocation, supervisor = supervisor)
@@ -65,3 +64,19 @@ doSupervision <- function(supervisor) {
   }
   return()
 }
+
+#' Check the parameters of RunModel methods
+#'
+#' Stop the execution if an error is detected.
+#'
+#' @param RunOptions a `GRiwrmRunOptions` object (See [CreateRunOptions.GRiwrmInputsModel])
+#' @param Param a [list] of [numeric] containing model parameters of each node of the network
+#'
+#' @return [NULL]
+#'
+checkRunModelParameters <- function(InputsModel, RunOptions, Param) {
+  if(!inherits(InputsModel, "GRiwrmInputsModel")) stop("`InputsModel` parameter must of class 'GRiwrmRunoptions' (See ?CreateRunOptions.GRiwrmInputsModel)")
+  if(!inherits(RunOptions, "GRiwrmRunOptions")) stop("Argument `RunOptions` parameter must of class 'GRiwrmRunOptions' (See ?CreateRunOptions.GRiwrmInputsModel)")
+  if(!is.list(Param) || !all(names(InputsModel) %in% names(Param))) stop("Argument `Param` must be a list with names equal to nodes IDs")
+  return()
+}

tests/testthat/test-RunModel.R

+context("RunModel.Supervisor")
+
+test_that("RunModelSupervisor with no regulation should returns same results as RunModel.GRiwrmInputsModel", {
+  data(Severn)
+  # Network configuration
+  nodes <- Severn$BasinsInfo[c(1,2,5), c("gauge_id", "downstream_id", "distance_downstream", "area")]
+  nodes$distance_downstream <- nodes$distance_downstream * 1000 # Conversion km -> m
+  nodes$model <- NA
+  nodes$model[1] <- "RunModel_GR4J"
+  griwrm <- GRiwrm(nodes, list(id = "gauge_id", down = "downstream_id", length = "distance_downstream"))
+  # InputsModel
+  DatesR <- Severn$BasinsObs[[1]]$DatesR
+  PrecipTot <- cbind(sapply(Severn$BasinsObs, function(x) {x$precipitation}))
+  PotEvapTot <- cbind(sapply(Severn$BasinsObs, function(x) {x$peti}))
+  Precip <- ConvertMeteoSD(griwrm, PrecipTot)
+  PotEvap <- ConvertMeteoSD(griwrm, PotEvapTot)
+  Qobs <- cbind(sapply(Severn$BasinsObs, function(x) {x$discharge_spec}))
+  InputsModel <- CreateInputsModel(griwrm, DatesR, Precip, PotEvap, Qobs)
+  # RunOptions
+  IndPeriod_Run <- seq(
+    length(InputsModel[[1]]$DatesR) - 365,
+    length(InputsModel[[1]]$DatesR)
+  )
+  IndPeriod_WarmUp = seq(IndPeriod_Run[1]-366,IndPeriod_Run[1]-1)
+  RunOptions <- CreateRunOptions(
+    InputsModel = InputsModel,
+    IndPeriod_WarmUp = IndPeriod_WarmUp,
+    IndPeriod_Run = IndPeriod_Run
+  )
+  Param <- list("54057" = c(0.727,  175.493,   -0.082,    0.029,    4.654))
+  OM_GriwrmInputs <- RunModel(
+    InputsModel,
+    RunOptions = RunOptions,
+    Param = Param
+  )
+  supervisor <- CreateSupervisor(InputsModel)
+  OM_Supervisor <- RunModel(
+    supervisor,
+    RunOptions = RunOptions,
+    Param = Param
+  )
+  expect_equal(OM_Supervisor[["54057"]]$Qsim, OM_GriwrmInputs[["54057"]]$Qsim)
+})