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

基于動態(tài)事件觸發(fā)機制的網(wǎng)絡(luò)化系統(tǒng)有限頻域故障檢測

2022-12-31 00:00:00朱淇姜順潘豐
計算機應(yīng)用研究 2022年9期

收稿日期:2022-02-28;修回日期:2022-04-08" 基金項目:國家自然科學(xué)基金資助項目(61403168,61876073)

作者簡介:朱淇(1997-),男,江蘇常州人,碩士,主要研究方向為網(wǎng)絡(luò)化系統(tǒng)的故障診斷;姜順(1981-),男(通信作者),河南信陽人,副教授,博士,主要研究方向為網(wǎng)絡(luò)化控制系統(tǒng)的分析與綜合(haveshun@jiangnan.edu.cn);潘豐(1963-),男,江蘇常熟人,教授,博士,主要研究方向為生產(chǎn)過程建模、優(yōu)化與控制.

摘 要:針對通信帶寬受限的網(wǎng)絡(luò)環(huán)境,引入一種基于動態(tài)事件觸發(fā)機制的數(shù)據(jù)傳輸策略,研究了一類非線性網(wǎng)絡(luò)化系統(tǒng)在隨機網(wǎng)絡(luò)攻擊下的有限頻域故障檢測問題。首先,在考慮故障靈敏性和擾動魯棒性的前提下,利用狀態(tài)增廣的方法將原系統(tǒng)的故障檢測問題轉(zhuǎn)換成H-/H∞濾波問題;然后,在考慮扇區(qū)有界非線性和隨機網(wǎng)絡(luò)攻擊的情況下,將故障的有限頻域特性考慮到H-性能指標(biāo)的設(shè)計中,并結(jié)合有限頻輸入特性,給出有限頻故障輸入下的故障檢測濾波器與動態(tài)事件觸發(fā)機制的聯(lián)合設(shè)計算法;最后,通過攪拌釜式反應(yīng)器系統(tǒng)的仿真算例驗證了該方法的有效性。

關(guān)鍵詞:故障檢測;網(wǎng)絡(luò)化系統(tǒng);動態(tài)事件觸發(fā)機制;有限頻域;網(wǎng)絡(luò)攻擊

中圖分類號:TP273"" 文獻標(biāo)志碼:A

文章編號:1001-3695(2022)09-030-2757-06

doi:10.19734/j.issn.1001-3695.2022.02.0081

Dynamic event-triggered fault detection for networked systems in finite-frequency domain

Zhu Qi1,2,Jiang Shun1,2,Pan Feng1,2

(1.College of Internet of Things Engineering,Jiangnan University,Wuxi Jiangsu 214122,China;2.Key Laboratory of Advanced Process Control of Light Industry for Ministry of Education,Wuxi Jiangsu 214122,China)

Abstract:This paper investigated the fault detection problem in finite-frequency domain for nonlinear networked system under stochastic cyber-attack.To save limited network resource,this paper introduced a novel dynamic event-triggered scheme.Firstly,under the consideration of fault sensitivity and disturbance robustness,this paper converted the addressed fault detection problem into an auxiliary H-/H∞ filtering problem by augmenting the states of the original system and the fault detection filter.Taking sector bounded nonlinearity and stochastic cyber-attacks into consideration,the design of H- performance index inclu-ded the frequency characteristics of fault signals.Combined with finite-frequency input characteristics,this paper proposed the joint design algorithm for fault detection filter and dynamic event-triggered scheme under the finite-frequency fault input.Finally,a simulation example of stirred tank reactor system verifies the effectiveness of the proposed method.

Key words:fault detection;networked system;dynamic event-triggered scheme;finite-frequency domain;cyber-attack

0 引言

網(wǎng)絡(luò)化系統(tǒng)是指被控對象和傳感器、控制器、執(zhí)行器以及其他系統(tǒng)組件之間通過共享通信網(wǎng)絡(luò)連接而形成的復(fù)雜控制系統(tǒng)。隨著通信技術(shù)和計算機網(wǎng)絡(luò)技術(shù)的高速發(fā)展,控制系統(tǒng)也向網(wǎng)絡(luò)化與智能化方向轉(zhuǎn)變,網(wǎng)絡(luò)化系統(tǒng)逐漸被廣泛地應(yīng)用于實際工業(yè)領(lǐng)域中,如移動無線傳感器網(wǎng)絡(luò)[1]、無人機[2]、遠程醫(yī)療[3]等。然而,共享通信網(wǎng)絡(luò)為網(wǎng)絡(luò)化系統(tǒng)帶來諸多便利的同時,也帶來了一系列挑戰(zhàn),如網(wǎng)絡(luò)時延[4]、數(shù)據(jù)丟包[5]、通信約束[6]等問題,甚至引起故障。近年來,隨著人們對控制系統(tǒng)可靠性與安全性要求的不斷提高,針對非理想環(huán)境下網(wǎng)絡(luò)化系統(tǒng)的故障檢測問題已成為控制領(lǐng)域的熱點研究問題[7,8]。

在另一個研究前沿,為了有效地利用系統(tǒng)中有限的通信資源,事件觸發(fā)機制應(yīng)運而生[9,10]。其主要策略是通過給出一個預(yù)先設(shè)定的事件觸發(fā)條件來判斷是否將測量輸出傳輸至濾波器,這是一種按需執(zhí)行的非等周期觸發(fā)方式。文獻[11,12] 在故障檢測方法中引入了靜態(tài)事件觸發(fā)機制,當(dāng)最新觸發(fā)時刻的測量值與當(dāng)前測量值之間的差值超出給定閾值時,更新數(shù)據(jù)傳輸?shù)挠|發(fā)時刻。然而,上述靜態(tài)事件觸發(fā)機制的低觸發(fā)參數(shù)靈活度限制了其適用范圍,因此有學(xué)者在此基礎(chǔ)上提出了動態(tài)事件觸發(fā)機制[13~15]。文獻[16] 考慮了一類連續(xù)時間非線性隨機系統(tǒng)的故障檢測問題,通過設(shè)計自適應(yīng)事件觸發(fā)控制率,使得觸發(fā)函數(shù)的閾值隨著測量輸出的變化動態(tài)調(diào)整。文獻[17] 針對離散時間網(wǎng)絡(luò)化系統(tǒng),在設(shè)計故障檢測觀測器的基礎(chǔ)上,構(gòu)建了基于集員估計的殘差評估動態(tài)閾值,進一步降低網(wǎng)絡(luò)資源占用且滿足更為靈活的系統(tǒng)設(shè)計要求。此外,在網(wǎng)絡(luò)信道傳輸數(shù)據(jù)過程中,系統(tǒng)可能會遭受網(wǎng)絡(luò)攻擊。與周期采樣機制相比,事件觸發(fā)機制下按需發(fā)送的數(shù)據(jù)包數(shù)量減少,但對系統(tǒng)性能的影響更大,若這些數(shù)據(jù)包被網(wǎng)絡(luò)攻擊竄改,則會破壞系統(tǒng)的穩(wěn)定性[18]。因此,在針對動態(tài)事件觸發(fā)機制下故障檢測的研究中設(shè)計抵御網(wǎng)絡(luò)攻擊的故障檢測策略至關(guān)重要。

需要注意的是,上述基于事件觸發(fā)機制的故障檢測方法只考慮了全頻域內(nèi)系統(tǒng)的整體性能,然而實際工程中的故障信號往往發(fā)生在有限頻率范圍內(nèi),例如:緩變故障處于低頻域,卡死型故障的頻率為零[9]。為了更好地刻畫有限頻域特性,文獻[19] 提出的廣義 Kalman-Yakubovich-Popov(KYP) 引理提供了一種處理故障信號頻率特性的有效方法。為了將廣義KYP引理的適用范圍推廣至非線性系統(tǒng),文獻[20] 研究了基于靜態(tài)事件觸發(fā)機制的有限頻域故障檢測問題,并將非線性誤差動態(tài)系統(tǒng)轉(zhuǎn)換成一類線性時變參數(shù)系統(tǒng)來處理。然而,目前基于動態(tài)事件觸發(fā)機制的非線性網(wǎng)絡(luò)化系統(tǒng)的有限頻域故障檢測問題尚未引起關(guān)注,更不用說考慮網(wǎng)絡(luò)攻擊的情況,這也是本文研究的主要動機和出發(fā)點。

基于以上討論,本文研究了隨機網(wǎng)絡(luò)攻擊下非線性網(wǎng)絡(luò)化系統(tǒng)的有限頻域故障檢測問題,主要貢獻體現(xiàn)在:a)所設(shè)計的故障檢測濾波器同時考慮H∞擾動魯棒性能和H-故障敏感性能;b)為了進一步節(jié)約網(wǎng)絡(luò)資源,采用一種具有更高靈活度參數(shù)的動態(tài)事件觸發(fā)機制;c)依據(jù)有限頻域不等式的時域解釋[21],直接在時域上設(shè)計有限頻域故障檢測濾波器。

1 問題描述

1.1 系統(tǒng)建模

考慮離散時間非線性網(wǎng)絡(luò)化系統(tǒng)如下:

x(k+1)=Ax(k)+Φ(k,x(k))+Bw(k)+Ff(k)

y(k)=Cx(k)+Dw(k)(1)

其中:x(k)∈Euclid Math TwoRApnx為系統(tǒng)的狀態(tài)向量;w(k)∈Euclid Math TwoRApnw為外部擾動且滿足2[0,∞)范數(shù)有界;f(k)∈Euclid Math TwoRApnf為待檢測的有限頻域故障信號;y(k)∈Euclid Math TwoRApny是測量輸出向量;系統(tǒng)參數(shù)矩陣A、B、F、C和D為已知的適維常數(shù)矩陣且(C,A)可觀測。非線性向量值函數(shù)Φ(k,x(k))滿足如下扇形有界條件:

[Φ(k,x(k))-Φ(k,(k))-S1(x(k)-(k))]T

×

[Φ(k,x(k))-Φ(k,(k))-S2(x(k)-(k))]lt;0(2)

其中:S1,S2∈Euclid Math TwoRApnx×nx為已知常數(shù)矩陣且Φ(k,0)=0。故障信號所處的頻率域可分為低頻段、中頻段與高頻段,統(tǒng)一描述為

Ωθ{θ∈Euclid Math TwoRAp|θ1≤θ≤θ2}(3)

當(dāng)0lt;θ2-θ1≤2π時,式(3)中的集合Ωθ表示中頻段,θ1和θ2分別為中頻段的下界和上界。

當(dāng)θ1=-θ2=-θl時,式(4)中的集合Ωθ表示低頻段,θl為低頻段的上界。

Ωθ{θ∈Euclid Math TwoRAp||θ|≤θl}(4)

當(dāng)θ1=θh,θ2=2π-θh時,式(5)中集合Ωθ表示高頻段,θh為高頻段的下界。

Ωθ{θ∈Euclid Math TwoRAp||θ|≥θh}(5)

定義1 考慮如下離散時間系統(tǒng):

η(k+1)=Aη(k)+Bv(k)(6)

若在輸入信號v(k)∈2[0,∞)作用下,系統(tǒng)狀態(tài)η(k)滿足:

∑∞k=0[η(k+1)-η(k)][η(k+1)-η(k)]T≤(2sinθl2)2∑∞k=0η(k)ηT(k)(7)

或ejθfc∑∞k=0[η(k+1)-ejθf1η(k)][η(k+1)-ejθf2η(k)]T≤0(8)

或∑∞k=0[η(k+1)-η(k)][η(k+1)-η(k)]T≥(2sinθh2)2∑∞k=0η(k)ηT(k)(9)

則稱v(k)∈2[0,∞)分別為對應(yīng)于[-θl,θl]的低頻輸入,對應(yīng)于[θ1,θ2]的中頻輸入或?qū)?yīng)于(-∞,-θh]∪[θh,+∞)的高頻輸入。式(7)~(9)中的θf表示故障信號f(k)的頻率,θf1和θf2分別表示故障頻率范圍的邊界值,θfp=(θf1+θf2)/2,θfc=(θf2-θf1)/2。

注釋1 定義1借鑒了文獻[21]對于有限頻域不等式的時域解釋。輸入信號的有限頻域特性表現(xiàn)為對系統(tǒng)狀態(tài)的驅(qū)動能力,如:頻域處于[-θl,θl]的低頻輸入v(k)∈2[0,∞)滿足式(7),意味著在輸入v(k)的作用下,系統(tǒng)狀態(tài)η(k)的相對變化速率低于2 sin(θl/2)。

1.2 動態(tài)事件觸發(fā)機制的建模

測量輸出y(k)通過帶寬有限的通信網(wǎng)絡(luò)進行傳輸,為了節(jié)省網(wǎng)絡(luò)資源,采用如下動態(tài)事件觸發(fā)機制:

ki+1=infk∈N+{kgt;ki|λ1(k)+σ1yT(k)y(k)-εeTy(k)Ωey(k)≤0}

(k+1)=λ2(k)+σ2yT(k)y(k)-εeTy(k)Ωey(k)(10)

