Этот пример показывает несколько методов идентификации, доступных в 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 обозначает оператор сдвига так, что A (q) y (t) действительно короче для y (t) - 1,5 y (t-1) + 0,7 y (t-2). Отображение m0 модели показывает, что это модель ARMAX.
Эта конкретная структура модели известна как модель ARMAX, где AR (Autoregressive) относится к A-полиному, MA (Скользящее среднее значение) к шуму C-полиному и X к входу «eXtra» B (q) u (t).
С точки зрения общих передаточных функций G и H, эта модель соответствует параметризации:
G(q) = B(q)/A(q) и H(q) = C(q)/A(q), с общими знаменателями
Генерируем входной сигнал u и моделируйте реакцию модели на эти входы. The idinput команда может использоваться, чтобы создать входной сигнал к модели и iddata команда упакует сигнал в подходящий формат. The 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 (s) и DEN (s) являются числителем и полиномами знаменателя, а 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 представлена как: A (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 структуры Output-Error (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 второго порядка и модель 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 с моделями пространства состояний и спектрального анализа. Для этого мы используем 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)
%Мы видим, что истинные, синие нули и полюсы находятся хорошо внутри зеленых областей неопределенности.