序列二次規(guī)劃法在多水源管網(wǎng)優(yōu)化調(diào)度中的應(yīng)用研究

申請免費試用、咨詢電話：400-8352-114

簡介： 供水管網(wǎng)優(yōu)化調(diào)度的一級優(yōu)化是一個非線性優(yōu)化問題，本文通過分析管網(wǎng)的水力關(guān)系，對管網(wǎng)水力關(guān)系進行合理的線性化，使目標(biāo)函數(shù)和約束條件顯式化，將問題轉(zhuǎn)化為序列二次規(guī)劃問題。在求解二次規(guī)劃問題中，考慮到大部分節(jié)點水頭的約束是非作用約束，利用線性化結(jié)果，將非作用約束從約束集中剔除，同時將齒行法的思想和水力學(xué)上的基本概念相結(jié)合，提出了一種適合本問題的修正齒行法，將二次規(guī)劃結(jié)果拉回到原約束面，保證了解的可行性。最后還初步分析了優(yōu)化計算的計算量。數(shù)值試驗表明本文的方法計算量小、效率高，結(jié)果可靠。
關(guān)鍵字：多水源管網(wǎng) 優(yōu)化 序列二次規(guī)劃 齒行法

(1.浙江大學(xué) 水工結(jié)構(gòu)與水環(huán)境研究所，浙江 杭州 310027；2.杭州市自來水總公司，浙江 杭州  310016)

1 問題的提出

    城市供水管網(wǎng)是城市的生命線之一。這一復(fù)雜的網(wǎng)絡(luò)系統(tǒng)，主要通過幾個供水泵站為城市血液提供能量，送至城市的各個角落。通過對供水泵站的優(yōu)化調(diào)度，可以降低企業(yè)的制水成本、使管網(wǎng)的供水壓力分布更合理，根據(jù)初步估計，對于一個日供水量為10萬噸的自來水公司，如果供水揚程降低1m，每年可以節(jié)電15萬kW·h；由于管道系統(tǒng)的滲漏與水頭有關(guān)，降低供水水頭也可以在一定程度上減少管網(wǎng)的滲漏；供水水頭的降低還可以減少爆管的風(fēng)險，這對于管網(wǎng)的管理有更深刻的意義。因此管網(wǎng)合理調(diào)度研究一直是供水企業(yè)一個重要課題，同時也是一個難題。
    供水管網(wǎng)運行的合理調(diào)度可以用一個最優(yōu)化問題來描述，管網(wǎng)的水力方程組是一組非線性方程，各水源水泵的開啟狀態(tài)作為離散變量，因此這是一個混合變量的非線性最優(yōu)化問題。由于離散變量與連續(xù)變量的同時存在，求解極為不方便，最常用的方法是將該問題分作兩級進行優(yōu)化：一級優(yōu)化是針對管網(wǎng)而言，目的在于求各水源的最佳供水量或最佳供水揚程；二級優(yōu)化是在一級最優(yōu)化的基礎(chǔ)上，根據(jù)水源的具體情況，確定滿意的水泵開啟方案和水泵的調(diào)速比。采用以上方法可以在一定程度上降低求解的困難，但一級優(yōu)化也是一個非線性的優(yōu)化問題，求解起來相當(dāng)麻煩，目前國內(nèi)外最常用的方法是廣義簡約梯度法。廣義簡約梯度法雖屬較優(yōu)秀的約束非線性規(guī)劃算法，根據(jù)作者在以往其他優(yōu)化應(yīng)用方面的研究，其重分析次數(shù)相當(dāng)多。在本優(yōu)化問題中水力計算是計算量的主體部分，由數(shù)值試驗的經(jīng)驗知，在目前中等配置的微機上完成一個2000個左右節(jié)點的供水管網(wǎng)，一次水力計算需要10s左右，如果采用廣義簡約梯度法，需要反復(fù)迭代計算，花費的時間是相當(dāng)可觀的。由此可見采用廣義簡約梯度法實現(xiàn)管網(wǎng)的在線優(yōu)化調(diào)度存在較大難度。
    本文針對一級優(yōu)化問題，采用序列二次規(guī)劃法進行求解。M.J.D.Powell所給出的序列二次規(guī)劃法實質(zhì)上是運用Kuhn?Tucker最優(yōu)化條件所形成的非線性方程進行迭代計算，而這一迭代過程恰好可以用求解一相應(yīng)的二次規(guī)劃問題替代，故原問題的求解過程轉(zhuǎn)化為求解一個二次規(guī)劃的序列。其中二次規(guī)劃問題的二次目標(biāo)函數(shù)是原問題Lagrange函數(shù)的二次展開式，包含了目標(biāo)與約束函數(shù)的二次信息。通常其二階導(dǎo)數(shù)矩陣由變尺度的思想通過先前迭代點的梯度信息逐步生成。序列二次規(guī)劃法綜合利用了K?Ｔ條件、變尺度、線性及二次近似等有效手段，在理論上是一種比廣義簡約梯度法優(yōu)秀的算法^[1]，但是它的迭代序列通常從不可行域逐步逼近可行域，需要在極限情況下才能完全達到約束要求，這顯然不利于盡快獲得可行的較優(yōu)解，故約束條件的妥善處理非常重要，本文將結(jié)構(gòu)優(yōu)化中齒行法的思想和水力學(xué)的基本概念相結(jié)合，提出了一種新的算法，可以方便地將迭代中的非可行點拉回到約束界面上，獲得了較高的計算效率，有助于實現(xiàn)管網(wǎng)的在線優(yōu)化調(diào)度。

2  供水優(yōu)化調(diào)度一級優(yōu)化的數(shù)學(xué)模型

    管網(wǎng)的運行調(diào)度一般以經(jīng)濟性作為目標(biāo)函數(shù)，與水源的供水量、供水水頭有關(guān)，據(jù)此可以建立供水管網(wǎng)的目標(biāo)函數(shù)：

min·FG(Q_s,H_s)

(1)

式中：FG為各水源的制水成本和供水的動力費用；Q_s、H_s為各水源的供水量和供水水頭。
    供水調(diào)度的主要約束條件有：管網(wǎng)的水力關(guān)系，各水源的水量和水壓的約束，管網(wǎng)中各節(jié)點的最小服務(wù)水頭。這些約束條件分別表示如下：
    管網(wǎng)水力關(guān)系

F(H_s,H_N,Q_N) =0

(2)

    各水源的供水水頭約束

H_smin≤H_s≤H_smax 

(3)

    各水源的供水量約束

Q_smin(H_s)_≤Q_s≤Q_smax(H_s) 

(4)

    管網(wǎng)各節(jié)點服務(wù)水頭約束

(5)

    其中：H_smax、H_smin分為水源的最大、最小供水水頭；Q_smin(H_s)、Q_smax(H_s)分為水源的最大、最小供水能力，通常水源的供水量的能力與供水水頭有關(guān)。H_N為管網(wǎng)中各節(jié)點的服務(wù)水頭；H_Nmax、H_Nmin分為管網(wǎng)中各節(jié)點的最大、最小服務(wù)水頭；Q_N為管網(wǎng)中各節(jié)點的節(jié)點流量。

3  模型的求解
    模型求解主要有2個難點：(1)約束條件太多，一個中等復(fù)雜的城市管網(wǎng)可能會有上千個約束；(2)目標(biāo)函數(shù)中各變量隱式相關(guān)——水源的供水水頭H_s和供水水量Q_s隱式相關(guān)。如果能對以上兩個方面進行適當(dāng)?shù)奶幚?，可以大大的降低難度，提高求解效率。針對以上兩點，本文從管網(wǎng)的水力條件出發(fā)，提出了一套求解方法：
    在一定負荷Ｑ_N下，將管網(wǎng)的水力計算公式(2)在H₀處作一階泰勒展開有：

(6)

