Шар и балка
В этом примере оценивается модель шара на балке.
Данные взяты из базы данных идентификации STADIUS.
using DelimitedFiles, Plots
using ControlSystemIdentification, ControlSystemsBase
## Шар и балка
url = "https://ftp.esat.kuleuven.be/pub/SISTA/data/mechanical/ballbeam.dat.gz"
zipfilename = "/tmp/bb.dat.gz"
path = Base.download(url, zipfilename)
run(`gunzip -f $path`)
data = readdlm(path[1:end-3])
u = data[:, 1]' # Угол наклона балки
y = data[:, 2]' # Положение шара
d = iddata(y, u, 0.1)InputOutput data of length 1000, 1 outputs, 1 inputs, Ts = 0.1Входом является угол наклона балки, а выходом — положение шара на ней. Это неустойчивый процесс (любому, кто пробовал проконтролировать его, знаком характерный звук падения никелевого шарика на пол).
Прежде чем оценивать модель, рассмотрим данные и функцию когерентности.
plot(
plot(d),
coherenceplot(d),
)
При очень низких и очень высоких частотах когерентность низкая. Так как процесс неустойчивый, данные собираются в замкнутом контуре и на входе не так много энергии DC. Поэтому можно предположить, что восстановить характеристики DC модели будет непросто.
Так как данные собираются в замкнутом контуре, мы используем метод идентификации, не подверженный смещению при наличии обратной связи. Выберем метод на основе погрешности прогнозирования (PEM). Так как процесс неустойчивый, укажем, что подпрограмма идентификации должна допускать неустойчивую модель, задав значение stable=false. Если этого не сделать, функция newpem будет пытаться стабилизировать оцениваемую неустойчивую модель.
Кроме того, мы разделим данные пополам и используем первую половину для оценки, а вторую для проверки.
dtrain = d[1:end÷2]
dval = d[end÷2:end]
# Разумной представляется модель порядка 2–3.
model,_ = newpem(dtrain, 3, stable=false)
predplot(model, dval, h=1)
predplot!(model, dval, h=10, ploty=false)
predplot!(model, dval, h=20, ploty=false)
На изображениях выше показан результат прогнозирования на шагов в будущее. Так как процесс является неустойчивым, модель также является неустойчивой и неправдоподобной [1]. Уже при прогнозировании на 20 шагов проявляется тенденция к неустойчивости.
Оцененные модели можно также визуализировать в частотной области. Показаны как модель, оцененная с помощью PEM, так и непараметрическая оценка методом Фурье (tfest, этот метод также используется для оценки модели шума).
w = exp10.(LinRange(-1.5, log10(pi/d.Ts), 200))
bodeplot(model.sys, w, lab="PEM", plotphase=false)
plot!(tfest(d))
Похоже, что модели не согласуются при низких частотах, что вполне ожидаемо с учетом изложенного выше.