V02.1_param_optim.Rmd 7.71 KB
Newer Older
1
---
2
title: "Plugging in new calibration algorithms in airGR"
3
4
5
6
author: "François Bourgin"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteEngine{knitr::rmarkdown}
7
  %\VignetteIndexEntry{Plugging in new calibration algorithms}
8
9
10
11
12
13
14
15
16
17
18
  %\VignetteEncoding{UTF-8}
---



```{r, warning=FALSE, include=FALSE, fig.keep='none', results='hide'}
library(airGR)
library(DEoptim)
library(hydroPSO)
library(Rmalschains)
# source("airGR.R")
19
set.seed(321)
20
load(system.file("vignettesData/vignetteParamOptim.rda", package = "airGR"))
21
22
23
24
25
26
27
28
```



# Introduction

## Scope

29
30
The Michel's calibration strategy (`Calibration_Michel()` function) is the calibration algorithm proposed in **airGR**. However, other optimization methods can be used in combination with **airGR**.
We show here how to use different R packages to perform parameter estimation.
31
32

In this vignette, we use the **GR4J** model to illustrate the different optimization strategies.
33
In particular, we assume that the R global environment contains input climate data, observed discharge and functions from the [Get Started](V01_get_started.html) vignette, as shown below.
34
Please note that the calibration period is defined in the `CreateRunOptions()` function .
35

36
37
38
39
<!-- ```{r, warning=FALSE, fig.keep='none', results='hide', fig.height=10, fig.width=10, eval=TRUE, echo=FALSE, message=FALSE} -->
<!-- example("Calibration_Michel", echo = FALSE, ask = FALSE) -->
<!-- ``` -->

40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
```{r, echo=TRUE, eval=FALSE}
example("Calibration_Michel")
```


Regarding the different optimization strategies presented here, we refer to each package for in-depth information about the description of the methods used.

Please note that this vignette is only for illustration purposes and does not provide any guidance about which optimization strategies is recommended for the family of the **GR** models.

## Definition of the necessary function

Parameter estimation can be performed by defining a function that takes a parameter set as input and returns the value of the performance criterion.
There are two important steps: the transformation of parameters to real space and the computation of the value of the performance criterion.
Here we choose to minimize the root mean square error.

55
The change of the repository from the "real" parameter space to a "transformed" space ensures homogeneity of displacement in the different dimensions of the parameter space during the step-by-step procedure of the calibration algorithm of the model.
56
57


58
59
```{r, warning=FALSE, results='hide', eval=FALSE}
OptimGR4J <- function(ParamOptim) {
60
  ## Transformation of the parameter set to real space
61
62
  RawParamOptim <- airGR::TransfoParam_GR4J(ParamIn = ParamOptim,
                                            Direction = "TR")
63
64
65
  ## Simulation given a parameter set
  OutputsModel <- airGR::RunModel_GR4J(InputsModel = InputsModel,
                                       RunOptions = RunOptions,
66
                                       Param = RawParamOptim)
67
68
69
70
71
72
73
74
75
76
  ## Computation of the value of the performance criteria
  OutputsCrit <- airGR::ErrorCrit_RMSE(InputsCrit = InputsCrit,
                                       OutputsModel = OutputsModel,
                                       verbose = FALSE)
  return(OutputsCrit$CritValue)
}
```


In addition, we need to define the lower and upper bounds of the four **GR4J** parameters in the transformed parameter space:
77
```{r, warning=FALSE, results='hide', eval=FALSE}
78
79
80
81
82
83
84
lowerGR4J <- rep(-9.99, times = 4)
upperGR4J <- rep(+9.99, times = 4)
```

# Local optimization

We start with a local optimization strategy by using the PORT routines (using the `nlminb()` of the `stats` package) and by setting a starting point in the transformed parameter space:
85
```{r, warning=FALSE, results='hide', eval=FALSE}
86
87
startGR4J <- c(4.1, 3.9, -0.9, -8.7)
optPORT <- stats::nlminb(start = startGR4J, 
88
89
90
91
92
93
94
95
96
                         objective = OptimGR4J,
                         lower = lowerGR4J, upper = upperGR4J,
                         control = list(trace = 1))
```
The RMSE value reaches a local minimum value after 35 iterations.