, 稱為敏度矩陣。如果用哈真-威廉公式表示管道的能量損失，用矩陣A、B可以分別表示為 , ，管網(wǎng)的水力學(xué)公式可以用式(7)表 達。A和B僅與管網(wǎng)中管道的水力坡度有關(guān)。當(dāng)任一水源的供水水頭發(fā)生變化，由于管網(wǎng)自身的調(diào)節(jié)作用，每根管道的水力坡度的變化幅度要比節(jié)點水頭變化小得多，A、B的變化都比較小。管網(wǎng)的水力計算公式(2)在H0附近可以線性近似為式(7)，且方程有足夠的精度(算例的數(shù)值計算結(jié)果參見附錄)。

(7)

在文獻^[2]中已證明B是正定對稱矩陣，其逆矩陣存在。

令  

則

(8)

    矩陣C的分量c_i,j反映了第j個水源對節(jié)點i的影響，矩陣C也稱為影響矩陣。如果管網(wǎng)中所有水源的供水水頭同步上升Δh，即ΔH_N＝[Δh1,Δh2,…,Δhs]^T，相當(dāng)于管網(wǎng)的參考水位提高了Δh。由式(8)知管網(wǎng)中任一點的水頭上升的水位 ，因此其中矩陣C的行向量的各分量之和必等于1，各管段的水力坡度不變。如果各水源的供水水頭和節(jié)點流量已知，可求得管網(wǎng)中的各節(jié)點的水頭，同樣可以求出各水源的供水量。水源泵站供水的動力能耗可以表示為Q_s·H_s/γ(γ為水源效率)，供水的動力費用與耗能成正比。水源供水量在H⁰附近可以線性近似為：Q_s=L_sH_s。水源的制水0費用(除動力費用)可以表示成R_sL_sH_s(R_s表示各水源的單位制水成本)，因此目標(biāo)函數(shù)在H⁰處可以近似用水源水頭的二次函數(shù)表示如下： 

(9)

各水源的供水量約束在H⁰處可以線性近似表示： 

(10)

    由于在管網(wǎng)中往往只是一部分的最不利節(jié)點違反約束，只要最不利的節(jié)點滿足了服務(wù)水頭的要求，其他節(jié)點也滿足了要求，因此可以將最不利的一些節(jié)點與水源節(jié)點的水頭關(guān)系從式(9)中的影響矩陣C中抽取出來，表示成矩陣G，管網(wǎng)節(jié)點水頭的約束方程(5)可以簡化表示如下： 

(11)

其中Ω為最不利節(jié)點的集合。
    通過上述方法，一級優(yōu)化模型在H⁰附近可以近似表示為線性約束的二次規(guī)劃問題： 

(12a)

_{s.t .    Hsmin≤Hs≤Hsmax} 

(12b)

_{Qsmin(H0s) +KsminHs≤LsHs≤Qsmax(H0s) +KsmaxHs} 

(12c)

(12d)

    在原優(yōu)化問題中，各水源的供水量、管網(wǎng)中節(jié)點的水頭是水源供水水頭的函數(shù)，是隱式關(guān)系，求解起來非常不方便。通過把管網(wǎng)水力關(guān)系式(2)線性化，消去原目標(biāo)函數(shù)(1)中的變量——水源供水量Q_s，可以把目標(biāo)函數(shù)表示僅含水源水頭變量的形式，將管網(wǎng)中各節(jié)點的水頭H_N表示成水源的供水水頭H_s的線性函數(shù)，只取其中最不利一部分作為每次優(yōu)化計算的約束條件，這樣大大地減少了約束條件。如果管網(wǎng)有上千個節(jié)點，只要保證最不利的10%左右節(jié)點滿足服務(wù)水頭約束，就能基本上保證每次優(yōu)化計算結(jié)果不會離約束邊界太遠，同時優(yōu)化計算的計算量成倍的減少。

    由于采用了線性近似的方法簡化約束條件和目標(biāo)函數(shù)，采用二次規(guī)劃法(QP法)優(yōu)化之后會導(dǎo)致結(jié)果越過實際約束邊界，其中主要是最不利點不滿足管網(wǎng)最小服務(wù)水頭的要求。在結(jié)構(gòu)優(yōu)化設(shè)計中經(jīng)常采用齒行法進行優(yōu)化迭代，其基本思想是在每次優(yōu)化迭代后，通過射線步(即將所有設(shè)計變量以同一倍數(shù)放大或縮小)將結(jié)果拉到最嚴格的約束邊界上。根據(jù)管網(wǎng)水力學(xué)，所有的水源的供水水頭同時都提高或降低相同的水位，使管網(wǎng)的最不利點的水位恰好處于約束邊界上，不會改變各個水源的供水關(guān)系。利用這一特性，可以構(gòu)造一修正的射線步，能夠方便的將中間優(yōu)化迭代點拉回到約束界面上(見圖1)。由于在這一修正的射線步中，每個水源提高的水位相同，因此變化后的值與變化前的值與坐標(biāo)軸正好成45°。上述優(yōu)化方法的流程圖如圖2所示。 

4  算 例 

    為了檢驗上述模型的可行性和可靠性，本文采用Fortran語言編制了優(yōu)化計算程序，算例中管網(wǎng)基本形狀如圖3所示。管網(wǎng)有兩個水源，12個節(jié)點，19根管道，管網(wǎng)節(jié)點的基本數(shù)據(jù)見表1。
    在計算中，不失一般性，以泵站供水的能耗作為目標(biāo)函數(shù)，即 (γi是水源泵站的工作效率)，如果考慮取水、凈化等制水成本的話，目標(biāo)增加一個線性項或一個二次項，對本文方法沒有影響。收斂準(zhǔn)則采用目標(biāo)值｜obj^k-obj^k+1｜/obj^k≤0.0001和前后兩次迭代的步長‖H^k+1-H^k‖≤0.1。為了測試計算方法的有效性和計算效率，本文對多種情況下進行了測試，試驗1和試驗2是比較在不同效率下，優(yōu)化獲得的目標(biāo)值是否能夠充分反映目標(biāo)函數(shù)對計算結(jié)果的影響。由于最大供水量和最大供水揚程往往是制約管網(wǎng)供水能力的約束條件，試驗3、試驗4是檢驗最大供水量約束與最大揚程是作用約束時計算結(jié)果是否可靠。從表2中的計算結(jié)果可以看出計算結(jié)果在以上情況下都能獲得較滿意的結(jié)果。
    優(yōu)化迭代計算量的分析：在計算過程中，計算量主要包括以下部分：1)水力計算；2)約束關(guān)系的線性化；3)二次規(guī)劃問題的求解。水力計算是求解非線性方程組，一般是通過線性化，把非線性方程轉(zhuǎn)化為求解線性方程組。在文獻^[2]中已經(jīng)證明該線性后方程組的系數(shù)矩陣主對角占優(yōu)，對稱正定，可以采用迭代法、高斯消去法、LU分解法或喬斯基分解法進行計算。數(shù)值試驗表明對于較大規(guī)模的管網(wǎng)采用迭代法具有一定的優(yōu)勢。關(guān)于約束關(guān)系的線性化，如果是采用LU分解法進行的水力計算，只需要進行少量的矩陣變換就可獲得約束關(guān)系的線性表示。

表1管網(wǎng)基本數(shù)據(jù)

編號

節(jié)點流量/(L/S)

地面標(biāo)高/m

最小服務(wù)水頭/m

1
2
3
4
5
6
7
8
9
10
11
12
水源1
水源2

36.2
36.8
82.5
36.4
48.7
81.5
198.7
66.1
50.6
43.2
105.8
35.5
—
—

34.0
35.0
36.0
35.0
37.0
41.0
40.0
43.0
42.0
45.0
46.0
48.0
33
30

24.0
24.0
24.0
24.0
24.0
24.0
24.0
24.0
24.0
24.0
24.0
24.0
—
—

    如果是采用迭代法求解線性方程組，那么還需要額外求解S+1個線性方程組(S表示水源數(shù))，由于水源數(shù)一般比較少，約束線性化的計算量相比進行一次水力計算來講要少得多。再有是二次規(guī)劃問題的求解，影響二次規(guī)劃問題的計算量主要是約束量，通過線性化之后，剔除了大部分非作用約束，即使對于上千個節(jié)點的管網(wǎng)，也可以把約束量控制在100個左右。在本文算例的迭代計算過程中，服務(wù)水頭約束條件的最大的違反值是個別最不利點比最小服務(wù)水頭約束小0.45m，這一誤差能夠滿足實際應(yīng)用的要求。二次規(guī)劃問題采用Lemke算法，該算法是在單純形法基礎(chǔ)上加以適當(dāng)修改來求解二次規(guī)劃的K-T點，其計算量與單純行法計算量相當(dāng)，所以二次規(guī)劃計算中的計算量也比較小。綜合以上的分析，可以估計每次尋優(yōu)迭代的計算一般不會超過2～3次水力計算，比常用的廣義簡約梯度法計算量要小得多。本文中的算例特意在配置較低的計算機(PI-120)上進行測試，在以上所有算例的計算時間均少于1s，由此可見本文方法的效率相當(dāng)高。

