В этом примере показано несколько методов идентификации, доступных в 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 (авторегрессия) относится к A-многочлену, MA (скользящая средняя) - к C-многочлену шума и X - к входу «eXtra» 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, U и E являются преобразованиями Лапласа выходных, входных и ошибочных сигналов (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),
или в длинных руках,
+ 1) + 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 второго порядка и модель Бокса-Дженкинса второго порядка. Модель ARMAX имеет те же шумовые характеристики, что и моделируемая модель m0 и модель Бокса-Дженкинса позволяет более общее описание шума.
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 показывает, что ошибка предсказания лишена какой-либо информации - она не автокоррелирована и также не коррелирована с входными данными. Это показывает, что дополнительные динамические элементы модели Бокса-Дженкинса (многочлены 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 с моделями state-space и 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)
%Мы видим, что истинные, синие нули и полюса находятся внутри зеленых областей неопределенности.