Interfacing HydroPso with GR2M

[Josh-kara]

BasinObs<-read.csv("Monthly Katar.csv",header=TRUE) InputsModel <- CreateInputsModel(FUN_MOD = RunModel_GR2M, DatesR = as.POSIXlt(BasinObs$DatesR), Precip = BasinObs$P, PotEvap = BasinObs$E) Ind_Run <- seq(which(format(as.POSIXlt(BasinObs$DatesR), format = "%m/%Y")=="01/1986"), which(format(as.POSIXlt(BasinObs$DatesR), format = "%m/%Y")=="11/1989"))

RunOptions <- CreateRunOptions(FUN_MOD = RunModel_GR2M, InputsModel = InputsModel, IndPeriod_Run = Ind_Run)

InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE, InputsModel = InputsModel, RunOptions = RunOptions, Obs = BasinObs$Qmm[Ind_Run])

CalibOptions <- CreateCalibOptions(FUN_MOD = RunModel_GR2M, FUN_CALIB = Calibration_Michel)

OutputsCalib <- Calibration_Michel(InputsModel = InputsModel, RunOptions = RunOptions, InputsCrit = InputsCrit, CalibOptions = CalibOptions, FUN_MOD = RunModel_GR2M, FUN_CRIT = ErrorCrit_NSE) OptimGR2M <- function(ParamOptim) {

Transformation of the parameter set to real space

RawParamOptim <- airGR::TransfoParam_GR2M(ParamIn = ParamOptim, Direction = "TR")

Simulation given a parameter set

OutputsModel <- airGR::RunModel_GR2M(InputsModel = InputsModel, RunOptions = RunOptions, Param = RawParamOptim)

Computation of the value of the performance criteria

OutputsCrit <- airGR::ErrorCritRMSE(InputsCrit = InputsCrit, OutputsModel = OutputsModel, verbose = FALSE) return(OutputsCrit$CritValue) } lowerGR2M <- rep(0, times = 2) upperGR2M <- rep(+999.99, times = 2) startGR2M <- c(500, 2) optPORT <- stats::nlminb(start = startGR2M, objective = OptimGR2M, lower = lowerGR2M, upper = upperGR2M, control = list(trace = 1)) startGR2M <- expand.grid(data.frame(CalibOptions$StartParamDistrib)) optPORT <- function(x) { opt <- stats::nlminb(start = x, objective = OptimGR2M, lower = lowerGR2M, upper = upperGR2M, control = list(trace = 1)) } listOptPORT <- apply(startGR2M, MARGIN = 1, FUN = optPORT_) parPORT <- t(sapply(listOptPORT, function(x) x$par)) objPORT <- sapply(listOptPORT, function(x) x$objective) resPORT <- data.frame(parPORT, RMSE = objPORT) summary(resPORT) library(airGR) library(DEoptim) library(hydroPSO) library(Rmalschains) optDE <- DEoptim::DEoptim(fn = OptimGR2M, lower = lowerGR2M, upper = upperGR2M, control = DEoptim::DEoptim.control(NP = 40, trace = 10)) optPSO <- hydroPSO::hydroPSO(fn = OptimGR2M, lower = lowerGR2M, upper = upperGR2M, control = list(write2disk = FALSE, verbose = FALSE)) optMALS <- Rmalschains::malschains(fn = OptimGR2M, lower = lowerGR2M, upper = upperGR2M, maxEvals = 2000) resGLOB <- data.frame(Algo = c("airGR", "PORT", "DE", "PSO", "MA-LS"), round(rbind( OutputsCalib$ParamFinalR
airGR::TransfoParam_GR2M(ParamIn = as.numeric(optDE$optim$bestmem), Direction = "TR"), airGR::TransfoParam_GR2M(ParamIn = as.numeric(optPSO$par) , Direction = "TR"), airGR::TransfoParam_GR2M(ParamIn = optMALS$sol , Direction = "TR")), digits = 3)) resGLOB

But I keep getting this error

Error in if (any(ProdStore < 0) | any(RoutStore < 0) | any(UH1 < 0) | : missing value where TRUE/FALSE needed Called from: CreateIniStates(FUN_MOD = RunModel_GR2M, InputsModel = InputsModel, ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L], ExpStore = NULL, UH1 = NULL, UH2 = NULL, GCemaNeigeLayers = NULL, eTGCemaNeigeLayers = NULL, verbose = FALSE) Browse[1]> Q

Can you please help?

[codingpanda19]

Hello, did you solve your problem?

[codingpanda19]

Hello, i'm not sure if i understand you well but if your issue is this one :

**> "Error in if (any(ProdStore < 0) | any(RoutStore < 0) | any(UH1 < 0) | :
> missing value where TRUE/FALSE needed
> Called from: CreateIniStates(FUN_MOD = RunModel_GR2M, InputsModel = InputsModel,
> ProdStore = RESULTS$StateEnd[1L], RoutStore = RESULTS$StateEnd[2L],
> ExpStore = NULL, UH1 = NULL, UH2 = NULL, GCemaNeigeLayers = NULL,
> eTGCemaNeigeLayers = NULL, verbose = FALSE)
> Browse[1]> Q"**

It means that probably you have very low value of Prodstore and Routstore, then you need to create your initialstates with function IniStates in airGR.

Your boundaries: **lowerGR2M <- rep(0, times = 2) upperGR2M <- rep(+999.99, times = 2) startGR2M <- c(500, 2)** are probably too much high, rember that those are Transformed values of parameter X1 (capacity of prod store) and X2 (capacity of routstore), you can try boundaries like this one:

**lowerGR2M <- rep(0, times = 2) upperGR2M <- rep(+9.99, times = 2)** Please let me know if anything...

[hzambran]

For the calibration of GR2M and other GR* models, you might want to have a look at: https://doi.org/10.5281/zenodo.3712382