999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

周期性 Gompertz差分模型的最優脈沖收獲策略

2011-01-01 00:00:00劉彥平,王萬雄,雒志學
經濟數學 2011年2期

摘 要 考慮一個具有周期性脈沖收獲的Gompertz 差分系統. 推導了保證種群系統持續生存、絕滅以及存在全局吸引的正脈沖周期解的充要條件. 以一個周期內持續產量最大化為管理目標 ,通過利用離散的 Pontryagin 最大值原理獲得了最優的脈沖收獲策略,推廣了現有的結論.

關鍵詞 Gompertz 差分模型; 周期解;全局吸引;離散 Pontryagin 最大值原理;脈沖收獲

中圖分類號 O175.1 文獻標識碼 A

Optimal Impulsive Harvesting Policy for a Periodic Gompertz Difference Model

LIU Yan-ping1,WANG Wan-xiong1,LUO Zhi-xue2

(1.College of Science, Gansu Agricultural University, Lanzhou 730070,China;

2.Department of Mathematics, Lanzhou Jiaotong University, Lanzhou 710069,China)

Abstract A Gompertz difference system with periodic impulsive harvesting was investigated.The sufficient and necessary conditions, which guarantee the permanence,extinction and existence of a globally attractive impulsive periodic solution, wereobtained. Choosing maximum sustainable annual yield as the management objective,the optimal impulsive harvesting policywas derived via the discrete maximum principle.The theorems generalize many existing conclusions.

Keywords Gompertz difference model;periodic solution; global attraction ; discrete Pontryagin maximum principle;impulsive harvest

1 引 言

可再生資源的最優管理一直是饒有興趣的研究熱點[1-7].實踐中,資源的管理是一個多目標決策問題.盡管對種群動態行為的研究相當完善,但面臨各種多目標決策問題時,對于如何設計科學的最優收獲方案仍然困難重重.目前,對可再生資源最優化管理的研究,主要集中在連續進行開發的方面,其理論成果相當豐富[1-3].

周期性 Logistic程描述的資源的最優收獲模型,不論是連續收獲[1-2],還是脈沖收獲[3-5],都已經得到廣泛研究.Gompert方程(t)=rxln K/x作為刻畫單種群發展的重要模型之一,也得到充分重視.就連續型 Gompertz統的最優化問題而言,無論是自治系統,還是非自治形式,都得到充分的研究[1,6-7].

現實中,很多資源往往在一定時期內的某些固定時刻脈沖地收獲,以漁業資源的捕獲為例,漁民不可能在全天或24小時的時間段內時時刻刻不停地捕撈.顯然,他們只是在其中特定的時間段或時刻點捕魚,并且以一網一網的方式進行.此外,漁類種群的生長期、成熟期等季節效應也決定著確定的捕魚期.顯而易見,這些捕獲方式都將致使漁類種群遭受瞬時的沖動擾動.因此,資源種群以多次脈沖的方式進行收獲顯得意義重大.人們在生產活動可以通過脈沖擾動的方法來建模.實際上,脈沖微分方程可對這類系統提供一種自然而又合理的刻畫.因此,脈沖微分方程理論已經被應用到諸如種群生態學、疾病的化療治療,種群動力學等方面.2007年,王等[7]研究了關于Gompertz脈沖收獲系統的最優化問題:

(t)=r(t)x(t)ln K(t)x(t),t≠τk, n∈N;

x(t+)=(1-Ek)x(t),t=τk,k,τk∈N.(1)

通過利用比較原理,證明了系統存在唯一全局吸引的脈沖周期解的條件,并由此利用最大值原理獲得了使管理目標最優的最優收獲努力度及相應的最優種群規模.很多學者認為,與連續時間的系統相比,由差分方程控制的離散時間模型對描述世代不重疊的生物種群更加合理[4].此外,被開發的資源種群往往受環境的季節性影響,它們具有某種周期性波動的特征.對王麗敏等研究的模型作一個合乎邏輯的改動,就是借助具有分段常數變元的微分方程將其進行離散化,建立一個相應的差分方程模型.假定系統(1)中種群的平均數量增長率以間隔相等的時間段變化,且對種群數量的測量也以間隔相等的時間段進行,為了模擬資源種群的動態變化規律,需要將脈沖微分系統(1)修改為:

(t)/x(t)=r([t])ln K(t)x(t),

其中,t≠0, 1, 2, ..., [t]表示t的整數部分,t∈(0,+

).設t∈[n,n+1), 則上述方程在每個區間[n,n+1), (n=0, 1, 2, ...)上成立.對該方程在 [n, t)上積分,并令t→n+1, 得

x(n+1)=x(n)exp r(n)ln K(n)x(n).

此式即為連續模型 (1)的離散形式.

本文的目的是研究具有脈沖收獲的周期性Gompertz差分系統的最優化問題.

x(n+1)=x(n)exp r(n)ln K(n)x(n), t≠τk;

x(t+)=(1-Ek)x(t),t=τk,k,τk∈N,(2)

其中,r(n)>0, K(n)>0,Ek, x(τk)表示種群在 t=τk時刻收獲量,Ek (0≤Ek<1)代表脈沖收獲努力度. 假定系統(2)是T周期的,即存在一個正整數T使得r(n+T)=r(n), K(n+T)=K(n).并且假設在一個周期T之內的時刻t=τk (k=1,2,...,q)進行q次收獲,且滿足

0

τk+q=τk+T 以及 Ek+q=Ek (k, n∈N).

2持續生存與絕滅

可再生種群資源的開發直接關系到它們的可持續發展水平,而資源的可持續發展意味著它們可被永久地利用.為了能在較高的生育水平與豐厚經濟利潤基礎之上,實現資源持續發展,首先應該保證它們能夠持續生存.很有必要先討論系統(2)中資源種群的持續生存與絕滅條件.

對系統(2)而言,注意到如果初始值x(0)>1,則對任意n∈N,總有x(n)>1.因而,令u(n)=ln x(n), 則系統(2)可被等價的轉化為

u(n+1)=(1-r(n))u(n)+r(n)ln K(n),t≠τk,

u(t+)=u(t)+ln (1-Ek), 

t=τk,  k,n,τk∈N.(3)

顯然,系統(2)中種群持續生存當且僅當系統(3)的解具有正的上,下界. 而關于系統的正周期解(2), 可給出結論:

定理1 假設(H):∣∏T-1i=0(1-r(i))∣<1成立.且∏T-1i=0(1-r(i))W(0,T)>0,則系統(2)滿足x(0)>1的一切解x(n)都漸近收斂于正周期xP(n);否則,若∏T-1i=0(1-r(i))W(0,T)<0,則對系統(2)具有初值x(0)>1的所有解x(n),n→

時,有 x(n)→0,其中

W(0,T)=∑T-1i=0∏ij=0(1-r(j))-1r(i)ln K(i)+

∑0≤τk

