CN111874267B - 基于粒子群算法的低轨卫星离轨控制方法及系统 - Google Patents

基于粒子群算法的低轨卫星离轨控制方法及系统 Download PDF

Info

Publication number
CN111874267B
CN111874267B CN202010629155.7A CN202010629155A CN111874267B CN 111874267 B CN111874267 B CN 111874267B CN 202010629155 A CN202010629155 A CN 202010629155A CN 111874267 B CN111874267 B CN 111874267B
Authority
CN
China
Prior art keywords
target
particle
fitness
orbit
population
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
CN202010629155.7A
Other languages
English (en)
Other versions
CN111874267A (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.)
Peoples Liberation Army Strategic Support Force Aerospace Engineering University
Original Assignee
Peoples Liberation Army Strategic Support Force Aerospace Engineering 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 Peoples Liberation Army Strategic Support Force Aerospace Engineering University filed Critical Peoples Liberation Army Strategic Support Force Aerospace Engineering University
Publication of CN111874267A publication Critical patent/CN111874267A/zh
Priority to US17/363,363 priority Critical patent/US11787567B2/en
Application granted granted Critical
Publication of CN111874267B publication Critical patent/CN111874267B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

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
    • B64G1/242—Orbits and trajectories
    • 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/26—Guiding or controlling apparatus, e.g. for attitude control using jets
    • 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/62—Systems for re-entry into the earth's atmosphere; Retarding or landing devices
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/004—Artificial life, i.e. computing arrangements simulating life
    • G06N3/006—Artificial life, i.e. computing arrangements simulating life based on simulated virtual individual or collective life forms, e.g. social simulations or particle swarm optimisation [PSO]
    • 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/247—Advanced control concepts for autonomous, robotic spacecraft, e.g. by using artificial intelligence, neural networks or autonomous agents
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/02—Neural networks
    • G06N3/08—Learning methods

Abstract

本申请公开了一种基于粒子群算法的低轨卫星离轨控制方法及系统，属于星座管理领域，能够解决在离轨控制方向缺乏算法实现的研究等问题。所述方法包括：根据粒子群中粒子的位置参数和摄动加速度，以及预设的拉格朗日乘子和惩罚因子，构建目标函数；基于目标函数计算粒子群中各粒子的目标函数值，根据计算结果确定每个粒子的目标适应度和粒子群的种群适应度，并更新各粒子的位置参数和速度，以及获取最大迭代次数时的种群最优适应度；根据种群最优适应度更新拉格朗日乘子和惩罚因子，并将目标偏差值与预设收敛条件比较，根据比较结果确定目标种群适应度和目标位置参数，得到目标离轨方式。本申请用于低轨卫星的离轨控制。

Description

基于粒子群算法的低轨卫星离轨控制方法及系统
技术领域
本申请涉及星座运行管理技术领域，尤其涉及一种基于粒子群算法的低轨卫星离轨控制方法及系统。
背景技术
随着人类航天发射任务的不断增加，低轨及中高轨轨道空间逐渐拥挤，使在轨运行的卫星碰撞概率逐渐增大。根据机构空间碎片协调委员会(IADC)编订的《IADC空间碎片减缓指南》，在自然轨道超过25年的情况下，低轨到寿卫星需要在25年内进入处置轨道进而在大气层内烧毁。同时，小卫星技术的发展和发射成本的降低使得低轨卫星星座中卫星数量剧增，OneWeb公司计划发射由2620颗卫星组成的星座、Samgsung公司计划发射由4600颗卫星组成的星座、Boeing计划发射2956颗卫星组成的星座和Starlink最终计划发射4.2万颗卫星组成星座，这些星座规模巨大，其卫星到寿后将会影响低轨的空间环境。因此，采取合理方式对低轨卫星星座中到寿卫星进行离轨控制变得尤为重要。
相关技术中，针对低地球轨道(low earth orbit，LEO)卫星，离轨方法一般分为电动力电缆离轨、大气阻力离轨和卫星自身推力器离轨，卫星自身推力器离轨相较于其他离轨方式具有离轨时间短、控制灵活等特点。电推力器作为一种典型的小推力器已经在深空探测、轨道维持和轨道机动中得到广泛的应用，Starlink星座中所有卫星均带有电推力器。Fromm等人分别仿真不同小推力器作用下卫星的离轨时间，并得出了适于仿真对象卫星离轨的小推力器推力范围，Huang等人基于两种不同的小推力离轨策略对单星和OneWeb星座离轨过程进行了仿真分析，并得出最优的单星及星座离轨策略。然而，相关技术中离轨控制研究主要集中在控制率推导方面，对算法实现少有提及，粒子群算法在各类连续空间优化问题、神经网络训练等领域中取得良好效果，但在离轨控制方面少有提及，一定程度上限制了卫星空间运动的进一步研究和准确控制。
发明内容
本申请目的在于提供一种基于粒子群算法的低轨卫星离轨控制方法及系统，进而至少在一定程度上克服了在离轨控制方向缺乏算法实现的研究等问题。
为实现上述目的，本申请采用的技术方案如下：
根据本申请的一个方面，提供一种基于粒子群算法的低轨卫星离轨控制方法，所述方法包括：
步骤S1、根据粒子群中粒子的位置参数和摄动加速度，以及预设的拉格朗日乘子和惩罚因子，构建目标函数；
步骤S2、基于所述目标函数，计算所述粒子群中各粒子的目标函数值，根据计算结果确定每个粒子的目标适应度和所述粒子群的种群适应度，并更新各粒子的位置参数和速度，以及获取最大迭代次数时的种群最优适应度和其对应的粒子的位置参数；
步骤S3、根据所述种群最优适应度更新所述拉格朗日乘子和惩罚因子，并获取当前迭代次数和基于所述种群最优适应度的目标偏差值，将所述当前迭代次数和所述目标偏差值与预设迭代条件和预设收敛条件进行比较，根据比较结果确定目标种群适应度和目标位置参数，得到目标离轨方式。
在本申请的一种示例性实施例中，所述根据粒子群中粒子的位置参数和摄动加速度，以及预设的拉格朗日乘子和惩罚因子，构建目标函数，包括：
步骤S11、获取在推力器推力的作用下，各轨道根数随时间的变化率，并根据所述变化率确定含有当前协状态变量的控制率；
步骤S12、计算所述推力器推力的推力矢量、燃料消耗率及摄动加速度；
步骤S13、根据所述当前协状态变量确定粒子位置向量，并将所述粒子位置向量作为粒子群中粒子的位置参数；
步骤S14、基于所述位置参数和摄动加速度，将所述控制率进行轨道积分，并根据积分结果、预设的拉格朗日乘子和惩罚因子，构建所述目标函数。
在本申请的一种示例性实施例中，所述将所述当前迭代次数和所述目标偏差值与预设迭代条件和预设收敛条件进行比较，根据比较结果确定目标种群适应度和目标位置参数，包括：
若所述当前迭代次数满足预设迭代条件且所述目标偏差值满足收敛条件，则将所述种群最优适应度确定为目标种群适应度，将所述目标种群适应度对应的位置参数确定为目标位置参数；
若所述当前迭代次数不满足所述预设迭代条件或所述目标偏差值不满足所述收敛条件，则基于更新后的拉格朗日乘子和惩罚因子更新目标函数，并继续执行步骤S14及后续步骤，直至所述当前迭代次数满足所述预设迭代条件，且所述目标偏差值满足所述收敛条件。
在本申请的一种示例性实施例中，步骤S11中，所述根据所述变化率确定含有当前协状态变量的控制率，包括：
基于所述变化率和所述当前协状态变量，通过公式(7)得到哈密尔顿函数：
Figure GDA0003164856570000031
其中，H为哈密尔顿函数值，λa、λe和λi为当前协状态变量，a、e和i分别为工作轨道的半长轴、偏心率和轨道倾角；
根据所述哈密尔顿函数对所述推力器推力对应的预设角度求导，得到所述含有当前协状态变量的控制率。
在本申请的一种示例性实施例中，步骤S12中，通过公式(11)得到推力矢量，通过公式(12)得到燃料消耗率，通过公式(14)得到摄动加速度；
Figure GDA0003164856570000032
Figure GDA0003164856570000033
Figure GDA0003164856570000041
其中，
Figure GDA0003164856570000042
其中，η为小推力器效率，P为小推力器功率，m为航天器质量，g为重力加速度，Isp为小推力器比冲，GMe为地球引力常数，值为3.986×1014，Re为地球半径，大小为6378137m，
Figure GDA0003164856570000043
其中CD,
Figure GDA0003164856570000044
为卫星阻力系数和面质比，
Figure GDA0003164856570000045
σ＝7.292×10-5为地球自转速度，ρ为卫星所在轨道高度的大气密度。
在本申请的一种示例性实施例中，步骤S14中，构建公式(16)所示的目标函数J：
Figure GDA0003164856570000046
其中，
Figure GDA0003164856570000047
为在第k次迭代的粒子i的控制率所产生的当前轨道根数与目标轨道根数的差的绝对值，λj,rj分别为拉格朗日乘子和惩罚因子，
Figure GDA0003164856570000048
为第k次迭代的粒子i的控制率的产生的离轨时间，是基于所述位置参数和摄动加速度，将所述控制率进行轨道积分得到的。
在本申请的一种示例性实施例中，在进行轨道积分的过程中，若存在轨道高度增加和/或偏心率偏离阈值的情况，则将得到的积分结果删除。
在本申请的一种示例性实施例中，在步骤S2中，通过公式(18)更新各粒子的位置参数和速度：
Figure GDA0003164856570000051
其中，
Figure GDA0003164856570000052
其中，某一粒子在第k次迭代中的位置和速度分别为
Figure GDA0003164856570000053
和
Figure GDA0003164856570000054
ω为惯性权重，
Figure GDA0003164856570000055
为个体在过去k次迭代过程中的历史最佳位置，
Figure GDA0003164856570000056
为粒子群在k轮迭代后种群历史最佳位置，m1为[0,1]之间的随机数，ρk为随机方向参数，ns,nf表示当前种群历史最佳位置
Figure GDA0003164856570000057
与上一次迭代
Figure GDA0003164856570000058
的关系，若两者相同，nf累加一次，否则ns累加一次，s,f分别表示ns,nf的阈值；
其中，在历次迭代中，将个体粒子与个体历史值对比，将粒子群与种群历史值对比，得到的最小适应度对应的粒子位置分别为
Figure GDA0003164856570000059
和
Figure GDA00031648565700000510
在本申请的一种示例性实施例中，在步骤S14中，通过公式(20)更新所述拉格朗日乘子和惩罚因子：
Figure GDA00031648565700000511
其中，
Figure GDA0003164856570000061
的值等于步骤S2得到的所述种群最优适应度
Figure GDA0003164856570000062
εf为约束阈值。
在本申请的一种示例性实施例中，所述轨道根数包括半长轴a、偏心率e、轨道倾角i、升交点赤经Ω、近地点角距ω和真近点角θ。
在本申请的一种示例性实施例中，所述目标适应度和所述种群适应度为最小目标函数值对应的适应度。
根据本申请的另一个方面，提供一种基于粒子群算法的低轨卫星离轨控制系统，所述系统包括：
构建单元，用于根据粒子群中粒子的位置参数和摄动加速度，以及预设的拉格朗日乘子和惩罚因子，构建目标函数；
更新单元，用于基于所述目标函数，计算所述粒子群中各粒子的目标函数值，根据计算结果确定每个粒子的目标适应度和所述粒子群的种群适应度，并更新各粒子的位置参数和速度，以及获取最大迭代次数时的种群最优适应度和其对应的粒子的位置参数；
处理单元，用于根据所述种群最优适应度更新所述拉格朗日乘子和惩罚因子，并获取当前迭代次数和基于所述种群最优适应度的目标偏差值，将所述当前迭代次数和所述目标偏差值与预设迭代条件和预设收敛条件进行比较，根据比较结果确定目标种群适应度和目标位置参数，得到目标离轨方式。
本申请的有益效果：
(1)本申请所提供的基于粒子群算法的低轨卫星离轨控制方法，相较于传统打靶法对初值不敏感，在给定搜索范围后就能进行搜索，大大减少了初值确定的困难。
(2)本申请所提供的基于粒子群算法的低轨卫星离轨控制方法，针对粒子群算法的位置、速度迭代公式和惯性权重迭代公式进行了改进，减少了算法陷入局部最优解的可能性，提高了算法的整体搜索能力和全局收敛性。
(3)本申请所提供的基于粒子群算法的低轨卫星离轨控制方法，采用增广拉格朗日方法引入了二次惩罚项，进而通过迭代校正拉格朗日乘子，避免了罚函数法因参数值过大引起的病态问题，并将一些约束通过轨道积分程序转化到适应度函数中，减少了目标函数的约束并提高了程序的计算速度。利用该方法可以获得使低轨卫星小推力离轨时间最短的最优控制率，具有对初值不敏感、收敛速度快的优点，为低轨卫星小推力离轨的工程应用提供了有效方案。
附图说明
此处的附图被并入说明书中并构成本说明书的一部分，示出了符合本申请的实施例，并与说明书一起用于解释本申请的原理。显而易见地，下面描述中的附图仅仅是本申请的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1是本申请示例性实施方式的基于粒子群算法的低轨卫星离轨控制方法流程图；
图2是本申请示例性实施方式的增广拉格朗日粒子群算法的外层程序框图；
图3是本申请示例性实施方式的粒子群算法子程序程序框图；
图4是本申请示例性实施方式的小推力器推力方向示意图；
图5是本申请示例性实施方式的半长轴在自然条件下随时间变化的结果；
图6是本申请示例性实施方式的第一种处置轨道程序迭代过程；
图7是本申请示例性实施方式的第二种处置轨道程序迭代过程；
图8是本申请示例性实施方式的最优控制率作用下卫星半长轴、偏心率、轨道倾角变化(第一种处置轨道)；
图9是本申请示例性实施方式的最优控制率作用下卫星半长轴、偏心率、轨道倾角变化(第二种处置轨道)。
具体实施方式
为了使本发明的目的、技术方案及有益效果更加清楚明白，下面结合附图及实施例，对本发明进行进一步详细说明。应当注意，此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。
本申请的示例实施方式中首先提供一种基于粒子群算法的低轨卫星离轨控制方法，首先用粒子的三个维度代表三个协状态变量，在可行域范围内，随机初始化种群中所有粒子的位置和速度，并分别将拉格朗日乘子和惩罚因子初始化；而后将目标函数作为粒子群算法的适应度函数，针对每个粒子计算适应度函数值，更新每个粒子的最优适应度和种群最优适应度，并更新粒子的位置和速度，直到达到最大迭代次数；最后更新增广拉格朗日函数参数，生成新的目标函数，重复上一步直到满足迭代次数和约束阈值要求，输出优化结果。该过程结合哈密尔顿函数和增广拉格朗日粒子群算法，在增广拉格朗日函数的作用下，粒子群算法经多次迭代后可得出最短离轨时间所对应的协状态变量，进而得出最优控制率。
具体而言，图1示出了本申请示例性实施方式的基于粒子群算法的低轨卫星离轨控制方法流程图，图2示出了本申请示例性实施方式的增广拉格朗日粒子群算法的外层程序框图，下面结合图1和图2对本申请实施方式的方法进行说明：
步骤S1、根据粒子群中粒子的位置参数和摄动加速度，以及预设的拉格朗日乘子和惩罚因子，构建目标函数；
具体的，包括：
步骤S11、获取在推力器推力的作用下，各轨道根数随时间的变化率，并根据所述变化率确定含有当前协状态变量的控制率；
在本实施方式中，轨道根数包括半长轴a、偏心率e、轨道倾角i、升交点赤经Ω、近地点角距ω和真近点角θ，通过如下步骤得到含有当前协状态变量的控制率：
步骤S111、获取在推力器推力的作用下，各轨道根数随时间的变化率；具体地，推力器推力对各轨道根数的影响的计算方法为：
Figure GDA0003164856570000081
Figure GDA0003164856570000082
Figure GDA0003164856570000091
Figure GDA0003164856570000092
Figure GDA0003164856570000093
Figure GDA0003164856570000094
其中，p＝a(1-e2),
Figure GDA0003164856570000095
μ＝3.986×1014，p,h,r,v,μ如上式所示，v为卫星速度，μ为地球引力常数，an,t,h为密切轨道径向，速度方向和角动量方向的摄动加速度，沿径向向外，沿速度方向和密切轨道角动量方向为正。
步骤S112、基于所述变化率和所述当前协状态变量，通过公式(7)得到哈密尔顿函数；
在本实施方式中，哈密尔顿函数的计算方法为：
Figure GDA0003164856570000096
其中，H为哈密尔顿函数值，λa、λe和λi为当前协状态变量，a、e和i分别为工作轨道的半长轴、偏心率和轨道倾角。
步骤S113、根据所述哈密尔顿函数对所述推力器推力对应的预设角度求导，得到所述含有当前协状态变量的控制率；
在本实施方式中，图4是出了本申请示例性实施方式的小推力器推力方向示意图，下面结合图4对本实施方式确定控制率的过程进行说明；其中，计算如图4所示角度α和β的含有当前协状态变量的控制率的过程为，将哈密尔顿函数分别对α和β求偏导，并联立方程，结果为：
Figure GDA0003164856570000101
Figure GDA0003164856570000102
步骤S12、计算推力器推力的推力矢量、燃料消耗率及摄动加速度；
在本实施方式中，计算图4所示小推力矢量和燃料消耗率，并计算地球非球形J2摄动和大气阻力摄动的摄动加速度，具体而言，
步骤S121、推力器的推力矢量计算方法为：
Figure GDA0003164856570000103
Figure GDA0003164856570000104
燃料消耗率的计算方法为：
Figure GDA0003164856570000105
其中，η为小推力器效率，P为小推力器功率，m为航天器质量，g为重力加速度，Isp为小推力器比冲。
步骤S122、计算地球非球形J2摄动和大气阻力摄动的摄动加速度；
摄动加速度计算方法为：
Figure GDA0003164856570000111
Figure GDA0003164856570000112
其中，GMe为地球引力常数，值为3.986×1014，Re为地球半径，大小为6378137m，
Figure GDA0003164856570000113
其中CD,
Figure GDA0003164856570000114
为卫星阻力系数和面质比，
Figure GDA0003164856570000115
σ＝7.292×10-5为地球自转速度，ρ为卫星所在轨道高度的大气密度。
步骤S13、根据所述当前协状态变量确定粒子位置向量，并将所述粒子位置向量作为粒子群中粒子的位置参数；
在本实施方式中，将协状态变量作为粒子位置向量，该粒子位置向量的计算方法为：
λi＝[λa λe λi] (15)。
步骤S14、基于所述位置参数和摄动加速度，将所述控制率进行轨道积分，并根据积分结果、预设的拉格朗日乘子和惩罚因子，构建目标函数；
在本实施方式中，构建公式(16)所示的目标函数J：
Figure GDA0003164856570000116
其中，
Figure GDA0003164856570000117
为在第k次迭代的粒子i的控制率所产生的当前轨道根数与目标轨道根数的差的绝对值，λj,rj分别为拉格朗日乘子和惩罚因子，
Figure GDA0003164856570000121
为第k次迭代的粒子i的控制率的产生的离轨时间，是基于所述位置参数和摄动加速度(由式(1)至(6)积分实现)，将控制率进行轨道积分得到的，λj,rj在子问题运算过程中保持不变。
需要说明的是，一些约束条件通过相应的处理可以在轨道积分过程中实施。由于粒子的初始位置值是随机的，某些粒子在轨道积分时可能会造成轨道高度增加和偏心率异常等异常状况，针对这一问题，可以在轨道积分过程中设置中断条件，并将离轨时间返回为无穷大值。这些粒子的适应度值也就变为无穷大，从而在适应度筛选中淘汰。
步骤S2、基于所述目标函数，计算所述粒子群中各粒子的目标函数值，根据计算结果确定每个粒子的目标适应度和所述粒子群的种群适应度，并更新各粒子的位置参数和速度，以及获取最大迭代次数时的种群最优适应度和其对应的粒子的位置参数；其中，种群最优适应度为在满足预设的最大迭代次数时输出的种群适应度。
在本实施方式中，目标适应度和种群适应度为最小目标函数值对应的适应度。
通过比较种群中各粒子的目标函数值(最小为最优)得出每个粒子的最优适应度
Figure GDA0003164856570000122
和种群最优适应度
Figure GDA0003164856570000123
并更新每个粒子的位置和速度，在满足预设迭代次数时输出子程序优化结果(参见图3所示)；
其中，通过公式(18)更新各粒子的位置参数和速度：
Figure GDA0003164856570000124
其中，
Figure GDA0003164856570000125
其中，某一粒子在第k次迭代中的位置和速度分别为
Figure GDA0003164856570000126
和
Figure GDA0003164856570000127
ω为惯性权重，
Figure GDA0003164856570000131
为个体在过去k次迭代过程中的历史最佳位置，
Figure GDA0003164856570000132
为粒子群在k轮迭代后种群历史最佳位置，m1为[0,1]之间的随机数，ρk为随机方向参数，ns,nf表示当前种群历史最佳位置
Figure GDA0003164856570000133
与上一次迭代
Figure GDA0003164856570000134
的关系，若两者相同，nf累加一次，否则ns累加一次，s,f分别表示ns,nf的阈值；
其中，在历次迭代中，将个体粒子与个体历史值对比，将粒子群与种群历史值对比，得到的最小适应度对应的粒子位置分别为
Figure GDA0003164856570000135
和
Figure GDA0003164856570000136
步骤S3、根据所述种群最优适应度更新所述拉格朗日乘子和惩罚因子，并获取当前迭代次数和基于所述种群最优适应度的目标偏差值，将所述当前迭代次数和所述目标偏差值与预设迭代条件和预设收敛条件进行比较，根据比较结果确定目标种群适应度和目标位置参数，得到目标离轨方式；其中，基于所述种群最优适应度的目标偏差值为步骤S2中获取的种群最优适应度对应的位置参数积分产生的半长轴或近地点与目标的偏差值。
在本实施方式中，若当前迭代次数满足预设迭代条件且所述目标偏差值满足收敛条件，则将种群最优适应度确定为目标种群适应度，将目标种群适应度对应的位置参数确定为目标位置参数；
若当前迭代次数不满足预设迭代条件或所述目标偏差值不满足收敛条件，则基于更新后的拉格朗日乘子和惩罚因子更新目标函数，并继续执行步骤S14及后续步骤，直至所述当前迭代次数满足所述预设迭代条件，且所述目标偏差值满足所述收敛条件。
其中，更新拉格朗日乘子
Figure GDA0003164856570000137
和惩罚因子
Figure GDA0003164856570000138
的计算方法为：
Figure GDA0003164856570000139
式中：
Figure GDA00031648565700001310
的值等于步骤S2得到的所述种群最优适应度
Figure GDA00031648565700001311
(为图3所示的子问题得到的最优解)，εf为约束阈值。
其中，若最优解
Figure GDA00031648565700001312
不满足约束条件，则返回步骤S14，利用更新后的拉格朗日乘子和惩罚因子计算得出新的目标函数，进而进行后续步骤，直到得出满足约束条件εf的结果。
本申请的基于粒子群算法的低轨卫星离轨控制方法，首先由哈密尔顿函数推导出带协状态变量的控制率；然后将协状态变量向量作为粒子群中粒子的位置参数，以当前目标函数运行粒子群算法子程序(图3)，并在每次迭代中更新每个粒子的最优适应度
Figure GDA0003164856570000141
种群最优适应度
Figure GDA0003164856570000142
和粒子的位置和速度，当满足最大迭代次数时，输出优化结果；最后依据子程序得出的优化结果，更新增广拉格朗日函数中的参数，进而得出新的目标函数用于下一次子程序迭代，当满足迭代次数和约束阈值要求时，输出优化结果。该方法可以得出使低轨卫星小推力离轨时间最短的最优控制率，具有对初值不敏感、收敛速度快的优点，为低轨卫星小推力离轨的工程应用提供了有效方案。
实施例1
下面结合本实施例对本申请的基于粒子群算法的低轨卫星离轨控制方法进行说明，本实施例中，通过以下步骤实现：步骤一，给定航天器自身参数和小推力器参数；步骤二，给定航天器工作轨道和两种处置轨道半长轴、偏心率和轨道倾角；步骤三，给定增广拉格朗日粒子群算法参数；步骤四，计算在无推力作用下离轨时间；步骤五，计算两种处置轨道的最优控制率、离轨时间和燃料消耗。
其中，航天器自身参数和小推力器参数见表1，航天器工作轨道和两种处置轨道半长轴、偏心率和轨道倾角见表2，增广拉格朗日粒子群算法参数见表3，两种处置轨道最优控制率、离轨时间和燃料消耗见表4。
表1航天器参数
参数 数值 参数 数值
卫星质量 1000.0kg 比冲 4660s
大气阻力系数 2 推力器功率 115W
面质比 0.02m<sup>2</sup>/kg 发动机推力 3mN
推力器效率 24.5％ 大气模型 SA76
表2工作轨道与处置轨道参数
Figure GDA0003164856570000151
表3增广拉格朗日粒子群算法参数
算法参数 数值
种群个数N 500
最大迭代次数k<sub>max</sub> 50
最大惯性权重ω<sub>max</sub> 0.9
最小惯性权重ω<sub>min</sub> 0.4
n<sub>s</sub>阈值s 15
n<sub>f</sub>阈值f 5
随机方向参数初值ρ<sub>0</sub> 1
初始惩罚因子r<sub>0</sub> 1000
约束阈值ε<sub>f</sub> 0.0001
表4两种处置轨道最优控制率、离轨时间和燃料消耗
Figure GDA0003164856570000152
Figure GDA0003164856570000161
图5示出了本申请实施例1中的无推力作用下半长轴随时间变化图，由图5可知，在自然条件下无推力离轨时间为140年，大大超出25年内离轨的要求。图6和图7分别示出了两种处置轨道迭代结果，由图6可得：过程中总共迭代60次，在第57次时满足半长轴偏差收敛条件；由图7可得：过程中总共迭代25次，在第12次时满足近地点偏差收敛条件。图8和图9分别示出了两种处置轨道半长轴、偏心率和轨道倾角随时间变化图，图8和图9可知：航天器均在本申请得出的最优控制率作用下达到了相应的目标轨道(处置轨道)，基于此，证明了本申请的算法得出的最优控制率的有效性。
本申请另一示例实施例提供一种基于粒子群算法的低轨卫星离轨控制系统，包括：构建单元，用于根据粒子群中粒子的位置参数和摄动加速度，以及预设的拉格朗日乘子和惩罚因子，构建目标函数；更新单元，用于基于所述目标函数，计算所述粒子群中各粒子的目标函数值，根据计算结果确定每个粒子的目标适应度和所述粒子群的种群适应度，并更新各粒子的位置参数和速度，以及获取最大迭代次数时的种群最优适应度和其对应的粒子的位置参数；处理单元，用于根据所述种群最优适应度更新所述拉格朗日乘子和惩罚因子，并获取当前迭代次数和基于所述种群最优适应度的目标偏差值，将所述当前迭代次数和所述目标偏差值与预设迭代条件和预设收敛条件进行比较，根据比较结果确定目标种群适应度和目标位置参数，得到目标离轨方式。
上述控制系统中各个单元的具体描述可以参考控制方法中对每个步骤的描述，在此不再赘述，上述控制系统可以实现与控制方法侧同样的功能。
以上，仅是本发明的一个实施例，并非对本发明做任何形式的限制，虽然本发明以较佳实施例揭示如上，然而并非用以限制本发明，任何熟悉本专业的技术人员，在不脱离本发明技术方案的范围内，利用上述揭示的技术内容做出些许的变动或修饰均等同于等效实施案例，均属于技术方案范围内。

Claims (11)

1.一种基于粒子群算法的低轨卫星离轨控制方法，其特征在于，包括：
步骤S1、根据粒子群中粒子的位置参数和摄动加速度，以及预设的拉格朗日乘子和惩罚因子，构建目标函数；
步骤S2、基于所述目标函数，计算所述粒子群中各粒子的目标函数值，根据计算结果确定每个粒子的目标适应度和所述粒子群的种群适应度，并更新各粒子的位置参数和速度，以及获取最大迭代次数时的种群最优适应度和其对应的粒子的位置参数；
步骤S3、根据所述种群最优适应度更新所述拉格朗日乘子和惩罚因子，并获取当前迭代次数和基于所述种群最优适应度的目标偏差值，将所述当前迭代次数和所述目标偏差值与预设迭代条件和预设收敛条件进行比较，根据比较结果确定目标种群适应度和目标位置参数，得到目标离轨方式；
其中，所述根据粒子群中粒子的位置参数和摄动加速度，以及预设的拉格朗日乘子和惩罚因子，构建目标函数，包括：
步骤S11、获取在推力器推力的作用下，各轨道根数随时间的变化率，并根据所述变化率确定含有当前协状态变量的控制率；
步骤S12、计算所述推力器推力的推力矢量、燃料消耗率及摄动加速度；
步骤S13、根据所述当前协状态变量确定粒子位置向量，并将所述粒子位置向量作为粒子群中粒子的位置参数；
步骤S14、基于所述位置参数和摄动加速度，将所述控制率进行轨道积分，并根据积分结果、预设的拉格朗日乘子和惩罚因子，构建所述目标函数。
2.根据权利要求1所述的方法，其特征在于，所述将所述当前迭代次数和所述目标偏差值与预设迭代条件和预设收敛条件进行比较，根据比较结果确定目标种群适应度和目标位置参数，包括：
若所述当前迭代次数满足预设迭代条件且所述目标偏差值满足收敛条件，则将所述种群最优适应度确定为目标种群适应度，将所述目标种群适应度对应的位置参数确定为目标位置参数；
若所述当前迭代次数不满足所述预设迭代条件或所述目标偏差值不满足所述收敛条件，则基于更新后的拉格朗日乘子和惩罚因子更新目标函数，并继续执行步骤S14及后续步骤，直至所述当前迭代次数满足所述预设迭代条件，且所述目标偏差值满足所述收敛条件。
3.根据权利要求1所述的方法，其特征在于，步骤S11中，所述根据所述变化率确定含有当前协状态变量的控制率，包括：
基于所述变化率和所述当前协状态变量，通过公式(7)得到哈密尔顿函数：
Figure FDA0003211063550000021
其中，H为哈密尔顿函数值，λa、λe和λi为当前协状态变量，a、e和i分别为工作轨道的半长轴、偏心率和轨道倾角；
根据所述哈密尔顿函数对所述推力器推力对应的预设角度求导，得到所述含有当前协状态变量的控制率。
4.根据权利要求1所述的方法，其特征在于，步骤S12中，通过公式(11)得到推力矢量，通过公式(12)得到燃料消耗率，通过公式(14)得到摄动加速度；
Figure FDA0003211063550000022
Figure FDA0003211063550000023
Figure FDA0003211063550000024
其中，
Figure FDA0003211063550000031
其中，η为小推力器效率，P为小推力器功率，m为航天器质量，g为重力加速度，Isp为小推力器比冲，GMe为地球引力常数，值为3.986×1014，Re为地球半径，大小为6378137m，
Figure FDA0003211063550000032
其中
Figure FDA0003211063550000033
为卫星阻力系数和面质比，
Figure FDA0003211063550000034
σ＝7.292×10-5为地球自转速度，ρ为卫星所在轨道高度的大气密度；an为密切轨道径向的摄动加速度，沿径向向外的方向为正；at为密切轨道速度方向的摄动加速度，沿速度方向的方向为正；ah为密切轨道角动量方向的摄动加速度，沿密切轨道角动量方向为正；p＝a(1-e2)，
Figure FDA0003211063550000035
v为卫星速度；μ＝3.986005×1014m3/s2为地球引力常数；a为轨道半长轴；e为偏心率；i为轨道倾角；w为近地点幅角；θ为真近点角。
5.根据权利要求1所述的方法，其特征在于，步骤S14中，构建公式(16)所示的目标函数J：
Figure FDA0003211063550000036
其中，
Figure FDA0003211063550000037
Figure FDA0003211063550000038
为在第k次迭代的粒子i的控制率所产生的当前轨道根数与目标轨道根数的差的绝对值，λj,rj分别为拉格朗日乘子和惩罚因子，
Figure FDA0003211063550000041
为第k次迭代的粒子i的控制率的产生的离轨时间，是基于所述位置参数和摄动加速度，将所述控制率进行轨道积分得到的。
6.根据权利要求5所述的方法，其特征在于，在进行轨道积分的过程中，若存在轨道高度增加和/或偏心率偏离阈值的情况，则将得到的积分结果删除。
7.根据权利要求1所述的方法，其特征在于，在步骤S2中，通过公式(18)更新各粒子的位置参数和速度：
Figure FDA0003211063550000042
其中，
Figure FDA0003211063550000043
其中，某一粒子在第k次迭代中的位置和速度分别为
Figure FDA0003211063550000044
和
Figure FDA0003211063550000045
ω为惯性权重，
Figure FDA0003211063550000046
为个体在过去k次迭代过程中的历史最佳位置，
Figure FDA0003211063550000047
为粒子群在k轮迭代后种群历史最佳位置，m1为[0，1]之间的随机数，ρk为随机方向参数，ns，nf表示当前种群历史最佳位置
Figure FDA0003211063550000048
与上一次迭代
Figure FDA0003211063550000049
的关系，若两者相同，nf累加一次，否则ns累加一次，s，f分别表示ns，nf的阈值；
其中，在历次迭代中，将个体粒子与个体历史值对比，将粒子群与种群历史值对比，得到的最小适应度对应的粒子位置分别为
Figure FDA00032110635500000410
和
Figure FDA00032110635500000411
k为迭代轮数；ρk+1为迭代k+1轮时随机方向参数。
8.根据权利要求7所述的方法，其特征在于，在步骤S14中，通过公式(20)更新所述拉格朗日乘子和惩罚因子：
Figure FDA0003211063550000051
其中，
Figure FDA0003211063550000052
的值等于步骤S2得到的所述种群最优适应度
Figure FDA0003211063550000053
εf为约束阈值；
Figure FDA0003211063550000054
为对应轨道根数迭代n+1次的拉格朗日乘子；
Figure FDA0003211063550000055
为对应轨道根数迭代n次的拉格朗日乘子；
Figure FDA0003211063550000056
为对应轨道根数迭代n+1次的惩罚因子；
Figure FDA0003211063550000057
为对应轨道根数迭代n次的惩罚因子；
Figure FDA0003211063550000058
为迭代n次的最优适应度粒子的控制率所产生的对应轨道根数与目标轨道根数的差的绝对值。
9.根据权利要求2～8任一项所述的方法，其特征在于，所述轨道根数包括半长轴a、偏心率e、轨道倾角i、升交点赤经Ω、近地点角距ω和真近点角θ。
10.根据权利要求1所述的方法，其特征在于，所述目标适应度和所述种群适应度为最小目标函数值对应的适应度。
11.一种基于粒子群算法的低轨卫星离轨控制系统，其特征在于，包括：
构建单元，用于根据粒子群中粒子的位置参数和摄动加速度，以及预设的拉格朗日乘子和惩罚因子，构建目标函数；
更新单元，用于基于所述目标函数，计算所述粒子群中各粒子的目标函数值，根据计算结果确定每个粒子的目标适应度和所述粒子群的种群适应度，并更新各粒子的位置参数和速度，以及获取最大迭代次数时的种群最优适应度和其对应的粒子的位置参数；
处理单元，用于根据所述种群最优适应度更新所述拉格朗日乘子和惩罚因子，并获取当前迭代次数和基于所述种群最优适应度的目标偏差值，将所述当前迭代次数和所述目标偏差值与预设迭代条件和预设收敛条件进行比较，根据比较结果确定目标种群适应度和目标位置参数，得到目标离轨方式；
其中，所述构建单元具体用于：
获取在推力器推力的作用下，各轨道根数随时间的变化率，并根据所述变化率确定含有当前协状态变量的控制率；
计算所述推力器推力的推力矢量、燃料消耗率及摄动加速度；
根据所述当前协状态变量确定粒子位置向量，并将所述粒子位置向量作为粒子群中粒子的位置参数；
基于所述位置参数和摄动加速度，将所述控制率进行轨道积分，并根据积分结果、预设的拉格朗日乘子和惩罚因子，构建所述目标函数。
CN202010629155.7A 2020-04-30 2020-07-01 基于粒子群算法的低轨卫星离轨控制方法及系统 Active CN111874267B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
US17/363,363 US11787567B2 (en) 2020-04-30 2021-06-30 Low-orbit satellite deorbit control method and system based on particle swarm algorithm

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
CN202010367858 2020-04-30
CN2020103678587 2020-04-30

Publications (2)

Publication Number Publication Date
CN111874267A CN111874267A (zh) 2020-11-03
CN111874267B true CN111874267B (zh) 2021-09-28

Family

ID=73150118

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010629155.7A Active CN111874267B (zh) 2020-04-30 2020-07-01 基于粒子群算法的低轨卫星离轨控制方法及系统

Country Status (2)

Country Link
US (1) US11787567B2 (zh)
CN (1) CN111874267B (zh)

Families Citing this family (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112415545A (zh) * 2020-11-23 2021-02-26 许昌学院 海上多gnss系统粒子群选星方法及系统
CN113343369B (zh) * 2021-08-06 2021-10-01 中国空气动力研究与发展中心设备设计与测试技术研究所 一种航天器气动融合轨道摄动分析方法
CN113867379B (zh) * 2021-10-08 2022-05-10 北京理工大学 一种棱锥型离轨帆的构形构建与姿态控制方法
CN114852375A (zh) * 2022-03-24 2022-08-05 北京控制工程研究所 一种编队卫星相对轨道变化估计方法、估计装置
CN114826376B (zh) * 2022-03-30 2023-11-10 大连大学 一种微纳卫星多目标优化资源分配方法
CN117361359B (zh) * 2022-06-21 2024-04-12 张家港市申联建设机械有限公司 一种塔机起升速度确定方法和系统
CN115347940B (zh) * 2022-08-18 2023-06-20 中国星网网络应用有限公司 带惩罚机制的卫星通讯时间窗规划方法及系统
CN115583367B (zh) * 2022-10-09 2023-05-30 哈尔滨工业大学 一种基于分布式序列凸优化的卫星集群重构控制方法
CN116341390B (zh) * 2023-05-11 2023-11-17 西安现代控制技术研究所 一种全局搜索快速收敛多约束弹道优化方法
CN116647270B (zh) * 2023-07-19 2023-09-29 中南大学 基于自适应粒子群算法的双层卫星网络分簇方法和系统
CN116756469B (zh) * 2023-08-22 2023-10-31 中之力搏建设工程有限公司 一种户外照明灯具优化管理系统

Family Cites Families (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103971174B (zh) * 2014-05-06 2017-04-12 大连理工大学 基于改进量子粒子群算法的水电站群优化调度方法
CN108710939A (zh) * 2018-04-25 2018-10-26 广西民族大学 一种多目标过程系统优化的粒子群算法求解方法
US20200019435A1 (en) * 2018-07-13 2020-01-16 Raytheon Company Dynamic optimizing task scheduling
CN109977576A (zh) * 2019-04-03 2019-07-05 北京理工大学 一种基于代理模型的卫星星座系统多学科设计优化方法
CN110083170A (zh) * 2019-04-11 2019-08-02 北京航空航天大学 一种使用固体微推力器进行轨道保持的优化控制方法
CN110769512B (zh) * 2019-10-31 2023-04-18 哈尔滨工业大学（深圳） 星载资源分配方法、装置、计算机设备及存储介质
CN111738395A (zh) * 2020-06-01 2020-10-02 大连理工大学 一种求解含等式约束优化问题的双适应度混沌粒子群算法

Also Published As

Publication number Publication date
US11787567B2 (en) 2023-10-17
CN111874267A (zh) 2020-11-03
US20220002006A1 (en) 2022-01-06

Similar Documents

Publication Publication Date Title
CN111874267B (zh) 基于粒子群算法的低轨卫星离轨控制方法及系统
Pontani et al. Particle swarm optimization applied to space trajectories
CN109870906B (zh) 一种基于bbo优化人工势场的高速旋翼飞行器路径规划方法
CN106021784A (zh) 一种基于两层优化策略的全轨迹优化设计方法
CN110059285B (zh) 考虑j2项影响的导弹自由段弹道偏差解析预报方法
Liu et al. Collision-free trajectory design for long-distance hopping transfer on asteroid surface using convex optimization
Zhu et al. Fast approximation of optimal perturbed long-duration impulsive transfers via artificial neural networks
Zandavi Surface-to-air missile path planning using genetic and PSO algorithms
Rao et al. Numerical optimization study of multiple‐pass aeroassisted orbital transfer
Burchett Fuzzy logic trajectory design and guidance for terminal area energy management
Izzo et al. Real-Time Optimal Guidance and Control for Interplanetary Transfers Using Deep Networks
Shornikov et al. Simulation of controlled motion in an irregular gravitational field for an electric propulsion spacecraft
CN113955153B (zh) 一种燃料最优的连续小推力轨道转移方法
Bau et al. DROMO: a new regularized orbital propagator
Colombo Optimal trajectory design for interception and deflection of Near Earth Objects
Wang et al. Satellite Constellation Reconfiguration Method Using Improved PSO Algorithm for Communication and Navigation
Dracopoulos et al. Neuro-genetic adaptive attitude control
Mandalia et al. High-Altitude Divert Architecture for Future Robotic and Human Mars Missions
Hamed et al. Optimized curvilinear potential field based multi-objective satellite collision avoidance maneuver
Zhu et al. The Intelligent Trajectory Optimization of Multistage Rocket with Gauss Pseudo-Spectral Method.
Cheng et al. Trajectory Correction of the Rocket for Aerodynamic Load Shedding Based on Deep Neural Network and the Chaotic Evolution Strategy with Covariance Matrix Adaptation
Pi et al. Reinforcement learning trajectory generation and control for aggressive perching on vertical walls with quadrotors
Dong et al. Trajectory online optimization for unmanned combat aerial vehicle using combined strategy
Dong et al. Knowledge-Driven Accurate Opponent Trajectory Prediction for Gun-Dominated Autonomous Air Combat
Ionasescu Orbit determination covariance analysis for the cassini solstice mission

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