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

Як можна описати динаміку біологічних систем? У кожний момент часу біологічна система має набір деяких характеристик. Наприклад, спостерігаючи за популяцією якогось виду, можна реєструвати її чисельність, площу займаної території, кількість доступного харчування, температуру довкілляі т. д. Протікання хімічної реакціїможна характеризувати концентраціями речовин, що беруть участь, тиском, температурою, рівнем кислотності середовища. Сукупність значень всіх параметрів, які дослідник вибрав для опису системи, є станом системи в кожний момент часу. При створенні моделі у зазначеній сукупності виділяють змінні та параметри. Змінні – це величини, зміни яких насамперед цікавить дослідника, параметри – умови «довкілля». Саме обраних змінних становлять рівняння, що відбивають закономірності зміни системи у часі. Наприклад, при створенні моделі зростання культури мікроорганізмів, як змінна зазвичай виступає її чисельність, а як параметр – швидкість розмноження. Можливо, суттєвою виявиться температура, при якій відбувається зростання, тоді цей показник також включається до моделі як параметр. А якщо, наприклад, рівень аерації завжди є достатнім і не впливає на ростові процеси, тоді його взагалі не включають у модель. Як правило, параметри залишаються незмінними під час експерименту, проте слід зазначити, що це не завжди так.

Описувати динаміку біологічної системи (тобто зміна її стану у часі) можна дискретними, і безперервними моделями. У дискретних моделях передбачається, що час є дискретну величину. Це відповідає реєстрації значень змінних через певні фіксовані інтервали часу (наприклад, раз на годину або на рік). У безперервних моделях біологічна змінна є безперервною функцією часу, що позначається, наприклад, 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) дійсними одного знака; дорівнює нулю. Ці випадки визначають тип поведінки рішення системи звичайних диференціальних рівнянь. Відповідні фазові портрети представлені таблиці1.1.

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

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

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

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

Побудова фазового портрета починаємо з побудови головних ізоклін(Ізокліна - лінія, на всьому протязі якої нахил фазової кривої (траєкторії), що визначається рівнянням, зберігає постійне значення). Для системи двох лінійних диференціальних рівнянь – це прямі, які проходять через початок координат. Рівняння ізокліни горизонтальних дотичних: . Рівняння ізокліни вертикальних дотичних: . Для подальшої побудови фазового портрета корисно побудувати ізокліну дотичних, що проходять під кутом. Для знаходження відповідного рівняння ізокліни необхідно вирішити рівняння . Можна знаходити і ізокліни щодо інших кутів, користуючись приблизними значеннями тангенсів кутів. У побудові фазового портрета може допомогти відповідь питанням, під яким кутом фазові траєкторії повинні перетинати координатні осі. Для цього в рівняння ізокліни підставляємо відповідні рівності (для визначення кута перетину з віссю OY) та (для визначення кута перетину з віссю 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б).

Дослідження стійкості стаціонарних станів нелінійних систем другого порядку

Нехай біологічна система описується системою двох автономних диференціальних рівнянь другого порядку загального вигляду:

Стаціонарні значення змінних системи визначаються з рівнянь алгебри:

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

Поза

Маємо, , , особлива точка груба. Характеристичні коріння системи першого наближення рівні , обидва дійсні і негативні, отже, в околиці нульової особливої ​​точки поведінка фазових траєкторій системи відповідатиме типу стійкий вузол.

Вступ

Оскільки концепція нелінійної динамічної системи досить багата, щоб охопити надзвичайно широке коло процесів, у яких майбутня поведінка системи визначається минулим, методи аналізу, розроблені в цій галузі, корисні у величезній різноманітності контекстів

Нелінійна динаміка входить у літературу принаймні трьома способами. По-перше, бувають випадки, коли експериментальні дані про зміну в часі однієї або декількох величин збираються і аналізуються з використанням методик, заснованих на нелінійній динамічній теорії, з мінімальними припущеннями щодо рівнянь, що лежать в основі, що керують процесом, який і створює ці дані. Тобто це випадок, при якому прагнуть знайти кореляції в даних, які можуть направити розробку математичної моделі, замість того, щоб спочатку вгадати модель, а потім порівняти її з даними.

По-друге, є випадки, коли нелінійна динамічна теорія може бути використана для твердження, що деяка спрощена модель має демонструвати важливі особливостіданої системи, з чого випливає, що описову модель можна побудувати та вивчити у широкому діапазоні параметрів. Часто це призводить до моделей, які поводяться якісно по-різному при різних параметрах і демонструють, що одна область показує поведінку, дуже схожу на поведінку, що спостерігається в реальній системі. У багатьох випадках поведінка моделі досить чутлива до змін параметрів, тому якщо параметри моделі можуть бути виміряні в реальній системі, модель демонструє реалістичну поведінку при цих значеннях, і можна бути впевненим, що модель охопила суттєві особливості система.

По-третє, трапляються випадки, коли модельні рівняння будуються на основі докладних описіввідомої фізики. Потім чисельні експерименти можуть дати інформацію про змінні, недоступні фізичним експериментам.

Спираючись на другий шлях, ця робота є розширенням моєї попередньої роботи «Нелінійна динамічна модель взаємозалежних галузей виробництва», а також іншої роботи (Dmitriev, 2015)

Усі необхідні визначення та інші теоретичні відомості, необхідні у роботі, з'являтимуться у першому розділі, у міру їхньої необхідності. Тут же буде наведено два визначення, що необхідні для розкриття самої теми дослідження.

Для початку дамо визначення системної динаміки. Відповідно до одного з визначень, системна динаміка - підхід імітаційного моделювання, який завдяки своїм методам та інструментам допомагає оцінити структуру складних системта їх динаміку (Штерман). Варто додати, що системна динаміка також і метод моделювання, який використовують з метою відтворення вірних (з точки зору точності) комп'ютерних моделей для складних систем заради їх майбутнього використання для того, щоб створити ефективнішу компанію/організацію, а також покращити методи взаємодії з даною системою. Переважно, необхідність у системній динаміці виникає при зіткненні з довгостроковими, стратегічними моделями, а також варто відзначити, що вона досить абстрактна.

Говорячи про нелінійну диференціальну динаміку, ми розглядатимемо нелінійну систему, яка за визначенням, є системою, в якій зміна результату не пропорційна зміні вхідних параметрів, і в якій функція описує залежність зміни в часі та положенні точки у просторі (Boeing, 2016).

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

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

Для досягнення поставленої мети було виділено такі завдання:

Визначення сталості системи.

Побудова фазових портретів.

Знаходження інтегральних траєкторій систем.

Побудова імітаційних моделей.

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

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

Що стосується причини вибору якісного, а не кількісного аналізу, то варто зазначити, що в переважній кількості випадків результати та висновки з якісного аналізу динамічних систем виявляються значнішими за результати їх кількісного аналізу. У такій ситуації доречно зазначити висловлювання В.П. Мілованова, у якому він стверджує, що зазвичай вважають, що результати, очікувані при застосуванні математичних методів для аналізу реальних об'єктів, повинні зводитися до чисельного результату. У цьому сенсі якісні методи покладається дещо інше завдання. У ньому акцентується увагу на досягненні результату, що описує якості системи, пошуку характерних особливостей всіх явищ в цілому, на прогнозування. Зрозуміло, важливо розуміти, як зміниться попит при зміні цін на певний вид товарів, але не варто забувати, що набагато важливіше розуміння, чи в таких умовах настане дефіцит чи надлишок цих товарів (Дмитрієв, 2016).

Об'єктом цього дослідження є нелінійна диференціальна та системна динаміка.

У такому разі предмет дослідження - опис процесу взаємодії між компаніями через нелінійну диференціальну та системну динаміку.

Говорячи про практичне застосування дослідження, варто відразу розбити його на дві частини. А саме на теоретичну, тобто якісний аналіз систем, та практичну, в якій буде розглянуто побудову імітаційних моделей.

Теоретична частина даного дослідження дає базові поняття та явища. У ній розглянуті прості диференціальні системи, як і в роботах багатьох інших авторів (Teschl, 2012; Nolte, 2015), але дозволяють описати взаємодію між компаніями. Грунтуючись на цьому надалі можна буде проводити більш поглиблені дослідження, або ж починати своє знайомство з тим, що являє собою якісний аналіз систем.

Практичну частину роботи можна використовуватиме створення системи підтримки прийняття рішень. Система підтримки прийняття рішень – автоматизована інформаційна система, Націлена на підтримку бізнесу або прийняття рішень в організації, що дозволяє вибрати між безліччю різних альтернатив (Keen, 1980). Нехай в даний момент моделі і не мають високої точності, але змінивши їх під конкретну компанію, можна досягти більш вірних результатів. Таким чином, при зміні в них різних параметрів і умов, здатних виникнути на ринку, можна отримати прогноз на майбутнє і прийняти заздалегідь вигідніше рішення.

1. Взаємодія підприємств за умов мутуалізму

У роботі будуть представлені двовимірні системи, які досить прості порівняно з системами вищого порядку, але при цьому дозволяють продемонструвати необхідні нам взаємини між організаціями.

Почати роботу варто з вибору виду взаємодії, яка надалі і буде описуватися, оскільки для кожного з видів системи, що їх описують, нехай небагато, але різні. На малюнку 1.1, зазначена модифікована під економічну взаємодію класифікація Юджима Одума для взаємодії популяцій (Одум, 1968), ґрунтуючись на якій надалі ми розглядатимемо взаємодію компаній.

Малюнок 1.1. Типи взаємодії між підприємствами

Грунтуючись на малюнку 1.1, виділимо 4 типи взаємодії і наведемо для кожного з них, що описує їхню систему рівнянь, засновану на моделі Мальтуса (Malthus, 1798). Відповідно до неї, швидкість зростання має пропорційну залежність від поточної чисельності виду, інакше кажучи, її можна описати наступним диференціальним рівнянням:

де a - певний параметр, що залежить від природного приросту популяції. Також варто додати, що в системах, що розглядаються далі, всі параметри, а також змінні приймають невід'ємні значення.

Виробництво сировини - виробництво продукції, що аналогічно моделі хижака-жертва. Модель хижак-жертва, також відома як модель Лотки-Вольтерри, - пара нелінійних диференціальних рівнянь першого порядку, що описують динаміку біологічної системи з двома видами, один із яких хижаки, а інший жертви (Llibre, 2007). Зміна чисельності цих видів описується наступною системою рівнянь:

(1.2)

де - характеризує зростання продукції першого підприємства без впливу другого (у разі моделі хижак-жертва, зростання популяції жертв без хижаків),

Характеризує зростання продукції другого підприємства без впливу першого (зростання популяції хижаків без жертв),

Характеризує зростання продукції першого підприємства з урахуванням впливом нього другого (зростання чисельності жертв при взаємодії з хижаками),

Характеризує зростання продукції другого підприємства, враховуючи вплив нього першого (зростання чисельності хижаків за її взаємодії з жертвами).

На одного, на хижака, як видно із системи, а також класифікації Одума, їхня взаємодія накладає сприятливий вплив. На іншого несприятливе. Якщо розглядати в економічних реаліях, то як видно на малюнку найпростішим аналогом є фірма-виробник та її постачальник ресурсів, які відповідають хижакові та жертві відповідно. Таким чином, без сировини випуск продукції експоненційно знижується.

Конкуренція - це суперництво між двома і більше (у разі ми розглядаємо двовимірні системи, тому беремо саме двовидову конкуренцію) видами, економічними групами за території, обмежені ресурси чи інші цінності (Elton, 1968). Зміни в кількості видів, або кількості продукції в нашому випадку, описується системою нижче:

(1.3)

В даному випадку види або компанії, що випускають один продукт, несприятливо впливають один на одного. Тобто, за відсутності конкурента, зростання продукції експоненційно зросте.

Тепер перейдемо до симбіотичного взаємодії, у якому обидва підприємства мають одне одного позитивний вплив. Спочатку розглянемо мутуалізм. Мутуалізм - тип відносин між різними видами, при якому кожен з них отримує вигоду від дії іншого, причому варто зазначити, що присутність партнера є обов'язковою умовою існування (Thompson, 2005). Такий тип відносин описується системою:

(1.4)

Оскільки взаємодія між компаніями необхідна їх існування, то відсутність продукту однієї компанії, випуск товарів інший експоненційно знижується. Таке можливо, коли компанії просто не мають інших альтернатив для закупівель.

Розглянемо ще один вид симбіотичної взаємодії, протокооперацію. Протокооперація схожа на мутуалізм з єдиним винятком, немає потреби в обов'язковому існуванні партнера, оскільки, наприклад, існують інші альтернативи. Оскільки вони схожі, то їх системи виглядають практично аналогічно один одному:

(1.5)

Таким чином, відсутність продукту однієї компанії не заважає зростанню продукту іншої.

Звісно, ​​крім перелічених у пунктах 3 і 4, можна назвати й інші види симбіотичних відносин: комменсалізм і аменсалізм (Hanski, 1999). Але вони не згадуватимуться надалі, оскільки в коменсалізмі одному з партнерів байдуже його взаємодія з іншим, а ми таки розглядаємо випадки, коли вплив є. А аменсалізм не розглядається, оскільки з економічної точки зору таких відносин, коли одному їхня взаємодія шкодить, а іншому байдуже, просто не може бути.

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

Ця глава присвячена взаємодії підприємств за умов мутуалізму. У ній будуть розглянуті дві системи, що є подальшим розвитком систем, заснованих на моделі Мальтуса, а саме системами з обмеженнями на збільшення продукції, що накладаються.

Динаміку пари, пов'язаної мутуалістичними відносинами, як було зазначено вище, у першому наближенні можна описати системою:

(1.6)

Можна зауважити, що при великій початковій кількості продукції система необмежено зростає, а при малому виробництво продукції падає. У цьому полягає некоректність білінійного опису ефекту, що виникає при мутуалізмі. Щоб спробувати виправити картину, введемо фактор, що нагадує насичення хижака, тобто фактор, який дозволить зменшити швидкість зростання продукції при його надлишку. У цьому випадку приходимо до наступної системи:

(1.7)

де - зростання виробництва продукту першої компанії при її взаємодії з другою з урахуванням насичення,

Зростання виробництва продукту другої компанії при її взаємодії з першою з урахуванням насичення,

Коефіцієнти насичення.

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

1.1 Стійкість систем у першому наближенні

