一、節水灌溉製度優化
在供水量不受限製的條件下,灌溉製度的設計方法已在《農田水利學》中介紹,這裏不需贅述。但是當水源不足,無法充分滿足作物生育期灌溉水量時,如何確定節水灌溉製度,便是本部分要論述的內容。下麵以缺水地區旱作物為討論對象,僅介紹確定節水灌溉製度的動態規劃法,其他方法如模糊優化法和存貯模型等,可參閱有關文獻。
在幹旱缺水地區,水資源十分短缺,這將不可避免地引起農作物不同程度的減產,減產的程度隨著不同作物、不同生育階段的缺水程度而異。對於這類地區,合理的灌溉應是在弄清各種作物在不同生長時期缺水減產情況的基礎上實行限額灌溉,尋求最優的灌溉製度及灌溉用水方案,使全生長期或全灌區的總增產值最大,或總的減產值最小。
某種農作物灌溉製度優化設計過程是個多階段決策過程,可以用動態規劃的數學模型。
(1)階段變量。以農作物的生長階段為階段變量,i=1,2,…,N,其編號與階段初編號一致。
(2)決策變量。各生長階段的灌水量di,m3/hm2,及實際騰發量ETai,m3/hm2,i=1,2,…,N。
(3)狀態變量。各階段初可用於分配的水量qi,m3/hm2,及計劃濕潤層內可供農作物利用的土壤水量Si,m3/hm2,它是土壤含水率的函數。
Si=10000γ/γ水H(β一βω)(565)
式中,γ,γ水-土壤幹容重,t/m3;H-計劃濕潤層深度,m;β-計劃濕潤層內土壤平均含水率,以占幹土重的百分數計;βω-土壤含水率下限,約大於凋萎係數,以占幹土重的百分數計。
(4)係統方程。係統方程是描述係統在運動過程中狀態轉移的方程,由於本係統有兩個狀態變量,係統方程也有兩個:
其一是水量分配方程,若對第i個生長階段采用決策di時,可表達為:
qi+1=qi-di(566)
式中,qi、qi+1分別為第i及第i十1階段初可用於分配的水量,m3/hm2。
其二是土壤計劃濕潤層內的水量平衡方程,可寫成
ETai+Ei+Si+1=Si+di+Pi+Ki(567)
式中,Ei-第i階段的滲漏量,m3/hm2;Pi-第i階段的有效降雨量,m3/hm2;Si,Si+1-分別為第i及第i+1階段初土壤中可供利用的水量,m3/hm2;Ki-第i階段的地下水補給量,m3/hm2。
(5)目標函數。采用Jensen提出的在供水不足條件下,水量和農作物實際產量的連乘模型,目標函數為單位麵積的產量最大。
F=max(Ya/Ym)=maxNi=1(ETa/ETm)λii(568)
式中,Ya-作物實際產量(kg/hm2);Ym-供水充分條件下作物最大產量(kg/hm2);ETm-最大騰發量,m3/hm2;λi-第i生育階段農作物對缺水的敏感性指數。
(6)約束條件。
決策約束:0≤di≤qi,i=1,2,…,N(569)
∑Ni=1di=Q(570)
ETmin,i≤ETai≤ETmi(571)
式中,Q-全生長期單位麵積上可供分配的水量,m3/hm2;ETmin,i-第i階段最小需水量,m3/hm2。
土壤含水率約束:
βω≤β≤βf(572)
式中,βf為田間持水率,以占幹土重的百分數計。
(7)初始條件。假定作物播種時的土壤含水率為已知,即
β1=β0(573)
則有S1=10000γH(β0-βω)(574)
設第1時段初可用於分配的水量為農作物全生長期可用於分配的水量,即
q1=Q(575)
本模型是一個具有兩個狀態變量及兩個決策變量的二維動態規劃問題,可用多維動態規劃的求解方法求解,如用逐次漸近法(DPSA)求解。
二、農作物間水量分配優化
在一個灌區內.往往種植有幾種農作物,在水源不足的情況下,如何安排每一種作物的種植麵積,才能使全灌區獲得最大的經濟效益,這就是一定水量在幾種作物之間的最優分配問題,也就是最優種植比例問題。同時,在作物種植麵積已定的情況下,各農作物之間在全生長期及生長期內的各個時段也存在對水量的競爭,因而也有水量的最優分配問題。
係統分析方法中的許多優化技術,都可用於求解上述兩類水量分配問題,下麵列舉的數例,隻是為了說明數學模型的建立及求解,以揭示水量分配問題的實質,至於采用何種優化技術,則不受此限。
(一)農作物最優種植比例模型
某地區擬種植三種作物(冬小麥、棉花、夏玉米),對於一個特定的年份,由於灌溉水量的限製,不能隨意安排三種作物的種植麵積,可建立如下的數學規劃模型,尋求一定水量在三種作物之間的最優分配。
1.目標函數
以灌區純收入最大作為擇優準則,可寫成:
F=max(B-C)(576)
式中,F-最大純收入,萬元;B-毛收入,萬元;C-灌溉費用及其他農業費用,萬元。
(1)毛收入。對於K種作物而言,其農業毛收入可用下式表示:
B=∑Ki=1(Pi+P′iLi)YiXiQi+∑Ki=1(Pi+P′iLi)YoiAoi(577)
即
B=∑Ki=1(Pi+P′iLi)(YiXiQi+YoiAoi)(577)'
式中,Pi-第i種作物的實物單價,元/kg;Pi'-第i種作物的秸稈單價,元/kg;Li-第i種作物單位麵積上秸稈產量與實物產量的比值;Yi-第i種作物單位麵積上灌水量,如Qi(m3/hm2)的實物產量,kg/hm2;Yoi-第i種作物不灌水條件下的單位麵積產量,kg/hm2;Aoi-第i種作物不灌水部分的種植麵積,畝;Xi-第i種作物分得的有效灌溉水量,m3;Qi-第i種作物單位麵積灌水量,m3/hm2;K-作物種類數。
(2)費用。費用一般包括兩項,其一為不變費用,即農業的基本投入,包括種子、化肥、農藥、中耕除草……等費用,它不隨灌水量而變化;其二是可變費用,如灌水費用、勞工費用以及由於灌溉而增加的農業費用等,可用下式表示:
C=∑Ki=1(PωiQi+VciYci+Fci+Foi)XiQi+∑Ki=1FoiAoi(578)
式中,Pωi-第i種作物每立方米灌水的水費(井灌條件下含灌水費用及抽水動力費),元/m3;Vci-第i種作物單位產量的費用,元/kg;Yci-第i種作物單位麵積上的產量kg/hm2;Fci-第i種作物單位麵積上由於灌溉而增加的農業費用,假定與灌溉增產量呈直線關係,元/畝;Foi-第i種作物單位麵積上的基本投入(不變費用),元/hm2。
其中
Pωi=Pωi1+Pωi2(579)
Fci=Ci(Yi-Yoi)(580)
Vci·Yci=Ki[Qi/50+Ei(Yi-Yoi)/250](581)
將(579)~(581)代入(578),則有:
C=∑ki=1{(Pωi1+Pωi2)Qi+Ki[Qi/50+Ei(Yi-Yoi)/250+Ci(Yi-Yoi)+Foi]XiQi}+∑ki=1FoiAoi(582)
式中,Pωi1、Pωi2-分別為第i種作物單位灌水量的費用和抽水電費,元/m3;Qi/50-單位麵積澆地用工數(這裏假定每畝地灌水50m3用一個工日),d;Ei-第i種作物收250kg產品所需工日數,d;Ki-單位工日價格,元/d;Ci-由農作物單位增產而增加的農業費用,調查求得。
2.約束條件
(1)水量平衡約束。
∑ki=1Xi=V(583)
式中,V-該係統可用於灌溉的有效水量,m3。
(2)政策約束。若政策要求滿足糧食生產能自給自足,如按人均口糧400kg計,則有:
∑ki=1(YiXiQi+YoiAoi)≥TPN×400·r1(584)
r1+r2+…+rk1=1(585)
式中,TPN-係統內人口數量;ri-第i種糧食作物的權重;K1-糧食作物種類。
若對於經濟作物,也應控製在一定範圍之內,即有:
∑Ki=K1+1XjQj+Aoj≤Aj(586)
式中,Xj-分配給第j種經濟作物的灌溉水量,m3;Qj-單位麵積上第j種經濟作物的灌溉水量,m3/hm2;Aj-第j種經濟作物的最大種植麵積,畝;Aoj-第j種經濟作物在不灌水條件下的種植麵積,hm2。
(3)前後期連種作物的麵積約束。對於接茬作物,後期的種植麵積不應大於前期的種植麵積,可以寫成:
∑K1T=1(XTQT+A0T)≤∑K1L=1(XLQL+A0L)(587)
式中,K2-接茬作物種類數;角標L和T分別表示前期和後期作物。
(4)耕地麵積約束。任何時候,各種作物同時種植麵積之和不應大於總耕地麵積,即:
∑K1i=1(XiQi+Aoi)+∑Kj=K1+1(XjQj+Aoj)≤A(588)
式中,A-係統總耕地麵積,hm2。
對於這樣一個多種農作物的優化配水模型,由於Xi、Xj及Qi、Qj均為變量,實質上是一個非線性規劃問題,可用非線性規劃模型求解方法求解,也可設法將其變換成一係列線性規劃問題,再用修正單純形法分別求解。
(二)灌溉水量在多種作物的全生長期及生長期內的最優分配
1.數學模型
假定在生長期開始,可供分配給幾種作物的總水量為V0,各種作物的種植麵積為Ak(k=1,2,…,n)。總的目標就是把V0分配給幾種作物的種植麵積上,使之產生最大的經濟效益。
首先求解全生長期水量分配問題。為了考慮灌水時間對各種作物產量的影響,可從灌溉製度模型中導出全生長期的水分生產函數。而各種作物每周的灌水計劃,可以利用灌溉製度模型,由該作物分配到的全生長期的供水量推算出來。這個灌水計劃在每一周都以各種作物爭水的形式在屏幕上顯示出來。如在某一周發生了多種作物爭水現象,可啟動第二級模型求解。在這種情況下,所有農作物的灌水計劃對於該周及後續周次都要進行修改,這個步驟逐周依次重複進行,直到全生長期期結束。因此,整個數學模型由三個模型組成:①單一作物灌溉製度模型,對於給定的全生長期灌溉供水量,使各種作物的產量最高(模型1);②第一級模型,把全生長期的可供水量分配到各種作物(模型2);③第二級模型,在生長期內的任一時段,把適宜水量分配給各種作物(模型3)。模型1給模型2及模型3提供了輸入資料。以下將介紹這三個模型的組成及功能。
(1)單一作物灌溉製度模型-模型1。
某一種農作物的灌溉製度優化設計問題,通常可用動態規劃模型求解,如本節第一部分所述。為了獲得每周的灌溉用水計劃,單一作物模型應分兩步求解。第一步,用動態規劃使水分生產函數最大化,以求得各個生長階段的水量分配,在每個生育階段開始,可供水量及土壤含水量是兩個狀態變量,水量分配要滿足土壤水平衡方程及灌溉係統所具有的約束條件;第二步,把分配給每個生長階段的水量,再以順序的方式分配到組成該階段的每一周。第二步分配的基礎是田間試驗資料,在生育期早期的若幹周內,僅靠土壤原有含水量就可得到較優的產量,而在作物生長臨界期,則應及時灌水。例如,如果揚花期的灌水推遲了,其結果是開花較少,這就導致該階段給定用水水平下的產量降低,這也就是農學家們推薦在作物臨界生長期開始時就進行灌溉的道理,如分蘖期、揚花期灌水等。
這裏所用到的水分生產函數是:
YaYm=Ni=1[1-λi(1-ETa/ETm)i](589)
ETai=∑Mij=1eaij(590)
ETmi=∑Mij=1emij(591)
在上述各式中,N-生長階段數;Ya-實際的作物產量;Ym-最大的(ETa=ETm)作物產量;λi-第i個生長階段作物對缺水的敏感性指標;eaij及emij-分別為第i階段第j周實際的及潛在騰發量;Mi-第i階段的周數。
emij可用一定的模式,按蒸發皿蒸發量資料估算,而eaij則應從土壤水平衡模型中估算。模型假定均勻的土壤儲水層深度(Zij)等於所考慮的該階段作物有效根係活動層深度。Zij的值可從一個線性的根係生長模型計算出來。某階段的期望降雨量(rij)及灌溉水深(Uij)歸並在一起,並在時段開始時輸入到土壤水庫中。有了這些資料以後,模型可以根據土壤含水量及土壤耗水量因子P計算實際騰發量。
如果以Sij表示時段(i,j)末可供利用的土壤含水量,以Sc表示單位土層深度的最大土壤含水量,則:
eaij=emijZijSij≥(1-P)ScZij
eaij=SijZijemij(1-P)ScZij,ZijSij<(1-P)ScZij(592)
而Sij可用土壤水平衡模型求得:
Sij=min(Sij-1Zij-1+rij+Uij+S0ΔZij-eaij)/Zij
Sc(593)
在(593)中
ΔZij=Zij-Zi,j-1(594)
S0-生育期開始時,可供利用的土壤含水量,ΔZij-由於該階段根係活動層增加而增加的土壤含水量增量,當入滲水量(rij+Uij)超過了田間持水率時、多餘水量應排出根係活動層。
每個時段開始時的灌水深度(Uij)應滿足下列約束條件。如果在不灌溉時(但考慮了期望降雨量)土壤中可供利用的水量足以維護潛在蒸發量速率直到該時段結束,則灌水量等於零。否則,灌水量要受土壤貯水能力的製約,各時段灌水量之和不應超過該階段所分配到的水量(Xi),或受該時段渠係配水能力(AWCij)的約束,這就是:
Uij=0ZijSij≥(1-P)ScZij+emij
Uij=minScZij-Si,j-1Zi,j-1-rij-S0ΔZij
(可供利用的土壤含水量)
Xi-∑j-1i=1Uit
(可供分配的供水量)
AWCij
(渠道供水能力)
ZijSij<P=ScZij+epij(595)
假定Q0為生長期開始時可供利用的水層深度,Xi是分配到每個生長階段的水量,則:
Q0≥∑Ni=1Xi(596)
Xi≥∑Mii=1Uij(597)
利用動態規劃法使(589)式最大化,其遞推方程如下:
fi(Q,S)=max[1-λi(1-ETa/ETm)i=f*i+1(Q-XiSi+1)(598)
0≤Xi≤Q,0≤Q≤Q0,i=N-1,N-2,…,1
Si+1=Si,Mi
fN(Q,S)=max[1-λNETa/ETm]N(599)
0≤XN≤Q,0≤Q≤Q0
對於特定的Q0及S0,使(598)及(599)式最大化,並滿足(590)~(597)式,可得到生長階段i的最優水量分配Xi*(i=1,2,…,N)。對於整個生長期,以Xi=Xi*運行土壤水平衡模型(592)~(595)式,可得到每周的最優配水量Uij,j=1,2,…,Mij,i=1,2,…,N。
(2)多種作物全生長期的水量最優分配-模型2。
1)各種作物全生長期的水分生產函數。模型1可以在一定範圍內(0≤Vk/Ak≤V0/Ak)對每種作物,對各種生育期灌水定額(Vkl/Ak,l=1,2,…,m)求解。在(598)中,使Q0=V0/Ak,用模型1的動態規劃程序,求解單一作物模型,可得到相應於一組生育期灌水量(X*kl=Vkl/Ak;l=1,2,…,m)的一組作物相對產量(Y*kl;l=1,2,…,m)。所得到的這些組數據(Y*kl,X*kl)之間的關係就構成該作物生育期的水分生產函數,它們可以是二次函數、指數函數或其他形式的函數,可通過曲線擬合求得。如采用指數函數,對每種作物k,可建立如下的相對產量(Yk)與生育期灌水量(Xk)的關係:
Yk-Y0k=(ak-Y0k)exp[-bk(X0k-Xk)/Xk],k=1,…,n(5100)