Понятие о качественном анализе динамических систем. Качественный анализ динамических систем. Построение фазовых портретов ДС. Построение матрицы Коши

Автоматика и телемеханика, Л- 1, 2007

РАС Б 02.70.-c, 47.ll.-j

© 2007 г. Ю.С. ПОПКОВ, д-р техн. наук (Институт системного анализа РАН, Москва)

КАЧЕСТВЕННЫЙ АНАЛИЗ ДИНАМИЧЕСКИХ СИСТЕМ С Вд-ЭНТРОПИЙНЫМ ОПЕРАТОРОМ

Предлагается метод исследования существования, единственности и локализации сингулярных точек рассматриваемого класса ДСЭО. Получены условия устойчивости "в малом" и "в большом". Приводятся примеры применения полученных условий.

1. Введение

Многие проблемы математического моделирования динамических процессов могут быть решены на основе концепции динамических систем с энтропийным оператором (ДСЭО) . ДСЭО представляет собой динамическую систему, в которой нелинейность описывается параметрической задачей максимизации энтропии. Феио-моиологически ДСЭО является моделью макросистемы с "медленным" самовоспроизведением и "быстрым" распределением ресурсов . Некоторые свойства ДСЭО исследовали в . Настоящая работа продолжает цикл исследований качественных свойств ДСЭО.

Рассматривается динамическая система с Вд-энтропийным оператором :

^ = £{х,у{х)), х е Еп:

у(х) = а^шах{Нв (у) | Ту = ц(х), у е Е^} > 0.

В этих выражениях:

С(х,у), ц(х) - непрерывно дифференцируемые вектор-функции;

Энтропия

(1.2) Нв (у) = уз 1п аз > 0, з = Т~т;

Т - (г х ш)-матрица с элементами ^ 0 имеет полный ранг, равный г;

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

(1.3) Q = {ц: 0 <оТ ^ ц ^ а+} С Е+,

где а- и а + - векторы из Е+, прпчем а- - вектор с малыми компонентами.

Воспользовавшись известным представлением энтропийного оператора через множители Лагранжа . преобразуем систему (1.1) к следующему виду:

- = £(х,у(г)), х е Кп, у(г) е К?, г е Ег+

Уз (г) = аз\\ ^, 3 = 1,т-

О(х,г) = Ту(г) = д(х),

где гк = ехр(-Ак) > 0 - экспоненциальные множители Лагранжа.

Наряду с ДСЭО общего вида (1.1) будем рассматривать, следуя классификации, приведенной в .

ДСЭО с сепарабельным потоком:

(1-5) ^ = I (х) + Ву(г),

где В (п х т)-матрица;

ДСЭО с мультипликативным потоком:

(1.6) ^ = х ® (а - х ® Шу(г)), аЬ

где Ш - (п х т)-матрица с неотрицательными элементами, а - вектор с положительными компонентами, ® - знак покоординатного умножения.

Задача данной работы состоит в исследовании существования, единственности и локализации сингулярных точек ДСЭО и их устойчивости.

2. Сингулярные точки

2.1. Существование

Рассмотрим систему (1.4). Сингулярные точки этой динамической системы определяются следующими уравнениями:

(2.1) С^(х,у(г))=0, г = ТП;

(2.2) уз (г) = а^ г^, 3 = Т^:

(2.3) вк (г) = ^ аз г^ = дк (х), к =1,г.

Рассмотрим вначале вспомогательную систему уравнений:

(2.4) С(д,г) = г, д е Я,

где множество Я определено равенством (1.3) и С(д,г) - вектор-функция с компонентами

(2.5) Ск(д,г) = - Ок(г), а- < дк < а+, к =1,г.

Уравнение (2.4) имеет единственное решение г* при каждом фиксированном векторе д, что следует из свойств Вд-энтропийного оператора (см. ).

Из определения компонент вектор-функции С(д,г) имеет место очевидная оценка:

(2.6) С(а+,г) < С(д, г) < С(а-,г), г в Е+. Рассмотрим два уравнения:

Обозначим решение первого уравнения через г+ и второго - через г-. Определим

(2.7) C (a+,z) = z, C(a

(2.8) zmaX = max z+, zmin = mm zk

и r-мерные векторы

(2.9) z {zmax, zmax}, z {zmin , zmin}.

Лемма 2.1. Для всех q G Q (1 . 3) решения z*(q) уравнения (2.4) принадлежат, вектор 1 юму отрезку

zmin < z*(q) < zmax,

где векторы zmin и zmax определяются выражениями (2.7)-(2.9).

Доказательство теоремы приведено в Приложении. Qq

ций qk(x) (1.3) для x G Rn, то имеет место

Следствие 2.1. Пусть выполнены условия леммы 2.1 и функции qk(x) удовлетворяют условиям (1.3) для вс ex x G Rn. Тогда для в сех x G Rm решен ия z* уравнения (2.3) принадлежат векторному отрезку

zmin < z* < zmax

Вернемся теперь к уравнениям (2.2). которые определяют компоненты вектор-функции y(z). Элементы ее якобиана имеют вид

(2.10) jb aj zk JJ & > 0

для всех z G R+ за исключением 0 и ж. Следовательно, вектор-функция y(z) строго монотонно-возрастающая. Согласно лемме 2.1 она ограничена снизу и сверху, т.е. для всех z G Rr (следовательно, для всех x G Rn) ее значения принадлежат множеству

(2.11) Y = {у: у- < y < y+},

где компоненты векторов yk, y+ определяются выражениями:

(2.12) yk = aj y+ = aj znlax, j = h™.

(2.13) bj = Y, tsj, 3 =1,

Рассмотрим первое уравнение в (2.1) и перепишем его в виде:

(2.14) L(x,y) = 0 для всех у е Y С Е^.

Это уравнение определяет зависимость переменной x от переменной у, принадле-Y

мы (1.4) сводится к существованию неявной функции x(y), определяемой уравнением (2.14).

Лемма 2.2. Пусть выполняются следующие условия:

а) вектор-функция L(x,y) непрерывна по совокупности переменных;

б) lim L(x, у) = ±<ж для любого фиксированного у е Y;

в) det J (x, у) = 0 для вс ex x е Еп для любого фиксиров анн ого у е Y.

Тогда существует единственная неявная функция x*(у), определенная на Y. В этой лемме J(x, у) - якобиан с элементами

(2.15) Ji,i (x,y) = --i, i,l = l,n.

Доказательство приведено в Приложении. Из приведенных лемм следует

Теорема 2.1. Пусть выполнены условия лемм 2.1 и 2.2. Тогда существует единственная сингулярная точка ДСЭО (1.4) и, соответственно, (1.1).

2.2. Локализация

Под исследованием локализации сингулярной точки понимается возможность установления интервала, в котором она находится. Задача эта весьма не простая, но для некоторого класса ДСЭО такой интервал можно установить.

Обратимся к первой группе уравнений в (2.1) и представим их в виде

(2.16) L(x,y)=0, у- й у й у+,

где у- и у+ определяются равенствами (2.12), (2.13).

Теорема 2.2. Пусть вектор-функция L(x,y) непрерывно дифференцируема и монотонно-возрастающая по обеим переменным, т.е.

-- > 0, -- > 0; i,l = 1, n; j = 1,m. dxi dyj

Тогда решение системы (2.16) по переменной x принадлежит интервалу (2.17) xmin й x й xmax,

а) векторы xmin, xmax имеют вид

Min = i x 1 xmax = r x т;

\xmin: . .., xminlxmax, . . ., xmax} :

xmin - ^Qin ^ ■ , xmax - ^QaX ^ ;

6) x- и x+ - компоненты решении следующих уравнений

(2.19) L(x,y-)=0, L(x,y+) = 0

с oo m в етственно.

Доказательство теоремы приведено в Приложении.

3. Устойчивость ДСЭО "в малом"

3.1. ДСЭО с сепарабельпым потоком Обратимся к уравнениям ДСЭО с сепарабельпым потоком, представив их в виде:

- = /(х) + Бу(г(х)), х е Кп аЬ

У- (г(Х)) = азП (Х)У33 , 3 = 1,"~ 8 = 1

0(х, г(х)) = Ту(г(х)) = д(х), г е Нг,.

Здесь значения компонент вектор-функции д(х) принадлежат множеству Q (1.3), (п х ш)-матрица Б имеет полный ранг, равный п (п < ш). Вектор-функция / (х) непрерывно дифференцируемая.

Пусть рассматриваемая система имеет сингулярную точку ж. Для исследования устойчивости этой сингулярной точки!"в малом" построим линеаризованную систему

где А - (п х п)-матрица, элементы которой вычислены в точке ж, и вектор £ = х - ж. Согласно первому уравнению в (3.1) матрица линеаризованной системы имеет

А = 7 (х) + БУг (ж)Их(х), х = г(х),

| 3 = 1,ш,к = 1,

I к = 1,г, I = 1,п

Из (3.1) определяются элементы матрицы Уг: ду.

"Ькз П " 8=1

3 , г8 х8 , 5 1,г.

Для определения элементов матрицы Zx обратимся к последней группе уравнений в (3.1). В показано, что данные уравнения определяют неявную вектор-функцию г(х), которая непрерывно-дифференцируема, если вектор-функция д(х) непрерывно дифференцируема. Якобиан Zx вектор-функции г(х) определяется уравнением

<Эг (z)Zx(Х) = Qx(Х),

вг (Х) = Т Уг (X),

ддк, -т- , -" -- к = 1,г, I = 1,п дх\

Из этого уравнения имеем (3.9) Zx(x) = в-1(z)Qx(x).

Подставляя этот результат в равенство (3.3). получим:

А = 1 (х) + Р (х), Р (х) = ВУг (г)[ТУг (г)]-1 Qx(x).

Таким образом, уравнение линеаризованной системы приобретает вид

(з.и) | = (j+р)е

Здесь элементы матриц J, Р вычислены в сингулярной точке. Достаточные условия устойчивости "в малом" ДСЭО (3.1) определяет следующая

Теорема 3.1. ДСЭО (3.1) имеет устойчивую "в малом" сингулярную точку x, если выполняются следующие условия:

а) матрицы J, Р (3.10) линеаршованной системы (3.11) имеют вещественные и различные собственные числа, причем матрица J имеет максимальное собственное число

Птах = max Пг > 0,

Wmax = max Ui < 0;

Umax + Птах <

Из этой теоремы и равенства (3.10) следует, что для сингулярных точек, для которых Qx(x) = 0 и (или) для X, = 0 щи tkj ^ 1 для всex k,j достаточные условия теоремы не выполняются.

3.2. ДСЭО с мультипликативным потоком Рассмотрим урвиения (1.6). представив их в виде:

X ® (a - x ® Wy(z(x))), x е Rn;

yj(z(x)) = aj ПZs(x)]isi" j = 1,m;

(ЗЛ2) yj (z(x)) = a^ <~"ts

Q(x, z(x)) = Ty(z(x)) = q(x), z е R++.

системы. Будем иметь:

(3.13) А = ^ [см] - 2ХШУх (г^х(х).

В этом выражении diag С] - диагональныя матрица с положительными элементами а1,..., ап, Уг, Zx - матрицы, определенные равенствами (3.4)-(3.7).

Представим матрицу A в виде

(3.14) A = diag+P (x),

(3.15) P (x) = -2xWYz (z)Zx(x).

Обозначим: maxi ai = nmax и wmax - максимадьное собственное число матрицы P(x) (3.15). Тогда теорема 3.1 справедлива и для ДСЭО (1.6). (3.12).

4. Устойчивость ДСЭО "в большом"

Обратимся к уравнениям ДЭСО (1.4), в которых значения компонент вектор-функции q(x) принадлежат множеству Q (1.3). В рассматриваемой системе существует сингулярная точка Z, которой соответствуют векторы z(x) = z ^ z- > 0 и

y(x) = y(z) = у > у- > 0.

Введем векторы отклонений £, C, П от сингулярной точки: (4.1) £ = x - x, (= у - у, п = z - z.

ЖЕЖЕРУН А.А., ПОКРОВСКИЙ А.В. - 2009 г.

1

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

двоичная динамическая система

динамическое свойство

качественный анализ

булевы ограничения

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

1. Biere A., Ganesh V., Grohe M., Nordstrom J., Williams R. Theory and Practice of SAT Solving. Dagstuhl Reports. 2015. vol. 5. no. 4. Р. 98–122.

2. Marin P., Pulina L., Giunchiglia E., Narizzano M., Tacchella A. Twelve Years of QBF Evaluations: QSAT Is PSPACE-Hard and It Shows. Fundam. Inform. 2016. vol. 149. Р. 133–58.

3. Бохман Д., Постхоф Х. Двоичные динамические системы. М.: Энергоатомиздат, 1986. 400 с.

4. Маслов С.Ю. Теория дедуктивных систем и ее применение. М.: Радио и связь, 1986. 133 с.

5. Jhala R., Majumdar R. Software model checking. ACM Computing Surveys. 2009. vol. 41. no. 4. Р. 21:1–21:54.

6. Васильев С.Н. Метод редукции и качественный анализ динамических систем. I–II // Известия РАН. Теория и системы управления. 2006. № 1. С. 21–29. № 2. С. 5–17.

7. DIMACS format [Электронный ресурс]. Режим доступа: http://www.cs.utexas.edu/users/moore/acl2/manuals/current/manual/index-seo.php/SATLINK____DIMACS (дата обращения: 24.07.2018).

8. QDIMACS standard [Электронный ресурс]. Режим доступа: http://qbflib.org/qdimacs.html (дата обращения: 24.07.2018).

9. Delgado-Eckert E., Reger J., Schmidt K. Discrete Time Systems with Event-Based Dynamics: Recent Developments in Analysis and Synthesis Methods. Mario Alberto Jordan (Ed.). Discrete Time Systems. InTech. 2011. Р. 447–476.

10. Васильев С.Н. Достижимость и связность в автоматной сети с общим правилом переключения состояний // Дифференциальные уравнения. 2002. Т. 38. № 11. С. 1533–1539.

11. Бычков И.В., Опарин Г.А., Богданова В.Г., Горский С.А., Пашинин А.А. Мультиагентная технология автоматизации параллельного решения булевых уравнений в распределенной вычислительной среде // Вычислительные технологии. 2016. Т. 21. № 3. С. 5–17.

12. Lonsing F., Biere A. DepQBF. A Dependency-Aware QBF solver. Journal on Satisfiability. Boolean Modeling and Computation. 2010. vol. 9. Р. 71–76.

13. Oparin G.A., Bogdanova V.G., Pashinin A.A., Gorsky S.A. Distributed Solvers of Applied Problems Based on Microservices and Agent Networks. Proc. Of the 41th Intern. Convention on Information and Communication Technology, Electronics and Microelectronics (MIPRO-2018). Р. 1643–1648.

14. Bogdanova V.G., Gorsky S.A. Scalable Parallel Solver of Boolean Satisfiability Problems. Proc. Of the 41th Intern. Convention on Information and Communication Technology. Electronics and Microelectronics (MIPRO-2018). Р. 244–249.

15. Bychkov I.V., Oparin G.A., Bogdanova V.G., Pashinin A.A. The Applied Problems Solving Technology Based on Distributed Computational Subject Domain Model: a Decentralized Approach // Параллельные вычислительные технологии XII международная конференция, ПаВТ’2018, г. Ростов-на-Дону, 2–6 апреля 2018 г. Короткие статьи и описания плакатов. Челябинск: Издательский центр ЮУрГУ, 2018. C. 34–48.

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

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

Для ДДС, функционирование которых рассматривается на конечном интервале времени, такие ограничения являются булевыми и записываются на языке булевых уравнений или булевых формул с кванторами. Первый тип ограничений приводит к необходимости решения SAT-задачи (задачи булевой выполнимости ); второй тип ограничений связан с решением задачи TQBF (проверки истинности квантифицированных булевых формул ). Первая задача является типичным представителем класса сложности NP, а вторая задача - класса сложности PSPACE. Как известно, PSPACE-полнота дискретной задачи дает более сильное свидетельство о ее труднорешаемости, чем NP-полнота. В силу этого сведение задачи качественного анализа ДДС к SAT-задаче более предпочтительно, чем сведение к задаче TQBF. В общем случае исследование не каждого свойства ДДС можно представить на языке булевых уравнений.

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

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

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

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

1) разработка ориентированного на специалиста по динамике систем логического языка спецификации динамических свойств;

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

