简介这是面向流行病学建模与仿真研究的MATLAB开源软件包DESSABNeT采用基于代理的模拟方法构建由个体代理组成的人工社会可对新冠、流感等经人际接触传播的传染病在特定城市或区域内的扩散过程进行推演。程序会结合真实人口统计、年龄分层、家庭及工作结构、每周社交频度等数据并为每个代理安排动态时间表以模拟隔离、封锁等社交限制对疫情的影响。资源共包含42个文件其中33个m文件为完整MATLAB源码覆盖人口数据生成、各类社交关系容器建立、动态感染模型、基本再生数计算、病例输入与隔离代码等模块另含7个mat预置仿真数据集、1个xlsx真实病例数据表及1个md使用说明整体大小仅2.7MB。目前已有712人浏览学习。借助这套代码读者能够复现基于悉尼实际人口数据的COVID-19传播仿真实验也可替换数据与参数开展其他城市或假想场景的模拟研究适合作为传染病建模课程设计、科研预研和系统开发参考。1. DESSABNeT把死亡通道从SEIR里拆出来的MATLAB开源仿真做传播病防控复盘时团队拿到过去60天的确诊与死亡序列想估算隔离措施的边际效果。手头现成的SEIR代码把康复和死亡合并在一个R里医疗负荷预测明显偏乐观决策层根本不敢用。DESSABNeT消毒网是一套系统开源的MATLAB项目就是冲这个痛点来的它在经典SEIR模型的基础上加入独立的死亡仓室D同时保留暴露仓室E让潜伏期、感染期、病死率各自拥有独立的可调参数。名字里的“消毒”也贴切——数据先进清洗流程再进模型网络推演。适合做曲线拟合、干预模拟和教学演示的开发同学不依赖额外工具箱从拿到源码到跑出第一张图半小时足够。2. 架构视角SEIRD微分方程与DESSABNeT的工程取舍2.1 SEIRD模型的数学骨架与每个仓室的物理含义SEIRD与经典SEIR最大的差别是把R仓室拆成了R康复和D死亡。这句话看起来只是多一个变量但做资源规划时意义完全不同康复者对应出院人群死亡者对应床位与防护资源的刚需两者混在一起模拟出的ICU压力会低一大截。模型总人口N S E I R D固定。状态随时间变化由下面这组常微分方程控制% core/seird_ode.m % 输入 x [S; E; I; R; D]p 为参数结构体t 为当前时刻 function dxdt seird_ode(t, x, p) S x(1); E x(2); I x(3); R x(4); D x(5); N p.N; infection p.beta * S * I / N; % 单位时间新增感染人数 dxdt zeros(5, 1); dxdt(1) -infection; % 易感者减少 dxdt(2) infection - p.sigma * E; % 潜伏者到期转感染 dxdt(3) p.sigma * E - p.gamma * I - p.mu * I; % 感染者分流 dxdt(4) p.gamma * I; % 康复累积单调不减 dxdt(5) p.mu * I; % 死亡累积单调不减 end感染项写成βSI/N而不是βSI原因是让传播率不随人口规模放大N越大每个感染者平均有效接触带来的新感染按比例摊薄。σ 1/潜伏期γ 1/感染期μ在这套代码里直接是单位时间死亡速率。demo默认参数见下表参数含义demo默认值单位调参方向beta有效接触率0.301/天增大则峰值提前、峰值变高sigma潜伏期倒数0.201/天增大则E仓室快速耗尽gamma感染期倒数0.101/天增大则康复速度加快mu日死亡速率0.021/天增大则D曲线更陡注意守恒关系SEIRD理论上不变。但ode45是数值积分浮点误差会让总量缓慢漂移连续模拟120天后可能偏零点几个百分点。正确做法是在后处理阶段对曲线做一次性归一化而不是在方程里强行把总量拉回N——改导数方向会破坏求解器的误差控制这一步是常见翻车点。2.2 目录结构与分层设计config / core / postprocess 各管一摊DESSABNeT按“参数—方程—求解—后处理”四段拆开目录结构如下DESSABNeT/ ├── main.m ├── config/ │ └── default_params.m ├── core/ │ ├── seird_ode.m │ ├── solve_seird.m │ └── fit_params.m ├── postprocess/ │ ├── plot_results.m │ └── export_csv.m ├── data/ │ └── demo_epidemic.csv └── tests/ └── run_tests.mconfig目录只负责返回一个结构体params包含N、beta、sigma、gamma、mu、初值I0和模拟天数core层的三个文件互不依赖main里的变量solve_seird只接收params和初值返回时间向量t与状态矩阵Xpostprocess拿到(t, X)后出图和导出CSV。这个分层最实用的收益是换数据不用动core换干预策略不用动postprocess后续加“隔离强度”、加“无症状仓室”都在core里加方程其它层基本不动。数据流一句话main读config → solve_seird调ode45 → seird_ode求导 → postprocess收图收表。所有中间状态都通过结构体传递避免全局变量污染工作区。为什么坚持用结构体而不是一堆全局变量因为MATLAB的ODE函数签名只要求f(t, x)多余参数用匿名函数包一层即可(t,x) seird_ode(t,x,p)闭包携带参数清晰且不会污染命令窗口。代价是参数结构体字段名一旦改所有引用位置都要同步所以default_params.m里每个字段我都会写注释改字段名先全局搜索。为什么不建议手写固定步长RK4ode45内置自适应步长对SEIRD这类不算刚性的系统精度和速度平衡最好固定步长RK4容易在峰值附近踩台阶出现锯齿振荡。这些都是代码评审时反复验证过的结论。3. 核心模块拆解从main到plot_results的一次完整数据流3.1 入口main.m参数装载与初值设定的约定main.m做三件事把路径挂进MATLAB、读默认参数、构造初值并调用求解。代码非常短% main.m clear; close all; clc; addpath(./core, ./config, ./postprocess); p default_params(); % 返回结构体模型参数集中在这 x0 [p.N - p.init_I, 0, p.init_I, 0, 0]; % S,E,I,R,D [t, X] solve_seird(p, x0, [0, p.days]); plot_results(t, X, p);初值x0的顺序必须与seird_ode里的x(1)~x(5)严格一致这是最容易踩的坑新手常把R和D的位置写反导致死亡曲线和康复曲线互串。demo默认init_I 1也就是第一天有1个感染者其余全部易感。如果改成init_I 0模型会永远停在无疫情状态因为感染项是S×I的乘积没有传染源就没有后续这不是bug是模型本身的机制约束。3.2 求解核心solve_seird.mode45与容差配置的意义solve_seird的代码只有几行但容差配置决定了曲线是否可信% core/solve_seird.m function [t, X] solve_seird(p, x0, tSpan) opts odeset(RelTol, 1e-6, AbsTol, 1e-7, ... NonNegative, 1:5); [t, X] ode45((t, x) seird_ode(t, x, p), tSpan, x0, opts); endRelTol和AbsTol控制局部误差。人口基数到一万量级时AbsTol设到1e-7能避免感染人数在峰值之后被计算误差抹成微小数甚至负数。NonNegative 1:5是五条仓室的“保险丝”强制所有状态不允许变负。但要注意NonNegative不是万能的它本质是让求解器在状态逼近零时缩小步长代价是计算变慢如果发现模拟半天跑不完优先检查是不是AbsTol设得太狠了。3.3 参数反演fit_params.m由观测序列反推传播参数拿到真实观测数据后要回答的问题是给定I和D的累计曲线能反推beta和mu吗DESSABNeT在fit_params里用的是最小二乘加fminsearch% core/fit_params.m function p fit_params(obs_t, obs_I, p0) theta0 [p0.beta, p0.sigma, p0.gamma, p0.mu]; opt optimset(Display, off); theta_opt fminsearch((th) residual(th, obs_t, obs_I, p0), theta0, opt); p p0; p.beta theta_opt(1); p.sigma theta_opt(2); p.gamma theta_opt(3); p.mu theta_opt(4); end function r residual(theta, obs_t, obs_I, p) p.beta theta(1); p.sigma theta(2); p.gamma theta(3); p.mu theta(4); [~, X] solve_seird(p, [p.N-1, 0, 1, 0, 0], [0, max(obs_t)]); I_sim interp1(t, X(:, 3), obs_t); % 插值到观测时间点 r sum((I_sim - obs_I).^2); end注意内层求解始终用I0 1的初值与外层观测的day0对齐。如果观测序列不是从第0天开始必须做前置对齐否则拟合出来的参数会把相位差硬塞给beta结果从物理含义上就没法解释。fminsearch是MATLAB基础函数不依赖优化工具箱代价是对初值敏感我一般会把theta0做2到3个随机扰动起点取残差最小的那组作为最终结果。3.4 后处理plot_results.m出图与导出CSVplot_results不仅画五条状态曲线还额外画“日新增感染”曲线作用是直接看峰值出现在第几天% postprocess/plot_results.m function plot_results(t, X, p) figure(Color, w); plot(t, X(:,1), b-, t, X(:,2), g-, t, X(:,3), r-, ... t, X(:,4), m-, t, X(:,5), k-, LineWidth, 1.5); legend(S, E, I, R, D); xlabel(天数); ylabel(人数); grid on; dI [0; diff(X(:,3))]; % 逐日新增感染 figure; bar(t, dI); title(日新增感染); end这里的dI直接用diff求相邻时刻的I差在平台期会因为数值锯齿产生负的小柱子画出来有点难看。我习惯改成滑动平均再做图这个细节在第5章会说。4. 运行与参数调优让demo数据变成可用的决策曲线4.1 环境准备与五分钟跑通运行门槛很低装好MATLABR2016b以后的版本都可以我在R2020b上实测过就能直接跑。把项目文件放到本地后进入DESSABNeT根目录在系统命令行执行cd DESSABNeT matlab -batch main如果MATLAB版本较老、不支持-batch参数直接在MATLAB命令窗口里输入main再回车效果一样。跑完后会弹出两个图窗同时在工作目录下生成一个CSV结果文件。第一次运行建议不要改任何参数用demo默认值验证环境。验证标准感染峰值大约出现在第30到40天之间D曲线持续上升但斜率在后期回落。环境验证通过后再进入调参环节。4.2 调参数看曲线beta与mu的肉眼敏感性试验调参与其说是“优化”不如先做“肉眼敏感性分析”。改config/default_params.m里的beta就能直观看到传播速度对曲线的影响% config/default_params.m p.N 10000; % 人口基数 p.beta 0.30; % 有效接触率 p.sigma 0.20; % 1/潜伏期 p.gamma 0.10; % 1/感染期 p.mu 0.02; % 单位时间死亡速率 p.init_I 1; % 初始感染人数 p.days 120; % 模拟天数把beta从0.3改成0.5感染项变大峰值明显提前、峰值变高改成0.1曲线平缓很多峰值后移120天可能还没到峰顶。这套试验跑一次只要一两分钟适合快速理解模型行为。常见误用是同时改好几个参数最后看不出曲线变化到底归因于谁。我一般一次只动一个参数记录三项指标峰值高度、峰值出现日、第120天时的D累计值。4.3 多情景对比把一次性模拟改成批量跑法真实决策场景通常不是跑一条曲线而是对比“不干预”“中等隔离”“强隔离”几种情景。DESSABNeT没有专门做批量入口但利用solve_seird的纯函数特性写个循环就能批量跑% scripts/run_scenarios.m betas [0.3, 0.2, 0.1]; % 三个干预强度 colors {r, b, k}; figure; hold on; for i 1:length(betas) p default_params(); p.beta betas(i); [t, X] solve_seird(p, [p.N-1, 0, 1, 0, 0], [0, p.days]); plot(t, X(:,3), colors{i}, LineWidth, 1.5); end legend(不干预, 中等隔离, 强隔离);这里有个MATLAB细节结构体是值类型copy-on-write机制保证每次循环里p default_params()拿到的都是全新实例修改beta不会影响上一轮结果。如果把结构体换成全局变量这个批量脚本就会互相污染跑出来的三条曲线全挤在一起。5. 常见问题与避坑指南5.1 曲线出现负值或NaN现象跑着跑着I(t)变成负数再往后输出全是NaN。 原因ode45默认容差只到1e-3级别在感染峰值附近步长跨过拐点误差把接近零的I压成负值负值再通过感染项扩散到其它仓室最终数值积分崩溃。 解决给odeset加上NonNegative 1:5并收紧AbsTol到1e-7如果参数组合下系统变得刚硬把求解器换成ode15s再试一次。记住NonNegative不是开关而是让求解器在边界自动缩小步长计算时间会翻倍但曲线可信度大幅提升。5.2 拟合发散或震荡现象fit_params跑出来的beta接近0.8重建曲线明显超调。 原因四个参数同时反演sigma、gamma、mu之间存在强相关目标函数面上有一整条“等效参数谷”很多组参数都能产生相似的I曲线。说白了多参数同时拟合在小样本下多少有点玄学成分。 解决改成两阶段拟合先固定潜伏期和感染期不参与估计只反演beta和mu得到稳定结果后再放宽到更多参数做微调。fminsearch对初值敏感每次从不同起点跑取残差最小的一组。5.3 R与D曲线在图上互串现象死亡曲线和康复曲线形状完全一样只是高度不同怎么看都不像真实数据。 原因初值x0 [N-1, 0, 1, 0, 0]中R和D的位置写反或者seird_ode里dxdt(4)和dxdt(5)的顺序与初值不一致。 解决先确认main.m里的初值顺序是S, E, I, R, D再看seird_ode的dxdt(4)gammaI、dxdt(5)muI两处严格对应。这个坑从头到尾检查一遍只要一分钟却能省掉后面一整轮重跑属于最值得记住的检查项。5.4 观测数据起始时刻对不齐现象用真实的累计I序列直接拟合残差很大且反推参数毫无意义。 原因模型day0是“出现第1个感染者”的时刻而采集序列通常从首例公开通报后才开始记录起点对不齐时所有相位差都被参数吸收。 解决拟合前先做对齐把观测序列前推到虚拟的I1时刻或者把模型初值时间轴平移让两者在同一个日历单位上。没有后悔药开局不看对齐后面拟合全是白跑。5.5 日新增柱状图出现负数现象bar(t, dI)里冒出一根向下的柱子。 原因数值解在平台期有微小锯齿diff相邻时刻一减负的“新增”就出现了。 解决不要直接用diff(X(:,3))改成滑动平均或用累计区间差值比如每7天求和做日平均画出来就是平滑的周均新增曲线。6. 二次开发把DESSABNeT接入自己的疫情数据6.1 数据清洗与接口改造想换自己的数据准备一个两列CSVday、new_I每日新增感染。清洗时注意三点累计曲线必须单调不减缺失日期用线性插值补齐如果采集序列是从某天开始的先用前置对齐把起点处理好。把数据读进来后转成累计曲线% 数据清洗示例 raw readtable(my_data.csv); day raw.day; dailyI raw.new_I; cumI cumsum(dailyI); % 每日新增转累计 obs_t day; obs_I cumI;6.2 加入时变β(t)干预项与模型验证最常见的一次开发是把固定β改成时变β(t)比如模拟第30天开始加强隔离β从0.3降到0.15。只需要改seird_ode的感染项function beta_t switch_beta(t, p) if t p.intervention_day beta_t p.beta * p.intervention_strength; else beta_t p.beta; end end % 方程里调用 infection beta_t * S * I / N验证方法先用fit_params拟合无干预段的参数设置intervention_day后再拟合干预段比较残差平方和与简洁程度指标确认加入干预项后拟合效果显著提升而不是纯粹多了一个自由度。从那以后我每次给模型喂数据前都先画一遍累计曲线做单调性检查改完参数强制走一遍残差核查省掉了无数次“看起来很正常、一上线就翻车”的深夜排查。希望帮到你。本文还有配套的精品资源点击获取