R/get_Reward.R

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

#' Get the reward matrix from simulations, mainly used in \code{Grid_Matrix}
#'
#' @param simulation_names
#' generated by \code{runWaterValuesSimulation}
#' @param pattern A pattern to identify simulations.
#' @param district_name Name of the district used to store output.
#' @param opts
#'   List of simulation parameters returned by the function
#'   \code{antaresRead::setSimulationPath}
#' @param correct_monotony Binary argument (default to false). True to correct monotony of rewards.
#' @param method_old If T, linear interpolation used between simulations reward, else smarter interpolation based on marginal prices
#' @param possible_controls If method_old=F, data.frame {week,u} of controls evaluated per week
#' @param mcyears Vector of years used to evaluate rewards
#' @param area Area used to calculate watervalues
#' @param pump_eff Pumping efficiency
#' @param district_balance Name of district used to evaluate controls on the stock
#' @param simulation_values Values generated by \code{runWaterValuesSimulation}
#' @param max_hydro data.frame {timeId,pump,turb} returned by the function \code{get_max_hydro}, should be hourly values
#' @param expansion Binary. True if mode expansion was used to run simulations
#' @param fictive_areas Vector of chr. Fictive areas used in simulation
#'
#' @return list containing a data.table {timeid,MCyear,simulation overall cost}, list of simulations names and list of simulations values
#' @export
#'



get_Reward <- function(simulation_values = NULL,simulation_names=NULL, pattern = NULL,
                       district_name = "water values district",
                       opts = antaresRead::simOptions(), correct_monotony = FALSE,
                       method_old = TRUE, possible_controls = NULL,
                       max_hydro, mcyears = "all",area=NULL,pump_eff=NULL,
                       district_balance="water values district",
                       expansion=F,fictive_areas=NULL) {
  assertthat::assert_that(class(opts) == "simOptions")
  assertthat::assert_that(district_name %in% antaresRead::getDistricts(opts=opts))
  studyPath <- opts$studyPath
  area = tolower(area)

  # just a test if there is a simulation done or not
  {if (is.null(simulation_names)) {
    if (is.null(pattern))
      stop("If 'simulation_names' is not provided, 'pattern' cannot be NULL.")
    simulation_names <- getSimulationNames(pattern = pattern, studyPath = studyPath)
  }}

  # this part prepare the environment of each simulation
  {
    opts_o <- lapply(
      X = simulation_names,
      FUN = function(i) {
        suppressWarnings({
          antaresRead::setSimulationPath(path = studyPath, simulation = i)
        })
      }
    )
  }

  # check that the MC years are in simulations
  for (o in 1:length(opts_o)){
    assertthat::assert_that(all(mcyears %in% opts_o[[o]]$mcYears),
                            msg="Those MC years didn't have been all simulated, check your simulation.")
  }

  if(method_old){

    #generate a table containing the year, the time id (IE week) and overall cost
    {reward <- lapply(
      X = opts_o,
      FUN = function(o) {
        res <- get_weekly_cost(district = district_name, mcyears = mcyears, opts = o,
                               fictive_areas = fictive_areas, expansion = expansion)
        res$simulation <- o$name
        res
      }
    )}


    reward <- rbindlist(reward)   #merge the all simulations tables together

    # Getting the controls applied in each simulation
    decisions <- simulation_values %>%
      dplyr::mutate(sim=as.double(stringr::str_extract(.data$sim, "\\d+$")))

    # Joining controls to rewards
    reward <- reward %>%
      dplyr::mutate(sim=as.double(stringr::str_extract(.data$simulation, "\\d+$")))
    if ("mcYear" %in% names(decisions)){
      reward <- reward %>%
        dplyr::left_join(decisions,by=c("sim","timeId"="week","mcYear"))
    } else {
      reward <- reward %>%
        dplyr::left_join(decisions,by=c("sim","timeId"="week"))
    }
    reward <- reward %>%
      dplyr::rename(reward="ov_cost",control="u")

    if (correct_monotony){
      cost <- reward
      # Getting possible controls
      U <- cost %>%
        dplyr::select("control") %>%
        dplyr::distinct() %>%
        dplyr::arrange()
      # Initialize reward
      cost <- cost %>% dplyr::mutate(min_previous_reward=.data$reward) %>%
        dplyr::arrange(.data$mcYear, .data$timeId, .data$control)
      for (u in U$control){# for each control, and each MC year,
        # get the minimum reward for all possible controls smaller than u
        cost[cost$control==u,'min_previous_reward'] <- cost %>% dplyr::filter(.data$control<=u) %>%
          dplyr::group_by(.data$mcYear, .data$timeId) %>%
          dplyr::mutate(min_previous_reward = min(.data$reward)) %>%
          dplyr::ungroup() %>%
          dplyr::filter(.data$control==u) %>%
          dplyr::select("min_previous_reward")
      }

      cost <- cost %>% dplyr::select(-c("reward")) %>%
        dplyr::rename("reward" = "min_previous_reward")
      # replace values
      reward <- cost[,c("timeId","mcYear","reward","simulation","sim","control")]
    }

    # Retrieving reward for control 0 for each MC year and each week
    #  and subtracting this value to all rewards with same year and same week
    reward <- dplyr::filter(reward,.data$control==0) %>%
      dplyr::select("mcYear","timeId","reward") %>%
      dplyr::right_join(reward,by=c("mcYear","timeId"),suffix=c("_0","")) %>%
      # dplyr::mutate(reward=.data$reward_0-.data$reward) %>%
      dplyr::mutate(reward=-.data$reward) %>%
      dplyr::select(-c("reward_0","simulation","sim"))
    reward <- as.data.table(reward)

    options("antares" = opts)
    # Prepare output
    output <- list()
    output$reward <- reward
    output$simulation_names <- simulation_names
    output$simulation_values <- simulation_values

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

    max_hydro <- dplyr::rename(max_hydro,"P_max"="pump","T_max"="turb")
    assertthat::assert_that(min(max_hydro$T_max)>0)

    # Creating possible controls if none
    if(is.null(possible_controls)){
      possible_controls <- simulation_values %>%
        dplyr::select("week","u")
    }

    # Transforming simulation values such that for each week there is a line
    # and for each simulation there is a column
    u <- simulation_values %>%
      dplyr::mutate(sim=as.double(stringr::str_extract(.data$sim,"\\d+$"))) %>%
      dplyr::group_by(.data$sim) %>%
      tidyr::nest() %>%
      tidyr::pivot_wider(names_from=.data$sim,values_from=.data$data)

    # Interpolate reward for each simulation
    {reward <- mapply(
      FUN = function(o,u) {
        res <- get_local_reward(o,possible_controls,max_hydro,area,mcyears,u[[1]],
                                  district_balance,pump_eff)
        res <- reward_offset(o,res, u[[1]],mcyears,district_name,fictive_areas=fictive_areas, expansion=expansion)
        res <- dplyr::mutate(res,simulation=o$name)
        res
      },
      o = opts_o,
      u = u,
      SIMPLIFY = F
    )}


    reward <- rbindlist(reward)   #merge the all simulations tables together
    local_reward <- reward

    assertthat::assert_that(sum(is.na(local_reward))==0,
                            msg="NaN values in local reward, something went wrong.")

    # Getting the minimum reward for each year, each week and each control (u)
    reward <- reward %>%
      dplyr::group_by(.data$mcYear,.data$week,.data$u) %>% dplyr::summarise(reward=min(.data$reward),.groups="drop")
    #Subtracting the reward corresponding to control 0 for each year and each week
    reward <- dplyr::filter(reward,.data$u==0) %>% dplyr::select("mcYear","week","reward") %>%
      dplyr::right_join(reward,by=c("mcYear","week"),suffix=c("_0","")) %>%
      # dplyr::mutate(reward=.data$reward-.data$reward_0) %>%
      dplyr::rename("timeId"="week","control"="u") %>%
      dplyr::select(-c("reward_0"))
    reward <- as.data.table(reward)

    options("antares" = opts)
    # Prepare output
    output <- list()
    output$reward <- reward
    output$local_reward <- local_reward
    output$simulation_names <- simulation_names
    output$simulation_values <- possible_controls
  }

  class(output) <- "Reward matrix , simulation names and values"

  assertthat::assert_that(sum(is.na(reward))==0,
                          msg="NaN values in reward, something went wrong.")

  return(output)

}

#' Calculate rewards for a simulation based on marginal prices, mainly used in \code{Get_Reward}
#'
#' @param opts List of simulation parameters returned by the function
#'   \code{antaresRead::setSimulationPath}
#' @param possible_controls data.frame {week,u} of controls evaluated per week
#' @param area_price Area used to evaluate marginal prices
#' @param mcyears Vector of years used to evaluate rewards
#' @param district_balance Name of district used to evaluate controls on the stock
#' @param pump_eff Pumping efficiency
#' @param max_hydro data.frame {timeId,pump,turb} returned by the function \code{get_max_hydro}, should be hourly values
#' @param u0 data.table {week,u} Constraint values per week used in the simulation
#'
#' @return a data.table {mcYear,week,u,reward}
#' @export
get_local_reward <- function(opts,possible_controls,max_hydro,area_price,mcyears,u0,
                             district_balance="water values district",pump_eff=1){
  # Get hourly marginal prices and energy pumped and generated for each hour
  price <- antaresRead::readAntares(areas=area_price,select=c("MRG. PRICE"),
                       opts=opts,mcYears = mcyears, timeStep = "hourly") %>%
    dplyr::select(-c("day","month","hour","area","time")) %>%
    dplyr::left_join(antaresRead::readAntares(districts=district_balance,select=c("BALANCE"),
                          opts=opts,mcYears = mcyears, timeStep = "hourly"),by=c("timeId","mcYear")) %>%
    dplyr::select(-c("day","month","hour","district","time")) %>%
    dplyr::mutate(week=(.data$timeId-1)%/%168+1) %>%
    dplyr::rename(price="MRG. PRICE",balance="BALANCE") %>%
    dplyr::left_join(max_hydro,by=c("timeId")) %>%
    dplyr::mutate(balance = dplyr::if_else(.data$balance<(-.data$T_max),-.data$T_max,.data$balance),
                  balance = dplyr::if_else(.data$balance>(.data$P_max),.data$P_max,.data$balance))

  price <- price %>%
    dplyr::cross_join(data.frame(pumping=c(T,F))) %>%
    dplyr::mutate(price = dplyr::if_else(.data$pumping,.data$price/pump_eff,.data$price)) %>%
    dplyr::mutate(gap_greater_control = dplyr::if_else(.data$pumping,
                                       dplyr::if_else(.data$balance>0,.data$balance*pump_eff,0),
                                       dplyr::if_else(.data$balance<0,.data$T_max+.data$balance,.data$T_max)),
                  gap_lower_control = dplyr::if_else(.data$pumping,
                                      dplyr::if_else(.data$balance>0,(.data$balance-.data$P_max)*pump_eff,(-.data$P_max)*pump_eff),
                                      dplyr::if_else(.data$balance<0,.data$balance,0)))
  greater_control_reward <- price %>%
    dplyr::group_by(.data$mcYear,.data$week) %>%
    dplyr::arrange(dplyr::desc(.data$price),.data$balance) %>%
    dplyr::mutate(reward=cumsum(.data$price*.data$gap_greater_control),
                  dif_vol=cumsum(.data$gap_greater_control)) %>%
    dplyr::select("mcYear","week","dif_vol","reward") %>%
    dplyr::ungroup()

  lower_control_reward <- price %>%
    dplyr::group_by(.data$mcYear,.data$week) %>%
    dplyr::arrange(.data$price,.data$balance) %>%
    dplyr::mutate(reward=cumsum(.data$price*.data$gap_lower_control),
                  dif_vol=cumsum(.data$gap_lower_control)) %>%
    dplyr::select("mcYear","week","dif_vol","reward") %>%
    dplyr::ungroup()


  control_reward <- rbind(greater_control_reward,lower_control_reward) %>%
    dplyr::distinct(.data$week,.data$mcYear,.data$dif_vol,.data$reward)


  if (!("mcYear" %in% names(possible_controls))){
    possible_controls <- possible_controls %>%
      dplyr::cross_join(data.frame(mcYear=mcyears))
  }

  u0 <- u0 %>%
    dplyr::rename("u0"="u")
  if ("mcYear" %in% names(u0)){
    control_reward <- control_reward  %>%
      dplyr::left_join(u0,by=c("week","mcYear"))
  } else {
    control_reward <- control_reward %>%
      dplyr::left_join(u0,by=c("week"))
  }

  control_reward <- control_reward %>%
    dplyr::mutate(dif_vol = .data$dif_vol + .data$u0)%>%
    dplyr::group_by(.data$mcYear,.data$week) %>%
    dplyr::arrange(.data$dif_vol) %>%
    dplyr::mutate(marg = (dplyr::lead(.data$reward)-.data$reward)/(dplyr::lead(.data$dif_vol)-.data$dif_vol))%>%
    dplyr::mutate(dif_vol_sup = dplyr::lead(.data$dif_vol),reward_sup=dplyr::lead(.data$reward)) %>%
    dplyr::rename(dif_vol_inf = "dif_vol", reward_inf = "reward") %>%
    tidyr::drop_na() %>%
    dplyr::ungroup()

  extreme_control <- possible_controls %>%
    dplyr::group_by(.data$week,.data$mcYear) %>%
    dplyr::summarise(min_vol=min(.data$u),
                     max_vol=max(.data$u),.groups="drop")

  control_reward <- control_reward %>%
    dplyr::left_join(extreme_control,by=c("week","mcYear")) %>%
    dplyr::group_by(.data$mcYear,.data$week) %>%
    dplyr::mutate(reward_inf = dplyr::if_else(.data$dif_vol_inf==min(.data$dif_vol_inf),
          .data$reward_sup - .data$marg* (.data$dif_vol_sup-.data$min_vol),
          .data$reward_inf),
          dif_vol_inf = dplyr::if_else(.data$dif_vol_inf==min(.data$dif_vol_inf),
                                       .data$min_vol,
                                       .data$dif_vol_inf),
          reward_sup = dplyr::if_else(.data$dif_vol_sup==max(.data$dif_vol_sup),
                                      .data$reward_inf + .data$marg* (.data$max_vol-.data$dif_vol_inf),
                                      .data$reward_sup),
          dif_vol_sup = dplyr::if_else(.data$dif_vol_sup==max(.data$dif_vol_sup),
                                       .data$max_vol,
                                       .data$dif_vol_sup))

  df_reward <- possible_controls %>%
    dplyr::left_join(control_reward, by=dplyr::join_by(x$week==y$week,
                                                       x$mcYear==y$mcYear,
                                                       x$u>=y$dif_vol_inf,
                                                       x$u<=y$dif_vol_sup)) %>%
    dplyr::mutate(reward = dplyr::if_else(.data$dif_vol_sup!=.data$dif_vol_inf,
      (.data$reward_sup-.data$reward_inf)/(.data$dif_vol_sup-
      .data$dif_vol_inf)*(.data$u-.data$dif_vol_inf)+.data$reward_inf,
      .data$reward_sup)) %>%
    dplyr::distinct(.data$week,.data$mcYear,.data$u,.data$reward)%>%
    dplyr::select("mcYear","week","u","reward") %>%
    dplyr::ungroup()

  return(df_reward)

}

