Обзор методов анализа моделей, систем и сигналов

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

Идентификация моделей

Рабочие файлы: [ABCD.vsm]

Идентификация компьютерных моделей
Та или иная алгоритмическая процедура, в результате выполнения которой моделирующая программа получает численные значения коэффициентов модели, структура которой и параметры были выбраны пользователем произвольно. Искомые коэффициенты могут быть представлены либо в форме ABCD-матриц, либо в форме коэффициентов полиномов числителя и знаменателя передаточной функции. Цель поиска коэффициентов состоит в предоставлении исходных данных для частотного, корневого и других видов анализа.

Аналитическая процедура приведения направленного графа сложна для программной реализации. По этой причине библиотеки анализа программ математического моделирования не отслеживают ни схему соединений блоков, определяемую пользователем, ни их параметры. Т.е. подобная информация для идентификации модели им не нужна (см. рис. 1).

gif-file, 2KB

Рис. 1

Но, какова бы ни была структура модели линейной системы, она может быть описана в матричной форме именуемой "пространством состояний":

$u′=Au+Bx$
$y=Cu+Dx$

Данная система уравнений для SISO-систем (модели с одним входом и выходом) имеет следующую масштабируемую топологию:

gif-file, 2KB

Легко увидеть, что в целях унификации вычислительных алгоритмов моделирующих программ можно объединить матрицы A, B, C и D, а так же матрицу $U$ с координатой $x$, и матрицу $U′$ с координатой $y$, как показано на рис. 2.

gif-file, 2KB

Рис. 2

Таким образом можно утверждать, что любую линейную модель можно преобразовать к универсальной блок-схеме изображенной на рис. 2. При этом в результате преобразования граф распространения сигнала в модели не претерпит изменений. Т.е. матричная блок-схема является универсальной формой отображения любых соединений блоков линейной модели, а не просто её эквивалентом.

gif-file, 2KB

Рис. 3

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

  1. Моделирующая программа для всех интеграторов модели (см. рис. 1) отключает вызов функции численного интегрирования, и, к освободившимся входам и выходам блоков 1/S, подключает библиотеку анализа. Дополнительными точками подключения библиотеки являются входная и выходная координаты модели. Библиотека анализа использует выходы блоков интеграторов для задания возмущений на граф модели, а их входы для фиксации отклика графа.
  2. В момент подключения (инициализации) библиотека анализа создает будущий (искомый) "образ модели" в виде обнуленной ABCD-матрицы.
  3. Далее, моделирующая программа инициирует процесс спец-симуляции модели длинной в $n+1$ шаг (на данном этапе функционирования, программы не используют графический интерфейс, т.е. скрывают свои действия от пользователя).
  4. На каждом шаге спец-симуляции, библиотека анализа поочередно устанавливает на выходе одного из интеграторов модели единицу (другие возмущения обнуляются), и фиксирует реакцию графа.
  5. Используя полученную совокупность отдельных реакций графа, библиотека анализа заполняет столбцы объединенной ABCD-матрицы в искомом "образе модели". Суть этой операции поясняет рис. 3.

Описанный сценарий легко поддается алгоритмизации и реализован в большинстве программ математического моделирования с поточной моделью управления. Следует отметить, что он может быть использован для идентификации мультичастотных дискретных, гибридных, а так же непрерывных систем со звеньями чистого запаздывания. В этом случае список блоков, к которым должна подключаться библиотека анализа расширяется регистрами задержки 1/Z и звеньями чистого запаздывания $e^{-τs}$. Но, составляя вектор задания, библиотека анализа должна хранить информацию о том, каким устройством формируется та или иная координата модели и его параметры ($T_ц$ для 1/Z и τ для $e^{-τs}$), что требуется для процедуры расчета частотных характеристик.

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

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

Билинейное преобразование

Рабочие файлы: [Аппроксиматоры]

Билинейное преобразование
Бинаправленная процедура пересчета коэффициентов линейной модели, используемая с целью нахождения соответствующего (дискретного или непрерывного) аналога, т.е. для трансформации моделей из $s$-домена в $z$-домен и обратно.

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

Если требуется перейти от непрерывного прототипа к дискретной модели, то процедура идентификации (в данном случае "перерассчитываемых") коэффициентов будет отличаться лишь тем, что библиотека анализа подключится не к выводам интеграторов модели, а к регистрам задержки, на которых эти (квазианалоговые) интеграторы реализованы. Следует уточнить, что в этой процедуре имеет смысл использовать лишь квазианалоговый интегратор первого порядка (на одном регистре), реализующий метод трапеций, поэтому библиотеки автоматически активируют этот метод интегрирования в настройках программ.

Если же требуется обратная трансформация модели, то процедура идентификации "перерассчитываемых" коэффициентов так же чуть модифицируется. Библиотека анализа замещает все регистры задержки модели их непрерывными аналогами – фазосдвигающими звеньями и подключается к интеграторам, на которых эти аппроксиматоры реализованы.

gif-file, 2KB

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

gif-file, 2KB или gif-file, 2KB .

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

gif-file, 2KB .

Разложение математического описания модели САР в степенной ряд

Рабочие файлы: [ABCD → с1с2с3]

Разложение математического описания модели в степенной ряд
Процедура разложения в степенной ряд системы дифференциальных уравнений САР решенной относительно требуемой координаты. Таким образом результат процедуры – первые коэффициенты разложения (коэффициенты ошибок, или коэффициенты отклика) – получается в результате трансформации либо ABCD-матрицы коэффициентов, либо коэффициентов числителя и знаменателя ПФ.

Известно, что например коэффициенты ошибок связывают коэффициенты разложения входного сигнала $g(t)$ в ряд Тейлора с соответствующими составляющими ошибки $x(t)$ в установившемся режиме движения:

$x(t)=c_0g(t)+c_1g′(t)/1!+c_2g′′(t)/2!+c_3 g′′′(t)/3!+…$

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

В общем виде математическое обоснование техники разложения математического описания модели САР в степенной ряд объемно и опирается на аппарат матричных исчислений (лишь матричное описание позволяет учесть любые возможные взаимосвязи между координатами модели). Поэтому опишем лишь алгоритм, согласно которому функционируют соответствующая процедура библиотеки анализа моделирующих программ:

gif-file, 2KB

  1. Пользователь указывает две координаты линейной непрерывной модели: вход и выход, – относительно которых требуется выполнить разложение и активирует соответствующую процедуру анализа.
  2. Библиотека анализа замещает все интеграторы модели специальной блок-схемой, изображенной на рисунке справа. К указанному пользователем входу библиотека подключает генератор импульса единичной амплитуды, и подготавливается к измерениям выходной координаты (начиная с этого шага алгоритма допустимо, чтобы программы не использовали графический интерфейс, т.е. скрывали технические подробности от пользователя).
  3. Далее, библиотека анализа инициирует процесс спец-симуляции модели. На каждом шаге выполняется измерение одного коэффициента разложения (в силу особенностей задающего сигнала значение выходной координаты равно коэффициенту разложения).

Достоинство описанного алгоритма заключается в том, что после упомянутой подстановки, величина шага симуляции моделирующей программы не имеет влияния на результат разложения. Однако процедура балансировки модели требует решения системы линейных алгебраических уравнений, порядок которой равен порядку исследуемой САР.

Если исходной информацией к разложению будут являться не ABCD-матрицы коэффициентов, а коэффициенты числителя и знаменателя ПФ разомкнутой САР, то балансировка модели в установившихся режимах движения может быть упрощена, поскольку итерационному решателю моделирующей программы достаточно прорешивать алгебраическое уравнение лишь первого порядка. Требуемые в этом случае трансформации исходной модели и результат разложения показаны на рисунке.

gif-file, 2KB

Частотный анализ моделей и систем

Частотный анализ модели или системы
Процесс идентификации частотной характеристики модели или системы сопровождаемый уточнением соответствующей частотной передаточной функции. Цель частотного анализа состоит в предоставлении исходных данных для решения задач оценки качества, коррекции, синтеза САР.

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

Вычислительные алгоритмы идентификации частотных характеристик моделей

Рабочие файлы: [fr_ABCD.mcd] [fr_ABCD_sz.mcd]

Главная особенность вычислительных алгоритмов идентификации ЧХ состоит в том, что они ни коим образом не анализируют входной и выходной сигналы, а поэтому могут применяться только к математическим моделям. Исходной информацией для расчета ЧХ являются коэффициенты модели, поэтому предварительная процедура её идентификации является обязательной.

Если программе известны коэффициенты числителя и знаменателя ПФ модели, то алгоритмизация процедуры расчета ЧХ не вызывает затруднений. Гораздо чаще программы математического моделирования располагают результатом идентификации в матричной ABCD-форме. Выведем формулу расчета ЧХ для данного случая:

gif-file, 2KB

ЧХ дискретной системы так же можно рассчитать с помощью представленной формулы, но нужно выполнить постановку: $jω←е^{jωT}$. В случае если требуется построить частотную характеристику для домена псевдочастоты λ (что может потребоваться при использовании устаревших методик), требуемая подстановка имеет вид: $jω←(2+jλT)/(2-jλT)$.

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

Во избежание необоснованных обвинений моделирующих программ, каждый специалист должен знать, что рассчитанная на основе известных коэффициентов ЧХ ни как не учитывает тех погрешностей, которые не может не вносить ЦВМ, согласно своей природе, в процесс симуляции моделей (см. рис.). Например, для перехода в частотную область используется подстановка идеального частотного оператора $s←jω$ (без вещественной составляющей), который не учитывает погрешностей дискретных квазианалогов интеграторов моделирующих программ (что, впрочем, методически верно).

gif-file, 2KB

На рисунке показаны переходные процессы системы, вызванные ненулевыми начальными условиями. Частотная характеристика разомкнутой системы очевидна (-40 дБ/дек. & -180°), и вычислительные алгоритмы частотного анализа ее подтверждают. Но симуляция движения координат модели выявляет несоответствия между временным и частотным доменами при переключении методов интегрирования (методу Эйлера с запаздыванием соответствует расходящийся переходный процесс; методу Эйлера с упреждением – сходящийся; методу трапеций – синусоида с постоянной амплитудой)

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

Процедуры идентификации ЧХ, основанные исключительно на спектральном или гармоническом анализе результатов измерений входного и выходного сигналов, имеют особую ценность в программах математического моделирования. Они обеспечивают наиболее эффективные и адекватные переходы от абстрактных математических моделей к проектируемым физическим прототипам и обратно, а так же являются одним из не многих инструментов оценки адекватности функционирования самих моделирующих программ в режиме симуляции.

Универсальных алгоритмов измерений ЧХ нет. Существуют три вида ограничений, которые порождают семейства методов измерений ЧХ, которые в той или иной степени могут быть оптимальными для исследования трех классов идентифицируемых линейных объектов (см. рис.). К упомянутым ограничениям относятся:

  1. Невозможность использования дельта-воздействий для объектов с эффектами насыщения.
  2. Затруднения в расширении диапазона частот свыше "трехдекадного барьера".
  3. Большие временные затраты на расширение динамического диапазона ниже уровня шума.

gif-file, 2KB

Если на вход системы подать весь спектр частот с единичными амплитудами, то определение частотной ПФ на основе измерительной информации упрощается:

$W(jω)=Y(jω)/X(jω)|_{X(jω)=1}=Y(jω)$.

Подобный единичный спектр имеет дельта-функция Дирака $δ(t)$. Реакция же систем на дельта-функцию называется функцией веса $w(t)$, поэтому частотную ПФ можно получить вычислением её Фурье-изображения:

$W(jω)=\FT{w(t)}$,

или же, для дискретных сигналов:

$W[jkω]=\DFT{w[n]}=_{n=0}^{N-1}∑w[n]е^{-jkωnT}$.

где: ω – частота первой гармоники в спектре сигнала длинной в $NT$ выборок; $k$ – порядковый номер гармоники (независимая переменная); $N$ – число выборок функции веса $w[n]$ (обычно кратно степени двойки); $k≤(N/2)$.

Подобный подход используется в простейших алгоритмах идентификации ЧХ и достаточно легко реализуется в любых математических программах:

gif-file, 2KB

Измерения частотных характеристик выполнены для двух дискретных фильтров с конечной импульсной характеристикой (FIR). Операция быстрого преобразования Фурье (нижний график) выполнена библиотечным блоком осциллограф программы VisSim (в свойствах блока осциллограф активирован режим вычисления БПФ для осциллограммы)

gif-file, 2KB

Измерения частотных характеристик выполнены для модели колебательного звена и инверсного фильтра Чебышева десятого порядка, т.е. для двух непрерывных систем с бесконечной импульсной характеристикой (IIR)

Сравнивая приведенные на рисунках ЧХ КИХ и БИХ-фильтров, легко понять суть "трехдекадного барьера" (присущего алгоритмам на базе БПФ), который проявляется при частотном анализе БИХ-фильтров, чьи ЧХ представляются в логарифмическом масштабе по оси частот. Очевидно, что, в этом случае, кратная двойке сетка анализируемых гармоник в выходном массиве процедуры БПФ ограничивает разрешение по частоте в НЧ диапазоне. Увеличение частотного диапазона на декаду требует увеличения времени симуляции в десять раз.

При измерении ЧХ реальных систем использование дельта-функции невозможно. Её замещают либо суперпозицией синусоид, либо белым шумом, либо другим сигналом конечной амплитуды. Основным методом расширения динамического диапазона ниже уровня шума является усреднение результатов повторных измерений. В случае использования алгоритмов измерения ЧХ на базе ДПФ для ослабления эффекта наложения частот используют методику взвешивания массива измерительной информации окнами Бартлета, Хэмминга, с хэннингом, Блэкмана, Хариса и т.д.

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

Идея алгоритмов идентификации ЧХ систем на основе технологий распознавания образов заключена в уточнении априорных предположений о порядке системы и порядке её астатизма. Согласно сценарию, алгоритм создает образ идентифицируемой системы (модель), и, в итерационном процессе, так подгоняет положение асимптот его частотной характеристики, дабы реакции идентифицируемой системы и образа на произвольные сигналы совпадали.

Простейшие алгоритмы данной группы обеспечивают погрешности идентификации не превышающие трети декады и 6 дБ, что вполне приемлемо при проектировании систем управления. В целях снижения погрешностей для идентификации выбирают фрагменты осциллограмм сигналов с наиболее широким спектром (переходные режимы функционирования).

Безусловно, группа измерительных алгоритмов в любых случаях более предпочтительна. Но в условиях, когда возможность вывода системы из производственного процесса для исследований отсутствует, или же невозможно подать на систему требуемые измерительные воздействия выбора нет – нужно использовать те сигналы, которые есть.