Документация Engee

Стеклоплавильная печь

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

Данные взяты из базы данных идентификации 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)
PJkyc3dusITRUiVwRD0CzligyERsHq1aslEsnXX3+9YsUKLsXV1fWrr75q3FYRmjpErgiGoPnJFfEsUx9YlkUIIYQ0ExmG0fp6XFdUKlVKSopCoXBycnJycvp3bSQ0PYhcEQwBkavXQgZCAoFAILRoiLEMgUAgEFo0ZCAkEAgEQouGDIQEAoFAaNGQgZBAIBAILRoyEBIIBAKhRUMGQgKBQCC0aP4PKiFldYbGdEEAAAAASUVORK5CYII=

Мы разделим данные пополам и используем первую половину для оценки, а вторую для проверки. Для регистрации прямой передачи на выход 4 в системе требуется параметр zeroD=false, иначе степень соответствия для выхода 4 никогда не будет достаточно высокой.

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

model = subspaceid(dtrain, 7, zeroD=false)

Мы можем посмотреть на матрицу в оцениваемой модели:

model.D
6×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)
4fydK0s14uqy8AAAAASUVORK5CYII=

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

Оцененную модель можно также визуализировать в частотной области.

w = exp10.(LinRange(-3, log10(pi/d.Ts), 200))
sigmaplot(model.sys, w, lab="MOESP")
7b8BwICDIARg4JEkqaOj85RphmQgJydnwYIFJiYm8fHx1tbWT994ypQpj5yWHFgikeiZXUokQSsQyBKcGgUAAKDQoPsEAAAAhQZBCAAAQKH9H8GUMFYv83hDAAAAAElFTkSuQmCC