考慮水庫(kù)來(lái)用水過(guò)程關(guān)聯(lián)性的多維隨機(jī)動(dòng)態(tài)規(guī)劃算法

郭愛(ài)軍，暢建霞※，王義民，黃 強(qiáng)，吳 彬，張 春

（1. 西安理工大學(xué)西北旱區(qū)生態(tài)水利國(guó)家重點(diǎn)實(shí)驗(yàn)室，西安 710048；2. 華能西藏雅魯藏布江水電開(kāi)發(fā)投資有限公司，拉薩 850000）

0 引 言

水庫(kù)是實(shí)現(xiàn)水資源時(shí)空均衡分配、解決水資源供需矛盾的關(guān)鍵工程措施。水庫(kù)興利調(diào)度作為水庫(kù)控制運(yùn)用的核心非工程措施，旨在根據(jù)水庫(kù)開(kāi)發(fā)任務(wù)和綜合利用要求，根據(jù)來(lái)、用水過(guò)程制定水庫(kù)蓄泄規(guī)則，最大限度滿足各個(gè)用水部門(mén)用水需求，實(shí)現(xiàn)水資源供給過(guò)程與需求過(guò)程的高度匹配。根據(jù)對(duì)水庫(kù)入庫(kù)徑流的處理方式，水庫(kù)調(diào)度可分為確定性調(diào)度與隨機(jī)調(diào)度，確定性調(diào)度認(rèn)為入庫(kù)徑流是已知和確定的，但往往受限于歷史有限的實(shí)測(cè)資料，難以全面反映徑流變化過(guò)程；隨機(jī)調(diào)度將徑流作為具有一定統(tǒng)計(jì)特性的隨機(jī)變量或隨機(jī)過(guò)程處理，能夠更好地體現(xiàn)徑流序列的不確定性。

作為求解水庫(kù)中長(zhǎng)期隨機(jī)優(yōu)化調(diào)度模型的算法之一，隨機(jī)動(dòng)態(tài)規(guī)劃算法自提出之后在流域水資源管理和水庫(kù)調(diào)度中被廣泛采用。黃強(qiáng)分別采用確定性、獨(dú)立隨機(jī)性、簡(jiǎn)單馬氏性過(guò)程與混合性過(guò)程描述了徑流過(guò)程，并改進(jìn)了有徑流預(yù)報(bào)下的隨機(jī)動(dòng)態(tài)規(guī)劃求解策略，應(yīng)用于以發(fā)電量最大為目標(biāo)的水庫(kù)優(yōu)化調(diào)度模型中；周惠成等以二灘水電站為例，提出了有、無(wú)預(yù)報(bào)與隨機(jī)過(guò)程相結(jié)合的徑流描述方法，建立了發(fā)電優(yōu)化調(diào)度模型，制定了電站優(yōu)化調(diào)度圖；盧迪等采用徑流預(yù)報(bào)概率修正先驗(yàn)概率來(lái)描述徑流的不確定性，建立了貝葉斯隨機(jī)動(dòng)態(tài)規(guī)劃模型，分析了有跨流域引水的碧流河水庫(kù)調(diào)度策略；崔遠(yuǎn)來(lái)等考慮降雨過(guò)程的隨機(jī)性，建立了解稻田高效節(jié)水優(yōu)化灌溉制度的隨機(jī)動(dòng)態(tài)規(guī)劃模型；Nop等基于降雨過(guò)程的馬爾可夫性，將隨機(jī)動(dòng)態(tài)規(guī)劃算法應(yīng)用于雨水收集系統(tǒng)中的水資源管理問(wèn)題；Cassagnole等分析了徑流預(yù)報(bào)質(zhì)量的不同對(duì)水庫(kù)電站效益的影響；Zhou等采用離散時(shí)間馬爾可夫鏈描述水庫(kù)入庫(kù)徑流過(guò)程，采用隨機(jī)動(dòng)態(tài)規(guī)劃算法以及機(jī)會(huì)約束模型推求了溪洛渡水電站的中期發(fā)電運(yùn)用計(jì)劃。

如上所述，隨機(jī)動(dòng)態(tài)規(guī)劃算法的核心之一是入庫(kù)徑流或降雨等水文過(guò)程的隨機(jī)性的描述。農(nóng)業(yè)灌溉用水是中國(guó)的主要用水部門(mén)，其需水量和需水過(guò)程受氣候條件的影響顯著，具有很強(qiáng)的隨機(jī)性；并且，農(nóng)業(yè)灌溉需水與入庫(kù)天然徑流同受降雨、氣溫等氣象因素影響，在流域氣象條件的相似性背景下，農(nóng)業(yè)灌溉需水與入庫(kù)天然徑流呈一定的負(fù)相關(guān)性，且流域范圍越小，負(fù)相關(guān)性越為顯著。如何在水庫(kù)隨機(jī)優(yōu)化調(diào)度中同時(shí)考慮需水與徑流二維過(guò)程的隨機(jī)性與關(guān)聯(lián)性，以往研究對(duì)此考慮不足。有限的研究也主要是從徑流的隨機(jī)性角度出發(fā)，忽略了農(nóng)業(yè)灌溉需水的隨機(jī)性及其與徑流之間的關(guān)聯(lián)性，使得隨機(jī)規(guī)劃理論模型與現(xiàn)實(shí)水文過(guò)程之間的匹配性不足。本研究結(jié)合馬爾可夫過(guò)程與Copula函數(shù)描述需水與徑流二維過(guò)程的隨機(jī)性與關(guān)聯(lián)性，基于隨機(jī)動(dòng)態(tài)規(guī)劃理論，提出考慮水庫(kù)來(lái)、需水過(guò)程關(guān)聯(lián)性與隨機(jī)性的水庫(kù)優(yōu)化調(diào)度算法，并進(jìn)行相關(guān)案例分析。該算法對(duì)于農(nóng)業(yè)水土資源隨機(jī)優(yōu)化配置等領(lǐng)域的優(yōu)化決策求解具有支撐作用。

1 水庫(kù)來(lái)、用水多維過(guò)程關(guān)聯(lián)性與隨機(jī)性

1.1 一維過(guò)程隨機(jī)性與關(guān)聯(lián)性

隨機(jī)動(dòng)態(tài)規(guī)劃是確定性動(dòng)態(tài)規(guī)劃與隨機(jī)過(guò)程二者的結(jié)合。相較于確定性動(dòng)態(tài)規(guī)劃，隨機(jī)動(dòng)態(tài)規(guī)劃的復(fù)雜性主要體現(xiàn)在隨機(jī)過(guò)程的處理方面，即固定階段（或時(shí)段）隨機(jī)變量的隨機(jī)性描述與階段（或時(shí)段）之間隨機(jī)變量的關(guān)聯(lián)性；求解步驟與思路方面，隨機(jī)動(dòng)態(tài)規(guī)劃與確定性動(dòng)態(tài)規(guī)劃一致，隨機(jī)動(dòng)態(tài)規(guī)劃具體求解步驟可參考文獻(xiàn)[8, 22]。

徑流或需水一維隨機(jī)過(guò)程的描述可從固定時(shí)段下變量的隨機(jī)性與相鄰時(shí)段間變量的關(guān)聯(lián)性兩個(gè)方面展開(kāi)。此處，固定時(shí)段指以水庫(kù)調(diào)度時(shí)間尺度為標(biāo)準(zhǔn)，如開(kāi)展以月為時(shí)間尺度的中長(zhǎng)期水庫(kù)調(diào)度時(shí)，固定時(shí)段則為一月、二月、三月、……、十二月，重點(diǎn)關(guān)注各個(gè)固定時(shí)段下變量的隨機(jī)性，如長(zhǎng)系列一月徑流的隨機(jī)性等；相鄰時(shí)段則以固定時(shí)段相鄰為具體情形，如一月與二月相鄰時(shí)段、二月與三月相鄰時(shí)段等，本文重點(diǎn)關(guān)注相鄰時(shí)段間隨機(jī)變量關(guān)聯(lián)性，采用狀態(tài)轉(zhuǎn)移概率矩陣來(lái)表征兩兩隨機(jī)變量不同狀態(tài)間的關(guān)聯(lián)特性。

1）固定時(shí)段下隨機(jī)性描述

徑流與需水的隨機(jī)性采用傳統(tǒng)單變量概率分布函數(shù)描述。本次研究選擇水文領(lǐng)域應(yīng)用較廣的伽馬分布（Gam）、對(duì)數(shù)正態(tài)分布（Logn）、廣義極值分布（Gev）、威布爾分布（Weibull）、正態(tài)分布（Norm）、指數(shù)分布（Exp）作為擬合固定時(shí)段下徑流或需水隨機(jī)變量的備選概率分布類(lèi)型，采用極大似然法（Maximum Likelihood Estimate，MLE）估計(jì)不同概率分布函數(shù)參數(shù)，選用王義民等提出的基于多評(píng)價(jià)準(zhǔn)則的綜合指數(shù)法選擇適宜概率分布函數(shù)?；诙嘣u(píng)價(jià)準(zhǔn)則的綜合指數(shù)包括修正的赤池信息準(zhǔn)則（Corrected Akaike’s Information Criterion）、離差絕對(duì)值和最小準(zhǔn)則、均方根誤差和相對(duì)離差平方和最小準(zhǔn)則，利用不同準(zhǔn)則所側(cè)重的評(píng)價(jià)方向，綜合優(yōu)選適宜概率分布函數(shù)。

2）相鄰時(shí)段間過(guò)程關(guān)聯(lián)性描述

隨機(jī)變量時(shí)段之間的關(guān)聯(lián)性影響隨機(jī)過(guò)程描述的方式。若相鄰時(shí)段之間隨機(jī)變量相互獨(dú)立，則可判斷該隨機(jī)過(guò)程為獨(dú)立過(guò)程，隨機(jī)變量不同狀態(tài)間的轉(zhuǎn)移概率為0；若相鄰時(shí)段之間的隨機(jī)變量相關(guān)性較強(qiáng)，難以忽視，即對(duì)該隨機(jī)過(guò)程的演變而言，在已知目前狀態(tài)的條件下與以往的狀態(tài)無(wú)關(guān)，該特性稱為馬爾可夫性，若隨機(jī)過(guò)程轉(zhuǎn)移至下一階段某一狀態(tài)時(shí)只與系統(tǒng)當(dāng)前狀態(tài)相關(guān)，稱之為一階馬爾可夫過(guò)程或簡(jiǎn)單馬爾可夫過(guò)程，若與其前面連續(xù)的個(gè)狀態(tài)相關(guān)，與個(gè)狀態(tài)以前的狀態(tài)無(wú)關(guān)，稱之為階馬爾可夫過(guò)程或復(fù)雜馬爾科夫過(guò)程。徑流與需水過(guò)程是典型的隨機(jī)過(guò)程，若相鄰時(shí)段之間或當(dāng)前時(shí)段與前面連續(xù)的個(gè)時(shí)段的變量具有顯著相關(guān)性，則可認(rèn)為徑流與需水過(guò)程是馬爾可夫過(guò)程；由于多時(shí)段相關(guān)的徑流與需水過(guò)程描述較為復(fù)雜，實(shí)際應(yīng)用中也較少，因此本次研究重點(diǎn)關(guān)注相鄰時(shí)段間隨機(jī)變量的狀態(tài)轉(zhuǎn)移過(guò)程，即未來(lái)時(shí)段徑流與需水狀態(tài)的演變只與當(dāng)前時(shí)段狀態(tài)有關(guān)。具體做法為，根據(jù)時(shí)段之間變量的相關(guān)性是否顯著，采用簡(jiǎn)單馬爾科夫過(guò)程（時(shí)段間變量相關(guān)性顯著）或獨(dú)立過(guò)程（時(shí)段之間變量相關(guān)性不顯著）描述徑流、需水過(guò)程。徑流、需水過(guò)程描述的關(guān)鍵是推求相鄰時(shí)段隨機(jī)變量不同狀態(tài)的條件轉(zhuǎn)移概率P。以相鄰時(shí)段徑流不同狀態(tài)轉(zhuǎn)移概率{Q()|Q()}計(jì)算為例，其公式如下：

式中為固定時(shí)段隨機(jī)變量離散總數(shù)；、為時(shí)段隨機(jī)變量第個(gè)狀態(tài)與時(shí)段+1隨機(jī)變量第個(gè)狀態(tài)；H(·)為時(shí)段與時(shí)段+1變量聯(lián)合分布函數(shù)；F(·)為時(shí)段變量概率分布函數(shù)；q、q為時(shí)段徑流的最小值與最大值，m/s；、分別為時(shí)段徑流在狀態(tài)下的下限與上限；、分別為時(shí)段+1徑流在狀態(tài)下的下限與上限。

變量狀態(tài)數(shù)量為時(shí)，全部狀態(tài)轉(zhuǎn)移及相應(yīng)的轉(zhuǎn)移概率便有種，則轉(zhuǎn)移概率矩陣P為

此處，轉(zhuǎn)移概率矩陣滿足以下條件：

式中表示隨機(jī)變量的第個(gè)離散狀態(tài)。

本次研究選用在水文水資源領(lǐng)域、經(jīng)濟(jì)領(lǐng)域應(yīng)用廣泛的Copula函數(shù)表征時(shí)段與時(shí)段+1變量聯(lián)合分布函數(shù)H(·)。Copula函數(shù)是一種描述多變量相關(guān)性的連接函數(shù)，其優(yōu)勢(shì)在于可靈活選用邊緣分布函數(shù)且聯(lián)合分布函數(shù)構(gòu)造簡(jiǎn)單、靈活，易于求解。需要注意的是，Copula函數(shù)僅用于擬合變量間相關(guān)性顯著的時(shí)段，若變量間相關(guān)性不顯著，則認(rèn)為變量彼此獨(dú)立。

Copula函數(shù)(·)表達(dá)形式如下：

式中(·)為邊緣概率分布函數(shù)，()、()分別為變量、的邊緣概率分布函數(shù)。Copula函數(shù)發(fā)展至今，已有多種類(lèi)型（如橢圓型、Archimedean型、二次型等）。其中，Archimedean型Copula函數(shù)是水文水資源領(lǐng)域應(yīng)用較為廣泛的類(lèi)型之一。本次研究選擇Archimedean型Copula函數(shù)中的Clayton Copula、Gumbel-Hougaard Copula、Frank Copula三種函數(shù)為備選類(lèi)型，采用極大似然函數(shù)估計(jì)Copula函數(shù)的參數(shù)，選用Cramer-von Mises統(tǒng)計(jì)量進(jìn)行擬合優(yōu)度檢驗(yàn)，進(jìn)而選擇最優(yōu)的Copula函數(shù)類(lèi)型。

Copula函數(shù)的參數(shù)可用于量化變量之間關(guān)聯(lián)程度，以Clayton Copula函數(shù)為例，參數(shù)→0時(shí)，表明變量與趨于獨(dú)立，參數(shù)→+∞時(shí)，表明變量與趨于完全相關(guān)，參數(shù)在（0，+∞）內(nèi)變化，與肯德?tīng)栔认嚓P(guān)系數(shù)的關(guān)系為=/(+2)。

1.2 徑流與需水多維過(guò)程關(guān)聯(lián)性與隨機(jī)性

