Оценка Гаммерштейна — Винера для нелинейной системы ременного привода

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

Система подробно описывается в этом отчете, а данные доступны по ссылке во фрагменте кода ниже.

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

Оценка модели Винера производится с помощью функции newpem; дополнительные сведения см. в разделе Identification of nonlinear models.

using DelimitedFiles, Plots
using ControlSystemIdentification, ControlSystemsBase

url = "http://www.it.uu.se/research/publications/reports/2017-024/CoupledElectricDrivesDataSetAndReferenceModels.zip"
zipfilename = "/tmp/bd.zip"
cd("/tmp")
path = Base.download(url, zipfilename)
run(`unzip -o $path`)
data = readdlm("/tmp/DATAUNIF.csv", ',')[2:end, 1:4]
iddatas = map(0:1) do ind
    u = data[:, 1 + ind]' .|> Float64 # Вход
    y = data[:, 3 + ind]' .|> Float64 # Выход
    iddata(y, u, 1/50)
end

plot(plot.(iddatas)...)

d = iddatas[1] # Для оценки используется один набор данных
coherenceplot(d)

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

output_nonlinearity = (y, p) -> y .= abs.(y)

nx = 3 # Порядок модели

results = map(1:40) do _ # Этот пример немного сложнее, поэтому мы попробуем сделать инициализацию более случайной.
    sysh, x0h, opt = newpem(d, nx; output_nonlinearity, show_trace=false, focus=:simulation)
    (; sysh, x0h, opt)
end;

(; sysh, x0h, opt) = argmin(r->r.opt.minimum, results) # Находим модель с наименьшей стоимостью

dv = iddatas[2] # Используем второй набор данных для проверки
yh = simulate(sysh, dv, x0h)
output_nonlinearity(yh, nothing) # Необходимо вручную применить нелинейность выхода к модели
plot(dv, lab=["Measured nonlinear output" "Input"], layout=(2,1), xlabel="Time")
plot!(dv.t, yh', lab="Simulation", sp=1, l=:dash)

bodeplot(sysh)

Если все прошло нормально, модель должна достаточно хорошо прогнозировать выход, а оцененная модель должна иметь резонансный пик на частоте примерно 20 рад/с (сравните с рис. 8 в отчете).

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