diff --git a/DESCRIPTION b/DESCRIPTION
index 9b38f1f6..1bfb3d1d 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -36,14 +36,17 @@ Imports:
     Matrix,
     BayesianTools,
     crayon,
+    DEoptim,
     doParallel,
     dplyr,
     ggplot2,
+    graDiEnt,
     lhs,
     lifecycle,
     lubridate,
     methods,
     nloptr,
+    parma,
     purrr,
     rlang,
     rstudioapi,
@@ -76,4 +79,4 @@ Encoding: UTF-8
 Language: en-US
 LazyData: true
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.3.1
+RoxygenNote: 7.3.3
diff --git a/NAMESPACE b/NAMESPACE
index 165dc471..3ac4784a 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -14,7 +14,10 @@ export(likelihood_log_ciidn_corr)
 export(plot_estimVSinit)
 export(plot_valuesVSit)
 export(plot_valuesVSit_2D)
+export(plot_valuesVSit_2D_go)
+export(plot_valuesVSit_go)
 export(test_wrapper)
+export(wrap_cmaes)
 importFrom(BayesianTools,MAP)
 importFrom(BayesianTools,applySettingsDefault)
 importFrom(BayesianTools,correlationPlot)
@@ -24,6 +27,7 @@ importFrom(BayesianTools,gelmanDiagnostics)
 importFrom(BayesianTools,getSample)
 importFrom(BayesianTools,marginalPlot)
 importFrom(BayesianTools,runMCMC)
+importFrom(DEoptim,DEoptim)
 importFrom(dplyr,"%>%")
 importFrom(dplyr,filter)
 importFrom(dplyr,select)
@@ -45,7 +49,12 @@ importFrom(ggplot2,scale_y_log10)
 importFrom(ggplot2,theme)
 importFrom(ggplot2,xlim)
 importFrom(ggplot2,ylim)
+importFrom(graDiEnt,GetAlgoParams)
+importFrom(graDiEnt,optim_SQGDE)
+importFrom(parma,cmaes)
+importFrom(parma,cmaes.control)
 importFrom(purrr,modify)
+importFrom(rgenoud,genoud)
 importFrom(rlang,.data)
 importFrom(stats,setNames)
 importFrom(tibble,tibble)
diff --git a/R/global_optim_function.R b/R/global_optim_function.R
new file mode 100644
index 00000000..e57e161c
--- /dev/null
+++ b/R/global_optim_function.R
@@ -0,0 +1,472 @@
+#' @title Summarizes results of global optim methods
+#'
+#' @inheritParams estim_param
+#'
+#' @param optim_results Results list returned by global optim method wrappers
+#'
+#' @return Prints results of global optim methods
+#'
+summary_global_optim <- function(optim_options, param_info, optim_results, out_dir) {
+  param_names <- get_params_names(param_info)
+  nb_params <- length(param_names)
+  final_values <- optim_results$final_values
+  min_crit_value <- optim_results$min_crit_value
+  est_values <- matrix(final_values, nrow = 1, dimnames = list(NULL, param_names))
+
+  cat(paste(
+    "\nList of observed variables used:",
+    paste(optim_results$obs_var_list, collapse = ", ")
+  ))
+  # Display of parameters values for the candidate which has the
+  # smallest criterion
+  for (ipar in 1:nb_params) {
+    cat(paste(
+      "\nEstimated value for", param_names[ipar], ": ",
+      format(final_values[[ipar]],
+        scientific = FALSE,
+        digits = 2, nsmall = 0
+      )
+    ))
+  }
+  cat(paste(
+    "\nMinimum value of the criterion:",
+    format(min_crit_value, scientific = FALSE, digits = 2, nsmall = 0)
+  ))
+  cat(paste(
+    "\nComplementary graphs and results can be found in ", out_dir,
+    "\n"
+  ))
+}
+
+#' @title Post-treat results of global optim methods
+#'
+#' @inheritParams estim_param
+#'
+#' @param optim_results Results list returned by global optim method wrappers
+#' @param crit_options List containing several arguments given to `estim_param`
+#'  function: `param_names`, `obs_list`, `model_function`,
+#'  `model_options`, `param_info`, `transform_obs`, `transform_sim`
+#' that must be passed to main_crit function by the methods wrappers.
+#'
+#' @return Updated results of global optim method
+#'
+
+post_treat_global_optim <- function(optim_options, param_info, optim_results,
+                                    crit_options) {
+  param_names <- get_params_names(param_info)
+  nb_params <- length(param_names)
+  info_crit_list <- crit_options$info_crit_list
+
+  # Recompute final value of minimized criterion
+  # (just to check it is correct and to get the observation list used)
+  info_final <- main_crit(
+    param_values = optim_results$final_values,
+    crit_options = c(crit_options, return_detailed_info = TRUE)
+  )
+  if (info_final$crit != optim_results$min_crit_value) {
+    stop(paste(
+      "Internal error: incoherent computation of minimum criterion value. \nValue obtained in method wrapper:",
+      optim_results$min_crit_value, "\nValue obtained afterwards:",
+      info_final$crit
+    ))
+  }
+  optim_results$forced_param_values <- info_final$forced_param_values
+  sapply(info_crit_list, function(x) {
+    final_info_crit <- x(
+      obs_list = info_final$obs_intersect,
+      crit = info_final$crit,
+      param_nb = nb_params
+    )
+    optim_results[x()$name] <<- final_info_crit
+  })
+  return(optim_results)
+}
+
+
+#' @title Generate plots for global optim methods
+#'
+#' @inheritParams estim_param
+#'
+#' @param optim_results Results list returned by global optim method wrappers
+#'
+#' @return Returns the list of plots + save them in a pdf file.
+#'
+#' @keywords internal
+#'
+plot_global_optim <- function(optim_options, param_info, optim_results, out_dir) {
+  est_values <- optim_results$est_values
+  trace_df <- optim_results$trace_df
+  p_all <- list()
+
+  # ValuesVSit
+
+  if (!is.null(trace_df)) {
+    tryCatch(
+      {
+        grDevices::pdf(
+          file = file.path(out_dir, "ValuesVSit_globaloptim.pdf"),
+          width = 9, height = 9
+        )
+      },
+      error = function(cond) {
+        filename <- paste0("ValuesVSit_new.pdf")
+        warning(
+          "Error trying to create ", out_dir,
+          "/ValuesVSit.pdf file. It is maybe opened in a pdf viewer and locked. It will be created under the name ",
+          filename
+        )
+        grDevices::pdf(
+          file = file.path(out_dir, filename),
+          width = 9, height = 9
+        )
+      }
+    )
+
+    # trace_df : params, ind, iter, crit, eval, rep, method
+    p <- plot_valuesVSit_go(
+      df = trace_df,
+      param_info = param_info,
+      iter_or_eval = "iter",
+      crit_log = FALSE,
+      ind_label = "begin"
+    )
+    for (plot in p) print(plot)
+    grDevices::dev.off()
+    p_all$valuesVSit <- p
+  }
+
+
+  #  ValuesVSit_2D
+
+  if (!is.null(trace_df)) {
+    df_2d <- trace_df
+
+    tryCatch(
+      {
+        grDevices::pdf(
+          file = file.path(out_dir, "ValuesVSit_2D.pdf"),
+          width = 9, height = 9
+        )
+      },
+      error = function(cond) {
+        filename <- paste0("ValuesVSit_2D_new.pdf")
+        warning(
+          "Error trying to create ", out_dir,
+          "/ValuesVSit_2D.pdf file. It is maybe opened in a pdf viewer and locked. It will be created under the name ",
+          filename
+        )
+        grDevices::pdf(
+          file = file.path(out_dir, filename),
+          width = 9, height = 9
+        )
+      }
+    )
+
+    p <- plot_valuesVSit_2D_go(
+      df = df_2d,
+      param_info = param_info,
+      iter_or_eval = "eval",
+      fill = "crit",
+      crit_log = FALSE,
+      lines = FALSE,
+      rep_label = "end"
+    )
+
+    for (plot in p) print(plot)
+    grDevices::dev.off()
+    p_all$valuesVSit_2D <- p
+  }
+
+  return(p_all)
+}
+
+
+#' @title Create plots of parameters and criterion values per iteration or
+#' evaluation number
+#'
+#' @param param_info List describing the parameters to estimate
+
+#' @param df Data.frame containing values of each individual and iteration parameters values (one column per
+#' estimated parameter), criterion (crit column), individual index (ind),
+#' iteration number (iter) and evaluation number (eval)
+#'  (similar to params_and_crit).
+#' See Details section for comments about the difference between evaluations
+#' and iterations.
+#' @param iter_or_eval Values of the x axis: "iter" for iteration number,
+#' "eval" for evaluation number
+#' @param crit_log If TRUE, consider criterion values in log scale
+#' @param ind_label Indicate if labels for the individual number must be
+#' plotted at both beginning and end of lines ("begin_end"),
+#' only at the beginning ("begin") or only at the end ("end")
+#'
+#' @return A named list containing one plot per parameter and a plot for the
+#' criterion.
+#'
+#' @details Evaluation means evaluation of the criterion from proposed values of
+#' the parameters by the parameter estimation algorithm.
+#' An iteration is reached when an evaluation lead to a better value of the
+#' criterion than the previously obtained values.
+#' There are thus more evaluations than iterations. The criterion decreases when
+#' iteration number increases while it is not the case when evaluation number
+#' increases.
+#'
+#' @importFrom ggplot2 ggplot aes_string theme element_text geom_point
+#' scale_color_gradient2 geom_line geom_label aes labs scale_y_log10
+#' @importFrom dplyr select filter %>%
+#'
+#' @export
+#'
+
+plot_valuesVSit_go <- function(df, param_info, iter_or_eval = c("iter", "eval"),
+                               crit_log = TRUE,
+                               ind_label = c("begin_end", "begin", "end")) {
+  param_names <- get_params_names(param_info)
+  bounds <- get_params_bounds(param_info)
+
+  lab <- "evaluations"
+  if (iter_or_eval[1] == "iter") {
+    df <- filter(df, !is.na(.data$iter))
+    lab <- "iterations"
+  }
+  have_crit <- "crit" %in% names(df) && !all(is.na(df$crit))
+  trans <- "identity"
+  mid <- NA_real_
+  if (have_crit) {
+    mid <- (max(df$crit) - min(df$crit)) / 2 + min(df$crit)
+    if (crit_log) {
+      if (all(df$crit > 0)) {
+        trans <- "log10"
+        mid <- (max(log10(df$crit)) -
+          +min(log10(df$crit))) / 2 + min(log10(df$crit))
+      } else {
+        warning("The criterion takes negative values, log transformation will not be done.")
+        crit_log <- FALSE
+      }
+    }
+  } else {
+    crit_log <- FALSE
+  }
+  tmp <- rbind(bounds$lb, bounds$ub, select(df, all_of(param_names)))
+  tmp[tmp == Inf | tmp == -Inf] <- NA
+  minvalue <- apply(tmp, 2, min, na.rm = TRUE)
+  maxvalue <- apply(tmp, 2, max, na.rm = TRUE)
+  minvalue <- minvalue - 0.05 * (maxvalue - minvalue)
+  maxvalue <- maxvalue + 0.05 * (maxvalue - minvalue)
+
+  p <- list()
+
+  for (param_name in param_names) {
+    p[[param_name]] <- ggplot(
+      df, aes(
+        x = !!rlang::sym(iter_or_eval[1]),
+        y = !!rlang::sym(param_name),
+        color = crit
+      )
+    ) +
+      labs(
+        title = paste0(
+          "Evolution of ", param_name,
+          " \n in function of the minimization ", lab
+        ),
+        y = param_name,
+        x = paste(lab, "number"),
+        fill = "Criterion"
+      ) +
+      theme(plot.title = element_text(hjust = 0.5)) +
+      geom_point(alpha = 0.5)
+
+    if (have_crit) {
+      p[[param_name]] <- p[[param_name]] +
+        scale_color_gradient2(
+          midpoint = mid, low = "blue", mid = "yellow",
+          high = "red", space = "Lab", trans = trans
+        )
+    }
+
+
+    p[[param_name]] <- p[[param_name]] +
+      ylim(minvalue[param_name], maxvalue[param_name])
+  }
+
+  df$ind <- as.factor(df$ind)
+  p[["criterion"]] <- ggplot(df, aes(
+    x = !!rlang::sym(iter_or_eval[1]),
+    y = crit,
+    color = crit
+  )) +
+    labs(
+      title = paste0(
+        "Evolution of the minimized criterion \n in function of the minimization ",
+        lab
+      ),
+      y = "Minimized criterion",
+      x = paste(lab, "number"),
+      color = "Criterion"
+    ) +
+    theme(plot.title = element_text(hjust = 0.5)) +
+    geom_point(alpha = 0.6)
+
+  if (have_crit) {
+    p[["criterion"]] <- p[["criterion"]] +
+      scale_color_gradient2(
+        midpoint = mid, low = "blue", mid = "yellow",
+        high = "red", space = "Lab", trans = trans
+      )
+  }
+
+  if (crit_log) {
+    p[["criterion"]] <- p[["criterion"]] + scale_y_log10()
+  }
+
+  if (crit_log) {
+    p[["criterion"]] <- p[["criterion"]] + scale_y_log10()
+  }
+
+  return(p)
+}
+
+
+#' @title Create 2D plots of parameters values evolution per iteration or
+#' evaluation number
+#'
+#' @inheritParams estim_param
+#' @param df Data.frame containing values of parameters (one column per
+#' estimated parameter), criterion (crit column), repetition number (rep),
+#' iteration number (iter) and evaluation number (eval)
+#'  (similar to params_and_crit).
+#' See Details section for comments about the difference between evaluations
+#' and iterations.
+#' @param iter_or_eval "iter" for plotting the values for each iteration,
+#' "eval" for plotting the values for each evaluation
+#' @param fill If "crit", colours the points and lines in function of the
+#' minimized criterion value, if "rep" colours in function of the
+#' repetition number.
+#' @param crit_log If TRUE, consider criterion values in log scale
+#' @param lines If TRUE add lines between points of a same repetition
+#' @param rep_label Indicate if labels for the repetition number must be plotted
+#' at both beginning and end of lines ("begin_end"), only at the beginning
+#' ("begin") or only at the end ("end")
+#'
+#' @return A list containing one plot per parameter pair.
+#'
+#' @details Evaluation means evaluation of the criterion from proposed values of
+#' the parameters by the parameter estimation algorithm.
+#' An iteration is reached when an evaluation lead to a better value of the
+#' criterion than the previously obtained values.
+#' There are thus more evaluations than iterations. The criterion decreases when
+#' iteration number increases while it is not the case when evaluation number
+#' increases.
+#'
+#' @importFrom ggplot2 ggplot aes_string theme element_text geom_point labs
+#' xlim ylim geom_path scale_y_log10
+#' @importFrom dplyr select filter %>%
+#'
+#' @export
+#'
+plot_valuesVSit_2D_go <- function(df, param_info, iter_or_eval = c("eval", "iter"),
+                                  fill = c("crit", "rep"), crit_log = TRUE,
+                                  lines = FALSE,
+                                  rep_label = c("begin_end", "begin", "end")) {
+  param_names <- get_params_names(param_info)
+  if (length(param_names) <= 1) {
+    return()
+  }
+  bounds <- get_params_bounds(param_info)
+
+  lab <- "evaluations"
+  if (iter_or_eval[1] == "iter") {
+    df <- filter(df, !is.na(.data$iter))
+    lab <- "iterations"
+  }
+  if ("rep" %in% names(df)) {
+    df$rep <- as.factor(df$rep)
+  }
+  has_crit <- "crit" %in% names(df)
+  trans <- "identity"
+  mid <- NA_real_
+  if (fill[1] == "crit" && has_crit) {
+    mid <- (max(df$crit) - min(df$crit)) / 2 + min(df$crit)
+    if (crit_log) {
+      if (all(df$crit > 0)) {
+        trans <- "log10"
+        mid <- (max(log10(df$crit)) -
+          min(log10(df$crit))) / 2 + min(log10(df$crit))
+      } else {
+        warning("The criterion takes negative values, log transformation will not be done.")
+        crit_log <- FALSE
+      }
+    }
+  } else {
+    crit_log <- FALSE
+  }
+  tmp <- rbind(bounds$lb, bounds$ub, select(df, all_of(param_names)))
+  # check ...
+  tmp[tmp == Inf | tmp == -Inf] <- NA
+  minvalue <- apply(tmp, 2, min, na.rm = TRUE)
+  maxvalue <- apply(tmp, 2, max, na.rm = TRUE)
+  minvalue <- minvalue - 0.05 * (maxvalue - minvalue)
+  maxvalue <- maxvalue + 0.05 * (maxvalue - minvalue)
+
+  p <- list()
+
+  df_pairs <- utils::combn(param_names, 2)
+
+  for (ipair in seq_len(ncol(df_pairs))) {
+    p[[ipair]] <- ggplot(df, aes_string(
+      x = df_pairs[1, ipair],
+      y = df_pairs[2, ipair], color = fill[1]
+    )) +
+      labs(
+        title = paste0(
+          "Evolution of ", df_pairs[1, ipair], " and ",
+          df_pairs[2, ipair], " \n in function of the minimization ",
+          lab
+        ),
+        y = paste("Estimated value for", df_pairs[2, ipair]),
+        x = paste("Estimated value for", df_pairs[1, ipair]),
+        fill = "Criterion"
+      ) +
+      theme(plot.title = element_text(hjust = 0.5)) +
+      geom_point(alpha = 0.5)
+
+    if (fill[1] == "crit" && has_crit) {
+      p[[ipair]] <- p[[ipair]] +
+        scale_color_gradient2(
+          midpoint = mid, low = "blue", mid = "yellow",
+          high = "red", space = "Lab", trans = trans
+        )
+    } else if (fill[1] == "eval") {
+      p[[ipair]] <- p[[ipair]] +
+        scale_color_gradient2(
+          low  = "lightblue",
+          high = "darkblue"
+        )
+    }
+
+    if (lines) {
+      for (irep in unique(df$rep)) {
+        p[[ipair]] <- p[[ipair]] +
+          geom_path(data = filter(df, rep == irep))
+        if (rep_label[1] == "begin_end" || rep_label[1] == "begin") {
+          p[[ipair]] <- p[[ipair]] +
+            geom_label(aes(label = rep),
+              data = filter(df, rep == irep) %>% filter(eval == min(.data$eval)),
+              size = 3
+            )
+        }
+        if (rep_label[1] == "begin_end" || rep_label[1] == "end") {
+          p[[ipair]] <- p[[ipair]] +
+            geom_label(aes(label = rep),
+              data = filter(df, rep == irep) %>% filter(eval == max(.data$eval)),
+              size = 3
+            )
+        }
+      }
+    }
+    p[[ipair]] <- p[[ipair]] +
+      xlim(minvalue[df_pairs[1, ipair]], maxvalue[df_pairs[1, ipair]]) +
+      ylim(minvalue[df_pairs[2, ipair]], maxvalue[df_pairs[2, ipair]])
+  }
+
+  return(p)
+}
diff --git a/R/multi_steps_functions.R b/R/multi_steps_functions.R
index 967a2520..61f741ce 100644
--- a/R/multi_steps_functions.R
+++ b/R/multi_steps_functions.R
@@ -1,7 +1,5 @@
 #' @title Summarizes results of multi-step procedure
 #'
