杭電自動(dòng)化控制系統(tǒng)仿真課程設(shè)計(jì)報(bào)告 終極版

本文格式為Word版，下載可任意編輯——杭電自動(dòng)化控制系統(tǒng)仿真課程設(shè)計(jì)報(bào)告終極版

控制系統(tǒng)仿真課程設(shè)計(jì)

（2023屆）

題目學(xué)院專業(yè)班級(jí)學(xué)號(hào)

報(bào)警算法的設(shè)計(jì)和應(yīng)用

自動(dòng)化自動(dòng)化

學(xué)生姓名指導(dǎo)教師完成日期

2023.9.14

控制系統(tǒng)仿真課程設(shè)計(jì)

本課程設(shè)計(jì)的目的著重于DCS中報(bào)警算法的基本原理與設(shè)計(jì)方法的綜合實(shí)際應(yīng)用。主要內(nèi)容包括報(bào)警器性能指標(biāo)MAR、FAR和AAD的求取方法，基于MAR和FAR的最優(yōu)報(bào)警器設(shè)計(jì)方法，直接門限報(bào)警器設(shè)計(jì)方法、滑動(dòng)平均濾波報(bào)警器設(shè)計(jì)方法、基于ROC曲線的最優(yōu)門限選取等Matlab仿真算法的設(shè)計(jì)。

通過本課程設(shè)計(jì)的實(shí)踐，把握自動(dòng)控制理論工程設(shè)計(jì)的基本方法和工具。

1工業(yè)報(bào)警器設(shè)計(jì)中的主要性能指標(biāo)

報(bào)警器對(duì)于保證現(xiàn)代工業(yè)設(shè)備安全有效的運(yùn)行是至關(guān)重要的，那么如何評(píng)價(jià)報(bào)警器的性能使之發(fā)揮不可替代的作用，成為信息與控制領(lǐng)域以及工業(yè)應(yīng)用領(lǐng)域中，研究者和工程技術(shù)人員普遍關(guān)注和研究的問題。尋常來說，報(bào)警器的性能可以用多種指標(biāo)來衡量，譬如“單位時(shí)間警報(bào)數(shù)量〞、“單位時(shí)間峰值警報(bào)數(shù)〞、“單位時(shí)間高（低）優(yōu)先級(jí)警報(bào)數(shù)〞，“警報(bào)識(shí)別率〞等等。其中，被公認(rèn)為最基本和重要的性能指標(biāo)就是其中的誤報(bào)率(FAR)、漏報(bào)率(MAR)、平均報(bào)警延遲(AAD)三特性能指標(biāo)。

1.1三個(gè)主要性能指標(biāo)(FAR/MAR/AAD)的概念與意義

首先令所研究的過程變量記為x，對(duì)它的觀測(cè)為一個(gè)離散的采樣信號(hào)x(t)，采樣周期為h，與其相關(guān)的閾值為xtp。報(bào)警器設(shè)計(jì)實(shí)際上可歸結(jié)為模式分類問題，即x具有“正常(Normal)〞和“異常(Abnormal)〞兩種模式，它們分別對(duì)應(yīng)設(shè)備“無故障〞和“故障〞兩種運(yùn)行狀態(tài)，x(t)經(jīng)過報(bào)警算法處理后會(huì)被映射到“警報(bào)(Alarm)〞和“未警報(bào)(No-alarm)〞兩種模式，亦即x(t)超過閾值xtp則發(fā)出警報(bào)，反之則不發(fā)出警報(bào)。然而，由于x(t)變化的不確定性或xtp的選取不當(dāng)，都會(huì)引起報(bào)警器給出兩種錯(cuò)誤的警報(bào)。一種是誤報(bào)，即在x(t)真實(shí)處于正常狀態(tài)時(shí)而報(bào)警器錯(cuò)誤的發(fā)出警報(bào)；二種是漏報(bào)，即在x(t)真實(shí)處于異常狀態(tài)時(shí)而報(bào)警器未曾發(fā)出警報(bào)。假設(shè)對(duì)于x的一個(gè)采樣序列為{x(0h),x(1h),x(2h),?}，在此次采樣過程中x經(jīng)歷了一次從正常狀態(tài)到異常狀態(tài)的變化，則可以給出一個(gè)混淆矩陣來描述x的兩個(gè)模式和報(bào)警器的兩個(gè)模式之間的關(guān)系，如表1所示。

表1報(bào)警器設(shè)計(jì)中的混淆矩陣

過程變量x的模式異常狀態(tài)正常狀態(tài)誤報(bào)（錯(cuò)誤）個(gè)數(shù)FA未警報(bào)（正確）個(gè)數(shù)TNA個(gè)報(bào)警器模式警報(bào)未警報(bào)警報(bào)（正確）個(gè)數(shù)TA漏報(bào)（錯(cuò)誤）個(gè)數(shù)MNA(A)(NA)當(dāng)xtp變化時(shí)該矩陣中的每個(gè)元素取值都會(huì)發(fā)生變化，并且它們之和為采樣序列{x(0h),x(1h),x(2h),?}的長(zhǎng)度。根據(jù)該矩陣可以給出誤報(bào)率(FAR)和漏報(bào)率(MAR)的定義如下：

FAR=(FA/(FA+TNA))?100%(1)MAR=(MNA/(MNA+TNA))?100%(2)

FAR和MAR是兩個(gè)最重要也是最基本的性能指標(biāo)。報(bào)警器設(shè)計(jì)的最直接也是最簡(jiǎn)單的標(biāo)準(zhǔn)就是找到一個(gè)最優(yōu)閾值，其對(duì)應(yīng)的FAR和MAR都取最小值。這一尋優(yōu)過程可以在接收操作特性(ROC)曲線上進(jìn)行，它描述了當(dāng)閾值取不同數(shù)值時(shí)，相應(yīng)FAR和MAR的變化狀況。最優(yōu)閾值尋常是指ROC曲線上離原點(diǎn)（FAR=0%,MAR=0%）最近的那個(gè)點(diǎn)所對(duì)應(yīng)的閾值。這里舉例說明求取最優(yōu)閾值的過程。

若令設(shè)備發(fā)生異常的時(shí)間為t0=1000h，發(fā)生警報(bào)的時(shí)間ta=1001h，則該組樣本序列下的延遲時(shí)間Td=ta-t0=1h。若有N組如此形式的采樣序列，則可以得到N個(gè)延遲時(shí)間Td1,Td2,?,TdN，那么平均報(bào)警延遲可定義為

AAD=(Td1+Td2+?+TdN)/N(3)

1.2三個(gè)主要性能指標(biāo)的概率表達(dá)式

當(dāng)過程變量x概率分布或密度函數(shù)確切已知的狀況下，不再需要基于樣本數(shù)據(jù)用式(1)-(3)給出三特性能指標(biāo)，而是可以給出FAR、MAR和AAD的概率表達(dá)式，這里將給予詳細(xì)介紹。假設(shè)過程變量x處于正常狀態(tài)的概率密度函數(shù)為p(x)，處于異常狀態(tài)的概率密度函數(shù)為q(x)。

在給定閾值xtp時(shí)，誤報(bào)率FAR定義為：

FAR??簡(jiǎn)記FAR為p1，令p2=1-p1，即p2??同理，漏報(bào)率MAR定義為：

xtp+?xtpp(x)dx(4)

-?p(x)dx，它表示正常狀態(tài)下無報(bào)警的概率。

MAR??q(x)dx(5)

-?xtp簡(jiǎn)記MAR為q1，令q2=1-q1，即

q2??+?xtpq(x)dx，它表示異常狀態(tài)下發(fā)生報(bào)警的概率。

