From 2d16783be45c2d96317609e749f8595a0ed0035a Mon Sep 17 00:00:00 2001
From: David <david.dorchies@inrae.fr>
Date: Wed, 26 Oct 2022 18:49:48 +0200
Subject: [PATCH] feat(RunModel.Supervision): Diversion node implementation
- Modify code to run Diversion nodes with Supervision
- Add vignette
Refs #95, #100
---
R/CreateGRiwrm.R | 8 +-
R/CreateInputsModel.GRiwrm.R | 3 +-
R/RunModel.InputsModel.R | 19 +-
R/RunModel.SD.R | 6 +-
R/RunModel.Supervisor.R | 89 +++++---
R/plot.GRiwrm.R | 6 +-
R/utils.R | 31 ++-
man/getNoSD_Ids.Rd | 4 +-
man/getSD_Ids.Rd | 4 +-
.../V06_Modelling_regulated_diversion.Rmd | 205 ++++++++++++++++++
10 files changed, 320 insertions(+), 55 deletions(-)
create mode 100644 vignettes/V06_Modelling_regulated_diversion.Rmd
diff --git a/R/CreateGRiwrm.R b/R/CreateGRiwrm.R
index 81c9c59..2c46b0e 100644
--- a/R/CreateGRiwrm.R
+++ b/R/CreateGRiwrm.R
@@ -256,15 +256,15 @@ getNodeProperties <- function(id, griwrm) {
upstreamIds <- griwrm$id[!griwrm$id %in% griwrm$down]
gaugedIds <- griwrm$id[!is.na(griwrm$model) & griwrm$model != "Ungauged"]
divertedIds <- griwrm$id[!is.na(griwrm$model) & griwrm$model == "Diversion"]
- p <- c(
+ p <- list(
position = ifelse(id %in% upstreamIds, "Upstream", "Intermediate"),
hydrology = ifelse(id %in% gaugedIds, "Gauged",
ifelse(is.na(griwrm$model[griwrm$id == id]),
"DirectInjection",
"Ungauged"))
)
- if (length(divertedIds) > 0) {
- if (id %in% divertedIds) p["diverted"] <- "Diversion"
- }
+ p$Upstream <- p$position == "Upstream"
+ p$DirectInjection = p$hydrology == "DirectInjection"
+ p$Diversion = id %in% divertedIds
return(p)
}
diff --git a/R/CreateInputsModel.GRiwrm.R b/R/CreateInputsModel.GRiwrm.R
index ec38978..6eb5dfb 100644
--- a/R/CreateInputsModel.GRiwrm.R
+++ b/R/CreateInputsModel.GRiwrm.R
@@ -245,10 +245,9 @@ CreateEmptyGRiwrmInputsModel <- function(griwrm) {
#' @return \emph{InputsModel} object for one.
#' @noRd
CreateOneGRiwrmInputsModel <- function(id, griwrm, ..., Qobs, Qmin) {
- hasDiversion <- "Diversion" %in% getNodeProperties(id, griwrm)
+ hasDiversion <- getNodeProperties(id, griwrm)$Diversion
if (hasDiversion) {
rowDiv <- which(griwrm$id == id & griwrm$model == "Diversion")
- hasDiversion <- TRUE
diversionOutlet <- griwrm$down[rowDiv]
griwrm <- griwrm[-rowDiv, ]
}
diff --git a/R/RunModel.InputsModel.R b/R/RunModel.InputsModel.R
index fd38a4a..12727eb 100644
--- a/R/RunModel.InputsModel.R
+++ b/R/RunModel.InputsModel.R
@@ -37,21 +37,30 @@ RunModel.InputsModel <- function(x, RunOptions, Param, FUN_MOD = NULL, ...) {
#' On a Diversion node, this function is called after `airGR::RunModel` to
#' divert a part of the flow to another node than the original downstream one.
#'
-#' @param InputsModel \[object of class \emph{InputsModel}\] see [airGR::CreateInputsModel] for details
+#' @param InputsModel \[object of class \emph{InputsModel}\] see
+#' [airGR::CreateInputsModel] for details
+#' @param RunOptions Same parameter as in [RunModel.GRiwrmInputsModel]
#' @param OutputsModel Output of [airGR::RunModel]
+#' @param updateQsim [logical] for updating Qsim after diversion in the output
#'
-#' @return Updated OutputsModel object with diversion:
-#' Qsim
+#' @return Updated `OutputsModel` object after diversion
#' @noRd
#'
-RunModel_Diversion <- function(InputsModel, RunOptions, OutputsModel) {
+RunModel_Diversion <- function(InputsModel,
+ RunOptions,
+ OutputsModel,
+ updateQsim = TRUE) {
OutputsModel$Qnat <- OutputsModel$Qsim
lQ <- calc_Qdiv(OutputsModel$Qsim_m3,
InputsModel$Qdiv[RunOptions$IndPeriod_Run],
InputsModel$Qmin[RunOptions$IndPeriod_Run])
+ #message(paste(InputsModel$Qdiv[RunOptions$IndPeriod_Run], lQ$Qdiv, lQ$Qsim, InputsModel$Qmin[RunOptions$IndPeriod_Run], sep = ", "))
OutputsModel$Qdiv_m3 <- lQ$Qdiv
OutputsModel$Qsim_m3 <- lQ$Qsim
- OutputsModel$Qsim <- OutputsModel$Qsim_m3 / sum(InputsModel$BasinAreas, na.rm = TRUE) / 1e3
+ if (updateQsim) {
+ OutputsModel$Qsim <-
+ OutputsModel$Qsim_m3 / sum(InputsModel$BasinAreas, na.rm = TRUE) / 1e3
+ }
if ("WarmUpQsim" %in% RunOptions$Outputs_Sim) {
lQ <- calc_Qdiv(OutputsModel$RunOptions$WarmUpQsim_m3,
InputsModel$Qdiv[RunOptions$IndPeriod_WarmUp],
diff --git a/R/RunModel.SD.R b/R/RunModel.SD.R
index 6d141d4..5498bcf 100644
--- a/R/RunModel.SD.R
+++ b/R/RunModel.SD.R
@@ -8,5 +8,9 @@
#' @export
#'
RunModel.SD <- function(x, RunOptions, Param, QcontribDown, ...) {
- airGR::RunModel_Lag(x, RunOptions = RunOptions, Param = Param[1], QcontribDown = QcontribDown)
+ OutputsModel <- airGR::RunModel_Lag(x,
+ RunOptions = RunOptions,
+ Param = Param[1],
+ QcontribDown = QcontribDown)
+ return(OutputsModel)
}
diff --git a/R/RunModel.Supervisor.R b/R/RunModel.Supervisor.R
index ed7a588..78f6f2b 100644
--- a/R/RunModel.Supervisor.R
+++ b/R/RunModel.Supervisor.R
@@ -33,32 +33,35 @@ RunModel.Supervisor <- function(x, RunOptions, Param, ...) {
# Run runoff model for each sub-basin
x$OutputsModel <- lapply(X = x$InputsModel, FUN = function(IM) {
- RunModel.GR(IM,
- RunOptions = RunOptions[[IM$id]],
- Param = Param[[IM$id]])
+ OM_GR <- RunModel.GR(IM,
+ RunOptions = RunOptions[[IM$id]],
+ Param = Param[[IM$id]])
+ if (IM$hasDiversion) OM_GR$Qnat <- OM_GR$Qsim
+ return(OM_GR)
})
- 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
- if(!is.null(x$InputsModel[[downId]])) {
- x$InputsModel[[downId]]$Qupstream[RunOptions[[downId]]$IndPeriod_Run, id] <-
- x$OutputsModel[[id]]$Qsim_m3
- }
+ class(x$OutputsModel) <- c("GRiwrmOutputsModel", class(x$OutputsModel))
+
+ # Copy simulated pure runoff flows (no SD nor Diversion nodes) to Qupstream
+ for(id in getNoSD_Ids(x$InputsModel, include_diversion = FALSE)) {
+ updateQupstream.Supervisor(x, id, IndPeriod_Run)
}
- # Save Qsim for step by step simulation
+ # Save OutputsModel for step by step simulation
QcontribDown <- do.call(
cbind,
lapply(x$OutputsModel, "[[", "Qsim")
)
-
Qsim_m3 <- do.call(
cbind,
lapply(x$OutputsModel, "[[", "Qsim_m3")
)
+ if (length(getDiversionRows(x$griwrm)) > 0) {
+ # Outputs of Diversion nodes
+ Qdiv_m3 <- Qsim_m3[, sv$griwrm$id[getDiversionRows(x$griwrm)], drop = FALSE] * NA
+ Qnat <- Qdiv_m3
+ }
- # Initialisation of model states by running the model with no supervision on warm-up period
+ # Initialization of model states by running the model with no supervision on warm-up period
RunOptionsWarmUp <- RunOptions
for(id in names(x$InputsModel)) {
RunOptionsWarmUp[[id]]$IndPeriod_Run <- RunOptionsWarmUp[[id]]$IndPeriod_WarmUp
@@ -96,23 +99,35 @@ RunModel.Supervisor <- function(x, RunOptions, Param, ...) {
doSupervision(x)
}
# Loop over sub-basin using SD model
- for(id in getSD_Ids(x$InputsModel)) {
- # Run the SD model for the sub-basin and one time step
- RunOptions[[id]]$IndPeriod_Run <- iTS
+ for(id in getSD_Ids(x$InputsModel, add_diversion = TRUE)) {
+ # Run model for the sub-basin and one time step
RunOptions[[id]]$IniStates <- serializeIniStates(x$OutputsModel[[id]]$StateEnd)
- x$OutputsModel[[id]] <- RunModel.SD(
- x$InputsModel[[id]],
- RunOptions = RunOptions[[id]],
- Param = Param[[id]],
- QcontribDown = QcontribDown[x$ts.index, id]
- )
- # Storing Qsim in the data.frame Qsim
+ RunOptions[[id]]$IndPeriod_Run <- iTS
+ if (RunOptions[[id]]$FeatFUN_MOD$IsSD) {
+ # Route upstream flows for SD nodes
+ x$OutputsModel[[id]] <- RunModel.SD(
+ x$InputsModel[[id]],
+ RunOptions = RunOptions[[id]],
+ Param = Param[[id]],
+ QcontribDown = QcontribDown[x$ts.index, id]
+ )
+ } else {
+ x$OutputsModel[[id]]$Qsim_m3 <- Qsim_m3[x$ts.index, id]
+ }
+ if (x$InputsModel[[id]]$hasDiversion) {
+ # Compute diverted and simulated flows on Diversion nodes
+ x$OutputsModel[[id]] <-
+ RunModel_Diversion(x$InputsModel[[id]],
+ RunOptions = RunOptions[[id]],
+ OutputsModel = x$OutputsModel[[id]])
+ }
+ # Storing Qsim_m3 and Qdiv_m3 data.frames
Qsim_m3[x$ts.index, id] <- x$OutputsModel[[id]]$Qsim_m3
- # Routing Qsim to Qupstream of downstream nodes
- if(!is.na(x$InputsModel[[id]]$down)) {
- x$InputsModel[[x$InputsModel[[id]]$down]]$Qupstream[iTS, id] <-
- x$OutputsModel[[id]]$Qsim_m3
+ if (x$InputsModel[[id]]$hasDiversion) {
+ Qdiv_m3[x$ts.index, id] <- x$OutputsModel[[id]]$Qdiv_m3
}
+ # Routing Qsim_m3 and Qdiv_m3 to Qupstream of downstream nodes
+ updateQupstream.Supervisor(x, id, iTS)
}
x$ts.previous <- x$ts.index
}
@@ -123,6 +138,9 @@ RunModel.Supervisor <- function(x, RunOptions, Param, ...) {
x$OutputsModel[[id]]$Qsim_m3 <- Qsim_m3[, id]
x$OutputsModel[[id]]$Qsim <-
Qsim_m3[, id] / sum(x$InputsModel[[id]]$BasinAreas, na.rm = TRUE) / 1e3
+ if (x$InputsModel[[id]]$hasDiversion) {
+ x$OutputsModel[[id]]$Qdiv_m3 <- Qdiv_m3[, id]
+ }
}
attr(x$OutputsModel, "Qm3s") <- OutputsModelQsim(x$InputsModel, x$OutputsModel, IndPeriod_Run)
@@ -131,3 +149,18 @@ RunModel.Supervisor <- function(x, RunOptions, Param, ...) {
return(x$OutputsModel)
}
+
+updateQupstream.Supervisor <- function(x, id, iTS) {
+ downId <- x$InputsModel[[id]]$down
+ if(!is.null(x$InputsModel[[downId]])) {
+ x$InputsModel[[downId]]$Qupstream[iTS, id] <-
+ x$OutputsModel[[id]]$Qsim_m3
+ }
+ if (x$InputsModel[[id]]$hasDiversion) {
+ divOutId <- x$InputsModel[[id]]$diversionOutlet
+ if (!is.na(divOutId)) {
+ x$InputsModel[[divOutId]]$Qupstream[iTS, id] <-
+ x$OutputsModel[[id]]$Qdiv_m3
+ }
+ }
+}
diff --git a/R/plot.GRiwrm.R b/R/plot.GRiwrm.R
index 88b49c9..aa2ee08 100644
--- a/R/plot.GRiwrm.R
+++ b/R/plot.GRiwrm.R
@@ -82,11 +82,11 @@ plot.GRiwrm <- function(x,
getNodeClass <- function(id, griwrm) {
props <- getNodeProperties(id, griwrm)
- if (props["hydrology"] == "DirectInjection") {
- nc <- props["hydrology"]
+ if (props$DirectInjection) {
+ nc <- "DirectInjection"
} else {
nc <- paste0(props["position"], props["hydrology"])
}
- if("diverted" %in% names(props)) nc <- paste0(nc, "Diversion")
+ if (props$Diversion) nc <- paste0(nc, "Diversion")
return(nc)
}
diff --git a/R/utils.R b/R/utils.R
index 286ee69..4e961f9 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -1,15 +1,16 @@
#' Function to obtain the ID of sub-basins using SD model
#'
#' @param InputsModel \[`GRiwrmInputsModel` object\]
+#' @param add_diversion [logical] for adding upstream nodes with diversion
#'
#' @return [character] IDs of the sub-basins using SD model
#' @export
-getSD_Ids <- function(InputsModel) {
+getSD_Ids <- function(InputsModel, add_diversion = FALSE) {
if (!inherits(InputsModel, "GRiwrmInputsModel")) {
stop("Argument `InputsModel` should be of class GRiwrmInputsModel")
}
bSDs <- sapply(InputsModel, function (IM) {
- inherits(IM, "SD")
+ inherits(IM, "SD") | IM$hasDiversion
})
names(InputsModel)[bSDs]
}
@@ -17,15 +18,16 @@ getSD_Ids <- function(InputsModel) {
#' Function to obtain the ID of sub-basins not using SD model
#'
#' @param InputsModel \[`GRiwrmInputsModel` object\]
+#' @param include_diversion [logical] for including diversion nodes
#'
#' @return [character] IDs of the sub-basins not using the SD model
#' @export
-getNoSD_Ids <- function(InputsModel) {
+getNoSD_Ids <- function(InputsModel, include_diversion = TRUE) {
if (!inherits(InputsModel, "GRiwrmInputsModel")) {
stop("Argument `InputsModel` should be of class GRiwrmInputsModel")
}
bSDs <- sapply(InputsModel, function (IM) {
- !inherits(IM, "SD")
+ !inherits(IM, "SD") & (include_diversion | !IM$hasDiversion)
})
names(InputsModel)[bSDs]
}
@@ -45,7 +47,7 @@ getDataFromLocation <- function(loc, sv) {
stop("Reaching output of other controller is not implemented yet")
} else {
if(sv$nodeProperties[[loc]]["hydrology"] != "DirectInjection") {
- if (sv$nodeProperties[[loc]]["position"] == "Upstream") {
+ if (sv$nodeProperties[[loc]]$Upstream) {
sv$OutputsModel[[loc]]$Qsim_m3[sv$ts.previous]
} else {
sv$OutputsModel[[loc]]$Qsim_m3
@@ -67,13 +69,22 @@ getDataFromLocation <- function(loc, sv) {
#' @noRd
setDataToLocation <- function(ctrlr, sv) {
l <- lapply(seq(length(ctrlr$Unames)), function(i) {
- node <- sv$griwrm4U$down[sv$griwrm4U$id == ctrlr$Unames[i]]
# limit U size to the number of simulation time steps of the current supervision time step
U <- ctrlr$U[seq.int(length(sv$ts.index)),i]
- # ! Qupstream contains warm up period and run period => the index is shifted
- if(!is.null(sv$InputsModel[[node]])) {
- sv$InputsModel[[node]]$Qupstream[sv$ts.index0 + sv$ts.index,
- ctrlr$Unames[i]] <- U
+
+ locU <- ctrlr$Unames[i]
+ if (sv$nodeProperties[[locU]]$DirectInjection) {
+ # Direct injection node => update Qusptream of downstream node
+ node <- sv$griwrm4U$down[sv$griwrm4U$id == locU]
+ # ! Qupstream contains warm up period and run period => the index is shifted
+ if(!is.null(sv$InputsModel[[node]])) {
+ sv$InputsModel[[node]]$Qupstream[sv$ts.index0 + sv$ts.index, locU] <- U
+ }
+ } else if (sv$nodeProperties[[locU]]$Diversion){
+ # Diversion node => update Qdiv with -U
+ sv$InputsModel[[locU]]$Qdiv[sv$ts.index0 + sv$ts.index] <- -U
+ } else {
+ stop("Node ", locU, " must be a Direct Injection or a Diversion node")
}
})
}
diff --git a/man/getNoSD_Ids.Rd b/man/getNoSD_Ids.Rd
index 5075069..644d49e 100644
--- a/man/getNoSD_Ids.Rd
+++ b/man/getNoSD_Ids.Rd
@@ -4,10 +4,12 @@
\alias{getNoSD_Ids}
\title{Function to obtain the ID of sub-basins not using SD model}
\usage{
-getNoSD_Ids(InputsModel)
+getNoSD_Ids(InputsModel, include_diversion = TRUE)
}
\arguments{
\item{InputsModel}{[\code{GRiwrmInputsModel} object]}
+
+\item{include_diversion}{\link{logical} for including diversion nodes}
}
\value{
\link{character} IDs of the sub-basins not using the SD model
diff --git a/man/getSD_Ids.Rd b/man/getSD_Ids.Rd
index 1a86310..a1987ec 100644
--- a/man/getSD_Ids.Rd
+++ b/man/getSD_Ids.Rd
@@ -4,10 +4,12 @@
\alias{getSD_Ids}
\title{Function to obtain the ID of sub-basins using SD model}
\usage{
-getSD_Ids(InputsModel)
+getSD_Ids(InputsModel, add_diversion = FALSE)
}
\arguments{
\item{InputsModel}{[\code{GRiwrmInputsModel} object]}
+
+\item{add_diversion}{\link{logical} for adding upstream nodes with diversion}
}
\value{
\link{character} IDs of the sub-basins using SD model
diff --git a/vignettes/V06_Modelling_regulated_diversion.Rmd b/vignettes/V06_Modelling_regulated_diversion.Rmd
new file mode 100644
index 0000000..f8f4c36
--- /dev/null
+++ b/vignettes/V06_Modelling_regulated_diversion.Rmd
@@ -0,0 +1,205 @@
+---
+title: "Severn_06: Modelling of a regulated diversion"
+output: rmarkdown::html_vignette
+vignette: >
+ %\VignetteIndexEntry{Severn_06: Modelling of a regulated diversion}
+ %\VignetteEngine{knitr::rmarkdown}
+ %\VignetteEncoding{UTF-8}
+---
+
+```{r, include = FALSE}
+knitr::opts_chunk$set(
+ collapse = TRUE,
+ comment = "#>",
+ fig.width = 6,
+ fig.asp = 0.68,
+ out.width = "70%",
+ fig.align = "center"
+)
+```
+
+```{r setup}
+library(airGRiwrm)
+data(Severn)
+```
+
+## The study case
+
+In this vignette, we are still using the study case of the vignette #1 and #2, and
+we are going to simulate a diversion at the node "54029" to provide flow to the node
+"54001" in order to support low flows in the Severn river.
+
+This objective is to maintain a minimum flow equals to the 3th decile of the flow
+at station "54001" without remaining less than the first decile of the flow downstream
+the station "54029".
+
+Let's compute the flow quantiles at each station in m3/s:
+
+```{r}
+Qobs <- cbind(sapply(Severn$BasinsObs, function(x) {x$discharge_spec}))
+Qobs <- Qobs[, Severn$BasinsInfo$gauge_id]
+Qobs_m3s <- t(apply(Qobs, 1, function(r) r * Severn$BasinsInfo$area * 1E3 / 86400))
+apply(Qobs_m3s[, c("54029", "54001")], 2, quantile, probs = seq(0,1,0.1), na.rm = TRUE)
+```
+
+The rule to apply is expressed as follow:
+
+> The diversion from "54029" is computed to maintain a minimum flow of 20 m3/s at
+the gauging station "54001". The diversion is allowed as far as the flow at the
+gauging station "54029" is above 2.6 m3/s.
+
+## Network
+
+```{r}
+
+nodes_div <- Severn$BasinsInfo[, c("gauge_id", "downstream_id", "distance_downstream", "area")]
+nodes_div$model <- "RunModel_GR4J"
+nodes_div <- rbind(nodes_div, data.frame(gauge_id = "54029",
+ downstream_id = "54001",
+ distance_downstream = 25,
+ model = "Diversion",
+ area = NA))
+renameCols <- list(id = "gauge_id", down = "downstream_id", length = "distance_downstream")
+griwrmV06 <- CreateGRiwrm(nodes_div, renameCols)
+plot(griwrmV06)
+```
+
+## GRiwrmInputsModel object
+
+We produce below the same operations as in the vignette "V02_Calibration_SD_model" to prepare the input data:
+
+```{r}
+data(Severn)
+nodes <- Severn$BasinsInfo[, c("gauge_id", "downstream_id", "distance_downstream", "area")]
+nodes$model <- "RunModel_GR4J"
+griwrm <- CreateGRiwrm(nodes, list(id = "gauge_id", down = "downstream_id", length = "distance_downstream"))
+BasinsObs <- Severn$BasinsObs
+DatesR <- BasinsObs[[1]]$DatesR
+PrecipTot <- cbind(sapply(BasinsObs, function(x) {x$precipitation}))
+PotEvapTot <- cbind(sapply(BasinsObs, function(x) {x$peti}))
+Precip <- ConvertMeteoSD(griwrm, PrecipTot)
+PotEvap <- ConvertMeteoSD(griwrm, PotEvapTot)
+```
+
+The main difference here is that we need to provide a diverted flow and a minimum
+remaining flow at station "54029". As, the diverted flow will be calculated during
+simulation, we provide a initial diverted flow equals to zero for all the time steps.
+
+```{r}
+Qdiv <- matrix(rep(0, length(DatesR)), ncol = 1)
+colnames(Qdiv) <- "54029"
+Qmin <- matrix(rep(2.6 * 86400, length(DatesR)), ncol = 1)
+colnames(Qmin) <- "54029"
+IM_div <- CreateInputsModel(griwrmV06, DatesR, Precip, PotEvap, Qobs = Qdiv, Qmin = Qmin)
+```
+
+## Implementation of the regulation controller
+
+### The supervisor
+
+The simulation is piloted through a `Supervisor` that can contain one or more
+`Controller`. The supervision time step will be here the same as the simulation
+time step: 1 day. Each day, the decision is taken for the current day from the
+measurement simulated the previous time step.
+
+```{r}
+sv <- CreateSupervisor(IM_div, TimeStep = 1L)
+```
+
+### The control logic function
+
+On a Diversion node, the minimum remaining flow is provided with the inputs and
+is automatically taken into account by the model. So the control logic function
+has only to take care about how much water to divert to the gauging station
+"54002".
+
+We need to enclose the logic function in a function factory providing the
+Supervisor in the environment of the function:
+
+```{r}
+#' @param sv the Supervisor environment
+logicFunFactory <- function(sv) {
+ #' @param Y Flow measured at "54002" the previous time step
+ function(Y) {
+ Qnat <- Y
+ # We need to remove the diverted flow to compute the natural flow at "54002"
+ lastU <- sv$controllers[[sv$controller.id]]$U
+ if (length(lastU) > 0) {
+ Qnat <- Y - lastU
+ }
+ return(-max(15 * 86400 - Qnat, 0))
+ }
+}
+```
+
+### The controller
+
+We declare the controller which defines where is the measurement `Y` , where to
+apply the decision `U` with which logic function:
+
+```{r}
+CreateController(sv,
+ ctrl.id = "Low flow support",
+ Y = "54001",
+ U = "54029",
+ FUN = logicFunFactory(sv))
+```
+
+
+## Running the simulation
+
+First we need to create a `GRiwrmRunOptions` object and load the parameters calibrated in the vignette "V02_Calibration_SD_model":
+
+```{r}
+# Running simulation between 2002 and 2005
+IndPeriod_Run <- seq(
+ which(DatesR == as.POSIXct("2002-01-01")),
+ which(DatesR == as.POSIXct("2005-01-01"))
+)
+IndPeriod_WarmUp = seq(IndPeriod_Run[1] - 366,IndPeriod_Run[1] - 1)
+RunOptions <- CreateRunOptions(IM_div,
+ IndPeriod_WarmUp = IndPeriod_WarmUp,
+ IndPeriod_Run = IndPeriod_Run)
+ParamV02 <- readRDS(system.file("vignettes", "ParamV02.RDS", package = "airGRiwrm"))
+```
+
+And we run the supervised model:
+
+```{r}
+OM_div <- RunModel(sv, RunOptions = RunOptions, Param = ParamV02)
+```
+
+To compare results with diversion and without diversion, we run the model without
+the supervision (remember that we have set the diverted flow at zero in the inputs):
+
+```{r}
+OM_nat <- RunModel(IM_div, RunOptions = RunOptions, Param = ParamV02)
+```
+
+
+## Exploring results
+
+Let's plot the diverted flow, and compare the flow at stations 54029 and 54001,
+with and without low-flow support at station 54001:
+
+```{r, fig.asp=1}
+dfQdiv <- data.frame(DatesR = OM_div[[1]]$DatesR,
+ Diverted_flow = OM_div$`54029`$Qdiv_m3 / 86400)
+
+oldpar <- par(mfrow=c(3,1), mar = c(2.5,4,1,1))
+plot.Qm3s(dfQdiv)
+
+# Plot natural and influenced flow at station "54001" and "54029"
+thresholds <- c("54029" = 2.6, "54001" = 15)
+lapply(c("54029", "54001"), function(id) {
+ df <- cbind(attr(OM_div, "Qm3s")[, c("DatesR", id)],
+ attr(OM_nat, "Qm3s")[, id])
+ names(df) <- c("DatesR",
+ paste(id, "with irrigation"),
+ paste(id, "natural flow"))
+ plot.Qm3s(df, ylim = c(0,50))
+ abline(h = thresholds[id], col = "green")
+})
+par(oldpar)
+```
+
--
GitLab