R/Grid_Matrix.R

In rte-antares-rpackage/antaresWaterValues: Generating water values for Antares

#' Calculate grid layer matrix of Bellman values and water values
#'
#' @param area An 'antares' area.
#' @param simulation_names Names of simulations to retrieve.
#' @param simulation_values Values for the simulation.
#' @param reward_db a table contains the reward values generated by the function \code{get_Reward}
#' If it's NULL it auto calculated.
#' @param inflow a table contains the hydro storage generated by the function \code{readAntares}
#' with th option \code{hydrostorage = TRUE}. If it's NULL it auto calculated.
#' @param nb_cycle Number of times to run the algorithm.
#' @param district_name Name of the district used to store output.
#' @param mcyears MC years to consider, by default all of them.
#' @param week_53 Water values for week 53, by default 0.
#' @param method Perform mean of grids algorithm or grid of means algorithm or
#'  grid of quantile algorithm.
#' @param states_step_ratio Discretization ratio to generate steps levels
#' between the reservoir capacity and zero . Defaults to 0.05
#' @param q_ratio from 0 to 1. the probability used in quantile method
#' to determine a bellman value which q_ratio all bellman values are equal or
#' less to it. (quantile(q_ratio))
#' @param test_week the week number u want to see it's calculation information
#' in the console
#' @param reservoir_capacity Reservoir capacity for the given area in MWh,
#'  if \code{NULL} (the default), value in Antares is used if available else
#'  a prompt ask the user the value to be used.
#' @param na_rm Remove NAs
#' @param only_input if TRUE skip bellman values calculation and return the input
#' @param until_convergence Boolean. TRUE to repeat cycle until convergence or
#'  attending the limit.
#' @param convergence_rate from 0 to 1. Define the convergence level from which
#'  we suppose that no need to continue another cycle..
#' @param convergence_criteria the value define convergence. if the difference
#' between two water values is less then this value those values are converged.
#' @param cycle_limit Define cycles limit when you are in the until_convergence mod.
#' @param pumping Boolean. True to take into account the pumping.
#' @param efficiency in [0,1]. efficient ratio of pumping.
#' @param opts
#'   List of simulation parameters returned by the function
#'   \code{antaresRead::setSimulationPath}
#' @param shiny Boolean. True to run the script in shiny mode.
#' @param stop_rate the percent from which the calculation stop. for example
#' \code{stop_rate=5} means the calculation stop if there is a week with less then
#' 5\% accessibles states.
#' @param debug_week the number of the week to open the process in debug mode
#' @param correct_concavity Binary argument (default to false). True to correct concavity of Bellman values.
#' @param correct_monotony_gain Binary argument (default to false). True to correct monotony of gains.
#' @param ... further arguments passed to or from other methods.
#' @param penalty_low Penalty for violating the bottom rule curve, comparable to the unsupplied energy cost
#' @param penalty_high Penalty for violating the top rule curve, comparable to the spilled energy cost
#' @param method_old_gain If T, linear interpolation used between simulations reward, else smarter interpolation based on marginal prices
#' @param hours_reward_calculation If method_old_gain=F, vector of hours used to evaluate costs/rewards of pumping/generating
#' @param controls_reward_calculation If method_old_gain=F, vector of controls evaluated
#' @param max_hydro_hourly Hourly maximum pumping and turbining powers
#' @param max_hydro_weekly Weekly maximum pumping and turbining powers
#' @param force_final_level Binary. Whether final level should be constrained
#' @param final_level Final level (in percent between 0 and 100) if final level is constrained but different from initial level
#' @param penalty_final_level_low Penalties for both bottom rule curve to constrain final level
#' @param penalty_final_level_high Penalties for top rule curve to constrain final level
#'
#' @return List of a data.frame with aggregated water values and
#' a data.frame of more detailed water values
#' @export
#'


  Grid_Matrix <- function(area, simulation_names,reward_db=NULL,inflow=NULL,
                             simulation_values = NULL, nb_cycle = 1L,
                             district_name = "water values district", mcyears = NULL,
                             week_53 = 0,
                             states_step_ratio = 0.01,
                             reservoir_capacity = NULL,
                             na_rm = FALSE,
                             method ,
                             only_input=FALSE,
                             q_ratio=0.5,
                             test_week=NULL,
                             opts = antaresRead::simOptions(),
                             shiny=F,
                             until_convergence=F,
                             convergence_rate=0.9,
                             convergence_criteria=1,
                             cycle_limit=10,
                             pumping=F,efficiency=1,stop_rate=0,
                             debug_week=54,
                             correct_concavity = FALSE,
                             correct_monotony_gain = FALSE,
                             penalty_low = 3000,
                             penalty_high = 3000,
                             method_old_gain = F,
                             hours_reward_calculation = c(seq.int(0,168,10),168),
                             controls_reward_calculation = NULL,
                          max_hydro_hourly=NULL,
                          max_hydro_weekly=NULL,
                          force_final_level = F,
                          final_level = NULL,
                          penalty_final_level_low = NULL,
                          penalty_final_level_high = NULL,
                        ...) {



  #----- shiny Loader

  # check the method chosen to calculate Bellman values is a valid method
  methods <- c("mean-grid","grid-mean","quantile")
  if (!method %in% methods){
    stop("Unknown method. available methods: ('mean-grid','grid-mean','quantile') ")
  }

  assertthat::assert_that(class(opts) == "simOptions")

  # Number of weeks
  n_week <- 52


  # Replacing NULL values
  if(is.null(reward_db)){
    if(!is.null(simulation_names)){
      if (is.null(simulation_values)) {
        simulation_values <- constraint_generator(area,length(simulation_names),pumping,opts, mcyears=mcyears)
        message(paste0("Using simulation_values: ", paste(simulation_values, collapse = ", ")))
      }}
    }else
      {simulation_names <- reward_db$simulation_names
      simulation_values <- reward_db$simulation_values
      reward_db <- as.data.table(reward_db$reward)
      }
  if (is.null(mcyears)) {
    mcyears <- opts$parameters$general$nbyears
    mcyears <- seq(0,mcyears)
  }



  # Getting reservoir capacity (ie maximum level)
  {if (is.null(reservoir_capacity)) {
    niveau_max <- get_reservoir_capacity(area = area)
    if (length(niveau_max) == 0) {
      ask_niveau_max <- "Failed to retrieve reservoir capacity from Antares, please specify value (in GWh, e.g. 10000 for France):\n"
      niveau_max <- readline(prompt = ask_niveau_max)
      niveau_max <- as.numeric(niveau_max)*1000
    }
    if (length(niveau_max) == 0)
      stop("Failed to retrieve reservoir capacity, please specify it explicitly with 'reservoir_capacity'.")
  } else {
    niveau_max <- reservoir_capacity
  }}


  # Checking states_step_ratio
  if ((!is.numeric(states_step_ratio))|(states_step_ratio<0))
    stop("Failed Not valid States_step_ratio, please change it between 0 and 1.")

  states_steps <- niveau_max*states_step_ratio


  # States matrix for all weeks plus an other week representing the end of the year
  states <- matrix( rep(seq(from = niveau_max, to = 0, by = -states_steps), n_week + 1), byrow = FALSE, ncol = n_week + 1)


  if(is.null(efficiency))
  {
    efficiency <- getPumpEfficiency(area = area,opts = opts)
  }

  # Getting inflow in the reservoir
  {
    inflow <- get_inflow(area=area, opts=opts,mcyears=mcyears)
    inflow[with(inflow, order(inflow$tsId, inflow$timeId)),]
    inflow <- inflow[, c("area", "tsId" , "timeId", "time", "hydroStorage")]
    # inflow[, timeId := gsub(pattern = "\\d{4}-w", replacement = "", x = time)]
    inflow[, "timeId" := as.numeric(inflow$timeId)]
    inflow <- inflow %>%
      dplyr::group_by(.data$area, .data$timeId, .data$tsId) %>%
      dplyr::mutate(hydroStorage = sum(.data$hydroStorage, na.rm = TRUE)) %>%
      as.data.table()

  }

  options("antares" = opts)

  # Getting maximum pumping and generating powers for each hour
  if (is.null(max_hydro_hourly)){
    max_hydro_hourly <- get_max_hydro(area,timeStep = "hourly")
  }


  # Calculating reward
  {
    if(is.null(reward_db))
    {
      if(is.null(controls_reward_calculation)&method_old_gain==F){
          controls_reward_calculation <- constraint_generator(area=area,nb_disc_stock=10,
                                                            pumping=pumping,
                                                            pumping_efficiency=efficiency,
                                                            opts=opts, mcyears=mcyears)
      }

      # Addig controls used is the simulation to the controls used to interpolate reward
      if (("mcYear" %in% names(simulation_values))&!("mcYear" %in% names(controls_reward_calculation))){
        controls_reward_calculation <- dplyr::cross_join(controls_reward_calculation,
                                                 data.frame(mcYear=mcyears))
      }
      controls_reward_calculation <- rbind(simulation_values,controls_reward_calculation) %>%
        dplyr::select(-c("sim")) %>%
        dplyr::distinct() %>%
        dplyr::arrange(.data$week,.data$u)

      # Checking the minimum number of controls to calculate reward
      if (pumping==T){
        if (method_old_gain==T){
          assertthat::assert_that(length(simulation_names)>=3,
                                  msg="If you have less than 3 simulations, you have to interpolate with marginal prices")
        } else {
          if (!("mcYear" %in% names(controls_reward_calculation))){
            n_distinct <- controls_reward_calculation %>%
              dplyr::group_by(.data$week)
          } else {
            n_distinct <- controls_reward_calculation %>%
              dplyr::group_by(.data$week,.data$mcYear)
          }
          nb_distinct <- n_distinct %>%
            dplyr::summarise(n=dplyr::n()) %>%
            dplyr::pull("n") %>%
            min()
          assertthat::assert_that(nb_distinct>=3,
                                  msg="You should have at least 3 different controls for each week.")
        }
      } else {
        if (method_old_gain==T){
          assertthat::assert_that(length(simulation_names)>=2,
                                  msg="If you have less than 3 simulations, you have to interpolate with marginal prices")
        } else {
          if (!("mcYear" %in% names(controls_reward_calculation))){
            n_distinct <- controls_reward_calculation %>%
              dplyr::group_by(.data$week)
          } else {
            n_distinct <- controls_reward_calculation %>%
              dplyr::group_by(.data$week,.data$mcYear)
          }
          nb_distinct <- n_distinct %>%
            dplyr::summarise(n=dplyr::n()) %>%
            dplyr::pull("n") %>%
            min()
          assertthat::assert_that(nb_distinct>=2,
                                  msg="You should have at least 2 different controls for each week.")
        }
      }

      # Calculate reward
      reward_db <- get_Reward(simulation_names = simulation_names, district_name = district_name,
                           opts = opts, correct_monotony = correct_monotony_gain,
                           method_old = method_old_gain,max_hydro=max_hydro_hourly,
                           hours=hours_reward_calculation,
                           possible_controls=controls_reward_calculation,
                           simulation_values = simulation_values, mcyears=mcyears,area=area,
                           district_balance=district_name, pump_eff = efficiency)

      # Retriving controls (u) for each week
      decision_space <- reward_db$simulation_values
      if ("sim" %in% names(decision_space)){
        decision_space <- dplyr::select(decision_space,-c("sim"))
      }
      decision_space <- round(decision_space)
      reward_db <- reward_db$reward


    } else {
      if ("sim" %in% names(simulation_values)){
        decision_space <- simulation_values %>% dplyr::select(-c("sim"))
      } else {
        decision_space <- simulation_values
      }
      decision_space <- round(decision_space)
    }

    reward_db <- reward_db[reward_db$timeId %in% seq_len(n_week)]}

  # Reservoir (rule curves)
  {
    reservoir <- readReservoirLevels(area, timeStep = "weekly",
                                     byReservoirCapacity = FALSE, opts = opts)
    vars <- c("level_low", "level_avg", "level_high")
    reservoir[,
              (vars) := lapply(.SD, function(x) {round(x * max(states))}),
              .SDcols = vars
    ]
  }

  if (force_final_level){
    final_level <- final_level*max(states)/100
  }

  # Prepare data.table watervalues that give Bellman value for each week, each MC year and each state
  watervalues <- data.table(expand.grid(weeks = seq_len(n_week+1), years = mcyears, statesid=seq_len(nrow(states))))

  # add states values
  {
    statesdt <- as.data.table(states)  #convert states matrix to data.table
    statesdt <- melt(data = statesdt, measure.vars = seq_len(ncol(states)), variable.name = "weeks", value.name = "states")
    statesdt[, "weeks" := as.numeric(gsub("V", "", statesdt$weeks))] #turn weeks to numbers V1==> 1
    statesdt[, "statesid" := seq_along(states), by = c("weeks")] # add id to refer to the state
    statesdt[, "states" := round(statesdt$states)]
  }

  # add states plus 1 (ie states for the following week)
  {
    statesplus1 <- copy(statesdt)
    statesplus1[, "weeks" := statesplus1$weeks - 1]
    statesplus1 <- statesplus1[, list(states_next = list(unlist(states))), by = c("weeks")]
    statesplus1 <- dplyr::left_join(x = statesdt, y = statesplus1, by = c("weeks"))
    watervalues <- dplyr::right_join(x = watervalues, y = statesplus1, by = c("weeks","statesid"))
  }

  # add inflow
  watervalues <- dplyr::left_join(x = watervalues, y = inflow[, list(weeks = inflow$timeId, years = inflow$tsId, hydroStorage=inflow$hydroStorage)], by = c("weeks", "years"))
  #at this point water values is the table containing (weeks,year,states,statesid;states_next,hydroStorage)

  #add reward
  watervalues <- dplyr::nest_join(x = watervalues, y = reward_db, by = c("weeks"="timeId","years"="mcYear"))


  #at this point we added the rewards for each weekly_amount

  # add rule curves
  watervalues <- dplyr::left_join(x = watervalues, y = reservoir[, list(weeks = reservoir$timeId,
                                                                        level_low = reservoir$level_low,
                                                                        level_high = reservoir$level_high)], by = "weeks")
  #here we added the lvl_high and low of the reservoir

  # add empty columns for future results
  watervalues$value_node <- NA_real_
  watervalues$transition <- NA_real_
  watervalues$transition_reward <- NA_real_
  watervalues$next_bellman_value <- NA_real_


  # get maximum generating and pumping energy for each week
  if (is.null(max_hydro_weekly)){
    max_hydro_weekly <- get_max_hydro(area,timeStep = "weekly")
  }
  E_max <-max_hydro_weekly$turb
  P_max <- max_hydro_weekly$pump*efficiency



  # prepare Bellman values for the end of the year (week 53)
  {
  if (length(week_53) == 1) week_53 <- rep_len(week_53, length(states))
  next_week_values <- week_53   # approximation to get initial bellman values from initial water values
  counter <- 0

  }

  # at this state, everything is ready to apply dynamic programming equation

  ####

  if (only_input) return(watervalues)

  if(!until_convergence){ # with a predefined number of cycles
    next_week_values <- rep_len(next_week_values, nrow(watervalues[watervalues$weeks==52]))

    for (n_cycl in seq_len(nb_cycle)) {

      cat("Calculating value nodes, cycle number:", n_cycl, "\n")

      pb <- utils::txtProgressBar(min = 0, max = 51, style = 3)

      watervalues[watervalues$weeks==53,"value_node" :=next_week_values]

      for (i in rev(seq_len(52))) { # rep(52:1, times = nb_cycle)


        temp <- watervalues[watervalues$weeks==i]

        if(debug_week==i)browser()

        # Bellman equation for week i
        temp <- Bellman(Data_week=temp,
                        next_week_values_l = next_week_values,
                        decision_space=dplyr::filter(decision_space,week==i),
                        E_max=E_max[i],
                        P_max=P_max[i],
                        states_steps=states_steps,
                        method=method,
                        mcyears = mcyears,
                        q_ratio= q_ratio,
                        counter = i,
                        niveau_max=niveau_max,
                        stop_rate=stop_rate,
                        penalty_level_low=if((i==52)&force_final_level&(n_cycl==1)){penalty_final_level_low}else{penalty_low},
                        penalty_level_high=if((i==52)&force_final_level&(n_cycl==1)){penalty_final_level_high}else{penalty_high},
                        lvl_high =if((i==52)&force_final_level&(n_cycl==1)){final_level}else{temp$level_high[1]},
                        lvl_low =if((i==52)&force_final_level&(n_cycl==1)){final_level}else{temp$level_low[1]}
                        )



        if(shiny&n_cycl==1&i==52){
          for (p in c("shinybusy")){
            if (!requireNamespace(p, quietly = TRUE)) {
              stop(
                paste0("Packageb", p, " must be installed to use this function."),
                call. = FALSE
              )
            }
          }
          shinybusy::show_modal_spinner(spin = "atom",color = "#0039f5")
        }

        # Correct concacity of Bellman values if needed
        if(correct_concavity){
          temp$value_node <- correct_concavity(temp,i:i)
        }

        # Write results for week i
        watervalues[watervalues$weeks==i,"value_node" :=temp$value_node]
        watervalues[watervalues$weeks==i,"transition" :=temp$transition]
        watervalues[watervalues$weeks==i,"transition_reward" :=temp$transition_reward]
        watervalues[watervalues$weeks==i,"next_bellman_value" :=temp$next_bellman_value]

        # Get Bellman values for week i that correspond to future Bellman values for week i-1
        next_week_values <- temp$value_node

        utils::setTxtProgressBar(pb = pb, value = 52 - i)

      }
      close(pb)
      next_week_values <- dplyr::filter(temp,.data$weeks==1)$value_node
      if(nrow(watervalues[is.na(watervalues$value_node)&(watervalues$weeks<=52)])>=1){
        message("Error in the calculation of Bellman values")
      }
      # Calculate water values by derivating Bellman values and applying penalties on rules curves for the current week
      value_nodes_dt <- build_data_watervalues(watervalues,statesdt,reservoir,penalty_high,penalty_low,
                                               force_final_level=if(n_cycl==1){force_final_level}else{F},penalty_final_level_low,penalty_final_level_high,final_level=final_level)

    }

  }else{ # calculating Bellman values until convergence
    next_week_values <- rep_len(next_week_values, nrow(watervalues[watervalues$weeks==52]))

    for (n_cycl in seq_len(cycle_limit)) {

      cat("Calculating value nodes, cycle number:", n_cycl, "\n")

      pb <- utils::txtProgressBar(min = 0, max = 51, style = 3)

      watervalues[watervalues$weeks==53,"value_node" :=next_week_values]

      for (i in rev(seq_len(52))) { # rep(52:1, times = nb_cycle)


        temp <- watervalues[watervalues$weeks==i]

        temp <- Bellman(Data_week=temp,
                        next_week_values_l = next_week_values,
                        decision_space=dplyr::filter(decision_space,week==i),
                        E_max=E_max[i],
                        P_max=P_max[i],
                        states_steps=states_steps,
                        method=method,
                        mcyears = mcyears,
                        q_ratio= q_ratio,
                        counter = i,
                        niveau_max=niveau_max,
                        stop_rate=stop_rate,
                        penalty_level_low=if((i==52)&force_final_level&(n_cycl==1)){penalty_final_level_low}else{penalty_low},
                        penalty_level_high=if((i==52)&force_final_level&(n_cycl==1)){penalty_final_level_high}else{penalty_high},
                        lvl_high =if((i==52)&force_final_level&(n_cycl==1)){final_level}else{temp$level_high[1]},
                        lvl_low =if((i==52)&force_final_level&(n_cycl==1)){final_level}else{temp$level_low[1]}
        )

        if(shiny&n_cycl==1&i==52){
          for (p in c("shinybusy")){
            if (!requireNamespace(p, quietly = TRUE)) {
              stop(
                paste0("Packageb", p, " must be installed to use this function."),
                call. = FALSE
              )
            }
          }
          shinybusy::show_modal_spinner(spin = "atom",color = "#0039f5")
        }

        if(correct_concavity){
          temp$value_node <- correct_concavity(temp,i:i)
        }

        watervalues[watervalues$weeks==i,"value_node" :=temp$value_node]
        watervalues[watervalues$weeks==i,"transition" :=temp$transition]
        watervalues[watervalues$weeks==i,"transition_reward" :=temp$transition_reward]
        watervalues[watervalues$weeks==i,"next_bellman_value" :=temp$next_bellman_value]



        next_week_values <- temp$value_node

        utils::setTxtProgressBar(pb = pb, value = 52 - i)

      }
      close(pb)
      next_week_values <- temp[temp$weeks==1]$value_node
      if(nrow(watervalues[is.na(watervalues$value_node)&(watervalues$weeks<=52)])>=1){
        message("Error in the calculation of Bellman values")
      }
      value_nodes_dt <- build_data_watervalues(watervalues,statesdt,reservoir,penalty_high,penalty_low,
                                               force_final_level=if(n_cycl==1){force_final_level}else{F},
                                               penalty_final_level_low,penalty_final_level_high,final_level=final_level)


      if(n_cycl>1){
        diff_vect <- last_wv -value_node_gen(watervalues,statesdt,reservoir)$vu
        convergence_value <- converged(diff_vect,conv=convergence_criteria)
        convergence_percent <- sprintf((convergence_value*100), fmt = '%#.2f')

        if(n_cycl>2){
        if(convergence_value>convergence_rate){
            cat(paste0("\033[0;42m", "        Cycle number:", n_cycl, ",    ==>",convergence_percent,"% Converged <==", "\033[0m \n"))
          break}


          if (identical(diff_vect,last_diff_vect)){
          cat(paste0("\033[0;43m", "        Cycle number:", n_cycl, ",    ==>", "Converged with:",convergence_percent,"% <==", "\033[0m \n"))
          break}


        if(convergence_value==last_conv){
          counter <- counter+1
          if(counter>1){
            cat(paste0("\033[0;43m", "        Cycle number:", n_cycl, ",    ==>", "Converged with:",convergence_percent,"% <==", "\033[0m \n"))
            break

          }
        }
        }


        cat(paste0("\033[0;41m", "        Cycle number:", n_cycl, ",    ==>",convergence_percent,"% Converged <==", "\033[0m \n"))
        last_diff_vect <- diff_vect
        last_conv <- convergence_value

      }

      last_wv <- value_node_gen(watervalues,statesdt,reservoir)$vu

      }
    }# end else


  if(shiny){
    for (p in c("shinybusy")){
      if (!requireNamespace(p, quietly = TRUE)) {
        stop(
          paste0("Packageb", p, " must be installed to use this function."),
          call. = FALSE
        )
      }
    }
    shinybusy::remove_modal_spinner()

  }

  # Prepare output
  result <- list()
  result$watervalues <- watervalues
  result$aggregated_results <- value_nodes_dt
  class(result) <- "detailled and aggregated results of watervalues calculation"
  return(result)
}