趙吉松，王江華，王泊喬，張金明，朱航標

(1.南京航空航天大學航天學院，南京 210016；2.北京機電工程研究所，北京 100074)

0 引 言

再入返回技術是載人航天飛行的關鍵和難點之一。與從近地軌道返回不同，月球探測器的返回速度高達11 km/s，火星探測器的返回速度高達14.5 km/s。對于如此高的返回速度，飛行器如果采用直接再入的方式返回地球，最大過載高達17以上，遠遠超出了人的承受能力。跳躍式再入是減小最大過載的一種有效途徑。如圖1所示，飛行器以較小的再入角進入大氣層，減速的同時利用大氣提供的升力，跳出大氣層，在大氣層外作一段橢圓軌道飛行，然后重新再入大氣層[1]。這種跳躍式再入方式通過兩次再入減速，延長了減速時間，從而能夠減小飛行器再入過程的最大過載和最大熱流等[2]。

圖1 跳躍式再入示意圖

再入軌跡優(yōu)化對于探月返回飛行器總體設計和制導控制系統(tǒng)設計具有重要意義[3-4]。與常規(guī)軌跡優(yōu)化問題相比，跳躍式再入軌跡優(yōu)化問題的特殊之處是首次再入段的微小軌跡誤差會在跳出大氣層的慣性飛行段被放大，進而給二次再入段軌跡帶來比較大的偏差。這種二次再入段對首次再入段誤差的高度敏感特性給軌跡優(yōu)化帶來了挑戰(zhàn)。目前的方法為了簡化，對首次再入段和二次再入段分開優(yōu)化[5-6]，但是這種解耦處理方法難以保證軌跡的整體最優(yōu)性。文獻[7]采用高斯偽譜法對探月飛船跳躍式再入軌跡的可達域進行了分析，但是沒有研究軌跡優(yōu)化結果與數(shù)值積分結果的誤差。此外，還有一些文獻從設計的角度對跳躍式再入軌跡的相關參數(shù)影響規(guī)律[8]和能量管理方法[9]等進行了研究，但是這些方法不能充分挖掘跳躍式再入軌跡的性能。

跳躍式再入軌跡優(yōu)化問題本質上屬于最優(yōu)控制問題，其求解方法主要分為間接法和直接法[10]。間接法借助變分法或者最大值原理，把泛函優(yōu)化問題轉化為兩點邊值問題求解；直接法通過對控制變量和/或狀態(tài)變量進行離散把泛函優(yōu)化轉化為非線性規(guī)劃(Nonlinear programming, NLP)問題，然后采用各種非線性規(guī)劃算法求解，比如基于序列二次規(guī)劃(Sequential quadratic programming，SQP)算法的SNOPT[11]和基于內(nèi)點法的IPOPT[12]。直接法中的配點法[10]由于不需要推導最優(yōu)性必要條件，并且對初值的敏感性較低，容易收斂，近年來得到廣泛研究和應用。此外，研究還表明，將配點法與網(wǎng)格細化技術相結合，能夠進一步提高配點法對復雜軌跡優(yōu)化問題的適應能力[13-17]。這種方法的原理是應用網(wǎng)格細化技術在優(yōu)化過程中根據(jù)軌跡的特點動態(tài)調(diào)整離散節(jié)點分布，從而采用較少的離散節(jié)點達到較高的優(yōu)化精度，降低了配點法的計算量。以文獻[17]為例，該研究基于稀疏配點法和網(wǎng)格細化技術對近地軌道返回的高超聲速滑翔再入軌跡進行了優(yōu)化，研究結果表明所述方法能夠快速優(yōu)化出一條嚴格滿足路徑約束和端點約束的再入軌跡。與高超聲速滑翔再入軌跡不同，探月返回飛行器的再入返回軌跡具有跳躍特性并且由于二次再入軌跡對首次再入段的誤差非常敏感，給軌跡優(yōu)化帶來挑戰(zhàn)。

本文在節(jié)點自適應稀疏配點法的基礎上，給出跳躍式再入軌跡的一種高精度優(yōu)化方法，以探月飛行器跳躍式再入返回軌跡為例，通過數(shù)值仿真檢驗了所建立的軌跡優(yōu)化方法的有效性。