Стійкість систем у першому наближенні розглядається у багатьох, як іноземних (Hairer, 1993; Bhatia, 2002; Khalil, 2001; Strogatz, 2001 та інші), так і російськомовних роботах (Ахромєєва, 1992; Беллман, 1954;, Демид; 1959 та інші), та її визначення є базовим кроком для аналізу процесів, що відбуваються в системі. Для цього виконаємо такі необхідні кроки:

Знайдемо рівноважні точки.

Знайдемо матрицю Якобі системи.

Знайдемо власні значення матриці Якобі.

Класифікуємо рівноважні точки з теореми Ляпунова.

Розглянувши кроки, варто докладніше зупинитися на їхнє роз'яснення, тому дам визначення та опишу методи, якими ми користуватимемося в кожному з цих кроків.

Перший крок – пошук рівноважних точок. Для їхнього знаходження прирівняємо кожну функцію до нуля. Тобто вирішимо систему:

де a і b маються на увазі всі параметри рівняння.

Наступний крок – пошук матриці Якобі. У нашому випадку це буде матриця 2 на 2 з першими похідними в деякій точці , як представлено нижче:

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

Де точка відповідає рівноважним точкам, знайденим у першому кроці.

Знайшовши і перейдемо до четвертого кроку і скористаємося наступними теоремами Ляпунова (Parks, 1992):

Теорема 1: Якщо всі корені характеристичного рівняння мають негативну дійсну частину, то рівноважна точка відповідна початковій та лінеаризованій системам – асимптотично стійка.

Теорема 2: Якщо хоча б один з коренів характеристичного рівняння має позитивну дійсну частину, то рівноважна точка, що відповідає початковій та лінеаризованій системам, - асимптотично нестійка.

Також, дивлячись на і можна і точніше визначити тип стійкості, ґрунтуючись на наведеному на рисунках 1.2 поділу (Lamar University).

Малюнок 1.2. Типи стійкості рівноважних точок

Розглянувши потрібні теоретичні відомості, перейдемо до аналізу систем.

Розглянемо систему без насичення:

Вона дуже проста та не підходить для практичного застосування, оскільки не має жодних обмежень. Але як перший приклад аналізу системи підходить для розгляду.

Спочатку знайдемо рівноважні точки, прирівнявши праві частини рівнянь до нуля. Таким чином, виявляємо дві рівноважні точки, назвемо їх A і B: .

Об'єднаємо крок із пошуком матриці Якобі, коренів характеристичного рівняння та визначенням типу стійкості. Оскільки вони елементарні, то одразу отримаємо відповідь:

1. У точці , , стійкий вузол.

У точці: , , сідло.

Як я вже писав, дана система занадто тривіальна, тому не потрібно ніяких пояснень.

Тепер проведемо аналіз системи з насичення:

(1.9)

Поява обмеження на взаємонасичення продукцією підприємствами наближає нас до реальної картини того, що відбувається, а також трохи ускладнює систему.

Як і раніше, прирівнюємо праві частини системи до нуля та вирішуємо отриману систему. Крапка залишилася без змін, а ось інша точка в даному випадку містить більше параметрів, ніж раніше: .

У цьому випадку матриця Якобі набуває такого вигляду:

Віднімемо від неї одиничну матрицю, помножену на і прирівняємо визначник отриманої матриці в точці А і B до нуля.

У точці аналогічну ранню картину:

Стійкий вузол.

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

2 Фазові портрети систем

Переважна більшість нелінійних динамічних моделей являють собою складні диференціальні рівняння, які або не вирішити, або це є деякою складністю. Прикладом може бути система з попереднього розділу. Незважаючи на простоту, знаходження типу стійкості в другій рівноважної точки було не легкою справою (нехай і не з математичної точки зору), а зі збільшенням параметрів, обмежень і рівнянь для збільшення кількості взаємодіючих підприємств, складність буде лише зростати. Звичайно, якщо параметри будуть являти собою числові вирази, то все стане неймовірно просто, але тоді аналіз певною мірою втратить будь-який сенс, адже в результаті, ми зможемо знайти рівноважні точки і дізнатися їх типи стійкості лише для конкретного випадку, а не загального.

У таких випадках, варто згадати про фазову площину і фазові портрети. У прикладній математиці, зокрема в контексті нелінійного системного аналізу, фазова площина є візуальним відображенням певних характеристик деяких видів диференціальних рівнянь (Nolte, 2015). Координатна площиназ осями значень будь-якої пари змінних, що характеризують стан системи – двовимірний випадок загального n-вимірного фазового простору.

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

Рішення диференціального рівняння є сімейством функцій. Графічно це можна побудувати у фазовій площині як двовимірне векторне поле. На площині малюються вектори, які представляють похідні у характерних точках за яким-небудь параметром, у разі по часу, тобто (). При достатній кількості цих стрілок в одній області можна візуалізувати поведінку системи та легко ідентифікувати граничні цикли (Boeing, 2016).

Векторне поле є фазовим портретом, конкретний шлях вздовж лінії потоку (тобто шлях, який завжди стосується векторів) є фазовим шляхом. Потоки у векторному полі вказують на зміну системи у часі, що описується диференціальним рівнянням (Jordan, 2007).

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

Отже, фазові площини корисні візуалізації поведінки фізичних систем. Зокрема, коливальних систем, таких як згадана вище модель хижак-жертва. У цих моделях фазові траєкторії можуть «закручуватися» у напрямі нуля, «виходити зі спіралі» в нескінченність або досягати нейтральної сталої ситуації, яка називається центрами. Це корисно при визначенні того, чи стабільна динаміка чи ні (Jordan, 2007).

Представлені в цьому розділі фазові портрети будуть побудовані за допомогою інструментів WolframAlpha, або наведені з інших джерел. Мальтузіанська модель зростання без насичення.

Побудуємо фазовий портрет першої системи з трьома набором параметрів, щоб порівняти їхню поведінку. Набір А ((1,1), (1,1)), який надалі називатиметься одиничним набором, набір B ((10,0.1), (2,2)), при виборі якого в системі спостерігається різкий спад виробництва продукції , і набір C ((1,10), (1,10)), при якому навпаки виникає різке і необмежене зростання. Варто відзначити, що значення по осях у всіх випадках будуть знаходитися в одних і тих же інтервалах від -10 до 10, для зручності порівняння фазових діаграм між собою. Звичайно, це не стосується якісного портрета системи, у якого осі безрозмірні.

Малюнок 1.3 Фазовий портрет із параметрами А

мутуалізм диференціальний граничний рівняння

На малюнку 1.3, наведених вище, продемонстровані фазові портрети системи при трьох зазначених наборах параметрів, а також фазовий портрет, що описує якісну поведінку системи. Не варто забувати, що найважливішою з практичної точки зору є перша чверть, оскільки кількість продукції, яка може бути лише невід'ємною, є наші осі.

На кожному з малюнків явно видно стійкість у рівноважній точці (0,0). І першому малюнку також помітно «сідло» у точці (1,1), інакше кажучи, якщо підставити значення набору параметрів у систему, то рівноважної точці У. При зміні меж побудови моделі, сідлова точка виявляється і інших фазових портретах.

Мальтузіанська модель зростання з насичення.

Побудуємо фазові діаграми для другої системи, в якій є насичення, з трьома новими наборами значень параметрів. Набір А, ((0.1,15,100), (0.1,15,100)), набір ((1,1,0.5), (1, 1,0.5)) і набір С ((20,1,100), (20,1,100) )).

Малюнок 1.4. Фазовий портрет із параметрами А

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

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

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

Інтегральні криві відомі під іншими назвами, залежно від природи та інтерпретації диференціального рівняння або векторного поля. У фізиці інтегральні криві для електричного поля або магнітного полявідомі як лінії поля, а інтегральні криві для поля швидкостей рідини відомі як лінії струму. У динамічних системах інтегральні криві для диференціального рівняння називають траєкторіями.

Малюнок 1.5. Інтегральні криві

Рішення будь-якої системи можна розглядати і як рівняння інтегральних кривих. Очевидно, що кожна фазова траєкторія - проекція деякої інтегральної кривої в просторі x,y,tна фазову поверхню.

Для побудови інтегральних кривих є кілька способів.

Один з них – метод ізоклін. Ізокліна - це крива, що проходить через точки, в яких нахил функції буде завжди одним і тим же, незалежно від початкових умов (Hanski, 1999).

Він часто використовується як графічний метод вирішення звичайних диференціальних рівнянь. Наприклад, у рівнянні виду y"= f(x, y) ізокліни є лініями на площині (x, y), отриманими прирівнюванням f (x, y) до константи. Це дає ряд ліній (для різних констант), вздовж яких криві рішення мають один і той же градієнт, обчислюючи цей градієнт для кожної ізокліни, поле нахилу можна візуалізувати, що дозволяє порівняно легко намалювати наближені криві рішення.

Малюнок 1.6. Метод ізоклін

Цей метод не вимагає обчислень на комп'ютері, і був дуже популярним раніше. Зараз існують програмні рішення, які на комп'ютерах побудують інтегральні криві гранично точно та швидко. Однак, навіть так метод ізоклін непогано зарекомендував себе як інструмент для вивчення поведінки рішень, оскільки дозволяє показати області типової поведінки інтегральних кривих.

Мальтузіанська модель зростання без насичення.

Почнемо з того, що, незважаючи на існування різних способів побудови, показати інтегральні криві системи рівнянь не так просто. Метод ізоклін, що згадується раніше, не підходить, оскільки він працює для диференціальних рівнянь першого порядку. А програмні засоби, що мають можливість побудови таких кривих, не знаходяться у відкритому доступі. Наприклад, Wolfram Mathematica, здатна цього, платна. Тому намагатимемося максимально використовувати можливості Wolfram Alpha, робота з якою описана в різних статтях та роботах (Orca, 2009). Навіть незважаючи на те, що картина буде явно не зовсім достовірною, але принаймні дозволить показати залежність у площинах (x, t), (y, t). Для початку розв'яжемо кожне з рівнянь щодо t. Тобто виведемо залежність кожної із змінних щодо часу. Для цієї системи отримуємо:

(1.10)

(1.11)

Рівняння симетричні, тому розглянемо лише одне їх, саме x(t). Нехай константа дорівнює 1. У разі скористаємося функцією побудови графіків.

Малюнок 1.7. Тривимірна модель для рівняння (1.10)

Мальтузіанська модель зростання з насичення.

Виконаємо аналогічні дії і для іншої моделі. Зрештою отримуємо два рівняння, що демонструють залежність змінних від часу.

(1.12)

(1.13)

Знову побудуємо тривимірну модель та лінії рівня.

Малюнок 1.8. Тривимірна модель для рівняння (1.12)

Оскільки значення змінних невід'ємні, то дроби при експоненті отримуємо від'ємне число. Таким чином згодом інтегральна крива зменшується.

Раніше було дано визначення системної динаміки для розуміння суті роботи, тепер зупинимося на цьому детальніше.

Системна динаміка - методологія та метод математичного моделюваннядля формування, розуміння та обговорення складних проблем, спочатку розроблена в 1950-х роках Джеєм Форрестером, та описана в його роботі (Forrester, 1961).

Системна динаміка одна із аспектів теорії систем як методу розуміння динамічного поведінки складних систем. Основою методу є визнання того, що структура будь-якої системи складається з численних відносин між її компонентами, які найчастіше настільки важливі визначення її поведінки, як і самі окремі компоненти. Прикладами є теорія хаосу та соціальна динаміка, описані в роботах різних авторів (Grebogi, 1987; Sontag, 1998; Кузнєцов, 2001; Табор, 2001). Також стверджується, що, оскільки у властивостях елементів часто не можуть бути знайдені властивості-цілого, у деяких випадках поведінка цілого не може бути пояснена з погляду поведінки частин.

Моделювання може по-справжньому показати всю практичну значимістьдинамічної системи. Хоча воно і можливе в електронних таблицях, існує безліч програмних пакетів, які були оптимізовані спеціально для цього.

Саме собою моделювання - це процес створення та аналізу прототипу фізичної моделі для прогнозування її продуктивності у світі. Імітаційне моделювання використовується, щоб допомогти проектувальникам та інженерам зрозуміти, за яких умов і в яких випадках процес може зазнати невдачі та які навантаження він може витримати (Хемді, 2007). Моделювання також може допомогти передбачити поведінку потоків рідини та інші фізичні явища. У моделі аналізується приблизні умови роботи за рахунок застосовуваних імітаційних програмних засобів (Строгалев, 2008).

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

Проте можливості моделювання не безмежні. Насамперед, це пов'язано з тим, що важко оцінити сферу застосування імітаційної моделі, зокрема, період часу, для якого прогноз може бути побудований з необхідною точністю (Law, 2006). Крім того, за своєю природою імітаційна модель прив'язана до конкретного об'єкта, а при спробі застосувати до іншого, навіть аналогічного йому об'єкта, вимагає радикального коригування або принаймні суттєвої модифікації.

Існує загальна причина існування обмежень на імітаційне моделювання. Побудова і чисельний розрахунок «точної» моделі успішний лише за існуванні кількісної теорії, тобто у тому разі, якщо всі рівняння відомі, а завдання зводиться лише вирішення даних рівнянь з певною точністю (Базыкин, 2003).

Але навіть незважаючи на це, імітаційне моделювання – чудовий засіб візуалізації динамічних процесів, що дозволяє, за більш-менш правильної моделі, приймати рішення, засновані на її результатах.

У цій роботі моделі систем будуть побудовані за допомогою засобів системної динаміки, які пропонуються програмою AnyLogic.

Мальтузіанська модель зростання без насичення/

Перед побудовою моделі необхідно розглянути елементи системної динаміки, якими ми користуватимемося, та зв'язати їх із нашою системою. Наступні визначення були взяті із довідкової інформації програми AnyLogic.

Накопичувач – основний елемент діаграм системної динаміки. Вони застосовуються уявлення об'єктів реального світу, у яких накопичуються деякі ресурси: гроші, речовини, чисельності груп людей, деякі матеріальні об'єкти тощо. Накопичувачі відбивають статичний стан моделируемой системи, які значення змінюються згодом у відповідність із існуючими у системі потоками. Звідси випливає, що динаміку системи задають потоки. Вхідні та вихідні з накопичувача потоки збільшують або зменшують значення накопичувача.

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

Поки накопичувачі визначають статичну частину системи, потоки визначають швидкість зміни значень накопичувачів, тобто у часі відбуваються зміни запасів і, таким чином, визначають динаміку системи.

Агент може містити змінні. Змінні зазвичай використовуються для моделювання змінних характеристик агента або зберігання результатів роботи моделі. Зазвичай динамічні змінні складаються з функції накопичувачів.

