利用扩展卡尔曼滤波算法和测量输出数据估计范德堡尔振荡器的状态。振荡器有两个状态和一个输出。

为振荡器创建一个扩展的卡尔曼滤波器对象。使用以前写的和保存的状态转换和测量功能，vdpStateFcn.m而且vdpMeasurementFcn.m．这些函数描述了具有非线性参数的范德堡尔振荡器的离散近似 $\mathrm{\mu }$等于1。该函数假定系统中存在附加的过程噪声和测量噪声。将两个状态的初始状态值指定为[1;0]。这是初始时刻状态值的猜测k，基于知识的系统输出直到时间k - 1， $$\underset{}{\overset{ˆ}{x}}［k|k-1］$$．

obj = extendedKalmanFilter(@vdpStateFcn，@vdpMeasurementFcn，[1;0]);

加载测量输出数据y从振荡器。在本例中，使用模拟静态数据进行说明。数据存储在vdp_data.mat文件。

负载vdp_data.maty

指定振荡器的过程噪声和测量噪声协方差。

obj。ProcessNoise = 0.01;obj。米easurementNoise = 0.16;

初始化数组以捕获估计的结果。

residBuf = [];xcorBuf = [];xpredBuf = [];

实现了扩展的卡尔曼滤波算法来估计振荡器的状态正确的而且预测命令。你首先纠正 $$\underset{}{\overset{ˆ}{x}}［k|k-1］$$及时使用测量k得到 $$\underset{}{\overset{ˆ}{x}}［k|k］$$．然后，预测下一个时间步的状态值 $$\underset{}{\overset{ˆ}{x}}［k+1|k］$$使用 $$\underset{}{\overset{ˆ}{x}}［k|k］$$，时间步长的状态估计k这是估计使用测量直到时间k．

要模拟实时数据测量，请每次一步使用测量数据。计算预测值和实际测量值之间的残差，以评估滤波器的性能和收敛程度。计算残差是一个可选步骤。当你使用剩余时，将命令放在正确的命令。如果预测值与测量值相匹配，则残差为零。

执行时间步骤的实时命令后，缓冲结果，以便在运行完成后绘制它们。

为k = 1:大小(y)[残差，残差方差]=残差(obj,y(k));[CorrectedState,CorrectedStateCovariance] = correct(obj,y(k));[PredictedState,PredictedStateCovariance] = predict(obj);residBuf(k，:) =残差;xcorBuf(k，:) = CorrectedState';xpredBuf(k，:) = PredictedState';结束

当你使用正确的命令,obj。状态而且obj。StateCovariance更新与修正的状态和状态估计误差协方差值的时间步长k，CorrectedState而且CorrectedStateCovariance．当你使用预测命令,obj。状态而且obj。StateCovariance是否用时间步长的预测值进行更新k + 1，PredictedState而且PredictedStateCovariance．当你使用剩余命令时，不进行任何修改obj属性。

在本例中，您使用了正确的之前预测因为初始状态值是 $$\underset{}{\overset{ˆ}{x}}［k|k-1］$$，初始时刻状态值的猜测k根据系统输出直到时间k - 1．如果初始状态值为 $$\underset{}{\overset{ˆ}{x}}［k-1|k-1］$$为前一时刻的值k - 1基于测量直到k - 1，然后使用预测命令。有关使用顺序的更多信息预测而且正确的,请参阅使用预测和正确的命令．

使用后校正值绘制估计状态。

plot(xcorBuf(:，1)， xcorBuf(:，2))“估计状态”）

画出实际测量值、修正后的估计测量值和残差。的测量函数vdpMeasurementFcn，测量为第一状态。

M = [y,xcorBuf(:，1)，residBuf];情节(M)网格在标题('实际和估计测量，剩余')传说(“测量”，“估计”，“残留”）

估算密切跟踪测量结果。在初始瞬态之后，残余量在整个运行过程中保持相对较小。

加载范德波尔ODE数据，并指定采样时间。

vdpODEdata.mat包含了具有初始条件的非线性参数mu=1的van der Pol ODE的模拟，使用ode45(2, 0)．真实状态随采样时间的变化而变化Dt = 0.05．

目录(fullfile (matlabroot,“例子”，“控制”，“主要”）)添加示例数据负载(“vdpODEdata.mat”，“xTrue”，“dt”) tSpan = 0:dt:5;

测量尺寸。对于本例，传感器测量带有标准偏差的高斯噪声的第一个状态0.04．

sqrtR = 0.04;yMeas = xTrue(:，1) + sqrtR*randn(nummel (tSpan)，1);

创建一个粒子滤波器，并设置状态转换和测量似然函数。

myPF = particleFilter(@vdpParticleFilterStateFcn，@ vdpmeasurementlikehoodfcn);

初始化粒子过滤器状态[2;0]用单位协方差，和使用1000粒子。

初始化(myPF, 1000,(2, 0),眼(2));

选择的意思是状态估计和系统的重采样方法。

myPF。StateEstimationMethod =“的意思是”；myPF。ResamplingMethod =“系统”；

估计状态使用正确的而且预测命令，并存储估计的状态。

xEst = 0 (size(xTrue));为k=1:size(xTrue,1) xEst(k，:) =正确的(myPF,yMeas(k));预测(myPF);结束

绘制结果，并比较估计状态和真实状态。

图(1)情节(xTrue (: 1), xTrue (:, 2),“x”x (: 1) x (:, 2),“罗”)传说(“真正的”，“估计”）

rmpath (fullfile (matlabroot,“例子”，“控制”，“主要”）)删除示例数据

考虑一个有输入的非线性系统u的状态x和测量y根据以下状态转换和测量方程演化:

过程噪声w系统的噪声是可加性的，而测量噪声是可加性的v非相加。

为系统创建状态转换函数和测量函数。使用附加输入指定函数u．

F = @(x,u)(根号(x+u));H = @(x,v,u)(x+2*u+v^2);

f而且h是匿名函数的函数句柄，分别存储状态转换和测量函数。在测量函数中，由于测量噪声是非加性的，v也指定为输入。请注意,v是否在附加输入之前指定为输入u．

创建一个扩展的卡尔曼滤波器对象，用于使用指定的函数估计非线性系统的状态。指定状态的初始值为1，测量噪声为非可加性。

obj = extendedKalmanFilter(f,h,1，“HasAdditiveMeasurementNoise”、假);

指定测量噪声协方差。

obj。米easurementNoise = 0.01;

方法可以估计系统的状态预测而且正确的命令。传递的值u来预测而且正确的，然后分别将它们传递给状态转换函数和测量函数。

用测量纠正状态估计y［k] =0.8和输入u［k] =0.2时间步长k．

正确的(obj, 0.8, 0.2)

预测给定的下一个时间步的状态u［k] =0.2．

预测(obj, 0.2)

检索错误，或者剩余在预测和测量之间。

[Residual, ResidualCovariance] = Residual (obj,0.8,0.2);