一种航天器轨道的并行自适应计算方法及系统

1.本发明涉及航空航天技术领域，具体为一种航天器轨道的并行自适应计算方法及系统。


背景技术：

2.在航空航天领域，许多问题都可以用非线性常微分方程的数学模型来描述。因此，在许多航空航天的工程任务中，非线性常微分方程的求解效率严重影响着任务的执行效果。近年来，一类结合了配点法和迭代法的非线性常微分方程求解方法被广泛应用到航天器轨道递推，轨道转移等航天工程任务中。配点迭代类方法具有比传统有限差分类方法更高的计算精度和更快的计算速度，但该方法的计算效率通常会受到计算参数的严重影响，选择合适的计算参数可以避免此类方法计算效率低下。但现有配点迭代类算法的自适应参数选择策略都是针对串行方法设计的，串行的配点迭代类方法因为串行计算节点不足会出现运算操作排队等待计算节点，从而造成计算效率降低的问题。


技术实现要素：

3.针对现有技术中存在的问题，本发明提供一种航天器轨道的并行自适应计算方法及系统，以提高航天器轨道的计算效率。
4.本发明是通过以下技术方案来实现：
5.一种航天器轨道的并行自适应计算方法，其特征在于，包括以下步骤：
6.步骤1，选定航天器轨道的总计算区间和总区间的配点数m，以及设定初始子区间的步长；
7.步骤2、采用并行的多个配点迭代算法分别对各配点时刻的初始状态进行迭代修正，并且各配点迭代算法采用不同的初始子区间步长，根据各配点迭代算法的收敛状态和迭代次数的关系，确定下一子区间的步长；
8.步骤3、根据确定的下一子区间的步长，重复步骤2，直至完成所有子区间状态的迭代修正，得到各子区间的状态修正值，进而得到航天器轨道。
9.优选的，步骤2中，采用gpu实现多个配点迭代算法的并行运行，具体如下：
10.在gpu上分配多个相互独立且并行的cuda流，各cuda流运行不同初始子区间步长的配点迭代算法。
11.优选的，所述配点迭代算法修正子区间状态的方法如下：
12.采用线性运算函数将配点迭代计算方法重写为误差迭代方程和误差迭代方程；
13.将误差迭代方程和状态迭代方程结合常矩阵对子区间中各配点时刻的初始状态进行迭代修正，直至状态修正值达到预定精度，得到该子区间的航天器轨道。
14.优选的，所述误差迭代方程的表达式如下：
[0015][0016]
状态迭代方程如下：
[0017][0018]
其中，为第n次迭代修正后航天器的位置及速度估计值，和分别为和中第d维状态对应的向量。
[0019]
优选的，步骤2中根据收敛状态和迭代次数的关系，确定下一子区间的步长的方法如下；
[0020]
s个配点迭代算法的步长分别为n1×
dt,n2×
dt,
……
，n
s-1
×
dt；
[0021]
其中，1＝n1＜n2＜
…
＜n
s-1
，s≥2；
[0022]
按照步长自小到大对各配点迭代算法进行排序，当排序中有且仅有前c个配点迭代算法收敛时，即0＜c≤s；
[0023]
获取各配点迭代算法的迭代次数分别为i1，i2……
，ic，根据迭代次数缩小下一子区间的步长。
[0024]
优选的，所述下一子区间的步长的计算方法如下：
[0025]
当时：
[0026]
dt＝nm×
dt，
[0027]
当时：
[0028]
dt＝dt/n
m+1
，
[0029]
优选的，当所有的配点迭代算法均未收敛，缩短当前子区间的步长，并采用修改后的步长对当前子区间的初始状态重新迭代修正。
[0030]
优选的，所述修改后的步长如下：
[0031]
dt＝dt/n
s-1
。
[0032]
上述航天器轨道的并行自适应计算方法的系统，包括
[0033]
初始化模块，用于选定航天器轨道的总计算区间和总区间的配点数m，以及设定初始子区间的步长；
[0034]
并行迭代模块，用于采用并行的多个配点迭代算法分别对各配点时刻的初始状态进行迭代修正，并且各配点迭代算法采用不同的初始子区间步长，根据各配点迭代算法的收敛状态和迭代次数的关系，确定下一子区间的步长；
[0035]
轨道模块，用于根据确定的下一子区间的步长，并结合并行的多个配点迭代算法完成所有子区间状态的迭代修正，输出航天器轨道。
[0036]
与现有技术相比，本发明具有以下有益的技术效果：
[0037]
本发明公开了一种航天器轨道的并行自适应计算方法，采用并行的多个配点迭代算法分别对各配点时刻的初始状态进行迭代修正，并且各配点迭代算法采用不同的初始子区间步长，根据各配点迭代算法的收敛状态和迭代次数的关系，确定下一子区间的步长，然后根据确定的步长重复上述过程完成航天器轨道的计算，本发明根据配点迭代类方法在不同计算步长下的收敛情况和迭代次数来调整计算的步长，使得在计算过程中子区间的步长
总是适当的，避免了配点迭代类方法因为计算步长过长或过短而导致计算精度或效率降低的问题，充分发挥了配点迭代类方法大步长计算的优势。
附图说明
[0038]
图1为本发明航天器轨道的自适应计算方法原理图；
[0039]
图2为本发明航天器轨道的自适应计算方法流程图；
[0040]
图3为本发明航天器轨道计算结果图；
[0041]
图4为本发明自适应并行计算方法所带来的加速比图。
具体实施方式
[0042]
下面结合附图对本发明做进一步的详细说明，所述是对本发明的解释而不是限定。
[0043]
如图1和图2所示，一种航天器轨道的并行自适应计算方法，包括如下步骤：
[0044]
步骤1，选定航天器轨道的总计算区间和总区间的配点数m，以及设定初始子区间的步长；
[0045]
根据经验选定固定配点数，基于gpu的并行计算方法的配点数通常取为32的整数倍，将总计算区间[t0,tf]分为若干个子区间，并设定初始子区间的步长为dt，即该子区间的长度。
[0046]
其中，t0为总计算区间的起始时间，tf为总计算区间的终止时间。
[0047]
在总计算区间[t0,tf]范围内根据公式选取m个cgl节点，其中τm＝-cos((m-1)π/(m-1)),m＝1,2,
…
,m。
[0048]
步骤2，采用并行的多个配点迭代算法分别对各配点时刻的初始状态进行迭代修正，并且各配点迭代算法采用不同的初始子区间步长，根据各配点迭代算法的收敛状态和迭代次数的关系，确定对下一子区间的步长，具体如下：
[0049]
s2.1、根据配点数确定线性运算函数；
[0050]
利用cublas库中的一般矩阵乘法子函数cublasdgemm，对角矩阵乘法子函数cublasddgmm和gpu核函数生成所需要的常矩阵q,p,h和这些常矩阵可预先生成并保存在gpu端备用；
[0051]
切比雪夫基函数微分常方阵q的表达式如下：
[0052][0053]
其中，φ(τ)＝[φ0(τ),φ1(τ),...,φn(τ)]为n+1个chebyshev基函数在τ处的值组成的向量；
[0054]
切比雪夫基函数积分常方阵p的表达式如下：
[0055][0056]
其中，φ0(τ)＝1,φ1(τ)＝cos(arccost),...,φn(τ)＝cos((n-1)arccost)；
[0057]
若τ＝
±
1，则1，则
[0058][0059]
忽略小量只保留常数的常雅可比矩阵的表达式如下：
[0060][0061]
其中，0
3m
×
3m
为3m
×
3m的全零矩阵，i
3m
×
3m
为3m
×
3m的单位矩阵；
[0062]
只和配点数有关的常矩阵h的表达式如下：
[0063]
h＝pt-tp，t＝dig(τ1,τ2,
…
,τm)
[0064]
s2.2、根据冷启动方案确定初始子区间的初始状态即令初始运动为速度为初速度的匀速直线运动，此过程利用gpu核函数实现，同时利用核函数得到初始状态的一阶导数。
[0065]
s2.3、采用线性运算函数将配点迭代计算方法重写为误差迭代方程和误差迭代方程。
[0066]
误差迭代方程如下：
[0067][0068]
状态迭代方程如下：
[0069][0070]
其中，为第n次迭代修正后航天器的位置及速度估计值，和分别为和中第d维状态对应的向量。
[0071]
s2.4、在gpu上分配多个相互独立且并行的cuda流，并对各cuda流设定不同的初始子区间的步长。
[0072]
s2.5、各cuda流根据各自子区间的步长，将误差迭代方程和状态迭代方程并结合常矩阵对子区间中各配点时刻的初始状态进行迭代修正；
[0073]
根据切比雪夫基函数微分常方阵q、切比雪夫基函数积分常方阵p、忽略小量只保留常数的常雅可比矩阵与配点数有关的常矩阵h以及初始子区间的初始状态结合配点迭代方法对初始子区间中各配点时刻的状态进行迭代修正，直至各cuda流均完成迭代，得到初始子区间的状态修正值。
[0074]
s2.6、根据各cuda流的配点迭代算法的收敛状态和迭代次数，确定各cuda流的下一子区间的步长。
[0075]
设共有s个cuda流，其计算步长分别为n1×
dt,n2×
dt,
……
，n
s-1
×
dt，其中1＝n1＜n2＜
…
＜n
s-1
，s≥2。设s个cuda流中有且仅有前c个cuda流收敛，0≤c≤s。若c＝0，则令步长dt＝dt/n
s-1
。若c≠0，则记各cuda流对应的迭代次数分别为i1，i2……
，ic，当时，令步长dt＝nm×
dt，当时，令步长dt＝dt/n
m+1
，其中，
[0076]
步骤3、重复步骤s2.5-2.6，直至完成总计算区间内各子区间均收敛。
[0077]
实施例1
[0078]
在gpu上分配两个相互独立的cuda流，cuda流1和cuda流2，两个cuda流上可以同时运行不同的并行程序，两个cuda流采用不同的初始子区间步长，采用配点迭代方法对初始子区间中各配点时刻的初始状态进行迭代修正。
[0079]
在本实施例中重点对根据当前子区间的迭代状态对下一子区间的步长设置方法进行阐述。
[0080]
s2.4、在cuda流1采用dt为初始子区间的步长，在cuda流2上采用2
×
dt为初始子区间的步长。
[0081]
s2.5、各cuda流根据各自初始子区间的步长，将误差迭代方程和状态迭代方程并结合常矩阵对初始子区间中各配点时刻的初始状态进行迭代修正；
[0082]
s2.6、根据依据cuda流1和cuda流2的配点迭代算法的收敛状态和迭代次数，确定各cuda流的下一子区间的步长的变化情况，具体如下：
[0083]
确定步长最大的cuda流2的迭代结果是否收敛，将该收敛结果记为第一判断结果；
[0084]
若第一判断结果为是，由于cuda流1的步长小于cuda流2，因此cuda流1必然收敛，记录cuda流2的计算结果，令ti＝ti+2
×
dt，ti为子区间的初始时刻，ti增加代表当前子区间迭代结束，开始下一子区间的计算，根据第一判断结果并结合各cuda流的迭代次数确定下一子区间的步长。
[0085]
记cuda流1的迭代次数为：i1；
[0086]
记cuda流2的迭代次数为：i2；
[0087]
判断i2是否大于2倍的i1，结果记为第二判断；
[0088]
若第二判断的结果为是，则令下一子区间的步长dt＝dt/2；
[0089]
若第二判断的结果为否，则令下一子区间的步长dt＝2
×
dt；
[0090]
若cuda流1计算收敛，且cuda流2计算未收敛，记录cuda流1的算结果，将步长增倍，然后进行下一区间的计算。
[0091]
若cuda流1和cuda流2均未收敛，将步长减半后在本区间重新进行计算。
[0092]
若第一判断的结果为否，判断cuda流1计算是否收敛，结果记为第三判断；
[0093]
若第三判断的结果为是，则记录cuda流1的计算结果，并令ti＝ti+dt，dt＝dt/2；
[0094]
若第三判断的结果为否，则令下一子区间的步长dt＝dt/2；
[0095]
步骤s3、根据步骤2得到的步长，并重复步骤s2.5-s2.6对下一子区间进行迭代计算，直至完成所有子区间的计算，得到航天器的轨道。
[0096]
当下一子区间的起始时刻ti等于总计算区间的终止时间，则表明总计算区间的轨道计算均已完成。
[0097]
实施例2
[0098]
一种航天器轨道的并行自适应计算方法，包括以下步骤：
[0099]
步骤1、对于卫星常见的轨道预测问题，根据运行的硬件支持选择固定配点数m，迭代精度误差限ε，允许的最大迭代次数，初始子区间步长dt和整体计算区间[t0,tf]，将总计算区间划分为多个连续的子区间，依据选定的固定配点数确定配点时刻tm，构造与配点数对应的常矩阵q,p,h和
[0100]
步骤2、在gpu上分配两个相互独立的cuda流，cuda流1和cuda流2，在cuda流1上以dt为计算步长，m为配点数进行迭代计算，在cuda流2上以2倍的dt为计算步长，m为配点数进行迭代计算。
[0101]
具体的计算过程为：
[0102]
利用gpu核函数生成初始子区间的初始状态由初始状态开始进行迭代计算。将状态值代入核函数中计算当前状态的一阶导数并利用公式和由当前状态值和一阶导数计算下一次迭代修正后的状态值
[0103]
是一个6m维的列向量。
[0104]
其中向量第1至第m个元素对应航天器在各配点处x轴方向的位置估计值，向量第m+1至第2m个元素对应航天器在各配点处y轴方向的位置估计值，向量第2m+1至第3m个元素对应航天器在各配点处z轴方向的位置估计值，向量第3m+1至第4m个元素对应航天器在各配点处x轴方向的速度估计值，向量第4m+1至第5m个元素对应航天器在各配点处y轴方向的速度估计值，向量第5m+1至第6m个元素对应航天器在各配点处z轴方向的速度估计值。
[0105]
同理，也是一个6m维的列向量。
[0106]
其中，向量第1至第3m个元素对应航天器在各配点处三个方向的速度估计值，向量第3m+1至第6m个元素对应航天器在各配点处各个方向的重力场估计值。
[0107]
依据式计算第n+1次与第n次迭代修正后航天器各配点处状态估计值差值向量的二范数e，若e＞ε，则进行下一次迭代计算，直到e≤ε为止，记录计算结果并退出迭代计算；在迭代计算过程中，若迭代次数大于所设置的最大迭代次数后仍未收敛，则表明此次迭代计算不收敛，退出迭代计算。
[0108]
步骤3、在每个子区间的迭代计算中，依据cuda流1和cuda流2的收敛情况和迭代次数的数量关系，判断确定下一子区间的步长应增倍还是减半。
[0109]
若cuda流1和cuda流2计算均收敛，记录cuda流2的计算结果，并结合cuda流1和cuda流2的迭代次数确定下一子区间的步长。
[0110]
若cuda流2的迭代次数i2大于2倍的cuda流2的迭代次数i1，则令步长dt＝dt/2，反之，则令步长dt＝2
×
dt。
[0111]
若cuda流1计算收敛，且cuda流2计算未收敛，记录cuda流1的算结果，令步长dt＝dt/2，然后进行下一区间的计算。
[0112]
若cuda流1和cuda流2均未收敛，令步长dt＝dt/2，对当前子区间的状态重新迭代
计算。
[0113]
该航天器轨道的并行自适应计算方法，选定固定配点数m,初始计算步长dt和总计算区间；分配cuda流1以dt为步长，cuda流2以2倍的dt为步长，使用并行的配点迭代类方法进行迭代计算；依据cuda流1和cuda流2的收敛情况和迭代次数的数量关系，判断计算步长应增倍还是减半；依据cuda流1和cuda流2的收敛情况，确定是否进行下一区间的迭代计算，直至子区间完成迭代。本发明使用gpu并行加速技术，实现了基于gpu的反馈加速picard迭代方法，通过不断迭代修正当前解使其逐渐逼近真实解，并将迭代计算量分摊到gpu上大量的计算节点上来加速计算，为航天器轨道预测提供了一种高效率高精度的并行计算方法。
[0114]
实施例3
[0115]
40阶egm2008地球重力场模型下近地卫星轨道递推问题的动力学方程为：
[0116][0117]
其中,r＝[x,y,z]
t
为近地卫星位置矢量,μ＝3986004.418
×
108m3/s2,表示地球引力常数,r＝||r||，表示地心与近地卫星质心间的距离,为轨道递推初始时刻，a为摄动项，本实施例中考虑40阶egm2008地球重力场模型摄动。
[0118]
本实施例采用初始计算条件如表1，计算参数如表2，计算得到的近地卫星运行轨道如图3所示，通过自适应优化计算参数后的计算精度如图4所示。
[0119]
表1实施例3采用的初始计算条件
[0120][0121][0122]
表2实施例3采用的计算参数
[0123]
航天器运行时间tf/s配点个数迭代精度误差限ε/m最大迭代次数80006410-6
80
[0124]
本发明公开了一种航天器轨道的并行自适应计算方法，所述计算方法包括如下步骤：选定固定配点数m,初始计算步长dt和总计算区间；分配cuda流1以dt为步长，cuda流2以2倍的dt为步长，使用并行的配点迭代类方法进行迭代计算；依据cuda流1和cuda流2的收敛情况和迭代次数的数量关系，判断计算步长应增倍还是减半；依据cuda流1和cuda流2的收敛情况，确定是否进行下一区间的迭代计算，直至总计算区间内各子区间均完成迭代计算，得打航天器的轨道。
[0125]
本发明使用不同cuda流以不同的计算参数进行迭代计算，并根据不同cuda流计算的收敛性和计算迭代次数的数量关系来确定配点迭代法合适的计算参数，为航天器轨道预测提供了一种可自适应选定计算参数的高效率并行计算方法。
[0126]
在另一实施例中，本发明还提供了一种航天器轨道的并行自适应计算方法的系统，包括初始化模块，并行迭代模块和轨道模块。
[0127]
初始化模块，用于选定航天器轨道的总计算区间和总区间的配点数m，以及设定初始子区间的步长；
[0128]
并行迭代模块，用于采用并行的多个配点迭代算法分别对各配点时刻的初始状态进行迭代修正，并且各配点迭代算法采用不同的初始子区间步长，根据各配点迭代算法的收敛状态和迭代次数的关系，确定下一子区间的步长；
[0129]
轨道模块，用于根据确定的下一子区间的步长，并结合并行的多个配点迭代算法完成所有子区间状态的迭代修正，输出航天器轨道。
[0130]
以上内容仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明权利要求书的保护范围之内。

