Merge pull request #45 from ucsusa/adv-nuclear-test

Final model update

.gitignore

@@ -169,6 +169,6 @@ data/*.csv
 time_series
 spatial_data
 networks
-notebooks
+experimental
 
 results

config.yml

@@ -1,18 +1,21 @@
 state_abbr: 'IL'
-version: "11.0"
-scenario: "multi-invest-optimization-test"
+version: "3.0"
+scenario: "illinois"
 
-solver: 'highs'  # 'cplex','highs','gurobi'
+solver: 'cplex'  # 'cplex','highs','gurobi'
 geo_res: 'rto'  # accepts: 'rto' or 'county'
 myopic: False
 multi_investment_periods: True
-time_res: 8  # 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
-load_growth: 0.00  # % annual growth
+total_demand: 140e6  # Annual MWh in the first year
+load_growth: 0.02  # % annual growth
 random_seed: 1234
 
+ptc_value: 30  # $/MWh, solar, pv, advanced nuclear
+itc_value: 0.4  # % tax credit for renewables and storage
+
 # model_years: [2025, 2030, 2035, 2040, 2045]
 model_years: [2030, 2035, 2040, 2045, 2050]
 # model_years: [2025]
@@ -64,6 +67,10 @@ rto_subba: ['0004','CE']  # MISO-Z4, ComEd
 region_names: ['MISO-Z4','ComEd']
 
 discount_rate: 0.07
+battery_wacc_nom: 0.055
+battery_wacc_real: 0.038  
+
+rate: 'WACC Real'  # accepts: 'WACC Nominal', 'discount_rate'
 
 growth_rates:  # MW/year
   Nuclear: "inf"
@@ -139,7 +146,13 @@ 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
+  current_capacity_factor: 0.382  # average annual capacity factor
+  capacity_factor: 
+    2030: 0.39
+    2035: 0.3975
+    2040: 0.405
+    2045: 0.4125
+    2050: 0.42
 
 nuclear_params:  # for existing reactors, 2023 NEI costs in context
   capacity_factor: 0.972

scripts/add_electricity.py

@@ -13,7 +13,8 @@
             "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'])
+wind_cf = float(snakemake.config['turbine_params']['current_capacity_factor'])
+wind_cf_list = snakemake.config['turbine_params']['capacity_factor']
 try: 
     retirements_df = pd.DataFrame(snakemake.config['retirements']).fillna(0)
 except:
@@ -23,6 +24,16 @@
     capacity_limits_df = pd.DataFrame(snakemake.config['capacity_max']).fillna(0)
 except:
     capacity_limits_df = None
+
+try:
+    itc_credit = float(snakemake.config['itc_value'])
+except:
+    itc_credit = 0
+
+try:
+    ptc_credit = float(snakemake.config['ptc_value'])
+except:
+    ptc_credit = 0
 
 BUILD_YEAR = 2025  # a universal build year place holder
 
@@ -56,15 +67,17 @@ def load_costs():
     r = float(snakemake.config['discount_rate'])
     lifetimes = snakemake.config['lifetime']
 
-    costs = costs.assign(rate=r, lifetime=20)
+    costs = costs.assign(discount_rate=r, lifetime=20)
 
     # assign lifetimes to technologies
     carriers = snakemake.config['atb_params']['carrier']
     for carrier in carriers:
         costs.loc[(carrier, slice(None), slice(None)),
                   'lifetime'] = float(lifetimes[carrier])
 
-    annuity_col = annuity(costs["rate"], costs["lifetime"])
+
+    rate_col = snakemake.config['rate']
+    annuity_col = annuity(costs[rate_col], costs["lifetime"])
 
     costs = costs.assign(
         capital_cost=(
@@ -101,7 +114,7 @@ def load_build_years():
     return build_years
 
 
-def linear_growth(init_value, start_year, growth_rate, end_year=2050):
+def linear_growth(init_value, start_year, growth_rate, end_year=2051):
     def model(x, init_val, start, rate):
         return rate * init_val * (x - start) + init_val
     years = np.arange(start_year, end_year, 1).astype('int')
@@ -225,12 +238,14 @@ def attach_renewables(
                         continue
                     name = f"{bus} {tech} EXIST"
                     extendable = False
+                    marginal_cost = item.marginal_cost
                 elif model_year:
                     build_year = model_year
                     p_nom = 0.0
                     name = f"{bus} {tech} {model_year}"
                     extendable = tech in snakemake.config['extendable_techs']
                     build_year = model_year
+                    marginal_cost = item.marginal_cost - ptc_credit
 
                 p_max_pu = re_profile[bus]
 
@@ -243,7 +258,7 @@ def attach_renewables(
                       p_nom_extendable=extendable,
                       carrier=carrier,
                       capital_cost=item.capital_cost,
-                      marginal_cost=item.marginal_cost,
+                      marginal_cost=marginal_cost,
                       lifetime=item.lifetime,
                       build_year=build_year)
     return
@@ -311,13 +326,20 @@ def attach_generators(
                 else:
                     p_max_pu = 1
                     p_min_pu = 0
+
+
+                # add ptc
+                if (tech in ['Biopower', 'NuclearSMR']) & ('EXIST' not in name):
+                    capital_cost = item.capital_cost * (1-itc_credit)
+                else:
+                    capital_cost = item.capital_cost
 
                 # 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
+                    cost_data = costs_ts.loc[str(pudl_year), carrier].values[:8760]
 
                     # get the VOM cost
                     cost_data = cost_data + item.VOM
@@ -335,7 +357,7 @@ def attach_generators(
                       p_min_pu=p_min_pu,
                       p_max_pu=p_max_pu,
                       carrier=carrier,
-                      capital_cost=item.capital_cost,
+                      capital_cost=capital_cost,
                       marginal_cost=marginal_cost,
                       lifetime=item.lifetime,
                       ramp_limit_down=ramp_limit_down,
@@ -384,7 +406,7 @@ def attach_storage(
                       p_nom_min=p_nom,
                       p_nom_extendable=extendable,
                       carrier=carrier,
-                      capital_cost=item.capital_cost,
+                      capital_cost=item.capital_cost*(1-itc_credit),
                       marginal_cost=item.marginal_cost,
                       lifetime=item.lifetime,
                       max_hours=float(tech.split(' ')[0].strip('Hr')),
@@ -521,9 +543,9 @@ def add_energy_max(n):
 
     # add new technology
     for year in model_years:
-        # current_costs = costs.xs((slice(None), slice(None), year))
-        # current_costs = current_costs.reset_index()
-        # current_costs.loc[current_costs['technology_alias']=='Solar', 'techdetail'] = 'Utility PV'
+        current_costs = costs.xs((slice(None), slice(None), year))
+        current_costs = current_costs.reset_index()
+        current_costs.loc[current_costs['technology_alias']=='Solar', 'techdetail'] = 'Utility PV'
         attach_renewables(n,
                           costs=current_costs,
                           model_year=year
@@ -563,7 +585,17 @@ def add_energy_max(n):
 
     # 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))
+
+    wind_vintage = wind_gen[wind_gen.str.contains('EXIST')] 
+
+    n.generators_t.p_max_pu.loc[:, wind_vintage] = ((n.generators_t.p_max_pu[wind_vintage] / 
+                                                    (n.generators_t.p_max_pu[wind_vintage].sum() / 
+                                                    (len(n.snapshots))) * wind_cf))
+
+    for year in model_years:
+        wind_vintage = wind_gen[wind_gen.str.contains(str(year))] 
+        n.generators_t.p_max_pu.loc[:, wind_vintage] = ((n.generators_t.p_max_pu[wind_vintage] / 
+                                                     (n.generators_t.p_max_pu[wind_vintage].sum() / 
+                                                      (len(n.snapshots))) * float(wind_cf_list[year])))
 
     n.export_to_netcdf(snakemake.output.elec_network)

scripts/nuclear_sensitivity.py

@@ -0,0 +1,33 @@
+import pypsa
+import numpy as np
+import matplotlib.pyplot as plt
+import pandas as pd
+import scipy as sp
+from glob import glob
+from tqdm import tqdm
+import seaborn as sb
+from pathlib import Path
+
+if __name__ == "__main__":
+    results_path = Path("../pypsa-illinois/results/advanced_nuclear_sensitivity/")
+    results_path.mkdir(parents=True, exist_ok=True)
+
+    n = pypsa.Network("../pypsa-illinois/results/advanced_nuclear_v1.1/networks/illinois_solved.nc")
+
+    delta = 0.25
+    costs = np.arange(1.5, 3+delta, delta)
+
+    smr_cost_2030 = n.generators.loc[n.generators.index.str.contains('SMR'), 'capital_cost'].unique()[0]
+
+    for cost in tqdm(costs):
+        n_test = n.copy()
+        n_test.generators.loc[n_test.generators.index.str.contains('SMR'), 'capital_cost'] = smr_cost_2030*cost
+
+        n_test.optimize(solver_name='cplex', multi_investment_periods=True)
+
+        scenario = f"cost-2023_growth_0.01_demand-1.36E+08_atb-Moderate-X-{cost}_v1.1"
+
+        network_folder = results_path/scenario
+        network_folder.mkdir(exist_ok=True)
+
+        n.export_to_netcdf(str(network_folder/"illinois_solved.nc"))

scripts/retrieve_costs.py

@@ -48,25 +48,65 @@
                                      'core_metric_variable'],
                               columns='core_metric_parameter',
                               values='value')
-    df_pivot = df_pivot.xs((atb_params['scenario'],
-                            atb_params['case'],
-                            str(atb_params['crp']),
-                            slice(None),
-                            slice(None),
-                            slice(None)
-                            ))
-    df_pivot = df_pivot.loc[(atb_params['carrier'],
+    # breakpoint()
+    df_pivot = df_pivot.xs((atb_params['scenario'], 
+                            atb_params['case'], 
+                            slice(None), 
+                            slice(None), 
+                            slice(None), 
+                            slice(None)))
+
+    wacc_data = df_pivot.loc[('*', 
+                              atb_params['carrier'], 
+                              slice(None), 
+                              slice(None)), 
+                             ['WACC Nominal', 'WACC Real']]
+    wacc_data.to_csv("data/wacc.csv")
+
+    # df_pivot = df_pivot.xs((atb_params['scenario'],
+    #                         atb_params['case'],
+    #                         str(atb_params['crp']),
+    #                         slice(None),
+    #                         slice(None),
+    #                         slice(None)
+    #                         ))
+
+    df_pivot = df_pivot.loc[(str(atb_params['crp']),
+                             atb_params['carrier'],
                              atb_params['technology'],
                              slice(None)),
                             ['CAPEX', 'Fixed O&M',
-                             'Variable O&M', 'Fuel']]\
-        .dropna(axis=0, how='all')\
-        .drop_duplicates(keep='first')
+                             'Variable O&M', 'Fuel']]
 
     df_pivot.rename(columns={'Fixed O&M': 'FOM',
                              'Variable O&M': 'VOM',
                              },
                     inplace=True)
+
+    df_pivot = df_pivot.reset_index()
+    wacc_data = wacc_data.reset_index()
+
+    wacc_data['crpyears'] = str(atb_params['crp'])
+
+    merged_data = df_pivot.merge(wacc_data, 
+                                 on=['technology_alias',
+                                     'core_metric_variable',
+                                     'crpyears'],
+                                 how='left')
+
+    merged_data = merged_data.drop(columns=['techdetail_y','crpyears'])
+    merged_data.rename(columns={'techdetail_x':'techdetail'}, inplace=True)
+    merged_data = merged_data.loc[~((merged_data['technology_alias']=='Wind') 
+                                    & (merged_data['techdetail'] == 'Utility PV'))]
+    merged_data = merged_data.loc[~((merged_data['technology_alias']=='Solar') 
+                                    & (merged_data['techdetail'] == 'Land-Based Wind'))]
+
+    df_pivot = merged_data.set_index(['technology_alias','techdetail','core_metric_variable'])
+
+
+    wacc_data = wacc_data.drop(columns=['crpyears', 
+                                        'techdetail']).set_index(['technology_alias',
+                                                                  'core_metric_variable'])
     # add natural gas and coal costs
 
     prices_url = "https://www.eia.gov/electricity/annual/html/epa_07_04.html"
@@ -111,22 +151,38 @@
                  'Fuel'] = coal_price * coal_heatrate
     # add petroleum
     df_t = df_pivot.T
-
     # [capital cost, fixed om cost, variable om cost, fuel cost]
     for year in range(2020, 2051, 1):
         # update costs for existing nuclear, these costs are from 2022
         df_t['Nuclear', 'LWR', year] = [
-            nuclear_cap_cost, nuclear_fixed_om, 0, nuclear_fuel]
+            nuclear_cap_cost, 
+            nuclear_fixed_om, 
+            0, 
+            nuclear_fuel, 
+            wacc_data.loc[('Nuclear',year), 'WACC Nominal'],
+            wacc_data.loc[('Nuclear',year), 'WACC Real'],]
 
         # these costs are from 2021 and should be updated to reflect inflation.
         # from here https://www.eia.gov/electricity/generatorcosts/
         df_t['Petroleum', 'Petroleum', year] = [1158,
                                                 27.94,
                                                 1.78,
                                                 (petroleum_price *
-                                                 petroleum_heatrate)]
+                                                 petroleum_heatrate),
+                                                float(snakemake.config['discount_rate']),
+                                                float(snakemake.config['discount_rate'])]
     df_pivot = df_t.T
 
+    # set values for the battery WACC
+    df_pivot.loc[('Batteries', 
+                '4Hr Battery Storage', 
+                slice(None)), 
+                ['WACC Real']] = float(snakemake.config['battery_wacc_real'])
+    df_pivot.loc[('Batteries', 
+                '4Hr Battery Storage', 
+                slice(None)), 
+                ['WACC Nominal']] = float(snakemake.config['battery_wacc_nom'])
+
     df_pivot.fillna(0., inplace=True)
 
     df_pivot.to_csv(snakemake.output.costs)

sensitivity-script.py

@@ -12,9 +12,9 @@
 
 # Define the variables and their possible values
 variable_ranges = {
-    'fuel_cost_year': [2018, 2019, 2020, 2021, 2022, 2023],
-    'load_growth': [0, 1],
-    'total_demand': [136e6, 185e6],
+    'fuel_cost_year': [2023],
+    'load_growth': [0.0],  #[0.01, 0.02, 0.025],
+    'total_demand': [140e6, 185e6],
     'atb_scenario': ['Moderate']  # 'Conservative', 'Advanced'
 }