[流体输配管网] 树状管网水力分析、牛顿-拉夫森算法/哈代-克罗斯算法求解环能量方程组、牛顿-拉夫森算法求解节点校正压力方程组 的 Matlab 实现与例题求解

1.版权声明

% 管网的水力学特性分析
% Analysis of Hydraulic Characteristics of Pipe Network
% 作者：中山大学 2019 级 廉
% Author: SYSU 2019 CheapMiao

% 参考书目：《给水排水管网系统 第四版》 刘遂庆 9787112254743 中国建筑工业出版社

% 如果原题目中的节点编号不是从 1 开始的连续整数，那么需要重新编号

% 管网恒定流方程组求解条件
% 1.节点流量或压力必须有一个已知
% 2.管网中至少有一个定压节点

% mathcad 的求解器确实方便……

2.函数库

function hf = HazenWilliamsFormula(q, L, D, Cw, n)
    % 海曾威廉公式 提一个 q 出来 hf=10.67/Cw^1.852/D^4.87*L*q*abs(q)^(n-1)
    % 管段流量 q 单位 m^3/s
    % 管段长度 L 单位 m
    % 管段直径 D 单位 m
    % 管材粗糙系数(海曾-威廉系数) Cw 
    
    hf = 10.67./Cw.^1.852./D.^4.87.*L.*q.*abs(q).^(n-1);
    % D = 0 的管段的水头损失定义为 0
    for i = 1:1:size(D,1)
        if abs(D(i)) < 1e-10
            hf(i) = 0;
        end
    end
end

function hf = ManningFormula(q, L, D, nM)
    % 曼宁公式 hf=10.29*n^2*q^2/D^5.333*L
    % 管段长度 L 单位 m
    % 管段直径 D 单位 m
    % 曼宁粗糙系数 nM
    
    hf = 10.29.*nM.^2.*q.^2./D.^(5.333).*L;
    % D = 0 的管段的水头损失定义为 0
    for i = 1:1:size(D,1)
        if abs(D(i)) < 1e-10
            hf(i) = 0;
        end
    end
end

function hp = HydraulicCharacteristicFormula(he, sp, qp, n)
    % 水力特性公式 提一个 q 出来 hp=he-sp*qp*abs(qp)^(n-1)
    % 水泵扬程 hp 单位 m
    % 水泵流量 qp 单位 m^3/s
    % 水泵静扬程 he 单位 m
    % 水泵内阻系数 sp
    % 与水头损失计算指数公式相同的指数 n
    
    hp = he-sp.*qp.*abs(qp).^(n-1);
end

function [vector_var, vector_var_subnum] = FindVars(vector)
    % vector 是一个包含数值和符号变量的符号类型的列向量
    
    % vector 之中的变量的记录数组
    vector_var = [];
    % vector 之中的变量的下标的记录数组
    vector_var_subnum = [];
    
    % 寻找 vector 之中的变量
    for i = 1:1:size(vector,1)
        % 如果可以 double，说明是赋的数值
        try
            % tmp 没用，只是接收值而已
            tmp = double(vector(i));
        % 如果不可以 double，说明是赋的变量，记录变量和下标
        catch ME
            vector_var = [vector_var;vector(i)];
            vector_var_subnum = [vector_var_subnum;i];
        end
    end
end

function [q, Q] = NodalFlowEqsSolver(A, q, Q)
    % 节点流量方程组求解器
    % 节点流量方程组 A*q+Q=0
    % 关联矩阵 A 已知
    % 管段流量 q 单位 m^3/s 部分为未知数
    % 节点流量 Q 单位 m^3/s 部分为未知数
    
    % 1.
    % 现在已经把节点流量方程组解出来了
    % 但是返回的解是一个结构体，不能通过 double 强制类型转换
    % 我还以为这个结构体可以像 lua 那样直接 table[key] 访问，但是不行
    
    % 这样的话，我想要自动地把求解结果放到一个向量里面就有点难办
    % 它 matlab 的 结构体居然不支持通过字符串键访问值
    % 必须要显式的圆点访问
    % 也就是说，我解出了 q，我还需要 q = [Sol.q1 Sol.q2 ... Sol.q9] 这样赋值
    % 这灵活性太低了……比如 q 的长度一变，就要重新改代码
    % 而且参与求解的变量里面还不单单有 q，修改代码的时候还要考虑 Q 的赋值
    % 那么接着就是要考虑 q 的位置和 Q 的位置，好麻烦
    
    % 我想着有没有一个办法让结构体不返回一个结构体，而是直接返回各个值
    % 官方的文档里面是有这样的例子的
    
    % % Sol 会返回 1*1 的结构体
    % Sol = solve(eq_NodalFlow,vars);
    % 
    % % Sol 会返回 1*1 的结构体
    % Sol = zeros(size(vars,1),1);
    % Sol = solve(eq_NodalFlow,vars);
    % 
    % % Sol 会返回 1*1 的结构体
    % Sol = sym(zeros(size(vars,1),1));
    % Sol = solve(eq_NodalFlow,vars);
    % 
    % % Sol 会返回 1*1 的结构体
    % syms Sol [size(vars,1),1];
    % Sol = solve(eq_NodalFlow,vars);
    
    % 但是后面我没试出来返回多个值的情况
    
    % 然后我就想着使用 eval
    % 一开始的想法是用 eval 返回一个变量名，如下代码所示
    
    % % 会报错：动态结构体引用的参数的计算结果必须为有效字段名称。
    % Sol = solve(eq_NodalFlow,vars);
    % Sol.(eval('q1'));
    
    % 这个代码的问题在于，eval 实际上是返回了一个变量，那么实际上传入 Sol 的 fieldname 是一个变量
    % 但是传入 Sol 的 fieldname 应该是一个常量文本
    
    % 因此接下来的方法就是要更进一步，把所有东西都放到 eval 里面！这就是！奥义·动态代码！
    
    % 首先是要动态生成文本
    
    % ['hello','world'] 自动拼接成一个字符串
    % ['Sol.',string(vars(i,1))] string(['Sol.',string(vars(i))]) 是一个 1*2 的 string 类型的数组，没有自动拼接
    % 所以后者需要使用 strcat
    
    % 然后要知道 q Q 里面哪些未知的符号变量
    % 由于传入的 q Q 都是 sym 类型的，所以不能通过 class 判断哪个是赋的数值
    % 但是可以通过类型转换的方式试探，如果能够 double 那就是赋的数值，如果不是，那就是变量
    % 这还是写振动力学求解 w 的时候带来的启发
    % 因此写了 FindVars 函数
    
    % 2.
    % 如果只有一个变量，那么 solve 可以直接解得一个数值
    % 因此需要 try 一下 Sol 是否只是一个数值，代码如下
    
    % % 如果 Sol 的类型为符号，说明只有一个变量
    % if class(Sol) == 'sym'
    %     % 转换成 double 类型
    %     try
    %         Sol = double(Sol);
    %     % 不能转换成 double 类型
    %     catch ME
    %         
    %     end
    % end
    
    % 对 例5.1 运行，报错为：
    % 矩阵维度必须一致。
    % 
    % 出错 AnalysisOfHydraulicCharacteristicsOfPipeNetwork>NodalFlowEqsSolver (第 284 行)
    %     if class(Sol) == 'sym'
    
    % debug 代码如下
    
    % class(Sol)
    % % 如果 Sol 的类型为符号，说明只有一个变量
    % if class(Sol) == 'sym'
    %     1
    %     % 转换成 double 类型
    %     try
    %         Sol = double(Sol);
    %     % 不能转换成 double 类型
    %     catch ME
    %         
    %     end
    % end
    
    % 实际没有输出 1，说明对于结构体，不适用 class(Sol) == 'sym'
    % 但是可以输出 ans = 'struct'，说明对于结构体可以使用 class
    % 那么还有哪里有问题呢……没搞懂
    % 只好把 if class(Sol) == 'sym' 放到 try 里面了，看上去很怪，但是没办法
    
    % 3.
    % 如果 Sol 中出现了 0-1 empty sym，说明该方程没有解
    
    % q Q 之中的变量的下标的记录数组
    [q_var, q_var_subnum] = FindVars(q);
    [Q_var, Q_var_subnum] = FindVars(Q);
    
    % 节点流量方程组
    eq_NodalFlow = A*q+Q == 0;
    % 待求解的变量
    vars = [q_var;Q_var];
    % 求解方程组
    Sol = solve(eq_NodalFlow,vars);
    
    % 转换成 double 类型
    try
        % 如果 Sol 的类型为符号，说明只有一个变量
        if class(Sol) == 'sym'
            Sol = double(Sol);
            % 把求解得到的变量从 Sol 之中拿出来
            if size(q_var_subnum,1) == 1
                q(q_var_subnum(1)) = Sol;
            elseif size(Q_var_subnum,1) == 1
                Q(Q_var_subnum(1)) = Sol;
            end
        end
    % 不能转换成 double 类型
    catch ME
        if class(Sol) == 'struct'
            % 把求解得到的变量从 Sol 结构体之中拿出来 
            for i = 1:1:size(q_var,1)
                % 检查解的存在性
                if size(eval(strcat('Sol.',string(q_var(i)))),1) == 0
                    error('该方程无解');
                end
                % 存在解，拿出
                q(q_var_subnum(i)) = eval(strcat('Sol.',string(q_var(i))));
            end
            for i = 1:1:size(Q_var,1)
                % 检查解的存在性
                if size(eval(strcat('Sol.',string(Q_var(i)))),1) == 0
                    error('该方程无解');
                end
                % 存在解，拿出
                Q(Q_var_subnum(i)) = eval(strcat('Sol.',string(Q_var(i))));
            end
        end
    end
    
    % 转换成 double 类型
    q = double(q);
    Q = double(Q);
end

function H = PressureDropEqsSolver(A, h, H)
    % 管段压降方程组 A'*H=h
    % 关联矩阵 A 
    % 节点压降 H 单位 m
    % 管段压降 h 单位 m 默认全部已知
    
    % 1.
    % 节点水头有时候解不出来，可能就是单纯的系数矩阵的秩与增广矩阵的秩不同
    % 这样的话就必须要手算了……
    % 严格来说节点流量方程似乎也会遇到这样情况？
    
    % 2.
    % 先根据已知变量的下标，把 A' 中已知为数值的那列去掉，保留其他，作为新的系数矩阵 co
    % 使用 A' 中已知为数值的那列与 H 相乘，移项到右边，最后得到新的列向量 b
    % 根据 x=Ab 可以把作为变量的 H 求出来，接下来就是把解放到 H 中
    
    % 3.
    % 似乎不管什么情况都可以使用 x=Ab 来解，这样的话，solve 似乎完全没有必要了
    % 我只是担心有的时候还是需要解更多奇怪的变量
    
    % q Q 之中的变量的下标的记录数组
    [H_var, H_var_subnum] = FindVars(H);
    
    % 对于 Ax=b
    % 构造系数矩阵 A
    for i = 1:1:size(H_var_subnum,1)
        co(:,i) = A(H_var_subnum(i),:);
    end
    % 列向量 b 初始化
    b = h;
    % 构造列向量 b
    for i = 1:1:size(A',2)
        % 如果在下表记录数组中找到了 i，说明 i 对应了变量，不变
        if find(H_var_subnum==i)
            % null
        % 如果在下表记录数组中没有找到 i，说明 i 对应了已知量，需要移项到右边  
        else
            % 需要去掉 A' 的第 i 列，即去掉 A 的第 i 行
            tmp = A(i,:);
            % 移项到右边
            for j = 1:1:size(tmp,2)
                b(j) = b(j) - tmp(j)*H(i);
            end
        end
    end
    % 构造增广矩阵
    augmented = [co,b];
    
    % 如果系数矩阵的秩与增广矩阵的秩相同，说明有解，可以直接用 solve 解
    if rank(co) == rank(augmented)
        % 管段压降方程组
        eq_PressureDrop = A'*H == h;
        % 待求解的变量
        vars = [H_var];
        % 求解方程组
        Sol = solve(eq_PressureDrop,vars);
        
        % 转换成 double 类型
        try
            % 如果 Sol 的类型为符号，说明只有一个变量
            if class(Sol) == 'sym'
                Sol = double(Sol);
                % 把求解得到的变量从 Sol 之中拿出来
                if size(H_var_subnum,1) == 1
                    H(H_var_subnum(1)) = Sol;
                end
            end
        % 不能转换成 double 类型
        catch ME
            if class(Sol) == 'struct'
                 % 把求解得到的变量从 Sol 结构体之中拿出来
                for i = 1:1:size(H_var,1)
                    % 检查解的存在性
                    if size(eval(strcat('Sol.',string(H_var(i)))),1) == 0
                        error('该方程无解');
                    end
                    % 存在解，拿出
                    H(H_var_subnum(i)) = eval(strcat('Sol.',string(H_var(i))));
                end
            end
        end
    % 如果系数矩阵的秩与增广矩阵的秩不同，说明无解，不能用 solve 解
    else
        % 解 Ax=b
        H_tmp = cob;
        % 把求解得到的变量放回 H
        for i = 1:1:size(H_tmp)
            H(H_var_subnum(i)) = H_tmp(i);
        end
    end
    
    % 转换成 double 类型
    H = double(H);
end

function sf = GetPipeFrictionLossCoefficient(L, D, Cw, varargin)
    % 通过指数形式的海曾威廉公式计算管段的摩阻系数 sf=10.67/Cw^1.852/D^4.87*L
    % 通过指数形式的曼宁公式计算管段的摩阻系数 sf=10.29*n^2/D^5.333*L
    % 海曾威廉公式 hf=10.67/Cw^1.852/D^4.87*L*q*abs(q)^(n-1)
    % 曼宁公式 hf=10.29*n^2*q^2/D^5.333*L
    % 水头损失公式的指数形式 hf=sf*q^n
    % 管段长度 L 单位 m
    % 管段直径 D 单位 m
    % 管材粗糙系数(海曾-威廉系数) Cw 
    % 曼宁粗糙系数 nM
    
    % 通过指数形式的曼宁公式计算管段的摩阻系数
    if size(varargin,1)
        if varargin{1} == 'Manning'
            n = Cw;
            sf = 10.29.*n.^2./D.^5.333.*L;
        end
    % 默认通过指数形式的海曾威廉公式计算管段的摩阻系数
    else
        sf = 10.67./Cw.^1.852./D.^4.87.*L;
    end
    % D = 0 的管段的水头损失定义为 0
    for i = 1:1:size(D,1)
        if abs(D(i)) < 1e-10
            sf(i) = 0;
        end
    end
end

function s = GetPipeResistanceCoefficient(L, D, sp, Cw, varargin)
    % 计算管段阻力系数 = 管段的摩阻系数 + 泵站的摩阻系数
    % 管段阻力系数 s
    % 管段长度 L 单位 m
    % 管段直径 D 单位 m
    % 水泵内阻系数(泵站的摩阻系数) sp
    % 管材粗糙系数(海曾-威廉系数) Cw
    % 曼宁粗糙系数 nM
    
    % 通过指数形式的曼宁公式计算管段的摩阻系数
    if size(varargin,1)
        if varargin{1} == 'Manning'
            sf = GetPipeFrictionLossCoefficient(L, D, Cw, 'Manning');
        end
    % 默认通过指数形式的海曾威廉公式计算管段的摩阻系数
    else
        sf = GetPipeFrictionLossCoefficient(L, D, Cw);
    end
    s = sf+sp;
end

function z = GetPipeDampingCoefficient(s, q, n)
    % 计算管段阻尼系数 z=n*s*abs(q)^(n-1);
    % 管段阻尼系数 z 
    % 管段阻力系数 s
    % 管段流量 q 单位 m^3/s
    % 与水头损失计算指数公式相同的指数 n
    
    z = n.*s.*abs(q).^(n-1);
end

function delta_h = GetWaterheadClosureError(Bf, h)
    % 计算水头闭合差 Δh=Bf*h
    % 水头闭合差 Δh 单位 m
    % 基本回路矩阵 Bf 顺时针为正
    % 管段压降 h 单位 m
    
    delta_h = Bf*h;
end

function bool = IsErrorAllowable(delta_h, eh) 
    % 例：判断各水头闭合差是否都小于最大允许闭合差
    % 水头闭合差 Δh 单位 m
    % 环水头闭合差的最大允许值 eh 单位 m
    
    bool = true;
    for i = 1:1:size(delta_h,1)
        if abs(delta_h(i)) > eh
            bool = false;
            return;
        end
    end
end

function F = GetLoopFlowCorrectionCoefficient(Bf, z)
    % 牛顿-拉夫森算法 计算与环校正流量的非线性函数等效的系数的矩阵 
    %     Σ(i∈Rl)(zi(0)), l = j,系数矩阵的对角元素
    % F = -zi(0),          i 为相邻环 l 和 j 的公共管段
    %     0，              l ≠ j 且环 l 和 j 不相邻
    % 其实这个书上的系数矩阵的计算公式有一点不对，就是计算公共管段的时候，假设了公共管段只有一段
    % 实际上有若干段是公共的完全是可能的嘛
    % 因此正确公式应该是
    %     Σ(i∈Rl)(zi(0)),   l = j,系数矩阵的对角元素
    % F = -Σ(i∈Rl)(zi(0)),  i 为相邻环 l 和 j 的公共管段
    %     0，                l ≠ j 且环 l 和 j 不相邻
    % 基本回路矩阵 Bf 顺时针为正
    % 管段阻尼系数 z
    
    % 基本回路矩阵的绝对值
    Bf_abs = abs(Bf);
    % 系数矩阵的对角线
    F_dia = Bf_abs*z;
    % 系数矩阵初始化
    F = diag(F_dia);
    % 获得系数矩阵的非对角元素
    for i = 1:1:size(F,1)
        for j = 1:1:size(F,2)
            % 非对角线
            if i ~= j
                % 第 i 个环和第 j 个环的公共边列向量
                pipe_shared = Bf_abs(i,:).*Bf_abs(j,:);
                % 系数矩阵的非对角元素
                F(i,j) = -pipe_shared*z;
            end
        end
    end
end

function q = CorrectQuantityOfFlow(Bf, q, delta_q)
    % 矫正各管段流量 qi(k+1)=qi(k)±Δql(k), i∈Rl
    % 基本回路矩阵 Bf 顺时针为正
    % 管段流量 q 单位 m^3/s
    % 管段矫正流量 delta_q 单位 m^3/s
    
    for i = 1:1:size(Bf,1)
        q = q + delta_q(i)*Bf(i,:)';
    end
end

function [hf, hp, h] = NewtonRaphsonProcess_LoopFlowCorrection(Bf, q, L, D, s, he, sp, Cw, n, eh, varargin)
    % 牛顿-拉夫森算法 环流量校正 判断各水头闭合差是否都小于最大允许闭合差，否则迭代
    % 基本回路矩阵 Bf 顺时针为正
    % 管段流量 q 单位 m^3/s
    % 管段长度 L 单位 m
    % 管段直径 D 单位 m
    % 管段阻力系数 s
    % 水泵静扬程 he 单位 m
    % 水泵内阻系数 sp
    % 管材粗糙系数(海曾-威廉系数) Cw 或曼宁公式系数 nM
    % 与水头损失计算指数公式相同的指数 n
    % 环水头闭合差的最大允许值 eh 单位 m
    
    % 牛顿-拉夫森算法 判断各水头闭合差是否都小于最大允许闭合差，否则迭代
    % do...while 型
    while 1
        % 通过指数形式的曼宁公式计算各管段水头损失
        if size(varargin,1)
            if varargin{1} == 'Manning'
                nM = Cw;
                hf = ManningFormula(q, L, D, nM);
            end
        % 默认通过指数形式的海曾威廉公式计算各管段水头损失
        else
            hf = HazenWilliamsFormula(q, L, D, Cw, n);
        end
        % 计算各管段泵站扬程
        hp = HydraulicCharacteristicFormula(he, sp, q, n);
        % 计算各管段压降
        h = hf - hp;
        % 计算管段阻尼系数 z
        z = GetPipeDampingCoefficient(s, q, n);
        % 计算与环校正流量的非线性函数等效的系数的矩阵
        F = GetLoopFlowCorrectionCoefficient(Bf, z);
        % 计算水头闭合差
        delta_h = GetWaterheadClosureError(Bf, h);
        
        % 求解基本环能量方程组
        delta_q = -Fdelta_h;
        % 矫正各管段流量（平差计算）
        q = CorrectQuantityOfFlow(Bf, q, delta_q);

        % 误差允许，结束迭代
        if IsErrorAllowable(delta_h, eh)
            break;
        end
    end
end

function [hf, hp, h] = HardyClausProcess_LoopFlowCorrection(Bf, q, L, D, s, he, sp, Cw, n, eh, varargin)
    % 哈代-克罗斯算法 环流量校正 判断各水头闭合差是否都小于最大允许闭合差，否则迭代
    % 基本回路矩阵 Bf 顺时针为正
    % 管段流量 q 单位 m^3/s
    % 管段长度 L 单位 m
    % 管段直径 D 单位 m
    % 管段阻力系数 s
    % 水泵静扬程 he 单位 m
    % 水泵内阻系数 sp
    % 管材粗糙系数(海曾-威廉系数) Cw 或曼宁公式系数 nM
    % 与水头损失计算指数公式相同的指数 n
    % 环水头闭合差的最大允许值 eh 单位 m
    
    % 哈代-克罗斯算法 判断各水头闭合差是否都小于最大允许闭合差，否则迭代
    % do...while 型
    while 1
        % 通过指数形式的曼宁公式计算各管段水头损失
        if size(varargin,1)
            if varargin{1} == 'Manning'
                nM = Cw;
                hf = ManningFormula(q, L, D, nM);
            end
        % 默认通过指数形式的海曾威廉公式计算各管段水头损失
        else
            hf = HazenWilliamsFormula(q, L, D, Cw, n);
        end
        % 计算各管段泵站扬程
        hp = HydraulicCharacteristicFormula(he, sp, q, n);
        % 计算各管段压降
        h = hf - hp;
        % 计算管段阻尼系数 z
        z = GetPipeDampingCoefficient(s, q, n);
        % 计算水头闭合差
        delta_h = GetWaterheadClosureError(Bf, h);
        % 基本回路矩阵的绝对值
        Bf_abs = abs(Bf);
        % 系数矩阵的对角线
        F_dia = Bf_abs*z;
        
        % 求解基本环能量方程组
        delta_q = -delta_h./F_dia;
        % 矫正各管段流量（平差计算）
        q = CorrectQuantityOfFlow(Bf, q, delta_q);
        
        % 误差允许，结束迭代
        if IsErrorAllowable(delta_h, eh) 
            break;
        end
    end    
end

function [hf, hp, h] = ImprovedHardyClausProcess_LoopFlowCorrection(Bf, q, L, D, s, he, sp, Cw, n, eh, varargin)
    % 改进哈代-克罗斯法 环流量校正
    % 基本回路矩阵 Bf 顺时针为正
    % 管段流量 q 单位 m^3/s
    % 管段长度 L 单位 m
    % 管段直径 D 单位 m
    % 管段阻力系数 s
    % 水泵静扬程 he 单位 m
    % 水泵内阻系数 sp
    % 管材粗糙系数(海曾-威廉系数) Cw 或曼宁公式系数 nM
    % 与水头损失计算指数公式相同的指数 n
    % 环水头闭合差的最大允许值 eh 单位 m
    
    % 1.
    % 改"各环同时平差"为"每次只平差一个环"
    % 各环同时平差时, 每个环平差时未考虑前面已平差环传来的闭合差, 不能将其传递出去, 故收敛速度慢
    % 且同时平差不能选择平差顺序, 不能避免对闭合差已经小于允许值的环进行多元的平差计算
    % 每次只平一个环, 平差后立即更新管段流量, 后面的平差利用新的管段流量计算闭合差, 则平差收敛速度得到很大提高
    % 平差时优先平差闭合差最大的环, 这样传递出去的闭合差最大, 计算效率最高
    % 平差接近最终解时管段流量变化量已经很小, 故在若干次平差后可将阻尼系数固定不变从而节省部分计算量
    % [1]张涛,蔡毅,李湧.环状管网自动平差研究[J].山西建筑,2008,(19):364-365.
    
    % 改进哈代-克罗斯算法 判断各水头闭合差是否都小于最大允许闭合差，否则迭代
    % do...while 型
    while 1
        % 通过指数形式的曼宁公式计算各管段水头损失
        if size(varargin,1)
            if varargin{1} == 'Manning'
                nM = Cw;
                hf = ManningFormula(q, L, D, nM);
            end
        % 默认通过指数形式的海曾威廉公式计算各管段水头损失
        else
            hf = HazenWilliamsFormula(q, L, D, Cw, n);
        end
        % 计算各管段泵站扬程
        hp = HydraulicCharacteristicFormula(he, sp, q, n);
        % 计算各管段压降
        h = hf - hp;
        % 计算管段阻尼系数 z
        z = GetPipeDampingCoefficient(s, q, n);
        % 计算水头闭合差
        delta_h = GetWaterheadClosureError(Bf, h);
        % 基本回路矩阵的绝对值
        Bf_abs = abs(Bf);
        % 系数矩阵的对角线
        F_dia = Bf_abs*z;
        % 寻找闭合差绝对值最大的环
        [~, delta_h_max_subnum] = max(abs(delta_h));
        
        % 求解基本环能量方程组 但只解闭合差绝对值最大的环
        delta_q = zeros(size(delta_h));
        delta_q(delta_h_max_subnum) = -delta_h(delta_h_max_subnum)./F_dia(delta_h_max_subnum);
        % 矫正各管段流量（平差计算）
        q = CorrectQuantityOfFlow(Bf, q, delta_q);
        
        % 误差允许，结束迭代
        if IsErrorAllowable(delta_h, eh) 
            break;
        end
    end
end

function c = GetNodalPressureCorrectionCoefficient(A, q, s, H_var_subnum, n)
    % 牛顿-拉夫森算法 计算与节点校正流量的非线性函数等效的系数的矩阵 
    %           Σ(k∈Si)(1/(n*sjk*|qjk(0)|^(n-1))), 系数矩阵的主对角元素, Sj——节点 j 的关联集
    % ∂G/∂ΔHk = -1/(n*sjk*|qjk(0)|^(n-1))         , 节点 j 和 k 衔接, 第 j 行第 k 列元素及其对称元素
    %           0                                 , 节点 j 和 k 不衔接, 第 j 行第 k 列元素及其对称元素
    % 关联集就是看与这个点相连的其他节点
    % 衔接就是只看 j 和 k 决定的这个管段
    % 关联矩阵 A
    % 管段流量 q 单位 m^3/s
    % 管段阻力系数 s
    
    % 1.
    % 噢！这个公式分母不是 q^n 而是 q^(n-1)
    % 所以公式不是
    % c(i,j) = c(i,j) + 1./n./s(col(k))./q(col(k))./abs(q(col(k)))^(n-1);
    % c(i,j) = -1./n./s(col(k))./q(col(k))./abs(q(col(k)))^(n-1);
    % 而是
    % c(i,j) = c(i,j) + 1./n./s(col(k))./abs(q(col(k)))^(n-1);
    % c(i,j) = -1./n./s(col(k))./abs(q(col(k)))^(n-1);
    
    % 系数矩阵初始化
    c = zeros(size(H_var_subnum,1),size(H_var_subnum,1));
    % 构造系数矩阵
    for i = 1:1:size(H_var_subnum,1)
        for j = 1:1:size(H_var_subnum,1)
            % 对角线元素，看关联集
            % H_var_subnum(i) 是由 H_var_subnum 决定的第几个 H 的序号
            if i == j
                % 从关联矩阵中第 H_var_subnum(i) 个节点对应的第 H_var_subnum(i) 行中找到含非零元素的列数 col
                % col 即表示与第 H_var_subnum(i) 个节点相连的管段编号
                col = find(A(H_var_subnum(i),:)); 
                % 对于每一个与第 i 个节点相连的管段
                for k = 1:1:size(col,2)
                    c(i,j) = c(i,j) + 1./n./s(col(k))./abs(q(col(k)))^(n-1);
                end
            % 非对角线元素，看衔接
            else
                % 从关联矩阵中第 H_var_subnum(i) 个节点对应的第 H_var_subnum(i) 行中找到含非零元素的列数 col
                % col 即表示与第 H_var_subnum(i) 个节点相连的管段编号
                col = find(A(H_var_subnum(i),:));
                % 对于每一个与第 H_var_subnum(i) 个节点相连的管段
                for k = 1:1:size(col,2)
                    % 如果在第 k 个与第 H_var_subnum(i) 个节点相连的管段中找到了非零元素
                    % 那么检查这个元素是否是第 j 个节点，即这个元素是否在第 j 行
                    if find(find(A(:,col(k)))==j)
                        c(i,j) = -1./n./s(col(k))./abs(q(col(k)))^(n-1);
                    end
                end    
            end
        end
    end
end

function [q, H] = NewtonRaphsonProcess_NodalPressureCorrection(A, H, H_var_subnum, Q, L, D, sp, Cw, n, eq, varargin)
    % 牛顿-拉夫森算法 节点压力校正 判断各水头闭合差是否都小于最大允许闭合差，否则迭代
    % 关联矩阵 A 
    % 节点压降 H 单位 m
    % 管段流量 q 单位 m^3/s
    % 水泵内阻系数 sp
    % 管段阻力系数 s
    % 管材粗糙系数(海曾-威廉系数) Cw 或曼宁公式系数 nM
    % 与水头损失计算指数公式相同的指数 n
    % 节点流量闭合差的最大允许值 eq 单位 L/s 或节点压力的最大允许值 ε
    
    % 1.
    % 节点压力平差计算步骤
    % 1.确定已知节点压力，设定位置压力节点的初始压力 Hi(0)
    % 2.用当前节点压力 Hi 计算各管段水头损失和管段流量 hj=HFj-HTj 和 qj=(hj/sj)^(1/n),(j=1,2,...,M)
    % 3.计算各节点流量闭合差 ΔQi=Σ(j∈Si)((aj)*qj)+Qi,(i=1,2,...,N)
    % aj = 1,当节点 i 是管段 j 的出点；aj = -1,当节点 i 是管段 j 的入点
    % 4.计算节点校正压力 
    % 牛顿-拉夫森算法 求解节点校正压力方程组，得到其系数矩阵
    %           Σ(k∈Si)(1/(n*sjk*|qjk(0)|^(n-1))), 系数矩阵的主对角元素, Sj——节点 j 的关联集
    % ∂G/∂ΔHk = -1/(n*sjk*|qjk(0)|^(n-1))         , 节点 j 和 k 衔接, 第 j 行第 k 列元素及其对称元素
    %           0                                 , 节点 j 和 k 不衔接, 第 j 行第 k 列元素及其对称元素
    % 节点压力平差法 ΔHi=-ΔQi/Σ(j∈Si)(cij),(i=1,2,...,N)
    % 5.如果任一节点校正压力 ΔHi(k)>ε，计算新的节点压力：Hi(k+1)=Hi(k)+ΔHi(k),(i=1,2,...,N,k 为迭代计算次数)，返回步骤 2
    % 6.如果 ΔHi(k)<ε，平差计算完成，求管段流量和节点自由压力
    
    % 2.
    % 本来打算通过如下代码传入计算系数矩阵的参数
    
    % % 向量初始化
    % q_var = zeros(size(H_var_subnum));
    % s_var = zeros(size(H_var_subnum));
    % % 构造
    % for i = 1:1:size(H_var_subnum,1)
    %     q_var(i) = q(H_var_subnum(i));
    %     s_var(i) = s(H_var_subnum(i));
    % end
    % 
    % % 计算与节点校正流量的非线性函数等效的系数的矩阵
    % c = GetNodalPressureCorrectionCoefficient(A, q_var, s_var, n)
    
    % 但是这时不正确的，因为砍掉的自由度是要砍掉节点，而这种做法砍掉了管段
    % 结果例如，(1)~(7) 的节点计算 c 时需要用到管段 [8]，但是却访问不到
    % 这个时候的 GetNodalPressureCorrectionCoefficient 也是根据已经消掉管段自由度来设计的，所以还要重新做一个
    % 比如向 GetNodalPressureCorrectionCoefficient 传入 H_var_subnum
    
    % 3.
    % 初始时水源节点 (1) 和 (5) 的泵站扬程未知，使得 he 未知，故不能使用如下代码在一开始确定一个 H 的初值
    
    % % 默认通过指数形式的海曾威廉公式计算各管段水头损失
    % hf = HazenWilliamsFormula(q, L, D, Cw, n);
    % % 计算各管段泵站扬程
    % hp = HydraulicCharacteristicFormula(he, sp, q, n);
    % % 计算各管段压降
    % h = hf - hp;
    % % 求解管段压降方程组
    % H = PressureDropEqsSolver(A, h, H)
    
    % 4.
    % 根据 qj=(hj/sj)^(1/n),(j=1,2,...,M) 先计算出 q 时，不能立即根据 Aq+Q=0 算出 Q，然后使用这个 Q 值计算 ΔQ
    % 也就是说，我将 Q 误认为一个具有初值，只受 ΔQ 影响的值
    % 实际上，不能在迭代的过程中通过 Aq+Q=0 算出 Q，因为这个时候还存在着节点流量闭合差
    % 就是说，Aq+Q=0不完全成立，误差是节点流量闭合差
    % 书本上，q 的取值就是取最后一次迭代时计算出的 q
    
    % 取第一次迭代时计算出的 q，根据 Aq+Q=0 计算出 Q
    
    % Q = 8×1    
    %    -0.1051
    %    -0.0221
    %     0.0372
    %     0.0346
    %    -0.0429
    %    -0.0182
    %     0.0903
    %     0.0262
        
    % 如下代码所示，取最后一次迭代时计算出的 q，根据 Aq+Q=0 计算出 Q
    
    % % 关联矩阵 A 
    % A = [
    %      1     0     0     0     0     0     0     0     0
    %     -1     1     0     0     1     0     0     0     0
    %      0    -1     1     0     0     1     0     0     0
    %      0     0    -1    -1     0     0     1     0     0
    %      0     0     0     1     0     0     0     0     0
    %      0     0     0     0    -1     0     0     1     0
    %      0     0     0     0     0    -1     0    -1     1
    %      0     0     0     0     0     0    -1     0    -1];
    % % 课本得出的管段流量结果
    % q_test = [0.0989+0.0989 0.0828 0.0113 0.0168+0.0168 0.0997 0.0203 0.0241 0.0645 0.0026]';
    % % 展示结果
    % disp(-A*q_test)
    
    % >>disp(-A*q_test)
    %    -0.1978
    %     0.0153
    %     0.0512
    %     0.0208
    %    -0.0336
    %     0.0352
    %     0.0822
    %     0.0267
    
    % 而节点流量初值为
    
    % Q = [
    %     -198.08
    %     14.55
    %     51.17
    %     20.77
    %     -33.42
    %     35.03
    %     82.33
    %     27.65]'/1000;
    
    % 这就验证了 q 的取值就是取最后一次迭代时计算出的 q

    % 5.
    % 一开始，第一次迭代计算出来的 ΔQ 和书本中给出的表格数据不符
    % 然后我就验算了一下书本上面是否是对的
    
    % 书本上 q1 = 0.1038 m^3/s q10 = 0.1038 m^3/s Q1 = -0.19808 m^3/s 根据公式 ΔQi=Σ(j∈Si)(qj)+Qi,(i=1,2,...,N) 得
    % ΔQ1 = q1 + q10 + Q1 = 0.1038+0.1038-0.19808 m^3/s = 0.0095 m^3/s
    % 书本表格数据为 ΔQ1 = 0.0095 m^3/s 与计算数据相符
    
    % 书本上 q1 = 0.1038 m^3/s q10 = 0.1038 m^3/s q2 = 0.0708 m^3/s q5 = 0.0564 m^3/s Q2 = 0.01455 m^3/s
    % 根据公式得 ΔQ2 = q1 + q10 + q2 + q5 + Q2 = 0.1038+0.1038+0.0708+0.0564+0.01455 m^3/s = 0.3493 m^3/s
    % 书本表格数据为 ΔQ2 = -0.0611 m^3/s 与计算数据不符
    
    % 我感觉这应该是还是需要在 qi 的前面加一个正负号的系数，书上写少了，比如
    % ΔQ2 = -q1 - q10 + q2 + q5 + Q2 = -0.1038-0.1038+0.0708+0.0564+0.01455 m^3/s = -0.0658 m^3/s
    % 这样就跟书本数据相差不大了
    
    % 因此 公式 ΔQi=Σ(j∈Si)(qj)+Qi,(i=1,2,...,N) 应该还是少了一个正负号，正确的为
    % ΔQi=Σ(j∈Si)((aj)*qj)+Qi,(i=1,2,...,N)
    % aj = 1,当节点 i 是管段 j 的出点；aj = -1,当节点 i 是管段 j 的入点
    
    % 6.
    % 在 gitee 上搜到了一个同样是 给水管网节点压力平差计算 的程序
    % https://gitee.com/q_rain/PipeNetModel/tree/master
    % 但是我看不懂他那个计算公式为什么那么写：
    % Qi = Σ(j∈Si)((aj)*qj)-0.5/dhfdq*H(k)
    % 当前管段的另一端节点的节点编号 k
    % 沿程水损对流量的导数 dhfdq = d(hf)/d(q) = n*sp*|q|^(n-1)

    % 7.
    % 又算错了，但是不知道原因，对比数据也不可靠的时候，还是只能自己手算了
    % 比如我第二个算错了
    % 已知在 delta_Q_var(i) = delta_Q_var(i) + Q(H_var_subnum(i)); 中
    % i=2时，-0.1348+0.0146=-0.1203
    % 但是实际上应为 -0.2056+0.0708+0.0564+0.01455=-0.0638 与书中数据表相近
    % 又发现 -0.2056+0.0708=-0.1348
    % 所以错误就是 delta_Q_var(i) = delta_Q_var(i) + node1*q(col(k)); 少执行了一次
    % 截到那一步，发现本来应该是 +0.0564 的 node1*q(col(k)) 现在算出来是 0
    % 那就说明这个时候算出来的 node1=0，就说明没有在关联矩阵里面找到管段端点
    % 要不是关联矩阵写错了，要不是找点找错了
    % 重新看了一遍关联矩阵，确认没有错，那就是找点找错了
    % 然后就发现是 node1 = A(H_var_subnum(i),k); 写错了，正确的应该是
    % node1 = A(H_var_subnum(i),col(k));
    
    % 通过指数形式的曼宁公式计算各管段水头损失
    if size(varargin,1)
        if varargin{1} == 'Manning'
            nM = Cw;
            % 计算管段阻力系数 s
            s = GetPipeResistanceCoefficient(L, D, sp, nM, 'Manning');
        else
            % 计算管段阻力系数 s
            s = GetPipeResistanceCoefficient(L, D, sp, Cw);
        end
    % 默认通过指数形式的海曾威廉公式计算各管段水头损失
    else
        % 计算管段阻力系数 s
        s = GetPipeResistanceCoefficient(L, D, sp, Cw);
    end
    
    delta_Q_var = zeros(size(H_var_subnum));
    H_var = zeros(size(H_var_subnum));
    % 构造
    for i = 1:1:size(H_var_subnum,1)
        H_var(i) = H(H_var_subnum(i));
    end
    
    % do...while 型
    while(1)
        % 2.用当前节点压力 Hi 计算各管段水头损失和管段流量 hj=HFj-HTj 和 qj=(hj/sj)^(1/n),(j=1,2,...,M)
        h = A'*H;
        q = (h./s).^(1/n);
        % 3.计算各节点流量闭合差 ΔQi=Σ(j∈Si)(qj)+Qi,(i=1,2,...,N)
        % 第 i 个闭合差对应第 H_var_subnum(i) 个节点
        for i = 1:1:size(H_var_subnum,1)
            col = find(A(H_var_subnum(i),:)); 
            % 对于每一个与第 H_var_subnum(i) 个节点相连的管段
            for k = 1:1:size(col,2)
                % node1 为查找结果的第 k 个管段中作为起点的第 H_var_subnum(i) 个节点的值
                % 即 node1 是查找结果的第 k 个管段中的一个端点的值
                node1 = A(H_var_subnum(i),col(k));
                % node1 是 q 前的系数
                delta_Q_var(i) = delta_Q_var(i) + node1*q(col(k));
            end
            delta_Q_var(i) = delta_Q_var(i) + Q(H_var_subnum(i));
        end

        % 计算与节点校正流量的非线性函数等效的系数的矩阵
        c = GetNodalPressureCorrectionCoefficient(A, q, s, H_var_subnum, n);
        % 4.计算节点校正压力
        delta_H_var = -cdelta_Q_var;
        % 通过自定义的节点压力校正公式计算节点压力
        if size(varargin,1)
            if class(varargin{2}) == 'function_handle'
                fun = varargin{2};
                % 节点压力校正公式
                for i = 1:1:size(H_var_subnum,1)
                    H(H_var_subnum(i)) = fun(H(H_var_subnum(i)),delta_H_var(i));
                end
            end
        % 默认节点压力校正公式计算节点压力
        else
            % 节点压力校正公式
            for i = 1:1:size(H_var_subnum,1)
                H(H_var_subnum(i)) = H(H_var_subnum(i)) + delta_H_var(i);
            end
        end
        
        % 误差允许，求解完毕
        if IsErrorAllowable(delta_Q_var,eq)
            return;
        end
    end
end

% 节点压力平差算法求解节点矫正压力方程组的函数就懒得写了，就也只是改一下 c 而已

3.树状管网水力分析

% 书本 例5.1 树状管网水力分析

clear;

% 已知条件

% 关联矩阵 A
A = [
    1	0	0	0	0	0	0	0	0
    -1	1	0	0	0	0	0	0	0
    0	-1	1	0	1	0	0	0	0
    0	0	-1	1	0	0	0	0	0
    0	0	0	-1	0	0	0	0	0
    0	0	0	0	-1	1	0	0	1
    0	0	0	0	0	-1	1	0	0
    0	0	0	0	0	0	-1	1	0
    0	0	0	0	0	0	0	-1	0
    0	0	0	0	0	0	0	0	-1];
% 管段流量 q 单位 L/s 部分为未知数
syms q [size(A,2) 1] real;
% 节点流量 Q 单位 L/s 部分为未知数
syms Q1 real;
Q = sym([Q1;5.37;16.11;7.16;4.48;30.74;7.52;7.07;3.67;11.63]);

% 管段长度 L 单位 m
L = [600;300;150;250;450;230;190;205;650];
% 管段直径 D 单位 m
D = [400;400;150;100;300;200;150;100;150]/1000;
% 水泵静扬程 he 单位 m
he = [42.6;0;0;0;0;0;0;0;0];
% 水泵内阻系数 sp
sp = [311.1;0;0;0;0;0;0;0;0];
% 水泵流量 qp 单位 m^3/s
% 等于管段流量 q

% 设计水头 H 单位 m
syms H [size(A,1) 1] real; 
H(1) = sym(7.8);
% 节点地面标高 H_att 单位 m 
H_att = [9.8;11.5;11.8;15.2;17.4;13.3;12.8;13.7;12.5;15];

% 海曾威廉系数 Cw
Cw = 100;
% 与水头损失计算指数公式相同的指数 n
n = 1.852;

% 求解节点流量方程组
[q, Q] = NodalFlowEqsSolver(A, q, Q);
% 管段流量 q 单位 m^3/s 全部已知
q = q/1000;
% 节点流量 Q 单位 m^3/s 全部已知
Q = Q/1000;
% 管段流速 v 单位 m/s
v = q./(pi*D.^2/4);
% 计算各管段水头损失
hf = HazenWilliamsFormula(q, L, D, Cw, n);
% 计算各管段泵站扬程
hp = HydraulicCharacteristicFormula(he, sp, q, n);
% 计算各管段压降
h = hf - hp;
% 求解管段压降方程组
H = PressureDropEqsSolver(A, h, H);
% 计算自由水压
H_free = H - H_att;

4.管网环方程组水力分析 牛顿-拉夫森算法求解环能量方程组

% 书本 例5.2 管网环方程组水力分析 牛顿-拉夫森算法求解环能量方程组

clear;

% 已知条件

% 题中节点编号 重新编号
% 2           1
% 3           2
% 4           3
% 6           4
% 7           5
% 8           6

% 题中管段编号 重新编号
% 2           1
% 3           2
% 5           3
% 6           4
% 7           5
% 8           6
% 9           7

% 关联矩阵 A 
A = [
     1     0     1     0     0     0     0
    -1     1     0     1     0     0     0
     0    -1     0     0     1     0     0
     0     0    -1     0     0     1     0
     0     0     0    -1     0    -1     1
     0     0     0     0    -1     0    -1];
% 基本回路矩阵 Bf 顺时针为正
Bf = [
     1     0    -1     1     0    -1     0
     0     1     0    -1     1     0    -1];
% 管段流量 q 单位 m^3/s 全部已知
q = [0.08990 0.00627 0.08990 0.03246 0.02265 0.05487 0.00500]';
% 节点流量 Q 单位 m^3/s 部分为未知数
syms Q6 real;
Q = sym([-0.17980 0.05117 -0.01638 0.03503 0.08233 Q6]');

% 管段长度 L 单位 m
L = [650 550 330 350 360 590 490]';
% 管段直径 D 单位 m
D = [0.300 0.200 0.300 0.200 0.20 0.300 0.100]';
% 水泵静扬程 he 单位 m
he = [0 0 0 0 0 0 0]';
% 水泵内阻系数 sp
sp = [0 0 0 0 0 0 0]';
% 水泵流量 qp 单位 m^3/s
% 等于管段流量 q

% 设计水头 H 单位 m
syms H [size(A,1) 1] real; 
H(6) = sym(41.50);
% 节点地面标高 H_att 单位 m 
H_att = [18.80 19.10 22.00 18.30 17.30 17.50]';

% 海曾威廉系数 Cw
Cw = 110;
% 与水头损失计算指数公式相同的指数 n
n = 1.852;
% 环水头闭合差的最大允许值 eh 单位 m
eh = 0.01;

% 求解节点流量方程组
[q, Q] = NodalFlowEqsSolver(A, q, Q);
% 计算管段阻力系数 s
s = GetPipeResistanceCoefficient(L, D, sp, Cw);
% 牛顿-拉夫森算法
[hf, hp, h] = NewtonRaphsonProcess_LoopFlowCorrection(Bf, q, L, D, s, he, sp, Cw, n, eh);
% 管段流速 v 单位 m/s
v = q./(pi*D.^2/4);
% 求解管段压降方程组
H = PressureDropEqsSolver(A, h, H);
% 计算自由水压
H_free = H - H_att;

5.管网环方程组水力分析 哈代-克罗斯算法求解环能量方程组

% 书本 例5.3 管网环方程组水力分析 哈代-克罗斯算法求解环能量方程组

clear;

% 已知条件

% 题中节点编号 重新编号
% 0           1
% 1           2
% 2           3
% 3           4
% 4           5
% 5           6
% 6           7
% 7           8
% 8           9

% 关联矩阵 A 
A = [
     0     0     0     0     0     0     0     0     0     1     1
     1     0     0     0     0     0     0     0     0    -1     0
    -1     1     0     0     1     0     0     0     0     0     0
     0    -1     1     0     0     1     0     0     0     0     0
     0     0    -1    -1     0     0     1     0     0     0     0
     0     0     0     1     0     0     0     0     0     0    -1
     0     0     0     0    -1     0     0     1     0     0     0
     0     0     0     0     0    -1     0    -1     1     0     0
     0     0     0     0     0     0    -1     0    -1     0     0];
% 基本回路矩阵 Bf 顺时针为正
Bf = [
     0     1     0     0    -1     1     0    -1     0     0     0
     0     0     1     0     0    -1     1     0    -1     0     0
     1     1     1    -1     0     0     0     0     0     1    -1];
% 管段流量 q 单位 m^3/s 全部已知
q = [198.08 81.17 10.00 33.42 102.36 20.00 22.65 67.33 5.00 198.08 33.42]'/1000;
% 节点流量 Q 单位 m^3/s 全部已知
% 不需要求解
% 管段长度 L 单位 m
L = [320 650 550 270 330 350 360 590 490 0 0]';
% 管段直径 D 单位 m
D = [2^(2/5.333)*300 300 200 2^(2/5.333)*200 300 200 200 300 100 0 0]'/1000;
% 水泵静扬程 he 单位 m
he = [48.96 0 0 0 0 0 0 0 0 12 48]';
% 水泵内阻系数 sp
sp = [138.5 0 0 0 0 0 0 0 0 0 0]';
% 水泵流量 qp 单位 m^3/s
% 等于管段流量 q

% 设计水头 H 单位 m
syms H [size(A,1) 1] real; 
H(1) = sym(0);
H(5) = sym(48.00);
% 节点地面标高 H_att 单位 m 
H_att = [0 13.60 18.80 19.10 22.00 32.20 18.30 17.30 17.50]';

% 海曾威廉系数 Cw
Cw = 110;
% 与水头损失计算指数公式相同的指数 n
n = 1.852;
% 环水头闭合差的最大允许值 eh 单位 m
eh = 0.01;

% 计算管段阻力系数 s
s = GetPipeResistanceCoefficient(L, D, sp, Cw);
s = 1.0e+05 *[
    0.0021
    0.0059
    0.0435
    0.0073
    0.0030
    0.0277
    0.0285
    0.0054
    1.5631
         0
         0];
% 哈代-克罗斯算法
[hf, hp, h] = HardyClausProcess_LoopFlowCorrection(Bf, q, L, D, s, he, sp, Cw, n, eh);
% 管段流速 v 单位 m/s
v = q./(pi*D.^2/4);
% 求解管段压降方程组
H = PressureDropEqsSolver(A, h, H);
% 计算自由水压
H_free = H - H_att;

6.习题1

% 书本 习题5.1 对于 例5.3
% 管段水头损失改用曼宁公式计算，取 n = 0.012
% 将水塔高度提高 0.5m，即节点 (5) 水头改为 48.50m，其余数据不变
% 使用改进哈代-克罗斯平差算法

clear;

% 已知条件

% 题中节点编号 重新编号
% 0           1
% 1           2
% 2           3
% 3           4
% 4           5
% 5           6
% 6           7
% 7           8
% 8           9

% 关联矩阵 A 
A = [
     0     0     0     0     0     0     0     0     0     1     1
     1     0     0     0     0     0     0     0     0    -1     0
    -1     1     0     0     1     0     0     0     0     0     0
     0    -1     1     0     0     1     0     0     0     0     0
     0     0    -1    -1     0     0     1     0     0     0     0
     0     0     0     1     0     0     0     0     0     0    -1
     0     0     0     0    -1     0     0     1     0     0     0
     0     0     0     0     0    -1     0    -1     1     0     0
     0     0     0     0     0     0    -1     0    -1     0     0];
% 基本回路矩阵 Bf 顺时针为正
Bf = [
     0     1     0     0    -1     1     0    -1     0     0     0
     0     0     1     0     0    -1     1     0    -1     0     0
     1     1     1    -1     0     0     0     0     0     1    -1];
% 管段流量 q 单位 m^3/s 全部已知
q = [198.08 81.17 10.00 33.42 102.36 20.00 22.65 67.33 5.00 198.08 33.42]'/1000;
% 节点流量 Q 单位 m^3/s 全部已知
% 不需要求解
% 管段长度 L 单位 m
L = [320 650 550 270 330 350 360 590 490 0 0]';
% 管段直径 D 单位 m
D = [2^(2/5.333)*300 300 200 2^(2/5.333)*200 300 200 200 300 100 0 0]'/1000;
% 水泵静扬程 he 单位 m 书本 P85 的图 5.7 对管段 [10]、[11] 的水泵静扬程的标注多了负号
he = [48.96 0 0 0 0 0 0 0 0 12 48]';
% 水泵内阻系数 sp
sp = [138.5 0 0 0 0 0 0 0 0 0 0]';
% 水泵流量 qp 单位 m^3/s
% 等于管段流量 q

% 设计水头 H 单位 m
syms H [size(A,1) 1] real; 
H(1) = sym(0);
H(6) = sym(48.5);
% 节点地面标高 H_att 单位 m 
H_att = [0 13.60 18.80 19.10 22.00 32.20 18.30 17.30 17.50]';

% 曼宁粗糙系数 nM
nM = 0.012;
% 与水头损失计算指数公式相同的指数 n
n = 1.852;
% 环水头闭合差的最大允许值 eh 单位 m
eh = 0.01;

% 计算管段阻力系数 s
s = GetPipeResistanceCoefficient(L, D, sp, nM, 'Manning');
% 改进的哈代-克罗斯算法
[hf, hp, h] = ImprovedHardyClausProcess_LoopFlowCorrection(Bf, q, L, D, s, he, sp, nM, n, eh, 'Manning');
% 管段流速 v 单位 m/s
v = q./(pi*D.^2/4);
% 求解管段压降方程组
H = PressureDropEqsSolver(A, h, H);
% 计算自由水压
H_free = H - H_att;

7.管网节点方程组水力分析 牛顿-拉夫森算法求解节点校正压力方程组

% 书本 例5.4 管网节点方程组水力分析 牛顿-拉夫森算法求解节点校正压力方程组

% 1.
% 书本 P87 的式 (5.35) 的关联矩阵的第一行应写为 1 0 1 0 0 0 0

% 2.
% 排掉了之前的 bug，计算的时候发现，q 算的不对
% 如果 L, D 都没错的话，那就只能是 sp 的问题了
% 因此我才发现题目没说明 sp，就真的是全为零，这书也太坑了吧
    
% 3.
% 不知道为什么，书 P93 的管网节点压力平差计算过程数据表中，没有把管段 [9],[10],[11] 的数据展示出来
% 本来就很难比对我的计算结果对不对，他不写的话，直接就根本比对不聊了
% 更别说我还是把并联管道化简到一根了，难顶

clear;

% 关联矩阵 A 
A = [
     1     0     0     0     0     0     0     0     0
    -1     1     0     0     1     0     0     0     0
     0    -1     1     0     0     1     0     0     0
     0     0    -1    -1     0     0     1     0     0
     0     0     0     1     0     0     0     0     0
     0     0     0     0    -1     0     0     1     0
     0     0     0     0     0    -1     0    -1     1
     0     0     0     0     0     0    -1     0    -1];
% 节点流量 Q 单位 m^3/s 初值已知
Q = [-198.08 14.55 51.17 20.77 -33.42 35.03 82.33 27.65]'/1000;
% 管段长度 L 单位 m
L = [320 650 550 270 330 350 360 590 490]';
% 管段直径 D 单位 m
D = [2^(2/5.333)*300 300 200 2^(2/5.333)*200 300 200 200 300 100]'/1000;
% 水泵静扬程 he 单位 m 未知
he = [0 0 0 0 0 0 0 0 0]';
% 水泵内阻系数 sp 题目虽然没有给出，但是根据题目解答给出的值，应该是和 例5.3 相同的数值
sp = [0 0 0 0 0 0 0 0 0]';
% 水泵流量 qp 单位 m^3/s
% 等于管段流量 q

% 节点水头初值 H 单位 m
H = [46 43 40 39 40 42 39 37.5]';
% 要求解的节点水头的下标
H_var_subnum = [1 2 3 4 5 6 7]';
% 节点压力的校正公式
H_CorrectionFun = @(H,delta_H) H+0.5*delta_H;
% 节点地面标高 H_att 单位 m 
H_att = [13.60 18.80 19.10 22.00 32.20 18.30 17.30 17.50]';

% 海曾威廉系数 Cw
Cw = 110;
% 与水头损失计算指数公式相同的指数 n
n = 1.852;
% 节点流量闭合差的最大允许值 eq 单位 L/s
% 取太小了矩阵就奇异了，没办法，就取大一点吧
eq = 0.01;

% 牛顿-拉夫森算法
[q, H] = NewtonRaphsonProcess_NodalPressureCorrection(A, H, H_var_subnum, Q, L, D, sp, Cw, n, eq, 'n', H_CorrectionFun);
% 管段流速 v 单位 m/s
v = q./(pi*D.^2/4);
% 计算自由水压
H_free = H - H_att;

8.习题2

% 书本 习题5.2 

% 管网恒定流方程组
% 1.质量守恒方程 Σ(j∈Ωi)(qij)+Qi=0
% 2.能量守恒方程 Σ(j∈Λi)(hi,j)=0
% 3.物理关系方程 hi,j=si,j*qi,j^α
% 管段流量 qi,j
% 管段水头损失 hi,j
% 节点流量 Qi
% 与节点 i 关联的节点集合 Ωi
% 环 i 的管段集合 Λi
% 管段摩阻 si,j 当节点 i 与节点 j 不相关联，si,j -> ∞
% 指数 α

% 当水塔高度未知时，应给定水塔供水流量
% 即已知节点流量，该节点为定流节点，通过水力分析可以求解出节点压力，从而确定水塔高度
% 当水塔高度已经确定时
% 即已知该节点水头，该节点为定压节点，通过水力分析可以求解出节点流量，从而确定水塔供水量

% 1.
% 如上所述，如果水塔高度要确定节点水头的话，那就是地面标高为零才能让我这么干了
% 所以我把地面标高都设为了 0

% 关联矩阵 A 
A = [
     1     0     0
     0     1     0
     0     0     1
    -1    -1    -1];
% 管段流量 q 单位 m^3/s 全部未知
q = sym(zeros(size(A,2),1));
% 节点流量 Q 单位 m^3/s 全部未知
Q = sym(zeros(size(A,1),1));
% 管段长度 L 单位 m
L = [200 150 300]';
% 管段直径 D 单位 m
D = [100 100 150]'/1000;
% 水泵静扬程 he 单位 m
he = [0 0 0]';
% 水泵内阻系数 sp
sp = [0 0 0]';
% 水泵流量 qp 单位 m^3/s
% 等于管段流量 q

% 设计水头 H 单位 m
% 每题有所不同
% 节点地面标高 H_att 单位 m 
H_att = [0 0 0 0]';

% 海曾威廉系数 Cw
Cw = 100;
% 与水头损失计算指数公式相同的指数 n
n = 1.852;
% 节点流量闭合差的最大允许值 eq 单位 L/s
eq = 0.01;

% 第一问

% 设计水头 H 单位 m
H = [25.20 24.70 28.30 20]';
% 求解管段压降方程组
h = A'*H;
% 计算管段阻力系数 s
s = GetPipeResistanceCoefficient(L, D, sp, Cw);
% 根据 qj=(hj/sj)^(1/n),(j=1,2,...,M) he 都为零
q = (h./s).^(1/n);
% 求解节点流量方程组
Q = -A*q;

% 第二问

% 水塔供水，说明用水点水头低于水塔水头
% 水塔进水，说明用水点水头高于水塔水头
% 水塔不供水也不进水，说明用水点水头与水塔水头相等

% 综上，第二问情境下，用水点与水塔水头中位数相等

% 设计水头 H 单位 m
H = [25.20 24.70 28.30 25.20]';
% 求解管段压降方程组
h = A'*H;
% 计算管段阻力系数 s
s = GetPipeResistanceCoefficient(L, D, sp, Cw);
% 根据 qj=(hj/sj)^(1/n),(j=1,2,...,M) he 都为零
q = (h./s).^(1/n);
% 求解节点流量方程组
Q = -A*q;

% 第三问

% 设计水头 H 初值 单位 m
H = [25.20 24.70 28.30 20]';
% 水头变量下标记录
H_var_subnum = [4]';
% 设置 Q 初值 单位 L/s
Q_initial = [-16 -16 -18 50]';
% do...while 型
while 1
    % 设置 Q 初值 单位 m^3/s
    Q = Q_initial/1000;
    % 牛顿-拉夫森算法
    [q, H] = NewtonRaphsonProcess_NodalPressureCorrection(A, H, H_var_subnum, Q, L, D, sp, Cw, n, eq);
    % 求解节点流量方程组
    Q = -A*q;
    % 解得用水点的用水量少于 50 L/s，那么给初值增加一点
    if Q(4)-50/1000 < -(1e-5)
        Q_initial(4) =  Q_initial(4) + 1;
    % 解得用水点的用水量少于 50 L/s，那么给初值减少一点
    elseif Q(4)-50/1000 > 1e-5
        Q_initial(4) =  Q_initial(4) - 1;
    % 解得用水点的用水量与 50 L/s 之间的误差在允许范围内，退出迭代
    else
        break;
    end
end