合工大周啸课题组

基于梯度算法的供水管网水力计算实现

论文介绍

供水管网水力分析（稳态流分析）是管网设计与优化的核心环节。其任务是在已知管网拓扑、管道参数、节点需水量和固定水头（水库、水塔）的条件下，求解所有管道的流量 $Q$ 和所有未知节点的测压水头 $H$。传统方法（如 Hardy‑Cross、牛顿法、线性理论、优化法）各有缺陷：或收敛慢、或依赖初值、或方程组规模过大、或数值实现复杂。

作者提出了一种梯度算法（Gradient Algorithm）：

  1. 在流量‑水头联合空间 $(Q, H)$ 中应用牛顿法，并借助 Content Model（内容模型） 的凸性证明了解的存在唯一性——这是算法无条件收敛的理论基石。
  2. 通过分块矩阵消元，将每步牛顿迭代压缩为只求解一个与未知节点水头数 $nn$ 同规模的线性系统，再通过简单回代得到新的管道流量。
  3. 该线性系统具有对称正定、稀疏、Stieltjes 矩阵的良好性质，可采用 不完全 Cholesky 分解 + 共轭梯度法（ICF/MCG） 高效求解，存储量仅与管道数 $np$ 成正比，适合在个人计算机上求解大型管网。

论文（核心）结构：

本文提出的梯度算法后来成为全球广泛使用的供水管网物理模拟软件 EPANET 的核心求解器（称为 Global Gradient Algorithm）。

1 Introduction

引言其实在讲四件事：

这就是整篇论文真正的核心出发点。

2 The Derivation of the Proposed Newton-Raphson Algorithm

2.1 问题定义

已知量

未知量

具体的形式定义

对于管段：

对于节点：

管段与节点的连接关系使用 $A_{12}(i,j)$ （行为管段、列为节点）来描述：

$$ A_{12}(i,j)=\begin{cases} 1 & \text{if flow of pipe i enters node j}\\ 0 & \text{if pipe i and node j are not connected}\\ -1 & \text{if flow of pipe i leaves node j} \end{cases} $$

另外，$A_{21}=A_{12}^T$. 另外还有 $A_{10}$ 是 管道与固定水头节点 的关联矩阵，形状为 $(np,no) $，其元素定义与 $A_{12}$ 完全相同，只是针对的是固定头节点（如水库、水塔）。

2.2 控制方程的矩阵形式

首先来看质量守恒：

$$ A_{21}Q=q $$

形状变换过程为：$[nn, np] × [np, 1] = [nn, 1]$. 其实也就是说，质量守恒的个数 = 未知水头节点数 $nn$.

再来看能量守恒：

$$ A_{12}H+F(Q)=-A_{10}H_0 $$

形状变换过程为：$[np, nn] × [nn, 1] = [np, 1] = [np, no] × [no, 1]$. 其实也就是说，能量守恒的个数 = 未知流量管段数 $np$.

2.3 一个简单的小 case

假设存在一个WDN系统：一个定水头源点 + 两个未知节点 + 两根管，网络如下：

$$ \text{Reservoir }(H_0) \xrightarrow{Q_1} \text{Node 1 }(H_1) \xrightarrow{Q_2} \text{Node 2 }(H_2) $$

设定：

所以，未知水头节点数 $nn=2$，定水头节点数 $no=1$，管段数 $np=2$.

先写质量守恒，对于 Node 1 来说，流入 Node 1 的是 $Q_1$，流出 Node 1 的是 $Q_2$，Node 1 还要消耗 $q_1$。因此，

$$ Q_1-Q_2=q_1 $$

对于 Node 2 来说，流入 Node 2 的是 $Q_2$，没有流出，Node 2 还要消耗 $q_2$。因此，

$$ Q_2=q_2 $$

$$ \begin{cases} Q_1-Q_2=q_1 \\ Q_2=q_2 \end{cases} $$

然后再另：

$$ Q= \begin{bmatrix} Q_1\\ Q_2 \end{bmatrix}, \qquad q= \begin{bmatrix} q_1\\ q_2 \end{bmatrix} $$

那么：

$$ A_{21}= \begin{bmatrix} 1 & -1\\ 0 & 1 \end{bmatrix} $$

注意：这里是 $A_{21}$，$A_{12}$ 中行为管段、列为节点，$A_{21}$ 则刚好相反。

所以：

$$ A_{21}Q= \begin{bmatrix} 1 & -1 \\ 0 & 1 \end{bmatrix} \begin{bmatrix} Q_1 \\ Q_2 \end{bmatrix} = \begin{bmatrix} Q_1-Q_2 \\ Q_2 \end{bmatrix} = \begin{bmatrix} q_1\\ q_2 \end{bmatrix} $$

再来写能量守恒，未知水头向量：

$$ H= \begin{bmatrix} H_1\\ H_2 \end{bmatrix} $$

定水头向量：

$$ H_0= \begin{bmatrix} H_0 \end{bmatrix} $$

水头损失向量：

$$ F(Q)= \begin{bmatrix} f_1(Q_1) \\ f_2(Q_2) \end{bmatrix} $$

那么可以写成：

$$ A_{12} = A_{21}^T = \begin{bmatrix} 1 & 0\\ -1 & 1 \end{bmatrix}, A_{10}= \begin{bmatrix} -1 \\ 0 \end{bmatrix} $$

于是：

$$ A_{12}H+F(Q)=-A_{10}H_0 $$

展开左边：

$$ A_{12}H= \begin{bmatrix} 1 & 0 \\ -1 & 1 \end{bmatrix} \begin{bmatrix} H_1 \\ H_2 \end{bmatrix} = \begin{bmatrix} H_1 \\ -H_1+H_2 \end{bmatrix} $$

右边：

$$ -A_{10}H_0 = - \begin{bmatrix} -1 \\ 0 \end{bmatrix} H_0 = \begin{bmatrix} H_0 \\ 0 \end{bmatrix} $$

所以总式子是：

$$ \begin{bmatrix} H_1 \\ -H_1+H_2 \end{bmatrix} + \begin{bmatrix} f_1(Q_1) \\ f_2(Q_2) \end{bmatrix} = \begin{bmatrix} H_0 \\ 0 \end{bmatrix} $$

这里可能部分读者会向我一样有一个小疑问：在这个 case 中，我们规定了管段流向，但是在真实的 WDN 计算过程中，我们是不知道具体的管段流向的，这个问题我们应该怎么去处理？关于这个问题，这里建立的“参考方向”其实只是为了建立矩阵和方程，不是说你已经知道水一定这么流。随后求解过程得到计算后的流量，如果为负，代表与规定方向相反；如果为正，则说明其与规定方向同向。

2.3 解的存在唯一性证明

凸函数与凸优化

一个优化问题由两部分组成：目标函数 + 约束条件。凸优化要求同时满足：

  1. 目标函数是凸函数

    • 直观理解：函数图像像一个“碗”（开口向上），任意两点连线都在图像上方。
    • 数学定义：$f(\theta x + (1-\theta)y) \le \theta f(x) + (1-\theta) f(y)$，对任意 $x,y$ 和 $\theta \in [0,1]$ 成立。

  2. 约束条件定义的可行域是凸集

    • 线性约束（如 $A_{21}Q = q$）定义的集合是凸集（超平面）。
    • 直观理解：集合内任意两点连线上的所有点仍在集合内。

当这两个条件都满足时，这个优化问题就是凸优化问题。

为什么凸优化问题的解（如果存在）一定是唯一的？

下面通过反证法去理解。假设存在两个不同的最优解 $x^*$ 和 $y^*$（$x^* \ne y^*$），且它们的目标函数值相等且最小：

$$ f(x^*) = f(y^*) = p^* \quad (\text{最小值}) $$

因为可行域是凸集，所以中点 $\frac{x^* + y^*}{2}$ 也在可行域内。

又因为 $f$ 是凸函数，有：

$$ f\left(\frac{x^* + y^*}{2}\right) \le \frac{f(x^*) + f(y^*)}{2} = p^* $$

这说明中点的函数值 不大于 最小值 $p^*$，但它又不可能小于最小值（因为 $p^*$ 是最小值），所以只能等于 $p^*$。

于是我们得到：$f(x^*) = f(y^*) = f(\frac{x^*+y^*}{2}) = p^*$，并且凸函数的不等式必须取等号。

对于严格凸函数，不等式是严格的（除非 $x^* = y^*$）：

$$ f\left(\frac{x^* + y^*}{2}\right) < \frac{f(x^*) + f(y^*)}{2} \quad \text{当 } x^* \ne y^* $$

这就产生了矛盾：左边 $< p^*$，而右边 $= p^*$，意味着存在一个可行点比最小值还小，不可能。

因此假设不成立，只能有 $x^* = y^*$，即解唯一。

凸函数的几何直觉：

后续论文中的 Content Model 目标函数是 严格凸 的（因为 $R_i>0, n_i>0$ 保证了每个管道的贡献项 $ \frac{R_i|Q_i|^{n_i+1}}{n_i+1} $ 是严格凸的，且线性项不破坏严格凸性），加上线性等式约束（凸集），所以是严格凸优化问题。 如果这样的问题有解，解必定唯一。

解的唯一性

由连续性方程和能量守恒方程所表示的系统，其解的个数可能依赖于各个 $f_i(Q_i)$ 的函数形状而不止一个。若所有 $f_i(Q_i)$ 都是单调递增函数，则可以证明该系统的解是存在且唯一的（Pilati and Todini, 1984）。

The system represented by equation (1) may have more than one solution depending upon the shape of $f_i(Q_i)$. If all $f_i(Q_i)$ are monotonically increasing functions, it can be proved that the solution of system (1) exists and is unique (Pilati and Todini, 1984).

如果流量越大，损失也越大，那么系统在物理上就是“正常阻尼型”的，数学上也更容易保证：

而对常见的管道阻力模型，这个条件通常是成立的。

Content Model 优化函数

论文将管网水力平衡问题重写为如下约束优化问题：

$$ \min C(Q) = \sum_{i=1}^{np} R_i \frac{|Q_i|^{n_i+1}}{n_i+1} + \sum_{j=1}^{no} H_{0j}\sum_{i=1}^{np} A_{01}(j,i)Q_i, \; \text{subject to} \sum_{i=1}^{np} A_{21}(j,i)Q_i-q_j=0,(j=1,\ldots,nn) $$

求解原方程组 ⟺ 求解 Content Model 的最优解（KKT 点）。

其中：

上式是对原始管网水力平衡问题的一种等价数学重写。通过构造一个全局目标函数 $C(Q)$，使其在连续约束条件下的一阶最优性条件能够对应原始的网络平衡关系。

目标函数第一项为：

$$ \sum_{i=1}^{np} R_i \frac{|Q_i|^{n_i+1}}{n_i+1} $$

它来源于单管损失函数：$f_i(Q_i)=R_i|Q_i|^{n_i-1}Q_i$ 对 $Q_i$ 的积分。即对每根管道构造单管项：$\Phi_i(Q_i)=\int_0^{Q_i} f_i(\xi),d\xi$，在 Hazen–Williams 形式下可得：$\Phi_i(Q_i)=R_i \frac{|Q_i|^{n_i+1}}{n_i+1}$。

因此，第一项代表了由各单管损失关系积分后得到的单管目标项之和。物理上可理解为管道中储存的“耗散能量”（或称为“内容”），整个管网这项求和，就是所有管道的总耗散内容。另外，后续对 $Q_i$ 求导时，可以重新恢复原来的单管损失函数 $f_i(Q_i)$。

目标函数第二项为：

$$ \sum_{j=1}^{no} H_{0j}\sum_{i=1}^{np} A_{01}(j,i)Q_i $$

该项可紧凑写为：$H_0^T A_{01} Q$，它表示为定水头边界节点通过与其相连的边界管段，对网络施加的边界驱动作用。

具体分析来看，可以这样解读：

  • $A_{01}(j,i)$ 表示固定水头节点 $j$ 与管道 $i$ 的连接关系（+1 流入，-1 流出）。
  • 所以 $\sum_i A_{01}(j,i) Q_i$ 是从固定水头节点 $j$ 流入管网的总流量。
  • 乘以该节点的水头 $H0_j$，就是该节点提供的势能输入功率（或称为“供应能”）。

因为固定水头节点的水头是常数，因此这一项是线性的。另外，在后续对 $Q$ 求导时，定水头边界 $H_0$ 能自然出现在一阶条件中，从而正确恢复包含边界条件的能量关系。

另外，此公式还有一个约束条件：$\sum_{i=1}^{np} A_{21}(j,i)Q_i-q_j=0\:(j=1,\ldots,nn)$. 这就是每一个未知水头节点上的连续方程。这个约束条件的加入，使得此公式成为了：在所有满足节点质量守恒的流量分布中，寻找对应的全局解。

接下来，作者进一步通过拉格朗日乘子法，将带连续约束的优化问题写成拉格朗日函数：

$$ \min L(Q,\lambda)= \sum_{i=1}^{np} R_i \frac{|Q_i|^{n_i+1}}{n_i+1} + \sum_{j=1}^{no} H_{0j}\sum_{i=1}^{np} A_{01}(j,i)Q_i+ \sum_{j=1}^{nn}\lambda_j\sum_{i=1}^{np}(A_{21}(j,i)Q_i-q_j) \qquad (4) $$

其中：

先来通过一个简单的概念来理解下什么是拉格朗日乘子。假设你想让两个数的乘积 $ x \times y $ 最大，但这两个数必须满足 $ x + y = 10 $（比如周长为定值，让面积最大）。

  • 没有约束时，$ x \times y $ 可以无限大。
  • 但有了 $ x + y = 10 $ 这条“路”，你能选的 $ x, y $ 就被限制住了。

拉格朗日的方法就是构造一个新函数：

$$ \mathcal{L}(x, y, \lambda) = x \times y + \lambda \times (10 - x - y) $$

然后把这个新函数当作没有约束来处理，对 $ x, y, \lambda $ 分别求偏导等于0。你会解出 $ x = 5, y = 5 $，这就是最优解。在计算过程中，这个 $ \lambda $ 就像是“如果我把约束放松一点（比如把 10 改成 11），最优值会变多少”的一个指示器。λ 的数值告诉你：约束对目标值的“紧箍”有多强。

这个 Eq.4 的作用是将 Eq.3 中的约束优化问题，改写为一个统一的拉格朗日函数形式，以便通过一阶条件同时恢复连续方程与能量方程。

Eq.4 与 Eq.3 的关系具体来说，Eq.3 为：

$$ \min C(Q) \quad \text{s.t. } A_{21}Q=q $$

这是一个带约束的最小化问题。

Eq.4 则是在此基础上引入拉格朗日乘子 $\lambda_j$，把约束吸收到总函数中，写成：

$$ L(Q,\lambda)=C(Q)+\lambda^T(A_{21}Q-q) $$

因此：

  • Eq.3 是约束优化问题
  • Eq.4 是对应的拉格朗日函数形式

Eq.4 可以看成由三部分组成：

但在本文推导中，$ \lambda $ 不是单纯的数学符号。后续通过一阶条件可以发现：$ \lambda $ 与未知节点水头 $H$ 具有对应关系，且最终会被识别为与节点水头同类的物理量。

公式 Eq.3 虽然已经把原问题写成了约束优化问题，但如果继续直接处理约束，会比较不方便。引入公式 Eq.4 后，可以统一通过偏导来处理：

因此，公式 Eq.4 的意义在于：把原先分开的“目标函数”和“约束条件”统一成一个总函数，为后续推导网络平衡方程的一阶条件形式提供基础。

对 $\lambda_j$ 求导的意义

对 Eq.4 中第 $j$ 个拉格朗日乘子求导，有：

$$ \frac{\partial L}{\partial \lambda_j} = \sum_{i=1}^{np}(A_{21}(j,i)Q_i-q_j) $$

令其为零：

$$ \sum_{i=1}^{np}(A_{21}(j,i)Q_i-q_j)=0 $$

这正好恢复第 $j$ 个未知水头节点上的连续方程。

因此，对乘子求导，对应恢复质量守恒。

对 $Q_i$ 求导的意义

对 Eq.4 中某个管段流量 $Q_i$ 求导时，会同时得到三类贡献：

  1. 来自第一项恢复单管损失函数：$f_i(Q_i)=R_i|Q_i|^{n_i-1}Q_i$
  2. 来自第二项恢复定水头边界项：$\sum_{j=1}^{no} H_{0j}A_{01}(j,i)$, 这是一个常数（因为 $H0_j$ 已知，$A_{01}(j,i)$ 已知），它代表了固定水头节点通过管道 $i$ 施加的“边界能量贡献”。
  3. 来自第三项恢复与该管段相连节点对应的乘子项：$\sum_{j=1}^{nn}\lambda_j A_{21}(j,i)$, 后面我们求出 $ \lambda $ 后，会发现这一项其实就是管道 $i$ 两端未知节点水头的组合。

最后，合成一阶条件（$\partial L / \partial Q_i = 0$）为：

$$ f_i(Q_i) \;+\; \sum_{j=1}^{no} H0_j A_{01}(j,i) \;+\; \sum_{j=1}^{nn} \lambda_j A_{21}(j,i) \;=\; 0 $$

这其实就是管段能量平衡方程。直观理解非常抽象，下面通过一个直观的例子讲解：

假设一根管道 $i$ 连接两个未知节点 $p$ 和 $q$，设参考方向设为从 $p$ 到 $q$，那么对于节点与管段的连接关系（行为节点、列为管段）：

$$ A_{21}(p,i) = -1 \quad (\text{管道离开节点 } p) $$

$$ A_{21}(q,i) = +1 \quad (\text{管道进入节点 } q) $$

$$ A_{21}(其他, i) = 0 $$

未知节点乘子项变为：

$$ \sum_j \lambda_j A_{21}(j,i) = \lambda_p \cdot (-1) + \lambda_q \cdot (+1) = \lambda_q - \lambda_p $$

如果 $\lambda = H$（不要急，先这么理解，后面会证明，让子弹再飞一会儿~），这就是 $H_q - H_p$（终点水头减起点水头）。因此方程变为：

$$ f_i(Q_i) + \text{(边界项)} + (H_q - H_p) = 0 $$

移项：

$$ H_p - H_q = f_i(Q_i) + \text{(边界项)} $$

这正是起点水头 − 终点水头 = 水头损失 + 边界贡献，如果没有边界项（管道不直接连水库），就是 $H_p - H_q = f_i(Q_i)$。

如果管道一端连着固定水头节点 （水源、水塔）$r$，设管道从固定节点 $r$ 到未知节点 $p$（参考方向 $r \to p$）。

  • 固定节点贡献：$\sum_j H0_j A_{01}(j,i) = H0_r \cdot A_{01}(r,i)$，因为管道从 $r$ 流出，所以 $A_{01}(r,i) = -1$（离开固定节点），该项 = $-H0_r$。
  • 未知节点贡献：只有节点 $p$ 有 $A_{21}(p,i) = +1$（流入），所以该项 = $+\lambda_p = +H_p$。

方程：

$$ f_i(Q_i) + (-H0_r) + H_p = 0 \quad\Rightarrow\quad H_p - H0_r = -f_i(Q_i) $$

即 $H0_r - H_p = f_i(Q_i)$，正确（水源到节点的水头差等于损失）。

总结

求导的物理意义：

求导对象得到的条件物理含义
对 $Q_i$$\partial L/\partial Q_i = 0$管道 $i$ 的能量平衡：两端水头差 = 水头损失 + 边界贡献
对 $\lambda_j$$\partial L/\partial \lambda_j = 0$节点 $j$ 的质量守恒：净流入流量 = 需水量

所以 Eq.4 的一阶条件，就是完整的水力方程组。拉格朗日函数就像一个“总能量表达式”，求极值的过程自然导出了能量守恒和质量守恒。这正是变分原理在管网问题中的应用。

2.4 等效性证明

在上一节的 Eq.3 中，已经构造出拉格朗日函数 $L(Q,\lambda)$。本小节的目标是：

通过对 $Q$ 和 $ \lambda $ 分别求偏导并令其为零，建立该约束优化问题的一阶条件系统，并进一步说明该系统与原始管网稳态水力平衡方程组等价。

论文中指出，当：

时，拉格朗日函数 $L$ 具有凸性，因此该问题的解存在且唯一，并满足极小值成立的充分条件。

根据上一节的结论，我们只需另 Eq (4) 一阶偏导等于零，即可得到其唯一解。

$$ \begin{cases} \frac{\partial L}{\partial Q_i}=0,\qquad i=1,\ldots,np \\ \frac{\partial L}{\partial \lambda_j}=0,\qquad j=1,\ldots,nn \end{cases} \qquad (5) $$

Eq.5 的物理与数学意义上一小节也进行了详细说明：

Eq.5 中的偏导条件可以统一写成矩阵形式：

$$ \begin{bmatrix} A_{11} & A_{12} \\ A_{21} & 0 \end{bmatrix} \begin{bmatrix} Q \\ \lambda \end{bmatrix} \begin{bmatrix} -A_{10}H_{0} \\ q \end{bmatrix} \qquad (6) $$

具体的求解过程是这样的。首先另 Eq.4 中的三部分分别记为：$L(Q,\lambda)=T_1+T_2+T_3$，即：

$$ T_1=\sum_{i=1}^{np} R_i \frac{|Q_i|^{n_i+1}}{n_i+1}, T_2=\sum_{j=1}^{no} H_{0j}\sum_{i=1}^{np} A_{01}(j,i)Q_i, T_3=\sum_{j=1}^{nn}\lambda_j\sum_{i=1}^{np}(A_{21}(j,i)Q_i-q_j) $$

首先计算对 $\lambda_j$ 的偏导，$T_1$ 与 $T_2$ 都不含 $\lambda_j$ 项，

$$ \frac{\partial L}{\partial \lambda_j} = \frac{\partial T_3}{\partial \lambda_j}= \sum_{i=1}^{np}A_{21}(j,i)Q_i-q_j $$

令其为零，得到 $\sum_{i=1}^{np}A_{21}(j,i)Q_i-q_j=0$，即：

$$ \sum_{i=1}^{np}A_{21}(j,i)Q_i=q_j, \: (j=1,\ldots,nn) $$

将这 $nn$ 条方程合在一起，写成矩阵形式就是：$A_{21}Q=q$，这正好就是 Eq.6 的下半部分，也就是连续性方程。

接下来对 $Q_i$ 求导。第一项 $T_1$ 对某个固定的 $Q_i$ 求导，只保留第 $i$ 项：

$$ \frac{\partial T_1}{\partial Q_i} = R_i |Q_i|^{n_i-1}Q_i $$

这一步就是前面一直说的：

$$ \frac{d}{dQ_i} \left( R_i\frac{|Q_i|^{n_i+1}}{n_i+1} \right) = f_i(Q_i) = R_i |Q_i|^{n_i-1}Q_i $$

第二项某个固定的 $Q_i$ 求导：

$$ \frac{\partial T_2}{\partial Q_i} = \sum_{j=1}^{no}H_{0j}A_{01}(j,i) $$

这是与第 $i$ 根管相关的定水头边界贡献。

第三项某个固定的 $Q_i$ 求导：$q_j$ 与 $Q_i$ 无关，所以只对含 $Q_i$ 的部分求导：

$$ \frac{\partial T_3}{\partial Q_i} = \sum_{j=1}^{nn}\lambda_j A_{21}(j,i) $$

最后，将这三项加起来：

$$ \frac{\partial L}{\partial Q_i} = \frac{\partial T_1}{\partial Q_i} + \frac{\partial T_2}{\partial Q_i} + \frac{\partial T_3}{\partial Q_i} = R_i|Q_i|^{n_i-1}Q_i + \sum_{j=1}^{no}H_{0j}A_{01}(j,i) + \sum_{j=1}^{nn}\lambda_jA_{21}(j,i) $$