1 跳躍式再入軌跡優(yōu)化問題

1.1 再入軌跡運動模型

將飛行器簡化為質點，那么描述飛行器質心運動的微分方程組為(忽略地球自轉影響)[10]

(1)

式中：r，θ，φ分別為飛行器在地心赤道坐標系的矢徑、經(jīng)度、緯度；v，ψ，γ分別為飛行器的速度、航向角和航跡角；g為重力加速度，g=μ/r2，μ為地球引力常數(shù)，μ=3.98603195×1014m3/s2。

空氣動力產(chǎn)生的加速度的沿飛行軌跡切向、法向和側向的三個分量as，an，aw分別為

式中：σ為速度傾側角，m為飛行器的質量。

升力和阻力的表達式分別為

式中：ρ為大氣密度，A為氣動參考面積，CL和CD分別為飛行器的升力系數(shù)和阻力系數(shù)。

1.2 再入初始條件

再入軌跡優(yōu)化問題的初始條件如下

(2)

式中：t0為初始時刻；h為飛行高度，h=r-Re，Re為地球半徑，Re= 6371.20 km；r0，θ0，φ0，v0，ψ0和γ0分別為相應的狀態(tài)變量的初值。

1.3 路徑約束

再入軌跡優(yōu)化問題的控制變量為速度傾側角，即u=σ(t)，其變化范圍受到如下限制

σmin≤σ(t)≤σmax

(3)

式中：σmin和σmax為速度傾側角的邊界。

飛行器再入返回過程的過載直接關系到宇航員的生命安全和舒適度，需要對過載進行限制

(4)

式中：nmax為再入過程允許的最大過載。

為了使得飛行器安全返回，需要對高超聲速再入過程的對流氣動加熱進行限制，即

(5)

此外，考慮到飛行器結構的承受能力，還需要對再入過程的動壓進行限制，即

(6)

式中：qmax為再入過程允許的動壓上限。

1.4 末端約束

當飛行器的速度降至預定速度時，阻力降落傘打開，飛行器進一步減速下降。為了安全著陸，還需要對開傘時的飛行高度進行約束。綜上所述，再入飛行器軌跡優(yōu)化問題的終端約束條件為

v(tf)=vf,h(tf)≥hf, min

(7)

式中：tf為再入軌跡的終端時刻，vf為開傘時飛行器的速度，hf, min允許的最低開傘高度。

1.5 目標函數(shù)

軌跡優(yōu)化的目標函數(shù)可根據(jù)實際設計需求選取。最大橫向航程是衡量飛行器再入機動能力的重要指標。本文選取最大橫向航程作為優(yōu)化指標。由于φ0= 0°，ψ0= 0°，因而目標函數(shù)可寫為

J=minφ(tf)

(8)

式中：φ(tf)<0，式(8)使橫向航程最大化。

綜上所述，跳躍式再入返回軌跡優(yōu)化問題描述為：求解最優(yōu)控制變量u(t)=σ(t)，使得目標函數(shù)最小化，并且滿足再入動力學方程組(1)，初始條件(2)，軌跡路徑約束(3)～(5)以及終端約束(7)。

2 跳躍式再入軌跡優(yōu)化方法

2.1 自適應稀疏局部配點法

(9)

狀態(tài)方程為

(10)

端點條件為

E(x(t0),t0,x(tf),tf)=0

(11)

路徑約束為

C(x(t),u(t),t)≤0,t∈[t0,tf]

(12)

本文采用局部配點法求解軌跡優(yōu)化問題。首先利用積分變換τ=(t-t0)/(tf-t0)將軌跡優(yōu)化問題(方程(9)～(12))的時間區(qū)間變換至歸一化的時間區(qū)間τ∈[0, 1]。假設單位區(qū)間[0, 1]上的N個離散節(jié)點為

G={τi:τi∈[0,1],i=0,1,…,N;τ0=0,τN=τf=1;τi<τi+1,i=0,1, …,N-1}

(13)

式中：τi稱為離散節(jié)點或網(wǎng)格節(jié)點，τi在[0, 1]上可以均勻分布，也可以非均勻分布。

記xi=x(τi)，ui=u(τi)，對于狀態(tài)方程(10)，基于q階Runge-Kutta(RK)方法的離散格式為

(14)

式中：Δt=tf-t0，hi=τi+1-τi，fij=f(xij,uij,τij；t0,tf)，xij，uij和τij為中間變量，xij由下式給出

(15)

式中：τij=τi+hiρj，uij=u(τij)，ρj,βj,αjl均為已知常數(shù)并且滿足0≤ρ1≤…≤ρq≤1。當αjl= 0(l≥j)時，為顯式格式，否則為隱式格式。采用類似的方法，可將目標函數(shù)離散化。常用的離散格式有梯形格式(q= 2)，Hermite-Simpson格式(q= 3，簡記為HS格式)，經(jīng)典四階Runge-Kutta格式(q= 4)等。

(16)

并且滿足如下約束

(17)

Ci=C(xi,ui,τi;t0,tf)≤0

(18)

(19)

E(x0,t0,xf,tf)=0

(20)

本研究采用HS格式，該格式需要用到區(qū)間中點的變量和函數(shù)值，為此將區(qū)間中點的控制量也作為優(yōu)化變量并且在區(qū)間中點添加路徑約束，即

(21)

(22)

約束條件為

(23)

Ci=C(xi,ui,τi;t0,tf)≤0

(24)

(25)

E(x0,t0,xf,tf)=0

(26)

其中

在數(shù)值優(yōu)化時，為了使軌跡優(yōu)化問題具有實際物理意義，還需要添加時間差約束

Δt=tf-t0>0

(27)

求解方程(22)～(27)所示的NLP問題即可得到軌跡優(yōu)化問題的最優(yōu)解。本文采用美國斯坦福大學基于SQP算法開發(fā)的SNOPT[11]求解NLP問題。

為了提高軌跡優(yōu)化的求解效率和精度，本文應用稀疏差分算法[18]提高NLP偏導數(shù)的計算效率，應用節(jié)點細化技術[19]動態(tài)調(diào)整離散節(jié)點分布，提高對跳躍式再入軌跡的適應能力。在采用SNOPT求解NLP時，為SNOPT提供NLP的一階偏導數(shù)矩陣(即雅克比矩陣，定義為目標函數(shù)和全部約束對全部自變量的偏導數(shù)矩陣)能夠顯著提高優(yōu)化的計算效率，但是偏導數(shù)矩陣的計算量比較大。稀疏差分法[18]通過分析偏導數(shù)矩陣的稀疏特性，將其中的占絕大多數(shù)的零元素識別出來，并且將其中的非零元素的值通過矩陣鏈式求導運算分解為原始軌跡優(yōu)化問題的偏導數(shù)，然后采用稀疏差分方法計算，從而大幅度減小偏導數(shù)的計算量。節(jié)點細化技術[19]的原理是數(shù)據(jù)壓縮，每次優(yōu)化出一條軌跡之后，根據(jù)軌跡的變化特性，在光滑區(qū)域對離散節(jié)點進行壓縮，減少節(jié)點數(shù)量，在非光滑區(qū)域插入一些新節(jié)點，使之分布更密，然后再次優(yōu)化軌跡，綜合效果是采用較少的離散節(jié)點高精度描述軌跡。稀疏差分法和節(jié)點細化技術都具有很強的通用性和魯棒性，對于本文的優(yōu)化問題，只需要設置用于判斷軌跡是否光滑的閾值參數(shù)即可(判斷每個節(jié)點處的控制量與采用其周圍節(jié)點的插值結果的差異是否超過閾值，若超過，則認為軌跡不光滑，否則認為軌跡光滑)。

2.2 跳躍式再入軌跡優(yōu)化方法

對于跳躍式再入軌跡，由于二次再入段對首次再入段的誤差非常敏感，若將首次再入段和二次再入段作為整體進行優(yōu)化會導致二次再入軌跡的優(yōu)化精度較低。本文在自適應稀疏局部配點法的基礎上，提出一種能夠高精度求解跳躍式再入軌跡優(yōu)化問題的方法。該方法的關鍵之處是從跳躍軌跡的最高點開始，對二次再入軌跡進行重新優(yōu)化。如圖2所示，首先，將跳躍式再入軌跡作為整體進行優(yōu)化，得到最優(yōu)傾側角隨時間變化曲線；然后，按照最優(yōu)傾側角積分再入動力學方程組得到最優(yōu)軌跡；當積分至跳躍軌跡的最高點時，以實際積分得到的狀態(tài)變量作為新的初始條件(其它條件保持不變)，對二次再入軌跡重新優(yōu)化，求解新的最優(yōu)傾側角?？紤]到優(yōu)化計算需要耗費一些時間，為了使得方法在制導跟蹤領域具有實用性，本文方法在解出二次優(yōu)化問題的最優(yōu)傾側角之前，仍然采用整體優(yōu)化得到的最優(yōu)傾側角作為二次再入軌跡的控制輸入。因為在跳躍軌跡的最高點附近大氣非常稀薄，氣動力對軌跡的影響非常小，所以這樣近似處理是合理的。

圖2 跳躍式再入軌跡優(yōu)化策略

本文方法的具體流程如下：

1)應用基于稀疏差分法和節(jié)點自適應細化技術的局部配點法求解由方程(1)～(8)描述的跳躍式再入軌跡優(yōu)化問題，得到最優(yōu)傾側角σ1(t)。

2)以σ1(t)為控制變量，采用四階Runge-Kutta方法積分由方程組(1)～(2)描述的微分方程組初值問題至跳躍軌跡的最高點(即最大高度處)，得到最高點的狀態(tài)變量x1,f以及對應的時刻t1,f。

3)以t1,f和x1,f作為新的初始條件，對由方程組(1)和(3)～(8)描述的再入軌跡優(yōu)化問題進行重新優(yōu)化(即二次優(yōu)化)，得到新的最優(yōu)傾側角σ2(t)和終端時刻t2,f，記下二次優(yōu)化的耗時為Δt2。

4)從跳躍軌跡的最高點開始，以x1,f為初始條件，以σ1(t)為控制變量，繼續(xù)積分再入動力學方程組(1)至t2,1=t1,f+ Δt2，得到狀態(tài)變量x2,1。

5)以t2,1和x2,1為初始條件，以σ2(t)為控制變量，積分再入動力學方程組(1)至終端時刻t2,f。

6)將在步驟2)、4)和5)中積分得到的三段軌跡(含各段采用的控制變量)組合在一起，即可得到最優(yōu)跳躍式再入軌跡和相應的控制變量。

相對于間接法和直接打靶法等其它方法，配點法的優(yōu)勢之一是對優(yōu)化初值的敏感性較低。盡管如此，良好的優(yōu)化初值仍然有利于加速軌跡優(yōu)化的收斂。本文在對跳躍式再入軌跡進行第一次優(yōu)化時，狀態(tài)變量的優(yōu)化初值選取初始狀態(tài)變量和終端狀態(tài)變量的連線(對于終端狀態(tài)給定的變量)或者初值狀態(tài)變量常值(對于終端狀態(tài)沒有給定的變量)，控制變量的優(yōu)化初值取σ(t)=-45°。在后續(xù)的節(jié)點細化以及二次優(yōu)化時，均利用前一步得到的最優(yōu)軌跡通過插值提供較為準確的優(yōu)化初值。

3 仿真結果

3.1 方法1(整體優(yōu)化)

為了對比，本文首先不考慮二次優(yōu)化，在整個時間區(qū)間對跳躍式再入軌跡進行優(yōu)化。采用本文3.1節(jié)介紹的自適應稀疏局部配點法求解該問題。圖3為優(yōu)化出的最優(yōu)控制變量隨時間變化曲線。圖4為優(yōu)化出的最優(yōu)高度和速度隨時間變化曲線。圖3和圖4中的圓圈表示離散最優(yōu)解(總共采用142個非均勻分布的離散節(jié)點)。圖3中的細實線為對離散控制變量進行插值得到的連續(xù)控制變量隨時間變化曲線。圖4中的細實線為根據(jù)圖3所示的最優(yōu)控制變量采用四階Runge-Kutta數(shù)值積分方法對再入動力學方程組進行數(shù)值積分得到的結果。

