Оценивающие простые модели из действительной лаборатории обрабатывают данные
Этот пример показывает, как разработать и анализировать простые модели из действительной лаборатории, обрабатывают данные. Мы запускаем с маленького описания процесса, изучаем, как импортировать данные к тулбоксу и предварительно обработать/обусловить его и затем продолжить систематически оценивать параметрические и непараметрические модели. Если модели были идентифицированы, мы сравниваем предполагаемые модели и также подтверждаем модель к фактическим выходным данным из эксперимента.
Системное описание
Это тематическое исследование касается данных, собранных от лабораторных весов "фен". (Преподаватель Процесса обратной связи PT326; См. также страницу 525 в Ljung, 1999). Процесс работает можно следующим образом: Воздух вентилируется через трубу и нагревается во входном отверстии. Температура воздуха измеряется термопарой при выходе. Вход является напряжением по нагревающемуся устройству, которое является только сеткой проводов резистора. Вывод является температурой воздуха выхода, представленной измеренным напряжением термопары.
Подготовка Данных для Анализа
Сначала мы загружаем данные ввода - вывода к MATLAB® Workspace.
load dryer2;
Векторный y2, вывод, содержит 1 000 измерений напряжения термопары, которое пропорционально температуре в воздушном потоке выхода. Векторный u2 содержит 1 000 точек входных данных, состоящих из напряжения, применился к нагревателю. Вход был сгенерирован как бинарная случайная последовательность, которая переключается от одного уровня до другого с вероятностью 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: Импульсный ответ модели FIR, оцененной с помощью 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
Report: [1x1 idresults.ssest]
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]
Чтобы выбрать значения матриц пространства состояний и их 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+15 *
0.1470 0.0957 0.0732
0.2935 0.1998 0.2061
1.2756 0.6196 0.2394
dB =
1.0e+13 *
0.3388
1.1938
4.0974
dC =
1.0e+14 *
2.2060 1.5977 0.3666
dD =
0
dK =
1.0e+13 *
1.3862
4.9801
7.2037
Неуверенность является довольно большой даже при том, что модель соответствует данным об оценке хорошо. Это вызвано тем, что модель сверхпараметризована, то есть, она имеет более свободные параметры, чем, что могло быть однозначно определено. Отклонение параметров в таких случаях не четко определено. Однако это не подразумевает, что модель ненадежна. Мы можем построить время - и частотная характеристика этого графика и просмотреть отклонение как области уверенности, как обсуждено затем.
Анализ предполагаемой модели
Диаграмма Боде сгенерированной модели может быть получена с помощью функции 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 = 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);
Команда spa производит модель IDFRD. Функция 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: Предвещайте ответы m1, m2 и m3, сравненного с непараметрической спектральной моделью gs оценки.
Частотные характеристики из этих трех моделей/методов очень близки. Это указывает, что этот ответ надежен.
Кроме того, годограф Найквиста может анализироваться с областями неуверенности, отмеченными на определенных частотах:
showConfidence(nyquistplot(m1,'b',m2,'r',m3,'c',gs,'g'),3)
Рисунок 13: годографы Найквиста моделей m1, m2, m3 и gs.
Непараметрическая модель gs показывает большую часть неуверенности в ответ.
Дополнительная информация
Для получения дополнительной информации об идентификации динамических систем с System Identification Toolbox посещают страницу информации о продукте System Identification Toolbox.