基于拟蒙特卡罗方法的供水管网抗震可靠性分析并行化研究

doi:10.3785/j.issn.1008-973X.2020.02.004

基于拟蒙特卡罗方法的供水管网抗震可靠性分析并行化研究

龙立^,, 郑山锁^,, 周炎, 贺金川, 孟宏立, 蔡永龙

Parallel study of seismic reliability analysis of water supply pipe network based on quasi-Monte Carlo method

LONG Li^,, ZHENG Shan-suo^,, ZHOU Yan, HE Jin-chuan, MENG Hong-li, CAI Yong-long

通讯作者: 郑山锁，男，教授. orcid.org/0000-0003-1171-4511. E-mail： zhengshansuo@263.net

收稿日期: 2019-07-30

Received: 2019-07-30

作者简介 About authors

龙立（1992—），男，博士生，从事城市生命线工程地震灾害损失评估.orcid.org/0000-0001-8429-5608.E-mail：longliforever@163.com , E-mail：longliforever@163.com

摘要

为了提高基于蒙特卡罗（Monte Carlo）方法的供水管网抗震可靠性分析效率，以低偏差Sobol点列替代伪随机数序列对供水管网节点和管段破坏概率进行抽样，结合宽度优先搜索算法，提出基于拟Monte Carlo方法和统一计算设备架构（CUDA）的供水管网抗震可靠性分析并行算法，并从内存、执行配置和指令等方面优化并行算法. 以某城市供水管网系统为例，对比串行和并行计算方法的精度及效率，分析Sobol点列和伪随机数序列对管网可靠性分析的影响. 结果表明，并行和串行方法计算结果的误差最大为0.52%，并行方法最高加速比为串行算法的96倍，在保证结果精度的同时大幅度提高计算效率. 基于Sobol点列进行1 000次并行模拟及基于伪随机数序列进行5 000次并行模拟，2种模拟结果与基于模糊数学法的解析值的最大误差分别为0.2%、0.4%，表明基于拟Monte Carlo的并行方法具有更高的精确度，更快的收敛速度.

关键词： 供水管网 ; 统一计算设备架构（CUDA） ; 宽度优先搜索 ; 并行计算 ; 网络可靠性分析 ; 拟Monte Carlo方法

Abstract

In order to improve the seismic reliability analysis efficiency of water supply pipe network based on Monte Carlo simulation, the failure probabilities of water supply pipe network nodes and pipes were sampled by using low discrepancy Sobol sequence instead of pseudo-random number sequence. Combined with the breadth-first search algorithm, a parallel algorithm for seismic reliability analysis of water supply pipe network based on quasi-Monte Carlo method and compute unified device architecture (CUDA) was proposed. The parallel algorithm was optimized from the aspects of memory, execution configuration and instructions. A city water supply pipe network was taken as the computational example, the accuracy and efficiency of serial and parallel computing methods were compared, and the influence of Sobol sequence and pseudo-random number sequence on the reliability analysis of pipe network was analyzed. Results show that the maximum error of the parallel and serial methods is 0.52%. The maximum acceleration ratio of the parallel method is 96 times that of the serial method, and the parallel method significantly improves the computational efficiency while ensuring the accuracy of results. 1 000 parallel simulations were performed based on Sobol sequences and 5 000 parallel simulations were performed based on pseudo-random number sequences, and the maximum errors between the two simulation results and the analytical value based on fuzzy mathematics were 0.2% and 0.4%, respectively. It indicates that the parallel method based on quasi-Monte Carlo has higher accuracy and faster convergence speed.

Keywords： water supply pipe network ; compute unified device architecture (CUDA) ; breadth-first search ; parallel computing ; network reliability analysis ; quasi-Monte Carlo method

PDF (795KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

龙立, 郑山锁, 周炎, 贺金川, 孟宏立, 蔡永龙. 基于拟蒙特卡罗方法的供水管网抗震可靠性分析并行化研究. 浙江大学学报(工学版)[J], 2020, 54(2): 241-247 doi:10.3785/j.issn.1008-973X.2020.02.004

LONG Li, ZHENG Shan-suo, ZHOU Yan, HE Jin-chuan, MENG Hong-li, CAI Yong-long. Parallel study of seismic reliability analysis of water supply pipe network based on quasi-Monte Carlo method. Journal of Zhejiang University(Engineering Science)[J], 2020, 54(2): 241-247 doi:10.3785/j.issn.1008-973X.2020.02.004

供水管网系统是生命线工程中不可或缺的基础设施，在近年来多次破坏性地震中，均遭到严重破坏，致使城市供水功能失效，严重影响人们的生活^[1-4]. 因此，有必要进行供水管网抗震可靠性分析，为供水管网震前加固和震后修复提供科学支撑. 目前，在供水管网抗震可靠性分析方面，国内外学者采用的主流方法有解析法、代理法和模拟法^[5]. 其中，模拟法分为蒙特卡罗（Monte Carlo）法^[6]、拟Monte Carlo法^[7-8]和重要度抽样法^[9]等. 相比于解析法和代理法，模拟法的优点是能考虑可靠性评估过程中多种不确定性，对管网的各种工况进行模拟，并回避网络可靠度分析中多项式复杂程度的非确定性问题，对问题的维数不敏感，无须考虑网络拓扑的复杂性，建模方便^{[5, 10]}. 其缺点也很明显，须多次运算才能获得较为精确和稳定的结果，往往会花费较长的计算时间.

近年来，随着图形处理单元（graphic processing units，GPU）技术的飞速发展及统一计算设备架构（compute unified device architecture，CUDA）的推出，基于CUDA的GPU并行计算技术已经获得广泛应用. 由于模拟法具有天然的并行性，每次模拟产生的结果相互独立，采用并行计算较适合. 季经纬等^[11]研究热辐射Monte Carlo法在CUDA平台上的应用，采用不同计算能力的GPU进行计算，运算时间相比单个中央处理器（central processing unit，CPU）的运算时间最高缩短33倍. 严立等^[12]基于CUDA Fortran语言对直接Monte Carlo法进行并行优化，在保证计算精度的情况下，有效提升了计算速度. 国外学者^[13-16]也将基于CUDA的并行模拟法应用到各自的研究领域，并取得了较好的加速效果. 目前，在供水管网抗震可靠性分析领域的应用还相对较少.

本研究以城市管网系统为研究对象，基于拟Monte Carlo方法，研究供水管网抗震可靠性分析在CUDA平台上的并行算法，并从内存、指令和执行配置方面对并行算法进行优化. 以某供水管网为例，对比串行和并行方法计算结果的精度及计算效率，并对Monte Carlo法与拟Monte Carlo法的收敛效率进行对比.

1. 基于拟Monte Carlo法的可靠性分析方法

供水管网系统是由供水源点、汇点及不同类型的节点单元和管线单元组成的复杂网络系统，本质上可简化为由点和边组成的网络拓扑结构，其数学模型可用“图”表示. “图”是由若干给定的顶点及连接两顶点的边构成的图形，常用来描述事物间的某种特定关系，其顶点用于代表事物，连接顶点的边则用于表示相应事物间具有的关系. 为了便于对供水管网进行抗震可靠性分析，须进行如下假设：1）只考虑水厂及泵站节点的破坏，其他节点及汇点不失效；2）节点及管段在地震作用下只有完好和破坏2种状态，分别用1和0表示，且节点与管段的状态是相互独立的；3）水的流通是单向的，只能从水源点、水厂、泵站等节点流到其他普通节点.

基于上述理论及假设，假设供水系统有 $m$个节点，其抗震可靠性分析可按以下步骤进行.

1）分别计算供水系统中各节点及管段在地震作用下的破坏概率 ${P_i}$、 ${P_{ij}}$，其中普通节点不失效，其破坏概率为0.

2）在各节点及管段上分别产生0~1之间的随机数 ${\delta _i}$、 ${R_{ij}}$，并与对应节点及管段的破坏概率进行比较，若 ${\delta _i} \geqslant {P_i}$，认为节点正常工作，反之认为节点失效；若 ${R_{ij}} \geqslant {P_{ij}}$，认为管段正常工作，反之认为管段失效. 从而产生关于节点及管段的2个矩阵N=[a_i]和L=[b_ij]：

(1) $ {a_i} = \left\{ {\begin{array}{*{20}{c}} {1,}&{{\delta _i} \geqslant {P_i};}\\ {0,}&{{\delta _i} < {P_i}.} \end{array}} \right. $

(2) $ {b_{ij}} = \left\{ {\begin{array}{*{20}{c}} {1,}&{{R_{ij}} \geqslant {P_{ij}};} \\ {0,}&{{R_{ij}} < {P_{ij}}.} \end{array}} \right. $

式中：i，j=1，2，···，m. 若 ${a_i}$=0，则更新矩阵L，将所有与节点i连接的管段功能状态设置为0. 根据矩阵L得到 $m \times m$阶的邻阶矩阵A=[a_ij].

3）以水源点为起点，采用宽度优先搜索算法遍历邻接矩阵A. 将水源点放入队列Q中，搜索所有与水源点直接相连的节点，并将相连节点放入队列Q中，删除当前队列Q中第1个节点，即水源点出列；以队列Q中的第1个顶点为起点继续搜索，重复此步骤直至队列Q为空，在此过程中采用向量V=[v_i]（i=1，2，···，m）记录节点的访问状态，对于节点i，若它在队列Q中出现过，则v_i=1，表示节点i与水源点连通，否则设为0. 对于多源点的管网系统，可增设1个虚拟源点，并假设虚拟源点到各水源点连接管段的破坏概率为0，即可将多源点系统转换为单源点管网系统进行分析.

4）将步骤3）访问结果V累加到向量T=[t_i]中：

(3) $ {t_i} = {t_i} + {v_i};\quad i=1,2,\cdots,m. $

5）重复步骤2）~4）k次，即进行k次模拟，得到水源点到各节点间的连通概率矩阵P_c=[P_c_i]：

(4) $ {P_{\rm{c}}}_i = {{{t_i}} / k};\quad i=1,2,\cdots,m. $

须说明的是，步骤2）中产生随机数的算法可分为伪随机数算法和准随机数算法，两者最根本的区别^[17]在于：前者保证随机数序列的随机性，后者保证随机数序列的确定性和均匀性. 准随机数算法牺牲序列随机性，以此提高序列的均匀性，力求任意大小的样本都能满足低差异性. 目前，常见的低偏差序列有Sobol序列^[18]、Halton序列^[19]、Niederreiter序列^[20]等. CUDA框架中的cuRAND函数库提供伪随机数生成器和准随机数生成器，可用于在主机和设备端生成随机数. Sobol序列是唯一被主机和设备端API支持的准随机数序列^[21]. 因此，本研究选用Sobol序列作为拟Monte Carlo模拟的随机数. 为了提高准随机数序列的随机性，一些学者采用加扰技术对准随机数序列进行随机化，使其在保证低偏差的同时具有一定的随机性.

2. 可靠性分析的并行算法设计

2.1. CUDA编程原理

CUDA是英伟达公司为通用计算图形处理单元开发的并行计算平台和编程模型. 借助CUDA，开发人员可以利用GPU强大的并行计算引擎提高计算应用程序的运算速度. 在GPU加速的应用程序中，计算复杂、逻辑性强的部分在CPU上运行，而数据并行、计算密集、逻辑简单的部分在GPU上处理.

从软件角度看，CUDA编程模型通过网格（grid）、线程块（block）形式组织线程（thread），每个网格可包含多个线程块，每个线程块可包含多个线程. 线程以一维、二维或三维的形式组织，每个线程都有唯一的索引. 通过CUDA编程模型，开发者可以方便地设置网格、线程块的大小和维度. CUDA并行程序由串行和并行两部分程序组成，串行部分程序在CPU上执行，并行部分程序在GPU上执行，运行在GPU上的程序称为核函数（kernel），一个应用程序可有多个核函数，核函数以网格形式组织，一个核函数只能在一个网格上执行. 从硬件角度看，CUDA建立在流多处理器的概念上. 流多处理器是由多个标量处理器组成的单指令多线程单元，这些标量处理器能够执行特定的任务. 标量处理器只能访问自己的寄存器，因此寄存器不能用于在线程之间共享数据，而其他类型的内存对所有标量处理器都是通用的. 此外，每个流多处理器还包括一组本地32位寄存器、共享存储器、常量缓存和纹理缓存. 在流多处理器中，指令调用采用流水线设计，在指令级层面最大化并行处理能力. 在CUDA设备上除了流多处理器，还包括全局内存、常量内存、纹理内存和本地内存等存储器. CPU只能访问设备上的全局内存、常量内存和纹理内存. GPU具有对本地内存和全局内存的读/写权限，但对常量内存和纹理内存只有读权限. CUDA设备上的存储空间示意图如图1所示.

图 1

新窗口打开| 下载原图ZIP| 生成PPT

图 1 CUDA设备上的内存空间

Fig.1 Memory space on CUDA device

当启动一个内核网格时，网格中的线程块分布在流多处理器中. 一旦线程块被调度到一个流多处理器上，线程块中的线程会被进一步划分为线程束. 一个线程束由32个连续的线程组成，同一线程束中所有线程共享从存储器中获取的单条指令. 流多处理器以线程束为基本单位来调度线程，并将线程分派给标量处理器来执行.

2.2. 基于CUDA的供水管网可靠性分析算法

由可靠性分析串行方法可知，该算法须进行 $k$次拟Monte Carlo模拟，随着模拟次数及邻接矩阵规模的增大，计算时间也显著增长. 每一次管网可靠性分析都是独立的，互不干扰，因此，其具有天然的可并行性. 每一次管网可靠性计算可用一个线程完成. 对于拥有成百上千个计算核心（core）的GPU而言，即使模拟上万次，每个计算核心运行数次即可完成. 基于CUDA的供水管网抗震可靠性并行计算流程图如图2所示. 具体计算步骤可分为4步. 1）在CPU中输入每个节点及管段破坏概率、模拟次数等数据，并通过cudaMalloc函数在GPU上分配内存，通过cudaMemcpy函数把数据从CPU拷贝到GPU；2）启动核函数，通过CUDA框架中提供的cuRAND函数库及加扰Sobol序列生成器，通过scramble_c参数，设置加扰序列的初始值，在每个线程中生成随机数，并与管段及节点破坏概率比较，生成邻阶矩阵A，通过宽度优先搜索获得节点访问状态向量V，将结果保存在一维指针中；3）在并行计算完成后，将访问状态矩阵V拷贝到CPU进行求和及均值运算，确定最终的连通概率 ${P_{\rm{c}}}$；4）使用cudaFree、cublasDestroy和cusparseDestroy等释放设备内存、函数库句柄等资源，同时释放CPU上分配的内存.

图 2

新窗口打开| 下载原图ZIP| 生成PPT

图 2 基于CUDA的管网抗震可靠性分析并行计算流程

Fig.2 Parallel computing flow diagram of seismic reliability analysis of pipe network based on CUDA

对于初步完成的CUDA程序，须从内存、执行配置和指令等方面进行性能优化，以实现在保证精度前提下计算时间的最小化. 优化具体流程如下.

1）内存优化. 使用尽可能多的快速内存以实现带宽最大化. 在设备的6种内存中，寄存器的存储速度最快，共享存储器次之，全局内存的存储速度最慢. 由于在供水管网中节点及管段的破坏概率在所有的线程中都是相同的，本研究程序采用共享内存进行存储. 在数据传输方面，应最小化主机和设备之间的数据传输. 中间数据结构应该在设备内存中创建，由设备操作并在不被主机映射或复制到主存的情况下被销毁. 由于拟Monte Carlo方法会用到大量的随机数，采用主机端生成随机数并传递到显存上的方法较耗时，本研究程序在设备端生成随机数，节省数据传输时间.

2）执行配置的优化. 流多处理器应该保持较高的占用率，因为线程指令是在CUDA中顺序执行的，在一个线程束繁忙时执行其他线程束是隐藏延迟的唯一方法；每个线程块的线程数最好是32的倍数. 由于城市区域管网规模是变化的，线程块和线程块中的线程数应该根据所占用的寄存器空间及共享内存空间动态分配，以达到最佳效率. 针对本研究案例，在程序计算时，将线程块及线程块中线程数分别设置为8、64.

3）指令优化. 为了最大限度地提高指令吞吐量，应该尽量使用简单的算术指令；尽量避免流控制指令（if、switch、do、for、while）的使用，从而避免在同一线程束中使用不同的执行路径；在误差可接受情况下，使用单精度浮点值. 在本研究程序中，节点和管段的破坏概率及产生的随机数均采用单精度浮点值.

3. 案例分析

选取文献[10]中展示的某城市区域供水管网，采用本研究所提出的计算方法计算其系统抗震可靠性. 该系统共有49个管线段、33个节点. 其中，节点包括1个二级泵站和高地水池，其他节点为需水节点，视为普通节点，不考虑其破坏情况，管网平面布置如图3所示. 文献[10]给出了各管线在7、8、9度地震作用下的破坏概率，但其假设节点不破坏，各节点基本完好的概率为1. 由于本研究重心是研究管网可靠性分析并行计算方法的加速效率，故不再详细研究水厂和泵站等节点的破坏概率，在计算节点失效矩阵N时，假设节点可靠度为0.98.

图 3

新窗口打开| 下载原图ZIP| 生成PPT

图 3 某城市区域管网平面布置图

Fig.3 Layout plan of pipe network in an urban area

3.1. 计算精度比较

采用串行和并行计算方法，进行10 000次模拟，串行方法采用Matlab语言编程实现. 在模拟计算时，所采用的CPU为Inter（R）Xeon（R）E5-2407，4核，主频2.20 GHZ，4.00 GB内存；显卡采用GeForceGT 610，2.00 GB显存. 分析结果如图4所示. 图中，N为节点号. 对比7、8、9度地震作用下串行和并行计算结果，最大误差为0.52%，证明本研究并行计算方法的准确性.

图 4

新窗口打开| 下载原图ZIP| 生成PPT

图 4 某城市区域供水管网在不同地震作用下的连通概率

Fig.4 Connection probabilities of water supply pipe network in an urban area under different earthquake actions

3.2. 并行效率分析

以7度地震作用下管网节点和管段的破坏概率作为输入，改变拟Monte Carlo模拟次数k，采用串行和并行计算方法分别进行1 000、3 000、5 000、10 000次模拟，串行、并行模拟时间T_s、T_p和加速比R_s如表1所示. 可以看出，采用并行计算方法加速效果显著，最高为串行计算方法的96倍，计算时间提升较明显，说明GPU适用于管网抗震可靠性分析.

表 1 串行和并行计算时间对比

Tab.1 Comparison of serial and parallel computing time

k	T_s/s	T_p/s	R_s
1 000	9.260	0.116	80
3 000	28.351	0.313	91
5 000	47.333	0.498	95
10 000	93.723	0.972	96

新窗口打开| 下载CSV

3.3. 伪随机数序列与Sobol序列精度分析

分别采用伪随机数生成器和准随机数生成器在设备端产生随机数序列. 如表2所示，比较生成490 000个32位伪随机数序列和Sobol点列的计算时间T. 可以看出，两者所花费的时间相差较小，表明采用不同随机数生成器，不会对可靠性计算整体耗时产生太大影响. 采用并行计算方法，分别采用基于伪随机数序列的Monte Carlo（MC）方法和基于Sobol序列的拟MC方法进行管网抗震可靠性分析，模拟次数分别为100、1 000、5 000次，分析结果如图5所示. 通过与解析值对比，发现随着模拟次数的增加，模拟结果逐渐向解析值逼近，须强调的是，本研究解析值采用模糊数学法计算得出，详细计算理论见文献[10]. 基于伪随机数序列进行5 000次模拟，模拟值与解析值的最大误差为0.4%；基于Sobol序列进行1 000次模拟，模拟结果较接近解析值，与解析值的最大误差仅为0.2%. 表明拟Monte Carlo法能有效提高收敛速度并获得更高精度结果，从而缩短管网抗震可靠性分析的求解时间.

表 2 抽样点生成计算时间

Tab.2 Computational time to generate sampling points

序列	T/s
伪随机数	0.258
Sobol序列	0.256

新窗口打开| 下载CSV

图 5

新窗口打开| 下载原图ZIP| 生成PPT

图 5 Monte Carlo与拟Monte Carlo在7度地震作用下的模拟结果对比

Fig.5 Comparison of Monte Carlo and quasi-Monte Carlo simulation results under action of seven degree earthquake

4. 结　论

（1）对某城市区域管网系统进行抗震可靠性分析，对比7、8、9度地震作用下串行和并行计算结果，最大误差为0.52%，证明本研究提出的基于拟Monte Carlo法和CUDA的供水管网抗震可靠性分析并行算法的准确性.

（2）以7度地震作用下管网节点和管段的破坏概率作为输入，进行1 000、3 000、5 000、10 000次并行和串行模拟，并行模拟方法的最高的加速比为串行方法的96倍，有效缩短了管网抗震可靠性分析时间.

（3）利用伪随机数序列和低偏差序列作为随机数输入进行管网抗震可靠性分析，模拟结果显示，基于伪随机数序列进行5 000次模拟，模拟值与解析值的最大误差为0.4%；基于Sobol序列进行1 000次模拟，模拟值与解析值的最大误差仅为0.2%，表明基于低偏差序列的抗震可靠性分析具有更好的收敛效率及精度.

本研究所提出的管网抗震可靠性并行分析方法具有通用性，还可应用于其他生命线系统（供气系统、道桥系统、电力系统等）的抗震可靠性分析.

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

杨丹, 郭恩栋, 王祥建, 等

汶川地震供水系统功能失效模式

[J]. 灾害学, 2010, 25 (Suppl.1): 350