Идентификация параметра динамической системы с помощью расширенного фильтра Калмана
Введение
Оценка состояния динамических систем по зашумлённым измерениям является фундаментальной задачей в теории управления, навигации и обработке сигналов. Для линейных систем с гауссовскими шумами оптимальным решением является фильтр Калмана, однако многие реальные системы требуют одновременной оценки как переменных состояния, так и неизвестных параметров модели.
В таких случаях применяют расширенный фильтр Калмана (РФК), основанный на локальной линеаризации нелинейных функций модели через разложение в ряд Тейлора первого порядка. В данной работе рассматривается задача идентификации параметра и оценки состояния апериодического звена первого порядка:
,
где — неизвестный параметр динамики, — известный коэффициент, — входное воздействие, — шум процесса. Для решения используется расширение вектора состояния за счёт включения в него неизвестного параметра .
Цель работы — продемонстрировать реализацию РФК для совместного оценивания, исследовать сходимость оценки параметра и проанализировать поведение ковариационной матрицы ошибки.
Необходимые библиотеки
В данном примере нам понадобятся следующие библиотеки:
- LinearAlgebra — для матричных операций;
- Statistics — для генерации шумов.
using LinearAlgebra, Statistics
Исходные параметры модели
Определим параметры системы:
T_sim = 200— количество шагов моделирования;k_u = 1.0— коэффициент усиления;T_const = 10.0— постоянная времени.
Вычислим истинные параметры: и , описывающие апериодическое звено первого порядка.
# Параметры модели
T_sim = 200;
k_u = 1.0;
T_const = 10.0;
a_true = 1.0 - 1.0/T_const;
b_true = k_u / T_const;
display(a_true)
Определим начальную оценку параметра a_init = 0.85 и входной сигнал u_d = 1.0. Также определим шумы через среднеквадратические отклонения: σ_v = √0.01 — шум процесса, σ_ζ = 0.001 — шум параметра (моделирует дрейф параметра), σ_w = √0.01 — шум измерений.
# Начальные условия
a_init = 0.85;
u_d = 1.0;
# Ковариации шумов
σ_v = sqrt(0.01);
σ_ζ = 0.001;
σ_w = sqrt(0.01);
Генерация данных
Создадим функцию generate_data(), которая имитирует поведение реальной системы. На каждом шаге генерируем шум процесса v и вычисляем следующее состояние по формуле:
Затем добавим шум измерений: .
Функция возвращает массивы истинных состояний и зашумленных измерений.
# Генерация истинного процесса
function generate_data(T_sim, σ_v, a_true, b_true, u_d, σ_w)
x_true = zeros(T_sim)
y_meas = zeros(T_sim)
x = 0.0
for n in 1:T_sim
v = σ_v * randn()
x = a_true * x + b_true * u_d + b_true * v
x_true[n] = x
y_meas[n] = x + σ_w * randn()
end
return x_true, y_meas
end
Функция расширенного фильтра Калмана
Создадим функцию ekf_filter(), реализующую расширенный фильтр Калмана для совместной оценки состояния x и параметра a. Инициализируем расширенный вектор состояния X = [x; a] = [0; 0.85] и ковариационную матрицу P = diag(1.0, 0.1).
На каждом шаге цикла:
-
Формирование матриц: — матрица перехода, — матрица входа, — матрица шума.
-
Прогноз состояния:
-
Линеаризация: вычисление матрицы Якоби , где учтена нелинейность из-за произведения .
-
Прогноз ковариации: , где .
-
Коррекция: вычисление расхождения ,
инновационной ковариации , коэффициента Калмана .
-
Обновление: ,
.
function ekf_filter(y_meas::Vector{Float64}, a_init::Float64, b_true::Float64, u_d::Float64, σ_v::Float64, σ_ζ::Float64, σ_w::Float64)
T_sim = length(y_meas)
# Расширенный вектор состояния: [x; a]
X = [0.0; a_init]
# Ковариационная матрица ошибки
P = diagm([1.0, 0.1])
# Матрица измерения
C = [1.0 0.0]
# Ковариационная матрица шума процесса
M = diagm([σ_v^2, σ_ζ^2])
# Единичная матрица 2x2
I2 = diagm(ones(2))
# Хранение результатов
x_est = zeros(T_sim)
a_est = zeros(T_sim)
P_hist = zeros(2, 2, T_sim)
for n in 1:T_sim
x_n = X[1]
a_n = X[2]
# Прогноз состояния
# Матрица перехода A = [a 0; 0 1]
A = [a_n 0.0;
0.0 1.0]
# Матрица управления B = [b; 0] (b - константа)
B = [b_true; 0.0]
X_pred = A * X + B * u_d
# Линеаризация
F = [a_n x_n;
0.0 1.0]
# Матрица шума процесса V = [b 0; 0 1]
V_mat = [b_true 0.0;
0.0 1.0]
# Прогноз ковариации
P_pred = F * P * F' + V_mat * M * V_mat'
# Коррекция
y_pred = (C * X_pred)[1]
innovation = y_meas[n] - y_pred
S = (C * P_pred * C')[1, 1] + σ_w^2
K = P_pred * C' / S
X = X_pred + K * innovation
P = (I2 - K * C) * P_pred
# Сохранение
x_est[n] = X[1]
a_est[n] = X[2]
P_hist[:, :, n] = P
end
return x_est, a_est, P_hist
end
Запуск моделирования
Запустим моделирование: сгенерируем данные и применим фильтр.
x_true, y_meas = generate_data(T_sim, σ_v, a_true, b_true, u_d, σ_w);
x_est, a_est, P_hist = ekf_filter(y_meas, a_init, b_true, u_d, σ_v, σ_ζ, σ_w);
Для визуального контроля, выведем первые и последние значения вектора оценки расширенного фильтра Калмана.
x_est
Визуализация результатов
Построим три графика. На первом покажем истинное состояние (синяя линия), зашумленные измерения (оранжевые точки) и оценку фильтра (зеленая линия). На втором — истинное значение параметра (черный пунктир) и его оценку (оранжевая линия), которая начинается с и сходится к истинному значению . На третьем — сходимость неопределённости оценки параметра .
p1 = plot(1:T_sim, x_true, label="Истинное состояние", lw=2)
plot!(p1, 1:T_sim, y_meas, label="Измерения", alpha=0.4, marker=(:circle, 2))
plot!(p1, 1:T_sim, x_est, label="Оценка РФК", lw=2)
xlabel!(""); ylabel!("x")
title!("Оценка состояния")
p2 = plot(1:T_sim, a_true*ones(T_sim), label="Истинный параметр a", lw=2, color=:black, ls=:dash)
plot!(p2, 1:T_sim, a_est, label="Оценка параметра a", lw=2)
xlabel!(""); ylabel!("a")
title!("Оценка параметра")
P_a = [P_hist[2, 2, n] for n in 1:T_sim]
p3 = plot(1:T_sim, P_a,
label="Дисперсия параметра a",
lw=2,
color=:purple,
fillrange=0,
fillalpha=0.2,
legend=:topright)
xlabel!("Шаг")
ylabel!("Дисперсия")
title!("Сходимость неопределённости оценки параметра a")
plot(p1, p2, p3, layout=(3,1), size=(800, 750))
Заключение
В ходе работы реализован расширенный фильтр Калмана для совместной оценки состояния и неизвестного параметра линейного динамического объекта. Моделирование показало:
-
Сходимость оценки параметра. Несмотря на намеренное отклонение начальной оценки, фильтр быстро скорректировал её, и после первых десятков итераций оценка стабилизировалась в окрестности истинного значения.
-
Сходимость неопределённости. Дисперсия ошибки оценивания параметра монотонно убывала, что свидетельствует о накоплении информации и повышении уверенности алгоритма в оценке.
-
Качество оценки состояния. Восстановленная траектория хорошо согласуется с истинной, несмотря на шумы процесса и измерений.
Важно подчеркнуть, что корректная работа РФК напрямую зависит от согласованности модели с реальным процессом, при несоответствии оценка параметра смещается.
Рассмотренный подход широко применяется в реальных задачах:
-
Адаптивное управление — оценка параметров объекта (постоянные времени, коэффициенты усиления) в реальном времени для построения адаптивных регуляторов.
-
Навигация и локализация — совместное уточнение координат, ориентации и калибровочных параметров датчиков в инерциальных системах, GPS/ГЛОНАСС.
-
Энергетика — оценка параметров электрических машин и состояния аккумуляторных батарей (заряд, ёмкость, износ).
-
Медицина и биотехнологии — оценка фармакокинетических параметров по данным измерений, обработка сигналов ЭЭГ и ЭКГ.
-
Метеорология — ассимиляция данных в численные модели погоды, где одновременно уточняются состояние атмосферы и параметры модели.
Таким образом, продемонстрированная схема совместного оценивания состояния и параметров на основе РФК является универсальным инструментом для работы с неполными или зашумлёнными данными о динамических системах.