圖3 最優(yōu)傾側角變化曲線(方法1)

圖4 最優(yōu)高度和速度變化曲線(方法1)

由圖4可知，在初次再入階段，數(shù)值積分結果與離散最優(yōu)解的差異微小，但是從飛行器跳躍出大氣層開始，數(shù)值積分結果與離散最優(yōu)解逐漸出現(xiàn)了明顯的差異。在軌跡終端，數(shù)值積分得到的終端高度為7361.8 m，終端速度為106.1 m/s，均不滿足開傘條件。數(shù)值積分結果與離散最優(yōu)解的高度差異高達-4768.3 m，速度差異高達-43.9 m/s，顯然超出了合理范圍。這種差異主要是因為在初次再入段，數(shù)值積分結果與離散最優(yōu)解存在微小的誤差(由于數(shù)值離散造成的，不可避免)，而跳躍式再入軌跡的跳躍段和二次再入段對首次再入段的誤差非常敏感，使得首次再入段的誤差沿軌跡累積并傳播，最終導致軌跡的終端誤差過大。對于常規(guī)軌跡優(yōu)化問題，數(shù)值積分結果與離散最優(yōu)解也存在誤差，只是誤差通常非常小，其影響可以忽略不計。理論上，增加離散節(jié)點的數(shù)量能夠減小這種誤差，但是實際上如果離散節(jié)點數(shù)量過多，NLP的規(guī)模過大，優(yōu)化算法的收斂性通常會變差，甚至無法收斂[21]。

3.2 方法2(整體優(yōu)化+二次優(yōu)化)

為了解決數(shù)值積分結果和離散最優(yōu)解差異過大問題，本文在跳躍式再入軌跡整體優(yōu)化結果的基礎上，從跳躍的最高點(對應的時間t1,f=711.2 s)開始，對二次再入軌跡重新優(yōu)化，具體方法參見2.2節(jié)。圖5給出優(yōu)化的控制變量隨時間變化曲線，圖6給出優(yōu)化的高度和速度隨時間變化曲線。其中，跳躍最高點之前的軌跡是通過對再入軌跡整體優(yōu)化得到(參見2.1節(jié))，跳躍最高點之后的軌跡為二次優(yōu)化的結果(本文考慮了二次優(yōu)化的耗時，因而連接點相對最高點略微右移)。圖5和圖6中的圓圈表示二次優(yōu)化得到的離散最優(yōu)解(采用84個非均勻分布的節(jié)點)。圖5中的細實線是對離散控制變量進行插值得到的連續(xù)控制變量隨時間變化曲線，圖6中的細實線是根據(jù)圖5所示的最優(yōu)控制變量對再入動力學方程組進行數(shù)值積分得到的結果。

對比圖5和前述圖3可知，對完整再入軌跡優(yōu)化和對二次再入軌跡重新優(yōu)化得到的二次再入控制變量差別不大。對比圖6和圖4可知，引入二次優(yōu)化之后，二次再入軌跡的數(shù)值積分結果與離散最優(yōu)解非常接近，終端高度誤差僅為2.12218 m，終端速度誤差僅為0.02580 m/s，很好地滿足了再入飛行器的開傘條件約束?？梢?，通過對二次再入段重新優(yōu)化顯著提高了跳躍式再入軌跡的優(yōu)化精度。

圖5 最優(yōu)控傾側角變化曲線(方法2)

圖6 最優(yōu)高度和速度變化曲線(方法2)

