Repository: PyPSA/PyPSA

Issue:
Improve n.statistics() defaults
Avoid confusion in df aggregation methods. Instead of defining a single default time aggregation method in class `StatisticsAccessor`, use defaults from each metric function. 

Discussion:

Patch:
diff --git a/doc/release_notes.rst b/doc/release_notes.rst
index c5177de2..3a525cb3 100644
--- a/doc/release_notes.rst
+++ b/doc/release_notes.rst
@@ -22,6 +22,7 @@ Upcoming Release
   to change the default snapshot weightings, and ``weightings_from_timedelta``
   to compute the weights if snapshots are of type ``pd.DatetimeIndex``.
 * The statistics function ``n.statistics()`` now also supports the calculation of the ``Market Value`` of components.
+* The statistics function allows now calculating energy balances ``n.statistics.energy_balance()`` and dispatch ``n.statistics.dispatch()``, as well as time series (e.g. ``n.statistics.curtailment(aggregate_time=False)``).
 * Bug fix in linearized unit commitment
 * PyPSA now supports quadratic marginal cost terms. A new column `marginal_cost_quadratic` was added to generators and links. The quadratic marginal cost is added to the objective function when calling ``n.optimize()``.
 * The minimum required version of ``linopy`` is now ``0.2.1``.
diff --git a/pypsa/statistics.py b/pypsa/statistics.py
index 8891c1aa..1009de8b 100644
--- a/pypsa/statistics.py
+++ b/pypsa/statistics.py
@@ -12,6 +12,7 @@ __copyright__ = (
     "MIT License"
 )
 
+import logging
 from functools import wraps
 
 import numpy as np
@@ -19,6 +20,8 @@ import pandas as pd
 
 from pypsa.descriptors import nominal_attrs
 
+logger = logging.getLogger(__name__)
+
 
 def get_carrier(n, c):
     """
@@ -34,17 +37,32 @@ def get_carrier(n, c):
     )
 
 
-def get_bus_and_carrier(n, c):
+def get_bus_and_carrier(n, c, port=""):
     """
     Get the buses and nice carrier names for a component.
     """
-    if "bus" not in n.df(c):
-        bus = "bus0"
+    if port == "":
+        if "bus" not in n.df(c):
+            bus = "bus0"
+        else:
+            bus = "bus"
     else:
-        bus = "bus"
+        bus = f"bus{port}"
     return [n.df(c)[bus].rename("bus"), get_carrier(n, c)]
 
 
+def get_carrier_and_bus_carrier(n, c, port=""):
+    """
+    Get component carrier and bus carrier in one combined DataFrame.
+
+    Used for MultiIndex in energy balance.
+    """
+    bus = f"bus{port}"
+    bus_and_carrier = pd.concat(get_bus_and_carrier(n, c, port), axis=1)
+    bus_carrier = n.df(c)[bus].map(n.buses.carrier).rename("bus_carrier")
+    return pd.concat([bus_and_carrier, bus_carrier], axis=1)
+
+
 def get_operation(n, c):
     """
     Get the operation time series of a component.
@@ -75,6 +93,8 @@ def aggregate_timeseries(df, weights, agg="sum"):
         return df.multiply(weights, axis=0).sum().div(weights.sum())
     elif agg == "sum":
         return weights @ df
+    elif not agg:
+        return df.T
     else:
         return df.agg(agg)
 
@@ -84,6 +104,7 @@ def aggregate_components(n, func, agg="sum", comps=None, groupby=None):
     Apply a function and group the result for a collection of components.
     """
     d = {}
+    kwargs = {}
     if comps is None:
         comps = n.branch_components | n.one_port_components
     if groupby is None:
@@ -95,11 +116,14 @@ def aggregate_components(n, func, agg="sum", comps=None, groupby=None):
             grouping = [n.df(c)[key] for key in groupby]
         elif isinstance(groupby, str):
             grouping = n.df(c)[groupby]
+        elif isinstance(groupby, dict):
+            grouping = None
+            kwargs = groupby
         else:
             ValueError(
-                f"Argument `groupby` must be a function, list or string, got {type(groupby)}"
+                f"Argument `groupby` must be a function, list, string or dict, got {type(groupby)}"
             )
-        d[c] = func(n, c).groupby(grouping).agg(agg)
+        d[c] = func(n, c).groupby(grouping, **kwargs).agg(agg)
     return pd.concat(d)
 
 
@@ -125,9 +149,9 @@ class StatisticsAccessor:
     def __call__(
         self,
         comps=None,
-        aggregate_time="mean",
         aggregate_groups="sum",
         groupby=get_carrier,
+        **kwargs,
     ):
         """
         Calculate statistical values for a network.
@@ -141,17 +165,13 @@ class StatisticsAccessor:
         comps: list-like
             Set of components to consider. Defaults to one-port and branch
             components.
-        aggregate_time : str, optional
-            Type of aggregation when aggregating time series.
-            Note that for {'mean', 'sum'} the time-series are aggregated
-            using snapshot weightings. The default is 'mean'.
         aggregate_groups : str, optional
             Type of aggregation when component groups. The default is 'sum'.
         groupby : callable, list, str, optional
             Specification how to group assets within one component class.
             If a function is passed, it should have the arguments network and
             component name. If a list is passed it should contain
-            column names of the static dataframe, same for a single string.
+            column names of the static DataFrame, same for a single string.
             Defaults to `get_carrier`.
 
         Returns
@@ -159,15 +179,19 @@ class StatisticsAccessor:
         df :
             pandas.DataFrame with columns given the different quantities.
         """
-        static_funcs = [
+        if "aggregate_time" in kwargs:
+            logger.warn(
+                "Argument 'aggregate_time' ignored in overview table. Falling back to individual function defaults."
+            )
+
+        funcs = [
             self.capex,
             self.optimal_capacity,
             self.installed_capacity,
-        ]
-        dynamic_funcs = [
             self.opex,
             self.supply,
             self.withdrawal,
+            self.dispatch,
             self.curtailment,
             self.capacity_factor,
             self.revenue,
@@ -175,30 +199,29 @@ class StatisticsAccessor:
         ]
         kwargs = dict(comps=comps, aggregate_groups=aggregate_groups, groupby=groupby)
         res = []
-        for func in static_funcs:
-            res.append(func(**kwargs))
-        for func in dynamic_funcs:
-            res.append(func(aggregate_time=aggregate_time, **kwargs))
+        for func in funcs:
+            df = func(**kwargs)
+            res.append(df.rename(df.attrs["name"]))
         return pd.concat(res, axis=1).sort_index(axis=0).sort_index(axis=1)
 
-    def get_carrier(self, c):
+    def get_carrier(self, n, c):
         """
         Get the buses and nice carrier names for a component.
         """
-        return get_carrier(self._parent, c)
+        return get_carrier(n, c)
 
-    def get_bus_and_carrier(self, c):
+    def get_bus_and_carrier(self, n, c):
         """
         Get the buses and nice carrier names for a component.
         """
-        return get_bus_and_carrier(self._parent, c)
+        return get_bus_and_carrier(n, c)
 
     def capex(self, comps=None, aggregate_groups="sum", groupby=None):
         """
-        Calculate the capital expenditure of the network.
+        Calculate the capital expenditure of the network in given currency.
 
         For information on the list of arguments, see the docs in
-        `Network.statistics` or `pypsa.statitics.StatisticsAccessor`.
+        `Network.statistics` or `pypsa.statistics.StatisticsAccessor`.
         """
         n = self._parent
 
@@ -209,15 +232,16 @@ class StatisticsAccessor:
         df = aggregate_components(
             n, func, comps=comps, agg=aggregate_groups, groupby=groupby
         )
-        df.attrs["unit"] = "€"
-        return df.rename("Capital Expenditure")
+        df.attrs["name"] = "Capital Expenditure"
+        df.attrs["unit"] = "currency"
+        return df
 
     def optimal_capacity(self, comps=None, aggregate_groups="sum", groupby=None):
         """
-        Calculate the optimal capacity of the network components.
+        Calculate the optimal capacity of the network components in MW.
 
         For information on the list of arguments, see the docs in
-        `Network.statistics` or `pypsa.statitics.StatisticsAccessor`.
+        `Network.statistics` or `pypsa.statistics.StatisticsAccessor`.
         """
         n = self._parent
 
@@ -228,15 +252,16 @@ class StatisticsAccessor:
         df = aggregate_components(
             n, func, comps=comps, agg=aggregate_groups, groupby=groupby
         )
+        df.attrs["name"] = "Optimal Capacity"
         df.attrs["unit"] = "MW"
-        return df.rename("Optimal Capacity")
+        return df
 
     def installed_capacity(self, comps=None, aggregate_groups="sum", groupby=None):
         """
-        Calculate the installed capacity of the network components.
+        Calculate the installed capacity of the network components in MW.
 
         For information on the list of arguments, see the docs in
