基于多目標(biāo)粒子群算法的停泵水錘防護(hù)優(yōu)化

劉亞萌，蔣 勁，李 婷，劉 承，陳英強

(武漢大學(xué)動力與機械學(xué)院 水力機械過渡過程教育部重點實驗室，武漢 430072)

水錘現(xiàn)象一直以來是泵站及輸水工程的常見問題。事故停泵造成的壓力波動會造成管道產(chǎn)生巨大的水錘壓力和汽化現(xiàn)象，嚴(yán)重時將導(dǎo)致管道破裂、設(shè)備損壞等事故，威脅泵站的運行安全[1,2]。通常需要增加必要的水錘防護(hù)措施，以保證泵站系統(tǒng)的安全運行，但由此會帶來工程投資的增加，因此為了解決安全和經(jīng)濟之間的矛盾，對其進(jìn)行優(yōu)化計算是十分必要的[3]。

目前實際工程和科學(xué)研究的對象大多是多目標(biāo)優(yōu)化問題，與單目標(biāo)優(yōu)化不同的是，多目標(biāo)優(yōu)化涉及的參數(shù)較多，目標(biāo)量不唯一，導(dǎo)致在算法選擇和計算工作量上有較大的難度[3-5]。本文計算了國內(nèi)某輸水工程的停泵過渡過程，以水泵出口閥門的關(guān)閥規(guī)律和單向調(diào)壓塔的直徑、初始液位作為設(shè)計變量，以系統(tǒng)最大、最小壓力以及單向調(diào)壓塔的有效體積作為目標(biāo)函數(shù)，建立數(shù)學(xué)模型并采用多目標(biāo)粒子群算法對該工程進(jìn)行優(yōu)化計算，得出了有效可靠且經(jīng)濟的水錘防護(hù)方案。

1 泵站水錘數(shù)值計算

國內(nèi)某輸水工程，全長2.4 km，其中引水隧洞長為347 m，直徑為2.2 m，泵站布置4臺臥式單級、雙吸離心泵，總裝機容量4×2 240 kW，水泵額定轉(zhuǎn)速1 480 r/min，轉(zhuǎn)動慣量為223.24 kg/m2，設(shè)計揚程為164.8 m，設(shè)計流量2.90 m3/s，整個輸水管線距離短、揚程高、流量大且管線布置起伏較大，駝峰點較多，具體的管線布置圖如圖1所示。

圖1 輸水工程管線縱剖面圖Fig.1 Longitudinal section of pipeline

1.1 水錘計算特征線法

泵站水錘的計算是針對整個水泵抽水裝置，包括管道內(nèi)點及與管道連接的泵裝置中的各部分。在水錘計算中，對于管道系統(tǒng)內(nèi)點的計算是求解水錘基本方程，即由運動方程和連續(xù)性方程組成的雙曲型偏微分方程組[6,7]。采用特征線方法將該偏微分方程組離散化，沿特征線方向?qū)⑺D(zhuǎn)換為水錘全微分方程：

(2)

由公式(1)和(2)進(jìn)行有限差分近似，可得編入計算機程序的相容性方程為：

(4)

式中：B=a/gA；R=fΔx/2gDA2。

1.2 單向調(diào)壓塔模型及工作原理

單向調(diào)壓塔與管道中間裝有止回閥，只允許塔中水單向流入主干管中，其構(gòu)造原理圖如圖2所示。水泵正常運行時，由于單向調(diào)壓塔內(nèi)水位低于管中正常工作水頭，止回閥處于關(guān)閉狀態(tài)，向主管道補給水體的流量為零，即Qp3=0；發(fā)生事故停泵后，當(dāng)管道中的壓力低于單向調(diào)壓塔預(yù)先設(shè)定的壓力值時，即Hp

圖2 單向調(diào)壓塔構(gòu)造原理圖Fig.2 Schematic diagram of one way Surge Tower

由連續(xù)性原理可知：

(5)

若任意時刻單向調(diào)壓塔內(nèi)水位HP3，可由此計算時段Δt初的水位H3和流出水體的體積對應(yīng)的水深，求得：

HP3=H3-0.5Δt(QP3+Q3)/Ast

(6)

主管道的相容性特征線方程為：

QP1=(Cp-Hp)/B

QP2=(Cp-Hm)/B

(7)

由公式(5)～(7)，可解出QP3的表達(dá)式為：

QP3=g(CdAp)2-{0.5 (B+Δt/Ast)+{0.25 (B+Δt/Ast)2-

[Cp+Cm+Δt/AstQ3]-2 (H3+ZT)/[g(CdAp)]2}1/2}

