Fix the consumption function

src/muse/quantities.py

@@ -289,54 +289,75 @@ def consumption(
     production: xr.DataArray,
     prices: Optional[xr.DataArray] = None,
     timeslice_level: Optional[str] = None,
-    **kwargs,
 ) -> xr.DataArray:
     """Commodity consumption when fulfilling the whole production.
 
-    Currently, the consumption is implemented for commodity_max == +infinity. If prices
-    are not given, then flexible consumption is *not* considered.
+    Firstly, the degree of technology activity is calculated (i.e. the amount of
+    technology flow required to meet the production). Then, the consumption of fixed
+    commodities is calculated in proportion to this activity.
+
+    In addition, if there are flexible inputs, then the single lowest-cost option is
+    selected (minimising price * quantity). If prices are not given, then flexible
+    consumption is *not* considered.
+
+    Arguments:
+        technologies: Dataset of technology parameters. Must contain `fixed_inputs`,
+            `flexible_inputs`, and `fixed_outputs`.
+        production: DataArray of production data. Must have "timeslice" and "commodity"
+            dimensions.
+        prices: DataArray of prices for each commodity. Must have "timeslice" and
+            "commodity" dimensions. If not given, then flexible inputs are not
+            considered.
+        timeslice_level: the desired timeslice level of the result (e.g. "hour", "day")
+
+    Return:
+        A data array containing the consumption of each commodity. Will have the same
+        dimensions as `production`.
+
     """
-    from muse.commodities import is_enduse, is_fuel
     from muse.utilities import filter_with_template
 
     params = filter_with_template(
-        technologies[["fixed_inputs", "flexible_inputs"]], production, **kwargs
+        technologies[["fixed_inputs", "flexible_inputs", "fixed_outputs"]],
+        production,
     )
 
-    # sum over end-use products, if the dimension exists in the input
-
-    comm_usage = technologies.comm_usage.sel(commodity=production.commodity)
+    # Calculate degree of technology activity
+    prod_amplitude = production_amplitude(
+        production, params, timeslice_level=timeslice_level
+    )
 
-    production = production.sel(commodity=is_enduse(comm_usage)).sum("commodity")
-    params_fuels = is_fuel(params.comm_usage)
-    consumption = production * broadcast_timeslice(
-        params.fixed_inputs.where(params_fuels, 0), level=timeslice_level
+    # Calculate consumption of fixed commodities
+    consumption_fixed = prod_amplitude * broadcast_timeslice(
+        params.fixed_inputs, level=timeslice_level
     )
+    assert all(consumption_fixed.commodity.values == production.commodity.values)
 
+    # If there are no flexible inputs, then we are done
+    if not (params.flexible_inputs > 0).any():
+        return consumption_fixed
+
+    # If prices are not given, then we can't consider flexible inputs, so just return
+    # the fixed consumption
     if prices is None:
-        return consumption
-
-    if not (params.flexible_inputs.sel(commodity=params_fuels) > 0).any():
-        return consumption
-
-    prices = filter_with_template(prices, production, installed_as_year=False, **kwargs)
-    # technology with flexible inputs
-    flexs = params.flexible_inputs.where(params_fuels, 0)
-    # cheapest fuel for each flexible technology
-    assert prices is not None
-    flexprice = [i for i in flexs.commodity.values if i in prices.commodity.values]
-    assert all(flexprice)
-    priceflex = prices.loc[dict(commodity=flexs.commodity)]
+        return consumption_fixed
+
+    # Flexible inputs
+    flexs = broadcast_timeslice(params.flexible_inputs, level=timeslice_level)
+
+    # Calculate the cheapest fuel for each flexible technology
+    priceflex = prices * flexs
     minprices = flexs.commodity[
         priceflex.where(flexs > 0, priceflex.max() + 1).argmin("commodity")
     ]
-    # add consumption from cheapest fuel
-    assert all(flexs.commodity.values == consumption.commodity.values)
+
+    # Consumption of flexible commodities
+    assert all(flexs.commodity.values == consumption_fixed.commodity.values)
     flex = flexs.where(
-        minprices == broadcast_timeslice(flexs.commodity, level=timeslice_level), 0
+        broadcast_timeslice(flexs.commodity, level=timeslice_level) == minprices, 0
     )
-    flex = flex / (flex > 0).sum("commodity").clip(min=1)
-    return consumption + flex * production
+    consumption_flex = flex * prod_amplitude
+    return consumption_fixed + consumption_flex
 
 
 def maximum_production(
@@ -366,8 +387,6 @@ def maximum_production(
         technologies: xr.Dataset describing the features of the technologies of
             interests.  It should contain `fixed_outputs` and `utilization_factor`. It's
             shape is matched to `capacity` using `muse.utilities.broadcast_techs`.
-        timeslices: xr.DataArray of the timeslicing scheme. Production data will be
-            returned in this format.
         filters: keyword arguments are used to filter down the capacity and
             technologies. Filters not relevant to the quantities of interest, i.e.
             filters that are not a dimension of `capacity` or `technologies`, are
@@ -531,4 +550,40 @@ def capacity_to_service_demand(
         level=timeslice_level,
     ) * broadcast_timeslice(technologies.utilization_factor, level=timeslice_level)
     capa_to_service_demand = demand / timeslice_outputs
-    return capa_to_service_demand.max(("commodity", "timeslice"))
+    return capa_to_service_demand.where(np.isfinite(capa_to_service_demand), 0).max(
+        ("commodity", "timeslice")
+    )
+
+
+def production_amplitude(
+    production: xr.DataArray,
+    technologies: xr.Dataset,
+    timeslice_level: Optional[str] = None,
+) -> xr.DataArray:
+    """Calculates the degree of technology activity based on production data.
+
+    We do this by dividing the production data by the output flow per unit of activity.
+    Taking the max of this across all commodities, we get the minimum units of
+    technology activity required to meet (at least) the specified production of all
+    commodities.
+
+    For example:
+    A technology has the following reaction: 1A -> 2B + 3C
+    If production is 4B & 6C, this is equal to a production amplitude of 2
+
+    Args:
+        production: DataArray with commodity-level production for a set of technologies.
+            Must have `timeslice` and `commodity` dimensions. May also have other
+            dimensions e.g. `region`, `year`, etc.
+        technologies: Dataset of technology parameters
+        timeslice_level: the desired timeslice level of the result (e.g. "hour", "day").
+            Must match the timeslice level of `production`
+
+    Returns:
+        DataArray with production amplitudes for each technology in each timeslice.
+        Will have the same dimensions as `production`, minus the `commodity` dimension.
+    """
+    return (
+        production
+        / broadcast_timeslice(technologies.fixed_outputs, level=timeslice_level)
+    ).max("commodity")

tests/test_objectives.py

@@ -1,14 +1,15 @@
 from pytest import fixture, mark
 
+YEAR = 2030
+
 
 @fixture
 def _technologies(technologies, retro_agent, search_space):
     techs = retro_agent.filter_input(
         technologies,
         technology=search_space.replacement,
-        year=retro_agent.forecast_year,
     ).drop_vars("technology")
-    return techs
+    return techs.sel(year=YEAR)
 
 
 @fixture
@@ -26,14 +27,14 @@ def _demand(demand_share, search_space):
 @fixture
 def _prices(retro_agent, agent_market):
     prices = retro_agent.filter_input(agent_market.prices)
-    return prices
+    return prices.sel(year=YEAR)
 
 
 def test_fixtures(_technologies, _demand, _prices):
     """Validating that the fixtures have appropriate dimensions."""
     assert set(_technologies.dims) == {"commodity", "replacement"}
     assert set(_demand.dims) == {"asset", "commodity", "timeslice"}
-    assert set(_prices.dims) == {"commodity", "timeslice", "year"}
+    assert set(_prices.dims) == {"commodity", "timeslice"}
 
 
 @mark.usefixtures("save_registries")
@@ -167,7 +168,7 @@ def test_emission_cost(_technologies, _demand, _prices):
     assert set(result.dims) == {"replacement", "asset", "timeslice"}
 
 
-def test_fuel_consumption(_technologies, _demand, _prices):
+def test_fuel_consumption_cost(_technologies, _demand, _prices):
     from muse.objectives import fuel_consumption_cost
 
     result = fuel_consumption_cost(_technologies, _demand, _prices)