Merge pull request #256 from PyPSA/fix-bug-myopic-co2

Fix bug myopic co2

Snakefile

@@ -464,12 +464,14 @@ rule prepare_sector_network:
         overrides="data/override_component_attrs",
         network=pypsaeur('networks/elec_s{simpl}_{clusters}_ec_lv{lv}_{opts}.nc'),
         energy_totals_name='resources/energy_totals.csv',
+        eurostat=input_eurostat,
         pop_weighted_energy_totals="resources/pop_weighted_energy_totals_s{simpl}_{clusters}.csv",
         transport_demand="resources/transport_demand_s{simpl}_{clusters}.csv",
         transport_data="resources/transport_data_s{simpl}_{clusters}.csv",
         avail_profile="resources/avail_profile_s{simpl}_{clusters}.csv",
         dsm_profile="resources/dsm_profile_s{simpl}_{clusters}.csv",
         co2_totals_name='resources/co2_totals.csv',
+        co2="data/eea/UNFCCC_v23.csv",
         biomass_potentials='resources/biomass_potentials_s{simpl}_{clusters}.csv',
         heat_profile="data/heat_load_profile_BDEW.csv",
         costs=CDIR + "costs_{planning_horizons}.csv",
@@ -568,7 +570,9 @@ rule plot_summary:
     input:
         costs=SDIR + '/csvs/costs.csv',
         energy=SDIR + '/csvs/energy.csv',
-        balances=SDIR + '/csvs/supply_energy.csv'
+        balances=SDIR + '/csvs/supply_energy.csv',
+        eurostat=input_eurostat,
+        country_codes='data/Country_codes.csv',
     output:
         costs=SDIR + '/graphs/costs.pdf',
         energy=SDIR + '/graphs/energy.pdf',

config.default.yaml

@@ -42,7 +42,12 @@ scenario:
   # decay with initial growth rate 0
   planning_horizons: # investment years for myopic and perfect; or costs year for overnight
     - 2030
-  # for example, set to [2020, 2030, 2040, 2050] for myopic foresight
+  # for example, set to
+  # - 2020
+  # - 2030
+  # - 2040
+  # - 2050
+  # for myopic foresight
 
 # CO2 budget as a fraction of 1990 emissions
 # this is over-ridden if CO2Lx is set in sector_opts
@@ -134,7 +139,8 @@ solar_thermal:
 
 # only relevant for foresight = myopic or perfect
 existing_capacities:
-  grouping_years: [1980, 1985, 1990, 1995, 2000, 2005, 2010, 2015, 2020, 2025, 2030]
+  grouping_years_power: [1980, 1985, 1990, 1995, 2000, 2005, 2010, 2015, 2020, 2025, 2030]
+  grouping_years_heat: [1980, 1985, 1990, 1995, 2000, 2005, 2010, 2015, 2019] # these should not extend 2020
   threshold_capacity: 10
   conventional_carriers:
     - lignite

scripts/add_existing_baseyear.py