表2 優(yōu)化約束條件與優(yōu)化結(jié)果

試驗1

試驗2

試驗3

試驗4

水源1

水源2

水源1

水源2

水源1

水源2

水源1

水源2

工作效率
最大揚程約束
最小揚程約束
最大供水量約束
最小供水量約束
供水水位
供水量

0.7
80
60
0.8
0.02
98.48
0.586

0.7
80
60
0.5
0.02
100.95
0.235

0.6
80
60
0.8
0.02
96.90
0.573

0.7
80
60
0.5
0.02
104.87
0.249

0.6
80
60
0.8
0.02
100.43
0.602

0.7
80
60
0.22
0.02
96.83
0.220

0.6
80
60
0.8
0.02
100.87
0.605

0.7
66
60
0.23
0.02
95.99
0.217

目標(biāo)值
迭代次數(shù)

78.83
11

87.76
7

88.66
4

88.66
4

 注：試驗1和試驗2水源水頭初值采用附表1中的狀態(tài)1，試驗1和試驗2水源水頭初值采用附表1中的狀態(tài)2。 

5 小 結(jié)
    供水管網(wǎng)優(yōu)化調(diào)度的一級優(yōu)化問題是一個復(fù)雜的非線性優(yōu)化問題，針對求解過程中的難點，本文通過分析管網(wǎng)的水力關(guān)系，對管網(wǎng)水力關(guān)系進行了合理的線性化，使目標(biāo)函數(shù)和約束條件顯式化，將問題轉(zhuǎn)化為序列二次規(guī)劃問題。在求解二次規(guī)劃問題中，同時考慮到大部分節(jié)點水頭的約束是非作用約束，利用線性化結(jié)果，將非作用約束從約束集中剔除，大大地減少了求解的工作量，并借鑒結(jié)構(gòu)優(yōu)化中齒行法的思想，結(jié)合水力學(xué)上的一些基本概念，提出了適合于本文的修正齒行法，將二次規(guī)劃結(jié)果拉回到約束面，保證了解的可行性。最后本文還簡要分析了優(yōu)化迭代的計算量，數(shù)值試驗表明采用本文的方法計算量小、效率高，結(jié)果可靠，為大規(guī)模管網(wǎng)的實時優(yōu)化調(diào)度提供了良好的基礎(chǔ)。

參 考 文 獻： 

    [1] 汪樹玉，楊德銓，劉國華，張科鋒.優(yōu)化原理、優(yōu)化方法與工程應(yīng)用[M].杭州：浙江大學(xué)出版社，1991，10：375-381，221-228.
    [2] 程偉平，包志仁.管網(wǎng)恒定流水力計算的最優(yōu)化方法[J].水利學(xué)報，2000，7：86-90.
    [3] 楊欽，嚴煦世.給水工程[M].北京：中國建筑工業(yè)出版社，1987.

附錄：

    本文在進行線性化推導(dǎo)過程中，認為矩陣C和D在計算過程中變化比較小，本文在計算中選取了供水水頭差異較大的情況比較了矩陣C和向量E的值(見附表1)，計算結(jié)果表明即使在2種供水狀態(tài)差異較大的情況下C和向量E的變化比較小。表中 ，與理論分析有一點差別，主要是由計算中截斷誤差引起。
    同時本文還比較了計算了在不同的比例負荷作用下，(即各節(jié)點的流量同時乘以一個相同的比例系數(shù))，此時敏度矩陣C矩陣仍然保持基本不變，矩陣假設(shè)基本合理。向量E的各系數(shù)與比例系數(shù)的1?852次方成正比，數(shù)字計算中間結(jié)果見附表2。 

附表１ 不同狀態(tài)下C、E矩陣的變化情況

狀態(tài)1
水源1供水水頭：100.227m
水源2供水水頭： 97.227m

狀態(tài)2
水源1供水水頭： 98.402m
水源2供水水頭：101.131m

1
2
3
4
5
6
7
8
9
10
11
12

C_i,1
0.90448
0.90326
0.66842
0.04095
0.94500
0.89536
0.66775
0.37115
0.91821
0.86529
0.73265
0.66456

C_i,2
0.03544
0.09666
0.33150
0.95897
0.05491
0.10456
0.33217
0.62877
0.08171
0.13463
0.26726
0.33534

E_i
-4.3374
-9.9901
-22.4775
-2.0184
-6.8753
-11.7309
-24.0387
-17.4159
-11.4742
-16.0033
-26.5912
-27.2211

C_i,1
0.96663
0.90943
0.69408
0.04145
0.94833
0.90198
0.69326
0.37706
0.92307
0.87392
0.07514
0.69799

C_i,2
0.03329
0.09048
0.37706
0.95837
0.05158
0.09793
0.30663
0.62293
0.07684
0.12598
0.24847
0.30189

E_i
-4.3379
-9.9913
-22.4826
-2.0183
-6.8760
-11.7323
-24.0439
-17.4166
-11.4752
-16.0051
-26.5951
-27.2264

附表2 不同比例負荷狀態(tài)下C、E矩陣的變化情況 

?狀態(tài)1
水源1供水水頭：100.227m，水源2供水水頭：97.227m，
負荷比Ｑ1/Q0=1.2，且(Q₁/Q₀)^1.852＝1.4017

狀態(tài)2
水源1供水水頭：100.227m，水源2供水水頭：97.227m，
負荷比Ｑ₂/Q₀=1.5，且(Q₂/Q₀)^1.852＝2.1189

1
2
3
4
5
6
7
8
9
10
11
12

C_i,1
0.96482
0.90425
0.67251
0.04103
0.94554
0.89642
0.67181
0.37203
0.91899
0.86668
0.73567
0.66970

C_i,2
0.03510
0.09567
0.32741
0.95889
0.05439
0.10350
0.32811
0.62789
0.08093
0.13324
0.26425
0.33022

E_i
-6.0807
-14.006
-31.5198
-2.8293
-9.6385
-16.4462
-33.7080
-24.4138
-16.0853
-22.4356
-37.2818
-38.1722

Ei/E_i^0*
1.4019
1.4020
1.4018
1.4019
1.4019
1.4020
1.4022
1.4018
1.4018
1.4019
1.4020
1.4048

C_i,1
0.96511
0.90506
0.67588
0.04109
0.94600
0.89727
0.67517
0.37278
0.91963
0.86782
0.73815
0.67400

C_i,2
0.034815
0.094856
0.32404
0.95883
0.053944
0.10263
0.32475
0.62714
0.080289
0.13210
0.26177
0.32592

E_i
-9.1932
-21.1757
-47.66.9
4.2773
-14.5721
-24.8649
-50.9687
-36.9088
-24.3185
-33.9200
-56.3678
-57.7025

Ei/E_i^0*
2.1195
2.1197
2.1204
2.1192
2.1151
2.1196
2.1203
2.1193
2.1194
2.1195
2.1198
2.1204

注：其中E⁰_i采用附表1中狀態(tài)1的計算值

    綜上所述，管網(wǎng)中各節(jié)點在比例負荷的情況下，管網(wǎng)中任一節(jié)點的水頭可以線性近似如下：其中：S為所有水源的集合，K_i為與管網(wǎng)總流量Q₀相對應(yīng)的系數(shù)，其中。 

收稿日期：2002-01-04
作者簡介：劉國華(1963-)，男，福建莆田人，教授，博導(dǎo)，主要從事水工結(jié)構(gòu)優(yōu)化與系統(tǒng)分析方面的研究。