其中:ki表示第i個觸發(fā)時刻且k0=0;(k)是柔性變量且(0)≥0;λ1、λ2、σ1、σ2和ε均為給定的非負實數(shù);y(k)為當(dāng)前測量輸出且ey(k)=(k)-y(k),表示上一觸發(fā)時刻的測量輸出(k)與當(dāng)前測量輸出y(k)之間的差值,(k)=y(ki),k∈[ki,ki+1),i,k∈Euclid Math TwoNAp。只有當(dāng)y(k)滿足觸發(fā)條件式(10)時,測量輸出數(shù)據(jù)才會被傳輸至濾波器的輸入端。圖1給出了動態(tài)事件觸發(fā)機制下采樣和釋放瞬間的時間序列,由此可見,傳感器發(fā)送數(shù)據(jù)是非等周期的。

注釋2 當(dāng)λ1趨于0時,觸發(fā)條件式(10)可改寫為

ki+1=infk∈N+{kgt;ki|σ1yT(k)y(k)-εeTy(k)Ωey(k)≤0}(11)

可知式(11)退化為靜態(tài)事件觸發(fā)機制。相比于文獻[11,12] 考慮的具有固定觸發(fā)閾值的靜態(tài)事件觸發(fā)機制,本文所考慮的動態(tài)事件觸發(fā)機制更具一般性與靈活性。

1.3 故障檢測濾波器的設(shè)計

為了檢測系統(tǒng)中的執(zhí)行器故障,設(shè)計如下形式的濾波器:

(k+1)=A(k)+Φ((k))+L[(k)-C(k)]

r(k)=(k)-C(k)(12)

其中:(k)∈Euclid Math TwoRApnx和r(k)∈Euclid Math TwoRApny為濾波器的狀態(tài)向量和殘差向量;L∈Euclid Math TwoRApnx×ny為待設(shè)計的濾波器增益矩陣。本文考慮網(wǎng)絡(luò)中存在虛假數(shù)據(jù)注入攻擊,其通過竄改數(shù)據(jù)包內(nèi)容而破壞數(shù)據(jù)的真實性和完整性。濾波器端在隨機網(wǎng)絡(luò)攻擊下的輸入信號(k)可以表示為

(k)=[1-α(k)]Cg(x(k))+α(k)(k)(13)

α(k)為服從Bernoulli分布的白噪聲序列。

P{α(k)=1}=E{α(k)}=

P{α(k)=0}=1-E{α(k)}=1-

V{α(k)}=E{(α(k)-)2}=(1-)(14)

其中:∈[0,1]是已知常數(shù);P{·}表示事件發(fā)生概率;E{·}表示數(shù)學(xué)期望;V{·}表示方差。當(dāng)α(k)=0時,表示信號在傳輸中可能受到網(wǎng)絡(luò)攻擊干擾,其概率為1-。當(dāng)α(k)=1時,表示信號在傳輸中沒有遭受網(wǎng)絡(luò)攻擊,系統(tǒng)正常工作的概率為。越小,網(wǎng)絡(luò)攻擊出現(xiàn)的概率越高。由于網(wǎng)絡(luò)攻擊的攻擊效力是非線性的,故假設(shè)g(x(k))是一個非線性函數(shù),且滿足下列二次約束條件:

gT(k,x(k))g(k,x(k))≤κ2xT(k)HTHx(k)(15)

其中:κgt;0為已知非線性函數(shù)g(k,x(k))的有界參數(shù);H為常數(shù)矩陣。定義ξ(k)=[xT(k) eT(k)]T、e(k)=x(k)-(k)和ΔΦ(k)=Φ(k,x(k))-Φ(k,(k)),并結(jié)合式(1)(12)和(13),通過狀態(tài)增廣可以得到如下濾波誤差系統(tǒng):

ξ(k+1)=(1-(k)2)ξ(k)+(1-(k)2)w(k)+

f(k)+(1-(k)2)ey(k)+3Φ(k,x(k))+

4ΔΦ(k)+(1+(k)2)g(k,x(k))

r(k)=(1+(k)2)ξ(k)+(1+(k)2)w(k)+

(5+(k)6)ey(k)+(3-(k)4)g(k,x(k))

y(k)=2ξ(k)+2w(k) (16)

其中:0代表適當(dāng)維數(shù)的零矩陣;I代表適當(dāng)維數(shù)的單位矩陣。其余參數(shù)矩陣為1=A0

(1-)LCA-LC,2=00

LC0,1=B

B-LD,2=0

LD,1=0

-L,2=0

L,3=I0,4=0

I,=FF,1=0(-1)LC,2=0LC,1=(-1)CC,2=C0,1=D,2=D,5=I,6=I,3=(1-)C,4=C,(k)=α(k)-。

本節(jié)的目的是設(shè)計動態(tài)事件觸發(fā)機制與故障檢測濾波器,使得濾波誤差系統(tǒng)式(16)均方漸近穩(wěn)定且滿足H-/H∞性能指標(biāo),即滿足下列要求:a)濾波誤差系統(tǒng)式(16)是均方漸近穩(wěn)定的;b)當(dāng)f(k)=0時,系統(tǒng)滿足如下H∞性能指標(biāo):

∑∞k=0‖r(k)‖2≤γ2∑∞k=0‖w(k)‖2(17)

并使得γ盡可能小,γ越小,擾動信號對殘差越魯棒;c)當(dāng)w(k)=0時,系統(tǒng)滿足如下H-性能指標(biāo):

∑∞k=0‖r(k)‖2≥β2∑∞k=0‖f(k)‖2(18)

并使得β盡可能大,β越大,故障信號對殘差越敏感。

1.4 殘差評估機制

結(jié)合故障檢測濾波器生成的殘差,選取殘差評價函數(shù)J(k)與閾值Jth來檢測故障。J(k)與Jth分別為

J(k)=E{∑ks=0rT(s)r(s)}1/2,Jth=supw∈2,f=0J(T)(19)

其中:T表示有限的評估時間長度。通過比較J(k)與Jth的大小來判斷故障是否發(fā)生:

J(k)gt;Jth故障報警J(k)≤Jth無故障(20)

2 主要結(jié)果

本章以線性矩陣不等式的形式給出了使得濾波誤差系統(tǒng)式(16)滿足H∞性能指標(biāo)γ與H-性能指標(biāo)β的充分條件,并通過求解凸優(yōu)化問題得到濾波器的參數(shù)。下面給出推導(dǎo)過程中將用到的引理。

引理1[22] 對于一個向量值函數(shù)Φ(k,x(k))和已知常數(shù)矩陣不等式S1,S2∈Euclid Math TwoRApn×n,若存在:

[Φ(k,x(k))-S1x(k)]T[Φ(k,x(k))-S2x(k)]lt;0(21)

則有式(22)成立。

[xT(k) ΦT(k,x(k))]12

*Ix(k)

Φ(k,x(k))lt;0(22)

其中:1=(ST1S2+ST2S1)/2;2=-(ST1+ST2)/2。

引理2[23] 對于矩陣A,Q=QT和Pgt;0,矩陣不等式-Q+ATPAlt;0成立當(dāng)且僅當(dāng)存在一個矩陣G,滿足:

-QATGGTAP-G-GTlt;0(23)

2.1 穩(wěn)定性及殘差對擾動的魯棒性條件

定理1 給定滿足0lt;λ1≤λ2、λ1+λ2lt;1、0≤σ1≤σ3和εgt;0的實數(shù)λ1、λ2、σ1、σ2和ε,如果存在正定對稱矩陣Pd,矩陣W1、X和濾波器參數(shù)矩陣L,使得如下矩陣不等式成立:

[Np,q]15×15lt;0(24)

則稱濾波誤差系統(tǒng)式(16)均方漸近穩(wěn)定,且滿足H∞性能指標(biāo)γ。其中,對稱矩陣Np,q中的非零項如下:N1,1=-Pd1+κ2HTH-S^1,N1,2=-Pd2,N1,6=-2,

N1,9=-ATW11-μ1(1-)CTX,N1,10=-ATW12-(1-)CTX,

N1,11=μ1CTX,N1,12=CTX,N1,13=CT,N1,14=-(1-)CT,N1,15=CT,N2,2=-Pd3-1,N2,7=2,N2,9=-μ1ATW11+μ1CTX,N2,10=-ATW+CTX,N2,14=CT,N3,3=-γ2I,N3,9=-BTW11-μ1BTX+μ1DTX,N3,10=-BTW12-BTW+DTX,N3,11=μ1DTX,N3,12=DTX,N3,13=DT,N3,14=DT,N3,15=DT,N4,4=-2εΩ,N4,9=μ1X,N4,10=X,N4,11=μ1X,N4,12=X,N4,14=I,N4,15=I,N5,5=(λ1+λ2-1)I,N6,6=N7,7=N8,8=-I,N6,9=-W11,N6,10=-W12,N7,9=-μ1W,N7,10=-W,N8,9=μ1(1-)CTX,N8,10=(1-)CTX,N8,11=-μ1CTX,N8,12=-CTX,N8,14=(1-)CT,N8,15=-CT,N9,9=N11,11=Pd1-W11-WT11,N9,10=N11,13=Pd2-W12-μ1WT,N10,10=N12,12=Pd3-W-WT,N13,13=-(σ1+σ2)-1I,N14,14=N15,15=-I。

證明 構(gòu)造如下Lyapunov-Krasovskii泛函

Vd(k)=ξT(k)Pdξ(k)+(k)(25)

沿系統(tǒng)式(16)的軌跡求偏差可知:

E{ΔVd(k)}=E{Vd(k+1)}-E{Vd(k)}=

ξT(k+1)Pdξ(k+1)-ξT(k)Pdξ(k)+(k+1)-(k)=

[ξT(k)(T1-(k)T2)+wT(k)(T1-(k)T2)+(k)+

eTy(k)(T1+(k)T2)+gT(x(k))(T1+(k)T2)]×

Pd[(1-(k)2)ξ(k)+(1-(k)2)w(k)+(k)+(1+(k)2)ey(k)+(1+(k)2)g(x(k))]-ξT(k)Pdξ(k)+(λ3-1)(k)+σ3yT(k)y(k)-eTy(k)Ωey(k)(26)

根據(jù)扇形有界條件式(2)和引理1可得:

[Φ(k,x(k))-S1x(k)]T[Φ(k,x(k))-S2x(k)]=

[Φ(k,x(k))-S17ξ(k)]T[Φ(k,x(k))-S27ξ(k)]=ηT(k)Λ1η(k)lt;0(27)

[ΔΦ(k)-S1e(k)]T[ΔΦ(k)-S2e(k)]=[ΔΦ(k)-S18ξ(k)]T[ΔΦ(k)-S28ξ(k)]=ηT(k)Λ2η(k)lt;0(28)

其中:ηT(k)=[ξT(k) wT(k) eTy(k) T(k) ΦT(k,x(k)) ΔΦT(k) gT(x(k))],(k)=1/2(k),S^1=(ST1S2+ST2S1)/2,S^2=-(ST1+ST2)/2。

Λ1=T717000T7200

000000

00000

0000

I00

00

0,

Λ2=T8180000T820

000000

00000

0000

000

I0

0

其中:7=[I 0],8=[0 I]。結(jié)合式(26)~(28),動態(tài)事件觸發(fā)機制式(10)和引理2可推導(dǎo)出:

E{ΔVd(k)}=E{ΔVd(k)}+gT(k,x(k))g(k,x(k))-gT(k,x(k))

g(k,x(k))≤

E{ΔVd(k)}+κ2ξT(k)HT0H0ξ(k)-gT(k,x(k))

g(k,x(k))-

ηT(k)Λ1η(k)-ηT(k)Λ2η(k)-ξT(k)Pdξ(k)+(λ1+

λ2-1)

(k)+(σ1+σ2)yT(k)y(k)-2εeTy(k)Ωey(k)=ηT(k)

[Ξ11-Ξ12Ξ22ΞT12-diag{0,-γ2I,0,0,0,0,0}-Ξ13Ξ33ΞT13]η(k)(29)

當(dāng)w(k)=0時,等式ηT(k)diag{0,-γ2I,0,0,0,0}η(k)=0成立。根據(jù)Schur補引理及式(29),可知Ξ11-Ξ12Ξ22ΞT12-Ξ13Ξ33ΞT13lt;0成立的充分條件為

Ξ11Ξ12Ξ13Ξ14Ξ2200Ξ330Ξ44lt;0(30)

其中:Ξ11=Ξ0000-T72-T820-γ2I00000

-2εΩ0000

(λ1+λ2-1)I000

-I00

-I0

-I

Ξ0=-Pd-T717-T818+κ2HT0H0

Ξ12=T1-T2

T1-T2

T1-T2

00

T30

T40

T1T2,

Ξ13=T2

T2

00000,Ξ14=T1T2T1T2

T5T6000000T3-T4

Ξ22=diag{-P-1d,-P-1d},Ξ33=-(σ1+σ2)-1IΞ44=diag{-I,-I},H0=[H 0],=α-(1-α-)

故有Ξ11-Ξ12Ξ22ΞT12-Ξ13Ξ33ΞT13-ηT(k)diag{0,-γ2I,0,0,0,0,0}η(k)lt;0成立。因此,濾波誤差系統(tǒng)式(16)均方漸近穩(wěn)定。為了證明濾波誤差系統(tǒng)式(16)在零初始條件下滿足期望的H∞性能指標(biāo)γ,考慮如下評估指標(biāo)函數(shù):

J∞=E{∑∞k=0rT(k)r(k)-γ2wT(k)w(k)}≤

E{∑∞k=0rT(k)r(k)-γ2wT(k)w(k)}-E{Vd(0)}+E{Vd(∞)}=

E{∑∞k=0rT(k)r(k)-γ2wT(k)w(k)+ΔVd(k)}=

ηT(k)[Ξ11-Ξ12Ξ22ΞT12-Ξ13Ξ33ΞT13-Ξ14Ξ44ΞT14]η(k)(31)

根據(jù)Schur補引理和式(30),可得Ξ11-Ξ12Ξ22ΞT12-Ξ13Ξ33ΞT13-Ξ14Ξ44ΞT14lt;0,則評估指標(biāo)函數(shù)滿足J∞lt;0,即

∑∞k=0rT(k)r(k)≤γ2∑∞k=0wT(k)w(k)(32)

因此,濾波誤差系統(tǒng)式(16)滿足期望的H-性能指標(biāo)β。結(jié)合Shur補引理和引理2可知,式(30)可改寫為如下矩陣不等式:

Ξ1112Ξ13Ξ14

2200

Ξ330

Ξ44lt;0(33)

其中,12=

WT11WT11WT110WT13WT14WT11-WT12-WT12-WT12000WT12,22=diag{Pd-W1-WT1,Pd-W1-WT1}。將Pd=Pd1Pd2

Pd3,W1=W11W12

μ1WW,X=LTW代入式(33),可得式(24)成立。定理證畢。

2.2 殘差對有限頻域故障的敏感性條件

定理2 給定滿足0lt;λ1≤λ2、λ1+λ2lt;1、0≤σ1≤σ3和εgt;0的實數(shù)λ1、λ2、σ1、σ2和ε,如果存在正定對稱矩陣Pf、Qf,使得如下不等式成立:

110341

-20-20002

I000000T

Ξf110341

-20-20002

I000000+

1050003

206000-4

0I00000T

Πf1050003

206000-4

0I00000+

Γ0000-T72-T820

000000

-2εΩ0000

(λ1+λ2-1)I000

-I00

-I0

-Ilt;0(34)

則系統(tǒng)式(16)滿足有限頻域H-性能指標(biāo)‖r(k)‖2≥β‖f(k)‖2。其中,Γ0=(σ1+σ2)T22+κ2HT0H0-HT21H2,Πf=-I00

-I0β2I,Ξf在不同頻域率下的取值如表1所示。

證明 首先考慮無外部擾動下濾波誤差系統(tǒng)式(16)的中頻段情況。分別對式(35)左乘和右乘矩陣η(k)=[ξT(k)fT(k) eTy(k) ΦT(k,x(k)) ΔΦT(k) gT(x(k)) T(k)]及其轉(zhuǎn)置,并結(jié)合跡運算關(guān)系式uTQv=tr(QvuT),可以推導(dǎo)出:

-ξT(k)Pfξ(k)+ξT(k+1)Pfξ(k+1)-rT(k)r(k)+β2fT(k)f(k)+tr [Qf·(e-jθfpξ(k+1)ξT(k)+ejθfpξ(k)ξT(k+1)-

2 cosθfcξ(k)ξT(k))]+

κ2ξT(k)HT0H0ξ(k)-gT(k,x(k))g(k,x(k))-ηT(k)

(Λ1+Λ2)η(k)+(λ1+λ2-1)(k)+(σ1+σ2)

yT(k)y(k)-2εeTy(k)Ωey(k)≤0(35)

表1 集合θ與矩陣Ξ在不同頻域的取值

Tab.1 θ and Ξ for different frequency ranges

頻域θΞf

低頻|θ|≤θlPf0Qf

Pf0

-Pf-2 cos(θfl)Qf

中頻θ1≤θ≤θ2Pf0ejθfpQf

Pf0

-Pf-2 cos(θfm)Qf

高頻|θ|≥θhPf0-Qf

Pf0

-Pf+2 cos(θfh)Qf

考慮到非線性項的扇形有界條件、網(wǎng)絡(luò)攻擊函數(shù)的二次約束條件以及動態(tài)事件觸發(fā)機制式(2),對于k∈[ki,ki+1),可推導(dǎo)出

ξT(k+1)Pfξ(k+1)-ξT(k)Pfξ(k)+(k+1)-(k)-

rT(k)r(k)+β2fT(k)f(k)+tr [Qf·(e-jθfpξ(k+1)ξT(k)+

ejθfpξ(k)ξT(k+1)-2 cos θfcξ(k)ξT(k))]≤0(36)

選取Lyapunov-Krasovskii泛函為Vf(k)=ξT(k)Pfξ(k)+(k),則有

ΔVf(k)-rT(k)r(k)+β2fT(k)f(k)+tr [Qf·(e-jθfpξ(k+1)·

ξT(k)+ejθfpξ(k)ξT(k+1)-2 cosθfcξ(k)ξT(k))]≤0(37)

在零初始條件下,對式(37)兩邊同時取k從0到∞進行累加可得:

∑∞k=0[-rT(k)r(k)+β2fT(k)f(k)]+tr(QfM)≤0(38)

其中:M:=∑∞k=0[-e-jθfpξ(k+1)ξT(k)-ejθfpξ(k)ξT(k+1)+2 cos θfcξ(k+1)ξT(k+1)]。根據(jù)歐拉公式,-M可以在計算后轉(zhuǎn)換為定義1中不等式(8)的左側(cè)項,故M是半正定的。由于Qfgt;0,當(dāng)式(38)成立時,式(38)左側(cè)的最后一項為非負項。因此,可以得到∑∞k=0rT(k)r(k)≥β2∑∞k=0fT(k)f(k),即定義1中的有限頻域H-性能指標(biāo)β得到滿足。

同理,令θf1=-θf2=-θfl可以得出執(zhí)行器故障發(fā)生在低頻段時的結(jié)果;令θf1=θfh和θf2=2π-θfh可得出執(zhí)行器故障發(fā)生在低頻段時的結(jié)果,其證明過程與上述中頻段類似。定理證畢。

注釋3 式(34)給出了有限頻域H-故障敏感性的設(shè)計條件,是廣義KYP引理在事件觸發(fā)機制與系統(tǒng)存在非線性項情況下的推廣。

2.3 濾波器參數(shù)計算

結(jié)合定理1、2,使用MATLAB中的YALMIP工具箱求解如下凸優(yōu)化問題,可以獲得最優(yōu)濾波器參數(shù)L,以及H-/H∞性能指標(biāo)γmin和βmax:

minPdgt;0,Pfgt;0,Qfgt;0γ-β s.t.式(24)(34)(39)

其中:濾波器增益矩陣可由L=W-TXT計算求得。

注釋4 為了降低事件觸發(fā)機制式(10)的數(shù)據(jù)發(fā)送率觸發(fā)閾值參數(shù),ε應(yīng)取較小值,然而ε取值較小可能導(dǎo)致定理1、2中的約束條件無解。因此,為了在保證故障檢測性能的前提下獲得較低數(shù)據(jù)發(fā)送率,提出濾波器式(12)與動態(tài)事件觸發(fā)機制式(10)的聯(lián)合設(shè)計算法如下:

算法1 有限頻故障輸入下的H-/H∞故障檢測濾波器與動態(tài)事件觸發(fā)機制的聯(lián)合設(shè)計

a)初始化事件觸發(fā)參數(shù)λ1、λ2、σ1、σ2、ε0和步長增量Δε。

b)若定理1、2中的線性矩陣不等式可解,保存ε=ε0、Ω、L,進入步驟c);否則,重置步驟a)中參數(shù)。

c)計算ε=ε+Δε,若定理1、2中條件可行且εgt;0,更新ε、Ω、L,重復(fù)步驟c);否則,進入步驟d)。

d)輸出動態(tài)事件觸發(fā)機制式(10)的參數(shù)(ε,Ω)、濾波器參數(shù)L以及H-/H∞性能指標(biāo)γmin和βmax。

e)選取殘差評價函數(shù)J(k)與閾值Jth,并通過式(21)進行故障檢測決策,判斷故障是否發(fā)生。

3 仿真實驗與結(jié)果分析

本章以攪拌釜式反應(yīng)器系統(tǒng)[24]為例來說明所提方法的有效性,其平衡方程為

dCAdt=FV(CA0-CA)-k1CAe-ERT

dCBdt=-FVCB+k1CAe-ERT-k2CBe-ERT

dTdt=FV(T0-T)+kωARρVCp(Tc-T)-k1CAΔHABR+k2CBΔHBCRρCpe-ERT(40)

其中:CA為反應(yīng)物濃度;CA0為進料濃度;CB為生成物濃度;CB為生成物濃度;T為反應(yīng)溫度;Tc為冷卻劑溫度;V為反應(yīng)器容積;F為容積流量;ρ為液體密度;E/R為反應(yīng)激活能;AR為熱交換系數(shù);ΔH為反應(yīng)熱。選取x(t)=[x1(t) x2(t) x3(t)]T作為狀態(tài)變量,其中x1(t)表示進料物A在t時刻的濃度,x2(t)表示出料物B在t時刻的濃度,x3(t)表示t時刻反應(yīng)室的溫度。圖2為攪拌釜式反應(yīng)器系統(tǒng)的示意圖。

取步長Δt=0.1 s對系統(tǒng)式(40)進行歐拉離散化處理,可得形如系統(tǒng)式(1)的狀態(tài)空間表達式:

x(k+1)=Ax(k)+Φ(k,x(k))+Bw(k)+Ff(k)

y(k)=Cx(k)+Dw(k)

其中:參數(shù)矩陣為A=0.387 20.022 20.018 3

0.244 40.389 70.000 7

-0.068 50.971 10.400 8,B=0.009 5

0.647 4

0.677 9,F(xiàn)=-0.15

0.07

-0.01,C=[1 0 0],D=0.2。

仿真過程中,假定執(zhí)行器故障發(fā)生在低頻域|θ|≤π/10,選取參數(shù)λ1=0.1,λ2=0.3,σ1=0.1,σ2=2,ε0=3,κ=1,非線性函數(shù)Φ(x(k))=0.5(S1+S2)x(k)+0.5(S2-S1)x(k)·sin(x(k)), sin(x(k)):=diag{sin(x1(k)),sin(x2(k)),sin(x3(k))},S1={0.1,0.15,0.2},S2={0.05,0.1,0.15},H=[1 0 0],隨機網(wǎng)絡(luò)攻擊的函數(shù)為g(x(k))=-tanh[0.3xT1(k) 0.2xT2(k) 0.1xT3(k)]T,系統(tǒng)的狀態(tài)初值與濾波器初值分別為x(0)=[0 -1 0.8]T,(0)=[0 0 0]T,外部擾動信號為

w(k)=0.2e-0.02ksin(0.8k) 0≤klt;500(41)

通過求解凸優(yōu)化問題式(40),可以得到最優(yōu)性能指標(biāo)γ=0.083 9與β=0.819 4,濾波器增益矩陣為L=[-0.252 9 0.021 7 -0.441 8]T。此外,為了驗證本文方法的有效性,將故障檢測結(jié)果與不考慮故障頻率特性的H∞濾波器(參考文獻[16])進行對比,選取與上述有限頻域方法相同的γ值,即γ=0.083 9,求得濾波器增益矩陣為Lfull=[-0.472 0 0.021 7 -0.079 5]T。

本章主要考慮兩種類型的執(zhí)行器故障,分別為突變故障和緩變故障。首先,假設(shè)執(zhí)行器在時刻130≤klt;300內(nèi)發(fā)生幅值為0.4的突變故障,表達式為

f(k)=0.4130≤klt;300

0其他(42)

圖3給出了殘差信號在有故障和無故障時的波動情況,對比可以看出在系統(tǒng)發(fā)生故障的時間段130≤klt;300內(nèi),殘差信號有較為明顯的波動。圖4表示在突變故障發(fā)生前后殘差評價函數(shù)J(k)隨時間變化的曲線。系統(tǒng)在k=130時刻發(fā)生故障后,J(k)逐漸超過閾值Jth,實線表示有限頻域故障檢測方法下的殘差評價函數(shù),虛線表示無故障時的殘差評價函數(shù),根據(jù)殘差評價機制式(19)(20)進行計算比較可得J(134)=0.364 9gt;Jth=0.364 3,即4個時間步長內(nèi)檢測出故障的發(fā)生。黑色點畫線表示全頻域故障檢測方法下的殘差評價函數(shù),6個時間步長后檢測出故障。如圖5所示,在500個評估時刻內(nèi),采樣的次數(shù)為500次,而傳感器發(fā)送數(shù)據(jù)的次數(shù)為217次,意味著數(shù)據(jù)傳輸率僅有43.4%,節(jié)約了56.6%的網(wǎng)絡(luò)資源。

其次,考慮執(zhí)行器在時刻k=130后發(fā)生緩變故障,表達式為

f(k)=0.005(k-130)130≤klt;200

0.35+0.03sin(0.3π(k-200))200≤klt;300

0其他(43)

