綜合利用水庫的優(yōu)化調(diào)度受多因素影響，如徑流，水庫特性、用水特性以及電站的機電特性等，其中徑流的影響較大。本文采用馬爾可夫單鏈彈性相關(guān)理論處理徑流，以供水流量為決策變量，在考慮有效雨量的基礎(chǔ)上建立了動態(tài)規(guī)劃數(shù)學(xué)模型，編制了結(jié)構(gòu)簡明，功能完善，便于操作使用的大型優(yōu)化調(diào)度計算程序，自動繪制出三維優(yōu)化調(diào)度圖，利用優(yōu)化調(diào)度圖進行綜合利用水庫調(diào)節(jié)計算，在幾乎不增加投資的條件下，產(chǎn)生了巨大的經(jīng)濟效益。經(jīng)實踐證明，本方法準(zhǔn)確可靠，適合于大、中、小型水庫，也適合于平原水庫、地下水庫；更適合于我國北方水資源緊缺地區(qū)使用。

1 采用離散的馬爾可夫隨機過程描述徑流

1.1 用馬爾可夫過程描述徑流

　　為了計算和應(yīng)用的方便，將時間序列離散化(即分為若干時段：月)，相鄰時段存在著依賴關(guān)系，以水庫來水的3個相鄰時段t₁、t₂、t₃間徑流關(guān)系進行分析。用X₁、X₂、X₃表示3個時段的徑流，三者之間的相關(guān)情況可分為2種情況：(1)直接相關(guān)。即不管X₂取值怎樣(或不計X₂取值的影響)的條件下，X₁與X₃相關(guān)，稱為偏相關(guān)，其相關(guān)程度用相關(guān)系數(shù)表征，可用數(shù)量表示為γ₁₃。(2)間接相關(guān)。即因存在著X1和X2、X2和X3之間的相鄰時段相關(guān)關(guān)系，故X₁的大小影響著X₂的大小，從而又影響著X₃的大小。這種相關(guān)是由中間量X₂傳遞的，不是直接的，因此叫間接相關(guān)。

1.2 計算相應(yīng)條件概率

　　當(dāng)一年分成K個時段(月)，每個時段的徑流以平均值來表示，記作Q_K(K=1,2,3,……,K)。 

　　應(yīng)用相關(guān)理論分析，可以確定相鄰時段徑流Q_K,Q_K-1(如圖1所示)的條件概率分布函Q_K,Q_K-1的條件概率分布函數(shù)示意 數(shù)F(Q_K/Q_K-1)。其條件概率分布是一個二維分布，用概率理論及水文統(tǒng)計原理來推 求徑流的條件概率計算式。

研究相鄰時段的徑流相關(guān)聯(lián)系時，應(yīng)用相關(guān)系數(shù)R及回歸方程式求得

隔時段相關(guān)系數(shù)則為：

　　式中：Q_1i,Q_2i,Q_3i為第i年相鄰時段的實測徑流值；為平均值；n為徑流實測系列年數(shù)。本時段徑流的相關(guān)關(guān)系，應(yīng)用相關(guān)中的直線相關(guān)，以自回歸線性公式來表示：

　　式中：σ_K,σ_K-1分別為時段t_k,t_k-1的徑流均方差；R₁為相鄰時段徑流之間的相關(guān)系數(shù)。 

　　相鄰時段徑流之間應(yīng)用自回歸線性相關(guān)時，其間隔時段的徑流對回歸線的偏離值即誤差的分布，經(jīng)剛性和彈性相關(guān)比較后，采用了彈性相關(guān)處理方法即偏態(tài)分布，按皮爾遜Ⅲ型曲線分布。相應(yīng)于條件概率的流量Q_PK��可由下式求得：

　　式中：條件變差系數(shù)，其中C_vk為變差系數(shù)。一年劃分為K個時段，每個時段的徑流劃分為M級(即M個狀態(tài))，則相鄰時段的轉(zhuǎn)移概率：P_kij(k=1,2,3,……,k;i,j=1,2,3,……,M)表示的含義是t_k-1時段徑流為狀態(tài)i時，t_k時段徑流為狀態(tài)j時的概率

而矩陣

則表示t_k-1時段到tk時段狀態(tài)的轉(zhuǎn)移概率矩陣，顯然，這個矩陣的每行各非負(fù)元素之和為1，即：

　　為了計算P_kij轉(zhuǎn)移概率的方便，取等分的10個概率5%，15%，……95%，這樣轉(zhuǎn)移概率的值都為0.1，則相應(yīng)的條件概率的流量Q_pi由式(4)即可求得。

2 動態(tài)規(guī)劃

　　動態(tài)規(guī)劃法是美國數(shù)學(xué)家貝爾曼提出的，是一種研究多階段決策過程的數(shù)學(xué)方法。近年來廣泛應(yīng)用于水資源規(guī)劃管理領(lǐng)域中

　　2.1 動態(tài)規(guī)劃數(shù)學(xué)模型

　　把徑流當(dāng)作隨機過程的水庫優(yōu)化調(diào)度圖的計算是一個多階段的隨機決策過程。它的計算模型如下。 

　　(1)階段：將水庫調(diào)度圖按月(或者旬)劃分成12個相互關(guān)連的階段(時段)，以便求解

　　(2)狀態(tài)：因相鄰兩個階段的入庫平均流量Q_t和Q_t+1之間有相關(guān)關(guān)系，以面臨時段初的庫水位和本時段預(yù)報徑流量Q_t為狀態(tài)變量S_t(Z_t-1，Q_t)

　　(3)決策：在時段狀態(tài)確定后，作一個相應(yīng)的決定，即面臨時段的供水量q_t，同時確定了時段末水位，進行狀態(tài)轉(zhuǎn)移。水庫水位分M級，故有M個狀態(tài)轉(zhuǎn)移，按0.618法在決策域內(nèi)優(yōu)選，對每一個狀態(tài)變量S_t要選擇一最優(yōu)供水量q_t，S_t～q_t關(guān)系曲線為時段t的調(diào)度線，決策域為(Q_Dmin,t;Q_xmax,t)

　　對決策變量供水量q_t進行所有狀態(tài)優(yōu)選計算時，還要進行庫水位限制的檢查判別，若時段末蓄水量V₂大于允許的最高蓄水位或限制水位，則在水庫蓄滿前供水量仍按q_t放水計算，當(dāng)水庫蓄滿后則按入庫水量供水。當(dāng)入庫水量大于電廠最大過水能力時，超過部分作為棄水

　　(4)狀態(tài)轉(zhuǎn)移：水庫狀態(tài)和調(diào)度圖形式有關(guān)，因考慮當(dāng)時入庫徑流和短期徑流因素，水庫調(diào)度中將一年劃分為K個時段，每個時段由時段初庫水位Z��初和時段流量Q_t組成水庫的運行狀態(tài)，而每一種狀態(tài)有一個相應(yīng)的決策變量供水流量q_t，用函數(shù)關(guān)系表示為：

t_k為時段數(shù)，每一個決策就有一個相應(yīng)的時段末庫水位，水庫進行了狀態(tài)轉(zhuǎn)移，若將水庫的水位劃分為Z級，徑流劃分為M級。一個時段的水庫面臨狀態(tài)有Z×M種，全年水庫運行狀態(tài)有K×Z×M種，水庫優(yōu)化調(diào)度圖就是對全年各種運行狀態(tài)作出相應(yīng)決策變量的關(guān)系圖。 

　　由式(7)可知，當(dāng)時段tk的初始庫水位和徑流量已定時，時段的最優(yōu)決策供水量是一個定值，因而下一時段tK+1的初始庫水位(即時段t_k末的水位)也就是一個確定值。由于下一時段t_K+1的徑流不是一個確定值，而是依時段tK的徑流Q_t變化的隨機值，其值由條件概率分布函數(shù)(彈性相關(guān))決策。因此，水庫在時段tK處于狀態(tài)i，而時　　段t_K+1處于狀態(tài)j的狀態(tài)轉(zhuǎn)移概率為P_kij，則有，而矩陣P_k=(P_kij)則表示從時段t_K到時段t_K+1的水庫狀態(tài)轉(zhuǎn)移概率矩陣，P_k完全由時段t_K的調(diào)度方式和徑流狀態(tài)轉(zhuǎn)移矩陣決定。經(jīng)過多年運行后，水庫的運行狀態(tài)達(dá)到一個穩(wěn)定的概率分布

　　(5)效益函數(shù)：水庫進行狀態(tài)轉(zhuǎn)移，伴隨著產(chǎn)生了效益函數(shù)(包括了工業(yè)用水、生活用水、灌溉用水、發(fā)電用水及三個保證率)

　　其中灌溉用水：因灌溉需水量每年、每月、每天都不相同，因此是隨機變量，極難編制計算機程序計算，故首次引入《農(nóng)田水利學(xué)》的“有效雨量”概念，使整個優(yōu)化計算大大簡化，完全解決了水量平衡問題，整個優(yōu)化計算，水量平衡達(dá)到100%

　　有效雨量的計算：從水庫灌區(qū)試驗站獲取資料Mij即從1952～1999年歷年(i=1952～1999，j為第i年各月(或旬))的灌溉定額(是由歷年灌溉試驗站實測作物需水量采用通用電算程序計算出的)，而M_max是48年中最枯水年的灌溉定額。M_max-M_ij=P_0ij,i=1952,…,1999,j=1,…,12,逐一列表進行計算。把每年每月的有效雨量加到每年每月的來水量Q_t中，因Mmax是常數(shù)，所以僅有隨機變量Mij。其數(shù)學(xué)表達(dá)式如下：C_ixj=A_ixj-B_ixj，即：

