一维圣维南方程差分数值算法中稀疏矩阵求解方法比较及优选研究

王浩骅，管光华*，肖昌诚,2

·区域农业水管理·

一维圣维南方程差分数值算法中稀疏矩阵求解方法比较及优选研究

王浩骅1，管光华1*，肖昌诚1,2

（1.武汉大学 水资源与水电工程科学国家重点实验室，武汉 430072；2.中国电建集团 华东勘测设计研究院有限公司，杭州 310014）

【】寻找高效、稳定的大型稀疏线性方程组求解算法以提高圣维南方程组求解速度。归纳了4种基于四点偏心格式的圣维南方程组求解算法并加以改进，并通过仿真试验，对比了不同算法的计算效率。计算断面数较少时（小于500），所有方法的运算时间基本一致；当计算断面数较大时（大于500），4种算法的速度较传统算法有了一定的提高，在计算断面数为1 520时，4种算法的计算速度是传统算法的4倍，在计算断面数为3040时，计算速度更是达到了10倍以上。改进后的高斯消元法和 PM 算法在计算断面数较大时计算速度较快，计算效率较高。该方法可应用于MPC控制、LQR控制等渠系自动化控制技术中，以提高仿真程序运算速度。

明渠一维非恒定流；圣维南方程组；大型稀疏矩阵；单个断面运算时间

0 引 言

【研究意义】为缓解水资源短缺，我国建立了许多大型调水工程，由于其调水距离长、输送水量大、沿线过水建筑物众多，控制和调度过程十分复杂，渠道在运行调度过程中不可避免地会出现非恒定流，而圣维南方程组是用来描述和求解非恒定流的重要途径。自圣维南方程组提出以来，通过长期的研究和实践已趋于完善[1]。随着大型调水工程的建设以及运行调度和控制过程的复杂化，如南水北调工程[2]、东深供水工程[3]、引黄济青工程[4]等，明渠一维非恒定流的求解过程要求更高，原有圣维南方程组的求解方法已经满足不了计算量和计算速度的要求。提高圣维南方程组的计算速度，有利于推进渠道运行调控的发展和推广圣维南方程组在其他工程领域的应用。【研究进展】目前对圣维南方程组求解方法的探讨较多，国内方面，陈大宏等[5]提出的DORA 算法、王泗远等[6]提出的不同浅水理论的数值计算方法为圣维南方程组的求解提供了不同的思路，而李文丰等[7]对具线性耗散摩阻力的 Cauchy 问题进的研究、聂大勇等[8]对具张弛项的圣维南方程组的弱解的研究丰富了圣维南方程组的应用范围。国外方面，Isaacson等[9]率先提出了显式和隐式的离散方法，为Preissmann 提出其经典的四点差分隐格式创造了前提条件。此后，Fread等[10]提出的压缩存贮消元法、Liggett等[11]提出的双追赶法使圣维南方程组的求解方法更加丰富。在求解大型稀疏矩阵的算法优化方面，郑经纬等[12]提出的基于CUDA的PCG算法大大提升了稀疏矩阵的计算效率，而潘少华等[13]为探究低秩稀疏矩阵优化问题提出了凸松弛和因子分解法，为低秩稀疏矩阵的求解提供了一种思路。目前，稀疏矩阵的求解存在的主要问题有：①如何平衡好大型稀疏矩阵的求解速度与其计算结果的收敛性、精度仍是稀疏矩阵求解的主要问题，即加快求解速度的同时保证必要的收敛和精度，使结果可用；②大型稀疏矩阵元素较少，其对应的方程组易呈病态，现有的矩阵压缩存储效率有待提高，预处理方法和优化模型算法有待改进；③作为大型复杂稀疏矩阵，圣维南方程组的求解时间往往过长，如何利用并行计算技术进行圣维南方程组的算法研究以提高计算速度显得愈发重要。【切入点】圣维南方程组的求解速度和求解算法仍需改进，提升求解速度。【拟解决的关键问题】本文为寻找高效、稳定的大型稀疏线性方程组的求解算法以提高圣维南方程组求解速度，归纳了并优化了4种基于四点偏心格式的圣维南方程组求解算法，分别是追赶法（chase）、循环递减算法（CR）、高斯消元法（GE）以及划分算法（PM）。同时，还比较了不同求解算法的计算量，并通过仿真实验，对比了不同算法的计算效率。

