CN107425548B - 一种插值h∞扩展卡尔曼滤波发电机动态状态估计方法 - Google Patents

一种插值h∞扩展卡尔曼滤波发电机动态状态估计方法 Download PDF

Info

Publication number
CN107425548B
CN107425548B CN201710811192.8A CN201710811192A CN107425548B CN 107425548 B CN107425548 B CN 107425548B CN 201710811192 A CN201710811192 A CN 201710811192A CN 107425548 B CN107425548 B CN 107425548B
Authority
CN
China
Prior art keywords
state
matrix
interpolation
extended kalman
value
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
Application number
CN201710811192.8A
Other languages
English (en)
Other versions
CN107425548A (zh
Inventor
孙永辉
艾蔓桐
王�义
张博文
王加强
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Hohai University HHU
Original Assignee
Hohai University HHU
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Hohai University HHU filed Critical Hohai University HHU
Priority to CN201710811192.8A priority Critical patent/CN107425548B/zh
Publication of CN107425548A publication Critical patent/CN107425548A/zh
Application granted granted Critical
Publication of CN107425548B publication Critical patent/CN107425548B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • HELECTRICITY
    • H02GENERATION; CONVERSION OR DISTRIBUTION OF ELECTRIC POWER
    • H02JCIRCUIT ARRANGEMENTS OR SYSTEMS FOR SUPPLYING OR DISTRIBUTING ELECTRIC POWER; SYSTEMS FOR STORING ELECTRIC ENERGY
    • H02J3/00Circuit arrangements for ac mains or ac distribution networks
    • H02J3/38Arrangements for parallely feeding a single network by two or more generators, converters or transformers
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01RMEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
    • G01R31/00Arrangements for testing electric properties; Arrangements for locating electric faults; Arrangements for electrical testing characterised by what is being tested not provided for elsewhere
    • G01R31/34Testing dynamo-electric machines
    • G01R31/343Testing dynamo-electric machines in operation

Abstract

本发明公开了一种插值H扩展卡尔曼滤波(interpolation H extended Kalman filter,IHEKF)发电机动态状态估计方法,该方法分为自适应插值、多步预测和修正两部分,首先,通过在多步预测和修正之前计算状态方程和量测方程的非线性指标,使用有限状态机确定内插因子,再根据内插因子使用插值法在两个实际量测值之间增加伪量测值;然后根据所增加的伪量测值进行多步预测和修正,运用扩展卡尔曼滤波的预测步得到状态预测值和状态预测误差协方差,并进一步在扩展卡尔曼滤波的修正步引入H对预测值进行修正得到机电暂态过程中发电机功角和电角速度的估计值和估计误差协方差。最后,算例分析结果表明,本发明所提方法可以准确应用于发电机动态状态估计,且具有较好的鲁棒性。

Description

一种插值H∞扩展卡尔曼滤波发电机动态状态估计方法
技术领域
本发明涉及一种插值H扩展卡尔曼滤波发电机动态状态估计方法,属于电力系统监测、分析和控制技术领域。
背景技术
电力系统动态状态估计对于监测和控制电力系统稳定性是非常重要的。动态状态估计能提供发电机、负载和控制器的预测以此来设计电力系统稳定器和稳压器,进而维持和增进暂态稳定性。随着相量测量单元(phasor measurement unit,PMU)的发展,精确跟踪电力系统机电暂态成为了可能。动态状态估计可以利用PMU测得的数据进行滤波并预测系统未来可能的变化,以此来制定相应的控制策略维持电网的安全稳定运行。
卡尔曼滤波作为状态估计一种有效方法,已经被广泛应用于线性系统中。对非线性系统来说,应用最多的是扩展卡尔曼滤波,在众多非线性系统状态估计的应用中,扩展卡尔曼滤波取得了较好的效果。但是,由于传统扩展卡尔曼滤波使用一阶泰勒展开式逼近非线性函数,在高度非线性系统中,估计效果会不理想,并且对于干扰的鲁棒性也较差。
本发明提出的一种插值H扩展卡尔曼滤波发电机动态状态估计方法,不仅减轻了非线性在估计精度上的消极影响,还增强了滤波对噪声的鲁棒性能。
发明内容
发明目的:本发明所要解决的技术问题是针对扩展卡尔曼滤波对于高度非线性系统滤波精度低鲁棒性差的不足,提出一种插值H扩展卡尔曼滤波发电机动态状态估计方法。
技术方案:一种插值H扩展卡尔曼滤波发电机动态状态估计方法,依次按以下步骤实现的:
(1)获得所需估计发电机机组的参数信息;
(2)程序初始化;
(3)利用H扩展卡尔曼滤波进行一次滤波;
(4)引入非线性指标:根据下式计算状态方程和量测方程的非线性指标:
Figure GDA0002427083220000021
Figure GDA0002427083220000022
Figure GDA0002427083220000023
Figure GDA0002427083220000024
式中,f(·)和h(·)分别为状态转换函数和量测函数,xk为k时刻的状态变量,δx为状态摄动,Qk和Rk分别为过程噪声协方差和量测噪声协方差,εf为受扰状态和相应线性逼近之间的差值,εh为量测值和相应线性逼近之间的差值,nf和nh分别为状态方程和量测方程的非线性指标;
(5)确定内插因子:采用有限状态机根据步骤(4)的状态方程和量测方程的非线性指标确定内插因子;
(6)根据步骤(5)的内插因子,使用插值法在两个实际量测值之间增加伪量测值;
(7)预测步:采用扩展卡尔曼滤波的预测步计算状态预测值和状态预测误差协方差;
(8)修正步:采用扩展卡尔曼滤波的修正步对步骤(7)的状态预测值进行修正,引入H修正步骤(7)的状态预测误差协方差;
(9)判断滤波次数是否达到时间间隔所插入的伪量测值数量加一,若是,则进入步骤(10),若否,则返回步骤(7);
(10)判断是否达到估计时间长度,若是,则输出结果,退出程序;若否,则返回步骤(4)继续。
步骤(1)中参数信息包括:时间惯性常数、阻尼系数、同步转速、额定功率和发电机总机组数。
步骤(2)中程序初始化包括:设定状态变量初始值、设定系统模型噪声方差矩阵、设定量测误差方差矩阵、设定滤波协方差初始值、设定估计时间长度、设定H参数、设定有限状态机。
有益效果:针对电力系统机电暂态过程的变化特性和发电机本身的非线性,本发明提出一种插值H扩展卡尔曼滤波发电机动态状态估计的算法。首先计算状态方程和量测方程的非线性指标,通过有限状态机确定内插因子,使用插值法在两个实际量测值之间增加伪量测量,然后再进行滤波,以此来减轻非线性在估计精度上的消极影响。通过在扩展卡尔曼滤波的修正步引入H计算估计误差协方差,以此提高了对噪声的鲁棒性能。经过IEEE标准算例测试结果表明,本发明提出方法的滤波精度和对噪声的鲁棒性能均优于扩展卡尔曼滤波器。
附图说明
图1为本发明方法流程图;
图2为以4个状态有限状态机为例的示意图;
图3为IEEE9为标准测试系统图;
图4为实施例采用本发明方法与EKF滤波结果对比图;其中,(a)为发电机1功角动态估计曲线,(b)为发电机1电角速度动态估计曲线;
图5为实施例噪声变化下本方法与EKF滤波结果对比图;其中,(a)为发电机1功角动态估计曲线,(b)为发电机1电角速度动态估计曲线,(c)为发电机2功角动态估计曲线,(d)为发电机2电角速度动态估计曲线,(e)为发电机 3功角动态估计曲线,(f)为发电机3电角速度动态估计曲线。
具体实施方式
下面结合具体实施例,进一步阐明本发明,应理解这些实施例仅用于说明本发明而不用于限制本发明的范围,在阅读了本发明之后,本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。
1.多步自适应插值扩展卡尔曼滤波
在扩展卡尔曼滤波中,通过使用一阶泰勒展开式逼近非线性函数,使得在高度非线性系统中,由于忽略高阶项引起的线性化误差会导致扩展卡尔曼滤波出现错误估计。为了克服这一缺点,本发明采用多步自适应插值法:首先计算状态方程和量测方程的非线性指标nf和nh,再根据非线性指标使用有限状态机确定内插因子,然后通过插值法在两个实际量测值之间增加伪量测值。在增加伪量测值之后,按照H扩展卡尔曼滤波进行状态估计。
非线性指标:
Figure GDA0002427083220000041
Figure GDA0002427083220000042
Figure GDA0002427083220000043
Figure GDA0002427083220000044
式中,f(·)和h(·)分别为状态转换函数和量测函数,xk为k时刻的状态变量,δx为状态摄动,Qk和Rk分别为过程噪声协方差和量测噪声协方差,εf为受扰状态和相应线性逼近之间的差值,εh为量测值和相应线性逼近之间的差值,nf和nh分别为状态方程和量测方程的非线性指标。
2.H扩展卡尔曼滤波
扩展卡尔曼滤波器在修正步引入H滤波器,修正预测误差协方差。
根据预测步得到的预测误差协方差矩阵Pkk-1结合量测方程的雅可比矩阵Hk和量测噪声协方差矩阵Rk,计算修正后的估计误差协方差Pkk
Figure GDA0002427083220000045
Figure GDA0002427083220000046
Figure GDA0002427083220000047
式中上标-1代表矩阵的逆,上标T代表矩阵的转置,λ是比1大的正常数在电力系统动态状态估计中建议值在1和10之间变化,γ是一个中间变量,max{·}为求矩阵中元素的最大值,eig(·)为求矩阵的特征值,Pkk-1为k时刻预测误差协方差矩阵,Hk为k时刻量测方程的雅可比矩阵,Rk为k时刻量测噪声协方差矩阵,Pkk为k时刻估计误差协方差矩阵,I为单位矩阵,inv(·)为求矩阵的逆。
多步自适应插值扩展卡尔曼滤波根据状态方程和量测方程的非线性指标和有限状态机得到插值因子,进而在两实际量测值之间通过插值法增加伪量测值进行扩展卡尔曼滤波,减轻了非线性的消极影响;而H扩展卡尔曼滤波通过在修正步引入H修正预测误差协方差矩阵,提高扩展卡尔曼滤波的鲁棒性能。本发明结合上述两者,对机电暂态过程中的发电机进行动态状态估计,具体步骤如下:
(1)首先获得所需估计发电机机组的参数信息。包括:时间惯性常数、阻尼系数、同步转速、额定功率和发电机总机组数等;
(2)程序初始化。包括:设定状态变量初始值、设定系统模型噪声方差矩阵、设定量测误差方差矩阵、设定滤波协方差初始值、设定估计时间长度、设定 H参数、设定有限状态机;
(3)利用H扩展卡尔曼滤波进行一次滤波;
(4)引入非线性指标:根据下式计算状态方程和量测方程的非线性指标:
Figure GDA0002427083220000051
Figure GDA0002427083220000052
Figure GDA0002427083220000053
Figure GDA0002427083220000054
式中,f(·)和h(·)分别为状态转换函数和量测函数,δx为状态摄动,Qk和Rk分别为过程噪声协方差和量测噪声协方差,εf为受扰状态和相应线性逼近之间的差值,εh为量测值和相应线性逼近之间的差值,nf和nh分别为状态方程和量测方程的非线性指标;
(5)确定内插因子:采用有限状态机根据步骤(4)的状态方程和量测方程的非线性指标确定内插因子;
(6)根据步骤(5)的内插因子,使用插值法在两个实际量测值之间增加伪量测值;
(7)预测步:采用扩展卡尔曼滤波的预测步计算状态预测值和状态预测误差协方差;
(8)修正步:采用扩展卡尔曼滤波的修正步对步骤(7)的状态预测值进行修正,引入H修正步骤(7)的状态预测误差协方差;
(9)判断滤波次数是否达到时间间隔所插入的伪量测值数量加一,若是,则进入步骤(10),若否,则返回步骤(7);
(10)判断是否达到估计时间长度,若是,则输出结果,退出程序;若否,则返回步骤(4)继续。
下面介绍本发明的一个算例:
本发明测试的算例为IEEE9节点标准系统。IEEE9节点量测数据由BPA仿真真值添加随机噪声得到,仿真时采用发电机经典二阶模型并将调速器的作用考虑在内,并假设在第40周波(1周波为0.02s,即电力系统运行周期)时,其中节点4-节点8支路首端发生三相金属性短路,第58周波时短路故障消失。
本发明是基于扩展卡尔曼滤波(extended Kalman filter,EKF)之上进行改进,所以选取扩展卡尔曼滤波与本发明的插值H扩展卡尔曼滤波(interpolation H extendedKalman filter,IHEKF)的动态状态估计算法进行性能比较。
为了让每个算法之间的估计结果比较更加明显,本文采用均方误差(meansquared error,MSE)作为指标进行算法性能间的对比,定义如下:
Figure GDA0002427083220000061
式中,
Figure GDA0002427083220000062
代表k时刻状态变量的滤波值,xk代表k时刻状态变量的真实值(BPA 数据),n为采样周期数。
图4为IEEE9节点系统中发电机1在EKF和本发明提出的IHEKF作用下的滤波曲线。可见在故障出现前,两种算法的滤波精度基本一致。但是当出现故障后,IHEKF和EKF都能够较好的跟踪状态变化。结合表1可以看出,IHEKF 的均方误差明显比EKF的要小很多,在跟踪精度上有所提高。
图5为IEEE9节点系统中发电机在噪声均值非零,方差逐步增大时EKF和本发明提出的IHEKF的滤波曲线。量测噪声的设置情况如下:Case1,发电机1 量测噪声均值为零,标准差为3σ0;Case2,发电机2量测噪声均值为E0,标准差为σ0;Case3,发电机3量测噪声均值为3E0,标准差为3σ0。结合表2可以看出,随着噪声的变化,IHEKF对噪声的鲁棒性能明显比EKF的要强。
表1不同算法下发电机估计结果的均方误差
Figure GDA0002427083220000071
表2不同噪声下二种算法的估计结果
Figure GDA0002427083220000072

Claims (5)

1.一种插值H扩展卡尔曼滤波发电机动态状态估计方法,其特征在于,所述方法是依次按以下步骤实现的:
(1)获得所需估计发电机机组的参数信息;
(2)程序初始化;
(3)利用H扩展卡尔曼滤波进行一次滤波;
(4)引入非线性指标;
(5)确定内插因子:采用有限状态机根据步骤(4)的状态方程和量测方程的非线性指标确定内插因子;
(6)根据步骤(5)的内插因子,使用插值法在两个实际量测值之间增加伪量测值;
(7)预测步:采用扩展卡尔曼滤波的预测步计算状态预测值和状态预测误差协方差;
(8)修正步:采用扩展卡尔曼滤波的修正步对步骤(7)的状态预测值进行修正,引入H修正步骤(7)的状态预测误差协方差;
(9)判断滤波次数是否达到时间间隔所插入的伪量测值数量加一,若是,则进入步骤(10),若否,则返回步骤(7);
(10)判断是否达到估计时间长度,若是,则输出结果,退出程序;若否,则返回步骤(4)继续。
2.如权利要求1所述的插值H扩展卡尔曼滤波发电机动态状态估计方法,其特征在于,根据下式计算状态方程和量测方程的非线性指标:
Figure FDA0002420985020000011
Figure FDA0002420985020000012
Figure FDA0002420985020000013
Figure FDA0002420985020000014
式中,f(·)和h(·)分别为状态转换函数和量测函数,δx为状态摄动,Qk和Rk分别为过程噪声协方差和量测噪声协方差,εf为受扰状态和相应线性逼近之间的差值,εh为量测值和相应线性逼近之间的差值,nf和nh分别为状态方程和量测方程的非线性指标。
3.如权利要求1所述的插值H扩展卡尔曼滤波发电机动态状态估计方法,其特征在于,步骤(1)中参数信息包括:时间惯性常数、阻尼系数、同步转速、额定功率和发电机总机组数。
4.如权利要求1所述的插值H扩展卡尔曼滤波发电机动态状态估计方法,其特征在于,步骤(2)中程序初始化包括:设定状态变量初始值、设定系统模型噪声方差矩阵、设定量测误差方差矩阵、设定滤波协方差初始值、设定估计时间长度、设定H参数、设定有限状态机。
5.如权利要求1所述的插值H扩展卡尔曼滤波发电机动态状态估计方法,其特征在于,根据预测步得到的预测误差协方差矩阵Pkk-1结合量测方程的雅可比矩阵Hk和量测噪声协方差矩阵Rk,计算修正后的估计误差协方差Pkk
Figure FDA0002420985020000021
Figure FDA0002420985020000022
Figure FDA0002420985020000023
式中上标-1代表矩阵的逆,上标T代表矩阵的转置,λ是比1大的正常数在电力系统动态状态估计中建议值在1和10之间变化,γ是一个中间变量,max{·}为求矩阵中元素的最大值,eig(·)为求矩阵的特征值,Pk|k-1为k时刻预测误差协方差矩阵,Hk为k时刻量测方程的雅可比矩阵,Rk为k时刻量测噪声协方差矩阵,Pk|k为k时刻估计误差协方差矩阵,I为单位矩阵,inv(·)为求矩阵的逆。
CN201710811192.8A 2017-09-11 2017-09-11 一种插值h∞扩展卡尔曼滤波发电机动态状态估计方法 Active CN107425548B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710811192.8A CN107425548B (zh) 2017-09-11 2017-09-11 一种插值h∞扩展卡尔曼滤波发电机动态状态估计方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710811192.8A CN107425548B (zh) 2017-09-11 2017-09-11 一种插值h∞扩展卡尔曼滤波发电机动态状态估计方法

Publications (2)

Publication Number Publication Date
CN107425548A CN107425548A (zh) 2017-12-01
CN107425548B true CN107425548B (zh) 2020-06-16

Family

ID=60432674

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710811192.8A Active CN107425548B (zh) 2017-09-11 2017-09-11 一种插值h∞扩展卡尔曼滤波发电机动态状态估计方法

Country Status (1)

Country Link
CN (1) CN107425548B (zh)

Families Citing this family (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107807278A (zh) * 2017-12-06 2018-03-16 河海大学 基于h∞扩展卡尔曼滤波的低频振荡信号参数辨识方法
CN108155648B (zh) * 2018-01-09 2021-01-05 河海大学 基于自适应h无穷扩展卡尔曼滤波的状态估计方法
CN108614804B (zh) * 2018-04-25 2022-09-02 中国人民解放军战略支援部队信息工程大学 基于信噪比检验的正则化卡尔曼滤波方法
CN109100649B (zh) * 2018-06-25 2020-10-16 南京南瑞继保电气有限公司 基于相量测量的发电机励磁系统及调速系统参数估计方法
CN109270455B (zh) * 2018-10-24 2021-02-02 郑州轻工业学院 基于弱敏集合卡尔曼滤波的感应电机状态监测方法
CN110957716A (zh) * 2018-12-29 2020-04-03 上海电机学院 一种基于扩展卡尔曼滤波的增量配电网转供能力优化方法
CN110112770A (zh) * 2019-04-17 2019-08-09 河海大学 一种基于自适应h∞容积卡尔曼滤波的发电机动态状态估计方法
CN110032812A (zh) * 2019-04-18 2019-07-19 河海大学 一种基于自适应容积卡尔曼滤波的动态状态估计方法
CN109918862A (zh) * 2019-04-28 2019-06-21 河海大学 一种基于鲁棒无迹h无穷滤波的发电机动态估计方法
CN110133400B (zh) * 2019-05-10 2021-07-09 青岛大学 一种融合递推状态估计的动态电力系统异常检测方法
CN111948534B (zh) * 2020-07-31 2023-05-05 华北电力科学研究院有限责任公司 发电机状态预警方法及系统
CN112083349B (zh) * 2020-08-01 2023-04-07 南通长江电器实业有限公司 一种永磁同步电机定子绕组匝间短路故障诊断方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6285971B1 (en) * 1997-08-22 2001-09-04 Voyan Technology Method for real-time nonlinear system state estimation and control
CN104777426A (zh) * 2015-04-17 2015-07-15 河海大学 一种基于无迹变换强跟踪的发电机动态状态估计方法
CN105403834A (zh) * 2015-12-22 2016-03-16 华北电力大学 一种发电机动态状态估计方法
CN106786561A (zh) * 2017-02-20 2017-05-31 河海大学 一种基于自适应卡尔曼滤波的低频振荡模态参数辨识方法
CN106844952A (zh) * 2017-01-20 2017-06-13 河海大学 基于无迹粒子滤波理论的发电机动态状态估计方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6285971B1 (en) * 1997-08-22 2001-09-04 Voyan Technology Method for real-time nonlinear system state estimation and control
CN104777426A (zh) * 2015-04-17 2015-07-15 河海大学 一种基于无迹变换强跟踪的发电机动态状态估计方法
CN105403834A (zh) * 2015-12-22 2016-03-16 华北电力大学 一种发电机动态状态估计方法
CN106844952A (zh) * 2017-01-20 2017-06-13 河海大学 基于无迹粒子滤波理论的发电机动态状态估计方法
CN106786561A (zh) * 2017-02-20 2017-05-31 河海大学 一种基于自适应卡尔曼滤波的低频振荡模态参数辨识方法

Also Published As

Publication number Publication date
CN107425548A (zh) 2017-12-01

Similar Documents

Publication Publication Date Title
CN107425548B (zh) 一种插值h∞扩展卡尔曼滤波发电机动态状态估计方法
CN107590317B (zh) 一种计及模型参数不确定性的发电机动态估计方法
CN108155648B (zh) 基于自适应h无穷扩展卡尔曼滤波的状态估计方法
Yin et al. Real-time implementation of fault-tolerant control systems with performance optimization
CN107797448B (zh) 采用扰动扩张补偿的电机位置离散重复控制方法
CN107544244B (zh) 基于椭圆吸引律和等效扰动扩张状态补偿的用于电机伺服系统的离散重复控制方法
CN102779238A (zh) 一种基于自适应卡尔曼滤波的无刷直流电机系统辨识方法
JP2013511703A5 (zh)
KR101394556B1 (ko) 영구자석 동기전동기의 회전자 위치센서 고장검출 장치 및 그 방법
US20150168925A1 (en) Model-based predictive controller with steady-state model adaptation
WO2016147722A1 (ja) 推定装置、推定方法およびプログラム
CN109188908B (zh) 基于指数型无切换吸引律的数字控制器设计方法
CN109586645B (zh) 一种永磁同步电机惯量识别方法及设备
Akhlaghi et al. Exploring adaptive interpolation to mitigate non-linear impact on estimating dynamic states
JP6314543B2 (ja) 蓄電池の充電状態推定方法及び装置
JP2014232649A (ja) 電池温度推定装置及び電池温度推定方法
US8541972B2 (en) Method for suppressing speed ripple by using torque compensator based on activation function
CN110907835A (zh) 一种具有噪声免疫特性的电池模型参数辨识与soc估计方法
CN104838322A (zh) 周期性扰动自动抑制设备
CN108776274B (zh) 一种基于自适应滤波的风电变流器故障诊断
CN112511056A (zh) 一种基于相量测量的鲁棒发电机动态状态估计方法
CN109240085B (zh) 非高斯系统动态数据校正与系统控制性能优化方法
Takemura et al. Proposal of novel rolling friction compensation with data-based friction model for ball screw driven stage
CN111340647A (zh) 一种配电网数据平差方法及系统
CN109768752B (zh) 一种基于多用途扰动观测器的永磁同步电机无差拍电流预测控制方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant