diff --git a/tests/testthat/test-vignettes.R b/tests/testthat/test-vignettes.R
index 21a980206fe54a43b5df67391dd7dd1d0d1490f6..f17ee71b0b3b6eb7b9e5b4a784c9e177d3ba5916 100644
--- a/tests/testthat/test-vignettes.R
+++ b/tests/testthat/test-vignettes.R
@@ -9,9 +9,21 @@ test_that("V01_get_started works", {
 
 test_that("V02.1_param_optim works", {
   skip_on_cran()
-  skip("hydroPSO not working presently")
   rm(list = ls())
+  load(system.file("vignettesData/vignetteParamOptim.rda", package = "airGR"))
+  rda_resGLOB <- resGLOB
+  rda_resPORT <- resPORT
   expect_true(RunVignetteChunks("V02.1_param_optim"))
+  expect_equal(summary(resGLOB), summary(rda_resGLOB), tolerance = 1E-7)
+  resGLOB <- data.frame(Algo = c("airGR", "PORT", "DE", "PSO", "MA-LS"),
+    round(rbind(
+      OutputsCalib$ParamFinalR                          ,
+      airGR::TransfoParam_GR4J(ParamIn = optPORT$par                    , Direction = "TR"),
+      airGR::TransfoParam_GR4J(ParamIn = as.numeric(optDE$optim$bestmem), Direction = "TR"),
+      rda_resGLOB[4, c("X1", "X2", "X3", "X4")],
+      airGR::TransfoParam_GR4J(ParamIn = optMALS$sol                    , Direction = "TR")),
+      digits = 3))
+  expect_equal(resGLOB[,-1], rda_resGLOB[,-1], tolerance = 1E-2) # High tolerance due to randomisation in optimisations
 })
 
 test_that("V02.2_param_mcmc works", {
diff --git a/vignettes/V02.1_param_optim.Rmd b/vignettes/V02.1_param_optim.Rmd
index 1531c1c408682a741ba8bcb6b642a256ad60d703..c23151e00f2a370c4ea2f3efe82063260cf77450 100644
--- a/vignettes/V02.1_param_optim.Rmd
+++ b/vignettes/V02.1_param_optim.Rmd
@@ -13,7 +13,7 @@ vignette: >
 ```{r, warning=FALSE, include=FALSE, fig.keep='none', results='hide'}
 library(airGR)
 library(DEoptim)
-library(hydroPSO)
+library(hydroPSO) # Needs R version >= 3.6 or latticeExtra <= 0.6-28 on R 3.5
 library(Rmalschains)
 # source("airGR.R")
 set.seed(321)
@@ -84,7 +84,7 @@ upperGR4J <- rep(+9.99, times = 4)
 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:
 ```{r, warning=FALSE, results='hide', eval=FALSE}
 startGR4J <- c(4.1, 3.9, -0.9, -8.7)
-optPORT <- stats::nlminb(start = startGR4J, 
+optPORT <- stats::nlminb(start = startGR4J,
                          objective = OptimGR4J,
                          lower = lowerGR4J, upper = upperGR4J,
                          control = list(trace = 1))
@@ -97,7 +97,7 @@ For each starting point, a local optimization is performed.
 ```{r, warning=FALSE, results='hide', eval=FALSE}
 startGR4J <- expand.grid(data.frame(CalibOptions$StartParamDistrib))
 optPORT_ <- function(x) {
-  opt <- stats::nlminb(start = x, 
+  opt <- stats::nlminb(start = x,
                        objective = OptimGR4J,
                        lower = lowerGR4J, upper = upperGR4J,
                        control = list(trace = 1))
@@ -138,7 +138,7 @@ optDE <- DEoptim::DEoptim(fn = OptimGR4J,
 
 
 ## Particle Swarm
-```{r, warning=FALSE, results='hide', message=FALSE, eval=FALSE}
+```{r, warning=FALSE, results='hide', message=FALSE, eval=FALSE, purl=FALSE}
 optPSO <- hydroPSO::hydroPSO(fn = OptimGR4J,
                              lower = lowerGR4J, upper = upperGR4J,
                              control = list(write2disk = FALSE, verbose = FALSE))
@@ -147,7 +147,7 @@ optPSO <- hydroPSO::hydroPSO(fn = OptimGR4J,
 ## MA-LS-Chains
 ```{r, warning=FALSE, results='hide', eval=FALSE}
 optMALS <- Rmalschains::malschains(fn = OptimGR4J,
-                                   lower = lowerGR4J, upper = upperGR4J, 
+                                   lower = lowerGR4J, upper = upperGR4J,
                                    maxEvals = 2000)
 ```
 
@@ -155,8 +155,8 @@ optMALS <- Rmalschains::malschains(fn = OptimGR4J,
 
 As it can be seen in the table below, the four additional optimization strategies tested lead to very close optima.
 
-```{r, warning=FALSE, echo=FALSE, eval=FALSE}
-resGLOB <- data.frame(Algo = c("airGR", "PORT", "DE", "PSO", "MA-LS"), 
+```{r, warning=FALSE, echo=FALSE, eval=FALSE, purl=FALSE}
+resGLOB <- data.frame(Algo = c("airGR", "PORT", "DE", "PSO", "MA-LS"),
                       round(rbind(
                         OutputsCalib$ParamFinalR                          ,
                         airGR::TransfoParam_GR4J(ParamIn = optPORT$par                    , Direction = "TR"),