Описание системы

Линейная система с поврежденным шумом может быть описана следующим образом:

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 automatically

m =
  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 poles

mtf =
 
  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),

или длинными руками,

Эта модель линейна в коэффициентах. Коэффициенты ${\mathit{a}}_1^{}$, ${\mathit{a}}_2^{}$, ..., ${\mathit{b}}_1^{}$, ${\mathit{b}}_2^{}$, .. может быть эффективно вычислено с использованием метода оценки методом наименьших квадратов. Оценка модели 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 второго порядка и модель Box-Jenkins второго порядка. Модель ARMAX имеет те же шумовые характеристики, что и моделируемая модель m0 и модель Box-Jenkins позволяет более общее описание шума.

am2 = armax(ze,[2 2 2 1])  % 2nd order ARMAX model

am2 =
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 model

bj2 =
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)

%

Мы видим, что истинные, синие нули и полюсы находятся хорошо внутри зеленых областей неопределенности.