发明名称 一种自动的步进式有序时段划分方法
摘要 本发明公开了一种自动的步进式有序时段划分方法,该方法可以将一个多操作步骤的间歇过程根据过程特性的变化自动划分为不同的子时段,每个时段具有相似的过程特性可以用同一模型来表示,而不同时段具有不同特性,需要建立不同的模型,大大简化了模型数量和复杂度。本发明直接采用在线监测指标作为时段划分的判断依据,改进了时段模型的精度,并大大提高了后续的过程在线监测性能。该方法简单易于实施,在注塑成型过程中得到了成功应用,不仅有利于对具体过程特性的了解,而且增强了实际在线过程监测的可靠性和可信度,有助于工业工程师对过程运行状态做出准确判断,及时发现故障,从而保证实际生产的安全可靠运行和产品的高质量追求。
申请公布号 CN103116306A 申请公布日期 2013.05.22
申请号 CN201310046432.1 申请日期 2013.02.05
申请人 浙江大学 发明人 赵春晖;李文卿
分类号 G05B19/048(2006.01)I 主分类号 G05B19/048(2006.01)I
代理机构 杭州求是专利事务所有限公司 33200 代理人 周烽
主权项 1.一种自动的步进式有序时段划分方法,其特征在于,该方法包括以下步骤:步骤1:获取过程分析数据:设一个间歇操作具有J个测量变量和K个采样点,则每一个测量批次可得到一个K×J的矩阵,重复I批次的测量步骤后,得到的数据可以表述为一个三维矩阵<img file="FDA00002822189400011.GIF" wi="296" he="62" />其中测量变量为温度、速度、压力、位移等批次运行过程中可被测量的状态参数;步骤2:数据预处理:将三维矩阵<img file="FDA00002822189400012.GIF" wi="55" he="58" />按照采集批次方向展开,即将一个操作批次内的各采样点上的变量按照时间顺序排开得到二维矩阵X(I×KJ),由K个时间片矩阵X<sub>k</sub>(I×J)组成,其中,下标k为时间指标;设二维矩阵X<sub>k</sub>内任意一点的变量为x<sub>k,i,j</sub>对该变量进行减均值、除以标准差的标准化处理,其中,下标i代表批次,j代表变量,标准化处理的计算公式如下:<maths num="0001"><![CDATA[<math><mrow><msub><mi>x</mi><mrow><mi>k</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>j</mi></mrow></msub><mo>=</mo><mfrac><mrow><msub><mi>x</mi><mrow><mi>k</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>j</mi></mrow></msub><mo>-</mo><msub><mover><mi>x</mi><mo>&OverBar;</mo></mover><mrow><mi>k</mi><mo>,</mo><mi>j</mi></mrow></msub></mrow><msub><mi>s</mi><mrow><mi>k</mi><mo>,</mo><mi>j</mi></mrow></msub></mfrac><mo>;</mo><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>1</mn><mo>)</mo></mrow></mrow></math>]]></maths>其中:k是时间片指标。<img file="FDA00002822189400014.GIF" wi="76" he="67" />是X<sub>k</sub>矩阵任一列的均值,s<sub>k,j</sub>是X<sub>k</sub>矩相应列的标准差,<maths num="0002"><![CDATA[<math><mrow><msub><mover><mi>x</mi><mo>&OverBar;</mo></mover><mrow><mi>k</mi><mo>,</mo><mi>j</mi></mrow></msub><mo>=</mo><mfrac><mn>1</mn><mi>I</mi></mfrac><munderover><mi>&Sigma;</mi><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>I</mi></munderover><msub><mi>x</mi><mrow><mi>k</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>j</mi></mrow></msub><mo>,</mo></mrow></math>]]></maths><maths num="0003"><![CDATA[<math><mrow><msub><mi>s</mi><mrow><mi>k</mi><mo>,</mo><mi>j</mi></mrow></msub><mo>=</mo><msqrt><munderover><mi>&Sigma;</mi><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>I</mi></munderover><msup><mrow><mo>(</mo><msub><mi>x</mi><mrow><mi>k</mi><mo>,</mo><mi>i</mi><mo>,</mo><mi>j</mi></mrow></msub><mo>-</mo><msub><mover><mi>x</mi><mo>&OverBar;</mo></mover><mrow><mi>k</mi><mo>,</mo><mi>j</mi></mrow></msub><mo>)</mo></mrow><mn>2</mn></msup><mo>/</mo><mrow><mo>(</mo><mi>I</mi><mo>-</mo><mn>1</mn><mo>)</mo></mrow></msqrt><mo>;</mo><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>2</mn><mo>)</mo></mrow></mrow></math>]]></maths>步骤3:时间片PCA建模,该步骤由以下子步骤来实现:(3.1)对步骤2标准化处理后的每一个时间片矩阵X<sub>k</sub>(I×J)执行PCA分解,建立时间片PCA模型,其中PCA分解公式如下:<maths num="0004"><![CDATA[<math><mrow><msub><mi>X</mi><mi>k</mi></msub><mo>=</mo><msub><mi>T</mi><mi>k</mi></msub><msup><msub><mi>P</mi><mi>k</mi></msub><mi>T</mi></msup><mo>=</mo><munderover><mi>&Sigma;</mi><mrow><mi>r</mi><mo>=</mo><mn>1</mn></mrow><mi>J</mi></munderover><msub><mi>t</mi><mrow><mi>k</mi><mo>,</mo><mi>r</mi></mrow></msub><msub><mi>p</mi><mrow><mi>k</mi><mo>,</mo><mi>r</mi></mrow></msub><mo>;</mo><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>3</mn><mo>)</mo></mrow></mrow></math>]]></maths>其中:t<sub>k,n</sub>为正交的主元向量,p<sub>k,n</sub>为正交归一化的负载向量,r表示不同的PCA分解方向,上标T表示矩阵的转置;T<sub>k</sub>(I×J)代表保留全部主元的得分矩阵,P<sub>k</sub>(J×J)代表对应的负载矩阵。(3.2)选取主元个数,将公式(3)重新表述成如下形式:<maths num="0005"><![CDATA[<math><mrow><msub><mi>X</mi><mi>k</mi></msub><mo>=</mo><msub><mi>T</mi><mi>k</mi></msub><msup><msub><mi>P</mi><mi>k</mi></msub><mi>T</mi></msup><mo>+</mo><msub><mi>E</mi><mi>k</mi></msub><mo>=</mo><munderover><mi>&Sigma;</mi><mrow><mi>r</mi><mo>=</mo><mn>1</mn></mrow><mi>R</mi></munderover><msub><mi>t</mi><mrow><mi>k</mi><mo>,</mo><mi>r</mi></mrow></msub><msub><mi>p</mi><mrow><mi>k</mi><mo>,</mo><mi>r</mi></mrow></msub><mo>+</mo><msub><mi>E</mi><mi>k</mi></msub><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>4</mn><mo>)</mo></mrow></mrow></math>]]></maths>其中:r表示不同的PCA分解方向;T<sub>k</sub>(I×R<sub>k</sub>)与P<sub>k</sub>(J×R<sub>k</sub>)代表负载矩阵分别为保留R<sub>k</sub>个主元后的得分矩阵和负载矩阵,E<sub>k</sub>为残差矩阵。通过上述变换,多向主元分析法模型将原始数据空间分解为主元空间和残差空间,主元空间内代表主要的系统过程波动信息;这里所保留的主元个数R<sub>k</sub>能够反映原过程中90%的过程波动信息。综合考虑全过程,选取出现次数最多的R<sub>k</sub>值统一为最终的建模主元个数R,即所有时间片模型保留相同的主元个数。(3.3)计算残差空间中各时间片k中对应各个批次的SPE指标:SPE<sub>k,i</sub>=e<sub>k,i</sub><sup>T</sup>e<sub>k,i</sub>      (5)其中,下标i表示时间片中的批次,e<sub>k,i</sub>是对应k时刻第i批次的残差列向量。根据相同时刻上不同批次的SPE值服从带权重系数的χ<sup>2</sup>分布,从而确定出每个时间点上的控制限Ctr<sub>k</sub>,它反应了时间片PCA模型的重构能力。步骤4:确定基于变量展开模型的SPE指标控制限:从间歇过程初始点开始,依次将下一个时间片与之前的时间片组合在一起并按变量方式展开得到时间块X<sub>v,k</sub>(IK×J),其中下标v代表变量展开方式。对新时间块矩阵进行PCA分析,提取出负载矩阵P<sub>v,k</sub>(J×R)。计算其SPE值并根据相同时刻上不同批次的SPE值服从带权重系数的χ<sup>2</sup>分布,从而确定出每个时间点上的控制限Ctr<sub>k</sub>。步骤5:确定第一时段划分点k<sup>*</sup>:比较在相同时间区域内的每个时间点上Ctr<sub>k</sub>和Ctr<sub>v,k</sub>的大小,如果发现连续三个样本呈现Ctr<sub>v,k</sub>>αCtr<sub>k</sub>,那么新加入的时间片对该时间块的PCA监测模型及相应的监测性能都有重大的影响。记加入新时间片前的时刻为k<sup>*</sup>。其中,α是依附于Ctr<sub>k</sub>的常数,称作缓和因子,它反映的是与时间片模型相比,时间块模型允许监测精度损失的程度。则k<sup>*</sup>时刻之前的时间片可认为是一个子时段。步骤6:过程分析数据更新,确定所有划分时段:根据步骤5中所获得的时刻k<sup>*</sup>的指示,移除第一个子时段,把余下的间歇过程数据作为新的输入数据带入到第5步中并重复上述步骤5-6,划分不同时间段,直到没有数据余留。步骤7:基于时段划分结果的统计建模:根据步骤6时段划分结果,每个时段内的时间片按照变量展开方式组合成子时段代表性建模数据组,X<sub>c</sub>(IK<sub>c</sub>×J),其中,下标c是时段指标。每个时段内相似的过程潜在特性可以通过对X<sub>c</sub>(IK<sub>c</sub>×J)实施PCA分解提取出来:T<sub>c</sub>=X<sub>c</sub>P<sub>c</sub><maths num="0006"><![CDATA[<math><mrow><msub><mover><mi>X</mi><mo>^</mo></mover><mi>c</mi></msub><mo>=</mo><msub><mi>T</mi><mi>c</mi></msub><msup><msub><mi>P</mi><mi>c</mi></msub><mi>T</mi></msup><mo>=</mo><msub><mi>X</mi><mi>c</mi></msub><msub><mi>P</mi><mi>c</mi></msub><msup><msub><mi>P</mi><mi>c</mi></msub><mi>T</mi></msup><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>6</mn><mo>)</mo></mrow></mrow></math>]]></maths><maths num="0007"><![CDATA[<math><mrow><msub><mi>E</mi><mi>c</mi></msub><mo>=</mo><msub><mi>X</mi><mi>c</mi></msub><mo>-</mo><msub><mover><mi>X</mi><mo>^</mo></mover><mi>c</mi></msub></mrow></math>]]></maths>其中,P<sub>c</sub>(J×R<sub>c</sub>)是该子时段的PCA负载矩阵,揭示了该子时段内的主要波动方向,R<sub>c</sub>代表了该子时段模型保留的PCA主元个数。T<sub>c</sub>(IK<sub>c</sub>×R<sub>c</sub>)代表了该子时段的主元得分矩阵。因此,该时段的主要波动通过T<sub>c</sub>P<sub>c</sub><sup>T</sup>表征,代表了该子时段内的主要波动信息;剩余为子时段残差矩阵,E<sub>c</sub>(IK<sub>c</sub>×J)代表了该子时段内的噪声信息。各时间片主元得分,T<sub>k</sub>(I×R<sub>c</sub>),可以很容易从子时段得分矩阵T<sub>c</sub>(IK<sub>c</sub>×R<sub>c</sub>)中根据对应的过程时间抽取获得,从而可以据此计算得出各采样时刻上的协方差关系S<sub>k</sub>。各时间片残差矩阵E<sub>k</sub>(I×J)亦可以从子时段残差矩阵E<sub>c</sub>(IK<sub>c</sub>×J)中对应获得。步骤8:计算实时监测统计指标:根据从公式(6)计算的结果中获取的主元得分时间片T<sub>k</sub>(I×R<sub>c</sub>)和残差矩阵时间片E<sub>k</sub>(I×J)可以在每个时刻计算两个监测统计指标:Hotelling-T<sup>2</sup>统计指标和SPE统计量。Hotelling-T<sup>2</sup>统计指标用来测量各采样时刻过程变量偏离正常工况下平均轨迹的距离,并且其在显著性水平α下的控制限计算为:<maths num="0008"><![CDATA[<math><mrow><msubsup><mi>T</mi><mrow><mi>i</mi><mo>,</mo><mi>k</mi></mrow><mn>2</mn></msubsup><mo>-</mo><msup><mrow><mo>(</mo><msub><mi>t</mi><mrow><mi>i</mi><mo>,</mo><mi>k</mi></mrow></msub><mo>-</mo><msub><mover><mi>t</mi><mo>&OverBar;</mo></mover><mi>k</mi></msub><mo>)</mo></mrow><mi>T</mi></msup><msubsup><mi>S</mi><mi>k</mi><mrow><mo>-</mo><mn>1</mn></mrow></msubsup><mrow><mo>(</mo><msub><mi>t</mi><mrow><mi>i</mi><mo>,</mo><mi>k</mi></mrow></msub><mo>-</mo><msub><mover><mi>t</mi><mo>&OverBar;</mo></mover><mi>k</mi></msub><mo>)</mo></mrow><mo>~</mo><mfrac><mrow><msub><mi>R</mi><mi>c</mi></msub><mrow><mo>(</mo><mi>I</mi><mo>-</mo><mn>1</mn><mo>)</mo></mrow></mrow><mrow><mi>I</mi><mo>-</mo><msub><mi>R</mi><mi>c</mi></msub></mrow></mfrac><msub><mi>F</mi><mrow><msub><mi>R</mi><mi>c</mi></msub><mo>,</mo><msub><mrow><mi>I</mi><mo>-</mo><mi>R</mi></mrow><mi>c</mi></msub><mo>,</mo><mi>&alpha;</mi></mrow></msub><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>7</mn><mo>)</mo></mrow></mrow></math>]]></maths>其中,R<sub>c</sub>是该时段PCA模型中所保留的主元个数;t<sub>i,k</sub>(R<sub>c</sub>×1)是第k时刻第i个批次的主元得分,即对应时间片得分矩阵T<sub>k</sub>(I×R<sub>c</sub>)的第i行;而<img file="FDA00002822189400032.GIF" wi="189" he="58" />是T<sub>k</sub>(I×R<sub>c</sub>)针对不同批次的均值向量,由于各时间片测量数据已经在数据预处理时中心化为零均值,这里的<img file="FDA00002822189400033.GIF" wi="190" he="58" />其实就是零向量。对于残差子空间,各个时刻不同批次的SPE统计量计算为:<maths num="0009"><![CDATA[<math><mrow><msub><mi>SPE</mi><mrow><mi>i</mi><mo>,</mo><mi>k</mi></mrow></msub><msup><mrow><mo>(</mo><msub><mi>x</mi><mrow><mi>i</mi><mo>,</mo><mi>k</mi></mrow></msub><mo>-</mo><msub><mover><mi>x</mi><mo>^</mo></mover><mrow><mi>i</mi><mo>,</mo><mi>k</mi></mrow></msub><mo>)</mo></mrow><mi>T</mi></msup><mrow><mo>(</mo><msub><mi>x</mi><mrow><mi>i</mi><mo>,</mo><mi>k</mi></mrow></msub><mo>-</mo><msub><mover><mi>x</mi><mo>^</mo></mover><mrow><mi>i</mi><mo>,</mo><mi>k</mi></mrow></msub><mo>)</mo></mrow><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>8</mn><mo>)</mo></mrow></mrow></math>]]></maths>其中,x<sub>i,k</sub>代表k时刻第i个批次的过程测量向量,<img file="FDA00002822189400035.GIF" wi="58" he="58" />则是由PCA模型重构的对应结果。式(8)中的计算结果在每个时刻都可以构成一个I×1向量[SPE<sub>1,k</sub>,SPE<sub>2,k</sub>,...,SPE<sub>I,k</sub>]<sup>T</sup>,并且该向量近似服从加权χ<sup>2</sup>分布,从而获得每个时刻SPE的监测控制限。步骤9:基于时段模型的在线过程监测:基于步骤6划分的时段、步骤7建立的时段模型监测系统和步骤8所得的两个监测统计量可以在线监测注塑成型等新运行间歇过程的状态。该步骤由以下子步骤来实现:(9.1)采集新测量数据及新测量数据预处理:在线监测时,采集到新的过程测量数据x<sub>new</sub>(J×1)后,其中,下标new代表新样本,J为测量变量,与步骤1中的测量变量相同;首先需要进行数据预处理。根据步骤2中获得的均值和标准差,根据过程时间的指示调用对应该时刻的均值和标准差对采集到的过程测量数据如公式(1)中所示进行标准化预处理。(9.2)计算新监测统计量:数据预处理后,根据公式(6)计算的PCA子时段模型,调用对应该新采样时刻所在时段的模型P<sub>c</sub>(J×R<sub>c</sub>),按照如下方式计算得到主元得分,估计残差及其对应的Hotelling-T<sup>2</sup>与SPE两个监测统计指标:t<sub>new</sub><sup>T</sup>=x<sub>new</sub><sup>T</sup>P<sub>c</sub><maths num="0010"><![CDATA[<math><mrow><msub><mover><mi>x</mi><mo>^</mo></mover><mi>new</mi></msub><mo>=</mo><msub><mi>P</mi><mi>c</mi></msub><msub><mi>t</mi><mi>new</mi></msub></mrow></math>]]></maths><maths num="0011"><![CDATA[<math><mrow><mo>(</mo><mn>9</mn><mo>)</mo></mrow></math>]]></maths><maths num="0012"><![CDATA[<math><mrow><msubsup><mi>T</mi><mi>new</mi><mn>2</mn></msubsup><mo>=</mo><msup><mrow><mo>(</mo><msub><mi>t</mi><mi>new</mi></msub><mo>-</mo><msub><mover><mi>t</mi><mo>&OverBar;</mo></mover><mi>k</mi></msub><mo>)</mo></mrow><mi>T</mi></msup><msup><msub><mi>S</mi><mi>k</mi></msub><mrow><mo>-</mo><mn>1</mn></mrow></msup><mrow><mo>(</mo><msub><mi>t</mi><mi>new</mi></msub><mo>-</mo><msub><mover><mi>t</mi><mo>&OverBar;</mo></mover><mi>k</mi></msub><mo>)</mo></mrow></mrow></math>]]></maths><maths num="0013"><![CDATA[<math><mrow><msub><mi>SPE</mi><mi>new</mi></msub><mo>=</mo><msup><mrow><mo>(</mo><msub><mi>x</mi><mi>new</mi></msub><mo>-</mo><msub><mover><mi>x</mi><mo>^</mo></mover><mi>new</mi></msub><mo>)</mo></mrow><mi>T</mi></msup><mrow><mo>(</mo><msub><mi>x</mi><mi>new</mi></msub><mo>-</mo><msub><mover><mi>x</mi><mo>^</mo></mover><mi>new</mi></msub><mo>)</mo></mrow></mrow></math>]]></maths>其中<img file="FDA00002822189400045.GIF" wi="85" he="66" />是新的过程测量数据,<img file="FDA00002822189400046.GIF" wi="46" he="66" />是前面根据训练数据获得的T<sub>k</sub>的均值向量,S<sub>k</sub>是X<sub>k</sub>的协方差矩阵。(9.3)在线判断过程运行状态:实时比较两个监测指标与其各自的统计控制限,若两个监测指标都位于统计控制限之内,表明过程运行正常;若有一个以上监测指标超出正常控制限,表明过程有异常状况发生。
地址 310058 浙江省杭州市西湖区余杭塘路866号