使用卡尔曼滤波方法对非线性系统进行参数辨识
使用卡尔曼滤波方法（特别是扩展卡尔曼滤波EKF）对非线性系统进行参数辨识
1. 理论基础：EKF用于参数辨识
核心思想
将系统参数视为状态变量的一部分，通过状态扩展方法进行联合估计：
扩展状态向量：x = [系统状态; 系统参数]
关键步骤：
- 状态扩展：将待辨识参数作为额外的状态变量
- 预测步骤：使用系统模型预测状态
- 更新步骤：使用观测值校正状态估计
- 参数更新：参数作为状态的一部分被同时估计
2. MATLAB实现：完整示例
示例系统：带未知参数的Van der Pol振荡器
Van der Pol方程：ẍ - μ(1-x²)ẋ + ω²x = 0
待辨识参数：μ 和 ω
%% EKF参数辨识：Van der Pol振荡器
clear; close all; clc;
%% 1. 系统参数设置
dt = 0.01; % 采样时间
T = 10; % 总时间
t = 0:dt:T; % 时间向量
N = length(t); % 数据点数
% 真实参数（待辨识）
mu_true = 1.5; % 非线性阻尼系数
omega_true = 3; % 自然频率
% 噪声设置
Q = diag([0.01, 0.01, 0.001, 0.001]); % 过程噪声协方差
R = 0.1; % 测量噪声协方差
%% 2. 生成模拟数据（真实系统）
fprintf('生成模拟数据...\n');
% 定义Van der Pol方程
vdp_ode = @(t, x) [x(2);
mu_true*(1-x(1)^2)*x(2) - omega_true^2*x(1)];
% 初始状态
x0 = [1; 0]; % [位置; 速度]
x_true = zeros(2, N);
x_true(:, 1) = x0;
% 数值积分（使用ode45）
[~, x_ode] = ode45(vdp_ode, t, x0);
x_true = x_ode';
% 添加测量噪声
y_meas = x_true(1, :) + sqrt(R)*randn(1, N);
%% 3. 扩展卡尔曼滤波（EKF）实现
fprintf('开始EKF参数辨识...\n');
% 扩展状态向量：[x1; x2; mu; omega]
x_est = zeros(4, N); % 状态估计
P_est = zeros(4, 4, N); % 误差协方差
% 初始猜测
x_est(:, 1) = [x0; 1.0; 2.5]; % 参数初始猜测
P_est(:, :, 1) = diag([0.5, 0.5, 0.5, 0.5]);
% EKF主循环
for k = 1:N-1
% ---------- 预测步骤 ----------
% 状态预测（使用欧拉法）
x1 = x_est(1, k);
x2 = x_est(2, k);
mu_est = x_est(3, k);
omega_est = x_est(4, k);
% 状态转移
x_pred = [x1 + dt*x2;
x2 + dt*(mu_est*(1-x1^2)*x2 - omega_est^2*x1);
mu_est;
omega_est];
% 计算雅可比矩阵（状态转移的偏导数）
F = zeros(4, 4);
F(1, 2) = 1;
F(2, 1) = -dt*(2*mu_est*x1*x2 + omega_est^2);
F(2, 2) = 1 + dt*mu_est*(1-x1^2);
F(2, 3) = dt*(1-x1^2)*x2;
F(2, 4) = -2*dt*omega_est*x1;
F(3, 3) = 1;
F(4, 4) = 1;
% 协方差预测
P_pred = F * P_est(:, :, k) * F' + Q;
% ---------- 更新步骤 ----------
% 观测矩阵（只能观测到位置x1）
H = [1, 0, 0, 0];
% 计算卡尔曼增益
S = H * P_pred * H' + R;
K = P_pred * H' / S;
% 状态更新
y_pred = H * x_pred;
x_est(:, k+1) = x_pred + K * (y_meas(k+1) - y_pred);
% 协方差更新
P_est(:, :, k+1) = (eye(4) - K * H) * P_pred;
% 确保参数为正（物理约束）
x_est(3, k+1) = max(x_est(3, k+1), 0.1); % mu > 0.1
x_est(4, k+1) = max(x_est(4, k+1), 0.1); % omega > 0.1
end
%% 4. 结果可视化
figure('Position', [100, 100, 1200, 800]);
% 子图1：状态估计对比
subplot(2, 3, [1, 2]);
plot(t, x_true(1, :), 'b-', 'LineWidth', 2, 'DisplayName', '真实位置');
hold on;
plot(t, x_est(1, :), 'r--', 'LineWidth', 1.5, 'DisplayName', '估计位置');
plot(t, y_meas, 'g.', 'MarkerSize', 8, 'DisplayName', '测量值');
xlabel('时间 (s)'); ylabel('位置 x_1');
title('状态估计对比');
legend('Location', 'best');
grid on;
% 子图2：参数估计过程
subplot(2, 3, 3);
plot(t, mu_true*ones(1, N), 'b-', 'LineWidth', 2, 'DisplayName', '真实值');
hold on;
plot(t, x_est(3, :), 'r--', 'LineWidth', 1.5, 'DisplayName', '估计值');
xlabel('时间 (s)'); ylabel('\mu');
title('参数 \mu 估计');
legend('Location', 'best');
grid on;
subplot(2, 3, 4);
plot(t, omega_true*ones(1, N), 'b-', 'LineWidth', 2, 'DisplayName', '真实值');
hold on;
plot(t, x_est(4, :), 'r--', 'LineWidth', 1.5, 'DisplayName', '估计值');
xlabel('时间 (s)'); ylabel('\omega');
title('参数 \omega 估计');
legend('Location', 'best');
grid on;
% 子图3：参数估计误差
subplot(2, 3, 5);
mu_error = abs(x_est(3, :) - mu_true);
omega_error = abs(x_est(4, :) - omega_true);
semilogy(t, mu_error, 'r-', 'LineWidth', 1.5, 'DisplayName', '\mu 误差');
hold on;
semilogy(t, omega_error, 'b-', 'LineWidth', 1.5, 'DisplayName', '\omega 误差');
xlabel('时间 (s)'); ylabel('绝对误差');
title('参数估计误差（对数坐标）');
legend('Location', 'best');
grid on;
% 子图4：收敛分析
subplot(2, 3, 6);
plot(t, P_est(3, 3, :), 'r-', 'LineWidth', 1.5, 'DisplayName', 'P_{\mu\mu}');
hold on;
plot(t, P_est(4, 4, :), 'b-', 'LineWidth', 1.5, 'DisplayName', 'P_{\omega\omega}');
xlabel('时间 (s)'); ylabel('协方差');
title('参数估计协方差');
legend('Location', 'best');
grid on;
%% 5. 性能评估
fprintf('\n========== 参数辨识结果 ==========\n');
fprintf('真实参数: μ = %.3f, ω = %.3f\n', mu_true, omega_true);
fprintf('最终估计: μ = %.3f, ω = %.3f\n', x_est(3, end), x_est(4, end));
fprintf('相对误差: μ = %.2f%%, ω = %.2f%%\n', ...
abs(x_est(3, end)-mu_true)/mu_true*100, ...
abs(x_est(4, end)-omega_true)/omega_true*100);
% 计算RMSE
mu_rmse = sqrt(mean((x_est(3, 100:end) - mu_true).^2));
omega_rmse = sqrt(mean((x_est(4, 100:end) - omega_true).^2));
fprintf('RMSE (去除初始瞬态): μ = %.4f, ω = %.4f\n', mu_rmse, omega_rmse);
3. 通用EKF参数辨识函数
下面提供一个可重用的EKF参数辨识函数：
function [x_est, P_est, params_est] = ekf_parameter_estimation(...
y_meas, dt, params_init, Q, R, sys_func, jac_func)
% EKF参数辨识通用函数
% 输入:
% y_meas: 测量数据
% dt: 采样时间
% params_init: 参数初始猜测
% Q: 过程噪声协方差
% R: 测量噪声协方差
% sys_func: 系统函数句柄 @(x, params)
% jac_func: 雅可比函数句柄 @(x, params)
% 输出:
% x_est: 状态估计
% P_est: 协方差估计
% params_est: 参数估计
N = length(y_meas);
n_states = length(params_init) + 2; % 假设2个系统状态
% 初始化
x_est = zeros(n_states, N);
P_est = zeros(n_states, n_states, N);
% 初始状态猜测
x_est(:, 1) = [0.1; 0; params_init']; % [x1; x2; params]
P_est(:, :, 1) = eye(n_states);
for k = 1:N-1
% 当前估计
x_current = x_est(:, k);
params_current = x_current(3:end);
% 预测步骤
x_pred = x_current;
x_pred(1:2) = sys_func(x_current(1:2), params_current, dt);
% 计算雅可比矩阵
F = jac_func(x_current(1:2), params_current, dt);
% 协方差预测
P_pred = F * P_est(:, :, k) * F' + Q;
% 更新步骤
H = [1, 0, zeros(1, length(params_current))]; % 观测位置
y_pred = H * x_pred;
% 卡尔曼增益
S = H * P_pred * H' + R;
K = P_pred * H' / S;
% 状态更新
x_est(:, k+1) = x_pred + K * (y_meas(k+1) - y_pred);
P_est(:, :, k+1) = (eye(n_states) - K * H) * P_pred;
end
params_est = x_est(3:end, :);
end
4. 应用示例：电机系统参数辨识
%% 示例：直流电机参数辨识
% 系统模型: J*θ̈ + B*θ̇ = K*V
% 待辨识参数: J(转动惯量), B(阻尼系数), K(转矩常数)
clear; close all; clc;
% 仿真参数
dt = 0.001;
T = 5;
t = 0:dt:T;
N = length(t);
% 真实参数
J_true = 0.01; % kg·m²
B_true = 0.1; % N·m·s/rad
K_true = 0.5; % N·m/V
% 输入信号（PWM电压）
V = 12 * (square(2*pi*2*t) + 1)/2; % 12V方波
% 模拟真实系统
theta = zeros(1, N);
omega = zeros(1, N);
for k = 2:N
alpha = (K_true*V(k) - B_true*omega(k-1)) / J_true;
omega(k) = omega(k-1) + alpha*dt;
theta(k) = theta(k-1) + omega(k)*dt;
end
% 添加噪声
y_meas = theta + 0.001*randn(1, N);
% 使用UKF进行参数辨识（更适合强非线性系统）
fprintf('使用UKF进行电机参数辨识...\n');
% UKF参数设置
alpha = 1e-3;
beta = 2;
kappa = 0;
% 状态维度：3个参数 + 2个状态
n = 5;
% UKF初始化
x_est = zeros(n, N);
x_est(:, 1) = [theta(1); omega(1); 0.015; 0.15; 0.6]; % 初始猜测
P = diag([0.1, 0.1, 0.005, 0.05, 0.1]); % 初始协方差
Q = diag([1e-6, 1e-6, 1e-6, 1e-6, 1e-6]); % 过程噪声
R = 1e-4; % 测量噪声
% UKF主循环
for k = 1:N-1
% 生成Sigma点
[chi, Wm, Wc] = ut(x_est(:, k), P, alpha, beta, kappa);
% 预测步骤
chi_pred = zeros(n, 2*n+1);
for i = 1:2*n+1
J = chi(3, i);
B = chi(4, i);
Kp = chi(5, i);
% 状态转移
alpha_m = (Kp*V(k) - B*chi(2, i)) / J;
chi_pred(1, i) = chi(1, i) + chi(2, i)*dt;
chi_pred(2, i) = chi(2, i) + alpha_m*dt;
chi_pred(3:5, i) = chi(3:5, i); % 参数保持不变
end
% 预测均值和协方差
x_pred = zeros(n, 1);
for i = 1:2*n+1
x_pred = x_pred + Wm(i) * chi_pred(:, i);
end
P_pred = zeros(n, n);
for i = 1:2*n+1
P_pred = P_pred + Wc(i) * (chi_pred(:, i) - x_pred) * (chi_pred(:, i) - x_pred)';
end
P_pred = P_pred + Q;
% 更新步骤
[chi_up, ~, ~] = ut(x_pred, P_pred, alpha, beta, kappa);
% 观测预测
z_pred = zeros(1, 2*n+1);
for i = 1:2*n+1
z_pred(i) = chi_up(1, i); % 观测位置
end
% 观测均值和协方差
z_mean = sum(Wm .* z_pred);
P_zz = 0;
for i = 1:2*n+1
P_zz = P_zz + Wc(i) * (z_pred(i) - z_mean) * (z_pred(i) - z_mean)';
end
P_zz = P_zz + R;
% 互协方差
P_xz = zeros(n, 1);
for i = 1:2*n+1
P_xz = P_xz + Wc(i) * (chi_up(:, i) - x_pred) * (z_pred(i) - z_mean)';
end
% 卡尔曼增益
K = P_xz / P_zz;
% 状态更新
x_est(:, k+1) = x_pred + K * (y_meas(k+1) - z_mean);
P = P_pred - K * P_zz * K';
end
% 结果显示
fprintf('\n========== 电机参数辨识结果 ==========\n');
fprintf('参数 真实值 估计值 误差\n');
fprintf('J(惯量) %.4f %.4f %.2f%%\n', ...
J_true, x_est(3, end), abs(x_est(3, end)-J_true)/J_true*100);
fprintf('B(阻尼) %.4f %.4f %.2f%%\n', ...
B_true, x_est(4, end), abs(x_est(4, end)-B_true)/B_true*100);
fprintf('K(转矩) %.4f %.4f %.2f%%\n', ...
K_true, x_est(5, end), abs(x_est(5, end)-K_true)/K_true*100);
参考代码 使用卡尔曼滤波方法对非线性系统进行参数辨识 www.youwenfan.com/contentcnp/97380.html
5. 实用建议与技巧
1. 算法选择指南
| 方法 | 适用场景 | 优缺点 |
|---|---|---|
| 扩展卡尔曼滤波(EKF) | 中度非线性系统，计算资源有限 | 实现简单，计算快，但线性化可能引入误差 |
| 无迹卡尔曼滤波(UKF) | 强非线性系统，需要更高精度 | 精度高，无需雅可比，但计算量较大 |
| 粒子滤波(PF) | 高度非线性、非高斯系统 | 最通用，但计算量大，可能发散 |
2. 参数辨识的关键要点
- 初始值选择：好的初始猜测能加速收敛
- 噪声协方差调整：需要根据实际情况调整Q和R
- 可观性分析：确保系统参数是可辨识的
- 数据预处理：滤波、去趋势等预处理能提高辨识精度
3. 常见问题与解决
% 问题1：滤波器发散
% 解决方案：增加过程噪声Q或使用自适应EKF
% 自适应噪声估计
if k > 10
innovation = y_meas(k) - y_pred;
R_adapt = 0.95*R_adapt + 0.05*innovation^2;
end
% 问题2：参数收敛慢
% 解决方案：增加参数的过程噪声（允许更快变化）
Q_params = diag([1e-4, 1e-4]); % 参数噪声比状态噪声大
6. 扩展功能：自适应EKF
% 自适应EKF实现框架
function [x_est, P_est, R_adapt] = adaptive_ekf(...)
% 初始化自适应噪声估计
R_adapt = R_init;
N_window = 10; % 自适应窗口大小
for k = 1:N-1
% ... EKF预测步骤 ...
% 自适应测量噪声估计
innovation = y_meas(k+1) - y_pred;
innovation_buffer(mod(k, N_window)+1) = innovation;
if k >= N_window
R_adapt = cov(innovation_buffer) + 1e-6*eye(size(R));
end
% 使用自适应R进行更新
S = H * P_pred * H' + R_adapt;
K = P_pred * H' / S;
% ... 状态更新 ...
end
end
浙公网安备 33010602011771号