Addresses overly optimistic wind capacity factor issue

Snakefile

@@ -21,9 +21,18 @@ rule retrieve_supply_regions:
 
 rule retrieve_costs:
     output: 
-        costs = "data/technology_costs.csv"
+        costs = "data/technology_costs.csv",
+        heatrates = "data/heatrates.csv"
     script: "scripts/retrieve_costs.py"
 
+rule retrieve_fuel_costs:
+    input:
+        heatrates = "data/heatrates.csv"
+    output: 
+        fuel_costs = "data/thermal_fuel_costs.csv",
+        fuel_cost_timeseries = "data/thermal_fuel_cost_ts.csv"
+    script: "scripts/retrieve_fuel_prices.py"
+
 rule retrieve_load:
     output:
         load = "data/time_series/load.csv"
@@ -72,6 +81,7 @@ rule add_electricity:
         load = "data/time_series/load.csv",
         generators = "data/aggregated_generators.csv",
         costs = "data/technology_costs.csv",
+        fuel_cost_timeseries = "data/thermal_fuel_cost_ts.csv",
         wind_profile = "data/time_series/wind.csv",
         solar_profile = "data/time_series/solar.csv",
         base_network = "data/networks/base_network.nc",
@@ -96,5 +106,6 @@ rule plot_results:
         dispatch_figure = f"results/{results_folder}/figures/illinois_dispatch.png",
         capacity_figure = f"results/{results_folder}/figures/illinois_capacity.png",
         emissions_figure = f"results/{results_folder}/figures/illinois_emissions.png",
-        active_units_figure = f"results/{results_folder}/figures/illinois_active_units.png"
+        active_units_figure = f"results/{results_folder}/figures/illinois_active_units.png",
+        monthly_generation_figure = f"results/{results_folder}/figures/illinois_monthly_generation.png"
     script: "scripts/plot_results.py"

config.yml

@@ -1,13 +1,13 @@
 state_abbr: 'IL'
-version: "5.0-test"
-scenario: "no_growth_no_export"
+version: "1.0-test"
+scenario: "wind_cf"
 
 solver: 'cplex'  # 'cplex','highs','gurobi'
 geo_res: 'rto'  # accepts: 'rto' or 'county'
-time_res: 4  # hours 
+time_res: 1  # hours 
 load_filter: 60e3  # MW, load above this level will be removed as outliers.
-# total_demand: 185e6  # Annual MWh in the first year
-total_demand: 136e6  # Annual MWh in the first year
+total_demand: 185e6  # Annual MWh in the first year
+# total_demand: 136e6  # Annual MWh in the first year
 load_growth: 0.00  # % annual growth
 load_share:
   MISO-Z4: 0.33
@@ -30,6 +30,7 @@ co2_limits:
 solar_years: [2016]  # NSRDB only accepts years 1998-2020
 wind_years: [2009]  # WTK only goes from 2009-2013
 load_year: 2019
+fuel_cost_year: 2019
 
 rto_subba: ['0004','CE']  # MISO-Z4, ComEd
 region_names: ['MISO-Z4','ComEd']
@@ -110,6 +111,7 @@ turbine_params:
   diameter: 103 # m
   rated_power: 2.75  # MW
   air_density: 1.225  # kg/m^3
+  capacity_factor: 0.35  # average annual capacity factor
 
 nuclear_params:  # for existing reactors, 2023 NEI costs in context
   capacity_factor: 0.972

functions/data_functions.py

@@ -0,0 +1,9 @@
+import pandas as pd
+
+
+def load_heatrates():
+
+    heatrate_url = "https://www.eia.gov/electricity/annual/html/epa_08_01.html"
+    heatrates = pd.read_html(heatrate_url)[1].set_index("Year")
+
+    return heatrates

scripts/add_electricity.py

@@ -13,6 +13,8 @@
 idx_opts = {"rto": "balancing_authority_code",
             "county": "county"}
 growth_rates = snakemake.config['growth_rates']
+pudl_year = int(snakemake.config['fuel_cost_year'])
+wind_cf = float(snakemake.config['turbine_params']['capacity_factor'])
 
 BUILD_YEAR = 2025  # a universal build year place holder
 
@@ -68,6 +70,15 @@ def load_costs():
     return costs
 
 
+def load_costs_ts():
+
+    fuel_costs = pd.read_csv(snakemake.input.fuel_cost_timeseries,
+                             parse_dates=True,
+                             index_col=['report_date'])
+
+    return fuel_costs
+
+
 def load_existing_generators():
     generators = pd.read_csv(snakemake.input.generators,
                              index_col=idx_opts[scale])
