陳紅光 王瓊雅 李曉寧 王中君 李晨陽

(1.東北農(nóng)業(yè)大學(xué)水利與建筑學(xué)院， 哈爾濱 150030； 2.東北農(nóng)業(yè)大學(xué)工程學(xué)院， 哈爾濱 150030)

0 引言

隨著社會(huì)經(jīng)濟(jì)和人口的快速增長(zhǎng)，水資源短缺問題日益突出。隨著糧食需求不斷增加，農(nóng)業(yè)作為用水大戶，水資源供需問題不斷加重[1-2]。合理配置農(nóng)業(yè)用水、提高農(nóng)業(yè)用水效率對(duì)緩解水資源供需矛盾和用水壓力具有十分重要的意義[3-5]。灌區(qū)作為農(nóng)業(yè)用水主體，其水資源的合理優(yōu)化配置對(duì)提高水資源利用率、促進(jìn)農(nóng)業(yè)水資源的可持續(xù)開發(fā)利用、權(quán)衡灌區(qū)水資源經(jīng)濟(jì)和風(fēng)險(xiǎn)之間的關(guān)系具有重要作用。

近年來，水資源優(yōu)化配置得到了廣泛發(fā)展，出現(xiàn)了多種優(yōu)化配置方法。婁帥等[6]基于免疫遺傳算法構(gòu)建多階段群決策優(yōu)化模型，解決漳河流域水資源優(yōu)化配置問題。孫冬營(yíng)等[7]運(yùn)用模糊聯(lián)盟合作博弈方法將水資源進(jìn)行兩次分配，得到流域水資源的合理配置方案。陳述等[8]將粒子群人工蜂群算法運(yùn)用到灌區(qū)渠-塘-田的水資源優(yōu)化配置中。上述研究基于確定的水資源系統(tǒng)進(jìn)行配置，但在實(shí)際中，灌區(qū)灌溉系統(tǒng)受到氣象條件、作物需水量以及政策變化等因素影響，是一個(gè)復(fù)雜的動(dòng)態(tài)配水系統(tǒng)[9]。傳統(tǒng)方法在處理系統(tǒng)中的不確定因素時(shí)，存在一定的局限性。因此，為了解決這些包含不確定因素的問題，提出了不確定優(yōu)化方法，其中兩階段隨機(jī)規(guī)劃模型得到了廣泛應(yīng)用[10]。HUANG等[11]將不確定性優(yōu)化和兩階段隨機(jī)規(guī)劃方法應(yīng)用到水資源系統(tǒng)管理中；LI等[12]將區(qū)間兩階段機(jī)會(huì)約束模型運(yùn)用到水資源規(guī)劃中；ZHANG等[13-14]運(yùn)用改進(jìn)后的區(qū)間兩階段隨機(jī)規(guī)劃算法，對(duì)三江平原水資源進(jìn)行優(yōu)化配置。區(qū)間兩階段優(yōu)化算法在大量學(xué)者[15-17]的努力下，逐漸得到完善，用以處理水資源配置過程中的不確定性問題。區(qū)間兩階段隨機(jī)效用模型處理不確定因素十分有效，但卻忽略了系統(tǒng)風(fēng)險(xiǎn)問題，模型結(jié)果不具有絕對(duì)可行性。而魯棒優(yōu)化方法則能夠在規(guī)劃過程中對(duì)風(fēng)險(xiǎn)進(jìn)行有效規(guī)避，并權(quán)衡系統(tǒng)中可變隨機(jī)值和追索成本的關(guān)系[18-20]。

本文針對(duì)灌區(qū)水資源配置模型中存在風(fēng)險(xiǎn)的特點(diǎn)，構(gòu)建多水源、多作物配水模型，以兩階段線性規(guī)劃為基礎(chǔ)，建立區(qū)間兩階段魯棒優(yōu)化模型(Interval-parameter two-stage robust stochastic programming model，ITRM)，采用概率密度函數(shù)、離散區(qū)間表示系統(tǒng)的不確定性，用魯棒系數(shù)表示系統(tǒng)的風(fēng)險(xiǎn)。在不同來水情況下，對(duì)各個(gè)灌區(qū)的不同作物進(jìn)行水資源的優(yōu)化配置，以期為灌區(qū)管理者提供風(fēng)險(xiǎn)可控、成本最優(yōu)的配水方案。

