主要内容

基于粒子滤波块的Simulink参数和状态估计金宝app

这个例子演示了在系统辨识工具箱™使用粒子滤波块。离散时间传递函数参数估计问题重新并递归地解决作为状态估计问题。

介绍

系统辨识工具箱有三个Simulink模块用于非线性状态估计:金宝app

  • 颗粒过滤器:实现了离散时间粒子滤波算法。

  • 扩展卡尔曼滤波:实现一阶离散时间扩展卡尔曼滤波算法。

  • 无迹卡尔曼滤波:实现离散时间无迹卡尔曼滤波算法。

这些块支持使用在不同的采样金宝app速率操作的多个传感器的状态估计。用于使用这些块的典型工作流程如下:

  1. 建模您的工厂,并使用MATLAB或Simulink的功能传感器的行为。金宝app

  2. 配置块的参数。

  3. 模拟滤波器和分析结果,以获得在过滤性能的信心。

  4. 在硬件上部署过滤器。您可以使用Simulink Coder™软件为这些过滤器生成代码。金宝app

本例使用粒子滤波块来演示此工作流程的前两个步骤。最后两个步骤在简要讨论下一步部分。本例中的目标是递归地估计离散时间传递函数(输出误差模型)的参数,其中模型参数在每个时间步骤中随着新信息的到来而更新。

如果你有兴趣在扩展卡尔曼滤波器,看到的例子“非线性系统的评估各国多,多速率传感器”。采用无迹卡尔曼滤波的如下类似的步骤,以扩展卡尔曼滤波。

示例文件夹添加到MATLAB路径。

让addpath(完整文件(matlabroot,“例子”'IDENT'“主要”));

植物建模

大多数状态估计算法依赖于描述从一个时间步到下一个设备状态演变的工厂模型(状态转换功能)。该函数通常表示为$ X [K + 1] = F(X [k]的,W [k]的,U [K])$其中x是状态,w是过程噪声,u是任选的另外的输入,例如系统输入或参数。颗粒过滤器块需要你在一个稍微不同的语法提供这个功能,$ X [k + 1] = f {pf} u (X [k], [k])美元.区别是:

  • 微粒过滤器的工作原理是以下许多态假说(颗粒)的轨迹,和块一次通过了所有状态的假设给你的函数。具体而言,如果你的状态向量x美元拥有N_s美元元素和你选择N_p美元使用粒子,X美元有尺寸美元[N_s \;N_p]美元其中每列是一个状态的假设。

  • 计算过程噪声的影响w美元对态假说$ X [K + 1] $在你的函数$ {F_ PF}(...)$.该块不就过程噪声的概率分布的任何假设w美元,也不需要w美元作为输入。

这个函数$ {F_ PF}(...)$可以是一个MATLAB函数与MATLAB编码器™或Simulink的功能块的限制规定。金宝app在创建$ {F_ PF}(...)$,您可以指定在颗粒过滤器块的函数名。

在此示例中,将重整离散时间传递函数参数估计问题作为一个状态估计问题。该传递函数可以表示离散时间过程的动力学,或者它可以表示加上的信号重构器如零阶保持一些连续时间动态。假设你有兴趣在估计第一阶离散时间传递函数的参数:

$ $ y [k] = \压裂{20问^ {1}}{1 - 0.7 - q ^ {1}} u e [k] + [k] $ $

这里美元$ y [k]为植物产量,$ U [K] $是工厂输入,$ E [k]的$是测量噪声,$ Q ^ { -  1} $是延时操作,使$ Q ^ { -  1} U [K] = U [T-1] $.参数化的传递函数$$ \压裂{N \; Q ^ { -  1}} {1 + d \; Q ^ { -  1}} $$, 在哪里$ n $$ d $是待估计的参数。通过选择状态向量,传递函数和参数可以多种方式表示为必要的状态空间形式。一个选择是$ X [K] = [\;Y [k]的;d [K];N [k]的\;] $其中,第二和第三状态表示参数估计值。然后传递函数可以等效写成

