基于粒子群優(yōu)化算法的管網(wǎng)摩阻系數(shù)的反分析p

1、 基于粒子群優(yōu)化算法的管道摩阻系數(shù)的反分析劉力鵬1，郄志紅1*（1 河北農(nóng)業(yè)大學(xué) 水利工程系，河北 保定 071001）摘要：管道摩阻系數(shù)是管網(wǎng)設(shè)計、運行調(diào)度計算的重要水力參數(shù)，管道摩阻系數(shù)的反分析計算，對管網(wǎng)的優(yōu)化調(diào)度、改擴(kuò)建設(shè)計等具有重要意義。將粒子群優(yōu)化算法（PSO）與節(jié)點水壓法水力計算相耦合，提出了基于粒子群優(yōu)化算法的管道摩阻系數(shù)反演新方法。通過對某管網(wǎng)進(jìn)行的實例分析，結(jié)果表明該方法計算速度快、易于收斂到全局最優(yōu)解。關(guān)鍵詞：管網(wǎng) 粒子群優(yōu)化算法 摩阻系數(shù) 節(jié)點水壓法Inverse analysis of frictional coefficients of pipelines base

2、d on PSOLIU Li-peng1 ，QIE Zhi-hong1* ，ZHANG Guo-zhan2(1 Department of Water Conservancy Engineering, Agriculture University of Hebei, Baoding 071001, Hebei, China2 Baoding drain the head office, Baoding 071001, Hebei, China)Abstract：Frictional coefficient is an important hydraulic parameter for pipe

3、 network design and scheduling calculation. The inverse analysis to frictional coefficient has much meaning to pipe network optimization and reconstruction. Combining the particle swarm optimization with node water head method, a new inverse analysis method to frictional coefficient is put forward.

4、The example analysis result shows that the method is rapid and easy to find the globally optimal solution.Keywords： Pipe network; particle swarm optimization (PSO); frictional coefficient; node water head method中圖分類號：TU991.61 文獻(xiàn)標(biāo)識碼：A1 引 言作者簡介：劉力鵬（1983 - ），男，河北石家莊人，碩士研究生；研究方向：水利信息管理。保定 河北農(nóng)業(yè)大學(xué)城鄉(xiāng)建設(shè)學(xué)院，0

5、71001。Email:liulipeng_* 通訊作者：郄志紅（1969 - ），男，河北徐水人，教授，河北農(nóng)業(yè)大學(xué)博士生導(dǎo)師；研究方向：主要從事節(jié)水灌溉，水資源及水信息學(xué)方面的教學(xué)和研究工作。 對于給水管網(wǎng)中的管道，隨著使用年限的增加，由于水質(zhì)的影響，導(dǎo)致管壁腐蝕和結(jié)垢，摩阻系數(shù)C隨之增加。有資料報導(dǎo)1：某條DN1200的鑄鐵管，僅使用兩年，就因管內(nèi)積垢，粗糙系數(shù)n值從0.013增加到0.023，管徑縮小到900mm。這一增一減，使得摩阻系數(shù)C增大到原來的13.3倍。即使是同一管徑、同一管材、同時埋設(shè)的管道，也因水質(zhì)情況與水管的工作條件不同，其C的變化也不同。因此，在管網(wǎng)的運行管理和改擴(kuò)建

6、設(shè)中，不能套用管道原設(shè)計的摩阻系數(shù)，而須通過對管網(wǎng)的監(jiān)測進(jìn)行反分析，從而推算出管道實際的摩阻系數(shù)。在管網(wǎng)反分析的實際問題中，觀測量與反演量往往呈復(fù)雜的非線性關(guān)系。目前，對于非線性模型參數(shù)估計問題2，有許多成熟的系統(tǒng)辨識方法，如最小二乘法3、極大似然估計法3、遺傳算法4等等。但是最小二乘法和極大似然估計法都是在梯度方向上尋優(yōu)的局部搜索方法，在某種情況下可能陷入局部極值。利用遺傳算法對參數(shù)進(jìn)行估計時，雖可避免陷入局部最優(yōu)解，但要涉及繁瑣的編碼、解碼過程，影響其效率。由Eberhart博士和Kennedy博士提出的粒子群優(yōu)化（Particle Swarm Optimization, PSO）算法5是

