Оцените модель без предварительной обработки и фильтрации

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

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)')

Из данных, модели и графиков ошибок можно увидеть, что:

Предварительная обработка данных

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

Сделайте копию оригинала 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 рад/с, где вы хотите захватить динамику, но величина данных отклика низкая.

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)');

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

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