Численное решение дифференциально-разностного уравнения первого порядка Текст научной статьи по специальности «Математика»

Системный анализ. Моделирование. Транспорт. Энергетика. Строительство

УДК 517.994 Мухачев Юрий Сергеевич,

к. ф-м. н., доцент кафедры «Теоретической механики и приборостроения», ИрГУПС

тел.: (3952)77-44-63, e-mail: MukhachevJYu_S@mail.ru

Рябов Евгений Валерьевич, научный сотрудник. НИИПФ ИГУ, тел.: (3952)33-21-46 Корчевин Евгений Николаевич, соискатель ИрГУПС, тел.: 89027611628

ЧИСЛЕННОЕ РЕШЕНИЕ ДИФФЕРЕНЦИАЛЬНО-РАЗНОСТНОГО УРАВНЕНИЯ ПЕРВОГО ПОРЯДКА

Yu. S. Mukhachev, E. V. Ryabov, E.N. Korchevin

THE NUMERICAL SOLUTION OF THE DIFFERENTIAL-DIFFERENCE EQUATION

OF THE FIRST ORDER

Аннотация. В статье предлагается метод численного решения характеристического квазиполинома. Предложен алгоритм вычисления корней квазиполинома, исследовано влияние параметров на расположение корней на комплексной плоскости.

Ключевые слова: автоматическое управление, запаздывание, корни квазиполинома.

Abstract. The article proposes a method of numerical solution of a characteristic quasipolynom. The algorithm for calculating the roots is proposed, the parameter influence on the location of the roots on the complex plane is investigated.

Keywords: automatic control, time lag, qua-sipolynomial roots.

Анализ переходных процессов и устойчивости систем автоматического управления (САУ) систем, содержащих в своем составе звенья чистого запаздывания, традиционно считается сложной проблемой теории автоматического управления. Основная сложность задачи связана с тем, что подобные системы описываются дифференциальными уравнениями с запаздывающим аргументом (дифференциально-разностными уравнениями). Характеристическое уравнение системы имеет вид квазиполинома, является трансцендентным и не имеет точного решения [1, 2, 3]. Ситуация значительно осложняется тем, что квазиполином имеет бесконечное количество корней.

Анализ устойчивости управления систем с запаздыванием проводится по частотным критериям Найквиста и Михайлова, для осуществления которых не требуется вычисления точных значений полюсов.

Исследование качества работы САУ в переходных и установившихся режимах являются самостоятельной проблемой, для решения которой необходимо знать точное положение полюсов на комплексной плоскости и оценить влияние параметров системы на их смещение. Приближенные методы анализа основаны на замене звена с запаздыванием обыкновенным апериодическим звеном второго порядка или последовательным соединением нескольких апериодических звеньев [2, 3].

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

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

Передаточная функция исходной системы имеет вид

W (p ) =

Ke

-pi

(1)

Тр + 1

где К - коэффициент передачи,

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

Woc ( p) =

Ke-

W ( p ) =_

1 + W( p) Tp + 1 + Ke-

(2)

Для обратного преобразования Лапласа функцию (2) необходимо представить через полюсы и вычеты в виде [4]

кс(р)=Е

с

(р - Рк):

(3)

То + 1 + Ке сою = 0, Тю - Ке 81пют = 0 .

(7)

(8)

Преобразуем формулу (8) с целью разделения переменных следующим образом:

Та 1

То + 1 + Кесою = 0 ^ соБют =--. (10)

Ке-°

После элементарных преобразований (10) получим выражение для ют в виде

Го + \'

ют = агссоБ -

Ке-

(11)

Преобразуем формулу (8) с целью разделения переменных

Тю - Ке 0 81пют = 0 ^ о = -Тю

т I Кз1пют

. (12)

Подставим выражение для переменной а из формулы (12) в формулу (11) и получим неявное уравнение относительно переменной ю в следующем виде:

со8ют + -

81п ют

ют

т , I т „ | , I 81Пют - + Ц- К| + 1п1-

Т IТ ) I ют

0. (13)

где С - вычеты, р^ - полюсы.

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

Гр + \ + Ке-р = 0. (4)

К сожалению, общего метода решения уравнений типа (4) не существует.

Предпримем попытку решения уравнения (4) приближенными численными методами. Рассмотрим квазиполином более подробно. Представим комплексную переменную р в виде

р = о + ]ю, (5)

где о - действительная часть, ю - мнимая часть, ] - мнимая единица.

Подставив выражение (5) в квазиполином (4) и выполнив ряд элементарных преобразований, получим следующее уравнение:

(То + 1 + Кесою)+ ](гю - Ке81пют)= 0. (6)

Комплексная величина равна нулю в том случае, если равны нулю действительная и мнимая величины, поэтому уравнение (6) разделяется на два уравнения вида

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

Численное решение уравнения (13) требует предварительной оценки возможного положения нулей, поэтому введем вспомогательную переменную

2 (ю) = СОБют +

БШют

ют

и 1пI ^К

1п

Б1П ют

(14)

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

Рис. 1. Вспомогательная функция /(ю)

