Диагностируйте идентификацию частотного диапазона моделей передаточной функции

Этот пример показывает, как выполнить и диагностировать идентификацию системы SISO с помощью данных частотной характеристики (FRD). Методы, объясненные здесь, могут также быть применены к моделям MIMO и данным ввода - вывода частотного диапазона.

Когда вы используете команду tfest, чтобы оценить модель передаточной функции SISO от данных частотной характеристики, алгоритм оценки минимизирует следующую потерю наименьших квадратов (стоимость) функция:

$\underset{G (ω)}{m i n i m i z e} \sum_{k = 1}^{N_{f}} {| W (ω_{k}) (G (ω_{k}) - f (ω_{k})) |}^{2}$

Здесь W является зависимым частотой весом, который вы задаете, G является передаточной функцией, которая должна быть оценена, f является измеренными данными частотной характеристики, и $w$ частота. Nf является частотами номера, на которых данные доступны. $G (w) - f (w)$ ошибка частотной характеристики.

В этом примере вы сначала оцениваете модель, не предварительно обрабатывая данные или с помощью опций оценки, чтобы задать фильтр взвешивания. Вы затем применяете эти методы поиска и устранения неисправностей, чтобы улучшить образцовую оценку.

Оцените модель, не предварительно обрабатывая и фильтруя

Загрузите измеренные непрерывно-разовые данные о частотной характеристике.

load troubleshooting_example_data Gfrd;

Gfrd является объектом idfrd, который хранит данные.

Оцените первоначальную модель передаточной функции с 11 полюсами и 10 нулями при помощи данных об оценке.

Gfit = tfest(Gfrd,11,10);

Постройте значение частотной характеристики предполагаемой модели и измеренных данных частотной характеристики.

bodemag(Gfrd,Gfit);
ylim([-100 0])
legend('Measured','Estimated')

Предполагаемая модель содержит побочную динамику. Алгоритм оценки пропускает долину на уровне 60 рад/с и резонирующий peaks после этого. Поэтому модель не является подходящим вариантом к данным.

Алгоритм минимизирует значение квадратичной невязки, $| W (w) (G (w) - f (w)) |^{2}$ , в функции потерь. Постройте это количество как функцию частоты. Эта диаграмма погрешностей обеспечивает представление, которого точки данных способствуют больше всего функции потерь, и так являются вероятными ограничивающими факторами во время оценки. График может помочь вам идентифицировать, почему существуют побочные или неполученные движущие силы.

Поскольку вы не задали зависимый частотой вес, $W (w)$ 1.

w = Gfrd.Frequency;
r1 = squeeze(freqresp(Gfit,w));
r2 = squeeze(freqresp(Gfrd,w));
fitError = r1-r2;
semilogx(w,abs(fitError).^2)
title('Weighted Estimation Error');
xlabel('Frequency (rad/s)');
ylabel('Magnitude (abs)')

От данных, модели и диаграмм погрешностей вы видите что:

Самые большие подходящие ошибки ниже 10 рад/с.
Алгоритм фокусируется на подборе кривой шумным высоким точкам данных значения ниже 10 рад/с, которые имеют большой вклад в функцию потерь оптимизации. В результате алгоритм присваивает побочные полюса и нули в эту область данных. Чтобы решить эту проблему, можно предварительно обработать данные, чтобы улучшить отношение сигнал-шум в этой области. Можно также использовать зависимые частотой веса, чтобы заставить алгоритм поместить меньше особого внимания на эту область.
Ниже приблизительно 40 рад/с большинство изменений в данных происходит из-за шума. Нет никаких значительных системных режимов (долины или peaks) в данных. Чтобы решить эту проблему, можно использовать фильтр скользящего среднего значения по данным, чтобы сглаживать измеренный ответ.
Алгоритм игнорирует долину приблизительно 60 рад/с и три слегка ослабленных резонирующих peaks, которые следуют за ним. Эти функции способствуют мало функции потерь, потому что подходящая ошибка является небольшой на этих частотах. Чтобы решить эту проблему, можно задать зависимые частотой веса, чтобы заставить алгоритм уделить больше внимания этим частотам.

Предварительно обработайте данные

Чтобы улучшить предполагаемое образцовое качество, предварительно обработайте данные. Для этого вы обрезаете низкие фрагменты сигнала к шуму данных ниже 1 рад/с и выше 2e4 рад/с, которые не интересны. Затем вы используете фильтр скользящего среднего значения, чтобы сглаживать данные в низкочастотной области высокого значения ниже 40 рад/с. На этих частотах данные имеют низкое отношение сигнал-шум, но имеют динамику, что вы интересуетесь получением. Не применяйте фильтр на частотах выше 40 рад/с, чтобы не сглаживать данные, где вы видите долину и три peaks, которые следуют за ним.

Сделайте копию исходного объекта данных idfrd.

GfrdProcessed = Gfrd;

Обрежьте данные ниже 1 рад/с и выше 2e4 рад/с.

GfrdProcessed = fselect(GfrdProcessed,1,2e4);

Примените фильтр скользящего среднего значения в центре с тремя точками, чтобы сгладить данные частотной характеристики ниже 40 рад/с, которые содержат побочную динамику. Данные об ответе хранятся в свойстве ResponseData объекта.

w = GfrdProcessed.Frequency;
f = squeeze(GfrdProcessed.ResponseData);
idx = w<40;
f(idx) = movmean(f(idx),3);

Здесь f(idx) является данными частотной характеристики на частотах меньше чем 40 рад/с.

Поместите отфильтрованные данные назад в исходный объект данных.

GfrdProcessed.ResponseData = f;

Отобразите исходные и предварительно обработанные данные на графике.

bodemag(Gfrd,GfrdProcessed);
ylim([-100 0]);
legend('Original data','Preprocessed data');

График показывает, что все желаемые движущие силы неповреждены после предварительной обработки.

Задайте фильтр взвешивания

Используйте низкий вес для низкочастотной области менее чем 10 рад/с, где побочные движущие силы существуют. Этот низкий вес и сглаживание, примененное ранее к этим данным, уменьшают шанс побочного peaks в предполагаемом образцовом ответе в этой области.

Weight = ones(size(f));
idx = w<10;
Weight(idx) = Weight(idx)/10;

Используйте высокий вес для данных в частотном диапазоне 40-6e3 rad/s, где вы хотите получить динамику, но значение данных об ответе является низким.

idx = w>40 & w<6e3;
Weight(idx) = Weight(idx)*30;

Задайте веса в опции WeightingFilter набора опции оценки.

tfestOpt = tfestOptions('WeightingFilter',Weight);

Обратите внимание на то, что Weight является пользовательским фильтром взвешивания. Можно также задать WeightingFilter как 'inv' или 'invsqrt' для данных частотной характеристики. Эти опции задают вес как $1 / | f (w) |$ и $1 / \sqrt{| f (w) |}$ , соответственно. Эти опции позволяют вам быстро протестировать эффект использования более высокого веса для низких областей значения данных. 'invsqrt' обычно является хорошим начальным выбором. Если эти веса не приводят к хорошим результатам оценки, можно обеспечить пользовательские веса как показано в этом примере.

Оцените модель Используя предварительно обработанные и фильтрованные данные

Оцените модель передаточной функции с 11 полюсами и 10 нулями с помощью предварительно обработанных данных и заданного фильтра взвешивания.

GfitNew = tfest(GfrdProcessed,11,10,tfestOpt);

Отобразите на графике исходные данные, начальный образцовый ответ и новый образцовый ответ.

bodemag(Gfrd,Gfit,GfitNew);
ylim([-100 0]);
legend('Original Data','Original Model','New Model');

Постройте ошибку оценки. Вычислите ошибку оценки включением фильтра взвешивания Weight, который вы использовали для оценки GfitNew.

w = GfrdProcessed.Frequency;
r1 = squeeze(freqresp(GfitNew,w));
r2 = squeeze(freqresp(GfrdProcessed,w));
fitErrorNew = Weight.*(r1-r2);
semilogx(w,abs(fitErrorNew).^2)
title('Weighted Estimation Error');
xlabel('Frequency (rad/s)');
ylabel('Magnitude (abs)');

Новая модель успешно получает всю системную динамику интереса.

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

Смотрите также

tfest | tfestOptions

Документация