當(dāng)過程變量的統(tǒng)計(jì)特性確切已知時(shí)，顯然報(bào)警延遲時(shí)間Td就成為是一個(gè)離散隨機(jī)變量，其取值屬于集合{mh|m=1,2,3,?}，那么在概率統(tǒng)計(jì)意義下，平均報(bào)警延遲AAD定義為Td的

[29]

期望值，AAD的計(jì)算往往假定監(jiān)控過程變量x(t)是獨(dú)立同分布的，在此假設(shè)下，AAD的具體計(jì)算公式為

AAD?E?Td?=?mh?P(Td=mh)=?mh?q1m?q2=h?q1/q2(6)

m=0m=0??(6)式中h為采樣間隔時(shí)間，為了便利計(jì)算，常假定其值為1s。

同樣的，在過程變量確切聽從一定概率分布的假設(shè)下，常用的濾波方法的FAR、MAR和AAD都可以根據(jù)過程變量的概率密度函數(shù)確切計(jì)算。

濾波方法對(duì)過程變量x做濾波處理，從而消除干擾信號(hào)對(duì)報(bào)警結(jié)果產(chǎn)生的影響。尋常假定濾波器的形式為關(guān)于監(jiān)控變量x的函數(shù)，記為y=f(x)，具體為有限記憶因果濾波器，即y(t)=f(x(t),x(t-1),?,x(t-n+1))，n為濾波器窗口長(zhǎng)度即濾波器階數(shù)。比較常用的有滑動(dòng)平均濾波、滑動(dòng)方差濾波器，分別如式(7)所示

y(t)?(x(t?n?1)???x(t?1)?x(t))n,t?n,n?1,?

y(t)為濾波后的變量，與其前n個(gè)時(shí)刻的x(t)取值有關(guān)。這里滑動(dòng)平均濾波方法經(jīng)常適用

于設(shè)備異常的發(fā)生改變過程變量均值的狀況，而滑動(dòng)方差濾波方法適用于設(shè)備異常的發(fā)生改變過程變量方差的狀況。

假定y(t)的概率密度函數(shù)為pY(y),qY(y)，其形式可根據(jù)x(t),x(t-1),?,x(t-n+1)的概率密度函數(shù)給出，那么，濾波報(bào)警器的誤報(bào)率定義為：

FAR?漏報(bào)率定義為：

?+?ytppY(y)dy(7)

MAR??qY(y)dy(8)

-?ytp2基于MAR和FAR的報(bào)警器最優(yōu)設(shè)計(jì)任務(wù)內(nèi)容

2.1直接門限法設(shè)計(jì)步驟

步驟一：設(shè)定x(t)在正常狀態(tài)和異常狀態(tài)下的統(tǒng)計(jì)分布

2??x(t)~N(?1,?1)正常狀況?2??x(t)~N(?2,?2)異常狀況

?其中均值μ1=0.XX，小數(shù)點(diǎn)后的XX為你的出生日期

?其中均值μ2=1+0.XX，標(biāo)準(zhǔn)差σ1=σ2=1.2

?例如李剛出生日期為1985.5.6日，則XX=06，μ1=0.06?μ2=1+0.X=1.06

要求：寫出針對(duì)自己出生日期的x(t)的分布形式

步驟二：在概率密度已知狀況下，設(shè)定所需遍歷的xtp的取值及個(gè)數(shù)

?設(shè)定xtp取值范圍為區(qū)間I=[μ1-3σ1,μ2+3σ2]?在區(qū)間I中每間隔0.1取一個(gè)點(diǎn)作為xtp的取值，共計(jì)可以得到Num=(μ2+3σ2)-(μ

1-3σ1)/0.1個(gè)xtp的取值點(diǎn)

要求：寫出針對(duì)自己的I，Num取值

步驟三：在概率密度已知狀況下，畫出關(guān)于FAR和MAR的ROC曲線并找到最優(yōu)的xtp取值。

?針對(duì)每一個(gè)xtp的取值，利用定積分式(4)-(5)計(jì)算FAR和MAR的取值?對(duì)于共計(jì)Num個(gè)xtp的取值，共計(jì)可以得到Num組FAR和MAR?畫出ROC曲線，從中找到最優(yōu)的xtp取值

