助推－滑翔導(dǎo)彈再入彈道快速優(yōu)化①

趙 欣，閆循良，張金生，王仕成，何安榮

(1.第二炮兵工程大學(xué)，西安 710025;2.第二炮兵駐211廠軍代室，北京 100084)

助推－滑翔導(dǎo)彈再入彈道快速優(yōu)化①

趙 欣1，閆循良1，張金生1，王仕成1，何安榮2

(1.第二炮兵工程大學(xué)，西安 710025;2.第二炮兵駐211廠軍代室，北京 100084)

基于LGL(Legendre-Gauss-Lobatto)偽譜法，研究了臨近空間助推-滑翔導(dǎo)彈再入段彈道快速優(yōu)化問題。首先，基于改進(jìn)的氣動(dòng)力模型建立了較為精確的再入數(shù)學(xué)模型;其次，針對(duì)該優(yōu)化問題在氣動(dòng)數(shù)據(jù)處理和優(yōu)化求解上存在的困難，基于LGL偽譜法系統(tǒng)地建立了再入最優(yōu)飛行彈道的求解步驟，為解決直接利用LGL偽譜法存在的困難，設(shè)計(jì)了一種基于LGL偽譜法的串行優(yōu)化求解策略;最后，分別采用積分推進(jìn)法和協(xié)狀態(tài)映射原理對(duì)優(yōu)化結(jié)果進(jìn)行了可行性和最優(yōu)性驗(yàn)證。仿真結(jié)果表明，本文的彈道優(yōu)化方法優(yōu)化1條再入彈道所用時(shí)間為3～4 s，計(jì)算效率較高，路徑約束和端點(diǎn)約束均得到很好滿足，算法求解精度較高，有效地實(shí)現(xiàn)了多約束多變量大型稀疏的再入彈道導(dǎo)彈快速優(yōu)化。

臨近空間;助推-滑翔導(dǎo)彈;再入彈道優(yōu)化;偽譜方法;優(yōu)化控制

0 引言

當(dāng)前，由于導(dǎo)彈防御系統(tǒng)彈道預(yù)測(cè)性不斷增強(qiáng)，彈道導(dǎo)彈所面臨的生存環(huán)境越來越惡劣。近些年，為提高導(dǎo)彈的戰(zhàn)場(chǎng)生存能力，各軍事強(qiáng)國競(jìng)相關(guān)注并研究中段突防變軌技術(shù)［1］，而助推-滑翔導(dǎo)彈正是在這一技術(shù)發(fā)展中應(yīng)運(yùn)而生的新概念武器［2］。助推-滑翔導(dǎo)彈采用升力體形，依靠氣動(dòng)力控制，可在臨近空間區(qū)域長時(shí)間高超聲速滑翔飛行，從而實(shí)現(xiàn)遠(yuǎn)程快速攻擊，其地位具有深遠(yuǎn)戰(zhàn)略意義。

在臨近空間助推-滑翔導(dǎo)彈的研究中，再入彈道優(yōu)化是系統(tǒng)總體優(yōu)化中的關(guān)鍵技術(shù)之一，它直接影響武器的整體突防性能。然而，由于這類導(dǎo)彈的再入彈道對(duì)控制變量高度敏感，且再入過程的非線性約束較強(qiáng)，致使彈道的可行域范圍較窄，這決定其再入彈道數(shù)值優(yōu)化存在一定困難［3］。目前，求解彈道優(yōu)化問題的主流方法有:基于極大值原理的間接法和基于非線性規(guī)劃理論的直接法［4］。間接法將最優(yōu)控制問題最終轉(zhuǎn)換為一個(gè)兩點(diǎn)邊值問題［5］，可以獲得精確的最優(yōu)解，但其要求求解狀態(tài)方程和協(xié)態(tài)方程，這對(duì)于復(fù)雜的非線性動(dòng)力學(xué)方程來說較為麻煩，計(jì)算量較大;同時(shí)，求解過程對(duì)共軛變量的初值高度敏感且難以估計(jì)，很難收斂，且算法魯棒性較差。直接法采用參數(shù)化方法將連續(xù)空間的最優(yōu)控制問題求解轉(zhuǎn)化為一個(gè)非線性規(guī)劃(Non Linear Programming，NLP)問題，通過數(shù)值求解非線性規(guī)劃問題來獲得最優(yōu)軌跡。與間接法相比，直接法具有以下優(yōu)點(diǎn):(1)算法簡單、易于理解;(2)對(duì)初始值要求不嚴(yán)格，收斂半徑大，魯棒性好;(3)易于處理包含路徑約束的最優(yōu)控制問題［6］。因此，直接法在解決實(shí)際問題的適應(yīng)性上更具有優(yōu)勢(shì)。

直接法又可分為2種基本類型:一是離散控制變量(又稱直接打靶法)［7］;二是同時(shí)離散控制變量和狀態(tài)變量(又稱直接配點(diǎn)法)［8］。前一種方法的主要不足點(diǎn)是節(jié)點(diǎn)數(shù)目不易確定，取得過少會(huì)導(dǎo)致所得到的最優(yōu)控制解不準(zhǔn)確，過多則會(huì)增大非線性規(guī)劃問題的規(guī)模，出現(xiàn)“維數(shù)災(zāi)難”;后一種方法中，一般的直接配點(diǎn)法由于采用分段多項(xiàng)式來近似狀態(tài)和控制變量時(shí)間歷程，使得設(shè)計(jì)變量數(shù)目龐大，且初值不易給定，控制曲線不光滑。

近些年發(fā)展起來的偽譜法(PSM)開始應(yīng)用于軌跡優(yōu)化問題，它也是配點(diǎn)法的一種，它采用全局插值多項(xiàng)式有限基在一系列離散點(diǎn)上近似狀態(tài)變量和控制變量，相對(duì)于一般的配點(diǎn)法，PSM能以較少的節(jié)點(diǎn)獲得較高的精度［9-10］。同時(shí)，采用PSM轉(zhuǎn)化得到的NLP問題的KKT(Karush-Kuhn-Tucker)條件與原最優(yōu)控制問題一階最優(yōu)必要條件的離散形式具有一致性［11］，彌補(bǔ)了一般直接法的不足。上述特點(diǎn)使得PSM成為目前求解復(fù)雜最優(yōu)控制問題最為有效的方法之一。

本文研究Legendre-Gauss-Lobatto(LGL)PSM在臨近空間助推-滑翔導(dǎo)彈再入彈道優(yōu)化中的應(yīng)用。首先，從工程應(yīng)用角度出發(fā)建立了與真實(shí)飛行環(huán)境比較相符的臨近空間助推-滑翔導(dǎo)彈以再入時(shí)間最短為目標(biāo)的多約束彈道優(yōu)化模型;在此基礎(chǔ)上，采用基于LGL PSM和序列二次規(guī)劃(SQP)算法的串行優(yōu)化求解策略處理問題;最后，采用解析法對(duì)優(yōu)化結(jié)果進(jìn)行了可行性和最優(yōu)性驗(yàn)證。通過仿真對(duì)上述方法進(jìn)行了數(shù)值驗(yàn)證。

1 再入彈道優(yōu)化問題數(shù)學(xué)模型

所研究問題的數(shù)學(xué)模型包括再入動(dòng)力學(xué)模型、氣動(dòng)力模型和大氣模型。

1.1 再入動(dòng)力學(xué)模型

基于相關(guān)假設(shè)，考慮地球自轉(zhuǎn)效應(yīng)，升力式再入飛行器三自由度無量綱動(dòng)力學(xué)方程為

式中r為再入飛行器質(zhì)心到地球球心的無量綱地心距;θ和φ分別為再入飛行器質(zhì)心位置的經(jīng)度和緯度;V為相對(duì)地球的無量綱速度;γ為彈道傾角;ψ為速度方位角;Ω為無量綱地球旋轉(zhuǎn)角速度;σ為傾側(cè)角;D和L分別為無量綱的阻力加速度和升力加速度。

1.2 改進(jìn)的氣動(dòng)力模型

目前，再入軌跡優(yōu)化中的氣動(dòng)力模型一般采用傳統(tǒng)簡化形式，即只將氣動(dòng)系數(shù)視為攻角函數(shù)［3］。然而，實(shí)際上氣動(dòng)系數(shù)受多種因素的影響，其中攻角α和飛行速度Vm是主要的2個(gè)因素，為進(jìn)一步減小因模型不精確引入的誤差，需給出更為精確的氣動(dòng)力模型。文獻(xiàn)［12］研究指出，升力系數(shù)與攻角近似呈線性關(guān)系，與馬赫數(shù)近似呈負(fù)指數(shù)變化關(guān)系;阻力系數(shù)與攻角近似呈二次函數(shù)關(guān)系，與馬赫數(shù)也近似呈負(fù)指數(shù)變化關(guān)系?；诖?，給出如下的改進(jìn)的氣動(dòng)力系數(shù)模型:

其中，d0、d1、d2、d3、l0、l1、l2、l3均為模型的待辨識(shí)參數(shù)，本文采用非線性最小二乘法進(jìn)行辨識(shí)。

1.3 大氣模型

大氣密度模型一般采用1976標(biāo)準(zhǔn)大氣模型，該模型為一分段函數(shù)，給出了大氣密度和聲速隨高度變化的函數(shù)。具體模型見文獻(xiàn)［13］。

2 再入彈道優(yōu)化模型

2.1 再入狀態(tài)變量和控制變量選擇

再入優(yōu)化問題狀態(tài)變量為

不考慮指令延遲，控制變量可以選擇為攻角α和傾側(cè)角σ，因此得到控制空間為

根據(jù)以上分析，得到優(yōu)化變量為(r，θ，φ，V，γ，ψ，α，σ)共8個(gè)參數(shù)。若設(shè)計(jì)攻角剖面為馬赫數(shù)的分段線性函數(shù)，則優(yōu)化控制量只有σ。因此，新的優(yōu)化變量為(r，θ，φ，V，γ，ψ，σ)，共 7 個(gè)參數(shù)。

2.2 再入彈道約束

再入飛行過程是一個(gè)高動(dòng)態(tài)的過程，需滿足一系列苛刻約束條件的限制，如過載限制、動(dòng)壓限制、熱流限制和初終端條件限制等。

(1)熱流約束

再入熱流約束一般是指對(duì)飛行器表面駐點(diǎn)處的熱流限制，其數(shù)學(xué)表達(dá)式為

駐點(diǎn)熱流計(jì)算式為

式中 ρ分別為無量綱大氣密度;RN為駐點(diǎn)曲率半徑;k、ρ0、Vc為常數(shù);v為飛行速度;n由邊界層的特性決定，一般情況下把駐點(diǎn)附近的邊界層按照層流處理，此時(shí)n=1/2;m則需根據(jù)氣體的粘性等因素來確定，本文取3.25;再入飛行熱流限制取=7 000 kW/m2。

(2)過載約束

最大總過載限制主要取決于再入飛行器結(jié)構(gòu)強(qiáng)度和機(jī)載設(shè)備的可承受過載范圍。再入飛行器總過載表達(dá)式滿足:

式中nmax為額定最大值，文中取20gn。

(3)動(dòng)壓約束

再入飛行過程中的動(dòng)壓約束一方面取決于隨動(dòng)壓增大而增大的氣動(dòng)控制鉸鏈力矩;另一方面取決于飛行器表面防熱材料的強(qiáng)度。動(dòng)壓約束為

其中qmax為最大動(dòng)壓，文中取qmax=3×105N/m2。

(4)端點(diǎn)條件

再入點(diǎn)初始時(shí)刻需滿足一定的條件，同時(shí)考慮到與末制導(dǎo)段交班的要求，再入終端狀態(tài)也應(yīng)滿足一定的條件，端點(diǎn)約束條件為

2.3 優(yōu)化目標(biāo)函數(shù)

對(duì)于導(dǎo)彈武器系統(tǒng)而言，再入點(diǎn)和目標(biāo)點(diǎn)確定后，飛行時(shí)間是一個(gè)重要指標(biāo)，即要在盡量短的時(shí)間內(nèi)攻擊目標(biāo)，提高武器的快速打擊性能。所以，選取到達(dá)目標(biāo)點(diǎn)飛行時(shí)間最短作優(yōu)化性能指標(biāo):

3 基于LGL偽譜法的再入段彈道優(yōu)化問題求解

3.1 彈道優(yōu)化問題描述

上述再入彈道優(yōu)化問題可歸納為以Bolza型為優(yōu)化目標(biāo)的非線性最優(yōu)控制問題:

式(17)～式(21)分別對(duì)應(yīng)彈道優(yōu)化問題的目標(biāo)函數(shù)、動(dòng)力學(xué)方程、邊界條件、路徑約束、狀態(tài)量和控制量自身范圍約束。

3.2 基于LGL偽譜法彈道優(yōu)化問題的轉(zhuǎn)化

(1)區(qū)間變換處理

軌跡優(yōu)化問題的時(shí)間區(qū)間為［t0，tf］，而LGL偽譜法通過τ∈［-1，1］上的拉格朗日多項(xiàng)式擬合控制量和狀態(tài)量。因此，需將時(shí)間區(qū)間轉(zhuǎn)換至［-1，1］，最終得到上述最優(yōu)控制問題的新形式為

(2)時(shí)間區(qū)間劃分及LGL點(diǎn)

設(shè)將時(shí)間歷程劃分為N個(gè)子區(qū)間，則有N+1個(gè)LGL 點(diǎn) τi(i=0，1，…，N)，其中 τ0=-1，τN= τf=1;中間節(jié)點(diǎn) τi(i=1，2，…，N-1)為N階 Legendre多項(xiàng)式的導(dǎo)數(shù)˙LN的零點(diǎn)。其中，LN為N階Legendre正交多項(xiàng)式。這些點(diǎn)在［-1，1］上并不是均勻分布的，而是兩端密集，中間稀疏，有利于提高擬合精度。

(3)控制變量和狀態(tài)變量近似

對(duì)狀態(tài)變量和控制變量基于LGL點(diǎn)進(jìn)行拉格朗日插值近似

式中 φi(·)為N階拉格朗日多項(xiàng)式。

在LGL點(diǎn)可得到以下關(guān)系式:

(4)動(dòng)力學(xué)微分方程約束處理

對(duì)連續(xù)微分約束用插值多項(xiàng)式的微分來代替，對(duì)狀態(tài)近似方程進(jìn)行微分，在LGL點(diǎn)進(jìn)行擬合得:

式中Dki為(N+1)×(N+1)階微分矩陣D的元素。

動(dòng)力學(xué)微分方程約束可近似轉(zhuǎn)換為代數(shù)約束，在LGL點(diǎn)擬合得到:

(5)性能指標(biāo)近似

采用Gauss-Lobatto積分準(zhǔn)則對(duì)性能指標(biāo)中的積分項(xiàng)進(jìn)行數(shù)值積分，得到離散的性能指標(biāo)函數(shù):

其中，LGL權(quán):

通過以上對(duì)最優(yōu)控制問題在LGL點(diǎn)近似，得到了以下離散化的非線性規(guī)劃問題，即尋找最優(yōu)參數(shù)X=(x0，x1，…，xN)，U=(u0，u1，…，uN)，t0，tf，使得性能指標(biāo)JN(X，U)最小化，同時(shí)滿足代數(shù)方程約束、端點(diǎn)約束、路徑等約束。其中，優(yōu)化變量數(shù)目為(n+m)(N+1)+2，約束數(shù)目為(n+c+m)(N+1)+q?？梢?，上述離散得到的問題為一大型稀疏NLP問題，對(duì)于這類問題的求解，本文采用序列二次規(guī)劃(SQP)算法，它被認(rèn)為是目前求NLP問題最有效的方法之一［14］。

3.3 基于LGL偽譜法的串行快速優(yōu)化策略

理論上，上述LGL偽譜法可直接求解本文研究的問題，然而實(shí)際應(yīng)用中存在以下因難:再入軌跡優(yōu)化由于狀態(tài)量和控制量數(shù)目較多，如果選取節(jié)點(diǎn)數(shù)目較多時(shí)，優(yōu)化設(shè)計(jì)變量的數(shù)目將會(huì)十分龐大，同時(shí)得到的約束方程數(shù)目也是非常巨大的，這會(huì)大大增加計(jì)算負(fù)擔(dān)。另外，對(duì)于大范圍再入問題，設(shè)計(jì)變量的初始猜測(cè)較為復(fù)雜，不恰當(dāng)?shù)某踔禃?huì)帶來計(jì)算代價(jià)的增加，甚至可能導(dǎo)致問題不收斂或收斂到不可行解。針對(duì)再入彈道優(yōu)化存在的問題，為了提高算法的計(jì)算速度及收斂性，本文采用基于LGL偽譜法的串行快速優(yōu)化策略來進(jìn)行最優(yōu)彈道的求解。

(1)設(shè)計(jì)初值生成器

利用LGL偽譜法能以較少節(jié)點(diǎn)獲得高精度的優(yōu)勢(shì)，先采用較少的節(jié)點(diǎn)，計(jì)算最優(yōu)軌跡和最優(yōu)控制，并以此為下一步精確計(jì)算的初值。

(2)采用從可行解到最優(yōu)解的串行優(yōu)化策略

串行優(yōu)化策略是指，不直接找尋滿足所有等式約束和不等式約束的解，而是先將等式約束轉(zhuǎn)換為目標(biāo)函數(shù)，將得到的最優(yōu)解x0作為初值，求解原非線性問題的解。

具體流程:首先，基于從可行解到最優(yōu)解的策略，由初值生成器獲得初值;再選取較多的節(jié)點(diǎn)，以求得高精度的最優(yōu)軌跡和最優(yōu)控制。

3.4 可行性和最優(yōu)性驗(yàn)證

為了證明極值解的合理性，研究中進(jìn)一步完成了優(yōu)化解的可行性和最優(yōu)性驗(yàn)證工作。

(1)基于積分推進(jìn)法的可行性驗(yàn)證

采用積分推進(jìn)法進(jìn)行極值解的可行性驗(yàn)證，具體方法:在保證所有約束滿足的前提下，利用所得到的控制量，積分推進(jìn)運(yùn)動(dòng)方程得到狀態(tài)量軌跡，通過比較優(yōu)化獲得的狀態(tài)量和積分推進(jìn)得到的狀態(tài)量軌跡，進(jìn)行解的可行性驗(yàn)證。如果二者匹配程度在可接受的誤差范圍內(nèi)，證明解是可行的。

(2)最優(yōu)性驗(yàn)證

采用Covector Mapping Principle(CMP)進(jìn)行一階最優(yōu)性驗(yàn)證，即根據(jù)龐特里亞金極小值原理，增廣哈米爾頓函數(shù) ˉH滿足 ?ˉH/?u=0，將其重新改寫，得到控制律的顯示表達(dá)，將解析得到的控制量與優(yōu)化得到的控制量進(jìn)行比較，如果二者是一致的，則一階最優(yōu)性是滿足的;同時(shí)，根據(jù)KKT定理可得到相應(yīng)的補(bǔ)充條件，判斷補(bǔ)充條件是否成立，最終實(shí)現(xiàn)結(jié)果的最優(yōu)性驗(yàn)證。

4 仿真算例

4.1 再入仿真條件

再入端點(diǎn)條件為

其中，考慮上述初值設(shè)定，對(duì)應(yīng)的狀態(tài)邊界條件為

對(duì)應(yīng)的各控制量邊界為

路徑約束已經(jīng)在前文中給出。若設(shè)計(jì)攻角剖面為馬赫數(shù)的分段線性函數(shù)，即

則優(yōu)化控制量只有σ，因此新的優(yōu)化變量為(r，θ，φ，γ，ψ，σ)，共7 個(gè)參數(shù)。

4.2 仿真結(jié)果及分析

采用初值生成技術(shù)和從可行解到最優(yōu)解的串行優(yōu)化技術(shù)進(jìn)行求解。初值生成器中取LGL節(jié)點(diǎn)數(shù)目為5，則離散化之后得到最終優(yōu)化變量數(shù)目為36，綜合考慮各種約束的總約束數(shù)目為60。將優(yōu)化結(jié)果作為初值進(jìn)行多節(jié)點(diǎn)優(yōu)化，選取多節(jié)點(diǎn)數(shù)目為30，則最終優(yōu)化變量數(shù)目為211，總約束數(shù)目為310。

