杭电自动化控制系统仿真课程设计报告-终极版

1、控制系统仿真课程设计（2012届）题 目报警算法的设计和应用学 院自动化专 业自动化班 级学 号学生姓名指导教师完成日期2015.9.14控制系统仿真课程设计本课程设计的目的着重于DCS中报警算法的基本原理与设计方法的综合实际应用。主要内容包括报警器性能指标MAR、FAR和AAD的求取方法，基于MAR和FAR的最优报警器设计方法，直接门限报警器设计方法、滑动平均滤波报警器设计方法、基于ROC曲线的最优门限选取等Matlab仿真算法的设计。通过本课程设计的实践，掌握自动控制理论工程设计的基本方法和工具。1 工业报警器设计中的主要性能指标报警器对于保证现代工业设备安全有效的运行是至关重要的，那么如

2、何评价报警器的性能使之发挥不可替代的作用，成为信息与控制领域以及工业应用领域中，研究者和工程技术人员普遍关注和研究的问题。通常来说，报警器的性能可以用多种指标来衡量，比如“单位时间警报数量”、“单位时间峰值警报数”、“单位时间高（低）优先级警报数”，“警报识别率”等等。其中，被公认为最基本和重要的性能指标就是其中的误报率(FAR)、漏报率(MAR)、平均报警延迟(AAD)三个性能指标。1.1 三个主要性能指标(FAR/MAR/AAD)的概念与意义首先令所研究的过程变量记为x，对它的观测为一个离散的采样信号x(t)，采样周期为h，与其相关的阈值为xtp。报警器设计实际上可归结为模式分类问题，即x

3、具有“正常(Normal)”和“异常(Abnormal)”两种模式，它们分别对应设备“无故障”和“故障”两种运行状态，x(t)经过报警算法处理后会被映射到“警报(Alarm)”和“未警报(No-alarm)”两种模式，亦即x(t)超过阈值xtp则发出警报，反之则不发出警报。然而，由于x(t)变化的不确定性或xtp的选取不当，都会引起报警器给出两种错误的警报。一种是误报，即在x(t)真实处于正常状态时而报警器错误的发出警报；二种是漏报，即在x(t)真实处于异常状态时而报警器未曾发出警报。假设对于x的一个采样序列为x(0h), x(1h), x(2h),L，在此次采样过程中x经历了一次从正常状态到

