張志國，余夢(mèng)倫

（北京宇航系統(tǒng)工程研究所，北京，100076）

論文與報(bào)告

基于解析法最優(yōu)軌跡的在線實(shí)時(shí)制導(dǎo)律研究

張志國，余夢(mèng)倫

（北京宇航系統(tǒng)工程研究所，北京，100076）

新型制導(dǎo)系統(tǒng)單純采用軌跡優(yōu)化的數(shù)值算法的優(yōu)點(diǎn)是精度高，但耗時(shí)長，算法穩(wěn)定性有待提高，隨著計(jì)算機(jī)技術(shù)的發(fā)展，耗時(shí)長的缺點(diǎn)能得到緩解，但仍然無法滿足實(shí)時(shí)動(dòng)態(tài)制導(dǎo)控制需求。數(shù)值法以前，解析法曾得到深入的研究，如果能結(jié)合二者的優(yōu)點(diǎn)，將能更好地達(dá)到提升飛行器實(shí)時(shí)高精度制導(dǎo)控制的效果。用解析法求解了軌跡優(yōu)化問題中的兩點(diǎn)邊值問題，并和數(shù)值打靶法對(duì)比分析了入軌精度和計(jì)算時(shí)間的差異，得出各自方法的優(yōu)缺點(diǎn)。

實(shí)時(shí)在線制導(dǎo)；解析法；軌跡優(yōu)化

0 引 言

飛行器軌跡優(yōu)化和制導(dǎo)控制一直是飛行器總體設(shè)計(jì)的主要研究方向。軌跡優(yōu)化研究中常用到的方法有間接法、直接法和混合法，目前都是依賴數(shù)值算法[1]。但隨著飛行器性能的日益提升，更高任務(wù)要求的不斷提出，對(duì)于飛行器軌跡設(shè)計(jì)和制導(dǎo)性能提出了更高的要求[1]。例如，飛行軌跡實(shí)時(shí)在線規(guī)劃問題[2]，預(yù)測(cè)制導(dǎo)問題[3]等，在實(shí)際飛行過程中需要迅速給出最優(yōu)控制策略，耗時(shí)長的數(shù)值算法是無法滿足實(shí)時(shí)動(dòng)態(tài)制導(dǎo)控制需求的。雖然隨著計(jì)算機(jī)技術(shù)的發(fā)展，數(shù)值算法耗時(shí)長的缺點(diǎn)能得到緩解，但仍然無法滿足實(shí)時(shí)動(dòng)態(tài)制導(dǎo)控制的需求，算法穩(wěn)定性和魯棒性也是需要注意的問題。解析法優(yōu)點(diǎn)是計(jì)算時(shí)間短、穩(wěn)定性高，但是精度高低取決于簡化程度、展開階數(shù)等因素。若將數(shù)值法和解析法相結(jié)合考慮，將能更好地達(dá)到提升飛行器實(shí)時(shí)高精度控制的效果。

Hu Steve S[4]在1966年給出了一套解析法求解最優(yōu)飛行導(dǎo)引方程的方法，但是并沒有給出具體的求解過程和精度分析。Gilchrist C和Morrow L[5]隨后對(duì)制導(dǎo)理論的解析進(jìn)行了總結(jié)，對(duì)于解析法求解兩點(diǎn)邊值問題給出了詳細(xì)的論述，重點(diǎn)分析了解析法展開到不同階數(shù)時(shí)協(xié)態(tài)變量相對(duì)標(biāo)準(zhǔn)值的偏差。解析求解方法后續(xù)在多項(xiàng)式制導(dǎo)和迭代制導(dǎo)的研究中被引用。隨著計(jì)算機(jī)技術(shù)的發(fā)展，數(shù)值方法獲得了迅速發(fā)展，如目前求解兩點(diǎn)邊值問題的打靶法。

本文中提出了求解兩點(diǎn)邊值問題的一種解析算法，并和現(xiàn)今數(shù)值打靶法獲得的結(jié)果進(jìn)行對(duì)比，分析比較2種方法的入軌精度和計(jì)算時(shí)間，得出各自方法的優(yōu)缺點(diǎn)，為后續(xù)火箭制導(dǎo)方法選擇和制導(dǎo)方案設(shè)計(jì)提供參考。

1 兩點(diǎn)邊值問題

考慮一般的飛行器軌跡優(yōu)化問題，從已知的位置速度到目標(biāo)圓軌道的時(shí)間最優(yōu)問題，根據(jù)極大值原理可以將問題轉(zhuǎn)化為兩點(diǎn)邊值問題。簡化的力學(xué)模型，飛行器受到推力F（大小和方向）和重力的作用。坐標(biāo)系選擇為Oxyz慣性坐標(biāo)系。其動(dòng)力學(xué)方程為

式中 x，y為位置坐標(biāo)；u，v分別為x方向和y方向的速度；gx，gy分別為重力加速度在x方向和y方向的分量；χ為控制方向角。

飛行器在慣性系下運(yùn)動(dòng)如圖1所示。

圖1 Oxyz慣性坐標(biāo)系

根據(jù)系統(tǒng)動(dòng)力學(xué)方程得到系統(tǒng)的Hamilton函數(shù)為

式中 λ1，λ2分別為速度變量u，v的Lagrange乘子（也稱協(xié)態(tài)變量）；λ3，λ4分別為位置變量x，y的Lagrange乘子。

由最優(yōu)控制理論可以得到系統(tǒng)的協(xié)態(tài)方程為

根據(jù)極大值原理，得到最優(yōu)控制條件：

從式（4）可以看出最優(yōu)控制方向矢量僅僅取決于和速度變量u，v對(duì)應(yīng)的協(xié)態(tài)變量λ1，λ2。協(xié)態(tài)變量歸一化條件：

橫截條件：

對(duì)于本文考慮的火箭上升進(jìn)入目標(biāo)圓軌道情況，末端條件僅需通過入軌位置約束F1、速度約束F2和方向約束F33個(gè)約束即可進(jìn)行描述：

式中 Rco，Vco分別為末端位置和速度大小；F3為入軌方向的切向入軌約束。式（1）～（7）構(gòu)成了標(biāo)準(zhǔn)的軌跡優(yōu)化兩點(diǎn)邊值問題，針對(duì)此類問題，求解的數(shù)值方法有很多，如間接法求解兩點(diǎn)邊值問題的打靶法[6]。

2 解析法求解最優(yōu)軌跡

求解最優(yōu)飛行軌跡的兩點(diǎn)邊值問題獲得控制量，可以轉(zhuǎn)化為求解和速度相對(duì)應(yīng)的協(xié)態(tài)變量，根據(jù)協(xié)態(tài)變量歸一化條件平方和等于1，不妨假設(shè)協(xié)態(tài)變量（控制量）的初值為

式中 x0為控制方向角的初值。記飛行器整個(gè)飛行時(shí)間為tco，起始時(shí)刻為t0，則剩余飛行時(shí)間Δt可以記為

解析法的關(guān)鍵是將當(dāng)前時(shí)刻的協(xié)態(tài)變量（控制量）寫成過程狀態(tài)量的函數(shù)，即Lagrange乘子的解可以寫成明確的解析表達(dá)式，通過狀態(tài)變量和飛行器任務(wù)參數(shù)（例如，質(zhì)量m、推力F、末端位置Rco和末端速度Vco等）顯示表示出來，就是所希望的自適應(yīng)導(dǎo)引函數(shù)：

式（10）表示了協(xié)態(tài)變量和狀態(tài)量參數(shù)間的函數(shù)關(guān)系，再根據(jù)協(xié)態(tài)變量和控制量的關(guān)系式（4），就可以得到控制變量X和狀態(tài)量參數(shù)間的函數(shù)關(guān)系式（11）。為了獲得式（10）的具體解析表達(dá)式，需要構(gòu)造4個(gè)未知Lagrange乘子的4個(gè)代數(shù)方程，乘子在方程中明確地出現(xiàn)，并且其系數(shù)根據(jù)方程中的變量確定，其中標(biāo)準(zhǔn)化方程式（5）和橫截條件方程式（6）已經(jīng)可以確定，另外2個(gè)可根據(jù)式（7）獲得。將式（7）關(guān)于剩余飛行時(shí)間Δt= tco- t0進(jìn)行Taylor級(jí)數(shù)展開，這3個(gè)末端條件根據(jù)初始條件展開為

解析法作為精確解的近似，級(jí)數(shù)展開僅取有限項(xiàng)數(shù)p，式（12）中Δt的系數(shù)在每一個(gè)當(dāng)前t0時(shí)刻得到，協(xié)態(tài)變量λ1～λ4將在級(jí)數(shù)展開過程中出現(xiàn)，但是當(dāng)Δt的值未知時(shí)，仍無求得協(xié)態(tài)變量的解析表達(dá)式。為了確定Δt，將式（12）改寫成有限級(jí)數(shù)或者多項(xiàng)式形式：

進(jìn)一步如果有：

式中 未知的Δt可通過已知的系數(shù)Bn表示（Bn中包含狀態(tài)量和協(xié)態(tài)量信息）；q為反轉(zhuǎn)求過程中的展開項(xiàng)數(shù)，從式（13）、式（14）過程中獲得Bn解析表達(dá)式的方法叫做級(jí)數(shù)反轉(zhuǎn)。這種方法程序技巧已在文獻(xiàn)[7]中給出。式（13）可以反轉(zhuǎn)求得的表達(dá)式（14），再代回到式（12）即可求得協(xié)態(tài)變量的解析方程式。

針對(duì)3個(gè)末端條件方程，通過大量計(jì)算[4]，采用基于末端速度約束F2來求解剩余飛行時(shí)間精度優(yōu)于其他2個(gè)方程。由F2得到剩余飛行時(shí)間后，代回另外2個(gè)末端約束F1和F3，綜合前面的式（5）、式（6），實(shí)現(xiàn)了求解4個(gè)協(xié)態(tài)變量的4個(gè)代數(shù)方程：

式中 i+j+k+I=0,1,2，...，表示協(xié)態(tài)變量在不同項(xiàng)級(jí)數(shù)展開過程中的不同階次。Rijkd=f(x,y,u,v,m,F,m,Rco,Vco)，是通過方程式（12）求得的協(xié)態(tài)變量解析方程式。隨著剩余飛行時(shí)間展開階數(shù)和末端約束展開階數(shù)的增加，方程的項(xiàng)數(shù)也將迅速增加，綜合考慮入軌精度要求和計(jì)算效率，本文考慮了展開1～4階的情況進(jìn)行計(jì)算仿真。

前面已經(jīng)得到了關(guān)于4個(gè)協(xié)態(tài)變量代數(shù)方程組，該方程組是非線性方程組，可通過Newton-Raphson等方法求解?？紤]到協(xié)態(tài)變量經(jīng)過歸一化處理都小于1，且飛行器在真空飛行段的速度θ傾角一般介于±30°之間，即x=90°-θ介于60～120°之間，且接近90°，因此λ2=cosx接近于0。此時(shí)非線性方程組的求解可以采用連續(xù)替代法[8]，經(jīng)過幾步迭代求解可以收斂。由式（15）中的第3式和第4式得到：

3 仿真分析

采用文獻(xiàn)[4]中的一組飛行參數(shù)，飛行器初始質(zhì)量150 t，推力889 kN，發(fā)動(dòng)機(jī)比沖420 s，初始徑向位置6555.9 km，水平飛行速度6.78 km/s。末端目標(biāo)圓軌道，軌道半徑為6565.7 km。分別采用解析法和數(shù)值法開展仿真實(shí)驗(yàn)。解析法每隔1 s進(jìn)行一次控制更新，至剩余飛行時(shí)間Δt小于10 s時(shí)控制量保持不變，直至剩余時(shí)間小于1×10-4s時(shí)停止運(yùn)算。數(shù)值法采用打靶法求解，積分步長設(shè)置為1 s，收斂精度設(shè)為1×10-6m。解析法分別采用剩余時(shí)間和末端條件約束的不同展開階數(shù)進(jìn)行計(jì)算，比較不同組合的入軌精度以及單步平均計(jì)算時(shí)間，并和數(shù)值法進(jìn)行比較，見表1。

表1 解析法和數(shù)值法精度分析和計(jì)算時(shí)間對(duì)比

從表1中可以看出，解析法隨著剩余時(shí)間和末端條件約束展開階數(shù)的增加，入軌精度提高，其中入軌速度誤差量級(jí)維持在10-3m/s，入軌位置誤差精度明顯提高。同時(shí)隨著解析法展開階數(shù)的增加，單步計(jì)算時(shí)間逐漸增加，這是由于隨階數(shù)增加解析模型復(fù)雜程度迅速增加。數(shù)值法計(jì)算采用VC++6.0語言編碼，解析法采用的是Matlab R2012語言編碼且計(jì)算時(shí)公式?jīng)]有充分化簡，計(jì)算時(shí)間仍有一定提升空間。綜合來看解析法的計(jì)算精度雖然沒有數(shù)值法的精度高，但是計(jì)算時(shí)間明顯快于數(shù)值法，解析法降低了一部分入軌精度帶來計(jì)算時(shí)間的節(jié)省，為未來實(shí)時(shí)動(dòng)態(tài)制導(dǎo)律設(shè)計(jì)提供可能。下面給出2種方法速度傾角隨時(shí)間變化曲線，速度傾角θ=90°-x。

圖2中最下方的曲線（粗實(shí)線）為數(shù)值法速度傾角隨時(shí)間變化曲線，末端有一段傾角保持常值的曲線是各個(gè)階次解析法速度傾角控制曲線。從圖2中可以看出，解析法能夠很好地逼近數(shù)值結(jié)果，并且隨著解析法階數(shù)增加，曲線逼近程度更高。此外解析法初始速度傾角會(huì)有小幅震蕩，這是由于解析法無法猜測(cè)初值，本文給出固定的初始速度傾角50°，但只需幾步運(yùn)算更新就可以實(shí)現(xiàn)傾角平穩(wěn)變化。

圖2 解析法和數(shù)值法速度傾角變化曲線

4 結(jié)束語

本文針對(duì)實(shí)時(shí)在線規(guī)劃制導(dǎo)控制需求，用解析法求解了軌跡優(yōu)化問題中的兩點(diǎn)邊值問題，給出了實(shí)時(shí)在線更新控制律的解析算法，并和數(shù)值打靶法對(duì)比分析了入軌精度和計(jì)算時(shí)間的差異。解析法的計(jì)算精度雖然沒有數(shù)值法的精度高，但是計(jì)算時(shí)間明顯快于數(shù)值法，隨著解析階數(shù)的增加，解析法入軌精度提高。

[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-），男，博士研究生，主要研究方向?yàn)轱w行器動(dòng)力學(xué)，制導(dǎo)與控制