-        `Network.statistics` or `pypsa.statitics.StatisticsAccessor`.
+        `Network.statistics` or `pypsa.statistics.StatisticsAccessor`.
         """
         n = self._parent
 
@@ -247,12 +272,13 @@ class StatisticsAccessor:
         df = aggregate_components(
             n, func, comps=comps, agg=aggregate_groups, groupby=groupby
         )
+        df.attrs["name"] = "Installed Capacity"
         df.attrs["unit"] = "MW"
-        return df.rename("Installed Capacity")
+        return df
 
     def expanded_capacity(self, comps=None, aggregate_groups="sum", groupby=None):
         """
-        Calculate the expanded capacity of the network components.
+        Calculate the expanded capacity of the network components in MW.
 
         For information on the list of arguments, see the docs in
         `Network.statistics` or `pypsa.statistics.StatisticsAccessor`.
@@ -262,8 +288,9 @@ class StatisticsAccessor:
         ) - self.installed_capacity(
             comps=comps, aggregate_groups=aggregate_groups, groupby=groupby
         )
+        df.attrs["name"] = "Expanded Capacity"
         df.attrs["unit"] = "MW"
-        return df.rename("Expanded Capacity")
+        return df
 
     def opex(
         self,
@@ -273,10 +300,17 @@ class StatisticsAccessor:
         groupby=None,
     ):
         """
-        Calculate the operational expenditure in the network.
+        Calculate the operational expenditure in the network in given currency.
 
         For information on the list of arguments, see the docs in
-        `Network.statistics` or `pypsa.statitics.StatisticsAccessor`.
+        `Network.statistics` or `pypsa.statistics.StatisticsAccessor`.
+
+        Parameters
+        ----------
+        aggregate_time : str, bool, optional
+            Type of aggregation when aggregating time series.
+            Note that for {'mean', 'sum'} the time series are aggregated
+            using snapshot weightings. With False the time series is given. Defaults to 'sum'.
         """
         n = self._parent
 
@@ -293,8 +327,9 @@ class StatisticsAccessor:
         df = aggregate_components(
             n, func, comps=comps, agg=aggregate_groups, groupby=groupby
         )
-        df.attrs["unit"] = "€"
-        return df.rename("Operational Expenditure")
+        df.attrs["name"] = "Operational Expenditure"
+        df.attrs["unit"] = "currency"
+        return df
 
     def supply(
         self,
@@ -304,7 +339,8 @@ class StatisticsAccessor:
         groupby=None,
     ):
         """
-        Calculate the supply of components in the network.
+        Calculate the supply of components in the network. Units depend on the
+        regarded bus carrier.
 
         For information on the list of arguments, see the docs in
         `Network.statistics` or `pypsa.statitics.StatisticsAccessor`.
@@ -324,8 +360,9 @@ class StatisticsAccessor:
         df = aggregate_components(
             n, func, comps=comps, agg=aggregate_groups, groupby=groupby
         )
-        df.attrs["unit"] = "MWh"
-        return df.rename("Supply")
+        df.attrs["name"] = "Supply"
+        df.attrs["unit"] = "carrier dependent"
+        return df
 
     def withdrawal(
         self,
@@ -335,7 +372,8 @@ class StatisticsAccessor:
         groupby=None,
     ):
         """
-        Calculate the withdrawal of components in the network.
+        Calculate the withdrawal of components in the network. Units depend on
+        the regarded bus carrier.
 
         For information on the list of arguments, see the docs in
         `Network.statistics` or `pypsa.statitics.StatisticsAccessor`.
@@ -355,8 +393,104 @@ class StatisticsAccessor:
         df = aggregate_components(
             n, func, comps=comps, agg=aggregate_groups, groupby=groupby
         )
