Transcription

1 博士研究生学位论文 题目 : 有限推力航天器最优轨道机 动控制 姓 名 : 岳新成 学 号 : 院 系 : 工学院 专 业 : 力学系统与控制 研究方向 : 航天器动力学与控制 导师姓名 : 黄琳院士杨莹副教授 二 一 年四月

2

3 版权声明 任何收存和保管本论文各种版本的单位和个人, 未经本论文作者同意, 不 得将本论文转借他人, 亦不得随意复制 抄录 拍照或以任何方式传播 否 则, 引起有碍作者著作权之问题, 将可能承担法律责任

4

5 中文摘要 摘 要 航天器轨道机动控制问题一直是航空航天领域的研究热点 研究轨道机动的方法主要有两类 : 一是脉冲法, 二是制导法 脉冲法由于假设卫星在瞬时就能获得所需速度增量而不引起卫星向径的突变而大大简化了问题, 因而得到广泛研究, 理论比较成熟 但当发动机推力比达到千秒级时, 轨道机动需要花费更多的时间, 就必须把推力作为连续函数来考虑 连续推力下的最优轨道转移问题的数值方法主要有直接法 间接法和混合法 本论文主要应用间接法研究连续推力航天器最优轨道机动问题 考虑航天器为常值推力发动机点质量模型, 在二体轨道基础上, 采用位置 速度和改进轨道根数作为状态变量描述了轨道动力学方程 轨道机动模式包括轨道转移 交会和拦截, 应用极大值原理设计了时间或能量最优控制律, 并进行了奇异性分析 根据极大值原理一阶必要条件, 将轨迹优化问题转化为两点边值问题, 然后在数值上应用打靶法来求解相应的边值问题 论文还研究了最优轨道规避问题, 在给定中间过渡轨道的条件下, 应用非线性规划方法研究了时间最优和燃料最优轨道规避问题 具体的研究工作和成果如下 : 1. 应用开普勒型动力学方程研究了时间最优轨道拦截问题 考虑控制输入分别在惯性坐标系和轨道坐标系下分解, 应用极大值原理设计最优控制律并进行了奇异性分析, 根据极大值原理一阶必要条件, 将轨道拦截问题转化为两点边值问题来求解 数值仿真结果表明最优状态解具振荡性, 且同样存在最大推力幅值与最小飞行时间乘积近似成常数关系的现象 2. 应用改进轨道根数描述的高斯型动力学方程, 在统一的框架下研究了时间最优轨道转移 交会和拦截问题 应用间接法设计了最优控制律, 极大值原理一阶条件将轨迹优化问题转化为两点边值问题 仿真结果表明对于相同的最大推力幅值, 轨道转移和交会问题具有几乎相同的最优轨迹 ; 调整最大推力幅值, 重复数值实验, 结果表明最大推力幅值与最小飞行时间乘积近似成常数关系的猜想同样适用于轨道拦截和交会问题 3. 研究了累积赤经最优轨道转移问题 将累积赤经引入系统作为新的自变量, 应用间接法研究考虑变质量和推力约束的最优累积赤经轨道转移问题 ; 借鉴平均化方法, 设计了近似最优控制律, 同样将最优轨道转移问题转化为两点边值问题 数值仿真结果表明, 最大推力幅值与最小累积赤经乘积近似成常数 I

6 中文摘要 关系 4. 研究了燃料最优轨道转移问题 基于改进轨道根数描述轨道动力学方程, 应用间接法研究燃料最优控制 由于燃料最优控制律存在奇异性, 不易求解, 转而应用 ϵ 算法将燃料最优和能量最优结合起来考虑 根据极小值原理一阶必要条件, 将轨迹优化问题转化为两点边值问题, 并在数值上对这三种形式的最优轨道转移做出比较, 结果表明能量最优 燃料最优及能量 - 燃料混合最优三种不同性能指标对应的最优轨迹几乎是相同的 5. 研究常值连续推力时间最优和燃料最优规避问题 在中间过渡轨道已知的条件下, 将规避问题划分为转移和返回两个阶段, 然后应用非线性规划方法统一来考虑时间或燃料最优规避问题 数值仿真结果表明, 对于同一最大发动机推力幅值, 时间最优情形与燃料最优情形相比较, 航天器完成规避任务历经的飞行时间短, 飞行轨迹的椭圆轨道半长轴变化较小, 但是绕飞地球的圈数多 数值仿真同时也表明, 在考虑时间最优或燃料最优轨道规避时, 航天器均一直以最大推力工作 关键词 : 有限推力航天器 ; 轨道机动 ; 轨道规避 ; 间接方法 ; 极大值原理 ; 两 点边值问题 ; 打靶法 II

7 On Finite-thrust Spacecraft Orbital Maneuver Optimal Control 英文摘要 Xincheng Yue (Mechanical System and Control) Directed by Professor Lin Huang and Associate Professor Ying Yang An important problem in astronautics is to transfer a satellite between elliptic orbits, which has been widely studied by many researchers in both the impulsive case and the continuous low-thrust case. Recently, much attention has been focused on the computation of the optimal control law for the case of the finite-thrust orbital transfer, which can be performed by minimizing the cost of the final time or the fuel-consumption under some additional constraints. The finite-thrust trajectory optimization methods fall mainly into three categories: indirect methods, direct methods and hybrid methods. In this dissertation, the finite-thrust optimal orbital maneuver problems are studied by using the indirect method, and the trajectory optimization problems are deduced into the two-point boundary-value problem by the Maximum Principle. The main contributions of this dissertation can be summarized as follows: 1.The finite-thrust time-optimal orbital intercept problem is studied firstly. The dynamics of the spacecraft is described by the so-called Keplerian equation using the position-speed variables in the Cartesian coordinate and in the orbital coordinate attached to the spacecraft. The optimal control laws are designed by the the Pontryagin Maximum Principle, and the singularity is also analyzed, then the trajectory optimization problem is converted into a two-point boundary-value problem. Two special examples are studied in numerical simulations, and the results show that the product of the minimum flight time and the maximal thrust is approximately a constant in these three cases of orbital maneuver. 2.The problems of indirect optimization of continuous finite-thrust minimum-time orbital maneuvers including the transfer, interception and rendezvous are studied. In virtue of the Pontryagin Maximum Principle, the trajectory optimization problem is converted into a two-point boundary-value problem, and the terminal constraints for these three orbital maneuver missions are also considered respectively. The single shooting method is used to solve the corresponding boundary-value problem, and the simulating results show that the product of the minimum flight time and the maximal III

8 英文摘要 thrust is approximately a constant in these three cases of orbital maneuver. 3.The continuous finite-thrust near minimum cumulative longitude orbital transfer problem using the indirect method is addressed. The movement of the satellite is described by the Gauss equation using the modified equinoctial elements and replacing the time as a system independent variable by the cumulative longitude. The Maximum Principle is adapted to design the optimal control in order to minimize the final cumulative longitude, and a two-point-boundary-value problem is derived from the orbital transfer problem. The single shooting method is applied in the numerical experiment, and the results of simulations demonstrate that the orbital transfer mission is achieved and the product of the maximal thrust and the minimum cumulative longitude is nearly a constant. 4.The problems of optimal continuous-thrust orbital transfer for minimizing energy, fuel and energy-fuel combination are studied by using indirect method. The fuel optimal orbital transfer problem is hard to solve because of the possible existence of the singularities in solution, but the minimum energy problem is relatively easy to deal with, so the ϵ-algorithm is applied to the fuel optimal problem by adding the energy term in the preference index. The trajectory optimization problem is converted into a two-point-boundary-value problem, and the shooting method is used to dispose to the problem in numerical computation. The results of the simulations demonstrate that the optimal trajectories are almost alike for the three optimal cases under the same thrust and transfer time. 5.The nonlinear programming method is used to study the continuous-thrust minimum-time and minimum-fuel-consumption orbital avoidance. The dynamics of the spacecraft is described by the modified elements. For the given intermediate parking orbit, the orbital avoidance mission is divided into transfer and round-trip phases. Considering the minimum-time or minimum-fuel-consumption performance, the optimal orbital avoidance problem is established. Furthermore, the nonlinear programming method is applied to the two phase orbital transfer problems. The results of simulation demonstrate that the minimum time/fuel-consumption orbital avoidance mission is well accomplished. The simulating results also show that the spacecraft flies around the earth with the maximum magnitude of the thrust at whole time history both for the minimum-time and minimum-fuel-consumption cases. Key Words: Finite-thrust Spacecraft; Orbital Maneuver; Orbital Avoidance;Indirect IV

9 英文摘要 Method; the Maximum Principle; Two-point boundary-value Problem; Shooting Method V

10

11 目 录 目 录 摘要 I ABSTRACT( 英文摘要 ) III 第一章 绪 论 研究背景 研究现状 轨道动力学描述 轨迹优化方法 轨道机动使命和要求 研究内容与结构安排 第二章 航天器轨道动力学基础 二体问题及轨道根数 开普勒型方程 控制在地心坐标系下的动力学方程 控制在轨道坐标系下的动力学方程 轨道根数描述高斯方程 极坐标方程 卫星质量方程 小结 第三章 时间最优轨道机动 Ⅰ: 开普勒型方程 时间最优轨道机动问题描述 最优轨迹设计 地心坐标系描述情形 轨道坐标系描述情形 数值仿真 算例 算例 VII

12 目 录 3.4 小结 第四章 时间最优轨道机动 Ⅱ: 轨道根数描述方程 问题描述 正则最优控制律设计 边值问题和横截条件 轨道转移 轨道拦截 轨道交会 数值仿真 小结 第五章 近似最优赤经轨道转移 问题描述 近似最优控制律设计 数值实验 小结 第六章 能量 / 燃料最优轨道转移 问题描述 轨道转移最优控制律设计 燃料最优 能量最优 ϵ 算子法求解燃料最优 数值计算 小结 第七章 连续推力最优轨道规避 问题描述 非线性规划求解最优规避问题 数值仿真 小论与展望 VIII

13 目 录 第八章结论与展望 参考文献 附录 A 符号与常数示例 附录 B 极大值原理 附录 C 两点边值问题与打靶法 个人简历 在学期间的研究成果 致谢 IX

14

15 表格 表格 3.1 最小飞行时间 时间最优轨道机动问题的数值条件 最小飞行时间 累积赤经最优转移数值边值 燃料 / 能量最优轨道转移数值边界条件 最优规避问题数值仿真边值条件 XI

16 插图 插图 2.1 轨道根数与坐标系简单示意图 球坐标系统示意图 时间最优轨道拦截 3D 轨迹,1N 时间最优轨道拦截状态变化轨迹,1N 时间最优轨道拦截最优控制轨迹,1N 最小飞行时间和最大推力幅值关系 时间最优轨道拦截 3D 轨迹,1N 时间最优轨道拦截状态变化轨迹,1N 时间最优轨道拦截最优控制轨迹,1N 时间最优轨道拦截 3D 轨迹,2N 时间最优轨道拦截状态变化轨迹,2N 时间最优轨道拦截最优控制轨迹,2N 状态最优轨迹,4N 轨道拦截 3D 最优轨迹,4N 轨道交会的 3D 最优轨迹,4N t fmin T max 近似常数关系 初始时刻 L 对轨道拦截时间的影响 累积赤经最优轨道转移 3D 最优轨迹,.8 N 累积赤经最优轨道转移状态的最优解,.8N 累积赤经最优轨道转移控制输入,.8 N L fmin T max 近似常数关系 不同推力幅值下 σ 变化轨迹 状态随时间变化轨迹 最优控制随时间变化轨迹 能量最优 3D 轨迹 燃料最优 3D 轨迹 XII

17 插图 6.5 终端时刻卫星质量与最大推力幅值的关系 时间最优规避 3D 最优轨迹,3N 最小燃料消耗规避 3D 最优轨迹,3N 时间最优规避状态最优变化轨迹,3N 最小燃料消耗规避状态最优变化轨迹,3N C.1 打靶法简单流程示意 XIII

18

19 第一章绪 论 第一章绪论 1.1 研究背景早在半个世纪前, 人类开始了空间技术的研究, 虽然这是一个相对较新的研究领域, 但是已经取得了巨大的成功 自从第一个人造卫星于 1957 年发射升空之后, 人类和平利用开发太空的活动, 诸如深空探测 月球探测 空间站等等, 都取得令人瞩目的成就, 也激励着众多科学家将兴趣聚焦在空间科学技术研究上 上世纪六十年代, 包括地球轨道间转移 地月转移 星际航行等在内的轨道机动任务的轨迹优化研究引起了众多科学家的研究兴趣 对于大多的空间使命, 诸如有效载荷和飞行时间等指标都是事先严格设计好的 携带有效载荷过多或者飞行时间过长都会降低效率, 或者完成一次任务需要实施更多次的空间机动飞行 因此, 空间任务的优化设计是不可避免的 当前, 航天器轨道控制技术主要基于航天器的推进系统, 因此根据推进方式的不同, 研究轨道机动的方法主要有两类 : 一是脉冲法, 二是制导法 [1, 2] 早期的航天器主要采用化学推进系统, 产生的推力较大, 可以近似看作是脉冲推力, 可以通过解析的办法来计算常值推力弧 脉冲法在理论上假设发动机推力充分大, 发动机的推力工作时间与轨道运行周期相比是相当短暂的, 航天器在瞬时就能获得所需速度增量而不引起卫星向径的突变 [1, 2] 此时发动机的推力比冲 I sp 约为 25~45s, 这一假设大大简化了问题的研究 早在上世纪五六十年代, 脉冲控制轨道机动包括轨道转移 拦截和交会问题已得到广泛深入的研究 [3] 早期的代表性工作是 Lawden 开展的, 他把燃料最优轨迹优化问题看做一个 Mayer 型问题, 应用变分法和拉格朗日乘子法寻求最优轨迹, 发展了一张 主矢量 理论, 脉冲机动的最优方向沿主矢量方向 [4] Lawden 的工作也证实了 Hohmann 关于最小燃料消耗脉冲轨道转移的结果 [1, 2] 对于轨道机动, 脉冲式大推力虽可使航天器在短时间内获得较大的速度增量, 但所需的推进剂多, 从而减少了有效载荷 随着空间探测任务的深入, 化学火箭发动机比冲小 效率低的缺点日益暴露, 不能满足未来任务的需要 与此同时, 以电推进 离子推进 太阳帆推进为代表的小推力推进技术受到了广泛的关注, 并迅速发展, 化学 - 电子混合推进方式同时也得到应用 [5], 1

20 1.2. 研究现状 越来越多的空间任务选择这种推力持续变轨的机动模式 1998 年发射的 Deep Space 1 第一次成功地使用太阳能电力推进系统作为星际航行的主要推进方式 [6];23 年 9 月, 欧洲第一个月球航天器 SMART-1 应用电推进系统使得航天器逃离地球引力并进入月球轨道 另一个小推力推进系统的重要例子是太阳帆, 航天器没有发动机, 直接靠太阳能产生动力来驱动 [7] 这些小推力推进系统的特点是, 推力小, 不能在短的轨道周期内完成轨道机动任务, 通常推力弧是连续的, 因而轨迹计算需要多次综合才能完成 小推力推进也常称作有限推力或连续推力, 此时, 发动机比冲 I sp 达到千秒级, 轨道机动需要花费更多的时间, 必须把推力作为连续函数来考虑 尽管表面上看, 小推力推进系统在技术上比较复杂, 需要的机动任务时间也较长, 但是小推力推进系统具有明显的有效载荷优势, 而且与传统的化学推进系统相比, 小推力推进系统能够提供更高级的使命设计任务和方法, 因而激起了众多学者的研究兴趣 本论文主要研究小推力时间最优或者能量最优轨道机动问题, 所指的小推力推进系统并不特指是电推进或者太阳帆推进, 涉及的小推力轨迹优化问题广泛的指非脉冲式使命设计并具有精确的微分运动方程描述 本论文的小推力轨道机动模式主要包括地球卫星之间的轨道转移 交会 拦截和规避问题 小推力最优轨道机动即轨迹优化和使命设计的研究极具挑战性, 所研究的问题括轨道动力学选择 轨迹优化方法 连续推力弧的计算 连续控制函数的选取以及不同的使命要求等等, 以下就这几方面当前的研究现状进行综述 1.2 研究现状关于有限推力最优轨道机动问题, 国内外已经开展了广泛深入的研究 Yang Gao 的博士论文应用了两个直接优化方法 ( 混合法和直接打靶法 ) 研究了小推力轨道转移问题 [8];James D.T. 研究了连续推力时间最优轨道转移问题, 应用欧拉 - 拉格朗日理论将轨迹优化问题转化为两点边值问题来求解, 并应用 Kustaanheimo-Stiefel(KS) 变换研究了最优初始协状态问题 [9];Mischa Kim 则基于数值和解析的方法研究了连续推力时间最优轨道转移问题, 应用解析和半解析技术提高了最优解寻优过程 [1]; 而 Troy Goodson 则主要应用间接法研究了有限推力燃料最优轨道转移问题, 并修改了 Patched Method 将其应用于最优轨道转移问题 [11] 国内方面, 哈尔滨工业大学任远的博士论文则是应用遗传算法研究了深空探测中的小推力轨道优化问题 [12] 2

21 第一章绪 论 轨道动力学描述 研究有限推力轨道机动问题, 一个重要内容是选择合适的轨道动力学描述 航天器轨道动力学的基础是二体轨道问题 [1, 2, 13, 14], 在地球平方反比力场中将航天器看做点质量体, 无控轨道运动满足二体开普勒方程 考虑外加控制输入或摄动的开普勒型动力学方程, 状态变量为位置和速度, 许多文献即采用该方程研究轨道优化问题 [15 21] 由于位置 速度分量描述的轨道动力学方程的解具有振荡性, 而轨道根数变化缓慢平滑, 数值上可以采用较大的积分步长, 而且用轨道根数描述的动力学方程更易于分析控制输入或者摄动对轨道的影响, 因而得到更广泛的应用, 如文献 [22 24] 采用经典的轨道根数描述方程研究了小推力最优轨道优化问题 但是, 由于经典的轨道根数描述的动力学方程在零偏心率和零倾角时存在奇异性, 一种称之为赤道轨道根数 (Equinoctial Elements) 的非奇异轨道根数被引入轨道机动研究 [14, 25, 26],Walker 等人提出的另一种改进的赤道轨道根数 (Modified Equinoctial Elements, MEE) [27, 28] 由于适用于所有二次曲线轨道, 因而用其描述的高斯型动力学方程 [29] 被更多学者采用 [3 33] Kechichian 也引入了一系列的 Equinoctial Elements 来研究小推力地球卫星间轨道转移问题 [34 37] 对于实际的轨道机动任务, 都是在特定的坐标系系统下实施, 笛卡尔惯性坐标系被自然地应用到位置和速度作为变量描述的开普勒型方程中 [15 21]; 但是为便于优化设计, 轨道根数描述的动力学方程多采用轨道固连坐标系, 在这种坐标系中, 通常控制推力分量有两种分解方式 : 切向 - 法向分解和径向 - 法向分解 [29, 3]; 球坐标或者二维极坐标系统也广泛应用于轨道转移研究中 [38 4] 轨迹优化方法 有限推力最优轨道机动问题, 在数值研究方法上归结为轨迹优化研究 轨迹优化数值方法按照传统的分类主要有间接法 直接法和混合法 [8 11,41,42], 简要分类如图 1.1 所示 [43] 关于轨迹优化的综述文献,Betts 在 1998 年发表的一篇总结早期的轨迹优化数值方法的文章被广泛引用 [41] Hull 从参数化变量类型 数值积分方法和积分阶次角度对参数化方法进行了分类和概述 [44] Ross 等人则从最优控制问题的转换方法 变量离散和算法三方面对轨迹优化方法进行了分类 [45] 文献 [42] 则从应用角度出发, 结合近年来轨迹优 3

22 1.2. 研究现状 化研究领域的进展, 阐述了不同分类标准下的常用轨迹优化方法以及高斯伪谱等新方法的基本原理 优缺点以及应用背景 国内学者如崔平远 [46 48] 唐国金团队 [49, 5] 杨涤团队 [51, 52] 等多采用直接法 混合法或直接搜索法对连续推力轨道机动或星际任务开展了广泛的研究 4

23 第一章绪 论 变分法 最速下降法 间接法 极大值原理 共轭梯度法 牛顿法 Patched 方法 边值打靶法 边界极值法 动态规划 最速下降法 共轭梯度法 轨迹优化方法 直接法 一般非线 性规划法 约束非线性规划 无约束非线性规划 变尺度法高斯 - 牛顿最小二乘法模式收索法单纯形法 SWIFT 法拉格朗日乘子法可行方向法投影梯度法 单纯形法 简单规划法 遗传算法 序贯二次规划法 随机收索法 模拟退火 微粒群法 神经网络 图 1.1 轨迹优化数值方法 5

24 1.2. 研究现状 间接法应用极大值原理 (Pontryagin Maximum Principle,PMP) 或变分法 [53 55] 推导最优控制的一阶必要条件, 将连续空间的最优控制问题转化为一个两点边值问题 (two-point boundary-value problem, TPBVP), 通过单值打靶法 多点打靶法等求解相应的边值问题来获得最优轨迹 间接法的优点, 一是解的精度较高, 收敛速度快, 最优解满足一阶最优必要条件, 便于推导 ; 其不足有三点 [41, 55]:1) 基于极大值原理推导一阶必要条件的过程比较复杂和烦琐 ;2) 求解两点边值问题的收敛域较小, 对初始猜测值比较敏感, 对抽象的 ( 无物理意义的 ) 协变量初始值猜测比较困难 ;3) 对于有路径约束的最优控制问题, 间接法处理比较困难 直接法采用参数化方法将连续空间的最优控制问题转化为一个非线性规划 (Nonlinear Programming,NLP) 问题, 通过数值求解非线性规划问题来获得最优轨迹 [56] 根据参数化方法不同, 直接法分为两种基本类型 [42]:1) 通过将控制变量离散化, 将连续空间的控制变量参数化, 而状态变量由数值积分获得 ;2) 同时将控制变量和状态变量离散化, 节点之间采用多项式表示状态变量随时间的变化关系, 连续动力学微分方程约束转化为代数约束 直接法的优点, 一是不需要推导最优条件 ; 二是收敛域相对间接法更宽广, 对初始估计精度要求不高, 不需要猜测协变量的初值, 不需要切换结构的先验信息 但是许多直接法不提供主矢量信息, 因此不能保证非线性规划解是原最优控制问题的解, 而将控制变量参数化的直接法容易收敛到局部最优解, 即依赖于初始猜测值 直接法还有一个显著缺点就是依赖节点的选取, 计算量大, 计算精度不高, 解不够光滑 混合法, 也引入协状态, 依据极大值原理得到正则方程组, 将问题转化为在终端约束条件下, 针对协状态初值及其他一些特定参数 ( 如发动机开关机时刻, 总飞行时间等 ) 的参数优化问题 间接法已成功应用于大量的小推力轨迹优化问题 Kechichian 采用了改进轨道根数描述的动力学方程研究了最优轨道转移问题, 应用极大值原理一阶必要条件给出协状态方程, 并推导出了协状态方程的具体表达式 [36, 37] C. H. Chuang 等人应用开普勒型方程, 并考虑了 J2 摄动和大气摄动, 应用间接法研究了燃料最优轨道转移问题, 并应用横截条件确定了终端边值约束 [19, 57, 58] Oberle, H.J. 应用极坐标方程, 研究了最小燃料消耗地球轨道 - 火星轨道转移问题, 应用极大值原理必要条件, 并证明了解的存在性 [59] Bonnard, B.,Caillau, JB, Gergaud, J. 等人所在的研究小组系统的应用间接法研究了时间最优 能量最优轨道转移问题, 用改进的轨道根数 MEE 描述卫 6

25 第一章绪 论 星动力学方程, 应用几何的观点分析了最优控制特性 [3 32, 6 63] Bonnard, B. 和 Chyba, M. 的专著 [64] 也详细的讨论了极大值原理的应用问题, 特别是奇异控制问题, 虽然该书并未强调是应用于航天器最优轨道机动问题, 但是书中的成果都是基于文献作者关于航天器最优轨道转移问题的研究成果, 完全适用于轨道机动问题的研究 间接法即是应用 Pontryagin 极值原理将轨迹优化问题转化为两点边值问题, 两点边值问题的求解有大量的文献可考, 应用在最优轨道机动上的数值方法主要有单值打靶法, 多点打靶法,Patched 方法和最小边界条件法 (Minimizing- Boundary-Condition method,mbcm) 单值打靶法首先确定打靶函数, 令其满足边值条件, 然后寻求合适的算法求解, 常用的有牛顿型算法 文献 [61] 介绍了单值打靶技术的应用, 文献 [3, 31, 65] 应用单值打靶法求解了时间最优轨道转移问题, 文献 [18, 32, 6] 研究了其在能量最优轨道转移问题上的应用 由于 Hamiltonian 系统的不稳定性, 某些情况单值打靶法是不充分的, 因而必须采用多点打靶法 应用多值打靶法必须预先知道最优控制的结构, 比较常用的算法是 H.J.Oberle 发展的 Boundsco 算法, 但这个算法需要猜测切换时间点 文献 [11, 57, 59, 66 68] 应用 Boundsco 算法研究了最优轨道机动问题 MBCM 实际上是单值打靶法的改进, 通过消除一个边值条件而增大了可行解集, 未知边界条件的个数并没改变 因为解集的扩大, 相应的边值问题易于求解, 一旦完成求解, 寻找与最终边值条件不符的解就可看作最小化问题, 计算梯度用于更新初始状态直至最后一个边界条件满足 在 MBCM 中, 最优控制的切换结构是明显的, 因为每一个积分步骤都将检验切换函数, 因而可以直接应用切换条件, 收敛域比单值打靶法大 文献 [11, 57] 应用 MBCM 研究了最优轨道转移问题 Goodson T. 也用 Patched 方法研究了能量最优轨道转移问题 应用打靶法求解两点边值问题由于需要猜测协状态初值而变得难于处理, 系统的收敛域对协状态的初始值非常敏感, 因此协状态的初值分析和设计也是最优控制领域当前的一个研究热点,David G. Hull 系统的总结了打靶法中初始协状态值的处理方法 [69] 常用的一种方法是规范变换法(Canonical transformations), 考虑不同坐标系或者变量之间的变换, 选择合适的坐标系或者变量来代替原协状态初始值 [7], 文献 [21, 71] 应用规范变换方法研究了最优轨道转移问题 文献 [72] 应用球坐标系统, 研究了最小燃料消耗最优轨道转移问题, 将协状态变量在初始时刻领域用一阶近似泰勒级数表示, 通过求解初始时刻领域内的方程来得到近似协变量初始猜测值 生成函数方法由于不需猜测协 7

26 1.2. 研究现状 状态的初始值而广泛应用于最优交会问题 [73 75] 生成函数有四种形式 : F 1 (x, x, t, t ), F 2 (x, p, t, t ), F 3 (p, x, t, t ), F 4 (p, p, t, t ) 最优交会初始状态和终端状态都是已知的, 这样就可以选取 F 1 为生成函数, 相 应有变换 : p = F 1(x, x, t, t ), p = F 1(x, x, t, t ) x x 其他生成函数可以通过 Legendre 变换得到 [74] 但是选取合适的生成函数也是 比较困难的事情, 而且由于存在多解而可能出现奇异 应用极大值原理设计最优控制律时, 另外一个问题是存在奇异控制 [76], 此时需要继续判定二阶条件, 分析切换函数的特性 文献 [62] 研究了时间最优 共面轨道转移情形, 证明了切换函数是连续可微的而且仅存在有限多个孤立零 点, 即不存在奇异弧而仅存在切换点, 而在任意两个切换点之间, 解是光滑连 续的 进而证明了在奇异点都是 π 奇异, 即在切换点控制输入方向旋转 18 度 文献 [3, 61] 应用幂零逼近的方法也同样证明了 π 奇异性 虽然文献 [62] 证明了 在近地点或远地点若存在奇异点则至多有三个连续切换点, 但是在通常情况下 是否存在奇异点, 以及若存在奇异点则切换时间的确定, 仍然是未解决的难 题, 在数值算法中则可以通过调整积分步长来处理奇异性 [77] 根据极大值原 理, 极值解满足一阶必要条件, 局部充分条件则需满足 Legendre-Clebsche 二阶 条件 2 H u 2 < Morse 定理建立起二阶条件与共轭点之间的联系 [64], 在第一个共轭点之 前的轨迹都是局部时间最优的 [64, 78, 79] Bonnard 等人提出了正则情况下 的 COTCOT 算法, 即检验 Jacobi 方程的解组成的 Jacobi 域的向量空间维数, 用于 检验共轭点, 并将其用于时间最优轨道转移 [78, 79], 但该算法需要预先知道轨 道转移时间, 而转移时间是我们需要优化的目标, 因而该算法在数值上并不能 准确检验二阶条件 而在奇异情况下, 单输入仿射系统同样可以检验共轭点, 但是多输入系统的二阶条件检验则是开问题 [64] 文献 [8] 则基于 Riccati 方法检 验终端时间自由最优控制问题的二阶充分条件, 并将其应用于时间最优地球 - 火 星轨道转移问题 而文献 [81, 82] 则考虑了单终端约束和多终端约束情况下的二 阶条件检测过程, 并将其应用于固定时间能量或燃料最优轨道转移问题 而对 8

27 第一章绪 论 于切换函数, 文献 [83] 假设应用同伦方法确定的切换函数不存在奇异弧即仅存在孤立的零点 ; 文献 [18] 则直接研究了燃料最优轨道转移, 假设切换函数仅存在孤立零点, 应用二分法检测切换点, 提高了打靶法的收敛性和运行效率 文献 [15] 则基于考虑 大气摄动的开普勒型动力学模型, 提出了一种基于振荡频率分析称之为采样法的切换算法, 能够提高运算速度和切换点的精确度 应用间接法研究轨迹优化问题还有一个难点就是路径约束或者状态约束, 虽然关于控制约束 状态等式或不等式约束问题, 理论上分析与设计已很完善 [55, 76], 但是在实际问题的数值处理中仍很困难, 如很多文献在研究最优轨道转移问题时都假设路径约束自然满足 [18,31,32,6,62] 文献 [21,66] 在研究返回式空间探测任务中, 应用开普勒型方程, 考虑能量最优指标, 并假定在目标轨道停留时间必须大于一定时间, 则将返回式问题分成转移段 停留段和返回转移段, 在目标轨道上的进入点和离开点约束则作为角点条件来处理, 应用极大值原理将问题转化成多点边值问题来求解 关于直接法, 一般是将轨迹优化问题转化为非线性规划问题 Betts 的专著系统的介绍了非线性规划方法的原理 方法和应用 [41], 他也应用非线性规划方法研究了最优轨道转移问题 [84 86] 文献 [33] 研究了时间轨道转移问题, 应用 direct-transcription 方法将连续最优控制问题转化为 NLP 问题 配点法 (collocation) 也应用于地球轨道转移任务中 [87, 88]; 文献 [89] 应用球坐标方程, 采用 sequential quadratic Programming(SQP) 和配点法研究了轨道转移问题 文献 [9] 则考虑了分段连续模型, 应用参数优化方法和二阶梯度方法研究了最小燃料消耗轨道转移问题 包括神经网络 模拟退火 进化算法 遗传算法等的随机直接收索方法也归入直接法范畴内 文献 [39, 91] 应用遗传算法 文献 [92] 应用进化算法研究了最优轨道转移问题 由于直接法和间接法各有优缺点, 将两个方法结合使用的混合方法也广泛应用于最优轨道转移中, 文献 [91] 应用开普勒型方程, 利用极大值原理将轨迹优化问题转化成两点边值问题, 然后应用遗传算法求解相应的边值问题, 而文献 [92] 则是用进化算法求解 文献 [93] 对球坐标系统描述的异面开普勒轨道转移问题, 应用极大值原理设计了最优控制律, 然后将时间区间分割, 节点作为内点约束, 再利用直接配点法来求解 轨道机动使命和要求 小推力 ( 或有限推力, 连续推力 ) 轨道机动任务包括轨道转移 交会和拦 9

28 1.2. 研究现状 截 最优轨道转移是研究的热点, 包括地球卫星轨道间转移 地球轨道 - 行星轨 道转移以及星际转移任务 文献 [15, 17 2, 57] 应用位置 - 速度描述的动力学方程 研究了地球卫星轨道转移问题, 而文献 [3 32, 36, 37, 6 63] 则用改进轨道根数 描述动力学方程研究了轨道优化问题 地球 - 行星轨道转移主要是研究地球 - 火 星轨道转移任务 [39, 59, 94, 95], 地球轨道 - 月球轨道转移任务也有研究 [96], 而 文献 [48, 97, 98] 则研究了星际航行问题 不同轨道机动任务的区别主要在终端 约束边值条件, 如用位置 - 速度描述, 在终端时刻位置和速度向量均确定则是轨 道交会问题, 而只有位置向量确定则是轨道拦截问题 因而, 不同的终端约束 条件可以表示不同的轨道机动形式, 如文献 [19,57] 在应用位置 - 速度描述的动力 学方程研究轨道转移任务时就是用角动量和偏心率向量来确定终端条件, 而文 献 [65] 则考虑了多种不同的终端约束条件来描述轨道转移任务 而关于轨道交会和拦截任务, 早期应用脉冲法开展了大量的研究 [3] 小推 力控制情形时, 相对轨道动力学方程广泛应用于轨道交会任务 [99, 1], 罗亚 中应用遗传算法 模拟退火算法等智能算法对相对动力学方程描述的轨道交会 问题开展研究 [11 13], 而文献 [14] 则用相对轨道根数推导出相对动力学方 程并将其应用于轨道交会问题中 Kechichian 发展了一套改进的轨道根数来描述 轨道运动, 应用间接法研究了最优轨道交会问题 [35, 15] 文献 [16] 在太阳系 坐标系统中用开普勒型方程研究了能量最优轨道交会问题, 考虑能量受限和比 冲受限, 并应用极大值原理设计了最优控制律 Park 则在太阳系坐标系统中应 用球坐标描述的动力学方程研究了时间 - 燃料混合最优轨道拦截和交会问题, 不 同的终端约束条件表示了轨道拦截和交会任务, 并应用横截条件推导出边值约 束 [98] 轨道优化的性能指标主要是要求飞行时间最短或者消耗燃料最小, 时间最 优轨道转移问题得到了广泛深入的研究 [3, 31, 33, 61 63, 17, 18] 而对于燃料 最优问题, 通常性能指标要求是初始质量最小或者终端质量最大, 对于常值 推力发动机, 性能指标则可描述为 J = t f u dt, 由于应用间接法研究燃料最优问题时存在奇异性, 通常在性能指标中引入能量指标 [32, 6] 或者其他惩罚 函数 [19], 在数值上再应用连续方法如同伦方法来研究 [32, 6] 也有文献将 时间和质量最优统一来考虑, 两者之间存在权重系数 [98], 但是由于量纲不统 一, 权重系数不易确定, 作者认为这种性能指标不是很合理, 文献也较少应用 这种指标形式 1

