三種群捕食系統(tǒng)動(dòng)力行為與最優(yōu)捕獲策略

馮慶紅, 陳榮軍, 馮艷青, 王忠英, 陳立范

(1. 常州工學(xué)院理學(xué)院, 江蘇 常州 213032; 2. 上海健康醫(yī)學(xué)院文理教學(xué)部, 上海 201318)

B-D功能反應(yīng)函數(shù)考慮了捕食者之間的相互競(jìng)爭(zhēng), 比Holling反應(yīng)函數(shù)[1]更加貼近實(shí)際．對(duì)帶B-D功能反應(yīng)函數(shù)的食餌捕食模型, 眾多學(xué)者分析了各種因素對(duì)系統(tǒng)的影響, 如擴(kuò)散[2-3]、食餌避難所[4-5]、脈沖反饋控制[6]、分支[7-8]、Allee效應(yīng)[9]、時(shí)滯[10-11]、 階段結(jié)構(gòu)[12-13]以及捕獲[14]等．由于捕獲問(wèn)題既關(guān)系到經(jīng)濟(jì)效益, 又與保護(hù)物種和生態(tài)環(huán)境有關(guān), 自然引起了生物數(shù)學(xué)學(xué)者的極大關(guān)注．最近, Haugen等[14]提出一種基于優(yōu)化的拖網(wǎng)作業(yè)規(guī)劃系統(tǒng), 建立了一個(gè)最優(yōu)控制問(wèn)題的公式, 為漁船和漁網(wǎng)自主捕魚(yú)找到了可行的軌跡; He等[15]建立了一類具全局反饋邊界條件和年齡結(jié)構(gòu)的非線性偏微分方程模型,并利用緊致性和最大化序列確定了最優(yōu)捕獲策略．受上述工作啟發(fā), 本文擬討論帶B-D功能反應(yīng)函數(shù)的三種群捕食系統(tǒng), 并對(duì)系統(tǒng)中2個(gè)競(jìng)爭(zhēng)捕食者種群同時(shí)進(jìn)行捕獲的最優(yōu)策略進(jìn)行分析．

1 三種群捕食系統(tǒng)

三種群捕食系統(tǒng)為

(1)

其中u,v,ω分別表示t時(shí)刻食餌種群和兩競(jìng)爭(zhēng)捕食者種群的密度;r表示食餌種群的內(nèi)稟增長(zhǎng)率;k表示食餌種群的最大環(huán)境容納量;αv/(a+bu+cv)表示B-D功能反應(yīng)函數(shù);d表示捕食者種群v的死亡率;γ表示捕食者種群ω的內(nèi)稟增長(zhǎng)率;μ表示捕食者種群ω和食餌種群u在平衡狀態(tài)下生物量之比;q2,q3分別表示對(duì)捕食者種群v和ω的可捕系數(shù);E2,E3分別表示相應(yīng)的捕獲努力量．假設(shè)所有參數(shù)為正數(shù), 顯然γ-q3E3>0．據(jù)生態(tài)學(xué)意義, 僅在區(qū)域{(u,v,ω)|u>0,v>0,ω>0}=R+內(nèi)對(duì)系統(tǒng)(1)進(jìn)行討論．

2 正平衡點(diǎn)的存在性與唯一性

定理1系統(tǒng)(1)存在唯一正的平衡點(diǎn)．

3 正平衡點(diǎn)的局部漸近穩(wěn)定性

定理2當(dāng)cβ>bα,r/k>bα,γb>μβ,c>1時(shí),正平衡點(diǎn)局部漸近穩(wěn)定．

證明 令f1(u,v,ω)=ru(1-u/k)-αuv/(a+bu+cv)-muω,f2(u,v,ω)=v[-d+βu/(a+bu+cv)]-q2E2v,f3(u,v,ω)=γω[1-ω/(μu)]-q3E3ω, 則系統(tǒng)(1)在正平衡點(diǎn)(u*,v*,ω*)處的Jacobian矩陣為

設(shè)對(duì)應(yīng)的特征值方程λ3+Dλ2+Fλ+G=0, 可知:G>0; 當(dāng)cβ>bα?xí)r,D>0; 且當(dāng)cβ>bα,r/k>bα,γb>μβ,c>1時(shí),DF-G<0. 根據(jù)Routh-Hurwitz判別法, 在上述條件下, 矩陣A的所有特征值都有負(fù)實(shí)部, 因此正平衡點(diǎn)局部漸近穩(wěn)定．

4 經(jīng)濟(jì)平衡點(diǎn)

假設(shè)兩競(jìng)爭(zhēng)捕食者種群v,ω的捕撈成本分別為常數(shù)c2,c3, 出售單價(jià)分別為常數(shù)p2,p3, 通過(guò)捕獲獲得的經(jīng)濟(jì)效益分別為L(zhǎng)2=(p2q2v-c2)E2和L3=(p3q3ω-c3)E3, 則總捕獲效益L=(p2q2v-c2)E2+(p3q3ω-c3)E3．顯然, 當(dāng)L2>0,L3>0時(shí), 捕獲效益大于捕獲費(fèi)用, 可以對(duì)兩競(jìng)爭(zhēng)捕食者種群進(jìn)行捕獲．

經(jīng)濟(jì)平衡點(diǎn)不同于生物平衡點(diǎn), 它是在經(jīng)濟(jì)利潤(rùn)被完全消耗時(shí)獲得的平衡點(diǎn)(u∞,v∞,ω∞,E2∞,E3∞), 滿足

解此方程組可得v∞=c2/(p2q2),ω∞=c3/(p3q3)．當(dāng){ar/k+rcc2/(kp2q2)-b[mc3/(p3q3)-r]}2>4brαc2/(kp2q2)時(shí),u∞存在, 且E2∞=βu∞/[q2(a+bu∞+cv∞)]-d/q2,E3∞=γ[1-ω∞/(μu∞)]/q3．由此可得如下定理．

定理3當(dāng){ar/k+rcc2/(kp2q2)-b[mc3/(p3q3)-r]}2>4brαc2/(kp2q2)時(shí), 系統(tǒng)的經(jīng)濟(jì)平衡點(diǎn)(u∞,v∞,ω∞,E2∞,E3∞) 存在．

5 最優(yōu)捕獲策略

最優(yōu)捕獲問(wèn)題可以轉(zhuǎn)化為如下的最優(yōu)控制問(wèn)題, 其目標(biāo)函數(shù)為max{J},使得

(2)

其中Ei(i=2,3)為捕獲控制變量, 且0≤Ei≤Emax．

定理4系統(tǒng) (1) 滿足條件{ar/k+rcc2/(kp2q2)-b[mc3/(p3q3)-r]}2>4brαc2/(kp2q2)時(shí), 存在唯一最優(yōu)平衡解(uδ,vδ,ωδ); 相應(yīng)地,對(duì)兩競(jìng)爭(zhēng)捕食者種群v,ω的最優(yōu)捕獲努力量分別為E2δ=βuδ/[q2(a+buδ+cvδ)]-d/q2和E3δ=γ[1-ωδ/(μuδ)]/q3．

證明 首先, 構(gòu)建此問(wèn)題的Hamilton函數(shù)H=e-δt[(p2q2v-c2)E2+(p3q3ω-c3)E3]+λ2{v[-d+βu/(a+bu+cv)]-q2E2v}+λ3{γω[1-ω/(μu)]-q3E3ω}, 其中λ1(t),λ2(t),λ3(t)為伴隨變量, 伴隨方程為

dλ1/dt=-?H/?u=-λ1[r(1-2u/k)-αv(a+cv)/(a+bu+cv)2-mω]-
λ2βv(a+cv)/(a+bu+cv)-λ3γω2/(μu2),
dλ2/dt=-?H/?v=-e-δtp2q2E2+λ1αu(a+bu)/(a+bu+cv)2-
λ2[-d+βu(a+bu)/(a+bu+cv)2-q2E2],
dλ3/dt=-?H/?ω=-e-δtp3q3E3+λ1mμ-λ3γ+2λ3γω/(μu)-q3E3．

考慮最優(yōu)平衡解, 由系統(tǒng)(2)得E2=βu/[q2(a+bu+cv)]-d/q2,E3=γ[1-ω/(μu)]/q3．

假設(shè)系統(tǒng)在Ei=0或Ei=Emax時(shí)不會(huì)出現(xiàn)平衡點(diǎn), 則根據(jù)Pantryagin最大值原理, 最優(yōu)捕獲控制量Ei(i=2,3)在極值點(diǎn)達(dá)到, 此時(shí)?H/?E=0, 則有λ2=-e-δt[p2-c2/(q2v)],λ3=-e-δt[p3-c3/(q3ω)]．伴隨方程可寫(xiě)成dλ1/dt-A1λ1=-B1e-δt, 得λ1=B1e-δt/(A1+δ), 其中

A1=ru*/k-bαu*v*/(a+bu*+cv*)2,

B1=βv*(a+cv*)[p2-c2/(q2v*)]/(a+bu*+cv*)+
γω*2[p3-c3/(q3ω*)][bαu*v*/(a+bu*+cv*)2-ru*/k]/(μu*2)．

同理,可得λ2=B2e-δt/(A2+δ),λ3=B3e-δt/(A3+δ), 其中A2=cβu*v*/(a+bu*+cv*)2,B2=p2[βu*/(a+bu*+cv*)-d]-B1αu*(a+bu*)/[(A1+δ)(a+bu*+cv*)2];A3=γω*/μu*,B3=p3γ[1-ω*/(μu*)]-B1mu*/(A1+δ)．于是得到奇異軌線為p2-c2/(q2v*)=B2/(A2+δ),p3-c3/(q3ω*)=B3/(A3+δ)．

令f(v*)=p2-c2/(q2v*)-B2/(A2+δ), 顯然limv*→0+f(v*)<0．若存在k>0使f(k)>0, 并且對(duì)任意v*>0, 顯然有f′(v*)>0．則根據(jù)零點(diǎn)定理, 在區(qū)間(0,k)內(nèi), 方程f(v*)=0存在唯一的正根v*=vδ．同理可得ω*=ωδ．再由式(2)中r(1-u*/k)-αv*/(a+bu*+cv*)-mω*=0, 可得唯一u*=uδ, 故E2δ=βuδ/[q2(a+buδ+cvδ)]-d/q2,E3δ=γ[1-ωδ/(μuδ)]/q3．