В этом примере показано, как разработать и анализировать простые модели из действительной лаборатории, обрабатывают данные. Мы запускаем с маленького описания процесса, изучаем, как импортировать данные к тулбоксу и предварительно обработать/обусловить его и затем продолжить систематически оценивать параметрические и непараметрические модели. Если модели были идентифицированы, мы сравниваем предполагаемые модели и также подтверждаем модель к фактическим выходным данным из эксперимента.
Это тематическое исследование касается данных, собранных от лабораторных весов "фен". (Преподаватель Процесса обратной связи 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.