CN114154253B - 考虑发动机关机和参数强非线性的连续推力轨迹优化方法 - Google Patents

考虑发动机关机和参数强非线性的连续推力轨迹优化方法 Download PDF

Info

Publication number
CN114154253B
CN114154253B CN202210123114.XA CN202210123114A CN114154253B CN 114154253 B CN114154253 B CN 114154253B CN 202210123114 A CN202210123114 A CN 202210123114A CN 114154253 B CN114154253 B CN 114154253B
Authority
CN
China
Prior art keywords
thrust
engine
shutdown
track
continuous
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
CN202210123114.XA
Other languages
English (en)
Other versions
CN114154253A (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.)
Beijing Institute of Technology BIT
Original Assignee
Beijing Institute of Technology BIT
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 Beijing Institute of Technology BIT filed Critical Beijing Institute of Technology BIT
Priority to CN202210123114.XA priority Critical patent/CN114154253B/zh
Publication of CN114154253A publication Critical patent/CN114154253A/zh
Application granted granted Critical
Publication of CN114154253B publication Critical patent/CN114154253B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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
    • 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
    • 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)
  • Computer Hardware Design (AREA)
  • General Engineering & Computer Science (AREA)
  • Evolutionary Computation (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • Automation & Control Theory (AREA)
  • Combined Controls Of Internal Combustion Engines (AREA)

Abstract

本发明公开的考虑发动机关机和参数强非线性的连续推力轨迹优化方法，属于航空航天技术领域。实现方法为：首先根据开关机约束和发动机类型给出强非线性推力和比冲变化模型；然后建立基于改进春分点根数考虑发动机关机和参数强非线性的转移轨迹动力学模型；之后根据任务需求，给出考虑发动机关机和参数强非线性的连续推力轨迹优化问题的约束和性能指标；然后通过线性化动力学和松弛非线性等式约束，将考虑发动机关机和参数强非线性的连续推力轨迹优化问题凸化；然后通过数值积分和逐次逼近快速求解考虑关机约束的变参数连续推力轨迹优化问题，得到最优转移轨迹及对应的最优控制方向，根据优化结果解决推力控制领域相关问题。本发明对近地轨道连续推力转移和深空轨迹连续推力转移均适用，适用范围广且鲁棒性强。

Description

考虑发动机关机和参数强非线性的连续推力轨迹优化方法
技术领域
本发明涉及一种考虑发动机关机和参数强非线性的连续推力轨迹优化方法，尤其适用于携带可变比冲和推力发动机探测器，同时具有定期关机约束的连续推力轨迹快速优化计算，属于航空航天技术领域。
背景技术
相比于传统的化学推进发动机，由于连续推力发动机的高比冲，探测器可以携带重量更少的燃料实现相同的速度增量，从而获得更加丰富的太空探测任务回报。日益复杂的探测任务对机载自主能力需求也日益增加，自主轨迹规划能力尤为关注。然而，对于搭载连续推力发动机的探测器而言，实现快速的在线轨迹规划具有很强的挑战性。除此之外，考虑发动机强非线性参数变化以及定期关机约束相比传统的常推力且无定期关机的情形更加符合实际任务需求，如果不考虑关机需求会导致无法准确对探测器进行定轨，末端误差大，而传统连续推力轨迹优化方法难以解决考虑上述两种非线性约束的轨迹优化问题。在已发展的关于连续推力轨迹优化先技术[1]（参见 Jiang F, Baoyin H, Li J. Practicaltechniques for low-thrust trajectory optimization with homotopic approach[J].Journal of Guidance, Control, and Dynamics, 2012, 35(1): 245-258.）提出基于同伦法的常推力连续推力转移轨迹优化方法，虽然该方法极大地提升了问题的收敛性，但是需要从多组不同协态变量初值的结果中筛选最优转移轨迹，因此其计算效率一般。
发明内容
为了解决不考虑发动机关机和参数强非线性导致的实际工程中无法准确定轨以及变参数发动机无法采用传统轨迹优化方法进行优化的问题，本发明的主要目的是提供一种考虑发动机关机和参数强非线性的连续推力轨迹优化方法，通过逐次凸规划实现考虑测控等限制需定期关机和携带类似于推力、比冲随日心距变化的离子发动机时的连续推力轨迹优化问题的高效求解，实现考虑发动机关机和参数强非线性的连续推力轨迹优化。
本发明的目的是通过下述技术方案实现的：
本发明公开的考虑发动机关机和参数强非线性的连续推力轨迹优化方法，首先根据开关机约束和发动机类型给出强非线性推力和比冲变化模型；然后建立基于改进春分点根数考虑发动机关机和参数强非线性的转移轨迹动力学模型；之后根据实际任务需求，给出考虑发动机关机和参数强非线性的连续推力轨迹优化问题的约束和性能指标；然后通过线性化动力学和松弛非线性等式约束，将考虑发动机关机和参数强非线性的连续推力轨迹优化问题凸化；然后通过数值积分和逐次逼近快速求解考虑关机约束的变参数连续推力轨迹优化问题，得到最优转移轨迹及对应的最优控制方向，即实现考虑发动机关机和参数强非线性的连续推力轨迹优化。
本发明公开的考虑发动机关机和参数强非线性的连续推力轨迹优化方法，包括如下步骤：
步骤一：根据关机约束和发动机类型给出强非线性推力和比冲变化模型。
强非线性参数发动机的推力和比冲随探测器帆板的输入功率变化，而帆板输入功率与日心距成反比：
Figure DEST_PATH_IMAGE001
（1）
其中r AU 是日心距离，单位为 AU，P 0为日心距为1AU时帆板的输入功率，P in 是日心距为r AU 时帆板输入功率；发动机最大推力T max和比冲随帆板输入功率的变化有两种模式，一种为多项式变化，一种是阶梯变化；将最大推力T max和比冲I sp 表示为输入功率的函数G和H，即
Figure DEST_PATH_IMAGE002
（2）
探测器进行转移时，需要定期对探测器进行定轨，此时发动机需定期关闭；考虑发动机开关机约束用标志符ut表示为
Figure DEST_PATH_IMAGE003
（3）
其中，t为飞行时间，t on 和t off 分别表示发动机的最大连续开机时长和执行定轨所需的单次最短关机时长；进而考虑关机约束的发动机最大推力T max为
Figure DEST_PATH_IMAGE004
(4)
步骤二：建立基于改进春分点根数考虑发动机关机和参数强非线性的转移轨迹动力学模型。
采用改进春分点根数
Figure DEST_PATH_IMAGE005
描述探测器的状态，其中
Figure DEST_PATH_IMAGE006
，
Figure DEST_PATH_IMAGE007
，
Figure DEST_PATH_IMAGE008
，
Figure DEST_PATH_IMAGE009
，
Figure DEST_PATH_IMAGE010
，
Figure DEST_PATH_IMAGE011
，a、e、i、ω、Ω和υ分别为探测器轨迹的半长轴、偏心 率、轨道倾角、近地点幅角、升交点赤经和真近点角；则连续推力探测器的轨道动力学模型 为
Figure DEST_PATH_IMAGE012
(5)
其中
Figure DEST_PATH_IMAGE013
(6)
且控制向量的系数矩阵为：
Figure DEST_PATH_IMAGE014
(7)
Figure DEST_PATH_IMAGE015
，
Figure DEST_PATH_IMAGE016
，
Figure DEST_PATH_IMAGE017
其中，I sp (r)和T(r)是与日心距相关的发动机强非线性比冲和推力大小；
Figure DEST_PATH_IMAGE018
=1 + f cos L + g sin L, s 2 = 1 + h 2 + k 2，μ是中心天体的引力常数，m是探测器质量，C为控制 向量，C=[fr,ft,fn,T]T，fr,ft,fn分别发动机径向、切向、法向推力大小。
步骤三：根据实际任务需求，给出考虑发动机关机和参数强非线性的连续推力轨迹优化问题的约束和性能指标。始末端探测器状态均由任务时间对应的星历给出。相应的改进春分点轨道根数为：
Figure DEST_PATH_IMAGE019
(8)
Figure DEST_PATH_IMAGE020
(9)
其中t 0 与t f 分别为始末时间；同时，初始质量m(t 0 )固定，末端质量m(t f )无约束；探测器的推力分量满足
Figure DEST_PATH_IMAGE021
(10)
且
Figure DEST_PATH_IMAGE022
(11)
期望得到燃料最优连续推力轨迹，故将优化问题的性能指标设定为：
Figure DEST_PATH_IMAGE023
(12)
目的为最大化探测器末端质量，即最小化燃料消耗。
综上所述，燃料最优连续推力转移轨迹优化问题P1总结为：
Figure DEST_PATH_IMAGE024
约束方程：式(8) ~ (11)。
步骤四：通过线性化动力学和松弛非线性等式约束，将考虑发动机关机和参数强非线性的连续推力轨迹优化问题凸化。
为了将探测器动力学凸化，将非线性动力学方程(5)基于小扰动的连续线性化方法近似。连续近似过程中第k次迭代存在一个解X k ；然后，第k+1次迭代过程中，在X k 存在的前提下将动力学方程线性化，主项H(X)在Xk附近线性化；故线性化后的动力学方程为：
Figure DEST_PATH_IMAGE025
(13)
其中
Figure DEST_PATH_IMAGE026
；状态向量系数矩阵为：
Figure DEST_PATH_IMAGE027
(14)
其中，
Figure DEST_PATH_IMAGE028
，
Figure DEST_PATH_IMAGE029
，
Figure DEST_PATH_IMAGE030
，
Figure DEST_PATH_IMAGE031
；
方程(10)中推力向量的非线性约束函数是非凸的；在运用凸优化方法前必须将其化为凸函数；故将式(10)中等号松弛为不等号，然后该约束就转化为凸约束，即：
Figure DEST_PATH_IMAGE032
(15)
从而，线性化后的动力学(13)、始末约束(8)(9)、凸约束(15)、性能指标(12)一起构成了凸子问题P2。
步骤五：通过数值积分和逐次逼近快速求解考虑关机约束的变参数连续推力轨迹优化问题，得到最优转移轨迹及对应的最优控制方向，即实现考虑发动机关机和参数强非线性的连续推力轨迹优化。采用梯形法对式(13)中的数值积分进行转化，从而将问题 P2转化为凸优化问题形式。将转移任务时间区间[t 0 ，t f ]划分为N+1个节点，第i个节点处探测器的状态量和控制量记为X i 、U i ，则变量集为[X 0, ..., X i, ..., X N ]，待求解控制集为[C 0, ..., C i , , ..., C N ]，则动力学积分(13)转化为
Figure DEST_PATH_IMAGE034
(16)
其中i = 0, 1, …, N-1，
Figure DEST_PATH_IMAGE036
Figure DEST_PATH_IMAGE038
Figure DEST_PATH_IMAGE040
。
从而，凸子问题P2转化为凸问题。
对凸问题进行连续逼近以便迭代求解得到最优的考虑关机约束的变参数连续推力转移轨迹，直至其解收敛于P1的解。
令k=0，给出初始状态向量的猜测值X0，由从初值X(t 0 )到终值X(t f )的直线获取；L(t)的猜测值选取从L(t 0 )到L(t f )＋2∏·Γ的一条直线，其中Γ为连续推力转移轨迹运行圈数的估计值；m(t)的猜测值选取m(t 0 )到m(t f )的一条直线。
对于第i+1次迭代，选取第k次迭代的解X k 作为状态向量的猜测初值，解对为
Figure DEST_PATH_IMAGE041
。
检查是否满足收敛条件：
Figure DEST_PATH_IMAGE042
(17)
其中，ρ为收敛精度要求；如果不满足式(17)，需要继续迭代求解，如果满足式(17)即得到问题P1的解X*＝X k+1、C*＝C k+1。
至此，得到考虑发动机定期关机约束和强非线性变化参数的连续推力轨迹优化问题的解，X*为最优转移轨迹，C*为对应的最优控制方向。
还包括步骤六：根据步骤一至步骤五优化得到的最优转移轨迹及对应的最优控制方向， 解决推力控制领域相关问题。
有益效果：
（1）本发明的一种考虑发动机关机和参数强非线性的连续推力轨迹优化方法，由于对发动机的开关机和强非线性参数进行建模并基于此进行连续推力轨迹优化问题的处理和求解，所以对实际工程中考虑定轨约束需定期关机以及发动机比冲和推力大小强非线性变化的连续推力场景适用性强，末端状态精度高。
（2）本发明的一种考虑发动机关机和参数强非线性的连续推力轨迹优化方法，通过将非线性连续推力转移轨迹优化问题转化为凸问题，并应用序列凸规划方法进行推力和比冲可变且定期关机的连续推力轨迹优化，计算效率高，通常只需要秒级运算，相较于应用较为广泛的同伦法效率显著增加。
（3）本发明的一种考虑发动机关机和参数强非线性的连续推力轨迹优化方法，优化求解过程具有普适性，因此对近地轨道和深空探测轨迹考虑发动机关机和强参数非线性的连续推力转移均适用，适用范围广。
（4）本发明的一种考虑发动机关机和参数强非线性的连续推力轨迹优化方法，相较于传统的连续推力轨迹优化初值猜测效率较低的问题，本发明初值猜测只需对始末状态线性离散作为猜测初值，因此鲁棒性强。
附图说明
图 1 为本发明公开的一种考虑发动机关机和参数强非线性的连续推力轨迹优化方法流程图；
图 2 为本发明公开的一种考虑发动机关机和参数强非线性的连续推力轨迹优化方法求解实施例1得到的最优连续推力转移轨迹，图2(a)为 XY 平面轨迹，图 2(b)为YZ平面轨迹；
图 3 为本发明公开的一种考虑发动机关机和参数强非线性的连续推力轨迹优化方法求解实施例1得到的最优推力大小变化曲线；
图 4 为本发明公开的一种考虑发动机关机和参数强非线性的连续推力轨迹优化方法求解实施例1得到的最优径向、切向、法向推力分量变化曲线；
图 5 为本发明公开的一种考虑发动机关机和参数强非线性的连续推力轨迹优化方法求解实施例2得到的最优连续推力转移轨迹；
图 6 为本发明公开的一种考虑发动机关机和参数强非线性的连续推力轨迹优化方法求解实施例2得到的最优推力大小变化曲线；
图 7 为本发明公开的一种考虑发动机关机和参数强非线性的连续推力轨迹优化方法求解实施例2得到的最优径向、切向、法向推力分量变化曲线。
具体实施方式
为了更好地说明本发明的目的和优点，下面结合具体实施示例对本发明做出详细解释。
实施例1：
选择目标天体为彗星67P，设计从地球至彗星67P的连续推力转移轨迹。探测器初 始质量m(t 0 )为1217.7 kg。考虑发动机定期关机约束为最大连续开机时长 ton = 14天，最 小单次关机时长toff = 1天。地球出发日期为2021年10月28日 0:0:0，探测器与彗星67P交 会日期为2028年10月28日0:0:0。根据星历可得地球出发时和与67P交会时探测器状态用改 进春分点根数
Figure 929605DEST_PATH_IMAGE005
表示如表1所示。
表1地球-67P转移初始、终端状态
改进春分点根数 P，m f g h k L，rad
初始值 1.495568e+11 -0.0038 0.0163 -2.4755e-5 2.2439e-6 13.1684
末端值 3.051940e+11 0.2917 0.5707 0.0394 0.0472 2.5061
如图1所示，本实施例提供一种考虑发动机关机和参数强非线性的连续推力轨迹优化方法，具体实现步骤如下：
步骤一：根据关机约束和发动机类型给出强非线性推力和比冲变化模型。 强非线性参数发动机的推力和比冲随探测器帆板的输入功率变化，而帆板输入功率一般与日心距成反比：
Figure 201186DEST_PATH_IMAGE001
(18)
其中r AU 是日心距离，单位为 AU，P 0为日心距为1AU时帆板的输入功率，P in 是日心距为r AU 时帆板输入功率；发动机最大推力T max和比冲随帆板输入功率的变化有两种模式，一种为多项式变化，一种是阶梯变化；将最大推力T max和比冲I sp 表示为输入功率的函数G和H，即
Figure 229185DEST_PATH_IMAGE002
(19)
探测器进行转移时，需要定期对探测器进行定轨，此时发动机需定期关闭；考虑发动机开关机约束用标志符ut表示为
Figure DEST_PATH_IMAGE043
(20)
其中，t为飞行时间，t on 和t off 分别表示发动机的最大连续开机时长和执行定轨所需的单次最短关机时长；进而考虑关机约束的发动机最大推力T max为
Figure 53178DEST_PATH_IMAGE004
(21)
采用如表 2 所示的发动机参数阶梯变化模型，每个输入功率级别对应一组最大推力和比冲。
表2 发动机参数阶梯变化表
参数 数值 数值 数值 数值 数值 数值 数值 数值 数值 数值 数值 数值
输入功率（kW） 4.64 4.16 3.34 2.74 2.32 2.08 1.86 1.67 1.37 1.14 0.93 0.57
最大推力（mN） 184.6 166.6 132.2 105.6 92.3 83.3 69.6 66.1 52.8 43.6 34.8 21.8
比冲（s） 3313 3293 3291 3300 3313 3293 2974 3291 3300 2188 2974 2188
选取探测器初始质量m(t 0 )=1217.7 kg。考虑发动机定期关机约束为最大连续开机时长ton = 14天，最小单次关机时长toff = 1天。地球出发日期为2021年10月28日0:0:0，探测器与彗星67P 交会日期为2028年10月28日0:0:0。
步骤二：建立基于改进春分点根数考虑发动机关机和参数强非线性的转移轨迹动力学模型。
采用改进春分点根数
Figure 898643DEST_PATH_IMAGE005
描述探测器的状态，其中
Figure 454259DEST_PATH_IMAGE006
，
Figure 407652DEST_PATH_IMAGE007
，
Figure 365113DEST_PATH_IMAGE008
，
Figure 725687DEST_PATH_IMAGE009
，
Figure 768598DEST_PATH_IMAGE010
，
Figure 797122DEST_PATH_IMAGE011
Figure DEST_PATH_IMAGE044
,a、e、i、ω、Ω和υ分别为探测器轨迹的 半长轴、偏心率、轨道倾角、近地点幅角、升交点赤经和真近点角；则连续推力探测器的轨道 动力学模型为
Figure DEST_PATH_IMAGE045
(22)
其中
Figure DEST_PATH_IMAGE046
(23)
且控制向量的系数矩阵为：
Figure DEST_PATH_IMAGE047
(24)
Figure DEST_PATH_IMAGE048
，
Figure DEST_PATH_IMAGE049
，
Figure 197575DEST_PATH_IMAGE017
其中，I sp (r)和T(r)是与日心距相关的发动机强非线性比冲和推力大小；
Figure DEST_PATH_IMAGE050
=1 + f cos L + g sin L, s 2 = 1 + h 2 + k 2，μ是中心天体的引力常数，m是探测器质量，C为控制 向量，C=[fr,ft,fn,T]T，fr,ft,fn分别发动机径向、切向、法向推力大小。
步骤三：根据实际任务需求，给出考虑发动机关机和参数强非线性的连续推力轨迹优化问题的约束和性能指标。
在本任务场景中，始末端探测器状态均由任务时间对应的地球和彗星67P星历给出。相应的改进春分点轨道根数为：
X(t0)＝[1.495568e＋11, －0.0038, 0.0163, －2.4755e－5, 2.2439e－6,13.1684] (25)
X(tf)＝[3.051940e＋11, 0.2917, 0.5707, 0.0394, 0.0472, 2.5061] (26)
其中t0与tf分别为始末时间。同时，初始质量m(t0)=1217.7kg，末端质量m(tf)无约束。探测器的推力分量满足
满足
Figure DEST_PATH_IMAGE051
(27)
且
Figure 964936DEST_PATH_IMAGE051
(28)
期望得到燃料最优连续推力轨迹，故将优化问题的性能指标设定为：
Figure DEST_PATH_IMAGE052
(29)
目的为最大化探测器末端质量，即最小化燃料消耗。
综上所述，燃料最优小推力转移轨迹优化问题（也被称为P1）可总结为：
Figure DEST_PATH_IMAGE053
约束方程：式(25) ~ (28)。
步骤四：通过线性化动力学和松弛非线性等式约束，将考虑发动机关机和参数强非线性的连续推力轨迹优化问题凸化。
为了将探测器动力学凸化，将非线性动力学方程(22)基于小扰动的连续线性化方法近似。连续近似过程中第k次迭代存在一个解X k ；然后，第k+1次迭代过程中，在X k 存在的前提下将动力学方程线性化，主项H(X)在Xk附近线性化；故线性化后的动力学方程为：
Figure DEST_PATH_IMAGE054
(30)
其中
Figure DEST_PATH_IMAGE055
；状态向量系数矩阵为：
Figure DEST_PATH_IMAGE057
(31)
其中，
Figure 692546DEST_PATH_IMAGE028
，
Figure 521831DEST_PATH_IMAGE029
，
Figure 532512DEST_PATH_IMAGE030
，
Figure 952998DEST_PATH_IMAGE031
；
方程(27)中推力向量的非线性约束函数是非凸的；在运用凸优化方法前必须将其化为凸函数；故将式(27)中等号松弛为不等号，然后该约束就转化为凸约束，即：
Figure DEST_PATH_IMAGE058
(32)
从而，线性化后的动力学(27)、始末约束(25)(26)、凸约束(32)、性能指标(29)一起构成了凸子问题P2。
步骤五：通过数值积分和逐次逼近快速求解考虑关机约束的变参数连续推力轨迹优化问题。采用梯形法对式(30)中的数值积分进行转化，从而将问题P2转化为凸优化问题形式。将转移任务时间区间[t 0 ，t f ]划分为N+1个节点，第i个节点处探测器的状态量和控制量记为X i 、U i ，则变量集为[X 0, ..., X i, ..., X N ]，待求解控制集为[C 0, ..., C i , , ..., C N ]，则动力学积分(30)转化为
Figure DEST_PATH_IMAGE059
(33)
其中i = 0, 1, …, N-1，
Figure DEST_PATH_IMAGE036A
Figure DEST_PATH_IMAGE060
Figure DEST_PATH_IMAGE061
，
从而，凸子问题P2转化为凸问题。
对凸问题进行连续逼近以便迭代求解得到最优的考虑关机约束的变参数连续推力转移轨迹，直至其解收敛于P1的解。
令k=0，给出初始状态向量的猜测值X0，由从初值X(t 0 )到终值X(t f )的直线获取；L(t)的猜测值选取从L(t 0 )到L(t f )＋2∏·Γ的一条直线，其中Γ为连续推力转移轨迹运行圈数的估计值；m(t)的猜测值选取m(t 0 )到m(t f )的一条直线。
i+1次迭代，选取第k次迭代的解次迭代，选取第k次迭代的解X k 作为状态向量的猜 测初值，解对为
Figure DEST_PATH_IMAGE062
。
检查是否满足收敛条件：
Figure 794003DEST_PATH_IMAGE042
(34)
其中，ρ为收敛精度要求；如果不满足式(34)，需要继续迭代求解，如果满足式(34)即得到问题P1的解X*＝X k+1、C*＝C k+1，即得到考虑发动机定期关机约束和强非线性变化参数的连续推力轨迹优化问题的解，X*为最优转移轨迹，C*为对应的最优控制方向。
经过优化，在16次迭代后满足收敛条件(34)，得到探测器在转移过程中燃料消耗为337.5kg，考虑定轨约束发动机定期关机时能够准确到达末端约束处的状态。最优连续推力转 移轨迹如图2所示，对应的最优推力幅值变化曲线如图3所示，其中阶梯状变化幅值与发动机推力特性相符，每条起于X轴的纵向竖线对应一次发动机开/关机，最优径向、切向、法向推力分量变化曲线如图4所示。
步骤六：根据步骤一至步骤五优化得到的最优转移轨迹及对应的最优控制方向，对实际工程中考虑定轨约束需定期关机以及发动机比冲和推力大小强非线性变化的连续推力场景适用性强，末端状态精度高。
实施例2：
选择目标天体为小行星2011UW158，设计从地球至2011UW158 的连续推力转移轨 迹。 探测器初始质量m(t 0 )为1217.7 kg。考虑发动机定期关机约束为最大连续开机时长ton = 7 天，最小单次关机时长toff = 1 天。地球出发日期为2021年10月28日0:0:0，探测器与 彗星67P交会日期为2026年10月28日0:0:0。根据星历可得地球出发时和与小行星 2011UW158交会时探测器状态用改进春分点根数
Figure DEST_PATH_IMAGE063
表示如表1所示。
表1 地球-2011UW158转移初始、终端状态
改进春分点根数 P，m f g h k L，rad
初始值 1.495568e+11 -0.0038 0.0163 -2.4755e-5 2.2439e-6 13.1684
末端值 2.081291e+11 0.1576 -0.3416 0.0110 -0.0384 8.1845
如图1所示，本实施例提供一种考虑发动机关机和参数强非线性的连续推力轨迹优化方法，具体实现步骤如下：
步骤一：根据关机约束和发动机类型给出强非线性推力和比冲变化模型。
强非线性参数发动机的推力和比冲随探测器帆板的输入功率变化，而帆板输入功率一般 与日心距成反比：
Figure 161399DEST_PATH_IMAGE001
(35)
其中r AU 是日心距离，单位为 AU，P 0为日心距为1AU时帆板的输入功率，P in 是日心距为r AU 时帆板输入功率；发动机最大推力T max和I sp 比冲随帆板输入功率的变化有两种模式，一种为多项式变化，一种是阶梯变化。可以将其统一表示为输入功率的函数G和H，即
Figure 416800DEST_PATH_IMAGE002
(36)
探测器进行转移时，需要定期对探测器进行定轨，此时发动机需定期关闭。考虑发动机开关机约束用标志符ut表示为
Figure 214379DEST_PATH_IMAGE043
(37)
其中，t为飞行时间，t on 和t off 分别表示发动机的最大连续开机时长和执行定轨所需的单次最短关机时长。进而考虑关机约束的发动机最大推力T max为
Figure 453600DEST_PATH_IMAGE004
(38)
采用如式(38)所示的发动机参数连续变化模型，发动机最大推力T max和燃料质量 秒流量
Figure DEST_PATH_IMAGE064
为
Figure DEST_PATH_IMAGE066
（39）
其中，g0=9.8m/s为海平面重力加速度，a0, a1, a2, a3, a4, b0, b1, b2, b3, b4为变参数发动机模型系数，其数值分别为
a0=-0.1974, a1=15.24614, a2=2.5679, a3=-0.7898, a4=0.0532,
b0=438.1, b1=819.7, b2=146.9, b3=-63.2, b4=4.8,
选取探测器初始质量m(t 0 )为1217.7kg。选取探测器初始质量ton=7天，最小单次关机时长toff = 1天。地球出发日期为2021年10月28日0:0:0，探测器与小行星 2011UW158交会日期为2026年10月28日0:0:0。
步骤二：建立基于改进春分点根数考虑发动机关机和参数强非线性的转移轨迹动力学模型。
采用改进春分点根数
Figure DEST_PATH_IMAGE067
描述探测器的状态，其中
Figure DEST_PATH_IMAGE068
,
Figure DEST_PATH_IMAGE069
,
Figure DEST_PATH_IMAGE070
,
Figure DEST_PATH_IMAGE071
，
Figure DEST_PATH_IMAGE072
，
Figure DEST_PATH_IMAGE073
，a、e、i、ω、Ω和υ分别为探测器轨迹的半长轴、偏心率、轨道倾角、近地点幅角、升交点赤经 和真近点角；则连续推力探测器的轨道动力学模型为
Figure DEST_PATH_IMAGE074
(40)
其中
Figure DEST_PATH_IMAGE075
(41)
且控制向量的系数矩阵为：
Figure 244926DEST_PATH_IMAGE047
(42)
Figure 230200DEST_PATH_IMAGE015
，
Figure 195750DEST_PATH_IMAGE016
，
Figure 922267DEST_PATH_IMAGE017
其中，I sp (r)和T(r)是与日心距相关的发动机强非线性比冲和推力大小；
Figure 362957DEST_PATH_IMAGE050
=1 + f cos L + g sin L, s 2 = 1 + h 2 + k 2，μ是中心天体的引力常数，m是探测器质量，C为控制 向量，C=[fr,ft,fn,T]T，fr,ft,fn分别发动机径向、切向、法向推力大小。
步骤三：根据实际任务需求，给出考虑发动机关机和参数强非线性的连续推力轨迹优化 问题的约束和性能指标。
X（t0）=[1.495568e+11, -0.0038, 0.0163, -2.4755e-5, 2.2439e-6, 13.1684] (43)
X（tf）=[2.081291e 11, 0.1576, 0.3416, 0.0110, -0.0384, 8.1845] (44)
其中t 0 与t f 分别为始末时间。同时，初始质量m(t0) = 1217.7kg，末端质量m(t f )无约束。探测器的推力分量满足
满足
Figure DEST_PATH_IMAGE076
(45)
且
Figure DEST_PATH_IMAGE077
（46）
期望得到燃料最优连续推力轨迹，故将优化问题的性能指标设定为：
Figure DEST_PATH_IMAGE078
（47）
目的为最大化探测器末端质量，即最小化燃料消耗。
综上所述，燃料最优小推力转移轨迹优化问题（也被称为P1）可总结为：
Figure DEST_PATH_IMAGE079
约束方程：式(43) ~ (46)。
步骤四：通过线性化动力学和松弛非线性等式约束，将连续推力轨迹优化问题凸化。
为了将探测器动力学凸化，将非线性动力学方程(40)基小扰动的连续线性化方法近似。连续近似过程中第k次迭代存在一个解X k ；然后，第k+1次迭代过程中，在X k 存在的前提下将动力学方程线性化，主项H(X)在Xk附近线性化；故线性化后的动力学方程为：
Figure DEST_PATH_IMAGE080
(48)
其中
Figure DEST_PATH_IMAGE081
；状态向量系数矩阵为：
Figure DEST_PATH_IMAGE082
(49)
其中，
Figure 105187DEST_PATH_IMAGE028
，
Figure 976060DEST_PATH_IMAGE029
，
Figure 268501DEST_PATH_IMAGE030
，
Figure 784321DEST_PATH_IMAGE031
；
方程(45)中推力向量的非线性约束函数是非凸的；在运用凸优化方法前必须将其化为凸函数；故将式(45)中等号松弛为不等号，然后该约束就转化为凸约束，即：
Figure DEST_PATH_IMAGE083
(50)
从而，线性化后的动力学(48)、始末约束(43)(44)、凸约束(50)、性能指标(47)一起构成了凸子问题P2。
步骤五：通过线性化动力学和松弛非线性等式约束，将考虑发动机关机和参数强非线性 的连续推力轨迹优化问题凸化。
采用梯形法对式(48)中数值积分进行转化，从而将问题P2转化为凸优化问题形式；将转移任务时间区间[t 0 ，t f ]划分为N+1个节点，第i个节点处探测器的状态量和控制量记为X i 、U i ，则变量集为[X 0, ..., X i, ..., X N ]，待求解控制集为[C 0, ..., C i , , ..., C N ]，则动力学积分(48)转化为
Figure DEST_PATH_IMAGE084
(51)
其中i = 0, 1, …, N-1，
Figure DEST_PATH_IMAGE036AA
Figure DEST_PATH_IMAGE038A
Figure DEST_PATH_IMAGE085
；
从而，凸子问题P2转化为凸问题。
对凸问题进行连续逼近以便迭代求解得到最优的考虑关机约束的变参数连续推力转移轨迹，直至其解收敛于P1的解。
令k=0，给出初始状态向量的猜测值X0，由从初值X(t 0 )到终值X(t f )的直线获取；L(t)的猜测值选取从L(t 0 )到L(t f )＋2∏·Γ的一条直线，其中Γ为连续推力转移轨迹运行圈数的估计值；m(t)的猜测值选取m(t 0 )到m(t f )的一条直线。
对于第i+1次迭代，选取第k次迭代的解X k 作为状态向量的猜测初值，解对为
Figure 908966DEST_PATH_IMAGE062
。
检查是否满足收敛条件：
Figure DEST_PATH_IMAGE086
(52)
其中，其中，ρ为收敛精度要求；如果不满足式(52)，需要继续迭代求解，如果满足式(52)即得到问题P1的解X*＝X k+1、C*＝C k+1，即得到考虑发动机定期关机约束和强非线性变化参数的连续推力轨迹优化问题的解，X*为最优转移轨迹，C*为对应的最优控制方向。
经过优化，在13次迭代后满足收敛条件(52)，得到探测器在转移过程中燃料消耗为 292.6kg，发动机开机时长约20208h，考虑定轨约束发动机定期关机时能够准确到达末端约束处的状态。最优连续推力转移轨迹如图5所示，对应的最优推力幅值变化曲线如图6所示，其中阶梯状变化幅值与发动机推力特性相符，每条起于X轴的纵向竖线对应一次发动机开/关机，最优径向、切向、法向推力分量变化曲线如图7所示。作为对比，若发动机推力恒为 90mN，比冲恒为3100s，则转移过程燃料消耗为293.5kg，发动机开机时长约38156h，燃耗相差不大，但是由于考虑发动机参数变化和关机时，转移过程中整体推力大小大于恒定推力场景，所以发动机开机时长更短。考虑发动机参数变化和关机的场景更加符合实际工程要求，定期的关机定轨使得探测器更精确到达目标位置。
步骤六：根据步骤一至步骤五优化得到的最优转移轨迹及对应的最优控制方向，对实际工程中考虑定轨约束需定期关机以及发动机比冲和推力大小强非线性变化的连续推力场景适用性强，末端状态精度高。
以上所述的具体描述，对发明的目的、技术方案和优点进行了进一步详细说明。需要理解的是，以上所述仅为本发明的具体实施示例，用于解释本发明，并不用于限定本发明的保护范围。凡在本发明的精神和原则之内所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (2)

1.考虑发动机关机和参数强非线性的连续推力轨迹优化方法，其特征在于：包括如下步骤，
步骤一：根据关机约束和发动机类型给出强非线性推力和比冲变化模型；
步骤一实现方法为，
强非线性参数发动机的推力和比冲随探测器帆板的输入功率变化，而帆板输入功率与日心距成反比：
Figure 391899DEST_PATH_IMAGE001
(1)
其中r AU 是日心距离，单位为 AU，P 0为日心距为1AU时帆板的输入功率，P in 是日心距为r AU 时帆板输入功率；发动机最大推力T max和比冲随帆板输入功率的变化有两种模式，一种为多项式变化，一种是阶梯变化；将最大推力T max和比冲I sp 表示为输入功率的函数G和H，即
Figure 578161DEST_PATH_IMAGE002
(2)
探测器进行转移时，需要定期对探测器进行定轨，此时发动机需定期关闭；考虑发动机开关机约束用标志符ut表示为
Figure 238949DEST_PATH_IMAGE003
（3）
其中，t为飞行时间，t on 和t off 分别表示发动机的最大连续开机时长和执行定轨所需的单次最短关机时长；进而考虑关机约束的发动机最大推力T max为
Figure 709114DEST_PATH_IMAGE004
(4)；
步骤二：建立基于改进春分点根数考虑发动机关机和参数强非线性的转移轨迹动力学模型；
步骤二实现方法为，
采用改进春分点根数
Figure 797155DEST_PATH_IMAGE005
描述探测器的状态，其中
Figure 329768DEST_PATH_IMAGE006
,
Figure 935193DEST_PATH_IMAGE007
,
Figure 869651DEST_PATH_IMAGE008
,
Figure 394173DEST_PATH_IMAGE009
，
Figure 538715DEST_PATH_IMAGE010
，
Figure 806885DEST_PATH_IMAGE011
，a、e、i、ω、Ω和υ分别为探测器轨迹的半长轴、偏心率、轨道倾角、近地点幅角、升交点赤经和真近点角；则连续推力探测器的轨道动力学模型为
Figure 595850DEST_PATH_IMAGE012
(5)
其中
Figure 432219DEST_PATH_IMAGE013
(6)
且控制向量的系数矩阵为：
Figure 939424DEST_PATH_IMAGE014
(7)
Figure 480126DEST_PATH_IMAGE015
,
Figure 779390DEST_PATH_IMAGE016
,
Figure 380135DEST_PATH_IMAGE017
其中，I sp (r)和T(r)是与日心距相关的发动机强非线性比冲和推力大小；
Figure 515581DEST_PATH_IMAGE018
=1 + f cosL + g sin L, s 2 = 1 + h 2 + K 2，μ是中心天体的引力常数，m是探测器质量，C为控制向量，C=[fr,ft,fn,T]T，fr,ft,fn分别发动机径向、切向、法向推力大小；
步骤三：根据实际任务需求，给出考虑发动机关机和参数强非线性的连续推力轨迹优化问题的约束和性能指标；
步骤三实现方法为，
始末端探测器状态均由任务时间对应的星历给出；相应的改进春分点轨道根数为：
Figure 859975DEST_PATH_IMAGE019
(8)
Figure 623532DEST_PATH_IMAGE020
(9)
其中t 0 与t f 分别为始末时间；同时，初始质量m(t 0 )固定，末端质量m(t f )无约束；探测器的推力分量满足
Figure 50971DEST_PATH_IMAGE021
(10)
且
Figure 267189DEST_PATH_IMAGE022
(11)
期望得到燃料最优连续推力轨迹，故将优化问题的性能指标设定为：
Figure 680852DEST_PATH_IMAGE023
(12)
目的为最大化探测器末端质量，即最小化燃料消耗；
综上所述，燃料最优连续推力转移道轨迹优化问题P1总结为：
Figure 174282DEST_PATH_IMAGE024
约束方程：式(8) ~ (11)；
步骤四：通过线性化动力学和松弛非线性等式约束，将考虑发动机关机和参数强非线性的连续推力轨迹优化问题凸化；
步骤四实现方法为，
为了将探测器动力学凸化，将非线性动力学方程(5)基于小扰动的连续线性化方法近似；连续近似过程中第k次迭代存在一个解X k ；然后，第k+1次迭代过程中，在X k 存在的前提下将动力学方程线性化，主项H(X)在Xk附近线性化；故线性化后的动力学方程为：
Figure 647988DEST_PATH_IMAGE025
(13)
其中
Figure 351502DEST_PATH_IMAGE026
；状态向量系数矩阵为：
Figure 162332DEST_PATH_IMAGE027
(14)
其中，
Figure 900481DEST_PATH_IMAGE028
,
Figure 545089DEST_PATH_IMAGE029
，
Figure 876844DEST_PATH_IMAGE030
,
Figure 366732DEST_PATH_IMAGE031
；
方程(10)中推力向量的非线性约束函数是非凸的；在运用凸优化方法前必须将其化为凸函数；故将式(10)中等号松弛为不等号，然后该约束就转化为凸约束，即：
Figure 84021DEST_PATH_IMAGE032
(15)
从而，线性化后的动力学(13)、始末约束(8)(9)、凸约束(15)、性能指标(12)一起构成了凸子问题P2；
步骤五：通过数值积分和逐次逼近快速求解考虑关机约束的变参数连续推力轨迹优化问题，得到最优转移轨迹及对应的最优控制方向，即实现考虑发动机关机和参数强非线性的连续推力轨迹优化；
步骤五实现方法为，
采用梯形法对式(13)中的数值积分进行转化，从而将问题P2转化为凸优化问题形式；将转移任务时间区间[t 0 ，t f ]划分为N+1个节点，第i个节点处探测器的状态量和控制量记为X i 、U i ，则变量集为[X 0, ..., X i, ..., X N ]，待求解控制集为[C 0, ..., C i , , ..., C N ]，则动力学积分(13)转化为
Figure 633951DEST_PATH_IMAGE033
(16)
其中i = 0, 1, …, N-1，
Figure 312057DEST_PATH_IMAGE034
，
Figure 746580DEST_PATH_IMAGE035
，
Figure 193742DEST_PATH_IMAGE036
，
Figure 648994DEST_PATH_IMAGE037
；
从而，凸子问题P2转化为凸问题；
对凸问题进行连续逼近以便迭代求解得到最优的考虑关机约束的变参数连续推力转移轨迹，直至其解收敛于P1的解；
令k=0，给出初始状态向量的猜测值X0，由从初值X(t 0 )到终值X(t f )的直线获取；L(t)的猜测值选取从L(t 0 )到L(t f )＋2π·Γ的一条直线，其中Γ为连续推力转移轨迹运行圈数的估计值；m(t)的猜测值选取m(t 0 )到m(t f )的一条直线；
对于第i+1次迭代，选取第k次迭代的解X k 作为状态向量的猜测初值，解对为
Figure 204609DEST_PATH_IMAGE038
；
检查是否满足收敛条件：
Figure 301878DEST_PATH_IMAGE039
(17)
其中，ρ为收敛精度要求；如果不满足式(17)，需要继续迭代求解，如果满足式(17)即得到问题P1的解X*＝X k+1、C*＝C k+1；
至此，得到考虑发动机定期关机约束和强非线性变化参数的连续推力轨迹优化问题的解，X*为最优转移轨迹，C*为对应的最优控制方向。
2.如权利要求1所述的考虑发动机关机和参数强非线性的连续推力轨迹优化方法，其特征在于：还包括步骤六，根据步骤一至步骤五优化得到的最优转移轨迹及对应的最优控制方向，解决推力控制领域相关问题。
CN202210123114.XA 2022-02-10 2022-02-10 考虑发动机关机和参数强非线性的连续推力轨迹优化方法 Active CN114154253B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210123114.XA CN114154253B (zh) 2022-02-10 2022-02-10 考虑发动机关机和参数强非线性的连续推力轨迹优化方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210123114.XA CN114154253B (zh) 2022-02-10 2022-02-10 考虑发动机关机和参数强非线性的连续推力轨迹优化方法

Publications (2)

Publication Number Publication Date
CN114154253A CN114154253A (zh) 2022-03-08
CN114154253B true CN114154253B (zh) 2022-05-10

Family

ID=80450135

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210123114.XA Active CN114154253B (zh) 2022-02-10 2022-02-10 考虑发动机关机和参数强非线性的连续推力轨迹优化方法

Country Status (1)

Country Link
CN (1) CN114154253B (zh)

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110736470A (zh) * 2019-11-06 2020-01-31 北京理工大学 一种不规则形状小天体附近连续推力轨道混合搜索方法
CN110806212A (zh) * 2019-11-12 2020-02-18 北京理工大学 基于逐次凸规划的小行星探测小推力转移轨迹优化方法
CN113111434A (zh) * 2021-03-30 2021-07-13 北京航空航天大学 一种基于凸混合整数规划的组合动力飞行器轨迹优化方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110736470A (zh) * 2019-11-06 2020-01-31 北京理工大学 一种不规则形状小天体附近连续推力轨道混合搜索方法
CN110806212A (zh) * 2019-11-12 2020-02-18 北京理工大学 基于逐次凸规划的小行星探测小推力转移轨迹优化方法
CN113111434A (zh) * 2021-03-30 2021-07-13 北京航空航天大学 一种基于凸混合整数规划的组合动力飞行器轨迹优化方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
Rapid planning for aerocapture trajectory via convex optimization;Hongwei Han 等;《Aerospace Science and Technology》;20190131;第84卷;第763-775页 *
基于同伦方法的地月系L2点小推力转移轨道优化;潘迅 等;《载人航天》;20190228;第25卷(第1期);第25-30页 *
小天体多目标多模式探测任务设计;黄江川 等;《中国科学: 物理学 力学 天文学》;20190617;第49卷(第8期);第1-11页 *

Also Published As

Publication number Publication date
CN114154253A (zh) 2022-03-08

Similar Documents

Publication Publication Date Title
Hu et al. Robust saturated finite time output feedback attitude stabilization for rigid spacecraft
Tang et al. Fuel-optimal low-thrust trajectory optimization using indirect method and successive convex programming
Huo et al. Electric sail thrust model from a geometrical perspective
Shamma et al. Gain-scheduled missile autopilot design using linear parameter varying transformations
Wu et al. Nonfragile output tracking control of hypersonic air-breathing vehicles with an LPV model
Lu et al. Flight control design for small-scale helicopter using disturbance-observer-based backstepping
CN110806212B (zh) 基于逐次凸规划的小行星探测小推力转移轨迹优化方法
Ren et al. Singular perturbation-based fault-tolerant control of the air-breathing hypersonic vehicle
Wang et al. Nonlinear hierarchy-structured predictive control design for a generic hypersonic vehicle
Kuipers et al. Adaptive control of an aeroelastic airbreathing hypersonic cruise vehicle
Wang et al. Sliding mode decoupling control of a generic hypersonic vehicle based on parametric commands
Huo et al. Altitude and velocity tracking control for an airbreathing hypersonic cruise vehicle
CN114154253B (zh) 考虑发动机关机和参数强非线性的连续推力轨迹优化方法
Gui et al. Robustness analysis and performance tuning for the quaternion proportional–derivative attitude controller
Zhang et al. On-line ascent phase trajectory optimal guidance algorithm based on pseudo-spectral method and sensitivity updates
CN109625332B (zh) 一种平动点轨道交会无需初始误差符号的预设性能控制方法
Zhou et al. Ascent trajectory optimization for air‐breathing vehicles in consideration of launch window
Martin et al. Velocity pointing error reduction for spinning, thrusting spacecraft via heuristic thrust profiles
Pontani et al. Neighboring optimal guidance and attitude control of low-thrust earth orbit transfers
Kluever et al. Trajectory-tracking guidance law for low-thrust earth-orbit transfers
Liu et al. Ascent trajectory optimization for air-breathing hypersonic vehicles based on IGS-MPSP
Zhang et al. Velocity-to-be-gained deorbit guidance law using state space perturbation method
Wu et al. Control of wing rock based on high order sliding mode disturbance observer
Wu et al. Symplectic transformation based analytical and numerical methods for linear quadratic control with hard terminal constraints
Li et al. A convex approach for trajectory optimization of super-synchronous-transfer-orbit large launch systems with time-position constraints

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