主要内容

基于粒子滤波块的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 Coder™限制的MATLAB函数，也可以是Simulink函数块。金宝app创建之后f {pf}(…)美元美元，在粒子过滤器块中指定函数名。

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

$ $ 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]美元．将传递函数参数化为$$ \frac{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] x_1[k] + x_3[k] u[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 %实现状态转换函数% xNext = [x (1) * (2) + x (3) * u;% x (2);% x (3)];%的向量化形式(所有粒子)。xNext = x;xNext (1) = bsxfun (@times x (1:) - x (2:)) + x (3:) * u;%添加一个小的过程噪声(相对于每个状态的预期大小)，以%增加粒子多样性xNext = xNext + bsxfun(@times，[1;1飞行;1 e 1], randn(大小(xNext)));结束

传感器建模

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

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

类型pfBlockMeasurementLikelihoodFcnExample；
function likelihood = pfblockmeasurementlikehoodfcnexample(粒子，测量)% pfblockmeasurementlikehoodfcnexample测量似然函数用于粒子过滤器% %测量是第一个状态% %输入:%粒子-一个numberofstats -by- numberofparticles矩阵% measurement -系统输出，一个标量% %输出:% likelihood -一个含有NumberOfParticles元素的向量，其中第n个%元素是第n个粒子的可能性%#编码原%预测测量yHat =粒子(1，:);根据实际测量值%与预测测量值%之间的误差计算每个粒子的似然性，假设误差分布在多元正态分布中，%平均值为零，方差为1。求对应的概率%密度函数e = bsxfun(@minus, yHat, measurement(:)');%误差数为1;Mu = 0;%平均西格玛=眼睛(数量的测量);%方差测量enterprod = dot((e-mu)， Sigma \ (e-mu)， 1);c = 1/√((2*pi)^numberOfMeasurements * det(Sigma));似然= c * exp(-0.5 * measurementErrorProd); end

滤波器结构

配置粒子过滤块进行估计。您可以指定状态转换和测量似然函数名称、粒子数以及这些粒子的初始分布。

在系统模型选项卡，指定以下参数:

状态转换

  1. 指定状态转换函数，pfBlockStateTransitionFcnExample,在函数．输入函数名后，单击应用，块检测到你的函数有一个额外的输入，你美元，并创建输入端口StateTransitionFcnInputs．您将系统输入连接到此端口。

初始化

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

  2. 指定高斯在分布从多元高斯分布中得到初始粒子集。然后指定[0;0;0]在的意思是因为你有三种状态，这是你的最佳猜测。指定Diag ([10 5 100])在协方差为第三个状态指定较大的方差(不确定性)，为前两个状态指定较小的方差。至关重要的是，这一初始粒子集的分布足够广泛(大方差)，以覆盖潜在的真实状态。

测量1

  1. 指定测量似然函数的名称，pfBlockMeasurementLikelihoodFcnExample,在函数．

样品时间

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

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

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

模拟及结果

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

open_system (“pfBlockExample”）

对系统进行仿真，比较估计参数和真实参数:

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

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

故障排除

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

粒子滤波器的故障排除通常通过查看粒子集及其权重来执行，可以通过选择来获得输出所有粒子而且输出所有权重在块输出，多速率块对话框的标签。

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

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

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

下一个步骤

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

  2. 生成代码:Particle Filter模块支持使用Simulink Coder™金宝app软件生成C和c++代码。金宝app向此块提供的函数必须符合MATLAB Coder™软件(如果您正在使用MATLAB函数对系统建模)和Simulink Coder软件(如果您正在使用Simulink函数块对系统建模)的限制。金宝app

总结

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

close_system (“pfBlockExample”, 0)

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

rmpath (fullfile (matlabroot,“例子”，“识别”，“主要”));

参考文献

西蒙，丹。最优状态估计:卡尔曼，H∞和非线性方法。约翰·威利父子出版社，2006年出版。

[2]杜塞，阿诺，亚当·m·约翰森。《粒子滤波和平滑教程:15年后》非线性滤波手册12.656-704(2009):3。

另请参阅

||

相关的话题

系统识别工具箱文档

金宝app

试试MATLAB、Sim金宝appulink和其他产品下载188bet金宝搏

现在就去审判吧