Гибкий манипулятор робота

В этом примере оценивается модель гибкого манипулятора робота.

Данные взяты из базы данных идентификации STADIUS.

using DelimitedFiles, Plots
using ControlSystemIdentification, ControlSystemsBase

url = "https://ftp.esat.kuleuven.be/pub/SISTA/data/mechanical/robot_arm.dat.gz"
zipfilename = "/tmp/flex.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.01) # Интервал выборки, не указанный для данных; предположительно 0,01

InputOutput data of length 1024, 1 outputs, 1 inputs, Ts = 0.01

Входом является крутящий момент двигателя, а выходом — ускорение манипулятора.

Прежде чем оценивать модель, рассмотрим данные и функцию когерентности.

xt = [2,5,10,20,50]
plot(
    plot(d),
    coherenceplot(d, xticks=(xt,xt), hz=true),
)

Когерентность низкая при высоких частотах и частотах от 2 до 6 Гц. Поэтому в данных диапазонах частот полагаться на оцененную модель следует с осторожностью. Причиной низкой когерентности может быть плохое отношение сигнала к шуму или наличие нелинейностей. Для систем с антирезонансами, таких как рассматриваемая, отношение сигнала к шуму часто является плохим на частотах режекции (ведь, согласно определению, частота режекции — это частота, на которой сигнал почти отсутствует). Спектральные плотности входа и выхода можно исследовать с помощью функции welchplot (см. также описание функции specplot).

welchplot(d, yticks=(xt,xt))

Неудивительно, что на входе мощность выше 20 Гц очень низкая. Именно это является причиной низкой когерентности в данном диапазоне. Обычно бывает полезно ограничить полосу частот сигнала возбуждения: на высоких частотах в механических конструкциях часто проявляются резонансы высокого порядка и нелинейное поведение, о чем свидетельствует спектральная плотность выхода на частоте примерно 34 Гц.

Кроме того, мы разделим данные пополам и используем первую половину для оценки, а вторую для проверки. Для оценки мы применим алгоритм идентификации на основе подпространства.

dtrain = d[1:end÷2]
dval = d[end÷2:end]

# Разумной представляется модель порядка 4, модель двойной массы. Мы оценим две модели: одну методом идентификации на основе подпространства, а другую методом на основе погрешности прогнозирования.
model_ss = subspaceid(dtrain, 4, focus=:prediction)
model_pem, _ = newpem(dtrain, 4; sys0 = model_ss, focus=:prediction)

predplot(model_ss, dval, h=1)
predplot!(model_ss, dval, h=5, ploty=false)
simplot!(model_ss, dval, ploty=false)

На изображениях выше показан результат прогнозирования на шагов в будущее.

Оцененные модели можно также визуализировать в частотной области. Показаны как модель, оцененная с помощью PEM, так и непараметрическая оценка методом Фурье (tfest, этот метод также используется для оценки модели шума).

w = exp10.(LinRange(-1, log10(pi/d.Ts), 200))
bodeplot(model_pem.sys, w, lab="PEM", plotphase=false, hz=true)
bodeplot!(model_ss.sys, w, lab="Subspace", plotphase=false, hz=true)
plot!(tfest(d), legend=:bottomleft, hz=true, xticks=(xt,xt))

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

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

using Statistics
hs = [1:40; 45:5:80]
perrs_pem = map(hs) do h
    yh = predict(model_pem, d; h)
    ControlSystemIdentification.rms(d.y - yh) |> mean
end
perrs_ss = map(hs) do h
    yh = predict(model_ss, d; h)
    ControlSystemIdentification.rms(d.y - yh) |> mean
end
serr_pem = ControlSystemIdentification.rms(d.y - simulate(model_pem, d)) |> mean
serr_ss = ControlSystemIdentification.rms(d.y - simulate(model_ss, d)) |> mean

plot(hs, perrs_pem, lab="Prediction errors PEM", xlabel="Prediction Horizon", ylabel="RMS error")
plot!(hs, perrs_ss, lab="Prediction errors Subspace")
hline!([serr_pem], lab="Simulation error PEM", l=:dash, c=1, ylims=(0, Inf))
hline!([serr_ss], lab="Simulation error Subspace", l=:dash, c=2, legend=:bottomright, ylims=(0, Inf))

Как видим, модель на основе погрешности прогнозирования немного лучше справляется с прогнозированием в несколько шагов (ведь именно для этого предназначен метод PEM), в то время как модель, идентифицированная с помощью subspaceid, лучше проявляется себя в разомкнутом контуре. Эффективность моделирования можно еще улучшить, задав параметр focus=:prediction при оценке моделей.