CN109002576B - 一种线性高阶比例制导系统脱靶量的幂级数解 - Google Patents

一种线性高阶比例制导系统脱靶量的幂级数解 Download PDF

Info

Publication number
CN109002576B
CN109002576B CN201810593916.0A CN201810593916A CN109002576B CN 109002576 B CN109002576 B CN 109002576B CN 201810593916 A CN201810593916 A CN 201810593916A CN 109002576 B CN109002576 B CN 109002576B
Authority
CN
China
Prior art keywords
power series
order
solution
guidance system
convergence
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.)
Active
Application number
CN201810593916.0A
Other languages
English (en)
Other versions
CN109002576A (zh
Inventor
陈万春
赵石磊
赫泰龙
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Beihang University
Original Assignee
Beihang University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Beihang University filed Critical Beihang University
Priority to CN201810593916.0A priority Critical patent/CN109002576B/zh
Publication of CN109002576A publication Critical patent/CN109002576A/zh
Application granted granted Critical
Publication of CN109002576B publication Critical patent/CN109002576B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00Computer-aided design [CAD]
    • G06F30/20Design optimisation, verification or simulation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Control Of Position, Course, Altitude, Or Attitude Of Moving Bodies (AREA)

Abstract

本发明公开一种线性高阶比例制导系统脱靶量的幂级数解,步骤一:高阶线性比例制导系统建模。步骤二:求解伴随系统的微分方程的幂级数解;包括求解幂级数解系数递推公式和幂级数收敛半径。步骤三:选择合适指数项衰减常数k。本发明优点在于:(1)推导了一般高阶比例制导系统脱靶量的幂级数解系数的递推关系;对于不同有效导引比、不同形式的高阶制导系统,幂级数解系数递推求解过程具有一致性。(2)求出了幂级数的收敛半径并给出了指数型衰减常数的选取方案;利用收敛幂级数的部分和可以得到脱靶量解析的、形式统一的逼近公式。(3)给出的幂级数解是脱靶量一种精确的解的形式,可用来研究脱靶量的性质以及寻找某些特殊条件下的解析解。

Description

一种线性高阶比例制导系统脱靶量的幂级数解
技术领域
本发明提供了一种线性高阶比例制导系统脱靶量的幂级数解,属于航天技术、武器技术领域。
背景技术
比例导引是最经典的制导律,由于其简洁、有效和易于物理实现,目前世界上战术导弹几乎都采用比例导引制导。导弹通常包含导引头、弹体环节、自动驾驶仪等系统,数学上最简化的导弹制导系统也要三阶微分方程来描述,很难得到这些高阶微分方程的解析解(有限项初等函数显示表示),所以一般采用计算机数值仿真的方法来研究比例导引制导系统。
脱靶量是衡量导弹拦截目标或者目标逃逸策略的最重要最直观的性能指标,也是设计分析制导系统的关键核心。利用伴随法进行数值仿真是评估导弹系统设计和求解脱靶量的通用方法。对于线性一阶比例导引制导系统,当有效导引比为正整数时,利用伴随方程可直接得到脱靶量的解析解;但是对于一般的高阶制导系统,并不存在脱靶量的解析解。
幂级数法是求解常微分方程组的有效手段,通常先假设微分方程的解为系数待定的收敛的幂级数,然后将其代入微分方程,利用幂级数恒等条件,得到待定系数序列的递推关系;通常还要求出该幂级数的收敛半径,以确定幂级数解的适用区间。幂级数法在求解特殊线性微分方程、非线性微分方程和实际工程问题中都有重要应用。但是利用幂级数法来研究比例制导系统的工作还很少。
发明内容
本发明的目的是提供一种线性高阶比例制导系统脱靶量的幂级数解,以填补现有技术中利用幂级数法来研究比例制导系统的空白。
本发明首先建立通用的线性高阶比例导引制导系统模型,由伴随构造得到相应的伴随系统,并对伴随系统的状态量和输出量(即脱靶量)进行无量纲化或归一化。假设伴随系统的状态量和脱靶量为幂级数和指数函数乘积的形式,得到幂级数系数的递推关系,并求出了幂级数解的收敛半径。由于脱靶量幂级数解是无穷级数表示的精确解,但实际应用中只能取级数的部分和来近似计算,还分析了参数对一般高阶制导系统幂级数解收敛速度的影响,并给出了参数的选择方案,从而得到了脱靶量的幂级数解。
本发明为一种线性高阶比例制导系统脱靶量的幂级数解,它包括以下三个步骤:
步骤一:线性高阶比例制导系统建模;
考虑一般的线化比例导引制导系统,研究目标阶跃机动和导弹初始瞄准误差角引起的脱靶量。伴随系统的微分方程为:
Figure BDA0001691516670000021
其中,N为有效导引比,t表示剩余飞行时间或总的飞行时间;z1、z2、zu和ζ是伴随系统的状态量,其初值分别为z1(0)=1、z2(0)=0、zu(0)=0和ζ(0)=0,其中z1的初值为1是由伴随系统的脉冲输入转化而来。伴随系统的输出量,即脱靶量为
Figure BDA0001691516670000022
其中,nT为目标机动水平大小,VM为导弹速度,θHE为导弹初始瞄准误差角;MnT表示目标阶跃机动引起的脱靶量,MHE表示导弹初始瞄准误差角引起的脱靶量。
G(s)为表征制导系统动态特性的稳定传递函数,包含导引头动力学、噪声滤波、飞控系统等环节;通常G(s)可以表示成如下形式
Figure BDA0001691516670000023
其中,Q为制导系统阶数,T为参考时间常数或总制导系统时间常数,λq(q=0,1,...Q)是多项式系数。
进一步,将伴随系统的微分方程进行无量纲化,以得到归一化的脱靶量便于应用。引入如下无量纲变量,将伴随系统的状态量、脱靶量以及时间变量转化为无量纲和归一化变量:
Figure BDA0001691516670000031
Figure BDA0001691516670000032
传递函数替换为:
Figure BDA0001691516670000033
因此,为了简化表达式,如无特别说明,后文中无量纲或归一化的变量仍使用原变量符号,方程的求解和结果讨论都是针对无量纲化的变量。
步骤二:求解伴随系统的微分方程的幂级数解;包括求解幂级数解系数递推公式和幂级数收敛半径。
1.幂级数解的系数递推公式
假设无量纲方程有如下形式的幂级数解
Figure BDA0001691516670000034
其中,an、bn、cn、dn分别为各级数的待定系数,e表示自然指数,参数k表示指数项衰减常数,可用来调节幂级数解整体收敛速度。注意ζ的级数解中含有tQ项,这是因为关于ζ动态是由传递函数G(s)来描述的,等价于如下微分方程
Figure BDA0001691516670000041
其中ζ(q)表示变量ζ的q阶导数。
利用关于时间多项式各次幂的系数相等,并结合伴随系统的微分方程的状态初值,可以得到如下递推关系
Figure BDA0001691516670000042
当n≥1时,
Figure BDA0001691516670000043
其中
Figure BDA0001691516670000044
Figure BDA0001691516670000045
这里Bn和Pn都是中间变量,用于简化表达式书写;A和C及其上下标表示排列数和组合数。cn和dn分别是归一化脱靶量MHE和MnT幂级数的系数。
2.幂级数解的收敛半径
首先求k=0时,幂级数
Figure BDA0001691516670000046
的收敛半径,其中an,0代表k=0时幂级数系数an的值。
记状态向量
Figure BDA0001691516670000051
则伴随系统的微分方程可以写成如下状态空间描述
Figure BDA0001691516670000052
其中
Figure BDA0001691516670000053
式中,OQ×1,OQ×Q,O1×Q分别代表Q行1列,Q行Q列,1行Q列的零矩阵。
Figure BDA0001691516670000054
以及状态初值
X(0)=[0 0 0 ... 0 1]T
显然,矩阵R的特征值只有0,任何正整数都不是R的特征值;函数矩阵A(t)是常值矩阵,在t=0处是解析的,而且其幂级数展开收敛半径为无穷大,则可得,上述微分方程的解可以表示为收敛半径为无穷大的幂级数
Figure BDA0001691516670000055
其中Xn是向量值级数系数,维数与X相同。,利用幂级数恒等条件,可以得到递推关系
Figure BDA0001691516670000056
式中,I是大小为(Q+1)×(Q+1)的单位矩阵。
由此系数序列Xn唯一确定。向量Xn中最后一个分量是状态z1幂级数的系数,通过递推关系消掉Xn其它分量,可以验证最后一个分量序列与幂级数递推关系在k=0时所确定的序列an,0相同,由此可得k=0时z1的级数
Figure BDA0001691516670000057
是收敛的且收敛半径无穷大。
其次,证明k为任意数时,幂级数是收敛的。设当
Figure BDA0001691516670000061
为任意数时,幂级数系数序列记为
Figure BDA0001691516670000062
Figure BDA0001691516670000063
级数
Figure BDA0001691516670000064
正是k=0时收敛的级数
Figure BDA0001691516670000065
与ekt的麦克劳林级数的柯西积:
Figure BDA0001691516670000066
级数
Figure BDA0001691516670000067
Figure BDA0001691516670000068
同理。所以,当
Figure BDA0001691516670000069
为任意数时,幂级数仍然收敛。事实上,e-kt麦克劳林级数收敛半径为无穷大,只要k取某一值幂级数解收敛,则k取其它值时相应的幂级数解仍收敛,且收敛半径相同。
步骤三:选择合适指数项衰减常数k;
前两个步骤得到了伴随系统的微分方程的幂级数解,解中包含一个指数项,其随时间的衰减速度由参数k决定,为了使幂级数解能够尽快收敛,需要进一步分析参数k对幂级数解收敛速度的影响,给出选取参数k的方案。
引入用来衡量幂级数解的收敛速度指标变量ncr,其意义是使得部分和逼近误差余式Rn小于指定精度ε的索引变量n的最小值,换句话说,至少需要ncr+1项部分和,才会使得逼近精确解的误差小于ε;ncr与k有关,我们将ncr表示成关于k的函数形式
Figure BDA00016915166700000610
从定义可以看出ncr越小越好,意味着收敛速度越快。使得ncr取值最小的k是最优的,记为
Figure BDA00016915166700000611
同时考虑幂级数精度、收敛速度和递推关系的简化,可以按照如下方案选取参数k:一阶制导系统选取k=1,欠阻尼二阶制导系统选取k=ξ/β,Q阶二项式系统选取k=Q,对于一般的高阶系统可选取
Figure BDA00016915166700000612
上述各式中,αi(i=1,2,...Q1)为一阶环节特征参数,ξj(j=1,2,...Q2)和βj(j=1,2,...Q2)为二阶环节特征参数。
而相应的部分和项数可由式(15)确定。
本发明的优点在于:
(1)推导了一般比例高阶制导系统脱靶量的幂级数解系数的递推关系;对于不同有效导引比、不同形式的高阶制导系统,幂级数解系数递推求解过程具有一致性。
(2)求出了幂级数的收敛半径并给出了指数型衰减常数的选取方案;利用收敛幂级数的部分和可以得到脱靶量解析的、形式统一的逼近公式。
(3)给出的幂级数解是脱靶量一种精确的解的形式,可以用来研究脱靶量的性质以及寻找某些特殊条件下的解析解。
附图说明
图1是本发明流程图。
图2是线性比例制导系统及其伴随系统框图。
图3a是k=3时含有一个二阶环节的五阶制导系统目标阶跃机动脱靶量的伴随仿真结果与幂级数解部分和对比图。
图3b是k=7时含有一个二阶环节的五阶制导系统目标阶跃机动脱靶量的伴随仿真结果与幂级数解部分和对比图。
图3c是k=10时含有一个二阶环节的五阶制导系统目标阶跃机动脱靶量的伴随仿真结果与幂级数解部分和对比图。
图4a是t=10,ε=10-4的一阶制导系统,不同的有效导引比N得到ncr(k)的曲线图。
图4b是N=4,t=10,ε=10-4的欠阻尼二阶制导系统,不同的阻尼比ξ得到ncr(k)的曲线图。
图4c是N=4,t=10,ε=10-4的欠阻尼二阶制导系统,不同的系统阶数Q得到ncr(k)的曲线图。
图4d是N=4,t=10,ε=10-4的一般的五阶制导系统,时间常数分布不同的传递函数G(s)得到ncr(k)的曲线。
图4e是N=4,t=10的含有一个二阶环节的五阶制导系统,不同误差精度ε得到ncr(k)的曲线。
图4f是N=4,t=10的含有一个二阶环节的五阶制导系统,在不同时间t处得到ncr(k)的曲线。
图5a是二项式形式的五阶制导系统目标阶跃机动脱靶量的伴随仿真结果与幂级数解部分和对比图。
图5b是二项式形式的五阶制导系统初始瞄准误差脱靶量的伴随仿真结果与幂级数解部分和对比图。
图5c是各时间常数不等的一阶环节乘积形式的五阶制导系统目标阶跃机动脱靶量的伴随仿真结果与幂级数解部分和对比图。
图5d是各时间常数不等的一阶环节乘积形式的五阶制导系统初始瞄准误差脱靶量的伴随仿真结果与幂级数解部分和对比图。
图5e是含有一个二阶环节的五阶制导系统目标阶跃机动脱靶量的伴随仿真结果与幂级数解部分和对比图。
图5f是含有一个二阶环节的五阶制导系统初始瞄准误差脱靶量的伴随仿真结果与幂级数解部分和对比图。
上述图中,涉及到的符号、代号说明如下:
N为有效导引比,伴随系统中t表示剩余飞行时间或总的飞行时间,nT为目标机动水平大小,VM为导弹速度,θHE为导弹初始瞄准误差角;MnT表示目标阶跃机动引起的脱靶量,MHE表示导弹初始瞄准误差角引起的脱靶量。z1、z2、zu和ζ为伴随系统状态。
具体实施方式
下面将结合附图和实施案例对本发明作进一步的详细说明。
本发明为一种线性高阶比例制导系统脱靶量的幂级数解,包括三个步骤,具体流程如图1所示,下面我们具体介绍上述三个步骤。
步骤一:线性高阶比例制导系统建模;考虑一般的线化比例制导系统,研究目标阶跃机动和导弹初始瞄准误差角引起的脱靶量。图2给出了原线化比例制导系统以及其伴随系统的框图,则可得到如下伴随系统的微分方程:
Figure BDA0001691516670000091
其中,N为有效导引比,t表示剩余飞行时间或总的飞行时间;z1、z2、zu和ζ是伴随系统的状态量,其初值分别为z1(0)=1、z2(0)=0、zu(0)=0和ζ(0)=0,其中z1的初值为1是由伴随系统的脉冲输入转化而来。伴随系统的输出量,即脱靶量为
Figure BDA0001691516670000092
其中,nT为目标机动水平大小,VM为导弹速度,θHE为导弹初始瞄准误差角;MHE是由于导弹初始航向误差造成的脱靶量,MnT是由于目标造成的脱靶量。
G(s)表示一般的稳定传递函数,包含导引头动力学、噪声滤波、飞控系统等环节;通常G(s)可以表示成如下形式
Figure BDA0001691516670000093
即G(s)是由Q1个一阶环节和Q2个二阶环节组成,αi(i=1,2,...Q1)为一阶环节特征参数,ξj(j=1,2,...Q2)和βj(j=1,2,...Q2)为二阶环节特征参数,均为无量纲的正数;不失一般性,令
Figure BDA0001691516670000094
T为参考时间常数或总制导系统时间常数,具有时间的量纲。为了推导一般形式的脱靶量的解,将G(s)关于s分母多项式展开得到
Figure BDA0001691516670000095
其中Q=Q1+2Q2;λq(q=1,2,...Q)是多项式系数,由式(3)中的αi、ξj和βj唯一确定,结合式(4)可以得到
λ0=1,λ1=1 (6)
进一步,将伴随系统的微分方程(1)和输出方程(2)进行无量纲化,以得到归一化的脱靶量便于应用。注意到脱靶量具有长度的量纲,由(2)可得z2具有时间的量纲,zu具有平方时间的量纲,z1和ζ无量纲;因此,引入如下无量纲变量,将时间变量、伴随系统的状态量和输出量(脱靶量)转化为无量纲和归一化变量:
Figure BDA0001691516670000101
Figure BDA0001691516670000102
经过简单求导运算可以得出,无量纲的伴随系统状态量关于无量纲的时间的微分方程与(1)相同,只需将传递函数替换为:
Figure BDA0001691516670000103
因此,为了简化表达式,如无特别说明,后文中无量纲或归一化的变量仍使用原变量符号,方程的求解和结果讨论都是针对无量纲化的变量。
步骤二:求解伴随系统的微分方程的幂级数解;包括求解幂级数解的系数递推公式和幂级数解的收敛半径。
1.幂级数解的系数递推公式
一般地,伴随系统的微分方程(1)不存在由有限项初等函数构成的解析解,但是当G(s)为一阶环节且有效导引比N为正整数时,利用Laplace变换可以得到方程(1)解析解(归一化后结果)为
Figure BDA0001691516670000111
注意式(10)中解析解都是e指数与多项式函数乘积的形式,受此启发,当N不是正整数或者制导系统的阶数是高阶时,我们可以求解和探讨无量纲方程(1)如下形式的幂级数解
Figure BDA0001691516670000112
其中an、bn、cn、dn分别为各级数的待定系数,参数k表示指数项衰减常数,可用来调节幂级数解整体收敛速度。注意ζ的级数解中含有tQ项,这是因为关于ζ动态是由传递函数(9)来描述的,等价于如下微分方程
Figure BDA0001691516670000113
其中ζ(q)表示变量ζ的q阶导数。对z1、ζ、z2和zu级数求导可得
Figure BDA0001691516670000121
其中A和C及其上下标分别表示排列数和组合数。
将式(11)及其相应的导数式(13)代入微分方程(1)和(12),等号两侧的指数部分消掉,利用关于时间多项式各次幂的系数相等,并结合伴随系统状态量的初值,可以得到如下递推关系
Figure BDA0001691516670000122
当n≥1时,
Figure BDA0001691516670000123
其中
Figure BDA0001691516670000124
Figure BDA0001691516670000125
这里Bn和Pn都是中间变量,用于简化表达式书写;A和C及其上下标仍然表示排列数和组合数。cn和dn分别是归一化脱靶量MHE和MnT幂级数的系数。
2.幂级数解的收敛半径
首先求k=0时,幂级数解的收敛半径。因为代表脱靶量的级数z2和zu可以由z1积分得到,故只需求解k=0时z1的级数
Figure BDA0001691516670000131
的收敛半径即可,其中an,0代表k=0时幂级数系数an的值。
记状态向量
X=[ζ ζ(1) ζ(2 ) … ζ(Q-1) z1]T (16)
则微分方程(1)和(11)可以写成如下状态空间描述
Figure BDA0001691516670000132
其中
Figure BDA0001691516670000133
式中,OQ×1,OQ×Q,O1×Q分别代表Q行1列,Q行Q列,1行Q列的零矩阵。
Figure BDA0001691516670000134
以及状态初值
X(0)=[0 0 0 … 0 1]T (20)
显然,矩阵R的特征值只有0,任何正整数都不是R的特征值;函数矩阵A(t)是常值矩阵,在t=0处是解析的,而且其幂级数展开收敛半径为无穷大,则可得,微分方程(17)的解可以表示为收敛半径为无穷大的幂级数
Figure BDA0001691516670000135
其中Xn是向量值级数系数,维数与X相同。将(21)代入(17),利用幂级数恒等条件,可以得到递推关系
Figure BDA0001691516670000141
式中,I是大小为(Q+1)×(Q+1)的单位矩阵。
由此系数序列Xn唯一确定。向量Xn中最后一个分量是状态z1幂级数的系数,通过递推关系消掉Xn其它分量,可以验证最后一个分量序列与递推关系(15)在k=0时所确定的序列an,0相同,由此可得k=0时z1的级数
Figure BDA0001691516670000142
是收敛的且收敛半径无穷大。
其次,证明k为任意数时,由递推关系(14)和(15)所确定的幂级数(11)是收敛的。设当
Figure BDA0001691516670000143
k为任意数时,由递推关系(14)和(15)得到的序列记为
Figure BDA0001691516670000144
Figure BDA0001691516670000145
容易验证
Figure BDA0001691516670000146
满足递推关系(14)和(15),再由该递推关系生成的序列是唯一的,所以式(23)恒成立。仔细观察发现级数
Figure BDA0001691516670000147
正是k=0时收敛的级数
Figure BDA0001691516670000148
与ekt的麦克劳林级数的柯西积:
Figure BDA0001691516670000149
级数
Figure BDA00016915166700001410
Figure BDA00016915166700001411
同理。所以,当
Figure BDA00016915166700001412
k为任意数时,幂级数(11)仍然收敛。事实上,e-kt麦克劳林级数收敛半径为无穷大,只要k取某一值幂级数解收敛,则k取其它值时相应的幂级数解仍收敛,且收敛半径相同。
步骤三:选择合适指数项衰减常数k。
前两个步骤得到了伴随系统的微分方程的幂级数解,解中包含一个指数项,其随时间的衰减速度由参数k决定,为了使幂级数解能够尽快收敛,需要进一步分析参数k对幂级数解收敛速度的影响,给出选取参数k的方案。
首先来看一下不同参数k对同一制导系统幂级数解的影响,这里选取含有一个二阶环节的五阶制导系统
Figure BDA0001691516670000151
式中,α1=0.1,α2=0.2,α3=0.56,ξ=0.7,β=0.1。
图3a、b、c给出了k分别取3、7和10时脱靶量幂级数解(部分和Sn)与伴随仿真结果,算例中N=4。
从图中可以看出,k=7时大约需要60项部分和(即S60)就可以很好地逼近精确解,而k=10时大约80项,k=3时则要超过150项,k=7时幂级数解收敛速度要快于k=3和k=10时。这可以从幂级数解zu的假设形式(13)来给出解释,具体地,zu是由常规幂级数和指数函数e-kt乘积组成。当k较大时,指数函数随着t的衰减速度很快,部分和Sn(t)中指数部分占主导(在n较小时),比如k=10时S60(t)的曲线在t=8处就衰减到0附近,导致偏离脱靶量精确解。当k较小时,指数函数随着t的衰减速度要小于脱靶量精确解的实际衰减速度,而多项式在t较大时随着t的增大是发散的,此时部分和Sn(t)中多项式部分在t较大处可能占主导,比如k=3时S80(t)、S100(t)、S120(t)的曲线随着t的增加严重偏离脱靶量精确解。
我们讨论幂级数部分和Sn来逼近脱靶量精确解的误差都是指截断误差,也就是余式Rn,理论上,Rn随着n的增加趋于0。为了进一步准确客观分析,引入用来衡量幂级数收敛速度指标变量ncr,其意义是使得部分和逼近误差余式Rn小于指定精度ε的索引变量n的最小值,换句话说,至少需要ncr+1项部分和,才会使得逼近精确解的误差小于ε;ncr与k有关,我们将ncr表示成关于k的函数形式
Figure BDA0001691516670000153
从定义可以看出ncr越小越好,意味着收敛速度越快。使得ncr取值最小的k,是最优的,记为
Figure BDA0001691516670000152
注意ncr和kopt与整个制导系统的参数和分析条件有关,包括制导系统阶数、时间常数分布、有效导引比N和分析时间t、误差精度ε等。图4a~图4f给出了不同参数和条件下幂级数解收敛速度指标ncr(k)的曲线,为了保证数值精度,采用四精度甚至更高精度浮点数运算,并且用级数前1001项(足够精确了)部分和S1000(t)作为脱靶量精确解来计算余式,即Rn(t)=S1000(t)-Sn(t)。
综上分析,同时考虑幂级数收敛速度和递推关系的简化,可以按照如下方案选取参数k:一阶制导系统选取k=1,欠阻尼二阶制导系统选取k=ξ/β,Q阶二项式系统选取k=Q,对于一般的高阶系统可选取
Figure BDA0001691516670000161
上述各式中,αi(i=1,2,...Q1)为一阶环节特征参数,ξj(j=1,2,...Q2)和βj(j=1,2,...Q2)为二阶环节特征参数。
而相应的部分和项数可由式(26)确定。通常取值k大一些更好,这时不同制导系统收敛速度几乎相同,而且数值运算稳定性更好,避免双精度浮点数运算舍入误差累积发散;另外,k较大时系数序列dn在n充分大时同号,而且dn/dn-1随着n的增大单调趋于0,则对于任意小的正数q,存在正整数L,使得当n>L时,有
Figure BDA0001691516670000162
利用此式我们可以得到简单的余式估计
Figure BDA0001691516670000163
这里t需要满足qt<1。考虑到收敛速度,k取值一般不超过10。
实施案例
为了验证此幂级数解的精度,研究五阶制导系统,将数值仿真得到的脱靶量与其幂级数解进行对比。比例导引有效导引比N=4,总飞行时间tf=10s。
传递函数G(s)分别选取二项式形式:
Figure BDA0001691516670000164
各时间常数不等的一阶环节乘积形式:
Figure BDA0001691516670000171
其中,α1=0.0667,α2=0.133,α3=0.2,α4=0.267,α5=0.333。
以及三个一阶环节和一个二阶环节乘积形式:
Figure BDA0001691516670000172
其中,α1=0.1,α2=0.2,α3=0.56,ξ=0.7,β=0.1。由以上传递函数可知,总时间常数均为T=1s。
对于传递函数G1(s),指数项衰减常数取k=5,目标阶跃机动引起的脱靶量和导弹初始瞄准误差引起的脱靶量的幂级数解部分和Sn(t)与伴随仿真结果的对比分别如如图5a、b所示。随着n增加,部分和Sn(t)都逐渐逼近脱靶量伴随仿真结果(或者精确解),S40已经与伴随结果相差无几。
对于传递函数G2(s),指数项衰减常数取k=15,目标阶跃机动引起的脱靶量和导弹初始瞄准误差引起的脱靶量的幂级数解部分和Sn(t)与伴随仿真结果的对比分别如如图5c、d所示。随着n增加,部分和Sn(t)都逐渐逼近脱靶量伴随仿真结果(或者精确解),S130已经与伴随结果相差无几。
对于传递函数G3(s),指数项衰减常数取k=10,目标阶跃机动引起的脱靶量和导弹初始瞄准误差引起的脱靶量的幂级数解部分和Sn(t)与伴随仿真结果的对比分别如如图5e、f所示。随着n增加,部分和Sn(t)都逐渐逼近脱靶量伴随仿真结果(或者精确解),S90已经与伴随结果相差无几。

Claims (1)

1.一种线性高阶比例制导系统脱靶量的幂级数解,其特征在于:其包括以下三个步骤:
步骤一:线性高阶比例制导系统建模;
伴随系统的微分方程为:
Figure FDA0003236231470000011
其中,N为有效导引比,t表示总的飞行时间;z1、z2、zu和ζ是伴随系统的状态量,
Figure FDA0003236231470000012
Figure FDA0003236231470000013
Figure FDA0003236231470000014
分别为z1、z2和zu对t的导数;z1、z2、zu和ζ的初值分别为z1(0)=1、z2(0)=0、zu(0)=0和ζ(0)=0,其中z1的初值为1是由伴随系统的脉冲输入转化而来;伴随系统的输出量,即脱靶量为:
Figure FDA0003236231470000015
其中,nT为目标机动水平大小,VM为导弹速度,θHE为导弹初始瞄准误差角;MnT表示目标阶跃机动引起的脱靶量,MHE表示导弹初始瞄准误差角引起的脱靶量;G(s)为表征制导系统动态特性的稳定传递函数,包含导引头动力学、噪声滤波以及飞控系统的动态特性;通常G(s)表示成如下形式:
Figure FDA0003236231470000016
其中,Q为制导系统阶数,T为总制导系统时间常数,λq是多项式系数,q=0,1,...,Q表示多项式的每一项的阶数;
为了便于幂级数的推导,进一步将式(1)所示的伴随系统的微分方程转化为无量纲化的微分方程;定义
Figure FDA0003236231470000017
Figure FDA0003236231470000018
分别为伴随系统的状态量z1、z2、zu、ζ及总飞行时间t的无量纲变量,
Figure FDA0003236231470000019
Figure FDA00032362314700000110
分别为MnT和MHE的无量纲变量,他们的表达式分别为:
Figure FDA0003236231470000021
Figure FDA0003236231470000022
经过简单求导运算得出,无量纲化的微分方程与(1)相同,只需将传递函数G(s)替换为无量纲化的传递函数
Figure FDA0003236231470000023
其表达式为:
Figure FDA0003236231470000024
步骤二:求解伴随系统的微分方程的幂级数解;包括求解幂级数解系数递推公式和幂级数收敛半径;
(一)幂级数解的系数递推公式
假设无量纲化的微分方程有如下形式的幂级数解:
Figure FDA0003236231470000025
其中,an、bn、cn、dn分别为各幂级数的待定系数,e表示自然指数,参数k表示指数项衰减常数,用来调节幂级数解整体收敛速度;
Figure FDA0003236231470000026
描述了状态量
Figure FDA0003236231470000027
的动态特性,其等价于如下微分方程:
Figure FDA0003236231470000028
式中
Figure FDA0003236231470000029
表示变量
Figure FDA00032362314700000210
的q阶导数;
Figure FDA00032362314700000211
Figure FDA00032362314700000212
的初始值;
利用关于时间多项式各次幂的系数相等,并结合伴随系统的微分方程的状态初值,得到如下递推关系:
Figure FDA0003236231470000031
当n≥1时,
Figure FDA0003236231470000032
其中,
Figure FDA0003236231470000033
Figure FDA0003236231470000034
这里Bn和Pn都是中间变量,用于简化表达式书写;A和C及其上下标表示排列数和组合数;cn和dn分别是无量纲化脱靶量
Figure FDA0003236231470000035
Figure FDA0003236231470000036
幂级数的系数;
(二)幂级数解的收敛半径
首先,求指数型衰减常数k=0时,幂级数
Figure FDA0003236231470000037
的收敛半径,其中an,0代表指数型衰减常数k=0时幂级数系数an的值;
记状态向量:
Figure FDA0003236231470000038
则伴随系统的无量纲化的微分方程表示为如下矩阵形式:
Figure FDA0003236231470000041
式中
Figure FDA0003236231470000042
为状态向量X对总飞行时间t的导数;矩阵R和A为如下所示的常值矩阵:
Figure FDA0003236231470000043
Figure FDA0003236231470000044
式中,OQ×1,OQ×Q,O1×Q分别代表Q行1列,Q行Q列,1行Q列的零矩阵;状态向量X的初值为:
X(0)=[0 0 0 … 0 1]T
显然,矩阵R的特征值只有0,任何正整数都不是R的特征值;矩阵A在t=0处是解析的,而且其幂级数展开收敛半径为无穷大,则式(11)所示的微分方程的解表示为如下收敛半径为无穷大的幂级数:
Figure FDA0003236231470000045
其中Xn是状态向量X的幂级数的系数,其为维数与X相同的向量;将式(12)代入式(11),利用幂级数恒等条件,得到递推关系:
Figure FDA0003236231470000046
式中,I是大小为(Q+1)×(Q+1)的单位矩阵;因此系数Xn是唯一的;向量Xn的最后一个分量是状态z1幂级数的系数,其与指数型衰减常数k=0时由幂级数递推关系(10)得到的z1幂级数的系数an,0相同,由此得到指数型衰减常数k=0时z1的级数
Figure FDA0003236231470000047
是收敛的且收敛半径无穷大;
其次,证明指数型衰减常数k为任意数时,式(7)所示的幂级数是收敛的;设当指数型衰减常数
Figure FDA0003236231470000051
Figure FDA0003236231470000052
为任意数时,由递推关系(9)和(10)得到的幂级数的系数为
Figure FDA0003236231470000053
Figure FDA0003236231470000054
容易验证:
Figure FDA0003236231470000055
而幂级数
Figure FDA0003236231470000056
正是幂级数
Figure FDA0003236231470000057
Figure FDA0003236231470000058
的麦克劳林级数的柯西积,即
Figure FDA0003236231470000059
幂级数
Figure FDA00032362314700000510
Figure FDA00032362314700000511
同理;所以,当指数型衰减常数
Figure FDA00032362314700000512
Figure FDA00032362314700000513
为任意数时,幂级数
Figure FDA00032362314700000514
Figure FDA00032362314700000515
仍然收敛;
步骤三:选择合适指数项衰减常数k;
定义Sn为式(7)所示幂级数前n项之和,Rn为逼近误差余式,即式(7)所示幂级数与Sn的差;引入用来衡量幂级数解的收敛速度指标变量ncr,其定义为使Rn小于指定精度ε的式(7)所示幂级数中索引变量n的最小值,表示成函数形式为:
Figure FDA00032362314700000517
从定义看出ncr越小越好,意味着收敛速度越快;使得收敛速度指标变量ncr取值最小的指数型衰减常数k是最优的,记为:
Figure FDA00032362314700000516
同时考虑幂级数精度、收敛速度和递推关系的简化,按照如下方案选取指数型衰减常数k:一阶制导系统选取k=1;欠阻尼二阶制导系统选取k=ξ/β,式中ξ和β为欠阻力二阶制导系统的特征参数;Q阶二项式系统选取k=Q;对于高阶系统选取:
Figure FDA0003236231470000061
式(18)中,Q1为一阶环节的个数,αi(i=1,2,...Q1)为一阶环节特征参数;Q2为二阶环节的个数,ξj(j=1,2,...Q2)和βj(j=1,2,...Q2)为二阶环节特征参数。
CN201810593916.0A 2018-06-11 2018-06-11 一种线性高阶比例制导系统脱靶量的幂级数解 Active CN109002576B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810593916.0A CN109002576B (zh) 2018-06-11 2018-06-11 一种线性高阶比例制导系统脱靶量的幂级数解

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810593916.0A CN109002576B (zh) 2018-06-11 2018-06-11 一种线性高阶比例制导系统脱靶量的幂级数解

Publications (2)

Publication Number Publication Date
CN109002576A CN109002576A (zh) 2018-12-14
CN109002576B true CN109002576B (zh) 2021-11-02

Family

ID=64601227

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810593916.0A Active CN109002576B (zh) 2018-06-11 2018-06-11 一种线性高阶比例制导系统脱靶量的幂级数解

Country Status (1)

Country Link
CN (1) CN109002576B (zh)

Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105659793B (zh) * 2009-09-25 2013-06-19 北京航空航天大学 一种基于标控脱靶量概念的防空导弹助推段制导方法
CN103728976A (zh) * 2013-12-30 2014-04-16 北京航空航天大学 一种基于广义标控脱靶量概念的多过程约束和多终端约束末制导律
CN103888100A (zh) * 2014-03-29 2014-06-25 北京航空航天大学 一种基于负熵的非高斯线性随机系统滤波方法
CN104266546A (zh) * 2014-09-22 2015-01-07 哈尔滨工业大学 一种基于视线的有限时间收敛主动防御制导控制方法
CN105549387A (zh) * 2015-12-07 2016-05-04 北京航空航天大学 一种助推段广义标控脱靶量解析制导方法
CN105807781A (zh) * 2014-12-31 2016-07-27 上海新跃仪表厂 一种基于比例导引的空间近距飞越末制导方法
CN106843272A (zh) * 2017-02-28 2017-06-13 北京航空航天大学 一种具有终端速度、弹道倾角和过载约束的显式制导律
CN107179021A (zh) * 2017-06-14 2017-09-19 北京理工大学 一种驾束制导体制下多弹协同高精度制导控制方法
CN108036676A (zh) * 2017-12-04 2018-05-15 北京航空航天大学 一种基于三维再入弹道解析解的全射向自主再入制导方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8933382B2 (en) * 2011-03-31 2015-01-13 Raytheon Company Guidance system and method for missile divert minimization

Patent Citations (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105659793B (zh) * 2009-09-25 2013-06-19 北京航空航天大学 一种基于标控脱靶量概念的防空导弹助推段制导方法
CN103728976A (zh) * 2013-12-30 2014-04-16 北京航空航天大学 一种基于广义标控脱靶量概念的多过程约束和多终端约束末制导律
CN103888100A (zh) * 2014-03-29 2014-06-25 北京航空航天大学 一种基于负熵的非高斯线性随机系统滤波方法
CN104266546A (zh) * 2014-09-22 2015-01-07 哈尔滨工业大学 一种基于视线的有限时间收敛主动防御制导控制方法
CN105807781A (zh) * 2014-12-31 2016-07-27 上海新跃仪表厂 一种基于比例导引的空间近距飞越末制导方法
CN105549387A (zh) * 2015-12-07 2016-05-04 北京航空航天大学 一种助推段广义标控脱靶量解析制导方法
CN106843272A (zh) * 2017-02-28 2017-06-13 北京航空航天大学 一种具有终端速度、弹道倾角和过载约束的显式制导律
CN107179021A (zh) * 2017-06-14 2017-09-19 北京理工大学 一种驾束制导体制下多弹协同高精度制导控制方法
CN108036676A (zh) * 2017-12-04 2018-05-15 北京航空航天大学 一种基于三维再入弹道解析解的全射向自主再入制导方法

Non-Patent Citations (7)

* Cited by examiner, † Cited by third party
Title
A new interpretation of adjoint method in linear time-varying system analysis;Tailong He等;《2017 IEEE International Conference on Cybernetics and Intelligent Systems (CIS) and IEEE Conference on Robotics, Automation and Mechatronics (RAM)》;20180101;第58-63页 *
Guidance law with circular no-fly zone constraint;webbin yu等;《Nonlinear Dynamics》;20140722;第1953–1971页 *
基于零控脱靶量的大气层外超远程拦截制导;陈峰等;《航空学报》;20090925;第30卷(第09期);第1583-1589页 *
基于预测脱靶量的远程拦截速度增益导引;陈峰等;《航空学报》;20081125;第29卷(第6期);第1665-1672页 *
导弹制导精度MATRIXx伴随分析系统;邹晖等;《飞行力学》;20011215;第19卷(第4期);第73-77页 *
扩展弹道成型制导系统脱靶量特性分析;王辉等;《红外与激光工程》;20130525;第42卷(第05期);第1322-1329页 *
考虑一阶驾驶仪动力学的扩展弹道成型制导系统解析研究;王辉等;《系统工程与电子技术》;20140307;第36卷(第03期);第509-518页 *

Also Published As

Publication number Publication date
CN109002576A (zh) 2018-12-14

Similar Documents

Publication Publication Date Title
Brezinski Computational aspects of linear control
CN109254533B (zh) 基于状态积分的梯度-修复算法的高超声速飞行器快速轨迹优化方法
Burstein Numerical methods in multidimensional shocked flows
Shinbrot A least squares curve fitting method with applications to the calculation of stability coefficients from transient-response data
Kim et al. Finite horizon integrated guidance and control for terminal homing in vertical plane
Wilson Monte-Carlo calculations for the lattice gauge theory
Qian et al. LPV/PI control for nonlinear aeroengine system based on guardian maps theory
CN109002576B (zh) 一种线性高阶比例制导系统脱靶量的幂级数解
CN102566427A (zh) 一种飞行器鲁棒控制方法
Biertümpfel et al. Time‐varying robustness analysis of launch vehicles under thrust perturbations
Gruca et al. Approximation of high-order systems by low-order models with delays
Meadows Solution of systems of linear ordinary differential equations with periodic coefficients
Lan et al. Finite difference based iterative learning control with initial state learning for a class of fractional order two‐dimensional continuous‐discrete linear systems
Napoli et al. A new collocation algorithm for solving even-order boundary value problems via a novel matrix method
Vanbeveren et al. On optimal and suboptimal actuator selection strategies
Cutrone et al. Transition prediction in hypersonic regime on complex geometries with RANS-based models
Feyel Evolutionary fixed-structure controller tuning against multiple requirements
Guang-Ren et al. Robust control system design using proportional plus partial derivative state feedback
Sarolea et al. Integral equation method in the theory of liquids
Cao et al. Stochastic stability of sigma-point unscented predictive filter
Drobyshevich et al. An algorithm of solution of parabolic equations with different time-steps in subdomains
Bajodek et al. Instability conditions for reaction-diffusion-ODE systems
Gray et al. Toward the numerical design of non linear feedback systems by Zakian's method of inequalities
Sias et al. Enhanced convergence of integral transform solution of ablation problems
Waghmare et al. Real Time Heat Exchanger Control Using Predictive Functional Control

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