圖7給出優(yōu)化的三維再入軌跡，其中圓圈為離散最優(yōu)解，細實線為數(shù)值積分結果。為了便于展示軌跡特征，圖7中還給出了飛行器在地球表面投影點的軌跡?？梢?，為了取得最大橫向航程，飛行器在首次再入過程中利用氣動力改變速度方向進行轉彎，在二次再入過程中進一步利用氣動力進行橫向機動。圖8給出再入過程中過載等路徑約束沿軌跡變化曲線?？梢?，最優(yōu)軌跡的過載、熱流和動壓都滿足路徑約束要求，其中過載在首次再入段的最低點附近達到上限，其它約束均未達到上限。

圖7 三維最優(yōu)跳躍式再入軌跡(方法2)

圖8 最優(yōu)軌跡的路徑約束(方法2)

表1給出方法1(整體優(yōu)化)和方法2(整體優(yōu)化+二次優(yōu)化)的優(yōu)化結果對比。對于方法2，表中還給出了不同的二次優(yōu)化耗時Δt2對數(shù)值積分終端誤差的影響。從表1中可以看出，通過對二次再入軌跡重新優(yōu)化，可以將終端時刻的速度誤差和高度誤差降低3～4個數(shù)量級，從而使得飛行器的開傘條件得到嚴格滿足。在目標函數(shù)方面，兩種方法解出的最大橫向航程幾乎沒有差別。需要強調(diào)的是，方法1的優(yōu)化結果對應的數(shù)值積分終端誤差過大、不滿足開傘條件，實際上并不是可行軌跡。

在本算例中，二次再入軌跡重新優(yōu)化消耗的時間為Δt2= 2.4 s。為了使方法具有實用性，本文在二次優(yōu)化得到新的控制變量之前，仍然采用整體優(yōu)化得到的控制變量作為輸入對再入動力學方程組進行積分。雖然二次優(yōu)化得出的控制變量和整體優(yōu)化得出的控制變量存在差異，但是仿真結果表明這樣處理對二次再入軌跡幾乎沒有影響。目前，國產(chǎn)宇航芯片的主頻可達到300 MHz[22]，而本文采用的計算平臺的處理器主頻為1.6 GHz，考慮到內(nèi)存和軟件等性能的差異，預計宇航級計算機的計算性能與本文采用的計算平臺相差10～20倍。進一步仿真研究表明，即使在計算能力較低的宇航級計算機上(以計算能力降低至1/20為例)，二次優(yōu)化的耗時會顯著增加，但是上述處理方式對二次再入軌跡仍然幾乎沒有影響，如表1所示。因為跳躍最高點的高度為117.8 km，大氣非常稀薄，氣動力非常小，因而在最高點附近采用不夠準確的控制變量積分再入動力學方程組對軌跡幾乎沒有影響。因此，本文提出的處理方法使得二次再入軌跡的優(yōu)化是一種準實時優(yōu)化，在再入軌跡的制導領域具有應用潛力。

表1 不同方法的優(yōu)化結果對比

4 結論

探月返回飛行器跳躍式再入返回軌跡優(yōu)化問題是復雜的多約束、非線性最優(yōu)控制問題，特別是二次再入軌跡對首次再入段的誤差非常敏感，給優(yōu)化帶來挑戰(zhàn)。本文提出一種高精度求解跳躍式再入軌跡優(yōu)化問題的方法。該方法首先應用節(jié)點自適應稀疏配點法對完整跳躍式再入軌跡進行優(yōu)化，然后根據(jù)優(yōu)化得到的控制變量積分再入動力學方程組，當積分到跳躍軌跡的最高點時，以積分得到的狀態(tài)變量值作為新的初始條件，對二次再入軌跡重新優(yōu)化。在求解二次再入軌跡期間(約2～3 s)，繼續(xù)基于整體優(yōu)化解出的最優(yōu)控制積分再入動力學方程組，直到優(yōu)化算法重新優(yōu)化出二次再入段的最優(yōu)控制，然后根據(jù)新的最優(yōu)控制積分再入動力學方程組至軌跡終端。仿真結果表明：與不采用二次優(yōu)化的方法相比，采用二次優(yōu)化能夠將跳躍式再入軌跡的終端高度和速度誤差降低3～4個數(shù)量級。此外，在跳躍的最高點附近進行的二次優(yōu)化是一種準實時優(yōu)化，在再入軌跡的精確制導領域具有應用潛力。