(8)

聯(lián)立公式(7)、(8)求得：

(9)

式中：Qp1、Qp2分別為Δt流進(jìn)、流出調(diào)壓塔的流量；Qp3、Q3分別為Δt時段末、初由調(diào)壓塔流向主管道的流量；Hp3、H3分別為Δt時段末、初調(diào)壓塔的水位；Ast為單向調(diào)壓塔的面積；Hp為Δt時段末的壓力水頭。

1.3 無防護(hù)措施事故停泵計算結(jié)果

泵站4臺機組同時正常運行時系統(tǒng)最大壓力出現(xiàn)在水泵出口，壓力值為197.5 m；最不利的運行工況即4臺機組同時事故停泵，無防護(hù)措施兩階段關(guān)閉閥門(4 s快關(guān)至0.3，40 s全關(guān))，最大壓力值為415.5 m，達(dá)到穩(wěn)態(tài)運行最大壓力的2.1倍，引水管及沿線多處出現(xiàn)-8 m的壓力，圖3為正常運行和事故停泵關(guān)閥的管道壓力包絡(luò)線。實踐證明不加防護(hù)措施僅控制泵出口的閥門不能改善水錘壓力過大和負(fù)壓情況，因此設(shè)置單向調(diào)壓塔進(jìn)行防護(hù)，但是單向調(diào)壓塔的尺寸以及閥門的關(guān)閉規(guī)律對系統(tǒng)的影響很大，為使整個系統(tǒng)安全運行并降低工程的投資成本，將對調(diào)壓塔尺寸和閥門的關(guān)閉規(guī)律進(jìn)行優(yōu)化。

圖3 正常運行和事故停泵關(guān)閥的管道壓力包絡(luò)線Fig.3 Pipe pressure envelope for normal operation and shutdown of pump

2 優(yōu)化模型建立

本文采用的多目標(biāo)粒子群優(yōu)化算法(Multi-Objective Particle Swarm Optimization,MOPSO)對停泵水錘計算進(jìn)行優(yōu)化，該算法將原來只能用于單目標(biāo)上的粒子群算法(PSO)經(jīng)過改進(jìn)應(yīng)用于多目標(biāo)優(yōu)化問題中,具體算法如下：

設(shè)粒子群算法中粒子本身找到的最優(yōu)解為個體極值(pBest)，整個群體所經(jīng)歷過的最好位置為全局極值(gBest)，粒子 的速度和位置將按如下公式進(jìn)行更新[11,12]：

Vi=wVi+c1rand() (pBest[i]-Xi)+

c2Rand() (pBest[g]-Xi)

(10)

Xi=Xi+Vid

(11)

式中：c1，c2稱為學(xué)習(xí)因子，均為常數(shù)；rand()和Rand()是[0，1]上的隨機數(shù)；w為慣性權(quán)重。

解空間的粒子根據(jù)式(10)和式(11)不斷調(diào)整自己的位置和速度。一個具有n個設(shè)計變量，m個目標(biāo)的多目標(biāo)優(yōu)化可以描述為：

miny=F(x)=[f1(x),f2(x),…,fm(x)]

(12)

s.t.gi(x)≤0

(13)

式中：x=(x1,x2,…,xm)∈X?Rn為n維決策向量，X為n維決策空間；y=(y1,y2,…,xn)∈Y?Rm為m維的目標(biāo)變量；Y為m目標(biāo)空間；gi(x)≤0是系統(tǒng)約束。

通常情況下，各個目標(biāo)函數(shù)之間可能相互制約甚至是矛盾沖突的，難以使每個目標(biāo)函數(shù)同時達(dá)到最優(yōu)，而多目標(biāo)優(yōu)化的最終目的是在各個目標(biāo)之間進(jìn)行協(xié)調(diào)權(quán)衡和折衷處理，尋找到相對較優(yōu)解[13,14]。多目標(biāo)粒子群優(yōu)化算法可按照如下公式對傳統(tǒng)粒子群算法采用變異策略：

vm=2 (r3-1)βVmax

(14)

xdi(t)=xdi(t)+vm

(15)

式中：vm是變異值；β∈[0,1]為變異系數(shù)，可調(diào)節(jié)變異程度；r3為在[0,1]范圍隨機變化的值；xdi表示第i個粒子的隨機選中的第d維。

經(jīng)過改進(jìn)后，多目標(biāo)粒子群算法可提高解的多樣性，增加粒子的全局搜索能力，且由變異產(chǎn)生的優(yōu)異粒子能夠?qū)ζ渌猱a(chǎn)生吸引，以此增強粒子對局部最優(yōu)的逃逸能力，避免其陷入局部最優(yōu)[7,15]。通常優(yōu)化計算的過程流程圖如圖4所示：首先需要建立模型，明確工程中的設(shè)計變量、優(yōu)化目標(biāo)和約束條件，進(jìn)而確定設(shè)計變量的范圍和參數(shù)初值，在優(yōu)化算法的驅(qū)動下進(jìn)行迭代計算，若目標(biāo)值的計算結(jié)果不滿足工程要求，需要重新給定初值進(jìn)行優(yōu)化，若滿足可輸出結(jié)果。

圖4 優(yōu)化流程圖Fig.4 Optimization flow chart

針對本文計算的系統(tǒng)，在泵站其他參數(shù)均不變的前提下，僅對約束變量進(jìn)行優(yōu)化，選取系統(tǒng)中對水錘防護(hù)影響較大的參數(shù)進(jìn)行約束，其中設(shè)計變量為：單向調(diào)壓塔的直徑、初始液位、閥門快關(guān)時間、快關(guān)百分比、慢關(guān)時間；目標(biāo)變量為：系統(tǒng)最大壓力、出水管最小壓力、引水管最小壓力和單向調(diào)壓塔的有效體積，約束的條件需滿足《泵站設(shè)計規(guī)范》以及本工程的設(shè)計要求，據(jù)此設(shè)計變量、目標(biāo)量和約束的設(shè)定具體如表1和表2所示。

表1 設(shè)計變量和約束條件的設(shè)定Tab.1 Design variables and constraints

3 優(yōu)化計算

表2 預(yù)期優(yōu)化目標(biāo)值Tab.2 expected optimization target value

注：D、H分別為單向調(diào)壓塔的直徑和初始液位。

參照表1和表2列出的設(shè)計變量、目標(biāo)變量和約束條件，在模型優(yōu)化計算前初步試算了20種方案，并提取出其中5種有代表性的結(jié)果，如表3所示。

從方案1、2、3可看出，在不改變調(diào)壓塔尺寸的情況下，增加快關(guān)速率和慢關(guān)時間，系統(tǒng)最大水錘壓力明顯減小且負(fù)壓減輕；對比方案3、4、5可得，關(guān)閥方案相同時，改變單向調(diào)壓塔的直徑和初始液位，調(diào)壓塔的有效體積越大，負(fù)壓情況改善明顯，但同時工程造價也就越高。這5種方案的計算結(jié)果雖能改善最大壓力和負(fù)壓，但均沒有達(dá)到工程要求，因此以方案3的參數(shù)作為初始值，并結(jié)合表1的約束和工程要求，設(shè)定各參數(shù)的范圍進(jìn)行優(yōu)化計算，其中：快關(guān)時間在4～6 s內(nèi)、慢關(guān)時間在40～60 s內(nèi)，慢關(guān)百分比在0.1～0.2內(nèi)(快關(guān)百分比=1-慢關(guān)百分比，故不再單獨設(shè)定)、單向調(diào)壓塔的直徑和初始液位分別在3～5 m和3～7 m內(nèi)。

表3 試算方案的參數(shù)設(shè)定及結(jié)果Tab.3 Parameter setting and results

采用多目標(biāo)粒子群算法進(jìn)行優(yōu)化后，選取快關(guān)時間和調(diào)壓塔初始液位2個設(shè)計變量的優(yōu)化歷史，如圖5所示，圖中兩條直線相交的點即為最優(yōu)的取值。在優(yōu)化過程中，可清晰地看出初始設(shè)計變量的取值遍布了整個變量可行域，說明優(yōu)化算法中的粒子具有全局搜索的能力，隨后逐漸縮小取值范圍，沿較優(yōu)的結(jié)果方向進(jìn)行計算，經(jīng)過600次的迭代后，得到了優(yōu)化目標(biāo)的全局最優(yōu)解。結(jié)果顯示在第392次迭代時獲得最優(yōu)規(guī)律：(水泵在第10 s發(fā)生事故停泵，故快關(guān)和慢關(guān)時間需減去10 s)快關(guān)時間為4.54 s，慢關(guān)百分比為0.10(即4.54 s關(guān)閥至0.1)，慢關(guān)時間為59.96 s，單向調(diào)壓塔的初始液位6.19 m，直徑3.18 m，在該工況下系統(tǒng)最大壓力為247 m，出水管最小壓力-1.53 m，引水管2點處的最小壓力-0.18 m，且單向調(diào)壓塔有效體積為49.16 m3相對較小，滿足運行可靠和經(jīng)濟的要求。

圖6和圖7是最大壓力和最小壓力與調(diào)壓塔尺寸及關(guān)閥參數(shù)的三維圖形。

圖6中最大壓力值與左圖的關(guān)閥參數(shù)在XZ、YZ軸上規(guī)律性較強，相比而言難以提取最大壓力與調(diào)壓塔尺寸之間的關(guān)系，說明在該泵站系統(tǒng)中慢關(guān)時間、關(guān)閥速率對最大壓力的貢獻(xiàn)率較大。分析原因可知，停泵后管道中的水流先由正向流動逐漸減小到零流速，在到達(dá)零流速之前快關(guān)閥門至某一角度，會造成很大的局部阻力，而閥門此時尚未全部關(guān)閉，因此水流在重力作用下倒流會使升壓減小；再逐漸慢關(guān)閥門，使流速變化的增量減小，由于壓力的升高與流速的變化成正比，從而可使管道的升壓限制在允許的范圍之內(nèi)。因此合理控制關(guān)閥速率和慢關(guān)時間，可達(dá)到減小水錘壓力升高的防護(hù)效果。

同理由圖7可知，單向調(diào)壓塔的直徑和初始液位對負(fù)壓的改善較為明顯，而最小壓力與關(guān)閥規(guī)律的曲面圖起伏較大，凸起點較多規(guī)律不明確。單向調(diào)壓塔須具有足夠的容積，以滿足主管道內(nèi)產(chǎn)生負(fù)壓及汽穴空腔時所需的補給水量，保證塔內(nèi)水體不被泄空；其次由于塔內(nèi)水體是靠重力作用向管道補水，故調(diào)壓塔的初始液位需滿足一定的要求，如果水壓不夠，向管道補水不及時則難以達(dá)到破壞水柱分離的效果，因此單向調(diào)壓塔的直徑和初始液位對負(fù)壓有著較大的影響。

圖5 快關(guān)時間和調(diào)壓塔初始液位的優(yōu)化過程Fig.5 The optimization process of the quick closing time and the initial level of the surge tank

圖6 最大壓力與關(guān)閥參數(shù)和調(diào)壓塔尺寸的三維圖Fig.6 3D diagram of maximum pressure between closing valve parameters and size of surge tank

圖7 最小壓力與單向調(diào)壓塔尺寸和關(guān)閥參數(shù)的三維圖Fig.7 3D diagram of the minimum pressure between closing valve parameters and size of surge tank

4 結(jié)果對比

將初始方案和優(yōu)化計算后的結(jié)果進(jìn)行對比，具體見表4，可知優(yōu)化后調(diào)壓塔的有效體積由87.96 m3降至49.16 m3，減小了44.1%；最大壓力值由276.71 m降至247 m，最小壓力值雖有初始的-0.61 m降為-1.53 m，但依然控制在了-2 m以內(nèi)，符合《泵站設(shè)計規(guī)范》以及本工程的設(shè)計要求，這也說明了多目標(biāo)優(yōu)化的解的特點，即該解對一個或幾個目標(biāo)函數(shù)不能進(jìn)一步改進(jìn)，且其他目標(biāo)值不至于劣化，即為尋求到的最優(yōu)解。

由于優(yōu)化方案的設(shè)計變量取值較為精確，實際操作中難以實現(xiàn)，故提出了建議值。建議方案與優(yōu)化計算的方案相比，在關(guān)閥規(guī)律和調(diào)壓塔參數(shù)上做了調(diào)整，更容易泵站的實際操作，且調(diào)整后建議方案的各個目標(biāo)量的結(jié)果均符合設(shè)計要求。

表4 計算結(jié)果對比 Tab.4 Comparison of calculation results

圖8和圖9分別為初始方案和優(yōu)化方案的壓力包絡(luò)線，以及二者調(diào)壓塔液位的變化曲線對比圖，由圖8可明顯看出水錘壓力有較大程度的減小，最大壓力由穩(wěn)態(tài)運行壓力的1.4倍降為1.25倍。圖9中優(yōu)化后單向調(diào)壓塔的初始液位由初始方案的7 m降為6.19 m，初始液位的降低將改變調(diào)壓塔的有效體積，進(jìn)而設(shè)計調(diào)壓塔的實際尺寸時也可適當(dāng)減小，降低了工程投資，達(dá)到了經(jīng)濟、可靠的目的。

圖8 初始方案和優(yōu)化方案的管道壓力包絡(luò)線對比圖Fig.8 Comparison of pipeline pressure envelope with initial scheme and optimization scheme

圖9 初始方案和優(yōu)化方案的調(diào)壓塔液位變化圖Fig.9 Comparison of liquid level of pressure regulating tower in initial scheme and optimization scheme

5 結(jié) 語

本文以國內(nèi)某輸水工程為例，采用多目標(biāo)粒子群算法進(jìn)行優(yōu)化計算，將閥門關(guān)閉規(guī)律和單向調(diào)壓塔初始液位和直徑作為設(shè)計變量，控制調(diào)壓塔體積和最大、最小壓力，尋找到最優(yōu)的方案，且得到以下結(jié)論。

(1)采用多目標(biāo)粒子群算法對本工程進(jìn)行優(yōu)化計算后，尋找到了折衷最優(yōu)解，單向調(diào)壓塔體積由87.96 m3降至49.16 m3，最大水錘壓力由穩(wěn)態(tài)壓力的1.4倍降至1.25倍。

(2)本文結(jié)合約束條件和試算方案中的較優(yōu)值，給定了設(shè)計變量的范圍，從而通過精細(xì)化計算得到更優(yōu)的結(jié)果，為優(yōu)化計算減少了工作量，提高了運算效率。

(3)多目標(biāo)粒子群算法的優(yōu)化結(jié)果與初始方案比較，各目標(biāo)值均有較大的改善，能夠在工程要求和經(jīng)濟投資中尋找到較優(yōu)解，說明該優(yōu)化算法能夠有效地計算水錘問題，且為其他工程實際中的多目標(biāo)優(yōu)化問題提供了解決思路。

□

[1] 王 娜. 泵站壓力管道的水錘研究[J]. 水利規(guī)劃與設(shè)計, 2016,(1):85-88.

[2] 王 圃, 許蘭森, 王 穎. 高揚程多起伏輸水管道水錘模擬及其防護(hù)[J]. 給水排水, 2014,(6):112-114.

[3] 肖曉偉, 肖 迪, 林錦國,等. 多目標(biāo)優(yōu)化問題的研究概述[J]. 計算機應(yīng)用研究, 2011,28(3):805-808.

[4] 謝 濤, 陳火旺. 多目標(biāo)優(yōu)化與決策問題的演化算法[J]. 中國工程科學(xué), 2002,4(2):59-68.

[5] DEB K,PRATAP A,AAGRWALI S,et al. A fastand elitist multiobjective genetic algorithm: NSGA-II[J]. IEEE Trans on Evolutionary Computation,2002,6(2):182-197.

[6] 蘭 剛, 蔣 勁, 李東東,等. 空氣罐在長距離輸水管線水錘防護(hù)中的應(yīng)用[J]. 水電站機電技術(shù), 2013,(4):42-44.

[7] 楊祖強, 羅 希, 王文文. 基于改進(jìn)型多目標(biāo)遺傳算法的含氣水錘防護(hù)方案優(yōu)化研究[J]. 水資源與水工程學(xué)報, 2014,25(1):200-204.

[8] 劉梅清, 劉志勇, 蔣 勁. 基于遺傳算法的單向調(diào)壓塔尺寸優(yōu)化研究[J]. 中國給水排水, 2008,24(23):56-60.

[9] 徐艷艷. 長距離高揚程多起伏輸水管道采用箱式雙向調(diào)壓塔等措施的水錘防護(hù)研究[D]. 西安：長安大學(xué), 2008.

[10] 劉光臨, 劉志勇, 王聽權(quán),等. 單向調(diào)壓塔水錘防護(hù)特性的研究[J]. 給水排水, 2002,28(2):82-85.

[11] SIERRA M R, COELLO C A C. Improving PSO-based multi-objective optimization using crowding, multation and Edominance[C]∥ Lecture Notes in Computer Scence. Proc of the 3rd International Conference on Evolutionary Multi-Criterion Optimization,2005:505-519.

[12] 劉衍民, 牛 奔, 趙慶禎. 多目標(biāo)優(yōu)化問題的粒子群算法仿真研究[J]. 計算機應(yīng)用研究, 2011,28(2):458-460.

[13] TSAI S J,SUN T Y,LIU C C. An improved multi-objective particle swarm optimizer for problems[J].Expert Systems with Applications,2010,37:872-886.

[14] 畢榮山, 楊 霞, 譚心舜,等. 基于動態(tài)Pareto解集的微粒群優(yōu)化算法及其在多目標(biāo)規(guī)劃中的應(yīng)用[J]. 計算機工程與應(yīng)用, 2004,40(32):85-88.

[15] 劉衍民, 牛 奔, 趙慶禎. 多目標(biāo)優(yōu)化問題的粒子群算法仿真研究[J]. 計算機應(yīng)用研究, 2011,28(2):458-460.