本发明涉及过程控制等自动化技术领域，特别是涉及一种适于系统辨识的历史数据段自动查找方法。
背景技术：
在流程工业中，动态过程的数学模型经常用于评价、诊断和提高控制系统的性能。闭环控制系统的性能对于石油、化工和发电等流程工业的安全高效运行有着非常重要的作用。工业调查表明很多工业控制系统面临着性能差的现象，这使得在产品质量、操作成本和过程安全等方面产生了非常严重的负面影响。
在火力发电厂中，给煤量、蒸汽量和给水量等主要过程变量在有功功率的过程中会产生大幅度变化，磨煤机、风机和水泵等相关设备会在多种运行工况下运行，因此这些控制系统的被控对象具有不同的动力学特性。
控制器参数和被控对象之间不匹配会导致控制系统的性能变差，因此在一种运行工况下设计的控制器参数可能不适用于另一个运行工况下的动态过程。为了改善控制系统的性能，需要对控制参数进行重新调整。系统辨识被广泛应用于流程工业中，用于获得被控对象的数学模型，进而根据数学模型重新调整控制参数。
通常可以通过特殊设计的辨识实验，用于辨识动态过程的数学模型。但是，这种系统辨识实验需要格外施加激励信号，会干扰正常的工业生产，所以这种实验在实际中很难进行。
综上所述，现有技术中对如何从控制系统大量的历史运行数据中找到适合系统辨识的数据段，尚缺乏有效的解决方案。
技术实现要素：
为了解决现有技术的不足，本发明提供了一种适于系统辨识的历史数据段自动查找方法，本发明一方面，克服了现有系统辨识实验因需要格外施加激励信号，而会干扰正常的工业生产，并且在实际中很难进行的难题；另一方面，本发明可以用在长时间的数据样本搜索，并且搜索的数据长度可以满足用户的要求。
一种适于系统辨识的历史数据段自动查找方法，包括以下步骤：
获得一系列离散的时间信号组成的数据样本，采用非参数化的概率分布变化检测方法获得数据样本中变化点的位置，找到概率分布存在变化的不同数据段；
判断每一个数据段是否经过一个明显的幅度变化或者为一个恒定值，将后者作为零假设，前者作为替代假设，从而形成假设检验，根据假设检验判断数据段是否能作为系统辨识数据段；
基于假设检验中指示序列都为零的情况下验证两个相邻数据段之间是否存在显著差异；
综合闭环控制系统的设定值获得可用于系统辨识的数据段。
进一步优选的技术方案，所述采用非参数化的概率分布变化检测方法获得数据样本中变化点的位置，具体为：
(1)初始化当前的数据段x(tk:tk+K+1)＝{x(tk),x(tk+1),...,x(tkK-1)},它包含整个数据样本x(1),x(2),Λ,x(N)；
(2)计算一个当前数据段x(tk:tk+K-1)内x(t)的一个相对位置，记作：
其中，t＝tk,tk+1,Λ,tk+K-1，
规定一个统计量U(t)＝U(t-1)+V(t)，初始值U(tk-1)＝0，规定一个假设检验为是否变化位置是时间点U(t)绝对值的最大值，记作：假设检验的P值记作：如果P的值小于第一类错误α的概率值，那么tmax就是变化点的位置，否则，tmax就不是变化点的位置；
(3)把当前数据段x(tk:tk+K-1)分x(tk:tmax)＝{x(tk,x(tk+1),...,x(tmax)}和x(tmax+1:tk+K-1)＝{x(tmax+1),x(tmax+2),...,x(tk+K-1)}两个数据段，tmax是变化的位置；
(4)对每个数据段重复(2)和(3)的步骤直到找不出变化的位置，数据段的长度不小于用户选择的参数K0，作为数据段的最小长度。
进一步优选的技术方案，如果一个数据段长度为K大于2·K0,那么该数据段应该分为K/K0段，经过运算，就能得到所有变化点的位置ti's(i＝1,2,...,I)。
进一步优选的技术方案，如果一个数据段x(ti:ti+1-1)的范围大于一个用户可选择的参数d0，d0为样本显著性变化最小的值，零假设被拒绝，也就是说前一种情况被接受，则x(ti:ti+1-1)可以作为系统辨识数据段。
进一步优选的技术方案，在一个数据段x(ti:ti+1-1)是否作为系统辨识数据段时用指示序列Ix(t)表示数据段能否用来系统辨识，若数据段能够用来系统辨识，则Ix(ti:ti+1-1)＝1；否则Ix(ti:ti+1-1)＝0。
进一步优选的技术方案，验证两个相邻数据段之间是否存在显著差异时，用假设检验来比较两个相邻数据段的平均值。
进一步优选的技术方案，用假设检验来比较两个相邻数据段的平均值时，对于标准高斯分布大于第一类错误类型的一个确定的概率γ，这里γ是一个很小的正数，和Ni分别是样本均值、标准差和数据段的长度，那么，它们的平均差异大于最小显著性变化值d0，如果第二个数据段有一个相对于第一个数据段明显的幅度变化，那么第二段数据的指标Ix的取值从0变为1。
进一步优选的技术方案，如果r(t)在第k个数据段没有变化，那么这个数据段不能用来系统辨识，此时对应的指标序列变为：I(tk,s:tk,e)＝0，用于系统辨识的数据段取值为I(t)＝1，在其它时间点I(t)＝0。
进一步优选的技术方案，本申请的上述适于系统辨识的历史数据段自动查找方法适用于多输入单输出的系统。
与现有技术相比，本发明的有益效果是：
本发明提出了一种从历史数据中搜索适合于系统辨识数据段的新方法，搜索条件是闭环控制系统中设定值和过程输出信号出现了明显的变化。一方面，克服了现有系统辨识实验因需要格外施加激励信号，而会干扰正常的工业生产，并且在实际中很难进行的难题；另一方面，本发明可以用在长时间的数据样本搜索，并且搜索的数据长度可以满足用户的要求。本发明在流程工业对控制回路进行评价、诊断和提高性能中是不可缺少的一部分。
附图说明
构成本申请的一部分的说明书附图用来提供对本申请的进一步理解，本申请的示意性实施例及其说明用于解释本申请，并不构成对本申请的不当限定。
图1为闭环控制系统结构示意图；
图2为本发明的系统辨识历史数据段自动查找方法流程图；
图3(a)至图3(c)分别为r(t),y(t),u(t)和d(t)在24小时内的趋势图；
图4为步骤二中的Iy(t)用虚线表示，步骤三中的Iy(t)用点表示；
图5(a)为y(t)的信号图；
图5(b)为u(t)的信号图，图5(a)与图5(b)用于找第一个变化点的位置；
图6(a)为y(t)∈[7405,8493]比较(实线)和(虚线)；
图6(b)为y(t)∈[8494,0721]比较(实线)和(虚线)；
图7(a)中r(t)实线表示，y(t)点表示，Ir(t)点画线；
图7(b)为u(t)信号图；
图7(c)为d(t)信号图。
具体实施方式
应该指出，以下详细说明都是例示性的，旨在对本申请提供进一步的说明。除非另有指明，本文使用的所有技术和科学术语具有与本申请所属技术领域的普通技术人员通常理解的相同含义。
需要注意的是，这里所使用的术语仅是为了描述具体实施方式，而非意图限制根据本申请的示例性实施方式。如在这里所使用的，除非上下文另外明确指出，否则单数形式也意图包括复数形式，此外，还应当理解的是，当在本说明书中使用术语“包含”和/或“包括”时，其指明存在特征、步骤、操作、器件、组件和/或它们的组合。
考虑到如图1描述的闭环控制系统，这里y和u表示动态过程G的输出和输入。r是y的设定值，d是影响闭环控制系统的外部干扰，v是一个不可测的噪声。控制器C通过不同的y和r改变u的值，q为延迟因子，例如q x(t)＝x(t-1)，如果动态过程G和干扰模型H是线性时不变的，那么一个多输入单输出的系统可表示为：
y(t)＝G(q-1)u(t)+H(q-1)d(t)
系统辨识的目的是从收集的历史数据段中，获得等于或者接近G的数学模型。因此，需要解决的问题是，给定的历史数据段r(t)、y(t)、u(t)其中t＝1,2,Λ,N,如何找出适合辨识G的数据段。本实施例子中收集的数据样本为工业现场的数据即可。
如图2所示，本发明的系统辨识历史数据段自动查找方法运行流程图，包括：
步骤一：找到概率分布存在变化的不同数据段：收集到的数据样本是x(1),x(2),Λ,x(N)等一系列离散的时间信号x(t)。现有变化点位置检测方法是基于先验知识的假设，这些假设在实际中很难验证是否成立。因此采用了一个非参数化的概率分布变化检测方法，获得变化点的位置。该步骤中的数据样本为泛指一般性信号，其中r(t)和y(t)为系统中的具体的实施例子信号。
该方法的具体步骤如下：
(1)初始化当前的数据段x(tk:tk+K+1)＝{x(tk),x(tk+1),...,x(tk+K-1)},它包含整个数据样本x(1),x(2),Λx(N)。
其中，tk表示第k次循环中当前数据段第一个样本时刻，K为当前数据段的样本个数，当k＝1时，tk＝1,即x(tk)＝x(1),x(tk+K-1)＝x(N)。
(2)计算一个当前数据段x(tk:tk+K-1)内x(t)的一个相对位置，记作：
其中，t＝tk,tk+1,...,tk+K-1。
规定一个统计量U(t)＝U(t-1)+V(t)，其中t＝tk,tk+1,...,tk+K-1。初始值U(tk-1)＝0。规定一个假设检验为是否变化位置是时间点U(t)绝对值的最大值，记作：假设检验的P值记作：
如果P的值小于第一类错误α的概率值，那么tmax就是变化点的位置，否则，tmax就不是变化点的位置，假设检验的P值即为第一类错误的概率。
(3)把当前数据段x(tk:tk+K-1)分x(tk:tmax)＝{x(tk,x(tk+1),...,x(tmax)}和x(tmax+1:tk+K-1)＝{x(tmax+1),x(tmax+2),...,x(tk+K-1)}两个数据段，tmax是变化的位置。
(4)对每个数据段重复(2)和(3)的步骤直到找不出变化的位置，数据段的长度不小于用户选择的参数K0，作为数据段的最小长度。另一方面，为了避免过长的数据段，如果一个数据段长度为K大于2·K0,那么该数据段应该分为K/K0段。经过步骤一的运算，就能得到所有变化点的位置ti's(i＝1,2,...,I)，I为变化点的个数。
步骤二：判断每一个数据段是否经过一个明显的幅度变化(即数据段的取值范围大于参数d0)或者为一个恒定值。将后者作为零假设，即数据段为一个恒定值，前者作为替代假设，即数据段存在幅度变化，从而形成假设检验。如果x(ti:ti+1-1)的范围大于一个用户可选择的参数d0(样本显著性变化最小的值)，零假设被拒绝，也就是说替代假设被接受。那么，x(ti:ti+1-1)可以作为系统辨识数据段。在这里，用指示序列Ix(t)表示数据段能否作为系统辨识的数据段，若数据段可以用来系统辨识，那么Ix(ti:ti+1-1)＝1；否则Ix(ti:ti+1-1)＝0。
在本文的零假设中，每个样本跨越均值的概率是相同的，因此跨越均值的两个相邻时刻差服从指数分布。将样本跨越均值的时刻记作tl's,l＝1,2,...,L，L为样本跨越均值时刻的个数,这些时刻由下式定义其中，ti≤tl≤ti+1-1。用Tl＝tl+1-tl表示跨越均值的两个相邻时刻差。定义DT是指数分布函数与估计分布函数的差值的最大绝对值，如下式：
式中，是随机变量T的指数分布函数，表达式如下：
其中，是估计分布函数，表达式如下：
其中，
如果DT大于相应的阈值，那么跨越均值的相邻时刻差就不是指数分布，所以零假设被拒绝。如果x(ti:ti+1-1)的范围大于一个用户可选择的参数d0(样本显著性变化最小的值)。那么，x(ti:ti+1-1)是可用作系统辨识的候选数据段。在这里，用指标序列Ix(t)表示，若数据段可以用来系统辨识，那么Ix(ti:ti+1-1)＝1；否则Ix(ti:ti+1-1)＝0。
步骤三：验证两个相邻数据段之间是否存在显著差异：这一步骤是在第二步中指示序列都为零的情况下进行的。用假设检验来比较两个相邻数据段的平均值。对于标准高斯分布大于第一类错误类型的一个确定的概率γ，这里γ是一个很小的正数。和Ni分别是样本均值，标准差和数据段的长度。那么，它们的平均差异大于最小显著性变化值d0。如果第二个数据段有一个相对于第一个数据段明显的幅度变化，那么第二段数据的指标Ix的取值从0变为1。
步骤四：综合r(t)和y(t)获得可用于系统辨识的数据段。定义一个指示序列I(t)，当Ir(t)或者Iy(t)等于1的时候，I(t)＝1，即：
I(t)＝Ir(t)YIy(t) (7)
其中，用tk,s和tk,e分别表示第k个数据段开始和结束的时间，是I(t)＝1连续的，即：但是，如果r(t)在第k个数据段没有变化，那么这个数据段不能用来系统辨识，此时对应的指标序列变为：I(tk,s:tk,e)＝0。可用于系统辨识的数据段取值为I(t)＝1，在其它时间点I(t)＝0。
需要说明的是，在该实施例子中，这里r(t)有更高的优先级，即r(t)不变，则整个数据段不能使用，不需要再考虑y(t)。
以下是本发明所述方法在具体示例中的应用。
这里给出一个工业应用案例，y(t)是锅炉的水位，r(t)是期望的水位值，u(t)是输出控制器，d(t)是锅炉的主蒸汽流量。y(t)的变化范围是-300mm到300mm，设置高低报警阈值为-100mm、100mm，u(t)是由两个水泵控制流入锅炉的进水量。图3(a)至图3(c)给出了r(t)、y(t)、u(t)和d(t)在24小时内的趋势图。
本方法找到了三个变化的数据段，[2990,9721]，[16549,24553]，[72296,77886]。取t∈[1,12000]的y(t)的时间趋势，图4中，y(t)取0到12000实线，Iy(t)是虚线，竖直的点画线是变化点的位置。图5(a)-图5(b)分别给出了y(t)和u(t)在第一个数据段的具体结果。
步骤一，找到数据段中变化点的位置：图4是t∈[1,12000]的y(t)的时间趋势并通过上述步骤一计算获得变化点的位置(竖直的点画线)，第一个变化点的位置是t＝7404，p＝5.5887-181。DT＝0.2688,阈值取0.2313。
步骤二，图6(a)比较公式9中的和y(t)∈[7405,8493]，DT取最大值0.2688，阈值为0.2313，指示序列Iy(t)在y(t)∈[7405,8493]取1，对比图6(b)是y(t)∈[8494,9721]内y(t)的其他数据，DT＝0.0880，它小于阈值0.0919。因此,在这个数据段内Iy(t)＝0。
步骤三，如图4中步骤二中的Iy(t)用虚线表示，步骤三中的Iy(t)用点表示，分别高于和低于Iy(t)的水平。数据段y(t)∈[2912,4782]是幅度比前一点时间大的数据段。所以另一个数据段是y(t)∈[7405,8493]。因此，步骤三中的Iy(t)比步骤二中的Iy(t)有更加宽的区间。
步骤四，图7(a)-图7(c)给出了两个指示序列Ir(t)和Iy(t)。很明显，Ir(t)围绕在r(t)发生斜坡变化的附近取值1。然而，由于变量d(t)，Iy(t)取1的范围比r(t)大。因此，有必要将Ir(t)和Iy(t)进行“或”运算得到最后的指示序列I(t)。
以上所述仅为本申请的优选实施例而已，并不用于限制本申请，对于本领域的技术人员来说，本申请可以有各种更改和变化。凡在本申请的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本申请的保护范围之内。