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

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

Введение

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

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

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

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

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

Цель работы

АНАЛИЗ ДИНАМИЧЕСКИХ СИСТЕМ С ЗАПАЗДЫВАНИЕМ

Блок «Идеальное транспортное запаздывание»

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



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

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



где:

T(L, s) – изображение по Лапласу сигнала на выходе из трубопровода;

T(0, s) – изображение по Лапласу сигнала на входе в трубопровод;

τ=L/u – постоянная запаздывания (время транспортировки).

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



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

Сформировать новую структурную схему и задать подписи блокам согласно рисунку (Рисунок 1).

Для формирования структурной схемы необходимо использовать блок «Ступенька» с вкладки «Источники», блок «Идеальное транспортное запаздывание» и два блока «Апериодическое звено 1-го порядка» с вкладки «Динамические», два блока «Мультиплексор» и два блока «Демультиплексор» с вкладки «Векторные», а так же блок «Временной график» с вкладки «Вывод данных».

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



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

Переместить курсор на кнопку «Скрипт» и выполнить нажатие левой кнопкой мыши. Откроется окно редактора «Скрипт страницы». Ввести текст, идентичный примеру (Рисунок 2). Константы «n1» и «n2» будут задавать количество последовательно соединенных инерционных звеньев первого порядка в двух параллельных цепях, аппроксимирующих свойства идеального запаздывающего звена. Закрыть окно «Скрипт страницы» нажатием на кнопку «Закрыть и применить».

Рисунок: Скрипт окна проекта.



В свойствах блока «Ступенька» задать значение свойства «Время срабатывания» равным «2», значение свойства «Начальное значение» равным «0», значение свойства «Конечное значение» равным «1».

В свойствах блока «Идеальное транспортное запаздывание» задать значение «Время запаздывания» равным «2», что означает что данный блок реализует постоянное запаздывание на две секунды. Открыть окно «Свойства» блока «Инерционное звено 1-го порядка» с подписью «8 последовательных звеньев» и заполнить его согласно рисунку (Рисунок 3).

Рисунок: Свойства блока «Инерционное звена 1-го порядка» с подписью «8 последовательных звеньев».



В поле «Формула» свойства «Коэффициенты усиления» введено «n1#1». Это означает, что в данной строке введен числовой вектор из n1 единиц. Можно было ввести данную строку и как вектор-строку. Символ «#» в диалоговых строках эквивалентен предлогу «по» ⇒ n1 элементов по 1.

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

В поле «Формула» свойства «Постоянные времени» задан вектор из n1 одинаковых постоянных времени, равных 2/n1 = 2/8 = 0.25 c.

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

Рисунок: Свойства блока «Инерционное звена 1-го порядка» с подписью «20 последовательных звеньев».



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

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

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



Поскольку алгоритм работы верхнего блока «Инерционное звено 1-го порядка» векторизован, то на вход блока должен поступать векторный сигнал, размерностью n1. Поэтому к скалярному сигналу от блока «Ступенька» необходимо добавить n1-1 сигналов, чтобы после блока «Мультиплексор» векторный сигнал имел размерность n1.

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

Данный пример отображает принцип работы сдвига сигналов в векторе. Необходимо рассмотреть реализацию сдвига, на примере выходного сигнала блока «Ступенька».

Сигнал от блока «Ступенька» поступает на первый элемент вектора входного порта, этот сигнал обрабатывается блоком «Инерционное звено 1-го порядка», первый сигнал вектора с помощью блока «Демультиплексор» подается на место второго элемента вектора в блоке «Мультиплексора», этот сигнал обрабатывается блоком «Инерционное звено 1-го порядка», второй сигнал вектора с помощью блока «Демультиплексор» подается на место третьего элемента в блоке «Мультиплексор» и т.д.

В итоге на втором выходном порте блока «Демультиплексор» будет сигнал, который n1 раз был обработан блоком «Инерционное звено 1-го порядка».

По аналогии с предыдущими блоками, задать свойства блоков «Мультиплексор» и «Демультиплексор» в цепи, реализующей двадцать последовательно соединенных звеньев.

Открыть окно «Параметры проекта» и задать параметры согласно рисунку (Рисунок 6).

Рисунок: Окно «Параметры проекта».



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

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

Рисунок: График переходных процессов.



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

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

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

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

Сформировать новую структурную схему согласно рисунку и задать подписи блокам (Рисунок 8).

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



Для формирования структурной схемы необходимо использовать блок «Полином n-й степени» с вкладки «Источники», блок «Размножитель» с вкладки «Векторные», блок «Язык программирования» с вкладки «Динамические», а также три блока «Фазовый портрет» с вкладки «Вывод данных».

В параметрах расчёта проекта необходимо установить параметры «Минимальный шаг интегрирования» и «Максимальный шаг интегрирования» равным «0.001», а параметр «Конечное время расчёта» равным «3».

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

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

Двойным нажатием левой кнопкой мыши по блоку «Язык программирования» вызвать окно редактора блока. Заполнить окно редактора согласно рисунку (Рисунок 9).

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

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



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

Lm(ω)=20⋅lg(A(ω)),

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

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

φ(ω)=arg(W(iω))

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

FНайквиста=F(Re(W(iω)), Im(W(iω)))

