Expanding Property Packages

[dimi3saad]

Hello!

I have been trying to develop my own property package (cubic EoS). I initially began following the IDAES documentation, only to find that a lot of it is very close to the generic BT one available. Therefore, I opted to modify the already existing BT one to include my components of interest (CO2, H2, etc.). I thus gathered all the required parameters for these species.

Initially I had errors with allowing CO2 and H2 to exist in anything other than the vapor phase. When I removed this, the error changed, and now the final error statement goes like:
ValueError: math domain error

I have a few questions:
1- Is modifying the generic property package a good idea, or is there a better way?
2- What could this error mean?

Thanks in advance for your help. I am willing and ready to share with you the files I am using.

[andrewlee94]

Hi @dimi3saad, a math domain error indicates that there is an issue evaluating on of the numerical expression in the model, such as a very large number or a division by zero. Given what you said above, my guess is that you tried to make CO2 or H2 a component in both phases and that the vapour pressure of those components is causing the problem. At the temperatures and pressures of interest for a benzene-toluene system, CO2 and H2 are both supercritical and cannot be modeled using Raoults' Law - this is why they were set as vapour phase only components.

If you want to model dissolution of thee species in the liquid phase, you would need to use Henry's Law instead. However, the existing property package do not support Henry's Law yet, so you would need to custom code this yourself for now (we are looking at adding Henry's Law to the generic property framework in the upcoming month).

If you do chose to write your own property package, I would suggest not starting with the generic framework - there is far to much complexity in there to make it easy to modify. Instead, you could look at the older Cubic EoS model (https://github.com/IDAES/idaes-pse/tree/main/idaes/generic_models/properties/cubic_eos) as a starting point.

[dimi3saad]

Thanks for your reply @andrewlee94 !

I have no interest in modeling CO2 and H2 as liquids (feed is at atmospheric temperature and pressure). Because of that, I have assigned in the package I developed that CO2 and H2 only exist in the vapor phase. Doing so however still leads to the math domain error.

Do you think it could be caused by something else in the package? Thank you in advance.

[andrewlee94]

We would need more information on the code and actual error to give you anything more specific. Given that this is a math domain error and not a solver error, it implies that this is occurring during evaluation of an expression during code execution. This could theoretically happen for any expression, so there is no way to narrow it down without more information.

Could you provide us with the full stack trace for the error message, as well as the files you are using?

[dimi3saad]

Definitely. I have attached the files here in a ZIP file. The .py file is the actual property package code, which should be put in the cubic EoS folder on your PC. The 2nd (.ipynb) is the one I was running, which is a modified version of the heater example.

Thank you so much. Please let me know if you need anything else.
CO2_propertypackage.zip

[andrewlee94]

@dimi3saad I think you sent the wrong Jupyter notebook; the one you attached appears to be the standard Heater example and uses the benzene-toluene property package instead of your custom one. Also, could you please include the full error traceback so that we can see what is going wrong for you and where - this can be important in case we cannot reproduce the error ourselves.

I had a quick look at your property package and did not see any obvious issues. However, I would suggest you consider using the newer Generic Properties Framework instead; it is much easier to use and hopefully more robust.

[dimi3saad]

@andrewlee94 apologies for that, I have 2 IDAES folders and I got mixed up again. Attached is the correct one now. I double checked this time. Regarding the traceback, here is what I get:
_**---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
in
1 # Initialize the flowsheet, and set the output at WARNING
----> 2 m.fs.heater.initialize(outlvl=idaeslog.WARNING)
3 # From the output it can be inferred that since there are no errors or warnings encountered during initialization, nothing is displayed

~/anaconda3/envs/epse/lib/python3.8/site-packages/idaes/core/unit_model.py in initialize(blk, state_args, outlvl, solver, optarg)
647 # ---------------------------------------------------------------------
648 # Initialize control volume block
--> 649 flags = blk.control_volume.initialize(
650 outlvl=outlvl,
651 optarg=optarg,

~/anaconda3/envs/epse/lib/python3.8/site-packages/idaes/core/control_volume0d.py in initialize(blk, state_args, outlvl, optarg, solver, hold_state)
1364
1365 # Initialize state blocks
-> 1366 in_flags = blk.properties_in.initialize(
1367 outlvl=outlvl,
1368 optarg=optarg,

~/anaconda3/envs/epse/lib/python3.8/site-packages/idaes/generic_models/properties/cubic_eos/cubic_prop_pack.py in initialize(blk, state_args, state_vars_fixed, hold_state, outlvl, solver, optarg)
434 if cons_count > 0:
435 with idaeslog.solver_log(solve_log, idaeslog.DEBUG) as slc:
--> 436 res = solve_indexed_blocks(opt, [blk], tee=slc.tee)
437 else:
438 res = ""

~/anaconda3/envs/epse/lib/python3.8/site-packages/idaes/core/util/initialization.py in solve_indexed_blocks(solver, blocks, **kwds)
204
205 # Solve temporary Block
--> 206 results = solver.solve(tmp, **kwds)
207
208 finally:

~/anaconda3/envs/epse/lib/python3.8/site-packages/pyomo/opt/base/solvers.py in solve(self, *args, **kwds)
573 initial_time = time.time()
574
--> 575 self._presolve(*args, **kwds)
576
577 presolve_completion_time = time.time()

~/anaconda3/envs/epse/lib/python3.8/site-packages/pyomo/opt/solver/shellcmd.py in _presolve(self, *args, **kwds)
196 self._define_signal_handlers = kwds.pop('use_signal_handling',None)
197
--> 198 OptSolver._presolve(self, *args, **kwds)
199
200 #

~/anaconda3/envs/epse/lib/python3.8/site-packages/pyomo/opt/base/solvers.py in _presolve(self, *args, **kwds)
670 write_start_time = time.time()
671 (self._problem_files, self._problem_format, self._smap_id) =
--> 672 self._convert_problem(args,
673 self._problem_format,
674 self._valid_problem_formats,

~/anaconda3/envs/epse/lib/python3.8/site-packages/pyomo/opt/base/solvers.py in _convert_problem(self, args, problem_format, valid_problem_formats, **kwds)
740 # Otherwise, we try to convert the object explicitly.
741 #
--> 742 return convert_problem(args,
743 problem_format,
744 valid_problem_formats,

~/anaconda3/envs/epse/lib/python3.8/site-packages/pyomo/opt/base/convert.py in convert_problem(args, target_problem_type, valid_problem_types, has_capability, **kwds)
103 tmpkw = kwds
104 tmpkw['capabilities'] = has_capability
--> 105 problem_files, symbol_map = converter.apply(*tmp, **tmpkw)
106 return problem_files, ptype, symbol_map
107

~/anaconda3/envs/epse/lib/python3.8/site-packages/pyomo/solvers/plugins/converter/model.py in apply(self, *args, **kwds)
185 else:
186 (problem_filename, symbol_map_id) =
--> 187 instance.write(
188 filename=problem_filename,
189 format=args[1],

~/anaconda3/envs/epse/lib/python3.8/site-packages/pyomo/core/base/block.py in write(self, filename, format, solver_capability, io_options)
1809 if solver_capability is None:
1810 def solver_capability(x): return True
-> 1811 (filename, smap) = problem_writer(self,
1812 filename,
1813 solver_capability,

pyomo/repn/plugins/ampl/ampl_.pyx in pyomo.repn.plugins.ampl.ampl_.ProblemWriter_nl.call()

pyomo/repn/plugins/ampl/ampl_.pyx in pyomo.repn.plugins.ampl.ampl_.ProblemWriter_nl.call()

pyomo/repn/plugins/ampl/ampl_.pyx in pyomo.repn.plugins.ampl.ampl_.ProblemWriter_nl.call()

pyomo/repn/plugins/ampl/ampl_.pyx in pyomo.repn.plugins.ampl.ampl_.ProblemWriter_nl._print_model_NL()

pyomo/repn/standard_repn.pyx in pyomo.repn.standard_repn.generate_standard_repn()

pyomo/repn/standard_repn.pyx in pyomo.repn.standard_repn._generate_standard_repn()

pyomo/repn/standard_repn.pyx in pyomo.repn.standard_repn._collect_sum()

pyomo/repn/standard_repn.pyx in pyomo.repn.standard_repn._collect_standard_repn()

pyomo/repn/standard_repn.pyx in pyomo.repn.standard_repn._collect_nonl()

pyomo/core/expr/numeric_expr.pyx in pyomo.core.expr.numeric_expr.UnaryFunctionExpression._apply_operation()

ValueError: math domain error**_

I will check the Generic Properties Framework as well. Thank you so much for your support!

heater_CO2.ipynb.zip

[andrewlee94]

Hi again,

I just had a look through the lines of code that were failing, and the problem lies in the Cubic EoS model. It was written before we included support for non-condensable components and assumes that all components exist in the liquid phase. Thus, it includes the non-condensable species when trying to calculate initial values for the bubble and dew points - as these are well above their critical points, the resulting vapour pressure is extremely large which results in the math domain errors.

So, what are the possible fixes:

1. The Generic Property Framework has addressed this issue and should work for this system.
2. You could try setting the A parameter in the Antoine equation for the non-condensable components to a very small number. This might be enough to make it work (although I cannot guarantee that). The Antoine equation is only used during initialisation, so it will not affect the final solution at all.

I will also talk to the other developers and see if we want to update the Cubic EoS model, although I think we are gradually trying to deprecate it.

[dimi3saad]

@andrewlee94 Thank you! This was very insightful. I will try these and hopefully everything works. Have a great day!