與突變故障類似,緩變故障同樣處于低頻域段。在系統(tǒng)發(fā)生故障的時間段130≤klt;300內(nèi),圖6中的殘差信號有明顯的波動。圖7表示在緩變故障發(fā)生前后殘差評價函數(shù)J(k)隨時間變化的曲線。系統(tǒng)在k=130時刻發(fā)生故障后,J(k)逐漸超過閾值Jth,紅色實線表示有限頻域故障檢測方法(H-/H∞濾波器) 下的殘差評價函數(shù),藍色虛線表示無故障時的殘差評價函數(shù)(見電子版),根據(jù)殘差評價機制式(19)(20)進行計算比較可得J(135)=0.954 7gt;Jth=0.954 2,即5個時間步長內(nèi)檢測出故障的發(fā)生。黑色點畫線表示全頻域故障檢測方法(H∞濾波器) 下的殘差評價函數(shù),12個時間步長后檢測出故障。結(jié)合圖4與7可知,相比于傳統(tǒng)的全頻域方法[16]采用的H∞濾波器,本文中有限頻域故障檢測方法選用的H-/H∞濾波器能夠充分考慮故障的頻率特性,所生成的殘差對故障有更大的敏感度且能更快速地檢測出故障的發(fā)生。

此外,圖8給出了緩變故障下的觸發(fā)時刻與觸發(fā)間隔。如圖8所示,500個評估時刻內(nèi),采樣的次數(shù)為500次,而傳感器發(fā)送數(shù)據(jù)的次數(shù)為208次,意味著數(shù)據(jù)傳輸率僅有41.6%,節(jié)約了58.4%的網(wǎng)絡(luò)資源。圖9、10給出在緩變故障的影響下,周期性事件觸發(fā)機制與靜態(tài)事件觸發(fā)機制下的觸發(fā)時刻與觸發(fā)間隔。對比圖8~10可知,事件觸發(fā)機制避免了每個采樣周期都傳輸測量輸出數(shù)據(jù),且動態(tài)事件觸發(fā)機制能進一步降低傳輸率。相應(yīng)地,表2給出了本文方法與不同傳輸機制下的故障檢測方法的數(shù)據(jù)傳輸率對比。可以看出,在突變與緩變故障情況下本文采用的動態(tài)事件觸發(fā)機制均能更顯著地節(jié)約網(wǎng)絡(luò)資源。

綜上所述,本文所設(shè)計的H-/H∞故障檢測濾波器可以將故障的有限頻域特性考慮到H-性能指標(biāo)的設(shè)計中,且快速地檢測網(wǎng)絡(luò)攻擊下攪拌釜式反應(yīng)器系統(tǒng)中故障的發(fā)生。表3清晰地呈現(xiàn)了本文仿真實驗數(shù)據(jù)。此外,引入的動態(tài)事件觸發(fā)傳輸策略可以進一步地減少輸出測量數(shù)據(jù)向濾波器的發(fā)送量,更顯著地節(jié)約網(wǎng)絡(luò)資源。

4 結(jié)束語

本文提出了一種動態(tài)事件觸發(fā)機制下非線性網(wǎng)絡(luò)系統(tǒng)的有限頻域故障檢測方法。在考慮扇區(qū)有界非線性和隨機網(wǎng)絡(luò)攻擊的情況下,將故障的有限頻域特性考慮到H-性能指標(biāo)的設(shè)計中,利用Lyapunov穩(wěn)定性理論和線性矩陣不等式方法得到故障檢測濾波器存在的充分條件,并使用凸優(yōu)化方法獲取最優(yōu)濾波器參數(shù)。仿真結(jié)果表明,與全頻域故障檢測方法相比,有限頻域方法對故障更加敏感,具有更快的檢測速度,且動態(tài)事件觸發(fā)機制能更顯著地節(jié)約網(wǎng)絡(luò)資源。此外,將本文所提出的故障檢測方法推廣至多智能體系統(tǒng)[25]與信息物理系統(tǒng)[26],是下一步需要研究的問題。

參考文獻:

[1]Zhang Wenan,Dong Hui,Guo Ge,et al.Distributed sampled-data filtering for sensor networks with nonuniform sampling periods[J].IEEE Trans on Industrial Informatics,2014,10(2):871-881.

[2]Cuenca A,Antunes D J,Castillo A,et al.Periodic event-triggered sampling and dual-rate control for a wireless networked control system with applications to UAVs[J].IEEE Trans on Industrial Electro-nics,2019,66(4):3157-3166.

[3]Meng Cai,Wang Tianmiao,Chou Wusheng,et al.Remote surgery case:robot-assisted teleneurosurgery[C]//Proc of IEEE International Conference on Robotics and Automation.2004:819-823.

[4]黃可望,劉艷,潘豐.時延網(wǎng)絡(luò)化控制系統(tǒng)的量化輸出反饋耗散控制[J].計算機應(yīng)用研究,2018,35(10):3053-3056.(Huang Kewang,Liu Yan,Pan Feng.Quantized output feedback dissipative control networked control systems with time-delay[J].Application Research of Computers,2018,35(10):3053-3056.)

[5]Yan Huaicheng,Qian Fengfeng,Yang Fuwen,et al.Filtering for nonlinear networked systems with randomly occurring distributed delays,missing measurements and sensor saturation[J].Information Sciences,2016,370-371:772-782.

[6]汪浩,姜順,潘豐.基于Round-Robin協(xié)議網(wǎng)絡(luò)化系統(tǒng)的故障檢測[J].信息與控制,2019,48(5):595-602.(Wang Hao,Jiang Shun,Pan Feng.Fault detection for networked control systems based on Round-Robin protocol[J].Information and Control,2019,48(5):595-602.)

[7]Li Hongyi,Gao Yabin,Shi Peng,et al.Observer-based fault detection for nonlinear systems with sensor fault and limited communication capacity[J].IEEE Trans on Automatic Control,2016,61(9):2745-2751.

[8]黎煊,吳曉蓓.具有長時延和丟包的網(wǎng)絡(luò)控制系統(tǒng)的故障檢測[J].計算機工程與應(yīng)用,2008,44(33):221-223.(Li Xuan,Wu Xiaobei.Fault detection for networked control systems with long time-delay and packet dropout[J].Computer Engineering and Applications,2008,44(33):221-223.)

[9]劉健辰.事件觸發(fā)通信機制下的有限頻故障檢測[J].信息與控制,2017,46(1):13-18.(Liu Jianchen.Fault detection in the finite frequency domain with an event-triggered communication scheme[J].Information and Control,2017,46(1):13-18.)

[10]李煒,陳文婧.事件觸發(fā)非均勻傳輸NCS少保守性魯棒主動優(yōu)化選擇容錯控制[J].計算機應(yīng)用研究,2019,36(7):2049-2053.(Li Wei,Chen Wenjing.Active optimization choice of NCS less conservative robust fault-tolerant controller under event trigger non-uniform transmission[J].Application Research of Computers,2019,36(7):2049-2053.)