Рисунок: Окно редактора блока «Язык программирования».



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

Для корректности вывода результатов расчета в векторной форме следует использовать блок «Размножитель».

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

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



Рисунок: График ФЧХ.



Рисунок: График ЛАХ.



На рисунках приведены результаты расчета годографов АФЧХ (годографов Найквиста) (Рисунок 10), фазовых частотных характеристик (ФЧХ) (Рисунок 11) и логарифмических амплитудных характеристик (ЛАХ) (Рисунок 12), соответственно. Чёрными линиями представлены характеристики идеального запаздывающего звена, синими линиями – для цепи из восьми звеньев, и красными – для цепи из двадцати звеньев.

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

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

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

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

Демонстрационно-ознакомительная задача в виде структурной схемы САР была рассмотрена в лабораторной работе №1 и представлена на рисунке (Рисунок 13).

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



Объект управления представлен в виде блока «Колебательное звено» с подписью «W2(s)». В блоке «Колебательное звено» свойство «Коэффициент усиления» равно «1», свойство «Постоянная времени» равно «1», свойство «Коэффициент демпфирования» равно «0.5».

Локальная обратная связь представлена в виде блока «Апериодическое звено 1-го порядка» с подписью «W3(s)». В блоке «Апериодическое звено 1-го порядка» свойство «Коэффициент усиления» равно «1», свойство «Постоянная времени» равно «5».

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

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

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



Этапы, которые необходимо выполнить:

  1. Сформировать математическую модель динамики САР с коэффициентом усиления интегрирующего регулятора равным «0.35».
  2. Определить критическое значение постоянной запаздывания τкрит в идеальном запаздывающем звене.
  3. Изменяя значение постоянной запаздывания в идеальном запаздывающем звене в пределах «0.1·τкрит» до «0.9·τкрит» (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 yвых(tk+1) τзапk+1
... ... ... ... ... ...
m tm θm ym yвых(tm) τзапm
... ... ... ... ... ...

При t≤tk значение yвых(t)=y0, а при t>tk значение определяется с использованием данных таблицы (Таблица 1) по алгоритму tm→θm→(θm-1)→y(θm-1)≡yвых(t). Последняя процедура (вычисление yвых(t)) проводится с использованием линейной интерполяции данных таблицы (Таблица 1).

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

  • при t≡tk значение τзап≡tk;
  • при 0<t<tk значение τзап(t) вычисляется по соотношению:

    τзап(tj)=tj+τ0зап⋅(1-θj)

  • при t>tk значение определяется с использованием алгоритма tm→θm(tm)→θm-1→t(θm-1)→τзап(tm)=tm-t(θm-1), причем вычисление t(θm-1) проводится с использованием линейной интерполяции данных таблицы (Таблица 1).

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

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

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

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

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

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

Необходимо сформировать новую структурную схему.

Для формирования структурной схемы необходимо использовать блок «Синусоида» и блок «Кусочно линейная» с вкладки «Источники», блок «Идеальное транспортное запаздывание» и блок «Переменное транспортное запаздывание» с вкладки «Динамические», блок «Запись в список сигналов» с вкладки «Сигналы», а так же блок «Временной график» с вкладки «Вывод данных».

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

Рисунок: Список сигналов проекта.



В свойствах блока «Запись в список сигналов» в поле «Значение» задать свойство «Имена сигналов» равным «U1».

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

В свойствах блока «Синусоида» задать значения свойства «Амплитуда» равным «1», свойства «Частота» равным «0.5», свойства «Сдвиг фазы» равным «0».

В свойствах блока «Кусочно линейная» задать значения свойства «Время» равным «[0, 5, 10, 20, 25, 40]», свойства «Значения функции» равным «[2, 2, 5, 5, 2, 2]».

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



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

  • 0…5 с мгновенное время запаздывания постоянно и равно 2 с;
  • 5…10 с мгновенное время запаздывания линейно растет от 2 с до 5 с;
  • 10…20 с мгновенное время запаздывания постоянно и равно 5 с;
  • 20…25 с мгновенное время запаздывания линейно убывает от 5 с до 2 с;
  • 25…40 с мгновенное время запаздывания постоянно и равно 2 с.

В параметрах расчёта проекта необходимо задать параметр «Время интегрирования» равный «40», параметры «Максимальный шаг интегрирования» и «Минимальный шаг интегрирования» равный «0.01».

Сохранить проект и запустить его на расчёт.

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

Рисунок: График переходных процессов.



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

САМОСТОЯТЕЛЬНОЕ ИССЛЕДОВАНИЕ ИЗВЕСТНЫХ ДИНАМИЧЕСКИХ ЗАДАЧ

МЕТОД СТРУКТУРНОГО МОДЕЛИРОВАНИЯ

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

  1. Уравнением Ван-дер-Поля:

    y``(t)+b⋅[y2(t)-a]⋅y`(t)+y(t)=0

    в диапазоне от t = 0 до t = 40 сек, если y(0) = 1, y`(0) = 0.

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

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

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

  2. Уравнением Матье:

    y``(t)+ε⋅y`(t)+p2⋅[1-2⋅μ⋅cos(2⋅ω⋅t)]⋅y(t)+β⋅y3(t)=0

    в диапазоне от 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)» и параметра «р» на характер движения системы.