CN110489791B - 一种脉冲速度增量改变轨道根数的高阶解析计算方法 - Google Patents

一种脉冲速度增量改变轨道根数的高阶解析计算方法 Download PDF

Info

Publication number
CN110489791B
CN110489791B CN201910621394.5A CN201910621394A CN110489791B CN 110489791 B CN110489791 B CN 110489791B CN 201910621394 A CN201910621394 A CN 201910621394A CN 110489791 B CN110489791 B CN 110489791B
Authority
CN
China
Prior art keywords
order
pulse
expression
track
speed increment
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
CN201910621394.5A
Other languages
English (en)
Other versions
CN110489791A (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.)
Northwestern Polytechnical University
Original Assignee
Northwestern Polytechnical University
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 Northwestern Polytechnical University filed Critical Northwestern Polytechnical University
Priority to CN201910621394.5A priority Critical patent/CN110489791B/zh
Publication of CN110489791A publication Critical patent/CN110489791A/zh
Application granted granted Critical
Publication of CN110489791B publication Critical patent/CN110489791B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10Complex mathematical operations
    • G06F17/11Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G06F17/13Differential equations

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • Operations Research (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Algebra (AREA)
  • Complex Calculations (AREA)

Abstract

本发明提供一种脉冲速度增量改变轨道根数的高阶解析计算方法,包括以下步骤:1,输入脉冲速度增量改变轨道根数的计算方法所需要达到的精度阶数N;2,采用微分形式的高斯变分方程通过链式法则依次得到轨道根数的1~N阶导数;3,采用泰勒展开法得到轨道根数变化关于脉冲速度增量的N阶解析函数;4,利用所建立的N阶解析函数计算给定脉冲速度增量下的轨道根数改变量,完成高阶解析计算。本发明的计算方法是完全解析的,因而计算效率高;且能够达到任意阶次,因而计算精度高。这在航天器轨道控制方面具有广泛的使用价值。

Description

一种脉冲速度增量改变轨道根数的高阶解析计算方法
技术领域
本发明涉及航空航天技术领域,尤其涉及一种脉冲速度增量改变轨道根数的高阶解析计算方法。
背景技术
航天器的轨道控制通常是通过推力器瞬间产生的等效脉冲速度增量实现的。在任意给定的脉冲速度增量下,航天器的轨道根数变化量可通过高斯变分方程来评估。现在被广泛使用的积分型高斯变分方程(IGVEs)是在微分型高斯变分方程(DGVEs)基础上通过一阶近似和脉冲速度增量等效替换得到的,因此也称作一阶IGVEs。值得指出的是,尽管DGVEs是精确的,但它只能适用于连续力控制模式,而不能适用于脉冲速度增量控制模式,因而在工程上使用较少。但能够适用于脉冲速度增量控制模式的一阶IGVEs的精度不高,只能用于评估短时间或较小脉冲速度增量控制下的轨道改变,这大大限制了工程实践的应用。为此,有学者提出了构造二阶IGVEs方程的方法。该方法的实质是对轨道积分与位置速度矢量之间的解析函数进行二阶泰勒级数展开。相比一阶 IGVEs,二阶IGVEs能够以更高精度评估脉冲速度增量产生的轨道根数变化。但用于求出该二阶IGVEs的方法具有明显的缺点,即轨道积分与位置速度矢量之间的解析函数是非线性的、隐形式的,非常不便于进行高阶解的进一步求解。这也是上述方法仅得到了二阶IGVEs而没有得到三阶及更高阶IGVEs的原因。考虑到航天器轨道精确控制对更高阶IGVEs的需求,有必要采用一种简易便捷且通用的方法,从而构造出解析的任意阶IGVEs。
发明内容
为解决现有技术中存在的问题,本发明的目的在于提供一种脉冲速度增量改变轨道根数的高阶解析计算方法,本发明能够快速精确地确定给定脉冲速度增量下航天器轨道根数的变化。
为实现上述目的,本发明采用以下技术手段:
一种脉冲速度增量改变轨道根数的高阶解析计算方法,包括以下步骤:
S1,获取用于计算脉冲速度增量改变轨道根数的解析解的阶数N;其中N 为任意正整数;
S2,根据微分形式的高斯变分方程通过链式法则依次得到轨道根数的1~N 阶导数;
S3,根据泰勒展开法得到轨道根数变化关于脉冲速度增量的N阶解析函数;
S4,根据建立的解析函数计算给定脉冲速度增量下的轨道根数改变量,完成高阶解析计算。
作为本发明的进一步改进,步骤S2具体为:
S2.1,给出轨道根数的微分形式的高斯变分方程;
S2.2,采用微积分的链式法则,对上述微分形式的高斯变分方程进一步求导。
作为本发明的进一步改进,步骤S2具体为:
S2.1,给出微分形式的高斯变分方程如下:
Figure BDA0002125636080000021
Figure BDA0002125636080000022
Figure BDA0002125636080000031
Figure BDA0002125636080000032
Figure BDA0002125636080000033
Figure BDA0002125636080000034
其中,等式左边的t为时间,a、e、i、Ω、ω、f为轨道根数,p=a(1-e2)、
Figure BDA0002125636080000035
为与轨道根数相关的两个轨道参数,μ为地球引力系数,Fr、Ft、Fn为连续控制力;具体地:
a、轨道半长轴;e、轨道偏心率;i、轨道倾角;Ω、升交点赤经;ω、纬度幅角;f、真近点角;p、轨道通径;h、轨道角动量;Fr、径向控制力;Ft、切向控制力;Fn、法向控制力;
S2.2,采用微积分的链式法则,对上述微分形式的高斯变分方程进一步求导:对任意六个轨道根数中的一个αk,若已得到其1~m-1阶的导数,则按照下述公式求解其第m阶的导数:
Figure BDA0002125636080000036
其中,
Figure BDA0002125636080000037
的表达式采用微分型高斯变分方程;
当m=N时,终止上述求导过程,否则继续。
作为本发明的进一步改进,步骤S3具体为:
S3.1,采用泰勒展开法得到轨道根数变化关于脉冲速度增量的N阶形式表达式;
S3.2,将表达式中的连续力和时间脉冲按照下述方式结合形成脉冲速度增量;
S3.3,将N阶形式表达式的连续力表达式变换为脉冲速度增量表达式,完成了N阶IGVEs的构建。
作为本发明的进一步改进,步骤S3具体为:
S3.1,采用泰勒展开法得到轨道根数变化关于脉冲速度增量的N阶形式表达式:
Figure BDA0002125636080000041
S3.2,将表达式中的连续力和时间脉冲按照下述方式结合形成脉冲速度增量:
Δvr=Fr·Δt,Δvt=Ft·Δt,Δvn=Fn·Δt (4)
S3.3,考虑公式(4)中的等价关系,将公式(3)中的连续力表达式变换为脉冲速度增量表达式,从而得到:
Figure BDA0002125636080000042
其中,ηkpq是经整理后得到的系数;完成了N阶IGVEs的构建。
与现有技术相比,本发明的有益效果是:
本发明先输入脉冲速度增量改变轨道根数的计算方法所需要达到的精度阶数,再采用微分形式的高斯变分方程通过链式法则依次得到1~N阶轨道根数导数;然后采用泰勒展开法得到轨道根数变化关于脉冲速度增量的N阶解析函数;最后利用所建立的解析函数计算给定脉冲速度增量下的轨道根数改变量,完成高阶解析计算。该方法能够逐次计算得到任意高阶的积分型高斯变分方程,因而可根据工程需要提高问题求解的精度。
优选地,所得结果具有解析表达式,因而计算代价小、计算效率高,适合在轨实时应用。
附图说明
图1是本发明的高阶解析计算方法的计算流程示意图;
图2是1阶、2阶、3阶IGVEs方法在10000次蒙特卡洛仿真中轨道半长轴改变量的计算误差统计结果。
图3是1阶、2阶、3阶IGVEs方法在10000次蒙特卡洛仿真中轨道偏心率改变量的计算误差统计结果。
图4是1阶、2阶、3阶IGVEs方法在10000次蒙特卡洛仿真中轨道倾角改变量的计算误差统计结果。
图5是1阶、2阶、3阶IGVEs方法在10000次蒙特卡洛仿真中升交点赤经改变量的计算误差统计结果。
图6是1阶、2阶、3阶IGVEs方法在10000次蒙特卡洛仿真中纬度幅角改变量的计算误差统计结果。
图7是1阶、2阶、3阶IGVEs方法在10000次蒙特卡洛仿真中真近点角改变量的计算误差统计结果。
具体实施方式
下面结合附图对本发明的结构和工作原理作进一步详细说明。
如图1所示,本发明提出的一种脉冲速度增量改变轨道根数的高阶解析计算方法,包括以下步骤:
S1,输入脉冲速度增量改变轨道根数的计算方法所需要达到的精度阶数N;其中N可以为1、2、3等任意正整数;
S2,采用微分形式的高斯变分方程通过链式法则依次得到轨道根数的1~N 阶导数;
S3,采用泰勒展开法得到轨道根数变化关于脉冲速度增量的N阶解析函数;
S4,利用所建立的解析函数计算给定脉冲速度增量下的轨道根数改变量,完成高阶解析计算。
优选的,S2具体包括以下步骤:
S2.1,给出微分形式的高斯变分方程如下:
Figure BDA0002125636080000061
Figure BDA0002125636080000062
Figure BDA0002125636080000063
Figure BDA0002125636080000064
Figure BDA0002125636080000065
Figure BDA0002125636080000066
其中,等式左边的t为时间,a、e、i、Ω、ω、f为轨道根数,p=a(1-e2)、
Figure BDA0002125636080000067
为与轨道根数相关的两个轨道参数,μ为地球引力系数,Fr、Ft、Fn为连续控制力。具体地:
a——轨道半长轴;
e——轨道偏心率;
i——轨道倾角;
Ω——升交点赤经;
ω——纬度幅角;
f——真近点角;
p——轨道通径;
h——轨道角动量;
Fr——径向控制力;
Ft——切向控制力;
Fn——法向控制力;
S2.2,采用微积分的链式法则,对上述微分形式的高斯变分方程进一步求导。即,对任意六个轨道根数中的一个αk,若已得到其1~m-1阶的导数,则按照下述公式求解其第m阶的导数:
Figure BDA0002125636080000071
其中,
Figure RE-GDA0002209124770000072
的表达式采用微分型高斯变分方程。
当m=N时,终止上述求导过程,否则继续。
优选的,S3具体包括以下步骤:
S3.1,采用泰勒展开法得到轨道根数变化关于脉冲速度增量的N阶形式表达式:
Figure BDA0002125636080000081
S3.2,将表达式中的连续力和时间脉冲按照下述方式结合形成脉冲速度增量:
Δvr=Fr·Δt,Δvt=Ft·Δt,Δvn=Fn·Δt (4)
S3.3,考虑公式(4)中的等价关系,将公式(3)中的连续力表达式变换为脉冲速度增量表达式,从而得到:
Figure BDA0002125636080000042
其中,ηkpq是经整理后得到的系数。这样能就完成了N阶IGVEs的构建。
下面列举一个具体实施例和附图,说明本发明的具体计算过程。
实施例1
如图1所示,本发明的计算方法包括:
S1,输入脉冲速度增量改变轨道根数的计算方法所需要达到的精度阶数 N=3;
S2,采用微分形式的高斯变分方程通过链式法则依次得到轨道根数的1、2、 3阶导数;具体地,
S2.1给出微分形式的高斯变分方程作为轨道根数的1阶导数:
Figure BDA0002125636080000083
Figure BDA0002125636080000084
Figure BDA0002125636080000091
Figure BDA0002125636080000092
Figure BDA0002125636080000093
Figure BDA0002125636080000094
S2.2对上述方程进行求导,采用链式法则,得到轨道根数的2阶导数:
Figure BDA0002125636080000095
Figure BDA0002125636080000096
Figure BDA0002125636080000097
Figure BDA0002125636080000098
Figure BDA0002125636080000101
Figure BDA0002125636080000102
S2.3对上述方程进行求导,利用链式法则,得到轨道根数的3阶导数:
Figure BDA0002125636080000103
Figure BDA0002125636080000111
Figure BDA0002125636080000112
Figure BDA0002125636080000113
Figure BDA0002125636080000121
Figure BDA0002125636080000122
S3,采用泰勒展开法得到轨道根数变化关于脉冲速度增量的3阶解析函数;具体地,
S3.1,采用泰勒展开法得到轨道根数变化关于脉冲速度增量的3阶形式表达式:
Figure BDA0002125636080000131
S3.2,将表达式中的连续力和时间脉冲按照下述方式结合形成脉冲速度增量:
Δvr=Fr·Δt,Δvt=Ft·Δt,Δvn=Fn·Δt (5)
S3.3,考虑公式(5)中的等价关系,将公式(4)中的连续力表达式变换为脉冲速度增量表达式,从而得到:
Figure BDA0002125636080000132
其中上标(1)、(2)、(3)分别表示1阶、2阶、3阶项,具体表达式为:
Figure BDA0002125636080000133
Figure BDA0002125636080000134
Figure BDA0002125636080000135
Figure BDA0002125636080000136
Figure BDA0002125636080000137
Figure BDA0002125636080000138
Figure BDA0002125636080000141
Figure BDA0002125636080000142
Figure BDA0002125636080000143
Figure BDA0002125636080000144
Figure BDA0002125636080000145
Figure BDA0002125636080000151
Figure BDA0002125636080000152
Figure BDA0002125636080000153
Figure BDA0002125636080000154
Figure BDA0002125636080000161
Figure BDA0002125636080000162
Figure BDA0002125636080000171
S4,利用所建立的解析函数计算给定脉冲速度增量下的轨道根数改变量,完成高阶解析计算。
结果验证,采用10000次蒙特卡洛仿真说明本发明方法的有效性。
假设参考轨道的六个轨道根数按照如下方式随机选取:
a∈[amin,amax],其中amin=RE+300km,amax=RE+20000km,RE=6378.14km;
e∈[emin,emax],其中emin=0.01,emax=1-amin/a;
i∈[0,π],Ω∈[0,2π],ω∈[0,2π],f∈[0,2π].
假设脉冲速度增量按照如下方式随机选取:
Δvr∈[-100,100]m/s、Δvt∈[-100,100]m/s、Δvn∈[-100,100]m/s。
则通过采用精确的非线性的位置速度到轨道根数的函数关系计算得到轨道根数的真实变化量,记作Δαtrue。按照本发明提供的1阶、2阶、3阶解,得到的相应轨道根数变化量记作Δαest。计算两者之间的偏差,并按照偏差量级进行统计,得到结果如图2~7所示。由图可知,对六个轨道根数中的任意一个,三阶解的精度总体上高于二阶解,而二阶解的精度总体上高于一阶解。
以上实施例仅用于说明本发明的技术方案而非对其限制,尽管参照上述实施例对本发明进行了详细说明,所属领域的普通技术人员依然可以对本发明的具体实施方案进行修改或者等同替换,而这些并未脱离本发明精神和范围的任何修改或者等同替换,其均在本发明的权利要求保护范围之内。

Claims (3)

1.一种脉冲速度增量改变轨道根数的高阶解析计算方法,其特征在于,包括以下步骤:
S1,获取用于计算脉冲速度增量改变轨道根数的解析解的阶数N;其中N为任意正整数;
S2,根据微分形式的高斯变分方程通过链式法则依次得到轨道根数的1~N阶导数;
S3,根据泰勒展开法得到轨道根数变化关于脉冲速度增量的N阶解析函数;
S4,根据建立的解析函数计算给定脉冲速度增量下的轨道根数改变量,完成高阶解析计算;
步骤S3具体为:
S3.1,采用泰勒展开法得到轨道根数变化关于脉冲速度增量的N阶形式表达式;
S3.2,将表达式中的连续力和时间脉冲按照下述方式结合形成脉冲速度增量;
S3.3,将N阶形式表达式的连续力表达式变换为脉冲速度增量表达式,完成了N阶IGVEs的构建;
具体为:
S3.1,采用泰勒展开法得到轨道根数变化关于脉冲速度增量的N阶形式表达式:
Figure FDA0002940589110000011
S3.2,将表达式中的连续力和时间脉冲按照下述方式结合形成脉冲速度增量:
Δvr=Fr·Δt, Δvt=Ft·Δt, Δvn=Fn·Δt (4)
S3.3,考虑公式(4)中的等价关系,将公式(3)中的连续力表达式变换为脉冲速度增量表达式,从而得到:
Figure FDA0002940589110000021
其中,ηkpq是经整理后得到的系数;完成了N阶IGVEs的构建。
2.根据权利要求1所述的脉冲速度增量改变轨道根数的高阶解析计算方法,其特征在于,步骤S2具体为:
S2.1,给出轨道根数的微分形式的高斯变分方程;
S2.2,采用微积分的链式法则,对上述微分形式的高斯变分方程进一步求导。
3.根据权利要求1或2所述的脉冲速度增量改变轨道根数的高阶解析计算方法,其特征在于,步骤S2具体为:
S2.1,给出微分形式的高斯变分方程如下:
Figure FDA0002940589110000022
Figure FDA0002940589110000023
Figure FDA0002940589110000024
Figure FDA0002940589110000025
Figure FDA0002940589110000026
Figure FDA0002940589110000027
其中,等式左边的t为时间,a、e、i、Ω、ω、f为轨道根数,p=a(1-e2)、
Figure FDA0002940589110000031
为与轨道根数相关的两个轨道参数,μ为地球引力系数,Fr、Ft、Fn为连续控制力;具体地:
a、轨道半长轴;e、轨道偏心率;i、轨道倾角;Ω、升交点赤经;ω、纬度幅角;f、真近点角;p、轨道通径;h、轨道角动量;Fr、径向控制力;Ft、切向控制力;Fn、法向控制力;
S2.2,采用微积分的链式法则,对上述微分形式的高斯变分方程进一步求导:对任意六个轨道根数中的一个αk,若已得到其1~m-1阶的导数,则按照下述公式求解其第m阶的导数:
Figure FDA0002940589110000032
其中,
Figure FDA0002940589110000033
的表达式采用微分型高斯变分方程;
当m=N时,终止上述求导过程,否则继续。
CN201910621394.5A 2019-07-10 2019-07-10 一种脉冲速度增量改变轨道根数的高阶解析计算方法 Active CN110489791B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910621394.5A CN110489791B (zh) 2019-07-10 2019-07-10 一种脉冲速度增量改变轨道根数的高阶解析计算方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910621394.5A CN110489791B (zh) 2019-07-10 2019-07-10 一种脉冲速度增量改变轨道根数的高阶解析计算方法

Publications (2)

Publication Number Publication Date
CN110489791A CN110489791A (zh) 2019-11-22
CN110489791B true CN110489791B (zh) 2021-06-08

Family

ID=68546913

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910621394.5A Active CN110489791B (zh) 2019-07-10 2019-07-10 一种脉冲速度增量改变轨道根数的高阶解析计算方法

Country Status (1)

Country Link
CN (1) CN110489791B (zh)

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106570285A (zh) * 2016-11-09 2017-04-19 中国人民解放军装备学院 一种基于状态传递矩阵解析解的J2摄动Lambert问题求解方法
CN108875175A (zh) * 2018-06-06 2018-11-23 北京航空航天大学 一种高阶非中心引力场下的不变相对轨道初值确定方法

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102230969B (zh) * 2011-03-22 2013-05-01 航天恒星科技有限公司 一种卫星星座星间链路的长时间自主维持方法
US20180241497A1 (en) * 2017-02-22 2018-08-23 Electronics And Telecommunications Research Instit Ute Planar electromagnetic wave generation apparatus for concentrating orbital angular momentum and method therefor

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106570285A (zh) * 2016-11-09 2017-04-19 中国人民解放军装备学院 一种基于状态传递矩阵解析解的J2摄动Lambert问题求解方法
CN108875175A (zh) * 2018-06-06 2018-11-23 北京航空航天大学 一种高阶非中心引力场下的不变相对轨道初值确定方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
Definition and analysis of properties of optimal three-impulse point-to-orbit transfers under time constraint;Ivashkin V V 等;《Acta Astronautica》;19810131;第8卷(第1期);第11-23页 *
轨道根数偏差修正的脉冲控制比较研究;任鑫东 等;《计算机仿真》;20101031;第36卷(第5期);第69-72页 *

Also Published As

Publication number Publication date
CN110489791A (zh) 2019-11-22

Similar Documents

Publication Publication Date Title
Valsecchi et al. Resonant returns to close approaches: Analytical theory
CN109255096B (zh) 一种基于微分代数的地球同步卫星轨道不确定演化方法
CN112257343B (zh) 一种高精度地面轨迹重复轨道优化方法及系统
CN101226062B (zh) 一种星上实时计算环月轨道的方法
CN110489791B (zh) 一种脉冲速度增量改变轨道根数的高阶解析计算方法
CN111044082A (zh) 一种基于星敏感器辅助的陀螺误差参数在轨快速标定方法
CN102288177B (zh) 一种基于角速率输出的捷联系统速度解算方法
CN105182906A (zh) 基于高阶s型运动轨迹的位置与速度控制方法
CN104316048A (zh) 一种普适性的脉冲星自主导航测量模型构建方法
CN101226206A (zh) 一种对加速度计零位偏差进行在轨标定的方法
CN107782307A (zh) 一种sins/dr组合导航系统里程计异常数据后处理方法
Zlenko A celestial-mechanical model for the tidal evolution of the Earth-Moon system treated as a double planet
CN105651287B (zh) 一种太阳耀斑到达时间差分测量及组合导航方法、系统
CN105334731A (zh) 一种远地点点火高精度解析轨道自主预报方法
CN107576441B (zh) 动量轮摩擦力矩在轨实时估计方法
Gautier et al. Joint stiffness identification from only motor force/torque data
CN111428912B (zh) 一种基于支持向量机的火星探测器轨道预测方法及系统
Verma et al. Prospects of dynamical determination of General Relativity parameter beta and solar quadrupole moment J2 with asteroid radar astronomy
CN105004338A (zh) 轨道动力学辅助动态脉冲星信号的处理方法
CN106325099A (zh) 一种基于伪相对运动的航天器实时轨道改进方法
CN112327333B (zh) 一种卫星位置计算方法及装置
CN102935898A (zh) 空间纯引力轨道万有引力摄动在轨飞行验证方法
CN111351483B (zh) 一种递归多子样大动态惯性导航方法
CN106153052A (zh) 一种脉冲星的脉冲模板的构建方法
CN116176869B (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