CN108171142B - 一种确定复杂工业过程中关键变量因果关系的方法 - Google Patents
一种确定复杂工业过程中关键变量因果关系的方法 Download PDFInfo
- Publication number
- CN108171142B CN108171142B CN201711427123.3A CN201711427123A CN108171142B CN 108171142 B CN108171142 B CN 108171142B CN 201711427123 A CN201711427123 A CN 201711427123A CN 108171142 B CN108171142 B CN 108171142B
- Authority
- CN
- China
- Prior art keywords
- key variables
- causal relationship
- time
- variables
- key
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 230000001364 causal effect Effects 0.000 title claims abstract description 101
- 238000000034 method Methods 0.000 title claims abstract description 90
- 238000004519 manufacturing process Methods 0.000 title claims abstract description 66
- 238000013507 mapping Methods 0.000 claims abstract description 61
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 23
- 238000001514 detection method Methods 0.000 claims abstract description 11
- 230000001419 dependent effect Effects 0.000 claims abstract description 8
- 230000000694 effects Effects 0.000 claims abstract description 7
- 238000000342 Monte Carlo simulation Methods 0.000 claims abstract description 6
- 230000008569 process Effects 0.000 claims description 38
- 238000004517 catalytic hydrocracking Methods 0.000 claims description 20
- 241000208125 Nicotiana Species 0.000 claims description 12
- 235000002637 Nicotiana tabacum Nutrition 0.000 claims description 12
- 238000004458 analytical method Methods 0.000 claims description 9
- 230000002457 bidirectional effect Effects 0.000 claims description 8
- 238000004364 calculation method Methods 0.000 claims description 8
- UFHFLCQGNIYNRP-UHFFFAOYSA-N Hydrogen Chemical compound [H][H] UFHFLCQGNIYNRP-UHFFFAOYSA-N 0.000 claims description 5
- 229910052739 hydrogen Inorganic materials 0.000 claims description 5
- 239000001257 hydrogen Substances 0.000 claims description 5
- 238000005504 petroleum refining Methods 0.000 claims description 5
- 238000001035 drying Methods 0.000 claims description 4
- 238000007781 pre-processing Methods 0.000 claims description 4
- 230000009466 transformation Effects 0.000 claims description 4
- FFBHFFJDDLITSX-UHFFFAOYSA-N benzyl N-[2-hydroxy-4-(3-oxomorpholin-4-yl)phenyl]carbamate Chemical compound OC1=C(NC(=O)OCC2=CC=CC=C2)C=CC(=C1)N1CCOCC1=O FFBHFFJDDLITSX-UHFFFAOYSA-N 0.000 claims description 3
- 230000001934 delay Effects 0.000 claims description 3
- 238000012614 Monte-Carlo sampling Methods 0.000 claims description 2
- 230000008901 benefit Effects 0.000 abstract description 6
- 230000002159 abnormal effect Effects 0.000 description 12
- 230000000875 corresponding effect Effects 0.000 description 10
- 238000010586 diagram Methods 0.000 description 6
- 230000008878 coupling Effects 0.000 description 5
- 238000010168 coupling process Methods 0.000 description 5
- 238000005859 coupling reaction Methods 0.000 description 5
- 239000000243 solution Substances 0.000 description 5
- 238000006243 chemical reaction Methods 0.000 description 4
- 238000012360 testing method Methods 0.000 description 4
- 230000008859 change Effects 0.000 description 3
- 238000003745 diagnosis Methods 0.000 description 3
- 230000005284 excitation Effects 0.000 description 3
- 238000003723 Smelting Methods 0.000 description 2
- 229910000831 Steel Inorganic materials 0.000 description 2
- 238000005314 correlation function Methods 0.000 description 2
- 238000002347 injection Methods 0.000 description 2
- 239000007924 injection Substances 0.000 description 2
- 239000010959 steel Substances 0.000 description 2
- 238000012546 transfer Methods 0.000 description 2
- 241001598984 Bromius obscurus Species 0.000 description 1
- 230000005856 abnormality Effects 0.000 description 1
- 230000009471 action Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 239000003054 catalyst Substances 0.000 description 1
- 235000019504 cigarettes Nutrition 0.000 description 1
- 230000002596 correlated effect Effects 0.000 description 1
- 238000005336 cracking Methods 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000005194 fractionation Methods 0.000 description 1
- 239000000295 fuel oil Substances 0.000 description 1
- 238000005984 hydrogenation reaction Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 230000010354 integration Effects 0.000 description 1
- 238000006317 isomerization reaction Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 238000005312 nonlinear dynamic Methods 0.000 description 1
- 239000003921 oil Substances 0.000 description 1
- 230000000704 physical effect Effects 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 239000002994 raw material Substances 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 238000007670 refining Methods 0.000 description 1
- 238000012827 research and development Methods 0.000 description 1
- 230000004044 response Effects 0.000 description 1
- 238000013215 result calculation Methods 0.000 description 1
- 238000005096 rolling process Methods 0.000 description 1
- 238000005070 sampling Methods 0.000 description 1
- 229920006395 saturated elastomer Polymers 0.000 description 1
- 238000000926 separation method Methods 0.000 description 1
- 238000011144 upstream manufacturing Methods 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
- 239000002699 waste material Substances 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
- G06F2218/02—Preprocessing
- G06F2218/04—Denoising
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Data Mining & Analysis (AREA)
- General Engineering & Computer Science (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Physics (AREA)
- Computational Mathematics (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Algebra (AREA)
- Probability & Statistics with Applications (AREA)
- Operations Research (AREA)
- Bioinformatics & Computational Biology (AREA)
- Evolutionary Biology (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Artificial Intelligence (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Signal Processing (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明公开了一种确定复杂工业过程中关键变量因果关系的方法，通过对工业过程中待确定因果关系的关键变量的历史数据采用伪近邻思想计算每个关键变量的最佳时序嵌入维度；针对两两关键变量，假定因果关系，以假定因变量的最佳时序嵌入维度为标准，构造两个关键变量的时序重构流形，利用收敛交叉映射算法计算两者间的收敛交叉映射能力；基于蒙特卡洛模拟确定CCM的能力判定阈值，以此判定关键变量间假定因果关系的正确性，从而构建工业过程中关键变量的初步因果关系网络；利用时滞检测方法修正初步因果关系网络，得到最终关键变量因果关系网络。本发明充分利用生产离线数据，对生产过程无干扰作用，提高了安全性和经济效益。
Description
技术领域
本发明涉及工业过程异常工况诊断与溯因分析中的关键变量因果关系分析等技术领域，具体涉一种确定复杂工业过程中关键变量因果关系的方法。
背景技术
在石油炼制、钢铁冶炼等大型复杂工业***，流程中运行变量较多，各变量之间耦合严重，集成度高。例如，加氢裂化过程是石油炼制流程中一个子过程，是氢气在较高压力和温度下经催化剂作用使重质油发生加氢、裂化和异构化反应转化为轻质油的加工过程，主要包括加氢精制反应、加氢裂化反应、高低压分离以及分馏***这四个重要部分。加氢裂化的工业过程作为这类典型的庞大而复杂的流程工业，具有流程长、变量多、工况多以及非线性的特点。当流程中一个关键变量出现异常时，可能会引起相关变量出现异常，从而导致整个流程出现波动，然而操作人员根据经验得到的结果往往无法及时并准确地定位到异常发生的原因，使工业过程长时间处于非优状态，造成资源能源的浪费。
在异常工况诊断研究中，有关学者提出了基于定性模型、基于定量模型及基于数据驱动的方法。其中，基于定量模型和定性模型是从对流程物理特性的基本理解发展而来，依靠模型的先验知识，但它只适用于机理明确的工业过程。目前，在流程工业中存有海量的历史数据，因而基于数据驱动的诊断方法对比其他两类方法具有明显的优势，而通过这类方法分析异常传播路径、定位异常原因的核心是寻找变量之间的因果关系，但是由于工业过程中变量具有非线性、耦合严重等特性，变量之间即使存在相关性也不能保证变量之间存在因果关系，给流程中变量间的因果关系识别带来了挑战。
针对上述挑战，目前已有不少学者从不同角度出发提出了一些解决方法。格兰杰提出的Granger因果关系检验主要适合于变量之间存在线性相关性，同时还要求测试的两个变量相互独立，但是该方法对于耦合动态***不适用，虽然一些学者改进并提出了非线性Granger因果关系分析方法，但是该方法在不可分离非线性动态***中测试失败。互相关函数(Cross Correlation Function)通过计算变量之间最大CCF值得到时间延迟来确定信号的因果关系，但是这个方法适用于线性***，对于非线性***计算结果不能保证准确性。Schreiber提出的传递熵(Transfer Entropy)算法来判断因果关系虽然对非线性和线性***都能适用，但是不仅对参数改变敏感，而且计算量大、实现相对困难，同时结果计算的准确性与概率密度函数的估计密切相关，因此要求数据具有时间稳定性。此外，还有学者提出通过输入激励信号来监测模型的变化，但是这类方法对于实际工业流程来说，不仅难以设计激励信号，而且加入的激励信号本身会影响工业流程的正常运行。
针对上述方法的局限性以及复杂工业过程的实际特点，亟需提出一种用于确定复杂工业过程中关键变量间因果关系的方法。考虑Sugihara等提出的收敛交叉映射(CCM)方法不仅适用于非线性、中强耦合的复杂***，而且计算方便参数选择量少的特点，本专利采用收敛交叉映射(CCM)的思想确定关键变量间的因果关系。然而由于CCM初始应用于生态***，并没有对时序嵌入维度的选择进行讨论，而且对于具有高维空间的复杂工业过程，依据人工经验选择嵌入维度和CCM能力判定阈值无法确保结果的准确性，并且很难实现自动监测。因此本专利通过改进CCM方法确定复杂工业过程中关键变量间因果关系的方法，最终形成关键变量的因果关系网。
发明内容
本发明的目的在于提供一种计算方便，参数选择量少，充分考虑了关键变量之间的耦合问题，能解决具有中强耦合的非线性复杂动态***中关键变量间的因果关系识别，为工业过程异常工况的溯因分析的方法。
为了实现上述发明目的，本发明提供技术方案如下：一种确定复杂工业过程中关键变量因果关系的方法，所述方法包含以下步骤：
步骤一：收集工业过程中关键变量的生产历史数据，对样本进行数据预处理：
选取待确定因果关系的n个关键变量，收集这些关键变量一段连续的时间序列的生产历史数据组成一个变量维度为n，样本长度为L的数据集作为样本；然后利用小波变换对所述的样本进行去噪处理；所述的n为大于1的自然数；
步骤二：确定重构流形的最佳时序嵌入维度E：
对步骤一预处理后的样本中每个关键变量，采用伪近方法计算其最佳时序嵌入维度E；即对所述的样本中的关键变量的第i个重构流形，设定嵌入维度E＝k的最近邻点的距离为Di(k)，嵌入维度为E＝k+1的最近邻点的距离为Di(k+1)，计算两者的几何平均，通过对数形式比较两者的值，使其结果具有鲁棒性；则有
其中，L表示时间序列的样本长度，τ为选取的时间间隔；当流形中的最近邻点的距离随着嵌入维度的增加不变时，即Di(k)无限接近于Di(k+1)，有G(k)的值无限接近于0，称为饱和状态；当G(k)的值不变时，选取最先进入饱和状态时刻对应的嵌入维度E为重构流形的最佳时序嵌入维度；
步骤三：计算两两关键变量间的收敛交叉映射能力：
对于两两关键变量，设定因果关系，并以设定因变量的最佳时序嵌入维度E为标准，构造两变量的重构流形，然后采用k近邻算法计算得到最近重构流形坐标，利用CCM算法计算两者间的收敛交叉映射能力；
步骤四：基于蒙特卡洛模拟方法确定收敛交叉映射能力的判定阈值，构建初步因果关系网络：
采用蒙特卡洛抽样，选取多对不相关的随机时间序列，分别计算它们在不同样本长度L下收敛交叉映射能力的均值和标准差；然后，曲线拟合得到与样本长度L相关的均值函数和标准差函数σr＝1.17L-0.46；接着，根据确定的样本长度L在拟合函数中找到对应的均值与标准差；以此按照3σ准则得到收敛交叉映射能力的判定阈值r0＝-3.8L-1.08+3.51L-0.46；若步骤三中计算的某两关键变量收敛交叉映射能力的收敛值超出该判定阈值，则认为设定的关键变量因果关系成立，遍历所有关键变量，判定所有关键变量两两间的因果关系，构建初步因果关系网络。
进一步地，还包括步骤五，采用时滞检测方法修正初步因果关系网络，得到最终的工业过程中关键变量因果关系网络：
首先，从步骤四中得到的初步因果关系网络中挑选出双向因果关系的关键变量：根据工业过程的工艺确定关键变量之间存在的时间延迟范围，设为Δt1，Δt2，其中Δt1≤0，Δt2≥0，设定时间序列{Xi}和{Xj}，利用CCM算法通过Xi(t)来反映Xj(t-Δt)，其中Δt∈[Δt1，Δt2]，计算不同时滞下的收敛交叉映射能力；比较不同时滞下的收敛交叉映射能力，得到的最大收敛交叉映射能力为最优结果，此最优结果对应的时滞认为是在流程中两关键变量之间的时间延迟；如果最优收敛交叉映射能力对应的时滞为正，则表示时间序列{Xi}的未来值能准确的反映Xj的过去值，因此存在从Xj到Xi方向上的因果关系，反之，时滞为负则表示不存在；
然后，在消除了假象双向因果关系的网络中挑选出存在相同原因且两者之间存在因果关系的关键变量，如果变量X对变量Z的影响先于X对Y的影响，此时表示Z的过去值能准确的预测Y的未来值，此时判断最优收敛交叉映射能力对应的时滞为负；从步骤四中得到的初步因果关系网络中经过以上分析后，对步骤四中得到的初步因果关系网络进行修正，得到最终的工业过程中关键变量因果关系网络。
进一步地，步骤三中所述的“利用CCM算法计算关键变量两两间的收敛交叉映射能力”的过程为：
对于两两含有时序的关键变量Xi和Xj，设定Xj到Xi存在因果关系；以假定因变量Xj的最佳时序嵌入维度为标准，构造关键变量Xi和Xj的重构流形；对于长度为L的样本时间序列{Xi}和{Xj}，其中i，j＝1，2，3，......，n且i≠j，则时间序列在t时刻的重构流形坐标分别为：
Mi(t)＝[xi(t)，xi(t-τ)，......，xi(t-(E-1)τ)]
Mj(t)＝[xj(t)，xj(t-τ)，......，xj(t-(E-1)τ)]
由上式可得时间序列{Xi}和{Xj}的重构流形分别为Mi＝{Mi(t)}，Mj＝{Mj(t)}；从流形Mi中找到距离Mi(t)最近的E+1个点Mi，t＝{Mi(t1)，Mi(t2)，Mi(t3)，......，Mi(tE+1)}，采用k近邻算法来寻找这E+1个点；根据时间序列{Xi}和{Xj}存在从Xj到Xi的因果关系的设定，通过Mi，t来估计时间序列{Xj}的值，得到估计的时间序列其具体计算公式如下：
计算估计时间序列与原时间序列{Xj}的相关系数得到收敛交叉映射能力ri→j如下：
其中，表示时间序列{Xj}的均值，表示估计时间序列的均值；根据CCM算法原理可得，如果存在从Xj到Xi的因果关系，则随着样本时间序列长度L的增加，相关系数ri→j将收敛于一个大于0的值，此时，称该相关系数为收敛交叉映射能力；如果收敛交叉映射能力收敛于大于0的值，则认为两个关键变量之间存在因果关系，否则不存在。
进一步地，步骤一中所述的“一段连续的时间序列”为无人工调整的一段时间。
进一步地，所述的L为大于或等于200的自然数。
进一步地，所述的工业过程为石油炼制的加氢裂化过程，所述的关键变量包括出口温度、入口温度、一床层顶温、二床层顶温、三床层顶温、冷氢注入一二床层量、冷氢注入二三床层量。
进一步地，所述的工业过程为烟草制丝的烘丝过程，所述的关键变量包括入口水分、热风风量、一区筒壁压力、二区筒壁压力、出口水分这五个关键变量。
与现有技术相比，本发明的有益效果在于：一是通过改进CCM的方法，能解决非线性中强耦合的复杂动态***中关键变量间的因果关系识别问题，为分析复杂工业过程中异常工况的溯因分析提供重要条件。二是采用伪近邻的思想确定CCM算法中的重要调优参数(时序嵌入维度)，解决了由于人工选择工业过程中多个关键变量的调优参数而导致算法准确性降低的问题。三是基于蒙特卡洛模拟方法确定CCM的判定阈值，减少因果分析中的人为判断。四是通过时滞检测消除假象双因果关系以及部分同源因果造成的伪因果关系，降低了工业过程关键变量因果关系网络的干扰项。依据本专利构建的工业过程中关键变量因果关系网络，对比现有方法，适用于复杂的工业实际过程，计算更加简便，充分利用工业过程的离线数据，且对工业过程无干扰作用，提高了安全性和经济效益，对回溯异常传播路径、定位异常原因效果良好。
附图说明
图1为本发明的一个实施例的算法流程图。
图2为本发明的一个实施例加氢裂化流程精制反应装置图。
图3为本发明的一个实施例嵌入维度曲线图。
图4为本发明的一个实施例交叉映射能力曲线。
图5为本发明的一个实施例初步因果关系网络图。
图6为本发明的一个实施例最终因果关系网络图。
图7为本发明的另一个实施例的嵌入维度曲线图。
图8为本发明的另一个实施例的交叉映射能力曲线.。
图9为本发明的另一个实施例的初步因果关系网络。
图10为本发明的另一个实施例的最终因果关系网络。
具体实施方式
为了更好地说明本发明的研发目的、技术方案及优势，以及充分公开本发明技术方案，下面对本发明的技术方案做进一步说明。
本发明提供的技术适用于石油炼制、烟草制丝、钢铁冶炼等复杂工业过程中关键变量间因果关系的确定。
一种确定复杂工业过程中关键变量因果关系的方法，包含以下步骤：
步骤一：收集工业过程关键变量的生产历史数据，对样本进行数据预处理：
选取待确定因果关系的n(n大于等于1)个关键变量，收集这些关键变量一段连续的时间序列的生产历史数据组成一个变量维度为n，样本长度为L的数据集作为样本；然后利用小波变换对所述的样本进行去噪处理。
为了使选取的样本更具稳定性，生产历史数选取相对稳定的、无人工调整的一段时间产生的数据。
步骤二：确定重构流形的最佳时序嵌入维度E：
对步骤一预处理后的样本中每个关键变量，采用伪近方法计算其最佳时序嵌入维度E；即对所述的样本中的关键变量的第i个重构流形，设定嵌入维度E＝k的最近邻点的距离为Di(k)，嵌入维度为E＝k+1的最近邻点的距离为Di(k+1)，计算两者的几何平均，通过对数形式比较两者的值，使其结果具有鲁棒性。则有
其中，L表示时间序列的样本长度，τ为选取的时间间隔；当流形中的最近邻点的距离随着嵌入维度的增加不变时，即Di(k)无限接近于Di(k+1)，有G(k)的值无限接近于0，称为饱和状态。当G(k)的值不变时，选取最先进入饱和状态时刻对应的嵌入维度E为重构流形的最佳时序嵌入维度。
步骤三：计算两两关键变量间的收敛交叉映射能力：
对于两两关键变量，设定因果关系，并以设定因变量的最佳时序嵌入维度E为标准，构造两变量的重构流形，然后采用k近邻算法计算得到最近重构流形坐标，利用CCM算法计算两者间的收敛交叉映射能力；
对于两两含有时序的关键变量Xi和Xj，设定Xj到Xi存在因果关系；以假定因变量Xj的最佳时序嵌入维度为标准，构造关键变量Xi和Xj的重构流形；对于长度为L的样本时间序列{Xi}和{Xj}，其中i，j＝1，2，3，......，n且i≠j，则时间序列在t时刻的重构流形坐标分别为：
Mi(t)＝[xi(t)，xi(t-τ)，......，xi(t-(E-1)τ)]
Mj(t)＝[xj(t)，xj(t-τ)，......，xj(t-(E-1)τ)]
由上式可得时间序列{Xi}和{Xj}的重构流形分别为Mi＝{Mi(t)}，Mj＝{Mj(t)}；从流形Mi中找到距离Mi(t)最近的E+1个点Mi，t＝{Mi(t1)，Mi(t2)，Mi(t3)，......，Mi(tE+1)}，采用k近邻算法来寻找这E+1个点；根据时间序列{Xi}和{Xj}存在从Xj到Xi的因果关系的设定，通过Mi，t来估计时间序列{Xj}的值，得到估计的时间序列其具体计算公式如下：
计算估计时间序列与原时间序列{Xj}的相关系数得到收敛交叉映射能力ri→j如下：
其中，表示时间序列{Xj}的均值，表示估计时间序列的均值。根据CCM算法原理可得，如果存在从Xj到Xi的因果关系，则随着样本时间序列长度L的增加，相关系数ri→j将收敛于一个大于0的值，此时，称该相关系数为收敛交叉映射能力；如果收敛交叉映射能力收敛于大于0的值，则认为两个关键变量之间存在因果关系，否则不存在。
步骤四：基于蒙特卡洛模拟方法确定收敛交叉映射能力的判定阈值，构建初步因果关系网络：
采用蒙特卡洛抽样，选取多对不相关的随机时间序列，分别计算它们在不同样本长度L下收敛交叉映射能力的均值和标准差；然后，曲线拟合得到与样本长度L相关的均值函数和标准差函数σr＝1.17L-0.46；接着，根据确定的样本长度L在拟合函数中找到对应的均值与标准差；以此按照3σ准则得到收敛交叉映射能力的判定阈值r0＝-3.8L-1.08+3.51L-0.46；若步骤三中计算的某两关键变量收敛交叉映射能力的收敛值超出该判定阈值，则认为设定的关键变量因果关系成立，遍历所有关键变量，判定所有关键变量两两间的因果关系，构建初步因果关系网络。
为了消除假象双因果关系以及部分同源因果造成的伪因果关系，降低了工业过程关键变量因果关系网络的干扰项，对前述技术方案的一种改进为：还包括：
步骤五，采用时滞检测方法修正初步因果关系网络，得到最终的工业过程中关键变量因果关系网络：
首先，从步骤四中得到的初步因果关系网络中挑选出双向因果关系的关键变量：根据工业过程的工艺确定关键变量之间存在的时间延迟范围，设为Δt1，Δt2，其中Δt1≤0，Δt2≥0，设定时间序列Xi和Xj，利用CCM算法通过Xi(t)来反映Xj(t-Δt)，其中Δt∈[Δt1，Δt2]，计算不同时滞下的收敛交叉映射能力；比较不同时滞下的收敛交叉映射能力，得到的最大收敛交叉映射能力为最优结果，此最优结果对应的时滞认为是在流程中两关键变量之间的时间延迟；如果最优收敛交叉映射能力对应的时滞为正，则表示时间序列{Xi}的未来值能准确的反映Xj的过去值，因此存在从Xj到Xi方向上的因果关系，反之，时滞为负则表示不存在；
然后，在消除了假象双向因果关系的网络中挑选出存在相同原因且两者之间存在因果关系的关键变量，如果变量X对变量Z的影响先于X对Y的影响，此时表示Z的过去值能准确的预测Y的未来值，此时判断最优收敛交叉映射能力对应的时滞为负；从步骤四中得到的初步因果关系网络中经过以上分析后，对步骤四中得到的初步因果关系网络进行修正，得到最终的工业过程中关键变量因果关系网络。
具体实施例1：
下面将结合我国某炼油厂加氢裂化过程的实际运行数据举例说明，针对本发明中提出的基于改进收敛交叉映射算法确定加氢裂化过程中关键变量因果关系，根据图1中流程图执行，并做出详细说明。应该强调的是，下述说明仅仅是示例性的，且下述中实例只是加氢裂化流程中一部分关键变量，而不是为了限制本发明的范围及其应用。所提供方法包括如下步骤：
步骤一，收集加氢裂化流程关键变量的生产历史数据，进行样本数据预处理。即，对于待确定因果关系的加氢裂化流程加氢精制反应部分7个关键变量，其中加氢精制反应装置简化示意图如图2所示，从历史数据库中选取同一工况下非人工干预的连续时间序列来组成数据集长度为1000的样本，并利用小波变换对样本进行去噪处理。表1为待测的加氢裂化流程中部分关键变量的描述。
表1关键变量描述
步骤二，对步骤一得到的预处理后的每个关键变量，采用伪近邻思想得到重构流形的最佳时序嵌入维度E。由于提取加氢裂化流程关键变量的时间序列本身是离散的点，采样间隔为每分钟，因此选择时间间隔为τ＝1。计算得到各关键变量的G(k)图形如图3所示。根据G(k)曲线，当其变化速率小于0.02时，认为G(k)的值趋于饱和得到最佳嵌入维度。表2为通过伪近邻计算方法得到各关键变量的重构流形的最佳嵌入维度。
表2关键变量最佳嵌入维度
步骤三，计算两两关键变量间的收敛交叉映射能力。对于两两关键变量，假定因果关系，并以假定因变量的最佳嵌入维度E为标准，构造两变量的重构流形，得到重构流形，采用k近邻算法得到最近重构流形坐标。利用CCM算法计算关键变量两两间的收敛交叉映射能力，在不同样本长度下计算收敛交叉映射能力得到的收敛曲线，然后对比曲线是否收敛得到收敛结果，图4给出两两关键变量间收敛的CCM结果曲线。需注意的是，由于计算的两关键变量之间可能存在嵌入维度不一致的情况，此时选择原因变量的最佳时序嵌入维度来重构流形。表3为关键变量两两之间计算得到的收敛交叉映射能力，其中行代表原因变量，列代表结果变量，不同关键变量之间的“/”表示计算的CCM值不收敛。
表3收敛交叉能力
步骤四，基于蒙特卡洛模拟方法确定收敛交叉映射能力的判定阈值，构建初步因果关系网络。采用假设检验的方式，使用多对不相关的随机时间序列，分别计算在不同样本长度下收敛交叉映射能力的均值和标准差，曲线拟合得到与样本长度相关的均值和标准差函数，并根据3σ准则得到收敛交叉映射能力的阈值r0＝-3.8L-1.08+3.51L-0.46。在实例中，样本长度L＝1000，因而得到收敛交叉映射能力的阈值为r0＝0.1441，判断表3中的数据，得到加氢裂化流程关键变量的初步因果关系网络如图5所示。
步骤五，采用时滞检测修正初步因果关系网络。在加氢裂化流程中，前后变量之间存在时间延迟，如果上游变量的当前值改变，经过一段时间后才能影响到下游变量的值。因此可通过时滞检测区分出由于强单向因果关系引起的同步现象以及部分同源因果引起的伪因果关系。根据加氢裂化流程工艺大致确定两关键变量之间可能存在的时间延迟范围，从初步因果关系网络中可得到存在双向因果关系的关键变量，经过时滞检测消除假象双向因果关系后得到进一步的因果关系网络，再从中挑选出存在相同原因且两者之间存在因果关系的关键变量，经过时滞检测后消去变量之间的伪因果关系。关键变量之间的时滞正负结果如表4所示。依据时滞检测结果修正初步因果关系网络后得到最终的因果关系网络如图6所示。
表4关键变量时滞检测结果
经过验证，得到的最终因果关系网络与加氢裂化流程工艺了解的经验一致，符合实际情况。该结果与人工通过流程经验知识判别因果关系相比，不仅具有快速性而且降低了操作人员的工作量。本专利完成的加氢裂化流程因果关系网络，对比现有方法，计算更加简便，充分利用加氢裂化流程的离线数据，且对加氢裂化流程无干扰作用，提高了安全性和经济效益，进而为回溯异常传播路径、定位异常原因打下良好基础。
具体实施例2：
再以烟草制丝生产为例，烟草制丝生产是典型的间歇生产过程，烘丝是最关键的一道工序，它的目标是把烟叶等原材料加工成符合卷烟质量标准以及下一阶段卷制生产要求的叶丝。该生产过程较为复杂，且叶丝受到水分、温度等方面因素的影响，因此，以烘丝过程为例，分析该过程中关键变量之间的因果关系，主要步骤如下：
步骤一，选取烘丝过程中入口水分、热风风量、一区筒壁压力、二区筒壁压力、出口水分这五个关键变量，样本的长度为300；
步骤二，确定重构流形的最佳时序嵌入维度E，各关键变量的G(k)图形如图7所示，因而得到各关键变量的最佳时序嵌入维度分别为7、6、6、6、5；
步骤三，计算得到关键变量两两间的收敛交叉映射能力，判断是否收敛，得到收敛的曲线图如图8所示，表5为关键变量两两之间计算得到的收敛交叉映射能力，其中行代表原因变量，列代表结果变量，不同关键变量之间的“/”表示计算的CCM值不收敛。
表5收敛交叉能力
步骤四，判定阈值判定得到初步因果关系网络图如图9所示；
步骤五，进行时滞检测结果，得到最终因果关系网络如图10所示，结果符合实际情况。
以上所述仅为本发明的实施例，并非因此限制本发明的专利范围，凡是利用本发明说明书内容所作的等效结构或等效流程变换，或直接或间接运用在其它相关的技术领域，均同理包括在本发明的保护范围内。
Claims (7)
1.一种确定复杂工业过程中关键变量因果关系的方法，其特征在于，所述方法包含以下步骤：
步骤一：收集工业过程中关键变量的生产历史数据，对样本进行数据预处理：
选取待确定因果关系的n个关键变量，收集这些关键变量一段连续的时间序列的生产历史数据组成一个变量维度为n，样本长度为L的数据集作为样本；然后利用小波变换对所述的样本进行去噪处理；所述的n为大于1的自然数；
步骤二：确定重构流形的最佳时序嵌入维度E：
对步骤一预处理后的样本中每个关键变量，采用伪近方法计算其最佳时序嵌入维度E；即对所述的样本中的关键变量的第i个重构流形，设定嵌入维度E＝k的最近邻点的距离为Di(k)，嵌入维度为E＝k+1的最近邻点的距离为Di(k+1)，计算两者的几何平均，通过对数形式比较两者的值，使其结果具有鲁棒性；则有
其中，L表示时间序列的样本长度，τ为选取的时间间隔；当流形中的最近邻点的距离随着嵌入维度的增加不变时，即Di(k)无限接近于Di(k+1)，有G(k)的值无限接近于0，称为饱和状态；当G(k)的值不变时，选取最先进入饱和状态时刻对应的嵌入维度E为重构流形的最佳时序嵌入维度；
步骤三：计算两两关键变量间的收敛交叉映射能力：
对于两两关键变量，设定因果关系，并以设定因变量的最佳时序嵌入维度E为标准，构造两变量的重构流形，然后采用k近邻算法计算得到最近重构流形坐标，利用CCM算法计算两者间的收敛交叉映射能力；
步骤四：基于蒙特卡洛模拟方法确定收敛交叉映射能力的判定阈值，构建初步因果关系网络：
采用蒙特卡洛抽样，选取多对不相关的随机时间序列，分别计算它们在不同样本长度L下收敛交叉映射能力的均值和标准差；然后，曲线拟合得到与样本长度L相关的均值函数和标准差函数σr＝1.17L-0.46；接着，根据确定的样本长度L在拟合函数中找到对应的均值与标准差；以此按照3σ准则得到收敛交叉映射能力的判定阈值r0＝-3.8L-1.08+3.51L-0.46；若步骤三中计算的某两关键变量收敛交叉映射能力的收敛值超出该判定阈值，则认为设定的关键变量因果关系成立，遍历所有关键变量，判定所有关键变量两两间的因果关系，构建初步因果关系网络。
2.根据权利要求1所述的确定复杂工业过程中关键变量因果关系的方法，其特征在于，还包括步骤五，采用时滞检测方法修正初步因果关系网络，得到最终的工业过程中关键变量因果关系网络：
首先，从步骤四中得到的初步因果关系网络中挑选出双向因果关系的关键变量：根据工业过程的工艺确定关键变量之间存在的时间延迟范围，设为Δt1，Δt2，其中Δt1≤0，Δt2≥0，设定时间序列{Xi}和{Xj}，利用CCM算法通过Xi(t)来反映Xj(t-Δt)，其中Δt∈[Δt1，Δt2]，计算不同时滞下的收敛交叉映射能力；比较不同时滞下的收敛交叉映射能力，得到的最大收敛交叉映射能力为最优结果，此最优结果对应的时滞认为是在流程中两关键变量之间的时间延迟；如果最优收敛交叉映射能力对应的时滞为正，则表示时间序列{Xi}的未来值能准确的反映Xj的过去值，因此存在从Xj到Xi方向上的因果关系，反之，时滞为负则表示不存在；
然后，在消除了假象双向因果关系的网络中挑选出存在相同原因且两者之间存在因果关系的关键变量，如果变量X对变量Z的影响先于X对Y的影响，此时表示Z的过去值能准确的预测Y的未来值，此时判断最优收敛交叉映射能力对应的时滞为负；从步骤四中得到的初步因果关系网络中经过以上分析后，对步骤四中得到的初步因果关系网络进行修正，得到最终的工业过程中关键变量因果关系网络。
3.根据权利要求1所述的确定复杂工业过程中关键变量因果关系的方法，其特征在于，步骤三中所述的“利用CCM算法计算关键变量两两间的收敛交叉映射能力”的过程为：
对于两两含有时序的关键变量Xi和Xj，设定Xj到Xi存在因果关系；以假定因变量Xj的最佳时序嵌入维度为标准，构造关键变量Xi和Xj的重构流形；对于长度为L的样本时间序列{Xi}和{Xj}，其中i，j＝1，2，3，......，n且i≠j，则时间序列在t时刻的重构流形坐标分别为：
Mi(t)＝[xi(t)，xi(t-τ)，......，xi(t-(E-1)τ)]
Mj(t)＝[xj(t)，xj(t-τ)，......，xj(t-(E-1)τ)]
由上式可得时间序列{Xi}和{Xj}的重构流形分别为Mi＝{Mi(t)}，Mj＝{Mj(t)}；从流形Mi中找到距离Mi(t)最近的E+1个点Mi，t＝{Mi(t1)，Mi(t2)，Mi(t3)，......，Mi(tE+1)}，采用k近邻算法来寻找这E+1个点；根据时间序列{Xi}和{Xj}存在从Xj到Xi的因果关系的设定，通过Mi，t来估计时间序列{Xj}的值，得到估计的时间序列其具体计算公式如下：
计算估计时间序列与原时间序列{Xj}的相关系数得到收敛交叉映射能力ri→j如下：
其中，表示时间序列{Xj}的均值，表示估计时间序列的均值；根据CCM算法原理可得，如果存在从Xj到Xi的因果关系，则随着样本时间序列长度L的增加，相关系数ri→j将收敛于一个大于0的值，此时，称该相关系数为收敛交叉映射能力；如果收敛交叉映射能力收敛于大于0的值，则认为两个关键变量之间存在因果关系，否则不存在。
4.根据权利要求1所述的确定复杂工业过程中关键变量因果关系的方法，其特征在于，步骤一中所述的“一段连续的时间序列”为无人工调整的一段时间。
5.根据权利要求1所述的确定复杂工业过程中关键变量因果关系的方法，其特征在于，所述的L为大于或等于200的自然数。
6.根据权利要求1所述的确定复杂工业过程中关键变量因果关系的方法，其特征在于，所述的工业过程为石油炼制的加氢裂化过程，所述的关键变量包括出口温度、入口温度、一床层顶温、二床层顶温、三床层顶温、冷氢注入一二床层量、冷氢注入二三床层量。
7.根据权利要求1所述的确定复杂工业过程中关键变量因果关系的方法，其特征在于，所述的工业过程为烟草制丝的烘丝过程，所述的关键变量包括入口水分、热风风量、一区筒壁压力、二区筒壁压力、出口水分这五个关键变量。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201711427123.3A CN108171142B (zh) | 2017-12-26 | 2017-12-26 | 一种确定复杂工业过程中关键变量因果关系的方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201711427123.3A CN108171142B (zh) | 2017-12-26 | 2017-12-26 | 一种确定复杂工业过程中关键变量因果关系的方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN108171142A CN108171142A (zh) | 2018-06-15 |
| CN108171142B true CN108171142B (zh) | 2019-02-12 |
Family
ID=62520499
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201711427123.3A Active CN108171142B (zh) | 2017-12-26 | 2017-12-26 | 一种确定复杂工业过程中关键变量因果关系的方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN108171142B (zh) |
Families Citing this family (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| TWI662424B (zh) * | 2017-12-25 | 2019-06-11 | 帆宣系統科技股份有限公司 | 領先輔助參數的選擇方法以及結合關鍵參數及領先輔助參數進行設備維護預診斷的方法 |
| CN109886409B (zh) * | 2019-02-15 | 2023-12-12 | 南京信息工程大学 | 一种多维时间序列的定量因果关系判定方法 |
| CN109814513B (zh) * | 2019-03-20 | 2020-12-29 | 广东辛孚科技有限公司 | 一种基于数据模型的催化裂化装置优化方法 |
| CN109816940B (zh) * | 2019-03-21 | 2023-05-09 | 北京天诚同创电气有限公司 | 污水处理厂的故障报警方法和装置 |
| CN111445674B (zh) * | 2020-04-08 | 2021-09-14 | 浙江浙能技术研究院有限公司 | 一种面向百万千瓦超超临界机组制粉***报警管理的因果网络构建方法 |
| CN111626099B (zh) * | 2020-04-10 | 2023-06-09 | 浙江大学 | 一种基于改进ccm的工业控制***多回路振荡因果关系分析方法 |
| CN112037106B (zh) * | 2020-08-07 | 2023-12-15 | 汉威科技集团股份有限公司 | 一种基于特征互相关性和概率密度的数据异常分析方法 |
| CN112101480B (zh) * | 2020-09-27 | 2022-08-05 | 西安交通大学 | 一种多变量聚类与融合的时间序列组合预测方法 |
| CN112765011B (zh) * | 2020-12-30 | 2023-10-10 | 上海昆涞生物科技有限公司 | 质控状态判定方法、装置及电子设备 |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105205224A (zh) * | 2015-08-28 | 2015-12-30 | 江南大学 | 基于模糊曲线分析的时间差高斯过程回归软测量建模方法 |
| CN106156434A (zh) * | 2016-07-11 | 2016-11-23 | 江南大学 | 基于局部时滞重构的滑动窗时间差‑高斯过程回归建模方法 |
Family Cites Families (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103020591B (zh) * | 2012-11-21 | 2015-09-16 | 燕山大学 | 一种基于因果网络分析的中等规模人群异常行为检测方法 |
| JP6354192B2 (ja) * | 2014-02-14 | 2018-07-11 | オムロン株式会社 | 因果ネットワーク生成システム |
| CN104504637A (zh) * | 2014-12-11 | 2015-04-08 | 广东工业大学 | 基于行为时间序列的社交网络因果关系发现算法 |
| CN105023044B (zh) * | 2015-07-21 | 2017-10-24 | 清华大学 | 基于大量时间序列的交通流因果关系挖掘方法 |
-
2017
- 2017-12-26 CN CN201711427123.3A patent/CN108171142B/zh active Active
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105205224A (zh) * | 2015-08-28 | 2015-12-30 | 江南大学 | 基于模糊曲线分析的时间差高斯过程回归软测量建模方法 |
| CN106156434A (zh) * | 2016-07-11 | 2016-11-23 | 江南大学 | 基于局部时滞重构的滑动窗时间差‑高斯过程回归建模方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN108171142A (zh) | 2018-06-15 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN108171142B (zh) | 一种确定复杂工业过程中关键变量因果关系的方法 | |
| CN112101480B (zh) | 一种多变量聚类与融合的时间序列组合预测方法 | |
| CN109992921B (zh) | 一种燃煤电厂锅炉热效率的在线软测量方法及*** | |
| CN104699077B (zh) | 一种基于嵌套迭代费舍尔判别分析的故障变量隔离方法 | |
| CN101520652B (zh) | 一种数控装备服役可靠性的评估方法 | |
| CN111046341A (zh) | 一种基于主成分分析的非常规天然气压裂效果评价及产能预测方法 | |
| CN105225979A (zh) | 一种半导体器件制程预测***和方法 | |
| CN112651166B (zh) | 脱硝***入口氮氧化物浓度预测方法、装置及脱硝*** | |
| CN105930629B (zh) | 一种基于海量运行数据的在线故障诊断方法 | |
| CN111275255A (zh) | 一种混凝土坝变形监测预报模型的构建方法 | |
| CN104865827B (zh) | 一种基于多工况模型的抽油机采油优化方法 | |
| CN106773693A (zh) | 一种工业控制多回路振荡行为稀疏因果分析方法 | |
| CN106156434A (zh) | 基于局部时滞重构的滑动窗时间差‑高斯过程回归建模方法 | |
| Zhong et al. | Multimode non‐Gaussian process monitoring based on local entropy independent component analysis | |
| CN103488561A (zh) | 一种在线升级主样本模型的kNN故障检测方法 | |
| CN108628288B (zh) | 一种用于一阶时滞组合积分控制***的性能评价方法 | |
| WO2019071438A1 (zh) | 对流程行业中的设备进行状态监测的方法、装置及介质 | |
| CN103926919B (zh) | 基于小波变换和Lasso函数的工业过程故障检测方法 | |
| CN117556366B (zh) | 基于数据筛选的数据异常检测***及方法 | |
| CN106647274B (zh) | 一种连续生产过程中运行工况稳态判别方法 | |
| CN102750445B (zh) | 基于复化Simpson公式改进多变量灰色模型的故障预测方法 | |
| CN103995985A (zh) | 基于Daubechies小波变换和弹性网的故障检测方法 | |
| CN108345214B (zh) | 一种基于替代数据法的工业过程非线性检测方法 | |
| CN107220705A (zh) | 常减压装置常压塔顶干点预测方法 | |
| CN102914970B (zh) | 工业过程控制***中不能直接测量的性能参数的在线数据驱动估计方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |