Анализ устойчивости и коррекция сар по частотным характеристикам и по полюсам

Лабораторная работа №4 по курсу «Управление в технических системах»

Введение

В лабораторных работах, выполненных Вами в прошлом семестре, были рассмотрены основные процедуры работы в SimInTech применительно к моделированию и анализу динамических процессов в линейных системах автоматического управления (САУ). Выполнив в прошлом семестре самостоятельно также и домашнее задание, Вы “закрепили” полученные знания.

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

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

Кроме того, существует значительное количество методов моделирования и анализа динамических систем в SimInTech, пока не известных Вам.

Поэтому лабораторный практикум настоящего семестра направлен, во-первых, на изучение методов моделирования и анализа нелинейных динамических систем и, во-вторых, на освоение Вами новых процедур работы в SimInTech.

Одна из задач настоящей лабораторной работы посвящена анализу динамических систем с запаздыванием, которые в Теории Управления обычно относят к классу особых динамических систем.

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

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

Цель работы

Анализ динамических систем с запаздыванием

Блок Идеальное запаздывающее звено

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

где T(x,t) – какая-то скалярная субстанция (например, температура потока), переносимая с постоянной скоростью u; х – продольная координата.

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

где: T(L,s) – изображение по Лапласу сигнала на выходе из трубопровода; T(0,s) – изображение по Лапласу сигнала на входе в трубопровода; τ = L/u – постоянная запаздывания (время транспортировки).

Часто передаточную функцию идеального запаздывающего звена аппроксимируют типовыми линейными звеньями, например, цепью из n последовательно соединенных апериодических звеньев 1-го порядка:

В учебной литературе нередко утверждается, что если n ~ 6…8, то этого достаточно для аппроксимации передаточной функции идеального запаздывающего звена. Покажем, что это не совсем так.

Используя полученный в прошлом семестре опыт работы в SimInTech, сформируйте “с чистого схемного окна” структурную схему, подобную рисунку (Рисунок 1).

На 1-ом этапе перенесите из “Линейки типовых блоков" в Схемное окно необходимые блоки, расположите их на требуемые места и соедините линиями связи.

Второй этап требует пояснений. Главная особенность структурной схемы на рисунке (Рисунок 1) – использование векторизованной обработки и передачи данных.

Переместите курсор на закладку Скрипт и выполните щелчок левой клавишей “мыши”: откроется окно Редактора Скрипта Проекта. Введите с клавиатуры текст, идентичный приведенному на рисунке (Рисунок 2).
const n1=8; const n2=20;
Константы n1 и n2 будут задавать количество последовательно соединенных апериодических звеньев 1-го порядка в двухпараллельных цепях, аппроксимирующих свойства идеального запаздывающего звена. Закройте Скрипт (нажатием кнопки Применить).

Рисунок 1. Структурная схема

Рисунок 2. Константы проекта

Откройте окно свойств блока Ступенька и введите в диалоговой строке параметры смещенного ступенчатого воздействия:Время срабатывания – 2, Начальное значение – 0, Конечное значение – 1. Введенное означает, что через 2 секунды после начала моделирования сигнал на выходе блока скачком изменится с 0 (нуля) до 1 (единицы).

Откройте окно свойств блока Идеальное запаздывание и введите в 1-ой строке (Время запаздывания) число 2 (два), что означает что данный блок реализует постоянное запаздывание 2 секунды.

Число, введенное во второй диалоговой строке задает начальный размер стека данных, в который будут записываться данные на входе блока после каждого шага интегрирования. Если стек заполнится полностью, то он будет увеличен до 1200, если снова заполнится – до 1400 и т.д. Выходной сигнал определяется линейной интерполяцией значений в стеке данных. Оставьте начальный размер стека (по умолчанию).

Откройте окно свойств верхнего блока Апериодическое звено 1-го порядка (8 последовательных звеньев) и заполните его так же, как это выполнено на рисунке (Рисунок 3).

Рисунок 3. Свойства «верхнего» Апериодического звена 1-го порядка

В 1-ой диалоговой строке (Коэффициенты усиления) введено n1#1. Это означает, что в данной строке введен числовой вектор из n1 (8) единиц (1). Можно было ввести данную строку и так (вектор-строка): [1, 1, 1, 1, 1, 1, 1, 1] (через запятую в квадратных скобках). Символ # в диалоговых строках эквивалентен предлогу “по” ⇒ n1 элементов по1.

В последней диалоговой строке (Начальные условия) аналогичным образом задан вектор из n1 (восьми) нулей.

В средней (во 2-ой) диалоговой строке задан вектор из n1 (восьми) одинаковых постоянных времени, равных 2/n1 = 2 c/8 = 0.25 c.

По аналогии с предыдущим заполните диалоговое окно для другого блока Апериодическое звено 1-го порядка (см. Рисунок 4). Очевидно, что данный блок предназначен для аппроксимации идеального запаздывающего звена цепью из 20-ти последовательно соединенных апериодических звеньев 1-го порядка.

Рисунок 4. «Нижнее» Апериодическое звено 1-го порядка

Рисунок 5. Свойства Демультиплексора

Откройте диалоговое окно блока Мультиплексор в цепи, реализующей 8 последовательно соединенных звеньев, и заполните диалоговую строку (Количество портов) во вкладке Свойства, указав там 2 порта.

Откройте диалоговое окно блока Демультиплексор и заполните его, как это выполнено на рисунке (Рисунок 5).

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

Поскольку алгоритм работы верхнего блока Апериодическое звено 1-го порядка (см. Рисунок 3) – векторизован, то на вход блока должен поступать векторный сигнал, размерностью n1 (8 элементов). Поэтому к скалярному сигналу от блока Ступенька необходимо добавить n1-1 (7) сигналов, чтобы после блока Мультиплексор векторный сигнал имел размерность n1.

Векторный сигнал, поступающий на 2-ой (нижний) порт блока Мультиплексор, сформирован из (n1-1) сигналов на 1-ом выходном порте блока Демультиплексор (см. Рисунок 1).

Фактически реализован сдвиг “жил” сигналов. Рассмотрим реализацию сдвига, “отталкиваясь” от сигнала блока Ступенька.

Сигнал от блока ступенька поступает на 1-ю “жилу” входного порта ⇒ далее “проход” через Апериодическое звено ⇒ далее сигнал 1-ой выходной “жилы” Демультиплексора подается на 2-ую входную “жилу” Мультиплексора ⇒ далее “проход” через Апериодическое звено ⇒ далее сигнал 2-ой выходной “жилы” Демультиплексора подается на 3-ю входную “жилу” Мультиплексора и т.д.

В итоге на втором выходном порте блока Демультиплексор будет сигнал, который n1 раз “прошел” через Апериодическое звено 1-го порядка.

По аналогии с рисунком (Рисунок 5) заполните диалоговые окна блоков Демультиплексор в цепи, аппроксимирующей звено идеального запаздывания 20-ю последовательно соединенными звеньями.

На этом формирование структурной схемы и ее параметров завершено.

Рисунок 6. Параметры расчета

Переместите курсор мыши на любое свободное место схемы проекта и нажмите правую клавишу мыши. В контекстном меню выберите пункт Параметры расчета и заполните диалоговое окно так же, как это выполнено на рисунке (Рисунок 6). Заполнив окно Параметры расчета, закройте его щелчком “мыши” по кнопке Ок. Запустите задачу на счет. Мгновенно на графике отобразятся результаты расчета. Используя процедуры редактирования окна графика, придайте ему вид, близкий рисунку (Рисунок 7), где линии: штриховая – блок идеальное запаздывание, пунктирная – цепь из 8 блоков, сплошная – из 20 блоков.

Рисунок 7. Сравнение переходных процессов

Сравнение графиков переходных процессов показывает, что даже при аппроксимации блока Идеальное запаздывание цепью из 20-ти последовательно соединенных звеньев “фронт” скачка существенно “размыт”, а при аппроксимации цепью из 8-ми блоков – тем более.

Резюме: сравнение данных результатов расчета переходных процессов показало, что вышеупомянутое утверждение о достаточности для аппроксимации цепи из 6…8 последовательно соединенных Апериодических звеньев 1-го порядка является фактически некорректным для входных воздействий типа “ступенька”. Дополним сравнение динамических свойств “классического” Идеального запаздывающего звена и его “аппроксиматоров” сопоставлением амплитудно-фазовых частотных характеристик.

Построение амплитудно-фазовых частотных характеристик для сопоставляемых звеньев

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

На первом этапе откройте новое схемное окно Схема автоматики,сформируйте блоки схемного окна как показано на рисунке 1.8 и соедините их линиями связи.

Рисунок 8.

Переместите курсор на блок Полином n-ой степени и задайте свойство Коэффициенты полинома равными [-2, 1], как показано на рисунке 1.9 – численно это будет равно диапазону изменения логарифма частоты lg(ω), поскольку мы установим параметры расчета соответствующим образом (конечное время счета = 3 с, шаг = 0.001 с). Таким образом, за 3 секунды счета будет сделано 3000 шагов и выход блока будет равномерно увеличиваться от -2 до +1.

Рисунок 9.

Входным сигналом для построения ЛАХ, ФЧХ и критерия Найквиста является десятичный логарифм частоты lg(ω), заданный блоком Полином n-ой степени.

Откройте блок Язык программирования. Входной сигнал задается ключевым словом input (строка input lgw;) что означает определение переменной lgw как входного сигнала в блок. Далее, исходя из формулы вычисления десятичного логарифма, находим значение собственно частоты колебаний звеньев ω=10lgw.

Следует задать в блоке Язык программирования значения полиномов передаточных функций по степеням «s» N(s) и L(s), K – общий коэффициент усиления звена (системы) = 1 – для каждой из функций сопоставляемых звеньев.

Определяем передаточную функцию следующим соотношением:

Логарифмическая амплитудно-частотная характеристика (или логарифмическая амплитудная характеристика, далее ЛАХ) определяется формулой:

где A(ω)=|W()|.

Фазочастотная характеристика (ФЧХ) определяется следующим соотношением:

Годограф Найквиста определяется функцией:

Для задания данных критериев анализе следует переместить курсор на блок Язык программирования, открыть его и набрать текст, представленный на рисунке (Рисунок 10).

Рисунок 10.

Проанализируйте представленный текст самостоятельно – вы должны (!?) понять, что означает каждая его строка и зачем она нужна.
Следует помнить, что в данном примере мы пользуемся векторными сигналами для расчета частотных характеристик. Поэтому программирование вывода результата расчета частотных критериев выполнено в векторной форме в блоке Язык программирования:
output Lm[3], fi[3], Re[3], Im[3];

Для корректности вывода результатов расчета в векторной форме следует использовать блок Размножитель. Зададим коэффициент размножения равным 3#1, как показано на рисунке (Рисунок 11), что означает 3 одинаковых векторных сигнала, которые подаются на входы в блок Язык программирования и блоки Фазовый портрет для ЛАХ и ФЧХ.

Рисунок 11. Свойство Размножителя

Рисунок 12. ФЧХ

Рисунок 13. Годографы Найквиста

Рисунок 14. ЛАХ

На рисунках (Рисунок 12-Рисунок 14) (в качестве “эталона” для Ваших графиков) приведены результаты расчета годографов АФЧХ (годографов Найквиста), фазовых частотных характеристик (ФЧХ) и логарифмических амплитудных характеристик (ЛАХ), соответственно. Чёрными линиями представлены характеристики Идеального запаздывающего звена, синими (?) линиями – для цепи из 8 звеньев, и красными (?) – для цепи из 20 блоков. Убедитесь в этом самостоятельно…

Анализ графиков частотных характеристик показывает, что в области низких частот (менее 1 Гц) аппроксимирующие цепи близки к Идеальному запаздывающему звену.

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

Резюме: вышеприведенное утверждение о достаточности для аппроксимации цепи из 6…8 последовательно соединенных Апериодических звеньев 1-го порядка является относительно корректным только для медленно изменяющихся входных воздействий.

Определение устойчивости линейных систем с запаздыванием

В прошлом семестре изучение основных процедур работы в SimInTech Вы проводили в рамках демонстрационно-ознакомительной задачи, в которой структурная схема САР имела вид, близкий рисунку (Рисунок 15).

Рисунок 15. Структурная схема САР

Объект управления с передаточной функцией W₂(s), соответствовал типовому звену (колебательному) с коэффициентом усиления k₂ = 1.0; постоянной времени T₂ = 1 c; коэффициентом демпфирования b = 0.5; начальные условия – нулевые.

Местная обратная связь с передаточной функцией W₃(s), соответствовала типовому звену – апериодическому 1-го порядка со свойствами: k₃ = 0.6; T₃ = 5 c.

Локальное сравнивающее устройство обеспечивало отрицательную обратную связь, т.е. “работало” в режиме обычного вычитания.

В ходе выполнения ознакомительной задачи Вы подобрали коэффициент усиления k₁ интегрирующего регулятора с передаточной функцией W₁(s) таким образом, что при подаче ступенчатого управляющего воздействия u(t) = 0.8·1(t) перерегулирование отсутствовало (т.е. ymax <= 0.8) и время переходного процесса не превышало 20 с. Значение коэффициента усиления k₁ интегрирующего регулятора оказалось равным 0.35.

Рисунок 16.

В настоящей лабораторной работе Вам предстоит скорректировать структурную схему САР, добавив в “прямую” цепь Идеальное запаздывающее звено. Структурная схема скорректированной САР должна иметь вид, близкий рисунку (Рисунок 16).

Этапы, которые Вы должны выполнить:
  1. Вы должны фактически снова сформировать математическую модель динамики “знакомой” САР (с найденным ранее “оптимальным” значением k₁ = 0.35).
  2. Определить критическое значение постоянной запаздывания τкрит в Идеальном запаздывающем звене.
  3. Варьируя постоянную запаздывания в Идеальном запаздывающем звене в пределах 0.1·τкрит …0.9·τкрит (4 значения) выполнить моделирование переходных процессов.
  4. Выполнить анализ полученных результатов.

Блок Переменное транспортное запаздывание

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

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

и основана на допущении о постоянстве линейной скорости переноса распадающейся субстанции в пределах участка для каждого момента времени при граничных условиях y(0,t) = f(t) и начальных условиях y(z,0) = y0(z), z ∊ [0,L]. В уравнении (*) y(t) – переносимая скалярная субстанция, u(t) – скорость переноса, L – длина участка переноса скалярной субстанции, z – пространственная (продольная) координата.

После ввода безразмерной пространственной координаты x = z/L и мгновенного времени переноса скалярной субстанции в пределах участка τ(t) = L/u(t) уравнение записывается как

а начальные условия принимают вид y(z,0)=y0(x), x∊ [0,1].

Вводя дополнительное дифференциальное уравнение для новой переменной θ

дифференциальное уравнение (*) принимает вид:
Используя преобразование Лапласа, получаем решение в виде уравнения:
где сомножитель y(0, θ-1) описывает составляющую, обусловленную только транспортным запаздыванием, а сомножитель exp(-λ·τзап(t)) описывает ослабление выходного сигнала блока, обусловленное только распадом субстанции за время ее пребывания в пределах участка транспортного запаздывания.

При вычислении yвых(t) в расчете используется запоминание текущих значений t, yвхθ(t) в стековой таблице (Таблица 1) и последующая обработка табличных данных.

Таблица 1.
Индекс записи Модельное время t θ(t) yвх(t) yвых(t) τзап(t)
0 0 0 y0 y0 τ0зап
1 t1 θ1 y1 y0 τ1зап
2 t2 θ2 y2 y0 τ2зап
y0
J tj θj yj y0 τjзап
y0
k-1 tk-1 θk-1<1.0 yk-1 y0 τk-1зап
k tk θk=1.0 yk y0 τkзап = tk
k+1 tk+1 θk+1<1.0 yk+1 y0 τk+1зап
m tm θm ym yвых(tm) τmзап
При ttk значение yвых(t) = y0, а при t > tk значение определяется с использованием данных таблицы (Таблица 1) по алгоритму tm → θmm-1) → y(θm-1) ≡ yвых(t). Последняя процедура (вычисление yвых(t)) проводится с использованием линейной интерполяции данных таблицы (Таблица 1).

Расчет фактического времени запаздывания τзап = t в блоке Переменное транспортное запаздывание при u(t) ≠ 0 проводится следующим образом:

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

На 1-ый входной порт подается сигнал, соответствующий значению скалярной субстанции на входе в участок транспортировки. На 2-ой входной порт подается сигнал, соответствующий значению мгновенного времени переноса скалярной субстанции в пределах участка транспортировки.

На 1-ом выходном порте формируется сигнал, соответствующий значению скалярной субстанции на выходе из участка транспортировки. На 2-ом выходном порте формируется сигнал, соответствующий значению времени пребывания “метки” скалярной субстанции в пределах участка транспортировки.

Блок Переменное транспортное запаздывание имеет 2 свойства (две диалоговые строки). Для работы блока необходимо задать:

По умолчанию блок Переменное транспортное запаздывание реализует алгоритм преобразования скалярного входного сигнала для нераспадающейся скалярной субстанции (λ=0).

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

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

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

Переместите курсор “мыши” на Панель инструментов и выберите в разделе меню Сервис пункт Сигналы…. Для задания глобального сигнала величины запаздывания следует нажать кнопку Добавить сигнал и ввести в колонках Имя сигнала(U1), Название (Время запаздывания), Режим (Вход), Значение (2), Способ расчета (Переменная), как показано на рисунке Рисунок 17.

Рисунок 17.

Закройте это диалоговое окно, выполнив щелчок “мышью” по кнопке Ок. Переместите курсор на блок Запись в список сигналов и откройте его диалоговое окно 2-х кратным щелчком “мыши”. Введите имя сигнала U1. Соедините блоки линиями связи, как это выполнено на рисунке (Рисунок 19).

Рисунок 18.

Задайте в свойствах блока Идеально запаздывающее звено значение U1 в качестве постоянной запаздывания. Отметим, что на самом деле свойство tau блока в процессе моделирования будет переменным, так как его значение численно равно фактическому времени запаздывания τзап(t) в блоке Переменное транспортное запаздывание.

Откройте свойства блока Синусоида и введите значение амплитуды (1), частоты (0.5) и сдвига фазы (0).

Откройте диалоговое окно свойств блока Произвольное кусочно-линейное воздействие и введите в первой строке [0, 5, 10, 20, 25, 40], а во второй диалоговой строке [2, 2, 5, 5, 2, 2]. Закройте это диалоговое окно.

Рисунок 19.

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

Откройте диалоговое окно Параметры расчета и введите: Время интегрирования40 с; Минимальный шаг интегрирования0.01 с; Максимальный шаг интегрирования0.01 с. Остальные параметры – по умолчанию.

Не забудьте сохранить проект на диск под оригинальным именем.

Выполните расчет переходного процесса (щелчок по кнопке Пуск). Если Вы выполните оформление графика, то его вид будет подобен рисунке (Рисунок 20). Данные расчета показывают, что блок Идеальное запаздывающее звено фактически реализовал математическую модель блока Переменное транспортное запаздывание.

Рисунок 20.

Убедитесь самостоятельно в том, что если мгновенное время запаздывания в блоке Переменное транспортное запаздывание постоянно, то блок фактически эквивалентен Идеальному запаздывающему звену.

Исследование известных динамических задач методами структурного моделирования

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

  1. Уравнением Ван-дер-Поля:
    в диапазоне от t = 0 до t = 40 сек, если y(0) = 1, y’(0) = 0.

    Используя типовые блоки библиотеки Данные (Временной график и Фазовый портрет) построить зависимости y(t) и траектории на фазовой плоскости (y, y’), если:

    • b = 1 и варьируемые значения а а = -1; 0; 1; 5;
    • a = 5 и варьируемые значения bb = 1; 2; 5; 10.

    Завершив моделирование, по виду переходных процессов сделайте вывод о роли значений a и b на характер движения системы.

  2. Уравнением Матье:
    в диапазоне от t = 0 до t = 200 сек, если y’(0) = 0, а y(0) = var.

    Используя типовые блоки библиотеки Данные (Временной график и Фазовый портрет) построить зависимости y(t) и траектории на фазовой плоскости (y, y’), если:

    • p = 1; e = 0.1; m = 0.2; b = 1; w = 1, а y(0) = 0.01; 0.1; 1.0.
    • y(0) = 0.5; e = 0.1; m = 0.2; b = 1; w = 1, а p = 0.5; 0.9; 0.95; 1.0; 1.05; 1.5.

    Завершив моделирование, по виду переходных процессов сделайте вывод о влиянии начальных условий y(0) и параметра р на характер движения системы.