基于解析法最优轨迹的在线实时制导律研究

张志国，余梦伦

（北京宇航系统工程研究所，北京，100076）

论文与报告

基于解析法最优轨迹的在线实时制导律研究

张志国，余梦伦

（北京宇航系统工程研究所，北京，100076）

新型制导系统单纯采用轨迹优化的数值算法的优点是精度高，但耗时长，算法稳定性有待提高，随着计算机技术的发展，耗时长的缺点能得到缓解，但仍然无法满足实时动态制导控制需求。数值法以前，解析法曾得到深入的研究，如果能结合二者的优点，将能更好地达到提升飞行器实时高精度制导控制的效果。用解析法求解了轨迹优化问题中的两点边值问题，并和数值打靶法对比分析了入轨精度和计算时间的差异，得出各自方法的优缺点。

实时在线制导；解析法；轨迹优化

0 引 言

飞行器轨迹优化和制导控制一直是飞行器总体设计的主要研究方向。轨迹优化研究中常用到的方法有间接法、直接法和混合法，目前都是依赖数值算法[1]。但随着飞行器性能的日益提升，更高任务要求的不断提出，对于飞行器轨迹设计和制导性能提出了更高的要求[1]。例如，飞行轨迹实时在线规划问题[2]，预测制导问题[3]等，在实际飞行过程中需要迅速给出最优控制策略，耗时长的数值算法是无法满足实时动态制导控制需求的。虽然随着计算机技术的发展，数值算法耗时长的缺点能得到缓解，但仍然无法满足实时动态制导控制的需求，算法稳定性和鲁棒性也是需要注意的问题。解析法优点是计算时间短、稳定性高，但是精度高低取决于简化程度、展开阶数等因素。若将数值法和解析法相结合考虑，将能更好地达到提升飞行器实时高精度控制的效果。

Hu Steve S[4]在1966年给出了一套解析法求解最优飞行导引方程的方法，但是并没有给出具体的求解过程和精度分析。Gilchrist C和Morrow L[5]随后对制导理论的解析进行了总结，对于解析法求解两点边值问题给出了详细的论述，重点分析了解析法展开到不同阶数时协态变量相对标准值的偏差。解析求解方法后续在多项式制导和迭代制导的研究中被引用。随着计算机技术的发展，数值方法获得了迅速发展，如目前求解两点边值问题的打靶法。

本文中提出了求解两点边值问题的一种解析算法，并和现今数值打靶法获得的结果进行对比，分析比较2种方法的入轨精度和计算时间，得出各自方法的优缺点，为后续火箭制导方法选择和制导方案设计提供参考。

1 两点边值问题

考虑一般的飞行器轨迹优化问题，从已知的位置速度到目标圆轨道的时间最优问题，根据极大值原理可以将问题转化为两点边值问题。简化的力学模型，飞行器受到推力F（大小和方向）和重力的作用。坐标系选择为Oxyz惯性坐标系。其动力学方程为

式中 x，y为位置坐标；u，v分别为x方向和y方向的速度；gx，gy分别为重力加速度在x方向和y方向的分量；χ为控制方向角。

飞行器在惯性系下运动如图1所示。

图1 Oxyz惯性坐标系

根据系统动力学方程得到系统的Hamilton函数为

式中 λ1，λ2分别为速度变量u，v的Lagrange乘子（也称协态变量）；λ3，λ4分别为位置变量x，y的Lagrange乘子。

由最优控制理论可以得到系统的协态方程为

根据极大值原理，得到最优控制条件：

从式（4）可以看出最优控制方向矢量仅仅取决于和速度变量u，v对应的协态变量λ1，λ2。协态变量归一化条件：

横截条件：

对于本文考虑的火箭上升进入目标圆轨道情况，末端条件仅需通过入轨位置约束F1、速度约束F2和方向约束F33个约束即可进行描述：

式中 Rco，Vco分别为末端位置和速度大小；F3为入轨方向的切向入轨约束。式（1）～（7）构成了标准的轨迹优化两点边值问题，针对此类问题，求解的数值方法有很多，如间接法求解两点边值问题的打靶法[6]。

2 解析法求解最优轨迹

求解最优飞行轨迹的两点边值问题获得控制量，可以转化为求解和速度相对应的协态变量，根据协态变量归一化条件平方和等于1，不妨假设协态变量（控制量）的初值为

式中 x0为控制方向角的初值。记飞行器整个飞行时间为tco，起始时刻为t0，则剩余飞行时间Δt可以记为

解析法的关键是将当前时刻的协态变量（控制量）写成过程状态量的函数，即Lagrange乘子的解可以写成明确的解析表达式，通过状态变量和飞行器任务参数（例如，质量m、推力F、末端位置Rco和末端速度Vco等）显示表示出来，就是所希望的自适应导引函数：

式（10）表示了协态变量和状态量参数间的函数关系，再根据协态变量和控制量的关系式（4），就可以得到控制变量X和状态量参数间的函数关系式（11）。为了获得式（10）的具体解析表达式，需要构造4个未知Lagrange乘子的4个代数方程，乘子在方程中明确地出现，并且其系数根据方程中的变量确定，其中标准化方程式（5）和横截条件方程式（6）已经可以确定，另外2个可根据式（7）获得。将式（7）关于剩余飞行时间Δt= tco- t0进行Taylor级数展开，这3个末端条件根据初始条件展开为

解析法作为精确解的近似，级数展开仅取有限项数p，式（12）中Δt的系数在每一个当前t0时刻得到，协态变量λ1～λ4将在级数展开过程中出现，但是当Δt的值未知时，仍无求得协态变量的解析表达式。为了确定Δt，将式（12）改写成有限级数或者多项式形式：

进一步如果有：

式中 未知的Δt可通过已知的系数Bn表示（Bn中包含状态量和协态量信息）；q为反转求过程中的展开项数，从式（13）、式（14）过程中获得Bn解析表达式的方法叫做级数反转。这种方法程序技巧已在文献[7]中给出。式（13）可以反转求得的表达式（14），再代回到式（12）即可求得协态变量的解析方程式。

针对3个末端条件方程，通过大量计算[4]，采用基于末端速度约束F2来求解剩余飞行时间精度优于其他2个方程。由F2得到剩余飞行时间后，代回另外2个末端约束F1和F3，综合前面的式（5）、式（6），实现了求解4个协态变量的4个代数方程：

式中 i+j+k+I=0,1,2，...，表示协态变量在不同项级数展开过程中的不同阶次。Rijkd=f(x,y,u,v,m,F,m,Rco,Vco)，是通过方程式（12）求得的协态变量解析方程式。随着剩余飞行时间展开阶数和末端约束展开阶数的增加，方程的项数也将迅速增加，综合考虑入轨精度要求和计算效率，本文考虑了展开1～4阶的情况进行计算仿真。

前面已经得到了关于4个协态变量代数方程组，该方程组是非线性方程组，可通过Newton-Raphson等方法求解。考虑到协态变量经过归一化处理都小于1，且飞行器在真空飞行段的速度θ倾角一般介于±30°之间，即x=90°-θ介于60～120°之间，且接近90°，因此λ2=cosx接近于0。此时非线性方程组的求解可以采用连续替代法[8]，经过几步迭代求解可以收敛。由式（15）中的第3式和第4式得到：

3 仿真分析

采用文献[4]中的一组飞行参数，飞行器初始质量150 t，推力889 kN，发动机比冲420 s，初始径向位置6555.9 km，水平飞行速度6.78 km/s。末端目标圆轨道，轨道半径为6565.7 km。分别采用解析法和数值法开展仿真实验。解析法每隔1 s进行一次控制更新，至剩余飞行时间Δt小于10 s时控制量保持不变，直至剩余时间小于1×10-4s时停止运算。数值法采用打靶法求解，积分步长设置为1 s，收敛精度设为1×10-6m。解析法分别采用剩余时间和末端条件约束的不同展开阶数进行计算，比较不同组合的入轨精度以及单步平均计算时间，并和数值法进行比较，见表1。

表1 解析法和数值法精度分析和计算时间对比

