Resolve "Feature request: feedback control"

Closes #19 (closed)

R/CreateInputsModel.GRiwrm.R

@@ -20,6 +20,7 @@ CreateInputsModel.GRiwrm <- function(x, DatesR, Precip, PotEvap, Qobs, ...) {
       id, x, DatesR,Precip[,id], PotEvap[,id], Qobs, ...
     )
   }
+  attr(InputsModel, "TimeStep") <- getModelTimeStep(InputsModel)
   return(InputsModel)
 }
 
@@ -92,3 +93,28 @@ CreateOneGRiwrmInputsModel <- function(id, griwrm, DatesR, Precip, PotEvap, Qobs
 
   return(InputsModel)
 }
+
+
+#' Check time steps of the model of all the nodes and return the time step in seconds
+#'
+#' This function is called inside [CreateInputsModel.GRiwrm] for defining the time step of the big model.
+#'
+#' @param InputsModel a `GRiwrmInputsModel`
+#'
+#' @return A [numeric] representing the time step in seconds
+#'
+getModelTimeStep <- function(InputsModel) {
+  TS <- sapply(InputsModel, function(x) {
+    if (inherits(x, "hourly")) {
+      TimeStep <- 60 * 60
+    } else if (inherits(x, "daily")) {
+      TimeStep <- 60 * 60 * 24
+    } else {
+      stop("All models should be at hourly or daily time step")
+    }
+  })
+  if(length(unique(TS)) != 1) {
+    stop("Time steps of the model of all nodes should be identical")
+  }
+  return(unique(TS))
+}

R/RunModel.GRiwrmInputsModel.R

