CN103400043B - 一种基于高精度多项式运算的复杂系统传递函数计算方法 - Google Patents

Abstract

本发明属于运载火箭控制技术领域，具体涉及一种基于高精度多项式运算的复杂系统传递函数计算方法。本发明的方法包括以下步骤：步骤1输入数据前处理；步骤2状态空间矩阵计算；步骤3采用FADEEVA法计算箭体的传递函数；步骤4扩充浮点数据信息，得到箭体加平台加陀螺的传递函数和箭体加平台加陀螺加舵机的传递函数；步骤5得到姿控系统开环传递函数。本发明突破了机器字长的限制保留有效数字，减小了计算舍入误差，提高了处理复杂系统时传递函数的求解精度。

Description

一种基于高精度多项式运算的复杂系统传递函数计算方法

技术领域

本发明属于运载火箭控制技术领域，具体涉及一种基于高精度多项式运算的复杂系统传递函数计算方法。

背景技术

分析运载火箭姿态控制系统的稳定性，目前通常有两种方法：其一是根轨迹法，即考察描述姿态控制系统的微分方程所对应的代数特征根的分布状况，若所有的根具有负实部，则判定整个系统是稳定的；其二是频率法，即通过求系统的开环幅相特性和幅相裕度来判定整个系统的稳定性。由于姿态控制系统中伺服机构的非线性及迟滞特性一般不可能用一个线性系统来描述，这就给用根轨迹法设计姿态控制系统带来了一定困难。因此，频率法是目前姿态控制系统设计中最常用的方法，而这其中最重要的就是动力学模型传递函数计算。

FADEEVA算法是一种常用的动力学模型传递函数计算方法，能够很好解决维数较低、简单的姿控系统传递函数计算问题，其计算精度和计算速度均能够满足设计要求。然而在处理大型捆绑运载火箭传递函数计算问题时，由于状态变量维数高（考虑刚、晃、弹的助推段模型状态变量高于200维），且状态矩阵病态严重（数值尺度相差几十次），受限于计算机字长（C语言双精度数有效数字通常为14～15位）在进行多项式运算时存在较大舍入误差，且随着多项式运算次数增加舍入误差也逐渐累积，并最终导致传递函数计算精度难以满足姿控系统分析设计需求。此外，FADEEVA算法易导致“维数灾难”，频率特性易发生畸变，亟需提出一套新的高精度传递函数计算方法。

发明内容

本发明需要解决的技术问题为：现有技术中的动力学模型传递函数计算方法在处理复杂系统时计算精度差，易导致“维数灾难”，频率特性易发生畸变。

本发明的技术方案如下所述：

一种基于高精度多项式运算的复杂系统传递函数计算方法，包括以下步骤：步骤1输入数据前处理；步骤2状态空间矩阵计算；步骤3采用FADEEVA法计算箭体的传递函数；步骤4扩充浮点数据信息，进行高精度多项式运算；将惯导平台通道和速率陀螺通道的两个传递函数相加，得到箭体加平台加陀螺的传递函数；将箭体加平台加陀螺的传递函数与舵机的传递函数相乘，得到箭体加平台加陀螺加舵机的传递函数；步骤5将箭体加平台加陀螺加舵机的传递函数与较正网络的传递函数相乘，得到姿控系统开环传递函数。

步骤2包括以下步骤：

火箭姿态控制系统状态方程表述为如下状态方程形式：

