基于Gauss伪谱方法的高超声速飞行器再入轨迹快速优化

第29卷第6期2008年11月　
宇　航　学　报
Journal of Astronautics
V ol.29N ovember
N o.6
2008
基于G auss 伪谱方法的高超声速飞行器
再入轨迹快速优化
雍恩米,唐国金,陈　磊
(国防科技大学航天与材料工程学院,长沙410073)
摘　要:基于一种求解最优控制问题的新方法—G auss 伪谱法(G auss Pseudospectral Method 2G P M ),研究了高超声速飞行器滑翔式再入的快速轨迹优化问题。

针对远程多约束条件下滑翔式再入轨迹优化问题的难点,提出了基于G P M 的串行分段优化策略,包括三个方面:(1)构造了设计变量初值生成器,获得近似最优解作为优化初值;(2)提出从可行解到最优解的串行优化策略;(3)引入平衡滑翔条件构造动态分段点,将再入轨迹分为初始下降段和滑翔段分别求解。

以某高超声速再入飞行器为对象进行轨迹优化计算,仿真结果验证了本文的轨迹优化方法具有较高的精度和计算效率。

关键词:轨迹优化;再入;G auss 伪谱方法
中图分类号:V412.44 文献标识码:A 文章编号:100021328(2008)0621766207
DOI :10.3873Πj.issn.100021328.2008.06.016
收稿日期:2007212225;　修回日期:2008205212
0　引言
近年来,由于远程快速物资运送或精确打击等需求,各类空天飞行的高超声速飞行器成为各国研究热点。

本文研究一类从轨道或亚轨道高度再入,
具有较大升阻比(大于一般的可重复使用运载器),依靠气动力控制实现远距离滑翔的再入飞行器的轨迹优化问题。

其中具有代表性的是2003年美国国防部和空军联合推出的“猎鹰”计划中的通用航空飞
行器(C omm on Aero Vehicle 2C AV )[1]。

此类飞行器的轨迹优化问题具有如下特点:(1)再入轨迹对控制变量高度敏感;
(2)再入过程的气动热、动压和过载等约束使
得再入轨迹的可行域限制在较为狭窄的范围内。

由此可见,高超声速飞行器再入轨迹数值优化存在一定困难。

间接法和直接法是目前求解轨迹优化问题的主流方法
[2]。

间接法将最优控制问题最终
转换为一个两点边值问题,但求解过程对共轭变量的初值高度敏感且难以估计,很难收敛。

同时再入过程中存在的各种约束,使间接法求解的推导过程更为复杂
[3,4]。

相对间接法,直接法在收敛的鲁棒性
和解决实际问题的适应性上具有优势[5]。

仅离散控
制变量和同时离散控制与状态变量是直接法的两种基本类型。

由于对控制变量的高度敏感,前一种方法用于本文的问题容易陷入局部解,甚至收敛到不可行解。

后一种方法又称为直接配点法,它将节点状态变量也作为设计变量
[6]。

一般的直接配点法采
用分段多项式来近似状态和控制变量时间历程,该方法的设计变量数目庞大,且不易给定恰当的初值。

最近发展起来的一种求解最优控制问题的G auss 伪谱方法是一种基于全局插值多项式的直接配点法,它相对于一般直接配点法的优势是可以用较少的节点获得较高的精度
[7-9]。

为实现远程多约束条件下高超声速飞行器再入轨迹快速优化,本文利用G auss 伪谱方法的优点,提出含初值估计的串行并分段的优化策略,即:
(1)构造设计变量初值生成器(Initial G uess G enerator 2IGG ),先以较少节点的G auss 伪谱方法计算近似最优解,通过插值获得更多节点对应的设计变量初值,再求多节点下的更精确的解;
(2)在求解初值时,采用了从可行解到最优解的串行优化策略;
(3)将再入轨迹分为初始下降段和滑翔段分别求解。

1　再入轨迹优化问题
1.1　无动力再入运动模型
(1)无量纲再入运动方程
考虑地球为旋转圆球时,高超声速飞行器无动力再入在极坐标系下的无量纲运动方程为[10]
d r
d t
=V sinγ(1a)
dθd t =V
cosγsinψ
r cosφ
(1b)
dφd t =
V cosγcosψ
r
(1c)
d V d t =-珚D-
sinγ
r2
+Ω2r cosφ·
(sinγcosφ-cosγsinφcosψ)(1d)
dγd t =
1
V
[珔L cosσ+
cosγ
r
(V2-1
r
)+
2ΩV cosφsinψ+Ω2r cosφ·
(cosγcosφ+sinγcosψsinφ)](1e)
dψd t =
1
V
珔L sinσ
cosγ
+V
2
r
cosγsinψtanφ-
2ΩV(tanγcosφcosψ-sinφ)+
Ω2r
cosγ
sinφcosφsinψ(1f)
其中无量纲地心距r、速度V、时间t和地球自转角速度Ω的无量纲化参数分别为地球平均半径R0、V c=g0R0、R0/g0和g0/R0,g0为海平面引力加速度。

θ、φ、γ、ψ分别为经度、纬度、航迹角和航向角。

而无量纲升力、阻力加速度珔L、珚D分别为:
珔L=ρ(VV
c )2SC
L
/(2mg0)
珚D=ρ(VV
c )2SC
D
/(2mg0)
其中ρ为大气密度,m、S分别为再入飞行器质量和气动参考面积。

(2)气动力模型
在高超声速飞行条件下,升力、阻力系数C
L、C D可近似表示为攻角的函数:
C L=C Lα·α
C D=C D0+CαD·α2
(2)
升力系数Cα
L
,阻力系数C D0、CαD由气动力数据
表拟合而成。

(3)再入攻角模型
对于大升阻比再入飞行器,热防护是攻角设计考虑的主要因素。

为结合再入制导问题,并充分利用已有研究成果,不直接设计离散点的攻角控制变量,而是事先确定攻角曲线。

由于高超声速再入飞行器的防热问题很严重,其再入初始阶段一般采用最大攻角飞行,才能满足热流约束。

而后减小攻角,并以最升阻比攻角飞行,以增大航程。

据此,本文设计攻角曲线为速度的分段线性函数,即:
α=
α
m ax
V1≤V<V e
a m ax L/D-αm ax
V2-V1
(V-V
1
)+αm ax V2≤V<V1
α
m ax L/D
V f≤V<V2
(3)
其中α
max
、α
max L/D
分别为最大飞行攻角和最大升阻
比对应的攻角。

V
e、V f为再入点速度和终端速度。

参数V
1、V2由航程要求决定。

1.2　约束条件
(1)弹道约束
为保证再入飞行器在结构和热防护上的可靠性,再入过程要求严格满足热流密度、动压和过载约束,即
Q=C1
R
ρ
ρ
s
n
V m≤
Q max(4) q=
1
2
ρ(VV
c
)2≤q
max
(5)
n=珔L2+珚D2=q C2D+C2L
S
mg
≤n
max
(6) 对高超声速飞行器,一般取n=0.5,m=3或3.15,其中C1为常数,R为飞行器头部曲率半径,
ρ
s为海平面大气密度,V c为环绕速度。

(2)终端约束
考虑到与末制导段交班的要求,再入终端状态应满足一定的约束。

终端状态等式约束一般为包括经纬度和地心距的位置约束。

另外,某些飞行器对终端航向角和速度大小也有一定要求。

约束条件见(7)、(8)式。

θ(t
f
)=θ
f
,φ(t f)=φf,r(t f)=r f(7) |Δψ(t f)|≤Δψf,V(t f)≥V f(8) 1.3　目标函数
对高超声速再入飞行器来说,防热问题是系统设计的难点之一。

再入轨迹优化的一个重要目标就是减小气动加热,以降低防热系统的重量,增加有效载荷的比重。

对驻点热流密度积分,得到再入过程的热载作为优化目标:
Q s=∫t f t0 Q d t(9)
7671
第6期雍恩米等:基于G auss伪谱方法的高超声速飞行器再入轨迹快速优化
另外,考虑到弹道平滑性,引入航迹角相关指标:
I=c·∫t f t0 γ2d t(10) 2　G auss伪谱方法
2.1　基本原理
G auss伪谱方法将状态变量和控制变量在一系列G auss点上离散,并以这些离散点为节点构造La2 grange插值多项式来逼近状态和控制变量。

通过对全局插值多项式求导来近似状态变量对时间的导数,从而将微分方程约束转换为一组代数约束。

性能指标中的积分项由G auss积分计算。

终端状态也由初始状态和对右函数的积分获得。

经上述变换,可将最优控制问题转化为具有一系列代数约束的参数优化问题,即非线性规划问题(N LP)。

2.2　连续最优控制问题的离散
(1)时域变换
采用G auss伪谱方法需要将轨迹优化问题的时间区间[t
,t f]转换到[-1,1],对时间变量t作变换:
τ=2t
t f-t0-
t f+t0
t f-t0
(11)
(2)全局插值多项式近似状态与控制变量
G auss伪谱法取K阶Legendre2G auss(LG)点以及τ
=-1作为节点,构成K+1个Lagrange插值多项式,并以此为基函数构造状态变量的近似表达式,即:
x(τ)≈X(τ)=6K i=0L i(τ)x(τi)(12)其中Lagrange插值基函数
L i(τ)=7K j=0,j≠iτ-τj
τ
i
-τj
且使节点处的近似状态与实际状态相等,即X
i
= X(τi)=x(τi)。

G auss伪谱方法中未用到控制变量的导数,因此任何满足性质
u(τk)=U k,(k=1,…,K)
的近似对一个N LP问题来说都是等价的。

为了形式上的统一,这里仍以LG点作为节点,以K个La2 grange插值多项式珘L i(τ)(i=1,…,K)为基函数来近似控制变量,即:
u(τ)≈U(τ)=6K i=1珘L i(τ)U(τi)(13) (3)动力学微分方程约束
对(12)式求微分有
x(τk)≈ X(τk)=6K i=0D ki(τk)X(τi)(14)其中微分矩阵D∈R K×(K+1)可以离线确定,即:
D ki=
L i(τk)
=
(1+τ
k
) P
K
(τ
k
)+P
K
(τ
k
)
(τ
k
-τi)[(1+τi) P K(τi)+P K(τi)]
　i≠k
(1+τ
i
)¨P
K
(τ
i
)+2 P
K
(τ
i
)
2[(1+τi) P K(τi)+P K(τi)]
i=k
(15)
从而动力学微分方程约束转换为代数约束:
6K i=0D ki X(τi)-t f-t02·
f(X(τk),U(τk),τk;t0,t f)=0(16)其中k=1,…,K。

(4)终端状态约束
状态变量的近似表达式未包括终端时刻节点,终端状态应满足动力学方程约束,即:
　x(τ
f
)=x(τ0)+∫1-1f(x(τ),u(τ),τ)dτ(17) 将终端约束条件离散并用G auss积分来近似,可得:
X(τf)=X(τ0)+
t f-t0
26w k f(X(τk),
U(τk),τ,t0,t f)(18) (5)性能指标函数的近似
将Bolza型性能指标函数中的积分项用G auss 积分来近似,得到G auss伪谱方法中的近似性能指标函数:
J=Φ(X0,t0,X f,t f)+
t f-t0
26
2
k=1
w k g(X k,U k,τk;t0,t f)(19) 根据上述数学变换,基于G auss伪谱法的最优
控制问题可描述为:求离散节点上的状态X
i
(i=0,
…,K)和控制变量U
k
(k=1,…,K),以及初末时
刻t
0、t f(若未给定),使得性能指标(19)最小,并满足动力学方程约束(16),终端状态约束(17),以及原最优控制问题的边界条件
φ(X
,t0,X f,t f)=0(20)和过程约束
C(X k,U K,τk;t0,t f)≤0(k=1,…,K)(21) 从而将原最优控制问题转化成一个一般非线性规划问题,即:
min
y∈R M
　F(y)
s.t.　g j(y)≥0,　j=1,2,…,p
8671宇航学报第29卷
h j (y )=0,　j =1,2,…,l (22)
其中y 为包含状态、控制变量和端点时间的设计变量。

3　求解优化问题的策略及流程
3.1　基于G PM 求轨迹优化问题的策略分析
从理论上分析,直接应用2.2节的G auss 伪谱法可以求解本文第1节的再入轨迹优化问题,然而在实际应用中存在以下困难:(1)当G auss 节点较多时,设计变量的数目会比较庞大。

