高阶线性比例制导系统脱靶量幂级数解

赫泰龙，陈万春，周浩

北京航空航天大学 宇航学院，北京 100083

比例导引是最经典的制导律，由于其简洁、有效和易于物理实现，目前世界上战术导弹几乎都采用比例导引制导[1-3]。脱靶量是设计分析制导系统的关键性能指标，通过研究线性化比例制导系统，Zarchan[1]提出了零控脱靶量的概念，用于解释比例导引，还用于推导和理解其他先进的制导律[3-6]。直接对制导系统微分方程进行数值仿真是求解脱靶量的通用方法，但是伴随法才是战术导弹制导系统最主要的分析设计工具[1,7-9]。对于一阶线性比例制导系统，当有效导引比为正整数时，利用伴随方程可直接得到脱靶量的解析解[1,10]；但是对于一般的高阶制导系统，并不存在脱靶量的解析解。另外，比例导引伴随系统的初始时刻是该微分方程的正则奇点[11-12]，在伴随仿真时需要在伴随时间上增加一个小量来避免奇异，不过大量的应用算例表明伴随仿真结果是可靠的[1]。

幂级数法在求解特殊线性微分方程、非线性微分方程和实际工程问题中都有重要应用[13-17]。但是利用幂级数法来研究比例导引制导系统的工作还很少，文献[11]给出了一阶制导系统幂级数解的简单示例；Holt[18]研究了一个特殊的三阶比例制导系统，该系统由3个不同带宽的单延迟环节表示，并且得到了脱靶量的幂级数解。文献[11,18]的工作都是针对特殊的制导系统，缺乏通用性,都没有给出幂级数解的收敛性证明；也没有研究幂级数部分和逼近精确解的速度[15,19]、精度和适用区间。

本文针对一般的高阶线性比例导引制导系统，推导得到了幂级数和指数函数e-kt乘积的形式的脱靶量的解，并证明了该幂级数解的收敛性。然后给出了一阶系统和高阶二项式系统幂级数解简化递推关系。考虑实际应用，本文还重点分析了参数k对幂级数解收敛速度的影响，同时对比伴随仿真结果，给出了参数k的选择方法。幂级数解为比例制导系统的设计评估提供了一种新的手段。

1 线性比例导引制导系统及脱靶量

考虑一般的线化比例导引制导系统，研究目标阶跃机动和导弹初始瞄准误差角引起的脱靶量。图1给出了原线化比例导引制导回路以及伴随系统框图，脱靶量则可以由伴随系统一次仿真得到[1]。N为有效导引比，δnT为目标阶跃机动转化而来得脉冲输入，δHE为导弹瞄准误差角初始状态转化而来的脉冲输入，VC为导弹和目标接近速度，nL为导弹活度的加速度，伴随系统中t为剩余飞行时间或总的飞行时间，nT为目标机动水平大小，VM为导弹速度，θHE为导弹初始瞄准误差角，MHE为导弹初始瞄准误差角引起的脱靶量，MnT为目标阶跃机动引起的脱靶量；记伴随系统状态分别为z1、z2、zu和ζ，则可得到伴随系统微分方程为

(1)

图1 线性比例导引制导系统及其伴随系统Fig.1 Linear proportional navigation guidance system and adjoint system

式中:z1的初值为1,是由伴随系统的脉冲输入转化而来。该伴随系统输出得到脱靶量为

(2)

G(s)表示一般的稳定传递函数，可能包含导引头动力学、噪声滤波、飞控系统等环节；通常G(s)可以表示为

(3)

即G(s)是由Q1个一阶环节和Q2个二阶环节组成，αi、ξj和βj为各环节特征参数系数，无量纲的正数；不失一般性，令

(4)

T为参考时间常数或总制导系统时间常数，具有时间的量纲。为了推导一般形式的脱靶量的解，将G(s)关于s分母多项式展开可得到

(5)

其中:Q=Q1+2Q2、λ0，λ1，…，λQ为多项式系数，由式(3)中的αi、ξj和βj唯一确定，结合式(4)可以得到

λ0=1,λ1=1

(6)

进一步，便于应用，将伴随式(1)、式(2)进行无量纲化：

(7)

(8)

经过简单求导运算可以得出，无量纲的伴随状态关于无量纲的时间的微分方程与式(1)相同，只需将传递函数替换为

(9)

因此，为了简化表达式，如无特别说明，后文中无量纲或归一化的变量仍使用原变量符号，方程的求解和结果讨论都是针对无量纲化的变量。

2 伴随方程的幂级数解

2.1 脱靶量幂级数的系数的递推关系

一般地，伴随系统微分式(1)不存在由有限项初等函数构成的解析解，但是当G(s)为一阶环节且有效导引比N为正整数时，利用Laplace变换可以得到式(1)解析解[1]，N等于3、4和5的结果(未进行归一化)如表1所示。

表1 N为正整数时一阶制导系统伴随方程的解析解

注意表1中解析解都是e指数与多项式函数乘积的形式，受此启发，当N不是正整数或者制导系统的阶数是高阶时，可以求解和探讨无量纲式(1)如下形式的幂级数解：

(10)

(11)

(12)

(13)

式中:an、bn、cn、dn分别为各级数的待定系数，参数k为指数项衰减常数，可用来调节幂级数解整体收敛速度，这些幂级数都假设为收敛的。注意ζ的级数解中含有tQ项，这是因为关于ζ动态是由传递函数式(9)来描述的，等价于如下微分方程：

(14)

式中:ζ(q)为变量ζ的q阶导数。对z1、ζ、z2和zu级数求导可得

(15)

式中:A和C及其上下标分别表示排列数和组合数。

将式(10)～式(13)及其相应的导数式(15)代入微分式(1)和式(14)，等号两侧的指数部分消掉，利用关于时间多项式各次幂的系数相等，并结合伴随方程的状态初值，可以得到如下递推关系：

(16)

当n≥1时

(17)

式中:

Bn和Pn都是中间变量，用于简化表达式书写；A和C及其上下标仍然表示排列数和组合数;cn和dn分别为归一化脱靶量MHE和MnT幂级数的系数。

2.2 脱靶量幂级数的收敛半径为无穷大

定理对任意参数k，脱靶量幂级数的收敛半径为无穷大。

证明由于e-kt在t=0处泰勒展开幂级数的收敛半径是无穷大，容易证明如果参数k=0时，由递推关系式(16)和式(17)所确定的幂级数式(10)～式(13)是收敛的，则k为任意数时，由同样递推关系所确定的幂级数仍然收敛，且收敛半径相同。

记状态向量

(18)

则微分方程式(1)和式(14)可以写成如下状态空间描述：

(19)

式中:

(20)

A(t)=

(21)

以及状态初值

(22)

显然，矩阵R的特征值只有0，任何正整数都不是R的特征值；函数矩阵A(t)是常值矩阵，在t=0处是解析的，而且其幂级数展开收敛半径为无穷大，则由文献[14]中定理6.1可得，微分方程式(19)的解可以表示为收敛半径为无穷大的幂级数：

(23)

式中:Xn为向量值级数系数，维数与X相同。将式(23)代入式(19)，利用幂级数恒等条件，可以得到递推关系为

(24)

事实上，t=0是微分式(1)的正则奇点，在t=0处式(1)的解是解析的，也就是在t=0的邻域内，式(1)的解可以表示为收敛的幂级数展开式，而且由于R和A(t)的特殊性，或者从递推式(24)也可以直接得到，幂级数解的收敛半径为无穷大。

为了便于研究脱靶量幂级数解的性质，本文还会用到级数的前n+1项部分和Sn与级数的n+1项以后的余式Rn。例如，当研究目标阶跃机动引起的脱靶量幂级数解式(13)时,有前面已经证明幂级数解收敛半径为无穷大，这意味着在我们关心的归一化的时间区间[0,tf]里，幂级数是一致收敛的，可以用部分和Sn来一致地逼近脱靶量真实的解，而且理论上可以以任意精度逼近，只要n足够大。

(25)

3 脱靶量幂级数解讨论分析

本节将针对一阶环节、高阶二项式环节以及一般的高阶环节等制导系统，研究其在特定参数k时的脱靶量幂级数解，并对比伴随仿真结果进行讨论分析。

3.1 一阶环节制导系统

考虑如下一阶滞后环节的制导系统：

(26)

这里G(s)是归一化后传递函数，可以用来代表简化的飞控系统或噪声滤波等环节；相应地，Q=1,λ0=1,λ1=1，代入一般的高阶系统幂级数解待定系数的递推关系式(16)和式(17)，可得如下递推关系:

(27)

当n≥1时

(28)

注意到bn递推中含有(1-k)，如果取参数k=1，则递推关系可以得到很大程度简化，而且经过简单的数列技巧变换，可以得出级数z2和zu的待定系数cn和dn的递推关系:

(29)

(30)

以及cn和dn的通项表达式分别为

(31)

(32)

显然，当N为大于2的正整数时成立

(33)

所以当N为正整数时，归一化脱靶量z2和zu的幂级数部分就转化为有限项的多项式，从而得到脱靶量的解析解；特别地，当取N=3,4,5时，就可以得到表1中所列出的结果，也就是用幂级数法同样可以得到解析解。

当N为正整数时，脱靶量幂级数解就是有限项的解析解，与伴随仿真结果重合，当N为非正整数时，得不到脱靶量的解析解。所以作为算例，这里给出N=3.5和N=4.5时脱靶量幂级数解部分和Sn与伴随仿真结果的对比，如图2所示。

图中横轴和纵轴分别表示归一化的时间和脱靶量。由于脱靶量解中指数部分随着时间快速衰减，而且主要研究导弹末段自寻的制导回路，所以重点关心归一化的时间在10以内的脱靶量情况。从图中可以看出，一阶制导系统脱靶量幂级数解收敛速度很快，取级数的前7项部分和就可以很好的逼近精确脱靶量，特别是在总飞行时间较小的时候。

参数k表示级数中指数部分衰减常数，对于一阶制导系统传递函数式(26)，其齐次解或脉冲响应输出的衰减常数正是1，仿真结果显示也符合脱靶量随飞行时间衰减的速率，因此选择参数k=1是合理的，而且还直接得到了待定系数cn和dn的通项，以及N为正整数时脱靶量的有限项解析表达式。这里只给出目标阶跃机动引起的脱靶量结果，导弹初始瞄准角误差引起脱靶量幂级数解性质类似，后文算例也只讨论目标阶跃机动脱靶量。

图2 一阶制导系统目标阶跃机动脱靶量的伴随仿真结果与幂级数解部分和Fig.2 Adjoint simulation results and partial sums of power series solutions for single-lag guidance system in step target maneuver

3.2 高阶二项式环节制导系统

考虑如下归一化的Q阶二项式制导系统传递函数：

(34)

是由Q个相同的一阶环节乘积得到，形式简单，常用于高阶制导系统性能的初步分析[1]。对比式(9)可以得到

参考3.1节一阶制导系统参数k的选取，考虑到Q阶二项式传递函数式(34)齐次响应指数衰减常数为Q，选取k=Q符合二项式制导系统脱靶量随飞行时间衰减的速率，而且经过排列组合等运算，可以将通用的待定系数的递推关系式(17)简化，即当n≥1时

(35)

(36)

(37)

进一步利用一些技巧变换，可以得到n≥1时，待定系数cn和dn的等效递推关系为

(38)

(39)

当然，在求解cn的前提下，利用式(35)中的递推关系求解dn计算量更少。

3.3 一般的高阶制导系统

对于一般的高阶制导系统，可以利用通用的待定系数的递推关系式(16)和式(17)，得到脱靶量幂级数；以五阶制导系统为例，传递函数G(s)分别选取二项式形式：

(40)

各时间常数不等的一阶环节乘积形式：

(41)

以及三个一阶环节和一个二阶环节乘积形式：

(42)

式中:G2(s)和G3(s)中各参数取值如表2所示。

参数k=9时，3种时间常数分布不同的五阶制导系统，目标阶跃机动引起的脱靶量幂级数解部分和Sn(t)与伴随仿真结果的对比如图3所示，算例中N=4。从图3中可以看出，随着n增加，部分和Sn都逐渐逼近脱靶量伴随仿真结果(或者精确解)，S70已经与伴随结果相差无几，而且直观上3个五阶制导系统逼近速度(级数收敛速度)差不多。另外，不同的时间常数分布的五阶制导系统脱靶量结果很相近[1]。

表2 G2(s)和G3(s)中各参数取值Table 2 Values of parameters in G2(s) and G3(s)

图3 不同时间常数分布的五阶制导系统目标阶跃机动脱靶量的伴随仿真结果与幂级数解部分和逼近Fig.3 Adjoint simulation results and partial sums approximation of normalized miss distance in target maneuver for fifth-order guidance systems with different configurations

4 参数k对幂级数解收敛速度的影响

本节将进一步分析参数k对幂级数解收敛速度的影响，给出选取参数k的方案。

首先来看一下不同参数k对同一制导系统幂级数解的影响，这里选取前面含有一个二阶环节的五阶制导系统G(s)=G3(s)。图4给出了k分别取3、7和10时脱靶量幂级数解(部分和Sn)与伴随仿真结果，算例中N=4。

图4 含有一个二阶环节的五阶制导系统目标阶跃机动脱靶量的伴随仿真结果与幂级数解部分和逼近Fig.4 Adjoint simulation results and partial sums approximation of normalized miss distance in target maneuver for fifth-order guidance system with a quadratic distribution

从图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。如图4(a)所示，随着n的增加，Sn逼近误差较小的区间逐渐扩大，这与普通函数的幂级数展开性质类似；但是当n增大到约150时，就无法通过增加部分和项数n来进一步改善逼近效果，甚至在t=8附近S150(t)的曲线呈现出近似随机振荡的特征，这是由于数值计算过程中舍入误差[20]的累积导致的，本文所有数值算例默认使用双精度浮点数，对于S150(t)中高阶多项式部分在t较大时求值运算(Horner方法[21])双精度是不够的。

为了进一步准确客观分析，引入用来衡量幂级数收敛速度指标变量ncr，其意义是使得部分和逼近误差余式Rn小于指定精度ε的索引变量n的最小值，换句话说，至少需要ncr+1项部分和，才会使得逼近精确解的误差小于ε；ncr与k有关，将ncr表示成关于k的函数形式[19]

ncr(k)≐minnRn(t)≤ε;k

(43)

从定义可以看出ncr越小越好，意味着收敛速度越快。使得ncr取值最小的k，是最优的，记为

(44)

注意ncr和kopt与整个制导系统的参数和分析条件有关，包括制导系统阶数、时间常数分布、有效导引比N和分析时间t、误差精度ε等。图5给出了不同参数和条件下幂级数解收敛速度指标ncr(k)的曲线，为了保证数值精度，采用四精度甚至更高精度浮点数运算，并且用级数前1 001项(足够精确了)部分和S1000(t)作为脱靶量精确解来计算余式，即Rn(t)=S1000(t)-Sn(t)。

图5(a)考虑高阶二项式制导系统，N=4，t=10，ε=10-4，针对不同的系统阶数Q得到ncr(k)的曲线，可以看出最佳参数kopt与阶数Q有关，但都在Q附近，选取参数k=Q是合理的。当k较小时，二项式系统阶数越高，幂级数解收敛速度越慢；当k较大时，不同阶数制导系统幂级数解收敛速度几乎相同。

图5(b)考虑一般的五阶制导系统，N=4，t=10，ε=10-4，针对时间常数分布不同的传递函数G(s)得到ncr(k)的曲线，可以看出最佳参数kopt与系统时间常数分布有关，分布范围越大(如G2(s))kopt越大，二项式分布kopt最小。当k较小时，时间常数散布越大(如G2(s))，ncr越大，幂级数解收敛速度越慢；当k较大时，不同时间常数分布的制导系统幂级数解收敛速度几乎相同。

图5(c)针对一般的含有一个二阶环节的五阶制导系统G(s)=G3(s)，N=4，t=10，得到不同误差精度ε时ncr(k)曲线，显然精度越高(即ε越小)，ncr越大，逼近所需要的部分和项数越多。不同误差精度ε，ncr(k)曲线形式相同，kopt也几乎不变。

图5 不同参数和条件下幂级数解收敛速度指标ncr(k)的曲线Fig.5 Plots of convergence rate index ncr(k) as a function of k with different system parameters and conditions

图5(d)针对一般的含有一个二阶环节的五阶制导系统G(s)=G3(s),N=4,ε=10-4，在不同时间t处得到ncr(k)的曲线，可以看出t越大，逼近到同样精度ε所需要的部分和项数越多；ncr(k)曲线形式相同，kopt也几乎不变。

综上分析，同时考虑幂级数收敛速度和递推关系的简化，可以按照如下方案选取参数k：一阶制导系统选取k=1，Q阶二项式系统选取k=Q，对于一般的高阶系统式(3)或式(9)可选取

(45)

而相应的部分和项数可由ncr(k)确定。通常k取值大一些更好，这时不同制导系统收敛速度几乎相同，而且数值运算稳定性更好，避免双精度浮点数运算舍入误差累积发散；考虑到收敛速度，k取值一般不超过10。

5 结 论

本文研究了线性比例导引制导系统脱靶量的幂级数解，通过理论论证和大量数值算例表明幂级数解是分析比例导引制导系统新的有效手段。

1) 推导了一般高阶比例制导系统脱靶量的幂级数解系数的递推关系，通用性强。

2) 严格证明了脱靶量幂级数解的收敛半径为无穷大。

3) 针对一阶环节和高阶二项式环节等制导系统，通过选取适当的参数k，得到了幂级数解系数简化的递推关系，并研究了不同参数和条件下幂级数解部分和的逼近性质。

4) 通过数值方法讨论分析了参数k对脱靶量幂级数解收敛速度的影响，给出了选取参数k以及部分和项数的方法，为幂级数解的实际应用奠定了基础。