1 明渠一维非恒定流及其数值解法

1.1 明渠一维非恒定流模型

圣维南方程组由反映质量守恒定律的连续性方程和反映动量守恒定律的运动方程组成。

1）连续性方程

根据质量守恒，流体系统的质量随时间的变化率为零，并运用输出方程，并考虑旁侧入流，则一维明渠非恒定渐变流的连续性微分方程[14]为：

式中：为水面宽（m）；为水位（m）；为时间（s）；为流量（m3/s）；为水深（m）；为断面的距离坐标（m）；为旁侧入流量（m3/(s·m)）。

2）运动方程

对于明渠的非恒定渐变流，通常取在水面，水面相对压强为0，考虑到流量和流速的关系，一维明渠非恒定渐变流的运动微分方程[14]为：

式中：为水面宽（m）；为水位（m）；为时间（s）；为流量（m3/s）；为断面的距离坐标（m）；g为重力加速度（m/s2）；为过水断面面积（m2）；为水流沿轴线方向的流速（m/s）。

1.2 有限差分法

有限差分法是求解圣维南方程组常用的数值计算方法之一。有限差分法又可分为显式差分法和隐式差分法，隐式差分法的基本思想是直接求解由内断面方程和边界方程联立组成的方程组。在所有隐式差分算法中，应用较多的是四点偏心格式。

四点偏心格式，也称Preissmann差分格式[15-16]，是针对矩形网格中间某点M将因变量的微分形式转化为差分形式。在距离步长上M点处于正中心，而在时间步长上，存在一个权重因子，偏向已知时层为，偏向未知时层为1-。

U、D、L、R 4 点的流量可以通过线性差值得到，由此可以得到 M 点流量和水位的差分形式，并代入圣维南方程组，可以得到圣维南方程组的差分方程[14]：

式中：=1,2,……,-1。

在实际计算过程中，若将全河段划分为个断面，则有-1个河段，每个断面有、共2个未知数，共2个未知数。而每个河段可建立2个方程，-1个河段一共可以建立2-2个方程。因此还需要根据上下游边界的实际情况补充2个方程，构成2个方程才能求解出所有未知数。上下游边界的条件可统一写为：

上述的方程组构成了大型稀疏矩阵。该方程组的矩阵形式为：

2 稀疏矩阵求解方法

2.1 基于四点偏心格式的算法优化

圣维南方程组的求解需要确定边界条件，而边界条件的存在可能使矩阵的某些主对角元素为0，这使得算法求解过程中会出错，因此需要进行处理。处理的方式有3种：

1)通过初等变换使该位置的元素不为0，这会增加算法的计算量，对于边界条件较少的矩阵可采用这种方法。

2)将为0的主对角元素用一个极小值来替换，这会引入一定的误差，其误差的大小可以通过极小值的大小进行调整。

3)将稀疏矩阵划分为多个子矩阵，然后通过初等变换进行消元处理，再采用并行技术求解。适用于三对角线性矩阵。

图1 稀疏矩阵算法优化流程

2.1.1 追赶法（chase）

如式（7）所示，四点偏心格式需要求解的系数矩阵是五对角系数矩阵的一种，首先要通过矩阵间的初等变换将其转化为三对角线性方程，再进行追赶法的操作，能简化很多步骤。本文采用的算例的边界条件较少，因此选择第一种方式消除主对角的0元素。

2.1.2 循环递减算法（CR算法）

该算法[17]主要分为2个步骤：向前约化和向后替代。前一步骤为消元，即通过矩阵变化消去所有下标为奇数的变量，保留下标为偶数的变量。此时，大型线性方程组的规模减半，并且消元之后构成新的线性方向组，其系数矩阵仍为原对角形式。按照这个方法执行下去，逐步消元递减，直到最后只剩下含2个变量的2个方程组时停止消元。对该二元线性方程组进行直接求解，得到2个未知量的解。后一步骤为回代求解，将求得的未知量带到消元过程的上一个方程组，求出上一方程组的未知量，如此往复求解计算，直到求解得到整个线性方程组的解。

与追赶法求解四点偏心格式相似，CR算法计算前需要将五对角线性方程组转化为三对角线性方程组。主对角0元素的存在也会影响CR算法的计算过程。通过初等变换使主对角元素为0项进行调整的方法并不适合本算法，因此，本文将为0的主对角元素用一个极小值来替换。

2.1.3 高斯消元法(GE)

以三对角线性方程组为例，其求解思路[18]可分为3个部分：首先是消去系数矩阵对角线以下的元素。将第2个方程乘以2(1)1得到的结果加到第2个方程上，这样，第2个方程就消去了系数2。同理，之后的每一个方程逐步消元递减，直到消去对角线以下的所有元素，接着消去系数矩阵对角线以上的元素。如式（8）所示，将第个方程乘以-b-1(a)-1得到的结果加到第-1 个方程上，这样，第-1个方程就消去了系数b-1。同理，之后的每一个方程逐步消元递减，直到消去对角线以下的所有元素，线性方程组变成：

最后求解这一线性方程组，就可以得到原线性方程组的解。与CR算法求解四点偏心格式相似，高斯消元法将为0的主对角元素用一个极小值来替换。

2.1.4 划分算法（PM算法）

对于矩阵系数为大型稀疏矩阵的方程组，可根据计算机并行的核数，进一步分解为周期块三对角线性方程组。假设计算机的核数为，可以得到下列方程组：

式中：D为m阶方阵，I为m阶单位方阵，H和G分别为m×m+1和m×m-1矩阵。

根据式（12）可以得到：

Preissmann 差分格式需要求解的大型稀疏矩阵属于五对角矩阵，其每一行或每一列最大非0元素的个数为4，属于元素较少的大型稀疏矩阵。以下是根据该类型矩阵的算法改进。

若B矩阵只有第一列元素非0的情况，则H也只有第一列元素非 0；同理，若C矩阵只有最后一列元素非0的情况，则G也只有最后一列元素非 0。以下就以这种情况为研究对象。

由于H只有第一列元素非0，G只有最后一列元素非2，式(13)可转化为：

四点偏心格式需要求解的大型稀疏矩阵，可以看作是一种周期块三对角线性方程组。通过矩阵分块发现，不同的计算机核数分解产生的B和C矩阵格式也有所差异，可能出现以下2种情况（以B矩阵为例）：

①其分解所得的矩阵H仍只有第一列元素非0，可直接求解；

②其分解所得的矩阵H第一列和第二元素非0，因此在进行矩阵分解前，需要进行矩阵行变换，消去部分项，将B矩阵第二列元素消去。

同理，C矩阵也按上述格式进行矩阵行变换。

2.2 对比分析

通过对四点偏心格式的不同求解算法的分析，可以得到不同算法的计算量与其系数矩阵阶次的关系，其中1/2。在计算机CPU运算过程中，乘法的计算速度比加法慢近10倍。当然，在不同的编译环境里采用不同的算法，计算速度是不同的[20]，比如C++中，加法一般比乘法快1倍，而刘东[21]提出的Booth 算法的乘法器可以大幅提高乘法的运算速度。因此本文主要比较算法的乘除法计算量。如表1和图2所示，随着矩阵阶次的增加，传统算法随阶次的平方增长，而4种算法的计算量呈线性增长，4种算法的计算量较传统算法大大减少了。当计算阶次取值较大时，4种算法中计算量较小的是循环递减法算法，较大的是划分算法。

每种算法的运算量的分析是基于理论基础的，而实际计算过程中，其数据传输过程和计算阶次也会对运算时间产生较大的影响。因此确定每种算法的计算速度还需具体实验进行验证。

表1 5种算法的计算量

Table 1 The amount of calculation of five algorithms

3 算例验证及结果分析

3.1 仿真渠段

漳河灌区以漳河水库为主要水源，是典型的“长藤结瓜”灌溉系统。灌区渠系分支众多，在整个灌区面积上形成了庞大复杂的灌溉渠网系统。灌区内有较多的渠道及渠系水工建筑物，是我国典型的综合型灌区。探索如何在这种庞大复杂的渠网系统中进行高效的运行调度显得至关重要，这也是本文仿真的目的之一。本文选取湖北省漳河灌区三干渠三分干部分渠段2004年8月21日的数据进行仿真研究。该渠段中渠系建筑物众多，包括大量的农耕桥、涵闸以及取水口等，为了仿真和计算方便，在建模过程中做了相应的简化：为这些渠系附属建筑物设置相应的水头损失系数，将其视为对应节点的水头损失。简化后渠段的几何参数和水力特性见图3和表2。其中渠首闸的闸门参数为：闸门宽度4m，闸门流量系数0.83，闸门最大开度3m。

仿真的边界条件和初始条件为：①边界条件：给定上游、下游流量边界。②初始条件：等体积运行，通过设计水深求各渠段体积，反算各渠段初始下游水深，得到初始水面线。

图3 仿真渠段示意

表2 仿真渠池的几何参数及水力特性

Table 2 Geometric parameters and hydraulic characteristics of simulation canals

3.2 算例设计

为了说明以上多种大型稀疏矩阵求解方法的计算速度，设计以下试验：

1）运算平台：Windows7 系统 Microsoft Visual C++ 2013 (C)。

2）运算环境：Matlab 程序包[22]。

3）计算断面数：分别为190、380、760、1 520、3 040、4 560、6 080。

4）渠道仿真参数：仿真时间48h；时间步长5min，空间步长随计算断面数变化。

对于同一渠道的非恒定流分析，计算精度会随着选取计算断面数的增加而增加；但并不是计算断面数越多越好，实际上计算断面数满足计算目的要求即可。在边际效应的作用下，即便断面数进一步增加，其计算精度并无实质性提高。本文为了简化算例，仅在计算对象不变的情况下增加断面数目即为考虑实际工程的规模可能越发巨大。

每次计算过程为整个渠道的非恒定流仿真过程，因此其运算时间不仅包括大型稀疏矩阵的求解，也包括其他部分的算法和数据处理。

为提高实验数据的准确度，在运行过程中关闭其他后台程序。同时，每种工况应执行3次，取3次实验平均值作为该工况的实验数据。

3.3 仿真计算及结果分析

3.3.1 运算时间

通过上述工况进行实验，得到5种算法在不同工况下的运算时间，如表3及图4所示。

表3 5 种算法在不同计算断面下的运算时间

注 计算断面数大于3 040时，采用传统算法计算时内存不足，故无数据。

图4 5 种算法运算时间与运算断面数的关系图

3.3.2 单个断面运算时间

为了更直观地反映5种算法的计算效率，即运算时间与计算断面数的关系，本文引入单个断面运算时间单：

式中：单为单个断面运算时间（s）；总为总运算时间（s）；计算断面数。

由图4及图5可知，计算断面数较少时（小于500），所有方法的运算时间基本一致，说明求解大型稀疏矩阵的时间占整个非恒定流问题的很小一部分，因此矩阵求解算法的优劣无法体现。当计算断面数较大时（大于500），4种算法的速度较传统算法有了一定的提高，在计算断面数为1 520时，4 种算法的计算速度是传统算法的4倍，在计算断面数为3 040时，计算速度更是达到了10倍以上。传统算法的单随断面数的增加而呈指数型增长，而其他4种算法的单随着计算断面的增加而减少，然后趋于稳定；其原因在于传统算法计算过程中未忽略0元素的计算，随着计算断面数的增大，其计算量呈指数型增长；其他4种算法忽略了大量0元素，只计算非0元素，计算量大大减少。

图 5 5 种算法运算的单个断面运算时间与运算断面的关系

3.3.3 误差分析

本文将4种算法的结果与传统算法进行误差对比。为了进行定量比较，本文选取渠道非恒定流计算过程的各个断面的最终水位作为比较值，分别计算4种算法与传统方法的水位之间的水位差值，并取其平均值作为绝对误差进行分析。

由表4和图7可知，4种算法的误差都较小，数量级均在-5以下，当断面较大时，PM算法和GE算法计算速度更快，误差更小。其中算法误差最小的是PM算法，算法误差最大的是CR算法。这是因为PM算法是通过划分并行来求解，相对其他3种算法而言优化处理过程不存在误差，所以误差最小。而CR算法不仅采用了“将为0的主对角元素用一个极小值来替换”的优化处理方法，而且其算法的步骤为消元加迭代的过程，会使误差反复累积，因此误差最大。综上所述，误差的结果较为合理。

表 4 4种算法的水位绝对误差平均值

Table 4 Mean absolute water level error of four algorithms m

4 讨 论

国内外学者在大型稀疏矩阵问题上做了大量的研究工作，并取得了较好的效果。然而，由于大型稀疏矩阵的复杂性，始终无法平衡好大型稀疏矩阵求解中求解速度与收敛性和精度的矛盾[15]，也无法彻底解决由于主对角0元素导致的病态方程组问题[16]。

基于此，本文在Preissmann 差分格式研究的基础上，提出了4种优化算法并进行了仿真验证。结果表明：①5种算法单个断面运算时间的差异体现了大型稀疏矩阵中0元素对计算量的巨大影响。通过正确的方式处理0元素，可以使计算量呈指数级下降。当计算断面数为1 520时，4种算法的计算速度是传统算法的4倍；当计算断面数为3 040时，新算法与传统算法的计算速度差距达到了10倍以上。②不同算法的优化处理方式差异是导致了其绝对误差产生差异的内因。PM算法是通过划分并行来求解，相对其他3种算法而言优化处理过程不存在误差，所以误差最小。而CR算法不仅采用了“将为0的主对角元素用一个极小值来替换”的优化处理方法，而且其算法的步骤为消元加迭代的过程，会使误差反复累积，因此误差最大，但其误差均不超过-5数量级。

未来在对圣维南方程组的求解研究应加强以下几个方面：①在优化圣维南方程组的算法过程中，应结合实际问题来考虑合理的误差允许值。目前的许多求解算法没有结合实际的问题背景，往往导致其无法很好的运用于解决实际问题。②建立大量的模型试验来验证仿真结果。圣维南方程组的求解目的是为了解决实际的水力学问题，仿真结果只有通过大量的模型试验验证，才有实际价值。③在圣维南方程组的求解中引入并行计算技术。目前已有很多关于并行计算技术在大型稀疏矩阵中应用的研究，但并行计算技术在圣维南方程组中的应用的相关研究仍较少。并行计算能够有效加快运算速度，提高运算能力，求解更大规模的问题。但其同时有其局限性：求解的问题具有一定的并行度且需要合理的并行算法。因此，并行计算技术在圣维南方程组中的应用也需要进一步研究。

本文的结果对提高大型渠道非恒定流仿真的运算速度具有参考价值，其方法可应用于MPC控制、LQR控制等渠系自动化控制技术中，以提高仿真程序运算速度。

5 结 论

1)当计算断面较多时，4种算法的计算速度优势得以体现，高斯消元法和PM算法的计算速度优势相对更大，且计算速度的差距随着断面数的增加不断扩大。

2)当计算断面数小于500左右时，5 种算法的计算速度基本一致；当断面数在500到1 500范围内，可选择除传统算法以外的其他4种算法进行计算；当断面数大于1 500时，选择高斯消元法和PM算法为佳。相较高斯消元法，PM算法编程复杂，但易实现并行计算。

3）4种算法的绝对误差都较小，最大绝对误差平均值也都不超过-5数量级，完全在可接受范围之内，并且其误差值并不会随着计算断面的增加而增大。

4）在计算断面较大时，推荐采用PM算法与高斯消元法；当考虑并行计算时，推荐采用PM算法。

[1] 聂大勇.一维圣维南方程组整体经典解[D].郑州:华北水利水电学院,2007.

NIE Dayong. Global Classical Solution of One-dimensionalSaint-Venant Equations [D]. Zhengzhou: North China University ofWater Resources and Electric Power, 2007.

[2] 黄会勇, 闫弈博, 高汉,等. 南水北调中线总干渠运行调度反馈控制方式研究[J]. 人民长江, 2014, 45(6): 56-59.

HUANG Huiyong, YAN Yibo, GAO Han, et al. Study of feedback control on operation and dispatch of main canal of middle route project of south-to-north water diversion[J]. Yangtze River, 2014, 45(6): 56-59.

[3] 曾庚运. 东深供水改造工程全线供水调度与管理自动化[J]. 中国农村水利水电, 2003(9): 1-5.

ZENG Gengyun. All-line water supply dispatching and the management automation in Dongshen water supply improvement project[J]. China Rural Water and Hydropower, 2003(9): 1-5.

[4] 韩延成, 赵洪丽, 张烨,等. 引黄济青工程调度运行自动控制系统研究[A]. 宫崇楠, 张金丽. 山东水利学会第十届优秀学术论文集[C]. 济南: 山东省科学技术协会, 2005.

HAN Yancheng, ZHAO Hongli, ZHANG Ye, et al. Research on the automatic control system for dispatching operation of the Yellow River Diversion Project[A]. Gong Chongnan, Zhang Jinli. The 10th Excellent Academic Papers of Shandong Water Conservancy Society[C]. Jinan: Shandong Association for Science and Technology, 2005.

[5] 陈大宏, 蓝霄峰, 杨小亭. 求解圣维南方程组的DORA算法[J]. 武汉大学学报(工学版), 2005, 38(5): 41-44.

CHEN Dahong, LAN Xiaofeng, YANG Xiaoting. DORA approach for solution of the Saint-Venant equations[J]. Engineering Journal of Wuhan University, 2005, 38(5): 41-44.

[6] 王泗远, 郭泽庆, 王高亚. 明渠非恒定流问题的数值计算方法: 以圣• 维南方程组的数值计算方法为例[J]. 中国新技术新产品, 2011 (10): 7.

WANG Siyuan, GUO Zeqing, WANG Gaoya. Numerical calculation method of non-constant flow in open channel: Taking the numerical calculation method of Saint Vinnan equations as an example [J]. China New Technologies and Products, 2011 (10): 7.

[7] 李文丰, 聂大勇. 具线性耗散摩阻力的圣维南方程组Cauchy问题探讨[J].黄河水利职业技术学院学报,2013, 25(1): 51-54.

LI Wenfeng, NIE Dayong. Cauchy Problem of Saint Venan Equations with Linear Dissipative Friction[J]. Journal of Yellow River Conservancy Technical Institute, 2013, 25 (1): 51-54.

[8] 聂大勇, 高杰, 陈征. 具张弛项的圣维南方程组的弱解[J]. 应用数学, 2016, 29(2): 381-387.

NIE Dayong, GAO Jie, CHEN Zheng. Weak solutions for saint-venant systems with relaxation[J]. Mathematica Applicata, 2016, 29(2): 381-387.

[9] ISAACSON E, STOKER J J, TROESCH A. Numerical Solution of Flood Prediction and River Regulation Problems[M]. New York: New York University Institute of Mathematical Science, 1956.

[10] FREAD D L, SHEDD R C, SMITH G F, et al. Modernization in the national weather service river and flood program[J]. Weather and Forecasting, 1995, 10(3): 477-484.

[11] LIGGETT J, CUNGE J. Numerical methods of solution of the unsteady flow equations[M]. Fort Collins, Colorado: Water Resources Publications, 1975.

[12] 郑经纬, 安雪晖, 黄绵松. 基于 CUDA 的大规模稀疏矩阵的 PCG 算法优化[J]. 清华大学学报(自然科学版), 2014, 54(8): 1 006-1 012.

ZHENG Jingwei, AN Xuehui, HUANG Miansong. CUDA-based PCG algorithm optimization for a large sparse matrix[J]. Journal of Tsinghua University (Science and Technology), 2014, 54(8): 1 006-1 012.

[13] 潘少华, 文再文. 低秩稀疏矩阵优化问题的模型与算法[J]. 运筹学学报, 2020, 24(3): 1-26.

PAN Shaohua, WEN Zaiwen. Models and algorithms for low-rank and sparse matrix optimization problems[J]. Operations Research Transactions, 2020, 24(3): 1-26.

[14] 赵昕, 张晓元, 赵明登. 水力学[M]. 北京: 中国电力出版社, 2009.

ZHAO Xin, ZHANG Xiaoyuan, ZHAO Mingdeng. Hydraulics[M]. Beijing: China Electric Power Press, 2009.

[15] 顾峰峰, 倪汉根. 四点时空偏心隐格式的改进求解[J]. 大连理工大学学报, 2007, 47(3): 419-423.

GU Fengfeng, NI Hangen. Two improved calculation methods of Preissmann four-point linear implicit scheme[J]. Journal of Dalian University of Technology, 2007, 47(3): 419-423.

[16] 葛华. Preissmann 四点时空偏心隐格式思想在二维数学模型中应用初探[J]. 科学技术与工程, 2007, 7(1): 91-93.

GE Hua. Application of the preissmann difference scheme to the 2-D model (preliminary study)[J]. Science Technology and Engineering, 2007, 7(1): 91-93.

[17] 唐光平. 基于三对角线性方程组的混合并行算法研究[D]. 长沙: 湖南大学, 2015, 2015.

TANG Guangping. Research on hybrid parallel algorithms based on tridiagonal linear equations [D]. Changsha: Hunan University, 2015.

[18] 彭朝英. 高斯消元法的改进及其在工程上的应用[J]. 邵阳学院学报(自然科学版), 2011, 8(2): 26-30.

PENG Chaoying. Improvement of the Gaussian elimination method and its application in engineering MCU[J]. Journal of Shaoyang University (Science and Technology), 2011, 8(2): 26-30.

[19] 朱成. 五对角线性方程组的并行求解算法的研究[D]. 长沙: 湖南大学, 2016.

ZHU Cheng. Research on Parallel Solving Algorithm of Pentadiagonal Linear Equations [D]. Changsha: Hunan University, 2016.

[20] 罗福强, 冯裕忠, 茹鹏. 计算机组成原理[M]. 北京: 清华大学出版社, 2011.

LUO Fuqiang, FENG Yuzhong, RU Peng. Principles of computer organization[M]. Beijing: Tsinghua University Press, 2011.

[21] 刘东. 采用 Booth 算法的16×16 并行乘法器设计[J]. 现代电子技术, 2003, 26(9): 21-22, 25.

LIU Dong. Design of parallel multiplier for borth algorithm[J]. Modern Electronics Technique, 2003, 26(9): 21-22, 25.

[22] 武汉大学. 输水渠道系统运行仿真与控制软件 V1.0: 中国, 2011SR034392[P]. 2011-06-03.

Wuhan University. Simulation and control software for water conveyance system operation V1.0: China, 2011SR034392[P]. 2011-06-03.

[23] BONDARENCO M, GAMAZO P, EZZATTI P. A comparison of various schemes for solving the transport equation in many-core platforms[J]. The Journal of Supercomputing, 2017, 73(1): 469-481.

[24] NAVABPOUR B, OSTAD ALI ASKARI K, ESLAMIAN S, et al. Comparison of solutions of Saint-Venant equations by characteristics and finite difference methods for unsteady flow analysis in open channel[J]. International Journal of Hydrology Science and Technology, 2018, 8(3): 229.

[25] 黄凯, 管光华, 刘大志,等. 串联渠系 PID 改进积分与微分环节仿真研究[J]. 灌溉排水学报, 2017, 36(2): 1-11.

HUANG Kai, GUAN Guanghua, LIU Dazhi, et al. Simulation for PID controller of modified integral and differential on tandem canal system[J]. Journal of Irrigation and Drainage, 2017, 36(2): 1-11.

[26] 李一鸣, 管光华, 陈琛,等. 渠道糙率影响因素分析及预测模型研究[J]. 灌溉排水学报, 2017, 36(S2): 155-161, 189.

LI Yiming, GUAN Guanghua, CHEN Chen, et al. Analysis of influencing factors of canal roughness and prediction model [J] .Journal of Irrigation and Drainage, 2017, 36 (S2): 155-161, 189.

[27] 张慧, 伍萍辉, 张馨,等. 非垂直拍摄获取细叶作物覆盖度优化算法研究[J]. 灌溉排水学报, 2019, 38(9): 55-62.

ZHANG Hui, WU Pinghui, ZHANG Xin, et al. Research on optimization algorithm of fine leaf crop coverage based on non-vertical shooting[J]. Journal of Irrigation and Drainage, 2019, 38(9): 55-62.

Comparison and Optimization of Sparse Matrix Solution Methods in One-dimensional Saint-venant Equation Difference Numerical Algorithm

WANG Haohua1, GUAN Guanghua1, XIAO Changcheng1, 2

(1. State Key Laboratory of Water Resources and Hydropower Engineering Science, Wuhan University, Wuhan 430072, China; 2. Power China HuaDong Engineering Corporation, Hangzhou 310014, China)

【】In order to alleviate the shortage of water resources, China has established many water transfer projects. Due to its long water transfer distance, large water delivery volume, and numerous water passing buildings along the line, the control process is very complicated. Unsteady flow will inevitably appear in the channel operation scheduling process, and the Saint-venant equations are an important way to describe and solve the unsteady flow.【】With the construction of large-scale water transfer projects and the complexity of the operation scheduling and control process, the traditional method of solving the sparse matrix of the original Saint-Venant equations has been unable to meet the requirements of calculation volume and calculation speed. In order to find efficient and stable algorithms for solving large and sparse linear equations to improve the speed of solving Saint-Venant equations.【】In this paper, four algorithms for solving Saint-Venant equations based on the four-point eccentric scheme are summarized and improved, and the calculation efficiency of different algorithms is compared through simulation experiments.【】From the simulation results, it can be seen that when the number of calculation sections is small (less than 500), the calculation time of all methods is basically the same; when the number of calculation sections is large (more than 500), the speed of the four algorithms is improved compared with the traditional algorithm.When the number of cross sections is 1 520, the calculation speed of the four algorithms is 4 times that of the traditional algorithm. When the number of cross sections is3040, the calculation speed is more than 10 times.【】The improved GE and PM algorithms have faster calculation speed and higher calculation efficiency when the number of cross sections is larger. The results of this paper have reference value for improving the calculation speed of large channel non-constant flow simulation. The method can be applied to canal system automatic control technology such as MPC control, LQR control point, etc., to improve the running speed of simulation program.

one dimensional unsteady flow in open-channel; Saint-venant equations; large sparse matrices; calculation time of a single section

TV133

A

10.13522/j.cnki.ggps.2020182

1672-3317（2021）03-0116-09

王浩骅, 管光华, 肖昌诚. 一维圣维南方程差分数值算法中稀疏矩阵求解方法比较及优选研究[J]. 灌溉排水学报, 2021, 40(3): 116-124.

WANG Haohua, GUAN Guanghua, XIAO Changcheng. Comparison and Optimization of Sparse Matrix Solution Methods in One-dimensional Saint-venant Equation Difference Numerical Algorithm[J]. Journal of Irrigation and Drainage, 2021, 40(3): 116-124.

2020-03-31

国家自然科学基金项目（51979202，51439006，51009108）；“十三五”国家重点研发项目（2016YFC0401810）

王浩骅（1996-），男，浙江台州人。硕士研究生，主要从事灌排自动化研究。E-mail: 2290586701@qq.com

管光华（1979-），男，江苏阜宁人。副教授，硕士生导师，主要从事灌排自动化及量水理论研究。E-mail: GGH@whu.edu.cn

责任编辑：白芳芳