CC BY
127
2010 Дискретные модели реальных процессов №4(10)
УДК 621.391.1:004.7
МЕТОД ПОСТРОЕНИЯ КЛЕТОЧНО-АВТОМАТНЫХ МОДЕЛЕЙ ПРОЦЕССОВ ФОРМИРОВАНИЯ УСТОЙЧИВЫХ СТРУКТУР1
О. Л. Бандман
Институт вычислительной математики и математической геофизики СО РАН,
г. Новосибирск, Россия
E-mail: bandman@ssd.sscc.ru, bandman@academ.org
Представлен метод построения клеточных автоматов (КА), моделирующих процессы самоорганизации при формирования устойчивых структур. Метод основан на параллельной композиции двух КА. Рассмотрено два случая: 1) когда один компонентный КА функционирует независимо, влияя на эволюцию второго, и 2) когда оба КА взаимодействуют на каждой итерации. Метод иллюстрируется результатами компьютерного моделирования двух самоорганизующихся систем: формирования устойчивых пространственных структур на подогретой поверхности и достижения баланса между хищником и жертвой.
Ключевые слова: математическое моделирование, клеточный автомат, диссипативные структуры, самоорганизация.
Введение
Интерес к математическому и компьютерному моделированию диссипативных процессов в неравновесных системах вызван, прежде всего, интенсификацией исследований химических, биологических, социологических и других явлений. Среди процессов с такими свойствами известны обычные, спиралевидные, концентрические, вращающиеся и стоячие волны (автоволны), а также колебательные и устойчивые пространственные структуры. Наиболее известны теоретические исследования диссипативных процессов [1], химических автоволн [2], биологических и экологических явлений [3]. Изучение устойчивых диссипативных структур началось с появления работы А. Тьюринга по морфологии [4]. Полный и хорошо написанный обзор исследований в этой области дан в [5]. Почти все известные работы основаны на математических моделях традиционного типа — системах дифференциальных уравнениях с частными производными. Поскольку исследуемые явления диссипативны и пространственно распределены, они относятся к процессам типа «реакция — диффузия», описываемым уравнениями с нелинейными членами, которые не имеют аналитических решений, а численные методы требует больших вычислительных ресурсов.
Однако, поскольку научный интерес к диссипативным структурам растет, поиск новых моделей интенсифицируется. В частности, появляются модели дискретного и стохастического типа, в которых химические преобразования и пространственная подвижность частиц вещества моделируются непосредственным (прямым) образом. К таким моделям относятся клеточные автоматы [6], клеточно-нейронные сети [7] и кинетические методы Монте-Карло [8]. В данной работе используется математическая модель, основанная на расширенном понятии клеточного автомата [9], отличающаяся
хРабота поддержана Программой фундаментальных исследований Президиума РАН №2-6 (2010) и Сибирским отделением РАН, Интеграционный проект 32 (2010).
от классического КА фон Неймана тем, что в ней допускается применение любого алфавита, любых функций переходов, любых режимов смены состояний клеток. Такая широкая трактовка КА позволяет строить математические описания пространственновременных процессов разного характера, в том числе и обладающих свойствами самоорганизации и формирования устойчивых образов (структур).
Использование КА-моделей для компьютерного моделирования пространственной динамики в химии и физике привлекательно по следующим причинам.
• КА моделирует дискретные изменения состояний в дискретных пространстве и времени, что позволяет непосредственным образом отображать в модели перемещения и трансформации реальных частиц или агентов.
• На функции переходов КА не накладывается никаких ограничений: они могут быть нелинейными, разрывными и вероятностными, что позволяет моделировать такие процессы, как фазовые переходы и химические превращения, а также учитывать условия, при которых допустимо изменение того или иного состояния.
• КА допускают моделирование не синхронизированных процессов (асинхронные КА), что соответствует естественному течению событий во всех реальных явлениях, в которых не введена искусственно синхронизация.
• КА моделируют объекты, которые в [6, 10] называются сложными системами (complex systems), так как имеют очень простые математические представления, но моделируют сложные пространственно-временные процессы. Отсюда вытекает простота программирования КА-моделей как в последовательном, так и в параллельном вариантах.
Исследование КА как дискретных моделей пространственной динамики началось в 80-х годах прошлого столетия работами [11, 12] и в настоящее время интенсивно развивается, охватывая все новые области задач и совершенствуя теорию и методологию [10, 13, 14]. К развитию методологии КА-моделирования и относится предлагаемая работа. В ней после введения даны формальные определения и пример простого КА (п. 1), моделирующего процесс образования устойчивых структур. В п. 2 представлен метод параллельной композиции КА [15], в п. 3 и 4 иллюстрируются два типовых случая его применения для построения КА-моделей этого класса.
1. Простые самоорганизующиеся клеточные автоматы
Атомарным понятием КА является клетка, которая характеризуется именем x Е X и состоянием а Е A, где X — множество имен клеток, которое обычно интерпретируется как множество координат точек в дискретном пространстве конечных размеров; A — алфавит состояний клеток. Состояние клетки x в момент времени t обозначается как vx(t). На множестве имен определен шаблон соседства T(x) = [pj(x) : j = = 0,... , q}, который определяет имена клеток, взаимодействующих с клеткой x. Если пространство моделирования представлено решеткой с прямоугольной, треугольной или квадратной формой клетки, то функции ipj (x) Е T(x) имеют вид сдвиговых: Pj (x) = x + lj, а шаблон состоит из q плотно примыкающих друг к другу клеток с условным «центром» в x. Шаблон связывает с каждой клеткой x Е X локальную конфигурацию V(x) = [vi,...,vq} таким образом, что каждая клетка (x) Е T(x) имеет состояние Vj Е V(x). Множество состояний всех клеток x Е X в момент t называется глобальной конфигурацией Q(t) = [vx(t) : x Е X}.
Функционирование КА задается локальным оператором над подстановками 0(01(x),... , 0ra(x)). Если локальный оператор задан одной подстановкой, то он назы-
вается простым. Подстановка имеет следующий вид:
0(x) : V(x) U V''(x) ^ V'(x), (1)
где V(x), V'(x) и V''(x) — базовая, новая и контекстная локальные конфигурации соответственно, причем V(x) и V'(x) связаны с одним и тем же шаблоном T(x). Подстановка называется вероятностной, если ее применение выполняется при условии rand < p, где rand — случайное число в интервале [0,1], p — заданная вероятность. Применение в к клетке x состоит в замене состояний vj Е V(x) на новые vj Е V'(x), которые являются значениями функций перехода вида
vj = fj(V(x) U V''(x)). (2)
Применение подстановки к клетке x на итерации t называется успешным, если
(V(x) U V''(x) С Q(t)) & (rand < p), (3)
даже если не произошло ни одного изменения состояний, т. е. fj (V(x)) = V(x). Если (3) не выполняется, то попытка применения подстановки безуспешна.
Функционирование КА подчиняется итерационному алгоритму, в котором итерацией считается применение локального оператора (успешное и безуспешное) ко всем клеткам x Е X. При этом происходит переход Q(t) ^ Q(t + 1). Последовательность П(0),... , Q(t),... , Q(i) называется эволюцией КА. Если после какой-то итерации t = i дальнейшее успешное применение локального оператора не изменяет глобальную конфигурацию или меняет ее не более чем в заданном числе клеток е, т. е.
|X| — |Q(T +1) П ft(T)| <е, (4)
то КА обладает свойством самоорганизации и относится к классу 2 по классификации Вольфрама [6], ас точки зрения качественной теории нелинейных систем эти КА характеризуются тем, что имеют устойчивый фокус. Свойство самоорганизации проявляется тогда, когда система состоит из множества элементов, взаимодействующих как с положительной, так и с отрицательной обратной связью (физическая интерпретация). Иными словами, взаимодействующие элементы являются по отношению друг к другу либо активаторами, либо ингибиторами (химическая интерпретация).
Порядок применений подстановок при переходе от Q(t) к Q(t + 1) называется режимом. Базовыми являются два режима: синхронный и асинхронный. При синхронном режиме итерация состоит из двух следующих шагов:
1) вычисляются значения функции переходов (2) для всех клеток, в которых (3) выполнено;
2) в базовых конфигурациях всех клеток производится смена состояний согласно (1).
Для синхронного режима важно, чтобы выполнялось условие корректности вычислений [13]
Vt = 1,... , £ (T(xj) П T(xj) = 0).
Очевидно, что оно всегда выполнено при | V'(x)| = 1.
При асинхронном режиме итерация состоит из последовательности |X| шагов. На каждом шаге
1) случайно выбирается имя клетки xj Е X;
2) вычисляются значения функций переходов для новой локальной конфигурации У'(ач) по (2);
3) в базовой конфигурации V(хг) производится смена состояний на соответствующие из V'(хг).
Среди простых самоорганизующихся КА наиболее известны два типа: 1) КА, моделирующие процессы, называемые «разделением фаз», и 2) КА с взвешенными шаблонами, моделирующие процессы формирования устойчивых структур [7]. Оба эти типа КА имеют устойчивые состояния как при синхронном, так и при асинхронном режимах. КА первого типа подробно описаны в [13, 14], их эволюция обладает свойством самоорганизации, но получающиеся в результате устойчивые структуры тривиальны — в них •цг(£) = 1 или •цг(£) = 0 для всех х Е X. Второй тип КА принадлежит к классу клеточно-нейронных сетей [7]. Использование в их функциях перехода взвешенных шаблонов позволяет получать устойчивые структуры в виде разного вида пятен и полос.
Пример 1 (рис. 1). Типичный представитель класса КА формирования устойчивых структур имеет А = {0,1}, X = {(г,') : г,' = 0,1,... , М}, а локальный оператор
опирается на шаблон
т (г') = {(г + ^ : = -3,..., ^..^ 3}. (5)
Этому шаблону ставится в соответствие весовая матрица Ш, элементы которой равны
ч Г а, если Ы ^ ^ ^ ^ 1, (6)
ф \ Ь иначе, (6)
где а > 0 соответствует активатору, а Ь < 0 — ингибитору. Локальный оператор 0(г,') меняет состояние только одной клетки, используя в функции переходов состояния всех клеток из т(г, '):
^(г') : {^+М+ : = ^..^ 0,..., 3} ^ , (7)
где
#=0 #=10 #=25 t=50
Рис. 1. Эволюция простого КА, моделирующего процесс формирования устойчивых структур типа «полоски»
Эволюции этого КА при работе в синхронном и асинхронном режимах имеют один и тот же поведенческий характер, причем асинхронный режим приводит к устойчивому состоянию за меньшее число итераций.
2. Параллельная композиция КА
Параллельная композиция [15] двух КА N1 = (А1,Х1, @1) и Н2 = (А2,Х2, @2) предполагает, что каждый КА функционирует на своем собственном клеточном поле, имеет свой собственный алфавит состояний, а между множествами их имен Х1 и Х2 существует взаимно-однозначное соответствие.
Композицию могут составлять только КА с одинаковыми режимами работы: либо оба синхронные, либо оба асинхронные. В синхронном случае смена состояний автоматов происходит в обоих КА одновременно, обозначая переход к следующей итерации. В асихронной композиции клетки каждого КА меняют свои состояния независимо от своего окружения, что дает возможность работать компонентным КА с произвольными скоростями. Однако «разноскоростные» композиции пока не исследованы. Принято считать, что на одном отрезке времени в каждом КА происходит одинаковое количество применений локальных операторов.
В [15] введено два типа параллельной композиции КА: однонаправленная и двунаправленная.
3. Однонаправленная параллельная композиция
В однонаправленной композиции один из КА, пусть Н1, является основным. Он имитирует исследуемый процесс. Второй КА Н2 играет роль контекста — он эволюционирует автономно, влияя на работу первого. В функции переходов локального оператора ©2 аргументами являются состояния клеток из П2, а в функции переходов @1 — состояния клеток из П2 и П1. Соответственно шаблоны локальных конфигураций принадлежат разным множествам имен:
Ух е Х1 (Т1(х) С Х2 и Х1), Ух е Х2 (Т2(х) С Х2). (9)
Однонаправленная композиция двух КА имеет эволюцию, формирующую устойчивые структуры, если этим свойством обладает основной КА. Контекстный КА играет роль управления, которое может динамически изменять параметры основного как во времени, так и в пространстве. Это дает возможность существенно увеличить разнообразие получаемых структур.
Пример 2. Примером однонаправленной асинхронной композиции может служить КА-модель коагуляции некоторого вещества, наличие которого в точке, соответствующей клетке х, имитируется состоянием их = 1, а отсутствие — состоянием их = 0. На процесс воздействует тепловое поле, конфигурация которого представляет собой холодный квадрат в центральной части и сильно нагретую область вокруг него (рис. 2). Температура в каждой точке имитируется осредненным значением состояния в соответствующей клетке. Коагуляцию моделирует асинхронный КА Ни = (Аи,Хи, @и), а изменения теплового поля — асинхронный КА диффузии N = (Аи,Хи, @и). При этом Аи = А^ = {0,1}, Х^ = Хи = {(г,') : г,' = 0,..., 300}.
КА коагуляции Ни принадлежит к классу асинхронных КА, моделирующих процесс формирования устойчивых образов, его локальный оператор такой же, как (7), но отличается значениями ингибиторов Ь, которые равны осредненным значениям состояний в соответствующих клетках П (¿), т. е.
Ь = И = (2г + 1)-1 Е vг+g,j+h, (10)
д,Н=-т
где г = 10 — радиус осреднения. N — асинхронный КА, моделирующий так называемую наивную диффузию [16]. Этот КА работает независимо от Ни, имитируя распро-
Рис. 2. Начальная глобальная конфигурация контекстного КА Ov(0)
странение тепла от горячей области (V = 1) к холодному квадрату в центре (V = 0) и к теплой области по краям, где плотность распределения единиц равна (V) =0,5 (рис. 3).
t=2
=20
t= 150
Рис. 3. Эволюция основного КА вида (7), в котором значения ингибиторов управляются распространяющимся теплом, моделируемым контекстным КА Hv с начальным значением Ov (0), показанным на рис. 2
Локальный оператор КА Nv равен 0v = {#v}, где
$v : {vi,j , v^1(i,j), v^2 (i,j), v^3 (i,j) , v^4 (i,j)} ^ {Vi,j , vtp 1 (i,j) , v<^2 (i,j), v<^3 (i,j) , v<^4 (i,j) }
с функцией перехода
v0 = v^fe(i,j), если 0,25k < rand < 0,25(k + 1);
v0, если 0,25k < rand < 0,25(k + 1),
Vu =
(12)
VVk (*,j)
иначе,
k = 1,..., 3.
В процессе эволюции контекстного КА плотность состояний в нем становится равномерной (температура выравнивается) и система приходит в устойчивое состояние, в котором основной массив имеет вид горизонтальных и вертикальных полосок (рис. 3).
4. Двунаправленная параллельная композиция
В двунаправленной композиции ^ и работают в тесном взаимодействии, причем каждый играет роль активатора для себя и ингибитора для другого, но ни один из них сам по себе не должен быть самоорганизующимся КА. В функциях перехода их
локальных операторов @1 и @2 контекстными локальными конфигурациями являются состояния клеток из П2 и П1. Соответственно для шаблонов локальных конфигураций выполняется
Ух е Х1 (Т1(х) С Х1 и Х2), Ух Е Х2 (Тг(х) С Х2 и Х1). (13)
Пример 3. Типичным примером двунаправленной композиции может служить задача установления баланса в системе типа «хищник — жертва». Задача интерпретируется следующим образом. На площади, представленной дискретным пространством, обитают два типа живых существ: хищник и жертва. Хищник питается жертвой. Если для хищника имеется достаточно пищи, то его плотность увеличивается с вероятностью, пропорциональной количеству сытых хищников (только сытые хищники размножаются). Если пищи не хватает, то плотность хищников уменьшается. Жертва всегда стремится размножиться. Обе популяции перемещаются (диффундируют), причем хищник более подвижен, поэтому характеризуется коэффициентом диффузии, много большим, чем жертва (^/^и = К, К > 10). Пусть N = (А^, Х^, ©^) моделирует поведение хищника, а Ни = (Аи,Хи, ©и) —жертвы, причем состояния клеток (г,'% Е Х^ обозначаются символами уга-, а состояния клеток (г,')и Е Хи — символами . Локальные операторы содержат по две подстановки: ©^ = {0^,0^}, ©и = {6и1, 6и2}. В обоих случаях первые подстановки 6^1 и 6и1 моделируют диффузию, используя КА из примера 1 с локальным оператором вида (11). Разница в коэффициентах диффузии реализуется путем разного количества применений 6^1 и 6и1 в течение одной итерации. Вторые подстановки моделируют взаимоотношения хищника и жертвы, которые зависят от состояний локальных конфигураций в обоих клеточных массивах:
: К,-} и (К"(г>') и К/(г')) ^ К.?},
6и2 : {иг,^'} и ' ) и ^ {<^},
где контексты ^^,/(г,') и УйТ(г,') содержат состояния клеток из
') = {(г + 0'' + ^ = -10^ .. ,-1,1,..., 10}
а новые состояния являются значениями вероятностных функций переходов
Г14)
u
0, если (ui,j) > (vi,j) & rand <pv^0,
1, если (vi,j) > (ui,j) & rand <pv^!,
0, если (ui,j) < (vi,j) & rand <pu^0,
1, если (ui,j) > (vi,j) & rand <ри^ь
f15)
где (у) —осредненное по (10) значение у.
Значения вероятностей в (15) определяются исходя из следующих соображений. Если (уга-) > (ига-), то плотность хищников уменьшится на (уга-) — (ига-) из-за нехватки пищи. Эта разность равна вероятности того, что у' ^ = 0. Если же (уга-) < (иг,^), то пищи достаточно всем, и хищник увеличивает свою плотность в соответствии с функцией размножения ^((уг,^)) = с■ (уг,^)(1 —(уг,^)) [17]. Аналогично, если (у^-) > (и^), то жертва поедается хищником, ее плотность уменьшается с вероятностью, равной плотности хищника. В противоположном случае остаток жертвы размножается с вероятностью,
равной (мг-) — (гг,-). Таким образом,
Р^0 = (^г-) — («г,-))/(^г,і), если (г^) > (Мі,-),
р^і = 0,5(гг,-)(1 — (гі,-}), если (гі,-) < (мг,-), (16)
Р«^о = К-), если (г*-) > (и*,-),
р«^1 = (Мг,і) — (г*,-), если (г*,-) < (мг,- ).
Композиция Н и Над функционирует в соответствии со следующей итерационной процедурой:
1) подстановка 0^1 применяется К раз, всякий раз к случайно выбранной клетке из ;
2) подстановка 0ад1 применяется К х раз, всякий раз к случайно выбранной клетке из Хад;
3) случайно выбирается клетка из , к ней применяется 6^2;
4) случайно выбирается клетка из Хи, к ней применяется 0и2.
Эти действия повторяются до тех пор, пока не достигается устойчивое состояние. На рис. 4 приведены три глобальных состояния в эволюции Н. В исходном состоянии (і = 0) плотность жертвы на всем пространстве равна = 0,5, а хищник плотно сосредоточен в средней полосе (гг-) = 1.
t=0 t=4
ф* ¿г ■«
г ж ■ тЛ
фШ 4
t=15 t=30
Рис. 4. Три глобальных состояния в эволюции КА ^v. При t > 30 изменений Qv и не наблюдается
Моделирование выполнялось на компьютере Intel® Core™ i7. Результаты показали, что представленная выше система «жертва — хищник» быстро самоорганизуется и очень устойчива. Все начальные конфигурации, имеющие в каждом КА даже самые малые плотности и по-разному распределенные, эволюционировали к одному и тому же устойчивому состоянию (рис. 4, t = 30).
Заключение
Предложен метод моделирования нелинейных процессов типа «реакция — диффузия» при помощи параллельной композиции клеточных автоматов. Рассмотрено два случая: однонаправленная параллельная композиция, которая позволяет ввести динамическое управление моделируемым процессом, и двунаправленная параллельная композиция, которая моделирует процессы формирования самоорганизующихся структур. Результаты позволяют надеяться, что предложенный подход будет полезен в создании систематической методологии синтеза самоорганизующихся КА.
ЛИТЕРАТУРА
1. Nicolis G., Prigogine I. Self-Organization in Nonequilibrium Systems. N.Y.: Wiley-Interscience, 1977. 236 p.
2. Жаботинский А. М. Периодические процессы окисления малоновой кислоты в растворе (исследование кинетики реакции Белоусова) // Биофизика. 1964. Т. 9. C. 306-310.
3. Deutsch A., Dormann S. Cellular Automaton Modeling of Biological Pattern Formarion. Berlin: Birkhauser, 2004. 330 p.
4. Turing A.M. The chemical basis of Morphogenesis // Phil. Trans. R. Soc. London. 1952. V. B 237. No. 641. P. 37-82.
5. Ванаг В. К. Диссипативные структуры в реакционно-диффузионных системах. М.; Ижевск: Ин-т компьютерных исследований, 2008. 300 c.
6. Wolfram S. A new kind of science. USA: Wolfram Media Inc., 2002. 1197p.
7. Chua L. CNN: A paradigm of complexity. Singapore: World Scientific, 2002. 320 p.
8. Elokhin V. I., Latkin E. I., Matveev A. V., Gorodetskii V. V. Application of Statistical Lattice Models to the Analysis of Oscillatory and Autowave Processes on the Reaction of Carbon Monoxide Oxidation over Platinum and Palladium Surfaces // Kinet. Catalys. 2003. V. 4. No. 5. P. 672-700.
9. Achasova S., Bandman O., Markova V., Piskunov S. Parallel Substitution Algorithm. Theory and Application. Singapore: World Scientific, 1994. 180 p.
10. Simulating Complex Systems by Cellular Automata / eds. A. G. Hoekstra, J. Kroc, P.M.A. Sloot. Berlin: Springer, 2010. 350p.
11. Toffolli T. Cellular Automata as an Alternative to (rather than Approximation of) Differential Equations in Modeling Physics // Physica D. 1984. V. 10. P. 117-127.
12. Wolfram S. Statistical mechanics of Cellular automata // Rev. Mod. Phys. 1993. V. 55. P. 607-640.
13. Бандман О. Л. Клеточно-автоматные модели пространственной динамики // Системная информатика: Сб. научн. тр. Новосибирск: Изд-во СО РАН, 2006. Вып. 10. C. 59-111.
14. Бандман О. Л. Дискретное моделирование физико-химических процессов // Прикладная дискретная математика. 2009. №3(5). C. 33-49.
15. Бандман О. Л. Методы композиции клеточных автоматов для моделирования пространственной динамики // Вестник Томского госуниверситета. Приложение. 2002. №9(1). С.188-192.
16. Toffolli T., Margolus N. Cellular Automata Machine. USA: MIT Press, 1987. 280 p.
17. Свирежев Ю. М. Нелинейные волны, диссипативные структуры и катастрофы в экологии. М.: Наука, 1987. 320 c.
