主要内容

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

此示例使用：

打开脚本

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

介绍

该控制系统工具箱具有非线性状态估计3个Simulink模块：金宝app

  • 粒子过滤器：实现离散时间粒子过滤器算法。

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

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

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

  1. 使用MATLAB或Simulink函数对设备和传感器行为建模。金宝app

  2. 配置块的参数。

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

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

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

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

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

addpath（完整文件（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}（...）$，可以在粒子过滤器块中指定函数名称。

在本例中，您将离散时间传递函数参数估计问题重新表述为状态估计问题。此传递函数可能表示离散时间过程的动力学，也可能表示与信号重构器（如零阶保持）耦合的某些连续时间动力学。假设您有兴趣估计一阶离散时间传递函数的参数：

$ $ 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=PFBlockStateTransitionFCN示例（x，u）%PFBlockStateTransfc粒子%过滤器的状态转换函数示例，用于估计%SISO、一阶、离散时间传递%函数模型%%的参数输入：%x-粒子、NumberOfStates by NumberOfParticles矩阵%u-系统输入、标量%%输出：%x下一个-预测的粒子，具有相同的维数作为输入，x%以矢量化形式（针对所有粒子）实现状态转换函数%xNext=[x（1）*x（2）+x（3）*u；%x（2）；%x（3）]；%，[1；1e-2；1e-1]，随机编号（尺寸（xNext））；结束

传感器建模

粒子过滤器块要求您提供一个测量似然函数，用于计算每个状态假设的似然（概率）。此函数具有以下形式：$ 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}（…）$，nd指定粒子过滤器块中的函数名。

类型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]）在里面协方差为第三种状态指定较大的偏差（不确定性），为前两种状态指定较小的偏差。至关重要的是，这组初始粒子的分布要足够宽（大的差异），以覆盖潜在的真实状态。

测量1

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

采样时间

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

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

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

仿真和结果

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

open_system (“pfBlockExample”）

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

SIM（“pfBlockExample”)开放式系统(“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,“例子”那“控制”那“主要”));

参考文献

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

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

另见

||

相关的话题

控制系统工具箱文档

金宝app

了解如何自动调整PID控制器增益

了解如何自动调整PID控制器增益

下载代码示例