采用聯立法求解大姿態終端約束的上升段軌跡優化

邱 豐，宋征宇

(北京航天自動控制研究所，北京100854)

采用聯立法求解大姿態終端約束的上升段軌跡優化

邱 豐，宋征宇

(北京航天自動控制研究所，北京100854)

針對恒定推力和未配置輔助動力系統的運載火箭，提出了一種基于聯立框架的直接法(“聯立法”)來求解帶大姿態終端約束的動態在線軌跡規劃問題，解決了傳統迭代制導方法中飛行軌跡線性規劃的不足，并可根據姿態跟蹤精度的需求來控制程序角速率，從而同時滿足各種約束條件。該聯立法結合了直接配點法和偽譜法的特點，采用有限元正交配置法將狀態變量和控制變量進行離散化;結果表明與偽譜法相比，計算效率得到了大幅提升。在此基礎上，提出了通過合理選擇初值、調整入軌精度偏差、自適應有限元分段等提升計算速度的措施，進一步提高了這一算法的實時性，為未來解決更復雜多樣的制導任務提供了基礎。

運載火箭;姿態終端約束;聯立法;軌跡規劃

0 引 言

近年來，自主控制的需求在航天應用中愈發增強。上升段飛行中，對于目標變更和推力下降等異常情況，自主規劃的適應性較跟蹤標準軌跡的攝動制導要高很多。自載人航天交會對接任務起，我國運載火箭首次采用了迭代制導方法［1］。這是一種在線動態自主軌跡規劃的方法，但由于飛行軌跡的控制是通過調整姿態角來實現的，使得入軌姿態的散布很大。

解決入軌姿態可以從總體方案著手，例如在火箭的末級增加輔助動力系統(也稱調姿系統)，或者依靠有效載荷自身進行姿態調整，但這些方案均額外增加了系統資源和復雜性。為此，文獻［2］提出了一種預測修正的迭代制導改進方法。但是，對于一些較為特殊的入軌姿態角，現有迭代制導中線性規劃的前提條件不再存在，導致入軌誤差很大。對于此類問題，即使傳統攝動制導也難以滿足精度要求(這將在后文進行分析)。

本文提出了一種基于聯立框架的直接法(下文簡稱“聯立法”)，將入軌時間最短作為性能指標，通過控制推力方向變化使得火箭以要求姿態進入預定軌道。該方法采用統一化的處理流程，非常適合計算機自動處理，在自主制導領域有較好的應用前景。

1 帶終端姿態約束的最優入軌問題

對于推力大小不可調節的運載火箭上升段實時軌跡優化問題，其本質上是帶有狀態約束和控制約束的最優控制問題，間接法和直接法是最常用的兩種方法，但間接法對于末端有姿態角約束且入軌時間浮動的情況無法完成求解［3］。而直接法將連續空間最優控制問題轉化為非線性規劃(NLP)問題，利用數值方法求解。文獻［4］利用偽譜法對固體運載火箭上升段軌跡優化進行研究，但計算耗時無法滿足實時控制要求;文獻［5］改善了計算效率，但收斂性還需進一步驗證。目前研究成果表明［6－8］，直接法的實時性和收斂性問題并沒有完全解決。

本文研究更為一般性的帶終端姿態約束的精確入軌問題，采用有限元正交配置法將狀態變量和控制變量進行離散化，實現對原命題的分段低階逼近處理，使原問題轉化為NLP問題，同時提出了提高計算效率的措施。

1.1 火箭質心運動模型

在不考慮氣動力影響的條件下，建立其相對于發射慣性坐標系(以下簡稱發慣系)的質心運動方程如下［9］:

式中:x，y，z為火箭當前點到發射點矢徑在發慣系上的投影，Vx，Vy，Vz為火箭速度矢量在發慣系上的投影，gx，gy，gz為引力加速度在發慣系上的投影，m0為飛行器初始質量，m·為秒耗量，P為推力，φ為俯仰角，ψ為偏航角。對于定推力發動機，秒耗量及推力設為定值。

由于狀態變量之間數量級相差較大，為了便于數值求解，需要對運動方程進行單位歸一化處理:

1.2 需要滿足的約束條件

有效載荷的運行軌道由五個軌道根數決定(如果考慮其某一時刻在軌道上的具體位置，則需要六個軌道根數)［10］。只有所處位置在軌道上的一點，并且速度大小及方向恰好是軌道上該點所需速度的情況下，才能實現對五個軌道根數的精確控制。

帶姿態約束的最優入軌問題中，不僅要控制火箭精確進入預定軌道，而且需滿足一定的入軌姿態。此外，由于火箭燃料有限，為了確保得到的解是可行的，還需要加入最大入軌時間的約束。

綜上所述，帶姿態約束的火箭最優入軌問題中要滿足的約束如表1所示。

表1 命題約束列表Table 1 Constraints list of the proposition

2 基于有限元正交配置的聯立法

2.1 有限元正交配置法

假定整個問題的時間區間為［t0，tf］，將上升段入軌問題描述為一般動態優化命題形式［11］:

式中:Φ(x(tf))是目標函數，x(t)∈Rnx為狀態變量組成的nx維向量;u(t)∈Rnu為控制變量組成的nu維向量;p∈Rnp是模型參數;f為微分方程組的右邊約束;gE為代數方程組;微分代數方程組中的狀態變量和控制變量的邊界約束均歸結到gI中; x0∈Rnx為狀態變量的初始狀態，N1(x(tf))為終端目標集約束。

傳統偽譜法采用全局插值多項式擬合，為了保證逼近精度要求多項式具有較高階次，離散化后NLP問題規模很大，導致了計算難度和求解時間的增加。為了在保證計算精度的同時降低離散化后命題的復雜程度，可將整個時域劃分為N個區間，每個區間稱為有限元，在每段有限元［ti－1，ti］上，采用正交多項式對狀態變量x及控制變量u進行逼近，采用拉格朗日插值多項式可使變量在配置點上的值正好等于其系數。這種方法稱之為有限元正交配置法［12］。

此外，在有限元端點上需通過連接方程來保證狀態變量的連續性;同時應滿足初始和終端約束條件。結合以上離散化后的條件，可以得到有限元正交配置法處理后的NLP問題。

2.2 聯立法對問題的重構

以本問題為例，將整個時域［t0，tf］分為N個有限元，每個有限元內用K階拉格朗日插值多項式對變量和約束進行逼近。采用Gauss-Radau點［13］作為配置點進行配置。

目標函數為入軌時間:

狀態變量為火箭的速度、位置:

控制變量為火箭的姿態程序角:

在推力大小不變的情況下，定點入軌是無法保證的。因此在求解過程中將真近點角f作為代數變量進行處理，通過五個軌道根數及f，可以確定入軌點速度、位置，從而形成六個等式約束方程。

(1)過程約束

火箭質心運動的狀態變量微分方程表示如下:

式中，i∈［0，N)，j∈［0，K)。

狀態變量的連續性方程表示為(以x為例):

實際飛行中，還應加入姿態角及姿態角速率的不等式約束(以俯仰角為例，ω為角速率上限):

(2)終端約束

狀態變量的初始約束和終端約束(以x為例):

目標入軌時間最大值的約束:

綜合式(7)～(11)可得帶終端姿態約束最優入軌問題離散化后的大規模非線性規劃命題。

由式(9)和(11)可以看出該非線性規劃命題中含有不等式約束，直接對其求解非常困難，且約束個數隨著有限元個數N和插值階數K增加而成倍增加，故采用內點算法［12］對命題進行求解。

至此，聯立法的計算過程可歸納為圖1所示，分為設計輸入、離散處理和實時規劃三部分。通過利用基于內點算法的NLP求解器迭代求解，將求得變量代入離散化后的命題得到約束殘差再返回求解器。通過迭代計算使殘差漸漸逼近零，求得的變量也越來越逼近最優解。在該計算過程中，終端姿態角約束被統一處理，并不區分姿態角的大小，體現了該方法的通用性。

2.3 聯立法的性能分析

在［t0，tf］內將整個時域分成N段，聯立法通過有限元將狀態變量x(t)分別定義到多個時段內。在有限元［ti－1，ti］上，由式(7)得到的近似最優解xk+1(t)和真實解x(t)之間的誤差邊界可以通過下式來描述［13］:

式中:C為只與階次K有關的常數，hi為第i個有限元的長度，為真實解K+1次導數在區間［ti－1，ti］上的無窮范數。因此，當有界時，可以得到:

也就是說當有限元的劃分無限加密的時候，聯立法得到的近似解無限逼近真實的解，故可以通過控制離散化精度來間接保證解的最優性。但工程實際應用中不能為了少量精度的提高而增加大量的有限元，在線規劃時可根據實際用例選取合適的有限元，以達到求解精度和運算時間的平衡。

3 仿真校驗

3.1 仿真場景

表2～表3以某火箭上升段飛行為例列出了仿真場景的各項參數。暫不考慮大氣影響，即從火箭二級飛行段開始。表4為終端姿態要求，其中終端俯仰角為－110°，正是這一特殊約束，導致傳統算法無法同時兼顧入軌精度和終端姿態的要求。

表2 初始參數Table 2 Initial parameters

表3 目標軌道參數Table 3 Target orbit parameters

表4 姿態角約束Table 4 Attitude constraints

平均取30段有限元和3階Radau配點進行離散化，初值取法如下:

1)各配點速度、位置初值均設為初始點值;

2)真近點角變量的初值取為0°;

3)目標函數(總飛行時間)的初值取標準彈道入軌時間376.4 s;

4)由于俯仰角變化大，其初值采用等差選取的方法，而偏航角變化不大，初值取初始狀態的值。

3.2 仿真結果與分析

本次數值實驗是基于 IntelRCoreTMi5－3470 CPU處理器、3.2 GHz主頻的計算機進行的，求解器為IPOPT(3.6.1版)［14］。由于決定運算速度的主要因素是處理器主頻，所以不同硬件平臺下的計算耗時可以根據CPU的主頻按比例近似折算。

3.2.1 與標準彈道進行對比分析

圖2給出了標準彈道設計與聯立法計算結果的飛行程序角對比:對于終端俯仰角－110°的約束條件，相當于火箭向下飛行，這對運載能力是有損失的，因此標準彈道設計時將姿態調整段壓縮在飛行的末斷，時間要盡可能短，俯仰角速率最大達到9.35°/s，如圖3所示。六自由度仿真結果表明，在該角速率條件下，姿態跟蹤誤差達到了2°，這反過來影響了入軌精度。這就是引言中所提采用攝動制導跟蹤標準軌跡飛行也難以滿足入軌精度的原因。

為此在聯立法計算中將角速率約束定為4°/s，終端姿態和入軌精度均能滿足要求，入軌時間為377.5 s，比標準彈道多 1.1 s(多消耗推進劑約118.6 kg，在安全余量范圍內)，如圖3所示。且軌跡可以實時動態調整，及時修正之前的誤差，即誤差不具有積累的效果。

從圖4速度和位置曲線可以看出，該飛行軌跡盡管前后段均近似線性，但整體難以用線性規劃來統一描述;而大角速率轉彎飛行的起點是浮動的，無法提前獲知，使得分段線性化處理也不可能，這是導致現有迭代制導無法滿足該算例入軌精度要求的主要原因。

3.2.2 與偽譜法求解結果進行對比分析

在同一運算平臺下，分別用全局30和300個Radau配置點的偽譜法對命題進行求解，得到的結果與30段3階Radau配置點的聯立法求解情況對比，如表5和圖5所示。

表5 聯立法與偽譜法結果對比Table 5 Comparison of the simultaneous method and the pseudo-spectral method

由圖5可見，求解時間近似的情況下，聯立法求解精度高于使用全局近似的偽譜法;而求解精度相似時，偽譜法的耗時遠遠超過聯立法。因此聯立法在求解效率上優于全局插值近似的偽譜法。

3.3 提升算法實時性的幾種手段

3.3.1 選取合理的初值

由于每個控制周期均迭代計算，即使在一個周期內受到干擾等影響，實際飛行軌跡也不會偏離過大。因此可采用當前點速度位置作為起始點初值，上一次規劃的相關參數作為其他狀態及控制變量的初值，以縮短優化時間，這對實時飛行的迭代計算意義十分重大。

針對本算例，將所有變量的初值均設置為已有優化結果得到的最優值附近，假設偏差為3%，則各個變量的初值取為:

計算結果對比如表6所示。

表6 不同初值結果對比Table 6 Comparison of different initial value

由表6可以看出，該方法可以大大縮短運算時間;且兩次計算之間的時間間隔越短，初值與最優解之間的偏差也越小。

3.3.2 選擇合適的入軌點精度偏差

聯立法是將軌道根數轉化得到的速度、位置作為終端等式約束條件，即要求入軌點的速度、位置與理論值的差值為0。而實際應用中，不同任務對不同軌道根數的精度要求也不同;并且在合理的偏差范圍內，即差值不為0也是可以接受的。因此可以對不同軌道根數區別處理以提升計算效率。

以本算例為例，在保持整體容差設置為1×10－8不變的情況下，將變量值與要求值的差除以10000，則相當于把六個約束放寬到1×10－4，結果對比如表7和表8所示。

表7 放寬誤差約束后計算效率對比Table 7 Comparison of computational efficiency after the relaxation of error constraints

從表7中看出，放寬誤差約束可以縮短優化時間，但對入軌精度的影響在可接受的范圍之內，如表8所示。因此，未來可以根據不同算例的特點對約束進行區別調整，在不影響精度的情況下進一步提高實時性。

表8 放寬誤差約束后入軌精度對比Table 8 Comparison of orbit injection accuracy after the relaxation of error constraints

3.3.3 自適應有限元分段

減少有限元的分段能有效減小離散后NLP命題的規模，但會降低解的精度和最優性。通過上述飛行曲線的觀察可以發現，控制變量的大幅變化基本集中在初始和終端位置，中間很長的時間變化平緩。根據這一特點，在變量快速變化的區間內加密有限元劃分，而在變化平緩的區域減少有限元分段，這種方法可稱作自適應有限元分段法。

為此，將整個飛行過程分為三個時間段:［0，20 s］、［20 s，350 s］和［350，tf］，三段內各取5個有限元，每個有限元內同樣取3階Radau點進行配置，結果如表9所示。

表9 不同有限元分段結果對比Table 9 Comparison of different finite element configurations

在實際使用中，可以利用起飛前離線規劃的軌跡，通過姿態角變化率的求解，調整有限元分段;在飛行過程中，也可以利用這一規律自適應地動態調整分段。

4 結 論

本文針對傳統迭代制導方法無法保證終端大姿態約束的問題，利用有限元正交配置進行低階逼近和聯立求解，與偽譜法等直接法相比計算效率得到大幅提高。并以某火箭上升段最優入軌任務為例，對算法的效果進行了仿真驗證。結果表明，入軌精度和終端大姿態約束均能滿足，求解時間也縮短到百毫秒級，為實時在線規劃打下了基礎。

［1］ 宋征宇.從準確、精確到精益求精——載人航天推動運載火箭制導方法的發展［J］.航天控制，2013，31(1):4－10.［Song Zheng-yu.From accurate，precise to perfect-manned space promotes the development of guidance method on launch vehicle［J］.Aerospace Control，2013，31(1):4－10.］

［2］ 呂新廣，宋征宇，鞏慶海.運載火箭軌跡預測制導方法研究［J］.載人航天，2013，19(6):24－28.［Lv Xin-guang;Song Zheng-yu;Gong Qing-hai.Research on trajectory prediction guidance for launch vehicle［J］.Manned Spaceflight，2013，19 (6):24－28.］

［3］ 宋征宇.信息集成在航天運輸控制系統中的應用研究［J］.宇航學報，2015，36(4)，365－374.［Song Zheng-yu.Application study ofinformation integration in aerospace transportation control system［J］.Journal of Astronautics，2015，36(4):365－374.］

［4］ 楊希祥，張為華.基于Gauss偽譜法的固體運載火箭上升段軌跡快速優化研究［J］.宇航學報，2011，32(1)，15－21.［Yang Xi-xiang，Zhang Wei-hua.Rapid optimization of ascent trajectory for solid launch vehicles based on gauss pseudo-spectral method［J］.Journal of Astronautics，2011，32(1):15－21.］

［5］ 閆循良，廖守億，張金生，等.基于節點改善策略的偽譜軌跡優化［J］.宇航學報，2013，34(7):891－900.［Yan Xun-liang，Liao Shou-yi，Zhang Jin-sheng，et al.Trajectory optimization using pseudospectral method based on a grid node refinement strategy［J］.Journal of Astronautics，2013，34(7):891－900］

［6］ Calise A J，Brqandt N.Generation of launch vehicle abort trajectories using a hybrid optimization method［J］.Journal of Guidance，Control，and Dynamics［J］.2004，27(6)，930－937.

［7］ Mease K D，Chen D T，Teufel P，et al.Reduced-order entry trajectory planning for acceleration guidance［J］.Journal of Guidance，Control，and Dynamics，2002，25(2):257－266.

［8］ 雍恩米，陳磊，唐國金，等.飛行器軌跡優化數值方法綜述［J］.宇航學報，2008，29(2):397－406.［Yong En-mi，Chen Lei，Tang Guo-jin，et al.A survey of numerical methods for trajectory optimization of spacecraft［J］.Journal of Astronautics，2008，29(2):397－406.］

［9］ 徐延萬.彈道導彈、運載火箭控制系統設計與分析［M］.北京:宇航出版社，1999:41－53.

［10］ 張蒿.固體運載火箭彈道設計與優化［J］.固體火箭技術，1997，20(1):1－5.［Zhang Song.Calculation and optimization of solid launch vehicle trajectory［J］.Journal of Solid Rocket Technology，1997，20(1):1－5.］

［11］ Betts J T.Survey of numerical methods for trajectory optimization［J］.Journal of Guidance，Control and Dynamics，1998，21 (2):193－207.

［12］ 邱豐，宋征宇.基于聯立框架的軌跡優化算法［J］.航天控制，2015，33(2):9－15.［Qiu Feng，Song Zheng-yu.A trajectory optimization algorithm based on simultaneous strategies［J］.Aerospace Control，2015，33(2):9－15.］

［13］ 陳偉鋒.大規模復雜過程系統的高性能優化理論與方法研究［D］.杭州:浙江大學，2011.［Chen Wei-feng.Research on high-performance optimization for large-scale complex process systems［D］.Hangzhou:Zhejiang University，2011.］

［14］ Wachter A，Biegler L T.On the implementation of an interiorpoint filterline search algorithm forlarge-scale nonlinear programming［J］.Mathematical Programming，2006，106(1):25－57.

通信地址:北京142信箱402分箱(100854)

電話:(010)88525917

E-mail:qfeng1991@126.com

宋征宇(1970－)，男，碩士，研究員，主要從事飛行器控制、制導與仿真，系統集成、智能自主控制技術等研究。本文通信作者。

通信地址:北京142信箱402分箱(100854)

電話:(010)68389311

E-mail:zycalt12@sina.com

(編輯:張宇平)

Large Terminal Attitude Constrained Trajectory Optimization of Ascent Stage Via Simultaneous Method

QIU Feng，SONG Zheng-yu
(Beijing Aerospace Automatic Control Institute，Beijing 100854，China)

For launch vehicles with constant thrust and without auxiliary propulsion system，a direct method on simultaneous framework(simultaneous method)is presented to solve the dynamic online trajectory planning problem with large terminal attitude constraint，making up for the deficiency of linear programming in the traditional iterative guidance method，meanwhile，the angular rate of the trajectory can be limited according to the requirement of the attitude tracking error，so as to satisfy all constraint conditions.By combing the characteristics of the direct collocation method and pseudo-spectral method，and based on orthogonal collocation on finite element，this paper adopts the simultaneous method and discretized state variables and control variables to solve this problem，which is much more efficient than the pseudo-spectral method.On this basis，the paper puts forward some measures to increase the calculation efficiency，such as initial value selection，orbit precision error adjustment，and adaptive finite element configuration，further improving the real-time performance and providing feasible solutions to more complex and diverse guidance tasks in the future.

Launch Vehicles;Terminal attitude constraint;Simultaneous method;Trajectory design

V4

A

1000-1328(2017)01-0018-08

10.3873/j.issn.1000-1328.2017.01.003

邱 豐(1991－)，女，碩士，主要從事導航、制導與控制。

2016-05-23;

2016-10-16