ABM-基于简单几何形状定义个体活动范围

模型假设与符号定义

• 假设个体在二维平面的一个正方形有界闭区域内运动，区域中没有个体的流入和流出。

• 假设初始时刻个体在区域内均匀分布。

• 假设区域的上边界与下边界连通，左边界与右边界连通，因此个体在撞上区域的边界时，会重新回到该区域内。

• 假设当易感个体与具有传染性的个体之间的距离小于某一个特定的阈值时，可能发生疾病的传播，传播的概率是一个固定的常数。

• 假设该疾病的潜伏期、传染期服从正态分布。

• 假设个体在每一个单位时间内会以一定的概率决定是否发生移动，并且等可能的北、东北、东、东南、南、西南、西、西北八个方向中随机选择一个进行移动，移动的距离服从正态分布。

• 假设个体在被感染之后获得终生免疫力，不会发生二次感染。

• 不考虑出生率和个体的自然死亡率。

个体属性和人群状态初始化

个体行为

个体位移更新规则

个体在每一个单位时间内会以一定的概率（moveProbability）决定是否发生移动，如果发生移动：

移动方向：等可能地从北、东北、东、东南、南、西南、西、西北八个方向中随机选择一个

移动距离：在均数是 moveAverage，标准差是 moveStd 的正态分布中随机抽样

接触传播

我们假设每个个体的活动范围是一个半径为 distanceContact 的圆，并且把两个个体的活动范围有交集视为接触，所以当易感个体与具有传染性的个体之间的距离小于某一个特定的阈值（2*distanceContact）时，认为该易感个体会以一定概率（$q$）被感染，如果被感染，该个体的疾病状态从易感变为感染但未具有传染性。

患病状态之间的转化

当感染但未具有传染性的个体在该状态下停留的时间大于等于该个体的潜伏期时，该个体状态变为感染且具有传染性；当感染且具有传染性的个体在该状态下停留的时间大于等于该个体的传染期时，该个体状态变为被移除，不会造成传播也不会被感染，且此状态不会再发生改变。

环境

个体的活动范围是二维平面上的一个正方形有界闭区域，因此我们把个体在笛卡尔坐标系中的横坐标限制在一个用行向量 $\mathbb{R}^{2}$（xlimit）表示的闭区间内，纵坐标限制在一个用行向量 $\mathbb{R}^{2}$（ylimit）表示的闭区间内。并且区域的上边界与下边界连通，左边界与右边界连通，因此个体在撞上区域的边界时，会重新回到该区域内。

代码

%函数
function [quadtreeCoordinates] = quadtreePartition(xlimit, ylimit, x, y, index, quadtreeCoordinates0, ending, distanceContact)
%QUADTREEPARTITION 使用四叉树分割优化接触检测
% input
%xlimit 表示在二维直角坐标系中横坐标范围的1*2的行向量，xlimit(2)大于等于xlimit(1)
%ylimit 表示在二维直角坐标系中纵坐标范围的1*2的行向量，ylimit(2)大于等于ylimit(1)
%x 每个个体的横坐标组成的列向量
%y 每个个体的纵坐标组成的列向量，x与y大小一致
%index 输入的坐标在最初的坐标列向量中的索引行向量
%quadtreeCoordinates0 quadtreeCoordinates的初始元胞数组
%ending 叶节点的最大对象数，用于表示进行分割的终止条件，整数标量，当某一区域内部的个体数量少于等于ending时，停止分割
%distanceContact 表示接触判定阈值的标量，当两个个体之间的距离小于2*distanceContact时，认为他们发生接触
% output
%quadtreeCoordinates N*1的元胞数组，四叉树有损地将二维空间中每个点的坐标保存为d个四分类变量的笛卡尔积，d>=1且d<=D，每一个元胞都包含一个1*d的表示笛卡尔积的行向量

quadtreeCoordinates = quadtreeCoordinates0;
%进行一次四叉树分割
%个体无法完美划分（以个体坐标为圆心，distanceContact为半径的圆与分割线相交）时，保留在当前节点中
idx0 = (abs(x - (xlimit(2)+xlimit(1))./2) < distanceContact) | (abs(y - (ylimit(2)+ylimit(1))./2) < distanceContact);
x(idx0) = [];
y(idx0) = [];
index(idx0) = [];
%个体划分到第一象限，坐标记为1
idx1 = (x <= xlimit(2) & x >= (xlimit(2)+xlimit(1))./2) & (y <= ylimit(2) & y >= (ylimit(2)+ylimit(1))./2);
for i = index(idx1)
    quadtreeCoordinates{i} = [quadtreeCoordinates{i}, 1];