$\left\{\begin{array}{c}\stackrel{\·}{x}=\mathrm{Ax}+\mathrm{Bu}\\ y=\mathrm{Cx}+\mathrm{Du}\end{array}\right.$

式中：

x=[x₁,x₂,...,x_m]^T为m维状态变量；

y=[y₁,y₂,...,y_m]^T为m维输出变量；

u=[u₁,u₂,...,u_r]^T为r维输入变量；

A、B、C、D为状态矩阵。

步骤3包括以下步骤：

纯箭体传递函数输入包括：合成摆角芯级发动机摆角和助推器发动机摆角；输出包括：刚体姿态角信号、和加速度信号、敏感元件测量的姿态角信号、和测量得到的加速度信号；

对步骤2所述火箭姿态控制系统状态方程做拉氏变换，得：

$\left\{\begin{array}{c}\mathrm{sx}\left(s\right)=\mathrm{Ax}\left(s\right)+\mathrm{Bu}\left(s\right)\\ y\left(s\right)=\mathrm{Cx}\left(s\right)+\mathrm{Du}\left(s\right)\end{array}\right.,$

y(s)=C(sI-A)^-1Bu(s)+Du(s)；

输出变量到输入变量的传递函数矩阵为：

$G\left(s\right)=\frac{y\left(s\right)}{u\left(s\right)}=C{(\mathrm{sI}-A)}^{-1}B+D;$

式中，矩阵G(s)的第(i,j)个元素为第i个输出变量到第j个输入变量的传递函数：

$g_{\mathrm{ij}}\left(s\right)=\frac{y_i\left(s\right)}{u_j\left(s\right)};$

FADEEVA法将预解矩阵(sI-A)^-1表示为：

${(\mathrm{sI}-A)}^{-1}=\frac{P\left(s\right)}{|\mathrm{sI}-A|}=\frac{P_0{s}^{n-1}+P_1{s}^{n-2}+...+P_{n-1}}{s^n+{\mathrm{\α}}_1{s}^{n-1}+...+{\mathrm{\α}}_n};$

则其分子的系数矩阵和分母的系数按如下FADEEVA递推公式计算：

$\left\{\begin{array}{cc}{P}_0=I,& {\mathrm{\α}}_1=-\mathrm{trA}\\ P_1={\mathrm{AP}}_0+{\mathrm{\α}}_{1}I,& {\mathrm{\α}}_2=-\frac{1}{2}\mathrm{tr}\left({\mathrm{AP}}_1\right)\\ P_2={\mathrm{AP}}_1+{\mathrm{\α}}_{2}I,& {\mathrm{\α}}_3=-\frac{1}{3}\mathrm{tr}\left({\mathrm{AP}}_2\right)\\ .& .\\ .& .\\ .& .\\ P_{n-1}={\mathrm{AP}}_{n-3}+{\mathrm{\α}}_{n-2}I,& {\mathrm{\α}}_{n-1}=-\frac{1}{n-1}\mathrm{tr}\left({\mathrm{AP}}_{n-2}\right)\\ P_{n-1}={\mathrm{AP}}_{n-2}+{\mathrm{\α}}_{n-1}I,& {\mathrm{\α}}_n=-\frac{1}{n}\mathrm{tr}\left({\mathrm{AP}}_{n-1}\right)\end{array}\right.$

其中tr(A)表示A的迹。

步骤4包括以下步骤：

箭体综合传递函数为：

其中，

W_PT(s)为惯性平台动态特性；

W_ST(s)为速率陀螺动态特性；

依次为俯仰通道静态增益和动态增益；

依次为芯一级俯仰通道静态增益和动态增益；

依次为助推器俯仰通道静态增益和动态增益；

W_SFxj(s)、W_SFzt(s)依次为芯一级发动机伺服机构和助推器伺服机构动态特性；

对浮点数据信息进行扩充，实现高精度多项式运算：

针对浮点数据定义如下数据结构形式：

结构成员的定义如下表：

成员名称	sign	pData	length	exp
					成员类型	char	char*	int	int
含义	符号	有效数字字符串	有效数字个数	指数

将惯导平台通道和速率陀螺通道的两个传递函数相加，得到箭体加平台加陀螺的传递函数；将箭体加平台加陀螺的传递函数与舵机的传递函数相乘，得到箭体加平台加陀螺加舵机的传递函数.

步骤5包括以下步骤：

开环系统传递函数如下式所示：

当考虑加表反馈时，开环系统传递函数为：

将箭体加平台加陀螺加舵机的传递函数与较正网络的传递函数相乘，得到姿控系统开环传递函数。

本发明的有益效果为：

（1）本发明的一种基于高精度多项式运算的复杂系统传递函数计算方法，突破了机器字长的限制保留有效数字，减小了计算舍入误差，提高了处理复杂系统时传递函数的求解精度；

（2）本发明的一种基于高精度多项式运算的复杂系统传递函数计算方法，实现了有效数字位数参数化配置机制，保证了算法对不同复杂程度问题的求解适应性；

（3）本发明的一种基于高精度多项式运算的复杂系统传递函数计算方法已成功应用于二代导航卫星发射技术的工程研制中，后续可进一步推广至我国新一代大型捆绑液体运载火箭的研制过程中。

附图说明

图1为本发明的一种基于高精度多项式运算的复杂系统传递函数计算方法流程图。

具体实施方式

下面结合附图和实施例对本发明的一种基于高精度多项式运算的复杂系统传递函数计算方法进行详细说明。

本发明的一种基于高精度多项式运算的复杂系统传递函数计算方法，包括以下步骤：

步骤1输入数据前处理；

步骤2状态空间矩阵计算；

步骤3采用传统FADEEVA法计算箭体的传递函数；

步骤4通过对浮点数据信息进行扩充，改变其算法内部处理机制，重新定义其多项式运算算法，实现高精度多项式运算；采用上述高精度方法将惯导平台通道和速率陀螺通道的两个传递函数相加，得到箭体加平台加陀螺的传递函数；采用高精度方法将箭体加平台加陀螺的传递函数与舵机的传递函数相乘，得到箭体加平台加陀螺加舵机的传递函数；

步骤5采用步骤4所述高精度方法将箭体加平台加陀螺加舵机的传递函数与较正网络的传递函数相乘，得到姿控系统开环传递函数。

具体而言，步骤1输入数据前处理

此步骤为本领域技术人员公知常识。

步骤2状态空间矩阵计算

火箭姿态控制系统状态方程表述为如下状态方程形式：

$\left\{\begin{array}{c}\stackrel{\·}{x}=\mathrm{Ax}+\mathrm{Bu}\\ y=\mathrm{Cx}+\mathrm{Du}\end{array}\right.$

式中：

x=[x₁,x₂,...,x_m]^T为m维状态变量；

y=[y₁,y₂,...,y_m]^T为m维输出变量；

u=[u₁,u₂,...,u_r]^T为r维输入变量；

A、B、C、D为状态矩阵。

步骤3采用传统FADEEVA法计算箭体的传递函数

纯箭体传递函数是以发动机摆角为输入量，以姿态角、角速度和加速度为输出量的传递函数。

其中，输入包括合成摆角芯级发动机摆角和助推器发动机摆角输出包括刚体姿态角信号 和加速度信号敏感元件测量的姿态角信号 和测量得到的加速度信号

对步骤2所述火箭姿态控制系统状态方程做拉氏变换，可得：

$\left\{\begin{array}{c}\mathrm{sx}\left(s\right)=\mathrm{Ax}\left(s\right)+\mathrm{Bu}\left(s\right)\\ y\left(s\right)=\mathrm{Cx}\left(s\right)+\mathrm{Du}\left(s\right)\end{array}\right.$

消去x(s)，可得：

y(s)=C(sI-A)^-1Bu(s)+Du(s)

由此可得输出变量到输入变量的传递函数矩阵

$G\left(s\right)=\frac{y\left(s\right)}{u\left(s\right)}=C{(\mathrm{sI}-A)}^{-1}B+D$

其中G(s)是m×r的矩阵，此时第(i,j)个元素即为第i个输出变量到第j个输入变量的传递函数：

$g_{\mathrm{ij}}\left(s\right)=\frac{y_i\left(s\right)}{u_j\left(s\right)}$

FADEEVA法将预解矩阵(sI-A)^-1表示为

${(\mathrm{sI}-A)}^{-1}=\frac{P\left(s\right)}{|\mathrm{sI}-A|}=\frac{P_0{s}^{n-1}+P_1{s}^{n-2}+...+P_{n-1}}{s^n+{\mathrm{\α}}_1{s}^{n-1}+...+{\mathrm{\α}}_n}$

则其分子的系数矩阵和分母的系数按如下FADEEVA递推公式计算：

$\left\{\begin{array}{cc}{P}_0=I,& {\mathrm{\α}}_1=-\mathrm{trA}\\ P_1={\mathrm{AP}}_0+{\mathrm{\α}}_{1}I,& {\mathrm{\α}}_2=-\frac{1}{2}\mathrm{tr}\left({\mathrm{AP}}_1\right)\\ P_2={\mathrm{AP}}_1+{\mathrm{\α}}_{2}I,& {\mathrm{\α}}_3=-\frac{1}{3}\mathrm{tr}\left({\mathrm{AP}}_2\right)\\ .& .\\ .& .\\ .& .\\ P_{n-1}={\mathrm{AP}}_{n-3}+{\mathrm{\α}}_{n-2}I,& {\mathrm{\α}}_{n-1}=-\frac{1}{n-1}\mathrm{tr}\left({\mathrm{AP}}_{n-2}\right)\\ P_{n-1}={\mathrm{AP}}_{n-2}+{\mathrm{\α}}_{n-1}I,& {\mathrm{\α}}_n=-\frac{1}{n}\mathrm{tr}\left({\mathrm{AP}}_{n-1}\right)\end{array}\right.$

其中tr(A)表示A的迹。

步骤4通过对浮点数据信息进行扩充，改变其算法内部处理机制，重新定义其多项式运算算法，实现高精度多项式运算；采用上述高精度方法将惯导平台通道和速率陀螺通道的两个传递函数相加，得到箭体加平台加陀螺的传递函数；采用高精度方法将箭体加平台加陀螺的传递函数与舵机的传递函数相乘，得到箭体加平台加陀螺加舵机的传递函数

箭体综合传递函数是考虑测量器件的动态特性，考虑动静态增益，考虑伺服机构特性，不考虑校正网络特性的情况下，整个箭体信号综合以后的传递函数，为单输入单输出传递函数。以其波特图为基础进行姿控系统校正网络设计，所以其计算精度直接影响了姿控系统的设计效果。

对于新一代大型运载火箭，助推发动机摆动参与火箭控制，则箭体综合传递函数为：

其中，

W_PT(s)为惯性平台动态特性；

W_ST(s)为速率陀螺动态特性；

依次为俯仰通道静态增益和动态增益；

依次为芯一级俯仰通道静态增益和动态增益；

依次为助推器俯仰通道静态增益和动态增益；

W_SFxj(s)、W_SFzt(s)依次为芯一级发动机伺服机构和助推器伺服机构动态特性。

考虑到在现有计算机环境下双精度浮点类型数据的有效数字为15位，难以保证传统FADEEVA算法在计算大型火箭姿态控制系统传递函数运算过程中的精度，本发明的方法对浮点数据信息进行扩充，改变其算法内部处理机制，重新定义其多项式运算算法，实现高精度多项式运算。

针对浮点数据定义如下数据结构形式：

结构成员的定义如下表：

成员名称	sign	pData	length	exp
					成员类型	char	char*	int	int
含义	符号	有效数字字符串	有效数字个数	指数

该方法的基本思路是对任意一个浮点类型数据进行重新定义，利用字符串格式数据pData最大限度保存其有效数字，避免算法运算时对字长以外的有效数字进行舍入，并重新定义浮点运算的运算法则，在所有运算过程结束后再将结构类型的运算结果转换为传统浮点类型数据，以备姿控系统设计分析使用。该方法虽然在一定程度上增加了计算量，却突破机器字长的限制最大限度的保留了计算过程中的有效数字。

鉴于改进后算法的计算量与有效数字位数直接相关，因此算法中采用有效数字位数参数化可配置的处理方式，可根据需要自行设定有效数字保留的最大位数，以取得计算效率和计算精度之间的平衡，保证了算法对不同问题的适应性。采用上述高精度方法将惯导平台通道和速率陀螺通道的两个传递函数相加，得到箭体加平台加陀螺的传递函数；采用高精度方法将箭体加平台加陀螺的传递函数与舵机的传递函数相乘，得到箭体加平台加陀螺加舵机的传递函数。

步骤5采用步骤4所述高精度方法将箭体加平台加陀螺加舵机的传递函数与较正网络的传递函数相乘，得到姿控系统开环传递函数

开环系统传递函数是考虑测量器件动态特性，考虑动静态增益，考虑伺服机构特性，考虑校正网络特性的情况下，整个开环系统的传递函数，为单输入单输出传递函数。开环系统传递函数判断姿控系统频域性能指标是否满足的依据。

开环系统传递函数如下式所示：

当考虑加表反馈时，开环系统传递函数为：

采用步骤4所述高精度方法将箭体加平台加陀螺加舵机的传递函数与较正网络的传递函数相乘，得到姿控系统开环传递函数。

Claims (3)

1.一种基于高精度多项式运算的复杂系统传递函数计算方法，其特征在于：包括以下步骤：

步骤1输入数据前处理；

步骤2状态空间矩阵计算；

步骤3采用FADEEVA法计算箭体的传递函数；

步骤4扩充浮点数据信息，进行高精度多项式运算；将惯导平台通道和速率陀螺通道的两个传递函数相加，得到箭体加平台加陀螺的传递函数；将箭体加平台加陀螺的传递函数与舵机的传递函数相乘，得到箭体加平台加陀螺加舵机的传递函数；

步骤5将箭体加平台加陀螺加舵机的传递函数与较正网络的传递函数相乘，得到姿控系统开环传递函数；

步骤2包括以下步骤：

火箭姿态控制系统状态方程表述为如下状态方程形式：

$\left\{\begin{array}{c}\stackrel{\·}{x}=Ax+Bu\\ y=Cx+Du\end{array}\right.$

式中：

x＝[x₁,x₂,...,x_m]^T为m维状态变量；

y＝[y₁,y₂,...,y_m]^T为m维输出变量；

u＝[u₁,u₂,...,u_r]^T为r维输入变量；

A、B、C、D为状态矩阵；

步骤3包括以下步骤：

纯箭体传递函数输入包括：合成摆角芯级发动机摆角和助推器发动机摆角输出包括：刚体姿态角信号和加速度信号敏感元件测量的姿态角信号和测量得到的加速度信号

对步骤2所述火箭姿态控制系统状态方程做拉氏变换，得：

$\left\{\begin{array}{c}sx\left(s\right)=Ax\left(s\right)+Bu\left(s\right)\\ y\left(s\right)=Cx\left(s\right)+Du\left(s\right)\end{array}\right.,$

y(s)＝C(sI-A)^-1Bu(s)+Du(s)；

输出变量到输入变量的传递函数矩阵为：

$G\left(s\right)=\frac{y\left(s\right)}{u\left(s\right)}=C{(sI-A)}^{-1}B+D;$

式中，矩阵G(s)的第(i,j)个元素为第i个输出变量到第j个输入变量的传递函数：

$g_{ij}\left(s\right)=\frac{y_i\left(s\right)}{u_j\left(s\right)};$

FADEEVA法将预解矩阵(sI-A)^-1表示为：

${(sI-A)}^{-1}=\frac{P\left(s\right)}{|sI-A|}=\frac{P_0{s}^{n-1}+P_1{s}^{n-2}+...+P_{n-1}}{s^n+{\mathrm{\α}}_1{s}^{n-1}+...+{\mathrm{\α}}_n};$

则其分子的系数矩阵和分母的系数按如下FADEEVA递推公式计算：

$\left\{\begin{array}{cc}{P}_0=I,& {\mathrm{\α}}_1=-trA\\ P_1={\mathrm{AP}}_0+{\mathrm{\α}}_{1}I,& {\mathrm{\α}}_2=-\frac{1}{2}tr\left({\mathrm{AP}}_1\right)\\ P_2={\mathrm{AP}}_1+{\mathrm{\α}}_{2}I,& {\mathrm{\α}}_3=-\frac{1}{3}tr\left({\mathrm{AP}}_2\right)\\ .& .\\ .& .\\ .& .\\ P_{n-1}={\mathrm{AP}}_{n-3}+{\mathrm{\α}}_{n-2}I,& {\mathrm{\α}}_{n-1}=-\frac{1}{n-1}tr\left({\mathrm{AP}}_{n-2}\right)\\ P_{n-1}={\mathrm{AP}}_{n-2}+{\mathrm{\α}}_{n-1}I,& {\mathrm{\α}}_n=-\frac{1}{n}tr\left({\mathrm{AP}}_{n-1}\right)\end{array}\right.$

其中tr(A)表示A的迹。

2.根据权利要求1所述的一种基于高精度多项式运算的复杂系统传递函数计算方法，其特征在于：步骤4包括以下步骤：

箭体综合传递函数为：

其中，

W_PT(s)为惯性平台动态特性；

W_ST(s)为速率陀螺动态特性；

依次为俯仰通道静态增益和动态增益；

依次为芯一级俯仰通道静态增益和动态增益；

依次为助推器俯仰通道静态增益和动态增益；

W_SFxj(s)、W_SFzt(s)依次为芯一级发动机伺服机构和助推器伺服机构动态特性；

对浮点数据信息进行扩充，实现高精度多项式运算：

针对浮点数据定义如下数据结构形式：

结构成员的定义如下：

成员名称：sign，成员类型：char，含义：符号；

成员名称：pData，成员类型：char*，含义：有效数字字符串；

成员名称：length，成员类型：int，含义：有效数字个数；

成员名称：exp，成员类型：int，含义：指数；

将惯导平台通道和速率陀螺通道的两个传递函数相加，得到箭体加平台加陀螺的传递函数；将箭体加平台加陀螺的传递函数与舵机的传递函数相乘，得到箭体加平台加陀螺加舵机的传递函数。

3.根据权利要求2所述的一种基于高精度多项式运算的复杂系统传递函数计算方法，其特征在于：步骤5包括以下步骤：

开环系统传递函数如下式所示：

当考虑加表反馈时，开环系统传递函数为：

将箭体加平台加陀螺加舵机的传递函数与较正网络的传递函数相乘，得到姿控系统开环传递函数。