張曉波，鄭雄偉，陳昌軍

（浙江省水利水電勘測設(shè)計(jì)院，浙江 杭州 310002）

隨著計(jì)算機(jī)技術(shù)的廣泛引用，水電站徑流調(diào)節(jié)計(jì)算方法已由傳統(tǒng)的列表試算法、圖解法轉(zhuǎn)向數(shù)值計(jì)算方法。大連理工大學(xué)于1981年首次提出采用改進(jìn)歐拉法進(jìn)行數(shù)值求解，受此啟發(fā)，筆者將水電站徑流調(diào)節(jié)計(jì)算的數(shù)值解法進(jìn)行拓展，分別采用二分法、改進(jìn)歐拉法、不動點(diǎn)迭代法、牛頓迭代法對水電站徑流調(diào)節(jié)進(jìn)行數(shù)值求解推導(dǎo)，并從迭代效率、計(jì)算精度、實(shí)用性等方面進(jìn)行比較，推薦實(shí)用、簡便的不動點(diǎn)迭代法應(yīng)用于徑流調(diào)節(jié)計(jì)算。

1 計(jì)算原理

水電站徑流調(diào)節(jié)基于以下基本方程[1]：

式（1）為水量平衡方程，Q(t)為入庫流量；q(t)為下泄流量；V(t)為水庫庫容。為便于公式表示，采用三次多項(xiàng)式擬合庫水位和庫容關(guān)系[1]，如式（3）：

需要指出，地形復(fù)雜的水庫其水位庫容曲線很難用式（3）完整地表達(dá),很多時(shí)候需采用分段函數(shù)擬合，程序計(jì)算應(yīng)采用插值法。

式（2）為水電站出力計(jì)算公式，N(t)為出力；Z上為庫水位；Z下為電站尾水位；k0為出力系數(shù)，k0=9.81η，η為機(jī)組出力效率系數(shù)，可取一定值；DH(t)為水頭損失，為便于公式表示，表達(dá)為q(t)的二次函數(shù)，如式（4）：

需要指出，實(shí)際發(fā)電中，機(jī)組出力效率是變化的，水頭損失也受出力不均勻的影響，式（4）不能反映這些影響，程序計(jì)算應(yīng)考慮這些具體問題。

電站尾水位也可擬合為q(t)的n次函數(shù)，如式（5）：

聯(lián)解式（1），（2）及相關(guān)公式，顯然難以直接求解。對該非線性方程，列舉以下方法求解之。

2 數(shù)值解法

2.1 改進(jìn)歐拉法

（1）數(shù)學(xué)原理。歐拉法迭代原理如圖1，對于方程 y=y(x)，可以從初始點(diǎn) P(x0，y0)出發(fā),沿該點(diǎn)切線方向推進(jìn)到x=x1上一點(diǎn)P1，然后再從P1依該點(diǎn)切線方向推進(jìn)到x=x2上一點(diǎn)P2依次類推[2]，記y′(xk)=f(xk，y(xk))，有求解關(guān)系：

歐拉法相當(dāng)粗略，多次推進(jìn)后，歐拉法計(jì)算的Pn顯著偏離實(shí)際曲線，因此提出了改進(jìn)歐拉法[2]：

（2）應(yīng)用求解。聯(lián)立式（1），（2），并代入輔助方程式（3），（4），（5），使得未知數(shù) q(t)和 Z上(t)顯現(xiàn)（以下 Z上(t)簡化為 Z(t)），式（1）可變形為：

圖1 歐拉迭代法

式（2）可變形為：

將式（10）代入式（9），可得：

式（8）為非線性的微分方程。大連理工大學(xué)陳守煜教授將式（11）看作常微分方程的初值問題求解[3]，采用改進(jìn)歐拉法求解，陳守煜教授對式（11）的改進(jìn)歐拉求解公式為：

應(yīng)用式（12），（13）可以對水電站徑流調(diào)節(jié)直接求解。q=f-1(Z)無法直接求解，需要迭代法求解。