@@ -131,7 +131,8 @@ def add_power_capacities_installed_before_baseyear(n, grouping_years, costs, bas
         'Oil': 'oil',
         'OCGT': 'OCGT',
         'CCGT': 'CCGT',
-        'Natural Gas': 'gas'
+        'Natural Gas': 'gas',
+        'Bioenergy': 'urban central solid biomass CHP',
     }
 
     fueltype_to_drop = [
@@ -139,7 +140,6 @@ def add_power_capacities_installed_before_baseyear(n, grouping_years, costs, bas
         'Wind',
         'Solar',
         'Geothermal',
-        'Bioenergy',
         'Waste',
         'Other',
         'CCGT, Thermal'
@@ -150,10 +150,29 @@ def add_power_capacities_installed_before_baseyear(n, grouping_years, costs, bas
         'Storage Technologies'
     ]
 
+    # drop unused fueltyps and technologies
     df_agg.drop(df_agg.index[df_agg.Fueltype.isin(fueltype_to_drop)], inplace=True)
     df_agg.drop(df_agg.index[df_agg.Technology.isin(technology_to_drop)], inplace=True)
     df_agg.Fueltype = df_agg.Fueltype.map(rename_fuel)
 
+    # Intermediate fix for DateIn & DateOut
+    # Fill missing DateIn
+    biomass_i = df_agg.loc[df_agg.Fueltype=='urban central solid biomass CHP'].index
+    mean = df_agg.loc[biomass_i, 'DateIn'].mean()
+    df_agg.loc[biomass_i, 'DateIn'] = df_agg.loc[biomass_i, 'DateIn'].fillna(int(mean))
+    # Fill missing DateOut
+    dateout = df_agg.loc[biomass_i, 'DateIn'] + snakemake.config['costs']['lifetime']
+    df_agg.loc[biomass_i, 'DateOut'] = df_agg.loc[biomass_i, 'DateOut'].fillna(dateout)
+
+
+    # drop assets which are already phased out / decomissioned
+    phased_out = df_agg[df_agg["DateOut"]<baseyear].index
+    df_agg.drop(phased_out, inplace=True)
+
+    # calculate remaining lifetime before phase-out (+1 because assumming
+    # phase out date at the end of the year)
+    df_agg["lifetime"] = df_agg.DateOut - df_agg.DateIn + 1
+
     # assign clustered bus
     busmap_s = pd.read_csv(snakemake.input.busmap_s, index_col=0).squeeze()
     busmap = pd.read_csv(snakemake.input.busmap, index_col=0).squeeze()
@@ -182,35 +201,52 @@ def add_power_capacities_installed_before_baseyear(n, grouping_years, costs, bas
         aggfunc='sum'
     )
 
+    lifetime = df_agg.pivot_table(
+        index=["grouping_year", 'Fueltype'],
+        columns='cluster_bus',
+        values='lifetime',
+        aggfunc='mean'  # currently taken mean for clustering lifetimes
+    )
+
     carrier = {
         "OCGT": "gas",
         "CCGT": "gas",
         "coal": "coal",
         "oil": "oil",
         "lignite": "lignite",
-        "nuclear": "uranium"
+        "nuclear": "uranium",
+        'urban central solid biomass CHP': "biomass",
     }
 
     for grouping_year, generator in df.index:
 
+
         # capacity is the capacity in MW at each node for this
         capacity = df.loc[grouping_year, generator]
         capacity = capacity[~capacity.isna()]
         capacity = capacity[capacity > snakemake.config['existing_capacities']['threshold_capacity']]
-
-        if generator in ['solar', 'onwind', 'offwind']:
-
         suffix = '-ac' if generator == 'offwind' else ''
-            name_suffix = f' {generator}{suffix}-{baseyear}'
+        name_suffix = f' {generator}{suffix}-{grouping_year}'
+        asset_i = capacity.index + name_suffix
+        if generator in ['solar', 'onwind', 'offwind']:
 
             # to consider electricity grid connection costs or a split between
             # solar utility and rooftop as well, rather take cost assumptions
             # from existing network than from the cost database
             capital_cost = n.generators.loc[n.generators.carrier==generator+suffix, "capital_cost"].mean()
 
+            # check if assets are already in network (e.g. for 2020)
+            already_build = n.generators.index.intersection(asset_i)
+            new_build = asset_i.difference(n.generators.index)
+
+            # this is for the year 2020
+            if not already_build.empty:
+                n.generators.loc[already_build, "p_nom_min"] = capacity.loc[already_build.str.replace(name_suffix, "")].values
+            new_capacity = capacity.loc[new_build.str.replace(name_suffix, "")]
+
             if 'm' in snakemake.wildcards.clusters:
 
-                for ind in capacity.index:
+                for ind in new_capacity.index:
 
                     # existing capacities are split evenly among regions in every country
                     inv_ind = [i for i in inv_busmap[ind]]
@@ -225,7 +261,7 @@ def add_power_capacities_installed_before_baseyear(n, grouping_years, costs, bas
                         [i + name_suffix for i in inv_ind],
                         bus=ind,
                         carrier=generator,
-                        p_nom=capacity[ind] / len(inv_ind), # split among regions in a country
+                        p_nom=new_capacity[ind] / len(inv_ind), # split among regions in a country
                         marginal_cost=costs.at[generator,'VOM'],
                         capital_cost=capital_cost,
                         efficiency=costs.at[generator, 'efficiency'],
@@ -236,14 +272,15 @@ def add_power_capacities_installed_before_baseyear(n, grouping_years, costs, bas
 
             else:
 
-                p_max_pu = n.generators_t.p_max_pu[capacity.index + name_suffix]
+                p_max_pu = n.generators_t.p_max_pu[capacity.index + f' {generator}{suffix}-{baseyear}']
 
+                if not new_build.empty:
                     n.madd("Generator",
-                    capacity.index,
+                        new_capacity.index,
-                    suffix=' ' + generator +"-"+ str(grouping_year),
+                        suffix=' ' + name_suffix,
-                    bus=capacity.index,
+                        bus=new_capacity.index,
                         carrier=generator,
-                    p_nom=capacity,
+                        p_nom=new_capacity,
                         marginal_cost=costs.at[generator, 'VOM'],
                         capital_cost=capital_cost,
                         efficiency=costs.at[generator, 'efficiency'],
@@ -257,20 +294,49 @@ def add_power_capacities_installed_before_baseyear(n, grouping_years, costs, bas
             if "EU" not in vars(spatial)[carrier[generator]].locations:
                 bus0 = bus0.intersection(capacity.index + " gas")
 
+            already_build = n.links.index.intersection(asset_i)
+            new_build = asset_i.difference(n.links.index)
+            lifetime_assets = lifetime.loc[grouping_year,generator].dropna()
+
+            # this is for the year 2020
+            if not already_build.empty:
+                n.links.loc[already_build, "p_nom_min"] = capacity.loc[already_build.str.replace(name_suffix, "")].values
+
+            if not new_build.empty:
+                new_capacity = capacity.loc[new_build.str.replace(name_suffix, "")]
+
+                if generator!="urban central solid biomass CHP":
                     n.madd("Link",
-                capacity.index,
+                        new_capacity.index,
-                suffix= " " + generator +"-" + str(grouping_year),
+                        suffix= name_suffix,
                         bus0=bus0,
-                bus1=capacity.index,
+                        bus1=new_capacity.index,
                         bus2="co2 atmosphere",
                         carrier=generator,
                         marginal_cost=costs.at[generator, 'efficiency'] * costs.at[generator, 'VOM'], #NB: VOM is per MWel
                         capital_cost=costs.at[generator, 'efficiency'] * costs.at[generator, 'fixed'], #NB: fixed cost is per MWel
-                p_nom=capacity / costs.at[generator, 'efficiency'],
+                        p_nom=new_capacity / costs.at[generator, 'efficiency'],
                         efficiency=costs.at[generator, 'efficiency'],
                         efficiency2=costs.at[carrier[generator], 'CO2 intensity'],
                         build_year=grouping_year,
-                lifetime=costs.at[generator, 'lifetime']
+                        lifetime=lifetime_assets.loc[new_capacity.index],
+                    )
+                else:
+                    key = 'central solid biomass CHP'
+                    n.madd("Link",
+                        new_capacity.index,
+                        suffix= name_suffix,
+                        bus0=spatial.biomass.df.loc[new_capacity.index]["nodes"].values,
+                        bus1=new_capacity.index,
+                        bus2=new_capacity.index + " urban central heat",
+                        carrier=generator,
+                        p_nom=new_capacity / costs.at[key, 'efficiency'],
+                        capital_cost=costs.at[key, 'fixed'] * costs.at[key, 'efficiency'],
+                        marginal_cost=costs.at[key, 'VOM'],
+                        efficiency=costs.at[key, 'efficiency'],
+                        build_year=grouping_year,
+                        efficiency2=costs.at[key, 'efficiency-heat'],
+                        lifetime=lifetime_assets.loc[new_capacity.index]
                     )
 
 
@@ -376,8 +442,8 @@ def add_heating_capacities_installed_before_baseyear(n, baseyear, grouping_years
         for i, grouping_year in enumerate(grouping_years):
 
             if int(grouping_year) + default_lifetime <= int(baseyear):
-                ratio = 0
+                continue
-            else:
+
             # installation is assumed to be linear for the past 25 years (default lifetime)
             ratio = (int(grouping_year) - int(grouping_years[i-1])) / default_lifetime
 
@@ -443,7 +509,7 @@ def add_heating_capacities_installed_before_baseyear(n, baseyear, grouping_years
             # delete links with p_nom=nan corresponding to extra nodes in country
             n.mremove("Link", [index for index in n.links.index.to_list() if str(grouping_year) in index and np.isnan(n.links.p_nom[index])])
 
-            # delete links if their lifetime is over and p_nom=0
+            # delete links with capacities below threshold
             threshold = snakemake.config['existing_capacities']['threshold_capacity']
             n.mremove("Link", [index for index in n.links.index.to_list() if str(grouping_year) in index and n.links.p_nom[index] < threshold])
 
@@ -454,11 +520,11 @@ if __name__ == "__main__":
         snakemake = mock_snakemake(
             'add_existing_baseyear',
             simpl='',
-            clusters="37",
+            clusters="45",
             lv=1.0,
             opts='',
-            sector_opts='168H-T-H-B-I-solar+p3-dist1',
+            sector_opts='365H-T-H-B-I-A-solar+p3-dist1',
-            planning_horizons=2020,
+            planning_horizons=2030,
         )
 
     logging.basicConfig(level=snakemake.config['logging_level'])
@@ -485,15 +551,17 @@ if __name__ == "__main__":
         snakemake.config['costs']['lifetime']
     )
 
-    grouping_years = snakemake.config['existing_capacities']['grouping_years']
+    grouping_years_power = snakemake.config['existing_capacities']['grouping_years_power']
-    add_power_capacities_installed_before_baseyear(n, grouping_years, costs, baseyear)
+    grouping_years_heat = snakemake.config['existing_capacities']['grouping_years_heat']
+    add_power_capacities_installed_before_baseyear(n, grouping_years_power, costs, baseyear)
 
     if "H" in opts:
         time_dep_hp_cop = options["time_dep_hp_cop"]
         ashp_cop = xr.open_dataarray(snakemake.input.cop_air_total).to_pandas().reindex(index=n.snapshots)
         gshp_cop = xr.open_dataarray(snakemake.input.cop_soil_total).to_pandas().reindex(index=n.snapshots)
         default_lifetime = snakemake.config['costs']['lifetime']
-        add_heating_capacities_installed_before_baseyear(n, baseyear, grouping_years, ashp_cop, gshp_cop, time_dep_hp_cop, costs, default_lifetime)
+        add_heating_capacities_installed_before_baseyear(n, baseyear, grouping_years_heat,
+                                                         ashp_cop, gshp_cop, time_dep_hp_cop, costs, default_lifetime)
 
     n.meta = dict(snakemake.config, **dict(wildcards=dict(snakemake.wildcards)))
     n.export_to_netcdf(snakemake.output[0])

scripts/build_energy_totals.py

@@ -127,17 +127,16 @@ to_ipcc = {
 }
 
 
-def build_eurostat(countries, year):
+def build_eurostat(input_eurostat, countries, report_year,  year):
     """Return multi-index for all countries' energy data in TWh/a."""
 
-    report_year = snakemake.config["energy"]["eurostat_report_year"]
     filenames = {
         2016: f"/{year}-Energy-Balances-June2016edition.xlsx",
         2017: f"/{year}-ENERGY-BALANCES-June2017edition.xlsx"
     }
 
     dfs = pd.read_excel(
-        snakemake.input.eurostat + filenames[report_year],
+        input_eurostat + filenames[report_year],
         sheet_name=None,
         skiprows=1,
         index_col=list(range(4)),
@@ -563,18 +562,18 @@ def build_energy_totals(countries, eurostat, swiss, idees):
     return df
 
 
-def build_eea_co2(year=1990):
+def build_eea_co2(input_co2, year=1990, emissions_scope="CO2"):
 
     # https://www.eea.europa.eu/data-and-maps/data/national-emissions-reported-to-the-unfccc-and-to-the-eu-greenhouse-gas-monitoring-mechanism-16
     # downloaded 201228 (modified by EEA last on 201221)
-    df = pd.read_csv(snakemake.input.co2, encoding="latin-1")
+    df = pd.read_csv(input_co2, encoding="latin-1")
 
     df.replace(dict(Year="1985-1987"), 1986, inplace=True)
     df.Year = df.Year.astype(int)
     index_col = ["Country_code", "Pollutant_name", "Year", "Sector_name"]
     df = df.set_index(index_col).sort_index()
 
-    emissions_scope = snakemake.config["energy"]["emissions"]
+    emissions_scope = emissions_scope
 
     cts = ["CH", "EUA", "NO"] + eu28_eea
 
@@ -611,9 +610,9 @@ def build_eea_co2(year=1990):
     return emissions / 1e3
 
 
-def build_eurostat_co2(countries, year=1990):
+def build_eurostat_co2(input_eurostat, countries, report_year, year=1990):
 
-    eurostat = build_eurostat(countries, year)
+    eurostat = build_eurostat(input_eurostat, countries, report_year, year)
 
     specific_emissions = pd.Series(index=eurostat.columns, dtype=float)
 
@@ -702,7 +701,9 @@ if __name__ == "__main__":
     idees_countries = countries.intersection(eu28)
 
     data_year = config["energy_totals_year"]
-    eurostat = build_eurostat(countries, data_year)
+    report_year = snakemake.config["energy"]["eurostat_report_year"]
+    input_eurostat = snakemake.input.eurostat
+    eurostat = build_eurostat(input_eurostat, countries, report_year, data_year)
     swiss = build_swiss(data_year)
     idees = build_idees(idees_countries, data_year)
 
@@ -710,8 +711,9 @@ if __name__ == "__main__":
     energy.to_csv(snakemake.output.energy_name)
 
     base_year_emissions = config["base_emissions_year"]
-    eea_co2 = build_eea_co2(base_year_emissions)
+    emissions_scope = snakemake.config["energy"]["emissions"]
-    eurostat_co2 = build_eurostat_co2(countries, base_year_emissions)
+    eea_co2 = build_eea_co2(snakemake.input.co2, base_year_emissions, emissions_scope)
+    eurostat_co2 = build_eurostat_co2(input_eurostat, countries, report_year, base_year_emissions)
 
     co2 = build_co2_totals(countries, eea_co2, eurostat_co2)
     co2.to_csv(snakemake.output.co2_name)

scripts/plot_summary.py

@@ -363,7 +363,7 @@ def historical_emissions(cts):
 
 
 
-def plot_carbon_budget_distribution():
+def plot_carbon_budget_distribution(input_eurostat):
     """
     Plot historical carbon emissions in the EU and decarbonization path
     """
@@ -385,9 +385,9 @@ def plot_carbon_budget_distribution():
     ax1.set_xlim([1990,snakemake.config['scenario']['planning_horizons'][-1]+1])
 
     path_cb = snakemake.config['results_dir'] + snakemake.config['run'] + '/csvs/'
-    countries=pd.read_csv(path_cb + 'countries.csv',  index_col=1)
+    countries = pd.read_csv(snakemake.input.country_codes, index_col=1)
     cts = countries.index.to_list()
-    e_1990 = co2_emissions_year(cts, opts, year=1990)
+    e_1990 = co2_emissions_year(cts, input_eurostat, opts, year=1990)
     CO2_CAP=pd.read_csv(path_cb + 'carbon_budget_distribution.csv',
                         index_col=0)
 
@@ -439,7 +439,6 @@ if __name__ == "__main__":
         from helper import mock_snakemake
         snakemake = mock_snakemake('plot_summary')
 
-    update_config_with_sector_opts(snakemake.config, snakemake.wildcards.sector_opts)
 
     n_header = 4
 
@@ -453,4 +452,4 @@ if __name__ == "__main__":
         opts=sector_opts.split('-')
         for o in opts:
             if "cb" in o:
-                plot_carbon_budget_distribution()
+                plot_carbon_budget_distribution(snakemake.input.eurostat)

scripts/prepare_sector_network.py

@@ -158,21 +158,22 @@ def get(item, investment_year=None):
         return item
 
 
-def co2_emissions_year(countries, opts, year):
+def co2_emissions_year(countries, input_eurostat, opts, emissions_scope, report_year, year):
     """
     Calculate CO2 emissions in one specific year (e.g. 1990 or 2018).
     """
-
+    emissions_scope = snakemake.config["energy"]["emissions"]
-    eea_co2 = build_eea_co2(year)
+    eea_co2 = build_eea_co2(snakemake.input.co2, year, emissions_scope)
 
     # TODO: read Eurostat data from year > 2014
     # this only affects the estimation of CO2 emissions for BA, RS, AL, ME, MK
+    report_year = snakemake.config["energy"]["eurostat_report_year"]
     if year > 2014:
-        eurostat_co2 = build_eurostat_co2(year=2014)
+        eurostat_co2 = build_eurostat_co2(input_eurostat, countries, report_year, year=2014)
     else:
-        eurostat_co2 = build_eurostat_co2(year)
+        eurostat_co2 = build_eurostat_co2(input_eurostat, countries, report_year, year)
 
-    co2_totals = build_co2_totals(eea_co2, eurostat_co2)
+    co2_totals = build_co2_totals(countries, eea_co2, eurostat_co2)
 
     sectors = emission_sectors_from_opts(opts)
 
@@ -185,7 +186,7 @@ def co2_emissions_year(countries, opts, year):
 
 
 # TODO: move to own rule with sector-opts wildcard?
-def build_carbon_budget(o, fn):
+def build_carbon_budget(o, input_eurostat, fn, emissions_scope, report_year):
     """
     Distribute carbon budget following beta or exponential transition path.
     """
@@ -202,10 +203,12 @@ def build_carbon_budget(o, fn):
 
     countries = n.buses.country.dropna().unique()
 
-    e_1990 = co2_emissions_year(countries, opts, year=1990)
+    e_1990 = co2_emissions_year(countries, input_eurostat, opts, emissions_scope,
+                                report_year, year=1990)
 
     #emissions at the beginning of the path (last year available 2018)
-    e_0 = co2_emissions_year(countries, opts, year=2018)
+    e_0 = co2_emissions_year(countries, input_eurostat, opts, emissions_scope,
+                             report_year,year=2018)
 
     planning_horizons = snakemake.config['scenario']['planning_horizons']
     t_0 = planning_horizons[0]
@@ -233,8 +236,9 @@ def build_carbon_budget(o, fn):
         co2_cap = pd.Series({t: exponential_decay(t) for t in planning_horizons}, name=o)
 
     # TODO log in Snakefile
-    if not os.path.exists(fn):
+    csvs_folder = fn.rsplit("/", 1)[0]
-        os.makedirs(fn)
+    if not os.path.exists(csvs_folder):
+        os.makedirs(csvs_folder)
     co2_cap.to_csv(fn, float_format='%.3f')
 
 
@@ -2423,7 +2427,7 @@ if __name__ == "__main__":
             opts="",
             clusters="37",
             lv=1.5,
-            sector_opts='Co2L0-168H-T-H-B-I-solar3-dist1',
+            sector_opts='cb40ex0-365H-T-H-B-I-A-solar+p3-dist1',
             planning_horizons="2020",
         )
 
@@ -2527,9 +2531,11 @@ if __name__ == "__main__":
         limit_type = "carbon budget"
         fn = snakemake.config['results_dir'] + snakemake.config['run'] + '/csvs/carbon_budget_distribution.csv'
         if not os.path.exists(fn):
-            build_carbon_budget(o, fn)
+            emissions_scope = snakemake.config["energy"]["emissions"]
+            report_year = snakemake.config["energy"]["eurostat_report_year"]
+            build_carbon_budget(o, snakemake.input.eurostat, fn, emissions_scope, report_year)
         co2_cap = pd.read_csv(fn, index_col=0).squeeze()
-        limit = co2_cap[investment_year]
+        limit = co2_cap.loc[investment_year]
         break
     for o in opts:
         if not "Co2L" in o: continue