-        df.attrs["unit"] = "MWh"
-        return df.rename("Withdrawal")
+        df.attrs["name"] = "Withdrawal"
+        df.attrs["unit"] = "carrier dependent"
+        return df
+
+    def dispatch(
+        self,
+        comps=None,
+        aggregate_time="sum",
+        aggregate_groups="sum",
+        groupby=None,
+    ):
+        """
+        Calculate the dispatch of components in the network. Units depend on
+        the regarded bus carrier.
+
+        For information on the list of arguments, see the docs in
+        `Network.statistics` or `pypsa.statistics.StatisticsAccessor`.
+
+        Parameters
+        ----------
+        aggregate_time : str, bool, optional
+            Type of aggregation when aggregating time series.
+            Note that for {'mean', 'sum'} the time series are aggregated to MWh
+            using snapshot weightings. With False the time series is given in MW. Defaults to 'sum'.
+        """
+        n = self._parent
+
+        @pass_empty_series_if_keyerror
+        def func(n, c):
+            if c in n.branch_components:
+                p = -n.pnl(c).p0
+            else:
+                p = n.pnl(c).p * n.df(c).sign
+            weights = get_weightings(n, c)
+            return aggregate_timeseries(p, weights, agg=aggregate_time)
+
+        df = aggregate_components(
+            n, func, comps=comps, agg=aggregate_groups, groupby=groupby
+        )
+        df.attrs["name"] = "Dispatch"
+        df.attrs["unit"] = "carrier dependent"
+        return df
+
+    def energy_balance(
+        self,
+        comps=None,
+        aggregate_time="sum",
+        aggregate_groups="sum",
+        aggregate_bus=True,
+    ):
+        """
+        Calculate the energy balance of components in the network. Positive
+        values represent a supply and negative a withdrawal. Units depend on
+        the regarded bus carrier.
+
+        For information on the list of arguments, see the docs in
+        `Network.statistics` or `pypsa.statistics.StatisticsAccessor`.
+
+        Additional parameter
+        ----------
+        aggregate_bus: bool, optional
+            Whether to obtain the nodal or carrier-wise energy balance. Default is True, corresponding to the carrier-wise balance.
+        aggregate_time : str, bool, optional
+            Type of aggregation when aggregating time series.
+            Note that for {'mean', 'sum'} the time series are aggregated to MWh
+            using snapshot weightings. With False the time series is given in MW. Defaults to 'sum'.
+        """
+        n = self._parent
+
+        @pass_empty_series_if_keyerror
+        def func(n, c):
+            sign = -1 if c in n.branch_components else n.df(c).get("sign", 1)
+            ports = [col[3:] for col in n.df(c).columns if col[:3] == "bus"]
+            p = list()
+            for port in ports:
+                mask = n.df(c)[f"bus{port}"] != ""
+                df = sign * n.pnl(c)[f"p{port}"].loc[:, mask]
+                index = get_carrier_and_bus_carrier(n, c, port=port)[mask]
+                df.columns = pd.MultiIndex.from_frame(index.reindex(df.columns))
+                p.append(df)
+            p = pd.concat(p, axis=1)
+            weights = get_weightings(n, c)
+            return aggregate_timeseries(p, weights, agg=aggregate_time)
+
+        groupby = ["carrier", "bus_carrier"]
+        if not aggregate_bus:
+            groupby.append("bus")
+
+        df = aggregate_components(
+            n,
+            func,
+            comps=comps,
+            agg=aggregate_groups,
+            groupby={"level": groupby},
+        )
+        df.attrs["name"] = "Energy Balance"
+        df.attrs["unit"] = "carrier dependent"
+        return df
 
     def curtailment(
         self,
@@ -366,13 +500,20 @@ class StatisticsAccessor:
         groupby=None,
     ):
         """
-        Calculate the curtailment of components in the network.
+        Calculate the curtailment of components in the network in MWh.
 
         The calculation only considers assets with a `p_max_pu` time
         series, which is used to quantify the available power potential.
 
         For information on the list of arguments, see the docs in
-        `Network.statistics` or `pypsa.statitics.StatisticsAccessor`.
+        `Network.statistics` or `pypsa.statistics.StatisticsAccessor`.
+
+        Parameters
+        ----------
+        aggregate_time : str, bool, optional
+            Type of aggregation when aggregating time series.
+            Note that for {'mean', 'sum'} the time series are aggregated to MWh
+            using snapshot weightings. With False the time series is given in MW. Defaults to 'sum'.
         """
         n = self._parent
 
@@ -385,8 +526,9 @@ class StatisticsAccessor:
         df = aggregate_components(
             n, func, comps=comps, agg=aggregate_groups, groupby=groupby
         )
+        df.attrs["name"] = "Curtailment"
         df.attrs["unit"] = "MWh"