end
%个体划分到第二象限，坐标记为2
idx2 = (x >= xlimit(1) & x <= (xlimit(2)+xlimit(1))./2) & (y <= ylimit(2) & y >= (ylimit(2)+ylimit(1))./2);
for i = index(idx2)
    quadtreeCoordinates{i} = [quadtreeCoordinates{i}, 2];
end
%个体划分到第三象限，坐标记为3
idx3 = (x >= xlimit(1) & x <= (xlimit(2)+xlimit(1))./2) & (y >= ylimit(1) & y <= (ylimit(2)+ylimit(1))./2);
for i = index(idx3)
    quadtreeCoordinates{i} = [quadtreeCoordinates{i}, 3];
end
%个体划分到第四象限，坐标记为4
idx4 = (x <= xlimit(2) & x >= (xlimit(2)+xlimit(1))./2) & (y >= ylimit(1) & y <= (ylimit(2)+ylimit(1))./2);
for i = index(idx4)
    quadtreeCoordinates{i} = [quadtreeCoordinates{i}, 4];
end

%递归：在函数的定义中调用自身
%对于四叉树这种具有明显递归性质的数据结构，递归算法可以使代码更加简洁，可读性也更好
%这是一个多路递归，该函数在递归调用时，同时调用了四个子问题
if sum(idx1) > ending %终止条件
    [quadtreeCoordinates] = quadtreePartition([(xlimit(2)+xlimit(1))./2, xlimit(2)], [(ylimit(2)+ylimit(1))./2, ylimit(2)], x(idx1), y(idx1), index(idx1), quadtreeCoordinates, ending, distanceContact);
end
if sum(idx2) > ending %终止条件
    [quadtreeCoordinates] = quadtreePartition([xlimit(1), (xlimit(2)+xlimit(1))./2], [(ylimit(2)+ylimit(1))./2, ylimit(2)], x(idx2), y(idx2), index(idx2), quadtreeCoordinates, ending, distanceContact);
end
if sum(idx3) > ending %终止条件
    [quadtreeCoordinates] = quadtreePartition([xlimit(1), (xlimit(2)+xlimit(1))./2], [ylimit(1), (ylimit(2)+ylimit(1))./2], x(idx3), y(idx3), index(idx3), quadtreeCoordinates, ending, distanceContact);
end
if sum(idx4) > ending %终止条件
    [quadtreeCoordinates] = quadtreePartition([(xlimit(2)+xlimit(1))./2, xlimit(2)], [ylimit(1), (ylimit(2)+ylimit(1))./2], x(idx4), y(idx4), index(idx4), quadtreeCoordinates, ending, distanceContact);
end

end

%脚本
clc %清理命令行窗口
clear %清理工作区变量
close all %关闭所有图窗

%% 在二维平面上设置简单几何形状的横纵坐标范围

% 正方形 四个端点坐标为(0,0)、(100,0)、(100,100)、(0,100)
xlimit = [0, 100];
ylimit = [0, 100];

% 长方形 四个端点坐标为(0,0)、(100,0)、(100,75)、(0,75)
%xlimit = [0, 100];
%ylimit = [0, 75];

% 平行四边形 四个端点坐标为(50,0)、(100,0)、(50,100)、(0,100)
%xlimit = [0, 100];
%x是从区间xlimit的均匀分布中得到的随机数
%if x <= 50
%   ylimit = [100-2.*x, 100];
%else
%   ylimit = [0, 200-2.*x];

% 圆 圆心(50,50)、半径为50
%xlimit = [0, 100];
%x是从区间xlimit的均匀分布中得到的随机数
%ylimit = [-sqrt(2500-(x-50).^2) + 50, sqrt(2500-(x-50).^2) + 50];

% 椭圆 圆心(50,25)、长轴为2*50=100、短轴为2*25=50、焦点连线与x轴平行
%xlimit = [0, 100];
%x是从区间xlimit的均匀分布中得到的随机数
%ylimit = [-sqrt(625-(x-50).^2./4) + 25, sqrt(625-(x-50).^2./4) + 25];

% 椭圆 圆心(25,50)、长轴为2*50=100、短轴为2*25=50、焦点连线与y轴平行
%xlimit = [0, 50];
%x是从区间xlimit的均匀分布中得到的随机数
%ylimit = [-sqrt(2500-(x-25).^2.*4) + 50, sqrt(2500-(x-25).^2.*4) + 50];

%% 将人群状态向量化，使用简单数组表示，向量默认为列向量

rng(626, 'twister'); %指定生成器算法和种子来重复生成随机数数组。每次使用相同算法和种子初始化生成器时，始终都可以获得相同的结果。
N = 200; %个体总数
Ts = 1; %时间步长为1天
simulationDuration = 100; %模拟总时长，单位是天
q = 0.5; %单次接触感染概率
distanceContact = 1; %接触判定阈值，当两个个体之间的距离小于2*distanceContact时，认为他们发生接触
initialCases = 20; %初始病例数
moveProbability = 0.7; %个体位置发生改变的概率
directionMatrix = [0,1; sqrt(2).*0.5,sqrt(2).*0.5; 1,0; sqrt(2).*0.5,-sqrt(2).*0.5; 0,-1; -sqrt(2).*0.5,-sqrt(2).*0.5; -1,0; -sqrt(2).*0.5,sqrt(2).*0.5];
%一个8*2的矩阵，每行表示二维平面的一个方向，且每个行向量的长度为1
moveAverage = 3;
moveStd = 0.5; %个体在一个时间步长内移动的距离服从均值为3，标准差为0.5的正态分布
stage = ones(N, 1, "int8"); %疾病状态向量，1-易感，2-暴露，3-传染，4-恢复
incubationPeriod = randn(N, 1, "single") .* 1 + 7; %潜伏期向量，假设潜伏期服从均值为7，标准差为1的正态分布，单位是天
incubationPeriod(incubationPeriod < 0) = 0;
infectiousPeriod = randn(N, 1, "single") .* 2 + 10; %传染期向量，假设传染期服从均值为10，标准差为2的正态分布，单位是天
infectiousPeriod(infectiousPeriod < 0) = 0;
durationSinceEnteringStage = zeros(N, 1, "single"); %进入各个状态的时间向量

%以正方形为例，在二维平面的正方形区域中随机生成每个个体的横纵坐标
x = rand(N, 1, "single") .* 100;
y = rand(N, 1, "single") .* 100;
%设置人群中的初始病例
stage(randperm(N, initialCases)) = 3;
%初始化人群在不同时刻下的状态矩阵
simulatedResults = zeros(simulationDuration+1, 4, "int16");
simulatedResults(1,:) = [N-initialCases, 0, initialCases, 0];
%可视化人群初始状态
figure
xline(100)
yline(100)
scatter(x(stage == 1), y(stage == 1), "MarkerEdgeColor", [0 0.4470 0.7410])
hold on
scatter(x(stage == 3), y(stage == 3), "MarkerEdgeColor", [0.9290 0.6940 0.1250])
hold off

%% 开始模拟