29 第一章绪 论 1.3 研究内容与结构安排本论文主要应用间接法研究轨道转移 交会和拦截问题, 考虑航天器为点质量模型和常值推力发动机, 分别对位置 - 速度和 MEE 描述的轨道动力学方程应用极大值原理设计时间最优控制律, 并进行了奇异性分析, 根据极大值原理一阶必要条件, 将轨迹优化问题转化为两点边值问题, 然后在数值上应用打靶法来求解相应的边值问题 将最终累积赤经作为新的系统自变量, 借鉴平均化方法研究了累积赤经最优轨道转移问题 对于能量最优轨道转移问题, 综合考虑燃料 能量以及能量 - 燃料混合最优问题, 应用极大值原理设计最优控制律, 数值上应用打靶法和 ε 算法来求解 论文还研究了最优轨道规避问题, 假设存在中间过渡轨道, 应用非线性规划方法研究了时间最优和燃料最优轨道规避问题 根据上述简单描述的研究内容, 论文在结构上安排如下 : 第二章, 首先讨论了航天器轨道动力学基础, 在二体轨道理论框架内, 给出了用位置 - 速度描述的开普勒型和用轨道根数描述的高斯型轨道动力学方程, 并简单讨论了球坐标描述的方程和航天器质量方程, 为后续章节的研究奠定基础 第三章, 在惯性坐标系和轨道坐标系下分别应用开普勒型方程研究时间最优轨道拦截问题, 应用极大值原理设计最优控制律并进行了奇异性分析, 将轨道拦截问题转化为两点边值问题来求解 数值仿真部分给出两个不同终端约束的轨道拦截任务 第四章则应用改进轨道根数描述的高斯型动力学方程, 在统一的框架下研究时间最优轨道转移 交会和拦截问题, 应用间接法设计了最优控制律, 极大值原理一阶条件将轨迹优化问题转化为两点边值问题, 然后应用单值打靶法求解相应的边值问题 第五章, 将累积赤经引入系统作为新的变自量, 应用间接法研究考虑变质量和具推力约束下的最优累积赤经轨道转移问题 ; 借鉴平均化方法, 设计了近似最优控制律, 同样将最优轨道转移问题转化为两点边值问题 第六章, 基于 MEE 描述轨道动力学方程, 应用间接法研究燃料最优 能量最优以及能量 - 燃料混合最优轨道转移, 分别设计了最优控制律, 尝试在数值上对这三种形式的最优轨道转移做比较并得出结论 第七章研究常值连续推力时间最优和燃料最优规避问题 在给定中间过渡 11

30 1.3. 研究内容与结构安排 轨道条件下, 将规避问题划分为转移和返回两个阶段, 然后应用非线性规划方法来研究时间和燃料最优规避问题 第八章则给出全文的结论, 并对未来的继续研究方向做了展望 在附录部分, 也简要介绍了极大值原理 两点边值问题和打靶方法 12

31 第二章航天器轨道动力学基础 第二章 航天器轨道动力学基础 这一章, 我们主要讨论航天器轨道动力学, 为后续章节的研究奠定基础 首先讨论二体问题及描述轨道运动的轨道根数, 然后讨论常用的三种描述连续推力轨道机动问题的轨道动力学方程, 即 : 一是在笛卡尔坐标系下用位置和速度表示的开普勒型方程, 二是在轨道坐标系下用改进赤道轨道根数表示的高斯方程, 三是在球坐标系下用球半径 幅角以及速度表示的动力学方程 2.1 二体问题及轨道根数轨道动力学, 是把航天器当做质点来研究其在外层空间随时间变化过程 轨道 在航天器轨道的分析问题中, 常假设其在地球中心力场中运动, 忽略其他各种摄动 ( 如地球非球形 密度分布不均匀引起的摄动, 太阳 月球引力摄动, 大气压力摄动等 ) 这种研究一个质点在另一个质点( 本论文分别指卫星和地球 ) 的万有引力作用下运动的问题, 即是二体问题, 这种卫星轨道则称为二体轨道, 它是轨道动力学的基础 运用牛顿关于物体线动量的变化率等于作用力 F = d dt (mv) 这一定律可以导出轨道动力学的基本方程 [1, 2, 14, 29, 11] r = µ r 2 r r (2.1) 其中 r 为卫星相对地球地心位置矢量,r = r 为其矢量幅值,µ = Gm e 为地 球引力常数 (G 为万有引力常数,m e 为地球质量, 且假设卫星质量远远小于地 球质量 ) 方程 2.1 称之为开普勒方程, 由此可以看出, 作用在卫星上的地球 中心引力仅与卫星质量成正比, 与卫星的地心距的平方成反比 ; 由引力产生的 加速度幅值 反 [29] µ r 2 与卫星质量无关, 引力加速度幅值与地心距单位矢量 r r 的方向相 动力学方程 2.1 可化成二体系统的基本动力学微分方程, 它是一个六阶的非 线性常微分方程组, 若要完全求解该方程组则必须找出含有六个相互独立的积 分常数的解 常选择一组意义明确且相互独立的六个积分常数代表运动轨道的 基本量, 称之为轨道根数, 用 σ 表示 13

32 2.1. 二体问题及轨道根数 对于椭圆轨道, 常用的经典轨道根数为 σ = (a, e, i, Ω, ω, θ) [1, 14]: a 椭圆轨道半长轴 e 偏心率 i 轨道倾角 Ω 升交点赤经 ω 近地点幅角 θ 真近地点角 其中 a 和 e 表示轨道的大小和形状,Ω 和 i 表示轨道面在空间的指向,ω 表示在轨 道面内近地点方向的指向,θ 表示过近地点时间, 也即表示卫星在轨道上的位 置 [1], 其间关系由图 2.1 简单所示 经典的轨道根数与航天器在地心惯性坐标系下的位置和速度矢量是可以相 互转换的, 用经典的轨道根数描述的位置和速度矢量为 [14] r = v = P 1+e cos θ P 1+e sin θ (cos Ω cos(ω + θ) sin Ω sin(ω + θ) cos i) (cos Ω cos(ω + θ) + cos Ω sin(ω + θ) cos i) P sin(ω + θ) sin i 1+e cos θ µ (cos Ω(sin θ + e sin ω) + sin Ω(cos θ + e cos ω) cos i) P µ (sin Ω(sin θ + e sin ω) cos Ω(cos θ + e cos ω) cos i) P 其中 P = a(1 e 2 ) 表示轨道的半通径 µ (cos θ + e cos ω) sin i P (2.2) 描述轨道运动的轨道根数有很多种形式, 如改进的经典轨道根数 庞加 莱根数 欧拉参数等等 [111] 由于经典轨道根数描述的动力学方程在偏心 率 e = 和 ( 或 ) 轨道倾角 i = 时存在奇异,Walker 提出了一种称为 Modified Equinoctial Elements(MEE) 的轨道根数 [27], 以该根数为状态的轨道动力学方 程得以广泛应用, 可以适用于除 i = 18 以外的任意圆锥曲线轨道 该轨道根 数与经典的轨道根数的关系为 14

33 第二章航天器轨道动力学基础 图 2.1 轨道根数与坐标系简单示意图 P = a(1 e 2 ) e x = e cos(ω + Ω) e y = e sin(ω + Ω) h x = tan( i ) cos Ω 2 h y = tan( i ) sin Ω 2 L = Ω + ω + θ (2.3) 其中第一个分量确定轨道的大小, 第二个和第三个分量确定轨道的形状 ( 圆, 椭圆, 抛物线, 双曲线 ), 第四个和第五个分量确定轨道倾角, 第六个分量确定卫星在轨道上的位置 前五个分量可看作慢变量, 第六个则是快变量 第一个分量也可考虑半长轴 a, 能量, 角动量 h = P/µ 等, 第四 五个分量也可考虑为 i cos Ω, i sin Ω 等, 第六个也可考虑偏赤经 Ω + ω + E 或平赤经 Ω + ω + M 同样, 改进的轨道根数 2.3 与惯性坐标系下的位置和速度分量的关系 15

34 2.2. 开普勒型方程 为 [14, 112] r 1 = P CW ((1 + h2 x h 2 y) cos L + 2h x h y sin L) r 2 = P CW ((1 h2 x + h 2 y) sin L + 2h x h y cos L) r 3 = P CW (2Z) v 1 = 1 µ C P (2h xh y (e x + cos L) (1 + h 2 x h 2 y)(e y + sin L)) v 2 = 1 µ C P ((1 h2 x + h 2 y)(e x + cos L) 2h x h y (e y + sin L)) v 3 = 1 µ C P (2h x(e x + cos L) + 2h y (e y + sin L)) (2.4a) (2.4b) (2.4c) (2.4d) (2.4e) (2.4f) 其中 (r 1, r 2, r 3 ), (v 1, v 2, v 3 ) 分别表示位置和速度矢量分量, 其他参数为 C = 1 + h 2 x + h 2 y, W = 1 + e x cos L + e y sin L 本文主要应用位置 速度变量和改进的赤道轨道根数描述的动力学方程, 接下来我们分别考虑其描述的系统方程 2.2 开普勒型方程考虑以位置 速度为状态变量的动力学方程, 在不考虑姿态时假设卫星为点质量, 在外力作用下的动力学方程为 : r = µ r r 3 + T m + T p (2.5) 其中 r 为卫星相对地球地心位置矢量,m 为卫星质量,µ 为地球引力常数 ;T 为发动机推力 ( 控制输入 ) 矢量, 并且是有界的满足 T T max ;T p 表示摄动力矢量 当 T = T p = 时方程 2.5 即为开普勒方程 2.1, 因此我们将以位置和速度为变量描述的方程称之为开普勒型方程 在本节, 我们不考虑摄动力的影响, 为书写方便, 以下省略 T p 控制在地心坐标系下的动力学方程 首先考虑控制输入在地心坐标系下分解的动力学方程 惯性坐标 系 (X,Y,Z) 由图 2.1 所示, 定义为 :X 轴方向从地心指向春分点,Y 轴方向位于 16

35 第二章航天器轨道动力学基础 赤道平面内与于 Y 方向正交,Z 为赤道面法向方向,(X,Y,Z) 构成右手直角坐标系 卫星在该坐标系下位置和速度矢量为 ] T r = [r 1 r 2 r 3 ] T v = [v 1 v 2 v 3 r = r = 2 r r r 2 3 则卫星的动力学方程为 ṙ = v v = µ r r 3 + u m (2.6) 其中控制约束需满足 u T max 控制在轨道坐标系下的动力学方程 上面我们考虑了推力在惯性坐标系下分解的动力学方程, 通常情况下, 发动机推力方向随轨道坐标系变化 在轨道坐标系下, 推力方向主要有两种分解形式 : 切向 - 法向分解和径向 - 法向分解 [29, 3] 切向- 法向标架 (T,C,N) 定义为, 在轨道坐标系中,T 方向沿轨道切向方向,C 在轨道面内与 T 方向垂直,N 方向为轨道面瞬时法向平行, 构成右手直角坐标系 ; 而径向 - 法向标架 (R,T,N) 如图 2.1 所示, 定义为, 在轨道坐标系中,R 方向沿卫星与地球连线径向方向,T 在轨道面内与 R 方向垂直,N 方向与瞬时轨道面法向平行, 构成右手直角坐标系 接下来我们在轨道坐标系下研究轨道动力学方程, 设 : 航天器的速度空间为 T M R 3 (M 为构型空间 ), 我们首先定义切丛空间 T R 3 的一组基向量, [ r ] T [ v = r 1 r 2 r 3 v 1 v 2 ] T v 3 则状态变量可写作 : [ ] ṙ = ṙ 1 ṙ 2 ṙ 3 r = ṙ r = v r [ ] v = v 1 v 2 v 3 v = v v 17

36 2.2. 开普勒型方程 令 由轨道动力学方程 2.5 则有 x = ẋ = ẋ = [ r [ ṙ [ ṙ ] v T R 3 ] ] [ r v ] v v ẋ = = [ ṙ [ v ] [ ] r v v µ r r 3 + u m = v r µ r r 3 ] [ ] r v v + u m v 令 则有轨道动力学方程 f = v r µ r r 3 v f = u m v ẋ = f + 1 m f a)(t,c,n) 型方程 假定推力不为零, 推力方向则可在切向 - 法向标架下分解为 f = f v = u tf t v + u nf n v + u cf c v = u tf t + u n f n + u c f c (2.7) 其中 : f t = v v f n = r v r v f c = f c f t = (r v) v (r v) v 18

37 第二章航天器轨道动力学基础 其中符号 表示向量的外积, 定义运算为 : r v = r 3 r 2 r 3 r 1 r 2 r 1 v 1 v 2 v 3 T = r 2 v 3 r 3 v 2 r 3 v 1 r 1 v 3 r 1 v 2 r 2 v 1 T 根据上述定义可知 (r v) v = r, v v v, v r = (r T v)v 1 v 2 r 1 (r T v)v 2 v 2 r 2 (r T v)v 3 v 2 r 3 T 因此方程 2.7 中描述运动动态的各分量为 f = 1 [ ] r 3 v 1 r 3 v 2 r 3 v 3 r 3 r 1 µ r 2 µ r 3 µ f t = 1 [ ] v v 1 v 2 v 3 1 [ ] f n = r r v 2 v 3 r 3 v 2 r 3 v 1 r 1 v 3 r 1 v 2 r 2 v 1 1 [ ] f c = (rt v)v v 2 r (r T v)v 1 v 2 r 1 (r T v)v 2 v 2 r 2 (r T v)v 3 v 2 r 3 则轨道动力学方程写成分量形式有 : ṙ 1 = v 1 ṙ 2 = v 2 ṙ 3 = v 3 v 1 = µr 1 r 3 + v 1 m v u t + (rt v)v 1 v 2 r 1 m (rt v)v v 2 r u c + r 2v 3 r 3 v 2 m r v u n v 2 = µr 2 r 3 + v 2 m v u t + (rt v)v 2 v 2 r 2 m (rt v)v v 2 r u c + r 3v 1 r 1 v 3 m r v u n v 3 = µr 3 r 3 + v 3 m v u t + (rt v)v 3 v 2 r 3 m (rt v)v v 2 r u c + r 1v 2 r 2 v 1 m r v u n 19

38 2.3. 轨道根数描述高斯方程 b)(r,t,n) 型方程 另一种重要分解为径向 - 法向分解 f = u r f r + u t f t + u n f n 其中 f r = r r v f n = r v r v v f t = f c f r f c f r 按照上述同样方法可得 : (r v) r = r, v r r, r v (r T v)r 1 r 2 T v 1 = (r T v)r 2 r 2 v 2 (r T v)r 3 r 2 v 3 则轨道动力学方程写成分量形式为 ṙ 1 = v 1 ṙ 2 = v 2 ṙ 3 = v 3 v 1 = µr 1 r 3 + r 1 m r u r + (rt v)r 1 r 2 v 1 m (rt v)r r 2 v u t + r 2v 3 r 3 v 2 m r v u n (2.8) v 2 = µr 2 r 3 + r 2 m r u r + (rt v)r 2 r 2 v 2 m (r T v)r r 2 v u t + r 3v 1 r 1 v 3 m r v u n v 3 = µr 3 qr 3 + r 3 m r u r + (rt v)r 3 r 2 v 3 m (rt v)r r 2 v u t + r 1v 2 r 2 v 1 m r v u n 2.3 轨道根数描述高斯方程 2 由于用位置 - 速度描述的轨道动力学方程具有振荡性, 因而常用轨道根数

39 第二章航天器轨道动力学基础 来描述卫星的运动 [32, 62] 实际应用中, 卫星轨道常采用六个轨道根数来描述 :a - 椭圆轨道半长轴 ;e- 偏心率 ;i- 轨道倾角 ;ω- 近地点幅角 ;Ω- 升交点赤经 ;M- 平近地点角 若航天器的作用力为保守力, 比如各种摄动力, 可以利用常数变易原理导出拉格朗日 (Lagrange) 型的摄动运动方程 [29] 而对于非保守力, 作用力在轨道坐标系 (R,T,N) 标架下分解, 则用经典轨道根数描述的高斯型轨道动力学方程为 [2] da dt = 2 n 1 e (e sin θu 2 r + (1 + e cos θ)u t ) de 1 e dt = 2 (sin θu r + (cos E + cos θ)u t ) na dω r sin (ω + θ) = dt na 2 1 e 2 sin i u n di r cos (ω + θ) = dt na 2 1 e u 2 n ( dω 1 e dt = 2 cos θu r e cos θ ) nae 1 + e cos θ sin θu t cos i dω dt dm = n 1 ( e2 ( 2er dt nae P cos θ)u r + (1 + r ) P ) sin θu t (2.9a) (2.9b) (2.9c) (2.9d) (2.9e) (2.9f) 方程 2.9 中,n = µ a 3 为真近点角的平均角速率,r = a(1 e2 ) 1+e cos θ 为地心距,E 为偏近点角, 偏近点角与真近点角的关系满足 而平近点角与偏近点角的关系满足 tan( θ ( ) 1/2 1 + e 2 ) = tan( E 1 e 2 ) M = E e sin E 由于经典轨道六根数描述的高斯行星方程 2.9 在偏心率 e = 和 ( 或 ) 轨道倾角 i = 时存在奇异, 为消除奇异性, 以改进的赤道轨道根数 2.3 为状态的轨道运动方程得以更加广泛的应用, 可以适用于除 i = 18 以外的任意圆锥曲线轨道 21

40 2.3. 轨道根数描述高斯方程 下面考虑改进轨道根数描述的动力学方程, 令状态为 [ x = a e x e y h x h y L ] ẋ = ȧ a + ė x + ė y + e x e ḣx + y h ḣy + x h L y L = ẋ x 则该系统在切向 - 法向标架下分解为 [3] ẋ = f + 1 m (u tf t + u c f c + u n f n ) 式中 : µ W 2 [ f = P P L = f t = 1 P W µ = 1 P [ 2W P η W µ (1 e 2 ) 2 f c = 1 ( P Dηx W η y W µ η = 1 P [ Dη x W η y W µ η f n = 1 W = 1 W P µ P [ µ µ W 2 P P ( 2W P η (1 e 2 ) 2 a + 2W η x η 2W η x η + 2W η y e x η 2W η y η + Dη x + W η y e x η ] x = f x ) e y ] x = f t x ) e y (2.1a) (2.1b) (2.1c) ] Dη x +W η y η x = f c x ( Ze y + Ze x + C cos L + C sin L + Z ) e x e y 2 h x 2 h y L Ze y Ze x C cos L 2 C sin L 2 Z ] x = f n x (2.1d) 22

41 第二章航天器轨道动力学基础 其中 P = a(1 e 2 ) η = (η x, η y ) = (e x + cos L, e y + sin L) W = 1 + e x cos l + e y sin l D = e x sin L e y cos L C = 1 + h 2 x + h 2 y Z = h x sin L h y cos L 按分量表示, 则描述轨道动力学方程为 da dt = P µ de x dt = 1 W de y dt = 1 W dh x dt = 1 W dh y dt = 1 W dl dt = 2P η (1 e 2 ) u 2 t ( 2W ηx P µ P µ P µ P µ η u t + Dη ) x W η y u c Ze y u n η ( 2W ηy η u t + Dη ) x + W η y u c + Ze x u n η C cos L u n 2 C sin L u n 2 µ W 2 P P + Zu n (2.11) 上述表式中应用的改进轨道根数为 (a, e x, e y, h x, h y, L), 现将第一个根数改为半 通径 P, 推力矢量在 (R,T,N) 轨道坐标系下分解, 轨道动力学方程可描述为 : ẋ = f + 1 m (u rf r + u t f t + u n f n ) 23

42 2.3. 轨道根数描述高斯方程 式中描述卫星运动的四个动态为 : µ W 2 [ f = P P L = ] µ W 2 P P x = f x f r = 1 ( P W sin L W cos L ) W µ e x e y = 1 P [ ] W sin L W cos L W µ x = f r x f t = 1 ( P 2P W µ P + (W cos L + η x) + (W sin L + η y ) ) e x e y = 1 P [ ] 2P (W cos L + η W µ x ) (W sin L + η y ) x = f t f n = 1 ( P Ze y + Ze x + C cos l + C sin L + Z W µ e x e y 2 h x 2 h y L = 1 W P [ µ Ze y Ze x C cos ll 2 C sin L 2 Z ] x = f n x x ) (2.12) 其中 W = 1 + e x cos L + e y sin L, Z = h x sin L h y cos L, η x = e x + cos L, η y = e y + sin L 和 C = 1 + h 2 x + h 2 y 动态 2.12 描述的轨道动力学方程, 漂移项 f 是周期的, 且在切空间内, 向 量 f, f r, f t, f n 张成的李括号 lie x (f, f r, f t, f n ) = T x M, M 为构型流形 因此, 根据仿射输入非线性系统可控性判断条件 [64], 系统保持可控性 本论文主要应用改进轨道根数描述的方程 2.12 来研究轨道机动问题 球坐 标系统 ( 也称为极坐标系统 ) 也被广泛应用于轨道机动, 特别是深空探测使 命, 如文献 [98] 应用球坐标方程研究了在太阳系系统中对那些对地球造成碰撞 威胁的流星或彗星实施拦截问题 下面简要给出极坐标描述的轨道运动方程 24

43 第二章航天器轨道动力学基础 图 2.2 球坐标系统示意图 2.4 极坐标方程 球坐标系统如图 2.2 所示, 运动方程则可描述为 [2, 98] ṙ = v r θ = v θ r cos φ φ = v φ r v r = v2 θ r + v2 φ r µ r + f T a r v θ = v rv θ r v φ = v rv φ r + v φv θ r tan φ + f T a θ v2 θ r tan φ + f T a φ (2.13) 其中 r 为航天器距坐标系系统中心的半径幅值,θ 和 φ 为两个位置角度, 如 图 2.2 所示,v r, v θ, v φ 则定义三个速度分量,f T 为发动机推力幅值,a r, a θ, a φ 表 25

44 2.5. 卫星质量方程 示推力方向在轨道坐标系下的分量,a r 沿半径方向,a θ 垂直于 a r 且平行于 XY 平面,a φ 则满足右手定则 2.5 卫星质量方程上面研究了卫星轨道动力学方程的三种描述形式, 通常情况下假设卫星为点质量, 并未考虑卫星质量变化 实际上, 对于轨道机动任务, 卫星携带的燃料是有限的, 燃料是受限的, 因此必须考虑卫星质量变化 对于连续小推力发动机, 本论文并不具体的去考虑采取何种推进方式, 若发动机推力为常值, 即比冲固定时, 卫星质量变化可描述为 [3] ṁ(t) = T. = βt (2.14) I sp g 其中 I sp 为发动机推力比冲,g 为当地重力加速度,T 为固定比冲时发动机推力且 T T max 若考虑变比冲发动机, 推力 T 和排气速度 v e = I s pg 的关系为 T = ṁ v e 而能量与排气速度的关系则为 [98] εp = 1 2 ṁ v2 e 其中 ε 为推进系统的效率 因此能量和推力的关系为 则卫星质量变化方程为 T = 2εP I sp g ṁ(t) = T. T 2 = I sp g 2εP. = 2εP (I sp g ) 小结 26 本章主要讨论了二体轨道问题的轨道根数和轨道动力学方程 在平方反比

45 第二章航天器轨道动力学基础 地球引力场中, 考虑有控外力作用, 常用位置 - 速度矢量和轨道根数作为变量, 描述航天器的运动 以位置 - 速度分量作为变量, 分别讨论了在地球惯性坐标系和轨道坐标系下的有控动力学方程 ; 应用改进的轨道根数, 分别讨论了轨道坐标系中常用两种分解方式下的动力学方程 ; 简单讨论了极坐标描述的动力学方程和卫星质量变化方程 在第三章, 将应用开普勒型方程, 考虑卫星质量变化, 研究时间最优轨道机动问题 27

46

47 第三章时间最优轨道机动 Ⅰ: 开普勒型方程 第三章 时间最优轨道机动 Ⅰ: 开普勒型方程 有限推力时间最优轨道机动问题是航空航天领域的研究热点, 轨道机动模式包括轨道转移 拦截与交会, 目前众多的研究集中在轨道转移方面 轨道动力学的基础是二体轨道动力学 [2, 14], 在地球平方反比力场中, 有控模型满足如下开普勒型方程 2.5 r = µ r r + T 3 m + T p (3.1) 上式描述的轨道动力学方程的状态变量是位置和速度矢量分量, 被广泛应用于轨道优化问题研究 文献 [57] 研究了燃料最优轨道转移问题, 考虑了 J2 摄动和大气摄动, 应用间接法将问题化成两点边值问题, 然后应用多点打靶法 BOUNDSCO 来求解 ; 而文献 [91] 则将直接法和间接法结合使用, 应用遗传算法求解相应的边值问题 文献 [15] 应用开普勒型方程研究了时间最优轨道转移问题, 应用极大值原理设计了最优控制律, 并提出了一种切换函数检验和分析的方法 不同轨道机动任务的区别主要在终端约束边值条件, 在终端时刻位置和速度向量均确定则是轨道交会问题, 而只有位置向量确定则是轨道拦截问题 因而, 不同的终端约束条件可以表示不同的轨道机动形式, 如文献 [19, 57] 用位置 - 速度表示的角动量和偏心率向量来确定终端条件 本章应用位置 - 速度描述的开普勒型方程研究时间最优轨道机动问题, 由于不同的轨道使命的区别主要在终端约束条件, 而且应用该方程研究轨道转移的文献较多, 本章则主要研究时间最优轨道拦截问题 另外, 已有的文献应用开普勒型方程, 控制输入都是在惯性坐标系下分解, 由于在轨道坐标系中, 更容易分析控制分量对轨道参数改变的影响, 因此本章主要应用轨道坐标系下描述的开普勒型方程研究轨道拦截问题 本章的主要贡献有两点 :1) 应用控制在轨道坐标系下分解的动力学模型, 根据极大值原理设计最优控制律, 分析了奇异性, 将时间最优轨道拦截问题化成两点边值问题 ;2) 考虑了一种新的轨道拦截问题的终端条件描述, 将目标轨道参数引入系统, 使轨道拦截问题更合理 3.1 时间最优轨道机动问题描述在二体轨道动力学基础上, 航天器 ( 本文主要是指地球卫星 ) 动力学方程可 29

48 3.1. 时间最优轨道机动问题描述 写成一仿射输入非线性系统 ẋ = f (x) + T max m 3 i=1 u if i (x) (3.2) 考虑常值推力固定比冲发动机, 其中 T max 为发动机最大推力幅值 同时, 卫星 质量变化满足方程 而推力方向作为控制则满足约束 u = ṁ = βt max u (3.3) u u u (3.4) 对于特定的轨道机动任务, 仍需满足初始和终端边值条件 Φ (t, x(t ), m(t )) = (3.5) Φ 1 (t f, x(t f ), m(t f )) = (3.6) 对于时间最优问题, 则需最小化性能指标 J = tf t dt (3.7) 方程 构成了时间最优轨道机动问题, 为简化研究问题, 我们考虑以 下假设条件 A1) 系统模型的状态始终满足路径约束 A = {(x, m) r > R e, m m s } 即卫星在椭圆域内飞行, 在地心坐标系下位置向量幅值始终大于地球半径 R e,m s 为不携带燃料时的卫星本体质量 A2) 初始轨道和目标轨道是相异的, 即最优控制是非空的 A3) 卫星在终端时质量是自由的 3

49 第三章时间最优轨道机动 Ⅰ: 开普勒型方程 3.2 最优轨迹设计 考虑以位置 速度为状态变量的动力学方程, 在不考虑姿态影响时假设卫 星为点质量体, 在外力作用下 ( 不包括其他摄动力 ) 的动力学方程为 : r = µ r r + T 3 m (3.8) 其中 r 为卫星相对地球地心位置矢量,m 为卫星质量,µ 为地球引力常数 ;T 为 发动机推力 ( 控制输入 ), 需满足约束 T T max 地心坐标系描述情形 首先考虑控制在地心惯性坐标系 XYZ 下分解情形, 令 T = T max u, 其中常数 T max > 为最大推力幅值, 控制约束变为 u 1, 则考虑卫星质量变化的的动力学状态方程为 ṙ = v v = µ r r + T max 3 m u (3.9) ṁ = βt max u 其中 x = (r, v, m) R 3 R 3 R + 为系统状态 考虑时间最优轨道转移, 问题满足的初始条件为 x(t ) = (r(t ), v(t ), m(t )) (3.1) 而对于轨道拦截问题, 则在终端时刻 r(t f ) 已知,v(t f ) 自由,m(t f ) m e 利用 Pontryagin 极大值原理, 得系统的 Hamiltonian 函数为 H = 1 + p T r v + p T v( r µ r 3 + T max m u) p mβt u (3.11) 其中 p = (p r, p v, p m ) R 3 R 3 R 为系统状态 (r, v, m) 对应的协态变量 当 p v 时, 对应的最优解 x = (r, v, m ) 称为正则最优解, 对应的控制 31

50 3.2. 最优轨迹设计 称为正则最优控制, 这时使 H 达极大值的最优推力取为 : u = p v p v (3.12) 上述分析中, 我们假定 p v, 若在某时间点 t i [t, t f ], i = 1, 2,..., N 有 p v =, 而 p v 时, 有 H u = pt vt max m = 在这样的点控制 u 将改变方向, 进行 18 度反向切换 ; 而当在一个测度不为 的时间段 [t 1, t 2 ] [t, t f ] 上 p v = 时, 相应的最优解 x = (r, v, m ) 称为奇异最优解, 对应的控制称为奇异最优控制, 由于此时有 ṗ v =, 而选择这种坐标系的情况下 ṗ v = p r, 这将导致在一个测度不为 的时间段 [t 1, t 2 ] [t, t f ] 上协态变量 p, 根据最大值原理, 这种情况是不可能的 所以在这种坐标系的选择下, 不存在奇异最优解和奇异控制 因此, u = 1 几乎处处满足, 系统的最优 Hamiltonian 函数为 : H = 1 + p T r v µ r 3 pt vr + ( p v m βp m)t max (3.13) 根据极大值原理一阶必要条件得正则状态方程 ṙ = v v = µ r r 3 + T max m u (3.14) ṁ = βt max 和协态方程 ṗ r = µ( p v r 3 3pT vr r 5 r) ṗ v = p r ṗ m = T m 2 p v (3.15) 32

51 第三章时间最优轨道机动 Ⅰ: 开普勒型方程 考虑横截条件, 则在终端时刻有 : p r (t f ) = [λ r1 λ r2 λ r3 ] T p v (t f ) = {, m(tf ) > m e p m (t f ) = λ m, m(t f ) = m e (3.16) H (t f ) = 1 + p T r (t f )v(t f ) βp m (t f )T max = 其中 λ i 为对应于终端约束的拉格朗日乘子 针对我们研究的问题, 假设当轨道 转移完成时, 燃料并未消耗完, 即满足假设条件 A3, 则横截条件取为 : p v (t f ) = p m (t f ) = (3.17) 而在终端时刻位置同时满足 r(t f ) = [r f 1 r f 2 r f 3 ] T 因此将地心惯性坐标下描述的轨道拦截问题转化为两点边值问题 轨道坐标系描述情形 考虑推力方向在轨道坐标下沿切向 - 法向 (R,T,N) 分解, 则有 u = u 1 f 1 + u 2 f 2 + u 3 f 3 33

52 3.2. 最优轨迹设计 其中 f 1, f 2, f 3 为轨道坐标系的基向量 : f 1 = 1 v f 2 = f 3 = v 1 v 2 v 3 1 (rt v)v v 2 r 1 r v r 2 v 3 r 3 v 2 r 3 v 1 r 1 v 3 r 1 v 2 r 2 v 1 (r T v)v 1 v 2 r 1 (r T v)v 2 v 2 r 2 (r T v)v 3 v 2 r 3 (3.18) 记 u = [u 1 u 2 u 3 ] T, 则 u 满足 u 1, 系统的状态方程为 : ṙ = v v = µ r r 3 + T max m (u 1f 1 + u 2 f 2 + u 3 f 3 ) ṁ = βt max u (3.19) 利用 Pontryagin 极大值原理, 得系统的 Hamiltonian 函数为 其中 H = 1 + H + 3 i=1 u ih i H = p T r v µ r 3 pt vr βp m T max u H i = T max m pt vf i, i = 1, 2, 3 (3.2) 令 Ψ = [H 1, H 2, H 3 ], 当 Ψ, 相应的最优解 x = (r, v, m ) 称为正则最优 解, 对应的控制称为正则最优控制, 这时使 H 达极大值的最优推力取为 : u = Ψ Ψ (3.21) 若假设不存在奇异性, u = 1 几乎处处满足, 则可得正则情况下系统的最 34

53 第三章时间最优轨道机动 Ⅰ: 开普勒型方程 优 Hamiltonian 函数为 : H = 1 + H + Ψ 对应的状态方程 ṙ = v 和协状态方程 v = µ r r 3 + ṁ = βt max u 同样考虑横截条件, 则在终端时刻有 : T 3 max i=1 pt vf i u i m( 3 i=1 (pt vf i u i ) 2 ) 1 2 ṗ r = H r ṗ v = H v ṗ m = T max m 2 Ψ (3.22) (3.23) p v (t f ) = p m (t f ) = (3.24) 终端时刻位置同时满足 r(t f ) = [r f 1 r f 2 r f 3 ] T, 因此将 (R,T,N) 轨道坐标系下描述 的轨道拦截问题转化为两点边值问题 在上面的分析中, 我们假定 Ψ, 现进行奇异性分析 若在某时间 点 t i [t, t f ], i = 1, 2,..., N 有 Ψ =, 而 Ψ 时, 有 H u = 在这样的点控制 u 将改变方向, 进行 18 度反向切换 ; 而当在一个测度不为 的 时间段 [t 1, t 2 ] [t, t f ] 上 Ψ = 时, 相应的最优解 x = (r, v, m ) 称为奇异最优 解, 对应的控制称为奇异最优控制 由坐标系的选取, 易知此时有 p v (t) =, 而 由 Ψ =, 可得方程 Ḣ i = {H i, H} = {H i, H } + 3 j i {H i, H j }u j =, i = 1, 2, 3 35

54 3.3. 数值仿真 或者 h 12 h 13 h 12 h 23 h 13 h 23 u 1 u 2 u 3 = h 1 h 2 h 3 其中 h ij 是由下式定义的 Poisson 括号 h ij = {H i, H j } = 3 k=1 H i v k H j p vk H i p vk H j v k, i =, 1, 2, 3, j = 1, 2, 3 当系数矩阵非奇异时, 可解出最优奇异控制 易证 Poisson 括号和李括号有如下关系 {H i, H j } = T max m pt v[f i, f j ] 若向量场 f 1, f 2, f 3 张成的分布 (r, v) = span(f 1, f 2, f 3 ) 是对和的, 则任意 的 [f i, f j ], i, j = 1, 2, 3, 可由 f 1, f 2, f 3 线性表出, 由此可知 {H i, H j } 可表示 成 H 1, H 2, H 3 线性的线性组合, 从而有 Ḣ i = {H i, H } =, i = 1, 2, 3 由此可得 p r =, 与控制在地心坐标系分解的情形相同, 这将导致在一个测度不为 的时间段 [t 1, t 2 ] [t, t f ] 上协态变量 p =, 根据最大值原理, 这种情况是不可能的 所以在这种坐标系的选择下, 当分布 (r, v) = span(f 1, f 2, f 3 ) 对和时不存在奇异最优解和奇异控制 注记 3.1: 在本章, 我们应用极大值原理设计最优控制律 若应用所谓的极小值原理, 则控制律和终端横截条件的符号则相反 [76] 但是, 不管是应用极大值还是极小值原理, 化成的两点边值问题的最优解是相同的 在后续章节中, 应用极大值原理研究轨迹优化问题, 都有类似的情况, 不再另行说明 3.3 数值仿真在 3.2 节, 我们考虑了以卫星位置和速度为变量的动力学方程, 对于时间最优轨道拦截问题, 应用极大值原理, 分别研究了惯性坐标系和轨道坐标系下 36

55 第三章时间最优轨道机动 Ⅰ: 开普勒型方程 最优推力矢量设计, 从而将时间最优轨道拦截问题转化为两点边值问题 在本 节, 则考虑应用单值打靶法求解相应的两点边值问题, 为便于问题的数值处 理, 转而考虑无量纲化系统 设 R e = Km 为地球半径,t e 为卫星绕半 径为 R e 的圆飞行一圈的时间, 则有 µ ω e = R 3 e, te = 2π re = 2π 3 ω e µ 定义无量纲量 : r = 1 R e r, v = t e R e v, t = t R e, µ = t2 e, µ = 4π 2 Re 3 本节的数值仿真采用的是 节轨道坐标系描述方程, 为书写方便, 以下则以 惯性坐标系下系统为例 描述卫星运动的原动力学方程可写为 d r d t = t e dr R e dt = t e v = t e R e v = v R e R e t e d v d t = t2 e dv r e dt = t2 e ( r µ R e r + T max 3 m u) = t2 e r µ r R + t2 e T max e 3 3 R e m u 即系统方程为 : d r d t = v d v d t = r µ r + t2 e T max 3 r e m u dm d t = βt maxt e u (3.25) 数为 对系统 3.25, 同样应用极大值原理设计最优控制律, 相应的 Hamiltonian 函 H = 1 + p r v + p v ( r µ r 3 + t2 e r e T max m u) p mβt max t e u (3.26) 37

56 3.3. 数值仿真 则根据极大值原理一阶必要条件, 系统的协状态方程为 ṗ r = H r ṗ v = H v ṗ m = H m (3.27) 对于时间最优轨道拦截问题, 初始边界条件为 Φ = (r() r, v() v, m() m ) = 而终端边值条件则为 Φ 1 = (r(t f ) r f, p v (t f ), p m (t f )) = 接下来应用单值打靶法求解化成的两点边值问题 算例 1 系统模型中的物理参数为 µ = Mm 3 /h 2 β = h/mm 而数值边界条件分别为 r = ( Mm,, ) v = (,.9969Mm/h,.2193Mm/h) m = 15kg r f = (3.9392Mm, Mm, Mm) 首先考虑最大推力幅值 T max = 1N, 图 3.1 所示为在地球惯性坐标系下的最优轨迹, 上图中的箭头表示发动机推力方向矢量, 左下图是最优轨迹在地球赤道平面上的投影, 而由右下图则可以看出飞行轨道的轨道倾角变化规律 38

57 第三章时间最优轨道机动 Ⅰ: 开普勒型方程 表 3.1 最小飞行时间 T max (N) t fmin (h) 图 3.2 所示为状态变量即位置 速度分量已经质量随时间变化轨迹, 任务完成时间为 h, 其中为图示方便,R = r, V = v 分别表示位置和速度向量的幅值 由图可以看出位置和速度分量随时间是周期振荡 图 3.3 所示的则是最优推力方向随时间变化轨迹, 由图可以看出满足 u = 1, 因而质量变化可以直接写出 m(t) = m βt max t, 由图 3.2 也可看出与时间成线性变化关系 图 3.1 时间最优轨道拦截 3D 轨迹,1N 图 3.5,3.6,3.7 所示的则是最大推力幅值为 T max = 1N 时相应的数值仿真结果, 飞行时间为 h 改变推力幅值大小, 重复数值实验, 飞行时间分别为相应的结果由图 3.4 表示, 可以看出最小飞行时间和最大推力幅值之间存在如下关系 t f min T max c 其中 c 为常数 39

58 3.3. 数值仿真 r 1 (Mm) r 2 (Mm) r 3 (Mm) t(h) t(h) t(h) R(Mm) t(h) 15 v 1 (Mm/h) v 2 (Mm/h) v 3 (Mm/h) t(h) t(h) t(h) V(Mm/h) t(h) m(kg) t(h) 图 3.2 时间最优轨道拦截状态变化轨迹,1N 算例 2 算例 1 考虑的只是单颗卫星, 给定初始和终端条件, 严格意义上讲并不是轨 道拦截问题, 仍是轨道转移问题 对于轨道拦截问题, 则需研究两个航天器的 相对运动, 终端约束条件不是一个确定值, 而是一个流形约束, 即 ϕ 1 (t f ) = r 1 (t f ) r f 1 = ϕ 2 (t f ) = r 2 (t f ) r f 2 = ϕ 3 (t f ) = r 3 (t f ) r f 3 = (3.28) 即要求在终端时刻两个航天器相对位置为, 而相对速度自由的 式 3.28 中 r f = (r f 1, r f 2, r f 3 ) 则需满足目标轨道约束, 是关于真近点角的变化的函数, 与飞行时间有关 现考虑用经典轨道根数 (a f, e f, i f, Ω f, ω f, θ f ) 描述目标轨道, 前五个参数描述目标轨道的大小 形状和轨道倾角, 是确定量, 而第六个参数真近点角 θ f 则描述目标航天器某一时刻在轨道上的位置, 随飞行时间变化 以下为书写方便, 省略下标 f, 目标航天器在终端时刻的位置可用轨道根数 4

59 第三章时间最优轨道机动 Ⅰ: 开普勒型方程 1 u t(h) 1 u t(h) 1 u t(h) 1 u t(h) 图 3.3 时间最优轨道拦截最优控制轨迹,1N 表示为 r f 1 = a(1 e2 ) (cos Ω cos(ω + θ) sin Ω sin(ω + θ) cos(i)) 1 + e cos θ r f 2 = a(1 e2 ) (cos Ω cos(ω + θ) + cos Ω sin(ω + θ) cos(i)) 1 + e sin θ r f 3 = a(1 e2 ) (sin(ω + θ) sin(i)) 1 + e cos θ (3.29) 其中的真近点角是关于时间变化的 µ θ = a 3 (1 e 2 ) (1 + e cos 3 θ)2 (3.3) 方程 3.3 难以写出 θ 关于时间 t 的解析表达式, 为便于数值处理, 将其引入系 统方程, 从而 θ 成为一新的系统变量 根据极大值原理横截条件可得以下终端约 41

60 3.3. 数值仿真 t fmin (h) /T max (1/N) 图 3.4 最小飞行时间和最大推力幅值关系 束条件 p r1 = p r2 = p r3 = 3 i=1 3 i=1 3 i=1 λ i ϕ i r 1 = λ 1 λ i ϕ i r 2 = λ 2 λ i ϕ i r 3 = λ 3 (3.31) 其中 λ i, i = 1, 2, 3 是 Lagrange 乘子, 而关于 θ 的横截条件则为 p θ = 3 i=1 λ i ϕ i θ = 3 i=1 λ i r f i θ (3.32) 由于终端时刻速度矢量和质量均是自由的, 因而由横截条件可知 p v = p m = (3.33) 42

61 第三章时间最优轨道机动 Ⅰ: 开普勒型方程 图 3.5 时间最优轨道拦截 3D 轨迹,1N Lagrange 乘子 λ i, i = 1, 2, 3 由方程组 3.31 确定, 将求解得到的 λ i 代人方程 3.32 可得到一个终端边值条件, 因此方程 3.28,3.33 和 3.32 组成了 7 个终端边值条件从而将问题转化为两点边值问题 数值仿真中, 初始数值条件为 (r 1, r 2, r 3 )() = (2.881Mm, Mm,.293Mm) (v 1, v 2, v 3 )() = ( 2.231Mm/h, 3.458Mm/h,.685Mm/h) m() = 15kg θ() = 1.472rad 而目标轨道参数为 (a f, e f, i f, Ω f, ω f ) = (32.378Mm,.5, 1.472rad,.5236rad,.5236rad) 以下是最大推力幅值 T max = 2N 的数值仿真结果, 图 3.8 是在惯性坐标系下 的最优飞行轨迹, 左上图中点划线表示的是目标轨道, 其上的实线部分则表示 任务完成时间内目标航天器的飞行轨迹, 由图可以看出, 拦截航天器的轨道倾 43

62 3.4. 小结 r 1 (Mm) r 2 (Mm) r 3 (Mm) t(h) R(Mm) t(h) t(h) t(h) 15 v 1 (Mm/h) v 2 (Mm/h) v 3 (Mm/h) t(h) t(h) t(h) V(Mm/h) t(h) m(kg) t(h) 图 3.6 时间最优轨道拦截状态变化轨迹,1N 角变化较小且比较缓慢, 在任务时间内大部分是实现轨道半长轴的增大, 由左下图也可明显看出 图 3.9 和图 3.1 则分别是状态和最优控制随时间变化轨迹, 与算例 1 的仿真结果作比较, 我们则可以看出位置和速度分量仍是振荡的, 但是在地球惯性坐标系下的绝对位置和速度则是近似线性变化的, 这与图 3.8 的结果是相对应的 从图 3.1 则可看出, 在任务总时间的前半段时间内, 控制分量在 u 1 方向几乎一直以最大值允许, 而在 u 2 方向的分量几乎为 ; 在任务即将完成时间段, 控制分量在 u 1, u 2 方向几乎为, 而在 u 3 方向则几乎一直以最大推力工作 这是因为在轨道坐标系下, 推力方向沿径向 - 法向标架分解,u 1 方向主要是控制轨道半长轴,u 2 主要控制轨道偏心率, 而 u 3 则主要控制轨道倾角的变化 这与图 3.8 的结果是相对应的 3.4 小结在本章, 考虑了用位置和速度作为变量描述的轨道动力学方程, 对于时间最优轨道拦截问题, 应用极大值原理, 分别研究了推力矢量在惯性坐标系和轨道坐标系下分解的最优控制律设计 在假设终端时刻卫星质量是自由的条件 44

63 第三章时间最优轨道机动 Ⅰ: 开普勒型方程 u t(h) 1 u t(h) 1 u t(h) 1 u t(h) 图 3.7 时间最优轨道拦截最优控制轨迹,1N 下, 依据极大值原理理论证明了不存在奇异控制情形, 最优控制几乎处处满足 u = 1 应用极大值原理一阶必要条件和横截条件, 将时间最优轨道拦截问题转化为两点边值问题 应用单值打靶法求解相应的边值问题, 数值仿真结果表明轨道拦截问题能够完成, 位置和速度分量随时间呈振荡变化 ; 改变推力幅值大小重复数值试验, 结果表明最大推力幅值和最小飞行时间的乘积近似成常数关系 45

64 3.4. 小结 图 3.8 时间最优轨道拦截 3D 轨迹,2N r 1 (Mm) r 2 (Mm) r 3 (Mm) m(kg) t(h) t(h) t(h) t(h) 15 R(Mm) 149 v 1 (Mm/h) v 2 (Mm/h) v 3 (Mm/h) V(Mm/h) t(h) t(h) t(h) t(h) t(h) 图 3.9 时间最优轨道拦截状态变化轨迹,2N 46

65 第三章时间最优轨道机动 Ⅰ: 开普勒型方程 1 u t(h) 1 u t(h) 1 u t(h) 1 u t(h) 图 3.1 时间最优轨道拦截最优控制轨迹,2N 47

66

67 第四章时间最优轨道机动 Ⅱ: 轨道根数描述方程 第四章 时间最优轨道机动 Ⅱ: 轨道根数描述方程 在第三章, 分别应用在地心坐标系和轨道坐标系下位置 - 速度表示的卫星轨道运动方程, 研究了时间最优轨道拦截问题, 根据极大值原理设计了最优控制律, 并证明了应用开普勒型方程不存在奇异最优解和奇异控制 根据数值仿真结果可以看出该方程解具有振荡性, 且假设的路径约束条件 A1 不易满足, 数值上处理比较困难 因此, 在本章应用更具普遍性的改进轨道根数描述轨道动力学方程来研究时间最优轨道机动问题, 并考虑 J2 摄动的影响 应用改进轨道根数描述的轨道动力学方程, 根据几何方法, 系统满足可控性条件 [61] 根据 Fillopov 定理可以证明最优解的存在性 [64], 但是最优控制是否是唯一和连续的仍然是开问题 [113] 文献 [3] 研究了单输入系统的可控性, 在切向 - 法向标架下, 沿 f n, f c 方向在椭圆域内不是无条件可控的, 而控制沿切向的单输入系统 ẋ = f + uf t 在椭圆域内是完全可控的 若记 p 为协状态,H = p, f, H t = p, f t, 则在 奇异面 Σ = {(p, x) : H = {H, H t } = } 上有控制输入 u s = {H, {H, H t }} {H t, {H, H t }} 式中 :{, } 为 Poisson 括号 而对于多输入仿射非线性系统 ẋ = f + m i=1 u if i 奇异控制则满足 {H i, H } + m j=1 {H i, H j }u j = 其中 [{H i, H j }] 1 i,j m 是对角元素为 的反对称矩阵, 求解控制输入 u i 需要分 m 是奇数和偶数两种情况 [64, 114] 众多学者应用改进轨道根数描述的动力学方程开展了时间最优轨道转移问题研究 [3, 31, 33, 61 63, 17, 18] 文献 [62] 研究了共面轨道转移情形, 证明了 49

68 4.1. 问题描述 切换函数是连续可微的而且仅存在有限多个孤立零点, 即不存在奇异弧而仅存在切换点, 而在任意两个切换点之间, 解是光滑连续的 进而证明了在奇异点都是 π 奇异, 即在切换点控制输入方向立刻旋转 18 度 文献 [115, 116] 应用单值打靶法研究了椭圆轨道向地球同步轨道转移问题, 数值计算结果表明最大推力幅值与最优轨道转移时间有如下关系 t f min T max c 这是数值计算的结果, 理论上仍是开问题 [113] 在第三章的讨论中, 可知轨道机动任务也包括轨道拦截和交会, 它们和轨道转移任务之间的区别主要在终端边值条件 根据作者目前掌握的文献来看, 应用改进轨道根数描述的高斯动力学方程研究时间最优轨道拦截和交会的文献很少见到 因此, 在本章, 在统一的框架下, 我们应用改进轨道根数描述的高斯动力学方程研究包括轨道转移 交会和拦截在内的时间最优轨道机动问题, 应用极大值原理, 将轨迹优化问题转化为两点边值问题, 同时分别讨论了轨道转移 交会和拦截问题的终端边值条件 本章的主要贡献在于 :1) 应用改进轨道根数描述的高斯型动力学方程, 间接优化方法用于研究时间最优轨道拦截和交会问题, 并详细讨论了各自终端条件的描述方法 2) 数值仿真结果表明, 最大推力幅值和最小飞行时间近似成常数关系的猜想同样适用于轨道拦截和交会问题 ; 同时也表明, 最大推力幅值固定时, 轨道交会和转移问题几乎具有相同的最优轨迹 4.1 问题描述考虑应用改进轨道根数描述卫星的轨道运动, 令 x = [P e x e y h x h y L] T, 则描述卫星运动的系统方程为 ẋ = f + T max m 3 u i f i + i=1 3 g i f i (4.1) i=1 其中 T max 为常值推力发动机最大推力幅值,u = [u 1 u 2 u 3 ] 为发动机推力 方向矢量 若推力在轨道坐标系下沿径向 - 法向 (R,T,N) 标架下分解, 则定义卫星 5

69 第四章时间最优轨道机动 Ⅱ: 轨道根数描述方程 动态的四个向量场如方程 (2.12) 描述, 写成向量形式为 µ [ ] T f = W 2 P P P [ f 1 = sin L cos L µ P 1 [ f 2 = µ W f 3 = P µ 1 W [ 2P W cos L + η x W sin L + η y ] T Ze y Ze x C 2 cos L C 2 sin L Z ] T ] T (4.2) 其中 µ 为地球引力常数, 其他参数则有如下形式 : W = 1 + e x cos L + e y sin L, η x = e x + cos L, η y = e y + sin L, Z = h x sin L h y cos L, C = 1 + h 2 x + h 2 y 方程 4.1 中的 g i 表示摄动作用力, 若是考虑地球扁率摄动的影响, 一阶扁率在径向 - 法向方向分解用经典轨道根数描述为 [2] g 1 = 3 2 J 2 用改进轨道根数描述则为 [117] µre 2 r 4 [1 3 sin2 i sin 2 (ω + f)] g 2 = 3 2 J µre 2 2 sin 2 i sin 2(ω + f) r 4 g 3 = 3 2 J µre 2 2 sin 2i sin 2 (ω + f) r 4 g 1 = 3 2 J µre 2 2 r [1 12(h x sin L h y cos L) 2 ] 4 (1 + h 2 x + h 2 y) 2 µre 2 (h x sin L h y cos L)(h x cos L + h y sin L) g 2 = 12J 2 r 4 (1 + h 2 x + h 2 y) 2 µr 2 (h e x sin L h y cos L)(1 h 2 x h 2 y) g 3 = 6J 2 r 4 (1 + h 2 x + h 2 y) 2 51

