EcoSimulation — 05. Імітація обмеження доступної кількості ресурсів


Підручник за курсом
«Імітаційне моделювання надорганізмових систем (з використанням LibreOffice Calc)»

Д.А. Шабанов

Врахування демографічної структури популяцій Імітація обмеження доступної кількості ресурсів Аналітичні вставки в імітаційній моделі
Імітаційне моделювання
надорганізмових систем-04 
Імітаційне моделювання надорганізмових систем-05 Імітаційне моделювання
надорганізмових систем-06

 

5. Імітація обмеження доступної кількості ресурсів

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

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

Якщо особини в популяції відрізняються одна від одної за потрібною для них кількістю ресурсів, нам доведеться інакше задавати обмеження, які середовище накладає на зростання чисельності популяції. Ми будемо задавати це обмеження не через місткість середовища K, задану так, що загальна чисельність особин після скорочення чисельності ωN не повинна перевищувати K (тобто ωN≤K), а як забезпеченість ресурсами (volume). У такому випадку (відповідно до прийнятих нами правил) нам знадобиться ще два позначення. Використовуємо букву D (demand) для позначення потреби в ресурсах, а букву U (use) – для позначення спожитої кількості ресурсів. Швидше за все, для новонароджених самиць і самців потреба в ресурсах має бути однаковою: df0=dm0; а крім цих величин, доведеться задати ще й df1, df2 тощо. Потреба в ресурсах може обчислюватися для чисельності модельної популяції на різних етапах: для αN або для ωN; відповідно до цього можна обчислювати αD або ωD.

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

Припустимо, в нашій моделі початкова чисельність популяції буде позначена як αN, альфа-чисельність. Після скорочення альфа-чисельності, викликаного різною виживаністю особин різних груп, ми перейдемо до ωN, омега-чисельності. Якщо оселище забезпечує усю потрібну кількість ресурсів для популяції (тобто ωDV), то споживання популяції відповідає її потребам (ωD=U). Якщо потреби популяції перевищують забезпеченість ресурсами, вона стикається з певними обмеженнями.

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

  • ωD≥V, тобто потреба в ресурсах вже наявних особин перевершує кількість ресурсів у середовищі. Приплоду немає.
  • ωD<V, но (ωD+nf0*df0+nm0*dm0)>V, тобто для наявних особин ресурсів достатньо, але якщо до наявних особин додати все можливе потомство, ресурсів не вистачатиме. У цьому випадку кількість потомства визначається як (V-ωD)/df0.
  • (ωD+nf0*df0+nm0*dm0)≤V, тобто ресурсів вистачає і для наявних особин, і для всіх їх можливих нащадків. В такому випадку приплід визначається кількістю потенційних батьків і їх плодючістю.

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

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

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

Ми розглядаємо популяцію або ГПС, геміклональну популяційну систему гібридогенного комплексу зелених жаб, як в роботі Кравченко, а може бути — гільдію, тобто сукупність різних популяцій, що використовують один і той же ресурс. У розглянутій сукупності представлені особини, що відносяться до різних форм, позначених g (genotypes). У попередньому прикладі форм дві: самиці (F) і самці (M), але, строго кажучи, їх може бути скільки завгодно; якщо їх кількість дорівнює k, різні генотипи — це g1, g2, g3 ... gk. Ці генотипи можуть належати до одного виду, як самиці та самці, або до різних видів, або до батьківських видів гібридогенного комплексу та до їх різноманітних гібридів. Представники кожного генотипу можуть бути представлені різними віковими групами. Якщо у системі, яку ми розглядаємо, представлено max_a вікових груп, починаючи з нульової (новонароджених), то генотип gg може бути представлений наступними віковими групами: gg0, gg1, gg2 ... ggmax_a.

Особливістю даного алгоритму є те, що смертність в ньому розділена на дві групи: неконкурентна смертність і конкурентна смертність. Мірою, що визначає неконкурентну смертність, є виживаність (s). Неконкурентна смертність не залежить від кількості ресурсів; конкурентна має місце тільки в тому випадку, коли спостерігається нестача ресурсів. Логічно неконкурентну смертність імітувати до того, як імітувати конкурентну смертність.

Для обговорюваного алгоритму серед інших початкових параметрів слід задавати значення виживаності для всіх груп: sg0, sg1, sg2 ... sgmax_a. Це можна зробити двома різними способами: таблицею та функцією. Перший шлях полягає в тому, щоб перерахувати кожне з значень sga і вказати у відповідному віконці моделі відповідну величину (в навчальній моделі, яку будуть робити студенти, найпростіше зробити саме так). Другий підхід полягає в тому, щоб задати певну функцію, що описує, як змінюється виживання представників кожного генотипу з віком. Конкретні значення sga будуть при цьому обчислюватися на підставі значення цієї функції.

Мірою, що визначає конкурентну смертність, є відносна конкурентоспроможність (с). Значення відносної конкурентоспроможності для всіх груп: сg0, сg1, сg2 ... сgmax_a задається так само, як і значення виживання.

Перетворення, які відбуваються з кожною групою особин gga, можна описати як перехід по ланцюжку станів αnga → βngaωnga. Приймемо, що альфа-чисельність — це чисельність на початку циклу, бета-чисельність — чисельність після неконкурентної смертності, а омега-чисельність — чисельність після конкурентного скорочення чисельності.

Отже, на початку кожного циклу чисельність кожної групи становить αnga. На початку циклу t початкова чисельність αtnga визначається остаточної чисельністю попереднього віку на циклі t-1: αtnga=ωt-1nga-1. У тому випадку, коли мова йде про потомство, що з'явилося на попередньому циклі роботи моделі (O), можна прийняти, що чисельність самиць і самців першого року життя приблизно дорівнює половині від загальної чисельності потомства минулого року. Обчислювати їх можна так: спочатку визначити чисельність самиць, а потім в залежності від неї визначити і чисельність самців: αtnf1~1/2*t-1O, а також αtnft-1O - αtnf1.

Потім імітується неконкурентна смертність, яка залежить від виживаності: βtnga = αtnga * sga.

На наступному кроці імітується конкурентне скорочення чисельності βtnf1ωtnf1. Звісно, найскладнішою частиною цієї моделі є алгоритм конкурентної смертності, обчислення омега-чисельності. Для того, щоб спростити розрахунки, ми можемо обчислювати проміжні величини: квотовані за конкурентоспроможністю чисельності всіх груп. Позначимо їх як гамма-чисельності, і задамо, що γtnf1 = βtnga * cga. Крім того, для подальших обчислень нам будуть потрібні ще дві величини, які необхідно розрахувати: сумарна потреба в ресурсах усіх особин до скорочення чисельності, βD = Σ(βtnga * dga), і сумарна потреба в ресурсах, що відповідає чисельності всіх квотованих груп, γD = Σ(γtnga * dga= Σ(βtng* cga * dga).

Після обчислення цих величин ми можемо встановити, за яким сценарієм буде проходити конкурентне скорочення чисельності. Позначимо сценарій як W (way — шлях, спосіб).

I сценарій: βD≤V; ресурсів всім вистачає, конкурентне скорочення чисельності не проводиться: ωnga = βnga.

II сценарій: βD>V і γD≥V; ресурсів не вистачає навіть після квотованого скорочення; в цьому випадку квотувану чисельність слід ще скоротити пропорційно V/γD: ωnga = γnga * V/γD. У разі, якщо після квотованого скорочення потреба усієї сукупності особин як раз дорівнює наявній кількості ресурсів, також можна застосовувати це формули (просто V/γD=1).

III сценарій: βD>V і γD<V; ресурсів не вистачає для всіх, але більше, ніж потрібно квотованим групам; в цьому випадку ωnga = βnga - (βngaγnga) * (βD-V)/(βD-γD).

Один з можливих варіантів розміщення такої моделі на аркуші LibreOffice Calc показаний на рис. 5.1. Чи треба копіювати цей варіант розміщення? Ні. Але моделі слід роботи так, щоб можна було зрозуміти, як вони побудовані. Не копіюйте цей варіант, зробіть краще!

Рис. 5.1. Загальний вигляд моделі з конкурентним скороченням чисельності при нестачі ресурсів. Показано формулу у комірці W26, яка описує конкурентне скорочення чисельності 

Природно, найменш інтуїтивно зрозумілим є третій сценарій. Щоб його зрозуміти, треба врахувати наступне. (βD-V)/(βD-γD)це відношення нестачі ресурсів для нескороченої чисельності до різниці в потреби в ресурсах між нескороченою (β) і квотованою (γ) чисельністю. Наприклад, бета-чисельність (нескорочена) потребує 500 одиниць ресурсу, гамма-чисельність (квотована) — 200 одиниць. Всього в наявності 400. Нестача ресурсів для бета-чисельності — 100 одиниць, різниця в потребах бета- і гамма-чисельностей — 300 одиниць. Відношення (βD-V)/(βD-γD) дорівнює 1/3. Значить, з бета-чисельності треба відняти третину різниці між бета- і гамма-чисельністю: 500 - 300/3 = 400.

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

1) чисельність до скорочення знижується до величини, яка відповідає забезпеченості ресурсами;

2) чисельність кожної групи до конкурентного скорочення βnga знижується до такої чисельності після конкурентного скорочення ωnga, що в кожній групі частка особин, які зберіглися після конкурентного скорочення, пропорційна конкурентоспроможності представників цієї групи: ωnga/βnga~cga. Після виконання цього алгоритму корисно обчислити остаточне споживання ресурсів: ωU = Σ(ωnga * dga). Якщо ωU відповідає V, то скорочення проведено правильно. Треба підкреслити, що відповідність ωU і V буде точним лише у тому разі, якщо обчислення проводилися без округлення, як це зроблено в моделі, показаної на рис. 5.1 (що, загалом-то, є далеким відходом від біологічного сенсу). В описаних розрахунках слід використовувати або скорочення вниз (в разі детермінованих моделей), або скорочення після додавання випадкового числа (в разі імовірнісних моделей).

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

5.3. Порівняння детерміністської та ймовірнісної моделей

Великий французький астроном та математик кінця XVIII – початку XIX століття П'єр-Симон Лаплас вважав, що світ детермінований своїм початковим станом. «Демон Лапласа» – уявна істота, яка знає координати, напрям руху та швідкість усіх частинок Всесвіту. За Лапласом, ця істота може розрахувати минулі та майбутні стани Всесвіту на будь-який час. Світ, за Лапласом, є детерміністським, тобто таким, зміни якого повністю визначаються його початковим станом. Важливим результатом, що був отриманий сучасною наукою, є розуміння того, що світ є не детерміністським, а є статистично-ймовірнісним. Особливо яскраво це проявляється у квантовій механіці. З початкових параметрів стану певної системи не можна отримати однозначний прогноз її динаміки; прогноз майбутнього може мати лише вигляд розподілу ймовірностей її майбутніх станів. 

Моделі, які ми будуємо, набагато простіші за Всесвіт (і саме тому їх використання має сенс). Більшість моделей, які ми будували, є детерміністськими. Детерміністські моделі переходять до певного, однозначного стану. Часто люди, знайомі з моделями лише поверхово, очікують, що модель дасть точний прогноз динаміки якогось природного процесу. Найчастіше це неможливо. На перебіг природних процесів в дійсності впливає величезна кількість випадкових подій і непередбачуваних факторів. Природно, в моделі не можна відбити усе різноманіття процесів, що йдуть в дійсності. Однак імітаційна модель може так чи інакше імітувати випадковий хід природних процесів. Такими були наші моделі з ймовірнисним округленням. Результат імовірнісного моделювання буде різним раз від разу. Проаналізувавши розподіл результатів моделювання, ми можемо встановити, до яких станів модель (при заданих початкових параметрах) приходить найчастіше, до яких – рідко, а до яких – не переходить ніколи. Ймовірнісні моделі породжують певний розподіл ймовірностей кінцевих станів. Ця інформація теж буває дуже цінною для розуміння природних процесів, які вивчаються з використанням імітаційного моделювання.

На цьому етапі засвоєння нашого курсу є гарна можливість порівняти ці два типи моделей. Скопіюйте рядки, на яких розміщується система перетворень (див. пункт 1.8) моделі, яку ви побудували на попередньому кроці (рис. 5.1). Вставте копію цих рядків нижче. Якщо все зроблено правильно, блок введення початкових параметрів моделі буде «працювати» на обидва варіанти системи перетворень. Тепер один варіантів системи перетворення, що дублюють один одного, можна змінити, щоб встановити, як впливає ця зміна на роботу моделі.

Залиште верхній варіант системи перетворень детерміністським, а нижній зробіть ймовірнісним. Як ми вже встановили (див. пункт 3.5), це можна зробити, додавши в формули округлення: =ROUNDDOWN("значення_комірки"+RAND();)}. Як ви зрозуміли, «значення_комірки» – це та формула, яка стояла в детерміністській моделі. Природно, розподіл усіх округлення слід додавати тільки в ті формули, які визначають чисельності особин і в яких можуть виникати нецілі значення; наприклад, в формулу для обчислення потреби в ресурсах ніяке округлення додавати не потрібно. Ймовірно, при визначенні чисельності особин першого року (в тому випадку, якщо на попередньому циклі визначено загальну кількість нащадків), досить округлити чисельність представників однієї статі, а потім відняти її з чисельності потомства і отримати чисельність представників протилежної статі.

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

Скоріше за все, результати роботи детерміністської та ймовірнісної моделей будуть відрізнятися. Якщо ми здійснимо перерахунок (хоча б за допомогою кнопки F9), результати роботи ймовірнісної моделі зміняться, а детерміністської – лишаться незмінними. Швидше за все, ми захочемо встановити розподіл результатів імовірнісного моделювання.

Рис. 5.2. Порівняння детерміністської та ймовірнісної моделей «з пташиного польоту» (один з можливих варіантів розміщення моделі на аркуші). При натисканні на цій картинці вона відкриється в новому вікні з більш високою роздільною здатністю; там її можна збільшити та роздивитися як слід. Показано формулу у комірці AS3; ця формула забезпечує запис результатів імітації при певному значенні лічильника. Позначення cnt відповідає показнику лічильника (AJ3), а probab – комірці, в якій міститься підсумок ймовірнісного моделювання (AJ10)

Щоб побудувати розподіл результатів моделювання, використовуємо лічильник ітерацій (ітерації – це повтори, в даному випадку – повторні прогони моделі). Для роботи з лічильником треба включити режим ітерацій. Це можна зробити (рис. 5.3) в меню «Засоби / Параметри / LibreOffice Calc / Обчислити».

Рис. 5.3. Режим ітерацій ввімкнено

Ввімкніть (якщо воно ще не ввімкнене) меню «Перегляд / Панелі інструментів / Елементи керування», робота з яким пояснювалася в пункті 4.1. Включіть «Режим розробки» і вставте в модель елемент управління «Лічильник». Цей елемент управління змінює з кожним клацанням по ньому значення в певній комірці, якою він керує. Надамо комірці, якою керує лічильник, ім'я cnt (counter).

Тепер нам знадобляться комірки, в яких ми будемо запам'ятовувати результати кожної ітерації. Для кожній з цих комірок можна було б вказати, якому значенню лічильника відповідає ця комірка (в одну комірку поміщається значення результату моделювання при показаннях лічильника 1, в іншу – при показаннях лічильника 2 і так далі). Однак найпростіше зробити ряд комірок зі значеннями ряду натуральних чисел (поставити в одну клітинку 1, в іншу – 2 і розтягнути ряд), а в комірки для запам'ятовування результатів моделювання вписати формулу з умовою рівності показань лічильника значенням комірки в цьому ж ряду. На рис. 5.4 показаний такий варіант: ряд з 20 комірок з номерами «обслуговує» 100 комірок для результатів. У першому стовпці комірок для результатів показання лічильника порівнюються зі значеннями комірок для номерів; у другому стовпці до значення комірок з номерами додається 20, в третьому – 40 і так далі.

Результати моделювання відображаються в двох комірках: одній для детерміністської моделі і другій для ймовірнісної. Для характеристики результату моделювання вибрано загальна кількість «статевозрілих» особин ( «дворічок» і «трирічок»). Результати детерміністського моделювання збираються в одній групі комірок, ймовірнісного – в інший. У кожній групі в комірках стоять формули, аналогічні тій, яку видно на рис. 5.2: =IF(cnt=0;"";IF(cnt=AQ4;probab;AR4)). Якщо в комірці лічильника (cnt) стоїть 0, всі комірки для запису результатів звільняються (вираз "" відповідає порожній комірці). Якщо значення лічильника відповідає значенню комірки з номером, в комірку для запису результатів поміщається результат ймовірного моделювання. Якщо значення лічильника якесь інше, значення комірки дорівнює самому собі, тобто не змінюється.

Значення результатів детерміністського моделювання залишаються тими й самими (чого і слід було очікувати). Результати ймовірнісного моделювання змінюються. У моделі на рис. 5.4 використаний наступний прийом. У стовпці AY відображається кількість різних результатів. Для цього використано формулу, що підраховує у певному діапазоні комірок кількість таких, що містять певне значення: =COUNTIF(діапазон;значення). Значення, що шукає ця функція, задані у стовпчику AX. У сьомому рядку цього стовпчика воно відповідає округленому вниз результату детерміністського моделювання, в шостий – на одиницю менше, в сьомий – на одиницю більше тощо. Отриманий розподіл показується на гістограмі.

Рис. 5.4. Лічильник ітерацій (фрагмент моделі, що показана на рис. 5.2). Зверніть увагу на формулу комірки AY9; вона підраховує, скільки разів в блоці з результатами ймовірнісного моделювання зустрічалося значення, що дорівнює числу в комірці AX9

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