Сравнение различных образцовых методов идентификации
Этот пример показывает несколько методов идентификации, доступных в 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".
Чтобы построить частотную характеристику, это обычно, чтобы использовать Диаграмму Боде с командой bodeplot или bode:
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)
%Мы видим, что истинные, синие нули и полюса хорошо в зеленых областях неуверенности.