留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

滑翔飞行器线性伪谱模型预测控制三维轨迹规划

孙建波 潘幸华 杨良 陈万春 赵育善

孙建波, 潘幸华, 杨良, 陈万春, 赵育善. 滑翔飞行器线性伪谱模型预测控制三维轨迹规划[J]. 红外与激光工程, 2020, 49(9): 20200279. doi: 10.3788/IRLA20200279
引用本文: 孙建波, 潘幸华, 杨良, 陈万春, 赵育善. 滑翔飞行器线性伪谱模型预测控制三维轨迹规划[J]. 红外与激光工程, 2020, 49(9): 20200279. doi: 10.3788/IRLA20200279
Sun Jianbo, Pan Xinghua, Yang Liang, Chen Wanchun, Zhao Yushan. 3D trajectory planning for gliding vehicle using linear pseudospectral model predictive control[J]. Infrared and Laser Engineering, 2020, 49(9): 20200279. doi: 10.3788/IRLA20200279
Citation: Sun Jianbo, Pan Xinghua, Yang Liang, Chen Wanchun, Zhao Yushan. 3D trajectory planning for gliding vehicle using linear pseudospectral model predictive control[J]. Infrared and Laser Engineering, 2020, 49(9): 20200279. doi: 10.3788/IRLA20200279

滑翔飞行器线性伪谱模型预测控制三维轨迹规划

doi: 10.3788/IRLA20200279
详细信息
    作者简介:

    孙建波(1990-),男,博士生,主要从事再入轨迹规划和制导方面的研究工作。Email:sunjianboshen@163.com

    赵育善(1957-),男,教授,博士生导师,博士,主要从事航天器动力学与控制方面的研究工作。Email:yszhao@buaa.edu.cn

    通讯作者: 潘幸华(1970−),男,研究员,博士生导师,硕士,主要从事飞行器总体设计方面的研究工作。Email:2631365575@qq.com
  • 中图分类号: V448.235

3D trajectory planning for gliding vehicle using linear pseudospectral model predictive control

  • 摘要: 对高升阻比滑翔飞行器,在线性伪谱模型预测控制基础上提出新的再入制导律,除了满足传统终端约束与路径约束,还能以特定航向角抵达终点。以高维多项式代理技术泛化升阻比,得到关于能量和攻角的表达式,攻角在线调节升阻比以增强规划能力。再入飞行分为下降段和滑翔段,下降段维持最大攻角和零倾侧角以限制热流率。滑翔段应用线性伪谱模型预测控制,用降阶动力学模型预测终端偏差,线性化模型获得误差传播方程。由于积分计算复杂,以高斯伪谱法获取控制量的修正值,修正攻角、倾侧角相关参数和两次倾侧反转的能量时刻消除终端偏差。方法简单易行,精度高,便于在线应用,仿真结果显示该方法能满足提出的规划需求。
  • 图  1  升阻比曲线

    Figure  1.  Lift to drag ratio curves

    图  2  高度速度曲线

    Figure  2.  Height-vs-velocity histories

    图  3  地面轨迹曲线

    Figure  3.  Groud track curves

    图  4  航向角曲线

    Figure  4.  Heading angle histories

    表  1  初始终端状态

    Table  1.   Initial and terminal states

    CaseInitial heading angle/(°)Terminal heading angle/(°)Terminal longitude/(°)Terminal latitude/(°)
    1 90 65 120 0
    2 90 70 120 0
    3 90 75 120 0
    4 110 135 110 0
    5 110 130 110 0
    6 110 125 110 0
    下载: 导出CSV

    表  2  控制变量

    Table  2.   Control variables

    CaseAngle of attack/(°)Bank angle/(°)Reversal energy
    1 9.091 0.6514 −4.559 5E7
    2 13.53 0.7008 −4.564 6E7
    3 18.55 0.7661 −4.572 3E7
    4 11.82 0.6876 −4.258 8E7
    5 15.81 0.7970 −4.223 8E7
    6 20.64 0.7661 −4.170 40E7
    下载: 导出CSV

    表  3  终端误差

    Table  3.   Terminal error

    CaseLongitude/(°)Latitude/(°)Heading angle/(°)
    CASE-1 2.92E-10 1.53E-6 2.42E-10
    CASE-2 2.07E-10 2.70E-7 3.39E-10
    CASE-3 7.37E-11 7.00E-7 5.48E-11
    CASE-4 4.31E-10 1.07E-7 1.29E-10
    CASE-5 8.99E-9 4.21E-8 1.63E-10
    CASE-6 5.39E-9 6.65E-7 5.43E-7
    下载: 导出CSV
  • [1] 邱家稳, 王强, 马继楠. 深空探测技术(特约)[J]. 红外与激光工程, 2020, 49(5): 20201001. doi:  20201001

    Qiu Jiawen, Wang Qiang, Ma Ji'nan. Deep space exploration technology (Invited) [J]. Infrared and Laser Engineering, 2020, 49(5): 20201001. (in Chinese) doi:  20201001
    [2] Harpold J C, Graves C A. Shuttle entry guidance [J]. Journal of Astronautical Sciences, 1979, 37(3): 239-268.
    [3] Shen Z, Lu P. On-board generation of three-dimensional constrained entry trajectories [J]. Journal of Guidance, Control and Dynamics, 2003, 26(1): 111-121. doi:  10.2514/2.5021
    [4] Shen Z, Lu P. Dynamic lateral guidance logic [J]. Journal of Guidance, Control and Dynamics, 2004, 27(6): 949-959. doi:  10.2514/1.8008
    [5] Yu W, Chen W. Entry guidance with real-time planning of reference based on analytical solutions [J]. Advances in Space Research, 2015, 55(9): 2325-2345. doi:  10.1016/j.asr.2015.02.002
    [6] Yang L, Liu X, Chen W, et al. Autonomous entry guidance using Linear Pseudospectral Model Predictive Control [J]. Aerospace Science and Technology, 2018, 80(9): 38-55.
    [7] 骞微著, 杨立保. 基于小波神经网络的光纤陀螺误差补偿方法[J]. 中国光学, 2018, 11(6): 1024-1031. doi:  10.3788/co.20181106.1024

    Qian Weizhu, Yang Libao. A fiber optic gyro error compensation method based on wavelet neural network [J]. Chinese Optics, 2018, 11(6): 1024-1031. (in Chinese) doi:  10.3788/co.20181106.1024
    [8] 邵会兵, 崔乃刚, 韦常柱. 滑翔导弹末段多约束智能弹道规划[J]. 光学 精密工程, 2019, 27(2): 410-420. doi:  10.3788/OPE.20192702.0410

    Shao Huibing, Cui Naigang, Wei Changzhu. Multi-constrained intelligent trajectory planning for gliding missiles [J]. Optics and Precision Engineering, 2019, 27(2): 410-420. (in Chinese) doi:  10.3788/OPE.20192702.0410
    [9] 刘蓉, 黄大庆, 姜定国. 高超声速飞行器的反步滑模神经网络控制系统[J]. 光学 精密工程, 2019, 27(11): 2392-2401. doi:  10.3788/OPE.20192711.2392

    Liu Rong, Huang Daqing, Jiang Dingguo. Backstepping sliding mode neural network control system for hypersonic vehicle [J]. Optics and Precision Engineering, 2019, 27(11): 2392-2401. (in Chinese) doi:  10.3788/OPE.20192711.2392
    [10] 魏海坤, 徐嗣鑫, 宋文忠. 神经网络的泛化理论和泛化方法[J]. 自动化学报, 2001, 27(6): 806-815.

    Wei Haikun, Xu Sixin, Song Wenzhong. Generalization theory and generalization methods for neural networks [J]. Acta Automatica Sinica, 2001, 27(6): 806-815. (in Chinese)
    [11] 李一涵, 胡海洋, 王强. 高超声速飞行器红外探测窗口辐射透射特性研究[J]. 红外与激光工程, 2020, 49(4): 0404002.

    Li Yihan, Hu Haiyang, Wang Qiang. Radiative transmission property of infrared window in hypersonic vehicle [J]. Infrared and Laser Engineering, 2020, 49(4): 0404002. (in Chinese)
    [12] 余晓娅, 刘立拓, 李瑞, 等. 高超声速再入试验的辐射光谱定量测量[J]. 中国光学, 2020, 13(2): 87-94.

    Yu Xiaoya, Liu Lituo, Li Rui, et al. Measurements of absolute radiative emissions for supersonic reentry [J]. Chinese Optics, 2020, 13(2): 87-94. (in Chinese)
    [13] Phillips T H. A common aero vehicle (CAV) model, description, and employment Guide[R]. US: Schafer Corp for Air Force Research Laboratory and Air Force Command, 2003.
  • [1] 李云辉, 霍炬, 杨明.  基于图模型的小尺寸飞行器地面试验中位姿估计方法 . 红外与激光工程, 2020, 49(4): 0413003-0413003-8. doi: 10.3788/IRLA202049.0413003
    [2] 孙浩添, 杜福嘉, 张志永.  基于模型预测控制的大口径快摆镜随动系统 . 红外与激光工程, 2020, 49(2): 0214001-0214001. doi: 10.3788/IRLA202049.0214001
    [3] 张雄雄, 宋言明, 孟凡勇, 孙广开, 祝连庆.  变体飞行器柔性复合蒙皮植入式光纤形状传感 . 红外与激光工程, 2019, 48(6): 622003-0622003(8). doi: 10.3788/IRLA201948.0622003
    [4] 高晔, 周军, 郭建国, 冀华.  红外成像制导导弹分布式协同制导律研究 . 红外与激光工程, 2019, 48(9): 904007-0904007(9). doi: 10.3788/IRLA201948.0904007
    [5] 刘显著, 王超, 江伦, 刘壮, 杨进华, 姜会林.  二维多项式位相光瞳滤波实现超分辨望远成像 . 红外与激光工程, 2018, 47(4): 418007-0418007(6). doi: 10.3788/IRLA201847.0418007
    [6] 王宏涛, 石德平, 刘恒军.  带热流密度限制的推力可控导弹三维制导律 . 红外与激光工程, 2018, 47(3): 304004-0304004(11). doi: 10.3788/IRLA201847.0304004
    [7] 史震, 何晨迪, 郑岩.  攻角约束下的二阶滑模控制器的协同制导律设计 . 红外与激光工程, 2018, 47(6): 617006-0617006(8). doi: 10.3788/IRLA201847.0617006
    [8] 张文渊, 兰志, 夏群利, 祁载康.  导引头隔离度对最优落角制导律制导性能影响研究 . 红外与激光工程, 2017, 46(3): 331001-0331001(6). doi: 10.3788/IRLA201746.0331001
    [9] 李召龙, 沈同圣, 娄树理.  基于多项式逼近的红外系统渐晕效应校正方法 . 红外与激光工程, 2016, 45(S1): 6-10. doi: 10.3788/IRLA201645.S104002
    [10] 王辉, 武涛, 杜运理, 林德福.  无动力学滞后的广义最优制导律解析研究 . 红外与激光工程, 2015, 44(1): 341-347.
    [11] 孙文卿, 陈磊, 李金鹏, 乌兰图雅, 何勇.  非圆孔径离散采样点正交多项式波前拟合 . 红外与激光工程, 2015, 44(3): 1068-1072.
    [12] 李艳辉, 厉明, 周凌, 张楠.  基于模型匹配的光电侦察无人机飞行控制器设计方法 . 红外与激光工程, 2015, 44(2): 693-698.
    [13] 王乃祥, 徐钰蕾, 史磊, 程志峰, 姚园.  高马赫飞行器迎风面与攻角对光学窗口周围流场的影响分析 . 红外与激光工程, 2015, 44(4): 1267-1272.
    [14] 梁冬生, 刘朝晖, 刘文, 袁辉, 刘夫成.  航空飞行器天文自主导航定位技术 . 红外与激光工程, 2014, 43(9): 3020-3025.
    [15] 肖松, 谭贤四, 王红, 李志淮.  变结构多模型临近空间高超声速飞行器跟踪算法 . 红外与激光工程, 2014, 43(7): 2362-2370.
    [16] 范国伟, 邢斯瑞.  挠性卫星的姿态机动滚动优化控制 . 红外与激光工程, 2014, 43(S1): 108-115.
    [17] 路远, 冯云松, 凌永顺, 乔亚.  飞行器尾焰红外辐射及其被动测距 . 红外与激光工程, 2013, 42(7): 1660-1664.
    [18] 徐顶国, 冯维林, 桑建华.  飞行器后机身蒙皮红外辐射特性 . 红外与激光工程, 2013, 42(1): 7-13.
    [19] 王奇涛, 佟首峰, 徐友会.  采用Zernike多项式对大气湍流相位屏的仿真和验证 . 红外与激光工程, 2013, 42(7): 1907-1911.
    [20] 张跃, 储海荣.  捷联式光学导引头特性与多维度最优制导律 . 红外与激光工程, 2013, 42(11): 2967-2973.
  • 加载中
图(4) / 表(3)
计量
  • 文章访问数:  26
  • HTML全文浏览量:  24
  • PDF下载量:  10
  • 被引次数: 0
出版历程
  • 收稿日期:  2020-04-09
  • 修回日期:  2020-05-12
  • 网络出版日期:  2020-09-22
  • 刊出日期:  2020-09-22

滑翔飞行器线性伪谱模型预测控制三维轨迹规划

doi: 10.3788/IRLA20200279
    作者简介:

    孙建波(1990-),男,博士生,主要从事再入轨迹规划和制导方面的研究工作。Email:sunjianboshen@163.com

    赵育善(1957-),男,教授,博士生导师,博士,主要从事航天器动力学与控制方面的研究工作。Email:yszhao@buaa.edu.cn

    通讯作者: 潘幸华(1970−),男,研究员,博士生导师,硕士,主要从事飞行器总体设计方面的研究工作。Email:2631365575@qq.com
  • 中图分类号: V448.235

摘要: 对高升阻比滑翔飞行器,在线性伪谱模型预测控制基础上提出新的再入制导律,除了满足传统终端约束与路径约束,还能以特定航向角抵达终点。以高维多项式代理技术泛化升阻比,得到关于能量和攻角的表达式,攻角在线调节升阻比以增强规划能力。再入飞行分为下降段和滑翔段,下降段维持最大攻角和零倾侧角以限制热流率。滑翔段应用线性伪谱模型预测控制,用降阶动力学模型预测终端偏差,线性化模型获得误差传播方程。由于积分计算复杂,以高斯伪谱法获取控制量的修正值,修正攻角、倾侧角相关参数和两次倾侧反转的能量时刻消除终端偏差。方法简单易行,精度高,便于在线应用,仿真结果显示该方法能满足提出的规划需求。

