implement better bubble point with lle

Author: yoelcortes

thermosteam/_stream.py

@@ -2338,7 +2338,7 @@ def bubble_point_at_T(self, T: Optional[float] = None, IDs: Optional[Sequence[st
         z = self.get_normalized_mol(bp.IDs)
         return bp(z, T=T or self.T)
 
-    def bubble_point_at_P(self, P: Optional[float] = None, IDs: Optional[Sequence[str]] = None, vlle=False):
+    def bubble_point_at_P(self, P: Optional[float] = None, IDs: Optional[Sequence[str]] = None, lle=False):
         """
         Return a BubblePointResults object with all data on the bubble point at constant pressure.
 
@@ -2356,27 +2356,9 @@ def bubble_point_at_P(self, P: Optional[float] = None, IDs: Optional[Sequence[st
         BubblePointValues(T=357.14, P=101325, IDs=('Water', 'Ethanol'), z=[0.836 0.164], y=[0.492 0.508])
 
         """
-        if vlle:
-            phase = self.phase
-            try:
-                # TODO: Optimize
-                vlle = self.vlle
-                vlle.top_chemical = self.chemicals.IDs[self.imol.data.sum(axis=0).argmax()]
-                if P is None: P = self.P
-                vlle.bubble_point_at_P(P)
-                BP = self.get_bubble_point(IDs)
-                x = self.imol['L', BP.IDs]
-                bp = BP(P=P, z=x / x.sum())
-                z = self.imol[bp.IDs]
-                z /= z.sum()
-                bp.z = z / z.sum() 
-            finally:
-                self.phase = phase
-            return bp
-        else:
-            BP = self.get_bubble_point(IDs)
-            z = self.get_normalized_mol(BP.IDs)
-            return BP(z, P=P or self.P)
+        BP = self.get_bubble_point(IDs)
+        z = self.get_normalized_mol(BP.IDs)
+        return BP(z, P=P or self.P, lle=lle)
 
     def dew_point_at_T(self, T: Optional[float] = None, IDs: Optional[Sequence[str]] = None):
         """

thermosteam/equilibrium/bubble_point.py

@@ -11,6 +11,7 @@
 import flexsolve as flx
 from .fugacity_coefficients import IdealFugacityCoefficients
 from .domain import vle_domain
+from .lle import solve_lle_mol
 from ..exceptions import InfeasibleRegion
 from .. import functional as fn
 from .._settings import settings
@@ -156,6 +157,18 @@ def _T_error(self, T, P, z_over_P, z_norm, y):
         y[:] = solve_y(y_phi, self.phi, T, P, y)
         return 1. - y.sum()
 
+    def _T_error_lle(self, T, P, z, y, x):
+        if T <= 0: raise InfeasibleRegion('negative temperature')
+        Psats = np.array([i(T) for i in self.Psats], dtype=float)
+        mol = solve_lle_mol(self.gamma, z, T, P, sample=x)
+        x[:] = mol / mol.sum() if mol.any() else z
+        y_phi =  (x / P
+                  * Psats
+                  * self.gamma(x, T, P) 
+                  * self.pcf(T, P, Psats))
+        y[:] = solve_y(y_phi, self.phi, T, P, y)
+        return 1. - y.sum()
+    
     def _P_error(self, P, T, z_Psat_gamma, Psats, y):
         if P <= 0: raise InfeasibleRegion('negative pressure')
         y_phi = z_Psat_gamma * self.pcf(T, P, Psats) / P
@@ -168,6 +181,14 @@ def _P_error_dep(self, P, T, z, Psats, z_Psats, y):
         y[:] = solve_y(y_phi, self.phi, T, P, y)
         return 1. - y.sum()
 
+    def _P_error_lle(self, P, T, z, Psats, y, x):
+        if P <= 0: raise InfeasibleRegion('negative pressure')
+        mol = solve_lle_mol(self.gamma, z, T, P, sample=x)
+        x[:] = mol / mol.sum() if mol.any() else z
+        y_phi = Psats * x * self.gamma(x, T, P) * self.pcf(T, P, Psats) / P
+        y[:] = solve_y(y_phi, self.phi, T, P, y)
+        return 1. - y.sum()
+    
     def _T_error_reactive(self, T, P, z, dz, y, x, liquid_conversion):
         if T <= 0: raise InfeasibleRegion('negative temperature')
         dz[:] = liquid_conversion(z, T, P, 'l')
@@ -217,21 +238,21 @@ def _Py_ideal(self, z_Psat_gamma_pcf):
         y = z_Psat_gamma_pcf / P
         return P, y
 
-    def __call__(self, z, *, T=None, P=None, liquid_conversion=None):
+    def __call__(self, z, *, T=None, P=None, liquid_conversion=None, lle=False):
         z = np.asarray(z, float)
         if T:
             if P: raise ValueError("may specify either T or P, not both")
-            P, *args = self.solve_Py(z, T, liquid_conversion)
+            P, *args = self.solve_Py(z, T, liquid_conversion, lle)
         elif P:
-            T, *args = self.solve_Ty(z, P, liquid_conversion)
+            T, *args = self.solve_Ty(z, P, liquid_conversion, lle)
         else:
             raise ValueError("must specify either T or P")
         if liquid_conversion:
             return ReactiveBubblePointValues(T, P, self.IDs, z, *args)
         else:
             return BubblePointValues(T, P, self.IDs, z, *args)
 
-    def solve_Ty(self, z, P, liquid_conversion=None, guess=None):
+    def solve_Ty(self, z, P, liquid_conversion=None, lle=False):
         """
         Bubble point at given composition and pressure.
 
@@ -270,11 +291,15 @@ def solve_Ty(self, z, P, liquid_conversion=None, guess=None):
             y = z.copy()
             return T, fn.normalize(y)
         elif liquid_conversion is None:
-            f = self._T_error
-            z_norm = z / z.sum()
-            z_over_P = z/P
+            z = z / z.sum()
+            z_over_P = z / P
             T_guess, y = self._Ty_ideal(z_over_P)
-            args = (P, z_over_P, z_norm, y)
+            if lle:
+                f = self._T_error_lle
+                args = (P, z, y, y.copy())
+            else:
+                f = self._T_error
+                args = (P, z_over_P, z, y)
             try:
                 T = flx.aitken_secant(f, T_guess, T_guess + 1e-3,
                                       self.T_tol, 5e-12, args,
@@ -311,7 +336,7 @@ def solve_Ty(self, z, P, liquid_conversion=None, guess=None):
                                          maxiter=self.maxiter)
             return T, dz, fn.normalize(y), x
 
-    def solve_Py(self, z, T, liquid_conversion=None):
+    def solve_Py(self, z, T, liquid_conversion=None, lle=False):
         """
         Bubble point at given composition and temperature.
 
@@ -353,17 +378,23 @@ def solve_Py(self, z, T, liquid_conversion=None):
             if T > self.Tmax: T = self.Tmax
             elif T < self.Tmin: T = self.Tmin
             Psats = np.array([i(T) for i in self.Psats])
-            z_norm = z / z.sum()
+            z = z / z.sum()
             if self.gamma.P_dependent:
                 f = self._P_error_dep
-                z_Psats = z_norm * Psats
+                z_Psats = z * Psats
                 P_guess, y = self._Py_ideal(z_Psats)
-                args = (T, z_norm, Psats, z_Psats, y)
+                args = (T, z, Psats, z_Psats, y)
             else:
-                z_Psat_gamma = z_norm * Psats * self.gamma(z_norm, T, 101325)
-                P_guess, y = self._Py_ideal(z_Psat_gamma)
-                f = self._P_error
-                args = (T, z_Psat_gamma, Psats, y)
+                if lle:
+                    z_Psat_gamma = z * Psats * self.gamma(z, T, 101325)
+                    P_guess, y = self._Py_ideal(z_Psat_gamma)
+                    f = self._P_error_lle
+                    args = (P, T, z, Psats, y, y.copy())
+                else:
+                    z_Psat_gamma = z * Psats * self.gamma(z, T, 101325)
+                    P_guess, y = self._Py_ideal(z_Psat_gamma)
+                    f = self._P_error
+                    args = (T, z_Psat_gamma, Psats, y)
             try:
                 P = flx.aitken_secant(f, P_guess, P_guess-1, self.P_tol, 1e-9,
                                       args, checkiter=False, maxiter=self.maxiter)

thermosteam/equilibrium/fugacities.py

@@ -56,7 +56,7 @@ def unweighted(self, x, T, P=101325.):
         return self.gamma(x, T, P) * self.pcf(T, P, Psats) * Psats
 
     def __call__(self, x, T, P=101325., reduce=False):
-        f_reduced = x * self.gamma(x, T)
+        f_reduced = x * self.gamma(x, T, P)
         if reduce: return f_reduced
         Psats = np.array([i.Psat(T) for i in self.chemicals], dtype=float)
         return f_reduced * Psats * self.pcf(T, P, Psats)

thermosteam/equilibrium/lle.py

@@ -12,10 +12,11 @@
 from .equilibrium import Equilibrium
 from ..exceptions import NoEquilibrium
 from .binary_phase_fraction import phase_fraction
-from .tangential_plane_stability import TangentPlaneStabilityAnalysis
+from .tangential_plane_stability import lle_tangential_plan_analysis
 from scipy.optimize import shgo, differential_evolution
 import flexsolve as flx
 import numpy as np
+from scipy.optimize import minimize
 
 __all__ = ('LLE', 'LLECache')
 
@@ -117,6 +118,34 @@ def pseudo_equilibrium(K, phi, z, T, n, f_gamma, gamma_args, inner_loop_options,
     else:
         return z
 
+# Light weight
+def solve_lle_mol(gamma, z, T, P, sample=None, phi=1):
+    n = z.size
+    stability = lle_tangential_plan_analysis(gamma, z, T, P, sample=sample)
+    if stability.unstable:
+        y = stability.candidate
+        y[y < 1e-64] = 1e-64
+        if stability.sample_unstable: 
+            phi = min(0.99 * (z / y).min(), phi)
+        else:
+            phi = 0.99 * (z / y).min()
+        x = z - phi * y
+        x /= x.sum()
+        K = gamma(y, T) / gamma(x, T)
+    else:
+        return z
+    K[K <= 0] = 1e-9
+    mol = pseudo_equilibrium(
+        K, phi, z, T, n, gamma.f, gamma.args, 
+        LLE.pseudo_equilibrium_inner_loop_options,
+        LLE.pseudo_equilibrium_outer_loop_options,
+    )
+    if mol.any():
+        return mol / mol.sum()
+    else:
+        return z
+    
+
 class LLE(Equilibrium, phases='lL'):
     """
     Create a LLE object that performs liquid-liquid equilibrium when called.
@@ -328,8 +357,7 @@ def solve_lle_liquid_mol(self, z, T, lle_chemicals, single_loop):
                 sample = x0
             else:
                 sample = None
-            TPSA = TangentPlaneStabilityAnalysis('lL', lle_chemicals, thermo=self.thermo)
-            stability = TPSA(z, T, 101325, sample=sample)
+            stability = lle_tangential_plan_analysis(gamma, z, T, 101325, sample=sample)
             if stability.unstable:
                 y = stability.candidate
                 y[y < 1e-64] = 1e-64

thermosteam/equilibrium/tangential_plane_stability.py

@@ -63,7 +63,6 @@ def edge_points_simplex_masked(z: np.ndarray,
     t = 0.5 * (1 - np.cos(np.pi * k / (points_per_edge - 1)))
 
     pts = []
-
     # edges only among the selected components
     for i, j in combinations(active, 2):
         for tau in t:
@@ -90,6 +89,75 @@ def edge_points_simplex_masked(z: np.ndarray,
 
     return np.array(pts)
 
+# Light weight
+def lle_tangential_plan_analysis(gamma, z, T, P, sample=None):
+    MW = np.array([i.MW for i in gamma.chemicals])
+    logfz = np.log(z * gamma(z, T, P) + 1e-30)
+    
+    def objective(w, T, P, logfz, softmax=False):
+        if softmax:
+            w = np.exp(w - np.max(w)) # Softmax for unconstrained optimization
+            w /= w.sum()
+        return np.dot(w, np.log(w * gamma(w, T, P) + 1e-30) - logfz)
+    
+    args = (T, P, logfz)
+    if sample is None:
+        best_val = np.inf
+        best_result = None
+    else:
+        best_result = sample
+        best_val = objective(sample, *args)
+        result = minimize(
+            objective, 
+            best_result, 
+            method="L-BFGS-B", 
+            options=dict(maxiter=20),
+            args=(*args, True),
+        )
+        value = result.fun
+        if value < best_val:
+            w = result.x
+            w = np.exp(w - np.max(w)) # Softmax for unconstrained optimization
+            w /= w.sum()
+            best_result = w
+            best_val = value
+        unstable = best_val < 0
+        if unstable:
+            return StabilityReport(
+                unstable=unstable,
+                candidate=sample,
+                tpd=best_val,
+                sample_unstable=True,
+            )
+    w = z * MW
+    w /= w.sum()
+    samples = edge_points_simplex_masked(w)
+    samples /= MW
+    samples /= samples.sum(axis=1, keepdims=True)
+    for sample in samples:
+        value = objective(sample, *args)
+        if value < best_val:
+            best_val = value
+            best_result = sample
+            result = minimize(
+                objective, 
+                best_result, 
+                method="L-BFGS-B", 
+                options=dict(maxiter=20),
+                args=(*args, True),
+            )
+            value = result.fun
+            if value < best_val:
+                w = result.x
+                w = np.exp(w - np.max(w)) # Softmax for unconstrained optimization
+                w /= w.sum()
+                best_result = w
+                best_val = value
+    return StabilityReport(
+        unstable=(best_val < 0),
+        candidate=best_result,
+        tpd=best_val
+    )
 
 class TangentPlaneStabilityAnalysis:
     __slots__ = ('fugacity_models', 'MW')
@@ -175,6 +243,6 @@ def __call__(self, z, T, P, reference_phase='l', potential_phase='L', sample=Non
 
 
 if __name__ == '__main__':
-    tmo.settings.set_thermo(['Water', 'Octanol'])
+    tmo.settings.set_thermo(['Water', 'Octanol', 'Ethanol'])
     TPSA = TangentPlaneStabilityAnalysis('lL', tmo.settings.chemicals)
-    report = TPSA(np.array([0.2, 0.8]), 298.15, 101325)
+    report = TPSA(np.array([0.25, 0.8, 0.05]), 298.15, 101325)