Стеклоплавильная печь
В этом примере оценивается модель стеклоплавильной печи.
Данные взяты из базы данных идентификации STADIUS.
using DelimitedFiles, Plots
using ControlSystemIdentification, ControlSystemsBase
url = "https://ftp.esat.kuleuven.be/pub/SISTA/data/process_industry/glassfurnace.dat.gz"
zipfilename = "/tmp/furnace.dat.gz"
path = Base.download(url, zipfilename)
run(`gunzip -f $path`)
data = readdlm(path[1:end-3])
u = data[:, 2:4]'
y = data[:, 5:10]'
d = iddata(y, u, 1)InputOutput data of length 1247, 6 outputs, 3 inputs, Ts = 1Входы включают два входа нагрева и один вход охлаждения, а 6 выходов — это датчики температуры в поперечном сечении печи.
Прежде чем оценивать модель, рассмотрим данные.
plot(d, layout=9)
Мы разделим данные пополам и используем первую половину для оценки, а вторую для проверки. Для регистрации прямой передачи на выход 4 в системе требуется параметр zeroD=false, иначе степень соответствия для выхода 4 никогда не будет достаточно высокой.
dtrain = d[1:2end÷3]
dval = d[2end÷3:end]
model = subspaceid(dtrain, 7, zeroD=false)Мы можем посмотреть на матрицу в оцениваемой модели:
model.D6×3 Matrix{Float64}:
-0.0162229 -0.0143999 0.214359
-0.00356019 -0.00725287 0.0991567
-0.00343267 -0.00755522 0.0986383
-0.0712456 -0.0359864 0.665711
-0.011395 -0.0131599 0.127809
-0.0103743 -0.00876221 0.133655Действительно, элемент (4,3) имеет достаточно большое значение.
Проверим модель путем прогнозирования на основе проверочных данных:
predplot(model, dval, h=1, layout=6)
predplot!(model, dval, h=10, ploty=false)
На изображениях выше показан результат прогнозирования на шагов в будущее.
Оцененную модель можно также визуализировать в частотной области.
w = exp10.(LinRange(-3, log10(pi/d.Ts), 200))
sigmaplot(model.sys, w, lab="MOESP")