使用元胞自动机模拟森林火灾

问题情境：在一片森林中（有的地方为空地），当火灾发生时火灾蔓延到的地方如果是树，那树将被烧毁转化为空地；如果本身是空地，那就仍为空地。

编程思路

• 初始化空间，分别用0，1，2表示空地、树与火灾
• 根据上一时刻空间中各个格点的状态，确定这一时刻空间状态：
  • 如果上一时刻为空地，仍为空地
  • 如果上一时刻为火灾，这一时刻为空地
  • 如果上一时刻为树，这一时刻是否为树、火灾取决于上一时刻该位置（上下左右）是否有火灾，若有，则转为火灾；若无，则仍为树

使用Numpy生成随机数、进行函数运算；使用matplotlib进行绘图呈现

Python代码

引入Python库

# 引入Python库
import numpy as np # 生成数组
import numpy.random as rand # 随机整数
import matplotlib.pyplot as plt # 绘图
from IPython.display import display, clear_output # 显示 
import time # 时间

设置模拟空间

# 设置模拟空间
def set_board(board_size=50, f_trees_start=0.5):
    '''
    创建初始空间

    输入:
        board_size: 空间长度
        f_trees_start: 一个元胞含树的概率

    输出一个2d numpy数组，其中的值是0，1，2中的一个（分别代表空地，树，火灾）
    '''

    # 所有元胞的初始值都为0
    game_board = np.zeros((board_size, board_size), dtype='int64')


    # 在每个元胞处进行伯努利随机数生成，如果随机数小于f_tree_start，就设置为1（即树）
    for i in range(board_size):
        for j in range(board_size):
            if rand.random() <= f_trees_start:
                game_board[i, j] = 1

    # 我们将最左侧的一列都设为2
    game_board[:, 0] = 2

    return game_board

生成50 * 50 空间

set_board()

array([[2, 0, 1, ..., 1, 1, 1],
       [2, 1, 1, ..., 1, 1, 1],
       [2, 0, 0, ..., 0, 0, 1],
       ...,
       [2, 1, 1, ..., 0, 0, 1],
       [2, 1, 0, ..., 1, 1, 1],
       [2, 0, 1, ..., 1, 1, 0]], dtype=int64)

绘制空间

def plotgrid(myarray=np.random.choice((0,1,2),size=(50,50),p=(0.55,0.4,0.05))):
    # 设置网格
    plt.tick_params(axis='both', which='both',
        bottom=False, top=False, left=False, right=False,
        labelbottom=False, labelleft=False)
    # 初始化行和列的值
    x_range = np.linspace(0, myarray.shape[1] - 1, myarray.shape[1])
    y_range = np.linspace(0, myarray.shape[0] - 1, myarray.shape[0])

    # 使用numpy创建格点
    x_indices, y_indices = np.meshgrid(x_range, y_range)

    # 确定树与火灾的位置
    tree_x = x_indices[myarray == 1];
    tree_y = y_indices[myarray == 1];
    fire_x = x_indices[myarray == 2];
    fire_y = y_indices[myarray == 2];

    # 绘制树和火的图像
    plt.plot(tree_x, myarray.shape[0] - tree_y - 1, 'gs', markersize=10)
    plt.plot(fire_x, myarray.shape[0] - fire_y - 1, 'rs', markersize=10)

    # 设置图片上下限
    plt.xlim([-1, myarray.shape[1]])
    plt.ylim([-1, myarray.shape[0]])

随机生成一个森林图片

plotgrid()
plt.show()

其中白色代表空地，红色代表火灾，绿色代表树。

更新空间

随着火灾的蔓延，树木、空地的数量会发生改变。

# 更新空间
def advance_board(game_board):
    '''
    根据给定值推进游戏进行
    输入: 初始空间
    输出: 进一步的游戏空间
    '''

    # 设置一个初始值全为0的空间
    new_board = np.zeros_like(game_board)

    # 遍历每个元胞，然后确定该如何做
    for i in range(len(game_board)):
        for j in range(len(game_board)):
            place = game_board[i, j]
            # 如果上一轮值为0，就仍为0
            if place == 0:
                new_board[i, j] = 0
            # 如果上一轮火灾，那这一轮就是空地
            if place == 2:
                new_board[i, j] = 0
            # 如果上一轮是树，那么就需要再具体确定这一轮的状态
            if place == 1:
                # 在新的空间中还是将元胞设为树
                new_board[i, j] = 1
                # 如果上一轮有一个邻居格点着火
                # 那现在这个节点就是火灾
                if i > 0:
                    if game_board[i - 1, j] == 2:
                        new_board[i, j] = 2
                # 检查火灾是不是发生在右侧
                if i < game_board.shape[0] - 1:
                    if game_board[i + 1, j] == 2:
                        new_board[i, j] = 2
                # 检验火灾是不是发生在上侧
                if j > 0:
                    if game_board[i, j - 1] == 2:
                        new_board[i, j] = 2
                # 检查火灾是不是发生在下侧
                if j < game_board.shape[1] - 1:
                    if game_board[i, j + 1] == 2:
                        new_board[i, j] = 2
    # 返回新空间
    return new_board

