Документация Engee
Notebook

Идентификация параметров модели в Engee

Открыть пример в Engee

Описание примера

Рассмотрим реализацию задачи параметрической идентификации модели на основании данных. Нам предстоит определить несколько неизменных во времени параметров некоторой модели, структура которой нам известна заранее.

В этом примере мы будем анализировать сигнал, порожденный динамической моделью Engee. Она находится в папке start/examples/data_analysis/parametric_identification.

Модель состоит из двух подсистем: одна порождает анализируемый сигнал, другая – эталонная модель (без шума).

image.png

В подсистеме ModelSystem, при помощи суммирования, объединяются пять сигналов:

  • постоянная компонента: число 5,
  • временная компонента: линия с наклоном 0.01,
  • первая гармоника: синусоида с амплитудой 1 и угловой частотой $2\pi \cdot 0.02$ рад/с,
  • менее заметная вторая гармоника: синусоида с амплитудой 0.6 и угловой частотой $2\pi \cdot 0.008$ рад/с,
  • белый шум с нулевым смещением и дисперсией 0.3.

image.png

В подсистеме ModelBlock мы сможем подставить идентифицированные параметры. Пока в ней заданы эталонные параметры и она генерирует эталонный сигнал. В ней модель описана в виде кода:

function (c::ModelBlock)(t::Real)
    T0 = 5
    T1 = 0.01
    A1, f1, phi1 = 1, 0.02, 0
    A2, f2, phi2 = 0.6, 0.008, 0
    
    signal = T0 + T1 * t +
             A1 * sin( 2*pi * f1 * t + phi1 ) .+
             A2 * sin( 2*pi * f2 * t + phi2 );
    return signal
end

Сохраним эталонные значения параметров:

In [ ]: 
# "Идеальные" параметры линейной части сигнала (смещение, наклон)
T0m = 5;
T1m = 0.01;
# "Идеальные" параметры первой синусоиды (амплитуда, частота, фаза)
A1m = 1;
f1m = 0.02;
phi1m = 0;
# "Идеальные" параметры второй синусоиды (амплитуда, частота, фаза)
A2m = 0.6;
f2m = 0.008;
phi2m = 0;

В результате запуска модели мы получаем 500 секунд сигнала. Вот график последних 200 секунд:

image_2.png

Блок ModelBlock сейчас инициализируется точно теми же параметрами, что и динамическая модель в подсистеме ModelSystem, поэтому график from_code (темно-синего цвета) полностью повторяет систематическую часть сигнала from_model (оранжевого цвета).

Какие данные сохраняет модель

Оба сигнала сохраняются в CSV таблице param_id_data.csv со следующими столбцами:

  • time – столбец времени,
  • m – столбец данных, порожденных динамической моделью (эталонные параметры и погрешность),
  • c – столбец данных, порожденных эталонной моделью (выполненной в виде кода).

Порядок действий в этом примере

Нам предстоит выполнить следующие действия:

  1. запустить динамическую модель чтобы синтезировать модельные данные
  2. подобрать параметры всех звеньев модели в правильном порядке (предполагаем, что сигнал позволяет угадать топологию модели, либо что она нам уже известна)
  3. проверить, что после вычета систематических составляющих из входного сигнала, получится нормально распределенный шум
  4. подставить параметры в код блока Engee Function и проверить, насколько модели соответствуют друг другу.

Синтез модельных данных

Для синтеза данных, которые мы будем анализировать, запустите следующую ячейку:

In [ ]: 
modelName = "param_identification_model"
#engee.close( modelName, force=true ) # Для повторного запуска модели
model = modelName in [m.name for m in engee.get_all_models()] ? engee.open( modelName ) : engee.load( "$(@__DIR__)/$modelName.engee");

results = engee.run( modelName, verbose=false )
Out[0]:
Dict{String, DataFrames.DataFrame} with 2 entries:
  "from_code"  => 5001×2 DataFrame…
  "from_model" => 5001×2 DataFrame…

Загрузка библиотек

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

In [ ]: 
using CSV, DataFrames, FFTW, Polynomials

Загрузка и подключение библиотеки для проверки гипотез (мы используем тестирование на нормальность):

In [ ]: 
#import Pkg; Pkg.add( "HypothesisTests", io=devnull );
using HypothesisTests

Загрузим данные

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

Мы будем визуализировать сигнал модели ModelSystem, который состоит из вектора времени и вычисленных значений:

In [ ]: 
t = results["from_model"].time;
signal = results["from_model"].value;

Также сохраним частоту дискретизации временного сигнала (величина, обратная расстоянию во времени между двумя соседними измерениями):

In [ ]: 
sample_rate = 1/(t[2] - t[1])
Out[0]:
10.0

Что мы будем идентифицировать

Эталонная модель ModelBlock инициализируется эталонными значениями параметров (теми же, что в модели ModelSystem).

image.png

Модель определяется 8 параметрами:

  • параметры линейного тренда: $T_0 + t \cdot T_1$
  • амплитуда, частота и фаза первой синусоиды: $A_1 sin (2 \pi f_1 t + \varphi_1)$
  • амплитуда, частота и фаза второй синусоиды: $A_2 sin (2 \pi f_2 t + \varphi_2)$

Параметры шума нас интересовать не будут, хотя их легко идентифицировать функциями mean и var.

Идентифицируем линейную составляющую

Первым делом мы идентифицируем и устраним из сигнала линейный тренд.

In [ ]: 
trend = Polynomials.fit( t, signal, 1 ) # Последний аргумент: степень полинома
Out[0]:
5.267818747854857 + 0.00894393330796865∙x

Сохраним параметры полинома в отдельные переменные

In [ ]: 
T0, T1 = trend[0], trend[1]
Out[0]:
(5.267818747854857, 0.00894393330796865)

Выведем найденный тренд на графике и покажем, как выглядят наши данные за вычетом этого тренда:

In [ ]: 
plot(
    plot( t, [signal trend.(t)], label=["Данные" "Тренд"] ),
    plot(t, signal - trend.(t), lc=:darkblue, label="Без тренда"),
    layout=(2,1)
)
Out[0]:

Отлично, сигнал стал похож на сочетание синусоид и шума.

Идентифицируем первую синусоиду

Чтобы определить параметры входящих в сигнал синусоид, мы построим спектр сигнала при помощи библиотеки FFTW.

Дискретизация которого будет пропорциональна количеству точек в исходном сигнале. То есть, если оставить исходную дискретизацию, мы можем получить низкую точность спектра. Чтобы избавиться от высокочастотной области, сократим количество точек в исходном сигнале, понизив частоту дискретизации исходного сигнала. Для этого мы задаем переменную downsample_factor. Значение 40 означает, что мы возьмем каждую 40-ю точку из исходного сигнала.

In [ ]: 
downsample_factor = 40;

# Дальше мы будем работать с сигналом пониженной дискретизации (signal_downsampled)
idx_downsample = 1:downsample_factor:length(t)
t_downsampled = t[idx_downsample]

# Дискретизация сигнала, из которого мы вычли линейный тренд
signal_downsampled = (signal - trend.(t))[idx_downsample];

Поскольку нас интересует только положительная часть спектра, используем функцию rfft (нам также не придется делать fftshift):

In [ ]: 
# Найдем спектр сигнала
F1 = rfft(signal_downsampled)

# Рассчитаем вектор частот найденного спектра
freqs = rfftfreq(length(t_downsampled), sample_rate/downsample_factor)

# Выведем график
plot( freqs, abs.(F1) )
Out[0]:

Найдем максимум на этом спектре и определим, какой частоте соответствует это значение.

In [ ]: 
_,idx = findmax( abs.(F1) )
print( "Частота наиболее заметной компоненты: ", freqs[idx] )

Частота наиболее заметной компоненты: 0.01984126984126984

Когда мы нашли максимальную точку на спектре, мы можем идентифицировать амплитуду частоту и фазу наиболее заметной гармонической компоненты. Теперь нам нужно сделать рекомбинацию – сгенерировать гармонический сигнал из комплексного числа, которое описывает его параметры.

In [ ]: 
# Находим параметры синусоиды
A1 = 2 * abs( F1[idx] ) ./ length(t_downsampled )
# Умножение на 2 происходит из-за того, что мощность сигнала в спектре
# делится между двумя пиками: положительным и отрицательным

f1 = abs( freqs[idx] ) # Частота синусоиды
phi1 = angle( F1[idx] ) + pi/2 # Фаза синусоиды

print("A1 = ", A1, "\nf1 = ", f1, "\nphi1 = ", phi1)

A1 = 0.9630993921284998
f1 = 0.01984126984126984
phi1 = 0.3848888013993568

Зададим первую синусоиду как функцию от времени и выведем график:

In [ ]: 
# Функция, задающая синусоиду с найденными параметрами
sinusoid1(tt) = A1 * sin.( 2*pi*f1*tt .+ phi1 );

# Выводим графики 
plot(
    plot( t, [signal .- trend.(t), sinusoid1.(t)], label=["Данные без лин.тренда" "Первая синусоида"]  ), 
    plot( t, signal .- trend.(t) .- sinusoid1(t), lc=:darkblue, label="Остаточный сигнал" ),
    layout=(2,1)
)
Out[0]:

Хорошо, теперь сигнал является сочетанием синусоиды и шума.

Идентифицируем вторую синусоиду

Повторим те же операции, что и для поиска первой синусоиды:

In [ ]: 
# Сигнал пониженной дискретизации (без тренда и первой синусоиды)
signal_downsampled2 = (signal .- trend.(t) .- sinusoid1(t))[idx_downsample]

# Находим спектр сигнала после вычета первых двух компонент
F2 = rfft(signal_downsampled2)

# Рассчитаем вектор частот найденного спектра
freqs = rfftfreq(length(t_downsampled), sample_rate/downsample_factor)

# Выведем график
plot( freqs, abs.(F2) )
Out[0]:

Найдем параметры второй синусоиды

In [ ]: 
# Находим максимум на спектре
_,idx2 = findmax( abs.(F2) )

# Определяем параметры синусоиды, соответствующие этому максимуму
A2 = 2 * abs( F2[idx2] ) ./ length( t_downsampled ) # Амплитуда
f2 = abs( freqs[idx2] ) # Частота
phi2 = angle( F2[idx2] ) + pi/2 # Фаза

print( "A2 = ", A2, "\nf2 = ", f2, "\nphi2 = ", phi2 )

A2 = 0.6193504244923119
f2 = 0.007936507936507936
phi2 = 0.13537312956162761

Зададим вторую синусоиду как функцию от времени и выведем график:

In [ ]: 
# Функция, задающая синусоиду с найденными параметрами
sinusoid2(tt) = A2 * sin.( 2*pi*f2*tt .+ phi2 );

# Выводим графики
plot(
    plot(t, [signal.-trend.(t).-sinusoid1.(t) sinusoid2.(t)], 
        label = ["Сигнал минус две компоненты" "Вторая синусоида"]),
    plot(t, signal.-trend.(t).-sinusoid1(t).-sinusoid2(t), lc=:darkblue,
        label = "Остаточный сигнал", legend=:outertopright),
    layout=(2,1)
)
Out[0]:

Как узнать, остались ли в этом сигнале существенные систематические составляющие? Один из способов – представить, что этот сигнал – реализация случайной величины и проверить его на нормальность распределения.

Анализ остаточного сигнала

Проверим гипотезу о нормальном распределении остаточной компоненты сигнала при помощи теста Харке-Бера (JB-теста):

In [ ]: 
h = JarqueBeraTest( signal .- trend.(t) .- sinusoid1(t) .- sinusoid2(t) )
Out[0]:
Jarque-Bera normality test
--------------------------
Population details:
    parameter of interest:   skewness and kurtosis
    value under h_0:         "0 and 3"
    point estimate:          "-0.010841427456395828 and 2.903158191134578"

Test summary:
    outcome with 95% confidence: fail to reject h_0
    one-sided p-value:           0.3584

Details:
    number of observations:         5001
    JB statistic:                   2.05218

Обратите внимание на результат (outcome) теста:

  • fail to reject h_0 означает отсутствие веских оснований отклонять гипотезу $H_0$: с уверенностью 95% мы утверждаем, что наблюдаемый сигнал порожден случайным процессом с нормальным распределением
  • reject h_0 означает, что уверенность в том, что наблюдаемый сигнал происходит из нормального распределения, не превышает 5%

Попробуйте исключить какую-нибудь составляющую из формулы – скорее всего, такой сигнал не пройдет тест Харке-Бера и вы получите результат reject h_0.

Анализ результатов

Полученный результат служит основанием для завершения поиска дополнительных систематических компонент.

Выведем финальный отчет и погрешность оценки параметров.

In [ ]: 
df = DataFrame( 
    ["T0"    T0    T0m    string(round(100*abs(T0 - T0m)/T0m, digits=2))*"%";
     "T1"    T1    T1m    string(round(100*abs(T1 - T1m)/T1m, digits=2))*"%";
     "A1"    A1    A1m    string(round(100*abs(A1 - A1m)/A1m, digits=2))*"%";
     "f1"    f1    f1m    string(round(100*abs(f1 - f1m)/f1m, digits=2))*"%";
     "phi1"  phi1  phi1m  string(round(abs(phi1 - phi1m),digits=2))*" рад";
     "A2"    A2    A2m    string(round(100*abs(A2 - A2m)/A2m, digits=2))*"%";
     "f2"    f2    f2m    string(round(100*abs(f2 - f2m)/f2m, digits=2))*"%";
     "phi2"  phi2  phi2m  string(round(abs(phi2 - phi2m),digits=2))*" рад"
    ],
    ["Параметр", "Идентифицировано", "Модельное значение", "Погрешность"]
 )
Out[0]:

8 rows × 4 columns

ПараметрИдентифицированоМодельное значениеПогрешность
AnyAnyAnyAny
1T05.2678255.36%
2T10.008943930.0110.56%
3A10.96309913.69%
4f10.01984130.020.79%
5phi10.38488900.38 рад
6A20.619350.63.23%
7f20.007936510.0080.79%
8phi20.13537300.14 рад

Создадим функцию, которая рассчитывает данные по модели с идентифицированными параметрами:

In [ ]: 
function ModelBlock(t)
    global T0, T1, A1, f1, phi1, A2, f2, phi2;

    signal = T0 .+ T1 .* t .+
          A1 .* sin.( 2*pi .* f1 .* t + phi1 ) .+
          A2 .* sin.( 2*pi .* f2 .* t + phi2);
    
    return signal
end;

Сравним результат выполнения идентифицированной модели и результат запуска ModelBlock:

In [ ]: 
plot( results["from_model"].time,
          [results["from_model"].value, results["from_code"].value, ModelBlock.(t) ],
      label=["Исходный сигнал" "Эталонная модель" "Идентифицированная модель"],
      linecolor=[:lightblue :red :green])
Out[0]:

В этом примере мы:

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

Блоки, использованные в примере