CN113298170B - 非稳定工况下能耗模式识别方法及系统 - Google Patents
非稳定工况下能耗模式识别方法及系统 Download PDFInfo
- Publication number
- CN113298170B CN113298170B CN202110613823.1A CN202110613823A CN113298170B CN 113298170 B CN113298170 B CN 113298170B CN 202110613823 A CN202110613823 A CN 202110613823A CN 113298170 B CN113298170 B CN 113298170B
- Authority
- CN
- China
- Prior art keywords
- stationary
- sequence
- energy consumption
- time
- sample
- 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
- 238000005265 energy consumption Methods 0.000 title claims abstract description 106
- 238000000034 method Methods 0.000 title claims abstract description 76
- 238000012544 monitoring process Methods 0.000 claims abstract description 52
- 230000008569 process Effects 0.000 claims abstract description 46
- 230000005654 stationary process Effects 0.000 claims abstract description 46
- 239000011159 matrix material Substances 0.000 claims abstract description 39
- 238000009776 industrial production Methods 0.000 claims abstract description 29
- 238000005259 measurement Methods 0.000 claims abstract description 29
- 238000012549 training Methods 0.000 claims abstract description 29
- 238000005520 cutting process Methods 0.000 claims abstract description 19
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 18
- 230000007704 transition Effects 0.000 claims abstract description 16
- 230000033228 biological regulation Effects 0.000 claims abstract description 13
- 238000007781 pre-processing Methods 0.000 claims abstract description 10
- 238000004519 manufacturing process Methods 0.000 claims abstract description 6
- 239000013598 vector Substances 0.000 claims description 25
- 238000010586 diagram Methods 0.000 claims description 18
- 238000003909 pattern recognition Methods 0.000 claims description 7
- 238000012545 processing Methods 0.000 claims description 6
- 238000004088 simulation Methods 0.000 claims description 6
- 238000003860 storage Methods 0.000 claims description 6
- 230000009471 action Effects 0.000 claims description 5
- 238000012952 Resampling Methods 0.000 claims description 4
- 238000010276 construction Methods 0.000 claims description 4
- 238000000605 extraction Methods 0.000 claims description 3
- 238000000513 principal component analysis Methods 0.000 claims description 3
- 238000007637 random forest analysis Methods 0.000 claims description 3
- 230000002776 aggregation Effects 0.000 claims description 2
- 238000004220 aggregation Methods 0.000 claims description 2
- 238000005354 coacervation Methods 0.000 claims 1
- 238000004458 analytical method Methods 0.000 description 10
- 238000004364 calculation method Methods 0.000 description 7
- 230000006641 stabilisation Effects 0.000 description 5
- 238000011105 stabilization Methods 0.000 description 5
- 230000008859 change Effects 0.000 description 3
- 238000000354 decomposition reaction Methods 0.000 description 3
- 238000009826 distribution Methods 0.000 description 3
- 230000006870 function Effects 0.000 description 3
- CURLTUGMZLYLDI-UHFFFAOYSA-N Carbon dioxide Chemical compound O=C=O CURLTUGMZLYLDI-UHFFFAOYSA-N 0.000 description 2
- 229910000831 Steel Inorganic materials 0.000 description 2
- 230000008901 benefit Effects 0.000 description 2
- 239000004568 cement Substances 0.000 description 2
- 238000012512 characterization method Methods 0.000 description 2
- 238000010924 continuous production Methods 0.000 description 2
- 238000007405 data analysis Methods 0.000 description 2
- 238000013507 mapping Methods 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 239000002994 raw material Substances 0.000 description 2
- 230000009467 reduction Effects 0.000 description 2
- 238000005070 sampling Methods 0.000 description 2
- 238000000638 solvent extraction Methods 0.000 description 2
- 239000010959 steel Substances 0.000 description 2
- 239000000126 substance Substances 0.000 description 2
- 230000002159 abnormal effect Effects 0.000 description 1
- 230000005856 abnormality Effects 0.000 description 1
- 229910002092 carbon dioxide Inorganic materials 0.000 description 1
- 239000001569 carbon dioxide Substances 0.000 description 1
- 230000008878 coupling Effects 0.000 description 1
- 238000010168 coupling process Methods 0.000 description 1
- 238000005859 coupling reaction Methods 0.000 description 1
- 238000013075 data extraction Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000018109 developmental process Effects 0.000 description 1
- 238000004134 energy conservation Methods 0.000 description 1
- 230000007613 environmental effect Effects 0.000 description 1
- 239000000446 fuel Substances 0.000 description 1
- 238000005272 metallurgy Methods 0.000 description 1
- 238000005065 mining Methods 0.000 description 1
- 238000010606 normalization Methods 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 238000000746 purification Methods 0.000 description 1
- 230000011218 segmentation Effects 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
- 230000001960 triggered effect Effects 0.000 description 1
- 238000013024 troubleshooting Methods 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
- 239000002699 waste material Substances 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/23—Clustering techniques
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/21—Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
- G06F18/213—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods
- G06F18/2135—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods based on approximation criteria, e.g. principal component analysis
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/21—Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
- G06F18/214—Generating training patterns; Bootstrap methods, e.g. bagging or boosting
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/24—Classification techniques
- G06F18/243—Classification techniques relating to the number of classes
- G06F18/24323—Tree-organised classifiers
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
- G06N3/00—Computing arrangements based on biological models
- G06N3/02—Neural networks
- G06N3/08—Learning methods
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02P—CLIMATE CHANGE MITIGATION TECHNOLOGIES IN THE PRODUCTION OR PROCESSING OF GOODS
- Y02P90/00—Enabling technologies with a potential contribution to greenhouse gas [GHG] emissions mitigation
- Y02P90/80—Management or planning
- Y02P90/84—Greenhouse gas [GHG] management systems
Landscapes
- Engineering & Computer Science (AREA)
- Data Mining & Analysis (AREA)
- Theoretical Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Engineering & Computer Science (AREA)
- Life Sciences & Earth Sciences (AREA)
- Artificial Intelligence (AREA)
- Computer Vision & Pattern Recognition (AREA)
- General Physics & Mathematics (AREA)
- Evolutionary Computation (AREA)
- Bioinformatics & Computational Biology (AREA)
- Evolutionary Biology (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Health & Medical Sciences (AREA)
- Biomedical Technology (AREA)
- Biophysics (AREA)
- Computational Linguistics (AREA)
- General Health & Medical Sciences (AREA)
- Molecular Biology (AREA)
- Computing Systems (AREA)
- Mathematical Physics (AREA)
- Software Systems (AREA)
- Testing And Monitoring For Control Systems (AREA)
Abstract
本发明实施例提供一种非稳定工况下能耗模式识别方法及一种非稳定工况下能耗模式识别系统，属于生产调控技术领域。所述方法包括：获取流程工业生产系统正常负荷工况下的时间序列数据；对时间序列数据进行预处理，并对预处理后的时间序列数据进行特征提取，获得工况数据集；通过预设切割算法对工况数据集进行切割，获得非平稳序列多变量时间序列样本集；根据非平稳序列多变量时间序列样本集构建多变量时间序列样本的相似性度量矩阵，并进行基于相似性度量矩阵的层次聚类；根据聚类结果进行非平稳过程能耗监测模型训练；采用训练获得的非平稳过程能耗监测模型进行流程工业生产系统非平稳过程监测。本发明方案对保证流程工业过程稳定快速过渡具有重要意义。
Description
技术领域
本发明涉及生产调控技术领域，具体地涉及方法一种非稳定工况下能耗模式识别方法及一种非稳定工况下能耗模式识别系统。
背景技术
流程行业包含水泥、钢铁、化工、造纸、冶金、电力等行业，不仅为机械、航空航天、军工、建筑提供原材料，也为国家发展提供电力等重要能源，在国民经济中占主导地位。其中钢铁、化工和水泥行业等行业为典型的高能耗和高二氧化碳排放的工业，在国家节能环保政策日益趋严的形势下，流程工业大型复杂设备的生产节能降耗是当前关注的热点问题。对于连续生产过程而言，设备耦合、系统规模庞大、生产连续、调控频繁的特性及高温、高压、高速旋转的工作环境会导致设备运行状态呈现变负荷、多工况、高能耗、和大范围非平稳特性，当出现外界扰动(原料/燃料成分变动)时，系统运行状态极不稳定，往往会带来设备振动增大、热量散失、空机运行甚至停机风险，这些都在无形中加大了系统的能耗浪费。然而由于操作人员的调控决策和实施水平的参差不齐，一方面，系统的高效低耗得不到保证，且依赖人工经验的调控过程且多聚焦于状态平稳、高产量、高品质的目标，多将低消耗高能效剥离出来且未给予足够重视，另一方面不当的调控策略会使系统从非平稳工况向平稳工况过渡时产生重大异常甚至导致故障停机。针对现有人工经验调试方法存在的诸多弊端，需要创造一种非稳定工况下能耗模式识别方法。
发明内容
本发明实施方式的目的是提供一种非稳定工况下能耗模式识别及一种非稳定工况下能耗模式识别系统，以至少解决现有操作人员的调控决策和实施水平的参差不齐，系统的高效低耗得不到保证甚至产生重大故障的问题。
为了实现上述目的，本发明第一方面提供一种非稳定工况下能耗模式识别方法，应用于流程工业生产系统在非稳定运行工况下的能耗状态模拟，所述方法包括：获取流程工业生产系统正常负荷工况下的时间序列数据；对所述时间序列数据进行预处理，并对预处理后的时间序列数据进行特征提取，获得工况数据集；通过预设切割算法对所述工况数据集进行切割，获得非平稳序列多变量时间序列样本集；根据所述非平稳序列多变量时间序列样本集构建多变量时间序列样本的相似性度量矩阵，并进行基于相似性度量矩阵的层次聚类；根据聚类结果进行非平稳过程能耗监测模型训练；采用训练获得的非平稳过程能耗监测模型进行流程工业生产系统非平稳过程监测。
可选的，所述对所述时间序列数据进行预处理，至少包括：对所述时间序列数据进行时序数据重采样、状态中间变量构建、数据小波降噪和时间序列数据标准化。
可选的，通过随机森林特征提取模型对预处理后的时间序列数据进行特征提取。
可选的，所述通过预设切割算法对所述工况数据集进行切割中，所述预设切割算法包括：通过预设关键时间点标识获取各预设关键时段；根据时间序列进行各预设关键时段的时间序列数据拼接，获得非平稳序列；其中，所述预设关键时间点包括：非平稳序列起点、前过渡期起点、调控组合动作起点、调控组合动作终点、后过渡期终点和非平稳序列终点；所述预设关键时段为分别对应每两个相邻预设关键时间点之间的时段。
可选的，所述根据所述非平稳序列多变量时间序列样本集构建多变量时间序列样本的相似性度量矩阵，并进行基于相似性度量矩阵的层次聚类，包括：分别计算所述非平稳序列多变量时间序列样本集中每两个样本之间的 DTW距离、时间序列统计特性和多变量时滞互相关系数；根据每两个样本之间的DTW距离、时间序列统计特性距离和多变量时滞互相关系数距离构建多变量时间序列样本的相似性度量矩阵；根据预设凝聚层次聚类算法将所述多变量时间序列样本的相似性度量矩阵聚类为多个簇；其中，所述多变量时间序列样本的相似性度量矩阵构建关系式为：
其中，SM(Zr-1,Zr)表示第r-1个非平稳序列多变量时间序列样本与第r 个非平稳序列多变量时间序列样本之间的相似距离；每两个非平稳序列多变量时间序列样本之间的相似距离计算公式为：
SM(Zx,Zy)＝SMDTW(Zx,Zy)+SMSC(Zx,Zy)+SMTLCC(Zx,Zy)
其中，SMDTW(Zx,Zy)表示样本x与样本y之间的DTW距离；
SMSC(Zx,Zy)表示样本x与样本y之间的时间序列统计特性距离；
SMTLCC(Zx,Zy)表示样本x与样本y之间的多变量时滞互相关系数距离。
可选的，所述非平稳序列多变量时间序列样本集中每两个样本之间的 DTW距离计算公式为：
SMDTW(Zx,Zy)＝w1*D(zx1,zy1)+w2*D(zx2,zy2)+…+wp*D(zxp,zyp)
其中，w1,w2,…,wp分别为基于第1,2，…，p个变量相对能耗变量的重要性权值；所述非平稳序列多变量时间序列样本集中每两个样本之间的时间序列统计特性距离计算公式为：
其中，fxjk表示样本x中第j个变量的第k个统计特性值；所述统计特征向量中的统计特性包括方差、均值、极差、偏度、峭度、熵和KL散度；所述非平稳序列多变量时间序列样本集中每两个样本之间的多变量时滞互相关系数距离计算公式为：
其中，TLCCxjq表示样本x中第j个变量和第q个变量之间的时滞互相关系数。
可选的，所述非平稳过程能耗模型训练的训练模型为核主成分分析模型 (KernalPrincipal Component Analysis,KPCA)；所述根据聚类结果进行非平稳过程能耗监测模型训练，包括：分别获取多个簇中每一个簇的非平稳序列多变量时间序列样本数据集，并根据所述非平稳序列多变量时间序列样本数据集构建对应簇的能耗模式的KPCA模型；根据获得的各KPCA模型 Hotelling-T2和SPE统计量的控制限形成非平稳过程能耗模式控制限向量集演化图；根据各KPCA模型结合所述非平稳过程能耗模式控制限向量集演化图获得非平稳过程能耗监测模型。
可选的，所述通过训练获得的非平稳过程能耗监测模型进行流程工业生产系统非平稳过程监测，包括：获取流程工业生产系统的运行参数；根据所述运行参数判断所述流程工业生产系统的当前能耗模态；根据所述流程工业生产系统的当前能耗模态选择对应的非平稳过程能耗监测模型；将所述流程工业生产系统的当前运行参数的演化路径与对应的非平稳过程能耗监测模型的演化路径进行对比，在对比结果满足报警条件时输出报警指令；执行所述报警指令，生成对应的报警信息。
本发明第二方面提供一种非稳定工况下能耗模式识别系统，应用于流程工业生产系统在非稳定运行工况下的能耗状态模拟，所述系统包括：采集单元，用于获取流程工业生产系统正常负荷工况下的时间序列数据；处理单元，用于：对所述时间序列数据进行预处理，并对预处理后的时间序列数据进行特征提取，获得工况数据集；通过预设切割算法对所述工况数据集进行切割，获得非平稳序列多变量时间序列样本集；根据所述非平稳序列多变量时间序列样本集构建多变量时间序列样本的相似性度量矩阵，并进行基于相似性度量矩阵的层次聚类；训练单元，用于根据聚类结果进行非平稳过程能耗监测模型训练；监测单元，用于采用训练获得的非平稳过程能耗监测模型进行流程工业生产系统非平稳过程监测。
另一方面，本发明提供一种计算机可读储存介质，该计算机可读存储介质上储存有指令，其在计算机上运行时使得计算机执行上述的非稳定工况下能耗模式识别方法。
通过上述技术方案，通过对既往历史运行工况数据的挖掘，考虑了能耗状态的关键变量，将具有相同潜在变化特性的能耗模式利用同一个模型进行表征，获得不同能耗模式的控制限向量演化图，在在线监测中，通过对比当前运行状态的统计量与不同能耗模式辨识模型的控制限，来判断系统从非平稳工况向平稳工况的过渡是否正常，该方法实时获取过程数据保证了判断的可靠性，对保证流程工业过程稳定快速过渡具有重要意义。
本发明实施方式的其它特征和优点将在随后的具体实施方式部分予以详细说明。
附图说明
附图是用来提供对本发明实施方式的进一步理解，并且构成说明书的一部分，与下面的具体实施方式一起用于解释本发明实施方式，但并不构成对本发明实施方式的限制。在附图中：
图1是本发明一种实施方式提供的非稳定工况下能耗模式识别方法的步骤流程图；
图2是本发明一种实施方式提供的非平稳序列划分示意图；
图3是本发明一种实施方式提供的多变量时间序列样本的相似性度量矩阵构建流程图；
图4是本发明一种实施方式提供的非平稳过程能耗监测模型训练步骤流程图；
图5是本发明一种实施方式提供的流程工业生产系统非平稳过程监测步骤流程图；
图6是本发明一种实施方式提供的非稳定工况下能耗模式识别系统的系统结构图。
附图标记说明
10-采集单元；20-处理单元；30-训练单元；40-监测单元。
具体实施方式
以下结合附图对本发明的具体实施方式进行详细说明。应当理解的是，此处所描述的具体实施方式仅用于说明和解释本发明，并不用于限制本发明。
图6是本发明一种实施方式提供的非稳定工况下能耗模式识别系统的系统结构图。如图6所示，本发明实施方式提供一种非稳定工况下能耗模式识别系统，所述系统包括：采集单元10，用于获取流程工业生产系统正常负荷工况下的时间序列数据；处理单元20，用于：对所述时间序列数据进行预处理，并对预处理后的时间序列数据进行特征提取，获得工况数据集；通过预设切割算法对所述工况数据集进行切割，获得非平稳序列多变量时间序列样本集；根据所述非平稳序列多变量时间序列样本集构建多变量时间序列样本的相似性度量矩阵，并进行基于相似性度量矩阵的层次聚类；训练单元30，用于根据聚类结果进行非平稳过程能耗监测模型训练；监测单元40，用于采用训练获得的非平稳过程能耗监测模型进行流程工业生产系统非平稳过程监测。
图1是本发明一种实施方式提供的一种非稳定工况下能耗模式识别方法的方法流程图。如图1所示，本发明实施方式提供一种非稳定工况下能耗模式识别方法，所述方法包括：
步骤S10：获取流程工业生产系统正常负荷工况下的时间序列数据。
具体的，进行流程工业生产系统非稳定工况具体能耗监测模型分析，就必须获取对应流程工业生产系统的运行参数，然后根据运行参数进行对应系统运行工况模拟，根据模拟情况获取整个运行过程中的能耗状态。因为非稳定工况与理论运行工况存在区别，且内外部不确定性干扰因素众多，所以直接获取非稳定状态下的能耗曲线，对于后续的参考意义也不会太大，因为两次随机情况的能耗曲线一致的概率非常低。但只要系统是在正常工作，那即使在非稳定工况下，最终的能耗曲线依然会回落到正常曲线上来，所以获取系统正常负荷工况下的运行数据，判断非稳定工况下的曲线趋势，若最终趋势回落到稳定状态，那便意味着系统由非稳定工况能够正常过渡到稳定工况，并未发生故障或重大异常。所以需要准确获取系统在正常负荷工况下的运行参数，优选的，在系统主控模块进行采集单元10构建，通过系统本身的运行参数存储库进行系统正常负荷下的时间序列数据提取。这些数据包括系统在正常负荷工况下的能耗情况和对应的时间身份信息，通过整理时间序列数据，可以获得系统在正常负荷工况下的能耗曲线。
步骤S20：对所述时间序列数据进行预处理，并对预处理后的时间序列数据进行特征提取，获得工况数据集。
具体的，时间序列数据包括很多设备身份信息和协议信息等无用信息，这些信息对于本发明提出的能耗模式识别无参考价值，且这些信息体量很大，若混合在原始数据中进行分析，会极大增加模型训练负担，造成训练干扰的同时还为隐形训练效率。所以在进行数据分析前，需要对原始的时间序列数据进行预处理，进行信息提纯。时序数据重采样即将时间序列由高频向低频转换，然后基于原始过程变量构建能耗状态表征变量M。对时域信号，选择小波基和分解层数，对信号进行小波分解，对分解得到的各层系数进行阈值处理，最后重构去噪后的小波系数获取去噪后的信号。由于变量间量纲不同，因此对降噪后数据进行标准化处理。这些预处理方法包括时序数据重采样、状态中间变量构建、数据小波降噪和时间序列数据标准化。通过预处理后的数据，无论是数据内容还是数据标准，均适用于后续的分析条件，便于后续分析使用。基于预处理后的数据，选择能耗状态表征变量M为结果变量，利用随机森林特征选择算法，计算各变量的重要性得分，根据得分排序分别赋予不同权重(w1,w2,…,wp)，获取能耗状态的工况特征集 (P1,P2,…,Pp)。此时获取的工况数据集为：
Zj(t)(j＝1,2,…,p；t＝1,2,…,m)
其中，p为变量数，m为观测数，t为索引每个时刻的观测值。一个多变量时间序列对象Zi包括p≥2的观测数，可以由以下m×p矩阵表示：
其中，zij(t)元素对应对象Zi(i＝1,2,…,n；n＝m×p)中变量j(j＝1,2,…,p) 在时刻t(t＝1,2,…,m)的一个观测值。列j包含了与变量j相关的时间序列。
步骤S30：通过预设切割算法对所述工况数据集进行切割，获得非平稳序列多变量时间序列样本集。
具体的，如图2，通过预设关键时间点标识获取各预设关键时段；根据时间序列进行各预设关键时段的时间序列数据拼接，获得非平稳序列；其中，所述预设关键时间点为6个，包括：非平稳序列起点、前过渡期起点、调控组合动作起点、调控组合动作终点、后过渡期终点和非平稳序列终点；所述预设关键时段为5个，分别对应每两个相邻预设关键时间点之间的时段。5 个关键时段即前稳定期、前过渡期、调控期、后过渡期、后稳定期，由5个关键时段共同组合成一个完整的非平稳序列。根据调控动作的6个关键时间点，对工况数据集Zj(t)(j＝1,2,…,p；t＝1,2,…,m)基于关键时间节点进行序列分割，获得S个非平稳序列多变量时间序列样本集(Z1,Z2,…Zs),非平稳序列多变量时间序列样本Zs，其表达式为：
其中，ms表示观测数，p表示变量数。
步骤S40：根据所述非平稳序列多变量时间序列样本集构建多变量时间序列样本的相似性度量矩阵，并进行基于相似性度量矩阵的层次聚类。
具体的，本发明通过对既往历史运行工况数据的挖掘，考虑了能耗状态的关键变量，将具有相同潜在变化特性的能耗模式利用同一个模型进行表征，获得不同能耗模式的控制限向量演化图，在在线监测中，通过对比当前运行状态的统计量与不同能耗模式辨识模型的控制限，来判断非平稳工况向平稳工况过渡是否正常。所以需要进行工况类别化，即筛选出具有相同潜在变化特性的能耗模式。具体的，如图3，包括以下步骤：
步骤S401：计算非平稳序列多变量时间序列样本中每两个样本之间的DTW距离。
具体的，对于两个不同的非平稳序列多变量时间序列样本Zx和Zy(分别由两个mx×p和my×p矩阵表示的多变量时间序列，变量数均为p，观测数分别为mx和my)，根据预设的DTW算法计算两个样本之间基于特征加权的多变量DTW距离，其计算公式为：
SMDTW(Zx,Zy)＝w1*D(zx1,zy1)+w2*D(zx2,zy2)+…+wp*D(zxp,zyp)
其中，w1,w2,…,wp分别为基于第1,2，…，p个变量相对能耗变量的重要性权值。
步骤S402：分析非平稳序列多变量时间序列样本中每两个样本之间的时间序列统计特性距离。
具体的，对于两个单变量时间序列Zxj和Zyj，获取他们的统计特性包括方差、均值、极差、偏度、峭度、熵、KL散度等数据分布特征，分别构成统计特性向量fxj(fxj1,…,fxj7)和统计特性向量fyj(fyj1,…,fyj7)。根据预设的分析公式进行Zxj和Zyj之前的时间序列统计特性相似距离分析，预设分析公式为：
则对于两个不同的非平稳序列多变量时间序列样本Zx和Zy(分别由两个m×p矩阵表示的多变量时间序列，变量数均为p)，根据预设的分析公式进行Zx和Zy之前的时间序列统计特性相似距离分析，预设分析公式为：
步骤S403：计算非平稳序列多变量时间序列样本中每两个样本之间的多变量时滞互相关系数距离。
具体的，确定时间滞后系数H，其中， -min(m1,m2,…,ms)＜H≤min(m1,m2,…,ms)；对于非平稳序列多变量时间序列样本Zx，其关键特征为(P1,P2,…Pp)，使用时滞互相关算法计算该样本的两个单变量时间序列Zxj(j＝1,2,…p)和Zxq(q＝1,2,…p)间的时滞互相关系数TLCCxjq。如果TLCCxjq＞0，则Zxj在引导Zxq，即前者先于后者变化；如果 TLCCxjq＜0，则Zxq在引导Zxj，即后者先于前者变化。通过计算，获得样本Zx的时滞相关系数集为：
(TLCCx12,TLCCx13,…,TLCCxjq,…,TLCCx(p-1)q)(j＝1,2,…,p；q＝1,2,…,p)
根据预设公式计算Zx和Zy之前的多变量时滞互相关系数相似距离，其计算公式为：
步骤S404：根据每两个样本之间的DTW距离、时间序列统计特性距离和多变量时滞互相关系数距离构建多变量时间序列样本的相似性度量矩阵。
具体的，根据如下公式计算Zx和Zy之前的相似距离：
SM(Zx,Zy)＝SMDTW(Zx,Zy)+SMSC(Zx,Zy)+SMTLCC(Zx,Zy)
构建S个非平稳序列多变量时间序列样本的相似性度量矩阵，该矩阵为对角线对0的上三角矩阵：
步骤S405：进行基于相似性度量矩阵的层次聚类。
在获取所述相似性度量矩阵后，使用凝聚层次聚类算法(ward linkage) 对非平稳序列多变量时间序列样本集聚类为G个簇，即分别选取两簇中所有节点的平均距离作为簇间距离，将具有最短距离的类别相合并，依次迭代，直至最终聚合成一类，聚类过程以可视化的聚类树呈现。对于新的样本点，根据如下公式计算所属簇：
其中，u是s和t组成的新的聚类，v是未使用的聚类，|v|、|s|和|t|分别是v、s和t中元素的个数，T＝|v|+|s|+|t|。根据以下轮廓系统公式确定G的最佳值：
其中，a表示平均簇内距离，即样本与同簇中其他样本点的距离均值。b 表示平均最近簇距离，即样本与其他邻近簇的所有样本点的距离均值。
步骤S50：根据聚类结果进行非平稳过程能耗监测模型训练。
具体的，对聚类结果中各簇非平稳序列多变量时间序列样本做加窗 KPCA建模，形成能耗模式控制限向量集演化图。由于非平稳过程的变量数据呈现出非线性、非高斯特性和动态特性，因此使用核主成分分析(Kernel Principal Component Analysis,KPCA)模型估计不同能耗模式的数据模式。分别获取G个簇的非平稳序列多变量时间序列样本数据集，分别构建G个能耗模式的KPCA模型，即对具有相同过程特性的能耗模式建立统一的辨识模型，对不同模式建立不同模型。如图4，单个能耗模式的离线KPCA建模过程具体包括：
步骤S501：基于前稳定期对非平稳序列多变量时间序列样本进行KPCA 建模窗口分割。
具体的，设定前稳定期(即P1至P2时间点区间)时间长度为h，样本Zs共分割为L个长度为h的窗口，其中：
每个窗口的数据表示为
(由h×p矩阵表示，变量数均为p，观测数分别为，l＝1,2,…,L)。针对各个窗口数据分别构建KPCA模型，鉴于
为前稳定期数据，因此将第一个窗口的KPCA模型记为参考模型KPCAref。
步骤S502：针对各个窗口的数据样本构建KPCA。
具体的，KPCA基于核函数原理，利用非线性函数φ将原始输入空间的非线性数据映射到高维空间F中，变成线性可分的数据，然后在高维空间中进行PCA变换。设标准化原始空间数据xn∈RJ(n＝1,2,…,N)，则xn在特征空间F对应的映射数据记为φ(xn)，即：
xn,xn,…,xn∈RJ→φ(x1),φ(x2),…φ(xn)∈RH
φ(X)的协方差矩阵表示为：
存在系数an(c＝1,2,…,N)使∑F特征向量ν能够线性表示，结合对应的特征值λ可得到：
其中，定义矩阵:
K∈RN×N
kij＝<φ(xi),φ(xj)>(i,j＝1,2,…,N)
其中，kij是核函数。计算映射数据φ(X)的非线性主成分得分，其计算公式为：
KPCA特征空间中核主成分模型内部波动与PCA类似，使用主成分子空间的Hotelling-T2统计量来衡量，KPCA中某一时刻采样点到模型空间距离与PCA类似，使用残差子空间的SPE统计量来衡量。然后计算所有建模数据的T2和SPE统计量，其中，Hotelling-T2统计量定义为：
T2＝[t1,t2,…,tA]Λ-1[t1,t2,…,tA]T
其中，Λ是得分向量相对应的协方差矩阵，T2统计量的控制限根据以下 F分布得到：
其中，N为模型中样本个数，A为保留的核主成分个数，a为显著性检验水平。则SPE统计量的定义为：
最后确定T2和SPE统计量的控制限，记为向量(T2(1),SPE(1))，也记为 (T2(ref),SPE(ref))。
步骤S503：绘制单个能耗模式控制限向量集演化图。
具体的，获取上述L个KPCA模型的T2和SPE统计量的控制限向量集， L个KPCA模型的T2和SPE统计量的控制限向量集记为：
((T2(1),SPE(1)),(T2(2),SPE(2)),…,(T2(l),SPE(l)))
根据预设公式计算相邻窗口的控制向量余弦距离，余弦值的正负和大小分别代表了相邻窗口控制限向量的变化方向和相似程度，预设的计算公式为：
依据cosθ绘制控制限向量集演化图，即表征了单个能耗模式的过程相关性特性。
步骤S504：绘制G个能耗模式控制限向量集演化图。
具体的，重复执行步骤S503，获得G个能耗模式控制限向量集演化图。 G个能耗模式的控制限向量集演化图分别描述了G种不同的非平稳过程相关性特性。
步骤S60：采用训练获得的非平稳过程能耗监测模型进行流程工业生产系统非平稳过程监测。
具体的，获得非平稳过程能耗监测模型后，后续便可直接使用该模型进行对应流程系统监测，即在后续系统运行过程中，实时采集系统的运行参数，根据实时的运行参数模拟出当前的系统运行工况，然后对应获得非平稳过程能耗监测模型，对比当前运行参数的演化路径和对应非平稳过程能耗模型演化路径，并在两个路径差距大于预设阈值时，输出报警指令；执行所述报警指令，生成对应的报警信息。在线监测时，对于新的过程数据，进行加窗判断当前状态属于哪个模态，从而使用对应的能耗模式监测模型监测并指导当前状态的过渡过程的异常与否。进行在线监测时，具体的，如图5，主要包括以下步骤：
步骤S601：获取在线监测数据。
具体的，设当前采样时刻t，取历史非平稳序列多变量时间序列样本中最短的前稳定期长度lmin，从第t-lmin+1个采样到第t个采样时刻的连续lmin个在线样本数据进行在线监测。
步骤S602：进行数据分析。
具体的，根据当前状态窗口数据
(由lmin×p矩阵表示的多变量时间序列数据)的统计量(T2(1),SPE(1))，与G个能耗模式的统计量(T2(ref),SPE(ref)) 进行对比，选择最为接近的能耗模式控制限向量集演化图作为参考演化图进行监测。
步骤S603：将实际演化图与参考演化图对比，判断是否符合正常的非稳定工况向稳定工况过渡的过程。
具体的，循环选取lmin长度的窗口数据计算统计量T2和SPE，将其与参考演化图对比，若控制限向量集与参考演化图演化轨迹一致，则表示过渡过程正常，若不一致，则表示出现了故障模态或新模态，则生成报警指令，并生成对应的报警信息，提醒相关人员进行故障排查。优选的，在另一种可能的实施方式中，为避免频繁报警，提高系统的稳定性能，在生成报警信息前，先生成人工审核的辅助信息，供人工进行工况审核，并在人员确定故障后再生成报警指令，若相关人员通过经验判断系统正常，则不触发报警，保持系统正常运行。
本发明实施方式还提供一种计算机可读储存介质，该计算机可读存储介质上储存有指令，其在计算机上运行时使得计算机执行上述的非稳定工况下能耗模式识别方法。
本领域技术人员可以理解实现上述实施方式的方法中的全部或部分步骤是可以通过程序来指令相关的硬件来完成，该程序存储在一个存储介质中，包括若干指令用以使得单片机、芯片或处理器(processor)执行本发明各个实施方式所述方法的全部或部分步骤。而前述的存储介质包括：U盘、移动硬盘、只读存储器(ROM，Read-Only Memory)、随机存取存储器(RAM， Random Access Memory)、磁碟或者光盘等各种可以存储程序代码的介质。
以上结合附图详细描述了本发明的可选实施方式，但是，本发明实施方式并不限于上述实施方式中的具体细节，在本发明实施方式的技术构思范围内，可以对本发明实施方式的技术方案进行多种简单变型，这些简单变型均属于本发明实施方式的保护范围。另外需要说明的是，在上述具体实施方式中所描述的各个具体技术特征，在不矛盾的情况下，可以通过任何合适的方式进行组合。为了避免不必要的重复，本发明实施方式对各种可能的组合方式不再另行说明。
此外，本发明的各种不同的实施方式之间也可以进行任意组合，只要其不违背本发明实施方式的思想，其同样应当视为本发明实施方式所公开的内容。
Claims (8)
1.一种非稳定工况下能耗模式识别方法，应用于流程工业生产系统在非稳定运行工况下的能耗状态模拟，其特征在于，所述方法包括：
获取流程工业生产系统正常负荷工况下的时间序列数据；
对所述时间序列数据进行预处理，并对预处理后的时间序列数据进行特征提取，获得工况数据集；
通过预设切割算法对所述工况数据集进行切割，获得非平稳序列多变量时间序列样本集；其中，所述预设切割算法包括：
通过预设关键时间点标识获取各预设关键时段；
根据时间序列进行各预设关键时段的时间序列数据拼接，获得非平稳序列；其中，所述预设关键时间点包括：
非平稳序列起点、前过渡期起点、调控组合动作起点、调控组合动作终点、后过渡期终点和非平稳序列终点；
所述预设关键时段为分别对应每两个相邻预设关键时间点之间的时段；
根据所述非平稳序列多变量时间序列样本集构建多变量时间序列样本的相似性度量矩阵，并进行基于相似性度量矩阵的层次聚类；包括：
分别计算所述非平稳序列多变量时间序列样本集中每两个样本之间的DTW距离、时间序列统计特性距离和多变量时滞互相关系数距离；
根据每两个样本之间的DTW距离、时间序列统计特性距离和多变量时滞互相关系数距离构建多变量时间序列样本的相似性度量矩阵；
根据预设凝聚层次聚类算法将所述多变量时间序列样本的相似性度量矩阵聚类为多个簇；其中，
所述多变量时间序列样本的相似性度量矩阵构建关系式为：
其中，SM(Zr-1,Zr)表示第r-1个非平稳序列多变量时间序列样本与第r个非平稳序列多变量时间序列样本之间的相似距离；
每两个非平稳序列多变量时间序列样本之间的相似距离计算公式为：
SM(Zx,Zy)＝SMDTW(Zx,Zy)+SMSC(Zx,Zy)+SMTLCC(Zx,Zy)
其中，SMDTW(Zx,Zy)表示样本x与样本y之间的DTW距离；
SMSC(Zx,Zy)表示样本x与样本y之间的时间序列统计特性距离；
SMTLCC(Zx,Zy)表示样本x与样本y之间的多变量时滞互相关系数距离；
根据聚类结果进行非平稳过程能耗监测模型训练；
采用训练获得的非平稳过程能耗监测模型进行流程工业生产系统非平稳过程监测。
2.根据权利要求1所述的方法，其特征在于，所述对所述时间序列数据进行预处理，至少包括：
对所述时间序列数据进行时序数据重采样、状态中间变量构建、数据小波降噪和时间序列数据标准化。
3.根据权利要求1所述的方法，其特征在于，通过随机森林特征提取模型对预处理后的时间序列数据进行特征提取。
4.根据权利要求1所述的方法，其特征在于，所述非平稳序列多变量时间序列样本集中每两个样本之间的DTW距离计算公式为：
SMDTW(Zx,Zy)＝w1*D(zx1,zy1)+w2*D(zx2,zy2)+…+wp*D(zxp,zyp)
其中，w1,w2,…,wp分别为基于第1,2，…，p个变量相对能耗变量的重要性权值；
所述非平稳序列多变量时间序列样本集中每两个样本之间的时间序列统计特性距离公式为：
其中，fxjk表示样本x的第j个变量的第k个统计特性值；统计特征向量中的统计特性包括方差、均值、极差、偏度、峭度、熵和KL散度；
所述非平稳序列多变量时间序列样本集中每两个样本之间的多变量时滞互相关系数距离计算公式为：
其中，TLCCxjq表示样本x的第j个变量和第q个变量之间的时滞互相关系数。
5.根据权利要求4所述的方法，其特征在于，所述非平稳过程能耗监测模型训练的训练模型为核主成分分析模型；
所述根据聚类结果进行非平稳过程能耗监测模型训练，包括：
分别获取多个簇中每一个簇的非平稳序列多变量时间序列样本数据集，并根据所述非平稳序列多变量时间序列样本数据集构建对应簇的能耗模式的KPCA模型；
根据获得的所有KPCA模型形成非平稳过程能耗模式控制限向量集演化图；
根据各KPCA模型和模型对应的非平稳过程能耗模式控制限向量集演化图获得非平稳过程能耗监测模型。
6.根据权利要求5所述的方法，其特征在于，所述通过训练获得的非平稳过程能耗监测模型进行流程工业生产系统非平稳过程监测，包括：
获取流程工业生产系统的运行参数；
根据所述运行参数判断所述流程工业生产系统的当前能耗模态；
根据所述流程工业生产系统的当前能耗模态选择对应的非平稳过程能耗监测模型；
将所述流程工业生产系统的当前运行参数的演化路径与对应的非平稳过程能耗监测模型的演化路径进行对比，在对比结果满足报警条件时输出报警指令；
执行所述报警指令，生成对应的报警信息。
7.一种非稳定工况下能耗模式识别系统，应用于流程工业生产系统在非稳定运行工况下的能耗状态模拟，其特征在于，所述系统包括：
采集单元，用于获取流程工业生产系统正常负荷工况下的时间序列数据；
处理单元，用于：
对所述时间序列数据进行预处理，并对预处理后的时间序列数据进行特征提取，获得工况数据集；
通过预设切割算法对所述工况数据集进行切割，获得非平稳序列多变量时间序列样本集；其中，所述预设切割算法包括：
通过预设关键时间点标识获取各预设关键时段；
根据时间序列进行各预设关键时段的时间序列数据拼接，获得非平稳序列；其中，所述预设关键时间点包括：
非平稳序列起点、前过渡期起点、调控组合动作起点、调控组合动作终点、后过渡期终点和非平稳序列终点；
所述预设关键时段为分别对应每两个相邻预设关键时间点之间的时段；
根据所述非平稳序列多变量时间序列样本集构建多变量时间序列样本的相似性度量矩阵，并进行基于相似性度量矩阵的层次聚类；包括：
分别计算所述非平稳序列多变量时间序列样本集中每两个样本之间的DTW距离、时间序列统计特性距离和多变量时滞互相关系数距离；
根据每两个样本之间的DTW距离、时间序列统计特性距离和多变量时滞互相关系数距离构建多变量时间序列样本的相似性度量矩阵；
根据预设凝聚层次聚类算法将所述多变量时间序列样本的相似性度量矩阵聚类为多个簇；其中，
所述多变量时间序列样本的相似性度量矩阵构建关系式为：
其中，SM(Zr-1,Zr)表示第r-1个非平稳序列多变量时间序列样本与第r个非平稳序列多变量时间序列样本之间的相似距离；
每两个非平稳序列多变量时间序列样本之间的相似距离计算公式为：
SM(Zx,Zy)＝SMDTW(Zx,Zy)+SMSC(Zx,Zy)+SMTLCC(Zx,Zy)
其中，SMDTW(Zx,Zy)表示样本x与样本y之间的DTW距离；
SMSC(Zx,Zy)表示样本x与样本y之间的时间序列统计特性距离；
SMTLCC(Zx,Zy)表示样本x与样本y之间的多变量时滞互相关系数距离；
训练单元，用于根据聚类结果进行非平稳过程能耗监测模型训练；
监测单元，用于采用训练获得的非平稳过程能耗监测模型进行流程工业生产系统非平稳过程监测。
8.一种计算机可读储存介质，该计算机可读存储介质上储存有指令，其在计算机上运行时使得计算机执行权利要求1至6中任一项权利要求所述的非稳定工况下能耗模式识别方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110613823.1A CN113298170B (zh) | 2021-06-02 | 2021-06-02 | 非稳定工况下能耗模式识别方法及系统 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110613823.1A CN113298170B (zh) | 2021-06-02 | 2021-06-02 | 非稳定工况下能耗模式识别方法及系统 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113298170A CN113298170A (zh) | 2021-08-24 |
| CN113298170B true CN113298170B (zh) | 2022-09-09 |
Family
ID=77326892
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110613823.1A Active CN113298170B (zh) | 2021-06-02 | 2021-06-02 | 非稳定工况下能耗模式识别方法及系统 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113298170B (zh) |
Families Citing this family (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114637263B (zh) * | 2022-03-15 | 2024-01-12 | 中国石油大学(北京) | 一种异常工况实时监测方法、装置、设备及存储介质 |
| CN114646477B (zh) * | 2022-04-01 | 2023-05-02 | 丰发伟业(北京)科技有限公司 | 一种汽车混合动力系统的测试方法及系统 |
| CN117421618B (zh) * | 2023-11-24 | 2024-04-05 | 上海东方低碳科技产业股份有限公司 | 一种建筑能耗监测方法及系统 |
| CN118366108A (zh) * | 2024-05-06 | 2024-07-19 | 嘉兴南湖区路空协同立体交通产业研究院 | 一种路空协同交通场景智能识别方法 |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108304567A (zh) * | 2018-02-09 | 2018-07-20 | 杭州安脉盛智能技术有限公司 | 高压变压器工况模式识别与数据分类方法及系统 |
| CN110879580A (zh) * | 2019-12-10 | 2020-03-13 | 浙江大学 | 一种面向大范围非平稳瞬变连续过程的分析和监测方法 |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US11264799B2 (en) * | 2018-11-07 | 2022-03-01 | State Grid of China, State Grid Jiangsu Electric Power Co Ltd | Submodular load clustering with robust principal component analysis |
| CN110909912B (zh) * | 2019-10-08 | 2022-08-09 | 西安交通大学 | 基于自适应误差反馈的园区电力系统净负荷组合预测方法 |
-
2021
- 2021-06-02 CN CN202110613823.1A patent/CN113298170B/zh active Active
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108304567A (zh) * | 2018-02-09 | 2018-07-20 | 杭州安脉盛智能技术有限公司 | 高压变压器工况模式识别与数据分类方法及系统 |
| CN110879580A (zh) * | 2019-12-10 | 2020-03-13 | 浙江大学 | 一种面向大范围非平稳瞬变连续过程的分析和监测方法 |
Non-Patent Citations (1)
| Title |
|---|
| 基于非参数核密度估计和改进谱多流形聚类的负荷曲线分类研究;高亚静等;《电网技术》;20171016(第05期);全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113298170A (zh) | 2021-08-24 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN113298170B (zh) | 非稳定工况下能耗模式识别方法及系统 | |
| CN111275288B (zh) | 基于XGBoost的多维数据异常检测方法与装置 | |
| CN110262450B (zh) | 面向汽轮机的多种故障特性协同分析的故障预测方法 | |
| CN111949012B (zh) | 一种基于双权重多邻域保持嵌入算法的间歇过程故障检测方法 | |
| AU2012284460A1 (en) | System of sequential kernel regression modeling for forecasting and prognostics | |
| EP2791745B1 (en) | A method of operating a process or machine | |
| CN117290800B (zh) | 一种基于超图注意力网络的时序异常检测方法及系统 | |
| CN106094749B (zh) | 基于改进核熵成分分析的非线性故障检测方法及应用 | |
| CN113095402A (zh) | 一种基于编码输入的生成对抗网络故障检测方法及系统 | |
| CN111340110A (zh) | 一种基于工业过程运行状态趋势分析的故障预警方法 | |
| Gao et al. | A novel fault detection model based on vector quantization sparse autoencoder for nonlinear complex systems | |
| CN112418577A (zh) | 一种工业产品生产过程的可视化监控方法及系统 | |
| CN111062848A (zh) | 一种消防工程异常状态监测的智能监控方法 | |
| Khan et al. | Advanced statistical and meta-heuristic based optimization fault diagnosis techniques in complex industrial processes: a comparative analysis | |
| CN112947649B (zh) | 一种基于互信息矩阵投影的多变量过程监控方法 | |
| CN108171271B (zh) | 一种设备劣化早期预警方法和系统 | |
| Li et al. | Data driven condition monitoring of wind power plants using cluster analysis | |
| CN110347138A (zh) | 自适应核主元分析的化工过程故障诊断方法 | |
| CN113985853B (zh) | 一种基于数据依赖核判别分析的工业过程可视化监控方法 | |
| Kang et al. | A data-driven method of health monitoring for spacecraft | |
| Yang et al. | A novel decentralized weighted ReliefF-PCA method for fault detection | |
| CN109523195B (zh) | 基于样本多模态分类的分层式非高斯算法的故障检测方法 | |
| Lu et al. | Visual monitoring of industrial operation states based on kernel fisher vector and self-organizing map networks | |
| Xuemin et al. | A fault detection method using multi-scale kernel principal component analysis | |
| Fingerhut et al. | Data-Driven Usage Profiling and Anomaly Detection in Support of Sustainable Machining Processes |
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 |