diff --git a/DESCRIPTION b/DESCRIPTION
index 432e108c82b68c478888467dc237dabcbc63d4dc..30551c58c7ed12d117d79b6a45150a8bf4c6e922 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -13,3 +13,9 @@ Encoding: UTF-8
 LazyData: true
 Roxygen: list(markdown = TRUE)
 RoxygenNote: 7.1.1
+Imports: 
+    dplyr,
+    jsonlite,
+    lubridate,
+    TSstudio
+VignetteBuilder: knitr
diff --git a/R/app_sys.R b/R/app_sys.R
new file mode 100644
index 0000000000000000000000000000000000000000..e19b22b703b589a56aa660bd2a2d9add869a631d
--- /dev/null
+++ b/R/app_sys.R
@@ -0,0 +1,9 @@
+#' Access files in the current app
+#'
+#' @param ... Character vector specifying directory and or file to
+#'     point to inside the current package.
+#'
+#' @noRd
+app_sys <- function(...){
+  system.file(..., package = "rvgest")
+}
diff --git a/R/get_objectives.R b/R/get_objectives.R
new file mode 100644
index 0000000000000000000000000000000000000000..cb2f1a9b784f3157f72a13e17deff78336b7971a
--- /dev/null
+++ b/R/get_objectives.R
@@ -0,0 +1,60 @@
+#' Load objective data for \code{\link{vgest_run}}
+#'
+#' @param objective_thresholds File location of threshold table
+#' @param objective_stations File location of objective station table
+#' @param lakes File location of lake table
+#'
+#' @return A dataframe containing one line by objective and the following columns:
+#' 
+#' - station  (character): identifier of the objective station
+#' - flood (boolean): TRUE for high flow mitigation objective, and FALSE for low flow support
+#' - level (character): code composed with "l" for low-flow and "h" for high flow followed by the number of the level 
+#' - threshold (numeric): value of the threshold in m3/s
+#' - lakes (dataframe): Dataframe containing lake details (name, min storage, max storage)
+#' 
+#' @export
+#'
+#' @examples
+#' # Get objectives stored in IRMaRA package
+#' df <- get_objectives()
+#' # Objectives at Paris
+#' df[df[,"station"] == "PARIS_05",]
+#' # Lake details concerning the 40th objective in the table
+#' df[57, "lakes"]
+get_objectives <- function(
+  objective_thresholds = app_sys("seine", "objective_thresholds.txt"),
+  objective_stations = app_sys("seine", "objective_stations.json"), 
+  lakes = app_sys("seine", "lakes.txt")
+) {
+  thresholds <- read.delim(objective_thresholds)
+  stations <- jsonlite::fromJSON(
+    readLines(objective_stations)
+  )
+  dfLakes <- read.delim(lakes)
+  row.names(dfLakes) <- dfLakes$name
+  bFirst = TRUE
+  for(iStation in 1:nrow(thresholds)) {
+    station <- thresholds[iStation, 1]
+    for(type in c("l", "h")) {
+      bFlood = (type == "h")
+      for(level in grep(paste0(type, ".*"), names(thresholds))) {
+        threshold <- thresholds[iStation, level]
+        df1 <- data.frame(
+          station = station,
+          flood = bFlood,
+          level = names(thresholds)[level],
+          threshold = threshold
+        )
+        df1$lakes <- list(dfLakes[stations[[station]]$lakes,])
+        if(bFirst) {
+          df <- df1
+          bFirst <- FALSE
+        } else {
+          df <- rbind(df, df1)
+        }
+      }
+    }
+  }
+  class(df) <- c("Objectives", class(df))
+  return(df)
+}
\ No newline at end of file
diff --git a/R/get_vobj_ts.R b/R/get_vobj_ts.R
new file mode 100644
index 0000000000000000000000000000000000000000..2da97414b405e8b1ab9cb6b8c0859dd349d72a12
--- /dev/null
+++ b/R/get_vobj_ts.R
@@ -0,0 +1,72 @@
+#' Method to transform VOBJ data into time series
+#'
+#' Transform VOBJ data read by [vgest_read], [vgest_read_one] or [vgest_read_all].
+#'
+#' @param x
+#'
+#' @return
+#' @rdname get_vobj_ts
+#' @export
+#'
+#' @examples
+#' \dontrun{
+#' dfTs <- get_vobj_ts(vgest_read_all(get_objectives()[1,], "./database"))
+#' }
+get_vobj_ts <- function(x) {
+  UseMethod("get_vobj_ts", x)
+}
+
+#' Transform a VOBJ data into time series
+#'
+#' @param x a dataframe provided by [vgest_read_one] or [vgest_read], or an item of the list provided by [vgest_read_all]
+#'
+#' @return A dataframe with one column `date` and one column `storage` containing the time series of the lake storage in m3.
+#' @rdname get_vobj_ts
+#' @export
+#'
+#' @examples
+#' \dontrun{
+#' # Get the first lake of the first objective
+#' dfTs <- get_vobj_ts(vgest_read_all(get_objectives()[1,], "./database")[[1]])
+#' # Can also be done by
+#' dfTs <- get_vobj_ts(vgest_read_one(1, get_objectives()[1,], "./database"))
+#' }
+get_vobj_ts.Vobj <- function(x) {
+  # Build the date vector
+  # Get year list
+  vYears <- as.numeric(
+    gsub("Y", "", names(x)[grep(pattern = "Y[0-9]{4}", names(x))])
+  )
+  # Combine dates dd/mm with years
+  vDates <- as.vector(sapply(vYears, function(x) {paste(x[,"DAY"], x, sep = "/")}))
+  # Convert dates to POSIX
+  vDates <- lubridate::dmy(vDates)
+  ts <- unlist(x[,grep(pattern = "Y[0-9]{4}", names(x))])
+  data.frame(date=vDates, storage=ts)
+}
+
+#' Transform list of VOBJ data into time series
+#'
+#' Use for binding results of several lakes for one objective at one station.
+#'
+#' @param x list of VOBJ provided by [vgest_read_all]
+#'
+#' @return Dataframe with the time series with a column `date` and one column by lake with the storage of the lake in m3.
+#' @rdname get_vobj_ts
+#' @export
+#'
+#' @examples
+#' \dontrun{
+#' dfTs <- get_vobj_ts(vgest_read_all(get_objectives()[1,], "./database"))
+#' }
+get_vobj_ts.ListVobj <-function(x) {
+  lVobjTs <- lapply(x, get_vobj_ts.Vobj)
+  df <- lVobjTs[[1]]
+  if(length(lVobjTs) > 1) {
+    for(i in 2:length(lVobjTs)) {
+      df <- cbind(df, lVobjTs[[i]]$storage)
+    }
+  }
+  names(df) <- c("date", names(x))
+  return(df)
+}
\ No newline at end of file
diff --git a/R/plot_isofrequency.R b/R/plot_isofrequency.R
new file mode 100644
index 0000000000000000000000000000000000000000..8b4d744a7c31241aee147046f1bbc2405cbf4a5b
--- /dev/null
+++ b/R/plot_isofrequency.R
@@ -0,0 +1,65 @@
+#' Plot one lake isofrequency annual curves for one objective at one station
+#'
+#' @param x [list] containing one [data.frame] per lake
+#' @param freq [numeric] vector of frequencies to plot
+#' @param result.dir [character] The path to the result database (default: "database")
+#'
+#' @return [NULL]
+#' @export
+#'
+plot_isofrequency <- function(x, freq, result.dir = "database") {
+  lObj <- vgest_read_all(x, result.dir)
+  nLakes <- nrow(x$lakes[[1]])
+  if(nLakes > 1) {
+    par(mfrow=c(ceiling(nLakes/2),2))
+  }
+  lapply(1:nLakes, function(i) {plot_isofrequency_lake(lObj[[i]], freq, x$lakes[[1]][i,], ceiling(i/2))})
+  sLowHigh <- c("low", "high")
+  mtext(
+    paste("Objective", x$threshold, "m3/s", sLowHigh[x$flood + 1], "flows threshold at", x$station),
+    side = 3, col = "black", line = -2, cex = 1.2, font = 2, outer = TRUE
+  )
+}
+
+#' Plot one lake isofrequency annual curves for one objective at one station
+#'
+#' This function is mainly called by [plot_isofrequency()].
+#'
+#' @param vObj dataframe produced by [vgest_read_one()] and stored in a list by [vgest_read_all()]
+#' @param frequencies vector of frequencies to plot
+#' @param lake lake data extract from column `lakes` of objective data given by [get_objectives()]
+#' @param top.margin top margin applied on the plot for the title
+#'
+#' @return [NULL]
+#' @export
+#'
+plot_isofrequency_lake <- function(vObj, frequencies, lake, top.margin) {
+  frequencies <- paste0("F", format(frequencies, digits = 5, nsmall = 5))
+  vObj[,-1] <- (vObj[,-1]) / 1E6 + lake$min
+  # Contrasted color palette adapted from https://sashat.me/2017/01/11/list-of-20-simple-distinct-colors/
+  sColors <- c("#003A80", "#f58231", "#3cb44b", "#e6194B", "#911eb4", "#ffe119", "#f032e6", "#9A6324", "#000075", "#808000", "#42d4f4", "#a9a9a9", "#bfef45", "#469990")
+  par(mgp=c(2,0.5,0), mar = c(2.5,3.5,5.5 - 2.2*top.margin,0.5))
+  plot(
+    rep(lake$min, 365),
+    type = "l", lwd = 1.5, col = "grey", lty = 2,
+    xaxt="n",
+    main = NULL,
+    ylab = "Reservoir Storage (mcm)",
+    xlab = NULL,
+    ylim = c(0, lake$max),
+    cex.lab = 0.6, cex.axis = 0.6
+  )
+  lines(rep(lake$max, 365), lwd = 1.5, col = "grey", lty = 2)
+  i <- 0
+  for(frequency in frequencies) {
+    i <- i + 1
+    lines(vObj[,frequency], lwd = 2, col = sColors[i])
+  }
+  mtext(
+    paste(lake$name, "reservoir (Active cap.", lake$max - lake$min, "mcm)"),
+    side = 3, font = 2, cex = 0.75
+  )
+  xtick <- c(0, 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31)
+  xtick <- sapply(1:13, function(x) {sum(xtick[1:x])})
+  axis(side=1, at=xtick, labels = unlist(strsplit("JFMAMJJASONDJ", NULL)), cex.axis = 0.6)
+}
diff --git a/R/save_isofrequency.R b/R/save_isofrequency.R
new file mode 100644
index 0000000000000000000000000000000000000000..d405e3f861d1bc4fb2867143a19e3bd12f250efa
--- /dev/null
+++ b/R/save_isofrequency.R
@@ -0,0 +1,15 @@
+save_isofrequency <- function(x, freq, result.dir = "database", output.dir = "output") {
+  dir.create(output.dir, showWarnings = FALSE)
+  file <- file.path(
+    output.dir, paste0(
+      paste("isofreq", x$station, x$level, x$threshold, sep = "_"),
+      ".png")
+  )
+  cat("Plotting", file, "...")
+  png(file = file,
+    width = 16, height = 10, units = "cm", res = 200
+  )
+  plot_isofrequency(x, freq, result.dir)
+  dev.off()
+  cat(" OK\n")
+}
\ No newline at end of file
diff --git a/R/vgest_cost.R b/R/vgest_cost.R
new file mode 100644
index 0000000000000000000000000000000000000000..4b1cc299e1cec439fb808f98031b7022ff5c352f
--- /dev/null
+++ b/R/vgest_cost.R
@@ -0,0 +1,98 @@
+#' Calculate the cost function for one simulation (One objective at one station)
+#'
+#' Uses ouput of PaChrono.txt or VObj\[1-4\].dat for the calculation.
+#'
+#' @param data Data source, see details
+#' @param objective one row of the dataframe given by [get_objectives()]
+#'
+#' @return the total cost for one objective at one station in m3/day
+#' @rdname vgest_cost
+#' @export
+#'
+vgest_cost <- function(data, objective) {
+  UseMethod("vgest_cost", data)
+}
+
+
+#' Calculate the total cost for all the lakes
+#'
+#' For each lake, it's the mean daily storage for a low-flow support objective and the mean daily available storage capacity for a high flow mitigation objective.
+#'
+#' @param Vobj A [matrix] or a [data.frame] with one column by lake, in the same order as `objective$lakes`
+#' @param objective one row of the dataframe given by [get_objectives()]
+#'
+#' @return total cost of all the lakes in m3/day
+#'
+vgest_cost_lakes <- function(Vobj, objective) {
+  lakes <- objective$lakes[[1]]
+  sum(
+    sapply(
+      1:nrow(lakes), function(i) {
+        V <- Vobj[,i]
+        if(objective$flood) {
+          V <- (lakes$max[i] - lakes$min[i]) * 1E6 - V
+        }
+        mean(V, na.rm = TRUE)
+      }
+    )
+  )
+}
+
+
+#' @param objective one row of the dataframe given by [get_objectives()]
+#'
+#' @rdname vgest_cost
+#' @export
+#'
+#' @examples
+#' \dontrun{
+#' # This should be done after the execution of vgest for the concerned objective
+#' objective <- get_objectives()[1,]
+#' lResultVobj <- vgest_read_all(objective)
+#' vgest_cost(lResultVobj, objective)
+#' }
+vgest_cost.ListVobj <- function(data, objective) {
+  lVobj <- lapply(data, function(x) {unlist(x[,grep(pattern = "Y[0-9]{4}", names(x))])})
+  Vobj <- matrix(unlist(lVobj), ncol = length(lVobj))
+  vgest_cost_lakes(Vobj, objective)
+}
+
+
+#' @param data the [data.frame] given by one item of the [list] provided by [vgest_read_chrono()]
+#'
+#' @return the total cost for one objective at one station in m3/day
+#' @rdname vgest_cost
+#' @export
+#'
+#' @examples
+#' \dontrun{
+#' # This should be done after the execution of vgest for the concerned objective
+#' objective <- get_objectives()[1,]
+#' lChrono <- vgest_read_chrono(objective)
+#' vgest_cost(lChrono[[1]], objective)
+#' }
+vgest_cost.Chrono <- function(data, objective) {
+  vgest_cost_lakes(as.data.frame(data[,grep("Vobj\\.[1-9]\\.", names(data))]), objective)
+}
+
+
+#'
+#' @param data the [list] given by [vgest_read_chrono()]
+#' @param objective [data.frame] given by [get_objectives()]
+#'
+#' @return A list with items named \[station\]_\[high/low\]_\[threshold\] containing the total cost for a list of objectives and stations in m3/day
+#' @rdname vgest_cost
+#' @export
+#'
+#' @examples
+#' \dontrun{
+#' # This should be done after the execution of vgest for the concerned objective
+#' objective <- get_objectives()[1,]
+#' lChronos <- vgest_read_chrono(objective, distributionType = 2)
+#' vgest_cost(lChronos, objective)
+#' }
+vgest_cost.ListChrono <- function(data, objective) {
+  lCosts <- lapply(1:length(data), function(i) {vgest_cost.Chrono (data[[i]], objective[i,])})
+  names(lCosts) <- names(data)
+  return(lCosts)
+}
\ No newline at end of file
diff --git a/R/vgest_read.R b/R/vgest_read.R
new file mode 100644
index 0000000000000000000000000000000000000000..17903dfece647705b83a7c52ab12bf802ce33119
--- /dev/null
+++ b/R/vgest_read.R
@@ -0,0 +1,90 @@
+#' Read result file `VOBJi.DAT` for one lake of one objective at one station
+#'
+#' This function is called by [vgest_read_one()], please run this one instead.
+#'
+#' @param file complete path to the `VOBJi.DAT` file
+#' @param bFlood boolean `TRUE` for high flow mitigation objective, `FALSE` for low flow support
+#'
+#' @return dataframe with the content of the `VOBJi.DAT` file
+#' @export
+#'
+#' @examples
+#' \dontrun{
+#' df <- vgest_read("database/PARIS_05/high_800/VOBJ1.dat")
+#' }
+vgest_read <- function(file, bFlood) {
+  # Managing headers
+  s <- readLines(file, n = 11)
+  headers <- unlist(strsplit(s[10], "\\s+"))
+  frequencies <- unlist(strsplit(s[11], "\\s+"))[-1]
+  if(bFlood) {
+    frequencies[-1] <- rev(frequencies[-1])
+  }
+  headers <- c(
+    "DAY",
+    paste0("Y", headers[2:(which(headers == "RANG-->")-1)]),
+    paste0("R", headers[which(headers == "RANG-->"):(which(headers == "RETOUR-->")-1)]),
+    paste0("F", frequencies)
+  )
+  # Managing column widths of fixed format
+  widths <- c(8, rep(16, length(headers) - 1))
+  # Read file in fixed format
+  df <- read.fwf(
+    file = file,
+    widths = widths,
+    header = FALSE,
+    skip = 11,
+    col.names = headers
+  )
+  # Remove redundancy columns
+  df$`FFREQUENCE...` <- NULL
+  df$`RRANG...` <- NULL
+  class(df) <- c("Vobj", class(df))
+  return(df)
+}
+
+#' Read result file `VOBJi.DAT` for one lake of one objective at one station
+#'
+#' This function is preferred to [vgest_read()] because it builds the path for the file to read.
+#'
+#' @param iLake Lake number for the current station
+#' @param x one line of the dataframe produced by [get_objectives()]
+#' @param result.dir path where results of VGEST runs are stored (See [vgest_run_store()])
+#'
+#' @return dataframe with the content of the `VOBJi.DAT` file
+#' @export
+#'
+#' @examples
+#' \dontrun{
+#' vgest_read_one(1, get_objectives()[1,], "./database")
+#' }
+vgest_read_one <- function(iLake, x, result.dir) {
+  sLowHigh <- c("low", "high")
+  file <- paste0(
+    file.path(
+      result.dir, x$station, paste0(sLowHigh[x$flood + 1], "_", x$threshold), "VOBJ"
+    ),
+    iLake,
+    ".dat"
+  )
+  vgest_read(file, x$flood)
+}
+
+#' Read all result files `VOBJi.DAT` for all the lakes of one objective at one station
+#'
+#' @param x one line of the dataframe produced by [get_objectives()]
+#' @param result.dir path where results of VGEST runs are stored (See [vgest_run_store()])
+#'
+#' @return list with dataframes produced by
+#' @export
+#'
+#' @examples
+#' \dontrun{
+#' vgest_read_all(get_objectives()[1,], "./database")
+#' }
+vgest_read_all <- function(x, result.dir = "database") {
+  lObj <- lapply(1:nrow(x$lakes[[1]]), vgest_read_one, x, result.dir)
+  names(lObj) <- x$lakes[[1]]$name
+  class(lObj) <- c("ListVobj", class(lObj))
+  return(lObj)
+}
\ No newline at end of file
diff --git a/R/vgest_read_chrono.R b/R/vgest_read_chrono.R
new file mode 100644
index 0000000000000000000000000000000000000000..462a682c7f696dc9cfe7bd479be30fc6dc0a6145
--- /dev/null
+++ b/R/vgest_read_chrono.R
@@ -0,0 +1,118 @@
+#' Read Chrono.txt or PaChrono.txt
+#'
+#' @description
+#' Read the time series provided by VGEST for the forward calculation (Chrono.txt) and backward calculation (PaChrono.txt)
+#'
+#' @details
+#' The format of the file is has follow. Headers are in line 53 followed by the complete time series which is backward in time for PaChrono.txt. The columns and respected widths are:
+#' - Downstream date (10)
+#' - QXsous (22)
+#'
+#' For each lake, we next have:
+#' - Upstream date (14)
+#' - 25 columns of flow and storage data (22)
+#' - if `DistributionType %in% c(4,5)`, 1 extra column (22)
+#' - 5 columns of codes (6)
+#'
+#' And then at the end:
+#' - 7 columns of flow (22)
+#' - 3 columns of codes (8)
+#' - 1 extra column of code for `DistributionType==5` (8)
+#'
+#' The file is saved in RDS format for quicker reading the next time.
+#'
+#' @param x file to read with txt extension
+#' @param ... further arguments passed to or from other methods
+#'
+#' @return
+#' @export
+#' @rdname vgest_read_chrono
+#'
+vgest_read_chrono <- function(x, ...) {
+  UseMethod("vgest_read_chrono", x)
+}
+
+
+#'
+#' @param nLakes number of lakes in the file
+#' @param distributionType Distribution type. See [vgest_write_batch] details
+#' @return [data.frame] with the content of the file
+#' @export
+#' @rdname vgest_read_chrono
+#'
+vgest_read_chrono.default <- function(x, nLakes, distributionType, ...) {
+  # Reading cached file if exists
+  rdsFile <- paste0(sub('\\..[^\\.]*$', '', x), ".rds")
+  if(file.exists(rdsFile)) {
+    df <- readRDS(rdsFile)
+  } else {
+    # Widths
+    width1Lake <- c(14, rep(22, 25))
+    if(distributionType %in% c(4,5)) { width1Lake <- c(width1Lake, 22)}
+    width1Lake <- c(width1Lake, rep(6,5))
+    widths <- c(10, 22, rep(width1Lake, nLakes), rep(22,7), rep(8,3))
+    if(distributionType == 5) { widths <- c(widths, 8) }
+    # Headers
+    if(distributionType %in% c(4,5)) {
+      headLine <- 54
+    } else {
+      headLine <- 53
+    }
+
+
+    s <- readLines(x, n = headLine)[headLine]
+    headers = trimws(unlist(read.fwf(textConnection(s), widths = widths, as.is = TRUE), use.names = FALSE))
+    # Read file in fixed format
+    df <- read.fwf(
+      file = x,
+      col.names = headers,
+      widths = widths,
+      skip = headLine,
+      strip.white = TRUE
+    )
+    df[,grep(pattern = "date", names(df))] <- lapply(df[,grep(pattern = "date", names(df))], as.Date, format = "%d/%m/%Y")
+    df <- dplyr::arrange(df, date.aval)
+    saveRDS(df, file = rdsFile)
+  }
+  class(df) <- c("Chrono", class(df))
+  return(df)
+}
+
+#'
+#' @param result.dir path for storing the result of vgest run. The result is stored in a subfolder named high or low (depending on \code{bFlood}) followed by the threshold
+#' @param distributionType Distribution type. See [vgest_write_batch] details
+#' @param backward boolean `TRUE` for reading "PaChrono.txt", `FALSE` for reading "Chrono.txt"
+#'
+#' @return A list with items named \[station\]_\[high/low\]_\[threshold\] containing [data.frame] with the content of each file
+#' @export
+#' @rdname vgest_read_chrono
+#'
+#' @examples
+#' \dontrun{
+#' objective <- get_objectives()[1,]
+#' distributionType <- 2
+#' vgest_run_store(objective,
+#'                 1, 1, "Q_NAT_1900-2009.txt",
+#'                 "01/01/1900", "31/12/2009", distributionType)
+#' lChrono <- vgest_read_chrono(objective, distributionType)
+#' }
+vgest_read_chrono.Objectives <- function(x, distributionType, result.dir = "database", backward = TRUE, ...) {
+  if(backward) {
+    fileName <- "PaChrono.txt"
+  } else {
+    fileName <- "Chrono.txt"
+  }
+  sLowHigh <- c("low", "high")
+  objectivePaths <- file.path(result.dir, x$station, paste0(sLowHigh[x$flood + 1], "_", x$threshold))
+  l <- lapply(1:nrow(x), function(i) {
+    cat("Treating folder", objectivePaths[i], "...\n")
+    vgest_read_chrono.default(
+      file.path(objectivePaths[i], fileName),
+      nLakes = nrow(x[i, "lakes"][[1]]),
+      distributionType = distributionType
+    )
+  })
+  names(l) <- paste(x$station, sLowHigh[x$flood + 1], x$threshold, sep = "_")
+  class(l) <- c("ListChrono", class(l))
+  return(l)
+}
diff --git a/R/vgest_read_qnat.R b/R/vgest_read_qnat.R
new file mode 100644
index 0000000000000000000000000000000000000000..097e760f16cba560cc6f92589513b5a41e358fb1
--- /dev/null
+++ b/R/vgest_read_qnat.R
@@ -0,0 +1,16 @@
+#' Read daily flow time series used by VGEST
+#' 
+#' @details The file should be an ASCII file with tab separators, the first column contains the date in format `YYYYMMDD`.
+#' The returned [data.frame] has the class attributes "Qnat" for further use with S3 methods.
+#'
+#' @param file Path to the file to read
+#'
+#' @return a [data.frame] with the content of the file
+#' @export
+#'
+vgest_read_qnat <- function(file) {
+  Qnat <- read.delim(file)
+  Qnat$Dates <- as.POSIXct(as.character(Qnat$Dates), format = "%Y%m%d", tz = "UTC")
+  class(Qnat) <- c("Qnat", class(Qnat))
+  return(Qnat)
+}
\ No newline at end of file
diff --git a/R/vgest_run.R b/R/vgest_run.R
new file mode 100644
index 0000000000000000000000000000000000000000..d7617cb5085d1b929230bcf58201281b16c84acb
--- /dev/null
+++ b/R/vgest_run.R
@@ -0,0 +1,23 @@
+#' Run VGEST and stop execution if an error is encountered during execution
+#'
+#' @details
+#' Delete the content of the `RESULTAT` folder before running
+#'
+#' @param vgest_location Location of VGEST installation (Default "../vgest")
+#'
+#' @return Nil
+#' @export
+#'
+#' @examples
+#' \dontrun{
+#' # Run vgest with the current configuration
+#' vgest_run()
+#' }
+vgest_run <- function(vgest_location = "../vgest") {
+  workingDir <- getwd()
+  setwd(vgest_location)
+  file.remove(list.files(path = "RESULTAT", full.names = TRUE))
+  iError <- system(command = "VGEST.exe")
+  setwd(workingDir)
+  if(iError != 0) stop("Running VGEST has returned an error", iError)
+}
\ No newline at end of file
diff --git a/R/vgest_run_store.R b/R/vgest_run_store.R
new file mode 100644
index 0000000000000000000000000000000000000000..99a4d6900c1cff78e5d8c12a436330b48306ff56
--- /dev/null
+++ b/R/vgest_run_store.R
@@ -0,0 +1,86 @@
+#' Prepare, run and store results of one or several VGEST instances
+#'
+#' @param ... other parameter passed to the method [vgest_run_store.character]
+#'
+#' @return [NULL]
+#' @export
+#' @rdname vgest_run_store
+#'
+vgest_run_store <- function(x, ...) {
+  UseMethod("vgest_run_store", x)
+}
+
+
+#'
+#' @inheritParams vgest_write_batch
+#' @param x See argument station
+#' @param result.dir path for storing the result of vgest run. The result is stored in a subfolder named high or low (depending on \code{bFlood}) followed by the threshold
+#'
+#' @export
+#' @rdname vgest_run_store
+#'
+vgest_run_store.character <- function(x, reservoirRuleSet, networkSet,
+                            Qfile, startDate, endDate,
+                            bFlood, threshold,
+                            distributionType, distributionOption = NULL,
+                            vgest_location = "../vgest",
+                            objective_file = "BATCH",
+                            formatResult = 1,
+                            result.dir = "database", ...) {
+  sLowHigh <- c("low", "high")
+  station <- x
+  cat("Run VGEST for configuration: ", station, Qfile, sLowHigh[bFlood + 1], threshold, "...")
+  vgest_write_batch(reservoirRuleSet, networkSet, Qfile, startDate, endDate,
+                    station, bFlood, threshold, distributionType, distributionOption,
+                    vgest_location, objective_file, formatResult)
+  vgest_run(vgest_location)
+  sTo <- file.path(result.dir, station, paste0(sLowHigh[bFlood + 1], "_", threshold))
+  cat(" Store files to", sTo)
+  dir.create(sTo, showWarnings = FALSE, recursive = TRUE)
+  file.remove(list.files(path = sTo, full.names = TRUE))
+  if(!all(file.copy(
+    list.files(file.path(vgest_location, "RESULTAT"), full.names = TRUE),
+    sTo, overwrite = TRUE
+  ))) stop("Error copying Vojb files from VGEST to Irmara")
+  cat(" - OK\n")
+}
+
+#' @param x row(s) of a [data.frame] provided by [get_objectives()]
+#'
+#' @export
+#' @rdname vgest_run_store
+#'
+#' @examples
+#' \dontrun{
+#' # Example with `vgest_run_store.Objectives`
+#' # Running vgest for:
+#' # - the first objective returned by `get_objectives()`
+#' # - the first configuration of reservoir rules
+#' # - the first configuration of network
+#' # - the naturalized hydrological flows of the file located in DONNEES/Q_NAT_1900-2009.txt
+#' # - doing the optimization on the period between 01/01/1900 and 31/12/2009
+#' # - a task distribution function of present volumes and maximum usable volume replenishment times from the start of time steps
+#' vgest_run_store(get_objectives()[1,],
+#'                 1, 1, "Q_NAT_1900-2009.txt",
+#'                 "01/01/1900", "31/12/2009", 2)
+#'
+#' # Example with `vgest_run_store.character`
+#' # Running vgest for:
+#' # - the first configuration of reservoir rules
+#' # - the first configuration of network
+#' # - the naturalised hydrological flows of the file located in DONNEES/Q_NAT_1900-2009.txt
+#' # - doing the optimisation on the period between 01/01/1900 and 31/12/2009
+#' # - an objective located at "PARIS_05"
+#' # - for a flood threshold of 800 m3/s
+#' # - a fixed task distribution between reservoirs with a relative change of -20%, -10%, +50%, -30%
+#' vgest_run_store("PARIS_05", 1, 1, "Q_NAT_1900-2009.txt", "01/01/1900", "31/12/2009",
+#'                   TRUE, 800, 1, c(-0.2, -0.1, 0.5, -0.3))
+#' }
+vgest_run_store.Objectives <- function(x, reservoirRuleSet, networkSet,
+                          Qfile, startDate, endDate, ...) {
+  nothing <- apply(x, 1, function(y) {
+    vgest_run_store(y$station, reservoirRuleSet, networkSet,
+                            Qfile, startDate, endDate,
+                            y$flood, y$threshold, ...)
+  })
+}
\ No newline at end of file
diff --git a/R/vgest_write_batch.R b/R/vgest_write_batch.R
new file mode 100644
index 0000000000000000000000000000000000000000..a39a2eeb09182df66aa23f2b0d04a60e7b367d15
--- /dev/null
+++ b/R/vgest_write_batch.R
@@ -0,0 +1,85 @@
+#' Write configuration files for running VGEST
+#'
+#' Write CHOIX.TXT, objective file and eventually REGLAGES.TXT in PARAMETR folder of VGEST.
+#'
+#' @details
+#' The format of the PARAMETR/CHOIX.TXT file is as follow:
+#'
+#' - line 1: name of the station where the flow target is located
+#' - line 2: rank of the set of parameters describing the reservoirs and their management rules (constraints and local instructions)
+#' - line 3: rank of the set of parameters describing the network
+#' - line 4: name of the flow data file (hydrological regime without reservoirs)
+#' - line 5: name of the file describing the annual objective hydrograph at the bottom station of the system
+#' - line 6: objective type: 0 for low flow support, 1 for high flow lamination
+#' - line 7: start of calculation period
+#' - line 8: end of calculation period
+#' - line 9: code indicating which sub programs should be run:
+#'   - 0: All sub programs
+#'   - 1: Only backward simulation
+#' - line 10: code for output format of volume results: 1 for absolute values in m3 , 2 for relative values with respect to Vtot or sum(Vtot)
+#' - line 11: code indicating how the flow to be stored is distributed between the reservoirs:
+#'   - 1=fixed;
+#'   - 2=function of present volumes and maximum usable volume replenishment times from the start of time steps;
+#'   - 3=aimed at balancing the filling rates at the end of the time step;
+#'   - 4=aiming to balance at the end of the time step the times Tpot of reservoir evolution towards extreme state (see line 11) with average inputs;
+#'   - 5=aiming to balance at the end of the time step the times Tpot of evolution of the reservoir towards extreme state (see line 11) with given quantity of the contributions of each quantity.
+#' - line 12 (used only if 4 or 5 on line 11): code indicating the nature of Tpot :
+#'   - 1: Tpot is the minimum potential duration Tpot1 of reconstitution of the maximum usable volume (obtaining V=Vtot or V=0, depending on the nature of the objective (support or rolling) and the direction of the calculations)
+#'   - 2: Tpot is the minimum potential duration Tpot2 of exhaustion of the usable volume (obtaining V=Vtot or V=0, depending on the nature of the objective (support or rolling) and the direction of the calculations)
+#'
+#' @param reservoirRuleSet rank of the set of parameters describing the reservoirs and their management rules (constraints and local instructions)
+#' @param networkSet rank of the set of parameters describing the network
+#' @param station name of the station where the flow target is located
+#' @param Qfile name of the flow data file (hydrological regime without reservoirs)
+#' @param bFlood boolean describing the objective type: \code{FALSE} for low flow support, \code{TRUE} for high flow lamination
+#' @param startDate start of calculation period in "dd/mm/yyyy" format
+#' @param endDate end of calculation period in "dd/mm/yyyy" format
+#' @param threshold value of the threshold in m3/s
+#' @param distributionType see in details the description of line 10 of CHOIX.TXT
+#' @param distributionOption for \code{distributionType=1}, this should contains the relative change in comparison with a repartition based on storage capacity of fixed repartition between reservoirs. If \code{distributionType=4} or \code{5}, see in details the description of line 11 of CHOIX.TXT (default \code{NULL})
+#' @param vgest_location location of the vgest software (default "../vgest")
+#' @param objective_file name of the file used for storing the threshold hydrograph
+#' @param formatResult code for output format of volume results: 1 for absolute values in m3 , 2 for relative values with respect to Vtot or sum(Vtot)
+#' @param subPrograms [integer] 0 for running all sub programs, 1 for running only backward simulation (default value)
+#'
+#' @return
+#' @export
+#'
+#' @examples
+#' \dontrun{
+#' # Preparing the run of vgest for:
+#' # - the first configuration of reservoir rules
+#' # - the first configuration of network
+#' # - the naturalised hydrological flows of the file located in DONNEES/Q_NAT_1900-2009.txt
+#' # - doing the optimisation on the period between 01/01/1900 and 31/12/2009
+#' # - an objective located at "PARIS_05"
+#' # - for a flood threshold of 800 m3/s
+#' # - a fixed task distribution between reservoirs with a relative change of -20%, -10%, +50%, -30%
+#' vgest_write_batch(1, 1, "Q_NAT_1900-2009.txt", "01/01/1900", "31/12/2009",
+#'                   "PARIS_05", TRUE, 800, 1, c(-0.2, -0.1, 0.5, -0.3))
+#' }
+vgest_write_batch <- function(reservoirRuleSet, networkSet,
+                              Qfile, startDate, endDate,
+                              station, bFlood, threshold,
+                              distributionType, distributionOption = NULL,
+                              vgest_location = "../vgest",
+                              objective_file = "BATCH",
+                              formatResult = 1,
+                              subPrograms = 2) {
+
+  if(distributionType %in% c(1, 4, 5) & is.null(distributionOption)) {
+    stop("distribution should be defined for distributionType = 1, 4 or 5")
+  }
+  sMode <- as.numeric(bFlood) # 0 for drought, 1 for flood
+  s <- c(station, reservoirRuleSet, networkSet,
+         Qfile, objective_file, sMode, startDate, endDate,
+         subPrograms, formatResult, distributionType)
+  if(distributionType == 1) {
+    writeLines(as.character(distributionOption), file.path(vgest_location, "PARAMETR", "REGLAGE.txt"))
+  } else if(distributionType %in% c(4,5)) {
+    s <- c(s,distributionOption)
+  }
+
+  writeLines(s, file.path(vgest_location, "PARAMETR", "CHOIX.txt"))
+  writeLines(paste("01/01", threshold), file.path(vgest_location, "OBJECTIF", objective_file))
+}
diff --git a/README.md b/README.md
index 3f557413c41d8912b66a87a769f3836d61f77598..6c82cdd86c06e43e4a55f48a032f9faef5378576 100644
--- a/README.md
+++ b/README.md
@@ -1,3 +1,52 @@
 # RVGEST
 