3) представление полученной модели в международном формате DIMACS или QDIMACS ;

4) выбор (разработка) эффективного параллельного (распределенного) решателя задачи выполнимости булевых ограничений (SAT или TQBF решателя);

5) разработка инструментальных средств создания программных сервисов;

6) разработка сервисов для качественного исследования разнообразных динамических свойств ДДС.

Целью настоящего исследования является решение только первых двух задач применительно к алгоритмизации качественных исследований автономных (без управляющих входов) синхронных ДДС. Такие системы в англоязычных публикациях принято называть синхронными булевыми сетями (Boolean network). Другие аспекты применения метода булевых ограничений (в том числе и для ДДС с управляющими входами) являются предметом следующих публикаций.

Математическая модель автономной ДДС

Пусть X = Bn (B = {0, 1} - множество двоичных векторов размерности n (пространство состояний ДДС). Через t∈T = {1,…,k} обозначим дискретное время (номер такта).

Для каждого состояния x0∈X, называемого начальным состоянием, определим траекторию x(t, x0) как конечную последовательность состояний x0, x1,…, xk из множества X. Далее будем рассматривать ДДС, в которых каждая пара смежных состояний xt, x(t - 1) (t∈T) траектории связана отношением

xt = F(xt - 1). (1)

Здесь F:X>X - векторная функция алгебры логики, называемая функцией переходов. Таким образом, для любых x0∈X система булевых уравнений (1) представляет модель динамики поведения траекторий ДДС в пространстве состояний X на конечном интервале времени T = {1, 2,…,k}. Здесь и далее величина k в определении множества T предполагается наперед заданной постоянной. Такое ограничение является вполне естественным. Дело в том, что при качественном анализе поведения траекторий ДДС практический интерес представляет вопрос о том, что можно сказать о выполнимости какого-либо динамического свойства при фиксированном, не слишком большом k. Выбор величины k в каждом конкретном случае осуществляется исходя из априорных сведений о длительности протекания процессов в моделируемой дискретной системе.

Известно , что система булевых уравнений (1) с начальным состоянием x0∈X для T = {1, 2,…,k} эквивалентна одному булевому уравнению вида

При k = 1 (рассматриваются только одношаговые переходы) уравнение (2) приобретает вид

(3)

Решения этого уравнения определяют направленный граф, состоящий из 2n вершин, отмеченных одним из 2n состояний множества X. Вершины x0 и x1 графа соединены дугой, направленной от состояния x0 к состоянию x1. Такой граф в теории двоичных автоматов называется диаграммой переходов. Представление поведения ДДС в виде диаграммы переходов весьма наглядно как при построении траекторий, так и исследовании их свойств, но практически реализуемо лишь для небольших размерностей n вектора состояния x∈X.

Языковые средства спецификации динамических свойств

Наиболее удобно задавать спецификацию динамического свойства на языке формальной логики. Следуя работе , обозначим через X0∈X, X1∈X, X*∈X - множества начальных, допустимых и целевых состояний.

Основными синтаксическими элементами логической формулы динамического свойства являются: 1) предметные переменные (компоненты векторов x0, x1,…, xk, время t); 2) ограниченные кванторы существования и всеобщности; 3) логические связки v, &; заключительные формулы. Заключительная формула представляет утверждение о принадлежности некоторых состояний множества траекторий x(t, x0) (x0∈X0) оценочным множествам X* и X1.

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

где A(y) - предикат, ограничивающий значение переменной y.

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

В дальнейшем будем предполагать, что элементы множеств X0, X1, X* определяются соответственно нулями следующих булевых уравнений

или характеристическими функциями этих множеств , .

С учетом ограничения на начальные состояния G0(x) = 0 наряду с уравнениями (2, 3) для сокращения записи будем использовать следующие булевы уравнения:

(4)

Предварительный качественный анализ автономных ДДС

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

Состояние x1 в (3) будем называть последователем состояния x0, а x0 - предшественником состояния x1. В автономной ДДС каждое состояние имеет только одного последователя, а число предшественников данного состояния может изменяться от нуля до 2n - 1. Все непосредственные предшественники x0 состояния s∈X являются нулями булева уравнения

Если уравнение (6) не имеет решений, то предшественники состояния s отсутствуют.

Равновесные состояния (если они существуют) являются решениями булева уравнения

Траектория x0, x1,…, xk называется циклом длины k, если состояния x0, x1,…, xk-1 попарно отличны друг от друга и xk = x0. Циклическая последовательность длины k (если она существует) является решением булева уравнения

где = 0 () - условия попарного различия множества состояний C цикла длины k. Если ни одно из состояний цикла не имеет предшественников, не принадлежащих множеству C, то такой цикл называется изолированным. Пусть элементы s множества C определяются решением булева уравнения Gc(s) = 0. Тогда несложно показать, что условие изолированности цикла эквивалентно отсутствию нулей у следующего булева уравнения:

Решения уравнения (7) (если они существуют) определяют состояния цикла, имеющие предшественников, не принадлежащих множеству C.

Так как равновесное состояние является циклом длины k = 1, то условие его изолированности аналогично условию изолированности с k ≥ 2 за тем отличием, что Gc(s) имеет вид полной дизъюнкции, определяющей это равновесное состояние.

Неизолированные равновесные состояния и циклы в дальнейшем будем называть аттракторами.

Спецификация динамических свойств типа достижимости

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

Определение этого свойства связано с заданием введенных ранее множеств X0, X*, X1 (соответствующих этим множествам булевых уравнений). Предполагается, что для множеств X0, X*, X1 выполняется ограничение

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

1. Основное свойство достижимости множества X* из множества X0 формулируется следующим образом: любая траектория, выпущенная из множества начальных состояний X0, достигает целевого множества X*. С использованием ограниченных кванторов существования и всеобщности, формула этого свойства имеет вид:

2. Свойство безопасности обеспечивает для любой траектории, выпущенной из X0, недостижимость множества X*:

3. Свойство одновременной достижимости. В ряде случаев может выставляться более «жесткое требование», которое состоит в том, чтобы каждая траектория достигала целевого множества ровно за k тактов (k∈T):

4. Свойство достижимости при фазовых ограничениях:

Это свойство гарантирует, что все траектории, выпускаемые из множества X0, до момента попадания в целевое множество X* находятся в множестве X1.

