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

控制系统仿真课程设计（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定义为：+()dtpxFARpxx(4)简记FAR为p1，令p2=1-p1，即2-()dtpxppxx，它表示正常状态下无报警的概率。同理，漏报率MAR定义为：-()dtpxMARqxx(5)简记MAR为q1，令q2=1-q1，即+2()dtpxqqxx，它表示异常状态下发生报警的概率。当过程变量的统计特性精确已知时，显然报警延迟时间Td就成为是一个离散随机变量，其取值属于集合{mh|m=1,2,3,}，那么在概率统计意义下，平均报警延迟AAD定义为Td的期望值，AAD的计算往往假定监控过程变量x(t)是独立同分布的[29]，在此假设下，AAD的具体计算公式为1212=0=0AAD=(=)==/mddmmETmhPTmhmhqqhqq(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)所示()((1)(1)()),,1,ytxtnxtxtntnny(t)为滤波后的变量，与其前n个时刻的x(t)取值有关。这里滑动平均滤波方法经常适用于设备异常的发生改变过程变量均值的情况，而滑动方差滤波方法适用于设备异常的发生改变过程变量方差的情况。假定y(t)的概率密度函数为pY(y),qY(y)，其形式可根据x(t),x(t-1),,x(t-n+1)的概率密度函数给出，那么，滤波报警器的误报率定义为：+()dtpYyFARpyy(7)漏报率定义为：-()dtpyYMARqyy(8)2基于MAR和FAR的报警器最优设计任务内容2.1直接门限法设计步骤步骤一：设定x(t)在正常状态和异常状态下的统计分布211222()~(,)()~(,)xtNxtN正常情况异常情况•其中均值μ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)在正常状态和异常状态下的统计分布21,1,22,2,()~(,)()~(,)xxxxxtNxtN正常情况异常情况•其中均值μ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阶滑动平均滤波器()((2)(1)())3,3,,1,ytxtxtxtntnn则由x(t)的分布可以得到y(t)的分布为21,1,1,1,1,1,022,2,2,2,2,2,0()~(,),=/3()~(,)=,=/3yyyxyxyyyxyxytNttytNtt正常情况异常情况要求：写出针对自己的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曲线Matla