水電站徑流調節的數值解法

張曉波，鄭雄偉，陳昌軍

（浙江省水利水電勘測設計院，浙江 杭州 310002）

隨著計算機技術的廣泛引用，水電站徑流調節計算方法已由傳統的列表試算法、圖解法轉向數值計算方法。大連理工大學于1981年首次提出采用改進歐拉法進行數值求解，受此啟發，筆者將水電站徑流調節計算的數值解法進行拓展，分別采用二分法、改進歐拉法、不動點迭代法、牛頓迭代法對水電站徑流調節進行數值求解推導，并從迭代效率、計算精度、實用性等方面進行比較，推薦實用、簡便的不動點迭代法應用于徑流調節計算。

1 計算原理

水電站徑流調節基于以下基本方程[1]：

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

需要指出，地形復雜的水庫其水位庫容曲線很難用式（3）完整地表達,很多時候需采用分段函數擬合，程序計算應采用插值法。

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

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

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

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

2 數值解法

2.1 改進歐拉法

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

歐拉法相當粗略，多次推進后，歐拉法計算的Pn顯著偏離實際曲線，因此提出了改進歐拉法[2]：

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

圖1 歐拉迭代法

式（2）可變形為：

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

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

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

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

2.2 二分法

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

圖2 二分迭代法

（2）應用求解。將式（1）離散化，得：

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

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

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

式（17）為僅包含未知量 qt＋1的非線形方程。因此，水電站徑流調節計算可轉化為非線性方程的求解問題。

為求解式（17）：F(qt＋1)=0。可以令 a0=0，b0=qm，qm為最大過機流量。

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

2.3 不動點迭代法

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

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

圖3 不動點迭代法

（2）應用求解。將式（15）等價改寫為：

由此可以建立迭代關系式（20），（21）：

（上標 k 為迭代次數，k=0，1，2，…）

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

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

于是對水電站徑流調節迭代公式有：

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

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

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

（2）應用求解。對于式（15）建立迭代關系為：

聯立式（15），（24），（25），（26）可求解。

圖4 牛頓迭代法

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

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

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

3 不動點迭代法與二分法比較

改進歐拉迭代法已在文獻[3]中與二分法迭代相比較，精度要低于二分法計算成果，而且發電流量仍需迭代，效率總體提高不大。牛頓迭代法由于要求導計算較為復雜，且對出力不均勻的水頭損失、出力系數隨裝機效率變化的情況難以用函數簡單考慮，故在實用中使用較為困難。文中以某水電站44年長系列逐月調節為例，選擇較為實用的改進的不動點迭代法和二分法迭代計算進行比較 。可以看出不動點迭代收斂速度約是二分法的3倍。

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

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

[2]李慶揚，王能超，易大義.數值分析[M].北京：清華大學出版社,2008.

[3]陳守煜.水電站水庫徑流調節的迭代數值計算[J].水利學報，1981，（2）：1-10.

[4]張卷美.一類不動點迭代法的求解[J].河南理工大學學報，2006，（2）：169-171.