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

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

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

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

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

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

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

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

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

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

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

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

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

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

,

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

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

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

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

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

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

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

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

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

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

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

,

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

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

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

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

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

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

,

,

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

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

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

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

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

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

(1.3*)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Вне

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

Введение

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

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

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

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

Опираясь на второй путь, данная работа является расширением моей предыдущей работы «Нелинейная динамическая модель взаимозависимых отраслей производства», а также другой работы (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; Демидович, 1967; Красовский, 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.

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

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

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

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

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

Связь - элемент системной динамики, использующийся для определения зависимости между элементами диаграммы потоков и накопителей.не создает связи автоматически, а принуждает пользователя явно рисовать их в графическом редакторе (однако стоит заметить, что 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). Stability theory of dynamical systems. 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 the 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.). Berlin: 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, and Yorke J. (1987). Chaos, Strange Attractors, and Fractal Basin Boundaries in 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.). John wiley.

15. Llibre J., Valls C. (2007). Global analytic first integrals for the real planar 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). Differential 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). Thirty years of Polynomial System Solving, and now? Journal of Symbolic Computation. 44 (3): 222-231.

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

25. Malthus T.R. (1798). An Essay on the Principle of Population, in Oxford World"s Classics reprint. p 61, end of Chapter VII

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;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Min = i x 1 xmax = r x т;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

| 3 = 1,ш,к = 1,

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

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

"Ькз П " 8=1

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

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

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

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

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

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

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

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

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

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

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

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

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

Птах = max Пг > 0,

Wmax = max Ui < 0;

Umax + Птах <

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

(3)

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

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

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

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

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

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

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

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

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

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

(4)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

(14)

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

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

Таблица 1

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

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

Свойство 1

Свойство 2

Свойство 3

Свойство 4

Свойство 5

Таблица 2

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Заключение

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

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

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

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

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

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

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

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

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

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

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