diff --git a/biosteam/facilities/hxn/_heat_exchanger_network.py b/biosteam/facilities/hxn/_heat_exchanger_network.py index 5edeaece..8e732a45 100644 --- a/biosteam/facilities/hxn/_heat_exchanger_network.py +++ b/biosteam/facilities/hxn/_heat_exchanger_network.py @@ -11,7 +11,7 @@ """ import biosteam as bst import numpy as np -from .hxn_synthesis import synthesize_network, StreamLifeCycle +from .hxn_synthesis import synthesize_network, StreamLifeCycle, plot_pinch_diagram from warnings import warn __all__ = ('HeatExchangerNetwork',) @@ -199,6 +199,8 @@ def _cost(self): self.force_ideal_thermo, self.avoid_recycle, self.sort_hus_by_T) new_HXs = HXs_hot_side + HXs_cold_side + self.new_HXs_hot_side = HXs_hot_side + self.new_HXs_cold_side = HXs_cold_side self.cold_indices = cold_indices self.original_heat_exchangers = hxs self.new_HXs = new_HXs @@ -370,6 +372,22 @@ def _get_stream_life_cycles(self): self.stream_life_cycles = stream_life_cycles return stream_life_cycles + def plot_pinch_diagram(self, file=None, **kwargs): + """ + Draw the pinch diagram of the synthesized network; see + :func:`~biosteam.facilities.hxn.hxn_synthesis.plot_pinch_diagram` + for the keyword arguments. Returns the matplotlib figure and axes. + """ + if not hasattr(self, 'new_HXs_hot_side'): + raise RuntimeError('simulate the heat exchanger network before ' + 'plotting its pinch diagram') + return plot_pinch_diagram( + self.stream_life_cycles, self.inlet_Ts, self.outlet_Ts, + self.new_HXs_hot_side, self.new_HXs_cold_side, + Qmin=self.Qmin, original_hxs=self.original_heat_exchangers, + file=file, **kwargs, + ) + def get_original_hxs_associated_with_streams(self): # pragma: no cover original_units = self.system.units original_heat_utils = self.original_heat_utils diff --git a/biosteam/facilities/hxn/hxn_synthesis.py b/biosteam/facilities/hxn/hxn_synthesis.py index 216485d8..1041c205 100644 --- a/biosteam/facilities/hxn/hxn_synthesis.py +++ b/biosteam/facilities/hxn/hxn_synthesis.py @@ -11,12 +11,13 @@ @author: sarangbhagwat """ from collections import namedtuple +import heapq import numpy as np import biosteam as bst from warnings import warn __all__ = ('StreamLifeCycle', 'ProblemTable', 'problem_table', - 'synthesize_network') + 'synthesize_network', 'plot_pinch_diagram') class LifeStage: @@ -671,3 +672,303 @@ def get_T_transient_hot_side(index): T_out_arr, pinch_T_arr, C_flow_vector, hx_utils_rearranged, streams_inlet, stream_HXs_dict,\ hot_indices, cold_indices + + +# Pinch diagram + +def _order_exchanger_columns(hxs, stream_life_cycles): + """ + Order heat exchangers left to right so that every stream meets its + exchangers in flow direction (cold streams flow left to right, hot streams + right to left). The per-stream stage orders define a precedence graph; + a topological sort (Kahn's algorithm, ties broken by the given order) + yields a consistent layout. Contradictory constraints, which would need a + stream to flow backwards, fall back to the given order. + """ + hxs = list(hxs) + position = {hx: i for i, hx in enumerate(hxs)} + successors = {hx: [] for hx in hxs} + N_predecessors = {hx: 0 for hx in hxs} + for life_cycle in stream_life_cycles: + stages = [i.unit for i in life_cycle.life_cycle if i.unit in position] + if not life_cycle.cold: stages.reverse() + for a, b in zip(stages, stages[1:]): + if b not in successors[a]: + successors[a].append(b) + N_predecessors[b] += 1 + ready = [position[hx] for hx in hxs if not N_predecessors[hx]] + heapq.heapify(ready) + ordered = [] + while ready: + hx = hxs[heapq.heappop(ready)] + ordered.append(hx) + for other in successors[hx]: + N_predecessors[other] -= 1 + if not N_predecessors[other]: heapq.heappush(ready, position[other]) + return ordered if len(ordered) == len(hxs) else hxs + +def _format_H(H): + mantissa, exponent = f'{H:.2e}'.split('e') + return f'{mantissa}E{int(exponent)}' + +def _auxiliary_name(unit): + """ + Return the (dotted) name of an auxiliary unit within its owner, e.g. + 'condenser' or 'evaporators[0].heat_exchanger', or None if the unit is + not auxiliary. + """ + owner = unit.owner + if owner is unit: return None + def search(parent, prefix): + for name, aux in parent.get_auxiliary_units_with_names(): + if aux is unit: return prefix + name + if hasattr(aux, 'get_auxiliary_units_with_names'): + found = search(aux, prefix + name + '.') + if found: return found + return search(owner, '') or unit.ID.lstrip('.') + +def _stream_label(hx, show_units, show_auxiliary_units, show_stream_IDs): + """ + Label of a stream from its original heat exchanger `hx`: + ' - ()', with each part + optional. + """ + parts = [] + if show_units: parts.append(hx.owner.ID) + if show_auxiliary_units: + auxname = _auxiliary_name(hx) + if auxname: parts.append(auxname) + label = ' - '.join(parts) + if show_stream_IDs: + ID = hx.ins[0].ID + if ID: label = f'{label} ({ID})' if label else ID + return label + +def plot_pinch_diagram(stream_life_cycles, inlet_Ts, outlet_Ts, + hot_side_HXs, cold_side_HXs, Qmin=1e-3, + original_hxs=None, show_units=True, + show_auxiliary_units=True, show_stream_IDs=True, + show_legend=True, ax=None, file=None, dpi=300): + """ + Draw a pinch diagram of a synthesized heat exchanger network: cold + streams (blue, flowing left to right) above hot streams (red, flowing + right to left), one vertical connector per process heat exchanger with + its duty, a dashed pinch line separating the cold-side from the + hot-side exchangers, and circles marking the utility exchangers that + bring each stream to its outlet temperature. + + Parameters + ---------- + stream_life_cycles : list[StreamLifeCycle] + One per stream, as built by HeatExchangerNetwork. + inlet_Ts, outlet_Ts : array-like + Stream inlet and outlet temperatures [K], indexed like the life cycles. + hot_side_HXs, cold_side_HXs : list[HXprocess] + Process exchangers above and below the pinch. + Qmin : float, optional + Utility exchangers with a duty at or below this [kJ/hr] are not marked. + original_hxs : list[Unit], optional + The original heat exchanger of each stream (indexed like the life + cycles). Required for the stream labels below. + show_units : bool, optional + Label each stream with the unit operation that owns its original heat + exchanger (the main unit for auxiliary exchangers). + show_auxiliary_units : bool, optional + Label each stream with the name of its original heat exchanger within + the main unit (e.g. 'condenser'), if it is an auxiliary unit. + show_stream_IDs : bool, optional + Label each stream with the ID of the original heat exchanger's inlet. + show_legend : bool, optional + Add a legend of the symbols below the diagram. + ax : matplotlib.axes.Axes, optional + Axes to draw on; a new figure is created if not given. + file : str, optional + If given, the figure is saved to this path. + dpi : int, optional + Resolution used when saving. + + Returns + ------- + fig : matplotlib.figure.Figure + ax : matplotlib.axes.Axes + + Notes + ----- + Temperatures are shown in degC and heat flows in kJ/hr at the inlet and + outlet of each stream. Exchanger columns on each side of the pinch are + ordered so that each stream meets them in flow direction whenever the + network allows it. Stream labels read ' - ()' + next to the stream index at the inlet. + + Examples + -------- + >>> import biosteam as bst + >>> bst.settings.set_thermo(['Water', 'Methanol', 'Glycerol']) + >>> feed1 = bst.Stream('feed1', flow=(8000, 100, 25)) + >>> feed2 = bst.Stream('feed2', flow=(10000, 1000, 10)) + >>> D1 = bst.ShortcutColumn('D1', ins=feed1, + ... outs=('distillate', 'bottoms_product'), + ... LHK=('Methanol', 'Water'), + ... y_top=0.99, x_bot=0.01, k=2, + ... is_divided=True) + >>> D1_H1 = bst.HXutility('D1_H1', ins = D1.outs[1], T = 300) + >>> D1_H2 = bst.HXutility('D1_H2', ins = D1.outs[0], T = 300) + >>> F1 = bst.Flash('F1', ins=feed2, + ... outs=('vapor', 'liquid'), V = 0.9, P = 101325) + >>> HXN = bst.HeatExchangerNetwork('HXN', T_min_app = 5.) + >>> sys = bst.System.from_units('sys', units=[D1, D1_H1, D1_H2, F1, HXN]) + >>> sys.simulate() + >>> fig, ax = HXN.plot_pinch_diagram() + >>> connectors = [i for i in ax.findobj() if (i.get_gid() or '').startswith('HX:')] + >>> len(connectors) == len(HXN.new_HXs) + True + >>> import matplotlib.pyplot as plt + >>> plt.close(fig) + + """ + import matplotlib.pyplot as plt + # Artists carry stable gids ('HX:', 'Util:', 'Label:') + # so the drawing can be checked structurally in tests. + show_labels = show_units or show_auxiliary_units or show_stream_IDs + if show_labels and original_hxs is None: + raise ValueError('original_hxs is required to label streams with ' + 'units, auxiliary units, or stream IDs') + cold_color, hot_color = '#2e6db4', '#d62728' + cold_bg, hot_bg = '#e6f0fa', '#fbe9e7' + process_hxs = set(hot_side_HXs) | set(cold_side_HXs) + # Stream index and stage of each side of every process exchanger, by identity + hx_streams = {hx: {} for hx in process_hxs} + for index, life_cycle in enumerate(stream_life_cycles): + for stage in life_cycle.life_cycle: + if stage.unit in hx_streams: + hx_streams[stage.unit][life_cycle.cold] = (index, stage) + cold_side_HXs = _order_exchanger_columns(cold_side_HXs, stream_life_cycles) + hot_side_HXs = _order_exchanger_columns(hot_side_HXs, stream_life_cycles) + columns = cold_side_HXs + hot_side_HXs + N_cs = len(cold_side_HXs) + N_columns = len(columns) + # x layout: 0 stream ends | 1 cold utilities | 2..N_cs+1 cold side | + # pinch | N_cs+2..N+1 hot side | N+2 hot utilities | N+3 stream ends + x_start, x_cold_util = 0., 1. + x_columns = {hx: 2. + i for i, hx in enumerate(columns)} + x_pinch = N_cs + 1.5 + x_hot_util = N_columns + 2. + x_end = N_columns + 3. + # y layout: cold streams on top, hot streams below, duty labels in between + cold_streams = [i for i, lc in enumerate(stream_life_cycles) if lc.cold] + hot_streams = [i for i, lc in enumerate(stream_life_cycles) if not lc.cold] + N_hot = len(hot_streams) + N_cold = len(cold_streams) + gap = 2.5 + y = {} + for k, i in enumerate(hot_streams): y[i] = N_hot - k + for k, i in enumerate(cold_streams): y[i] = N_hot + gap + N_cold - k + y_label = N_hot + (gap + 1.) / 2. + y_top = N_hot + gap + N_cold + 1. + y_bottom = 0. + if ax is None: + fig, ax = plt.subplots( + figsize=(max(6., 0.75 * (N_columns + 4) + 3.), 0.4 * y_top + 1.) + ) + else: + fig = ax.figure + # Background and pinch line + x_min, x_max = x_start - 1.8, x_end + 1.8 + ax.axvspan(x_min, x_pinch, color=cold_bg, lw=0, zorder=0) + ax.axvspan(x_pinch, x_max, color=hot_bg, lw=0, zorder=0) + ax.axvline(x_pinch, color='k', ls='--', lw=1, zorder=1) + ax.text(x_min + 0.2, y_bottom + 0.1, 'Cold side', color=cold_color, + weight='bold', ha='left', va='bottom') + ax.text(x_max - 0.2, y_bottom + 0.1, 'Hot side', color=hot_color, + weight='bold', ha='right', va='bottom') + # Column headers + header_kwargs = dict(ha='center', va='bottom', weight='bold', fontsize=8) + for x_T, x_H in ((x_start - 1.3, x_start - 0.6), (x_end + 0.6, x_end + 1.3)): + ax.text(x_T, y_top, 'T\n[°C]', **header_kwargs) + ax.text(x_H, y_top, 'H\n[kJ·h$^{-1}$]', **header_kwargs) + ax.text(x_start - 0.3, y_label, 'ΔH\n[kJ·h$^{-1}$]', + ha='right', va='center', weight='bold', fontsize=8) + # Streams + value_kwargs = dict(ha='center', va='center', fontsize=8) + for index, life_cycle in enumerate(stream_life_cycles): + cold = life_cycle.cold + color = cold_color if cold else hot_color + yi = y[index] + stages = life_cycle.life_cycle # never empty: each stream has a utility stage + H_in = stages[0].H_in + H_out = stages[-1].H_out + T_in = inlet_Ts[index] - 273.15 + T_out = outlet_Ts[index] - 273.15 + # T is the outer column on the left and the inner column on the right + T_left, H_left, T_right, H_right = ( + (T_in, H_in, T_out, H_out) if cold else (T_out, H_out, T_in, H_in) + ) + x_in, x_out, sign = (x_start, x_end, 1) if cold else (x_end, x_start, -1) + ax.annotate('', xy=(x_out, yi), xytext=(x_in, yi), + arrowprops=dict(arrowstyle='-|>', color=color, lw=1.2, + shrinkA=0, shrinkB=0), zorder=2) + ax.text(x_start - 1.3, yi, f'{T_left:.1f}', color=color, **value_kwargs) + ax.text(x_start - 0.6, yi, _format_H(H_left), color=color, **value_kwargs) + ax.text(x_end + 0.6, yi, f'{T_right:.1f}', color=color, **value_kwargs) + ax.text(x_end + 1.3, yi, _format_H(H_right), color=color, **value_kwargs) + # Index and label share a baseline above the stream, clear of the + # exchanger circles + y_text = yi + 0.25 + ax.text(x_in + sign * 0.3, y_text, str(index), color=color, + ha='center', va='baseline', weight='bold', fontsize=9) + if show_labels: + label = _stream_label(original_hxs[index], show_units, + show_auxiliary_units, show_stream_IDs) + # the smaller label reads as centered with the index when its + # baseline is slightly higher + ax.text(x_in + sign * 0.6, y_text + 0.08, label, color=color, + ha='left' if cold else 'right', va='baseline', fontsize=7, + zorder=6, gid=f'Label:{index}', + bbox=dict(boxstyle='square,pad=0.15', fc='w', ec='none')) + # Utility exchangers: a cold stream ends in a hot utility (red), a + # hot stream in a cold utility (blue) + x_util, util_color = (x_hot_util, hot_color) if cold else (x_cold_util, cold_color) + for stage in stages: + unit = stage.unit + if unit in process_hxs: continue + if abs(stage.H_out - stage.H_in) <= Qmin: continue + ax.plot([x_util], [yi], 'o', mfc='w', mec=util_color, mew=1.2, + ms=6, zorder=4, gid='Util:' + unit.ID) + # Process exchangers + for hx in columns: + streams = hx_streams[hx] + if len(streams) != 2: + warn(f'{hx.ID} is not in exactly one hot and one cold stream ' + 'life cycle; it is not drawn', RuntimeWarning) + continue + (i_cold, stage_cold), (i_hot, stage_hot) = streams[True], streams[False] + x = x_columns[hx] + Q = abs(stage_hot.H_in - stage_hot.H_out) + ax.plot([x, x], [y[i_hot], y[i_cold]], '-o', color='k', mfc='w', + mew=1.2, ms=6, lw=1.2, zorder=3, gid='HX:' + hx.ID) + ax.text(x, y_label, _format_H(Q), rotation=90, ha='center', + va='center', fontsize=8, zorder=5, + bbox=dict(boxstyle='square,pad=0.25', fc='w', ec='k', lw=0.8)) + if show_legend: + from matplotlib.lines import Line2D + handles = [ + Line2D([], [], color=cold_color, lw=1.2, marker='>', markevery=[-1], + ms=5, label='Cold stream'), + Line2D([], [], color=hot_color, lw=1.2, marker='<', markevery=[0], + ms=5, label='Hot stream'), + Line2D([], [], color='k', lw=1.2, marker='o', mfc='w', mew=1.2, + ms=6, label='Process heat exchange'), + Line2D([], [], ls='', marker='o', mfc='w', mec=hot_color, mew=1.2, + ms=6, label='Hot utility'), + Line2D([], [], ls='', marker='o', mfc='w', mec=cold_color, mew=1.2, + ms=6, label='Cold utility'), + Line2D([], [], color='k', ls='--', lw=1, label='Pinch'), + ] + ax.legend(handles=handles, loc='upper center', bbox_to_anchor=(0.5, 0.), + ncol=3, fontsize=7, frameon=False, handlelength=2.5, + columnspacing=1.5) + ax.set_xlim(x_min, x_max) + ax.set_ylim(y_bottom, y_top + 1.2) + ax.set_axis_off() + if file: fig.savefig(file, dpi=dpi, bbox_inches='tight') + return fig, ax diff --git a/tests/test_hxn.py b/tests/test_hxn.py index 5520b6fe..214372d5 100644 --- a/tests/test_hxn.py +++ b/tests/test_hxn.py @@ -253,6 +253,136 @@ def test_synthetic_network_reaches_MER(): # respects T_min_app; with non-equilibrium inlets the table is conservative assert actual_heat >= table.hot_util_load * (1 - 1e-3) +# --- pinch diagram ----------------------------------------------------------- + +class _FakeStage: + def __init__(self, unit): self.unit = unit + +class _FakeLifeCycle: + def __init__(self, units, cold=True): + self.cold = cold + self.life_cycle = [_FakeStage(u) for u in units] + +def test_pinch_diagram_column_order_follows_stream_direction(): + from biosteam.facilities.hxn.hxn_synthesis import _order_exchanger_columns + h1, h2, h3 = 'h1', 'h2', 'h3' + # stream A visits h2 then h1; stream B visits h1 then h3 -> h2, h1, h3 + cycles = [_FakeLifeCycle([h2, h1]), _FakeLifeCycle([h1, h3])] + assert _order_exchanger_columns([h1, h2, h3], cycles) == [h2, h1, h3] + # exchangers not in the requested subset are ignored, order is stable + assert _order_exchanger_columns([h3, h1], cycles) == [h1, h3] + # a hot stream flows right to left, so its stage order is reversed: + # hot stream visits h1 then h3 -> h3 left of h1 + cycles = [_FakeLifeCycle([h2, h1]), _FakeLifeCycle([h1, h3], cold=False)] + assert _order_exchanger_columns([h1, h2, h3], cycles) == [h2, h3, h1] + # contradictory constraints (a cycle) fall back to the given order + cycles = [_FakeLifeCycle([h1, h2]), _FakeLifeCycle([h2, h1])] + assert _order_exchanger_columns([h2, h1], cycles) == [h2, h1] + +def _gid_artists(ax, prefix): + return [a for a in ax.findobj() if (a.get_gid() or '').startswith(prefix)] + +def test_pinch_diagram_doctest_system(): + import matplotlib + matplotlib.use('Agg') + import matplotlib.pyplot as plt + sys, HXN, feed = build_system() + sys.simulate() + assert HXN.new_HXs_hot_side + HXN.new_HXs_cold_side == HXN.new_HXs + fig, ax = HXN.plot_pinch_diagram() + try: + # one connector per process exchanger + connectors = _gid_artists(ax, 'HX:') + assert {a.get_gid() for a in connectors} == {'HX:' + hx.ID for hx in HXN.new_HXs} + # one utility marker per utility exchanger with a duty above Qmin + utils = _gid_artists(ax, 'Util:') + expected = {'Util:' + hx.ID for hx in HXN.new_HX_utils + if abs(hx.outs[0].H - hx.ins[0].H) > HXN.Qmin} + assert {a.get_gid() for a in utils} == expected + # one row per stream, with inlet temperatures in degC + texts = {t.get_text() for t in ax.texts} + for T in HXN.inlet_Ts: assert f'{T - 273.15:.1f}' in texts + for T in HXN.outlet_Ts: assert f'{T - 273.15:.1f}' in texts + for i in range(len(HXN.inlet_Ts)): assert str(i) in texts + finally: + plt.close(fig) + +def test_pinch_diagram_stream_labels(): + import matplotlib + matplotlib.use('Agg') + import matplotlib.pyplot as plt + from biosteam.facilities.hxn.hxn_synthesis import _auxiliary_name, _stream_label + sys, HXN, feed = build_system() + sys.simulate() + D1 = bst.main_flowsheet.unit.D0 + D1_H1 = bst.main_flowsheet.unit.D0_H1 + assert _auxiliary_name(D1.condenser) == 'condenser' + assert _auxiliary_name(D1_H1) is None + # label composition + assert _stream_label(D1.condenser, True, True, False) == 'D0 - condenser' + assert _stream_label(D1.condenser, True, False, False) == 'D0' + assert _stream_label(D1.condenser, False, True, False) == 'condenser' + assert _stream_label(D1_H1, True, True, True) == 'D0_H1 (' + D1_H1.ins[0].ID + ')' + assert _stream_label(D1_H1, False, False, True) == D1_H1.ins[0].ID + assert _stream_label(D1_H1, False, False, False) == '' + # an unnamed inlet adds nothing + assert _stream_label(D1.reboiler, True, True, True) == 'D0 - reboiler' + # every stream gets a label on the figure; toggles remove them + fig, ax = HXN.plot_pinch_diagram() + try: + labels = {a.get_gid(): a.get_text() for a in _gid_artists(ax, 'Label:')} + hxs = HXN.original_heat_exchangers + assert labels == { + f'Label:{i}': _stream_label(hx, True, True, True) + for i, hx in enumerate(hxs) + } + assert any(text.startswith('D0 - condenser') for text in labels.values()) + assert any(text == 'F1 - heat_exchanger (feed_flash)' for text in labels.values()) + finally: + plt.close(fig) + fig, ax = HXN.plot_pinch_diagram(show_units=False, show_auxiliary_units=False, + show_stream_IDs=False) + try: + assert not _gid_artists(ax, 'Label:') + finally: + plt.close(fig) + +def test_pinch_diagram_requires_simulation(): + sys, HXN, feed = build_system() + with pytest.raises(RuntimeError, match='simulate'): + HXN.plot_pinch_diagram() + +def test_pinch_diagram_legend(): + import matplotlib + matplotlib.use('Agg') + import matplotlib.pyplot as plt + sys, HXN, feed = build_system() + sys.simulate() + fig, ax = HXN.plot_pinch_diagram() + try: + legend = ax.get_legend() + assert legend is not None + labels = [t.get_text() for t in legend.get_texts()] + assert labels == ['Cold stream', 'Hot stream', 'Process heat exchange', + 'Hot utility', 'Cold utility', 'Pinch'] + # utility markers are colored by utility type, consistent with the legend + handles = dict(zip(labels, legend.legend_handles)) + hot_util_color = handles['Hot utility'].get_markeredgecolor() + cold_util_color = handles['Cold utility'].get_markeredgecolor() + assert hot_util_color != cold_util_color + heaters = {hx.ID for hx in HXN.new_HX_utils if hx.outs[0].H > hx.ins[0].H} + for artist in _gid_artists(ax, 'Util:'): + heater = artist.get_gid()[len('Util:'):] in heaters + expected = hot_util_color if heater else cold_util_color + assert artist.get_markeredgecolor() == expected, artist.get_gid() + finally: + plt.close(fig) + fig, ax = HXN.plot_pinch_diagram(show_legend=False) + try: + assert ax.get_legend() is None + finally: + plt.close(fig) + if __name__ == '__main__': test_cache_network_matches_fresh_synthesis() test_cache_network_perturbed_feed() @@ -265,3 +395,8 @@ def test_synthetic_network_reaches_MER(): test_problem_table_non_monotone_stream_is_point_load() test_problem_table_point_load_cannot_heat_above_itself() test_synthetic_network_reaches_MER() + test_pinch_diagram_column_order_follows_stream_direction() + test_pinch_diagram_doctest_system() + test_pinch_diagram_stream_labels() + test_pinch_diagram_legend() + test_pinch_diagram_requires_simulation()