# -*- coding: utf-8 -*-
# hensmith: Heat Exchanger Network Synthesis, Modeling, Integration,
# Thermodynamics, and Heuristics
# Copyright (C) 2020-, Sarang Bhagwat <sarangbhagwat.developer@gmail.com>
#
# This module is under the UIUC open-source license. See
# github.com/BioSTEAMDevelopmentGroup/hensmith/blob/master/LICENSE.txt
# for license details.
"""
Pinch analysis and heat exchanger network synthesis: the problem table
(`problem_table`), the synthesis of a network at minimum energy requirement
(`synthesize_network`, on the planner of `hensmith._planner`), unsplit by
default and with stream splits where needed if `stream_splitting` (on
`hensmith._splitting`; realized splits are `StreamSplit`), stream life
cycles (`StreamLifeCycle`) and pinch diagrams (`plot_pinch_diagram`).
"""
from collections import namedtuple
from collections.abc import Iterable, Mapping, Sequence
import heapq
import math
import re
import numpy as np
import biosteam as bst
from warnings import warn
from ._curves import (StreamCurve, stream_curves, _end_state, _T_EQ, _T_SIDE,
_copy, _point_load_inlet)
from ._planner import (plan_network, _approach_shifts,
_approach_contributions, _contribution_value)
from ._splitting import Split, _ISO_TOL, _same_network
__all__ = ('StreamLifeCycle', 'ProblemTable', 'problem_table',
'synthesize_network', 'plot_pinch_diagram')
#: IDs of the synthesized exchangers: process exchangers above the pinch
#: ``HX_<cold>_<hot>_hs`` and below it ``HX_<hot>_<cold>_cs`` (the first
#: number is the stream at port 0), with ``_<n>`` for the n-th exchanger of a
#: repeated pair; utility exchangers ``Util_<index>_hs|cs``.
_PROCESS_ID = re.compile(r'^HX_(\d+)_(\d+)_(hs|cs)(?:_\d+)?$')
_UTILITY_ID = re.compile(r'^Util_(\d+)_(hs|cs)$')
def _stream_ports(unit):
"""Stream index at each inlet port of a synthesized exchanger, parsed
from its ID (``(a, b)`` for ``HX_<a>_<b>_...``, ``(a,)`` for
``Util_<a>_...``), or None if the ID is not one of the synthesizer's."""
ID = unit.ID
match = _PROCESS_ID.match(ID)
if match: return int(match.group(1)), int(match.group(2))
match = _UTILITY_ID.match(ID)
if match: return (int(match.group(1)),)
return None
#: An inlet port, ``unit.ins[index]`` (see `StreamLifeCycle.entry`).
_Port = namedtuple('_Port', ('unit', 'index'))
def _flow_order(trunk, splits, branches):
"""
A life cycle with splits in flow order: ``(None, stage)`` for every
stage of `trunk` (the stages on the whole flow, in flow order, the
utility last) and ``(split, branches[k])`` for split k of `splits` (in
flow order; `branches[k]` holds the stages of each of its branches),
placed after ``split.position`` trunk process stages and before the
utility.
"""
flow, k = [], 0
for t, stage in enumerate(trunk):
utility = isinstance(stage.unit, bst.HXutility)
while k < len(splits) and (utility or splits[k].position <= t):
flow.append((splits[k], branches[k]))
k += 1
flow.append((None, stage))
flow.extend(zip(splits[k:], branches[k:]))
return flow
def _split_flow(stages, splits, flow_order, index):
"""
Tag the `stages` of stream `index` on the branches of its `splits` (in
flow order) with their branch and fraction, sort the others (the
trunk) by `flow_order` with the utility last, and return
`_flow_order`'s flow of them.
"""
tags = {id(hx): (k, b, f) for k, split in enumerate(splits)
for b, (f, hxs) in enumerate(zip(split.fractions, split.branches))
for hx in hxs}
tagged, trunk = {}, []
for stage in stages:
if id(stage.unit) in tags: tagged[id(stage.unit)] = stage
else: trunk.append(stage)
branches = []
for split in splits:
branches.append([])
for hxs in split.branches:
branch = []
for hx in hxs:
stage = tagged.get(id(hx))
if stage is None:
raise ValueError(f'{hx.ID} of split {split!r} does '
f'not carry stream {index}')
k, b, stage.fraction = tags[id(hx)]
stage.branch = (k, b)
branch.append(stage)
branches[-1].append(branch)
trunk.sort(key=flow_order)
utilities = [s for s in trunk if isinstance(s.unit, bst.HXutility)]
trunk = [s for s in trunk if not isinstance(s.unit, bst.HXutility)]
return _flow_order(trunk + utilities, splits, branches)
[docs]
class LifeStage:
"""
One stage of a stream's passage through the synthesized network: the
heat exchanger it passes through and which of that exchanger's
inlet/outlet pairs carries the stream.
Parameters
----------
unit : HXprocess or HXutility
Heat exchanger of the synthesized network.
index : int
Position of the stream in `unit.ins` / `unit.outs` (0 or 1 for an
`HXprocess`; always 0 for an `HXutility`).
branch : tuple[int, int], optional
``(k, b)`` for a stage on branch b of the stream's split k (see
`StreamLifeCycle.splits`); None (the default) for a stage on the
whole flow.
fraction : float, optional
Fraction of the stream's flow through this stage: that of its
branch, 1 (the default) on the whole flow.
Attributes
----------
s_in : Stream
`unit.ins[index]`, the stream entering this stage.
s_out : Stream
`unit.outs[index]`, the stream leaving this stage.
H_in : float
Enthalpy of `s_in` [kJ/hr], read from the stream when accessed
(on a branch, the branch's: `fraction` of the whole flow's).
H_out : float
Enthalpy of `s_out` [kJ/hr], read from the stream when accessed.
"""
def __init__(self, unit, index, branch=None, fraction=1.):
self.unit = unit
self.index = index
self.branch = branch
self.fraction = fraction
@property
def s_in(self): return self.unit.ins[self.index]
@property
def s_out(self): return self.unit.outs[self.index]
@property
def H_in(self): return self.s_in.H
@property
def H_out(self): return self.s_out.H
def _branch_info(self):
branch = self.branch
if branch is None: return ''
return f", branch {branch}, fraction {self.fraction:.4g}"
def _info(self, N_tabs=1):
tabs = N_tabs*'\t'
return (f"{type(self).__name__}: {self.unit.ID}{self._branch_info()}\n"
+ tabs + f"H_in = {self.H_in:.3g} kJ/hr\n"
+ tabs + f"H_out = {self.H_out:.3g} kJ/hr")
def __repr__(self):
return (f"<{type(self).__name__}: {repr(self.unit)}{self._branch_info()}, H_in = {round(self.H_in, 4):.3g} kJ/hr, H_out = {round(self.H_out, 4):.3g} kJ/hr>")
def show(self):
print(self._info())
_ipython_display_ = show
[docs]
class StreamLifeCycle:
"""
The ordered sequence of heat exchangers one process stream passes through
in a synthesized heat exchanger network, from its inlet state to its
final utility exchanger.
Streams are numbered by their position in the rearranged utility list
of `synthesize_network` (cold streams first, then hot streams); the
network's stream copies and exchanger IDs embed that index
(``s_<index>__<exchanger ID>`` for the inlet streams of the exchangers,
``HX_<hot>_<cold>_cs`` / ``HX_<cold>_<hot>_hs`` for process exchangers,
with a suffix ``_<n>`` for the n-th exchanger of a repeated pair,
``Util_<index>_cs`` / ``Util_<index>_hs`` for utility exchangers), which
is how the life cycle is recovered from the exchangers: the IDs are
parsed (the first number is the stream at port 0, the second the stream
at port 1), so the stream indices are matched exactly and never as
substrings of other indices or of rewired stream IDs.
Parameters
----------
index : int
Stream index in the synthesized network.
cold : bool
True for a heated (cold) stream, False for a cooled (hot) stream.
Attributes
----------
index : int
Stream index in the synthesized network.
cold : bool
Whether the stream is a heated (cold) stream.
name : str
``'s_<index>'``, the prefix of the stream's copies in the network.
life_cycle : list[LifeStage] or None
Stages in flow order, set by `get_life_cycle`; None until then.
Where the stream splits, the stages of every branch follow those
before the split, branch by branch (see `get_life_cycle`).
splits : list[StreamSplit]
The stream's splits in flow order, set by `get_life_cycle`; empty
for an unsplit stream. Split k's branch stages have ``branch ==
(k, b)``.
entry : tuple[Unit, int] or None
``(unit, index)``, a named tuple: the inlet port `unit.ins[index]`
where the whole stream enters the network, i.e. the first
splitter of its first split if it splits at its inlet, else its
first stage's. None until `get_life_cycle` runs, and for a life
cycle whose `life_cycle` is assigned directly.
Notes
-----
`HeatExchangerNetwork` builds one life cycle per stream after synthesis
and stores them in `HeatExchangerNetwork.stream_life_cycles`, aligned
with `HeatExchangerNetwork.original_heat_exchangers`;
`plot_pinch_diagram` draws them.
"""
# Defaults for an instance without them, e.g. one unpickled from before
# stream splitting: an unsplit stream with no entry port recorded.
splits = ()
entry = None
def __init__(self, index, cold):
self.index = index
self.name = 's_%s'%index
self.cold = cold
self.life_cycle = None
self.splits = []
self.entry = None
@property
def H_in(self):
"""Enthalpy of the stream at `entry` [kJ/hr], read from the stream
when accessed: the whole flow's inlet, equal to the first stage's
`H_in` unless the stream splits at its inlet. Without an `entry`
(a `life_cycle` assigned directly), the first stage's `H_in`."""
entry = self.entry
if entry is None: return self.life_cycle[0].H_in
unit, index = entry
return unit.ins[index].H
[docs]
def get_relevant_units(self, index, new_HXs, new_HX_utils):
"""
Return the process and utility exchangers (two lists) that carry
stream `index`: those whose parsed ID (``HX_<a>_<b>_<hs|cs>[_<n>]``
or ``Util_<a>_<hs|cs>``) names it; for any other ID, those whose ID
contains ``_<index>_``.
"""
def relevant(hx):
ports = _stream_ports(hx)
if ports is None: return '_%s_'%index in hx.ID
return index in ports
new_HXs_relevant = [hx for hx in new_HXs if relevant(hx)]
new_HX_utils_relevant = [hx for hx in new_HX_utils if relevant(hx)]
return new_HXs_relevant, new_HX_utils_relevant
[docs]
def get_life_cycle(self, new_HXs, new_HX_utils, splits=None):
"""
Build and return the list of `LifeStage` objects for this stream.
Parameters
----------
new_HXs : list[HXprocess]
Process exchangers of the synthesized network.
new_HX_utils : list[HXutility]
Utility exchangers of the synthesized network.
splits : list[StreamSplit], optional
Splits of the network (``info['splits']`` of
`synthesize_network`); those of other streams are ignored.
Without any split of this stream, the life cycle is built
exactly as for an unsplit network.
Returns
-------
list[LifeStage]
One stage per port that carries this stream: for an exchanger
ID of the synthesizer (see the class notes) the port its ID
assigns to `index` (0 or 1 for a process exchanger, 0 for a
utility exchanger); for any other ID, every port among 0 and 1
(0 for a utility) whose inlet ID contains ``'s_<index>_'``.
Sorted in flow direction: by inlet enthalpy, ascending for a
cold stream and descending for a hot one; ties (zero-duty
stages only) put the stream's first side of the pinch first
(cold-side stages for a cold stream, hot-side stages for a hot
one) and the utility last. Also stored as `life_cycle`.
With splits, only the trunk stages (those on the whole flow,
which pass the stream's enthalpies in order) are sorted so;
the exchangers of split k's branches (``StreamSplit.branches``,
by identity) follow the first ``StreamSplit.position`` trunk
process stages, branch by branch and each branch in flow
order, as stages with ``branch == (k, b)`` and the branch's
`fraction`; the utility stays last. Enthalpies of different
branches are never compared. The stream's splits are stored,
in flow order, as `splits`, and its entry port as `entry`.
Raises
------
ValueError
If an exchanger of one of the stream's splits does not carry
the stream.
"""
index = self.index
name = self.name
cold = self.cold
new_HXs_relevant, new_HX_utils_relevant =\
self.get_relevant_units(index, new_HXs, new_HX_utils)
life_cycle = []
for units, N_ports in ((new_HXs_relevant, 2), (new_HX_utils_relevant, 1)):
for unit in units:
ports = _stream_ports(unit)
if ports is None:
life_cycle.extend([LifeStage(unit, k) for k in range(N_ports)
if name + '_' in unit.ins[k].ID])
else:
life_cycle.extend([LifeStage(unit, k) for k, i in enumerate(ports)
if i == index])
sign = 1. if cold else -1.
first_side = '_cs' if cold else '_hs'
def flow_order(stage):
ID = stage.unit.ID
if isinstance(stage.unit, bst.HXutility): rank = 2
else: rank = 0 if first_side in ID else 1
return (sign * stage.H_in, rank)
splits = [split for split in splits or () if split.stream == index]
if splits:
# flow order: by position, then (consecutive splits) by the
# enthalpy where they split
splits.sort(key=lambda split: (split.position,
sign * split.H_split))
flow = _split_flow(life_cycle, splits, flow_order, index)
life_cycle = [stage for split, part in flow
for stage in ([part] if split is None
else [s for b in part for s in b])]
else:
life_cycle.sort(key=flow_order)
flow = [(None, stage) for stage in life_cycle[:1]]
self.life_cycle = life_cycle
self.splits = splits
if flow:
split, part = flow[0]
self.entry = (_Port(part.unit, part.index) if split is None
else _Port(split.splitters[0], 0))
else:
self.entry = None
return life_cycle
def _flow(self):
"""The life cycle in flow order as ``(None, stage)`` for every trunk
stage and ``(split, branches)`` for every split, `branches` being
the stages of each of its branches (see `_flow_order`)."""
splits = self.splits
branches = [[[] for hxs in split.branches] for split in splits]
trunk = []
for stage in self.life_cycle:
branch = stage.branch
if branch is None: trunk.append(stage)
else: branches[branch[0]][branch[1]].append(stage)
return _flow_order(trunk, splits, branches)
[docs]
def connections(self):
"""
Yield the stream's connections in flow order, as ``(up_unit,
up_port, down_unit, down_port)``: ``up_unit.outs[up_port]`` feeds
``down_unit.ins[down_port]``. Wiring every connection of every
life cycle joins the network.
For an unsplit stream, these are its consecutive stages. At a
split: the stage before it feeds the first splitter (port 0);
splitter c's second outlet feeds splitter c + 1; the outlet of
branch b (``StreamSplit.outlet``) feeds the branch's first stage,
or inlet b of the mixer for a branch without exchangers (a
bypass); each branch stage feeds the next, the last one mixer
inlet b; the mixer's outlet feeds the next stage (or the next
split's first splitter).
"""
up = None # (unit, port) feeding the next node
for split, part in self._flow():
if split is None:
if up is not None: yield (*up, part.unit, part.index)
up = (part.unit, part.index)
continue
splitters = split.splitters
if up is not None: yield (*up, splitters[0], 0)
for a, b in zip(splitters, splitters[1:]): yield (a, 1, b, 0)
for b, stages in enumerate(part):
end = split.outlet(b)
for stage in stages:
yield (*end, stage.unit, stage.index)
end = (stage.unit, stage.index)
yield (*end, split.mixer, b)
up = (split.mixer, 0)
[docs]
def stage_pairs(self):
"""
Yield the pairs ``(a, b)`` of stages where the stream flows from
`a` into `b` directly or through splitters and mixers only: for an
unsplit stream, its consecutive stages. The stage before a split
precedes the first stage of every branch, and the last stage of
every branch the stage after the split; a bypass carries the stages
before its split past it. Stages on sibling branches make no pair.
"""
frontier = [] # the stages whose outflow reaches the next node
for split, part in self._flow():
if split is None:
for stage in frontier: yield stage, part
frontier = [part]
continue
ends = []
for stages in part:
if stages:
for stage in frontier: yield stage, stages[0]
yield from zip(stages, stages[1:])
last = [stages[-1]]
else: # a bypass
last = frontier
ends.extend([stage for stage in last
if not any(stage is end for end in ends)])
frontier = ends
def __repr__(self):
life_cycle = self.life_cycle
cold = self.cold
if not self.life_cycle:
return 'Not initialized; run StreamLifeCycle.get_life_cycle or\
HX_Network.get_stream_life_cycles first.'
else:
index = self.index
name = 'Stream_%s'%index
strtype = 'cold' if cold else 'hot'
rep = ''
for LifeStage in life_cycle:
line = '\t\t' + repr(LifeStage) + '\n'
rep += line
rep = '<StreamLifeCycle: ' + name + ', ' + strtype + '\n\tlife_cycle = [\n' + rep[:-1] + '\n\t]'
for k, split in enumerate(self.splits):
fractions = ', '.join(f'{f:.4g}' for f in split.fractions)
rep += (f'\n\tsplit {k}: {len(split.fractions)} branches '
f'({fractions})')
return rep + '>'
[docs]
def show(self):
"""Print the life cycle, one stage per line."""
info = repr(self).replace('[', '').replace(']', '').replace('life_cycle =', 'life_cycle:')
print(info[1:-1])
_ipython_display_ = show
ProblemTable = namedtuple(
'ProblemTable',
['Ts', 'interval_H', 'point_H', 'residual',
'hot_util_load', 'cold_util_load', 'pinch_T']
)
ProblemTable.__doc__ = """
Result of `problem_table`: the temperature-interval heat cascade of a set of
process streams on the *shifted* temperature scale (hot streams shifted down
by the minimum approach temperature, cold streams unshifted; with approach
contributions every stream by its own shift, see `problem_table`).
Attributes
----------
Ts : numpy.ndarray
Shifted grid temperatures [K], descending: every shifted breakpoint of
every stream's temperature-enthalpy curve, i.e. its end temperatures,
the phase boundaries inside its range (a pure component's saturation
temperature, a mixture's bubble and dew points), the samples of its
two-phase glides and of its curved single-phase stretches
(temperature-dependent heat capacity), and the outlet temperature of
every point-load stream (see `point_H`). Temperatures closer than
1e-9 K are one grid point.
interval_H : numpy.ndarray
(N streams x n-1 intervals) heat contributed by each stream to each
open interval (Ts[k], Ts[k+1]) [kJ/hr]: positive for hot streams (heat
released), negative for cold streams (heat required); zero outside the
stream's own temperature range.
point_H : numpy.ndarray
(N x n) heat contributed *at* each grid temperature [kJ/hr], with the
same sign convention: the jump of the stream's curve there, i.e. a pure
component's latent heat at its (shifted) saturation temperature, the
enthalpy by which a non-equilibrium end state departs from equilibrium
at its own end temperature, and the whole duty of a point-load stream
at its shifted outlet temperature (isothermal streams and streams whose
outlet temperature moves against their duty).
residual : numpy.ndarray
(n,) heat cascaded *leaving* each grid temperature, after its point
loads, when no hot utility is supplied [kJ/hr]; negative where that
cascade is infeasible.
hot_util_load : float
Minimum hot utility target [kJ/hr] (zero for a threshold problem).
cold_util_load : float
Minimum cold utility target [kJ/hr].
pinch_T : float
Shifted grid temperature of the pinch [K] (``Ts[0]`` for a threshold
problem). Stream j's pinch temperature is ``pinch_T + s_j``, its shift:
``pinch_T + T_min_app`` (hot) and ``pinch_T`` (cold) by default,
``pinch_T + T_min_app/2 + c_j`` (hot) and ``pinch_T + T_min_app/2 -
c_j`` (cold) with approach contributions.
See Also
--------
problem_table : builds the table and documents the cascade.
"""
[docs]
def problem_table(streams_inlet, streams_quenched, is_hot, T_min_app,
curves=None, approach_contributions=None):
"""
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]: the required approach of every
pair of streams unless `approach_contributions` are given.
curves : list, optional
Prebuilt temperature-enthalpy curves of the same streams, in the
same order (e.g. shared with the network synthesis); built here if
not given.
approach_contributions : Sequence[float or None], optional
Approach temperature contribution of each stream [K], in the order
of `streams_inlet`: None or nan (unspecified: ``T_min_app / 2``), or
a finite number > 0 (`T_min_app` must then be finite and > 0 too).
A hot and a cold stream then need the sum of their contributions
between them (see Notes). Defaults to None (every pair needs
`T_min_app`).
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
-----
Each stream is described by a piecewise-linear temperature-enthalpy
curve built once from a handful of flashes: breakpoints at its end
temperatures and at every phase boundary inside its range (a pure
component's saturation temperature, a mixture's bubble and dew points);
a flat (isothermal) segment for a pure component's latent heat, between
its saturated-liquid and saturated-vapor enthalpies; samples of each
mixture glide (binaries traced along their bubble-point curve); and
interior breakpoints of each curved single-phase stretch
(temperature-dependent Cp), dense enough that linear interpolation is
within 0.002 K of the true curve. Single-phase stretches are evaluated
with their phases fixed (no flash), so a grid point exactly at a
saturation temperature is never ambiguous. Hot streams are shifted down
by `T_min_app`; cold streams are not (see below for approach
contributions). The grid is the union of all
shifted breakpoints, so every point at which any stream's curve bends
or jumps is a grid point, every stream is within 0.002 K of linear
between grid points, and the grid minimum of the cascade is the true
one to within ``0.002 K * sum(CP)``.
At a grid temperature a stream contributes the jump of its curve there
as a point load (`point_H`), and between two grid temperatures the heat
of its curve in that open interval (`interval_H`). Inside its own
temperature range a stream is taken at equilibrium with its enthalpy
clipped to its real range, so a non-equilibrium end state (e.g. a
superheated liquid from a non-rigorous HXutility) can never inflate the
duty; what it departs from equilibrium at its own end temperature is a
point load there. Every stream's contributions therefore telescope
exactly to ``sign * |H_out - H_in|``. Streams whose outlet temperature
does not move with their duty (isothermal, or 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
first location the pinch. With the per-stream identity above,
hot_util_load - cold_util_load equals the net heating demand.
With `approach_contributions` ``c_j`` (Kemp 2007), hot streams are
shifted down by ``T_min_app/2 + c_j`` and cold streams by
``T_min_app/2 - c_j``, so that a hot and a cold stream at the same
shifted temperature are ``c_h + c_c`` apart: every match must keep the
sum of its two streams' contributions. An unspecified contribution is
``T_min_app/2``; with none specified, the table is the default one, bit
for bit.
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. (The grid
between the ends holds the breakpoints that follow the curvature of
liquid water's enthalpy.)
>>> import biosteam as bst
>>> from hensmith.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[[0, -1]] # shifted grid ends (hot streams 5 K down)
array([395., 295.])
>>> round(table.hot_util_load, 3)
0.0
>>> round(table.cold_util_load, -1)
1445550.0
>>> table.pinch_T
395.0
If the cold stream contributes 10 K and the hot one the default 2.5 K,
the cold stream needs 12.5 K against the hot one, which a 400 K inlet
cannot give against a 390 K outlet, so hot utility appears. The hot
stream's pinch temperature is ``395 + 5 = 400`` K and the cold one's
``395 - 7.5 = 387.5`` K:
>>> table = problem_table([hot_in, cold_in], [hot_out, cold_out],
... [True, False], 5., approach_contributions=[None, 10.])
>>> table.Ts[[0, -1]]
array([397.5, 295. ])
>>> round(table.hot_util_load, 3)
171724.348
>>> table.pinch_T
395.0
"""
return _problem_table(streams_inlet, streams_quenched, is_hot, T_min_app,
curves, approach_contributions)[0]
def _problem_table(streams_inlet, streams_quenched, is_hot, T_min_app,
curves=None, approach_contributions=None):
"""
Return the `ProblemTable` of `problem_table` together with the stream
curves it was built from and its grid enthalpies, as ``(table, curves,
grid)``.
`grid` is a dict with
* 'Ts': the table's shifted grid temperatures (descending, n);
* 'shift': each stream's shift (`T_min_app` for hot streams, 0 for
cold ones by default; ``T_min_app/2 + c_j`` and ``T_min_app/2 - c_j``
with `approach_contributions`; N), so stream j's real temperature at
grid index k is ``Ts[k] + shift[j]``;
* 'contributions': each stream's approach temperature contribution
[K], unspecified ones (all by default) at ``T_min_app / 2`` (N);
* 'Hl', 'Hr': (N x n) each stream's enthalpy at every grid temperature,
left (low-enthalpy) and right (high-enthalpy) limit; they differ only
where the stream's curve has a flat (point load) at that grid
temperature, and are H_hi above and H_lo below the stream's range;
* 'k_hi', 'k_lo': (N,) grid indices of each stream's own T_hi and T_lo
breakpoints (``k_hi <= k_lo``; equal for a point-load stream), so its
range is selected by position, not by comparing shifted floats.
The table is exactly ``point_H = sign * (Hr - Hl)`` and ``interval_H =
sign * (Hl[:, :-1] - Hr[:, 1:])`` (sign +1 hot, -1 cold): piecewise-
linear curves through the knots (Ts[k] + shift[j], Hl[j, k]) and
(Ts[k] + shift[j], Hr[j, k]) reproduce the table's cascade exactly.
"""
N = len(streams_inlet)
is_hot = np.asarray(is_hot, dtype=bool)
sign = np.where(is_hot, 1., -1.)
shift = _approach_shifts(is_hot, T_min_app, approach_contributions)
if curves is None:
curves = stream_curves(streams_inlet, streams_quenched, is_hot)
elif len(curves) != N:
raise ValueError(f'{len(curves)} curves given for {N} streams')
H_in = np.array([c.H_in for c in curves])
H_out = np.array([c.H_out for c in curves])
shifted = [c.T - shift[j] for j, c in enumerate(curves)]
sizes = [x.size for x in shifted]
values, inverse = np.unique(np.concatenate(shifted), return_inverse=True)
# merge grid temperatures closer than 1e-9 K (ascending clusters)
cluster = np.concatenate([[0], np.cumsum(np.diff(values) > _T_EQ)])
n = int(cluster[-1]) + 1
# the highest member of each cluster represents it (selected explicitly:
# NumPy does not specify which value a repeated fancy index keeps)
top = values[np.flatnonzero(np.append(np.diff(cluster) > 0, True))]
Ts = top[::-1].copy()
position = (n - 1) - cluster[inverse] # descending grid index of each breakpoint
offsets = np.concatenate([[0], np.cumsum(sizes)])
Hl = np.empty((N, n))
Hr = np.empty((N, n))
k_hi = np.empty(N, dtype=int)
k_lo = np.empty(N, dtype=int)
for j, c in enumerate(curves):
index = position[offsets[j]:offsets[j + 1]]
Hl[j], Hr[j] = c.grid_limits(Ts, shift[j], index)
k_hi[j] = index[-1]
k_lo[j] = index[0]
point_H = sign[:, None] * (Hr - Hl)
interval_H = sign[:, None] * (Hl[:, :-1] - Hr[:, 1:])
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.
table = ProblemTable(Ts, interval_H, point_H, residual,
hot_util_load, cold_util_load, Ts[k_pinch])
grid = dict(Ts=Ts, shift=shift, Hl=Hl, Hr=Hr, k_hi=k_hi, k_lo=k_lo,
contributions=_approach_contributions(
N, T_min_app, approach_contributions))
return table, curves, grid
def _pinch_cut(table):
"""
Return which side of the pinch the point loads *at* the pinch
temperature belong to: 'below' if the zero-heat-flow cut is the flow
arriving at `pinch_T` (before its point loads), 'above' if it is the
flow leaving it (after them); always 'below' for a threshold problem
(no hot utility: everything lies below the pinch at ``Ts[0]``). Split
every stream with ``side = 'right' if cut == 'below' else 'left'``
(see `pinch_state`) to agree with the table: the heat above the split
is then exactly the hot utility target and the heat below it the cold
one. (`synthesize_network` does not cut streams at the pinch, even
with `stream_splitting`, whose splits are parallel branches: the
planner finds the same cut in its own cascade, ``plan.cut``, which
reproduces the table's.)
"""
if table.hot_util_load == 0.: return 'below'
k = int(np.flatnonzero(table.Ts == table.pinch_T)[0])
point_total = table.point_H.sum(axis=0)
arriving = table.residual - point_total
return 'below' if arriving[k] <= table.residual[k] else 'above'
def _parent(unit, owners=None):
"""The unit one step up the owner chain of `unit`: ``unit._owner`` or,
without one, the unit that holds `unit` by `owners` (see
`_held_owners`), or None for a top unit."""
owner = getattr(unit, '_owner', None)
if owner is None and owners:
held = owners.get(id(unit))
if held is not None and held[0] is unit: owner = held[1]
return owner
def _top(unit, owners=None):
"""The top owner of `unit` along `_parent` (`unit` itself without an
owner)."""
while (owner := _parent(unit, owners)) is not None: unit = owner
return unit
def _held_owners(units, held=None):
"""
The owners that the units at hand give the auxiliary units they hold
without setting their owner: a dict from the id of such an auxiliary to
``(auxiliary, holder)``, for each `biosteam.Unit` without an ``_owner``
among the `auxiliary_units` of a unit of `units` or of one on their
owner chains (the first holder wins; none that would close a loop).
`biosteam.System.from_units` sets the owners of the auxiliaries that
exist when it is called, but some units build theirs anew without one:
`biosteam.MultistageCompressor` rebuilds its compressors and
intercoolers whenever `pr` or `n_stages` changes, so the keys of its
ID cover them only by this relation (see `_parent`). `held`, such a
dict taken earlier (a snapshot, which keeps its objects alive), comes
first: its auxiliaries keep their holders even once a holder has
rebuilt them (no longer holds them), and the units at hand add the
others, along its owners. Then a unit without an ``_owner`` that none
holds as an auxiliary, but whose heat utility is among the
`heat_utilities` of units on those chains (an exchanger kept outside
the `auxiliary_unit_names`, see `_heat_utility_units`), is held by the
lowest of them (an owner lists the heat utilities of its auxiliaries
too), none that would close a loop. Nothing is mutated.
"""
owners = dict(held) if held else {}
chains = _owner_chains(units, owners)
def closes_loop(aux, unit):
above = unit
while above is not None and above is not aux:
above = _parent(above, owners)
return above is not None
for unit in chains:
for aux in getattr(unit, 'auxiliary_units', ()):
if (not isinstance(aux, bst.Unit)
or getattr(aux, '_owner', None) is not None
or id(aux) in owners):
continue
if not closes_loop(aux, unit): owners[id(aux)] = (aux, unit)
by_heat_utilities = {}
for unit in chains:
for hx in _heat_utility_units([unit]):
if getattr(hx, '_owner', None) is None and id(hx) not in owners:
by_heat_utilities.setdefault(id(hx), (hx, []))[1].append(unit)
for hx, holders in by_heat_utilities.values():
holders = [unit for unit in holders if not closes_loop(hx, unit)]
# the lowest holder: one on the chain of no other holder
lowest = [unit for unit in holders
if not any(above is unit for other in holders
if other is not unit
for above in _owner_chains([other], owners)[1:])]
if lowest: owners[id(hx)] = (hx, lowest[0])
return owners
def _unit_name(unit, owners=None):
"""``(owner.ID, unit.ID)``, with the top owner (``(ID, ID)`` for a unit
without an owner; owners along `_parent`): the unit's tuple key."""
return _top(unit, owners).ID, unit.ID
def _id_path(unit, owners=None):
"""The IDs of `unit` and of its owners in turn (along `_parent`),
``(unit.ID, owner.ID, ..., top.ID)``, which tell apart auxiliary units
of the same ID under different owners."""
path = []
while unit is not None:
path.append(unit.ID)
unit = _parent(unit, owners)
return tuple(path)
def _unit_message_name(unit, owners=None):
"""How messages name a unit: its tuple key (`_unit_name`), followed by
its owners below the top one when there are any (the tuple key of an
auxiliary of an auxiliary is shared by its namesakes)."""
path = _id_path(unit, owners)
label = repr(_unit_name(unit, owners))
if len(path) <= 2: return label
return f"{label} (of {' of '.join(map(repr, path[1:-1]))})"
def _names_unit(key, unit, live=(), owners=None):
"""
Whether the approach contribution key `key` (of a valid type, see
`_check_contribution_key`) names `unit`, with owners along `_parent`
(`owners`, see `_held_owners`). A unit key names itself; a unit key
that is not a unit at hand (its id not in `live`, the ids of the units
at hand: a stale auxiliary that its owner rebuilt) names the unit of
the same `_id_path` instead. A tuple key names every unit of that
`_unit_name` (every auxiliary of that ID under that top owner, at any
depth) and a str key the unit of that ID without an owner.
"""
if isinstance(key, bst.Unit):
if key is unit: return True
return (id(key) not in live
and _id_path(key, owners) == _id_path(unit, owners))
if isinstance(key, tuple): return key == _unit_name(unit, owners)
return key == unit.ID and _parent(unit, owners) is None
def _contribution_key_level(key, unit, live=(), owners=None):
"""The number of steps along the owner chain of `unit` (`unit`, its
owner, ..., along `_parent` with `owners`) to the first unit that `key`
names (see `_names_unit`, which takes `live`), or None."""
level = 0
while unit is not None:
if _names_unit(key, unit, live, owners): return level
unit = _parent(unit, owners)
level += 1
return None
def _check_contribution_key(key):
"""Raise a TypeError unless `key` is a unit, a str or a pair of str."""
if isinstance(key, bst.HeatUtility):
raise TypeError('a heat utility cannot be an approach_contributions '
'key (biosteam replaces the heat utilities of a unit '
'at every simulation): use its unit (hu.unit), the '
"unit's ID or (owner.ID, unit.ID), or a callable (or, "
'outside a HeatExchangerNetwork, a sequence aligned '
'with the heat utilities)')
if isinstance(key, (bst.Unit, str)): return
if (isinstance(key, tuple) and len(key) == 2
and all(isinstance(i, str) for i in key)): return
raise TypeError('approach_contributions keys must be units, unit IDs '
'(str) or (owner.ID, unit.ID) tuples, not '
f'{type(key).__name__!r} ({key!r})')
def _owner_chains(units, owners=None):
"""Every unit on the owner chains of `units` (along `_parent` with
`owners`), once each."""
found = {}
for unit in units:
while unit is not None and id(unit) not in found:
found[id(unit)] = unit
unit = _parent(unit, owners)
return list(found.values())
def _nested_auxiliaries(unit):
"""
The `biosteam.Unit` objects that `unit` holds as auxiliary units, at
any depth, in the order of `biosteam.Unit.nested_auxiliary_units`.
Unlike that property, this reads the auxiliaries of units only (an
object that is not one, such as the `biosteam.VacuumSystem` of a flash
drum under vacuum, has none, and the property fails on it one level
down), and it takes each unit once, never `unit` itself, and walks
into each unit at most once (so units that hold each other end it).
"""
nested = []
found = {id(unit)} # the units listed (and `unit`)
walked = set() # the units `walk` went into
def add(aux):
if isinstance(aux, bst.Unit) and id(aux) not in found:
found.add(id(aux))
nested.append(aux)
def walk(holder):
# as `nested_auxiliary_units`: each attribute of `holder` (a unit or
# a list of them) before the auxiliaries of its units, and the
# auxiliaries of those that are units each before its own
walked.add(id(holder))
for name in getattr(holder, 'auxiliary_unit_names', ()):
held = getattr(holder, name, None)
if held is None: continue
group = list(held) if isinstance(held, Iterable) else [held]
for aux in group: add(aux)
for aux in group:
if not isinstance(aux, bst.Unit): continue
for sub in getattr(aux, 'auxiliary_units', ()):
add(sub)
if isinstance(sub, bst.Unit) and id(sub) not in walked:
walk(sub)
walk(unit)
return nested
def _heat_utility_units(units):
"""
The `biosteam.Unit` objects that the heat utilities of `units` belong
to (``hu.unit``), each once, other than the units of `units`: the
exchangers that a unit holds through its `heat_utilities` only, outside
its `auxiliary_unit_names`, as a unit may keep an exchanger as a plain
attribute and take its heat utilities in its `_design`
(``self.heat_utilities = self.hx.heat_utilities``; biosteam rejects
only the heat utilities of auxiliaries added again). A
`HeatExchangerNetwork` takes those heat utilities as any other
(`biosteam.process_tools.heat_exchanger_utilities_from_units` reads
`unit.heat_utilities`), so such exchangers are among its original
heat exchangers. An owner lists the heat utilities of its auxiliaries
too, whose units are its nested auxiliaries.
"""
found = {id(unit) for unit in units}
held = []
for unit in units:
for hu in getattr(unit, 'heat_utilities', None) or ():
hx = getattr(hu, 'unit', None)
if isinstance(hx, bst.Unit) and id(hx) not in found:
found.add(id(hx))
held.append(hx)
return held
def _with_auxiliaries(units):
"""Each unit of `units` followed by its nested auxiliary units (those
that are `biosteam.Unit` objects, see `_nested_auxiliaries`) and by the
exchangers that it or they hold through their heat utilities only (see
`_heat_utility_units`): the units that keys may name, as
`HeatExchangerNetwork` takes those available to it, whether or not they
hold a heat utility, and among which its original exchangers are found
while no unit has rebuilt them (see
`hensmith.targeting._rebuilt_label`)."""
candidates = []
for unit in units:
group = [unit, *_nested_auxiliaries(unit)]
candidates.extend(group)
candidates.extend(_heat_utility_units(group))
return candidates
def _key_universe(units, held=None):
"""The units that keys may name, given the units at hand `units`:
``(universe, live, owners)``, the owner chains of `units` (see
`_owner_chains`) along the owners they give the auxiliaries they hold
(`owners`, see `_held_owners`, after those of the snapshot `held`), and
the ids of the universe (`live`)."""
owners = _held_owners(units, held)
universe = _owner_chains(units, owners)
return universe, {id(u) for u in universe}, owners
def _auxiliary_attributes(key, universe, owners=None):
"""
The attributes that the tuple key `key`, of the ID that
`biosteam.Unit.auxiliary` gives an auxiliary, may name: ``(owner.ID,
'name[k]')`` the k-th unit of a list under the attribute `name`, and
``(owner.ID, name)`` a single unit there, with `name` one of the
`auxiliary_unit_names` of a unit under the top owner ``owner.ID``
(along `_parent` with `owners`). Yields ``(unit, name, k, held)`` for
each such unit, with k None for the single form and `held` what the
unit holds there (None if missing).
"""
owner_ID, ID = key
match = re.fullmatch(r'(.+)\[(0|[1-9]\d*)\]', ID)
name, k = (match[1], int(match[2])) if match else (ID, None)
for unit in universe:
if (_top(unit, owners).ID == owner_ID
and name in getattr(unit, 'auxiliary_unit_names', ())):
yield unit, name, k, getattr(unit, name, None)
def _auxiliary_slots(key, universe, owners=None):
"""
Where the units of `universe` hold the auxiliary unit of the tuple key
`key` by the ID that `biosteam.Unit.auxiliary` gives an auxiliary (see
`_auxiliary_attributes`). Yields ``(unit, held)`` for each such unit
whose attribute takes that form (a list for ``'name[k]'``) or holds
nothing (None or missing), with `held` the `biosteam.Unit` there, or
None if there is none (nothing, an object that is not a unit, or a
list too short).
"""
for unit, name, k, held in _auxiliary_attributes(key, universe, owners):
if held is not None:
if isinstance(held, list) != (k is not None): continue
if k is not None: held = held[k] if k < len(held) else None
yield unit, held if isinstance(held, bst.Unit) else None
def _other_form_hints(key, universe, owners=None):
"""
Hints for the tuple key `key` of the ID that `biosteam.Unit.auxiliary`
gives an auxiliary in the other form than the attribute it names takes
(see `_auxiliary_slots`, which skips such attributes): for ``(owner.ID,
name)`` where a unit under that top owner holds a list there, the keys
of its units (``"E1.evaporators is a list: ('E1', 'evaporators[0]'),
..."``), and for ``(owner.ID, 'name[k]')`` where it holds one unit,
that unit's key (``"D1.condenser is ('D1', 'condenser')"``).
"""
hints = []
for unit, name, k, held in _auxiliary_attributes(key, universe, owners):
if held is None or isinstance(held, list) == (k is not None): continue
if k is not None:
if isinstance(held, bst.Unit):
hints.append(f'; {unit.ID}.{name} is '
f'{_unit_name(held, owners)!r}')
continue
keys = [repr(_unit_name(aux, owners)) for aux in held
if isinstance(aux, bst.Unit)]
if keys:
shown = ', '.join(keys[:3]) + (', ...' if len(keys) > 3 else '')
hints.append(f'; {unit.ID}.{name} is a list: {shown}')
return hints
def _stage_tickets(compressor=None):
"""The ticket names that a `biosteam.MultistageCompressor` gives its
stages in their IDs (``f'{K.ID}_{ticket}{n}'``): those of the
`biosteam.IsentropicCompressor` and `biosteam.HXutility` stages it
builds from `pr` and `n_stages`, and those of the stages that
`compressor` holds now (given as `compressors` and `hxs`)."""
tickets = {getattr(getattr(bst, 'IsentropicCompressor', None),
'ticket_name', None), bst.HXutility.ticket_name}
if compressor is not None:
for stages in (getattr(compressor, 'compressors', None),
getattr(compressor, 'hxs', None)):
tickets.update(getattr(stage, 'ticket_name', None)
for stage in stages or ())
return {ticket for ticket in tickets if isinstance(ticket, str) and ticket}
def _is_stage_ID(ID, compressor_ID, tickets):
"""Whether `ID` is ``f'{compressor_ID}_{ticket}{n}'`` for a ticket name
of `tickets` and a stage number n >= 1 (see `_stage_tickets`)."""
prefix = f'{compressor_ID}_'
return ID.startswith(prefix) and any(
re.fullmatch(re.escape(ticket) + r'[1-9]\d*', ID[len(prefix):])
for ticket in tickets)
def _names_unbuilt_stage(key, universe, owners=None):
"""
Whether the tuple key `key` gives the ID of a stage that a
`biosteam.MultistageCompressor` K of `universe` under the top owner
``key[0]`` (along `_parent` with `owners`) may build:
``f'{K.ID}_{ticket}{n}'`` for any stage number n >= 1 (see
`_stage_tickets`). K sets its `auxiliary_unit_names` to the IDs of the
stages it builds whenever `pr` or `n_stages` changes, so ``('K1',
'K1_H3')`` names an intercooler of K1 with 3 stages or more and no unit
with fewer.
"""
MultistageCompressor = getattr(bst, 'MultistageCompressor', None)
if MultistageCompressor is None: return False
owner_ID, ID = key
return any(isinstance(unit, MultistageCompressor)
and _is_stage_ID(ID, unit.ID, _stage_tickets(unit))
and _top(unit, owners).ID == owner_ID
for unit in universe)
def _names_unbuilt_auxiliary(key, universe, live, owners=None):
"""
Whether the key `key` (of a valid type, see `_check_contribution_key`),
which names no unit of `universe` (see `_names_unit`; `live` holds the
ids of `universe`, `owners` the owners of `_parent`), names an
auxiliary unit that a unit of `universe` may build but does not hold
now, as a `biosteam.HXutilities` builds one exchanger per agent it
needs: an old auxiliary object whose top owner is in `universe`
(rebuilt with fewer auxiliaries), or a tuple key of the ID that biosteam
gives an auxiliary (see `_auxiliary_slots`) where a unit of `universe`
holds no unit, or holds the unit of that ID: ``('X1',
'heat_exchangers[2]')`` while the list ``X1.heat_exchangers`` is
shorter, say. Or a tuple key of the ID that a
`biosteam.MultistageCompressor` of `universe` gives a stage, of any
stage number (see `_names_unbuilt_stage`): ``('K1', 'K1_H3')`` while
K1 has 2 stages. Such a key is valid but names nothing. A slot that
holds a unit of another ID (``T1.heat_exchanger`` of a
`biosteam.EnzymeTreatment` T1 is ``('T1', '')``) does not make its key
valid (see `_unknown_key_hints`).
"""
if isinstance(key, bst.Unit):
top = _top(key, owners)
return top is not key and id(top) in live
if not isinstance(key, tuple): return False
return (any(held is None or held.ID == key[1]
for _, held in _auxiliary_slots(key, universe, owners))
or _names_unbuilt_stage(key, universe, owners))
def _unknown_keys(keys, universe, owners=None):
"""The keys of `keys` (units, str or pairs of str) that name neither a
unit of `universe` (see `_names_unit`, with owners along `_parent`)
nor an auxiliary unit that one of them may build (see
`_names_unbuilt_auxiliary`)."""
live = {id(u) for u in universe}
return [key for key in keys
if not any(_names_unit(key, u, live, owners) for u in universe)
and not _names_unbuilt_auxiliary(key, universe, live, owners)]
def _unknown_key_hints(keys, universe, owners=None):
"""
Hints for the unknown keys `keys` (see `_unknown_keys`) that name a
unit of `universe` (owners along `_parent` with `owners`) by another
key: for a str key, each auxiliary unit of that ID (``"'condenser' is
('D1', 'condenser')"``), and for a tuple key of the ID that biosteam
gives an auxiliary (see `_auxiliary_slots`), each unit held there under
another ID (``"T1.heat_exchanger is ('T1', '')"``), or the keys of the
units held there when the key takes the other form (a list's name
alone, or an index into a single one, see `_other_form_hints`). A tuple
key ``(owner_ID, ID)`` of a unit of `universe` without an owner gets that
unit's key (``"('K1', 'K1_H1') is 'K1_H1' (...)"``), and a str key,
the stages without an owner that a `biosteam.MultistageCompressor` of
that ID gives its IDs (see `_is_stage_ID`), such as the intercoolers
of the heat utilities of a compressor K1 outside a HeatExchangerNetwork
(which, with the `key_units` of `hensmith.PinchTargets.from_heat_utilities`,
alone takes K1 as their owner, see `_held_owners`). Returns them
as ``'; <hint>'`` each, joined ('' if none).
"""
hints = {}
tickets = _stage_tickets()
for key in keys:
if isinstance(key, str):
stages = []
for unit in universe:
if unit.ID == key and _parent(unit, owners) is not None:
hints[f'; {key!r} is {_unit_name(unit, owners)!r}'] = None
elif (_parent(unit, owners) is None
and _is_stage_ID(unit.ID, key, tickets)):
stages.append(unit.ID)
if stages:
one = len(stages) == 1
keyed = ('its ID is its key' if one
else 'their IDs are their keys')
hints[f"; {', '.join(map(repr, stages))} "
f"{'has' if one else 'have'} no owner here, so {keyed}"
f' (a MultistageCompressor {key!r} owns the stages it '
'builds only among the units available to a '
'HeatExchangerNetwork, or the key_units of '
'PinchTargets.from_heat_utilities)'] = None
elif isinstance(key, tuple):
for unit, held in _auxiliary_slots(key, universe, owners):
if held is not None and held.ID != key[1]:
hints[f'; {unit.ID}.{key[1]} is '
f'{_unit_name(held, owners)!r}'] = None
for hint in _other_form_hints(key, universe, owners):
hints[hint] = None
if key[0] != key[1] and any(
unit.ID == key[1] and _parent(unit, owners) is None
for unit in universe):
hints[f'; {key!r} is {key[1]!r} (a unit without an owner '
'here)'] = None
return ''.join(hints)
def _mapping_contributions(hus, mapping, units, held=None):
"""
`_stream_contributions` of a mapping from keys naming units to
contributions: checks the types of keys and values, then that every key
names a unit of the universe (the owner chains of the units of `hus`
and of `units`, the units at hand of `_names_unit`, along the owners
they give the auxiliaries they hold without one, after those of the
snapshot `held`, see `_key_universe`) or an auxiliary unit that one of them may build but does not hold now
(see `_names_unbuilt_auxiliary`: such a key names nothing), then that
no two keys of the same level for a unit of the universe give it
different values (None and nan are equal); then each heat utility takes
the value of the key of the lowest level on its unit's owner chain
(None without one).
"""
items = []
for key, value in mapping.items():
_check_contribution_key(key)
items.append((key, _contribution_value(value, f'key {key!r}')))
universe, live, owners = _key_universe([hu.unit for hu in hus]
+ list(units or ()), held)
unknown = _unknown_keys([key for key, _ in items], universe, owners)
if unknown:
raise ValueError('approach_contributions keys that name no unit: '
f"{', '.join(map(repr, unknown))} (a key is a unit, "
'the ID of a unit without an owner, or (owner.ID, '
'unit.ID)'
f'{_unknown_key_hints(unknown, universe, owners)})')
_check_key_levels(items, universe, live, owners, 'approach_contributions')
return [_nearest_key_value(items, hu.unit, live, owners) for hu in hus]
def _check_key_levels(items, universe, live, owners, what):
"""Raise a ValueError if two keys of `items` (``(key, value)`` pairs
of unit keys) of the same level for a unit of `universe` (see
`_contribution_key_level`, which takes `live` and `owners`) give it
different values; `what` names the keys (``'h'``, say)."""
for unit in universe:
first = {} # level -> (key, value)
for key, value in items:
level = _contribution_key_level(key, unit, live, owners)
if level is None: continue
other, other_value = first.setdefault(level, (key, value))
if other_value != value:
raise ValueError(
f'{what} keys {other!r} and {key!r} give unit '
f'{_unit_message_name(unit, owners)} different values '
f'at the same level: {other_value!r} and {value!r}')
def _nearest_key_value(items, unit, live, owners, default=None):
"""The value of the key of `items` (``(key, value)`` pairs of unit
keys) of the lowest level on the owner chain of `unit` (see
`_contribution_key_level`, which takes `live` and `owners`; the first
of that level), or `default` if none names it."""
best = None # (level, value) of the nearest key
for key, value in items:
level = _contribution_key_level(key, unit, live, owners)
if level is not None and (best is None or level < best[0]):
best = (level, value)
return default if best is None else best[1]
def _stream_contributions(hus, approach_contributions, units=None,
held=None):
"""
Resolve `approach_contributions` to one approach temperature
contribution per heat utility of `hus`, in the given order (zero-duty
ones included): a list of float, or None where unspecified (nan becomes
None; not filled with ``T_min_app / 2``). Accepts the forms of
`synthesize_network`'s parameter: None (all unspecified), a mapping
whose keys name units (see `_mapping_contributions`), a sequence
aligned with `hus` (not a str or bytes), or a callable of a heat
utility (called once per heat utility, in order). `units` adds units
(with their owner chains) to those that keys may name and that are
checked for conflicting keys, whether or not they hold a heat utility
of `hus`, and `held` gives the holders of auxiliaries taken earlier
(see `_held_owners`).
"""
N = len(hus)
x = approach_contributions
if x is None: return [None] * N
if isinstance(x, Mapping):
return _mapping_contributions(hus, x, units, held)
if isinstance(x, (Sequence, np.ndarray)) and not isinstance(x, (str, bytes)):
values = list(x)
if len(values) != N:
raise ValueError(f'{len(values)} approach contributions given for '
f'{N} heat utilities')
return [_contribution_value(v, f'heat utility {k}')
for k, v in enumerate(values)]
if callable(x):
# labeled with the owners the keys of a mapping would take
owners = _held_owners([hu.unit for hu in hus] + list(units or ()),
held)
return [_contribution_value(
x(hu), f'unit {_unit_message_name(hu.unit, owners)}')
for hu in hus]
raise TypeError('approach_contributions must be a mapping, a sequence or '
f'a callable, not {type(x).__name__!r}')
def _explicit_contributions(hus, approach_contributions, T_min_app,
units=None, held=None):
"""
`_stream_contributions` of `hus` (with `units` and `held`), or None if
nothing is specified (`approach_contributions` None, or every value
None or nan: the default path). With a value, `T_min_app` must be
finite and > 0.
"""
if approach_contributions is None: return None
values = _stream_contributions(hus, approach_contributions, units, held)
if all(v is None for v in values): return None
_approach_contributions(len(values), T_min_app, values) # checks T_min_app
return values
def _arrange_heat_utilities(hus, sort_hus_by_T):
"""
The heat utilities of `hus` in stream order, ``(hus_rearranged,
N_heating)``: the heating ones (``hu.duty > 0``, cold streams), then the
cooling ones; zero-duty ones are dropped. With `sort_hus_by_T`, heating
utilities are sorted by descending and cooling ones by ascending inlet
temperature of their unit.
"""
hus_heating = [hu for hu in hus if hu.duty > 0]
hus_cooling = [hu for hu in hus if hu.duty < 0]
if sort_hus_by_T:
hus_heating.sort(key=lambda i: i.unit.ins[0].T, reverse=True)
hus_cooling.sort(key=lambda i: i.unit.ins[0].T)
return hus_heating + hus_cooling, len(hus_heating)
def _process_streams(hxs, force_ideal_thermo):
"""
The process streams of the exchangers `hxs`, ``(streams_inlet,
streams_quenched)``: unregistered copies (see `hensmith._curves._copy`)
of each inlet and of each outlet re-flashed at its own enthalpy and
pressure, with ideal thermodynamics if `force_ideal_thermo`.
"""
if force_ideal_thermo:
streams_inlet = [_copy(hx.ins[0], hx.ins[0].thermo.ideal()) for hx in hxs]
streams_quenched = [_copy(hx.outs[0], hx.outs[0].thermo.ideal()) for hx in hxs]
else:
streams_inlet = [_copy(hx.ins[0]) for hx in hxs]
streams_quenched = [_copy(hx.outs[0]) for hx in hxs]
for i in streams_quenched: i.vle(H=i.H, P=i.P)
return streams_inlet, streams_quenched
def _pinch_analysis(hus, T_min_app=10, force_ideal_thermo=False,
sort_hus_by_T=False, approach_contributions=None):
"""
The first step of `synthesize_network`: prepare the process streams
behind `hus` and run the problem table on them. Returns the 12 values
of `temperature_interval_pinch_analysis` (which wraps this function)
and then the table's ``table, curves, grid`` (see `_problem_table`),
on which the network is planned, so that the curves are built only
once. `approach_contributions` (see `_stream_contributions`) is resolved
first, before any stream is copied or flashed.
"""
values = _explicit_contributions(hus, approach_contributions, T_min_app)
hx_utils_rearranged, N_heating = _arrange_heat_utilities(hus, sort_hus_by_T)
hxs = [hu.unit for hu in hx_utils_rearranged]
# unregistered copies; the inlets are registered under their own IDs
# below
streams_inlet, streams_quenched = _process_streams(hxs, force_ideal_thermo)
for i in range(len(streams_inlet)):
stream = streams_inlet[i]
ID = 'Util_%s'%i
stream.ID = 's_%s__%s'%(i,ID)
cold_indices = list(range(N_heating))
hot_indices = list(range(N_heating, len(hxs)))
indices = cold_indices + hot_indices
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
if values is None:
contributions = None
else: # in stream order
value_of = dict(zip(map(id, hus), values))
contributions = [value_of[id(hu)] for hu in hx_utils_rearranged]
table, curves, grid = _problem_table(streams_inlet, streams_quenched,
is_hot, T_min_app, None,
contributions)
hot_util_load = table.hot_util_load
cold_util_load = table.cold_util_load
pinch_T = table.pinch_T # shifted scale
shift = grid['shift']
# Per-stream pinch temperature, for information only (returned as
# `HeatExchangerNetwork.pinch_Ts`; `load_duties` splits a stream there):
# the network is planned on the curves (see `synthesize_network`), not
# on these temperatures. Stream i's pinch is the shifted pinch plus its
# shift, T_p (pinch_T for a cold and pinch_T + T_min_app for a hot
# stream by default). A stream already entirely on one side of the
# process pinch (T_in past T_p) is not split; its pinch_T is its own
# T_in. So is a non-monotone stream (T_out on the wrong side of T_in for
# its duty, e.g. a cold stream whose VLE outlet ends up cooler than it
# entered: a point load at T_out), whose whole duty `load_duties` then
# puts on a single side (hot side for a cold stream, cold side for a hot
# one).
pinch_T_arr = []
for i in cold_indices:
T_p = pinch_T + shift[i]
if T_in_arr[i] > T_p or T_in_arr[i] > T_out_arr[i]:
pinch_T_arr.append(T_in_arr[i])
elif T_out_arr[i] < T_p:
pinch_T_arr.append(T_out_arr[i])
else:
pinch_T_arr.append(T_p)
for i in hot_indices:
T_p = pinch_T + shift[i]
if T_in_arr[i] < T_p or T_in_arr[i] < T_out_arr[i]:
pinch_T_arr.append(T_in_arr[i])
elif T_out_arr[i] > T_p:
pinch_T_arr.append(T_out_arr[i])
else:
pinch_T_arr.append(T_p)
pinch_T_arr = np.array(pinch_T_arr)
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, table, curves, grid
def temperature_interval_pinch_analysis(hus,
T_min_app=10,
force_ideal_thermo=False,
sort_hus_by_T=False,
approach_contributions=None):
"""
Prepare the process streams behind `hus` and run the problem table on
them, as `synthesize_network` does first; a standalone pinch analysis
(the network itself is planned on the table's stream curves, which this
function does not return).
Heating utilities (``hu.duty > 0``, cold streams) come first, then
cooling utilities; zero-duty utilities are dropped. Each stream is a
copy of its exchanger's inlet (renamed ``s_<index>__Util_<index>``) and
of its outlet re-flashed at its own enthalpy. `approach_contributions`
takes the forms of `synthesize_network`'s parameter.
Returns
-------
pinch_T_arr : numpy.ndarray
Per-stream pinch temperature (see `synthesize_network`).
hot_util_load, cold_util_load : float
MER targets of `problem_table` [kJ/hr].
T_in_arr, T_out_arr : numpy.ndarray
Inlet and quenched outlet temperatures [K].
hxs : list[Unit]
The original heat exchangers, in stream order.
hot_indices, cold_indices, indices : list[int]
Stream indices of the hot streams, the cold streams, and all
(cold first).
streams_inlet, hx_utils_rearranged, streams_quenched : list
Inlet copies, heat utilities and quenched outlet copies, in stream
order.
"""
return _pinch_analysis(hus, T_min_app, force_ideal_thermo,
sort_hus_by_T, approach_contributions)[:12]
def pinch_state(stream_in, stream_out, T_pinch, side=None, curve=None):
"""
Return a copy of the stream in the state it has when it crosses the
pinch, with enthalpy guaranteed to lie within [min(H_in, H_out),
max(H_in, H_out)].
Parameters
----------
stream_in, stream_out : Stream
The stream's real end states (the outlet quenched to equilibrium at
its own enthalpy).
T_pinch : float
The stream's (real, unshifted) pinch temperature [K].
side : str, optional
Branch of a flat (isothermal) segment of the stream at `T_pinch`
(latent heat or a non-equilibrium end jump): 'right' takes its
high-enthalpy end (the load at the pinch goes below the pinch),
'left' its low-enthalpy end (the load goes above). Pass ``side =
'right' if cut == 'below' else 'left'``, with ``cut`` the pinch cut
of the table (`_pinch_cut`), to split every stream consistently
with the targets. If not given, an end temperature returns that
end's state (below) and elsewhere hot streams take 'right' and cold
streams 'left'.
curve : StreamCurve, optional
Prebuilt curve of the stream (e.g. from `_problem_table`); built
here if needed and not given.
Notes
-----
The state comes from the stream's temperature-enthalpy curve (see
`problem_table`), so it is deterministic: a pinch at the stream's own
saturation temperature is resolved by `side`, never by whatever phase
split a previous flash left, and a pinch inside a glide gets the
table's enthalpy. Inside the stream's temperature range the enthalpy is
clipped to its real range, as in the table: a non-equilibrium inlet
(e.g. a superheated liquid from a non-rigorous HXutility) has less
enthalpy than the equilibrium fluid at the pinch, so the state returned
is instead the equilibrium state at the nearer end enthalpy. Without
`side`, `T_pinch` equal to an end temperature returns that end's state
(the equilibrium state at the end enthalpy, see `_end_state`), because
flashing a non-equilibrium inlet at its own temperature does not
reproduce `H_in`. Using the *equilibrium* state at the end enthalpy,
rather than the stream as given, keeps the synthesizer consistent with
the problem table: the heat is offered at the temperature the
equilibrium model says it is available, not at a fictitious one.
Either way the hot-side and cold-side loads split `|H_in - H_out|`
exactly, and the state never carries heat the real stream does not
have.
A standalone analysis helper (see also `load_duties`): the network
synthesis does not cut streams at a pinch temperature (its only
splits, with `stream_splitting`, are parallel branches), it plans on
the curves themselves (see `synthesize_network`).
"""
if side is None:
T_lo, T_hi = sorted((stream_in.T, stream_out.T))
if T_pinch == stream_in.T: return _end_state(stream_in, T_lo, T_hi)
if T_pinch == stream_out.T: return _end_state(stream_out, T_lo, T_hi)
if curve is None: curve = StreamCurve(stream_in, stream_out)
if side is None: side = 'right' if curve.is_hot else 'left'
return curve.state_at_T(T_pinch, side)
def load_duties(streams, streams_quenched, pinch_T_arr, T_out_arr, indices,
is_cold, Q_hot_side, Q_cold_side):
"""
Fill `Q_hot_side` and `Q_cold_side` with each stream's duty above and
below its pinch temperature, ``[kind, duty]`` with kind 'heat' (cold
streams) or 'cool' (hot streams) and duties below 0.01 kJ/hr set to 0,
from the stream's `pinch_state` at ``pinch_T_arr[index]`` (e.g. from
`temperature_interval_pinch_analysis`). A standalone analysis helper,
like `pinch_state`: `synthesize_network` does not use it.
"""
for index in indices:
H_in = streams[index].H
H_out = streams_quenched[index].H
H_pinch = pinch_state(streams[index], streams_quenched[index],
pinch_T_arr[index]).H
if not is_cold(index):
dH1 = H_in - H_pinch
dH2 = H_pinch - H_out
if abs(dH1)<0.01: dH1 = 0
if abs(dH2)<0.01: dH2 = 0
Q_hot_side[index] = ['cool', dH1]
Q_cold_side[index] = ['cool', dH2]
else:
dH1 = H_out - H_pinch
dH2 = H_pinch - H_in
if abs(dH1)<0.01: dH1 = 0
if abs(dH2)<0.01: dH2 = 0
Q_hot_side[index] = ['heat', dH1]
Q_cold_side[index] = ['heat', dH2]
# %% Network synthesis
#: Exact-state approach acceptance [K]. Synthesized process exchangers get
#: ``HXprocess(dT=required - _APPROACH_TOL)``, with `required` their pair's
#: required approach, as a guard only: the planner and the exactness check
#: enforce the required approach on the exact states.
_APPROACH_TOL = 1e-6
#: A simulated process-exchanger duty that differs from the planned one by
#: more than this times the two streams' total duties is a deviation.
_DUTY_TOL = 1e-6
#: Status 'mer' needs the utilities of the realized network (from the
#: simulated duties) within this times the total stream duty of the targets:
#: the accuracy of the enthalpy flashes that realize the plan (the plan
#: itself reaches the targets within `hensmith._planner._MER_TOL`).
_ACHIEVED_TOL = 1e-6
#: Rounds of exact-state verification and local knot refinement.
_MAX_REFINE = 3
#: Extra rounds with stream splitting after the last refine round, each run
#: only if a violating exchanger lies on a split side whose network
#: (`hensmith._splitting._same_network`) is newly excluded. The retry
#: re-plans that side without its networks excluded so far: another network
#: if the side has one left, else the same one on the refined knots (a side's
#: last candidate is never excluded), after which no network is new.
_MAX_SPLIT_RETRY = 2
#: A split's mixer outlet further than this from the planned temperature [K]
#: is a deviation: 133 times the offset of a branch-flow PH flash (7.5e-10
#: K) and 10 times below the approach guard `_APPROACH_TOL`. It catches
#: thermosteam's PH flash disagreeing with the stream's curve (risk R-2),
#: measured at 0.8-2.4 K inside random water/methanol glides.
_MIX_T_TOL = 1e-7
def _grid_knots(curves, grid):
"""
Knots ``(T, H - H_lo)`` of every stream on the problem-table grid, the
planner's model of the stream: at every grid temperature inside the
stream's own range, its left and right enthalpy limits there (two knots
where the curve has a flat). Every breakpoint of every curve is a grid
point and the table is linear between grid points, so the planner's
cascade on these knots is exactly the table's. Enthalpies are relative
to the stream's H_lo, so that large absolute enthalpies cost the planner
no precision.
"""
Ts, shift, Hl, Hr = grid['Ts'], grid['shift'], grid['Hl'], grid['Hr']
knots = []
for j, curve in enumerate(curves):
k = np.arange(grid['k_lo'][j], grid['k_hi'][j] - 1, -1) # ascending T
T = np.repeat(Ts[k] + shift[j], 2)
H = np.column_stack((Hl[j, k], Hr[j, k])).ravel() - curve.H_lo
keep = np.ones(T.size, dtype=bool)
keep[1::2] = Hr[j, k] != Hl[j, k]
knots.append((T[keep], H[keep]))
return knots
def _knot_T(knots, H, hot):
"""
Temperature of a knot curve at enthalpies `H`, linear between knots; a
vertical stretch (a clipped non-equilibrium end) counts at its lowest
temperature for a hot stream and at its highest for a cold one, as in
the planner.
"""
T, Hk = knots
rises = np.diff(Hk) > 0.
if hot: keep = np.concatenate(([True], rises))
else: keep = np.concatenate((rises, [True]))
return np.interp(H, Hk[keep], T[keep])
def _curve_tol_T(curve):
"""Largest distance [K] between a stream's linearized curve (and so its
grid knots) and its exact states."""
if not curve.monotone: return 0.
return max(curve.tol_T, curve.glide_error)
def _interval_min(f, a, fa, b, fb):
"""
Minimum of a smooth function `f` on [a, b], given its end values.
A minimum inside the interval shows at an end as a slope that points
into it (f falls from `a`, or rises into `b`); both one-sided slopes are
taken from a step of 1e-6 of the interval. Without such a slope the
minimum is an end (the smooth pieces of temperature-enthalpy curves
bend one way over a knot interval). Otherwise a scan of the interval
brackets the dip and golden-section search narrows it to 1e-7 of the
interval. Returns the smallest value found.
"""
eps = 1e-6 * (b - a)
fa_ = f(a + eps)
fb_ = f(b - eps)
if fa_ >= fa and fb_ >= fb: return min(fa, fb)
xs = [a, a + eps, *np.linspace(a, b, 9)[1:-1].tolist(), b - eps, b]
fs = [fa, fa_, *[f(x) for x in xs[2:-2]], fb_, fb]
k = int(np.argmin(fs))
if k in (0, len(xs) - 1): return fs[k]
lo, x, hi, fx = xs[k - 1], xs[k], xs[k + 1], fs[k]
g = 0.5 * (3. - 5.**0.5) # golden section
xtol = 1e-7 * (b - a)
while hi - lo > xtol:
u = x + g * (hi - x) if hi - x > x - lo else x - g * (x - lo)
fu = f(u)
if fu < fx:
if u > x: lo = x
else: hi = x
x, fx = u, fu
elif u > x: hi = u
else: lo = u
return fx
def _exchanger_approach(curves, knots, h, c, H_hot_in, H_cold_in, Q,
T_min_app, fh=1., fc=1.):
"""
Exact-state check of one counter-current exchanger of duty `Q` in which
stream `h` enters hot at `H_hot_in` and stream `c` enters cold at
`H_cold_in` (enthalpies relative to each stream's H_lo); `T_min_app` is
the pair's required approach.
With stream splitting, `fh` and `fc` are the flow fractions of the hot
and the cold branch in the exchanger: a branch at fraction f has its
parent's states at f times the parent's enthalpy flow, so the duty
moves its parent-equivalent enthalpy by ``Q / f``. The duty position q
in [0, Q] (from the hot inlet and the cold outlet end) lies at
``H_hot_in - q / fh`` and ``H_cold_out - q / fc``, a knot H at
``(H_hot_in - H) fh`` or ``(H_cold_out - H) fc``; the states returned
are parent-equivalent (what `_refine_knots` takes). At f = 1 every
expression is the unsplit one, bit for bit (division and
multiplication by 1. are exact).
The positions are the ends and every knot and curve breakpoint inside
the exchanger. Between two consecutive positions both knot curves are
linear, so the planned approach (on the knots) is too, and each stream's
exact states are within `_curve_tol_T` of its knots: an interval whose
two ends have a planned approach of at least `T_min_app` plus that
margin is feasible on the exact states. Every other interval is
checked with `StreamCurve.T_exact` at its ends and midpoint and, where
the exact approach is not linear (curved single-phase stretches and
glides), searched for a minimum inside it (`_interval_min`).
Returns the smallest approach found [K] (exact where evaluated) and the
exact states ``(T_hot, H_hot, T_cold, H_cold)`` wherever the exact
approach is below the pair's required approach less the tolerance,
``T_min_app - _APPROACH_TOL``. With ``T_min_app = inf`` every interval
is checked, so the approach returned is the exact minimum over the
exchanger and the states are all those evaluated.
This is the only exact internal-approach check of the synthesis (the
exchangers themselves, `HXprocess`, check their two terminals only,
which misses an internal pinch at a phase change).
"""
hot, cold = curves[h], curves[c]
H_hot_out, H_cold_out = H_hot_in - Q / fh, H_cold_in + Q / fc
qs = [0., Q]
for H in (knots[h][1], hot.H - hot.H_lo):
qs.extend((H_hot_in - H[(H > H_hot_out) & (H < H_hot_in)]) * fh)
for H in (knots[c][1], cold.H - cold.H_lo):
qs.extend((H_cold_out - H[(H > H_cold_in) & (H < H_cold_out)]) * fc)
qs = np.unique(np.clip(qs, 0., Q))
planned = (_knot_T(knots[h], H_hot_in - qs / fh, True)
- _knot_T(knots[c], H_cold_out - qs / fc, False))
near = planned < T_min_app + _curve_tol_T(hot) + _curve_tol_T(cold) + 1e-9
limit = T_min_app - _APPROACH_TOL
points = []
states = {}
def exact(q):
if q not in states:
qh, qc = q / fh, q / fc
T_hot = hot.T_exact(hot.H_lo + H_hot_in - qh, 'low')
T_cold = cold.T_exact(cold.H_lo + H_cold_out - qc, 'high')
states[q] = dT = T_hot - T_cold
if dT < limit:
points.append((T_hot, H_hot_in - qh, T_cold, H_cold_out - qc))
return states[q]
worst = float(planned[~near].min()) if not near.all() else np.inf
if qs.size == 1:
if near[0]: worst = min(worst, exact(qs[0]))
return worst, points
for k in np.flatnonzero(near[:-1] | near[1:]):
a, b = qs[k], qs[k + 1]
fa, fm, fb = exact(a), exact(0.5 * (a + b)), exact(b)
worst = min(worst, fa, fm, fb)
if abs(fa + fb - 2. * fm) > 1e-9: # not linear: look for a dip
worst = min(worst, _interval_min(exact, a, fa, b, fb))
return worst, points
def _required_approach(required, h, c):
"""The required approach [K] of hot stream `h` against cold stream `c`:
`required` itself if it is a scalar, else ``required[h] - required[c]``
(`required` the streams' shifts, ``grid['shift']``)."""
if np.ndim(required) == 0: return required
return float(required[h] - required[c])
def _check_required(required, knots):
"""Raise a ValueError unless `required` is a scalar or has one shift per
stream of `knots`."""
if np.ndim(required) and len(required) != len(knots):
raise ValueError(f'{len(required)} shifts given for {len(knots)} '
'streams')
def _exact_approach(plan, duties, ends, curves, knots, required):
"""
`_exchanger_approach` of every exchanger in `duties` (duty by index into
``plan.exchangers``) at the enthalpies of the walk `ends` (see `_walk`),
with the flow fractions of its branches (1 on a trunk), against its
pair's required approach: `required` is the scalar pair requirement, or
the per-stream shift vector ``grid['shift']`` (never
``grid['contributions']``; see `_required_approach`).
Returns the smallest approach [K], the violating exact states by stream,
``{stream: [(T_exact, H - H_lo), ...]}``, and the violating exchangers.
"""
_check_required(required, knots)
worst = np.inf
violations = {}
bad = []
for n, Q in duties.items():
e = plan.exchangers[n]
h, c = e.hot, e.cold
approach, points = _exchanger_approach(
curves, knots, h, c, ends[n, h][0], ends[n, c][0], Q,
_required_approach(required, h, c), e.hot_frac, e.cold_frac
)
worst = min(worst, approach)
if points: bad.append(n)
for T_hot, H_hot, T_cold, H_cold in points:
violations.setdefault(h, []).append((T_hot, H_hot))
violations.setdefault(c, []).append((T_cold, H_cold))
return worst, violations, bad
def _approach_margin(plan, duties, ends, curves, knots, required):
"""
The smallest approach less its pair's required approach [K] over the
exchangers in `duties`, as `_exact_approach` evaluates it (`required`
as there).
"""
margin = np.inf
for n, Q in duties.items():
e = plan.exchangers[n]
h, c = e.hot, e.cold
req = _required_approach(required, h, c)
approach = _exchanger_approach(
curves, knots, h, c, ends[n, h][0], ends[n, c][0], Q, req,
e.hot_frac, e.cold_frac
)[0]
margin = min(margin, approach - req)
return margin
def _shrink(curves, knots, h, c, H_hot_in, H_cold_in, Q, T_min_app,
fh=1., fc=1.):
"""
Largest duty ``Q' <= Q`` (to 1e-9 of Q) at which the exchanger of
`_exchanger_approach` (branches at flow fractions `fh` and `fc`) keeps
the pair's required approach ``T_min_app`` less `_APPROACH_TOL` on the
exact states. With both inlets
fixed, a smaller duty lowers the cold stream's enthalpy (so its
temperature) at every position and shortens the exchanger, so the
approach can only grow: the feasible duties form an interval [0, Q']
and bisection finds its end. On branches too: duty position q lies at
``H_cold_in + (Q - q) / fc`` on the cold branch, which falls with Q,
and at ``H_hot_in - q / fh`` on the hot one, whatever Q.
"""
def ok(x):
return not _exchanger_approach(curves, knots, h, c, H_hot_in,
H_cold_in, x, T_min_app, fh, fc)[1]
lo, hi = 0., Q
# a first guess from the local heat capacity flow rates saves most of
# the bisection: the violations are within the chord error of the knots
approach, _ = _exchanger_approach(curves, knots, h, c, H_hot_in,
H_cold_in, Q, T_min_app, fh, fc)
CP = 0. # of the branches: f times the parent's
for j, f in ((h, fh), (c, fc)):
T, H = knots[j]
dT = np.diff(T)
dH = np.diff(H)
slopes = dH[dT > 0.] / dT[dT > 0.]
if slopes.size: CP = max(CP, f * float(slopes.max()))
guess = Q - 2. * (T_min_app - approach) * CP
if 0. < guess < Q and ok(guess): lo = guess
while hi - lo > 1e-9 * Q:
mid = 0.5 * (lo + hi)
if ok(mid): lo = mid
else: hi = mid
return lo
def _repair(plan, duties, knots, is_hot, curves, required):
"""
Shrink every exchanger that falls short of its pair's required approach
less `_APPROACH_TOL` on the exact states (`required` as in
`_exact_approach`) to its largest feasible duty (`_shrink`); the rest of
its duty goes to the utilities. Shrinking a match moves the later stages
of both its streams toward their inlets, which never reduces another
exchanger's approach (the curves are monotone), so one pass suffices;
the loop only guards against rounding. A branch exchanger shrinks with
its flow fractions; its later branch stages and the mix (at the split
enthalpy -/+ the surviving branch duties, see `_walk`) move toward the
inlet too. Returns the new duties and the changes, ``[(n, Q_before,
Q_after)]``.
"""
duties = dict(duties)
changes = []
for _ in range(len(duties) + 1):
ends = _walk(plan, duties, knots, is_hot)[0]
bad = _exact_approach(plan, duties, ends, curves, knots, required)[2]
if not bad: break
for n in bad:
e = plan.exchangers[n]
ends = _walk(plan, duties, knots, is_hot)[0]
Q = _shrink(curves, knots, e.hot, e.cold, ends[n, e.hot][0],
ends[n, e.cold][0], duties[n],
_required_approach(required, e.hot, e.cold),
e.hot_frac, e.cold_frac)
changes.append((n, duties[n], Q))
duties[n] = Q
return duties, changes
def _refine_knots(knots, violations):
"""
Return `knots` with the exact states of `violations` inserted (or, at an
existing knot, i.e. a chord point inside a glide, corrected), clamped
between the neighbouring knot temperatures so that every curve stays
monotone.
"""
knots = list(knots)
for j, points in violations.items():
T, H = knots[j]
T, H = T.tolist(), H.tolist()
for T_new, H_new in sorted(points, key=lambda p: p[1]):
i = int(np.searchsorted(H, H_new))
if i < len(H) and H[i] == H_new:
lo = T[i - 1] if i > 0 else T_new
hi = T[i + 1] if i + 1 < len(T) else T_new
T[i] = min(max(T_new, lo), hi)
else:
lo = T[i - 1] if i > 0 else T_new
hi = T[i] if i < len(T) else T_new
T.insert(i, min(max(T_new, lo), hi))
H.insert(i, H_new)
knots[j] = (np.array(T), np.array(H))
return knots
def _enthalpy_limit(curve, s_in, H, hot):
"""
`H` if `HXprocess` can take it as the enthalpy limit of the stream of
`curve` entering an exchanger in state `s_in`, else None (the other
stream's limit then sets the duty).
`HXprocess` flashes the stream to the limit and rejects an equilibrium
state on the wrong side of the inlet temperature (a heated stream
colder than its inlet). On a monotone curve that happens only strictly
inside a non-equilibrium end jump (see `StreamCurve.jumps`). A
point-load stream's (non-monotone curve's) equilibrium states are not
ordered with its real inlet temperature: e.g. a liquid fed above its
bubble point and boiled to its dew point, colder than its feed, has
every state past its real inlet on the wrong side, its outlet included.
It enters its first exchanger at equilibrium at its inlet enthalpy
instead (see `_first_inlet`), from which its states are ordered, unless
that flash failed. So the equilibrium state at `H` is compared with the
inlet directly.
"""
if curve.monotone:
tol = curve.tol_H
inside = any(H_a + tol < H < H_b - tol for H_a, H_b in curve.jumps)
return None if inside else H
try:
T = curve.state_at_H(H).T
except Exception:
return None
past = T <= s_in.T + _T_SIDE if hot else T >= s_in.T - _T_SIDE
return H if past else None
def _first_inlet(stream, point_load, T_point, hot):
"""
Bring `stream`, a copy of a process stream's real inlet, in place to
the state in which the stream enters its first process exchanger: the
real inlet, except that a point-load stream (non-monotone
`StreamCurve`, whose whole duty is planned at its outlet temperature
`T_point`; `hot` if it is cooled) enters at equilibrium at its inlet
enthalpy and pressure, which lies on the plan's side of `T_point`
(see `hensmith._curves._point_load_inlet`; kept as given if that flash
fails). The enthalpy is the same either way, so no balance changes.
`HXprocess` judges a match by the inlet temperatures (the hotter inlet
is the hot stream, no heat moves unless they are more than `dT` apart,
and the partner's outlet is capped at the inlet temperature -/+ `dT`),
so a point-load stream's real inlet, on the wrong side of `T_point` by
definition, would make it refuse or cut short a match that the plan
keeps the required approach for at `T_point`. E.g. a reboiler fed as a liquid
above its boiling point enters as the vapor-liquid mixture it flashes
to, and a vapor fed below its dew point (e.g. the ideal-thermo copy of
a saturated vapor whose ideal dew point is higher) as the mixture it
partially condenses to, hotter than its feed.
"""
if point_load: stream.copy_like(_point_load_inlet(stream, T_point, hot))
class _RealizationError(Exception):
"""An exchanger of the plan could not be simulated."""
def __init__(self, n, ID, error):
super().__init__(n, ID, error)
self.n, self.ID, self.error = n, ID, error
class _SplitRealizationError(Exception):
"""A splitter or the mixer of split `k` (index into ``plan.splits``)
could not be simulated; `branch_exchangers` holds its branch exchangers
as ``(index into plan.exchangers, ID)``."""
def __init__(self, k, branch_exchangers, error):
super().__init__(k, branch_exchangers, error)
self.k, self.branch_exchangers = k, branch_exchangers
self.error = error
[docs]
class StreamSplit:
"""
A split of one process stream into parallel branches that re-join in
a mixer, as realized by `synthesize_network` (``info['splits']``).
Attributes
----------
stream : int
The stream's index.
side : {'above', 'below'}
The side of the pinch.
index : int
Ordinal of the split among the stream's splits on its side, from 1
in flow order (its IDs end in ``_<index>`` from 2 on).
fractions : tuple[float]
Flow fraction of every branch; they sum to 1 within round-off.
splitters : list[biosteam.Splitter]
The splitter chain, IDs ``Split_<stream>_<hs|cs>[_<index>]``
(element 1), then that plus ``_b<c>`` (element c): element c sends
``f_c / (f_c + ... + f_n)`` of its feed to its first outlet (branch
c - 1) and the rest, by its second outlet, to element c + 1; the
last element's second outlet is the last branch.
mixer : biosteam.Mixer
The rigorous mixer where the branches re-join (inlet b is branch
b), ID ``Mix_<stream>_<hs|cs>[_<index>]``.
branches : list[list[HXprocess]]
The process exchangers of every branch, in flow order. A branch
whose exchangers were all dropped is empty: its splitter outlet
feeds the mixer directly (a bypass).
position : int
Number of the stream's trunk process exchangers before the split.
isothermal : bool
Every branch ends at the stream's state at `H_mix`: the branches
re-join at one temperature.
H_split, H_mix : float
Planned full-flow enthalpies of the stream where it splits and
where it re-joins [kJ/hr].
"""
__slots__ = ('stream', 'side', 'index', 'fractions', 'splitters',
'mixer', 'branches', 'position', 'isothermal', 'H_split',
'H_mix')
def __init__(self, stream, side, index, fractions, splitters, mixer,
branches, position, isothermal, H_split, H_mix):
self.stream, self.side, self.index = stream, side, index
self.fractions, self.splitters, self.mixer = (fractions, splitters,
mixer)
self.branches, self.position = branches, position
self.isothermal, self.H_split, self.H_mix = isothermal, H_split, H_mix
def outlet(self, b):
"""``(splitter, port)`` where branch `b` leaves the chain."""
return ((self.splitters[b], 0) if b < len(self.splitters)
else (self.splitters[-1], 1))
def __repr__(self):
fractions = ', '.join(f'{f:.4g}' for f in self.fractions)
return (f'<StreamSplit {self.splitters[0].ID}: stream {self.stream} '
f'{self.side}, {len(self.fractions)} branches ({fractions})>')
def _walk(plan, duties, knots, is_hot):
"""
Enthalpies (relative to each stream's H_lo) at which every stream enters
and leaves each of its exchangers in `duties` (duty by index into
``plan.exchangers``), walking it in flow order from its inlet
(``plan.stages``), and at which it enters its utility; plus the pair
index of every exchanger (its rank among the exchangers of the same
(side, hot, cold) pair, in the order the hot stream meets them). A
smaller duty (a dropped or shrunk match) shifts the later stages of both
streams toward their inlets.
With stream splitting (``plan.splits``), every stream is walked along
``plan.paths``: a split's branches start at the split enthalpy and a
branch exchanger of flow fraction f moves the branch by its duty over
f, so its `ends` are parent-equivalent (the enthalpies of the full flow
in the branch's state); the stream re-joins at the split enthalpy -/+
the sum of the surviving branch duties, which the fractions leave free
of round-off. A dropped branch exchanger shortens its branch (a branch
left without exchangers bypasses), and later stages move toward the
inlet by the duty lost (`_split_nodes` gives the split enthalpies).
"""
ends = {}
last = []
if not plan.splits:
for j, hot in enumerate(is_hot):
H = knots[j][1][-1] if hot else 0.
for n in plan.stages[j]:
if n not in duties: continue
Q = duties[n]
H_next = H - Q if hot else H + Q
ends[n, j] = (H, H_next)
H = H_next
last.append(H)
else:
for j, hot in enumerate(is_hot):
H = knots[j][1][-1] if hot else 0.
for item in plan.paths[j]:
if isinstance(item, Split):
H = _walk_split(item, duties, ends, H, hot)
continue
if item not in duties: continue
Q = duties[item]
H_next = H - Q if hot else H + Q
ends[item, j] = (H, H_next)
H = H_next
last.append(H)
count = {}
pair_index = {}
for j, hot in enumerate(is_hot):
if not hot: continue
for n in plan.stages[j]:
if n not in duties: continue
e = plan.exchangers[n]
key = (e.side, e.hot, e.cold)
count[key] = pair_index[n] = count.get(key, 0) + 1
return ends, last, pair_index
def _branch_duty(split, duties):
"""Sum (exactly rounded) of the duties of `split`'s branch exchangers
in `duties`."""
return math.fsum([duties[n] for ns in split.branches for n in ns
if n in duties])
def _walk_split(split, duties, ends, H, hot):
"""
`_walk` over `split` (a `hensmith._splitting.Split` of a stream that is
cooled if `hot`) from the split enthalpy `H`: fill the `ends` of its
branch exchangers in `duties`, each branch from `H` by duty over
fraction, and return the enthalpy where the branches re-join.
"""
j = split.stream
for f, ns in zip(split.fractions, split.branches):
Hb = H
for n in ns:
if n not in duties: continue
Q = duties[n] / f
H_next = Hb - Q if hot else Hb + Q
ends[n, j] = (Hb, H_next)
Hb = H_next
D = _branch_duty(split, duties)
return H - D if hot else H + D
def _split_nodes(plan, duties, ends):
"""
Where the streams of a plan with splits (``plan.splits``) split and
re-join, for the exchangers in `duties` and their walk `ends` (see
`_walk`, whose arithmetic this repeats).
Returns
-------
first_nodes : dict[int, tuple[int]]
For every stream, the first exchanger of every branch of the split
that is its first node (the first item of ``plan.paths[j]`` with an
exchanger in `duties`), in branch order: they all take the stream's
real inlet (`_realize`), as that split's splitter chain does
(`_realize_splits`). Empty if the first node is not a split.
split_ends : list[tuple or None]
For every split of ``plan.splits``, ``(H_split, [H_end, ...],
H_mix)``: the enthalpies (relative to the stream's H_lo) where it
splits, where each branch ends (a parent-equivalent enthalpy;
``H_split`` for a branch without exchangers, which bypasses) and
where the branches re-join. None for a split whose branches have
no exchanger left: it is not realized, the stream passes it
unchanged.
first_splits : dict[int, int]
For every stream whose first node is a split, that split's index
into ``plan.splits`` (the one rule for both of the above).
"""
split_ends = []
for split in plan.splits:
j = split.stream
live = [[n for n in ns if n in duties] for ns in split.branches]
n = next((ns[0] for ns in live if ns), None)
if n is None:
split_ends.append(None)
continue
H = ends[n, j][0]
D = _branch_duty(split, duties)
H_mix = H - D if plan.exchangers[n].hot == j else H + D
split_ends.append((H, [ends[ns[-1], j][1] if ns else H
for ns in live], H_mix))
order = {id(split): k for k, split in enumerate(plan.splits)}
first_nodes, first_splits = {}, {}
for j, path in plan.paths.items():
first_nodes[j] = ()
for item in path:
if not isinstance(item, Split):
if item in duties: break
continue
k = order[id(item)]
if split_ends[k] is not None:
first_nodes[j] = tuple(
next(n for n in ns if n in duties)
for ns in item.branches if any(n in duties for n in ns))
first_splits[j] = k
break
return first_nodes, split_ends, first_splits
def _discard(units):
"""Remove units, and the streams connected to them, from the registry of
the active flowsheet (after a failed realization)."""
for unit in units:
for s in (*unit.ins, *unit.outs): bst.main_flowsheet.stream.discard(s)
bst.main_flowsheet.unit.discard(unit)
def _realize(plan, duties, curves, knots, streams_inlet, is_hot, required):
"""
Build and run one plain `HXprocess` per exchanger in `duties` (duty by
index into ``plan.exchangers``), in plan order, each with ``dT`` its
pair's required approach less `_APPROACH_TOL` (`required` as in
`_exact_approach`). Raise
`_RealizationError` (after discarding the units built so far) if one of
them cannot be simulated. Returns ``(units, first, last)``: the units by
exchanger index, each stream's first exchanger (None if it has none)
and the enthalpy at which it enters its utility (relative to its H_lo).
A branch exchanger of a split stream (flow fraction f) runs f of the
stream's flow in the states of the full stream: its inlet is the full
stream's state at the branch's (parent-equivalent) inlet enthalpy,
scaled by f, and its enthalpy limit f times the full stream's at the
planned outlet (`_enthalpy_limit` judges the full-flow state). The
first exchanger of every branch of a split at a stream's inlet takes
the real inlet (`_split_nodes`).
"""
_check_required(required, knots)
ends, last, pair_index = _walk(plan, duties, knots, is_hot)
first = [next((n for n in plan.stages[j] if n in duties), None)
for j in range(len(is_hot))]
first_nodes = _split_nodes(plan, duties, ends)[0] if plan.splits else {}
units = {}
for n in sorted(duties):
e = plan.exchangers[n]
h, c = e.hot, e.cold
dT = _required_approach(required, h, c) - _APPROACH_TOL
suffix = '' if pair_index[n] == 1 else f'_{pair_index[n]}'
if e.side == 'above':
ID = f'HX_{c}_{h}_hs{suffix}'
ports = (c, h)
fractions = (e.cold_frac, e.hot_frac)
else:
ID = f'HX_{h}_{c}_cs{suffix}'
ports = (h, c)
fractions = (e.hot_frac, e.cold_frac)
ins, outs, H_lims = [], [], []
for j, f in zip(ports, fractions):
curve = curves[j]
H_in, H_out = ends[n, j]
if first[j] == n or n in first_nodes.get(j, ()):
s = _copy(streams_inlet[j]) # the real inlet state
_first_inlet(s, not curve.monotone, curve.T_out, is_hot[j])
else:
s = curve.state_at_H(curve.H_lo + H_in)
s.ID = f's_{j}__{ID}'
# a planned outlet whose equilibrium state is not past the inlet
# (inside a non-equilibrium end jump, or a point load's): leave
# it to the other stream's limit
H_lim = _enthalpy_limit(curve, s, curve.H_lo + H_out, is_hot[j])
if f != 1.: # a branch: f of the flow, in the same states
s.scale(f)
if H_lim is not None: H_lim *= f
ins.append(s)
outs.append(s.copy(f'{ID}__s_{j}'))
H_lims.append(H_lim)
hx = bst.HXprocess(ID=ID, ins=ins, outs=outs, H_lim0=H_lims[0],
H_lim1=H_lims[1], dT=dT, thermo=ins[0].thermo)
units[n] = hx
try:
hx._run()
except Exception as error:
_discard(units.values())
raise _RealizationError(n, ID, error)
return units, first, last
def _realize_splits(plan, duties, units, curves, knots, streams_inlet,
is_hot):
"""
Build and run, once each, the splitter chain and the rigorous mixer of
every split of ``plan.splits`` that keeps an exchanger in `duties`
(duty by index into ``plan.exchangers``), after `_realize` built its
exchangers `units` (by exchanger index). A split without one is not
realized: the stream passes it as a trunk.
The splitter chain's feed is the stream's real inlet if the split is
its first node (as its branches' first exchangers, `_split_nodes`),
else its state at the split enthalpy. The mixer's inlet b carries
fraction f_b of the flow: all in the stream's state at the mix
enthalpy if every branch ends there (an isothermal re-join), else each
in the state at its branch's end (parent-equivalent, see `_walk`). Its
outlet starts at the planned state (the full stream at the mix
enthalpy), which only speeds up the flash: where thermosteam's PH flash
disagrees with the stream's curve (inside some two-phase glides), the
outlet lands off the plan whatever its start, and the check on T
reports it. Such a mixer is reported, not undone: the PH flash at that
enthalpy lands there from every start state (measured), so a process
exchanger of the stream ending there would too, split or not.
Returns
-------
splits : list[StreamSplit]
The realized splits, in the order of ``plan.splits``.
deviations : list[dict]
The mixers whose outlet is off the plan, each as ``dict(ID,
T_plan, T, H_plan, H)``: ``|H - H_plan|`` above `_DUTY_TOL` times
the stream's duty, ``H_plan`` being the sum of its inlets'
enthalpies (the binding check for real thermo, whose rigorous
outlet is H of the converged T), or ``|T - T_plan|`` above
`_MIX_T_TOL`, ``T_plan`` being the planned state's temperature.
Raises
------
_SplitRealizationError
If a splitter or a mixer cannot be simulated, after discarding the
units built for the splits and `units`.
"""
ends = _walk(plan, duties, knots, is_hot)[0]
_, split_ends, first_splits = _split_nodes(plan, duties, ends)
tolQ = plan.info['tolQ']
# each realized split's ordinal (per stream and side, from 1 in flow
# order) and the stream's trunk exchangers before it
order = {id(split): k for k, split in enumerate(plan.splits)}
where = {}
for path in plan.paths.values():
count, trunk = {}, 0
for item in path:
if not isinstance(item, Split):
trunk += item in duties
continue
k = order[id(item)]
if split_ends[k] is None: continue
count[item.side] = count.get(item.side, 0) + 1
where[k] = (count[item.side], trunk)
splits, deviations, built = [], [], []
def run(unit, k, live):
try:
unit._run()
except Exception as error:
branch_exchangers = [(n, units[n].ID) for ns in live for n in ns]
_discard([*built, *units.values()])
raise _SplitRealizationError(k, branch_exchangers, error)
for k, split in enumerate(plan.splits):
if split_ends[k] is None: continue
H_split, H_ends, H_mix = split_ends[k]
j, fractions = split.stream, split.fractions
curve = curves[j]
index, position = where[k]
base = f"Split_{j}_{'hs' if split.side == 'above' else 'cs'}"
if index > 1: base += f'_{index}'
mix_ID = 'Mix' + base[len('Split'):]
live = [[n for n in ns if n in duties] for ns in split.branches]
if first_splits.get(j) == k: # the stream's first node
s = _copy(streams_inlet[j]) # the real inlet state
_first_inlet(s, not curve.monotone, curve.T_out, is_hot[j])
else:
s = curve.state_at_H(curve.H_lo + H_split)
s.ID = f's_{j}__{base}'
splitters = []
for c in range(1, len(fractions)):
ID = base if c == 1 else f'{base}_b{c}'
# tail sums: no cancellation, and the last element's rest is
# exactly its complement
r = fractions[c - 1] / math.fsum(fractions[c - 1:])
splitter = bst.Splitter(ID, ins=s, outs=(s.copy(f'{ID}__s_{j}_0'),
s.copy(f'{ID}__s_{j}_1')),
split=r, thermo=s.thermo)
built.append(splitter)
splitters.append(splitter)
run(splitter, k, live)
s = splitter.outs[1]
isothermal = all(abs(H - H_mix) <= _ISO_TOL * tolQ for H in H_ends)
ins = []
for b, (f, H) in enumerate(zip(fractions, H_ends)):
s = curve.state_at_H(curve.H_lo + (H_mix if isothermal else H))
s.ID = f's_{j}_{b}__{mix_ID}'
s.scale(f)
ins.append(s)
outlet = curve.state_at_H(curve.H_lo + H_mix) # the planned state
outlet.ID = f'{mix_ID}__s_{j}'
T_plan = outlet.T
mixer = bst.Mixer(mix_ID, ins=ins, outs=outlet, rigorous=True,
thermo=outlet.thermo)
built.append(mixer)
run(mixer, k, live)
out = mixer.outs[0]
H_plan = math.fsum([s.H for s in mixer.ins])
if (abs(out.H - H_plan) > _DUTY_TOL * abs(curve.H_out - curve.H_in)
or abs(out.T - T_plan) > _MIX_T_TOL):
deviations.append(dict(ID=mix_ID, T_plan=T_plan, T=out.T,
H_plan=H_plan, H=out.H))
splits.append(StreamSplit(
j, split.side, index, tuple(fractions), splitters, mixer,
[[units[n] for n in ns] for ns in live], position, isothermal,
curve.H_lo + H_split, curve.H_lo + H_mix))
return splits, deviations
[docs]
def synthesize_network(hus, T_min_app=5., Qmin=1e-3, force_ideal_thermo=False,
avoid_recycle=False, sort_hus_by_T=False, info=None,
stream_splitting=False, approach_contributions=None):
"""
Synthesize a heat exchanger network for the process streams behind a
set of utility heat exchangers: pinch analysis (`problem_table`), then a
pinch-outward plan that reaches the minimum energy requirement (MER)
targets whenever the search finds an unsplit network that does,
realized with one `HXprocess` per match and one rigorous `HXutility`
per stream. By default no stream is split; with `stream_splitting`, a
side of the pinch that no unsplit network serves at MER is planned with
stream splits and reaches the targets.
Parameters
----------
hus : list[HeatUtility]
One heat utility per process stream; `hu.unit` is the original heat
exchanger (its `ins[0]`/`outs[0]` are the stream's end states) and
the sign of `hu.duty` marks the stream: positive = heated (cold
stream), negative = cooled (hot stream); zero-duty utilities are
dropped. Heating utilities are placed before cooling utilities;
within each group the given order is kept unless `sort_hus_by_T`.
All returned per-stream arrays and lists are indexed in that
rearranged order (the stream index), which also breaks ties in the
planner's search.
T_min_app : float, optional
Minimum approach temperature [K]: kept on the exact stream states
at both ends of and everywhere inside every process exchanger, and
used to shift hot streams in the problem table. With
`approach_contributions`, the default approach: each stream without
a contribution contributes ``T_min_app / 2``. Defaults to 5.
Qmin : float, optional
Planned exchangers with a duty below this [kJ/hr] are dropped and
their duty left to the utilities (a large value can cost MER).
Defaults to 1e-3.
force_ideal_thermo : bool, optional
Analyze copies of the streams with ideal thermodynamics
(``thermo.ideal()``); the synthesized exchangers inherit that
thermo. Defaults to False.
avoid_recycle : bool, optional
Never match the same (hot, cold) pair twice anywhere, so no two
exchangers connect the same pair of streams (a second exchanger
between them can form a recycle loop in the network). This
forbids the repeated matches that some unsplit MER networks need.
Defaults to False.
sort_hus_by_T : bool, optional
Sort the heating utilities by inlet temperature, descending, and the
cooling utilities ascending, before analysis. Defaults to False.
info : dict, optional
If given, filled with the synthesis report: 'status' ('mer' if the
realized network's utilities equal the targets, else
'best_effort'), 'Q_hot_target' and 'Q_cold_target' (the problem
table's targets), 'Q_hot_plan' and 'Q_cold_plan' (the planned
utilities), 'Q_hot' and 'Q_cold' (the utilities of the realized
network, from the simulated exchanger duties), 'penalty'
(``Q_hot_plan - Q_hot_target``), 'sides' (per side of the pinch:
status, method, work, proof of a needed split, gaps, units),
'plan_targets' (the planner's own cascade in the first round:
Q_hot, Q_cold, pinch_T on the shifted scale, cut),
'refine_rounds', 'min_approach' (the smallest approach inside any
process exchanger [K]; exact on the stream states wherever it is
within the curves' linearization tolerance of its pair's required
approach, else from the knots; None without process exchangers),
'min_approach_margin' (the smallest approach less its pair's
required approach over the process exchangers [K], as evaluated by
the exact check: at least -1e-6 up to round-off; None when
'min_approach' is None), 'approach_contributions' (each stream's
approach temperature contribution as used [K], unspecified ones at
``T_min_app / 2``, in stream order: that of the returned arrays,
not of `hus`), 'deviations' (exchangers whose simulated duty
differs from the plan), 'qmin_dropped' (matches dropped by
`Qmin`), 'repaired' (matches shrunk to keep their pair's required
approach on the exact states, see Notes), 'dropped' (matches that
could not be simulated;
normally empty) and 'point_loads' (the indices of the streams
whose outlet temperature does not move with their duty, e.g. an
isothermal condenser or a reboiler fed as a liquid above its
boiling point, so that their whole duty is a point load at the
outlet temperature; each enters its first process exchanger at
equilibrium at its inlet enthalpy, see Notes). With
`stream_splitting` also 'stream_splitting' (True), 'splits' (the
realized splits, a list of `StreamSplit`) and 'split_deviations'
(mixers whose outlet is off the planned state; normally empty);
each side in 'sides' then also has 'split' (None for a side that
did not try to split, else the chosen candidate and its network
signature, both None if no candidate was found, see
`hensmith._planner.Plan`) and its method is 'split-<candidate>'
when it splits. Required with `stream_splitting`, which raises a
ValueError without it (the splitters and mixers are returned only
there).
stream_splitting : bool, optional
Allow a process stream to be split into parallel branches that
re-join. A side of the pinch that no unsplit network serves at the
minimum energy requirement (a pinch-rule proof, or an unsplit
search that leaves a utility penalty) is planned with splits
instead, and then reaches the MER targets exactly on the planner's
knots (see Notes, "Stream splitting"). Sides that an unsplit
network serves are never split, so a problem that needs no split
gets the same network as with the default. The branches are
realized with `biosteam.Splitter` chains and rigorous
`biosteam.Mixer` units, which are adiabatic and cost nothing; their
structure is reported in ``info['splits']`` (`synthesize_network`)
and ``synthesis_info['splits']`` (`HeatExchangerNetwork`). With
`avoid_recycle`, a split that would repeat a stream pair is not
used, and MER is then not guaranteed. Defaults to False.
approach_contributions : Mapping, Sequence or callable, optional
Approach temperature contribution of each process stream [K]: a
match between hot stream h and cold stream c must keep ``c_h +
c_c`` instead of `T_min_app` (see Notes, "Approach
contributions"). A value is None or nan (unspecified: ``T_min_app
/ 2``) or a finite number > 0 (`T_min_app` must then be finite and
> 0 too). Given as one of:
- a mapping ``{key: value}`` whose keys name units: a unit (that
very object; a stale auxiliary that its owner rebuilt, no longer
among the units of `hus` and their owners, names the one of the
same ID under owners of the same IDs instead, if any),
``(owner.ID, unit.ID)`` with the top owner (the unit's own ID
twice for a unit without an owner; it names every auxiliary of
that ID under that owner, at any depth, so the heaters of two
auxiliary flash drums are told apart by keying the drums or the
heater objects) or the ID of a unit without an owner (a bare
auxiliary ID such as ``'condenser'`` names nothing). Only the
units of `hus` and their owners are at hand, so an auxiliary that
its unit holds without setting its owner is a unit without an
owner here: an intercooler that a `biosteam.MultistageCompressor`
K1 built anew is ``'K1_H1'``, and ``'K1'`` raises (a
`HeatExchangerNetwork` that K1 is available to takes K1 as its
owner instead). A heat utility takes the value of
the key nearest to its unit along the unit's chain of owners: a
unit's own key beats its owner's, which beats the owner's owner's
(e.g. ``'D1'`` covers the condenser and the reboiler of column
D1, and ``('D1', 'condenser')`` overrides it for the condenser).
A None or nan value wins at its level like any other. Heat
utility keys are rejected (BioSTEAM replaces them at every
simulation). A key that names neither the unit of a heat utility
of `hus` nor one of its owners, nor an auxiliary that one of them
may build but does not hold now (such a key has no effect: e.g.
``('X1', 'heat_exchangers[2]')``, biosteam's ID of the third
exchanger of a `biosteam.HXutilities` X1, while X1 needs fewer,
or an old exchanger of X1 that none took the place of), or two
keys of the same level for a unit with different values, raise a
ValueError;
- a sequence (or `numpy.ndarray`; not a str or bytes) with one
value per heat utility, aligned with `hus` as given (before the
heating utilities are placed first; zero-duty utilities
included); a wrong length raises a ValueError;
- a callable ``f(hu)`` returning the value of a heat utility,
called once per heat utility in the order of `hus`.
Anything else raises a TypeError. Checked and resolved before any
stream is copied or flashed. If every value is unspecified, the
result is the default one. The contributions as used are reported
in ``info['approach_contributions']`` in stream order, not in the
order of `hus`. Defaults to None (every stream contributes
``T_min_app / 2``).
Returns
-------
HXs_hot_side : list[HXprocess]
Process exchangers of the hot-side (above-pinch) design, in plan
order (from the pinch outward), IDs ``HX_<cold>_<hot>_hs``;
``ins``/``outs`` [0] is the cold stream and [1] the hot stream.
HXs_cold_side : list[HXprocess]
Process exchangers of the cold-side (below-pinch) design, in plan
order, IDs ``HX_<hot>_<cold>_cs``; ``ins``/``outs`` [0] is the hot
stream and [1] the cold stream. The n-th exchanger (n >= 2) of the
same pair on the same side, counted in the order the hot stream
meets them, gets the suffix ``_<n>`` (e.g. ``HX_3_2_cs_2``).
new_HX_utils : list[HXutility]
One rigorous utility exchanger per stream (possibly of zero duty)
bringing it from its last process exchanger to its outlet enthalpy,
IDs ``Util_<index>_cs`` (hot streams) / ``Util_<index>_hs`` (cold
streams); listed hot streams first.
hxs : list[Unit]
The original heat exchangers, in stream order.
T_in_arr, T_out_arr : numpy.ndarray
Inlet and (quenched) outlet temperatures of each stream [K].
pinch_T_arr : numpy.ndarray
Per-stream pinch temperature [K] (informational): the process pinch
on the stream's own scale when the stream crosses it
(`ProblemTable.pinch_T` plus the stream's shift: `T_min_app` for a
hot stream, 0 for a cold one by default, ``T_min_app/2 + c_j`` and
``T_min_app/2 - c_j`` with approach contributions); the inlet
temperature of a stream whose inlet already
lies past the pinch in its direction of flow, or that is isothermal
or non-monotone; the outlet temperature of a stream that ends
before reaching the pinch.
C_flow_vector : numpy.ndarray
Heat capacity flow rate of each process stream, ``|H_out - H_in| /
|T_in - T_out|`` [kJ/hr/K] from its inlet and quenched outlet (the
temperature difference is replaced by 1e-12 for an isothermal
stream, which therefore ranks as a very large flow rate).
hx_utils_rearranged : list[HeatUtility]
The heat utilities of `hus` in stream order.
streams_inlet : list[Stream]
One copy of each stream's inlet, in stream order, as prepared for
the analysis (ideal-thermo copies if `force_ideal_thermo`). The
network works on further copies, so these keep their inlet state.
stream_HXs_dict : dict[int, list[Unit]]
For each stream index, its process exchangers in flow order, then
its utility exchanger. With splits, the order is topological: the
exchangers before a split, those of its branches (branch by
branch, each in flow order), then those after it.
hot_indices, cold_indices : list[int]
Stream indices of the hot and cold streams.
Notes
-----
*Curves and targets.* Every stream's outlet is quenched to equilibrium
at its own enthalpy and described by a piecewise-linear temperature-
enthalpy curve (see `problem_table`), built once. The problem table
[Kemp07]_ on the union of all breakpoints gives the targets, the pinch
and the side of the pinch that point loads at the pinch temperature
belong to. The planner models each stream by its knots on that grid,
which reproduces the table's cascade exactly.
*Planner.* Each stream is cut at the pinch; above it the hot streams,
below it the cold streams must be served completely by process matches
("musts"), while the partners ("flexes") leave any remainder to a
utility at their far end. A depth-first search builds each side from
the pinch outward, one match at a time. Every match keeps its pair's
required approach (`T_min_app` by default) at every knot (so internal
pinches, e.g. a condensing vapor against a boiling mixture, are
respected), and every step keeps the problem table of the remaining
problem feasible (remaining problem analysis
[Smith05]_, as a closed-form bound on the duty). At every level where
that table is tight, the pinch design rules [LH83]_ (see also
[Seider17]_, Chapter 9; number and heat-capacity-flow rules,
generalized to isothermal segments, which can serve several partners
in series) must hold; at the pinch itself a violation proves that MER
needs stream splitting.
Candidate duties are the largest feasible one and a finite set of
events (a stream ticked off, a partner saved for another stream, a
switch of partner, a return). The same pair may be matched repeatedly,
which emulates a split by series alternation. Budgets are counted in
deterministic work units, so results do not depend on machine speed.
A branch and bound then reduces the number of exchangers. A side that
is proven to need splits, or whose search runs out of budget, gets a
best-effort plan (unless `stream_splitting`, see "Stream splitting"
below): heat a must cannot place is moved to its pinch end,
where it crosses the pinch at the cost of an equal amount of extra hot
and cold utility (the penalty), minimized by greedy dives and a
bisection of these gaps. See `hensmith._planner` for the details.
*Realization.* Each stream is walked in flow order from its inlet: a
hot stream through its hot-side matches from its inlet end, then its
cold-side matches, then its cooler; a cold stream through its
cold-side matches, then its hot-side matches, then its heater.
Each match becomes a plain `HXprocess` whose two enthalpy limits
(`H_lim0`, `H_lim1`) are the planned outlet enthalpies, so that its
duty reproduces the plan; its `dT` is its pair's required approach
(`T_min_app` by default) less 1e-6 K, only a guard against rounding
(the approach is enforced by the plan and the check below). A planned
outlet whose equilibrium state is not past the stream's state at the
exchanger inlet cannot be a limit (`HXprocess`
rejects it): one strictly inside a non-equilibrium end jump (see
`StreamCurve.jumps`), or on a point-load stream colder (hotter) than
its inlet state when heated (cooled). That stream's limit is left out
and the other stream's sets the duty; a match in which neither stream
can take a limit runs to its `dT` guard (a deviation if its duty
differs). A stream's first exchanger gets its real inlet, except a
point-load stream (whose outlet temperature does not move with its
duty, so the plan places its whole duty there, a temperature its real
inlet lies beyond): it enters at equilibrium at its inlet enthalpy,
which lies on the plan's side of its outlet temperature (a reboiler
fed as a liquid above its boiling point enters as the mixture it
flashes to; a vapor fed below its dew point, e.g. under
`force_ideal_thermo`, as the mixture it partially condenses to),
because `HXprocess` compares the inlet temperatures with `dT` and
would refuse or cut short a match planned there. Later exchangers get
the stream's exact state at the planned enthalpy.
Every exchanger is simulated once; one whose duty differs from the
plan is reported in `info`. Each stream ends in one rigorous
`HXutility` to its outlet enthalpy; an `AssertionError` is raised if
it does not reproduce the quenched outlet within tolerance.
*Exactness.* The knots are exact at grid points but chords in between
(at most 0.002 K off inside glides and curved single-phase stretches).
Every planned exchanger is therefore checked on the exact stream states
wherever its planned approach is within that margin of its pair's
required approach. Where a MER plan falls short by more than 1e-6 K, the
exact states are inserted as knots and the network is planned again (at
most three rounds). A best-effort plan (or a MER plan still short after
the last round) instead has each violating match shrunk to the largest
duty that keeps the approach, the rest going to the utilities: with its
inlets fixed a smaller duty can only raise a match's approach, and the
later stages of both streams move toward their inlets, which never
reduces another match's approach. Constant heat capacity streams never
need either step.
*Stream splitting.* With `stream_splitting`, a side that no unsplit
network serves at MER is planned with splits (see
`hensmith._splitting`). A branch exchanger of flow fraction f runs f of
its stream's flow in the full stream's states: its inlet is the
stream's state at the branch's inlet enthalpy scaled by f, and its
enthalpy limit f times the full stream's; the exact check follows each
branch on its parent's curve with enthalpy steps of duty over f. Each
split becomes a chain of `biosteam.Splitter` units (one fewer than its
branches, fed with the real inlet if the split is the stream's first
node) and a rigorous `biosteam.Mixer`, whose outlet starts at the
planned state; a mixer outlet off that state is reported in
``info['split_deviations']``. Each refine round re-plans a split
side's previous pick first and keeps it while it plans the same
network (its signature, with the fractions the refined knots move). If
exchangers on a split side still violate after the last refine round,
that network is excluded and the side re-planned (at most two more
rounds): with another network if the side has one left, else with the
same one on the refined knots (a side's last candidate is never
excluded), which ends the retries. If that does worse, the best MER
plan of the rounds is restored. A split whose splitters or mixer
cannot be simulated becomes a trunk (its branch exchangers are
dropped, as in ``info['dropped']``).
*Guarantees and limits.* The utilities are never below the targets.
Every process exchanger keeps its pair's required approach (`T_min_app`
by default) less 1e-6 K on the exact states at its ends and at every
checked position inside it, and the
heat balance closes on every stream. 'mer' is reported only if the
realized network reaches the targets. The result is deterministic.
Completeness is empirical, not proven: every pruning test is a
necessary condition, so a missed MER network can only come from the
finite set of candidate duties, the caps on repeated pairs or the work
budgets; the planner reached MER on every unsplit-feasible problem of a
certified benchmark of about 1,700 problems with 2-40 streams. Networks
whose match order is cyclic cannot be represented. An unsplit MER
network can need many exchangers (series alternation approaches a
split only in the limit); MER always takes precedence over the number
of units. By default, problems that need stream splits get a
best-effort network whose penalty is small but not minimal in general,
and a side that needs splits without a pinch-rule proof spends its
whole MER budget before the best-effort step.
With `stream_splitting` (and without `avoid_recycle`), every side that
no unsplit network serves has a split plan at MER on the planner's
knots: the vertical core of `hensmith._splitting` always yields one
(Theorem M there), up to the part of the root deficit that the
cascade's own tolerances already absorbed into the targets (at most
about 1e-9 of the total duty). The remaining gap is the realization.
With constant heat capacities the knots are exact, and the realized
network reaches the targets to round-off. With real thermodynamics
the knots are chords, so a split plan can fall short of the required
approach on the exact states by up to the chords' tolerance; the refine
rounds close that as they do for unsplit plans, with up to two more rounds
that exclude the violating network (see "Stream splitting" above), and
only an exchanger still short after them is shrunk (reported in
'repaired', with status 'best_effort'). A rigorous mixer's outlet
enthalpy differs from the sum of its inlets by at most its flash
residual, which the stream's utility takes. The test suite reaches
'mer' with splits on all 38 problems of its corpus that provably need
them (25 with constant heat capacities, 13 with real thermodynamics)
and on the three such regression systems. The number of exchangers is
minimized only among the candidate split plans a side generates, and
splitters and mixers are not costed.
*Approach contributions.* With `approach_contributions`, each stream j
contributes ``c_j`` (``T_min_app / 2`` unless given), and a match
between hot stream h and cold stream c must keep ``c_h + c_c`` instead
of `T_min_app` [Kemp07]_: e.g. a small contribution for a condensing or
boiling stream (high film coefficients) and a large one for a gas.
Stream j is shifted down by ``s_j = T_min_app/2 + c_j`` if hot and
``T_min_app/2 - c_j`` if cold, so that ``T*_h >= T*_c`` if and only if
``T_h - T_c >= c_h + c_c``: the problem table, the targets and the
planner work on one shared shifted scale, and every guarantee above
carries over. A pair's required approach, ``s_h - s_c``, is enforced on
the exact states by the exactness check and the repair, is each
exchanger's `dT` guard (less 1e-6 K), and ``info['min_approach_margin']``
reports the smallest margin over it. Stream j's pinch temperature is
the shifted pinch plus ``s_j``. Only the process matches follow the
contributions: each stream's `HXutility` gets its agent from BioSTEAM's
own rule. With every contribution unspecified, the shifts are the
default ones bit for bit (``x/2 + x/2 == x``), and so is the result.
Thermosteam's TP flashes fail silently inside the glides of some
mixtures (e.g. water-ethanol with 20-50 % ethanol); an exchanger
simulated there can deviate from its plan (reported in
`info['deviations']`), and a mixer's outlet can land off its planned
state (reported in ``info['split_deviations']``).
`HeatExchangerNetwork` calls this function, rewires each stream's
stages in series (the branches of a split in parallel, from its
splitter chain to its mixer), converges the network as a `System` and
costs it.
Examples
--------
Problem r002: one hot stream against three cold ones (heat capacity
flow rates in kW/K, temperatures 300 K above those of the classic
problem, T_min_app = 10 K). Its only unsplit MER network matches the
hot stream twice with the same cold stream. A constant heat capacity
pseudo-component makes 1000 kmol/hr of fluid per kW/K:
>>> import biosteam as bst, thermosteam as tmo
>>> from hensmith.hxn_synthesis import synthesize_network
>>> Fluid = tmo.Chemical('Fluid', search_db=False, phase='l', MW=1.,
... Cn=3.6, default=True)
>>> bst.settings.set_thermo([Fluid], cache=True)
>>> def process_stream(ID, T_in, T_out, CP):
... inlet = bst.Stream(ID + '_in', Fluid=1000. * CP, T=T_in,
... units='kmol/hr')
... hx = bst.HXutility(ID, ins=inlet, T=T_out, rigorous=False)
... hx.simulate()
... return hx
>>> units = [process_stream('C1', 440., 470., 1.),
... process_stream('C2', 400., 420., 1.),
... process_stream('C3', 420., 550., 2.),
... process_stream('H1', 520., 350., 4.)]
>>> hus = [hx.heat_utilities[0] for hx in units]
>>> info = {}
>>> result = synthesize_network(hus, T_min_app=10., info=info)
>>> HXs_hot_side, HXs_cold_side, new_HX_utils = result[:3]
>>> for hx in HXs_hot_side + HXs_cold_side:
... print(hx.ID, round(hx.Q / 3600., 6), 'kW')
HX_3_2_cs 160.0 kW
HX_3_0_cs 30.0 kW
HX_3_2_cs_2 20.0 kW
HX_3_1_cs 20.0 kW
>>> info['status']
'mer'
The utilities equal the MER targets, 80 kW of heating and 450 kW of
cooling:
>>> duties = [(hx.outs[0].H - hx.ins[0].H) / 3600. for hx in new_HX_utils]
>>> round(sum(Q for Q in duties if Q > 0), 6), round(-sum(Q for Q in duties if Q < 0), 6)
(80.0, 450.0)
>>> round(info['Q_hot_target'] / 3600., 6), round(info['Q_cold_target'] / 3600., 6)
(80.0, 450.0)
With an approach temperature contribution of 2 K for the hot stream
(keyed by its unit) and the default ``T_min_app / 2 = 5`` K for the
cold streams, every match needs 7 K instead of 10 K, so the targets
fall to those of ``T_min_app = 7`` K, and every exchanger's `dT`
guard is 7 K less 1e-6 K:
>>> info = {}
>>> result = synthesize_network(hus, T_min_app=10., info=info,
... approach_contributions={units[3]: 2.})
>>> info['status'], round(info['Q_hot_target'] / 3600., 6), round(info['Q_cold_target'] / 3600., 6)
('mer', 74.0, 444.0)
>>> info['approach_contributions']
array([5., 5., 5., 2.])
>>> sorted({round(hx.dT, 6) for hx in result[0] + result[1]})
[6.999999]
References
----------
.. [LH83] Linnhoff, B., & Hindmarsh, E. (1983). The pinch design method
for heat exchanger networks. Chemical Engineering Science, 38(5),
745-763.
.. [Kemp07] Kemp, I. C. (2007). Pinch Analysis and Process Integration
(2nd ed.). Butterworth-Heinemann.
.. [Smith05] Smith, R. (2005). Chemical Process Design and Integration.
Wiley.
.. [Seider17] Seider, W. D., Lewin, D. R., Seader, J. D., Widagdo, S.,
Gani, R., & Ng, M. K. (2017). Product and Process Design Principles.
Wiley. Heat Exchanger Networks (Chapter 9).
"""
if stream_splitting and info is None:
raise ValueError('stream_splitting=True needs an info dict: the '
'splitters and mixers are returned in '
'info["splits"]')
# resolved before any stream is copied or flashed (None: the default)
values = _explicit_contributions(hus, approach_contributions, T_min_app)
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, table, curves, grid = \
_pinch_analysis(hus, T_min_app, force_ideal_thermo, sort_hus_by_T,
values)
# every pair's required approach is shift[h] - shift[c]; by default it
# is the scalar T_min_app itself, passed as given (not the shift vector,
# whose float() pair difference in `_required_approach` would promote a
# numpy float32 T_min_app to float64 and change the guards), so the
# default checks, repair and realization are unchanged to the bit; the
# planner gets the contributions only if any is given, likewise
shift = T_min_app if values is None else grid['shift']
contributions = None if values is None else grid['contributions']
N = len(hxs)
is_hot = [False] * N
for i in hot_indices: is_hot[i] = True
H_in_arr = np.array([c.H_in for c in curves])
H_out_arr = np.array([c.H_out for c in curves])
dTs = np.abs(T_in_arr - T_out_arr)
C_flow_vector = np.abs(H_out_arr - H_in_arr) / np.maximum(dTs, 1e-12)
duty = np.abs(H_out_arr - H_in_arr)
scale = float(duty.sum())
# Plan on the grid knots and check the plan on the exact states. Where it
# falls short, a MER plan is planned again on knots refined with the
# exact states (at most _MAX_REFINE rounds); a best-effort plan, or a MER
# plan still short after the last round, has its violating matches
# shrunk locally instead (a re-plan of a best-effort side repeats its
# whole search to recover a duty of the order of the chord error).
# With stream splitting, every round re-plans a split side's previous
# pick first (stickiness); after the last refine round, up to
# _MAX_SPLIT_RETRY more rounds run while violating exchangers lie on
# split sides whose networks can still be excluded (_same_network), and
# if they end worse, the best MER plan of the rounds is restored.
knots = _grid_knots(curves, grid)
exclude, prefer, best = {}, {}, None
n_rounds = _MAX_REFINE + (_MAX_SPLIT_RETRY if stream_splitting else 0)
for refine_round in range(n_rounds + 1):
plan = plan_network(knots, is_hot, T_min_app,
avoid_recycle=avoid_recycle, Qmin=Qmin,
**(dict(stream_splitting=True,
_split_exclude=exclude,
_split_prefer=prefer)
if stream_splitting else {}),
**({} if contributions is None else
dict(approach_contributions=contributions)))
if not refine_round: plan_targets = plan.info['cascade']
duties = {n: e.Q for n, e in enumerate(plan.exchangers)}
ends = _walk(plan, duties, knots, is_hot)[0]
min_approach, violations, bad = _exact_approach(
plan, duties, ends, curves, knots, shift
)
if stream_splitting and plan.status == 'mer' and (
best is None or len(bad) < best[0]):
best = (len(bad), plan, knots, duties, min_approach, bad)
if (not violations or plan.status != 'mer'
or refine_round == n_rounds): break
if refine_round >= _MAX_REFINE:
# a retry (stream splitting only): exclude the networks of the
# split sides of the violating exchangers; if none is new, stop
# where the default loop stops
grew = False
for n in bad:
side = plan.exchangers[n].side
split = plan.info['sides'][side]['split']
if split is None or split['candidate'] is None: continue
signature = split['signature']
if not any(_same_network(signature, excluded)
for excluded in exclude.get(side, ())):
exclude.setdefault(side, set()).add(signature)
grew = True
if not grew: break
if stream_splitting:
prefer = {}
for side, side_info in plan.info['sides'].items():
split = side_info['split']
if split is not None and split['candidate'] is not None:
prefer[side] = (split['candidate'], split['signature'])
knots = _refine_knots(knots, violations)
if (refine_round > _MAX_REFINE and best is not None
and (plan.status != 'mer' or len(bad) > best[0])):
# the retries did worse: back to the best MER plan (the fewest
# violating exchangers, the earliest round), with its knots
_, plan, knots, duties, min_approach, bad = best
repaired = []
if bad:
duties, changes = _repair(plan, duties, knots, is_hot, curves,
shift)
for n, Q_before, Q_after in changes:
e = plan.exchangers[n]
repaired.append(dict(side=e.side, hot=e.hot, cold=e.cold,
Q_plan=Q_before, Q=Q_after))
if Q_after < Qmin: del duties[n]
ends = _walk(plan, duties, knots, is_hot)[0]
min_approach = _exact_approach(plan, duties, ends, curves, knots,
shift)[0]
# the smallest approach less its pair's required approach: by default
# every pair's is T_min_app (rounding is monotone, so this is the
# smallest per-exchanger margin); else measured per pair, on fresh walk
# ends (`ends` is stale after a best plan was restored)
if not duties:
min_approach_margin = None
elif contributions is None:
min_approach_margin = min_approach - T_min_app
else:
min_approach_margin = _approach_margin(
plan, duties, _walk(plan, duties, knots, is_hot)[0], curves,
knots, shift)
# Realize; a match that cannot be simulated is dropped (its duty goes to
# the utilities: removing a match never reduces another's approach), and
# a split whose splitters or mixer cannot be simulated becomes a trunk
# (its branch exchangers are dropped). Each failure removes at least one
# exchanger, so the loop ends.
dropped = []
while True:
try:
units, first, last = _realize(plan, duties, curves, knots,
streams_inlet, is_hot, shift)
splits, split_deviations = (
_realize_splits(plan, duties, units, curves, knots,
streams_inlet, is_hot)
if plan.splits else ([], [])
)
except _RealizationError as failure:
dropped.append(dict(ID=failure.ID, Q=duties.pop(failure.n),
error=repr(failure.error)))
except _SplitRealizationError as failure:
for n, ID in failure.branch_exchangers:
if n in duties:
dropped.append(dict(ID=ID, Q=duties.pop(n),
error=repr(failure.error)))
else:
break
HXs_hot_side = [units[n] for n in sorted(units)
if plan.exchangers[n].side == 'above']
HXs_cold_side = [units[n] for n in sorted(units)
if plan.exchangers[n].side == 'below']
deviations = []
for n, hx in units.items():
e = plan.exchangers[n]
if abs(hx.Q - duties[n]) > _DUTY_TOL * (duty[e.hot] + duty[e.cold]):
deviations.append(dict(ID=hx.ID, Q_plan=duties[n], Q=hx.Q))
# One rigorous utility per stream, hot streams first
stream_HXs_dict = {i: [units[n] for n in plan.stages[i] if n in units]
for i in indices}
new_HX_utils = []
for i in hot_indices + cold_indices:
hot = is_hot[i]
curve = curves[i]
ID = 'Util_%s_cs'%i if hot else 'Util_%s_hs'%i
if first[i] is None:
s = _copy(streams_inlet[i])
else:
s = curve.state_at_H(curve.H_lo + last[i])
s.ID = 's_%s__%s'%(i, ID)
outlet = s.copy('%s__s_%s'%(ID, i))
new_HX_util = bst.units.HXutility(ID=ID, ins=s, outs=outlet,
H=H_out_arr[i], rigorous=True,
thermo=s.thermo)
new_HX_util._run()
s_out = new_HX_util.outs[0]
atol_T = 5. if 's' in hxs[i].outs[0].phases else 0.001
if hot:
np.testing.assert_allclose(s_out.H, H_out_arr[i], rtol=5e-3, atol=1.)
np.testing.assert_allclose(s_out.T, T_out_arr[i], rtol=5e-3, atol=atol_T)
else:
np.testing.assert_allclose(s_out.H, H_out_arr[i], rtol=1e-2, atol=1.)
np.testing.assert_allclose(s_out.T, T_out_arr[i], rtol=5e-2, atol=atol_T)
new_HX_utils.append(new_HX_util)
stream_HXs_dict[i].append(new_HX_util)
if info is not None:
# utilities of the plan and of the realized network: each stream's
# duty less what its process exchangers transfer
planned = {n: duties[n] for n in units}
simulated = {n: hx.Q for n, hx in units.items()}
loads = []
for Q_of in (planned, simulated):
Q_hot = Q_cold = 0.
for i in indices:
remaining = duty[i] - sum(Q_of[n] for n in plan.stages[i]
if n in units)
if is_hot[i]: Q_cold += remaining
else: Q_hot += remaining
loads.append((Q_hot, Q_cold))
(Q_hot_plan, Q_cold_plan), (Q_hot, Q_cold) = loads
tol = _ACHIEVED_TOL * scale
mer = (plan.status == 'mer'
and abs(Q_hot - table.hot_util_load) <= tol
and abs(Q_cold - table.cold_util_load) <= tol)
info.update(
status='mer' if mer else 'best_effort',
Q_hot_target=table.hot_util_load,
Q_cold_target=table.cold_util_load,
Q_hot_plan=Q_hot_plan, Q_cold_plan=Q_cold_plan,
Q_hot=Q_hot, Q_cold=Q_cold,
penalty=Q_hot_plan - table.hot_util_load,
sides=plan.info['sides'],
plan_targets=dict(Q_hot=plan_targets['Qh'],
Q_cold=plan_targets['Qc'],
pinch_T=plan_targets['pinch_T'],
cut=plan_targets['cut']),
refine_rounds=refine_round,
min_approach=min_approach if duties else None,
min_approach_margin=min_approach_margin if duties else None,
approach_contributions=grid['contributions'].copy(),
deviations=deviations,
qmin_dropped=plan.info['qmin_dropped'],
dropped=dropped, repaired=repaired,
point_loads=[i for i in indices if not curves[i].monotone],
)
if stream_splitting:
info.update(stream_splitting=True, splits=splits,
split_deviations=split_deviations)
return HXs_hot_side, HXs_cold_side, new_HX_utils, hxs, T_in_arr,\
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.
Exchangers left out of `hxs` pass the precedence on: two requested
exchangers of a stream are ordered whenever the stream flows from one
to the other. A life cycle with splits orders its exchangers by
`StreamLifeCycle.stage_pairs`, so the exchangers of sibling branches
are not ordered by that stream.
"""
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:
if getattr(life_cycle, 'splits', None):
pairs = _split_precedence(life_cycle, position)
else:
stages = [i.unit for i in life_cycle.life_cycle if i.unit in position]
if not life_cycle.cold: stages.reverse()
pairs = zip(stages, stages[1:])
for a, b in pairs:
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 _split_precedence(life_cycle, requested):
"""
Pairs ``(a, b)`` of `requested` exchangers of a life cycle with splits
where column `a` goes left of column `b`: the stream flows from `a` to
`b` (from `b` to `a` for a hot stream), through `stage_pairs` and the
stages not requested. Each requested stage is paired with the nearest
requested stages downstream, which gives the same precedence as their
transitive closure.
"""
successors = {}
for a, b in life_cycle.stage_pairs():
successors.setdefault(id(a), []).append(b)
pairs = []
for stage in life_cycle.life_cycle:
if stage.unit not in requested: continue
seen, stack = set(), list(successors.get(id(stage), ()))
while stack:
other = stack.pop()
if id(other) in seen: continue
seen.add(id(other))
if other.unit in requested: pairs.append((stage.unit, other.unit))
else: stack.extend(successors.get(id(other), ()))
if not life_cycle.cold: pairs = [(b, a) for a, b in pairs]
return pairs
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`:
'<owner> - <auxiliary name> (<inlet stream ID>)', 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
[docs]
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 (of its whole flow, `StreamLifeCycle.H_in` and
its utility's outlet). Exchanger columns on each side of the pinch are
ordered so that each stream meets them in flow direction whenever the
network allows it; the exchangers of a split stream's branches are
ordinary columns, with their branch duties, and sibling branches are
not ordered by that stream (splits are not drawn). Stream labels read
'<unit> - <auxiliary> (<stream>)' 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:<ID>', 'Util:<ID>', 'Label:<index>',
# 'H_in:<index>', 'H_out:<index>') 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
# the whole stream's inlet, also where its first stage is a branch
H_in = life_cycle.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)
)
gid_left, gid_right = ('H_in:', 'H_out:') if cold else ('H_out:', '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,
gid=gid_left + str(index), **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,
gid=gid_right + str(index), **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
