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

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

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

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

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

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

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

В подсистеме 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

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

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

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

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

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

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

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

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

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

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

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

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

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( model )

Dict{String, DataFrames.DataFrame} with 2 entries:
  "from_code"  => 5001×2 DataFrame…
  "from_model" => 5001×2 DataFrame…

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

Подключение библиотеки для вывода графиков:

using Plots

Выбор движка для динамической отрисовки графиков:

plotly()

[ Info: Precompiling PlotlyKaleido [f2990250-8cf9-495f-b13a-cce12b45703c]

Plots.PlotlyBackend()

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

using CSV, DataFrames, FFTW, Polynomials

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

import Pkg; Pkg.add( "HypothesisTests", io=devnull );
using HypothesisTests

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

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

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

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

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

sample_rate = 1/(t[2] - t[1])

10.0

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

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

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

параметры линейного тренда:
амплитуда, частота и фаза первой синусоиды:
амплитуда, частота и фаза второй синусоиды:

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

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

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

trend = Polynomials.fit( t, signal, 1 ) # Последний аргумент: степень полинома

Polynomial(5.267818747854857 + 0.00894393330796865*x)

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

T0, T1 = trend[0], trend[1]

(5.267818747854857, 0.00894393330796865)

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

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

interactive-scripts/images/data_analysis_param_identification_example/a24c404df226b877f73b473bbe8ce1bdd275c7b8

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

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

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

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

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):

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

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

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

interactive-scripts/images/data_analysis_param_identification_example/92ac506f93d4d3c97988cad8f999acf5ca866117

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

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

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

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

# Находим параметры синусоиды
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

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

# Функция, задающая синусоиду с найденными параметрами
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)
)

interactive-scripts/images/data_analysis_param_identification_example/4db8730c4ea19abe496332fa7e5a068371df764d

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

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

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

# Сигнал пониженной дискретизации (без тренда и первой синусоиды)
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) )

interactive-scripts/images/data_analysis_param_identification_example/9309314531d0a6fb60529aefc5e9066688c65337

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

# Находим максимум на спектре
_,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

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

# Функция, задающая синусоиду с найденными параметрами
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)
)

interactive-scripts/images/data_analysis_param_identification_example/f97d2bff9268c3df43cf45138405f27318611fe7

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

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

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

h = JarqueBeraTest( signal .- trend.(t) .- sinusoid1(t) .- sinusoid2(t) )

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 означает отсутствие веских оснований отклонять гипотезу : с уверенностью 95% мы утверждаем, что наблюдаемый сигнал порожден случайным процессом с нормальным распределением
reject h_0 означает, что уверенность в том, что наблюдаемый сигнал происходит из нормального распределения, не превышает 5%

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

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

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

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

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))*" рад"
    ],
    ["Параметр", "Идентифицировано", "Модельное значение", "Погрешность"]
 )

8×4 DataFrame
 Row │ Параметр  Идентифицировано  Модельное значение  Погрешность
     │ Any       Any               Any                 Any
─────┼─────────────────────────────────────────────────────────────
   1 │ T0        5.26782           5                   5.36%
   2 │ T1        0.00894393        0.01                10.56%
   3 │ A1        0.963099          1                   3.69%
   4 │ f1        0.0198413         0.02                0.79%
   5 │ phi1      0.384889          0                   0.38 рад
   6 │ A2        0.61935           0.6                 3.23%
   7 │ f2        0.00793651        0.008               0.79%
   8 │ phi2      0.135373          0                   0.14 рад

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

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:

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

interactive-scripts/images/data_analysis_param_identification_example/d686cde6b17940c90a79a3ea6abfb4facf831f1f

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

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

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

Engee Function