Сравнение различных методов идентификации модели
Этот пример показывает несколько методов идентификации, доступных в System Identification Toolbox™. Мы начинаем путем симуляции экспериментальных данных и используем несколько методов оценки, чтобы оценить модели из данных. Следующие стандартные программы оценки проиллюстрированы в этом примере: spa, ssest, tfest, arx, oe, armax и bj.
Системное описание
Шум повредил линейную систему, может быть описан:
y = G u + H e
где y выход и u вход, в то время как e обозначает неизмеренный (белый) источник шума. G передаточная функция системы и H дает описание шумового воздействия.
Существует много способов оценить G и H. По существу они соответствуют различным способам параметрировать эти функции.
Определение модели
System Identification Toolbox предоставляет пользователям опцию симуляции данных, как был бы получен из физического процесса. Давайте симулируем данные из модели IDPOLY с данным набором коэффициентов.
B = [0 1 0.5]; A = [1 -1.5 0.7]; m0 = idpoly(A,B,[1 -1 0.2],'Ts',0.25,'Variable','q^-1'); % The sample time is 0.25 s.
Третий аргумент [1 -1 0.2] дает характеристику воздействий, которые действуют на систему. Выполнение этих команд производит следующее дискретное время idpoly модель:
m0
m0 =
Discrete-time ARMAX model: A(q)y(t) = B(q)u(t) + C(q)e(t)
A(q) = 1 - 1.5 q^-1 + 0.7 q^-2
B(q) = q^-1 + 0.5 q^-2
C(q) = 1 - q^-1 + 0.2 q^-2
Sample time: 0.25 seconds
Parameterization:
Polynomial orders: na=2 nb=2 nc=2 nk=1
Number of free coefficients: 6
Use "polydata", "getpvec", "getcov" for parameters and their uncertainties.
Status:
Created by direct construction or transformation. Not estimated.
Здесь q обозначает оператор сдвига так, чтобы (q) y (t) действительно являлся сокращением от y (t) - 1,5 года (t-1) + 0,7 года (t-2). Отображение модели m0 показывает, что это - модель ARMAX.
Эта конкретная структура модели известна как модель ARMAX, где (Авторегрессивный) AR относится к A-полиному, MA (Скользящее среднее значение) к шумовому C-полиному и X к "дополнительному" входу B (q) u (t).
В терминах общих передаточных функций G и H, эта модель соответствует параметризации:
G(q) = B(q)/A(q) и H(q) = C(q)/A(q), с общими знаменателями
Симуляция модели
Мы генерируем входной сигнал u и симулируйте ответ модели к этим входным параметрам. idinput команда может использоваться, чтобы создать входной сигнал к модели и iddata команда группирует сигнал в подходящий формат. sim команда может затем использоваться, чтобы симулировать выход как показано ниже:
prevRng = rng(12,'v5normal'); u = idinput(350,'rbs'); %Generates a random binary signal of length 350 u = iddata([],u,0.25); %Creates an IDDATA object. The sample time is 0.25 sec. y = sim(m0,u,'noise') %Simulates the model's response for this
y =
Time domain data set with 350 samples.
Sample time: 0.25 seconds
Outputs Unit (if specified)
y1
%input with added white Gaussian noise e %according to the model description m0 rng(prevRng);
Один аспект, чтобы отметить то, что сигналы u и y оба объекты IDDATA. Затем мы собираем эти данные ввода - вывода, чтобы сформировать один объект iddata.
z = [y,u];
Чтобы получить информацию об объекте данных (который теперь включает обоих выборки входных и выходных данных), только введите его имя в командном окне MATLAB®:
z
z =
Time domain data set with 350 samples.
Sample time: 0.25 seconds
Outputs Unit (if specified)
y1
Inputs Unit (if specified)
u1
Построить первые 100 значений входа u и выход y, используйте plot команда на объекте iddata:
h = gcf; set(h,'DefaultLegendLocation','best') h.Position = [100 100 780 520]; plot(z(1:100));
Это - хорошая практика, чтобы использовать только фрагмент данных в целях оценки, ze и сохраните другую часть, чтобы подтвердить предполагаемые модели позже:
ze = z(1:200); zv = z(201:350);
Выполнение спектрального анализа
Теперь, когда мы получили симулированные данные, мы можем оценить модели и сделать сравнения. Давайте запустимся со спектрального анализа. Мы вызываем spa команда, которая возвращает оценку спектрального анализа функции частоты и шумового спектра.
GS = spa(ze);
Результатом спектрального анализа является комплексная функция частоты: частотная характеристика. Это упаковано в объект под названием IDFRD объект (Идентифицированные Данные о Частотной характеристике):
GS
GS = IDFRD model. Contains Frequency Response Data for 1 output(s) and 1 input(s), and the spectra for disturbances at the outputs. Response data and disturbance spectra are available at 128 frequency points, ranging from 0.09817 rad/s to 12.57 rad/s. Sample time: 0.25 seconds Output channels: 'y1' Input channels: 'u1' Status: Estimated using SPA on time domain data "ze".
Чтобы построить частотную характеристику, это обычно, чтобы использовать Диаграмму Боде с bode или bodeplot команда:
clf
h = bodeplot(GS); % bodeplot returns a plot handle, which bode does not
ax = axis; axis([0.1 10 ax(3:4)])
Оценка GS сомнительно, поскольку это формируется из поврежденных данных шума. Метод спектрального анализа обеспечивает также свою собственную оценку этой неопределенности. Чтобы отобразить неопределенность (говорят, например, 3 стандартных отклонения), мы можем использовать showConfidence команда на графике обрабатывает h возвращенный предыдущим bodeplot команда.
showConfidence(h,3)
В графике говорится, что несмотря на то, что номинальная оценка частотной характеристики (синяя линия) не обязательно точна, существует вероятность на 99,7% (3 стандартных отклонения нормального распределения), что истинный ответ в теневой (голубой) области.
Оценка параметрических моделей в пространстве состояний
Затем давайте оценим значение по умолчанию (без нас задающий конкретную структуру модели) линейная модель. Это будет вычислено как модель в пространстве состояний с помощью ошибочного метода предсказания. Мы используем ssest функция в этом случае:
m = ssest(ze) % The order of the model will be chosen automaticallym =
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
x1 -0.007167 1.743
x2 -2.181 -1.337
B =
u1
x1 0.09388
x2 0.2408
C =
x1 x2
y1 47.34 -14.4
D =
u1
y1 0
K =
y1
x1 0.04108
x2 -0.03751
Parameterization:
FREE form (all coefficients in A, B, C free).
Feedthrough: none
Disturbance component: estimate
Number of free coefficients: 10
Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties.
Status:
Estimated using SSEST on time domain data "ze".
Fit to estimation data: 75.08% (prediction focus)
FPE: 0.9835, MSE: 0.9262
В этой точке матрицы модели в пространстве состояний не говорят нам очень. Мы могли сравнить частотную характеристику модели m с тем из спектрального анализа оценивают GS. Обратите внимание на то, что GS модель дискретного времени в то время как m модель непрерывного времени. Возможно использовать ssest вычислить модели дискретного времени также.
bodeplot(m,GS) ax = axis; axis([0.1 10 ax(3:4)])
Мы отмечаем, что эти две частотных характеристики близки, несмотря на то, что они были оценены совсем другими способами.
Оценка простых передаточных функций
Существует большое разнообразие линейных структур модели, все соответствующие линейному различию или дифференциальным уравнениям, описывающим отношение между u и y. Отличные структуры все соответствуют различным способам моделирования шумового влияния. Самыми простыми из этих моделей являются модели ARX и передаточная функция. Передаточная функция принимает форму:
Y(s) = NUM(s)/DEN(s) U(s) + E(s)
где NUM и DEN являются полиномом числителя и полиномом знаменателя, и Y (s), U (s) и E (s) являются Преобразования Лапласа выхода, входных сигналов и сигналов ошибки (y (t), u (t) и e (t)) соответственно. NUM и DEN могут быть параметрированы их порядками, которые представляют количество нулей и полюсов. Для данного количества полюсов и нулей, мы можем оценить передаточную функцию с помощью tfest команда:
mtf = tfest(ze, 2, 2) % transfer function with 2 zeros and 2 polesmtf =
From input "u1" to output "y1":
-0.05428 s^2 - 0.02386 s + 29.6
-------------------------------
s^2 + 1.361 s + 4.102
Continuous-time identified transfer function.
Parameterization:
Number of poles: 2 Number of zeros: 2
Number of free coefficients: 5
Use "tfdata", "getpvec", "getcov" for parameters and their uncertainties.
Status:
Estimated using TFEST on time domain data "ze".
Fit to estimation data: 71.69%
FPE: 1.282, MSE: 1.195
Модель AN ARX представлена как: (q) y (t) = B (q) u (t) + e (t),
или обычным письмом,
Эта модель линейна в коэффициентах. Коэффициенты , , ..., , .. может быть вычислен эффективно с помощью метода оценки методом наименьших квадратов. Оценка модели ARX с предписанной структурой - два полюса, один нуль и одна задержка на входе получена как ([na nb nk] = [2 2 1]):
mx = arx(ze,[2 2 1])
mx =
Discrete-time ARX model: A(z)y(t) = B(z)u(t) + e(t)
A(z) = 1 - 1.32 z^-1 + 0.5393 z^-2
B(z) = 0.9817 z^-1 + 0.4049 z^-2
Sample time: 0.25 seconds
Parameterization:
Polynomial orders: na=2 nb=2 nk=1
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: 68.18% (prediction focus)
FPE: 1.603, MSE: 1.509
Сравнение эффективности предполагаемых моделей
Модели (mtf, mx) лучше или хуже, чем модель в пространстве состояний по умолчанию m? Один способ узнать состоит в том, чтобы симулировать каждую модель (свободный шум) с входом из данных о валидации zv и сравните соответствующие симулированные выходные параметры с измеренным выходом в наборе zv:
compare(zv,m,mtf,mx)
Подгонка является процентом выходного изменения, которое объяснено моделью. Очевидно m лучшая модель, несмотря на то, что mtf приближается. Передаточная функция дискретного времени может быть оценкой с помощью любого tfest или oe команда. Первый создает модель IDTF, в то время как последний создает модель IDPOLY структуры Ошибки на выходе (y (t) = B (q)/F (q) u (t) + e (t)), но они математически эквивалентны.
md1 = tfest(ze,2,2,'Ts',0.25) % two poles and 2 zeros (counted as roots of polynomials in z^-1)
md1 = From input "u1" to output "y1": 0.8383 z^-1 + 0.7199 z^-2 ---------------------------- 1 - 1.497 z^-1 + 0.7099 z^-2 Sample time: 0.25 seconds Discrete-time identified transfer function. Parameterization: Number of poles: 2 Number of zeros: 2 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: 71.46% FPE: 1.264, MSE: 1.214
md2 = oe(ze,[2 2 1])
md2 =
Discrete-time OE model: y(t) = [B(z)/F(z)]u(t) + e(t)
B(z) = 0.8383 z^-1 + 0.7199 z^-2
F(z) = 1 - 1.497 z^-1 + 0.7099 z^-2
Sample time: 0.25 seconds
Parameterization:
Polynomial orders: nb=2 nf=2 nk=1
Number of free coefficients: 4
Use "polydata", "getpvec", "getcov" for parameters and their uncertainties.
Status:
Estimated using OE on time domain data "ze".
Fit to estimation data: 71.46%
FPE: 1.264, MSE: 1.214
compare(zv, md1, md2)
Модели md1 и md2 поставьте идентичные подгонки к данным.
Остаточный анализ
Дальнейший способ получить сведения о качестве модели состоит в том, чтобы вычислить "остаточные значения": т.е. та часть e в y = Gu + He это не могло быть объяснено моделью. Идеально, они должны быть некоррелироваными с входом и также взаимно некоррелироваными. Остаточные значения вычисляются, и их свойства корреляции отображены resid команда. Эта функция может использоваться, чтобы оценить остатки оба во временном и частотном диапазоне. Сначала давайте получим остаточные значения для модели Output-Error во временном интервале:
resid(zv,md2) % plots the result of residual analysis
Мы видим, что взаимная корреляция между остаточными значениями и ввела, находится в области доверия, указывая, что нет никакой значительной корреляции. Оценка динамики G должен затем быть рассмотрен как соответствующий. Однако (автоматическая) корреляция e является значительным, таким образом, e не может быть рассмотрен как белый шум. Это означает что шумовая модель H не соответствует.
Оценка ARMAX и моделей поля-Jenkins
Давайте теперь вычислим модель ARMAX второго порядка и модель Box-Jenkins второго порядка. Модель ARMAX имеет те же шумовые характеристики как симулированная модель m0 и модель Box-Jenkins позволяет более общее шумовое описание.
am2 = armax(ze,[2 2 2 1]) % 2nd order ARMAX modelam2 =
Discrete-time ARMAX model: A(z)y(t) = B(z)u(t) + C(z)e(t)
A(z) = 1 - 1.516 z^-1 + 0.7145 z^-2
B(z) = 0.982 z^-1 + 0.5091 z^-2
C(z) = 1 - 0.9762 z^-1 + 0.218 z^-2
Sample time: 0.25 seconds
Parameterization:
Polynomial orders: na=2 nb=2 nc=2 nk=1
Number of free coefficients: 6
Use "polydata", "getpvec", "getcov" for parameters and their uncertainties.
Status:
Estimated using ARMAX on time domain data "ze".
Fit to estimation data: 75.08% (prediction focus)
FPE: 0.9837, MSE: 0.9264
bj2 = bj(ze,[2 2 2 2 1]) % 2nd order BOX-JENKINS modelbj2 =
Discrete-time BJ model: y(t) = [B(z)/F(z)]u(t) + [C(z)/D(z)]e(t)
B(z) = 0.9922 z^-1 + 0.4701 z^-2
C(z) = 1 - 0.6283 z^-1 - 0.1221 z^-2
D(z) = 1 - 1.221 z^-1 + 0.3798 z^-2
F(z) = 1 - 1.522 z^-1 + 0.7243 z^-2
Sample time: 0.25 seconds
Parameterization:
Polynomial orders: nb=2 nc=2 nd=2 nf=2 nk=1
Number of free coefficients: 8
Use "polydata", "getpvec", "getcov" for parameters and their uncertainties.
Status:
Estimated using BJ on time domain data "ze".
Fit to estimation data: 75.47% (prediction focus)
FPE: 0.9722, MSE: 0.8974
Сравнение предполагаемых моделей - симуляция и поведение предсказания
Теперь, когда мы оценили столько различных моделей, давайте сделаем некоторые сравнения модели. Это может быть сделано путем сравнения симулированных выходных параметров как прежде:
clf compare(zv,am2,md2,bj2,mx,mtf,m)
Мы можем также выдержать сравнение, как хорошо различные предполагаемые модели могут предсказать выход, скажем, 1 шаг вперед:
compare(zv,am2,md2,bj2,mx,mtf,m,1)
Остаточный анализ модели bj2 показывает, что ошибка предсказания лишена любой информации - это не автокоррелируется и также некоррелированое с входом. Это показывает, что дополнительные динамические элементы модели Box-Jenkins (C и полиномы D) смогли получить шумовую динамику.
resid(zv,bj2)
Сравнение функций частоты
Для того, чтобы выдержать сравнение, частота функционирует для сгенерированных моделей, мы используем снова bodeplot команда:
clf opt = bodeoptions; opt.PhaseMatching = 'on'; w = linspace(0.1,4*pi,100); bodeplot(GS,m,mx,mtf,md2,am2,bj2,w,opt); legend({'SPA','SS','ARX','TF','OE','ARMAX','BJ'},'Location','West');
Шумовые спектры предполагаемых моделей могут также анализироваться. Например, здесь мы сравниваем шумовые спектры ARMAX и модели Box-Jenkins с пространством состояний и модели Spectral Analysis. Для этого мы используем spectrum команда:
spectrum(GS,m,mx,am2,bj2,w) legend('SPA','SS','ARX','ARMAX','BJ');
Сравнение предполагаемых моделей с истинной системой
Здесь мы подтверждаем предполагаемые модели против истинной системы m0 то, что мы раньше симулировали входные и выходные данные. Давайте сравним функции частоты модели ARMAX с истинной системой.
bode(m0,am2,w)
Ответы практически совпадают. Шумовые спектры могут также быть сравнены.
spectrum(m0,am2,w)
Давайте также исследуем диаграмму нулей и полюсов:
h = iopzplot(m0,am2);
Это видно, что полюса и нули истинной (синей) системы и (зеленая) модель ARMAX очень близки.
Мы можем также оценить неопределенность в нулях и полюсах. Чтобы построить области доверия вокруг предполагаемых полюсов и нулей, соответствующих 3 стандартным отклонениям, используйте showConfidence или включите "характеристику" области Доверия от контекста графика (щелкают правой кнопкой) по меню.
showConfidence(h,3)
%Мы видим, что истинные, синие нули и полюса хорошо в зеленых областях неопределенности.