R/Bellman.R

In dhia-gharsallaoui/watervalues: Generating water values for Antares

  #' Compute Bellman values at step i from step i+1
  #'
  #' @param Data_week A "data.table" generated in Grid_Matrix code
  #' that contains:
  #'  * states Numeric. All the water values that can be set, listed in
  #'   decreasing order.
  #'  * value_inflow Numeric. Inflow values for each Monte-Carlo year.
  #'  * Rewards for each simulation value and each Monte-Carlo year.
  #'  * level_high Numeric. Highest possible reservoir value.
  #'  * level_low Numeric. Lowest possible reservoir value.
  #'  * states_next List of vectors enumerating all reachable states
  #' @param decision_space Simulation constraints values
  #' @param next_week_values_l Numeric. Bellman values at step i+1.
  #' @param niveau_max Numeric of length 1. Reservoir capacity.
  #' @param E_max Numeric of length 1. Maximum energy that can be generated by
  #'   hydro storage over one step of time.
  #' @param P_max Numeric of length 1. Maximum energy that can be pumped to
  #' reservoir over one step of time.
  #' @param method Character. Perform mean of grids algorithm or grid of means algorithm or
  #'  grid of quantile algorithm.
  #' @param na_rm Boolean. Remove NAs
  #' @param max_mcyear Numeric of length 1. Number of Monte-Carlo year
  #'  used in simulations
  #' @param q_ratio numeric in [0,1]. the probability used in quantile algorithm.
  #' @param print_test Boolean. print Bellman values.
  #' @param test_week Numeric of length 1. number of the week to print in test.
  #' @param counter Numeric of length 1. number of the week in calculation.
  #' @param correct_outliers If TRUE, outliers in Bellman values are replaced by spline
  #'   interpolations. Defaults to FALSE.
  #' @param inaccessible_states Numeric in [0,1]. Tolerance of inaccessible states.
  #' For example if equal to 0.9 we delete the state if this states is inaccessible by 90\% of scenarios.
  #' @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 debugger_feas open debug mode in case there is an error of no accessible states
  #' @param ... further arguments passed to or from other methods.
  #' @return a \code{data.table} like Data_week with the Bellman values
  #' @importFrom stats ave quantile
  #' @export


  Bellman <- function(Data_week,next_week_values_l,decision_space,E_max,P_max=0,
                      niveau_max,method,na_rm=TRUE,max_mcyear,print_test=FALSE,
                      correct_outliers=FALSE,q_ratio=0.75,test_week,counter,inaccessible_states=1,
                      stop_rate=5,debugger_feas=F,...){




    decision_space <- unlist(decision_space, use.names = FALSE)
    decision_space <- round(decision_space)
    alpha <- getOption(x = "watervalues.alpha", default = 0.0001)
    decimals <- getOption(x = "watervalues.decimals", default = 3)

    {
    for (i in (1:nrow(Data_week))){
      #init
      {
        states <- Data_week$states[i]
        statesID <- Data_week$statesid[i]
        level_high <- Data_week$level_high[i]
        level_low <- Data_week$level_low[i]
        value_inflow <- Data_week$hydroStorage[i]
        reward <- Data_week$reward[[i]]
        states_next <- Data_week$states_next[[i]]
        value_reward <- unlist(reward, use.names = FALSE)
        states_next <- unlist(states_next, use.names = FALSE)


        if (i%% (nrow(Data_week)/max_mcyear) ==1){
          next_week_values <- next_week_values_l[i:(i+(nrow(Data_week)/max_mcyear)-1)]
        }



        }



    #eliminate non accessible states
    if (states > round(level_high, decimals) - alpha) {
      Data_week$value_node[i] <- -Inf
      next
    }
    if (states < round(level_low, decimals) + alpha) {
      Data_week$value_node[i] <- -Inf
      next
    }


    # max possible decision --------

    largest_decisions <- largest_decisions(states,value_inflow,niveau_max,E_max,P_max)

    largest_turb <-largest_decisions$largest_turb

    largest_pump <- largest_decisions$largest_pump

    if (largest_turb < largest_pump) {
      Data_week$value_node[i] <- -Inf
      message <- sprintf("Week %d in states %d even using the max power of turbining we can't avoid to cross the reservoir capacity",counter,statesID)
      message(message)
      next
    }

    # the decisions that respect the max possible decision used in simulation constraints

    decisions_current <- check_largest_decisions(decision_space,largest_decisions,alpha)


    # Possible next states
    next_states <- states_next[states_next >= (states - E_max + value_inflow-alpha) & states_next <= (states + P_max + value_inflow + alpha) ]


    # Turbaned energy per transition

    turbined_energy <- turbined_energy(states,next_states,value_inflow,decisions_current,largest_decisions)

    #add element to decisions_current to cover all needed information in the future
    decisions_cover <-decisions_cover(turbined_energy,decisions_current)

    # List of accessible Rewards

    step_reward <- accessible_rewards(decisions_cover,decision_space,value_reward)

    provisional_steps <- step_reward$steps
    provisional_reward_line <- step_reward$rewards


    decisions <- generate_decisions(turbined_energy,decisions_cover,E_max,P_max)



    decision_rewards <- generate_decisions_rewards(decisions,step_reward,alpha)


    # respect the guide graph constraints

    decisions <-  guide_cs_check(decisions,states,value_inflow,level_high,level_low,alpha)


    # Bellman calculator
    Bellman <- bellman_calculator(decisions,next_week_values,decision_rewards,states,value_inflow,niveau_max,states_next,alpha,na_rm)

    Bellman_values <- Bellman$Bellman_values



    max_Bell <- -Inf
    ld <-length(decisions)
    sorted_Bellman <- sort(Bellman_values)
    lb <- length(sorted_Bellman)
    t <- 1
    while (t<ld){
      max_bell <- sorted_Bellman[lb-t]
      if(is.finite(suppressWarnings(min(Bellman$next_bellman_value[which(Bellman_values==max_bell)])))){
        max_Bell <- max_bell
        break
      }

      t <- t+1
    }

    Data_week$value_node[i] <- max_Bell
    if(length(decisions)>0){
      Data_week$transition[i] <- min(decisions[which(Bellman_values==max_Bell)])
      Data_week$transition_reward[i] <- min(Bellman$transition_reward[which(Bellman_values==max_Bell)])
      Data_week$next_bellman_value[i] <- min(Bellman$next_bellman_value[which(Bellman_values==max_Bell)])
      if(is.na(min(Bellman$next_bellman_value[which(Bellman_values==max_Bell)])))browser()
    }

  #----- little test -----

    if(is.numeric(test_week)){
     if (counter==test_week){
     print(sprintf("******Week %d *********",counter))
     print(sprintf("State: %d",states))
     print(sprintf("StateID: %d",statesID))
     print(sprintf("Decisions :  "))
     print(decisions)
     print(sprintf("Rewards : "))
     print(Bellman_values)
     print(sprintf("BELLMAN >>>>> %f",Data_week$value_node[i]))
     writeLines("#######--------------------###########\n")
   }}


    }
    }


  # test feasible week
    feasible_test_week(Data_week$value_node,counter,stop_rate,debug_feas=debugger_feas)

  # test scenarios
    Data_week <- scanarios_check(Data_week,counter)

  # inaccessible criteria
    mcyears <- length(unique(Data_week$years))
    Data_week[,accessibility_percent:=accessibility/mcyears]
    Data_week[accessibility_percent<inaccessible_states,value_node:=NaN]
  #------ mean-grid method---------

  if (method == "mean-grid") {
    if (correct_outliers) {
      Data_week[, value_node := correct_outliers(value_node), by = years]
    }
      return(Data_week)
  }

  #------ grid-mean method---------

  if(method=="grid-mean"){
    if (correct_outliers) {
      Data_week[, value_node := correct_outliers(value_node)]
    }
      Data_week$value_node <- stats::ave(Data_week$value_node, Data_week$statesid, FUN=function(x) mean_or_inf(x,inaccessible_states))

    return(Data_week)
  }

  if (method=="quantile"){
    if (correct_outliers) {
      Data_week[, value_node := correct_outliers(value_node)]
    }

    Data_week$value_node <- stats::ave(Data_week$value_node, Data_week$statesid, FUN=function(x) stats::quantile(x, q_ratio,na.rm =T))


    return(Data_week)
  }




     return(Data_week)

    }







  #' Compute Bellman values at step i from step i+1 (version used in parallel computing)
  #'
  #' @param Data_week A "data.table" generated in Grid_Matrix code
  #' that contains:
  #'  * states Numeric. All the water values that can be set, listed in
  #'   decreasing order.
  #'  * value_inflow Numeric. Inflow values for each Monte-Carlo year.
  #'  * Rewards for each simulation value and each Monte-Carlo year.
  #'  * level_high Numeric. Highest possible reservoir value.
  #'  * level_low Numeric. Lowest possible reservoir value.
  #'  * states_next List of vectors enumerating all reachable states
  #' @param i the row number of the point to calculate her bellman value.
  #' @param decision_space Simulation constraints values
  #' @param next_week_values_l Numeric. Bellman values at step i+1.
  #' @param niveau_max Numeric of length 1. Reservoir capacity.
  #' @param E_max Numeric of length 1. Maximum energy that can be generated by
  #'   hydro storage over one step of time.
  #' @param P_max Numeric of length 1. Maximum energy that can be pumped to
  #' reservoir over one step of time.
  #' @param na_rm Boolean. Remove NAs
  #' @param max_mcyear Numeric of length 1. Number of Monte-Carlo year
  #'  used in simulations
  #' @param alpha maximum tolerance.
  #' @param decimals number of decimals.
  #' @return a \code{data.table} like Data_week with the Bellman values
  #' @importFrom stats ave quantile
  #' @export

  bellman_parallel_value <- function(Data_week,i,next_week_values_l,max_mcyear,
                                     E_max,P_max,niveau_max,na_rm
                                     ,decision_space,alpha,decimals){


  #init
  {
    states <- Data_week$states[i]
    statesID <- Data_week$statesid[i]
    level_high <- Data_week$level_high[i]
    level_low <- Data_week$level_low[i]
    value_inflow <- Data_week$hydroStorage[i]
    reward <- Data_week$reward[[i]]
    states_next <- Data_week$states_next[[i]]


    value_reward <- unlist(reward, use.names = FALSE)
    states_next <- unlist(states_next, use.names = FALSE)


    next_week_values <- next_week_value(next_week_values_l,max_mcyear,i)

  }



  #eliminate non accessible states
  if (states > round(level_high, decimals) - alpha) {
    return( -Inf)

  }
  if (states < round(level_low, decimals) + alpha) {
    return(-Inf)

  }

    # max possible decision --------

    largest_decisions <- largest_decisions(states,value_inflow,niveau_max,E_max,P_max)

    largest_turb <-largest_decisions$largest_turb

    largest_pump <- largest_decisions$largest_pump



    # the decisions that respect the max possible decision used in simulation constraints

    decisions_current <- check_largest_decisions(decision_space,largest_decisions,alpha)


    # Possible next states
    next_states <- states_next[states_next >= (states - E_max + value_inflow) & states_next <= (states + P_max + value_inflow + alpha) ]



    # Turbaned energy per transition

    turbined_energy <- turbined_energy(states,next_states,value_inflow,decisions_current,largest_decisions)

    #add element to decisions_current to cover all needed information in the future
    decisions_cover <-decisions_cover(turbined_energy,decisions_current)

    # List of accessible Rewards

    step_reward <- accessible_rewards(decisions_cover,decision_space,value_reward)

    provisional_steps <- step_reward$steps
    provisional_reward_line <- step_reward$rewards


    decisions <- generate_decisions(turbined_energy,decisions_cover,E_max,P_max)



    decision_rewards <- generate_decisions_rewards(decisions,step_reward,alpha)


    # respect the guide graph constraints

    decisions <-  guide_cs_check(decisions,states,value_inflow,level_high,level_low,alpha)


    # Bellman calculator
    Bellman_values <- bellman_calculator(decisions,next_week_values,decision_rewards,states,value_inflow,niveau_max,states_next,alpha,na_rm)







  return(max(Bellman_values, na.rm = TRUE))

  }