主要内容

本页采用了机器翻译。点击此处可查看英文原文。

generateJacobianFcn

使用自动微分生成扩展卡尔曼滤波器的 MATLAB 雅可比函数

自 R2023a 起

说明

fcnStateJac = generateJacobianFcn(obj,'state',Us1,...,Usn) 利用自动微分技术,为扩展卡尔曼滤波器 (EKF) 生成状态转移雅可比矩阵函数。

该函数会在当前文件夹中生成两个 MATLAB® 函数文件:

  • stateTransitionJacobianFcn.m - 生成的状态转移雅可比函数

  • stateTransitionJacobianFcnAD.m - 一个利用自动微分生成状态转移雅可比矩阵的辅助函数

fcnStateJac 是一个匿名函数的句柄,该函数使用输入对象 extendedKalmanFilter (obj)、传递给 objpredict 函数的附加参量 (Us1,...,Usn),以及内部用于计算状态转移雅可比矩阵的常量来调用 stateTransitionJacobianFcn.m

要在 EKF 对象中使用此雅可比矩阵,请在该对象的 StateTransitionJacobianFcn 属性中指定 fcnStateJac。例如:

obj.StateTransitionJacobianFcn = fcnStateJac;

示例

fcnMeasurementJac = generateJacobianFcn(obj,'measurement',Um1,...,Umn) 利用自动微分技术,为扩展卡尔曼滤波器 (EKF) 生成测量雅可比矩阵函数。

该函数会在当前文件夹中生成两个 MATLAB 函数文件:

  • measurementJacobianFcn.m - 生成的测量雅可比函数

  • measurementJacobianFcnAD.m - 一个利用自动微分生成测量雅可比矩阵的辅助函数

fcnMeasurementJac 是一个匿名函数的句柄,该函数使用输入对象 extendedKalmanFilter (obj)、传递给 objcorrect 函数的附加参量 (Um1,...,Umn),以及内部用于计算测量雅可比矩阵的常量来调用 measurementJacobianFcn.m

要在 EKF 对象中使用此雅可比矩阵,请在该对象的 MeasurementJacobianFcn 属性中指定 fcnMeasurementJac。例如:

obj.MeasurementJacobianFcn = fcnMeasurementJac;

示例

[___,constants] = generateJacobianFcn(___) 还会返回用于计算雅可比函数的常量。您可以返回任意一个雅可比函数的常数。

[___] = generateJacobianFcn(___,FileName=filename) 指定了生成的雅可比矩阵函数文件的名称以及生成这些文件的文件夹位置。

示例

示例

全部折叠

为一个具有两个状态和一个输出的范德波尔振荡器创建一个扩展卡尔曼滤波器 (EKF) 对象。使用之前编写并保存的状态转移和测量函数 vdpStateFcn.mvdpMeasurementFcn.m。将这两个状态的初始状态值指定为 [2;0]

obj = extendedKalmanFilter(@vdpStateFcn,@vdpMeasurementFcn,[2;0]);

扩展卡尔曼滤波算法利用状态转移函数和测量函数的雅可比矩阵进行状态估计。您可以编写并保存雅可比函数,并将它们作为函数句柄提供给 EKF 对象。对于此对象,请使用之前编写并保存的函数 vdpStateJacobianFcn.mvdpMeasurementJacobianFcn.m

obj.StateTransitionJacobianFcn = @vdpStateJacobianFcn;
obj.MeasurementJacobianFcn = @vdpMeasurementJacobianFcn;

如果无法获得雅可比函数,可以使用自动微分来生成它们。

obj.StateTransitionJacobianFcn = generateJacobianFcn(obj,'state');
obj.MeasurementJacobianFcn = generateJacobianFcn(obj,'measurement');

如果您未指定函数的雅可比矩阵,软件将对其进行数值计算。这种数值计算可能会导致处理时间增加,并造成状态估计的数值不准确。

考虑一个输入为 u 的非线性系统,其状态 x 和测量值 y 分别按照以下状态转移方程和测量方程演化:

x[k]=x[k-1]+u[k-1]+w[k-1]

y[k]=x[k]+2*u[k]+v[k]2

系统的过程噪声 w 具有加性,而测量噪声 v 则不具有加性。

为该系统创建状态转移函数和测量函数。请使用附加输入 u 来指定函数。

f = @(x,u)(sqrt(x+u));
h = @(x,v,u)(x+2*u+v^2);

fh 分别是存储状态转移函数和测量函数的匿名函数的句柄。在测量函数中,由于测量噪声是非加性的,因此还将 v 指定为输入。请注意,v 被指定为在额外输入 u 之前的输入。

创建一个扩展卡尔曼滤波器 (EKF) 对象,用于利用状态转移函数和测量函数估计非线性系统的状态。将状态的初始值指定为 1,并将测量噪声指定为非加性。

obj = extendedKalmanFilter(f,h,1,"HasAdditiveMeasurementNoise",false);

使用自动微分生成 obj 的状态转移和测量雅可比函数。由于状态转移函数和测量函数都带有额外的输入 u,请向 generateJacobianFcn 传递一个与 u 类型和大小相同的任意值。

obj.StateTransitionJacobianFcn = generateJacobianFcn(obj,"state",0.2);
obj.MeasurementJacobianFcn = generateJacobianFcn(obj,"measurement",0.2);

使用自动微分技术生成 EKF 对象的状态转移和测量雅可比函数。将雅可比函数文件保存到非默认位置。

为一个具有两个状态和一个输出的范德波尔振荡器创建一个扩展卡尔曼滤波器 (EKF) 对象。使用之前编写并保存的状态转移和测量函数 vdpStateFcn.mvdpMeasurementFcn.m。将这两个状态的初始状态值指定为 [2;0]

obj = extendedKalmanFilter(@vdpStateFcn,@vdpMeasurementFcn,[2;0]);

创建一个文件夹,用于生成雅可比函数文件。将该文件夹添加到路径中。

% Change to a folder where you have write permission
CWD = pwd;
c = onCleanup(@()cd(CWD));
cd(tempdir)

% create a folder of the desired name
folder = 'jacobians';
[status,msg,msgID] = mkdir(folder);
% add the new folder to MATLAB path
addpath(folder)

将雅可比函数文件生成到新文件夹中。显示生成的文件。

fcnStateJac = generateJacobianFcn(obj,"state",...
    FileName=fullfile(folder,'sJac'));
fcnMeasurementJac = generateJacobianFcn(obj,"measurement",...
    FileName=fullfile(folder,'mJac'));

{dir(fullfile(folder,"*.m")).name}'
ans = 4×1 cell
    {'mJac.m'  }
    {'mJacAD.m'}
    {'sJac.m'  }
    {'sJacAD.m'}

将生成的函数句柄添加到 EKF 对象中。

obj.StateTransitionJacobianFcn = @fcnStateJac;
obj.MeasurementJacobianFcn = @fcnMeasurementJac;
delete(c); % restore working folder

输入参数

全部折叠

扩展卡尔曼滤波器,指定为一个 extendedKalmanFilter 对象。

状态转移函数的其他参量。状态转移函数 fobjStateTransitionFcn 属性中进行了定义。请指定与 obj 函数中的 predict 函数所使用的相同的附加参量。输入参量可以是任何类型。

测量函数的其他参量。测量函数 hobjMeasurementFcn 属性中进行了定义。请指定与 obj 函数中的 correct 函数所使用的相同的附加参量。输入参量可以是任何类型。

雅可比函数的文件名,指定为字符向量。该函数将此命名约定应用于生成的雅可比函数文件:

  • filename.m

  • filenameAD.m

如果 filename 包含绝对或相对文件夹路径,则 generateJacobianFcn 会将文件生成到该文件夹中。如果 filename 无法找到指定的路径,则 generateJacobianFcn 会返回一个错误。

如果 generateJacobianFcn 生成文件的文件夹不在 MATLAB 的路径上,则 generateJacobianFcn 将雅可比函数句柄作为 [] 返回。如果出现此问题,请按照以下步骤生成函数句柄:

注意

要使用此过程,您必须在调用 generateJacobianFcn 时返回 constants 参量。

  1. 将该文件夹添加到 MATLAB 搜索路径中,或者将生成的函数文件移动到 MATLAB 搜索路径中的某个文件夹内。

  2. 获取 filename 的名称部分,用于函数句柄名称 FcnName

    [~,FcnName] = fileparts(filename);
  3. 使用此名称生成函数句柄,其中 constants 是您之前对 generateJacobianFcn 调用的第二个参量。

    fcnJac = @(varargin)FcnName(varargin{:},constants);

然后,您可以在 obj 的相应雅可比函数属性中指定 fcnJac 的句柄。

有关使用 MATLAB 搜索路径的详细信息,请参阅什么是 MATLAB 搜索路径?

示例: FileName='measjac' 会在当前文件夹中生成 measjac.mmeasjacAD.m 文件。

示例: FileName='AD/measjac' 会在当前文件夹的 AD 子文件夹中生成 measjac.mmeasjacAD.m 文件。

示例: FileName='C:/AD/measjac'C:/AD 文件夹中生成 ekfjac.mmeasjacAD.m 文件。

数据类型: char

输出参量

全部折叠

用于计算雅可比矩阵的其他常量,以长度为 N 的常量值元胞数组形式返回,其中 N 为常量的个数。如果雅可比函数不需要常数,则 constants 是一个空元胞数组。

如果指定了 filename,而 filename 中指定的文件夹不在 MATLAB 的路径中,则可以使用这些常量手动构建函数句柄。

状态转移函数的雅可比矩阵,以函数句柄的形式返回。

  • 如果 obj.HasAdditiveProcessNoise = true,那么 fcnStateJac 的签名如下:

    dx = fcnStateJac(obj,Us1,...,Usn,constants)
  • 如果 obj.HasAdditiveProcessNoise = false,那么 fcnStateJac 的签名如下:

    [dx,dw] = fcnStateJac(obj,w,Us1,...,Usn,constants)

在以下函数签名中:

  • dx 是预测状态相对于前一状态的雅可比矩阵。

  • dw 是预测状态相对于过程噪声元素的雅可比矩阵。

  • w 是过程噪声变量。

  • Us1,...,Usnobjpredict 函数所使用的附加参量。

  • constants 是用于计算雅可比矩阵的附加常数。

测量函数的雅可比矩阵,以函数句柄的形式返回。

  • 如果 obj.HasAdditiveMeasurementNoise = true,那么 fcnMeasurementJac 的签名如下:

    dy = fcnMeasurementJac(obj,Um1,...,Umn,constants)
  • 如果 obj.HasAdditiveMeasurementNoise = false,那么 fcnMeasurementJac 的签名如下:

    [dy,dv] = fcnMeasurementJac(obj,v,Um1,...,Umn,constants)

在以下函数签名中:

  • dy 是测量函数相对于态的雅可比矩阵。

  • dv 是测量函数相对于测量噪声的雅可比矩阵。

  • v 是测量噪声变量。

  • Um1,...,Umnobjcorrect 函数所使用的附加参量。

  • constants 是用于计算雅可比矩阵的附加常数。

限制

  • 目前,自动微分仅支持有限的数学运算集,相关内容详见优化变量和表达式支持的运算 (Optimization Toolbox)。如果原始的状态转移或测量函数使用了列表中未包含的操作或函数,或者包含 if-else 语句或循环,则 generateJacobianFcn 将因错误而终止。

  • 要生成雅可比函数,请勿在原始函数中预分配任何优化变量。例如,假设您试图从包含以下代码的函数中生成雅可比矩阵。

    dxdt = zeros(2,1); 
    dxdt(1) = x(1)*x(2);
    dxdt(2) = x(1)/x(2); 
    这段代码会导致以下错误。
    Unable to perform assignment because value of type 
    'optim.problemdef.OptimizationExpression' 
    is not convertible to 'double'.
    请改用这段代码。
    dxdt = [x(1)*x(2); x(1)/x(2)];

  • 建议将 obj 中的状态转移和测量函数指定为当前文件夹或 MATLAB 路径下某个文件夹中的文件。虽然在生成雅可比矩阵函数时支持本地函数的句柄,但在生成 C/C++ 部署代码时则不支持。有关局部函数的信息,请参阅局部函数

版本历史记录

在 R2023a 中推出