Агент може мати параметри. Параметри часто використовуються для представлення деяких характеристик змодельованого об'єкта. Вони корисні, коли екземпляри об'єктів мають однакову поведінку, описану у класі, але відрізняються деякими значеннями параметрів. Між змінними та параметрами існує явна різниця. Змінна є станом моделі і може змінюватися під час моделювання. Параметр зазвичай використовується для статичного опису об'єктів. Під час одного прогону моделі параметр зазвичай є константою і змінюється тільки тоді, коли потрібно переналаштувати поведінку моделі.

Зв'язок - елемент системної динаміки, що використовується для визначення залежності між елементами діаграми потоків і накопичувачів. Як приклад, якщо якийсь елемент A згадується в рівнянні або початковому значенні елемента B, спочатку необхідно з'єднати ці елементи зв'язком, що йде від A до B, і тільки потім ввести вираз у властивостях B.

Існують і деякі інші елементи системної динаміки, але вони не будуть задіяні в ході роботи, тому опустимо їх.

Для початку розглянемо з чого складатиметься модель системи (1.4).

По-перше, відразу відзначаємо два накопичувачі, які будуть містити в собі значення кількості продукції кожного з підприємств.

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

По-третє, переходимо до змінних та параметрів. Змінних лише дві. X та Y, відповідальні за зростання продукції. А також у нас є чотири параметри.

По-четверте, щодо зв'язків, кожен із потоків повинен бути пов'язаний зі змінними та параметрами, що входять до рівняння потоку, а також обидві змінні повинні мати зв'язок із накопичувачами для зміни значення з часом.

Детальний опис побудови моделі, як приклад роботи в середовищі моделювання AnyLogic, залишимо для наступної системи, оскільки вона дещо складніша і в ній використовується більше параметрів, і одразу перейдемо до розгляду готового варіанта системи.

Нижче на малюнку 1.9 представлена ​​побудована модель:

Малюнок 1.9. Модель системної динаміки для системи (1.4)

Усі елементи системної динаміки відповідають описаному вище, тобто. два накопичувачі, чотири потоки (два вхідних, два вихідних), чотири параметри, дві динамічні змінні, та необхідні зв'язки.

На малюнку видно, що більше продукції, тим сильніше її зростання, що призводить до різкого зростання кількості товарів, що відповідає нашій системі. Але, як вже було раніше сказано, відсутність обмеження на це зростання не дозволяють застосовувати дану модель на практиці.

Мальтузіанська модель зростання з насичення/

Розглядаючи цю систему, докладніше зупинимося на побудові моделі.

Першим кроком додаємо два накопичувачі, назвемо їх X_stock та Y_stock. Кожному з них поставимо початкове значення рівне 1. Зазначимо, що відсутність потоків в класично заданому рівняннянакопичувача нічого немає.

Малюнок 1.10. Побудова моделі системи (1.9)

Наступний крок – додавання потоків. Побудуємо для кожного накопичувача вхідний та вихідний потік за допомогою графічного редактора. Не можна забувати, що один із країв потоку повинен перебувати в накопичувачі, інакше вони не будуть пов'язані.

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

Третім кроком у нас є додавання шести параметрів та двох динамічних змінних. Дамо кожному елементу ім'я у відповідність до його буквеним виразом у системі, а також задаємо початкові значення параметрів наступним чином: e1=e2=1, a12=a21=3, n1=n2=0,2.

Усі елементи рівнянь присутні, залишилося лише написати рівняння потоків, але цього спочатку необхідно додати зв'язки між елементами. Наприклад, вихідний потік, відповідальний за доданок , може бути пов'язані з e1 і x. А кожна динамічна змінна повинна бути пов'язана з відповідним їй накопичувачем (X_stock x, Y_stock y). Створення зв'язків відбувається аналогічно додавання потоків.

Після створення необхідних зв'язків можна переходити до написання рівнянь потоків, що демонструється правому малюнку. Звичайно, можна піти і в зворотному порядку, але при існуванні зв'язків, під час написання рівнянь з'являються підказки для встановлення потрібних параметрів/змінних, що полегшує завдання в складних моделях.

Після виконання всіх кроків можна запускати імітаційну модель і подивитися на її результат.

Розглянувши системи нелінійних диференціальних рівнянь взаємодії підприємств за умов мутуалізму, можна зробити кілька висновків.

Існує два стани системи: різке необмежене зростання, або прагнення кількості продукції до нуля. Який із двох станів прийме система залежить від параметрів.

Жодна із запропонованих моделей, у тому числі модель з урахуванням насичення, не підходить для практичного застосування через відсутність ненульового стійкого положення, а також причин, описаних у пункті 1.

У разі спроби подальшого дослідження даного типу симбіотичної взаємодії для створення моделі, що застосовується компаніями на практиці, необхідне подальше ускладнення системи та введення нових параметрів. Наприклад, Базикін у своїй книзі наводить приклад динаміки двох мутуалістичних популяцій із запровадженням додаткового чинника внутрішньовидової конкуренції. За рахунок чого система набуває вигляду:

(1.15)

І в такому разі з'являється ненульове стійке положення системи, відокремлене від нульового «сідлом», що наближає її до реальної картини того, що відбувається.

2. Взаємодія підприємств в умовах протокооперації

Всі основні теоретичні відомості були представлені в попередньому розділі, тому при аналізі моделей, що розглядаються в цьому розділі, теорія здебільшого буде опущена, за винятком кількох моментів, з якими ми не стикалися у попередньому розділі, а також можливе скорочення обчислень. Розглянута в цьому розділі модель взаємодії організацій в умовах протокооперації, що складається з двох рівнянь, заснованих на мальтузіанської моделі, виглядає як система (1.5). Проаналізовані в попередньому розділі системи показали, що для їх максимального наближення до діючих моделей необхідне ускладнення систем. Виходячи з даних висновків, відразу додамо на модель обмеження на зростання. На відміну від попереднього типу взаємодії, коли зростання, що не залежить від іншої компанії, негативний, у цьому випадку всі знаки позитивні, а отже, маємо постійне зростання. Уникаючи недоліків, описаних раніше, постараємося обмежити його логістичним рівнянням, також відомим як рівняння Ферхюльста (Gershenfeld, 1999), що має такий вигляд:

, (2.1)

де P – чисельність популяції, r – параметр, що показує швидкість зростання, K – параметр, що відповідає за максимально можливу чисельність популяції. Тобто з часом чисельність популяції (у разі продукції), йтиме до певного параметра До.

Дане рівняння допоможе стримати нестримне зростання продукції, яке ми спостерігали раніше. Таким чином система набуває наступного вигляду:

(2.2)

Не варто забувати, що обсяг товару, що зберігається на складі для кожної компанії різний, тому параметри, що обмежують зростання різні. Назвемо цю систему «», і надалі використовуватимемо цю назву, коли її розглядатимемо.

Другою системою, яку ми розглядатимемо, є подальший розвитокмоделі з обмеженням Ферхюльста Як і в попередньому розділі введемо обмеження на насичення, тоді система набуде вигляду:

(2.3)

Тепер кожен із доданків має власне обмеження, тому вже без подальшого аналізу можна помітити, що необмеженого зростання, як у моделях попереднього розділу, не буде. А оскільки кожен із доданків демонструє позитивне зростання, то й кількість продукції не впаде в нуль. Назвемо цю модель "модель протокооперації з двома обмеженнями".

Дані дві моделі розглядаються в різних джерелах про біологічні популяції. Тепер спробуємо дещо розширити системи. Для цього розглянемо наступний рисунок.

На малюнку продемонстровано приклад процесів двох компаній: сталеливарної та вугільної промисловості. В обох підприємствах є зростання продукції, що не залежить від іншої, а також є зростання продукції, яке виходить завдяки їхній взаємодії. Це ми вже враховували у ранніх моделях. Тепер варто звернути увагу, що фірми як виробляють продукцію, вони її ще й продають, наприклад, ринку чи взаємодіє з нею підприємства. Тобто. виходячи з логічних висновків, існує необхідність негативного зростання компаній за рахунок продажу продукції (на малюнку за це відповідають параметри ?1 і ?2), а також за рахунок передачі частини продукції іншому підприємству. Раніше ми враховували це лише з позитивним знаком у іншої компанії, але не розглядали те, що у першого підприємства під час передачі продукції її кількість зменшується. У такому разі отримуємо систему:

(2.4)

І якщо про доданок можна сказати, що якби в попередніх моделях було зазначено, що , характеризують природний приріст, а параметр може бути негативним, то різниці практично немає, то про доданок такого сказати не можна. До того ж надалі, розглядаючи подібну систему із запровадженим нею обмеженням, правильніше використовувати саме доданки позитивного і негативного зростання, оскільки у разі ними можуть накладатися різні обмеження, що неможливо для природного приросту. Назвемо її "розширена модель протокооперації".

І нарешті, четвертою моделлю, що розглядається, є розширена модель протокооперації з раніше згаданим логістичним обмеженням на зростання. І система для цієї моделі така:

, (2.5)

де - приріст продукції першого підприємства, який не залежить від другого, з урахуванням логістичного обмеження, - приріст продукції першої компанії, що залежить від другої, з урахуванням логістичного обмеження; - приріст продукції другого підприємства, що не залежить від першого, з урахуванням логістичного обмеження, - приріст продукції другої компанії, що залежить від першої, з урахуванням логістичного обмеження; - споживання товарів першого підприємства, не пов'язане з іншим; - споживання товарів другого підприємства, не пов'язане з іншим; - споживання товарів першої галузі другою галуззю; - споживання товарів другої галузі. першою галуззю.

Надалі ця модель позначатиметься як «розширена модель протооперації з логістичним обмеженням».

1 Стійкість систем у першому наближенні

Модель протокооперації з обмеженням Ферхюльста

Методи аналізу стійкості системи були зазначені в аналогічному розділі попереднього розділу. Насамперед знаходимо рівноважні точки. Одна з них, як завжди, нульова. Інша є точку з координатами .

Для нульової точки λ1 = , λ2 = , оскільки обидва параметри невід'ємні, отримуємо нестійкий вузол.

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

Аналіз даної системи складніший за попередню в силу того, що додається фактор насичення, таким чином з'являються нові параметри, а при знаходженні рівноважних точок доведеться вирішувати не лінійне, а білінійне рівняння через змінну в знаменнику. Тому, як і попередньому випадку, залишимо визначення типу стійкості на фазові діаграми.

Незважаючи на появу нових параметрів, Якобіан у нульовій точці, так само, як і коріння характеристичного рівняння, виглядає аналогічно до попередньої моделі. Таким чином, у нульовій точці є нестійкий вузол.

Перейдемо до розширених моделей. Перша їх не містить жодних обмежень і набуває вигляду системи (2.4)

Зробимо заміну змінних, , і . Нова система:

(2.6)

У такому разі отримуємо дві рівноважні точки, точка А(0,0), В(). Точка лежить у першій чверті, оскільки змінні мають невід'ємне значення.

Для рівноважної точки А отримуємо:

. - Нестійкий вузол,

. - Сідло,

. - Сідло,

. - стійкий вузол,

У точці B коріння характеристичного рівняння є комплексними числами: λ1 = , λ2 = . Ми не можемо визначити тип стійкості покладаючись на теореми Ляпунова, тому проведемо чисельне моделювання, яке не покаже всі можливі стани, але дозволить дізнатися хоча б деякі з них.

Малюнок 2.2. Чисельне моделювання пошуку типу стійкості

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

Не вдаючись у подробиці обчислень, приходимо до наступних рівноважних точок. Точка А(0,0) та точка В з наступними координатами:

(), де a =

Для точки А визначення типу стійкості - тривіальна задача. Коріння характеристичного рівняння таке λ1 = , λ2 = . Таким чином отримуємо чотири варіанти:

1. λ1 > 0, λ2 > 0 - нестійкий вузол.

2. λ1< 0, λ2 >0 – сідло.

3. λ1 ​​> 0, λ2< 0 - седло.

4. λ1< 0, λ2 < 0 - устойчивый узел.

Говорячи про точку В, варто погодитися, що підстановка скорочень у вираз для неї ускладнить роботу з Якобіаном і знаходженням коренів характеристичного рівняння. Наприклад, після спроби їх пошуку за допомогою обчислювальних засобів WolframAlpha, виведення значення коренів зайняв близько п'яти рядків, що не дозволяє працювати з ними буквально. Звичайно, за наявності вже наявних параметрів, можна швидко знайти рівноважну точку, але це окремий випадок, оскільки ми знайдемо стан рівноваги, якщо воно є, лише для даних параметрів, що не підходить для системи підтримки прийняття рішень, для якої модель і планується створюватися .

Через складність роботи з корінням характеристичного рівняння побудуємо взаємне розташування нуль-ізоклін за аналогією з розібраною в роботі Базикиною системою (Базикін, 2003). Це дозволить нам розглянути можливі стани системи, і надалі при побудові фазових портретів виявити рівноважні точки та типи їхньої стійкості.

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

(2.7)

Таким чином, ізокліни мають вигляд парабол.

Малюнок 2.3. Можливий варіант розташування нуль-ізоклін

Усього можливі чотири випадки їхнього взаємного розташування за кількістю загальних точок між параболами. Для кожного з них є свої набори параметрів, а значить і фазові портрети системи.

2 Фазові портрети систем

Побудуємо фазовий портрет системи, за умови, що а інші параметри дорівнюють 1. В даному випадку достатньо і одного набору змінних, оскільки якісна не зміниться.

Як видно з наведених нижче малюнків, нульова точка - нестійкий вузол, а друга точка, якщо підставити числові значення параметрів, то отримаємо (-1.5, -1.5) - сідло.

Малюнок 2.4. Фазовий портрет для системи (2.2)

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

Модель протокооперації із двома обмеженнями.

У цій системі є додатковий стримуючий фактор, тому фазові діаграми повинні відрізнятися від попереднього випадку, що і видно на малюнку. Нульова точка також - нестійкий вузол, але у цій системі з'являється стійке становище, саме стійкий вузол. При даних параметрах координати (5.5,5.5), він представлений малюнку.

Малюнок 2.5. Фазовий портрет для системи (2.3)

Таким чином, обмеження на кожен доданок дозволило отримати стійке положення системи.

Розширена модель протокооперації.

Побудуємо фазові портрети для розширеної моделі, але одразу використовуючи її модифікований вигляд:

Розглянемо чотири набори параметрів, причому таких, щоб розглянути всі випадки з нульовою рівноважною точкою, а також продемонструвати фазові діаграми чисельного моделювання, що використовується для рівномірної ненульової точки: набір А(1,0.5,0,5) відповідає стану , набір (1,0.5,-0.5) відповідає набір С(-1,0.5,0,5) а набір D(-1,0.5,-0,5) , тобто стійкого вузла в нульовій точці. Перші два набори продемонструють фазові портрети для параметрів, які ми розглядали при чисельному моделюванні.

Малюнок 2.6. Фазовий портрет системи (2.4) з параметрами А-D.

На малюнках необхідно звернути увагу на точки (-1,2) та (1,-2) відповідно, у них виникає «сідло». Для більш детального уявлення, малюнку представлений інший масштаб малюнка з сідловою точкою (1,-2). На малюнку в точках (1,2) та (-1,-2) видно стійкий центр. Що стосується нульової точки, то починаючи з малюнка по малюнок на фазових діаграмах чітко помітний нестійкий вузол, сідло, сідло та стійкий вузол.

Розширена модель протокооперації із логістичним обмеженням.

Як і в попередній моделі, продемонструємо фазові портрети для чотирьох випадків нульової точки, а також постараємося відзначити і ненульові рішення на цих діаграмах. Для цього візьмемо наступні набори параметрів із параметрами, вказаними в наступному порядку (): А(2,1,2,1), Б(2,1,1,2), С(1,2,2,1) та Д (1,2,1,2). Інші параметри для всіх наборів будуть наступними: , .

На поданих нижче малюнках можна спостерігати чотири рівноважні стани нульової точки, описаних у попередньому розділі для даної динамічної системи. А також на малюнках є стійке положення точки з однією ненульовою координатою.

Малюнок 2.7. Фазовий портрет системи (2.5) з параметрами А-B

3 Інтегральні траєкторії систем

Модель протокооперації з обмеженням Ферхюльста

Як і в попередньому розділі розв'яжемо кожне з диференціальних рівнянь окремо і явно виразимо залежність змінних від часового параметра.

(2.8)

(2.9)

З отриманих рівнянь видно, що значення кожної зі змінних зростає, що демонструється на тривимірної моделі нижче.

Малюнок 2.8. Тривимірна модель для рівняння (2.8)

Даний вид графіка на початку чимось нагадує тривимірне зображення мальтузіанської моделі без насичення, що розглядається в главі 1, оскільки має аналогічне їй швидке зростання, але надалі можна помітити зниження швидкості зростання через досягнення обмеження на обсяг продукції. Таким чином, підсумковий зовнішній вигляд інтегральних кривих схожий з графіком логістичного рівняння, яке було використане, щоб обмежити один із доданків.

Модель протокооперації із двома обмеженнями.

Вирішуємо кожне з рівнянь за допомогою Wolfram Alpha. Таким чином, залежність функції x(t) зводиться до наступного виду:

(2.10)

Для другої функції ситуація аналогічна, тому опустимо її розв'язання. Чисельні значення виникли через заміну властивостей деякими відповідними їм значеннями, що не впливає на якісну поведінку інтегральних кривих. На наведених нижче малюнках помітно використання обмежень на зростання, оскільки згодом експоненційне зростання перетворюється на логарифмічний.

Малюнок 2.9. Тривимірна модель для рівняння (2.10)

Розширена модель протокооперації

Майже аналогічно моделям при мутуалізмі. Єдина різниця у швидшому щодо тих моделей зростанні, що видно з наведених нижче рівнянь (якщо подивитися на ступінь експоненти) та графіків. Інтегральна крива повинна набувати вигляду експоненти.

(2.11)

(2.12)

Розширена модель протокооперації з логістичним обмеженням

Залежність x(t) виглядає так:

Без графіка складно оцінити, поведінка функції, тому скориставшись відомими нам засобами, побудуємо його.

Малюнок 2.10 Тривимірна модель рівняння

Значення функції зменшується при не малих значеннях іншої змінної, що пов'язано з відсутністю обмежень на негативний білінійний доданок, і є очевидним результатом

4 Системна динаміка взаємодіючих компаній

Модель протокооперації з обмеженням Ферхюльста.

Побудуємо систему (2.2). Використовуючи вже відомі нам інструменти, будуємо імітаційну модель. Цього разу на відміну від мутуалістичних моделей, у моделі буде логістичне обмеження.

Малюнок 2.11. Модель системної динаміки для системи (2.2)

Запустимо модель. У цій моделі варто відзначити той факт, що зростання взаємозв'язку нічим не обмежене, а зростання продукції без впливу іншого має специфічне обмеження. Якщо подивитися на вираз логістичної функції, можна помітити, що у разі, коли змінна (кількість товарів) перевищує максимально можливий обсяг зберігання, доданок стає негативним. У разі, коли існує лише логістична функція, таке неможливе, але при додатковому завжди позитивному факторі зростання таке можливо. І зараз важливо зрозуміти, що логістична функція впорається із ситуацією не надто швидкого зростання кількості продукції, наприклад, лінійного. Звернімо увагу на малюнки нижче.

Малюнок 2.12. Приклад роботи моделі системної динаміки системи (2.2)

На лівому малюнку показано 5 крок роботи програми відповідної запропонованої моделі. Але зараз варто звернути увагу на правий малюнок.

По-перше, для одного з вхідних потоків для Y_stock видалено зв'язок з х, виражена в доданку . Це зроблено для того, щоб показати різницю в роботі моделі при лінійному позитивному потоці, і білінійному зростанні, який представлений для X_stock. При лінійних необмежених потоках після перевищення параметра До система в якийсь момент приходить до рівноваги (у даній моделі рівноважний стан - 200 тисяч одиниць товару). Але набагато раніше білінійне зростання призводить до різкого зростання кількості товару, що переходить у нескінченність. Якщо ж залишити і правий і лівий постійно позитивні потоки білінійними, то вже приблизно на 20-30 кроці, значення накопичувача приходить до різниці двох нескінченностей.

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

Модель протокооперації із двома обмеженнями.

З'ясувавши недоліки минулої моделі і ввівши обмеження на другий доданок чинником насичення, побудуємо і запустимо нову модель.

Малюнок 2.13. Модель системної динаміки та приклад її роботи для системи (2.3)

Ця модель, зрештою, приносить довгоочікувані результати. Вийшло обмежити зростання значень накопичувача. Як очевидно з правого малюнка обох підприємств рівновагу досягається при невеликому перевищенні обсягу зберігання.

Розширена модель протокооперації.

При розгляді системної динаміки цієї моделі буде продемонстровано можливості програмного середовища AnyLogic для яскравої візуалізації моделей. Усі попередні моделі були побудовані лише з використанням елементів системної динаміки. Тому самі моделі виглядали непомітно, вони не дозволяли відстежити динаміку зміни кількості продукції в часі та змінювати параметри під час роботи програми. При роботі з цією та наступною моделями постараємося скористатися ширшим спектром можливостей програми для зміни трьох зазначених вище недоліків.

По-перше, у програмі поряд із розділом «системна динаміка» у програмі також присутні розділи «картинки», «3D-об'єкти», що дозволяють урізноманітнити модель, що корисно за її подальшої презентації, оскільки робить вигляд моделі «приємнішою».

По-друге, для відстеження динаміки зміни значень моделі існує розділ статистика, який дозволяє додавати діаграми та різні інструменти збору даних, пов'язуючи їх з моделлю.

По-третє, для зміни параметрів та інших об'єктів під час виконання моделі є розділ «елементи керування». Об'єкти цього розділу дозволяють змінювати параметри під час роботи моделі (приклад, «бігунок»), вибирати різні стани об'єкта (приклад, «перемикач») і виконувати інші дії, що змінюють спочатку задані дані під час роботи.

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

Розширена модель протокооперації із логістичним обмеженням.

Використовуючи готову попередню модель, додамо до неї параметри з логістичного рівняння для обмеження зростання.

Опустимо побудову моделі, оскільки на попередніх п'яти моделях, представлених у роботі, вже було продемонстровано всі необхідні інструменти та принципи роботи з ними. Варто лише відзначити, що її поведінка подібна до моделі протокооперації з обмеженням Ферхюльста. Тобто. відсутність насичення заважає її практичному застосуванню.

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

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

Відповідні моделі повинні мати обмеження на кожному з доданків, оскільки в іншому випадку різке зростання білінійних множників «руйнує» всю імітаційну модель.

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

Висновок

В результаті проведеного дослідження було проведено аналіз шести систем, що описують динаміку виробництва продукції підприємствами, що взаємно впливають одна на одну. У результаті рівноважні точки і типи їх стійкості були визначені одним із наступних способів: аналітично, або завдяки побудованим фазовим портретам у випадках, коли аналітичне рішення з будь-яких причин неможливо. Для кожної системи були побудовані фазові діаграми, а також побудовані тривимірні моделі, на яких при проектуванні можливо отримати інтегральні криві в площинах (x, t), (y, t). Після використання середовища моделювання AnyLogic були побудовані всі моделі і розглянуті варіанти їх поведінки при певних параметрах.

Після проведення аналізу систем і побудови їх імітаційних моделей стає очевидним, що дані моделі можуть розглядатися лише як навчальні, або ж для опису макроскопічних систем, але ніяк не як система підтримки прийняття рішень для окремих компаній, через свою низьку точність і в деяких місцях не зовсім достовірне уявлення процесів, що відбуваються. Але також не варто забувати, що якою б вірною не була описуюча модель динамічна система у кожної компанії/організації/галузі свої власні процеси та обмеження, таким чином створити і описати загальну модель неможливо. У кожному конкретному випадку вона видозмінюватиметься: ускладнюватиметься або навпаки спрощуватиметься для подальшої роботи.

Роблячи висновок з висновків до кожного розділу, варто загострити увагу на виявленому факті, що введення обмежень на кожне з доданків рівняння хоч і ускладнює систему, але дозволяє виявити стійкі положення системи, а також наблизити її до того, що відбувається насправді. І варто відзначити, що моделі протокооперації більше підходять для вивчення, оскільки мають ненульові стійкі положення на відміну від двох мутуалістичних моделей, які ми розглянули.

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

Література

1. Bhatia Nam Parshad; Szegх Giorgio P. (2002). Стабільність теорії динамічних систем. Springer.

2. Blanchard P.; Devaney, R. L.; Hall, G. R. (2006). Differential Equations. London: Thompson. pp. 96-111.

Boeing, G. (2016). Visual Analysis of Nonlinear Dynamical Systems: Chaos, Fractals, Self-Similarity and Limits of Prediction. Systems. 4(4): 37.

4. Campbell, David K. (2004). Nonlinear physics: Fresh breather. Nature. 432 (7016): 455-456.

Elton C.S. (1968) reprint. Animal ecology. Great Britain: William Clowes and Sons Ltd.

7. Forrester Jay W. (1961). Industrial Dynamics. MIT Press

8. Gandolfo, Giancarlo (1996). Economic Dynamics (Third ed.). Берлін: Springer. pp. 407-428.

9. Gershenfeld Neil A. (1999). The Nature of Mathematical Modeling. Cambridge, UK: Cambridge University Press.

10. Goodman M. (1989). Study Notes in System Dynamics. Pegasus.

Grebogi C, Ott E, та Yorke J. (1987). Chaos, Strange Attractors, і Фрактальні Basin Boundaries в Nonlinear Dynamics. Science 238 (4827), pp 632-638.

12. Hairer Ernst; Nørsett Syvert Paul; Wanner, Gerhard (1993), Solving ordinary differential equations I: Nonstiff problems, Berlin, New York

Hanski I. (1999) Metapopulation Ecology. Oxford University Press, Oxford, pp. 43-46.

Hughes-Hallett Deborah; McCallum, William G.; Gleason, Andrew M. (2013). Calculus: Single and Multivariable (6 ed.). Джон wiley.

15. Libro J., Valls C. (2007). Global analytic перші integrals для реального планування Lotka-Volterra system, J. Math. Phys.

16. Jordan D.W.; Smith P. (2007). Non-Linear Ordinary Differential Equations: Introduction for Scientists and Engineers (4th ed.). Oxford University Press.

Khalil Hassan K. (2001). Nonlinear Systems. Prentice Hall.

Lamar University, Online Math Notes - Phase Plane, P. Dawkins.

Lamar University, Online Math Notes - Systems of Differential Equations, P. Dawkins.

Lang Serge (1972). Різні manifolds. Reading, Mass.-London-Don Mills, Ont.: Addison-Wesley Publishing Co., Inc.

Law Averill M. (2006). Simulation Modeling and Analysis with Expertfit Software. McGraw-Hill Science.

Lazard D. (2009). Тридцяти років з Polynomial System Solving, and now? Journal of Symbolic Computation. 44 (3): 222-231.

24. Lewis Mark D. (2000). The Promise of Dynamic Systems Approaches для Integrated Account of Human Development. Child Development. 71 (1): 36-43.

25. Malthus T.R. (1798). An Essay на Принципі популяції, в Оксфордському світі з Classics reprint.

26. Morecroft John (2007). Strategic Modelling and Business Dynamics: A Feedback Systems Approach. John Wiley & Sons.

27. Nolte D.D. (2015), Introduction to Modern Dynamics: Chaos, Networks, Space and Time, Oxford University Press.

Автоматика та телемеханіка, Л-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) = ±<ж для любого фиксированного у е 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 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 роки з QBF Evaluations: QSAT є 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 з Event-Based Dynamics: Recent Developments in Analysis and Synthesis Methods. Маріо 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 Базовані на Microservices and Agent Networks. Proc. З 41th Intern. Конвенція про інформаційну і комунікаційну технологію, електронні та мікроелектроніки (MIPRO-2018). Р. 1643-1648.

14. Богданова В.Г., Горський С.А. Scalable Parallel Solver of Boolean Satisfiability Problems. Proc. З 41th Intern. Convention on Information and Communication Technology. Electronics and Microelectronics (MIPRO-2018). Р. 244-249.

15. Бичков I.V., Oparin G.A., Bogdanova V.G., Pashinin A.A. The Applied Problems Solving Technology Based on Distributed Computational Subject 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 - множини початкових, допустимих та цільових станів.

p align="justify"> Основними синтаксичними елементами логічної формули динамічної властивості є: 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*) у форматі КНФ:

Використовуючи закон Де-Моргана, знайдемо заперечення функції 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). Пропонуємо до вашої уваги журнали, що видаються у видавництві «Академія Природознавства»