　　式中C_ij為i年系列j時段(月)的有效雨量，a^ij為i年系列j時段農(nóng)作物需水量(j可按日計算后歸納成各農(nóng)作物生長期所需水量，再換算成月)。b_ij為i年系列j時段各類農(nóng)作物綜合灌溉水量。 

　　(6)目標(biāo)函數(shù)：根據(jù)水庫水資源不足的具體情況，擬定在滿足生活用水和工業(yè)用水保證率的條件下，盡量滿足農(nóng)業(yè)用水。目標(biāo)函數(shù)可表示為：滿足用水量保證率條件下供水量最大。目標(biāo)函數(shù)計算可用下列分段線性函數(shù)求得：

　　式中：q_t為水庫供水量，Q_Dmin為系統(tǒng)供水下限，即保證城市生活用水和工業(yè)用水的下限；Q_xmin為農(nóng)業(yè)保證供水量與Q_Dmin之和；Q_Dmax為電廠的最大過水能力；Qxmax為農(nóng)業(yè)供水量上限與Q_Dmin之和；C_E為發(fā)電專用水量小于Q_xmin時的折算系數(shù)，C_A為供水量小于Q_xmin時的折算系數(shù)，在計算中，可先任意假設(shè)C_A、C_E,C_A、C_E與Q_xmin的保證率成正比。給定一個C_A、C_E就可遞推得出一張優(yōu)化調(diào)度圖，用水庫多年入庫流量資料按調(diào)度圖進行歷時操作計算，若計算結(jié)果所得保證率低于要求的保證率，則修改C_A、C_E重新遞推計算(一般遞推2～3次即可)，求得另一優(yōu)化調(diào)度圖，再進行歷時操作，直至所得保證率符合要求為止。即經(jīng)過試算選擇滿足保證率要求的C_A、C_E值。 

　　2.2 動態(tài)規(guī)劃遞推方程 以q_t為t階段的決策變量，S_t(Z_t-1,Q_t)為t階段的狀態(tài)變量，則其逆時序動態(tài)規(guī)劃最優(yōu)遞推方程為：

　　式中：F_t(S_t,q_t)代表水庫從時刻t處于狀態(tài)St出發(fā)至水庫運行終了時刻N(計算周期末)的目標(biāo)函數(shù)值；ft(St,qt)代表時刻t水庫處于狀態(tài)S_t取供水量q_t時面臨時段效益期望值；F_t+1(S_t+1)代表水庫從時刻t+1處于S_t+1(j狀態(tài))出發(fā)至?xí)r刻期間各時段均采用最優(yōu)決策時所得的效益期望值；Qt表示計算中t時段所用的入庫徑流序列；p_i,j為t時刻采取q_t決策，系統(tǒng)由第t階段的第i種狀態(tài)S_t轉(zhuǎn)移為第t+1階段的第j種狀態(tài)S_t+1時的條件概率，F(xiàn)_t+1相應(yīng)S_t+1狀態(tài)最優(yōu)決策的效益。 

　　遞推方程的約束條件如下：①庫水位約束V_min,t≤V_t≤V_max,t，即各時段的庫水位不低于死水位V_min,t，也不能超過該時段允許的最高蓄水位V_max,t。②水量平衡約束V_t+1=V_t+(Q_t-q_t)·Δ_t-y_t-E_t,式中V_t+1、V_t代表時段t末、初的蓄水量；Q_t、q_t代表t時段平均入庫徑流量和供水量；yt為棄水量，Et為水庫蒸發(fā)滲漏損失。③供水約束和輸水能力約束Q_Dmax,t≥q_t≥Q_Dmin,t。t時段內(nèi)供水量不能超過水輪機的最大過水能力Q_Dmax,t，也不能小于下限Q_Dmin,t

　　2.3 動態(tài)規(guī)劃遞推計算 采取逆時序逐時段動態(tài)規(guī)劃遞推計算，即每時段對所有狀態(tài)逐一地優(yōu)選對應(yīng)的最優(yōu)決策。對時段的多個入庫流量代表值所產(chǎn)生的效益期望值。優(yōu)選方法采用0.618法，規(guī)定搜索點為20個

　　2.4 優(yōu)化調(diào)度圖 Howard��用Z變換方法證明式(10)隨年數(shù)t增加計算是收斂的，進行遞推計算采取逆時序遞推，即從N時段開始遞推到1時段，只要知道F_N(S_N)即可按式(10)遞推計算。開始可取庫水位(庫容)～蓄水量關(guān)系曲線作為初始遞推線F_N(S_N)。當(dāng)對第一個時段的所有狀態(tài)優(yōu)選出最優(yōu)決策后，即可往前遞推一個時段。當(dāng)?shù)谝荒曛饌�時段全部遞推計算完畢后，還要進行第二年周期的遞推計算，是因為初始遞推F_N(S_N)是任意假設(shè)的，故第一年周期遞推所得的策略并非穩(wěn)定的最優(yōu)策略，必需繼續(xù)遞推至各時段的遞推線均收斂為止，這時所得的策略才是穩(wěn)定的最優(yōu)策略。遞推線收斂的準(zhǔn)則是：前后兩年周期中同一時段的遞推線相差小于規(guī)定的相對誤差ε即：

　　式中：F_t(S_i)ⁿ代表第n年時段t遞推線上相應(yīng)于狀態(tài)Si的未來效益值；F_t(S_i)⁽ⁿ⁺¹⁾則是第n+1年時段t遞推線上同一狀態(tài)Si相應(yīng)的未來效益值，ε取0.001。一般最多遞推兩年就可以收斂，即可得出12時段或36個時段(旬)的最優(yōu)調(diào)度線。這時各時段的最優(yōu)決策構(gòu)成一個最優(yōu)策略，即為優(yōu)化調(diào)度圖。顯然，因考慮月(或旬)、相隔月(旬)的相關(guān)，即多用了一項概率預(yù)報，則相應(yīng)增加了經(jīng)濟效益。由于采用了馬爾可夫單鏈彈性相關(guān)理論對徑流進行處理，使水庫調(diào)度圖從二維坐標(biāo)變成三維坐標(biāo)，形成空間水庫優(yōu)化調(diào)度圖，再由調(diào)度圖換成一組以Qt為參數(shù)的方程，遞推線也由一條變成一組，即優(yōu)化調(diào)度線由一條線變成一組，形成一族調(diào)度曲線圖，為便于實際調(diào)度時使用。 

　　2.5 動態(tài)規(guī)劃計算程序 動態(tài)規(guī)劃的計算是一個非常復(fù)雜的過程，不同的規(guī)劃問題，要用不同的計算程序。我們根據(jù)最優(yōu)化(opt)問題的數(shù)學(xué)模型^［2］，用VISUL C編制了計算程序，用遞推方程找出最優(yōu)解。該程序在PⅡ微機上調(diào)試成功，經(jīng)實踐證明其具有功能強大，使用方便，運行速度快等特點，并能自動繪出三維空間水庫優(yōu)化調(diào)度圖及帶有一組參數(shù)的調(diào)度曲線圖。

3 應(yīng)用示例

　　本方法已應(yīng)用于山東沐浴、跋山和黃前等幾個大中型水庫，都取得理想效果。僅以沐浴水庫多目標(biāo)優(yōu)化調(diào)度的應(yīng)用情況來說明。

　　應(yīng)用滿足用水保證率條件下供水量最大為目標(biāo)函數(shù)合理地解決3個保證率的計算問題；建立了動態(tài)規(guī)劃數(shù)學(xué)模型^［5］，利用其優(yōu)化調(diào)度程序計算，計算結(jié)果理想，輸出了大量的表格，(如表1所示，限于篇幅，僅列一小部分)，并自動繪出了水庫優(yōu)化調(diào)度空間圖及多族調(diào)度曲線圖(如圖3、4所示)。利用優(yōu)化調(diào)度圖進行綜合調(diào)節(jié)計算，在幾乎不增加投資的情況下，增加了巨大的經(jīng)濟效益。

    依據(jù)制定的水庫優(yōu)化調(diào)度圖即馬爾可夫調(diào)度圖，對1952～1999年共48年水文年度的徑流資料進行長系列操作計算，計算結(jié)果表明，綜合利用水庫優(yōu)化調(diào)度后，工業(yè)用水保證率為95%，生活用水保證率為97%，灌溉保證率為80.5%；多年平均年發(fā)電量為384.7萬kW·h。灌溉保證率較常規(guī)調(diào)節(jié)計算的保證率75%增加到80.5%。如維持常規(guī)計算的灌溉保證率75%，則灌溉面積可從0.97萬hm²擴灌到1萬hm²。原沐浴水電站設(shè)計書的多年平均年發(fā)電量為311.3萬kW·h，優(yōu)化調(diào)度后年發(fā)電量凈增73萬kW·h，增加發(fā)電量24%。常規(guī)水量平衡48年總棄水量為40102.27萬m³，優(yōu)化調(diào)度后棄水量大大減少，僅棄水2335.14萬m³。

4  結(jié)語