CN104794101A - 一种分数阶非线性系统状态估计方法 - Google Patents

Abstract

本发明公开了一种新的分数阶非线性系统状态估计方法。首先，给出状态预测值和预测误差协方差初始值。接着，利用近似处理方法求解得到非高斯Lévy噪声下状态近似值和量测近似值，并由此计算系统噪声协方差和量测噪声协方差。然后，利用量测噪声协方差计算最优滤波增益，结合已得到的最优滤波增益计算状态估计值和估计误差协方差。最后，利用当前时刻状态估计值更新下一时刻状态预测值，利用当前时刻估计误差协方差和系统噪声协方差更新下一时刻状态预测协方差。本发明由于解决了带有非高斯Lévy噪声的分数阶非线性系统状态估计问题，从而拓展了分数阶理论的应用范围，而且本发明易于与已有的状态估计软件相结合。

Description

一种分数阶非线性系统状态估计方法

技术领域

本发明涉及一种新的分数阶非线性系统状态估计方法，属于系统分析与处理技术领域。

背景技术

系统分析与处理，旨在研究特定系统结构中各部分(各子系统)的相互作用，系统的对外接口与界面，以及系统整体的行为、功能和局限，从而为系统未来的变迁与有关决策提供参考和依据，其目标之一在于改善决策过程及系统性能，以达到系统的整体最优。在系统分析与处理领域，状态估计起着至关重要的作用。状态估计是根据可获取的量测数据估算动态系统内部状态的方法，由于对系统的输入和输出进行量测而得到的数据只能反映系统的外部特性，而系统的动态规律需要用内部状态变量来描述，因此状态估计对于了解和控制一个系统具有重要意义。

卡尔曼滤波作为一种状态估计的方法被广泛应用于各类系统中。传统的卡尔曼滤波只限于整数阶系统，直到分数阶微积分有了很大的发展，分数阶卡尔曼滤波才慢慢地被学者们研究并逐步得到广泛应用。对非线性系统来说，应用最多的是扩展卡尔曼滤波，不管是整数阶系统还是分数阶系统，扩展卡尔曼滤波都得到广泛应用，但是传统的扩展卡尔曼滤波要求系统噪声和量测噪声均为高斯白噪声，这些过于理想化的要求在现实在很难得到满足。因此，研究非高斯白噪声下的扩展卡尔曼滤波问题具有很重要的理论与现实意义。

发明内容

发明目的：基于以上分析，本发明采用分数阶系统理论，提出一种新的分数阶非线性系统状态估计方法，以期提高状态估计的估计精度。

由于非高斯噪声普遍存在于实际非线性系统中，那么利用传统的扩展卡尔曼滤波方法对这些系统进行状态估计时不能得出理想结果。本发明中所提到的一种新的非线性系统状态估计方法通过近似处理得到状态量和量测量的近似值，并由此推导得出系统噪声协方差和量测噪声协方差，再根据量测噪声协方差和已给定的状态预测误差协方差计算最优滤波增益，从而得到状态估计值和估计误差协方差。该方法与常规动态状态估计程序相结合，能很好地处理非高斯噪声情况下非线性系统的状态估计问题。

技术方案：一种分数阶非线性系统状态估计方法，所述方法是在计算机中依次按以下步骤实现的：

(1)、初始化，包括：设定状态预测量的初值和预测误差协方差的初值；

(2)、对非高斯Lévy噪声进行近似处理，推导得出状态近似值和量测输出近似值，计算步骤为：

