CN103400043A - 一种基于高精度多项式运算的复杂系统传递函数计算方法 - Google Patents
一种基于高精度多项式运算的复杂系统传递函数计算方法 Download PDFInfo
- Publication number
- CN103400043A CN103400043A CN2013103492643A CN201310349264A CN103400043A CN 103400043 A CN103400043 A CN 103400043A CN 2013103492643 A CN2013103492643 A CN 2013103492643A CN 201310349264 A CN201310349264 A CN 201310349264A CN 103400043 A CN103400043 A CN 103400043A
- Authority
- CN
- China
- Prior art keywords
- transport function
- add
- platform
- gyro
- rocket body
- 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
Images
Abstract
本发明属于运载火箭控制技术领域,具体涉及一种基于高精度多项式运算的复杂系统传递函数计算方法。本发明的方法包括以下步骤:步骤1输入数据前处理;步骤2状态空间矩阵计算;步骤3采用FADEEVA法计算箭体的传递函数;步骤4扩充浮点数据信息,得到箭体加平台加陀螺的传递函数和箭体加平台加陀螺加舵机的传递函数;步骤5得到姿控系统开环传递函数。本发明突破了机器字长的限制保留有效数字,减小了计算舍入误差,提高了处理复杂系统时传递函数的求解精度。
Description
技术领域
本发明属于运载火箭控制技术领域,具体涉及一种基于高精度多项式运算的复杂系统传递函数计算方法。
背景技术
分析运载火箭姿态控制系统的稳定性,目前通常有两种方法:其一是根轨迹法,即考察描述姿态控制系统的微分方程所对应的代数特征根的分布状况,若所有的根具有负实部,则判定整个系统是稳定的;其二是频率法,即通过求系统的开环幅相特性和幅相裕度来判定整个系统的稳定性。由于姿态控制系统中伺服机构的非线性及迟滞特性一般不可能用一个线性系统来描述,这就给用根轨迹法设计姿态控制系统带来了一定困难。因此,频率法是目前姿态控制系统设计中最常用的方法,而这其中最重要的就是动力学模型传递函数计算。
FADEEVA算法是一种常用的动力学模型传递函数计算方法,能够很好解决维数较低、简单的姿控系统传递函数计算问题,其计算精度和计算速度均能够满足设计要求。然而在处理大型捆绑运载火箭传递函数计算问题时,由于状态变量维数高(考虑刚、晃、弹的助推段模型状态变量高于200维),且状态矩阵病态严重(数值尺度相差几十次),受限于计算机字长(C语言双精度数有效数字通常为14~15位)在进行多项式运算时存在较大舍入误差,且随着多项式运算次数增加舍入误差也逐渐累积,并最终导致传递函数计算精度难以满足姿控系统分析设计需求。此外,FADEEVA算法易导致“维数灾难”,频率特性易发生畸变,亟需提出一套新的高精度传递函数计算方法。
发明内容
本发明需要解决的技术问题为:现有技术中的动力学模型传递函数计算方法在处理复杂系统时计算精度差,易导致“维数灾难”,频率特性易发生畸变。
本发明的技术方案如下所述:
一种基于高精度多项式运算的复杂系统传递函数计算方法,包括以下步骤:步骤1输入数据前处理;步骤2状态空间矩阵计算;步骤3采用FADEEVA法计算箭体的传递函数;步骤4扩充浮点数据信息,进行高精度多项式运算;将惯导平台通道和速率陀螺通道的两个传递函数相加,得到箭体加平台加陀螺的传递函数;将箭体加平台加陀螺的传递函数与舵机的传递函数相乘,得到箭体加平台加陀螺加舵机的传递函数;步骤5将箭体加平台加陀螺加舵机的传递函数与较正网络的传递函数相乘,得到姿控系统开环传递函数。
步骤2包括以下步骤:
火箭姿态控制系统状态方程表述为如下状态方程形式:
式中:
x=[x1,x2,...,xm]T为m维状态变量;
y=[y1,y2,...,ym]T为m维输出变量;
u=[u1,u2,...,ur]T为r维输入变量;
A、B、C、D为状态矩阵。
步骤3包括以下步骤:
对步骤2所述火箭姿态控制系统状态方程做拉氏变换,得:
y(s)=C(sI-A)-1Bu(s)+Du(s);
输出变量到输入变量的传递函数矩阵为:
式中,矩阵G(s)的第(i,j)个元素为第i个输出变量到第j个输入变量的传递函数:
FADEEVA法将预解矩阵(sI-A)-1表示为:
则其分子的系数矩阵和分母的系数按如下FADEEVA递推公式计算:
其中tr(A)表示A的迹。
步骤4包括以下步骤:
箭体综合传递函数为:
其中,
WPT(s)为惯性平台动态特性;
WST(s)为速率陀螺动态特性;
WSFxj(s)、WSFzt(s)依次为芯一级发动机伺服机构和助推器伺服机构动态特性;
对浮点数据信息进行扩充,实现高精度多项式运算:
针对浮点数据定义如下数据结构形式:
结构成员的定义如下表:
成员名称 | sign | pData | length | exp |
成员类型 | char | char* | int | int |
含义 | 符号 | 有效数字字符串 | 有效数字个数 | 指数 |
将惯导平台通道和速率陀螺通道的两个传递函数相加,得到箭体加平台加陀螺的传递函数;将箭体加平台加陀螺的传递函数与舵机的传递函数相乘,得到箭体加平台加陀螺加舵机的传递函数.
步骤5包括以下步骤:
开环系统传递函数如下式所示:
当考虑加表反馈时,开环系统传递函数为:
将箭体加平台加陀螺加舵机的传递函数与较正网络的传递函数相乘,得到姿控系统开环传递函数。
本发明的有益效果为:
(1)本发明的一种基于高精度多项式运算的复杂系统传递函数计算方法,突破了机器字长的限制保留有效数字,减小了计算舍入误差,提高了处理复杂系统时传递函数的求解精度;
(2)本发明的一种基于高精度多项式运算的复杂系统传递函数计算方法,实现了有效数字位数参数化配置机制,保证了算法对不同复杂程度问题的求解适应性;
(3)本发明的一种基于高精度多项式运算的复杂系统传递函数计算方法已成功应用于二代导航卫星发射技术的工程研制中,后续可进一步推广至我国新一代大型捆绑液体运载火箭的研制过程中。
附图说明
图1为本发明的一种基于高精度多项式运算的复杂系统传递函数计算方法流程图。
具体实施方式
下面结合附图和实施例对本发明的一种基于高精度多项式运算的复杂系统传递函数计算方法进行详细说明。
本发明的一种基于高精度多项式运算的复杂系统传递函数计算方法,包括以下步骤:
步骤1输入数据前处理;
步骤2状态空间矩阵计算;
步骤3采用传统FADEEVA法计算箭体的传递函数;
步骤4通过对浮点数据信息进行扩充,改变其算法内部处理机制,重新定义其多项式运算算法,实现高精度多项式运算;采用上述高精度方法将惯导平台通道和速率陀螺通道的两个传递函数相加,得到箭体加平台加陀螺的传递函数;采用高精度方法将箭体加平台加陀螺的传递函数与舵机的传递函数相乘,得到箭体加平台加陀螺加舵机的传递函数;
步骤5采用步骤4所述高精度方法将箭体加平台加陀螺加舵机的传递函数与较正网络的传递函数相乘,得到姿控系统开环传递函数。
具体而言,步骤1输入数据前处理
此步骤为本领域技术人员公知常识。
步骤2状态空间矩阵计算
火箭姿态控制系统状态方程表述为如下状态方程形式:
式中:
x=[x1,x2,...,xm]T为m维状态变量;
y=[y1,y2,...,ym]T为m维输出变量;
u=[u1,u2,...,ur]T为r维输入变量;
A、B、C、D为状态矩阵。
步骤3采用传统FADEEVA法计算箭体的传递函数
纯箭体传递函数是以发动机摆角为输入量,以姿态角、角速度和加速度为输出量的传递函数。
对步骤2所述火箭姿态控制系统状态方程做拉氏变换,可得:
消去x(s),可得:
y(s)=C(sI-A)-1Bu(s)+Du(s)
由此可得输出变量到输入变量的传递函数矩阵
其中G(s)是m×r的矩阵,此时第(i,j)个元素即为第i个输出变量到第j个输入变量的传递函数:
FADEEVA法将预解矩阵(sI-A)-1表示为
则其分子的系数矩阵和分母的系数按如下FADEEVA递推公式计算:
其中tr(A)表示A的迹。
步骤4通过对浮点数据信息进行扩充,改变其算法内部处理机制,重新定义其多项式运算算法,实现高精度多项式运算;采用上述高精度方法将惯导平台通道和速率陀螺通道的两个传递函数相加,得到箭体加平台加陀螺的传递函数;采用高精度方法将箭体加平台加陀螺的传递函数与舵机的传递函数相乘,得到箭体加平台加陀螺加舵机的传递函数
箭体综合传递函数是考虑测量器件的动态特性,考虑动静态增益,考虑伺服机构特性,不考虑校正网络特性的情况下,整个箭体信号综合以后的传递函数,为单输入单输出传递函数。以其波特图为基础进行姿控系统校正网络设计,所以其计算精度直接影响了姿控系统的设计效果。
对于新一代大型运载火箭,助推发动机摆动参与火箭控制,则箭体综合传递函数为:
其中,
WPT(s)为惯性平台动态特性;
WST(s)为速率陀螺动态特性;
依次为俯仰通道静态增益和动态增益;
WSFxj(s)、WSFzt(s)依次为芯一级发动机伺服机构和助推器伺服机构动态特性。
考虑到在现有计算机环境下双精度浮点类型数据的有效数字为15位,难以保证传统FADEEVA算法在计算大型火箭姿态控制系统传递函数运算过程中的精度,本发明的方法对浮点数据信息进行扩充,改变其算法内部处理机制,重新定义其多项式运算算法,实现高精度多项式运算。
针对浮点数据定义如下数据结构形式:
结构成员的定义如下表:
成员名称 | sign | pData | length | exp |
成员类型 | char | char* | int | int |
含义 | 符号 | 有效数字字符串 | 有效数字个数 | 指数 |
该方法的基本思路是对任意一个浮点类型数据进行重新定义,利用字符串格式数据pData最大限度保存其有效数字,避免算法运算时对字长以外的有效数字进行舍入,并重新定义浮点运算的运算法则,在所有运算过程结束后再将结构类型的运算结果转换为传统浮点类型数据,以备姿控系统设计分析使用。该方法虽然在一定程度上增加了计算量,却突破机器字长的限制最大限度的保留了计算过程中的有效数字。
鉴于改进后算法的计算量与有效数字位数直接相关,因此算法中采用有效数字位数参数化可配置的处理方式,可根据需要自行设定有效数字保留的最大位数,以取得计算效率和计算精度之间的平衡,保证了算法对不同问题的适应性。采用上述高精度方法将惯导平台通道和速率陀螺通道的两个传递函数相加,得到箭体加平台加陀螺的传递函数;采用高精度方法将箭体加平台加陀螺的传递函数与舵机的传递函数相乘,得到箭体加平台加陀螺加舵机的传递函数。
步骤5采用步骤4所述高精度方法将箭体加平台加陀螺加舵机的传递函数与较正网络的传递函数相乘,得到姿控系统开环传递函数
开环系统传递函数是考虑测量器件动态特性,考虑动静态增益,考虑伺服机构特性,考虑校正网络特性的情况下,整个开环系统的传递函数,为单输入单输出传递函数。开环系统传递函数判断姿控系统频域性能指标是否满足的依据。
开环系统传递函数如下式所示:
当考虑加表反馈时,开环系统传递函数为:
采用步骤4所述高精度方法将箭体加平台加陀螺加舵机的传递函数与较正网络的传递函数相乘,得到姿控系统开环传递函数。
Claims (5)
1.一种基于高精度多项式运算的复杂系统传递函数计算方法,其特征在于:包括以下步骤:
步骤1输入数据前处理;
步骤2状态空间矩阵计算;
步骤3采用FADEEVA法计算箭体的传递函数;
步骤4扩充浮点数据信息,进行高精度多项式运算;将惯导平台通道和速率陀螺通道的两个传递函数相加,得到箭体加平台加陀螺的传递函数;将箭体加平台加陀螺的传递函数与舵机的传递函数相乘,得到箭体加平台加陀螺加舵机的传递函数;
步骤5将箭体加平台加陀螺加舵机的传递函数与较正网络的传递函数相乘,得到姿控系统开环传递函数。
3.根据权利要求2所述的一种基于高精度多项式运算的复杂系统传递函 数计算方法,其特征在于:步骤3包括以下步骤:
对步骤2所述火箭姿态控制系统状态方程做拉氏变换,得:
y(s)=C(sI-A)-1Bu(s)+Du(s);
输出变量到输入变量的传递函数矩阵为:
式中,矩阵G(s)的第(i,j)个元素为第i个输出变量到第j个输入变量的传递函数:
FADEEVA法将预解矩阵(sI-A)-1表示为:
则其分子的系数矩阵和分母的系数按如下FADEEVA递推公式计算:
其中tr(A)表示A的迹。
4.根据权利要求3所述的一种基于高精度多项式运算的复杂系统传递函数计算方法,其特征在于:步骤4包括以下步骤:
箭体综合传递函数为:
其中,
WPT(s)为惯性平台动态特性;
WST(s)为速率陀螺动态特性;
WSFxj(s)、WSFzt(s)依次为芯一级发动机伺服机构和助推器伺服机构动态特性;
对浮点数据信息进行扩充,实现高精度多项式运算:
针对浮点数据定义如下数据结构形式:
结构成员的定义如下表:
将惯导平台通道和速率陀螺通道的两个传递函数相加,得到箭体加平台加陀螺的传递函数;将箭体加平台加陀螺的传递函数与舵机的传递函数相乘,得到箭体加平台加陀螺加舵机的传递函数。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201310349264.3A CN103400043B (zh) | 2013-08-12 | 2013-08-12 | 一种基于高精度多项式运算的复杂系统传递函数计算方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201310349264.3A CN103400043B (zh) | 2013-08-12 | 2013-08-12 | 一种基于高精度多项式运算的复杂系统传递函数计算方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN103400043A true CN103400043A (zh) | 2013-11-20 |
CN103400043B CN103400043B (zh) | 2016-08-31 |
Family
ID=49563670
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201310349264.3A Expired - Fee Related CN103400043B (zh) | 2013-08-12 | 2013-08-12 | 一种基于高精度多项式运算的复杂系统传递函数计算方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN103400043B (zh) |
Cited By (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN105973516A (zh) * | 2015-12-11 | 2016-09-28 | 北京强度环境研究所 | 一种用于识别固体火箭发动机的脉动推力的方法 |
Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US7989743B2 (en) * | 2006-03-07 | 2011-08-02 | Raytheon Company | System and method for attitude control of a flight vehicle using pitch-over thrusters and application to an active protection system |
CN202453698U (zh) * | 2011-11-15 | 2012-09-26 | 北京宇航系统工程研究所 | 一种运载火箭动力测控系统 |
-
2013
- 2013-08-12 CN CN201310349264.3A patent/CN103400043B/zh not_active Expired - Fee Related
Patent Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US7989743B2 (en) * | 2006-03-07 | 2011-08-02 | Raytheon Company | System and method for attitude control of a flight vehicle using pitch-over thrusters and application to an active protection system |
CN202453698U (zh) * | 2011-11-15 | 2012-09-26 | 北京宇航系统工程研究所 | 一种运载火箭动力测控系统 |
Non-Patent Citations (2)
Title |
---|
张予器: "超高精度浮点运算的关键技术研究", 《中国优秀硕士学位论文全文数据库 信息科技辑》, no. 11, 15 November 2006 (2006-11-15) * |
杨云飞,张立洲: "运载火箭控制系统频域分析软件开发", 《计算机仿真》, vol. 23, no. 9, 30 September 2006 (2006-09-30), pages 15 - 18 * |
Cited By (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN105973516A (zh) * | 2015-12-11 | 2016-09-28 | 北京强度环境研究所 | 一种用于识别固体火箭发动机的脉动推力的方法 |
CN105973516B (zh) * | 2015-12-11 | 2019-04-05 | 北京强度环境研究所 | 一种用于识别固体火箭发动机的脉动推力的方法 |
Also Published As
Publication number | Publication date |
---|---|
CN103400043B (zh) | 2016-08-31 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
Antonov et al. | Unscented Kalman filter for vehicle state estimation | |
Silva et al. | Development of reduced-order models for aeroelastic analysis and flutter prediction using the CFL3Dv6. 0 code | |
Niswonger et al. | MODFLOW-NWT, a Newton formulation for MODFLOW-2005 | |
CN103207568B (zh) | 一种抗舵机饱和的船舶航向自适应控制方法 | |
Guruswamy | A review of numerical fluids/structures interface methods for computations using high-fidelity equations | |
CN106778012A (zh) | 一种小天体附着探测下降轨迹优化方法 | |
Forti et al. | Efficient geometrical parametrisation techniques of interfaces for reduced-order modelling: application to fluid–structure interaction coupling problems | |
CN105136145A (zh) | 一种基于卡尔曼滤波的四旋翼无人机姿态数据融合的方法 | |
Zimmermann et al. | Improved extrapolation of steady turbulent aerodynamics using a non-linear POD-based reduced order model | |
CN102980580A (zh) | 基于张量积多胞鲁棒h2滤波的无陀螺卫星姿态确定方法 | |
CN103994698A (zh) | 基于过载与角速度测量的导弹俯仰通道简单滑模控制方法 | |
CN105955031A (zh) | 非线性预测控制的fpga硬件加速控制器及其加速实现方法 | |
CN105203110A (zh) | 一种基于大气阻力模型补偿的低轨卫星轨道预报方法 | |
Jacobson et al. | Flutter Analysis with Stabilized Finite Elements based on the Linearized Frequency-domain Approach | |
CN115496006A (zh) | 一种适用于高超声速飞行器的高精度数值模拟方法 | |
Kadioglu et al. | A fourth-order auxiliary variable projection method for zero-Mach number gas dynamics | |
Pan et al. | Accurate real-time truck simulation via semirecursive formulation and Adams–Bashforth–Moulton algorithm | |
CN112683274A (zh) | 一种基于无迹卡尔曼滤波的无人机组合导航方法和系统 | |
CN104504189A (zh) | 随机激励下大规模结构设计方法 | |
CN106354013A (zh) | 攻角的线性自抗扰控制方法 | |
CN105425589A (zh) | 提高航天器惯性参数辨识精度的输入信号设计方法 | |
CN103400043A (zh) | 一种基于高精度多项式运算的复杂系统传递函数计算方法 | |
Huang et al. | Aeroelastic simulation using CFD/CSD coupling based on precise integration method | |
Haug | An index 0 differential-algebraic equation formulation for multibody dynamics: Nonholonomic constraints | |
CN103425834A (zh) | 一种柔性材料的形变仿真方法和装置 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
C14 | Grant of patent or utility model | ||
GR01 | Patent grant | ||
CF01 | Termination of patent right due to non-payment of annual fee | ||
CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20160831 |