CN110222393A - 基于细粒度风电发电状态划分的风机叶片结冰异常监测方法 - Google Patents
基于细粒度风电发电状态划分的风机叶片结冰异常监测方法 Download PDFInfo
- Publication number
- CN110222393A CN110222393A CN201910448367.2A CN201910448367A CN110222393A CN 110222393 A CN110222393 A CN 110222393A CN 201910448367 A CN201910448367 A CN 201910448367A CN 110222393 A CN110222393 A CN 110222393A
- Authority
- CN
- China
- Prior art keywords
- state
- wind speed
- sheet
- monitoring
- wind
- 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.)
- Granted
Links
- 238000012544 monitoring process Methods 0.000 title claims abstract description 103
- 238000000034 method Methods 0.000 title claims abstract description 51
- 230000005856 abnormality Effects 0.000 title abstract description 4
- 230000005611 electricity Effects 0.000 title abstract 3
- 230000002159 abnormal effect Effects 0.000 claims abstract description 26
- 239000011159 matrix material Substances 0.000 claims description 51
- 238000009826 distribution Methods 0.000 claims description 44
- 230000003068 static effect Effects 0.000 claims description 27
- 238000010248 power generation Methods 0.000 claims description 20
- 238000004364 calculation method Methods 0.000 claims description 17
- 238000000354 decomposition reaction Methods 0.000 claims description 12
- 230000002087 whitening effect Effects 0.000 claims description 6
- 238000006243 chemical reaction Methods 0.000 claims description 3
- 238000012545 processing Methods 0.000 claims description 3
- 238000002360 preparation method Methods 0.000 claims 1
- 238000012423 maintenance Methods 0.000 abstract description 6
- 238000001514 detection method Methods 0.000 abstract description 4
- 230000000694 effects Effects 0.000 abstract description 3
- 238000000605 extraction Methods 0.000 abstract description 2
- 230000000875 corresponding effect Effects 0.000 description 12
- 230000001133 acceleration Effects 0.000 description 3
- 230000002596 correlated effect Effects 0.000 description 3
- 238000011161 development Methods 0.000 description 3
- 238000004458 analytical method Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 2
- 230000000717 retained effect Effects 0.000 description 2
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 2
- 238000012512 characterization method Methods 0.000 description 1
- 230000009194 climbing Effects 0.000 description 1
- 238000005336 cracking Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 230000008014 freezing Effects 0.000 description 1
- 238000007710 freezing Methods 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000000638 solvent extraction Methods 0.000 description 1
- 238000006467 substitution reaction 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
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Optimization (AREA)
- Mathematical Analysis (AREA)
- Computational Mathematics (AREA)
- General Engineering & Computer Science (AREA)
- Data Mining & Analysis (AREA)
- Software Systems (AREA)
- Databases & Information Systems (AREA)
- Algebra (AREA)
- Computing Systems (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- Wind Motors (AREA)
Abstract
本发明公开了一种基于细粒度风电发电状态划分的风机叶片结冰异常检测方法。在揭示风机变工况运行时各变量关系的分段性特点后,针对性地提出了一种基于慢特征提取的细粒度状态划分建模方法。针对风力发电的动态变化特性,提出了一种动静态协同的监测方法。利用风场SCADA系统采集到的数据,用以上方法对风机建立各子状态的监测模型,并离线了验证本文提出的方法检测风机出力异常的效果。充分利用了风机运行时数据的动态特性,有效地提升了检测效果,有助于风场维护人员对叶片结冰的情况进行及时的诊断以及处理,从而保障了风力发电机组正常稳定的运行,同时提升了人员财产的安全保障系数。
Description
技术领域
本发明属于风力发电过程监测领域,特别是涉及一种基于细粒度风电发电状态划分的风 机叶片结冰异常检测方法。
背景技术
据行业统计,2018年1-9月,全国新增风电并网容量1261万千瓦,到9月底累计 风电并网容量达到1.76亿千瓦;1-9月,全国风电发电量2676亿千瓦时,同比增长26%; 平均利用小时数1565小时,同比增加178小时;1-9月,全国弃风电量222亿千瓦时,同 比减少74亿千瓦时。
与此同时,经济性仍是制约风电发展的重要因素。与传统的化石能源电力相比,风电的 发电成本仍比较高,补贴需求和政策依赖性较强,行业发展受政策变动影响较大。而在风电 项目的开发过程中,风机能否在运转时期发挥最佳性能是衡量风场投资成败的关键因素之一。 风机极易出现一些利用普通手段难以监测的细微故障,这些故障的发生往往会导致风力发电 机无法达到额定的工作点,例如风机叶片结冰。叶片结冰会导致风机产生无法并网的电能, 造成了资源的浪费甚至对电网的稳定运行造成干扰。结冰程度较为严重时会导致冰块脱落, 这往往会威胁到风场人员财产安全。同时因为风机结构众多,测点数很多,导致数据量十分 庞大,不利于对此类故障的分析以及监测。
风力发电产业的工作环境复杂多变,并且由于其需要并入电网运行,因而对发电机运行 的稳定性要求较高,并且需要及时的人员维护。但是其发生故障的机率较高,若不及时检测 并排除故障就会导致风场无法正常运行,并消耗大量资金用于维修。例如:2016年2月16 日和2月20日,大唐河北乌登山风场和山西偏关水泉风场接连发生风机倒塔,造成设别损坏 的事件。河北乌登山风场110号风机倒塔事故的直接原因是叶片质量出现问题,在运行中开 裂,气流不平稳而引起风机的剧烈摇晃,山西偏关水泉风场项目14号风机倒塔事故的直接原 因则是风机振动值严重超标,造成法兰疲劳开裂导致的风机倒塔。
现阶段,暂时没有人针对风机运行变工况的特点,以及数据的静态与动态特性,来研究 风机的状态监测,判断其结冰故障。当前对风机结冰等故障的判断,主要依靠风电操作员的 经验,或者通过二维数据特征,通常为即风电功率曲线,来简单进行判断。但是这样做出的 结冰故障判断通常只能在严重结冰判断出来,无法及时在结冰早期就做出判断。这可能是由 于先前的方法无法有效利用风机中大量其他测点带来的完整信息,且没有考虑风机运行变工 况的特点,以及数据的静态与动态特性,从而无法对风机的整体状态形成完成、准确、及时 的监测。
发明内容
本发明的目的在于针对现有技术的不足,提供一种基于细粒度风电发电状态划分的风机 叶片结冰异常检测方法。
本发明的目的通过以下技术方案实现:一种基于细粒度风电发电状态划分的风机叶片结 冰异常检测方法,本方法的基本分析和建模单元是风速片数据矩阵,输出是沿速度方向的状 态分类,并且不同的子状态沿速度方向是连续的。通过顺序地添加新的风速片以迭代更新可 变展开建模单元,基于其与参考分布相似度的变化来检查所得模型的相似特性的变化,若特 性变化超过设定阈值,则认为不属于同一子状态。对于得出的各状态分类,可根据其统计监 测量的数据分布得监测置信限。该算法实际上是通过考察模型的改变其对监测性能的影响来 识别不同状态的。该方法包括两个部分:
(1)基于细粒度风电发电状态划分获得离线子状态,包括以下子步骤:
(1.1)获得风场正常工作的I条历史数据并将其描述为一个二维矩阵X(I×J)。矩阵中的J个 可测变量为运行过程中被测量到的J个状态参数。
(1.2)准备风速片数据矩阵:将步骤(1.1)描述的二维矩阵X(I×J)中的风速变量为参考进 行离散化处理,得到K个原始风速片。每个风速片中有Ik(k=1,2,...,K)个样本以及J个可测 变量,其中下标k表示原始风速片的序号,并有离散化后第k个原始风速片记为 Xk。对每一个原始风速片Xk进行各自单独的数据标准化处理,得到标准化后的风速片
其中,表示Xk的均值,D(Xk)表示Xk的方差。
(1.3)对风速片执行SFA算法并获得初始风速片模型:对的协方差矩阵进行奇异值分 解得到标准正交基所张成的矩阵Uk:
其中,Λk是以对的协方差矩阵进行奇异值分解得到的所有标准化特征值为主对角线 的矩阵。
进一步得到白化矩阵Qk:Qk=Λk -1/2Uk T
再对风速片进行白化,得到白化后的风速片Zk:
对Zk的一阶差分的协方差矩阵进行奇异值分解得到标准正交基所张成的矩阵Pk:
其中,Ωk是以对的协方差矩阵进行奇异值分解得到的所有标准化特征值为主对角线 的矩阵。
进一步得投影矩阵Wk:Wk=PkQk
得到的Ωk中标准化特征值数量为N,将特征值由小到大排序,保留前M个特征,得Ωk,M, 剩余的为Ωk,Me,Me=N-M。对应的,保留按特征值大小排序的M个标准正交基,张成矩阵Pk,M。剩余的标准正交基张成矩阵Pk,Me。由于各风速片变量维度为J,其特征值数量相同,N=J。
M决定方法如下:统计所有风速片低于平均变换速度的成分数量,采用出现次数最多的 成分的个数作为M。
把Pk,M、Pk,Me、Ωk,M、Ωk,Me作为监测模型,将风速片投影到特征空间,得风速片在特征空间的投影Sk,分为两部分计算:
Sk=[Sk,d Sk,e]
Sk,d=[Sk,1,Sk,2,...,Sk,M]=Pk,MZk
Sk,e=[Sk,1,Sk,2,...,Sk,Me]=Pk,MeZk
其中,Sk,d包含M个慢特征[Sk,1,Sk,2,...,Sk,M]。Sk,e包含Me个剩余特征[Sk,1,Sk,2,...,Sk,Me]。
计算Sk,d中每条向量的T2监测统计值
其中,表示第k个风速片下第j条数据的T2监测统计值,sk,d,j表示Sk,d的第j条向量。
统计离散化概率分布Prk,j,计算Sk,d的监测统计值与标准卡方分布的相似度Simk:
在标准卡方分布x轴上的与之间,取n个离散化段,n=60;统计 在个离散化段中出现的概率Prk,j(j=1,2,...,n),得其概率分布,并统计标准卡方分布在 离散化段的概率分布为Prk,jr(j=1,2,...,n)。
Simk的计算公式如下:
(1.4)将作为状态片Yk,nu的初始风速片:Yk,nu表示为Yk,nu(Iu×J),由当前 状态片中包含的风速片展开而成;Iu表示第u个风速片的样本数量;k为当前合成风速片中 初始风速片序数,最初等于1;nu为当前状态片内除初始风速片外的风速片个数,初始时, nu=0。
向状态片Yk,nu中逐一添加风速片:将nu加1,即nu=nu+1,把风速片添加到状态 片Yk,nu中,得到状态片Yk,nu(Iu×J),将状态片Yk,nu按照(1.3)逐步计算得到状 态片的Zk,nu、Ωk,nu,M、Ωk,nu,Me、Pk,nu,M和Pk,nu,Me,将Pk,nu,M、Pk,nu,Me作为状态片的监测模型,并代入状态片Yk,nu中的所有风速片中,计算获得每个风速片在新状态片Yk,nu中Sk,nu,u(u=0,1,...,nu),其中, Sk,nu,u=[Sk,nu,u,dSk,nu,u,e](u=0,1,...,nu),Sk,nu,u,d=Pk,nu,MZk+u,Sk,nu,u,e=Pk,nu,MeZk+u;Sk,nu,u表示使 用状态片Yk,nu当前的监测模型对第k+u个风速片在特征空间的投影,表示第k+u个 风速片,Zk+u表示第k+u个风速片经过如步骤(1.3)中所述白化处理后得到的矩阵。
进一步计算Sk,nu,u,d(u=0,1,...,nu)每条向量的监测统计值统计计算获得每个风速片的监测统计值在状态片Yk,nu中 与标准卡方分布的相似度Simk,nu,u(u=0,1,...,nu)。
(1.5)将每个风速片在状态片Yk,nu中的Simk,nu,u(u=0,1,...,nu)与初始风速片 的Simk进行比较:若有连续三个u满足则把当前的nu作为nu*, 继续执行步骤(1.6);若不满足该条件,则重复执行步骤(1.4)。
其中,α是一个附加于Simk的常数,称为松弛因子。
(1.6)将序号为k至k+nu*-1的状态片作为一个离线子状态;假设目前得到的是第a个离线 子状态Ya=Yk,nu*-1,则该离线子状态的监测模型为Pa,M=Pv,nu*-1,M和Pa,Me=Pv,nu*-1,Me,且该子状 态的Za=Zk,nu*-1、Ωa,M=Ωk,nu*-1,M、Ωa,Me=Ωk,nu*-1,Me,进一步得Sa,d=Pa,MZa,Sa,e=Pa,MeZa。 计算该离线子状态Yk,nu*-1的四个监测统计量并根据高斯核密度估计确定 这四个监测统计量对应的置信限(置信度α一般取0.95)ST1a,ST2a,DY1a,DY2a。
四个监测统计量的计算与服从分布如下:
Ta 2=sa,d,i Tsa,d,i,服从分布:
Ta,e 2=sa,e,i Tsa,e,i,服从分布:
Sa 2=sa,d,i TΩa,M -1sa,d,i,服从分布:
Sa,e 2=sa,e,i TΩa,Me -1sa,e,i,服从分布:
其中,sa,d,i为Sa,d中的第i条向量,sa,e,i为Sa,e中的第i条向量。表示自由度为p的卡方 分布,Fp,q表示自由度为p,q的方差比分布。
(1.7)将风速片作为新状态片的初始风速片,进行下一个离线子状态划分;新状态片 Yk,nu中:k=k+nu*,nu归零。递归地重复步骤(1.2)~(1.6)获得全部离线子状态。
(2)获取风场工作实时监测数据并利用步骤(1)状态划分的离线子状态对风场工作实时数 据划分状态,获得风机实时运行状态;具体包括如下步骤:
(2.1)获得数据并计算监测统计量:获取风场工作实时监测数据,包括当前时间及其之前的 Win条时序数据,每条数据包括J个状态参数。将监测数据标准化后根据风速划分成与步骤(1.7) 中获得的离线子状态相对应的λ个实时监测子状态,形成实时监测状态片Ymv(v=1,2,...,λ)。 根据步骤(1.3)中的相关公式计算得到Ymv对应的Zmv(v=1,2,...,λ),并根据各自对应的离线 子状态的Pa,M、Pa,Me、Ωa,M、Ωa,Me,计算得到Smv,d=Pa,MZmv、Smv,e=Pa,MeZmv,进一步按 照步骤(1.5)的相关公式计算得到当前实时监测状态片Ymv(v=1,2,...,λ)的Tv 2(v=1,2,...,λ)、 Tv,e 2(v=1,2,...,λ)、Sv 2(v=1,2,...,λ)、Sv,e 2(v=1,2,...,λ)监测统计量。
(2.2)比较实时监测统计量与离线监测阈值:根据对应的离线子状态中的控制限ST1a,ST2a, DY1a,DY2a判断实时监测子状态的监测统计量Tv 2(v=1,2,...,λ)、Tv,e 2(v=1,2,...,λ)、 Sv 2(v=1,2,...,λ)、Sv,e 2(v=1,2,...,λ)中每个数据是否超限,分别统计本次监测中各监测统计量 对应的超限个数TC_N、TeC_N、SC_N、SeC_N,及超限比例TC、TeC、SC、SeC, 其中TC、TeC作为静态指标,SC、SeC作为动态指标。
(2.3)判断风机实时运行状态:当指标超过预设的比例,这里采用10%,即判定指标异常。 如果动态指标正常,静态指标异常,则判断风机工况发生变动,但工作状态正常,判断其状 态为State1;如果动态指标正常,静态指标正常,则判断风机工况不变动,且工作状态正常, 判断其状态为State2;如果动态指标异常,静态指标也异常,则判断风机工况发生变动,且 工作状态异常,判断其状态为State3;如果动态指标异常,静态指标正常,则判断风机工况 不变动,但工作状态异常,判断其状态为State4。
本发明的有益效果在于:本发明方法利用风速作为主要关注变量,在揭示风机变工况运 行时的各变量关系分段性的特点和规律后,通过不同风速将风电数据划分为不同状态,详细 的划分了风力发电机组在不同风速下的数据特性,然后根据得到的不同的数据特性,建立不 同状态下的监测模型,提高了风电数据利用效率。该方法创新的提出了以风速为划分依据且 后续动静态协同监控的想法。有效利用了风机中大量其他测点带来的完整信息,且考虑风机 运行变工况的特点,以及数据的静态与动态特性,细化的区分了不同风速下的风机运行情况, 并动静协同,从而对风机的整体状态形成完整、准确、及时的监测。从而保证了风力发电机 组的安全可靠的运行。
附图说明:
图1是一种基于细粒度风电发电状态划分的风机叶片结冰异常检测方法的流程图,(a) 为离线建模过程流程图,(b)为在线监测流程图;
图2是本发明方法的细粒度状态划分展示图,为功率曲线的分类展示;
图3是本发明方法的在线监测图示,(a)正常状态段监测示意图,(b)异常状态段监测 示意图。
具体实施方式
针对风机有大量测点、变工况运行、数据的静态动态特性不一的情况,提出了一种细粒 度风电发电状态划分和风机叶片结冰异常监测方法。该方法在揭示风机变工况运行时的各变 量关系的分段性特点后,针对性地提出了一种基于慢特征提取的细粒度风电发电状态划分和 动静协同风机叶片结冰异常监测方法,实现离线的风机出力性能评价与刻画。其特点在于把 时序数据转化为风速片数据,然后提取慢特征进行状态划分建模。针对风力发电的动态变化 特点,提出了一种动静态协同的监测方法,对风机运行状态进行在线监测,通过判断风机当 前运行异常的情况,可以描述风机出力性能下降的情况。其特点在于监测时将数据对应到之 前划分的子模型中去,同时考虑静态与动态监测统计量来进行监测。下面结合附图及具体实 例,对本发明作进一步详细说明。
本发明以开源数据集,天池大赛的风机结冰故障预测数据集为例。包括如下变量:风速 (wind_speed)、发电机转速(generator_speed)、网侧有功功率/kw(power)、对风角/°(wind_direction)、25秒平均风向角(wind_direction_mean)、偏航位置(yaw_position)、偏 航速度(yaw_speed)、叶片1角度(pitch1_angle)、叶片2角度(pitch2_angle)、叶片3角度 (pitch3_angle)、变桨电机1温度(pitch1_moto_tmp)、变桨电机2温度(pitch2_moto_tmp)、 变桨电机3温度(pitch3_moto_tmp)、x方向加速度(acc_x)、y方向加速度(acc_y)、环境 温度(environment_tmp)、机舱温度(int_tmp)、ng5 1温度(pitch1_ng5_tmp)、ng5 2温度 (pitch2_ng5_tmp)、ng5 3温度(pitch3_ng5_tmp)、ng5 1充电器直流电流(pitch1_ng5_DC)、 ng5 2充电器直流电流(pitch2_ng5_DC)、ng5 3充电器直流电流(pitch3_ng5_DC)。
由先验知识,结合皮尔逊相关系数,可以判定pitch1_angle叶片1角度、pitch2_angle叶 片2角度、pitch3_angle叶片3角度是强相关的,可以保留其均值为pitch_angle_mean,并删 去原来的这三个变量。
pitch1_speed叶片1速度、pitch2_speed叶片2速度、pitch3_speed叶片3速度也为强相关, 可以保留其均值为pitch_speed_mean,并删去原来的这三个变量。
pitch1_moto_tmp变桨电机1温度、pitch2_moto_tmp变桨电机2温度、pitch3_moto_tmp 变桨电机3温度也为强相关,可以保留其均值为pitch_moto_tmp_mean,并删去原来的这三个 变量。
共有20个变量。因此,本实施例中提到J为20,历史数据I为350255。
应该理解,本发明不止局限于上述实例的风力发电过程,凡是熟悉本领域的技术人员在 不违背本发明的前提下还可以做出等同变型或替换,这些等同的变型或替换均包含在本申请 权利要求所限定的范围内。
如图1所示,本发明是一种基于细粒度风电发电状态划分的风机叶片结冰异常检测方法, 包括以下两个部分:
本发明的目的通过以下技术方案实现:一种基于细粒度风电发电状态划分的风机叶片结 冰异常检测方法,该方法包括两个部分:
(1)基于细粒度风电发电状态划分获得离线子状态,包括以下子步骤:
(1.1)获得风场正常工作的I条历史数据并将其描述为一个二维矩阵X(I×J)。矩阵中的J个 可测变量为运行过程中被测量到的J个状态参数。
(1.2)准备风速片数据矩阵:将步骤(1.1)描述的二维矩阵X(I×J)中的风速变量为参考进 行离散化处理,得到K个原始风速片。每个风速片中有Ik(k=1,2,...,K)个样本以及J个可测 变量,其中下标k表示原始风速片的序号,并有离散化后第k个原始风速片记为 Xk。对每一个原始风速片Xk进行各自单独的数据标准化处理,得到标准化后的风速片
其中,表示Xk的均值,D(Xk)表示Xk的方差。
(1.3)对风速片执行SFA算法并获得初始风速片模型:对的协方差矩阵进行奇异值分 解得到标准正交基所张成的矩阵Uk:
其中,Λk是以对的协方差矩阵进行奇异值分解得到的所有标准化特征值为主对角线 的矩阵。
进一步得到白化矩阵Qk:Qk=Λk -1/2Uk T
再对风速片进行白化,得到白化后的风速片Zk:
对Zk的一阶差分的协方差矩阵进行奇异值分解得到标准正交基所张成的矩阵Pk:
其中,Ωk是以对的协方差矩阵进行奇异值分解得到的所有标准化特征值为主对角线 的矩阵。
进一步得投影矩阵Wk:Wk=PkQk
得到的Ωk中标准化特征值数量为N,将特征值由小到大排序,保留前M个特征,得Ωk,M, 剩余的为Ωk,Me,Me=N-M。对应的,保留按特征值大小排序的M个标准正交基,张成矩阵Pk,M。剩余的标准正交基张成矩阵Pk,Me。由于各风速片变量维度为J,其特征值数量相同,N=J。
M决定方法如下:统计所有风速片低于平均变换速度的成分数量,采用出现次数最多的 成分的个数作为M。
把Pk,M、Pk,Me、Ωk,M、Ωk,Me作为监测模型,将风速片投影到特征空间,得风速片在特征空间的投影Sk,分为两部分计算:
Sk=[Sk,d Sk,e]
Sk,d=[Sk,1,Sk,2,...,Sk,M]=Pk,MZk
Sk,e=[Sk,1,Sk,2,...,Sk,Me]=Pk,MeZk
其中,Sk,d包含M个慢特征[Sk,1,Sk,2,...,Sk,M]。Sk,e包含Me个剩余特征[Sk,1,Sk,2,...,Sk,Me]。
计算Sk,d中每条向量的T2监测统计值
其中,表示第k个风速片下第j条数据的T2监测统计值,sk,d,j表示Sk,d的第j条向量。
统计离散化概率分布Prk,j,计算Sk,d的监测统计值与标准卡方分布的相似度Simk:
在标准卡方分布x轴上的与之间,取n个离散化段,n=60;统计 在个离散化段中出现的概率Prk,j(j=1,2,...,n),得其概率分布,并统计标准卡方分布在 离散化段的概率分布为Prk,jr(j=1,2,...,n)。
Simk的计算公式如下:
(1.4)将作为状态片Yk,nu的初始风速片:Yk,nu表示为Yk,nu(Iu×J),由当前 状态片中包含的风速片展开而成;Iu表示第u个风速片的样本数量;k为当前合成风速片中 初始风速片序数,最初等于1;nu为当前状态片内除初始风速片外的风速片个数,初始时, nu=0。
向状态片Yk,nu中逐一添加风速片:将nu加1,即nu=nu+1,把风速片添加到状态 片Yk,nu中,得到状态片Yk,nu(Iu×J),将状态片Yk,nu按照(1.3)逐步计算得到状 态片的Zk,nu、Ωk,nu,M、Ωk,nu,Me、Pk,nu,M和Pk,nu,Me,将Pk,nu,M、Pk,nu,Me作为状态片的监测模型,并代入状态片Yk,nu中的所有风速片中,计算获得每个风速片在新状态片Yk,nu中Sk,nu,u(u=0,1,...,nu),其中, Sk,nu,u=[Sk,nu,u,dSk,nu,u,e](u=0,1,...,nu),Sk,nu,u,d=Pk,nu,MZk+u,Sk,nu,u,e=Pk,nu,MeZk+u;Sk,nu,u表示使 用状态片Yk,nu当前的监测模型对第k+u个风速片在特征空间的投影,表示第k+u个 风速片,Zk+u表示第k+u个风速片经过如步骤(1.3)中所述白化处理后得到的矩阵。
进一步计算Sk,nu,u,d(u=0,1,...,nu)每条向量的监测统计值统计计算获得每个风速片的监测统计值在状态片Yk,nu中 与标准卡方分布的相似度Simk,nu,u(u=0,1,...,nu)。
(1.5)将每个风速片在状态片Yk,nu中的Simk,nu,u(u=0,1,...,nu)与初始风速片 的Simk进行比较:若有连续三个u满足则把当前的nu作为nu*, 继续执行步骤(1.6);若不满足该条件,则重复执行步骤(1.4)。
其中,α是一个附加于Simk的常数,称为松弛因子。
(1.6)将序号为k至k+nu*-1的状态片作为一个离线子状态;假设目前得到的是第a个离线 子状态Ya=Yk,nu*-1,则该离线子状态的监测模型为Pa,M=Pv,nu*-1,M和Pa,Me=Pv,nu*-1,Me,且该子状 态的Za=Zk,nu*-1、Ωa,M=Ωk,nu*-1,M、Ωa,Me=Ωk,nu*-1,Me,进一步得Sa,d=Pa,MZa,Sa,e=Pa,MeZa。 计算该离线子状态Yk,nu*-1的四个监测统计量并根据高斯核密度估计确定 这四个监测统计量对应的置信限(置信度α一般取0.95)ST1a,ST2a,DY1a,DY2a。
四个监测统计量的计算与服从分布如下:
Ta 2=sa,d,i Tsa,d,i,服从分布:
Ta,e 2=sa,e,i Tsa,e,i,服从分布:
Sa 2=sa,d,i TΩa,M -1sa,d,i,服从分布:
Sa,e 2=sa,e,i TΩa,Me -1sa,e,i,服从分布:
其中,sa,d,i为Sa,d中的第i条向量,sa,e,i为Sa,e中的第i条向量。表示自由度为p的卡方 分布,Fp,q表示自由度为p,q的方差比分布。
(1.7)将风速片作为新状态片的初始风速片,进行下一个离线子状态划分;新状态片 Yk,nu中:k=k+nu*,nu归零。递归地重复步骤(1.2)~(1.6)获得全部离线子状态。
本实施例中,分为了7个子状态,如图2所示。
(2)获取风场工作实时监测数据并利用步骤(1)状态划分的离线子状态对风场工作实时数 据划分状态,获得风机实时运行状态;具体包括如下步骤:
(2.1)获得数据并计算监测统计量:获取风场工作实时监测数据,包括当前时间及其之前的 Win条时序数据,每条数据包括J个状态参数。将监测数据标准化后根据风速划分成与步骤(1.7) 中获得的离线子状态相对应的λ个实时监测子状态,形成实时监测状态片Ymv(v=1,2,...,λ)。 根据步骤(1.3)中的相关公式计算得到Ymv对应的Zmv(v=1,2,...,λ),并根据各自对应的离线 子状态的Pa,M、Pa,Me、Ωa,M、Ωa,Me,计算得到Smv,d=Pa,MZmv、Smv,e=Pa,MeZmv,进一步按 照步骤(1.5)的相关公式计算得到当前实时监测状态片Ymv(v=1,2,...,λ)的Tv 2(v=1,2,...,λ)、 Tv,e 2(v=1,2,...,λ)、Sv 2(v=1,2,...,λ)、Sv,e 2(v=1,2,...,λ)监测统计量。
(2.2)比较实时监测统计量与离线监测阈值:根据对应的离线子状态中的控制限ST1a,ST2a, DY1a,DY2a判断实时监测子状态的监测统计量Tv 2(v=1,2,...,λ)、Tv,e 2(v=1,2,...,λ)、 Sv 2(v=1,2,...,λ)、Sv,e 2(v=1,2,...,λ)中每个数据是否超限,分别统计本次监测中各监测统计量 对应的超限个数TC_N、TeC_N、SC_N、SeC_N,及超限比例TC、TeC、SC、SeC, 其中TC、TeC作为静态指标,SC、SeC作为动态指标。
(2.3)判断风机实时运行状态:当指标超过预设的比例,这里采用10%,即判定指标异常。 如果动态指标正常,静态指标异常,则判断风机工况发生变动,但工作状态正常,判断其状 态为State1;如果动态指标正常,静态指标正常,则判断风机工况不变动,且工作状态正常, 判断其状态为State2;如果动态指标异常,静态指标也异常,则判断风机工况发生变动,且 工作状态异常,判断其状态为State3;如果动态指标异常,静态指标正常,则判断风机工况 不变动,但工作状态异常,判断其状态为State4。
(2.4)判断得到结果后,即可指导风机的运行与故障维修。
选取风机正常工作的某段时间段和风机异常工作,即叶片结冰的某段时间段,对监测效 果进行验证。若State3到达高位,对应严重异常,若State4到达高位,对应微小异常。其中 Abnormal到达高位,对应人工标注的结冰异常。如图3(a)所示,正常段监测时会出现少量 误报,在可接受范围之内;如图3(b)所示,异常段监测时,能够提前预测结冰异常的出现, 对风机的整体状态形成了完整、准确、及时的监测。有助于指导风机的故障维修,从而保障 风力发电机组安全可靠运行。
本方法创新的提出了以风速为划分依据且后续动静态协同监控的想法。有效利用了风机 中大量其他测点带来的完整信息,且考虑风机运行变工况的特点,以及数据的静态与动态特 性,细化的区分了不同风速下的风机运行情况,并动静协同,从而对风机的整体状态形成完 整、准确、及时的监测。从而保证了风力发电机组的安全可靠的运行。
Claims (1)
1.基于细粒度风电发电状态划分的风机叶片结冰异常检测方法,其特征在于,包括如下步骤:
(1)基于细粒度风电发电状态划分获得离线子状态,包括以下子步骤:
(1.1)获得风场正常工作的I条历史数据并将其描述为一个二维矩阵X(I×J)。矩阵中的J个可测变量为运行过程中被测量到的J个状态参数。
(1.2)准备风速片数据矩阵:将步骤(1.1)描述的二维矩阵X(I×J)中的风速变量为参考进行离散化处理,得到K个原始风速片。每个风速片中有Ik(k=1,2,...,K)个样本以及J个可测变量,其中下标k表示原始风速片的序号,并有离散化后第k个原始风速片记为Xk。对每一个原始风速片Xk进行各自单独的数据标准化处理,得到标准化后的风速片
其中,表示Xk的均值,D(Xk)表示Xk的方差。
(1.3)对风速片执行SFA算法并获得初始风速片模型:对的协方差矩阵进行奇异值分解得到标准正交基所张成的矩阵Uk:
其中,Λk是以对的协方差矩阵进行奇异值分解得到的所有标准化特征值为主对角线的矩阵。
进一步得到白化矩阵Qk:Qk=Λk -1/2Uk T
再对风速片进行白化,得到白化后的风速片Zk:
对Zk的一阶差分的协方差矩阵进行奇异值分解得到标准正交基所张成的矩阵Pk:
其中,Ωk是以对的协方差矩阵进行奇异值分解得到的所有标准化特征值为主对角线的矩阵。
进一步得投影矩阵Wk:Wk=PkQk
得到的Ωk中标准化特征值数量为N,将特征值由小到大排序,保留前M个特征,得Ωk,M,剩余的为Ωk,Me,Me=N-M。对应的,保留按特征值大小排序的M个标准正交基,张成矩阵Pk,M。剩余的标准正交基张成矩阵Pk,Me。由于各风速片变量维度为J,其特征值数量相同,N=J。
M决定方法如下:统计所有风速片低于平均变换速度的成分数量,采用出现次数最多的成分的个数作为M。
把Pk,M、Pk,Me、Ωk,M、Ωk,Me作为监测模型,将风速片投影到特征空间,得风速片在特征空间的投影Sk,分为两部分计算:
Sk=[Sk,d Sk,e]
Sk,d=[Sk,1,Sk,2,...,Sk,M]=Pk,MZk
Sk,e=[Sk,1,Sk,2,...,Sk,Me]=Pk,MeZk
其中,Sk,d包含M个慢特征[Sk,1,Sk,2,...,Sk,M]。Sk,e包含Me个剩余特征[Sk,1,Sk,2,...,Sk,Me]。
计算Sk,d中每条向量的T2监测统计值
其中,表示第k个风速片下第j条数据的T2监测统计值,sk,d,j表示Sk,d的第j条向量。
统计离散化概率分布Prk,j,计算Sk,d的监测统计值与标准卡方分布的相似度Simk:
在标准卡方分布x轴上的与之间,取n个离散化段,n=60;统计在个离散化段中出现的概率Prk,j(j=1,2,...,n),得其概率分布,并统计标准卡方分布在离散化段的概率分布为Prk,jr(j=1,2,...,n)。
Simk的计算公式如下:
(1.4)将作为状态片Yk,nu的初始风速片:Yk,nu表示为Yk,nu(Iu×J),由当前状态片中包含的风速片展开而成;Iu表示第u个风速片的样本数量;k为当前合成风速片中初始风速片序数,最初等于1;nu为当前状态片内除初始风速片外的风速片个数,初始时,nu=0。
向状态片Yk,nu中逐一添加风速片:将nu加1,即nu=nu+1,把风速片添加到状态片Yk,nu中,得到状态片Yk,nu(Iu×J),将状态片Yk,nu按照(1.3)逐步计算得到状态片的Zk,nu、Ωk,nu,M、Ωk,nu,Me、Pk,nu,M和Pk,nu,Me,将Pk,nu,M、Pk,nu,Me作为状态片的监测模型,并代入状态片Yk,nu中的所有风速片中,计算获得每个风速片在新状态片Yk,nu中Sk,nu,u(u=0,1,...,nu),其中,Sk,nu,u=[Sk,nu,u,dSk,nu,u,e](u=0,1,...,nu),Sk,nu,u,d=Pk,nu,MZk+u,Sk,nu,u,e=Pk,nu,MeZk+u;Sk,nu,u表示使用状态片Yk,nu当前的监测模型对第k+u个风速片在特征空间的投影,表示第k+u个风速片,Zk+u表示第k+u个风速片经过如步骤(1.3)中所述白化处理后得到的矩阵。
进一步计算Sk,nu,u,d(u=0,1,...,nu)每条向量的监测统计值统计计算获得每个风速片的监测统计值在状态片Yk,nu中与标准卡方分布的相似度Simk,nu,u(u=0,1,...,nu)。
(1.5)将每个风速片在状态片Yk,nu中的Simk,nu,u(u=0,1,...,nu)与初始风速片的Simk进行比较:若有连续三个u满足则把当前的nu作为nu*,继续执行步骤(1.6);若不满足该条件,则重复执行步骤(1.4)。
其中,α是一个附加于Simk的常数,称为松弛因子。
(1.6)将序号为k至k+nu*-1的状态片作为一个离线子状态;假设目前得到的是第a个离线子状态则该离线子状态的监测模型为和且该子状态的进一步得Sa,d=Pa,MZa,Sa,e=Pa,MeZa。计算该离线子状态的四个监测统计量并根据高斯核密度估计确定这四个监测统计量对应的置信限(置信度α一般取0.95)ST1a,ST2a,DY1a,DY2a。
四个监测统计量的计算与服从分布如下:
Ta 2=sa,d,i Tsa,d,i,服从分布:
Ta,e 2=sa,e,i Tsa,e,i,服从分布:
Sa 2=sa,d,i TΩa,M -1sa,d,i,服从分布:
Sa,e 2=sa,e,i TΩa,Me -1sa,e,i,服从分布:
其中,sa,d,i为Sa,d中的第i条向量,sa,e,i为Sa,e中的第i条向量。表示自由度为p的卡方分布,Fp,q表示自由度为p,q的方差比分布。
(1.7)将风速片作为新状态片的初始风速片,进行下一个离线子状态划分;新状态片Yk,nu中:k=k+nu*,nu归零。递归地重复步骤(1.2)~(1.6)获得全部离线子状态。
(2)获取风场工作实时监测数据并利用步骤(1)状态划分的离线子状态对风场工作实时数据划分状态,获得风机实时运行状态;具体包括如下步骤:
(2.1)获得数据并计算监测统计量:获取风场工作实时监测数据,包括当前时间及其之前的Win条时序数据,每条数据包括J个状态参数。将监测数据标准化后根据风速划分成与步骤(1.7)中获得的离线子状态相对应的λ个实时监测子状态,形成实时监测状态片Ymv(v=1,2,...,λ)。根据步骤(1.3)中的相关公式计算得到Ymv对应的Zmv(v=1,2,...,λ),并根据各自对应的离线子状态的Pa,M、Pa,Me、Ωa,M、Ωa,Me,计算得到Smv,d=Pa,MZmv、Smv,e=Pa,MeZmv,进一步按照步骤(1.5)的相关公式计算得到当前实时监测状态片Ymv(v=1,2,...,λ)的Tv 2(v=1,2,...,λ)、Tv,e 2(v=1,2,...,λ)、Sv 2(v=1,2,...,λ)、Sv,e 2(v=1,2,...,λ)监测统计量。
(2.2)比较实时监测统计量与离线监测阈值:根据对应的离线子状态中的控制限ST1a,ST2a,DY1a,DY2a判断实时监测子状态的监测统计量Tv 2(v=1,2,...,λ)、Tv,e 2(v=1,2,...,λ)、Sv 2(v=1,2,...,λ)、Sv,e 2(v=1,2,...,λ)中每个数据是否超限,分别统计本次监测中各监测统计量对应的超限个数TC_N、TeC_N、SC_N、SeC_N,及超限比例TC、TeC、SC、SeC,其中TC、TeC作为静态指标,SC、SeC作为动态指标。
(2.3)判断风机实时运行状态:当指标超过预设的比例,这里采用10%,即判定指标异常。如果动态指标正常,静态指标异常,则判断风机工况发生变动,但工作状态正常,判断其状态为State1;如果动态指标正常,静态指标正常,则判断风机工况不变动,且工作状态正常,判断其状态为State2;如果动态指标异常,静态指标也异常,则判断风机工况发生变动,且工作状态异常,判断其状态为State3;如果动态指标异常,静态指标正常,则判断风机工况不变动,但工作状态异常,判断其状态为State4。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201910448367.2A CN110222393B (zh) | 2019-05-28 | 2019-05-28 | 基于细粒度风电发电状态划分的风机叶片结冰异常监测方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201910448367.2A CN110222393B (zh) | 2019-05-28 | 2019-05-28 | 基于细粒度风电发电状态划分的风机叶片结冰异常监测方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN110222393A true CN110222393A (zh) | 2019-09-10 |
CN110222393B CN110222393B (zh) | 2020-12-18 |
Family
ID=67818163
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201910448367.2A Active CN110222393B (zh) | 2019-05-28 | 2019-05-28 | 基于细粒度风电发电状态划分的风机叶片结冰异常监测方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN110222393B (zh) |
Cited By (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN111524502A (zh) * | 2020-05-27 | 2020-08-11 | 科大讯飞股份有限公司 | 一种语种检测方法、装置、设备及存储介质 |
CN112651648A (zh) * | 2020-12-30 | 2021-04-13 | 广东电网有限责任公司电力科学研究院 | 一种基于数据白化的配电网可靠性指标预处理方法 |
CN114753980A (zh) * | 2022-04-29 | 2022-07-15 | 南京国电南自维美德自动化有限公司 | 一种风机叶片结冰监测方法及系统 |
Citations (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN104537260A (zh) * | 2015-01-14 | 2015-04-22 | 清华大学 | 基于缓慢特征回归的动态软测量方法和系统 |
WO2016093373A1 (ja) * | 2016-02-01 | 2016-06-16 | 株式会社小松製作所 | 作業機械の制御システム、作業機械、及び作業機械の管理システム |
CN106647718A (zh) * | 2017-01-20 | 2017-05-10 | 中国石油大学(华东) | 基于贝叶斯核慢特征分析的非线性工业过程故障检测方法 |
CN108803531A (zh) * | 2018-07-17 | 2018-11-13 | 浙江大学 | 基于动静特征协同分析和有序时段划分的闭环系统过程监测方法 |
CN109143995A (zh) * | 2018-07-13 | 2019-01-04 | 浙江大学 | 一种基于质量相关慢特征充分分解的闭环系统精细运行状态监测方法 |
CN109189020A (zh) * | 2018-09-11 | 2019-01-11 | 浙江大学 | 一种基于动静态特性协同分析的大型燃煤机组燃烧系统综合监测方法 |
CN109471420A (zh) * | 2018-09-21 | 2019-03-15 | 浙江大学 | 基于cva-sfa的智能电厂大型燃煤发电机组空气预热器控制性能监测方法 |
-
2019
- 2019-05-28 CN CN201910448367.2A patent/CN110222393B/zh active Active
Patent Citations (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN104537260A (zh) * | 2015-01-14 | 2015-04-22 | 清华大学 | 基于缓慢特征回归的动态软测量方法和系统 |
WO2016093373A1 (ja) * | 2016-02-01 | 2016-06-16 | 株式会社小松製作所 | 作業機械の制御システム、作業機械、及び作業機械の管理システム |
CN106647718A (zh) * | 2017-01-20 | 2017-05-10 | 中国石油大学(华东) | 基于贝叶斯核慢特征分析的非线性工业过程故障检测方法 |
CN109143995A (zh) * | 2018-07-13 | 2019-01-04 | 浙江大学 | 一种基于质量相关慢特征充分分解的闭环系统精细运行状态监测方法 |
CN108803531A (zh) * | 2018-07-17 | 2018-11-13 | 浙江大学 | 基于动静特征协同分析和有序时段划分的闭环系统过程监测方法 |
CN109189020A (zh) * | 2018-09-11 | 2019-01-11 | 浙江大学 | 一种基于动静态特性协同分析的大型燃煤机组燃烧系统综合监测方法 |
CN109471420A (zh) * | 2018-09-21 | 2019-03-15 | 浙江大学 | 基于cva-sfa的智能电厂大型燃煤发电机组空气预热器控制性能监测方法 |
Non-Patent Citations (5)
Title |
---|
CHUNHUI ZHAO ET AL.: ""A sparse dissimilarity analysis algorithm for incipient fault isolation with no priori fault information"", 《CONTROL ENGINEERING PRACTICE》 * |
WEI QIAO ET AL.: ""A Survey on Wind Turbine Condition Monitoring and Fault Diagnosis−Part II Signals and Signal Processing Methods"", 《IEEE TRANSACTIONS ON INDUSTRIAL ELECTRONICS》 * |
YAN QIN ET AL.: ""Comprehensive process decomposition for closed-loop process monitoring with quality-relevant slow feature analysis"", 《JOURNAL OF PROCESS CONTROL》 * |
张启亮: ""基于SCADA数据特征的风电机组叶桨结冰辨识研究"", 《中国优秀硕士学位论文全文数据库工程科技辑》 * |
王伟 等: ""基于PCA的卷烟制丝过程监测与故障诊断"", 《控制工程》 * |
Cited By (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN111524502A (zh) * | 2020-05-27 | 2020-08-11 | 科大讯飞股份有限公司 | 一种语种检测方法、装置、设备及存储介质 |
CN111524502B (zh) * | 2020-05-27 | 2024-04-30 | 科大讯飞股份有限公司 | 一种语种检测方法、装置、设备及存储介质 |
CN112651648A (zh) * | 2020-12-30 | 2021-04-13 | 广东电网有限责任公司电力科学研究院 | 一种基于数据白化的配电网可靠性指标预处理方法 |
CN112651648B (zh) * | 2020-12-30 | 2023-06-06 | 广东电网有限责任公司电力科学研究院 | 一种基于数据白化的配电网可靠性指标预处理方法 |
CN114753980A (zh) * | 2022-04-29 | 2022-07-15 | 南京国电南自维美德自动化有限公司 | 一种风机叶片结冰监测方法及系统 |
CN114753980B (zh) * | 2022-04-29 | 2024-06-04 | 南京国电南自维美德自动化有限公司 | 一种风机叶片结冰监测方法及系统 |
Also Published As
Publication number | Publication date |
---|---|
CN110222393B (zh) | 2020-12-18 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN110222393B (zh) | 基于细粒度风电发电状态划分的风机叶片结冰异常监测方法 | |
CN105134510A (zh) | 一种风力发电机组变桨系统的状态监测和故障诊断方法 | |
Stephen et al. | A copula model of wind turbine performance | |
CN109751206B (zh) | 风机叶片结冰故障预测方法、装置及存储介质 | |
CN110905732B (zh) | 风电机组风轮不平衡的辨识方法、系统及储存介质 | |
CN110761958B (zh) | 风力发电机组的叶片失速诊断方法及装置 | |
CN103925155A (zh) | 一种风电机组输出功率异常的自适应检测方法 | |
CN109184821B (zh) | 一种大型发电机组汽轮机的闭环信息分析的在线监测方法 | |
CN109595130A (zh) | 一种风机叶片结冰故障预测方法及系统 | |
CN108760037B (zh) | 一种基于频谱分析的风力发电机叶片结构损伤检测方法 | |
Wu et al. | Variable forgetting factor identification algorithm for fault diagnosis of wind turbines | |
CN110608133B (zh) | 一种海上风力发电控制系统及方法 | |
CN111794921B (zh) | 一种基于迁移成分分析的陆上风电机组叶片结冰诊断方法 | |
CN115450864A (zh) | 一种基于合成少数类样本的风电机组叶片结冰诊断方法 | |
CN114781259A (zh) | 一种风力发电机诊断方法及系统 | |
CN114542402A (zh) | 一种基于多参量分析的风电叶片故障类型在线诊断方法及系统 | |
Peng et al. | Short-Term Prediction of Generator Blade Ice Fault Based on Multi-AN | |
CN113848347A (zh) | 一种风力发电机测风仪健康状态检测方法 | |
CN108988319B (zh) | 一种基于深度前馈神经网络和数值积分灵敏度的快速紧急控制方法 | |
Bi et al. | Applying instantaneous SCADA data to artificial intelligence based power curve monitoring and WTG fault forecasting | |
Kadali et al. | Evaluation of Energy in Wind Turbine System Using Probability Distribution | |
Liu et al. | Fuzzy fmea of floating wind turbine based on related weights and topsis theory | |
Qin et al. | The Rapid Computational Model for Power Output of Ice-covered Wind Turbines Based on Airfoil Aerodynamic Parameters | |
Ye et al. | Using SCADA data fusion by swarm intelligence for wind turbine condition monitoring | |
CN111219287A (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 |