设LG 点的个数为K ,则设计变量数目为6(K +1)+K +1,当阶次大于8时,设计变量就会超过50个。

那么,对新型飞行器的轨迹优化问题,
设定设计变量的初值的工作会比较繁杂,且不恰当的初值可能使问题收敛到不可行解;(2)G auss 伪谱法将动力学微分方程约束转换为代数约束,相应的优化问题中等式约束为6N 个,同时再入轨迹优化问题还存在较为苛刻的不等式约束和终端状态约束,在如此多约束条件下,采用一般的优化算法很难直接有效地找到可行解。

另外,大升阻比升力体飞行器再入时在初始下降段和滑翔段呈现出不同的特性,初始下降段轨迹主要由再入条件决定,控制变量对该段轨迹影响较小,滑翔段一般保持拟平衡滑翔状态以有效利用升力增大航程。

充分理解再入轨迹的特点,可为轨迹优化问题带来方便。

根据以上高超声速再入轨迹优化存在问题和特点,本文提出以下基于G auss 伪谱方法求解最优轨迹的策略:
(1)构造设计变量初值生成器(IGG )利用G auss 伪谱法能以较少节点获得较高精度的优势,先采用较少的节点,例如K =5,计算近似的最优轨迹和控制变量,并以此为下一步精确计算的初值。

(2)采用从可行解到最优解的串行优化策略设原问题的等式约束同(22)式,串行优化策略则是指,不直接寻找满足所有等式和不等式约束的解,而是先将等式约束转换为目标函数,即
min F (y )=sqrt (
6
l
i =1
h i (y )2
)
s.t.　g j (y )≥0,　j =1,2,…,p
(23)
获得的最优解y 0。

再以y 0为初值,求解原最优控制问题(22)的解。

(3)将再入轨迹分为初始下降段和拟平衡滑翔段分别优化,分段点动态选取,即:
取再入轨迹在高度-速度(r -V )平面的斜率
d r /d V 满足平衡滑翔条件以及最大热流密度约束(再入初期一般动压和过载较小,可以满足约束,因此不考虑)时刻为分段转折点,即满足条件
d r d V
3DOF
-
d r d V
QEGC
<δ
Q Trans ≤ Q max
(24)其中
d r d V
3DOF
可由运动方程得到近似表达式
d r
d V
=V sin γ
-D -sin γ/r
2
求解。

d r
d V
QEGC
为平衡滑翔轨迹在
高度-速度剖面对应的斜率,给定一个速度V ,可由平衡滑翔条件求出对应r ,再由有限差分求解d r d V
QEGC。

3.2　优化流程
本文采用G PM 作为每阶段优化的基本方法。

根据上述优化策略,优化流程分两步,首先由IGG 获得初值,并采用从可行解到最优解的策略。

第二步选取足够多的节点(K =10～20),以求得具有一定精度的最优轨迹。

同时为充分利用再入轨迹不同阶段特性,将再入轨迹分为初始下降段和滑翔段分别优化。

图1描述了本文提出的分段串行优化策略的基本流程。

图1　基于G auss 伪谱方法的再入轨迹优化策略
Fig.1　The Reentry trajectory optimization strategy via G P M
9
671第6期雍恩米等:基于G auss 伪谱方法的高超声速飞行器再入轨迹快速优化
4　数值计算实例
以远程高超声速滑翔式再入飞行器为仿真对象,气动参数采用美国波音公司98年设计的带控制翼的锥形体再入机动飞行器C AV 的相关参数拟合
得到
[11],其最大升阻比为2.4,气动参考面积为
0.35m 2
,质量为907kg ,最大飞行攻角和最大升阻比
攻角分别为20°和10°。

再入初始弹道参数与热流密度、动压、过载和侧倾角变化率约束见表1。