優(yōu)化需要的目標(biāo)函數(shù)及約束的梯度函數(shù)(即雅可比陣)采用解析方法提供。最終，優(yōu)化指標(biāo)達(dá)到全局最小值為972.93 s，狀態(tài)量及控制量優(yōu)化結(jié)果如圖1～圖3所示，各路徑約束的約束變量曲線如圖4所示。

圖1 位置變化曲線Fig.1 Curves of position vs time

圖2 速度相關(guān)狀態(tài)量變化曲線Fig.2 Curves of state variables correlated with velocity vs time

圖3 最優(yōu)控制傾側(cè)角變化曲線Fig.3 Curves of optimal bank angle vs time

圖4 內(nèi)點(diǎn)路徑約束曲線Fig.4 Curves of interior point constraints

優(yōu)化計(jì)算在CPU為E7500/2.93 GHz，2G內(nèi)存的臺(tái)式計(jì)算機(jī)上實(shí)現(xiàn)，操作系統(tǒng)為Windows XP。采用MATLAB7.7.0(R2008b)語言進(jìn)行代碼編寫，MATLAB環(huán)境下進(jìn)行仿真，SQP算法實(shí)現(xiàn)參數(shù)優(yōu)化計(jì)算。算例中再入彈道優(yōu)化的時(shí)間包括初值生成時(shí)間和精確計(jì)算時(shí)間兩部分，總耗時(shí)為3.72 s;直接采用3.2節(jié)中的LGL偽譜法，優(yōu)化計(jì)算初值根據(jù)再入初始條件和終端條件線性插值簡單給出，LGL點(diǎn)個(gè)數(shù)取30，優(yōu)化計(jì)算所需時(shí)間為58.64 s，說明設(shè)計(jì)的快速優(yōu)化策略大大提高了優(yōu)化計(jì)算速度。另外，相對(duì)于傳統(tǒng)優(yōu)化算法，由于LGL偽譜法不包含積分計(jì)算，且離散得到的NLP問題具有典型的稀疏結(jié)構(gòu)。所以，從這兩點(diǎn)也可說明這種方法具有較高的計(jì)算效率。

另外，相關(guān)文獻(xiàn)研究指出，若在LINUX環(huán)境下，采用C或Fortran語言編寫代碼，并將代碼進(jìn)一步優(yōu)化，可使得計(jì)算速度提高至少100倍［15］?？梢?，在計(jì)算速度上，利用LGL偽譜法求解助推-滑翔導(dǎo)彈再入彈道優(yōu)化問題，計(jì)算效率較高。

由圖1～圖4的仿真結(jié)果可看出，所得結(jié)果滿足所有端點(diǎn)約束條件，同時(shí)滿足各種路徑約束，本文建立的優(yōu)化方法能以較少的節(jié)點(diǎn)和計(jì)算時(shí)間獲得狀態(tài)量的優(yōu)化解和優(yōu)化的控制量。

進(jìn)一步進(jìn)行可行性驗(yàn)證。將得到的最優(yōu)控制量σ*(t)進(jìn)行插值，得到整個(gè)時(shí)間區(qū)間上的σinterp1(t)，采用MATLAB的ode45龍哥庫塔積分器，積分推進(jìn)運(yùn)動(dòng)方程可得到狀態(tài)量推進(jìn)軌跡，將此結(jié)果與本文得到的優(yōu)化狀態(tài)量進(jìn)行對(duì)比，得到圖5～圖7所示結(jié)果。

由圖5～圖7的對(duì)比結(jié)果可看出，優(yōu)化結(jié)果與推進(jìn)結(jié)果吻合較好，雖然也存在一定的誤差，但對(duì)于本文所研究的大范圍的變軌問題而言，其相對(duì)誤差(定義為絕對(duì)誤差/變量值)是較小的，數(shù)量級(jí)在10-4～10-5左右，說明計(jì)算精度較高。即在理想的飛行條件下，臨近空間助推-滑翔導(dǎo)彈按照LGL偽譜法獲得的最優(yōu)制導(dǎo)指令進(jìn)行飛行，可從預(yù)定初始點(diǎn)到目標(biāo)點(diǎn)并實(shí)現(xiàn)交接班，從而驗(yàn)證了極值解的可行性。