令其为零，得：

$$ R_i|Q_i|^{n_i-1}Q_i + \sum_{j=1}^{no}H_{0j}A_{01}(j,i) + \sum_{j=1}^{nn}\lambda_jA_{21}(j,i) =0 $$

这就是对每根管 $i=1,\dots,np$ 的一组方程。

接下来就是如何将上面这 $np$ 条标量方程写成矩阵的形式。

  • 首先是第一项 $A_{11}Q$，我们定义了 $A_{11}$, 且可以发现 $A_{11}Q$ 的第 $i$ 个分量正好就是：$R_i|Q_i|^{n_i-1}Q_i$, 这正对应上面偏导方程里的第一项。
  • 然后是第二项 $A_{12}\lambda$，根据 $A_{12}=A_{21}^T$，可得 $A_{12}\lambda$ 的第 $i$ 个分量是：$\sum_{j=1}^{nn}A_{12}(i,j)\lambda_j$, 而因为 $A_{12}(i,j)=A_{21}(j,i)$，所以这就等于：$\sum_{j=1}^{nn}A_{21}(j,i)\lambda_j$, 这正好就是对 $Q_i$ 求导得到的第三项：$\sum_{j=1}^{nn}\lambda_jA_{21}(j,i)$.
  • 第三项，$(A_{10}H_0)_i$ 就是 $\sum_{j=1}^{no}A_{10}(i,j)H_{0j} = \sum_{j=1}^{no}A_{01}(j,i)H_{0j}$.

这个式子可以按上下两部分理解。上半部分：

$$ A_{11}Q + A_{12}\lambda = -A_{10}H_{0} $$

其中，矩阵 $A_{11}$ 定义为：

$$ A_{11}= \begin{bmatrix} R_1|Q_1|^{n_1-1} & & \\ & R_2|Q_2|^{n_2-1} & \\ & & \ddots \\ & & & R_{np}|Q_{np}|^{n_{np}-1} \end{bmatrix} $$

其每一个对角元都对应一根管段的当前等效损失系数：$R_i|Q_i|^{n_i-1}$. 因此：

矩阵 $A_{11}$ 可以理解为全网各管段在当前流量状态下的等效损失系数矩阵。

之所以是对角矩阵，是因为每一根管段的损失函数只依赖于自身流量 $Q_i$，并不直接依赖于其他管段的流量。

同时注意到：

$$ A_{11}Q = \begin{bmatrix} R_1|Q_1|^{n_1-1}Q_1 \\ R_2|Q_2|^{n_2-1}Q_2 \\ \vdots\\ R_{np}|Q_{np}|^{n_{np}-1}Q_{np} \end{bmatrix} = F(Q) $$

这说明：

矩阵 $A_{11}Q$ 恰好就是原始方程中的管段损失向量 $F(Q)$。

回到 Eq.6 中的上半部分，它来源于：

$$ \frac{\partial L}{\partial Q}=0 $$

表示管段上的损失项、节点乘子项与定水头边界项之间的平衡关系，对应于后续识别出的能量方程。

下半部分：

$$ A_{21}Q = q $$

它来源于：

$$ \frac{\partial L}{\partial \lambda}=0 $$

表示未知水头节点上的连续方程，即质量守恒。

因此，公式(6)已经把：

统一写进了一个块矩阵方程中。

最后，通过将 Eq.6 与原始方程 Eq.1 比较，可以直接赋予拉格朗日乘子 $ \lambda $ 物理意义。原始能量关系为：

$$ A_{12}H + F(Q) = -A_{10}H_{0} $$

而 Eq.6 上半部分为：$A_{11}Q + A_{12}\lambda = -A_{10}H_{0}$，由于已经有：$A_{11}Q = F(Q)$. 于是可得：$\lambda = H$.

也就是说：

拉格朗日乘子 $ \lambda $ 在物理上正对应于未知节点水头 $H$。

这一步非常关键，因为它说明 Eq.4 中引入的乘子并不是抽象的数学量，而是与真实水力物理量直接对应的。

在识别出 $\lambda = H$ 后， Eq.6 可进一步写成：

$$ \begin{bmatrix} A_{11} & A_{12} \\ A_{21} & 0 \end{bmatrix} \begin{bmatrix} Q \\ H \end{bmatrix} = \begin{bmatrix} - A_{10}H_{0} \\ q \end{bmatrix} \qquad (8) $$

所以，Eq.8 与最初的原始水力平衡方程组完全等价。

这说明：

从 Content Model 出发构造的优化表达、拉格朗日函数以及一阶条件系统，并不是另一套独立模型，而是原始稳态管网水力平衡问题的等价数学表达。

2.5 牛顿法（Newton–Raphson 迭代）的直接应用

在前面的内容中，已经通过拉格朗日函数的一阶条件，将原始管网稳态水力平衡问题重写为统一的矩阵形式 Eq.8.

但该系统本质上仍是非线性的，原因在于矩阵 $A_{11}$ 的对角元依赖于当前流量 $Q$：$A_{11}=\mathrm{diag}\left(R_i|Q_i|^{n_i-1}\right)$.

因此，本部分的任务是：

在统一方程 Eq.8 的基础上，引入 Newton–Raphson 方法，对非线性系统进行线性化，并建立每一步迭代的修正方程。

为什么要使用 Newton–Raphson 方法

由于矩阵 $A_{11}$ 随流量 $Q$ 而变化，公式 Eq.8 不能被看成普通的线性方程组，因此无法直接一次求解。

为求得其解，论文采用 Newton–Raphson 迭代思想，即：

一个关于数值稳定性的小问题：为什么要避免 $A_{11}$ 奇异：当某根管道两端水头相同、从而流量为零时，矩阵 $A_{11}$ 的某些对角元可能为零，从而导致矩阵奇异或数值不稳定。为了避免这一问题，作者指出：可以为 $A_{11}$ 的对角元设置一个下界，使其不至于变成零。

Newton–Raphson 修正方程

根据上面的守恒矩阵方程 Eq (8)，另

原始能量方程可写为：

$$ A_{11}Q + A_{12}H + A_{10}H_0 = 0 $$

如果说在当前迭代点 $(Q^k,H^k)$ 处，该条件尚未满足，则左侧会出现一个偏差：

$$ dE = A_{11}Q^k + A_{12}H^k + A_{10}H_0 $$

因此，$dE$ 表示：当前迭代点下能量方程的残差。

另外，根据原始连续方程为：$A_{21}Q = q$, 若当前流量 $Q^k$ 还未满足该条件，则有残差：

$$ dq = A_{21}Q^k - q $$

因此，$dq$ 表示：当前迭代点下连续方程的残差。

$dE$ 和 $dq$ 并不是新的物理变量，而是当前猜测解距离真实解的偏差量。在 Newton–Raphson 迭代中，其目标就是：通过不断求解修正量 $dQ$ 和 $dH$，使 $dE$ 和 $dq$ 逐步减小到零。

上面得到的 Eq (8) 为：

$$ \begin{bmatrix} A_{11} & A_{12} \\ A_{21} & 0 \end{bmatrix} \begin{bmatrix} Q \\ H \end{bmatrix} = \begin{bmatrix} -A_{10} H_0 \\ q \end{bmatrix} $$

但这里有一个微妙之处：$A_{11}$ 本身依赖 $Q$，所以这实际上是一个非线性系统：

$$ F(Q, H) = \begin{bmatrix} A_{11}(Q) Q + A_{12} H + A_{10} H_0 \\[2mm] A_{21} Q - q \end{bmatrix} = 0 $$

更清晰地，记：

$$ E(Q, H) = A_{11}(Q) Q + A_{12} H + A_{10} H_0 \quad (\text{能量方程残差}) $$

$$ M(Q) = A_{21} Q - q \quad (\text{连续性残差}) $$

我们要解 $E=0$ 和 $M=0$。

牛顿法用于求解 $F(x)=0$ 时，在当前点 $x^k$ 处做一阶泰勒展开：

$$ F(x^k + \Delta x) \approx F(x^k) + J(x^k) \Delta x $$

设我们希望 $F(x^k + \Delta x) = 0$，则：

$$ J(x^k) \Delta x = -F(x^k) $$

其中 $J(x^k)$ 是雅可比矩阵（偏导数矩阵）。

应用到我们的 $F(Q, H)$ 中，令：

$$ x = \begin{bmatrix} Q \\ H \end{bmatrix}, \quad \Delta x = \begin{bmatrix} dQ \\ dH \end{bmatrix}, \quad F = \begin{bmatrix} E(Q,H) \\ M(Q) \end{bmatrix} $$

然后计算 $J(x^k)$ 的分块。首先是 $E$ 关于 $Q$ 的偏导：

$$ \frac{\partial E}{\partial Q} = \frac{\partial}{\partial Q} \left( A_{11}(Q) Q \right) + \frac{\partial}{\partial Q} (A_{12} H) + 0 $$

关键是第一项。已知 $A_{11}(Q)$ 是对角矩阵，对角元为 $R_i |Q_i|^{n_i-1}$。 所以 $A_{11}(Q) Q$ 的第 $i$ 个分量为 $R_i |Q_i|^{n_i-1} Q_i = f_i(Q_i)$。 因此：

$$ \frac{\partial f_i}{\partial Q_i} = R_i \cdot n_i |Q_i|^{n_i-1} = n_i \cdot (A_{11})_{ii} $$

于是：

$$ \frac{\partial E}{\partial Q} = \operatorname{diag}(n_i) \cdot A_{11}(Q) = N A_{11}(Q) $$

于是雅可比矩阵为：

$$ J(x^k) = \begin{bmatrix} N A_{11}(Q^k) & A_{12} \\ A_{21} & 0 \end{bmatrix} $$

其中， 矩阵 $N$ 的含义为：

$$ N= \begin{bmatrix} n_1 & & \\ & n_2 & \\ & & \ddots \\ & & & n_{np} \end{bmatrix} $$

其为一个 $(np,np) $ 的对角矩阵，其中每个对角元为对应管段损失函数中的指数 $n_i$。其来源于在 Newton 线性化过程中，需要对单管损失函数 $f_i(Q_i)=R_i|Q_i|^{n_i-1}Q_i$ 关于 $Q_i$ 求导。其导数中会自然带出一个系数 $n_i$，因此在矩阵形式下就表现为 $N A_{11}$。

右端项 $-F(x^k)$ 为：

$$ -F(x^k) = -\begin{bmatrix} A_{11}(Q^k) Q^k + A_{12} H^k + A_{10} H_0 \\ A_{21} Q^k - q \end{bmatrix} = \begin{bmatrix} - (A_{11} Q^k + A_{12} H^k + A_{10} H_0) \\ - (A_{21} Q^k - q) \end{bmatrix} $$

为了与论文符号一致，通常定义：

$$ dE = A_{11} Q^k + A_{12} H^k + A_{10} H_0 \quad (\text{能量残差}) $$

$$ dq = A_{21} Q^k - q \quad (\text{连续性残差}) $$

则：

$$ -F(x^k) = \begin{bmatrix} -dE \\ -dq \end{bmatrix} $$

牛顿修正方程为：

$$ J(x^k) \begin{bmatrix} dQ \\ dH \end{bmatrix} = -F(x^k) $$

代入得：

$$ \begin{bmatrix} N A_{11} & A_{12} \\ A_{21} & 0 \end{bmatrix} \begin{bmatrix} dQ \\ dH \end{bmatrix} = \begin{bmatrix} -dE \\ -dq \end{bmatrix} $$

但是论文中的式 (9) 右端是 $\begin{bmatrix} dE \\ dq \end{bmatrix}$，而不是负号。 这是符号习惯的差异：有些作者把修正量定义为 $x^{k+1} = x^k + \Delta x$，则 $J \Delta x = -F$； 论文中可能隐含地将修正量定义为 $dQ = Q^{k+1} - Q^k$ 或 $dQ = Q^k - Q^{k+1}$？检查论文上下文，他们后来写 $dH = H^k - H^{k+1}$，所以修正量实际是负增量，从而右端变为 $+dE, +dq$。这种符号选择不影响最终迭代公式的数值结果，只是形式不同。为了与论文原文一致，我们就接受式 (9) 的形式，并记住它在数值实现时等价于标准的牛顿步。

该式表示：在当前迭代点处，通过求解一个线性修正系统，来确定如何调整 $Q$ 和 $H$，从而逐步逼近原非线性系统的真实解。

2.6 引入 $D$ 简化矩阵

定义：

$$ N A_{11} = D^{-1} \quad \Rightarrow \quad D = (N A_{11})^{-1} $$

因为 $N$ 和 $A_{11}$ 都是对角矩阵，所以 $D$ 也是对角矩阵，容易计算。

于是系统矩阵 Eq (9) 更新为：

$$ \begin{bmatrix} D^{-1} & A_{12} \\ A_{21} & 0 \end{bmatrix} \begin{bmatrix} dQ \\ dH \end{bmatrix} = \begin{bmatrix} dE \\ dq \end{bmatrix} $$

理论上，我们可以直接求解这个系统得到 $dQ$ 和 $dH$。但矩阵大小为 $(np+nn)\times(np+np)$，对于大型管网（np 可达数千甚至数万），直接求解成本很高。

2.7 分块矩阵求逆

注意到 $dQ$ 和 $dH$ 之间通过 $A_{21}dQ = dq$ 以及 $D^{-1}dQ + A_{12}dH = dE$ 耦合。如果能先解出 $dH$，再回代求 $dQ$，就能把问题规模从 $np+nn$ 降到 $nn$。这等价于消去 $dQ$，只保留 $dH$ 的方程。而消去过程可以通过分块求逆系统完成。

分块矩阵求逆的一般公式，对于形如：

$$ M = \begin{bmatrix} P & Q \\ R & 0 \end{bmatrix} $$

的块矩阵，若 $P$ 和 $R P^{-1} Q$ 可逆，则其逆矩阵为：

$$ M^{-1} = \begin{bmatrix} P^{-1} - P^{-1}Q (R P^{-1} Q)^{-1} R P^{-1} & P^{-1}Q (R P^{-1} Q)^{-1} \\ (R P^{-1} Q)^{-1} R P^{-1} & - (R P^{-1} Q)^{-1} \end{bmatrix} $$

这个公式可以通过块消元推导。

将上述分块矩阵求逆的技巧应用到我们的矩阵中，这里 $P = D^{-1}$（对角，可逆），$Q = A_{12}$，$R = A_{21}$，右下角为 $0$。

首先计算 Schur complement：

$$ S = - R P^{-1} Q = - A_{21} D A_{12} $$

因为 $P^{-1} = D$，所以：

$$ S = - A_{21} D A_{12} $$

注意 $A_{21} D A_{12}$ 是对称正定的（论文中称为 Stieltjes 矩阵），可逆。

代入分块求逆公式得到：

$$ M^{-1} = \begin{bmatrix} D - D A_{12} (A_{21} D A_{12})^{-1} A_{21} D & D A_{12} (A_{21} D A_{12})^{-1} \\ (A_{21} D A_{12})^{-1} A_{21} D & - (A_{21} D A_{12})^{-1} \end{bmatrix} = \begin{bmatrix} B_{11} & B_{12} \\ B_{21} & B_{22} \end{bmatrix} $$

其中每个子块的尺寸：

那么：

$$ \begin{bmatrix} dQ \\ dH \end{bmatrix} = \begin{bmatrix} B_{11} & B_{12} \\ B_{21} & B_{22} \end{bmatrix} \begin{bmatrix} dE \\ dq \end{bmatrix} = \begin{bmatrix} B_{11} dE + B_{12} dq \\ B_{21} dE + B_{22} dq \end{bmatrix} $$

2.8 消去 $dQ$，得到 $dH$ 的表达式

由上可得：

$$ dH = B_{21} dE + B_{22} dq $$