English Abstract

    • 滑翔飞行器在再入过程中要经历复杂的环境条件,因此再入制导律的设计至关重要[1]。再入制导方法通常可分为两类:标称轨迹制导和预测校正制导。经典的标称轨迹制导算法为航天飞机再入制导律[2],Shen和Lu利用拟平衡滑翔假设将路径约束简化为倾侧角约束,完成了三维再入轨迹在线快速规划[3],随后为了改进侧向制导的精度,Shen采用横程参数替代航天飞机制导算法的航向瞄准误差,提出了新的倾侧角反转判定方法[4]。尽管标称轨迹制导算法得到广泛应用,但是该类算法太依赖于参考轨迹的设计,为了提高再入制导算法的灵活性和鲁棒性,研究者们转向预测校正制导方法。该算法不依赖预先设计的轨迹,而是通过在线预测终端误差修正控制参数完成任务。

      目前的预测校正算法需要利用数值计算的方式求解控制参数的梯度信息,计算效能消耗随控制参数数量成倍增加,因此该类制导算法的控制参数通常很少,这样只能满足部分终端约束。此外,基于航向角误差走廊的传统侧向制导方法,不仅减弱了高升阻比飞行器的横向机动能力,同时横向误差较大,倾侧角反转次数较多。Yu和Chen利用弹道阻尼和谱分解方法,基于解析的在线轨迹规划提出了一种预测制导算法[5],Yang利用线性伪谱方法预测终端误差设计了模型预测制导方法[6]

      前述研究中攻角都是事先规划,没有在线调整能力,终端位置通过倾侧角幅值调整和符号反转完成,横向机动的调节能力还可增强。此外,为了实现一些特殊任务,还需考虑更多的终端约束,例如为了实现特殊的攻击角度,需要考虑特定的航向角约束。因此在Yang提出的线性伪谱再入制导方法的基础上,借鉴如今在导航、制导和控制里得到广泛运用的神经网络技术[7-9],从神经网络泛化的角度出发[10],通过引入高维多项式代理技术,研究将攻角纳入规划控制变量的再入轨迹规划方法。该方法能供快速规划所需要的再入弹道,并通过调整攻角达到对横向机动轨迹进行调整的目的,从而保证能够获得所需的终端航向角约束,开展了多个不同状态的算例测试,得到了满意的仿真结果。

    • 基于圆球模型的质点再入动力学方程为:

      $$\begin{split} &\dot r = v\sin \gamma {\rm{ ;}}\dot \theta = \dfrac{{v\cos \gamma \sin \psi }}{{r\cos \phi }};\dot \phi = \dfrac{{v\cos \gamma \cos \psi }}{r} \\ & \dot v = - \dfrac{D}{m} - g\sin \gamma + {\omega ^2}r\cos \phi \left( \begin{array}{l} \sin \gamma \cos \phi \\ - \cos \gamma \cos \psi \sin \phi \\ \end{array} \right){\rm{ }} \\ & \dot \gamma = \dfrac{{L\cos \sigma }}{{mv}} - \left( {\dfrac{g}{v} - \dfrac{v}{r}} \right)\cos \gamma + 2\omega \sin \psi \cos \phi + \\ & {\rm{ }}\dfrac{{{\omega ^2}r}}{v}\cos \phi \left( {\cos \gamma \cos \phi + \sin \gamma \cos \psi \sin \phi } \right){\rm{ }} \\ & \dot \psi = \dfrac{{L\sin \sigma }}{{vm\cos \gamma }} + \dfrac{v}{r}\cos \gamma \sin \psi \tan \phi - \\ & {\rm{ }}2\omega \left( {\tan \gamma \cos \psi \cos \phi - \sin \phi } \right) + \\ & {\rm{ }}\dfrac{{{\omega ^2}r}}{{v\cos \gamma }}\sin \psi \sin \phi \cos \phi {\rm{ }} \\ \end{split} $$ (1)

      式中:$r = {R_e} + h$h为飞行器质心高度,${R_e}$为地球半径,6 378 245 m;θϕ分别为经度和纬度;v为飞行器相对地球的速度值,m/s;γ为航迹角,表示飞行器相对地球的速度与当地水平面的夹角;ψ为航向角,表示飞行器相对地球的速度在当地水平面的投影与正北方向的夹角,顺时针为正;m为飞行器质量,kg;$g = \mu /{r^2}$表示当地重力加速度,m/s2μ为地球重力常数;σ为倾侧角,表示飞行器沿速度方向旋转的角度;ω为地球自转角速度,7.292 1×10−5 rad/s;LD分别表示飞行器受的升力和阻力:

      $$L = \frac{1}{2}\rho {v^2}{C_l}{S_{ref}};D = \frac{1}{2}\rho {v^2}{C_d}{S_{ref}}$$ (2)

      $\rho = {\rho _0}\exp \left( { - h/H} \right)$为大气密度,${\rho _0}$为海平面标准大气压,1.225 kg/m3H为大气密度常数,7200;Sref为飞行器特征面积;ClCd分别表示飞行器的升力、阻力系数,是关于马赫数和攻角的函数。

    • 路径约束包括热流密度[11-12]、动压和过载约束:

      $$\begin{array}{l} \dot Q = {K_Q}{\left( \rho \right)^{0.5}}{V^{3.15}} \leqslant {{\dot Q}_{\max }} \\ q = 0.5\rho {V^2} \leqslant {q_{\max }} \\ n = \dfrac{{\sqrt {{L^2} + {D^2}} }}{{m{g_0}}} \leqslant {n_{\max }} \\ \end{array} $$ (3)

      式中:$\dot Q$为热流密度;KQ为常数;$q$为动压;$n$为过载;${g_0}$为海平面重力加速度;${\dot Q_{\max }}$${q_{\max }}$${n_{\max }}$分别为飞行器能承受的最大热流、动压和过载。

      终端状态约束如下,新的航向角约束是为了实现特殊的飞行任务:

      $$\left\{ \begin{array}{l} h\left( {{E_f}} \right) = {h_f},v\left( {{E_f}} \right) = {v_f},\sigma \left( {{E_f}} \right) = 0 \\ \lambda \left( {{E_f}} \right) = {\lambda _f},\phi \left( {{E_f}} \right) = {\phi _f},\psi \left( {{E_f}} \right) = {\psi _f} \\ \end{array} \right.$$ (4)

      Ef为终端能量,$E = - \mu /r + {v^2}/2$;其余以f为下标的状态量为终端状态,终端倾侧角为后面的飞行保留调整余量。

    • 以CAV-H作为研究模型,质量为907 kg,特征面积为0.483 9 m2,最大升阻比接近3.5[13]ClCd用高阶多项式拟合,攻角的调节范围是5 ~ 20°。最大热流、动压和过载分别为500 W/cm2、100 kPa和2。攻角为分段函数:

      $$\left\{ \begin{array}{l} \alpha = {\alpha _{\rm{b}}},E > {E_{{\rm{mid}}}} \\ \alpha = f\left( E \right)\left( {{\alpha _{\rm{b}}} - {\alpha _f}} \right) + {\alpha _f},{E_f} \leqslant E \leqslant {E_{{\rm{mid}}}} \\ \end{array} \right.$$ (5)
      $$f\left( E \right) = \left( {\frac{{2\left( {E - {E_f}} \right)}}{{\left( {{E_{{\rm{mid}}}} - {E_f}} \right)}} - \frac{{{{\left( {E - {E_f}} \right)}^2}}}{{{{\left( {{E_{{\rm{mid}}}} - {E_f}} \right)}^2}}}} \right)$$ (6)

      式中:Emid以高度35 km和速度3 500 m/s计算,接近${E_f}$${\alpha _{\rm{b}}}$为控制参数;为了充分发挥攻角的调节能力,${\alpha _f}$为5°。

    • 针对不同的常值攻角,结合能量方程和平衡滑翔条件,联立求解可获得高度速度曲线和升阻比能量曲线。${\sigma _{\rm{0}}} = {30^ \circ }$,倾侧角留有调节余量。

      $$\left\{ \begin{array}{l} E{\rm{ = }} - \mu {\rm{/}}r + {v^2}/2 \\ 0 = L\cos {\sigma _{\rm{0}}}{\rm{/}}\left( {mv} \right) - \left( {g{\rm{/}}v - v/r} \right) \\ \end{array} \right.$$ (7)

      分别对攻角在${10^ \circ }$${15^ \circ }$${20^ \circ }$的情况进行仿真计算,如图1图2所示。

      图  1  升阻比曲线

      Figure 1.  Lift to drag ratio curves

      图  2  高度速度曲线

      Figure 2.  Height-vs-velocity histories

      升阻比可表示为以能量E和攻角α为自变量的函数,采用高维多项式拟合获得满足精度的解析函数:

      $$F\left( {\alpha ,E} \right) = {\left[ {\begin{array}{*{20}{c}} 1 \\ \alpha \\ {{\alpha ^2}} \end{array}} \right]^{\rm{T}}}\left[ {\begin{array}{*{20}{c}} {{p_{11}}}&{{p_{12}}}& \cdots &{{p_{16}}} \\ {{p_{21}}}&{{p_{22}}}& \cdots &{{p_{26}}} \\ {{p_{31}}}&{{p_{32}}}& \cdots &{{p_{36}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} 1 \\ \vdots \\ {{E^5}} \end{array}} \end{array}} \right]$$ (8)
    • 高升阻比飞行器再入飞行分为下降段和滑翔段。每个阶段的制导策略在后文详述。

    • 下降段从初始高度开始,到飞行器能够产生足够的升力保持滑翔时结束。此阶段飞行器速度很快,大气密度急速增加,为了限制最大热流,选择最大容许攻角和零倾侧角的控制指令。

      $${\alpha _{{\rm{cmd}}}} = {\alpha _{\max }};{\sigma _{{\rm{cmd}}}} = 0$$ (9)

      $\dot h = v\sin \gamma = 0$时,下降段结束,滑翔段开始。

    • 滑翔段以多段线性伪谱法更新倾侧反转点和控制量,引导飞行器飞向交班点,并为后续飞行保留足够的能量。

    • 根据平衡滑翔假设简化动力学方程得到如下方程,航向角方程中第一项的符号随倾侧角反转而改变。记${{x}} = {\left[ {\theta ,\phi ,\varphi } \right]^{\rm{T}}}$$\dot {{x}} = {{f}}$,根据符号的变化可将降阶动力学方程分为三段,分别记作${{{f}}_1}$${{{f}}_2}$${{{f}}_3}$

      $$\left\{ \begin{array}{l} \dfrac{{\rm{d}\theta }}{{{\rm{d}}E}}{\rm{ = }} - \dfrac{{\sin \psi }}{{{\varTheta }\left( {E,\psi ,\phi } \right)r\cos \phi }}F\left( {\alpha ,E} \right)u;{\rm{ }} \\ \dfrac{{\rm{d}\phi }}{{{\rm{d}}E}}{\rm{ = }} - \dfrac{{\cos \psi }}{{{\varTheta }\left( {E,\psi ,\phi } \right)r}}F\left( {\alpha ,E} \right)u; \\ \dfrac{{\rm{d}\psi }}{{{\rm{d}}E}}{\rm{ = }} \mp \dfrac{1}{{2\left( {E + \mu /r} \right)}}F\left( {\alpha ,E} \right)\sqrt {1 - {u^2}} - \\ \left( \begin{array}{l} \dfrac{{\sin \psi \tan \phi }}{{{\varTheta }\left( {E,\psi ,\phi } \right)r}}{\rm{ + }}\dfrac{{2{\omega _e}\sin \phi }}{{{\varTheta }\left( {E,\psi ,\phi } \right)\sqrt {2\left( {E + \mu /r} \right)} }}{\rm{ + }} \\ \dfrac{{\omega _e^2r\sin \psi \sin \phi \cos \phi }}{{2\left( {E + \mu /r} \right){\varTheta }\left( {E,\psi ,\phi } \right)}} \\ \end{array} \right)F\left( {\alpha ,E} \right)u \\ \end{array} \right.$$ (10)
      $$\begin{split} {\varTheta }\left( {E,\psi ,\phi } \right) =& - \left( {\dfrac{{2E}}{r} + \frac{\mu }{{{r^2}}}} \right) - \omega _e^2r{\cos ^2}\phi - \\ &2\sqrt {2\left( {E + \dfrac{\mu }{r}} \right)} {\omega _e}\sin \psi \cos \phi \\ \end{split} $$ (11)

      其中$u = \cos \sigma $,升阻比为F$h$的取值为当前高度与${h_f}$的平均值。${\varTheta }$关于$\phi $$\psi $的偏导数如下:

      $$\begin{array}{l} \dfrac{{\rm{d}{\varTheta }}}{{\rm{d}\phi }} = 2\sqrt {2\left( {E + \dfrac{\mu }{r}} \right)} {\omega _e}\sin \psi \sin \phi + 2\omega _e^2r\sin \phi \cos \phi \\ \dfrac{{\rm{d}{\varTheta }}}{{\rm{d}\psi }} = - 2\sqrt {2\left( {E + \dfrac{\mu }{r}} \right)} {\omega _e}\cos \psi \cos \phi \\ \end{array} $$ (12)
    • 为了既能消除终端误差又能保持较大的横向机动能力,设置两个倾侧反转点Ere1Ere2。先令${E_{re2}} = {E_{mid}}$,计算Ere1;第一次反转完成后,再计算Ere2$u$用分段函数来表征:

      $$\left\{ \begin{array}{l} u = {u_{\rm{b}}},E > {E_{{\rm{mid}}}} \\ u = f\left( E \right)\left( {{u_{\rm{b}}} - 1} \right) + 1,{E_f} \leqslant E \leqslant {E_{{\rm{mid}}}} \\ \end{array} \right.$$ (13)

      因此全段的控制就可以用参数${u_{\rm{b}}}$${\alpha _{\rm{b}}}$${E_{re{\rm{1}}}}$${E_{re2}}$来表征。

    • 对于降阶动力学系统,给定初始控制参数,可以通过数值积分预测飞行轨迹,从而得到相对终端状态约束的偏差$\delta {x_f} = x({E_f}) - {x_f}$。将降阶动力学方程在标称轨迹附近泰勒展开,忽略高阶项,令${U_{\rm{b}}} = {\left[ {{u_{\rm{b}}},{\alpha _{\rm{b}}}} \right]^{\rm{T}}}$得到以状态偏差为自变量的线性动力学方程如下:

      $$\left\{ \begin{array}{l} \delta {\dot{ x}} = {{{A}}_1}\left( E \right)\delta {{x}} + {{{B}}_1}\left( E \right)\delta {U_{\rm{b}}},E > {E_{re1}} \\ \delta {\dot{ x}} = {{{A}}_2}\left( E \right)\delta {{x}} + {{{B}}_2}\left( E \right)\delta {U_{\rm{b}}},{E_{re2}} \leqslant E \leqslant {E_{re1}} \\ \delta {\dot{ x}} = {{{A}}_3}\left( E \right)\delta {{x}} + {{{B}}_3}\left( E \right)\delta {U_{\rm{b}}},{E_f} \leqslant E \leqslant {E_{re2}} \\ \end{array} \right.$$ (14)

      其中${{x}} = {{{x}}_{ref}} - \delta {{x}}$${U_{\rm{b}}} = {U_{ref}} - \delta {U_{\rm{b}}}$

      $${{{B}}_i}\left( E \right) = \left\{ \begin{array}{l} {{B}}\left( E \right),E > {E_{{\rm{mid}}}}\\ {{B}}\left( E \right)f\left( E \right)\left( {\begin{array}{*{20}{c}} 1&0\\ 0&1 \end{array}} \right),\;{E_f} \le E \le {E_{{\rm{mid}}}} \end{array} \right.$$ (15)

      ${{{A}}_i}\left( E \right)$是一个3×3的矩阵,${{B}}\left( E \right)$是一个3×2的矩阵。与参考文献[6]中系数不同,此处的系数考虑了攻角变化,${{{A}}_i}\left( E \right)$${{B}}\left( E \right)$的表达式表示为:

      $${{{A}}_i}\left( E \right) = \left[ {\begin{array}{*{20}{c}} 0&{{a_{12}}}&{{a_{13}}}\\ 0&{{a_{22}}}&{{a_{23}}}\\ 0&{{a_{32}}}&{{a_{33}}} \end{array}} \right],\;{{B}}\left( E \right) = \left[ {\begin{array}{*{20}{c}} {{b_{{\rm{1}}1}}}&{{b_{{\rm{12}}}}}\\ {{b_{2{\rm{1}}}}}&{{b_{{\rm{22}}}}}\\ {{b_{3{\rm{1}}}}}&{{b_{3{\rm{2}}}}} \end{array}} \right]$$ (16)
      $${a_{12}} = \dfrac{{\sin \psi }}{r}\left( {\dfrac{{\partial {\varTheta }}}{{\partial \phi }}/\left( {{{\varTheta }^2}\cos \phi } \right) - \dfrac{{\sin \phi }}{{{\varTheta }{{\cos }^2}\phi }}} \right)F\left( \alpha \right)u$$ (17)
      $${a_{13}} = - \dfrac{1}{{r\cos \phi }}\left( {\dfrac{{\cos \psi }}{{\varTheta }} - \dfrac{{\partial {\varTheta } }}{{\partial \psi }}\dfrac{{\sin \psi }}{{{{\varTheta } ^2}}}} \right)F\left( \alpha \right)u$$ (18)
      $${b_{{\rm{1}}1}} = - \dfrac{{\sin \psi }}{{{\varTheta }r\cos \phi }}F\left( \alpha \right)$$ (19)
      $${b_{{\rm{12}}}} = - \frac{{\sin \psi }}{{{\varTheta } r\cos \phi }}\frac{{\partial F}}{{\partial \alpha }}u$$ (20)
      $${a_{22}} = \frac{{\cos \psi }}{{r{{\varTheta }^2}}}\frac{{\partial {\varTheta } }}{{\partial \phi }}F\left( \alpha \right)u$$ (21)
      $${a_{23}} = \dfrac{1}{r}\left( {\dfrac{{\sin \psi }}{{\varTheta }} + \dfrac{{\cos \psi }}{{{{\varTheta }^2}}}\dfrac{{\partial {\varTheta } }}{{\partial \psi }}} \right)F\left( \alpha \right)u$$ (22)
      $${b_{2{\rm{1}}}} = - \dfrac{{\cos \psi }}{{{\varTheta }r}}F\left( \alpha \right)$$ (23)
      $${b_{2{\rm{2}}}} = - \dfrac{{\cos \psi }}{{{\varTheta }r}}\frac{{\partial F}}{{\partial \alpha }}u$$ (24)
      $${a_{32}} = \left( \begin{array}{l} \left( \begin{array}{l} \dfrac{{\sin \psi }}{{r\cos \phi }} + \dfrac{{2{\omega _e}}}{{\sqrt {2\left( {E + \mu {\rm{/}}r} \right)} }} + \\ \dfrac{1}{2}\dfrac{{\omega _e^2r\sin \psi \cos \phi }}{{\left( {E + \mu {\rm{/}}r} \right)}} \end{array} \right)\dfrac{{\sin \phi }}{{{{\varTheta } ^2}}}\dfrac{{\partial {\varTheta }}}{{\partial \phi }}\\ - \left( \begin{array}{l} \dfrac{{\sin \psi }}{{r{{\cos }^2}\phi }}{\rm{ + }}\dfrac{{2{\omega _e}\cos \phi }}{{\sqrt {2\left( {E + \mu {\rm{/}}r} \right)} }}、{\rm{ + }}\\ \dfrac{1}{2}\dfrac{{\omega _e^2r\sin \psi \cos 2\phi }}{{\left( {E + \mu {\rm{/}}r} \right)}} \end{array} \right)\dfrac{1}{{\varTheta }} \end{array} \right)F\left( \alpha \right)u$$ (25)
      $${a_{33}} = \left( \begin{array}{l} \left( \begin{array}{l} \dfrac{{\sin \psi }}{{\cos \phi r}} + \dfrac{{2{\omega _e}}}{{\sqrt {{\rm{2}}\left( {E + \mu {\rm{/}}r} \right)} }} \\ \dfrac{{\omega _e^2r\cos \phi \sin \psi }}{{2\left( {E + \mu {\rm{/}}r} \right)}} \\ \end{array} \right)\dfrac{{\sin \phi }}{{{{\varTheta } ^2}}}\dfrac{{\partial {\varTheta }}}{{\partial \psi }} \\ - \left( {\dfrac{1}{{r\cos \phi }} + \dfrac{{\omega _e^2r\cos \phi }}{{2\left( {E + \mu {\rm{/}}r} \right)}}} \right)\dfrac{{\sin \phi \cos \psi }}{{\varTheta } } \\ \end{array} \right)F\left( \alpha \right)u$$ (26)
      $$\begin{array}{l} {b_{3{\rm{1}}}} = \pm \dfrac{{F\left( \alpha \right)u}}{{2\left( {E + \mu /r} \right)\sqrt {{\rm{1}} - {u^2}} }} - \\ {\rm{ }}\left( \begin{array}{l} \dfrac{{\sin \psi \tan \phi }}{{{\varTheta }r}}{\rm{ + }}\dfrac{{2{\omega _e}\sin \phi }}{{{\varTheta }\sqrt {2\left( {E + \mu {\rm{/}}r} \right)} }} \\ {\rm{ + }}\dfrac{{\omega _e^2r\sin \psi \sin \phi \cos \phi }}{{2\left( {E + \mu /r} \right){\varTheta }}} \\ \end{array} \right)F\left( \alpha \right) \\ \end{array} $$ (27)
      $$\begin{array}{l} {b_{3{\rm{2}}}} = \mp \dfrac{{\sqrt {1 - {u^2}} }}{{2\left( {E + \mu /r} \right)}}\dfrac{{\partial F}}{{\partial \alpha }}\sqrt {{\rm{1}} - {u^2}} -\\ \left( \begin{array}{l} \dfrac{{\sin \psi \tan \phi }}{{{\varTheta }r}}{\rm{ + }}\dfrac{{2{\omega _e}\sin \phi }}{{{\varTheta }\sqrt {2\left( {E + \mu {\rm{/}}r} \right)} }} \\ {\rm{ + }}\dfrac{{\omega _e^2r\sin \psi \sin \phi \cos \phi }}{{2\left( {E + \mu /r} \right){\varTheta }}} \\ \end{array} \right)\dfrac{{\partial F}}{{\partial \alpha }}u \\ \end{array} $$ (28)
    • 终端状态可以表示为控制参数的函数:

      $$\begin{array}{l} {{x}}\left( {{E_f}} \right) = {{g}}\left( {{U_{\rm{b}}},{E_{re1}},{E_{re2}}} \right) = {{{x}}_0} + \\ \displaystyle\int_{{E_0}}^{{E_{re1}}} {{{{f}}_1}\left( {{{x}},{U_{\rm{b}}}} \right){\rm{d}}E + } \int_{{E_{re1}}}^{{E_{re2}}} {{{{f}}_2}\left( {{{x}},{U_{\rm{b}}}} \right){\rm{d}}E + } \int_{{E_{re2}}}^{{E_f}} {{{{f}}_3}\left( {{{x}},{U_{\rm{b}}}} \right){\rm{d}}E} \\ \end{array} $$ (29)

      $\delta U = {\left[ {\delta {U_{\rm{b}}},\delta {E_{re1}},\delta {E_{re2}}} \right]^{\rm{T}}}$,假设$\delta U$为小量,将终端状态在标称轨迹附近泰勒展开并忽略高阶项可得:

      $$\delta {{{x}}_f} = \frac{{\partial {{g}}}}{{\partial U}}\delta U$$ (30)
      $$\dfrac{{\partial {{g}}}}{{\partial U}} = \left[ {\begin{array}{*{20}{c}} {\dfrac{{\partial {g_\theta }}}{{\partial {u_{\rm{b}}}}}}&{\dfrac{{\partial {g_\theta }}}{{\partial {\alpha _{\rm{b}}}}}}&{\dfrac{{\partial {g_\theta }}}{{\partial {E_{re1}}}}}&{\dfrac{{\partial {g_\theta }}}{{\partial {E_{re2}}}}} \\ {\dfrac {\partial {g_\phi }}{{\partial{u_{\rm{b}}}}}}&{\dfrac{{\partial {g_\phi }}}{{\partial {\alpha _{\rm{b}}}}}}&{\dfrac{{\partial {g_\phi }}}{{\partial {E_{re1}}}}}&{\dfrac{{\partial {g_\phi }}}{{\partial {E_{re2}}}}} \\ {\dfrac{{\partial {g_\psi }}}{{\partial {u_{\rm{b}}}}}}&{\dfrac{{\partial {g_\psi }}}{{\partial {\alpha _{\rm{b}}}}}}&{\dfrac{{\partial {g_\psi }}}{{\partial {E_{re1}}}}}&{\dfrac{{\partial {g_\psi }}}{{\partial {E_{re2}}}}} \end{array}} \right]$$ (31)
      $$\delta U = {\left( {\frac{{\partial {{g}}}}{{\partial U}}} \right)^{ - 1}}\delta {{{x}}_f}$$ (32)

      根据公式(29)直接计算控制量的修正值需要积分运算,很难得到解析解。用拉格朗日插值多项式来拟合状态量和控制量,可将微分动力学方程转化为一组线性代数方程,从而简化计算过程[6]。结合偏差微分方程(14)与高斯伪谱法得到终端偏差的表达式为:

      $$\begin{split} \delta {{{x}}_f} =& {{K}}_3^x{{K}}_2^x{{K}}_1^x\delta {x_0} + {{K}}_3^x{{K}}_2^x\delta {{{x}}_{re1}} + {{K}}_3^x\delta {{{x}}_{re2}}+\\ & \left( {{{K}}_3^x{{K}}_2^x{{K}}_1^u + {{K}}_3^x{{K}}_2^u + {{K}}_3^u} \right)\delta {U_{\rm{b}}} \\ \end{split} $$ (33)

      $\delta {{{x}}_0}$为初始状态偏差,轨迹规划从当前状态开始,因此$\delta {{{x}}_0} = 0$。关于控制参量的偏导数为:

      $${\left[ {\begin{array}{*{20}{c}} {\dfrac{{\partial {{g}}}}{{\partial {u_{\rm{b}}}}}}&{\dfrac{{\partial {{g}}}}{{\partial {\alpha _{\rm{b}}}}}} \end{array}} \right]^{\rm{T}}} = {{K}}_3^x{{K}}_2^x{{K}}_1^u + {{K}}_3^x{{K}}_2^u + {{K}}_3^u$$ (34)

      $\delta {E_{re1}}$$\delta {E_{re2}}$分别引起${E_{re1}}$${E_{re2}}$时刻的状态偏差$\delta {{{x}}_{re1}}$$\delta {{{x}}_{re2}}$,结合公式(29)由变分法得到:

      $$\delta {{{x}}_{re1}} = \left[ {{{{f}}_1}\left( {{U_{\rm{b}}},{E_{re1}},{E_{re2}}} \right)} \right.\left. { - {{{f}}_2}\left( {{U_{\rm{b}}},{E_{re1}},{E_{re2}}} \right)} \right]\delta {E_{re1}}$$ (35)
      $$\delta {{{x}}_{re2}} = \left[ {{{{f}}_2}\left( {{U_{\rm{b}}},{E_{re1}},{E_{re2}}} \right)} \right.\left. { - {{{f}}_3}\left( {{U_{\rm{b}}},{E_{re1}},{E_{re2}}} \right)} \right]\delta {E_{re2}}$$ (36)
      $$\frac{{\partial {{g}}}}{{\partial {E_{re1}}}} = {{K}}_3^x{{K}}_2^x\left[ {{{{f}}_1}\left( {{U_{\rm{b}}},{E_{re1}},{E_{re2}}} \right)} \right.\left. { - {{{f}}_2}\left( {{U_{\rm{b}}},{E_{re1}},{E_{re2}}} \right)} \right]$$ (37)
      $$\frac{{\partial {{g}}}}{{\partial {E_{re2}}}} = {{K}}_3^x\left[ {{{{f}}_2}\left( {{U_{\rm{b}}},{E_{re1}},{E_{re2}}} \right)} \right.\left. { - {{{f}}_3}\left( {{U_{\rm{b}}},{E_{re1}},{E_{re2}}} \right)} \right]$$ (38)
      $$\begin{array}{l} {{K}}_i^x = \left\{ {{{1}} - {\rm{d}}{E_i}{{{1}}^{\rm{T}}}{{{A}}_i}{{{W}}_i}{{\left[ {{{{D}}_{2:n}} - {\rm{d}}{E_i}{{{A}}_i}} \right]}^{ - 1}}{{{D}}_1}} \right\} \\ {{K}}_i^u = \left\{ \begin{array}{l} {\rm{d}}{E_i}^2{{{1}}^{\rm{T}}}{{{A}}_i}{{{W}}_i}{\left[ {{{{D}}_{2:n}} - {\rm{d}}{E_i}{{{A}}_i}} \right]^{ - 1}}{{{B}}_i} \\ + {\rm{d}}{E_i}{{B}}_i^{\rm{T}}{{{W}}_i}{{1}} \\ \end{array} \right\} \\ \end{array} $$ (39)
      $${{{A}}_i} = \left[ {\begin{array}{*{20}{c}} {{{{A}}_i}\left( {{E_1}} \right)}& \cdots &0 \\ \vdots & \ddots & \vdots \\ 0& \cdots &{{{{A}}_i}\left( {{E_n}} \right)} \end{array}} \right],{{{B}}_i} = \left[ {\begin{array}{*{20}{c}} {{{{B}}_i}\left( {{E_1}} \right)} \\ \vdots \\ {{{{B}}_i}\left( {{E_n}} \right)} \end{array}} \right]$$ (40)

      ${E_{i{\rm{ - 1}}}}$${E_i}$为该段初始与终端能量,${\rm{d}}{E_i}{\rm{ = }}\left( {{E_i} - {E_{i - 1}}} \right)/2$${E_1}$~${E_n}$为配点处能量,${{{W}}_i}$为高斯积分权重矩阵,${{{D}}_1},{{{D}}_{2:n}}$为高斯微分近似矩阵。控制修正量为:

      $$\left[ {\begin{array}{*{20}{c}} {\delta {u_{\rm{b}}}} \\ {\delta {\alpha _{\rm{b}}}} \\ {\delta {E_{rei}}} \end{array}} \right] = {\left[ {\begin{array}{*{20}{c}} {\dfrac{{\partial {g_\theta }}}{{\delta {u_{\rm{b}}}}}}&{\dfrac{{\partial {g_\theta }}}{{\partial {\alpha _{\rm{b}}}}}}&{\dfrac{{\partial {g_\theta }}}{{\partial {E_{rei}}}}} \\ {\dfrac{{\partial {g_\phi }}}{{\delta {u_{\rm{b}}}}}}&{\dfrac{{\partial {g_\phi }}}{{\partial {\alpha _{\rm{b}}}}}}&{\dfrac{{\partial {g_\phi }}}{{\partial {E_{rei}}}}} \\ {\dfrac{{\partial {g_\psi }}}{{\delta {u_{\rm{b}}}}}}&{\dfrac{{\partial {g_\psi }}}{{\partial {\alpha _{\rm{b}}}}}}&{\dfrac{{\partial {g_\psi }}}{{\partial {E_{rei}}}}} \end{array}} \right]^{ - 1}}\left[ {\begin{array}{*{20}{c}} {\delta {\theta _f}} \\ {\delta {\psi _f}} \end{array}} \right]$$ (41)

      其中i=1,2,迭代更新控制参数直至满足终端约束,调节${E_{re1}}$时,${E_{re2}} = {E_{{\rm{mid}}}}$

    • 路径约束可以转化为对倾侧角大小的约束:

      $$\begin{array}{l} \sigma \left( v \right) \leqslant {\sigma _{\max }}\left( v \right) = \\ \arccos \left( {\dfrac{{m\left( {g\left( {{r_{\min }}} \right) - {v^2}{\rm{/}}{r_{\min }}} \right)}}{{{k_L}L\left( {{r_{\min }}} \right)}} + {k_{cst}}\left( {\dfrac{{{\rm{d}}r}}{{{\rm{d}}v}} - {{\left. {\dfrac{{{\rm{d}}r}}{{{\rm{d}}v}}} \right|}_{\dot Q,q,n}}} \right)} \right) \\ \end{array} $$ (42)

      式中:${r_{\min }}$为给定速度下满足约束的最小高度;${k_L} = \tilde L/L$$L$为根据模型计算出的升力,$\tilde L$为测量得到的升力;${\rm{d}}r/{\rm{d}}v$$\left. {\left( {{\rm{d}}r/{\rm{d}}v} \right)} \right|\dot Q,q,n$分别为再入动力学和过程约束方程中高度对速度的导数;kcst是一个常数,当存在约束没有满足时大于零,其余情况下等于零[6]

    • 新的轨迹规划方法考虑了包含地球自转项在内的多种非线性项,通过将攻角引入规划变量,极大地增强了飞行器的调节能力,下面结合不同的初始和终端条件,开展仿真验证。初始状态为:高度55 km,速度 7 000 m/s;终端约束为:高度30 km,速度2 500 m/s,初始航向角、终端位置和航向角根据任务需求具体确定,其具体取值参考表1

      表 1  初始终端状态

      Table 1.  Initial and terminal states

      CaseInitial heading angle/(°)Terminal heading angle/(°)Terminal longitude/(°)Terminal latitude/(°)
      1 90 65 120 0
      2 90 70 120 0
      3 90 75 120 0
      4 110 135 110 0
      5 110 130 110 0
      6 110 125 110 0

      值得注意的是CASE-1,2,3选择初始倾侧角符号为1,而CASE-4,5,6,选择初始倾侧角符号为−1。所有算例的计算时间均小于0.1 s,运行环境为频率3.3 GHz、内存16 G计算机下的MATLAB2016b。图3为地面轨迹曲线,两类算例具有不同的初始航向角和飞行方向,但均满足终端位置要求,飞行路径具有较大的横向机动轨迹,通过调节攻角能够对横向轨迹进行一定的调整。图4为航向角曲线,两类算例的初始航向角相同,经过两次倾侧反转都满足所需要的终端航向角,从航向角曲线也能清晰的看出各算例具有不同的倾侧反转时刻。

      图  3  地面轨迹曲线

      Figure 3.  Groud track curves

      图  4  航向角曲线

      Figure 4.  Heading angle histories

      计算所获得控制变量如表2所示,攻角都在5 ~ 20°内变化,且根据终端航向角的需求不同,将攻角的调节能力完全发挥出来,航向角和第一次反转时刻也都确定。

      表 2  控制变量

      Table 2.  Control variables

      CaseAngle of attack/(°)Bank angle/(°)Reversal energy
      1 9.091 0.6514 −4.559 5E7
      2 13.53 0.7008 −4.564 6E7
      3 18.55 0.7661 −4.572 3E7
      4 11.82 0.6876 −4.258 8E7
      5 15.81 0.7970 −4.223 8E7
      6 20.64 0.7661 −4.170 40E7

      终端误差如表3所示,可见该算法具有很高的计算精度。综合而言,所提出的方法将攻角纳入控制量,能有效的调节横向机动轨迹,在保证终端位置的条件下,有效调节终端航向角。

      表 3  终端误差

      Table 3.  Terminal error

      CaseLongitude/(°)Latitude/(°)Heading angle/(°)
      CASE-1 2.92E-10 1.53E-6 2.42E-10
      CASE-2 2.07E-10 2.70E-7 3.39E-10
      CASE-3 7.37E-11 7.00E-7 5.48E-11
      CASE-4 4.31E-10 1.07E-7 1.29E-10
      CASE-5 8.99E-9 4.21E-8 1.63E-10
      CASE-6 5.39E-9 6.65E-7 5.43E-7
    • 为了满足特定的终端航向角约束,基于线性伪谱模型预测控制技术研究了将攻角纳入控制变量的三维再入轨迹规划算法。该方法首先采用高维多项式代理技术,并结合拟平衡滑翔条件泛化升阻比,获得升阻比关于攻角和能量的泛化函数。随后对再入动力学方程进行降阶处理,保留地球自转等非线性影响因素,并对降阶后的动力学方程进行线性化处理获得误差传播方程,推导了在伪谱离散条件下关于终端位置和航向角偏差的攻角、倾侧角和倾侧角时刻解析修正关系,最终实现轨迹的在线解算。同时由于攻角参与调节,极大增强了再入轨迹的横向机动调整能力。几类典型算的仿真结果表明:该方法不仅能够保证终端和航向角,还具有很快的求解速度和很高的计算精度,适合在线应用。

参考文献 (13)

目录

    /

    返回文章
    返回