Валидация
Для помощи в валидации оцененных моделей доступен ряд функций. Покажем их на примере.
Сгенерируем тестовые данные:
using ControlSystemIdentification, ControlSystemsBase, Random
using ControlSystemIdentification: newpem
Random.seed!(1)
T = 200
nx = 2
nu = 1
ny = 1
x0 = randn(nx)
σy = 0.5
sim(sys,u) = lsim(sys, u, 1:T)[1]
sys = tf(1, [1, 2*0.1, 0.1])
sysn = tf(σy, [1, 2*0.1, 0.3])
# Обучающие данные
u = randn(nu,T)
y = sim(sys, u)
yn = y + sim(sysn, randn(size(u)))
dn = iddata(yn, u, 1)
# Проверочные данные
uv = randn(nu, T)
yv = sim(sys, uv)
ynv = yv + sim(sysn, randn(size(uv)))
dv = iddata(yv, uv, 1)
dnv = iddata(ynv, uv, 1)InputOutput data of length 200, 1 outputs, 1 inputs, Ts = 1Затем выполним подгонку нескольких моделей:
res = [newpem(dn, nx, focus=:prediction) for nx = [2,3,4]];Iter Function value Gradient norm
0 8.214820e+01 4.786930e+02
* time: 6.890296936035156e-5
50 1.367893e+01 1.347632e+02
* time: 0.03708004951477051
Iter Function value Gradient norm
0 6.941423e+01 6.057682e+02
* time: 7.796287536621094e-5
50 1.905270e+01 8.542130e+02
* time: 0.07334589958190918
Iter Function value Gradient norm
0 5.781936e+01 6.010229e+02
* time: 8.106231689453125e-5
50 7.805465e+00 4.607881e+01
* time: 0.06663203239440918
100 7.798826e+00 5.230949e+01
* time: 0.12460589408874512После подгонки моделей проверим результаты с помощью проверочных данных и функций simplot и predplot (сравните их с функцией Matlab sys.id compare):
using Plots
ω = exp10.(range(-2, stop=log10(pi), length=150))
fig = plot(layout=4, size=(1000,600))
for i in eachindex(res)
sysh, x0h, opt = res[i]
simplot!( sysh, dnv, x0h; sp=1, ploty=false)
predplot!(sysh, dnv, x0h; sp=2, ploty=false)
end
plot!(dnv.y' .* [1 1], lab="y", l=(:dash, :black), sp=[1 2])
bodeplot!((getindex.(res,1)), ω, link = :none, balance=false, plotphase=false, subplot=3, title="Process", linewidth=2*[4 3 2 1])
bodeplot!(innovation_form.(getindex.(res,1)), ω, link = :none, balance=false, plotphase=false, subplot=4, linewidth=2*[4 3 2 1])
bodeplot!(sys, ω, link = :none, balance=false, plotphase=false, subplot=3, lab="True", l=(:black, :dash), legend = :bottomleft, title="System model")
bodeplot!(innovation_form(ss(sys),syse=ss(sysn)), ω, link = :none, balance=false, plotphase=false, subplot=4, lab="True", l=(:black, :dash), ylims=(0.1, 100), legend = :bottomleft, title="Noise model")
На рисунке результат моделирования сравнивается с истинной моделью вверху слева и с прогнозом вверху справа. Модели системы и модели шума визуализируются на графиках внизу. Все модели достаточно хорошо отражают динамику системы, но испытывают некоторые трудности с представлением коэффициента усиления динамики шума. У истинной системы четыре полюса (два в рамках процесса и еще два в рамках процесса шума), но иногда более простая модель может работать лучше.
Прогнозные модели можно также оценивать путем h-шагового прогнозирования, где h означает «горизонт».
figh = plot()
for i in eachindex(res)
sysh, x0h, opt = res[i]
predplot!(sysh, dnv, x0h, ploty=false, h=5)
end
plot!(dnv.y', lab="y", l=(:dash, :black))
figh
Обычно желательно проверять оцененную модель с горизонтом прогнозирования более единицы. В частности, может быть полезно проверить эффективность горизонта прогнозирования, приблизительно соответствующего главной временной постоянной процесса.
См. также описание функций simulate, predplot, simplot и coherenceplot.
Предикторы различной длины
При увеличении горизонта прогнозирования сопоставление приближается к соответствующему сопоставлению моделируемой системы, в то время как сопоставление становится все меньше.
using LinearAlgebra
G = c2d(DemoSystems.resonant(), 0.1)
K = kalman(G, I(G.nx), I(G.ny))
sys = add_input(G, K, I(G.ny)) # Строим модель инноваций с входами u и e
T = 10000
u = randn(G.nu, T)
e = 0.1randn(G.ny, T)
y = lsim(sys, [u; e]).y
d = iddata(y, u, G.Ts)
Gh,_ = newpem(d, G.nx, zeroD=true)
# Создаем предикторы с различными горизонтами
p1 = observer_predictor(Gh)
p2 = observer_predictor(Gh, h=2)
p10 = observer_predictor(Gh, h=10)
p100 = observer_predictor(Gh, h=100)
bodeplot([p1, p2, p10, p100], plotphase=false, lab=["1" "" "2" "" "10" "" "100" ""])
bodeplot!(sys, ticks=:default, plotphase=false, l=(:black, :dash), lab=["sim" ""], title=["From u" "From y"])
Погрешность прогнозирования как функция от горизонта прогнозирования приближается к погрешности моделирования.
using Statistics
hs = [1:40; 45:5:80]
perrs = map(hs) do h
yh = predict(Gh, d; h)
ControlSystemIdentification.rms(d.y - yh) |> mean
end
serr = ControlSystemIdentification.rms(d.y - simulate(Gh, d)) |> mean
plot(hs, perrs, lab="Prediction errors", xlabel="Horizon", ylabel="RMS error")
hline!([serr], lab="Simulation error", l=:dash, legend=:bottomright, ylims=(0, Inf))
API валидации
#
StatsAPI.predict — Function
predict(ARX::TransferFunction, d::InputOutputData)Прогнозирование на шаг вперед для процесса ARX. Длина возвращаемого прогнозирования равна length(d) - max(na, nb).
Пример:
julia> predict(tf(1, [1, -1], 1), iddata(1:10, 1:10))
9-element Vector{Int64}:
2
4
6
8
10
12
14
16
18predict(sys, d::AbstractIdData, args...)
predict(sys, y, u, x0 = nothing)См. также описание predplot.
yh = predict(ar::TransferFunction, y)Прогнозирование авторегрессионной модели
#
ControlSystemIdentification.fpe — Function
fpe(e, d::Int)Критерий итоговой ошибки прогноза Акаике (FPE) для выбора порядка модели.
e — ошибки прогнозирования, а d — количество оцениваемых параметров.