Invest module on flexibility

NAMESPACE

@@ -1,10 +1,12 @@
 # Generated by roxygen2: do not edit by hand
 
 export(Grid_Matrix)
+export(MultiStock_H2_Investment_reward_compute_once)
 export(calculateBellmanWithIterativeSimulations)
 export(calculateBellmanWithIterativeSimulationsMultiStock)
 export(changeHydroManagement)
 export(constraint_generator)
+export(edit_study_with_results)
 export(getBellmanValuesFromOneSimulationMultistock)
 export(getBellmanValuesSequentialMultiStockWithPlaia)
 export(getBellmanValuesWithPlaia)
@@ -23,6 +25,7 @@ export(plot_reward_mc)
 export(plot_reward_variation)
 export(plot_reward_variation_mc)
 export(readReservoirLevels)
+export(remove_candidate_cluster)
 export(remove_out)
 export(reward_offset)
 export(runWaterValuesSimulation)

R/Bellman.R

@@ -8,8 +8,7 @@
   #'  * 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 reward Data.frame reward generated by \code{get_Reward()}
   #' @param next_week_values_l Numeric. Bellman values at step i+1.
   #' @param E_max Numeric of length 1. Maximum energy that can be generated by
   #'   hydro storage over one step of time.
@@ -27,47 +26,50 @@
   #'
   #' @return a \code{data.table} like Data_week with the Bellman values
   #' @keywords internal
-  Bellman <- function(Data_week,next_week_values_l,decision_space,E_max,P_max=0,
+  Bellman <- function(Data_week,next_week_values_l,reward,E_max,P_max=0,
                       mcyears,cvar_value=1,niveau_max,
                       penalty_level_low,penalty_level_high,
                       lvl_high,lvl_low,overflow_cost,
                       next_state=NULL){
 
     # Getting all possible transitions between a state for the current week and a state for the next week
-    decision_space <- dplyr::select(decision_space,-c("week"))
-    if (!("mcYear") %in% names(decision_space)){
-      decision_space <- decision_space %>%
-        dplyr::cross_join(data.frame(mcYear=mcyears))
-    }
+    decision_space <- dplyr::select(reward,-c("timeId","reward"))
 
     if (is.null(next_state)){
       next_state = Data_week$states
     }
 
-    # Get interpolation function of rewards for each possible transition for each MC year
-    f_reward_year <- get_reward_interpolation(Data_week)
-
     # Build a data.table from Data_week that list for each state and each MC year, the possible transitions
     df_SDP <- build_all_possible_decisions(Data_week,decision_space,
                                            lvl_high,lvl_low,E_max,P_max,
                                            next_week_values_l,niveau_max,overflow_cost,
                                            next_states = next_state)
 
+    # Reward interpolation
+    setDT(df_SDP)
+    setDT(reward)
+    interp_reward <- reward[
+      , .(interp_fun = list({
+        ctrl <- sort(unique(control))
+        rew  <- reward[order(control)][match(ctrl, sort(control))]
+        approxfun(ctrl, rew)
+      })),
+      by = mcYear
+    ]
+    df_SDP <- interp_reward[df_SDP, on = .(mcYear = years)][, reward := interp_fun[[1]](control), by = mcYear][, interp_fun := NULL]
+
     # For each transition (control), find the associated reward and for each next state,
     # calculate penalties for violating rule curves. Then, find for each MC year and each state,
     # the maximum sum of reward, next bellman value and penalties
     df_SDP <- df_SDP %>%
-      dplyr::mutate(gain=mapply(function(y,x)f_reward_year[[which(y==mcyears)]](x), df_SDP$years, df_SDP$control),
-             penalty_low = dplyr::if_else(.data$next_state<=lvl_low,penalty_level_low*(.data$next_state-lvl_low),0),
+      dplyr::mutate(penalty_low = dplyr::if_else(.data$next_state<=lvl_low,penalty_level_low*(.data$next_state-lvl_low),0),
              penalty_high = dplyr::if_else(.data$next_state>=lvl_high,penalty_level_high*(lvl_high-.data$next_state),0),
-             sum=.data$gain+.data$next_value+.data$penalty_low+.data$penalty_high) %>%
-      dplyr::group_by(.data$years,.data$states) %>%
+             sum=.data$reward+.data$next_value+.data$penalty_low+.data$penalty_high) %>%
+      dplyr::group_by(.data$mcYear,.data$states) %>%
       dplyr::filter(.data$sum==max(.data$sum)) %>%
       dplyr::slice_max(.data$next_state, with_ties = F) %>%
-      dplyr::select(-c("value_node","transition","transition_reward",
-                "next_bellman_value")) %>%
-      dplyr::rename("value_node"="sum","transition"="control","transition_reward"="gain",
-             "next_bellman_value"="next_value")
+      dplyr::rename("value_node"="sum","transition"="control","transition_reward"="reward",
+             "next_bellman_value"="next_value","years"="mcYear")
 
     assertthat::assert_that(nrow(df_SDP)==nrow(Data_week),msg=paste0("Problem with Bellman"))
 

R/Bellman_subfunctions.R

@@ -1,21 +1,3 @@
-
-#' Create approximation of reward function for each scenario, used in \code{Bellman}
-#'
-#' @param Data_week Data frame generated in \code{Grid_Matrix} code containing
-#' reward database for each scenario (we suppose there is only one week at a time)
-#'
-#' @return List of \code{stats::approxfun} for each scenario
-#' @keywords internal
-get_reward_interpolation <- function(Data_week){
-
-  reward <- dplyr::distinct(Data_week[,c('years','reward_db')])
-
-  f_reward_year <- sapply(reward$reward_db,FUN = function(df) stats::approxfun(df$control,df$reward))
-
-  return(f_reward_year)
-}
-
-
 #' Create a data frame with all possible transition and associated next state for
 #' each scenario and each state for one particular week, used in \code{Bellman}.
 #' For each next state, Bellman value is calculated.
@@ -61,7 +43,6 @@ build_all_possible_decisions <- function(Data_week,decision_space,
 
   possible_control <- Data_week  %>%
     dplyr::right_join(decision_space,by=c("years"="mcYear"), relationship="many-to-many") %>%
-    dplyr::rename("control"="u") %>%
     dplyr::mutate(next_state = -.data$control+.data$states+.data$hydroStorage) %>%
     rbind(possible_control)
 
@@ -78,7 +59,8 @@ build_all_possible_decisions <- function(Data_week,decision_space,
   df_SDP <- df_SDP %>%
     dplyr::filter(.data$next_state>=0, .data$next_state<=niveau_max)%>%
     dplyr::mutate(next_value = f_next_value(.data$next_state) + .data$overflow*overflow_cost) %>%
-    dplyr::select(-c("overflow"))
+    dplyr::select(c("years","states","next_state","control","hydroStorage","next_value"))
 
   return(df_SDP)
 }
+

R/Grid_Matrix.R

@@ -28,6 +28,7 @@
 #' @param correct_concavity Binary argument (default to \code{FALSE}). \code{TRUE} to correct concavity of Bellman values.
 #' @param reservoir_capacity Double. Reservoir capacity for the given area in MWh given by \code{get_reservoir_capacity()}.
 #' @param inflow Generated by \code{get_inflow()}
+#' @param plot_watervalues Binary.
 #' @returns
 #' \item{watervalues}{A \code{dplyr::tibble()} with multiple columns and detailed results.}
 #' \item{aggregated_results}{A \code{dplyr::tibble()} with multiple columns and summarized results.}
@@ -57,7 +58,8 @@ Grid_Matrix <- function(area,
                         force_final_level = F,
                         final_level = NULL,
                         penalty_final_level_low = NULL,
-                        penalty_final_level_high = NULL) {
+                        penalty_final_level_high = NULL,
+                        plot_watervalues = TRUE) {
 
 
   area = tolower(area)
@@ -104,8 +106,7 @@ Grid_Matrix <- function(area,
     inflow <- get_inflow(area=area, opts=opts,mcyears=mcyears)
   }
 
-  decision_space = reward_db$decision_space
-  reward_db = reward_db$reward
+  reward = reward_db$reward
 
   # Reservoir (rule curves)
   {
@@ -130,24 +131,12 @@ Grid_Matrix <- function(area,
     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
-  }
-
-  # 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"))
+    watervalues <- dplyr::right_join(x = watervalues, y = statesdt, 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 water values is the table containing (weeks,year,states,statesid,hydroStorage)
 
   #at this point we added the rewards for each weekly_amount
 
@@ -211,7 +200,7 @@ Grid_Matrix <- function(area,
         # 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),
+                        reward=dplyr::filter(reward,timeId==i),
                         E_max=E_max[i],
                         P_max=P_max[i],
                         mcyears = mcyears,
@@ -261,7 +250,8 @@ Grid_Matrix <- function(area,
       }
       # 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)
+                                               force_final_level=if(n_cycl==1){force_final_level}else{F},penalty_final_level_low,penalty_final_level_high,final_level=final_level,
+                                               plot_watervalues = plot_watervalues)
 
     }
 
@@ -283,7 +273,7 @@ Grid_Matrix <- function(area,
 
         temp <- Bellman(Data_week=temp,
                         next_week_values_l = next_week_values,
-                        decision_space=dplyr::filter(decision_space,week==i),
+                        reward=dplyr::filter(reward,timeId==i),
                         E_max=E_max[i],
                         P_max=P_max[i],
                         mcyears = mcyears,
@@ -331,7 +321,8 @@ Grid_Matrix <- function(area,
       }
       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)
+                                               penalty_final_level_low,penalty_final_level_high,final_level=final_level,
+                                               plot_watervalues = plot_watervalues)
 
 
       if(n_cycl>1){