圖5 位置比較Fig.5 Comparison of positions

圖6 速度及速度方位角比較Fig.6 Comparison of velocity and velocity azimuth

進(jìn)一步分析偏差存在的原因主要有:

(1)偽譜法離散所選取的LGL節(jié)點(diǎn)數(shù)目為30，這個(gè)點(diǎn)的數(shù)目是較少的，要進(jìn)一步提高精度，必須增加節(jié)點(diǎn)數(shù);

(2)由于本文選擇傾側(cè)角而不是傾側(cè)角速率作為控制量，這導(dǎo)致無法對(duì)傾側(cè)角變化率進(jìn)行限制，致使傾側(cè)角的變化較為劇烈，而本文采用樣條曲線對(duì)控制量進(jìn)行擬合(見圖8)，插值所得控制量并不能與優(yōu)化結(jié)果很好地吻合，致使?fàn)顟B(tài)量推進(jìn)結(jié)果并不能完全真正地反映最優(yōu)控制量的作用。

接下來進(jìn)行最優(yōu)性驗(yàn)證。本文利用解析法驗(yàn)證問題的一階最優(yōu)性，驗(yàn)證原則:Hamiltonian函數(shù)沿最優(yōu)軌線保持為零。采用CMP對(duì)協(xié)狀態(tài)進(jìn)行估計(jì)，最終得到該問題的Hamiltonian，如圖9所示?？梢姡ㄟ^CMP及解析驗(yàn)證法得到的Hamiltonian在零附近震蕩，滿足一階最優(yōu)性必要條件。

圖7 彈道傾角比較Fig.7 Comparison of trajectory obliquity

圖8 控制量優(yōu)化結(jié)果與插值結(jié)果的比較Fig.8 Comparison of optimal result and interpolation result of the controller

圖9 最優(yōu)控制的HamiltonianFig.9 Hamiltonian of optimal control

同時(shí)，給出控制量傾側(cè)角對(duì)應(yīng)的拉格朗日乘子變化曲線，如圖10所示。可知，乘子μσ滿足KKT定理。因此，極值解的最優(yōu)性得到了驗(yàn)證。

圖10 傾側(cè)角與對(duì)應(yīng)的KKT乘子曲線Fig.10 Curves of bank angle and KKT multiplier

5 結(jié)論

(1)針對(duì)臨近空間助推-滑翔導(dǎo)彈再入彈道快速優(yōu)化問題，建立了多參數(shù)、多約束的再入彈道優(yōu)化模型。

(2)利用LGL偽譜法和SQP相結(jié)合的方法詳細(xì)建立了再入最優(yōu)飛行彈道的求解步驟，針對(duì)直接利用該方法的困難，設(shè)計(jì)了含初值生成器串行優(yōu)化求解策略，有效解決了助推-滑翔導(dǎo)彈再入段彈道快速優(yōu)化問題。仿真結(jié)果表明，對(duì)于再入彈道的優(yōu)化求解，LGL偽譜法具有較高的計(jì)算效率和計(jì)算精度，能夠?qū)崿F(xiàn)多約束、多變量、大型稀疏的再入軌跡快速優(yōu)化。

(3)完成了極值解的有效性檢驗(yàn)，計(jì)算分析表明，優(yōu)化結(jié)果滿足極小值原理對(duì)應(yīng)的最優(yōu)性必要條件。

(4)本文的研究成果能提供實(shí)時(shí)或近實(shí)時(shí)、開環(huán)最優(yōu)解，可為在線制導(dǎo)、反饋控制的設(shè)計(jì)提供基本保證。

［1］徐瑋，孫丕忠，夏智勛.助推-滑翔導(dǎo)彈總體一體化優(yōu)化設(shè)計(jì)［J］.固體火箭技術(shù)，2008，31(4):317-320.

［2］李瑜，楊志紅，崔乃剛.洲際助推-滑翔導(dǎo)彈全程突防彈道優(yōu)化［J］.固體火箭技術(shù)，2010，33(2):125-130.

