From 9b9090c76cb1fbd601626342ce569d87e490d9d0 Mon Sep 17 00:00:00 2001 From: chrstphtrs Date: Wed, 18 Oct 2023 16:59:49 +0200 Subject: [PATCH 01/23] add option for additional national carbon budget constraints --- config/config.default.yaml | 17 +++++++ rules/solve_myopic.smk | 2 + scripts/solve_network.py | 99 ++++++++++++++++++++++++++++++++++++++ 3 files changed, 118 insertions(+) diff --git a/config/config.default.yaml b/config/config.default.yaml index 7dc0cf76f..325bbbaa2 100644 --- a/config/config.default.yaml +++ b/config/config.default.yaml @@ -84,6 +84,22 @@ co2_budget: 2045: 0.032 2050: 0.000 +co2_budget_national: + 2030: + 'DE': 0.350 + 'AT': 0.450 + 'BE': 0.450 + 'CH': 0.450 + 'CZ': 0.450 + 'DK': 0.450 + 'FR': 0.450 + 'GB': 0.450 + 'LU': 0.450 + 'NL': 0.450 + 'NO': 0.450 + 'PL': 0.450 + 'SE': 0.450 + # docs in https://pypsa-eur.readthedocs.io/en/latest/configuration.html#electricity electricity: voltages: [220., 300., 380.] @@ -454,6 +470,7 @@ sector: hydrogen_turbine: false SMR: true SMR_cc: true + co2_budget_national: false regional_co2_sequestration_potential: enable: false attribute: 'conservative estimate Mt' diff --git a/rules/solve_myopic.smk b/rules/solve_myopic.smk index 8a93d24a8..06fd9b79a 100644 --- a/rules/solve_myopic.smk +++ b/rules/solve_myopic.smk @@ -88,11 +88,13 @@ rule solve_sector_network_myopic: co2_sequestration_potential=config["sector"].get( "co2_sequestration_potential", 200 ), + countries=config["countries"], input: network=RESULTS + "prenetworks-brownfield/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_{planning_horizons}.nc", costs="data/costs_{planning_horizons}.csv", config=RESULTS + "config.yaml", + co2_totals_name=RESOURCES + "co2_totals.csv", output: RESULTS + "postnetworks/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_{planning_horizons}.nc", diff --git a/scripts/solve_network.py b/scripts/solve_network.py index 224d47147..f5dd79e0b 100644 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -41,6 +41,8 @@ pypsa.pf.logger.setLevel(logging.WARNING) from pypsa.descriptors import get_switchable_as_dense as get_as_dense +from prepare_sector_network import emission_sectors_from_opts + def add_land_use_constraint(n, planning_horizons, config): if "m" in snakemake.wildcards.clusters: @@ -762,6 +764,92 @@ def add_pipe_retrofit_constraint(n): n.model.add_constraints(lhs == rhs, name="Link-pipe_retrofit") +def add_co2limit_country(n, config, limit_countries, nyears=1.0): + """ + Add a set of emissions limit constraints for specified countries. + + The countries and emissions limits are specified in the config file entry 'co2_budget_country_{investment_year}'. + + Parameters + ---------- + n : pypsa.Network + config : dict + limit_countries : dict + nyears: float, optional + Used to scale the emissions constraint to the number of snapshots of the base network. + """ + logger.info(f"Adding CO2 budget limit for each country as per unit of 1990 levels") + + # TODO: n.config (submodule) vs snakemake.config (main module, overwrite/overwritten config)? + # countries = config.countries + # print(config) + countries = ['AT', 'BE', 'CH', 'CZ', 'DE', 'DK', 'FR', 'GB', 'LU', 'NL', 'NO', 'PL', 'SE'] + + # TODO: import function from prepare_sector_network? Move to common place? + sectors = emission_sectors_from_opts(opts) + + # convert Mt to tCO2 + co2_totals = 1e6 * pd.read_csv(snakemake.input.co2_totals_name, index_col=0) + + co2_limit_countries = co2_totals.loc[countries, sectors].sum(axis=1) + co2_limit_countries = co2_limit_countries.loc[co2_limit_countries.index.isin(limit_countries.keys())] + + co2_limit_countries *= co2_limit_countries.index.map(limit_countries) * nyears + + p = n.model["Link-p"] # dimension: (time, component) + + # NB: Most country-specific links retain their locational information in bus1 (except for DAC, where it is in bus2) + country = n.links.bus1.map(n.buses.location).map(n.buses.country) + country_DAC = ( + n.links[n.links.carrier == "DAC"] + .bus2.map(n.buses.location) + .map(n.buses.country) + ) + country[country_DAC.index] = country_DAC + + lhs = [] + for port in [col[3:] for col in n.links if col.startswith("bus")]: + if port == str(0): + efficiency = ( + n.links["efficiency"].apply(lambda x: 1.0).rename("efficiency0") + ) + elif port == str(1): + efficiency = n.links["efficiency"].rename("efficiency1") + else: + efficiency = n.links[f"efficiency{port}"] + mask = n.links[f"bus{port}"].map(n.buses.carrier).eq("co2") + + idx = n.links[mask].index + + grouping = country.loc[idx] + + if not grouping.isnull().all(): + expr = ( + (p.loc[:, idx] * efficiency[idx]) + .groupby(grouping, axis=1) + .sum() + .sum(dims="snapshot") + ) + lhs.append(expr) + + lhs = sum(lhs) # dimension: (country) + lhs = lhs.rename({list(lhs.dims.keys())[0]: "country"}) + rhs = pd.Series(co2_limit_countries) # dimension: (country) + + for ct in lhs.indexes["country"]: + n.model.add_constraints( + lhs.loc[ct] <= rhs[ct], + name=f"GlobalConstraint-co2_limit_per_country{ct}", + ) + n.add( + "GlobalConstraint", + f"co2_limit_per_country{ct}", + constant=rhs[ct], + sense="<=", + type="", + ) + + def extra_functionality(n, snapshots): """ Collects supplementary constraints which will be passed to @@ -792,6 +880,17 @@ def extra_functionality(n, snapshots): add_carbon_budget_constraint(n, snapshots) add_retrofit_gas_boiler_constraint(n, snapshots) + if n.config["sector"]["co2_budget_national"]: + # prepare co2 constraint + nhours = n.snapshot_weightings.generators.sum() + nyears = nhours / 8760 + investment_year = int(snakemake.wildcards.planning_horizons[-4:]) + limit_countries = snakemake.config["co2_budget_national"][investment_year] + + # add co2 constraint for each country + logger.info(f"Add CO2 limit for each country") + add_co2limit_country(n, config, limit_countries, nyears) + def solve_network(n, config, solving, opts="", **kwargs): set_of_options = solving["solver"]["options"] From a35f5479aedd933773634a931e407d0535a6da64 Mon Sep 17 00:00:00 2001 From: chrstphtrs Date: Tue, 24 Oct 2023 14:06:17 +0200 Subject: [PATCH 02/23] add links instead of equal-and-opposite fuel/emissions load pairs for land transport oil (ICEs), naphtha for industry and kerosene for aviation (before summed as 'oil'), shipping oil, shipping methanol, agriculture machinery oil --- scripts/prepare_sector_network.py | 246 +++++++++++++++++++----------- 1 file changed, 160 insertions(+), 86 deletions(-) diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py index ee2f0e3c6..989bdb780 100644 --- a/scripts/prepare_sector_network.py +++ b/scripts/prepare_sector_network.py @@ -135,6 +135,7 @@ def define_spatial(nodes, options): spatial.oil = SimpleNamespace() spatial.oil.nodes = ["EU oil"] spatial.oil.locations = ["EU"] + spatial.oil.land_transport = nodes + " land transport oil" # uranium spatial.uranium = SimpleNamespace() @@ -1467,8 +1468,8 @@ def add_land_transport(n, costs): n.madd( "Bus", nodes, - location=nodes, suffix=" EV battery", + location=nodes, carrier="Li ion", unit="MWh_el", ) @@ -1568,28 +1569,31 @@ def add_land_transport(n, costs): ice_efficiency = options["transport_internal_combustion_efficiency"] n.madd( - "Load", - nodes, - suffix=" land transport oil", - bus=spatial.oil.nodes, + "Bus", + spatial.oil.land_transport, + location=nodes, carrier="land transport oil", - p_set=ice_share / ice_efficiency * transport[nodes], + unit="land transport", ) - co2 = ( - ice_share - / ice_efficiency - * transport[nodes].sum().sum() - / nhours - * costs.at["oil", "CO2 intensity"] + n.madd( + "Load", + spatial.oil.land_transport, + bus=spatial.oil.land_transport, + carrier="land transport oil", + p_set=ice_share / ice_efficiency * transport[nodes].rename(columns=lambda x: x + " land transport oil"), ) - n.add( - "Load", - "land transport oil emissions", - bus="co2 atmosphere", - carrier="land transport oil emissions", - p_set=-co2, + n.madd( + "Link", + spatial.oil.land_transport, + bus0=spatial.oil.nodes, + bus1=spatial.oil.land_transport, + bus2="co2 atmosphere", + carrier="land transport oil", + efficiency=ice_efficiency, + efficiency2=costs.at["oil", "CO2 intensity"], + p_nom_extendable=True, ) @@ -2611,46 +2615,36 @@ def add_industry(n, costs): ) p_set_methanol = shipping_methanol_share * p_set.sum() * efficiency - n.madd( - "Load", - spatial.methanol.nodes, - suffix=" shipping methanol", - bus=spatial.methanol.nodes, + n.add( + "Bus", + "EU shipping methanol", + location="EU", carrier="shipping methanol", - p_set=p_set_methanol, + unit="MWh_LHV", ) - # CO2 intensity methanol based on stoichiometric calculation with 22.7 GJ/t methanol (32 g/mol), CO2 (44 g/mol), 277.78 MWh/TJ = 0.218 t/MWh - co2 = p_set_methanol / options["MWh_MeOH_per_tCO2"] - n.add( "Load", - "shipping methanol emissions", - bus="co2 atmosphere", - carrier="shipping methanol emissions", - p_set=-co2, + "shipping methanol", + bus="EU shipping methanol", + carrier="shipping methanol", + p_set=p_set_methanol, ) - if shipping_oil_share: - p_set_oil = shipping_oil_share * p_set.sum() + if len(spatial.methanol.nodes) == 1: + link_names = ["EU shipping methanol"] + else: + link_names = nodes + " shipping methanol" n.madd( - "Load", - spatial.oil.nodes, - suffix=" shipping oil", - bus=spatial.oil.nodes, - carrier="shipping oil", - p_set=p_set_oil, - ) - - co2 = p_set_oil * costs.at["oil", "CO2 intensity"] - - n.add( - "Load", - "shipping oil emissions", - bus="co2 atmosphere", - carrier="shipping oil emissions", - p_set=-co2, + "Link", + link_names, + bus0=spatial.methanol.nodes, + bus1="EU shipping methanol", + bus2="co2 atmosphere", + carrier="shipping methanol", + p_nom_extendable=True, + efficiency2=1 / options["MWh_MeOH_per_tCO2"], # CO2 intensity methanol based on stoichiometric calculation with 22.7 GJ/t methanol (32 g/mol), CO2 (44 g/mol), 277.78 MWh/TJ = 0.218 t/MWh ) if "oil" not in n.buses.carrier.unique(): @@ -2683,6 +2677,41 @@ def add_industry(n, costs): marginal_cost=costs.at["oil", "fuel"], ) + if shipping_oil_share: + p_set_oil = shipping_oil_share * p_set.sum() + + n.add( + "Bus", + "EU shipping oil", + location="EU", + carrier="shipping oil", + unit="MWh_LHV", + ) + + n.add( + "Load", + "shipping oil", + bus="EU shipping oil", + carrier="shipping oil", + p_set=p_set_oil, + ) + + if len(spatial.oil.nodes) == 1: + link_names = ["EU shipping oil"] + else: + link_names = nodes + " shipping oil" + + n.madd( + "Link", + link_names, + bus0=spatial.oil.nodes, + bus1="EU shipping oil", + bus2="co2 atmosphere", + carrier="shipping oil", + p_nom_extendable=True, + efficiency2=costs.at["oil", "CO2 intensity"], + ) + if options["oil_boilers"]: nodes_heat = create_nodes_for_heat_sector()[0] @@ -2724,19 +2753,49 @@ def add_industry(n, costs): lifetime=costs.at["Fischer-Tropsch", "lifetime"], ) + # naphtha demand_factor = options.get("HVC_demand_factor", 1) - p_set = demand_factor * industrial_demand.loc[nodes, "naphtha"].sum() / nhours if demand_factor != 1: logger.warning(f"Changing HVC demand by {demand_factor*100-100:+.2f}%.") - n.madd( + # NB: CO2 gets released again to atmosphere when plastics decay + # except for the process emissions when naphtha is used for petrochemicals, which can be captured with other industry process emissions + # convert process emissions from feedstock from MtCO2 to energy demand + p_set = demand_factor * (industrial_demand.loc[nodes, "naphtha"] - industrial_demand.loc[nodes, "process emission from feedstock"] / costs.at["oil", "CO2 intensity"]).sum() / nhours + + n.add( + "Bus", + "EU naphtha for industry", + location="EU", + carrier="naphtha for industry", + unit="MWh_LHV", + ) + + n.add( "Load", - ["naphtha for industry"], - bus=spatial.oil.nodes, + "naphtha for industry", + bus="EU naphtha for industry", carrier="naphtha for industry", p_set=p_set, ) + if len(spatial.oil.nodes) == 1: + link_names = ["EU naphtha for industry"] + else: + link_names = nodes + " naphtha for industry" + + n.madd( + "Link", + link_names, + bus0=spatial.oil.nodes, + bus1="EU naphtha for industry", + bus2="co2 atmosphere", + carrier="naphtha for industry", + p_nom_extendable=True, + efficiency2=costs.at["oil", "CO2 intensity"], + ) + + # aviation demand_factor = options.get("aviation_demand_factor", 1) all_aviation = ["total international aviation", "total domestic aviation"] p_set = ( @@ -2748,29 +2807,36 @@ def add_industry(n, costs): if demand_factor != 1: logger.warning(f"Changing aviation demand by {demand_factor*100-100:+.2f}%.") - n.madd( + n.add( + "Bus", + "EU kerosene for aviation", + location="EU", + carrier="kerosene for aviation", + unit="MWh_LHV", + ) + + n.add( "Load", - ["kerosene for aviation"], - bus=spatial.oil.nodes, + "kerosene for aviation", + bus="EU kerosene for aviation", carrier="kerosene for aviation", p_set=p_set, ) - # NB: CO2 gets released again to atmosphere when plastics decay or kerosene is burned - # except for the process emissions when naphtha is used for petrochemicals, which can be captured with other industry process emissions - # tco2 per hour - co2_release = ["naphtha for industry", "kerosene for aviation"] - co2 = ( - n.loads.loc[co2_release, "p_set"].sum() * costs.at["oil", "CO2 intensity"] - - industrial_demand.loc[nodes, "process emission from feedstock"].sum() / nhours - ) + if len(spatial.oil.nodes) == 1: + link_names = ["EU kerosene for aviation"] + else: + link_names = nodes + " kerosene for aviation" - n.add( - "Load", - "oil emissions", - bus="co2 atmosphere", - carrier="oil emissions", - p_set=-co2, + n.madd( + "Link", + link_names, + bus0=spatial.oil.nodes, + bus1="EU kerosene for aviation", + bus2="co2 atmosphere", + carrier="kerosene for aviation", + p_nom_extendable=True, + efficiency2=costs.at["oil", "CO2 intensity"], ) # TODO simplify bus expression @@ -3018,28 +3084,36 @@ def add_agriculture(n, costs): ) if oil_share > 0: - n.madd( + n.add( + "Bus", + "EU agriculture machinery oil", + location="EU", + carrier="agriculture machinery oil", + unit="MWh_LHV", + ) + + n.add( "Load", - ["agriculture machinery oil"], - bus=spatial.oil.nodes, + "agriculture machinery oil", + bus="EU agriculture machinery oil", carrier="agriculture machinery oil", p_set=oil_share * machinery_nodal_energy.sum() * 1e6 / nhours, ) - co2 = ( - oil_share - * machinery_nodal_energy.sum() - * 1e6 - / nhours - * costs.at["oil", "CO2 intensity"] - ) + if len(spatial.oil.nodes) == 1: + link_names = ["EU agriculture machinery oil"] + else: + link_names = nodes + " agriculture machinery oil" - n.add( - "Load", - "agriculture machinery oil emissions", - bus="co2 atmosphere", - carrier="agriculture machinery oil emissions", - p_set=-co2, + n.madd( + "Link", + link_names, + bus0=spatial.oil.nodes, + bus1="EU agriculture machinery oil", + bus2="co2 atmosphere", + carrier="agriculture machinery oil", + p_nom_extendable=True, + efficiency2=costs.at["oil", "CO2 intensity"], ) From 94afba7c5d195b2cf6d7d17016e78040c4659440 Mon Sep 17 00:00:00 2001 From: chrstphtrs Date: Tue, 24 Oct 2023 16:39:33 +0200 Subject: [PATCH 03/23] add coal tech_color to config --- config/config.default.yaml | 1 + 1 file changed, 1 insertion(+) diff --git a/config/config.default.yaml b/config/config.default.yaml index 325bbbaa2..cafb9d1d1 100644 --- a/config/config.default.yaml +++ b/config/config.default.yaml @@ -816,6 +816,7 @@ plotting: Coal: '#545454' coal: '#545454' Coal marginal: '#545454' + coal for industry: '#343434' solid: '#545454' Lignite: '#826837' lignite: '#826837' From 7cb677d0e6056f52560e6ddb53e88f76861d8fc2 Mon Sep 17 00:00:00 2001 From: chrstphtrs Date: Tue, 24 Oct 2023 16:39:58 +0200 Subject: [PATCH 04/23] clean up function add_co2limit_country --- scripts/solve_network.py | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/scripts/solve_network.py b/scripts/solve_network.py index f5dd79e0b..b372b366e 100644 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -764,7 +764,7 @@ def add_pipe_retrofit_constraint(n): n.model.add_constraints(lhs == rhs, name="Link-pipe_retrofit") -def add_co2limit_country(n, config, limit_countries, nyears=1.0): +def add_co2limit_country(n, limit_countries, nyears=1.0): """ Add a set of emissions limit constraints for specified countries. @@ -780,10 +780,7 @@ def add_co2limit_country(n, config, limit_countries, nyears=1.0): """ logger.info(f"Adding CO2 budget limit for each country as per unit of 1990 levels") - # TODO: n.config (submodule) vs snakemake.config (main module, overwrite/overwritten config)? - # countries = config.countries - # print(config) - countries = ['AT', 'BE', 'CH', 'CZ', 'DE', 'DK', 'FR', 'GB', 'LU', 'NL', 'NO', 'PL', 'SE'] + countries = n.config["countries"] # TODO: import function from prepare_sector_network? Move to common place? sectors = emission_sectors_from_opts(opts) @@ -814,7 +811,7 @@ def add_co2limit_country(n, config, limit_countries, nyears=1.0): n.links["efficiency"].apply(lambda x: 1.0).rename("efficiency0") ) elif port == str(1): - efficiency = n.links["efficiency"].rename("efficiency1") + efficiency = n.links["efficiency"] else: efficiency = n.links[f"efficiency{port}"] mask = n.links[f"bus{port}"].map(n.buses.carrier).eq("co2") @@ -889,7 +886,7 @@ def extra_functionality(n, snapshots): # add co2 constraint for each country logger.info(f"Add CO2 limit for each country") - add_co2limit_country(n, config, limit_countries, nyears) + add_co2limit_country(n, limit_countries, nyears) def solve_network(n, config, solving, opts="", **kwargs): From e2b2eafbc12e17254a0a517c87ae76bf08962585 Mon Sep 17 00:00:00 2001 From: chrstphtrs Date: Tue, 24 Oct 2023 16:46:58 +0200 Subject: [PATCH 05/23] add geographical resolution to oil and methanol for options['co2_budget_national'] to include all necessary links in national co2 budget constraints --- scripts/add_existing_baseyear.py | 2 +- scripts/prepare_sector_network.py | 187 +++++++++++++++++------------- 2 files changed, 106 insertions(+), 83 deletions(-) diff --git a/scripts/add_existing_baseyear.py b/scripts/add_existing_baseyear.py index 1474b0042..7ddc6b1da 100644 --- a/scripts/add_existing_baseyear.py +++ b/scripts/add_existing_baseyear.py @@ -303,7 +303,7 @@ def add_power_capacities_installed_before_baseyear(n, grouping_years, costs, bas else: bus0 = vars(spatial)[carrier[generator]].nodes if "EU" not in vars(spatial)[carrier[generator]].locations: - bus0 = bus0.intersection(capacity.index + " gas") + bus0 = bus0.intersection(capacity.index + " " + carrier[generator]) # check for missing bus missing_bus = pd.Index(bus0).difference(n.buses.index) diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py index 989bdb780..34bfdce72 100644 --- a/scripts/prepare_sector_network.py +++ b/scripts/prepare_sector_network.py @@ -128,13 +128,33 @@ def define_spatial(nodes, options): # methanol spatial.methanol = SimpleNamespace() - spatial.methanol.nodes = ["EU methanol"] - spatial.methanol.locations = ["EU"] + + if options["co2_budget_national"]: + spatial.methanol.nodes = nodes + " methanol" + spatial.methanol.locations = nodes + spatial.methanol.shipping = nodes + " shipping methanol" + else: + spatial.methanol.nodes = ["EU methanol"] + spatial.methanol.locations = ["EU"] + spatial.methanol.shipping = ["EU shipping methanol"] # oil spatial.oil = SimpleNamespace() - spatial.oil.nodes = ["EU oil"] - spatial.oil.locations = ["EU"] + + if options["co2_budget_national"]: + spatial.oil.nodes = nodes + " oil" + spatial.oil.locations = nodes + spatial.oil.naphtha = nodes + " naphtha for industry" + spatial.oil.kerosene = nodes + " kerosene for aviation" + spatial.oil.shipping = nodes + " shipping oil" + spatial.oil.agriculture_machinery = nodes + " agriculture machinery oil" + else: + spatial.oil.nodes = ["EU oil"] + spatial.oil.locations = ["EU"] + spatial.oil.naphtha = ["EU naphtha for industry"] + spatial.oil.kerosene = ["EU kerosene for aviation"] + spatial.oil.shipping = ["EU shipping oil"] + spatial.oil.agriculture_machinery = ["EU agriculture machinery oil"] spatial.oil.land_transport = nodes + " land transport oil" # uranium @@ -2613,34 +2633,34 @@ def add_industry(n, costs): efficiency = ( options["shipping_oil_efficiency"] / options["shipping_methanol_efficiency"] ) - p_set_methanol = shipping_methanol_share * p_set.sum() * efficiency - n.add( + # need to aggregate potentials if methanol not nodally resolved + if options["co2_budget_national"]: + p_set_methanol = shipping_methanol_share * p_set * efficiency + else: + p_set_methanol = shipping_methanol_share * p_set.sum() * efficiency + + n.madd( "Bus", - "EU shipping methanol", - location="EU", + spatial.methanol.shipping, + location=spatial.methanol.locations, carrier="shipping methanol", unit="MWh_LHV", ) - n.add( + n.madd( "Load", - "shipping methanol", - bus="EU shipping methanol", + spatial.methanol.shipping, + bus=spatial.methanol.shipping, carrier="shipping methanol", p_set=p_set_methanol, ) - if len(spatial.methanol.nodes) == 1: - link_names = ["EU shipping methanol"] - else: - link_names = nodes + " shipping methanol" - n.madd( "Link", - link_names, + spatial.methanol.shipping, bus0=spatial.methanol.nodes, - bus1="EU shipping methanol", + bus1=spatial.methanol.shipping, bus2="co2 atmosphere", carrier="shipping methanol", p_nom_extendable=True, @@ -2678,34 +2698,33 @@ def add_industry(n, costs): ) if shipping_oil_share: - p_set_oil = shipping_oil_share * p_set.sum() + # need to aggregate potentials if oil not nodally resolved + if options["co2_budget_national"]: + p_set_oil = shipping_oil_share * p_set + else: + p_set_oil = shipping_oil_share * p_set.sum() - n.add( + n.madd( "Bus", - "EU shipping oil", - location="EU", + spatial.oil.shipping, + location=spatial.oil.locations, carrier="shipping oil", unit="MWh_LHV", ) - n.add( + n.madd( "Load", - "shipping oil", - bus="EU shipping oil", + spatial.oil.shipping, + bus=spatial.oil.shipping, carrier="shipping oil", p_set=p_set_oil, ) - if len(spatial.oil.nodes) == 1: - link_names = ["EU shipping oil"] - else: - link_names = nodes + " shipping oil" - n.madd( "Link", - link_names, + spatial.oil.shipping, bus0=spatial.oil.nodes, - bus1="EU shipping oil", + bus1=spatial.oil.shipping, bus2="co2 atmosphere", carrier="shipping oil", p_nom_extendable=True, @@ -2761,34 +2780,33 @@ def add_industry(n, costs): # NB: CO2 gets released again to atmosphere when plastics decay # except for the process emissions when naphtha is used for petrochemicals, which can be captured with other industry process emissions # convert process emissions from feedstock from MtCO2 to energy demand - p_set = demand_factor * (industrial_demand.loc[nodes, "naphtha"] - industrial_demand.loc[nodes, "process emission from feedstock"] / costs.at["oil", "CO2 intensity"]).sum() / nhours + # need to aggregate potentials if oil not nodally resolved + if options["co2_budget_national"]: + p_set = demand_factor * (industrial_demand.loc[nodes, "naphtha"] - industrial_demand.loc[nodes, "process emission from feedstock"] / costs.at["oil", "CO2 intensity"]) / nhours + else: + p_set = demand_factor * (industrial_demand.loc[nodes, "naphtha"] - industrial_demand.loc[nodes, "process emission from feedstock"] / costs.at["oil", "CO2 intensity"]).sum() / nhours - n.add( + n.madd( "Bus", - "EU naphtha for industry", - location="EU", + spatial.oil.naphtha, + location=spatial.oil.locations, carrier="naphtha for industry", unit="MWh_LHV", ) - n.add( + n.madd( "Load", - "naphtha for industry", - bus="EU naphtha for industry", + spatial.oil.naphtha, + bus=spatial.oil.naphtha, carrier="naphtha for industry", p_set=p_set, ) - if len(spatial.oil.nodes) == 1: - link_names = ["EU naphtha for industry"] - else: - link_names = nodes + " naphtha for industry" - n.madd( "Link", - link_names, + spatial.oil.naphtha, bus0=spatial.oil.nodes, - bus1="EU naphtha for industry", + bus1=spatial.oil.naphtha, bus2="co2 atmosphere", carrier="naphtha for industry", p_nom_extendable=True, @@ -2797,42 +2815,47 @@ def add_industry(n, costs): # aviation demand_factor = options.get("aviation_demand_factor", 1) - all_aviation = ["total international aviation", "total domestic aviation"] - p_set = ( - demand_factor - * pop_weighted_energy_totals.loc[nodes, all_aviation].sum(axis=1).sum() - * 1e6 - / nhours - ) if demand_factor != 1: logger.warning(f"Changing aviation demand by {demand_factor*100-100:+.2f}%.") - n.add( + all_aviation = ["total international aviation", "total domestic aviation"] + # need to aggregate potentials if oil not nodally resolved + if options["co2_budget_national"]: + p_set = ( + demand_factor + * pop_weighted_energy_totals.loc[nodes, all_aviation].sum(axis=1) + * 1e6 + / nhours + ) + else: + p_set = ( + demand_factor + * pop_weighted_energy_totals.loc[nodes, all_aviation].sum(axis=1).sum() + * 1e6 + / nhours + ) + + n.madd( "Bus", - "EU kerosene for aviation", - location="EU", + spatial.oil.kerosene, + location=spatial.oil.locations, carrier="kerosene for aviation", unit="MWh_LHV", ) - n.add( + n.madd( "Load", - "kerosene for aviation", - bus="EU kerosene for aviation", + spatial.oil.kerosene, + bus=spatial.oil.kerosene, carrier="kerosene for aviation", p_set=p_set, ) - if len(spatial.oil.nodes) == 1: - link_names = ["EU kerosene for aviation"] - else: - link_names = nodes + " kerosene for aviation" - n.madd( "Link", - link_names, + spatial.oil.kerosene, bus0=spatial.oil.nodes, - bus1="EU kerosene for aviation", + bus1=spatial.oil.kerosene, bus2="co2 atmosphere", carrier="kerosene for aviation", p_nom_extendable=True, @@ -3062,7 +3085,7 @@ def add_agriculture(n, costs): machinery_nodal_energy = pop_weighted_energy_totals.loc[ nodes, "total agriculture machinery" - ] + ] * 1e6 if electric_share > 0: efficiency_gain = ( @@ -3079,37 +3102,37 @@ def add_agriculture(n, costs): p_set=electric_share / efficiency_gain * machinery_nodal_energy - * 1e6 / nhours, ) if oil_share > 0: - n.add( + # need to aggregate potentials if oil not nodally resolved + if options["co2_budget_national"]: + p_set = oil_share * machinery_nodal_energy / nhours + else: + p_set = oil_share * machinery_nodal_energy.sum() / nhours + + n.madd( "Bus", - "EU agriculture machinery oil", - location="EU", + spatial.oil.agriculture_machinery, + location=spatial.oil.locations, carrier="agriculture machinery oil", unit="MWh_LHV", ) - n.add( + n.madd( "Load", - "agriculture machinery oil", - bus="EU agriculture machinery oil", + spatial.oil.agriculture_machinery, + bus=spatial.oil.agriculture_machinery, carrier="agriculture machinery oil", - p_set=oil_share * machinery_nodal_energy.sum() * 1e6 / nhours, + p_set=p_set, ) - if len(spatial.oil.nodes) == 1: - link_names = ["EU agriculture machinery oil"] - else: - link_names = nodes + " agriculture machinery oil" - n.madd( "Link", - link_names, + spatial.oil.agriculture_machinery, bus0=spatial.oil.nodes, - bus1="EU agriculture machinery oil", + bus1=spatial.oil.agriculture_machinery, bus2="co2 atmosphere", carrier="agriculture machinery oil", p_nom_extendable=True, From 2ad9ca8f7b10155d2b8f738d11a4948ac5f17fb1 Mon Sep 17 00:00:00 2001 From: chrstphtrs Date: Thu, 26 Oct 2023 11:17:57 +0200 Subject: [PATCH 06/23] add regionalised oil load for process emissions from naphtha as feedstock --- scripts/prepare_sector_network.py | 30 +++++++++++++++++++++++++++--- 1 file changed, 27 insertions(+), 3 deletions(-) diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py index 34bfdce72..548301063 100644 --- a/scripts/prepare_sector_network.py +++ b/scripts/prepare_sector_network.py @@ -2782,9 +2782,25 @@ def add_industry(n, costs): # convert process emissions from feedstock from MtCO2 to energy demand # need to aggregate potentials if oil not nodally resolved if options["co2_budget_national"]: - p_set = demand_factor * (industrial_demand.loc[nodes, "naphtha"] - industrial_demand.loc[nodes, "process emission from feedstock"] / costs.at["oil", "CO2 intensity"]) / nhours + p_set_plastics = demand_factor * (industrial_demand.loc[nodes, "naphtha"] - industrial_demand.loc[nodes, "process emission from feedstock"] / costs.at["oil", "CO2 intensity"]) / nhours else: - p_set = demand_factor * (industrial_demand.loc[nodes, "naphtha"] - industrial_demand.loc[nodes, "process emission from feedstock"] / costs.at["oil", "CO2 intensity"]).sum() / nhours + p_set_plastics = demand_factor * (industrial_demand.loc[nodes, "naphtha"] - industrial_demand.loc[nodes, "process emission from feedstock"] / costs.at["oil", "CO2 intensity"]).sum() / nhours + + if options["co2_budget_national"]: + p_set_process_emissions = ( + demand_factor + * (industrial_demand.loc[nodes, "process emission from feedstock"] + / costs.at["oil", "CO2 intensity"]) + / nhours + ) + else: + p_set_process_emissions = ( + demand_factor + * (industrial_demand.loc[nodes, "process emission from feedstock"] + / costs.at["oil", "CO2 intensity"] + ).sum() + / nhours + ) n.madd( "Bus", @@ -2799,7 +2815,15 @@ def add_industry(n, costs): spatial.oil.naphtha, bus=spatial.oil.naphtha, carrier="naphtha for industry", - p_set=p_set, + p_set=p_set_plastics, + ) + + n.madd( + "Load", + ["naphtha for industry into process emissions from feedstock"], + bus=spatial.oil.nodes, + carrier="naphtha for industry", + p_set=p_set_process_emissions, ) n.madd( From 82ac430fd92f2724918ff0e25568fdd57b75a9a5 Mon Sep 17 00:00:00 2001 From: chrstphtrs Date: Wed, 8 Nov 2023 09:57:24 +0100 Subject: [PATCH 07/23] fix spatial resolution for solid biomass links and naphtha oil loads under 'co2_spatial: true' flag --- scripts/prepare_sector_network.py | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py index 548301063..a5ca8941f 100644 --- a/scripts/prepare_sector_network.py +++ b/scripts/prepare_sector_network.py @@ -145,6 +145,7 @@ def define_spatial(nodes, options): spatial.oil.nodes = nodes + " oil" spatial.oil.locations = nodes spatial.oil.naphtha = nodes + " naphtha for industry" + spatial.oil.naphtha_process_emissions = nodes + " naphtha process emissions" spatial.oil.kerosene = nodes + " kerosene for aviation" spatial.oil.shipping = nodes + " shipping oil" spatial.oil.agriculture_machinery = nodes + " agriculture machinery oil" @@ -152,6 +153,7 @@ def define_spatial(nodes, options): spatial.oil.nodes = ["EU oil"] spatial.oil.locations = ["EU"] spatial.oil.naphtha = ["EU naphtha for industry"] + spatial.oil.naphtha_process_emissions = "EU naphtha process emissions" spatial.oil.kerosene = ["EU kerosene for aviation"] spatial.oil.shipping = ["EU shipping oil"] spatial.oil.agriculture_machinery = ["EU agriculture machinery oil"] @@ -2443,9 +2445,14 @@ def add_industry(n, costs): efficiency=1.0, ) + if len(spatial.biomass.industry_cc)<=1 and len(spatial.co2.nodes)>1: + link_names = nodes + " " + spatial.biomass.industry_cc + else: + link_names = spatial.biomass.industry_cc + n.madd( "Link", - spatial.biomass.industry_cc, + link_names, bus0=spatial.biomass.nodes, bus1=spatial.biomass.industry, bus2="co2 atmosphere", @@ -2820,7 +2827,7 @@ def add_industry(n, costs): n.madd( "Load", - ["naphtha for industry into process emissions from feedstock"], + spatial.oil.naphtha_process_emissions, bus=spatial.oil.nodes, carrier="naphtha for industry", p_set=p_set_process_emissions, From d9ec127f996f854cc775cdfcc6db18cd26cf3ea5 Mon Sep 17 00:00:00 2001 From: chrstphtrs Date: Tue, 21 Nov 2023 14:55:32 +0100 Subject: [PATCH 08/23] Add process emissions to country emissions constraint, fix snapshot weighting --- scripts/solve_network.py | 19 +++++++++++++++++-- 1 file changed, 17 insertions(+), 2 deletions(-) diff --git a/scripts/solve_network.py b/scripts/solve_network.py index b372b366e..e2edb2eb4 100644 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -795,7 +795,7 @@ def add_co2limit_country(n, limit_countries, nyears=1.0): p = n.model["Link-p"] # dimension: (time, component) - # NB: Most country-specific links retain their locational information in bus1 (except for DAC, where it is in bus2) + # NB: Most country-specific links retain their locational information in bus1 (except for DAC, where it is in bus2, and process emissions, where it is in bus0) country = n.links.bus1.map(n.buses.location).map(n.buses.country) country_DAC = ( n.links[n.links.carrier == "DAC"] @@ -803,6 +803,12 @@ def add_co2limit_country(n, limit_countries, nyears=1.0): .map(n.buses.country) ) country[country_DAC.index] = country_DAC + country_process_emissions = ( + n.links[n.links.carrier.str.contains("process emissions")] + .bus0.map(n.buses.location) + .map(n.buses.country) + ) + country[country_process_emissions.index] = country_process_emissions lhs = [] for port in [col[3:] for col in n.links if col.startswith("bus")]: @@ -818,13 +824,18 @@ def add_co2limit_country(n, limit_countries, nyears=1.0): idx = n.links[mask].index + international = n.links.carrier.map( + lambda x: 0.4 if x in ["kerosene for aviation", "shipping oil"] else 1.0 + ) grouping = country.loc[idx] if not grouping.isnull().all(): expr = ( - (p.loc[:, idx] * efficiency[idx]) + ((p.loc[:, idx] * efficiency[idx] * international[idx]) .groupby(grouping, axis=1) .sum() + *n.snapshot_weightings.generators + ) .sum(dims="snapshot") ) lhs.append(expr) @@ -935,6 +946,10 @@ def solve_network(n, config, solving, opts="", **kwargs): f"Solving status '{status}' with termination condition '{condition}'" ) if "infeasible" in condition: + m = n.model + labels = m.compute_infeasibilities() + print(labels) + m.print_infeasibilities() raise RuntimeError("Solving status 'infeasible'") return n From e8324b9c2788339837caf8718898b14f879ea4d8 Mon Sep 17 00:00:00 2001 From: lisazeyen <35347358+lisazeyen@users.noreply.github.com> Date: Fri, 24 Nov 2023 09:58:24 +0100 Subject: [PATCH 09/23] fix bug when oil copper plated --- scripts/prepare_sector_network.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py index a5ca8941f..81e4d6e30 100644 --- a/scripts/prepare_sector_network.py +++ b/scripts/prepare_sector_network.py @@ -153,7 +153,7 @@ def define_spatial(nodes, options): spatial.oil.nodes = ["EU oil"] spatial.oil.locations = ["EU"] spatial.oil.naphtha = ["EU naphtha for industry"] - spatial.oil.naphtha_process_emissions = "EU naphtha process emissions" + spatial.oil.naphtha_process_emissions = ["EU naphtha process emissions"] spatial.oil.kerosene = ["EU kerosene for aviation"] spatial.oil.shipping = ["EU shipping oil"] spatial.oil.agriculture_machinery = ["EU agriculture machinery oil"] From 3ff925e797574afc11193d7f63316fdbdde03e12 Mon Sep 17 00:00:00 2001 From: lisazeyen <35347358+lisazeyen@users.noreply.github.com> Date: Fri, 24 Nov 2023 10:00:07 +0100 Subject: [PATCH 10/23] add load shedding for all energy carriers --- scripts/solve_network.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/solve_network.py b/scripts/solve_network.py index e2edb2eb4..97c78dad5 100644 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -352,7 +352,7 @@ def prepare_network( # http://journal.frontiersin.org/article/10.3389/fenrg.2015.00055/full # TODO: retrieve color and nice name from config n.add("Carrier", "load", color="#dd2e23", nice_name="Load shedding") - buses_i = n.buses.query("carrier == 'AC'").index + buses_i = n.buses.index if not np.isscalar(load_shedding): # TODO: do not scale via sign attribute (use Eur/MWh instead of Eur/kWh) load_shedding = 1e2 # Eur/kWh From cea62de438b7c358bf23ad306dccb300b13beab7 Mon Sep 17 00:00:00 2001 From: Tom Brown Date: Mon, 4 Dec 2023 16:46:11 +0100 Subject: [PATCH 11/23] solve_network: quick fix so duals can be read from CO2 constrain --- scripts/solve_network.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/scripts/solve_network.py b/scripts/solve_network.py index 97c78dad5..2413f4c99 100644 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -841,10 +841,10 @@ def add_co2limit_country(n, limit_countries, nyears=1.0): lhs.append(expr) lhs = sum(lhs) # dimension: (country) - lhs = lhs.rename({list(lhs.dims.keys())[0]: "country"}) + lhs = lhs.rename({list(lhs.dims.keys())[0]: "snapshot"}) rhs = pd.Series(co2_limit_countries) # dimension: (country) - for ct in lhs.indexes["country"]: + for ct in lhs.indexes["snapshot"]: n.model.add_constraints( lhs.loc[ct] <= rhs[ct], name=f"GlobalConstraint-co2_limit_per_country{ct}", From 66178a5a27625b7055d029403c66bd7ac6df1da5 Mon Sep 17 00:00:00 2001 From: Tom Brown Date: Mon, 4 Dec 2023 16:46:45 +0100 Subject: [PATCH 12/23] solve_network: fix sign for country CO2 when bus0=atmosphere So that DAC extracts CO2 rather than pumping into air; for p>0, link withdraws from bus0, but injects into bus1/2/3, so you have to take account of this sign difference- --- scripts/solve_network.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/solve_network.py b/scripts/solve_network.py index 2413f4c99..53170da9a 100644 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -814,7 +814,7 @@ def add_co2limit_country(n, limit_countries, nyears=1.0): for port in [col[3:] for col in n.links if col.startswith("bus")]: if port == str(0): efficiency = ( - n.links["efficiency"].apply(lambda x: 1.0).rename("efficiency0") + n.links["efficiency"].apply(lambda x: -1.0).rename("efficiency0") ) elif port == str(1): efficiency = n.links["efficiency"] From bbf9ca2d9be0af6fe80ffcc667556405ab0bddbc Mon Sep 17 00:00:00 2001 From: Tom Brown Date: Fri, 8 Dec 2023 11:58:28 +0100 Subject: [PATCH 13/23] bug fix: naming of p_set when co2_national is True Without this naming fix, the p_set is a NaN once added --- scripts/prepare_sector_network.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py index 81e4d6e30..606e17b31 100644 --- a/scripts/prepare_sector_network.py +++ b/scripts/prepare_sector_network.py @@ -2643,7 +2643,7 @@ def add_industry(n, costs): # need to aggregate potentials if methanol not nodally resolved if options["co2_budget_national"]: - p_set_methanol = shipping_methanol_share * p_set * efficiency + p_set_methanol = shipping_methanol_share * p_set.rename(lambda x : x + " shipping methanol") * efficiency else: p_set_methanol = shipping_methanol_share * p_set.sum() * efficiency @@ -2707,7 +2707,7 @@ def add_industry(n, costs): if shipping_oil_share: # need to aggregate potentials if oil not nodally resolved if options["co2_budget_national"]: - p_set_oil = shipping_oil_share * p_set + p_set_oil = shipping_oil_share * p_set.rename(lambda x: x + " shipping oil") else: p_set_oil = shipping_oil_share * p_set.sum() @@ -2789,15 +2789,16 @@ def add_industry(n, costs): # convert process emissions from feedstock from MtCO2 to energy demand # need to aggregate potentials if oil not nodally resolved if options["co2_budget_national"]: - p_set_plastics = demand_factor * (industrial_demand.loc[nodes, "naphtha"] - industrial_demand.loc[nodes, "process emission from feedstock"] / costs.at["oil", "CO2 intensity"]) / nhours + p_set_plastics = demand_factor * (industrial_demand.loc[nodes, "naphtha"] - industrial_demand.loc[nodes, "process emission from feedstock"] / costs.at["oil", "CO2 intensity"]).rename(lambda x: x + " naphtha for industry") / nhours else: p_set_plastics = demand_factor * (industrial_demand.loc[nodes, "naphtha"] - industrial_demand.loc[nodes, "process emission from feedstock"] / costs.at["oil", "CO2 intensity"]).sum() / nhours + if options["co2_budget_national"]: p_set_process_emissions = ( demand_factor * (industrial_demand.loc[nodes, "process emission from feedstock"] - / costs.at["oil", "CO2 intensity"]) + / costs.at["oil", "CO2 intensity"]).rename(lambda x: x + " naphtha process emissions") / nhours ) else: @@ -2857,7 +2858,7 @@ def add_industry(n, costs): * pop_weighted_energy_totals.loc[nodes, all_aviation].sum(axis=1) * 1e6 / nhours - ) + ).rename(lambda x: x + " kerosene for aviation") else: p_set = ( demand_factor @@ -3139,7 +3140,7 @@ def add_agriculture(n, costs): if oil_share > 0: # need to aggregate potentials if oil not nodally resolved if options["co2_budget_national"]: - p_set = oil_share * machinery_nodal_energy / nhours + p_set = oil_share * machinery_nodal_energy.rename(lambda x: x + " agriculture machinery oil") / nhours else: p_set = oil_share * machinery_nodal_energy.sum() / nhours From 2d323d1b879751bc96303bd6c6a54fda2c90eccb Mon Sep 17 00:00:00 2001 From: Tom Brown Date: Fri, 8 Dec 2023 12:27:07 +0100 Subject: [PATCH 14/23] bug fix: ICE efficiency for land transport was applied twice This was overestimating ICE oil demand by factor 1/0.3. --- scripts/prepare_sector_network.py | 1 - 1 file changed, 1 deletion(-) diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py index 606e17b31..342a6b151 100644 --- a/scripts/prepare_sector_network.py +++ b/scripts/prepare_sector_network.py @@ -1613,7 +1613,6 @@ def add_land_transport(n, costs): bus1=spatial.oil.land_transport, bus2="co2 atmosphere", carrier="land transport oil", - efficiency=ice_efficiency, efficiency2=costs.at["oil", "CO2 intensity"], p_nom_extendable=True, ) From 00e86e6435816fe007fb25d62de90cf58fbc01c4 Mon Sep 17 00:00:00 2001 From: Tom Brown Date: Fri, 8 Dec 2023 13:28:08 +0100 Subject: [PATCH 15/23] bug fix: route process emissions from steam cracker to correct bus Now naphtha demand causes process emissions from steak crackers to route to process emissions bus, then rest of CO2 goes to atmosphere. --- scripts/prepare_sector_network.py | 46 +++++++------------------------ 1 file changed, 10 insertions(+), 36 deletions(-) diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py index 342a6b151..8e995dd65 100644 --- a/scripts/prepare_sector_network.py +++ b/scripts/prepare_sector_network.py @@ -145,7 +145,6 @@ def define_spatial(nodes, options): spatial.oil.nodes = nodes + " oil" spatial.oil.locations = nodes spatial.oil.naphtha = nodes + " naphtha for industry" - spatial.oil.naphtha_process_emissions = nodes + " naphtha process emissions" spatial.oil.kerosene = nodes + " kerosene for aviation" spatial.oil.shipping = nodes + " shipping oil" spatial.oil.agriculture_machinery = nodes + " agriculture machinery oil" @@ -153,7 +152,6 @@ def define_spatial(nodes, options): spatial.oil.nodes = ["EU oil"] spatial.oil.locations = ["EU"] spatial.oil.naphtha = ["EU naphtha for industry"] - spatial.oil.naphtha_process_emissions = ["EU naphtha process emissions"] spatial.oil.kerosene = ["EU kerosene for aviation"] spatial.oil.shipping = ["EU shipping oil"] spatial.oil.agriculture_machinery = ["EU agriculture machinery oil"] @@ -2783,31 +2781,10 @@ def add_industry(n, costs): if demand_factor != 1: logger.warning(f"Changing HVC demand by {demand_factor*100-100:+.2f}%.") - # NB: CO2 gets released again to atmosphere when plastics decay - # except for the process emissions when naphtha is used for petrochemicals, which can be captured with other industry process emissions - # convert process emissions from feedstock from MtCO2 to energy demand - # need to aggregate potentials if oil not nodally resolved - if options["co2_budget_national"]: - p_set_plastics = demand_factor * (industrial_demand.loc[nodes, "naphtha"] - industrial_demand.loc[nodes, "process emission from feedstock"] / costs.at["oil", "CO2 intensity"]).rename(lambda x: x + " naphtha for industry") / nhours - else: - p_set_plastics = demand_factor * (industrial_demand.loc[nodes, "naphtha"] - industrial_demand.loc[nodes, "process emission from feedstock"] / costs.at["oil", "CO2 intensity"]).sum() / nhours - - if options["co2_budget_national"]: - p_set_process_emissions = ( - demand_factor - * (industrial_demand.loc[nodes, "process emission from feedstock"] - / costs.at["oil", "CO2 intensity"]).rename(lambda x: x + " naphtha process emissions") - / nhours - ) + p_set_plastics = demand_factor * industrial_demand.loc[nodes, "naphtha"].rename(lambda x: x + " naphtha for industry") / nhours else: - p_set_process_emissions = ( - demand_factor - * (industrial_demand.loc[nodes, "process emission from feedstock"] - / costs.at["oil", "CO2 intensity"] - ).sum() - / nhours - ) + p_set_plastics = demand_factor * industrial_demand.loc[nodes, "naphtha"].sum() / nhours n.madd( "Bus", @@ -2825,13 +2802,10 @@ def add_industry(n, costs): p_set=p_set_plastics, ) - n.madd( - "Load", - spatial.oil.naphtha_process_emissions, - bus=spatial.oil.nodes, - carrier="naphtha for industry", - p_set=p_set_process_emissions, - ) + # some CO2 from naphtha are process emissions from steam cracker + # rest of CO2 released to atmosphere either in waste-to-energy or decay + process_co2_per_naphtha = industrial_demand.loc[nodes, "process emission from feedstock"].sum() / industrial_demand.loc[nodes, "naphtha"].sum() + emitted_co2_per_naphtha = costs.at["oil", "CO2 intensity"] - process_co2_per_naphtha n.madd( "Link", @@ -2839,9 +2813,11 @@ def add_industry(n, costs): bus0=spatial.oil.nodes, bus1=spatial.oil.naphtha, bus2="co2 atmosphere", + bus3=spatial.co2.process_emissions, carrier="naphtha for industry", p_nom_extendable=True, - efficiency2=costs.at["oil", "CO2 intensity"], + efficiency2=emitted_co2_per_naphtha, + efficiency3=process_co2_per_naphtha, ) # aviation @@ -2941,7 +2917,7 @@ def add_industry(n, costs): unit="t_co2", ) - sel = ["process emission", "process emission from feedstock"] + sel = ["process emission"] if options["co2_spatial"] or options["co2network"]: p_set = ( -industrial_demand.loc[nodes, sel] @@ -2952,8 +2928,6 @@ def add_industry(n, costs): else: p_set = -industrial_demand.loc[nodes, sel].sum(axis=1).sum() / nhours - # this should be process emissions fossil+feedstock - # then need load on atmosphere for feedstock emissions that are currently going to atmosphere via Link Fischer-Tropsch demand n.madd( "Load", spatial.co2.process_emissions, From 326ed63329d55d5a84f9230840161cdb3673e27a Mon Sep 17 00:00:00 2001 From: Tom Brown Date: Fri, 8 Dec 2023 17:53:28 +0100 Subject: [PATCH 16/23] add_brownfield: disable grid expansion if LV already hit Numerical problems were causing infeasibilities otherwise --- scripts/add_brownfield.py | 28 ++++++++++++++++++++++++++++ 1 file changed, 28 insertions(+) diff --git a/scripts/add_brownfield.py b/scripts/add_brownfield.py index 741025801..fb1453fd7 100644 --- a/scripts/add_brownfield.py +++ b/scripts/add_brownfield.py @@ -119,6 +119,32 @@ def add_brownfield(n, n_p, year): n.links.loc[new_pipes, "p_nom"] = 0.0 n.links.loc[new_pipes, "p_nom_min"] = 0.0 +def disable_grid_expansion_if_LV_limit_hit(n): + if not "lv_limit" in n.global_constraints.index: + return + + #calculate minimum LV + attr = "nom_min" + dc = n.links.index[n.links.carrier == "DC"] + tot = (n.lines["s_" + attr]*n.lines["length"]).sum() + (n.links.loc[dc,"p_" + attr]*n.links.loc[dc,"length"]).sum() + + diff = n.global_constraints.at["lv_limit","constant"]-tot + + #allow small numerical differences + limit = 1 + + if diff < limit: + logger.info(f"LV is already reached (gap {diff}), disabling expansion and LV limit") + expandable_acs = n.lines.index[n.lines.s_nom_extendable] + n.lines.loc[expandable_acs,"s_nom_extendable"] = False + n.lines.loc[expandable_acs,"s_nom"] = n.lines.loc[expandable_acs,"s_nom_min"] + + expandable_dcs = n.links.index[n.links.p_nom_extendable & (n.links.carrier == "DC")] + n.links.loc[expandable_dcs,"p_nom_extendable"] = False + n.links.loc[expandable_dcs,"p_nom"] = n.links.loc[expandable_dcs,"p_nom_min"] + + n.global_constraints.drop("lv_limit", + inplace=True) if __name__ == "__main__": if "snakemake" not in globals(): @@ -150,5 +176,7 @@ def add_brownfield(n, n_p, year): add_brownfield(n, n_p, year) + disable_grid_expansion_if_LV_limit_hit(n) + n.meta = dict(snakemake.config, **dict(wildcards=dict(snakemake.wildcards))) n.export_to_netcdf(snakemake.output[0]) From 830019a6e5d5ced3403bd4d5d9e28e2d66fe621b Mon Sep 17 00:00:00 2001 From: Tom Brown Date: Fri, 15 Dec 2023 09:50:47 +0100 Subject: [PATCH 17/23] add rule that allows cost data to be modified --- rules/build_sector.smk | 4 ++-- rules/retrieve.smk | 14 ++++++++++++++ scripts/modify_cost_data.py | 12 ++++++++++++ 3 files changed, 28 insertions(+), 2 deletions(-) create mode 100644 scripts/modify_cost_data.py diff --git a/rules/build_sector.smk b/rules/build_sector.smk index 5a9e8646d..596c0305d 100644 --- a/rules/build_sector.smk +++ b/rules/build_sector.smk @@ -743,9 +743,9 @@ rule prepare_sector_network: else RESOURCES + "biomass_potentials_s{simpl}_{clusters}_{planning_horizons}.csv", heat_profile="data/heat_load_profile_BDEW.csv", - costs="data/costs_{}.csv".format(config["costs"]["year"]) + costs="data/costs_{}-modified.csv".format(config["costs"]["year"]) if config["foresight"] == "overnight" - else "data/costs_{planning_horizons}.csv", + else "data/costs_{planning_horizons}-modified.csv", profile_offwind_ac=RESOURCES + "profile_offwind-ac.nc", profile_offwind_dc=RESOURCES + "profile_offwind-dc.nc", h2_cavern=RESOURCES + "salt_cavern_potentials_s{simpl}_{clusters}.csv", diff --git a/rules/retrieve.smk b/rules/retrieve.smk index b830be251..18b424ff7 100644 --- a/rules/retrieve.smk +++ b/rules/retrieve.smk @@ -259,3 +259,17 @@ if config["enable"]["retrieve"]: "../envs/environment.yaml" script: "../scripts/retrieve_monthly_fuel_prices.py" + + +rule modify_cost_data: + input: + costs="data/costs_{year}.csv", + output: + "data/costs_{year}-modified.csv" + log: + LOGS + "modify_cost_data_{year}.log", + resources: + mem_mb=1000, + retries: 2 + script: + "../scripts/modify_cost_data.py" diff --git a/scripts/modify_cost_data.py b/scripts/modify_cost_data.py new file mode 100644 index 000000000..3e1f12f47 --- /dev/null +++ b/scripts/modify_cost_data.py @@ -0,0 +1,12 @@ + +import pandas as pd + +costs = pd.read_csv(snakemake.input.costs, index_col=[0, 1]).sort_index() + +if "modifications" in snakemake.input.keys(): + modifications = pd.read_csv(snakemake.input.modifications, index_col=[0, 1]).sort_index() + costs.loc[modifications.index] = modifications + print(modifications) + print( costs.loc[modifications.index]) + +costs.to_csv(snakemake.output[0]) From c5a123b4f443a29e8f0f7446ed45ab48230ed1a9 Mon Sep 17 00:00:00 2001 From: Tom Brown Date: Fri, 15 Dec 2023 14:57:03 +0100 Subject: [PATCH 18/23] allow additional functionality for solving to be added by file To add this, overwrite the rule with a new argument: snakemake.input.additional_functionality --- scripts/solve_network.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/scripts/solve_network.py b/scripts/solve_network.py index 53170da9a..dce63efe0 100644 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -44,6 +44,7 @@ from prepare_sector_network import emission_sectors_from_opts + def add_land_use_constraint(n, planning_horizons, config): if "m" in snakemake.wildcards.clusters: _add_land_use_constraint_m(n, planning_horizons, config) @@ -899,6 +900,13 @@ def extra_functionality(n, snapshots): logger.info(f"Add CO2 limit for each country") add_co2limit_country(n, limit_countries, nyears) + if "additional_functionality" in snakemake.input.keys(): + import importlib, os, sys + sys.path.append(os.path.dirname(snakemake.input.additional_functionality)) + additional_functionality = importlib.import_module(os.path.splitext(os.path.basename(snakemake.input.additional_functionality))[0]) + + additional_functionality.additional_functionality(n, snapshots, snakemake.wildcards.planning_horizons) + def solve_network(n, config, solving, opts="", **kwargs): set_of_options = solving["solver"]["options"] From 1a7f093e037ed177468c163863a7bbd929d322c3 Mon Sep 17 00:00:00 2001 From: Tom Brown Date: Fri, 15 Dec 2023 17:18:36 +0100 Subject: [PATCH 19/23] solve: pass wildcards and config to additional_functionality --- scripts/solve_network.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/solve_network.py b/scripts/solve_network.py index dce63efe0..6f88b904a 100644 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -905,7 +905,7 @@ def extra_functionality(n, snapshots): sys.path.append(os.path.dirname(snakemake.input.additional_functionality)) additional_functionality = importlib.import_module(os.path.splitext(os.path.basename(snakemake.input.additional_functionality))[0]) - additional_functionality.additional_functionality(n, snapshots, snakemake.wildcards.planning_horizons) + additional_functionality.additional_functionality(n, snapshots, snakemake.wildcards, config) def solve_network(n, config, solving, opts="", **kwargs): From b3753d73d75eedaddb8da9f6dcfc4bcf8793831a Mon Sep 17 00:00:00 2001 From: Tom Brown Date: Wed, 20 Dec 2023 09:22:40 +0100 Subject: [PATCH 20/23] undo addition of script to allow cost modifications This undoes commit 830019a6e5d5ced3403bd4d5d9e28e2d66fe621b. Reason: this was introduced for the PyPSA-Ariadne derivative, but can be handled more elegantly within the derivative repository. --- rules/build_sector.smk | 4 ++-- rules/retrieve.smk | 14 -------------- scripts/modify_cost_data.py | 12 ------------ 3 files changed, 2 insertions(+), 28 deletions(-) delete mode 100644 scripts/modify_cost_data.py diff --git a/rules/build_sector.smk b/rules/build_sector.smk index 596c0305d..5a9e8646d 100644 --- a/rules/build_sector.smk +++ b/rules/build_sector.smk @@ -743,9 +743,9 @@ rule prepare_sector_network: else RESOURCES + "biomass_potentials_s{simpl}_{clusters}_{planning_horizons}.csv", heat_profile="data/heat_load_profile_BDEW.csv", - costs="data/costs_{}-modified.csv".format(config["costs"]["year"]) + costs="data/costs_{}.csv".format(config["costs"]["year"]) if config["foresight"] == "overnight" - else "data/costs_{planning_horizons}-modified.csv", + else "data/costs_{planning_horizons}.csv", profile_offwind_ac=RESOURCES + "profile_offwind-ac.nc", profile_offwind_dc=RESOURCES + "profile_offwind-dc.nc", h2_cavern=RESOURCES + "salt_cavern_potentials_s{simpl}_{clusters}.csv", diff --git a/rules/retrieve.smk b/rules/retrieve.smk index 18b424ff7..b830be251 100644 --- a/rules/retrieve.smk +++ b/rules/retrieve.smk @@ -259,17 +259,3 @@ if config["enable"]["retrieve"]: "../envs/environment.yaml" script: "../scripts/retrieve_monthly_fuel_prices.py" - - -rule modify_cost_data: - input: - costs="data/costs_{year}.csv", - output: - "data/costs_{year}-modified.csv" - log: - LOGS + "modify_cost_data_{year}.log", - resources: - mem_mb=1000, - retries: 2 - script: - "../scripts/modify_cost_data.py" diff --git a/scripts/modify_cost_data.py b/scripts/modify_cost_data.py deleted file mode 100644 index 3e1f12f47..000000000 --- a/scripts/modify_cost_data.py +++ /dev/null @@ -1,12 +0,0 @@ - -import pandas as pd - -costs = pd.read_csv(snakemake.input.costs, index_col=[0, 1]).sort_index() - -if "modifications" in snakemake.input.keys(): - modifications = pd.read_csv(snakemake.input.modifications, index_col=[0, 1]).sort_index() - costs.loc[modifications.index] = modifications - print(modifications) - print( costs.loc[modifications.index]) - -costs.to_csv(snakemake.output[0]) From 8a55a55d20215dfe32cc49bae5de3e0f05411f1f Mon Sep 17 00:00:00 2001 From: Tom Brown Date: Thu, 21 Dec 2023 16:08:43 +0100 Subject: [PATCH 21/23] copperplate oil/methanol supply; allow demand to be regional Force a single supply bus for oil/methanol (until we allow oil/methanol transport). Introduce new config switches "regional_oil/methanol_demand" that allow demand to be regionalised. This is important if regional CO2 budgets need to be enforced. --- config/config.default.yaml | 2 + scripts/prepare_sector_network.py | 103 ++++++++++++++++-------------- 2 files changed, 57 insertions(+), 48 deletions(-) diff --git a/config/config.default.yaml b/config/config.default.yaml index cafb9d1d1..c1e7ed0f9 100644 --- a/config/config.default.yaml +++ b/config/config.default.yaml @@ -471,6 +471,8 @@ sector: SMR: true SMR_cc: true co2_budget_national: false + regional_methanol_demand: false #set to true if regional CO2 constraints needed + regional_oil_demand: false #set to true if regional CO2 constraints needed regional_co2_sequestration_potential: enable: false attribute: 'conservative estimate Mt' diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py index 8e995dd65..b5a0c0d5d 100644 --- a/scripts/prepare_sector_network.py +++ b/scripts/prepare_sector_network.py @@ -127,35 +127,42 @@ def define_spatial(nodes, options): spatial.h2.locations = nodes # methanol + + #beware: unlike other carriers, uses locations rather than locations+carriername + #this allows to avoid separation between nodes and locations + spatial.methanol = SimpleNamespace() - if options["co2_budget_national"]: - spatial.methanol.nodes = nodes + " methanol" - spatial.methanol.locations = nodes + spatial.methanol.nodes = ["EU methanol"] + spatial.methanol.locations = ["EU"] + + if options["regional_methanol_demand"]: + spatial.methanol.demand_locations = nodes spatial.methanol.shipping = nodes + " shipping methanol" else: - spatial.methanol.nodes = ["EU methanol"] - spatial.methanol.locations = ["EU"] + spatial.methanol.demand_locations = ["EU"] spatial.methanol.shipping = ["EU shipping methanol"] # oil spatial.oil = SimpleNamespace() - if options["co2_budget_national"]: - spatial.oil.nodes = nodes + " oil" - spatial.oil.locations = nodes + spatial.oil.nodes = ["EU oil"] + spatial.oil.locations = ["EU"] + + if options["regional_oil_demand"]: + spatial.oil.demand_locations = nodes spatial.oil.naphtha = nodes + " naphtha for industry" spatial.oil.kerosene = nodes + " kerosene for aviation" spatial.oil.shipping = nodes + " shipping oil" spatial.oil.agriculture_machinery = nodes + " agriculture machinery oil" + spatial.oil.land_transport = nodes + " land transport oil" else: - spatial.oil.nodes = ["EU oil"] - spatial.oil.locations = ["EU"] + spatial.oil.demand_locations = ["EU"] spatial.oil.naphtha = ["EU naphtha for industry"] spatial.oil.kerosene = ["EU kerosene for aviation"] spatial.oil.shipping = ["EU shipping oil"] spatial.oil.agriculture_machinery = ["EU agriculture machinery oil"] - spatial.oil.land_transport = nodes + " land transport oil" + spatial.oil.land_transport = ["EU land transport oil"] # uranium spatial.uranium = SimpleNamespace() @@ -1588,10 +1595,15 @@ def add_land_transport(n, costs): ice_efficiency = options["transport_internal_combustion_efficiency"] + p_set_land_transport_oil = ice_share / ice_efficiency * transport[nodes].rename(columns=lambda x: x + " land transport oil") + + if not options["regional_oil_demand"]: + p_set_land_transport_oil = p_set_land_transport_oil.sum(axis=1).to_frame(name="EU land transport oil") + n.madd( "Bus", spatial.oil.land_transport, - location=nodes, + location=spatial.oil.demand_locations, carrier="land transport oil", unit="land transport", ) @@ -1601,7 +1613,7 @@ def add_land_transport(n, costs): spatial.oil.land_transport, bus=spatial.oil.land_transport, carrier="land transport oil", - p_set=ice_share / ice_efficiency * transport[nodes].rename(columns=lambda x: x + " land transport oil"), + p_set=p_set_land_transport_oil, ) n.madd( @@ -2638,16 +2650,15 @@ def add_industry(n, costs): options["shipping_oil_efficiency"] / options["shipping_methanol_efficiency"] ) - # need to aggregate potentials if methanol not nodally resolved - if options["co2_budget_national"]: - p_set_methanol = shipping_methanol_share * p_set.rename(lambda x : x + " shipping methanol") * efficiency - else: - p_set_methanol = shipping_methanol_share * p_set.sum() * efficiency + p_set_methanol = shipping_methanol_share * p_set.rename(lambda x : x + " shipping methanol") * efficiency + + if not options["regional_methanol_demand"]: + p_set_methanol = p_set_methanol.sum() n.madd( "Bus", spatial.methanol.shipping, - location=spatial.methanol.locations, + location=spatial.methanol.demand_locations, carrier="shipping methanol", unit="MWh_LHV", ) @@ -2684,7 +2695,8 @@ def add_industry(n, costs): # could correct to e.g. 0.001 EUR/kWh * annuity and O&M n.madd( "Store", - [oil_bus + " Store" for oil_bus in spatial.oil.nodes], + spatial.oil.nodes, + suffix=" Store", bus=spatial.oil.nodes, e_nom_extendable=True, e_cyclic=True, @@ -2702,16 +2714,16 @@ def add_industry(n, costs): ) if shipping_oil_share: - # need to aggregate potentials if oil not nodally resolved - if options["co2_budget_national"]: - p_set_oil = shipping_oil_share * p_set.rename(lambda x: x + " shipping oil") - else: - p_set_oil = shipping_oil_share * p_set.sum() + + p_set_oil = shipping_oil_share * p_set.rename(lambda x: x + " shipping oil") + + if not options["regional_oil_demand"]: + p_set_oil = p_set_oil.sum() n.madd( "Bus", spatial.oil.shipping, - location=spatial.oil.locations, + location=spatial.oil.demand_locations, carrier="shipping oil", unit="MWh_LHV", ) @@ -2781,15 +2793,15 @@ def add_industry(n, costs): if demand_factor != 1: logger.warning(f"Changing HVC demand by {demand_factor*100-100:+.2f}%.") - if options["co2_budget_national"]: - p_set_plastics = demand_factor * industrial_demand.loc[nodes, "naphtha"].rename(lambda x: x + " naphtha for industry") / nhours - else: - p_set_plastics = demand_factor * industrial_demand.loc[nodes, "naphtha"].sum() / nhours + p_set_plastics = demand_factor * industrial_demand.loc[nodes, "naphtha"].rename(lambda x: x + " naphtha for industry") / nhours + + if not options["regional_oil_demand"]: + p_set_plastics = p_set_plastics.sum() n.madd( "Bus", spatial.oil.naphtha, - location=spatial.oil.locations, + location=spatial.oil.demand_locations, carrier="naphtha for industry", unit="MWh_LHV", ) @@ -2826,26 +2838,21 @@ def add_industry(n, costs): logger.warning(f"Changing aviation demand by {demand_factor*100-100:+.2f}%.") all_aviation = ["total international aviation", "total domestic aviation"] - # need to aggregate potentials if oil not nodally resolved - if options["co2_budget_national"]: - p_set = ( + + p_set = ( demand_factor * pop_weighted_energy_totals.loc[nodes, all_aviation].sum(axis=1) * 1e6 / nhours ).rename(lambda x: x + " kerosene for aviation") - else: - p_set = ( - demand_factor - * pop_weighted_energy_totals.loc[nodes, all_aviation].sum(axis=1).sum() - * 1e6 - / nhours - ) + + if not options["regional_oil_demand"]: + p_set = p_set.sum() n.madd( "Bus", spatial.oil.kerosene, - location=spatial.oil.locations, + location=spatial.oil.demand_locations, carrier="kerosene for aviation", unit="MWh_LHV", ) @@ -3111,16 +3118,16 @@ def add_agriculture(n, costs): ) if oil_share > 0: - # need to aggregate potentials if oil not nodally resolved - if options["co2_budget_national"]: - p_set = oil_share * machinery_nodal_energy.rename(lambda x: x + " agriculture machinery oil") / nhours - else: - p_set = oil_share * machinery_nodal_energy.sum() / nhours + + p_set = oil_share * machinery_nodal_energy.rename(lambda x: x + " agriculture machinery oil") / nhours + + if not options["regional_oil_demand"]: + p_set = p_set.sum() n.madd( "Bus", spatial.oil.agriculture_machinery, - location=spatial.oil.locations, + location=spatial.oil.demand_locations, carrier="agriculture machinery oil", unit="MWh_LHV", ) From 1b569dde1bcbcd32175d41b0ba3ed265e76a1aad Mon Sep 17 00:00:00 2001 From: Tom Brown Date: Tue, 2 Jan 2024 16:02:10 +0100 Subject: [PATCH 22/23] move code for national CO2 budgets out of extra_functionality This can be added by derived workflows like PyPSA-Eur via additional_functionality. Changed additional_functionality to pass snakemake rather than wildcards and config separately. This gives maximal flexibility. --- config/config.default.yaml | 17 ------ scripts/solve_network.py | 110 +------------------------------------ 2 files changed, 1 insertion(+), 126 deletions(-) diff --git a/config/config.default.yaml b/config/config.default.yaml index c1e7ed0f9..6d2ebd9f9 100644 --- a/config/config.default.yaml +++ b/config/config.default.yaml @@ -84,22 +84,6 @@ co2_budget: 2045: 0.032 2050: 0.000 -co2_budget_national: - 2030: - 'DE': 0.350 - 'AT': 0.450 - 'BE': 0.450 - 'CH': 0.450 - 'CZ': 0.450 - 'DK': 0.450 - 'FR': 0.450 - 'GB': 0.450 - 'LU': 0.450 - 'NL': 0.450 - 'NO': 0.450 - 'PL': 0.450 - 'SE': 0.450 - # docs in https://pypsa-eur.readthedocs.io/en/latest/configuration.html#electricity electricity: voltages: [220., 300., 380.] @@ -470,7 +454,6 @@ sector: hydrogen_turbine: false SMR: true SMR_cc: true - co2_budget_national: false regional_methanol_demand: false #set to true if regional CO2 constraints needed regional_oil_demand: false #set to true if regional CO2 constraints needed regional_co2_sequestration_potential: diff --git a/scripts/solve_network.py b/scripts/solve_network.py index 6f88b904a..433b175b3 100644 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -41,9 +41,6 @@ pypsa.pf.logger.setLevel(logging.WARNING) from pypsa.descriptors import get_switchable_as_dense as get_as_dense -from prepare_sector_network import emission_sectors_from_opts - - def add_land_use_constraint(n, planning_horizons, config): if "m" in snakemake.wildcards.clusters: @@ -765,100 +762,6 @@ def add_pipe_retrofit_constraint(n): n.model.add_constraints(lhs == rhs, name="Link-pipe_retrofit") -def add_co2limit_country(n, limit_countries, nyears=1.0): - """ - Add a set of emissions limit constraints for specified countries. - - The countries and emissions limits are specified in the config file entry 'co2_budget_country_{investment_year}'. - - Parameters - ---------- - n : pypsa.Network - config : dict - limit_countries : dict - nyears: float, optional - Used to scale the emissions constraint to the number of snapshots of the base network. - """ - logger.info(f"Adding CO2 budget limit for each country as per unit of 1990 levels") - - countries = n.config["countries"] - - # TODO: import function from prepare_sector_network? Move to common place? - sectors = emission_sectors_from_opts(opts) - - # convert Mt to tCO2 - co2_totals = 1e6 * pd.read_csv(snakemake.input.co2_totals_name, index_col=0) - - co2_limit_countries = co2_totals.loc[countries, sectors].sum(axis=1) - co2_limit_countries = co2_limit_countries.loc[co2_limit_countries.index.isin(limit_countries.keys())] - - co2_limit_countries *= co2_limit_countries.index.map(limit_countries) * nyears - - p = n.model["Link-p"] # dimension: (time, component) - - # NB: Most country-specific links retain their locational information in bus1 (except for DAC, where it is in bus2, and process emissions, where it is in bus0) - country = n.links.bus1.map(n.buses.location).map(n.buses.country) - country_DAC = ( - n.links[n.links.carrier == "DAC"] - .bus2.map(n.buses.location) - .map(n.buses.country) - ) - country[country_DAC.index] = country_DAC - country_process_emissions = ( - n.links[n.links.carrier.str.contains("process emissions")] - .bus0.map(n.buses.location) - .map(n.buses.country) - ) - country[country_process_emissions.index] = country_process_emissions - - lhs = [] - for port in [col[3:] for col in n.links if col.startswith("bus")]: - if port == str(0): - efficiency = ( - n.links["efficiency"].apply(lambda x: -1.0).rename("efficiency0") - ) - elif port == str(1): - efficiency = n.links["efficiency"] - else: - efficiency = n.links[f"efficiency{port}"] - mask = n.links[f"bus{port}"].map(n.buses.carrier).eq("co2") - - idx = n.links[mask].index - - international = n.links.carrier.map( - lambda x: 0.4 if x in ["kerosene for aviation", "shipping oil"] else 1.0 - ) - grouping = country.loc[idx] - - if not grouping.isnull().all(): - expr = ( - ((p.loc[:, idx] * efficiency[idx] * international[idx]) - .groupby(grouping, axis=1) - .sum() - *n.snapshot_weightings.generators - ) - .sum(dims="snapshot") - ) - lhs.append(expr) - - lhs = sum(lhs) # dimension: (country) - lhs = lhs.rename({list(lhs.dims.keys())[0]: "snapshot"}) - rhs = pd.Series(co2_limit_countries) # dimension: (country) - - for ct in lhs.indexes["snapshot"]: - n.model.add_constraints( - lhs.loc[ct] <= rhs[ct], - name=f"GlobalConstraint-co2_limit_per_country{ct}", - ) - n.add( - "GlobalConstraint", - f"co2_limit_per_country{ct}", - constant=rhs[ct], - sense="<=", - type="", - ) - - def extra_functionality(n, snapshots): """ Collects supplementary constraints which will be passed to @@ -889,23 +792,12 @@ def extra_functionality(n, snapshots): add_carbon_budget_constraint(n, snapshots) add_retrofit_gas_boiler_constraint(n, snapshots) - if n.config["sector"]["co2_budget_national"]: - # prepare co2 constraint - nhours = n.snapshot_weightings.generators.sum() - nyears = nhours / 8760 - investment_year = int(snakemake.wildcards.planning_horizons[-4:]) - limit_countries = snakemake.config["co2_budget_national"][investment_year] - - # add co2 constraint for each country - logger.info(f"Add CO2 limit for each country") - add_co2limit_country(n, limit_countries, nyears) - if "additional_functionality" in snakemake.input.keys(): import importlib, os, sys sys.path.append(os.path.dirname(snakemake.input.additional_functionality)) additional_functionality = importlib.import_module(os.path.splitext(os.path.basename(snakemake.input.additional_functionality))[0]) - additional_functionality.additional_functionality(n, snapshots, snakemake.wildcards, config) + additional_functionality.additional_functionality(n, snapshots, snakemake) def solve_network(n, config, solving, opts="", **kwargs): From f494dd85b969f9491a2d9bf81ea98008452440a7 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 2 Jan 2024 15:21:49 +0000 Subject: [PATCH 23/23] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- scripts/add_brownfield.py | 31 ++++++++------ scripts/prepare_sector_network.py | 69 ++++++++++++++++++++----------- scripts/solve_network.py | 11 ++++- 3 files changed, 73 insertions(+), 38 deletions(-) diff --git a/scripts/add_brownfield.py b/scripts/add_brownfield.py index fb1453fd7..ffdaf46be 100644 --- a/scripts/add_brownfield.py +++ b/scripts/add_brownfield.py @@ -119,32 +119,39 @@ def add_brownfield(n, n_p, year): n.links.loc[new_pipes, "p_nom"] = 0.0 n.links.loc[new_pipes, "p_nom_min"] = 0.0 + def disable_grid_expansion_if_LV_limit_hit(n): if not "lv_limit" in n.global_constraints.index: return - #calculate minimum LV + # calculate minimum LV attr = "nom_min" dc = n.links.index[n.links.carrier == "DC"] - tot = (n.lines["s_" + attr]*n.lines["length"]).sum() + (n.links.loc[dc,"p_" + attr]*n.links.loc[dc,"length"]).sum() + tot = (n.lines["s_" + attr] * n.lines["length"]).sum() + ( + n.links.loc[dc, "p_" + attr] * n.links.loc[dc, "length"] + ).sum() - diff = n.global_constraints.at["lv_limit","constant"]-tot + diff = n.global_constraints.at["lv_limit", "constant"] - tot - #allow small numerical differences + # allow small numerical differences limit = 1 if diff < limit: - logger.info(f"LV is already reached (gap {diff}), disabling expansion and LV limit") + logger.info( + f"LV is already reached (gap {diff}), disabling expansion and LV limit" + ) expandable_acs = n.lines.index[n.lines.s_nom_extendable] - n.lines.loc[expandable_acs,"s_nom_extendable"] = False - n.lines.loc[expandable_acs,"s_nom"] = n.lines.loc[expandable_acs,"s_nom_min"] + n.lines.loc[expandable_acs, "s_nom_extendable"] = False + n.lines.loc[expandable_acs, "s_nom"] = n.lines.loc[expandable_acs, "s_nom_min"] + + expandable_dcs = n.links.index[ + n.links.p_nom_extendable & (n.links.carrier == "DC") + ] + n.links.loc[expandable_dcs, "p_nom_extendable"] = False + n.links.loc[expandable_dcs, "p_nom"] = n.links.loc[expandable_dcs, "p_nom_min"] - expandable_dcs = n.links.index[n.links.p_nom_extendable & (n.links.carrier == "DC")] - n.links.loc[expandable_dcs,"p_nom_extendable"] = False - n.links.loc[expandable_dcs,"p_nom"] = n.links.loc[expandable_dcs,"p_nom_min"] + n.global_constraints.drop("lv_limit", inplace=True) - n.global_constraints.drop("lv_limit", - inplace=True) if __name__ == "__main__": if "snakemake" not in globals(): diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py index b5a0c0d5d..f1ddce2d6 100644 --- a/scripts/prepare_sector_network.py +++ b/scripts/prepare_sector_network.py @@ -128,8 +128,8 @@ def define_spatial(nodes, options): # methanol - #beware: unlike other carriers, uses locations rather than locations+carriername - #this allows to avoid separation between nodes and locations + # beware: unlike other carriers, uses locations rather than locations+carriername + # this allows to avoid separation between nodes and locations spatial.methanol = SimpleNamespace() @@ -1595,10 +1595,16 @@ def add_land_transport(n, costs): ice_efficiency = options["transport_internal_combustion_efficiency"] - p_set_land_transport_oil = ice_share / ice_efficiency * transport[nodes].rename(columns=lambda x: x + " land transport oil") + p_set_land_transport_oil = ( + ice_share + / ice_efficiency + * transport[nodes].rename(columns=lambda x: x + " land transport oil") + ) if not options["regional_oil_demand"]: - p_set_land_transport_oil = p_set_land_transport_oil.sum(axis=1).to_frame(name="EU land transport oil") + p_set_land_transport_oil = p_set_land_transport_oil.sum(axis=1).to_frame( + name="EU land transport oil" + ) n.madd( "Bus", @@ -2454,7 +2460,7 @@ def add_industry(n, costs): efficiency=1.0, ) - if len(spatial.biomass.industry_cc)<=1 and len(spatial.co2.nodes)>1: + if len(spatial.biomass.industry_cc) <= 1 and len(spatial.co2.nodes) > 1: link_names = nodes + " " + spatial.biomass.industry_cc else: link_names = spatial.biomass.industry_cc @@ -2650,7 +2656,11 @@ def add_industry(n, costs): options["shipping_oil_efficiency"] / options["shipping_methanol_efficiency"] ) - p_set_methanol = shipping_methanol_share * p_set.rename(lambda x : x + " shipping methanol") * efficiency + p_set_methanol = ( + shipping_methanol_share + * p_set.rename(lambda x: x + " shipping methanol") + * efficiency + ) if not options["regional_methanol_demand"]: p_set_methanol = p_set_methanol.sum() @@ -2679,7 +2689,10 @@ def add_industry(n, costs): bus2="co2 atmosphere", carrier="shipping methanol", p_nom_extendable=True, - efficiency2=1 / options["MWh_MeOH_per_tCO2"], # CO2 intensity methanol based on stoichiometric calculation with 22.7 GJ/t methanol (32 g/mol), CO2 (44 g/mol), 277.78 MWh/TJ = 0.218 t/MWh + efficiency2=1 + / options[ + "MWh_MeOH_per_tCO2" + ], # CO2 intensity methanol based on stoichiometric calculation with 22.7 GJ/t methanol (32 g/mol), CO2 (44 g/mol), 277.78 MWh/TJ = 0.218 t/MWh ) if "oil" not in n.buses.carrier.unique(): @@ -2714,7 +2727,6 @@ def add_industry(n, costs): ) if shipping_oil_share: - p_set_oil = shipping_oil_share * p_set.rename(lambda x: x + " shipping oil") if not options["regional_oil_demand"]: @@ -2793,7 +2805,13 @@ def add_industry(n, costs): if demand_factor != 1: logger.warning(f"Changing HVC demand by {demand_factor*100-100:+.2f}%.") - p_set_plastics = demand_factor * industrial_demand.loc[nodes, "naphtha"].rename(lambda x: x + " naphtha for industry") / nhours + p_set_plastics = ( + demand_factor + * industrial_demand.loc[nodes, "naphtha"].rename( + lambda x: x + " naphtha for industry" + ) + / nhours + ) if not options["regional_oil_demand"]: p_set_plastics = p_set_plastics.sum() @@ -2816,7 +2834,10 @@ def add_industry(n, costs): # some CO2 from naphtha are process emissions from steam cracker # rest of CO2 released to atmosphere either in waste-to-energy or decay - process_co2_per_naphtha = industrial_demand.loc[nodes, "process emission from feedstock"].sum() / industrial_demand.loc[nodes, "naphtha"].sum() + process_co2_per_naphtha = ( + industrial_demand.loc[nodes, "process emission from feedstock"].sum() + / industrial_demand.loc[nodes, "naphtha"].sum() + ) emitted_co2_per_naphtha = costs.at["oil", "CO2 intensity"] - process_co2_per_naphtha n.madd( @@ -2840,11 +2861,11 @@ def add_industry(n, costs): all_aviation = ["total international aviation", "total domestic aviation"] p_set = ( - demand_factor - * pop_weighted_energy_totals.loc[nodes, all_aviation].sum(axis=1) - * 1e6 - / nhours - ).rename(lambda x: x + " kerosene for aviation") + demand_factor + * pop_weighted_energy_totals.loc[nodes, all_aviation].sum(axis=1) + * 1e6 + / nhours + ).rename(lambda x: x + " kerosene for aviation") if not options["regional_oil_demand"]: p_set = p_set.sum() @@ -3095,9 +3116,9 @@ def add_agriculture(n, costs): f"Total agriculture machinery shares sum up to {total_share:.2%}, corresponding to increased or decreased demand assumptions." ) - machinery_nodal_energy = pop_weighted_energy_totals.loc[ - nodes, "total agriculture machinery" - ] * 1e6 + machinery_nodal_energy = ( + pop_weighted_energy_totals.loc[nodes, "total agriculture machinery"] * 1e6 + ) if electric_share > 0: efficiency_gain = ( @@ -3111,15 +3132,15 @@ def add_agriculture(n, costs): suffix=" agriculture machinery electric", bus=nodes, carrier="agriculture machinery electric", - p_set=electric_share - / efficiency_gain - * machinery_nodal_energy - / nhours, + p_set=electric_share / efficiency_gain * machinery_nodal_energy / nhours, ) if oil_share > 0: - - p_set = oil_share * machinery_nodal_energy.rename(lambda x: x + " agriculture machinery oil") / nhours + p_set = ( + oil_share + * machinery_nodal_energy.rename(lambda x: x + " agriculture machinery oil") + / nhours + ) if not options["regional_oil_demand"]: p_set = p_set.sum() diff --git a/scripts/solve_network.py b/scripts/solve_network.py index 433b175b3..4bdbb543f 100644 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -793,9 +793,16 @@ def extra_functionality(n, snapshots): add_retrofit_gas_boiler_constraint(n, snapshots) if "additional_functionality" in snakemake.input.keys(): - import importlib, os, sys + import importlib + import os + import sys + sys.path.append(os.path.dirname(snakemake.input.additional_functionality)) - additional_functionality = importlib.import_module(os.path.splitext(os.path.basename(snakemake.input.additional_functionality))[0]) + additional_functionality = importlib.import_module( + os.path.splitext( + os.path.basename(snakemake.input.additional_functionality) + )[0] + ) additional_functionality.additional_functionality(n, snapshots, snakemake)