Идентификация частотного диапазона: оценка моделей с использованием данных частотного диапазона

Этот пример показывает, как оценить модели, используя данные частотного диапазона. Оценка и валидация моделей, использующих данные частотного диапазона, работают так же, как и с данными временного интервала. Это обеспечивает большую гибкость в оценке и анализе моделей с использованием временного и частотного диапазона, а также спектральных (FRF) данных. Вы можете одновременно оценить модели с помощью данных в обеих областях, сравнить и объединить эти модели. Модель, оцененная с использованием временного интервала данных, может быть подтверждена с помощью спектральных данных или наоборот.

Данные частотного диапазона не могут использоваться для оценки или валидации нелинейных моделей.

Введение

Экспериментальные данные частотного диапазона распространены во многих приложениях. Возможно, данные были собраны как данные частотной характеристики (частотные функции: FRF) из процесса с помощью частотного анализатора. Может также быть, что более практично работать с преобразованиями Фурье входного и выходного сигнала (FFT данных временной области), например, для обработки периодических или ограниченных по полосам данных. (Ограниченный полосой непрерывный сигнал времени не имеет частотных составляющих выше частоты Найквиста). В System Identification Toolbox данные ввода-вывода частотного диапазона представлены так же, как и данные временной области, то есть использование iddata объекты. Свойство 'Domain' объекта должно быть установлено на 'Frequency'. Данные частотной характеристики представлены как комплексные векторы или как векторы амплитуда/фаза как функция от частоты. Объекты IDFRD в тулбоксе используются для инкапсуляции FRF, где пользователь задает комплексные данные отклика и вектор частоты. Такие объекты IDDATA или IDFRD (а также объекты FRD Control System Toolbox) могут использоваться бесшовно с любой стандартной программой оценки (такой как procest, tfest и т.д.).

Проверка Частотного диапазона данных

Начнем с загрузки некоторых данных частотного диапазона:

load demofr

Этот MAT-файл содержит данные частотной характеристики на частотах W, с амплитудной характеристикой AMP и фазового отклика PHA. Давайте сначала посмотрим на данные:

subplot(211), loglog(W,AMP),title('Amplitude Response')
subplot(212), semilogx(W,PHA),title('Phase Response')

Теперь эти экспериментальные данные будут сохранены как объект IDFRD. Первое преобразование амплитуды и фазы в комплексный оцененный ответ:

zfr = AMP.*exp(1i*PHA*pi/180);
Ts = 0.1;
gfr = idfrd(zfr,W,Ts);

Ts - шаг расчета базовых данных. Если данные соответствуют непрерывному времени, для примера, поскольку вход был ограничен полосой, используйте Ts = 0.

Примечание. Если у вас есть Control System Toolbox™, можно использовать объект FRD вместо объекта IDFRD. IDFRD имеет опции для получения дополнительной информации, такие как спектры нарушений порядка и меры неопределенности, которые недоступны в объектах FRD.

Объект IDFRD gfr теперь содержит данные, и они могут быть построены и проанализированы по-разному. Чтобы просмотреть данные, мы можем использовать plot или bode:

clf
bode(gfr), legend('gfr')

Оценка моделей с использованием данных частотной характеристики (FRF)

Для оценки моделей теперь можно использовать gfr как набор данных со всеми командами тулбокса прозрачным способом. Единственным ограничением является то, что модели шума не могут быть построены. Это означает, что для полиномиальных моделей применяется только OE (модели выходной ошибки), а для моделей пространства состояний необходимо исправить K = 0.

m1 = oe(gfr,[2 2 1]) % Discrete-time Output error (transfer function) model
ms = ssest(gfr) % Continuous-time state-space model with default choice of order
mproc = procest(gfr,'P2UDZ') % 2nd-order, continuous-time model with underdamped poles
compare(gfr,m1,ms,mproc)
L = findobj(gcf,'type','legend');
L.Location = 'southwest'; % move legend to non-overlapping location
m1 =
Discrete-time OE model: y(t) = [B(z)/F(z)]u(t) + e(t)
  B(z) = 0.9986 z^-1 + 0.4968 z^-2                   
                                                     
  F(z) = 1 - 1.499 z^-1 + 0.6998 z^-2                
                                                     
Sample time: 0.1 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 frequency response data "gfr".
Fit to estimation data: 88.04%                      
FPE: 0.2512, MSE: 0.2492                            

ms =
  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  -1.785   3.097
   x2  -6.835  -1.785
 
  B = 
          u1
   x1  -4.15
   x2  27.17
 
  C = 
           x1      x2
   y1    1.97  0.3947
 
  D = 
       u1
   y1   0
 
  K = 
       y1
   x1   0
   x2   0
 
Parameterization:
   FREE form (all coefficients in A, B, C free).
   Feedthrough: none
   Disturbance component: none
   Number of free coefficients: 8
   Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                                
Estimated using SSEST on frequency response data "gfr".
Fit to estimation data: 88.04%                         
FPE: 0.2512, MSE: 0.2492                               

mproc =
Process model with transfer function:            
                     1+Tz*s                      
  G(s) = Kp * ---------------------- * exp(-Td*s)
              1+2*Zeta*Tw*s+(Tw*s)^2             
                                                 
         Kp = 7.4619                             
         Tw = 0.20245                            
       Zeta = 0.36242                            
         Td = 0                                  
         Tz = 0.013617                           
                                                 
Parameterization:
    {'P2DUZ'}
   Number of free coefficients: 5
   Use "getpvec", "getcov" for parameters and their uncertainties.

Status:                                                  
Estimated using PROCEST on frequency response data "gfr".
Fit to estimation data: 88.03%                           
FPE: 0.2517, MSE: 0.2492                                 

Как показано выше, различные типы линейных моделей могут быть оценены как в непрерывной, так и в дискретное время областях, используя спектральные данные. Эти модели могут быть проверены с помощью данных во временной области. Набор данных ввода-вывода временной области ztime, для примера, собирается из той же системы и может использоваться для валидации m1, ms и mproc:

compare(ztime,m1,ms,mproc) %validation in a different domain

Мы также можем рассмотреть невязки, чтобы подтвердить качество модели с помощью данных валидации ztime. Как наблюдалось, невязки почти белые:

resid(ztime,mproc) % Residuals plot

Данные конденсации с использованием SPAFDR

Важной причиной для работы с данными частотной характеристики является то, что легко конденсировать информацию с небольшими потерями. Команда SPAFDR позволяет вам вычислять сглаженные данные отклика на ограниченных частотах, например, с логарифмическим интервалом. Вот пример, где gfr данные сконденсированы до 100 логарифмически разнесенных значений частоты. При подобном методе также исходные данные о временной интервал могут быть конденсированы:

sgfr = spafdr(gfr) % spectral estimation with frequency-dependent resolution
sz = spafdr(ztime); % spectral estimation using time-domain data
clf
bode(gfr,sgfr,sz)
axis([pi/100 10*pi, -272 105])
legend('gfr (raw data)','sgfr','sz','location','southwest')
sgfr =
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 100 frequency points, ranging from 0.03142 rad/s to 31.42 rad/s.
 
Sample time: 0.1 seconds
Output channels: 'y1'
Input channels: 'u1'
Status:                                                 
Estimated using SPAFDR on frequency response data "gfr".

Эти Диаграммы Боде показывают, что информация в сглаженных данных была хорошо обеспечена. Теперь эти записи данных со 100 точками могут очень хорошо использоваться для оценки модели. Для примера:

msm = oe(sgfr,[2 2 1]);
compare(ztime,msm,m1) % msm has the same accuracy as M1 (based on 1000 points)

Оценка с использованием данных ввода-вывода частотного диапазона

Может быть, что измерения доступны как преобразования Фурье входных и выходных входов. Такие данные частотного диапазона от системы даны как сигналы Y и U. На логарифмических графиках они выглядят как

Wfd = (0:500)'*10*pi/500;
subplot(211),loglog(Wfd,abs(Y)),title('The amplitude of the output')
subplot(212),loglog(Wfd,abs(U)),title('The amplitude of the input')

Данные частотной характеристики по существу являются отношением между Y и U. Чтобы собрать частотным диапазоном данные как объект IDDATA, выполните следующие действия:

ZFD = iddata(Y, U, 'Ts', 0.1, 'Frequency', Wfd)
ZFD =

Frequency domain data set with responses at 501 frequencies.
Frequency range: 0 to 31.416 rad/seconds
Sample time: 0.1 seconds                                                                              
                                                                                                      
Outputs      Unit (if specified)                                                                      
   y1                                                                                                 
                                                                                                      
Inputs       Unit (if specified)                                                                      
   u1                                                                                                 
                                                                                                      

Теперь снова набор данных частотного диапазона ZFD могут использоваться в качестве данных во всех стандартных программах оценки, так же как данные временному интервалу и данные частотной характеристики:

mf = ssest(ZFD)   % SSEST picks best order in 1:10 range when called this way
mfr = ssregest(ZFD) % an alternative regularized reduction based state-space estimator
clf
compare(ztime,mf,mfr,m1)
mf =
  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   -1.78   3.095
   x2  -6.812   -1.78
 
  B = 
          u1
   x1   1.32
   x2  28.61
 
  C = 
              x1         x2
   y1          2  6.418e-08
 
  D = 
       u1
   y1   0
 
  K = 
       y1
   x1   0
   x2   0
 
Parameterization:
   FREE form (all coefficients in A, B, C free).
   Feedthrough: none
   Disturbance component: none
   Number of free coefficients: 8
   Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                              
Estimated using SSEST on frequency domain data "ZFD".
Fit to estimation data: 97.21%                       
FPE: 0.04288, MSE: 0.04186                           

mfr =
  Discrete-time identified state-space model:
    x(t+Ts) = A x(t) + B u(t) + K e(t)
       y(t) = C x(t) + D u(t) + e(t)
 
  A = 
               x1         x2         x3         x4         x5         x6
   x1       0.701   -0.05307    -0.4345   0.006642   -0.08085    -0.1158
   x2     -0.4539   -0.09623    -0.3629    -0.2113    -0.2219     0.3536
   x3      0.1719    -0.3996     0.4385    0.01558     0.1768     0.2467
   x4      0.1821    -0.3465    -0.3292    -0.1357     0.2578    -0.2483
   x5    -0.09939     -0.338     0.1236     0.2537     -0.387    0.05591
   x6    -0.06004      0.226    0.04117     0.6891    -0.0873     0.2818
   x7     -0.1056     0.1859    0.04223      0.629     0.1968    -0.7077
   x8     0.05337     0.1948    0.06355    0.09052     0.4216     0.3997
   x9     0.01696    0.05961    0.04891   -0.01251    0.03521     0.3876
   x10    0.01727     0.1232    0.03586     0.1187     0.1738    0.05051
 
               x7         x8         x9        x10
   x1     -0.3087   0.007547   -0.02011     0.1469
   x2    -0.01728   -0.04967    -0.1144   -0.03883
   x3     0.07107    -0.2977      0.129    0.06179
   x4    -0.08461   -0.03541   -0.06711    -0.1759
   x5     -0.5324     0.1778     0.1114  2.119e-05
   x6       0.155    -0.5047     -0.285    -0.3976
   x7     -0.2406     0.5628    0.09159     0.4845
   x8     -0.5674    -0.3337     -0.105    -0.1995
   x9      0.2024     0.4718   -0.01861    -0.5838
   x10   -0.01243    -0.2968     0.6808    -0.7726
 
  B = 
              u1
   x1   -0.09482
   x2     0.7665
   x3     -1.036
   x4     -0.162
   x5    -0.2926
   x6     0.0805
   x7    -0.1277
   x8   -0.02441
   x9    -0.0288
   x10  -0.03776
 
  C = 
            x1       x2       x3       x4       x5       x6       x7
   y1    1.956  -0.4539   -1.805    1.356   0.3662    1.691   0.8489
 
            x8       x9      x10
   y1    0.148   0.5203   0.4013
 
  D = 
       u1
   y1   0
 
  K = 
                y1
   x1       0.1063
   x2    -0.007395
   x3       0.0426
   x4    -0.004966
   x5     -0.01505
   x6    -0.005483
   x7   -0.0004218
   x8    -0.001044
   x9    -0.003066
   x10   -0.002273
 
Sample time: 0.1 seconds
  
Parameterization:
   FREE form (all coefficients in A, B, C free).
   Feedthrough: none
   Disturbance component: estimate
   Number of free coefficients: 130
   Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                                 
Estimated using SSREGEST on frequency domain data "ZFD".
Fit to estimation data: 76.61% (prediction focus)       
FPE: 3.448, MSE: 2.938                                  

Преобразования между представлениями данных (время - частота)

Наборы входно-выходных данных временных и частотных диапазонов могут быть преобразованы в любую область с помощью FFT и IFFT. Эти команды адаптируются к объектам IDDATA:

dataf = fft(ztime)
datat = ifft(dataf)
dataf =

Frequency domain data set with responses at 501 frequencies.
Frequency range: 0 to 31.416 rad/seconds
Sample time: 0.1 seconds                                                                              
                                                                                                      
Outputs      Unit (if specified)                                                                      
   y1                                                                                                 
                                                                                                      
Inputs       Unit (if specified)                                                                      
   u1                                                                                                 
                                                                                                      

datat =

Time domain data set with 1000 samples.
Sample time: 0.1 seconds                
                                        
Outputs      Unit (if specified)        
   y1                                   
                                        
Inputs       Unit (if specified)        
   u1                                   
                                        

Входно-выходные данные временных и частотных диапазонов могут быть преобразованы в данные частотной характеристики SPAFDR, SPA и ETFE:

g1 = spafdr(ztime)
g2 = spafdr(ZFD);
clf;
bode(g1,g2)
g1 =
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 100 frequency points, ranging from 0.06283 rad/s to 31.42 rad/s.
 
Sample time: 0.1 seconds
Output channels: 'y1'
Input channels: 'u1'
Status:                                            
Estimated using SPAFDR on time domain data "ztime".

Данные частотной характеристики могут также быть преобразованы в более сглаженные данные (меньше разрешения и меньше данных) SPAFDR и SPA;

g3 = spafdr(gfr);

Данные частотной характеристики могут быть преобразованы во входно-выходные сигналы частотного диапазона командой IDDATA:

gfd = iddata(g3)
plot(gfd)
gfd =

Frequency domain data set with responses at 100 frequencies.
Frequency range: 0.031416 to 31.416 rad/seconds
Sample time: 0.1 seconds                                                                                     
                                                                                                             
Outputs      Unit (if specified)                                                                             
   y1                                                                                                        
                                                                                                             
Inputs       Unit (if specified)                                                                             
   u1                                                                                                        
                                                                                                             

Использование данных частотного диапазона в непрерывном времени для оценки моделей в непрерывном времени

Временные интервалы данные могут, естественно, храниться и рассматриваться как дискретное время, выборочные данные. Данные частотного диапазона имеют то преимущество, что непрерывные временные данные могут быть представлены правильно. Предположим, что базовые непрерывные сигналы времени не имеют информации о частоте выше частоты Найквиста, например, потому что они дискретизированы быстро, или вход не имеет частотной составляющей выше частоты Найквиста, и что данные были собраны из установившегося эксперимента. Тогда дискретные преобразования Фурье (DFT) данных также являются преобразованиями Фурье непрерывных временных сигналов на выбранных частотах. Поэтому они могут использоваться, чтобы непосредственно соответствовать непрерывным временным моделям.

Это проиллюстрировано на следующем примере.

Рассмотрим непрерывную систему времени:

$$ G(s) = \frac{1}{s^2+s+1} $$

m0 = idpoly(1,1,1,1,[1 1 1],'ts',0)
m0 =
Continuous-time OE model: y(t) = [B(s)/F(s)]u(t) + e(t)
  B(s) = 1                                             
                                                       
  F(s) = s^2 + s + 1                                   
                                                       
Parameterization:
   Polynomial orders:   nb=1   nf=2   nk=0
   Number of free coefficients: 3
   Use "polydata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                                         
Created by direct construction or transformation. Not estimated.

Загрузка данных, полученных в результате статической симуляции этой системы с помощью периодических входов. Собранные данные были преобразованы в частотный диапазон и сохранены в файле CTFDData.mat.

load CTFDData.mat dataf % load continuous-time frequency-domain data.

Посмотрите на данные:

plot(dataf)
set(gca,'XLim',[0.1 10])

Использование dataf для оценки будет по умолчанию давать непрерывные модели времени: Пространство состояний:

m4 = ssest(dataf,2); %Second order continuous-time model

Для полиномиальной модели с nb = 2 коэффициент и nf = 2 числителя использование расчетных коэффициентов знаменателя:

nb = 2;
nf = 2;
m5 = oe(dataf,[nb nf])
m5 =
Continuous-time OE model: y(t) = [B(s)/F(s)]u(t) + e(t)
  B(s) = -0.01831 s + 1.009                            
                                                       
  F(s) = s^2 + 1.001 s + 0.997                         
                                                       
Parameterization:
   Polynomial orders:   nb=2   nf=2   nk=0
   Number of free coefficients: 4
   Use "polydata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                             
Estimated using OE on frequency domain data "dataf".
Fit to estimation data: 70.15%                      
FPE: 0.00491, MSE: 0.00468                          

Сравнение переходных характеристик с неопределенностью истинной системы m0 и модели m4 и m5. Доверительные интервалы показаны с закрашенными фигурами на графике.

clf
h = stepplot(m0,m4,m5);
showConfidence(h,1)
legend('show','location','southeast')

Хотя в этом случае не было необходимости, обычно рекомендуется фокусировать подгонку на ограниченной полосе частот (нижние частоты фильтруют данные) при оценке с использованием непрерывных временных данных. Система имеет пропускную способность около 3 рад/с и возбуждалась синусоидами до 6,2 рад/с. Разумная область значений для фокусировки подгонки равна [0 7] рад/с:

m6 = ssest(dataf,2,ssestOptions('WeightingFilter',[0 7])) % state space model
m6 =
  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.5011    1.001
   x2  -0.7446  -0.5011
 
  B = 
             u1
   x1  -0.01706
   x2     1.016
 
  C = 
               x1          x2
   y1       1.001  -0.0005347
 
  D = 
       u1
   y1   0
 
  K = 
       y1
   x1   0
   x2   0
 
Parameterization:
   FREE form (all coefficients in A, B, C free).
   Feedthrough: none
   Disturbance component: none
   Number of free coefficients: 8
   Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                                
Estimated using SSEST on frequency domain data "dataf".
Fit to estimation data: 87.03% (data prefiltered)      
FPE: 0.004832, MSE: 0.003631                           
m7 = oe(dataf,[1 2],oeOptions('WeightingFilter',[0 7])) % polynomial model of Output Error structure
m7 =
Continuous-time OE model: y(t) = [B(s)/F(s)]u(t) + e(t)
  B(s) = 0.9861                                        
                                                       
  F(s) = s^2 + 0.9498 s + 0.9704                       
                                                       
Parameterization:
   Polynomial orders:   nb=1   nf=2   nk=0
   Number of free coefficients: 3
   Use "polydata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                             
Estimated using OE on frequency domain data "dataf".
Fit to estimation data: 86.81% (data prefiltered)   
FPE: 0.004902, MSE: 0.003752                        
opt = procestOptions('SearchMethod','lsqnonlin',...
   'WeightingFilter',[0 7]); % Requires Optimization Toolbox(TM)
m8 = procest(dataf,'P2UZ',opt)  % process model with underdamped poles
m8 =
Process model with transfer function:
                     1+Tz*s          
  G(s) = Kp * ---------------------- 
              1+2*Zeta*Tw*s+(Tw*s)^2 
                                     
         Kp = 1.0124                 
         Tw = 1.0019                 
       Zeta = 0.5021                 
         Tz = -0.017474              
                                     
Parameterization:
    {'P2UZ'}
   Number of free coefficients: 4
   Use "getpvec", "getcov" for parameters and their uncertainties.

Status:                                                  
Estimated using PROCEST on frequency domain data "dataf".
Fit to estimation data: 87.03% (data prefiltered)        
FPE: 0.004832, MSE: 0.003631                             
opt = tfestOptions('SearchMethod','lsqnonlin',...
   'WeightingFilter',[0 7]); % Requires Optimization Toolbox
m9 = tfest(dataf,2,opt) % transfer function with 2 poles
m9 =
 
  From input "u1" to output "y1":
  -0.01662 s + 1.007
  ------------------
   s^2 + s + 0.995
 
Continuous-time identified transfer function.

Parameterization:
   Number of poles: 2   Number of zeros: 1
   Number of free coefficients: 4
   Use "tfdata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                                
Estimated using TFEST on frequency domain data "dataf".
Fit to estimation data: 87.03% (data prefiltered)      
FPE: 0.00491, MSE: 0.003629                            
h = stepplot(m0,m6,m7,m8,m9);
showConfidence(h,1)
legend('show')

Заключения

Мы видели, как время, частота и спектральные данные могут беспрепятственно использоваться, чтобы оценить множество линейных моделей как в непрерывной, так и в дискретной временной областях. Модели могут быть проверены и сравнены в областях, отличных от тех, в которых они были оценены. Форматы данных (время, частота и спектр) являются взаимоконвертируемыми, используя такие методы, как fft, ifft, spafdr и spa. Кроме того, прямая оценка в непрерывном времени достижима при помощи tfest, ssest и procest стандартные программы оценки. Плавное использование данных в любой области для оценки и анализа является важной функцией System Identification Toolbox.

См. также

| | |

Похожие темы

Для просмотра документации необходимо авторизоваться на сайте