@@ -14,21 +14,21 @@ RunModel.GRiwrmInputsModel <- function(x, RunOptions, Param, ...) {
   OutputsModel <- list()
   class(OutputsModel) <- c("GRiwrmOutputsModel", class(OutputsModel))
 
-  for(IM in x) {
-    message("RunModel.GRiwrmInputsModel: Treating sub-basin ", IM$id, "...")
+  for(id in names(x)) {
+    message("RunModel.GRiwrmInputsModel: Treating sub-basin ", x[[id]]$id, "...")
 
-    # Update x$Qupstream with simulated upstream flows
-    if(any(IM$UpstreamIsRunoff)) {
-      IM <- UpdateQsimUpstream(IM, RunOptions[[IM$id]]$IndPeriod_Run, OutputsModel)
+    # Update x[[id]]$Qupstream with simulated upstream flows
+    if(any(x[[id]]$UpstreamIsRunoff)) {
+      x[[id]] <- UpdateQsimUpstream(x[[id]], RunOptions[[id]]$IndPeriod_Run, OutputsModel)
     }
 
     # Run the model for the sub-basin
-    OutputsModel[[IM$id]] <- RunModel.InputsModel(
-      IM,
-      RunOptions = RunOptions[[IM$id]],
-      Param = Param[[IM$id]]
+    OutputsModel[[id]] <- RunModel.InputsModel(
+      x[[id]],
+      RunOptions = RunOptions[[id]],
+      Param = Param[[id]]
     )
-
   }
+  attr(OutputsModel, "Qm3s") <- OutputsModelQsim(x, OutputsModel, RunOptions[[1]]$IndPeriod_Run)
   return(OutputsModel)
 }

R/RunModel.Supervisor.R

@@ -44,7 +44,6 @@ RunModel.Supervisor <- function(x, RunOptions, Param, ...) {
     # Regulation occurs from second time step
     if(iTS > RunOptions[[1]]$IndPeriod_Run[1]) {
       doSupervision(x)
-      message("Supervision done")
     }
 
     # Loop over sub-basin using SD model
@@ -71,5 +70,6 @@ RunModel.Supervisor <- function(x, RunOptions, Param, ...) {
   for(id in getSD_Ids(x$InputsModel)) {
     x$OutputsModel[[id]]$Qsim <- Qsim[[id]]
   }
+  attr(x$OutputsModel, "Qm3s") <- OutputsModelQsim(x$InputsModel, x$OutputsModel)
   return(x$OutputsModel)
 }

R/utils.R

@@ -55,12 +55,10 @@ getDataFromLocation <- function(loc, sv) {
 #'
 #' @return [NULL]
 setDataToLocation <- function(control, sv) {
-  message("setDataToLocation[", control[1], "] <- ", control[2])
   node <- sv$griwrm$down[sv$griwrm$id == control[1]]
   # ! Qupstream contains warm up period and run period => the index is shifted
   sv$InputsModel[[node]]$Qupstream[sv$ts.index0 + sv$ts.index, control[1]] <-
     as.numeric(control[2])
-  message("setDataToLocation[", control[1], "] <- ", control[2])
 }
 
 
@@ -99,3 +97,55 @@ checkRunModelParameters <- function(InputsModel, RunOptions, Param) {
   if(!is.list(Param) || !all(names(InputsModel) %in% names(Param))) stop("Argument `Param` must be a list with names equal to nodes IDs")
   return()
 }
+
+
+#' Create a data.frame with simulated flows at each nodes of the [GRiwrm] object
+#'
+#' @details
+#' This function can only be called inside [RunModel.GRiwrmInputsmodel] or [RunModel.Supervisor]
+#' because it needs a `GRiwrmInputsModel` object internally modified by these functions
+#' (`Qupstream` updated with simulated flows).
+#'
+#' @param InputsModel a `GRiwrmInputsModel` object created by [CreateInputsModel.GRiwrm]
+#' @param OutputsModel a `GRiwrmOutputsModel` object created by [RunModel.GRiwrmInputsmodel] or [RunModel.Supervisor]
+#' @param IndPeriod_Run an [integer] vector (See [airGR::CreateRunOptions])
+#'
+#' @return a [data.frame] containing the simulated flows (in m3/time step) structured with the following columns:
+#' - 'DatesR' containing the timestamps of the time series
+#' - one column by node with the simulated flows
+#'
+OutputsModelQsim <- function(InputsModel, OutputsModel, IndPeriod_Run) {
+  griwrm <- attr(InputsModel, "GRiwrm")
+  # Get simulated flow for each node
+  # Flow for each node is available in InputsModel$Qupstream except for the downstream node
+  upperNodes <- griwrm$id[!is.na(griwrm$down)]
+  lQsim <- lapply(
+    upperNodes,
+    function(x, griwrm, IndPeriod_Run) {
+      node <- griwrm$down[griwrm$id == x]
+      InputsModel[[node]]$Qupstream[IndPeriod_Run, x]
+    },
+    griwrm = griwrm, IndPeriod_Run = IndPeriod_Run
+  )
+  names(lQsim) <- upperNodes
+  # Flow of the downstream node is only available in OutputsModel[[node]]$Qsim
+  downNode <- names(InputsModel)[length(InputsModel)]
+  lQsim[[downNode]] <- OutputsModel[[downNode]]$Qsim
+
+  # Conversion to m3/s
+  lQsim <- lapply(
+    names(lQsim),
+    function(x) {
+      i <- which(griwrm$id == x)
+      if(is.na(griwrm$area[i])) { # m3/time step => m3/s
+        return(lQsim[[x]] / attr(InputsModel, "TimeStep"))
+      } else { # mm/time step => m3/s
+        return(lQsim[[x]] * griwrm$area[i] * 1E3 / attr(InputsModel, "TimeStep"))
+      }
+    }
+  )
+  names(lQsim) <- c(upperNodes, downNode)
+  dfQsim <- cbind(data.frame(DatesR = as.POSIXct(InputsModel[[1]]$DatesR[IndPeriod_Run])),
+                  do.call(cbind,lQsim))
+  return(dfQsim)
+}

tests/testthat/test-RunModel.R

@@ -66,9 +66,9 @@ test_that("RunModelSupervisor with two regulations that cancel each other out sh
   InputsModel <- CreateInputsModel(griwrm2, DatesR, Precip, PotEvap, Qobs2)
   sv <- CreateSupervisor(InputsModel)
   # Function to withdraw half of the measured flow
-  fWithdrawal <- function(y) { message("fWithdrawal(", y, ")"); -y/2 }
+  fWithdrawal <- function(y) { -y/2 }
   # Function to release half of the the measured flow
-  fRelease <- function(y) { message("fWithdrawal(", y, ")"); y/2 }
+  fRelease <- function(y) { y/2 }
   # Controller that withdraw half of the flow measured at node "54002" at location "R1"
   createController(sv, "Withdrawal", Y = c("54002"), U = c("R1"), FUN = fWithdrawal)
   # Controller that release half of the flow measured at node "54002" at location "R2"