-#' @inheritParams estim_param
-#'
 #' @param results_multi_step Results of the multi_step procedure as returned by post_treat_multi_step
 #'
 #' @param path_results Folder path where results of the multi-step optimization process can be found
@@ -45,8 +43,6 @@ summary_multi_step <- function(results_multi_step, path_results) {
 
 #' @title Post-treat results of multi-step procedure
 #'
-#' @inheritParams estim_param
-#'
 #' @param optim_results_list List of results returned for each step of the parameter estimation process
 #'
 #' @return List of estimated and forced parameters values
diff --git a/R/optim_switch.R b/R/optim_switch.R
index 2e19e32b..dd05c51b 100644
--- a/R/optim_switch.R
+++ b/R/optim_switch.R
@@ -16,14 +16,19 @@
 #'
 
 optim_switch <- function(...) {
+  message("optim_switch nargs = ", nargs())
+  message("names(arguments) = ", paste(names(list(...)), collapse = ", "))
+
   # Store and save results at the end of the current optimization step
+  res <- list()
+  flag_error <- FALSE
   if (nargs() > 2) {
-    # Initialize res
-    res <- list()
-    flag_error <- FALSE
-
     on.exit({
-      res$forced_param_values <- crit_options$forced_param_values
+      ####
+      if (!is.null(crit_options)) {
+        res$forced_param_values <- crit_options$forced_param_values
+      }
+      #####
       if (!is.null(res$final_values)) {
         res$forced_param_values <- compute_eq_const(
           res$forced_param_values,
@@ -68,7 +73,8 @@ optim_switch <- function(...) {
         }
       }
 
-      if (!is.null(crit_options$out_dir) & length(res) > 0) {
+      # if (!is.null(crit_options$out_dir) & length(res) > 0) {
+      if (!is.null(crit_options) && !is.null(crit_options$out_dir) && length(res) > 0) {
         save(res, file = file.path(
           crit_options$out_dir,
           "optim_results.Rdata"
@@ -91,6 +97,19 @@ optim_switch <- function(...) {
   optim_method <- arguments$optim_method
   optim_options <- arguments$optim_options
   wrap_args <- within(arguments, rm("optim_method"))
+  ##########################
+  #
+  # param_info   <- NULL
+  # crit_options <- NULL
+
+  # if ("param_info" %in% names(arguments)) {
+  #   param_info <- arguments$param_info
+  # }
+  # if ("crit_options" %in% names(arguments)) {
+  #  crit_options <- arguments$crit_options
+  # }
+  ##########################
+
   if (nargs() > 2) {
     param_info <- arguments$param_info
     crit_options <- arguments$crit_options
@@ -132,7 +151,8 @@ optim_switch <- function(...) {
         res$obs_situation_list <- .croptEnv$obs_situation_list
         res$plots <- plot_bayesian(
           optim_options = optim_options,
-          param_info = param_info, optim_results = res,
+          param_info = param_info,
+          optim_results = res,
           out_dir = crit_options$out_dir
         )
         summary_bayesian(
@@ -141,37 +161,122 @@ optim_switch <- function(...) {
           out_dir = crit_options$out_dir
         )
       }
-    } else if (optim_method == "optim") {
-      res <- do.call(wrap_optim, wrap_args)
+    } else if (optim_method == "DEoptim") {
+      res <- do.call(wrap_DEoptim, wrap_args)
       if (nargs() > 2) {
         res$obs_var_list <- .croptEnv$obs_var_list
         res$obs_situation_list <- .croptEnv$obs_situation_list
         if (arguments$crit_options$info_level >= 1) {
           res$params_and_crit <- dplyr::bind_rows(.croptEnv$params_and_crit)
         }
-        res <- post_treat_frequentist(
+        res <- post_treat_global_optim(
           optim_options = optim_options,
           param_info = param_info,
           optim_results = res,
           crit_options = crit_options
         )
-        res$plots <- plot_frequentist(
+        res$plots <- plot_global_optim(
           optim_options = optim_options,
-          param_info = param_info, optim_results = res,
+          param_info = param_info,
+          optim_results = res,
           out_dir = crit_options$out_dir
         )
-        summary_frequentist(
+        summary_global_optim(
           optim_options = optim_options, param_info = param_info,
           optim_results = res,
           out_dir = crit_options$out_dir
         )
       }
+    } else if (optim_method == "graDiEnt") {
+      res <- do.call(wrap_graDiEnt, wrap_args)
+      if (nargs() > 2) {
+        res$obs_var_list <- .croptEnv$obs_var_list
+        res$obs_situation_list <- .croptEnv$obs_situation_list
+        if (arguments$crit_options$info_level >= 1) {
+          res$params_and_crit <- dplyr::bind_rows(.croptEnv$params_and_crit)
+        }
+        res <- post_treat_global_optim(
+          optim_options = optim_options,
+          param_info = param_info,
+          optim_results = res,
+          crit_options = crit_options
+        )
+        res$plots <- plot_global_optim(
+          optim_options = optim_options,
+          param_info = param_info,
+          optim_results = res,
+          out_dir = crit_options$out_dir
+        )
+        summary_global_optim(
+          optim_options = optim_options,
+          param_info = param_info,
+          optim_results = res,
+          out_dir = crit_options$out_dir
+        )
+      }
+    }else if (optim_method == "genoud") {
+      res <- do.call(wrap_rgenoud, wrap_args)
+      if (nargs() > 2) {
+        res$obs_var_list <- .croptEnv$obs_var_list
+        res$obs_situation_list <- .croptEnv$obs_situation_list
+        if (arguments$crit_options$info_level >= 1) {
+          res$params_and_crit <- dplyr::bind_rows(.croptEnv$params_and_crit)
+        }
+        res <- post_treat_global_optim(
+          optim_options = optim_options,
+          param_info = param_info,
+          optim_results = res,
+          crit_options = crit_options
+        )
+        res$plots <- plot_global_optim(
+          optim_options = optim_options,
+          param_info = param_info,
+          optim_results = res,
+          out_dir = crit_options$out_dir
+        )
+        summary_global_optim(
+          optim_options = optim_options,
+          param_info = param_info,
+          optim_results = res,
+          out_dir = crit_options$out_dir
+        )
+      }
+    } else if (optim_method == "cmaes") {
+      wrap_args$optim_options <- optim_options
+      res <- do.call(wrap_cmaes, wrap_args)
+
+      if (nargs() > 2) {
+        res$obs_var_list <- .croptEnv$obs_var_list
+        res$obs_situation_list <- .croptEnv$obs_situation_list
+        if (arguments$crit_options$info_level >= 1) {
+          res$params_and_crit <- dplyr::bind_rows(.croptEnv$params_and_crit)
+        }
+        res <- post_treat_global_optim(
+          optim_options = optim_options,
+          param_info = param_info,
+          optim_results = res,
+          crit_options = crit_options
+        )
+        res$plots <- plot_global_optim(
+          optim_options = optim_options,
+          param_info = param_info,
+          optim_results = res,
+          out_dir = crit_options$out_dir
+        )
+        summary_global_optim(
+          optim_options = optim_options,
+          param_info = param_info,
+          optim_results = res,
+          out_dir = crit_options$out_dir
+        )
+      }
     } else {
       flag_unknown_method <- TRUE
     },
     error = function(cond) {
-      warning(cond)
-      flag_error <<- TRUE
+      message("=== ERROR CAUGHT IN optim_switch ===")
+      traceback(20)
+      stop(cond)
     }
   )
 
@@ -179,7 +284,7 @@ optim_switch <- function(...) {
     flag_error <- TRUE
     stop(paste0(
       "Unknown method ", optim_method,
-      ", please choose between nloptr.simplex, BayesianTools.dreamzs and optim."
+      ", please choose between nloptr.simplex, BayesianTools.dreamzs, DEoptim and graDiEnt."
     ))
   }
 
diff --git a/R/wrap_DE.R b/R/wrap_DE.R
new file mode 100644
index 00000000..a21d6482
--- /dev/null
+++ b/R/wrap_DE.R
@@ -0,0 +1,137 @@
+#' @title A wrapper for DEoptim package
+#'
+#' @inheritParams optim_switch
+#'
+#' @return prints, graphs and a list containing:
+#' `final_values`, the vector of estimated values for optimized parameters
+#' for the repetition that lead to the lowest value of the criterion
+#' `init_values`, the vector of initial values for optimized parameters
+#' `est_values`, data frame of values of all population
+#' `ind_min_crit`, the index of the repetition that lead to the lowest value
+#' of the criterion
+#' `DE`, the data structure returned by DE
+#' `crit_values` the vector of criterion values for the all population
+#' @importFrom DEoptim DEoptim
+#' @keywords internal
+
+wrap_DEoptim <- function(optim_options, param_info, crit_options) {
+  if (is.null((ranseed <- optim_options$ranseed))) {
+    ranseed <- NULL
+  }
+  if (is.null((reltol <- optim_options$reltol))) {
+    optim_options$reltol <- 1e-10
+  }
+  if (is.null((NP <- optim_options$NP))) {
+    stop(
+      "Optim_options$NP is mandatory, please define it (e.g., 10 * number of parameters)"
+    )
+  }
+  if (is.null((trace <- optim_options$trace))) {
+    optim_options$trace <- FALSE
+  }
+  optim_options$ranseed <- NULL # ranseed is set to NULL because DEoptim.controll doesn't regonise it
+  optim_options$storepopfrom <- 1
+  optim_options$storepopfreq <- 1
+  control <- do.call(DEoptim::DEoptim.control, optim_options)
+  algorithm <- "DEoptim"
+
+  # return requested information if only optim_options is given in argument
+  if (nargs() == 1 & methods::hasArg(optim_options)) {
+    return(list(
+      package = "DEoptim", family = "Global",
+      method = "DEoptim", init_values_nb = control$NP
+    ))
+  }
+  param_names <- get_params_names(param_info)
+  nb_params <- length(param_names)
+
+  crit_options$tot_max_eval <- control$NP * control$itermax
+  bounds <- get_params_bounds(param_info)
+  init_values <- get_init_values(param_info)
+  control$initialpop <- as.matrix(init_values)
+
+
+  start_time <- Sys.time()
+
+  .trace_env <- new.env(parent = emptyenv())
+  .trace_env$x_list <- list()
+  .trace_env$crit_list <- list()
+  .trace_env$k <- 0L
+
+  if (!is.null(ranseed)) set.seed(ranseed)
+  fn_de <- function(x) {
+    val <- main_crit(x, crit_options = crit_options)
+    .trace_env$k <- .trace_env$k + 1L
+    .trace_env$x_list[[.trace_env$k]] <- as.numeric(x)
+    .trace_env$crit_list[[.trace_env$k]] <- val
+    return(val)
+  }
+  DE <- tryCatch(
+    DEoptim::DEoptim(
+      fn = fn_de,
+      lower = bounds$lb,
+      upper = bounds$ub,
+      control = control
+    ),
+    error = function(e) {
+      warning(sprintf("DEoptim failed: %s", e$message))
+      NULL
+    }
+  )
+  elapsed <- Sys.time() - start_time
+
+  # Verify criterion value
+  if (is.na(DE$optim$bestval)) {
+    stop(
+      "The value of the criterion is NA."
+    )
+  }
+
+
+  # Get the estimated values  ==> matrix of NP * nb_params
+  est_values <- NULL
+  est_values <- as.matrix(DE$member$pop)
+  colnames(est_values) <- param_names
+  if (ncol(est_values) == nb_params) {
+    colnames(est_values) <- param_names
+  }
+
+  # Final solution
+  final_values <- DE$optim$bestmem
+  names(final_values) <- param_names
+
+  trace_df <- NULL
+  if (!is.null(DE$member$storepop)) {
+    storepop <- DE$member$storepop
+    itermax <- length(storepop)
+    NP <- nrow(storepop[[1]])
+
+    pop_all <- do.call(rbind, lapply(seq_len(itermax), function(it) {
+      m <- as.data.frame(storepop[[it]])
+      colnames(m) <- param_names
+      m$ind <- seq_len(NP)
+      m$iter <- it
+      m
+    }))
+    crit_vec <- as.numeric(unlist(.trace_env$crit_list))
+    n <- min(nrow(pop_all), length(crit_vec))
+    pop_all <- pop_all[seq_len(n), , drop = FALSE]
+    pop_all$crit <- crit_vec[seq_len(n)]
+
+    pop_all$eval <- seq_len(n)
+    pop_all$rep <- 1L
+    pop_all$method <- "DEoptim"
+    trace_df <- pop_all
+  }
+
+  res <- list(
+    final_values = final_values,
+    init_values = init_values,
+    est_values = est_values,
+    min_crit_value = DE$optim$bestval,
+    DE = DE,
+    trace_df = trace_df
+  )
+
+  return(res)
+}
diff --git a/R/wrap_cmaes.R b/R/wrap_cmaes.R
new file mode 100644
index 00000000..0fb278cf
--- /dev/null
+++ b/R/wrap_cmaes.R
@@ -0,0 +1,214 @@
+#' Wrap CMA-ES using parma
+#'
+#' @inheritParams optim_switch
+#'
+#' @return list with final_values, init_values, est_values, GE, min_crit_value
+#' @importFrom parma cmaes cmaes.control
+#' @keywords internal
+#' @export
+wrap_cmaes <- function(optim_options, param_info = NULL, crit_options) {
+  if (!requireNamespace("parma", quietly = TRUE)) {
+    stop("the package 'parma' (not installed) becouse it required a dependencies that doesn't exist in the cluster")
+  }
+  if (is.null((ranseed <- optim_options$ranseed))) {
+    ranseed <- NULL
+  }
+  n_params <- optim_options$n_params
+  if (is.null(n_params)) {
+    n_params <- length(param_info)
+  }
+  if (is.null((trace <- optim_options$return_trace))) {
+    optim_options$return_trace <- TRUE
+  }
+  optim_options$ranseed <- NULL # ranseed is set to NULL because getAlgoParams doesn't regonise it
+  algorithm <- "cmaes"
+
+  if (is.null(param_info) || is.null(crit_options)) {
+    return(list(
+      package = "parma", family = "Global",
+      method = "cmaes", init_values_nb = init_nb <- if (!is.null(optim_options$ctrl$options$PopSize)) optim_options$ctrl$options$PopSize else NA_integer_
+    ))
+  }
+  param_names <- get_params_names(param_info)
+  nb_params <- length(param_names)
+  bounds <- get_params_bounds(param_info)
+  init_values <- get_init_values(param_info)
+
+  control_params <- optim_options$ctrl
+
+
+  ctrl <- parma::cmaes.control()
+
+
+  # if ctrl = list(options=..., CMA=...) existe
+  if (!is.null(control_params) && is.list(control_params) &&
+    (!is.null(control_params$options) || !is.null(control_params$CMA))) {
+    if (!is.null(control_params$options)) {
+      for (nm in names(control_params$options)) {
+        ctrl$options[[nm]] <- control_params$options[[nm]]
+      }
+    }
+
+    if (!is.null(control_params$CMA)) {
+      for (nm in names(control_params$CMA)) {
+        ctrl$CMA[[nm]] <- control_params$CMA[[nm]]
+      }
+    }
+
+    # si l'utilisateur a fourni directement une liste d'options (sans options/CMA)
+  } else if (!is.null(control_params) && is.list(control_params)) {
+    for (nm in names(control_params)) {
+      ctrl$options[[nm]] <- control_params[[nm]]
+    }
+  }
+
+  if (!is.null(ctrl$options$PopSize) && !is.null(ctrl$options$MaxIter)) {
+    crit_options$tot_max_eval <- as.integer(ctrl$options$PopSize * ctrl$options$MaxIter)
+  }
+
+  start_time <- Sys.time()
+
+  if (!is.null(ranseed)) set.seed(ranseed)
+  trace_env <- new.env(parent = emptyenv())
+  trace_env$eval <- 0L
+  trace_env$crit <- numeric(0)
+  trace_env$X <- list()
+  keep_trace <- isTRUE(optim_options$return_trace)
+  ObjFun <- function(x) {
+    x <- as.numeric(x)
+    names(x) <- param_names
+    val <- main_crit(x, crit_options = crit_options)
+
+    if (keep_trace) {
+      trace_env$eval <- trace_env$eval + 1L
+      trace_env$crit[trace_env$eval] <- val
+      trace_env$X[[trace_env$eval]] <- x
+    }
+    return(val)
+  }
+  insigma <- optim_options$insigma
+  if (is.null(insigma)) {
+    # default value if it is not provided
+    insigma <- rep(0.3, nb_params)
+  }
+  if (length(insigma) == 1L) insigma <- rep(insigma, nb_params)
+
+  # because parma::cmaes expects a numeric vector for 'pars' (length = nb_params)
+  x0 <- init_values
+
+  # init_values can be a matrix/data.frame with several rows (nb_values x nb_params) our case data frame
+  if (is.data.frame(x0)) x0 <- as.matrix(x0)
+
+  if (is.matrix(x0)) {
+    # take the first row as starting point
+    x0 <- x0[1, , drop = TRUE]
+  }
+
+  x0 <- as.numeric(x0)
+  names(x0) <- param_names
+
+
+  if (length(x0) != nb_params) {
+    x0 <- (bounds$lb + bounds$ub) / 2
+    names(x0) <- param_names
+  }
+
+  parma_res <- tryCatch(
+    parma::cmaes(
+      pars    = x0,
+      fun     = ObjFun,
+      lower   = bounds$lb,
+      upper   = bounds$ub,
+      insigma = insigma,
+      ctrl    = ctrl
+    ),
+    error = function(e) {
+      warning(sprintf("parma::cmaes failed: %s", e$message))
+      NULL
+    }
+  )
+  if (isTRUE(optim_options$debug)) {
+    message("parma_res fields: ", paste(names(parma_res), collapse = ", "))
+  }
+
+
+  elapsed <- Sys.time() - start_time
+
+  if (is.null(parma_res)) {
+    stop("parma::cmaes failed (object NULL).")
+  }
+
+  sol <- NULL
+  if (!is.null(parma_res$bestever$x)) {
+    sol <- parma_res$bestever$x
+  } else if (!is.null(parma_res$par)) {
+    sol <- parma_res$par
+  } else if (!is.null(parma_res$xbest)) {
+    sol <- parma_res$xbest
+  }
+  if (is.null(sol)) {
+    stop("Impossible to find a solution in the object sent by parma::cmaes().")
+  }
+  final_values <- as.numeric(sol)
+  names(final_values) <- param_names
+
+  # for post_treat_global_optim()
+  min_crit_value <- if (!is.null(parma_res$bestever$f)) {
+    parma_res$bestever$f
+  } else {
+    main_crit(final_values, crit_options)
+  }
+
+
+  cmaes_bestever_f <- if (!is.null(parma_res$bestever$f)) parma_res$bestever$f else NA_real_
+  cmaes_objective <- if (!is.null(parma_res$objective)) parma_res$objective else NA_real_
+
+
+  # Getting the estimated values
+  est_values <- matrix(
+    final_values,
+    nrow = 1L,
+    dimnames = list(NULL, param_names)
+  )
+  colnames(est_values) <- param_names
+  if (ncol(est_values) == nb_params) {
+    colnames(est_values) <- param_names
+  }
+
+  # Final solution
+  trace_df <- NULL
+  if (keep_trace && trace_env$eval > 0L) {
+    Xmat <- do.call(rbind, trace_env$X)
+    colnames(Xmat) <- param_names
+
+    trace_df <- as.data.frame(Xmat)
+    trace_df$crit <- trace_env$crit
+    trace_df$eval <- seq_len(trace_env$eval)
+    lambda <- ctrl$options$PopSize
+    if (is.null(lambda) || is.na(lambda)) {
+      lambda <- 30L
+    }
+    lambda <- as.integer(lambda)
+
+    trace_df$iter <- ((trace_df$eval - 1L) %/% lambda) + 1L # génération CMA-ES
+    trace_df$ind <- ((trace_df$eval - 1L) %% lambda) + 1L # individu dans la pop            # single trajectory
+    trace_df$rep <- 1L
+    trace_df$method <- "cmaes"
+  }
+
+  res <- list(
+    final_values = final_values,
+    init_values = init_values,
+    est_values = est_values,
+    min_crit_value = min_crit_value,
+    CMAES = parma_res,
+    trace_df = trace_df,
+    cmaes_bestever_f = cmaes_bestever_f,
+    cmaes_objective = cmaes_objective,
+    counteval = if (!is.null(parma_res$counteval)) parma_res$counteval else NA_integer_,
+    stopflag = if (!is.null(parma_res$stopflag)) parma_res$stopflag else NA_character_,
+    out_dir = crit_options$out_dir
+  )
+
+  return(res)
+}
diff --git a/R/wrap_graDiEnt.R b/R/wrap_graDiEnt.R
new file mode 100644
index 00000000..282652f4
--- /dev/null
+++ b/R/wrap_graDiEnt.R
@@ -0,0 +1,161 @@
+#' @title A wrapper for graDiEnt package
+#'
+#' @inheritParams optim_switch
+#'
+#' @return list with final_values, init_values, est_values, GE, min_crit_value
+#' @importFrom graDiEnt optim_SQGDE GetAlgoParams
+#' @keywords internal
+#'
+wrap_graDiEnt <- function(optim_options, param_info, crit_options) {
+  if (is.null((ranseed <- optim_options$ranseed))) {
+    ranseed <- NULL
+  }
+  if (is.null((n_params <- optim_options$n_params))) {
+    stop(
+      "Optim_options$n_params is mandatory, please define it (e.g., 10 * number of parameters)"
+    )
+  }
+  if (is.null((trace <- optim_options$return_trace))) {
+    optim_options$return_trace <- TRUE
+  }
+  optim_options$ranseed <- NULL # ranseed is set to NULL because getAlgoParams doesn't regonise it
+  algorithm <- "graDiEnt"
+
+  # return requested information if only optim_options is given in argument
+  if (nargs() == 1 && methods::hasArg(optim_options)) {
+    init_nb <- if (!is.null(optim_options$n_particles)) {
+      optim_options$n_particles
+    } else {
+      NA_integer_
+    }
+    return(list(
+      package = "graDiEnt", family = "Global",
+      method = "SQGDE", init_values_nb = init_nb
+    ))
+  }
+  param_names <- get_params_names(param_info)
+  nb_params <- length(param_names)
+  bounds <- get_params_bounds(param_info)
+  init_values <- get_init_values(param_info)
+
+  control_params <- do.call(graDiEnt::GetAlgoParams, optim_options)
+  crit_options$tot_max_eval <- control_params$n_particles * control_params$n_iter
+
+  start_time <- Sys.time()
+
+  .trace_env <- new.env(parent = emptyenv())
+  .trace_env$par_list <- list()
+  .trace_env$crit_list <- list()
+  .trace_env$k <- 0L
+
+  if (!is.null(ranseed)) set.seed(ranseed)
+  range_bounds <- bounds$ub - bounds$lb
+  ObjFun <- function(u) {
+    u <- as.numeric(u)
+    names(u) <- param_names
+    x <- bounds$lb + u * range_bounds
+    names(x) <- param_names
+    val <- main_crit(x, crit_options = crit_options)
+    .trace_env$k <- .trace_env$k + 1L
+    .trace_env$par_list[[.trace_env$k]] <- x
+    .trace_env$crit_list[[.trace_env$k]] <- val
+    return(val)
+  }
+  SQGDE <- tryCatch(
+    graDiEnt::optim_SQGDE(
+      ObjFun         = ObjFun,
+      control_params = control_params
+    ),
+    error = function(e) {
+      warning(sprintf("GraDiEnt failed: %s", e$message))
+      NULL
+    }
+  )
+  elapsed <- Sys.time() - start_time
+
+  # Verify criterion value
+  if (is.null(SQGDE) || is.null(SQGDE$solution)) {
+    stop(
+      "The value of the criterion is NA."
+    )
+  }
+
+  u_sol <- SQGDE$solution
+  names(u_sol) <- param_names
+
+  final_values <- bounds$lb + u_sol * range_bounds
+  names(final_values) <- param_names
+  min_crit_value <- SQGDE$solution_weight
+  if (is.null(min_crit_value)) {
+    min_crit_value <- main_crit(final_values, crit_options)
+  }
+  # Get the estimated values  ==> matrix of 1 * nb_params
+  est_values <- matrix(
+    final_values,
+    nrow = 1L,
+    dimnames = list(NULL, param_names)
+  )
+  colnames(est_values) <- param_names
+  if (ncol(est_values) == nb_params) {
+    colnames(est_values) <- param_names
+  }
+
+  # Final solution
+  if (!is.null(SQGDE$particles_trace)) {
+    tr <- SQGDE$particles_trace
+    if (length(dim(tr)) != 3) {
+      warning("Unexpected dimension for SQGDE$particles_trace; est_values set to final_values only.")
+      est_values <- matrix(
+        final_values,
+        nrow = 1L,
+        dimnames = list(NULL, param_names)
+      )
+    } else {
+      n_it <- dim(tr)[1]
+      last_pop <- tr[n_it, , ]
+      est_u <- as.matrix(last_pop)
+      colnames(est_u) <- param_names
+
+      # transformation
+      est_values <- sweep(est_u, 2, range_bounds, `*`)
+      est_values <- sweep(est_values, 2, bounds$lb, `+`)
+      colnames(est_values) <- param_names
+    }
+  } else {
+    est_values <- matrix(
+      final_values,
+      nrow = 1L,
+      dimnames = list(NULL, param_names)
+    )
+  }
+  trace_df <- NULL
+  if (.trace_env$k > 0L) {
+    n_pop <- control_params$n_particles
+    k <- .trace_env$k
+    pars_mat <- do.call(rbind, lapply(.trace_env$par_list, function(v) as.numeric(v)))
+    colnames(pars_mat) <- param_names
+    crit_vec <- as.numeric(unlist(.trace_env$crit_list))
+
+    n_it_est <- ceiling(k / n_pop)
+    iter <- rep(seq_len(n_it_est), each = n_pop)[seq_len(k)]
+    ind <- rep(seq_len(n_pop), times = n_it_est)[seq_len(k)]
+
+    trace_df <- as.data.frame(pars_mat)
+    trace_df$ind <- ind
+    trace_df$iter <- iter
+    trace_df$crit <- crit_vec
+    trace_df$eval <- seq_len(k)
+    trace_df$method <- "graDiEnt"
+    trace_df$rep <- 1L
+  }
+  res <- list(
+    final_values = final_values,
+    init_values = init_values,
+    est_values = est_values,
+    min_crit_value = min_crit_value,
+    SQGDE = SQGDE,
+    trace_df = trace_df
+  )
+
+  return(res)
+}
diff --git a/R/wrap_parma.R b/R/wrap_parma.R
new file mode 100644
index 00000000..43920b06
--- /dev/null
+++ b/R/wrap_parma.R
@@ -0,0 +1,136 @@
+#' @title A wrapper for parma::cmaes package (CMA-ES)
+#'
+#' @inheritParams optim_switch
+#'
+#' @return list with final_values, init_values, est_values, GE, min_crit_value
+#' @importFrom parma cmaes cmaes.control
+#' @keywords internal
+#'
+wrap_cmaes <- function(optim_options, param_info, crit_options) {
+  if(is.null((ranseed <- optim_options$ranseed))) {
+    ranseed = NULL
+  }
+  if (is.null((n_params <- optim_options$n_params))) {
+    stop(
+      "Optim_options$n_params is mandatory, please define it (e.g., 10 * number of parameters)"
+    )
+  }
+  if (is.null((trace <- optim_options$return_trace))) {
+    optim_options$return_trace <- TRUE
+  }
+  optim_options$ranseed <- NULL # ranseed is set to NULL because getAlgoParams doesn't regonise it
+  algorithm <- "parma_cmaes"
+
+  # return requested information if only optim_options is given in argument
+  if (nargs() == 1 && methods::hasArg(optim_options)) {
+    init_nb <- if (!is.null(optim_options$n_particles)) {
+      optim_options$n_particles
+    } else {
+      NA_integer_
+    }
+    return(list(
+      package = "parma", family = "Global",
+      method = "cmaes", init_values_nb = init_nb
+    ))
+  }
+  param_names <- get_params_names(param_info)
+  nb_params <- length(param_names)
+  bounds <- get_params_bounds(param_info)
+  init_values <- get_init_values(param_info)
+
+  control_params <- optim_options$ctrl
+  if (is.null(control_params )) {
+    control_params <- parma::cmaes.control()
+   } else{
+     if (is.null(control_params$options) && is.null(control_params$CMA)) {
+       ctrl <- parma::cmaes.control(options = control_params, CMA = list())
+       }else {
+         ctrl <- do.call(parma::cmaes.control, control_params)
+       }
+   }
+  if (!is.null(ctrl$options$PopSize) && !is.null(ctrl$options$MaxIter)){
+    crit_options$tot_max_eval <- as.integer(ctrl$options$PopSize * ctrl$options$MaxIter)
+    }
+
+  start_time <- Sys.time()
+
+  if (!is.null(ranseed)) set.seed(ranseed)
+  ObjFun <- function(x){
+    x <- as.numeric(x)
+    names(x) <- param_names
+    main_crit(x, crit_options = crit_options)
+    }
+  insigma <- optim_options$insigma
+  if (is.null(insigma)) {
+  # default value if it is not provided
+  insigma <- rep(0.3, nb_params)
+  }
+  if (length(insigma) == 1L) insigma <- rep(insigma, nb_params)
+  parma_res <- tryCatch(
+    parma::cmaes(
+    pars    = init_values,
+    fun     = ObjFun,
+    lower   = bounds$lb,
+    upper   = bounds$ub,
+    insigma = insigma,
+    ctrl    = ctrl
+    ),
+    error = function(e) {
+    warning(sprintf("parma::cmaes failed: %s", e$message))
+     NULL
+    }
+    )
+  elapsed <- Sys.time() - start_time
+
+  if (is.null(parma_res)) {
+    stop("parma::cmaes failed (object NULL).")
+  }
+
+  sol <- NULL
+  if (!is.null(parma_res$par)) sol <- parma_res$par
+  if (is.null(sol) && !is.null(parma_res$bestever$x))
+    sol <- parma_res$bestever$x
+  if (is.null(sol) && !is.null(parma_res$xbest))
+    sol <- parma_res$xbest
+  if (is.null(sol)) {
+    stop("Impossible to find a solution in the object sent by parma::cmaes().")
+    }
+
+  final_values <- as.numeric(sol)
+  names(final_values) <- param_names
+
+  min_crit_value <- NULL
+  if (!is.null(parma_res$value))
+    min_crit_value <- parma_res$value
+  if (is.null(min_crit_value) && !is.null(parma_res$bestever$f))
+    min_crit_value <- parma_res$bestever$f
+  if (is.null(min_crit_value)) {
+    min_crit_value <- main_crit(final_values, crit_options = crit_options)
+  }
+
+  # Get the estimated values  ==> matrix of 1 * nb_params
+  est_values <- matrix(
+    final_values,
+    nrow = 1L,
+    dimnames = list(NULL, param_names)
+  )
+  colnames(est_values) <- param_names
+  if (ncol(est_values) == nb_params) {
+    colnames(est_values) <- param_names
+  }
+
+  # Final solution
+  trace_df <- NULL
+
+  res <- list(
+    final_values = final_values,
+    init_values = init_values,
+    est_values = est_values,
+    min_crit_value = min_crit_value,
+    CMAES = parma_res,
+    trace_df = trace_df
+  )
+
+  return(res)
+}
+
diff --git a/R/wrap_rgenoud.R b/R/wrap_rgenoud.R
new file mode 100644
index 00000000..c01efda1
--- /dev/null
+++ b/R/wrap_rgenoud.R
@@ -0,0 +1,218 @@
+#' @title A wrapper for rgenoud package (single criterion)
+#' @inheritParams optim_switch
+#' @importFrom rgenoud genoud
+#' @return A list containing:
+#'   `final_values`: vector of estimated parameter values
+#'   corresponding to the best solution.
+#'   `init_values`: vector or matrix of initial parameter values.
+#'   `est_values`: matrix of estimated values (final solution).
+#'   `min_crit_value`: minimum value of the criterion.
+#'   `trace_df`: data.frame containing parameters, criterion,
+#'   iteration, evaluation and individual indices .
+#'   `obs_var_list`: list of observed variables used by the criterion.
+#'   `elapsed`: computation time.
+#'   `genoud.pro`: an output file called located in the temporary
+#'   directory provided by tempdir, containing solutions in each
+#'   generation (criterion and parameter values )
+#'
+#' @keywords internal
+
+wrap_rgenoud <- function(optim_options, param_info, crit_options) {
+  if (nargs() == 1 && methods::hasArg(optim_options)) {
+    init_TP <- if (!is.null(optim_options$pop.size))
+      optim_options$pop.size
+    else
+      NA_integer_
+    return(list(
+      package = "rgenoud", family = "Global",
+      method = "genoud", init_values_nb = init_TP
+    ))
+    }
+  if (is.null(optim_options$pop.size)) {
+    stop("optim_options$pop.size is mandatory, please define it (e.g., 10 * number of parameters)")
+    }
+
+  ranseed <- optim_options$ranseed
+  optim_options$ranseed <- NULL
+
+  param_names <- get_params_names(param_info)
+  nb_params <- length(param_names)
+
+  bounds <- get_params_bounds(param_info)
+  init_values <- get_init_values(param_info)
+
+  Domains <- cbind(bounds$lb, bounds$ub)
+
+  # read final population from genoud.pro to construct est_value
+  read_last_population_from_pro <- function(project_path, nb_params, param_names) {
+    if (is.null(project_path) || !file.exists(project_path))
+      return(NULL)
+    L <- readLines(project_path, warn = FALSE)
+    gen_idx <- grep("^Generation:\\s*[0-9]+", L)
+    if (length(gen_idx) == 0)
+      return(NULL)
+    start <- gen_idx[length(gen_idx)]   # last generation
+    block <- L[start:length(L)]
+    rows <- grep("^\\s*[0-9]+\\s+", block, value = TRUE)
+    if (length(rows) == 0)
+      return(NULL)
+    tok <- strsplit(trimws(rows), "\\s+")
+    mat <- do.call(rbind, lapply(tok, function(x) as.numeric(x)))
+    if (ncol(mat) < (2 + nb_params))  #(je vérifie si ce n est pas 3+nb_params)
+      return(NULL)
+    crit <- mat[,2]
+    pars <- mat[, 3:(2+nb_params), drop = FALSE]  # ici aussi je dois vérifier!!
+    colnames(pars) <- param_names
+    list(pars = pars, crit = crit)
+    }
+
+  starting.values <- init_values
+  if (is.vector(starting.values)) {
+    starting.values <- matrix(starting.values, nrow = 1)
+    }
+  if (is.data.frame(starting.values)) {
+    starting.values <- as.matrix(starting.values)
+    }
+  if (!is.null(ranseed)) {
+    set.seed(ranseed)
+    if (is.null(optim_options$unif.seed))
+      optim_options$unif.seed <- sample.int(.Machine$integer.max, 1)
+    if (is.null(optim_options$int.seed))
+      optim_options$int.seed  <- sample.int(.Machine$integer.max, 1)
+    }
+  # project file
+  if (is.null(optim_options$project.path)) {
+    optim_options$project.path <- file.path(
+      crit_options$out_dir,
+      "genoud.pro"
+    )
+    }
+  # tot_max_eval : genoud could do more (local search!!)
+  crit_options$tot_max_eval <- optim_options$pop.size * optim_options$max.generations
+
+  # trace env
+  .trace_env <- new.env(parent = emptyenv())
+  .trace_env$x_list <- list()
+  .trace_env$crit_list <- list()
+  .trace_env$k <- 0
+
+  fn_gen <- function(x) {
+    val <- main_crit(x, crit_options = crit_options)
+    .trace_env$k <- .trace_env$k + 1
+    .trace_env$x_list[[.trace_env$k]] <- as.numeric(x)
+    .trace_env$crit_list[[.trace_env$k]] <- val
+    val
+    }
+
+  start_time <- Sys.time()
+
+  GENOUD <- tryCatch(
+    do.call(
+      rgenoud::genoud,
+      c(
+        list(
+          fn = fn_gen,
+          nvars = nb_params,
+          max = FALSE,
+          Domains = Domains,
+          starting.values = starting.values
+          ),
+        optim_options
+        )
+      ),
+    error = function(e) {
+      warning(sprintf("genoud failed: %s", e$message))
+      NULL
+      }
+    )
+
+  elapsed <- Sys.time() - start_time
+
+  if (is.null(GENOUD)) {
+    return(NULL)
+  }
+
+  # final solution
+  final_values <- as.numeric(GENOUD$par)
+  names(final_values) <- param_names
+  min_crit_value <- GENOUD$value
+
+  # est_values : the best individual / final population en amont !!
+  est_values <- matrix(final_values, nrow = 1, dimnames = list(NULL, param_names))
+  # trace_df reconstruction from saved file
+  trace_df <- NULL
+  if (.trace_env$k > 0L) {
+    X <- do.call(rbind, .trace_env$x_list)     # k x nb_params
+    critv <- as.numeric(unlist(.trace_env$crit_list))
+    eval <- seq_along(critv)
+    popsize_used <- optim_options$pop.size
+
+
+    # iter = batch (pas génération genoud ), ind = index dans le batch
+    iter <- ((eval - 1L) %/% popsize_used) + 1L  #à ver
+    ind  <- ((eval - 1L) %%  popsize_used) + 1L
+
+    colnames(X) <- param_names
+    trace_df <- data.frame(
+      iter = as.integer(iter),
+      ind  = as.integer(ind),
+      as.data.frame(X),
+      crit = critv,
+      eval = as.integer(eval),
+      rep = 1L,
+      method = "rgenoud"
+    )
+    # est_values from genoud.pro (final gener)
+    pro_pop <- read_last_population_from_pro(
+      project_path = optim_options$project.path,
+      nb_params = nb_params,
+      param_names = param_names
+      )
+    if (!is.null(pro_pop)) {
+      est_values <- as.matrix(pro_pop$pars)
+      }
+    }
+  #obs_var_list
+  info_final <- tryCatch(
+    main_crit(
+      param_values = final_values,
+      crit_options = c(crit_options, return_detailed_info = TRUE)
+    ),
+    error = function(e) NULL
+    )
+
+  obs_var_list <- NULL
+  if (!is.null(info_final)) {
+    if (!is.null(info_final$obs_var_list))
+      obs_var_list <- info_final$obs_var_list
+    else if (!is.null(info_final$obs_intersect))
+      obs_var_list <- names(info_final$obs_intersect)
+    }
+  # Stop info
+  stop_reason <- "HARD MAXIMUM GENERATION LIMIT HIT"
+  if (!is.null(GENOUD$hard.generation.limit) && isTRUE(GENOUD$hard.generation.limit)) {
+    stop_reason <- "hard maximum generation limit hit"
+    }
+  if (!is.null(GENOUD$message) && nzchar(GENOUD$message)) {
+    stop_reason <- GENOUD$message
+    }
+  generations_run <- NA_integer_
+  if (!is.null(GENOUD$generations))
+    generations_run <- as.integer(GENOUD$generations)
+
+  res <- list(
+    final_values = final_values,
+    init_values = init_values,
+    est_values = est_values,
+    min_crit_value = min_crit_value,
+    GENOUD = GENOUD,
+    trace_df = trace_df,
+    obs_var_list = obs_var_list,
+    elapsed = elapsed,
+    stop_reason = stop_reason,
+    generations_run = generations_run,
+    project_path = optim_options$project.path
+  )
+
+  return(res)
+}
diff --git a/man/plot_global_optim.Rd b/man/plot_global_optim.Rd
new file mode 100644
index 00000000..2f2d32f5
--- /dev/null
+++ b/man/plot_global_optim.Rd
@@ -0,0 +1,53 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/global_optim_function.R
+\name{plot_global_optim}
+\alias{plot_global_optim}
+\title{Generate plots for global optim methods}
+\usage{
+plot_global_optim(optim_options, param_info, optim_results, out_dir)
+}
+\arguments{
+\item{optim_options}{List of options of the parameter estimation method, containing:
+\itemize{
+\item \code{ranseed} Set random seed so that each execution of estim_param give the same
+results when using the same seed. If you want randomization, set it to NULL,
+otherwise set it to a number of your choice (e.g. 1234) (optional, default to NULL, which means random seed)
+\item specific options depending on the method used. Click on the links to see examples with the \href{https://sticsrpacks.github.io/CroptimizR/articles/Parameter_estimation_simple_case.html}{simplex}
+and \href{https://sticsrpacks.github.io/CroptimizR/articles/Parameter_estimation_DREAM.html}{DreamZS} methods.
+\item \code{out_dir} \ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#deprecated}{\figure{lifecycle-deprecated.svg}{options: alt='[Deprecated]'}}}{\strong{[Deprecated]}} Definition of \code{out_dir} in \code{optim_options} is no longer supported, use the new argument \code{out_dir} of \code{estim_param} instead.
+}}
+
+\item{param_info}{Information on the parameters to estimate.
+Either
+a list containing:
+\itemize{
+\item \code{ub} and \code{lb}, named vectors of upper and lower bounds (-Inf and Inf can be used if init_values is provided),
+\item \code{default}, named vectors of default values (optional, corresponding parameters are set to these values when the parameter is part of the \code{candidate_param} list and when it is not estimated ; these values are also used as first initial values when the parameters are estimated)
+\item \code{init_values}, a data.frame containing initial
+values to test for the parameters (optional, if not provided, or if less values
+than number of repetitions of the minimization are provided, the, or part
+of the, initial values will be randomly generated using LHS sampling within
+parameter bounds).
+}
+
+or a named list containing for each parameter:
+\itemize{
+\item \code{sit_list}, list the groups of situations for which the current estimated
+parameter must take different values (see \href{https://sticsrpacks.github.io/CroptimizR/articles/Parameter_estimation_Specific_and_Varietal.html}{here}
+for an example),
+\item \code{ub} and \code{lb}, vectors of upper and lower bounds (one value per group),
+\item \code{init_values}, the list of initial values per group  (data.frame, one column per group, optional).
+\item \code{default}, vector of default values per group (optional, the parameter is set to its default value when it is part of the \code{candidate_param} list and when it is not estimated ; the default value is also used as first initial value when the parameter is estimated)
+}}
+
+\item{optim_results}{Results list returned by global optim method wrappers}
+
+\item{out_dir}{Path to the directory where the optimization results will be written. (optional, default to \code{getwd()})}
+}
+\value{
+Returns the list of plots + save them in a pdf file.
+}
+\description{
+Generate plots for global optim methods
+}
+\keyword{internal}
diff --git a/man/plot_valuesVSit_2D_go.Rd b/man/plot_valuesVSit_2D_go.Rd
new file mode 100644
index 00000000..0ee03ced
--- /dev/null
+++ b/man/plot_valuesVSit_2D_go.Rd
@@ -0,0 +1,79 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/global_optim_function.R
+\name{plot_valuesVSit_2D_go}
+\alias{plot_valuesVSit_2D_go}
+\title{Create 2D plots of parameters values evolution per iteration or
+evaluation number}
+\usage{
+plot_valuesVSit_2D_go(
+  df,
+  param_info,
+  iter_or_eval = c("eval", "iter"),
+  fill = c("crit", "rep"),
+  crit_log = TRUE,
+  lines = FALSE,
+  rep_label = c("begin_end", "begin", "end")
+)
+}
+\arguments{
+\item{df}{Data.frame containing values of parameters (one column per
+estimated parameter), criterion (crit column), repetition number (rep),
+iteration number (iter) and evaluation number (eval)
+(similar to params_and_crit).
+See Details section for comments about the difference between evaluations
+and iterations.}
+
+\item{param_info}{Information on the parameters to estimate.
+Either
+a list containing:
+\itemize{
+\item \code{ub} and \code{lb}, named vectors of upper and lower bounds (-Inf and Inf can be used if init_values is provided),
+\item \code{default}, named vectors of default values (optional, corresponding parameters are set to these values when the parameter is part of the \code{candidate_param} list and when it is not estimated ; these values are also used as first initial values when the parameters are estimated)
+\item \code{init_values}, a data.frame containing initial
+values to test for the parameters (optional, if not provided, or if less values
+than number of repetitions of the minimization are provided, the, or part
+of the, initial values will be randomly generated using LHS sampling within
+parameter bounds).
+}
+
+or a named list containing for each parameter:
+\itemize{
+\item \code{sit_list}, list the groups of situations for which the current estimated
+parameter must take different values (see \href{https://sticsrpacks.github.io/CroptimizR/articles/Parameter_estimation_Specific_and_Varietal.html}{here}
+for an example),
+\item \code{ub} and \code{lb}, vectors of upper and lower bounds (one value per group),
+\item \code{init_values}, the list of initial values per group  (data.frame, one column per group, optional).
+\item \code{default}, vector of default values per group (optional, the parameter is set to its default value when it is part of the \code{candidate_param} list and when it is not estimated ; the default value is also used as first initial value when the parameter is estimated)
+}}
+
+\item{iter_or_eval}{"iter" for plotting the values for each iteration,
+"eval" for plotting the values for each evaluation}
+
+\item{fill}{If "crit", colours the points and lines in function of the
+minimized criterion value, if "rep" colours in function of the
+repetition number.}
+
+\item{crit_log}{If TRUE, consider criterion values in log scale}
+
+\item{lines}{If TRUE add lines between points of a same repetition}
+
+\item{rep_label}{Indicate if labels for the repetition number must be plotted
+at both beginning and end of lines ("begin_end"), only at the beginning
+("begin") or only at the end ("end")}
+}
+\value{
+A list containing one plot per parameter pair.
+}
+\description{
+Create 2D plots of parameters values evolution per iteration or
+evaluation number
+}
+\details{
+Evaluation means evaluation of the criterion from proposed values of
+the parameters by the parameter estimation algorithm.
+An iteration is reached when an evaluation lead to a better value of the
+criterion than the previously obtained values.
+There are thus more evaluations than iterations. The criterion decreases when
+iteration number increases while it is not the case when evaluation number
+increases.
+}
diff --git a/man/plot_valuesVSit_go.Rd b/man/plot_valuesVSit_go.Rd
new file mode 100644
index 00000000..7fa7e99b
--- /dev/null
+++ b/man/plot_valuesVSit_go.Rd
@@ -0,0 +1,51 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/global_optim_function.R
+\name{plot_valuesVSit_go}
+\alias{plot_valuesVSit_go}
+\title{Create plots of parameters and criterion values per iteration or
+evaluation number}
+\usage{
+plot_valuesVSit_go(
+  df,
+  param_info,
+  iter_or_eval = c("iter", "eval"),
+  crit_log = TRUE,
+  ind_label = c("begin_end", "begin", "end")
+)
+}
+\arguments{
+\item{df}{Data.frame containing values of each individual and iteration parameters values (one column per
+estimated parameter), criterion (crit column), individual index (ind),
+iteration number (iter) and evaluation number (eval)
+(similar to params_and_crit).
+See Details section for comments about the difference between evaluations
+and iterations.}
+
+\item{param_info}{List describing the parameters to estimate}
+
+\item{iter_or_eval}{Values of the x axis: "iter" for iteration number,
+"eval" for evaluation number}
+
+\item{crit_log}{If TRUE, consider criterion values in log scale}
+
+\item{ind_label}{Indicate if labels for the individual number must be
+plotted at both beginning and end of lines ("begin_end"),
+only at the beginning ("begin") or only at the end ("end")}
+}
+\value{
+A named list containing one plot per parameter and a plot for the
+criterion.
+}
+\description{
+Create plots of parameters and criterion values per iteration or
+evaluation number
+}
+\details{
+Evaluation means evaluation of the criterion from proposed values of
+the parameters by the parameter estimation algorithm.
+An iteration is reached when an evaluation lead to a better value of the
+criterion than the previously obtained values.
+There are thus more evaluations than iterations. The criterion decreases when
+iteration number increases while it is not the case when evaluation number
+increases.
+}
diff --git a/man/post_treat_global_optim.Rd b/man/post_treat_global_optim.Rd
new file mode 100644
index 00000000..16900d51
--- /dev/null
+++ b/man/post_treat_global_optim.Rd
@@ -0,0 +1,55 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/global_optim_function.R
+\name{post_treat_global_optim}
+\alias{post_treat_global_optim}
+\title{Post-treat results of global optim methods}
+\usage{
+post_treat_global_optim(optim_options, param_info, optim_results, crit_options)
+}
+\arguments{
+\item{optim_options}{List of options of the parameter estimation method, containing:
+\itemize{
+\item \code{ranseed} Set random seed so that each execution of estim_param give the same
+results when using the same seed. If you want randomization, set it to NULL,
+otherwise set it to a number of your choice (e.g. 1234) (optional, default to NULL, which means random seed)
+\item specific options depending on the method used. Click on the links to see examples with the \href{https://sticsrpacks.github.io/CroptimizR/articles/Parameter_estimation_simple_case.html}{simplex}
+and \href{https://sticsrpacks.github.io/CroptimizR/articles/Parameter_estimation_DREAM.html}{DreamZS} methods.
+\item \code{out_dir} \ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#deprecated}{\figure{lifecycle-deprecated.svg}{options: alt='[Deprecated]'}}}{\strong{[Deprecated]}} Definition of \code{out_dir} in \code{optim_options} is no longer supported, use the new argument \code{out_dir} of \code{estim_param} instead.
+}}
+
+\item{param_info}{Information on the parameters to estimate.
+Either
+a list containing:
+\itemize{
+\item \code{ub} and \code{lb}, named vectors of upper and lower bounds (-Inf and Inf can be used if init_values is provided),
+\item \code{default}, named vectors of default values (optional, corresponding parameters are set to these values when the parameter is part of the \code{candidate_param} list and when it is not estimated ; these values are also used as first initial values when the parameters are estimated)
+\item \code{init_values}, a data.frame containing initial
+values to test for the parameters (optional, if not provided, or if less values
+than number of repetitions of the minimization are provided, the, or part
+of the, initial values will be randomly generated using LHS sampling within
+parameter bounds).
+}
+
+or a named list containing for each parameter:
+\itemize{
+\item \code{sit_list}, list the groups of situations for which the current estimated
+parameter must take different values (see \href{https://sticsrpacks.github.io/CroptimizR/articles/Parameter_estimation_Specific_and_Varietal.html}{here}
+for an example),
+\item \code{ub} and \code{lb}, vectors of upper and lower bounds (one value per group),
+\item \code{init_values}, the list of initial values per group  (data.frame, one column per group, optional).
+\item \code{default}, vector of default values per group (optional, the parameter is set to its default value when it is part of the \code{candidate_param} list and when it is not estimated ; the default value is also used as first initial value when the parameter is estimated)
+}}
+
+\item{optim_results}{Results list returned by global optim method wrappers}
+
+\item{crit_options}{List containing several arguments given to \code{estim_param}
+function: \code{param_names}, \code{obs_list}, \code{model_function},
+\code{model_options}, \code{param_info}, \code{transform_obs}, \code{transform_sim}
+that must be passed to main_crit function by the methods wrappers.}
+}
+\value{
+Updated results of global optim method
+}
+\description{
+Post-treat results of global optim methods
+}
diff --git a/man/summary_global_optim.Rd b/man/summary_global_optim.Rd
new file mode 100644
index 00000000..6c3f65a8
--- /dev/null
+++ b/man/summary_global_optim.Rd
@@ -0,0 +1,52 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/global_optim_function.R
+\name{summary_global_optim}
+\alias{summary_global_optim}
+\title{Summarizes results of global optim methods}
+\usage{
+summary_global_optim(optim_options, param_info, optim_results, out_dir)
+}
+\arguments{
+\item{optim_options}{List of options of the parameter estimation method, containing:
+\itemize{
+\item \code{ranseed} Set random seed so that each execution of estim_param give the same
+results when using the same seed. If you want randomization, set it to NULL,
+otherwise set it to a number of your choice (e.g. 1234) (optional, default to NULL, which means random seed)
+\item specific options depending on the method used. Click on the links to see examples with the \href{https://sticsrpacks.github.io/CroptimizR/articles/Parameter_estimation_simple_case.html}{simplex}
+and \href{https://sticsrpacks.github.io/CroptimizR/articles/Parameter_estimation_DREAM.html}{DreamZS} methods.
+\item \code{out_dir} \ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#deprecated}{\figure{lifecycle-deprecated.svg}{options: alt='[Deprecated]'}}}{\strong{[Deprecated]}} Definition of \code{out_dir} in \code{optim_options} is no longer supported, use the new argument \code{out_dir} of \code{estim_param} instead.
+}}
+
+\item{param_info}{Information on the parameters to estimate.
+Either
+a list containing:
+\itemize{
+\item \code{ub} and \code{lb}, named vectors of upper and lower bounds (-Inf and Inf can be used if init_values is provided),
+\item \code{default}, named vectors of default values (optional, corresponding parameters are set to these values when the parameter is part of the \code{candidate_param} list and when it is not estimated ; these values are also used as first initial values when the parameters are estimated)
+\item \code{init_values}, a data.frame containing initial
+values to test for the parameters (optional, if not provided, or if less values
+than number of repetitions of the minimization are provided, the, or part
+of the, initial values will be randomly generated using LHS sampling within
+parameter bounds).
+}
+
+or a named list containing for each parameter:
+\itemize{
+\item \code{sit_list}, list the groups of situations for which the current estimated
+parameter must take different values (see \href{https://sticsrpacks.github.io/CroptimizR/articles/Parameter_estimation_Specific_and_Varietal.html}{here}
+for an example),
+\item \code{ub} and \code{lb}, vectors of upper and lower bounds (one value per group),
+\item \code{init_values}, the list of initial values per group  (data.frame, one column per group, optional).
+\item \code{default}, vector of default values per group (optional, the parameter is set to its default value when it is part of the \code{candidate_param} list and when it is not estimated ; the default value is also used as first initial value when the parameter is estimated)
+}}
+
+\item{optim_results}{Results list returned by global optim method wrappers}
+
+\item{out_dir}{Path to the directory where the optimization results will be written. (optional, default to \code{getwd()})}
+}
+\value{
+Prints results of global optim methods
+}
+\description{
+Summarizes results of global optim methods
+}
diff --git a/man/wrap_DEoptim.Rd b/man/wrap_DEoptim.Rd
new file mode 100644
index 00000000..34e31468
--- /dev/null
+++ b/man/wrap_DEoptim.Rd
@@ -0,0 +1,33 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/wrap_DE.R
+\name{wrap_DEoptim}
+\alias{wrap_DEoptim}
+\title{A wrapper for DEoptim package}
+\usage{
+wrap_DEoptim(optim_options, param_info, crit_options)
+}
+\arguments{
+\item{optim_options}{see description in estim_param}
+
+\item{param_info}{see description in estim_param}
+
+\item{crit_options}{List containing several arguments given to \code{estim_param}
+function: \code{param_names}, \code{obs_list}, \code{crit_function}, \code{model_function},
+\code{model_options}, \code{param_info}, \code{transform_obs}, \code{transform_sim}, \code{out_dir}
+that must be passed to main_crit function by the methods wrappers.}
+}
+\value{
+prints, graphs and a list containing:
+\code{final_values}, the vector of estimated values for optimized parameters
+for the repetition that lead to the lowest value of the criterion
+\code{init_values}, the vector of initial values for optimized parameters
+\code{est_values}, data frame of values of all population
+\code{ind_min_crit}, the index of the repetition that lead to the lowest value
+of the criterion
+\code{DE}, the data structure returned by DE
+\code{crit_values} the vector of criterion values for the all population
+}
+\description{
+A wrapper for DEoptim package
+}
+\keyword{internal}
diff --git a/man/wrap_cmaes.Rd b/man/wrap_cmaes.Rd
new file mode 100644
index 00000000..6e8cb8a2
--- /dev/null
+++ b/man/wrap_cmaes.Rd
@@ -0,0 +1,31 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/wrap_cmaes.R, R/wrap_parma.R
+\name{wrap_cmaes}
+\alias{wrap_cmaes}
+\title{Wrap CMA-ES using parma}
+\usage{
+wrap_cmaes(optim_options, param_info, crit_options)
+
+wrap_cmaes(optim_options, param_info, crit_options)
+}
+\arguments{
+\item{optim_options}{see description in estim_param}
+
+\item{param_info}{see description in estim_param}
+
+\item{crit_options}{List containing several arguments given to \code{estim_param}
+function: \code{param_names}, \code{obs_list}, \code{crit_function}, \code{model_function},
+\code{model_options}, \code{param_info}, \code{transform_obs}, \code{transform_sim}, \code{out_dir}
+that must be passed to main_crit function by the methods wrappers.}
+}
+\value{
+list with final_values, init_values, est_values, GE, min_crit_value
+
+list with final_values, init_values, est_values, GE, min_crit_value
+}
+\description{
+Wrap CMA-ES using parma
+
+A wrapper for parma::cmaes package (CMA-ES)
+}
+\keyword{internal}
diff --git a/man/wrap_graDiEnt.Rd b/man/wrap_graDiEnt.Rd
new file mode 100644
index 00000000..b0e3224f
--- /dev/null
+++ b/man/wrap_graDiEnt.Rd
@@ -0,0 +1,25 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/wrap_graDiEnt.R
+\name{wrap_graDiEnt}
+\alias{wrap_graDiEnt}
+\title{A wrapper for graDiEnt package}
+\usage{
+wrap_graDiEnt(optim_options, param_info, crit_options)
+}
+\arguments{
+\item{optim_options}{see description in estim_param}
+
+\item{param_info}{see description in estim_param}
+
+\item{crit_options}{List containing several arguments given to \code{estim_param}
+function: \code{param_names}, \code{obs_list}, \code{crit_function}, \code{model_function},
+\code{model_options}, \code{param_info}, \code{transform_obs}, \code{transform_sim}, \code{out_dir}
+that must be passed to main_crit function by the methods wrappers.}
+}
+\value{
+list with final_values, init_values, est_values, GE, min_crit_value
+}
+\description{
+A wrapper for graDiEnt package
+}
+\keyword{internal}
diff --git a/man/wrap_rgenoud.Rd b/man/wrap_rgenoud.Rd
new file mode 100644
index 00000000..46329331
--- /dev/null
+++ b/man/wrap_rgenoud.Rd
@@ -0,0 +1,37 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/wrap_rgenoud.R
+\name{wrap_rgenoud}
+\alias{wrap_rgenoud}
+\title{A wrapper for rgenoud package (single criterion)}
+\usage{
+wrap_rgenoud(optim_options, param_info, crit_options)
+}
+\arguments{
+\item{optim_options}{see description in estim_param}
+
+\item{param_info}{see description in estim_param}
+
+\item{crit_options}{List containing several arguments given to \code{estim_param}
+function: \code{param_names}, \code{obs_list}, \code{crit_function}, \code{model_function},
+\code{model_options}, \code{param_info}, \code{transform_obs}, \code{transform_sim}, \code{out_dir}
+that must be passed to main_crit function by the methods wrappers.}
+}
+\value{
+A list containing:
+\code{final_values}: vector of estimated parameter values
+corresponding to the best solution.
+\code{init_values}: vector or matrix of initial parameter values.
+\code{est_values}: matrix of estimated values (final solution).
+\code{min_crit_value}: minimum value of the criterion.
+\code{trace_df}: data.frame containing parameters, criterion,
+iteration, evaluation and individual indices .
+\code{obs_var_list}: list of observed variables used by the criterion.
+\code{elapsed}: computation time.
+\code{genoud.pro}: an output file called located in the temporary
+directory provided by tempdir, containing solutions in each
+generation (criterion and parameter values )
+}
+\description{
+A wrapper for rgenoud package (single criterion)
+}
+\keyword{internal}
diff --git a/tests/testthat/test_estim_param_DEoptim.R b/tests/testthat/test_estim_param_DEoptim.R
new file mode 100644
index 00000000..a653cdf3
--- /dev/null
+++ b/tests/testthat/test_estim_param_DEoptim.R
@@ -0,0 +1,178 @@
+context("Test the estim_param function using a toy model")
+library(tibble)
+
+# Define a toy model and its wrapper
+toymodel <- function(nb_days, rB = 0.1, Bmin = 0.1, Bmax = 10, h = 0.5, Bini = 0.01) {
+  # Simulate biomass and yield with a simple logistic model
+  #
+  # Arguments
+  #   - rB: growth rate
+  #   - Bmin: minimum biomass
+  #   - Bmax: maximum biomass
+  #   - h: harvest rate
+  #   - Bini: initial biomass
+  #
+  # Value
+  #   A list containing the simulated biomass and yield
+  #
+  biom <- numeric(nb_days)
+  biom[1] <- Bini
+  for (i in 2:nb_days) {
+    biom[i] <- biom[i - 1] + rB * (biom[i - 1] + Bmin) * (1 - biom[i - 1] / Bmax)
+  }
+  yield <- h * biom
+  return(list(biom = biom, yield = yield))
+}
+
+toymodel_wrapper <- function(param_values = NULL, situation,
+                             model_options, ...) {
+  # A wrapper for toymodel
+  #
+  # Arguments
+  #   - param_values: (optional) a named vector or a tibble containing the
+  #     values of the toymodel parameters to force.
+  #   - situation: Vector of situations names for which results must be
+  #     returned.
+  #     In this case, the names of the situations are coded as "year_suffix"
+  #   - model_options: a tibble containing the begin and end date of simulation
+  #     for each situation.
+  #   - ...: mandatory since CroptimizR will give additional arguments not used
+  #     here
+  #
+  # Value:
+  #   A named list of tibble per situation.
+  #   Each tibble contains columns:
+  #      - Date (POSIXct dates of simulated results),
+  #      - One column per simulated variable (biomass and yield)
+  #
+  results <- list(
+    sim_list = setNames(
+      vector("list", length(situation)),
+      nm = situation
+    ),
+    error = FALSE
+  )
+  attr(results$sim_list, "class") <- "cropr_simulation"
+
+  # Remove dummy in param_values, just here for specific tests
+  if (!is.null(param_values) && "dummy" %in% names(param_values)) {
+    param_values <- param_values[!names(param_values) %in% "dummy"]
+  }
+
+  for (sit in situation) {
+    # Retrieve begin and end date from situation name
+    begin_date <- dplyr::filter(model_options, situation == sit) %>%
+      dplyr::select(begin_date)
+    end_date <- dplyr::filter(model_options, situation == sit) %>%
+      dplyr::select(end_date)
+
+    date_sequence <- seq(from = begin_date[[1]], to = end_date[[1]], by = "day")
+
+    if (!all(names(param_values) %in% c(
+      "rB", "Bmin", "Bmax", "h", "Bini"
+    ))) {
+      warning(
+        paste(
+          "Unknown parameters in param_values:",
+          paste(names(param_values), collapse = ",")
+        )
+      )
+      results$error <- TRUE
+      return(results)
+    }
+
+    # Call the toymodel with varying arguments depending on what is given in
+    # param_values
+    res <- do.call(
+      "toymodel",
+      c(nb_days = length(date_sequence), as.list(param_values))
+    )
+
+    # Fill the result variable
+    results$sim_list[[sit]] <-
+      dplyr::tibble(
+        Date = date_sequence,
+        biomass = res$biom,
+        yield = res$yield
+      )
+  }
+  return(results)
+}
+
+# To make these function accessible to the test environment
+.GlobalEnv$toymodel <- toymodel
+.GlobalEnv$toymodel_wrapper <- toymodel_wrapper
+
+# Use setup() to define shared data for all tests in this file
+setup({
+  # These variables will be available in all tests
+  .GlobalEnv$model_options <- tibble::tibble(
+    situation = c("sit1_2000", "sit1_2001", "sit2_2003", "sit2_2004"),
+    begin_date = as.Date(c("2000-05-01", "2001-05-12", "2003-05-05", "2004-05-15")),
+    end_date = as.Date(c("2000-11-05", "2001-11-20", "2003-11-15", "2004-11-18"))
+  )
+  .GlobalEnv$param_true_values <- c(rB = 0.08, h = 0.55, Bmax = 7)
+
+  # Generate synthetic observations
+  tmp <- toymodel_wrapper(
+    situation = c("sit1_2000", "sit1_2001", "sit2_2003", "sit2_2004"),
+    param_values = param_true_values,
+    model_options = model_options
+  )
+
+  .GlobalEnv$obs_synth <- lapply(tmp$sim_list, function(x) {
+    return(dplyr::filter(x, dplyr::row_number() %% 10 == 0))
+  })
+})
+
+# For cleaning up after tests
+teardown({
+  if (exists("model_options", envir = .GlobalEnv)) rm(model_options, envir = .GlobalEnv)
+  if (exists("param_true_values", envir = .GlobalEnv)) rm(param_true_values, envir = .GlobalEnv)
+  if (exists("obs_synth", envir = .GlobalEnv)) rm(obs_synth, envir = .GlobalEnv)
+})
+
+
+# Test estim_param with DEoptim
+test_that("estim_param 1 step OLS criterion", {
+  param_info <- list(
+    rB = list(lb = 0, ub = 1, init_values = 0.2),
+    h = list(lb = 0, ub = 1, init_values = 0.2),
+    Bmax = list(lb = 5, ub = 15, init_values = 6)
+  )
+
+  optim_options <- list(
+    ranseed = 1234,
+    itermax = 50,
+    reltol = 3.929623e-04,
+    steptol = 4,
+    NP = 30
+  )
+
+  res <- estim_param(
+    obs_list = obs_synth,
+    model_function = toymodel_wrapper,
+    model_options = model_options,
+    optim_method = "DEoptim",
+    crit_function = crit_ols,
+    optim_options = optim_options,
+    param_info = param_info,
+    obs_var = c("biomass", "yield"),
+    situation = c("sit1_2000", "sit1_2001", "sit2_2003"),
+    out_dir = tempdir()
+  )
+  str(res$trace_df)
+
+  expect_equal(res$final_values[["rB"]],
+    param_true_values[["rB"]],
+    tolerance = param_true_values[["rB"]] * 1e-2
+  )
+  expect_equal(res$final_values[["h"]],
+    param_true_values[["h"]],
+    tolerance = param_true_values[["h"]] * 1e-2
+  )
+  expect_equal(res$final_values[["Bmax"]],
+    param_true_values[["Bmax"]],
+    tolerance = param_true_values[["Bmax"]] * 1e-2
+  )
+})
diff --git a/tests/testthat/test_estim_param_GEoptim.R b/tests/testthat/test_estim_param_GEoptim.R
new file mode 100644
index 00000000..b0db6540
--- /dev/null
+++ b/tests/testthat/test_estim_param_GEoptim.R
@@ -0,0 +1,180 @@
+context("Test the estim_param function using a toy model")
+library(tibble)
+
+# Define a toy model and its wrapper
+toymodel <- function(nb_days, rB = 0.1, Bmin = 0.1, Bmax = 10, h = 0.5, Bini = 0.01) {
+  # Simulate biomass and yield with a simple logistic model
+  #
+  # Arguments
+  #   - rB: growth rate
+  #   - Bmin: minimum biomass
+  #   - Bmax: maximum biomass
+  #   - h: harvest rate
+  #   - Bini: initial biomass
+  #
+  # Value
+  #   A list containing the simulated biomass and yield
+  #
+  biom <- numeric(nb_days)
+  biom[1] <- Bini
+  for (i in 2:nb_days) {
+    biom[i] <- biom[i - 1] + rB * (biom[i - 1] + Bmin) * (1 - biom[i - 1] / Bmax)
+  }
+  yield <- h * biom
+  return(list(biom = biom, yield = yield))
+}
+
+toymodel_wrapper <- function(param_values = NULL, situation,
+                             model_options, ...) {
+  # A wrapper for toymodel
+  #
+  # Arguments
+  #   - param_values: (optional) a named vector or a tibble containing the
+  #     values of the toymodel parameters to force.
+  #   - situation: Vector of situations names for which results must be
+  #     returned.
+  #     In this case, the names of the situations are coded as "year_suffix"
+  #   - model_options: a tibble containing the begin and end date of simulation
+  #     for each situation.
+  #   - ...: mandatory since CroptimizR will give additional arguments not used
+  #     here
+  #
+  # Value:
+  #   A named list of tibble per situation.
+  #   Each tibble contains columns:
+  #      - Date (POSIXct dates of simulated results),
+  #      - One column per simulated variable (biomass and yield)
+  #
+  results <- list(
+    sim_list = setNames(
+      vector("list", length(situation)),
+      nm = situation
+    ),
+    error = FALSE
+  )
+  attr(results$sim_list, "class") <- "cropr_simulation"
+
+  # Remove dummy in param_values, just here for specific tests
+  if (!is.null(param_values) && "dummy" %in% names(param_values)) {
+    param_values <- param_values[!names(param_values) %in% "dummy"]
+  }
+
+  for (sit in situation) {
+    # Retrieve begin and end date from situation name
+    begin_date <- dplyr::filter(model_options, situation == sit) %>%
+      dplyr::select(begin_date)
+    end_date <- dplyr::filter(model_options, situation == sit) %>%
+      dplyr::select(end_date)
+
+    date_sequence <- seq(from = begin_date[[1]], to = end_date[[1]], by = "day")
+
+    if (!all(names(param_values) %in% c(
+      "rB", "Bmin", "Bmax", "h", "Bini"
+    ))) {
+      warning(
+        paste(
+          "Unknown parameters in param_values:",
+          paste(names(param_values), collapse = ",")
+        )
+      )
+      results$error <- TRUE
+      return(results)
+    }
+
+    # Call the toymodel with varying arguments depending on what is given in
+    # param_values
+    res <- do.call(
+      "toymodel",
+      c(nb_days = length(date_sequence), as.list(param_values))
+    )
+
+    # Fill the result variable
+    results$sim_list[[sit]] <-
+      dplyr::tibble(
+        Date = date_sequence,
+        biomass = res$biom,
+        yield = res$yield
+      )
+  }
+  return(results)
+}
+
+# To make these function accessible to the test environment
+.GlobalEnv$toymodel <- toymodel
+.GlobalEnv$toymodel_wrapper <- toymodel_wrapper
+
+# Use setup() to define shared data for all tests in this file
+setup({
+  # These variables will be available in all tests
+  .GlobalEnv$model_options <- tibble::tibble(
+    situation = c("sit1_2000", "sit1_2001", "sit2_2003", "sit2_2004"),
+    begin_date = as.Date(c("2000-05-01", "2001-05-12", "2003-05-05", "2004-05-15")),
+    end_date = as.Date(c("2000-11-05", "2001-11-20", "2003-11-15", "2004-11-18"))
+  )
+  .GlobalEnv$param_true_values <- c(rB = 0.08, h = 0.55, Bmax = 7)
+
+  # Generate synthetic observations
+  tmp <- toymodel_wrapper(
+    situation = c("sit1_2000", "sit1_2001", "sit2_2003", "sit2_2004"),
+    param_values = param_true_values,
+    model_options = model_options
+  )
+
+  .GlobalEnv$obs_synth <- lapply(tmp$sim_list, function(x) {
+    return(dplyr::filter(x, dplyr::row_number() %% 10 == 0))
+  })
+})
+
+# For cleaning up after tests
+teardown({
+  if (exists("model_options", envir = .GlobalEnv)) rm(model_options, envir = .GlobalEnv)
+  if (exists("param_true_values", envir = .GlobalEnv)) rm(param_true_values, envir = .GlobalEnv)
+  if (exists("obs_synth", envir = .GlobalEnv)) rm(obs_synth, envir = .GlobalEnv)
+})
+
+
+
+# Test estim_param with DEoptim
+test_that("estim_param 1 step OLS criterion", {
+  param_info <- list(
+    rB = list(lb = 0, ub = 1),
+    h = list(lb = 0, ub = 1),
+    Bmax = list(lb = 5, ub = 15)
+  )
+
+  optim_options <- list(
+    ranseed = 1234,
+    n_params = length(param_info),
+    n_iter = 50,
+    n_particles = 30,
+    n_diff = 2,
+    stop_check = 3,
+    return_trace = TRUE
+  )
+
+  res <- estim_param(
+    obs_list = obs_synth,
+    model_function = toymodel_wrapper,
+    model_options = model_options,
+    optim_method = "graDiEnt",
+    crit_function = crit_ols,
+    optim_options = optim_options,
+    param_info = param_info,
+    obs_var = c("biomass", "yield"),
+    situation = c("sit1_2000", "sit1_2001", "sit2_2003"),
+    out_dir = tempdir()
+  )
+
+  expect_equal(res$final_values[["rB"]],
+    param_true_values[["rB"]],
+    tolerance = param_true_values[["rB"]] * 1e-2
+  )
+  expect_equal(res$final_values[["h"]],
+    param_true_values[["h"]],
+    tolerance = param_true_values[["h"]] * 1e-2
+  )
+  expect_equal(res$final_values[["Bmax"]],
+    param_true_values[["Bmax"]],
+    tolerance = param_true_values[["Bmax"]] * 1e-2
+  )
+})
diff --git a/tests/testthat/test_estim_param_cmaes.R b/tests/testthat/test_estim_param_cmaes.R
new file mode 100644
index 00000000..d8c560a4
--- /dev/null
+++ b/tests/testthat/test_estim_param_cmaes.R
@@ -0,0 +1,179 @@
+context("Test the estim_param function using a toy model")
+library(tibble)
+
+# Define a toy model and its wrapper
+toymodel <- function(nb_days, rB = 0.1, Bmin = 0.1, Bmax = 10, h = 0.5, Bini = 0.01) {
+  # Simulate biomass and yield with a simple logistic model
+  #
+  # Arguments
+  #   - rB: growth rate
+  #   - Bmin: minimum biomass
+  #   - Bmax: maximum biomass
+  #   - h: harvest rate
+  #   - Bini: initial biomass
+  #
+  # Value
+  #   A list containing the simulated biomass and yield
+  #
+  biom <- numeric(nb_days)
+  biom[1] <- Bini
+  for (i in 2:nb_days) {
+    biom[i] <- biom[i - 1] + rB * (biom[i - 1] + Bmin) * (1 - biom[i - 1] / Bmax)
+  }
+  yield <- h * biom
+  return(list(biom = biom, yield = yield))
+}
+
+toymodel_wrapper <- function(param_values = NULL, situation,
+                             model_options, ...) {
+  # A wrapper for toymodel
+  #
+  # Arguments
+  #   - param_values: (optional) a named vector or a tibble containing the
+  #     values of the toymodel parameters to force.
+  #   - situation: Vector of situations names for which results must be
+  #     returned.
+  #     In this case, the names of the situations are coded as "year_suffix"
+  #   - model_options: a tibble containing the begin and end date of simulation
+  #     for each situation.
+  #   - ...: mandatory since CroptimizR will give additional arguments not used
+  #     here
+  #
+  # Value:
+  #   A named list of tibble per situation.
+  #   Each tibble contains columns:
+  #      - Date (POSIXct dates of simulated results),
+  #      - One column per simulated variable (biomass and yield)
+  #
+  results <- list(
+    sim_list = setNames(
+      vector("list", length(situation)),
+      nm = situation
+    ),
+    error = FALSE
+  )
+  attr(results$sim_list, "class") <- "cropr_simulation"
+
+  # Remove dummy in param_values, just here for specific tests
+  if (!is.null(param_values) && "dummy" %in% names(param_values)) {
+    param_values <- param_values[!names(param_values) %in% "dummy"]
+  }
+
+  for (sit in situation) {
+    # Retrieve begin and end date from situation name
+    begin_date <- dplyr::filter(model_options, situation == sit) %>%
+      dplyr::select(begin_date)
+    end_date <- dplyr::filter(model_options, situation == sit) %>%
+      dplyr::select(end_date)
+
+    date_sequence <- seq(from = begin_date[[1]], to = end_date[[1]], by = "day")
+
+    if (!all(names(param_values) %in% c(
+      "rB", "Bmin", "Bmax", "h", "Bini"
+    ))) {
+      warning(
+        paste(
+          "Unknown parameters in param_values:",
+          paste(names(param_values), collapse = ",")
+        )
+      )
+      results$error <- TRUE
+      return(results)
+    }
+
+    # Call the toymodel with varying arguments depending on what is given in
+    # param_values
+    res <- do.call(
+      "toymodel",
+      c(nb_days = length(date_sequence), as.list(param_values))
+    )
+
+    # Fill the result variable
+    results$sim_list[[sit]] <-
+      dplyr::tibble(
+        Date = date_sequence,
+        biomass = res$biom,
+        yield = res$yield
+      )
+  }
+  return(results)
+}
+
+# To make these function accessible to the test environment
+.GlobalEnv$toymodel <- toymodel
+.GlobalEnv$toymodel_wrapper <- toymodel_wrapper
+
+# Use setup() to define shared data for all tests in this file
+setup({
+  # These variables will be available in all tests
+  .GlobalEnv$model_options <- tibble::tibble(
+    situation = c("sit1_2000", "sit1_2001", "sit2_2003", "sit2_2004"),
+    begin_date = as.Date(c("2000-05-01", "2001-05-12", "2003-05-05", "2004-05-15")),
+    end_date = as.Date(c("2000-11-05", "2001-11-20", "2003-11-15", "2004-11-18"))
+  )
+  .GlobalEnv$param_true_values <- c(rB = 0.08, h = 0.55, Bmax = 7)
+
+  # Generate synthetic observations
+  tmp <- toymodel_wrapper(
+    situation = c("sit1_2000", "sit1_2001", "sit2_2003", "sit2_2004"),
+    param_values = param_true_values,
+    model_options = model_options
+  )
+
+  .GlobalEnv$obs_synth <- lapply(tmp$sim_list, function(x) {
+    return(dplyr::filter(x, dplyr::row_number() %% 10 == 0))
+  })
+})
+
+# For cleaning up after tests
+teardown({
+  if (exists("model_options", envir = .GlobalEnv)) rm(model_options, envir = .GlobalEnv)
+  if (exists("param_true_values", envir = .GlobalEnv)) rm(param_true_values, envir = .GlobalEnv)
+  if (exists("obs_synth", envir = .GlobalEnv)) rm(obs_synth, envir = .GlobalEnv)
+})
+
+
+# Test estim_param with CMA-ES (parma)
+test_that("estim_param 1 step OLS criterion", {
+  param_info <- list(
+    rB = list(lb = 0, ub = 1),
+    h = list(lb = 0, ub = 1),
+    Bmax = list(lb = 5, ub = 15)
+  )
+
+  optim_options <- list(
+    ranseed = 1234,
+    return_trace = TRUE,
+    insigma = 0.2,
+    ctrl = list(
+      options = list(MaxIter = 200, PopSize = 30, Trace = 0),
+      CMA = list(active = 1)
+    )
+  )
+
+  res <- estim_param(
+    obs_list = obs_synth,
+    model_function = toymodel_wrapper,
+    model_options = model_options,
+    optim_method = "cmaes",
+    crit_function = crit_ols,
+    optim_options = optim_options,
+    param_info = param_info,
+    obs_var = c("biomass", "yield"),
+    situation = c("sit1_2000", "sit1_2001", "sit2_2003"),
+    out_dir = tempdir()
+  )
+
+  expect_equal(res$final_values[["rB"]],
+    param_true_values[["rB"]],
+    tolerance = param_true_values[["rB"]] * 1e-2
+  )
+  expect_equal(res$final_values[["h"]],
+    param_true_values[["h"]],
+    tolerance = param_true_values[["h"]] * 1e-2
+  )
+  expect_equal(res$final_values[["Bmax"]],
+    param_true_values[["Bmax"]],
+    tolerance = param_true_values[["Bmax"]] * 1e-2
+  )
+})
diff --git a/tests/testthat/test_estim_param_rgenoud.R b/tests/testthat/test_estim_param_rgenoud.R
new file mode 100644
index 00000000..86359997
--- /dev/null
+++ b/tests/testthat/test_estim_param_rgenoud.R
@@ -0,0 +1,194 @@
+context("Test the estim_param function using a toy model")
+library(tibble)
+
+# Define a toy model and its wrapper
+toymodel <- function(nb_days, rB = 0.1, Bmin = 0.1, Bmax = 10, h = 0.5, Bini = 0.01) {
+  # Simulate biomass and yield with a simple logistic model
+  #
+  # Arguments
+  #   - rB: growth rate
+  #   - Bmin: minimum biomass
+  #   - Bmax: maximum biomass
+  #   - h: harvest rate
+  #   - Bini: initial biomass
+  #
+  # Value
+  #   A list containing the simulated biomass and yield
+  #
+  biom <- numeric(nb_days)
+  biom[1] <- Bini
+  for (i in 2:nb_days) {
+    biom[i] <- biom[i - 1] + rB * (biom[i - 1] + Bmin) * (1 - biom[i - 1] / Bmax)
+  }
+  yield <- h * biom
+  return(list(biom = biom, yield = yield))
+}
+clip <- function(x, lb, ub) pmin(pmax(x, lb), ub)
+
+toymodel_wrapper <- function(param_values = NULL, situation,
+                             model_options, ...) {
+  # A wrapper for toymodel
+  #
+  # Arguments
+  #   - param_values: (optional) a named vector or a tibble containing the
+  #     values of the toymodel parameters to force.
+  #   - situation: Vector of situations names for which results must be
+  #     returned.
+  #     In this case, the names of the situations are coded as "year_suffix"
+  #   - model_options: a tibble containing the begin and end date of simulation
+  #     for each situation.
+  #   - ...: mandatory since CroptimizR will give additional arguments not used
+  #     here
+  #
+  # Value:
+  #   A named list of tibble per situation.
+  #   Each tibble contains columns:
+  #      - Date (POSIXct dates of simulated results),
+  #      - One column per simulated variable (biomass and yield)
+  #
+  results <- list(
+    sim_list = setNames(
+      vector("list", length(situation)),
+      nm = situation
+    ),
+    error = FALSE
+  )
+  attr(results$sim_list, "class") <- "cropr_simulation"
+
+  # Remove dummy in param_values, just here for specific tests
+  if (!is.null(param_values) && "dummy" %in% names(param_values)) {
+    param_values <- param_values[!names(param_values) %in% "dummy"]
+  }
+
+  for (sit in situation) {
+    # Retrieve begin and end date from situation name
+    begin_date <- dplyr::filter(model_options, situation == sit) %>%
+      dplyr::select(begin_date)
+    end_date <- dplyr::filter(model_options, situation == sit) %>%
+      dplyr::select(end_date)
+
+    date_sequence <- seq(from = begin_date[[1]], to = end_date[[1]], by = "day")
+
+    if (!all(names(param_values) %in% c(
+      "rB", "Bmin", "Bmax", "h", "Bini"
+    ))) {
+      warning(
+        paste(
+          "Unknown parameters in param_values:",
+          paste(names(param_values), collapse = ",")
+        )
+      )
+      results$error <- TRUE
+      return(results)
+    }
+
+    # --- force parameters into bounds (important for Nelder-Mead) ---
+    if (!is.null(param_values)) {
+      if ("rB" %in% names(param_values))   param_values["rB"]   <- clip(param_values["rB"],   0, 1)
+      if ("h" %in% names(param_values))    param_values["h"]    <- clip(param_values["h"],    0, 1)
+      if ("Bmax" %in% names(param_values)) param_values["Bmax"] <- clip(param_values["Bmax"], 5, 15)
+    }
+    # ---------------------------------------------------------------
+
+    res <- do.call(
+      "toymodel",
+      c(nb_days = length(date_sequence), as.list(param_values))
+    )
+
+    # Fill the result variable
+    results$sim_list[[sit]] <-
+      dplyr::tibble(
+        Date = date_sequence,
+        biomass = res$biom,
+        yield = res$yield
+      )
+  }
+  return(results)
+}
+
+# To make these function accessible to the test environment
+.GlobalEnv$toymodel <- toymodel
+.GlobalEnv$toymodel_wrapper <- toymodel_wrapper
+
+# Use setup() to define shared data for all tests in this file
+setup({
+  # These variables will be available in all tests
+  .GlobalEnv$model_options <- tibble::tibble(
+    situation = c("sit1_2000", "sit1_2001", "sit2_2003", "sit2_2004"),
+    begin_date = as.Date(c("2000-05-01", "2001-05-12", "2003-05-05", "2004-05-15")),
+    end_date = as.Date(c("2000-11-05", "2001-11-20", "2003-11-15", "2004-11-18"))
+  )
+  .GlobalEnv$param_true_values <- c(rB = 0.08, h = 0.55, Bmax = 7)
+
+  # Generate synthetic observations
+  tmp <- toymodel_wrapper(
+    situation = c("sit1_2000", "sit1_2001", "sit2_2003", "sit2_2004"),
+    param_values = param_true_values,
+    model_options = model_options
+  )
+
+  .GlobalEnv$obs_synth <- lapply(tmp$sim_list, function(x) {
+    return(dplyr::filter(x, dplyr::row_number() %% 10 == 0))
+  })
+})
+
+# For cleaning up after tests
+teardown({
+  if (exists("model_options", envir = .GlobalEnv)) rm(model_options, envir = .GlobalEnv)
+  if (exists("param_true_values", envir = .GlobalEnv)) rm(param_true_values, envir = .GlobalEnv)
+  if (exists("obs_synth", envir = .GlobalEnv)) rm(obs_synth, envir = .GlobalEnv)
+})
+
+
+# Test estim_param with DEoptim
+test_that("estim_param 1 step OLS criterion", {
+  param_info <- list(
+    rB = list(lb = 0, ub = 1, init_values = 0.2),
+    h = list(lb = 0, ub = 1, init_values = 0.2),
+    Bmax = list(lb = 5, ub = 15, init_values = 6)
+  )
+
+
+  optim_options <- list(
+    ranseed = 1234,
+    pop.size = 30,
+    max.generations = 20,
+    wait.generations = 30,
+    BFGS = TRUE,
+    optim.method = "Nelder-Mead",
+    boundary.enforcement = 1,
+    gradient.check = FALSE,
+    control = list(maxit = 1, reltol = 1e-06),
+    solution.tolerance = 1e-3,
+    hard.generation.limit = FALSE,
+    print.level = 3,
+    P9 = 0
+  )
+
+  res <- estim_param(
+    obs_list = obs_synth,
+    model_function = toymodel_wrapper,
+    model_options = model_options,
+    optim_method = "genoud",
+    crit_function = crit_ols,
+    optim_options = optim_options,
+    param_info = param_info,
+    obs_var = c("biomass", "yield"),
+    situation = c("sit1_2000", "sit1_2001", "sit2_2003"),
+    out_dir = tempdir()
+  )
+  str(res$trace_df)
+
+  expect_equal(res$final_values[["rB"]],
+               param_true_values[["rB"]],
+               tolerance = param_true_values[["rB"]] * 1e-2
+  )
+  expect_equal(res$final_values[["h"]],
+               param_true_values[["h"]],
+               tolerance = param_true_values[["h"]] * 1e-2
+  )
+  expect_equal(res$final_values[["Bmax"]],
+               param_true_values[["Bmax"]],
+               tolerance = param_true_values[["Bmax"]] * 1e-2
+  )
+})