CN109086247B - 基于双时间尺度无迹卡尔曼滤波的系统故障参数估计方法 - Google Patents

基于双时间尺度无迹卡尔曼滤波的系统故障参数估计方法 Download PDF

Info

Publication number
CN109086247B
CN109086247B CN201811093396.3A CN201811093396A CN109086247B CN 109086247 B CN109086247 B CN 109086247B CN 201811093396 A CN201811093396 A CN 201811093396A CN 109086247 B CN109086247 B CN 109086247B
Authority
CN
China
Prior art keywords
system state
time
sampling
variable
observation
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
CN201811093396.3A
Other languages
English (en)
Other versions
CN109086247A (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.)
Hefei University of Technology
Original Assignee
Hefei University of Technology
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 Hefei University of Technology filed Critical Hefei University of Technology
Priority to CN201811093396.3A priority Critical patent/CN109086247B/zh
Publication of CN109086247A publication Critical patent/CN109086247A/zh
Application granted granted Critical
Publication of CN109086247B publication Critical patent/CN109086247B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/17—Function evaluation by approximation methods, e.g. inter- or extrapolation, smoothing, least mean square method
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q10/00—Administration; Management
    • G06Q10/06—Resources, workflows, human or project management; Enterprise or organisation planning; Enterprise or organisation modelling
    • G06Q10/063—Operations research, analysis or management
    • G06Q10/0639—Performance analysis of employees; Performance analysis of enterprise or organisation operations

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • Business, Economics & Management (AREA)
  • Human Resources & Organizations (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Operations Research (AREA)
  • Strategic Management (AREA)
  • Development Economics (AREA)
  • Economics (AREA)
  • Entrepreneurship & Innovation (AREA)
  • Educational Administration (AREA)
  • Marketing (AREA)
  • Game Theory and Decision Science (AREA)
  • Quality & Reliability (AREA)
  • Tourism & Hospitality (AREA)
  • General Business, Economics & Management (AREA)
  • Computing Systems (AREA)
  • Feedback Control In General (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了基于双时间尺度无迹卡尔曼滤波的系统故障参数估计方法，包括：对电动代步车系统的系统状态空间方程进行离散化处理，得到时间离散的系统状态空间方程；将待估计参数添加到系统状态变量中，得到增广后的系统状态空间方程，所述待估计参数为故障集合中的参数；采用具有双时间尺度的无迹卡尔曼滤波算法对待估计参数与原系统状态变量进行联合估计，在宏观尺度下，仅对增广前的原系统状态变量采用无迹卡尔曼滤波进行估计，在微观尺度下，对增广后的系统状态变量采用无迹卡尔曼滤波进行估计。本发明解决了故障诊断实时性差、故障源定位不及时的问题，减少了故障诊断时的计算量，提高了运算效率。

Description

基于双时间尺度无迹卡尔曼滤波的系统故障参数估计方法
技术领域
本发明涉及故障参数估计领域，尤其是基于双时间尺度无迹卡尔曼滤波的系统故障参数估计方法。
背景技术
在世界各国普遍出现人口老龄化问题以及残疾人口不断增加的背景下，专为出行不便的老年人及残疾人士提供出行便利的电动代步车产品，近年来在全球范围内受到了广泛关注。电动代步车的普及也带来了一些安全问题，作为一种专为老年人以及残疾人士设计的代步工具，一旦在道路行驶过程中发生故障，将可能造成严重的后果，因此开发一种实时的电动代步车故障诊断系统迫在眉睫。
目前，针对非线性系统故障参数估计的方法有扩展卡尔曼滤波法、粒子滤波法以及无迹卡尔曼滤波法等方法。但在一些实时性要求较高的场景下，这些传统的方法未必能满足要求，造成了故障诊断实时性差、故障源定位不及时的问题。
发明内容
为了克服上述现有技术的缺陷，本发明提供基于双时间尺度无迹卡尔曼滤波的系统故障参数估计方法，解决了故障诊断实时性差、故障源定位不及时的问题，减少了故障诊断时的计算量，提高了运算效率。
为实现上述目的，本发明采用以下技术方案，包括：
基于双时间尺度无迹卡尔曼滤波的系统故障参数估计方法，包括以下步骤：
S1，对电动代步车系统的连续的系统状态空间方程进行离散化处理，得到时间离散的系统状态空间方程；
S2，将L2-L1个待估计参数添加到电动代步车系统的L1个原系统状态变量中，得到增广后的系统状态空间方程；
其中，所述待估计参数为故障集合中的参数；增广前的系统状态空间方程的维度为L1，增广前的系统状态变量为L1个；增广后的系统状态空间方程的维度为L2，增广后的系统状态变量为L2个；
S3，采用具有双时间尺度的无迹卡尔曼滤波算法对L2-L1个待估计参数与L1个原系统状态变量进行联合估计，得到参数估计值。
步骤S3中，所述双时间尺度包括宏观尺度和微观尺度；在宏观尺度下，仅对增广前的L1个系统状态变量采用无迹卡尔曼滤波进行估计；在微观尺度下，对增广后的L2个系统状态变量采用无迹卡尔曼滤波进行估计。
步骤S1中，所述电动代步车系统的连续的系统状态空间方程包括状态方程和观测方程，如公式(1)所示：
Figure BDA0001804933460000021
其中，t表示实际的连续时间，x(t)表示t时刻下的系统状态变量，即增广前的原系统状态变量；u(t)表示t时刻下的系统输入；W表示过程噪声；
Figure BDA0001804933460000022
表示t时刻下的系统状态变量的一阶导数；y(t)为t时刻下的系统输出，即观测变量；V表示观测噪声；f(·)和g(·)分别表示状态函数和观测函数，即分别对应为状态方程和观测方程；x(·)表示原系统状态变量；y(·)表示观测变量；
采用前向差分法或后向差分法对所述连续的系统状态空间方程进行离散化，得到时间离散的系统状态空间方程，所述时间离散的系统状态空间方程，如公式(2)所示：
Figure BDA0001804933460000023
其中，k表示第k个采样点所对应的时刻，即表示第k个采样时刻；k+1表示第k+1个采样点所对应的采样时刻，即表示第k+1个采样时刻；x(k)表示第k个采样时刻的系统状态变量；u(k)表示第k个采样时刻的系统输入；x(k+1)表示第k+1个采样时刻的系统状态变量；y(k)为第k个采样采样时刻的系统输出。
步骤S2中，所述增广后的系统状态空间方程，如公式(3)所示：
Figure BDA0001804933460000031
其中，x1(k)表示第k个采样时刻的待估计参数，即为电动代步车系统的新添加的系统状态变量；x1(·)表示待估计参数，即新添加的系统状态变量。
步骤S3中，在宏观尺度下，仅对增广前的L1个系统状态变量进行计算，并对系统的协方差矩阵的部分行列、卡尔曼增益矩阵的部分行列以及观测变量的部分行列进行更新，所述部分行列为增广前的L1个系统状态变量所占的维度L1；在微观尺度下，对增广后的L2个系统状态变量进行计算，并对系统的协方差矩阵的全部行列、卡尔曼增益矩阵的全部行列以及观测变量的全部行列进行更新，所述全部行列为增广后的L2个系统状态变量所占的总维度L2；所述协方差矩阵为系统状态变量在相邻的两个采样点所对应的采样时刻之间的协方差，用于描述系统状态变量在相邻的两个采样时刻之间的递推依赖关系；所述卡尔曼增益矩阵为系统在卡尔曼滤波递推过程中从某个采样时刻到下一采样时刻的增益矩阵。
利用无迹卡尔曼滤波算法进行计算并更新系统的协方差矩阵、卡尔曼增益矩阵、观测变量；所述无迹卡尔曼滤波算法，包括以下步骤：
S31，对系统状态变量进行无迹变换，由无迹变换获得一组采样点，以及每个采样点所对应的权值，并将该组采样点称为Sigma点集；
S32，分别计算Sigma点集中的每个采样点的一步预测值；
S33，根据每个采样点所对应的权值和步骤S32中得到每个采样点的一步预测值，进行加权求和的计算，得到系统状态变量的一步预测值，再根据系统状态变量的一步预测值得到系统的协方差矩阵；
S34，对系统状态变量的一步预测值进行无迹变换，产生一组新的采样点，即新的Sigma点集，以及每个新的采样点对应的权值；
S35，分别将由步骤S34得到的每个新的采样点代入观测方程，分别得到每个新的采样点的观测变量的预测量；
S36，根据步骤S35得到每个新的采样点的观测变量的预测量，进行加权求和计算，得到观测变量的均值；根据观测变量的均值得到观测变量的协方差
Figure BDA0001804933460000041
系统状态变量和观测变量之间的协方差
Figure BDA0001804933460000042
S37，计算并更新系统的卡尔曼增益矩阵；
S38，计算并更新系统状态变量和系统的协方差矩阵。
本发明的优点在于：
(1)本发明采用的双时间尺度无迹卡尔曼滤波算法是一个针对非线性系统的算法，以无迹卡尔曼滤波算法为基本框架，对于一步预测方程采用了无迹变换来处理均值和协方差的非线性传递问题。它对非线性函数的概率密度分布进行近似，用一系列确定样本来逼近状态的后验概率密度，而不是对非线性函数进行近似，从而避免了线性化误差的引入，因此对于非线性分布的统计量有较高的精度。
(2)本发明采用了将原系统状态变量和待估计参数进行联合估计的方法，能够精确地估计出故障集合中参数值的大小，从而实现故障源的准确定位。
(3)本发明基于双时间尺度，由于系统的待估计参数的变化速率远远小于系统的状态变量的变化速率，采用了基于无迹卡尔曼滤波算法框架的宏观、微观的双时间尺度算法，降低了部分情形下矩阵的运算维度，从而改进了计算效率，提升了故障参数估计的实时性。
附图说明
图1为本发明的方法流程图。
图2为利用传统的无迹卡尔曼滤波的故障参数估计方法进行的参数估计结果。
图3(1)为利用本发明的基于双时间尺度无迹卡尔曼滤波的系统故障参数估计方法并以宏观步长为2进行的参数估计结果。
图3(2)为利用本发明的基于双时间尺度无迹卡尔曼滤波的系统故障参数估计方法并以宏观步长为5进行的参数估计结果。
图3(3)为利用本发明的基于双时间尺度无迹卡尔曼滤波的系统故障参数估计方法并以宏观步长为10进行的参数估计结果。
图3(4)为利用本发明的基于双时间尺度无迹卡尔曼滤波的系统故障参数估计方法并以宏观步长为15进行的参数估计结果。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
由图1所示，基于双时间尺度无迹卡尔曼滤波的系统故障参数估计方法，包括以下步骤：
S1，对电动代步车系统的连续的系统状态空间方程进行离散化处理，得到时间离散的系统状态空间方程；
S2，将L2-L1个待估计参数添加到电动代步车系统的L1个原系统状态变量中，得到增广后的系统状态空间方程；
其中，所述待估计参数为故障诊断中可能会发生的故障的参数，即故障集合中的参数；增广前的系统状态空间方程的维度为L1，增广前的系统状态变量为L1个；增广后的系统状态空间方程的维度为L2，增广后的系统状态变量为L2个；
S3，采用具有双时间尺度的无迹卡尔曼滤波算法对L2-L1个待估计参数与L1个原系统状态变量进行联合估计，得到参数估计值。
步骤S1中，所述电动代步车系统的连续的系统状态空间方程包括状态方程和观测方程，如下公式所示：
Figure BDA0001804933460000051
其中，t表示实际时间，x(t)表示t时刻下的系统状态变量，即增广前的原系统状态变量；u(t)表示t时刻下的系统输入；W表示过程噪声；
Figure BDA0001804933460000052
表示t时刻下的系统状态变量的一阶导数；y(t)为t时刻下的系统输出，即观测变量；V表示观测噪声；f(·)和g(·)分别表示状态函数和观测函数，即分别对应为状态方程和观测方程；x(·)表示原系统状态变量；y(·)表示观测变量；
采用前向差分法或后向差分法对所述连续的系统状态空间方程进行离散化，得到时间离散的系统状态空间方程，所述时间离散的系统状态空间方程，如下公式所示：
Figure BDA0001804933460000061
其中，k表示第k个采样点所对应的时刻，即表示第k个采样时刻；k+1表示第k+1个采样点所对应的采样时刻，即表示第k+1个采样时刻；x(k)表示第k个采样时刻的系统状态变量；u(k)表示第k个采样时刻的系统输入；x(k+1)表示第k+1个采样时刻的系统状态变量；y(k)为第k个采样采样时刻的系统输出。
本实施例中，所述电动代步车系统状态空间方程包括状态方程和观测方程，分别如下公式(1-1)和公式(1-2)所示：
Figure BDA0001804933460000062
Figure BDA0001804933460000063
其中，系统状态变量
Figure BDA0001804933460000064
观测变量
Figure BDA0001804933460000065
系统状态变量中θ1，
Figure BDA0001804933460000066
s，
Figure BDA0001804933460000067
θ2，
Figure BDA0001804933460000068
分别表示后轮的角位移、后轮的角速度、后轮的位移、车身的线速度、前轮的角位移、前轮的角速度；观测变量中
Figure BDA0001804933460000069
表示后轮的角速度，
Figure BDA00018049334600000610
表示车身的线速度，
Figure BDA00018049334600000611
表示前轮处的角速度，分别由安装在对应位置的传感器测得。uin表示输入电压；k1表示电压到电流的转换比；k2表示电流到转矩的转换比；k3表示转矩到角速度的转换比；k4表示车轮半径；J1表示电机的转动惯量；J2表示后轮的转动惯量；J3表示前轮的转动惯量；kf和Fu分别表示电机处的粘性摩擦和库伦摩擦；kf1和Fu1分别表示后轮处的粘性摩擦和库伦摩擦；kf2和Fu2分别表示前轮处的粘性摩擦和库伦摩擦；C1和C2分别表示后轮到车身以及车身到前轮的传动轴的刚度；m表示车重。
采用前向差分法或后向差分法对所述连续的系统状态空间方程进行离散化，得到时间离散的系统状态空间方程，所述时间离散的系统状态空间方程的状态方程和观测方程，分别如下公式(2-1)和公式(2-2)所示：
Figure BDA0001804933460000071
Figure BDA0001804933460000072
其中，ts表示采样时间间隔；k表示离散状态下的第k个采样点所对应的采样时刻，即表示第k个采样时刻。
步骤S2中，所述增广后的系统状态空间方程，如下公式所示：
Figure BDA0001804933460000073
其中，x1(k)表示第k个采样时刻的待估计参数，即为电动代步车系统的新添加的系统状态变量；x1(·)表示待估计参数，即新添加的系统状态变量。
若待估计参数包含在原系统状态空间方程中，即属于原系统状态变量中，则需要将新添加的待估计参数代替原系统状态空间方程中对应的系统状态变量。
本实施例中，以前轮传感器发生故障为例，β3为前轮处传感器的有效因子，选取x7＝β3添加到系统状态变量中，β3为待估计参数，即为新添加的系统状态变量；由于待估计参数β3不属于系统的原状态变量中，故在系统状态方程中不做替换，新添加的系统状态变量x7＝β3影响了前轮传感器的测量值，故在观测方程中乘以该新添加的系统状态变量，由此，得到增广后的系统状态空间方程的状态方程和观测方程，分别如下公式(3-1)和公式(3-2)所示：
Figure BDA0001804933460000081
Figure BDA0001804933460000082
步骤S3中，所述双时间尺度包括宏观尺度和微观尺度，所述宏观尺度下的采样时间间隔为大间隔，所述微观尺度下的采样时间间隔为小间隔。由于系统的待估计参数的变化速率远远小于系统的状态变量的变化速率，故在宏观尺度下，仅对增广前的L1个系统状态变量系统采用无迹卡尔曼滤波进行估计，对系统的协方差矩阵的部分行列、卡尔曼增益矩阵的部分行列以及观测变量的部分行列进行更新，所述部分行列为增广前的L1个原系统状态变量所占的维度L1，本实施例中，增广前的系统状态变量为X＝[x1,x2,x3,x4,x5,x6]；在微观尺度下，对增广后的L2个系统状态变量采用无迹卡尔曼滤波进行估计，即对增广前的L1个系统状态变量和L2-L1个待估计参数均采用无迹卡尔曼滤波进行估计，对系统的协方差矩阵的全部行列、卡尔曼增益矩阵的全部行列以及观测变量的全部行列进行更新，所述全部行列为增广后的L2个系统状态变量所占的总维度L2，本实施例中，增广后的系统状态变量为X＝[x1,x2,x3,x4,x5,x6,x7]。
所述协方差矩阵为系统状态变量在相邻的两个采样点所对应的采样时刻之间的协方差，用于描述系统状态变量在相邻的两个采样时刻之间的递推依赖关系；所述卡尔曼增益矩阵为系统在卡尔曼滤波递推过程中从某个采样时刻到下一采样时刻的增益矩阵。
所述无迹卡尔曼滤波为摒弃了对非线性函数进行线性化的做法，采用卡尔曼线性滤波框架，对于一步预测方程，使用无迹变换UT来处理均值和协方差的非线性传递问题。
无迹变换UT的实现方法为：在原状态分布中按某一规则选取一些采样点，使这些采样点的均值和协方差等于原状态分布的均值和协方差；将这些点代入非线性函数中，相应得到非线性函数值点集，通过这些点集求取变换后的均值和方差。这样得到的非线性变换后的均值和协方差精度最少具有2阶精度，对于高斯分布，可达到3阶精度。其采样点的选择是基于先验均值和先验协方差矩阵的平方根的相关列实现的。
无迹变换UT的基本原理为：若非线性变换y＝f(x)，状态x为n维随机变量，并且已知其均值
Figure BDA0001804933460000091
和方差P，则可通过无迹变换UT的方式得到2n+1个Sigma点即采样点，以及每个采样点相应的权值ω，以计算y的统计特征，包括以下步骤：
S11，计算2n+1个Sigma点，即采样点，这里的n指的是状态x的维数；
Figure BDA0001804933460000092
其中，
Figure BDA0001804933460000093
Figure BDA0001804933460000094
表示矩阵方根的第i列；
S12，计算这些采样点相应的权值。
Figure BDA0001804933460000095
其中，下标m表示均值，下标c表示协方差，上标表示第几个采样点；参数λ是一个缩放比例系数，用以降低总的预测误差，λ＝a2(n+κ)-n；α的选取控制了采样点的分布状态；κ为待选参数，其具体值虽然没有界限，但通常应确保矩阵(n+λ)P为半正定矩阵；待选参数β是一个非负的权系数，β≥0，用以合并方程中高阶项的动差，用以把高阶项的影响包括在内。
本实施例中，对于公式(2)或公式(3)所描述的系统的状态空间方程，其中，公式(2)为增广前的状态状态空间方程，公式(3)为增光后的系统状态空间方程，在宏观尺度下，根据公式(2)进行系统更新的相关计算，在微观尺度下，根据公式(3)进行系统更新的相关计算；具体的，系统状态变量X在第k个采样时刻的无迹卡尔曼滤波算法基本步骤如下：
S21，利用式(4)和(5)获得2n+1个采样点X(i)(k|k)即Sigma点集，以及每个采样点X(i)(k|k)对应的权值；上标(i)表示第几个采样点，i＝1，2，…，2n+1；
Figure BDA0001804933460000101
其中，k表示第k个采样时刻；
Figure BDA0001804933460000102
表示系统状态变量的均值；P(k|k)表示系统状态变量的方差；λ表示缩放比例系数；n表示系统状态变量X的维度，本实施例中，增广前的系统状态变量维度为6，增广后的系统状态变量维度为7；
S22，分别计算每个采样点X(i)(k|k)的一步预测值X(i)(k+1|k)，计算方式为：
X(i)(k+1|k)＝f[k,X(i)(k|k)]
其中，f(·)表示状态方程；
S23，根据步骤S22中得到每个采样点X(i)(k|k)的一步预测值X(i)(k+1|k)和每个采样点X(i)(k|k)对应的权值ω(i)，对该2n+1个采样点的一步预测值进行加权求和，得到系统状态变量的一步预测值
Figure BDA0001804933460000103
再根据系统状态变量的一步预测值
Figure BDA0001804933460000104
得到系统的协方差矩阵P(k+1|k)，具体计算方式为：
Figure BDA0001804933460000105
Figure BDA0001804933460000106
其中，权值ω(i)通过式(5)得到；
S24，根据系统状态变量的一步预测值
Figure BDA0001804933460000107
使用无迹变换UT，产生2n+1个新的采样点X(i)(k+1|k)，即新的Sigma点集，以及每个新的采样点对应的权值；上标(i)表示第几个采样点，i＝1，2，…，2n+1；2n+1个新的采样点X(i)(k+1|k)的计算方式为：
Figure BDA0001804933460000111
S25，分别将由步骤S24得到的2n+1个新的采样点X(i)(k+1|k)代入观测方程，分别得到此2n+1个新的采样点X(i)(k+1|k)的观测变量的预测量Y(i)(k+1|k)，计算方式为：
Y(i)(k+1|k)＝g[X(i)(k+1|k)]
其中，g(·)表示观测方程；
S26，根据步骤S25得到此2n+1个新的采样点X(i)(k+1|k)的观测变量的预测量Y(i)(k+1|k)，对此2n+1个新的采样点X(i)(k+1|k)的观测变量的预测量Y(i)(k+1|k)加权求和，得到观测变量的均值
Figure BDA00018049334600001111
根据观测变量的均值得到观测变量的协方差
Figure BDA0001804933460000112
系统状态变量和观测变量之间的协方差
Figure BDA0001804933460000113
具体计算方式为为：
Figure BDA0001804933460000114
Figure BDA0001804933460000115
Figure BDA0001804933460000116
S27，计算并更新系统的卡尔曼增益矩阵K(k+1)，计算方式为；
Figure BDA0001804933460000117
S28，计算并更新系统状态变量
Figure BDA0001804933460000118
和系统的协方差矩阵P(k+1|k+1)，计算方式为：
Figure BDA0001804933460000119
Figure BDA00018049334600001110
其中，Y(k+1)表示经无迹卡尔曼滤波算法得到的观测变量的值；KT(k+1)表示卡尔曼增益矩阵的转置；
本发明在无迹卡尔曼滤波算法的基础上，根据系统的待估计参数的变化速度远远小于系统的状态变量的变化速度的特点，增加了双时间尺度的概念，并设置一个较大的时间间隔为宏观尺度和一个较小的时间间隔为微观尺度。设置双时间尺度的目的在于，在宏观尺度下仅对系统的状态变量、协方差矩阵、卡尔曼增益矩阵以及观测量的部分行列进行更新，所述观测量的部分行列为系统的状态变量所占的维度；在微观尺度下对系统的状态变量、系统的待估计参数协方差矩阵、卡尔曼增益矩阵以及观测量的全部行列进行更新，所述观测量的全部行列为增广后的系统的状态变量和系统的待估计参数的总维度，以减少不必要的计算量，从而提升运算效率。
本实施例中，基于双时间尺度无迹卡尔曼滤波算法的部分代码如下所示：
Figure BDA0001804933460000121
Figure BDA0001804933460000131
如图2和图3所示，图2为利用传统的无迹卡尔曼滤波的故障参数估计方法进行的仿真实验；图3为利用本发明的基于双时间尺度无迹卡尔曼滤波的系统故障参数估计方法并分别以宏观步长为2、5、10、15进行的仿真实验；其中，真实值为仿真实验中注入的前轮传感器故障，β3为前轮处传感器的有效因子，真实值的曲线即为β3的变化曲线。
实验结果如下表所示：
UKF step＝2 step＝5 step＝10 step＝15
12.97s 11.60s 11.52s 11.48s 11.41s
其中，第一列表示利用传统的无迹卡尔曼滤波的故障参数估计方法进行故障参数估计时所需的时间，第二、三、四、五列分别表示利用本发明的基于双时间尺度无迹卡尔曼滤波的系统故障参数估计方法进行仿真实验，并分别以宏观步长为2、5、10、15进行故障参数估计时所需的时间。由此可知，本发明的基于双时间尺度无迹卡尔曼滤波的系统故障参数估计方法有效的提高了故障诊断的速率，且宏观步长取值越大，估计出故障参数所需的时间越短。
以上仅为本发明创造的较佳实施例而已，并不用以限制本发明创造，凡在本发明创造的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明创造的保护范围之内。

Claims (4)

1.基于双时间尺度无迹卡尔曼滤波的系统故障参数估计方法，其特征在于，包括以下步骤：
S1，对电动代步车系统的连续的系统状态空间方程进行离散化处理，得到时间离散的系统状态空间方程；
S2，将L2-L1个待估计参数添加到电动代步车系统的L1个原系统状态变量中，得到增广后的系统状态空间方程；
其中，所述待估计参数为故障集合中的参数；增广前的系统状态空间方程的维度为L1，增广前的系统状态变量为L1个；增广后的系统状态空间方程的维度为L2，增广后的系统状态变量为L2个；
S3，采用具有双时间尺度的无迹卡尔曼滤波算法对L2-L1个待估计参数与L1个原系统状态变量进行联合估计，得到参数估计值；
步骤S2中，所述增广后的系统状态空间方程，如公式(3)所示：
Figure FDA0003822702980000011
其中，x1(k)表示第k个采样时刻的待估计参数，即为电动代步车系统的新添加的系统状态变量；x1(·)表示待估计参数，即新添加的系统状态变量；
步骤S3中，在宏观尺度下，仅对增广前的L1个系统状态变量进行计算，并对系统的协方差矩阵的部分行列、卡尔曼增益矩阵的部分行列以及观测变量的部分行列进行更新，所述部分行列为增广前的L1个系统状态变量所占的维度L1；在微观尺度下，对增广后的L2个系统状态变量进行计算，并对系统的协方差矩阵的全部行列、卡尔曼增益矩阵的全部行列以及观测变量的全部行列进行更新，所述全部行列为增广后的L2个系统状态变量所占的总维度L2；所述协方差矩阵为系统状态变量在相邻的两个采样点所对应的采样时刻之间的协方差，用于描述系统状态变量在相邻的两个采样时刻之间的递推依赖关系；所述卡尔曼增益矩阵为系统在卡尔曼滤波递推过程中从某个采样时刻到下一采样时刻的增益矩阵。
2.根据权利要求1所述的基于双时间尺度无迹卡尔曼滤波的系统故障参数估计方法，其特征在于，步骤S3中，所述双时间尺度包括宏观尺度和微观尺度；在宏观尺度下，仅对增广前的L1个系统状态变量采用无迹卡尔曼滤波进行估计；在微观尺度下，对增广后的L2个系统状态变量采用无迹卡尔曼滤波进行估计。
3.根据权利要求2所述的基于双时间尺度无迹卡尔曼滤波的系统故障参数估计方法，其特征在于，步骤S1中，所述电动代步车系统的连续的系统状态空间方程包括状态方程和观测方程，如公式(1)所示：
Figure FDA0003822702980000021
其中，t表示实际的连续时间，x(t)表示t时刻下的系统状态变量，即增广前的原系统状态变量；u(t)表示t时刻下的系统输入；W表示过程噪声；
Figure FDA0003822702980000022
表示t时刻下的系统状态变量的一阶导数；y(t)为t时刻下的系统输出，即观测变量；V表示观测噪声；f(·)和g(·)分别表示状态函数和观测函数，即分别对应为状态方程和观测方程；x(·)表示原系统状态变量；y(·)表示观测变量；
采用前向差分法或后向差分法对所述连续的系统状态空间方程进行离散化，得到时间离散的系统状态空间方程，所述时间离散的系统状态空间方程，如公式(2)所示：
Figure FDA0003822702980000023
其中，k表示第k个采样点所对应的时刻，即表示第k个采样时刻；k+1表示第k+1个采样点所对应的采样时刻，即表示第k+1个采样时刻；x(k)表示第k个采样时刻的系统状态变量；u(k)表示第k个采样时刻的系统输入；x(k+1)表示第k+1个采样时刻的系统状态变量；y(k)为第k个采样采样时刻的系统输出。
4.根据权利要求1所述的基于双时间尺度无迹卡尔曼滤波的系统故障参数估计方法，其特征在于，利用无迹卡尔曼滤波算法进行计算并更新系统的协方差矩阵、卡尔曼增益矩阵、观测变量；所述无迹卡尔曼滤波算法，包括以下步骤：
S31，对系统状态变量进行无迹变换，由无迹变换获得一组采样点，以及每个采样点所对应的权值，并将该组采样点称为Sigma点集；
S32，分别计算Sigma点集中的每个采样点的一步预测值；
S33，根据每个采样点所对应的权值和步骤S32中得到每个采样点的一步预测值，进行加权求和的计算，得到系统状态变量的一步预测值，再根据系统状态变量的一步预测值得到系统的协方差矩阵；
S34，对系统状态变量的一步预测值进行无迹变换，产生一组新的采样点，即新的Sigma点集，以及每个新的采样点对应的权值；
S35，分别将由步骤S34得到的每个新的采样点代入观测方程，分别得到每个新的采样点的观测变量的预测量；
S36，根据步骤S35得到每个新的采样点的观测变量的预测量，进行加权求和计算，得到观测变量的均值；根据观测变量的均值得到观测变量的协方差和系统状态变量和观测变量之间的协方差；
S37，计算并更新系统的卡尔曼增益矩阵；
S38，计算并更新系统状态变量和系统的协方差矩阵。
CN201811093396.3A 2018-09-19 2018-09-19 基于双时间尺度无迹卡尔曼滤波的系统故障参数估计方法 Active CN109086247B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201811093396.3A CN109086247B (zh) 2018-09-19 2018-09-19 基于双时间尺度无迹卡尔曼滤波的系统故障参数估计方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201811093396.3A CN109086247B (zh) 2018-09-19 2018-09-19 基于双时间尺度无迹卡尔曼滤波的系统故障参数估计方法

Publications (2)

Publication Number Publication Date
CN109086247A CN109086247A (zh) 2018-12-25
CN109086247B true CN109086247B (zh) 2022-10-18

Family

ID=64842168

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201811093396.3A Active CN109086247B (zh) 2018-09-19 2018-09-19 基于双时间尺度无迹卡尔曼滤波的系统故障参数估计方法

Country Status (1)

Country Link
CN (1) CN109086247B (zh)

Families Citing this family (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109645995B (zh) * 2019-01-16 2021-09-07 杭州电子科技大学 基于肌电模型和无迹卡尔曼滤波的关节运动估计方法
CN110109049B (zh) * 2019-03-27 2021-04-20 北京邮电大学 用于大规模天线角度估计的无迹卡尔曼滤波方法及装置
CN110263924B (zh) * 2019-06-19 2021-08-17 北京计算机技术及应用研究所 一种神经元群模型的参数和状态估计方法
CN110161394A (zh) * 2019-07-04 2019-08-23 苏州妙益科技股份有限公司 一种基于无迹卡尔曼滤波的绝缘检测方法
CN110768600A (zh) * 2019-11-08 2020-02-07 江苏科技大学 一种pmsm无速度传感器转子检测方法

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101615794A (zh) * 2009-08-05 2009-12-30 河海大学 基于无迹变换卡尔曼滤波的电力系统动态状态估计方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8639476B2 (en) * 2008-04-22 2014-01-28 The United States Of America As Represented By The Secretary Of The Navy Process for estimation of ballistic missile boost state

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101615794A (zh) * 2009-08-05 2009-12-30 河海大学 基于无迹变换卡尔曼滤波的电力系统动态状态估计方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
基于改进强跟踪滤波的广义系统传感器故障诊断及隔离;梁天添等;《中国惯性技术学报》;20180815(第04期);全文 *
基于自适应无迹卡尔曼滤波的分布式驱动电动汽车车辆状态参数估计;王震坡等;《北京理工大学学报》;20180715(第07期);全文 *

Also Published As

Publication number Publication date
CN109086247A (zh) 2018-12-25

Similar Documents

Publication Publication Date Title
CN109086247B (zh) 基于双时间尺度无迹卡尔曼滤波的系统故障参数估计方法
CN108304612B (zh) 基于噪声补偿的迭代平方根ckf的汽车雷达目标跟踪方法
CN109472418B (zh) 基于卡尔曼滤波的机动目标状态预测优化方法
CN103927436A (zh) 一种自适应高阶容积卡尔曼滤波方法
CN108363824B (zh) 一种基于pso优化ukf的ehb系统压力估计方法
CN111090945B (zh) 一种针对切换系统的执行器和传感器故障估计设计方法
Lang A nonparametric polynomial identification algorithm for the Hammerstein system
CN105277896A (zh) 基于elm-mukf的锂电池剩余寿命预测方法
CN104462015B (zh) 处理非高斯Lévy噪声的分数阶线性离散系统状态更新方法
CN106896324A (zh) 一种soc估计方法
Wang et al. Comparison of Kalman Filter-based state of charge estimation strategies for Li-Ion batteries
CN106154168A (zh) 数据驱动的动力电池荷电状态估计方法
CN111942399A (zh) 一种基于无迹卡尔曼滤波的车速估算方法及系统
CN104794101A (zh) 一种分数阶非线性系统状态估计方法
CN116861155A (zh) 一种非线性模型下锂电池的遗忘因子梯度下降算法
CN111709147A (zh) 一种基于水文失事机理的设计洪水地区组成计算方法
CN105677936B (zh) 机电复合传动系统需求转矩的自适应递归多步预测方法
WO2021210526A1 (ja) 残容量推定装置、モデル生成装置、残容量推定方法、モデル生成方法、及びプログラム
CN104820788A (zh) 计及Lévy噪声的分数阶扩展卡尔曼滤波方法
CN113608442A (zh) 基于特征函数的非线性状态模型系统的状态估计方法
CN110470925B (zh) 一种基于可拓关联函数的电驱动力总成可靠度测试方法
CN113468471A (zh) 基于张量型加权Schatten-p范数的交通数据修复方法
CN111950123A (zh) 一种陀螺仪误差系数曲线拟合预测方法及系统
Öztürk et al. Mathematical estimation of expenditures in the health sector in turkey with grey modeling
CN114815785B (zh) 一种基于有限时间观测器的非线性系统执行器鲁棒故障估计方法

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