diff --git a/R/CreateInputsModel.GRiwrm.R b/R/CreateInputsModel.GRiwrm.R
index 2d8ad3bacab51ea6d3186409929bc39a01b2166d..04ab3f8dab0185524b8849d5e8c4b92e8dc9f990 100644
--- a/R/CreateInputsModel.GRiwrm.R
+++ b/R/CreateInputsModel.GRiwrm.R
@@ -11,8 +11,8 @@
#' @export
CreateInputsModel.GRiwrm <- function(x, DatesR, Precip, PotEvap, Qobs, ...) {
- InputsModel <- CreateEmptyGRiwrmInputsModel()
- Qobs[is.na(Qobs)] <- -99 # airGRCreateInputsModel doesn't accept NA values
+ InputsModel <- CreateEmptyGRiwrmInputsModel(x)
+ Qobs[is.na(Qobs)] <- -99 # airGR::CreateInputsModel doesn't accept NA values
for(id in getNodeRanking(x)) {
message("CreateInputsModel.GRiwrm: Treating sub-basin ", id, "...")
@@ -27,9 +27,10 @@ CreateInputsModel.GRiwrm <- function(x, DatesR, Precip, PotEvap, Qobs, ...) {
#' Create an empty InputsModel object for **airGRiwrm** nodes
#'
#' @return \emph{GRiwrmInputsModel} empty object
-CreateEmptyGRiwrmInputsModel <- function() {
+CreateEmptyGRiwrmInputsModel <- function(griwrm) {
InputsModel <- list()
class(InputsModel) <- c("GRiwrmInputsModel", class(InputsModel))
+ attr(InputsModel, "GRiwrm") <- griwrm
return(InputsModel)
}
diff --git a/R/CreateSupervisor.R b/R/CreateSupervisor.R
index 07ab00bb7e0aeb07c75a5853887f0badf29dd793..0d698a50d13423efb58cb38658b48120c8bd3d5c 100644
--- a/R/CreateSupervisor.R
+++ b/R/CreateSupervisor.R
@@ -25,7 +25,7 @@
#' 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
+ # Create Supervisor environment from the parent of GlobalEnv
e <- new.env(parent = parent.env(globalenv()))
class(e) <- c("Supervisor", class(e))
@@ -35,8 +35,9 @@ CreateSupervisor <- function(InputsModel) {
# Add pointer to itself in order to assign variable from function environment
e$supervisor <- e
- # Copy of the InputsModel
+ # Copy of InputsModel, griwrm and prepare OutputsModel
e$InputsModel <- InputsModel
+ e$griwrm <- attr(InputsModel, "GRiwrm")
e$OutputsModel <- list()
# Controller list
diff --git a/R/RunModel.Supervisor.R b/R/RunModel.Supervisor.R
index 9f35149d53bc13077d5731405d97d0a61fa2e195..6a10e2bec83f977eca69adeb7ee868964540727b 100644
--- a/R/RunModel.Supervisor.R
+++ b/R/RunModel.Supervisor.R
@@ -9,9 +9,7 @@
#' @export
RunModel.Supervisor <- function(x, RunOptions, Param, ...) {
- x$ts.index0 <- sapply(RunOptions, function(x) {
- x$IndPeriod_Run[1] - 1
- })
+ x$ts.index0 <- RunOptions[[1]]$IndPeriod_Run[1] - 1
# Run runoff model for each sub-basin
x$OutputsModel <- lapply(X = x$InputsModel, FUN = function(IM) {
@@ -20,6 +18,13 @@ RunModel.Supervisor <- function(x, RunOptions, Param, ...) {
Param = Param[[IM$id]])
})
class(x$OutputsModel) <- append(class(x$OutputsModel), "GRiwrmOutputsModel")
+ # Copy simulated pure runoff flows (no SD nodes) to Qupstream in downstream SD nodes
+ for(id in getNoSD_Ids(x$InputsModel)) {
+ downId <- x$InputsModel[[id]]$down
+ x$InputsModel[[downId]]$Qupstream[RunOptions[[downId]]$IndPeriod_Run, id] <-
+ x$OutputsModel[[id]]$Qsim
+ }
+
# Save Qsim for step by step simulation
Qsim <- lapply(x$OutputsModel, function(OM) {
OM$Qsim
@@ -34,19 +39,17 @@ RunModel.Supervisor <- function(x, RunOptions, Param, ...) {
# 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.index <- iTS - x$ts.index0
x$ts.date <- x$InputsModel[[1]]$DatesR[iTS]
- doSupervision(x)
+ # 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
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)
@@ -54,9 +57,15 @@ RunModel.Supervisor <- function(x, RunOptions, Param, ...) {
x$InputsModel[[id]],
RunOptions = RunOptions[[id]],
Param = Param[[id]],
- QsimDown = Qsim[[id]][iTS - x$ts.index0]
+ QsimDown = Qsim[[id]][x$ts.index]
)
- Qsim[[id]][iTS - x$ts.index0] <- x$OutputsModel[[id]]$Qsim
+ # Storing Qsim in the data.frame Qsim
+ Qsim[[id]][x$ts.index] <- x$OutputsModel[[id]]$Qsim
+ # Routing Qsim to the downstream node
+ if(!is.na(x$InputsModel[[id]]$down)) {
+ x$InputsModel[[x$InputsModel[[id]]$down]]$Qupstream[iTS, i] <-
+ x$OutputsModel[[id]]$Qsim
+ }
}
}
for(id in getSD_Ids(x$InputsModel)) {
diff --git a/R/RunModel_Lag.R b/R/RunModel_Lag.R
index b625ac1f6b98a679d714a4224ed580ce44a3a9aa..398478a045c1fa0634ac6a4fcb1e2a0a046105a5 100644
--- a/R/RunModel_Lag.R
+++ b/R/RunModel_Lag.R
@@ -78,20 +78,17 @@ RunModel_Lag <- function(InputsModel, RunOptions, Param) {
#message("Initstates: ",paste(IniStates, collapse = ", "))
for (upstream_basin in seq_len(NbUpBasins)) {
- Qupstream <- InputsModel$Qupstream[RunOptions$IndPeriod_Run, upstream_basin]
+ Qupstream <- c(IniStates[[upstream_basin]],
+ 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)])
- #message("Qupstream1: ", paste(Qupstream1, collapse = ", "))
- #message("Qupstream2: ", paste(Qupstream2, collapse = ", "))
-
- OutputsModel$Qsim <- OutputsModel$Qsim +
- Qupstream1 * HUTRANS[1, upstream_basin] +
- Qupstream2 * HUTRANS[2, upstream_basin]
+ #message("Qupstream[", upstream_basin, "]: ", paste(Qupstream, collapse = ", "))
+ OutputsModel$Qsim <-
+ OutputsModel$Qsim +
+ Qupstream[2:(1 + LengthTs)] * HUTRANS[1, upstream_basin] +
+ Qupstream[1:LengthTs] * HUTRANS[2, upstream_basin]
}
# Warning for negative flows
if (any(OutputsModel$Qsim < 0)) {
@@ -105,12 +102,7 @@ RunModel_Lag <- function(InputsModel, RunOptions, Param) {
if ("StateEnd" %in% RunOptions$Outputs_Sim) {
OutputsModel$StateEnd$SD <- lapply(seq(NbUpBasins), function(x) {
LengthTs <- tail(RunOptions$IndPeriod_Run,1)
- Qupstream <- InputsModel$Qupstream[(LengthTs - floor(PT[x])):LengthTs, x]
- if (!is.na(InputsModel$BasinAreas[x])) {
- # Upstream flow with area needs to be converted to m3 by time step
- Qupstream <- Qupstream * InputsModel$BasinAreas[x] * 1e3
- }
- return(Qupstream)
+ InputsModel$Qupstream[(LengthTs - floor(PT[x])):LengthTs, x]
})
#message("StateEnd: ",paste(OutputsModel$StateEnd$SD, collapse = ", "))
}
diff --git a/R/utils.R b/R/utils.R
index 4cf16933935889bd01ffdb97c01e305beb789d47..5f5650640d91b4b7127bfe621cc23427f6f607d0 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -14,35 +14,53 @@ getSD_Ids <- function(InputsModel) {
names(InputsModel)[bSDs]
}
+#' Id of sub-basins not using SD model
+#'
+#' @param InputsModel `GRiwrmInputsModel` object
+#'
+#' @return [character] IDs of the sub-basins not using SD model
+#'
+getNoSD_Ids <- function(InputsModel) {
+ if (!inherits(InputsModel, "GRiwrmInputsModel")) {
+ stop("Argument `InputsModel` should be of class GRiwrmInputsModel")
+ }
+ bSDs <- sapply(InputsModel, function (IM) {
+ !inherits(IM, "SD")
+ })
+ names(InputsModel)[bSDs]
+}
+
#' Retrieve data in the model for the current time steps
#'
#' This function should be call inside a Supervisor
#'
#' @param loc location of the data
-#' @param supervisor `Supervisor` (See [CreateSupervisor])
+#' @param sv a `Supervisor` (See [CreateSupervisor])
#'
#' @return [numeric] retrieved data at the location
-getDataFromLocation <- function(loc, supervisor) {
- if (grep("\\[[0-9]+\\]$", loc)) {
+getDataFromLocation <- function(loc, sv) {
+ if (length(grep("\\[[0-9]+\\]$", loc)) > 0) {
stop("Reaching output of other controller is not implemented yet")
} else {
- supervisor$OutputsModel[[loc]]$Qsim[supervisor$ts.index - 1]
+ sv$InputsModel[[sv$griwrm$down[sv$griwrm$id == loc]]]$Qupstream[sv$ts.index0 + sv$ts.index - 1, loc]
}
}
#' Write data to model input for the current time step
#'
-#' @param control [list] A row of the `U` [data.frame] from a `Controller`
-#' @param supervisor `Supervisor` (See [CreateSupervisor])
+#' @param control [vector] A row of the `U` [data.frame] from a `Controller`
+#' @param sv `Supervisor` (See [CreateSupervisor])
#'
#' @return [NULL]
-setDataToLocation <- function(control, supervisor) {
- node <- supervisor$InputsModel[[control$loc]]$down
+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
- supervisor$InputsModel[[node]]$Qupstream[supervisor$ts.index0[node] + supervisor$ts.index, control$loc] <-
- control$v
+ sv$InputsModel[[node]]$Qupstream[sv$ts.index0 + sv$ts.index, control[1]] <-
+ as.numeric(control[2])
+ message("setDataToLocation[", control[1], "] <- ", control[2])
}
@@ -55,12 +73,12 @@ doSupervision <- function(supervisor) {
for (id in names(supervisor$controllers)) {
# Read Y from locations in the model
supervisor$controllers[[id]]$Y$v <-
- sapply(supervisor$controllers[[id]]$Y$loc, getDataFromLocation, supervisor = supervisor)
+ sapply(supervisor$controllers[[id]]$Y$loc, getDataFromLocation, sv = supervisor)
# Run logic
supervisor$controllers[[id]]$U$v <-
sapply(supervisor$controllers[[id]]$Y$v, supervisor$controllers[[id]]$FUN)
# Write U to locations in the model
- sapply(supervisor$controllers[[id]]$U, setDataToLocation, supervisor = supervisor)
+ apply(supervisor$controllers[[id]]$U, 1, setDataToLocation, sv = supervisor)
}
return()
}
diff --git a/tests/testthat/test-RunModel.R b/tests/testthat/test-RunModel.R
index e32f9311fe752831531848b490d6672ad5f30e1d..ac972f392fdd38029f8b43aac025fcfea7b1487d 100644
--- a/tests/testthat/test-RunModel.R
+++ b/tests/testthat/test-RunModel.R
@@ -1,44 +1,83 @@
context("RunModel.Supervisor")
+# Load data
+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
+nTS <- 365
+IndPeriod_Run <- seq(
+ length(InputsModel[[1]]$DatesR) - nTS + 1,
+ 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
+)
+
+# RunModel.GRiwrmInputsModel
+Param <- list("54057" = c(0.727, 175.493, -0.082, 0.029, 4.654))
+OM_GriwrmInputs <- RunModel(
+ InputsModel,
+ RunOptions = RunOptions,
+ Param = Param
+)
+
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
- nTS <- 365
- IndPeriod_Run <- seq(
- length(InputsModel[[1]]$DatesR) - nTS + 1,
- 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,
+ sv <- CreateSupervisor(InputsModel)
+ OM_Supervisor <- RunModel(
+ sv,
RunOptions = RunOptions,
Param = Param
)
- supervisor <- CreateSupervisor(InputsModel)
+ expect_equal(OM_Supervisor[["54057"]]$Qsim, OM_GriwrmInputs[["54057"]]$Qsim)
+})
+
+test_that("RunModelSupervisor with two regulations that cancel each other out should returns same results as RunModel.GRiwrmInputsModel", {
+ # Add 2 nodes to the network
+ griwrm2 <- rbind(griwrm,
+ data.frame(
+ id = c("R1", "R2"),
+ down = "54057",
+ length = 100000,
+ area = NA,
+ model = NA
+ ))
+ # Add Qobs for the 2 new nodes
+ Qobs2 <- cbind(Qobs, matrix(data = rep(0, 2*nrow(Qobs)), ncol = 2))
+ colnames(Qobs2) <- c(colnames(Qobs2)[1:6], "R1", "R2")
+ 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 }
+ # Function to release half of the the measured flow
+ fRelease <- function(y) { message("fWithdrawal(", 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"
+ createController(sv, "Release", Y = c("54002"), U = c("R2"), FUN = fRelease)
OM_Supervisor <- RunModel(
- supervisor,
+ sv,
RunOptions = RunOptions,
Param = Param
)
expect_equal(OM_Supervisor[["54057"]]$Qsim, OM_GriwrmInputs[["54057"]]$Qsim)
})
+