ln xP(n)=∏n-1i=0(1-r(i))[ln x*+

∑n-1i=0∏ij=0(1-r(j))-1r(i)ln K(i)+

∑0≤τk

證明 根據線性差分方程理論,對任意x(0)>1,無脈沖收獲的差分系統(3)可以表示為

u(n)=∏n-1i=0(1-r(i))[u(0)+

∑n-1i=0(∏ij=0(1-r(j)))-1r(i)ln K(i)]. (4)

考慮脈沖時,則對任意n∈N(μT≤n<(μ+1)T),可分別取τk, T+τk, ...,μT+τq (k=1,2,...,q)作為新的初始條件.于是由式(4)得

u(τ1)=∏τ1-1i=0(1-r(i))[u(0)+

∑τ1-1i=0∏ij=0(1-r(j))-1r(i)ln K(i)],

u(τ2)=∏τ2-1i=τ1(1-r(i))[u(τ1)+

∑τ2-1i=τ1∏ij=τ1(1-r(j))-1#8226;

r(i)ln K(i)+ln (1-E1)],

u(τ3)=∏τ3-1i=τ2(1-r(i))[u(τ2)+

∑τ3-1i=τ2∏ij=τ2(1-r(j))-1#8226;r(i)ln K(i)+

ln (1-E2)],

……

u(n)=∏n-1i=μT(1-r(i))[u(μT)+

∑n-1i=μT∏ij=μT(1-r(j))-1r(i)ln K(i)].

將它們加起來,經直接計算可得

u(n)=∏n-1i=0(1-r(i))[u(0)+

∑n-1i=0∏ij=0(1-r(j))-1r(i)ln K(i)+

∑0≤τk

同時,由u(n)的周期性可知u(0)=u(T), 即

u(0)=[1-∏T-1i=0(1-r(i))]-1∏T-1i=0(1-r(i))W(0,T). (6)

結合式(5)~(6),則脈沖系統(2)周期解xP(n)表示為

ln xP(n)=∏n-1i=0(1-r(i))[ln x*+

∑n-1i=0∏ij=0(1-r(j))-1×r(i)ln K(i)+

∑0≤τk

其中,

ln x*=u(0).不難驗證式(7)所表示的周期解滿足

∑T-1n=0r(n)ln K(n)+∑qk=1ln (1-Ek)

=∑T-1n=0r(n)ln x(n).

假定W(0,T)>0,則對系統(3)滿足u(0)>0的解為u(n),都有

∣u(n)-ln xP(n)∣

=∏n-1i=0(1-r(i))∣u(0)-ln  x*∣.

對n≥1,不失一般性,可假設存在一個整數p使得n∈(pT,(p+1)T], 由此推出

∏n-1i=0(1-r(i))=[∏T-1i=0(1-r(i))]p∏n-1i=pT(1-r(i)).

若∏T-1i=0(1-r(i))W(0,T)>0,

則由(H)知n→

時,有∏n-1i=0(1-r(i))→0.從而,系統(1)具有初值u(0)>0的一切解u(n)都漸近收斂于周期解xP(n).同時,若

∏T-1i=0(1-r(i))W(0,T)<0,則由公式只含有限多項知, 一定存在足夠小的整數δ>0使得

∏T-1i=0(1-r(i))W(0,T)<-δ.(8)

為證資源種群在該條件下滅絕:lim n→

x(n)=0.需證lim n→

u(n)=-

.令

M=sup 1≤n≤T∣∏n-1i=0(1-r(i))∣,

N=sup 1≤n≤TM[∑n-1i=0∣∏ij=0(1-r(j))-1r(i)ln K(i)∣

+∑0≤τk

則對任意n∈(μT,(μ+1)T],由式(5)得

u(n)=∏μT-1i=0(1-r(i))∏n-1i=μT(1-r(i))u(0)+

∏μT-1i=0(1-r(i))[∑μT-1i=0∏ij=0(1-r(j))-1r(i)ln K(i)+

∑0

τk<μT-1∏τk-1j=0(1-r(j))-1ln (1-Ek)]+

∏n-1i=μT(1-r(i))[∑n-1i=μT∏ij=0(1-r(j))-1r(i)ln K(i)+

∑μT≤τk

∣∏T-1i=0(1-r(i))∣μMu(0)-μδ+N. (9)

由式(8)知lim n→

u(n)=-

.因此,對所有x(0)>1,有lim n→

x(n)=0.因為u(n)=ln x(n),故結論顯然.

令δ=min n∈[0,T]xp(n),γ=max n∈[0,T]xp(n), 則δ, γ均為正數,且δ<γ.因此,由定理1 知,對充分小的ε0>0,必定存在>0使得由式(5)所表示的任意解滿足:δ-ε0≤x(n)≤γ+ε0, n>. 即: 如下關于系統種群的持續生存結論成立.

推論1 假設(H):∣∏T-1i=0(1-r(i))∣<1成立,并且∏T-1i=0(1-r(i))W(0,T)>0,則系統種群能夠持續生存.即,在一個周期T內進行q次脈沖收獲的情形下,資源種群能夠得以可持續性發展.

由于實際生產活動周期與環境變動的周期間往往存在差異,一種更為現實且自然的考慮是討論比T周期系統(2)與式(3)更一般的情形.即假定生態環境具有周期T1,而人們生產活動的周期是T2,且滿足T1≠T2.于是有

r(n+T1)=r(n), K(n+T1)=K(n) (n∈N);

τk+q=τk+T2, Ek+q=Ek (k∈N).(10)

顯而易見,此時系統(2)的動態特征取決于兩類周期T1與T2之間的具體聯系.不妨設p1T1=p2T2=T0,其中p1, p2為互質的正整數,則上述周期性Gompertz脈沖差分系統的最優化問題轉化為在一個周期T0內進行p2×q次脈沖收獲的情形.對應于定理1,有結論

定理2 若(H):∣∏T0-1i=0(1-r(i))∣<1成立.且∏T0-1i=0(1-r(i))(0,T0)>0,則系統(2)滿足條件(10)的一切解都漸近收斂于脈沖周期解xP(n);否則,若∏T0-1i=0r(i))(0,T0)<0,則對于系統(2)滿足x(0)>1的任意解x(n),當n→

時,x(n)→0,其中

(0,T0)=p1∑T1-1i=0∏ij=0(1-r(j))-1r(i)ln K(i)+

∑0≤τk

證明注意到p1T1=p2T2=T0.由此可知

W(0,T0)=∏T0-1i=0(1-r(i))[∑T0-1i=0∏ij=0(1-

r(j))-1r(i)ln K(i)+∑0≤τk

r(j))-1ln (1-Ek)]

=∏T0-1i=0(1-r(i))[∑p1T1-1i=0∏ij=0(1-r(j))-1r(i)ln K(i)

∑0≤τk

=(0,T0).

從而,可將滿足條件(10)的非自治系統(2)視作T0周期系統,由定理1知該結論成立.

3 最優脈沖收獲策略

這部分的主要目的是推導系統(2)的最優脈沖收獲策略,這種收獲方案可使資源種群在一個周期內的持續產量最大化.為此,定義容許集合:

D={Ek≥0:∣∏T-1i=0(1-r(i))∣<1,且

∏T-1i=0(1-r(i))W(0,T)>0}以及容許控制集

S={Ek∈D:Ek+q=Ek,0≤Ek<1, 

k=0,1,2,…}.

假設存在一個脈沖序列τk∈N(k=0,1,2,…),使得系統這些脈沖時刻遭受一定擾動致使資源種群數量瞬間降低,則對任意n>m0,由于系統(3)有初值u(m0)>0的解u(n)可以表示為

u(n)=∏n-1i=m0(1-r(i))[u(m0)+

∑n-1i=m0∏ij=m0(1-r(j))-1×r(i)ln K(i)+

∑m0≤τk

由式(11)知,對任意n∈(τk,τk+1],有

u(n)=∏n-1i=τk(1-r(i))[u(τk)+

∑n-1i=τk∏ij=τk(1-r(j))-1r(i)ln K(i)]+

∑τk≤τs

這意味著

u(τk+1)=∏τk+1-1i=τk(1-r(i))[u(τk)+

∑τk+1-1i=τk∏ij=τk(1-r(j))-1r(i)ln K(i)+

ln (1-Ek)].(12)

記D(τk+1)=∏τk+1-1i=τk(1-r(i)),且

H(τk+1)=∑τk+1-1i=τk∏ij=τk(1-r(j))-1r(i)ln K(i).

則式(12)可化為

u(τk+1)=D(τk+1)u(τk)+D(τk+1)H(τk+1)+

D(τk+1)ln (1-Ek). (13)

Maximize YEkqk=1=∑qk=1Ekxn(τk), (14)

Subject to:

ln xn(τk+1)=D(τk+1)ln [Ek)xn(τk)]+

D(τk+1)H(τk+1),

Ek∈ S (k=1,2,…,q).

為了獲得更一般最優化問題的最優脈沖收獲策略,以便將其應用到其他類似的問題,在這里不妨先討論其狀態變量的約束條件可化為一般形式的線性差分方程的最優化問題.即

Maximize YEkqk=1=∑qk=1Ekxn(τk),(15)

Subject to:

G(xn(τk+1))=D(τk+1)G((1-Ek)xn(τk))+

D(τk+1)H(τk+1),Ek∈S(k=1,2,…,q),

其中,G(#8226;)是嚴格單調的可微函數, 即G(#8226;)存在反函數,用G-1(#8226;)表示.為了方便, 后文通過xn(τk)來表示資源種群在n+τk刻的種群規模.于是由離Pontryagin最大值原理知.

定理 3 設E*kqk=1,{x*n(τk)}qk=1分別為最優脈沖收獲努力度及其相應最優種群規模. 若

1x*n(τk)G-1(1D(τk+1)[G(xn(τk+1))-

D(τk+1)H(τk+1)])<1.

則最優脈沖收獲努力度為

E*k=1-1x*n(τk)G-1(1D(τk+1)[G(x*n(τk+1))-

D(τk+1)H(τk+1)]).(16)

相應最優種群規模滿足:

G′(x*n(τk+1))=D(τk+1)G′((1-E*k)x*n(τk)).(17)

最大化年持續產量為

YE*kqk=1=∑qk=1[x*n(τk)-

G-1(1D(τk+1)[G(x*n(τk+1))-

D(τk+1)H(τk+1)])]. (18)

證明 構造如下Hamilton函數

H(xn(τk), Ek,λk, τk)=Ekxn(τk)+

λk+1[G-1(D(τk+1)G((1-Ek)xn(τk))+

D(τk+1)H(τk+1))-xn(τk+1)],(19)

其中,λk表示伴隨變量.假設{x*n(τk)∣k=1,2,...,q}是對應于上述控制變量(E*1,E*2,…,E*q)的種群規模.由離散 Pontryagin最大值原理知, 使得問題最優的必要性條件為

Δλk=-Hxn(τk),HEk=0.

即:

Δλk=-Ek-λk+1#8226;

[D(τk+1)(1-Ek)G′(Z)G′((1-Ek)xn(τk))-1],

xn(τk)-λk+1D(τk+1)xn(τk)G′(Z)G′((1-Ek)xn(τk))=0,(20)

其中

Z=G-1(D(τk+1)G((1-Ek)xn(τk))+

D(τk+1)H(τk+1))

即:Z=xn(τk+1).由上述λk=1, (k=1,2,…,q)與式(20)中第二式得

G′(xn(τk+1)=D(τk+1)G′((1-Ek)xn(τk)). (21)

因此,結合式(15)可求出最優種群規模{x*n(τk)}qk=1.再由式(15)式得最優脈沖收獲努力度 E*kqk=1,由式(17)給出.最后,{x*n(τk)}qk=1與E*kqk=1的結果將推出表達式(18).證畢.

若選取G(x)=ln x, (x>0),則系統(15)退化為式(14),從而上述式(21)變為

(1-Ek)xn(τk)xn(τk+1)=D(τk+1)(k=1,2,…,q).

結合等式

ln(xn(τk+1))=D(τk+1)ln ((1-Ek)xn(τk))+D(τk+1)H(τk+1)

便可求得相應的最優種群規模

x*n(τk)=[D(τk)eH(τk)]D(τk)1-D(τk).(22)

因此,作為前面定理3的直接應用,有結論:

定理 4設E*kqk=1,{x*n(τk)}qk=1分別為最優脈沖收獲努力度及其相應最優種群規模.若

D(τk+1)11-D(τk)D(τk)-D(τk)1-D(τk)×

exp {[D(τk+1)H(τk+1)1-D(τk+1)-D(τk)H(τk)1-D(τk)]}<1,

則最優脈沖收獲努力度為

E*k=1-D(τk+1)11-D(τk)D(τk)-D(τk)1-D(τk)×

exp {[D(τk+1)H(τk+1)1-D(τk+1)-D(τk)H(τk)1-D(τk)]}.

相應最優種群規模由關系式(22)給定.而最大化的年持續產量可表示為

YE*kqk=1=∑qk=1exp {D(τk)H(τk)1-D(τk)}#8226;

[D(τk)D(τk)1-D(τk)-D(τk)11-D(τk)].

注1若D(τk+1)=exp {-∫τk+1τkr(s)ds}, 且

H(τk+1)=∫τk+1τkexp {-∫τk+1sr(s)ln K(s)ds},則定理4 恰好為文獻[6]與[7]的主要結論.所以,本文的結果推廣了文獻[6-7]的結論.

若令G(x)=1/xθ (x>0),并考慮到具周期系數的Gilpin-Ayala 脈沖系統具有如下脈沖解

1xθ(τk+1)=D(τk+1)((1-Ek)x(τk))θ+D(τk+1)H(τk+1),(23)

其中,H(τk+1)=∫τk+1τkexp {θ∫sτkr(ι)dι}θr(s)Kθ(s)ds,

D(τk+1)=exp {-θ∫τk+1τkr(s)ds},

則周期性Gilpin-Ayala 脈沖系統的脈沖收獲問題

Maximize YEkqk=1=∑qk=1Ekx(τk),(24)

Subject to:

(t)=r(t)x(t)[1-(x(t)K(t))θ],t≠τk,

x(t+)=(1-Ek)x(t),t=τk

可化為一般性最優化系統(15)的形式,其中r(t), K(t)正的T周期函數,θ>0為常數.于是定理3 中的最優化條件(21)可化為

(1-Ek)xn(τk)xn(τk+1)=D1/(1+θ)(τk+1)(k=1,2,…,q).

結合式(23), 便可求得最優種群規模

x*n(τk)=[1-D1/(1+θ)(τk)D(τk)H(τk)]1θ.(25)

于是,作為前面定理3的直接應用, 有

定理5 設E*kqk=1與{x*n(τk)}qk=1分別為最優化問題(24)的最優脈沖收獲努力度及其最優種群規模. 若

D(τk)Dθ/(1+θ)(τk+1)1-D1/(1+θ)(τk+1)1-D1/(1+θ)(τk)H(τk)H(τk+1)<1,

則最優脈沖收獲努力度為

E*k=1-D(τk)Dθ/(1+θ)(τk+1)1-D1/(1+θ)(τk+1)1-D1/(1+θ)(τk)H(τk)H(τk+1).

相應最優種群規模由關系式(25)給定.

最大化的年持續產量可表示為

YE*kqk=1=∑qk=1{(1-D1/(1+θ)(τk)D(τk)H(τk))1θ-

(1-D1/(1+θ)(τk)D(τk)H(τk))1-θθ1-D1/(1+θ)(τk+1)Dθ/(1+θ)(τk+1)H(τk+1)}.

注2 若θ=1,則該脈沖系統恰好是具有脈沖收獲的周期性logistic系統的脈沖收獲問題,并且定理5為文獻[4-6]的主要結果.

注3 根據線性微分方程的疊加原理,并利用完全類似的方法,也可獲得具有脈沖收獲與臨界退償效應的周期性logistic系統的優化問題:

MaximizeYEkqk=1=∑qk=1Ekx(τk),

Subject to:

(t)=r(t)x(t)(x(t)K0(t)-1)(1-x(t)K(t)),

t≠τk,

x(t+)=(1-Ek)x(t),t=τk.

注4 利用L'Hospital求導法則, 還可以證明:

在上述諸最優化問題中,隨著脈沖收獲間隔趨于零,脈沖收獲方式之下的的最優種群規模將無限趨近連續收獲方式下的的最優種群規模,從而二者間相應的最大化年持續產量YE*kqk=1也有類似的聯系.本文將它解釋為最優脈沖收獲策略為理論最優化結果的實際應用,而連續方式的收獲策略為脈沖收獲策略的理想狀態.

參考文獻

[1] C CLARK. Mathematical Bioeconomics: The Optimal Manage-ment of Renewable Resources[M].New York:Wiley, 1992.

[2] MENG FAN, KE WANG. Optimal harvesting policy for single population with periodic coefficients[J]. Math Biosci,1998,152(2):165-177.

[3] Xiaoying ZHANG,Zhisheng SHUAI, Ke WANG.Optimal harvesting policy for single population[J].Nonlinear Anal(RWA),2003,4(4): 639-651.

[4] Sanyi TANG , R CHEKE, Yanni XIAO .Optimal impulsive harvesting on nonautonomous Beverton-Holt difference equations[J]. Nonlinear Analysis,2006,65(12): 2311-2341.

[5] Yanni XIAO , Daizhan ChENG , Huasu Qin. Optimal impulsive control in periodic ecosystem [J]. Sysems and Control Letters,2006,55(7):558-565.

[6] Lingzhen DONG,LansunChEN , LihuaSun. Optimal harvesting policies for periodic Gompertz systems[J].Nonlinear Analysis (RWA),2007,8(2):572-578.

[7] Limin WANG , Yanshun TAN.Optimal impulsive control policy in periodic Gompertz ecosystem(in chinese)[J]. J. Sys Sci and Math Scis,2007,27(4): 520-528.

注:本文中所涉及到的圖表、注解、公式等內容請以PDF格式閱讀原文

主站蜘蛛池模板: 伊人色综合久久天天| 永久免费av网站可以直接看的| 日韩在线成年视频人网站观看| 日本国产精品一区久久久| 日韩美女福利视频| 无码一区18禁| 欧美视频在线第一页| 免费人成网站在线高清| 欧美第一页在线| 久久无码av一区二区三区| 丰满少妇αⅴ无码区| 欧美成人a∨视频免费观看| 香蕉视频国产精品人| 国产精彩视频在线观看| 波多野结衣一二三| 欧美精品v| 伊人欧美在线| 无码精品福利一区二区三区| 欧美日韩免费观看| 996免费视频国产在线播放| 91精品综合| 国产精品久线在线观看| 中文字幕在线观看日本| 亚洲男女在线| 亚洲综合国产一区二区三区| a国产精品| 亚洲综合极品香蕉久久网| 中文字幕无码电影| 国产精品久久久精品三级| 国产午夜精品一区二区三区软件| 久久精品娱乐亚洲领先| 亚洲无码在线午夜电影| 日韩在线1| 天天躁夜夜躁狠狠躁躁88| 欧美激情首页| 国产精品yjizz视频网一二区| 天天婬欲婬香婬色婬视频播放| 毛片手机在线看| 福利片91| 国产福利免费观看| hezyo加勒比一区二区三区| 国产高清无码麻豆精品| 欧美亚洲综合免费精品高清在线观看| 免费可以看的无遮挡av无码 | 精品国产免费观看一区| 日韩在线永久免费播放| 偷拍久久网| 欧美精品v欧洲精品| 人妻中文字幕无码久久一区| 97视频在线观看免费视频| 国产精品男人的天堂| jizz亚洲高清在线观看| 露脸真实国语乱在线观看| 最新亚洲av女人的天堂| 亚洲经典在线中文字幕| 亚洲无码高清一区| 69av在线| 成人噜噜噜视频在线观看| 91欧美亚洲国产五月天| 久久久久久久久久国产精品| 亚洲日本在线免费观看| 国产在线视频福利资源站| 国产尤物视频在线| 久久99国产视频| 日本一区中文字幕最新在线| 色欲色欲久久综合网| 久久国产乱子伦视频无卡顿| 无码精品国产VA在线观看DVD | 亚洲乱码精品久久久久..| 激情综合网址| 亚洲午夜片| 中文字幕天无码久久精品视频免费| 国产区网址| 日本亚洲国产一区二区三区| 青青青国产视频| 欧美一级99在线观看国产| 午夜精品一区二区蜜桃| 国产无码精品在线| 天天躁日日躁狠狠躁中文字幕| 久久99国产综合精品1| AⅤ色综合久久天堂AV色综合 | 欧美成人综合视频|