基于解析法最優軌跡的在線實時制導律研究

張志國，余夢倫

（北京宇航系統工程研究所，北京，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-），男，博士研究生，主要研究方向為飛行器動力學，制導與控制