1、 )及应用图论中的最短路经法,及动态规鲻、线性规财等等。管网的非线性规划模型比较真实完整地反映了管网优化闯题的实质,但其隶解复杂。特别是大型管网问题中变量多计算困难,而且不能得到理论上的全局最优解。我国也有一些城市对管网作了现状核算和水力分析计算,为管网的扩建和运行管理提供了依据,并取的了可贵的经验,收到了良好的经济效果。本文选题结合导师科研项目,密切联系生产实际情况,对大塑给水管网的计算问题作些探讨。并应用近代数学发展的成果,结合前入的研究成果,在方法上作些补充和改进,意在促进这一领域的研究成果向实际应用方面转化,使计算机技术能更好地在大型给水管网方面投入实际应用,为国民经济建设服务,同时也
2、是为大型给水管网的优化计算打些基础。管网水力分析计算有许多算法,我校的杨钦教授、陈霖庆教授编制了一系列的水力计算程序,如解节点压力法或解校正节点压力法,解连枝流量法或解回路校正流量法,及同时解连枝和树枝流量法等。其中解节点压力法输入数据较少,系数矩阵为一正定对称的稀疏矩阵,在计算机上的求解要优于其它方法。本文选择这种方法结合实际的管网水力计算,就其在大型管网上的实际应用中的一些问题如存贮策略、机时节省数据的输入和结果输出等方面作些探讨。二、大塑给水管网的水力分析、给水管网的水力计算模式给水管网的水力分析计算就是要求得满足水力学连续性方程、能量方程及阻耗方程的一组水力平衡解。连续性方程是流向任一
3、节点的流量必须等于漉离该节点的漉量,以保持水流的连续性。即 叭 可 易 卜式中 一一 节 , 舭的 节点 流量 【一一联接节点和的管段流量(设漉出节点为正,漉入为负)一一与节点相联接的节点数。能量方程就是管网中消耗于每个环内的管段水头损失之代数和等于零(设顺时针方向为正,反之为负)。即 式中。 :, , ()一一环内第条管段的水头损失;一一环内的管段数;一一管网的环数(基本环)。管段的阻耗公式可表示为 式中 一一管段的摩阻系数;余同上。因此,能量方程(式()又可表示为 ()终合运用以上方程,可得到一些不同的算法。如解环方程法、解节点压力法和解管段流量法等等。由于阻耗公式是非线性的,决定了管网水
4、力计算要对其进行线性化处理后进而迭代求解。这是一个反复求解线性方程组的过程,本文采用直接线性化方法直接求解节点压力。这种方法的系数矩阵为正定对称阵性态好,而且计算中不用该变步长,人工干预少使用方便。 给水管网是一个管道网络系统,可用图论理论和方法对其进行分析(详见参考)。根据管网的线性图论模墅,管网的连续性方程可写为。 : )或写成分块形式:广 广 式中,一一管网的衔接矩阵( );一一未知压力点;一一已知压力点,一一节点流量向量(),一一管段流量向量()一一管网的节点数;一一管网的管段数。能量方程的矩阵形式为:或: 广, , 广 , 式中。一一节点压向量;一一管段水头损失向量;一一加压泵站的加
5、压向量。压降方程为: 设:式中。 趴 ( 即视与为线性关系,在迭代过程中逐次予以修正。由式()、()和式()整理可得: 十 将上式代入连续性方程(式 )得: 或 () 一 一 设: 一 一 蛳则得: 解上线性方程组即可解出给水管网的节点压力,进而求出管段流量相应的一组解,由于初始的为的函数,故式蛐需要经过多次迭代反复求解。即先视 为定值(也就是先假定初始流量),形威系数矩阵和列向量,然后求解 ,解得节点压力,由此得到一组新的管段流量,再计算新的 ,重新求解线性方程直到满足精度要求为止。计算步骤如下:、假设 , 、形成矩阵 ;、形成列向量一 一;、求解 ;、计算管段水头损失;()、计算管段流量(
6、);了、判断若: 一 列计算终止否则、 ( );(、计算-8 ;、转继续计算。程序框图见四给水管网的存贮策略)由上节知给水管网的水力计算就是反复求解的元素为: 的过程,系数矩阵 ( 火火大( ) 节点与节点相关联 节点与节点不相关联 一 由此可知的每一行(列)的非零元素个数为该行(列)对应节点的度(与该点相关联的节点数),而一般给水管网平均每个节点所联的节点数不超过个,也就是说系数矩阵中每行(列)的非零元平均不超过个。同时,矩阵中的零元个数大致同节点数成比例,管网规模越大,系数矩阵的零元个数越多,也就是稀疏程度越大。个节点的管网零元约占娜,个节点的管网则为郏,个节点的管网为 斯。一般称零元占优
7、势(蛔以上)的矩阵为稀疏矩阵,而相应的线性方程组称为稀疏的线性方程组,管网水力计算的线性方程组恰是这样的一种稀疏线性方程组。用经典的高斯消元法求解阶线性方程逐需要的计算机存贮量与的平方威正比,长运算数(乘除运算)大致与成正比,如在每秒万次的计算机上解 阶线性方程逐仅需要个存贮单元和少于秒机器运行时问解阶问题需要。,个存贮单元和秒的运行时问,若同样的机器解,阶问题 则需要,、内存单元和大约,)秒的机器运行时间,折合约,时,贮存量和机时的增加都是惊人的,随营管网规模的增加,这类问题也很突出 如本文的算洌节点数在以上。对于大型的稀疏线性方程组的求解不能直接应用通常适于满阵的高斯消去法或分解法。其原因
8、是:由于矩阵阶致较大,若存贮整个矩阵将占存贮量太大,甚至于存不下,因此需要特殊的存贮方案,一般的是仅存非零元或尽可能少的存零元 在高斯消去法或分解法中仅对非零元进行操作;更为重要是要设法在消去过程中尽可能的保持矩阵的稀疏性。这类问题表明,要恨据稀疏矩阵考虑相应的存贮办法和改进措施以节省存贮量和计算机时。在考虑稀疏矩阵存贮结构的同时必须考虑消去法,存贮方案不同相应的算法选择亦不同。从某种意义上讲,设计好稀疏矩阵的存贮结构是稀疏线性方程组求解算法的关键之一。稀疏矩阵的常用存贮方法有:线性表链表按带宽存贮线性表是将非零元按列(或行)存在一个一维数组中,为了给出非零元的列号(或行号),再引入数组,它记
9、录数组中相应元素的行号,另有一数纽记录每列第一个非零元的位置。这种存贮方法所需的存贮量为:(为非零元个数,为矩阵阶数),它的优点是存贮量小,结构简单;缺点是不便于插入和删除,如插入一个非零元时位于它下面的所有非零元都必须向下移动一个位置,浪费机时。单链表克服了上述线性表的插入不便,即在线性表的基础上增加了一个链指针效组,()表示第非零元的下一个非零元的位置,同时设置效组记录各列(或行)非零元的个数,其优点是便于插入和删除。但增加了存贮量,而且此法有个缺点是访问链表中任一元表,均需要从一列到另一个非零元开始,期序地找下去,而且均需要问接访问(通过的访问),所以重访时间较线性表慢。带宽存贮分为一维
10、压缩存贮和两维存贮,其中一维变带宽压缩存贮是矩阵中每列(或行)的第一个非零元开始存连同带宽内的零元一起存入一个一组数组中,另外采用一个指标数组来记录对角元所在中的位置,这种方法适于结构对称的线性方程组,而且采用了变带宽较二维存贮节省单元,这一方法也在管网的水力计算中得以应用,但尚有一点需要进一步讨论的是在半带宽内尚有大量的零元要占内存并且参与计算,而且大型管网问题节点标号的编射尚无方法可有效地减少这部分零元,所以这一部分的浪费在大型管网的计算中更是不容忽视的。除以上三种方法外还有双链表,和分块存贮等方法,这些方法更适于结构不对称的稀疏矩阵的存贮问题,考虑到管网的计算特点这里不作介绍。对应于以上
11、存贮方法可派生出许多线性方程组的解法,大型稀疏线性方程组有许多解法,其取舍条件是保持利用系数矩阵的稀疏性,一般的标准是在分解过程中的局部填入量找长运算数最小,填入元即为在矩阵分解过程中原来零元的位置可产生非零元。一个稀疏矩阵在高斯消去或分解过程中,原来的非零元可能变成非零元。如在分解中, 一 :, 易 钆 , 一 一 卜 虹 一 ,式中 为分解的下三角阵元素 为分解的对角阵元素计算中 等于零时若某一 和 不等于零则使得 不等于零产生填入元,入图所示:广 分艉前为零元,不等于,而 一 而使得 不等于零即为填入元。长运算数即指乘除运算次数,局部长运算算数最少,印指消去过程中每一步。在所剩余的矩阵中
12、选择局部操作数最少的对元为主元。上述两个原则最有代表住的两个算法是廷尼一袄克算法和马尔科威算法。但就给水管网计算的性质来说是求解正定对称的稀疏线性方程组,以上两个准则则等价于非零元个数最少,如采用半带宽一维压缩存贮则为半带宽内的非零元最少,本文拟采用压缩一维存贮分解和结合重序过程的作法,以求在所占内存和计算机时间上有所改善。本文所采用重序算法()援为基于图论的理论。在分解过程中每次选取使半带宽最小的对角元为主元。这是一个模拟分解定序的过程,也相当于对管网的节点标号重新排序,不过这个过程是在对系数矩阵进行数值分解的外部进行的。在介绍重序算法之前有必要介绍几个图论的知识。距离:在连通图( )中两个
13、节点,弘连接如的最短路径的长就称作两点问的距离,记作(,)。偏心率:一个节点的侗心率是该点与图内所有点的距离的最大值,记作()即:()皿(如) 直径:连通图( )的直径为图内所有节点偏心率的最大值,记作()即:():() 图的分层结构:结点集合的一个划分,记作 即: , , )() ()( ) : () 使, ()为节点集合的邻集,即与内节点和相关联的所有节点的集合。由上定义可知每个,使连通图(,)产生一个包含两个连通图的截图()和( 、),即为连通圈的隔离集,简称隔集。数:( )为分层结构的长度。分层结构的宽度定义为:( ) 对于节点,根在的分层结构是: ()(),(), ,()其 中 ,
14、( 卜虹 ( )本文所用方法是按照某一节点的分层结构对每层内层所有节点 ()按其度大小顺序进行顺序编号。对于一般的线性图来说,很多顺序算法算法的效率与 ”起始节点 ”的选择关系极大,大量研究证明最好选用外围节点即():()的节点为 ”初始节点 ”,但目前尚无有效算法能找出一殷图的外围节点,等人给出了一个有效的算法能从线性图中找到偏心率接近图的直径的节点,即伪外围节点,算法如下:找一个有最小度数的节点;产生根在的分层结构, ()(),(), ,()();将()(即最后一层)中的节点按度数的增加顺序分类;按度数增加顺序依次从()中选出 ()()产生根在分层结构,若()()则将赋给转步骤(),若所有
15、()()则为伪外围节点。这个算法的缺点是集()()可能很大,在第()步骤需要执行 ()()次,计算量很大,对大墅管网问题,与实际求解运行时间相比显得得不偿失(这里是个用机时换空问的问题),可喜的是城市一般不是特别狭长或形状奇特的图形,而是方园在一定范围内较均匀的图形,特考虑省去寻找伪外围节点的过程,而是在管网图形的长的方向外围选一度数较小的节点代替,这样虽可能使结果内存贮有所增加但省去了许多机时。根据本文算例结果来看,内存的增加很少,结果是令人满意的。在本文的算例中节点数为 ,经重序过程优化后一维存贮敷组下界为 , ,与有关资料介绍的算法所得结果相近。如有资料介绍: 时,占 ,单元 时占, 单
16、元相比之下结果也是满意的,就是说至少对本算例来说因伪外围结点选择而造成的内存增加量不明显,这个过程在西门子机上所用机时接近秒,占总计算时问的左右,从这点上看再耗费机时去减少为数不多的内存可省量也是不合适的。正如前面所介绍的,结构对称稀疏矩阵线性方程组的非零元最少也等价于长运算敏最少。该过程对管网水力计算的速度也有可观的改进,如本文算倒中手工编于未经重序过程计算时间为 秒左右,重序后计算时问则在 秒以内,重序过程所用计时在 秒以内。节省机时 知左右。由此可见本文所编过程是可行的,而且效果良好。重序过程算法如下:()在管网选一长向外围度数小的节点,存入;()寻找与相关联之节点,送效组;()对内诸元
17、素(节点)按其度数大小顺序重新标记。存入若为空则结束,否则,转步骤()继续运行。程序框图见四计算机重序过程即可以节省存贮单元又可以简化管网节点的编号过程,对于较小的管网手工编号可以取得较好的结果。但对大型管网则不容易作到,而且编号过程也较费时、费力,而使用本程序则无需考虑其它因素,只需标注节点号即可,这在多方案计算,已知压力改变和插入节点时尤为方便,水力计算方法决定了已知压力是在编号的末尾,手工编号在已知压力点是变动时要对节点编号作相应改动,而重序过程的使用则省去了这一部分的工作,对插入节点也只要顺序编号即可。、 大型管同实例水力计算中遇到的问题水力计算的收敛性问题已有专门的论证,下面只就大型
18、管网计算中收敛方面的一些情况及采取的措旌作些讨论。由于管网规模的扩大,计算量相应增加,尤其在求解 方程组的分解过程中,计算量增加的更显著。这样舍入误差带来的影响也突出出来。在水力分析计算中,如前所述是反复求线性方程组 的过程,精度控制为相邻两次计算管段流量差的最大值,所以对求解线性方程组的精度要求不能只满足于其解值满足工程上的需要即可,而且要考虑值的误差对管段流量的影响。工程上一般达到 米足可满足要求,但这个精度可能造成管段流量的较大误差,不满足一般的升秒的要求。同时流量精度的达到也是对前面线性化假设的一种检验,娶从算法完整性的角度来讨论求解线性方程组精度要求问题。对求解线性方程组的台入误差分
19、析,可采用向后误差分析法,其基本思想就是计算过程中舍入误差的影响归结为原始数据变化的影响,即找出原始数据的某种变化使其对最后结果的影响同求解过程中的舍入误差等效。这样计算所得近 似解将严格满足下面的方程式:(十 )式中: 一一系数矩阵,一一列向量,一一解向量, 找到这样的 和 问题就转化为摄动分析问题了。这里只引用结论。对于平方根法解的精度估计如下月 一()( 一)(式中()一一矩阵的条件效;一一矩阵阶乘;-25一一计算机字长-8一一常系致;-5,) 一 一一向量的无穷范数; 上式表明解的精度除与矩阵性态(条件数()有关外,与矩阵的阶数及所用机器字长有关,那么对于不同阶数的问题,若对解要求保持
20、同一精度,可考虑提高步长。由阻耗公式 考虑舍入误差,十 (十) 一 一一一一一一一一一一 一一 一 ) (由此可得 十 ( )( ) ( ) ()() () ()。 ()( ( )十 () 一一 一一( ) , () 若 升,秒,管段最大流量为, 升,秒的话相应的对 求解中的精度要保持 若压力在 米以内, 的精度则要求小于 米而不是工程上的米,这是算法整体性的要求。对一般规格的管同的计算中,大多可满足 对 的精度要求,不出现含入误差对结果的影响,大型问题的阶段增大则可能使得求出的值的误差结果造成 () () , ,的精度降不下来而影响收毁。计算中的表现为前面几次迭代 的下降正常计算所得值也属正
21、常范围,前后两次的值的差值不大于 ,可 总在一定的范围摆动收敛不下来。这时即可采用双精度的办法,来消除舍入误差的影响,以保证收敛性。此法在实际计苒中得到预明效果, 管网计算的一些辅助手段一一图形处理计算机科学的飞速发展为人们提供了越来越有力的工具,它不仅能完成大量的繁杂的科学计算,而且还具有逻辑判断和图形处理方面的功能。计算机的图形功能的开发应用可使得计算的结果图形化满足工程上的应用和习惯。过去常用的管网计算程序的汁算结果是用打印报告的数据表格形式,读来不直观 不符台工程上的习惯,而且将计算结果填到图上的工作不仅烦琐而且易出现错误,管网舰漠越大这一问题越突出。本文拟就给水管网计算图形的数据输入
22、(拓扑关系,几何尺寸及有关参数),数据文件的生成及计算结果的输出过程的计算机自动完成作了些探索与尝试,力求充分地开发计算机的功能,简化一些繁索的敏据整理与输入过程,和计算结果的图形输出,以便于对结果进行分析。满足工程上的应用需要。给水管网计算需要输入的数据有节点流量、管径、管长、管段的粗键系数及拓扑关系等等。常用的方法是将管同的计算图形中的管段和节点进行标号,然后再整理成数据表的形式,最后从终端键盘输入。解节点方程一般需输入以下几个数组: () ()()一一管段上游节点号;一一管段下游节点号;一一管段管径; ()一一管长;()()一管段的租糙系数;一一节点流量;这些数组的输入工作的一部分,可利
23、用计算机的图处理功能进行简化转为计算机完成。如管网的拓扑关系,即)、()数组和直线管段的管长等,而且现有计算机的功能和软硬件配置也具有开发这种过程的条件。本文所编过程是在 一 软件包支持下,在 机上实现的,程序用 语言编制。 一 是一种目前在国内外微墅计算机系统中最普遍应用的一种高性能作图软件包。它具有适应性好、实用性强、适应性广等特点,它的图形文件经过适当的格式转换可以被 等其他程序或高级语言调用,本次使用的是 一 版本,配有 一 语君。 一 是一种嵌入在 一 中的 语言,编程功能强,利于用户开发图形功能。 语言是一类表处理语言,具有变量运算、递归、循环等功能,使用比较方变。本文所编过程的硬
24、件要求除主锣外,还需配置一台数字化仪用于图形的输入,和一台绘图仪用于图形的输出。该方法的主簧过程是将管网图形在数字化仪上 ”描 ”入计算机,同时自动生成 , 数组及直线段管长,配合 ”描 ”入过程同时在键盘上按描入次序输入管径,管段粗糙系数及节点流量。这样处理简化了数据整理过程和键盘输入量,而且将图形送入了计算机,以后的计算结果的输出只需将相应的数据填入图中即可。上述过程完成后,整个管网计算的已知数据(包括数据文件和图形文件)等已进入计算机,可对其进行水利计算。计算结果可生成供 程序调用的数据文件,也可以生成一串绘制结果的 语句,然后在 一 下启动输出过程将计算结果的数据填入管网图上。图形的输
25、出过程:)将管网图定位于数字化仪上。逐条描入臂程,同时按提示从键盘上输入值及相关点的节点滤量 过程框图见四中 过程框图。着使用中数字化仪太小,描入图形有困难,可输入节点坐标值,当然这样输入量有所增加,输出的结果无影响。本过程的例图谤见附录。本文的绘图程序是利用 的西线和图块功能,用 语言编铹。图块是由一组图形实体组合起来的复杂物体,一旦这种组合完成以后,这些图块便给定了一个图块名,可以在绘图过程中按需要在任意点处插入。图块的属性是随图块每次插入而变化的文字信息,它可以显示在图形中。具体的作法是,先将节点和管段的有关数据的图上格式射成块,并将欲显示的数据定义为属性,以便图中以文字的形式显示。这样
26、可以将一个块配以一组属性值,依次在插入点处绘钼出该点计算结果或已知效据。也就是说依次对节点按坐标值插入节点数据块,并输出该点的结果数据作为属性的响应,然后再依次对管段按中点坐标和倾斜角度插入譬段数据块,以管段的结果数据作为属性响应,这样就完成了计算结果的图形绘制工作。该图形可在计算机显示器上显示,亦可由绘图机绘出。在数据输入阶段,以上过程是同 ”描入过程相结合同时进行的,数据输入结束,上面标有已知数据的管网图形亦即生成。限于篇幅,这里不再作详细的介绍。程序框图见四中 过程框图。给水管网图形文件的生成使得给水管网的数据整理、校核工作带来了许多方便,尤其是大型臂网计算,数据近一千组(如本文实倒),
27、效据整理、输入工作很繁琐,在管段,节点编号及管段和节点数据的输入中难免出现错误。使用本文所编的输入过程可减少一些中问环节,减少出错机会,而且校核很方便。如节点标号及管段上、下游节点编号数组,一旦出现错误则造成管网的拓扑关系变化,管网的联接关系混乱使计算得出错误结果或引起不收敛,而且较难发现和查寻,若利用图形功能与原图对照则一目了然。而且在图上校核数据也较方便,符合工程设计的习惯。对计算结果的分析,借助于图形也较过去的常用方法来的方便和明了。在本文的工程实例中已体会到了这一点。三 给水管网的现状核算为了合理的进行管网的改建和扩建工程及确定合理的调度运行方案,应首先对现有管网进行复核计算。通过复核
28、,我们可以摸清管网的现状作为改建、扩建及调度的依据。给水管网是一个大型的网络系统,工程设计计算考虑的只是干管(即管径大于某一尺寸的管段),用水力学的一些方法,将千管沿线支管分配出的流量折算成节点流量丽形成管网的计算简图。这个过程有一定的近似,存在一些误差,尽管对管网的计算可求得精确的数学解,但不可能与实际运行状况完全一致,更何况管网的一些参数,特别是粗糙系数、水泵特性参数、节点流量等的不准确是难以避免的,都可能造成计算结果与实测值不吻合的情况但只要处理得当,数据相对准确,这些误差造成的影响可以控制在一定的范围内,满足工程上的需要。管网的现状核算就是要通过一些算法来调整管网的参数,使得管网计算模
29、型较真实地反映管网的实际状况。给水管网现状核算的一般作法:给水管网现状核算大致可分为三个部分:)现场数据调查和测定)进行水力平钠计算)比较铡压点的实测值和计算结果的差异若不满足精度要求,列)依一定的算法调整节点和管段摩阻,使压力的实测值和计算值的相差保持在精度要求范围内。以上算法的主要依据是:一般的讲,计算结果与实际测定结果的偏差主要是由于节点流量不准确和管段摩阻(它包括管径和粗糙系数不准确两者的影响)的不准确造成的。不同的调整方法形成了不同的算法,主要有灵敏服分析法和最速下降法。灵敏度分析法是求出节点的压力对节点流量和管段摩阻的灵敏度,然后选择这个系数较大的节点或管段进行调整,以使得测压点处
30、压差迅速而有效地消除。最速下将法认为:为消除测压点处实测压力与计算压力之闻的差值主要是要调整节点流量,而忽略管段摩阻的影响。另外还有的算法在灵敏度分析的基础上,结合优化理论,优化调整过程使得调整量最小的算法。以上算法的共同基础是基于现场同步铡流测压结果来进行调整的。可是由于种种原因所限不能对管网进行实测时,复核工作显得难以进行。自来水公司长年运行经验可对管网的现状有个总体上的认识,历年的测压资料(不是同步测值),抄表簿上的水量资料, 工况企业的分布和用水量情况以及调度数据等均可作为复核的参考资料,这些数据虽不能满足现状核算模式的数据要求,但可以依此数据作出初步的假定,然后计算根据结果的情况对原
31、始假定数据进行调整,可以说是一个试算的过程来完成。我们即采用了这种方法同协作单位的同志作对本文的算例作了管网核算的尝试,这是个繁琐的过程,好在协作单位的同志对管网的现状情况非常熟悉,使得这个过程得以顺利完成。这种作法,从算法上讲比较原始,而且所得结果也较难用精确的数量关系来评价,但对这种无力进行现场同步测流、测压的特定条件下,还是有一定意义的。就本文算例来说,结果较为满意。其计算结果经有关方面认定,认为基本可以反映了实际管网状况,可以作为下一步规化的依据。四变量说明与程序框图(一)、变量说明瑚一一管网节点数;一管网管段数;一一未知压力点敦;一一加压泵站个数; 一一 水源个数; 一一已知压力点个
32、数; () 一一 管段上游节点编号数组; ()() ()()() ()一一管段下游节点编号数组;一一管段管径数组;一一管段管长数组;一一管段粗糙系数投组;一一管段摩阻数组;一一管段水头损失数组; () 一一管段流量致组; 【) 一一管段流量暂存敦组;)()()一一节点漉量效组;一一列向量;一一带点压力数组; ()一一 加压泵站所在管段编号数组; () ()一一加压泵站的静扬程;一一加压泵站的摩阻系数; () () () ()一一一一一一一一加压泵站的内部阻力系数水源所节点编号;水源供水量;节点新编号; ()一一 ()一一已知压力节点标号;已知压力点压力值;水力分析程序幅圈 了口过程报圈 弋:!罗 弋:罗 ,雨念过程怄豳查寻口 的缡号()查寻口的翁号(寓 捕捉管殴的上潍节点口 坐标存入(),()键入()扛(捕捉管段的上游啭点口 口坐标有:入(),()礁入(扛()珂