（3）解法評價(jià)。改進(jìn)歐拉法使得所求的末時(shí)段的Z(t)顯式表示，無需迭代，方法簡便。但是，正是由于改進(jìn)歐拉法是非線性方程的二階近似，只有二階精度。對于在步長△t較大時(shí)（如逐月徑流調(diào)節(jié)），誤差較大，陳守煜教授在其文獻(xiàn)[3]中的示例，已經(jīng)說明了該法精度不及傳統(tǒng)的迭代法，因?yàn)閭鹘y(tǒng)的迭代法可控制誤差精度,歐拉法正是以犧牲精度來實(shí)現(xiàn)其簡單的計(jì)算方法。而且，發(fā)電流量q(t)仍需迭代求解，總體上效率提高有限。

2.2 二分法

（1）數(shù)學(xué)原理。二分法是解非線形方程最常用的一種方法。其原理如圖2。在有根區(qū)間[a，b]，取中點(diǎn)x0=(a+b)/2，假設(shè)中點(diǎn)x0不是f(x)的零點(diǎn)，進(jìn)行根的搜索,若f(x0)與f(a)同號,則所求根x*必在x0的右側(cè)，令 a1=x0，b1=b；否則 x*必在 x0的左側(cè),這時(shí)令 a1=a，b1=x0，…，如此迭代下去，最終逼近 x*[2]。

圖2 二分迭代法

（2）應(yīng)用求解。將式（1）離散化，得：

將式（3）代入式（14）可得：

將式（2），（4），（5）整理得：

將式（16）代入式（15），可得：

式（17）為僅包含未知量 qt＋1的非線形方程。因此，水電站徑流調(diào)節(jié)計(jì)算可轉(zhuǎn)化為非線性方程的求解問題。

為求解式（17）：F(qt＋1)=0?？梢粤?a0=0，b0=qm，qm為最大過機(jī)流量。

（3）解法評價(jià)。二分法算法簡單，只要給定合理的根區(qū)間，迭代且總是收斂的，且容易控制精度。缺點(diǎn)是由于尋根區(qū)間范圍較大，收斂速度較慢。在實(shí)際應(yīng)用計(jì)算中，除了采用式（17）的流量作為迭代因子外，也可采用水位作為迭代因子，迭代精度亦可采用流量、水位、出力任一指標(biāo)控制，只需應(yīng)用其原理即可。

2.3 不動點(diǎn)迭代法

（1）數(shù)學(xué)原理。不動點(diǎn)迭代的基本思想是：對于 f(x?)=0，可以改成等價(jià)形式 x=φ(x)。若要求 x*滿足 f(x?)=0，則 x?=φ(x*)；反之亦然。稱 x*為函數(shù) φ(x)的一個(gè)不動點(diǎn)。建立迭代關(guān)系：xk＋1=φ(xk)，k=0，1，2，…。

其幾何意義如圖3，方程x=φ(x)的求根問題就是確定直線y=x于曲線y=φ(x)的交點(diǎn)P?。給定初值 x0，求得 y0=φ(x0)，再令 x1=φ(x0)，求得 y1=φ(x1)，再令 x2=φ(x1)……，直至|xk＋1-xk|＜ε[2]。

圖3 不動點(diǎn)迭代法

（2）應(yīng)用求解。將式（15）等價(jià)改寫為：

由此可以建立迭代關(guān)系式（20），（21）：

（上標(biāo) k 為迭代次數(shù)，k=0，1，2，…）

（3）解法評價(jià)。同二分法相比，不動點(diǎn)迭代法具有方法簡便，收斂速度快的優(yōu)勢。本文在第3節(jié)給予實(shí)例說明。需要指出的是：不動點(diǎn)迭代法要求有收斂條件：對于x=φ(x)，存在x*的某個(gè)鄰域連續(xù) R:|x-x*|＜δ，且|φ′(x)|＜1，則不動點(diǎn)迭代局部收斂。

水電站徑流調(diào)節(jié)的迭代函數(shù)式（19）求導(dǎo)較為復(fù)雜，難以直接判斷，因此在直接采用式（20），（21）迭代時(shí)有時(shí)會遇到不收斂的情況。將x=φ(x)化為等價(jià)形式：x=(1-ω)x＋ωφ(x)[4]，ω 為松弛因子。于是有迭代公式：

于是對水電站徑流調(diào)節(jié)迭代公式有：

可以證明，取ω=0.5時(shí)，可使得式（23）迭代無條件收斂[4]。

因此，采用改進(jìn)的不動點(diǎn)迭代收斂公式式（23），可使得迭代條件同二分法一樣簡單，但收斂要快很多，同時(shí)可控制計(jì)算精度。與二分法類似，在實(shí)際應(yīng)用計(jì)算中，除了采用式（23）的流量作為迭代因子外，也可采用水位作為迭代因子，迭代精度亦可采用流量、水位、出力任一指標(biāo)控制，只需應(yīng)用其原理即可。2.4牛頓迭代法

（1）數(shù)學(xué)原理。牛頓迭代法的基本思想為：方程f(x)=0的根x?可解釋為曲線y=f(x)與x軸的交點(diǎn)的橫坐標(biāo)。設(shè)xk是根x?的某個(gè)近似值，過點(diǎn)Pk=xk，f(xk))作 y=f(x)的切線，與 x 軸相交于 xk＋1，如此逐步迭代，逼近x?[2]。其迭代方程為：

（2）應(yīng)用求解。對于式（15）建立迭代關(guān)系為：

聯(lián)立式（15），（24），（25），（26）可求解。

圖4 牛頓迭代法

需要指出的是：牛頓迭代法收斂速度很快，但對初值要求很高，若初值偏離真實(shí)解很遠(yuǎn)，牛頓法可能發(fā)散。因此，為保證收斂，在式（21）增加下山因子 λ(0＜λ≤1)，式（24）改寫為：

迭代計(jì)算可從λ=1開始，逐次將λ減半，直至|F(xk＋1)|＜|F(xk)|，可保證迭代收斂。

（3）解法評價(jià)。牛頓迭代法的優(yōu)點(diǎn)是收斂快，可以說是幾種數(shù)值迭代方法中最快的。但是，牛頓法需要求導(dǎo) F′(qt＋1)，如式（25），（26），顯然較為復(fù)雜。實(shí)際上，式（25）對于水庫地形復(fù)雜，庫容難以用多項(xiàng)式函數(shù)表示的情況難以適應(yīng)，式（26）對于機(jī)組效率k0非恒定（與發(fā)電流量有關(guān)）、水頭損失隨出力不均勻的影響等問題難以用函數(shù)形式直接求導(dǎo)，故牛頓迭代法實(shí)際應(yīng)用較為困難。

3 不動點(diǎn)迭代法與二分法比較

改進(jìn)歐拉迭代法已在文獻(xiàn)[3]中與二分法迭代相比較，精度要低于二分法計(jì)算成果，而且發(fā)電流量仍需迭代，效率總體提高不大。牛頓迭代法由于要求導(dǎo)計(jì)算較為復(fù)雜，且對出力不均勻的水頭損失、出力系數(shù)隨裝機(jī)效率變化的情況難以用函數(shù)簡單考慮，故在實(shí)用中使用較為困難。文中以某水電站44年長系列逐月調(diào)節(jié)為例，選擇較為實(shí)用的改進(jìn)的不動點(diǎn)迭代法和二分法迭代計(jì)算進(jìn)行比較 ?？梢钥闯霾粍狱c(diǎn)迭代收斂速度約是二分法的3倍。

在大中型水庫、電站中，常需要推求多條調(diào)度線，以充分提高水庫電站的調(diào)節(jié)能力；在新建工程設(shè)計(jì)方案比較中，有些電站需要逐日調(diào)節(jié)以比較其效益差別；甚至有的水電站需要根據(jù)來水預(yù)報(bào)進(jìn)行實(shí)時(shí)調(diào)度，采用實(shí)用的不動點(diǎn)迭代法，可有效提高計(jì)算效率。

[1]葉秉如.水利計(jì)算及水資源規(guī)劃[M].北京：中國水利水電出版社,2003.

[2]李慶揚(yáng)，王能超，易大義.數(shù)值分析[M].北京：清華大學(xué)出版社,2008.

[3]陳守煜.水電站水庫徑流調(diào)節(jié)的迭代數(shù)值計(jì)算[J].水利學(xué)報(bào)，1981，（2）：1-10.

[4]張卷美.一類不動點(diǎn)迭代法的求解[J].河南理工大學(xué)學(xué)報(bào)，2006，（2）：169-171.