参数和状态估计在Simulink中使用粒子滤波模块万博1manbetx

这个示例使用:

打开脚本

这个例子演示了粒子滤波块在控制系统工具箱中的使用。将离散时间传递函数参数估计问题转化为状态估计问题递推求解。

介绍

该控制系统工具箱具有非线性状态估计3个Simulink模块：万博1manbetx

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

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

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

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

  1. 使用MATLAB或Simulink函数对植物和传感器的行为建模。万博1manbetx

  2. 配置块的参数。

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

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

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

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

将示例文件文件夹添加到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的功能块的限制规定。万博1manbetx在创建$ {F_ PF}（...）$，在粒子滤波块中指定函数名。

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

$$ Y [k]的= \压裂{20Q ^ { - 1}} {1-0.7q ^ { - 1}} U [K] + E [k]的$$

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

$$ x[k+1] = \左[ 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)));结束

传感器建模

粒子滤波块要求您提供一个测量似然函数，计算每个状态假设的似然(概率)。这个函数有这种形式$ 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
函数= pfBlockMeasurementLikelihoodFcnExample可能性(粒子,测量)% pfBlockMeasurementLikelihoodFcnExample测量粒子滤波% %的似然函数是第一个国家输入% %:% - % NumberOfStates-by-NumberOfParticles矩阵测量粒子系统输出,一个标量% %输出:%似然——一个带有粒子数元素的向量，其第n个元素是第n个粒子的似然值%#codegen %预测测量yHat =粒子(1，:);根据实际测量值和预测测量值之间的误差%计算每个粒子的似然值%假设误差按多元正态分布，%的均值为0，方差为1。求对应的概率%密度函数e = bsxfun(@ -， yHat, measurement(:)');测量误差% = 1;μ= 0;平均西格玛=眼(测量数);%方差测量enterrorprod = dot((e-mu)， Sigma \ (e-mu)， 1);c = 1/sqrt((2*pi)^测量数* det(Sigma));似然= c * exp(-0.5 * measure rementerrorprod); end

过滤建设

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

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

状态转换

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

初始化

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

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

测量1

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

采样时间

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

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

默认情况下，该块仅输出状态假设，可以通过似然性加权的平均值。要看到所有的粒子，重量，还是要选择提取状态估计的不同方法，请在选项块输出,多重速率的选项卡。

仿真和结果

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

open_system (“pfBlockExample”）

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

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

该图显示了真正的分子和分母的参数，以及它们的粒子滤波估计。估计大约收敛于真实值后的10个时间步。收敛，即使初始状态的猜测是远离真实值获得。

故障排除

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

粒子滤波的故障排除通常是通过观察粒子的集合及其权值来完成的，通过选择权值可以得到粒子的权值输出所有颗粒和输出所有的权重在块输出,多重速率的块对话框的选项卡。

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

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

如果状态估计精度精细最初，但随着时间的推移不断恶化，该问题可以是粒子退化或粒子贫瘠[2]。当颗粒被分布过于广泛，而粒子贫瘠的发生是因为重采样之后结块在一起的粒子发生粒子退化。粒子退化导致颗粒贫困直接采样的结果。添加在本例中使用的状态转换函数人工过程噪声是一个实用的方法。有上解决这些问题，并根据您的应用程序更系统的方法可能提供了大量馆藏文献的。[1]，[2]是两个引用，可以是有帮助的。

下一个步骤

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

  2. 生成代码：在颗粒过滤器块支持利用Simulink编码器™软件C和C ++代码生成。万博1manbetx万博1manbetx到该块中，所提供的功能，必须符合MATLAB编码器™软件的限制，（如果你使用MATLAB功能，您的系统模型）和Simulink编码器软件（如果你使用的Simulink功能块到您的系统模型）。万博1manbetx

摘要

此示例示出了如何使用在控制系统工具箱颗粒过滤器块。您估计递归离散时间传递函数，其中的参数在每个时间步长作为新信息到达更新的参数。

close_system (“pfBlockExample”，0）

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

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

参考

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

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

另请参阅

||

相关话题

控制系统工具箱文档
万博1manbetx
学习如何自动调整PID控制器增益

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

下载代码示例