В этом примере показано, как оценивать модели с использованием данных частотной области. Оценка и валидация моделей с использованием данных частотной области работают так же, как и с данными временной области. Это обеспечивает большую гибкость в оценке и анализе моделей с использованием данных временной и частотной областей, а также спектральных (FRF) данных. Можно одновременно оценивать модели, используя данные в обеих областях, сравнивать и объединять эти модели. Модель, оцененная с использованием данных временной области, может быть проверена с использованием спектральных данных или наоборот.
Данные частотной области не могут использоваться для оценки или проверки нелинейных моделей.
Экспериментальные данные частотной области распространены во многих приложениях. Возможно, данные собирались как данные частотной характеристики (частотные функции: FRF) из процесса с помощью частотного анализатора. Также может оказаться более практичным работать с преобразованиями Фурье входных и выходных данных (БПФ данных временной области), например, обрабатывать периодические данные или данные с ограниченной полосой. (Непрерывный сигнал времени с ограниченной полосой частот не имеет частотных составляющих выше частоты Найквиста). В панели инструментов идентификации системы данные ввода-вывода в частотной области представлены так же, как и данные во временной области, т.е. с использованием iddata объекты. Свойство Domain объекта должно иметь значение Frequency. Данные частотной характеристики представлены как комплексные векторы или как векторы величины/фазы как функция частоты. Объекты IDFRD на панели инструментов используются для инкапсуляции FRF, где пользователь задает сложные данные отклика и частотный вектор. Такие объекты IDDATA или IDFRD (а также объекты FRD панели инструментов системы управления) могут использоваться с любой подпрограммой оценки (например, 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.
Примечание.При наличии Toolbox™ системы управления вместо объекта IDFRD можно использовать объект FRD. IDFRD имеет варианты дополнительной информации, такие как спектры возмущений и меры неопределенности, которые недоступны в объектах FRD.
Объект IDFRD gfr теперь содержит данные, которые могут быть нанесены на график и проанализированы различными способами. Для просмотра данных мы можем использовать plot или bode:
clf
bode(gfr), legend('gfr')
Для оценки моделей теперь можно использовать 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 позволяет вычислять сглаженные данные отклика на ограниченных частотах, например, с логарифмическим интервалом. Вот пример, где 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
Наборы данных ввода-вывода во временной и частотной областях могут быть преобразованы в любую область с использованием БПФ и ОБПФ. Эти команды адаптированы к объектам 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) данных также являются преобразованиями Фурье непрерывных временных сигналов на выбранных частотах. Поэтому они могут использоваться для непосредственного соответствия непрерывным моделям времени.
Это будет проиллюстрировано следующим примером.
Рассмотрим систему непрерывного времени:
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 для оценки будет по умолчанию давать непрерывные модели времени: State-space:
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 процедуры оценки. Бесшовное использование данных в любой области для оценки и анализа является важной особенностью инструментария идентификации системы.