$$ x[k+1] = \left[
c \开始{数组}{}& # xA;-x_2[k] [k] + x_3[k] [k] [k] [k];x_2 [k] \ \ & # xA;x_3 [k] \ \ & # xA;结束\{数组}& # xA;\] $ $

测量噪声项$ E [k]的$在传感器建模处理。在本例中,执行上述表达式的MATLAB函数,在用于计算效率的矢量形式:

类型pfBlockStateTransitionFcnExample
功能xNext = pfBlockStateTransitionFcnExample(X,U)%pfBlockStateTransitionFcnExample状态转换为粒子%滤光函数,为SISO,一阶,离散时间传递%函数模型%%输入估计的%参数:%X  - 颗粒,NumberOfStates逐-NumberOfParticles矩阵%U  - 系统输入,一个标量%%输出:%xNext  - 预计颗粒,具有相同的尺寸,输入x%实施状态转换函数%xNext = [X(1)* X(2)+ X(3)* U;%×(2);%×(3)];%在向量化形式(对于所有颗粒)。xNext = X;xNext(1,:) = bsxfun(@倍,X(1,:), -  X(2,:))+ X(3,:)* U;%添加一个小的过程噪声(相对于每一个状态的预期大小),以%增加粒子多样性xNext = xNext + bsxfun(@times,[1; 1E-2; 1E-1],randn(大小(xNext)));结尾

传感器建模

颗粒过滤器块需要你提供,其计算每个状态假设的可能性(概率)的测量似然函数。该函数的形式$ L [K] = {H_ PF}(X [K],Y [k]的,U [K])$$ L [k]的$是一个N_p美元元素向量,其中N_p美元是您选择的粒子数。$ M ^ {第} $元素$ L [k]的$可能性是$ M ^ {第} $粒子(列)中美元$ X [k]美元$ y [k]是传感器测量。$ U [K] $是一个可选的输入参数,它可以从状态转换函数的输入是不同的。

所述传感器测量在这个例子中,第一状态。本示例假定实际和预测测量值之间的误差是根据高斯分布分散的,但是任何任意的概率分布或一些其他方法可以被用来计算似然性。创建$ {H_ PF}(...)$,并在Particle Filter块中指定函数名。

类型pfBlockMeasurementLikelihoodFcnExample;
函数= pfBlockMeasurementLikelihoodFcnExample可能性(粒子,测量)% pfBlockMeasurementLikelihoodFcnExample测量粒子滤波% %的似然函数是第一个国家输入% %:% - % NumberOfStates-by-NumberOfParticles矩阵测量粒子系统输出,一个标量% %输出:% likelihood -一个具有NumberOfParticles元素的向量,其第n %元素是第n个粒子的可能性%#codegen % Predicted measurement yHat = particles(1,:);假设误差分布在每个多元正态分布上,%均值为0,方差为1。计算相应的概率%密度函数e = bsxfun(@minus, yHat, measurement(:)');% Error numberOfMeasurements = 1;μ= 0;%均值=眼(测量数);% Variance measurementErrorProd = dot((e-mu), Sigma \ (e-mu), 1);c = 1/√((2*pi)^numberOfMeasurements * det(Sigma));可能性= c * exp(-0.5 * measureenterrorprod); end

滤波器结构

配置估计颗粒过滤器块。指定的状态转换和测量似然函数名,粒子数,并且这些粒子的初始分布。

系统模型块对话框的选项卡中,指定以下参数:

状态转换

  1. 指定状态转换功能,pfBlockStateTransitionFcnExample,在函数.当您输入函数名并单击时申请,块检测到您的函数有一个额外的输入,你美元,并创建输入端口StateTransitionFcnInputs.您将系统输入连接到这个端口。

初始化

  1. 指定10000粒子数.较高的粒子数目通常对应于更好的估计,在增加的计算成本。

  2. 指定高斯分配从多变量高斯分布得到初始的一组颗粒。然后指定[0;0;0]的意思是因为你有三种状态,这是你最好的猜测。指定诊断([10 5 100])协方差以指定您的猜测第三状态的一大变化(不确定性),以及更小的方差为第2位。至关重要的是,这个初始组颗粒散布足够宽(大方差)以覆盖潜在的真实状态。

测量1

  1. 指定测量可能性函数的名称,pfBlockMeasurementLikelihoodFcnExample,在函数

采样时间

  1. 在块对话框的底部,在进入1采样时间.如果有中间状态转换和测量似然函数不同的采样时间,或者如果你有具有不同采样时间的多个传感器,这些可以在所配置的块输出,多重速率的标签。

颗粒过滤器包括删除与低似然性粒子和幼苗用的那些具有较高似然性的新的粒子。这是由下该选项控制重采样组。本示例使用默认设置。

默认情况下,区块只输出状态假设的平均值,并根据它们的可能性加权。要查看所有的粒子、权重,或选择提取状态估计的不同方法,请查看块输出,多重速率的标签。

仿真和结果

对于一个简单的测试,真正的工厂模型模拟与白噪声输入。的输入和从植物中的噪声测量被馈送到粒子滤波块。以下Simulink模型表示金宝app此设置。

open_system (“pfBlockExample”

模拟系统,比较估计的和真实的参数:

SIM(“pfBlockExample”)open_system(“pfBlockExample /参数估计”

图中显示了真实的分子和分母参数,以及它们的粒子滤波估计。经过10个时间步长,估计值近似收敛于真值。即使初始状态猜想与真实值相差甚远,也能获得收敛性。

故障排除

一些潜在的实施问题和故障排除的想法都在这里上市,情况如预期的应用程序中的微粒过滤器不执行。

粒子滤波器的故障排除通常是通过查看粒子集和它们的权重来执行的,这些权重可以通过选择得到输出所有粒子输出所有的权重块输出,多重速率的块对话框的选项卡。

第一个完整性检查是确保状态转换和度量可能性函数能够很好地捕获系统的行为。如果您的系统有一个仿真模型(因此可以访问仿真中的真实状态),那么可以尝试使用真实状态初始化过滤器。然后,您可以验证状态转移函数是否准确地计算了真实状态的时间传播,以及测量似然函数是否为这些粒子计算了很高的可能性。

初始粒子集很重要。在你的模拟开始时,确保至少有一些粒子有很高的可能性。如果真实状态在状态假设的初始传播范围之外,估计可能是不准确的,甚至会出现分歧。

如果状态估计精度最初很好,但随着时间的推移而恶化,问题可能是粒子退化或粒子枯竭[2]。当粒子分布太广时,粒子会发生简并;当粒子重新采样后,粒子会聚集在一起,粒子会发生贫化。由于直接重采样,粒子退化导致粒子贫化。在本例中使用的状态转移函数中加入人工过程噪声是一种实用的方法。有大量关于解决这些问题的文献,并且基于您的应用程序可以使用更系统的方法。[1]和[2]是两个有用的参考文献。

下一步

  1. 验证状态估计:一旦滤波器在模拟中按预期执行,通常使用广泛的蒙特卡罗模拟进一步验证性能。有关更多信息,请参见验证在线状态估计在Simulink金宝app.您可以使用选项下随机性在颗粒过滤器块对话组以促进这些模拟。

  2. 生成代码:在颗粒过滤器块支持利用Simulink编码器™软件C和C ++代码生成。金宝app金宝app到该块中,所提供的功能,必须符合MATLAB编码器™软件的限制,(如果你使用MATLAB功能,您的系统模型)和Simulink编码器软件(如果你使用的Simulink功能块到您的系统模型)。金宝app

总结

这个例子展示了如何在系统识别工具箱中使用粒子过滤器块。你递归地估计一个离散时间传递函数的参数,其中参数在每个时间步随着新信息的到来而更新。

close_system (“pfBlockExample”,0)

从MATLAB路径中删除示例文件夹。

rmpath (fullfile (matlabroot,“例子”'IDENT'“主要”));

参考文献

[1]西蒙,丹。最优状态估计:卡尔曼,H无穷大,和非线性方法。John Wiley&Sons,2006年。

[2]杜塞,阿诺,和Adam M.约翰森。“粒子的指南过滤和平滑:十五年后”。非线性滤波12.656-704手册(2009):3。

也可以看看

||

相关的话题