5. Свойство притяжения. Пусть X* является аттрактором. Тогда логическая формула свойства притяжения совпадает с формулой основного свойства достижимости:

т.е. для каждой траектории, выпущенной из множества X0, существует момент времени t∈T, начиная с которого траектория не выходит за пределы множества X*. Множество X0 в этом случае принадлежит части области притяжения множества X*(X0∈Xa, где Xа - полная область притяжения (бассейн) аттрактора).

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

Приведем два свойства, проверка которых приводит к необходимости решения задачи TQBF.

6. Свойство связности целевого множества:

т.е. существует начальное состояние x0∈X0 такое, что каждое целевое состояние x*⊆X* в некоторый момент t∈T достижимо, что означает существование соответствующей этому состоянию траектории, такой, что все целевые состояния x*∈X* принадлежат данной траектории.

7. Свойство тотальной достижимости множества X* из X0:

т.е. каждое целевое состояние является достижимым из X0.

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

Для свойств (1-5) проверка их выполнимости сводится к поиску нулей булева уравнения, технология формирования которого носит стандартизованный характер и подробно рассматривается только для основного свойства достижимости. Свойства (6, 7) приводят к задаче проверки истинности квантифицированной булевой формулы.

1. Основное свойство достижимости. Его логическая формула имеет вид

С учетом (4) запишем формулу (8) в виде

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

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

Формула (10) истинна тогда и только тогда, когда тождественно истинно подкванторное выражение, т.е.

Тождественная истинность импликации означает, что булева функция является логическим следствием функции , т.е. любая траектория с начальным состоянием x0∈X0 достигает целевого множества X*.

Выполнимость тождества (11) эквивалентна отсутствию нулей у булева уравнения

При получении (12) мы избавились от импликации и заменили ϕ*(x0, x1,..., xk) на . Если уравнение (12) имеет хотя бы одно решение, то свойство достижимости не имеет места. Такое решение представляет (в определенном смысле) контрпример для проверяемого свойства и может помочь исследователю выявить причину возникновения ошибки.

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

2. Свойство безопасности

3. Свойство одновременной достижимости

4. Свойство достижимости при фазовых ограничениях

5. Свойство притяжения. Выполнимость этого свойства проверяется в два этапа. На первом этапе выясняется, является ли множество X* аттрактором. Если ответ положительный, то на втором этапе проверяется основное свойство достижимости. Если X* достижимо из X0, то все условия свойства притяжения выполнены.

6. Свойство связности

7. Свойство тотальной достижимости`

Для свойств (6, 7) скалярная форма равенства двух булевых векторов xt = x* имеет вид

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

Обозначим через x0∈X = B3 начальное состояние модели (13). Пусть T = {1, 2}. Выпишем требуемые для спецификации свойств функции одношагового и двухшагового переходов модели (13):

(14)

где знаком «.» обозначена операция конъюнкции.

Для проверки выполнимости каждого свойства задаются начальное (X0) и целевое (X*) множества, определяемые нулями уравнений G0(x) = 0, G*(x) = 0 или характеристическими функциями этих множеств (см. п. 2). В качестве SAT-решателя используется решатель инструментального комплекса (ИК) РЕБУС , а TQBF-решателя - DepQBF . Кодировка переменных в булевых моделях рассматриваемых ниже свойств для этих решателей приведена в табл. 1, булевы модели этих свойств в форматах DIMACS и QDIMACS расположены в табл. 2.

Таблица 1

Кодировка переменных

Номер переменной в булевой модели

Свойство 1

Свойство 2

Свойство 3

Свойство 4

Свойство 5

Таблица 2

Булевы модели свойств

Свойство 1	Свойство 2	Свойство 3	Свойство 4 (А)	Свойство 4 (B)	Свойство 5
					e 1 2 3 4 5 6 7 8 9 0 4 -5 -6 7 -8 -9 -10 11 12 0 4 5 6 -7 8 9 10 -11 -12 0

1. Основное свойство достижимости (k = 2). Пусть X0 = {x∈X: x1 = 0}, X*={x∈X: x1 = 1}. Начальное и целевое множества определяются соответственно уравнениями G0(x) = x1 = 0 и . Булево уравнение (12) в этом случае приобретает вид

где функция ϕ(x0, x1, x2) определена в (14). Решатель ИК РЕБУС выдает ответ «unsat» (уравнение не имеет нулей), таким образом свойство достижимости X* из X0 выполняется, что наглядно видно из следующей диаграммы переходов, приведенной на рисунке.

2. Циклы длины k = 2. Циклическая последовательность длины 2 (если она существует) является решением булева уравнения

Функция имеет вид

Выражение R(x0, x1) при нахождении цикла не включалось в уравнение, так как циклы длины k = 1 (равновесные состояния) в модели (13) отсутствуют. С помощью решателя ИК РЕБУС получены два ответа (в выходном формате DIMACS): 1 2 3 4 5 -6 0 и 1 2 -3 4 5 6 0, соответствующие циклическим последовательностям (рисунок): {(1 1 1), (1 1 0)} и {(1 1 0), (1 1 1)}. Множества состояний обоих циклов совпадают, что означает наличие в модели (13) одного цикла длины k = 2.

Диаграмма переходов системы (13)

3. Свойство изолированности цикла. Если элементы s множества состояний C цикла длины k = 2 определяются решением булева уравнения Gc(s) = 0, то условие изолированности цикла эквивалентно отсутствию нулей у следующего булева уравнения:

Так как C = {(1 1 1), (1 1 0)}, имеем

Для этого уравнения решатель ИК РЕБУС находит два решения: -1 2 3 4 5 -6 0 и -1 2 -3 4 5 -6 0 (в двоичном представлении согласно кодировке переменных в табл. 1 это пары состояний (0 1 1), (1 1 0) и {(0 1 0), (1 1 0)). Таким образом, состояние цикла (1 1 0) имеет два предшественника, (0 1 1) и (0 1 0), не принадлежащих множеству состояний цикла. Это означает, что свойство изолированности цикла не выполняется, т.е. данный цикл является аттрактором.

4. Свойство притяжения. Пусть X* = C является аттрактором. Логическая формула свойства притяжения совпадает с формулой основного свойства достижимости

а соответствующее булево уравнение для нашего случая имеет вид

Выпишем функции G0(x0), ϕ(x0, x1, x2) и . Функция ϕ(x0, x1, x2) приведена в (14). Для X* = C выражение равно . Рассмотрим два варианта задания множества начальных состояний X0, для случаев выполнения (А) и невыполнения (B) свойства притяжения за k = 2 тактов.

A. Пусть . Тогда

В этом случае для булевого уравнения (15) выдается ответ «unsat». Свойство притяжения для заданного множества X0 выполняется.

B. Пусть . Тогда

В этом случае ИК РЕБУС для уравнения (15) находит решение: 1 -2 3 4 -5 -6 -7 8 9 0, которое соответствует траектории {(1 0 1),(1 0 0),(0 1 1)}. Эта траектория с начальным состоянием x0 = (1 0 1) за два такта не достигает множества X* = C, что означает невыполнимость свойства притяжения для заданного X0.

5. Свойство связности. Логическая формула свойства связности имеет вид следующего высказывания:

Для k = 2 ϕ*(x0, x1, x2) = G0(x0)∨ϕ(x0, x1, x2), где функция ϕ(x0, x1, x2) приведена в (14). В качестве начального выберем состояние (1 0 1). Тогда . Пусть целевое множество X* = {(0 1 1), (1 0 0)}. В этом случае функция G*(x*) имеет вид

Запишем G*(x*) в формате КНФ:

Используя закон Де-Моргана, найдем отрицание функции ϕ*(x0, x1, x2). Подставив в (16) все полученные функции и с учетом кодировки булевых переменных (табл. 1), получим булеву модель в формате QDIMACS (табл. 2). Решатель DepQBF выдает ответ «sat», что означает истинность высказывания (16). Свойство связности для заданных X0, X*, T = {1, 2} выполнено.

Заключение

К основным достоинствам метода булевых ограничений в качественном исследовании ДДС относятся:

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

2. По формуле свойства и уравнениям динамики автоматически выполняется построение соответствующего булева уравнения или квантифицированной булевой формулы.

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

4. Проблема сокращения перебора в той или иной мере решается разработчиками этих решателей и экранирована от специалистов по качественному анализу ДДС.

5. Обеспечивается возможность решения задачи качественного анализа ДДС для больших размерностей вектора состояния n на достаточно продолжительном промежутке времени T. По числу состояний метод булевых ограничений количественно соизмерим с методом model checking. В силу того, что в последние годы наблюдается существенный рост производительности специализированных алгоритмов решения SAT и TQBF-задач, общее количество переменных в булевой модели свойства для современных решателей может измеряться тысячами.

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

Следует отметить, что метод булевых ограничений является достаточно общим методом качественного анализа ДДС на конечном интервале времени. Он применим не только к автономным системам, но и к системам с управляющими входами, к системам с глубиной памяти больше единицы, к ДДС общего вида, когда функция переходов неразрешима относительно состояния xt и имеет вид F(xt, xt-1) = 0. Для ДДС со входами особое значение приобретает свойство управляемости и его различные вариации. Кроме задач анализа ДДС метод булевых ограничений применим к задачам синтеза обратной связи (статической или динамической, по состоянию или по входу), обеспечивающих в синтезируемой системе выполнение требуемого динамического свойства.

Исследование выполнено при поддержке РФФИ, проект № 18-07-00596/18.

Библиографическая ссылка

Опарин Г.А., Богданова В.Г., Пашинин А.А. МЕТОД БУЛЕВЫХ ОГРАНИЧЕНИЙ В КАЧЕСТВЕННОМ АНАЛИЗЕ ДВОИЧНЫХ ДИНАМИЧЕСКИХ СИСТЕМ // Международный журнал прикладных и фундаментальных исследований. – 2018. – № 9. – С. 19-29;
URL: https://applied-research.ru/ru/article/view?id=12381 (дата обращения: 18.03.2020). Предлагаем вашему вниманию журналы, издающиеся в издательстве «Академия Естествознания»

Отправить свою хорошую работу в базу знаний просто. Используйте форму, расположенную ниже

Студенты, аспиранты, молодые ученые, использующие базу знаний в своей учебе и работе, будут вам очень благодарны.

Размещено на http://www.allbest.ru/

Задание

управление автоматический найквист частотный

Произвести анализ динамических свойств системы автоматического управления, заданной структурной схемой, представленной на рисунке 1, включающей следующие этапы:

Выбор и обоснование методов исследования, построение математической модели САУ;

Расчетная часть, включающая математическое моделирование САУ на ЭВМ;

Анализ устойчивости математической модели объекта управления и САУ;

Исследование устойчивости математической модели объекта управления и САУ.

Структурная схема исследуемой САУ, где, передаточные функции объекта управления (ОУ), исполнительного механизма (ИМ), датчика (Д) и корректирующего устройства (КУ)

Значения коэффициентов К1, К2, К3, К4, Т1, Т2, Т3 и Т4 приведены в таблице1.

Вариант задания на курсовую работу

	Параметры

Введение

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

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

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

1 . Анализ устойчивости САУ по критерию Найквиста

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

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

Система, устойчивая в разомкнутом состоянии, сохранит устойчивость и после её замыкания отрицательной обратной связью, если годограф КЧХ в разомкнутом состоянии W(jщ) не охватывает в комплексной плоскости точку с координатами (-1; j0).

В приведенной формулировке критерия Найквиста считается, что годограф КЧХ W(jщ) «не охватывает» точку (-1; j0), если равен нулю общий угол поворота вектора проведенный из указанной точки к годографу W(jщ) при изменении частоты от щ=0 до щ > ?.

Если годограф КЧХ W(jщ) при некоторой частоте называемой критической частотой щк, проходит через точку (-1; j0), то переходный процесс в замкнутой системе представляет собой незатухающие колебания с частотой щк, т.е. система оказывается на границе устойчивости выраженные следующим образом:

Здесь W(p) - передаточная функция разомкнутой САУ. Предположим, что разомкнутая система устойчива. Тогда для устойчивости замкнутой САУ необходимо и достаточно, чтобы годограф амплитудно-фазовой характеристики W(jw) разомкнутой системы (указанная характеристика получается из W(p) заменой p=jw) не охватывал точку с координатами (-1, j0). Частота, на которой |W(jw)| = 1, называется частотой среза (w ср).

Для оценки насколько далеко от границы устойчивости находится система, вводятся понятие запасов устойчивости. Запас устойчивости по амплитуде (модулю) указывает, во сколько раз необходимо изменить длину радиуса-вектора годографа АФХ, чтобы, не меняя фазового сдвига, вывести систему на границу устойчивости. Для абсолютно устойчивых систем запас устойчивости по модулю DК вычисляется по формуле:

где частота w 0 определяется из соотношения arg W(jw 0) = - 180 0 .

Запас устойчивости по амплитуде DК вычисляется и по формуле:

DК = 1 - К 180 ;

где К 180 - значение коэффициента передачи при фазовом сдвиге -180°.

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

Запас устойчивости по фазе Dj вычисляется по формуле:

Dj = 180° - j К=1 ;

где j К=1 - значение фазового сдвига при коэффициенте передачи К = 1;

Величина Dj = 180 0 + arg W (j; w ср) определяет запас устойчивости по фазе. Из критерия Найквиста следует, что устойчивая в разомкнутом состоянии САУ будет устойчивой и в замкнутом состоянии, если сдвиг по фазе на частоте среза не достигает - 180°. Выполнение этого условия можно проверить, построив логарифмические частотные характеристики разомкнутой САУ.

2. Исследование устойчивости САУ по критерию Найквиста

Исследование устойчивости по критерию Найквиста путем анализа АФЧХ при разомкнутой САУ. Для этого разрываем систему как показано на структурной схеме исследуемой САУ:

Структурная схема исследуемой САУ

Ниже представлены передаточные функции объекта управления (ОУ), исполнительного механизма (ИМ), датчика (Д) и корректирующего устройства (КУ):

Значения коэффициентов по заданию следующие:

К1 =1,0; К2 = 0,2; К3 = 2; К4 = 1,0; Т1 = 0,4; Т2 = 0,2; Т3 = 0,07; Т4 = 0,4.

Произведем расчет передаточной функции после разрыва системы:

W(р) = W ку (р) Ч W им (р) ЧW оу (р) ЧW д (р);

W(р) = Ч Ч Ч

Подставив заданные коэффициенты в функцию получим:

Анализируя данную функцию в программе математического моделирования («МАТLАВ»), получим годограф амплитудно-фазочастотной характеристики (АФЧХ) разомкнутой САУ на комплексной плоскости, приведенную на рисунке.

Годограф АФЧХ разомкнутой САУ на комплексной плоскости.

Исследование устойчивости САУ по АФЧХ

Вычисляем коэффициент передачи при фазовом сдвиге -180°, К 180 = 0,0395.

Запас устойчивости по амплитуде DК по формуле:

DК = 1 - К 180 = 1 - 0,0395= 0,9605; где К 180 = 0,0395.

Определим запас по фазе Dj:

запас устойчивости по фазе Dj определяется по формуле: Dj = 180° - j К=1 ; где j К=1 - значение фазового сдвига при коэффициенте передачи К = 1. Но так как, j К=1 в нашем случае не наблюдается, (амплитуда всегда меньше единицы), то исследуемая система устойчива при любом значении фазового сдвига (САУ устойчива на всем диапазоне частот).

Исследование устойчивости САУ по логарифмическим характеристикам

Логарифмическая амплитудно-частотная характеристика разомкнутой САУ

Логарифмическая фазочастотная характеристика разомкнутой САУ

Используя программу математического моделирования («МАТLАВ»), получим логарифмические характеристики исследуемой САУ, которые представлены на рисунке 4 (логарифмическая амплитудно-частотная характеристика) и рисунке 5 (логарифмическая фазочастотная характеристика), где;

L(w) = 20lg|W (j; w) |).

Логарифмический критерий устойчивости САУ представляет собой выражение критерия Найквиста в логарифмической форме.

Для нахождения из значения фазового сдвига 180° (рисунок 5) проводим горизонтальную линию до пересечения с ЛФЧХ, с этой точки пересечения проводим вертикальную линию до пересечения с ЛФЧХ (рисунок 4). Получаем значение коэффициента передачи при фазовом сдвиге 180°:

20lgК 180 ° = - 28,05862;

при этом К 180 ° = 0,0395 (DК" = 28,05862).

Запас устойчивости по амплитуде находится продолжением вертикальной линии до значения 20lgК 180 ° = 0.

Для нахождения запаса устойчивости по фазе, пропускается горизонтальная линия по линии 20lgК 180 ° = 0 до пересечения с ЛАЧХ и пропускается из этой точки вертикальная линия до пересечения с ЛФЧХ. При этом разница между найденным значением фазового сдвига и фазовым сдвигом равным 180° и будет запасом устойчивости по фазе.

Dj = 180° - j К;

Dj = 180° - 0 = 180°.

где: j К - найденное значение фазового сдвига;

Так как ЛФЧХ исследуемой САУ лежит ниже линии 20lgК 180 ° = 0, поэтому САУ будет иметь запас устойчивости по фазе при любом значении фазового сдвига от нуля до 180°.

Вывод: проанализировав ЛАЧХ и ЛФЧХ, следует что исследуемая САУ устойчива на всем диапазоне частот.

Заключение

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

Библиография

1. И.Ф. Бородин, Ю.А. Судник. Автоматизация технологических процессов. Учебник для вузов. Москва. «Колос», 2004.

2. В.С. Гутников. Интегральная электроника в измерительных устройствах. «Энергоатомиздат». Ленинградское отделение, 1988.

3. Н.Н. Иващенко. Автоматическое регулирование. Теория и элементы систем. Москва. «Машиностроение», 1978.

Размещено на Allbest.ru

...

Подобные документы

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

курсовая работа , добавлен 21.02.2016


Исследование системы управления частотой вращения двигателя с корректирующей цепью и без нее. Оценка устойчивости системы по критериям Гурвица, Михайлова и Найквиста. Построение логарифмических амплитудно-частотной и фазово-частотной характеристик.

курсовая работа , добавлен 22.03.2015


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

курсовая работа , добавлен 24.03.2013


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

курсовая работа , добавлен 13.09.2010


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

курсовая работа , добавлен 12.08.2014


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

лабораторная работа , добавлен 15.12.2010


Построение элементарной схемы и исследование принципа работы системы автоматического управления, ее значение в реализации способа поднастройки системы СПИД. Основные элементы системы и их взаимосвязь. Анализ устойчивости контура и его оптимальных частот.

контрольная работа , добавлен 12.09.2009


Определение передаточной функции разомкнутой системы, стандартной формы ее записи и степени астатизма. Исследование амплитудно-фазовой, вещественной и мнимой частотных характеристик. Построение годографа АФЧХ. Алгебраические критерии Рауса и Гурвица.

курсовая работа , добавлен 09.05.2011


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

дипломная работа , добавлен 19.01.2017


Функциональная схема системы автоматического регулирования температуры приточного воздуха в картофелехранилище. Определение закона регулирования системы. Анализ устойчивости по критериям Гурвица и Найквиста. Качество управления по переходным функциям.

Введение 4

Априорный анализ динамических систем 5

Прохождение случайного сигнала через линейную систему 5

Эволюция фазового вектора системы 7

Эволюция ковариационной матрицы фазового вектора системы 8

Статистическая линеаризация 8

Первый способ 9

Второй способ 10

Вычисление коэффициентов линеаризации 10

Неоднозначность в нелинейных звеньях 14

Нелинейное звено, охваченное обратной связью 15

Моделирование случайных процессов 16

Формирующий фильтр 16

Моделирование белого шума 17

Оценивание статистических характеристик динамических систем методом Монте-Карло 18

Точность оценок 18

Нестационарные динамические системы 20

Стационарные динамические системы 21

Апостериорный анализ динамических систем 22

Фильтр Калмана 22

Модель движения 22

Модель измерений 23

Коррекция 23

Прогноз 23

Оценивание 23

Использование калмановской фильтрации в нелинейных задачах 25

Метод наименьших квадратов 27

Построение оценок 27

Прогноз 29

Использование метода наименьших квадратов в нелинейных задачах 29

Построение матрицы Коши 30

Моделирование измерений 30

Численные методы 31

Специальные функции 31

Моделирование случайных величин 31

Равномерно распределенные случайные величины 31

Гауссовские случайные величины 32

Случайные векторы 33

Интеграл вероятностей 34

Полиномы Чебышева 36

Интегрирование обыкновенных дифференциальных уравнений 36

Методы Рунге-Кутты 36

Точность результатов численного интегрирования 37

Вложенный метод Дормана-Принса 5(4) порядка 37

Многошаговые методы 39

Методы Адамса 39

Интегрирование уравнений с запаздывающим аргументом 40

Сравнение вычислительных качеств методов 40

Задача Аренсторфа 40

Эллиптические функции Якоби 41

Задача двух тел 41

Уравнение Ван-дер-Поля 42

«Брюсселятор» 42

Уравнение Лагранжа для висячей струны 42

«Плеяды» 42

Оформление пояснительной записки 43

Титульный лист 43

Раздел «Введение» 44

Раздел «Теория» 44

Раздел «Алгоритм» 44

Раздел «Программа» 45

Раздел «Результаты» 45

Раздел «Выводы» 45

Раздел «Список использованных источников» 45

Приложения 45

Литература 47

Введение

В настоящем учебном пособии содержатся методические указания к выполнению заданий курсовых проектов и к проведению практических занятий по курсу «Основы статистической динамики».

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

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

Статистическая линеаризация

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

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

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

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

Величина выбирается исходя из условия равенства математических ожиданий нелинейного и линеаризованного сигналов и носит название статистической средней характеристики эквивалентного звена:

,

где – плотность распределения входного сигнала .

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

– математическое ожидание входного сигнала ;
– статистический коэффициент усиления эквивалентного звена по средней составляющей .

Т.о. эквивалентная зависимость в данном случае приобретает вид:

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

Первый способ

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

,

где – дисперсия входного случайного воздействия.

Знак в выражении для определяется характером зависимости в окрестности значения аргумента . Если возрастает, то , а если убывает, то .

Второй способ

Значение по второму способу выбирается из условия минимизации средней квадратической ошибки линеаризации:

Окончательное соотношение для вычисления коэффициента по второму способу имеет вид:

.

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

Формирующий фильтр

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

при одинаковых степенях .

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

Например, спектральная плотность процесса , подлежащего моделированию имеет вид:

,

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

Очевидно, что числитель и знаменатель искомой передаточной функций должны иметь порядки 1 и 2 (в самом деле, будучи возведенной в квадрат по модулю передаточная функция образует частное полиномов 2-й и 4-й степеней)

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

,

а квадрат ее модуля:

Приравняем полученные соотношения:

Вынесем за скобку и в правой части равенства, приравнивая тем самым коэффициенты при нулевых степенях :

,

откуда с очевидностью вытекают следующие равенства:

; ; ; .

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

Моделирование белого шума

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

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

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

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

где – эмпирически определяемая граница полосы пропускания динамической системы.

Точность оценок

Оценки математического ожидания

и дисперсии

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

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

,

где
– истинное значение математического ожидания случайной величины ,
– среднеквадратическое отклонение случайной величины ,
– интеграл вероятностей.

На основе приведенного выше соотношения величина может быть определена следующим образом:

,

где – функция, обратная по отношению к интегралу вероятностей .

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

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

.

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

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

с точностью до величины , которая за неимением более точной информации может быть приблизительно определена из соотношения:

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

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

Например, для гауссовского закона распределения случайная величина

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

распределена по закону также с степенью свободы.

Фильтр Калмана

Модель движения

Как известно, Фильтр Калмана предназначен для оценивания вектора состояния линейной динамической системы, модель эволюции которого может быть записана в виде:

где
– матрица Коши, определяющая изменение вектора состояния системы в ее собственном движении (без управляющих и шумовых воздействий) от момента времени к моменту времени ;
– вектор вынуждающих неслучайных воздействий на систему (например, управляющих воздействий) в момент времени ;
– матрица влияния вынуждающих воздействий в момент времени на вектор состояния системы в момент времени ;
– вектор случайных независимых центрированных воздействий на систему в момент времени ;
– матрица влияния случайных воздействий в момент времени на вектор состояния системы в момент времени .

Модель измерений

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

где – матрица, связывающая векторы состояния и измерений в один и тот же момент времени .

Коррекция

Основу Фильтра Калмана составляют соотношения коррекции, являющиеся результатом минимизации следа ковариационной матрицы апостериорной плотности распределения линейной (по вектору измерений ) оценки вектора состояния системы :

Прогноз

Дополняя соотношения коррекции соотношениями прогноза, основанными на линейных свойствах модели эволюции системы:

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

Оценивание

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

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

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

Поясним алгоритм калмановской фильтрации при помощи рисунка.

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

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

В осях координат , (в канале измерений) в моменты времени и изображены результаты измерений и :

,

где
– истинное значение вектора измерений в момент времени ;
– вектор ошибок измерений, реализовавшихся в момент времени .

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

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

Метод наименьших квадратов

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

Построение оценок

Для случая линейной модели равноточных измерений:

имеем следующий алгоритм оценивания фазового вектора:

.

Для случая неравноточных измерений вводится в рассмотрение матрица , содержащая на диагонали весовые коэффициенты. С учетом весовых коэффициентов предыдущее соотношение примет вид:

.

Если в качестве весовой использовать матрицу, обратную к ковариационной матрице ошибок измерений , то с учетом того обстоятельства, что получим:

.

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

На рисунке показано некоторое возможное взаимное расположение моментов времени, к которым отнесены измерения и момента времени, к которому отнесен вектор оцениваемых параметров.

Для каждого вектора справедливо следующее соотношение:

, при .

Таким образом, в результирующем соотношении метода наименьших квадратов вектор и матрица имеют следующую структуру:

; .

где
– определяет неслучайное вынуждающее воздействие на систему;
– определяет случайное воздействие на систему.

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

где – ковариационная матрица вектора .

Построение матрицы Коши

В задачах построения оценок методами статистической обработки измерений часто встречается задача построения матрицы Коши. Эта матрица связывает фазовые векторы системы, отнесенные к разным моментам времени, в собственном движении.

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

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

; .

Моделирование измерений

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

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

Численные методы

Специальные функции

Случайные векторы

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

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

При сумма асимптотического ряда становится практически равной 1.

КИНЕТИКА БИОЛОГИЧЕСКИХ ПРОЦЕССОВ

Как можно описать динамику биологических систем? В каждый момент времени биологическая система обладает набором некоторых характеристик. Например, наблюдая за популяцией какого-то вида, можно регистрировать ее численность, площадь занимаемой территории, количество доступного питания, температуру окружающей среды и т. д. Протекание химической реакции можно характеризовать концентрациями участвующих веществ, давлением, температурой, уровнем кислотности среды. Совокупность значений всех характеристик, которые исследователь выбрал для описания системы, является состоянием системы в каждый момент времени. При создании модели в указанной совокупности выделяют переменные и параметры. Переменные – это те величины, изменения которых в первую очередь интересует исследователя, параметры – условия «внешней среды». Именно для выбранных переменных составляют уравнения, отражающие закономерности изменения системы во времени. Например, при создании модели роста культуры микроорганизмов, в качестве переменной обычно выступает ее численность, а в качестве параметра – скорость размножения. Возможно, существенной окажется температура, при которой происходит рост, тогда этот показатель также включается в модель в качестве параметра. А если, например, уровень аэрации всегда является достаточным и не оказывает никакого влияния на ростовые процессы, тогда его вообще не включают в модель. Как правило, параметры остаются неизменными во время эксперимента, однако стоит отметить, что это не всегда так.

Описывать динамику биологической системы (то есть изменение ее состояния во времени) можно как дискретными, так и непрерывными моделями. В дискретных моделях предполагается, что время представляет собой дискретную величину. Это соответствует регистрации значений переменных через определенные фиксированные интервалы времени (например, раз в час или раз в год). В непрерывных моделях биологическая переменная является непрерывной функцией времени, обозначаемой, например, x (t ).

Часто большое значение имеют начальные условия модели – состояние исследуемой характеристики в начальный момент времени, т.е. при t = 0.

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

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

Если правая часть дифференциального уравнения вида явно не зависит от времени, т.е. справедливо:

то такое уравнение называется автономным (система, описываемая таким уравнением, называется автономной ). Состояние автономных систем в каждый момент времени характеризуется одной единственной величиной - значением переменной x в данный момент времени t .

Зададимся вопросом: пусть дано дифференциальное уравнение для x (t ), можно ли найти все функции x (t ), удовлетворяющие этому уравнению? Или: если известно начальное значение некоторой переменной (например, начальный размер популяции, концентрация вещества, электропроводность среды и т.п.) и имеется информация о характере изменения этой переменной, то можно ли предсказать, каким будет ее значение во все последующие моменты времени? Ответ на поставленный вопрос звучит следующим образом: если заданы начальные условия при и для уравнения выполнены условия теоремы Коши (функция , заданная в некоторой области, и ее частная производная непрерывны в этой области), то имеется единственное решение уравнения, удовлетворяющее заданным начальным условиям. (Напомним, что любая непрерывная функция, которая удовлетворяет дифференциальному уравнению, называется решением этого уравнения.) Это означает, что мы можем однозначно предсказывать поведение биологической системы, если известны характеристики ее начального состояния, и уравнение модели удовлетворяет условиям теоремы Коши.

Стационарное состояние. Устойчивость

Будем рассматривать автономное дифференциальное уравнение

В стационарном состоянии значения переменных в системе не меняются со временем, то есть скорость изменения значений переменных равна 0: . Если левая часть уравнения (1.2) равна нулю, то и правая тоже равна нулю: . Корни этого алгебраического уравнения являются стационарными состояниями дифференциального уравнения (1.2).

Пример1.1: Найдите стационарные состояния уравнения .

Решение : Перенесем слагаемое, не содержащее производную, в правую часть равенства: . По определению в стационарном состоянии выполняется равенство: . Значит, должно выполняться равенство . Решаем уравнение:

,

Итак, уравнение имеет 3 стационарных состояния: , .

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

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

Состояние равновесия уравнения устойчиво по Ляпунову, если для любого всегда можно найти такое , что если , то для всех .

Существует аналитический метод исследования устойчивости стационарного состояния – метод Ляпунова. Для его обоснования напомним формулу Тейлора .

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

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

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

Пример1.2: Разложите функцию в ряд Тейлора в окрестности точки до 4 порядка.

Решение: Запишем ряд Тейлора до 4-го порядка в общем виде:

Найдем производные заданной функции в точке :

,

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

Аналитический метод исследования устойчивости стационарного состояния (метод Ляпунова ) состоит в следующем. Пусть – стационарное состояние уравнения . Зададим небольшое отклонение переменной x от ее стационарного значения: , где . Подставим выражение для точки x в исходное уравнение: . Левая часть уравнения примет вид: , поскольку в стационарном состоянии скорость изменения значения переменой равна нулю: . Правую часть разложим в ряд Тейлора в окрестности стационарного состояния, учитывая, что , оставим только линейный член в правой части уравнения:

Получили линеаризованное уравнение или уравнение первого приближения . Величина есть некоторая постоянная величина, обозначим ее a : . Общее решение линеаризованного уравнения имеет вид: . Это выражение описывает закон, по которому будет изменяться во времени заданное нами отклонение от стационарного состояния . Отклонение будет со временем затухать, т.е. при , если показатель степени в экспоненте будет отрицательным, т.е. . По определению стационарное состояние будет устойчивым . Если же , то с увеличением времени отклонение будет только увеличиваться, стационарное состояние – неустойчивое . В случае, когда уравнение первого приближения ответа на вопрос об устойчивости стационарного состояния дать не может. Необходимо рассматривать члены более высокого порядка в разложении в ряд Тейлора.

Кроме аналитического метода исследования устойчивости стационарного состояния существует и графический.

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

Решение: Найдем особые точки:

,

,

Строим график функции (рис.1.1).

Рис. 1.1. График функции (пример 1.3).

Определим по графику устойчиво ли каждое из найденных стационарных состояний. Зададим небольшое отклонение изображающей точки от особой точки влево: . В точке с координатой функция принимает положительное значение: или . Последнее неравенство означает, что со временем координата должна увеличиваться, то есть изображающая точка должна возвратиться к точке . Теперь зададим небольшое отклонение изображающей точки от особой точки вправо: . В этой области функция сохраняет положительное значение, следовательно, со временем координата x также увеличивается, то есть изображающая точка будет отдаляться от точки . Таким образом, малое отклонение выводят систему из стационарного состояния, следовательно, по определению особая точка неустойчива. Аналогичные рассуждения приводят к тому, что любое отклонение от особой точки со временем затухает, стационарное состояние устойчиво. Отклонение изображающей точки в любом направлении от стационарного состояния приводит к ее удалению от точки , это неустойчивое стационарное состояние.

Решение системы линейных дифференциальных уравнений

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

Анализ системы уравнений начинается с нахождения стационарных состояний. У систем вида (1.3) особая точка единственна, ее координаты – (0,0). Исключение составляет вырожденный случай, когда уравнения можно представить в виде:

(1.3*)

В этом случае все пары , удовлетворяющие соотношению , являются стационарными точками системы (1.3*). В частности, точка (0,0) также является стационарной для системы (1.3*). На фазовой плоскости в данном случае имеем прямую с коэффициентом наклона , проходящую через начало координат, каждая точка которой является особой точкой системы (1.3*) (см. таблицу1.1, пункт 6).

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

Общее решение системы двух линейных уравнений имеет вид:

Характеристические числа могут быть выражены через коэффициенты линейных уравнений следующим образом:

Характеристические числа могут быть 1) действительными разных знаков, 2) действительными одного знака, 3) комплексно сопряженными, а также, в вырожденных случаях, 4) чисто мнимыми, 5) действительными совпадающими, 6) действительными, один из которых (или оба) равен нулю. Эти случаи определяют тип поведения решения системы обыкновенных дифференциальных уравнений. Соответствующие фазовые портреты представлены в таблице1.1.

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

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

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

Рассмотрим систему линейных дифференциальных уравнений:

Построение фазового портрета начинаем с построения главных изоклин (изоклина – линия, на всём протяжении которой наклон фазовой кривой (траектории), определяемый уравнением, сохраняет постоянное значение). Для системы двух линейных дифференциальных уравнений – это всегда прямые, проходящие через начало координат. Уравнение изоклины горизонтальных касательных : . Уравнение изоклины вертикальных касательных : . Для дальнейшего построения фазового портрета полезно построить изоклину касательных, проходящих под углом . Для нахождения соответствующего уравнения изоклины необходимо решить уравнение . Можно находить и изоклины касательных других углов, пользуясь приблизительными значениями тангенсов углов. В построении фазового портрета также может помочь ответ на вопрос, под каким углом фазовые траектории должны пересекать координатные оси. Для этого в уравнение изоклины подставляем соответствующие равенства (для определения угла пересечения с оcью OY) и (для определения угла пересечения с оcью OХ).

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

Постройте фазовый и кинетический портрет системы.

Решение: Координаты особой точки – (0,0). Коэффициенты линейных уравнений равны: , , , . Определим тип стационарного состояния (см. раздел о характеристических числах):

Таким образом, характеристические корни являются мнимыми: , следовательно, особая точка рассматриваемой линейной системы имеет тип центр (рис. 1.2а).

Уравнение изоклины горизонтальных касательных: , уравнение изоклины вертикальных касательных: . Под углом в 45 ° траектории системы пересекают прямую .

После построения фазового портрета необходимо определить направление движения по найденным траекториям. Это можно сделать следующим образом. Возьмем произвольную точку на любой траектории. Например, на изоклине горизонтальных касательных (1,1). Подставим координаты этой точки в систему уравнений. Получим выражения для скоростей изменения переменных x , y в этой точке:

Полученные значения показывают, что скорость изменения переменной x – отрицательная, то есть ее значение должно уменьшаться, а переменная y не изменяется. Отмечаем полученное направление стрелкой. Таким образом, в рассматриваемом примере движение по фазовым траекториям направлено против часовой стрелки. Подставляя в систему координаты разных точек можно получить «карту» направлений скоростей, так называемое векторное поле .

Рис 1.2. Фазовый (а) и кинетический (б) портрет системы, пример 1.4

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

Построить кинетический портрет системы означает построить графики зависимости величин переменных x , y от времени. По фазовому портрету можно построить кинетический и наоборот. Одной фазовой траектории соответствует одна пара кинетических кривых. Выберем на фазовом портрете произвольную точку на произвольной фазовой траектории. Это начальная точка, соответствующая моменту времени . В зависимости от направления движения в рассматриваемой системе значения переменных x , y либо уменьшаются, либо увеличиваются. Пусть координаты начальной точки – (1,1). Согласно построенному фазовому портрету, стартуя из этой точки, мы должны двигаться против часовой стрелки, координаты x и y при этом будут уменьшаться. С течением времени координата x проходит через 0, значение y при этом остается положительным. Далее координаты x и y продолжают уменьшаться, координата y проходит через 0 (значение x при этом отрицательно). Величина x достигает минимального значения на изоклине вертикальных касательных, затем начинает увеличиваться. Величина y своего минимального значения достигает на изоклине горизонтальных касательных (значение x в этот момент времени отрицательно). Далее и величина x , и величина y увеличиваются, возвращаясь к начальным значениям (рис. 1.2б).

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

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

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

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

Вне

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