从表1中可以看出，解析法随着剩余时间和末端条件约束展开阶数的增加，入轨精度提高，其中入轨速度误差量级维持在10-3m/s，入轨位置误差精度明显提高。同时随着解析法展开阶数的增加，单步计算时间逐渐增加，这是由于随阶数增加解析模型复杂程度迅速增加。数值法计算采用VC++6.0语言编码，解析法采用的是Matlab R2012语言编码且计算时公式没有充分化简，计算时间仍有一定提升空间。综合来看解析法的计算精度虽然没有数值法的精度高，但是计算时间明显快于数值法，解析法降低了一部分入轨精度带来计算时间的节省，为未来实时动态制导律设计提供可能。下面给出2种方法速度倾角随时间变化曲线，速度倾角θ=90°-x。

图2中最下方的曲线（粗实线）为数值法速度倾角随时间变化曲线，末端有一段倾角保持常值的曲线是各个阶次解析法速度倾角控制曲线。从图2中可以看出，解析法能够很好地逼近数值结果，并且随着解析法阶数增加，曲线逼近程度更高。此外解析法初始速度倾角会有小幅震荡，这是由于解析法无法猜测初值，本文给出固定的初始速度倾角50°，但只需几步运算更新就可以实现倾角平稳变化。

图2 解析法和数值法速度倾角变化曲线

4 结束语

本文针对实时在线规划制导控制需求，用解析法求解了轨迹优化问题中的两点边值问题，给出了实时在线更新控制律的解析算法，并和数值打靶法对比分析了入轨精度和计算时间的差异。解析法的计算精度虽然没有数值法的精度高，但是计算时间明显快于数值法，随着解析阶数的增加，解析法入轨精度提高。

[1] Lu P, Pan B. Highly constrained optimal launch ascent guidance[J]. Journal of Guidance Control & Dynamics, 2012, 33(2): 404-414.

[2] Johnson E N, Calise A J, Curry M D. Adaptive guidance and control for autonomous hypersonic vehicles[J]. Journal of Guidance Control & Dynamics, 2006, 29(3): 725-736.

[3] Lu P. Predictor-corrector entry guidance for low lifting vehicles[C]. Hilton Head: AIAA Guidance, Navigation and Control Conference, 2007.

[4] Gilchrist C A, Morrow L. An analytical approach to solution of two point boundary condition problems in optimal guidance[R]. Summary report, TM-292-6-038, 1966.

[5] Hu S S, Thompson M L. A direct and analytical solution for space flight guidance function[C]. New York: AIAA Paper, Presented to 3rd Aerospace Science Meeting of AIAA, 1966.

[6] Haberkorn, T P, Gergaud M J. Low-thrust minimum-fuel orbital transfer: a homotopic approach[J]. Journal of Guidance Control & Dynamics, 2004, 27(6): 1046-1060.

[7] Thorne J D. Series reversion/inversion of Lambert’s time function. Master Degree[R], Air Force Institute of Technology Air University, Captain USAF, 1989.

[8] Thompson M L, Gray W J. A method for the solution of systems of nonlinear algebraic equations having non-numerical coefficients[R]. NSL Research Analysis Section Tech Memo No. 141, 1965.

Online Real-time Guidance Law Using an Analytical Optimal-trajectory Method

Zhang Zhi-guo, Yu Meng-lun
(Beijing Institute of Astronautical System Engineering, Beijing, 100076）

The numerical guidance methods of optimal trajectory are high accuracy, but time-consuming and algorithm-unstable. With the development of computer technology, time-consuming weakness can be alleviated, but still unable to meet the real-time dynamic guidance and control requirements. Analytical methods have been in-depth study before the instance of numerical methods. If combines the advantages of these two methods, it will be able to achieve the real-time control with high accuracy. Firstly, an analytical method is solved with the two-point boundary-value problem of optimal trajectory, then analyze the differences of in-orbit precision and computing time between the analytical and numerical methods. Finally, the advantages and disadvantages of each method are given.

Online real-time guidance; Analytical method; Trajectory optimization

V448.1

A

1004-7182(2016)05-0058-04

10.7654/j.issn.1004-7182.20160513

2015-11-19；

2016-07-13

张志国（1991-），男，博士研究生，主要研究方向为飞行器动力学，制导与控制