Add Wind Battery Paper Results

paper/Wind_battery_paper_2024/pt_annual_elec_rev.json

@@ -0,0 +1 @@
+{"2": [24198264.200666964, 28863440.94444351, 33270001.148331646, 37553948.20013763, 41757285.92109379, 45790390.649100706, 49742041.05346024, 53616531.2151789, 57441243.264177024, 61146055.04620315], "4": [26484807.39413671, 33330521.03068218, 39886499.960852034, 46259267.61980139, 52469690.18084988, 58458257.11011698, 64057677.62524079, 70061680.93902785, 75315434.46978994, 79977239.79768948], "6": [27736015.65145083, 35871389.24923533, 43775243.70044307, 51444679.76847142, 58846054.069521256, 66119202.922327384, 72330428.278673, 78710547.73478292, 84648645.79854389, 90705881.93816274], "8": [28447709.64532416, 37280504.29448119, 45827714.21682936, 54142071.67840834, 62288769.2426817, 69835501.0690678, 76488082.21983999, 83465276.09713222, 89747857.13171645, 95542419.0763853], "10": [28966069.086988572, 38311889.645500205, 47368466.726949155, 56190622.607311904, 64383686.673088975, 71788787.76691541, 79741318.3425242, 86209873.1380863, 93760443.52576275, 100385375.49069096]}

paper/Wind_battery_paper_2024/utils.py

