Моделирование системы "хищник-жертва" с использованием пакета SimPy

Модель Хищник-Жертва и SimPy

Простейшая модель взаимодействия двух популяций — "хищника" и "жертвы" была предложена в начале XX века американским ученым Альфредом Лоткой и является базовой математической моделью взаимодействия двух популяций, позволяющий проследить такие важные явления как периодические колебания, стабилизация и\или неограниченный рост численности.

Математически модель формулируется при помощи дифференциальных уравнений:
$$ \begin{array}{l} \dot{x} = \alpha\cdot x - \gamma_x\cdot x \cdot y \\ \dot{y} = \beta\cdot y + \gamma_y\cdot x\cdot y\end{array},$$

где \(\dot{x}, \dot{y}\) — скорости изменения численности видов "x" и "y" соответственно.

Представленная модель состоит из следующих положений: 1) в отсутствии фактора взаимодействия двух видов ("хищника" и "жертвы"), численность каждого из них удовлетворяет закону экспоненциального роста; 2) определенная часть биомассы (регулируется коэффициентом \(\gamma_x\)) "жертвы" (вид "x") поедается "хищником" (вид "y") исходя из представлений о случайной их встрече: чем больше "хищника" и "жертвы", тем больше биомассы "жертвы" поедается "хищником" (т.е. используется прямая пропорциональность);  3) часть биомассы (регулируется коэффициентом \(\gamma_y\)) съеденной "жертвы"  способствует увеличению численности/биомассы "хищника".

Исследование и решение данной системы не представляет особого труда. В контексте освоения системы для имитационного моделирования SimPy, больший интерес был вызван реализацией данной модели  исходя из ее базовых представлений (положения 1-3) в виде системы конкурирующий процессов: 1) рождения "жертвы"; 2) рождения "хищника"; 3) естественной смерти "жертвы"; 4) естественной смерти "хищника"; 5) поедания "жертвы" "хищником".

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

class Plot(object):
    '''
    env -- среда моделирования
    xinit, yinit -- начальные значения численности жертвы и хищника
    соответственно.
    '''

    def __init__(self, env, xinit, yinit):
        self.x = xinit
        self.y = yinit
        self.env = env

    # сохраняет текущие значения численности жертвы и хищника в глобальных 
    # переменных datax, datay, time
    def print_plot(self):
        global datax, datay, time
        print 'x={0}, y={1}, time={2}'.format(self.x, self.y, self.env.now)
        datax.append(self.x)
        datay.append(self.y)
        time.append(self.env.now)

    # Процесс, ответственный за рождение жертвы
    def bornx(self, alphax):
        while True:
            self.print_plot()
            amount = self.x * alphax
            self.x += amount
            yield self.env.timeout(1)

    # Процесс, ответственный за естественную смерть жертвы
    def deadx(self, betax):
        while True:
            amount = self.x * betax
            self.x -= amount
            yield self.env.timeout(1)

    # Процесс, ответственный за рождение хищника
    def borny(self, alphay):
        while True:
            amount = self.y * alphay
            self.y += amount
            yield self.env.timeout(1)

    # Процесс, ответственный за естественную смерть хищника
    def deady(self, betay):
        while True:
            amount = self.y * betay
            self.y -= amount
            yield self.env.timeout(1)

    # Процесс, ответственный за поедание хищником жертвы
    def yeatx(self, gammax, gammay):
        while True:
            amountx = self.y * self.x * gammax
            amounty = self.y * self.x * gammay
            self.x -= amountx
            self.y += amounty
            yield self.env.timeout(1)

Представитель класса Plot хранит в себе численности жертвы (аттрибут x), численность хищника (аттрибут y) и аттрибут модельной среды env (для доступа к текущему времени моделирования и отсчета шагов при помощи метода timeout).

Методы класса Plot — bornx, borny, deadx, deady, yeatx — отождествляются с процессами рождения жертвы, рождения хищника, естественной смерти жертвы, естественной смерти хищника, поедания хищником жертвы соответственно. Каждый из этих процессов выполняется ровно 1 единицу модельного времени, таким образом, формируемую модель можно рассматривать как дискретный аналог исходной непрерывной модели; правильней, однако, представлять данную модель как самостоятельную, построенную на базе непрерывных представлений. Это связанно с тем, что все перечисленные процессы будут выполняться при моделировании конкурентным образом, и каждый из них будет изменять состояние модели. 

Рассмотрим оставшуюся часть реализации модели:

# ------------- Параметры модели -----------------
alphax = 0.3  # рождаемость жертвы
alphay = 0.04  # рождаемость хищника
betax = 0.005  # естественная смертность жертвы
betay = 0.08  # естественная смертность хищника
gammax = 0.03  # коэффициент интенсивности поедания жертвы хищником
gammay = 0.001  # коэффициент перехода биомассы жертвы в биомассу хищника,
                # за счет поедания
# ------------------------------------------------


# Задаем начальные численности жертвы (30) и хищника (8)
# Единицы измерения здесь не важны, но лучше представлять их
# в данной модели как непрерывные показатели биомассы
p = Plot(env, 30, 8)


# Добавление конкурирующих процессов гибели, рождения, поедания
# в модельную среду
env.process(p.bornx(alphax))
env.process(p.borny(alphay))
env.process(p.deadx(betax))
env.process(p.deady(betay))
env.process(p.yeatx(gammax, gammay))

# выполнение имитационного моделирования до определенного
# времени, в данном случае от 0 до 400
# каждый процесс (гибели, рождения, поедания) занимает 1 единицу
# времени (это описывает строка self.env.timeout(1) в каждом
# процессе)  
env.run(until=400)

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

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

этой ссылке. (4,5 KB)

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

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

Изменение численности хищника и жертвы

Квазипериодический процесс

blog comments powered by Disqus