Испаритель

В этом примере оценивается модель четырехступенчатого испарителя, предназначенного для сокращения содержания воды в продукте, например молоке. Входов три: поток подачи, поток пара на первую ступень испарителя и поток охлаждающей воды. Выходов также три: содержание сухого вещества, поток и температура итогового продукта.

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

Zhu Y., Van Overschee P., De Moor B., Ljung L., Comparison of three classes of identification methods. Proc. of SYSID '94,

using DelimitedFiles, Plots
using ControlSystemIdentification, ControlSystemsBase

url = "https://ftp.esat.kuleuven.be/pub/SISTA/data/process_industry/evaporator.dat.gz"
zipfilename = "/tmp/evaporator.dat.gz"
path = Base.download(url, zipfilename)
run(`gunzip -f $path`)
data = readdlm(path[1:end-3])
# Входы:
# 	u1: поток подачи на первую ступень испарителя
# 	u2: поток пара на первую ступень испарителя
# 	u3: поток охлаждающей воды
# Выходы:
# 	y1: содержание сухого вещества
# 	y2: поток итогового продукта
# 	y3: температура итогового продукта
u = data[:, 1:3]'
y = data[:, 4:6]'
d = iddata(y, u, 1)

InputOutput data of length 6305, 3 outputs, 3 inputs, Ts = 1

Входы включают два входа нагрева и один вход охлаждения, а 6 выходов — это датчики температуры в поперечном сечении печи.

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

plot(d, layout=6)

Мы разделим данные пополам и используем первую половину для оценки, а вторую для проверки. Разумной представляется модель приблизительно 8-го порядка (в оригинальной работе используются порядки 6—13). Для регистрации прямой передачи в системе требуется параметр zeroD=false, иначе степень соответствия никогда не будет достаточно высокой.

dtrain = d[1:3300] # первый эксперимент завершается через 3300 секунд
dval = d[3301:end]

model,_ = newpem(dtrain, 8, zeroD=false)

Iter     Function value   Gradient norm
     0     7.896730e+02     2.039439e+02
 * time: 0.00010085105895996094
    50     7.576979e+02     4.949210e+02
 * time: 6.055285930633545
   100     7.371685e+02     7.425409e+01
 * time: 11.377646923065186
   150     7.314176e+02     9.001736e+01
 * time: 16.527439832687378
   200     7.298205e+02     2.484199e+01
 * time: 21.66542100906372
   250     7.295992e+02     5.194412e+00
 * time: 26.784355878829956
   300     7.295811e+02     9.931482e+00
 * time: 31.901824951171875
   350     7.295732e+02     2.194529e+01
 * time: 37.025290966033936
   400     7.295718e+02     7.056210e-01
 * time: 42.14753603935242
   450     7.295717e+02     4.192907e-05
 * time: 47.282536029815674

predplot(model, dval, h=1, layout=d.ny)
predplot!(model, dval, h=5, ploty=false)

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

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

w = exp10.(LinRange(-2, log10(pi/d.Ts), 200))
sigmaplot(model.sys, w, lab="PEM", plotphase=false)

Сравним эффективность прогнозирования с полученной в оригинальной работе:

ys = predict(model, dval, h=5)
ControlSystemIdentification.mse(dval.y-ys)

3×1 Matrix{Float64}:
 0.05761868851045837
 0.1563550175486351
 0.019269556556454893

Авторы получили следующие погрешности: [0,24, 0,39, 0,14]