CC BY
29
Серия «Математика»
2010. Т. 3, № 3. С. 59-66
Онлайн-доступ к журналу: http://isu.ru/izvestia
УДК 518.977, 51-76
Исследование модели динамики популяции методами теории оптимального управления
А. В. Букина
Иркутский государственный университет Ю. С. Букин
Лимнологический институт СО РАН
Аннотация. Рассматривается модель динамики популяции, учитывающая конкуренцию за ресурсы среды обитания между особями с различными адаптивными характеристиками. С помощью теории оптимального управления исследуется зависимость распределения особей по признакам от соотношения интенсивности конкуренции и меры рассеяния ресурсов по фракциям. Для соответствующей интегро-дифференциальной задачи оптимального управления выводится необходимое условие оптимальности и на его основе строится численный алгоритм решения прикладной задачи.
Ключевые слова: динамика популяции; задача оптимального управления интегро-дифференциальной системой; дифференциальный принцип максимума; метод условного градиента; овражный метод.
Пусть х(в,Ь) - плотность особей популяции в момент времени Ь с адаптивной характеристикой в, например, с размером клюва в, который соответствует размеру потребляемого ореха. Плотность ресурсов, доступных для особей с признаком в, задается функцией К (в). Уравнение динамики популяции имеет вид [6, 7]:
Возможные значения адаптивного признака находятся в интервале 5 = [зо,Зі], времени - в интервале Т = [і0,£і]. При этом предполагается, что начальное распределение численности особей популяции известно
1. Введение
у с (s,ox(£,t)d£)
S
(l.l)
и определено условием х(в,Ь0) = х0(в). Особи рождаются с интенсивностью г. Функция С (в, £) определяет конкуренцию между особями за ресурсы среды. Конкуренция предполагает, что ресурс, доступный для особей с признаком в, доступен и для особей с признаком £, только в меньшей мере.
В самом типичном случае функции К (в), С (в,£) являются гауссовскими и имеют следующий вид:
К
(в - 8)21 е (в - О
К(в) = —ехр —’ С(в’^ = —ехр ак\2п 1 2а к 1 ас V 2п
21
2аС
где К - общее количество ресурсов, 8 - наиболее распространенный ресурс, а к - мера рассеяния ресурсов по фракциям, е - количество ресурсов, потребляемое одной особью, ас характеризует убыль конкуренции по мере увеличения модуля разности между в и £.
Соотношение параметров ак и ас определяет характер распределения особей по признакам. И, в частности, может вызывать разделение популяции на два признаковых кластера, что служит основой симпат-рического видообразования (видообразования без географической изоляции). Согласно работам [6, 7] такое разветвление перемешанных особей происходит при ас < а к. Содержательно это примерно означает, что высокая конкуренция в точке наиболее распространенного ресурса не компенсируется соответствующим количеством ресурсов. Использование методов теории оптимального управления позволяет определить зависимость любого распределения особей по признакам и соотношения параметров ак и ас при любом временном периоде.
Задача будет состоять в нахождении параметров ас и а к, при которых в фиксированный момент времени плотность особей популяции максимально близка к какому-либо заданному распределению х(в):
3(а) = J(х(в, ^) — х(в))2йв ^ шш,
(ак, ас) € Q = [ако, аК1] х [асо, аС1]-
2. Задача оптимального управления
Опишем формальную постановку задачи. Будем искать управления, доставляющие минимум функционалу
3(и) = J <р(х(в,Ь1),в)йв + JJ Ф(х,у,и, в,Ь)йвйЬ, (2.1)
5 П
определенному на решениях интегро-дифференциальной задачи вида
Xt = f (x, y, u, s, t), y(s, t) = j g(x(€, t),u((, t),s, £, t)d£,
S
x(s,t0) = x0(s). (2.2)
При этом допустимые управлениях u = u(s, t) принадлежат множе-
ству
V = {u € L^(n) : u(s, t) € U}. (2.3)
Здесь (s,t) € П = S x T = (s0,si) x (t0,ti), x(s,t) € Rn, y(s,t) € Rm,
U c Rr.
На параметры задачи (2.1)—(2.3) наложим следующие ограничения: функции ф = ф(х, s), Ф = Ф(х, y, u, s, t), вектор-функции f = f (x,y,u, s,t), g = g(x,u, s,£,t) непрерывны по совокупности своих аргументов и имеют непрерывные производные по переменным x и y, x0(s) - измеримая существенно ограниченная в S вектор-функция. Также будем считать, что вектор-функции f и g удовлетворяют условию Липшица по переменным x, y и x в областях Rn x Rm и Rn соответственно.
Для интегральной формы задачи (2.2) классическим методом последовательных приближений можно доказать существование и единственность решения в пространстве L(П) для любого допустимого управления. А также получить оценки приращений решений, отвечающие двум допустимым управлениям u и u:
НДу^-ОН < (\\Айд(х,и,8,С,1)\\ + ||Дж(£,г)||)^ (2.4)
где константа М > 0 не зависит от входных данных х0, д и f.
Путем анализа формулы приращения целевого функционала на двух допустимых процессах (и, х,у) и (и = и + Ди, X = х + Дх, у = у + Ду) можно получить [3] следующее представление
Д7(и) = — JJ ДйН(ф,в,х,у,и,з,1)й8в;Ь + п(Ди), (2.5)
п
где Н(ф, в, х, у, и, в, £) - аналог функции Понтрягина вида
Н(ф,в,х,у,и,в,і) = і), /(х,у,и,в,і)) +
+ / {в(£,і),д(х,и,£, в,і)) — Ф(х,у,и, в, і),
я
ф,в - решения сопряженной задачи
фі = —Нх(ф, в, х, у, и, в, і), ф(в, Іі) = —ірх(х(в, І1), в),
в(в,і) = Ну(ф,в,х,у,и, в, і), (в, і) Є П, (2.6)
П(Ди) - остаток формулы приращения вида
я
— Л ({ДиНх, Дх) + {ДиНу, Ду) + он (|| (Дх, Ду)||) (ів(1і.
п
В рассматриваемой прикладной задаче все компоненты управления являются параметрами, функции /,д, Ф дифференцируемы по и, множество и выпукло. Поэтому рассмотрим формулу приращения функционала (2.5) на классической вариации управления, линеаризовав ее по и. С учетом оценок (2.4) получим дифференциальный принцип максимума.
Пусть процесс (и*,х*,у*) оптимален, ф*,в* - соответствующее ему решение сопряженной задачи (2.6). Тогда почти в каждой точке (в,Ь) € П справедливо
Данное необходимое условие оптимальности служит основой для построения градиентных методов численного поиска локальных решений. Выпишем основные конструкции метода условного градиента.
п
п
Н(ф, в, х, у, а, в, і)
я
У в(і,і)ехр
фг = —фт( 1 - у
ак\[2к (в — в)2'
ехр
<7С^/~2П .
к
в(е, Ь)ехр
27к
(е — в)2
272
ф(в, ^1) = —2(х(в,£1) — х(в)), в(в, ^) = —фгх7К^2Пехр ——2^• (2.7)
к
27к
Пусть задано начальное управление 70 = (7°к, 7°) € Q, положим к = 0. Из интегро-дифференциальных задач (2.2), (2.7) построим соответствующие значения фазовых и сопряженных траекторий хк, ук, фк, вк. Определим вспомогательное управление 7к как решение задачи
У3(7) = —JJ На(ф,в,х,у,7,в,1)йв<И. п
Подсчитаем величину Ш(7к) = (7к),7к — 7к^ < 0. Если Ш(7к) = 0,
то управление 7к удовлетворяет дифференциальному принципу максимума. В случае Ш(7к) < 0 зададим а-параметрическое семейство управлений 7к(а) = 7к + а(7к — 7к), а € [0,1], положим 7к+1 = 7к(ак). Существует несколько способов задания положительной величины ак [1]. Так как параметры задачи вместе с производными по переменным состояния и управления удовлетворяют условию Липшица по этим переменным, целевой функционал ограничен снизу, данный алгоритм является релаксационным и сходится по невязке дифференциального принципа максимума.
Реализация описанного алгоритма и серия экспериментов при разных начальных управлениях показали, что целевая функция имеет овражный характер. То есть существует некая зависимость параметров 7с и 7к , отклонение от которой приводит к резкому увеличению значений функционала. Для эффективного поиска точек на «дне оврага» был
применен один из вариантов «овражного метода» [1]. Пусть имеются к_____________1 к
две точки 7к 1, 7к на дне оврага, полученные несколькими итерациями метода условного градиента. Из точки
7к+1 = 7к — (7к — 7к_1)\\7к — 7к_1 Ц-1 Н в1дп(<1 (7к) — 3(7к_1))
совершается спуск градиентным методом и находится следующая точка на дне оврага 7к+1. Величина овражного шага Н подбирается эмпирически в соответствии с кривизной дна.
3. Результаты
Для установления зависимости точек, полученных овражным методом, применялся метод наименьших квадратов. Достоверной оказалась аппроксимация линейной функцией. Применительно к исследуемой модели данный результат означает, что всякое распределение численности особей по признакам определяется соответствующей линейной зависимостью параметров распределения ресурсов и интенсивности конкуренции.
Рассмотренный подход был опробован с использованием следующих значений параметров модели
Т = [0,100], 5 = [0.5,1.5], х(в, 0) = 5, г = 1, К = 100, 8 = 1,
Q = [0.2,1] х [0.2,1],
К
х(в) =
2 * 0.2 * V2п
г (в — в')2 г (в — в")2,
ехР[— ъ . п о2 ] + ехР[— ъ . п о2 ]
2 * 0.22
2 0.22
где значения в', в'' выбирались равномерно удаляющимися от 1 до концов интервала Б. Такое задание х(в) основывается на результатах, опубликованных в работах [2]—[5].
Приведем полученные соотношения 7с и 7к для четырех случаев:
1) для в'=1, в''=1 (Рис. 1а) получили зависимость 7к = 0.837с + 0.11,
2) в'=0.7, в''=1.3 (рис. 1Ь): 7к = 0.887с + 0.15,
3) в'=0.6, в''=1.4 (рис. 1с): 7К = 0.927С + 0.145,
4) в'=0.5, в''=1.5 (рис. 1ё): 7К = 0.97С + 0.17.
с й
Рис. 1. Распределение численности особей по признакам. -ление х(в),-----полученное распределение х(8,Ь{).
заданное распреде-
Если полученное распределение на рис. 1с охарактеризовать уже как формирование двух видов, то можно говорить о том, что при парамет-
X
Б
Ь
X
200
150
10(
50
рах 7с и 7к из Q, удовлетворяющих соотношению 7к > 0.927с + 0.145, происходит видообразование.
Неточное соответствие полученного распределения х(в,^) фактическому х(в), особенно в случае Ь, означает, что модель динамики популяции (1.1) нуждается в корректировке.
Список литературы
1. Васильев Ф. П. Численные методы решения экстремальных задач / Ф. П. Васильев. - М. : Наука, 1980. - 520 с.
2. Ситникова Т. Я. О глубоководных «карликах» и «гигантах» среди байкальских эндемичных гастропод / Т. Я. Ситникова, М. Н. Шимараев // Журн. общ. биологии. - 2001. - T. 62, вып. 3. - С. 226-238.
3. Терлецкий В. А. Вариационный принцип максимума в задаче оптимального управления интегро-дифференциальной системой / В. А. Терлецкий, А. В. Букина // Вестн. Бурят. ун-та. Сер. 13, Математика и информатика. - 2008. -Вып. 9. - С. 52-55.
4. Тихонова Е. Н. Морфологический анализ байкальских амфипод Pallasea Cancellus из реки Ангары / Е. Н. Тихонова, Р. М. Камалтынов // Бюл. ВСНЦ СО РАМН. - 2007. - Вып. 1. - С. 108-112.
5. Широкая А. А. Распределение моллюсков семейства ACROLOXIDAE (GASTROPODA, PULMONATA) в озере Байкал / А. А. Широкая, Н. В. Максимова, Т. Я. Ситникова // Зоол. журн. - 2008. - T. 87, вып. 5. - С. 532-546.
6. Dieckmann U. On the origin of species by sympatric speciation / U. Dieckmann, M. Doebeli // Nature. - 1999. - N 400. - P. 354-357.
7. Doebeli M. Evolutionary branching and sympatric speciation caused by different types of ecological interactions / M. Doebeli, U. Dieckmann // Am. Nat. - 2006. -Vol. 156. - P. 77-101.
A. V. Bukina, Yu. S. Bukin
Investigation of population dynamics model by the methods of optimal control theory
Abstract. We consider population dynamics model that takes into account resources competition between individuals with different adaptive characteristics. With the help of optimal control theory the dependence of distribution of individuals by trait on the correlation of competition intensity and the measure of resource distribution is investigated. For corresponding integro-differential optimal control problem we establish the necessary optimality condition and on its basis we construct a numerical solution method for the applied problem.
Keywords: population dynamics; optimal control of integro-differential system; differential principle of maximum; conditional gradient method; ravine method
Букина Анна Викторовна, младший научный сотрудник, Иркутский государственный университет, 664003, Иркутск, ул. К. Маркса, 1 тел.: 8-964-5-43-58-78 (annabukina@mai. ru)
Букин Юрий Сергеевич, кандидат биологических наук, старший научный сотрудник лаборатории геносистематики, Лимнологический институт СО РАН, 664033, Иркутск, 33, А/я 278, ул. Улан-Баторская, 3, тел. (3952) 42-29-23, ([email protected])
Bukina Anna, Irkutsk State University, 1, K. Marks St., Irkutsk, 664003, junior research scientist, Phone: 8-964-5-43-58-78, ([email protected]) Bukin Yury, Limnological Institute of SB RAS, 3 Ulan-Batorskaya Str. P.O. Box 278, Irkutsk, 664033, senior research scientist, Phone: (3952) 4229-23, ([email protected])