技术特征：
1.一种航天器轨道的并行自适应计算方法，其特征在于，包括以下步骤：步骤1，选定航天器轨道的总计算区间和总区间的配点数m，以及设定初始子区间的步长；步骤2、采用并行的多个配点迭代算法分别对各配点时刻的初始状态进行迭代修正，并且各配点迭代算法采用不同的初始子区间步长，根据各配点迭代算法的收敛状态和迭代次数的关系，确定下一子区间的步长；步骤3、根据确定的下一子区间的步长，重复步骤2，直至完成所有子区间状态的迭代修正，得到各子区间的状态修正值，进而得到航天器轨道。2.根据权利要求1所述的一种航天器轨道的并行自适应计算方法，其特征在于，步骤2中，采用gpu实现多个配点迭代算法的并行运行，具体如下：在gpu上分配多个相互独立且并行的cuda流，各cuda流运行不同初始子区间步长的配点迭代算法。3.根据权利要求1所述的一种航天器轨道的并行自适应计算方法，其特征在于，所述配点迭代算法修正子区间状态的方法如下：采用线性运算函数将配点迭代计算方法重写为误差迭代方程和误差迭代方程；将误差迭代方程和状态迭代方程结合常矩阵对子区间中各配点时刻的初始状态进行迭代修正，直至状态修正值达到预定精度，得到该子区间的航天器轨道。4.根据权利要求3所述的一种航天器轨道的并行自适应计算方法，其特征在于，所述误差迭代方程的表达式如下：状态迭代方程如下：其中，为第n次迭代修正后航天器的位置及速度估计值，和分别为和中第d维状态对应的向量。5.根据权利要求1所述的一种航天器轨道的并行自适应计算方法，其特征在于，步骤2中根据收敛状态和迭代次数的关系，确定下一子区间的步长的方法如下；s个配点迭代算法的步长分别为n1×
dt,n2×
dt,
……
，n
s-1
×
dt；其中，1＝n1＜n2＜
…
＜n
s-1
，s≥2；按照步长自小到大对各配点迭代算法进行排序，当排序中有且仅有前c个配点迭代算法收敛时，即0＜c≤s；获取各配点迭代算法的迭代次数分别为i1，i2……
，i
c
，根据迭代次数缩小下一子区间的步长。6.根据权利要求5所述的一种航天器轨道的并行自适应计算方法，其特征在于，所述下一子区间的步长的计算方法如下：当时：dt＝n
m
×
dt，
当时：dt＝dt/n
m+1
，7.根据权利要求5所述的一种航天器轨道的并行自适应计算方法，其特征在于，当所有的配点迭代算法均未收敛，缩短当前子区间的步长，并采用修改后的步长对当前子区间的初始状态重新迭代修正。8.根据权利要求7所述的一种航天器轨道的并行自适应计算方法，其特征在于，所述修改后的步长如下：dt＝dt/n
s-1
。9.一种用于执行权利要求1-8任一项所述航天器轨道的并行自适应计算方法的系统，其特征在于，包括初始化模块，用于选定航天器轨道的总计算区间和总区间的配点数m，以及设定初始子区间的步长；并行迭代模块，用于采用并行的多个配点迭代算法分别对各配点时刻的初始状态进行迭代修正，并且各配点迭代算法采用不同的初始子区间步长，根据各配点迭代算法的收敛状态和迭代次数的关系，确定下一子区间的步长；轨道模块，用于根据确定的下一子区间的步长，并结合并行的多个配点迭代算法完成所有子区间状态的迭代修正，输出航天器轨道。

技术总结
本发明公开了一种航天器轨道的并行自适应计算方法及系统，采用并行的多个配点迭代算法分别对各配点时刻的初始状态进行迭代修正，并且各配点迭代算法采用不同的初始子区间步长，根据各配点迭代算法的收敛状态和迭代次数的关系，确定下一子区间的步长；并重复迭代过程直至完成所有子区间状态的迭代修正，得到航天器轨道；本发明根据配点迭代类方法在不同计算步长下的收敛情况和迭代次数来调整计算的步长，使得在计算过程中子区间的步长总是适当的，避免了配点迭代类方法因为计算步长过长或过短而导致计算精度或效率降低的问题，充分发挥了配点迭代类方法大步长计算的优势。挥了配点迭代类方法大步长计算的优势。挥了配点迭代类方法大步长计算的优势。


技术研发人员：代洪华 王昌涛 张哲 汪雪川 岳晓奎 刘传凯 孙军
受保护的技术使用者：西北工业大学
技术研发日：2023.06.29
技术公布日：2023/9/20