#' Modify local reward to take into account overall cost of the simulation, mainly used in \code{Get_Reward}
#'
#' @param opts List of simulation parameters returned by the function
#'   \code{antaresRead::setSimulationPath}
#' @param df_reward data.table computed by the function \code{get_local_reward}
#' @param u0 data.table {week,u} Constraint values per week used in the simulation, empty list if none
#' @param mcyears Vector of years used to evaluate rewards
#' @param district_cost Name of district used to evaluate overall cost
#' @param expansion Binary. True if mode expansion was used to run simulations
#' @param fictive_areas Vector of chr. Fictive areas used in simulation
#'
#' @return a data.table {mcYear,week,u,reward}
#' @export
reward_offset <- function(opts, df_reward, u0=c(),mcyears,district_cost= "water values district",
                          fictive_areas=NULL, expansion = F){
  cost <- get_weekly_cost(district = district_cost, mcyears = mcyears, opts = opts,
                          fictive_areas = fictive_areas, expansion=expansion) %>%
    dplyr::rename(week="timeId") %>%
    dplyr::select("mcYear","week","ov_cost") %>%
    as.data.frame()
  if (sum(is.na(u0))>=1){
    u0 <- c()
  }
  if (length(u0)>0){
    u0 <- u0 %>%
      dplyr::rename("u0"="u")
    if ("mcYear" %in% names(u0)){
      df_reward <- df_reward %>%
        dplyr::left_join(u0,by=c("week","mcYear"))
    } else {
      df_reward <- df_reward %>%
        dplyr::left_join(u0,by=c("week"))
    }
    df_reward <- df_reward %>%
      dplyr::left_join(dplyr::select(dplyr::filter(df_reward,.data$u==u0),
                       "mcYear","week","reward"),
                by=c("mcYear","week"),suffix=c("","_0")) %>%
      dplyr::mutate(reward = .data$reward-.data$reward_0) %>%
      dplyr::select(-c("reward_0","u0"))
  }
  df_reward <- df_reward %>%
    dplyr::left_join(cost,by=c("mcYear","week")) %>%
    dplyr::mutate(reward = .data$reward-.data$ov_cost) %>%
    dplyr::select(-c("ov_cost"))
  return(df_reward)
}