diff --git a/Snakefile b/Snakefile index c0cbd90..49954a2 100644 --- a/Snakefile +++ b/Snakefile @@ -28,7 +28,8 @@ rule retrieve_load: rule retrieve_existing_generators: output: generators = "data/existing_generators.csv", - aggregated = "data/aggregated_generators.csv" + aggregated = "data/aggregated_generators.csv", + build_year = "data/build_year_generators.csv" script: "scripts/retrieve_generators.py" rule retrieve_renewable_profiles: @@ -39,6 +40,11 @@ rule retrieve_renewable_profiles: solar = "data/time_series/solar.csv" script: "scripts/retrieve_renewables.py" +rule retrieve_co2_emissions: + output: + emissions = "data/technology_emissions.csv" + script: "scripts/retrieve_emissions.py" + rule build_topology: input: supply_regions = "data/spatial_data/supply_regions.shp" @@ -63,17 +69,27 @@ rule add_electricity: costs = "data/technology_costs.csv", wind_profile = "data/time_series/wind.csv", solar_profile = "data/time_series/solar.csv", - base_network = "data/networks/base_network.nc" + base_network = "data/networks/base_network.nc", + emissions = "data/technology_emissions.csv", + build_years = "data/build_year_generators.csv" output: elec_network = "data/networks/electricity_network.nc", final_costs = "data/final_costs.csv" script: "scripts/add_electricity.py" - rule solve_network: input: elec_network = "data/networks/electricity_network.nc" output: - solved_network = "results/networks/illinois_solved.nc", - dispatch_figure = "results/figures/illinois_dispatch.png" - script: "scripts/solve_network.py" \ No newline at end of file + solved_network = "results/networks/illinois_solved.nc" + script: "scripts/solve_network.py" + +rule plot_results: + input: + solved_network = "results/networks/illinois_solved.nc" + output: + dispatch_figure = "results/figures/illinois_dispatch.png", + capacity_figure = "results/figures/illinois_capacity.png", + emissions_figure = "results/figures/illinois_emissions.png", + active_units_figure = "results/figures/illinois_active_units.png" + script: "scripts/plot_results.py" \ No newline at end of file diff --git a/config.yml b/config.yml index 4c448bb..7b49e92 100644 --- a/config.yml +++ b/config.yml @@ -1,12 +1,27 @@ state_abbr: 'IL' 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: 136e6 # Annual MWh in the first year +load_growth: 0.00 # % annual growth +load_share: + MISO-Z4: 0.33 + ComEd: 0.67 -model_year: 2025 +model_years: [2025, 2030, 2035, 2040, 2045] +plot_year: 2045 -solar_year: 2019 # NSRDB only accepts years 1998-2020 -wind_year: 2011 # WTK only goes from 2009-2013 +co2_limits: + 2020: 31.25 + 2025: 25 # MT + 2030: 18.75 + 2035: 12.5 + 2040: 6.25 + 2045: 0.0 + +solar_years: [2016,2017,2018,2019,2020] # NSRDB only accepts years 1998-2020 +wind_years: [2009,2010,2011,2012,2013] # WTK only goes from 2009-2013 load_year: 2019 rto_subba: ['0004','CE'] # MISO-Z4, ComEd @@ -14,14 +29,25 @@ region_names: ['MISO-Z4','ComEd'] discount_rate: 0.07 +growth_rates: # MW/year + Nuclear: 1e3 + Biomass: 1e3 + Natural Gas: 10e3 + Coal: 10e3 + Petroleum: 10e3 + Batteries: 5e3 + Solar: 5e3 + Wind: 5e3 + lifetime: Nuclear: 80 Biomass: 40 Natural Gas: 40 - Coal: 40 + Coal: 50 Solar: 20 Wind: 20 Batteries: 15 + Petroleum: 40 atb_params: atb_year: 2022 # the ATB publication year @@ -56,11 +82,10 @@ atb_params: lines: # this is an assumption v_nom: 400. s_nom: 0. - x: 1 # reactance - r: 1 # resistance + x: 0.001 # reactance + r: 0.001 # resistance s_nom_extendable: True - turbine_params: cut_in: 3.0 # m/s cut_out: 25.0 # m/s @@ -69,35 +94,25 @@ turbine_params: rated_power: 2.75 # MW air_density: 1.225 # kg/m^3 -# extendable_techs: [ -# 'CCAvgCF', -# 'CTAvgCF', -# 'Biopower', -# 'Land-Based Wind', -# 'Utility PV', -# 'IGCCAvgCF', -# 'LWR', -# 'NuclearSMR', -# '10Hr Battery Storage', -# '8Hr Battery Storage', -# '6Hr Battery Storage', -# '4Hr Battery Storage', -# '2Hr Battery Storage' -# ] +nuclear_params: # for existing reactors, 2023 NEI costs in context + capacity_factor: 0.972 + capital_cost: 6.88 # $/MWh + fixed_om: 18.68 # $/MWh + fuel: 5.37 # $/MWh extendable_techs: [ # 'CCAvgCF', # 'CTAvgCF', - # 'Biopower', - # 'Land-Based Wind', - # 'Utility PV', # 'IGCCAvgCF', + # 'Biopower', + 'Land-Based Wind', + 'Utility PV', # 'LWR', # 'NuclearSMR', # '10Hr Battery Storage', # '8Hr Battery Storage', # '6Hr Battery Storage', - # '4Hr Battery Storage', + '4Hr Battery Storage', # '2Hr Battery Storage' ] diff --git a/dag.png b/dag.png index 15b7047..89892c2 100644 Binary files a/dag.png and b/dag.png differ diff --git a/functions/nrel_data_api.py b/functions/nrel_data_api.py index f1aab75..9f6467e 100644 --- a/functions/nrel_data_api.py +++ b/functions/nrel_data_api.py @@ -21,7 +21,7 @@ 'year':2019, 'leap_day':'false', 'selector':'POINT', - 'utc':'false', + 'utc':'true', 'interval':'60', 'attr_list':['ghi']} diff --git a/scripts/add_electricity.py b/scripts/add_electricity.py index a60ecf2..c4fb0a2 100644 --- a/scripts/add_electricity.py +++ b/scripts/add_electricity.py @@ -6,19 +6,30 @@ from tqdm import tqdm import pypsa -model_year = snakemake.config['model_year'] +version = 22 +model_years = np.array(snakemake.config['model_years']).astype('int') +resolution = int(snakemake.config['time_res']) +scale = snakemake.config['geo_res'] +idx_opts = {"rto":"balancing_authority_code", + "county":"county"} +growth_rates = snakemake.config['growth_rates'] + +BUILD_YEAR = 2025 # a universal build year place holder + def annuity(r, n): return r / (1-1/(1+r)**n) def create_snapshots(): - snapshots = pd.date_range(start=f"{model_year}-01-01", - end=f"{model_year+1}-01-01", - freq='1h', - inclusive='left') + timestamps = pd.DatetimeIndex([]) + for year in model_years: + period = pd.date_range(start=f"{year}-01-01", + freq=f"{resolution}h", + periods=8760/resolution) + timestamps = timestamps.append(period) - return snapshots + return timestamps def load_costs(): @@ -48,22 +59,64 @@ def load_costs(): def load_existing_generators(): generators = pd.read_csv(snakemake.input.generators, - index_col='balancing_authority_code') + index_col=idx_opts[scale]) return generators + +def load_build_years(): + build_years = pd.read_csv(snakemake.input.build_years, + index_col=idx_opts[scale]) + return build_years + + +def linear_growth(init_value, start_year, growth_rate, end_year=2050): + 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') + growth_data = model(years, init_value, start_year, growth_rate) + + growth_df = pd.DataFrame({'demand':growth_data}) + growth_df.index = pd.date_range(start=str(start_year), + periods=(end_year-start_year), + freq='YE') -# attach components + + return growth_df + + +# attach components def attach_load(n): - load = pd.read_csv(snakemake.input.load, parse_dates=True, index_col="Interval End") + initial_demand = float(snakemake.config['total_demand']) - load = load.loc[str(snakemake.config['load_year'])] + if len(model_years) == 1: + demand = [initial_demand] + else: + rate = float(snakemake.config['load_growth']) + start_year = model_years[0] + end_year = model_years[-1] + N_years = end_year-start_year + demand = linear_growth(initial_demand, + model_years[0], + rate + ) + demand = demand.loc[demand.index.year.isin(model_years)].values.flatten() + + + load = pd.read_csv(snakemake.input.load, parse_dates=True, index_col="Interval End") - snapshots = create_snapshots() + start = int(snakemake.config['load_year']) + n_model_years = len(snakemake.config['model_years']) + end = start + n_model_years - 1 + load = load.loc[str(start):str(end)][:int(n_model_years*8760)] + # normalize the load data - load.set_index(snapshots, inplace=True) + for d,y in zip(demand, load.index.year.unique()): + load.loc[str(y)] = load.loc[str(y)].div(load.loc[str(y)].sum(axis=0).sum(),axis=1)*d load = load.resample(f"{snakemake.config['time_res']}h").mean() + load = load[:len(n.snapshots)] + load.set_index(n.snapshots, inplace=True) print('Adding loads to model') for bus in tqdm(n.buses.index): @@ -76,79 +129,194 @@ def attach_load(n): return -def attach_generators(n, costs, generators): - costs = costs.xs((slice(None), slice(None), model_year)) - - available_carriers = costs.index.get_level_values('technology_alias').unique().to_list() - - snapshots = create_snapshots() - wind_profile = pd.read_csv(snakemake.input.wind_profile, parse_dates=True, index_col=0) - wind_profile.set_index(snapshots, inplace=True) - wind_profile = wind_profile.resample(f"{snakemake.config['time_res']}h").mean() +def load_emissions(): + df = pd.read_csv(snakemake.input.emissions, index_col=0) - solar_profile = pd.read_csv(snakemake.input.solar_profile, parse_dates=True, index_col=0) - solar_profile.set_index(snapshots, inplace=True) - solar_profile = solar_profile.resample(f"{snakemake.config['time_res']}h").mean() + return df + + +def add_carriers(n, costs, emissions): + available_carriers = list(costs['technology_alias'].unique()) - # add carriers for carrier in available_carriers: + # emissions data + try: + co2_emissions = emissions.at[carrier,'tCO2 per MWh'] + except KeyError: + co2_emissions = 0.00 + n.add(class_name="Carrier", name=carrier, - color=snakemake.config['carrier_colors'][carrier]) - - # flatten index - costs = costs.reset_index() - costs.loc[costs['technology_alias']=='Solar', 'techdetail'] = 'Utility PV' + co2_emissions=co2_emissions, + color=snakemake.config['carrier_colors'][carrier], + max_growth=float(growth_rates[carrier])) + return + + +def load_re_profile(n, carrier='Solar'): + file_name = {'Solar':snakemake.input.solar_profile, + 'Wind':snakemake.input.wind_profile} + re_profile = pd.read_csv(file_name[carrier], parse_dates=True, index_col=0) + re_profile = re_profile.resample(f"{snakemake.config['time_res']}h").mean().dropna(axis=0) + re_profile.set_index(n.snapshots, inplace=True) + + return re_profile + + + +def attach_renewables(n, costs, generators=None, build_years=None, model_year=None): + carriers = ['Wind','Solar'] + for bus in n.buses.index: for item in costs.itertuples(): tech = item.techdetail carrier = item.technology_alias - capital_cost = item.capital_cost - marginal_cost = item.marginal_cost - lifetime = item.lifetime - - if carrier == 'Batteries': - class_name = "StorageUnit" - max_hours = float(tech.split(' ')[0].strip('Hr')) - cyclic_state_of_charge=True - else: - class_name = "Generator" - max_hours = 0.0 - cyclic_state_of_charge=False - - # existing capacity - if tech in generators.loc[bus].index: - p_nom = generators.at[bus, tech] - else: - p_nom = 0.0 - - # renewables - if carrier == 'Wind': - p_max_pu = p_min_pu = wind_profile[bus] - elif carrier == 'Solar': - p_max_pu = p_min_pu = solar_profile[bus] - else: - p_max_pu = 1 - p_min_pu = 0 - - extendable = tech in snakemake.config['extendable_techs'] - - n.add(class_name=class_name, - name=f"{bus} {tech}", - bus=bus, - p_nom=p_nom, - p_nom_min=p_nom, - p_min_pu=p_min_pu, - p_max_pu=p_max_pu, - p_nom_extendable=extendable, - carrier=carrier, - capital_cost=capital_cost, - marginal_cost=marginal_cost, - lifetime=lifetime, - max_hours=max_hours, - cyclic_state_of_charge=cyclic_state_of_charge) + if carrier in carriers: + re_profile = load_re_profile(n, carrier=carrier) + # existing capacity + if isinstance(generators, pd.DataFrame) and isinstance(build_years, pd.DataFrame): + try: + p_nom = generators.at[bus, tech] + build_year = build_years.at[bus, tech] + except: # the technology is not included in the list of existing generators + continue + name = f"{bus} {tech} EXIST" + extendable = False + 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 + + p_max_pu = re_profile[bus] + + + n.add(class_name="Generator", + name=name, + bus=bus, + p_nom=p_nom, + p_nom_min=p_nom, + p_max_pu=p_max_pu, + p_nom_extendable=extendable, + carrier=carrier, + capital_cost=item.capital_cost, + marginal_cost=item.marginal_cost, + lifetime=item.lifetime, + build_year=build_year) + return + + +def attach_generators(n, costs, generators=None, build_years=None, model_year=None): + carriers = ['Natural Gas','Biomass','Coal','Nuclear','Petroleum'] + for bus in n.buses.index: + for item in costs.itertuples(): + tech = item.techdetail + carrier = item.technology_alias + if carrier in carriers: + # existing capacity + if isinstance(generators, pd.DataFrame) and isinstance(build_years, pd.DataFrame): + try: + p_nom = generators.at[bus, tech] + build_year = build_years.at[bus, tech] + except: # the technology is not included in the list of existing generators + continue + name = f"{bus} {tech} EXIST" + extendable = False + 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 + + # ramp limits + if tech in ['LWR']: + ramp_limit_up = 0.01*resolution + ramp_limit_down = 0.01*resolution + p_nom_min = p_nom + elif tech in ['IGCCAvgCF', 'Biopower']: + ramp_limit_up = 0.1*resolution + ramp_limit_down = 0.1*resolution + p_nom_min = 0 + elif tech in ['CCAvgCF','CTAvgCF']: + ramp_limit_up = min(0.6*resolution,1.0) + ramp_limit_down = min(0.6*resolution,1.0) + p_nom_min = 0 + elif tech in ['Petroleum']: + ramp_limit_up = 1.0 + ramp_limit_down = 1.0 + p_nom_min = 0 + else: + ramp_limit_up = 1.0 + ramp_limit_down = 1.0 + p_nom_min = p_nom + + # minimum/maximum power output + if tech == 'LWR': + p_min_pu = 0.95 + p_max_pu = 1.0 + else: + p_max_pu = 1 + p_min_pu = 0 + + n.add(class_name="Generator", + name=name, + bus=bus, + p_nom=p_nom, + p_nom_min=p_nom_min, + p_nom_extendable=extendable, + p_min_pu = p_min_pu, + p_max_pu = p_max_pu, + carrier=carrier, + capital_cost=item.capital_cost, + marginal_cost=item.marginal_cost, + lifetime=item.lifetime, + ramp_limit_down = ramp_limit_down, + ramp_limit_up = ramp_limit_up, + build_year=build_year, + ) return - + + +def attach_storage(n, costs, generators=None, build_years=None, model_year=None): + carriers = ["Batteries"] + for bus in n.buses.index: + for item in costs.itertuples(): + tech = item.techdetail + carrier = item.technology_alias + if carrier in carriers: + # existing capacity + if isinstance(generators, pd.DataFrame) and isinstance(build_years, pd.DataFrame): + try: + p_nom = generators.at[bus, tech] + build_year = build_years.at[bus, tech] + except: # the technology is not included in the list of existing generators + continue + name = f"{bus} {tech} EXIST" + extendable = False + 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 + + n.add(class_name="StorageUnit", + name=name, + bus=bus, + p_nom=p_nom, + p_nom_min=p_nom, + p_nom_extendable=extendable, + carrier=carrier, + capital_cost=item.capital_cost, + marginal_cost=item.marginal_cost, + lifetime=item.lifetime, + max_hours=float(tech.split(' ')[0].strip('Hr')), + cyclic_state_of_charge=False, + build_year=build_year) + + return if __name__ == "__main__": @@ -157,11 +325,64 @@ def attach_generators(n, costs, generators): 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() + current_costs.loc[current_costs['technology_alias']=='Solar', 'techdetail'] = 'Utility PV' generators = load_existing_generators() + build_years = load_build_years() + emissions = load_emissions() attach_load(n) - attach_generators(n, costs=costs, generators=generators) + add_carriers(n, + costs=current_costs, + emissions=emissions) + + # add existing technology + attach_renewables(n, + costs=current_costs, + generators=generators, + build_years=build_years + ) + attach_generators(n, + costs=current_costs, + generators=generators, + build_years=build_years + ) + attach_storage(n, + costs=current_costs, + generators=generators, + build_years=build_years + ) + + # 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' + attach_renewables(n, + costs=current_costs, + model_year=year + ) + attach_generators(n, + costs=current_costs, + model_year=year + ) + attach_storage(n, + costs=current_costs, + model_year=year + ) + + # add co2 constraint + for y in model_years: + emissions_dict = snakemake.config['co2_limits'] + + n.add(class_name="GlobalConstraint", + name=f"CO2 Limit {y}", + carrier_attribute='co2_emissions', + sense="<=", + investment_period=y, + constant=float(emissions_dict[y])*1e6) n.export_to_netcdf(snakemake.output.elec_network) diff --git a/scripts/build_base_network.py b/scripts/build_base_network.py index 95ff7b8..f1a7112 100644 --- a/scripts/build_base_network.py +++ b/scripts/build_base_network.py @@ -4,6 +4,27 @@ import pypsa +model_years = np.array(snakemake.config['model_years']).astype('int') +resolution = int(snakemake.config['time_res']) + +def add_snapshots(network): + + n_hours = 8760 + + snapshots = pd.DatetimeIndex([]) + for year in model_years: + period = pd.date_range(start=f"{year}-01-01", + freq=f"{resolution}h", + periods=n_hours / resolution) + snapshots = snapshots.append(period) + + network.snapshots = pd.MultiIndex.from_arrays([snapshots.year, snapshots]) + network.investment_periods = model_years + + network.snapshot_weightings.loc[:,:] = resolution + + return + def base_network(): n = pypsa.Network() n.name = 'PyPSA-Illinois' @@ -14,14 +35,7 @@ def base_network(): line_config = snakemake.config['lines'] v_nom = line_config['v_nom'] - resolution = int(snakemake.config['time_res']) - model_year = snakemake.config['model_year'] - n.set_snapshots(pd.date_range(start=f"{model_year}-01-01", - end=f"{model_year+1}-01-01", - inclusive='left', - freq=f"{resolution}h")) - n.snapshot_weightings.loc[:,:] = resolution n.import_components_from_dataframe(buses, 'Bus') n.import_components_from_dataframe(lines, 'Line') @@ -29,4 +43,22 @@ def base_network(): if __name__ == "__main__": n = base_network() + add_snapshots(n) + + year_diff = np.diff(model_years) + if len(year_diff) == 0: + pad = 1 + else: + pad = year_diff[-1] + + n.investment_period_weightings["years"] = list(year_diff) + [pad] + + r = 0.01 + T = 0 + for period, nyears in n.investment_period_weightings.years.items(): + discounts = [(1 / (1 + r) ** t) for t in range(T, T + nyears)] + n.investment_period_weightings.at[period, "objective"] = sum(discounts) + T += nyears + n.investment_period_weightings + n.export_to_netcdf(snakemake.output.network) \ No newline at end of file diff --git a/scripts/plot_results.py b/scripts/plot_results.py new file mode 100644 index 0000000..7a634dc --- /dev/null +++ b/scripts/plot_results.py @@ -0,0 +1,149 @@ +import pypsa +import pandas as pd +import matplotlib.pyplot as plt +import numpy as np + +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 = p_by_carrier[['Nuclear', + 'Coal', + 'Natural Gas', + 'Biomass', + 'Petroleum', + 'Solar', + 'Wind']] + + 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) + + fig, ax = plt.subplots(figsize=(12, 6)) + + color = p_by_carrier.columns.map(n.carriers.color) + + p_by_carrier.where(p_by_carrier > 0).loc[time].plot.area( + ax=ax, + linewidth=0, + color=color, + ) + + charge = p_by_carrier.where(p_by_carrier < 0).dropna(how="all", axis=1).loc[time] + + if not charge.empty: + charge.plot.area( + ax=ax, + linewidth=0, + color=charge.columns.map(n.carriers.color), + ) + + 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) + plt.tight_layout() + + return fig, ax + + +def plot_capacity(n): + + fig, ax = plt.subplots(figsize=(12,8)) + + df = pd.concat([n.generators[['p_nom', 'p_nom_opt']], + n.storage_units[['p_nom','p_nom_opt']]])\ + .replace(0, np.nan)\ + .dropna(axis=0, how='all') + + df.plot.bar(ax=ax) + plt.tight_layout() + + return fig,ax + + +def plot_emissions(n, time_res): + + emissions_data = n.carriers\ + .reset_index()\ + .sort_values(by='Carrier')\ + .set_index('Carrier')[['co2_emissions']]\ + .drop('Batteries') + + emissions = n.generators_t.p * n.generators.carrier.map(n.carriers.co2_emissions) + 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_year = (p_by_carrier.groupby(p_by_carrier.index.get_level_values('period')).sum()) + + annual_emissions = (p_by_carrier_year * emissions_data['co2_emissions']).sum(axis=1).to_frame() + + annual_emissions = annual_emissions * time_res / 1e6 + + annual_emissions.columns = ['Mtonnes CO2/year'] + + fig, ax = plt.subplots(figsize=(12,8)) + + annual_emissions.plot.bar(ax=ax) + + ax.set_ylabel("Mtonnes CO2/year") + + return fig, ax + + +def plot_active_units(n): + + fig, ax = plt.subplots(figsize=(12,8)) + c = "StorageUnit" + df = pd.concat( + { + period: n.get_active_assets(c, period) * n.df(c).p_nom_opt + for period in n.investment_periods + }, + axis=1, + ) + df = df.groupby(n.storage_units.carrier).sum() + + c = "Generator" + df2 = pd.concat( + { + period: n.get_active_assets(c, period) * n.df(c).p_nom_opt + for period in n.investment_periods + }, + axis=1, + ) + df2 = df2.groupby(n.generators.carrier).sum() + df = pd.concat([df, df2]) + df.T.plot.bar( + ax=ax, + stacked=True, + edgecolor="white", + width=1, + ylabel="Capacity", + xlabel="Investment Period", + rot=0, + figsize=(10, 5), + ) + plt.tight_layout() + + return fig, ax + +if __name__ == "__main__": + + n = pypsa.Network(snakemake.input.solved_network) + + fig, ax = plot_dispatch(n, year=int(snakemake.config['plot_year'])) + + plt.savefig(snakemake.output.dispatch_figure) + + fig, ax = plot_capacity(n) + plt.savefig(snakemake.output.capacity_figure) + + fig, ax = plot_emissions(n, float(snakemake.config['time_res'])) + plt.savefig(snakemake.output.emissions_figure) + + fig, ax = plot_active_units(n) + plt.savefig(snakemake.output.active_units_figure) \ No newline at end of file diff --git a/scripts/retrieve_costs.py b/scripts/retrieve_costs.py index d5b1d21..f4a93e1 100644 --- a/scripts/retrieve_costs.py +++ b/scripts/retrieve_costs.py @@ -1,6 +1,21 @@ import pandas as pd from nrelpy.atb import ATBe +n_illinois_reactors = 11 +total_lwr_capacity = 12415.1 +average_reactor_size = total_lwr_capacity/n_illinois_reactors # MW +nuclear_params = snakemake.config['nuclear_params'] +capacity_factor = float(nuclear_params['capacity_factor']) +capital_cost = float(nuclear_params['capital_cost']) # $/MWh +nuclear_fuel = float(nuclear_params['fuel']) # $/MWh +fixed_om = float(nuclear_params['fixed_om']) # $/MWh + +# convert capital and fom to $/MW/yr +operating_hours = 8760*capacity_factor +nuclear_cap_cost = capital_cost*operating_hours/1e3 # $/kW/yr +nuclear_fixed_om = fixed_om*operating_hours/1e3 # $/kW/yr + + if __name__ == "__main__": atb_params = snakemake.config['atb_params'] @@ -85,11 +100,13 @@ # add petroleum df_t = df_pivot.T - # these costs are from 2021 and should be updated to reflect inflation. - df_t['Petroleum','Petroleum',snakemake.config['model_year']] = [1158, - 27.94, - 1.78, - petroleum_price*petroleum_heatrate] # from here https://www.eia.gov/electricity/generatorcosts/ + # [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] + + # these costs are from 2021 and should be updated to reflect inflation. + df_t['Petroleum','Petroleum',year] = [1158, 27.94, 1.78, petroleum_price*petroleum_heatrate] # from here https://www.eia.gov/electricity/generatorcosts/ df_pivot = df_t.T df_pivot.fillna(0., inplace=True) diff --git a/scripts/retrieve_emissions.py b/scripts/retrieve_emissions.py new file mode 100644 index 0000000..ba4da42 --- /dev/null +++ b/scripts/retrieve_emissions.py @@ -0,0 +1,18 @@ +import pandas as pd +from unyt import g, kg, tonne, pound, kWh, MWh, GWh + + +if __name__ == "__main__": + + dfs = pd.read_html("https://www.eia.gov/tools/faqs/faq.php?id=74&t=11") + df = dfs[0].head(3) + df = df.replace({'Natural gas':'Natural Gas'}) + df.columns = df.columns.droplevel() + df = df.set_index(df.iloc[:,0]) + df = df[['pounds per kWh']] + + df = df.iloc[:,0].astype('float').apply(lambda x: (x*(pound/kWh)).to(tonne/MWh)).to_frame() + + df.index.name = '' + df.columns = ['tCO2 per MWh'] + df.to_csv(snakemake.output.emissions) \ No newline at end of file diff --git a/scripts/retrieve_generators.py b/scripts/retrieve_generators.py index be47b5f..a663854 100644 --- a/scripts/retrieve_generators.py +++ b/scripts/retrieve_generators.py @@ -126,10 +126,20 @@ gdf.to_csv(snakemake.output.generators) + gdf['build_year'] = pd.to_datetime(gdf['operating-year-month']).dt.year + scale = snakemake.config['geo_res'] idx_opts = {"rto":"balancing_authority_code", "county":"county"} + + gen_build_year = gdf.pivot_table(index=idx_opts[scale], + columns='technology', + values='build_year', + aggfunc='mean').dropna(axis=1).astype('int') + + gen_build_year.to_csv(snakemake.output.build_year) + gen_agg = gdf.pivot_table(index=idx_opts[scale], columns='technology', values='nameplate-capacity-mw', diff --git a/scripts/retrieve_load.py b/scripts/retrieve_load.py index 7a47997..d114645 100644 --- a/scripts/retrieve_load.py +++ b/scripts/retrieve_load.py @@ -30,6 +30,7 @@ verbose=True) try: + # demand['Interval End'] = demand['Interval End'].dt.tz_convert('US/Central') demand['Interval End'] = demand['Interval End'].dt.tz_localize(None) except: pass @@ -40,7 +41,18 @@ columns='Subregion', values='MW') - demand_pivot[demand_pivot>60e3] = np.nan # this should be replaced with a function to remove outliers. e.g., 2 SD away from mean. + + # calculate outliers with IQR + # threshold = 20 + # q1 = demand_pivot.quantile(0.05) + # q3 = demand_pivot.quantile(0.95) + # IQR = q3-q1 + # upper = (q3+threshold*IQR) + # lower = (q1-threshold*IQR) + # demand_pivot = demand_pivot[~((demand_pivotupper))] + + + demand_pivot[demand_pivot>float(snakemake.config['load_filter'])] = np.nan # this should be replaced with a function to remove outliers demand_pivot = demand_pivot.interpolate("linear") demand_pivot.rename(columns=dict(zip(snakemake.config['rto_subba'], diff --git a/scripts/retrieve_renewables.py b/scripts/retrieve_renewables.py index 1e92c7c..a1043ce 100644 --- a/scripts/retrieve_renewables.py +++ b/scripts/retrieve_renewables.py @@ -3,12 +3,13 @@ import pandas as pd import geopandas as gpd from tqdm import tqdm -from unyt import m, s, MW, W, kg +from unyt import m, s, MW, W, kg, g sys.path.append("functions") from nrel_data_api import parameters, make_csv_url +model_years = np.array(snakemake.config['model_years']).astype('int') def handle_datetime(dataframe): """ @@ -21,10 +22,14 @@ def handle_datetime(dataframe): A pandas dataframe. """ frame = dataframe.copy() - model_year = snakemake.config['model_year'] - timestamps = pd.date_range(f"{model_year}-01-01",f"{model_year+1}-01-01", - freq='1h', inclusive='left') - + timestamps = pd.DatetimeIndex([]) + for year in model_years: + period = pd.date_range(start=f"{year}-01-01", + freq=f"1h", + periods=8760) + timestamps = timestamps.append(period) + + frame.index = timestamps try: frame.set_index(timestamps,inplace=True) frame.drop(columns=['Year','Month','Day','Hour','Minute'],inplace=True) @@ -43,24 +48,32 @@ def retrieve_solar_timeseries(region): region : :class:`gpd.GeoDataFrame` A geopandas dataframe containing modeled bus regions. """ + parameters['attr_list'] = ['ghi'] - parameters['year'] = int(snakemake.config['solar_year']) - pbar = tqdm(region[['name','x','y']].values, position=0, leave=True) - frames = [] - for n, i, j in pbar: - pbar.set_description(f"Processing {n}") - parameters['lon'] = i - parameters['lat'] = j - URL = make_csv_url(parameters=parameters, - kind='solar') - df = pd.read_csv(URL, skiprows=2) - df.rename(columns={'GHI':f"{n}"}, inplace=True) - df = handle_datetime(df) - frames.append(df) - - solar_df = pd.concat(frames, axis=1) - - return solar_df + outer_pbar = tqdm(snakemake.config['solar_years'], position=0, leave=True) + all_frames = [] + for year in outer_pbar: + outer_pbar.set_description(f"Processing {year}") + parameters['year'] = int(year) + frames = [] + inner_pbar = tqdm(region[['name','x','y']].values, position=1, leave=True) + for n, i, j in inner_pbar: + inner_pbar.set_description(f"Processing {n}") + parameters['lon'] = i + parameters['lat'] = j + URL = make_csv_url(parameters=parameters, + kind='solar') + df = pd.read_csv(URL, skiprows=2)[:8760] + df.rename(columns={'GHI':f"{n}"}, inplace=True) + frames.append(df) + + solar_df = pd.concat(frames, axis=1) + + all_frames.append(solar_df) + full_df = pd.concat(all_frames, axis=0) + full_df = handle_datetime(full_df) + + return full_df def retrieve_wind_timeseries(region): @@ -74,23 +87,31 @@ def retrieve_wind_timeseries(region): """ wind_attr = ['windspeed_80m'] parameters['attr_list'] = wind_attr - parameters['year'] = int(snakemake.config['wind_year']) - pbar = tqdm(region[['name','x','y']].values, position=0, leave=True) - frames = [] - for n, i, j in pbar: - pbar.set_description(f"Processing {n}") - parameters['lon'] = i - parameters['lat'] = j - URL = make_csv_url(parameters=parameters, - kind='wind') - df = pd.read_csv(URL, skiprows=1) - df.rename(columns={'wind speed at 80m (m/s)':f"{n}"}, inplace=True) - df = handle_datetime(df) - frames.append(df) - - wind_df = pd.concat(frames, axis=1) - - return wind_df + + outer_pbar = tqdm(snakemake.config['wind_years'], leave=True, position=0) + all_frames = [] + for year in outer_pbar: + outer_pbar.set_description(f"Processing {year}") + + parameters['year'] = int(year) + pbar = tqdm(region[['name','x','y']].values, position=1, leave=True) + frames = [] + for n, i, j in pbar: + pbar.set_description(f"Processing {n}") + parameters['lon'] = i + parameters['lat'] = j + URL = make_csv_url(parameters=parameters, + kind='wind') + df = pd.read_csv(URL, skiprows=1)[:8760] + df.rename(columns={'wind speed at 80m (m/s)':f"{n}"}, inplace=True) + frames.append(df) + + wind_df = pd.concat(frames, axis=1) + all_frames.append(wind_df) + full_df = pd.concat(all_frames, axis=0) + full_df = handle_datetime(full_df) + + return full_df def turbine_power(v): diff --git a/scripts/solve_network.py b/scripts/solve_network.py index f1fa404..7dd7c2e 100644 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -3,58 +3,14 @@ import pypsa import matplotlib.pyplot as plt -def plot_dispatch(n, time="2025-07"): - p_by_carrier = n.generators_t.p.groupby(n.generators.carrier, axis=1).sum().div(1e3) - - p_by_carrier = p_by_carrier[['Nuclear', - 'Coal', - 'Natural Gas', - 'Biomass', - 'Petroleum', - 'Solar', - 'Wind']] - 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) - - fig, ax = plt.subplots(figsize=(6, 3)) - - color = p_by_carrier.columns.map(n.carriers.color) - - p_by_carrier.where(p_by_carrier > 0).loc[time].plot.area( - ax=ax, - linewidth=0, - color=color, - ) - - charge = p_by_carrier.where(p_by_carrier < 0).dropna(how="all", axis=1).loc[time] - - if not charge.empty: - charge.plot.area( - ax=ax, - linewidth=0, - color=charge.columns.map(n.carriers.color), - ) - - n.loads_t.p_set.sum(axis=1).loc[time].div(1e3).plot(ax=ax, c="k") +if __name__ == "__main__": - ax.legend(loc=(1.05, 0)) - ax.set_ylabel("GW") - ax.set_ylim(n.generators_t.p.min().min(), n.loads_t.p_set.sum(axis=1).max()/1e3+2.5) - plt.tight_layout() - - return fig, ax + n = pypsa.Network(snakemake.input.elec_network) + try: + n.optimize(solver_name='cplex') + except: + n.optimize(solver_name='highs') -if __name__ == "__main__": - - n = pypsa.Network(snakemake.input.elec_network) - - n.optimize(solver_name='highs') - - n.export_to_netcdf(snakemake.output.solved_network) - - fig, ax = plot_dispatch(n) - - plt.savefig(snakemake.output.dispatch_figure) \ No newline at end of file + n.export_to_netcdf(snakemake.output.solved_network) \ No newline at end of file