Skip to content

Add carbon emissions to statistics module#1149

Open
gincrement wants to merge 19 commits intoPyPSA:masterfrom
gincrement:carbon_emission_stats
Open

Add carbon emissions to statistics module#1149
gincrement wants to merge 19 commits intoPyPSA:masterfrom
gincrement:carbon_emission_stats

Conversation

@gincrement
Copy link
Copy Markdown
Contributor

@gincrement gincrement commented Feb 15, 2025 •
edited by lkstrp
Loading

Closes #520

Changes proposed in this Pull Request

Checklist

  • Code changes are sufficiently documented; i.e. new functions contain docstrings and further explanations may be given in doc.
  • Unit tests for new features were added (if applicable).
  • A note for the release notes doc/release_notes.rst of the upcoming release is included.
  • I consent to the release of this PR's code under the MIT license.

@lkstrp lkstrp changed the title Add carbon emissions to statistics module (Issue #520) Add carbon emissions to statistics module Feb 27, 2025
lkstrp
lkstrp requested changes Feb 27, 2025
View reviewed changes
Copy link
Copy Markdown
Member

@lkstrp lkstrp left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for another PR @gincrement !
Could you have a look on all failing checks?

@fneum fneum requested a review from Copilot March 3, 2025 09:30
Copilot AI reviewed Mar 3, 2025
View reviewed changes
Copy link
Copy Markdown

Copilot AI left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

PR Overview

This PR adds a new carbon emissions calculation to the statistics module.

  • Added a new function carbon_emissions in the statistics expressions module.
  • Updated the overall statistics aggregation list to include carbon_emissions.
  • Extended the test suite and example notebook to demonstrate the new functionality.

Reviewed Changes

File Description
pypsa/statistics/expressions.py Adds the carbon_emissions function and refines port_efficiency.
test/test_statistics.py Introduces a test for the new carbon emissions calculation.
examples/notebooks/statistics.ipynb Updates example notebook to include carbon emissions example output.

Copilot reviewed 3 out of 3 changed files in this pull request and generated 2 comments.

Comments suppressed due to low confidence (1)

pypsa/statistics/expressions.py:1098

  • The use of 'n.df(c)' appears non-standard compared to other parts of the module where static data is accessed. Verify that 'n.df(c)' is the intended accessor for carrier information, or consider using a consistent method such as 'n.static(c)'.
* n.df(c).carrier.map(n.carriers.co2_emissions)

FabianHofmann
FabianHofmann requested changes Mar 3, 2025
View reviewed changes
Copy link
Copy Markdown
Contributor

@FabianHofmann FabianHofmann left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

nice idea to add this!

fneum
fneum requested changes Mar 4, 2025
View reviewed changes
Copy link
Copy Markdown
Member

@fneum fneum left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks mostly good apart from the comments raised already and one further question from me:

Is it really necessary to use the absolute value function for get_operation() and efficiency? Imagine the case where you're modelling an integrated BECCS power plant (without modelling the CO2 carrier specifically). Would the function then show net-negative emissions? Perhaps add a test for that case.

lkstrp
lkstrp reviewed Mar 6, 2025
View reviewed changes
FabianHofmann
FabianHofmann reviewed Mar 6, 2025
View reviewed changes
@lkstrp lkstrp mentioned this pull request Mar 24, 2025
Merged
21 tasks
@lkstrp lkstrp added this to the v0.35 milestone Apr 22, 2025
@FabianHofmann
Copy link
Copy Markdown
Contributor

FabianHofmann commented May 5, 2025

after notifying @gincrement that changing the port_efficiency function would lead to quite some problems in the other statistic functions, I added the efficiency calculation to the carbon_emission function directly. Note that the current implementation does only support generators, it does not include Stores and StorageUnits as imposed by the corresponding constraint. I honestly don't know if we can/should include them. probably it would be nice and add them but I had to think about it more...

fneum
fneum requested changes Jun 6, 2025
View reviewed changes
Copy link
Copy Markdown
Member

@fneum fneum left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think it should match the corresponding global constraint:

https://pypsa.readthedocs.io/en/latest/user-guide/optimal-power-flow.html#primary-energy

e.g. looking at the net supply of the store/storage unit between initial and final state of charge.

Just returning the carbon emissions of the generators could be confusuing and would not be aligned with the global constraint. That should be avoided.

Even though the statistics function prints a warning, I think the stores and storage units implementation needs to be included before merging.

@lkstrp lkstrp removed this from the v0.35 milestone Jun 23, 2025
@gincrement
Update expressions.py
c57f6e6 
removed abs from efficiency and get_operation and increased coverate to all technologies
@gincrement
Update test_statistics.py
9b906a5 
added co2 carrier, bus and load to test negative and positive flows
fneum
fneum requested changes Dec 18, 2025
View reviewed changes
Copy link
Copy Markdown
Member

@fneum fneum left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Another unit test would be good where emissions come from the net difference in state of charge of a Store / StorageUnit. Then, I'm fine for this to be merged.

See

PyPSA/pypsa/optimization/global_constraints.py

Lines 272 to 480 in 5031132

def define_primary_energy_limit(n: Network, sns: pd.Index) -> None:
"""Define primary energy constraints.
It limits the byproducts of primary energy sources (defined by carriers) such
as CO2.
Parameters
----------
n : pypsa.Network
The network to apply constraints to.
sns : list-like
Set of snapshots to which the constraint should be applied.
"""
m = n.model
weightings = n.snapshot_weightings.loc[sns]
glcs = n.c.global_constraints.static.query('type == "primary_energy"')
if n._multi_invest:
period_weighting = n.investment_period_weightings.years[sns.unique("period")]
weightings = weightings.mul(period_weighting, level=0, axis=0)
period_last_sns = pd.MultiIndex.from_frame(
sns.to_frame(index=False).groupby("period").timestep.last().reset_index()
)
storage_weightings = (
pd.Series(1, n.snapshots).mul(period_weighting).loc[period_last_sns]
)
unique_names = glcs.index.unique("name")
for name in unique_names:
if n.has_scenarios:
glc_group = glcs.xs(name, level="name")
scenarios = glc_group.index.get_level_values("scenario")
else:
glc_group = glcs.loc[name]
scenarios = [slice(None)]
expressions = []
for scenario in scenarios:
glc = glc_group.loc[scenario]
if isnan(glc.investment_period):
sns_sel = slice(None)
elif glc.investment_period in sns.unique("period"):
sns_sel = sns.get_loc(glc.investment_period)
else:
continue
lhs = []
emissions = n.c.carriers.static[glc.carrier_attribute][
lambda ds: ds != 0
].loc[scenario]
if emissions.empty:
continue
# Determine investment period for active asset filtering
period = glc.investment_period if n._multi_invest else None
# generators
gens = n.c.generators.static[
n.c.generators.static.carrier.isin(emissions.index)
]
gens = n.c.generators.filter_by_active_assets(gens, period)
if not gens.empty:
gens = gens.loc[scenario]
efficiency = (
n.c.generators._as_dynamic("efficiency")
.loc[:, scenario]
.loc[sns[sns_sel], gens.index]
)
em_pu = gens.carrier.map(emissions) / efficiency
em_pu = em_pu.multiply(weightings.generators[sns_sel], axis=0)
p = m["Generator-p"].sel(name=gens.index, snapshot=sns[sns_sel])
if n.has_scenarios:
p = p.sel(scenario=scenario, drop=True)
expr = (p * em_pu).sum()
lhs.append(expr)
# storage units
cond = "carrier in @emissions.index and not cyclic_state_of_charge"
sus = n.c.storage_units.static.query(cond)
sus = n.c.storage_units.filter_by_active_assets(sus, period)
if not sus.empty:
sus = sus.loc[scenario]
em_pu = sus.carrier.map(emissions)
soc = m["StorageUnit-state_of_charge"].sel(
name=sus.index, snapshot=sns[sns_sel]
)
if n._multi_invest:
sus_continuous = sus.query("not state_of_charge_initial_per_period")
if not sus_continuous.empty and period_weighting.ne(1).any():
msg = (
"Found non-cyclic storage units with associated carrier emissions "
"and continuous depletion over multiple investment periods "
"combined with investment period year weightings != 1. "
"The primary energy constraint will be inconsistent. "
"Please consider setting `state_of_charge_initial_per_period` to True, "
"using equal period weightings or a cyclic storage unit instead."
)
raise NotImplementedError(msg)
if not sus_continuous.empty and period_weighting.eq(1).all():
soc_final = (
soc.sel(name=sus_continuous.index)
.ffill("snapshot")
.isel(snapshot=-1)
)
if n.has_scenarios:
soc_final = soc_final.sel(scenario=scenario, drop=True)
lhs.append(
(soc_final * -em_pu).sum()
+ em_pu @ sus_continuous.state_of_charge_initial
)
sus_per_period = sus.query("state_of_charge_initial_per_period")
if not sus_per_period.empty:
soc_final = soc.loc[period_last_sns, sus_per_period.index]
if n.has_scenarios:
soc_final = soc_final.sel(scenario=scenario, drop=True)
soc_delta = -soc_final + sus_per_period.state_of_charge_initial
lhs.append((soc_delta * storage_weightings * em_pu).sum())
else:
soc_final = soc.ffill("snapshot").isel(snapshot=-1)
if n.has_scenarios:
soc_final = soc_final.sel(scenario=scenario, drop=True)
lhs.append(
(soc_final * -em_pu).sum() + em_pu @ sus.state_of_charge_initial
)
# stores
stores = n.c.stores.static.query(
"carrier in @emissions.index and not e_cyclic"
)
stores = n.c.stores.filter_by_active_assets(stores, period)
if not stores.empty:
stores = stores.loc[scenario]
em_pu = stores.carrier.map(emissions)
e = m["Store-e"].sel(name=stores.index, snapshot=sns[sns_sel])
if n._multi_invest:
stores_continuous = stores.query("not e_initial_per_period")
if not stores_continuous.empty and period_weighting.ne(1).any():
msg = (
"Found non-cyclic stores with associated carrier emissions "
"and continuous depletion over multiple investment periods "
"combined with investment period year weightings != 1. "
"The primary energy constraint will be inconsistent. "
"Please consider setting `e_initial_per_period` to True, "
"using equal period weightings or a cyclic store instead."
)
raise NotImplementedError(msg)
if not stores_continuous.empty and period_weighting.eq(1).all():
e_final = (
e.sel(name=stores_continuous.index)
.ffill("snapshot")
.isel(snapshot=-1)
)
if n.has_scenarios:
e_final = e_final.sel(scenario=scenario, drop=True)
lhs.append(
(e_final * -em_pu).sum()
+ em_pu @ stores_continuous.e_initial
)
stores_per_period = stores.query("e_initial_per_period")
if not stores_per_period.empty:
e_final = e.loc[period_last_sns, stores_per_period.index]
if n.has_scenarios:
e_final = e_final.sel(scenario=scenario, drop=True)
e_delta = -e_final + stores_per_period.e_initial
lhs.append((e_delta * storage_weightings * em_pu).sum())
else:
e_final = e.ffill("snapshot").isel(snapshot=-1)
if n.has_scenarios:
e_final = e_final.sel(scenario=scenario, drop=True)
lhs.append((e_final * -em_pu).sum() + em_pu @ stores.e_initial)
if not lhs:
continue
lhs = merge(lhs)
expressions.append(lhs)
if not expressions:
continue
if n.has_scenarios:
expression = merge(expressions, dim="scenario").assign_coords(
scenario=scenarios
)
else:
expression = expressions[0]
m.add_constraints(
expression,
glc_group.sense,
glc_group.constant,
name=f"GlobalConstraint-{name}",
)

test/test_statistics.py Outdated
Comment on lines +308 to +315
# add some co2 details
n.add("Carrier", "co2", co2_emissions=1)
n.add("Bus", "co2bus", carrier="co2", color="darkblue")
n.add("Load", "co2", bus="co2bus", carrier="co2")
n.dynamic('Load').p_set["co2"] = 0.0
n.loads_t.p_set.at[n.loads_t.p_set.index[-1], 'co2'] = 10
n.add("Generator", "co2purchase", bus="co2bus", carrier="co2", p_nom_extendable=True, capital_cost=10000, marginal_price=1, p_nom_max=2, efficiency=-1)
n.add("StorageUnit", "CO2storage", bus="co2bus", carrier="co2", p_nom_extendable=True, capital_cost=10000, standing_loss=0.1)
Copy link
Copy Markdown
Member

@fneum fneum Dec 18, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This should not be necessary. The ac_dc_network alone should produce a valid n.statistics.carbon_emissions() output, since one of its carriers ('gas') has specific emissions.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Reviewers

@FabianHofmann FabianHofmann FabianHofmann requested changes

@lkstrp lkstrp lkstrp requested changes

Copilot code review Copilot Copilot left review comments

@fneum fneum fneum requested changes

Requested changes must be addressed to merge this pull request.

Assignees

No one assigned

Labels

None yet

Projects

None yet

Milestone

No milestone

Development

Successfully merging this pull request may close these issues.

Add carbon intensity to statistics module

5 participants

@gincrement @FabianHofmann @fneum @lkstrp