Add new Complementarity formualtion for VLE with cubic EoSs #1397
Conversation
andrewlee94
added
enhancement
| def identify_VL_component_list(blk, phase_pair): | ||
| """ | ||
| Identify liquid and vapor phases and which components are in VL equilibrium | ||
|
|
||
| Args: | ||
| blk: StateBlock of interest | ||
| phase_pair: 2-tuple of Phases in equilibrium | ||
|
|
||
| Returns: | ||
| Lists of component names for: | ||
| * liquid Phase object | ||
| * vapor Phase object | ||
| * components using Raoult's Law | ||
| * components using Henry's Law | ||
| * liquid only components, | ||
| * vapor only components | ||
|
|
||
| """ |
There was a problem hiding this comment.
How does this function handle systems with multiple liquid phases? How does this function handle things with solid phases? Now that it's no longer private, we need to consider these questions.
There was a problem hiding this comment.
This is only for VLE systems, so is not expected to work in those cases. I have added an error for cases where a non-VL pair is passed as the argument.
| # Use lowest component temperature_crit as starting point | ||
| # Starting high and moving down generally works better, | ||
| # as it under-predicts next step due to exponential form of | ||
| # Psat. | ||
| # Subtract 1 to avoid potential singularities at Tcrit | ||
| Tbub0 = ( | ||
| min(blk.params.get_component(j).temperature_crit.value for j in raoult_comps) | ||
| - 1 |
There was a problem hiding this comment.
Should we convert 1 from Kelvin to the temperature_units of params?
There was a problem hiding this comment.
In this case, 1 is just an arbitrary value to move away from the exact critical value as some correlations are singular if T=Tc[j].
There was a problem hiding this comment.
Yeah, but in edge cases we could end up with a negative temperature (for example, if the user decided to work with temperature reduced by Tc).
| if l_phase is None or v_phase is None: | ||
| raise ConfigurationError( | ||
| f"{b.params.name} Generic Property Package phase pair {phase_pair[0]}-{phase_pair[1]} " | ||
| "was set to use Smooth VLE formulation, however this is not a vapor-liquid pair." | ||
| ) |
There was a problem hiding this comment.
This error should probably be thrown by identify_VL_component_list.
There was a problem hiding this comment.
A catch has been added to the utility method, but this can be retained as well.
| s = Var( | ||
| b.params.phase_list, | ||
| initialize=0.0, | ||
| bounds=(0, None), | ||
| doc="Slack variable for equilibrium temperature", | ||
| units=pyunits.dimensionless, | ||
| ) | ||
| b.add_component("s" + suffix, s) | ||
|
|
||
| # Equilibrium temperature | ||
| def rule_teq(b): | ||
| if b.params.get_phase(phase_pair[0]).is_vapor_phase(): | ||
| vapor_phase = phase_pair[0] | ||
| liquid_phase = phase_pair[1] | ||
| else: | ||
| vapor_phase = phase_pair[1] | ||
| liquid_phase = phase_pair[0] | ||
| return ( | ||
| b._teq[phase_pair] | ||
| - b.temperature | ||
| - s[vapor_phase] * t_units | ||
| + s[liquid_phase] * t_units | ||
| == 0 |
There was a problem hiding this comment.
Why make the slack variable s dimensionless when we're multiplying it by t_units here? That also implies it should be scaled.
| eps = Param( | ||
| default=1e-04, | ||
| mutable=True, | ||
| doc="Smoothing parameter for complementarities", | ||
| units=f_units, | ||
| ) | ||
| b.add_component("eps" + suffix, eps) |
There was a problem hiding this comment.
Does eps need to be adjusted by the user?
There was a problem hiding this comment.
Yes, users should adjust eps as necessary and this will be in the docs (the old formulation was the same).
| def rule_temperature_slack_complementarity(b, p): | ||
| flow_phase = b.flow_mol_phase[p] | ||
|
|
||
| return smooth_min(s[p] * f_units, flow_phase, eps) == 0 |
There was a problem hiding this comment.
Okay, we're treating s as having different units depending on the context. That's not ideal for scaling.
There was a problem hiding this comment.
Agreed - if you have any ideas on how to handle this better it would be appreciated. For now I have just focused on getting the code running.
There was a problem hiding this comment.
My instinct here is to just give s[p] units of temperature per flowrate, then scale this equation like temperature. However, it needs to be scaled with a "temperature difference" scaling factor, not an "absolute temperature" scaling factor.
There was a problem hiding this comment.
Based on some testing, scaling of s, g and epsilon depends on the larger of T and flow, thus it is hard to know what to do with these. I am inclined to leave these as unit-less for that reason, and any scaling routine can handle the necessary logic behind the scenes (as I doubt users will care enough to understand what is happening).
idaes/models/properties/modular_properties/phase_equil/smooth_VLE_2.py
Outdated
Show resolved
Hide resolved
andrewlee94
changed the title
Codecov ReportAttention: Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #1397 +/- ##
==========================================
- Coverage 77.63% 77.52% -0.12%
==========================================
Files 391 392 +1
Lines 64391 64510 +119
Branches 14264 14268 +4
==========================================
+ Hits 49989 50010 +21
- Misses 11830 11935 +105
+ Partials 2572 2565 -7 ☔ View full report in Codecov by Sentry. |
|
@viiibhav @jghouse88 Not sure if either of you still receive these, but input and review on this PR would be welcome. |
idaes/models/properties/modular_properties/base/tests/test_generic_property_integration.py
Outdated
Show resolved
Hide resolved
Will take a look this week if that works. The best test would be to see if the paper examples work with this pr and see if we get similar results. |
|
@jghouse88 I have tested this with the existing tests in IDAES, swapping the new formulation for the old one. All but two tests pass with the new formulation, but the two that fail appear to need further investigation - the temperature sweep test is failing to converge to the correct solution (it gets stuck somewhere). |
Codecov ReportAttention: Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #1397 +/- ##
==========================================
+ Coverage 77.80% 77.83% +0.02%
==========================================
Files 393 394 +1
Lines 64797 64926 +129
Branches 14390 14398 +8
==========================================
+ Hits 50413 50532 +119
- Misses 11794 11811 +17
+ Partials 2590 2583 -7 ☔ View full report in Codecov by Sentry. |
docs/explanations/components/property_package/general/pe/smooth_vle2.rst
Outdated
Show resolved
Hide resolved
docs/explanations/components/property_package/general/pe/smooth_vle2.rst
Outdated
Show resolved
Hide resolved
docs/explanations/components/property_package/general/pe/smooth_vle2.rst
Outdated
Show resolved
Hide resolved
|
|
||
| Each complementarity requires a smoothing parameter, named :math:`\epsilon_T` and :math:`\epsilon_Z` for the temperature and cubic root constraints respectively. Within the IDAES model, these are rendered as ``eps_t_phase1_phase2`` and ``eps_z_phase1_phase2``, where ``phase1`` and ``phase2`` are the names assigned to the liquid and vapor phases in the property package (order will depend on the order these are declared). | ||
|
|
||
| The tractability of the VLE problem depends heavily upon the values chosen for :math:`\epsilon_T` and :math:`\epsilon_Z`, with larger values resulting in smoother transitions at the phase boundaries (and thus increased tractability) at the expense of decreased accuracy near these points. It is recommended that users employ a 2-stage approach to solving these problems, starting with a larger value of :math:`\epsilon_T` and :math:`\epsilon_Z` initially to determine which region the solution lies in, followed by a second solve using smaller values to refine the solution. |
There was a problem hiding this comment.
"It is recommended that users employ a 2-stage approach to solving these problems, starting with a larger value of :math:\epsilon_T and :math:\epsilon_Z initially to determine which region the solution lies in, followed by a second solve using smaller values to refine the solution."
How is this going to translate to a typical user of IDAES? Are the smoothing parameters global or stuck on state blocks? How can the user find them?
There was a problem hiding this comment.
They are local, as they potentially need to vary unit-to-unit. Their names are show in the docs.
idaes/models/properties/modular_properties/phase_equil/__init__.py
Outdated
Show resolved
Hide resolved
idaes/models/properties/modular_properties/phase_equil/smooth_VLE_2.py
Outdated
Show resolved
Hide resolved
idaes/models/properties/modular_properties/phase_equil/smooth_VLE_2.py
Outdated
Show resolved
Hide resolved
| s = Var( | ||
| vl_phase_set, | ||
| initialize=EPS_INIT, | ||
| bounds=(0, None), | ||
| doc="Slack variable for equilibrium temperature", | ||
| units=t_units, | ||
| ) |
There was a problem hiding this comment.
Do we have any sort of scaling for these slack variables, or is the intention to let the autoscaler do it? This is a case where the results from the autoscaler is going to be highly dependent on initialization, and since the slack variable will frequently be near zero, you may not get any scaling at all.
| try: | ||
| teq_cons = getattr(b, "_teq_constraint" + suffix) | ||
| # pylint: disable-next=protected-access | ||
| iscale.set_scaling_factor(b._teq[phase_pair], sf_T) | ||
| iscale.constraint_scaling_transform(teq_cons, sf_T, overwrite=False) |
There was a problem hiding this comment.
I guess we're not scaling the slacks then. If we're giving them temperature units, we might as well scale them like temperature. It would be a good first guess, at any rate, and I think the autoscalers will faceplant.
There was a problem hiding this comment.
I am not making any changes to the legacy scaling methods at this point - this code was just copied over from the old implementation. We can think about what is needed for this in the new tools once we get there.
There was a problem hiding this comment.
Very well, but I do expect this case will need heuristic scaling.
| if self.temperature.value is not None: | ||
| t_value = value(self.temperature) | ||
| else: | ||
| t_value = None |
There was a problem hiding this comment.
If this is None, does the model break or is there a catch for this case?
There was a problem hiding this comment.
This just affects whether we can get initial guesses for Teq. Whether the model will solve or not is a separate issue at this point; the goal here is for the code to work without encountering a terminal error.
| "temperature_ref": (298.15, pyunits.K), | ||
| "parameter_data": { | ||
| "PR_kappa": { | ||
| ("foo", "bar"): 0.000, |
There was a problem hiding this comment.
| m = ConcreteModel() | ||
| m.params = DummyParameterBlock( | ||
| components={ | ||
| "a": {}, |
There was a problem hiding this comment.
Are these supposed to be empty?
There was a problem hiding this comment.
Yes - this is a unit test. We need to define the components so the code runs but don't care about any actual properties.
Replaces #977
Summary/Motivation:
This PR adds an implementation of the improved Smooth VLE formulation originally proposed in #977. This implementation preserves the old implementation for backward compatibility.
Changes proposed in this PR:
SmoothVLE2to implement improved formulationLegal Acknowledgement
By contributing to this software project, I agree to the following terms and conditions for my contribution: