В этом примере показано, как выполнить и устранить неполадки при идентификации системы SISO с использованием данных частотного отклика (FRD). Способы, описанные здесь, также могут быть применены к моделям MIMO и данным ввода-вывода частотной области.
При использовании tfest команда для оценки модели передаточной функции SISO из данных частотно-отклика, алгоритм оценки минимизирует следующую функцию потерь наименьших квадратов (стоимость):
(ωk)) |2
Здесь W - определяемый частотно-зависимый вес, G - передаточная функция, которая должна быть оценена, f - измеренные частотно-характеристические данные, а - частота. Nf - это число частот, на которых доступны данные. (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 рад/с и резонансные пики после этого. Поэтому модель не соответствует данным.
Алгоритм минимизирует величину ошибки в квадрате | w)) | 2 в функции потерь. Постройте график этой величины как функции частоты. Этот график ошибок дает представление о том, какие точки данных вносят наибольший вклад в функцию потерь, и также вероятные ограничивающие факторы во время оценки. Сюжет может помочь вам определить, почему существуют ложные или некаптурированные динамики.
Поскольку не указан частотно-зависимый вес, ) равно 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 рад/с, большинство вариаций данных связано с шумом. В данных отсутствуют значимые системные режимы (впадины или пики). Для решения этой проблемы можно использовать фильтр скользящего среднего по данным для сглаживания измеренного отклика.
Алгоритм игнорирует долину около 60 рад/с и три слегка затухающих резонансных пика, которые следуют за ней. Эти функции мало способствуют функции потерь, поскольку ошибка подгонки мала на этих частотах. Для решения этой проблемы можно указать зависящие от частоты весовые коэффициенты, чтобы алгоритм уделял больше внимания этим частотам.
Для улучшения расчетного качества модели выполните предварительную обработку данных. Для этого необходимо обрезать малошумящие части данных ниже 1 рад/с и выше 2е4 рад/с, которые не интересны. Затем используется фильтр скользящего среднего для сглаживания данных в низкочастотной области высокого значения ниже 40 рад/с. На этих частотах данные имеют низкое отношение сигнал/шум, но имеют динамику, которую вы заинтересованы в захвате. Не применяйте фильтр на частотах выше 40 рад/с, чтобы избежать сглаживания данных там, где вы видите долину и три следующих за ней пика.
Создание копии оригинала idfrd объект данных.
GfrdProcessed = Gfrd;
Усечение данных ниже 1 рад/с и выше 2е4 рад/с.
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 рад/с, где существует ложная динамика. Этот низкий вес и сглаживание, примененное ранее к этим данным, уменьшают вероятность ложных пиков в оценочной реакции модели в этой области.
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/| ) | (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)');
Новая модель успешно фиксирует всю интересующую системную динамику.
График взвешенных ошибок можно использовать для дальнейшего устранения неполадок, если первоначальный выбор весов не дает удовлетворительного результата.