@@ -0,0 +1,309 @@
+import copy
+import os
+import pathlib
+import pandas as pd
+import numpy as np
+import re
+import json
+
+'''
+prescient_output_to_df, get_gdf and summarize_results are modified from Ben's fork. 
+Run summarize_results will give us the csv file of the simulation including DA dispatch, DA LMP, RT dispatch and RT LMP (like the sweep results)
+
+summarize_revenue is modified from Xian's fork. This function return a dictionary of DA dispatch, DA revenue, RT dispatch and RT revenue and total revenue.
+
+summarize_H2_revenue is to calculate the total hydrogen revenue.
+'''
+
+
+def prescient_output_to_df(file_name):
+    '''Helper for loading data from Prescient output csv.
+        Combines Datetimes into single column.
+    '''
+    df = pd.read_csv(file_name)
+    df['Datetime'] = \
+        pd.to_datetime(df['Date']) + \
+        pd.to_timedelta(df['Hour'], 'hour') + \
+        pd.to_timedelta(df['Minute'], 'minute')
+    df.drop(columns=['Date','Hour','Minute'], inplace=True)
+    # put 'Datetime' in front
+    cols = df.columns.tolist()
+    cols = cols[-1:]+cols[:-1]
+
+    return df[cols]
+
+
+def get_gdf(directory, generator_file_name, generator_name, dispatch_name):
+    gdf = prescient_output_to_df(os.path.join(directory, generator_file_name))
+    gdf = gdf[gdf["Generator"] == generator_name][["Datetime", dispatch_name, dispatch_name + " DA"]]
+    gdf.set_index("Datetime", inplace=True)
+    gdf.rename(columns={ dispatch_name : generator_name + " Dispatch", dispatch_name + " DA" : generator_name + " Dispatch DA"}, inplace=True)
+
+    return gdf
+
+
+def summarize_results(base_directory, generator_name, bus_name, output_directory):
+    """
+    Summarize Prescient runs for a single generator
+    Args:
+        base_directory (str) : the base directory name (without index)
+        generator_name (str) : The generator name to get the dispatch for. Looks in thermal_gens.csv and then renewable_gens.csv.
+        bus_name (str) : The bus to get the LMPs for.
+        output_directory (str) : The location to write the summary files to.
+    Returns:
+        None
+    """
+
+    # param_file = os.path.join(output_directory, "sweep_parameters.csv")
+
+    # figure out if renewable or thermal or virtual
+    generator_file_names = ("thermal_detail.csv", "renewables_detail.csv", "virtual_detail.csv")
+    dispatch_name_map = { "thermal_detail.csv" : "Dispatch",
+                          "renewables_detail.csv" : "Output",
+                          "virtual_detail.csv" : "Output",
+                        }
+
+    def _get_gen_df(generator_name):
+        for generator_file_name in generator_file_names:
+            gdf = pd.read_csv(os.path.join(base_directory, generator_file_name))["Generator"]
+            if generator_name in gdf.unique():
+                return generator_file_name
+        else: # no break
+            raise RuntimeError("Could not find output for generator "+generator_name)
+
+    generator_file_name = _get_gen_df(generator_name)
+    dispatch_name = dispatch_name_map[generator_file_name]
+
+
+    # if not os.path.isfile(os.path.join(directory, "overall_simulation_output.csv")):
+    #     raise Exception(f"For index {idx}, the simulation did not complete!")
+
+    gdf = get_gdf(base_directory, generator_file_name, generator_name, dispatch_name)
+    df_list = [gdf]
+    RT_names = [gdf.columns[0]]
+    DA_names = [gdf.columns[1]]
+
+
+    bdf = prescient_output_to_df(os.path.join(base_directory, "bus_detail.csv"))
+    bdf = bdf[bdf["Bus"] == bus_name][["Datetime","LMP","LMP DA"]]
+    bdf.set_index("Datetime", inplace=True)
+    df_list.append(bdf)
+    RT_names.append(bdf.columns[0])
+    DA_names.append(bdf.columns[1])
+
+    odf = pd.concat(df_list, axis=1)[[*RT_names,*DA_names]]
+    lmp_array = odf["LMP"].to_numpy()
+    odf["LMP"] = np.clip(lmp_array, 0, 500)
+    odf.to_csv(os.path.join(output_directory, f"{base_directory + '_' + generator_name}.csv"))
+
+
+def summarize_revenue(sim_id, result_dir, gen_detail, bus_name, gen_name, cap_rt_lmp = False):
+
+    '''
+    Summary the total DA and RT dispatch and revenue.
+    Args:
+        sim_id(int): simulation id
+        result_dir(str): the result directory name
+        gen_detail(str): generator results file
+        bus_name(str): the bus to get the LMPs for.
+        gen_name(str): name of the generator
+        cap_rt_lmp: if we are going to cap the rt_lmp.
+    Returns:
+        summary(dict): the total DA and RT dispatch and revenue information
+    '''
+
+    df = pd.read_csv(os.path.join(result_dir, gen_detail))
+    df = df.loc[df["Generator"] == gen_name]
+    df["Time Index"] = range(len(df))
+    df.rename(columns={"Output": "Dispatch", "Output DA": "Dispatch DA"}, inplace=True)
+
+    bus_df = pd.read_csv(os.path.join(result_dir, "bus_detail.csv"))
+    bus_df = bus_df.loc[bus_df["Bus"] == bus_name]
+    bus_df["Time Index"] = range(len(bus_df))
+
+    df = df.merge(bus_df, how="left", left_on="Time Index", right_on="Time Index")
+
+    if cap_rt_lmp == True:
+        lmp_array = df["LMP"].to_numpy()
+        df["LMP"] = pd.DataFrame(np.clip(lmp_array, 0, 500))
+
+    df["Revenue DA"] = df["Dispatch DA"] * df["LMP DA"]
+    df["Revenue RT"] = (df["Dispatch"] - df["Dispatch DA"]) * df["LMP"]
+    df["Total Revenue"] = df["Revenue DA"] + df["Revenue RT"]
+
+    df = df[["Dispatch", "Dispatch DA", "Revenue DA", "Revenue RT", "Total Revenue"]]
+
+    summary = df.sum().to_dict()
+    summary["sim_id"] = sim_id
+
+    return summary
+
+def summarize_rt_revenue(sim_id, result_dir, gen_detail, bus_name, gen_name, cap_rt_lmp = False):
+
+    '''
+    Summary the total DA and RT dispatch and revenue.
+    Args:
+        sim_id(int): simulation id
+        result_dir(str): the result directory name
+        gen_detail(str): generator results file
+        bus_name(str): the bus to get the LMPs for.
+        gen_name(str): name of the generator
+        cap_rt_lmp: if we are going to cap the rt_lmp.
+    Returns:
+        summary(dict): the total DA and RT dispatch and revenue information
+    '''
+
+    df = pd.read_csv(os.path.join(result_dir, gen_detail))
+    df = df.loc[df["Generator"] == gen_name]
+    df["Time Index"] = range(len(df))
+    df.rename(columns={"Output": "Dispatch", "Output DA": "Dispatch DA"}, inplace=True)
+
+    bus_df = pd.read_csv(os.path.join(result_dir, "bus_detail.csv"))
+    bus_df = bus_df.loc[bus_df["Bus"] == bus_name]
+    bus_df["Time Index"] = range(len(bus_df))
+
+    df = df.merge(bus_df, how="left", left_on="Time Index", right_on="Time Index")
+
+    if cap_rt_lmp == True:
+        lmp_array = df["LMP"].to_numpy()
+        df["LMP"] = pd.DataFrame(np.clip(lmp_array, 0, 500))
+
+    df["Total Revenue"] = df["Dispatch"] * df["LMP"]
+
+    df = df[["Dispatch", "Dispatch DA", "Total Revenue"]]
+
+    summary = df.sum().to_dict()
+    summary["sim_id"] = sim_id
+
+    return summary
+
+
+def summarize_in_df(result_dir, gen_detail, bus_name, gen_name):
+
+    '''
+    Summary the hourly DA and RT dispatch and LMP information.
+    Args:
+        result_dir(str): the result directory name
+        gen_detail(str): generator results file
+        bus_name(str) the bus to get the LMPs for.
+        gen_name(str): name of the generator
+    Returns:
+        df(dataframe): the hourly DA and RT dispatch and revenue information
+    '''
+
+    df = pd.read_csv(os.path.join(result_dir, gen_detail))
+    df = df.loc[df["Generator"] == gen_name]
+    df["Time Index"] = range(len(df))
+    df.rename(columns={"Output": "Dispatch", "Output DA": "Dispatch DA"}, inplace=True)
+
+    bus_df = pd.read_csv(os.path.join(result_dir, "bus_detail.csv"))
+    bus_df = bus_df.loc[bus_df["Bus"] == bus_name]
+    bus_df["Time Index"] = range(len(bus_df))
+
+    df = df.merge(bus_df, how="left", left_on="Time Index", right_on="Time Index")
+
+    df["Revenue DA"] = df["Dispatch DA"] * df["LMP DA"]
+    df["Revenue RT"] = (df["Dispatch"] - df["Dispatch DA"]) * df["LMP"]
+    df["Total Revenue"] = df["Revenue DA"] + df["Revenue RT"]
+
+    df = df[["Dispatch", "Dispatch DA", "Revenue DA", "Revenue RT", "Total Revenue"]]
+
+    return df
+
+
+def summarize_H2_revenue(df, PEM_size, H2_price, gen_name):
+
+    '''
+    Summary the H2 revenue. 
+    Args:
+        df(dataframe): dataframe for the hourly dispatch information (from summarize_in_df)
+        PEM_size(float): the size of the PEM, MW.
+        H2_price(float): Hydrogen price, $/kg
+        gen_name(str): name of the generator
+    Returns:
+        (dict): the DA and RT dispatch and revenue information
+    '''
+    # read the wind data
+    df_wind = pd.read_csv("Real_Time_wind_hourly.csv")[gen_name]
+
+    # calculate excess elec (can be used for PEM or curtailed)
+    rt_dispatch = df["Dispatch"]
+    df_ph = df_wind - rt_dispatch
+    excess_elec = df_ph.to_numpy()
+
+    # Find and remove Nan indexes.
+    nan_idx = np.isnan(excess_elec)
+    excess_elec_without_nan = excess_elec[~nan_idx]
+
+    # calculate the electricity that PEM used to produce H2
+    pem_elec = np.clip(excess_elec_without_nan, 0, PEM_size)
+    # calculate H2 revenue
+    eh_rate = 54.953 # kWh/kg
+    h2_revenue = sum(pem_elec)/eh_rate*H2_price*1000
+
+    df_H2 = {}
+    df_H2["Total PEM electricity"] = np.sum(pem_elec)
+    df_H2["Total H2 revenue"] = h2_revenue
+    return df_H2
+
+
+def calculate_NPV(annual_revenue, wind_size, battery_size, duration = 4, year= 2023, scenario="moderate", discount_rate = 0.05, OM_cost = False, extant_wind = True):
+    discount_rate = discount_rate               # discount rate
+    N = 30                                      # years
+    PA = ((1+discount_rate)**N - 1)/(discount_rate*(1+discount_rate)**N) 
+
+    with open("wind_battery_cost_parameter.json", "rb") as f:
+        price_dict = json.load(f)
+
+    # the default argment in the list is [2hr, 4hr, 6hr, 8hr, 10hr]
+    arg_duration = int(duration/2 - 1)
+
+    batt_op_cost = price_dict["battery"]["fixed_om"][scenario][str(year)][arg_duration]
+    batt_cap_cost_kw = price_dict["battery"]["batt_cap_cost_param"][scenario][str(year)][0]
+    batt_cap_cost_kwh = price_dict["battery"]["batt_cap_cost_param"][scenario][str(year)][1]
+
+    wind_cap_cost = price_dict["wind"]["capital"][scenario][str(year)][0]
+    wind_op_cost = price_dict["wind"]["fixed_om"][scenario][str(year)][0]
+
+    capital_cost_wind = wind_cap_cost*wind_size*1000
+    capital_cost_battery = batt_cap_cost_kw*battery_size*1000 + batt_cap_cost_kwh*battery_size*duration*1000
+    op_cost = wind_op_cost*wind_size*1000 + battery_size*batt_op_cost*1000
+
+    if OM_cost:
+        if extant_wind:
+            NPV = (annual_revenue - op_cost)*PA - capital_cost_battery
+        else:
+            NPV = (annual_revenue - op_cost)*PA - capital_cost_wind - capital_cost_battery
+    else:
+        if extant_wind:
+            NPV = annual_revenue * PA - capital_cost_battery
+        else:    
+            NPV = annual_revenue * PA - capital_cost_wind - capital_cost_battery
+
+    return NPV
+
+
+def main():
+    gen_detail = "thermal_detail.csv"
+    generator_name = "303_WIND_1"
+    bus_name = "Caesar"
+    folder_path = "new_wind_battery_sweep_sb_small"
+    # get the list of folder names that stores the simulation results
+    file_list = os.listdir(folder_path)
+    # a dictionary that is going to save the simulation revenue results
+    result_dict = {}
+    for name in file_list:
+        # split the name and take the last value, which is the sim_id.
+        sim_path = os.path.join(folder_path, name)
+        sim_id = int(re.split("_", name)[-1])
+        res = summarize_revenue(sim_id, sim_path, gen_detail, bus_name, generator_name, cap_rt_lmp = True)
+        result_dict[sim_id] = res
+
+    with open (folder_path + 'cap_rt_lmp_summary', "w") as f:
+        json.dump(result_dict, f)
+
+    return result_dict
+
+
+if __name__ == "__main__":
+    main()

