主要内容

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

这个示例使用:

打开脚本

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

介绍

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

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

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

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

这些模块支持使用在不同采样金宝app率下工作的多个传感器进行状态估计。使用这些块的典型工作流如下所示：

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

  2. 配置块的参数。

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

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

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

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

将示例文件夹添加到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]美元．参数化传递函数为$$\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] [k] + x_3[k] [k] [k] [k];x_2 [k] \ \ & # xA;x_3 [k] \ \ & # xA;结束\{数组}& # xA;\] $ $

测量噪声项$e[k]$是在传感器建模中处理的。在这个例子中，你在一个MATLAB函数中实现了上面的表达式，以矢量形式计算效率:

类型PfBlockStateTransitionFcExample
函数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]，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矩阵测量粒子系统输出,一个标量% %输出:% 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. 指定状态转移函数，PfBlockStateTransitionFcExample,在函数．当您输入函数名并单击时申请，该块检测到您的函数有一个额外的输入，你美元，并创建输入端口StateTransitionFcnInputs．您将系统输入连接到这个端口。

初始化

  1. 指定10000在粒子数．粒子数越高，估计越准确，计算成本也就越高。

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

测量1

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

采样时间

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

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

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

仿真和结果

对于一个简单的测试，使用白噪声输入模拟真实的电厂模型。来自设备的输入和噪声测量被馈送到粒子滤波器块。以下Simulink模型表示此设置：金宝app

open_system (“pfBlockExample”)

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

模拟(“pfBlockExample”) open_system (“pfBlockExample/Parameter Estimates”)

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

故障排除

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

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

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

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

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

下一个步骤

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

  2. 生成代码：粒子过滤器块支持C和C++代码生成，使用Simulink Coder™ 金宝app软件为该块提供的函数必须符合MATLAB Coder的限制™ 软件（如果使用MAT金宝appLAB函数对系统建模）和Simulink编码器软件（如果使用Simulink功能块对系统建模）。

总结

此示例演示了如何使用“控制系统工具箱”中的粒子过滤器块。您递归地估计离散时间传递函数的参数，当新信息到达时，参数会在每个时间步更新。

close_system (“pfBlockExample”, 0)

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

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

参考文献

[1] 西蒙，丹。最优状态估计：Kalman、H无穷大和非线性方法。约翰·威利父子公司，2006年。

[2] 《粒子滤波与平滑教程：十五年后》，《非线性滤波手册》12.656-704（2009）：3。

另请参阅

||

相关的话题

控制系统工具箱文档

金宝app

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

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

下载代码示例