7、一種有效的隨機(jī)全局優(yōu)化技術(shù),已經(jīng)被證明是一種可靠的優(yōu)化方法6。PSO算法同遺傳算法類似，是一種基于迭代的優(yōu)化工具，系統(tǒng)初始化為一組隨機(jī)解，通過某種方式迭代尋找最優(yōu)解。但PSO算法并沒有遺傳算法中的交叉以及變異算子，編碼方式也比遺傳算法簡單7。本文引入粒子群優(yōu)化算法來對給水管網(wǎng)進(jìn)行反分析。2 粒子群算法2.1 粒子群算法介紹PSO算法是模擬鳥群捕食行為的一種方法。設(shè)想這樣的一個場景：一群鳥在隨機(jī)搜索食物，在這個區(qū)域里只有一塊食物，所有的鳥都不知道食物在那里，但是他們知道當(dāng)前的位置離食物還有多遠(yuǎn)。那么找到食物的最優(yōu)策略是什么呢？最簡單有效的就是搜尋目前離食物最近的鳥的周圍區(qū)域。在PSO算法中，優(yōu)化

8、問題的解是搜索空間中的一只鳥，可以稱之為“粒子”。把優(yōu)化問題的解初始化為一群隨機(jī)粒子（隨機(jī)解），然后通過迭代找到最優(yōu)解。在每一次迭代中，粒子通過跟蹤兩個“極值”來更新自己。第一個極值是粒子本身截止到目前所找到的最優(yōu)解，這個解叫做個體極值pbest ；另一個極值是整個種群截止到目前所找到的最優(yōu)解，這個解叫做全局極值gbest。在找到這兩個最優(yōu)值后，粒子根據(jù)、兩個公式來更新自己的速度和位置5： 式中V 是粒子的速度；present 是當(dāng)前粒子的位置；pbest 和gbest如前定義；rand()是介于（0，1）之間的隨機(jī)數(shù)；c1、c2是學(xué)習(xí)因子，通常c1= c2= 2。2.2 PSO算法求解最優(yōu)化

9、問題的算法流程 初始化粒子群，包括群體規(guī)模N，每個粒子的位置Xi和速度Vi8 。 計算每個粒子的適應(yīng)度值Fit i 。 對每個粒子，用它的適應(yīng)度值Fit i 和個體極值pbest i 相比較，如果Fit i pbest i ，則用Fit i 替換掉 pbest i 。 對每個粒子，用它的適應(yīng)度值Fit i 和全局極值gbest相比較，如果Fit i gbest，則用Fit i 替換gbest，按照一定更新概率或根據(jù)式(1)、(2)分別更新粒子的速度Vi和位置Xi，或隨機(jī)初始化更新。 如果滿足結(jié)束條件(誤差足夠好或到達(dá)最大循環(huán)次數(shù))退出，否則回到。3 管道摩阻系數(shù)的反分析在摩阻系數(shù)的反演中，觀測

10、量（某些節(jié)點的水壓）與反演量（摩阻系數(shù)）間呈非線性關(guān)系，一般采用非線性優(yōu)化的方法搜索滿足目標(biāo)函數(shù)的最優(yōu)解，此時的反演便轉(zhuǎn)化為最優(yōu)化求解問題。在摩阻系數(shù)的反分析過程中，需要調(diào)用管網(wǎng)的正向計算。一些研究結(jié)果證明1，管網(wǎng)正向計算中的環(huán)流量法的流量優(yōu)化分配問題是一個凹規(guī)劃問題，而節(jié)點水壓法就沒有這一問題，而且此方法比環(huán)流量法收斂速度快、精度高、存貯量小，并且不用考慮初始流量的分配。本文采用節(jié)點水壓法作為管網(wǎng)正向計算的方法。3.1 節(jié)點水壓法基本公式1：海曾-威廉斯(Hazen-Williams)公式： （i = 1,2,M）式中，=/，是節(jié)點i、j之間的管段流量，是節(jié)點i的節(jié)點流量閉合差，、D分別是管

11、長（m）、管徑（m）,摩阻系數(shù)C與管材有關(guān)，=0.54，M是節(jié)點數(shù)（不包括水塔節(jié)點），是節(jié)點流量（m3/s）,、是各節(jié)點的水壓值，SGN為符號函數(shù)：基本步驟： 給定水壓初值，由公式(4)計算（i = 1,2,M）。若所有節(jié)點i均有（為所需精度，可取 = 0.001 m3/s），則為所求；否則，令（k表示迭代次數(shù)，即初始化迭代次數(shù)）。 由式計算 （i ，j=1,2,M）, 系為在 = 處的值，按式進(jìn)行計算形成Jacobi矩陣； 解下列式線性方程組，求出（i = 1,2,M） 式中稱為校正水頭或水頭校正值； 按式求出的新的近似值（i = 1,2,M） 由式（4）計算（i = 1,2,M），若所有節(jié)