220.5

注：在ROC曲線上尋優(yōu)的準(zhǔn)則為R=(FAR+MAR)取最小值，那個(gè)最小的R值對(duì)應(yīng)的記為最優(yōu)的xtp

要求：畫出ROC曲線

步驟四：利用步驟一中設(shè)定的x(t)的兩個(gè)概率密度函數(shù)，產(chǎn)生正常狀況1000個(gè)點(diǎn)，異常狀況1000個(gè)點(diǎn)，組成采樣序列x(1),x(2),?,x(2000)要求：畫出序列x(t)的圖形

步驟五：依照步驟二中設(shè)定的xtp的取值個(gè)數(shù)，將每個(gè)xtp與x(t)，t=1,2,?,2000,比較，

依照直接門限法的規(guī)則，計(jì)算FAR和MAR，式(1)-(2)，畫出ROC曲線，計(jì)算最優(yōu)的xtp取值，以及相應(yīng)的MAR和FAR

220.5

注：在ROC曲線上尋優(yōu)的準(zhǔn)則為R=(FAR+MAR)取最小值，那個(gè)最小的R值對(duì)應(yīng)的記為最優(yōu)的xtp

要求：畫出ROC曲線

步驟六：比較步驟三和步驟五中計(jì)算出的最優(yōu)xtp取值是否一樣，說明原因？

在參數(shù)mulx1=0.12時(shí)，步驟三xotp的值是0.700,步驟五中是0.800（并隨機(jī)變化），它們是不一樣的。步驟三中采用的方法是理論計(jì)算，理論上說會(huì)比步驟五中的隨機(jī)模擬要確切，但它的步長(zhǎng)值設(shè)為0.1，而差值0.8-0.7=0.1。所以步長(zhǎng)大了。步長(zhǎng)越小越確切，改為0.01后，步驟三輸出xotp=0.700。所以最優(yōu)值為0.700是比較適合的。2.2滑動(dòng)平均濾波法設(shè)計(jì)步驟

步驟一：設(shè)定x(t)在正常狀態(tài)和異常狀態(tài)下的統(tǒng)計(jì)分布

2??x(t)~N(?1,x,?1,x)正常狀況?2??x(t)~N(?2,x,?2,x)異常狀況

?其中均值μ1,x=0.XX，小數(shù)點(diǎn)后的XX為你的出生日期

?其中均值μ2,x=1+0.XX，標(biāo)準(zhǔn)差σ1,x=σ2,x=1.2

?例如李剛出生日期為1985.5.6日，則XX=06，μ1,x=0.06?μ2,x=1+0.X=1.06

要求：寫出針對(duì)自己出生日期的x(t)的分布形式，同2.1中的步驟一

步驟二：設(shè)計(jì)與x(t)對(duì)應(yīng)的3階滑動(dòng)平均濾波器

y(t)?(x(t?2)?x(t?1)?x(t))3,n?3,t?n,n?1,?

則由x(t)的分布可以得到y(tǒng)(t)的分布為

2??y(t)~N(?1,y,?1,y)?1,y??1,x,?1,y=?1,x/3t?t0正常狀況?2??y(t)~N(?2,y,?2,y)?2,y=?2,x,?2,y=?2,x/3t?t0異常狀況要求：寫出針對(duì)自己的y(t)的分布形式

步驟三：在概率密度已知狀況下，設(shè)定所需遍歷的ytp的取值及個(gè)數(shù)

?設(shè)定ytp取值范圍為區(qū)間I=[μ1,y-3σ1,y,μ2,y+3σ2,y]

?在區(qū)間I中每間隔0.1取一個(gè)點(diǎn)作為ytp的取值，共計(jì)可以得到Num=(μ2,y+3σ

2,y)-(μ1,y-3σ1,y)/0.1個(gè)x