主要内容

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

这个示例使用:

打开脚本

这个例子演示了在控制系统工具箱™中粒子过滤器块的使用。将离散时间传递函数参数估计问题重新表述为状态估计问题并递归求解。

介绍

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

  • 粒子滤波:实现离散时间粒子滤波算法。

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

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

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

  1. 使用MATLAB或Simulink函数建模您的植物和传感器行为。金宝app

  2. 配置块的参数。

  3. 模拟滤波器并分析结果，以获得对滤波器性能的置信度。

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

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

如果您对扩展卡尔曼滤波器感兴趣，请参阅示例“用多个、多速率传感器估计非线性系统的状态”。Unscented卡尔曼滤波器的使用遵循与扩展卡尔曼滤波器相似的步骤。

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

目录(fullfile (matlabroot,“例子”，“控制”，“主要”））；

植物建模

大多数状态估计算法依赖于描述植物状态从一个时间步骤到下一个时间步骤的植物模型(状态转移函数)。这个函数通常表示为$ 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}(…)美元美元，则在Particle Filter块中指定函数名。

在这个例子中，您将离散时间传递函数参数估计问题重新表述为状态估计问题。这个传递函数可以表示离散时间过程的动力学，也可以表示与信号重构器(如零阶保持器)耦合的连续时间动力学。假设你有兴趣估计一阶离散时间传递函数的参数:

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

在这里美元$ y [k]为植物产量，u [k]美元为植物的输入，美元$ e [k]为测量噪声，美元问^ {1}$时延算子是这样的吗美元问^ {1}u [k] = [t - 1]美元．参数化传递函数为$ $ \压裂{n \;问^ {1}}{1 + d \;问^ {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粒子状态转换函数%过滤器,估算参数%的输出,一阶,离散时间传递% % %函数模型输入:% x -粒子,NumberOfStates-by-NumberOfParticles矩阵% u -系统输入,一个标量% %输出:实现状态转移函数% xNext = [x(1)*x(2) + x(3)*u;% x (2);% x (3)];%在矢量形式(为所有粒子)。xNext = x;xNext (1) = bsxfun (@times x (1:) - x (2:)) + x (3:) * u;添加一个小的过程噪声(相对于每个状态的预期大小)，以增加粒子多样性1飞行;1 e 1], randn(大小(xNext))); end

传感器建模

Particle Filter块要求您提供一个度量可能性函数，用于计算每个状态假设的可能性(概率)。这个函数有它的形式$ L [k] = h_ {pf} (X [k], [k], u [k])美元．L [k]美元是一个N_p美元元素的向量,N_p美元是你选择的粒子数。$ m ^ {th} $元素L [k]美元可能性是$ m ^ {th} $粒子(列)美元$ X [k]．美元$ y [k]为传感器测量。u [k]美元是一个可选输入参数，它可以不同于状态转移函数的输入。

在本例中，传感器测量第一种状态。这个例子假设实际测量和预测测量之间的误差是按照高斯分布分布的，但任何任意概率分布或其他方法都可以用来计算可能性。您创建美元$ h_ {pf}(…)，并在粒子过滤器块中指定函数名。

类型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]在的意思是因为你有三个状态，这是你的最佳猜测。指定诊断接头([100]5)在协方差为第三种状态指定较大的方差(不确定性)，为前两种状态指定较小的方差。至关重要的是，这一初始粒子集的分布范围足够广(大的方差)，以覆盖潜在的真实状态。

测量1

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

样品时间

  1. 在区块对话框的底部，输入1样品时间．如果在状态转移和测量似然函数之间有不同的采样时间，或者如果你有多个不同采样时间的传感器，这些可以在块输出,多重速率的选项卡。

粒子过滤器包括删除概率低的粒子，并使用概率高的粒子播种新的粒子。控件下的选项控制重采样组。本示例使用默认设置。

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

仿真和结果

为了简单的测试，在白噪声输入下对真植物模型进行了模拟。输入和噪声测量从工厂馈送粒子滤波器块。下面的Simulink模型表金宝app示了这个设置:

open_system (“pfBlockExample”）

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

sim卡(“pfBlockExample”) open_system (pfBlockExample /参数估计的）

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

故障排除

这里列出了一些潜在的实现问题和故障排除思路，以防粒子过滤器没有按照您的应用程序的预期执行。

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

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

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

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

下一个步骤

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

  2. 生成代码:粒子过滤器块支持使用Simulink Coder™软件生成C和c++代码金宝app。金宝app您提供给此块的函数必须符合MATLAB Coder™软件(如果您使用MATLAB函数来建模您的系统)和Simulink Coder软件(如果您使用Simulink Function块来建模您的系统)的限制。金宝app

总结

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

close_system (“pfBlockExample”, 0)

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

rmpath (fullfile (matlabroot,“例子”，“控制”，“主要”））；

参考文献

[1]西门,丹。最优状态估计:卡尔曼、H∞和非线性方法。John Wiley & Sons, 2006。

Doucet, Arnaud和Adam M. Johansen。“粒子滤波和平滑教程:15年后。”非线性滤波手册12.656-704(2009):3。

另请参阅

||

相关的话题

控制系统工具箱文档

金宝app

学习如何自动调整PID控制器增益

学习如何自动调整PID控制器增益

下载代码示例