From 2d5f77547ca723a3c72ac077f39a7ca6249daf1b Mon Sep 17 00:00:00 2001 From: sarangbhagwat Date: Sat, 22 Aug 2026 16:10:07 -0700 Subject: [PATCH 1/5] add failing energy-consistency tests for the HXN problem table Co-Authored-By: Claude Fable 5 --- tests/test_hxn.py | 140 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 140 insertions(+) diff --git a/tests/test_hxn.py b/tests/test_hxn.py index 758a22e9..eb1cc28e 100644 --- a/tests/test_hxn.py +++ b/tests/test_hxn.py @@ -10,6 +10,7 @@ Tests for the heat exchanger network facility. """ import warnings +import pytest import biosteam as bst import numpy as np from numpy.testing import assert_allclose @@ -90,8 +91,147 @@ def test_energy_balance_error_contributions_ignored_none(): assert len(errors) == N assert HXN.ignored is None +# --------------------------------------------------------------------------- +# Problem-table (pinch) analysis +# --------------------------------------------------------------------------- + +from biosteam.facilities.hxn.hxn_synthesis import ( + temperature_interval_pinch_analysis, problem_table, +) + +def utility_hx(ID, T, P, phase, T_out, **flow): + """A simulated HXutility acting as one process stream (kmol/hr flows).""" + s = bst.Stream(ID + '_in', T=T, P=P, phase=phase, units='kmol/hr', **flow) + hx = bst.HXutility(ID, ins=s, T=T_out, + rigorous=(phase == 'g')) # liquid streams stay liquid (report's case) + hx.simulate() + return hx + +def synthetic_units(): + """Report's 4-stream case: two cold, one condensing hot, one hot liquid.""" + bst.settings.set_thermo(['Water', 'Ethanol'], cache=True) + bst.main_flowsheet.set_flowsheet('test_hxn_synthetic') + return [utility_hx('C1', 300., 101325., 'l', 390., Water=2000.), + utility_hx('C2', 310., 101325., 'l', 345., Water=500., Ethanol=500.), + utility_hx('H1', 352., 101325., 'g', 340., Ethanol=300.), + utility_hx('H2', 420., 5e5, 'l', 320., Water=800.)] + +def heat_utilities(units): + hus = [hx.heat_utilities[0] for hx in units] + hus.sort(key=lambda hu: hu.duty) + return hus + +def pinch_streams(hus): + """Inlet/quenched-outlet stream copies exactly as the synthesizer prepares them.""" + streams_inlet = [hu.unit.ins[0].copy() for hu in hus] + streams_quenched = [hu.unit.outs[0].copy() for hu in hus] + for s in streams_quenched: s.vle(H=s.H, P=s.P) + is_hot = [hu.duty < 0 for hu in hus] + return streams_inlet, streams_quenched, is_hot + +def assert_energy_consistent(hus, T_min_app): + """Invariants of a correct problem table, independent of the synthesizer.""" + unit_duties = np.array([hu.unit_duty for hu in hus]) + table = problem_table(*pinch_streams(hus), T_min_app) + # (a) each stream's grid contributions telescope to its real duty + # (hot: +|dH|, cold: -dH) + per_stream = table.interval_H.sum(axis=1) + table.point_H.sum(axis=1) + assert_allclose(per_stream, -unit_duties, rtol=1e-9) + # (b) targets are non-negative and hot - cold is the net demand + assert table.hot_util_load >= 0. and table.cold_util_load >= 0. + assert_allclose(table.hot_util_load - table.cold_util_load, + unit_duties.sum(), rtol=1e-9) + # (c) a target can never exceed the un-integrated load + assert table.hot_util_load <= unit_duties[unit_duties > 0].sum() * (1 + 1e-9) + assert table.cold_util_load <= -unit_duties[unit_duties < 0].sum() * (1 + 1e-9) + # (d) the public wrapper reports the same targets + pinch_T_arr, hot, cold, *_ = temperature_interval_pinch_analysis(hus, T_min_app) + assert_allclose([hot, cold], [table.hot_util_load, table.cold_util_load], rtol=1e-12) + return table + +def test_problem_table_energy_consistency_doctest_system(): + sys, HXN, feed = build_system() + sys.simulate() + hus = HXN._get_original_heat_utilties() + hus.sort(key=lambda hu: hu.duty) + assert_energy_consistent(hus, 5.) + +@pytest.mark.parametrize('T_min_app', [5., 10., 20.]) +def test_problem_table_energy_consistency_synthetic(T_min_app): + hus = heat_utilities(synthetic_units()) + table = assert_energy_consistent(hus, T_min_app) + # the condensing ethanol stream (352 K in, ~351.4 K dew point) must keep + # its latent heat: hot target well below the un-integrated heating load + heating = sum(hu.unit_duty for hu in hus if hu.unit_duty > 0) + assert table.hot_util_load < 0.5 * heating + +def test_problem_table_two_streams_closed_form(): + bst.settings.set_thermo(['Water', 'Ethanol'], cache=True) + bst.main_flowsheet.set_flowsheet('test_hxn_two_streams') + hot = utility_hx('Hw', 400., 5e5, 'l', 300., Water=1000.) + # Case 1 (threshold problem): the hot stream covers every interval of the + # smaller cold stream -> no hot utility, cold utility = net surplus. + cold = utility_hx('Cs', 300., 5e5, 'l', 390., Water=900.) + hus = heat_utilities([hot, cold]) + table = problem_table(*pinch_streams(hus), 5.) + Q_hot = -hot.heat_utilities[0].unit_duty + Q_cold = cold.heat_utilities[0].unit_duty + assert table.hot_util_load == 0. + assert_allclose(table.cold_util_load, Q_hot - Q_cold, rtol=1e-9) + assert table.pinch_T == table.Ts[0] # no pinch: everything sits below it + # Case 2: the larger cold stream is short of heat everywhere; the cascade + # minimum is at the cold inlet (300 K on the shifted scale) and the hot + # utility is the cold duty minus what the hot stream gives down to 305 K. + cold = utility_hx('Cb', 300., 5e5, 'l', 390., Water=1200.) + hus = heat_utilities([hot, cold]) + table = problem_table(*pinch_streams(hus), 5.) + s = hot.ins[0].copy(); s.vle(T=305., P=s.P) + Q_hot_above_pinch = hot.ins[0].H - s.H + Q_cold = cold.heat_utilities[0].unit_duty + assert_allclose(table.hot_util_load, Q_cold - Q_hot_above_pinch, rtol=1e-9) + assert table.pinch_T == 300. + # remaining hot-stream heat below the pinch leaves as cold utility + assert_allclose(table.cold_util_load, s.H - hot.outs[0].H, rtol=1e-9) + +def test_problem_table_non_monotone_stream_is_point_load(): + # A heated stream whose outlet is colder than its inlet (e.g. a column + # reboiler outlet at VLE): treated as an isothermal load at T_out. + bst.settings.set_thermo(['Water', 'Ethanol'], cache=True) + a = bst.Stream('a', Water=100., T=372., P=101325., phase='l', units='kmol/hr') + b = bst.Stream('b', Water=100., T=371., P=101325., phase='g', units='kmol/hr') + dH = b.H - a.H + assert dH > 0 + table = problem_table([a], [b], [False], 5.) + assert_allclose(table.point_H, [[-dH]]) + assert table.interval_H.shape == (1, 0) + assert_allclose([table.hot_util_load, table.cold_util_load, table.pinch_T], + [dH, 0., 371.]) + table = problem_table([b], [a], [True], 5.) + assert_allclose([table.hot_util_load, table.cold_util_load, table.pinch_T], + [0., dH, 372. - 5.]) # outlet T 372 K, shifted by T_min_app + +def test_synthetic_network_reaches_MER(): + units = synthetic_units() + HXN = bst.HeatExchangerNetwork('HXN', T_min_app=5.) + sys = bst.System.from_units('sys_synthetic', units=[*units, HXN]) + sys.simulate() + hus = heat_utilities(units) + table = problem_table(*pinch_streams(hus), 5.) + actual_heat = sum(hu.unit_duty for hx in HXN.new_HX_utils + for hu in hx.heat_utilities if hu.unit_duty > 0) + # the greedy heuristic reaches the (corrected) MER on this case; it can + # never legitimately beat it + assert_allclose(actual_heat, table.hot_util_load, rtol=1e-2) + assert actual_heat >= table.hot_util_load * (1 - 1e-3) + if __name__ == '__main__': test_cache_network_matches_fresh_synthesis() test_cache_network_perturbed_feed() test_cache_network_duplicate_IDs() test_energy_balance_error_contributions_ignored_none() + test_problem_table_energy_consistency_doctest_system() + for T_min_app in (5., 10., 20.): + test_problem_table_energy_consistency_synthetic(T_min_app) + test_problem_table_two_streams_closed_form() + test_problem_table_non_monotone_stream_is_point_load() + test_synthetic_network_reaches_MER() From 08e873c57a710be85540a613ae56837bb7b3b188 Mon Sep 17 00:00:00 2001 From: sarangbhagwat Date: Sat, 22 Aug 2026 16:17:02 -0700 Subject: [PATCH 2/5] fix HXN problem table: evaluate streams at real T, clip to [H_in, H_out], point loads for isothermal streams temperature_interval_pinch_analysis built the heat cascade from interval enthalpies that were (1) evaluated by flashing hot streams at the *shifted* temperature T - T_min_app, so a stream entering within T_min_app of its dew point was already liquid and lost its latent heat, (2) never clipped to the stream's actual [H_in, H_out], so non-equilibrium column outlets inflated interval duties, and (3) assigned with the wrong sign to streams whose outlet temperature moves against their duty. Targets were impossible (hot target above the un-integrated heating load on the class doctest; hot - cold != net demand on sugarcane). The cascade is now built by problem_table(): each monotone stream is walked once down the shifted grid with a warm-started copy flashed at the real temperature, exact at its own end points and clipped in between, so its contributions telescope exactly to its duty; isothermal and non-monotone streams are point loads at their outlet temperature. Hence hot_util_load - cold_util_load == sum(unit_duty) and both targets are non-negative by construction. Signature and return value of temperature_interval_pinch_analysis are unchanged. Validated by tests/test_hxn.py: per-stream and net energy identities on the doctest system and a 4-stream case at dT = 5/10/20 K, two-stream closed forms, point-load handling, and the synthetic network reaching MER. Co-Authored-By: Claude Fable 5 --- biosteam/facilities/hxn/hxn_synthesis.py | 188 ++++++++++++++++------- 1 file changed, 133 insertions(+), 55 deletions(-) diff --git a/biosteam/facilities/hxn/hxn_synthesis.py b/biosteam/facilities/hxn/hxn_synthesis.py index 5d5ff658..f2d4bb56 100644 --- a/biosteam/facilities/hxn/hxn_synthesis.py +++ b/biosteam/facilities/hxn/hxn_synthesis.py @@ -10,6 +10,7 @@ @author: sarangbhagwat """ +from collections import namedtuple import numpy as np import biosteam as bst from warnings import warn @@ -124,8 +125,124 @@ 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): + """ + 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); 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. + """ + 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) + stream = stream_in.copy() + for k, T in enumerate(Ts): + if T == T_hi: + Hs[k] = H_top + elif T == T_lo: + Hs[k] = H_bottom + else: + 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!r} at {T_real:.2f} K " + f"({error!r}); interpolating enthalpy linearly in " + "temperature for the problem table", RuntimeWarning) + 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), `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 + minimum 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. + """ + 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]) + 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]) + residual = np.cumsum( + point_H.sum(axis=0) + np.concatenate([[0.], interval_H.sum(axis=0)]) + ) + k_pinch = int(np.argmin(residual)) + scale = np.abs(H_out - H_in).sum() + if -residual[k_pinch] <= 1e-9 * scale: # threshold problem: no hot utility + hot_util_load = 0. + k_pinch = 0 + else: + hot_util_load = -residual[k_pinch] + cold_util_load = residual[-1] + hot_util_load + 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 +267,38 @@ 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 the stream is split between the + # hot-side and cold-side designs. Non-monotone streams (T_out against + # the duty) are not split: their pinch is the inlet T, so load_duties + # puts the whole duty on the hot side, as before this fix. 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 From bcc56277ea5db0377c2cbada390538102d1219b2 Mon Sep 17 00:00:00 2001 From: sarangbhagwat Date: Sat, 22 Aug 2026 16:25:50 -0700 Subject: [PATCH 3/5] harden HXN problem table: endpoint enthalpies by index, traceable VLE warning, clean restart after failed flash Three review fixes to the problem_table helpers added on this branch: 1. _stream_H_at_boundaries compared grid temperatures to T_hi/T_lo with float equality to detect a stream's own end points; this only worked because Ts was built from the same floats. Assign Hs[0]/Hs[-1] by position instead (the caller always passes T_hi first, T_lo last, and monotone streams have T_hi > T_lo strictly so the slice has >= 2 entries, asserted), and loop only over the interior boundaries. 2. The VLE-failure warning printed {stream!r} for an anonymous inlet copy, giving no way to trace which stream failed. Pass the inlet stream's own ID (set by the wrapper to s___Util_) into the helper and name it in the warning. 3. On a failed flash, the warm-started copy was left in whatever state the failed VLE call put it in and reused for the next boundary. Re-copy stream_in in the except branch so the next boundary restarts clean. Also rewrote the comment above the pinch_T_arr loop in temperature_interval_pinch_analysis, which claimed non-monotone streams get pinch = T_in "as before this fix" -- not true in general (e.g. T_out < T_in <= pinch_cold_stream_T for a cold stream gave T_out/pinch under the old rule). The comment now states the actual rule (streams already on one side of the pinch, including non-monotone ones, are not split; load_duties assigns their whole duty to one side) without claiming equivalence with prior behavior. Covering tests unchanged (tests/test_hxn.py + hxn doctests): 12 passed in 13.36s. Co-Authored-By: Claude Fable 5 --- biosteam/facilities/hxn/hxn_synthesis.py | 76 +++++++++++++++--------- 1 file changed, 48 insertions(+), 28 deletions(-) diff --git a/biosteam/facilities/hxn/hxn_synthesis.py b/biosteam/facilities/hxn/hxn_synthesis.py index f2d4bb56..6df32c10 100644 --- a/biosteam/facilities/hxn/hxn_synthesis.py +++ b/biosteam/facilities/hxn/hxn_synthesis.py @@ -131,38 +131,50 @@ def get_sorted_life_cycle(self): '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): +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); 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. + 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, T in enumerate(Ts): - if T == T_hi: - Hs[k] = H_top - elif T == T_lo: - Hs[k] = H_bottom - else: - 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!r} at {T_real:.2f} K " - f"({error!r}); interpolating enthalpy linearly in " - "temperature for the problem table", RuntimeWarning) - H = H_lo + (H_hi - H_lo) * (T - T_lo) / (T_hi - T_lo) - Hs[k] = min(max(H, H_lo), H_hi) + 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): @@ -222,7 +234,8 @@ def problem_table(streams_inlet, streams_quenched, is_hot, T_min_app): 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]) + 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]) @@ -279,10 +292,17 @@ def temperature_interval_pinch_analysis(hus, 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 - # Per-stream pinch temperature: where the stream is split between the - # hot-side and cold-side designs. Non-monotone streams (T_out against - # the duty) are not split: their pinch is the inlet T, so load_duties - # puts the whole duty on the hot side, as before this fix. + # 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 or T_in_arr[i] > T_out_arr[i]: From 446e7e1dbb3167601d982cd45a638818d7919919 Mon Sep 17 00:00:00 2001 From: sarangbhagwat Date: Sat, 22 Aug 2026 18:05:04 -0700 Subject: [PATCH 4/5] fix HXN problem table cascade: check heat arriving at a boundary before its point loads The cascade residual[k] = sum(point_H[:, :k+1]) + sum(interval_H[:, :k]) is the heat *leaving* boundary Ts[k], after that boundary's point loads. Testing only this for non-negativity lets a hot point load at Ts[k] (an isothermal condenser, a phase-changing HXutility) mask a deficit in the cold interval (Ts[k-1], Ts[k]) directly above it - but a source at Ts[k] cannot serve a sink above Ts[k]. Example: a condensing stream at 400 K (shifted 395 K) against a cold stream heated 392 -> 398 K returned a hot utility target of zero; the correct target is the cold 395-398 K segment. The target is now the minimum over both the arriving flow (residual - point_total) and the leaving flow at every boundary, which is the standard problem-table treatment of point loads. Also: cold_util_load is clamped at zero in the threshold branch (residual[-1] could be negative by a rounding-level amount, and the tests assert non-negativity); problem_table/ProblemTable are exported in __all__ with an Examples doctest (two-stream threshold case). Note on synthesis behaviour (unchanged by this commit, introduced with the wrapper rewrite): non-monotone streams (outlet temperature moving against the duty) get pinch_T = T_in in pinch_T_arr, so load_duties places their whole duty on the hot side. The old code gave a meaningless negative dH2 for such streams when T_in was on the wrong side of the pinch. The table counts them as point loads at T_out; reconciling that with the synthesis heuristic is a separate item. Validation: regression test test_problem_table_point_load_cannot_heat_above_itself; tests/test_hxn.py + hxn doctests 14 passed (class doctest unchanged); canonical suite 71 failed / 476 passed / 62 skipped with the same pre-existing failure set as the baseline; tests/test_biorefineries.py 1 failed (test_cornstover, baseline) / 4 passed. Co-Authored-By: Claude Fable 5 --- biosteam/facilities/hxn/hxn_synthesis.py | 61 +++++++++++++++++++----- tests/test_hxn.py | 30 ++++++++++++ 2 files changed, 80 insertions(+), 11 deletions(-) diff --git a/biosteam/facilities/hxn/hxn_synthesis.py b/biosteam/facilities/hxn/hxn_synthesis.py index 6df32c10..b2770f33 100644 --- a/biosteam/facilities/hxn/hxn_synthesis.py +++ b/biosteam/facilities/hxn/hxn_synthesis.py @@ -15,7 +15,8 @@ import biosteam as bst from warnings import warn -__all__ = ('StreamLifeCycle', 'synthesize_network') +__all__ = ('StreamLifeCycle', 'ProblemTable', 'problem_table', + 'synthesize_network') class LifeStage: @@ -197,7 +198,8 @@ def problem_table(streams_inlet, streams_quenched, is_hot, T_min_app): 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), `hot_util_load`, + 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 @@ -210,10 +212,39 @@ def problem_table(streams_inlet, streams_quenched, is_hot, T_min_app): 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 - minimum 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. + 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, 3) + 1445547.086 + >>> table.pinch_T + 395.0 """ N = len(streams_inlet) is_hot = np.asarray(is_hot, dtype=bool) @@ -240,17 +271,25 @@ def problem_table(streams_inlet, streams_quenched, is_hot, T_min_app): 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_H.sum(axis=0) + np.concatenate([[0.], interval_H.sum(axis=0)]) + point_total + np.concatenate([[0.], interval_H.sum(axis=0)]) ) - k_pinch = int(np.argmin(residual)) + # 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 -residual[k_pinch] <= 1e-9 * scale: # threshold problem: no hot utility + if -flow[k_pinch] <= 1e-9 * scale: # threshold problem: no hot utility hot_util_load = 0. k_pinch = 0 else: - hot_util_load = -residual[k_pinch] - cold_util_load = residual[-1] + hot_util_load + hot_util_load = -flow[k_pinch] + # clamp: in the threshold branch residual[-1] may be negative by a + # rounding-level amount, and a negative cold utility is meaningless + cold_util_load = max(0., residual[-1] + hot_util_load) return ProblemTable(Ts, interval_H, point_H, residual, hot_util_load, cold_util_load, Ts[k_pinch]) diff --git a/tests/test_hxn.py b/tests/test_hxn.py index eb1cc28e..5520b6fe 100644 --- a/tests/test_hxn.py +++ b/tests/test_hxn.py @@ -210,6 +210,32 @@ def test_problem_table_non_monotone_stream_is_point_load(): assert_allclose([table.hot_util_load, table.cold_util_load, table.pinch_T], [0., dH, 372. - 5.]) # outlet T 372 K, shifted by T_min_app +def test_problem_table_point_load_cannot_heat_above_itself(): + # A source at shifted temperature T cannot serve sinks above T: an + # isothermal condensing hot stream at 400 K (shifted to 395 K) must not + # cover the 395-398 K segment of a cold stream that runs 392 -> 398 K. + bst.settings.set_thermo(['Water', 'Ethanol'], cache=True) + hot_in = bst.Stream('hv', Water=100., T=400., P=101325., phase='g', + units='kmol/hr') + hot_out = bst.Stream('hl', Water=100., T=400., P=101325., phase='l', + units='kmol/hr') + P_hot = hot_in.H - hot_out.H + assert P_hot > 0 + cold_in = bst.Stream('cl', Water=1000., T=392., P=5e5, phase='l', + units='kmol/hr') + cold_out = cold_in.copy('cl_out'); cold_out.vle(T=398., P=5e5) + s = cold_in.copy(); s.vle(T=395., P=5e5) + H_392, H_395, H_398 = cold_in.H, s.H, cold_out.H + assert P_hot > H_398 - H_392 # more than enough heat overall + table = problem_table([hot_in, cold_in], [hot_out, cold_out], + [True, False], 5.) + assert_allclose(table.Ts, [398., 395., 392.]) + # heat arriving at 395 K, before the point load there, is short by the + # 395-398 K segment of the cold stream + assert_allclose(table.hot_util_load, H_398 - H_395, rtol=1e-9) + assert table.pinch_T == 395. + assert_allclose(table.cold_util_load, P_hot - (H_395 - H_392), rtol=1e-9) + def test_synthetic_network_reaches_MER(): units = synthetic_units() HXN = bst.HeatExchangerNetwork('HXN', T_min_app=5.) @@ -222,6 +248,9 @@ def test_synthetic_network_reaches_MER(): # the greedy heuristic reaches the (corrected) MER on this case; it can # never legitimately beat it assert_allclose(actual_heat, table.hot_util_load, rtol=1e-2) + # the lower bound is exact here only because every synthetic inlet is an + # equilibrium state (so the clipped table is exact) and the synthesizer + # respects T_min_app; with non-equilibrium inlets the table is conservative assert actual_heat >= table.hot_util_load * (1 - 1e-3) if __name__ == '__main__': @@ -234,4 +263,5 @@ def test_synthetic_network_reaches_MER(): test_problem_table_energy_consistency_synthetic(T_min_app) test_problem_table_two_streams_closed_form() test_problem_table_non_monotone_stream_is_point_load() + test_problem_table_point_load_cannot_heat_above_itself() test_synthetic_network_reaches_MER() From bbea7059803161da2b505d017d64efbc6d077d96 Mon Sep 17 00:00:00 2001 From: sarangbhagwat Date: Sat, 22 Aug 2026 18:12:34 -0700 Subject: [PATCH 5/5] keep HXN problem-table energy identity exact when rounding makes cold utility negative The threshold branch could leave cold_util_load = residual[-1] negative by up to 1e-9 * scale; clamping it to zero alone would break hot_util_load - cold_util_load == sum(unit_duty) in relative terms (the identity the tests assert with rtol=1e-9). Absorb the rounding into hot_util_load instead, so both loads are non-negative and the identity is exact. Also round the problem_table doctest's cold utility to 10 kJ/hr so it does not demand ten significant digits from a VLE enthalpy across Python versions. tests/test_hxn.py + hxn doctests: 14 passed. Co-Authored-By: Claude Fable 5 --- biosteam/facilities/hxn/hxn_synthesis.py | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/biosteam/facilities/hxn/hxn_synthesis.py b/biosteam/facilities/hxn/hxn_synthesis.py index b2770f33..216485d8 100644 --- a/biosteam/facilities/hxn/hxn_synthesis.py +++ b/biosteam/facilities/hxn/hxn_synthesis.py @@ -241,8 +241,8 @@ def problem_table(streams_inlet, streams_quenched, is_hot, T_min_app): array([395., 390., 300., 295.]) >>> round(table.hot_util_load, 3) 0.0 - >>> round(table.cold_util_load, 3) - 1445547.086 + >>> round(table.cold_util_load, -1) + 1445550.0 >>> table.pinch_T 395.0 """ @@ -287,9 +287,13 @@ def problem_table(streams_inlet, streams_quenched, is_hot, T_min_app): k_pinch = 0 else: hot_util_load = -flow[k_pinch] - # clamp: in the threshold branch residual[-1] may be negative by a - # rounding-level amount, and a negative cold utility is meaningless - cold_util_load = max(0., residual[-1] + hot_util_load) + 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])