Coverage for backend/ahuora-compounds/ahuora_property_packages/humid_air/HumidAirSurrogate.py: 19%
191 statements
« prev ^ index » next coverage.py v7.10.7, created at 2026-07-22 05:22 +0000
1# Import Python libraries
2import logging
3import os
5# Import Pyomo libraries
6from pyomo.environ import (
7 Constraint,
8 Param,
9 Reals,
10 value,
11 Var,
12 NonNegativeReals,
13 units,
14)
16# Import IDAES cores
17from idaes.core import (
18 declare_process_block_class,
19 PhysicalParameterBlock,
20 StateBlockData,
21 StateBlock,
22 MaterialBalanceType,
23 EnergyBalanceType,
24 VaporPhase,
25 Component,
26)
27from idaes.core.util.model_statistics import degrees_of_freedom
29from idaes.core.surrogate.surrogate_block import SurrogateBlock
30from idaes.core.surrogate.pysmo_surrogate import PysmoSurrogate
31import idaes.logger as idaeslog
33from pyomo.environ import Block
34from pyomo.core.base.expression import ScalarExpression, Expression, _GeneralExpressionData, ExpressionData
35from pyomo.core.base.var import ScalarVar, _GeneralVarData, VarData, IndexedVar
38# Some more information about this module
39__author__ = "Stephen Burroughs"
42# Set up logger
43_log = logging.getLogger(__name__)
45class _StateBlock(StateBlock):
46 """
47 This Class contains methods which should be applied to Property Blocks as a
48 whole, rather than individual elements of indexed Property Blocks.
49 """
51 def initialize(
52 blk,
53 state_args=None,
54 hold_state=False,
55 outlvl=1,
56 state_vars_fixed=False,
57 solver="ipopt",
58 optarg={"tol": 1e-8},
59 ):
60 """
61 Initialisation routine for property package.
63 Keyword Arguments:
64 flow_mol : value at which to initialize component flows
65 (default=None)
66 pressure : value at which to initialize pressure (default=None)
67 temperature : value at which to initialize temperature
68 mole_flow_frac: value at which to initialise the molar flow fraction
69 (default=None)
70 outlvl : sets output level of initialisation routine
72 * 0 = no output (default)
73 * 1 = return solver state for each step in routine
74 * 2 = include solver output information (tee=True)
75 state_vars_fixed: Flag to denote if state vars have already been
76 fixed.
77 - True - states have already been fixed by the
78 control volume 1D. Control volume 0D
79 does not fix the state vars, so will
80 be False if this state block is used
81 with 0D blocks.
82 - False - states have not been fixed. The state
83 block will deal with fixing/unfixing.
84 optarg : solver options dictionary object (default=None)
85 solver : str indicating which solver to use during
86 initialization (default = 'ipopt')
87 hold_state : flag indicating whether the initialization routine
88 should unfix any state variables fixed during
89 initialization (default=False).
90 - True - states variables are not unfixed, and
91 a dict of returned containing flags for
92 which states were fixed during
93 initialization.
94 - False - state variables are unfixed after
95 initialization by calling the
96 release_state method
98 Returns:
99 If hold_states is True, returns a dict containing flags for
100 which states were fixed during initialization.
101 """
103 if state_vars_fixed is False:
104 # Fix state variables if not already fixed
105 Fcflag = {}
106 Pflag = {}
107 Tflag = {}
108 # Fmflag = {}
110 for k in blk.keys():
111 # if blk[k].mole_frac_comp["water"].fixed is True:
112 # Fmflag[k] = True
113 # else:
114 # Fmflag[k] = False
115 # if state_args is None:
116 # blk[k].mole_frac_comp["water"].fix()
117 # else:
118 # blk[k].mole_frac_comp["water"].fix(state_args["flow_mol_water"])
120 if blk[k].flow_mol.fixed is True:
121 Fcflag[k] = True
122 else:
123 Fcflag[k] = False
124 if state_args is None:
125 blk[k].flow_mol.fix()
126 else:
127 blk[k].flow_mol.fix(state_args["flow_mol"])
129 if blk[k].pressure.fixed is True:
130 Pflag[k] = True
131 else:
132 Pflag[k] = False
133 if state_args is None:
134 blk[k].pressure.fix()
135 else:
136 blk[k].pressure.fix(state_args["pressure"])
138 if blk[k].temperature.fixed is True:
139 Tflag[k] = True
140 else:
141 Tflag[k] = False
142 if state_args is None:
143 blk[k].temperature.fix()
144 else:
145 blk[k].temperature.fix(state_args["temperature"])
147 # If input block, return flags, else release state
148 flags = {"Fcflag": Fcflag, "Pflag": Pflag, "Tflag": Tflag}
150 else:
151 # Check when the state vars are fixed already result in dof 0
152 for k in blk.keys():
153 if degrees_of_freedom(blk[k]) != 0:
154 raise Exception(
155 "State vars fixed but degrees of freedom "
156 "for state block is not zero during "
157 "initialization."
158 )
160 if state_vars_fixed is False:
161 if hold_state is True:
162 return flags
163 else:
164 blk.release_state(flags)
166 def release_state(blk, flags, outlvl=0):
167 """
168 Method to release state variables fixed during initialisation.
170 Keyword Arguments:
171 flags : dict containing information of which state variables
172 were fixed during initialization, and should now be
173 unfixed. This dict is returned by initialize if
174 hold_state=True.
175 outlvl : sets output level of of logging
176 """
177 if flags is None:
178 return
180 # Unfix state variables
181 for k in blk.keys():
182 if flags["Fcflag"][k] is False:
183 blk[k].flow_mol.unfix()
184 if flags["Pflag"][k] is False:
185 blk[k].pressure.unfix()
186 if flags["Tflag"][k] is False:
187 blk[k].temperature.unfix()
188 # if flags["Fmflag"][k] is False:
189 # blk[k].mole_frac_comp["water"].unfix()
191 if outlvl > 0:
192 if outlvl > 0:
193 _log.info("{} State Released.".format(blk.name))
195class _StateBlockWrapper(_StateBlock):
197 def initialize(blk, *args, **kwargs):
198 for v, k in blk.items():
199 print(k)
200 k.constraints.deactivate()
201 return _StateBlock.initialize(blk, *args, **kwargs)
203 def release_state(blk, flags, outlvl=idaeslog.NOTSET):
204 _StateBlock.release_state(blk, flags, outlvl)
206 for v, k in blk.items():
207 k.constraints.activate()
209@declare_process_block_class("HAirStateBlock", block_class=_StateBlockWrapper)
210class HAirStateBlockData(StateBlockData):
211 """
212 An example property package for ideal gas properties with Gibbs energy
213 """
215 def build(self):
216 """
217 Callable method for Block construction
218 """
219 super(HAirStateBlockData, self).build()
220 self.constraints = Block()
221 self._make_state_vars()
222 if self.config.defined_state is False:
223 self.sum_mole_frac_out = Constraint(
224 expr = 1.0 == sum(self.mole_frac_comp[i] for i in self.component_list)
225 )
227 def constrain_component(blk, component: Var | Expression, value: float) -> Constraint | Var | None:
228 """
229 Constrain a component to a value
230 """
231 if isinstance(component, ScalarExpression):
232 c = Constraint(expr=component == value)
233 c.defining_state_var = True
234 blk.constraints.add_component(component.local_name, c)
235 return c
236 elif type(component) in (ScalarVar, _GeneralVarData, VarData, IndexedVar):
237 component.fix(value)
238 return component
239 elif type(component) in (_GeneralExpressionData, ExpressionData):
240 # allowed, but we don't need to fix it (eg. mole_frac_comp in helmholtz)
241 return None
242 else:
243 raise Exception(
244 f"Component {component} is not a Var or Expression: {type(component)}"
245 )
247 def _make_state_vars(self):
249 self.flow_mol = Var(
250 domain=NonNegativeReals,
251 initialize=1.0,
252 units=units.mol / units.s,
253 doc="Total molar flowrate [mol/s]",
254 )
255 self.pressure = Var(
256 domain=NonNegativeReals,
257 initialize=95000,
258 bounds=(10000, 900000),
259 units=units.Pa,
260 doc="State pressure [Pa]",
261 )
263 self.temperature = Var(
264 domain=NonNegativeReals,
265 initialize=350,
266 bounds=(193.15, 273.15+400),
267 units=units.K,
268 doc="Dry bulb temperature [K]",
269 )
271 self.temperature_wet_bulb = Var(
272 domain=NonNegativeReals,
273 bounds=(130, 350 + 273.15),
274 units=units.K,
275 doc="Wet bulb temperature [K]",
276 )
278 self.mole_frac_comp = Var(
279 self.params.component_list,
280 initialize = 1/len(self.params.component_list),
281 domain=Reals,
282 bounds = (0,1),
283 units=units.dimensionless
284 )
286 self.relative_humidity = Var(
287 domain=Reals,
288 initialize=0.3,
289 bounds = (0,1)
290 )
292 self.enth_mol = Var(
293 domain=Reals,
294 initialize=300,
295 units = units.J / units.mol,
296 doc = "Enthalpy [J/mol]"
297 )
299 self.entr_mol = Var(
300 domain=Reals,
301 initialize=40,
302 units = units.J / units.mol / units.K,
303 doc = "Entropy [J/mol/K]"
304 )
306 self.vol_mol = Var(
307 domain = NonNegativeReals,
308 initialize=40,
309 units = units.m**3 / units.mol
310 )
312 inputs = [self.temperature, self.pressure, self.mole_frac_comp["water"], self.mole_frac_comp["air"]]
313 outputs = [self.relative_humidity, self.temperature_wet_bulb, self.enth_mol, self.entr_mol, self.vol_mol]
314 script_dir = os.path.dirname(__file__)
315 self.pysmo_surrogate = PysmoSurrogate.load_from_file(
316 os.path.join(script_dir,"pysmo_humid_air.json")
317 )
318 self.surrogate = SurrogateBlock()
319 self.surrogate.build_model(
320 self.pysmo_surrogate,
321 input_vars=inputs,
322 output_vars=outputs,
323 )
325 def _vol_mass(self):
326 def _vol_mass_rule(b):
327 return b.vol_mol / sum(
328 b.mole_frac_comp[i]
329 * (1/b.params.mw_comp[i])
330 for i in b.params.component_list
331 )
332 self.vol_mass = Expression(rule=_vol_mass_rule)
334 def _enth_mass(self):
335 def enth_mass_rule(b):
336 return sum(b.enth_mass_comp[i] for i in b.params.component_list)
337 self.enth_mass = Expression (rule=enth_mass_rule)
339 def _enth_mass_comp(self):
340 def _rule_enth_mass_comp(b, i):
341 return b.enth_mol_comp[i] / b.params.mw_comp[i]
342 self.enth_mass_comp = Expression(
343 self.params.component_list,
344 rule=_rule_enth_mass_comp,
345 )
347 def _enth_mol_comp(self):###check this
348 def _rule_enth_mol_comp(b, i):
349 return b.enth_mol * b.mole_frac_comp[i]
350 self.enth_mol_comp = Expression(
351 self.params.component_list,
352 rule=_rule_enth_mol_comp,
353 )
355 def _entr_mass(self):
356 def entr_mass_rule(b):
357 return sum (b.entr_mass_comp[i] for i in b.params.component_list)
358 self.entr_mass = Expression (rule=entr_mass_rule)
360 def _entr_mass_comp(self):
361 def _rule_entr_mass_comp(b, i):
362 return b.entr_mol_comp[i] /b.params.mw_comp[i]
363 self.entr_mass_comp = Expression(
364 self.params.component_list,
365 rule=_rule_entr_mass_comp,
366 )
368 def _entr_mol_comp(self):
369 def _rule_entr_mol_comp(b, i):
370 return b.entr_mol * b.mole_frac_comp[i]
371 self.entr_mol_comp = Expression(
372 self.params.component_list,
373 rule=_rule_entr_mol_comp,
374 )
376 def _flow_mass(self):
377 def flow_mass_rule(b):
378 return sum(b.flow_mass_comp[i] for i in b.params.component_list)
379 self.flow_mass = Expression(rule = flow_mass_rule)
381 def _flow_mass_comp(self):
382 def _rule_flow_mass_comp(b, i):
383 return b.flow_mol_comp[i] * b.params.mw_comp[i]
384 self.flow_mass_comp = Expression(
385 self.params.component_list,
386 rule=_rule_flow_mass_comp,
387 )
389 def _flow_mol_comp(self):
390 def _rule_flow_mol_comp(b, i):
391 return b.mole_frac_comp[i] * b.flow_mol
392 self.flow_mol_comp = Expression(
393 self.params.component_list,
394 rule=_rule_flow_mol_comp,
395 )
397 def _flow_vol(self):
398 def _rule_flow_vol(b):
399 return b.vol_mol*b.flow_mol
400 self.flow_vol = Expression(rule=_rule_flow_vol)
402 def _mass_frac_comp(self):
403 def _mass_frac_comp_rule(b, i):
404 return b.flow_mol_comp[i] * b.params.mw_comp[i]
405 self.mass_frac_comp = Expression ( self.params.component_list, rule = _mass_frac_comp_rule)
407 def _total_energy_flow(self):
408 def _rule_total_energy_flow(b):
409 return b.flow_mass * b.enth_mass
410 self.total_energy_flow = Expression( rule = _rule_total_energy_flow)
412 # def _vapor_frac(self):
413 # self.vapor_frac = Var(
414 # domain = NonNegativeReals,
415 # initialize=1.0,
416 # )
418 def get_material_flow_terms(self, p, c):
419 return self.mole_frac_comp[c]*self.flow_mol
421 def get_enthalpy_flow_terms(self , p):
422 return self.flow_mol * self.enth_mol
424 def default_material_balance_type(self):
425 return MaterialBalanceType.componentTotal
427 def default_energy_balance_type(self):
428 return EnergyBalanceType.enthalpyTotal
431 def define_state_vars(self):
432 return {
433 "flow_mol": self.flow_mol,
434 "temperature": self.temperature,
435 "pressure": self.pressure,
436 "mole_frac_comp": self.mole_frac_comp
437 }
440 def model_check(blk):
441 """
442 Model checks for property block
443 """
444 # Check temperature bounds
445 if value(blk.temperature) < blk.temperature.lb:
446 _log.error("{} temperature set below lower bound.".format(blk.name))
447 if value(blk.temperature) > blk.temperature.ub:
448 _log.error("{} temperature set above upper bound.".format(blk.name))
450 # Check pressure bounds
451 if value(blk.pressure) < blk.pressure.lb:
452 _log.error("{} Pressure set below lower bound.".format(blk.name))
453 if value(blk.pressure) > blk.pressure.ub:
454 _log.error("{} Pressure set above upper bound.".format(blk.name))
457@declare_process_block_class("HAirParameterBlock")
458class PhysicalParameterData(PhysicalParameterBlock):
459 """
460 Property Parameter Block Class
462 Contains parameters and indexing sets associated with properties for
463 supercritical Humid Air
465 """
467 def build(self):
468 """
469 Callable method for Block construction.
470 """
471 super(PhysicalParameterData, self).build()
473 self._state_block_class = HAirStateBlock # noqa: F821
474 # List of valid phases in property package
475 self.Vap = VaporPhase()
477 # Component list - a list of component identifiers
478 self.water = Component()
479 self.air = Component()
480 mw_comp_dict = {"water":0.01801528, "air": 0.02896}
481 self.mw_comp = Param(
482 self.component_list,
483 mutable=False,
484 initialize=mw_comp_dict,
485 doc="Molecular weights of components [kg/mol]",
486 units=units.kg / units.mol,
487 )
489 @classmethod
490 def define_metadata(cls, obj):
491 obj.add_properties(
492 {
493 "flow_mol": {"method": None, "units": units.mol / units.s},
494 "flow_mol_comp": {"method": "_flow_mol_comp"},
495 "flow_mass": {"method": "_flow_mass", "units": units.kg / units.s},
496 "flow_mass_comp": {"method": "_flow_mass_comp"},
497 "flow_vol": {"method":"_flow_vol", "units": units.m**3 / units.s},
498 "pressure": {"method": None, "units": units.Pa},
499 "vol_mol" : {"method": None, "units": units.m**3 / units.mol},
500 "vol_mass": {"method": "_vol_mass", "units": units.m**3 / units.kg},
501 "enth_mol": {"method": None, "units": units.J / units.mol},
502 "enth_mol_comp": {"method": "_enth_mol_comp"},
503 "enth_mass": {"method": "_enth_mass", "units": units.J/units.kg},
504 "enth_mass_comp": {"method": "_enth_mass_comp"},
505 "mole_frac_comp": {"method": "_mole_frac_comp"},
506 "mass_frac_comp": {"method": "_mass_frac_comp"},
507 "entr_mol": {"method": None, "units": units.J / units.mol / units.K},
508 "entr_mol_comp": {"method": "_entr_mol_comp"},
509 "entr_mass": {"method": "_entr_mass", "units": units.J / units.kg / units.K},
510 "entr_mass_comp": {"method": "_entr_mass_comp"},
511 "temperature": {"method": None, "units": units.K},
512 "total_energy_flow": {"method": "_total_energy_flow", "units": units.kW},
513 # "vapor_frac": {"method": "_vapor_frac"}
515 }
516 )
518 obj.define_custom_properties(
519 {
520 "temperature_wet_bulb": {"method": None},
521 "relative_humidity": {"method": None},
522 }
523 )
525 obj.add_default_units(
526 {
527 "time": units.s,
528 "length": units.m,
529 "mass": units.kg,
530 "amount": units.mol,
531 "temperature": units.K,
532 }
533 )