figure
xline(100)
yline(100)
for i = 1:simulationDuration
  
    %进行四叉树分割，寻找与传染者接触的易感者，并判断其是否发生感染
    quadtreeCoordinates0 = cell(N,1);
    [quadtreeCoordinates] = quadtreePartition(xlimit, ylimit, x, y, 1:N, quadtreeCoordinates0, 10, distanceContact);        
    caseIndex = stage == 3;
    CaseQuadtreeCoordinates = quadtreeCoordinates(caseIndex);
    CaseXCoordinate = x(caseIndex);
    CaseYCoordinate = y(caseIndex);
    for j = 1:numel(CaseQuadtreeCoordinates)        
        CaseNeighborIndex = false(N, 1);
        currentCase = CaseQuadtreeCoordinates{j};
        %从四叉树根部开始搜索可能与传染者接触的个体
        if isempty(currentCase)
            CaseNeighborIndex = true(N,1);
        else
            for k = 1:N
                if isempty(quadtreeCoordinates{k})
                    CaseNeighborIndex(k) = true;
                elseif numel(quadtreeCoordinates{k}) >= numel(currentCase)
                    if isequal(currentCase, quadtreeCoordinates{k}(1:numel(currentCase)))
                        CaseNeighborIndex(k) = true;
                    end
                else
                    for m = 1:(numel(currentCase)-1)
                        if isequal(currentCase(1:m), quadtreeCoordinates{k})
                            CaseNeighborIndex(k) = true;
                            break
                        end
                    end
                end
            end
        end
        %判断是否发生感染
        CaseNeighborIndex1 = find(CaseNeighborIndex);
        newCaseIndex = (sqrt((x(CaseNeighborIndex) - CaseXCoordinate(j)).^2 + (y(CaseNeighborIndex) - CaseYCoordinate(j)).^2) < 2*distanceContact) & (rand(sum(CaseNeighborIndex),1) <= q) & (stage(CaseNeighborIndex) == 1);
        stage(CaseNeighborIndex1(newCaseIndex)) = 2;
        durationSinceEnteringStage(CaseNeighborIndex1(newCaseIndex)) = 0;
    end
    
    %更新个体的位置
    moveIndex = find(rand(N,1) <= moveProbability);
    directionIndex = randi(8, [1, numel(moveIndex)]); %每个个体的移动方向有8种可能：1-北、2-东北、3-东、4-东南、5-南、6-西南、7-西、8-西北
    x(moveIndex) = x(moveIndex) + (randn(numel(moveIndex),1).*moveStd + moveAverage) .* directionMatrix(directionIndex, 1);
    y(moveIndex) = y(moveIndex) + (randn(numel(moveIndex),1).*moveStd + moveAverage) .* directionMatrix(directionIndex, 2);
    temp = x(moveIndex);
    temp(temp > xlimit(2)) = xlimit(1) + temp(temp > xlimit(2)) - xlimit(2);
    temp(temp < xlimit(1)) = xlimit(2) + temp(temp < xlimit(1)) - xlimit(1);
    x(moveIndex) = temp;
    temp = y(moveIndex);
    temp(temp > ylimit(2)) = ylimit(1) + temp(temp > ylimit(2)) - ylimit(2);
    temp(temp < ylimit(1)) = ylimit(2) + temp(temp < ylimit(1)) - ylimit(1);
    y(moveIndex) = temp;

    %更新自然史状态
    idx1 = (stage == 2) & (durationSinceEnteringStage >= incubationPeriod);
    stage(idx1) = 3;
    durationSinceEnteringStage(idx1) = durationSinceEnteringStage(idx1) - incubationPeriod(idx1);
    idx2 = (stage == 3) & (durationSinceEnteringStage >= infectiousPeriod);
    stage(idx2) = 4;
    durationSinceEnteringStage(idx2) = durationSinceEnteringStage(idx2) - infectiousPeriod(idx2);

    %时间增加
    durationSinceEnteringStage = durationSinceEnteringStage + Ts;

    %记录此时人群的状态并可视化
    simulatedResults(i+1,:) = [sum(stage == 1), sum(stage == 2), sum(stage == 3), sum(stage == 4)];
    scatter(x(stage == 1), y(stage == 1))
    hold on
    scatter(x(stage == 2), y(stage == 2))
    scatter(x(stage == 3), y(stage == 3))
    scatter(x(stage == 4), y(stage == 4))
    hold off
end
figure
plot(simulatedResults)
legend(["S", "E", "I", "R"])

模拟设置

模拟结果与解释

模拟用时

迭代次数 100，总时长 4.636141s，单次迭代用时 0.04636s

模拟结果

散点图中，每个圆圈代表一个个体，不同颜色代表不同的疾病状态，蓝色代表易感，橙色代表感染但未具有传染性，黄色代表感染且具有传染性，紫色代表被移除，不会造成传播也不会被感染。

t=0

t=31

t=61

t=100

每天的新感染人数，在第20天左右达到峰值

每天的新恢复人数

各类人群随时间的变化趋势

从上图可以看出两类感染人群的数量一直维持在0到20之间，在第80天之后感染者逐渐消失，最终这次疫情造成77人感染（包括20名初始病例）。该设置条件下疾病并没有引起广泛的传播。

参数变化对于模拟结果的影响

传染期较长时，在封闭空间中，每个感染者有更多的时间将疾病传播给其他人，这会导致疾病传播的范围更广，感染人数的峰值更高，并且疫情持续的时间可能也会更长。相反，传染期较短时，疾病传播的机会相对减少，感染人数的增长速度会比较慢，疫情可能在较短时间内得到控制。

潜伏期长意味着感染者在较长时间内可能没有传染性，发病人数的变化曲线会出现多个峰，每一代感染者之间有明显的时间间隔。而在潜伏期短的情况则只会出现一个峰。

较小的时间步长可以更精细地模拟疾病传播的过程，能够更准确地捕捉到人群中短时间内的频繁接触和疾病传播情况，但这也会增加计算量和模拟的复杂性。大的时间步长会简化模拟过程，但可能会丢失一些细节，低估疾病的传播速度和范围。

总人数较多时，疾病传播的潜在机会就会增加。即使疾病的传染性不是很强，由于人口密度大，人与人之间接触的可能性高，疾病传播的范围也可能很广。感染人数的增长可能会呈现指数型增长，疫情可能会比较严重，而且持续时间可能较长。总人数较少时，疾病传播的范围相对有限。因为接触的人数有限，感染人数的增长可能会比较缓慢。