工程模型参数的识别

示例的描述

让我们考虑基于数据的模型参数化识别问题的实现。我们必须确定某个模型的几个时不变参数，其结构是我们事先知道的。

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

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

在子系统中 ModelSystem 使用求和，五个信号被组合:

*常量组件：编号5,
*时间分量：斜率为0.01的线,
*一次谐波：振幅为1和角频率的正弦波高兴/高兴,
*不太明显的二次谐波：振幅为0.6和角频率的正弦波高兴/高兴,
*具有零偏移和0.3方差的白噪声。

在子系统中 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);
返回信号

结束

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

Pkg.add(["CSV", "HypothesisTests"])

# "Идеальные" параметры линейной части сигнала (смещение, наклон)
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功能块代码并检查模型之间的匹配程度。

模型数据的合成

要综合我们将要分析的数据，请运行以下单元格:

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 )

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

加载库

连接用于下载数据、处理表、分析谱的库以及用于多项式回归的库。

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 ) # Последний аргумент: степень полинома

将多项式的参数保存在单独的变量中

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

很好，信号看起来像正弦波和噪声的组合。

我们确定了第一个正弦曲线

为了确定信号中包含的正弦曲线的参数，我们将使用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) )

让我们找到这个频谱上的最大值，并确定这个值对应于哪个频率。

_,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)
)

好了，现在信号是正弦波和噪声的组合。

识别第二正弦

让我们重复搜索第一个正弦曲线的相同操作。:

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

让我们找到第二个正弦曲线的参数

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

我如何发现这个信号中是否仍然有重要的系统成分？一种方法是想象这个信号是一个随机变量的实现，并检查它的分布的正态性。

残差信号分析

让我们使用[testхарке-Бера]测试关于信号残余分量正态分布的假设（https://ru.wikipedia.org/wiki/Тест_Харке_—_Бера )(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%

尝试从公式中排除某些成分-最有可能的是，这样的信号将无法通过[测试Фарке-Бера]（https://ru.wikipedia.org/wiki/Тест_Харке_—_Бера ），你会得到结果 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))*" рад"
    ],
    ["Параметр", "Идентифицировано", "Модельное значение", "Погрешность"]
 )

让我们创建一个函数，从具有识别参数的模型计算数据。:

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

在这个例子中，我们:
*实现了一个相当简单的参数识别程序,
*已确定所有参数的识别错误,
*我们确保具有识别参数的模型的图与参考模型的图非常相似。

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