pyRAVT.py

#!/usr/bin/env python
# coding: utf-8

# In[ ]:


# Script Name: pyRAVT_v0.1.0
# Author: Jo Ann Rañola
# Description: This Python script, pyRAVT, is an implementation of the Resource Adequacy Viewer Tool (RAVT), originally developed at EPRI (see https://ravt.epri.com/). pyRAVT calculates various Resource Adequacy (RA) metrics.
# Date Created: 2025-09-09


# In[105]:


import os
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
import numpy as np
import plotly.graph_objects as go
import plotly.io as pio
import datetime


# In[106]:


# Input the number of samples per weather year
no_samples = 150

# Input the number of weather years, start year, and end year
no_weather_years = 36
start_weather_year = 2012
end_weather_year = 2047

# Input the width and height of the figure
figure_width = 20
figure_height = 10
# figure_width = 12.5
# figure_height = 5.9

# Input the number of columns in LOLH/LOLD/EUE heatmap, scatterplot of event duration vs max USE, event duration/max USE/total USE histogram (all zones)
no_columns = 3
# no_columns = 1

# Input the CVaR alpha level (e.g., 0.05 for 5% CVaR, top 5% worst scenarios)
CVaR_alpha = 0.05
# CVaR_alpha = 1.00 # All scenarios


# In[107]:


# Set the directory
directory = os.getcwd()

# Create a directory named 'plot'
os.makedirs('output', exist_ok=True)

# Create a directory named 'plot' under 'output' folder
os.makedirs('output/plot', exist_ok=True)


# In[108]:


# Function to create datetime from Year, Month, Day, Hour columns
def create_datetime(row):
    return pd.to_datetime(f"{row['Year']}-{row['Month']}-{row['Day']} {row['Hour']}:00:00")

# Function to create date from Year, Month, Day columns
def create_date(row):
    return pd.to_datetime(f"{row['Year']}-{row['Month']}-{row['Day']}").date()

# Function to calculate the difference between two rows if they have the same 'Weather_Year' and 'Sample_Number' values
def calculate_difference(df):
    # Sort the DataFrame by 'Date', 'Weather_Year' and 'Sample_Number'
    df = df.sort_values(by=['Weather_Year', 'Sample_Number', 'Scenario_Name', 'Zone_Name', 'Datetime']) 
    # Calculate the difference using groupby and diff
    df['Difference'] = df.groupby(['Weather_Year', 'Sample_Number', 'Scenario_Name', 'Zone_Name'])['Datetime'].diff()    
    return df

# Function to calculate the difference between two rows if they have the same 'Weather_Year' and 'Sample_Number' values - all zones
def calculate_difference_all_zones(df):
    # Sort the DataFrame by 'Date', 'Weather_Year' and 'Sample_Number'
    df = df.sort_values(by=['Weather_Year', 'Sample_Number', 'Scenario_Name', 'Datetime']) 
    # Calculate the difference using groupby and diff
    df['Difference'] = df.groupby(['Weather_Year', 'Sample_Number', 'Scenario_Name'])['Datetime'].diff()    
    return df

# Function to add 'Event' column based on the specified conditions
def add_event_column(df):
    df['Event'] = None
    for i in range(len(df)):
        if i == 0:
            df.at[i, 'Event'] = 1
        else:
            if df.at[i, 'Difference'] == 3600:
                df.at[i, 'Event'] = df.at[i-1, 'Event']
            else:
                df.at[i, 'Event'] = df.at[i-1, 'Event'] + 1
    return df

# Function to create the heatmap of LOLH by hour and month
def heatmap_LOLH(data, **kwargs):
    pivot_data = data.pivot_table(values='LOLH', index='Hour', columns='Month', aggfunc='sum')
    pivot_data = pivot_data.reindex(index=range(24), columns=range(1, 13), fill_value=0)
    month_labels = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec']
    ax = sns.heatmap(pivot_data, annot=False, fmt=".1f", cmap="coolwarm", vmin=0, **kwargs)
    ax.set_xticklabels(month_labels)  # Set month names
    ax.invert_yaxis()  # Invert the y-axis
    for ax in g.axes.flat:
        for label in ax.get_xticklabels():
            label.set_rotation(90)

# Function to create the heatmap of LOLD by weather year and month
def heatmap_LOLD(data, **kwargs):
    pivot_data = data.pivot_table(values='LOLD', index='Weather_Year', columns='Month', aggfunc='sum')
    pivot_data = pivot_data.reindex(columns=range(1, 13), fill_value=0)
    month_labels = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec']
    ax = sns.heatmap(pivot_data, annot=False, fmt=".1f", cmap="coolwarm", vmin=0, **kwargs)
    ax.set_xticklabels(month_labels)  # Set month names
    for ax in g.axes.flat:
        for label in ax.get_xticklabels():
            label.set_rotation(90)

# Function to create the heatmap of EUE by hour and month
def heatmap_EUE(data, **kwargs):
    pivot_data = data.pivot_table(values='EUE', index='Hour', columns='Month', aggfunc='sum')
    pivot_data = pivot_data.reindex(index=range(24), columns=range(1, 13), fill_value=0)
    month_labels = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec']
    ax = sns.heatmap(pivot_data, annot=False, fmt=".1f", cmap="coolwarm", vmin=0, **kwargs)
    ax.set_xticklabels(month_labels)  # Set month names
    ax.invert_yaxis()  # Invert the y-axis
    for ax in g.axes.flat:
        for label in ax.get_xticklabels():
            label.set_rotation(90)

# Define the LOLH stacked bar plot function
def stacked_bar_plot_LOLH(data, **kwargs):
    pivot_data = data.pivot_table(values='LOLH', index='Month', columns='Weather_Year', aggfunc='sum').fillna(0)
    cmap = plt.get_cmap('viridis')  # Choose a colormap
    colors = [cmap(i / len(pivot_data.columns)) for i in range(len(pivot_data.columns))]
    ax = pivot_data.plot(kind='bar', stacked=True, ax=plt.gca(), color=colors)
    months_present = pivot_data.index.tolist()
    month_labels_full = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec']
    month_labels = [month_labels_full[m-1] for m in months_present]
    ax.set_xticklabels(month_labels, rotation=0)
    for ax in g.axes.flat:
        for label in ax.get_xticklabels():
            label.set_rotation(90)

# Define the LOLD stacked bar plot function
def stacked_bar_plot_LOLD(data, **kwargs):
    pivot_data = data.pivot_table(values='LOLD', index='Month', columns='Weather_Year', aggfunc='sum').fillna(0)
    cmap = plt.get_cmap('viridis')  # Choose a colormap
    colors = [cmap(i / len(pivot_data.columns)) for i in range(len(pivot_data.columns))]
    ax = pivot_data.plot(kind='bar', stacked=True, ax=plt.gca(), color=colors)
    months_present = pivot_data.index.tolist()
    month_labels_full = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec']
    month_labels = [month_labels_full[m-1] for m in months_present]
    ax.set_xticklabels(month_labels, rotation=0)
    for ax in g.axes.flat:
        for label in ax.get_xticklabels():
            label.set_rotation(90)

# Define the EUE stacked bar plot function
def stacked_bar_plot_EUE(data, **kwargs):
    pivot_data = data.pivot_table(values='EUE', index='Month', columns='Weather_Year', aggfunc='sum').fillna(0)
    cmap = plt.get_cmap('viridis')  # Choose a colormap
    colors = [cmap(i / len(pivot_data.columns)) for i in range(len(pivot_data.columns))]
    ax = pivot_data.plot(kind='bar', stacked=True, ax=plt.gca(), color=colors)
    months_present = pivot_data.index.tolist()
    month_labels_full = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec']
    month_labels = [month_labels_full[m-1] for m in months_present]
    ax.set_xticklabels(month_labels, rotation=0)
    for ax in g.axes.flat:
        for label in ax.get_xticklabels():
            label.set_rotation(90)

# Define the LOLH stacked bar plot function by weather year
def stacked_bar_plot_LOLH_wy(data, **kwargs):
    pivot_data = data.pivot_table(values='LOLH', index='Weather_Year', columns='Month', aggfunc='sum').fillna(0)
    cmap = plt.get_cmap('viridis')  # Choose a colormap
    colors = [cmap(i / len(pivot_data.columns)) for i in range(len(pivot_data.columns))]
    ax = pivot_data.plot(kind='bar', stacked=True, ax=plt.gca(), color=colors)
    handles, labels = g.axes[0, 0].get_legend_handles_labels()
    month_labels_full = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec']
    labels = [month_labels_full[int(l)-1] if l.isdigit() and 1 <= int(l) <= 12 else l for l in labels]
    g.fig.legend(handles, labels, title='Month', loc='upper right')
    for ax in g.axes.flat:
        for label in ax.get_xticklabels():
            label.set_rotation(90)

# Define the LOLD stacked bar plot function by weather year
def stacked_bar_plot_LOLD_wy(data, **kwargs):
    pivot_data = data.pivot_table(values='LOLD', index='Weather_Year', columns='Month', aggfunc='sum').fillna(0)
    cmap = plt.get_cmap('viridis')  # Choose a colormap
    colors = [cmap(i / len(pivot_data.columns)) for i in range(len(pivot_data.columns))]
    ax = pivot_data.plot(kind='bar', stacked=True, ax=plt.gca(), color=colors)
    handles, labels = g.axes[0, 0].get_legend_handles_labels()
    month_labels_full = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec']
    labels = [month_labels_full[int(l)-1] if l.isdigit() and 1 <= int(l) <= 12 else l for l in labels]
    g.fig.legend(handles, labels, title='Month', loc='upper right')
    for ax in g.axes.flat:
        for label in ax.get_xticklabels():
            label.set_rotation(90)

# Define the EUE stacked bar plot function by weather year
def stacked_bar_plot_EUE_wy(data, **kwargs):
    pivot_data = data.pivot_table(values='EUE', index='Weather_Year', columns='Month', aggfunc='sum').fillna(0)
    cmap = plt.get_cmap('viridis')  # Choose a colormap
    colors = [cmap(i / len(pivot_data.columns)) for i in range(len(pivot_data.columns))]
    ax = pivot_data.plot(kind='bar', stacked=True, ax=plt.gca(), color=colors)
    handles, labels = g.axes[0, 0].get_legend_handles_labels()
    month_labels_full = ['Jan', 'Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov', 'Dec']
    labels = [month_labels_full[int(l)-1] if l.isdigit() and 1 <= int(l) <= 12 else l for l in labels]
    g.fig.legend(handles, labels, title='Month', loc='upper right')
    for ax in g.axes.flat:
        for label in ax.get_xticklabels():
            label.set_rotation(90)

# Function to create the heatmap of LOLEv by weather year and sample number
def heatmap_LOLEv(data, **kwargs):
    pivot_data = data.pivot_table(values='LOLEv', index='Weather_Year', columns='Sample_Number', aggfunc='sum')
    ax = sns.heatmap(pivot_data, annot=False, fmt=".1f", cmap="coolwarm", vmin=0, **kwargs)
    for ax in g.axes.flat:
        for label in ax.get_xticklabels():
            label.set_rotation(90)

# Function to create area plots
def area_plot(data, **kwargs):
    plt.fill_between(data['Hour'], data['Unserved_Energy_MWh'], alpha=0.5)
    plt.xlim(0, 23)  # Set x-axis limits

# Calculate rolling average of LOLH or EUE by Scenario_Name with window size based on number of rows
def get_window_size(df, scenario_col='Scenario_Name'):
    min_window = 1
    window_sizes = {}
    for scenario in df[scenario_col].unique():
        n_rows = df[df[scenario_col] == scenario].shape[0]
        window_sizes[scenario] = max(int(n_rows), min_window)
    return window_sizes


# In[109]:


# Read the CSV file into a DataFrame
data_use = pd.read_csv(directory + '/input/unserved_energy.csv')

# Calculate the hourly weight
no_simulations = no_samples * no_weather_years
use_wgt = 1 / no_simulations

# Add the hourly weight
data_use['LOLH'] = use_wgt


# In[110]:


# Calculate the LOLH
LOLH = data_use.groupby(['Year', 'Scenario_Name', 'Zone_Name']).agg({
    'LOLH': 'sum'
}).reset_index()


# In[111]:


# Calculate the EUE
data_use['EUE'] = data_use['LOLH'] * data_use['Unserved_Energy_MWh']
EUE = data_use.groupby(['Year', 'Scenario_Name', 'Zone_Name']).agg({
    'EUE': 'sum'
}).reset_index()


# In[112]:


# Calculate the LOLD
data_use['Date'] = data_use.apply(create_date, axis=1)
data_use_day = data_use.groupby(['Year', 'Scenario_Name', 'Weather_Year', 'Sample_Number', 'Zone_Name', 'Date']).agg({
    'Unserved_Energy_MWh': 'sum'
}).reset_index()
data_use_day['Date'] = pd.to_datetime(data_use_day['Date'])
data_use_day['Month'] = data_use_day['Date'].dt.month
data_use_day['LOLD'] = use_wgt
LOLD = data_use_day.groupby(['Year', 'Scenario_Name', 'Zone_Name']).agg({
    'LOLD': 'sum'
}).reset_index()


# In[113]:


# Calculate the LOLEv
data_use['Datetime'] = data_use.apply(create_datetime, axis=1)
data_use = calculate_difference(data_use)
data_use['Difference'] = pd.to_timedelta(data_use['Difference'])
data_use['Difference'] = data_use['Difference'].dt.total_seconds()
data_use = add_event_column(data_use)
data_use_event = data_use.groupby(['Year', 'Scenario_Name', 'Weather_Year', 'Sample_Number', 'Zone_Name', 'Event']).agg({
    'Unserved_Energy_MWh': 'sum'
}).reset_index()
data_use_event['LOLEv'] = use_wgt
LOLEv = data_use_event.groupby(['Year', 'Scenario_Name', 'Zone_Name']).agg({
    'LOLEv': 'sum'
}).reset_index()


# In[114]:


# Merge the all the metrics dataframes on 'Scenario_Name' and 'Zone_Name' columns
merged_df = pd.merge(LOLH, LOLD, on=['Year', 'Scenario_Name', 'Zone_Name'])
merged_df['LOLP'] = merged_df['LOLD']/3.65
merged_df = pd.merge(merged_df, LOLEv, on=['Year', 'Scenario_Name', 'Zone_Name'])
merged_df = pd.merge(merged_df, EUE, on=['Year', 'Scenario_Name', 'Zone_Name'])

# Write to CSV file
merged_df.to_csv(directory + '/output/metrics_by_zone.csv', index=False)


# In[115]:


# Calculate metrics (all zones combined)
data_use_all_zones = data_use.groupby(['Scenario_Name', 'Weather_Year', 'Sample_Number', 'Year', 'Month', 'Day', 'Hour', 'Date', 'Datetime']).agg({
    'Unserved_Energy_MWh': 'sum'
}).reset_index()

data_use_all_zones['LOLH'] = use_wgt

data_use_all_zones['EUE'] = data_use_all_zones['LOLH'] * data_use_all_zones['Unserved_Energy_MWh']

LOLH_all_zones = data_use_all_zones.groupby(['Year', 'Scenario_Name']).agg({
    'LOLH': 'sum'
}).reset_index()

EUE_all_zones = data_use_all_zones.groupby(['Year', 'Scenario_Name']).agg({
    'EUE': 'sum'
}).reset_index()

data_use_day_all_zones = data_use_all_zones.groupby(['Year', 'Scenario_Name', 'Weather_Year', 'Sample_Number', 'Date']).agg({
    'Unserved_Energy_MWh': 'sum'
}).reset_index()
data_use_day_all_zones['LOLD'] = use_wgt

LOLD_all_zones = data_use_day_all_zones.groupby(['Year', 'Scenario_Name']).agg({
    'LOLD': 'sum'
}).reset_index()

data_use_all_zones = calculate_difference_all_zones(data_use_all_zones)
data_use_all_zones['Difference'] = pd.to_timedelta(data_use_all_zones['Difference'])
data_use_all_zones['Difference'] = data_use_all_zones['Difference'].dt.total_seconds()
data_use_all_zones = add_event_column(data_use_all_zones)
data_use_event_all_zones = data_use_all_zones.groupby(['Year', 'Scenario_Name', 'Weather_Year', 'Sample_Number', 'Event']).agg({
    'Unserved_Energy_MWh': 'sum'
}).reset_index()
data_use_event_all_zones['LOLEv'] = use_wgt

LOLEv_all_zones = data_use_event_all_zones.groupby(['Year', 'Scenario_Name']).agg({
    'LOLEv': 'sum'
}).reset_index()

# Merge the all the metrics dataframes on 'Scenario_Name' and 'Zone_Name' columns
merged_df_all_zones = pd.merge(LOLH_all_zones, LOLD_all_zones, on=['Year', 'Scenario_Name'])
merged_df_all_zones['LOLP'] = merged_df_all_zones['LOLD']/3.65
merged_df_all_zones = pd.merge(merged_df_all_zones, LOLEv_all_zones, on=['Year', 'Scenario_Name'])
merged_df_all_zones = pd.merge(merged_df_all_zones, EUE_all_zones, on=['Year', 'Scenario_Name'])

# Write to CSV file
merged_df_all_zones.to_csv(directory + '/output/metrics_all_zones.csv', index=False)


# In[116]:


# Calculate the individual events characteristics (all zones combined)
data_use_all_zones_x = data_use_all_zones.groupby(['Year', 'Scenario_Name', 'Weather_Year', 'Sample_Number', 'Event']).agg({
    'Unserved_Energy_MWh': [('USE_MWh', 'sum'), ('Max_USE_MWh', 'max')],
    'Event': [('Event_Duration', 'count'), ('Event_No', 'mean')]
}).reset_index()
# data_use_all_zones_x.to_csv(directory + '/output/data_use_all_zones_x.csv', index=False)


# In[117]:


# Calculate CVaR

# Number of simulations to consider for CVaR
CVaR_no_simulations = int(round(no_simulations * CVaR_alpha))

# Create the all simulations dataframe
scenario_names = data_use_all_zones['Scenario_Name'].unique()
weather_years = np.arange(start_weather_year, end_weather_year + 1)  
sample_numbers = np.arange(1, no_samples + 1)

rows = []
for scenario in scenario_names:
    sim_num = 1
    for wy in weather_years:
        for sn in sample_numbers:
            rows.append([scenario, wy, sn, sim_num])
            sim_num += 1

df_sim = pd.DataFrame(rows, columns=['Scenario_Name', 'Weather_Year', 'Sample_Number', 'Simulation_Number'])

# Calculate LOLH and EUE by simulation number
data_use_all_zones_sim = data_use_all_zones.groupby(['Scenario_Name', 'Weather_Year', 'Sample_Number']).agg({
    'LOLH': 'sum',
    'EUE': 'sum'
}).reset_index()

data_use_all_zones_sim['LOLH'] = data_use_all_zones_sim['LOLH'] * no_simulations
data_use_all_zones_sim['EUE'] = data_use_all_zones_sim['EUE'] * no_simulations

# Merge dataframes 
CVaR_merged_df = pd.merge(df_sim, data_use_all_zones_sim, how = 'left', indicator = False)
CVaR_merged_df['LOLH'] = CVaR_merged_df['LOLH'].fillna(0)
CVaR_merged_df['EUE'] = CVaR_merged_df['EUE'].fillna(0)

# Sort LOLH from highest to lowest and group by Scenario_Name
sorted_df_LOLH = (
    CVaR_merged_df
    .sort_values(['Scenario_Name', 'LOLH'], ascending=[True, False])
    .copy()
)

# Replace Simulation_Number so that each Scenario_Name starts at 1 and ends at no_simulations
sorted_df_LOLH['Simulation_Number'] = (
    sorted_df_LOLH.groupby('Scenario_Name').cumcount() + 1
)

# Sort EUE from highest to lowest and group by Scenario_Name
sorted_df_EUE = (
    CVaR_merged_df
    .sort_values(['Scenario_Name', 'EUE'], ascending=[True, False])
    .copy()
)

# Replace Simulation_Number so that each Scenario_Name starts at 1 and ends at no_simulations
sorted_df_EUE['Simulation_Number'] = (
    sorted_df_EUE.groupby('Scenario_Name').cumcount() + 1
)

# Calculate rolling average of LOLH by Scenario_Name with window size based on number of rows
window_sizes = get_window_size(sorted_df_LOLH)

sorted_df_LOLH['LOLH_rolling_avg'] = (
    sorted_df_LOLH
    .groupby('Scenario_Name')['LOLH']
    .apply(lambda x: x.rolling(window=window_sizes[x.name], min_periods=1).mean())
    .reset_index(level=0, drop=True)
)

# Calculate rolling average of EUE by Scenario_Name with window size based on number of rows
window_sizes = get_window_size(sorted_df_EUE)

sorted_df_EUE['EUE_rolling_avg'] = (
    sorted_df_EUE
    .groupby('Scenario_Name')['EUE']
    .apply(lambda x: x.rolling(window=window_sizes[x.name], min_periods=1).mean())
    .reset_index(level=0, drop=True)
)

