做传染病拟合的时候最烦的事情就是SEIR模型里那几个参数传播率β、潜伏期转阳率σ、康复率γ随便哪个手调都能折腾一下午。调到后来曲线还是对不齐尤其是峰值位置和达峰时间差个两三天后面的预测就完全没法看。后来我把哈里斯鹰优化算法Harris Hawks Optimization, HHO接到SEIR模型上做参数自动搜索用真实数据跑了几轮拟合效果比人工试凑稳定得多。这篇就把HHO-SEIR的实现思路、Matlab代码结构和踩过的坑一起整理出来给正在折腾传染病模型参数拟合的朋友一个可以直接抄作业的参考。这套方案的核心思路其实很简单把“找参数”变成“找最优解”SEIR模型负责生成预测曲线HHO算法负责搜索使得预测曲线与真实数据误差最小的参数组合。整个过程不需要梯度信息不需要人工反复试凑属于典型的黑箱优化。对于传染病的传播动力学这类非线性模型来说这种方案非常实用——你不需要对模型本身做任何线性化近似也不用担心参数之间存在耦合导致手调顾此失彼。文章后面会把SEIR模型、HHO算法原理、代码实现逐层拆开讲适合有一定Matlab基础、正在做传染病建模或者需要做参数拟合的研究生和工程师参考。1. 为什么SEIR模型的参数这么难拟合1.1 SEIR模型的基本结构SEIR模型把人群分成四个仓室易感者Susceptible、潜伏者Exposed、感染者Infectious和康复者Recovered。与更简单的SIR模型相比SEIR多了一个潜伏期状态E——个体被感染之后并不会立刻具备传染力而是经过一段潜伏期后才进入感染状态。这个设计对新冠这类具有明显潜伏期的传染病尤其重要。四个仓室之间的流转关系可以用微分方程组描述dS/dt -β·S·I/N dE/dt β·S·I/N - σ·E dI/dt σ·E - γ·I dR/dt γ·I其中β是有效接触感染率σ是潜伏期转阳率潜伏期的倒数γ是康复率感染周期的倒数N是总人口数。核心参数就是这三个再加上初始状态S0、E0、I0、R0就可以完整描述一条流行曲线。1.2 手调参数的痛点刚接触SEIR的人最喜欢干的事就是打开Excel或者Matlab写个循环一组一组试β和γ。这种做法的核心问题在于三个参数之间不是独立的。β增大会让感染高峰前移、峰值变高γ减小则会让感染持续期拉长、峰值后移σ会影响曲线的前段爬坡速率。调整其中一个参数往往需要连带调整另外两个才能保持拟合效果。参数之间有耦合关系时人工试凑就变成了一个高维搜索问题靠眼睛判断误差、靠手动微调参数效率极低。另一个痛点是缺乏统一的误差标准。手调的时候这次看峰值对不对下次看总感染人数对不对再下次看曲线末端的dying tail对不对——评价标准一直在变参数自然很难调到一个真正均衡的位置。而优化算法解决的就是这个问题先把“拟合得好不好”量化成一个目标函数然后让算法在这个目标函数的指引下自动搜索。1.3 传统参数拟合方法的局限传统做法里最小二乘曲线拟合和梯度下降是比较常见的两种。最小二乘的问题在于SEIR模型输出不是参数的简单线性组合用lsqcurvefit这类工具时非常依赖初始值初值给得不好直接收敛到局部最优。梯度下降则需要模型的导数信息对于微分方程组来说梯度需要通过伴随方法或数值差分计算实现复杂度立刻上升一个量级而且SEIR这种多峰曲线一旦陷入局部谷底就很难跳出来。这一步让我下定决心换元启发式算法HHO不依赖梯度对目标函数的形态没有要求全局搜索能力强参数设置也不算复杂。用一句话说就是把麻烦事交给随机搜索策略把判断力留给目标函数。2. 哈里斯鹰算法HHO的核心机制2.1 算法灵感来源哈里斯鹰算法的灵感来自哈里斯鹰在沙漠中的群体捕猎行为。这种鹰有个非常有意思的特点它们会以群体方式围猎兔子多只鹰轮流从不同方向惊扰猎物让兔子不断改变逃跑方向、消耗体力最后在合适的时机俯冲抓捕。整个过程从宏观上看就是一个搜索与开发的交替策略和优化算法的核心矛盾——既要广撒网探索可行域又要在有希望的区域重点挖——有着天然的对应关系。算法里把每只鹰的位置对应到一组待优化参数把兔子猎物的位置对应到当前找到的最优解把围猎过程的每一次位置更新对应到一次迭代搜索。这个比喻虽然听起来偏生物学但数学建模之后就是一个标准的启发式优化器。2.2 从探索到开发的转换机制HHO的迭代逻辑大致分为三个阶段探索阶段、探索到开发的转换阶段、开发阶段。转换阶段是关键它由逃逸能量E控制计算公式是E 2·E0·(1 - t/T)其中E0在每次迭代开始时取(-1,1)之间的随机值t是当前迭代次数T是最大迭代次数。随着迭代推进E的绝对值整体从接近2逐渐衰减到接近0算法也从探索行为过渡到开发行为。探索阶段发生在|E|≥1时。这时鹰的位置更新有两种策略一是随机落在种群中个体的附近二是向当前最优个体靠拢但加入随机偏移。这个机制保证了前期搜索范围足够大不容易错过全局最优解。当|E|1时进入开发阶段这一步又分成软围攻、硬围攻、带渐进式快速俯冲的软围攻和硬围攻四种。核心区别是兔子的逃逸能量大小用随机数r和E判断和鹰的俯冲策略不同。渐进式快速俯冲用到了Levy飞行可以让鹰在局部区域外快速跳跃探索避免在一个位置反复试探的效率浪费。2.3 为什么HHO适合SEIR参数优化SEIR参数优化的目标函数是由微分方程数值解生成的形状复杂而且很可能存在多个局部最优谷底。HHO对这个问题的适配性来自两点第一探索阶段维持了种群多样性即便初期收敛到了某个局部区域Levy飞行和随机机制依然能带来跳出局部的能力第二它属于无导数优化目标函数只需要能算出一个数值就行不管内部是微分方程还是查表对优化器来说都一视同仁。实测试验中我对同一个数据集分别跑了遗传算法GA、粒子群PSO和HHO三类算法种群规模和迭代次数相同的情况下HHO的收敛速度和最终拟合误差综合表现最好。当然这不能算严格对比实验但至少说明HHO在这个问题上是靠谱的。3. 基于HHO-SEIR的整体实现方案3.1 整体流程设计整套代码的逻辑划分为四个模块数据准备模块、SEIR模型求解模块、目标函数模块、HHO优化器模块。串联流程是先读入每日新增确诊数据将其处理成累计感染人数或每日新增序列然后定义SEIR微分方程函数交给ode45求解再把求解结果与真实数据对比计算误差作为HHO的适应度函数最后由HHO迭代搜索参数组合输出最优参数和相关评估指标。% 主流程示意 % data 真实每日新增或累计数据 % param_bounds [beta下限 beta上限; sigma下限 sigma上限; gamma下限 gamma上限] % MaxIter 200; NPop 20; [best_param, best_fitness, convergence] HHO_SEIR(seir_objective, param_bounds, NPop, MaxIter, data);3.2 SEIR模型求解模块模型求解我用的是Matlab内置的ode45自适应步长的RK45算法处理SEIR这类非刚性问题足够了。定义方程组时需要注意输入格式function dydt seir_ode(t, y, beta, sigma, gamma, N) S y(1); E y(2); I y(3); R y(4); dSdt -beta * S * I / N; dEdt beta * S * I / N - sigma * E; dIdt sigma * E - gamma * I; dRdt gamma * I; dydt [dSdt; dEdt; dIdt; dRdt]; end求解时用odeset设置非负约束避免数值振荡导致S或I出现负值。这点在实际运行中很容易被忽略尤其是初始感染人数很小或参数范围边界比较极端时负值会让目标函数突然跳到一个非常大的误差值。3.3 目标函数的设计策略目标函数是整个优化过程的核心它决定了算法往哪个方向走。我试过两种方案用每日新增感染人数作为拟合目标以及用累计感染人数作为拟合目标。两者的差异很大每日新增对模型的峰值位置和宽度非常敏感拟合难度高但信息量大累计数据更平滑拟合相对容易但容易掩盖内部的动态特征。实际采用哪种取决于手头有什么数据。如果只有累计确诊序列就用累计数据做拟合如果有每日新增序列建议优先用每日新增做目标函数。目标函数我一般写成归一化的均方根误差RMSE并对数量级差异做了缩放处理function err seir_objective(param, data, N, tspan) beta param(1); sigma param(2); gamma param(3); y0 [N - data(1) - 1; 0; data(1); 0]; [~, Y] ode45((t,y) seir_ode(t, y, beta, sigma, gamma, N), tspan, y0, opts); % 取每日新增I_new sigma * E采样到与data等长的序列 pred_new diff(Y(:,2)); % 不断累计E的变化和各天新增的对应关系需按实际模型推导 % 简化示意实际需要按每日间隔提取或对E做差分 err sqrt(mean((pred_new - data).^2)) / (max(data) - min(data) eps); end这里有个细节值得展开SEIR的观测量不直接等于I现存感染而往往是每日新增确诊也就是每天从E转到I的人数。所以目标函数里取的不是I本身而是σ·E在每步的累积变化。这个映射关系如果搞错了拟合结果会完全乱套。3.4 HHO优化器主体实现HHO优化器的主体可以按位置更新公式来实现。核心逻辑包括种群初始化、逃逸能量计算、位置更新策略选择以及适应度排序。方法步骤可以这样组织在参数边界内随机初始化种群每个个体的位置就是一组[β, σ, γ]。计算每个个体的适应度记录当前最优位置兔子位置。对每次迭代先更新逃逸能量E再逐个更新每个鹰的位置。对更新后的位置做边界吸收处理避免参数跑到无意义区域。重新计算适应度、更新最优解。记录每次迭代的最优适应度用于后续绘制收敛曲线。HHO位置更新的Matlab示意为function [new_pos] HHO_update(pos, rabbit, E, r, lb, ub) if abs(E) 1 % 探索策略 if r 0.5 new_pos pos - rand() * abs(pos - 2 * r * rabbit); else idx randi(size(pos,1)); new_pos (pos(idx) - rand() * abs(pos(idx) - 2*r*rabbit)); end else % 开发策略软围攻 / 硬围攻 / 渐进式俯冲 % ... 省略若干公式分支 ... end % 边界修正 new_pos max(new_pos, lb); new_pos min(new_pos, ub); end完整的Matlab代码量不算大核心逻辑加注释大概三四百行就能跑通。我在后面的章节会专门讲怎么一步步把整个流程跑起来以及实际运行中会遇到的报错和坑。4. 实操过程从数据准备到结果评估4.1 准备真实数据这里用公开的每日新增确诊数据作为示例。数据格式建议为一行一列或一行一日具体无所谓关键是日期序列要与模型采样的时间点对齐。在导入Matlab时把每日新增序列存成列向量data_new同时生成对应的时间向量tspan 1:length(data_new)时间单位通常为天。这里有个常见困惑SEIR的时间单位到底怎么取。如果传染病数据是逐日统计的时间单位就是天σ的单位就是1/天。对应到潜伏期的含义如果平均潜伏期5天那σ约等于1/50.2。用这些生理学上的先验信息来设置参数范围比盲目给一个很大的边界更高效。同理γ约等于1/感染周期如果平均感染周期为7天γ就在0.14附近。4.2 分步运行的步骤第一步先做数据初始化与参数边界设置边界区间可以是β∈[0.1, 1.5]、σ∈[0.05, 0.6]、γ∈[0.01, 0.5]具体按疫情特征和已有研究调整。第二步定义SEIR微分方程组并封装成函数。第三步写目标函数将模拟输出与真实数据对齐。第四步运行HHO主循环种群规模20、最大迭代次数200通常就有不错效果。第五步输出最优参数重新调用SEIR模型生成拟合曲线与真实数据叠加绘图对比。% 运行示例 data load(covid_daily_new.mat); % 真实每日新增数据 N 10000000; % 区域总人口 param_bounds [0.01 2; 0.05 0.6; 0.01 0.5]; [best_param, best_fit, curve] run_HHO_SEIR(data, N, param_bounds, 20, 200); fprintf(beta%.4f, sigma%.4f, gamma%.4f\n, best_param);如果只用累计感染数据那数据格式上要先做一次累加转换data_cum cumsum(data_new);然后把目标函数里的拟合对象从差分序列改成累计序列。两种方式我都实测过用累计序列拟合时HHO收敛得更快因为目标函数曲面更平滑用每日新增序列拟合时结果更能反映流行过程的峰值特征但需要更多的迭代次数才能稳定。4.3 如何评价拟合效果拟合效果不能只盯着一张图看。我习惯同时看三个指标RMSE误差值、决定系数R²、以及基础再生数R0。RMSE看绝对误差R²看曲线形状的整体贴合度R0 β/γ是传染病模型里最关心的中间量。R0估出来如果明显不合理即使拟合曲线很贴合参数也大概率不可靠——这可能意味着模型结构有误或者数据本身存在偏倚。我在这几个指标上遇到过一个有意思的情况某组参数组合的RMSE非常小但R0高达6.8明显超出了该传染病的现实范围。检查后发现是优化器跑到边界区域、用极端参数硬凑曲线导致的。这种情况需要回到参数边界设置把R0的合理范围折算进去压缩搜索空间。4.4 收敛性与稳定性验证HHO毕竟是随机优化算法每次运行时结果会有小幅波动。为了确保结果可复现我在主循环里加了随机数种子设置并建议多做几次独立运行检查最优参数的标准差。如果标准差偏大说明算法没有稳定收敛到同一个区域这时应该增大最大迭代次数或者调整逃逸能量公式中的随机分布方式。另外一个值得记录的经验是HHO收敛曲线下降速度很快前50代基本完成主要下降后150代只是微调。所以如果你的最大迭代次数只有几十次很可能还没进入充分的开发阶段就停下来了结果会出现明显偏差。我测试过50代和200代的对比后者的RMSE降低了约18%这个提升主要来自开发阶段的精细搜索。5. 常见问题与排查技巧实录5.1 ode45求解失败或结果异常SEIR模型求解整体不复杂但有个容易出问题的地方当感染人数很小比如个位数而总人口有百万甚至千万时S·I/N的量级很小微分方程数值解的绝对容差设置不当会导致求解器认为已经收敛实际却产生了数值偏差。解决办法是用odeset把相对容差RelTol设到1e-6同时启用非负约束NonNegative。如果运行时间过长再把MaxStep适当调大可以显著提速。5.2 优化过程陷入局部最优HHO有较强的全局搜索能力但不能保证100%收敛到全局最优。如果发现每次运行得到的结果差异大、或者拟合曲线在真实数据峰值附近偏移明显优先检查三件事种群规模是否太小建议20到30、最大迭代次数是否不足建议200到500、参数边界是否过宽导致搜索空间过于稀疏。边界过宽的问题经常被忽略比如把β的范围设成0到10可行域里大量区域对应的是完全没有流行病学意义的参数组合算法要浪费大量迭代在这些无用区域上。5.3 参数范围设置的经验法则参数范围设置的依据应该是传染病本身的生物学特征。以新冠为例平均潜伏期约5天左右σ的典型值≈0.2感染周期约7到10天γ的典型值≈0.1到0.15基本再生数R0在2到3之间β的典型值就可以从R0×γ来推算大约在0.2到0.45。这样设置边界后HHO搜索到的参数天然就在一个合理区间内后期人工验证也容易通过。5.4 Matlab运行环境的几个注意点代码在不同Matlab版本下的兼容性基本没问题ode45和基本的矩阵操作几十年没变过。真正容易踩坑的是路径和脚本命名不要把脚本命名为seir.m又同时定义一个seir_ode.m函数容易出现调用混淆也不要用中文路径。另外如果数据量大且迭代次数多Matlab内存会持续增长建议每50次迭代就用clear清理一次临时变量或者把收敛曲线数据预先分配为定长数组而不是动态增长。5.5 输出和绘图技巧最终绘图建议把三张图放到一起真实数据与拟合曲线的对比图、HHO收敛曲线、三个参数的迭代变化轨迹。收敛曲线能直观显示优化过程的健康程度参数轨迹能看出算法是否在早期就跳出了局部区域。绘图时注意线条粗细和图例清晰度set(gca,FontSize,12)这类细节虽然小但能让结果图直接能用在论文汇报里。6. 这套方案的适用边界与扩展方向HHO-SEIR这套组合的适用场景我认为核心是“参数不可直接观测但模型结构相对明确”的问题。除了传染病参数拟合它还可以直接迁移到其他仓室模型SIR、SEIRD、SEIRV只需要改微分方程和目标函数。甚至SIR模型的参数估计也是同一套流程只是少了一个σ。我个人的体会是优化算法本身不是万能的模型选错了参数拟合得再完美也没有意义。SEIR模型有一个隐含假设人群均匀混合、参数在整个流行过程中不变。现实中隔离措施、口罩、疫苗都会让β随时间变化。如果想提高拟合效果可以考虑把β改成随时间分段变化的函数再把分段点也作为优化变量一并搜索。这个思路我试过参数数量从3个涨到5个HHO依然能处理拟合效果明显提升。最后再分享一个小技巧跑HHO的时候把每次迭代的最优参数都保存下来不要只存最后一个值。优化后期如果有某次迭代突然跳出一个拟合误差更小的参数往往不是因为算法更聪明而是因为搜索过程中探索到一个之前没经过的优良区域。把整个参数轨迹留下来你就能在结果分析时看出哪些区域更有可能藏有真实参数这比只看最后输出的一组参数要有价值得多。