1, Introduction to cellular automata

Development of cellular automata
The original cellular automata was proposed by von Neumann in the 1950s to simulate the self replication of biological cells However, it has not been valued by the academic community
In 1970, after John Horton Conway of Cambridge University designed a computer game "life game", cellular automata attracted the attention of scientists

In 1983, S.Wolfram published a series of papers The models generated by 256 rules of elementary cellular machines are deeply studied, and their evolution behavior is described by entropy. Cellular automata are divided into stationary type, periodic type, chaotic type and complex type

2. Preliminary understanding of cellular automata
Cellular automata (CA) is a method used to simulate local rules and local connections. A typical cellular automata is defined on a grid. The grid at each point represents a cell and a finite state. Change rules apply to each cell and run simultaneously. Typical change rules depend on the state of the cell and its state (4 or 8) the status of neighbors.

The change rule of 3-cell & cell state
Typical change rules depend on the state of the cell and the state of its (4 or 8) neighbors.

Application of 4-cellular automata
Cellular automata has been applied to physical simulation, biological simulation and other fields.

matlab programming of 5-cell cellular automata
Combined with the above, we can understand that cellular automata simulation needs to understand three points. One is cell. In matlab, it can be understood as a square block composed of one or more points in the matrix. Generally, we use a point in the matrix to represent a cell. The second is the change rule, which determines the state of the cell at the next moment. The third is the state of cells. The state of cells is user-defined and usually opposite, such as the living state or death state of organisms, red light or green light, obstacles or no obstacles at this point, etc.

6 one dimensional cellular automata -- traffic rules
definition:
6.1 cells are distributed on one-dimensional linear grid
6.2 cells have only two states: car and empty

7 two dimensional cellular automata -- life game
definition:
7.1 cells are distributed on two-dimensional square grid
7.2 cells have only two states of life and death

The cellular state is determined by the surrounding eight neighbors
Rules:

Skeleton: death; Smiling face: survival
If there are three smiling faces around, the middle becomes a smiling face
Less than two smiling faces or more than three, the middle becomes death.

8 what is cellular automata
Discrete system: cell is defined in finite time and space, and the state of cell is finite
Dynamic system: the behavior of cellular automata has dynamic characteristics
Simplicity and complexity: cellular automata uses simple rules to control interacting cells to simulate a complex world

9 constituent elements

(1) Cell

Cell is the basic unit of cellular automata:
State: each cell has the function of memory storage state
Discrete: in simple cases, there are only two possible states of cells; In complex cases, cells have many states
Update: the state of the cell is constantly updated according to the dynamic rules
(2) Mesh (Lattice)
Different dimensional grid

Common 2D mesh

(3) Neighbors

(4) Boundary

Reflective: a state with itself as a boundary
Absorption type: regardless of the boundary (the car disappears when it reaches the boundary)

(5) Rule (state transition function)
Definition: determine the dynamic function of the cell state at the next time according to the current state of the cell and its neighbors. In short, it is a state transition function
Classification:
Summation: the state of a cell depends on and only depends on the current state of all its neighbors and its own current state
Legal type: sum type rules are legal type rules However, if the rules of cellular automata are limited to summation, cellular automata will have limitations
(6) Forest fire

Green: trees; Red: fire; Black: open space.
Three state cycle transitions:
Tree: when there is fire around or struck by lightning, it becomes fire.
Open space: change into trees with probability p
Rational analysis: red is fire; Ash is open space; Green is a tree

The density sum of the three states of the cell is 1

The density of fire into open space is equal to the density of open space into trees (newly grown trees are equal to burned trees)

f is the probability of lightning: far less than the probability of tree generation; T s m a x T_{smax}T smax
It's a time scale for a large group of trees to be burned
Program implementation
Periodic boundary conditions
Buy

The number is number
Building neighbor matrix

The number in the above matrix corresponds to the upper neighbor number of the same position number of the original matrix, one by one
Similarly:

(7) Traffic concept
Distance and density

Flow equation

conservation equation

Spatiotemporal trajectory (horizontal axis is space and vertical axis is time)

The intersection of the red line and the blue line indicates the location of each time car.
If it is a vertical line, it indicates the corresponding time of the vehicle in this position

Macro continuous model:

Most common rules:

The red bar indicates that the speed is full.

1 acceleration rule: no more than v m a x (2 grids / s) v#{max} (2 grids / s) V
max (2 grids / s)
2. Collision prevention: do not exceed the vehicle distance

Theoretical analysis:

Result analysis: density and flow

The first figure: the abscissa is the normalized density, and the ordinate is the traffic flow. The second figure: theoretical values and CA results

Result analysis: spatiotemporal trajectory

The dark area in the middle is the area of traffic jam.

2, Partial source code

function [ ] = main(~)
%Draw the motion diagram
A=zeros(113,26);
M=ceil(sqrt(113^2+26^2));
A1=zeros(29,22);
A1(:,[1 22])=M;
A1([1 29],:)=M;
A1([4 5],[10 11 12 13])=M;%Give the number of obstacles on the podium
A1(8:2:24,[4 5 6 10 11 12 13 17 18 19])=M;%Give the obstacle amount of the table
A1(9:2:25,[5 11 12 18])=A(9:2:25,[5 11 12 18]);%Give the distance from the middle seat of the adjacent seat to the door,It's removed here+1 Situation
A1(5,1)=0;%Give me a place to go out
A([1,29,57,85],5:26)=M;
A(113,:)=M;
A(:,[1,5,26])=M;
A(1:29,5:26)=A1;
A(29:57,5:26)=A1;
A(57:85,5:26)=A1;
A(85:113,5:26)=A1;%So far, the shape of the classroom, the distribution of obstacles and the position of the door are given
C=zeros(113,26);
C1=zeros(29,22);
C1(9:2:25,[4 5 6 10 11 12 13 17 18 19])=M;

C(85:113,5:26)=C1;%So far, the personnel distribution of the whole floor is given, as shown in the matrix C As shown in, the initial situation is given at this time, and everyone is sitting in the seat
for g=1:340;%200 This time is sure to end. You can try it while，I'm too lazy here
imshow(max(A,C)==0,'InitialMagnification','fit')%Take the largest as the distribution of personnel in the classroom, and give the classroom outline and obstacle image
%The next two parameters are to resize the image
pause(0.25);
C=run(C);
end
function [D]=choic(D)
A=zeros(29,22);%Generate a zero matrix with 29 rows and 22 columns, and the processed zero matrix A The matrix will represent the desire of each cell for the door
B=ones(25,22);%A matrix with 25 rows and 22 columns at each position of 1 is generated for processing
A(5:29,:)=sqrt((cumsum(B,1)-1).^2+(cumsum(B,2)-1).^2);%Use calculation skills to calculate the distance between each cell and the door
C=ones(4,22);%Generate a matrix of 4 rows and 22 columns with each position of 1
A(1:4,:)=sqrt((cumsum(C,1)-5).^2+(cumsum(C,2)-1).^2);%Use calculation skills to calculate the distance between each cell and the door
M=ceil(sqrt(113^2+26^2));%Take the upward rounding of the diagonal of the whole floor as the farthest distance. Here is also a fine adjustment to spell the two programs together after writing the corridor. Otherwise, take the diagonal of the room upward rounding
A(:,[1 22])=M;
A([1 29],:)=M;%Give the amount of obstruction on the wall
A([4 5],[10 11 12 13])=M;%Give the number of obstacles on the podium
A(8:2:24,[4 5 6 10 11 12 13 17 18 19])=M;%Give the obstacle amount of the table
A(9:2:25,[5 11 12 18])=A(9:2:25,[5 11 12 18]);%Give the distance from the middle seat of the adjacent seat to the door,It's removed here+1 Situation
A(5,1)=0;%Give the location of the door
%D=zeros(29,22);%Another table is given to show that in classrooms of the same size
%D(9:2:25,[4 5 6 10 11 12 13 17 18 19])=M;%Give the location of the people in the classroom
%The matrix is given here D It is convenient to call when verifying the correctness of program writing
X=zeros(29,22); %Give the target matrix, that is, the choice made by people and the position to be moved in the next moment
Y=zeros(29,22); %Give out what you want to move to X The same goal in( x,y)The location of the person,The value above is x*sqrt(2)+y*sqrt(3)Indicates that the target is coming( x,y)go
Z=rand(29,22);%Generate a random number table for final decision-making and judgment, and pay attention to the position of this table
H=zeros(29,22);%Generating matrix H In order to be the position of the person who finally succeeded in grabbing the position
%A(5,1)=0;%Give initial value
%D(5,1)=0;%give D The location of the door
%X(5,1)=0;
%Y(5,1)=0;
%Z(5,1)=0;
%These quantities are to make the people at the door of the room go out when considering a room alone. If considering the whole floor, the work of making the people in the classroom go out is placed in run In the function
for x=2:28;%give x Cycle range of
for y=2:21;%give y Cycle range of
if D(x,y)==M %If the center of the nine palaces is human, the following operations can be carried out
E=max(A,D);%take A and D The maximum of the corresponding value in indicates that the person is also in the classroom at the moment, and the value of each cell except his own position indicates the distance from the person to the door
F=E((x-1):1:(x+1),(y-1):(y+1));%Take out to( x,y)Centered Jiugong grid
F(2,2)=A(x,y);%The center of the Jiugong grid has an impact on others but not on yourself, so change back to the original value
G=sort(F);%Sort columns
b=G(1,:);%Take row vector
d=sort(b);%Sort row vectors
if length(find(F==d(1)))==1 %First consider the case where the minimum value found is 1
[r,c]=find(F==d(1)); %Find the location of the minimum value
if r==2&&c==2 %The minimum value is just in the center
X(x,y)=D(x,y);%The objective matrix is itself and remains unchanged
Y(x,y)=x*sqrt(2)+y*sqrt(3);%The arrival target matrix is given( x,y)Value at position of
else p=x-2+r; %When the minimum value is not in the center, the minimum value position found is the target position
q=y-2+c;
if p~=1
X(p,q)=M;%Give the objective matrix
Y(x,y)=p*sqrt(2)+q*sqrt(3);%The arrival target matrix is given( x,y)Value at position of
else X(x,y)=D(x,y);
Y(x,y)=x*sqrt(2)+y*sqrt(3);
end
end
else length(find(F==d(1)))==2; %Consider the simultaneous occurrence of two minimum values
[r,c]=find(F==d(1));
s=rand(1);
if s>0.5
p=x-2+r(1);
q=y-2+c(1);
if p~=1;
X(p,q)=M;
Y(x,y)=p*sqrt(2)+q*sqrt(3);
else X(x,y)=D(x,y);
Y(x,y)=x*sqrt(2)+y*sqrt(3);
end
else
p=x-2+r(2);
q=y-2+c(2);
if p~=1;
X(p,q)=M;
Y(x,y)=p*sqrt(2)+q*sqrt(3);
else X(x,y)=D(x,y);
Y(x,y)=x*sqrt(2)+y*sqrt(3);
end
end
end
end 
end
end
for x=2:28;%give x Cycle range of
for y=2:22;%give y Cycle range of
if X(x,y)>0;
Y1=Y((x-1):1:(x+1),(y-1):(y+1));%To target( x,y)Give the nine palace grid
Y2=Z((x-1):1:(x+1),(y-1):(y+1));%The random number of the Jiugong lattice is given as the subsequent judgment
w=x*sqrt(2)+y*sqrt(3);
%(Y1==w)*Y2;%Find the random number corresponding to the person with the center as the goal from the nine house lattice
t=max(max((Y1==w).*Y2));%Find the maximum value of the corresponding position of the person with the center as the target
[r1,c1]=find((Y1==w).*Y2==t);%Find the location of the maximum value

end
end
%To (5),1)How should the target be discussed
for x=2:28;%give x Cycle range of
for y=1;%give y Cycle range of
if X(x,y)>0;

[r1,c1]=find((Y1==w).*Y2==t);%Find the location of the maximum value
H(x-2+r1,c1)=M;%Will grab the center position in H Change the position value in to M
end 
end
end
D=D+X-H;%A very simple formula is applied to give the distribution of people in the classroom in the next time step
end

3, Operation results

4, matlab version and references

1 matlab version
2014a

2 references
[1] Steamed stuffed bun Yang, Yu Jizhou, Yang Shan Intelligent optimization algorithm and its MATLAB example (2nd Edition) [M]. Electronic Industry Press, 2016
[2] Zhang Yan, Wu Shuigen MATLAB optimization algorithm source code [M] Tsinghua University Press, 2017
[3][mathematical modeling] cellular automata. Blogger: binary artificial intelligence