Идентификация параметра динамической системы с помощью расширенного фильтра Калмана

Введение

Оценка состояния динамических систем по зашумлённым измерениям является фундаментальной задачей в теории управления, навигации и обработке сигналов. Для линейных систем с гауссовскими шумами оптимальным решением является фильтр Калмана, однако многие реальные системы требуют одновременной оценки как переменных состояния, так и неизвестных параметров модели.

В таких случаях применяют расширенный фильтр Калмана (РФК), основанный на локальной линеаризации нелинейных функций модели через разложение в ряд Тейлора первого порядка. В данной работе рассматривается задача идентификации параметра и оценки состояния апериодического звена первого порядка:

,

где — неизвестный параметр динамики, — известный коэффициент, — входное воздействие, — шум процесса. Для решения используется расширение вектора состояния за счёт включения в него неизвестного параметра .

Цель работы — продемонстрировать реализацию РФК для совместного оценивания, исследовать сходимость оценки параметра и проанализировать поведение ковариационной матрицы ошибки.

Необходимые библиотеки

В данном примере нам понадобятся следующие библиотеки:

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)

0.9

Определим начальную оценку параметра 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

generate_data (generic function with 1 method)

Функция расширенного фильтра Калмана

Создадим функцию 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

ekf_filter (generic function with 1 method)

Запуск моделирования

Запустим моделирование: сгенерируем данные и применим фильтр.

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

200-element Vector{Float64}:
 0.14310081619978648
 0.23915330583528285
 0.20957303927390009
 0.26607767068000626
 0.3884568888337616
 0.5100903333114909
 0.5344283013016132
 0.5692036761074832
 0.6351143666385379
 0.6772130891330526
 0.7201782083990114
 0.700702597158715
 0.7237308724987693
 ⋮
 1.007264754012645
 1.011002756933732
 0.9996417060121864
 1.007551076756899
 1.0070630128722802
 1.011825792285365
 1.011063881429158
 0.9942177135507816
 0.991925762549218
 1.0136964386882734
 1.0224545667226046
 1.0082160811729803

Визуализация результатов

Построим три графика. На первом покажем истинное состояние (синяя линия), зашумленные измерения (оранжевые точки) и оценку фильтра (зеленая линия). На втором — истинное значение параметра (черный пунктир) и его оценку (оранжевая линия), которая начинается с и сходится к истинному значению . На третьем — сходимость неопределённости оценки параметра .

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/ГЛОНАСС.
Энергетика — оценка параметров электрических машин и состояния аккумуляторных батарей (заряд, ёмкость, износ).
Медицина и биотехнологии — оценка фармакокинетических параметров по данным измерений, обработка сигналов ЭЭГ и ЭКГ.
Метеорология — ассимиляция данных в численные модели погоды, где одновременно уточняются состояние атмосферы и параметры модели.

Таким образом, продемонстрированная схема совместного оценивания состояния и параметров на основе РФК является универсальным инструментом для работы с неполными или зашумлёнными данными о динамических системах.

Расширенный фильтр Калмана