${\stackrel{\&OverBar;}{x}}_{k+1}=\left\{\begin{array}{cc}{\mathrm{\&Omega;}}_1+{\mathrm{\&delta;}}_1\&CenterDot;\mathrm{sign}\left({\mathrm{\&Omega;}}_2\right)& \mathrm{if}\left|{\mathrm{\&Omega;}}_2\right|\&GreaterEqual;{\mathrm{\&delta;}}_1\\ {\stackrel{\&OverBar;}{x}}_{k+1}& \mathrm{if}\left|{\mathrm{\&Omega;}}_2\right|<{\mathrm{\&delta;}}_1\end{array}\right.$

${\stackrel{\&OverBar;}{y}}_k=\left\{\begin{array}{cc}h\left(\widetilde{\stackrel{\&OverBar;}{x}}\right)+{\mathrm{\&delta;}}_2\&CenterDot;\mathrm{sign}({\stackrel{\&OverBar;}{y}}_k-h\left(\widetilde{\stackrel{\&OverBar;}{x}}\right))& \mathrm{if}|{\stackrel{\&OverBar;}{y}}_k-h\left(\widetilde{\stackrel{\&OverBar;}{x}}\right)|\&GreaterEqual;{\mathrm{\&delta;}}_2\\ {\stackrel{\&OverBar;}{y}}_k& \mathrm{if}|{\stackrel{\&OverBar;}{y}}_k-h\left(\widetilde{\stackrel{\&OverBar;}{x}}\right)|<{\mathrm{\&delta;}}_2\end{array}\right.$

式中

${\mathrm{\&Omega;}}_1=f({\widetilde{\stackrel{\&OverBar;}{x}}}_k,u_k)-\underset{j=1}{\overset{k+1}{\mathrm{\&Sigma;}}}{(-1)}^j{\mathrm{\&gamma;}}_j{\widetilde{\stackrel{\&OverBar;}{x}}}_{k+1-j}$

${\mathrm{\&Omega;}}_2={\stackrel{\&OverBar;}{x}}_{k+1}-f({\widetilde{\stackrel{\&OverBar;}{x}}}_k,u_k)+\underset{j=1}{\overset{k+1}{\mathrm{\&Sigma;}}}{(-1)}^j{\mathrm{\&gamma;}}_j{\widetilde{\stackrel{\&OverBar;}{x}}}_{k+1-j}$

表示含有Lévy噪声的状态向量，下标k+1表示第k+1时刻，δ₁表示所取的第一个阈值，sign(x)表示符号函数，如果x>0，返回1；如果x<0，返回-1。|x|表示取x的绝对值，表示量测输出向量，表示非线性函数，表示第k时刻状态预测向量，δ₂表示所取的第二个阈值，在Ω₁中，表示非线性函数，u_k表示第k时刻控制输入向量，γ_j表示

${\mathrm{\&gamma;}}_j=\mathrm{diag}\left(\right[\left(\begin{array}{c}{\mathrm{\&alpha;}}_1\\ j\end{array}\right)\left(\begin{array}{c}{\mathrm{\&alpha;}}_2\\ j\end{array}\right)\&CenterDot;\&CenterDot;\&CenterDot;\left(\begin{array}{c}{\mathrm{\&alpha;}}_N\\ j\end{array}\right)\left]\right),$ 其中 $\left(\begin{array}{c}\mathrm{\&alpha;}\\ j\end{array}\right)=\left\{\begin{array}{cc}1& j=0\\ \frac{\mathrm{\&alpha;}(\mathrm{\&alpha;}-1)\&CenterDot;\&CenterDot;\&CenterDot;(\mathrm{\&alpha;}-j+1)}{j!}& j>0\end{array}\right.$

α_N表示第N个分数阶阶次值，j＝1,2,…,k+1。

(3)、利用已得到的状态近似值和量测输出近似值，计算当前时刻系统噪声协方差和量测噪声协方差，计算步骤为：

$\begin{array}{c}{\stackrel{\&OverBar;}{Q}}_k=({\stackrel{\&OverBar;}{x}}_{k+1}-f({\widehat{\stackrel{\&OverBar;}{x}}}_k,u_k)+\underset{j=1}{\overset{k+1}{\mathrm{\&Sigma;}}}{(-1)}^j{\mathrm{\&gamma;}}_j{\widehat{\stackrel{\&OverBar;}{x}}}_{k+1-j}){({\stackrel{\&OverBar;}{x}}_{k+1}-f({\widehat{\stackrel{\&OverBar;}{x}}}_k,u_k)+\underset{j=1}{\overset{k+1}{\mathrm{\&Sigma;}}}{(-1)}^j{\mathrm{\&gamma;}}_j{\widehat{\stackrel{\&OverBar;}{x}}}_{k+1-j})}^T\\ +(F_k+{\mathrm{\&gamma;}}_1){\widehat{\stackrel{\&OverBar;}{p}}}_k{(F_k+{\mathrm{\&gamma;}}_2)}^T+\underset{j=2}{\overset{k+1}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_2{\widehat{\stackrel{\&OverBar;}{p}}}_{k+1-j}{\mathrm{\&gamma;}}_j^T\end{array}$

${\stackrel{\&OverBar;}{R}}_k=[{\stackrel{\&OverBar;}{y}}_k-h\left({\widetilde{\stackrel{\&OverBar;}{x}}}_k\right)]{[{\stackrel{\&OverBar;}{y}}_k-h\left({\widetilde{\stackrel{\&OverBar;}{x}}}_k\right)]}^T+H_k{\widetilde{\stackrel{\&OverBar;}{p}}}_k{H}_k^T$

式中，表示非线性函数，表示第k时刻状态估计向量，上标T表示转置，表示非线性函数在处的雅克比矩阵，表示第k时刻状态估计误差协方差，表示非线性函数处的雅克比矩阵，表示第k时刻状态预测误差协方差。

(4)、利用当前时刻量测噪声协方差和给定的当前时刻预测误差协方差计算当前时刻最优滤波增益，计算步骤为：

${\stackrel{\&OverBar;}{K}}_k={\widetilde{\stackrel{\&OverBar;}{p}}}_k{h}_k^T{(H_k{\widetilde{\stackrel{\&OverBar;}{p}}}_k{H}_k^T+{\stackrel{\&OverBar;}{R}}_k)}^{-1}$

(5)、结合已得到的当前时刻最优滤波增益和给定的当前时刻状态预测值计算当前时刻状态估计值，计算步骤为：

${\stackrel{\&OverBar;}{x}}_k={\widetilde{\stackrel{\&OverBar;}{x}}}_k+{\stackrel{\&OverBar;}{K}}_k[{\stackrel{\&OverBar;}{y}}_k-h\left({\widetilde{\stackrel{\&OverBar;}{x}}}_k\right)]$

(6)、利用当前时刻最优滤波增益和给定的当前时刻预测误差协方差计算当前时刻估计误差协方差，计算步骤为：

${\widehat{\stackrel{\&OverBar;}{p}}}_k=(I-{\stackrel{\&OverBar;}{K}}_k{H}_k){\widetilde{\stackrel{\&OverBar;}{p}}}_k$

(7)、利用当前时刻状态估计值更新下一时刻状态预测值，计算步骤为：

$\left\{\begin{array}{c}{\mathrm{\&Delta;}}^{\mathrm{\&alpha;}}{\widetilde{\stackrel{\&OverBar;}{x}}}_{k+1}=f({\widehat{\stackrel{\&OverBar;}{x}}}_k,u_k)\\ {\widetilde{\stackrel{\&OverBar;}{x}}}_{k+1}={\mathrm{\&Delta;}}^{\mathrm{\&alpha;}}{\widetilde{\stackrel{\&OverBar;}{x}}}_{k+1}-\underset{j=1}{\overset{k+1}{\mathrm{\&Sigma;}}}{(-1)}^j{\mathrm{\&gamma;}}_j{\widehat{\stackrel{\&OverBar;}{x}}}_{k+1-j}\end{array}\right.$

式中Δ表示分数阶算子，上标α表示分数阶阶次向量。

(8)、利用当前时刻估计误差协方差和系统噪声协方差更新下一时刻状态预测协方差，计算步骤为：

${\widetilde{\stackrel{\&OverBar;}{p}}}_{k+1}=(F_k+{\mathrm{\&gamma;}}_1){\widehat{\stackrel{\&OverBar;}{p}}}_k{(F_k+{\mathrm{\&gamma;}}_1)}^T+{\stackrel{\&OverBar;}{Q}}_k+\underset{j=2}{\overset{k+1}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_j{\widehat{\stackrel{\&OverBar;}{p}}}_{k+1-j}{\mathrm{\&gamma;}}_j^T$

(9)、判断k+1是否大于等于步长L，如果是，结束计算，否则返回步骤(2)进行下一次估计。

传统的扩展卡尔曼滤波算法处理的是整数阶系统，并且要求系统噪声和量测噪声均为高斯白噪声，这意味着系统的阶次必须是整数，噪声的概率分布是正态分布，而且满足它的二阶矩不相关。这些理想化的定义在实际系统中很难得到满足，比如一些粘弹性结构，有耗网络和扩散波等，这些本身就有着分数阶的性质，用整数阶模型很难对这些系统进行建模。另一方面，实际系统中的噪声大多为非高斯白噪声，这样的系统处理起来具有一定的复杂性和挑战性。

本发明提出一种新的分数阶非线性系统状态估计方法，在传统卡尔曼滤波的基础上运用分数阶系统理论，通过非高斯Lévy噪声替代传统的高斯白噪声，利用近似代替方法分解得到Lévy噪声的分解值，重新计算系统噪声和量测噪声的协方差，进而得到最优滤波增益，最后计算状态估计值和估计误差协方差。本发明由于结合分数阶系统理论和非高斯Lévy噪声，能有效地解决整数阶系统的局限性问题和非高斯噪声下的状态估计问题。

附图说明

图1为本发明实施例的方法流程图；

图2为本发明提出的状态估计方法得到的真实值与估计值；

图3为真实值与估计值之间的误差。

具体实施方式

下面结合具体实施例，进一步阐明本发明，应理解这些实施例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。

如图1所示，分数阶非线性系统状态估计方法，包括如下步骤：

(1)、初始化。包括：设定状态预测量的初值和预测误差协方差的初值。

(2)、对非高斯Lévy噪声进行近似处理，推导得出状态近似值和量测输出近似值。

(3)利用已得到的状态近似值和量测输出近似值，计算当前时刻系统噪声协方差和量测噪声协方差。

(4)、利用当前时刻量测噪声协方差和给定的当前时刻预测误差协方差计算当前时刻最优滤波增益。

(5)、结合已得到的当前时刻最优滤波增益和给定的当前时刻状态预测值计算当前时刻状态估计值。

(6)、利用当前时刻最优滤波增益和给定的当前时刻预测误差协方差计算当前时刻估计误差协方差。

(7)、利用当前时刻状态估计值更新下一时刻状态预测值。

(8)、利用当前时刻估计误差协方差和系统噪声协方差更新下一时刻状态预测协方差。

(9)、判断k+1是否大于等于步长L，如果是，结束计算，否则返回步骤(2)进行下一次估计。

目前分数阶非线性系统状态估计方法主要采用的是20世纪70年代初由R.E.Larson和Debs提出的卡尔曼滤波算法扩展而来的分数阶卡扩展尔曼滤波算法。

卡尔曼滤波是一种递推回归方法，它的基本思想是根据新量测的数据和由前面一步量测数据计算而得的估计值，再来推算出新的估计值，即

新估计值＝旧估计值十修正值

分数阶扩展卡尔曼滤波的本质是求出k时刻状态量真值x_k的最优估计值估计的准则是以状态量的估计误差协方差阵最小为目标函数，即：

$\min {\hat{P}}_k=\min E\left[{\hat{e}}_k{\hat{e}}_k^T\right]$

式中E为求期望函数，x_k为k时刻状态量的真值，为k时刻状态量的估计值，上标T表示转置。

由于实际中的一些非线性系统具有分数阶性质，并且存在非高斯白噪声，因此传统扩展卡尔曼滤波算法在应用上受到了很大的限制。分数阶扩展卡尔曼滤波算法具有全局相关性，能较好地体现系统函数发展的历史依赖过程，因而能处理具有分数阶性质的非线性系统，非高斯Lévy噪声通过近似分解可以得到具有高斯白噪声的性质。对于具有非高斯Lévy噪声的分数阶非线性系统，模型描述可以表示为：

$\left\{\begin{array}{c}{\mathrm{\&Delta;}}^{\mathrm{\&alpha;}}{\stackrel{\&OverBar;}{x}}_{k+1}=f({\stackrel{\&OverBar;}{x}}_k,u_k)+{\stackrel{\&OverBar;}{w}}_k\\ {\stackrel{\&OverBar;}{x}}_{k+1}={\mathrm{\&Delta;}}^{\mathrm{\&alpha;}}{\stackrel{\&OverBar;}{x}}_{k+1}-\underset{j=1}{\overset{k+1}{\mathrm{\&Sigma;}}}{(-1)}^j{\mathrm{\&gamma;}}_j{\stackrel{\&OverBar;}{x}}_{k+1-j}\\ {\stackrel{\&OverBar;}{y}}_k=h\left({\stackrel{\&OverBar;}{x}}_k\right)+{\stackrel{\&OverBar;}{v}}_k\end{array}\right.$

式中Δ表示分数阶算子，上标α表示分数阶阶次向量，和分别表示非线性函数，表示含有Lévy噪声的状态向量，下标k表示为第k时刻，u_k表示控制输入向量，为量测输出向量，和均为非高斯Lévy噪声，两者协方差矩阵分别表示为和γ_j表示

${\mathrm{\&gamma;}}_j=\mathrm{diag}\left(\right[\left(\begin{array}{c}{\mathrm{\&alpha;}}_1\\ j\end{array}\right)\left(\begin{array}{c}{\mathrm{\&alpha;}}_2\\ j\end{array}\right)\&CenterDot;\&CenterDot;\&CenterDot;\left(\begin{array}{c}{\mathrm{\&alpha;}}_N\\ j\end{array}\right)\left]\right),$ 其中 $\left(\begin{array}{c}\mathrm{\&alpha;}\\ j\end{array}\right)=\left\{\begin{array}{cc}1& j=0\\ \frac{\mathrm{\&alpha;}(\mathrm{\&alpha;}-1)\&CenterDot;\&CenterDot;\&CenterDot;(\mathrm{\&alpha;}-j+1)}{j!}& j>0\end{array}\right.$

式中α_N表示第N个分数阶阶次值，j＝1,2,…,k+1。

在以上所示分数阶非线性系统模型中，系统噪声和量测噪声均为非高斯Lévy噪声，由于其具有无限方差特性，因此不能用来直接计算协方差矩阵，必须用近似处理的方法得出具有高斯白噪声性质的近似值才可以用来计算协方差矩阵。近似处理后状态向量和量测输出向量可以表示为：

${\stackrel{\&OverBar;}{x}}_{k+1}=\left\{\begin{array}{cc}{\mathrm{\&Omega;}}_1+{\mathrm{\&delta;}}_1\&CenterDot;\mathrm{sign}\left({\mathrm{\&Omega;}}_2\right)& \mathrm{if}\left|{\mathrm{\&Omega;}}_2\right|\&GreaterEqual;{\mathrm{\&delta;}}_1\\ {\stackrel{\&OverBar;}{x}}_{k+1}& \mathrm{if}\left|{\mathrm{\&Omega;}}_2\right|<{\mathrm{\&delta;}}_1\end{array}\right.$

${\stackrel{\&OverBar;}{y}}_k=\left\{\begin{array}{cc}h\left(\widetilde{\stackrel{\&OverBar;}{x}}\right)+{\mathrm{\&delta;}}_2\&CenterDot;\mathrm{sign}({\stackrel{\&OverBar;}{y}}_k-h\left(\widetilde{\stackrel{\&OverBar;}{x}}\right))& \mathrm{if}|{\stackrel{\&OverBar;}{y}}_k-h\left(\widetilde{\stackrel{\&OverBar;}{x}}\right)|\&GreaterEqual;{\mathrm{\&delta;}}_2\\ {\stackrel{\&OverBar;}{y}}_k& \mathrm{if}|{\stackrel{\&OverBar;}{y}}_k-h\left(\widetilde{\stackrel{\&OverBar;}{x}}\right)|<{\mathrm{\&delta;}}_2\end{array}\right.$

式中

${\mathrm{\&Omega;}}_1=f({\widetilde{\stackrel{\&OverBar;}{x}}}_k,u_k)-\underset{j=1}{\overset{k+1}{\mathrm{\&Sigma;}}}{(-1)}^j{\mathrm{\&gamma;}}_j{\widetilde{\stackrel{\&OverBar;}{x}}}_{k+1-j}$

${\mathrm{\&Omega;}}_2={\stackrel{\&OverBar;}{x}}_{k+1}-f({\widetilde{\stackrel{\&OverBar;}{x}}}_k,u_k)+\underset{j=1}{\overset{k+1}{\mathrm{\&Sigma;}}}{(-1)}^j{\mathrm{\&gamma;}}_j{\widetilde{\stackrel{\&OverBar;}{x}}}_{k+1-j}$

表示含有Lévy噪声的状态向量，下标k+1表示第k+1时刻，δ₁表示所取的第一个阈值，sign(x)表示符号函数，如果x>0，返回1；如果x<0，返回-1。|x|表示取x的绝对值，表示量测输出向量，表示非线性函数，表示第k时刻状态预测向量，δ₂表示所取的第二个阈值，在Ω₁中，表示非线性函数，u_k表示第k时刻控制输入向量。

第k时刻系统噪声协方差和量测噪声协方差可分别表示为：

${\stackrel{\&OverBar;}{Q}}_k=E\left[{\stackrel{\&OverBar;}{w}}_k{\stackrel{\&OverBar;}{w}}_k^T\right]$

${\stackrel{\&OverBar;}{R}}_k=E\left[{\stackrel{\&OverBar;}{v}}_k{\stackrel{\&OverBar;}{v}}_k^T\right]$

利用给定模型中的表达式进行代换，可得到：

$\left\{\begin{array}{c}{\stackrel{\&OverBar;}{w}}_k={\stackrel{\&OverBar;}{x}}_{k+1}-f({\stackrel{\&OverBar;}{x}}_k,u_k)+\underset{j=1}{\overset{k+1}{\mathrm{\&Sigma;}}}{(-1)}^j{\mathrm{\&gamma;}}_j{\stackrel{\&OverBar;}{x}}_{k+1-j}\\ {\stackrel{\&OverBar;}{v}}_k={\stackrel{\&OverBar;}{y}}_k-h\left({\stackrel{\&OverBar;}{x}}_k\right)\end{array}\right.$

将上述两式分别代入系统和量测噪声协方差阵可得：

$\begin{array}{c}{\stackrel{\&OverBar;}{Q}}_k=({\stackrel{\&OverBar;}{x}}_{k+1}-f({\widehat{\stackrel{\&OverBar;}{x}}}_k,u_k)+\underset{j=1}{\overset{k+1}{\mathrm{\&Sigma;}}}{(-1)}^j{\mathrm{\&gamma;}}_j{\widehat{\stackrel{\&OverBar;}{x}}}_{k+1-j}){({\stackrel{\&OverBar;}{x}}_{k+1}-f({\widehat{\stackrel{\&OverBar;}{x}}}_k,u_k)+\underset{j=1}{\overset{k+1}{\mathrm{\&Sigma;}}}{(-1)}^j{\mathrm{\&gamma;}}_j{\widehat{\stackrel{\&OverBar;}{x}}}_{k+1-j})}^T\\ +(F_k+{\mathrm{\&gamma;}}_1){\widehat{\stackrel{\&OverBar;}{p}}}_k{(F_k+{\mathrm{\&gamma;}}_2)}^T+\underset{j=2}{\overset{k+1}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_2{\widehat{\stackrel{\&OverBar;}{p}}}_{k+1-j}{\mathrm{\&gamma;}}_j^T\end{array}$

${\stackrel{\&OverBar;}{R}}_k=[{\stackrel{\&OverBar;}{y}}_k-h\left({\widetilde{\stackrel{\&OverBar;}{x}}}_k\right)]{[{\stackrel{\&OverBar;}{y}}_k-h\left({\widetilde{\stackrel{\&OverBar;}{x}}}_k\right)]}^T+H_k{\widetilde{\stackrel{\&OverBar;}{p}}}_k{H}_k^T$

式中，表示非线性函数，表示第k时刻状态估计向量，上标T表示转置，表示非线性函数在处的雅克比矩阵，表示第k时刻状态估计误差协方差，表示非线性函数在处的雅克比矩阵，表示第k时刻状态预测误差协方差。

由已得到的k时刻的量测噪声协方差，可得k时刻最优滤波增益为：

${\stackrel{\&OverBar;}{K}}_k={\widetilde{\stackrel{\&OverBar;}{p}}}_k{h}_k^T{(H_k{\widetilde{\stackrel{\&OverBar;}{p}}}_k{H}_k^T+{\stackrel{\&OverBar;}{R}}_k)}^{-1}$

系统在k时刻的估计状态量和估计误差协方差可分别表示为：

${\widehat{\stackrel{\&OverBar;}{x}}}_k=E\left[{\stackrel{\&OverBar;}{x}}_k\right|{\stackrel{\&OverBar;}{y}}_k^{*}]$

${\widehat{\stackrel{\&OverBar;}{P}}}_k=\mathrm{cov}\left[{\stackrel{\&OverBar;}{x}}_k\right|{\stackrel{\&OverBar;}{y}}_k^{*}]$

其中，表示k-1时刻之前所有控制量和量测量的累计值，cov表示取协方差函数。

由k时候变到k+1时刻时，k+1时刻状态预测量和预测误差协方差分别为：

${\widetilde{\stackrel{\&OverBar;}{x}}}_{k+1}=E\left[{\stackrel{\&OverBar;}{x}}_{k+1}\right|{\stackrel{\&OverBar;}{y}}_k^{*}]$

${\widetilde{\stackrel{\&OverBar;}{P}}}_{k+1}=\mathrm{cov}\left[{\stackrel{\&OverBar;}{x}}_{k+1}\right|{\stackrel{\&OverBar;}{y}}_k^{*}]$

状态量估计步骤总结如下：

步骤1：在当前控制量和量测累计值基础上，计算当前时刻的状态估计量；

步骤2：在当前控制量和量测累计值基础上，计算当前时刻估计误差协方差；

步骤3：在当前控制量和量测累计值基础上，利用下一时刻的状态估计量更新状态预测量；

步骤4：在当前控制量和量测累计值基础上，利用下一时刻的估计误差协方差更新预测误差协方差。

步骤5：在步长范围内，利用下一时刻状态预测量更新状态估计量，利用预测误差协方差更新估计误差协方差。

根据以上步骤，k时刻状态量的估计值为：

$\begin{array}{c}{\widehat{\stackrel{\&OverBar;}{x}}}_k=E\left[{\stackrel{\&OverBar;}{x}}_k\right|{\stackrel{\&OverBar;}{y}}_k^{*}]\\ =E[f({\stackrel{\&OverBar;}{x}}_{k-1},u_{k-1})+{\stackrel{\&OverBar;}{w}}_{k-1}-\underset{j=1}{\overset{k}{\mathrm{\&Sigma;}}}{(-1)}^j{\mathrm{\&gamma;}}_j{\stackrel{\&OverBar;}{x}}_{k-1}|{\stackrel{\&OverBar;}{y}}_k^{*}]\\ ={\widetilde{\stackrel{\&OverBar;}{x}}}_k+{\stackrel{\&OverBar;}{K}}_k[{\stackrel{\&OverBar;}{y}}_k-h\left({\widetilde{\stackrel{\&OverBar;}{x}}}_k\right)]\end{array}$

式中，表示第k时刻状态的预测量，表示第k时刻最优滤波增益。

k时刻的状态估计误差协方差为：

$\begin{array}{c}{\widehat{\stackrel{\&OverBar;}{P}}}_k=\mathrm{cov}\left[{\stackrel{\&OverBar;}{x}}_k\right|{\stackrel{\&OverBar;}{y}}_k^{*}]\\ =\mathrm{cov}[f({\stackrel{\&OverBar;}{x}}_{k-1},u_{k-1})+{\stackrel{\&OverBar;}{w}}_{k-1}-\underset{j=1}{\overset{k}{\mathrm{\&Sigma;}}}{(-1)}^j{\mathrm{\&gamma;}}_j{\stackrel{\&OverBar;}{x}}_{k-j}|{\stackrel{\&OverBar;}{y}}_k^{*}]\\ =(I-{\stackrel{\&OverBar;}{K}}_k{H}_k){\widetilde{\stackrel{\&OverBar;}{p}}}_k\end{array}$

式中H_k表示非线性函数在处的雅克比矩阵，表示k时刻状态预测误差协方差。

接下来利用状态估计量更新状态预测量，计算步骤为：

$\begin{array}{c}{\widetilde{\stackrel{\&OverBar;}{x}}}_{k+1}=E\left[{\stackrel{\&OverBar;}{x}}_{k+1}\right|{\stackrel{\&OverBar;}{y}}_k^{*}]\\ =E[f({\stackrel{\&OverBar;}{x}}_k,u_k)+{\stackrel{\&OverBar;}{w}}_k-\underset{j=1}{\overset{k}{\mathrm{\&Sigma;}}}{(-1)}^j{\mathrm{\&gamma;}}_j{\stackrel{\&OverBar;}{x}}_{k+1-j}|{\stackrel{\&OverBar;}{y}}_k^{*}]\\ ={f({\widehat{\stackrel{\&OverBar;}{x}}}_k,u_k)-\underset{j=1}{\overset{k+1}{\mathrm{\&Sigma;}}}{(-1)}^j{\mathrm{\&gamma;}}_j\widehat{\stackrel{\&OverBar;}{x}}}_{k+1-j}\end{array}$

其中，表示第k时刻状态估计量，u_k表示第k时刻控制输入向量。

最后，利用估计误差协方差更新预测误差协方差，计算步骤为：

$\begin{array}{c}{\widetilde{\stackrel{\&OverBar;}{P}}}_{k+1}=\mathrm{cov}\left[{\stackrel{\&OverBar;}{x}}_{k+1}\right|{\stackrel{\&OverBar;}{y}}_k^{*}]\\ =\mathrm{cov}[f({\stackrel{\&OverBar;}{x}}_k,u_k)+{\stackrel{\&OverBar;}{w}}_k-\underset{j=1}{\overset{k+1}{\mathrm{\&Sigma;}}}{(-1)}^j{\mathrm{\&gamma;}}_j{\stackrel{\&OverBar;}{x}}_{k+1-j}|{\stackrel{\&OverBar;}{y}}_k^{*}]\\ ={(F_k+{\mathrm{\&gamma;}}_1)\widehat{\stackrel{\&OverBar;}{p}}}_k{(F_k+{\mathrm{\&gamma;}}_1)}^T+{\stackrel{\&OverBar;}{Q}}_k+\underset{j=2}{\overset{k+1}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_j{\widehat{\stackrel{\&OverBar;}{p}}}_{k+1-j}{\mathrm{\&gamma;}}_j^T\end{array}$

其中

$\begin{array}{c}{\stackrel{\&OverBar;}{Q}}_k=({\stackrel{\&OverBar;}{x}}_{k+1}-f({\widehat{\stackrel{\&OverBar;}{x}}}_k,u_k)+\underset{j=1}{\overset{k+1}{\mathrm{\&Sigma;}}}{(-1)}^j{\mathrm{\&gamma;}}_j{\widehat{\stackrel{\&OverBar;}{x}}}_{k+1-j}){({\stackrel{\&OverBar;}{x}}_{k+1}-f({\widehat{\stackrel{\&OverBar;}{x}}}_k,u_k)+\underset{j=1}{\overset{k+1}{\mathrm{\&Sigma;}}}{(-1)}^j{\mathrm{\&gamma;}}_j{\widehat{\stackrel{\&OverBar;}{x}}}_{k+1-j})}^T\\ +(F_k+{\mathrm{\&gamma;}}_1){\widehat{\stackrel{\&OverBar;}{p}}}_k{(F_k+{\mathrm{\&gamma;}}_2)}^T+\underset{j=2}{\overset{k+1}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_2{\widehat{\stackrel{\&OverBar;}{p}}}_{k+1-j}{\mathrm{\&gamma;}}_j^T\end{array}$

F_k表示非线性函数在处的雅克比矩阵，表示第k时刻估计误差协方差。

新的分数阶非线性系统状态估计方法计算公式如下：

a.预测步：

$\left\{\begin{array}{c}{\mathrm{\&Delta;}}^{\mathrm{\&alpha;}}{\widetilde{\stackrel{\&OverBar;}{x}}}_{k+1}=f({\widehat{\stackrel{\&OverBar;}{x}}}_k,u_k)\\ {\widetilde{\stackrel{\&OverBar;}{x}}}_{k+1}={\mathrm{\&Delta;}}^{\mathrm{\&alpha;}}{\widetilde{\stackrel{\&OverBar;}{x}}}_{k+1}-\underset{j=1}{\overset{k+1}{\mathrm{\&Sigma;}}}{(-1)}^j{\mathrm{\&gamma;}}_j{\widehat{\stackrel{\&OverBar;}{x}}}_{k+1-j}\\ {\widetilde{\stackrel{\&OverBar;}{p}}}_{k+1}=(F_k+{\mathrm{\&gamma;}}_1){\widehat{\stackrel{\&OverBar;}{p}}}_k{(F_k+{\mathrm{\&gamma;}}_1)}^T+{\stackrel{\&OverBar;}{Q}}_k+\underset{j=2}{\overset{k+1}{\mathrm{\&Sigma;}}}{\widehat{\stackrel{\&OverBar;}{p}}}_{k+1-j}{\mathrm{\&gamma;}}_j^T\end{array}\right.$

b.估计步：

$\left\{\begin{array}{c}{\stackrel{\&OverBar;}{K}}_k={\widetilde{\stackrel{\&OverBar;}{p}}}_k{H}_k^T{(H_k{\widetilde{\stackrel{\&OverBar;}{p}}}_k{H}_k^T+{\stackrel{\&OverBar;}{R}}_k)}^{-1}\\ {\widehat{\stackrel{\&OverBar;}{x}}}_k={\widetilde{\stackrel{\&OverBar;}{x}}}_k+{\stackrel{\&OverBar;}{K}}_k[{\stackrel{\&OverBar;}{y}}_k-h\left({\widetilde{\stackrel{\&OverBar;}{x}}}_k\right)]\\ {\widehat{\stackrel{\&OverBar;}{p}}}_k=(I-{\stackrel{\&OverBar;}{K}}_k{H}_k){\widetilde{\stackrel{\&OverBar;}{p}}}_k\end{array}\right.$

式中：表示系统噪声协方差，表示量测噪声协方差，并且有

$\begin{array}{c}{\stackrel{\&OverBar;}{Q}}_k=({\stackrel{\&OverBar;}{x}}_{k+1}-f({\widehat{\stackrel{\&OverBar;}{x}}}_k,u_k)+\underset{j=1}{\overset{k+1}{\mathrm{\&Sigma;}}}{(-1)}^j{\mathrm{\&gamma;}}_j{\widehat{\stackrel{\&OverBar;}{x}}}_{k+1-j}){({\stackrel{\&OverBar;}{x}}_{k+1}-f({\widehat{\stackrel{\&OverBar;}{x}}}_k,u_k)+\underset{j=1}{\overset{k+1}{\mathrm{\&Sigma;}}}{(-1)}^j{\mathrm{\&gamma;}}_j{\widehat{\stackrel{\&OverBar;}{x}}}_{k+1-j})}^T\\ +(F_k+{\mathrm{\&gamma;}}_1){\widehat{\stackrel{\&OverBar;}{p}}}_k{(F_k+{\mathrm{\&gamma;}}_2)}^T+\underset{j=2}{\overset{k+1}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_2{\widehat{\stackrel{\&OverBar;}{p}}}_{k+1-j}{\mathrm{\&gamma;}}_j^T\end{array}$

${\stackrel{\&OverBar;}{R}}_k=[{\stackrel{\&OverBar;}{y}}_k-h\left({\widetilde{\stackrel{\&OverBar;}{x}}}_k\right)]{[{\stackrel{\&OverBar;}{y}}_k-h\left({\widetilde{\stackrel{\&OverBar;}{x}}}_k\right)]}^T+H_k{\widetilde{\stackrel{\&OverBar;}{p}}}_k{H}_k^T$

下面介绍本发明的一个实施例：

考虑非高斯Lévy噪声的分数阶非线性系统模型：

$\left\{\begin{array}{c}{\mathrm{\&Delta;}}^{\mathrm{\&alpha;}}{\stackrel{\&OverBar;}{x}}_{k+1}=f({\stackrel{\&OverBar;}{x}}_k,u_k)+{\stackrel{\&OverBar;}{w}}_k\\ {\stackrel{\&OverBar;}{x}}_{k+1}={\mathrm{\&Delta;}}^{\mathrm{\&alpha;}}{\stackrel{\&OverBar;}{x}}_{k+1}-\underset{j=1}{\overset{k+1}{\mathrm{\&Sigma;}}}{(-1)}^j{\mathrm{\&gamma;}}_j{\stackrel{\&OverBar;}{x}}_{k+1-j}\\ {\stackrel{\&OverBar;}{y}}_k=h\left({\stackrel{\&OverBar;}{x}}_k\right)+{\stackrel{\&OverBar;}{v}}_k\end{array}\right.$

式中

$f({\stackrel{\&OverBar;}{x}}_k,u_k)=\left[\begin{array}{c}f({\stackrel{\&OverBar;}{x}}_{1,k},{\stackrel{\&OverBar;}{u}}_{1,k})\\ f({\stackrel{\&OverBar;}{x}}_{2,k},{\stackrel{\&OverBar;}{u}}_{2,k})\\ f({\stackrel{\&OverBar;}{x}}_{3,k},{\stackrel{\&OverBar;}{u}}_{3,k})\end{array}\right]=\left[\begin{array}{c}{\stackrel{\&OverBar;}{x}}_{1.k}\cos \left({\stackrel{\&OverBar;}{x}}_{2,k}\right)\\ -0.6\sin \left({\stackrel{\&OverBar;}{x}}_{1,k}\right)+{\sin }^2\left({\stackrel{\&OverBar;}{x}}_{2,k}\right)-0.1x_{3,k}\\ 0.4{\stackrel{\&OverBar;}{x}}_{1,k}{\cos }^2\left({\stackrel{\&OverBar;}{x}}_{3,k}\right)\end{array}\right]+c_k{u}_k$

$h\left({\stackrel{\&OverBar;}{x}}_k\right)=\left[\begin{array}{c}0.1\sin \left({\stackrel{\&OverBar;}{x}}_{1,k}\right)-0.2\cos \left({\stackrel{\&OverBar;}{x}}_{2,k}\right)\\ 0.15{\stackrel{\&OverBar;}{x}}_{1,k}\cos \left({\stackrel{\&OverBar;}{x}}_{3,k}\right)\end{array}\right]$

c_k＝[0.2 0.2 0.2]^T

u_k＝1

α＝[0.3 1.2 0.8]^T

运用本发明所提到的状态更新方法进行仿真，结果如图2所示；真实值与估计值之间的比较如图3所示。

Claims (1)

1.一种分数阶非线性系统状态估计方法，其特征在于，包括以下步骤：

(1)、初始化，包括：设定状态预测量的初值和预测误差协方差的初值；

(2)、对非高斯Lévy噪声进行近似处理，推导得出状态近似值和量测输出近似值，计算步骤为：

${\stackrel{\&OverBar;}{x}}_{k+1}=\left\{\begin{array}{cc}{\mathrm{\&Omega;}}_1+{\mathrm{\&delta;}}_1\&CenterDot;\mathrm{sign}\left({\mathrm{\&Omega;}}_2\right)& \mathrm{if}\left|{\mathrm{\&Omega;}}_2\right|\&GreaterEqual;{\mathrm{\&delta;}}_1\\ {\stackrel{\&OverBar;}{x}}_{k+1}& \mathrm{if}\left|{\mathrm{\&Omega;}}_2\right|<{\mathrm{\&delta;}}_1\end{array}\right.$

${\stackrel{\&OverBar;}{y}}_k=\left\{\begin{array}{cc}h\left(\widetilde{\stackrel{\&OverBar;}{x}}\right)+{\mathrm{\&delta;}}_2\&CenterDot;\mathrm{sign}({\stackrel{\&OverBar;}{y}}_k-h\left(\widetilde{\stackrel{\&OverBar;}{x}}\right))& \mathrm{if}|{\stackrel{\&OverBar;}{y}}_k-h\left(\widetilde{\stackrel{\&OverBar;}{x}}\right)|{\&le;\mathrm{\&delta;}}_2\\ {\stackrel{\&OverBar;}{y}}_k& \mathrm{if}|{\stackrel{\&OverBar;}{y}}_k-h\left(\widetilde{\stackrel{\&OverBar;}{x}}\right)|<{\mathrm{\&delta;}}_2\end{array}\right.$

式中

${\mathrm{\&Omega;}}_1=f({\widetilde{\stackrel{\&OverBar;}{x}}}_k,u_k)-\underset{j=1}{\overset{k+1}{\mathrm{\&Sigma;}}}{(-1)}^j{\mathrm{\&gamma;}}_j{\widehat{\stackrel{\&OverBar;}{x}}}_{k+1-j}$

${\mathrm{\&Omega;}}_2={\stackrel{\&OverBar;}{x}}_{k+1}-f({\widetilde{\stackrel{\&OverBar;}{x}}}_k,u_k)+\underset{j=1}{\overset{k+1}{\mathrm{\&Sigma;}}}{(-1)}^j{\mathrm{\&gamma;}}_j{\widetilde{\stackrel{\&OverBar;}{x}}}_{k+1-j}$

表示含有Lévy噪声的状态向量，下标k+1表示第k+1时刻，δ₁表示所取的第一个阈值，sign(x)表示符号函数，如果x>0，返回1；如果x<0，返回-1；|x|表示取x的绝对值，表示量测输出向量，表示非线性函数，表示第k时刻状态预测向量，δ₂表示所取的第二个阈值，在Ω₁中，表示非线性函数，u_k表示第k时刻控制输入向量，γ_j表示

${\mathrm{\&gamma;}}_j=\mathrm{diag}(\left[\begin{array}{cccc}\left(\begin{array}{c}{\mathrm{\&alpha;}}_1\\ j\end{array}\right)& \left(\begin{array}{c}{\mathrm{\&alpha;}}_2\\ j\end{array}\right)& ...& \left(\begin{array}{c}{\mathrm{\&alpha;}}_N\\ j\end{array}\right)\end{array}\right]),$ 其中 $\left(\begin{array}{c}\mathrm{\&alpha;}\\ j\end{array}\right)=\left\{\begin{array}{cc}1& j=0\\ \frac{\mathrm{\&alpha;}(\mathrm{\&alpha;}-1)...(\mathrm{\&alpha;}-j+1)}{j!}& j>0\end{array}\right.$

α_N表示第N个分数阶阶次值，j＝1,2,…,k+1。

(3)、利用已得到的状态近似值和量测输出近似值，计算当前时刻系统噪声协方差和量测噪声协方差，计算步骤为：

$\begin{array}{c}{\stackrel{\&OverBar;}{Q}}_k=({\stackrel{\&OverBar;}{x}}_{k+1}-f({\widehat{\stackrel{\&OverBar;}{x}}}_k,u_k)+\underset{j=1}{\overset{k+1}{\mathrm{\&Sigma;}}}{(-1)}^j{\mathrm{\&gamma;}}_j{\widehat{\stackrel{\&OverBar;}{x}}}_{k+1-j}){({\stackrel{\&OverBar;}{x}}_{k+1}-f({\widehat{\stackrel{\&OverBar;}{x}}}_k,u_k)+\underset{j=1}{\overset{k+1}{\mathrm{\&Sigma;}}}{(-1)}^j{\mathrm{\&gamma;}}_j{\widehat{\stackrel{\&OverBar;}{x}}}_{k+1-j})}^T\\ +(F_k+{\mathrm{\&gamma;}}_1){\widehat{\stackrel{\&OverBar;}{p}}}_k{(F_k+{\mathrm{\&gamma;}}_1)}^T+\underset{j=2}{\overset{k+1}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_j{\widehat{\stackrel{\&OverBar;}{p}}}_{k+1-j}{{\mathrm{\&gamma;}}^T}_j\end{array}$

${\stackrel{\&OverBar;}{R}}_k=[{\stackrel{\&OverBar;}{y}}_k-h\left({\widetilde{\stackrel{\&OverBar;}{x}}}_k\right)]{[{\stackrel{\&OverBar;}{y}}_k-h\left({\stackrel{}{\stackrel{\&OverBar;}{x}}}_k\right)]}^T+H_k{\widetilde{\stackrel{\&OverBar;}{p}}}_k{H}_k^T$

式中，表示非线性函数，表示第k时刻状态估计向量，上标T表示转置，表示非线性函数在处的雅克比矩阵，表示第k时刻状态估计误差协方差，表示非线性函数在处的雅克比矩阵，表示第k时刻状态预测误差协方差；

(4)、利用当前时刻量测噪声协方差和给定的当前时刻预测误差协方差计算当前时刻最优滤波增益，计算步骤为：

${\stackrel{\&OverBar;}{K}}_k={\widetilde{\stackrel{\&OverBar;}{p}}}_k{H}_k^T{(H_k{\widetilde{\stackrel{\&OverBar;}{p}}}_k{H}_k^T+{\stackrel{\&OverBar;}{R}}_k)}^{-1}$

(5)、结合已得到的当前时刻最优滤波增益和给定的当前时刻状态预测值计算当前时刻状态估计值，计算步骤为：

${\widehat{\stackrel{\&OverBar;}{x}}}_k={\widetilde{\stackrel{\&OverBar;}{x}}}_k+{\stackrel{\&OverBar;}{K}}_k[{\stackrel{\&OverBar;}{y}}_k-h\left({\widetilde{\stackrel{\&OverBar;}{x}}}_k\right)]$

(6)、利用当前时刻最优滤波增益和给定的当前时刻预测误差协方差计算当前时刻估计误差协方差，计算步骤为：

${\widehat{\stackrel{\&OverBar;}{p}}}_k=(I-{\stackrel{\&OverBar;}{K}}_k{H}_k){\widetilde{\stackrel{\&OverBar;}{p}}}_k$

(7)、利用当前时刻状态估计值更新下一时刻状态预测值，计算步骤为：

$\left\{\begin{array}{c}{\mathrm{\&Delta;}}^a{\widetilde{\stackrel{\&OverBar;}{x}}}_{k+1}=f({\widehat{\stackrel{\&OverBar;}{x}}}_k,u_k)\\ {\widetilde{\stackrel{\&OverBar;}{x}}}_{k+1}={\mathrm{\&Delta;}}^{a{\widetilde{\stackrel{\&OverBar;}{x}}}_{k+1}}-\underset{j=1}{\overset{k+1}{\mathrm{\&Sigma;}}}{(-1)}^j{\mathrm{\&gamma;}}_j{\widehat{\stackrel{\&OverBar;}{x}}}_{k+1-j}\end{array}\right.$

式中Δ表示分数阶算子，上标α表示分数阶阶次向量。

(8)、利用当前时刻估计误差协方差和系统噪声协方差更新下一时刻状态预测协方差，计算步骤为：

${\widetilde{\stackrel{\&OverBar;}{p}}}_{k+1}=(F_k+{\mathrm{\&gamma;}}_1){\widehat{\stackrel{\&OverBar;}{p}}}_k{(F_k+{\mathrm{\&gamma;}}_1)}^T+{\stackrel{\&OverBar;}{Q}}_k+\underset{j=2}{\overset{k+1}{\mathrm{\&Sigma;}}}{\mathrm{\&gamma;}}_j{\widehat{\stackrel{\&OverBar;}{p}}}_{k+1-j}{\mathrm{\&gamma;}}_j^T$

(9)、判断k+1是否大于等于步长L，如果是，结束计算，否则返回步骤(2)进行下一次估计。