一种解析的小推力圆轨道异面交会优化方法与流程

1.本发明属于航天导航控制技术领域，具体涉及一种解析的小推力圆轨道异面交会优化方法。


背景技术：

2.在空间操作或碎片清除任务中，航天器需要以最省燃料的方式交会目标并开展各类操作。电推进相比传统化学推进具有高比冲的优势，但推力较小，轨道转移时间较长，优化计算存在一定困难。间接法优化方式由于仅能满足最优性必要条件，因此求解容易陷入局部最优，且反复积分轨道耗时很长。现有的参数优化方法同样存在积分耗时，求解变量过多的问题，不容易获得最优交会轨迹。


技术实现要素：

3.本发明的目的是提供一种解析的小推力圆轨道异面交会优化方法，解决了现有技术中存在的参数优化方法积分耗时、求解变量过多的问题，且不容易获得最优交会轨迹的问题。
4.本发明所采用的技术方案是，一种解析的小推力圆轨道异面交会优化方法，具体按照以下步骤实施：
5.步骤1、首先获得输入的航天器初始轨道和需要交会的目标轨道，以及指定的交会时刻，计算航天器无推力情况下在给定的交会时刻，与需要交会的目标轨道的轨道根数差；
6.步骤2、设计一种推力的最优开关机策略，在整个交会飞行时间内分为三个阶段：任务起始阶段的定时开关机，任务中间阶段不开机，临近交会时刻前的定时开关机；定义一个推力开机系数表示第一、三阶段中推力开机时间与第一、三阶段本身时长的比值；得到一种描述小推力圆轨道异面交会轨迹的简化模型；
7.步骤3、根据动力学规律和步骤二定义的开关机策略及模型，建立固定推力开机系数，以步骤二中定义的第一阶段推力产生的半长轴、倾角、升交点赤经改变量为优化变量，以小推力等效速度增量最小为目标的等式约束优化模型；通过对推力开机系数在0到1之间的离散化并分别求解对应的等式约束优化模型，寻优得到最优的推力开机系数和半长轴、倾角、升交点赤经改变量；
8.步骤4、将最优的推力开机系数和半长轴、倾角、升交点赤经改变量转化为小推力实际控制率并代入轨道预报模型中，预报计算交会时刻的航天器轨道根数。得到交会时刻相对于需要交会的目标轨道的相位差；
9.步骤5、根据相位差计算对半长轴改变量的修正量，然后固定半长轴改变量为修正后的值，重新求解步骤3中的等式约束优化模型，得到新的倾角、升交点赤经改变量；
10.步骤6、返回步骤4重新预报交会时刻的轨道根数以及相对目标轨道的相位差，然后判断相位差是否满足给定阈值范围，如果不满足则进入步骤5 继续求解新的半长轴、倾角、升交点赤经改变量并预报相位差；直到相位差小于阈值则输出对应的小推力控制率和
交会转移轨道。
11.本发明的特点还在于，
12.步骤1具体按照以下步骤实施：
13.设输入的初始航天器轨道根数和交会目标轨道根数分别为 [a0,e0,i0,ω0,ω0,m0]和[af,ef,if,ωf,ωf,mf]，其中a为半长轴，e为偏心率，i为倾角，ω为近地点幅角，ω为升交点赤经，m为平近点角，下标0表示初始轨道，f表示目标轨道。则计算二者的差值表示为 [
△
a0,
△
i0,
△
ω0,
△
u0]＝[a
f-a0,i
f-i0,ω
f-ω0,(ωf+m
f-ω
0-m0)]。
[0014]
步骤2具体按照以下步骤实施：
[0015]
步骤2.1、定时开关机是指每圈开机时长与单圈周期的比值固定为一个待优化变量k
thrust
(取值在0到1之间)，推力方向仅包含切向分量和法向分量，且推力方向的切向分量保持不变，法向分量每半圈切换一次符号，
[0016]
u,u
on
,u
off
分别是半圈推力的等效作用中间点相位、开机点相位和关机点相位，有如下关系：设α为推力方向与切向的夹角，则推力f的切向分量为f cosα，法向分量在上半圈为f sinα，下半圈为-f sinα，并假设航天器周期和质量近似保持不变，分别记为常值t和m，则每圈推力对轨道根数的累计变化量写为：
[0017][0018]
其中
△at
表示t时间内推力引起半长轴的改变量，
△it
表示t时间内推力引起倾角的改变量，
△
ω
t
表示t时间内推力引起升交点赤经的改变量。
△vt
,
△
vn分别为每圈推力产生的切向和法向速度增量，a0,i0分别为航天器初始半长轴和倾角，v为初始轨道平均速度。
[0019]
记表示法向速度增量累积的弧段损失系数，并设定义的交会起始阶段定时开关机的总时长为
△
t1，等效中间点相位为u1，临近交会阶段定时开关机的总时长为
△
t2，等效中间点相位为u2，则根据式(2)得
△
t1、
△
t2对应的轨道根数改变量为：
[0020][0021]
其中
△
a1,
△
i1,
△
ω1分别表示交会起始阶段推力对半长轴、倾角和升交点赤经的改变量；
△
a2,
△
i2,
△
ω2分别表示临近交会阶段推力对半长轴、倾角和升交点赤经的改变量；这样就得到了设计的推力开关机策略与轨道根数改变量之间的关系。这个模型将用于下一节的优化模型中，其中
△
a1,
△
i1,
△
ω1，
△
a2,
△
i2,
△
ω2以及k
thrust
，α1，
△
t1，u1，α2，
△
t2，u2均为待优化计算的数值。
[0022]
步骤3具体按照以下步骤实施：
[0023]
步骤3.1、在给定k
thrust
时，为了满足交会需要，根据步骤3计算的轨道根数改变量
△
a1,
△
i1,
△
ω1，
△
a2,
△
i2,
△
ω2满足如下约束：
[0024][0025]
其中，为第一阶段开关机结束后的航天器升交点赤经相对目标轨道升交点赤经的漂移率，j2,re,μ分别为地球非球形摄动j2项、地球半径和引力常数。
[0026]
步骤3.2、对于固定的k
thrust
，建立以
△
a1,
△
i1,
△
ω1,
△
ω2为变量，以推力开机总时间最小为目标的等式约束优化模型：
[0027][0028]
其中参数均已在前面定义。结合式(4)，定义等式约束函数g为：
[0029][0030]
其中为初始轨道对应的升交点赤经漂移率。
[0031]
步骤3.3、令y＝
△
i1,z1＝
△
ω1，z2＝
△
ω2，则根据极小值原理，最优的x,y,z1,z2应满足如下关系：
[0032][0033]
其中λ表示约束乘子，分别表示
△
t1,
△
t2对x,y,z1,z2的导数：
[0034][0035]
求解该方程，式(7)中k
thrust
是固定值，为了得到最优的k
thrust
，分别取k
thrust
为0到1之间等间隔的50个数并求解式(7)，取f最小对应的k
thrust
及x,y,z1,z2为最优解。由式(3)计算得推力开关机控制变量
△
t1,
△
t2,α1,α2,u1,u2。
[0036]
步骤4具体按照以下步骤实施：
[0037]
步骤4.1、将
△
t1,
△
t2,α1,α2,u1,u2代入式(9)对应的小推力精确动力学方程，预报交会时刻轨道根数，动力学方程为：
[0038][0039]
其中，r，v分别为航天器位置速度矢量，可以和轨道根数互转；a
p
为摄动加速度之和，e
t
,en分别为轨道坐标系下的切向和法向单位矢量，计算方法为常识；c为随时间变化的推力开机和关机标志，0为关机，1为开机；α为随时间变化的推力偏航角。c和α的具体取值由
△
t1,
△
t2,α1,α2,u1,u2这几个参数决定，具体关系式如下分为两种情况。
[0040]
当时间t≤t0+
△
t1时：
[0041][0042]
当t≥tf‑△
t2时：
[0043][0044]
式中变量均已定义；
[0045]
步骤4.2、使用数值积分方法根据航天器初始轨道根数和式(9)计算航天器在交会时刻的位置速度以及轨道根数[a
p
,e
p
,i
p
,ω
p
,ω
p
,m
p
]，计算与目标轨道根数的相位差
△up
＝ωf+m
f-ω
p-m
p
。
[0046]
步骤5具体按照以下步骤实施：
[0047]
根据
△up
计算对
△
a1的修正量
[0048][0049]
其中μ为地球引力常数。然后固定重新求解优化模型(6) 中剩余三个变量y,z1,z2；然后根据式(3)就能重新计算参数
△
t1,
△
t2,α1,α2,u1,u2的取值。
[0050]
步骤6具体按照以下步骤实施：
[0051]
将步骤6获得的
△
t1,
△
t2,α1,α2,u1,u2代回到步骤5中，重新预报航天器轨道根数并计算
△up
。
[0052]
判断
△up
是否小于给定阈值，是则输出
△
t1,
△
t2,α1,α2,u1,u2以及对应的预报轨道根数过程中获得的航天器飞行中随时间变化的位置速度以及轨道根数，也就是小推力最优交会轨迹。否则返回步骤6继续对
△
a1进行修正，并继续计算新的
△up
，直到
△up
小于给定阈值。
[0053]
本发明的有益效果是，一种解析的小推力圆轨道异面交会优化方法，首先计算航天器无机动时与目标在交会时刻的轨道根数差；然后设计了一种解析的小推力开关机模式，建立小推力开机参数的等式约束优化模型；通过求解非线性方程组并进行迭代和局部寻优，获得小推力交会的最优开机参数和转移轨迹。本发明公开的方法解决了求解长时间小推力轨道交会耗时长，计算量大的问题，适用圆轨道异面交会的轨迹优化，求解效率显著提高。
附图说明
[0054]
图1是推力开关机策略示意图；
[0055]
图2是本发明方法的流程图；
[0056]
图3是交会飞行过程中推力随时间的变化图；
[0057]
图4是交会过程半长轴差随时间的变化图；
[0058]
图5是交会过程倾角差随时间的变化图；
[0059]
图6是交会过程升交点赤经差随时间的变化图；
[0060]
图7是交会过程相位差变随时间的化图。
具体实施方式
[0061]
下面结合附图和具体实施方式对本发明进行详细说明。
[0062]
针对多圈小推力轨道转移的特点，建立了一种参数较少，使用等式约束优化模型即可求解的快速优化方法，解决了长时间轨道交会的快速求解问题。为了实现上述技术目的，本发明所采用的技术方案是：一种解析的小推力圆轨道异面交会优化方法，即给定航天器初始轨道根数[a0,e0,i0,ω0,ω0,m0] 以及目标轨道[af,ef,if,ωf,ωf,mf]，给定交会时刻tf和时长
△
t，固定推力为f 的圆轨道交会，获得最优的推力开关机控制率和转移轨迹(交会过程中的位置速度变化)的过程。
[0063]
一种解析的小推力圆轨道异面交会优化方法，结合图1～2，具体按照以下步骤实施：
[0064]
步骤1、首先获得输入的航天器初始轨道和需要交会的目标轨道，以及指定的交会时刻，计算航天器无推力情况下在给定的交会时刻，与需要交会的目标轨道的轨道根数差；
[0065]
步骤1具体按照以下步骤实施：
[0066]
设输入的初始航天器轨道根数和交会目标轨道根数分别为 [a0,e0,i0,ω0,ω0,m0]和[af,ef,if,ωf,ωf,mf]，其中a为半长轴，e为偏心率，i为倾角，ω为近地点幅角，ω为升交点赤经，m为平近点角，下标0表示初始轨道，f表示目标轨道。则计算二者的差值表示为[
△
a0,
△
i0,
△
ω0,
△
u0]＝[a
f-a0,i
f-i0,ω
f-ω0,(ωf+m
f-ω
0-m0)]。
[0067]
步骤2、设计一种推力的最优开关机策略，在整个交会飞行时间内分为三个阶段：任务起始阶段的定时开关机，任务中间阶段不开机，临近交会时刻前的定时开关机；定义一
个推力开机系数表示第一、三阶段中推力开机时间与第一、三阶段本身时长的比值；得到一种描述小推力圆轨道异面交会轨迹的简化模型；
[0068]
步骤2具体按照以下步骤实施：
[0069]
步骤2.1、定时开关机是指每圈开机时长与单圈周期的比值固定为一个待优化变量k
thrust
(取值在0到1之间)，推力方向仅包含切向分量和法向分量，且推力方向的切向分量保持不变，法向分量每半圈切换一次符号，如附图1所示。附图1中u,u
on
,u
off
分别是半圈推力的等效作用中间点相位、开机点相位和关机点相位，有如下关系：设α为推力方向与切向的夹角，则推力f的切向分量为f cosα，法向分量在上半圈为f sinα，下半圈为-f sinα，并假设航天器周期和质量近似保持不变，分别记为常值t 和m，则每圈推力对轨道根数的累计变化量写为：
[0070][0071]
其中
△at
表示t时间内推力引起半长轴的改变量，
△it
表示t时间内推力引起倾角的改变量，
△
ω
t
表示t时间内推力引起升交点赤经的改变量。
△vt
,
△
vn分别为每圈推力产生的切向和法向速度增量，a0,i0分别为航天器初始半长轴和倾角，v为初始轨道平均速度。
[0072]
记表示法向速度增量累积的弧段损失系数，并设定义的交会起始阶段定时开关机的总时长为
△
t1，等效中间点相位为u1，临近交会阶段定时开关机的总时长为
△
t2，等效中间点相位为u2，则根据式(2)得
△
t1、
△
t2对应的轨道根数改变量为：
[0073][0074]
其中
△
a1,
△
i1,
△
ω1分别表示交会起始阶段推力对半长轴、倾角和升交点赤经的改变量；
△
a2,
△
i2,
△
ω2分别表示临近交会阶段推力对半长轴、倾角和升交点赤经的改变量；这样就得到了设计的推力开关机策略与轨道根数改变量之间的关系。这个模型将用于下一节的优化模型中，其中
△
a1,
△
i1,
△
ω1，
△
a2,
△
i2,
△
ω2以及k
thrust
，α1，
△
t1，u1，α2，
△
t2，u2均为待优化计算的数值。
[0075]
步骤3、根据动力学规律和步骤二定义的开关机策略及模型，建立固定推力开机系数，以步骤二中定义的第一阶段推力产生的半长轴、倾角、升交点赤经改变量为优化变量，以小推力等效速度增量最小为目标的等式约束优化模型；通过对推力开机系数在0到1之间的离散化并分别求解对应的等式约束优化模型，寻优得到最优的推力开机系数和半长轴、倾角、升交点赤经改变量；
[0076]
步骤3具体按照以下步骤实施：
[0077]
步骤3.1、在给定k
thrust
时，为了满足交会需要，根据步骤3计算的轨道根数改变量
△
a1,
△
i1,
△
ω1，
△
a2,
△
i2,
△
ω2满足如下约束：
[0078][0079]
其中，为第一阶段开关机结束后的航天器升交点赤经相对目标轨道升交点赤经的漂移率，j2,re,μ分别为地球非球形摄动j2项、地球半径和引力常数。
[0080]
步骤3.2、对于固定的k
thrust
，建立以
△
a1,
△
i1,
△
ω1,
△
ω2为变量，以推力开机总时间最小为目标的等式约束优化模型：
[0081]
[0082]
其中参数均已在前面定义。结合式(4)，定义等式约束函数g为：
[0083][0084]
其中为初始轨道对应的升交点赤经漂移率。
[0085]
步骤3.3、令y＝
△
i1,z1＝
△
ω1，z2＝
△
ω2，则根据极小值原理，最优的x,y,z1,z2应满足如下关系：
[0086][0087]
其中λ表示约束乘子，分别表示
△
t1,
△
t2对x,y,z1,z2的导数：
[0088][0089]
求解该方程，式(7)中k
thrust
是固定值，为了得到最优的k
thrust
，分别取k
thrust
为0到1之间等间隔的50个数并求解式(7)，取f最小对应的k
thrust
及x,y,z1,z2为最优解。由式(3)计算得推力开关机控制变量
△
t1,
△
t2,α1,α2,u1,u2。
[0090]
步骤4、将最优的推力开机系数和半长轴、倾角、升交点赤经改变量转化为小推力实际控制率并代入轨道预报模型中，预报计算交会时刻的航天器轨道根数。得到交会时刻相对于需要交会的目标轨道的相位差；
[0091]
步骤4具体按照以下步骤实施：
[0092]
步骤4.1、将
△
t1,
△
t2,α1,α2,u1,u2代入式(9)对应的小推力精确动力学方程，预报交会时刻轨道根数，动力学方程为：
[0093][0094]
其中，r，v分别为航天器位置速度矢量，可以和轨道根数互转；a
p
为摄动加速度之和，e
t
,en分别为轨道坐标系下的切向和法向单位矢量，计算方法为常识；c为随时间变化的推力开机和关机标志，0为关机，1为开机；α为随时间变化的推力偏航角。c和α的具体取值由
△
t1,
△
t2,α1,α2,u1,u2这几个参数决定，具体关系式如下分为两种情况。
[0095]
当时间t≤t0+
△
t1时：
[0096][0097]
当t≥tf‑△
t2时：
[0098][0099]
式中变量均已定义；
[0100]
步骤4.2、使用数值积分方法根据航天器初始轨道根数和式(9)计算航天器在交会时刻的位置速度以及轨道根数[a
p
,e
p
,i
p
,ω
p
,ω
p
,m
p
]，计算与目标轨道根数的相位差
△up
＝ωf+m
f-ω
p-m
p
。
[0101]
步骤5、根据相位差计算对半长轴改变量的修正量，然后固定半长轴改变量为修正后的值，重新求解步骤3中的等式约束优化模型，得到新的倾角、升交点赤经改变量；
[0102]
步骤5具体按照以下步骤实施：
[0103]
根据
△up
计算对
△
a1的修正量
[0104][0105]
其中μ为地球引力常数。然后固定重新求解优化模型(6) 中剩余三个变量y,z1,z2；然后根据式(3)就能重新计算参数
△
t1,
△
t2,α1,α2,u1,u2的取值。
[0106]
步骤6、返回步骤4重新预报交会时刻的轨道根数以及相对目标轨道的相位差，然后判断相位差是否满足给定阈值范围，如果不满足则进入步骤5 继续求解新的半长轴、倾角、升交点赤经改变量并预报相位差；直到相位差小于阈值则输出对应的小推力控制率和交会转移轨道。
[0107]
步骤6具体按照以下步骤实施：
[0108]
将步骤6获得的
△
t1,
△
t2,α1,α2,u1,u2代回到步骤5中，重新预报航天器轨道根数并计算
△up
。
[0109]
判断
△up
是否小于给定阈值，是则输出
△
t1,
△
t2,α1,α2,u1,u2以及对应的预报轨道根数过程中获得的航天器飞行中随时间变化的位置速度以及轨道根数，也就是小推力最优交会轨迹。否则返回步骤6继续对
△
a1进行修正，并继续计算新的
△up
，直到
△up
小于给定阈值。
[0110]
实施例
[0111]
设航天器初始质量为500kg，比冲为1000s，推力为0.5n。要求在14天后交会目标轨道。本星轨道和目标轨道根数如下表
[0112]
表1历元为初始时刻的本星和目标轨道
[0113] a(km)ei(度)ω(度)ω(度)m(度)本星7100098000目标7140098.1100
[0114]
求解步骤如下(1)首先计算交会时刻本星和目标轨道差为 [
△
a0,
△ex0
,
△ey0
,
△
i0,
△
ω0,
△
u0]＝[40km，0，0，0.1度，1度，0，0]。(2)将k
thrust
离散为1到0之间间隔0.02的51个数，依次求解式(7)方程组，得到k
thrust
＝0.38 时等效速度增量93.78m/s，为最优值。(3)计算得到u1,α1,u2,α2,
△
t1,
△
t2分别为[0.7931，2.0756，2.0050，0.9559，1.3019天，1.5547天]。(4)预报计算得交会时刻轨道根数与目标的轨道根数存在相位差
△up
＝-0.8076。(5)按照式(12)计算修正后的x＝x+
△a1c
，重新求解式(7)，得到使
△up
＝0的 [u1,α1,u2,α2,
△
t1,
△
t2]＝[0.7961，2.0255，2.0141，0.9181，1.3072天，1.5300天]。
[0115]
(6)按照式(9)输出最优解对应的轨道根数变化和推力开关机详细信息。其中推力大小和方向的变化趋势见附图3，航天器相对目标的轨道根数差的变化见附图4～7。与其他方法比较，本发明优化耗时极大减小(本发明中方法耗时仅为1秒以内，传统的间接法约20分钟以上，且不能稳定收敛到全局最优解)；二者燃料消耗仅相差5％以内；并且，本发明对应的推力方向角α在区间内为常值，相比间接法计算出的推力方向角是连续变化的，更容易工程实现。
[0116]
文献[1]罗亚中,黄岸毅,李恒年,伍升钢,张进,杨震.一种j2摄动下长时间轨道交会最优速度增量快速估计方法[p].湖南省：cn110789739a,2020-02-14。