［3］陳小慶，侯中喜，劉建霞.高超聲速滑翔式飛行器再入軌跡多目標(biāo)多約束優(yōu)化［J］.國防科技大學(xué)學(xué)報(bào)，2009，31(6):77-83.

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

［5］Gergaud J，Haberkorn T.Orbital transfer:some links between the low thrust and the impulse cases［J］.Acta Astronautica，2007，60:649-657.

［6］Hui T L，Ping Y J，Qun F，et al.Reentry skipping trajectory optimization usingdirectparameteroptimization method［C］//14th AIAA/AHI Space Planes and Hypersonic Systems and Technologies Conference，2006.

［7］Gao Y，Li W Q.Systematic direct approach for optimizing continuous-thrust earth-orbit transfers［J］.Journal of Aeronautics，2009，22:56-69.

［8］Desai P N，Conway B A.Six-degree-of-freedom trajectory optimization using a two-timescale collocation architecture［J］.Journal of Guidance，Control，and Dynamics，2008，31(5):1308-1315.

［9］Geoffrey T Huntington，David Benson，Anil V Rao.A comparison of accuracy and computational efficiency of three pseudospectral methods［R］.AIAA 2007-6405.

［10］Bemon A，Thorvaldsen T，Rao V.Direct tmjectory optimization and costae estimation via an orthogonal collocation method［J］.Journal of Guidance，Control，and Dynamics，2006，29(6):1435-1440.

［11］楊希祥，張為華.基于Gauss偽譜法的固體運(yùn)載火箭上升段軌跡快速優(yōu)化研究［J］.宇航學(xué)報(bào)，2011，32(1):15-21.

［12］孫勇，段廣仁，張卯瑞，等.高超聲速飛行器再入過程改進(jìn)氣動(dòng)系統(tǒng)模型［J］.系統(tǒng)工程與電子技術(shù)，2011，33(1):134-137.

［13］Etkin B，Reid L D.Dynamics of flight-stability and control［M］.New York:John Wiley and Sons lnc，1995.

［14］Philip E G.SNOPT:an SQP algorithm for large-scale constrained optimization［J］.SIAM Review，2005，47(1):99-131.

［15］閆循良.基于空間發(fā)射的組合機(jī)動(dòng)路徑規(guī)劃研究［D］.西安:西北工業(yè)大學(xué)，2011.

Rapid re-entry trajectory optimization for boost-glide missile

ZHAO Xin1，YAN Xun-liang1，ZHANG Jin-sheng1，WANG Shi-cheng1，HE An-rong2
(1.The Second Artillery Engineering University，Xi＇an 710025，China;2.The Military Deputy office of the Second Artillery in 211 Factory，Beijing 100084，China)

Rapid re-entry trajectory optimization of near space boost-glide missile was studied via the Legendre-Gauss-Lobatto(LGL)pseudospectral method.Firstly，an accurate mathematics model in re-entry phase was established based on an improved aerodynamic model.Aiming at the difficulties of the optimization problem in processing of aerodynamic data as well as optimization solving，the steps based on LGL pseudospectral method were listed systemically for solving the optimal re-entry flight trajectory.Then a serial strategy based on LGL pseudospectral method was presented to deal with the difficulties in optimization using the LGL pseudospectral method directly.Finally，the feasibility and the optimality of optimal result were validated by integral propagation method and covector mapping principle，respectively.Simulation results show that computational time required for optimizing one reentry trajectory is 3 ～4 s，and thus the computational efficiency is high.Path constrains and boundary constrains are well satisfied，and the precision of this trajectory optimization method is high.The rapid re-entry trajectory optimization with characteristics of multiple constraints，multiple variables and large sparsity is achieved.

near space;boost-glide missile;re-entry trajectory optimization;pseudospectral method;optimal control

V448.2

A

1006-2793(2012)04-0427-07

2012-04-09;

2012-07-09。

國家自然科學(xué)基金項(xiàng)目(60874093);863項(xiàng)目(2011AA7053016)。

趙欣(1984—)，男，博士，研究方向?yàn)轱w行動(dòng)力學(xué)及控制。E-mail:zhaoxin20062111@163.com

(編輯:崔賢彬)