與一維隨機(jī)過(guò)程不同的，多維隨機(jī)過(guò)程考慮特定時(shí)段下多維變量之間的獨(dú)立關(guān)系或相關(guān)關(guān)系；相鄰時(shí)段間，考慮多維變量時(shí)段第個(gè)狀態(tài)內(nèi)部關(guān)聯(lián)性的基礎(chǔ)上向時(shí)段+1第個(gè)狀態(tài)轉(zhuǎn)移的過(guò)程，與1.1節(jié)單一變量的狀態(tài)轉(zhuǎn)移形成區(qū)別。

1）固定時(shí)段下徑流與需水多維變量隨機(jī)性描述

固定時(shí)段下多維變量隨機(jī)性的描述選用多維聯(lián)合分布函數(shù)。由于Copula函數(shù)不需提前假定變量之間的獨(dú)立或相關(guān)關(guān)系，通過(guò)該函數(shù)的參數(shù)值大小可反映多維變量之間的相依關(guān)系或獨(dú)立關(guān)系，故本研究選用Copula函數(shù)構(gòu)造多維變量的聯(lián)合分布函數(shù)，以此描述徑流（，m3/s）與需水（，m3/s）多維變量的隨機(jī)性。

式中(,)為徑流與灌溉需水的聯(lián)合概率，()、()為徑流與需水的邊緣分布函數(shù)，通過(guò)基于多評(píng)價(jià)準(zhǔn)則的綜合指數(shù)法從備選概率分布中優(yōu)選所得；q、w分別為長(zhǎng)系列情況下月的徑流與需水；(x< q≤x)表示q處于[x,x]區(qū)間的概率，此處[x,x]表示的是公式（2）對(duì)徑流離散的個(gè)狀態(tài)中的某個(gè)狀態(tài)的上下限；(x< w≤x)表示w處于[x,x]區(qū)間的概率，[x,x]表示的是對(duì)需水離散的個(gè)狀態(tài)中的某個(gè)狀態(tài)的上下限；θ、θ分別為邊際分布函數(shù)()、()的參數(shù)集。

若固定時(shí)段下徑流與需水不相關(guān)，則

2）相鄰時(shí)段間徑流與需水過(guò)程關(guān)聯(lián)性描述