-        return df.rename("Curtailment")
+        return df
 
     def capacity_factor(
         self,
@@ -399,7 +541,14 @@ class StatisticsAccessor:
         Calculate the capacity factor of components in the network.
 
         For information on the list of arguments, see the docs in
-        `Network.statistics` or `pypsa.statitics.StatisticsAccessor`.
+        `Network.statistics` or `pypsa.statistics.StatisticsAccessor`.
+
+        Parameters
+        ----------
+        aggregate_time : str, bool, optional
+            Type of aggregation when aggregating time series.
+            Note that for {'mean', 'sum'} the time series are aggregated to
+            using snapshot weightings. With False the time series is given. Defaults to 'mean'.
         """
         n = self._parent
 
@@ -416,8 +565,9 @@ class StatisticsAccessor:
             comps=comps, aggregate_groups=aggregate_groups, groupby=groupby
         )
         df = df.div(capacity, fill_value=np.nan)
+        df.attrs["name"] = "Capacity Factor"
         df.attrs["unit"] = "p.u."
-        return df.rename("Capacity Factor")
+        return df
 
     def revenue(
         self,
@@ -427,10 +577,17 @@ class StatisticsAccessor:
         groupby=None,
     ):
         """
-        Calculate the revenue of components in the network.
+        Calculate the revenue of components in the network in given currency.
 
         For information on the list of arguments, see the docs in
-        `Network.statistics` or `pypsa.statitics.StatisticsAccessor`.
+        `Network.statistics` or `pypsa.statistics.StatisticsAccessor`.
+
+        Parameters
+        ----------
+        aggregate_time : str, bool, optional
+            Type of aggregation when aggregating time series.
+            Note that for {'mean', 'sum'} the time series are aggregated to
+            using snapshot weightings. With False the time series is given. Defaults to 'sum'.
         """
         n = self._parent
 
@@ -452,8 +609,9 @@ class StatisticsAccessor:
         df = aggregate_components(
             n, func, comps=comps, agg=aggregate_groups, groupby=groupby
         )
-        df.attrs["unit"] = "€"
-        return df.rename("Revenue")
+        df.attrs["name"] = "Revenue"
+        df.attrs["unit"] = "currency"
+        return df
 
     def market_value(
         self,
@@ -463,10 +621,18 @@ class StatisticsAccessor:
         groupby=None,
     ):
         """
-        Calculate the market value of components in the network.
+        Calculate the market value of components in the network in given
+        currency/MWh.
 
         For information on the list of arguments, see the docs in
-        `Network.statistics` or `pypsa.statitics.StatisticsAccessor`.
+        `Network.statistics` or `pypsa.statistics.StatisticsAccessor`.
+
+        Parameters
+        ----------
+        aggregate_time : str, bool, optional
+            Type of aggregation when aggregating time series.
+            Note that for {'mean', 'sum'} the time series are aggregated to
+            using snapshot weightings. With False the time series is given. Defaults to 'mean'.
         """
         kwargs = dict(
             comps=comps,
@@ -474,6 +640,7 @@ class StatisticsAccessor:
             aggregate_groups=aggregate_groups,
             groupby=groupby,
         )
-        df = self.revenue(**kwargs) / self.supply(**kwargs)
-        df.attrs["unit"] = "€ / MWh"
-        return df.rename("Market Value")
+        df = self.revenue(**kwargs) / self.dispatch(**kwargs)
+        df.attrs["name"] = "Market Value"
+        df.attrs["unit"] = "currency / MWh"
+        return df


Test patch:
diff --git a/test/test_lopf_global_constraints.py b/test/test_lopf_global_constraints.py
index c37ed76e..ee9858b1 100644
--- a/test/test_lopf_global_constraints.py
+++ b/test/test_lopf_global_constraints.py
@@ -27,7 +27,7 @@ def test_operational_limit_ac_dc_meshed(ac_dc_network, api):
     )
 
     optimize(n, api)
-    assert n.statistics.supply().loc[:, "gas"].sum().round(3) == limit
+    assert n.statistics.dispatch().loc[:, "gas"].sum().round(3) == limit
 
 
 @pytest.mark.parametrize("api", APIS)
diff --git a/test/test_statistics.py b/test/test_statistics.py
index c74a5fc9..cb9638f4 100644
--- a/test/test_statistics.py
+++ b/test/test_statistics.py
@@ -36,6 +36,17 @@ def test_default_solved(ac_dc_network_r):
     df = ac_dc_network_r.statistics.opex()
     assert not df.empty
 
+    df = ac_dc_network_r.statistics.energy_balance()
+    assert not df.empty
+    assert (
+        round(
+            df.groupby(level="bus_carrier").sum().sum()
+            / df.where(lambda x: x > 0).groupby(level="bus_carrier").sum().sum(),
+            3,
+        )
+        == 0
+    )
+
 
 def test_per_bus_carrier_unsolved(ac_dc_network):
     df = ac_dc_network.statistics(groupby=get_bus_and_carrier)