# Remove rows with Simulation_Number greater than CVaR_no_simulations
CVaR_sorted_df_LOLH = sorted_df_LOLH[sorted_df_LOLH['Simulation_Number'] <= CVaR_no_simulations].copy()
CVaR_sorted_df_EUE = sorted_df_EUE[sorted_df_EUE['Simulation_Number'] <= CVaR_no_simulations].copy()

# Filter sorted_df_LOLH where Simulation_Number equals CVaR_no_simulations
CVaR_LOLH = sorted_df_LOLH[sorted_df_LOLH['Simulation_Number'] == CVaR_no_simulations]
CVaR_LOLH = CVaR_LOLH.drop(columns=['LOLH', 'EUE'])

# Filter sorted_df_EUE where Simulation_Number equals CVaR_no_simulations
CVaR_EUE = sorted_df_EUE[sorted_df_EUE['Simulation_Number'] == CVaR_no_simulations]
CVaR_EUE = CVaR_EUE.drop(columns=['LOLH', 'EUE'])

# Write to CSV file
CVaR_LOLH.to_csv(directory + '/output/CVaR_LOLH.csv', index=False)
CVaR_EUE.to_csv(directory + '/output/CVaR_EUE.csv', index=False)


# In[118]:


# Plot scatterplot of event duration vs max USE with bubble size as total USE
df_scatter = data_use_all_zones_x.copy()
df_scatter = df_scatter.droplevel(1, axis=1) if isinstance(df_scatter.columns, pd.MultiIndex) else df_scatter
df_scatter['Event_Duration'] = data_use_all_zones_x['Event']['Event_Duration']
df_scatter['Max_USE_MWh'] = data_use_all_zones_x['Unserved_Energy_MWh']['Max_USE_MWh']
df_scatter['USE_MWh'] = data_use_all_zones_x['Unserved_Energy_MWh']['USE_MWh']
df_scatter['Scenario_Name'] = data_use_all_zones_x['Scenario_Name']
g = sns.FacetGrid(df_scatter, col="Scenario_Name", col_wrap = no_columns, sharex=True, sharey=True, height=5)
g.map_dataframe(
    sns.scatterplot,
    x='Event_Duration',
    y='Max_USE_MWh',
    size='USE_MWh',
    sizes=(20, 400),
    alpha=0.7,
    legend=False
)
g.set_axis_labels('Event Duration (hours)', 'Event Maximum Unserved Energy (MWh)')
g.set_titles(col_template="{col_name}")
g.fig.suptitle('Event Duration vs Max USE (Bubble Size = Total USE)')
g.fig.tight_layout(rect=[0, 0, 0.85, 0.95])
g.fig.set_size_inches(figure_width, figure_height)
plt.savefig(f"{directory}/output/plot/LOLEv_Duration_vs_MaxUSE_Bubble.svg")
# plt.close()


# In[119]:


# Plot histogram of event duration
df_hist = data_use_all_zones_x.copy()
df_hist = df_hist.droplevel(1, axis=1) if isinstance(df_hist.columns, pd.MultiIndex) else df_hist
df_hist['Event_Duration'] = data_use_all_zones_x['Event']['Event_Duration']
df_hist['Scenario_Name'] = data_use_all_zones_x['Scenario_Name']
g = sns.FacetGrid(df_hist, col="Scenario_Name", col_wrap = no_columns, sharex=True, sharey=True)
g.map(plt.hist, "Event_Duration", bins=30, color='skyblue', edgecolor='black')
g.set_axis_labels('Event Duration (hours)', 'Frequency')
g.set_titles(col_template="{col_name}")
g.fig.suptitle('Histogram of Event Duration')
g.fig.tight_layout(rect=[0, 0, 0.85, 1])
g.fig.set_size_inches(figure_width, figure_height)
plt.savefig(f"{directory}/output/plot/LOLEv_Duration.svg")
# plt.close()


# In[120]:


# Plot histogram of event max USE
df_hist = data_use_all_zones_x.copy()
df_hist = df_hist.droplevel(1, axis=1) if isinstance(df_hist.columns, pd.MultiIndex) else df_hist
df_hist['Event_Max_USE'] = data_use_all_zones_x['Unserved_Energy_MWh']['Max_USE_MWh']
df_hist['Scenario_Name'] = data_use_all_zones_x['Scenario_Name']
g = sns.FacetGrid(df_hist, col="Scenario_Name", col_wrap = no_columns, sharex=True, sharey=True)
g.map(plt.hist, "Event_Max_USE", bins=30, color='skyblue', edgecolor='black')
g.set_axis_labels('Event Maximum Unserved Energy (MWh)', 'Frequency')
g.set_titles(col_template="{col_name}")
g.fig.suptitle('Histogram of Event Maximum Unserved Energy')
g.fig.tight_layout(rect=[0, 0, 0.85, 1])
g.fig.set_size_inches(figure_width, figure_height)
plt.savefig(f"{directory}/output/plot/LOLEv_Max_USE.svg")
# plt.close()


# In[121]:


# Plot histogram of event total USE
df_hist = data_use_all_zones_x.copy()
df_hist = df_hist.droplevel(1, axis=1) if isinstance(df_hist.columns, pd.MultiIndex) else df_hist
df_hist['Event_Total_USE'] = data_use_all_zones_x['Unserved_Energy_MWh']['USE_MWh']
df_hist['Scenario_Name'] = data_use_all_zones_x['Scenario_Name']
g = sns.FacetGrid(df_hist, col="Scenario_Name", col_wrap = no_columns, sharex=True, sharey=True)
g.map(plt.hist, "Event_Total_USE", bins=30, color='skyblue', edgecolor='black')
g.set_axis_labels('Event Total Unserved Energy (MWh)', 'Frequency')
g.set_titles(col_template="{col_name}")
g.fig.suptitle('Histogram of Event Total Unserved Energy')
g.fig.tight_layout(rect=[0, 0, 0.85, 1])
g.fig.set_size_inches(figure_width, figure_height)
plt.savefig(f"{directory}/output/plot/LOLEv_USE.svg")
# plt.close()


# In[122]:


# Plot LOLH heatmap (all zones)
g = sns.FacetGrid(data_use, col="Scenario_Name", col_wrap = no_columns)
g.map_dataframe(heatmap_LOLH)
g.set_titles(row_template="{row_name}", col_template="{col_name}")
g.set_axis_labels("Month", "Hour")
g.fig.suptitle("LOLH Heatmap (All Zones)")
g.fig.tight_layout(rect=[0, 0, 0.85, 1])
g.fig.set_size_inches(figure_width, figure_height)
plt.savefig(directory + '/output/plot/LOLH_Heatmap.svg')
# plt.close()


# In[123]:


# Plot LOLD heatmap (all zones)
g = sns.FacetGrid(data_use_day, col="Scenario_Name", col_wrap = no_columns)
g.map_dataframe(heatmap_LOLD)
g.set_titles(row_template="{row_name}", col_template="{col_name}")
g.set_axis_labels("Month", "Weather Year")
g.fig.suptitle("LOLD Heatmap (All Zones)")
g.fig.tight_layout(rect=[0, 0, 0.85, 1])  
g.fig.set_size_inches(figure_width, figure_height)
plt.savefig(directory + '/output/plot/LOLD_Heatmap.svg')
# plt.close()


# In[124]:


# Plot EUE heatmap (all zones)
g = sns.FacetGrid(data_use, col="Scenario_Name", col_wrap = no_columns)
g.map_dataframe(heatmap_EUE)
g.set_titles(row_template="{row_name}", col_template="{col_name}")
g.set_axis_labels("Month", "Hour")
g.fig.suptitle("EUE Heatmap (All Zones)")
g.fig.tight_layout(rect=[0, 0, 0.85, 1])
g.fig.set_size_inches(figure_width, figure_height)
plt.savefig(directory + '/output/plot/EUE_Heatmap.svg')
# plt.close()


# In[125]:


# Plot LOLH heatmap by zone
g = sns.FacetGrid(data_use, col="Scenario_Name", row="Zone_Name")
g.map_dataframe(heatmap_LOLH)
g.set_titles(row_template="{row_name}", col_template="{col_name}")
g.set_axis_labels("Month", "Hour")
g.fig.suptitle("LOLH Heatmap by Zone")
g.fig.tight_layout(rect=[0, 0, 0.85, 1])
g.fig.set_size_inches(figure_width, figure_height)
plt.savefig(directory + '/output/plot/LOLH_Heatmap_Zone.svg')
# plt.close()


# In[126]:


# Plot LOLD heatmap by zone
g = sns.FacetGrid(data_use_day, col="Scenario_Name", row="Zone_Name")
g.map_dataframe(heatmap_LOLD)
g.set_titles(row_template="{row_name}", col_template="{col_name}")
g.set_axis_labels("Month", "Weather Year")
g.fig.suptitle("LOLD Heatmap")
g.fig.tight_layout(rect=[0, 0, 0.85, 1])  
g.fig.set_size_inches(figure_width, figure_height)
plt.savefig(directory + '/output/plot/LOLD_Heatmap_Zone.svg')
# plt.close()


# In[127]:


# Plot EUE heatmap by zone
g = sns.FacetGrid(data_use, col="Scenario_Name", row="Zone_Name")
g.map_dataframe(heatmap_EUE)
g.set_titles(row_template="{row_name}", col_template="{col_name}")
g.set_axis_labels("Month", "Hour")
g.fig.suptitle("EUE Heatmap by Zone")
g.fig.tight_layout(rect=[0, 0, 0.85, 1])
g.fig.set_size_inches(figure_width, figure_height)
plt.savefig(directory + '/output/plot/EUE_Heatmap_Zone.svg')
# plt.close()


# In[128]:


# Plot LOLH by Month
g = sns.FacetGrid(data_use, col='Scenario_Name', margin_titles=True)
g.map_dataframe(stacked_bar_plot_LOLH)
g.set_axis_labels('Month', 'LOLH (hours/year)')
g.set_titles(col_template='{col_name}')
g.fig.suptitle('LOLH by Month')
handles, labels = g.axes[0, 0].get_legend_handles_labels()
g.fig.legend(handles, labels, title='Weather Year', loc='upper right')
g.fig.tight_layout(rect=[0, 0, 0.85, 1])
g.fig.set_size_inches(figure_width, figure_height)
plt.savefig(directory + '/output/plot/LOLH_Month.svg')
# plt.close()


# In[129]:


# Plot LOLD by Month
g = sns.FacetGrid(data_use_day, col='Scenario_Name', margin_titles=True)
g.map_dataframe(stacked_bar_plot_LOLD)
g.set_axis_labels('Month', 'LOLD (days/year)')
g.set_titles(col_template='{col_name}')
g.fig.suptitle('LOLD by Month')
handles, labels = g.axes[0, 0].get_legend_handles_labels()
g.fig.legend(handles, labels, title='Weather Year', loc='upper right')
g.fig.tight_layout(rect=[0, 0, 0.85, 1]) 
g.fig.set_size_inches(figure_width, figure_height)
plt.savefig(directory + '/output/plot/LOLD_Month.svg')
# plt.close()


# In[130]:


# Plot EUE by Month
g = sns.FacetGrid(data_use, col='Scenario_Name', margin_titles=True)
g.map_dataframe(stacked_bar_plot_EUE)
g.set_axis_labels('Month', 'EUE (MWh/year)')
g.set_titles(col_template='{col_name}')
g.fig.suptitle('EUE by Month')
handles, labels = g.axes[0, 0].get_legend_handles_labels()
g.fig.legend(handles, labels, title='Weather Year', loc='upper right')
g.fig.tight_layout(rect=[0, 0, 0.85, 1])
g.fig.set_size_inches(figure_width, figure_height)
plt.savefig(directory + '/output/plot/EUE_Month.svg')
# plt.close()


# In[131]:


# Plot LOLH by Weather Year
g = sns.FacetGrid(data_use, col='Scenario_Name', margin_titles=True)
g.map_dataframe(stacked_bar_plot_LOLH_wy)
g.set_axis_labels('Weather Year', 'LOLH (hours/year)')
g.set_titles(col_template='{col_name}')
g.fig.suptitle('LOLH by Weather Year')
# Adjust layout to prevent cutting off labels
g.fig.tight_layout(rect=[0, 0, 0.85, 1]) 
g.fig.set_size_inches(figure_width, figure_height)
plt.savefig(directory + '/output/plot/LOLH_WY.svg')
# plt.close()


# In[132]:


# Plot LOLD by Weather Year
g = sns.FacetGrid(data_use_day, col='Scenario_Name', margin_titles=True)
g.map_dataframe(stacked_bar_plot_LOLD_wy)
g.set_axis_labels('Weather Year', 'LOLD (days/year)')
g.set_titles(col_template='{col_name}')
g.fig.suptitle('LOLD by Weather Year')
g.fig.tight_layout(rect=[0, 0, 0.85, 1])
g.fig.set_size_inches(figure_width, figure_height)
plt.savefig(directory + '/output/plot/LOLD_WY.svg')
# plt.close()


# In[133]:


# Plot EUE by Weather Year
g = sns.FacetGrid(data_use, col='Scenario_Name', margin_titles=True)
g.map_dataframe(stacked_bar_plot_EUE_wy)
g.set_axis_labels('Weather Year', 'EUE (MWh/year)')
g.set_titles(col_template='{col_name}')
g.fig.suptitle('EUE by Weather Year')
g.fig.tight_layout(rect=[0, 0, 0.85, 1])
g.fig.set_size_inches(figure_width, figure_height)
plt.savefig(directory + '/output/plot/EUE_WY.svg')
# plt.close()


# In[134]:


# Plot LOLH_rolling_avg, faceted by Scenario_Name (all simulations)
g = sns.FacetGrid(sorted_df_LOLH, col="Scenario_Name", col_wrap=no_columns, sharex=False, sharey=True, height=5)
g.map_dataframe(sns.lineplot, x="Simulation_Number", y="LOLH_rolling_avg")

# Add label to the last value in each facet
for ax, scenario in zip(g.axes.flat, sorted_df_LOLH['Scenario_Name'].unique()):
    df_scenario = sorted_df_LOLH[sorted_df_LOLH['Scenario_Name'] == scenario]
    last_row = df_scenario[df_scenario['Simulation_Number'] == df_scenario['Simulation_Number'].max()]
    if not last_row.empty:
        x = last_row['Simulation_Number'].values[0]
        y = last_row['LOLH_rolling_avg'].values[0]
        ax.text(x, y, f"{y:.2f}", color="red", fontsize=12, ha="right", va="bottom")

g.set_axis_labels("Simulation", "LOLH (hours/year)")
g.set_titles(col_template="{col_name}")
# g.fig.suptitle("LOLH Rolling Average by Scenario", y=1.05)
g.fig.suptitle(f"LOLH Rolling Average by Scenario (All Simulations)")
g.fig.tight_layout(rect=[0, 0, 0.85, 0.95])
g.fig.set_size_inches(figure_width, figure_height)
plt.savefig(f"{directory}/output/plot/Rolling_LOLH_by_Scenario.svg")
# plt.close()


# In[ ]:


# Plot LOLH_rolling_avg, faceted by Scenario_Name (CVaR only)
g = sns.FacetGrid(CVaR_sorted_df_LOLH, col="Scenario_Name", col_wrap=no_columns, sharex=False, sharey=True, height=5)
g.map_dataframe(sns.lineplot, x="Simulation_Number", y="LOLH_rolling_avg")

# Add label to the last value in each facet
for ax, scenario in zip(g.axes.flat, CVaR_sorted_df_LOLH['Scenario_Name'].unique()):
    df_scenario = CVaR_sorted_df_LOLH[CVaR_sorted_df_LOLH['Scenario_Name'] == scenario]
    last_row = df_scenario[df_scenario['Simulation_Number'] == df_scenario['Simulation_Number'].max()]
    if not last_row.empty:
        x = last_row['Simulation_Number'].values[0]
        y = last_row['LOLH_rolling_avg'].values[0]
        ax.text(x, y, f"{y:.2f}", color="red", fontsize=12, ha="right", va="bottom")

g.set_axis_labels("Simulation", "LOLH (hours/year)")
g.set_titles(col_template="{col_name}")
# g.fig.suptitle("LOLH Rolling Average by Scenario", y=1.05)
g.fig.suptitle(f"LOLH Rolling Average by Scenario | {CVaR_alpha} CVaR")
g.fig.tight_layout(rect=[0, 0, 0.85, 0.95])
g.fig.set_size_inches(figure_width, figure_height)
plt.savefig(f"{directory}/output/plot/Rolling_LOLH_by_Scenario_CVaR.svg")
# plt.close()


# In[4]:


# Plot EUE_rolling_avg, faceted by Scenario_Name (all simulations)
g = sns.FacetGrid(sorted_df_EUE, col="Scenario_Name", col_wrap=no_columns, sharex=False, sharey=True, height=5)
g.map_dataframe(sns.lineplot, x="Simulation_Number", y="EUE_rolling_avg")

# Add label to the last value in each facet
for ax, scenario in zip(g.axes.flat, sorted_df_EUE['Scenario_Name'].unique()):
    df_scenario = sorted_df_EUE[sorted_df_EUE['Scenario_Name'] == scenario]
    last_row = df_scenario[df_scenario['Simulation_Number'] == df_scenario['Simulation_Number'].max()]
    if not last_row.empty:
        x = last_row['Simulation_Number'].values[0]
        y = last_row['EUE_rolling_avg'].values[0]
        ax.text(x, y, f"{y:.2f}", color="red", fontsize=12, ha="right", va="bottom")

g.set_axis_labels("Simulation", "EUE (MWh/year)")
g.set_titles(col_template="{col_name}")
# g.fig.suptitle("LOLH Rolling Average by Scenario", y=1.05)
g.fig.suptitle(f"EUE Rolling Average by Scenario (All Simulations)")
g.fig.tight_layout(rect=[0, 0, 0.85, 0.95])
g.fig.set_size_inches(figure_width, figure_height)
plt.savefig(f"{directory}/output/plot/Rolling_EUE_by_Scenario.svg")
# plt.close()


# In[3]:


# Plot EUE_rolling_avg, faceted by Scenario_Name
g = sns.FacetGrid(CVaR_sorted_df_EUE, col="Scenario_Name", col_wrap=no_columns, sharex=False, sharey=True, height=5)
g.map_dataframe(sns.lineplot, x="Simulation_Number", y="EUE_rolling_avg")

# Add label to the last value in each facet
for ax, scenario in zip(g.axes.flat, CVaR_sorted_df_EUE['Scenario_Name'].unique()):
    df_scenario = CVaR_sorted_df_EUE[CVaR_sorted_df_EUE['Scenario_Name'] == scenario]
    last_row = df_scenario[df_scenario['Simulation_Number'] == df_scenario['Simulation_Number'].max()]
    if not last_row.empty:
        x = last_row['Simulation_Number'].values[0]
        y = last_row['EUE_rolling_avg'].values[0]
        ax.text(x, y, f"{y:.2f}", color="red", fontsize=12, ha="right", va="bottom")

g.set_axis_labels("Simulation", "EUE (MWh/year)")
g.set_titles(col_template="{col_name}")
# g.fig.suptitle("LOLH Rolling Average by Scenario", y=1.05)
g.fig.suptitle(f"EUE Rolling Average by Scenario | {CVaR_alpha} CVaR")
g.fig.tight_layout(rect=[0, 0, 0.85, 0.95])
g.fig.set_size_inches(figure_width, figure_height)
plt.savefig(f"{directory}/output/plot/Rolling_EUE_by_Scenario_CVaR.svg")
# plt.close()


# In[138]:


# Create an HTML file
with open(directory + '/output/combined_plots.html', 'w') as html_file:
    # Write the basic structure of the HTML file with a title
    html_file.write('<!DOCTYPE html>\n<html>\n<head>\n<title>RA Metric Results</title>\n')
    html_file.write('<style>\n')
    html_file.write('img { max-width: 100%; height: auto; display: block; margin: 10px auto; border: 2px solid black; }\n')
    html_file.write('div { text-align: center; margin: 20px; }\n')
    html_file.write('body { margin: 20px; }\n')
    html_file.write('</style>\n</head>\n<body>\n')
    html_file.write('<h1 style="text-align:center;">RA Metric Results</h1>\n')

    # Add datetime
    now = datetime.datetime.now().strftime('%Y-%m-%d %H:%M:%S')
    html_file.write(f'<p style="text-align:center;"><em>Generated on: {now}</em></p>\n')

    # Iterate over all files in the plot directory
    for filename in os.listdir(os.path.join(directory, 'output', 'plot')):
        if filename.endswith('.svg'):
            caption = os.path.splitext(filename)[0].replace('_', ' ')
            html_file.write(f'<div>\n')
            html_file.write(f'<img src="plot/{filename}" alt="{filename}">\n')
            # html_file.write(f'<h3 style="text-align:center;">{caption}</h3>\n')
            html_file.write(f'</div>\n')

    html_file.write('</body>\n</html>')