表1　初末条件和过程约束
T able 1　Initial conditions and path constraints
V (m/s )
h (km )
γ(°)θ(°)φ(°)ψ(°)初始条件7200100-2
1602558
终端条件≥1000
20
258
38
Q max (K
W/m 2)q max (kPa )
n max
σmax (deg/s )
过程约束
1000
5006
45
优化目标为再入过程热载最小,并同时考虑但
到平滑的因素,取综合性能指标:
J =
∫
t
f
t
Q d t +c ·∫
t
f
t
γ2d t
(25)
初值生成时取G auss 节点个数K =5,求最优
轨迹时初始下降段和滑翔段分别取K =5和K =12。

优化计算在CPU 为2.8G H z/Pentimu
I V ,操作系统为
Windows XP
的微机上实现,Matlab 环境下编程,S QP 算法实现参数优化计算。

算例中再入轨迹
优化的时间开销约为4-5分钟左右。

根据笔者的
经验[12]
,如果在C ++环境下编写同样数学模型的程序,运行时间应在10秒以内,相较于文献[13]的结果有所改进。

在Matlab 环境下,一般直接法求解一个二维的再入轨迹优化问题需要的机时约为30分钟[14]。

可见,G auss 伪谱方法相对一般直接法在计算速度上有很大优势。

优化计算结果见图2。

主要轨迹参数给出了四种情况进行比较,“G PM ”为本文提出策略下最终
771宇航学报第29卷
图2　最优轨迹计算结果Fig.2　Results of optimal trajectory
优化结果,G auss节点的轨迹参数是获得的设计变量值,近似解由IGG得到的初值参考轨迹,数值积分解是将伪谱方法得到的控制变量带入原运动方程积分得到的轨迹。

结果表明,近似解已经比较接近最优解,可以为进一步优化提供有效初值。

本文的G PM 得到的最优解与数值积分结果基本一致,具有较高的精度。

最优控制变量变化平缓,见图2(d)。

最优轨迹满足各种弹道约束,如图2(e)(f)所示。

5　结论
本文基于G auss伪谱法求解高超声速无动力再入飞行器的快速轨迹优化问题。

针对直接应用该方法的困难,提出了包括设计变量初值生成器、从可行解到最优解的串行优化以及将轨迹分段的优化的策略。

仿真计算结果表明,相较于一般的直接法,在该优化策下应用G auss伪谱法可以快速有效地求解高超声速飞行器多约束条件下的再入最优轨迹。

本文研究的成果可进一步应用到再入制导问题中。