@@ -235,7 +246,8 @@ def attach_generators(
         costs,
         generators=None,
         build_years=None,
-        model_year=None):
+        model_year=None,
+        costs_ts=None):
     carriers = ['Natural Gas', 'Biomass', 'Coal', 'Nuclear', 'Petroleum']
     for bus in n.buses.index:
         for item in costs.itertuples():
@@ -292,6 +304,20 @@ def attach_generators(
                     p_max_pu = 1
                     p_min_pu = 0
 
+                # time series marginal costs
+                if ((tech in ['CTAvgCF', 'CCAvgCF', 'IGCCAvgCF'])
+                        and isinstance(costs_ts, pd.DataFrame)):
+
+                    # select year to replicate
+                    cost_data = costs_ts.loc[str(pudl_year), carrier].values
+
+                    # get the VOM cost
+                    cost_data = cost_data + item.VOM
+
+                    marginal_cost = np.tile(cost_data, len(model_years))
+                else:
+                    marginal_cost = item.marginal_cost
+
                 n.add(class_name="Generator",
                       name=name,
                       bus=bus,
@@ -302,7 +328,7 @@ def attach_generators(
                       p_max_pu=p_max_pu,
                       carrier=carrier,
                       capital_cost=item.capital_cost,
-                      marginal_cost=item.marginal_cost,
+                      marginal_cost=marginal_cost,
                       lifetime=item.lifetime,
                       ramp_limit_down=ramp_limit_down,
                       ramp_limit_up=ramp_limit_up,
@@ -364,7 +390,6 @@ def attach_storage(
 
     n = pypsa.Network(snakemake.input.base_network)
     costs = load_costs()
-
     costs.to_csv("data/final_costs.csv")
     current_costs = costs.xs((slice(None), slice(None), 2020))
     current_costs = current_costs.reset_index()
@@ -374,6 +399,7 @@ def attach_storage(
     generators = load_existing_generators()
     build_years = load_build_years()
     emissions = load_emissions()
+    costs_ts = load_costs_ts()
 
     attach_load(n)
     add_carriers(n,
@@ -389,7 +415,8 @@ def attach_storage(
     attach_generators(n,
                       costs=current_costs,
                       generators=generators,
-                      build_years=build_years
+                      build_years=build_years,
+                      costs_ts=costs_ts
                       )
     attach_storage(n,
                    costs=current_costs,
@@ -408,7 +435,8 @@ def attach_storage(
                           )
         attach_generators(n,
                           costs=current_costs,
-                          model_year=year
+                          model_year=year,
+                          costs_ts=costs_ts
                           )
         attach_storage(n,
                        costs=current_costs,
@@ -430,4 +458,9 @@ def attach_storage(
         except (AttributeError, KeyError, TypeError):
             pass
 
+    # modify wind capacity factor
+    wind_gen = n.generators[n.generators.carrier == 'Wind'].index
+    n.generators_t.p_max_pu.loc[:, wind_gen] = ((n.generators_t.p_max_pu[wind_gen] / (
+        n.generators_t.p_max_pu[wind_gen].sum() / len(n.snapshots)) * wind_cf))
+
     n.export_to_netcdf(snakemake.output.elec_network)

scripts/plot_results.py

@@ -3,26 +3,41 @@
 import matplotlib.pyplot as plt
 import numpy as np
 
+TECH_ORDER = ['Nuclear',
+              'Coal',
+              'Natural Gas',
+              'Biomass',
+              'Petroleum',
+              'Solar',
+              'Wind']
+
+
+def power_by_carrier(n):
+    p_by_carrier = n.generators_t.p.T.groupby(
+        n.generators.carrier).sum().T
+
+    return p_by_carrier
+
 
 def plot_dispatch(n, year=2025, month=7):
 
     time = (year, f'{year}-0{month}')
-    p_by_carrier = n.generators_t.p.groupby(
-        n.generators.carrier, axis=1).sum().div(1e3)
+    p_by_carrier = power_by_carrier(n).div(1e3)
 
-    p_by_carrier = p_by_carrier[['Nuclear',
-                                 'Coal',
-                                 'Natural Gas',
-                                 'Biomass',
-                                 'Petroleum',
-                                 'Solar',
-                                 'Wind']]
+    p_by_carrier = p_by_carrier[TECH_ORDER]
 
     if not n.storage_units.empty:
         sto = n.storage_units_t.p.T.groupby(
             n.storage_units.carrier).sum().T.div(1e3)
         p_by_carrier = pd.concat([p_by_carrier, sto], axis=1)
 
+    # y-limits
+    y_min = -n.storage_units_t.p_store.max().max() / 1e3
+    y_max = n.loads_t.p_set.sum(axis=1).max() / 1e3
+    margin = 0.1
+    y_low = (1 + margin) * y_min
+    y_high = (1 + margin) * y_max
+
     fig, ax = plt.subplots(figsize=(12, 6))
 
     color = p_by_carrier.columns.map(n.carriers.color)
@@ -31,6 +46,7 @@ def plot_dispatch(n, year=2025, month=7):
         ax=ax,
         linewidth=0,
         color=color,
+        ylim=(y_low - margin, y_high + margin)
     )
 
     charge = p_by_carrier.where(
@@ -43,14 +59,13 @@ def plot_dispatch(n, year=2025, month=7):
             ax=ax,
             linewidth=0,
             color=charge.columns.map(n.carriers.color),
+            ylim=(y_low - margin, y_high + margin)
         )
 
     n.loads_t.p_set.sum(axis=1).loc[time].div(1e3).plot(ax=ax, c="k")
 
     ax.legend(loc=(1.05, 0))
-    ax.set_ylabel("GW")
-    ax.set_ylim(-n.storage_units_t.p_store.max().max() / 1e3 -
-                2.5, n.loads_t.p_set.sum(axis=1).max() / 1e3 + 2.5)
+    ax.set_ylabel("GW", fontsize=16)
     plt.tight_layout()
 
     return fig, ax
@@ -84,7 +99,7 @@ def plot_emissions(n, time_res):
     total_emissions = n.snapshot_weightings.generators @ emissions.sum(
         axis=1).div(1e6)
 
-    p_by_carrier = n.generators_t.p.groupby(n.generators.carrier, axis=1).sum()
+    p_by_carrier = power_by_carrier(n)
 
     p_by_carrier_year = (
         p_by_carrier.groupby(
@@ -149,6 +164,25 @@ def plot_active_units(n):
     return fig, ax
 
 
+def plot_monthly_generation(n, time_res):
+    p_by_carrier = power_by_carrier(n) * time_res
+
+    p_by_carrier = p_by_carrier.resample('ME', level='timestep').sum()
+
+    p_by_carrier = p_by_carrier[TECH_ORDER]
+
+    color = p_by_carrier.columns.map(n.carriers.color)
+
+    fig, ax = plt.subplots(figsize=(12, 8))
+    p_by_carrier.plot.area(ax=ax,
+                           color=color,
+                           fontsize=16)
+
+    ax.set_xlabel('')
+    ax.set_ylabel('Generation [MWh]', fontsize=18)
+    return fig, ax
+
+
 if __name__ == "__main__":
 
     n = pypsa.Network(snakemake.input.solved_network)
@@ -160,8 +194,12 @@ def plot_active_units(n):
     fig, ax = plot_capacity(n)
     plt.savefig(snakemake.output.capacity_figure)
 
-    fig, ax = plot_emissions(n, float(snakemake.config['time_res']))
+    time_res = float(snakemake.config['time_res'])
+    fig, ax = plot_emissions(n, time_res)
     plt.savefig(snakemake.output.emissions_figure)
 
     fig, ax = plot_active_units(n)
     plt.savefig(snakemake.output.active_units_figure)
+
+    fig, ax = plot_monthly_generation(n, time_res)
+    plt.savefig(snakemake.output.monthly_generation_figure)

scripts/retrieve_costs.py

@@ -1,5 +1,11 @@
+import numpy as np
 import pandas as pd
 from nrelpy.atb import ATBe
+import sys
+
+sys.path.append("functions")
+from data_functions import load_heatrates
+
 
 n_illinois_reactors = 11
 total_lwr_capacity = 12415.1
@@ -67,8 +73,7 @@
     fuel_prices = pd.read_html(prices_url)[1].set_index(
         pd.Series(range(2012, 2023))).T
 
-    heatrate_url = "https://www.eia.gov/electricity/annual/html/epa_08_01.html"
-    heatrates = pd.read_html(heatrate_url)[1].set_index("Year")
+    heatrates = load_heatrates()
 
     price_col = 'Average Cost (Dollars per MMBtu)'
 
@@ -92,7 +97,6 @@
 
     # converts btu/kwh to mmbtu/mwh
     ng_heatrate = heatrates.at[atb_params['atb_year'], 'Natural Gas'] / 1e3
-    # ng_heatrate =1  # converts btu/kwh to mmbtu/mwh
     coal_heatrate = heatrates.at[atb_params['atb_year'], 'Coal'] / 1e3
     petroleum_heatrate = heatrates.at[atb_params['atb_year'],
                                       'Petroleum'] / 1e3
@@ -126,3 +130,4 @@
     df_pivot.fillna(0., inplace=True)
 
     df_pivot.to_csv(snakemake.output.costs)
+    heatrates.to_csv(snakemake.output.heatrates)