fix HXN pinch analysis: energy-consistent problem table

biosteam/facilities/hxn/hxn_synthesis.py

@@ -10,11 +10,13 @@
 
 @author: sarangbhagwat
 """
+from collections import namedtuple
 import numpy as np
 import biosteam as bst
 from warnings import warn
 
-__all__ = ('StreamLifeCycle', 'synthesize_network')
+__all__ = ('StreamLifeCycle', 'ProblemTable', 'problem_table',
+           'synthesize_network')
 
 class LifeStage:
 
@@ -124,8 +126,179 @@ def get_sorted_life_cycle(self):
         return self.life_cycle
 
 
-def temperature_interval_pinch_analysis(hus, 
-                                        T_min_app=10, 
+ProblemTable = namedtuple(
+    'ProblemTable',
+    ['Ts', 'interval_H', 'point_H', 'residual',
+     'hot_util_load', 'cold_util_load', 'pinch_T']
+)
+
+def _stream_H_at_boundaries(stream_in, H_in, H_out, T_lo, T_hi, Ts, shift,
+                            stream_label):
+    """
+    Enthalpies [kJ/hr] of one monotone stream at the grid boundaries
+    `Ts` (shifted scale, descending, all within [T_lo, T_hi]).
+
+    Exact at the stream's own end points (H_in/H_out as given) by
+    *position*: `Ts[0]` and `Ts[-1]` are the stream's own T_hi/T_lo (every
+    monotone stream has `T_hi > T_lo` strictly, so `Ts` always has at least
+    these two entries) and are assigned H_in/H_out directly, without a float
+    comparison. In between, the inlet copy is flashed at the *real*
+    temperature `T + shift` and the result is clipped to
+    [min(H_in, H_out), max(H_in, H_out)] so that a non-equilibrium outlet
+    (e.g. a column reboiler/condenser product) can never inflate an
+    interval. A single copy is walked down the grid so each VLE is
+    warm-started from the previous boundary; `stream_label` (the inlet
+    stream's own ID) identifies the stream in the VLE-failure warning.
+    """
+    assert Ts.size >= 2, (
+        "boundary grid for a monotone stream must include both its own "
+        "end points"
+    )
+    H_lo, H_hi = sorted((H_in, H_out))
+    H_top, H_bottom = (H_in, H_out) if H_in > H_out else (H_out, H_in)
+    Hs = np.empty(Ts.size)
+    Hs[0] = H_top
+    Hs[-1] = H_bottom
+    stream = stream_in.copy()
+    for k in range(1, Ts.size - 1):
+        T = Ts[k]
+        T_real = T + shift
+        try:
+            stream.vle(T=T_real, P=stream.P)
+            H = stream.H
+        except Exception as error:
+            warn(f"could not solve VLE for stream {stream_label!r} at "
+                 f"{T_real:.2f} K ({error!r}); interpolating enthalpy "
+                 "linearly in temperature for the problem table",
+                 RuntimeWarning)
+            # restart the warm start from a clean copy so the failed flash
+            # does not leave `stream` in a bad state for the next boundary
+            stream = stream_in.copy()
+            H = H_lo + (H_hi - H_lo) * (T - T_lo) / (T_hi - T_lo)
+        Hs[k] = min(max(H, H_lo), H_hi)
+    return Hs
+
+def problem_table(streams_inlet, streams_quenched, is_hot, T_min_app):
+    """
+    Energy-consistent problem table (temperature-interval heat cascade).
+
+    Parameters
+    ----------
+    streams_inlet : list[Stream]
+        Inlet stream of each utility heat exchanger.
+    streams_quenched : list[Stream]
+        Corresponding outlet streams, re-flashed at their enthalpy.
+    is_hot : Sequence[bool]
+        True where the stream is cooled.
+    T_min_app : float
+        Minimum approach temperature [K].
+
+    Returns
+    -------
+    ProblemTable
+        Grid temperatures `Ts` (shifted scale, descending), per-stream
+        `interval_H` (N x n-1) and `point_H` (N x n) contributions (+ for
+        hot, - for cold), the cascade `residual` (n) *leaving* each
+        boundary (i.e. after its point loads), `hot_util_load`,
+        `cold_util_load` and the shifted-scale `pinch_T`.
+
+    Notes
+    -----
+    Hot streams are shifted down by `T_min_app`; cold streams are not. For
+    monotone streams the contribution to interval (Ts[k], Ts[k+1]) is
+    sign * (H(Ts[k]) - H(Ts[k+1])) with H evaluated at the real temperature
+    and clipped to [H_in, H_out], so every stream's contributions telescope
+    exactly to sign * |H_out - H_in|. Isothermal streams, and streams whose
+    outlet temperature moves against their duty (a heated stream that exits
+    colder than it entered, e.g. a reboiler outlet at VLE), are point loads
+    at their outlet temperature. The cascade starting from zero hot utility
+    is residual[k] = sum(point_H[:, :k+1]) + sum(interval_H[:, :k]), the
+    heat *leaving* boundary Ts[k]. Feasibility must also hold for the heat
+    *arriving* at Ts[k] before its point loads are applied,
+    arriving[k] = residual[k] - sum(point_H[:, k]), because a source at
+    Ts[k] cannot serve a sink above Ts[k]. The minimum over both flows,
+    min(residual, arriving), fixes the hot utility target,
+    `residual[-1] + hot_util_load` the cold one, and its location the
+    pinch. With the per-stream identity above, hot_util_load -
+    cold_util_load equals the net heating demand.
+
+    Examples
+    --------
+    A threshold problem: 1000 kmol/hr of water cooled 400 -> 300 K supplies
+    every interval of 900 kmol/hr of water heated 300 -> 390 K, so no hot
+    utility is needed and the surplus leaves as cold utility.
+
+    >>> import biosteam as bst
+    >>> from biosteam.facilities.hxn.hxn_synthesis import problem_table
+    >>> bst.settings.set_thermo(['Water'])
+    >>> hot_in = bst.Stream(Water=1000., T=400., P=5e5, phase='l', units='kmol/hr')
+    >>> hot_out = hot_in.copy(); hot_out.vle(T=300., P=5e5)
+    >>> cold_in = bst.Stream(Water=900., T=300., P=5e5, phase='l', units='kmol/hr')
+    >>> cold_out = cold_in.copy(); cold_out.vle(T=390., P=5e5)
+    >>> table = problem_table([hot_in, cold_in], [hot_out, cold_out],
+    ...                       [True, False], 5.)
+    >>> table.Ts
+    array([395., 390., 300., 295.])
+    >>> round(table.hot_util_load, 3)
+    0.0
+    >>> round(table.cold_util_load, -1)
+    1445550.0
+    >>> table.pinch_T
+    395.0
+    """
+    N = len(streams_inlet)
+    is_hot = np.asarray(is_hot, dtype=bool)
+    sign = np.where(is_hot, 1., -1.)
+    shift = np.where(is_hot, T_min_app, 0.)
+    T_in = np.array([s.T for s in streams_inlet])
+    T_out = np.array([s.T for s in streams_quenched])
+    H_in = np.array([s.H for s in streams_inlet])
+    H_out = np.array([s.H for s in streams_quenched])
+    monotone = (sign * (T_in - T_out)) > 0.
+    T_hi = np.where(monotone, np.maximum(T_in, T_out), T_out) - shift
+    T_lo = np.where(monotone, np.minimum(T_in, T_out), T_out) - shift
+    Ts = np.unique(np.concatenate([T_hi, T_lo]))[::-1]
+    n = Ts.size
+    interval_H = np.zeros((N, n - 1))
+    point_H = np.zeros((N, n))
+    for j in range(N):
+        if monotone[j]:
+            idx = np.flatnonzero((Ts <= T_hi[j]) & (Ts >= T_lo[j]))
+            Hs = _stream_H_at_boundaries(streams_inlet[j], H_in[j], H_out[j],
+                                         T_lo[j], T_hi[j], Ts[idx], shift[j],
+                                         streams_inlet[j].ID)
+            interval_H[j, idx[:-1]] = sign[j] * (Hs[:-1] - Hs[1:])
+        else:
+            k = np.searchsorted(-Ts, -T_hi[j])
+            point_H[j, k] = sign[j] * abs(H_out[j] - H_in[j])
+    point_total = point_H.sum(axis=0)
+    residual = np.cumsum(
+        point_total + np.concatenate([[0.], interval_H.sum(axis=0)])
+    )
+    # heat arriving at each boundary, before that boundary's point loads:
+    # a point source at Ts[k] cannot serve sinks above Ts[k], so the cascade
+    # must be non-negative both before and after the point loads
+    arriving = residual - point_total
+    flow = np.minimum(residual, arriving)
+    k_pinch = int(np.argmin(flow))
+    scale = np.abs(H_out - H_in).sum()
+    if -flow[k_pinch] <= 1e-9 * scale:  # threshold problem: no hot utility
+        hot_util_load = 0.
+        k_pinch = 0
+    else:
+        hot_util_load = -flow[k_pinch]
+    cold_util_load = residual[-1] + hot_util_load
+    if cold_util_load < 0.:
+        # only reachable in the threshold branch, by at most 1e-9 * scale:
+        # absorb the rounding into the hot utility so that
+        # hot_util_load - cold_util_load == sum(unit_duty) stays exact
+        hot_util_load -= cold_util_load
+        cold_util_load = 0.
+    return ProblemTable(Ts, interval_H, point_H, residual,
+                        hot_util_load, cold_util_load, Ts[k_pinch])
+
+def temperature_interval_pinch_analysis(hus,
+                                        T_min_app=10,
                                         force_ideal_thermo=False,
                                         sort_hus_by_T=False):
     hx_utils = hus
@@ -150,77 +323,45 @@ def temperature_interval_pinch_analysis(hus,
         ID = 'Util_%s'%i
         stream.ID = 's_%s__%s'%(i,ID)
     N_heating = len(hus_heating)
-    is_cold_stream_index = lambda x: x < N_heating
-    T_in_arr = np.array([stream.T for stream in streams_inlet])
-    T_out_arr = np.array([i.T for i in streams_quenched])
-    adj_T_in_arr = T_in_arr.copy()
-    # adj_T_in_arr[:N_heating] -= T_min_app
-    adj_T_in_arr[N_heating:] -= T_min_app
-    adj_T_out_arr = T_out_arr.copy()
-    # adj_T_out_arr[:N_heating] -= T_min_app
-    adj_T_out_arr[N_heating:] -= T_min_app
-    T_changes_tuples = list(zip(adj_T_in_arr, adj_T_out_arr))
-    all_Ts_descending = [*adj_T_in_arr, *adj_T_out_arr]
-    all_Ts_descending.sort(reverse=True)
-    stream_indices_for_T_intervals =\
-        {(all_Ts_descending[i], all_Ts_descending[i+1]):[]\
-            for i in range(len(all_Ts_descending)-1)}
-    H_for_T_intervals = dict.fromkeys(stream_indices_for_T_intervals, 0)
     cold_indices = list(range(N_heating))
     hot_indices = list(range(N_heating, len(hxs)))
     indices = cold_indices + hot_indices
-    for i in range(len(all_Ts_descending)-1):
-        T_start = all_Ts_descending[i]
-        T_end = all_Ts_descending[i+1]
-        for stream_index in indices:
-            T1, T2 = T_changes_tuples[stream_index]
-            if (T1 >= T_start and T2 <= T_end) or (T2 >= T_start and T1 <= T_end):
-                multiplier = -1 if is_cold_stream_index(stream_index) else 1
-                stream = streams_inlet[stream_index].copy()
-                if stream.T != T_start: stream.vle(T = T_start, P = stream.P)
-                H1 = stream.H
-                try:
-                    stream.vle(T = T_end, P = stream.P)
-                except:
-                    warn(f"could not solve VLE for {repr(stream)} at {repr(hxs[stream_index].owner)}", RuntimeWarning)
-                H2 = stream.H
-                H = multiplier*(H1 - H2)
-                H_for_T_intervals[(T_start, T_end)] += H
-
-    res_H_vector = []
-    prev_res_H = 0
-    for interval, H in H_for_T_intervals.items():
-        res_H_vector.append(prev_res_H + H)
-        prev_res_H = res_H_vector[len(res_H_vector)-1]
-    hot_util_load = - min(res_H_vector)
-    # assert hot_util_load>= 0, 'Hot utility load is negative'
-    if not hot_util_load>=0:
-        warn(f"Hot utility load is negative: {hot_util_load}", RuntimeWarning)
-    # print(hot_util_load)
-    # the lower temperature of the temperature interval for which the res_H is minimum
-    pinch_cold_stream_T = all_Ts_descending[res_H_vector.index(-hot_util_load)+1]
+    T_in_arr = np.array([stream.T for stream in streams_inlet])
+    T_out_arr = np.array([i.T for i in streams_quenched])
+    is_hot = np.zeros(len(hxs), dtype=bool)
+    is_hot[hot_indices] = True
+    table = problem_table(streams_inlet, streams_quenched, is_hot, T_min_app)
+    hot_util_load = table.hot_util_load
+    cold_util_load = table.cold_util_load
+    pinch_cold_stream_T = table.pinch_T
     pinch_hot_stream_T = pinch_cold_stream_T + T_min_app
-    cold_util_load = res_H_vector[len(res_H_vector)-1] + hot_util_load
-    # assert cold_util_load>=0, 'Cold utility load is negative'
-    if not cold_util_load>=0:
-        warn(f"Cold utility load is positive: {cold_util_load}", RuntimeWarning)
+    # Per-stream pinch temperature: where each stream is split between the
+    # hot-side and cold-side network designs. A stream already entirely on
+    # one side of the process pinch (T_in past pinch_cold_stream_T for a
+    # cold stream, or past pinch_hot_stream_T for a hot stream) is not
+    # split; its pinch_T is its own T_in. This clause also catches
+    # non-monotone streams (T_out on the wrong side of T_in for their duty,
+    # e.g. a cold stream whose VLE outlet ends up cooler than it entered):
+    # rather than split their problem_table point-load duty across the
+    # cascade, they get pinch_T = T_in too, so load_duties assigns their
+    # whole duty to a single side (Q_hot_side for a cold stream,
+    # Q_cold_side for a hot one).
     pinch_T_arr = []
     for i in cold_indices:
-        if T_in_arr[i] > pinch_cold_stream_T:
+        if T_in_arr[i] > pinch_cold_stream_T or T_in_arr[i] > T_out_arr[i]:
             pinch_T_arr.append(T_in_arr[i])
         elif T_out_arr[i] < pinch_cold_stream_T:
             pinch_T_arr.append(T_out_arr[i])
         else:
             pinch_T_arr.append(pinch_cold_stream_T)
     for i in hot_indices:
-        if T_in_arr[i] < pinch_hot_stream_T:
+        if T_in_arr[i] < pinch_hot_stream_T or T_in_arr[i] < T_out_arr[i]:
             pinch_T_arr.append(T_in_arr[i])
         elif T_out_arr[i] > pinch_hot_stream_T:
             pinch_T_arr.append(T_out_arr[i])
         else:
             pinch_T_arr.append(pinch_hot_stream_T)
     pinch_T_arr = np.array(pinch_T_arr)
-    # print(pinch_T_arr, hot_util_load, cold_util_load,)
     return pinch_T_arr, hot_util_load, cold_util_load, T_in_arr, T_out_arr,\
            hxs, hot_indices, cold_indices, indices, streams_inlet, hx_utils_rearranged, \
            streams_quenched