Support for drawing availability profiles from historical data (#1372)

Author: anamileva

data_toolkit/project/availability/outages/create_availability_iteration_input_csvs.py

@@ -1,4 +1,5 @@
 # Copyright 2016-2025 Blue Marble Analytics LLC.
+# Copyright 2026 Sylvan Energy Analytics LLC.
 #
 # Licensed under the Apache License, Version 2.0 (the "License");
 # you may not use this file except in compliance with the License.
@@ -74,14 +75,24 @@ def parse_arguments(args):
 
     parser.add_argument("-db", "--database", default="../../io.db")
     parser.add_argument(
-        "-csv",
-        "--input_csv",
+        "-o_csv",
+        "--outage_params_input_csv",
         default=None,
         help="""Path to the unit availability  params CSV file to load into the 
         raw_data_unit_availability_params table in the database. If not 
         specified, data will be assumed to have been already loaded into the 
         database.""",
     )
+    parser.add_argument(
+        "-hist_csv",
+        "--historical_availability_csv",
+        default=None,
+        help="""Path to the historical availability data CSV file to load for the
+        historical_year outage model. Expected columns: year, month, day_of_month,
+        hour_of_day, unit, derate. Each unit can have multiple years of hourly
+        derate data. If not specified, model will not be able to use historical_year
+        outage_model.""",
+    )
     parser.add_argument("-stage", "--stage_id", default=1, help="Defaults to 1.")
     parser.add_argument("-n_iter", "--n_iterations")
     parser.add_argument(
@@ -188,6 +199,7 @@ def get_weighted_availability_adjustment(
     project_iteration_seed,
     max_integer_for_unit_outage_seeding,
     hyb_stor_seed_unit_increment,
+    historical_data=None,
 ):
     project_outage_adjustment = []
     project_hyb_stor_outage_adjustment = []
@@ -227,6 +239,8 @@ def get_weighted_availability_adjustment(
             mttr=unit_mttr,
             n_units=n_units,
             unit_seed=unit_seed,
+            historical_data=historical_data,
+            unit=unit,
         )
 
         # Get the project outage
@@ -246,6 +260,8 @@ def get_weighted_availability_adjustment(
                     if unit_seed is None
                     else unit_seed + hyb_stor_seed_unit_increment
                 ),
+                historical_data=historical_data,
+                unit=unit,
             )
             project_hyb_stor_outage_adjustment.append(
                 unit_outage_adjustment * unit_weight
@@ -277,8 +293,10 @@ def simulate_project_availability(
     study_year,
     filepath,
     print_ones,
+    historical_data=None,
 ):
 
+    # print(f"Simulating project {project}, iteration {iteration_n}")
     stage_tmp_dict = get_temporal_structure(study_year)
 
     # No stage simulation at this point; assume single stage
@@ -291,6 +309,7 @@ def simulate_project_availability(
         project_iteration_seed=project_iteration_seed,
         max_integer_for_unit_outage_seeding=max_integer_for_unit_outage_seeding,
         hyb_stor_seed_unit_increment=hyb_stor_seed_unit_increment,
+        historical_data=historical_data,
     )
 
     export_df = pd.DataFrame(
@@ -344,16 +363,23 @@ def simulate_unit_outages(
     unit_seed,
     dt=1,
     starting_outage_states=None,
+    historical_data=None,
+    unit=None,
 ):
     """
-    outage_model: ["Derate", "MC_independent", "MC_sequential"]
+    outage_model: ["Derate", "MC_independent", "MC_sequential", "historical_year"]
     FOR: numpy array with the length of the simulation window and the FOR as
         value; note this can vary by timepoint
     N_units: integer, number of units modeled
     starting_outage_states: array with the starting outage state (1/0) for
         each of the N units
     dt: outage timestep length
+    historical_data: dict with unit names as keys and DataFrames containing
+        historical availability data with columns: year, month, day_of_month,
+        hour_of_day, unit, derate (for historical_year model)
+    unit: unit name (for historical_year model)
     """
+    # print(f"Simulating unit outages for {unit}... Outage model: {outage_model}")
     # Seed the simulation if requested
     if unit_seed is not None:
         np.random.seed(unit_seed)
@@ -411,6 +437,42 @@ def simulate_unit_outages(
             availability[t, :] = avail_tmp
             avail_last = avail_tmp
 
+    elif outage_model == "historical_year":
+        # Sample a random year from historical data for this unit
+        if historical_data is None or unit not in historical_data:
+            raise ValueError(
+                f"Historical data not provided for unit {unit}. "
+                "Please provide historical_availability_csv."
+            )
+
+        unit_hist_data = historical_data[unit]
+
+        # Get list of unique years in the historical data
+        available_years = unit_hist_data["year"].unique()
+
+        if len(available_years) == 0:
+            raise ValueError(f"No historical years found for unit {unit}.")
+
+        # Randomly select a year
+        selected_year = np.random.choice(available_years)
+
+        # Get the derate values for the selected year
+        # Data is already sorted by year, month, day_of_month, hour_of_day
+        year_data = unit_hist_data[unit_hist_data["year"] == selected_year]
+        derate_values = year_data["value"].values
+
+        # Check if we have the right number of hours
+        if len(derate_values) != len(for_array):
+            raise ValueError(
+                f"Historical data for unit {unit}, year {selected_year} "
+                f"has {len(derate_values)} hours, but expected {len(for_array)}."
+            )
+
+        # Create availability array - replicate the same derate pattern for each
+        # of the n_units (this represents multiple identical units with the same
+        # historical outage pattern)
+        availability = np.outer(derate_values, np.ones(n_units))
+
     else:
         availability = np.ones([len(for_array), n_units])
 
@@ -478,13 +540,28 @@ def main(args=None):
     conn = connect_to_database(parsed_args.database)
 
     # ### Load data from CSV
-    if parsed_args.input_csv is not None:
+    if parsed_args.outage_params_input_csv is not None:
         read_and_import_csv(
             conn=conn,
-            f_path=parsed_args.input_csv,
+            f_path=parsed_args.outage_params_input_csv,
             table="raw_data_unit_availability_params",
         )
 
+    # Load historical availability data if provided
+    # Not loaded into the database for now
+    historical_data = None
+    if parsed_args.historical_availability_csv is not None:
+        # Read the CSV file directly into a DataFrame
+        hist_df = pd.read_csv(parsed_args.historical_availability_csv)
+
+        # Expected columns: year, month, day_of_month, hour_of_day, unit, derate
+        # Create an hour index for ordering (1-8760 or 1-8784 for leap years)
+        # Sort by temporal columns to ensure proper ordering
+        hist_df = hist_df.sort_values(["year", "month", "day_of_month", "hour_of_day"])
+
+        # Group by unit for unit-level access
+        historical_data = {unit: group for unit, group in hist_df.groupby("unit")}
+
     # Make out directory if it doesn't exist
     if not os.path.exists(parsed_args.output_directory):
         os.makedirs(parsed_args.output_directory)