Coverage for backend/idaes_service/solver/custom/thermal_utility_systems/simple_heat_pump.py: 60%
157 statements
« prev ^ index » next coverage.py v7.10.7, created at 2025-11-06 23:27 +0000
1#################################################################################
2# The Institute for the Design of Advanced Energy Systems Integrated Platform
3# Framework (IDAES IP) was produced under the DOE Institute for the
4# Design of Advanced Energy Systems (IDAES).
5#
6# Copyright (c) 2018-2024 by the software owners: The Regents of the
7# University of California, through Lawrence Berkeley National Laboratory,
8# National Technology & Engineering Solutions of Sandia, LLC, Carnegie Mellon
9# University, West Virginia University Research Corporation, et al.
10# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md
11# for full copyright and license information.
12#################################################################################
13"""
14Heat Exchanger Models.
15"""
17__author__ = "Team Ahuora"
21# Import Pyomo libraries
22from pyomo.environ import (
23 Block,
24 Var,
25 Param,
26 log,
27 Reference,
28 PositiveReals,
29 ExternalFunction,
30 units as pyunits,
31 check_optimal_termination,
32)
33from pyomo.common.config import ConfigBlock, ConfigValue, In
35# Import IDAES cores
36from idaes.core import (
37 declare_process_block_class,
38 UnitModelBlockData,
39)
41import idaes.logger as idaeslog
42from idaes.core.util.functions import functions_lib
43from idaes.core.util.tables import create_stream_table_dataframe
44from idaes.models.unit_models.heater import (
45 _make_heater_config_block,
46 _make_heater_control_volume,
47)
49from idaes.core.util.misc import add_object_reference
50from idaes.core.util import scaling as iscale
51from idaes.core.solvers import get_solver
52from idaes.core.util.exceptions import ConfigurationError, InitializationError
53from idaes.core.initialization import SingleControlVolumeUnitInitializer
55_log = idaeslog.getLogger(__name__)
58class SimpleHeatPumpInitializer(SingleControlVolumeUnitInitializer):
59 """
60 Initializer for 0D Heat Exchanger units.
62 """
64 def initialization_routine(
65 self,
66 model: Block,
67 plugin_initializer_args: dict = None,
68 copy_inlet_state: bool = False,
69 duty=1000 * pyunits.W,
70 ):
71 """
72 Common initialization routine for 0D Heat Exchangers.
74 This routine starts by initializing the hot and cold side properties. Next, the heat
75 transfer between the two sides is fixed to an initial guess for the heat duty (provided by the duty
76 argument), the associated constraint is deactivated, and the model is then solved. Finally, the heat
77 duty is unfixed and the heat transfer constraint reactivated followed by a final solve of the model.
79 Args:
80 model: Pyomo Block to be initialized
81 plugin_initializer_args: dict-of-dicts containing arguments to be passed to plug-in Initializers.
82 Keys should be submodel components.
83 copy_inlet_state: bool (default=False). Whether to copy inlet state to other states or not
84 (0-D control volumes only). Copying will generally be faster, but inlet states may not contain
85 all properties required elsewhere.
86 duty: initial guess for heat duty to assist with initialization. Can be a Pyomo expression with units.
88 Returns:
89 Pyomo solver results object
90 """
91 return super(SingleControlVolumeUnitInitializer, self).initialization_routine(
92 model=model,
93 plugin_initializer_args=plugin_initializer_args,
94 copy_inlet_state=copy_inlet_state,
95 duty=duty,
96 )
98 def initialize_main_model(
99 self,
100 model: Block,
101 copy_inlet_state: bool = False,
102 duty=1000 * pyunits.W,
103 ):
104 """
105 Initialization routine for main 0D HX models.
107 Args:
108 model: Pyomo Block to be initialized.
109 copy_inlet_state: bool (default=False). Whether to copy inlet state to other states or not
110 (0-D control volumes only). Copying will generally be faster, but inlet states may not contain
111 all properties required elsewhere.
112 duty: initial guess for heat duty to assist with initialization, default = 1000 W. Can
113 be a Pyomo expression with units.
115 Returns:
116 Pyomo solver results object.
118 """
119 # Get loggers
120 init_log = idaeslog.getInitLogger(
121 model.name, self.get_output_level(), tag="unit"
122 )
123 solve_log = idaeslog.getSolveLogger(
124 model.name, self.get_output_level(), tag="unit"
125 )
127 # Create solver
128 solver = self._get_solver()
130 self.initialize_control_volume(model.source, copy_inlet_state)
131 init_log.info_high("Initialization Step 1a (heat output (sink)) Complete.")
133 self.initialize_control_volume(model.sink, copy_inlet_state)
134 init_log.info_high("Initialization Step 1b (cold side) Complete.")
135 # ---------------------------------------------------------------------
136 # Solve unit without heat transfer equation
137 model.heat_transfer_equation.deactivate()
139 # Check to see if heat duty is fixed
140 # We will assume that if the first point is fixed, it is fixed at all points
141 if not model.sink.heat[model.flowsheet().time.first()].fixed:
142 cs_fixed = False
144 model.sink.heat.fix(duty)
145 for i in model.source.heat:
146 model.source.heat[i].set_value(-duty)
147 else:
148 cs_fixed = True
149 for i in model.source.heat:
150 model.source.heat[i].set_value(model.sink.heat[i])
152 with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
153 res = solver.solve(model, tee=slc.tee)
154 init_log.info_high("Initialization Step 2 {}.".format(idaeslog.condition(res)))
155 if not cs_fixed:
156 model.sink.heat.unfix()
157 model.heat_transfer_equation.activate()
158 # ---------------------------------------------------------------------
159 # Solve unit
160 with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
161 res = solver.solve(model, tee=slc.tee)
162 init_log.info("Initialization Completed, {}".format(idaeslog.condition(res)))
164 return res
167def _make_heat_pump_config(config):
168 """
169 Declare configuration options for HeatExchangerData block.
170 """
171 config.declare(
172 "source",
173 ConfigBlock(
174 description="Config block for heat output (sink)",
175 doc="""A config block used to construct the heat output (sink) control volume.
176This config can be given by the heat output (sink) name instead of source.""",
177 ),
178 )
179 config.declare(
180 "sink",
181 ConfigBlock(
182 description="Config block for cold side",
183 doc="""A config block used to construct the cold side control volume.
184This config can be given by the cold side name instead of sink.""",
185 ),
186 )
187 _make_heater_config_block(config.source)
188 _make_heater_config_block(config.sink)
191@declare_process_block_class("SimpleHeatPump", doc="Simple 0D heat pump model.")
192class SimpleHeatPumpData(UnitModelBlockData):
193 """
194 Simple 0D heat pump unit.
195 Unit model to transfer heat from one material to another.
196 """
198 default_initializer = SimpleHeatPumpInitializer
200 CONFIG = UnitModelBlockData.CONFIG(implicit=True)
201 _make_heat_pump_config(CONFIG)
203 def build(self):
204 """
205 Building model
207 Args:
208 None
209 Returns:
210 None
211 """
212 ########################################################################
213 # Call UnitModel.build to setup dynamics and configure #
214 ########################################################################
215 super().build()
216 config = self.config
218 ########################################################################
219 # Add control volumes #
220 ########################################################################
221 source = _make_heater_control_volume(
222 self,
223 "source",
224 config.source,
225 dynamic=config.dynamic,
226 has_holdup=config.has_holdup,
227 )
228 sink = _make_heater_control_volume(
229 self,
230 "sink",
231 config.sink,
232 dynamic=config.dynamic,
233 has_holdup=config.has_holdup,
234 )
236 ########################################################################
237 # Add variables #
238 ########################################################################
239 # Use heat output (sink) units as basis
240 s1_metadata = self.source.config.property_package.get_metadata()
242 q_units = s1_metadata.get_derived_units("power")
243 temp_units = s1_metadata.get_derived_units("temperature")
245 self.work_mechanical = Var(
246 self.flowsheet().time,
247 domain=PositiveReals,
248 initialize=100.0,
249 doc="Mechanical work input",
250 units=q_units,
251 )
252 self.coefficient_of_performance = Var(
253 domain=PositiveReals,
254 initialize=2.0,
255 doc="Coefficient of performance",
256 units=pyunits.dimensionless,
257 )
258 self.efficiency = Var(
259 domain=PositiveReals,
260 initialize=0.5,
261 doc="Efficiency (Carnot = 1)",
262 units=pyunits.dimensionless,
263 )
264 self.delta_temperature_lift = Var(
265 self.flowsheet().time,
266 domain=PositiveReals,
267 initialize=10.0,
268 doc="Temperature lift between source inlet and sink outlet",
269 units=temp_units,
270 )
272 self.approach_temperature = Var(
273 domain=PositiveReals,
274 initialize=10.0,
275 doc="Approach temperature of refrigerant",
276 units=temp_units,
277 )
279 self.heat_duty = Reference(sink.heat)
280 ########################################################################
281 # Add ports #
282 ########################################################################
283 self.add_inlet_port(name="source_inlet", block=source, doc="heat output (sink) inlet")
284 self.add_inlet_port(
285 name="sink_inlet",
286 block=sink,
287 doc="Cold side inlet",
288 )
289 self.add_outlet_port(
290 name="source_outlet", block=source, doc="heat output (sink) outlet"
291 )
292 self.add_outlet_port(
293 name="sink_outlet",
294 block=sink,
295 doc="Cold side outlet",
296 )
298 ########################################################################
299 # Add temperature lift constraints #
300 ########################################################################
302 @self.Constraint(self.flowsheet().time)
303 def delta_temperature_lift_equation(b, t):
304 # Refrigerant saturation levels (K)
305 T_cond = b.sink.properties_out[0].temperature + b.approach_temperature
306 T_evap = b.source.properties_out[0].temperature - b.approach_temperature
307 # Positive temperature lift
308 return b.delta_temperature_lift[t] == T_cond - T_evap
310 ########################################################################
311 # Add a unit level energy balance #
312 ########################################################################
313 @self.Constraint(self.flowsheet().time)
314 def unit_heat_balance(b, t):
315 # The duty of the sink = the source + work input
316 return 0 == (
317 source.heat[t] + pyunits.convert(sink.heat[t], to_units=q_units) - b.work_mechanical[t]
318 )
321 ########################################################################
322 # Add Heat transfer equation #
323 ########################################################################
325 @self.Constraint(self.flowsheet().time)
326 def heat_transfer_equation(b, t):
327 # This equation defines the heat duty using the cop of heating.
328 return sink.heat[t] == b.coefficient_of_performance * b.work_mechanical[t]
330 ########################################################################
331 # Add COP Equation #
332 ########################################################################
333 @self.Constraint()
334 def cop_equation(b):
335 # Refrigerant saturation levels (K)
336 T_cond = b.sink.properties_out[0].temperature + b.approach_temperature
337 T_evap = b.source.properties_out[0].temperature - b.approach_temperature
338 # Carnot CoP with efficiency
339 return b.coefficient_of_performance == (T_cond / (T_cond - T_evap)) * b.efficiency
342 ########################################################################
343 # Add symbols for LaTeX equation rendering #
344 ########################################################################
345 self.work_mechanical.latex_symbol = "W_{mech}"
346 self.coefficient_of_performance.latex_symbol = "COP"
347 self.efficiency.latex_symbol = "\\eta"
348 self.heat_duty.latex_symbol = "Q_{HP}"
349 self.delta_temperature_lift.latex_symbol = "\\Delta T_{lift}"
352 def initialize_build(
353 self,
354 state_args_1=None,
355 state_args_2=None,
356 outlvl=idaeslog.NOTSET,
357 solver=None,
358 optarg=None,
359 duty=None,
360 ):
361 """
362 Heat exchanger initialization method.
364 Args:
365 state_args_1 : a dict of arguments to be passed to the property
366 initialization for the heat output (sink) (see documentation of the specific
367 property package) (default = {}).
368 state_args_2 : a dict of arguments to be passed to the property
369 initialization for the cold side (see documentation of the specific
370 property package) (default = {}).
371 outlvl : sets output level of initialization routine
372 optarg : solver options dictionary object (default=None, use
373 default solver options)
374 solver : str indicating which solver to use during
375 initialization (default = None, use default solver)
376 duty : an initial guess for the amount of heat transferred. This
377 should be a tuple in the form (value, units), (default
378 = (1000 J/s))
380 Returns:
381 None
383 """
384 # Set solver options
385 init_log = idaeslog.getInitLogger(self.name, outlvl, tag="unit")
386 solve_log = idaeslog.getSolveLogger(self.name, outlvl, tag="unit")
388 # Create solver
389 opt = get_solver(solver, optarg)
391 flags1 = self.source.initialize(
392 outlvl=outlvl, optarg=optarg, solver=solver, state_args=state_args_1
393 )
395 init_log.info_high("Initialization Step 1a (heat output (sink)) Complete.")
397 flags2 = self.sink.initialize(
398 outlvl=outlvl, optarg=optarg, solver=solver, state_args=state_args_2
399 )
401 init_log.info_high("Initialization Step 1b (cold side) Complete.")
402 # ---------------------------------------------------------------------
403 # Solve unit without heat transfer equation
404 self.heat_transfer_equation.deactivate()
406 # Get side 1 and side 2 heat units, and convert duty as needed
407 s1_units = self.source.heat.get_units()
408 s2_units = self.sink.heat.get_units()
410 # Check to see if heat duty is fixed
411 # WE will assume that if the first point is fixed, it is fixed at all points
412 if not self.sink.heat[self.flowsheet().time.first()].fixed: 412 ↛ 440line 412 didn't jump to line 440 because the condition on line 412 was always true
413 cs_fixed = False
414 if duty is None: 414 ↛ 429line 414 didn't jump to line 429 because the condition on line 414 was always true
415 # Assume 1000 J/s and check for unitless properties
416 if s1_units is None and s2_units is None: 416 ↛ 418line 416 didn't jump to line 418 because the condition on line 416 was never true
417 # Backwards compatibility for unitless properties
418 s1_duty = -1000
419 s2_duty = 1000
420 else:
421 s1_duty = pyunits.convert_value(
422 -1000, from_units=pyunits.W, to_units=s1_units
423 )
424 s2_duty = pyunits.convert_value(
425 1000, from_units=pyunits.W, to_units=s2_units
426 )
427 else:
428 # Duty provided with explicit units
429 s1_duty = -pyunits.convert_value(
430 duty[0], from_units=duty[1], to_units=s1_units
431 )
432 s2_duty = pyunits.convert_value(
433 duty[0], from_units=duty[1], to_units=s2_units
434 )
436 self.sink.heat.fix(s2_duty)
437 for i in self.source.heat:
438 self.source.heat[i].value = s1_duty
439 else:
440 cs_fixed = True
441 for i in self.source.heat:
442 self.source.heat[i].set_value(self.sink.heat[i])
444 with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
445 res = opt.solve(self, tee=slc.tee)
446 init_log.info_high("Initialization Step 2 {}.".format(idaeslog.condition(res)))
447 if not cs_fixed: 447 ↛ 449line 447 didn't jump to line 449 because the condition on line 447 was always true
448 self.sink.heat.unfix()
449 self.heat_transfer_equation.activate()
450 # ---------------------------------------------------------------------
451 # Solve unit
452 with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
453 res = opt.solve(self, tee=slc.tee)
454 init_log.info_high("Initialization Step 3 {}.".format(idaeslog.condition(res)))
455 # ---------------------------------------------------------------------
456 # Release Inlet state
457 self.source.release_state(flags1, outlvl=outlvl)
458 self.sink.release_state(flags2, outlvl=outlvl)
460 init_log.info("Initialization Completed, {}".format(idaeslog.condition(res)))
462 if not check_optimal_termination(res): 462 ↛ 463line 462 didn't jump to line 463 because the condition on line 462 was never true
463 raise InitializationError(
464 f"{self.name} failed to initialize successfully. Please check "
465 f"the output logs for more information."
466 )
468 def _get_performance_contents(self, time_point=0):
469 # this shows as the performance table when calling the report method
470 var_dict = {
472 "Coefficient of Performance": self.coefficient_of_performance,
473 "Efficiency": self.efficiency,
474 "Source Heat Duty": self.source.heat[time_point],
475 "Mechanical Work Input": self.work_mechanical[time_point],
476 "Sink Heat Duty": self.sink.heat[time_point],
477 "Temperature Lift": self.delta_temperature_lift[time_point],
478 }
480 expr_dict = {}
482 return {"vars": var_dict, "exprs": expr_dict}
484 def _get_stream_table_contents(self, time_point=0):
485 # this shows as t he stream table when calling the report method
486 # Get names for hot and cold sides
487 return create_stream_table_dataframe(
488 {
489 f"Source Inlet": self.source_inlet,
490 f"Source Outlet": self.source_outlet,
491 f"Sink Inlet": self.sink_inlet,
492 f"Sink Outlet": self.sink_outlet,
493 },
494 time_point=time_point,
495 )
497 def calculate_scaling_factors(self):
498 super().calculate_scaling_factors()
499 # TODO: Review this code to check it makes sense.
501 # Scaling for heat pump variables
502 # Mechanical work input: typical values vary, set default scaling to 0.01
503 sf_work = dict(
504 zip(
505 self.work_mechanical.keys(),
506 [
507 iscale.get_scaling_factor(v, default=0.01)
508 for v in self.work_mechanical.values()
509 ],
510 )
511 )
512 # COP and efficiency are dimensionless, usually between 1 and 10, default scaling 0.1
513 sf_cop = iscale.get_scaling_factor(self.coefficient_of_performance, default=0.1)
514 sf_eff = iscale.get_scaling_factor(self.efficiency, default=0.1)
516 # Delta Ts: typical range 1-100, default scaling 0.1
517 sf_dT_in = dict(
518 zip(
519 self.delta_temperature_in.keys(),
520 [
521 iscale.get_scaling_factor(v, default=0.1)
522 for v in self.delta_temperature_in.values()
523 ],
524 )
525 )
526 sf_dT_out = dict(
527 zip(
528 self.delta_temperature_out.keys(),
529 [
530 iscale.get_scaling_factor(v, default=0.1)
531 for v in self.delta_temperature_out.values()
532 ],
533 )
534 )
536 # Heat duty: depends on process, default scaling 0.01
537 sf_q = dict(
538 zip(
539 self.heat_duty.keys(),
540 [
541 iscale.get_scaling_factor(v, default=0.01)
542 for v in self.heat_duty.values()
543 ],
544 )
545 )
547 # Apply scaling to constraints
548 for t, c in self.heat_transfer_equation.items():
549 iscale.constraint_scaling_transform(
550 c, sf_cop * sf_work[t], overwrite=False
551 )
553 for t, c in self.unit_heat_balance.items():
554 iscale.constraint_scaling_transform(
555 c, sf_q[t], overwrite=False
556 )
558 for t, c in self.delta_temperature_in_equation.items():
559 iscale.constraint_scaling_transform(c, sf_dT_in[t], overwrite=False)
561 for t, c in self.delta_temperature_out_equation.items():
562 iscale.constraint_scaling_transform(c, sf_dT_out[t], overwrite=False)
564 # COP equation scaling
565 iscale.constraint_scaling_transform(
566 self.cop_equation, sf_cop, overwrite=False
567 )