[11]Ren Weijian,Sun Shibo,Hou Nan,et al.Event-triggered non-fragilefault detection for discrete time-delayed nonlinear systems with channel fadings[J].Journal of the Franklin Institute,2018,355(1):436-457.

[12]李艷輝,陶瑩瑩.事件觸發(fā)機制下分布時滯網(wǎng)絡(luò)化控制系統(tǒng)H∞故障檢測[J].控制與決策,2019,35(12):3059-3065.(Li Yanhui,Tao Yingying.Event-triggered H∞ fault detection for networked control systems with distributed delays[J].Control and Decision,2019,35(12):3059-3065.)

[13]Mousavi S H,Marquez H J.Integral-based event triggering controller design for stochastic LTI systems via convex optimization[J].International Journal of Control,2016,89(7):1416-1427.

[14]Dolk V S,Borgers D P,Heemels W P M H.Output-based and decentralized dynamic event-triggered control with guaranteed gain perfor-mance and Zeno-freeness[J].IEEE Trans on Automatic Control,2017,62(1):34-49.

[15]Liu Yan,Yang Guanghong,Li Xiaojian.Event-triggered fault detection observer design for affine fuzzy systems[J].Neurocomputing,2017,267:564-571.

[16]Ning Zhaoke,Wang Tong,Song Xiaona,et al.Fault detection of nonlinear stochastic systems via a dynamic event-triggered strategy[J].Signal Processing,2020,167:107283.

[17]Wang Xudong,F(xiàn)ei Zhongyang,Yan Huangcheng,et al.Dynamic event-triggered fault detection via zonotopic residual evaluation and its application to vehicle lateral dynamics[J].IEEE Trans on Industrial Informatics,2020,16(11):6952-6961.

[18]Hu Songlin,Yue Dong,Xie Xiangpeng,et al.Resilient event-triggered controller synthesis of networked control systems under periodic DoS jamming attacks[J].IEEE Trans on Cybernetics,2019,49(12):4271-4281.

[19]Iwasaki T,Hara S.Generalized KYP lemma:unified frequency domain inequalities with design applications[J].IEEE Trans on Automatic Control,2005,50(1):41-59.

[20]Gu Ying,Yang Guanghong.Event-triggered fault detection for discrete-time Lipschitz nonlinear networked systems in finite-frequency domain[J].Neurocomputing,2017,260:245-256.

[21]Iwasaki T,Hara S,F(xiàn)radkov A L.Time domain interpretations of frequency domain inequalities on(SEMI) finite ranges[J].Systems amp; Control Letters,2005,54(7):681-691.

[22]Gu Ying,Yang Guanghong.Reliable filter design for a class of discrete-time nonlinear systems with time-varying delay[J].Optimal Control Applications and Methods,2009;31(4):303-322.

[23]De Oliveira M C,Geromel J C,Bernussou J.Extended and norm cha-racterizations and controller parametrizations for discrete-time systems[J].International Journal of Control,2002,75(9):666-679.

[24]Chen Bo,Hu Guoqiang,Ho D W C,et al.Distributed robust fusion estimation with application to state monitoring systems[J].IEEE Trans on Systems,Man,and Cybernetics:Systems,2017,47(11):2994-3005.

[25]李富強,豆根生,鄭寶周.基于事件觸發(fā)機制的多智能體網(wǎng)絡(luò)平均一致性研究[J].計算機應(yīng)用研究,2017,34(3):665-670.(Li Fuqiang,Dou Gensheng,Zheng Baozhou.Study on average consistency of multi-agent networks based on event-triggered mechanism[J].Application Research of Computers,2017,34(3):665-670.)

[26]Gu Zhou,Zhou Xiaohong,Zhang Tao,et al.Event-triggered filter design for nonlinear cyber-physical systems subject to deception attacks[J].ISA Trans,2020,104:130-137.

主站蜘蛛池模板: 黄色一级视频欧美| 欧美中文字幕无线码视频| 狠狠亚洲婷婷综合色香| 亚洲精品手机在线| 999国产精品永久免费视频精品久久 | 幺女国产一级毛片| 国产va欧美va在线观看| 婷婷成人综合| 激情无码视频在线看| 伊人91在线| 国产精品福利一区二区久久| 亚洲成aⅴ人在线观看| 国产国模一区二区三区四区| 亚洲欧美成人| 亚洲成人在线网| 国产一二三区视频| 亚洲VA中文字幕| 最新痴汉在线无码AV| 亚洲午夜国产片在线观看| 亚洲精品日产精品乱码不卡| 日本三级欧美三级| 91成人免费观看在线观看| 国产精品自在在线午夜区app| 亚洲精品无码人妻无码| 成人午夜福利视频| 秘书高跟黑色丝袜国产91在线| 亚洲中文精品久久久久久不卡| 色欲色欲久久综合网| 亚洲综合精品第一页| 欧美激情福利| 午夜国产精品视频| 国产欧美网站| 激情无码视频在线看| 亚洲乱强伦| jizz国产在线| 亚洲国产第一区二区香蕉| 一区二区三区精品视频在线观看| 99久久国产精品无码| 一级毛片在线免费视频| 日本高清免费一本在线观看| 亚洲中文字幕手机在线第一页| 亚洲成人精品久久| 国产在线视频欧美亚综合| 中文字幕亚洲无线码一区女同| 91九色国产porny| 无码日韩人妻精品久久蜜桃| 欧美国产日韩在线观看| 精品国产一二三区| 国产精品久久久久久久伊一| 日韩精品视频久久| 国产免费久久精品99re不卡| 九色综合伊人久久富二代| 色综合成人| 天堂av综合网| 免费激情网址| 丁香亚洲综合五月天婷婷| 色成人综合| 国产精品人成在线播放| 亚洲首页在线观看| 欧美成人免费一区在线播放| 2022精品国偷自产免费观看| 久久久久青草线综合超碰| 国产日韩精品一区在线不卡| 亚洲无码在线午夜电影| 国产日韩精品一区在线不卡| 亚洲欧美在线精品一区二区| 亚洲欧美国产高清va在线播放| 亚洲第一区欧美国产综合| 久久精品国产精品一区二区| 国产不卡在线看| 亚洲床戏一区| 久久精品波多野结衣| 中国国语毛片免费观看视频| 国产成人精品免费视频大全五级| 国产欧美视频在线观看| 国内精品视频区在线2021| 8090成人午夜精品| 国产成人无码AV在线播放动漫 | 日韩免费毛片| 久久先锋资源| 欧美日韩精品在线播放| 少妇精品在线|