Фен для волос

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

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

Развернутую форму этого примера, в которой для оцениваемой системы проектируется MPC-регулятор, см. здесь, а также в этом видео.

using DelimitedFiles, Plots
using ControlSystemIdentification, ControlSystemsBase

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

InputOutput data of length 1000, 1 outputs, 1 inputs, Ts = 0.08

Входом является напряжение, подаваемое на нагревательный элемент, а выходом — температура воздуха.

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

plot(
    plot(d),
    coherenceplot(d),
)

Когерентность выглядит неплохо до 10 рад/с; выше этой частоты модель становится недостоверной. Мы видим, что данные не отцентрированы относительно нуля, и поэтому, прежде чем продолжить, удаляем среднее значение (при этом не забудьте, что следует сделать с рабочей точкой). Если посмотреть на импульсную характеристику:

d = detrend(d)
impulseestplot(d, 40, σ=3)

то можно увидеть, что имеется входная задержка в три шага.

Перед оценкой моделей разделим данные на оценочный и проверочный наборы. После этого оценим две разные модели: одну методом на основе погрешности прогнозирования (newpem), а другую стандартным методом наименьших квадратов (arx). При оценке модели ARX мы сообщаем средству оценки известную входную задержку в 3 выборки. Затем мы проверяем оцененные модели, используя их для прогнозирования и моделирования на основе проверочных данных.

dtrain = d[1:end÷2]
dval = d[end÷2:end]
## Разумной представляется модель порядка 3.
model_pem,_ = newpem(dtrain, 3)
model_arx = arx(dtrain, 2, 2, inputdelay=3)

predplot(model_pem, dval, sysname="PEM")
predplot!(model_arx, dval, ploty=false, sysname="ARX")
simplot!(model_pem, dval, ploty=false, sysname="PEM")
simplot!(model_arx, dval, ploty=false, sysname="ARX")

Модели приблизительно сопоставимы: модель PEM немного лучше показывает себя при прогнозировании, а модель ARX — при моделировании. Обратите внимание, что функция newpem принимает именованный аргумент focus, которому можно присвоить значение focus = :simulation, чтобы повысить эффективность моделирования.

Наконец, мы сравниваем модели в частотной области с непараметрической оценкой передаточной функции, полученной методами Фурье (tfest):

w = exp10.(LinRange(-1, log10(pi/d.Ts), 200))
bodeplot(model_pem.sys, w, lab="PEM", plotphase=false)
bodeplot!(model_arx, w, lab="ARX", plotphase=false)
plot!(tfest(d))

Параметрические модели весьма схожи и хорошо согласуются с непараметрической оценкой. Кроме того, мы видим, что непараметрическая оценка становится довольно зашумленной при частоте выше 10 рад/с, чего и следовало ожидать исходя из функции когерентности.

Импульсные характеристики оцененной модели можно сравнить с импульсной характеристикой, оцененной непосредственно на основе приведенных выше данных.

impulseestplot(d, 40, σ=3, lab="Data", seriestype=:steppost)
plot!(impulse(model_pem, 3), lab="PEM")
plot!(impulse(model_arx, 3), lab="ARX")

Импульсная характеристика модели ARX для первых трех выборок строго нулевая, так как при оценке этой модели было задано значение inputdelay=3. Для модели PEM этот параметр не был задан, но он тем не менее был определен на основе данных.

На последнем шаге проверки мы проводим остаточный анализ. Если модель извлекла всю полезную информацию, имеющуюся в данных, остатки должны образовывать последовательность белого шума и между входом и остатками не должно быть корреляции. Для этого имеется функция residualplot.

residualplot(model_pem, dval, lab="PEM")
residualplot!(model_arx, dval, lab="ARX")

Как видим, имеется небольшая корреляция слева в остатках: черные пунктирные линии показывают уровни значимости 95 %. Если степень соответствия модели высокая (то есть остатки небольшие), такая небольшая корреляция обычно не должна вызывать беспокойства, поэтому на этом мы сочтем задачу выполненной.