4、异常状态的变化，则可以给出一个混淆矩阵来描述x的两个模式和报警器的两个模式之间的关系，如表1所示。表1 报警器设计中的混淆矩阵过程变量x的模式异常状态正常状态报警器模式警报(A)警报（正确）个数TA误报（错误）个数FA未警报(NA)漏报（错误）个数MNA 未警报（正确）个数TNA个当xtp变化时该矩阵中的每个元素取值都会发生变化，并且它们之和为采样序列x(0h), x(1h), x(2h),L的长度。根据该矩阵可以给出误报率(FAR)和漏报率(MAR)的定义如下：FAR=(FA/(FA+TNA)100% (1)MAR=(MNA/(MNA+TNA) 100% (2)FAR和MAR是两个最重要也是

5、最基本的性能指标。报警器设计的最直接也是最简单的标准就是找到一个最优阈值，其对应的FAR和MAR都取最小值。这一寻优过程可以在接收操作特性(ROC)曲线上进行，它描述了当阈值取不同数值时，相应FAR和MAR的变化情况。最优阈值通常是指ROC曲线上离原点（FAR=0%,MAR=0%）最近的那个点所对应的阈值。这里举例说明求取最优阈值的过程。若令设备发生异常的时间为t0=1000h ，发生警报的时间ta=1001h，则该组样本序列下的延迟时间Td=ta-t0=1h。若有N组如此形式的采样序列，则可以得到N个延迟时间Td1,Td2, L,TdN，那么平均报警延迟可定义为AAD=(Td1+ Td2+L

6、+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就

7、成为是一个离散随机变量，其取值属于集合mh|m=1,2,3,L，那么在概率统计意义下，平均报警延迟AAD定义为Td的期望值，AAD的计算往往假定监控过程变量x(t)是独立同分布的29，在此假设下，AAD的具体计算公式为 (6)(6)式中h为采样间隔时间，为了方便计算，常假定其值为1s。同样的，在过程变量精确服从一定概率分布的假设下，常用的滤波方法的FAR、MAR和AAD都可以根据过程变量的概率密度函数精确计算。滤波方法对过程变量x做滤波处理，从而消除干扰信号对报警结果产生的影响。通常假定滤波器的形式为关于监控变量x的函数，记为y=f (x)，具体为有限记忆因果滤波器，即y(t)=f(x(t),

8、x(t-1) ,L, x(t-n+1)，n为滤波器窗口长度即滤波器阶数。比较常用的有滑动平均滤波、滑动方差滤波器，分别如式(7)所示y(t)为滤波后的变量，与其前n个时刻的x(t)取值有关。这里滑动平均滤波方法经常适用于设备异常的发生改变过程变量均值的情况，而滑动方差滤波方法适用于设备异常的发生改变过程变量方差的情况。 假定y(t)的概率密度函数为pY(y), qY(y)，其形式可根据x(t),x(t-1),L, x(t-n+1)的概率密度函数给出，那么，滤波报警器的误报率定义为： (7)漏报率定义为： (8)2 基于MAR和FAR的报警器最优设计任务内容2.1 直接门限法设计步骤步骤一：设定

9、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-31, 2+32 在区间I中每间隔0.1取一个点作为xtp的取值，共计可以得到 Num=(2+32)-(1-31)/0.1个xtp的取值点要求：写出针对自己的I，Num取值步骤三：在概率密度已知情况下，画出关于FAR和MA

10、R的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=

11、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。所以步长大了。

12、步长越小越精确，改为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)的分布形式步骤三

13、：在概率密度已知情况下，设定所需遍历的ytp的取值及个数 设定ytp取值范围为区间I=1,y-31,y, 2,y+32,y 在区间I中每间隔0.1取一个点作为ytp的取值，共计可以得到 Num=(2,y+32,y)-( 1,y-31,y)/0.1个xtp的取值点要求：写出针对自己的I，Num取值步骤四：在概率密度已知情况下，画出关于FAR和MAR的ROC曲线并找到最优的ytp取值。 针对每一个ytp的取值，利用定积分式(4)-(5)计算FAR和MAR的取值 对于共计Num个ytp的取值，共计可以得到Num组FAR和MAR 画出ROC曲线，从中找到最优的xtp取值注：在ROC曲线上寻优的准则为R

14、=(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(

15、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是比较合适的。实验图像：Fig

16、ure1：定积分法求滑动平均滤波法的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的随机数组。%=滑动平均滤波法设计=%clcclear %清空所有历史程序的记录%步骤一：设定x(t)在正常状态和异常状态下的统计分布,即x(t)的正态分布参数m

17、ux1=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/(30.5); %等于原来的值除以根号下3sigmay2=sigmax2/(30.5); %等于原来的值除以根号下3%步骤三：在概率密度已知情况下，设

18、定所需遍历的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+32)-(1-31)/0.1个xtp的取值点%步骤四：在概率密度已知情况下，画出关于FAR和MAR的ROC曲线并找到最优的ytp取值for i=1:Num %对每一个ytp数组中的值循环 FAR(i)=1-cdf(norm,ytp(i),mu

19、y1,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)的距离ends,sn=sort(R); %将R数组取值从小到大进行排序生成序列s，sn为s中元素在原来R中的位置yotp=ytp(sn(1); %sn(1)为yotp在ytp

20、向量中的位置，是最优的ytp取值display(yotp); %向控制台输出yotp值oFAR=FAR(sn(1); %计算最优阈值下的最优FARoMAR=MAR(sn(1); %计算最优阈值下的最优MAR%=%=定积分法求滑动平均滤波法的ROC曲线=%=%figure (1) %画出ROC曲线图plot(FAR,MAR,linewidth,2) %设定ROC曲线的X和Y轴分别为FAR和MARtitle(定积分法求滑动平均滤波法的ROC曲线) %加注图的标题text(FAR(sn(1),MAR(sn(1), leftarrow 最优阈值点,FontSize,12) %标注最优点xlabel(y

21、(t)的FAR); %加注图X坐标名称FARylabel(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)图像=%=%fi

22、gure (2)plot(x,linewidth,1)title(x(t)的采样数据) %加注图的标题xlabel(t); %加注图X坐标名称tylabel(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);for i=3:length(x) %循环赋值 y(i)=(x(i)+x(i-1)+

23、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和FARn_FAR=zeros(1,Num); %设定对于每个ytp取值，产生的误报警的个数n_FAR，并初始化计数值为0n_MAR=zeros(1,Num); %设定对于每个ytp取值

24、，产生的漏报警的个数n_MAR，并初始化计数值为0for i=1:Num %对每个ytp的值，都对正常的1000个y(t)的值比较一遍 for j=1:1000 if y(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;添加于此，也可以endfor i=1:Num %对每个ytp的值，都对正异常的1000个y(t)的值比较一遍 for j=1001:200

25、0 if y(j)ytp(i) %当y处于异常状态，但是y取值没有超过ytp（漏报） 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下的漏报率for i=1:Num R1(i)=(FAR1(i)2+(MAR1(i)2)0.5; %计算ROC曲线上每个点和原点(FAR=0，MAR=0)的距离,即综

26、合值ends1,sn1=sort(R1); %将R1取值从小到大进行排序生成序列s1，sn1()返回s1中元素在原来R中的位置；yotp1=ytp(sn1(1); %sn1(1)为yotp在ytp向量中的位置，也是进行排序后R数组的最小值display(yotp1); %向控制台输出yotp1值oFAR1=FAR1(sn1(1); %计算最优阈值下的最优FAR1oMAR1=MAR1(sn1(1); %计算最优阈值下的最优MAR1%=%=利用y(t)的具体数值得到ROC图像=%=%figure (4) %画出ROC曲线图plot(FAR1,MAR1,linewidth,2) %设定ROC曲线的X

27、和Y轴分别为FAR1和MAR1title(数值比较法求滑动平均滤波法ROC曲线) %加注图的标题text(FAR1(sn1(1),MAR1(sn1(1), leftarrow 最优阈值点,FontSize,12)xlabel(利用y(t)数值的FAR1); %加注图X坐标名称FAR1ylabel(利用y(t)数值的MAR1); %加注图Y坐标名称MAR1对实验结果的解释：从数学公式中可以看到，滑动滤波的标准差是原来的1/(30.5)，从而造成图像上的变化。变得又高又瘦，重叠部分小了，这就说明FAR和MAR相对减小，效果加强，如图：对比可以看到X和Y都比直接门限法要小。最优化的yotp比xotp也要小，上面已说明数值。3实践总结经过这两个星期的短学期，学到了很多新知识，对MATLAB的使用，对报警器系统更加熟悉。MATLAB