17.01.2018 08:35
938MeV 8 1900 Источник
В практиктических приложениях ТАУ часто необходимо точно и качественно идентифицировать объект управления. В этой статье речь пойдет об идентификации объекта управления частотным методом. Данный метод применим, когда есть возможность физически протестировать объект управления синусоидальным входным воздействиямем, изменяя частоту в широком диапазоне. Если это условие соблюдено, то результат, как правило, оправдывает самые оптимистичные ожидания.
Полюса передаточной функции

Что необходимо знать


Начнем с краткой теоретической справки. Для понимания материала статьи читателю необходимо иметь представление о следущих вещах:

  • Преобразование лапласа
  • Линейные динамические системы
  • Характеристическое уравнение
  • Передаточная функция
  • Дискретное преобразование Фурье
  • Характеристика Боде: ЛАЧХ и ЛФЧХ

Отклик на синусоидальное воздействие


Как известно, отклик динамической системы на синусоидальное воздействие — есть синусоида с той же частотой, но отличной амплитудой и фазой. Именно эти две характеристики: амплитуда и фаза образуют график Боде, то есть ЛАЧХ и ЛФЧХ. По сути, задача идентификации динамической системы сводится к экспериментальному нахождению этих двух графиков.

Для примера рассмотрим уравнение гармонического осциллятора с ненулевой правой частью как образец линейной динамической системы второго порядка

$\ddot{x}(t)+ 2\gamma\dot{x}(t)+\omega_n^2 x(t) = u(t)$


Обозначим преобразование Лапласа произвольной функции $f(t)$ через $\hat{f}(s) = \mathcal{L}\lbrace f(t) \rbrace$. Применим преобразование Лапласа к обеим частям уравнения

$s^2\,\hat{x}(s)+ 2\gamma\, s\, \hat{x}(s)+\omega_n^2 \, \hat{x}(s) = \hat{u}(s)$


Тогда характеристическое уравнение динамической системы будет

$\lambda^2 + 2\gamma\, \lambda+\omega_n^2 = 0$


А передаточная функция

$H(s) = \frac{\hat{x}(s)}{\hat{u}(s)} = \frac{1}{s^2 + 2\gamma\, s+\omega_n^2}$


Искомые характеристики ЛАЧХ и ЛФЧХ получаются заменой $s \rightarrow j\omega$ и взятием модуля и аргумента функции $H(j\omega)$ соответственно

$\begin{split} a(\omega) &= \left|H(j\omega)\right|\\ \varphi(\omega) &= \arg{H(j\omega)} \end{split}$


Здесь $a(\omega)$ — амплитуда (Magnitude), а $\varphi(\omega)$ — фаза (Phase) соответствующей компоненты на частоте $\omega$. В результате получится нечто похожее на это
Характеристика Боде: ЛАЧХ и ЛФЧХ
Характеристика Боде: ЛАЧХ и ЛФЧХ

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

Подавая синусоидальное воздействие на вход, мы будем наблюдать примерно такую картину
Синусоидальное входное воздействие и отклик системы
Синусоидальное входное воздействие и отклик системы

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

Дискретное преобразование Фурье


Допустим, мы измерили два графика, входное воздействие $u(t)$ и отклик системы $x(t)$, для какой-то заданной частоты входного воздействия $f$. Необходимо, чтобы оба измерения были сделаны синхронно и с одинаковым периодом дискретизации $T$, который должен быть известен достаточно точно. Таким образом, мы имеем два набора дискретных значений $u_n = u(t_n)$ и $x_n = x(t_n)$, где $n = 0,1,\dots,N$, а $u_n$ и $x_n$ — значения $u(t)$ и $x(t)$ в соответствующие дискретные моменты времени $t_n$. Следует заметить, что полное время измерения должно быть достаточно большим, чтобы захватить хотя бы несколько (я бы рекомендовал 3 или больше) периодов колебаний.

К дискретным наборам $\lbrace u_n\rbrace$ и $\lbrace x_n\rbrace$ можем применить дискретное преобразование Фурье. Дискретное преобразование Фурье переводит вышеуказанные сигналы из временной области в частотную, то есть

$\begin{split} \lbrace u_n \rbrace &\rightarrow \lbrace U_k \rbrace \\ \lbrace x_n \rbrace &\rightarrow \lbrace X_k \rbrace \end{split}$


где $k = 0,1,\dots,N/2$, а $U_k$ и $X_k$ — соответствующие комплексные амплитуды $k$-й гармоники.

Применение метода


Теперь применим дискретное преобразование Фурье к нашим двум сигналам. На рисунке ниже показаны графики амплитуд $|U_k|$ и $|X_k|$
Амплитудный спектр входного и выходного сигнала
Амплитудный спектр входного и выходного сигнала

График также получен при обработке экспериментальных данных. Как можно видеть, на графике выделяется острый пик на частоте $f = 0.4$ Hz. Это «несущая» частота входного воздействия, то есть частота, на которой происходило возбуждение объекта управления. На обоих графиках, входном и выходном, для данной гармоники наблюдается пик. Из двух значений комплексных амплитуд $U_k$ и $X_k$ на данной частоте получим значение передаточной функции $H_k$.

Как мы помним,

$H(s) = \frac{\hat{x}(s)}{\hat{u}(s)}$


В случае же дискретных $u_n$ и $x_n$ имеем

$H_k = \frac{X_k}{U_k}$


Тогда, на графики ЛАЧХ и ЛФЧХ можем нанести точку:

$\begin{split} a_k &= \left|H_k\right|\\ \varphi_k &= \arg{H_k} \end{split}$


Какую же частоту имеет гармоника с номером $k$? Отвечаем: частота гармоники дается формулой

$f_k = \frac{k}{NT}$


Можно также записать

$\omega_k = \frac{2 \pi k}{NT}$


Выражение же для, собственно, гармоники следующее

$h_k(t) = a_k \sin\left(\omega_k t + \varphi_k\right)$


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

Проделав такой эксперимент и обработав данные, получаем таблицу частотных характеристик

Таблица частотных характеристик

Можно сразу построить измеренные графики ЛАЧХ и ЛФЧХ. Выглядит это примерно вот так

ЛАЧХ и ЛФЧХ, построенные по экспериментальным данным
ЛАЧХ и ЛФЧХ, построенные по экспериментальным данным

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

Теперь нужно воспользоваться MATLAB, а конкретнее System Identification Toolbox. В составе этого Toolbox'а есть System Identification App — интерактивное приложение, которое может идентифицировать вашу систему по частотным данным. К слову, есть и другие опции, как то, идентификация напрямую по временным измерениям.

Для корректной идентификации необходимо знать порядок системы. Здесь нам поможет график ЛФЧХ. Чтобы узнать порядок системы, посмотрите на график ЛФЧХ и прикиньте, на сколько раз по 90 градусов отстает фаза на высокой частоте. Количество раз по 90 градусов и будет порядок (знаменателя) вашей системы.

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

tf1 =
 
  From input "u1" to output "y1":
   -97.64 s + 1.063e04
  ---------------------
  s^2 + 1.547 s + 176.7
 
Name: tf1
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 response data "h".
Fit to estimation data: 92.26% (stability enforced)  
FPE: 120.6, MSE: 112.3

Согласно отчету об идентификации, входные данные удовлетворяются найденной моделью на 92.26%. Имеем следующую передаточную функцию физической системы

$H(s) = \frac{-97.64\,s+1.063\times10^4}{s^2+1.547\,s+176.7}$


Теперь у вас есть объект типа IDTF c именем tf1, с которым можно делать все то же самое, что и с любой другой LTI System в MATLAB. Кроме того объект содержит в себе информацию о неопределенности внутренних параметров, и если мы построим Bode Plot, то на графике можно вызвать отображение доверительных интервалов. В настройках можно указать количество стандартных отклонений.

Идентифицированная модель системы с доверительными интервалами
Идентифицированная модель системы с доверительными интервалами

Чтобы проверить правильность идентификации можно совместить экспериментальные точки с графиками частотных характеристик идентифицированной модели

Совмещенные графики для ЛАЧХ и ЛФЧХ
Совмещенные графики для ЛАЧХ и ЛФЧХ

Заключение


Использование данного метода существенно облегчает процесс проектирования системы управления. Автор статьи успешно разработал и воплотил в железе LQG котроллер для гидравлической цепи специального назначения. Кроме того, данный метод можно использовать для оценки эффективности системы управления, при условии, что вам физически доступно входное воздействие на объект управления со стороны нежелательных возмущений.

Комментарии (8)


  1. Arastas
    17.01.2018 18:41
    #10621612

    Ого, месяц ТАУ продолжается!


    Теперь серьёзнее.
    У вас на рисунке с данными экспериментов фаза начинается от -90, а у ЛАХ наклон на низких частотах. А в той части, где вы совмещаете эксперимент с моделью, там фаза от нуля идёт, и у ЛАХ нет наклона. Куда интегратор дели?


    Чтобы узнать порядок системы, посмотрите на график ЛФЧХ и прикиньте, на сколько раз по 90 градусов отстает фаза на высокой частоте. Количество раз по 90 градусов и будет порядок (знаменателя) вашей системы.

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


    Кстати, не обязательно строить по точкам. Как вариант, можно подать достаточно частотно богатый входной сигнал, и построить частотную характеристику сразу в диапазоне (tfestimate в Matlab).


  1. 938MeV Автор
    17.01.2018 21:27
    #10621896

    Каюсь. Это два разных эксперимента. График привел просто показать как пример. Там наклон был оттого что это система второго порядка помноженная на 1/s. Грубо говоря, мотор там стоял еще.


    Пожалуй, да. Про нули передаточной функции я как-то не подумал. А как тогда оценить порядок системы?


    Про "богатый" сигнал. А вы представляете себе как это будет выглядеть? Когда объект управления механический — это просто будет жуть и страх! Все будет ходить ходуном, хаотично трястись и с грохотом разваливаться на части.
    А если вы имеете ввиду ступеньку на вход подать. Так тоже можно. Но, на мой взгляд, менее надежно. Это хорошо работает для систем второго порядка. А что если порядок системы будет выше, да еще и нули в произвольных местах. Тут уже точность определения по ступеньке будет гораздо хуже.


    1. Arastas
      17.01.2018 23:11
      #10621988

      А как тогда оценить порядок системы?

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


    1. Arastas
      17.01.2018 23:17
      #10622000

      Про "богатый" сигнал. А вы представляете себе как это будет выглядеть? Когда объект управления механический — это просто будет жуть и страх! Все будет ходить ходуном, хаотично трястись и с грохотом разваливаться на части.

      Прекрасно представляю, ничего там особо страшного не происходит. :) И на гидропривода подавал шум, и на механику, как весом с десяток киллограм, так и в несколько тонн. Просто надо иметь базовое представление о физике процесса и разумно выбирать форму сигнала и частотный состав.
      На самом деле, если вам шум не нравится, то есть более приятные варианты — синусоида с нарастающей частотой. Кроме того, что гладкая и хорошая, она ещё позволяет в одном эксперименте оценить целый диапазон частот. Достаточно распространённый сигнал на практике. Лучше, чем ступенька. :)


  1. Arastas
    17.01.2018 23:16
    #10621998

    //промахнулся


    1. 938MeV Автор
      18.01.2018 11:21
      #10622544

      А как анализировать сигнал с нарастающей частотой? Подскажите?


      1. Arastas
        18.01.2018 13:05
        #10622826

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


  1. SNEURO
    18.01.2018 11:19
    #10622534

    Спасибо вам за статью.
    Встречал людей, которые частными методами идентифицировали крупные, инерционные тепловые объекты. Частоты там приходилось выбирать совсем низкие. Процесс идентифицировали занимал не один час. Интересно, стоит ли оно затраченного времени.
    Подскажите, есть ли преимущество у частотного метода перед step response, например?