以此次研究涉及的徑流與需水多維隨機(jī)過(guò)程為例，時(shí)段間水資源供需狀態(tài)轉(zhuǎn)移概率P{[(Q(),W()] |[(Q(),W()]}計(jì)算方法如下：

記時(shí)段徑流處于狀態(tài)（即公式（2）中對(duì)徑流離散的個(gè)狀態(tài)中的狀態(tài)，該狀態(tài)下徑流的上下限區(qū)間為[x,x]，公式可表示為Q() =x<q≤x）的概率為(Q())，其公式為

故，在時(shí)段徑流處于狀態(tài)條件下，時(shí)段+1徑流處于狀態(tài)（即公式（2）中對(duì)徑流離散的個(gè)狀態(tài)中的狀態(tài)，該狀態(tài)下徑流的上下限區(qū)間為[z,z]，公式可表示為Q() =z<q≤z）的概率為(Q()) =(Q()|Q())，其公式為

同理，時(shí)段需水處于狀態(tài)（即W() =xw≤x）的概率(W())為

對(duì)于上述線性思維模式影響下所造成的諸多歷史事實(shí)的長(zhǎng)期被“遮蔽”和歪曲，學(xué)界也不乏反思，如王瑞來(lái)《宋代士大夫主流精神論—以范仲淹為中心的考察》一文：

故，時(shí)段需水處于狀態(tài)條件下（即(W())），時(shí)段+1需水處于狀態(tài)（即W() =zw≤z）的概率為(W()) =(W()|W()) ，其公式為

考慮到變量與的關(guān)聯(lián)性，故時(shí)段+1多維變量處于狀態(tài)[Q(),W()]的聯(lián)合概率[Q(),W()]為

式中[Q(),W()]表示的是在時(shí)段多維變量處于狀態(tài)[(Q(),W()]條件下，時(shí)段+1多維變量處于狀態(tài)[Q(),W()]的條件概率P{[(Q(),W()] |[(Q(),W()]}），即

其中，式（9）與式（11）基于Copula函數(shù)的條件概率公式（式（1））計(jì)算；式（12）采用基于Copula函數(shù)的聯(lián)合概率計(jì)算。

若多維過(guò)程（,）彼此獨(dú)立，則

多維過(guò)程的一步狀態(tài)轉(zhuǎn)移矩陣為P一般式可寫(xiě)為

2 考慮多維過(guò)程關(guān)聯(lián)性與隨機(jī)性的隨機(jī)動(dòng)態(tài)規(guī)劃算法

式中f(V, Q())表示-1時(shí)段入庫(kù)徑流為Q()、水庫(kù)庫(kù)容為V情形下從時(shí)段到調(diào)度期末的最優(yōu)余留期效益；{Q()|Q()}表示時(shí)段-1時(shí)入庫(kù)徑流為Q()條件下時(shí)段入庫(kù)徑流為Q()的概率；O(V,Q())是當(dāng)前時(shí)段入庫(kù)徑流為Q()時(shí)的優(yōu)化調(diào)度目標(biāo)值；f(V, Q())表示時(shí)段入庫(kù)徑流為Q()、水庫(kù)庫(kù)容為V情形下從+1時(shí)段到調(diào)度期末的最優(yōu)余留期效益。

類(lèi)比可得多維隨機(jī)過(guò)程（以徑流與需水二維隨機(jī)過(guò)程為例）下，隨機(jī)動(dòng)態(tài)規(guī)劃遞推公式如下：

式中f(V,(Q(),W())) 表示-1時(shí)段入庫(kù)徑流為Q()、需水為W()、水庫(kù)庫(kù)容為V情形下從時(shí)段到調(diào)度期末的最優(yōu)余留期效益；{[Q(), W()]|[Q(),W()]}表示時(shí)段-1時(shí)入庫(kù)徑流與需水為[Q(),W()]條件下時(shí)段入庫(kù)徑流與需水為[Q(), W()]狀態(tài)的概率。

一維隨機(jī)過(guò)程與多維隨機(jī)過(guò)程時(shí)段間的狀態(tài)轉(zhuǎn)移示意圖如圖1所示。

圖1 一維與多維隨機(jī)過(guò)程下的狀態(tài)轉(zhuǎn)移過(guò)程 Fig.1 State transition process of one-dimensional and multi-dimensional stochastic process

3 實(shí)例研究

本文以涇河流域東莊水庫(kù)為例，水庫(kù)總庫(kù)容32.76億m，調(diào)節(jié)庫(kù)容5.78億m，為不完全年調(diào)節(jié)水庫(kù)，死水位與正常蓄水位分別為756、789 m。水庫(kù)供水對(duì)象主要為城鎮(zhèn)生活用水、工業(yè)用水與灌區(qū)灌溉用水。其中，城鎮(zhèn)生活與工業(yè)用水多年平均用水量2.13億m，年內(nèi)變化較為均勻，且壩前取水占65%，供水時(shí)優(yōu)先滿足二者用水要求；徑流序列采用1960—2013年還原后的天然入庫(kù)徑流序列；灌區(qū)灌溉用水量計(jì)算方面，首先采用國(guó)際糧農(nóng)組織（FAO）推薦的Penman-Monteith方法和作物系數(shù)法計(jì)算作物需水量和凈灌溉需水量，然后不同作物種植面積乘以對(duì)應(yīng)作物凈灌溉需水量并除以灌區(qū)灌溉水有效利用系數(shù)（以0.65計(jì)算）。研究以農(nóng)業(yè)灌溉用水缺水量最小為目標(biāo)，約束條件包括水量平衡約束、水庫(kù)水位約束、流量約束及變量非負(fù)約束。

3.1 數(shù)學(xué)模型

1）目標(biāo)函數(shù)

2）約束條件

① 水量平衡約束

式中V表示時(shí)段末水庫(kù)庫(kù)容，m；q表示時(shí)段水庫(kù)來(lái)水流量，m/s；q表示時(shí)段水庫(kù)出庫(kù)流量，m/s；q表示時(shí)段水庫(kù)上游取水流量，m/s；q表示時(shí)段水庫(kù)發(fā)電流量，m/s；q表示時(shí)段水庫(kù)棄水流量，m/s；q表示時(shí)段水庫(kù)下游河道生態(tài)流量，m/s；q表示時(shí)段水庫(kù)的灌區(qū)供水流量，m/s；q表示時(shí)段水庫(kù)其他供水流量（如工業(yè)及生活），m/s。

② 水位約束

式中Z為時(shí)段末水庫(kù)允許消落的最低水位，m；Z為時(shí)段末水庫(kù)允許達(dá)到的最高水位，m。該最高水位在汛期為防洪限制水位，在非汛期則為正常蓄水位。、Z分別為水庫(kù)整個(gè)調(diào)度期初、末水位，Z為水庫(kù)的起調(diào)水位（通常為死水位），m。

③ 流量約束

式中q、q為時(shí)段水庫(kù)允許的出庫(kù)流量最大值及最小值，m/s；q、q為時(shí)段水庫(kù)發(fā)電流量最大值及最小值，m/s；q為時(shí)段水庫(kù)下游河道所需最小生態(tài)流量，m/s。

④ 變量非負(fù)約束。調(diào)度模型中所有的變量均大于等于0。

3.2 徑流與需水過(guò)程隨機(jī)性與關(guān)聯(lián)性結(jié)果分析

采用Gam、Logn、Gev、Weibull、Norm、Exp擬合1960—2013年各月的徑流與需水?dāng)?shù)據(jù)，應(yīng)用基于多評(píng)價(jià)準(zhǔn)則的綜合指數(shù)法選擇適宜的概率分布函數(shù)，該方法的統(tǒng)計(jì)指標(biāo)值越小表明該概率分布函數(shù)擬合隨機(jī)變量的效果越優(yōu)。基于多評(píng)價(jià)準(zhǔn)則的綜合指數(shù)法的統(tǒng)計(jì)指標(biāo)結(jié)果見(jiàn)表1。由表可知，表征徑流隨機(jī)性的適宜概率分布函數(shù)不同月份之間差異較大，整體來(lái)看，Logn分布與Weibull分布函數(shù)選擇較多；需水方面，大多數(shù)月份下，Gev分布函數(shù)是描述需水隨機(jī)性的適宜函數(shù)。

表1 基于多評(píng)價(jià)準(zhǔn)則的綜合指數(shù)法統(tǒng)計(jì)指標(biāo) Table 1 Statistical index obtained by comprehensive index method based on multiple evaluation rules

采用肯德?tīng)栔认嚓P(guān)系數(shù)描述徑流、需水過(guò)程時(shí)段之間的關(guān)聯(lián)性，結(jié)果見(jiàn)圖2。由圖可知，徑流方面，除卻6 —7月、7—8月外，其余月份之間徑流過(guò)程呈顯著性相關(guān)；灌溉需水方面，1—2月、5—6月、7—8月以及12—1月四個(gè)相鄰時(shí)段灌溉需水相關(guān)性顯著，其他時(shí)段間，灌溉需水相關(guān)性不顯著。故，構(gòu)建多維隨機(jī)動(dòng)態(tài)規(guī)劃模型時(shí)，根據(jù)相鄰時(shí)段間徑流與需水的相關(guān)程度統(tǒng)計(jì)特性，需選擇性地采用適宜Copula函數(shù)表征一維變量時(shí)段之間關(guān)聯(lián)性。

圖2 徑流與需水隨機(jī)過(guò)程時(shí)段之間關(guān)聯(lián)性分析。 Fig.2 Relevancy analysis of runoff and water demand process at adjacent time

采用Cramer-von Mises方法選擇適宜的Copula函數(shù)，以此表征徑流、需水時(shí)段之間關(guān)聯(lián)性，結(jié)果表明：徑流方面，依照?qǐng)D2中所述時(shí)序，1—2月、2—3月、3—4月、4—5月、5—6月、8—9月、9—10月、10—11月、11 —12月、12—1月時(shí)段間徑流的適宜Copula函數(shù)依次為Frank、Frank、Gumbel-Hougaard、Gumbel-Hougaard、Frank、Clayton、 Frank、Gumbel-Hougaard、Gumbel-Hougaard與Frank；灌溉需水方面，1—2月、5 —6月、7—8月、12—1月時(shí)段間灌溉需水的適宜Copula函數(shù)依次為Gumbel-Hougaard、Gumbel-Hougaard、Clayton、Frank。采用如上Copula函數(shù)推求時(shí)段之間變量狀態(tài)轉(zhuǎn)移概率，節(jié)選5—6月為例，結(jié)果如圖3。以圖3a為例，5月徑流為7 m/s左右，轉(zhuǎn)移至6月時(shí)，徑流為6 m/s的概率最大，78 m/s時(shí)的可能性最??；同樣，5月徑流為200 m/s左右，轉(zhuǎn)移至6月時(shí)，徑流為78 m/s的概率最大，小于78 m/s情形下的概率逐漸減小。通過(guò)查閱該圖，實(shí)際生產(chǎn)管理人員可獲取不同月份之間徑流狀態(tài)轉(zhuǎn)移概率，推算在當(dāng)前月份的徑流狀態(tài)下，下月流域發(fā)生不同狀態(tài)徑流的可能性。

圖3 5—6月徑流與需水狀態(tài)轉(zhuǎn)移概率矩陣 Fig.3 State transitional probability matrix of runoff and water demand from May to June

采用肯德?tīng)栔认嚓P(guān)系數(shù)診斷二維隨機(jī)變量徑流與需水各個(gè)時(shí)段相關(guān)性，結(jié)果表明4月、5月、6月、7月、9月、10月徑流與需水呈顯著性相關(guān)，其余月份二者間相關(guān)關(guān)系不顯著。在此基礎(chǔ)上，選擇適宜Copula函數(shù)描述相關(guān)性顯著月份的徑流與需水聯(lián)合概率分布，相關(guān)性不顯著月份徑流與需水的聯(lián)合概率采用邊緣概率分布函數(shù)相乘方式獲得。結(jié)果表明，由于徑流與需水之間的負(fù)相關(guān)性，F(xiàn)rank Copula函數(shù)是擬合各月徑流與需水相依關(guān)系的最優(yōu)Copula函數(shù)。以5月份為例，采用Frank Copula函數(shù)表征徑流與需水二維變量的聯(lián)合分布，結(jié)果見(jiàn)圖4。以圖4中的累積概率分布函數(shù)（CDF，Cumulative Distribution Function）曲線為0.1為例，由于徑流與需水之間的負(fù)相關(guān)關(guān)系，需水大、徑流小的組合與需水小、徑流大的組合位于同一累積概率分布曲線上。

圖4 5月徑流與需水二維聯(lián)合概率分布 Fig.4 Two dimensional joint probability distribution of runoff and water demand in May

二維隨機(jī)變量時(shí)段之間關(guān)聯(lián)性通過(guò)二維組合過(guò)程的狀態(tài)轉(zhuǎn)移矩陣描述。圖5展示了采用公式（14）、公式（15）計(jì)算得到的在考慮關(guān)聯(lián)性（相關(guān)）與不考慮關(guān)聯(lián)性（獨(dú)立）兩種情況下，5月徑流與需水二維組合狀態(tài)[,] = [72.47, 24.70]、[,] = [17.23, 42.04]向6月份100種[,]組合狀態(tài)轉(zhuǎn)移的概率。此處，100即是公式（2）中對(duì)徑流與需水進(jìn)行離散后的狀態(tài)總數(shù)，圖5中的1至100的狀態(tài)數(shù)是分別將徑流與需水離散為10種狀態(tài)，然后兩者之間互相組合，共形成100種方案。其中，的10種變化方案為5.55、12.20、14.55、20.29、26.49、34.78、48.19、78.81、161.22、374.65 m/s，的10種變化方案為0.04、1.76、2.77、5.36、7.96、10.95、14.83、21.00、30.20、78.22 m/s，以為基準(zhǔn)，由大至小輪次組合，形成圖5中的6月份[,]100種方案。

圖5 考慮與不考慮徑流與需水相關(guān)性時(shí)狀態(tài)的轉(zhuǎn)移概率 Fig.5 State transition probability between runoff and water demand with and without considering relevancy

由圖可知，考慮徑流與需水關(guān)聯(lián)性與否兩種情形下，二維隨機(jī)過(guò)程狀態(tài)轉(zhuǎn)移概率不同，且不同轉(zhuǎn)移狀態(tài)下兩種情形計(jì)算所得的轉(zhuǎn)移概率大小不同，即考慮相關(guān)時(shí)的轉(zhuǎn)移概率并不是全大于或全小于獨(dú)立時(shí)的轉(zhuǎn)移概率；此外，圖5中兩種情形下轉(zhuǎn)移概率最大差別僅為0.04，差別相對(duì)較小，發(fā)生在5月的二維隨機(jī)過(guò)程[,] = [17.23, 42.04]向6月的二維隨機(jī)過(guò)程[12.20, 1.76]轉(zhuǎn)移過(guò)程中。結(jié)合隨機(jī)動(dòng)態(tài)規(guī)劃遞推公式（式（17））可知，考慮變量關(guān)聯(lián)性與否將影響最終最優(yōu)水位決策的制定，具體影響程度取決于不同情形下?tīng)顟B(tài)轉(zhuǎn)移概率的差異程度。

3.3 水庫(kù)優(yōu)化調(diào)度結(jié)果分析

為了對(duì)比分析考慮與不考慮徑流、需水多維變量關(guān)聯(lián)性對(duì)水庫(kù)優(yōu)化調(diào)度結(jié)果的影響，本次研究分別采用所提出的考慮徑流與需水關(guān)聯(lián)性、不考慮徑流與需水關(guān)聯(lián)性下的兩種算法求解優(yōu)化調(diào)度模型，推求了水庫(kù)優(yōu)化調(diào)度規(guī)則。兩種情形下水庫(kù)調(diào)度規(guī)則差別較小，圖6展示了考慮與不考慮徑流與需水關(guān)聯(lián)性時(shí)6月份的水庫(kù)水位最優(yōu)策略。[,]=[10.33, 13.93]情況下，考慮徑流與需水關(guān)聯(lián)性與不考慮二者關(guān)聯(lián)性水庫(kù)水位決策有一定的差異，不考慮二者關(guān)聯(lián)性（圖中帶正方形線條）時(shí)水庫(kù)蓄存水位比考慮二者關(guān)聯(lián)性（圖中帶六角形線條）時(shí)水庫(kù)蓄存水位較低；其余[,]情況下，考慮與不考慮二者關(guān)聯(lián)性時(shí)，水庫(kù)水位最優(yōu)決策過(guò)程無(wú)變化，故在圖6中重合。

圖6 考慮徑流與需水關(guān)聯(lián)性的六月份水庫(kù)最優(yōu)決策 Fig.6 Optimal reservoir operations in June considering the relevancy between runoff and water demand

將長(zhǎng)系列的來(lái)水與需水?dāng)?shù)據(jù)輸入水庫(kù)優(yōu)化調(diào)度規(guī)則，發(fā)現(xiàn)當(dāng)考慮變量關(guān)聯(lián)性時(shí)，長(zhǎng)系列下缺水總量為10.37 億m，不考慮時(shí)，缺水總量為10.49億m，考慮多維變量關(guān)聯(lián)性較不考慮時(shí)水庫(kù)調(diào)度的優(yōu)化目標(biāo)（即缺水總量）更優(yōu)，缺水總量相對(duì)差距0.12億m，差別較小，差距量占缺水總量（10.37億m3）的1.16%。該結(jié)論與黃強(qiáng)的研究結(jié)論一致（其發(fā)現(xiàn)采用獨(dú)立隨機(jī)過(guò)程或考慮相鄰時(shí)段徑流相關(guān)關(guān)系的簡(jiǎn)單馬氏過(guò)程對(duì)水庫(kù)優(yōu)化調(diào)度成果的影響不大），其原因主要為：考慮多維變量關(guān)聯(lián)性影響時(shí)，其核心是對(duì)隨機(jī)動(dòng)態(tài)規(guī)劃算法的遞推公式（公式17）中的狀態(tài)轉(zhuǎn)移概率產(chǎn)生影響，考慮二者關(guān)聯(lián)性影響時(shí)，狀態(tài)轉(zhuǎn)移概率{[Q(), W()]| [Q(),W()]}采用公式（12）計(jì)算，不考慮二者關(guān)聯(lián)性影響（即徑流與需水呈獨(dú)立關(guān)系）時(shí)，狀態(tài)轉(zhuǎn)移概率采用公式（14）計(jì)算，將狀態(tài)轉(zhuǎn)移概率{·}與時(shí)段優(yōu)化調(diào)度結(jié)果O(·)與余留期效益f(·)之和相乘，進(jìn)而遞推即可得到水庫(kù)調(diào)度優(yōu)化目標(biāo)結(jié)果。通過(guò)圖5所展示的轉(zhuǎn)移概率結(jié)果對(duì)比發(fā)現(xiàn)，考慮二者關(guān)聯(lián)性影響與不考慮二者關(guān)聯(lián)性影響時(shí)，狀態(tài)轉(zhuǎn)移概率{·}的差別較小，因此導(dǎo)致水庫(kù)調(diào)度優(yōu)化調(diào)度規(guī)則差別較小，因此將長(zhǎng)系列徑流與需水實(shí)測(cè)資料輸入該優(yōu)化調(diào)度規(guī)則時(shí)所得到的優(yōu)化調(diào)度結(jié)果差異較小。此外，從理論上講，徑流與需水同受氣象因素影響，且由于流域內(nèi)氣象條件的相似性，徑流與需水二者之間的相關(guān)性切實(shí)存在，難以忽視，故此，建議采用考慮多維變量關(guān)聯(lián)性的多維隨機(jī)動(dòng)態(tài)規(guī)劃算法求解含農(nóng)業(yè)灌溉用水的水庫(kù)優(yōu)化調(diào)度問(wèn)題。此外，若水庫(kù)以工業(yè)、生活用水等為目標(biāo)時(shí)，來(lái)水與工業(yè)、生活用水之間的關(guān)聯(lián)性較差，可采用本文中的公式（15）計(jì)算來(lái)水與工業(yè)、生活用水之間的狀態(tài)轉(zhuǎn)移概率，進(jìn)而應(yīng)用本文多維動(dòng)態(tài)規(guī)劃算法進(jìn)行求解。

4 結(jié) 論

水庫(kù)來(lái)水以及灌溉用水二維過(guò)程的隨機(jī)性與關(guān)聯(lián)性切實(shí)存在。本文利用Copula函數(shù)在表征多維隨機(jī)變量的非線性相關(guān)與隨機(jī)性方面的優(yōu)勢(shì)，描述同類(lèi)變量不同時(shí)段之間的隨機(jī)性以及不同變量相同時(shí)段之間的關(guān)聯(lián)性，在此基礎(chǔ)上，與傳統(tǒng)隨機(jī)動(dòng)態(tài)規(guī)劃算法結(jié)合，提出了考慮多維變量關(guān)聯(lián)性與不考慮多維變量關(guān)聯(lián)性的多維隨機(jī)動(dòng)態(tài)規(guī)劃算法。

以東莊水庫(kù)優(yōu)化調(diào)度為例，結(jié)果表明，采用考慮水庫(kù)來(lái)用水過(guò)程關(guān)聯(lián)性的多維隨機(jī)動(dòng)態(tài)規(guī)劃算法求解優(yōu)化調(diào)度模型，缺水總量為10.37億m，采用不考慮二者關(guān)聯(lián)性的多維隨機(jī)動(dòng)態(tài)規(guī)劃算法求解時(shí)，缺水總量為10.49 億m，缺水量相對(duì)降低0.12億m，考慮來(lái)用水過(guò)程關(guān)聯(lián)性的結(jié)果更優(yōu)，但差距較小，相對(duì)提升百分比為1.16%；考慮來(lái)用水過(guò)程關(guān)聯(lián)性時(shí)的狀態(tài)轉(zhuǎn)移概率較不考慮二者關(guān)聯(lián)性時(shí)的狀態(tài)轉(zhuǎn)移概率差異較小是造成優(yōu)化調(diào)度結(jié)果差異較小的主要原因?？紤]來(lái)用水過(guò)程關(guān)聯(lián)性的多維隨機(jī)動(dòng)態(tài)規(guī)劃算法理論上更符合水資源供需過(guò)程天然存在的隨機(jī)性與關(guān)聯(lián)性特征，本文建議采用考慮來(lái)用水過(guò)程關(guān)聯(lián)性的多維隨機(jī)動(dòng)態(tài)規(guī)劃算法求解水庫(kù)優(yōu)化調(diào)度模型。

若水庫(kù)以工業(yè)與生活供水為目標(biāo)，該類(lèi)用水過(guò)程與水庫(kù)來(lái)水過(guò)程之間并無(wú)機(jī)理上的關(guān)聯(lián)性，建議采用本文所提的不考慮隨機(jī)變量間關(guān)聯(lián)性的多維隨機(jī)動(dòng)態(tài)規(guī)劃算法求解。此外，本文所提的隨機(jī)動(dòng)態(tài)規(guī)劃算法在其他領(lǐng)域的隨機(jī)優(yōu)化決策問(wèn)題方面（如基于水庫(kù)調(diào)度的水風(fēng)光多能互補(bǔ)等）同樣適用。未來(lái)，如何擴(kuò)展二維隨機(jī)動(dòng)態(tài)規(guī)劃算法至更高維情形，是該項(xiàng)研究的重點(diǎn)所在。