基于梯度算法的供水管网水力计算实现
论文介绍论文题目：A Gradient Algorithm for the Analysis of Pipe Networks作者：E. Todini, S. Pilati (1988)供水管网水力分析（稳态流分析）是管网设计与优化的核心环节。其任务是在已知管网拓扑、管道参数、节点需水量和固定水头（水库、水塔）的条件下，求解所有管道的流量 $Q$ 和所有未知节点的测压水头 $H$。传统方法（如 Hardy‑Cross、牛顿法、线性理论、优化法）各有缺陷：或收敛慢、或依赖初值、或方程组规模过大、或数值实现复杂。作者提出了一种梯度算法（Gradient Algorithm）：在流量‑水头联合空间 $(Q, H)$ 中应用牛顿法，并借助 Content Model（内容模型） 的凸性证明了解的存在唯一性——这是算法无条件收敛的理论基石。通过分块矩阵消元，将每步牛顿迭代压缩为只求解一个与未知节点水头数 $nn$ 同规模的线性系统，再通过简单回代得到新的管道流量。该线性系统具有对称正定、稀疏、Stieltjes 矩阵的良好性质，可采用 不完全 Cholesky 分解 + 共轭梯度法（ICF/MCG） 高效求解，存储量仅与管道数 $np$ 成正比，适合在个人计算机上求解大型管网。论文（核心）结构：第二章：推导节点梯度法的完整迭代公式。第三章：从同一修正方程导出线性理论和环梯度法，展示三种方法的统一性。第四章：从输入、初值、规模、求解效率四个维度比较，证明节点梯度法最优。第五章：介绍 ICF/MCG 算法，实现节点法的高效求解。本文提出的梯度算法后来成为全球广泛使用的供水管网物理模拟软件 EPANET 的核心求解器（称为 Global Gradient Algorithm）。1 Introduction引言其实在讲四件事：管网问题分两步：一是设计与优化，二是水力分析。本文主要研究第二个，即稳态管网水力分析。稳态分析要解什么：在稳态（管网达到稳定运行的）条件下，求：每条管段的流量以及每个节点的测压水头。旧方法有什么问题：作者回顾了四类方法，Hardy Cross 一类的局部梯度法：可能发散，收敛慢牛顿法：更好，但仍依赖初值和网络复杂度线性理论法：更稳，但方程组更大优化法：理论漂亮，但数值实现较重作者的方法新在哪里：作者的方法想把两类优点结合起来，从优化思路那里获得“解存在且唯一”的理论基础从牛顿法那里获得高效迭代求解能力最终再把问题压缩成一个只和未知节点水头数有关的线性系统这就是整篇论文真正的核心出发点。2 The Derivation of the Proposed Newton-Raphson Algorithm2.1 问题定义已知量管网拓扑（管道连接关系）每根管道的特性：$R_i, n_i$（Hazen‑Williams 或 Darcy‑Weisbach 参数）每个节点的需水量 $q_j$（正值表示流出）固定水头节点（水库、水塔）的水头值 $H0$（$no$ 个）未知量每根管道的流量 $Q_i$（$np$ 个）每个未知节点的水头 $H_j$（$nn$ 个）具体的形式定义对于管段：$Q^T=[Q_1,Q_2,\ldots,Q_{np}]$，形状为 $(1,np)$，表示每一根管道中的流量。$F^T(Q)=[f_1,f_2,\ldots,f_{np}]$，形状为 $(1,np)$，用于表示各管道中的水头损失规律，其中每一项 $f_i$ 都是对应流量 $Q_i$ 的函数，即 $f_i=f_i(Q_i)$。若采用 Hazen–Williams 水头损失函数，则对于每一根管道可写为：$f_i(Q_i)=R_{i}\left|Q_{i}\right|^{n_{i}-1} Q_{i}$.对于节点：$q^T=[q_1,q_2,\ldots,q_{nn}]$，形状为 $(1,nn)$，表示各未知节点的需水量。$H^T=[H_1,H_2,\ldots,H_{nn}]$，形状为 $(1,nn)$，表示各未知节点的水头。$H_0^T=[H_{0_{1}},H_{0_{2}},\ldots,H_{0_{no}}]$，形状为 $(1,no)$，表示各定水头节点的水头。管段与节点的连接关系使用 $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 控制方程的矩阵形式质量守恒：作用在节点上，流进来的水和流出去的水必须平衡（节点 $q$ 与管段 $Q$ 的关系）。能量守恒：作用在管段上，管段两端的水头差必须和这条管道上的能量损失相匹配（管段多个属性与节点 $H$ 的关系）。首先来看质量守恒：$$ 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) $$设定：源点是定水头节点，水头已知为 $H_0$Node 1 和 Node 2 是未知水头节点两根管道：Pipe 1：源点到 Node 1，流量 $Q_1$Pipe 2：Node 1 到 Node 2，流量 $Q_2$节点需水：Node 1 需水 $q_1$Node 2 需水 $q_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} $$第一行对应 Node 1：Pipe 1 流入，记 $+1$；Pipe 2 流出，记 $-1$；第二行对应 Node 2：Pipe 2 流入，记 $+1$.注意：这里是 $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 解的存在唯一性证明凸函数与凸优化一个优化问题由两部分组成：目标函数 + 约束条件。凸优化要求同时满足：目标函数是凸函数直观理解：函数图像像一个“碗”（开口向上），任意两点连线都在图像上方。数学定义：$f(\theta x + (1-\theta)y) \le \theta f(x) + (1-\theta) f(y)$，对任意 $x,y$ 和 $\theta \in [0,1]$ 成立。约束条件定义的可行域是凸集线性约束（如 $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 点）。其中：$Q$：全局管段流量向量$np$：未知流量管段数$nn$：未知水头节点数$no$：定水头节点数$R_i$：第 $i$ 根管道的阻力参数$H_{0j}$：第 $j$ 个定水头节点的已知水头$A_{21}$：未知水头节点—管段关联矩阵$A_{01}$：定水头节点—管段关联矩阵$q_j$：第 $j$ 个未知水头节点的需水量上式是对原始管网水力平衡问题的一种等价数学重写。通过构造一个全局目标函数 $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) $$其中：$\lambda=[\lambda_1,\lambda_2,\dots,\lambda_{nn}]^T$：与节点连续方程对应的拉格朗日乘子向量。先来通过一个简单的概念来理解下什么是拉格朗日乘子。假设你想让两个数的乘积 $ 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 可以看成由三部分组成：单管损失积分项：$\sum_{i=1}^{np} R_i \frac{|Q_i|^{n_i+1}}{n_i+1}$，保证后续对 $Q_i$ 求导时，可以恢复第 $i$ 根管段的损失函数 $f_i(Q_i)$；定水头边界项：$\sum_{j=1}^{no} H_{0j}\sum_{i=1}^{np} A_{01}(j,i)Q_i$，保证后续对 $Q_i$ 求导时，定水头边界 $H_0$ 能进入一阶条件；拉格朗日乘子约束项：$\sum_{j=1}^{nn}\lambda_j\sum_{i=1}^{np}(A_{21}(j,i)Q_i-q_j)$，将原来“写在外面”的节点连续约束吸收到统一函数中，并使节点信息能够在后续对 $Q$ 求导时进入一阶条件。但在本文推导中，$ \lambda $ 不是单纯的数学符号。后续通过一阶条件可以发现：$ \lambda $ 与未知节点水头 $H$ 具有对应关系，且最终会被识别为与节点水头同类的物理量。公式 Eq.3 虽然已经把原问题写成了约束优化问题，但如果继续直接处理约束，会比较不方便。引入公式 Eq.4 后，可以统一通过偏导来处理：对 $ \lambda $ 求导，恢复连续方程对 $Q$ 求导，恢复能量关系因此，公式 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$ 求导时，会同时得到三类贡献：来自第一项恢复单管损失函数：$f_i(Q_i)=R_i|Q_i|^{n_i-1}Q_i$来自第二项恢复定水头边界项：$\sum_{j=1}^{no} H_{0j}A_{01}(j,i)$, 这是一个常数（因为 $H0_j$ 已知，$A_{01}(j,i)$ 已知），它代表了固定水头节点通过管道 $i$ 施加的“边界能量贡献”。来自第三项恢复与该管段相连节点对应的乘子项：$\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 $ 分别求偏导并令其为零，建立该约束优化问题的一阶条件系统，并进一步说明该系统与原始管网稳态水力平衡方程组等价。论文中指出，当：各阻力参数 $R_i > 0$各指数 $n_i > 0$时，拉格朗日函数 $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 的物理与数学意义上一小节也进行了详细说明：对 $Q_i$ 求偏导：对应第 $i$ 根管段上的平衡条件，最终可恢复能量关系。对 $\lambda_j$ 求偏导：对应第 $j$ 个未知水头节点上的约束条件，最终可恢复连续方程（质量守恒）。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 迭代思想，即：先给定当前迭代点 $(Q^k,H^k)$计算此时原方程的残差在该点附近对原系统进行线性化解出流量修正量 $dQ$ 与水头修正量 $dH$迭代更新，直到残差趋于零一个关于数值稳定性的小问题：为什么要避免 $A_{11}$ 奇异：当某根管道两端水头相同、从而流量为零时，矩阵 $A_{11}$ 的某些对角元可能为零，从而导致矩阵奇异或数值不稳定。为了避免这一问题，作者指出：可以为 $A_{11}$ 的对角元设置一个下界，使其不至于变成零。Newton–Raphson 修正方程根据上面的守恒矩阵方程 Eq (8)，另$Q^k$：第 $k$ 次迭代时的流量向量$H^k$：第 $k$ 次迭代时的未知节点水头向量原始能量方程可写为：$$ 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) $$对 $E$ 关于 $H$ 的偏导：$\frac{\partial E}{\partial H} = A_{12}$（因为 $A_{12} H$ 对 $H$ 是线性的，且 $A_{11} Q$ 和 $A_{10} H_0$ 不含 $H$）对 $M$ 关于 $Q$ 的偏导：$\frac{\partial M}{\partial Q} = A_{21}$（因为 $A_{21} Q - q$ 对 $Q$ 是线性的）对 $M$ 关于 $H$ 的偏导：$\frac{\partial M}{\partial H} = 0$于是雅可比矩阵为：$$ 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} $$其中每个子块的尺寸：$B_{11}$：$np \times np$$B_{12}$：$np \times nn$$B_{21}$：$nn \times np$$B_{22}$：$nn \times nn$那么：$$ \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 Procedures3.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} A_{12} = 0$$M_{31} A_{10} = 0$（因为固定水头节点不在任何环上，或通过 dummy 管道处理后成立）所以：$$ 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} $$缺点：需要初始流量 $Q^0$ 满足 $A_{21} Q^0 = q$。丢失固定水头信息 $H_0$，解不唯一。3.5 结论三种方法每步收敛速度相同（都是牛顿法的二阶收敛），但：节点法最易用（任意初值，无需环信息）线性理论最笨重（规模大，非对称）环法规模最小，但初值要求苛刻4 Pros and Cons of the Different Formulations4.1 输入要求 (Input requirements)节点梯度法 (Nodal)：只需提供管道表（起止节点、长度、直径等）和节点表（需水量、固定水头标记等）。不需要定义环或环关联矩阵。打分 1（最简）。线性理论 (Linear T.)：需要节点数据和管道数据，额外需要环关联矩阵 $M_{13}$ 或能自动生成它的信息。打分 2。环梯度 (Loop)：同样需要环关联矩阵。此外，当有多个固定水头节点时，还需要添加 dummy 管道（人工虚拟管道）使公式 $n\ell = np - nn$ 成立，这增加了输入工作量。打分 2。注意：虽然自动生成环矩阵在技术上可行，但论文认为相对于节点法，这仍然是额外的复杂度。4.2 初始解要求 (Initial solution requirements)节点梯度法：可以从任意初始流量开始（例如全零），算法能收敛。打分 1。线性理论：同样任意初始流量，因为它的方程中已经包含了连续性残差 $dq$，会自动修正。打分 1。环梯度：必须提供一个初始流量 $Q^0$ 满足节点连续性方程 $A_{21} Q^0 = q$。这通常需要预先求解一个线性系统（论文式29），多了一步准备工作。打分 2。4.3 线性系统的规模 (Size of the system of linear equations)论文比较的是每步迭代中需要求解的线性方程组的矩阵维数：节点梯度：矩阵为 $A_{21} D A_{12}$，大小 $nn \times nn$（未知水头节点数）。打分 2（中等）。线性理论：矩阵为 $\begin{bmatrix} M_{31} D^{-1} \\ A_{21} \end{bmatrix}$ 之类（实际推导中需要求解一个 $np \times np$ 系统，非对称），规模 $np$（管道数）。在管网中 $np$ 通常远大于 $nn$，所以打分 3（最差）。环梯度：矩阵为 $M_{31} A_{11} M_{13}$，大小 $n\ell \times n\ell$，$n\ell = np - nn$（对于单水源）。通常 $n\ell$ 小于 $nn$ 也远小于 $np$，所以打分 1（最优）。注意：矩阵规模小不一定求解快，因为还涉及稀疏性和条件数。所以下一项更重要。4.4 线性系统的求解效率 (Efficient solution)节点梯度：矩阵 $A_{21} D A_{12}$ 对称正定，稀疏，非零元素约 $np + nn$（每个管道贡献两个非零元，加上对角）。可使用共轭梯度法（CG）或不完全 Cholesky 预处理（ICF/MCG），迭代次数少，存储低。打分 1。线性理论：矩阵非对称，难以用 CG，一般需直接高斯消去或非对称迭代法，填充严重，运算量约 $O(np^2)$ 甚至更高。打分 3。环梯度：矩阵 $M_{31} A_{11} M_{13}$ 对称正定，但通常比节点矩阵更稠密（因为环与环共享管道导致填充）。尽管维数小，但实际迭代求解时可能需更多次迭代或更多存储。论文引用多个作者的数值实验表明：节点法实际求解时间通常更短。打分 2。4.5 综合打分表及结论三种方法来自同一 Newton–Raphson 修正方程，只是在不同的空间投影。论文给出了一个评分表：方法INPUTINITIAL SOLUTIONSIZEEFFICIENT SOLUTIONSummaryNODAL1121✅ 输入简单，任意初始解，矩阵对称正定稀疏，求解最快LINEAR T.2133❌ 输入稍复杂，任意初始解，矩阵非对称规模最大，最慢LOOP2212❌ 输入复杂（需环矩阵），需初始解满足连续性，矩阵规模最小但较稠密，求解速度中等结论：节点梯度法总分最低（5分），最适宜用于大型管网，尤其是微机（PC）。4.6 补充说明：为什么环梯度法矩阵规模小但求解效率不是最好？填充：环矩阵 $M_{31} A_{11} M_{13}$ 中，每个环与环之间如果有共享管道，就会产生非零元。对于实际管网，环往往相互重叠，导致该矩阵的密度远大于节点矩阵。节点矩阵的每一行对应一个节点，非零元只来自与该节点相连的管道（通常度数不高），因此非常稀疏。迭代次数：CG 的收敛速度取决于矩阵的条件数。虽然两种矩阵都是对称正定，但环矩阵的条件数可能更差（尤其当环的阻抗差异大时），导致需要更多迭代。直接法：如果使用直接法（如 Cholesky 分解），节点矩阵的填充也较少，而环矩阵的填充更多，计算量更大。因此，虽然环梯度法的矩阵维数小，但实际计算成本并不占优。5．The ICF/MCG Solution for the System of Linear Equations5.1 节点梯度法需要求解的线性系统在节点梯度法的每一步迭代 Eq (18) 中，都需要求解一个形如$$ A x = b $$的线性方程组。这里具体对应的是：$A = A_{21} D A_{12} = A_{21} N^{-1} A_{11}^{-1} A_{12}$，对称正定矩阵，大小 $nn \times nn$（未知水头节点数）$x = H^{k+1}$，未知节点水头向量$b = - \left[ A_{21} N^{-1} (Q^k + A_{11}^{-1} A_{10} H_0) + (q - A_{21} Q^k) \right]$，右端项论文强调：矩阵 $A$ 具有三个重要性质：对称：$A^T = A$正定：所有特征值 > 0Stieltjes 类型：非对角元 ≤ 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_k$（矩阵‑向量乘积）需要向量内积（$r_k^T r_k$，$p_k^T A p_k$ 等）需要简单的向量加法与数乘完全不需要知道 $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}$。例子：如果 $\kappa = 100$，则 $\frac{\sqrt{100}-1}{\sqrt{100}+1} = \frac{10-1}{10+1} = \frac{9}{11} \approx 0.818$。每步误差减少约 18%。要减少到 1e-6，大约需要 $k \approx \frac{\ln(1e-6)}{\ln(0.818)} \approx \frac{-13.8}{-0.20} \approx 69$ 步。如果 $\kappa = 10^6$，则 $\frac{1000-1}{1000+1} \approx 0.998$，每步只减少 0.2%，需要上千步。如果 $\kappa = 1$，一步就收敛（直接解）。在管网问题中，$A = A_{21} D A_{12}$ 的条件数通常不会极大（尤其是经过适当标度），因此实际收敛步数远小于 $nn$（比如几十步）。论文中使用的 不完全 Cholesky 预处理 进一步将条件数降低到接近 1，从而只需很少的 CG 迭代（比如 5～10 步）即可达到工程精度。为什么“只计算 $A p$”能高效收敛？不存储 $A$：内存需求从 $O(nn^2)$ 降到 $O(np)$（因为只需要存储管道拓扑和 $D$ 对角元）。计算 $A p$ 快速：利用稀疏性，计算 $A p$ 只需遍历所有管道一次，耗时 $O(np)$，而不是 $O(nn^2)$。迭代次数少：即使每次迭代 $O(np)$，乘以几十步后总工作量远小于直接法（比如 Cholesky 分解的 $O(nn^3)$）。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) 算法：对 $A$ 计算不完全 Cholesky 分解：$A \approx L L^T$（只计算与 $A$ 相同稀疏模式的元素）。使用预处理共轭梯度法（PCG）求解 $A x = b$，每次迭代中需要求解形如 $L L^T z = r$ 的三角方程组（容易求解，因为 $L$ 是下三角稀疏）。迭代直到残差足够小。优点：存储量小：只需存储 $A$ 的非零元（约 $np + nn$ 个）和 $L$ 的非零元（与 $A$ 相同数量）。计算量小：每次 CG 迭代只需 $O(np)$ 次操作（因为矩阵-向量乘利用了稀疏性）。收敛快：预处理后条件数大大改善，通常 10～20 次迭代即可达到工程精度。
❓ 需求分析首先展示下我的问题，以便分析本篇博客中的哪些步骤适用于你的需要。已知条件：有某地区管网的.inp文件，并且有该地区水源节点历史两个月的余氯时间序列数据。求解结果：得到该管网结构中所有节点在历史两个月中的余氯值，更准确的结果如下表展示所示：Time节点 01节点 02……节点 m时刻 012.564.28……2.56时刻 023.564.56……2.56…………………………时刻 n3.443.89……2.56🫡 基础知识什么是WNTR？简单来说：WNTR（Water Network Tool for Resilience）是由美国国家标准与技术研究院（NIST）开发的一套用于分析水分配系统的开源工具，用于模拟和评估给水管网系统的性能和恢复力。主要作用包括：模拟水分配系统、恢复力分析、情景分析，以及优化调度等。简单来说，WNTR就是基于Python的高级版EPANET。EPANET中的Patterns是干什么的？在EPANET中，"Patterns"（模式）是用来表示随时间变化的水需求、压力或其他参数的变化模式。模式是一个数组或一系列因子，这些因子用于调整基础需求或其他输入参数，以模拟系统在不同时间段的行为。应用场景：水需求模式：模式通常用于表示水需求的日变化。例如，居民区的用水需求在早晨和傍晚会有高峰，而在夜间会较低。通过设置不同时间段的需求因子，可以精确模拟真实的用水情况。泵操作模式：模式还可以用于控制泵的操作。例如，可以设置泵在特定时间段的运行速度或开启/关闭状态，以模拟不同的运行策略。压力调节模式：在某些情况下，模式也可以用于调节系统中的压力，模拟压力管理策略。举个例子：假设有一个居民区的用水模式，表示一天24小时内的需求变化：时间步长（小时）: 1 模式因子: [0.5, 0.6, 0.7, 0.8, 1.0, 1.2, 1.5, 1.8, 2.0, 1.8, 1.6, 1.4, 1.2, 1.0, 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.3, 0.4, 0.5]在这个模式中，0点到1点的需求因子是0.5，表示这段时间的需求是基础需求的50%。早晨6点到9点的需求因子是1.2到2.0，表示用水需求较高，为基础需求的120%到200%。🕹️ 具体步骤简要介绍在开始细致讲解前，先对整体的操作步骤进行一个大致总结，以便各位观众按需研究：以文本形式打开.inp文件，将管网水源节点的余氯时间序列添加为一个Pattern；在EPANET中打开这个.inp文件，将该Pattern与水源节点相关联，并进行相关设置；使用WNTR完成水力模拟并保存所有节点的水质数据。管网数据文件格式的调整在.inp文件中，可以将管网数据以文本的形式进行编辑，非常方便。那么第一步需要进行的操作就是：如果你的管网文件是.net格式，需要将其调整为.inp格式。这里有个注意点，将.net文件另存为.inp文件不能直接使用 File → Save As... 命令来操作，否则生成的文件会显示打开错误：正确的操作步骤为：File → Export → Network，这样导出的.inp文件才是完整并且可直接进行文本编辑的。将时间序列设置为Pattern完成上面的.inp文件格式转换后，就可以直接使用文本编辑工具（Notepad、VS Code等）打开并进行编辑了。在文件中直接搜索[PATTERNS]，定位到Pattern的存储位置，会有如下的字样：[PATTERNS] ;ID Multipliers ; 1 0.781 0.752 0.738 0.706 0.704 0.73 然后即可在后面进行Pattern的添加，例如：[PATTERNS] ;ID Multipliers ; 1 0.781 0.752 0.738 0.706 0.704 0.73 ; demo 1 0.986378898 0.974703668 0.976000916 0.980541283 0.97340642上面添加的Pattern名为demo，可以使用EPANET打开该文件进行查看：这里有个小问题，如果我的水源节点历史时间序列长度为2个月（60天）、每5分钟一个数据，那么一共会有60*24*(60/5)=17280条数据，也就是说我设置的Pattern长度也必须为17280。同时，还需要进行对应的时间设置：另外还需要将水源节点的余氯值与这个Pattern相关联：关于节点的两个属性：Initial Quality代表了水力模型起始运行那一刻的水质；Source Quality代表了水力模拟过程中，水源水质的变化情况。为什么要在EPANET中进行这些设置呢？就是因为我是用WNTR倒腾了几周还是会有各种各样的问题，一直没有解决，然后才采用了这种方法。在EPANET中设置的属性是会保存在.inp文件中的，但水力模拟的结果或者过程信息却是不会保存在文件中的。在WNTR中进行水力模拟其实，按照上面的步骤进行操作后，我们已经可以在EPANET中进行水力模拟，得到想要的结果了。但是问题就在于：在 Report → Table 的导出选项中，每次仅可以选择一个节点或者管道的属性进行导出：这如果一个网络拓扑结果中有成千上万个节点，那我得导出到猴年马月呢~😟因此，在这里我们借助于WNTR进行水力模拟，其可以得到Numpy数据格式的文件，非常方便进行数据分析与处理。使用WNTR进行水力模拟就非常简单了，直接上代码：# 导入WNTR包 import wntr # 加载EPANET模型 inp_file = r"C:\Users\myxc\Desktop\Anytown-2.inp" wn = wntr.network.WaterNetworkModel(inp_file) # 运行水力模拟(使用EPANET内核) sim = wntr.sim.EpanetSimulator(wn) results = sim.run_sim() # 获取所有节点的余氯数据 cl_concentration = results.node['quality'] * 1000 cl_concentration请注意：在WNTR中设置Pattern或者水源节点的一些具体属性会非常麻烦，但是应该也是可以进行的。本文作者折腾了很长时间，还是不能完美实现自己的需求。因此，才采用了在EPANET中进行属性设置、在WNTR中进行水力模拟的方法。其中，cl_concentration返回结果便是管网中所有节点在任意时刻的余氯监测值（模拟）：但是目前还有个疑问我没有搞明白：可以看到我上面对返回的数据矩阵进行了 * 1000操作，WNTR返回的数据单位需要 * 1000才与EPANET的相等。😕暂时还没找到原因~