简介：DESSABNeT 是一个基于 MATLAB 编写的开源传染病传播模拟软件包，采用基于代理的模型（ABM）构建人工社会，以真实人口统计数据为输入，可模拟任何城市或区域中传染病的传播动态，最初专门用于澳大利亚城市的 COVID-19 疫情建模。面向流行病学研究者、公共卫生决策者及计算社科领域的 MATLAB 用户，可用于评估社交限制措施、模拟不同疾病爆发场景。资源包共 42 个文件，以 33 个 .m 格式 MATLAB 脚本为核心，覆盖人口设置、社交网络构建、每周时间表配置、疾病传播动态模块及社交限制阶段模拟；另有 7 个 .mat 预生成数据文件、1 个 Excel 数据表和 1 个 README 说明文档，压缩包整体 2.7MB，按模块分目录组织。已有 712 人学习下载。读者可获得完整的 ABM 建模框架代码，包括悉尼人口结构设置、家庭/工作/学校网络构建、感染病例注入与传播计算、R0 矩阵计算、隔离与社交限制模拟等，并提供可直接运行的示例数据，便于移植到其他传染病或城市场景开展仿真实验。1. 这份MATLAB模型，为什么值得你先跑一遍再说做传染病仿真的人多数有一个共识：SEIR已经不够用了。真实世界每一轮疫情都会出现死亡病例，而标准SEIR只把人群分成易感、暴露、感染、恢复四类，死亡和康复只占一个R槽，模拟出来反而比实际结果乐观不少。这套名为DESSABNeT的SEIRD模型MATLAB代码，正好补上了这段缺口——它把死亡（D）单独拆出来，同时把消毒网（Disinfection Network）的影响以参数形式接入传播链路，跑出来的曲线更贴合实际决策场景。适合谁用？一类是做疾控策略模拟的研究生和工程师，另一类是想在MATLAB里快速做多情形对比的仿真开发者。不需要额外装任何工具包，打开主脚本就能看到感染峰值从哪一天开始放缓、消毒强度变化会带来多大差别。2. 从SEIR到SEIRD：模型演进与技术栈选型2.1 标准SEIR丢掉的死亡D，在这里单独成仓先厘清一个概念：SEIRD模型里五个字母分别代表Susceptible（易感）、Exposed（暴露）、Infectious（感染）、Recovered（恢复）、Deceased（死亡）。传统的SEIR把R视为“移除”，康复和死亡混在一起，这会带来两个问题：第一，死亡病例无法参与后续决策策略评估，尤其是在计算床位占用率、病亡率时完全失真；第二，当恢复率高而死亡少时还好，一旦疫情严重，R的累加值会掩盖真实致死趋势。DESSABNeT在这套实现里把D配成独立仓室，并用“消毒网”参数把外部干预动作嵌进微分方程组。简单说，消毒不再只是模型外部的一个注释，而是直接作用在传染系数β上的削减因子。细看代码能发现，传染率的大小除了受基础传播系数影响，还乘了一个随时段改变的消毒强度序列，这样每一次消毒动作都会被模型记录成传染曲线上的一个拐点。2.2 技术栈与仓库结构：两个文件夹分工十分明确DESSABNeT的MATLAB工程结构做得很轻，没有复杂的包管理器依赖，核心就两个目录：sim/和plot/。前者放模型主体函数，后者放可视化脚本，职责清楚，不会出现“改了绘图却要连带翻模型代码”的情况。文件夹大致是这样：DESSABNeT/ ├── sim/ │ ├── main_seird.m # 主运行脚本，定义参数 │ ├── seird_ode.m # SEIRD微分方程核心函数 │ ├── disinfection.m # 消毒策略生成函数 │ └── contact_matrix.m # 人群接触矩阵构建（可选） ├── plot/ │ ├── plot_curves.m # 绘制五仓室曲线 │ ├── plot_peak.m # 感染峰值对比绘图 │ └── plot_heatmap.m # 消毒强度-峰值热力图 └── data/ └── baseline_params.mat # 默认参数文件main_seird.m是入口，所有参数集中在文件头部，适合做批量仿真；seird_ode.m是模型本体，没有任何界面代码，方便二次开发时抽出来单测。这套结构对于研究者来说很友好，因为它遵循的是一条黄金法则：参数、方程、可视化三者分离。如果某天要把模型迁移到Python或R，只需要重写seird_ode.m，其他部分都可以原样保留。3. 核心模块拆解：SEIRD微分主循环与DESSABNeT消毒策略接入3.1 主模型函数：五仓室ODE的底层逻辑打开seird_ode.m，核心代码并不长，但每一行的物理含义都需要说清楚。这个文件是整套模型的运算心脏，MATLAB的ode45会反复调用它来计算每一步的仓室变化量。function dydt = seird_ode(t, y, beta0, sigma, gamma, delta, dis_amp) % y 列向量: [S; E; I; R; D] S = y(1); E = y(2); I = y(3); R = y(4); D = y(5); N = S + E + I + R + D; % 总人口，保持恒定 % 消毒网作用：传染率随时间被削减 beta = beta0 * (1 - dis_amp * exp(-0.1 * t)); dS = -beta * S * I / N; dE = beta * S * I / N - sigma * E; dI = sigma * E - (gamma + delta) * I; dR = gamma * I; dD = delta * I; dydt = [dS; dE; dI; dR; dD]; end这段代码最值得细看的是beta那一行。beta0是不做任何干预时的基础传染率，dis_amp是消毒强度系数，exp(-0.1*t)表示消毒效果随时间的衰减。这种设计隐含了一个公共卫生常识：消毒的边际效应是逐日递减的，不可能连续三十天保持同样的消杀效率。如果直接把消毒设为常量，模拟后期会出现曲线异常偏低、与真实疫情趋势脱节的问题。参数sigma代表潜伏期结束进入感染期的速率，gamma是恢复速率，delta是因病死亡速率。需要特别留意的是，gamma + delta共同决定感染者仓室的迁出速度，两者此消彼长会直接影响感染峰值的形状。如果调高delta而不动gamma，I仓室会更快清空，但D仓室会累积出一个令人警醒的死亡曲线。3.2 消毒策略接入：disinfection函数怎么生成干预计划disinfection.m的作用是生成一段可读的消毒强度时间表，再把这张表喂给主模型。单独拆这个函数出来是有好处的，因为实际场景中消毒策略不是连续函数，而是分段操作——今天全城消杀，明天可能只有重点区域消毒。function dis_table = disinfection(base_strength, schedule, tspan) % schedule: [start_day, end_day, strength_scale] % 返回与 tspan 等长的消毒强度向量 dis_table = zeros(length(tspan), 1); for i = 1:size(schedule, 1) start_day = schedule(i, 1); end_day = schedule(i, 2); scale = schedule(i, 3); idx = tspan >= start_day & tspan <= end_day; dis_table(idx) = base_strength * scale; end end这个函数的输入schedule是一个矩阵，每一行代表一次消毒行动。比如[10, 20, 1.5]表示第10天到第20天之间执行1.5倍强度的消杀作业。输出向量与仿真时间轴等长，方便直接点乘到传染率上。我一般会把schedule从Excel里读进来，这样非编程背景的同事也能直接修改干预计划，不用动代码。二次开发时也可以把一个真实的消毒日志转成这个格式——只要日期和强度对齐就行。contact_matrix.m是另一个可选模块。当总人口N被分成多个子群时（比如学生、在职人员、老年人），每个子群的接触频率不同，用单一β值就不够了。这时传入一个接触矩阵C，传染率会变成beta0 * C * S * I / N，更适合做异质性人群的传播模拟。默认工程只跑单一人群，但结构上已经预留了这个扩展点。3.3 主入口main_seird：串联全流程的关键脚本main_seird.m是直接运行的文件，它的任务是把参数定义、ODE求解、结果返回整合成一条流水线。很多人在这个文件上翻车，因为改完参数立刻运行，结果还在看上一次的，原因是MATLAB的工作区缓存没有清理。clear; clc; % 基础参数 N = 100000; % 总人口 I0 = 10; % 初始感染者 E0 = 0; % 初始暴露者 sigma = 1/5.2; % 潜伏期约5.2天 gamma = 1/10; % 恢复期约10天 delta = 0.015; % 死亡率 beta0 = 0.32; % 基础传染率 % 消毒策略 schedule = [5, 15, 1.2; 20, 30, 0.8]; dis_amp = 0.25; % 时间跨度 tspan = 0:0.5:100; y0 = [N - I0; E0; I0; 0; 0]; % 求解 [t, y] = ode45(@(t, y) seird_ode(t, y, beta0, sigma, gamma, delta, dis_amp), ... tspan, y0); % 留存结果 result.t = t; result.y = y; save('data/run_result.mat', 'result');注意这里tspan用步长0.5而不是默认的整数步长，原因在于ode45是自适应步长算法，输出节点不一定均匀，如果后面要做傅里叶分析或频谱分析，不均匀时间轴会带来很多额外麻烦。统一成0.5天一步，既保持曲线光滑，又方便和真实日报告数据做对齐。参数I0设置为10，意味着第一天就有10个感染者进场。实际使用中应该把I0设成1或者更小，因为社区传播的起点往往很低，设置过大会让峰值前置、显得疫情爆发太快。4. 从克隆到跑通：运行配置与三大参数组的调优方法4.1 环境准备与两种运行方式MATLAB环境基本没有额外要求，只要是R2016b之后的版本都能直接跑。仓库不依赖Simulink，也不依赖任何第三方工具箱，普通的MATLAB基础授权就够。运行方式有两种。第一种是图形界面操作：打开main_seird.m，在编辑器中直接按F5运行，结果保存在data/run_result.mat，然后运行plot/plot_curves.m绘图。第二种是命令行批处理，适合批量跑参数扫描：matlab -batch "run('main_seird.m'); plot_curves"用-batch模式可以避开MATLAB启动GUI的额外开销，连续跑几十组参数时体感快很多。如果机器上装了并行计算工具箱，还能用parfor把多组参数同时送进求解器。需要注意的一点是，-batch模式下工作区不保留变量，必须在脚本内部显式用save把结果写到磁盘，否则跑完什么都没留下。4.2 三大参数组：传播强度、病程时长、消毒强度这组参数可以划分为三组理解，每组对应模型行为的一个方面：表格：核心参数调节方向参数组参数名含义调高后的影响传播强度beta0基础传染率峰值更高、峰值时间提前病程时长sigma,gamma潜伏期/恢复期时长病程越长，感染平台期越久消毒干预dis_amp,schedule消毒强度与时间表峰值下降，曲线出现明显转折beta0的取值可以参考基本再生数R0来推算，近似公式是R0 = beta0 / (gamma + delta)。比如希望模拟R0=3.5的传染病，假设gamma=0.1,delta=0.015，那么beta0应该设在3.5 * 0.115 ≈ 0.4附近。这里特别提醒，算出来的值只是一个参考起点，实际拟合时还需要根据疫情数据反向校准。sigma写的是1/5.2，意思是潜伏期平均5.2天，这对应指数分布的速率参数。如果改成1/3就表示潜伏期缩短到3天，感染峰值会明显左移。调整时注意不要把恢复期和潜伏期设成同一数量级，否则E和I仓室会同时膨胀，曲线之间互相遮掩，后续可视化很难解读。4.3 验证模型行为：曲线不收敛时怎么看模型跑通之后，先别急着做策略对比，第一步是检查五条曲线是否符合常识。正常情况下，S曲线缓慢下降，E和I先后出现一个单峰，R曲线单调上升，D曲线缓慢累积。如果看到I曲线剧烈震荡或者出现负值，说明数值求解出了问题。绝大多数震荡来自ode45对刚性方程的不适应。SEIRD模型在参数极端时（比如beta0很大、delta接近0）会呈现刚性特征，此时换成ode15s求解器往往能立刻解决：[t, y] = ode15s(@(t, y) seird_ode(t, y, beta0, sigma, gamma, delta, dis_amp), ... tspan, y0);另外一个常见验证手段是计算总人口守恒。每行S+E+I+R+D理论上恒等于初始总人口N，如果数值上出现漂移，大概率是因为某个仓室的更新率在积分中丢失。可以在求解结束后跑一行sum(y, 2)查看随时间的变化，偏差超过总人口的1%就得回头查微分方程里的正负号。plot_curves.m会自动把五条曲线画在同一张图上，但默认的纵轴范围可能让D和R曲线因为数值太小而压扁。我会习惯性把D曲线单独画一个对数坐标的子图，否则D的累积趋势很难看清。5. 避坑指南：MATLAB模型调试中的5个实际问题5.1 仿真结果全是NaN现象：跑完main_seird.m后，y矩阵里出现大量NaN。原因：初始感染者I0设为0，导致感染项beta * S * I / N在第一步就是0，所有仓室保持初始值不变，但数值求解器的自适应步长尝试把时间步调大后，某个仓室计算出现除以零的情况，触发NaN连锁扩散。解决：I0至少设为1。还要检查beta0是否过大（比如超过1.5），在ode45默认容差下容易数值爆炸。把beta0回退到0.3附近，NaN现象基本就会消失。5.2 改了参数但输出结果完全不变现象：在编辑器中修改beta0后运行，绘图结果和上一次一模一样。原因：MATLAB的脚本运行默认会沿用上一次工作区中残留的变量。如果某次运行后main_seird.m里的赋值语句被跳过（比如代码报错后继续执行），旧参数就留在工作区里被重复使用。解决：在每个主脚本开头都写clear; clc;，并且使用run('main_seird.m')强制刷新。注意clear会清空所有变量，如果依赖外部导入的Excel数据，需要放在clear之后重新读取。5.3 消毒强度调得很大，曲线却几乎没变化现象：把dis_amp从0.2调到0.9，感染峰值只下降了不到5%。原因：消毒效果被指数衰减项exp(-0.1*t)主导，当t超过30天后衰减系数已经降到0.05以内，消毒强度再大也推不动曲线。解决：缩短消毒作用窗口，把衰减系数从0.1改到0.05，或者直接把指数衰减换成矩形窗函数，这样消毒强度在作用窗口内保持恒定，效果才会真实反映在曲线上。还有一种做法是把dis_amp从标量改成向量，配合schedule使用分段常数控制。5.4 绘图时D曲线被压成一条直线现象：plot_curves.m生成的图中，D和R曲线贴在横轴上，毫无辨识度。原因：D曲线数值比I曲线小两个数量级，线性坐标下几乎不可见。这是SEIRD可视化的通病，不是模型出错。解决：用subplot把D单独画到第二个子图，并给D曲线一个独立的纵轴。也可以改用yyaxis left和yyaxis right把I和D分别放在左右两侧。示例如下：figure; yyaxis left; plot(t, y(:,3), 'LineWidth', 2); ylabel('Infectious'); yyaxis right; plot(t, y(:,5), 'LineWidth', 2); ylabel('Deceased');5.5 多次运行，每次结果都有微小差异现象：同样参数运行两次，峰值时间有时差一天，曲线尾部有细微抖动。原因：如果代码里用了rand或randn生成初始扰动，而随机种子没有固定，结果自然无法复现。这个问题在接触矩阵版中尤其常见。解决：在main_seird.m开头固定随机种子：rng(2025);这一行加上之后，任何随机因素都会被锁死，多轮实验才能公平对比。6. 怎么改出第二个模型：二次开发的三个入手方向第一个入手方向是把单一人群改成多群体模型。现在contact_matrix.m只是预留了接口，没有真正启用。要改成多群体，需要给每个子群单独分配一组beta和gamma。以学生和老年人两个子群为例，学生群体的接触频率高但病程短，老年人群体反之。将接触矩阵写成2x2的数值表，每个元素表示两个群体的接触次数，然后主模型中的传染项从beta * S * I / N改成beta * C(i,j) * S_i * I_j / N，需要同时更新所有仓室的下标引用。完成后可以观察到一个典型现象：感染的峰值不再是单峰，而是出现拖尾的双峰形态。第二个方向是引入时变参数。真实疫情中，防控措施的强度是逐周调整的，而不是一沉不变。可以把beta0从一个标量改成向量，向量的每个元素对应一天的传染率。比如：beta_daily = linspace(0.32, 0.2, 101);用interp1在求解过程中动态插值。这个改法特别适合复盘某次特定疫情，把每日真实传染率曲线直接喂进模型，再对比模型输出的感染人数和真实报告数据，使用最小二乘法拟合出最优参数组合。第三个方向是把输出从曲线改成指标表。现在的输出只是mat文件，真正写报告时还需要峰值日期、累计病例、总死亡数、第90天残余感染数五个指标。我习惯在主脚本末尾加一段汇总代码：peak_I = max(y(:,3)); peak_day = t(y(:,3) == peak_I); total_R = y(end, 4); total_D = y(end, 5); fprintf('峰值感染数: %.0f 第%.1f天\n', peak_I, peak_day); fprintf('累计恢复: %.0f 累计死亡: %.0f\n', total_R, total_D);这段代码让每次仿真的结果一眼就能读取，不用再打开绘图窗口肉眼找峰值。跑完一组参数扫描后，把这些指标收集成一个表格，消毒策略的优化问题就变成了一张可以直观对比的Excel表格。最后说一个我反复踩过的坑。以前我总喜欢在seird_ode.m里直接改参数做试验，结果每次临时改完忘记还原，后续仿真全都被带偏。后来我强制规定：方程文件只保留模型结构，所有参数必须由main_seird.m传入，任何实验性修改一律写在主脚本里，跑完立刻save。从那以后，每次做参数扫描，我都在运行前先clear清场、固定随机种子、校验N的守恒，确认全绿才去看曲线。这套流程虽然多花十秒，但能省掉后面至少两小时的排查时间。希望这篇解析能帮你顺利跑通这套SEIRD模型，并把消毒网策略真正嵌进你的仿真流程里。本文还有配套的精品资源，点击获取
简介：DESSABNeT 是一个基于 MATLAB 编写的开源传染病传播模拟软件包，采用基于代理的模型（ABM）构建人工社会，以真实人口统计数据为输入，可模拟任何城市或区域中传染病的传播动态，最初专门用于澳大利亚城市的 COVID-19 疫情建模。面向流行病学研究者、公共卫生决策者及计算社科领域的 MATLAB 用户，可用于评估社交限制措施、模拟不同疾病爆发场景。资源包共 42 个文件，以 33 个 .m 格式 MATLAB 脚本为核心，覆盖人口设置、社交网络构建、每周时间表配置、疾病传播动态模块及社交限制阶段模拟；另有 7 个 .mat 预生成数据文件、1 个 Excel 数据表和 1 个 README 说明文档，压缩包整体 2.7MB，按模块分目录组织。已有 712 人学习下载。读者可获得完整的 ABM 建模框架代码，包括悉尼人口结构设置、家庭/工作/学校网络构建、感染病例注入与传播计算、R0 矩阵计算、隔离与社交限制模拟等，并提供可直接运行的示例数据，便于移植到其他传染病或城市场景开展仿真实验。1. 这份MATLAB模型，为什么值得你先跑一遍再说做传染病仿真的人多数有一个共识：SEIR已经不够用了。真实世界每一轮疫情都会出现死亡病例，而标准SEIR只把人群分成易感、暴露、感染、恢复四类，死亡和康复只占一个R槽，模拟出来反而比实际结果乐观不少。这套名为DESSABNeT的SEIRD模型MATLAB代码，正好补上了这段缺口——它把死亡（D）单独拆出来，同时把消毒网（Disinfection Network）的影响以参数形式接入传播链路，跑出来的曲线更贴合实际决策场景。适合谁用？一类是做疾控策略模拟的研究生和工程师，另一类是想在MATLAB里快速做多情形对比的仿真开发者。不需要额外装任何工具包，打开主脚本就能看到感染峰值从哪一天开始放缓、消毒强度变化会带来多大差别。2. 从SEIR到SEIRD：模型演进与技术栈选型2.1 标准SEIR丢掉的死亡D，在这里单独成仓先厘清一个概念：SEIRD模型里五个字母分别代表Susceptible（易感）、Exposed（暴露）、Infectious（感染）、Recovered（恢复）、Deceased（死亡）。传统的SEIR把R视为“移除”，康复和死亡混在一起，这会带来两个问题：第一，死亡病例无法参与后续决策策略评估，尤其是在计算床位占用率、病亡率时完全失真；第二，当恢复率高而死亡少时还好，一旦疫情严重，R的累加值会掩盖真实致死趋势。DESSABNeT在这套实现里把D配成独立仓室，并用“消毒网”参数把外部干预动作嵌进微分方程组。简单说，消毒不再只是模型外部的一个注释，而是直接作用在传染系数β上的削减因子。细看代码能发现，传染率的大小除了受基础传播系数影响，还乘了一个随时段改变的消毒强度序列，这样每一次消毒动作都会被模型记录成传染曲线上的一个拐点。2.2 技术栈与仓库结构：两个文件夹分工十分明确DESSABNeT的MATLAB工程结构做得很轻，没有复杂的包管理器依赖，核心就两个目录：sim/和plot/。前者放模型主体函数，后者放可视化脚本，职责清楚，不会出现“改了绘图却要连带翻模型代码”的情况。文件夹大致是这样：DESSABNeT/ ├── sim/ │ ├── main_seird.m # 主运行脚本，定义参数 │ ├── seird_ode.m # SEIRD微分方程核心函数 │ ├── disinfection.m # 消毒策略生成函数 │ └── contact_matrix.m # 人群接触矩阵构建（可选） ├── plot/ │ ├── plot_curves.m # 绘制五仓室曲线 │ ├── plot_peak.m # 感染峰值对比绘图 │ └── plot_heatmap.m # 消毒强度-峰值热力图 └── data/ └── baseline_params.mat # 默认参数文件main_seird.m是入口，所有参数集中在文件头部，适合做批量仿真；seird_ode.m是模型本体，没有任何界面代码，方便二次开发时抽出来单测。这套结构对于研究者来说很友好，因为它遵循的是一条黄金法则：参数、方程、可视化三者分离。如果某天要把模型迁移到Python或R，只需要重写seird_ode.m，其他部分都可以原样保留。3. 核心模块拆解：SEIRD微分主循环与DESSABNeT消毒策略接入3.1 主模型函数：五仓室ODE的底层逻辑打开seird_ode.m，核心代码并不长，但每一行的物理含义都需要说清楚。这个文件是整套模型的运算心脏，MATLAB的ode45会反复调用它来计算每一步的仓室变化量。function dydt = seird_ode(t, y, beta0, sigma, gamma, delta, dis_amp) % y 列向量: [S; E; I; R; D] S = y(1); E = y(2); I = y(3); R = y(4); D = y(5); N = S + E + I + R + D; % 总人口，保持恒定 % 消毒网作用：传染率随时间被削减 beta = beta0 * (1 - dis_amp * exp(-0.1 * t)); dS = -beta * S * I / N; dE = beta * S * I / N - sigma * E; dI = sigma * E - (gamma + delta) * I; dR = gamma * I; dD = delta * I; dydt = [dS; dE; dI; dR; dD]; end这段代码最值得细看的是beta那一行。beta0是不做任何干预时的基础传染率，dis_amp是消毒强度系数，exp(-0.1*t)表示消毒效果随时间的衰减。这种设计隐含了一个公共卫生常识：消毒的边际效应是逐日递减的，不可能连续三十天保持同样的消杀效率。如果直接把消毒设为常量，模拟后期会出现曲线异常偏低、与真实疫情趋势脱节的问题。参数sigma代表潜伏期结束进入感染期的速率，gamma是恢复速率，delta是因病死亡速率。需要特别留意的是，gamma + delta共同决定感染者仓室的迁出速度，两者此消彼长会直接影响感染峰值的形状。如果调高delta而不动gamma，I仓室会更快清空，但D仓室会累积出一个令人警醒的死亡曲线。3.2 消毒策略接入：disinfection函数怎么生成干预计划disinfection.m的作用是生成一段可读的消毒强度时间表，再把这张表喂给主模型。单独拆这个函数出来是有好处的，因为实际场景中消毒策略不是连续函数，而是分段操作——今天全城消杀，明天可能只有重点区域消毒。function dis_table = disinfection(base_strength, schedule, tspan) % schedule: [start_day, end_day, strength_scale] % 返回与 tspan 等长的消毒强度向量 dis_table = zeros(length(tspan), 1); for i = 1:size(schedule, 1) start_day = schedule(i, 1); end_day = schedule(i, 2); scale = schedule(i, 3); idx = tspan >= start_day & tspan <= end_day; dis_table(idx) = base_strength * scale; end end这个函数的输入schedule是一个矩阵，每一行代表一次消毒行动。比如[10, 20, 1.5]表示第10天到第20天之间执行1.5倍强度的消杀作业。输出向量与仿真时间轴等长，方便直接点乘到传染率上。我一般会把schedule从Excel里读进来，这样非编程背景的同事也能直接修改干预计划，不用动代码。二次开发时也可以把一个真实的消毒日志转成这个格式——只要日期和强度对齐就行。contact_matrix.m是另一个可选模块。当总人口N被分成多个子群时（比如学生、在职人员、老年人），每个子群的接触频率不同，用单一β值就不够了。这时传入一个接触矩阵C，传染率会变成beta0 * C * S * I / N，更适合做异质性人群的传播模拟。默认工程只跑单一人群，但结构上已经预留了这个扩展点。3.3 主入口main_seird：串联全流程的关键脚本main_seird.m是直接运行的文件，它的任务是把参数定义、ODE求解、结果返回整合成一条流水线。很多人在这个文件上翻车，因为改完参数立刻运行，结果还在看上一次的，原因是MATLAB的工作区缓存没有清理。clear; clc; % 基础参数 N = 100000; % 总人口 I0 = 10; % 初始感染者 E0 = 0; % 初始暴露者 sigma = 1/5.2; % 潜伏期约5.2天 gamma = 1/10; % 恢复期约10天 delta = 0.015; % 死亡率 beta0 = 0.32; % 基础传染率 % 消毒策略 schedule = [5, 15, 1.2; 20, 30, 0.8]; dis_amp = 0.25; % 时间跨度 tspan = 0:0.5:100; y0 = [N - I0; E0; I0; 0; 0]; % 求解 [t, y] = ode45(@(t, y) seird_ode(t, y, beta0, sigma, gamma, delta, dis_amp), ... tspan, y0); % 留存结果 result.t = t; result.y = y; save('data/run_result.mat', 'result');注意这里tspan用步长0.5而不是默认的整数步长，原因在于ode45是自适应步长算法，输出节点不一定均匀，如果后面要做傅里叶分析或频谱分析，不均匀时间轴会带来很多额外麻烦。统一成0.5天一步，既保持曲线光滑，又方便和真实日报告数据做对齐。参数I0设置为10，意味着第一天就有10个感染者进场。实际使用中应该把I0设成1或者更小，因为社区传播的起点往往很低，设置过大会让峰值前置、显得疫情爆发太快。4. 从克隆到跑通：运行配置与三大参数组的调优方法4.1 环境准备与两种运行方式MATLAB环境基本没有额外要求，只要是R2016b之后的版本都能直接跑。仓库不依赖Simulink，也不依赖任何第三方工具箱，普通的MATLAB基础授权就够。运行方式有两种。第一种是图形界面操作：打开main_seird.m，在编辑器中直接按F5运行，结果保存在data/run_result.mat，然后运行plot/plot_curves.m绘图。第二种是命令行批处理，适合批量跑参数扫描：matlab -batch "run('main_seird.m'); plot_curves"用-batch模式可以避开MATLAB启动GUI的额外开销，连续跑几十组参数时体感快很多。如果机器上装了并行计算工具箱，还能用parfor把多组参数同时送进求解器。需要注意的一点是，-batch模式下工作区不保留变量，必须在脚本内部显式用save把结果写到磁盘，否则跑完什么都没留下。4.2 三大参数组：传播强度、病程时长、消毒强度这组参数可以划分为三组理解，每组对应模型行为的一个方面：表格：核心参数调节方向参数组参数名含义调高后的影响传播强度beta0基础传染率峰值更高、峰值时间提前病程时长sigma,gamma潜伏期/恢复期时长病程越长，感染平台期越久消毒干预dis_amp,schedule消毒强度与时间表峰值下降，曲线出现明显转折beta0的取值可以参考基本再生数R0来推算，近似公式是R0 = beta0 / (gamma + delta)。比如希望模拟R0=3.5的传染病，假设gamma=0.1,delta=0.015，那么beta0应该设在3.5 * 0.115 ≈ 0.4附近。这里特别提醒，算出来的值只是一个参考起点，实际拟合时还需要根据疫情数据反向校准。sigma写的是1/5.2，意思是潜伏期平均5.2天，这对应指数分布的速率参数。如果改成1/3就表示潜伏期缩短到3天，感染峰值会明显左移。调整时注意不要把恢复期和潜伏期设成同一数量级，否则E和I仓室会同时膨胀，曲线之间互相遮掩，后续可视化很难解读。4.3 验证模型行为：曲线不收敛时怎么看模型跑通之后，先别急着做策略对比，第一步是检查五条曲线是否符合常识。正常情况下，S曲线缓慢下降，E和I先后出现一个单峰，R曲线单调上升，D曲线缓慢累积。如果看到I曲线剧烈震荡或者出现负值，说明数值求解出了问题。绝大多数震荡来自ode45对刚性方程的不适应。SEIRD模型在参数极端时（比如beta0很大、delta接近0）会呈现刚性特征，此时换成ode15s求解器往往能立刻解决：[t, y] = ode15s(@(t, y) seird_ode(t, y, beta0, sigma, gamma, delta, dis_amp), ... tspan, y0);另外一个常见验证手段是计算总人口守恒。每行S+E+I+R+D理论上恒等于初始总人口N，如果数值上出现漂移，大概率是因为某个仓室的更新率在积分中丢失。可以在求解结束后跑一行sum(y, 2)查看随时间的变化，偏差超过总人口的1%就得回头查微分方程里的正负号。plot_curves.m会自动把五条曲线画在同一张图上，但默认的纵轴范围可能让D和R曲线因为数值太小而压扁。我会习惯性把D曲线单独画一个对数坐标的子图，否则D的累积趋势很难看清。5. 避坑指南：MATLAB模型调试中的5个实际问题5.1 仿真结果全是NaN现象：跑完main_seird.m后，y矩阵里出现大量NaN。原因：初始感染者I0设为0，导致感染项beta * S * I / N在第一步就是0，所有仓室保持初始值不变，但数值求解器的自适应步长尝试把时间步调大后，某个仓室计算出现除以零的情况，触发NaN连锁扩散。解决：I0至少设为1。还要检查beta0是否过大（比如超过1.5），在ode45默认容差下容易数值爆炸。把beta0回退到0.3附近，NaN现象基本就会消失。5.2 改了参数但输出结果完全不变现象：在编辑器中修改beta0后运行，绘图结果和上一次一模一样。原因：MATLAB的脚本运行默认会沿用上一次工作区中残留的变量。如果某次运行后main_seird.m里的赋值语句被跳过（比如代码报错后继续执行），旧参数就留在工作区里被重复使用。解决：在每个主脚本开头都写clear; clc;，并且使用run('main_seird.m')强制刷新。注意clear会清空所有变量，如果依赖外部导入的Excel数据，需要放在clear之后重新读取。5.3 消毒强度调得很大，曲线却几乎没变化现象：把dis_amp从0.2调到0.9，感染峰值只下降了不到5%。原因：消毒效果被指数衰减项exp(-0.1*t)主导，当t超过30天后衰减系数已经降到0.05以内，消毒强度再大也推不动曲线。解决：缩短消毒作用窗口，把衰减系数从0.1改到0.05，或者直接把指数衰减换成矩形窗函数，这样消毒强度在作用窗口内保持恒定，效果才会真实反映在曲线上。还有一种做法是把dis_amp从标量改成向量，配合schedule使用分段常数控制。5.4 绘图时D曲线被压成一条直线现象：plot_curves.m生成的图中，D和R曲线贴在横轴上，毫无辨识度。原因：D曲线数值比I曲线小两个数量级，线性坐标下几乎不可见。这是SEIRD可视化的通病，不是模型出错。解决：用subplot把D单独画到第二个子图，并给D曲线一个独立的纵轴。也可以改用yyaxis left和yyaxis right把I和D分别放在左右两侧。示例如下：figure; yyaxis left; plot(t, y(:,3), 'LineWidth', 2); ylabel('Infectious'); yyaxis right; plot(t, y(:,5), 'LineWidth', 2); ylabel('Deceased');5.5 多次运行，每次结果都有微小差异现象：同样参数运行两次，峰值时间有时差一天，曲线尾部有细微抖动。原因：如果代码里用了rand或randn生成初始扰动，而随机种子没有固定，结果自然无法复现。这个问题在接触矩阵版中尤其常见。解决：在main_seird.m开头固定随机种子：rng(2025);这一行加上之后，任何随机因素都会被锁死，多轮实验才能公平对比。6. 怎么改出第二个模型：二次开发的三个入手方向第一个入手方向是把单一人群改成多群体模型。现在contact_matrix.m只是预留了接口，没有真正启用。要改成多群体，需要给每个子群单独分配一组beta和gamma。以学生和老年人两个子群为例，学生群体的接触频率高但病程短，老年人群体反之。将接触矩阵写成2x2的数值表，每个元素表示两个群体的接触次数，然后主模型中的传染项从beta * S * I / N改成beta * C(i,j) * S_i * I_j / N，需要同时更新所有仓室的下标引用。完成后可以观察到一个典型现象：感染的峰值不再是单峰，而是出现拖尾的双峰形态。第二个方向是引入时变参数。真实疫情中，防控措施的强度是逐周调整的，而不是一沉不变。可以把beta0从一个标量改成向量，向量的每个元素对应一天的传染率。比如：beta_daily = linspace(0.32, 0.2, 101);用interp1在求解过程中动态插值。这个改法特别适合复盘某次特定疫情，把每日真实传染率曲线直接喂进模型，再对比模型输出的感染人数和真实报告数据，使用最小二乘法拟合出最优参数组合。第三个方向是把输出从曲线改成指标表。现在的输出只是mat文件，真正写报告时还需要峰值日期、累计病例、总死亡数、第90天残余感染数五个指标。我习惯在主脚本末尾加一段汇总代码：peak_I = max(y(:,3)); peak_day = t(y(:,3) == peak_I); total_R = y(end, 4); total_D = y(end, 5); fprintf('峰值感染数: %.0f 第%.1f天\n', peak_I, peak_day); fprintf('累计恢复: %.0f 累计死亡: %.0f\n', total_R, total_D);这段代码让每次仿真的结果一眼就能读取，不用再打开绘图窗口肉眼找峰值。跑完一组参数扫描后，把这些指标收集成一个表格，消毒策略的优化问题就变成了一张可以直观对比的Excel表格。最后说一个我反复踩过的坑。以前我总喜欢在seird_ode.m里直接改参数做试验，结果每次临时改完忘记还原，后续仿真全都被带偏。后来我强制规定：方程文件只保留模型结构，所有参数必须由main_seird.m传入，任何实验性修改一律写在主脚本里，跑完立刻save。从那以后，每次做参数扫描，我都在运行前先clear清场、固定随机种子、校验N的守恒，确认全绿才去看曲线。这套流程虽然多花十秒，但能省掉后面至少两小时的排查时间。希望这篇解析能帮你顺利跑通这套SEIRD模型，并把消毒网策略真正嵌进你的仿真流程里。本文还有配套的精品资源，点击获取