计算状态

计算森林中树和空地的比例，如果连续两次数目、空地、火灾所占比例不发生改变，则认为状态稳定，程序终止。

# 计算状态
def calc_stats(game_board, start_board):
    '''
    计算空间中树与空地的比例

    输入: 状态空间

    输出: 空地和树的比例
    '''

    # 用numpy 空地数量占比
    frac_empty = np.sum(game_board[game_board == 0])/(game_board.shape[0]*game_board.shape[1])

    # 树的数量占比
    frac_tree = np.sum(game_board[game_board == 1]) /(game_board.shape[0]*game_board.shape[1])

    return frac_empty, frac_tree

主程序

开始模拟火灾。火灾从森林西侧（左）开始，逐步蔓延。

# 开始模拟

start_board = set_board() # 设置初始空间
f_trees_start = 0.7 # 设置树的比例
board_size = 50 # 设置空间大小

# 将图片大小设为(10,10)
fig = plt.figure(figsize=(10, 10))

# 初始化模拟空间
game_board = set_board(board_size=board_size, f_trees_start=f_trees_start)

# 绘制初始空间
plotgrid(game_board)

# 设置终止条件
on_fire = True

frac_empty_last,frac_trees_last = 0,0
frac_fire_last = 0
""# 开始循环
while on_fire == True:

    # 更新空间
    game_board = advance_board(game_board)

    # 呈现图片
    plotgrid(game_board)
    time.sleep(0.1)  # 停顿
    clear_output(wait=True) # 清空图片
    display(fig) # 呈现图片
    fig.clear() # 清空图片

    # 获得空地与树的数量
    frac_empty, frac_trees = calc_stats(game_board, start_board)
    frac_fire = 1 - frac_empty - frac_trees
    print(frac_empty,frac_trees)
    # 如果火灾，树，空地不发生改变，则终止
    if frac_empty_last == frac_empty and frac_trees_last==frac_trees_last and frac_fire == frac_fire_last:
        on_fire = False
    frac_empty_last,frac_trees_last,frac_fire_last = frac_empty, frac_trees,frac_fire
# 关闭绘图
plt.close()

初始森林图像（森林比例为0.7）如下：

演变过程：

更改数目密度(比例为0.5），更稀疏的森林，火灾蔓延范围更少：

更稠密森林(比例为0.9），蔓延范围更广：

当然，大家也可以尝试调整其他参数，观察火灾的变化情况。

# 引入Python库
import numpy as np # 生成数组
import numpy.random as rand # 随机整数
import matplotlib.pyplot as plt # 绘图
from IPython.display import display, clear_output # 显示 
import time # 时间

# 设置模拟空间
def set_board(board_size=50, f_trees_start=0.5):
    '''
    创建初始空间

    输入:
        board_size: 空间长度
        f_trees_start: 一个元胞含树的概率

    输出一个2d numpy数组，其中的值是0，1，2中的一个（分别代表空地，树，火灾）
    '''

    # 所有元胞的初始值都为0
    game_board = np.zeros((board_size, board_size), dtype='int64')


    # 在每个元胞处进行伯努利随机数生成，如果随机数小于f_tree_start，就设置为1（即树）
    for i in range(board_size):
        for j in range(board_size):
            if rand.random() <= f_trees_start:
                game_board[i, j] = 1

    # 我们将最左侧的一列都设为2
    game_board[:, 0] = 2

    return game_board

set_board()

array([[2, 0, 1, ..., 1, 1, 1],
       [2, 1, 1, ..., 1, 1, 1],
       [2, 0, 0, ..., 0, 0, 1],
       ...,
       [2, 1, 1, ..., 0, 0, 1],
       [2, 1, 0, ..., 1, 1, 1],
       [2, 0, 1, ..., 1, 1, 0]], dtype=int64)

def plotgrid(myarray=np.random.choice((0,1,2),size=(50,50),p=(0.55,0.4,0.05))):
    # 初始化行和列的值
    plt.tick_params(axis='both', which='both',
        bottom=False, top=False, left=False, right=False,
        labelbottom=False, labelleft=False)
    x_range = np.linspace(0, myarray.shape[1] - 1, myarray.shape[1])
    y_range = np.linspace(0, myarray.shape[0] - 1, myarray.shape[0])

    # 使用numpy创建格点
    x_indices, y_indices = np.meshgrid(x_range, y_range)

    # 确定树与火灾的位置
    tree_x = x_indices[myarray == 1];
    tree_y = y_indices[myarray == 1];
    fire_x = x_indices[myarray == 2];
    fire_y = y_indices[myarray == 2];

    # 绘制树和火的图像
    plt.plot(tree_x, myarray.shape[0] - tree_y - 1, 'gs', markersize=10)
    plt.plot(fire_x, myarray.shape[0] - fire_y - 1, 'rs', markersize=10)

    # 设置图片上下限
    plt.xlim([-1, myarray.shape[1]])
    plt.ylim([-1, myarray.shape[0]])
    # 显示网格

plotgrid()
plt.show()

# 更新空间
def advance_board(game_board):
    '''
    根据给定值推进游戏进行
    输入: 初始空间
    输出: 进一步的游戏空间
    '''

    # 设置一个初始值全为0的空间
    new_board = np.zeros_like(game_board)

    # 遍历每个元胞，然后确定该如何做
    for i in range(len(game_board)):
        for j in range(len(game_board)):
            place = game_board[i, j]
            # 如果上一轮值为0，就仍为0
            if place == 0:
                new_board[i, j] = 0
            # 如果上一轮火灾，那这一轮就是空地
            if place == 2:
                new_board[i, j] = 0
            # 如果上一轮是树，那么就需要再具体确定这一轮的状态
            if place == 1:
                # 在新的空间中还是将元胞设为树
                new_board[i, j] = 1
                # 如果上一轮有一个邻居格点着火
                # 那现在这个节点就是火灾
                if i > 0:
                    if game_board[i - 1, j] == 2:
                        new_board[i, j] = 2
                # 检查火灾是不是发生在右侧
                if i < game_board.shape[0] - 1:
                    if game_board[i + 1, j] == 2:
                        new_board[i, j] = 2
                # 检验火灾是不是发生在上侧
                if j > 0:
                    if game_board[i, j - 1] == 2:
                        new_board[i, j] = 2
                # 检查火灾是不是发生在下侧
                if j < game_board.shape[1] - 1:
                    if game_board[i, j + 1] == 2:
                        new_board[i, j] = 2
    # 返回新空间
    return new_board

# 计算状态
def calc_stats(game_board, start_board):
    '''
    计算空间中树与空地的比例

    输入: 状态空间

    输出: 空地和树的比例
    '''

    # 用numpy 空地数量占比
    frac_empty = np.sum(game_board[game_board == 0])/(game_board.shape[0]*game_board.shape[1])

    # 树的数量占比
    frac_tree = np.sum(game_board[game_board == 1]) /(game_board.shape[0]*game_board.shape[1])

    return frac_empty, frac_tree

# 开始模拟

start_board = set_board() # 设置初始空间
f_trees_start = 0.9 # 设置树的比例
board_size = 50 # 设置空间大小

# 将图片大小设为(10,10)
fig = plt.figure(figsize=(10, 10))

# 初始化模拟空间
game_board = set_board(board_size=board_size, f_trees_start=f_trees_start)

# 绘制初始空间
plotgrid(game_board)

# 设置终止条件
on_fire = True

frac_empty_last,frac_trees_last = 0,0
frac_fire_last = 0
""# 开始循环
n = 0
while on_fire == True:

    # 更新空间
    game_board = advance_board(game_board)

    # 呈现图片
    plotgrid(game_board)

    clear_output(wait=True) # 清空图片
    display(fig) # 呈现图片
    time.sleep(0.1)  # 停顿
    if n % 10 == 0:
        time.sleep(3)
    n += 1
    fig.clear() # 清空图片

    # 获得空地与树的数量
    frac_empty, frac_trees = calc_stats(game_board, start_board)
    frac_fire = 1 - frac_empty - frac_trees
    print(frac_empty,frac_trees)
    # 如果火灾，树，空地不发生改变，则终止
    if frac_empty_last == frac_empty and frac_trees_last==frac_trees_last and frac_fire == frac_fire_last:
        on_fire = False
    frac_empty_last,frac_trees_last,frac_fire_last = frac_empty, frac_trees,frac_fire
# 关闭绘图
plt.close()

0.0 0.0004

game_board = set_board(board_size=board_size, f_trees_start=f_trees_start)

# 绘制初始空间
plotgrid(game_board)
plt.show()