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 hours If method_old=F, vector of hours used to evaluate costs/rewards of pumping/generating
#' @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
#'
#' @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, hours=0:168, possible_controls = NULL,
                       max_hydro, mcyears = "all",area=NULL,pump_eff=NULL,
                       district_balance="water values district") {
  assertthat::assert_that(class(opts) == "simOptions")
  assertthat::assert_that(district_name %in% antaresRead::getDistricts(opts=opts))
  studyPath <- opts$studyPath

  # 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)
        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::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)){
      nb_hours <- length(hours)
      possible_controls <- data.frame(expand.grid(h=1:(2*nb_hours-1),week=1:52)) %>%
        dplyr::left_join(get_max_hydro(area,opts,"weekly"),by=c("week"="timeId")) %>%
        dplyr::mutate(i=dplyr::if_else(.data$h>nb_hours,2*nb_hours-.data$h,.data$h)) %>%
        dplyr::mutate(u=dplyr::if_else(.data$h>nb_hours,.data$turb/168*(168-hours[.data$i]),.data$pump/168*(hours[.data$i]-168)*pump_eff)) %>%
        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) {
        if (min(max_hydro$P_max)>0){
          res <- get_local_reward(o,hours,possible_controls,max_hydro,area,mcyears,
                                  district_balance,pump_eff)
        } else {
          res <- get_local_reward_turb(o,possible_controls,max_hydro,area,mcyears,
                                       district_balance)
        }
        res <- reward_offset(o,res, u,mcyears,district_name)
        res
      },
      o = opts_o,
      u = u,
      SIMPLIFY = F
    )}


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

    # 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$simulation_names <- simulation_names
    output$simulation_values <- possible_controls
  }

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

  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 hours vector of hours used to evaluate costs/rewards of pumping/generating
#' @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
#'
#' @return a data.table {mcYear,week,u,reward}
#' @export
get_local_reward <- function(opts,hours,possible_controls,max_hydro,area_price,mcyears,
                             district_balance="water values district",pump_eff=1){
  # Transform hours in negative hours
  hours <- unique(c(-hours, hours))

  # 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")


  # Evaluate how much reward one can make by generating at the maximum power for each hour of the week
  price_turb_more <- price %>%
    dplyr::group_by(.data$mcYear,.data$week) %>% dplyr::arrange(dplyr::desc(.data$price),.data$balance) %>%
    dplyr::left_join(max_hydro,by=c("timeId")) %>%
    dplyr::mutate(reward_turb=cumsum(price*dplyr::if_else(.data$balance<0,(.data$T_max+.data$balance),.data$T_max)),
                  vol_turb=cumsum(dplyr::if_else(.data$balance<0,(.data$T_max+.data$balance),.data$T_max)),
                  hour_turb=.data$vol_turb/.data$T_max) %>%
    dplyr::select("mcYear","week","hour_turb","reward_turb") %>% dplyr::ungroup()
  # Evaluate how much reward one can loose by not generating for each hour of the week
  price_turb_less <- price %>%
    dplyr::group_by(.data$mcYear,.data$week) %>%
    dplyr::left_join(max_hydro,by=c("timeId")) %>%
    dplyr::arrange(.data$price,dplyr::desc(.data$balance)) %>%
    dplyr::mutate(reward_turb=cumsum(price*dplyr::if_else(.data$balance<0,.data$balance,0L)),
                  vol_turb=cumsum(dplyr::if_else(.data$balance<0,.data$balance,0L)),
                  hour_turb=.data$vol_turb/.data$T_max) %>%
    dplyr::select("mcYear","week","hour_turb","reward_turb") %>% dplyr::ungroup()

  price_turb <- rbind(price_turb_less,price_turb_more) %>%
    dplyr::distinct(.data$mcYear,.data$week,.data$hour_turb,.data$reward_turb)

  # Transforming price into intervals of hours and associated reward
  price_turb_int <- price_turb %>%
    dplyr::group_by(.data$mcYear,.data$week) %>%
    dplyr::arrange(.data$hour_turb) %>%
    dplyr::mutate(hour_turb_inf=round(.data$hour_turb,5),reward_turb_inf=.data$reward_turb,
           hour_turb_sup=round(dplyr::lead(.data$hour_turb),5),reward_turb_sup=dplyr::lead(.data$reward_turb)) %>%
    dplyr::select(-c("hour_turb","reward_turb")) %>%
    tidyr::drop_na()

  # Evaluate how much has been generated in the simulation
  hour_turb_0 <- price_turb %>% dplyr::group_by(.data$mcYear,.data$week) %>%
    dplyr::summarise(hour_turb_0=-round(min(.data$hour_turb),5),.groups="drop")

  # Evaluate for each hour h of hours the variation of gain if the number of generating hours
  # varies by h
  hour_turb <- data.frame(tidyr::expand_grid(mcYear=mcyears,week=1:52,hour_turb=hours)) %>%
    rbind(dplyr::select(dplyr::mutate(hour_turb_0,hour_turb=round(-.data$hour_turb_0,5)),-c("hour_turb_0"))) %>%
    rbind(dplyr::select(dplyr::mutate(hour_turb_0,hour_turb=round(168-.data$hour_turb_0,5)),-c("hour_turb_0"))) %>%
    dplyr::left_join(hour_turb_0,by=c("mcYear","week")) %>%
    dplyr::filter(.data$hour_turb+.data$hour_turb_0>=0,.data$hour_turb+.data$hour_turb_0<=168) %>%
    dplyr::select(-c("hour_turb_0")) %>%
    dplyr::left_join(price_turb_int, by=dplyr::join_by(x$mcYear==y$mcYear,
                                                x$week==y$week,
                                                x$hour_turb>=y$hour_turb_inf,
                                                x$hour_turb<=y$hour_turb_sup)) %>%
    dplyr::distinct(.data$mcYear,.data$week,.data$hour_turb,.keep_all=T) %>%
    dplyr::mutate(reward_turb=dplyr::if_else(.data$hour_turb_inf!=.data$hour_turb_sup,
                               .data$reward_turb_inf+(.data$reward_turb_sup-.data$reward_turb_inf)/(.data$hour_turb_sup-.data$hour_turb_inf)*(.data$hour_turb-.data$hour_turb_inf),
                               .data$reward_turb_inf)) %>%
    dplyr::select(-c("hour_turb_inf","hour_turb_sup","reward_turb_inf","reward_turb_sup"))

  # Transforming into interval of hours and associated reward
  price_turb_int <- hour_turb %>%
    dplyr::group_by(.data$mcYear,.data$week) %>%
    dplyr::arrange(.data$hour_turb) %>%
    dplyr::mutate(hour_turb_inf=round(.data$hour_turb,5),reward_turb_inf=.data$reward_turb,
           hour_turb_sup=round(dplyr::lead(.data$hour_turb),5),reward_turb_sup=dplyr::lead(.data$reward_turb)) %>%
    dplyr::select(-c("hour_turb","reward_turb")) %>%
    tidyr::drop_na()

  # We next do exactly the same for pumping
  price_pump_more <- price %>%
    dplyr::group_by(.data$mcYear,.data$week) %>%
    dplyr::left_join(max_hydro,by=c("timeId")) %>%
    dplyr::arrange(.data$price,dplyr::desc(.data$balance)) %>%
    dplyr::mutate(cost_pump=cumsum(price*dplyr::if_else(.data$balance>0,(.data$P_max-.data$balance),.data$P_max)),
                  vol_pump=cumsum(dplyr::if_else(.data$balance>0,(.data$P_max-.data$balance),.data$P_max)),
                  hour_pump=.data$vol_pump/.data$P_max) %>%
    dplyr::select("mcYear","week","hour_pump","cost_pump") %>% dplyr::ungroup()
  price_pump_less <- price %>%
    dplyr::group_by(.data$mcYear,.data$week) %>%
    dplyr::left_join(max_hydro,by=c("timeId")) %>%
    dplyr::arrange(dplyr::desc(.data$price),.data$balance) %>%
    dplyr::mutate(cost_pump=cumsum(price*dplyr::if_else(.data$balance>0,-.data$balance,0L)),
                  vol_pump=cumsum(dplyr::if_else(.data$balance>0,-.data$balance,0L)),
                  hour_pump=.data$vol_pump/.data$P_max) %>%
    dplyr::select("mcYear","week","hour_pump","cost_pump") %>% dplyr::ungroup()
  price_pump <- rbind(price_pump_less,price_pump_more) %>%
    dplyr::distinct(.data$mcYear,.data$week,.data$hour_pump,.data$cost_pump)

  price_pump_int <- price_pump %>%
    dplyr::group_by(.data$mcYear,.data$week) %>%
    dplyr::arrange(.data$hour_pump) %>%
    dplyr::mutate(hour_pump_inf=round(.data$hour_pump,5),cost_pump_inf=.data$cost_pump,
           hour_pump_sup=round(dplyr::lead(.data$hour_pump),5),cost_pump_sup=dplyr::lead(.data$cost_pump)) %>%
    dplyr::select(-c("hour_pump","cost_pump")) %>%
    tidyr::drop_na()

  hour_pump_0 <- price_pump_int %>% dplyr::group_by(.data$mcYear,.data$week) %>%
    dplyr::summarise(hour_pump_0=-round(min(.data$hour_pump_inf),5),.groups="drop")

  hour_pump <- data.frame(tidyr::expand_grid(mcYear=mcyears,week=1:52,hour_pump=hours)) %>%
    rbind(dplyr::select(dplyr::mutate(hour_pump_0,hour_pump=round(-.data$hour_pump_0,5)),-c("hour_pump_0"))) %>%
    rbind(dplyr::select(dplyr::mutate(hour_pump_0,hour_pump=round(168-.data$hour_pump_0,5)),-c("hour_pump_0"))) %>%
    dplyr::left_join(hour_pump_0,by=c("mcYear","week")) %>%
    dplyr::filter(.data$hour_pump+.data$hour_pump_0>=0,.data$hour_pump+.data$hour_pump_0<=168) %>%
    dplyr::select(-c("hour_pump_0")) %>%
    dplyr::left_join(price_pump_int, by=dplyr::join_by(x$mcYear==y$mcYear,
                                                x$week==y$week,
                                                x$hour_pump>=y$hour_pump_inf,
                                                x$hour_pump<=y$hour_pump_sup)) %>%
    dplyr::distinct(.data$mcYear,.data$week,.data$hour_pump,.keep_all=T) %>%
    dplyr::mutate(cost_pump=dplyr::if_else(.data$hour_pump_inf!=.data$hour_pump_sup,
                             .data$cost_pump_inf+(.data$cost_pump_sup-.data$cost_pump_inf)/(.data$hour_pump_sup-.data$hour_pump_inf)*(.data$hour_pump-.data$hour_pump_inf),
                             .data$cost_pump_inf)) %>%
    dplyr::select(-c("hour_pump_inf","hour_pump_sup","cost_pump_inf","cost_pump_sup"))

  price_pump_int <- hour_pump %>%
    dplyr::group_by(.data$mcYear,.data$week) %>%
    dplyr::arrange(.data$hour_pump) %>%
    dplyr::mutate(hour_pump_inf=round(.data$hour_pump,5),cost_pump_inf=.data$cost_pump,
           hour_pump_sup=round(dplyr::lead(.data$hour_pump),5),cost_pump_sup=dplyr::lead(.data$cost_pump)) %>%
    dplyr::select(-c("hour_pump","cost_pump")) %>%
    tidyr::drop_na()

  # Maximum generating and pumping powers per hour
  max_hydro <- max_hydro %>%
    dplyr::mutate(week=(.data$timeId-1)%/%168+1) %>%
    dplyr::group_by(week) %>%
    dplyr::summarise(P_max=mean(.data$P_max),T_max=mean(.data$T_max),.groups = "drop")

  # Initialize a data.frame with one line per MC year per control
  if ("mcYear" %in% names(possible_controls)){
    df_reward <- possible_controls
    assertthat::assert_that(identical(unique(possible_controls$mcYear),mcyears))
  } else {
    df_reward <- data.frame(tidyr::expand_grid(mcYear=mcyears,possible_controls))
  }
  df_reward <- df_reward %>%
    dplyr::left_join(max_hydro,by=c("week")) %>%
    dplyr::left_join(hour_turb_0,by=c("mcYear","week")) %>%
    dplyr::left_join(hour_pump_0,by=c("mcYear","week"))

  # Evaluate extreme generating and pumping, ie for each control, the adequate amount of hours
  # pumping and generating such that the control corresponded and the total number of hours is 168
  df_reward_extreme <- df_reward %>%
    dplyr::mutate(hour_turb_exact=round((.data$u+168*.data$P_max*pump_eff)/(.data$T_max+.data$P_max*pump_eff)-.data$hour_turb_0,5),
           hour_pump_exact=168-.data$hour_turb_exact-.data$hour_pump_0-.data$hour_turb_0) %>%
    dplyr::left_join(price_turb_int, by=dplyr::join_by(x$mcYear==y$mcYear,
                                                x$week==y$week,
                                                x$hour_turb_exact>=y$hour_turb_inf,
                                                x$hour_turb_exact<=y$hour_turb_sup)) %>%
    dplyr::mutate(reward_turb=dplyr::if_else(.data$hour_turb_inf!=.data$hour_turb_sup,
                               .data$reward_turb_inf+(.data$reward_turb_sup-.data$reward_turb_inf)/(.data$hour_turb_sup-.data$hour_turb_inf)*(.data$hour_turb_exact-.data$hour_turb_inf),
                               .data$reward_turb_inf),
           reward_turb=round(.data$reward_turb,5)) %>%
    dplyr::select(-c("hour_turb_inf","hour_turb_sup","reward_turb_inf","reward_turb_sup")) %>%
    dplyr::rename(hour_turb="hour_turb_exact") %>%
    dplyr::distinct(.data$mcYear,.data$week,.data$u,.keep_all = T)

  # Calculate from hour_turb how much pumping is necessary to correspond to control and calculate the reward
  df_reward_turb <- dplyr::full_join(hour_turb, df_reward, by=c("mcYear","week"),relationship = "many-to-many") %>%
    dplyr::mutate(hour_pump_exact=round((-.data$u+.data$T_max*(.data$hour_turb_0+.data$hour_turb))/pump_eff/.data$P_max-.data$hour_pump_0,5)) %>%
    dplyr::filter((.data$hour_pump_exact+.data$hour_pump_0>=0)&(.data$hour_turb+.data$hour_turb_0+.data$hour_pump_exact+.data$hour_pump_0<=168)) %>%
    rbind(df_reward_extreme) %>%
    dplyr::distinct(.data$mcYear,.data$week,.data$hour_turb,.data$u,.keep_all = T) %>%
    dplyr::left_join(price_pump_int, by=dplyr::join_by(x$mcYear==y$mcYear,
                                                x$week==y$week,
                                                x$hour_pump_exact>=y$hour_pump_inf,
                                                x$hour_pump_exact<=y$hour_pump_sup)) %>%
    dplyr::mutate(cost_pump=dplyr::if_else(.data$hour_pump_inf!=.data$hour_pump_sup,
                             .data$cost_pump_inf+(.data$cost_pump_sup-.data$cost_pump_inf)/(.data$hour_pump_sup-.data$hour_pump_inf)*(.data$hour_pump_exact-.data$hour_pump_inf),
                             .data$cost_pump_inf),
           reward = .data$reward_turb-.data$cost_pump) %>%
    dplyr::select("mcYear","week","u","reward")

  # Same from hour_pump
  df_reward_pump <- dplyr::full_join(hour_pump, df_reward, by=c("mcYear","week"),relationship = "many-to-many") %>%
    dplyr::mutate(hour_turb_exact=round((.data$u+.data$P_max*pump_eff*(.data$hour_pump_0+.data$hour_pump))/.data$T_max-.data$hour_turb_0,5)) %>%
    dplyr::filter((.data$hour_turb_exact+.data$hour_turb_0>=0)&(.data$hour_turb_exact+.data$hour_turb_0+.data$hour_pump+.data$hour_pump_0<=168)) %>%
    dplyr::distinct(.data$mcYear,.data$week,.data$hour_pump,.data$u,.keep_all = T) %>%
    dplyr::left_join(price_turb_int, by=dplyr::join_by(x$mcYear==y$mcYear,
                                                x$week==y$week,
                                                x$hour_turb_exact>=y$hour_turb_inf,
                                                x$hour_turb_exact<=y$hour_turb_sup)) %>%
    dplyr::mutate(reward_turb=dplyr::if_else(.data$hour_turb_inf!=.data$hour_turb_sup,
                               .data$reward_turb_inf+(.data$reward_turb_sup-.data$reward_turb_inf)/(.data$hour_turb_sup-.data$hour_turb_inf)*(.data$hour_turb_exact-.data$hour_turb_inf),
                               .data$reward_turb_inf),
           reward_turb=round(.data$reward_turb,5),
           reward = .data$reward_turb-.data$cost_pump) %>%
    dplyr::select("mcYear","week","u","reward")

  # Bind the data.frame and calculate the best combination of generating and pumping hours for each control
  df_reward <- rbind(df_reward_pump,df_reward_turb) %>%
    dplyr::group_by(.data$mcYear,.data$week,.data$u) %>%
    dplyr::slice_max(.data$reward,with_ties = F) %>% dplyr::ungroup()

  return(df_reward)

}

#' Calculate rewards for a simulation based on marginal prices for a stock
#' with only generating power, 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 max_hydro data.frame {timeId,pump,turb} returned by the function \code{get_max_hydro}, should be hourly values
#'
#' @return a data.table {mcYear,week,u,reward}
#' @export
get_local_reward_turb <- function(opts,possible_controls,max_hydro,area_price,mcyears,
                                  district_balance="water values district"){
  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")


  price_turb_more <- price %>%
    dplyr::group_by(.data$mcYear,.data$week) %>%
    dplyr::left_join(max_hydro,by=c("timeId")) %>%
    dplyr::arrange(dplyr::desc(.data$price),.data$balance) %>%
    dplyr::mutate(reward_turb=cumsum(.data$price*dplyr::if_else(.data$balance<0,(.data$T_max+.data$balance),.data$T_max)),
                  vol_turb=cumsum(dplyr::if_else(.data$balance<0,(.data$T_max+.data$balance),.data$T_max)),
                  hour_turb=.data$vol_turb/.data$T_max) %>%
    dplyr::select("mcYear","week","hour_turb","reward_turb") %>% dplyr::ungroup()
  price_turb_less <- price %>%
    dplyr::group_by(.data$mcYear,.data$week) %>%
    dplyr::left_join(max_hydro,by=c("timeId")) %>%
    dplyr::arrange(.data$price,dplyr::desc(.data$balance)) %>%
    dplyr::mutate(reward_turb=cumsum(price*dplyr::if_else(.data$balance<0,.data$balance,0L)),
                  vol_turb=cumsum(dplyr::if_else(.data$balance<0,.data$balance,0L)),
                  hour_turb=.data$vol_turb/.data$T_max) %>%
    dplyr::select("mcYear","week","hour_turb","reward_turb") %>% dplyr::ungroup()
  price_turb <- rbind(price_turb_less,price_turb_more) %>%
    dplyr::distinct(.data$mcYear,.data$week,.data$hour_turb,.data$reward_turb)

  hour_turb_0 <- price_turb %>% dplyr::group_by(.data$mcYear,.data$week) %>%
    dplyr::summarise(hour_turb_0=-min(.data$hour_turb),.groups="drop")

  price_turb_int <- price_turb %>%
    dplyr::group_by(.data$mcYear,.data$week) %>%
    dplyr::arrange(.data$hour_turb) %>%
    dplyr::mutate(hour_turb_inf=round(.data$hour_turb,5),reward_turb_inf=.data$reward_turb,
           hour_turb_sup=round(dplyr::lead(.data$hour_turb),5),reward_turb_sup=dplyr::lead(.data$reward_turb)) %>%
    dplyr::select(-c("hour_turb","reward_turb")) %>%
    tidyr::drop_na()

  # hour_turb <- data.frame(tidyr::expand_grid(mcYear=mcyears,week=1:52,hour_turb=hours)) %>%
  #   dplyr::left_join(price_turb_int, by=dplyr::join_by(mcYear,week,hour_turb>=hour_turb_inf,
  #                                        hour_turb<=hour_turb_sup)) %>%
  #   dplyr::mutate(reward_turb=dplyr::if_else(hour_turb_inf!=hour_turb_sup,
  #                         reward_turb_inf+(reward_turb_sup-reward_turb_inf)/(hour_turb_sup-hour_turb_inf)*(hour_turb-hour_turb_inf),
  #                         reward_turb_inf)) %>%
  #   dplyr::select(-c(hour_turb_inf,hour_turb_sup,reward_turb_inf,reward_turb_sup))


  # Maximum generating and pumping powers per hour
  max_hydro <- max_hydro %>%
    dplyr::mutate(week=(.data$timeId-1)%/%168+1) %>%
    dplyr::group_by(week) %>%
    dplyr::summarise(P_max=mean(.data$P_max),T_max=mean(.data$T_max),.groups = "drop")

  if ("mcYear" %in% names(possible_controls)){
    df_reward <- possible_controls
    assertthat::assert_that(identical(unique(possible_controls$mcYear),mcyears))
  } else {
    df_reward <- data.frame(tidyr::expand_grid(mcYear=mcyears,possible_controls))
  }
  df_reward <- df_reward %>%
    dplyr::left_join(max_hydro,by=c("week")) %>%
    dplyr::left_join(hour_turb_0,by=c("mcYear","week"))

  df_reward <- df_reward %>%
    dplyr::mutate(hour_turb_exact=round(.data$u/.data$T_max-.data$hour_turb_0,5)) %>%
    dplyr::left_join(price_turb_int, by=dplyr::join_by(x$mcYear==y$mcYear,
                                                x$week==y$week,
                                                x$hour_turb_exact>=y$hour_turb_inf,
                                                x$hour_turb_exact<=y$hour_turb_sup)) %>%
    dplyr::mutate(reward=dplyr::if_else(.data$hour_turb_inf!=.data$hour_turb_sup,
                          .data$reward_turb_inf+(.data$reward_turb_sup-.data$reward_turb_inf)/(.data$hour_turb_sup-.data$hour_turb_inf)*(.data$hour_turb_exact-.data$hour_turb_inf),
                          .data$reward_turb_inf),
           reward=round(.data$reward,5)) %>%
    dplyr::select(-c("hour_turb_inf","hour_turb_sup","reward_turb_inf","reward_turb_sup")) %>%
    dplyr::rename(hour_turb="hour_turb_exact") %>%
    dplyr::select("mcYear","week","u","reward") %>%
    dplyr::distinct(.data$mcYear,.data$week,.data$u,.keep_all = T)

  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
#'
#' @return a data.table {mcYear,week,u,reward}
#' @export
reward_offset <- function(opts, df_reward, u0=c(),mcyears,district_cost= "water values district"){
  cost <- get_weekly_cost(district = district_cost, mcyears = mcyears, opts = opts) %>%
    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[[1]] %>%
      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)
}