1 模型建立與求解

1.1 模型建立

灌區(qū)多水源優(yōu)化配置是一個(gè)復(fù)雜的系統(tǒng)過程，供水情況受來水量等因素影響，存在一定的隨機(jī)特點(diǎn)[21]。且多水源聯(lián)合優(yōu)化調(diào)度目的是使水資源得到合理的分配，滿足灌區(qū)用水需求，同時(shí)盡量減少用水成本。因此，本文以灌區(qū)農(nóng)作物需水量為決策變量，引入用水成本、缺水懲罰系數(shù)，并分兩個(gè)階段確定灌區(qū)水資源的最優(yōu)配置。在保證灌區(qū)水資源承載力的前提下，以正常水平年作物需水量為依據(jù)，確定作物預(yù)先配水目標(biāo)，并將其作為第1階段的決策變量；由于受灌區(qū)初始儲(chǔ)水量、天然來水量、水源蒸發(fā)量等因素的影響，供水量可能小于供水目標(biāo)，這就需要決策者進(jìn)行調(diào)整，減少灌區(qū)供水量或者外調(diào)水源進(jìn)行補(bǔ)水，而減少供水量會(huì)影響作物產(chǎn)量，外調(diào)水源則會(huì)增加用水成本，都會(huì)產(chǎn)生經(jīng)濟(jì)懲罰，為降低用水成本，需要對(duì)第1階段的配水量進(jìn)行調(diào)整，將缺水量作為第2階段的決策因素。故兩階段隨機(jī)規(guī)劃模型為

(1)

式中Cij——水源i向作物j的輸水成本，元/m3

Wij——水源i向作物j的預(yù)先目標(biāo)配水量，m3

E——隨機(jī)變量的期望值

Dij——水源i未滿足作物j預(yù)先目標(biāo)配水量時(shí)的缺水懲罰系數(shù)，元/m3

Sij——水源i未滿足作物j預(yù)先目標(biāo)配水量時(shí)的缺水量，m3

i——水源，i為1、2分別表示地表水和地下水

j——作物，j為1、2、3分別表示水稻、玉米、大豆

f——系統(tǒng)總成本，元

(2)

式中Sijh——來水量水平為h時(shí)，i水源向j作物供水，未滿足預(yù)先目標(biāo)配水量的缺水量，m3

已有的研究較少考慮系統(tǒng)風(fēng)險(xiǎn)問題，無法保證模型最優(yōu)解的絕對(duì)可行性。魯棒優(yōu)化方法將風(fēng)險(xiǎn)體現(xiàn)在函數(shù)當(dāng)中，對(duì)風(fēng)險(xiǎn)進(jìn)行有效測(cè)評(píng)，并在規(guī)劃過程中規(guī)避風(fēng)險(xiǎn)，平衡系統(tǒng)的成本和風(fēng)險(xiǎn)的關(guān)系，可以有效增加模型求解的可行性和系統(tǒng)的穩(wěn)定性[22]。因此，本文在兩階段隨機(jī)規(guī)劃模型的基礎(chǔ)上引入魯棒優(yōu)化方法，建立兩階段隨機(jī)魯棒優(yōu)化模型。模型表述為

(3)

式中ρ——魯棒系數(shù)

(4)

式中θih——松弛變量

θih主要作用是保證模型的穩(wěn)定性和非負(fù)性。作物預(yù)先目標(biāo)配水量(需水量)Wij具有不確定性，且作物的價(jià)格和產(chǎn)量變動(dòng)導(dǎo)致懲罰系數(shù)Dij也不確定。為表示不確定性，引入?yún)^(qū)間參數(shù)來表示不確定性參數(shù)。“+”表示參數(shù)的上限，“-”表示參數(shù)的下限，建立區(qū)間兩階段魯棒優(yōu)化模型(ITRM)

(5)

(1)需水量約束

(6)

式中Wijmax——作物j正常生長(zhǎng)的最大需水量

Wijmin——作物j正常生長(zhǎng)的最小需水量

(2)水源可用水量約束

地表水約束為

(7)

地下水約束為

(8)

式中Qij——灌區(qū)初期水源i的儲(chǔ)水量，m3

Qim——灌區(qū)末期水源i的蓄水量，m3

不錯(cuò)，《罕哈冉惠傳》中出現(xiàn)了不少反映佛教思想的內(nèi)容。如學(xué)者們指出的，當(dāng)哈冉惠和他的兩個(gè)兄弟遭蟒古斯暗害，誤食了蟒古斯投下的毒藥后，哈冉惠變成了長(zhǎng)有九十五顆頭的大黑蟒古斯，他的胞弟烏蘭岱·莫日根變成了一尊石人，烏蘭岱的坐騎變成了一尊石馬，他的義兄弟吉爾吉斯·賽因·貝托爾則變成了一頭黃色的野豬。是兩位公主派出了維蘭·索龍嘎的兒子去向菩薩求救，菩薩將三件寶物——萬能的金套繩、能起死回生的仙丹、智慧的金盤賜給了小童，小童用此三寶使三位英雄得救，恢復(fù)了原貌。菩薩，為眾人所知，當(dāng)然屬于重要的佛的形象，然而當(dāng)菩薩向小童授予三件寶物之前，卻說出下面一段頗令人費(fèi)解的話：

qih——水源i的來水量，m3

Qsi——灌區(qū)水源i的蒸發(fā)量，m3

Qimin——水源i最低蓄水量，m3

(3)追索變量約束

追索變量約束為

(9)

(4)水源最大供水量約束

水源最大供水量約束為

(10)

式中Wimax——水源i最大可供水量，m3

(5)非負(fù)約束

非負(fù)約束為

(11)

1.2 模型求解

將ITRM分為兩個(gè)子模型進(jìn)行求解，對(duì)應(yīng)的目標(biāo)函數(shù)下限子模型為

(12)

約束條件為

(13)

(14)

約束條件為

(15)

(16)

(17)

zij=zijopt(?i,j)

(18)

最優(yōu)配水量為

(19)

式中Aij——水源i向作物j的最優(yōu)配水量

其算法流程如圖1所示。

圖1 區(qū)間兩階段魯棒優(yōu)化方法算法流程圖Fig.1 Algorithm flow chart of interval-parameter two-stage robust stochastic programming method

2 案例研究

2.1 區(qū)域概況

三江平原地區(qū)牡丹江灌區(qū)位于黑龍江省東南部(131°13′～133°37′E，44°57′～46°12′N)，耕地面積為3.07×109m2，所轄12個(gè)現(xiàn)代化農(nóng)場(chǎng)。灌區(qū)內(nèi)土壤肥沃，水資源豐富，主要農(nóng)作物為水稻、大豆、玉米，以井灌為主。然而，近些年隨著作物種植面積的不斷增加，作物需水量增大，導(dǎo)致灌區(qū)用水量壓力增大，水資源危機(jī)加重。實(shí)際配置水資源量如果小于農(nóng)作物的需水量,供水不足將造成農(nóng)作物減產(chǎn)，增大懲罰成本，若滿足農(nóng)作物需水量，則又增加用水成本和系統(tǒng)供水風(fēng)險(xiǎn)。因此，本文運(yùn)用ITRM有效的平衡系統(tǒng)成本和風(fēng)險(xiǎn)，解決灌區(qū)水資源分配問題。

2.2 參數(shù)確定

選取水稻、大豆、玉米3種作物為研究對(duì)象，根據(jù)《牡丹江統(tǒng)計(jì)年鑒》[24]、《黑龍江墾區(qū)統(tǒng)計(jì)年鑒》[25]等相關(guān)資料可以得到3種作物的灌區(qū)種植面積、單位面積灌溉用水區(qū)間值，作物灌溉面積比例參照文獻(xiàn)[27]，水稻、大豆、玉米分別為100%、10%、10%，根據(jù)文獻(xiàn)[24-25]和當(dāng)?shù)囟嗄杲y(tǒng)計(jì)數(shù)據(jù)，得到作物單位面積最小灌溉量和最大灌溉量(表1)。以灌區(qū)2011年3種作物種植面積數(shù)據(jù)為已知條件，將作物單位面積最小(最大)用水量乘以作物種植面積、作物灌溉面積比例，得到作物最小(最大)需水量，查閱文獻(xiàn)[24-25]得到灌區(qū)水量分配的相關(guān)經(jīng)濟(jì)系數(shù)(表2)。根據(jù)《黑龍江省農(nóng)墾水利志》[26](1987—2015年)對(duì)灌區(qū)多年天然降水量的統(tǒng)計(jì)結(jié)果，可知1992、2002、2008、2011年灌區(qū)的降水量為400～450 mm，1987、1990、2009、2013年降水量為600 mm以上，其他年份降水量在450～600 mm之間。把來水量水平分為高、中、低3個(gè)水平，由統(tǒng)計(jì)結(jié)果可知，中流量年份比高降流量和低流量年份多，且高流量和低流量出現(xiàn)概率相近，故假設(shè)預(yù)測(cè)年份的降水量水平高、中、低出現(xiàn)的概率為0.2、0.6、0.2。查閱文獻(xiàn)[27]獲得地下水和地表水的供水比例為4∶1。其中地下水的主要補(bǔ)給來源為側(cè)向徑流補(bǔ)給[28]，2011年全年牡丹江地區(qū)地表水過境水量為2.482×1010m3，但其提水量2.3×107m3，僅占過境水量的0.093%。根據(jù)《牡丹江統(tǒng)計(jì)年鑒》[24]、《黑龍江墾區(qū)統(tǒng)計(jì)年鑒》[25]以及多年數(shù)據(jù)統(tǒng)計(jì)分析，得到灌區(qū)在不同來水量水平下的可供水量區(qū)間值(表3)，灌區(qū)地表水最大可供水量為6.2×108m3，最小蓄水量為1.1×108m3，地下水最大可供水量為9.7×108m3，最小儲(chǔ)蓄量為2.55×108m3。

表1 灌區(qū)作物種植面積及單位面積灌溉用水量Tab.1 Irrigation crops area and water per unit in irrigation district

表2 灌區(qū)作物最大最小需水量、預(yù)先配水量和經(jīng)濟(jì)系數(shù)Tab.2 Irrigation corps maximum and minimum water distribution, target value of crops water distribution in advance and economic parameters

表3 不同來水量水平下農(nóng)作物灌溉水源可用水量Tab.3 Water available for crops irrigation under different water inflow levels

2.3 模型計(jì)算結(jié)果與分析

圖2為灌區(qū)用于水稻、玉米、大豆的最優(yōu)目標(biāo)配水量，其中地表水最優(yōu)目標(biāo)配水量為6.39×107m3，地下水為2.696 8×108m3，在中水平流量下灌區(qū)可用的地表水、地下水資源量依次為[2.25×108,2.42×108] m3、[3.25×108，3.47×108] m3，此時(shí)配水量充足，不存在缺水。而在低流量情況下，灌區(qū)的地下水可用水量為[2.58×108,2.76×108] m3，此時(shí)可分配水量幾乎達(dá)到上限，作物供水存在缺失。可能的原因是:①水資源不合理開發(fā)，過度利用地下水。②缺少合理的規(guī)劃和利用。

最優(yōu)目標(biāo)配水量由模型第1階段得到，最優(yōu)分配水量則由式(19)確定。魯棒優(yōu)化方法旨在第2階段中對(duì)成本期望進(jìn)行追索，并對(duì)系統(tǒng)的風(fēng)險(xiǎn)進(jìn)行評(píng)估。表4中給出了ρ取不同值時(shí)，對(duì)應(yīng)的缺水量和最優(yōu)配水量。當(dāng)ρ=0時(shí)，模型為普通區(qū)間兩階段隨機(jī)規(guī)劃模型，代表了決策者對(duì)風(fēng)險(xiǎn)持中立態(tài)度，不考慮成本可變性。由表4可知：①對(duì)于水稻，地表水的決策變量為0，當(dāng)ρ=0時(shí)，在低、中、高3種來水水平下，其缺水量分別為[3.68×106，5.68×106] m3、[1.53×106，3.53×106] m3、[2.10×105，1.21×106] m3，因此，最優(yōu)配水量為[4.419×107，4.619×107] m3、[4.634×107，4.834×107] m3、[4.866×107，4.966×107] m3。由表3可知，水稻的地表水最小需水量為[4.321×107，4.444×107] m3,可滿足水稻的供水需求。對(duì)于玉米，地表水的決策變量為1，當(dāng)ρ=0時(shí)不產(chǎn)生缺水量，最優(yōu)配水量為6.03×106m3。在供水量充足的情況下，系統(tǒng)更偏好避免懲罰風(fēng)險(xiǎn)，所以將目標(biāo)配水量上限作為最優(yōu)目標(biāo)配水量。②隨著ρ的增大，缺水量不斷增加。例如，在低水平流量下，當(dāng)ρ=0.4，1，2時(shí)，水稻的地下水缺水量為[2.179×107，2.788×107] m3、[3.305×107，3.877×107] m3、[3.628×107，4.356×107] m3，最優(yōu)配水量為[1.716 1×108，1.767 0×108] m3、[1.607 2×108，1.664 4×108] m3、[1.559 3×108，1.632 1×108] m3，玉米的地表水缺水量為0 m3、[8.3×105，9.7×105] m3、[1.02×106，1.34×106] m3，最優(yōu)配水量為6.03×106m3、[5.06×106，5.20×106] m3、[4.69×106，5.01×106] m3。低水平時(shí)系統(tǒng)來水量較少，通過減少配水量來調(diào)節(jié)系統(tǒng)的最優(yōu)配水量，即，魯棒系數(shù)越大，系統(tǒng)的穩(wěn)定性越強(qiáng)。模型增加了作物缺水量，降低了作物的配水量，從而使系統(tǒng)的穩(wěn)定性增加。③當(dāng)ρ逐漸增大時(shí)，作物的缺水量隨之增大，但當(dāng)ρ≥3之后，缺水量幾乎不再變化。例如，當(dāng)ρ=2,3,5時(shí)，在中水平流量下，大豆的地表水缺水量為[8.80×105，1.19×106] m3、[1.32×106，2.58×106] m3、[1.36×106，2.60×106] m3，最優(yōu)配水量為[6.81×106，7.12×106] m3、[5.42×106，6.68×106] m3、[5.46×106，6.64×106] m3，玉米的地下水缺水量為[6.23×106，6.88×106] m3、[7.22×106，8.67×106] m3、[7.22×106，8.67×106] m3，最優(yōu)配水量為[1.723×107，1.788×107] m3、[1.544×107，1.689×107] m3、[1.544×107，1.689×107] m3，缺水量的增加會(huì)使系統(tǒng)穩(wěn)定性增加，但為了保證作物正常生長(zhǎng)的最小需水量，缺水量不再增大，此時(shí)系統(tǒng)趨于穩(wěn)定。

表4 不同ρ取值的作物缺水量和最優(yōu)配置水量Tab.4 Crops water shortage and water optimal allocation of different values of ρ

續(xù)表4

表5 決策變量值Tab.5 Value of decision variable

圖2 作物最優(yōu)目標(biāo)配水量Fig.2 Optimized allocation targets for different crops

最優(yōu)系統(tǒng)成本區(qū)間如圖3、4所示，系統(tǒng)成本隨著魯棒系數(shù)的變化而變化，呈遞增趨勢(shì)。如圖3所示，當(dāng)ρ=0時(shí)，最優(yōu)系統(tǒng)總成本為[1.104 32×109,2.049 95×109]元，當(dāng)ρ=1時(shí)，模型最優(yōu)系統(tǒng)總成本為[1.331 55×109，2.235 76×109]元。當(dāng)ρ=5時(shí)，最優(yōu)系統(tǒng)總成本為[1.943 77×109，2.657 69×109]

元。在低流量水平下(圖4a)，最優(yōu)系統(tǒng)成本在[6.835 5×108，9.286 7×108]元(ρ=5)和[3.693 9×108，6.962 5×108]元(ρ=0)之間變化。在中流量水平下(圖4b)，最優(yōu)系統(tǒng)成本在[6.496 7×108，9.130 5×108]元(ρ=5)和[3.683 3×108，6.864 5×108]元(ρ=0)之間變化。在高流量水平下(圖4c)，系統(tǒng)成本在[6.105 5×108，8.160 2×108]元(ρ=5)和[3.666 1×108，6.672 5×108]元(ρ=0)之間變化。對(duì)比圖3、4，隨著水資源最優(yōu)分配量的變化，系統(tǒng)成本呈現(xiàn)一定的變化規(guī)律：①魯棒系數(shù)增加，引起系統(tǒng)成本增大，當(dāng)ρ≥3之后，成本幾乎不變，說明系統(tǒng)已經(jīng)趨于穩(wěn)定。②隨著魯棒系數(shù)增加，成本的上下限差值變小，系統(tǒng)的穩(wěn)定性增加，經(jīng)濟(jì)性和穩(wěn)定性得到了較好的平衡。③較高的成本對(duì)應(yīng)較高的缺水水平。當(dāng)可用水量較高時(shí)，決策者可利用的水資源量也會(huì)較多，如果實(shí)際分配量較少，則會(huì)產(chǎn)生較高的系統(tǒng)風(fēng)險(xiǎn)和較多的懲罰成本；相反，如果可用水量較少，則決策者需要減少實(shí)際供水量，采取保守決策，降低系統(tǒng)風(fēng)險(xiǎn)增加穩(wěn)定性。決策者可以根據(jù)系統(tǒng)分析結(jié)果，針對(duì)灌區(qū)實(shí)際情況，制定風(fēng)險(xiǎn)和經(jīng)濟(jì)相協(xié)調(diào)的水資源配置策略。

圖3 最優(yōu)系統(tǒng)總成本Fig.3 Optimized net system cost

圖4 不同流量水平下最優(yōu)系統(tǒng)成本Fig.4 Optimized net system cost under different water inflow levels

3 結(jié)論

(1)針對(duì)水資源分配過程中存在風(fēng)險(xiǎn)的問題，將魯棒優(yōu)化方法與兩階段規(guī)劃方法耦合，建立了區(qū)間兩階段魯棒優(yōu)化模型，以三江平原牡丹江灌區(qū)農(nóng)業(yè)水資源配置為例進(jìn)行了研究。模型結(jié)果表明，系統(tǒng)總成本隨著魯棒系數(shù)的變化有一定的規(guī)律。當(dāng)模型不考慮系統(tǒng)風(fēng)險(xiǎn)時(shí)，即ρ=0時(shí)，系統(tǒng)成本在[1.104 32×109,2.049 95×109]元之間變化，隨著魯棒系數(shù)的增大，模型的缺水量增加，使得系統(tǒng)的穩(wěn)定性增強(qiáng)、成本增加，當(dāng)ρ=0.4、1、2時(shí)，系統(tǒng)成本分別在[1.263 25×109，2.185 67×109]元、[1.331 55×109,2.235 76×109]元、[1.608 79×109,2.415 52×109]元之間變化，但當(dāng)ρ=3、5時(shí)，系統(tǒng)缺水量不再增大，此時(shí)系統(tǒng)達(dá)到穩(wěn)定狀態(tài)，成本在[1.903 27×109,2.634 75×109]元、[1.943 77×109,2.657 69×109]元之間變化。由此可知，在不同來水量水平下，通過增加魯棒系數(shù)增加模型結(jié)果的可行性，對(duì)灌區(qū)的用水成本、系統(tǒng)的穩(wěn)定性和系統(tǒng)的風(fēng)險(xiǎn)三者之間進(jìn)行了充分的權(quán)衡，使配置方案更具有實(shí)際操作性和靈活性。

(2)與傳統(tǒng)區(qū)間兩階段隨機(jī)規(guī)劃模型相比，區(qū)間兩階段魯棒優(yōu)化方法不但可以有效地解決不確定條件下的隨機(jī)問題和區(qū)間問題，魯棒優(yōu)化方法的加入可以捕獲規(guī)劃過程中產(chǎn)生的風(fēng)險(xiǎn)問題，避免優(yōu)化結(jié)果出現(xiàn)高風(fēng)險(xiǎn)狀態(tài)，彌補(bǔ)了模型存在風(fēng)險(xiǎn)的缺陷，增強(qiáng)了系統(tǒng)的穩(wěn)定性。將區(qū)間兩階段魯棒優(yōu)化方法應(yīng)用到灌區(qū)水資源優(yōu)化配置中，驗(yàn)證了模型的應(yīng)用性和有效性。模型結(jié)果表明，通過魯棒系數(shù)的變化，生成一系列對(duì)應(yīng)的不同風(fēng)險(xiǎn)水平、不同情境的可行性方案，顯示了系統(tǒng)經(jīng)濟(jì)性和系統(tǒng)穩(wěn)定性之間的權(quán)衡。