CN107861517A - 基于线性伪谱的跳跃式再入飞行器在线弹道规划制导方法 - Google Patents
基于线性伪谱的跳跃式再入飞行器在线弹道规划制导方法 Download PDFInfo
- Publication number
- CN107861517A CN107861517A CN201711058719.0A CN201711058719A CN107861517A CN 107861517 A CN107861517 A CN 107861517A CN 201711058719 A CN201711058719 A CN 201711058719A CN 107861517 A CN107861517 A CN 107861517A
- Authority
- CN
- China
- Prior art keywords
- mrow
- msub
- mtd
- mtr
- msubsup
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 28
- 238000001228 spectrum Methods 0.000 title claims abstract description 14
- 238000012937 correction Methods 0.000 claims abstract description 29
- 230000010354 integration Effects 0.000 claims abstract description 22
- 230000010006 flight Effects 0.000 claims abstract description 9
- 230000011218 segmentation Effects 0.000 claims abstract description 8
- 239000011159 matrix material Substances 0.000 claims description 34
- 230000006870 function Effects 0.000 claims description 15
- 238000000926 separation method Methods 0.000 claims description 8
- 230000033001 locomotion Effects 0.000 claims description 5
- 238000005094 computer simulation Methods 0.000 claims description 3
- 238000012546 transfer Methods 0.000 claims description 3
- 230000000149 penetrating effect Effects 0.000 claims description 2
- 206010048669 Terminal state Diseases 0.000 claims 1
- 235000013399 edible fruits Nutrition 0.000 claims 1
- 238000004364 calculation method Methods 0.000 abstract description 7
- 230000001360 synchronised effect Effects 0.000 abstract description 2
- 241000218228 Humulus Species 0.000 description 14
- 238000004422 calculation algorithm Methods 0.000 description 6
- 238000004611 spectroscopical analysis Methods 0.000 description 6
- 238000013461 design Methods 0.000 description 4
- 238000010586 diagram Methods 0.000 description 3
- 238000012986 modification Methods 0.000 description 3
- 230000004048 modification Effects 0.000 description 3
- 238000010606 normalization Methods 0.000 description 3
- 239000007787 solid Substances 0.000 description 3
- 238000013459 approach Methods 0.000 description 2
- 125000004122 cyclic group Chemical group 0.000 description 2
- 238000009795 derivation Methods 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 239000003550 marker Substances 0.000 description 2
- NCGICGYLBXGBGN-UHFFFAOYSA-N 3-morpholin-4-yl-1-oxa-3-azonia-2-azanidacyclopent-3-en-5-imine;hydrochloride Chemical compound Cl.[N-]1OC(=N)C=[N+]1N1CCOCC1 NCGICGYLBXGBGN-UHFFFAOYSA-N 0.000 description 1
- 241000208340 Araliaceae Species 0.000 description 1
- PEDCQBHIVMGVHV-UHFFFAOYSA-N Glycerine Chemical compound OCC(O)CO PEDCQBHIVMGVHV-UHFFFAOYSA-N 0.000 description 1
- 244000025221 Humulus lupulus Species 0.000 description 1
- 235000005035 Panax pseudoginseng ssp. pseudoginseng Nutrition 0.000 description 1
- 235000003140 Panax quinquefolius Nutrition 0.000 description 1
- 230000001133 acceleration Effects 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 210000000709 aorta Anatomy 0.000 description 1
- 239000002775 capsule Substances 0.000 description 1
- 230000019771 cognition Effects 0.000 description 1
- 238000012938 design process Methods 0.000 description 1
- 239000010432 diamond Substances 0.000 description 1
- 238000009826 distribution Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000010304 firing Methods 0.000 description 1
- 239000012530 fluid Substances 0.000 description 1
- 235000008434 ginseng Nutrition 0.000 description 1
- 230000005484 gravity Effects 0.000 description 1
- 230000009191 jumping Effects 0.000 description 1
- 238000013507 mapping Methods 0.000 description 1
- 230000013011 mating Effects 0.000 description 1
- 239000000203 mixture Substances 0.000 description 1
- 239000002245 particle Substances 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 230000000750 progressive effect Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G05—CONTROLLING; REGULATING
- G05D—SYSTEMS FOR CONTROLLING OR REGULATING NON-ELECTRIC VARIABLES
- G05D1/00—Control of position, course, altitude or attitude of land, water, air or space vehicles, e.g. using automatic pilots
- G05D1/10—Simultaneous control of position or course in three dimensions
- G05D1/107—Simultaneous control of position or course in three dimensions specially adapted for missiles
-
- F—MECHANICAL ENGINEERING; LIGHTING; HEATING; WEAPONS; BLASTING
- F42—AMMUNITION; BLASTING
- F42B—EXPLOSIVE CHARGES, e.g. FOR BLASTING, FIREWORKS, AMMUNITION
- F42B15/00—Self-propelled projectiles or missiles, e.g. rockets; Guided missiles
- F42B15/01—Arrangements thereon for guidance or control
Landscapes
- Engineering & Computer Science (AREA)
- Aviation & Aerospace Engineering (AREA)
- General Engineering & Computer Science (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Automation & Control Theory (AREA)
- Chemical & Material Sciences (AREA)
- Combustion & Propulsion (AREA)
- Control Of Position, Course, Altitude, Or Attitude Of Moving Bodies (AREA)
- Other Investigation Or Analysis Of Materials By Electrical Means (AREA)
Abstract
本发明提供了一种基于线性伪谱的跳跃式再入飞行器在线弹道规划制导方法，包括如下步骤：初始化，标控飞行段，跳跃段倾侧反转判断，跳跃段预测弹道积分，精度判定，跳跃段倾侧反转时刻修正，重新建立参数控制，跳跃段预测弹道积分，精度判定，跳跃段控制参数修正，Kepler轨道飞行。本发明通过迭代更新倾侧角模值参数进而修正终端射程偏差，当分段条件不由时间决定时，引入控制补偿，为末段制导提供了良好的初始条件。而且在一个制导周期内同步更新倾侧角模值控制规律以及倾侧角反转时间，进而同时确保终端横程与纵程精度，得益于控制量修正值与终端偏差之间的解析关系，每一次控制量的更新都不需要迭代积分，单次计算耗时仅为2ms，非常适合在线应用。
Description
技术领域
本发明涉及跳跃式再入飞行器制导领域，更具体的说是涉及一种基于线性伪谱的跳跃式再入飞行器在线弹道规划制导方法。
背景技术
随着嫦娥计划的不断进展，探月返回轨道的设计已经成为研究的热点。区别于一般的高超声速再入飞行器，以美国的阿波罗探月返回舱为例，这一类飞行器再入速度大、初始能量高，通过设计跳跃式再入弹道，可以在确保热流约束与过载约束的前提下，有效地降低飞行器能量，使其二次再入时能够直接进入末制导段。加之这类飞行器普遍采用低升阻比的气动设计，再入过程中对弹道偏差的修正能力较弱。所以，这类低升阻比跳跃式再入飞行器的制导律，在设计过程中将更具有挑战性。
目前，这类制导律的设计，普遍采用基于弹道积分的预测校正方法。这类方法的基本思想是，设置初始的倾侧角模值控制规律和反转时刻，通过弹道积分获得终端射程偏差和横程偏差；而后修改模值规律与反转时刻，再次通过弹道积分评判终端偏差；如此往复，直至终端偏差达到精度要求，这一过程中往往采用牛顿迭代。需要强调的是，经过往复地循环积分后，只能得到当前时刻下的控制修正。实际飞行中，由于气动参数偏差以及大气模型的偏差，制导律需要能够在线修正弹道积累误差。而循环积分对弹载计算机要求高，计算耗时长，很不利于在线应用。这是目前方法的缺点之一。
其次，无论是通过迭代积分修正模值控制规律，以使二次再入点对应的剩余射程满足终端约束；还是通过迭代积分确定偏射角度进而修正倾侧翻转时刻，以使飞行器经过Kepler惯性飞行后，能够抵消地球自转带来的横程偏差。这种预测校正方法采用类似于试凑的循环积分，没有触及问题本质，它不能够以数学形式描述终端偏差与控制修正量之间的对应关系。
因此，如何提供一种保证终端纵向射程约束，并有效控制横向落点偏差的再入制导方法是本领域技术人员亟需解决的问题。
发明内容
有鉴于此，本发明利用线性伪谱的公式推导，获得了终端脱靶量与控制量修正值之间的解析关系，区别于以往将横程与纵程分开制导的方式，本发明可以使用同一套制导方法，在一个制导周期内同步更新倾侧角模值控制规律以及倾侧角反转时间，进而同时确保终端横程与纵程精度。所以本发明不仅在理论层面统一了横程制导与纵程制导的设计思想，同时得益于解析公式的获得，该方法也非常适用于实际工程的在线应用。
为了达到上述目的，本发明采用如下技术方案：
一种基于线性伪谱的跳跃式再入飞行器在线弹道规划制导方法，包括如下步骤：
S1：初始化：设置初始、终端的计算仿真参数，包括终端纵向射程偏差δS和横程落点偏差δχ的精度要求，通过离线弹道优化与参数化处理，将倾侧角变化规律表示为倾侧反转时刻tre和倾侧角初始模值σ1的多段参数函数，并将tre和σ1作为标准控制；
S2：标控飞行段：当飞行器距地面高度大于100km时，按标准控制飞行；当高度低于100km时，进入步骤(3)；
S3：跳跃段倾侧反转判断：按标准控制飞行，根据当前时间与倾侧角反转时刻的对应关系，判断是否进行倾侧反转，如果当前时刻未达到反转时刻时，不进行反转，进入步骤S4，如果当前时刻大于或等于反转时刻，且不是最后一个反转点时，重新建立非线性参数控制问题，进入步骤S4，当作为最后一个反转点时，进入步骤S7；
S4：跳跃段预测弹道积分：使用当前时刻的状态量作为积分初值，所述的标准控制作为控制输入，通过弹道积分获得终端状态偏差δxf，以及全局多段弹道信息：状态量序列Xk、控制量序列Uk；
S5：精度判定：由终端状态偏差δxf求解终端纵向射程偏差δS和横程落点偏差δχ，当二者都满足步骤S1设置的精度要求时，保持标准控制不变，进入步骤S3，否则，进入步骤S6；
S6：跳跃段倾侧反转时刻修正：基于步骤S4获得的多段预测弹道积分，获得终端偏差与控制修正量之间的解析关系，并更新标准控制，即倾侧反转时刻tre和倾侧角初始模值σ1，
考虑到弹道积分的初值即为当前实际运动状态，所以初始状态偏差始终为零，终端状态偏差最终表示为
其中，K矩阵表示误差传递矩阵，上标数字分别表示第1段、第2段和第3段，下标x和u分别表示相应的终端状态偏差是由状态量偏差或控制量偏差导致；
二次再入点是以高度重新达到81.5km确定的，不由飞行时间决定；考虑到二次再入点的高度约束表示为
结合标准轨迹上基准位置处的导数信息，终端状态偏差修正为
δx0为初始状态偏差，hf表示由等时弹道积分获得的终端高度，即二次再入点的高度，终端状态偏差进一步表示为
其中，Y5为时间修正矩阵，用以描述终端时刻的修正关系；
同理，不仅仅第三段的终端状态是以高度81.5km作为分界点，跳跃式再入轨迹的第二段与第三段的分段条件也是以高度81.5km确定的，对终端状态偏差做出修正
其中，时间修正矩阵Y4用以描述二三段分界点处的修正关系；
进入步骤S7；
S7：重新建立参数控制：仅将倾侧角模值参数作为标准控制，用以修正纵向射程偏差，当飞行器高度再次超过81.5km，此时认为飞行器脱离大气层，不受气动力的影响，进入步骤S11，否则，进入步骤S8；
S8：跳跃段预测弹道积分：在当前标准控制的作用下，通过弹道积分获得终端状态偏差δxf，以及全局多段弹道信息：状态量序列Xk、控制量序列Uk；
S9：精度判定：根据终端状态偏差δxf求解终端射程偏差δS，当δS满足步骤S1设置的精度要求时，保持标准控制不变，进入步骤S7，否则，终端偏差过大，进入步骤S10；
S10：跳跃段控制参数修正：基于步骤S8获得的多段预测弹道积分，计算得到终端射程偏差与倾侧角模值参数修正值之间的解析关系，并进一步更新标准控制，进入步骤S7；
S11：Kepler轨道飞行：飞行器沿标准Kepler轨道飞行，直至二次再入时距地面高度小于等于81.5km，Kepler段飞行结束，此时再入飞行器完成跳跃段制导飞行，仿真结束。
优选的，在上述基于线性伪谱的跳跃式再入飞行器在线弹道规划制导方法中，所述步骤S6中：
综合考虑倾侧角初始模值σ1和倾侧角反转时刻tre的影响，终端偏差可以进一步整理为
其中，Sx2维的系数矩阵M表征终端偏差与控制量修正值之间的解析关系。
选取两组终端约束函数，用于求解控制量的修正值
y1＝cos(S)＝sinφsinΦ+cosΦcosφcos(Θ-θ)
其中，S为飞行器纵向射程，θ和φ分别表示飞行器所在位置的经度和纬度，其单位是弧度rad；Θ和Φ分别表示飞行器着陆点的经度和纬度；将飞行器质心、着陆点、地球圆心三点围成的射面与飞行器当地北向的夹角定义为视线角Ψ，ψ表示速度矢量在当地水平面的投影与正北方向的夹角，顺时针方向为正，二者的单位是弧度；对上式两侧求微分
并定义系数矩阵Z1
所以，标准控制的修正量可以解析地表示为
标准控制迭代更新为
其中，下标k表示第k次迭代。
优选的，在上述基于线性伪谱的跳跃式再入飞行器在线弹道规划制导方法中，所述步骤S10中：
飞行器处于第二段与第三段的飞行过程中，不考虑倾侧反转时刻，终端偏差表示为
其中，Sx1维的系数矩阵N为是
仅将纵向射程作为终端约束函数
y1＝cos(S)＝sinφsinΦ+cosΦcosφcos(Θ-θ)
标准控制的修正量解析地表示为
其中，定义系数矩阵Z2
迭代更新标准控制为
σ1|k+1＝σ1|k-δσ1
其中，下标k表示第k次迭代。
经由上述的技术方案可知，与现有技术相比，本发明公开提供了一种基于线性伪谱的跳跃式再入飞行器在线弹道规划制导方法，通过迭代更新倾侧角模值参数进而修正终端射程偏差，当分段条件不由时间决定时，引入了控制补偿，为末段制导提供了良好的初始条件。而且，得益于控制量修正值与终端偏差之间的解析关系，每一次控制量的更新都不需要迭代积分，单次计算耗时仅为2ms，非常适合在线应用。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据提供的附图获得其他的附图。
图1附图为本发明基于线性伪谱的跳跃式再入飞行器在线弹道规划制导方法的流程示意图；
图2附图为本发明跳跃式再入弹道轨迹分段示意图；
图3附图为本发明倾侧角模值随飞行时间的参数化曲线；
图4附图为本发明制导算法与初始偏差的对比图；
图5附图为本发明状态量对时间的变分关系图；
图6附图为本发明引入时间补偿的高度随射程变化曲线；
图7附图为本发明全程制导的高度随射程变化曲线；
图8附图为本发明全程制导的横程偏差随时间变化曲线。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
跳跃式再入轨迹如图2所示，以距地面高度81.5km为界，脱离大气层点和二次再入点可以将全程分为三部分，即跳跃段、Kepler段和末段。本发明只考虑跳跃段和Kepler段的运动状态，其中将高度100km作为启控点，100km以上时，由于不具备气动修正能力，飞行器按照标准控制飞行。实施例中只设置一个倾侧反转点，在图2中，倾侧反转时刻记为tre。
下面将结合实例对本发明作进一步的详细说明。
一种基于线性伪谱的跳跃式再入飞行器在线弹道规划制导方法，具体步骤步骤如下：
S1:初始化：设置初始、终端的计算仿真参数，包括终端纵向射程偏差δS和横程落点偏差δχ的精度要求，通过离线弹道优化与参数化处理，将倾侧角变化规律表示为倾侧反转时刻tre和倾侧角初始模值σ1的多段参数函数，并将tre和σ1作为标准控制。
即初始化建模，具体步骤如下：
1)再入动力学方程
在圆球，考虑地球自转的假设条件下，三自由度质点再入动力学方程可以表示为
其中，所有状态量都是无量纲量，所有动力学方程都是对无量纲时间τ的导数。时间的归一化常数是地球半径R0＝6378135m，地球表面的重力加速度地球引力常数μ＝3.986005×1014m3/s2；r表示飞行器质心距地球圆心的地心距离，长度单位的归一化常量是R0；θ和φ分别表示飞行器所在位置的经度和纬度，其单位是弧度rad，不需要归一化；V表示飞行器相对于地球的速度大小，速度单位的归一化常量是γ表示飞行器速度矢量与当地水平面的夹角，称为弹道倾角，ψ表示速度矢量在当地水平面的投影与正北方向的夹角，顺时针方向为正，二者的单位是弧度，不需要归一化；σ表示飞行器的倾侧角，单位弧度，不需要归一化；Ω表示地球自转角速度的无量纲量，Ω＝0.058798，角速度单位的归一化常量是
L和D是飞行器所受升力和阻力的归一化值，其表达式为
其中，空气密度的求解公式与MATLAB的aero工具包保持一致，当距地面高度h大于11km时，ρ＝0.3639·exp[(11000-h)·1.577·10-4]；m表示飞行器质量；Sref表示飞行器的参考面积；升力系数Cl和阻力系数Cd都可以表示为攻角和马赫数的函数。需要说明的是，对于Kepler段的无控飞行，只需将方程(1.1)中的气动项置零，
此时与矢量方程
理论上保持一致，其中r表示飞行器在Kepler段的位置矢量。
2)飞行器模型
美国的猎户座载人飞船是一种典型的低升阻比再入飞行器，这种飞船的底面直径为5m，参考面积Sref＝19.635m2，总重8382kg。飞行器的攻角曲线随马赫数变化，攻角在跳跃段保持在160.2度，此时的升阻比约为0.289。
3)任务分析
文中将美国爱德华兹空军基地设为飞船的预定着落点，其经纬度信息可以表示为
Θ＝242.12° Φ＝34.905° (1.3)
剩余射程sto-go，即当前位置(θ,φ)与着陆点(Θ,Φ)之间的大圆距离可以表示为
cos(sto-go)＝sinφsinΦ+cosΦcosφcos(Θ-θ) (1.4)
将飞行器质心、着陆点、地球圆心三点围成的射面与飞行器当地北向的夹角定义为视线角Ψ，顺时针为正，表达式为
则飞行器的终端横程偏差可以定义为
χ＝sin-1[sinsto-gosin(ψ-Ψ)] (1.6)
需要说明的是，以上计算过程中，全部使用弧度单位，sto-go和χ都是无量纲量。
4)控制指令参数化
由方程(1.1)可知，在攻角规律确定的情况下，飞行器的再入轨迹由倾侧角唯一确定。倾侧角模值的变化规律如图3所示，从再入飞行器跳跃段的起点到Kepler段的终点，tre表示倾侧角反转时刻，tk代表脱离大气层的时刻，即Kepler段起点，tf代表二次再入时刻，即Kepler段的终点。th作为匀速下降段和保持段的分界点，虽然在Kepler段倾侧角的变化不影响飞行轨迹，但为了保持控制量的连续，倾侧角的模值在飞行全程都被设计为时间的参数函数。
当飞行时间介于th与tf之间时，倾侧角保持常值σf＝70deg；当飞行时间小于th时，倾侧角模值可以表示为
其中，σ1表示倾侧角初始模值，t表示当前飞行时间。σ1与th的选取方式并不唯一，二者需要通过离线弹道优化，合理组合以使飞行器的纵程满足终端约束。本发明中，σ1＝83.2°，th近似地选为tk与tf的中间时刻如此既可以确保倾侧角保持段不消失以满足终端70度的约束，也可以确保飞行器在进入Kepler飞行前保持控制规律的唯一性。
综上，由倾侧反转时刻tre和倾侧角模值分段点th，倾侧角可以被分为三段，
其中，倾侧角在第二段与第三段的符号相同，并于第一段相反，与±表示倾侧角符号的改变。
至此，倾侧角已被表示为倾侧反转时刻tre和倾侧角初始模值σ1的分段参数函数，所以将标准控制选为tre和σ1。
S2：标控飞行段：当飞行器距地面高度大于100km时，按标准控制飞行；当高度低于100km时，进入步骤S3；
具体地，该标准控制是指步骤S1中离线获得的倾侧反转时刻tre和倾侧角初始模值σ1。
将起控点高度定在100km的原因主要有两点：其一，距地面高度大于100km时，空气稀薄，大气对飞行器的减速效果不明显，这将导致飞行器速度随时间的变化曲线整体不具备较好的多项式拟合特性，将严重影响伪谱法的求解精度，而从100km以后可以获得较高的拟合精度；其二，空气稀薄的环境下，飞行器既不会受到气流扰动，也基本不具备气动修正能力。所以飞行器在100km以上时按标准控制飞行，不引入控制量的修正。
S3:跳跃段倾侧反转判断：按标准控制飞行，根据当前时间与倾侧角反转时刻的对应关系，判断是否进行倾侧反转，如果当前时刻未达到反转时刻时，不进行反转，进入步骤S4，如果当前时刻大于或等于反转时刻，且不是最后一个反转点时，重新建立非线性参数控制问题，进入步骤S4，当作为最后一个反转点时，进入步骤S8；
S4：跳跃段预测弹道积分：使用当前时刻的状态量作为积分初值，倾侧反转时刻tre和倾侧角初始模值σ1作为标准控制，通过弹道积分可以获得终端状态偏差δxf，以及全局多段弹道信息：状态量序列Xk、控制量序列Uk；
S5：精度判定：由终端状态偏差δxf求解终端横程偏差δχ，当δχ满足步骤S1设置的精度要求时，保持倾侧角反转时刻不变，进入步骤S3，如果不满足精度要求，进入步骤S6；
S6：跳跃段倾侧反转时刻修正：基于步骤S4获得的多段预测弹道积分，结合参数化控制和多段线性伪谱法，并利用等高度条件修正终端时间偏差，结合变分原理推导倾侧点时间改变量对终端状态偏差的影响。最终获得终端偏差与控制修正量之间的解析关系，并更新标准控制，即倾侧反转时刻tre和倾侧角初始模值σ1，进入步骤S7；
具体地，此步骤中首先将控制变量参数化处理，如此有利于降低控制实现的难度并提升控制规律的认知。此外，将运动方程线性化并结合线性伪谱法，可以获得倾侧反转时间的解析修正解。
1)多段非线性方程的线性化
考虑带有终端约束的再入动力学方程，能够表示为一般的非线性动力学方程的形式：
其中，状态量x＝[r θ φ V γ ψ]T。将方程(1.9)在预测弹道周围泰勒展开，忽略二阶以上高阶项，可以获得一组以状态偏差δx为自变量的误差传播动力学方程
其中，x＝xref-δx，u＝uref-δu。系数矩阵A是6×6的矩阵，系数矩阵B是6×1的向量，控制量偏差δu是标量。经过推导，矩阵A和矩阵B中各元素的表达式为
其中
A43＝-Ω2r(sinγsin2φ+cosγcosψcos2φ)
其中，考虑到全程飞行的马赫数在20以上，此时升阻力系数随马赫数变化不大，可以忽略CL、CD对地心距离r、飞行速度V的偏导。将涉及到的计算公式整理如下
考虑到控制量已参数化为时间t的连续函数，结合公式，控制量偏差δu可以表示为
则误差传播动力学方程可以表示为
其中，系数矩阵为6×1的向量，且
其中，表示向量的张量积。因为δσ1表示倾侧角初始模值的改变，所以受倾侧反转的影响，与的符号相反。
2)含有分段时间修正的多段线性伪谱修正算法
线性伪谱修正算法的核心思想是，通过将非线性问题在预测弹道周围拟线性化处理，得到以状态量偏差作为自变量的线性微分方程。得益于本发明的控制变量已经参数化处理，结合Gauss伪谱法，利用正交配点以及拉格朗日插值多项式，误差传播方程的终端偏差可以表示为初始偏差和控制量修正的解析函数。
本文所要处理的跳跃式再入飞行器横向制导问题是一个典型的多段非线性问题，前文已经论述，作为分段点的倾侧反转时刻，其值的选择将直接影响终端横向偏差。所以，本文将结合线性伪谱法，利用变分原理，推导终端偏差与分段时刻的解析关系，并时时修正倾侧反转点，实现在线制导。
因为Gauss伪谱中所使用的Legendre-Gauss正交节点分布在[-1,+1]，所以首先应把时间区间[t0,tf]映射到区间[-1,+1]，区间变换公式如下
则误差传播动力学方程可以表示为
选取N个Legendre-Gauss正交节点(τ1,τ2,...,τN)，并将端点值表示为τ0＝-1以及τf＝1，则N阶Lagrange插值多项式可以表示为
任意时刻的状态量可以通过插值公式拟合
其中，xi表示第i个插值点处的状态量xi＝x(τi)。
将上式两侧对τ求一阶导，可以得到状态量微分的插值拟合公式
构建N×N+1阶的微分逼近矩阵D，D矩阵各元素的具体表达形式为
将公式(1.18)带入方程(1.14)中，则线性动力学方程转化为一组基于LG配点的代数约束
其中，上标i(i＝1,2,3)表示第i段，假设状态量有s维，则有
[t1 t2 t3]T＝[tre tk tf]T
以式(1.19)为例，结合Gauss积分公式，第一段飞行过程的终端偏差可以由初始状态偏差配点处状态偏差δx1，以及初始倾侧角的模值偏差δσ1表示
其中，s×Ns阶系数矩阵W1为
将式(1.23)代入式(1.19)，消去配点处的状态偏差δx1
类似地，可以将第二段的终端偏差第三段的终端偏差表示为
定义误差传递矩阵
其中，上标i表示第i段。
考虑到弹道积分的初值即为当前实际运动状态，所以初始状态偏差始终为零。结合公式(1.24)和公式(1.25)终端状态偏差最终可以表示为
需要特别指出的是，式(1.27)中所求解的终端状态偏差是指，由当前初始状态在标准控制的作用下，经过同等时间后，所获得的终端状态值x(tf)与终端状态约束xf之间的偏差，即
δxf＝x(tf)-xf (1.28)
一般性的，所求解的δxf可直接用于以时间为分段标志的多段问题。但是，对于以其他状态量为分段标志的问题，比如本发明中所阐述的跳跃式再入弹道，如图4所示，其二次再入点是以距地面高度等于81.5km确定的，此时的终端偏差应该表示为
δxhf＝x(hf)-xhf (1.29)
δxhf可由δxf修正后得到，为阐述二者间的关系，先从一个简单算例开始。如图4所示，黑色实线表示标准轨迹，距地面高度81.5km处的二次再入点所对应的剩余射程为2000km，严格满足终端射程约束。假设飞行器再次跃出大气层，即在Kepler起点处，有模值为20m/s的速度扰动，所对应的偏差轨迹如虚线所示。飞行器经历与标准轨迹相同的飞行时间后，对应的等时积分终端值表示为红色加号。同时，结合公式
可以通过初始速度扰动解析地求出终端偏差，将矩阵计算的结果表示为黑色圈状。可以看出二者近似重合，这表明公式(1.30)所描述的矩阵计算方法是行之有效的。
但是，二次再入点是以高度重新达到81.5km确定的，不由飞行时间决定。沿着红色虚线所代表的偏差轨迹，找到高度等于81.5km的点并以红色方块状表示。该点与基准位置之间的偏差才是方程(1.30)中所描述的终端偏差考虑到二次再入点的高度约束可以表示为
结合标准轨迹上基准位置处的导数信息，对公式(1.30)作如下修正
引入时间修正矩阵Y5用以描述终端时刻的修正关系
方程(1.30)可以进一步表示为
图4中，方程(1.33)修正矩阵的计算结果被表示为黑色菱形，结果表明其与红色实心方块所表示的等高积分结果基本一致。所以方程(1.33)中的修正思想完全适用于求解以高度为分界点的终端偏差。
事实上，不仅仅第三段的终端状态是以高度81.5km作为分界点，如图2所示，跳跃式再入轨迹的第二段与第三段的分段条件也是以高度81.5km确定的，所以，类似地，对公式(1.27)做出修正
其中，时间修正矩阵Y4用以描述二三段分界点处的修正关系
标准控制由倾侧角初始模值σ1和倾侧角反转时刻tre两个参数组成
其中，上述公式仅仅描述了反转时刻不改变的情况。进一步地，为了考虑δtre对终端偏差的影响，还需要考虑状态量对时间的变分。如图5所示，δx(tre)是δtre所引起的状态量的改变，可以由倾侧反转前后的微分方程表示
δx(tre)＝[f1(x(tre),u,tre)-f2(x(tre),u,tre)]δtre (1.35)
式中f1(x(tre),u,tre),f2(x(tre),u,tre)分别表示方程组(1.1)在与的取值。
在本发明中，方程(1.35)的关系可以具体表示为
其中，的第7个分量的符号，与第一段倾侧角符号相同。
结合方程(1.34)，并考虑对高度的修正，定义Sx2维的系数矩阵M
则全程终端偏差可以表示为
选取两组终端约束函数，用于求解控制量的修正值
对上式两侧求微分
将公式(1.40)整理为
其中
定义系数矩阵Z1
则
所以，标准控制的修正量可以解析地表示为
综上，通过步骤S6，标准控制可以迭代更新
其中，下标k表示第k次迭代。
S7：重新建立参数控制问题，仅将倾侧角模值参数作为标准控制，用以修正纵向射程偏差。当飞行器高度再次超过81.5km，此时认为飞行器脱离大气层，不受气动力的影响，进入步骤S11，否则，进入步骤S8；
S8：跳跃段预测弹道积分：使用当前时刻的状态量作为积分初值，在标准控制的作用下，通过弹道积分可以获得终端状态偏差δxf，以及全局多段弹道信息：状态量序列Xk、控制量序列Uk；
S9：精度判定：根据终端状态偏差δxf求解终端射程偏差δS，当δS满足步骤S1设置的精度要求时，保持标准控制不变，进入步骤S7，否则，终端偏差过大，进入步骤S10；
S10：跳跃段控制参数修正：基于步骤S8获得的多段预测弹道积分，结合参数化控制、多段伪谱法，获得终端射程偏差与倾侧角模值参数修正值之间的解析关系，并更新标准控制，进入步骤S7；
步骤S10具体如下：
区别于步骤S6的是，飞行器此时处于第二段与第三段的飞行过程中，不再需要考虑倾侧反转，所以有别于公式(1.27)，此时的终端偏差可以表示为
其中，Sx1维的系数矩阵N是
为求解δσ1，仅将射程作为终端约束函数
y1＝cos(S)＝sinφsinΦ+cosΦcosφcos(Θ-θ) (1.46)
将式(1.46)两侧求微分，并进一步表示为
定义系数矩阵Z2
则
Z2δσ1＝δy1
所以，标准控制的修正量可以解析地表示为
综上，通过步骤S10标准控制可以迭代更新
σ1|k+1＝σ1|k-δσ1 (1.50)
其中，下标k表示第k次迭代。
S11：Kepler轨道飞行：飞行器沿标准Kepler轨道飞行，直至二次再入时距地面高度小于等于81.5km，Kepler段飞行结束。此时再入飞行器完成跳跃段制导飞行，仿真结束。
下面举例说明本发明执导方法终端横向偏差的控制精度及迭代计算速度。
算例1
此处将提供一个通过迭代更新倾侧角模值参数进而修正终端射程偏差的算例，用以阐明当分段条件不由时间决定时，引入控制补偿的必要性，并进一步阐述算法的实施过程。
如图6中虚线所示，在距地面80km的起始点，如果飞行器按照初始倾侧角模值参数σ0＝83.7°继续飞行，在经历以81.5km作为高度分界线的Kepler轨道飞行后，与2000km剩余射程之间存在724.87km的终端纵向射程偏差。在不考虑倾侧角反转时刻以及横程偏差的前提下，结合前文的推导，终端偏差δxf可以由倾侧角初始模值的修正值δσ1解析表示
上式中，有必要再次强调的是，由于跳跃式再入弹道Kepler段的起点与终点的划分完全由高度确定，必须引入时间修正矩阵Y4与Y5，才能确保积分终端满足81.5km的高度分界条件。
为求解δσ1，仅将射程作为终端约束函数
y1＝cos(S)＝sinφsinΦ+cosΦcosφcos(Θ-θ) (2.1)
将上式两侧求微分，并进一步表示为
定义系数矩阵Z2
则标准控制的修正量可以解析地表示为
通过上式求得δσ1＝0.178°，更新控制后的修正轨迹如图6中实线所示，仅经过一次修正，终端射程偏差就缩减为58.2km。如果不引入时间修正矩阵Y4与Y5，经过迭代修正后的积分弹道如图6中虚点线所示，这也从反面再次验证了本发明方法的准确性。
算例2
此处将提供一个低升阻比跳跃式再入飞行器的全程制导算例，旨在通过分析飞行器二次再入点的终端射程偏差和横程偏差，检验本发明所提出方法的准确性和实用性。算例的仿真条件按照步骤1中的描述设置，其中气动数据按照美国猎户座载人飞船选取，飞行器的初始状态如下表所示
初始的倾侧角模值参数和倾侧角反转时间分别设为83.7°和175s，如图7和图8所示，如果不采取控制修正，飞行器在第二次再入时的射程偏差和横程偏差分别为770.3km和161.7km，严重影响后续的末段制导。采用本发明提出的低升阻比跳跃式再入飞行器制导方法后，终端射程偏差和横程偏差分别减小为2.2km和3.2km，为后续的末制导段创造了良好的初始条件。
将制导算法在使用i7-7700k处理器的个人电脑中安装的MATLAB2016a中运行，全程共执行四次修正计算，其中倾侧翻转后修正一次，每次的计算时间以及修正指令如下表所示
可以看出，由于本发明提出的算法使用了梯度信息，所以只需要四次修正，就能使二次再入点的终端偏差小于5km，为末段制导提供了良好的初始条件。而且，得益于控制量修正值与终端偏差之间的解析关系，每一次控制量的更新都不需要迭代积分，单次计算耗时仅为2ms，非常适合在线应用。
本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的装置而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。
对所公开的实施例的上述说明，使本领域专业技术人员能够实现或使用本发明。对这些实施例的多种修改对本领域的专业技术人员来说将是显而易见的，本文中所定义的一般原理可以在不脱离本发明的精神或范围的情况下，在其它实施例中实现。因此，本发明将不会被限制于本文所示的这些实施例，而是要符合与本文所公开的原理和新颖特点相一致的最宽的范围。
Claims (3)
1.一种基于线性伪谱的跳跃式再入飞行器在线弹道规划制导方法，其特征在于，包括如下步骤：
S1：初始化：设置初始、终端的计算仿真参数，包括终端纵向射程偏差δS和横程落点偏差δχ的精度要求，通过离线弹道优化与参数化处理，将倾侧角变化规律表示为倾侧反转时刻tre和倾侧角初始模值σ1的多段参数函数，并将tre和σ1作为标准控制；
S2：标控飞行段：当飞行器距地面高度大于100km时，按标准控制飞行；当高度低于100km时，进入步骤(3)；
S3：跳跃段倾侧反转判断：按标准控制飞行，根据当前时间与倾侧角反转时刻的对应关系，判断是否进行倾侧反转，如果当前时刻未达到反转时刻时，不进行反转，进入步骤S4，如果当前时刻大于或等于反转时刻，且不是最后一个反转点时，重新建立非线性参数控制问题，进入步骤S4，当作为最后一个反转点时，进入步骤S7；
S4：跳跃段预测弹道积分：使用当前时刻的状态量作为积分初值，所述的标准控制作为控制输入，通过弹道积分获得终端状态偏差δxf，以及全局多段弹道信息：状态量序列Xk、控制量序列Uk；
S5：精度判定：由终端状态偏差δxf求解终端纵向射程偏差δS和横程落点偏差δχ，当二者都满足步骤S1设置的精度要求时，保持标准控制不变，进入步骤S3，否则，进入步骤S6；
S6：跳跃段倾侧反转时刻修正：基于步骤S4获得的多段预测弹道积分，获得终端偏差与控制修正量之间的解析关系，并更新标准控制，即倾侧反转时刻tre和倾侧角初始模值σ1，
考虑到弹道积分的初值即为当前实际运动状态，所以初始状态偏差始终为零，终端状态偏差最终表示为
<mrow>
<msub>
<mi>&delta;x</mi>
<mi>f</mi>
</msub>
<mo>=</mo>
<msubsup>
<mi>&delta;x</mi>
<mi>f</mi>
<mn>3</mn>
</msubsup>
<mo>=</mo>
<msubsup>
<mi>K</mi>
<mi>x</mi>
<mn>3</mn>
</msubsup>
<mrow>
<mo>(</mo>
<msubsup>
<mi>K</mi>
<mi>x</mi>
<mn>2</mn>
</msubsup>
<msubsup>
<mi>K</mi>
<mi>u</mi>
<mn>1</mn>
</msubsup>
<mo>+</mo>
<msubsup>
<mi>K</mi>
<mi>u</mi>
<mn>2</mn>
</msubsup>
<mo>)</mo>
</mrow>
<msub>
<mi>&delta;&sigma;</mi>
<mn>1</mn>
</msub>
</mrow>
其中，K矩阵表示误差传递矩阵，上标数字分别表示第1段、第2段和第3段，下标x和u分别表示相应的终端状态偏差是由状态量偏差或控制量偏差导致；
二次再入点是以高度重新达到81.5km确定的，不由飞行时间决定；考虑到二次再入点的高度约束表示为
<mrow>
<msub>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<msub>
<mi>&sigma;</mi>
<mn>1</mn>
</msub>
</mtd>
</mtr>
<mtr>
<mtd>
<msub>
<mi>t</mi>
<mrow>
<mi>r</mi>
<mi>e</mi>
</mrow>
</msub>
</mtd>
</mtr>
</mtable>
</mfenced>
<mrow>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
</msub>
<mo>=</mo>
<msub>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<msub>
<mi>&sigma;</mi>
<mn>1</mn>
</msub>
</mtd>
</mtr>
<mtr>
<mtd>
<msub>
<mi>t</mi>
<mrow>
<mi>r</mi>
<mi>e</mi>
</mrow>
</msub>
</mtd>
</mtr>
</mtable>
</mfenced>
<mi>k</mi>
</msub>
<mo>-</mo>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mi>&delta;</mi>
<msub>
<mi>&sigma;</mi>
<mn>1</mn>
</msub>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<msub>
<mi>&delta;t</mi>
<mrow>
<mi>r</mi>
<mi>e</mi>
</mrow>
</msub>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
</mrow>
结合标准轨迹上基准位置处的导数信息，终端状态偏差修正为
<mrow>
<msubsup>
<mi>&delta;x</mi>
<mi>f</mi>
<mn>3</mn>
</msubsup>
<mo>=</mo>
<msubsup>
<mi>K</mi>
<mi>x</mi>
<mn>3</mn>
</msubsup>
<msubsup>
<mi>&delta;x</mi>
<mn>0</mn>
<mn>3</mn>
</msubsup>
<mo>-</mo>
<msubsup>
<mover>
<mi>x</mi>
<mo>&CenterDot;</mo>
</mover>
<mi>f</mi>
<mn>3</mn>
</msubsup>
<mfrac>
<mrow>
<msubsup>
<mi>&delta;h</mi>
<mi>f</mi>
<mn>3</mn>
</msubsup>
</mrow>
<msubsup>
<mover>
<mi>h</mi>
<mo>&CenterDot;</mo>
</mover>
<mi>f</mi>
<mn>3</mn>
</msubsup>
</mfrac>
</mrow>
δx0为初始状态偏差，hf表示由等时弹道积分获得的终端高度，即二次再入点的高度，终端状态偏差进一步表示为
<mrow>
<msubsup>
<mi>&delta;x</mi>
<mi>f</mi>
<mn>3</mn>
</msubsup>
<mo>=</mo>
<msub>
<mi>Y</mi>
<mn>5</mn>
</msub>
<msubsup>
<mi>K</mi>
<mi>x</mi>
<mn>3</mn>
</msubsup>
<msubsup>
<mi>&delta;x</mi>
<mn>0</mn>
<mn>3</mn>
</msubsup>
</mrow>
其中，Y5为时间修正矩阵，用以描述终端时刻的修正关系；
同理，不仅仅第三段的终端状态是以高度81.5km作为分界点，跳跃式再入轨迹的第二段与第三段的分段条件也是以高度81.5km确定的，对终端状态偏差做出修正
<mrow>
<msub>
<mi>&delta;x</mi>
<mi>f</mi>
</msub>
<mo>=</mo>
<msubsup>
<mi>&delta;x</mi>
<mi>f</mi>
<mn>3</mn>
</msubsup>
<mo>=</mo>
<msub>
<mi>Y</mi>
<mn>5</mn>
</msub>
<msubsup>
<mi>K</mi>
<mi>x</mi>
<mn>3</mn>
</msubsup>
<msub>
<mi>Y</mi>
<mn>4</mn>
</msub>
<mrow>
<mo>(</mo>
<msubsup>
<mi>K</mi>
<mi>x</mi>
<mn>2</mn>
</msubsup>
<msubsup>
<mi>K</mi>
<mi>u</mi>
<mn>1</mn>
</msubsup>
<mo>+</mo>
<msubsup>
<mi>K</mi>
<mi>u</mi>
<mn>2</mn>
</msubsup>
<mo>)</mo>
</mrow>
<msub>
<mi>&delta;&sigma;</mi>
<mn>1</mn>
</msub>
</mrow>
其中，时间修正矩阵Y4用以描述二三段分界点处的修正关系；
进入步骤S7；
S7：重新建立参数控制：仅将倾侧角模值参数作为标准控制，用以修正纵向射程偏差，当飞行器高度再次超过81.5km，此时认为飞行器脱离大气层，不受气动力的影响，进入步骤S11，否则，进入步骤S8；
S8：跳跃段预测弹道积分：在当前标准控制的作用下，通过弹道积分获得终端状态偏差δxf，以及全局多段弹道信息：状态量序列Xk、控制量序列Uk；
S9：精度判定：根据终端状态偏差δxf求解终端射程偏差δS，当δS满足步骤S1设置的精度要求时，保持标准控制不变，进入步骤S7，否则，终端偏差过大，进入步骤S10；
S10：跳跃段控制参数修正：基于步骤S8获得的多段预测弹道积分，计算得到终端射程偏差与倾侧角模值参数修正值之间的解析关系，并进一步更新标准控制，进入步骤S7；
S11：Kepler轨道飞行：飞行器沿标准Kepler轨道飞行，直至二次再入时距地面高度小于等于81.5km，Kepler段飞行结束，此时再入飞行器完成跳跃段制导飞行，仿真结束。
2.根据权利要求1所述的基于线性伪谱的跳跃式再入飞行器在线弹道规划制导方法，其特征在于，所述步骤S6中：
综合考虑倾侧角初始模值σ1和倾侧角反转时刻tre的影响，终端偏差可以进一步整理为
<mrow>
<msub>
<mi>&delta;x</mi>
<mi>f</mi>
</msub>
<mo>=</mo>
<mi>M</mi>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mi>&delta;</mi>
<msub>
<mi>&sigma;</mi>
<mn>1</mn>
</msub>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<msub>
<mi>&delta;t</mi>
<mrow>
<mi>r</mi>
<mi>e</mi>
</mrow>
</msub>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
</mrow>
其中，Sx2维的系数矩阵M表征终端偏差与控制量修正值之间的解析关系。
选取两组终端约束函数，用于求解控制量的修正值
y1＝cos(S)＝sinφsinΦ+cosΦcosφcos(Θ-θ)
<mrow>
<msub>
<mi>y</mi>
<mn>2</mn>
</msub>
<mo>=</mo>
<mi>t</mi>
<mi>a</mi>
<mi>n</mi>
<mi>&Psi;</mi>
<mo>-</mo>
<mi>t</mi>
<mi>a</mi>
<mi>n</mi>
<mi>&psi;</mi>
<mo>=</mo>
<mfrac>
<mrow>
<mi>s</mi>
<mi>i</mi>
<mi>n</mi>
<mrow>
<mo>(</mo>
<mi>&Theta;</mi>
<mo>-</mo>
<mi>&theta;</mi>
<mo>)</mo>
</mrow>
</mrow>
<mrow>
<mi>c</mi>
<mi>o</mi>
<mi>s</mi>
<mi>&phi;</mi>
<mi>t</mi>
<mi>a</mi>
<mi>n</mi>
<mi>&Phi;</mi>
<mo>-</mo>
<mi>s</mi>
<mi>i</mi>
<mi>n</mi>
<mi>&phi;</mi>
<mi>c</mi>
<mi>o</mi>
<mi>s</mi>
<mrow>
<mo>(</mo>
<mi>&Theta;</mi>
<mo>-</mo>
<mi>&theta;</mi>
<mo>)</mo>
</mrow>
</mrow>
</mfrac>
<mo>-</mo>
<mi>t</mi>
<mi>a</mi>
<mi>n</mi>
<mi>&psi;</mi>
</mrow>
其中，S为飞行器纵向射程，θ和φ分别表示飞行器所在位置的经度和纬度，其单位是弧度rad；Θ和Φ分别表示飞行器着陆点的经度和纬度；将飞行器质心、着陆点、地球圆心三点围成的射面与飞行器当地北向的夹角定义为视线角Ψ，ψ表示速度矢量在当地水平面的投影与正北方向的夹角，顺时针方向为正，二者的单位是弧度；对上式两侧求微分
<mrow>
<msub>
<mi>&delta;y</mi>
<mn>1</mn>
</msub>
<mo>=</mo>
<mfrac>
<mrow>
<mo>&part;</mo>
<msub>
<mi>y</mi>
<mn>1</mn>
</msub>
</mrow>
<mrow>
<mo>&part;</mo>
<mi>&theta;</mi>
</mrow>
</mfrac>
<mi>&delta;</mi>
<mi>&theta;</mi>
<mo>+</mo>
<mfrac>
<mrow>
<mo>&part;</mo>
<msub>
<mi>y</mi>
<mn>1</mn>
</msub>
</mrow>
<mrow>
<mo>&part;</mo>
<mi>&phi;</mi>
</mrow>
</mfrac>
<mi>&delta;</mi>
<mi>&phi;</mi>
</mrow>
<mrow>
<msub>
<mi>&delta;y</mi>
<mn>2</mn>
</msub>
<mo>=</mo>
<mfrac>
<mrow>
<mo>&part;</mo>
<msub>
<mi>y</mi>
<mn>2</mn>
</msub>
</mrow>
<mrow>
<mo>&part;</mo>
<mi>&theta;</mi>
</mrow>
</mfrac>
<mi>&delta;</mi>
<mi>&theta;</mi>
<mo>+</mo>
<mfrac>
<mrow>
<mo>&part;</mo>
<msub>
<mi>y</mi>
<mn>2</mn>
</msub>
</mrow>
<mrow>
<mo>&part;</mo>
<mi>&phi;</mi>
</mrow>
</mfrac>
<mi>&delta;</mi>
<mi>&phi;</mi>
<mo>+</mo>
<mfrac>
<mrow>
<mo>&part;</mo>
<msub>
<mi>y</mi>
<mn>2</mn>
</msub>
</mrow>
<mrow>
<mo>&part;</mo>
<mi>&psi;</mi>
</mrow>
</mfrac>
<mi>&delta;</mi>
<mi>&psi;</mi>
</mrow>
并定义系数矩阵Z1
<mrow>
<msub>
<mi>Z</mi>
<mn>1</mn>
</msub>
<mo>=</mo>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mfrac>
<mrow>
<mo>&part;</mo>
<msub>
<mi>y</mi>
<mn>1</mn>
</msub>
</mrow>
<mrow>
<mo>&part;</mo>
<mi>&theta;</mi>
</mrow>
</mfrac>
</mtd>
<mtd>
<mfrac>
<mrow>
<mo>&part;</mo>
<msub>
<mi>y</mi>
<mn>1</mn>
</msub>
</mrow>
<mrow>
<mo>&part;</mo>
<mi>&phi;</mi>
</mrow>
</mfrac>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
</mtr>
<mtr>
<mtd>
<mfrac>
<mrow>
<mo>&part;</mo>
<msub>
<mi>y</mi>
<mn>2</mn>
</msub>
</mrow>
<mrow>
<mo>&part;</mo>
<mi>&theta;</mi>
</mrow>
</mfrac>
</mtd>
<mtd>
<mfrac>
<mrow>
<mo>&part;</mo>
<msub>
<mi>y</mi>
<mn>2</mn>
</msub>
</mrow>
<mrow>
<mo>&part;</mo>
<mi>&phi;</mi>
</mrow>
</mfrac>
</mtd>
<mtd>
<mfrac>
<mrow>
<mo>&part;</mo>
<msub>
<mi>y</mi>
<mn>2</mn>
</msub>
</mrow>
<mrow>
<mo>&part;</mo>
<mi>&psi;</mi>
</mrow>
</mfrac>
</mtd>
</mtr>
</mtable>
</mfenced>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<mn>1</mn>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
</mtr>
<mtr>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<mn>1</mn>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
</mtr>
<mtr>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<mn>1</mn>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
</mtr>
</mtable>
</mfenced>
<mi>M</mi>
</mrow>
所以，标准控制的修正量可以解析地表示为
<mrow>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mi>&delta;</mi>
<msub>
<mi>&sigma;</mi>
<mn>1</mn>
</msub>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<msub>
<mi>&delta;t</mi>
<mrow>
<mi>r</mi>
<mi>e</mi>
</mrow>
</msub>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>=</mo>
<msubsup>
<mi>Z</mi>
<mn>1</mn>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msubsup>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mi>&delta;</mi>
<msub>
<mi>y</mi>
<mn>1</mn>
</msub>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<msub>
<mi>&delta;y</mi>
<mn>2</mn>
</msub>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
</mrow>
标准控制迭代更新为
<mrow>
<msub>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<msub>
<mi>&sigma;</mi>
<mn>1</mn>
</msub>
</mtd>
</mtr>
<mtr>
<mtd>
<msub>
<mi>t</mi>
<mrow>
<mi>r</mi>
<mi>e</mi>
</mrow>
</msub>
</mtd>
</mtr>
</mtable>
</mfenced>
<mrow>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
</msub>
<mo>=</mo>
<msub>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<msub>
<mi>&sigma;</mi>
<mn>1</mn>
</msub>
</mtd>
</mtr>
<mtr>
<mtd>
<msub>
<mi>t</mi>
<mrow>
<mi>r</mi>
<mi>e</mi>
</mrow>
</msub>
</mtd>
</mtr>
</mtable>
</mfenced>
<mi>k</mi>
</msub>
<mo>-</mo>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mi>&delta;</mi>
<msub>
<mi>&sigma;</mi>
<mn>1</mn>
</msub>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<msub>
<mi>&delta;t</mi>
<mrow>
<mi>r</mi>
<mi>e</mi>
</mrow>
</msub>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
</mrow>
其中，下标k表示第k次迭代。
3.根据权利要求1所述的基于线性伪谱的跳跃式再入飞行器在线弹道规划制导方法，其特征在于，所述步骤S10中：
飞行器处于第二段与第三段的飞行过程中，不考虑倾侧反转时刻，终端偏差表示为
<mrow>
<msub>
<mi>&delta;x</mi>
<mi>f</mi>
</msub>
<mo>=</mo>
<msubsup>
<mi>&delta;x</mi>
<mi>f</mi>
<mn>3</mn>
</msubsup>
<mo>=</mo>
<msub>
<mi>N</mi>
<mrow>
<mi>s</mi>
<mo>&times;</mo>
<mn>1</mn>
</mrow>
</msub>
<msub>
<mi>&delta;&sigma;</mi>
<mn>1</mn>
</msub>
</mrow>
其中，Sx1维的系数矩阵N为是
<mrow>
<msub>
<mi>N</mi>
<mrow>
<mi>s</mi>
<mo>&times;</mo>
<mn>1</mn>
</mrow>
</msub>
<mo>=</mo>
<msub>
<mi>Y</mi>
<mn>5</mn>
</msub>
<msubsup>
<mi>K</mi>
<mi>x</mi>
<mn>3</mn>
</msubsup>
<msub>
<mi>Y</mi>
<mn>4</mn>
</msub>
<msubsup>
<mi>K</mi>
<mi>u</mi>
<mn>2</mn>
</msubsup>
<msub>
<mi>&delta;&sigma;</mi>
<mn>1</mn>
</msub>
</mrow>
仅将纵向射程作为终端约束函数
y1＝cos(S)＝sinφsinΦ+cosΦcosφcos(Θ-θ)
标准控制的修正量解析地表示为
<mrow>
<msub>
<mi>&delta;&sigma;</mi>
<mn>1</mn>
</msub>
<mo>=</mo>
<msubsup>
<mi>Z</mi>
<mn>2</mn>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msubsup>
<msub>
<mi>&delta;y</mi>
<mn>1</mn>
</msub>
</mrow>
其中，定义系数矩阵Z2
<mrow>
<msub>
<mi>Z</mi>
<mn>2</mn>
</msub>
<mo>=</mo>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mfrac>
<mrow>
<mo>&part;</mo>
<msub>
<mi>y</mi>
<mn>1</mn>
</msub>
</mrow>
<mrow>
<mo>&part;</mo>
<mi>&theta;</mi>
</mrow>
</mfrac>
</mtd>
<mtd>
<mfrac>
<mrow>
<mo>&part;</mo>
<msub>
<mi>y</mi>
<mn>1</mn>
</msub>
</mrow>
<mrow>
<mo>&part;</mo>
<mi>&phi;</mi>
</mrow>
</mfrac>
</mtd>
</mtr>
</mtable>
</mfenced>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<mn>1</mn>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
</mtr>
<mtr>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<mn>1</mn>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
</mtr>
</mtable>
</mfenced>
<mi>N</mi>
</mrow>
迭代更新标准控制为
σ1|k+1＝σ1|k-δσ1
其中，下标k表示第k次迭代。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201711058719.0A CN107861517B (zh) | 2017-11-01 | 2017-11-01 | 基于线性伪谱的跳跃式再入飞行器在线弹道规划制导方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201711058719.0A CN107861517B (zh) | 2017-11-01 | 2017-11-01 | 基于线性伪谱的跳跃式再入飞行器在线弹道规划制导方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN107861517A true CN107861517A (zh) | 2018-03-30 |
| CN107861517B CN107861517B (zh) | 2020-06-16 |
Family
ID=61696428
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201711058719.0A Active CN107861517B (zh) | 2017-11-01 | 2017-11-01 | 基于线性伪谱的跳跃式再入飞行器在线弹道规划制导方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107861517B (zh) |
Cited By (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109062241A (zh) * | 2018-09-29 | 2018-12-21 | 北京航空航天大学 | 基于线性伪谱模型预测控制的自主全射向再入制导方法 |
| CN109091109A (zh) * | 2018-07-02 | 2018-12-28 | 南京大学 | 基于全矩阵滤波和时间反转算子的优化型光声断层成像的图像重构方法 |
| CN109270960A (zh) * | 2018-12-05 | 2019-01-25 | 中南大学 | 基于Radau伪谱法的在线最优反馈再入制导方法 |
| CN110425943A (zh) * | 2019-08-06 | 2019-11-08 | 西北工业大学 | 面向变质心飞行器的工程化再入制导方法 |
| CN110955859A (zh) * | 2019-11-15 | 2020-04-03 | 北京航空航天大学 | 一种基于预定飞行轨的实际飞行轨道气动参数的计算方法 |
| CN111351401A (zh) * | 2018-12-21 | 2020-06-30 | 北京理工大学 | 应用于捷联导引头制导飞行器的防侧偏制导方法 |
| CN112034879A (zh) * | 2019-12-16 | 2020-12-04 | 哈尔滨工业大学 | 一种基于高度-射程比的标准轨迹跟踪制导方法 |
| CN112124631A (zh) * | 2020-09-11 | 2020-12-25 | 中国运载火箭技术研究院 | 自适应分档控制方法、控制装置及存储介质 |
| CN112861253A (zh) * | 2020-12-25 | 2021-05-28 | 航天科工微电子系统研究院有限公司 | 一种异质多飞行器复杂构型协同弹道规划方法 |
| CN113311871A (zh) * | 2021-06-02 | 2021-08-27 | 上海机电工程研究所 | 跳跃-滑翔导弹的制导律优化方法及系统 |
| CN115407664A (zh) * | 2022-11-01 | 2022-11-29 | 北京航天自动控制研究所 | 一种基于神经网络训练的非程序制导方法 |
Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6629668B1 (en) * | 2002-02-04 | 2003-10-07 | The United States Of America As Represented By The Secretary Of The Army | Jump correcting projectile system |
| CN103863579A (zh) * | 2014-03-31 | 2014-06-18 | 北京控制工程研究所 | 一种深空探测返回过程的预测校正制导方法 |
| CN103921957A (zh) * | 2014-01-16 | 2014-07-16 | 中国人民解放军国防科学技术大学 | 一种探月飞船跳跃式再入的跃起点能量管理方法 |
| CN104850129A (zh) * | 2014-12-19 | 2015-08-19 | 北京控制工程研究所 | 一种跳跃式再入的射向预偏置横向制导方法 |
| CN104843197A (zh) * | 2014-12-19 | 2015-08-19 | 北京控制工程研究所 | 一种跳跃式再入的双环制导方法 |
| CN106586033A (zh) * | 2017-01-13 | 2017-04-26 | 北京航空航天大学 | 自适应分段的多段线性伪谱广义标控脱靶量再入制导方法 |
| CN105222648B (zh) * | 2015-09-11 | 2017-04-26 | 北京航空航天大学 | 一种线性伪谱广义标控脱靶量制导控制方法 |
-
2017
- 2017-11-01 CN CN201711058719.0A patent/CN107861517B/zh active Active
Patent Citations (7)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US6629668B1 (en) * | 2002-02-04 | 2003-10-07 | The United States Of America As Represented By The Secretary Of The Army | Jump correcting projectile system |
| CN103921957A (zh) * | 2014-01-16 | 2014-07-16 | 中国人民解放军国防科学技术大学 | 一种探月飞船跳跃式再入的跃起点能量管理方法 |
| CN103863579A (zh) * | 2014-03-31 | 2014-06-18 | 北京控制工程研究所 | 一种深空探测返回过程的预测校正制导方法 |
| CN104850129A (zh) * | 2014-12-19 | 2015-08-19 | 北京控制工程研究所 | 一种跳跃式再入的射向预偏置横向制导方法 |
| CN104843197A (zh) * | 2014-12-19 | 2015-08-19 | 北京控制工程研究所 | 一种跳跃式再入的双环制导方法 |
| CN105222648B (zh) * | 2015-09-11 | 2017-04-26 | 北京航空航天大学 | 一种线性伪谱广义标控脱靶量制导控制方法 |
| CN106586033A (zh) * | 2017-01-13 | 2017-04-26 | 北京航空航天大学 | 自适应分段的多段线性伪谱广义标控脱靶量再入制导方法 |
Non-Patent Citations (4)
| Title |
|---|
| 曾亮: "探月飞船跳跃式返回再入制导方法研究", 《中国优秀硕士学位论文全文数据库》 * |
| 杨良等: "多约束在线高斯伪谱末制导方法", 《弹道学报》 * |
| 郭敏文 等: "基于改进Gauss伪谱法的探月返回器跳跃式再入轨迹优化设计", 《航天控制》 * |
| 陈万春 等: "Longitudinal characteristics of steady glide trajectory for hypersonic vehicle", 《国际控制自动信息科学会议(ICCAIS)》 * |
Cited By (18)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109091109B (zh) * | 2018-07-02 | 2021-04-20 | 南京大学 | 基于全矩阵滤波和时间反转算子的优化型光声断层成像的图像重构方法 |
| CN109091109A (zh) * | 2018-07-02 | 2018-12-28 | 南京大学 | 基于全矩阵滤波和时间反转算子的优化型光声断层成像的图像重构方法 |
| CN109062241A (zh) * | 2018-09-29 | 2018-12-21 | 北京航空航天大学 | 基于线性伪谱模型预测控制的自主全射向再入制导方法 |
| CN109270960A (zh) * | 2018-12-05 | 2019-01-25 | 中南大学 | 基于Radau伪谱法的在线最优反馈再入制导方法 |
| CN111351401A (zh) * | 2018-12-21 | 2020-06-30 | 北京理工大学 | 应用于捷联导引头制导飞行器的防侧偏制导方法 |
| CN111351401B (zh) * | 2018-12-21 | 2022-12-23 | 北京理工大学 | 应用于捷联导引头制导飞行器的防侧偏制导方法 |
| CN110425943A (zh) * | 2019-08-06 | 2019-11-08 | 西北工业大学 | 面向变质心飞行器的工程化再入制导方法 |
| CN110425943B (zh) * | 2019-08-06 | 2021-05-07 | 西北工业大学 | 面向变质心飞行器的工程化再入制导方法 |
| CN110955859A (zh) * | 2019-11-15 | 2020-04-03 | 北京航空航天大学 | 一种基于预定飞行轨的实际飞行轨道气动参数的计算方法 |
| CN110955859B (zh) * | 2019-11-15 | 2021-05-14 | 北京航空航天大学 | 一种基于预定飞行轨的实际飞行轨道气动参数的计算方法 |
| CN112034879A (zh) * | 2019-12-16 | 2020-12-04 | 哈尔滨工业大学 | 一种基于高度-射程比的标准轨迹跟踪制导方法 |
| CN112124631A (zh) * | 2020-09-11 | 2020-12-25 | 中国运载火箭技术研究院 | 自适应分档控制方法、控制装置及存储介质 |
| CN112861253A (zh) * | 2020-12-25 | 2021-05-28 | 航天科工微电子系统研究院有限公司 | 一种异质多飞行器复杂构型协同弹道规划方法 |
| CN112861253B (zh) * | 2020-12-25 | 2023-08-22 | 航天科工微电子系统研究院有限公司 | 一种异质多飞行器复杂构型协同弹道规划方法 |
| CN113311871A (zh) * | 2021-06-02 | 2021-08-27 | 上海机电工程研究所 | 跳跃-滑翔导弹的制导律优化方法及系统 |
| CN113311871B (zh) * | 2021-06-02 | 2022-12-13 | 上海机电工程研究所 | 跳跃-滑翔导弹的制导律优化方法及系统 |
| CN115407664A (zh) * | 2022-11-01 | 2022-11-29 | 北京航天自动控制研究所 | 一种基于神经网络训练的非程序制导方法 |
| CN115407664B (zh) * | 2022-11-01 | 2023-02-03 | 北京航天自动控制研究所 | 一种基于神经网络训练的非程序制导方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN107861517B (zh) | 2020-06-16 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107861517A (zh) | 基于线性伪谱的跳跃式再入飞行器在线弹道规划制导方法 | |
| Shim | Hierarchical flight control system synthesis for rotorcraft-based unmanned aerial vehicles | |
| CN111306989B (zh) | 一种基于平稳滑翔弹道解析解的高超声速再入制导方法 | |
| CN106371312B (zh) | 基于模糊控制器的升力式再入预测-校正制导方法 | |
| CN103488814B (zh) | 一种适用于再入飞行器姿态控制的闭环仿真系统 | |
| CN103838914B (zh) | 一种高超声速飞行器滑翔段弹道解析求解方法 | |
| CN110015446B (zh) | 一种半解析的火星进入制导方法 | |
| CN109062241B (zh) | 基于线性伪谱模型预测控制的自主全射向再入制导方法 | |
| Xie et al. | Highly constrained entry trajectory generation | |
| CN107992074A (zh) | 一种基于飞行路径角规划的再入轨迹设计方法 | |
| De Marco et al. | A general solution to the aircraft trim problem | |
| CN107121929B (zh) | 基于线性协方差模型预测控制的鲁棒再入制导方法 | |
| Zhu et al. | Impact time and angle control guidance independent of time-to-go prediction | |
| CN103853050A (zh) | 一种四旋翼飞行器的pid优化控制方法 | |
| CN109446582B (zh) | 一种考虑地球自转的高精度降阶平稳滑翔动力学建模方法 | |
| Bogdanov et al. | State-dependent Riccati equation control of a small unmanned helicopter | |
| CN105865455A (zh) | 一种利用gps与加速度计计算飞行器姿态角的方法 | |
| CN110908407A (zh) | 一种rlv再入热流率跟踪的改进预测制导方法 | |
| Yu et al. | Design and Implementation of a Hardware‐in‐the‐Loop Simulation System for a Tilt Trirotor UAV | |
| CN105718660B (zh) | 临近空间大范围机动弹道三维包络计算方法 | |
| Biannic et al. | Flare control law design via multi-channel ℋ∞ synthesis: Illustration on a freely available nonlinear aircraft benchmark | |
| CN107796401B (zh) | 跳跃式再入飞行器线性伪谱参数修正横向制导方法 | |
| CN104536448A (zh) | 一种基于Backstepping法的无人机姿态系统控制方法 | |
| Lugo-Cárdenas et al. | The MAV3DSim: A simulation platform for research, education and validation of UAV controllers | |
| Liang et al. | Evolved atmospheric entry corridor with safety factor |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |