CN114169072A - 一种复杂约束下航天器连续推力轨道转移控制及优化方法 - Google Patents
一种复杂约束下航天器连续推力轨道转移控制及优化方法 Download PDFInfo
- Publication number
- CN114169072A CN114169072A CN202111473972.9A CN202111473972A CN114169072A CN 114169072 A CN114169072 A CN 114169072A CN 202111473972 A CN202111473972 A CN 202111473972A CN 114169072 A CN114169072 A CN 114169072A
- Authority
- CN
- China
- Prior art keywords
- spacecraft
- orbit
- transfer
- constraint
- constraints
- 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.)
- Granted
Links
- 238000012546 transfer Methods 0.000 title claims abstract description 115
- 238000005457 optimization Methods 0.000 title claims abstract description 69
- 238000000034 method Methods 0.000 title claims abstract description 64
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 56
- PEDCQBHIVMGVHV-UHFFFAOYSA-N Glycerine Chemical compound OCC(O)CO PEDCQBHIVMGVHV-UHFFFAOYSA-N 0.000 claims abstract description 25
- 238000011217 control strategy Methods 0.000 claims abstract description 22
- 238000005259 measurement Methods 0.000 claims abstract description 22
- 241001183967 Isodon Species 0.000 claims abstract description 4
- 238000004364 calculation method Methods 0.000 claims description 12
- 239000000446 fuel Substances 0.000 claims description 9
- 230000002068 genetic effect Effects 0.000 claims description 5
- 238000013507 mapping Methods 0.000 claims description 5
- 238000006243 chemical reaction Methods 0.000 claims description 3
- 230000000694 effects Effects 0.000 claims description 3
- 239000002245 particle Substances 0.000 claims description 3
- 230000006870 function Effects 0.000 description 21
- 238000004891 communication Methods 0.000 description 9
- 230000008569 process Effects 0.000 description 9
- 230000001133 acceleration Effects 0.000 description 6
- 238000013461 design Methods 0.000 description 6
- 238000004088 simulation Methods 0.000 description 6
- 238000004458 analytical method Methods 0.000 description 5
- 230000008859 change Effects 0.000 description 5
- 238000009795 derivation Methods 0.000 description 3
- 230000005484 gravity Effects 0.000 description 3
- 238000011160 research Methods 0.000 description 3
- 238000011161 development Methods 0.000 description 2
- 230000010354 integration Effects 0.000 description 2
- 239000011159 matrix material Substances 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 230000003287 optical effect Effects 0.000 description 2
- 230000007704 transition Effects 0.000 description 2
- 241000959638 Geastrum saccatum Species 0.000 description 1
- 241000283973 Oryctolagus cuniculus Species 0.000 description 1
- 230000009471 action Effects 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 238000013473 artificial intelligence Methods 0.000 description 1
- 238000013528 artificial neural network Methods 0.000 description 1
- 230000001174 ascending effect Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 238000005286 illumination Methods 0.000 description 1
- 238000012933 kinetic analysis Methods 0.000 description 1
- 238000010801 machine learning Methods 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 230000002035 prolonged effect Effects 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
- 238000004800 variational method Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/10—Geometric CAD
- G06F30/15—Vehicle, aircraft or watercraft design
-
- 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
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/27—Design optimisation, verification or simulation using machine learning, e.g. artificial intelligence, neural networks, support vector machines [SVM] or training a model
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/04—Constraint-based CAD
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Geometry (AREA)
- General Physics & Mathematics (AREA)
- Evolutionary Computation (AREA)
- Remote Sensing (AREA)
- Aviation & Aerospace Engineering (AREA)
- Computer Hardware Design (AREA)
- General Engineering & Computer Science (AREA)
- Software Systems (AREA)
- Medical Informatics (AREA)
- Automation & Control Theory (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Computational Mathematics (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Pure & Applied Mathematics (AREA)
- Chemical & Material Sciences (AREA)
- Combustion & Propulsion (AREA)
- Radar, Positioning & Navigation (AREA)
- Artificial Intelligence (AREA)
- Control Of Position, Course, Altitude, Or Attitude Of Moving Bodies (AREA)
Abstract
本发明涉及一种复杂约束下航天器连续推力轨道转移控制及优化方法,根据航天器轨道动力学特性,设计由航天器当前轨道位置和可调参数集确定的轨道转移控制策略,得到轨道转移控制模型;结合地球自转模型和日、月星历模型,引入航天器在实际任务中需要面对的开关机约束,得到带有复杂约束的轨道转移控制模型;所述开关机约束包括地影区关机约束、月影区关机约束、地面测控站可见时开机约束及日凌现象发生时关机约束;将轨道转移策略中的可调参数集作为输入,建立控制性能指标作为输出,得到带有复杂约束的轨道转移控制优化模型;利用智能优化算法输入的可调参数集进行参数优化,最终完成复杂约束下航天器连续推力轨道转移控制及优化工作。
Description
技术领域
本发明涉及一种复杂约束下航天器连续推力轨道转移控制及优化方法,属于航天器控制技术领域。
背景技术
在当前航天任务中,采用连续推力对航天器进行轨道控制成为一种发展的趋势,尤其是电推进等新型推进技术的发展使得,可以大幅减少燃料质量,有效降低卫星总重或提高有效载荷的质量比,延长任务的寿命,降低发射成本。因此关于连续推力作用下的航天器轨道转移最优控制方法的研究受到了研究学者的广泛关注,数值求解方法可以分为直接法、间接法及两者结合的混合法。直接法将轨道转移最优控制问题进行状态变量和控制变量的离散参数化,从而转化为参数优化问题,进一步通过非线性规划算法求解。间接法利用变分法和庞特里亚金极值原理推导最优性必要条件,将最优控制问题转化为两点边值问题,通过对协态变量的初值进行搜索进而得到最优控制律。求解边值问题协态变量初值可以采用打靶法,但是打靶函数对于初值非常敏感,不容易收敛。因此同伦法、智能优化算法、UKF参数估计法[8]等手段均被用于对协态变量初值进行搜索。混合法通常假定了推力的函数结构,推导最优控制的必要条件后求解。通过对控制结构进行假设和对问题进行简化,可以得到一些轨迹优化的解析求解方法。田百义等针对GTO-GEO轨道转移问题,采用了轨道面内面外解耦控制的策略,在面内保持远地点高度不变条件下解析求解径向和法向推力夹角,面外控制依靠面内外推力夹角控制,控制角设计为纬度幅角的正弦函数,在升降交点处取最大值。朱政帆等提出了一种利用人工智能方法进行轨道转移设计的构想,自动进行初解、延拓和拼接三个步骤,发挥出机器学习高效计算得优越性。但是,现有研究中的轨道控制方法与轨道动力学模型相关性很高,且需要复杂的公式推导。
求解轨道转移问题后得到的最优控制律往往是一组复杂的数值曲线,需要设计制导控制律,使航天器实现对最优轨迹的跟踪。Yang等于李雅普诺夫理论和神经网络方法实现了考虑地球J2摄动与地影区影响的轨道转移制导控制。张冉等提出一种多项式曲线拟合的方法,逼近轨道转移的数值解,建立了一种分段常值推力跟踪参考轨道的闭环制导策略。但是,采用制导律设计需要有参考轨迹进行跟踪,不便于直接使用。
直接法的优点是不需要推导复杂的最优性条件,在状态和控制参数化时可以考虑路径约束,缺点是难以保证控制的最优性,并且需要给出控制变量的初值猜测。间接法中两点边值问题的协态变量没有物理意义,且对协态初值非常敏感,使得初值猜测的难度很大,尤其是在初末状态变量差异和转移时间增大的情况下,更加难以得到收敛的初值猜测值,并且采用间接法进行最优性必要条件推导时难以处理路径约束。同时,由于轨道转移的解析解往往难以得到,无论间接法还是直接法求解得到的控制律都是数值解,是一组自变量的数值曲线,难以直接在工程控制中应用,往往需要设计额外的制导律进行轨迹跟踪。
发明内容
本发明的技术解决问题:克服现有技术的不足,提供一种复杂约束下航天器连续推力轨道转移控制及优化方法,实现任意两条开普勒椭圆轨道间的轨道转移控制,并且在转移过程中能够满足开关机的路径约束,同时可以避免复杂的制导律设计。
本发明技术解决方案:一种复杂约束下航天器连续推力轨道转移控制及优化方法,其特点在于包括以下步骤:
步骤1:基于航天器连续推力动力学模型,提出一种由航天器当前轨道特征点邻域和可调参数集确定的轨道转移控制策略,依靠常幅值径向、切向和法向分量的发动机推力实现轨道转移;
步骤2:基于步骤1中的轨道转移控制策略,依据航天器在实际任务中需要面对的开关机约束,建立多种航天器开关机约束条件,将轨道转移控制策略和多种航天器开关机约束条件结合,得到带有复杂约束的轨道转移控制模型;所述多种航天器开关机约束包括地影区内的关机约束、月影区内的关机约束、地面测控站可见区开机约束及日凌区关机约束;
步骤3:根据步骤2得到的带有复杂约束的轨道转移控制模型,将步骤1中的可调参数集构造为状态变量输入,建立控制性能指标函数作为输出,构建状态变量和控制性能指标间的函数映射,得到带有复杂约束的轨道转移优化问题。利用智能优化算法求解带有复杂约束的轨道转移优化问题,得到使性能指标函数值最小的状态变量,最终完成复杂约束下航天器连续推力轨道转移控制及优化工作。
所述步骤1中,航天器轨道转移控制策略如下:
式中,ustr为航天器发动机推力的方向矢量,Xc和XT分别为当前轨道和目标轨道的改进春分点根数,R为当前航天器位置,通过轨道要素转换方法可以由Xc计算得到,以函数R=g(Xc)描述;
f(Xc,XT)=[ur,ut,uh]T,其中:
式中的下标r、t和h分别表示径向、切向和法向,则ur、ut和uh分别是ustr沿径向、切向和法向的分量。定义控制参数向量u0=[ur0,ut0,uh0],在D1和D2中进行面内椭圆调整,只依靠径向和切向的推力,其中径向推力用于调整近地点幅角,切向推力用于调整轨道半长轴;在D3和D4中进行轨道面调整,只依靠法向推力。ur、ut和uh的取值方法具体给出如下:
(1)若R∈D1(ξ1,δ1),若远地点高度低,则flag1=1,否则flag1=-1;若ωc<ωT,则flag4=1,否则flag4=-1,有ut=flag1*ut0,ur=flag4*ur0,flag1和flag3为符号标记;
(2)若R∈D2(ξ2,δ2),若近地点高度低,则flag2=1,否则flag2=-1;若ωc<ωT,则flag4=1,否则flag4=-1,有ut=flag2*ut0,ur=flag4*ur0,flag2为符号标记;
(3)若R∈D3(ξ3,δ3),flag3=1,uh=flag3*uh0,flag3为符号标记;
(4)若R∈D4(ξ4,δ4),flag3=-1,uh=flag3*uh0;
其中,{D1(ξ1,δ1),D2(ξ2,δ2),D3(ξ3,δ3),D4(ξ4,δ4)}为航天器当前轨道椭圆上的四个弧段,分别由四个特征点和邻域角定义,四个特征点分别是当前轨道的近地点、远地点、下穿越点和上穿越点,下穿越点为航天器由上向下经过目标轨道面的位置,反之则是上穿越点,ξ为特征向量,由地心指向特征点;δ为邻域角,用于限制对应弧段上的点与特征点关于地心的夹角。定义一个邻域角向量Δa=[δ1,δ2,δ3,δ4],由于近地点和远地点连线,上穿越点和下穿越点连线均经过地心,因此有0°<δi<90°,i=1,2,3,4。可调参数集包括Δa=[δ1,δ2,δ3,δ4]和u0=[ur0,ut0,uh0]共七个参数。
所述步骤2具体实现如下:
依据步骤1中的方法可以得到航天器由初始轨道向目标轨道转移的控制律,在这一步中将以航天器的当前的位置和时刻为输入,建立发动机的开关机约束,只有满足开机条件下推力才会以步骤1给出的值输出,否则将置0,具体实现如下:
(1)由当前航天器轨道要素计算航天器在地心J2000坐标系中的位置和速度;
(2)根据地球自转模型和星历表计算得到当前时刻计算日、月和测控站在地心J2000坐标系中的位置;
(3)建立对日观测坐标系Oxhyhzh(Sh),将航天器位置和月球位置转换到Oxhyhzh(Sh)中,由柱形阴影区定义得到地影区和月影区的范围,得到地影区和月影区的关机约束,判断航天器是否处于地影区或月影区内;
(4)在地球固连坐标系Oxeyeze(Se)中以地球圆球模型计算给定观测站的坐标RO,e,得到地面测控站可见区开机约束,并将其转换到J2000坐标系中得到RO,i,计算航天器相对观测站的仰角,若仰角ε0>5°则航天器在观测站的可见区内;
(5)在J2000坐标系中根据航天器、观测站和太阳的位置得到日凌区关机约束,计算航天器-观测站-太阳的夹角,若夹角小于1.5°则航天器进入日凌区;
(6)将四个开关机约束作为一个判断条件P,若四个约束均满足,即航天器不在地影区和月影区内,地面测控站对航天器可见,不发生日凌现象,则P=1,否则P=0;
得到的带有复杂约束的轨道转移控制模型为:
所述步骤3具体实现如下:
依据步骤1和2中的方法可以得到带有复杂约束的轨道转移控制律,在这一步建立合适的性能函数,对轨道转移过程进行优化,具体实现如下:
(1)将可调参数集中的七个参数设为待优化函数的输入,构成一个7维状态变量S=[u0,Δa];
(3)根据每一个状态变量S,在初始轨道和目标轨道确定后,可以计算得到一个性能指标函数J,因此得到了一个状态变量S到性能指标函数J的映射,即:
J=F(S);
(4)智能优化算法包括遗传算法、粒子群算法、差分进化算法,从中选取一种算法,设置优化算法的迭代次数及关于状态变量S的搜索空间OS,得到带有复杂约束的轨道转移优化问题为:
min:J=F(S),S∈OS
式中min表示让性能指标函数J的值最小;
(5)利用智能优化算法求解带有复杂约束的轨道转移优化问题,得到使性能指标函数J最小的状态变量S*。
本发明提出的轨道转移控制及优化方法相比现有方法的优点在于:
(1)现有的轨道转移控制方法侧重于最优性分析,需要进行复杂的模型推导,且无法考虑转移过程中的开关机约束,轨道控制律求解算法的收敛性严重依赖于数学规划算法的性能,当初末轨道差异较大时任意出现算法无法收敛的情况,并且算法计算时间较长且难以控制。本发明根据航天器轨道转移的目标及工程约束条件,设计出一种仅依赖于航天器位置信息的轨道控制及优化方法,在考虑影响航天器发动机开关机的地影区、月影区、观测站通信以及日凌约束条件下,实现了连续推力轨道转移,并结合高性能的智能优化算法,进一步实现控制律优化,从而完成航天器转移轨迹优化。同时,本发明中设计的控制算法对于航天器的轨道动力学模型、约束条件模型、发动机推力配置等没有限制,可以根据实际的工程任务情况进行二次开发来细化和完善轨控模型,保证所生成的轨控策略符合各种任务场景中的约束条件。
(2)本发明以航天器的当前轨道与目标轨道的差异确定控制律,在转移过程中考虑多种开关机约束,得到满足多种开关机约束的轨道转移控制律。
(3)本发明轨道转移控制律受控制向量u0和邻域角向量Δa影响,给定轨道转移性能指标后,利用先进的智能优化算法可以对S=[u0,Δa]进行优化,得到更符合任务需求的轨道转移控制律,实现对满足多种开关机约束的轨道转移控制律的进一步优化。
(4)本发明中径向、切向和法向的推力输出只有{负,零,正}三种状态,适合于实际工程中实现。
(5)本发明中整个轨道转移控制及优化方法中的轨道动力学模型和约束模型与控制策略无关,因此允许将对应的计算模型替换为类似的其它模型,如航天器二体问题动力学模型可以增加高精度摄动力计算模块;日、月位置计算可以由高精度的JPL星历表模型改为简化计算的星历模型来提高程序计算效率;地球圆球模型替换为更精确的椭球模型;允许增加更多的开关机约束。在本发明所提出的轨道转移控制策略基础上,可以进一步根据实际工程任务的需要进行模型修改,使控制律在满足更多任务约束的条件下实现航天器轨道转移。
附图说明
图1为对日观测坐标系与J2000坐标系的关系;
图2为航天器变轨控制流程图;
图3为四个特征点邻域示意图;
图4为本发明的轨道转移控制及优化方法流程图;
图5为转移轨迹及相关状态变量时间历程,(a)算例1,(b)算例2,(c)算例3;
图6为燃料消耗优化曲线;(a)算例1,(b)算例2,(c)算例3。
具体实施方式
下面结合附图及实施例对本发明进行详细说明。
如图1所示,为了便于计算日、月、地、星的光路,定义对日观测坐标系Oxhyhzh(Sh),由地心J2000坐标系旋转得到。Oxhyhzh(Sh)的坐标原点为地心,xh轴由地心指向太阳,整个坐标系由J2000坐标系先绕zi轴旋转αh,再绕yi轴旋转-δh得到,其中αh和δh分别为太阳赤经和太阳赤纬。
将图2的变轨控制流程作为一个控制求解模块,进而可以将整套算法的求解流程在图4中给出。采用的智能优化算法可以选择当前研究中对函数最优解搜索效率高的现有算法,并根据具体的问题进行适当的改进,从而可以起到更好的优化效果。
如图2所示,本发明针对多发动机开关机约束下的航天器连续推力轨道转移问题,提出一种复杂约束下航天器连续推力轨道转移控制及优化方法。在轨道转移过程中,将发动机推力沿为用于调整轨道椭圆形状的密切轨道面内的径向和切向,以及用于调整轨道面方位的垂直于密切轨道面的法向,设计符合轨道动力学特性的控制策略,构建参数化轨道转移模型,并利用智能优化算法进行求解。同时考虑地影、月影、日凌及航天器与测控站通讯对发动机开关机的路径约束,使得轨道转移控制律更符合实际工程需求。
如图2所示,给定航天器轨道转移任务的初末轨道后,人为地设置好控制向量和邻域角向量地参数后,可以生成轨道控制律。控制律由航天器的实时位置与特征点邻域的关系确定,确定控制输出后还需要判断是否满足地影、月影、测控站可见、日凌等开关机约束,若不满足开机约束,则控制输出会被设置为0。确定控制状态后可以代入轨道动力学方程进行数值积分,得到下一时刻的航天器轨道位置,迭代计算直到航天器到达目标轨道或者判断轨道转移失败。航天器到达目标轨道后,可以得到航天器转移路线、转移代价等信息。
如图3所示,航天器的四个特征点是当前椭圆轨道上的点,其对应邻域为椭圆弧段,满足弧段内的点与特征点关于地心的夹角小于对应的邻域角。由图可知,特征点3和4所在的特征交线所在的直线即为当前轨道法向量和目标轨道法向量所在平面过地心的垂线。
如图4所示,确定轨道转移任务的信息后,将控制向量和邻域角向量作为待优化参数,利用智能优化算法对其进行调参,得到满足在“轨道转移任务代价评估”中最优的参数值,并得到其对应的轨道转移控制律,具体实现如下:
1.连续推力轨道动力学模型
在分析连续推力轨道转移问题时,转移时间通常很长,尤其是以电推进方式进行的小推力轨道转移,需要的飞行圈数可达几百圈。因此为了提高积分精度,往往采用轨道要素进行动力学分析。同时考虑到常用的分析场景为GEO转移问题,Kepler六要素在轨道倾角为0时与目标轨道为圆轨道时会出现奇点,因此采用改进春分点轨道要素X=[p,ex,ey,hx,hy,L]T作为状态量,其与Kepler轨道要素的关系为:
式中:p为椭圆轨道半通径,a为椭圆轨道半长轴,e为椭圆偏心率,i为轨道倾角,Ω为升交点赤经,ω为近地点幅角,θ为真近点角,μ为地球引力常数。
在地心J2000坐标系Oxiyizi(Si)中,航天器的动力学方程为:
式中:
其中的辅助变量定义在式(4)中。式(3)中的aacce为除了中心引力加速度外的加速度向量,[ar au ah]T为aacce以地心轨道坐标系(地心轨道坐标系Ox0y0z0原点在地球中心,x0轴由地心指向航天器,z0轴与航天器动量矩向量方向相同,y0轴与z0轴、x0轴满足右手准则,在轨道平面内与矢径方向相垂直)为基准的向量形式,aeng为航天器发动机推力加速度向量,adis为摄动力加速度向量。
将航天器质量设为m,Tmax为航天器发动机推力上限,比冲为常值Isp,重力加速度取地球表面赤道处有g=9.7804m/s2,则引入质量微分方程后的动力学模型为:
式中u为控制矢量,方向为发动机推力方向,幅值为输出推力与最大推力Tmax的比率。
2.航天器发动机开关机约束
实际航天任务中,考虑到航天器轨道环境、观测、通信与安全性需求,往往要求发动机只能在有光照、与地面站可通信的条件下才能开机,这就为轨道转移问题中的控制变量带来了复杂的路径约束,并且约束不仅与自身轨道位置相关,还与当前的时刻有关,使得理想中的最优控制解难以求解。
本发明中考虑的开关机约束包括地影、月影、日凌及航天器与测控站通讯四种,认为航天器发动机只能在非地影区、非月影区、与地面测控站可通信,且需避开“日凌中断通信”的条件下才允许开机。日、月在地心J2000坐标系中的位置向量RS,i和RM,i可由JPL的星历表DE430计算得到。
2.1地影区约束
当航天器运动到地球相对于太阳得背面时,无法接收到太阳光照射,即进入了地影区。常用得地影模型分为柱形地影模型和锥形地影模型,其中柱形模型认为日地距离远大于太阳和地球半径,可近似将到达地球的太阳光视为平行光。锥形模型则考虑了日地半径差与日地距离,对地影区进行更精确的划分,可分为本影区、半影区、光照区和伪影区。
本发明中按照柱形地影模型计算地影区,并且将地球视为半径为赤道半径rE=6378.137km的球体。由于地球是椭球,最大半径为赤道半径,因此这样计算得到的地影区比实际地影区保守。
为了便于计算日、月、地、星的光路,定义对日观测坐标系Oxhyhzh(Sh),由地心J2000坐标系旋转得到,关系如图1所示。Oxhyhzh(Sh)的坐标原点为地心,xh轴由地心指向太阳,整个坐标系由J2000坐标系先绕zi轴旋转αh,再绕yi轴旋转-δh得到,其中αh和δh分别为太阳赤经和太阳赤纬,可由RS,i计算得到[13]。由Oxiyizi(Si)到Oxhyhzh(Sh)的坐标转换矩阵为:
航天器在Oxiyizi(Si)中的坐标Ri=[xi,yi,zi]T可以根据当前时刻的改进春分点要素计算得到,则其在Oxhyhzh(Sh)中坐标可由Rh=[xh,yh,zh]T=LhiRi计算得到。当航天器进入地影区,在Oxhyhzh(Sh)中观察,显然满足式(7)条件。
月影区与地影区类似,为太阳光受到月球遮挡无法照射到航天器的区域。将月球视为半径rM=1738km的球体,同样在Oxhyhzh(Sh)中采用柱形模型来计算月影区约束。月心在Oxhyhzh(Sh)中坐标由RM,h=[xM,h,yM,h,zM,h]T=LhiRM,i计算,由于地星距离远小于地月距离,则航天器进入月影区的条件可由式(8)判断。
2.3观测站通信约束
为了保证航天器在轨运行的安全性,在任务中通常要求航天器必须处于观测站可见范围内才能进行轨道机动。观测站在地球固连坐标系Oxeyeze(Se)中的位置由地理纬度、经度及当地海拔确定,采用地球球模型计算观测站当前时刻在Oxeyeze(Se)中坐标RO,e [13],认为地球是匀速自转,可以得到Oxeyeze(Se)到Oxiyizi(Si)的转换矩阵Lie,则有RO,i=LieRO,e。为实现测控站对航天器的观测,需满足仰角ε0>5°,因此满足测控站通信约束的条件为:
地面测控网中包含N个测控站,则需至少有一个测控站可以对航天器进行观测。
2.4日凌约束
当航天器运动到接近测控站和太阳间连线的位置时,测控站与航天器间的通信会受到太阳噪声的严重干扰,即出现日凌现象,此时会要求航天器停止进行轨道机动。考虑到测控站天线的直径和工作波长,太阳的等视圆面,因此设不发生日凌的条件为航天器-观测站-太阳的夹角大于1.5°,即:
将四个开关机约束作为一个判断条件P,若四个约束均满足,即航天器不在地影区和月影区内,地面测控站对航天器可见,不发生日凌现象,则P=1,否则P=0;
3.轨道转移控制策略及优化问题模型
通过对式(2)、(3)的动力学模型分析,可知径向推力和切向推力只改变轨道面内椭圆的形状,轨道面的方位只能由法向推力控制,因此将初末轨道的差异分为轨道面内椭圆差异和轨道面方位差异进行分析。椭圆差异由控制形状的半长轴和偏心率,及控制拱线相位的近地点幅角体现。为了直观表现椭圆形状,用近地点和远地点的地星距更方便调整控制策略。
根据轨道动力学特性,切向推力比径向推力改变半长轴大小的效率更高,径向推力更适合改变拱线方向,法向推力在当前轨道面与目标轨道面的交线处的控制效率最高。依此特性,脉冲轨道在近地点施加切向脉冲调整远地点高度,在远地点调整近地点高度,在经过目标轨道面处施加法向脉冲调整轨道面。然而,连续推力无法瞬时改变速度矢量来调整轨道,因此在这些特征点的邻域内施加推力可以有更高的变轨效率。轨道转移的结果由推力的幅值与邻域大小的设置决定,所以给出控制策略如图2所示。
为了便于工程中实现,三个方向的推力输出只有{负,零,正}三种状态,即:
式中[ur0,ut0,uh0]定义为控制向量u0,需要在确定初末轨道信息后,开始变轨前确定,可以通过智能优化算法根据变轨需求进行参数优化。特征点的邻域由其对应的特征向量(由地心指向特征点)与航天器位置向量的邻域角确定,航天器轨道共有四个特征点,分别是其密切轨道的近地点、远地点、下穿越点和上穿越点,对应四个邻域{D1(ξ1,δ1),D2(ξ2,δ2),D3(ξ3,δ3),D4(ξ4,δ4)},其中ξ为特征向量,δ为邻域角,如图3所示。定义一个邻域角向量Δa=[δ1,δ2,δ3,δ4],同样可以通过智能优化算法根据变轨需求进行参数优化。下穿越点为航天器由上向下经过目标轨道面的位置,反之则是上穿越点。给定当前密切轨道法向量nosc和目标轨道法向量nT后,下穿越点的特征向量可由式(12)得到。由航天器和特征点邻域的位置关系可以得到如下控制策略:
式中,ustr为航天器发动机推力的方向矢量,Xc和XT分别为当前轨道和目标轨道的改进春分点根数,R为当前航天器位置,通过轨道要素转换方法可以由Xc计算得到,以函数R=g(Xc)描述。f(Xc,XT)=[ur,ut,uh]T中各控制分量的确定方法在表1中给出。结合轨道转移过程中的开关机约束,得到带有复杂开关机约束的轨道转移控制策略为:
径向推力适合于进行拱线调相,若目标轨道是圆轨道,则无需进行调相,只依靠切向推力和法向推力就能实现轨道转移。当目标轨道为一般的椭圆时,可以根据航天器当前密切轨道与目标轨道的近地点幅角差进行反馈控制,在D1(ξ1,δ1)和D2(ξ2,δ2)内控制切向推力,配合轨道面和轨道形状调整,实现航天器到目标轨道的转移。
表1变轨控制策略
*表1中ωc和ωT分别为当前轨道和目标轨道的近地点幅角
四个特征邻域中,D1和D2中进行面内椭圆进行调整,D3和D4中进行轨道面调整,两组邻域可能会出现重合。依靠式(11)中三个方向控制量的模均不大于保证了||u||≤1。但是这种各分量上界平均分配往往并不合适,当初末轨道面差异较小时,对uh0的上界可以适当减小;当面内椭圆形状差异较小时,对ut0的上界可以适当减小;当椭圆相位差异较小时,对ur0的上界可以适当减小。考虑到最优控制中Bang-Bang控制的特性,期望控制输出在0和最大值间切换,因此将依据控制策略给出控制量定义为ustr,实际的控制为
式(13)保证了发动机在四个特征邻域内开机时有||u||=1。
智能优化算法当前发展迅速,尤其是元启发式算法,通过引入随机搜索机制,极大地提升了对复杂优化模型的最优解的搜索效率。常用的智能优化算法包括遗传算法、粒子群算法、差分进化算法,在实际使用中选取合适的优化算法对优化问题求解。根据前面的轨道转移策略,设7维状态变量S=[u0,Δa],将其作为优化算法的输入参数,输出的性能指标函数则根据任务需求设置,对于时间最优问题,有对于燃料最优问题,有对于能量最优问题,有根据每一个状态变量S,在初始轨道和目标轨道确定后,可以计算得到一个性能指标函数J,因此得到了一个状态变量S到性能指标函数J的映射,即J=F(S)。基于现有的智能优化算法,设置优化算法的迭代次数及关于状态变量S的搜索空间OS,得到带有复杂约束的轨道转移优化问题为:
min:J=F(S),S∈OS
最终可以利用智能优化算法求解带有复杂约束的轨道转移优化问题,得到使性能指标函数J最小的状态变量S*。
4.仿真案例
针对大轨道倾角转移,椭圆形状差异明显及椭圆相位差异三种情况,设置三组仿真案例,初末轨道信息在表2中给出。其余仿真参数设置为:最大推力Tmax=10N,初始质量m0=1500kg,海平面重力加速度g=9.7804m/s2,比冲Isp=3600s,地球引力常数μ=3.986004415×1014m2/s3,不考虑摄动力影响,即adis=0,发射时刻设置为2021年7月1日12时。案例中采用差分进化算法进行参数优化,性能函数采用的是燃料最优问题,即定义一组标称参数为用于体现差分进化算法的参数优化效果,将此时三组案例中轨道转移燃料消耗值作为参考值。仿真结果在表3中给出。
表2初末轨道状态参数
*表中目标轨道参数的L无需给出
表3仿真结果
通过表3的结果可以看出,本发明所提出的轨道控制及优化可以实现多约束条件下的航天器连续推力轨道转移,并能利用智能优化算法对参数进行调节,实现推力控制的优化。
5.仿真结果
根据表3中的结果可以生成航天器轨道转移的轨迹如图5所示。观察可知,根据所设计的轨道控制算法,航天器能够实现在复杂约束下的轨迹转移,可满足实际工程中轨道转移任务的需求。三组算例中燃料消耗随迭代次数的变化曲线如图6所示。
图5给出了三个算例的轨道转移结果,包括转移轨迹、轨道要素和质量的时间曲线、控制输出的时间曲线。由图可知,本发明所提出的轨道控制策略可以实现给定初末轨道间的转移任务,所适用的任务类型范围广,包括三个算例,(a)算例1,(b)算例2,(c)算例3,算例1的降低轨道倾角,算例2的同时进行轨道倾角降低和轨道圆化,算例3的同时进行轨道倾角降低和椭圆拱线旋转。
图6给出了采用差分进化算法对轨道转移过程中的燃料消耗进行优化的结果,依靠控制向量和邻域角向量进行参数调节来实现。由图可知,在三个算例即(a)算例1,(b)算例2,(c)算例3中进行参数优化可以有效地改善轨道转移的代价消耗。
提供以上实施例仅仅是为了描述本发明的目的,而并非要限制本发明的范围。本发明的范围由所附权利要求限定。不脱离本发明的精神和原理而做出的各种等同替换和修改,均应涵盖在本发明的范围之内。
Claims (4)
1.一种复杂约束下航天器连续推力轨道转移控制及优化方法,其特征在于,包括以下步骤:
步骤1:基于航天器连续推力动力学模型,提出一种由航天器当前轨道特征点邻域和可调参数集确定的轨道转移控制策略,依靠常幅值径向、切向和法向分量的发动机推力实现轨道转移;
步骤2:基于步骤1中的轨道转移控制策略,依据航天器在实际任务中需要面对的开关机约束,建立多种航天器开关机约束条件,将轨道转移控制策略和多种航天器开关机约束条件结合,得到带有复杂约束的轨道转移控制模型;所述多种航天器开关机约束包括地影区内的关机约束、月影区内的关机约束、地面测控站可见区开机约束及日凌区关机约束;
步骤3:根据步骤2得到的带有复杂约束的轨道转移控制模型,将步骤1中的可调参数集构造为状态变量输入,建立控制性能指标函数作为输出,构建状态变量和控制性能指标间的函数映射,得到带有复杂约束的轨道转移优化问题;利用智能优化算法求解带有复杂约束的轨道转移优化问题,得到使性能指标函数值最小的状态变量,最终完成复杂约束下航天器连续推力轨道转移控制及优化工作。
2.根据权利要求1所述的复杂约束下航天器连续推力轨道转移控制及优化方法,其特征在于:所述步骤1中,轨道转移控制策略如下:
式中,ustr为航天器发动机推力的方向矢量,Xc和XT分别为当前轨道和目标轨道的改进春分点根数,R为当前航天器位置,通过轨道要素转换方法由Xc计算得到,以函数R=g(Xc)描述;
f(Xc,XT)=[ur,ut,uh]T,其中:
式中下标r、t和h分别表示径向、切向和法向,则ur、ut和uh分别是ustr沿径向、切向和法向的分量;定义控制参数向量u0=[ur0,ut0,uh0],在D1和D2中进行面内椭圆调整,只依靠径向和切向的推力,其中径向推力用于调整近地点幅角,切向推力用于调整轨道半长轴;在D3和D4中进行轨道面调整,只依靠法向推力;ur、ut和uh的取值方法具体给出如下:
(1)若R∈D1(ξ1,δ1),若远地点高度低,则flag1=1,否则flag1=-1;若ωc<ωT,则flag4=1,否则flag4=-1,有ut=flag1*ut0,ur=flag4*ur0,flag1和flag3为符号标记;
(2)若R∈D2(ξ2,δ2),若近地点高度低,则flag2=1,否则flag2=-1;若ωc<ωT,则flag4=1,否则flag4=-1,有ut=flag2*ut0,ur=flag4*ur0,flag2为符号标记;
(3)若R∈D3(ξ3,δ3),flag3=1,uh=flag3*uh0,flag3为符号标记;
(4)若R∈D4(ξ4,δ4),flag3=-1,uh=flag3*uh0;
其中,{D1(ξ1,δ1),D2(ξ2,δ2),D3(ξ3,δ3),D4(ξ4,δ4)}为航天器当前轨道椭圆上的四个弧段,分别由四个特征点和邻域角定义,四个特征点分别是当前轨道的近地点、远地点、下穿越点和上穿越点,下穿越点为航天器由上向下经过目标轨道面的位置,反之则是上穿越点,ξ为特征向量,由地心指向特征点;δ为邻域角,用于限制对应弧段上的点与特征点关于地心的夹角;定义一个邻域角向量Δa=[δ1,δ2,δ3,δ4],由于近地点和远地点连线,上穿越点和下穿越点连线均经过地心,因此有0°<δi<90°,i=1,2,3,4;可调参数集包括Δa=[δ1,δ2,δ3,δ4]和u0=[ur0,ut0,uh0]共七个参数,直接影响轨道转移的效果。
3.根据权利要求1所述的复杂约束下航天器连续推力轨道转移控制及优化方法,其特征在于:所述步骤2具体实现如下:
(1)由当前航天器轨道要素计算航天器在地心J2000坐标系中的位置和速度;
(2)根据地球自转模型和星历表计算得到当前时刻计算日、月和测控站在地心J2000坐标系中的位置;
(3)建立对日观测坐标系Oxhyhzh(Sh),将航天器位置和月球位置转换到Oxhyhzh(Sh)中,由柱形阴影区定义得到地影区和月影区的范围,得到地影区和月影区的关机约束,判断航天器是否处于地影区或月影区内;
(4)在地球固连坐标系Oxeyeze(Se)中以地球圆球模型计算给定观测站的坐标RO,e,得到地面测控站可见区开机约束,并将其转换到J2000坐标系中得到RO,i,计算航天器相对观测站的仰角,若仰角ε0>5°则航天器在观测站的可见区内;
(5)在J2000坐标系中根据航天器、观测站和太阳的位置得到日凌区关机约束,计算航天器-观测站-太阳的夹角,若夹角小于1.5°则航天器进入日凌区;
(6)将四个开关机约束作为一个判断条件P,若四个约束均满足,即航天器不在地影区和月影区内,地面测控站对航天器可见,不发生日凌现象,则P=1,否则P=0;
因此,带有复杂约束的轨道转移控制模型为:
4.根据权利要求1所述的复杂约束下航天器连续推力轨道转移控制及优化方法,其特征在于:所述步骤3具体实现如下:
(1)将可调参数集中的七个参数设为待优化函数的输入,构成一个7维状态变量S=[u0,Δa];
(3)根据每一个状态变量S,在初始轨道和目标轨道确定后,可以计算得到一个性能指标函数J,因此得到了一个状态变量S到性能指标函数J的映射,即:
J=F(S);
(4)智能优化算法包括遗传算法、粒子群算法、差分进化算法,从中选取一种算法,设置优化算法的迭代次数及关于状态变量S的搜索空间OS,得到带有复杂约束的轨道转移优化问题为:
min:J=F(S),S∈OS;
式中min表示让性能指标函数J的值最小;
(5)利用智能优化算法求解带有复杂约束的轨道转移优化问题,得到使性能指标函数J最小的状态变量S*。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202111473972.9A CN114169072B (zh) | 2021-12-02 | 2021-12-02 | 一种复杂约束下航天器连续推力轨道转移控制及优化方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202111473972.9A CN114169072B (zh) | 2021-12-02 | 2021-12-02 | 一种复杂约束下航天器连续推力轨道转移控制及优化方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN114169072A true CN114169072A (zh) | 2022-03-11 |
CN114169072B CN114169072B (zh) | 2024-06-11 |
Family
ID=80483275
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202111473972.9A Active CN114169072B (zh) | 2021-12-02 | 2021-12-02 | 一种复杂约束下航天器连续推力轨道转移控制及优化方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN114169072B (zh) |
Cited By (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN114701670A (zh) * | 2022-03-24 | 2022-07-05 | 北京航空航天大学 | 一种面向非合作目标的防碰撞悬停轨道控制方法和系统 |
CN114722711A (zh) * | 2022-04-11 | 2022-07-08 | 中国科学院空天信息创新研究院 | 一种受摄条件下的兰伯特转移轨道确定方法和系统 |
CN114955009A (zh) * | 2022-05-20 | 2022-08-30 | 北京航天飞行控制中心 | 一种地球轨道交会对接点的选取方法、系统、介质及设备 |
CN115196046A (zh) * | 2022-09-19 | 2022-10-18 | 航天东方红卫星有限公司 | 一种太阳同步轨道卫星超寿运行轨控策略确定方法 |
CN115314101A (zh) * | 2022-04-12 | 2022-11-08 | 中国人民解放军战略支援部队航天工程大学 | 一种基于并行计算的低轨通信卫星星座快速建模方法 |
Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN106202640A (zh) * | 2016-06-28 | 2016-12-07 | 西北工业大学 | 日‑地三体引力场中的晕轨道航天器偏置轨道设计方法 |
US20170297746A1 (en) * | 2015-11-20 | 2017-10-19 | Thales | Orbit transfer method for a spacecraft using a continuous or quasi-continuous thrust and embedded driving system for implementing such a method |
CN109375648A (zh) * | 2018-12-07 | 2019-02-22 | 北京理工大学 | 一种多约束条件下椭圆轨道卫星编队构形初始化方法 |
CN110276159A (zh) * | 2019-07-01 | 2019-09-24 | 北京理工大学 | 一种基于多模型融合的卫星系统多学科优化方法 |
-
2021
- 2021-12-02 CN CN202111473972.9A patent/CN114169072B/zh active Active
Patent Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20170297746A1 (en) * | 2015-11-20 | 2017-10-19 | Thales | Orbit transfer method for a spacecraft using a continuous or quasi-continuous thrust and embedded driving system for implementing such a method |
CN106202640A (zh) * | 2016-06-28 | 2016-12-07 | 西北工业大学 | 日‑地三体引力场中的晕轨道航天器偏置轨道设计方法 |
CN109375648A (zh) * | 2018-12-07 | 2019-02-22 | 北京理工大学 | 一种多约束条件下椭圆轨道卫星编队构形初始化方法 |
CN110276159A (zh) * | 2019-07-01 | 2019-09-24 | 北京理工大学 | 一种基于多模型融合的卫星系统多学科优化方法 |
Non-Patent Citations (1)
Title |
---|
孟雅哲: "航天器燃耗最优轨道直接/间接混合法延拓求解", 航空学报, vol. 38, no. 1, 25 January 2017 (2017-01-25), pages 264 - 285 * |
Cited By (8)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN114701670A (zh) * | 2022-03-24 | 2022-07-05 | 北京航空航天大学 | 一种面向非合作目标的防碰撞悬停轨道控制方法和系统 |
CN114722711A (zh) * | 2022-04-11 | 2022-07-08 | 中国科学院空天信息创新研究院 | 一种受摄条件下的兰伯特转移轨道确定方法和系统 |
CN114722711B (zh) * | 2022-04-11 | 2022-12-06 | 中国科学院空天信息创新研究院 | 一种受摄条件下的兰伯特转移轨道确定方法和系统 |
CN115314101A (zh) * | 2022-04-12 | 2022-11-08 | 中国人民解放军战略支援部队航天工程大学 | 一种基于并行计算的低轨通信卫星星座快速建模方法 |
CN115314101B (zh) * | 2022-04-12 | 2023-08-01 | 中国人民解放军战略支援部队航天工程大学 | 一种基于并行计算的低轨通信卫星星座快速建模方法 |
CN114955009A (zh) * | 2022-05-20 | 2022-08-30 | 北京航天飞行控制中心 | 一种地球轨道交会对接点的选取方法、系统、介质及设备 |
CN115196046A (zh) * | 2022-09-19 | 2022-10-18 | 航天东方红卫星有限公司 | 一种太阳同步轨道卫星超寿运行轨控策略确定方法 |
CN115196046B (zh) * | 2022-09-19 | 2022-12-13 | 航天东方红卫星有限公司 | 一种太阳同步轨道卫星超寿运行轨控策略确定方法 |
Also Published As
Publication number | Publication date |
---|---|
CN114169072B (zh) | 2024-06-11 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN114169072A (zh) | 一种复杂约束下航天器连续推力轨道转移控制及优化方法 | |
CN109592079A (zh) | 一种限定时间的航天器共面交会变轨策略确定方法 | |
CN110986974A (zh) | 面向复杂动力学环境的多航天器任务智能规划与控制方法 | |
CN105912819B (zh) | 一种地月l1拉格朗日点转移轨道的快速设计方法 | |
CN112016187B (zh) | 一种基于混合动力的近地小行星交会任务轨道优化方法 | |
Baig et al. | Light-levitated geostationary cylindrical orbits are feasible | |
Lücking et al. | Electrochromic orbit control for smart-dust devices | |
CN110031003B (zh) | 一种火箭上面级最优可达轨道快速规划计算方法 | |
RU2608186C2 (ru) | Способ и система для управления группой, по меньшей мере, из двух спутников, выполненных с возможностью обеспечения обслуживания | |
CN109240340B (zh) | 一种基于拟周期轨道的洛伦兹力多星编队构型方法 | |
CN112629543A (zh) | 一种大椭圆轨道及小倾角圆轨道的轨道规划方法 | |
CN115373423B (zh) | 一种用于商业卫星的编队捕获方法 | |
CN112572835A (zh) | 一种具有姿态切换的卫星在轨角动量管理及控制方法 | |
CN110723315B (zh) | 一种天体表面弹道式飞行探测的轨迹生成方法 | |
CN115258196A (zh) | 低轨卫星星座组网电推进变轨策略优化方法和系统 | |
CN106446313B (zh) | 基于极地悬停卫星轨道的系统设计方法 | |
CN110096746A (zh) | 一种卫星集群初始轨道设计方法及装置 | |
CN113093246A (zh) | 地面多目标点成像快速判定及任务参数计算方法 | |
Inalhan et al. | Precise formation flying control of multiple spacecraft using carrier-phase differential GPS | |
Shirazi | Multi-objective optimization of orbit transfer trajectory using imperialist competitive algorithm | |
Mingjian et al. | Flight trajectory optimization of sun-tracking solar aircraft under the constraint of mission region | |
Kolawole et al. | Backtracking search algorithm for non-aligned thrust optimization for satellite formation | |
CN115892519A (zh) | 一种用于近距离航天器轨道脉冲博弈的航天器控制方法 | |
Petropoulos | Some analytic integrals of the averaged variational equations for a thrusting spacecraft | |
CN114384803B (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 |