12、點i均有，則為所求；否則令，返回步驟。3.2 管道摩阻系數(shù)的反演反演步驟： 確定PSO算法的運行參數(shù)，包括最大迭代次數(shù)、速度和位置的運動范圍、慣性權(quán)重等。 初始化微粒種群，包括群體規(guī)模，每個微粒的速度Vi和位置Xi（位置Xi與摩阻系數(shù)c相對應(yīng)）。 通過節(jié)點水壓法計算管網(wǎng)各觀測節(jié)點的水壓值，并計算觀測節(jié)點處水壓觀測值與計算水壓值的二乘誤差，將二乘誤差值作為個體微粒的適應(yīng)度函數(shù)值Fiti，評價各微粒的適應(yīng)度函數(shù)值Fiti。 對于每個微粒，將其適應(yīng)度函數(shù)值Fiti與歷史最好位置pbesti相比較，若當(dāng)前適應(yīng)度函數(shù)值Fiti較優(yōu)，則用適應(yīng)度函數(shù)值Fiti來更新pbesti。 將每個微粒的適應(yīng)度值Fit

13、i與群體經(jīng)歷過的歷史最佳位置gbest相比較,如果當(dāng)前群體中最好的適應(yīng)度值Fiti較好,則將其作為新的gbest. 根據(jù)（1）、(2)兩式分別更新每個微粒的速度Vi和位置Xi，檢驗gbest是否滿足最小誤差，若滿足，則微粒的位置Xi即為摩阻系數(shù)c的反演值；若不滿足，則回到。摩阻系數(shù)反演流程圖如圖1所示。確定PSO算法的運行參數(shù)隨機(jī)初始化群體，包括位置和速度把各觀測節(jié)點水壓值與相應(yīng)的計算節(jié)點水壓值之差的平方和設(shè)為個體適應(yīng)值Fiti對于每個微粒，其適應(yīng)度值Fiti與其經(jīng)歷的最好位置pbesti相比較，將較優(yōu)者作為pbesti每個微粒的適應(yīng)度值Fit i 與整個群體所經(jīng)歷的最好位置gbest進(jìn)行比較

14、，將較優(yōu)者作為gbest更新微粒的速度和位置滿足迭代終止條件嗎？NY得到摩阻系數(shù)c的反演值結(jié) 束圖1 摩阻系數(shù)的反演流程Fig. 1 Calculating flow in inverse analysis of frictional coefficient4 實例如圖2所示，為某示例給水管網(wǎng)，其中管道分為兩類：鋼管（摩阻系數(shù)設(shè)計采用值c=120）、混凝土管(摩阻系數(shù)設(shè)計采用值c=110)，各節(jié)點出流量 (m3/s)、管長(m)、管徑D(mm)，兩個泵站水源，一個水塔水源。泵站水源的水頭方程分別為（Q為21號當(dāng)量水管的管段流量）,（Q為22號當(dāng)量水管的管段流量）,水塔水面標(biāo)高為74m，圖中,分

15、別表示水源1，2，3，節(jié)點的數(shù)字為節(jié)點標(biāo)號。4.1 基于PSO的管道摩阻系數(shù)反分析計算 圖2 各管段尺寸及各節(jié)點出流量Fig. 2 Pipelines arrangement and discharged flux of nodes圖3 各管段編號及計算模型圖Fig. 3 Serial numbers of pipelines and calculation model管網(wǎng)共有22個管段（其中，18號管段是混凝土管；920號管段是鋼管；21、22號管段為當(dāng)量水管），14個節(jié)點，如圖3所示。在第1，3，10節(jié)點處分別設(shè)立三個水壓監(jiān)測點，測得兩年以后三個監(jiān)測點的水壓值（不考慮管網(wǎng)的鋪設(shè)年代對管道摩

16、阻系數(shù)的影響），見表1所示。把水壓監(jiān)測值作為觀測量對管網(wǎng)進(jìn)行反分析計算 ，得到由PSO反演出的摩阻系數(shù)反演值。摩阻系數(shù)反演值的準(zhǔn)確度由反演精度評價函數(shù)evl來控制，evl按式計算。 式中為監(jiān)測節(jié)點的計算水壓值，為監(jiān)測水壓值。在節(jié)點水壓法正向計算中，節(jié)點流量允許閉合差設(shè)定為0.001m3/s，步長因子設(shè)定為0.58，最大迭代次數(shù)設(shè)定為500；在PSO反分析計算中，將種群數(shù)設(shè)定為15，最大迭代次數(shù)設(shè)定為500，初始權(quán)重設(shè)為0.95，權(quán)重以線性方式變化，終止迭代誤差（精度）設(shè)為0.001。摩阻系數(shù)反演的收斂情況如圖3所示。計算結(jié)果見表1、2。表1 各測點觀測水壓值及摩阻系數(shù)反演結(jié)果Tab.1 Obs

