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