代入 $B_{21}, B_{22}$ 表达式：

$$ dH = (A_{21} D A_{12})^{-1} A_{21} D \, dE \;-\; (A_{21} D A_{12})^{-1} dq $$

再代入 $dE$ 和 $dq$ 的定义：

$$ dE = A_{11} Q^k + A_{12} H^k + A_{10} H_0, \quad dq = A_{21} Q^k - q $$

经过整理（合并同类项）就得到论文中的 Eq (15):

$$ dH = (A_{21} D A_{12})^{-1} \bigl[ A_{21} D (A_{11} Q^k + A_{10} H_0) + (q - A_{21} Q^k) \bigr] $$

注意 $dH = H^k - H^{k+1}$，所以 $H^{k+1} = H^k - dH$。

2.9 最终的牛顿迭代格式

根据 Eq (9)：

$$ \begin{cases} N A_{11} dQ + A_{12} dH = dE &\\ A_{21} dQ = dq \end{cases} $$

其中：

$$ dE = A_{11} Q^k + A_{12} H^k + A_{10} H_0, \quad dq = A_{21} Q^k - q $$

且 $dQ = Q^{k+1} - Q^k$，$dH = H^{k+1} - H^k$（注意：论文后来定义为 $dH = H^k - H^{k+1}$，符号不同，我们先用标准定义，最后调整）。

对等式：

$$ N A_{11} dQ = dE - A_{12} dH $$

左乘 $(N A_{11})^{-1} = D$：

$$ dQ = D (dE - A_{12} dH) $$

然后将其带入等式：$A_{21} dQ = dq$ 得：

$$ A_{21} D (dE - A_{12} dH) = dq $$

$$ A_{21} D dE - A_{21} D A_{12} dH = dq $$

然后代入 $dE$ 和 $dq$ 的定义：

$$ A_{21} D (A_{11} Q^k + A_{12} H^k + A_{10} H_0) - A_{21} D A_{12} dH = A_{21} Q^k - q $$

注意到 $D A_{11} = (N A_{11})^{-1} A_{11} = N^{-1}$，所以：

$$ A_{21} N^{-1} Q^k + A_{21} D A_{12} H^k + A_{21} D A_{10} H_0 - A_{21} D A_{12} dH = A_{21} Q^k - q $$

然后将含 $H^k$ 和 $dH$ 的项移到一边

$$ A_{21} D A_{12} H^k - A_{21} D A_{12} dH = A_{21} Q^k - q - A_{21} N^{-1} Q^k - A_{21} D A_{10} H_0 $$

左边提取 $A_{21} D A_{12}$：

$$ A_{21} D A_{12} (H^k - dH) = A_{21} Q^k - q - A_{21} N^{-1} Q^k - A_{21} D A_{10} H_0 $$

若定义 $dH = H^k - H^{k+1}$，则 $H^k - dH = H^{k+1}$。代入：

$$ A_{21} D A_{12} H^{k+1} = A_{21} Q^k - q - A_{21} N^{-1} Q^k - A_{21} D A_{10} H_0 $$

整理右端项：

$$ A_{21} D A_{12} H^{k+1} = - \left[ A_{21} N^{-1} Q^k + A_{21} D A_{10} H_0 + q - A_{21} Q^k \right] $$

然后左乘 $(A_{21} D A_{12})^{-1}$，即可解出 $H^{k+1}$：

$$ H^{k+1} = - (A_{21} D A_{12})^{-1} \left[ A_{21} N^{-1} Q^k + A_{21} D A_{10} H_0 + q - A_{21} Q^k \right] $$

只是右边多了一个负号？检查：你的第一个公式右边没有负号，这里推导出来有负号。差别来自 $dH$ 的定义符号。我们后面统一。

然后继续将 $D = A_{11}^{-1} N^{-1}$ 代入矩阵部分：

$$ A_{21} D A_{12} = A_{21} A_{11}^{-1} N^{-1} A_{12} $$

其逆为：

$$ (A_{21} D A_{12})^{-1} = (A_{21} A_{11}^{-1} N^{-1} A_{12})^{-1} $$

由于 $A_{11}^{-1}$ 和 $N^{-1}$ 都是对角阵且可交换，可写成 $(A_{21} N^{-1} A_{11}^{-1} A_{12})^{-1}$。另外，右端项中的 $A_{21} D A_{10} H_0 = A_{21} A_{11}^{-1} N^{-1} A_{10} H_0$，而论文中写的是 $A_{21} N^{-1} (A_{11}^{-1} A_{10} H_0)$，两者相同。

最终得到论文中的 Eq (18)：

$$ H^{k+1} = - (A_{21} N^{-1} A_{11}^{-1} A_{12})^{-1} \left[ A_{21} N^{-1} (Q^k + A_{11}^{-1} A_{10} H_0) + (q - A_{21} Q^k) \right] $$

接下来开始推导 $Q^{k+1}$，还是从：

$$ N A_{11} dQ + A_{12} dH = dE $$

开始。首先，将 $dQ$ 和 $dH$ 的表达式代入：

$$ N A_{11} (Q^k - Q^{k+1}) + A_{12} (H^k - H^{k+1}) = dE $$

将已知的 $dE = A_{11} Q^k + A_{12} H^k + A_{10} H_0$，代入上式：

$$ N A_{11} (Q^k - Q^{k+1}) + A_{12} (H^k - H^{k+1}) = A_{11} Q^k + A_{12} H^k + A_{10} H_0 $$

展开左边：

$$ N A_{11} Q^k - N A_{11} Q^{k+1} + A_{12} H^k - A_{12} H^{k+1} = A_{11} Q^k + A_{12} H^k + A_{10} H_0 $$

移项，把含 $Q^{k+1}, H^{k+1}$ 的项留在左边，其余移到右边：

$$ - N A_{11} Q^{k+1} - A_{12} H^{k+1} = A_{11} Q^k + A_{12} H^k + A_{10} H_0 - N A_{11} Q^k - A_{12} H^k $$

右边化简：$A_{12} H^k - A_{12} H^k = 0$，于是：

$$ - N A_{11} Q^{k+1} - A_{12} H^{k+1} = A_{11} Q^k - N A_{11} Q^k + A_{10} H_0 $$

$$ - N A_{11} Q^{k+1} - A_{12} H^{k+1} = (I - N) A_{11} Q^k + A_{10} H_0 $$

两边乘以 $-1$：

$$ N A_{11} Q^{k+1} + A_{12} H^{k+1} = (N - I) A_{11} Q^k - A_{10} H_0 $$

$$ N A_{11} Q^{k+1} = (N - I) A_{11} Q^k - A_{10} H_0 - A_{12} H^{k+1} $$

左乘 $(N A_{11})^{-1} = A_{11}^{-1} N^{-1}$（对角矩阵可交换）：

$$ Q^{k+1} = A_{11}^{-1} N^{-1} (N - I) A_{11} Q^k - A_{11}^{-1} N^{-1} (A_{10} H_0 + A_{12} H^{k+1}) $$

注意 $A_{11}^{-1} N^{-1} (N - I) A_{11} = A_{11}^{-1} (N^{-1} N - N^{-1}) A_{11} = A_{11}^{-1} (I - N^{-1}) A_{11}$。
由于 $A_{11}$ 是对角矩阵，与 $I - N^{-1}$ 可交换，所以：

$$ A_{11}^{-1} (I - N^{-1}) A_{11} = I - N^{-1} $$

因此：

$$ Q^{k+1} = (I - N^{-1}) Q^k - N^{-1} A_{11}^{-1} (A_{12} H^{k+1} + A_{10} H_0) $$

3 A Unifying Assessment of the Loop Gradient Procedures

3.1 环关联矩阵的定义

在管网分析中，独立环的数量为：

$$ n\ell = np - nn \quad (\text{当 } no=1 \text{ 时，否则需加 dummy 管道}) $$

以一个田字格（3x3 节点，4 个小方格）为例：

1---2---3
|   |   |
4---5---6
|   |   |
7---8---9
  • 总节点数 =9，管道数 =12
  • 小环有 4 个，还有 1 个大环
  • 但是 独立环数 = 管道数 - 节点数 + 连通子图数 = 12 - 9 + 1 = 4。正好等于小方格的数量。也就是说，这四个小方格环是独立环，所有其他环（如外围大环）都可以由这四个小环组合得到（大环 = 四个小环的对称差）。

结论：

  • 环：任意一个闭合回路。
  • 独立环：一组最少的环，使得所有其他环都可以由它们组合得到。其数量由公式决定。

环关联矩阵 $ \mathbf{M13} $ 的大小为 $ np \times n\ell $，其元素定义如下：

$$ \mathbf{M13}(i,k) = \begin{cases} +1 & \text{若管道 } i \text{ 属于环 } k \text{，且参考方向与环方向一致} \\ -1 & \text{若管道 } i \text{ 属于环 } k \text{，且参考方向与环方向相反} \\ 0 & \text{若管道 } i \text{ 不属于环 } k \end{cases} $$

同时定义 $ \mathbf{M31} = \mathbf{M13}^{\mathrm{T}} $。

3.2 拓扑性质

论文给出两个关键性质（式 20）：

$$ \mathbf{M31} \mathbf{A12} = \mathbf{0}, \quad \mathbf{A21} \mathbf{M13} = \mathbf{0} $$

$$ \mathbf{M31} \mathbf{A10} = \mathbf{0}, \quad \mathbf{A01} \mathbf{M13} = \mathbf{0} $$

这些性质来源于其物理意义。 乘积 $\mathbf{M31} \mathbf{A12}$ 的第 $(k,j)$ 个元素为：

$$ \sum_{i=1}^{np} M31(k,i) \cdot A12(i,j) $$

即：环 $k$ 中所有管道对节点 $j$ 的“带符号贡献”之和。这个说法其实对应的是另一个常用的环性质：

$$ \sum_{\text{pipe in loop}} \text{head loss} = 0 $$

而水头损失可以写成 $H_{\text{from}} - H_{\text{to}}$。当我们将这个求和写成矩阵形式时，确实会用到 $\mathbf{M31} \mathbf{A12}$，但那是乘以水头向量 $H$ 之后才出现的：

$$ \mathbf{M31} (\mathbf{A12} H) = 0 $$

因为 $\mathbf{A12} H$ 给出了每个管道两端的水头差（有符号），然后环矩阵把这些水头差加起来，正好是环一周的总水头变化，当然为零。 但 $\mathbf{M31} \mathbf{A12} = \mathbf{0}$ 本身并不需要 $H$，它是在说：即使你把 $H$ 替换成任意向量，只要 $H$ 是节点上的值，环上水头差之和总为零。这其实是同一个拓扑性质的表现。

3.3 线性理论法

前面第二章中，已经详细讲解了节点梯度法。接下来，我们会继续推导另外两种方法：

方法操作求解规模特点
节点梯度法直接解 $dH$（第二章）$nn \times nn$对称正定，任意初值
线性理论法左乘 $\begin{bmatrix} M_{31} & 0 \\ 0 & I \end{bmatrix}$，消去 $dH$$np \times np$非对称，规模最大
环梯度法再令 $dQ = M_{13} dQ_\ell$，消去 $dH$ 和节点流量$n\ell \times n\ell$对称正定，但需初值满足连续性

原始修正方程 Eq (9)（已用 $D$ 简化）：

$$ \begin{bmatrix} D^{-1} & A12 \\ A21 & 0 \end{bmatrix} \begin{bmatrix} dQ \\ dH \end{bmatrix} = \begin{bmatrix} dE \\ dq \end{bmatrix} $$

其中 $dE, dq$ 见 Eq (10)。现在，构造一个变换矩阵（大小为 $np \times (np+nn)$）：

$$ \begin{bmatrix} \mathbf{M31} & \mathbf{0} \\ \mathbf{0} & \mathbf{I}_{nn} \end{bmatrix} $$

左乘该矩阵于修正方程两边，得到：

$$ \begin{bmatrix} \mathbf{M31} D^{-1} & \mathbf{M31} A12 \\ A21 & 0 \end{bmatrix} \begin{bmatrix} dQ \\ dH \end{bmatrix} = \begin{bmatrix} \mathbf{M31} dE \\ dq \end{bmatrix} $$

利用性质 $\mathbf{M31} A12 = 0$，左上角变为 $\mathbf{M31} D^{-1}$，左下角仍为 $A21$，右上角为 0。于是方程组变成：

$$ \begin{cases} \mathbf{M31} D^{-1} dQ = \mathbf{M31} dE \\ A21 dQ = dq \end{cases} $$

这是一个关于 $dQ$ 的线性系统（共 $np$ 个方程），不再含有 $dH$。这就是线性理论法的核心：直接求解管道流量的修正量 $dQ$，不需要先求水头。 写成矩阵形式：

$$ \begin{bmatrix} \mathbf{M31} D^{-1} \\ A21 \end{bmatrix} dQ = \begin{bmatrix} \mathbf{M31} dE \\ dq \end{bmatrix} $$

通常需要求解一个 $np \times np$ 的线性系统（虽然不对称），这就是论文中 Eq (22) 所表达的内容。

3.4 环梯度法

如果我们进一步利用环的流量修正量 $dQ_\ell$（大小为 $n\ell$）与管道流量修正量的关系：

$$ dQ = \mathbf{M13} \, dQ_\ell $$

即认为管道流量的修正量可以表示为环修正量的线性组合（每个环一个修正量，沿环方向叠加）。 将 $dQ = M13 \, dQ_\ell$ 代入线性理论的第一式 $\mathbf{M31} D^{-1} dQ = \mathbf{M31} dE$，得到：

$$ \mathbf{M31} D^{-1} \mathbf{M13} \, dQ_\ell = \mathbf{M31} dE $$

记 $\mathbf{L} = \mathbf{M31} D^{-1} \mathbf{M13}$，这是一个 $n\ell \times n\ell$ 的对称正定矩阵（因为 $D^{-1}$ 正定）。于是：

$$ dQ_\ell = \mathbf{L}^{-1} (\mathbf{M31} dE) $$

由定义：

$$ dE = A_{11} Q^k + A_{12} H^k + A_{10} H_0 $$

左乘 $M_{31}$：

$$ M_{31} dE = M_{31} A_{11} Q^k + M_{31} A_{12} H^k + M_{31} A_{10} H_0 $$

利用拓扑性质：

所以：

$$ M_{31} dE = M_{31} A_{11} Q^k $$

这个化简非常关键：它消去了水头 $H^k$ 和固定水头 $H_0$。因此，环梯度法的每步迭代不需要计算水头，只依赖当前流量。因此，可以得到：

$$ (M_{31} D^{-1} M_{13}) \, dQ_\ell = M_{31} A_{11} Q^k $$

注意 $D^{-1} = N A_{11}$，所以：

$$ M_{31} D^{-1} M_{13} = M_{31} N A_{11} M_{13} $$

因此：

$$ dQ_\ell = \left( M_{31} A_{11} M_{13} \right)^{-1} M_{31} A_{11} Q^k $$

因为 $dQ = Q^k - Q^{k+1} = M_{13} dQ_\ell$，所以：

$$ Q^{k+1} = Q^k - M_{13} \, dQ_\ell = Q^k - M_{13} \left( M_{31} A_{11} M_{13} \right)^{-1} M_{31} A_{11} Q^k $$

这就是论文的Eq (25)，环梯度法的迭代公式为：

$$ Q^{k+1} = Q^k - \mathbf{M13} \left( \mathbf{M31} \mathbf{A11} \mathbf{M13} \right)^{-1} \mathbf{M31} \mathbf{A11} Q^k \tag{25} $$

缺点：

3.5 结论

三种方法每步收敛速度相同（都是牛顿法的二阶收敛），但：

4 Pros and Cons of the Different Formulations

4.1 输入要求 (Input requirements)

注意：虽然自动生成环矩阵在技术上可行，但论文认为相对于节点法，这仍然是额外的复杂度。

4.2 初始解要求 (Initial solution requirements)

4.3 线性系统的规模 (Size of the system of linear equations)

论文比较的是每步迭代中需要求解的线性方程组的矩阵维数：

注意：矩阵规模小不一定求解快，因为还涉及稀疏性和条件数。所以下一项更重要。

4.4 线性系统的求解效率 (Efficient solution)

4.5 综合打分表及结论

三种方法来自同一 Newton–Raphson 修正方程，只是在不同的空间投影。论文给出了一个评分表：

方法INPUTINITIAL SOLUTIONSIZEEFFICIENT SOLUTIONSummary
NODAL1121✅ 输入简单，任意初始解，矩阵对称正定稀疏，求解最快
LINEAR T.2133❌ 输入稍复杂，任意初始解，矩阵非对称规模最大，最慢
LOOP2212❌ 输入复杂（需环矩阵），需初始解满足连续性，矩阵规模最小但较稠密，求解速度中等

结论：节点梯度法总分最低（5分），最适宜用于大型管网，尤其是微机（PC）。

4.6 补充说明：为什么环梯度法矩阵规模小但求解效率不是最好？

因此，虽然环梯度法的矩阵维数小，但实际计算成本并不占优。

5．The ICF/MCG Solution for the System of Linear Equations

5.1 节点梯度法需要求解的线性系统

在节点梯度法的每一步迭代 Eq (18) 中，都需要求解一个形如

$$ A x = b $$

的线性方程组。这里具体对应的是：

论文强调：矩阵 $A$ 具有三个重要性质：

  1. 对称：$A^T = A$
  2. 正定：所有特征值 > 0
  3. Stieltjes 类型：非对角元 ≤ 0，对角元 > 0，且是 M-矩阵的一种

这些性质使得可以使用共轭梯度法（Conjugate Gradient, CG）来迭代求解，并且可以结合不完全 Cholesky 预处理加速收敛。

5.2 为什么不用直接法？

直接法（如高斯消元、Cholesky 分解）虽然理论上精确，但对于稀疏矩阵，分解过程中会产生填充（fill-in）：原来为零的位置变成非零，导致存储量和计算量急剧增加。对于大型管网（例如 $nn=1000$），填充后可能接近稠密矩阵，需要 $O(nn^3)$ 运算，不可接受。

因此，论文采用迭代法，特别是 CG 法，只利用矩阵的非零元进行矩阵-向量乘法，存储量小。

5.3 共轭梯度法（CG）的基本思想

对于线性方程组：

$$ A x = b $$

其中 $A$ 是 对称正定（SPD）矩阵（比如节点梯度法中的 $A_{21} D A_{12}$），$b$ 已知，$x$ 未知（比如 $H^{k+1}$）。

CG 是一种迭代法：从一个初始猜测 $x_0$ 开始，逐步生成序列 $x_1, x_2, \dots$，使得 $x_k$ 快速逼近精确解 $x^* = A^{-1}b$。

核心思想：极小化一个二次函数

对于 SPD 矩阵 $A$，方程 $A x = b$ 等价于极小化下面的二次函数：

$$ \phi(x) = \frac{1}{2} x^T A x - b^T x $$

因为梯度 $\nabla \phi(x) = A x - b$，令其为零即得原方程。

CG 在每一步沿着某个搜索方向 $p_k$ 前进，选择步长 $\alpha_k$ 使得 $\phi(x_k + \alpha p_k)$ 最小。关键技巧是：这些搜索方向是 A-共轭 的，即 $p_i^T A p_j = 0$ 对 $i \ne j$。这个性质保证了每一步的更新不会破坏前面已经获得的最优性，因此理论上最多 $n$ 步（$n$ 是未知数个数）就能到达精确解。

CG 每步只需要计算 $A p_k$

CG 的迭代公式（标准形式）如下所示。首先，初始化：$x_0$，$r_0 = b - A x_0$，$p_0 = r_0$，重复直到收敛：

$$ \begin{aligned} \alpha_k &= \frac{r_k^T r_k}{p_k^T A p_k} \\[2mm] x_{k+1} &= x_k + \alpha_k p_k \\[2mm] r_{k+1} &= r_k - \alpha_k A p_k \\[2mm] \beta_k &= \frac{r_{k+1}^T r_{k+1}}{r_k^T r_k} \\[2mm] p_{k+1} &= r_{k+1} + \beta_k p_k \end{aligned} $$

从公式可见：

完全不需要知道 $A$ 的单个元素，只需要能够给出输入向量 $p$ 对应的输出 $A p$。这就是“矩阵 $A$ 不需要显式存储”的含义。

如何通过“只计算 $A p$”实现收敛？

CG 本质上是通过 Krylov 子空间 方法迭代求解。它生成的解 $x_k$ 满足：

$$ x_k \in x_0 + \mathcal{K}_k(A, r_0) $$

其中 Krylov 子空间 $\mathcal{K}_k(A, r_0) = \text{span}\{r_0, A r_0, A^2 r_0, \dots, A^{k-1} r_0\}$。

但是，CG 不直接计算这些高次幂，而是通过递推产生共轭方向。每一步只要知道 $A p_k$，就能隐式地拓展 Krylov 子空间。当 $k = n$ 时，Krylov 子空间已经包含了所有可能的搜索方向，解必然精确。这就是为什么理论上最多 $n$ 步收敛。

实际收敛速度：为什么通常远少于 $n$ 步？

CG 的收敛速度取决于矩阵 $A$ 的 条件数 $\kappa(A) = \frac{\lambda_{\max}}{\lambda_{\min}}$（最大特征值除以最小特征值）。在 $A$‑范数下，误差满足：

$$ \|x_k - x^*\|_A \le 2 \left( \frac{\sqrt{\kappa} - 1}{\sqrt{\kappa} + 1} \right)^k \|x_0 - x^*\|_A $$

其中 $\|z\|_A = \sqrt{z^T A z}$。

例子：

在管网问题中，$A = A_{21} D A_{12}$ 的条件数通常不会极大（尤其是经过适当标度），因此实际收敛步数远小于 $nn$（比如几十步）。论文中使用的 不完全 Cholesky 预处理 进一步将条件数降低到接近 1，从而只需很少的 CG 迭代（比如 5～10 步）即可达到工程精度。

为什么“只计算 $A p$”能高效收敛？

5.4 预处理（Preconditioning）：为什么需要不完全 Cholesky？

CG 的收敛速度取决于矩阵 $A$ 的条件数 $\kappa(A) = \lambda_{\max}/\lambda_{\min}$。如果条件数很大，CG 需要很多步。

预处理就是找一个可逆矩阵 $P$，使得 $P^{-1} A$ 的条件数更接近 1，从而加速收敛。常用方法是选择一个容易求逆的矩阵 $M \approx A$，然后求解等价系统 $M^{-1} A x = M^{-1} b$（实际实现更精巧）。

不完全 Cholesky 分解（ICF）：对 $A$ 做近似 Cholesky 分解：

$$ A \approx L L^T $$

其中 $L$ 是下三角矩阵，与 $A$ 具有相同的非零模式（即不允许填充）。这样分解得到的 $L$ 比完全 Cholesky 的 $L$ 稀疏得多。然后取 $M = L L^T$ 作为预处理子。

为什么叫“不完全”：因为完全 Cholesky 分解会产生填充（原来为零的位置变成非零），而不完全分解禁止填充，只允许在原来非零的位置上计算，因此分解结果只是 $A$ 的近似。

5.5 ICF/MCG 算法步骤（简略）

论文引用的 D. Kershaw (1978) 算法：

  1. 对 $A$ 计算不完全 Cholesky 分解：$A \approx L L^T$（只计算与 $A$ 相同稀疏模式的元素）。
  2. 使用预处理共轭梯度法（PCG）求解 $A x = b$，每次迭代中需要求解形如 $L L^T z = r$ 的三角方程组（容易求解，因为 $L$ 是下三角稀疏）。
  3. 迭代直到残差足够小。

优点：

当前页面是本站的「Google AMP」版。查看和发表评论请点击：完整版 »