From 3575dc2b3bc6cb07f5ecefa8134fba111f55a9a7 Mon Sep 17 00:00:00 2001 From: John Jasa Date: Thu, 13 Aug 2026 22:34:06 -0600 Subject: [PATCH 1/5] First pass at hetereo-commodity handling --- CHANGELOG.md | 1 + .../system_level/cost_minimization_control.py | 55 +-- .../system_level/demand_following_control.py | 61 +--- .../profit_maximization_control.py | 70 +--- .../system_level/system_level_control_base.py | 316 ++++++++++++++++++ .../system_level/test/test_slc_controllers.py | 310 +++++++++++++++++ 6 files changed, 655 insertions(+), 158 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index b5e276524..39fc54ee8 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,7 @@ - Enable `PySAMWindPlantPerformanceModel` to accept more than 300 turbines by overriding the default maximum in the PySAM model. [PR 831](https://github.com/NatLabRockies/H2Integrate/pull/831) - Add `PySAMWavePerformanceModel` and `WaveResource` to wrap PySAM MhkWave as an H2I performance model, replacing the HOPP wave module in example 09. [PR 825](https://github.com/NatLabRockies/H2Integrate/pull/825) - Replace HOPP with native H2I wind, solar, and battery models in example 11. Adds `percent_load_missed` and `curtailment_percent` outputs to `DemandComponentBase`, allows zero capacity in wind/solar/battery performance models. [PR 826](https://github.com/NatLabRockies/H2Integrate/pull/826) +- Add heterogeneous-commodity system-level control that translates demand for one commodity into upstream set-points across converters using static per-technology conversion ratios defined in the tech config. [PR TBD](https://github.com/NatLabRockies/H2Integrate/pull/TBD) ## 0.9 [August 10, 2026] diff --git a/h2integrate/control/control_strategies/system_level/cost_minimization_control.py b/h2integrate/control/control_strategies/system_level/cost_minimization_control.py index 869976d50..66e5771c5 100644 --- a/h2integrate/control/control_strategies/system_level/cost_minimization_control.py +++ b/h2integrate/control/control_strategies/system_level/cost_minimization_control.py @@ -1,5 +1,3 @@ -import numpy as np - from h2integrate.control.control_strategies.system_level.system_level_control_base import ( SystemLevelControlBase, ) @@ -42,53 +40,8 @@ def setup(self): self._setup_marginal_costs() def compute(self, inputs, outputs): - demand = inputs[self.demand_input_name].copy() - - # 1. Fixed techs: always produce, subtract from demand - for fixed_tech in self.fixed_techs: - commodity_from_tech = self._get_commodity_for_tech(fixed_tech) - if self.commodity in commodity_from_tech: - demand = self._subtract_fixed(fixed_tech, demand, self.commodity, inputs) - - # 2. Flexible techs: full production - for flexible_tech in self.flexible_techs: - commodity_from_tech = self._get_commodity_for_tech(flexible_tech) - if self.commodity in commodity_from_tech: - demand = self._subtract_flexible( - flexible_tech, demand, self.commodity, inputs, outputs - ) - - # 3. Storage dispatch - # number of storage components that produce the demanded commodity - n_storage = len( - [s for s in self.storage_techs if self.commodity in self._get_commodity_for_tech(s)] - ) - for storage_tech in self.storage_techs: - commodity_from_tech = self._get_commodity_for_tech(storage_tech) - if self.commodity in commodity_from_tech: - demand = self._dispatch_storage( - storage_tech, demand / n_storage, self.commodity, inputs, outputs - ) - - # 4. Merit-order dispatch: cheapest dispatchable first - remaining = np.maximum(demand, 0.0) - - marginal_costs = self._compute_marginal_costs(inputs) - - # Merit order: sort by mean marginal cost (cheapest first) - mean_costs = np.array([mc.mean() for mc in marginal_costs]) - dispatch_order = np.argsort(mean_costs) - - # Initialize all dispatchable set-point outputs to zero - for set_point_name in self.dispatchable_set_point_names: - outputs[set_point_name] = np.zeros(self.n_timesteps) - - # Dispatch in merit order - for idx in dispatch_order: - set_point_name = self.dispatchable_set_point_names[idx] - rated_name = self.dispatchable_rated_names[idx] - rated = inputs[rated_name] + self._run_dispatch(inputs, outputs) - dispatch = np.minimum(remaining, rated) - outputs[set_point_name] = dispatch - remaining -= dispatch + def _dispatch_dispatchables(self, commodity, remaining_demand, inputs, outputs): + """Merit-order dispatch: cheapest dispatchables first, up to rated capacity.""" + return self._merit_order_dispatch(commodity, remaining_demand, inputs, outputs) diff --git a/h2integrate/control/control_strategies/system_level/demand_following_control.py b/h2integrate/control/control_strategies/system_level/demand_following_control.py index bb83712d0..a61d604a2 100644 --- a/h2integrate/control/control_strategies/system_level/demand_following_control.py +++ b/h2integrate/control/control_strategies/system_level/demand_following_control.py @@ -1,5 +1,3 @@ -import numpy as np - from h2integrate.control.control_strategies.system_level.system_level_control_base import ( SystemLevelControlBase, ) @@ -29,57 +27,14 @@ class DemandFollowingControl(SystemLevelControlBase): storage. The remaining demand (floored at zero) is split **evenly** across all dispatchable techs that produce the demanded commodity (each receives ``remaining_demand / n_dispatchable``). + + For heterogeneous-commodity systems, this four-step order is applied at + each commodity level: after the demand commodity is dispatched, the demand + placed on upstream converters (for example a hydrogen electrolyzer feeding + an ammonia synloop) is translated into demand for their input commodity via + the configured conversion ratios and the same dispatch is repeated. See + ``SystemLevelControlBase._run_dispatch``. """ def compute(self, inputs, outputs): - commodity = self.commodity - demand = inputs[self.demand_input_name].copy() - - # 1. Fixed techs: always produce, subtract from demand - for fixed_tech in self.fixed_techs: - commodity_from_tech = self._get_commodity_for_tech(fixed_tech) - for tech_commodity in commodity_from_tech: - if tech_commodity == commodity: - demand = self._subtract_fixed(fixed_tech, demand, commodity, inputs) - - # 2. Flexible techs: operate at full production - for flexible_tech in self.flexible_techs: - commodity_from_tech = self._get_commodity_for_tech(flexible_tech) - for tech_commodity in commodity_from_tech: - if tech_commodity == commodity: - demand = self._subtract_flexible( - flexible_tech, demand, commodity, inputs, outputs - ) - else: - if f"{flexible_tech}_rated_{tech_commodity}_production" in inputs: - # set the per-tech set-point as the rated production - outputs[f"{flexible_tech}_{tech_commodity}_set_point"] = inputs[ - f"{flexible_tech}_rated_{tech_commodity}_production" - ] * np.ones(self.n_timesteps) - - # 3. Storage dispatch - # number of storage components that produce the demanded commodity - n_storage = len( - [s for s in self.storage_techs if commodity in self._get_commodity_for_tech(s)] - ) - for storage_tech in self.storage_techs: - commodity_from_tech = self._get_commodity_for_tech(storage_tech) - if commodity in commodity_from_tech: - demand = self._dispatch_storage( - storage_tech, demand / n_storage, commodity, inputs, outputs - ) - - # 4. Dispatchable techs - remaining_demand = np.maximum(demand, 0.0) - - # calculate the number of dispatchable technologies that - # produce the demanded commodity - n_dispatchable = len( - [s for s in self.dispatchable_techs if commodity in self._get_commodity_for_tech(s)] - ) - for dispatchable_tech in self.dispatchable_techs: - commodity_from_tech = self._get_commodity_for_tech(dispatchable_tech) - if commodity in commodity_from_tech: - outputs[f"{dispatchable_tech}_{commodity}_set_point"] = ( - remaining_demand / n_dispatchable - ) + self._run_dispatch(inputs, outputs) diff --git a/h2integrate/control/control_strategies/system_level/profit_maximization_control.py b/h2integrate/control/control_strategies/system_level/profit_maximization_control.py index 9a8eb0a42..7e492d8a5 100644 --- a/h2integrate/control/control_strategies/system_level/profit_maximization_control.py +++ b/h2integrate/control/control_strategies/system_level/profit_maximization_control.py @@ -1,4 +1,3 @@ -import numpy as np from attrs import field, define from h2integrate.core.utilities import BaseConfig @@ -74,7 +73,7 @@ def _resolve_sell_price(self, config): price = group.get("model_inputs", {}).get("commodity_sell_price", None) if price is None: raise ValueError( - f"Finance group '{raw}' does not contain " f"model_inputs.commodity_sell_price." + f"Finance group '{raw}' does not contain model_inputs.commodity_sell_price." ) return price return raw @@ -101,57 +100,20 @@ def setup(self): self._setup_marginal_costs() def compute(self, inputs, outputs): - demand = inputs[self.demand_input_name].copy() - sell_price = inputs["commodity_sell_price"] # shape (n_timesteps,) - - # 1. Fixed techs: always produce, subtract from demand - for fixed_tech in self.fixed_techs: - commodity_from_tech = self._get_commodity_for_tech(fixed_tech) - if self.commodity in commodity_from_tech: - demand = self._subtract_fixed(fixed_tech, demand, self.commodity, inputs) - - # 2. Flexible techs: full production (always profitable) - for flexible_tech in self.flexible_techs: - commodity_from_tech = self._get_commodity_for_tech(flexible_tech) - if self.commodity in commodity_from_tech: - demand = self._subtract_flexible( - flexible_tech, demand, self.commodity, inputs, outputs - ) - - # 3. Storage dispatch - # number of storage components that produce the demanded commodity - n_storage = len( - [s for s in self.storage_techs if self.commodity in self._get_commodity_for_tech(s)] - ) - for storage_tech in self.storage_techs: - commodity_from_tech = self._get_commodity_for_tech(storage_tech) - if self.commodity in commodity_from_tech: - demand = self._dispatch_storage( - storage_tech, demand / n_storage, self.commodity, inputs, outputs - ) - - # 4. Profit-driven merit-order dispatch - remaining = np.maximum(demand, 0.0) - - marginal_costs = self._compute_marginal_costs(inputs) + self._run_dispatch(inputs, outputs) - # Merit order: sort by mean marginal cost (cheapest first) - mean_costs = np.array([mc.mean() for mc in marginal_costs]) - dispatch_order = np.argsort(mean_costs) + def _dispatch_dispatchables(self, commodity, remaining_demand, inputs, outputs): + """Merit-order dispatch, profit-gated for the primary demand commodity. - # Initialize all dispatchable set-point outputs to zero - for set_point_name in self.dispatchable_set_point_names: - outputs[set_point_name] = np.zeros(self.n_timesteps) - - # Dispatch only where profitable (element-wise comparison) - for idx in dispatch_order: - mc = marginal_costs[idx] # per-timestep array - profitable = mc < sell_price # boolean mask per timestep - - set_point_name = self.dispatchable_set_point_names[idx] - rated_name = self.dispatchable_rated_names[idx] - rated = inputs[rated_name] - - dispatch = np.where(profitable, np.minimum(remaining, rated), 0.0) - outputs[set_point_name] = dispatch - remaining -= dispatch + The demand commodity is dispatched only at timesteps where a + technology's marginal cost is below the sell price. Upstream (derived) + commodities are dispatched at minimum cost to supply the profitable + primary production, since the profitability decision has already been + made at the demand-commodity level. + """ + if commodity == self.commodity: + sell_price = inputs["commodity_sell_price"] + return self._merit_order_dispatch( + commodity, remaining_demand, inputs, outputs, sell_price=sell_price + ) + return self._merit_order_dispatch(commodity, remaining_demand, inputs, outputs) diff --git a/h2integrate/control/control_strategies/system_level/system_level_control_base.py b/h2integrate/control/control_strategies/system_level/system_level_control_base.py index 362b551ed..7e3d1e2a5 100644 --- a/h2integrate/control/control_strategies/system_level/system_level_control_base.py +++ b/h2integrate/control/control_strategies/system_level/system_level_control_base.py @@ -1,3 +1,6 @@ +import warnings +from collections import defaultdict + import numpy as np import networkx as nx import openmdao.api as om @@ -202,6 +205,12 @@ def setup(self): self._setup_tech_category("storage", self.storage_techs) self._setup_feedstock_category(self.feedstock_comps) + # Detect commodity converters and load their static conversion ratios + # (Phase 1: static ratios from tech_config). Enables backward demand + # propagation across converter boundaries for heterogeneous-commodity + # systems. + self._build_conversion_ratios() + def _setup_commodity( self, tech_name, @@ -588,6 +597,313 @@ def _get_commodity_for_tech(self, tech_name): return tech_commodities + # ------------------------------------------------------------------ + # Heterogeneous-commodity dispatch (backward demand propagation) + # ------------------------------------------------------------------ + + def _build_conversion_ratios(self): + """Detect commodity converters and read their static conversion ratios. + + A converter is a controller-managed technology that produces a commodity + it does not itself consume (for example an electrolyzer: + electricity -> hydrogen, or an ammonia synloop: hydrogen -> ammonia). + Consumed commodities are read directly from each technology's incoming + graph edges and produced commodities from ``techs_to_commodities``; this + direct-edge definition is robust for converter chains (A -> B -> C, + where B and C both convert). + + Populates three attributes used by ``_run_dispatch`` to translate demand + for one commodity into demand for an upstream (input) commodity: + + - ``self._converters``: set of ``(in_commodity, tech_name, + out_commodity)`` tuples (empty for single-commodity systems). + - ``self.conversion_ratios``: mapping of ``(tech_name, in_commodity, + out_commodity)`` to a float ratio read from the tech config at + ``technologies..model_inputs.control_parameters. + conversion_ratios._per_``. + - ``self._missing_ratio_warned``: set used to emit the "missing ratio" + warning at most once per converter. + + Ratios are interpreted such that ``input_rate = output_rate * ratio`` + in each commodity's rate units (for example ``51 kWh/kg`` translates a + hydrogen production rate in ``kg/h`` into an electricity demand in + ``kW``). Converters whose input commodity has no controller-managed + producer (e.g. a feedstock-supplied nitrogen stream) do not + require a ratio. + """ + + def _as_list(commodity): + if commodity is None: + return [] + if isinstance(commodity, str): + return [commodity] + return list(commodity) + + # Detect converters: techs that produce a commodity they do not consume + self._converters = set() + for tech_name in self.input_techs: + produced_commodities = set(self._get_commodity_for_tech(tech_name)) + if not produced_commodities: + continue + + consumed_commodities = set() + for _src, _dst, edge_commodity in self.technology_graph.in_edges( + tech_name, data="commodity" + ): + consumed_commodities.update(_as_list(edge_commodity)) + + for out_commodity in produced_commodities - consumed_commodities: + for in_commodity in consumed_commodities - produced_commodities: + self._converters.add((in_commodity, tech_name, out_commodity)) + + # Read each converter's static input-per-output ratio from the tech config + self.conversion_ratios = {} + self._missing_ratio_warned = set() + technologies = self.options["tech_config"].get("technologies", {}) + for in_commodity, tech_name, out_commodity in self._converters: + control_params = ( + technologies.get(tech_name, {}) + .get("model_inputs", {}) + .get("control_parameters", {}) + ) + ratios = control_params.get("conversion_ratios", {}) + key = f"{in_commodity}_per_{out_commodity}" + if key in ratios: + self.conversion_ratios[(tech_name, in_commodity, out_commodity)] = float( + ratios[key] + ) + + def _run_dispatch(self, inputs, outputs): + """Dispatch every commodity level required to meet the demand profile. + + The demand commodity is dispatched first. For each converter that + produces a just-dispatched commodity, the converter's committed output + set-point is translated into demand for its input commodity (via its + static conversion ratio) and accumulated. Commodities are processed in + topological order of this demand-flow graph so that all contributions + to an input commodity are gathered before it is dispatched. + + Subclasses implement the strategy-specific dispatchable step by + overriding ``_dispatch_dispatchables``; the fixed, flexible, and + storage steps are shared across all strategies. + """ + # Flexible techs can only curtail, so they always run at their rated + # production for every commodity they produce. Commanding them here + # guarantees they are set even when their commodity is not explicitly + # dispatched below (commodity-level dispatch may re-issue the value). + for tech_name in self.flexible_techs: + for commodity in self._get_commodity_for_tech(tech_name): + rated_name = f"{tech_name}_rated_{commodity}_production" + set_point_name = f"{tech_name}_{commodity}_set_point" + if rated_name in inputs and set_point_name in outputs: + outputs[set_point_name] = inputs[rated_name] * np.ones(self.n_timesteps) + + converters = getattr(self, "_converters", None) or set() + + # Single-commodity system: dispatch the demand commodity and return + if not converters: + self._dispatch_commodity( + self.commodity, inputs[self.demand_input_name].copy(), inputs, outputs + ) + return + + # Build the demand-flow graph (out_commodity -> in_commodity) and group + # converters by the commodity they output + demand_flow = nx.DiGraph() + demand_flow.add_node(self.commodity) + converters_by_output = defaultdict(list) + for in_commodity, tech_name, out_commodity in converters: + demand_flow.add_edge(out_commodity, in_commodity) + converters_by_output[out_commodity].append((in_commodity, tech_name, out_commodity)) + + try: + commodity_order = list(nx.topological_sort(demand_flow)) + except nx.NetworkXUnfeasible: + # Commodity cycle (unusual): process the demand commodity first, + # then the remaining commodities in arbitrary order + commodity_order = [self.commodity] + [ + c for c in demand_flow.nodes if c != self.commodity + ] + + derived_demand = {self.commodity: inputs[self.demand_input_name].copy()} + + for commodity in commodity_order: + demand = derived_demand.get(commodity) + if demand is None: + continue + + self._dispatch_commodity(commodity, demand.copy(), inputs, outputs) + + # Translate committed converter output into upstream input demand + for in_commodity, tech_name, out_commodity in converters_by_output.get(commodity, []): + self._accumulate_derived_demand( + derived_demand, tech_name, in_commodity, out_commodity, inputs, outputs + ) + + def _dispatch_commodity(self, commodity, demand, inputs, outputs): + """Dispatch all technologies producing ``commodity`` to meet ``demand``. + + Applies the shared four-step priority order (fixed, flexible, storage, + dispatchable) for a single commodity. Only technologies that produce + ``commodity`` participate. The dispatchable step is delegated to + ``_dispatch_dispatchables`` so cost-aware subclasses can override it. + + Args: + commodity (str): Commodity to dispatch (demand or a derived input). + demand (np.ndarray): Demand profile for ``commodity`` (may be mutated). + inputs: OpenMDAO inputs. + outputs: OpenMDAO outputs. + + Returns: + np.ndarray: Remaining unmet demand after dispatchables. + """ + # 1. Fixed techs: always produce, subtract from demand + for fixed_tech in self.fixed_techs: + if commodity in self._get_commodity_for_tech(fixed_tech): + demand = self._subtract_fixed(fixed_tech, demand, commodity, inputs) + + # 2. Flexible techs: run at rated production, subtract from demand + for flexible_tech in self.flexible_techs: + if commodity in self._get_commodity_for_tech(flexible_tech): + updated = self._subtract_flexible(flexible_tech, demand, commodity, inputs, outputs) + if updated is not None: + demand = updated + + # 3. Storage techs: split residual demand evenly and dispatch + n_storage = len( + [s for s in self.storage_techs if commodity in self._get_commodity_for_tech(s)] + ) + for storage_tech in self.storage_techs: + if commodity in self._get_commodity_for_tech(storage_tech): + updated = self._dispatch_storage( + storage_tech, demand / n_storage, commodity, inputs, outputs + ) + if updated is not None: + demand = updated + + # 4. Dispatchable techs: strategy-specific (subclass hook) + remaining = np.maximum(demand, 0.0) + return self._dispatch_dispatchables(commodity, remaining, inputs, outputs) + + def _dispatch_dispatchables(self, commodity, remaining_demand, inputs, outputs): + """Split remaining demand evenly across dispatchables producing ``commodity``. + + This is the default (cost-agnostic) dispatchable step used by + ``DemandFollowingControl``. Cost-aware controllers override this method + to apply merit-order or profit-aware dispatch. + + Returns: + np.ndarray: Remaining unmet demand (zeros under the even-split + assumption that dispatchables absorb their full share). + """ + dispatchables = [ + t for t in self.dispatchable_techs if commodity in self._get_commodity_for_tech(t) + ] + n_dispatchable = len(dispatchables) + if n_dispatchable == 0: + return remaining_demand + + for tech_name in dispatchables: + outputs[f"{tech_name}_{commodity}_set_point"] = remaining_demand / n_dispatchable + + return np.zeros(self.n_timesteps) + + def _merit_order_dispatch(self, commodity, remaining_demand, inputs, outputs, sell_price=None): + """Dispatch dispatchables producing ``commodity`` in ascending marginal-cost order. + + Each technology is dispatched up to its rated production until demand + is met. When ``sell_price`` is provided, a technology is only + dispatched at timesteps where its marginal cost is below the sell price + (profit gating); demand may then go unmet. + + Marginal costs come from ``_compute_marginal_costs`` and are aligned + with ``self.dispatchable_techs``. Requires ``_setup_marginal_costs`` to + have been called in the subclass ``setup``. + + Returns: + np.ndarray: Remaining unmet demand after dispatch. + """ + dispatchables = [ + t for t in self.dispatchable_techs if commodity in self._get_commodity_for_tech(t) + ] + + # Initialize set-points for these dispatchables to zero + for tech_name in dispatchables: + outputs[f"{tech_name}_{commodity}_set_point"] = np.zeros(self.n_timesteps) + + if not dispatchables: + return remaining_demand + + marginal_cost_by_tech = dict( + zip(self.dispatchable_techs, self._compute_marginal_costs(inputs)) + ) + + # Merit order: cheapest mean marginal cost first + dispatch_order = sorted(dispatchables, key=lambda t: marginal_cost_by_tech[t].mean()) + + remaining = np.array(remaining_demand, dtype=float) + for tech_name in dispatch_order: + rated = inputs[f"{tech_name}_rated_{commodity}_production"] + if sell_price is not None: + profitable = marginal_cost_by_tech[tech_name] < sell_price + dispatch = np.where(profitable, np.minimum(remaining, rated), 0.0) + else: + dispatch = np.minimum(remaining, rated) + outputs[f"{tech_name}_{commodity}_set_point"] = dispatch + remaining = remaining - dispatch + + return remaining + + def _accumulate_derived_demand( + self, derived_demand, tech_name, in_commodity, out_commodity, inputs, outputs + ): + """Add a converter's induced input-commodity demand to ``derived_demand``. + + The converter's committed output (its ``{tech}_{out}_set_point``) is + multiplied by the input-per-output conversion ratio to obtain the demand + it places on its input commodity. Converters whose input commodity has + no controller-managed producer (for example a feedstock-supplied stream) + are skipped. + + Backward propagation is opt-in: if no conversion ratio is configured for + a converter whose input commodity does have controllable producers, the + propagation is skipped (upstream techs keep their default dispatch, the + legacy behavior) and a one-time warning is emitted so the missing ratio + is discoverable when heterogeneous-commodity control is intended. + """ + has_producers = any( + in_commodity in self._get_commodity_for_tech(t) for t in self.input_techs + ) + if not has_producers: + return + + key = (tech_name, in_commodity, out_commodity) + if key not in self.conversion_ratios: + if key not in self._missing_ratio_warned: + self._missing_ratio_warned.add(key) + warnings.warn( + f"No conversion ratio defined for converter '{tech_name}' " + f"('{in_commodity}' -> '{out_commodity}'); system-level control will " + f"not translate '{out_commodity}' demand into '{in_commodity}' demand, " + f"and upstream '{in_commodity}' technologies keep their default " + f"dispatch. Define technologies.{tech_name}.model_inputs." + f"control_parameters.conversion_ratios.{in_commodity}_per_{out_commodity} " + f"in the tech config to enable heterogeneous-commodity control.", + stacklevel=2, + ) + return + + # ``key`` is guaranteed present here (checked above); apply the static + # input-per-output ratio to the converter's committed output + ratio = np.full(self.n_timesteps, self.conversion_ratios[key]) + set_point = np.asarray(outputs[f"{tech_name}_{out_commodity}_set_point"], dtype=float) + contribution = np.maximum(set_point, 0.0) * ratio + + if in_commodity in derived_demand: + derived_demand[in_commodity] = derived_demand[in_commodity] + contribution + else: + derived_demand[in_commodity] = contribution + # ------------------------------------------------------------------ # Marginal-cost helpers for cost-aware controllers # ------------------------------------------------------------------ diff --git a/h2integrate/control/control_strategies/system_level/test/test_slc_controllers.py b/h2integrate/control/control_strategies/system_level/test/test_slc_controllers.py index cfd46a14c..24c58a625 100644 --- a/h2integrate/control/control_strategies/system_level/test/test_slc_controllers.py +++ b/h2integrate/control/control_strategies/system_level/test/test_slc_controllers.py @@ -722,3 +722,313 @@ def test_feedstock_no_feedstock_raises(self): with pytest.raises(ValueError, match="at least one feedstock"): _build_problem(CostMinimizationControl, plant_config, slc_topology, demand=50000) + + +# --------------------------------------------------------------------------- +# Heterogeneous-commodity dispatch (backward demand propagation) +# --------------------------------------------------------------------------- +def _tech_config_with_ratios(ratios_by_tech): + """Build a tech_config carrying static conversion ratios. + + Args: + ratios_by_tech (dict): Mapping of ``tech_name`` to a dict of + ``{"_per_": ratio}`` entries. + + Returns: + dict: A ``tech_config`` with the nested ``model_inputs. + control_parameters.conversion_ratios`` structure the base class reads. + """ + return { + "technologies": { + tech: {"model_inputs": {"control_parameters": {"conversion_ratios": ratios}}} + for tech, ratios in ratios_by_tech.items() + } + } + + +def _build_hetero_problem( + slc_cls, + plant_config, + slc_topology, + tech_config, + demand, + upstream_out=None, + commodity_units=None, +): + """Build an SLC problem wiring non-demand commodity outputs via IVCs. + + Every ``(tech, commodity)`` output whose commodity differs from the demand + commodity is fed by an ``IndepVarComp`` and connected into the controller, + mirroring the real plant connections so ``units_by_conn`` inputs resolve. + + Args: + slc_cls: Controller class to instantiate. + plant_config (dict): Plant config. + slc_topology (dict): SLC topology. + tech_config (dict): Tech config (carries conversion ratios). + demand (float | array): Demand-commodity demand profile. + upstream_out (dict): Optional ``{(tech, commodity): value}`` outputs for + the wired IVCs (defaults to zero). + commodity_units (dict): Optional ``{commodity: units}`` for the IVCs. + + Returns: + om.Problem: A setup problem with the demand value applied. + """ + upstream_out = upstream_out or {} + commodity_units = commodity_units or {} + n_timesteps = plant_config["plant"]["simulation"]["n_timesteps"] + demand_commodity = slc_topology["demand_commodity"] + + prob = om.Problem() + + ivc_connections = [] + for i, (tech, commodity) in enumerate(sorted(slc_topology["tech_to_commodity"])): + if commodity == demand_commodity: + continue + val = upstream_out.get((tech, commodity), 0.0) + val = np.full(n_timesteps, val) if np.isscalar(val) else np.asarray(val, dtype=float) + ivc = prob.model.add_subsystem(f"src{i}", om.IndepVarComp()) + ivc.add_output(f"{tech}_{commodity}_out", val=val, units=commodity_units.get(commodity)) + ivc_connections.append((f"src{i}.{tech}_{commodity}_out", f"slc.{tech}_{commodity}_out")) + + prob.model.add_subsystem( + "slc", + slc_cls( + driver_config={}, + plant_config=plant_config, + tech_config=tech_config, + slc_topology=slc_topology, + ), + ) + for src, dst in ivc_connections: + prob.model.connect(src, dst) + + prob.setup() + prob.set_val(f"slc.{demand_commodity}_demand", demand) + return prob + + +@pytest.mark.unit +class TestHeterogeneousCommodityControl: + """Backward demand propagation across commodity converters.""" + + def test_detect_converters_chain(self): + """Converters are detected across a grid -> electrolyzer -> synloop chain.""" + tech_connections = [ + ["grid", "electrolyzer", "electricity", "cable"], + ["electrolyzer", "synloop", "hydrogen", "pipe"], + ["synloop", "demand", "ammonia", "pipe"], + ] + plant_config = _build_plant_config(tech_connections) + tech_graph = _build_technology_graph(tech_connections) + classifiers = _build_tech_control_classifiers( + dispatchable=["grid", "electrolyzer", "synloop"] + ) + slc_topology = _build_slc_topology( + tech_graph, classifiers, demand_commodity="ammonia", demand_commodity_rate_units="kg/h" + ) + tech_config = _tech_config_with_ratios( + { + "electrolyzer": {"electricity_per_hydrogen": 51.0}, + "synloop": {"hydrogen_per_ammonia": 0.18}, + } + ) + prob = _build_hetero_problem( + DemandFollowingControl, + plant_config, + slc_topology, + tech_config, + demand=100.0, + commodity_units={"electricity": "kW", "hydrogen": "kg/h"}, + ) + assert prob.model.slc._converters == { + ("electricity", "electrolyzer", "hydrogen"), + ("hydrogen", "synloop", "ammonia"), + } + + def test_single_converter_static_propagation(self): + """A single converter translates hydrogen demand into electricity demand.""" + tech_connections = [ + ["grid", "electrolyzer", "electricity", "cable"], + ["electrolyzer", "demand", "hydrogen", "pipe"], + ] + plant_config = _build_plant_config(tech_connections) + tech_graph = _build_technology_graph(tech_connections) + classifiers = _build_tech_control_classifiers(dispatchable=["grid", "electrolyzer"]) + slc_topology = _build_slc_topology( + tech_graph, + classifiers, + demand_commodity="hydrogen", + demand_commodity_rate_units="kg/h", + ) + tech_config = _tech_config_with_ratios({"electrolyzer": {"electricity_per_hydrogen": 51.0}}) + prob = _build_hetero_problem( + DemandFollowingControl, + plant_config, + slc_topology, + tech_config, + demand=100.0, + commodity_units={"electricity": "kW"}, + ) + prob.set_val("slc.grid_rated_electricity_production", 1e9) + prob.set_val("slc.electrolyzer_rated_hydrogen_production", 1e6) + prob.run_model() + + np.testing.assert_allclose(prob.get_val("slc.electrolyzer_hydrogen_set_point"), 100.0) + np.testing.assert_allclose(prob.get_val("slc.grid_electricity_set_point"), 100.0 * 51.0) + + def test_chained_converter_static_propagation(self): + """Demand propagates through a two-converter chain to the electricity source.""" + tech_connections = [ + ["grid", "electrolyzer", "electricity", "cable"], + ["electrolyzer", "synloop", "hydrogen", "pipe"], + ["synloop", "demand", "ammonia", "pipe"], + ] + plant_config = _build_plant_config(tech_connections) + tech_graph = _build_technology_graph(tech_connections) + classifiers = _build_tech_control_classifiers( + dispatchable=["grid", "electrolyzer", "synloop"] + ) + slc_topology = _build_slc_topology( + tech_graph, + classifiers, + demand_commodity="ammonia", + demand_commodity_rate_units="kg/h", + ) + tech_config = _tech_config_with_ratios( + { + "electrolyzer": {"electricity_per_hydrogen": 51.0}, + "synloop": {"hydrogen_per_ammonia": 0.18}, + } + ) + prob = _build_hetero_problem( + DemandFollowingControl, + plant_config, + slc_topology, + tech_config, + demand=100.0, + commodity_units={"electricity": "kW", "hydrogen": "kg/h"}, + ) + prob.set_val("slc.grid_rated_electricity_production", 1e9) + prob.set_val("slc.electrolyzer_rated_hydrogen_production", 1e6) + prob.set_val("slc.synloop_rated_ammonia_production", 1e6) + prob.run_model() + + np.testing.assert_allclose(prob.get_val("slc.synloop_ammonia_set_point"), 100.0) + np.testing.assert_allclose( + prob.get_val("slc.electrolyzer_hydrogen_set_point"), 100.0 * 0.18 + ) + np.testing.assert_allclose( + prob.get_val("slc.grid_electricity_set_point"), 100.0 * 0.18 * 51.0 + ) + + def test_derived_demand_reuses_flexible_and_dispatchable(self): + """Derived electricity demand flows through the shared flexible/dispatchable steps.""" + tech_connections = [ + ["wind", "electrolyzer", "electricity", "cable"], + ["grid", "electrolyzer", "electricity", "cable"], + ["electrolyzer", "demand", "hydrogen", "pipe"], + ] + plant_config = _build_plant_config(tech_connections) + tech_graph = _build_technology_graph(tech_connections) + classifiers = _build_tech_control_classifiers( + flexible=["wind"], dispatchable=["grid", "electrolyzer"] + ) + slc_topology = _build_slc_topology( + tech_graph, + classifiers, + demand_commodity="hydrogen", + demand_commodity_rate_units="kg/h", + ) + tech_config = _tech_config_with_ratios({"electrolyzer": {"electricity_per_hydrogen": 50.0}}) + prob = _build_hetero_problem( + DemandFollowingControl, + plant_config, + slc_topology, + tech_config, + demand=100.0, + upstream_out={("wind", "electricity"): 1000.0}, + commodity_units={"electricity": "kW"}, + ) + prob.set_val("slc.wind_rated_electricity_production", 2000.0) + prob.set_val("slc.grid_rated_electricity_production", 1e9) + prob.set_val("slc.electrolyzer_rated_hydrogen_production", 1e6) + prob.run_model() + + # Derived electricity demand = 100 kg/h * 50 kWh/kg = 5000 kW + # Flexible wind runs at rated (2000 kW), curtailing 5000 - 1000 = 4000 kW of demand + np.testing.assert_allclose(prob.get_val("slc.electrolyzer_hydrogen_set_point"), 100.0) + np.testing.assert_allclose(prob.get_val("slc.wind_electricity_set_point"), 2000.0) + np.testing.assert_allclose(prob.get_val("slc.grid_electricity_set_point"), 4000.0) + + def test_missing_ratio_warns_and_keeps_legacy_dispatch(self): + """A converter without a ratio warns once and skips propagation.""" + tech_connections = [ + ["grid", "electrolyzer", "electricity", "cable"], + ["electrolyzer", "demand", "hydrogen", "pipe"], + ] + plant_config = _build_plant_config(tech_connections) + tech_graph = _build_technology_graph(tech_connections) + classifiers = _build_tech_control_classifiers(dispatchable=["grid", "electrolyzer"]) + slc_topology = _build_slc_topology( + tech_graph, + classifiers, + demand_commodity="hydrogen", + demand_commodity_rate_units="kg/h", + ) + prob = _build_hetero_problem( + DemandFollowingControl, + plant_config, + slc_topology, + tech_config={}, + demand=100.0, + commodity_units={"electricity": "kW"}, + ) + prob.set_val("slc.grid_rated_electricity_production", 1e9) + prob.set_val("slc.electrolyzer_rated_hydrogen_production", 1e6) + + with pytest.warns(UserWarning, match="No conversion ratio"): + prob.run_model() + + # Hydrogen demand is still met; electricity is not driven (keeps default set-point) + np.testing.assert_allclose(prob.get_val("slc.electrolyzer_hydrogen_set_point"), 100.0) + np.testing.assert_allclose(prob.get_val("slc.grid_electricity_set_point"), 1.0) + + def test_cost_min_merit_order_at_derived_level(self): + """Merit order applies to derived electricity demand under cost minimization.""" + tech_connections = [ + ["cheap", "electrolyzer", "electricity", "cable"], + ["expensive", "electrolyzer", "electricity", "cable"], + ["electrolyzer", "demand", "hydrogen", "pipe"], + ] + plant_config = _build_plant_config( + tech_connections, cost_per_tech={"cheap": 0.03, "expensive": 0.08} + ) + tech_graph = _build_technology_graph(tech_connections) + classifiers = _build_tech_control_classifiers( + dispatchable=["cheap", "expensive", "electrolyzer"] + ) + slc_topology = _build_slc_topology( + tech_graph, + classifiers, + demand_commodity="hydrogen", + demand_commodity_rate_units="kg/h", + ) + tech_config = _tech_config_with_ratios({"electrolyzer": {"electricity_per_hydrogen": 50.0}}) + prob = _build_hetero_problem( + CostMinimizationControl, + plant_config, + slc_topology, + tech_config, + demand=100.0, + commodity_units={"electricity": "kW"}, + ) + prob.set_val("slc.cheap_rated_electricity_production", 3000.0) + prob.set_val("slc.expensive_rated_electricity_production", 5000.0) + prob.set_val("slc.electrolyzer_rated_hydrogen_production", 1e6) + prob.run_model() + + # Derived electricity demand = 5000 kW; cheapest tech fills first + np.testing.assert_allclose(prob.get_val("slc.electrolyzer_hydrogen_set_point"), 100.0) + np.testing.assert_allclose(prob.get_val("slc.cheap_electricity_set_point"), 3000.0) + np.testing.assert_allclose(prob.get_val("slc.expensive_electricity_set_point"), 2000.0) From 8fd0c38d715346f69827ac63555d428fe66d6f83 Mon Sep 17 00:00:00 2001 From: John Jasa Date: Thu, 13 Aug 2026 23:10:47 -0600 Subject: [PATCH 2/5] Extending to timeseries for hetereo commodities --- CHANGELOG.md | 1 + .../system_level/system_level_control_base.py | 235 ++++++++++++++---- .../system_level/test/test_slc_controllers.py | 106 ++++++++ .../system_level/test/test_slc_examples.py | 2 +- h2integrate/core/h2integrate_model.py | 42 ++++ 5 files changed, 331 insertions(+), 55 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 39fc54ee8..5338b2ec6 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,7 @@ - Add `PySAMWavePerformanceModel` and `WaveResource` to wrap PySAM MhkWave as an H2I performance model, replacing the HOPP wave module in example 09. [PR 825](https://github.com/NatLabRockies/H2Integrate/pull/825) - Replace HOPP with native H2I wind, solar, and battery models in example 11. Adds `percent_load_missed` and `curtailment_percent` outputs to `DemandComponentBase`, allows zero capacity in wind/solar/battery performance models. [PR 826](https://github.com/NatLabRockies/H2Integrate/pull/826) - Add heterogeneous-commodity system-level control that translates demand for one commodity into upstream set-points across converters using static per-technology conversion ratios defined in the tech config. [PR TBD](https://github.com/NatLabRockies/H2Integrate/pull/TBD) + - Extend the heterogeneous-commodity control to prefer measured conversion ratios computed from each converter's consumed and produced streams per timestep, falling back to the static ratio when a measurement is unavailable. [PR TBD](https://github.com/NatLabRockies/H2Integrate/pull/TBD) ## 0.9 [August 10, 2026] diff --git a/h2integrate/control/control_strategies/system_level/system_level_control_base.py b/h2integrate/control/control_strategies/system_level/system_level_control_base.py index 7e3d1e2a5..9dce96b74 100644 --- a/h2integrate/control/control_strategies/system_level/system_level_control_base.py +++ b/h2integrate/control/control_strategies/system_level/system_level_control_base.py @@ -84,6 +84,58 @@ def _get_buy_price_default_and_shape(tech_config, tech_name, n_timesteps, plant_ return 0.0, n_timesteps +def detect_commodity_converters(technology_graph, input_techs, produced_by_tech): + """Detect commodity converters from the technology graph. + + A converter is a controller-managed technology that produces a commodity it + does not itself consume (for example an electrolyzer: electricity -> + hydrogen, or an ammonia synloop: hydrogen -> ammonia). Consumed commodities + are read directly from each technology's incoming graph edges and produced + commodities from ``produced_by_tech``; this direct-edge definition is robust + for converter chains (A -> B -> C, where B and C both convert). + + This module-level helper is shared by the controller component (which calls + it with its own graph and classification) and ``H2IntegrateModel`` (which + calls it during classification so the same converter set drives the + consumption-signal connections). + + Args: + technology_graph (nx.DiGraph): Directed technology graph with a + ``commodity`` attribute on each edge. + input_techs (Iterable[str]): Controller-managed technologies that may + produce a commodity (fixed, flexible, dispatchable, storage). + produced_by_tech (Mapping[str, Iterable[str]]): Mapping of technology + name to the commodities it produces. + + Returns: + set[tuple[str, str, str]]: ``(in_commodity, tech_name, out_commodity)`` + tuples, one per detected conversion. + """ + + def _as_list(commodity): + if commodity is None: + return [] + if isinstance(commodity, str): + return [commodity] + return list(commodity) + + converters = set() + for tech_name in input_techs: + produced_commodities = set(produced_by_tech.get(tech_name, ())) + if not produced_commodities: + continue + + consumed_commodities = set() + for _src, _dst, edge_commodity in technology_graph.in_edges(tech_name, data="commodity"): + consumed_commodities.update(_as_list(edge_commodity)) + + for out_commodity in produced_commodities - consumed_commodities: + for in_commodity in consumed_commodities - produced_commodities: + converters.add((in_commodity, tech_name, out_commodity)) + + return converters + + class SystemLevelControlBase(om.ExplicitComponent): """Base class for system-level controllers. @@ -602,59 +654,50 @@ def _get_commodity_for_tech(self, tech_name): # ------------------------------------------------------------------ def _build_conversion_ratios(self): - """Detect commodity converters and read their static conversion ratios. + """Detect commodity converters and read their conversion ratios. A converter is a controller-managed technology that produces a commodity it does not itself consume (for example an electrolyzer: electricity -> hydrogen, or an ammonia synloop: hydrogen -> ammonia). - Consumed commodities are read directly from each technology's incoming - graph edges and produced commodities from ``techs_to_commodities``; this - direct-edge definition is robust for converter chains (A -> B -> C, - where B and C both convert). + Converters are taken from ``slc_topology["converters"]`` when + ``H2IntegrateModel`` provides them, and otherwise detected here with + ``detect_commodity_converters`` so the component remains usable + standalone (for example in unit tests). - Populates three attributes used by ``_run_dispatch`` to translate demand + Populates the attributes used by ``_run_dispatch`` to translate demand for one commodity into demand for an upstream (input) commodity: - ``self._converters``: set of ``(in_commodity, tech_name, out_commodity)`` tuples (empty for single-commodity systems). - ``self.conversion_ratios``: mapping of ``(tech_name, in_commodity, - out_commodity)`` to a float ratio read from the tech config at + out_commodity)`` to a float static ratio read from the tech config at ``technologies..model_inputs.control_parameters. conversion_ratios._per_``. + - ``self._converter_consumed_names``: mapping of the same key to the + ``{tech}_{in_commodity}_consumed`` input registered for the dynamic + (measured) ratio path. - ``self._missing_ratio_warned``: set used to emit the "missing ratio" warning at most once per converter. - Ratios are interpreted such that ``input_rate = output_rate * ratio`` - in each commodity's rate units (for example ``51 kWh/kg`` translates a - hydrogen production rate in ``kg/h`` into an electricity demand in - ``kW``). Converters whose input commodity has no controller-managed - producer (e.g. a feedstock-supplied nitrogen stream) do not - require a ratio. + Static ratios are interpreted such that ``input_rate = output_rate * + ratio`` in each commodity's rate units (for example ``51 kWh/kg`` + translates a hydrogen production rate in ``kg/h`` into an electricity + demand in ``kW``). When a converter reports its consumption via a + ``{in_commodity}_consumed`` output (wired by ``H2IntegrateModel``), the + ratio is measured per timestep as ``consumed / produced`` and the static + ratio is used only as the zero-output fallback. Converters whose input + commodity has no controller-managed producer (for example a + feedstock-supplied nitrogen stream) do not require a ratio. """ - - def _as_list(commodity): - if commodity is None: - return [] - if isinstance(commodity, str): - return [commodity] - return list(commodity) - - # Detect converters: techs that produce a commodity they do not consume - self._converters = set() - for tech_name in self.input_techs: - produced_commodities = set(self._get_commodity_for_tech(tech_name)) - if not produced_commodities: - continue - - consumed_commodities = set() - for _src, _dst, edge_commodity in self.technology_graph.in_edges( - tech_name, data="commodity" - ): - consumed_commodities.update(_as_list(edge_commodity)) - - for out_commodity in produced_commodities - consumed_commodities: - for in_commodity in consumed_commodities - produced_commodities: - self._converters.add((in_commodity, tech_name, out_commodity)) + converters = self.options["slc_topology"].get("converters") + if converters is None: + produced_by_tech = defaultdict(set) + for tech_name, commodity in self.techs_to_commodities: + produced_by_tech[tech_name].add(commodity) + converters = detect_commodity_converters( + self.technology_graph, self.input_techs, produced_by_tech + ) + self._converters = set(converters) # Read each converter's static input-per-output ratio from the tech config self.conversion_ratios = {} @@ -673,15 +716,50 @@ def _as_list(commodity): ratios[key] ) + # Register a measured-consumption input for every converter whose input + # commodity has a controller-managed producer. H2IntegrateModel wires + # the converter's ``{in_commodity}_consumed`` output to this input so + # the ratio can be measured per timestep. The NaN default marks the + # input as unconnected (dynamic ratio unavailable) when running the + # component standalone + self._converter_consumed_names = {} + for in_commodity, tech_name, out_commodity in self._converters: + has_producers = any( + in_commodity in self._get_commodity_for_tech(t) for t in self.input_techs + ) + if not has_producers: + continue + + if in_commodity in self.commodities_to_units: + unit_kwargs = {"units": self.commodities_to_units[in_commodity]} + elif in_commodity in self.commodities_to_ref_var: + unit_kwargs = { + "units": None, + "copy_units": self.commodities_to_ref_var[in_commodity], + } + else: + unit_kwargs = {"units": None} + + consumed_name = f"{tech_name}_{in_commodity}_consumed" + self.add_input( + consumed_name, + val=np.full(self.n_timesteps, np.nan), + shape=self.n_timesteps, + desc=f"Measured {in_commodity} consumed by converter {tech_name}", + **unit_kwargs, + ) + self._converter_consumed_names[(tech_name, in_commodity, out_commodity)] = consumed_name + def _run_dispatch(self, inputs, outputs): """Dispatch every commodity level required to meet the demand profile. The demand commodity is dispatched first. For each converter that produces a just-dispatched commodity, the converter's committed output set-point is translated into demand for its input commodity (via its - static conversion ratio) and accumulated. Commodities are processed in - topological order of this demand-flow graph so that all contributions - to an input commodity are gathered before it is dispatched. + conversion ratio, measured when available and otherwise static) and + accumulated. Commodities are processed in topological order of this + demand-flow graph so that all contributions to an input commodity are + gathered before it is dispatched. Subclasses implement the strategy-specific dispatchable step by overriding ``_dispatch_dispatchables``; the fixed, flexible, and @@ -865,11 +943,12 @@ def _accumulate_derived_demand( no controller-managed producer (for example a feedstock-supplied stream) are skipped. - Backward propagation is opt-in: if no conversion ratio is configured for - a converter whose input commodity does have controllable producers, the - propagation is skipped (upstream techs keep their default dispatch, the - legacy behavior) and a one-time warning is emitted so the missing ratio - is discoverable when heterogeneous-commodity control is intended. + Backward propagation is opt-in: if a converter whose input commodity has + controllable producers exposes neither a static ratio nor a connected + consumption signal, the propagation is skipped (upstream techs keep their + default dispatch, the legacy behavior) and a one-time warning is emitted + so the missing ratio is discoverable when heterogeneous-commodity control + is intended. """ has_producers = any( in_commodity in self._get_commodity_for_tech(t) for t in self.input_techs @@ -877,25 +956,24 @@ def _accumulate_derived_demand( if not has_producers: return - key = (tech_name, in_commodity, out_commodity) - if key not in self.conversion_ratios: + ratio = self._conversion_ratio(tech_name, in_commodity, out_commodity, inputs) + if ratio is None: + key = (tech_name, in_commodity, out_commodity) if key not in self._missing_ratio_warned: self._missing_ratio_warned.add(key) warnings.warn( - f"No conversion ratio defined for converter '{tech_name}' " + f"No conversion ratio available for converter '{tech_name}' " f"('{in_commodity}' -> '{out_commodity}'); system-level control will " f"not translate '{out_commodity}' demand into '{in_commodity}' demand, " f"and upstream '{in_commodity}' technologies keep their default " - f"dispatch. Define technologies.{tech_name}.model_inputs." - f"control_parameters.conversion_ratios.{in_commodity}_per_{out_commodity} " - f"in the tech config to enable heterogeneous-commodity control.", + f"dispatch. Connect the converter's '{in_commodity}_consumed' output or " + f"define technologies.{tech_name}.model_inputs.control_parameters." + f"conversion_ratios.{in_commodity}_per_{out_commodity} in the tech config " + f"to enable heterogeneous-commodity control.", stacklevel=2, ) return - # ``key`` is guaranteed present here (checked above); apply the static - # input-per-output ratio to the converter's committed output - ratio = np.full(self.n_timesteps, self.conversion_ratios[key]) set_point = np.asarray(outputs[f"{tech_name}_{out_commodity}_set_point"], dtype=float) contribution = np.maximum(set_point, 0.0) * ratio @@ -904,6 +982,55 @@ def _accumulate_derived_demand( else: derived_demand[in_commodity] = contribution + def _conversion_ratio(self, tech_name, in_commodity, out_commodity, inputs): + """Return a per-timestep input-per-output conversion ratio, or ``None``. + + Ratio precedence: + + 1. Dynamic (measured): where the converter reports its consumption via a + connected ``{tech}_{in_commodity}_consumed`` input, the ratio is + ``consumed / produced`` at every timestep with nonzero production. + This is the general, nonlinear-aware path; the plant solver resolves + the resulting feedback without an analytic Jacobian. + 2. Static (fallback): the constant ratio from the tech config, used at + zero-output timesteps and when no consumption is measured. When no + static ratio is configured, the mean measured ratio is used as the + zero-output fallback. + + Returns ``None`` when neither a static ratio nor a connected consumption + signal is available, signalling the caller to skip propagation. + """ + key = (tech_name, in_commodity, out_commodity) + static = self.conversion_ratios.get(key) + + consumed = None + consumed_name = self._converter_consumed_names.get(key) + if consumed_name is not None: + measured = np.asarray(inputs[consumed_name], dtype=float) + # A registered but unconnected input keeps its NaN sentinel; treat it + # as "no measurement" so the static ratio (or the warning) applies + if np.any(np.isfinite(measured)): + consumed = measured + + if static is None and consumed is None: + return None + + if consumed is None: + return np.full(self.n_timesteps, float(static)) + + produced = np.asarray(inputs[f"{tech_name}_{out_commodity}_out"], dtype=float) + valid = np.isfinite(consumed) & (produced != 0.0) + dynamic = np.divide(consumed, produced, out=np.zeros_like(produced), where=valid) + + if static is not None: + nominal = float(static) + elif np.any(valid): + nominal = float(dynamic[valid].mean()) + else: + nominal = 0.0 + + return np.where(valid, dynamic, nominal) + # ------------------------------------------------------------------ # Marginal-cost helpers for cost-aware controllers # ------------------------------------------------------------------ diff --git a/h2integrate/control/control_strategies/system_level/test/test_slc_controllers.py b/h2integrate/control/control_strategies/system_level/test/test_slc_controllers.py index 24c58a625..cf84fed9e 100644 --- a/h2integrate/control/control_strategies/system_level/test/test_slc_controllers.py +++ b/h2integrate/control/control_strategies/system_level/test/test_slc_controllers.py @@ -1,5 +1,7 @@ """Unit tests for system-level control base class and all controller strategies.""" +import warnings + import numpy as np import pytest import networkx as nx @@ -1032,3 +1034,107 @@ def test_cost_min_merit_order_at_derived_level(self): np.testing.assert_allclose(prob.get_val("slc.electrolyzer_hydrogen_set_point"), 100.0) np.testing.assert_allclose(prob.get_val("slc.cheap_electricity_set_point"), 3000.0) np.testing.assert_allclose(prob.get_val("slc.expensive_electricity_set_point"), 2000.0) + + def test_dynamic_ratio_overrides_static(self): + """Measured consumption drives the ratio, overriding the static tech-config value.""" + tech_connections = [ + ["grid", "electrolyzer", "electricity", "cable"], + ["electrolyzer", "demand", "hydrogen", "pipe"], + ] + plant_config = _build_plant_config(tech_connections) + tech_graph = _build_technology_graph(tech_connections) + classifiers = _build_tech_control_classifiers(dispatchable=["grid", "electrolyzer"]) + slc_topology = _build_slc_topology( + tech_graph, + classifiers, + demand_commodity="hydrogen", + demand_commodity_rate_units="kg/h", + ) + # Static ratio is 50, but the measured ratio (5100 / 100 = 51) should win. + tech_config = _tech_config_with_ratios({"electrolyzer": {"electricity_per_hydrogen": 50.0}}) + prob = _build_hetero_problem( + DemandFollowingControl, + plant_config, + slc_topology, + tech_config, + demand=100.0, + commodity_units={"electricity": "kW"}, + ) + prob.set_val("slc.grid_rated_electricity_production", 1e9) + prob.set_val("slc.electrolyzer_rated_hydrogen_production", 1e6) + prob.set_val("slc.electrolyzer_hydrogen_out", 100.0) + prob.set_val("slc.electrolyzer_electricity_consumed", 5100.0) + prob.run_model() + + np.testing.assert_allclose(prob.get_val("slc.grid_electricity_set_point"), 100.0 * 51.0) + + def test_dynamic_ratio_time_varying_with_zero_output_fallback(self): + """Per-timestep measured ratios apply; a zero-output timestep falls back to static.""" + tech_connections = [ + ["grid", "electrolyzer", "electricity", "cable"], + ["electrolyzer", "demand", "hydrogen", "pipe"], + ] + plant_config = _build_plant_config(tech_connections) + tech_graph = _build_technology_graph(tech_connections) + classifiers = _build_tech_control_classifiers(dispatchable=["grid", "electrolyzer"]) + slc_topology = _build_slc_topology( + tech_graph, + classifiers, + demand_commodity="hydrogen", + demand_commodity_rate_units="kg/h", + ) + tech_config = _tech_config_with_ratios({"electrolyzer": {"electricity_per_hydrogen": 50.0}}) + prob = _build_hetero_problem( + DemandFollowingControl, + plant_config, + slc_topology, + tech_config, + demand=100.0, + commodity_units={"electricity": "kW"}, + ) + prob.set_val("slc.grid_rated_electricity_production", 1e9) + prob.set_val("slc.electrolyzer_rated_hydrogen_production", 1e6) + # Third timestep produces no hydrogen, so its ratio falls back to the static 50. + prob.set_val("slc.electrolyzer_hydrogen_out", [100.0, 100.0, 0.0, 100.0]) + prob.set_val("slc.electrolyzer_electricity_consumed", [5100.0, 4000.0, 9999.0, 6000.0]) + prob.run_model() + + # ratio = [51, 40, 50 (fallback), 60]; derived electricity = 100 * ratio + np.testing.assert_allclose( + prob.get_val("slc.grid_electricity_set_point"), [5100.0, 4000.0, 5000.0, 6000.0] + ) + + def test_dynamic_ratio_without_static_does_not_warn(self): + """A connected consumption signal enables propagation with no static ratio or warning.""" + tech_connections = [ + ["grid", "electrolyzer", "electricity", "cable"], + ["electrolyzer", "demand", "hydrogen", "pipe"], + ] + plant_config = _build_plant_config(tech_connections) + tech_graph = _build_technology_graph(tech_connections) + classifiers = _build_tech_control_classifiers(dispatchable=["grid", "electrolyzer"]) + slc_topology = _build_slc_topology( + tech_graph, + classifiers, + demand_commodity="hydrogen", + demand_commodity_rate_units="kg/h", + ) + prob = _build_hetero_problem( + DemandFollowingControl, + plant_config, + slc_topology, + tech_config={}, + demand=100.0, + commodity_units={"electricity": "kW"}, + ) + prob.set_val("slc.grid_rated_electricity_production", 1e9) + prob.set_val("slc.electrolyzer_rated_hydrogen_production", 1e6) + prob.set_val("slc.electrolyzer_hydrogen_out", 100.0) + prob.set_val("slc.electrolyzer_electricity_consumed", 5100.0) + + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + prob.run_model() + + assert not any("No conversion ratio" in str(w.message) for w in caught) + np.testing.assert_allclose(prob.get_val("slc.grid_electricity_set_point"), 100.0 * 51.0) diff --git a/h2integrate/control/control_strategies/system_level/test/test_slc_examples.py b/h2integrate/control/control_strategies/system_level/test/test_slc_examples.py index e97e63bb2..89696019b 100644 --- a/h2integrate/control/control_strategies/system_level/test/test_slc_examples.py +++ b/h2integrate/control/control_strategies/system_level/test/test_slc_examples.py @@ -187,7 +187,7 @@ def test_slc_yes_hydrogen(subtests, temp_copy_of_example): pytest.approx( model.prob.get_val("finance_subgroup_hydrogen.LCOH", units="USD/kg"), rel=1e-6 ) - == 14.878096642042243 + == 9.46456262909423 ) diff --git a/h2integrate/core/h2integrate_model.py b/h2integrate/core/h2integrate_model.py index 91acb6e9d..459c22337 100644 --- a/h2integrate/core/h2integrate_model.py +++ b/h2integrate/core/h2integrate_model.py @@ -1,6 +1,7 @@ import re import importlib.util from enum import IntEnum +from collections import defaultdict import numpy as np import networkx as nx @@ -20,6 +21,7 @@ SLCSolverOptionsConfig, ) from h2integrate.control.control_strategies.system_level.system_level_control_base import ( + detect_commodity_converters, _get_tech_buy_price_input_name, ) @@ -632,6 +634,22 @@ def _classify_slc_technologies(self): slc_topology["tech_control_classifiers"] = upstream_tech_control_classifiers + # Detect commodity converters (e.g. electrolyzer: electricity -> hydrogen) + # so the same converter set drives backward demand propagation in the + # controller and the consumption-signal connections made below + converter_input_classifiers = {"fixed", "flexible", "dispatchable", "storage"} + converter_input_techs = { + tech + for tech, classifier in upstream_tech_control_classifiers.items() + if classifier in converter_input_classifiers + } + produced_by_tech = defaultdict(set) + for tech, commodity in tech_to_commodity: + produced_by_tech[tech].add(commodity) + slc_topology["converters"] = detect_commodity_converters( + upstream_tech_graph, converter_input_techs, produced_by_tech + ) + return slc_topology def add_system_level_controller(self, slc_topology): @@ -806,6 +824,30 @@ def add_system_level_controller(self, slc_topology): f"{tech_name}.{commodity}_set_point", ) + # --- Step 3b: Connect converter consumption signals --------------- + # For each detected converter whose input commodity has a + # controller-managed producer, wire the converter's measured + # ``{in_commodity}_consumed`` output to the controller so it can derive + # a per-timestep conversion ratio (heterogeneous-commodity control). + converters = slc_topology.get("converters", set()) + classifiers = slc_topology["tech_control_classifiers"] + converter_input_classifiers = {"fixed", "flexible", "dispatchable", "storage"} + produced_by_tech = defaultdict(set) + for tech_name, commodity in slc_topology["tech_to_commodity"]: + produced_by_tech[tech_name].add(commodity) + for in_commodity, converter_tech, _out_commodity in converters: + has_producers = any( + in_commodity in produced_by_tech[tech] + for tech, classifier in classifiers.items() + if classifier in converter_input_classifiers + ) + if not has_producers: + continue + self.plant.connect( + f"{converter_tech}.{in_commodity}_consumed", + f"system_level_controller.{converter_tech}_{in_commodity}_consumed", + ) + # --- Step 4: Connect marginal-cost inputs (cost-aware strategies) - if strategy_name in ("CostMinimizationControl", "ProfitMaximizationControl"): cost_per_tech = plant_slc_config.get("control_parameters", {}).get("cost_per_tech", {}) From 3c22f9a3b2b0109e3d1e915c44cf0eafeee84451 Mon Sep 17 00:00:00 2001 From: John Jasa Date: Fri, 14 Aug 2026 00:53:27 -0600 Subject: [PATCH 3/5] Expanding hetereo complexity --- CHANGELOG.md | 1 + .../slc_demand_following.md | 6 + .../system_level_control_base.md | 22 +- .../driver_config.yaml | 4 + .../heterogeneous_commodity.yaml | 4 + .../heterogeneous_commodity/plant_config.yaml | 105 ++++++++ .../run_heterogeneous.py | 252 ++++++++++++++++++ .../heterogeneous_commodity/tech_config.yaml | 241 +++++++++++++++++ .../system_level/test/test_slc_examples.py | 70 +++++ 9 files changed, 704 insertions(+), 1 deletion(-) create mode 100644 examples/35_system_level_control/heterogeneous_commodity/driver_config.yaml create mode 100644 examples/35_system_level_control/heterogeneous_commodity/heterogeneous_commodity.yaml create mode 100644 examples/35_system_level_control/heterogeneous_commodity/plant_config.yaml create mode 100644 examples/35_system_level_control/heterogeneous_commodity/run_heterogeneous.py create mode 100644 examples/35_system_level_control/heterogeneous_commodity/tech_config.yaml diff --git a/CHANGELOG.md b/CHANGELOG.md index 5338b2ec6..8ecf923e3 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,7 @@ - Replace HOPP with native H2I wind, solar, and battery models in example 11. Adds `percent_load_missed` and `curtailment_percent` outputs to `DemandComponentBase`, allows zero capacity in wind/solar/battery performance models. [PR 826](https://github.com/NatLabRockies/H2Integrate/pull/826) - Add heterogeneous-commodity system-level control that translates demand for one commodity into upstream set-points across converters using static per-technology conversion ratios defined in the tech config. [PR TBD](https://github.com/NatLabRockies/H2Integrate/pull/TBD) - Extend the heterogeneous-commodity control to prefer measured conversion ratios computed from each converter's consumed and produced streams per timestep, falling back to the static ratio when a measurement is unavailable. [PR TBD](https://github.com/NatLabRockies/H2Integrate/pull/TBD) + - Add the `35_system_level_control/heterogeneous_commodity` example, which serves an ammonia demand from a wind, battery, grid, electrolyzer, hydrogen-storage, and ammonia synthesis loop chain to demonstrate demand propagating from ammonia to hydrogen to electricity. The example uses profit-maximizing control so that wind always runs, the battery charges on wind surplus and discharges to cover deficits, and the grid is only dispatched to backfill the electricity that wind and the battery cannot supply, and it generates dispatch and dynamic conversion-ratio figures. [PR TBD](https://github.com/NatLabRockies/H2Integrate/pull/TBD) ## 0.9 [August 10, 2026] diff --git a/docs/control/system_level_control/slc_demand_following.md b/docs/control/system_level_control/slc_demand_following.md index 3dbd9b3b6..140bef989 100644 --- a/docs/control/system_level_control/slc_demand_following.md +++ b/docs/control/system_level_control/slc_demand_following.md @@ -92,6 +92,12 @@ The inputs for technologies classified as `feedstock` are: The `DemandFollowingControl` controller can be used in hybrid systems where technologies produce different commodities. For example, in a system where an electrolyzer produces hydrogen and the demand commodity is hydrogen, the controller can set the electricity-generating *curtailable* technologies' set-points to meet the hydrogen demand. +The controller can also propagate demand backward across converters, which are technologies whose output commodity differs from an input commodity they consume. +When a converter has a controllable producer for one of its input commodities, the demand for the converter's output commodity is translated into a demand for that input commodity using a conversion ratio. +The controller prefers a measured conversion ratio, computed per timestep from the converter's consumed and produced streams, and falls back to a static ratio supplied under `technologies..model_inputs.control_parameters.conversion_ratios` when a measurement is unavailable. +This propagation can chain across multiple converters. For example, an ammonia demand is translated into a hydrogen demand across an ammonia synthesis loop, and that hydrogen demand is then translated into an electricity demand across an electrolyzer, which the curtailable and dispatchable electricity technologies serve. +The `examples/35_system_level_control/heterogeneous_commodity` example demonstrates this end to end with a wind, battery, grid, electrolyzer, hydrogen-storage, and ammonia synthesis loop chain served against an ammonia demand. The example uses profit-maximizing control, so wind always runs, the battery charges on wind surplus and discharges to cover deficits, and the grid is only dispatched to backfill the electricity that wind and the battery cannot supply. Its runner script writes dispatch and dynamic conversion-ratio figures to an `outputs` folder. + This framework provides a starting point for hybrid energy system control but is intended to be extended with more sophisticated strategies for complex multi-commodity systems. ## Limitations diff --git a/docs/control/system_level_control/system_level_control_base.md b/docs/control/system_level_control/system_level_control_base.md index 7319cad36..471265408 100644 --- a/docs/control/system_level_control/system_level_control_base.md +++ b/docs/control/system_level_control/system_level_control_base.md @@ -18,7 +18,27 @@ Setup I/O for SLC controllers. - `_setup_tech_category()` - `_setup_feedstock_category()` - `find_converter_techs()` - - Note: this method is currently is not used but will be used for heterogeneous commodity systems. + +Heterogeneous-commodity conversion. These methods let a controller translate demand +for one commodity into demand for the upstream commodities that a converter consumes +to produce it (for example, translating ammonia demand into hydrogen demand across a +synthesis loop, and then into electricity demand across an electrolyzer). +- `detect_commodity_converters()` + - Module-level helper that identifies converter technologies (technologies whose + output commodity differs from an input commodity) that have a controllable + producer for that input commodity, so demand can propagate backward across them. +- `_build_conversion_ratios()` + - Registers each converter's `{commodity}_consumed` inputs and reads any static + conversion ratios supplied under + `technologies..model_inputs.control_parameters.conversion_ratios`. +- `_conversion_ratio()` + - Returns the per-timestep conversion ratio for a converter. It prefers a measured + ratio computed from the converter's consumed and produced streams, and falls back + to the static ratio from the technology config when a measurement is unavailable. +- `_accumulate_derived_demand()` + - Adds the demand derived through a converter (output demand times conversion ratio) + to the upstream commodity's demand, warning only when neither a measured nor a + static ratio is available. Functions for controlling components based on assigned control classifier. - `_subtract_curtailable()` diff --git a/examples/35_system_level_control/heterogeneous_commodity/driver_config.yaml b/examples/35_system_level_control/heterogeneous_commodity/driver_config.yaml new file mode 100644 index 000000000..0a1a05894 --- /dev/null +++ b/examples/35_system_level_control/heterogeneous_commodity/driver_config.yaml @@ -0,0 +1,4 @@ +name: driver_config +description: This analysis runs a wind and battery powered green ammonia plant under system-level control +general: + folder_output: outputs diff --git a/examples/35_system_level_control/heterogeneous_commodity/heterogeneous_commodity.yaml b/examples/35_system_level_control/heterogeneous_commodity/heterogeneous_commodity.yaml new file mode 100644 index 000000000..e09f3dcda --- /dev/null +++ b/examples/35_system_level_control/heterogeneous_commodity/heterogeneous_commodity.yaml @@ -0,0 +1,4 @@ +name: H2Integrate_config +driver_config: driver_config.yaml +plant_config: plant_config.yaml +technology_config: tech_config.yaml diff --git a/examples/35_system_level_control/heterogeneous_commodity/plant_config.yaml b/examples/35_system_level_control/heterogeneous_commodity/plant_config.yaml new file mode 100644 index 000000000..1d11f8ac6 --- /dev/null +++ b/examples/35_system_level_control/heterogeneous_commodity/plant_config.yaml @@ -0,0 +1,105 @@ +name: plant_config +description: This plant is located in Texas, USA. +sites: + site: + latitude: 30.6617 + longitude: -101.7096 + resources: + wind_resource: + resource_model: WTKNLRDeveloperAPIWindResource + resource_parameters: + resource_year: 2013 +# array of arrays containing left-to-right technology interconnections +technology_interconnections: + # wind output available for battery charging (electricity_in) + - [wind, battery, electricity, cable] + # wind to the combined electricity output + - [wind, elec_combiner, electricity, cable] + # battery net output to the combined electricity output + - [battery, elec_combiner, electricity, cable] + # dispatchable grid backup to the combined electricity output + - [grid, elec_combiner, electricity, cable] + # combined electricity to the electrolyzer + - [elec_combiner, electrolyzer, electricity, cable] + # electrolyzer to hydrogen storage + - [electrolyzer, h2_storage, hydrogen, pipe] + # combine the hydrogen streams from the electrolyzer and storage + - [electrolyzer, h2_combiner, hydrogen, pipe] + - [h2_storage, h2_combiner, hydrogen, pipe] + # hydrogen supply to the ammonia synthesis loop + - [h2_combiner, ammonia, hydrogen, pipe] + # nitrogen and electricity feedstocks for the synthesis loop + - [n2_feedstock, ammonia, nitrogen, pipe] + - [electricity_feedstock, ammonia, electricity, cable] + # ammonia supply to the ammonia demand + - [ammonia, ammonia_load_demand, ammonia, pipe] +resource_to_tech_connections: + # connect the wind resource to the wind technology + - [site.wind_resource, wind, wind_resource_data] +plant: + plant_life: 30 + simulation: + n_timesteps: 8760 + dt: 3600 +system_level_control: + control_strategy: ProfitMaximizationControl + demand_component: ammonia_load_demand + # Profit-maximizing dispatch. The controller honors each technology's control + # classification, which produces the desired behavior for this plant: + # - wind is flexible, so it always runs at its available (curtailable) output; + # - the battery is storage, so it charges on surplus and discharges to cover deficits; + # - the grid is a dispatchable source, so merit-order dispatch only calls on it to + # backfill the electricity that wind and the battery cannot supply. + # Ammonia is produced at every timestep where its sell price exceeds its marginal + # cost, so the whole chain runs whenever ammonia production is profitable. + control_parameters: + commodity_sell_price: 0.60 # USD/kg ammonia sell price (above the marginal cost so demand is served) + cost_per_tech: + # Ammonia marginal cost is taken from its upstream nitrogen and direct-electricity + # feedstocks; the grid, wind, and battery are not feedstocks and are not summed here. + ammonia: feedstock + solver_options: + solver_name: gauss_seidel + max_iter: 50 + convergence_tolerance: 1.0e-6 +finance_parameters: + finance_groups: + profast_lco: + finance_model: ProFastLCO + model_inputs: + params: + analysis_start_year: 2032 + installation_time: 36 # months + inflation_rate: 0.0 # 0 for nominal analysis + discount_rate: 0.09 + debt_equity_ratio: 2.62 + property_tax_and_insurance: 0.03 + total_income_tax_rate: 0.257 + capital_gains_tax_rate: 0.15 + sales_tax_rate: 0.07375 + debt_interest_rate: 0.07 + debt_type: Revolving debt + loan_period_if_used: 0 + cash_onhand_months: 1 + admin_expense: 0.00 + capital_items: + depr_type: MACRS + depr_period: 5 + refurb: [0.] + finance_subgroups: + ammonia: + commodity: ammonia + commodity_stream: ammonia + finance_groups: [profast_lco] + technologies: + - wind + - battery + - grid + - electrolyzer + - h2_storage + - n2_feedstock + - electricity_feedstock + - ammonia + cost_adjustment_parameters: + cost_year_adjustment_inflation: 0.025 # used to adjust modeled costs to target_dollar_year + target_dollar_year: 2022 diff --git a/examples/35_system_level_control/heterogeneous_commodity/run_heterogeneous.py b/examples/35_system_level_control/heterogeneous_commodity/run_heterogeneous.py new file mode 100644 index 000000000..dbc427373 --- /dev/null +++ b/examples/35_system_level_control/heterogeneous_commodity/run_heterogeneous.py @@ -0,0 +1,252 @@ +import os +from pathlib import Path + +import numpy as np +import matplotlib + + +# Use a non-interactive backend so the figures render when the script is run headless. +matplotlib.use("Agg") +import matplotlib.pyplot as plt + +from h2integrate import EXAMPLE_DIR +from h2integrate.core.h2integrate_model import H2IntegrateModel + + +EXAMPLE_FOLDER = EXAMPLE_DIR / "35_system_level_control" / "heterogeneous_commodity" +os.chdir(EXAMPLE_FOLDER) + +################################## +# Create an H2I model with an ammonia demand served by a heterogeneous commodity chain. +# The profit-maximizing system-level controller translates ammonia demand backward into +# hydrogen demand (across the synthesis loop) and then into electricity demand (across the +# electrolyzer). Wind runs whenever it is available, the battery charges on wind surplus and +# discharges to cover short deficits, and the grid is only dispatched to backfill the +# electricity that wind and the battery cannot supply. +h2i = H2IntegrateModel("heterogeneous_commodity.yaml") + +h2i.setup() + +# Run the model +h2i.run() + +# Post-process the results +h2i.post_process() + + +################################## +# Plot the resulting dispatch, backward demand propagation, and the dynamic +# conversion ratios that the controller uses to translate demand across commodities. +def _get(name, units=None): + """Return a flattened numpy array for a promoted model output.""" + return np.asarray(h2i.prob.get_val(name, units=units)).flatten() + + +def make_plots(figure_dir): + figure_dir = Path(figure_dir) + figure_dir.mkdir(parents=True, exist_ok=True) + + n_timesteps = h2i.prob.get_val("ammonia.ammonia_out", units="kg/h").size + hours = np.arange(n_timesteps) + + # A representative week (7 days) chosen to show variable wind and grid backup. + week_start = 3000 + week = slice(week_start, week_start + 168) + week_hours = hours[week] - week_start + + # Electricity streams (MW). + wind_mw = _get("wind.electricity_out", "MW") + grid_mw = _get("grid.electricity_out", "MW") + battery_mw = _get("battery.electricity_out", "MW") + electrolyzer_load_mw = _get("electrolyzer.electricity_consumed", "MW") + synloop_load_mw = _get("ammonia.electricity_consumed", "MW") + + # Hydrogen streams (kg/h). + h2_produced = _get("electrolyzer.hydrogen_out", "kg/h") + h2_to_synloop = _get("ammonia.hydrogen_consumed", "kg/h") + + # Ammonia streams (kg/h). + ammonia_out = _get("ammonia.ammonia_out", "kg/h") + ammonia_demand = _get("ammonia_load_demand.ammonia_demand_out", "kg/h") + + # Controller set points that show backward propagation of demand. + ammonia_set_point = _get("system_level_controller.ammonia_ammonia_set_point", "kg/h") + hydrogen_set_point = _get("system_level_controller.electrolyzer_hydrogen_set_point", "kg/h") + total_electricity_load_mw = electrolyzer_load_mw + synloop_load_mw + + # Measured (dynamic) conversion ratios the controller derives per timestep. + electrolyzer_load_kw = _get("electrolyzer.electricity_consumed", "kW") + synloop_load_kw = _get("ammonia.electricity_consumed", "kW") + with np.errstate(divide="ignore", invalid="ignore"): + electricity_per_hydrogen = np.where( + h2_produced > 1e-6, electrolyzer_load_kw / h2_produced, np.nan + ) + hydrogen_per_ammonia = np.where(ammonia_out > 1e-6, h2_to_synloop / ammonia_out, np.nan) + electricity_per_ammonia = np.where( + ammonia_out > 1e-6, synloop_load_kw / ammonia_out, np.nan + ) + + # ----------------------------------------------------------------- + # Figure 1: the commodity cascade over a representative week. + # Electricity (wind + grid) drives hydrogen, which drives ammonia. + # ----------------------------------------------------------------- + fig, axes = plt.subplots(3, 1, figsize=(10, 9), sharex=True) + + axes[0].stackplot( + week_hours, + wind_mw[week], + np.clip(battery_mw[week], 0.0, None), + grid_mw[week], + labels=["Wind", "Battery discharge", "Grid (firm)"], + colors=["#4C9F70", "#F2C14E", "#8C8C8C"], + ) + axes[0].plot( + week_hours, + total_electricity_load_mw[week], + color="black", + lw=1.5, + label="Electrolyzer + synloop load", + ) + axes[0].set_ylabel("Electricity (MW)") + axes[0].set_title("Electricity supply: wind prioritized, grid provides firm backup") + axes[0].legend(loc="upper right", ncol=2, fontsize=8) + + axes[1].plot(week_hours, h2_produced[week], color="#3B7DD8", label="Electrolyzer output") + axes[1].plot( + week_hours, h2_to_synloop[week], color="#D8663B", ls="--", label="Synloop consumption" + ) + axes[1].set_ylabel("Hydrogen (kg/h)") + axes[1].set_title("Hydrogen: electrolyzer output feeds the ammonia synthesis loop") + axes[1].legend(loc="upper right", fontsize=8) + + axes[2].plot(week_hours, ammonia_out[week], color="#7B4FA3", label="Ammonia production") + axes[2].plot(week_hours, ammonia_demand[week], color="black", ls="--", label="Ammonia demand") + axes[2].set_ylabel("Ammonia (kg/h)") + axes[2].set_xlabel("Hour of representative week") + axes[2].set_title("Ammonia: production tracks the demand that drives the whole chain") + axes[2].legend(loc="lower right", fontsize=8) + + fig.suptitle("Heterogeneous-commodity dispatch: electricity -> hydrogen -> ammonia", y=0.995) + fig.tight_layout() + fig.savefig(figure_dir / "dispatch_cascade.png", dpi=150) + plt.close(fig) + + # ----------------------------------------------------------------- + # Figure 2: backward demand propagation. A firm ammonia demand is + # translated into a firm hydrogen demand and then a firm electricity + # load. The demands are near-constant, so we show their magnitudes and + # annotate the conversion ratios that connect them. + # ----------------------------------------------------------------- + fig, ax_left = plt.subplots(figsize=(10, 5)) + ax_left.plot( + week_hours, ammonia_set_point[week], color="#7B4FA3", label="Ammonia set point (kg/h)" + ) + ax_left.plot( + week_hours, hydrogen_set_point[week], color="#D8663B", label="Hydrogen set point (kg/h)" + ) + ax_left.set_xlabel("Hour of representative week") + ax_left.set_ylabel("Commodity set point (kg/h)") + ax_left.set_ylim(0.0, 1.15 * float(ammonia_set_point[week].max())) + + ax_right = ax_left.twinx() + ax_right.plot( + week_hours, + total_electricity_load_mw[week], + color="#3B7DD8", + ls="--", + label="Derived electricity load (MW)", + ) + ax_right.set_ylabel("Electricity load (MW)") + ax_right.set_ylim(0.0, 1.15 * float(total_electricity_load_mw[week].max())) + + # Annotate the conversion ratios that link the three firm demand levels. + mean_hydrogen_per_ammonia = float(np.nanmean(hydrogen_per_ammonia)) + mean_electricity_per_ammonia_mw = float(total_electricity_load_mw.mean() / ammonia_out.mean()) + ax_left.annotate( + f"x {mean_hydrogen_per_ammonia:0.3f} kg H2 / kg NH3", + xy=(0.5, 0.62), + xycoords="axes fraction", + color="#D8663B", + fontsize=9, + ha="center", + ) + ax_left.annotate( + f"x {mean_electricity_per_ammonia_mw * 1000:0.2f} kWh / kg NH3 (total)", + xy=(0.5, 0.12), + xycoords="axes fraction", + color="#3B7DD8", + fontsize=9, + ha="center", + ) + + lines_left, labels_left = ax_left.get_legend_handles_labels() + lines_right, labels_right = ax_right.get_legend_handles_labels() + ax_left.legend( + lines_left + lines_right, labels_left + labels_right, loc="center right", fontsize=8 + ) + ax_left.set_title("Backward demand propagation: ammonia -> hydrogen -> electricity") + fig.tight_layout() + fig.savefig(figure_dir / "demand_propagation.png", dpi=150) + plt.close(fig) + + # ----------------------------------------------------------------- + # Figure 3: dynamic conversion ratios across the year. The controller + # prefers these measured ratios over the static seed values. The + # electrolyzer ratio drifts up as the stack degrades over the year. + # ----------------------------------------------------------------- + fig, axes = plt.subplots(3, 1, figsize=(10, 9), sharex=True) + + axes[0].plot(hours, electricity_per_hydrogen, color="#3B7DD8", lw=0.8) + axes[0].axhline(51.0, color="black", ls=":", lw=1.0, label="Static seed (51 kWh/kg)") + axes[0].set_ylabel("kWh / kg H2") + axes[0].set_title( + "Electrolyzer measured ratio (electricity per hydrogen) drifts up with degradation" + ) + axes[0].legend(loc="upper left", fontsize=8) + + axes[1].plot(hours, hydrogen_per_ammonia, color="#D8663B", lw=0.8) + axes[1].axhline(0.2, color="black", ls=":", lw=1.0, label="Static seed (0.2 kg/kg)") + axes[1].set_ylabel("kg H2 / kg NH3") + axes[1].set_title("Synloop measured ratio (hydrogen per ammonia)") + axes[1].legend(loc="upper left", fontsize=8) + + axes[2].plot(hours, electricity_per_ammonia, color="#7B4FA3", lw=0.8) + axes[2].axhline(0.530645243, color="black", ls=":", lw=1.0, label="Static seed (0.53 kWh/kg)") + axes[2].set_ylabel("kWh / kg NH3") + axes[2].set_xlabel("Hour of year") + axes[2].set_title("Synloop measured ratio (electricity per ammonia)") + axes[2].legend(loc="upper left", fontsize=8) + + fig.suptitle("Dynamic conversion ratios used to translate demand across converters", y=0.995) + fig.tight_layout() + fig.savefig(figure_dir / "conversion_ratios.png", dpi=150) + plt.close(fig) + + # ----------------------------------------------------------------- + # Figure 4: annual electricity source mix and the resulting LCOA. + # ----------------------------------------------------------------- + wind_energy = wind_mw.sum() + grid_energy = grid_mw.sum() + battery_energy = np.clip(battery_mw, 0.0, None).sum() + lcoa = float(h2i.prob.get_val("finance_subgroup_ammonia.LCOA", units="USD/kg")[0]) + + fig, ax = plt.subplots(figsize=(6, 6)) + ax.pie( + [wind_energy, battery_energy, grid_energy], + labels=["Wind", "Battery", "Grid (firm)"], + colors=["#4C9F70", "#F2C14E", "#8C8C8C"], + autopct="%1.1f%%", + startangle=90, + ) + ax.set_title( + f"Annual electricity supplied to the chain\nAmmonia LCOA = ${lcoa:0.2f}/kg", fontsize=11 + ) + fig.tight_layout() + fig.savefig(figure_dir / "electricity_source_mix.png", dpi=150) + plt.close(fig) + + return figure_dir + + +figures = make_plots(EXAMPLE_FOLDER / "outputs") +print(f"Saved dispatch and conversion-ratio figures to {figures}") diff --git a/examples/35_system_level_control/heterogeneous_commodity/tech_config.yaml b/examples/35_system_level_control/heterogeneous_commodity/tech_config.yaml new file mode 100644 index 000000000..827d29326 --- /dev/null +++ b/examples/35_system_level_control/heterogeneous_commodity/tech_config.yaml @@ -0,0 +1,241 @@ +name: technology_config +description: Wind and battery power an electrolyzer that feeds an ammonia synthesis loop, with the system-level controller + translating ammonia demand backward into hydrogen and electricity demand. +technologies: + wind: + performance_model: + model: PYSAMWindPlantPerformanceModel + cost_model: + model: ATBWindPlantCostModel + model_inputs: + performance_parameters: + num_turbines: 20 + turbine_rating_kw: 6000 + hub_height: 115 + rotor_diameter: 170 + create_model_from: default + config_name: WindPowerSingleOwner + pysam_options: + Farm: + wind_farm_wake_model: 0 + Losses: + ops_strategies_loss: 10.0 + layout: + layout_mode: basicgrid + layout_options: + row_D_spacing: 5.0 + turbine_D_spacing: 5.0 + rotation_angle_deg: 0.0 + row_phase_offset: 0.0 + layout_shape: square + cost_parameters: + capex_per_kW: 1300 + opex_per_kW_per_year: 39 + cost_year: 2022 + battery: + performance_model: + model: StoragePerformanceModel + cost_model: + model: GenericStorageCostModel + model_inputs: + shared_parameters: + commodity: electricity + commodity_rate_units: kW + max_charge_rate: 10000 # kW (10 MW) + max_capacity: 40000 # kWh (40 MWh, 4-hour duration) + init_soc_fraction: 0.5 + max_soc_fraction: 1.0 + min_soc_fraction: 0.1 + performance_parameters: + round_trip_efficiency: 0.90 + demand_profile: 10000 # kW, required by storage base config + cost_parameters: + cost_year: 2022 + capacity_capex: 310 # $/kWh + charge_capex: 311 # $/kW + opex_fraction: 0.025 + elec_combiner: + performance_model: + model: GenericCombinerPerformanceModel + model_inputs: + performance_parameters: + commodity: electricity + commodity_rate_units: kW + in_streams: 3 + electrolyzer: + performance_model: + model: ECOElectrolyzerPerformanceModel + cost_model: + model: SingliticoCostModel + model_inputs: + shared_parameters: + location: onshore + electrolyzer_capex: 1295 # $/kW overnight installed capital costs for a 1 MW system in 2022 USD/kW + performance_parameters: + size_mode: normal + n_clusters: 15 # 15 x 3 MW = 45 MW, sized so the load sits within the wind's reach + cluster_rating_MW: 3 + eol_eff_percent_loss: 10 # eol defined as x% change in efficiency from bol + uptime_hours_until_eol: 80000. # number of 'on' hours until electrolyzer reaches eol + include_degradation_penalty: true # include degradation + turndown_ratio: 0.1 # turndown_ratio = minimum_cluster_power/cluster_rating_MW + financial_parameters: + capital_items: + depr_period: 7 + replacement_cost_percent: 0.15 # percent of capex - H2A default case + # Static seed ratio for the electricity -> hydrogen converter. The controller + # prefers the measured ratio (electricity_consumed / hydrogen_out) once the + # solver has values and uses this seed on iteration zero and for zero-output + # timesteps. + control_parameters: + conversion_ratios: + electricity_per_hydrogen: 51.0 # kWh electricity per kg hydrogen + h2_storage: + performance_model: + model: StoragePerformanceModel + cost_model: + model: GenericStorageCostModel + model_inputs: + shared_parameters: + commodity: hydrogen + commodity_rate_units: kg/h + max_charge_rate: 1200.0 # kg/time step + max_capacity: 3700.0 # kg + performance_parameters: + max_soc_fraction: 1.0 # fraction (0-1) + min_soc_fraction: 0.1 # fraction (0-1) + init_soc_fraction: 0.1 # fraction (0-1) + max_discharge_rate: 1200.0 # kg/time step + charge_efficiency: 1.0 # fraction (0-1) + discharge_efficiency: 1.0 # fraction (0-1) + demand_profile: 4000.0 # kg/h (see commodity_rate_units) + cost_parameters: + cost_year: 2022 + capacity_capex: 100 # $/kg + charge_capex: 100 # $/kg/h + opex_fraction: 0.025 + h2_combiner: + performance_model: + model: GenericCombinerPerformanceModel + model_inputs: + performance_parameters: + commodity: hydrogen + commodity_rate_units: kg/h + n2_feedstock: + performance_model: + model: FeedstockPerformanceModel + cost_model: + model: FeedstockCostModel + model_inputs: + shared_parameters: + commodity: nitrogen + commodity_rate_units: t/h + performance_parameters: + rated_capacity: 6.0 # metric tonnes of N2/hour (oversized so ammonia is hydrogen-limited) + cost_parameters: + cost_year: 2022 + price: 5.0 + annual_cost: 0. + start_up_cost: 0.0 + electricity_feedstock: + performance_model: + model: FeedstockPerformanceModel + cost_model: + model: FeedstockCostModel + model_inputs: + shared_parameters: + commodity: electricity + commodity_rate_units: MW + performance_parameters: + rated_capacity: 5.0 # MW (covers the synloop's direct electricity draw) + cost_parameters: + cost_year: 2022 + price: 25.0 # USD/(MW*h) == 0.025 USD/kWh, matching the grid buy price + annual_cost: 0. + start_up_cost: 0.0 + grid: + performance_model: + model: GridPerformanceModel + cost_model: + model: GridCostModel + model_inputs: + shared_parameters: + interconnection_size: 50000.0 # kW (50 MW) firm backup for the electrolyzer + cost_parameters: + cost_year: 2022 + electricity_buy_price: 0.025 # USD/kWh firm grid power + interconnection_capex_per_kw: 100.0 + interconnection_opex_per_kw: 5.0 + fixed_interconnection_cost: 0.0 + ammonia: + performance_model: + model: AmmoniaSynLoopPerformanceModel + cost_model: + model: AmmoniaSynLoopCostModel + model_inputs: + shared_parameters: + production_capacity: 4000.0 # kg ammonia per hour + catalyst_consumption_rate: 0.000091295354067341 + catalyst_replacement_interval: 3 + performance_parameters: + size_mode: normal + capacity_factor: 0.9 + energy_demand: 0.530645243 # kWh electricity per kg ammonia + heat_output: 0.8299956 + feed_gas_t: 25.8 + feed_gas_p: 20 + feed_gas_x_n2: 0.25 + feed_gas_x_h2: 0.75 + feed_gas_mass_ratio: 1.13 + purge_gas_t: 7.5 + purge_gas_p: 275 + purge_gas_x_n2: 0.26 + purge_gas_x_h2: 0.68 + purge_gas_x_ar: 0.02 + purge_gas_x_nh3: 0.04 + purge_gas_mass_ratio: 0.07 + # Dynamics disabled for a well-behaved system-level control fixed-point loop. + turndown_ratio: 0.0 + ramp_up_rate_fraction: 1.0 + ramp_down_rate_fraction: 1.0 + include_cold_start: false + include_warm_start: false + # Static seed ratios for the two controllable inputs of the synloop + # (hydrogen and electricity). Nitrogen is supplied by a feedstock and is not + # controller-managed, so it needs no ratio. + control_parameters: + conversion_ratios: + hydrogen_per_ammonia: 0.2 # kg hydrogen per kg ammonia + electricity_per_ammonia: 0.530645243 # kWh electricity per kg ammonia + cost_parameters: + baseline_capacity: 52777.6 + base_cost_year: 2016 + capex_scaling_exponent: 0.6 + labor_scaling_exponent: 0.25 + asu_capex_base: 236920646 + synloop_capex_base: 302460908 + heat_capex_base: 7069100 + cool_capex_base: 4799200 + other_eqpt_capex_base: 0 + land_capex_base: 4112701.84103543 + deprec_noneq_capex_rate: 0.42 + labor_rate_base: 57 + num_workers_base: 50 + hours_yr: 2080 + gen_admin: 0.2 + prop_tax_ins: 0.02 + maint_rep: 0.005 + oxygen_byproduct_rate: 0.29405077250145 + water_consumption_rate: 0.049236824 + rebuild_cost_base: 0 + cooling_water_cost_base: 0.000113349938601175 + catalyst_cost_base: 23.19977341 + oxygen_price_base: 0.0285210891617726 + ammonia_load_demand: + performance_model: + model: GenericDemandComponent + model_inputs: + performance_parameters: + commodity: ammonia + commodity_rate_units: kg/h + demand_profile: 4000.0 # constant ammonia demand (kg/h) diff --git a/h2integrate/control/control_strategies/system_level/test/test_slc_examples.py b/h2integrate/control/control_strategies/system_level/test/test_slc_examples.py index 89696019b..0f7524fac 100644 --- a/h2integrate/control/control_strategies/system_level/test/test_slc_examples.py +++ b/h2integrate/control/control_strategies/system_level/test/test_slc_examples.py @@ -403,3 +403,73 @@ def test_slc_upstream_demand(subtests, temp_copy_of_example): # check that no hydrogen systems are in model.prob.get_val(slc_h2s_output_var, units="kg/h") assert f"Variable '{slc_h2s_output_var}' not found. " in str(excinfo.value) + + +@pytest.mark.integration +@pytest.mark.parametrize( + "example_folder,resource_example_folder", + [("35_system_level_control/heterogeneous_commodity", None)], +) +def test_slc_heterogeneous_commodity(subtests, temp_copy_of_example): + example_folder = temp_copy_of_example + + model = H2IntegrateModel(example_folder / "heterogeneous_commodity.yaml") + + model.run() + + with subtests.test("Ammonia set point follows demand"): + assert np.all( + model.prob.get_val("system_level_controller.ammonia_ammonia_set_point", units="kg/h") + == model.prob.get_val("ammonia_load_demand.ammonia_demand_out", units="kg/h") + ) + + with subtests.test("Hydrogen demand propagated from ammonia"): + # The converter's ammonia set point drives a nonzero electrolyzer hydrogen set point. + assert ( + model.prob.get_val( + "system_level_controller.electrolyzer_hydrogen_set_point", units="kg/h" + ).sum() + > 0.0 + ) + + with subtests.test("Electricity demand propagated to wind and battery"): + # Hydrogen demand across the electrolyzer propagates to an electricity demand, so the + # electricity producers appear in the controller as set-point outputs. + assert ( + model.prob.get_val( + "system_level_controller.wind_electricity_set_point", units="kW" + ).sum() + > 0.0 + ) + assert ( + model.prob.get_val( + "system_level_controller.battery_electricity_set_point", units="kW" + ).size + > 0 + ) + + with subtests.test("Ammonia production does not exceed demand"): + assert ( + model.prob.get_val("ammonia.ammonia_out", units="kg/h").sum() + <= model.prob.get_val("ammonia_load_demand.ammonia_demand_out", units="kg/h").sum() + + 1e-6 + ) + + with subtests.test("Ammonia annual production"): + assert ( + pytest.approx(35036000.0, rel=1e-5) + == model.prob.get_val("ammonia.ammonia_out", units="kg/h").sum() + ) + + with subtests.test("Electrolyzer hydrogen set point total"): + assert ( + pytest.approx(7030464.51294698, rel=1e-5) + == model.prob.get_val( + "system_level_controller.electrolyzer_hydrogen_set_point", units="kg/h" + ).sum() + ) + + with subtests.test("LCOA"): + assert pytest.approx(2.0024798555505585, rel=1e-5) == model.prob.get_val( + "finance_subgroup_ammonia.LCOA", units="USD/kg" + ) From 0306807456668d8fd3f39768c85d57a97795e9d4 Mon Sep 17 00:00:00 2001 From: John Jasa Date: Fri, 14 Aug 2026 09:42:00 -0600 Subject: [PATCH 4/5] Reworking where the commodities are set --- .../heterogeneous_commodity/plant_config.yaml | 11 +++ .../run_heterogeneous.py | 2 +- .../heterogeneous_commodity/tech_config.yaml | 16 +---- .../profit_maximization_control.py | 1 + .../system_level/system_level_control_base.py | 27 ++++---- .../system_level/test/test_slc_controllers.py | 67 +++++++++++-------- .../system_level/test/test_slc_examples.py | 2 +- 7 files changed, 69 insertions(+), 57 deletions(-) diff --git a/examples/35_system_level_control/heterogeneous_commodity/plant_config.yaml b/examples/35_system_level_control/heterogeneous_commodity/plant_config.yaml index 1d11f8ac6..ebfae07de 100644 --- a/examples/35_system_level_control/heterogeneous_commodity/plant_config.yaml +++ b/examples/35_system_level_control/heterogeneous_commodity/plant_config.yaml @@ -58,6 +58,17 @@ system_level_control: # Ammonia marginal cost is taken from its upstream nitrogen and direct-electricity # feedstocks; the grid, wind, and battery are not feedstocks and are not summed here. ammonia: feedstock + # Static seed ratios for the commodity converters, keyed by technology. The + # controller prefers the measured ratio (input_consumed / output_produced) + # once the solver has values and uses these seeds on iteration zero and for + # zero-output timesteps. Feedstock-supplied inputs (for example nitrogen into + # the synloop) are not controller-managed and need no ratio. + conversion_parameters: + electrolyzer: + electricity_per_hydrogen: 51.0 # kWh electricity per kg hydrogen + ammonia: + hydrogen_per_ammonia: 0.2 # kg hydrogen per kg ammonia + electricity_per_ammonia: 0.530645243 # kWh electricity per kg ammonia solver_options: solver_name: gauss_seidel max_iter: 50 diff --git a/examples/35_system_level_control/heterogeneous_commodity/run_heterogeneous.py b/examples/35_system_level_control/heterogeneous_commodity/run_heterogeneous.py index dbc427373..b2331df29 100644 --- a/examples/35_system_level_control/heterogeneous_commodity/run_heterogeneous.py +++ b/examples/35_system_level_control/heterogeneous_commodity/run_heterogeneous.py @@ -197,7 +197,7 @@ def make_plots(figure_dir): fig, axes = plt.subplots(3, 1, figsize=(10, 9), sharex=True) axes[0].plot(hours, electricity_per_hydrogen, color="#3B7DD8", lw=0.8) - axes[0].axhline(51.0, color="black", ls=":", lw=1.0, label="Static seed (51 kWh/kg)") + # axes[0].axhline(51.0, color="black", ls=":", lw=1.0, label="Static seed (51 kWh/kg)") axes[0].set_ylabel("kWh / kg H2") axes[0].set_title( "Electrolyzer measured ratio (electricity per hydrogen) drifts up with degradation" diff --git a/examples/35_system_level_control/heterogeneous_commodity/tech_config.yaml b/examples/35_system_level_control/heterogeneous_commodity/tech_config.yaml index 827d29326..0274db901 100644 --- a/examples/35_system_level_control/heterogeneous_commodity/tech_config.yaml +++ b/examples/35_system_level_control/heterogeneous_commodity/tech_config.yaml @@ -75,7 +75,7 @@ technologies: size_mode: normal n_clusters: 15 # 15 x 3 MW = 45 MW, sized so the load sits within the wind's reach cluster_rating_MW: 3 - eol_eff_percent_loss: 10 # eol defined as x% change in efficiency from bol + eol_eff_percent_loss: 30 # big degradation to demonstrate the changing conversion ratios uptime_hours_until_eol: 80000. # number of 'on' hours until electrolyzer reaches eol include_degradation_penalty: true # include degradation turndown_ratio: 0.1 # turndown_ratio = minimum_cluster_power/cluster_rating_MW @@ -83,13 +83,6 @@ technologies: capital_items: depr_period: 7 replacement_cost_percent: 0.15 # percent of capex - H2A default case - # Static seed ratio for the electricity -> hydrogen converter. The controller - # prefers the measured ratio (electricity_consumed / hydrogen_out) once the - # solver has values and uses this seed on iteration zero and for zero-output - # timesteps. - control_parameters: - conversion_ratios: - electricity_per_hydrogen: 51.0 # kWh electricity per kg hydrogen h2_storage: performance_model: model: StoragePerformanceModel @@ -200,13 +193,6 @@ technologies: ramp_down_rate_fraction: 1.0 include_cold_start: false include_warm_start: false - # Static seed ratios for the two controllable inputs of the synloop - # (hydrogen and electricity). Nitrogen is supplied by a feedstock and is not - # controller-managed, so it needs no ratio. - control_parameters: - conversion_ratios: - hydrogen_per_ammonia: 0.2 # kg hydrogen per kg ammonia - electricity_per_ammonia: 0.530645243 # kWh electricity per kg ammonia cost_parameters: baseline_capacity: 52777.6 base_cost_year: 2016 diff --git a/h2integrate/control/control_strategies/system_level/profit_maximization_control.py b/h2integrate/control/control_strategies/system_level/profit_maximization_control.py index 7e492d8a5..14db24ea2 100644 --- a/h2integrate/control/control_strategies/system_level/profit_maximization_control.py +++ b/h2integrate/control/control_strategies/system_level/profit_maximization_control.py @@ -10,6 +10,7 @@ class ProfitMaximizationControlConfig(BaseConfig): commodity_sell_price: float = field(default=0.0) cost_per_tech: dict = field(default={}) + conversion_parameters: dict = field(default={}) class ProfitMaximizationControl(SystemLevelControlBase): diff --git a/h2integrate/control/control_strategies/system_level/system_level_control_base.py b/h2integrate/control/control_strategies/system_level/system_level_control_base.py index 9dce96b74..b32b6704f 100644 --- a/h2integrate/control/control_strategies/system_level/system_level_control_base.py +++ b/h2integrate/control/control_strategies/system_level/system_level_control_base.py @@ -670,9 +670,9 @@ def _build_conversion_ratios(self): - ``self._converters``: set of ``(in_commodity, tech_name, out_commodity)`` tuples (empty for single-commodity systems). - ``self.conversion_ratios``: mapping of ``(tech_name, in_commodity, - out_commodity)`` to a float static ratio read from the tech config at - ``technologies..model_inputs.control_parameters. - conversion_ratios._per_``. + out_commodity)`` to a float static ratio read from the plant config at + ``system_level_control.control_parameters.conversion_parameters. + ._per_``. - ``self._converter_consumed_names``: mapping of the same key to the ``{tech}_{in_commodity}_consumed`` input registered for the dynamic (measured) ratio path. @@ -699,17 +699,18 @@ def _build_conversion_ratios(self): ) self._converters = set(converters) - # Read each converter's static input-per-output ratio from the tech config + # Read each converter's static input-per-output ratio from the plant + # config's system-level-control ``conversion_parameters`` block, keyed by + # tech name with ``_per_`` entries. self.conversion_ratios = {} self._missing_ratio_warned = set() - technologies = self.options["tech_config"].get("technologies", {}) + conversion_parameters = ( + self.options["plant_config"]["system_level_control"] + .get("control_parameters", {}) + .get("conversion_parameters", {}) + ) for in_commodity, tech_name, out_commodity in self._converters: - control_params = ( - technologies.get(tech_name, {}) - .get("model_inputs", {}) - .get("control_parameters", {}) - ) - ratios = control_params.get("conversion_ratios", {}) + ratios = conversion_parameters.get(tech_name, {}) key = f"{in_commodity}_per_{out_commodity}" if key in ratios: self.conversion_ratios[(tech_name, in_commodity, out_commodity)] = float( @@ -967,8 +968,8 @@ def _accumulate_derived_demand( f"not translate '{out_commodity}' demand into '{in_commodity}' demand, " f"and upstream '{in_commodity}' technologies keep their default " f"dispatch. Connect the converter's '{in_commodity}_consumed' output or " - f"define technologies.{tech_name}.model_inputs.control_parameters." - f"conversion_ratios.{in_commodity}_per_{out_commodity} in the tech config " + f"define system_level_control.control_parameters.conversion_parameters." + f"{tech_name}.{in_commodity}_per_{out_commodity} in the plant config " f"to enable heterogeneous-commodity control.", stacklevel=2, ) diff --git a/h2integrate/control/control_strategies/system_level/test/test_slc_controllers.py b/h2integrate/control/control_strategies/system_level/test/test_slc_controllers.py index cf84fed9e..5393193ea 100644 --- a/h2integrate/control/control_strategies/system_level/test/test_slc_controllers.py +++ b/h2integrate/control/control_strategies/system_level/test/test_slc_controllers.py @@ -729,23 +729,24 @@ def test_feedstock_no_feedstock_raises(self): # --------------------------------------------------------------------------- # Heterogeneous-commodity dispatch (backward demand propagation) # --------------------------------------------------------------------------- -def _tech_config_with_ratios(ratios_by_tech): - """Build a tech_config carrying static conversion ratios. +def _plant_config_with_conversion_parameters(plant_config, conversion_parameters): + """Attach static conversion ratios to a plant config. Args: - ratios_by_tech (dict): Mapping of ``tech_name`` to a dict of + plant_config (dict): Plant config to extend in place. + conversion_parameters (dict): Mapping of ``tech_name`` to a dict of ``{"_per_": ratio}`` entries. Returns: - dict: A ``tech_config`` with the nested ``model_inputs. - control_parameters.conversion_ratios`` structure the base class reads. + dict: The same ``plant_config`` with the ``system_level_control. + control_parameters.conversion_parameters`` structure the base class + reads. """ - return { - "technologies": { - tech: {"model_inputs": {"control_parameters": {"conversion_ratios": ratios}}} - for tech, ratios in ratios_by_tech.items() - } - } + control_parameters = plant_config.setdefault("system_level_control", {}).setdefault( + "control_parameters", {} + ) + control_parameters["conversion_parameters"] = conversion_parameters + return plant_config def _build_hetero_problem( @@ -829,17 +830,18 @@ def test_detect_converters_chain(self): slc_topology = _build_slc_topology( tech_graph, classifiers, demand_commodity="ammonia", demand_commodity_rate_units="kg/h" ) - tech_config = _tech_config_with_ratios( + plant_config = _plant_config_with_conversion_parameters( + plant_config, { "electrolyzer": {"electricity_per_hydrogen": 51.0}, "synloop": {"hydrogen_per_ammonia": 0.18}, - } + }, ) prob = _build_hetero_problem( DemandFollowingControl, plant_config, slc_topology, - tech_config, + {}, demand=100.0, commodity_units={"electricity": "kW", "hydrogen": "kg/h"}, ) @@ -863,12 +865,14 @@ def test_single_converter_static_propagation(self): demand_commodity="hydrogen", demand_commodity_rate_units="kg/h", ) - tech_config = _tech_config_with_ratios({"electrolyzer": {"electricity_per_hydrogen": 51.0}}) + plant_config = _plant_config_with_conversion_parameters( + plant_config, {"electrolyzer": {"electricity_per_hydrogen": 51.0}} + ) prob = _build_hetero_problem( DemandFollowingControl, plant_config, slc_topology, - tech_config, + {}, demand=100.0, commodity_units={"electricity": "kW"}, ) @@ -897,17 +901,18 @@ def test_chained_converter_static_propagation(self): demand_commodity="ammonia", demand_commodity_rate_units="kg/h", ) - tech_config = _tech_config_with_ratios( + plant_config = _plant_config_with_conversion_parameters( + plant_config, { "electrolyzer": {"electricity_per_hydrogen": 51.0}, "synloop": {"hydrogen_per_ammonia": 0.18}, - } + }, ) prob = _build_hetero_problem( DemandFollowingControl, plant_config, slc_topology, - tech_config, + {}, demand=100.0, commodity_units={"electricity": "kW", "hydrogen": "kg/h"}, ) @@ -942,12 +947,14 @@ def test_derived_demand_reuses_flexible_and_dispatchable(self): demand_commodity="hydrogen", demand_commodity_rate_units="kg/h", ) - tech_config = _tech_config_with_ratios({"electrolyzer": {"electricity_per_hydrogen": 50.0}}) + plant_config = _plant_config_with_conversion_parameters( + plant_config, {"electrolyzer": {"electricity_per_hydrogen": 50.0}} + ) prob = _build_hetero_problem( DemandFollowingControl, plant_config, slc_topology, - tech_config, + {}, demand=100.0, upstream_out={("wind", "electricity"): 1000.0}, commodity_units={"electricity": "kW"}, @@ -1016,12 +1023,14 @@ def test_cost_min_merit_order_at_derived_level(self): demand_commodity="hydrogen", demand_commodity_rate_units="kg/h", ) - tech_config = _tech_config_with_ratios({"electrolyzer": {"electricity_per_hydrogen": 50.0}}) + plant_config = _plant_config_with_conversion_parameters( + plant_config, {"electrolyzer": {"electricity_per_hydrogen": 50.0}} + ) prob = _build_hetero_problem( CostMinimizationControl, plant_config, slc_topology, - tech_config, + {}, demand=100.0, commodity_units={"electricity": "kW"}, ) @@ -1051,12 +1060,14 @@ def test_dynamic_ratio_overrides_static(self): demand_commodity_rate_units="kg/h", ) # Static ratio is 50, but the measured ratio (5100 / 100 = 51) should win. - tech_config = _tech_config_with_ratios({"electrolyzer": {"electricity_per_hydrogen": 50.0}}) + plant_config = _plant_config_with_conversion_parameters( + plant_config, {"electrolyzer": {"electricity_per_hydrogen": 50.0}} + ) prob = _build_hetero_problem( DemandFollowingControl, plant_config, slc_topology, - tech_config, + {}, demand=100.0, commodity_units={"electricity": "kW"}, ) @@ -1083,12 +1094,14 @@ def test_dynamic_ratio_time_varying_with_zero_output_fallback(self): demand_commodity="hydrogen", demand_commodity_rate_units="kg/h", ) - tech_config = _tech_config_with_ratios({"electrolyzer": {"electricity_per_hydrogen": 50.0}}) + plant_config = _plant_config_with_conversion_parameters( + plant_config, {"electrolyzer": {"electricity_per_hydrogen": 50.0}} + ) prob = _build_hetero_problem( DemandFollowingControl, plant_config, slc_topology, - tech_config, + {}, demand=100.0, commodity_units={"electricity": "kW"}, ) diff --git a/h2integrate/control/control_strategies/system_level/test/test_slc_examples.py b/h2integrate/control/control_strategies/system_level/test/test_slc_examples.py index 0f7524fac..eac91a2e1 100644 --- a/h2integrate/control/control_strategies/system_level/test/test_slc_examples.py +++ b/h2integrate/control/control_strategies/system_level/test/test_slc_examples.py @@ -470,6 +470,6 @@ def test_slc_heterogeneous_commodity(subtests, temp_copy_of_example): ) with subtests.test("LCOA"): - assert pytest.approx(2.0024798555505585, rel=1e-5) == model.prob.get_val( + assert pytest.approx(2.00231927, rel=1e-5) == model.prob.get_val( "finance_subgroup_ammonia.LCOA", units="USD/kg" ) From b0c082da25ede9a8f21e437161de366c224ed483 Mon Sep 17 00:00:00 2001 From: John Jasa Date: Fri, 14 Aug 2026 10:40:58 -0600 Subject: [PATCH 5/5] Adding dynamics and figure plotting for hetereo example --- CHANGELOG.md | 1 + .../run_heterogeneous.py | 423 ++++++++++++++---- .../heterogeneous_commodity/tech_config.yaml | 9 +- 3 files changed, 352 insertions(+), 81 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 8ecf923e3..3d046e6ab 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,7 @@ - Add heterogeneous-commodity system-level control that translates demand for one commodity into upstream set-points across converters using static per-technology conversion ratios defined in the tech config. [PR TBD](https://github.com/NatLabRockies/H2Integrate/pull/TBD) - Extend the heterogeneous-commodity control to prefer measured conversion ratios computed from each converter's consumed and produced streams per timestep, falling back to the static ratio when a measurement is unavailable. [PR TBD](https://github.com/NatLabRockies/H2Integrate/pull/TBD) - Add the `35_system_level_control/heterogeneous_commodity` example, which serves an ammonia demand from a wind, battery, grid, electrolyzer, hydrogen-storage, and ammonia synthesis loop chain to demonstrate demand propagating from ammonia to hydrogen to electricity. The example uses profit-maximizing control so that wind always runs, the battery charges on wind surplus and discharges to cover deficits, and the grid is only dispatched to backfill the electricity that wind and the battery cannot supply, and it generates dispatch and dynamic conversion-ratio figures. [PR TBD](https://github.com/NatLabRockies/H2Integrate/pull/TBD) + - Add a system block diagram and a conversion-ratio chain figure to the heterogeneous-commodity example that show the technologies, their control classifications, the per-converter conversion ratios, and how those ratios multiply to propagate ammonia demand back into hydrogen and electricity demand. [PR TBD](https://github.com/NatLabRockies/H2Integrate/pull/TBD) ## 0.9 [August 10, 2026] diff --git a/examples/35_system_level_control/heterogeneous_commodity/run_heterogeneous.py b/examples/35_system_level_control/heterogeneous_commodity/run_heterogeneous.py index b2331df29..aca73dd63 100644 --- a/examples/35_system_level_control/heterogeneous_commodity/run_heterogeneous.py +++ b/examples/35_system_level_control/heterogeneous_commodity/run_heterogeneous.py @@ -1,6 +1,7 @@ import os from pathlib import Path +import yaml import numpy as np import matplotlib @@ -8,11 +9,70 @@ # Use a non-interactive backend so the figures render when the script is run headless. matplotlib.use("Agg") import matplotlib.pyplot as plt +from matplotlib.patches import Patch, FancyBboxPatch from h2integrate import EXAMPLE_DIR from h2integrate.core.h2integrate_model import H2IntegrateModel +# Technology fill colors keyed by control classification and arrow colors keyed by commodity. +# These are shared by the block diagrams so the same legend applies everywhere. +_CLASS_COLORS = { + "flexible": "#4C9F70", + "storage": "#F2C14E", + "dispatchable": "#9AA0A6", + "feedstock": "#6FA8DC", + "combiner": "#D9D9D9", + "demand": "#B39DDB", +} +_COMMODITY_COLORS = { + "electricity": "#3B7DD8", + "hydrogen": "#D8663B", + "ammonia": "#7B4FA3", + "nitrogen": "#4C9F70", +} + + +def _box(ax, xy, text, facecolor, width=12.0, height=8.0, fontsize=9, text_color="black"): + """Draw a rounded technology box centered at ``xy`` and return its center.""" + x, y = xy + patch = FancyBboxPatch( + (x - width / 2, y - height / 2), + width, + height, + boxstyle="round,pad=0.2,rounding_size=1.2", + linewidth=1.3, + edgecolor="#2F2F2F", + facecolor=facecolor, + zorder=3, + ) + ax.add_patch(patch) + ax.text(x, y, text, ha="center", va="center", fontsize=fontsize, color=text_color, zorder=5) + return x, y + + +def _flow_arrow(ax, start, end, color, label=None, lw=2.2, rad=0.0, label_dy=1.6): + """Draw a commodity flow arrow from ``start`` to ``end`` with an optional label.""" + ax.annotate( + "", + xy=end, + xytext=start, + arrowprops={ + "arrowstyle": "-|>", + "color": color, + "lw": lw, + "shrinkA": 6, + "shrinkB": 6, + "connectionstyle": f"arc3,rad={rad}", + }, + zorder=2, + ) + if label: + mx = (start[0] + end[0]) / 2 + my = (start[1] + end[1]) / 2 + label_dy + ax.text(mx, my, label, ha="center", va="center", fontsize=7.5, color=color, zorder=6) + + EXAMPLE_FOLDER = EXAMPLE_DIR / "35_system_level_control" / "heterogeneous_commodity" os.chdir(EXAMPLE_FOLDER) @@ -70,7 +130,7 @@ def make_plots(figure_dir): ammonia_demand = _get("ammonia_load_demand.ammonia_demand_out", "kg/h") # Controller set points that show backward propagation of demand. - ammonia_set_point = _get("system_level_controller.ammonia_ammonia_set_point", "kg/h") + _get("system_level_controller.ammonia_ammonia_set_point", "kg/h") hydrogen_set_point = _get("system_level_controller.electrolyzer_hydrogen_set_point", "kg/h") total_electricity_load_mw = electrolyzer_load_mw + synloop_load_mw @@ -132,72 +192,13 @@ def make_plots(figure_dir): plt.close(fig) # ----------------------------------------------------------------- - # Figure 2: backward demand propagation. A firm ammonia demand is - # translated into a firm hydrogen demand and then a firm electricity - # load. The demands are near-constant, so we show their magnitudes and - # annotate the conversion ratios that connect them. - # ----------------------------------------------------------------- - fig, ax_left = plt.subplots(figsize=(10, 5)) - ax_left.plot( - week_hours, ammonia_set_point[week], color="#7B4FA3", label="Ammonia set point (kg/h)" - ) - ax_left.plot( - week_hours, hydrogen_set_point[week], color="#D8663B", label="Hydrogen set point (kg/h)" - ) - ax_left.set_xlabel("Hour of representative week") - ax_left.set_ylabel("Commodity set point (kg/h)") - ax_left.set_ylim(0.0, 1.15 * float(ammonia_set_point[week].max())) - - ax_right = ax_left.twinx() - ax_right.plot( - week_hours, - total_electricity_load_mw[week], - color="#3B7DD8", - ls="--", - label="Derived electricity load (MW)", - ) - ax_right.set_ylabel("Electricity load (MW)") - ax_right.set_ylim(0.0, 1.15 * float(total_electricity_load_mw[week].max())) - - # Annotate the conversion ratios that link the three firm demand levels. - mean_hydrogen_per_ammonia = float(np.nanmean(hydrogen_per_ammonia)) - mean_electricity_per_ammonia_mw = float(total_electricity_load_mw.mean() / ammonia_out.mean()) - ax_left.annotate( - f"x {mean_hydrogen_per_ammonia:0.3f} kg H2 / kg NH3", - xy=(0.5, 0.62), - xycoords="axes fraction", - color="#D8663B", - fontsize=9, - ha="center", - ) - ax_left.annotate( - f"x {mean_electricity_per_ammonia_mw * 1000:0.2f} kWh / kg NH3 (total)", - xy=(0.5, 0.12), - xycoords="axes fraction", - color="#3B7DD8", - fontsize=9, - ha="center", - ) - - lines_left, labels_left = ax_left.get_legend_handles_labels() - lines_right, labels_right = ax_right.get_legend_handles_labels() - ax_left.legend( - lines_left + lines_right, labels_left + labels_right, loc="center right", fontsize=8 - ) - ax_left.set_title("Backward demand propagation: ammonia -> hydrogen -> electricity") - fig.tight_layout() - fig.savefig(figure_dir / "demand_propagation.png", dpi=150) - plt.close(fig) - - # ----------------------------------------------------------------- - # Figure 3: dynamic conversion ratios across the year. The controller + # Figure 2: dynamic conversion ratios across the year. The controller # prefers these measured ratios over the static seed values. The # electrolyzer ratio drifts up as the stack degrades over the year. # ----------------------------------------------------------------- fig, axes = plt.subplots(3, 1, figsize=(10, 9), sharex=True) axes[0].plot(hours, electricity_per_hydrogen, color="#3B7DD8", lw=0.8) - # axes[0].axhline(51.0, color="black", ls=":", lw=1.0, label="Static seed (51 kWh/kg)") axes[0].set_ylabel("kWh / kg H2") axes[0].set_title( "Electrolyzer measured ratio (electricity per hydrogen) drifts up with degradation" @@ -223,30 +224,298 @@ def make_plots(figure_dir): plt.close(fig) # ----------------------------------------------------------------- - # Figure 4: annual electricity source mix and the resulting LCOA. + # Figure 3: system block diagram. Technologies are colored by their + # control classification, arrows show the commodity flows, and the two + # converters (electrolyzer and ammonia synloop) carry the conversion + # ratios that the controller multiplies to propagate demand upstream. # ----------------------------------------------------------------- - wind_energy = wind_mw.sum() - grid_energy = grid_mw.sum() - battery_energy = np.clip(battery_mw, 0.0, None).sum() - lcoa = float(h2i.prob.get_val("finance_subgroup_ammonia.LCOA", units="USD/kg")[0]) - - fig, ax = plt.subplots(figsize=(6, 6)) - ax.pie( - [wind_energy, battery_energy, grid_energy], - labels=["Wind", "Battery", "Grid (firm)"], - colors=["#4C9F70", "#F2C14E", "#8C8C8C"], - autopct="%1.1f%%", - startangle=90, + # The static seed ratios now live in plant_config under + # system_level_control -> control_parameters -> conversion_parameters, + # so read them straight from there to keep the diagram in sync with the config. + with Path("plant_config.yaml").open() as f: + plant_cfg = yaml.safe_load(f) + conv = plant_cfg["system_level_control"]["control_parameters"]["conversion_parameters"] + seed_elec_per_h2 = float(conv["electrolyzer"]["electricity_per_hydrogen"]) + seed_h2_per_nh3 = float(conv["ammonia"]["hydrogen_per_ammonia"]) + seed_elec_per_nh3 = float(conv["ammonia"]["electricity_per_ammonia"]) + + # Measured (dynamic) ratios that the controller actually used, averaged over the year. + m_elec_per_h2 = float(np.nanmean(electricity_per_hydrogen)) + m_h2_per_nh3 = float(np.nanmean(hydrogen_per_ammonia)) + m_elec_per_nh3 = float(np.nanmean(electricity_per_ammonia)) + + # Representative demand magnitudes used to annotate the propagation band. + nh3_rate = float(ammonia_demand.mean()) + h2_rate = float(hydrogen_set_point.mean()) + electrolyzer_rate_mw = float(electrolyzer_load_mw.mean()) + synloop_rate_mw = float(synloop_load_mw.mean()) + + fig, ax = plt.subplots(figsize=(16, 9)) + ax.set_xlim(0, 100) + ax.set_ylim(0, 100) + ax.axis("off") + + color_class = _CLASS_COLORS + color_flow = _COMMODITY_COLORS + + # Technology boxes, laid out left to right along the forward commodity flow. + _box(ax, (9, 86), "Wind\n(flexible)\n120 MW", color_class["flexible"]) + _box(ax, (9, 68), "Battery\n(storage)\n10 MW / 40 MWh", color_class["storage"]) + _box(ax, (9, 50), "Grid\n(dispatchable)\n50 MW", color_class["dispatchable"]) + _box(ax, (26, 68), "Electricity\ncombiner", color_class["combiner"], width=11, height=7) + _box( + ax, + (44, 68), + "Electrolyzer\n(dispatchable)\nconverter: elec -> H2", + color_class["dispatchable"], + width=15, + height=11, + ) + _box(ax, (62, 88), "H2 storage\n(storage)", color_class["storage"], width=11, height=7) + _box(ax, (62, 68), "H2\ncombiner", color_class["combiner"], width=11, height=7) + _box(ax, (80, 90), "N2 feedstock\n(feedstock)", color_class["feedstock"], width=12, height=7) + _box( + ax, + (80, 68), + "Ammonia synloop\n(dispatchable)\nconverter: H2 -> NH3\n(+ direct electricity)", + color_class["dispatchable"], + width=16, + height=13, + ) + _box( + ax, + (80, 48), + "Electricity feedstock\n(feedstock)", + color_class["feedstock"], + width=13, + height=7, + ) + _box(ax, (95, 68), "Ammonia\ndemand\n4000 kg/h", color_class["demand"], width=10, height=9) + + # Electricity flows (blue). + _flow_arrow(ax, (14.5, 84), (20.5, 70), color_flow["electricity"]) + _flow_arrow(ax, (14.5, 68), (20.5, 68), color_flow["electricity"], label="electricity") + _flow_arrow(ax, (14.5, 52), (20.5, 66), color_flow["electricity"]) + _flow_arrow(ax, (31.5, 68), (36.5, 68), color_flow["electricity"], label="electricity") + _flow_arrow(ax, (80, 51.5), (80, 61.5), color_flow["electricity"], label="direct", label_dy=0) + _flow_arrow(ax, (5, 82), (5, 72), color_flow["electricity"], lw=1.6, label="charge", label_dy=0) + + # Hydrogen flows (orange). + _flow_arrow(ax, (51.5, 71), (56.5, 86), color_flow["hydrogen"], rad=0.25) + _flow_arrow(ax, (51.5, 68), (56.5, 68), color_flow["hydrogen"], label="hydrogen") + _flow_arrow(ax, (62, 84.5), (62, 71.5), color_flow["hydrogen"]) + _flow_arrow(ax, (67.5, 68), (72, 68), color_flow["hydrogen"], label="hydrogen") + + # Nitrogen feedstock and the final ammonia product. + _flow_arrow(ax, (80, 86.5), (80, 74.5), color_flow["nitrogen"], label="nitrogen", label_dy=0) + _flow_arrow(ax, (88, 68), (90, 68), color_flow["ammonia"], label="ammonia") + + # Conversion-ratio callouts on the two converters (seed and measured values). + ax.text( + 44, + 60.5, + f"x {seed_elec_per_h2:0.1f} kWh/kg H2 (seed)\nmeasured ~{m_elec_per_h2:0.1f}", + ha="center", + va="top", + fontsize=7.5, + color=color_flow["electricity"], + ) + ax.text( + 80, + 60.0, + f"x {seed_h2_per_nh3:0.2f} kg H2/kg NH3 (seed)\n" + f"+ {seed_elec_per_nh3:0.3f} kWh/kg NH3 direct", + ha="center", + va="top", + fontsize=7.5, + color=color_flow["hydrogen"], + ) + + # Backward demand-propagation band along the bottom. + band = FancyBboxPatch( + (4, 6), + 92, + 22, + boxstyle="round,pad=0.4,rounding_size=1.5", + linewidth=1.0, + edgecolor="#999999", + facecolor="#F5F5F5", + zorder=1, + ) + ax.add_patch(band) + ax.text( + 50, + 25.5, + "System-level controller: backward demand propagation " + "(measured consumed/produced ratios override the plant_config seeds)", + ha="center", + va="center", + fontsize=9.5, + weight="bold", + ) + _box(ax, (16, 15), f"NH3 demand\n{nh3_rate:,.0f} kg/h", "#EDE3F7", width=15, height=8) + _box(ax, (42, 15), f"H2 demand\n{h2_rate:,.0f} kg/h", "#FBE3D6", width=15, height=8) + _box( + ax, + (72, 19), + f"Electrolyzer load\n{electrolyzer_rate_mw:,.1f} MW", + "#DCE8F8", + width=17, + height=7, + ) + _box(ax, (72, 10), f"Synloop direct\n{synloop_rate_mw:,.1f} MW", "#DCE8F8", width=17, height=7) + _flow_arrow( + ax, + (23.5, 15), + (34.5, 15), + color_flow["hydrogen"], + label=f"x {seed_h2_per_nh3:0.2f}", + label_dy=1.4, + ) + _flow_arrow( + ax, + (49.5, 15), + (63.5, 19), + color_flow["electricity"], + label=f"x {seed_elec_per_h2:0.1f}", + label_dy=1.4, + ) + _flow_arrow( + ax, + (23.5, 13), + (63.5, 10), + color_flow["electricity"], + rad=0.1, + label=f"x {seed_elec_per_nh3:0.3f} (direct)", + label_dy=-1.8, + ) + + legend_handles = [ + Patch( + facecolor=color_class["flexible"], edgecolor="#2F2F2F", label="flexible (curtailable)" + ), + Patch(facecolor=color_class["storage"], edgecolor="#2F2F2F", label="storage"), + Patch(facecolor=color_class["dispatchable"], edgecolor="#2F2F2F", label="dispatchable"), + Patch(facecolor=color_class["feedstock"], edgecolor="#2F2F2F", label="feedstock"), + Patch(facecolor=color_class["combiner"], edgecolor="#2F2F2F", label="combiner"), + Patch(facecolor=color_class["demand"], edgecolor="#2F2F2F", label="demand"), + ] + ax.legend( + handles=legend_handles, + loc="upper center", + ncol=6, + fontsize=8, + frameon=False, + bbox_to_anchor=(0.5, 1.03), ) ax.set_title( - f"Annual electricity supplied to the chain\nAmmonia LCOA = ${lcoa:0.2f}/kg", fontsize=11 + "Heterogeneous-commodity system: technologies, conversion ratios, and demand propagation", + fontsize=13, + pad=22, + ) + fig.tight_layout() + fig.savefig(figure_dir / "system_block_diagram.png", dpi=150) + plt.close(fig) + + # ----------------------------------------------------------------- + # Figure 4: how the conversion ratios multiply along the chain to set + # the electricity intensity of one kilogram of ammonia. Seeds come from + # plant_config; the measured values are what the controller actually used. + # ----------------------------------------------------------------- + elec_via_electrolysis = m_h2_per_nh3 * m_elec_per_h2 # kWh electricity per kg NH3 through H2 + total_elec_per_nh3 = elec_via_electrolysis + m_elec_per_nh3 + + fig, (ax_chain, ax_bar) = plt.subplots( + 2, 1, figsize=(12, 8), gridspec_kw={"height_ratios": [3, 1]} + ) + ax_chain.set_xlim(0, 100) + ax_chain.set_ylim(0, 100) + ax_chain.axis("off") + + _box( + ax_chain, (12, 60), "1 kg\nammonia", color_class["demand"], width=14, height=16, fontsize=10 ) + _box(ax_chain, (45, 78), f"{m_h2_per_nh3:0.3f} kg\nhydrogen", "#FBE3D6", width=15, height=15) + _box( + ax_chain, + (80, 78), + f"{elec_via_electrolysis:0.2f} kWh\n(electrolysis)", + "#DCE8F8", + width=17, + height=15, + ) + _box( + ax_chain, + (80, 40), + f"{m_elec_per_nh3:0.3f} kWh\n(synloop direct)", + "#DCE8F8", + width=17, + height=15, + ) + + _flow_arrow( + ax_chain, + (19, 64), + (37.5, 78), + color_flow["hydrogen"], + label=f"x {m_h2_per_nh3:0.3f}\n(H2 per NH3)", + label_dy=3.5, + ) + _flow_arrow( + ax_chain, + (52.5, 78), + (71.5, 78), + color_flow["electricity"], + label=f"x {m_elec_per_h2:0.1f}\n(elec per H2)", + label_dy=3.5, + ) + _flow_arrow( + ax_chain, + (19, 56), + (71.5, 40), + color_flow["electricity"], + rad=0.1, + label=f"x {m_elec_per_nh3:0.3f} (direct elec per NH3)", + label_dy=-3.5, + ) + + ax_chain.text( + 50, + 10, + f"Total electricity intensity = {elec_via_electrolysis:0.2f} + {m_elec_per_nh3:0.3f}" + f" = {total_elec_per_nh3:0.2f} kWh per kg NH3", + ha="center", + va="center", + fontsize=9.5, + ) + ax_chain.set_title( + "Conversion ratios multiply along the chain: NH3 -> H2 -> electricity", fontsize=12 + ) + + ax_bar.barh( + [0], [elec_via_electrolysis], color=color_flow["electricity"], label="via electrolysis" + ) + ax_bar.barh( + [0], + [m_elec_per_nh3], + left=[elec_via_electrolysis], + color="#9AA0A6", + label="synloop direct", + ) + ax_bar.set_yticks([]) + ax_bar.set_xlabel("Electricity per kg ammonia (kWh/kg)") + ax_bar.legend(loc="lower right", fontsize=8, ncol=2) + ax_bar.set_xlim(0, total_elec_per_nh3 * 1.15) + ax_bar.text( + total_elec_per_nh3, 0, f" {total_elec_per_nh3:0.2f} kWh/kg", va="center", fontsize=9 + ) + fig.tight_layout() - fig.savefig(figure_dir / "electricity_source_mix.png", dpi=150) + fig.savefig(figure_dir / "conversion_ratio_chain.png", dpi=150) plt.close(fig) return figure_dir figures = make_plots(EXAMPLE_FOLDER / "outputs") -print(f"Saved dispatch and conversion-ratio figures to {figures}") +print(f"Saved dispatch, block-diagram, and conversion-ratio figures to {figures}") diff --git a/examples/35_system_level_control/heterogeneous_commodity/tech_config.yaml b/examples/35_system_level_control/heterogeneous_commodity/tech_config.yaml index 0274db901..60d37c10f 100644 --- a/examples/35_system_level_control/heterogeneous_commodity/tech_config.yaml +++ b/examples/35_system_level_control/heterogeneous_commodity/tech_config.yaml @@ -79,6 +79,8 @@ technologies: uptime_hours_until_eol: 80000. # number of 'on' hours until electrolyzer reaches eol include_degradation_penalty: true # include degradation turndown_ratio: 0.1 # turndown_ratio = minimum_cluster_power/cluster_rating_MW + ramp_up_rate_fraction: 0.2 # fraction of electrolyzer rating per timestep + ramp_down_rate_fraction: 0.2 # fraction of electrolyzer rating per timestep financial_parameters: capital_items: depr_period: 7 @@ -187,10 +189,9 @@ technologies: purge_gas_x_ar: 0.02 purge_gas_x_nh3: 0.04 purge_gas_mass_ratio: 0.07 - # Dynamics disabled for a well-behaved system-level control fixed-point loop. - turndown_ratio: 0.0 - ramp_up_rate_fraction: 1.0 - ramp_down_rate_fraction: 1.0 + turndown_ratio: 0.2 + ramp_up_rate_fraction: 0.2 + ramp_down_rate_fraction: 0.2 include_cold_start: false include_warm_start: false cost_parameters: