Частотная идентификация динамических систем

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

Результатом идентификации может быть импульсная или переходная характеристика объекта или спектральные характеристики. Эти характеристики могут использоваться в дальнейшем для структурной и параметрической идентификации.

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

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

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

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

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

Рассмотрим частотный метод идентификации, предложенный А.Г. Александровым [19]. Частотный метод относится к активным методам идентификации. Рассматривается передаточная функция объекта

                    .               (2.13)

В частотном методе идентификации для определения параметров передаточной функции (2.13) используется испытательный сигнал вида

                                     .                                (2.14)

На выходе объекта (2.13) измеряется сигнал

                          .                     (2.15)

Заменяя в (2.13) s = jw, можно получить комплексный коэффициент передачи

             ,        (2.16)

, (2.17)

Результат деления в операторах суммирования (2.17) округляется до ближайшего снизу целого числа.

При частотном методе идентификации эксперимент состоит в том, что тестовый сигнал вида (2.14) подается n раз с частотами w ₁, w ₂,…, w _n. Каждый раз по результатам эксперимента вычисляются числа l _r(w _r), v _r(w _r), r =1,2,…, n. Далее на основе (2.17) получается система 2 n линейных алгебраических уравнений для определения неизвестных параметров a _i, b _j:

            .       (2.18)

Пусть вначале на вход подан сигнал

                                     .                                (2.19)

После завершения переходных процессов выходной сигнал будет равен

                         .                    (2.20)

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

   . (2.21)

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

. (2.22)

Подавая далее на вход объекта сигнал  определяем l ₂ и v ₂ и т.д.

Рассмотрим пример частотной идентификации объекта, заданного следующей передаточной функцией:

                                  .                             (2.23)

Для проведения численного эксперимента соберем схему в Simulink, представленную на рис. 2.3.

Рисунок 2.3 – Схема частотной идентификации

В схеме используются блоки:

– генератор гармонического сигнала Sin Wave;

– передаточная функция Transfer Fnc;

– осциллограф Scope.

В блоке передаточная функция Transfer Fnc задаем параметры передаточной функции в соответствии с выражением (2.3), как показано на рис. 2.4. В первой строке Numerator coefficients задаются коэффициенты числителя начиная со старшего. Во второй строке Denominator coefficients задаются коэффициенты знаменателя передаточной функции, также начиная со старшего коэффициента.

В блоке генератора гармонического сигнала Sin Wave задаются (см. рис. 2.5): амплитуда Amplitude, равная 1; смещение синусоиды относительно нуля Bias, равное 0; частота Frequency, равная 3 рад/с; начальная фаза Phase, равная 0; интервал дискретизации Sample time, для непрерывного моделирования равный 0.

Рисунок 2.4 – Задание параметров блока Transfer Fnc

Рисунок 2.5 – Задание параметров блока Sin Wave

В блоке осциллографа Scope необходимо в меню Parameters, расположенном на его панели быстрого доступа, установить число осей, равное 2 и убрать ограничение на число отображаемых точек, убрав метку из окна Limit data points to last (по умолчанию равно 5 000), как показано на рис. 2.6.

Среди общих настроек моделирования требуется определить общее время моделирования и максимальный шаг интегрирования.

Рисунок 2.6 – Настройка осциллографа Scope

Время моделирования задается в окне модели Simulink, расположенном на панели быстрого доступа. По умолчанию это время равно 10. Требуемое врем моделирования вычисляется частоте задающего сигнала на входе. В рассматриваемом примере для ПФ (2.4) находим нули и полюса (корни полиномов числителя и знаменателя):

                          .                      (2.24)

Для нахождения нулей нужно приравнять к нулю числитель ПФ и решить полученной линейное алгебраическое уравнение. Для нахождения полюсов нужно приравнять к нулю знаменатель ПФ и решить полученной линейное алгебраическое уравнение. Для систем высокого порядка можно использовать численные процедуры Matlab. Например, оператор ltiview позволяет построить распределение нулей и полюсов линейной системы автоматически.

Кроме того учтем, что частота гармонического сигнала равна 3 рад/с, что соответствует корню

                                        .                                   (2.25)

Отметим, что задаваемая частота не должна быть слишком высокой, т.е. корень (2.25) не должен слишком сильно превышать значения корней (2.24). Кроме того, нежелательно, чтобы он совпадал с корнями (2.24).

По выражению (2.25) вычисляем время моделирования:

                     с.               (2.26)

Если после запуска моделирования графики на осциллографе получаются негладкими, то необходимо задать максимальный шаг моделирования в меню Simulation->Configuratio Parameters, как представлено на рис. 2.7, где задан максимальный шаг интегрирования 0.01 с.

Рисунок 2.7 – Задание максимального шага интегрирования

Далее схема, представленная на рис. 2.3, модифицируется с целью автоматизации вычисления экспериментальных коэффициентов. Модификация схемы моделирования представлена на рис. 2.8.

Дополнительная часть схемы вычисляет коэффициенты преобразования Фурье.

Для первого эксперимента настраиваем параметры следующим образом:

1) в блоке генератора Sin Wave 2 устанавливаем амплитуду, равную 1, частоту равную 3 и фазу, равную 0;

2) в генераторе Constant задаем 0;

Рисунок 2.8 – Схема эксперимента частотной идентификации

3) в блоке Switch задаем параметр Threshold (время переключения), равный:

            с.      (2.27)

Такие параметры означают, что на интегратор 9 периодов синусоиды будет поступать 0, а затем один период будет поступать выходной сигнал.

4) В блоке Gain задаем коэффициент усиления, равный:

                       с.                 (2.28)

В этом случае при моделировании вычисляется выражением вида:

                           с.                     (2.29)

5) Запустив на моделирование схему получим на осциллографе Scope 2 в конечный момент моделирования значение коэффициента .

Далее задаем в генераторе Sin Wave 2 начальную фазу, равную -1.57 рад моделируем снова схему и получаем значение коэффициента . Так как осуществлен сдвиг фазы в генераторе Sin Wave 2 на -1.57 рад, то данное значение вычислено по формуле

                          с.                    (2.30)

Далее меняет частоту в генераторах на 1.5 рад/с, пересчитываем коэффициент усиления (2.28) при новой частоте и повторяем моделирование. При этом время моделирования сохраняем прежнее, а время переключения переключателя вычисляется заново:

 с. (2.31)

Тогда значения коэффициентов для новой частоты равны: ; .

Отметим, что можно измерять амплитуду и фазу выходного сигнала графически и определять коэффициенты  в соответствии с выражениями

                         .                    (2.32)

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

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

Представим искомую ПФ в виде

                                 .                            (2.33)

Подставим в ПФ (2.11х) замену , в результате чего получим комплексный коэффициент передачи:

                .           (2.34)

Теперь воспользуемся уравнениями, связывающими искомые коэффициенты ПФ (2.34) с экспериментальными данными (2.28) и (2.30):

            .       (2.35)

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

                       .                  (2.36)

Решим полученную систему с помощью Matlab. Для этого преобразуем ее к матричной форме:

                 .            (2.37)

Программа решения системы (2.37) представлена ниже.

%---------------------------------------------------

w1=3; w2=1.5;

l1=-0.0077; v1=0.4382;

l2=0.2832; v2=0.8967;

A=[l1 w1*v1 -1 0; v1 -w1*l1 0 w1; l2 w2*v2 -1 0; v2 -w2*l2 0 w2];

B=[w1^2*l1; w1^2*v1; w2^2*l2; w2^2*v2];

resh1=A^(-1)*B

%----------------------------------------------------

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

                                 .                            (2.38)

Сравнивая (2.36) и (2.23), видим, что погрешность определения коэффициентов передаточной функции в данном примере не превысила 5 %.

 

2. 5. Проектное задание 6

На основе примера, приведенного выше провести определение параметров передаточной функции вида (2.33) для параметров, заданных вариантом из табл. 2.4.

Таблица 2.4 – Варианты заданий для частотной идентификации объекта

№
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

Продолжение таблицы 2.4

17
18
19
20
21
22
23
24
25
26
27
28
29
30

По результатам идентификации вычислить погрешность по каждому коэффициенту.