杭电自动化控制系统仿真课程设计报告 终极版.docx
《杭电自动化控制系统仿真课程设计报告 终极版.docx》由会员分享,可在线阅读,更多相关《杭电自动化控制系统仿真课程设计报告 终极版.docx(14页珍藏版)》请在冰豆网上搜索。
杭电自动化控制系统仿真课程设计报告终极版
控制系统仿真课程设计
(2012届)
题目
报警算法的设计和应用
学院
自动化
专业
自动化
班级
学号
学生姓名
指导教师
完成日期
2015.9.14
控制系统仿真课程设计
本课程设计的目的着重于DCS中报警算法的基本原理与设计方法的综合实际应用。
主要内容包括报警器性能指标MAR、FAR和AAD的求取方法,基于MAR和FAR的最优报警器设计方法,直接门限报警器设计方法、滑动平均滤波报警器设计方法、基于ROC曲线的最优门限选取等Matlab仿真算法的设计。
通过本课程设计的实践,掌握自动控制理论工程设计的基本方法和工具。
1工业报警器设计中的主要性能指标
报警器对于保证现代工业设备安全有效的运行是至关重要的,那么如何评价报警器的性能使之发挥不可替代的作用,成为信息与控制领域以及工业应用领域中,研究者和工程技术人员普遍关注和研究的问题。
通常来说,报警器的性能可以用多种指标来衡量,比如“单位时间警报数量”、“单位时间峰值警报数”、“单位时间高(低)优先级警报数”,“警报识别率”等等。
其中,被公认为最基本和重要的性能指标就是其中的误报率(FAR)、漏报率(MAR)、平均报警延迟(AAD)三个性能指标。
1.1三个主要性能指标(FAR/MAR/AAD)的概念与意义
首先令所研究的过程变量记为x,对它的观测为一个离散的采样信号x(t),采样周期为h,与其相关的阈值为xtp。
报警器设计实际上可归结为模式分类问题,即x具有“正常(Normal)”和“异常(Abnormal)”两种模式,它们分别对应设备“无故障”和“故障”两种运行状态,x(t)经过报警算法处理后会被映射到“警报(Alarm)”和“未警报(No-alarm)”两种模式,亦即x(t)超过阈值xtp则发出警报,反之则不发出警报。
然而,由于x(t)变化的不确定性或xtp的选取不当,都会引起报警器给出两种错误的警报。
一种是误报,即在x(t)真实处于正常状态时而报警器错误的发出警报;二种是漏报,即在x(t)真实处于异常状态时而报警器未曾发出警报。
假设对于x的一个采样序列为{x(0h),x(1h),x(2h),},在此次采样过程中x经历了一次从正常状态到异常状态的变化,则可以给出一个混淆矩阵来描述x的两个模式和报警器的两个模式之间的关系,如表1所示。
表1报警器设计中的混淆矩阵
过程变量x的模式
异常状态
正常状态
报警器模式
警报(A)
警报(正确)
个数TA
误报(错误)
个数FA
未警报(NA)
漏报(错误)
个数MNA
未警报(正确)
个数TNA个
当xtp变化时该矩阵中的每个元素取值都会发生变化,并且它们之和为采样序列{x(0h),x(1h),x(2h),}的长度。
根据该矩阵可以给出误报率(FAR)和漏报率(MAR)的定义如下:
FAR=(FA/(FA+TNA))100%
(1)
MAR=(MNA/(MNA+TNA))100%
(2)
FAR和MAR是两个最重要也是最基本的性能指标。
报警器设计的最直接也是最简单的标准就是找到一个最优阈值,其对应的FAR和MAR都取最小值。
这一寻优过程可以在接收操作特性(ROC)曲线上进行,它描述了当阈值取不同数值时,相应FAR和MAR的变化情况。
最优阈值通常是指ROC曲线上离原点(FAR=0%,MAR=0%)最近的那个点所对应的阈值。
这里举例说明求取最优阈值的过程。
若令设备发生异常的时间为t0=1000h,发生警报的时间ta=1001h,则该组样本序列下的延迟时间Td=ta-t0=1h。
若有N组如此形式的采样序列,则可以得到N个延迟时间Td1,Td2,,TdN,那么平均报警延迟可定义为
AAD=(Td1+Td2++TdN)/N(3)
1.2三个主要性能指标的概率表达式
当过程变量x概率分布或密度函数精确已知的情况下,不再需要基于样本数据用式
(1)-(3)给出三个性能指标,而是可以给出FAR、MAR和AAD的概率表达式,这里将给予详细介绍。
假设过程变量x处于正常状态的概率密度函数为p(x),处于异常状态的概率密度函数为q(x)。
在给定阈值xtp时,误报率FAR定义为:
(4)
简记FAR为p1,令p2=1-p1,即
,它表示正常状态下无报警的概率。
同理,漏报率MAR定义为:
(5)
简记MAR为q1,令q2=1-q1,即
,它表示异常状态下发生报警的概率。
当过程变量的统计特性精确已知时,显然报警延迟时间Td就成为是一个离散随机变量,其取值属于集合{mh|m=1,2,3,},那么在概率统计意义下,平均报警延迟AAD定义为Td的期望值,AAD的计算往往假定监控过程变量x(t)是独立同分布的[29],在此假设下,AAD的具体计算公式为
(6)
(6)式中h为采样间隔时间,为了方便计算,常假定其值为1s。
同样的,在过程变量精确服从一定概率分布的假设下,常用的滤波方法的FAR、MAR和AAD都可以根据过程变量的概率密度函数精确计算。
滤波方法对过程变量x做滤波处理,从而消除干扰信号对报警结果产生的影响。
通常假定滤波器的形式为关于监控变量x的函数,记为y=f(x),具体为有限记忆因果滤波器,即y(t)=f(x(t),x(t-1),,x(t-n+1)),n为滤波器窗口长度即滤波器阶数。
比较常用的有滑动平均滤波、滑动方差滤波器,分别如式(7)所示
y(t)为滤波后的变量,与其前n个时刻的x(t)取值有关。
这里滑动平均滤波方法经常适用于设备异常的发生改变过程变量均值的情况,而滑动方差滤波方法适用于设备异常的发生改变过程变量方差的情况。
假定y(t)的概率密度函数为pY(y),qY(y),其形式可根据x(t),x(t-1),,x(t-n+1)的概率密度函数给出,那么,滤波报警器的误报率定义为:
(7)
漏报率定义为:
(8)
2基于MAR和FAR的报警器最优设计任务内容
2.1直接门限法设计步骤
步骤一:
设定x(t)在正常状态和异常状态下的统计分布
•其中均值μ1=0.XX,小数点后的XX为你的出生日期
•其中均值μ2=1+0.XX,标准差σ1=σ2=1.2
•例如李刚出生日期为1985.5.6日,则XX=06,μ1=0.06
•μ2=1+0.X=1.06
要求:
写出针对自己出生日期的x(t)的分布形式
步骤二:
在概率密度已知情况下,设定所需遍历的xtp的取值及个数
•设定xtp取值范围为区间I=[μ1-3σ1,μ2+3σ2]
•在区间I中每间隔0.1取一个点作为xtp的取值,共计可以得到Num=(μ2+3σ2)-(μ1-3σ1)/0.1个xtp的取值点
要求:
写出针对自己的I,Num取值
步骤三:
在概率密度已知情况下,画出关于FAR和MAR的ROC曲线并找到最优的xtp取值。
•针对每一个xtp的取值,利用定积分式(4)-(5)计算FAR和MAR的取值
•对于共计Num个xtp的取值,共计可以得到Num组FAR和MAR
•画出ROC曲线,从中找到最优的xtp取值
注:
在ROC曲线上寻优的准则为R=(FAR2+MAR2)0.5取最小值,那个最小的R值对应的记为最优的xtp
要求:
画出ROC曲线
步骤四:
利用步骤一中设定的x(t)的两个概率密度函数,产生正常情况1000个点,异常情况1000个点,组成采样序列x
(1),x
(2),…,x(2000)
要求:
画出序列x(t)的图形
步骤五:
按照步骤二中设定的xtp的取值个数,将每个xtp与x(t),t=1,2,…,2000,比较,按照直接门限法的规则,计算FAR和MAR,式
(1)-
(2),画出ROC曲线,计算最优的xtp取值,以及相应的MAR和FAR
注:
在ROC曲线上寻优的准则为R=(FAR2+MAR2)0.5取最小值,那个最小的R值对应的记为最优的xtp
要求:
画出ROC曲线
步骤六:
比较步骤三和步骤五中计算出的最优xtp取值是否一样,说明原因?
在参数mulx1=0.12时,步骤三xotp的值是0.700,步骤五中是0.800(并随机变化),它们是不一样的。
步骤三中采用的方法是理论计算,理论上说会比步骤五中的随机模拟要精确,但它的步长值设为0.1,而差值0.8-0.7=0.1。
所以步长大了。
步长越小越精确,改为0.01后,步骤三输出xotp=0.700。
所以最优值为0.700是比较合适的。
2.2滑动平均滤波法设计步骤
步骤一:
设定x(t)在正常状态和异常状态下的统计分布
•其中均值μ1,x=0.XX,小数点后的XX为你的出生日期
•其中均值μ2,x=1+0.XX,标准差σ1,x=σ2,x=1.2
•例如李刚出生日期为1985.5.6日,则XX=06,μ1,x=0.06
•μ2,x=1+0.X=1.06
要求:
写出针对自己出生日期的x(t)的分布形式,同2.1中的步骤一
步骤二:
设计与x(t)对应的3阶滑动平均滤波器
则由x(t)的分布可以得到y(t)的分布为
要求:
写出针对自己的y(t)的分布形式
步骤三:
在概率密度已知情况下,设定所需遍历的ytp的取值及个数
•设定ytp取值范围为区间I=[μ1,y-3σ1,y,μ2,y+3σ2,y]
•在区间I中每间隔0.1取一个点作为ytp的取值,共计可以得到Num=(μ2,y+3σ2,y)-(μ1,y-3σ1,y)/0.1个xtp的取值点
要求:
写出针对自己的I,Num取值
步骤四:
在概率密度已知情况下,画出关于FAR和MAR的ROC曲线并找到最优的ytp取值。
•针对每一个ytp的取值,利用定积分式(4)-(5)计算FAR和MAR的取值
•对于共计Num个ytp的取值,共计可以得到Num组FAR和MAR
•画出ROC曲线,从中找到最优的xtp取值
注:
在ROC曲线上寻优的准则为R=(FAR2+MAR2)0.5取最小值,那个最小的R值对应的记为最优的ytp
要求:
画出ROC曲线
步骤五:
利用下式中的两个概率密度函数,产生正常情况1000个点,异常情况1000个点,组成采样序列x
(1),x
(2),…,x(2000)
要求:
画出序列x(t)的图形
步骤六:
由步骤五中产生的x(t)序列产生y(t)序列
要求:
画出序列y(t)的图形
步骤七:
按照步骤三中设定的ytp的取值个数,将每个ytp与y(t),t=1,2,,…,2000,比较,按照超限即报警的规则,计算FAR和MAR,式
(1)-
(2),画出ROC曲线,计算最优的ytp取值,以及相应的MAR和FAR,
注:
y
(1)=x
(1),y
(2)=x
(2)亦即头两个y的值还取x的取值,在ROC曲线上寻优的准则为R=(FAR2+MAR2)0.5取最小值,那个最小的R值对应的记为最优的ytp
要求:
画出ROC曲线
步骤八:
比较步骤四和步骤七中计算出的最优ytp取值是否一样,说明原因?
在参数muly1=0.12时,步骤四输出yotp=0.6415,步骤七输出0.5415(随机变化),它们不一样。
和上面一样,步骤四中采用的方法是理论计算,理论上说会比步骤七中的随机模拟要精确,但它的步长值设为0.1,而差值也为0.1。
所以步长大了。
步长越小越精确,改为0.01后,步骤三输出yotp=0.6215。
所以最优值为0.6215是比较合适的。
实验图像:
Figure1:
定积分法求滑动平均滤波法的ROC曲线
Figure2:
x(t)的采样数据
Figure3:
基于x(t)数值的y(t)的数据图
Figure4:
数值比较法求滑动平均滤波法ROC曲线
Matlab代码及其解释:
相关函数说明:
CDF(NAME,X,A,B):
NAME代表概率分布曲线类型,X代表积分界限,A,B是相应分布曲线的参数;NORMRND(MU,SIGMA,M,N...):
在以参数MU,SIGMA为参数的概率分布曲线中,随机产生M-N的随机数组。
%==========================滑动平均滤波法设计=============================%
clc
clear%%%清空所有历史程序的记录
%%%步骤一:
设定x(t)在正常状态和异常状态下的统计分布,即x(t)的正态分布参数
mux1=0.12;%%%正常状态的均值mux1=0.12;
mux2=1+mux1;%%%异常状态的均值mux2=1.12;
sigmax1=1.2;%%%正常状态的标准差sigmax1=1.2;
sigmax2=sigmax1;%%%异常状态的标准差sigmax2=1.2;
%%%步骤二:
设计与x(t)对应的3阶滑动平均滤波器,即y(t)随x(t)的正态分布参数;
muy1=mux1;%%%和原来相同
muy2=mux2;%%%和原来相同
sigmay1=sigmax1/(3^0.5);%%%等于原来的值除以根号下3
sigmay2=sigmax2/(3^0.5);%%%等于原来的值除以根号下3
%%%步骤三:
在概率密度已知情况下,设定所需遍历的ytp的取值及个数
%%%ytp为把测量的y(t)转换为警报值的阈值
I_l=muy1-3*sigmay1;%%%设定xtp取值范围为区间I的左端点
I_r=muy2+3*sigmay2;%%%设定xtp取值范围为区间I的右端点
ytp=I_l:
0.1:
I_r;%%%在区间I中每间隔0.1取一个点作为xtp的取值
Num=length(ytp);%%%Num=(μ2+3σ2)-(μ1-3σ1)/0.1个xtp的取值点
%%%步骤四:
在概率密度已知情况下,画出关于FAR和MAR的ROC曲线并找到最优的ytp取值
fori=1:
Num%%%对每一个ytp数组中的值循环
FAR(i)=1-cdf('norm',ytp(i),muy1,sigmay1);%%%利用函数指令cdf计算FAR的定积分,cdf()函数计算从负无穷到ytp(i)下曲线和横坐标轴围成的面积,norm代指正态分布函数
MAR(i)=cdf('norm',ytp(i),muy2,sigmay2);%%%利用函数指令cdf计算MAR的定积分,异常分布曲线中漏掉的部分,即MAR值
R(i)=((FAR(i))^2+(MAR(i))^2)^0.5;%%%计算ROC曲线上每个点和原点(FAR=0,MAR=0)的距离
end
[s,sn]=sort(R);%%%将R数组取值从小到大进行排序生成序列s,sn为s中元素在原来R中的位置
yotp=ytp(sn
(1));%%%sn
(1)为yotp在ytp向量中的位置,是最优的ytp取值
display(yotp);%%%向控制台输出yotp值
oFAR=FAR(sn
(1));%%%计算最优阈值下的最优FAR
oMAR=MAR(sn
(1));%%%计算最优阈值下的最优MAR
%===================================%
%==定积分法求滑动平均滤波法的ROC曲线==%
%===================================%
figure
(1)%%%画出ROC曲线图
plot(FAR,MAR,'linewidth',2)%%%设定ROC曲线的X和Y轴分别为FAR和MAR
title('定积分法求滑动平均滤波法的ROC曲线')%%%加注图的标题
text(FAR(sn
(1)),MAR(sn
(1)),'\leftarrow最优阈值点','FontSize',12)%%%标注最优点
xlabel('y(t)的FAR');%%%加注图X坐标名称FAR
ylabel('y(t)的MAR');%%%加注图Y坐标名称FAR
%%%步骤五:
利用x(t)的概率密度函数,产生正常情况1000个点,异常情况1000个点,组成采样序列,并画出x(t)图形
x_nor=normrnd(mux1,sigmax1,1,1000);%%%产生正常状态下的1000个点(normal),normrnd()函数返回对应正态分布的随机点
x_abnor=normrnd(mux2,sigmax2,1,1000);%%%产生异常状态下的1000个点(abnormal)
x=[x_nor,x_abnor];
%===================================%
%====利用正态分布随机值的x(t)图像====%
%===================================%
figure
(2)
plot(x,'linewidth',1)
title('x(t)的采样数据')%%%%加注图的标题
xlabel('t');%%%加注图X坐标名称t
ylabel('x(t)');%%%加注图Y坐标名称x(t)
%%%步骤六:
由步骤五中产生的x(t)序列产生y(t)序列,并画出y(t)图形;
%%%%其中y(t)=((x(t)+x(t-1)+x(t-2))/3;其中n=3,t=n,n+1,n+2......
y=zeros(1,length(x));%%%初始化y数组的2000个值
y
(1)=x
(1);%%%y(t)的头几个值和x(t)相同
y(20)=x
(2);
fori=3:
length(x)%%%循环赋值
y(i)=(x(i)+x(i-1)+x(i-2))/3;
end
%===================================%
%===利用随机x(t)的值产生的y(t)图像===%
%===================================%
figure(3)
plot(y,'linewidth',1)
title('基于x(t)数值的y(t)的数据图')
xlabel('t')
ylabel('y(t)')
%%%步骤七:
按照步骤三中设定的ytp的取值个数,将每个ytp与y(t),t=1,2,,…,2000比较
%%%按照超限即报警的规则,计算FAR和MAR,画出ROC曲线,计算最优的ytp取值,以及相应的MAR和FAR
n_FAR=zeros(1,Num);%%%设定对于每个ytp取值,产生的误报警的个数n_FAR,并初始化计数值为0
n_MAR=zeros(1,Num);%%%设定对于每个ytp取值,产生的漏报警的个数n_MAR,并初始化计数值为0
fori=1:
Num%%%对每个ytp的值,都对正常的1000个y(t)的值比较一遍
forj=1:
1000
ify(j)>=ytp(i)%%%当y处于正常状态,但是y取值超过ytp(误报),对每个ytp的值,取1000个正常点比较
n_FAR(i)=n_FAR(i)+1;%%%计算在ytp(i)这个阈值下,误报的个数,保存对应ytp的值的误报个数
end
end%%%FAR1(i)=n_FAR(i)/1000;添加于此,也可以
end
fori=1:
Num%%%对每个ytp的值,都对正异常的1000个y(t)的值比较一遍
forj=1001:
2000
ify(j)n_MAR(i)=n_MAR(i)+1;%%%计算在ytp(i)这个阈值下,漏报的个数
end
end%%%MAR1(i)=n_MAR(i)/1000添加与此,也可以
end
FAR1=n_FAR/1000;%%%对数组中的每个值都除以1000,即计算每个ytp下的误报率
MAR1=n_MAR/1000;%%%对数组中的每个值都除以1000,即计算每个ytp下的漏报率
fori=1:
Num
R1(i)=((FAR1(i))^2+(MAR1(i))^2)^0.5;%%%计算ROC曲线上每个点和原点(FAR=0,MAR=0)的距离,即综合值
end
[s1,sn1]=sort(R1);%%%将R1取值从小到大进行排序生成序列s1,sn1()返回s1中元素在原来R中的位置;
yotp1=ytp(sn1
(1));%%%sn1
(1)为yotp在ytp向量中的位置,也是进行排序后R数组的最小值
display(yotp1);%%%向控制台输出yotp1值
oFAR1=FAR1(sn1
(1));%%%计算最优阈值下的最优FAR1
oMAR1=MAR1(sn1
(1));%%%计算最优阈值下的最优MAR1
%===================================%
%===利用y(t)的具体数值得到ROC图像====%
%===================================%
figure(4)%%%画出ROC曲线图
plot(FAR1,MAR1,'linewidth',2)%%%设定ROC曲线的X和Y轴分别为FAR1和MAR1
title('数值比较法求滑动平均滤波法ROC曲线')%%%加注图的标题
text(FAR1(sn1
(1)),MAR1(sn1
(1)),'\leftarrow最优阈值点','FontSize',12)
xlabel('利用y(t)数值的FAR1');%%%加注图X坐标名称FAR1
ylabel('利用y(t)数值的MAR1');%%%加注图Y坐标名称MAR1
对实验结果的解释:
从数学公式中可以看到,滑动滤波的标准差是原来的1/(3^0.5),从而造成图像上的变化。
变得又高又瘦,重叠部分小了,这就说明FAR和MAR相对减小,效果加强,如图:
对比可以看到X和Y都比直接门限法要小。
最优化的yotp比xotp也要小,上面已说明数值。
3实践总结
经过这两个星期的短学期,学到了很多新知识,对MATLAB的使用,对报警器系统更加熟悉。
MATLAB是一个很好的数学处理软件,很简单的方式处理数学问题,代码通俗易懂,是很好的编程平台。
并且在理解报警系统中的随机性问题中,很好的理解了概率问题。
概率是个体和总体的关系,个体身上体现了总体的各种倾向,总体又是由相互独立的个体构成的。
而在处理中加入动态的思想是非常契合实际的,因为在实际中理想的数学模型是不存在的,会有各种预见和不能预见的因素影响实际结果,而动态滤波,在很大程度上克服了这个问题,可以说是以变化的方法来应对概率中的、现实中的变化特性,可以在其他问题,尤其是工程问题中加以借鉴。
实际实验中无法区分测的量的正常异常,可以通过故障实验的方法来的到准确的数据,即观察,凭借人的经验去判定。
如果要完全的自动化,可以对实际中产生的数据进行实时的均值统计,均值曲线的拐点意味着数值规律的变化,比如从正态分布到其它的分布,然后调用相应的处理函数即可。
4致谢
非常感谢***、***这两位老师的帮助和耐心的指导。
5参考文献
[1]宋晓静,***,文成林.一种基于证据理论的工业报警器设计方法[J].杭州电子科技大学学报.2011,31(6):
135-138.
[2]XuJ,WangJ,IzadiI,etal.PerformanceAssessmentandDesignforUnivariateAlarmSystemsBasedonFAR,MAR,andAAD[J].IEEETran