We can also try a multi-start approach to test the consistency of the local optimization.
Here we use the same grid used for the filtering step of the Michel's calibration strategy (`Calibration_Michel()` function).
For each starting point, a local optimization is performed.
97
```{r, warning=FALSE, results='hide', eval=FALSE}
98
99
100
101
102
103
104
startGR4J <- expand.grid(data.frame(CalibOptions$StartParamDistrib))
optPORT_ <- function(x) {
  opt <- stats::nlminb(start = x, 
                       objective = OptimGR4J,
                       lower = lowerGR4J, upper = upperGR4J,
                       control = list(trace = 1))
}
105
listOptPORT <- apply(startGR4J, MARGIN = 1, FUN = optPORT_)
106
107
108
```

We can then extract the best parameter sets and the value of the performance criteria:
109
110
111
```{r, warning=FALSE, results='hide', eval=FALSE}
parPORT <- t(sapply(listOptPORT, function(x) x$par))
objPORT <- sapply(listOptPORT, function(x) x$objective)
112
resPORT <- data.frame(parPORT, RMSE = objPORT)
113
114
115
116
```

As can be seen below, the optimum performance criterion values (column *objective*) can differ from the global optimum value in many cases, resulting in various parameter sets.
```{r, warning=FALSE}
117
summary(resPORT)
118
119
120
121
122
123
124
125
126
127
```

The existence of several local minima illustrates the importance of defining an appropriate starting point or of using a multi-start strategy or a global optimization strategy.


# Global optimization

Global optimization is most often used when facing a complex response surface, with multiple local mimina.
Here we use the following R implementation of some popular strategies:

128
129
130
* [DEoptim: differential evolution](https://cran.r-project.org/package=DEoptim)
* [hydroPSO: particle swarm](https://cran.r-project.org/package=hydroPSO)
* [Rmalschains: memetic algorithms](https://cran.r-project.org/package=Rmalschains)
131
132

## Differential Evolution
133
```{r, warning=FALSE, results='hide', eval=FALSE}
134
135
136
137
138
139
140
optDE <- DEoptim::DEoptim(fn = OptimGR4J,
                          lower = lowerGR4J, upper = upperGR4J,
                          control = DEoptim::DEoptim.control(NP = 40, trace = 10))
```


## Particle Swarm
141
```{r, warning=FALSE, results='hide', message=FALSE, eval=FALSE}
142
143
144
145
146
147
optPSO <- hydroPSO::hydroPSO(fn = OptimGR4J,
                             lower = lowerGR4J, upper = upperGR4J,
                             control = list(write2disk = FALSE, verbose = FALSE))
```

## MA-LS-Chains
148
```{r, warning=FALSE, results='hide', eval=FALSE}
149
150
151
152
153
154
155
optMALS <- Rmalschains::malschains(fn = OptimGR4J,
                                   lower = lowerGR4J, upper = upperGR4J, 
                                   maxEvals = 2000)
```

# Results

156
As it can be seen in the table below, the four additional optimization strategies tested lead to very close optima.
157

158
```{r, warning=FALSE, echo=FALSE, eval=FALSE}
159
resGLOB <- data.frame(Algo = c("airGR", "PORT", "DE", "PSO", "MA-LS"), 
160
161
162
163
164
165
166
167
                      round(rbind(
                        OutputsCalib$ParamFinalR                          ,
                        airGR::TransfoParam_GR4J(ParamIn = optPORT$par                    , Direction = "TR"),
                        airGR::TransfoParam_GR4J(ParamIn = as.numeric(optDE$optim$bestmem), Direction = "TR"),
                        airGR::TransfoParam_GR4J(ParamIn = as.numeric(optPSO$par)         , Direction = "TR"),
                        airGR::TransfoParam_GR4J(ParamIn = optMALS$sol                    , Direction = "TR")),
                        digits = 3))
```
168
```{r, warning=FALSE, echo=FALSE}
169
resGLOB
170
171
172
173
174
175
```

<!-- This is an expected result because the response surface for quadratic performance criteria of the **GR4J** model is generally sufficiently well defined in the transformed parameter space to allow using a local optimization strategy instead of a more time consuming global optimization strategy. -->