Оценка простых моделей на основе реальных данных лабораторного процесса
Этот пример показывает, как разработать и проанализировать простые модели из реальных данных лабораторного процесса. Мы начинаем с небольшого описания процесса, узнаем, как импортировать данные в тулбокс и обработаем/обусловим его, а затем систематически продолжаем оценивать параметрические и непараметрические модели. Как только модели были идентифицированы, мы сравниваем предполагаемые модели, а также подтверждаем модель с фактическими выходными данными эксперимента.
Описание системы
Это тематическое исследование касается данных, собранных с лабораторной шкалы «фен». (Feedback's Process Trainer PT326; См. также стр. 525 в Ljung, 1999). Процесс работает следующим образом: Воздух раздувается через трубку и нагревается на входном отверстии. Температура воздуха измеряется термопарой на выходе. Входным входом является напряжение на нагревательном устройстве, которое является просто mesh резисторных проводов. Выходным сигналом является температура воздуха на выходе, представленная измеренным напряжением термопары.
Настройка данных для анализа
Сначала загружаем данные ввода-вывода в Рабочую область MATLAB ®.
load dryer2;
Векторные y2, выход, содержит 1000 измерений напряжения термопары, пропорционального температуре в выходном воздушном потоке. Векторные u2 содержит 1000 входных точек данных, состоящих из напряжения, приложенного к нагревателю. Вход был сгенерирован как двоичная случайная последовательность, которая переключается с одного уровня на другой с вероятностью 0,2. Значение шага расчета составляет 0.08 секунд.
Следующим шагом является настройка данных как объекта iddata
dry = iddata(y2,u2,0.08);
Чтобы получить информацию о данных, просто введите имя iddata объект в командном окне MATLAB:
dry
dry =
Time domain data set with 1000 samples.
Sample time: 0.08 seconds
Outputs Unit (if specified)
y1
Inputs Unit (if specified)
u1
Чтобы просмотреть свойства вышеуказанного объекта iddata, используйте get команда:
get(dry)
ans =
struct with fields:
Domain: 'Time'
Name: ''
OutputData: [1000x1 double]
y: 'Same as OutputData'
OutputName: {'y1'}
OutputUnit: {''}
InputData: [1000x1 double]
u: 'Same as InputData'
InputName: {'u1'}
InputUnit: {''}
Period: Inf
InterSample: 'zoh'
Ts: 0.0800
Tstart: []
SamplingInstants: [1000x0 double]
TimeUnit: 'seconds'
ExperimentName: 'Exp1'
Notes: {}
UserData: []
Для лучшего ведения бухгалтерского учета рекомендуется давать имена входным и выходным каналам и модулям времени. Эти имена будут распространяться на протяжении всего анализа этого объекта iddata:
dry.InputName = 'Heater Voltage'; dry.OutputName = 'Thermocouple Voltage'; dry.TimeUnit = 'seconds'; dry.InputUnit = 'V'; dry.OutputUnit = 'V';
Теперь, когда у нас готов набор данных, мы выбираем первые 300 точек данных для оценки модели.
ze = dry(1:300)
ze =
Time domain data set with 300 samples.
Sample time: 0.08 seconds
Outputs Unit (if specified)
Thermocouple Voltage V
Inputs Unit (if specified)
Heater Voltage V
Предварительная обработка данных
Постройте график интервала от выборки 200 до 300:
plot(ze(200:300));
Фигура 1: Снимок измеренных данных о фене.
Из приведенного выше графика можно заметить, что данные не являются средними нулями. Поэтому давайте удалим постоянные уровни и сделаем данные нулевыми средними.
ze = detrend(ze);
Тот же набор данных после того, как он был детрендирован:
plot(ze(200:300)) %show samples from 200 to 300 of detrended data
Фигура 2: Детрендированные данные оценки.
Оценка непараметрических и параметрических моделей
Теперь, когда набор данных был детрендирован и нет явных выбросов, давайте сначала оценим импульсную характеристику системы путем корреляционного анализа, чтобы получить некоторое представление о временных константах и тому подобное:
clf mi = impulseest(ze); % non-parametric (FIR) model showConfidence(impulseplot(mi),3); %impulse response with 3 standard %deviations confidence region
Фигура 3: Импульсная характеристика модели конечной импульсной характеристики, оцененная с использованием ze.
Затененная область отмечает 99,7% доверительный интервал. Существует временная задержка (потеря времени) из 3 выборок, прежде чем выход ответит на вход (значительный выход за пределами доверительного интервала).
Самый простой способ начать работу с параметрической стандартной программой оценки - создать модель пространства состояний, где порядок модели определяется автоматически, используя метод ошибки предсказания. Давайте оценим модель, используя ssest метод оценки:
m1 = ssest(ze);
m1 является идентифицированной в непрерывном времени моделью пространства состояний, представленной idss объект. Алгоритм оценки выбирает 3 в качестве оптимального порядка модели. Чтобы просмотреть свойства предполагаемой модели, достаточно ввести имя модели в командном окне:
m1
m1 =
Continuous-time identified state-space model:
dx/dt = A x(t) + B u(t) + K e(t)
y(t) = C x(t) + D u(t) + e(t)
A =
x1 x2 x3
x1 -0.4839 -2.011 2.092
x2 3.321 -1.913 5.998
x3 1.623 -17.01 -15.61
B =
Heater Volta
x1 -0.05753
x2 0.02004
x3 1.377
C =
x1 x2 x3
Thermocouple -14.07 0.07729 0.04252
D =
Heater Volta
Thermocouple 0
K =
Thermocouple
x1 -0.9457
x2 -0.02097
x3 2.102
Parameterization:
FREE form (all coefficients in A, B, C free).
Feedthrough: none
Disturbance component: estimate
Number of free coefficients: 18
Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties.
Status:
Estimated using SSEST on time domain data "ze".
Fit to estimation data: 95.32% (prediction focus)
FPE: 0.001621, MSE: 0.001526
Отображение предполагает, что модель является свободной формой (все записи матриц A, B и C были рассмотрены как свободные параметры) и что предполагаемая модель довольно хорошо подходит для данных (более 90% подгонки). Чтобы извлечь свойства этой модели, например, чтобы получить A матрица дискретного объекта пространства состояний, сгенерированного выше, мы можем использовать оператор точки:
A = m1.a;
Смотрите пример «Данные и объекты модели в System Identification Toolbox» для получения дополнительной информации об объектах модели. Чтобы узнать, какие свойства объекта модели можно извлечь, используйте get команда:
get(m1)
A: [3x3 double]
B: [3x1 double]
C: [-14.0706 0.0773 0.0425]
D: 0
K: [3x1 double]
StateName: {3x1 cell}
StateUnit: {3x1 cell}
Structure: [1x1 pmodel.ss]
NoiseVariance: 1.2587e-04
InputDelay: 0
OutputDelay: 0
Ts: 0
TimeUnit: 'seconds'
InputName: {'Heater Voltage'}
InputUnit: {'V'}
InputGroup: [1x1 struct]
OutputName: {'Thermocouple Voltage'}
OutputUnit: {'V'}
OutputGroup: [1x1 struct]
Notes: [0x1 string]
UserData: []
Name: ''
SamplingGrid: [1x1 struct]
Report: [1x1 idresults.ssest]
Чтобы получить значения матриц пространства состояний и их 1 стандартная неопределенность отклонения, используйте idssdata команда:
[A,B,C,D,K,~,dA,dB,dC,dD,dK] = idssdata(m1)
A =
-0.4839 -2.0112 2.0917
3.3205 -1.9135 5.9981
1.6235 -17.0096 -15.6070
B =
-0.0575
0.0200
1.3770
C =
-14.0706 0.0773 0.0425
D =
0
K =
-0.9457
-0.0210
2.1019
dA =
1.0e+14 *
1.1413 1.1945 0.8367
2.0899 1.4613 1.2117
9.8253 7.1186 1.8997
dB =
1.0e+13 *
0.4013
1.5604
4.5776
dC =
1.0e+14 *
2.1049 1.3852 0.4309
dD =
0
dK =
1.0e+13 *
1.3260
4.5015
7.0925
Неопределенности довольно большие, хотя модель хорошо соответствует данным оценки. Это потому, что модель чрезмерно параметризована, то есть она имеет больше свободных параметров, чем то, что можно однозначно идентифицировать. Отклонение параметров в таких случаях не четко определена. Однако это не означает, что модель ненадежна. Мы можем построить график временной и частотной характеристики этого графика и просмотреть отклонение как доверительные области, как обсуждается далее.
Анализ предполагаемой модели
Диаграмма Боде сгенерированной модели может быть получен с помощью bode функция, как показано ниже:
h = bodeplot(m1);
Фигура 4: Диаграмма Боде предполагаемой модели.
Щелкните правой кнопкой мыши график и выберите «Характеристики» - > «Доверительная область». Или используйте showConfidence команда для просмотра отклонения отклика.
showConfidence(h,3) % 3 standard deviation (99.7%) confidence region
Фигура 5: Диаграмма Боде с 3 стандартным отклонением доверия области.
Точно так же мы можем сгенерировать график шага и связанную с ним 3 стандартную доверительную область отклонения. Мы можем сравнить отклики и связанные отклонения параметрической модели m1 с таковой непараметрической модели mi:
showConfidence(stepplot(m1,'b',mi,'r',3),3)
Фигура 6: Шаговый график m1 моделей и mi с доверительными областями.
Можно также рассмотреть годограф Найквиста и отметить области неопределенности на определенных частотах эллипсами, соответствующими 3 стандартным отклонениям:
Opt = nyquistoptions;
Opt.ShowFullContour = 'off';
showConfidence(nyquistplot(m1,Opt),3)
Фигура 7: Годограф Найквиста предполагаемой модели, показывающей области неопределенности на определенных частотах.
Графики отклика показывают, что предполагаемая модель m1 довольно надежно.
Оценка моделей с предписанной структурой
System Identification Toolbox также может использоваться, чтобы получить модель с предписанной структурой. Для примера модель разностного уравнения с 2 полюсами, 1 нулями и 3 выборка задержками может быть получена с помощью arx функция, как показано ниже:
m2 = arx(ze,[2 2 3]);
Чтобы просмотреть модель, введите имя модели в командном окне.
m2
m2 =
Discrete-time ARX model: A(z)y(t) = B(z)u(t) + e(t)
A(z) = 1 - 1.274 z^-1 + 0.3942 z^-2
B(z) = 0.06679 z^-3 + 0.04429 z^-4
Sample time: 0.08 seconds
Parameterization:
Polynomial orders: na=2 nb=2 nk=3
Number of free coefficients: 4
Use "polydata", "getpvec", "getcov" for parameters and their uncertainties.
Status:
Estimated using ARX on time domain data "ze".
Fit to estimation data: 95.08% (prediction focus)
FPE: 0.001756, MSE: 0.001687
Непрерывная передаточная функция времени с 2 полюсами, одним нулем и 0,2 секундной задержкой переноса может быть оценена с помощью tfest команда:
m3 = tfest(ze, 2, 1, 0.2)
m3 =
From input "Heater Voltage" to output "Thermocouple Voltage":
1.183 s + 26.55
exp(-0.2*s) * ---------------------
s^2 + 11.61 s + 28.63
Continuous-time identified transfer function.
Parameterization:
Number of poles: 2 Number of zeros: 1
Number of free coefficients: 4
Use "tfdata", "getpvec", "getcov" for parameters and their uncertainties.
Status:
Estimated using TFEST on time domain data "ze".
Fit to estimation data: 88.79%
FPE: 0.009126, MSE: 0.008768
Валидация предполагаемой модели к экспериментальному выходу
Насколько хороша предполагаемая модель? Один из способов выяснить это - симулировать его и сравнить результат модели с измеренным выходом. Выберите фрагмент исходных данных, которые не использовались при создании модели, скажем, из выборок от 800 до 900. После предварительной обработки данных валидации мы используем compare функция, показанная ниже, чтобы просмотреть качество предсказания:
zv = dry(800:900); % select an independent data set for validation zv = detrend(zv); % preprocess the validation data set(gcf,'DefaultLegendLocation','best') compare(zv,m1); % perform comparison of simulated output
Фигура 8: Симулированный отклик модели по сравнению с выходными данными валидации.
Здесь можно заметить, что соглашение очень хорошее. Показанное значение «Fit» вычисляется как:
Fit = 100*(1 - norm(yh - y)/norm(y-mean(y)))
где y - измеренный выход (=|zv.y|), и yh - выход модели m1.
Сравнение предполагаемых моделей
Чтобы сравнить эффективность моделей, которые мы оценили, например m1, m2 и m3 с данными валидации zv, мы можем снова использовать compare команда:
compare(zv,m1,'b',m2,'r',m3,'c');
Фигура. 9. Сравнение характеристик моделей m1, m2, m3 на наборе данных валидации ze.
Диаграммы нулей и полюсов для моделей можно получить с помощью iopzplot:
h = iopzplot(m1,'b',m2,'r',m3,'c');
Фигура 10: Поляки и нули моделей m1, m2 и m3.
Неопределенности в полюсах и нулях также могут быть получены. В следующем операторе '3' относится к количеству стандартных отклонений.
showConfidence(h,3);
Фигура 11: Карта ноль полюсов с областями неопределенности.
Частотные функции выше, которые получаются из моделей, могут быть сравнены с той, которая получена с помощью непараметрического метода спектрального анализа (spa):
gs = spa(ze);
The spa команда создает модель IDFRD. The bode функция может снова использоваться для сравнения с передаточными функциями полученных моделей.
w = linspace(0.4,pi/m2.Ts,200); opt = bodeoptions; opt.PhaseMatching = 'on'; bodeplot(m1,'b',m2,'r',m3,'c',gs,'g',w,opt); legend('m1','m2','m3','gs')
Фигура 12: Bode-ответы m1, m2 и m3 по сравнению с непараметрической моделью спектральной оценки gs.
Частотные характеристики трех моделей/методов очень близки. Это указывает, что этот ответ надежен.
Кроме того, годограф Найквиста может быть проанализирован с областями неопределенности, отмеченными на определенных частотах:
showConfidence(nyquistplot(m1,'b',m2,'r',m3,'c',gs,'g'),3)
Фигура 13: Годографы Найквиста моделей m1, m2, m3 и gs.
Непараметрическая модель gs показывает наибольшую неопределенность в ответ.