Sugarscape Model («сахарная модель»)
Имеется некоторая окружающая среда — двумерный ландшафт, где разбросан сахар. Туда же помещены агенты-жуки, которые ползают по сахарному ландшафту по простым правилам: агенту нужно есть сахар, и он перемещается туда, где сахара больше.
- Для каждого агента указывается количество сахара, которое он должен потребить в единицу времени (метаболизм).
- У каждого агента имеется запас сахара, который он может переносить с собой из клетки в клетку (переменная величина).
- Агент может «видеть» и соответственно перемещаться в соседние клетки. Он осматривает видимые им клетки и выбирает пустую с наибольшим количеством сахара. Если таковых несколько, то выбирается ближайшая.
- Агент передвигается в нее и получает весь находящийся там сахар. Таким образом, количество сахара у агента складывается из того, что он имел, и полученного в новой клетке за вычетом потребления в данном временном интервале.
- Агент погибает, если ему нечего есть.
Цикл работы «Сахарного мира»:
Если в клетке есть агент, то
Запомнить текущие расположение и количество сахара
Найти среди соседних ячеек ячейку с наибольшим количеством сахара
Если количество сахара в новой ячейке больше, чем в текущей, то перейти в новую ячейку
Конец «если»
Программа
function sugarscape
% Модель "Сахарный мир" (Sugarscape, by Epstein J. and Axtell R.
% "Growing Artifical Societies: Social Scienses from the Bottom Up"
% Washington DC, 1996
% Правило для ландшафта: Ginf (с.26); правило для агентов: M (с.25).
% Параметры модели:
nstep = 20; % число шагов
size = 50; % size x size - размер мира (четное число)
maxmetabolism = 4; % макс. метаболизм (потребление сахара за шаг)
maxvision = 6; % макс. число клеток, на которые агент способен "видеть" (< size)
maxsugar = 20; % макс. объем сахара в клетке
% Задаем ландшафт
s = initSugarScape(size, maxsugar);
% Задаем расположение и параметры агентов
a_str = initAgents(size, s, maxvision, maxmetabolism);
% Цикл моделирования
for step = 1:nstep;
% Отобразим ландшафт и агентов
displayScape(a_str, s, size)
% Применим правила к каждому агенту
for i = 1:size
for j = 1:size
% i = randperm(size) агентов можно выбирать случайным образом
if a_str(i,j).active == 1 % В данной клетке есть агент?
% Запомним положение агента и количество сахара в текущей клетке
temps = s(i,j);
tempi = i;
tempj = j;
% Ищем наилучшую ячейку среди тех, что агент может "видеть"
for k = a_str(i,j).vision : -1 : 1
[temps,tempi,tempj] = see(i,j,k,a_str,s,size,temps,tempi,tempj);
end
% Агент движется в лучшую ячейку, обновляя свой запас сахара
a_str = moveAgent(a_str, s, i, j, temps, tempi, tempj);
end % if
end % for j
end % for i
end % for step
function s = initSugarScape(size, maxsugar)
% Задание распределения сахара с двумя пиками
% s1 - распределение сахара с северо-восточным пиком
x = -ceil(0.75*size) : size-ceil(0.75*size)-1;
y = -ceil(0.25*size) : size-ceil(0.25*size)-1;
[X,Y] = meshgrid(x,y);
s1 = maxsugar./(abs(X) + abs(Y));
s1(X == 0 & Y == 0) = maxsugar;
% s2 - распределение сахара с юго-западным пиком
s2 = s1';
% Распределение с двумя пиками
s = s1 + s2;
function a_str = initAgents(size, s, maxvision, maxmetabolism)
% Задание расположения и параметров агентов
for i = 1:size;
for j = 1:size;
if rand < 0.2 % поместить в клетку агента и задать его параметры
a_str(i,j).active = 1;
a_str(i,j).metabolism = ceil(rand * maxmetabolism);
a_str(i,j).vision = ceil(rand * maxvision);
a_str(i,j).wealth = s(i,j); % запас сахара данного агента
else % оставить клетку пустой
a_str(i,j).active = 0;
a_str(i,j).metabolism = 0;
a_str(i,j).vision = 0;
a_str(i,j).wealth = 0;
end
end
end
function displayScape(a_str, s, size)
% Отображение ландшафта и агентов
a = zeros(size);
for i = 1:size;
for j = 1:size;
if (a_str(i,j).active == 1)
a(i,j) = a_str(i,j).active;
end
end
end
imagesc(s);
alpha(0.7); hold on
spy(a,'ko'); axis('square','off'); hold off
drawnow
function [temps, tempi, tempj] = see(i,j,k,a_str,s,size,temps,tempi,tempj)
% Поиск наилучшей ячейки в поле зрения агента
south = [i+k size i+k-size j i+k j];
north = [k-i -1 i-k+size j i-k j];
east = [j+k size i j+k-size i j+k];
west = [k-j -1 i j-k+size i j-k];
c{1} = south; c{2} = north; c{3} = east; c{4} = west;
for m = 1:4
if (c{m}(1) > c{m}(2))
u = c{m}(3);
v = c{m}(4);
[temps, tempi, tempj] = neighbor(u,v,a_str,s,temps,tempi,tempj);
else
u = c{m}(5);
v = c{m}(6);
[temps, tempi, tempj] = neighbor(u,v,a_str,s,temps,tempi,tempj);
end
end
function [temps, tempi, tempj] = neighbor(u,v,a_str,s,temps,tempi,tempj)
% Если ячейка пуста и количество сахара в ней больше, чем в текущей,
% то запоминаем эту ячейку в качестве текущей
if (a_str(u,v).active == 0)
if (s(u,v) >= temps)
temps = s(u,v);
tempi = u;
tempj = v;
end
end
function a_str = moveAgent(a_str, s, i, j, temps, tempi, tempj)
% Перемещение агентов
if (temps > s(i,j))
% Агент переходит в соседнюю ячейку с наибольшим количеством сахара,
a_str(tempi,tempj) = a_str(i,j);
% освобождает ранее занятую клетку
a_str(i,j).active = 0;
a_str(i,j).vision = 0;
a_str(i,j).metabolism = 0;
a_str(i,j).wealth = 0;
% и обновляет свой запас сахара
a_str(tempi,tempj).wealth = ...
a_str(tempi,tempj).wealth + temps - a_str(tempi,tempj).metabolism;
% Если wealth <= 0, клетка считается пустой
if (a_str(tempi,tempj).wealth <= 0)
a_str(tempi,tempj).active = 0;
a_str(tempi,tempj).vision = 0;
a_str(tempi,tempj).metabolism = 0;
a_str(tempi,tempj).wealth = 0;
end
else
% Агент остается на месте и обновляет запас
a_str(i,j).wealth = a_str(i,j).wealth + temps - a_str(i,j).metabolism;
if (a_str(i,j).wealth <= 0)
a_str(i,j).active = 0;
a_str(i,j).vision = 0;
a_str(i,j).metabolism = 0;
a_str(i,j).wealth = 0;
end
end
Игра "Жизнь"
Правила описаны в Википедии.
clear all; clc
n = 64; % размеры КА
% создадим массивы
z = zeros(n,n);
c = z; % клетки КА
sum = z; % длЯ хранениЯ суммы живых клеток окрестности
% случайнаЯ начальнаЯ конфигурациЯ
%c = (rand(n,n))<.5 ;
% резонансное возбуждение
c(n/2-1:n/2+1,n/2) = 1;
c(n/2-1,n/2+1) = 1;
c(n/2,n/2-1) = 1;
% создадим картинку и нарисуем ее
imh = image(cat(3,c,c,z));
set(imh, 'erasemode', 'none') % рисование без морганий
axis equal
axis tight
% внутренние точки КА
x = 2:n-1;
y = 2:n-1;
% Основной цикл
ptime = 0.2; % времЯ задержки
for t=1:1000
% сумма живых клеток в окрестности данной
sum(x,y) = c(x,y-1) + c(x,y+1) + ...
c(x-1, y) + c(x+1,y) + ...
c(x-1,y-1) + c(x-1,y+1) + ...
c(3:n,y-1) + c(x+1,y+1);
% а вот и правила игры
c = (sum==3) | (sum==2 & c);
% рисуем картинку
set(imh, 'cdata', cat(3,c,c,z) )
pause(ptime) % задержка, если хотим рассмотреть подробности
% обновим счетчик шагов и покажем его значение на экране
drawnow % нужно чтобы рисовалось немедленно
end
Хижник-жертва
Основная программа:
% predator_prey.m
clear; clc
% Задаем начальные условия
t0 = 0;
tfinal = 15;
y0 = [20 20]';
% Численное интегрирование ДУ
tfinal = tfinal*(1+eps);
[t,y] = ode23(@lotka,[t0 tfinal],y0);
subplot(1,2,1)
plot(t,y)
title('Временная диаграмма')
subplot(1,2,2)
plot(y(:,1),y(:,2))
title('Фазовый портрет')
Функция в правой части уравнений:
% lotka.m
function yp = lotka(t,y)
% модель Лотка-Вольтерра системы "хищник-жертва"
yp(1) = 1*y(1) - .01*y(1)*y(2);
yp(2) =-1*y(2) + .02*y(1)*y(2);
yp = yp';
Пространственная модель "хищник-жертва"
Здесь описание получилось подлиннее. Приведем здесь только финальную программу. Ее и остальные программы можно скачать в архиве.
% serengeti4.m
function serengeti4
clc
nstep= 150; % Число итераций
size = 10; % Размер игрового мира (size x size)
grad = 1.5; % Градиент изменения растительности
dv = 0.1; % Прирост растительности за одну итерацию
% параметр = [жертва,хищник]
emax = [20,30]; % Максимальный запас энергии одной особи
eat = [1,0]; % Масса корма, съедаемого за одну итерацию
addenrg = [10,10]; % Прирост энергии за счет поедания корма
delenrg = [5,3]; % Расход энергии особью за одну итерацию
minenrg = [1,2]; % Минимальный запас энергии, необходимый для жизни
minrenrg = [2,5]; % Минимальный запас энергии, необходимый для воспроизводства
crmin = [10,40]; % Нижняя граница счетчика воспроизводства
crmax = [20,50]; % Верхняя граница счетчика воспроизводства
cpr1 = 0.5; % Вероятность поймать единичную особь жертвы
v = initVegetation(size, grad);
vmax = v;
% Задаем расположение и параметры агентов
[prey, pred, maxnumprey, maxnumpred] = initAgents(size, emax);
for step = 1:nstep
v(v < vmax) = v(v < vmax) + dv; % Прирост травы
todel = [];
r = countPopulation(size,pred);
for k = 1:length(prey) % Травоядные
prey(k).cr = prey(k).cr + 1;
if r(prey(k).i,prey(k).j) > 0 % В текущей ячейке хищник?
[inew, jnew] = selectNeighbour(size,r,prey(k).i,prey(k).j,2);
elseif (v(prey(k).i,prey(k).j) - eat(1)) > 0 % Корм в ячейке есть?
v(prey(k).i,prey(k).j) = v(prey(k).i,prey(k).j) - eat(1);
inew = prey(k).i; jnew = prey(k).j;
prey(k).energy = prey(k).energy + addenrg(1);
else
[inew, jnew] = selectNeighbour(size,v,prey(k).i,prey(k).j,1);
end
prey(k).i = inew;
prey(k).j = jnew;
prey(k).energy = prey(k).energy - delenrg(1);
if prey(k).energy < minenrg(1)
todel = [todel,k];
end
if prey(k).cr > crmax(1) % Счетчик воспроизводства превысил верхнюю границу?
prey(k).cr = 0;
elseif ( prey(k).cr >= crmin(1) && prey(k).cr <= crmax(1) ) && ...
prey(k).energy >= minrenrg(1)
[prey,maxnumprey] = reproduceAgents(prey,k,maxnumprey,emax(1));
prey(k).cr = 0;
end
end
prey(todel) = [];
todelp = [];
s = countPopulation(size,prey);
for p = 1:length(pred) % Хищники
pred(p).cr = pred(p).cr + 1;
if (s(pred(p).i,pred(p).j) - eat(2)) > 0 % Корм в ячейке есть?
cpr = cpr1; % Вероятность поимки
numdelprey = ceil( cpr*s(pred(p).i,pred(p).j) ); % Количество пойманных жертв
if numdelprey > 0
prey = delPrey(prey,pred(p).i,pred(p).j,numdelprey); % Уменьшение численности жертв
s = countPopulation(size,prey);
pred(p).energy = pred(p).energy + (numdelprey*addenrg(2) - delenrg(2));
end
else
[inew, jnew] = selectNeighbour(size,s,pred(p).i,pred(p).j,1);
pred(p).i = inew;
pred(p).j = jnew;
pred(p).energy = pred(p).energy - delenrg(2);
end
if pred(p).energy < minenrg(2)
todelp = [todelp,p];
end
if pred(p).cr > crmax(2) % Счетчик воспроизводства превысил верхнюю границу?
pred(p).cr = 0;
elseif ( pred(p).cr >= crmin(2) && pred(p).cr <= crmax(2) ) && ...
pred(p).energy >= minrenrg(2)
[pred,maxnumpred] = reproduceAgents(pred,p,maxnumpred,emax(2));
pred(p).cr = 0;
end
end
pred(todelp) = [];
mycmap = summer; colormap(mycmap(end:-1:1,:))
imagesc( v );
alpha(0.7); hold on
title( ['Число шагов = ', num2str(step)] )
[s,r] = displayAgents(size,prey,pred,1,0.2);
end
function x = initVegetation(s, A)
% Задаем растительный покров x:
% x = initVegetation(s, A)
% s - размер игрового мира
% A - градиент изменения растительности
[ X, Y ] = meshgrid( 1:s, 1:s );
x = A*X;
function [a,b,maxk,maxp] = initAgents(size, maxenergy)
% Задание расположения и параметров агентов
k = 1; p = 1;
for i = 1:size;
for j = 1:size;
if rand < 0.2 % поместить в клетку жертву и задать ее параметры
a(k).i = i;
a(k).j = j;
a(k).energy = maxenergy(1);
a(k).number = k;
a(k).cr = 0;
k = k + 1;
end
if rand < 0.05 % поместить в клетку хищника и задать его параметры
b(p).i = i;
b(p).j = j;
b(p).energy = maxenergy(2);
b(p).number = p;
b(p).cr = 0;
p = p + 1;
end
end
end
maxk = k - 1; maxp = p - 1;
function [in,jn] = selectNeighbour(size,a,i,j,key)
temp =-111*ones(3); % -111 означает, что клетка не принадлежит окрестности
I = 0;
for i1 = i-1:i+1
I = I + 1; J = 0;
for j1 = j-1:j+1
J = J + 1;
if i1 > 0 && j1 > 0 && i1 < size+1 && j1 < size+1
temp(I,J) = a(i1,j1);
end
end
end
switch key
case 1 % Поиск максимального значения
[u,v] = find( temp == max(max(temp)) & temp ~= -111 );
case 2 % Поиск пустой ячейки
[u,v] = find( temp == 0 );
end
n = length(u);
in = i + u( randperm(n) ) - 2;
jn = j + v( randperm(n) ) - 2;
in = in(1); jn = jn(1);
function [s,r] = displayAgents(size,a,b,txt,ptime)
s = countPopulation(size,a);
r = countPopulation(size,b);
spy(s,'w.',20);
spy(r,'r.',20); axis('square','off'); hold off
if txt ~= 0
for i = 1:size
for j = 1:size
if s(i,j) ~= 0, text( j,i,num2str(s(i,j)) ), end
if r(i,j) ~= 0, text( j,i,num2str(r(i,j)) ), end
end
end
end
pause(ptime)
function [a,maxk] = reproduceAgents(a,k,maxk,maxenergy)
m = length(a) + 1;
maxk = maxk + 1;
a(m).i = a(k).i;
a(m).j = a(k).j;
a(m).energy = maxenergy;
a(m).number = maxk;
a(m).cr = 0;
function x = countPopulation(size,a)
x = zeros(size);
n = length(a);
for k = 1:n
x(a(k).i,a(k).j) = x(a(k).i,a(k).j) + 1;
end
function a = delPrey(a,i,j,ndel)
n = length(a);
d = [];
for k = 1:n
if (a(k).i == i) && (a(k).j == j)
d = [d,k];
end
end
d = d( randperm(length(d)) );
d = sort(d,'descend');
for p = 1:ndel
a(d(p)) = [];
end
Комментарии
comments powered by Disqus