70 4.1. 问题描述 其中 r = P 1 + e x cos L + e y sin L 令 f = f + 3 g i f i i=1 则系统可描述为 ẋ = f + T max m 3 i=1 u if i (4.3) 考虑卫星质量变化方程为 ṁ = βt max u (4.4) 而控制输入需满足约束 u 1 对于轨道机动任务, 仍需满足边值约束 Φ(, x(), m()) = Φ(t f, x(t f ), m(t f )) = 至此, 我们完成了轨道机动问题的简单描述, 为便于应用极大值原理设计时间最优控制律, 控制约束改写成 u = 1, T T max 对于固定比冲常值推力发动机, 控制则变成推力方向矢量和推力幅值大小 u c = T u, T T max, u = 1 (4.5) 52

71 第四章时间最优轨道机动 Ⅱ: 轨道根数描述方程 综合考虑, 时间最优轨道机动问题可描述为 min J = t f dt 性能指标 ẋ = f (x) + T 3 m i=1 u if i (x) 卫星动态 ṁ = βt 质量变化 (T P ) t Φ(, x(), m()) = 初始条件 (4.6) Φ(t f, x(t f ), m(t f )) = u = 1, T T max 终端条件 控制约束 为简化研究问题, 我们考虑以下假设条件 A1) 系统模型的状态始终满足路径约束 A = {(x, m) P >, (e x, e y ) < 1, m m s } 即卫星在椭圆域内飞行, 在地心坐标系下位置向量幅值始终大于地球半径,m s 为不携带燃料时的卫星本体质量 A2) 初始轨道和目标轨道是相异的, 即最优控制是非空的 A3) 卫星在终端时质量是自由的 4.2 正则最优控制律设计 这一节, 我们研究时间最优轨道机动 (T P ) t 的最优控制问题, 应用极大值 原理, 相应的 Hamilontian 函数为 H = p + H + T m 3 i=1 u ih i βp m T + ν( u 2 1) (4.7) 在式 4.7 中, p 是一个非正的常数,H i (i =, 1, 2, 3) 是 Hamilontian 提升 p, f i,p 是状态 x 对应的协变量,p m 是质量 m 对应的协变量,ν 则是 Lagrange 乘子 则正则情况下,p 是负常数且被整定取为 1 定义 Ψ = 53

72 4.2. 正则最优控制律设计 (H 1, H 2, H 3 ), 而最优控制则由下式确定 H = u u = 1 2 H u 2 (4.8) 当任何时候 Ψ 时, 求解方程组 可得到 ν 和 u 定义切换函数为 ν = T 2m Ψ u = Ψ Ψ S = Ψ m βp m (4.9) 因此, 可将与推力有关的项合并起来, 从而 Hamiltonian 函数可以重新写成 H = p + H + S T (4.1) 根据极大值原理一阶必要条件, 描述协变量 p m 的动态方程为 ṗ m = T m 2 Ψ 而根据假设 A2, 由极大值原理横截条件可知 p m (t f ) =, 由上式我们可以看 出 p m 是非正的增函数且趋近于 p m (t f ), 因此切换函数 S 是非负的, 当推力幅值 最大时 Hamiltonian 函数也取最大 根据以上分析, 有以下命题 : 命题 4.1: 当假设 A1 和 A2 满足, 若 Ψ = (H 1, H 2, H 3 ) 不为 时, 则时间最优轨道机 动问题 (T P ) t 的最优控制为 u = Ψ Ψ, T = T max (4.11) 54 上面我们考虑了正则最优控制律设计, 接下来考虑奇异情况

73 第四章时间最优轨道机动 Ⅱ: 轨道根数描述方程 令 (x, m, p, p m, u) 是系统极值解, 极值的分类则依赖于轨迹与切换面 Ψ = 的接触阶数 当 Ψ, 控制 u 是光滑的且有式 4.11 形式时, 则说极值解是 阶的, 若 Ψ 时则是奇异的 [3] Ψ 是一个在时域内与奇异弧和切换点相关的函数, 文献 [62] 已经证明了 Ψ 是连续可微的, 几何分析也表明在集合 {Ψ = } 内仅存在有限多个切换点 因此有下述命题 : 命题 4.2: 对于固定比冲发动机, 设发动机推力最大幅值为 T max, 若在路径约束安全区域 A 内存在最优轨迹 (x, m), 则时间最优轨道机动问题 (T P ) t 相应的最优控制几乎处处满足 u = 1, T = T max 注记 4.1: 在文献 [3] 和 [62] 中, 作者研究了时间最优轨道转移问题, 得到了类似命题 4.1 和 4.2 的结论 应用 Pontryagin 极大值原理, 轨道拦截或交会问题与轨道转移问题的主要区别在于终端条件, 因此命题 4.1 和 4.2 同样适用于轨道拦截和交会问题 式函数 根据命题 4.2, T = T max 几乎处处满足, 因此卫星质量可以表示成时间的显 m(t) = m βt max t (4.12) 其中 m 是卫星在任务初始时刻的质量 进一步, 为简化问题, 我们令 τ = t t f, 从而将飞行时间整定到区间 [, 1], 最终飞行时间 t f 则作为系统一个新的常数变 量 因此, 时间最优轨道机动问题 (T P ) t 可转化为如下 Mayer 型最优性能指标的 形式 (T P ) τ min J = t f (1) ( dx dτ = t ) f f (x) + Tmax 3 m(t f τ) i=1 u if i (x), τ [, 1] dt f dτ = Φ(, x(), t f ()) =, Φ(1, x(1), t f (1)) = u = 1 (4.13) 若 (x, t f, u) 是问题 (T P ) τ 的解, 则有绝对连续协变量 p 和 p tf 分别对应于状 态 x 和 t f, 因此根据极大值原理一阶必要条件,(x, t f, p, p tf ) 是下述两点边值问 55

74 4.3. 边值问题和横截条件 题的解 ẋ = H p t f = p tf ṗ = H x = H t f (4.14a) (4.14b) (4.14c) (4.14d) 同时满足初始条件 Φ() = (P () P, e x () e x, e y () e y, h x () h x, h y () h y, L() L, p tf () p t f ) = (4.15) 和终端边值条件 Φ(1) = Φ(1, x(1), t f (1), p(1), p tf (1)) = (4.16) 对于轨道转移 拦截和交会问题, 假设满足同样的初始边值条件, 而终端边值 条件则需要由具体的任务来确定 式 4.14 中, 新的 Hamilontian 函数具有如下形 式 ( H = p, t f f (x) + T ) max m(t f τ) (u 1f 1 (x) + u 2f 2 (x) + u 3f 3 (x)) (4.17) 其中最优控制 u = (u 1, u 2, u 3) 仍然是由 4.11 定义的函数 至此, 将连续推力时 间最优轨道机动问题转化为两点边值问题, 下一节我们将讨论不同轨道机动使 命的终端边值条件 4.3 边值问题和横截条件 应用极大值原理, 我们将时间最优轨道机动问题转化为两点边值问题 作 为问题出发点, 我们首先讨论时间最优轨道转移问题, 在此基础上, 在统一的 框架下研究时间最优轨道转移 拦截与交会问题 首先给出轨道机动方式的简 单描述 : 56

75 第四章时间最优轨道机动 Ⅱ: 轨道根数描述方程 轨道转移 : 从初始轨道运动到目标轨道, 初始轨道和目标轨道已知, 但航天器到达目标轨道的位置不确定 轨道拦截 : 初始轨道和航天器在初始轨道上的位置已知, 目标轨道已知, 要求在终端时刻转移轨道上的航天器与目标轨道上的拦截目标的相对位置为, 相对速度自由 轨道交会 : 初始轨道和航天器在初始轨道上的位置已知, 目标轨道已知, 要求在终端时刻转移轨道上的航天器与目标轨道上的交会目标的相对位置为, 相对速度也为 根据上面对轨道机动方式的定义, 可知在初始条件相同的情况下, 轨道转移 拦截与交会的差别主要是终端条件 令 (P f, e f x, e f y, h f x, h f y, L f ) 为描述目标轨道根数, 接下来, 我们分别讨论三种机动方式的终端条件 轨道转移 轨道转移问题得到了广泛研究, 其终端目标轨道满足条件 [31, 17, 18] (P (1) P f, e x (1) e f x, e y (1) e f y, h x (1) h f x, h y (1) h f y) = 和 (p L (1) p f L, p t f (1) p f t f ) = 卫星在目标轨道上的位置不确定, 因而终端时 刻变量 L 是自由的, 根据极大值原理横截条件可知 p f L = 而系统的性能指标为 Mayer 型, 最小飞行时间为 t f, 同样根据横截条件可知 p f t f = 1 因此, 时间最优轨道转移问题满足终端边界条件 Φ(1) = (P (1) P f, e x (1) e f x, e y (1) e f y, h x (1) h f x, h y (1) h f y, p L (1), p tf (1) + 1) = (4.18) 注记 4.2: 事实上, 若终端时刻累积赤经 L 也是固定的, 则航天器在目标轨道的位置 速度和绕地球飞行圈数也是固定的 若同时飞行时间 t f 也预先确定, 这 57

76 4.3. 边值问题和横截条件 种情况则可以看做是轨道交会问题 [6] 但是在时间最优轨道机动问题中, 赤经 L 与时间 t 一样也是快变量, 因此最终赤经也与最大推力幅值相关, 并不能随意选取 已有文献研究表明, 最小累积赤经和最大推力幅值的乘积近似成常数关系 [32, 112] 因此, 对于固定终端累积赤经的时间最优轨道转移问题, 并不能简单的看作轨道交会问题 轨道拦截 对于轨道拦截问题, 要求拦截星与拦截目标在终端时刻相对位置为, 而相 对速度自由 假设目标星在固定轨道做无控运动, 初始位置给定 因此, 拦截 星的终端约束不是一个固定轨道而是一个子流形, 满足如下条件 : ϕ 1. = r1 (1) r f 1 = (4.19a) ϕ 2. = r2 (1) r f 2 = (4.19b) ϕ 3. = r3 (1) r f 3 = (4.19c) 其中 (r 1, r 2, r 3 ) 表示拦截星在惯性坐标系下的位置, 可用改进的轨道根数来描述 ;(r f 1, r f 2, r f 3 ) 是目标星在地球惯性坐标系下的位置, 可将 (P f, e f x, e f y, h f x, h f y, L f ) 代人公式 (2.4a-2.4c) 得到 然而, 给定初始位置, 若飞行时间未知, 则目标星在自身轨道上的终端位置也是未知的, 位置主要由赤经决定, 因此将描述目标星在目标轨道上位置的变量 L f 的无控方程 dl f dτ = t f µ P f (1 + e f x cos L f + e f y sin L f ) 2 P f (4.2) 引入系统作为一个方程, 则系统状态变量则为 (x, t f, L f ) R 8 R + R +, 而将 L f () L f = 引入方程 4.15 作为一个初始边界条件 根据极大值原理的横截条件, 由终端约束 4.19, 可以得到关于 x 的协状态 p = [p P, p ex, p ey, p hx, p hy, p L ] T 的终端边值条件 : ϕ 1 p P (1) = (λ 1 P (1) + λ ϕ 2 2 P (1) + λ ϕ 3 3 P (1) ) ϕ 1 p ex (1) = (λ 1 e x (1) + λ ϕ 2 2 e x (1) + λ ϕ 3 3 e x (1) ) (4.21a) (4.21b) 58

77 第四章时间最优轨道机动 Ⅱ: 轨道根数描述方程 ϕ 1 p ey (1) = (λ 1 e y (1) + λ ϕ 2 2 e y (1) + λ ϕ 3 3 e y (1) ) ϕ 1 p hx (1) = (λ 1 h x (1) + λ ϕ 2 2 h x (1) + λ ϕ 3 3 h x (1) ) ϕ 1 p hy (1) = (λ 1 h y (1) + λ ϕ 2 2 h y (1) + λ ϕ 3 3 h y (1) ) ϕ 1 p L (1) = (λ 1 L(1) + λ ϕ 2 2 L(1) + λ ϕ 3 3 L(1) ) r f 2 r f 3 r f 1 p L f (1) = λ 1 L f (1) + λ 2 L f (1) + λ 3 L f (1) (4.21c) (4.21d) (4.21e) (4.21f) (4.21g) 其中 λ i (i = 1, 2, 3) 是 Lagrange 乘子. 从非线性代数方程组 4.21 中任选 3 个方程求解出 λ i (i = 1, 2, 3), 再代入方程组 4.21 另外剩余 4 个方程则可得 4 个边界条件 因此这 4 个边界条件, 方程 4.19 再加上 p f t f = 1 组成了时间最优拦截问题的 8 个终端边值条件 注记 4.3: 终端边值 4.21 也可根据坐标变换得到 令 x M, y = f(x) N, 考虑坐标变换,y 表示成 x 根据极大值原理可知协变量分别为 p x T x M, p y T f(x)n 定义运算 f x : T xm T f(x) N f 若 q x T x M 则可导出 q x x T f(x) M 而 p x, q x =, 因而有 p y, f x q x = ( 则有 f x) T f(x) N Tx M 由此可以导出 ( ) f p y, q x = x ( ) f p y = p x x 即 p y = (( ) ) 1 f p x (4.22) x 59

78 4.3. 边值问题和横截条件 考虑位置和速度 (r, v) 表示的状态变量与 MEE 之间的转化, 对于拦截问题则有速度对应的协状态终端约束为 p v (t f ) =, 因此可根据公式 4.22 式求得 3 个对应的轨道根数表示的协状态终端约束条件, 避免了求取 Lagrange 乘子 λ i (i = 1, 2, 3), 再加上 3 个终端位置约束, 则构成了 6 个终端条件 但是为将上述轨道拦截化成两点边值问题, 仍需另外寻找其他合适的终端边值条件 轨道交会 与轨道拦截问题不同, 轨道交会问题则要求交会星与交会目标在终端时刻 相对位置和相对速度均为, 因此终端约束则需满足如下条件 :. ϕ 1 = r1 (1) r f 1 = (4.23a). ϕ 2 = r2 (1) r f 2 = (4.23b). ϕ 3 = r3 (1) r f 3 = (4.23c). ϕ 4 = v1 (1) v f 1 = (4.23d). ϕ 5 = v2 (1) v f 2 = (4.23e). ϕ 6 = v3 (1) v f 3 = (4.23f) 同样, 式中 (r f 1, r f 2, r f 3 ) 和 (v f 1, v f 2, v f 3 ) 表示目标星在惯性坐标系下的终端位置和速度分量, 可将 (P f, e f x, e f y, h f x, h f y, L f ) 代人公式 ( 2.4a-2.4f) 得到 同样, 将描述 L f 的无控方程 4.2 引入系统作为系统一个方程, 则系统状态变量则为 (x, t f, L f ) R 8 R + R +, 而将 L f () L f = 引入方程 4.15 作为一个初始边界条件 根据极大值原理的横截条件, 由终端约束 4.23, 可以得到如下终端边值条件 : p(1) = 6 i=1 λ i ϕ i x(1) (4.24) 于 p L f 求解上非线性代数方程组求得拉格朗日乘子 λ i, i = 1,, 6, 然后代入关 的协状态终端条件 r f 2 r f 3 r f 1 p L f (1) = λ 1 L f (1) + λ 2 L f (1) + λ 3 L f (1) + λ 4 L f (1) + λ 5 L f (1) + λ 6 L f (1) (4.25) 6 v f 1 v f 2 v f 3

79 第四章时间最优轨道机动 Ⅱ: 轨道根数描述方程 表 4.1 时间最优轨道机动问题的数值条件 变量 初始条件 终端条件 P 15.6Mm 32.1Mm e x e y..3 h x.612. h y..1 L π 自由 L f 1.57 自由 可以得到一个终端边值条件 因此方程 4.23, 4.25 和 p tf (1) = 1. 组成了时间最优交会问题的 8 个终端边值条件 4.4 数值仿真在 4.2 节, 应用极大值原理设计了时间最优控制律, 而 4.3 节则讨论了不同轨道机动使命的边值条件, 从而将时间最优轨道机动问题转化为两点边值问题 在这一节, 我们则应用打靶法求解相应的两点边值问题 以轨道转移问题为例, 两点边值问题等价于寻找所谓的打靶函数的零点, 即 : 寻找 (p, t f ) R6 R 使得下式成立 b(φ (x, t f, p, p t f )) = (4.26) 其中的边界函数 b 由式 (4.18) 确定 在数值仿真中, 系统的物理常数 µ, β 与第三章数值仿真算例相同, 初始质量为 m = 15kg, 后续章节仍采用该数值, 不再另行说明 数值边值条件则如表 4.1 所示, 表中的单位 1Mm = 1 6 m 首先在同一幅图中给出时间最优轨道转移 拦截与交会问题最优状态变化轨迹, 如图 4.1 所示 图示给出了六个轨道根数 P, e x, e y, h x, h y, L 和卫星质量 m 随时间变化轨迹, 横坐标表示时间历程, 实线表示轨道转移, -. 线表示轨道交会, 而 - - 则表示轨道拦截, 图中水平直线则表示目标的轨道根数值 对于同样的最大推力 4N, 轨道转移任务的最小飞行时间是 h, 交会任务的最小飞行是 h, 而轨道拦截任务的最小飞行时间是 h 61

80 4.4. 数值仿真 e x P(Mm) time(h) time(h).5 transfer L h x rendezvous intercept time(h) time(h).2 e y h y time(h) time(h) m(kg) time(h) 图 4.1 状态最优轨迹,4N 由图 4.1 可以看出, 状态变量随时间的变化轨迹比较平滑, 也可以看出最优交会问题和转移问题几乎具有相同的最优轨迹, 状态以几乎相同的轨迹到达目标轨道期望值 然而, 轨道拦截问题的最终轨道参数值与目标值相差较大, 由图 4.2 地球坐标系下三维空间轨迹图也可以看出, 拦截星的最后运行轨道几乎垂直于目标星轨道 图 4.2 所示的是时间最优轨道拦截在地球惯性坐标系下 3D 空间的最优变化轨迹, 箭头代表推力的实时方向 实线表示初始轨道, 点划线则表示目标星的固定运行轨道 ; * 表示拦截星的初始点, o 则表示拦截星与目标星的遭遇点 同样, 图 4.3 所示的是时间最优轨道交会问题的 3D 最优轨迹, 由图 8 图 9 可以看出, 轨道拦截情形时的拦截星的运行轨道倾角与轨道交会情形相比则变化较小 由于轨道转移和交会几乎有相同的最优轨迹, 轨道转移问题的 3D 最优轨迹不再给出图示 调整最大推力幅值大小, 重复上述数值实验, 三种轨道机动方式的最小飞行时间总结于表 4.2, 最小飞行时间和最大推力之间的关系则由图 4.4 所示 由图 4.4 左图可以看出, 对于固定最大推力幅值, 轨道转移和轨道交会时的最小 62

81 第四章时间最优轨道机动 Ⅱ: 轨道根数描述方程 图 4.2 轨道拦截 3D 最优轨迹,4N 图 4.3 轨道交会的 3D 最优轨迹,4N 63

82 4.4. 数值仿真 表 4.2 最小飞行时间 T max (N) 轨道转移轨道交会轨道拦截 t fmin (h) b t fmin (h) b t fmin (h) b E E E E E E E E E E E E E E E E E E E E E E E E E E E E-13 飞行时间几乎相同, 而轨道拦截情形的最小飞行时间相比则小 一些学者研究 时间最优轨道转移问题时, 发现最大推力幅值和最小时间的乘积近似成常数关 系 :t fmin T max c [31], 由图 4.4 也可以看出, 该结论同样适用于轨道拦截和 交会问题 然而, 从左图的局部放大图和右图, 我们也可以看出, 这个关系不 适用于大推力情况 随着推力幅值的增大, 到某一值 ( 如 6N) 时, 最小飞行 时间并不随最大推力增大而显著减小 因而, 自然地可以猜测, 脉冲推力可以 看作是某种情况下连续推力的极限形式, 但是需要理论进一步的证明 另外需 要说明的是, 对于 15kg 的卫星产生 6N 的推力并不现实, 在这里只是为了数 值仿真的需要 对于轨道拦截问题, 我们改变拦截星在初始轨道上的出发位置, 即分别选 择 L =, π/2, π, 3π/2, 重复上述数值实验, 相应的最小飞行时间和最大推力 的关系如图 4.5 所示 由图可以看出, 对于相同的推力幅值, 越早实施轨道拦 截, 所需时间越少 ; 最小时间和最大推力的乘积值 c 与 L 的选取近似成线性增 长方式, 对于本算例则可简单归为 c = 6 + 8L π 64

83 第四章时间最优轨道机动 Ⅱ: 轨道根数描述方程 25 transfer 25 transfer rendezvous rendezvous intercept intercept t fmin (h) t fmin (h) /T max (1/N) T max (N) 图 4.4 t fmin T max 近似常数关系 4.5 小结采用改进轨道根数描述卫星轨道运动, 应用间接法研究了常值连续推力时间最优轨道转移 拦截与交会问题 应用极大值原理设计最优控制律, 理论分析表明, 卫星一直以最大推力方式运行, 根据极大值原理一阶必要条件将轨迹优化问题转化为两点边值问题 轨道转移 拦截与交会方式的主要差别在于终端约束条件, 根据横截条件, 分别研究了三种机动方式的终端边值条件 应用单值打靶法分别研究了相应的两点边值问题, 数值仿真结果表明 : (1) 轨道转移和轨道交会使命对于同一最大推力几乎有相同的最优轨迹 ; (2) 对于相同的最大推力幅值, 轨道拦截所需的时间要比轨道转移所需时间短 ; (3) 最小时间与最大推力幅值乘积近似成常数关系, 这一数值结论同样适用于轨道拦截与交会问题 ; (4) 对于轨道拦截问题, 越早实施变轨机动, 所需任务时间越少 但是上述结论只是数值仿真得到结果, 并未得到理论严格的证明, 需要进一步 65

84 4.5. 小结 t fmin (h) 5 L = L =π/2 L =π L =1.5π /T max(1/n) 图 4.5 初始时刻 L 对轨道拦截时间的影响 的开展工作, 完善上述结论 66

85 第五章近似最优赤经轨道转移 第五章 近似最优赤经轨道转移 在第四章, 我们应用间接法研究了时间最优轨道机动问题, 数值结果如图 4.1 所示表明变量 L 随时间成线性增长 ; 而对于轨道拦截问题, 改变拦截星在初始轨道上的出发位置, 即选择不同的初始赤经, 数值试验也表明最小时间和最大推力的乘积值 c 与 L 的选取近似成线性增长方式 这一章, 我们应用间接法研究累积赤经最优轨道转移问题 对于连续推力最小燃料消耗轨道转移问题, 由于赤经 L 和时间 t 一样都是快变量, 因此该问题根据最终时间 t f 和最终赤经 L f 划分为四种情形 :1) t f 和 L f 均自由 ;2) t f 自由,L f 固定 ;3) t f 固定,L f 自由 ;4) t f 和 L f 均固定 事实上,L f 固定意味着卫星在目标轨道的位置和绕地球飞行圈数已确定, 该问题则对应着轨道交会问题 J.Gergaud 应用同伦算法研究了 t f 固定 L f 自由以及 t f 自由 L f 固定燃料最优轨道转移问题 [32] 在 t f 和 L f 均自由情形, 虽没有理论证明, 但是事实表明该问题无解 [32, 112] 在最小燃料消耗问题上, 最大推力幅值给定时, 固定飞行时间 t f 可通过 t fmin (h) T max (N) c 来确定 [115,116] 同样为保证任务可行性, 最终 L f 也不能随意选取, 也必须严格大于最小转移赤经 因此研究最小累积赤经轨道转移问题非常有意义, 该问题在文献 [32,112,118] 简单的提到, 并给出了最小赤经和最大推力的乘积近似成常数的结论, 但是文献作者并没有给出解决该问题的实际方法 在文献 [118, 119] 中 平均近似方法被用来研究连续推力能量最优轨道转移问题, 用累积赤经 L 来重新参数化系统, 仿射控制系统则转换成无漂移项的次黎曼系统, 但是作者并没有考虑卫星质量的变化 而且, 文献 [118] 仅考虑了推力约束下的共面轨道转移, 文献 [119] 虽考虑了非共面转移, 但却没有考虑输入约束 本章应用间接最优方法研究了连续推力最小累积赤经轨道转移问题, 其中 L 代替时间作为系统的自变量, 对于平均化近似系统, 应用极大值原理设计最优控制律, 并将问题转化为两点边值问题 主要贡献有两点, 一是系统研究了连续推力最小赤经轨道转移问题, 包括问题的描述和控制律设计 ; 二是用 L 将系统参数化, 采用近似平均化方法研究了变质量和带有推力约束的非共面轨道转移问题 67

86 5.1. 问题描述 5.1 问题描述 采用改进轨道根数描述轨道动力学方程, 定义状态变量 x = [P e x e y h x h y ] T, 卫星发动机产生的推力在轨道坐标系 (R,T,N) 标架下 分解, 则卫星轨道转移可由下式描述 : 3 ẋ = 1 m u if i (L, x) i=1 L = g (L, x) + 1 (5.1) m g 1(L, x)u 3 其中 和 f 1 = f 2 = f 3 = P µ P µ P µ 1 (W sin L W cos L ) W e x e y 1 W (2P P + (W cos L + η x) e x + (W sin L + η y ) e y ) 1 W ( Ze y + Ze x + C e x e y 2 cos L + C h x 2 sin L ) h y µ g = P P g 1 = µ W 2 P Z W L L (5.2) 其中的参数为 W = 1 + e x cos L + e y sin L, η x = e x + cos L, η y = e y + sin L, 考虑卫星质量变化, 与发动机推力幅值成比例 Z = h x sin L h y cos L, C = 1 + h 2 x + h 2 y ṁ = β u (5.3) 因此描述卫星运动的状态实际为 (x, L, m) M R + R + 考虑任务, 卫星的 68

87 第五章近似最优赤经轨道转移 飞行轨迹必须满足安全区域 [62] A = {(x, m) P >, (e x, e y ) < 1, m m s } (5.4) 即卫星在椭圆域内飞行, 在地心坐标系下位置向量幅值始终大于地球半径,m s 为无燃料的卫星质量 对于轨道转移问题, 描述卫星初始轨道和目标轨道的边界约束为 x() = x, L() = L, m() = m, Φ(x(t f )) = (5.5) 而控制输入则需满足 u T max (5.6) 其中 T max 为发动机最大推力幅值 接下来, 设 u = T max ū 和 ϵ = T max /m 由于变量 L 和时间 t 一样是快变 量, 终端时间 t f 不固定, 因此需重新参数化系统为 : dx dl = ϵ g + ϵg 1 ū 3 3 i=1 ūif i (5.7) 卫星质量变化为 : dm dl = βmϵ ū (5.8) g + ϵg 1 ū 3 因此最小累积赤经轨道转移任务的性能指标为 : J = Lf L dl (5.9) 轨道转移的任务就是寻找绝对连续的状态 (x, m) 和在有界控制输入 ū 使得指标 5.9 最小化, 同时需要满足动态 (5.7),(5.8), 路径约束 (5.4), 边值条件 (5.5) 和控制约束 ū 近似最优控制律设计这一节, 我们研究最小累积赤经轨道转移的最优控制问题 应用极大值原 69

88 5.2. 近似最优控制律设计 理, 最优轨迹就是极值轨迹, 相应的 Hamiltonian(H) 函数为 : ( 3 ) ϵ H(L, x, m, p, p, p m, ū) = p + ū i H i βmp m ū g + ϵg 1 ū 3 i=1 (5.1) 在 (5.1) 中, p, H i = p, f i 是 Hamiltonian 提升, p 是对应 x 的协状态, 而 p m 是 m 对应的协状态. 正常情况下,p 是负的且整定为 1 为简化问题, 我们假 设终端质量是自由的, 因此约束 m m s 不起作用 令 σ = ϵg 1 ū 3 /g, 通过简单的数值计算和量纲匹配, 我们可以确定 σ < 1, 因此我们把 H 函数关于 σ 用泰勒级数展开有 : ( 3 ) H(L, x, m, p, p, p m, ū) = 1 + ϵ 1 ū i H i βmp m ū g 1 + ϵg 1 ū 3 /g i=1 ( = 1 + ϵ ( 3 ) 1 σ + σ2 ) ū i H i βmp m ū g i=1 ( 3 ) = 1 + ϵ ū i H i βmp m ū + O(σ) g i=1 (5.11) 忽略 O(σ) 项, 则 H 函数可描述为 H(L, x, m, p, p m, ū) 1 + T max mg ( 3 ) ū i H i βmp m ū i=1 (5.12) 命题 5.1: 在安全区域内, 沿最优解有, (Ⅰ) 质量对应的协变量 p m 是非正的增函数 ; (Ⅱ) 令 Ψ = (H 1, H 2, H 3 ), 当 Ψ 不为零时, 最优控制输入为 ū = Ψ Ψ (5.13) 证明 : 应用柯西 - 苏瓦兹不等式, 式 5.12 可写为 7 H(L, x, m, p, p m, ū) 1 + T max mg ( Ψ ū βmp m ū ) = 1 + T max mg ( Ψ βmp m ) ū

89 第五章近似最优赤经轨道转移 当 ū = αψ, α 时等号成立, 因此使得 H 函数取极大的解可写成 ū = α(y(l)) Ψ, α(y(l)) [, 1] Ψ 其中 y(l) 是系统的最优解 因为在路径约束内,g 大于, 因此定义切换函数 S = Ψ βmp m 则关于质量的协变量方程为 p m = T max m 2 g 3 ū i H i i=1 根据上式可知 p m 是增函数且趋近于 p m (L f ), 由于假设终端质量自由, 由 横截条件可知 p m (L f ) =, 因此 p m 是非正的 因为 p m 是非正的, 因此当 Ψ 不 为 时, 切换函数 S 严格大于, 因而 α(y(l)), 论断 (Ⅱ) 成立 令 (x, p, ū) 是最优解, 正则极值的分类则根据切换面 {Ψ = } 的接触阶数, 当 Ψ 不为, 由 (5.13) 式确定的控制输入 ū 是光滑的, 极值则称为 阶的 ; 若 Ψ = 则是奇异的 [3] 接下来, 我们讨论 Ψ 的性质 命题 5.2: 沿最优轨迹,Ψ 是连续可微的且仅有有限个孤立零点 证明 : 令 C 表示 Ψ 的零点集合 C = {L [L, L f ] Ψ = } 令 k(x) = T max mg, 因此相应的 H 函数为 : ( 3 ) H = k(x) ū i H i βmp m ū i=1 71

90 5.2. 近似最优控制律设计 此时的切换面函数 Ψ = (H 1, H 2, H 3 ) 是绝对连续的, 而且 Ψ = ( Ḣ 1, Ḣ2, Ḣ3), 而 { Ḣ i = {H i, H} = H i, k ( 3 )} ū i H i βmp m ū i=1 ( 3 ) = ū i H i βmp m ū {H i, k(x)} + k(x) j i ūj {H i, H j } i=1 令 3 g ij = ū i H i {H i, k(x)} + k(x) j i ūj {H i, H j }, L [L, L f ]\ C i=1 g ij =, L C 在切换点之外面, g ij 是连续的. 若令 L C 是切换点, 则有当 L L 时 g ij Ψ {H i, k(x)} + j i k(x) {H i, H j } 又因为 则 [f i, f j ] V ect({f 1, f 2, f 3 }), 1 i < j 3 {H j, H i }( L) = H [fj,f i ]( L) = p( L), [f j, f i ](x( L)) = 因此 Ψ 是连续的,Ψ 属于 C 1 接下来证明 Ψ 仅有有限个孤立零点 因为 L 是切换点, 则 Ψ( L) =, 那么有 H i ( L) =, 而 Ḣ i ( L) = βmp m ū f T i k x 根据命题 5.1, 可知 p m 是非正增函数, 因此 H i ( L) = 不可能同时满足. 此, Ψ( L) 而且每个切换点在区间 [L, L f ] 是孤立的, 也就是说 Ψ 仅有有 限个零点, 而不存在零值区间 因 根据上述命题的结论, 可知最优轨迹是零阶极值, 因而有下述命题 命题 5.3: 对固定推力 T max, 若最优解 (x, m) 在路径约束的内部, 任何最优控 制 ū 几乎处处满足 ū = 1 72

91 第五章近似最优赤经轨道转移 5.3 数值实验 为简化问题, 我们将问题整定在区间 [, 1], 取 L =, 令 τ = L L f, 则自由的 终端值 L f 作为一个新的系统变量引入系统, 可描述为 令 x = (x, m, L f ), p = (p, p m, p Lf ), 新的 H 函数 dl f dτ = (5.14) ( 3 ) ϵ H(τ, x, p, v) = L f ū i H i βmp m ū g + ϵg 1 ū 3 i=1 (5.15) 根据极大值原理的一阶必要条件, 最优轨迹 ( x, p, v) 满足方程 d x dτ = H ( x, p, ū) p d p dτ H = ( x, p, ū) x (5.16) 其中 ū 满足方程 (5.13), 目标是最小化 L f, 因此连续推力最小累积赤经轨道转移 问题转化为两点边值问题 系统由 (5.16) 确定, 初始和终端边界条件满足 x() = (P, e x, e y, h x, h y, m ) Φ( x) = (P P f, e x e f x, e y e f y, h x h f x, h y h f y) = 由横截条件可知另需满足边值条件 p Lf () = p Lf (f) = p m (f) = 数值仿真算例的边值条件值如表 5.1 所示 接下来, 应用单值打靶法求解两 点边值问题, 打靶函数为 S( x) = (P P f, e x e f x, e y e f y, h x h f x, h y h f y) 73

92 5.3. 数值实验 表 5.1 累积赤经最优转移数值边值 变量 初始条件 终端条件 P 15.6Mm 32.1Mm e x e y..3 h x.22. h y..1 数值仿真结果如图所示, 图 5.1 所示是最大推力.8 牛 N 对应的最优轨迹, 上图显示的是 (r 1, r 2, r 3 )3D 轨迹, 箭头表示推力的大小和方向 ; 左下图是 (r 1, r 2 ) 平面中偏心率的变化轨迹, 右下图则是 (r 2, r 3 ) 平面轨道倾角的变化轨迹 图中符号 * 表示起始点, O 则代表终端点 对于.8N 推力, 可以看出卫星完成任务使命需要绕飞地球约 88 圈 图 5.1 累积赤经最优轨道转移 3D 最优轨迹,.8 N 图 5.2 显示的是状态最优变化轨迹, 图 5.3 则是最优控制输入的结构变化 由图可以看出, 状态变量的变化平滑, 累积赤经 L 与独立变量 τ 近似成线性关系 由图 5.2, 可以看出六个变量 (P, e x, e y, h x, h y, m) 都达到了期望的目标 我 74

93 第五章近似最优赤经轨道转移 们也可以由图 5.3 看出 ū 几乎处处等于 1, 而 Ψ 在本算例中没有零点 P τ 1 m τ.4 e x τ.5 h x τ.1 e y h y τ τ L τ 图 5.2 累积赤经最优轨道转移状态的最优解,.8N 调整最大推力值, 重复数值实验, 最小累积赤经 L f 与最大推力幅值 T max 的关系如图 15 所示, 可以看出 L fmin T max 近似成常数 相同的结论在文献 [32, 112, 118] 中已提到过, 而在本文中为方便设定 L() = 在不同推力下的 σ 的变化轨迹如图 5.5 所示, 在数值上证实了 σ < 1 的假设条件 5.4 小结本节研究了连续推力最小累积赤经轨道转移问题 累积赤经 L 代替时间 t 作为系统的自变量, 应用极大值原理将轨道转移问题转化为两点边值问题, 数值实验结果表明单值打靶法能够有效解决该问题 本节的主要贡献就是应用近似最优方法系统研究了累积赤经最优轨道转移问题 在数值计算中, 我们同样得出 L fmin T max 近似成常数关系的结论, 但是并没有得到数学上的严格证明 另一个问题就是, 在数值实验中, 我们并未发现存在 Ψ = 的切换点或者奇异弧, 但是我们依然不知道是否会在其他算例中出现 Ψ = 的情况, 数值仿真发现的结论和问题将给理论研究提供了有益的启示 75

94 5.4. 小结 1 1 u 1 ψ τ τ 2 u 2 ψ τ τ 2 u 3 ψ τ τ 2 u.8 Ψ τ τ 图 5.3 累积赤经最优轨道转移控制输入,.8 N Lfmin /Tmax 图 5.4 L fmin T max 近似常数关系 76

95 第五章近似最优赤经轨道转移 图 5.5 不同推力幅值下 σ 变化轨迹 77

96

97 第六章能量 / 燃料最优轨道转移 第六章 能量 / 燃料最优轨道转移 在第三和第四章, 我们研究了时间最优轨道机动问题, 本章则应用改进轨 道根数描述的动力学方程研究燃料最优轨道转移问题 问题同样满足可控性条 件, 文献 [32, 65] 则证明了以位置 速度表示的开普勒型方程在共面情况下轨 道转移时间固定时燃料最优解的存在性 当转移时间自由且终端赤经 L 固定时 问题仍可解, 但当终端 L 也自由时虽未得到理论证明, 数值计算却表明问题无 解 [32, 83, 83] 对于固定时间轨道转移任务, 通常性能指标要求消耗燃料最少, 即终端 时刻卫星总质量最大, 但是由于奇异弧的控制不易确定, 燃料最优不易求 解 [32] 文献 [18] 则直接研究了燃料最优轨道转移, 假设不存在奇异弧即切换 函数仅存在孤立零点, 应用二分法检测切换点, 提高了打靶法的收敛性和运行 时间 文献 [15] 则基于考虑地球扁率摄动 大气摄动的开普勒型动力学模型, 提出了一种基于振荡频率分析称之为采样法的切换算法, 能够提高运算速度和 切换点的精确度, 同时滑行段的切换时间可以作为优化参数 文献 [77, 112, 118, 119] 研究了一种次黎曼 (SR) 问题, 形式为 ẋ = m i=1 u if i min tf u dt tf 根据 Maupertuis 原理,SR 问题等价于最小能量 u 2 dt 问题, 其中时间 t f 是固定的 [15-18] 卫星轨道动力学方程是典型的带漂移项的仿射输入非线性系 统, 可以考虑设法消除漂移项从而将燃料最优问题转化为能量最优问题, 文 献 [65, 77, 83, 118, 119] 采用平均化方法研究了卫星质量不变条件下的能量最优 轨道转移问题 文献 [119] 基于推力在切向 - 法向分解下的共面高斯方程, 将平 赤经作为自变量, 从而 Hamiltonian 函数关于 L 是周期的, 平均化 Hamiltonian 函 数, 将问题转化为 SR 问题 ; 对极值解进行了几何分析, 证明了共面转移的平均 化系统是完全可积的 文献 [118] 则基于改进轨道根数 (MEE) 描述方程考虑了 共面和异面转移的平均化方法, 根据黎曼测地线的存在与完备性, 给出了二维 椭圆域任意两点存在能量最优轨迹的条件 文献 [77] 则分别分析了推力在单方 向 切向 - 法向 (R,T,N) 分解 径向 - 法向 (T,C,N) 分解等不同情况的平均化系统 另一种研究燃料最优轨道转移的方法是在性能指标中引入能量指标从而消 79

98 6.1. 问题描述 除切换函数的奇异性, 常用的方法是同伦方法 文献 [32] 研究了基于 MEE 表示的轨道动力学方程燃料最优的同伦方法, 文献 [6] 则研究了开普勒方程形式 由上述分析可知, 能量最优情形相对容易求解, 应用同伦方法求解燃料最优的一个关键就是寻求同伦算子 λ 由 向 1 过渡时的零路径跟随 (Zero PathFollowing), 文献 [32] 介绍了离散连续法 分段线性法和预测校正法三种不同的零路径跟随方法 对于能量最优, 当推力小于 5N 时, 单值打靶法也难于求解, 文献 [32, 112] 引入了一种初始化同伦方法 应用间接法研究基于 MEE 描述动力学方程的能量最优轨道转移问题, 作者目前并未见相关文献 本章仍基于 MEE 描述的高斯型动力学方程, 应用间接法研究燃料最优 能量最优以及能量 - 燃料混合最优轨道转移, 分别设计了最优控制律, 并尝试在数值上对这三种形式的最优轨道转移做比较 6.1 问题描述有限推力下最小燃料消耗轨道转移问题模型可描述如下 : min J = m(t f ) t f u dt ẋ = f (x) + Tmax 3 m i=1 u if i (x) ṁ = βt max u x() = x, m() = m, Φ(x(t f )) = u 1 其中 T max 为发动机推力的最大幅值, 控制 u = (u 1, u 2, u 3 ) 为推力在轨道坐标系中方向 卫星的轨道动力学方程的状态仍采用 MEE 形式, 定义状态为 x = [P e x e y h x h y L] T 和质量 m 若推力在轨道坐标系下沿径向- 法向 (R,T,N) 标架下分解, 则定义卫星动态的四个向量场 f i (x) 仍由方程 4.2 描述 对于轨道转移任务, 要求初始轨道和目标轨道是不同的, 因而最优控制是非空的 问题满足的初值边界条件可用 MEE 描述为 (x, m ) = (P, e x, e y, h x, h y, L, m ) R 7 8

99 第六章能量 / 燃料最优轨道转移 而终端边界条件则为 Φ(x) = (P P f, e x e f x, e y e f y, h x h f x, h y h f y) = 为使问题便于求解, 我们仍做以下假设 : A1) 系统模型的状态始终满足路径约束 A = {(x, m) P >, (e x, e y ) < 1, m m s } 即卫星在椭圆域内飞行, 在地心坐标系下位置向量幅值始终大于地球半径,m s 为不携带燃料时的卫星本体质量 A2) 最终飞行时间 t f 要严格大于最短轨道转移时间 t fmin A3) 卫星在终端时的质量是自由的 至此, 我们完成了燃料最优轨道转移最优控制问题的描述, 下一节则应用极小值原理设计最优控制律 6.2 轨道转移最优控制律设计 燃料最优 燃料最优的性能指标取为 Lagrange 型, 应用极小值原理, 系 统 Hamiltonian(H) 函数为 H = (p βt max p m ) u + H + T max m 3 i=1 u ih i (6.1) 其中 p 为大于 的常数, 通常取为 1;H i = p, f i, i =,, 3 为 Hamiltonian 提 升,p, p m 分别为状态 x, m 对应的协状态 根据极值原理可知 p(t), p m (t) 不同时 为 应用 Cauchy-Schwarz 不等式, 令 Ψ = (H 1, H 2, H 3 ), 取 p = 1, 则 6.1 式有 H (1 βt max p m ) u + H + T max Ψ u (6.2) m 则当 u = αψ, α 时 6.2 取等号 因此令 z = (x, m, p, p m ), 当 Ψ = 81

100 6.2. 轨道转移最优控制律设计 (H 1, H 2, H 3 ) 不为 时, 最小化 H 函数的解可写为 定义切换函数 则当 Ψ 时, 最优控制为 u = α(z) Ψ, α(z) [, 1] (6.3) Ψ S(t, z) = 1 βt max p m T max m Ψ, S(t, z) < Ψ u (t) = α Ψ, S(t, z) = Ψ, S(t, z) > Ψ (6.4) (6.5) 其中 α [, 1] 而当 Ψ = 时, 则有 [32] u(t) S(, 1), 1 βt max p m < u(t) B(, 1), 1 βt max p m = u(t) =, 1 βt max p m > (6.6) 由上述分析可知, 最优控制函数是由 Bang-Bang 弧和奇异子弧 (Singular subarc) 所组成 上述问题中当 Ψ, S(t, z) = 时的奇异控制无法确定, 即使假设 Ψ, 即不存在奇异弧时, 在数值上应用打靶法求解燃料最优轨道转移仍然是非常困难 能量最优 数为 T 能量最优的性能指标为 u 2 dt, 同样根据极值原理, 系统 Hamiltonian 函 H = u 2 βt max p m u + H + T max m 当 Ψ = (H 1, H 2, H 3 ) 不为 时, 定义切换函数 3 i=1 u ih i (6.7) α(t, z) = βt maxp m + T max m 2 Ψ (6.8) 82

101 第六章能量 / 燃料最优轨道转移 则最优控制为 Ψ, α(t, z) > 1 Ψ u (t) = α(t, z) Ψ, α(t, z) [, 1] Ψ, α(t, z) < (6.9) 当 Ψ = 时, 由于无需处理奇异性, 问题相对燃料最优情况容易求解 由极 大值原理的一阶必要条件有 p m = Tmax Ψ, 则可知 p m 2 m 是减函数, 而根据假设 A3 由横截条件则有 p m (t f ) =, 因而 p m 是非负的, 因而切换函数 α(t, z) 也是 非负的, 则式 6.9 中的最优控制不会出现为 的情况 ϵ 算子法求解燃料最优 由于存在奇异性, 直接应用单值打靶法求解燃料最优情形比较困难, 那么 自然想到要设法寻找合适的算法消除奇异性 而能量最优控制情形不存在奇异 性, 因此可以考虑在燃料性能指标中引入能量指标, 从而消除奇异性 可以考 虑 ϵ 算子法 [53, 12], 即在性能指标的被积函数中增加 ϵ k u T u, 从而将系统的性 能指标修改成 J(u, ϵ k ) = 同样根据极值原理, 定义切换函数 T ( u + ϵ k u 2 )dt (6.1) α(t, z) = βt maxp m + Tmax m Ψ 1 2ε k (6.11) 则最优控制为 Ψ, α(t, z) > 1 Ψ u (t) = α(t, z) Ψ, α(t, z) [, 1] Ψ, α(t, z) < (6.12) 文献 [53] 证明了 ϵ 算法的收敛性, 同时也表明当 lim T k ϵ k = 时, 将收敛 于 u dt 的最小值, 实际上当 ϵ k 降低到足够小的数值时, 可以得到原奇异解的一个相当好的近似 若取 ϵ k = 1, 就是燃料 - 能量混合最优 6.3 数值计算 6.2 节应用讨论了燃料 能量最优轨道转移问题的最优控制设计问题, 应用 83

102 6.3. 数值计算 表 6.1 燃料 / 能量最优轨道转移数值边界条件 变量 初始条件 终端条件 P 15.6Mm 32.1Mm e x e y..3 h x.612. h y..1 L π 自由 m 15kg 自由 极小值原理一阶必要条件, 系统的状态和协状态方程为 ẋ = p H (t, x, m, p, p m, u ) ṁ = pm H (t, x, m, p, p m, u ) ṗ = x H (t, x, m, p, p m, u ) ṗ m = m H (t, x, m, p, p m, u ) (6.13) 而在数值算例中, 我们考虑目标轨道的 L f 自由情形, 初始轨道和目标轨道的轨道根数如表 6.1 所示 根据假设 A3, 终端卫星质量自由, 而终端 L f 也是自由的, 因此根据极值原理的横截条件, 质量和 L 对应的协变量仍需满足终端约束条件 p m (t f ) = p L (t f ) = (6.14) 至此, 根据极值原理, 系统方程 6.13, 初始边值条件和终端边界条件 6.1 以及横 截条件 6.14 将轨迹优化问题转化为两点边值问题 下面则应用单值打靶方法求 解相应的两点边值问题, 打靶映射为 [61] S(t, p ) = exp x (t, p ) x f (6.15) 84

103 第六章能量 / 燃料最优轨道转移 具体到本问题则为 S(t f, p ) = ( P (t f ) P f, e x (t f ) e f x, e y (t f ) e f y, h x h f x, h y (t f ) h f y, p m (t f ), p L (t f ) ) 而对于最短飞行时间, 第四章的数值计算结果表明最大推力与最优轨道转移时间有如下关系 t fmin (h) T max (N) c (6.16) 因此对于固定时间燃料或能量最优轨道转移, 假设 A2 可以由下式来满足 t f = c tf t fmin (6.17) 其中 c tf 为大于 1 常数 对于表 6.1 表示的边值条件, 根据数值仿真结果式 6.16 中 的常数值为 8 求解 6.15 式的经典工具是牛顿型算法, 本文应用的非线性求解器是 More JJ 等人开发的 HYBRD/HYBRJ [121], 这个求解器使得两点边值问题的变得相对易 于求解 通过求解打靶函数确定 p, 然后求解初值问题 而对于 ϵ 算法, 我们则 采用文献 [32] 提到的预测 校正零路径跟随方法,ϵ k : 1 变化, 从而寻求燃 料最优控制 当 ϵ k = 1 时, 对应的为燃料 能量混合最优情形, 而当 ϵ k = 时则 对应燃料最优情形 以下数值仿真都是最大推力幅值 T max = 5N 时的结果, 根据式 6.16 可知最短 飞行时间约为 16h, 式 6.17 中的 c tf 取为 1.2 图示中 - 线表示的是能量最优控制 情形, -. 表示的是能量 - 燃料混合最优控制情形, 而 则表示的应用 ϵ 算法求解 的燃料最优控制情形 图 6.1 所示的是状态随时间变化的最优轨迹, 由图可以看 出状态变化比较平滑, 同时也可看出三种情况下的最优轨迹几乎是重合的, 那 么可以猜测在给定初始轨道和目标轨道 给定轨道转移时间的条件下, 无论是 燃料最优还是能量最优抑或燃料 - 能量混合最优, 轨道转移的结果都是一样的 但是从图示也可看出, 终端 L 存在较大差别, 即在目标轨道上终端位置是不同 的, 图示的 L 整定为 mod(l, 2π) π 从而可以看出绕飞地球的圈数约为 圈 图 6.2 为最优控制随时间的变化规律, 由图可以看出, 卫星运行的初始几圈 1 85

104 6.3. 数值计算 4 state vs time,t max = 5 N, t f = 192h 15 P/Mm 2 m/kg time/h time/h.1 e x.5 h x time/h time/h.2 e y h y time/h time/h L/rad time/h 图 6.1 状态随时间变化轨迹 控制规律比较相似, 而 u 的变化轨迹则显示能量最优和能量 - 燃料混合最优情形都未出现为 的时刻, 最优控制是连续的 在能量最优轨道转移情形中, 我们已经证明了不会出现 u = 的情况 ; 但是在能量 - 燃料混合最优控制情形, 切换函数为 α(t, z) = βt maxp m + T max m Ψ, 并未证明是非负的 而在燃料最优控制情 2 形, 未出现奇异弧, 即仅存在有限个切换点, 最优控制是典型 Bang-Bang 控制 同时也可看出我们基于的假设 Ψ = 是成立的, 但是在理论上仍是未解决的开问题 图 6.3 图 6.4 则分别表示是能量最优控制情形和燃料最优控制情形卫星在地心坐标系下的 3D 飞行轨迹, 燃料 - 能量混合最优则与能量最优情形相似, 在此未另行给出 图中红色的箭头代表推力方向, 由图可以看出能量最优或能量 - 燃料混合最优情形发动机一直工作, 而在能量最优情形则存在零推力阶段 改变最大推力幅值重复上述数值实验, 结果如图 6.5 所示, 能量最优 燃料最优及能量 - 燃料最优三种情形下轨道转移任务完成时卫星的质量和发动机最大推力之间的关系, 由图可以看出三种情况下的最终质量几乎相同, 而对于不同的最优轨道转移任务终端质量也几乎是恒定的 可以猜测给定初始轨道和目标轨道, 发动机最大推力的大小只改变轨道转移的圈数, 即推力大则圈数少 轨 86

105 第六章能量 / 燃料最优轨道转移 u 1 u 2 u 3 u 1 Control vs time, T max = 5 N, t f = 192h time(h) time(h) time(h) time(h) 1 Ψ time(h) 图 6.2 最优控制随时间变化轨迹 道转移时间短, 而推力小则圈数多 轨道转移时间长, 并不影响燃料的消耗, 即 T max u t f 是恒定的 6.4 小结本文基于间接法分别研究了时间固定下的连续推力能量最优 燃料最优以及燃料 - 能量混合最优轨道转移控制 燃料最优因存在奇异性而不易数值求解, 转而研究能量最优和能量 - 燃料混合最优, 通过在性能指标中引入能量指标消除奇异性来研究燃料最优控制 应用极值原理将轨迹优化问题转化为两点边值问题, 采用单值打靶法求解相应的两点边值问题, 数值结果表明基于的假设 Ψ = 是成立的 ; 能量最优或能量燃料混合最优情形最优控制是连续的, 而燃料最优情形是典型的 bang-bang 控制 ; 给定初始和目标轨道, 轨道转移任务消耗的燃料是确定的, 推力大小只改变轨道转移的圈数 ; 能量最优, 燃料最优及能量 - 燃料混合最优三种不同性能指标对应的最优轨迹几乎是相同的 两点边值问题的数值求解仍是未很好解决的难题, 由于数值计算受到计算精度 算法收敛域等影响, 上述数值得到的结果并未得到理论证明, 结论是否具普遍性仍需进一步研究 87

106 6.4. 小结 图 6.3 能量最优 3D 轨迹 图 6.4 燃料最优 3D 轨迹 88