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

这个示例使用:

打开脚本

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

简介

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

粒子滤波:实现离散时间粒子滤波算法。
扩展卡尔曼滤波:实现一阶离散时间扩展卡尔曼滤波算法。
无气味卡尔曼滤波:实现离散时间无气味卡尔曼滤波算法。

这些块支持使用在不同采样率下运行的多个传感器进行状态估计。使用这些块的典型工作流程如下:

使用MATLAB或Simulink函数为植物和传感器行为建模。
配置块的参数。
模拟滤波器并分析结果，以获得对滤波器性能的信心。
在硬件上部署过滤器。您可以使用Simulink Coder™软件为这些过滤器生成代码。

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

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

将示例文件夹添加到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])美元．的差异是:

粒子滤波器通过跟踪许多状态假设(粒子)的轨迹工作，并且块将所有的状态假设立即传递给你的函数。具体来说，如果你的状态向量有元素，你选择使用粒子,有尺寸 $美元[N_s \;N_p]美元$ 每一列都是一个状态假设。
计算过程噪声的影响关于状态假设在你的函数．该块对过程噪声的概率分布不做任何假设，不需要作为输入。

这个函数 f {pf}(…)美元美元可以是一个符合MATLAB Coder™的MATLAB函数，也可以是一个Simulink函数块。在您创建，在粒子筛选块中指定函数名。

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

$$ 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] 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], [k], u [k])美元$ ． L [k]美元是一个 N_p美元元素的向量,是你选择的粒子数。 $m ^ {th}$ 元素的可能性是 $m ^ {th}$ 粒子(列) 美元$ X [k] ．美元$ y [k] 是传感器测量。 u [k]美元是一个可选的输入参数，它可以不同于状态转换函数的输入。

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

类型pfBlockMeasurementLikelihoodFcnExample

pfblockmeasurementlikehoodfcnexample(粒子，测量)% pfblockmeasurementlikehoodfcnexample测量似是函数粒子过滤器% %测量是第一状态% %输入:%粒子-一个NumberOfStates-by-NumberOfParticles矩阵%测量-系统输出，一个标量% %输出:% likelihood -一个含有NumberOfParticles元素的向量，其第n %元素是第n个粒子的可能性%#codegen %预测测量yHat =粒子(1，:);假设误差按多元正态分布分布，均值为0，方差为1。计算相应的概率%密度函数e = bsxfun(@minus, yHat, measurement(:)');%误差数= 1;μ= 0;%平均西格玛=眼睛(numberOfMeasurements);% Variance measurementErrorProd = dot((e-mu)， Sigma \ (e-mu)， 1);c = 1/√((2*pi)^numberOfMeasurements * det(Sigma));likelihood = c * exp(-0.5 * measurementErrorProd); end

滤波器结构

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

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

状态转换

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

初始化

指定10000在粒子数．较高的粒子数通常对应更好的估计，增加计算成本。
指定高斯在分布从多元高斯分布中得到初始粒子集。然后指定[0;0;0]在的意思是因为你有三种状态，这是你最好的猜测。指定诊断接头([100]5)在协方差为第三种状态指定较大的方差(不确定性)，为前两种状态指定较小的方差。至关重要的是，这一初始粒子集要传播得足够广(大的方差)，以覆盖潜在的真实状态。

测量1

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

样品时间

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

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

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

仿真和结果

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

open_system (“pfBlockExample”）

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

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

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

故障排除

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

粒子滤波器的故障排除通常通过查看粒子集及其权重来完成，这些权重可以通过选择获得输出所有粒子而且输出所有重量在块输出,多重速率的选项卡的块对话框。

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

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

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

下一个步骤

验证状态估计:一旦滤波器在模拟中按预期执行，通常使用广泛的蒙特卡洛模拟进一步验证性能。有关更多信息，请参见在Simulink中验证在线状态估计．您可以使用下面的选项随机性组在粒子过滤器块对话框，以促进这些模拟。
生成代码:粒子过滤器块支持使用Simulink Coder™软件生成C和c++代码。您提供给此块的函数必须符合MATLAB Coder™软件(如果您使用MATLAB函数对您的系统建模)和Simulink Coder软件(如果您使用Simulink函数块对您的系统建模)的限制。

总结

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

close_system (“pfBlockExample”, 0)

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

rmpath (fullfile (matlabroot,“例子”，“控制”，“主要”));

参考文献

[1]西门,丹。最优状态估计:卡尔曼，H∞和非线性方法。约翰·威利父子公司，2006年。

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

另请参阅

扩展卡尔曼滤波器|无味卡尔曼滤波|粒子滤波