Алгоритм поиска корней состоит из следующих действий:

а) задается область значений аргумента ют, в которой необходимо найти корни;

б) производится предварительное вычисление значений 2(ю) с грубым шагом по координате ют;

в) выделяется интервал, ограниченный двумя соседними значениями ют, внутри которого переменная 2(ю) меняет знак;

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

д) в найденном узком интервале находится приближенное значение корня методом секущей;

к

Системный анализ. Моделирование. Транспорт. Энергетика. Строительство

е) найденные значения корней на оси юх используются для вычисления переменной а по формуле (12);

ж) для полученных приближенных значений производится расчет абсолютной и относительной погрешностей вычисленных значений полюса;

з) полученные значения а и ю объединяются в комплексную переменную р = а + , которая изображается на комплексной плоскости;

и) численные значения полюсов и относительных погрешностей выводятся на монитор.

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

1) основную программу, задающую исходные значения и обрабатывающую результаты численного решения,

2) подпрограмму непосредственного вычисления корней квазиполинома.

На рис. 2 показан пример поиска решений характеристического квазиполинома для постоянной времени Т = 3,0, времени запаздывания п =1,0, коэффициента передачи К = 6,0.

-2 -1 а, а.е.

Рис. 2. Пример расчета полюсов квазиполинома

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

На рис. 3 показана зависимость вспомогательной функции (18) от юх, рассчитанная при разных значениях коэффициента усиления К. Из рисунка видно, что зависимость от коэффициента усиления имеет две принципиально разных области, разделенные критическим значением ККР = 0,79.

1.5

0.5

-0.5

-1.5

-1

\1к=3,о1

/ / \

/ ик=0,79| \

/ \ Ь{К=0 5] 1

у /// 1 1

V у

-10

0

о, а.е.

10

15

Рис. 3. Вспомогательная функция при коэффициенте усиления выше критического (К=3), ниже критического (К=0,5) и равном критическому (К=0,79)

При значениях выше критического коэффициента, т. е. К > ККР, в интервале юх ={-п, +п} существует два комплексно сопряженных полюса, ниже критического коэффициента при К < ККР в этом интервале вообще не существует полюсов. При точном равенстве К = ККР существует один действительный полюс.

Величину критического коэффициента усиления ККР можно оценить по формуле (13), положив юх =0, тогда получим выражение

Т

1 + -- 1п

Т

\хККР

= 0,

(15)

из которого следует решение относительно ККР в виде

К

1+1+1п 1

т т

(16)

При условиях для Т = 3 и п = 1 формула (16) дает точное значение для критического коэффициента усиления ККР = 0,79079. Зависимость положения полюсов на комплексной плоскости и при варьировании коэффициента усиления выше и ниже критического усиления показана на рис. 4.

Из рисунка 4 видно, что принципиальным моментом является наличие критического коэффициента ККР.

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

На рисунке 6 показана зависимость положения полюсов от времени запаздывания. Из рисунка видно, что время задержки оказывает значитель-

е

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

Рис. 4. Влияние коэффициента усиления на расположение полюсов

Рис. 5. Влияние постоянной времени на расположение полюсов

60

40

20

-20

-40

-60.

- т-0,5 - \

-8 -6 -4-2 0 2

о, а.е.

Рис. 6. Влияние времени задержки на расположение полюсов

На рисунке 3 показана зависимость вспомогательной функции (18) от ют, рассчитанная при разных значениях коэффициента усиления К. Из рисунка видно, что зависимость от коэффициента усиления имеет две принципиально разных области, разделенные критическим значением ККР = 0,79.

Заключение

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

2. Зависимость расположения полюсов от коэффициента усиления указывает на наличие критического коэффициента усиления:

- выше критического коэффициента наблюдается пара комплексно сопряженных полюсов;

- при критическом коэффициенте наблюдается один действительный полюс в отрицательной полуплоскости;

- если коэффициент усиления ниже критического, то низкочастотные полюса при ют ={-п , +п} отсутствуют.

3. Зависимость расположения полюсов от постоянной времени указывает, что постоянная времени существенно влияет только на параметр затухания, почти не изменяя частоту колебаний.

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

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

БИБЛИОГРАФИЧЕСКИЙ СПИСОК

1. Методы классической и современной теории автоматического управления : учеб. : в 5 т. / под ред. К.А. Пупкова, Н.Д. Егупова. - 2-е изд., перераб. и доп. М. : Изд-во МГТУ им. Н. Э. Баумана, 2004. Т. 1 : Математические модели, динамические характеристики и анализ систем автоматического управления. 635 с.

2. Ерофеев А. А. Теория автоматического управления : учеб. для вузов. 2-е изд., перераб. и доп. СПб. : Политехника, 2005. 302 с.

3. Бесекерский В. А., Попов Е. П. Теория систем автоматического управления. Изд. 4-е, пере-раб. и доп. СПб. : Профессия, 2004. 752 с. (Специалист).

4. Беллман Р., Кук К. Дифференциально-разностные уравнения : пер. с англ. М. : МИР, 1967. 548 с.