CN107792396A - 发动机工作时干扰力矩实时估计方法 - Google Patents
发动机工作时干扰力矩实时估计方法 Download PDFInfo
- Publication number
- CN107792396A CN107792396A CN201710876529.3A CN201710876529A CN107792396A CN 107792396 A CN107792396 A CN 107792396A CN 201710876529 A CN201710876529 A CN 201710876529A CN 107792396 A CN107792396 A CN 107792396A
- Authority
- CN
- China
- Prior art keywords
- disturbance torque
- thruster
- integration
- data
- time
- 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
Classifications
-
- B—PERFORMING OPERATIONS; TRANSPORTING
- B64—AIRCRAFT; AVIATION; COSMONAUTICS
- B64G—COSMONAUTICS; VEHICLES OR EQUIPMENT THEREFOR
- B64G1/00—Cosmonautic vehicles
- B64G1/22—Parts of, or equipment specially adapted for fitting in or to, cosmonautic vehicles
- B64G1/24—Guiding or controlling apparatus, e.g. for attitude control
- B64G1/242—Orbits and trajectories
-
- 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/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
- G06F17/13—Differential equations
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Optimization (AREA)
- Mathematical Analysis (AREA)
- Computational Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Remote Sensing (AREA)
- Theoretical Computer Science (AREA)
- Aviation & Aerospace Engineering (AREA)
- Radar, Positioning & Navigation (AREA)
- Combustion & Propulsion (AREA)
- Operations Research (AREA)
- Algebra (AREA)
- Chemical & Material Sciences (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Control Of Position, Course, Altitude, Or Attitude Of Moving Bodies (AREA)
Abstract
本发明公开了一种发动机工作时干扰力矩实时估计方法,其包括下列步骤:步骤一,角速度与推力器计时数据采集,高频率采集星上角速度,慢传推力器工作计时、喷气量数据;步骤二,迭代时间选取,迭代步长表示两次积分求和的时间间隔,决定了干扰力矩计算时数据输出频率;步骤三,积分求和;步骤四,干扰力矩解算,对卫星姿态动力学方程等式两边各项分别积分,利用上面步骤中求得的各项积分结果,反解得到干扰力矩。本发明提出的方法简单,且可实现多用途复用,可在高轨道卫星转移段变轨过程、深空探测飞行器轨道控制或者其他轨道控制过程使用,也可以在其他航天器中推广使用。
Description
技术领域
本发明涉及一种干扰力矩实时估计方法,特别是涉及一种发动机工作时干扰力矩实时估计方法。
背景技术
高轨道卫星在转移轨道段需利用490N发动机进行轨道控制。由于490N发动机自身结构以及在地面调试校准过程中不可能做到完全精准无误,发动机工作时会产生偏离理论推力方向的推力。如果卫星质心在490N发动机推力矢量方向上,推力不会对姿态控制产生干扰力矩;但是由于490N发动机自身偏斜以及卫星多次变轨过程中燃料消耗引起质心会不断变化,490N发动机推力不可能始终通过卫星质心,必然会产生干扰力矩,引起用于姿态控制的推力器工作多消耗燃料,甚至卫星姿态失稳。因此对490N发动机工作时干扰力矩估计具有重要意义。
另外,卫星在轨运行期间会受到空间干扰力矩的影响,包括外界环境引起的干扰力矩以及自身转动部件产生的干扰力矩的影响。为了减小干扰力矩对姿态的影响,一般在地面对干扰力矩进行建模,然后在系统设计时对干扰力矩进行补偿。然而由于在轨空间环境、490N发动机实际工作等状态与地面建模存在差异,造成干扰力矩模型存在一定的误差和不确定性,这样对姿态控制精度、干扰力矩补偿效果等与理想值有一定误差。因此,因此通过在轨星上实时数据对干扰力矩进行实时估计是十分必要的。
目前国内外都在对卫星在轨干扰力矩的模型辨识进行研究。但是较多集中在卫星在轨稳态运行阶段,对于高轨道卫星、深空探测卫星等大推力发动机在变轨过程干扰力矩估计还没有发现相关的报道。
发明内容
本发明所要解决的技术问题是提供一种发动机工作时干扰力矩实时估计方法,其在不增加卫星数据采集前提下,利用卫星的实时角速度和推力器工作计时等数据,估计490N发动机工作时干扰力矩信息,通过引入耦合项积分以及必要的近似假设,保证估计准确性的同时,使过程简单易行,本发明提出的方法简单,且可实现多用途复用,可在高轨道卫星转移段变轨过程、深空探测飞行器轨道控制或者其他轨道控制过程使用,也可以在其他航天器中推广使用。
本发明是通过下述技术方案来解决上述技术问题的:一种发动机工作时干扰力矩实时估计方法,其包括下列步骤:
步骤一,角速度与推力器计时数据采集,高频率采集星上角速度,慢传推力器工作计时、喷气量数据;
步骤二,迭代时间选取,迭代步长表示两次积分求和的时间间隔,决定了干扰力矩计算时数据输出频率;
步骤三,积分求和;
步骤四,干扰力矩解算,对卫星姿态动力学方程等式两边各项分别积分,利用上面步骤中求得的各项积分结果,反解得到干扰力矩。
优选地,所述步骤三包括下列步骤:
步骤五,求取动量矩的积分;
步骤六,求取推力器工作控制力矩积分;
步骤七,求取干扰力矩积分;
步骤八,求取其它坐标轴与本轴的耦合项积分。
优选地,所述步骤三对于动力学中各项分别根据各自特点进行数值积分,动量矩的积分等于积分时长内角动量的变化,推力器工作控制力矩只存在推力器工作时间,因此积分时长中只有一段时间有作用;干扰力矩一直存在,整个积分时长内都由作用,耦合项积分采用每一时刻的数据乘以步长叠加得到。
优选地,所述步骤七求取干扰力矩积分过程中考虑了耦合项的影响。
优选地,所述步骤一基于卫星姿态动力学原理,利用星上角速度测量数据和推力器工作计时数据估计干扰力矩,角速度数据和推力器计时数据更新快慢不同,并根据数据更新频率确定迭代步长以及积分时长。
优选地,所述步骤四通过对卫星姿态动力学方程中各项的积分来反解干扰力矩,而不是直接用存在角速度差分的原方程反解干扰力矩。
本发明的积极进步效果在于:本发明发动机工作时干扰力矩实时估计方法在不增加卫星数据采集前提下,利用卫星的实时角速度和推力器工作计时等数据,估计490N发动机工作时干扰力矩信息,通过引入耦合项积分以及必要的近似假设,保证估计准确性的同时,使过程简单易行,本发明提出的方法简单,且可实现多用途复用,可在高轨道卫星转移段变轨过程、深空探测飞行器轨道控制或者其他轨道控制过程使用,也可以在其他航天器中推广使用。
附图说明
图1为本发明发动机工作时干扰力矩实时估计方法的流程图。
图2为发动机工作时对卫星本体+X轴的干扰力矩的示意图。
图3为发动机工作时对卫星本体+Y轴的干扰力矩的示意图。
图4发动机工作时对卫星本体+Z轴的干扰力矩的示意图。
具体实施方式
下面结合附图给出本发明较佳实施例,以详细说明本发明的技术方案。
如图1所示,本发明发动机工作时干扰力矩实时估计方法包括下列步骤:
步骤一,角速度与推力器计时数据采集,高频率采集星上角速度,慢传推力器工作计时、喷气量数据;
步骤二,迭代时间选取,迭代步长表示两次积分求和的时间间隔,决定了干扰力矩计算时数据输出频率;
步骤三,积分求和;
步骤四,干扰力矩解算,对卫星姿态动力学方程等式两边各项分别积分,利用上面步骤中求得的各项积分结果,反解得到干扰力矩。
步骤三包括下列步骤:
步骤五,求取动量矩的积分;
步骤六,求取推力器工作控制力矩积分;
步骤七,求取干扰力矩积分;
步骤八,求取其它坐标轴与本轴的耦合项积分。
步骤三对于动力学中各项分别根据各自特点进行数值积分,动量矩的积分等于积分时长内角动量的变化,推力器工作控制力矩只存在推力器工作时间,因此积分时长中只有一段时间有作用;干扰力矩一直存在,整个积分时长内都由作用,耦合项积分采用每一时刻的数据乘以步长叠加得到。
步骤七求取干扰力矩积分过程中考虑了耦合项的影响。
步骤一基于卫星姿态动力学原理,利用星上角速度测量数据和推力器工作计时数据估计干扰力矩,角速度数据和推力器计时数据更新快慢不同,并根据数据更新频率确定迭代步长以及积分时长。
步骤四通过对卫星姿态动力学方程中各项的积分来反解干扰力矩,而不是直接用存在角速度差分的原方程反解干扰力矩。
本发明发动机工作时干扰力矩实时估计方法基本原理为卫星姿态动力学方程(其中,H为卫星角动量,ω为卫星星体角速度,Mc为卫星控制力矩,Md为卫星受到的干扰力矩),根据星上角速度、姿态控制推力器工作计时数据实时解算干扰力矩。
姿态动力学公式可以进一步表示为对上式积分反解得到干扰力矩。采用这样积分求和方法的好处是对数据进行了滤波处理,尤其是角速度数据。否则对角速度求差分,会引入了很大噪声。具体过程包括以下环节:
步骤一,角速度与推力器计时数据采集。
星上角速度反应卫星姿态变化,一般实时性要求较高,采集频率较高;推力器工作计时、喷气量可以慢传。星上角速度测量数据更新周期为T1;推力器计时更新周期为T2(T2是T1的整数倍)。这样推力器的推力器计时数据在T2时间内被采集m=T2/T1次。
步骤二,迭代时间选取。
迭代步长T3表示两次积分求和的时间间隔,决定了干扰力矩计算时数据输出频率,不小于角速度测量数据更新周期T1,取为nT1,(1≤n≤m,n为常整数)。由于推力器计时更新周期为T2,为保证数据有效,积分时长T4不小于推力器计时更新周期T2。
步骤三,积分求和。步骤三包括下列步骤:
步骤五,动量矩的积分,即对姿态动力学方程中动量矩项积分。动量矩的积分可以用星体角动量变化来表示,积分结束时刻角动量减去积分初始时刻角动量就是动量矩的积分,如下式(1):
上式中,分别表示三轴角速度的微分,即三轴角加速度。为t0+T4时刻x轴角速度,为t0时刻x轴角速度,ωx(i+m)为采集数据中第i+m个x轴角速度,ωxi为采集数据中第i个x轴角速度,另外两个方向类同。
步骤六,推力器工作控制力矩Mc积分,仅在推力器工作时有控制力矩作用,因此需要分别求出产生某一轴正负力矩的推力器工作时长。推力器控制力矩为常值,实际可以计算得到。对各个推力产生的控制力矩与各自工作时长的乘积进行叠加,以实现推力器工作控制力矩积分,如下式(2):
上式中,Mcx为推力器产生的x轴控制力矩,McxP为推力器在x轴产生的正向力矩,McxN推力器在x轴产生的负向力矩,TxP(i+m)、TxPi分别为产生x轴正向力矩的推力器第(i+m)个和第i个喷气计时数据,TxN(i+m)、TxNi分别为产生x轴负向力矩的推力器第(i+m)个和第i个喷气计时数据。另外两轴类同。
步骤七,干扰力矩Md积分,干扰力矩在整个过程中一直存在,而且变化缓慢,认为在积分时长内是不变化的,因此干扰力矩积分用干扰力矩与积分时长的乘积得到,如下式(3):
上式中,Mdx,Mdy,Mdz上式分别表示490N发动机工作时的卫星三轴干扰力矩。
步骤八,其它坐标轴与本轴的耦合项积分。
对X方向来说,耦合项为-(Jyy-Jzz)ωyωz,对其积分如下式(4):
对Y方向来说,耦合项为-(Jzz-Jxx)ωzωx,对其积分如下式(5):
对Z方向来说,耦合项为-(Jxx-Jyy)ωxωy,对其积分如下式(6):
步骤四,干扰力矩解算。
对卫星姿态动力学方程等式两边各项分别积分,分解可得如下列公式(7)和(8)所示:
利用上面步骤中求得的各项积分结果,反解得到三轴干扰力矩。
根据数据需求确定得到的迭代步长T3,每隔T3时长取一个点作为本次积分的t0时刻,然后重复上面步骤三与步骤四,得到连续干扰力矩输出。
下面给出采用上述方法对某型号高轨卫星变轨490N发动机工作过程中干扰力矩估计的实例。实际参数如下:某次变轨末期,该卫星惯量Jxx、Jyy、Jzz分别为6719.04、6086.74、7512.04,单位为kg.m2。姿态控制推力器产生的x轴控制力矩McxP=McxN=18.9Nm,产生的y轴控制力矩McyP=McyN=19.75Nm,产生的z轴控制力矩MczP=MczN=15.4Nm。角速度更新周期T1为0.5s,推力器工作计时更新周期T2为16s。为保证数据量适中和数据有效,迭代步长T3和积分时长T4均取16s。依据上述方法估计得到某段时间490N工作时对卫星本体+X轴、+Y轴、+Z轴产生的干扰力矩散点图如图2至图4所示。
本发明由于采取上述方法,在不增加卫星数据采集前提下,利用卫星的实时角速度和推力器工作计时等数据,估计490N发动机工作时干扰力矩信息,通过引入耦合项积分以及必要的近似假设,保证估计准确性的同时,使过程简单易行,本发明提出的方法简单,且可实现多用途复用,可在高轨道卫星转移段变轨过程、深空探测飞行器轨道控制或者其他轨道控制过程使用,也可以在其他航天器中推广使用。本发明的结果可用于高轨道卫星、深空探测飞行器大推力发动机变轨过程中姿态动力学模型修正,并用于姿态控制律的设计。本发明主要用于490N发动机。本发明利用卫星在轨实时角速度、姿态控制推力器的工作计时数据,可以快速地估计出卫星受到的干扰力矩,进而实现精确建模与控制。
以上所述的具体实施例,对本发明的解决的技术问题、技术方案和有益效果进行了进一步详细说明,所应理解的是,以上所述仅为本发明的具体实施例而已,并不用于限制本发明,凡在本发明的精神和原则之内,所做的任何修改、等同替换、改进等,均应包含在本发明的保护范围之内。
Claims (6)
1.一种发动机工作时干扰力矩实时估计方法,其特征在于,其包括下列步骤:
步骤一,角速度与推力器计时数据采集,高频率采集星上角速度,慢传推力器工作计时、喷气量数据;
步骤二,迭代时间选取,迭代步长表示两次积分求和的时间间隔,决定了干扰力矩计算时数据输出频率;
步骤三,积分求和;
步骤四,干扰力矩解算,对卫星姿态动力学方程等式两边各项分别积分,利用上面步骤中求得的各项积分结果,反解得到干扰力矩。
2.如权利要求1所述的发动机工作时干扰力矩实时估计方法,其特征在于,所述步骤三包括下列步骤:
步骤五,求取动量矩的积分;
步骤六,求取推力器工作控制力矩积分;
步骤七,求取干扰力矩积分;
步骤八,求取其它坐标轴与本轴的耦合项积分。
3.如权利要求2所述的发动机工作时干扰力矩实时估计方法,其特征在于,所述步骤三对于动力学中各项分别根据各自特点进行数值积分,动量矩的积分等于积分时长内角动量的变化,推力器工作控制力矩只存在推力器工作时间,因此积分时长中只有一段时间有作用;干扰力矩一直存在,整个积分时长内都由作用,耦合项积分采用每一时刻的数据乘以步长叠加得到。
4.如权利要求2所述的发动机工作时干扰力矩实时估计方法,其特征在于,所述步骤七求取干扰力矩积分过程中考虑了耦合项的影响。
5.如权利要求1所述的发动机工作时干扰力矩实时估计方法,其特征在于,所述步骤一基于卫星姿态动力学原理,利用星上角速度测量数据和推力器工作计时数据估计干扰力矩,角速度数据和推力器计时数据更新快慢不同,并根据数据更新频率确定迭代步长以及积分时长。
6.如权利要求1所述的发动机工作时干扰力矩实时估计方法,其特征在于,所述步骤四通过对卫星姿态动力学方程中各项的积分来反解干扰力矩,而不是直接用存在角速度差分的原方程反解干扰力矩。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201710876529.3A CN107792396B (zh) | 2017-09-25 | 2017-09-25 | 发动机工作时干扰力矩实时估计方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201710876529.3A CN107792396B (zh) | 2017-09-25 | 2017-09-25 | 发动机工作时干扰力矩实时估计方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN107792396A true CN107792396A (zh) | 2018-03-13 |
CN107792396B CN107792396B (zh) | 2021-02-19 |
Family
ID=61531985
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201710876529.3A Active CN107792396B (zh) | 2017-09-25 | 2017-09-25 | 发动机工作时干扰力矩实时估计方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN107792396B (zh) |
Cited By (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN108820257A (zh) * | 2018-06-25 | 2018-11-16 | 上海卫星工程研究所 | 一种卫星用步进电机机构的力矩校核方法 |
CN109018442A (zh) * | 2018-06-15 | 2018-12-18 | 上海卫星工程研究所 | 新型低成本卫星三轴姿态分时解耦高复用喷气控制方法 |
CN110362112A (zh) * | 2019-07-22 | 2019-10-22 | 江南机电设计研究所 | 一种抑制发动机干扰的引入方法 |
CN111026142A (zh) * | 2019-12-11 | 2020-04-17 | 北京控制工程研究所 | 一种大干扰和小惯量情况下的快速姿态机动方法及系统 |
CN117184456A (zh) * | 2023-11-08 | 2023-12-08 | 北京控制工程研究所 | 轨控发动机干扰力矩的估计方法、装置、设备及介质 |
Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20080251646A1 (en) * | 2007-04-13 | 2008-10-16 | The Boeing Company | Singularity Escape and Avoidance Using A Virtual Array Rotation |
CN103303495A (zh) * | 2013-04-11 | 2013-09-18 | 北京控制工程研究所 | 一种动力下降过程干扰力矩的估计方法 |
CN105116906A (zh) * | 2015-07-17 | 2015-12-02 | 中国空间技术研究院 | 基于向量理论的航天器变轨发动机干扰力矩计算方法 |
-
2017
- 2017-09-25 CN CN201710876529.3A patent/CN107792396B/zh active Active
Patent Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20080251646A1 (en) * | 2007-04-13 | 2008-10-16 | The Boeing Company | Singularity Escape and Avoidance Using A Virtual Array Rotation |
CN103303495A (zh) * | 2013-04-11 | 2013-09-18 | 北京控制工程研究所 | 一种动力下降过程干扰力矩的估计方法 |
CN105116906A (zh) * | 2015-07-17 | 2015-12-02 | 中国空间技术研究院 | 基于向量理论的航天器变轨发动机干扰力矩计算方法 |
Non-Patent Citations (1)
Title |
---|
梁彤, 张奕群等: ""一种空间飞行器轨控发动机干扰力矩的测试方法"", 《现代防御技术》 * |
Cited By (8)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN109018442A (zh) * | 2018-06-15 | 2018-12-18 | 上海卫星工程研究所 | 新型低成本卫星三轴姿态分时解耦高复用喷气控制方法 |
CN108820257A (zh) * | 2018-06-25 | 2018-11-16 | 上海卫星工程研究所 | 一种卫星用步进电机机构的力矩校核方法 |
CN108820257B (zh) * | 2018-06-25 | 2020-05-29 | 上海卫星工程研究所 | 一种卫星用步进电机机构的力矩校核方法 |
CN110362112A (zh) * | 2019-07-22 | 2019-10-22 | 江南机电设计研究所 | 一种抑制发动机干扰的引入方法 |
CN111026142A (zh) * | 2019-12-11 | 2020-04-17 | 北京控制工程研究所 | 一种大干扰和小惯量情况下的快速姿态机动方法及系统 |
CN111026142B (zh) * | 2019-12-11 | 2023-04-14 | 北京控制工程研究所 | 一种大干扰和小惯量情况下的快速姿态机动方法及系统 |
CN117184456A (zh) * | 2023-11-08 | 2023-12-08 | 北京控制工程研究所 | 轨控发动机干扰力矩的估计方法、装置、设备及介质 |
CN117184456B (zh) * | 2023-11-08 | 2024-01-30 | 北京控制工程研究所 | 轨控发动机干扰力矩的估计方法、装置、设备及介质 |
Also Published As
Publication number | Publication date |
---|---|
CN107792396B (zh) | 2021-02-19 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN107792396A (zh) | 发动机工作时干扰力矩实时估计方法 | |
CN102175246B (zh) | 一种x脉冲星探测器等效器的航天器导航系统 | |
Gair et al. | Forced motion near black holes | |
Cho et al. | Explicit solution to the full nonlinear problem for satellite formation-keeping | |
CN103674032B (zh) | 融合脉冲星辐射矢量和计时观测的卫星自主导航系统及方法 | |
CN103112603B (zh) | 欠驱动高速自旋卫星建立正常姿态的方法 | |
CN105843237A (zh) | 一种用于抑制柔性振动的航天器姿态参考指令生成方法 | |
CN109885075A (zh) | 一种基于t-s模糊建模的卫星姿态抗干扰容错控制方法 | |
CN102073280A (zh) | 一种复杂挠性航天器模糊奇异摄动建模与姿态控制方法 | |
CN104142686A (zh) | 一种卫星自主编队飞行控制方法 | |
CN103970964A (zh) | 一种挠性卫星模态参数在轨辨识方法 | |
CN104848862B (zh) | 一种环火探测器精密同步定位守时方法及系统 | |
CN104656447A (zh) | 一种航天器抗干扰姿态跟踪的微分几何非线性控制方法 | |
CN107402516B (zh) | 基于联合执行机构的递阶饱和模糊pd姿态控制方法 | |
CN105912007A (zh) | 空间机械臂抗干扰姿态稳定的微分几何非线性控制方法 | |
CN105930305B (zh) | 一种三脉冲交会接近制导方法 | |
CN105759827A (zh) | 一种抑制不期望柔性振动的航天器姿态控制系统 | |
CN104316048A (zh) | 一种普适性的脉冲星自主导航测量模型构建方法 | |
CN106248300A (zh) | 基于成对推力器连续工作的卫星质心位置测量方法 | |
CN110134137A (zh) | 基于扩张状态观测器的航天器姿态跟踪控制方法 | |
CN104192322A (zh) | 一种行星动力下降段轨迹在线生成的抗干扰制导控制方法 | |
Babcock | CubeSat attitude determination via Kalman filtering of magnetometer and solar cell data | |
CN107421533A (zh) | 一种深空探测器x射线脉冲星toa/dtoa组合导航方法 | |
CN103678897B (zh) | 一种基于凯恩方程的飞轮隔振平台专用动力学建模方法 | |
CN110955255B (zh) | 基于cmg的高精度轨控姿态维持方法、系统及介质 |
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 |