CN115258196A - 低轨卫星星座组网电推进变轨策略优化方法和系统 - Google Patents

低轨卫星星座组网电推进变轨策略优化方法和系统 Download PDF

Info

Publication number
CN115258196A
CN115258196A CN202210857788.2A CN202210857788A CN115258196A CN 115258196 A CN115258196 A CN 115258196A CN 202210857788 A CN202210857788 A CN 202210857788A CN 115258196 A CN115258196 A CN 115258196A
Authority
CN
China
Prior art keywords
orbit
thrust
electric propulsion
module
networking
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.)
Pending
Application number
CN202210857788.2A
Other languages
English (en)
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.)
Shanghai Institute of Satellite Engineering
Original Assignee
Shanghai Institute of Satellite Engineering
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 Shanghai Institute of Satellite Engineering filed Critical Shanghai Institute of Satellite Engineering
Priority to CN202210857788.2A priority Critical patent/CN115258196A/zh
Publication of CN115258196A publication Critical patent/CN115258196A/zh
Pending legal-status Critical Current

Links

Images

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
    • 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/244—Spacecraft control systems
    • B64G1/245—Attitude control algorithms for spacecraft attitude control
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02T—CLIMATE CHANGE MITIGATION TECHNOLOGIES RELATED TO TRANSPORTATION
    • Y02T90/00—Enabling technologies or technologies with a potential or indirect contribution to GHG emissions mitigation

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Chemical & Material Sciences (AREA)
  • Combustion & Propulsion (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • Automation & Control Theory (AREA)
  • Control Of Position, Course, Altitude, Or Attitude Of Moving Bodies (AREA)

Abstract

本发明提供了一种低轨卫星星座组网电推进变轨策略优化方法和系统，包括：步骤1：根据电推进的推力大小，在仅有开和关两种状态的工程约束的基础上，建立时间最优的星座组网电推进转移轨道优化问题；步骤2：对组网过程进行阶段划分，并将推力大小及方向进行分阶段参数化处理；步骤3：确定待求解变量，将星座组网电推进转移轨道优化问题转换为控制变量受约束的打靶问题；步骤4：对打靶问题的初值进行猜测，采用牛顿迭代算法进行求解，最终得到星座组网电推进转移轨道优化问题的最优解。本发明可以保证转移轨道满足预先设定的组网约束，对姿控系统要求低，结构简单，求解效率高，有利于星上实现，可以有效实现轨道根数的联合控制。

Description

低轨卫星星座组网电推进变轨策略优化方法和系统
技术领域
本发明涉及航天器轨道设计与优化技术领域，具体地，涉及一种低轨卫星星座组网电推进变轨策略优化方法和系统。
背景技术
近年来，高覆盖率、短重访、甚至是全天时持续覆盖的低轨星座相继被提出，成为当前的研究热点。此类低轨星座的规模可达数十颗到数万颗，占据大量的空间轨位资源，且通常配置电推进系统，采用全电推的方式完成组网机动、构型保持、离轨等轨控任务。低轨卫星星座组网机动是指在星箭分离后，卫星从初始轨道变轨爬升至目标轨道的转移过程。与一般的单星任务相比，星座组网对轨道六根数均提出了设计约束，其转移轨道的设计也存在很大的差异。因此，全电推模式下的组网转移轨道设计与优化成了低轨星座发展的关键技术之一。
电推进轨道优化是一个经典的航天动力学设计问题，有很多的学者对其进行了研究，例如小行星探测转移轨道设计、人工太阳同步轨道设计、GTO-GEO多圈小推力转移、椭圆轨道卫星对地覆盖的电推进控制等，而对于低轨星座组网转移轨道设计的公开文献仍比较少。对于低轨卫星星座，为了降低星座中工作的卫星发生碰撞的概率，卫星的工作轨道常选择小偏心率冻结的近圆轨道。因此，低轨星座组网需要解决多个轨道根数联合控制问题，确保在变轨完成后相位差、半长轴、偏心率及近地点幅角均满足要求，且满足一定的性能指标。针对该问题，本专利考虑电推进推力大小固定，且仅有开和关两种状态的工程约束，构建时间最优的星座组网电推进转移轨道优化问题，提出一种推力参数化方法，将原问题转化为一个控制变量受约束的打靶问题，并用非线性规划算法对其进行求解。该方法可以实现相位差、半长轴、偏心率及近地点幅角四者的联合控制，适用于目标轨道是小偏心率冻结的低轨星座任务。最后，采用地球21阶次引力场模型对该方法的有效性进行了校验。
目前，人们对航天器变轨策略进行了一些研究，经检索，与变轨策略设计相关的主要有：
专利《一种静止轨道卫星变轨策略优化方法》(公开号：CN102424116A)克服现有技术的不足，提供了一种静止轨道卫星变轨策略优化方法，合理确定变轨策略设计的各个约束条件，以减少变轨策略设计过程中的人工干预以及计算时间和计算量。该专利关注的静止轨道卫星的变轨策略，模型与优化方法与低轨卫星组网存在很大的差异，其相关设计方法无法借用。
专利《基于粒子群算法的GEO卫星变轨策略计算方法、系统及介质》(公开号：CN108216687A)也是针对静止轨道卫星变轨，将点火时刻和点火方向作为优化变量，将推进剂消耗作为目标函数，通过设定初始粒子种群并按算法进行进化计算，更快速的获得期望解，提高计算效率。该专利同样关注的静止轨道卫星的变轨策略，模型与优化方法与低轨卫星组网存在很大的差异，其相关设计方法无法借用。
专利《基于异面变轨策略的卫星星座重构方法、设备及存储设备》(公开号：CN107885917A)提供了基于异面变轨策略的卫星星座重构方法、设备及存储设备，通过卫星机动变轨的方式，优化卫星的轨道参数，改变卫星星座的空间构型，对卫星星座进行重构，从而满足应急任务对卫星网络组网的性能需求。该专利关注的星座重构方法，而本专利则关注卫星发射入轨初期的组网机动，设计对象不同。
专利《一种卫星的自主变轨方法》(公开号：CN101219713)克服现有技术的不足，提出了一种星上自主进行的姿态确定，姿态机动以及自主轨控的自主变轨方法，解决了深空探测卫星的变轨问题，保证变轨准确、可靠、自动地完成。该专利关注深空探测航天器变轨的自主性，与本专利研究对象不同。
论文《连续小推力条件下星座轨道机动方法研究》(中国空间科学技术，2020年第01期：1-16)研究了一箭多星星座任务的组网轨道机动问题，采用分时升轨机动实现了相位差和半长轴的联合调整，但并未涉及半长轴、偏心率及近地点幅角三者的联合调整，因而无法直接适用于目标轨道为小偏心率冻结的近圆轨道组网任务。而本专利主要关注的是半长轴、偏心率及近地点幅角三者的联合调整。
专利文献CN105799954A(申请号：CN201410845134.3)公开了一种天基分散部署微纳载荷的模块化飞行器及其变轨制导方法，所述飞行器包括：承载与推进一体化单元、测量与控制一体化单元、轻小型级间适配分离装置；可以实现远程自主快速变轨到目标轨道。然而该专利仅对飞行器的变轨进行了研究，与本专利研究对象不同。
发明内容
针对现有技术中的缺陷，本发明的目的是提供一种低轨卫星星座组网电推进变轨策略优化方法和系统。
根据本发明提供的低轨卫星星座组网电推进变轨策略优化方法，包括：
步骤1：根据电推进的推力大小，在仅有开和关两种状态的工程约束的基础上，建立时间最优的星座组网电推进转移轨道优化问题；
步骤2：对组网过程进行阶段划分，并将推力大小及方向进行分阶段参数化处理；
步骤3：确定待求解变量，将星座组网电推进转移轨道优化问题转换为控制变量受约束的打靶问题；
步骤4：对打靶问题的初值进行猜测，采用牛顿迭代算法进行求解，最终得到星座组网电推进转移轨道优化问题的最优解。
优选的，所述步骤1包括：
步骤1.1：在地心惯性参考坐标系下，建立卫星的轨道动力学方程：
Figure BDA0003756199530000031
Figure BDA0003756199530000032
Figure BDA0003756199530000033
其中，r和v分表为位置和速度矢量；fP为包含J2、大气阻力的摄动加速度；μ为地球引力常数；r为地球半径；ge为地球海平面重力加速度；m为航天器质量；Isp为推力器比冲；α为推力方向的单位矢量；Tmax为最大推力；u为实际推力与最大推力的比值；
在仅有开和关两种状态的工程约束的基础上，u的取值为0或者1；
推力方向的单位矢量在惯性坐标系下进一步表示为：
Figure BDA0003756199530000034
其中，
Figure BDA0003756199530000035
为推力方位角，δ为推力高度角；
令卫星的瞬时轨道根数为：
kosc＝[a e i Ω ω M]
其中，a为半长轴，e为偏心率，i为轨道倾角，Ω为升交点赤经，ω为近地点幅角，M为平近点角；
令
Figure BDA0003756199530000036
表示卫星的平均轨道根数，x＝[r v]为直角坐标参数，则kosc、kave及x三者可互相转换；
步骤1.2：令初始时刻为t0，末端时刻为tf，则初末时刻卫星面内的组网约束表示为：
r(t0)＝r0，v(t0)＝v0，m(t0)＝m0
Figure BDA0003756199530000041
Figure BDA0003756199530000042
式中，
Figure BDA0003756199530000043
及
Figure BDA0003756199530000044
为常值，表示末端时刻期望的半长轴、偏心率及近地点幅角，
Figure BDA0003756199530000045
表示沿迹角，定义为：
Figure BDA0003756199530000046
在末端时刻，期望的
Figure BDA0003756199530000047
根据如下公式确定：
Figure BDA0003756199530000048
式中，
Figure BDA0003756199530000049
表示t0时刻的初始，
Figure BDA00037561995300000410
表示沿迹角变化率；
步骤1.3：令控制变量和状态变量：
u＝[u α]T
y＝[r v m]T
则轨道动力学方程表示为：
Figure BDA00037561995300000411
相应的边界约束表示为：
bini[y(t0)]＝y(t0)-y0＝0
Figure BDA00037561995300000412
其中，y0＝[r0 v0 m0]T
组网过程电推进转移轨道问题P为寻找最优控制变量：
Figure BDA00037561995300000413
同时满足：
Figure BDA00037561995300000414
bini[y(t0)]＝0，bter[y(tf)]＝0
u＝0或者u＝1，αTα＝1。
优选的，所述步骤2包括：
步骤2.1：根据最快入轨的性能指标要求和组网轨道的运动特性，将u和α在时间域上按时刻t1、t2及t3进行划分，在推力为全开的基础上，将转移轨道按时间顺序分为巡航阶段、升轨阶段及偏心率矢量调整阶段；
步骤2.2：将各阶段的推力大小和方向设定为：
Figure BDA0003756199530000051
Figure BDA0003756199530000052
式中，αa和αeω分别为升轨阶段和偏心率矢量调整阶段的推力方向单位矢量；
步骤2.3：为满足最快入轨的性能指标要求，在升轨阶段和偏心率矢量调整阶段，将推力方向分别设定沿着切向方向和轨道面内垂直拱线方向，表示为：
Figure BDA0003756199530000053
αeω＝(cosσ)ip+(sinσ)jp
其中，σ表示推力方向与拱线的夹角，且：
Figure BDA0003756199530000054
Figure BDA0003756199530000055
其中，ip和jp的计算是基于平均轨道参数进行的，确保推力方向平缓变化。
优选的，所述步骤3包括：
步骤3.1：对于给定的初始状态t0和y0，定义待求解的变量为：
z＝[Δt1 Δt2 Δt3 σ]T
其中，Δt1＝t1-t0，Δt2＝t2-t1，Δt3＝t3-t2；
步骤3.2：将问题P转换为如下打靶问题：
Figure BDA0003756199530000061
优选的，所述步骤4包括：
步骤4.1：对Δt1、Δt2、Δt3及σ进行猜测，确定迭代变量初值；
步骤4.2：采用龙格库塔法计算代价函数F(z)，采用牛顿迭代算法进行迭代计算，直至||F(z)||收敛至小于指定阈值。
根据本发明提供的低轨卫星星座组网电推进变轨策略优化系统，包括：
模块M1：根据电推进的推力大小，在仅有开和关两种状态的工程约束的基础上，建立时间最优的星座组网电推进转移轨道优化问题；
模块M2：对组网过程进行阶段划分，并将推力大小及方向进行分阶段参数化处理；
模块M3：确定待求解变量，将星座组网电推进转移轨道优化问题转换为控制变量受约束的打靶问题；
模块M4：对打靶问题的初值进行猜测，采用牛顿迭代算法进行求解，最终得到星座组网电推进转移轨道优化问题的最优解。
优选的，所述模块M1包括：
模块M1.1：在地心惯性参考坐标系下，建立卫星的轨道动力学方程：
Figure BDA0003756199530000062
Figure BDA0003756199530000063
Figure BDA0003756199530000064
其中，r和v分表为位置和速度矢量；fP为包含J2、大气阻力的摄动加速度；μ为地球引力常数；r为地球半径；ge为地球海平面重力加速度；m为航天器质量；Isp为推力器比冲；α为推力方向的单位矢量；Tmax为最大推力；u为实际推力与最大推力的比值；
在仅有开和关两种状态的工程约束的基础上，u的取值为0或者1；
推力方向的单位矢量在惯性坐标系下进一步表示为：
Figure BDA0003756199530000066
其中，
Figure BDA0003756199530000065
为推力方位角，δ为推力高度角；
令卫星的瞬时轨道根数为：
kosc＝[a e i Ω ω M]
其中，a为半长轴，e为偏心率，i为轨道倾角，Ω为升交点赤经，ω为近地点幅角，M为平近点角；
令
Figure BDA0003756199530000071
表示卫星的平均轨道根数，x＝[r v]为直角坐标参数，则kosc、kave及x三者可互相转换；
模块M1.2：令初始时刻为t0，末端时刻为tf，则初末时刻卫星面内的组网约束表示为：
r(t0)＝r0，v(t0)＝v0，m(t0)＝m0
Figure BDA0003756199530000072
Figure BDA0003756199530000073
式中，
Figure BDA0003756199530000074
及
Figure BDA0003756199530000075
为常值，表示末端时刻期望的半长轴、偏心率及近地点幅角，
Figure BDA0003756199530000076
表示沿迹角，定义为：
Figure BDA0003756199530000077
在末端时刻，期望的
Figure BDA0003756199530000078
根据如下公式确定：
Figure BDA0003756199530000079
式中，
Figure BDA00037561995300000710
表示t0时刻的初始，
Figure BDA00037561995300000711
表示沿迹角变化率；
模块M1.3：令控制变量和状态变量：
u＝[u α]T
y＝[r v m]T
则轨道动力学方程表示为：
Figure BDA00037561995300000712
相应的边界约束表示为：
bini[y(t0)]＝y(t0)-y0＝0
Figure BDA0003756199530000081
其中，y0＝[r0 v0 m0]T
组网过程电推进转移轨道问题P为寻找最优控制变量：
Figure BDA0003756199530000082
同时满足：
Figure BDA0003756199530000083
bini[y(t0)]＝0，bter[y(tf)]＝0
u＝0或者u＝1，αTα＝1。
优选的，所述模块M2包括：
模块M2.1：根据最快入轨的性能指标要求和组网轨道的运动特性，将u和α在时间域上按时刻t1、t2及t3进行划分，在推力为全开的基础上，将转移轨道按时间顺序分为巡航阶段、升轨阶段及偏心率矢量调整阶段；
模块M2.2：将各阶段的推力大小和方向设定为：
Figure BDA0003756199530000084
Figure BDA0003756199530000085
式中，αa和αeω分别为升轨阶段和偏心率矢量调整阶段的推力方向单位矢量；
模块M2.3：为满足最快入轨的性能指标要求，在升轨阶段和偏心率矢量调整阶段，将推力方向分别设定沿着切向方向和轨道面内垂直拱线方向，表示为：
Figure BDA0003756199530000086
αeω＝(cosσ)ip+(sinσ)jp
其中，σ表示推力方向与拱线的夹角，且：
Figure BDA0003756199530000091
Figure BDA0003756199530000092
其中，ip和jp的计算是基于平均轨道参数进行的，确保推力方向平缓变化。
优选的，所述模块M3包括：
模块M3.1：对于给定的初始状态t0和y0，定义待求解的变量为：
z＝[Δt1 Δt2 Δt3 σ]T
其中，Δt1＝t1-t0，Δt2＝t2-t1，Δt3＝t3-t2；
模块M3.2：将问题P转换为如下打靶问题：
Figure BDA0003756199530000093
优选的，所述模块M4包括：
模块M4.1：对Δt1、Δt2、Δt3及σ进行猜测，确定迭代变量初值；
模块M4.2：采用龙格库塔法计算代价函数F(z)，采用牛顿迭代算法进行迭代计算，直至||F(z)||收敛至小于指定阈值。
与现有技术相比，本发明具有如下的有益效果：
本发明中的方法考虑电推进推力大小固定，仅有开和关两种状态，根据最快入轨要求，充分利于组网轨道的运动特性，提出一种推力参数化方法，可以保证转移轨道满足预先设定的组网约束，对姿控系统要求低，且推力参数仅为4个，结构简单，求解效率高，有利于星上实现，可以有效实现轨道根数的联合控制。
附图说明
通过阅读参照以下附图对非限制性实施例所作的详细描述，本发明的其它特征、目的和优点将会变得更明显：
图1为本发明方法流程图；
图2为具体实施方式中推力剖面；
图3为具体实施方式中半长轴及沿迹角偏差；
图4为具体实施方式中偏心率矢量及推力大小；
图5为具体实施方式中惯性坐标系下推力方向角及燃料消耗。
具体实施方式
下面结合具体实施例对本发明进行详细说明。以下实施例将有助于本领域的技术人员进一步理解本发明，但不以任何形式限制本发明。应当指出的是，对本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变化和改进。这些都属于本发明的保护范围。
实施例：
本发明提供一种低轨卫星星座组网电推进变轨策略优化算法，考虑电推进推力大小固定，仅有开和关两种状态，根据最快入轨要求，充分利于组网轨道的运动特性，提出一种推力参数化方法，可以保证转移轨道满足预先设定的组网约束，对姿控系统要求低，且推力参数仅为4个，结构简单，求解效率高，有利于星上实现，可以有效实现轨道根数的联合控制。
如图1，本发明提供了一种低轨卫星星座组网电推进变轨策略优化算法，包括如下步骤：
步骤1：考虑电推进推力大小固定，仅有开和关两种状态的工程约束，建立时间最优的星座组网电推进转移轨道优化问题。
其中，具体过程包括：
步骤101：在地心惯性参考坐标系下，建立卫星的轨道动力学方程为：
Figure BDA0003756199530000101
Figure BDA0003756199530000102
其中，r和v分表为位置和速度矢量，fP为包含J2、大气阻力等摄动加速度，μ为地球引力常数，r为地球半径，ge为地球海平面重力加速度，m为航天器质量，Isp为推力器比冲，α为推力方向的单位矢量，Tmax为最大推力，u为实际推力与最大推力的比值。本文考虑发动机推力不可连续调整，仅有开和关两种状态的情况，因此u的取值为0或者1。推力方向的单位矢量在惯性坐标系下可进一步表示为：
Figure BDA0003756199530000103
其中，
Figure BDA0003756199530000111
为推力方位角，δ为推力高度角。
令卫星的瞬时轨道根数为：
kosc＝[a e i Ω ω M]
其中，a为半长轴，e为偏心率，i为轨道倾角，Ω为升交点赤经，ω为近地点幅角，M为平近点角。令
Figure BDA0003756199530000112
表示卫星的平均轨道根数，x＝[r v]为直角坐标参数，则kosc、kave及x三者可互相转换。
步骤102：令初始时刻为t0，末端时刻为tf，则初末时刻卫星面内的组网约束可以表示为：
r(t0)＝r0，v(t0)＝v0，m(t0)＝m0
Figure BDA00037561995300001113
Figure BDA0003756199530000113
式中，
Figure BDA0003756199530000114
及
Figure BDA0003756199530000115
为常值，表示末端时刻期望的半长轴、偏心率及近地点幅角，
Figure BDA0003756199530000116
表示沿迹角，定义为：
Figure BDA0003756199530000117
在末端时刻，期望的
Figure BDA0003756199530000118
可根据如下式子确定：
Figure BDA0003756199530000119
式中，
Figure BDA00037561995300001110
表示t0时刻的初始，
Figure BDA00037561995300001111
表示沿迹角变化率。
步骤103：为便于阐述，令控制变量和状态变量：
u＝[u α]T
y＝[r v m]T
则轨道动力学方程可以表示为：
Figure BDA00037561995300001112
相应的边界约束表示为：
bini[y(t0)]＝y(t0)-y0＝0
Figure BDA0003756199530000121
其中，y0＝[r0 v0 m0]T。
组网过程电推进转移轨道问题P可以阐述如下：寻找最优控制变量；
Figure BDA0003756199530000122
同时满足：
Figure BDA0003756199530000123
bini[y(t0)]＝0，bter[y(tf)]＝0
u＝0或者1，αTα＝1
步骤2：对组网过程进行阶段划分，并将推力大小及方向分阶段参数化处理。
其中，步骤2具体包括：
步骤201：根据最快入轨的性能指标要求，充分考虑组网轨道的运动特性，将u和α在时间域上按时刻t1、t2及t3进行划分，在推力基本为全开的基础上，将转移轨道按时间顺序分为巡航阶段、升轨阶段及偏心率矢量调整阶段，如图2。
步骤202：将各阶段的推力大小和方向设定为：
Figure BDA0003756199530000124
Figure BDA0003756199530000125
式中，αa和αeω分别为升轨阶段和偏心率矢量调整阶段的推力方向单位矢量。
步骤203：为满足最快入轨的性能指标要求，在升轨阶段和偏心率矢量调整阶段，将推力方向分别设定沿着切向方向和轨道面内垂直拱线方向，具体表示为：
Figure BDA0003756199530000126
αeω＝(cosσ)ip+(sinσ)jp
其中，σ表示推力方向与拱线的夹角，且：
Figure BDA0003756199530000131
Figure BDA0003756199530000132
需要指出的是，为确保推力方向平缓变化，ip和jp的计算是基于平均轨道参数进行的。
步骤3：确定待求解变量，建立一个控制变量受约束的打靶问题。
该步骤的具体过程为：
步骤301：经过步骤2的转换，对于给定的初始状态t0和y0，定义待求解的变量为
z＝[Δt1 Δt2 Δt3 σ]T
其中，Δt1＝t1-t0，Δt2＝t2-t1，Δt3＝t3-t2。
步骤302：将问题P转换为如下打靶问题：
Figure BDA0003756199530000133
可以看出，上述方程组中，未知变量的个数与约束方程数量相等，因此方程的解是完全确定的。
步骤4：对打靶问题的初值进行猜测，采用牛顿迭代算法进行求解。
该步骤的具体过程为：
步骤401：对Δt1、Δt2、Δt3及σ进行猜测，确定迭代变量初值；
步骤402：采用龙格-库塔数值计算算法计算代价函数F(z)，采用牛顿迭代算法进行迭代计算，直至||F(z)||收敛至小于指定阈值。
以下为一种低轨卫星星座组网电推进变轨策略优化算法的数值仿真验证。
假设卫星采用偏低发射入轨，星间分离后卫星的初始轨道高度为1000km，目标轨道高度为1200km。假设卫星的初始质量为170kg，推力大小为20mN，比冲为1600s。在表中，将目标轨道的近地点幅角和偏心率选择为90deg和0.001是为了实现目标轨道的冻结。考虑目标轨道的倾角接近90deg，其轨道冻结受地球高阶引力场影响大，为此在仿真中，将考虑地球21阶次的引力场模型。由于初始轨道偏心率接近0，近地点幅角无法定义，为此在后续分析中采用偏心率矢量的两个分量来描述，即：
ex＝ecosω，ey＝esinω
在目标轨道，ex和ey的目标值为0和0.001。在仿真计算结果和分析中，若无特别说明，则所有的轨道根数均是指平均轨道根数。
根据以上计算条件，进行数值仿真计算，结果如图3、图4、图5所示。
根据本发明提供的低轨卫星星座组网电推进变轨策略优化系统，包括：模块M1：根据电推进的推力大小，在仅有开和关两种状态的工程约束的基础上，建立时间最优的星座组网电推进转移轨道优化问题；模块M2：对组网过程进行阶段划分，并将推力大小及方向进行分阶段参数化处理；模块M3：确定待求解变量，将星座组网电推进转移轨道优化问题转换为控制变量受约束的打靶问题；模块M4：对打靶问题的初值进行猜测，采用牛顿迭代算法进行求解，最终得到星座组网电推进转移轨道优化问题的最优解。
所述模块M1包括：模块M1.1：在地心惯性参考坐标系下，建立卫星的轨道动力学方程：
Figure BDA0003756199530000141
Figure BDA0003756199530000142
Figure BDA0003756199530000143
其中，r和v分表为位置和速度矢量；fP为包含J2、大气阻力的摄动加速度；μ为地球引力常数；r为地球半径；ge为地球海平面重力加速度；m为航天器质量；Isp为推力器比冲；α为推力方向的单位矢量；Tmax为最大推力；u为实际推力与最大推力的比值；在仅有开和关两种状态的工程约束的基础上，u的取值为0或者1；推力方向的单位矢量在惯性坐标系下进一步表示为：
Figure BDA0003756199530000144
其中，
Figure BDA0003756199530000145
为推力方位角，δ为推力高度角；令卫星的瞬时轨道根数为：kosc＝[a e i Ω ω M]，其中，a为半长轴，e为偏心率，i为轨道倾角，Ω为升交点赤经，ω为近地点幅角，M为平近点角；令
Figure BDA0003756199530000146
表示卫星的平均轨道根数，x＝[r v]为直角坐标参数，则kosc、kave及x三者可互相转换；
模块M1.2：令初始时刻为t0，末端时刻为tf，则初末时刻卫星面内的组网约束表示为：r(t0)＝r0，v(t0)＝v0，m(t0)＝m0
Figure BDA0003756199530000147
Figure BDA0003756199530000151
式中，
Figure BDA0003756199530000152
及
Figure BDA0003756199530000153
为常值，表示末端时刻期望的半长轴、偏心率及近地点幅角，
Figure BDA0003756199530000154
表示沿迹角，定义为：
Figure BDA0003756199530000155
在末端时刻，期望的
Figure BDA0003756199530000156
根据如下公式确定：
Figure BDA0003756199530000157
式中，
Figure BDA0003756199530000158
表示t0时刻的初始，
Figure BDA0003756199530000159
表示沿迹角变化率；
模块M1.3：令控制变量和状态变量：u＝[u α]T，y＝[r v m]T
则轨道动力学方程表示为：
Figure BDA00037561995300001510
相应的边界约束表示为：bini[y(t0)]＝y(t0)-y0＝0
Figure BDA00037561995300001511
其中，y0＝[r0 v0 m0]T
组网过程电推进转移轨道问题P为寻找最优控制变量：
Figure BDA00037561995300001512
同时满足：
Figure BDA00037561995300001513
bini[y(t0)]＝0，bter[y(tf)]＝0，u＝0或者u＝1，αTα＝1。
所述模块M2包括：模块M2.1：根据最快入轨的性能指标要求和组网轨道的运动特性，将u和α在时间域上按时刻t1、t2及t3进行划分，在推力为全开的基础上，将转移轨道按时间顺序分为巡航阶段、升轨阶段及偏心率矢量调整阶段；模块M2.2：将各阶段的推力大小和方向设定为：
Figure BDA00037561995300001514
式中，αa和αeω分别为升轨阶段和偏心率矢量调整阶段的推力方向单位矢量；模块M2.3：为满足最快入轨的性能指标要求，在升轨阶段和偏心率矢量调整阶段，将推力方向分别设定沿着切向方向和轨道面内垂直拱线方向，表示为：
Figure BDA00037561995300001515
αeω＝(cosσ)ip+(sinσ)jp，其中，σ表示推力方向与拱线的夹角，且：
Figure BDA0003756199530000161
其中，ip和jp的计算是基于平均轨道参数进行的，确保推力方向平缓变化。
所述模块M3包括：模块M3.1：对于给定的初始状态t0和y0，定义待求解的变量为：z＝[Δt1 Δt2 Δt3 σ]T，其中，Δt1＝t1-t0，Δt2＝t2-t1，Δt3＝t3-t2；
模块M3.2：将问题P转换为如下打靶问题：
Figure BDA0003756199530000162
所述模块M4包括：模块M4.1：对Δt1、Δt2、Δt3及σ进行猜测，确定迭代变量初值；模块M4.2：采用龙格库塔法计算代价函数F(z)，采用牛顿迭代算法进行迭代计算，直至||F(z)||收敛至小于指定阈值。
本领域技术人员知道，除了以纯计算机可读程序代码方式实现本发明提供的系统、装置及其各个模块以外，完全可以通过将方法步骤进行逻辑编程来使得本发明提供的系统、装置及其各个模块以逻辑门、开关、专用集成电路、可编程逻辑控制器以及嵌入式微控制器等的形式来实现相同程序。所以，本发明提供的系统、装置及其各个模块可以被认为是一种硬件部件，而对其内包括的用于实现各种程序的模块也可以视为硬件部件内的结构；也可以将用于实现各种功能的模块视为既可以是实现方法的软件程序又可以是硬件部件内的结构。
以上对本发明的具体实施例进行了描述。需要理解的是，本发明并不局限于上述特定实施方式，本领域技术人员可以在权利要求的范围内做出各种变化或修改，这并不影响本发明的实质内容。在不冲突的情况下，本申请的实施例和实施例中的特征可以任意相互组合。

Claims (10)

1.一种低轨卫星星座组网电推进变轨策略优化方法，其特征在于，包括：
步骤1：根据电推进的推力大小，在仅有开和关两种状态的工程约束的基础上，建立时间最优的星座组网电推进转移轨道优化问题；
步骤2：对组网过程进行阶段划分，并将推力大小及方向进行分阶段参数化处理；
步骤3：确定待求解变量，将星座组网电推进转移轨道优化问题转换为控制变量受约束的打靶问题；
步骤4：对打靶问题的初值进行猜测，采用牛顿迭代算法进行求解，最终得到星座组网电推进转移轨道优化问题的最优解。
2.根据权利要求1所述的低轨卫星星座组网电推进变轨策略优化方法，其特征在于，所述步骤1包括：
步骤1.1：在地心惯性参考坐标系下，建立卫星的轨道动力学方程：
Figure FDA0003756199520000011
Figure FDA0003756199520000012
Figure FDA0003756199520000013
其中，r和v分表为位置和速度矢量；fP为包含J2、大气阻力的摄动加速度；μ为地球引力常数；r为地球半径；ge为地球海平面重力加速度；m为航天器质量；Isp为推力器比冲；α为推力方向的单位矢量；Tmax为最大推力；u为实际推力与最大推力的比值；
在仅有开和关两种状态的工程约束的基础上，u的取值为0或者1；
推力方向的单位矢量在惯性坐标系下进一步表示为：
Figure FDA0003756199520000014
其中，
Figure FDA0003756199520000015
为推力方位角，δ为推力高度角；
令卫星的瞬时轨道根数为：
kosc＝[a e i Ω ω M]
其中，a为半长轴，e为偏心率，i为轨道倾角，Ω为升交点赤经，ω为近地点幅角，M为平近点角；
令
Figure FDA0003756199520000016
表示卫星的平均轨道根数，x＝[r v]为直角坐标参数，则kosc、kave及x三者可互相转换；
步骤1.2：令初始时刻为t0，末端时刻为tf，则初末时刻卫星面内的组网约束表示为：
r(t0)＝r0，v(t0)＝v0，m(t0)＝m0
Figure FDA0003756199520000021
Figure FDA0003756199520000022
式中，
Figure FDA0003756199520000023
及
Figure FDA0003756199520000024
为常值，表示末端时刻期望的半长轴、偏心率及近地点幅角，
Figure FDA0003756199520000025
表示沿迹角，定义为：
Figure FDA0003756199520000026
在末端时刻，期望的
Figure FDA0003756199520000027
根据如下公式确定：
Figure FDA0003756199520000028
式中，
Figure FDA0003756199520000029
表示t0时刻的初始，
Figure FDA00037561995200000210
表示沿迹角变化率；
步骤1.3：令控制变量和状态变量：
u＝[u α]T
y＝[r v m]T
则轨道动力学方程表示为：
Figure FDA00037561995200000211
相应的边界约束表示为：
bini[y(t0)]＝y(t0)-y0＝0
Figure FDA00037561995200000212
其中，y0＝[r0 v0 m0]T
组网过程电推进转移轨道问题P为寻找最优控制变量：
Figure FDA00037561995200000213
同时满足：
Figure FDA0003756199520000031
bini[y(t0)]＝0，bter[y(tf)]＝0
u＝0或者u＝1，αTα＝1。
3.根据权利要求2所述的低轨卫星星座组网电推进变轨策略优化方法，其特征在于，所述步骤2包括：
步骤2.1：根据最快入轨的性能指标要求和组网轨道的运动特性，将u和α在时间域上按时刻t1、t2及t3进行划分，在推力为全开的基础上，将转移轨道按时间顺序分为巡航阶段、升轨阶段及偏心率矢量调整阶段；
步骤2.2：将各阶段的推力大小和方向设定为：
Figure FDA0003756199520000032
Figure FDA0003756199520000033
式中，αa和αeω分别为升轨阶段和偏心率矢量调整阶段的推力方向单位矢量；
步骤2.3：为满足最快入轨的性能指标要求，在升轨阶段和偏心率矢量调整阶段，将推力方向分别设定沿着切向方向和轨道面内垂直拱线方向，表示为：
Figure FDA0003756199520000034
αeω＝(cosσ)ip+(sinσ)jp
其中，σ表示推力方向与拱线的夹角，且：
Figure FDA0003756199520000035
Figure FDA0003756199520000036
其中，ip和jp的计算是基于平均轨道参数进行的，确保推力方向平缓变化。
4.根据权利要求3所述的低轨卫星星座组网电推进变轨策略优化方法，其特征在于，所述步骤3包括：
步骤3.1：对于给定的初始状态t0和y0，定义待求解的变量为：
z＝[Δt1 Δt2 Δt3 σ]T
其中，Δt1＝t1-t0，Δt2＝t2-t1，Δt3＝t3-t2；
步骤3.2：将问题P转换为如下打靶问题：
Figure FDA0003756199520000041
5.根据权利要求4所述的低轨卫星星座组网电推进变轨策略优化方法，其特征在于，所述步骤4包括：
步骤4.1：对Δt1、Δt2、Δt3及σ进行猜测，确定迭代变量初值；
步骤4.2：采用龙格库塔法计算代价函数F(z)，采用牛顿迭代算法进行迭代计算，直至||F(z)||收敛至小于指定阈值。
6.一种低轨卫星星座组网电推进变轨策略优化系统，其特征在于，包括：
模块M1：根据电推进的推力大小，在仅有开和关两种状态的工程约束的基础上，建立时间最优的星座组网电推进转移轨道优化问题；
模块M2：对组网过程进行阶段划分，并将推力大小及方向进行分阶段参数化处理；
模块M3：确定待求解变量，将星座组网电推进转移轨道优化问题转换为控制变量受约束的打靶问题；
模块M4：对打靶问题的初值进行猜测，采用牛顿迭代算法进行求解，最终得到星座组网电推进转移轨道优化问题的最优解。
7.根据权利要求6所述的低轨卫星星座组网电推进变轨策略优化系统，其特征在于，所述模块M1包括：
模块M1.1：在地心惯性参考坐标系下，建立卫星的轨道动力学方程：
Figure FDA0003756199520000042
Figure FDA0003756199520000043
Figure FDA0003756199520000044
其中，r和v分表为位置和速度矢量；fP为包含J2、大气阻力的摄动加速度；μ为地球引力常数；r为地球半径；ge为地球海平面重力加速度；m为航天器质量；Isp为推力器比冲；α为推力方向的单位矢量；Tmax为最大推力；u为实际推力与最大推力的比值；
在仅有开和关两种状态的工程约束的基础上，u的取值为0或者1；
推力方向的单位矢量在惯性坐标系下进一步表示为：
Figure FDA0003756199520000051
其中，
Figure FDA0003756199520000052
为推力方位角，δ为推力高度角；
令卫星的瞬时轨道根数为：
kosc＝[a e i Ω ω M]
其中，a为半长轴，e为偏心率，i为轨道倾角，Ω为升交点赤经，ω为近地点幅角，M为平近点角；
令
Figure FDA0003756199520000053
表示卫星的平均轨道根数，x＝[r v]为直角坐标参数，则kosc、kave及x三者可互相转换；
模块M1.2：令初始时刻为t0，末端时刻为tf，则初末时刻卫星面内的组网约束表示为：
r(t0)＝r0，v(t0)＝v0，m(t0)＝m0
Figure FDA0003756199520000054
Figure FDA0003756199520000055
式中，
Figure FDA0003756199520000056
及
Figure FDA0003756199520000057
为常值，表示末端时刻期望的半长轴、偏心率及近地点幅角，
Figure FDA0003756199520000058
表示沿迹角，定义为：
Figure FDA0003756199520000059
在末端时刻，期望的
Figure FDA00037561995200000510
根据如下公式确定：
Figure FDA00037561995200000511
式中，
Figure FDA00037561995200000512
表示t0时刻的初始，
Figure FDA00037561995200000513
表示沿迹角变化率；
模块M1.3：令控制变量和状态变量：
u＝[u α]T
y＝[r v m]T
则轨道动力学方程表示为：
Figure FDA0003756199520000061
相应的边界约束表示为：
bini[y(t0)]＝y(t0)-y0＝0
Figure FDA0003756199520000062
其中，y0＝[r0 v0 m0]T
组网过程电推进转移轨道问题P为寻找最优控制变量：
Figure FDA0003756199520000063
同时满足：
Figure FDA0003756199520000064
bini[y(t0)]＝0，bter[y(tf)]＝0
u＝0或者u＝1，αTα＝1。
8.根据权利要求7所述的低轨卫星星座组网电推进变轨策略优化系统，其特征在于，所述模块M2包括：
模块M2.1：根据最快入轨的性能指标要求和组网轨道的运动特性，将u和α在时间域上按时刻t1、t2及t3进行划分，在推力为全开的基础上，将转移轨道按时间顺序分为巡航阶段、升轨阶段及偏心率矢量调整阶段；
模块M2.2：将各阶段的推力大小和方向设定为：
Figure FDA0003756199520000065
Figure FDA0003756199520000066
式中，αa和αeω分别为升轨阶段和偏心率矢量调整阶段的推力方向单位矢量；
模块M2.3：为满足最快入轨的性能指标要求，在升轨阶段和偏心率矢量调整阶段，将推力方向分别设定沿着切向方向和轨道面内垂直拱线方向，表示为：
Figure FDA0003756199520000071
αeω＝(cosσ)ip+(sinσ)jp
其中，σ表示推力方向与拱线的夹角，且：
Figure FDA0003756199520000072
Figure FDA0003756199520000073
其中，ip和jp的计算是基于平均轨道参数进行的，确保推力方向平缓变化。
9.根据权利要求8所述的低轨卫星星座组网电推进变轨策略优化系统，其特征在于，所述模块M3包括：
模块M3.1：对于给定的初始状态t0和y0，定义待求解的变量为：
z＝[Δt1 Δt2 Δt3 σ]T
其中，Δt1＝t1-t0，Δt2＝t2-t1，Δt3＝t3-t2；
模块M3.2：将问题P转换为如下打靶问题：
Figure FDA0003756199520000074
10.根据权利要求9所述的低轨卫星星座组网电推进变轨策略优化系统，其特征在于，所述模块M4包括：
模块M4.1：对Δt1、Δt2、Δt3及σ进行猜测，确定迭代变量初值；
模块M4.2：采用龙格库塔法计算代价函数F(z)，采用牛顿迭代算法进行迭代计算，直至||F(z)||收敛至小于指定阈值。
CN202210857788.2A 2022-07-20 2022-07-20 低轨卫星星座组网电推进变轨策略优化方法和系统 Pending CN115258196A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210857788.2A CN115258196A (zh) 2022-07-20 2022-07-20 低轨卫星星座组网电推进变轨策略优化方法和系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210857788.2A CN115258196A (zh) 2022-07-20 2022-07-20 低轨卫星星座组网电推进变轨策略优化方法和系统

Publications (1)

Publication Number Publication Date
CN115258196A true CN115258196A (zh) 2022-11-01

Family

ID=83767093

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210857788.2A Pending CN115258196A (zh) 2022-07-20 2022-07-20 低轨卫星星座组网电推进变轨策略优化方法和系统

Country Status (1)

Country Link
CN (1) CN115258196A (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116520711A (zh) * 2023-07-03 2023-08-01 中国西安卫星测控中心 一种电推卫星walker星座组网控制筹划方法
CN116692034A (zh) * 2023-08-07 2023-09-05 北京航天驭星科技有限公司 时间最优相位差调整方法、系统、电子设备和介质

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116520711A (zh) * 2023-07-03 2023-08-01 中国西安卫星测控中心 一种电推卫星walker星座组网控制筹划方法
CN116520711B (zh) * 2023-07-03 2023-10-13 中国西安卫星测控中心 一种电推卫星walker星座组网控制筹划方法
CN116692034A (zh) * 2023-08-07 2023-09-05 北京航天驭星科技有限公司 时间最优相位差调整方法、系统、电子设备和介质
CN116692034B (zh) * 2023-08-07 2023-09-29 北京航天驭星科技有限公司 时间最优相位差调整方法、系统、电子设备和介质

Similar Documents

Publication Publication Date Title
CN115258196A (zh) 低轨卫星星座组网电推进变轨策略优化方法和系统
Baig et al. Light-levitated geostationary cylindrical orbits are feasible
CN114460952B (zh) 椭圆轨道伴飞构型初始化的双星协同变轨方法及系统
Zhang et al. Autonomous guidance for rendezvous phasing based on special-point-based maneuvers
CN115755598A (zh) 一种智能航天器集群分布式模型预测路径规划方法
Benedikter et al. Convex optimization of launch vehicle ascent trajectory with heat-flux and splash-down constraints
Zhao et al. A novel two-level optimization strategy for multi-debris active removal mission in LEO
Kluever Optimal Earth-Moon trajectories using combined chemical-electric propulsion
Edlerman et al. Analysis of nanosatellite formation establishment and maintenance using electric propulsion
Smet et al. Systematic exploration of solar gravity driven orbital transfers in the Martian system using artificial neural networks
Zhu et al. The Intelligent Trajectory Optimization of Multistage Rocket with Gauss Pseudo-Spectral Method.
Ortega Fuzzy logic techniques for rendezvous and docking of two geostationary satellites
. Brusch Constrained impulsive trajectory optimization for orbit-to-orbit transfer
Schäff et al. Advanced Electric Orbit-Raising Optimization for Operational Purpose
Chernenkii Study of gas giant satellites system
Shirobokov et al. Parametric analysis of low-thrust lunar transfers with resonant encounters
Pagano et al. Global launcher trajectory optimization for lunar base settlement
Ozimek Low-thrust trajectory design and optimization of lunar south pole coverage missions
Sin Trajectory optimization and control of small spacecraft constellations
Vaz et al. Sub-optimal orbital maneuvers for artificial satellites
CN113779788B (zh) 一种大入轨偏差情况下组合体分离轨道确定方法和系统
Marrdonny et al. A guidance algorithm for launch to equatorial orbit
Pergola Semianalytic Approach for Optimal Configuration of Electric Propulsion Spacecraft
Jesick et al. Automated design of optimal finite thrust orbit insertion with ballistic flyby constraints
Sosa et al. Optimized Launch to Rendezvous Trajectories for a Mars Ascent Vehicle Using Direct Collocation

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