17、ervational hydraulic pressure values of nodes and inverse analysis calculation results of frictional coefficients節(jié)點號節(jié)點觀測水壓值/m混凝土管摩阻系數(shù)反演值鋼管摩阻系數(shù)反演值終止迭代時的evl值1 77.06132.034348153.0924220.0009463 74.8710 73.39圖4 摩阻系數(shù)收斂曲線Fig. 4 Convergence curve of frictional coefficient表2 在管網(wǎng)正向計算中與反演結(jié)果相對應(yīng)的各節(jié)點流量閉合差Tab.2

18、Closed error of flux correlated with reverse calculation result of nodes節(jié)點號節(jié)點流量閉合差節(jié)點號節(jié)點流量閉合差節(jié)點號節(jié)點流量閉合差10.0008460.000543110.00009620.0003270.000196120.00004630.00035980.000144130.00010740.00080890.000199140.00051950.000991100.000262把表1中的摩阻系數(shù)反演值與管道摩阻系數(shù)設(shè)計采用值進(jìn)行比較可知，管網(wǎng)的摩阻系數(shù)變大，說明管網(wǎng)經(jīng)兩年運行，管道內(nèi)壁可能產(chǎn)生了沉積物和水垢，在一

19、定程度上影響了管道的配水能力。由反演收斂圖4可知，PSO算法收斂速度很快，收斂精度高。4.2 基于PSO的摩阻系數(shù)反演結(jié)果的穩(wěn)定性分析由反演結(jié)果分析可知反演結(jié)果是正確的，但是管網(wǎng)的實際運行會受到各種因素的影響，節(jié)點水壓的監(jiān)測難免會出現(xiàn)誤差，因此，管網(wǎng)反演模型的穩(wěn)定性問題是一個非常重要的問題?，F(xiàn)刨除觀測節(jié)點水壓時產(chǎn)生的誤差，對基于PSO算法的反演模型引入噪聲誤差，由噪聲的加入對反演模型的影響，來判斷基于PSO算法的反演模型的穩(wěn)定性。將表1中的三個節(jié)點觀測水壓值設(shè)為三個節(jié)點的真實水壓值（不含誤差），則在這三個水壓值之上隨機(jī)加入0.03m的噪聲（共有7種組合），計算結(jié)果見表3。 表3 反演模型加入噪

20、聲后的摩阻系數(shù)Tab.3 Frictional coefficients result with noise組合數(shù)混凝土管鋼管加噪聲后計算 出的摩阻系數(shù)摩阻系數(shù)誤差加噪聲后計算 出的摩阻系數(shù)摩阻系數(shù)誤差1132.034211-0.000137153.0926050.0001832132.0755620.041214153.7181930.6257713132.0758060.041458153.7182610.6258394132.3482670.313919152.400039-0.6923835132.3484650.314117152.400032-0.692396132.3795480

21、.3452153.016304-0.0761187132.3796850.345337153.016327-0.076095圖5 摩阻系數(shù)比較曲線Fig. 5 Comparison curve of frictional coefficient由表3和圖5可知，加噪聲前后，兩種類型管段的摩阻系數(shù)的總體變化值很小，最大差值不超過1，說明基于PSO算法的反分析計算模型的穩(wěn)定性是好的。5 結(jié)束語本文應(yīng)用PSO算法對管道摩阻系數(shù)的反演進(jìn)行了探討，并以實例進(jìn)行了驗證，分析了基于PSO的反演模型的穩(wěn)定性。分析結(jié)果證明了基于PSO的管道摩阻系數(shù)反演計算是可行、簡單且穩(wěn)定的，為給水管網(wǎng)的反分析計算提供了一種新的方法。應(yīng)該指出的是：實際管網(wǎng)是不斷更新改造的，不同埋設(shè)年代的管道由于其工作的起始時間不同，其粗糙程度亦各有不同。反演計算時，應(yīng)該將不同管材、不同年代的管道的粗糙系數(shù)以不同的反演變量代表，這樣不同時間反演量的個數(shù)是不同的，需要根據(jù)具體情況設(shè)定反演問題和模型，而且隨著反演量個數(shù)的增加，監(jiān)測量的個數(shù)也應(yīng)相應(yīng)增加，必須保證獨立的觀測量的個數(shù)反演量的個數(shù)，否則問題將是不適定的。此外由于觀測誤差、測點位置等都對反演結(jié)果有較大影響，這些問題有待深入研究。參考文獻(xiàn)1 許仕榮,邱振華.給水管網(wǎng)的計算理論與電算應(yīng)用M.長沙:湖