District Heating Network Hydraulic Flow¶
Author: Ruizhi Yu
This example solves the hydraulic subproblem of a district heating system (DHS): given each node’s mass injection demand, compute the steady-state pipe mass flow distribution so that mass is conserved at every node and the pressure drop around every loop is zero.
The example is written in two equivalent forms — one element-wise and
one using Mat_Mul — so it serves both as a worked walkthrough of the
matrix-vector calculus pipeline and as a regression test against
SolUtil’s DhsFlow on the
well-known Barry Island 35-node network.
Governing equations¶
Let \(n_\text{node}\) be the number of nodes, \(n_\text{pipe}\) the number of pipes, and \(n_\text{loop} = n_\text{pipe} - (n_\text{node} - 1)\) the number of independent loops. For Barry Island, the counts are \(n_\text{node} = 35\), \(n_\text{pipe} = 35\), \(n_\text{loop} = 1\).
Mass continuity¶
Let \(V \in \mathbb{R}^{n_\text{node} \times n_\text{pipe}}\) be the signed node–pipe incidence matrix:
Then mass balance at every node reads \(V\,m = m^\text{inj}\), where \(m \in \mathbb{R}^{n_\text{pipe}}\) is the pipe-flow vector and \(m^\text{inj} \in \mathbb{R}^{n_\text{node}}\) is the signed node injection (negative at sources/slacks, positive at loads, zero at intermediate nodes).
For a connected network the rows of \(V\) satisfy \(\mathbf{1}^\top V = 0\), so one equation is redundant. We drop the slack row and keep \(n_\text{node} - 1\) independent balance equations.
Loop pressure¶
Let \(L \in \mathbb{R}^{n_\text{loop} \times n_\text{pipe}}\) be the loop–pipe incidence matrix: \(L_{\ell j} = \pm 1\) if pipe \(j\) belongs to loop \(\ell\), with sign giving the reference direction. For each loop, the head loss around the closed contour is zero:
where \(K \in \mathbb{R}^{n_\text{pipe}}\) is the per-pipe quadratic resistance coefficient and \(|\cdot|\) is element-wise absolute value. The \(m \odot |m| = m^2 \cdot \operatorname{sign}(m)\) form handles signed flow directions correctly.
Mat_Mul form¶
Putting the two together, the Mat_Mul model has one vector equation of length \(n_\text{node} - 1\) and one vector equation of length \(n_\text{loop}\):
def build_matmul_model(hc: dict, m_init: np.ndarray, m_inj: np.ndarray):
"""Build the hydraulic subproblem using matrix-vector form.
The two governing equations are a single vector equation each:
V_ns @ m - m_inj_ns = 0 (one vector eqn of length n_node-1)
L @ (K ⊙ m ⊙ |m|) = 0 (one vector eqn of length n_loop)
where ``V_ns`` is ``V`` with the slack-node row removed and
``m_inj_ns`` the matching slice of the injection vector. ``L`` is
the loop incidence built from ``hc['pinloop']``.
The Mat_Mul form produces analytical Jacobian blocks of the shapes
Solverz's matrix-calculus engine recognises:
- ``∂(V_ns @ m) / ∂m = V_ns`` — constant-matrix derivative.
- ``∂(L @ (K ⊙ m ⊙ |m|)) / ∂m = L · (diag(K ⊙ |m|) + diag(K ⊙ m ⊙
sign(m)))`` — a *mutable matrix* Jacobian block (depends on m),
triggering the ``Diag(v)@L`` scatter-add code path.
"""
model = Model()
model.m = Var('m', m_init)
# --- Signed node-pipe incidence V (sparse), one slack row removed ---
n_node = hc['n_node']
n_pipe = hc['n_pipe']
V_dense = np.zeros((n_node, n_pipe))
for edge in hc['G'].edges(data=True):
fnode, tnode, data = edge
pipe = data['idx']
V_dense[tnode, pipe] = 1 # flows INTO tnode
V_dense[fnode, pipe] = -1 # flows OUT of fnode
# Drop exactly one row — see :func:`build_elementwise_model` for the
# rationale. Multi-slack systems drop only the first slack row.
slack_nodes = sorted(hc['slack_node'].tolist())
skip_node = slack_nodes[0] if slack_nodes else 0
non_slack_rows = [n for n in range(n_node) if n != skip_node]
V_ns = csc_array(V_dense[non_slack_rows, :])
model.V_ns = Param('V_ns', V_ns, dim=2, sparse=True)
model.m_inj_ns = Param('m_inj_ns', m_inj[non_slack_rows])
model.K = Param('K', hc['K'])
# --- Equation 1: mass continuity on non-slack nodes ---
model.mass_balance = Eqn(
'mass_balance',
Mat_Mul(model.V_ns, model.m) - model.m_inj_ns)
# --- Equation 2: loop pressure drop, one row per independent loop ---
# Radial networks contribute no loop equations at all — the balance
# system alone is square in that case.
pinloop = np.atleast_2d(np.asarray(hc['pinloop']))
if pinloop.shape[1] != n_pipe and pinloop.shape[0] == n_pipe:
pinloop = pinloop.T
# Strip any all-zero rows (radial sub-loops) before passing to the
# model so the resulting L has exactly the loops we actually need.
nontrivial_rows = [i for i in range(pinloop.shape[0])
if np.any(pinloop[i] != 0)]
if nontrivial_rows:
L_sparse = csc_array(pinloop[nontrivial_rows].astype(np.float64))
model.L = Param('L', L_sparse, dim=2, sparse=True)
model.loop_pressure = Eqn(
'loop_pressure',
Mat_Mul(model.L, model.K * model.m * Abs(model.m)))
return model
The symbolic Jacobian derived by the matrix-calculus engine is:
The first block is a constant matrix Jacobian (just \(V_\text{ns}\)). The second is a mutable matrix Jacobian — it depends on \(m\) — and Solverz decomposes it into scatter-add loops at code-generation time.
Element-wise form¶
For comparison, the same system written with one scalar equation per node and per loop:
def build_elementwise_model(hc: dict, m_init: np.ndarray, m_inj: np.ndarray):
"""Build the hydraulic subproblem using one scalar ``Eqn`` per node /
per loop.
Parameters
----------
hc : dict
DHS case dictionary (from ``SolUtil.sysparser.load_hs``) with
fields ``n_node``, ``n_pipe``, ``G`` (networkx DiGraph), ``K``,
``pinloop``, ``slack_node``.
m_init : ndarray (n_pipe,)
Initial guess for the pipe mass flow vector.
m_inj : ndarray (n_node,)
Node injection vector (signed). ``V @ m`` should equal ``m_inj``.
Returns
-------
Model
A Solverz ``Model`` ready for ``create_instance``.
"""
model = Model()
model.m = Var('m', m_init)
model.K = Param('K', hc['K'])
model.m_inj = Param('m_inj', m_inj)
slack_nodes = sorted(hc['slack_node'].tolist())
n_node = hc['n_node']
n_pipe = hc['n_pipe']
# Mass continuity: we drop exactly ONE redundant row. For a
# connected graph with all node injections prescribed the balance
# system ``V m = m_inj`` has a single linearly dependent row
# (summing every row gives ``0 = 0``), no matter how many slack
# nodes are declared. Dropping one slack row is sufficient and
# necessary; dropping one row per slack would make the system
# underdetermined when there are multiple slacks.
skip_node = slack_nodes[0] if slack_nodes else None
for node in range(n_node):
if node == skip_node:
continue
rhs = -model.m_inj[node]
for edge in hc['G'].in_edges(node, data=True):
pipe = edge[2]['idx']
rhs = rhs + model.m[pipe]
for edge in hc['G'].out_edges(node, data=True):
pipe = edge[2]['idx']
rhs = rhs - model.m[pipe]
model.__dict__[f'mass_balance_{node}'] = Eqn(
f'mass_balance_{node}', rhs)
# Loop pressure: one scalar equation per independent loop.
# ``hc['pinloop']`` can be a 1-D row vector (single loop) or a 2-D
# array shaped ``(n_loop, n_pipe)``. Radial networks
# (``n_loop == 0``) get no loop equation — the balance rows alone
# are sufficient.
pinloop = np.atleast_2d(np.asarray(hc['pinloop']))
if pinloop.shape[1] != n_pipe and pinloop.shape[0] == n_pipe:
pinloop = pinloop.T
for loop_idx in range(pinloop.shape[0]):
loop_row = pinloop[loop_idx]
if not np.any(loop_row != 0):
# Purely radial row — skip, no pressure constraint.
continue
rhs = 0
for j in range(n_pipe):
coeff = int(loop_row[j])
if coeff != 0:
rhs = rhs + coeff * model.K[j] * model.m[j] * Abs(model.m[j])
model.__dict__[f'loop_pressure_{loop_idx}'] = Eqn(
f'loop_pressure_{loop_idx}', rhs)
return model
Both formulations produce numerically identical pipe mass flows.
Regression test¶
The test test_heat_flow_hydraulic_matmul_regression solves the
hydraulic system via all four code paths through Solverz
(inline/module × element-wise/Mat_Mul) and checks that each one
matches the ground truth from SolUtil’s DhsFlow:
def test_heat_flow_hydraulic_matmul_regression(datadir):
"""Full regression: all four solver paths must agree with SolUtil."""
df, m_inj = _load_ground_truth(datadir)
hc = df.hc
# A flat-but-non-zero initial guess so Newton has something to do.
m_init = 0.5 * np.ones(hc['n_pipe'])
# -------- 1. Element-wise, inline mode --------
mdl_ew = build_elementwise_model(hc, m_init, m_inj)
m_ew_inline, _, _ = _solve(mdl_ew)
np.testing.assert_allclose(
m_ew_inline, df.m, rtol=1e-6, atol=1e-8,
err_msg='element-wise inline disagrees with SolUtil ground truth')
# -------- 2. Mat_Mul, inline mode --------
mdl_mm = build_matmul_model(hc, m_init, m_inj)
m_mm_inline, spf_mm, y0_mm = _solve(mdl_mm)
np.testing.assert_allclose(
m_mm_inline, df.m, rtol=1e-6, atol=1e-8,
err_msg='Mat_Mul inline disagrees with SolUtil ground truth')
# Mat_Mul and element-wise inline must agree to machine precision —
# they solve the same system and start from the same y0.
np.testing.assert_allclose(
m_mm_inline, m_ew_inline, rtol=1e-10, atol=1e-12,
err_msg='Mat_Mul inline diverges from element-wise inline')
# -------- 3. Element-wise, module printer (jit=True) --------
tmp_dir = os.path.join(
os.path.dirname(os.path.abspath(__file__)), '_hf_generated')
if os.path.exists(tmp_dir):
shutil.rmtree(tmp_dir)
mdl_ew_mod = build_elementwise_model(hc, m_init, m_inj)
m_ew_mod = _solve_module(mdl_ew_mod, 'heat_flow_ew_mod', tmp_dir, jit=True)
np.testing.assert_allclose(
m_ew_mod, df.m, rtol=1e-6, atol=1e-8,
err_msg='element-wise module printer disagrees with SolUtil')
# -------- 4. Mat_Mul, module printer (jit=True) --------
mdl_mm_mod = build_matmul_model(hc, m_init, m_inj)
m_mm_mod = _solve_module(mdl_mm_mod, 'heat_flow_mm_mod', tmp_dir, jit=True)
np.testing.assert_allclose(
m_mm_mod, df.m, rtol=1e-6, atol=1e-8,
err_msg='Mat_Mul module printer disagrees with SolUtil')
np.testing.assert_allclose(
m_mm_mod, m_mm_inline, rtol=1e-10, atol=1e-12,
err_msg='Mat_Mul module printer diverges from its inline twin')
shutil.rmtree(tmp_dir)
The second test test_heat_flow_stepwise_jacobian_match is the
strict check that caught the v0.8.0 bug where the module printer
froze its mutable-matrix Jacobian blocks at their initial values. It
drives Newton iteration with the inline model and, at each step,
asserts that the module printer’s \(F\) and \(J\) match the inline
values to machine precision. If the mutable-matrix block assembly
path ever regresses, this test fails immediately on step 0 or 1
instead of only degrading the final solution’s accuracy.
References¶
The Mat_Mul formulation of DHS hydraulic and thermal equations follows:
Yu, Gu, Lu, Yao, Zhang, Ding, Lu. Non-Iterative Calculation of Quasi-Dynamic Energy Flow in the Heat and Electricity Integrated Energy Systems. IEEE Transactions on Power Systems, 38(5), 4148–4162 (2023). doi:10.1109/TPWRS.2022.3210167
equations (11)–(17). SolUtil’s DhsFlow implements the equivalent
element-wise form used as ground truth.