paper/Wind_battery_paper_2024/wind_battery_cost_parameter.json

@@ -0,0 +1 @@
+{"battery": {"capital": {"conservative": {"2023": [809.46, 1353.058, 1978.254, 2562.651, 3147.048], "2030": [681.622, 1173.726, 1665.829, 2157.932, 2650.036], "2040": [681.622, 1173.726, 1665.829, 2157.932, 2650.036], "2050": [681.622, 1173.726, 1665.829, 2157.932, 2650.036]}, "moderate": {"2023": [746.037, 1255.708, 1765.38, 2275.051, 2784.723], "2030": [579.552, 894.692, 1209.833, 1524.973, 1840.113], "2040": [507.455, 783.3, 1059.145, 1334.991, 1610.836], "2050": [435.108, 671.464, 907.819, 1144.174, 1380.53]}, "advanced": {"2023": [720.301, 1240.329, 1760.357, 2280.384, 2800.412], "2030": [394.893, 679.989, 965.085, 1250.181, 1535.277], "2040": [345.531, 594.99, 844.449, 1093.908, 1343.367], "2050": [296.169, 509.991, 723.814, 937.636, 1151.458]}}, "fixed_om": {"conservative": {"2023": [20.237, 34.846, 49.456, 64.066, 78.676], "2030": [17.041, 29.343, 41.646, 53.948, 66.251], "2040": [17.041, 29.343, 41.646, 53.948, 66.251], "2050": [17.041, 29.343, 41.646, 53.948, 66.251]}, "moderate": {"2023": [18.651, 31.39, 44.009, 57.01, 69.618], "2030": [14.489, 22.367, 30.246, 38.124, 46.003], "2040": [12.686, 19.583, 26.479, 33.375, 40.271], "2050": [10.878, 16.787, 22.695, 28.604, 34.513]}, "advanced": {"2023": [18.008, 31.008, 44.009, 57.01, 70.01], "2030": [9.872, 17, 24.127, 31.255, 38.382], "2040": [8.638, 14.875, 21.111, 27.348, 33.584], "2050": [7.404, 12.75, 18.095, 23.441, 28.786]}}, "batt_cap_cost_param": {"conservative": {"2023": [204.66349791171595, 294.2384502770722], "2030": [189.51879999994418, 246.05170000000768], "2040": [189.51879999994418, 246.05170000000768], "2050": [189.51879999994418, 246.05170000000768]}, "moderate": {"2023": [236.365, 254.835], "2030": [264.41169999997675, 157.5701500000031], "2040": [231.60950000000406, 137.92264999999932], "2050": [198.75279999998963, 118.17770000000148]}, "advanced": {"2023": [200.27350000000357, 260.0138499999994], "2030": [109.79700000000014, 142.54799999999997], "2040": [96.07200000000005, 124.72949999999999], "2050": [82.34670000003243, 106.91114999999567]}}}, "wind": {"capital": {"conservative": {"2023": [1323.4], "2030": [1000], "2040": [950], "2050": [900]}, "moderate": {"2023": [1308], "2030": [950], "2040": [855], "2050": [760]}, "advanced": {"2023": [1233.4], "2030": [700], "2040": [612.5], "2050": [525]}}, "fixed_om": {"conservative": {"2023": [43], "2030": [43], "2040": [42.025], "2050": [41.05]}, "moderate": {"2023": [41.78], "2030": [38.95], "2040": [36.029], "2050": [33.108]}, "advanced": {"2023": [40.414], "2030": [34.38], "2040": [29.223], "2050": [24.066]}}}}