From d36b248acf69bd280bfdc70d92e384ad88c45fe5 Mon Sep 17 00:00:00 2001 From: Tom Bland Date: Fri, 13 Dec 2024 15:42:56 +0000 Subject: [PATCH 1/4] Copy changes to quantities module from #556 --- src/muse/quantities.py | 119 +++++++++++++++++++-------- tests/test_quantities.py | 172 +++++---------------------------------- 2 files changed, 107 insertions(+), 184 deletions(-) diff --git a/src/muse/quantities.py b/src/muse/quantities.py index 27ccfba3c..74e85e317 100644 --- a/src/muse/quantities.py +++ b/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,38 @@ 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. + + 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. + """ + assert set(technologies.dims).issubset(set(production.dims)) + + return ( + production + / broadcast_timeslice(technologies.fixed_outputs, level=timeslice_level) + ).max("commodity") diff --git a/tests/test_quantities.py b/tests/test_quantities.py index 9f89f16b4..ed5e09525 100644 --- a/tests/test_quantities.py +++ b/tests/test_quantities.py @@ -9,16 +9,13 @@ def production( technologies: xr.Dataset, capacity: xr.DataArray, timeslice ) -> xr.DataArray: - from numpy.random import random - from muse.timeslices import broadcast_timeslice, distribute_timeslice - comms = xr.DataArray( - random(len(technologies.commodity)), - coords={"commodity": technologies.commodity}, - dims="commodity", + return ( + broadcast_timeslice(capacity) + * distribute_timeslice(technologies.fixed_outputs) + * broadcast_timeslice(technologies.utilization_factor) ) - return broadcast_timeslice(capacity) * distribute_timeslice(comms) def test_supply_enduse(technologies, capacity, timeslice): @@ -116,109 +113,21 @@ def test_decommissioning_demand(technologies, capacity, timeslice): ).values == approx(ufac * fouts * (current - forecast)) -def test_consumption_no_flex(technologies, production, market): - from muse.commodities import is_enduse, is_fuel +def test_consumption(technologies, production, market): from muse.quantities import consumption - from muse.timeslices import broadcast_timeslice - - fins = ( - technologies.fixed_inputs.where(is_fuel(technologies.comm_usage), 0) - .interp(year=sorted(set(production.installed.values)), method="slinear") - .sel( - technology=production.technology, - region=production.region, - year=production.installed, - ) - ) - services = technologies.commodity.sel(commodity=is_enduse(technologies.comm_usage)) - expected = ( - (production.rename(commodity="comm_in") * broadcast_timeslice(fins)) - .sel(comm_in=production.commodity.isin(services).rename(commodity="comm_in")) - .sum("comm_in") - ) - actual = consumption(technologies, production) - assert set(actual.dims) == set(expected.dims) - assert actual.values == approx(expected.values) + # Prices not provided, so flexible inputs are ignored + consump = consumption(technologies, production) + assert set(production.dims) == set(consump.dims) + # Prices provided, but no flexible inputs -> should be the same as above technologies.flexible_inputs[:] = 0 - actual = consumption(technologies, production, market.prices) - assert actual.values == approx(expected.values) - - -def test_consumption_with_flex(technologies, production, market, timeslice): - from itertools import product + consump2 = consumption(technologies, production, market.prices) + assert consump2.values == approx(consump.values) - from muse.commodities import is_enduse, is_fuel - from muse.quantities import consumption - from muse.timeslices import broadcast_timeslice, distribute_timeslice - - techs = technologies.copy() - techs.fixed_inputs[:] = 0 - techs.flexible_inputs[:] = 0 - consumables = is_fuel(techs.comm_usage) - while (techs.flexible_inputs.sel(commodity=consumables) == 0).all(): - techs.flexible_inputs[:] = ( - np.random.randint(0, 2, techs.flexible_inputs.shape) != 0 - ) - techs.flexible_inputs[{"commodity": ~consumables}] = 0 - - def one_dim(dimension): - from numpy import arange - from numpy.random import shuffle - - data = arange(len(dimension), dtype="int") - shuffle(data) - return xr.DataArray(data, coords=dimension.coords, dims=dimension.dims) - - year = one_dim(production.year) - asset = one_dim(production.asset) - region = one_dim(market.region) - timeslice = one_dim(market.timeslice) - commodity = one_dim(market.commodity) - - prices = ( - timeslice - + broadcast_timeslice(commodity) - + broadcast_timeslice(year) * broadcast_timeslice(region) - ) - assert set(prices.dims) == set(market.prices.dims) - noenduse = ~is_enduse(techs.comm_usage) - production = distribute_timeslice(asset * year + commodity) - production.loc[{"commodity": noenduse}] = 0 - - actual = consumption(technologies, production, prices) - assert set(actual.dims) == {"year", "timeslice", "asset", "commodity"} - assert (year.year == actual.year).all() - assert (timeslice.timeslice == actual.timeslice).all() - assert (asset.asset == actual.asset).all() - assert (commodity.commodity == actual.commodity).all() - - fuels = techs.commodity.loc[{"commodity": consumables}].values - dims = ("timeslice", "asset", "year") - allprods = list(product(*(actual[u] for u in dims))) - allprods = [ - allprods[i] for i in np.random.choice(range(len(allprods)), 50, replace=False) - ] - for ts, asset, year in allprods: - flexs = techs.flexible_inputs.sel( - region=asset.region, technology=asset.technology - ).interp(year=asset.installed, method="slinear") - comm_prices = prices.sel(region=asset.region, year=year, timeslice=ts) - comm_prices = [int(p) for p, f in zip(comm_prices, flexs) if f > 0] - min_price = min(comm_prices) if comm_prices else None - ncomms = max(len([u for u in comm_prices if u == min_price]), 1) - for comm in fuels: - current_price = prices.sel( - region=asset.region, year=year, timeslice=ts, commodity=comm - ) - coords = dict(timeslice=ts, year=year, asset=asset, commodity=comm) - if current_price != min_price: - assert actual.sel(coords).values == approx(0) - continue - prod = production.sel(asset=asset, year=year).sum("commodity") - expected = prod.sel(timeslice=ts) / ncomms * flexs.sel(commodity=comm) - assert expected.values == approx(actual.sel(coords).values) + # Flexible inputs considered + consump3 = consumption(technologies, production, market.prices) + assert set(production.dims) == set(consump3.dims) def test_production_aggregate_asset_view( @@ -307,52 +216,6 @@ def test_capacity_in_use(production: xr.DataArray, technologies: xr.Dataset): assert capa.values == approx(prod / fout / ufac) -def test_supply_cost(production: xr.DataArray, timeslice: xr.Dataset): - from numpy import average - from numpy.random import random - - from muse.costs import supply_cost - - timeslice = timeslice.timeslice - production = production.sel(year=production.year.min(), drop=True) - # no zero production, because it does not sit well with np.average - production[:] = random(production.shape) - lcoe = xr.DataArray( - random((len(production.asset), len(timeslice))), - coords={"timeslice": timeslice, "asset": production.asset}, - dims=("asset", "timeslice"), - ) - - production, lcoe = xr.broadcast(production, lcoe) - actual = supply_cost(production, lcoe, asset_dim="asset") - for region in set(production.region.values): - expected = average( - lcoe.sel(asset=production.region == region), - weights=production.sel(asset=production.region == region), - axis=production.get_axis_num("asset"), - ) - - assert actual.sel(region=region).values == approx(expected) - - -def test_supply_cost_zero_prod(production: xr.DataArray, timeslice: xr.Dataset): - from numpy.random import randn - - from muse.costs import supply_cost - - timeslice = timeslice.timeslice - production = production.sel(year=production.year.min(), drop=True) - production[:] = 0 - lcoe = xr.DataArray( - randn(len(production.asset), len(timeslice)), - coords={"timeslice": timeslice, "asset": production.asset}, - dims=("asset", "timeslice"), - ) - production, lcoe = xr.broadcast(production, lcoe) - actual = supply_cost(production, lcoe, asset_dim="asset") - assert actual.values == approx(0e0) - - def test_emission(production: xr.DataArray, technologies: xr.Dataset): from muse.commodities import is_enduse, is_pollutant from muse.quantities import emission @@ -423,3 +286,10 @@ def test_supply_capped_by_min_service(technologies, capacity, timeslice): ) assert (spl == approx(demand)).all() assert (spl <= minprod).all() + + +def test_production_amplitude(production, technologies): + from muse.quantities import production_amplitude + + result = production_amplitude(production, technologies) + assert set(result.dims) == set(production.dims) - {"commodity"} From bfe81eea5eb127e675aac1c714ea22926d87fc50 Mon Sep 17 00:00:00 2001 From: Tom Bland Date: Fri, 13 Dec 2024 15:51:54 +0000 Subject: [PATCH 2/4] Fix tests --- src/muse/quantities.py | 2 +- tests/test_objectives.py | 11 ++++++----- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/src/muse/quantities.py b/src/muse/quantities.py index 74e85e317..85fbfc21d 100644 --- a/src/muse/quantities.py +++ b/src/muse/quantities.py @@ -579,7 +579,7 @@ def production_amplitude( DataArray with production amplitudes for each technology in each timeslice. Will have the same dimensions as `production`, minus the `commodity` dimension. """ - assert set(technologies.dims).issubset(set(production.dims)) + # assert set(technologies.dims).issubset(set(production.dims)) return ( production diff --git a/tests/test_objectives.py b/tests/test_objectives.py index 593877f4d..5727df2d2 100644 --- a/tests/test_objectives.py +++ b/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) From fb2eba15b78aac464dc2e6f9104aaf30a8ef8a57 Mon Sep 17 00:00:00 2001 From: Tom Bland Date: Fri, 13 Dec 2024 15:57:31 +0000 Subject: [PATCH 3/4] Remove line --- src/muse/quantities.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/muse/quantities.py b/src/muse/quantities.py index 85fbfc21d..5d424b266 100644 --- a/src/muse/quantities.py +++ b/src/muse/quantities.py @@ -579,8 +579,6 @@ def production_amplitude( DataArray with production amplitudes for each technology in each timeslice. Will have the same dimensions as `production`, minus the `commodity` dimension. """ - # assert set(technologies.dims).issubset(set(production.dims)) - return ( production / broadcast_timeslice(technologies.fixed_outputs, level=timeslice_level) From 60c6bdcea0f8679be44e679f445cffc42ea9507d Mon Sep 17 00:00:00 2001 From: Tom Bland Date: Tue, 17 Dec 2024 10:41:07 +0000 Subject: [PATCH 4/4] Add simple example to docstring --- src/muse/quantities.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/muse/quantities.py b/src/muse/quantities.py index 5d424b266..23ce7e409 100644 --- a/src/muse/quantities.py +++ b/src/muse/quantities.py @@ -567,6 +567,10 @@ def production_amplitude( 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