R/math_functions.R

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

#' Calculate and plot watervalues with \code{value_node_gen}. Penalties taken into account.
#' Used in \code{Grid_Matrix}
#'
#' @param watervalues Data frame generated in \code{Grid_Matrix}
#' @param statesdt Data frame of possible states, generated in \code{Grid_Matrix}
#' @param reservoir Data frame describing rule curves, generated by \code{readReservoirLevels}
#' @param penalty_level_high Penalty for violating the bottom rule curve
#' @param penalty_level_low Penalty for violating the top rule curve
#' @param force_final_level Binary. Whether final level should be constrained
#' @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
#' @param final_level Double. Final level to consider
#'
#' @return Data frame with water value (vu) for each week (weeks) and each state (states).
#' vu_pen corresponds to water value without penalties
build_data_watervalues <- function(watervalues,statesdt,reservoir,
                                   penalty_level_high,penalty_level_low,
                                   force_final_level=F,penalty_final_level_low=0,
                                   penalty_final_level_high=0,final_level=NULL){
  # calculate derivative of Bellman values (ie watervalues not yet penalized)
  value_nodes_dt <- value_node_gen(watervalues,statesdt,reservoir)

  # removing artificial week 53
  value_nodes_dt <- value_nodes_dt[value_nodes_dt$weeks!=1,]

  #add penalties
  if (!force_final_level){
    value_nodes_dt <- value_nodes_dt %>%
      dplyr::mutate(value_nodes_dt,
                    vu_pen=dplyr::case_when(states>.data$level_high ~ vu - penalty_level_high,
                                            states<.data$level_low ~ vu + penalty_level_low,
                                            TRUE ~ vu)) %>%
      dplyr::mutate(force_final_level=FALSE) %>%
      dplyr::mutate(penalty_low=penalty_level_low,
                    penalty_high=penalty_level_high) %>%
      dplyr::rename(vu="vu_pen",vu_pen="vu")
  } else {
    value_nodes_dt <- value_nodes_dt %>%
      dplyr::mutate(level_high=dplyr::if_else(.data$weeks!=53,
                                              .data$level_high,final_level)) %>%
      dplyr::mutate(level_low=dplyr::if_else(.data$weeks!=53,
                                             .data$level_low,final_level)) %>%
      dplyr::mutate(vu_pen=dplyr::if_else(.data$weeks!=53,
                                          dplyr::case_when(states>.data$level_high ~ vu - penalty_level_high,
                                            states<.data$level_low ~ vu + penalty_level_low,
                                            TRUE ~ vu),
                                          dplyr::case_when(states>.data$level_high ~ vu - penalty_final_level_high,
                                                           states<.data$level_low ~ vu + penalty_final_level_low,
                                                           TRUE ~ vu))) %>%
      dplyr::mutate(force_final_level=TRUE) %>%
      dplyr::mutate(penalty_low=dplyr::if_else(.data$weeks!=53,
                                               penalty_level_low,penalty_final_level_low),
                    penalty_high=dplyr::if_else(.data$weeks!=53,
                                                penalty_level_high,penalty_final_level_high)) %>%
      dplyr::rename(vu="vu_pen",vu_pen="vu")
  }

  value_nodes_dt <- value_nodes_dt %>%
    dplyr::mutate(weeks=.data$weeks-1)

  # plotting
  print(waterValuesViz(value_nodes_dt))
  return(value_nodes_dt)

}


#' Calculate water values from Bellman values, used in \code{Grid_Matrix} and
#' \code{build_data_watervalues}
#'
#' @param watervalues an intermediate result in Grid_Matrix contains the bellman values
#' @param statesdt an intermediate result in Grid_Matrix contains the states dicretization
#' @param reservoir an intermediate result in Grid_Matrix contains the reservoir levels
value_node_gen <- function(watervalues,statesdt,reservoir){
  # Calculating mean of Bellman values for each week and each state
  value_nodes_dt <- watervalues %>%
    dplyr::group_by(.data$weeks, .data$statesid) %>%
    dplyr::summarise(value_node=mean_finite(.data$value_node),.groups = "drop") %>%
    dplyr::mutate(value_node=dplyr::if_else(!is.finite(.data$value_node),
                                            NaN,.data$value_node)) %>%
    dplyr::left_join(statesdt,by=c("statesid","weeks")) %>%
    as.data.table()

  # Adding rule curves for the beginning of the week and then calculating derivative of Bellman values
  names(reservoir)[1] <- "weeks"
  value_nodes_dt <- dplyr::mutate(value_nodes_dt, beg_week=dplyr::if_else(.data$weeks>1,.data$weeks-1,52))
  value_nodes_dt <- dplyr::left_join(value_nodes_dt,reservoir,by=c("beg_week"="weeks")) %>%
    dplyr::select(-c("beg_week","level_avg")) %>%
    dplyr::arrange(.data$weeks, -.data$statesid) %>%
    dplyr::group_by(.data$weeks) %>%
    dplyr::mutate(value_node_dif=.data$value_node-dplyr::lag(.data$value_node),# Delta of Bellman values
                  states_dif=.data$states-dplyr::lag(.data$states), # Delta of states
                  vu=.data$value_node_dif /.data$states_dif, # Ratio
                  vu=round(.data$vu,2)) %>%
    as.data.table()
  return(value_nodes_dt)
}



#' test a difference vector convergence, used in \code{Grid_Matrix}
#' @param diff_vect is a vector of water values differences
#' @param conv is the value from which the difference become converged

converged <- function(diff_vect,conv=1){
  t <- abs(diff_vect)
  t2 <- is.nan(t)
  nan_values <- length(t2[t2==T])
  numeric_values <- length(t)-nan_values

  converged_values <- length(t[t<conv]) - nan_values
  converge_percent <-  converged_values/numeric_values
  return(converge_percent)
  }

#----- Mean of finite values
#' Calculate the mean of finite values.
#' Return \code{-Inf} if all \code{-Inf}.
#' @param x numeric vector whose mean is wanted.

mean_finite <- function(x) {
  if (all(!is.finite(x))) {
    -Inf
  } else {
    mean(x[is.finite(x)])
  }
}

#------- states to percent -----
#' Convert Reservoir levels from MWh to percent of reservoir.
#' @param data  A data.table contains the statesid and states columns
#' @param states_step_ratio percent step between two successive levels.

states_to_percent <- function(data,states_step_ratio=0.01){

  # rescale levels to round percentages ranging from 0 to 100
  states_ref <- data[, .SD[1], by = c("statesid"), .SDcols = "states"]
  states_ref[, "states_percent" := 100*states_ref$states/max(states_ref$states)]

  interv <- seq(from=0,to=100,by=round(100*states_step_ratio))
  nearest_states <- states_ref$statesid[sapply(interv, function(x) which.min(abs(x - states_ref$states_percent)))]

  states_ref_0_100 <- data.table(
    states_round_percent = interv,
    statesid = nearest_states
  )

  res <- CJ(weeks = unique(data$weeks), states_round_percent = interv)

  res <- res %>%
    dplyr::left_join(dplyr::select(states_ref_0_100,
                                   c("states_round_percent",
                                     "statesid")),by=c("states_round_percent")) %>%
    dplyr::left_join(dplyr::select(data,
                                   c("weeks", "statesid",
                                     "value_node","value_node_dif","vu")),
                                   by=c("weeks", "statesid")) %>%
    as.data.table()

  return(res)
  }


#' Correct concavity of Bellman values for the given weeks to have nice monotony for water values,
#' used in \code{Grid_Matrix}
#'
#' @param weeks Weeks for which we want to correct concavity
#' @param df_value_node DataFrame containing bellman values
#'
#' @return vector of corrected Bellman values
correct_concavity <- function(df_value_node, weeks){

  for (s in weeks){
    # Getting values for the week s
    df_week <- df_value_node[df_value_node$weeks==s,c("weeks","states","value_node","years")]
    # Removing NaN and infinite values
    df_week <- dplyr::filter(df_week,!is.na(df_week$value_node))
    df_week <- dplyr::filter(df_week,is.finite(df_week$value_node))
    df_week <- unique(df_week)
    df_week <- dplyr::arrange(df_week,.data$states,.data$years)
    n <- nrow(df_week)
    # Initialize corrected values
    df_week$new_value <- df_week$value_node
    # Getting all possible states
    states <- dplyr::distinct(df_week,states) %>% dplyr::pull(states)
    for (x in states){
      # Getting state x and its Bellman value, then calculating slope between point x
      # and all other points
      df_week <- df_week %>% dplyr::filter(states==x) %>%
        dplyr::rename(value_x="new_value",states_x="states") %>%
        dplyr::select("years", "states_x", "value_x") %>%
        dplyr::right_join(df_week, by=c("years")) %>%
        dplyr::mutate(coef=(.data$new_value-.data$value_x)/(states-.data$states_x)) %>%
        dplyr::mutate(coef=dplyr::if_else(.data$coef<0,0,.data$coef))
      # Keeping only states bigger than x then for each year, getting the maximum slope m and
      # state y associated with this slope. Then for each state between x and y, correct the value
      # such as the function became concave thanks to slope m
      df_week <- df_week %>% dplyr::filter(states>x) %>% dplyr::group_by(.data$years) %>%
        dplyr::slice(which.max(.data$coef)) %>%
        dplyr::transmute(m=.data$coef,states_y=.data$states, years=.data$years) %>%
        dplyr::right_join(df_week,by="years") %>%
        dplyr::mutate(new_value=dplyr::if_else((states>.data$states_x)&(states<=.data$states_y),
                                 .data$m*(states-.data$states_x)+.data$value_x,.data$new_value)) %>%
        dplyr::select("years", "weeks", "states", "value_node", "new_value")
    }
    # Replacing values
    df_value_node[df_value_node$weeks==s,"new_value"]  <- dplyr::left_join(df_value_node[df_value_node$weeks==s,c("weeks","states","years")],
                                                                    df_week[,c("weeks","states","new_value","years")],
                                                                    by=c("weeks","states","years"))$new_value
  }
  return(df_value_node$new_value)
}