-R package wrapper for VGEST
\ No newline at end of file
+R package wrapper for VGEST
+
+## Installation
+
+### Requirements
+
+The third party program 'VGEST' is required to perform the optimisations. Clone or download 'VGEST' sources and install it at the same folder level as the rvgest folder ('rvgest' and 'vgest' folders should be located in the same folder). Follow the instruction at https://gitlab.irstea.fr/in-wop/vgest for installation and compilation. Then, prepare the model for the Seine Basin case by following the steps of the "get started" example.
+
+Finally, copy a naturalised flow database in the `vgest/DONNEES` folder.
+
+### Local installation steps
+
+Installation of the package from source should be done in 3 steps:
+
+- download sources
+- run `roxygen` for generating `NAMESPACE` file and documentation from sources
+- install the package
+
+### Download sources
+
+First possibility, clone the repository with git (recommended):
+
+If you have configured a SSH key for gitlab (See: https://docs.gitlab.com/ee/gitlab-basics/create-your-ssh-keys.html): `git clone git@gitlab-ssh.irstea.fr:in-wop/rvgest.git`
+
+Otherwise, use HTTPS: `git clone https://gitlab.irstea.fr/in-wop/rvgest.git`
+
+Second possibility, download the source code at https://gitlab.irstea.fr/in-wop/rvgest/-/archive/master/rvgest-master.zip
+
+### Build and install the package
+
+Open the project file `rvgest.Rproj` in RStudio and type:
+
+```
+roxygen2::roxygenise()
+remotes::install_local()
+```
+
+## Usage
+
+### Running VGEST and plot outputs
+
+See the dedicated vignette "run_vgest_basics" in the `vignettes` folder of the package.
+
+### Documentation
+
+See all available functions and their documentation by typing `help(package = "rvgest")` in a R console.
+
+## Code of Conduct
+
+Please note that the rvgest project is released with a [Contributor Code of Conduct](https://contributor-covenant.org/version/2/0/CODE_OF_CONDUCT.html). By contributing to this project, you agree to abide by its terms.
\ No newline at end of file
diff --git a/inst/seine/lakes.txt b/inst/seine/lakes.txt
new file mode 100644
index 0000000000000000000000000000000000000000..d8408ab584b7d8e24a401a213b515c2d5246b71f
--- /dev/null
+++ b/inst/seine/lakes.txt
@@ -0,0 +1,5 @@
+name	min	max
+YONNE	8.5	82.5
+SEINE	6.6	219.5
+AUBE	2.3	183.5
+MARNE	10	364.5
diff --git a/inst/seine/objective_description.txt b/inst/seine/objective_description.txt
new file mode 100644
index 0000000000000000000000000000000000000000..775297b10dc79596ee095f5e336f072b16079b45
--- /dev/null
+++ b/inst/seine/objective_description.txt
@@ -0,0 +1,9 @@
+code	type	name	description
+l1	low	vigilance	sensibilisation for water usage reduction
+l2	low	alert	restriction of 30% of non productive water withdrawals
+l3	low	reinforced alert	restriction of 50% of water withdrawals excluding drinking water supply
+l4	low	crisis	only the drinking water supply and respect for biological life
+are insured
+h1	high	yellow	localized overflow flooding
+h2	high	orange	flood causing major overflows
+h3	high	red	major flood, direct and widespread threat
diff --git a/inst/seine/objective_stations.json b/inst/seine/objective_stations.json
new file mode 100644
index 0000000000000000000000000000000000000000..3e11b82a93c6e793901eef82c33b0233cf94c6ea
--- /dev/null
+++ b/inst/seine/objective_stations.json
@@ -0,0 +1,29 @@
+{
+  "ALFOR_16": {
+    "lakes": ["YONNE",  "SEINE", "AUBE"]
+  },
+  "ARCIS_24": {
+    "lakes": ["AUBE"]
+  },
+  "CHALO_21": {
+    "lakes": ["MARNE"]
+  },
+  "COURL_23": {
+    "lakes": ["YONNE"]
+  },
+  "GURGY_02": {
+    "lakes": ["YONNE"]
+  },
+  "MERY-_22": {
+    "lakes": ["SEINE"]
+  },
+  "NOGEN_13": {
+    "lakes": ["SEINE", "AUBE"]
+  },
+  "NOISI_17": {
+    "lakes": ["MARNE"]
+  },
+  "PARIS_05": {
+    "lakes": ["YONNE",  "SEINE", "AUBE", "MARNE"]
+  }
+}
diff --git a/inst/seine/objective_thresholds.txt b/inst/seine/objective_thresholds.txt
new file mode 100644
index 0000000000000000000000000000000000000000..7bdde573c2b704eb808399000815076890b812d5
--- /dev/null
+++ b/inst/seine/objective_thresholds.txt
@@ -0,0 +1,10 @@
+station	l1	l2	l3	l4	h1	h2	h3
+ARCIS_24	6.3	5	4	3.5	110	260	400
+MERY-_22	7.3	5	4	3.5	140	170	400
+NOGEN_13	25	20	17	16	180	280	420
+GURGY_02	14	12.5	11	9.2	220	340	400
+COURL_23	23	16	13	11	550	700	900
+ALFOR_16	64	48	41	36	850	1200	1400
+CHALO_21	12	11	9	8	330	520	700
+NOISI_17	32	23	20	17	350	500	650
+PARIS_05	81	60	51	45	950	1600	2000
diff --git a/rvgest.Rproj b/rvgest.Rproj
index 69fafd4b6dddad27500cfc67efb9fb16e86a96bd..cba1b6b7a0a292098af0e0b1ce41b7846a59bd5e 100644
--- a/rvgest.Rproj
+++ b/rvgest.Rproj
@@ -14,7 +14,6 @@ LaTeX: pdfLaTeX
 
 AutoAppendNewline: Yes
 StripTrailingWhitespace: Yes
-LineEndingConversion: Posix
 
 BuildType: Package
 PackageUseDevtools: Yes
diff --git a/vignettes/.gitignore b/vignettes/.gitignore
new file mode 100644
index 0000000000000000000000000000000000000000..097b241637da023174b0f2e3715bd0291d9ded37
--- /dev/null
+++ b/vignettes/.gitignore
@@ -0,0 +1,2 @@
+*.html
+*.R
diff --git a/vignettes/run_vgest_basics.Rmd b/vignettes/run_vgest_basics.Rmd
new file mode 100644
index 0000000000000000000000000000000000000000..dc6119b58e9d360aed9825b0387e6837c81f6303
--- /dev/null
+++ b/vignettes/run_vgest_basics.Rmd
@@ -0,0 +1,82 @@
+---
+title: "Running VGEST and plot outputs"
+output: rmarkdown::html_vignette
+vignette: >
+  %\VignetteIndexEntry{Running VGEST and plot outputs}
+  %\VignetteEngine{knitr::rmarkdown}
+  %\VignetteEncoding{UTF-8}
+---
+
+```{r, include = FALSE}
+knitr::opts_chunk$set(
+  collapse = TRUE,
+  comment = "#>"
+)
+```
+
+```{r setup}
+library(irmara)
+```
+
+## Prerequisites
+
+Download the vgest repository or clone it at https://gitlab.irstea.fr/in-wop/vgest. Follow the instruction for compile it and follow the steps of the "get started" example.
+
+Copy a naturalised flow database in the `vgest/DONNEES` folder.
+
+In the following functions called, the location of VGEST is defined in the parameter `vgest_location`. By default, its value is :
+
+```{r}
+vgest_location = "../../vgest"
+```
+
+
+This correspond to an installation of VGEST located at the same folder level as the irmara folder (they should appear together in the same folder).
+
+## Get the list of objectives
+
+An objective is a combination of one downstream location, one type of objective (low/high flow), and a threshold.
+
+These information are provided in the data of the Irmara package.
+
+```{r}
+objectives <- get_objectives()
+objectives
+```
+
+For selecting the first high flow threshold at Paris:
+
+```{r}
+objParisH1 <- objectives[objectives$station == "PARIS_05" & objectives$level == "h1",]
+objParisH1
+```
+The column `lakes` contains a list containing a data.frame detailing information about the lakes upstream the concerned station:
+
+```{r}
+objParisH1$lakes[[1]]
+```
+
+## Running vgest for one or several objectives
+
+`vgest_run_store` is a function that run a configuration of VGEST for one or several objectives and store the result folder of each instance of VGEST in a folder. Type `?vgest_run_store` in the console for the help of the function.
+
+The instruction below run VGEST for:
+
+* The first level of high flow objective at Paris
+* The reservoir rule set number 1 (all local constraints)
+* The network set number 1
+* The naturalized flows provided in `Q_NAT_1900-2009.txt` starting from 1900 up to 2009
+* The type of distribution task number 2 which function of present volumes and maximum usable volume replenishment times from the start of time steps
+
+```{r}
+vgest_run_store(x = objParisH1, 
+                reservoirRuleSet = 1,
+                networkSet = 1, 
+                Qfile = "Q_NAT_1900-2009.txt", 
+                startDate = "01/01/1900",
+                endDate = "31/12/2009", 
+                distributionType = 2, 
+                vgest_location = vgest_location, 
+                result.dir = "../database")
+```
+
diff --git a/vignettes/v01_plotting_time_series.Rmd b/vignettes/v01_plotting_time_series.Rmd
new file mode 100644
index 0000000000000000000000000000000000000000..d641f364ce4067e469e3486004a7e510da5c0ee9
--- /dev/null
+++ b/vignettes/v01_plotting_time_series.Rmd
@@ -0,0 +1,51 @@
+---
+title: "v01_plotting_time_series"
+output: rmarkdown::html_vignette
+vignette: >
+  %\VignetteIndexEntry{v01_plotting_time_series}
+  %\VignetteEngine{knitr::rmarkdown}
+  %\VignetteEncoding{UTF-8}
+---
+
+```{r, include = FALSE}
+knitr::opts_chunk$set(
+  collapse = TRUE,
+  comment = "#>"
+)
+```
+
+```{r setup}
+library(irmara)
+```
+
+Run VGEST for the first objectif high flow mitigation at the station "PARIS_05" (950 m3/s). 
+
+The reservoir rule set corresponds to taking into account the local reserved flow at the intakes of the lakes (and Yonne lake return flow) and a maximum outflow on the Yonne lake not limited by hydropower. "Reference" (high flow) local constraints are not taken into account.
+
+Distribution task used is function of present volumes and maximum usable volume replenishment times from the start of time steps.
+
+```{r}
+vgest_run_store(get_objectives()[61,], 
+                reservoirRuleSet = 5, 
+                networkSet = 1, 
+                Qfile = "Q_NAT_1900-2009.txt", 
+                start = "01/01/1900", 
+                end = "31/12/2009", 
+                distributionType = 2,
+                vgest_location = "../../vgest",
+                result.dir = "../database"
+                )
+```
+
+Get the time series from VOBJ files.
+
+```{r}
+dfTs <- get_v.obj_ts(vgest_read_all(get_objectives()[61,], "../database"))
+```
+
+Plot time series. Use the zoom function to explore it.
+
+```{r}
+TSstudio::ts_plot(dfTs, type = "multiple")
+```
+