参考文献:
[1]　G eorge R.The comm on aero vehicle2space delivery system of the fu2
ture the comm on aero vehicle2space delivery system of the future[C]ΠΠ
AIAA S pace T echnology C on ference&Exposition,Albuquerque,
NM,1999:1-11.
[2]　Betts J T.Survey of numerical methods for trajectory optim ization
[J].Journal of G uidance,C ontrol and Dynam ics,1998,21(2):
193-206.
[3]　Istratie V.Optimal profound entry into atm osphere with m inimum heat
and constraints[C]ΠΠAIAA Atm ospheric Flight M echanics C on ference
and Exhibit.M ontreal,Canada,2001.
[4]　E lnagar J,K azem i M A.Pseudospectral chebyshev optimal control of
constrained nonlinear dynam icla systems[J].C om putational Optim i2
zation and Applications,1998(11):195-217.
[5]　Hui T L,Ping YJ,Qun F,et al.Reentry skipping trajectory opti2
m ization using direct parameter optim ization method[C]ΠΠ14th AIAA/
AHI S pace Planes and Hypers onic Systems and T echnologies C on fer2
ence,2006.
[6]　R obert T.T rajectory optim ization for a fixed2trim reentry vehicle us2
ing direct collocation and nonlinear programm ing[C]ΠΠAIAA G uid2
ance,Navigation,and C ontrol C on ference&Exhibit.Denver,US A,
2000.
[7]　Huntington G T.Advancement and analysis of a gauss pseudospectral
transcription[D].Cambridge,M A:M assachusetts Institute of techn2
ology,2007.
[8]　Bens on A,Thorvaldsen T,Rao V.Direct trajectory optim ization and
costae estimation via an orthog onal collocation method[J].Journal of
G uidance,C ontrol and Dynam ics,2006,29(6):1435-1440.
[9]　罗建军,王明光,袁建平.基于伪光谱方法的月球软着陆轨道
快速优化[J].宇航学报,2007,28(5):1119-1122.[LUO Jian2
jun,W ANG M ing2guang,Y UAN Jian2ping.The legendre
pseudospectral method for rapid lunar s oft2landing trajectory optim iza2
tion[J].Journal of Astronautics,2007,28(5):1119-1122(in Chi2
nese).]
[10]　Vinh N X,Busemann A,Culp R D.Hypers onic and Planetary Entry
Flight M echanics[M].Ann Arbor,MI:Univ.of M ichigan Press,
1980:26-27.
[11]　C orporation T P.A comm on aero vehicle m odel,description,and
em ployment guide[E B/O L].w w w.dtic.m il/matris/sbir/sbir041/
srch/af031a.doc,2003
[12]　雍恩米,唐国金,陈磊.高超声速无动力远程滑翔飞行器多
约束条件下的轨迹快速生成[J].宇航学报,2008,29(1):46-
52.[Y ONG En2m i,T ANG G uo2jin,Chen Lei.Rapid trajectory
planning for hypers onic unpowered long2range reentry vehicles with
1771
第6期雍恩米等:基于G auss伪谱方法的高超声速飞行器再入轨迹快速优化
multi2constraints[J].Journal of Astronautics,2008,29(1):46-52
(in Chinese).]
[13]　王明光,袁建平,罗建军.R LV再入轨迹机载快速优化[J].宇
航学报,2005,26(3):253-256.[W ANGM ing2guang,Y UAN Jian2
ping,LUO Jian2jun.Onbord optim ization of three2dimensional con2
strained reentry trajectory for R LV[J].Journal of Astronautics,
2005,26(3):253-256(in Chinese).]
[14]　雍恩米,唐国金,陈磊.助推-滑翔式导弹中段弹道方案的初
步分析[J].国防科技大学学报,2006,28(6):6-10.[Y ONG
En2m i,T ANG G uo2jin,Chen Lei.Schematic study of m id2course
trajectories for the boost2glide m issiles[J].Journal of National Uni2versity of Defense T echnology,2006,28(6):6-10(in Chinese).
]
作者简介:雍恩米(1979-),女,博士生,飞
行器设计专业,研究方向为飞行器再入轨
迹优化与制导。

通信地址:湖南长沙国防科技大学航天与
材料工程学院108教研室(410073)
电话:(0731)4576009
E2mail:w fynm@
R apid Trajectory Optimization for H ypersonic R eentry V ehicle
Via G auss Pseudospectral Method
Y ONG En2mi,T ANG G uo2jin,CHE N Lei
(C ollege of Aerospace and M aterial Engineering,National University of Defense T echnology,Changsha410073,China)
Abstract:T rajectory optimization of a hypers onic glide2reentry vehicle is studied via a new optimal control method2G auss Pseudospectral Method(G P M).A pipeling2segmenting strategy is proposed to deal with the difficulties in optimization of long2rang trajectory with several path constraints.The strategy includes:(1)The Initial G uess G enerator(IGG)is constructed to present a property initial value with a near optimal s olution;(2)A pipeling optimization framew ork is adopted to obtain the optimal s olution from a feasible s olution;(3)The whole reentry trajectory is dynamically divided into initial descending phase and glide phase by im posing the equilibrium glide condition.The precision and efficiency of this trajectory optimization method are dem onstrated by applying it to a hypers onic reentry vehicle.
K ey w ords:T rajectory optimization;Reentry;G auss pseudospectral method
(上接第1705页)
Numerical Simulation on C avitation of Submarine Launched
Missile’s Surface at Large Angles of Attack
QUAN X iao2bo,WEI Hai2peng,K ONG De2cai,LI Y an
(Beijing Institute of Astronautical System Engineering,Beijng100076,China)
Abstract:Cavitation on the missile’s sur face has great effect on the forces and mechanical environment of submarine launched missile m oving underwater with a high speed.The numerical com putations for cavitation at large angles of attack are made.And its precision is dem onstrated with the com paris ons between the calculation results and experiments.The distribution of pressure on sur face and the shape of cavitation at different attack angles is obtained.The in fluence of attack angle and cavitation number on the hydrodynamic force is given.
K ey w ords:Submarine launched missile;Cavitation flows;Numerical simulation
2771宇航学报第29卷。