Engee 文档
Notebook

Engee 中模型参数的识别

示例说明

让我们来考虑基于数据的模型参数识别任务的实现。我们必须识别某个模型的几个时变参数，我们事先知道该模型的结构。

在本例中，我们将分析动态模型 Engee 生成的信号。它位于start/examples/data_analysis/parametric_identification 文件夹中。

该模型由两个子系统组成：一个生成分析信号，另一个是参考模型（无噪声）。

image.png

在子系统ModelSystem 中，五个信号通过求和的方式进行组合：

  • 常数分量：数字 5、
  • 时间分量：斜率为 0.01 的直线、
  • 第一次谐波：振幅为 1、角频率为$2\pi \cdot 0.02$ rad/s 的正弦波、
  • 不太突出的二次谐波：振幅为 0.6、角频率为$2\pi \cdot 0.008$ rad/s 的正弦波、
  • 白噪声，偏移为零，色散为 0.3。

image.png

在子系统ModelBlock 中，我们可以替换已确定的参数。只要在其中设置参考参数并生成参考信号即可。在该子系统中，模型被描述为一个代码：

``朱莉娅 函数 (c::ModelBlock)(t::Real) T0 = 5 T1 = 0.01 A1、f1、phi1 = 1、0.02、0 A2、f2、phi2 = 0.6、0.008、0

信号 = T0 + T1 * t +
         A1 * sin( 2*pi * f1 * t + phi1 ) .+
         A2 * sin( 2*pi * f2 * t + phi2 )；
返回信号

结束

让我们保存参数的参考值：

In [ ]: 
Pkg.add(["CSV", "HypothesisTests"])
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 功能块代码中，并检查模型是否相互匹配。

合成模型数据

运行以下单元格，合成我们要分析的数据：

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

我们怎样才能知道这个信号中是否还存在重要的系统成分呢？一种方法是将该信号想象成随机变量的实现，并测试其分布的正态性。

分析残差信号

让我们使用 Harker-Bera 检验 （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

试着从公式中排除某些成分--这样的信号很可能无法通过 Harkke-Bera 检验，你将得到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]:

在这个例子中，我们

  • 实施了一个相当简单的参数识别程序、
  • 确定所有参数的识别误差、
  • 验证带识别参数的模型图与参考模型图非常相似。

示例中使用的块