Coverage for backend/ahuora-builder/src/ahuora_builder/custom/thermal_utility_systems/header.py: 83%
328 statements
« prev ^ index » next coverage.py v7.10.7, created at 2026-07-22 05:22 +0000
1"""Header unit model.
2Implementation by Tim and Keegan, taken from Ahuora-UnitOperations repo"""
4from typing import List, Optional, Iterable
6from pyomo.common.config import Bool, ConfigBlock, ConfigValue, In
7import pyomo.environ as pyo
8from pyomo.environ import (
9 Constraint,
10 Expression,
11 Param,
12 PositiveReals,
13 RangeSet,
14 Suffix,
15 Var,
16 check_optimal_termination,
17 value,
18 units as pyunits,
19)
20from pyomo.core.base.reference import Reference
22from idaes.core import StateBlock, UnitModelBlockData, declare_process_block_class, useDefault
23from idaes.core.initialization import ModularInitializerBase
24from idaes.core.solvers import get_solver
25from idaes.core.util import scaling as iscale
26from idaes.core.scaling import CustomScalerBase, ConstraintScalingScheme
27from idaes.core.util.config import is_physical_parameter_block
28from idaes.core.util.exceptions import InitializationError
29from idaes.core.util.math import smooth_max, smooth_min
30from idaes.core.util.model_diagnostics import DiagnosticsToolbox
31from idaes.core.util.model_statistics import degrees_of_freedom, report_statistics
32from idaes.core.util.tables import create_stream_table_dataframe
34import idaes.logger as idaeslog
36_log = idaeslog.getLogger(__name__)
38__author__ = "Ahuora Centre for Smart Energy Systems, University of Waikato, New Zealand"
41def _build_config(config: ConfigBlock) -> None:
42 """Declare configuration options for the Header unit.
44 Declares property package references and integer counts for inlets and outlets.
46 Args:
47 config (ConfigBlock): The mutable configuration block to populate.
49 Raises:
50 ValueError: If invalid option values are provided by the caller (via IDAES).
51 """
53 config.declare(
54 "property_package",
55 ConfigValue(
56 default=useDefault,
57 domain=is_physical_parameter_block,
58 description="Property package to use for control volume",
59 ),
60 )
61 config.declare(
62 "property_package_args",
63 ConfigBlock(
64 implicit=True,
65 description="Arguments to use for constructing property packages",
66 ),
67 )
68 config.declare(
69 "num_inlets",
70 ConfigValue(
71 default=1,
72 domain=In(list(range(1, 100))),
73 description="Number of utility providers at inlets.",
74 ),
75 )
76 config.declare(
77 "num_outlets",
78 ConfigValue(
79 default=1,
80 domain=In(list(range(0, 100))),
81 description=(
82 "Number of utility users at outlets. Excludes outlets "
83 "associated with condensate and vent flows."
84 ),
85 ),
86 )
87 config.declare(
88 "is_liquid_header",
89 ConfigValue(
90 default=False,
91 domain=Bool,
92 description="Flag for selecting liquid or vapour (including steam and other gases).",
93 ),
94 )
97@declare_process_block_class("simple_header")
98class SimpleHeaderData(UnitModelBlockData):
99 """Thermal utility header unit operation."""
101 CONFIG = UnitModelBlockData.CONFIG()
102 _build_config(CONFIG)
104 def build(self) -> None:
105 """Build the unit model structure and equations."""
106 super().build()
107 units_meta = self.config.property_package.get_metadata().get_derived_units
109 if self.config.num_inlets < 1: 109 ↛ 110line 109 didn't jump to line 110 because the condition on line 109 was never true
110 raise ValueError("Header requires at least one provider (num_inlets >= 1).")
111 if self.config.num_outlets < 1: 111 ↛ 112line 111 didn't jump to line 112 because the condition on line 111 was never true
112 raise ValueError("Header requires at least one user (num_outlets >= 1).")
114 self.inlet_list = [f"inlet_{i+1}" for i in range(self.config.num_inlets)]
115 self.outlet_list = [f"outlet_{i+1}" for i in range(self.config.num_outlets)] + [
116 "outlet_condensate",
117 "outlet_vent",
118 ]
120 """
121 1. Build inlet, outlet and internal state blocks and associate
122 with ports (where applicable)
123 """
124 self.inlet_blocks = self._build_state_blocks(
125 stream_name_list=self.inlet_list,
126 has_phase_equilibrium=False,
127 is_defined_state=True,
128 is_build_port=True,
129 )
132 self._outlet_supply_blocks = self._build_state_blocks(
133 stream_name_list=[f"outlet_{i+1}" for i in range(self.config.num_outlets)],
134 has_phase_equilibrium=False,
135 is_defined_state=False,
136 is_build_port=True,
137 )
138 self._outlet_vent_blocks = self._build_state_blocks(
139 stream_name_list=["outlet_condensate", "outlet_vent"],
140 has_phase_equilibrium=False,
141 is_defined_state=False,
142 is_build_port=True,
143 )
144 self.outlet_blocks = self._outlet_supply_blocks + self._outlet_vent_blocks
145 self.internal_blocks = self._build_state_blocks(
146 stream_name_list=["mixed"],
147 has_phase_equilibrium=True,
148 is_defined_state=False,
149 is_build_port=False,
150 )
151 for sb in (self.inlet_blocks + self.internal_blocks + self.outlet_blocks):
152 for t in sb:
153 sb[t].flow_mol.setlb(0.0)
155 """
156 2. Create parameters, variables, references and expressions
157 """
158 # References
159 self.total_flow_mol = Reference(self.mixed_state[:].flow_mol)
160 self.total_flow_mass = Reference(self.mixed_state[:].flow_mass)
161 self.pressure = Reference(self.mixed_state[:].pressure)
162 self.temperature = Reference(self.mixed_state[:].temperature)
163 self.enth_mol = Reference(self.mixed_state[:].enth_mol)
164 self.enth_mass = Reference(self.mixed_state[:].enth_mass)
165 self.vapor_frac = Reference(self.mixed_state[:].vapor_frac)
167 # State variables
168 self.heat_loss = Var(
169 self.flowsheet().time,
170 initialize=0.0,
171 bounds=(0, None),
172 units=units_meta("power"),
173 doc="Heat loss",
174 )
175 self.pressure_loss = Var(
176 self.flowsheet().time,
177 initialize=0,
178 bounds=(0, None),
179 units=units_meta("pressure"),
180 doc="Pressure loss.",
181 )
183 # Non-normal state variables, other parameters and expression
184 if self.config.is_liquid_header:
185 self._liq_out_enth_mol = Var(
186 self.flowsheet().time,
187 initialize=42.0 * 18 * pyunits.J / pyunits.mol,
188 units=units_meta("energy") / units_meta("amount"),
189 doc="Molar enthalpy of the liquid outlets",
190 )
191 else:
192 self._vap_out_enth_mol = Var(
193 self.flowsheet().time,
194 initialize=2700.0 * 18 * pyunits.J / pyunits.mol,
195 units=units_meta("energy") / units_meta("amount"),
196 doc="Molar enthalpy of the vapour outlets",
197 )
198 inlet_idx = RangeSet(len(self.inlet_blocks))
199 self._minimum_pressure = Var(
200 self.flowsheet().time,
201 inlet_idx,
202 doc="Variable for calculating minimum inlet pressure",
203 units=units_meta("pressure"),
204 )
205 self._eps_pressure = Param(
206 mutable=True,
207 initialize=1e-3,
208 domain=PositiveReals,
209 doc="Smoothing term for minimum inlet pressure",
210 units=units_meta("pressure"),
211 )
213 # Expressions
214 self.degree_of_superheat = Expression(
215 self.flowsheet().time,
216 rule=lambda b, t: (
217 b.temperature[t] - b.outlet_condensate_state[t].temperature
218 ),
219 )
220 self._partial_total_flow_mol = Expression(
221 self.flowsheet().time,
222 rule=lambda b, t: (
223 sum(o[t].flow_mol for o in (b.inlet_blocks + b._outlet_supply_blocks))
224 ),
225 )
226 self.balance_flow_mol = Expression(
227 self.flowsheet().time,
228 rule=lambda b, t: (
229 sum(i[t].flow_mol for i in b.inlet_blocks)
230 -
231 sum(
232 o[t].flow_mol
233 for o in (
234 b._outlet_supply_blocks
235 +
236 [
237 b.outlet_vent_state
238 if self.config.is_liquid_header
239 else b.outlet_condensate_state,
240 ]
241 )
242 )
243 ),
244 doc="Flow imbalance between inlets and outlets; positive if in excess of supply to outlets.",
245 )
246 self.makeup_flow_mol = Expression(
247 self.flowsheet().time,
248 rule=lambda b, t: (
249 (
250 b.outlet_condensate_state[t].flow_mol
251 if self.config.is_liquid_header
252 else b.outlet_vent_state[t].flow_mol
253 )
254 -
255 b.balance_flow_mol[t]
256 ),
257 )
259 """
260 3. Declare constraints to define mass, energy, and momentum balances,
261 unit operation performance and other constraint
262 """
263 # a) Material balance equations
264 @self.Constraint(self.flowsheet().time, doc="Mixed state material balance")
265 def mixed_state_material_balance(b, t):
266 return b.mixed_state[t].flow_mol == sum(i[t].flow_mol for i in b.inlet_blocks)
268 eps_smooth = 1e-5 # smoothing parameter; smaller = closer to exact max, larger = smoother
269 eps_div = 1e-6 # small number to prevent division by zero
270 if self.config.is_liquid_header:
271 # Assigns excess liquid flow to outlet_condensate
272 @self.Constraint(self.flowsheet().time, doc="Condensate flow balance.")
273 def condensate_flow_balance(b, t):
274 return (
275 b.outlet_condensate_state[t].flow_mol
276 ==
277 smooth_max(
278 b.balance_flow_mol[t] / (b._partial_total_flow_mol[t] + eps_div*pyunits.mol / pyunits.s),
279 0.0,
280 eps_smooth,
281 ) * (b._partial_total_flow_mol[t] + eps_div*pyunits.mol / pyunits.s)
282 )
284 # Removes any gas/vapour from a liquid header
285 @self.Constraint(self.flowsheet().time, doc="Vent flow balance.")
286 def vent_flow_balance(b, t):
287 return b.outlet_vent_state[t].flow_mol == b.mixed_state[t].flow_mol * b.mixed_state[t].vapor_frac
289 else:
290 # Assigns excess steam/vapour flow to outlet_vent
291 @self.Constraint(self.flowsheet().time, doc="Vent flow balance.")
292 def vent_flow_balance(b, t):
293 return (
294 b.outlet_vent_state[t].flow_mol
295 ==
296 smooth_max(
297 b.balance_flow_mol[t] / (b._partial_total_flow_mol[t] + eps_div*pyunits.mol / pyunits.s),
298 0.0,
299 eps_smooth,
300 ) * (b._partial_total_flow_mol[t] + eps_div*pyunits.mol / pyunits.s)
301 )
303 # Removes any condensate/liquid from a steam/gas header
304 @self.Constraint(self.flowsheet().time, doc="Condensate flow balance.")
305 def condensate_flow_balance(b, t):
306 return b.outlet_condensate_state[t].flow_mol == b.mixed_state[t].flow_mol * (1 - b.mixed_state[t].vapor_frac)
308 # b) Energy balance equations
309 @self.Constraint(self.flowsheet().time, doc="Energy balance for inlets to mixed state including heat loss")
310 def inlets_to_mixed_state_energy_balance(b, t):
311 return (
312 b.mixed_state[t].flow_mol * b.mixed_state[t].enth_mol
313 + b.heat_loss[t]
314 == sum(i[t].flow_mol * i[t].enth_mol for i in b.inlet_blocks)
315 )
317 @self.Constraint(self.flowsheet().time, doc="Energy balance for mixed state to outlets")
318 def mixed_state_to_outlets_energy_balance(b, t):
319 return (
320 b.mixed_state[t].enth_mol * sum(o[t].flow_mol for o in b.outlet_blocks)
321 ==
322 sum(o[t].flow_mol * o[t].enth_mol for o in b.outlet_blocks)
323 )
325 if self.config.is_liquid_header:
326 self._outlet_blocks_exc_vent = self._outlet_supply_blocks + [self.outlet_condensate_state]
327 @self.Constraint(self.flowsheet().time, self._outlet_blocks_exc_vent,
328 doc="All liquid outlets (incl. condensate) share a common liquid enthalpy",
329 )
330 def molar_enthalpy_equality_eqn(b, t, o):
331 return o[t].enth_mol == b._liq_out_enth_mol[t]
332 else:
333 self._outlet_blocks_exc_condensate = self._outlet_supply_blocks + [self.outlet_vent_state]
334 @self.Constraint(self.flowsheet().time, self._outlet_blocks_exc_condensate,
335 doc="All vapour outlets (incl. vent) share a common vapour enthalpy",
336 )
337 def molar_enthalpy_equality_eqn(b, t, o):
338 return o[t].enth_mol == b._vap_out_enth_mol[t]
340 # c) Momentum balance equations
341 @self.Constraint(self.flowsheet().time, inlet_idx,
342 doc="Calculation for minimum inlet pressure",
343 )
344 def minimum_pressure_constraint(b, t, i):
345 if i == inlet_idx.first():
346 return b._minimum_pressure[t, i] == b.inlet_blocks[i - 1][t].pressure
347 else:
348 return (
349 b._minimum_pressure[t, i]
350 ==
351 smooth_min(
352 b._minimum_pressure[t, i - 1],
353 b.inlet_blocks[i - 1][t].pressure,
354 b._eps_pressure,
355 )
356 )
357 # Set mixed pressure to minimum inlet pressure minus any pressure loss
358 @self.Constraint(self.flowsheet().time, doc="Pressure equality constraint from minimum inlet to mixed state")
359 def mixture_pressure(b, t):
360 return b.mixed_state[t].pressure == (
361 b._minimum_pressure[t, inlet_idx.last()] - b.pressure_loss[t]
362 )
363 # Set outlet pressures to mixed pressure
364 @self.Constraint(self.flowsheet().time, self.outlet_blocks,
365 doc="Pressure equality constraint from mixed state to outlets",
366 )
367 def pressure_equality_eqn(b, t, o):
368 return b.mixed_state[t].pressure == o[t].pressure
370 # d) Additional constraints for vapour/liquid separation and outlet splits
371 if self.config.is_liquid_header:
372 @self.Constraint(self.flowsheet().time, doc="Vent vapour fraction.")
373 def vent_vapour_fraction(b, t):
374 return b.outlet_vent_state[t].vapor_frac == 1 # 1 - 1e-6
375 else:
376 @self.Constraint(self.flowsheet().time, doc="Condensate vapour fraction.")
377 def condensate_vapour_fraction(b, t):
378 return b.outlet_condensate_state[t].vapor_frac == 0 # 1e-6
380 """
381 4. Other model components and references
382 """
383 self.scaling_factor = Suffix(direction=Suffix.EXPORT)
384 self.outlet_idx = pyo.Set(initialize=self.outlet_list)
385 # Map each (t, o) to the outlet state's flow var
386 ref_map = {}
387 for o in self.outlet_list:
388 if o != "vent": 388 ↛ 387line 388 didn't jump to line 387 because the condition on line 388 was always true
389 outlet_state_block = getattr(self, f"{o}_state")
390 for t in self.flowsheet().time:
391 ref_map[(t, o)] = outlet_state_block[t].flow_mol
393 self.split_flow = Reference(ref_map)
396 def initialize_build(self, state_args=None, outlvl=idaeslog.NOTSET, solver=None, optarg=None):
397 """
398 General wrapper for template initialization routines
400 Keyword Arguments:
401 state_args : a dict of arguments to be passed to the property
402 package(s) to provide an initial state for
403 initialization (see documentation of the specific
404 property package) (default = {}).
405 outlvl : sets output level of initialization routine
406 optarg : solver options dictionary object (default=None)
407 solver : str indicating which solver to use during
408 initialization (default = None)
410 Returns: None
411 """
412 init_log = idaeslog.getInitLogger(self.name, outlvl, tag="unit")
413 solve_log = idaeslog.getSolveLogger(self.name, outlvl, tag="unit")
414 t0 = self.flowsheet().config.time.first()
416 opt = get_solver(solver, optarg)
418 pp = self.mixed_state[t0].params
419 state_args = {} if state_args is None else dict(state_args)
421 shared_state_args = {
422 key: state_args[key]
423 for key in ("flow_mol", "pressure", "enth_mol", "temperature")
424 if key in state_args
425 }
427 def _state_seed_args(stream_name):
428 local_args = dict(shared_state_args)
429 local_args.update(state_args.get(stream_name, {}))
430 return local_args
432 inlet_state_args = {name: _state_seed_args(name) for name in self.inlet_list}
433 outlet_state_args = {name: _state_seed_args(name) for name in self.outlet_list}
434 mixed_state_args = _state_seed_args("mixed")
436 # Helper functions
437 def _value_or_none(obj): # returns the value of a Var or Param if it is fixed, otherwise returns None
438 return value(obj, exception=False)
440 def _pick_seed(*candidates): # loops through different options of seeding and picks the first (in order of preference)
441 for candidate in candidates: 441 ↛ 444line 441 didn't jump to line 444 because the loop on line 441 didn't complete
442 if candidate is not None:
443 return candidate
444 return None
446 def _enthalpy_from_tp(temperature, pressure):
447 if temperature is None or pressure is None: 447 ↛ 448line 447 didn't jump to line 448 because the condition on line 447 was never true
448 return None
449 return value(pp.htpx(temperature * pyunits.K, pressure * pyunits.Pa))
452 def _safe_nonnegative(val, fallback=0.0):
453 if val is None:
454 return fallback
455 return max(val, 0.0)
457 #-----------------------------------------------------
458 # 1. Build state seeds. Prefer explicit state_args, then connected/fixed
459 # values already present on the unit, then infer from local specs, then
460 # fall back to generic guesses.
462 # Inlets flows ranking
463 # 1. Explicit state_args
464 # 2. A fixed inlet flow on this unit
465 # 3. A propagated upstream inlet flow
466 # 4. Back-calculate from fixed outlet flow(s) + 10% divided evenly among inlets. The 10% is for venting
467 # 5. Nominal fallback
469 # First get an estimate of the total inlet based on outlet flows if they were fixed
470 fixed_outlet_flows = []
471 for outlet_sb in self.outlet_blocks:
472 outlet_flow_fixed = (
473 _value_or_none(outlet_sb[t0].flow_mol)
474 if outlet_sb[t0].flow_mol.fixed
475 else None
476 )
477 if outlet_flow_fixed is not None:
478 fixed_outlet_flows.append(outlet_flow_fixed)
480 total_fixed_outlet_flow = (
481 sum(fixed_outlet_flows) if fixed_outlet_flows else None
482 )
483 nominal_total_flow = _pick_seed(total_fixed_outlet_flow, 500.0)
484 nominal_inlet_flow = nominal_total_flow / max(len(self.inlet_blocks), 1)
486 # Now build the inlet seeds using the ranking above
487 inlet_seeds = {}
488 seeded_inlet_flows = []
489 for inlet_name, inlet_sb in zip(self.inlet_list, self.inlet_blocks):
490 inlet_has_source = len(list(getattr(self, inlet_name).sources())) > 0
491 local_args = inlet_state_args[inlet_name]
493 flow_from_fixed = (
494 _value_or_none(inlet_sb[t0].flow_mol)
495 if inlet_sb[t0].flow_mol.fixed
496 else None
497 )
498 flow_from_upstream = (
499 _value_or_none(inlet_sb[t0].flow_mol)
500 if inlet_has_source
501 else None
502 )
503 flow_from_fixed_outlets = (
504 total_fixed_outlet_flow / max(len(self.inlet_blocks), 1) * 1.1
505 if total_fixed_outlet_flow is not None
506 else None
507 )
509 f_in = _pick_seed(
510 local_args.get("flow_mol"),
511 flow_from_fixed,
512 flow_from_upstream,
513 flow_from_fixed_outlets,
514 nominal_inlet_flow,
515 )
517 inlet_seeds[inlet_name] = {"flow_mol": f_in}
518 seeded_inlet_flows.append(f_in)
520 # Pass these seeded inlet flows to the mixed state and outlet states
521 f_mixed = sum(seeded_inlet_flows)
523 # Outlet flow ranking:
524 # 1. Fixed outlet flow on this unit
525 # 2. Even split of seeded inlet flow (minus 10% for venting) divided among numbered outlets, vent gets 10%
527 # First split the total flow into outlet and vent portions. The outlet portion is then divided evenly among the numbered outlets, and the vent portion is assigned to the vent outlet. This is just a starting guess to help initialization, and will be overridden if there are fixed outlet flows that provide a better seed.
528 numbered_outlet_flow = 0.9 * f_mixed / max(self.config.num_outlets, 1)
529 vent_flow = 0.1 * f_mixed
531 # Now build the outlet seeds using the ranking above
532 outlet_seeds = {}
533 seeded_outlet_flows = []
534 for outlet_name, outlet_sb in zip(self.outlet_list, self.outlet_blocks):
535 flow_from_fixed = (
536 _value_or_none(outlet_sb[t0].flow_mol)
537 if outlet_sb[t0].flow_mol.fixed
538 else None
539 )
541 if outlet_name.startswith("outlet_") and outlet_name not in (
542 "outlet_condensate",
543 "outlet_vent",
544 ):
545 default_flow = numbered_outlet_flow
546 elif outlet_name == "outlet_vent":
547 default_flow = vent_flow
548 else:
549 default_flow = 0.0
551 f_out = _pick_seed(
552 flow_from_fixed,
553 default_flow,
554 )
556 outlet_seeds[outlet_name] = {"flow_mol": f_out}
557 seeded_outlet_flows.append(f_out)
559 # Inlet pressure ranking
560 # 1. Explicit state_args
561 # 2. A propagated upstream inlet pressure
562 # 3. A fixed inlet pressure on this unit
563 # 4. Assume all inlet pressures equal to the outlet pressure + pressure loss
564 # 5. Nominal fallback of 10 bar
566 # First determine the outlet pressure if it were fixed
567 pressure_loss = _pick_seed(_value_or_none(self.pressure_loss[t0]), 0.0)
568 fixed_outlet_pressure = None
569 for outlet_sb in self.outlet_blocks:
570 outlet_pressure_fixed = (
571 _value_or_none(outlet_sb[t0].pressure)
572 if outlet_sb[t0].pressure.fixed
573 else None
574 )
575 if outlet_pressure_fixed is not None: 575 ↛ 576line 575 didn't jump to line 576 because the condition on line 575 was never true
576 fixed_outlet_pressure = outlet_pressure_fixed
577 break
579 inlet_pressure_from_outlet = (
580 fixed_outlet_pressure + pressure_loss
581 if fixed_outlet_pressure is not None
582 else None
583 )
585 # Now build the inlet pressure seeds using the ranking above
586 seeded_inlet_pressures = []
587 for inlet_name, inlet_sb in zip(self.inlet_list, self.inlet_blocks):
588 inlet_has_source = len(list(getattr(self, inlet_name).sources())) > 0
589 local_args = inlet_state_args[inlet_name]
591 pressure_from_upstream = (
592 _value_or_none(inlet_sb[t0].pressure)
593 if inlet_has_source
594 else None
595 )
596 pressure_from_fixed = (
597 _value_or_none(inlet_sb[t0].pressure)
598 if inlet_sb[t0].pressure.fixed
599 else None
600 )
602 p_in = _pick_seed(
603 local_args.get("pressure"),
604 pressure_from_upstream,
605 pressure_from_fixed,
606 inlet_pressure_from_outlet,
607 10e5,
608 )
610 inlet_seeds[inlet_name]["pressure"] = p_in
611 seeded_inlet_pressures.append(p_in)
613 # Pass the seeded mixed pressure to the outlet states
614 p_mixed = min(seeded_inlet_pressures) - pressure_loss
615 for outlet_name, outlet_sb in zip(self.outlet_list, self.outlet_blocks):
616 outlet_seeds[outlet_name]["pressure"] = p_mixed
618 # Inlet enthalpy ranking
619 # 1. Explicit state_args
620 # 2. A propagated upstream inlet enthalpy
621 # 3. Fixed inlet enthalpy on this unit
622 # 4. Back calculate from fixed outlet enthalpies and seeded inlet flows with user heat loss
623 # 5. Assume inlet enthalpies all equal to the saturation enthalpy + 10 C superheat # TODO: adapt this for water headers
625 # For rank 4, first determine the outlet enthalpies if they were fixed, and calculate the inlet enthalpy with heat loss
626 heat_loss = _pick_seed(_value_or_none(self.heat_loss[t0]), 0.0)
627 inlet_enthalpy_from_outlets = None
628 if sum(seeded_inlet_flows) > 0:
629 outlet_energy_terms = []
630 all_outlet_enthalpies_fixed = True
631 for outlet_name, outlet_sb in zip(self.outlet_list, self.outlet_blocks): 631 ↛ 647line 631 didn't jump to line 647 because the loop on line 631 didn't complete
632 outlet_enthalpy_fixed = (
633 _value_or_none(outlet_sb[t0].enth_mol)
634 if outlet_sb[t0].enth_mol.fixed
635 else None
636 )
637 seeded_outlet_flow = outlet_seeds[outlet_name].get("flow_mol")
639 if outlet_enthalpy_fixed is None or seeded_outlet_flow is None: 639 ↛ 643line 639 didn't jump to line 643 because the condition on line 639 was always true
640 all_outlet_enthalpies_fixed = False
641 break
643 outlet_energy_terms.append(
644 seeded_outlet_flow * outlet_enthalpy_fixed
645 )
647 if all_outlet_enthalpies_fixed: 647 ↛ 648line 647 didn't jump to line 648 because the condition on line 647 was never true
648 inlet_enthalpy_from_outlets = (
649 sum(outlet_energy_terms) + heat_loss
650 ) / (sum(seeded_inlet_flows) + 1e-6)
652 seeded_inlet_enthalpies = []
653 for inlet_name, inlet_sb in zip(self.inlet_list, self.inlet_blocks):
654 inlet_has_source = len(list(getattr(self, inlet_name).sources())) > 0
655 local_args = inlet_state_args[inlet_name]
656 p_in = inlet_seeds[inlet_name]["pressure"]
657 flow_mol = inlet_seeds[inlet_name]["flow_mol"]
658 enthalpy_from_upstream = (
659 _value_or_none(inlet_sb[t0].enth_mol)
660 if inlet_has_source
661 else None
662 )
663 enthalpy_from_fixed = (
664 _value_or_none(inlet_sb[t0].enth_mol)
665 if inlet_sb[t0].enth_mol.fixed
666 else None
667 )
669 inlet_sb[t0].pressure.set_value(p_in)
670 temperature_sat = _value_or_none(
671 getattr(inlet_sb[t0], "temperature_sat", None)
672 )
674 h_in = _pick_seed(
675 local_args.get("enth_mol"),
676 enthalpy_from_upstream,
677 enthalpy_from_fixed,
678 inlet_enthalpy_from_outlets,
679 _enthalpy_from_tp(
680 temperature_sat + 10.0,
681 p_in,
682 ) if temperature_sat is not None else None,
683 )
685 inlet_seeds[inlet_name]["enth_mol"] = h_in
686 seeded_inlet_enthalpies.append(h_in*flow_mol)
688 # Seed inlet enthalpy to mixed and outlet states
689 h_mixed = (sum(seeded_inlet_enthalpies) - pyo.value(self.heat_loss[t0])) / (sum(seeded_inlet_flows) + 1e-6)
690 for outlet_name, outlet_sb in zip(self.outlet_list, self.outlet_blocks):
691 outlet_seeds[outlet_name]["enth_mol"] = h_mixed
693 # TODO: do better seeding of condensate enthalpy
694 #-----------------------------------------------------
695 # 2. Initialize state blocks using explicitly seeded values for all state variables
696 for inlet_name, inlet_sb in zip(self.inlet_list, self.inlet_blocks):
697 inlet_sb.initialize(
698 solver=solver,
699 optarg=optarg,
700 outlvl=outlvl,
701 state_args=inlet_seeds[inlet_name],
702 )
703 init_log.info_high(f"{inlet_name} state initialization complete")
705 flags_mixed = self.mixed_state.initialize(
706 solver=solver,
707 optarg=optarg,
708 outlvl=outlvl,
709 state_args={
710 "flow_mol": _pick_seed(mixed_state_args.get("flow_mol"), f_mixed),
711 "pressure": _pick_seed(mixed_state_args.get("pressure"), p_mixed),
712 "enth_mol": _pick_seed(mixed_state_args.get("enth_mol"), h_mixed),
713 },
714 hold_state=True,
715 )
716 init_log.info_high("mixed state initialization complete")
718 held_outlet_states = {}
719 for outlet_name, outlet_sb in zip(self.outlet_list, self.outlet_blocks):
720 held_outlet_states[outlet_name] = outlet_sb.initialize(
721 solver=solver,
722 optarg=optarg,
723 outlvl=outlvl,
724 state_args=outlet_seeds[outlet_name],
725 hold_state=True,
726 )
727 init_log.info_high(f"{outlet_name} state initialization complete")
728 # #-----------------------------------------------------
729 # # 3. Deactivate constraints not specified in Step 2 then solve first pass
730 relaxed_eqns = [
731 self.mixed_state_material_balance,
732 self.condensate_flow_balance,
733 self.vent_flow_balance,
734 self.inlets_to_mixed_state_energy_balance,
735 self.mixed_state_to_outlets_energy_balance,
736 self.molar_enthalpy_equality_eqn,
737 self.minimum_pressure_constraint,
738 self.mixture_pressure,
739 self.pressure_equality_eqn,
740 ]
741 if self.config.is_liquid_header:
742 relaxed_eqns.append(self.vent_vapour_fraction)
743 else:
744 relaxed_eqns.append(self.condensate_vapour_fraction)
746 for con in relaxed_eqns:
747 con.deactivate()
749 report_statistics(self)
751 dt = DiagnosticsToolbox(self)
752 dt.report_structural_issues()
753 dt.display_underconstrained_set()
755 active_constraints = sum(
756 1 for _ in self.component_data_objects(Constraint, active=True, descend_into=True)
757 )
759 if active_constraints > 0: 759 ↛ 760line 759 didn't jump to line 760 because the condition on line 759 was never true
760 with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
761 res = opt.solve(self, tee=slc.tee)
762 if not check_optimal_termination(res):
763 dt.report_numerical_issues()
764 raise InitializationError(f"{self.name} failed relaxed initialization")
765 else:
766 init_log.info_high(
767 "Relaxed initialization pass skipped: no active constraints remained after seeding."
768 )
770 # Restore full model
771 self.mixed_state.release_state(flags_mixed)
772 for outlet_name, flags in held_outlet_states.items():
773 getattr(self, f"{outlet_name}_state").release_state(flags)
775 for con in relaxed_eqns:
776 con.activate()
778 dof = degrees_of_freedom(self)
779 # NOTE: Seems header mostly has >0 DoF in tests so check is skipped, but should be revisited to debug it
780 # if dof != 0:
781 # raise InitializationError(
782 # f"{self.name} degrees of freedom were not 0 before final solve. DoF = {dof}"
783 # )
785 with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
786 res = opt.solve(self, tee=slc.tee)
787 if not check_optimal_termination(res): 787 ↛ 788line 787 didn't jump to line 788 because the condition on line 787 was never true
788 raise InitializationError(f"{self.name} failed final initialization")
790 init_log.info(f"Initialization complete: {idaeslog.condition(res)}")
792 def _get_performance_contents(self, time_point=0, is_full_report=True):
793 """Collect performance results for reporting.
795 Args:
796 time_point (int | float): Time index at which to report values.
797 is_full_report (bool): Flag for full or partial performance report.
799 Returns:
800 dict: A report of internal unit model results.
801 """
802 if is_full_report:
803 var_dict = {
804 "Heat Loss": self.heat_loss[time_point],
805 "Pressure Drop": self.pressure_loss[time_point],
806 "Mass Flow": self.mixed_state[time_point].flow_mass,
807 "Molar Flow": self.mixed_state[time_point].flow_mol,
808 "Balance Flow": self.balance_flow_mol[time_point],
809 "Pressure": self.mixed_state[time_point].pressure,
810 "Temperature": self.mixed_state[time_point].temperature,
811 "Degree of Superheat": self.degree_of_superheat[time_point],
812 "Vapour Fraction": self.mixed_state[time_point].vapor_frac,
813 "Mass Specific Enthalpy": self.mixed_state[time_point].enth_mass,
814 "Molar Specific Enthalpy": self.mixed_state[time_point].enth_mol,
815 }
816 else:
817 var_dict = {
818 "Balance Flow": self.balance_flow_mol[time_point],
819 "Pressure": self.mixed_state[time_point].pressure,
820 "Temperature": self.mixed_state[time_point].temperature,
821 "Degree of Superheat": self.degree_of_superheat[time_point],
822 }
824 return {"vars": var_dict}
826 def diagnose(self) -> list[tuple[pyo.Component, str]]:
827 """Report common header formulation issues for flowsheet diagnostics."""
828 problems = []
830 for time in self.flowsheet().time:
831 balance_flow = value(self.balance_flow_mol[time], exception=False)
832 condensate_flow = value(
833 self.outlet_condensate_state[time].flow_mol,
834 exception=False,
835 )
836 total_flow = value(self.total_flow_mol[time], exception=False)
838 if balance_flow is not None and abs(balance_flow) > 1e-6:
839 problems.append(
840 (
841 self.balance_flow_mol[time],
842 f"Balance flow is {balance_flow:.6g} mol/s. "
843 "This means the header is adding flow from nowhere "
844 "to make the material balance close. Replace balance flow "
845 "with an inlet flow variable to specify the makeup.",
846 )
847 )
849 if (
850 not self.config.is_liquid_header
851 and condensate_flow is not None
852 and total_flow is not None
853 ):
854 condensate_tolerance = max(abs(total_flow) * 1e-3, 1e-6)
855 if (
856 total_flow > 1e-6
857 and abs(condensate_flow - total_flow) <= condensate_tolerance
858 ):
859 problems.append(
860 (
861 self.total_flow_mol[time],
862 f"Condensate flow ({condensate_flow:.6g} mol/s) "
863 f"is approximately equal to total header flow "
864 f"({total_flow:.6g} mol/s). This means everything "
865 "entering the steam header is liquid after mixing, "
866 "so the steam header will fail to solve. Check the "
867 "inlet conditions and ensure there is sufficient "
868 "vapour entering the header.",
869 )
870 )
872 return problems
875 # -----------------------------------------------------------------
876 # Common utilities
877 # -----------------------------------------------------------------
879 def calculate_scaling_factors(self):
880 super().calculate_scaling_factors()
881 iscale.set_scaling_factor(self.heat_loss, 1e-6)
882 iscale.set_scaling_factor(self.pressure_loss, 1e-6)
883 iscale.set_scaling_factor(self.balance_flow_mol, 1e-3)
884 iscale.set_scaling_factor(self._partial_total_flow_mol, 1e-3)
885 if hasattr(self, "_vap_out_enth_mol"):
886 iscale.set_scaling_factor(self._vap_out_enth_mol, 1e-6)
889 def _build_state_blocks(
890 self,
891 stream_name_list: Iterable[str],
892 has_phase_equilibrium: bool,
893 is_defined_state: Optional[bool] = False,
894 is_build_port: Optional[bool] = False,
895 ) -> List[StateBlock]:
896 blocks: List[StateBlock] = []
898 base_args = dict(self.config.property_package_args)
899 base_args["has_phase_equilibrium"] = has_phase_equilibrium
900 base_args["defined_state"] = is_defined_state
902 for stream_name in stream_name_list:
903 args = dict(base_args)
904 args["doc"] = f"Thermophysical properties at {stream_name}"
905 sb = self.config.property_package.build_state_block(self.flowsheet().time, **args)
906 setattr(self, f"{stream_name}_state", sb)
907 blocks.append(sb)
909 if is_build_port: # No port is needed for intermediate/internal state blocks
910 self.add_port(name=stream_name, block=sb)
912 return blocks
915 def _get_stream_table_contents(self, time_point=0):
916 io_dict = {name: getattr(self, name) for name in [*self.inlet_blocks, *self.outlet_blocks, *self.internal_blocks]}
917 return create_stream_table_dataframe(io_dict, time_point=time_point)
920 @staticmethod
921 def ahuora_metadata():
922 from ahuora_unit_ops.json_config import (
923 JsonAdapterArgConfig,
924 JsonFrontendConfig,
925 JsonGraphicObjectConfig,
926 JsonIdaesAdapterConfig,
927 JsonPortAdapterConfig,
928 JsonPortConfig,
929 JsonPropertyConfig,
930 JsonPropertySetGroupConfig,
931 JsonSchemaPortMappingConfig,
932 JsonUnitOpConfig,
933 )
935 return JsonUnitOpConfig(
936 key='simple_header',
937 objectType='simple_header',
938 enumMember='SimpleHeader',
939 displayType='Simple Header',
940 displayName='Simple Header',
941 categoryPath=['chemical', 'mixer'],
942 ports={
943 'outlet_condensate': JsonPortConfig(
944 displayName='Condensate',
945 type='outlet',
946 streamType='stream',
947 many=False,
948 default=1,
949 minimum=1,
950 makeStream=True,
951 streamOffset=3,
952 streamName='Condensate S',
953 ),
954 'outlet_vent': JsonPortConfig(
955 displayName='Vent Outlet',
956 type='outlet',
957 streamType='stream',
958 many=False,
959 default=1,
960 minimum=1,
961 makeStream=True,
962 streamOffset=3,
963 streamName='Vent S',
964 ),
965 'inlet': JsonPortConfig(
966 displayName='Inlet',
967 type='inlet',
968 streamType='stream',
969 many=True,
970 default=1,
971 minimum=1,
972 makeStream=True,
973 streamOffset=3,
974 streamName='S',
975 ),
976 'outlet': JsonPortConfig(
977 displayName='Outlet',
978 type='outlet',
979 streamType='stream',
980 many=True,
981 default=2,
982 minimum=1,
983 makeStream=True,
984 streamOffset=3,
985 streamName='S',
986 ),
987 },
988 propertyPackagePorts={
989 '': ['inlet', 'outlet', 'outlet_vent', 'outlet_condensate'],
990 },
991 graphicObject=JsonGraphicObjectConfig(
992 kind='explicit',
993 width=25,
994 height=800,
995 autoHeight=True,
996 ),
997 indexSets=['splitter_fraction'],
998 properties={
999 'split_flow': JsonPropertyConfig(
1000 propertySetGroup='default',
1001 displayName='Flow',
1002 indexSets=['splitter_fraction'],
1003 sumToOne=False,
1004 value=None,
1005 unit=None,
1006 unitType='molarflow',
1007 description=None,
1008 type='numeric',
1009 many=False,
1010 default=1,
1011 options={},
1012 hasTimeIndex=True,
1013 ),
1014 'outlet_vent.flow_mol': JsonPropertyConfig(
1015 propertySetGroup='default',
1016 displayName='Vent Flow',
1017 indexSets=None,
1018 sumToOne=False,
1019 value=None,
1020 unit=None,
1021 unitType='molarflow',
1022 description=None,
1023 type='numeric',
1024 many=False,
1025 default=1,
1026 options={},
1027 hasTimeIndex=True,
1028 ),
1029 'heat_loss': JsonPropertyConfig(
1030 propertySetGroup='default',
1031 displayName='Heat Loss',
1032 indexSets=None,
1033 sumToOne=False,
1034 value=0,
1035 unit=None,
1036 unitType='heatflow',
1037 description=None,
1038 type='numeric',
1039 many=False,
1040 default=1,
1041 options={},
1042 hasTimeIndex=True,
1043 ),
1044 'pressure_loss': JsonPropertyConfig(
1045 propertySetGroup='default',
1046 displayName='Pressure Drop',
1047 indexSets=None,
1048 sumToOne=False,
1049 value=0,
1050 unit=None,
1051 unitType='pressure',
1052 description=None,
1053 type='numeric',
1054 many=False,
1055 default=1,
1056 options={},
1057 hasTimeIndex=True,
1058 ),
1059 'total_flow_mass': JsonPropertyConfig(
1060 propertySetGroup='default',
1061 displayName='Mass Flow',
1062 indexSets=None,
1063 sumToOne=False,
1064 value=None,
1065 unit=None,
1066 unitType='massflow',
1067 description=None,
1068 type='numeric',
1069 many=False,
1070 default=1,
1071 options={},
1072 hasTimeIndex=True,
1073 ),
1074 'total_flow_mol': JsonPropertyConfig(
1075 propertySetGroup='default',
1076 displayName='Molar Flow',
1077 indexSets=None,
1078 sumToOne=False,
1079 value=None,
1080 unit=None,
1081 unitType='molarflow',
1082 description=None,
1083 type='numeric',
1084 many=False,
1085 default=1,
1086 options={},
1087 hasTimeIndex=True,
1088 ),
1089 'balance_flow_mol': JsonPropertyConfig(
1090 propertySetGroup='default',
1091 displayName='Balance Flow',
1092 indexSets=None,
1093 sumToOne=False,
1094 value=None,
1095 unit=None,
1096 unitType='molarflow',
1097 description=None,
1098 type='numeric',
1099 many=False,
1100 default=1,
1101 options={},
1102 hasTimeIndex=True,
1103 ),
1104 'temperature': JsonPropertyConfig(
1105 propertySetGroup='default',
1106 displayName='Temperature',
1107 indexSets=None,
1108 sumToOne=False,
1109 value=None,
1110 unit=None,
1111 unitType='temperature',
1112 description=None,
1113 type='numeric',
1114 many=False,
1115 default=1,
1116 options={},
1117 hasTimeIndex=True,
1118 ),
1119 'degree_of_superheat': JsonPropertyConfig(
1120 propertySetGroup='default',
1121 displayName='Degree of Superheat',
1122 indexSets=None,
1123 sumToOne=False,
1124 value=None,
1125 unit=None,
1126 unitType='deltaTemperature',
1127 description=None,
1128 type='numeric',
1129 many=False,
1130 default=1,
1131 options={},
1132 hasTimeIndex=True,
1133 ),
1134 'pressure': JsonPropertyConfig(
1135 propertySetGroup='default',
1136 displayName='Pressure',
1137 indexSets=None,
1138 sumToOne=False,
1139 value=None,
1140 unit=None,
1141 unitType='pressure',
1142 description=None,
1143 type='numeric',
1144 many=False,
1145 default=1,
1146 options={},
1147 hasTimeIndex=True,
1148 ),
1149 'vapor_frac': JsonPropertyConfig(
1150 propertySetGroup='default',
1151 displayName='Vapor Fraction',
1152 indexSets=None,
1153 sumToOne=False,
1154 value=None,
1155 unit=None,
1156 unitType='ratio',
1157 description=None,
1158 type='numeric',
1159 many=False,
1160 default=1,
1161 options={},
1162 hasTimeIndex=True,
1163 ),
1164 'enth_mass': JsonPropertyConfig(
1165 propertySetGroup='default',
1166 displayName='Mass Specific Enthalpy',
1167 indexSets=None,
1168 sumToOne=False,
1169 value=None,
1170 unit=None,
1171 unitType='massEnthalpy',
1172 description=None,
1173 type='numeric',
1174 many=False,
1175 default=1,
1176 options={},
1177 hasTimeIndex=True,
1178 ),
1179 'enth_mol': JsonPropertyConfig(
1180 propertySetGroup='default',
1181 displayName='Molar Specific Enthalpy',
1182 indexSets=None,
1183 sumToOne=False,
1184 value=None,
1185 unit=None,
1186 unitType='molarEnthalpy',
1187 description=None,
1188 type='numeric',
1189 many=False,
1190 default=1,
1191 options={},
1192 hasTimeIndex=True,
1193 ),
1194 },
1195 propertySetGroups={
1196 'default': JsonPropertySetGroupConfig(
1197 type='stateVars',
1198 displayName='Properties',
1199 stateVars=['split_flow', 'heat_loss', 'pressure_loss'],
1200 toggle=None,
1201 ),
1202 },
1203 keyProperties=['pressure', 'temperature', 'degree_of_superheat', 'enth_mass', 'vapor_frac'],
1204 splitterFractionName='Outlet',
1205 idaesAdapter=JsonIdaesAdapterConfig(
1206 constructor='ahuora_builder.custom.thermal_utility_systems.header.simple_header',
1207 args={
1208 'property_package': JsonAdapterArgConfig(
1209 kind='property_package',
1210 label=None,
1211 ),
1212 'num_inlets': JsonAdapterArgConfig(
1213 kind='port_count',
1214 port='inlet',
1215 ),
1216 'num_outlets': JsonAdapterArgConfig(
1217 kind='port_count',
1218 port='outlet',
1219 ),
1220 },
1221 ports=JsonPortAdapterConfig(
1222 kind='schema_ports',
1223 mappings=[
1224 JsonSchemaPortMappingConfig(
1225 type='group',
1226 port='inlet',
1227 outputKeyTemplate='inlet_{index}',
1228 inlet=True,
1229 connectedOnly=True,
1230 sortByIndex=True,
1231 ),
1232 JsonSchemaPortMappingConfig(
1233 type='group',
1234 port='outlet',
1235 outputKeyTemplate='outlet_{index}',
1236 inlet=False,
1237 connectedOnly=True,
1238 sortByIndex=True,
1239 ),
1240 JsonSchemaPortMappingConfig(
1241 type='fixed',
1242 port='outlet_vent',
1243 outputKey='outlet_vent',
1244 inlet=False,
1245 optional=True,
1246 connectedOnly=True,
1247 ),
1248 JsonSchemaPortMappingConfig(
1249 type='fixed',
1250 port='outlet_condensate',
1251 outputKey='outlet_condensate',
1252 inlet=False,
1253 optional=True,
1254 connectedOnly=True,
1255 ),
1256 ],
1257 ),
1258 properties=None,
1259 ),
1260 frontend=JsonFrontendConfig(
1261 showInPanel=True,
1262 variant=None,
1263 ),
1264 )