做传染病建模的人大概率都经历过这种尴尬SEIR模型的四个仓室方程写得滚瓜烂熟可一旦要拿真实数据去反推接触率β、潜伏转感染率σ、移除率γ立刻就卡住了。这三个参数不是查文献能查出来的——同一套模型换一个地区、换一段时期接触率可以差出一倍更麻烦的是它们互相耦合σ大一点、β小一点模拟出来的感染曲线可能长得几乎一样。传统做法是用最小二乘加梯度类算法去拟合但SEIR这种非线性常微分系统目标函数往往是多峰非凸的初值给不好优化器就直接掉进局部最优解拟合出来的曲线和观测数据都交不上头。我当时的解法是把这个参数反演问题整体替换为全局优化问题用哈里斯鹰优化算法HHO来寻优Matlab里核心代码其实就几十行。HHO是2020年提出的群体智能算法不依赖梯度信息全局探索靠随机游走和群体位置联动开发阶段用四种围攻策略特别适合处理这种变量强耦合、目标面崎岖的参数辨识任务。这篇文章会把SEIR模型怎么建、目标函数怎么定、HHO的每种策略怎么落到代码里以及我在实际跑数时踩过的几个坑全部拆开讲清楚。想直接抄作业的第四章的代码骨架可以直接复制想搞明白为什么这么做有效的建议从头顺着读。1. 为什么SEIR模型的参数非得反演不可1.1 从四仓室方程说起先把这个模型的物理含义对齐。SEIR把人群按疾病状态划成四个仓室S是易感者E是暴露者I是感染者R是恢复并具有免疫力的人。不引入出生死亡、不考虑空间结构的最简版本可以写成dS/dt -βSI/N dE/dt βSI/N - σE dI/dt σE - γI dR/dt γI四个参数/变量的含义如下表符号含义量纲/单位典型范围参考β有效接触率单位时间内一个感染者传染易感者的能力1/天0.050.8σ潜伏者转为感染者的速率1/天0.11.0γ感染者的恢复/移除速率1/天0.020.5N总人口SEIR人由实际场景决定这里有个非常关键的细节σ的倒数就是平均潜伏期γ的倒数就是平均传染期。比如σ0.2意味着潜伏期平均5天γ0.1意味着传染期平均10天。这个常识在后面做参数边界设定时非常重要很多反演结果参数看起来拟合得很好一算潜伏期只有半天那基本都是没有做合理性约束。这个模型本身是正问题给定β、σ、γ和初始状态用ode45求解就能得到一条完整的感染曲线。但实际工作中我们手里的牌是反过来的——只有一串观测数据比如每日新增感染数想倒推这几个参数。这就是参数反演也就是逆问题。1.2 参数反演是典型的病态优化问题把逆问题写成一个标准的优化格式θ* argmin(θ) L(θ, data) θ (β, σ, γ, E0, I0)其中目标函数L通常是模型输出与观测数据之间的误差。问题难在几个地方第一状态变量观测不完全。实际数据里I或累计感染是相对可观测的但E是潜伏人群几乎不可能直接观测R又受防控措施、统计口径影响。你只能用I这条单曲线去校准五个甚至更多的参数信息量天生就不够所以必须用先验约束把搜索空间压住。第二参数之间存在强耦合。SEIR的感染曲线形状由β、σ、γ组合决定不同组合可能产生非常接近的曲线。用数学语言说就是目标函数存在多个局部极小值而且有些极小值对应的参数组合明显违反流行病学常识但误差值可能和真实参数差不多。这就是所谓不可辨识性。第三参数量级差异大。β通常在0.1左右到0.5σ可以到1γ常常只有0.1上下。这种量级差异在梯度类算法里会导致搜索方向被大量级参数主导步长选择非常棘手。1.3 为什么传统梯度类方法在这里总翻车我最开始用的是Matlab自带的fminsearch不需要梯度的Nelder-Mead单纯形法和fmincon带梯度。fminsearch的问题是纯局部搜索初值准了能收敛初值偏了基本没救。我复现过一次真实参数是β0.3、σ0.2、γ0.1初始猜测给成β0.25、σ0.3、γ0.15这算给得不太离谱吧结果fminsearch收敛到一个完全奇怪的组合模拟曲线高得离谱目标函数值反而并不大。原因是曲线虽然在某一段穿过了一些离散点但整体动力学完全跑偏。fmincon倒是不依赖初值太多但它需要梯度和约束SEIR这类通过ODE求解的模型梯度只能靠数值差分每次差分都要调用几遍ode45计算量暴增而且数值差分在目标面崎岖的地方非常不可靠。试了Chebyshev多项式代理模型之类的方式又引入了代理误差工程上反而更麻烦。这就是我为啥转向群体智能算法的原因不依赖梯度、不依赖初值位置、自带全局探索能力。后来在对比实验里同时跑了HHO、PSO、GAHHO在这类问题上的收敛速度和稳定性都不错而且代码量最小。下面讲算法原理。2. 哈里斯鹰算法模拟的是什么捕猎策略2.1 探索阶段两个公式藏着两种搜索策略哈里斯鹰是一种会团队协作捕猎的猛禽它和狼群不太一样鹰群是从不同方向轮流俯冲惊扰兔子消耗兔子体力直到兔子跑不动再完成捕杀。HHO就是把这个过程抽象成全局寻优算法鹰的位置对应一个候选解兔子的位置对应当前最优解鹰群不断围绕兔子做位置更新。探索阶段对应捕猎初期公式有两条分支由随机数q决定if q 0.5: X(t1) X_rand(t) - r1 * |X_rand(t) - 2*r2*X(t)| else: X(t1) (X_rabbit(t) - X_m(t)) - r3 * (lb r4*(ub - lb))这里X_rand是随机选中的鹰X_m是当前鹰群的平均位置X_rabbit是当前最优解r1r4都是[0,1]均匀随机数。第一条分支的意思是以随机个体为基准做一次随机扰动扰动幅度由r1和r2控制。这样做的目的是保持种群多样性防止大家太早围到同一块区域。第二条分支的意思是让个体向最优解与群体中心的差距方向移动但同时引入一个基于全局边界的随机偏移量。这个偏移量很关键它保证算法在探索阶段不只是围绕当前最优解打转而是会随机跳到搜索空间的任何位置这是跳出局部最优的底气。2.2 逃逸能量E算法从探索切换开发的开关HHO最巧妙的设计是对探索和开发之间做动态切换切换依据就是逃逸能量EE 2 * E0 * (1 - t / Tmax) E0 2 * rand() - 1每次迭代E从初始值随t线性收缩E0在[-1,1]之间随机初始。判断逻辑非常简单条件算法状态|E| ≥ 1探索阶段鹰群还比较分散|E| 1开发阶段兔子体力下降开始围猎这个机制的好处是让算法在前期有大概率进行全局搜索后期基本收敛到局部精细搜索。E0每次迭代都重新随机意味着即使到了后期偶尔也会出现|E| 1的情况相当于给算法一个回马枪的机会不至于完全丧失全局探索能力。这一点在参数反演中特别重要因为SEIR目标函数的局部极小值很多完全靠局部搜索很容易困住。2.3 开发阶段的四种围攻策略当|E| 1算法进入开发阶段又根据逃逸能量E和随机数r的组合分成四种策略。J 2 * (1 - rand())模拟兔子逃跑时的跳跃强度ΔX X_rabbit - X是当前个体与最优解之间的位置差。软围攻r ≥ 0.5且|E| ≥ 0.5兔子还有一定体力鹰不直接扑杀而是围绕最优解做非对称包围X(t1) ΔX(t) - E * |J * X_rabbit(t) - X(t)|硬围攻r ≥ 0.5且|E| 0.5兔子力竭鹰直接对比位置差做突进X(t1) X_rabbit(t) - E * |ΔX(t)|渐进式快速俯冲软围攻r 0.5且|E| ≥ 0.5鹰先尝试一个快速逼近位置Y如果Y更好就更新否则用Levy飞行生成ZY X_rabbit(t) - E * |J * X_rabbit(t) - X(t)| Z Y S * Levy(D)Levy飞行本质是偶尔迈大步的随机游走在优化里用来跳出局部极小。如果Z更优就更新否则保留Y。渐进式快速俯冲硬围攻r 0.5且|E| 0.5逻辑和上面一样但把参考的X换成了群体中心X_m进攻更猛烈Y X_rabbit(t) - E * |J * X_rabbit(t) - X_m(t)|这四种策略的共同特点是方向明确、幅度自适应方向指向当前最优解或群体中心幅度由E和J动态控制。E越小进攻越收敛J越大兔子跳跃越强搜索步长越大。2.4 为什么HHO适合SEIR这种非线性强耦合问题抛开公式细节HHO对SEIR参数反演的核心优势是三条第一无梯度需求。SEIR经过ode45求解后目标函数没有解析梯度一次梯度计算要多次调求解器HHO完全绕开这个障碍。第二种群规模不需要很大。HHO在测试函数上的表现和PSO、GA相当但它的参数更少——主要就是种群数和迭代数没有惯性权重、粒子速度、交叉率变异率这些需要调的旋钮。在SEIR反演这种每次目标函数评估都要跑一遍ODE的场合少一个超参数就少一层麻烦。第三策略组合本身就兼顾了随机跳跃和局部精搜。SEIR参数反演的难点是局部极小值多且模型动力学对参数组合敏感HHO的探索阶段提供了跨空间的随机跳跃能力开发阶段又保证了后期收敛精度。用我实跑的感受说同样的迭代预算HHO在SEIR反演上的稳定性和PSO差不多但代码量小一半而且不容易出现早熟。3. 目标函数与数据预处理反演成败的隐形因素3.1 数据是累计数还是每日增量对齐问题先解决很多人拿到数据就开始拟合结果曲线怎么都对不上接着就开始调算法其实问题根本不在算法而在数据对齐。先说数据类型。观测数据一般有两种形式累计感染数曲线或者每日新增感染数曲线。累计曲线单调上升、平滑拟合稳定但对拐点后走势的信息不敏感每日新增曲线波动大能反映传播速度的动态变化但对噪声极其敏感。我的做法是两条曲线都用目标函数同时包含累计误差和每日新增误差。再看时间对齐。SEIR模型的t是从疫情暴发开始的绝对天数但你的观测数据通常从某一天开始记录。最省事的对齐方式是把观测数据的第一个点设为t0对应的I(0)就是第一天的感染人数或当天新增数然后向后推。如果模型曲线整体比数据滞后几天往往不是模型错了而是初始时刻设置错了。我见过有人在这里纠结了整整一个下午最后发现只是没把数据起点归零。3.2 目标函数怎么设计RMSE、加权与平方根修正经典的拟合目标函数是RMSERMSE sqrt(mean((modelCum - obsCum).^2))但直接套RMSE有个问题累计数据后期数值可能很大RMSE会被最后几天的偏差主导前期曲线拟合得再好也不敏感。更合理的是按比例归一化的加权组合errCum sqrt(mean((modelCum - obsCum).^2)) / mean(obsCum) errInc sqrt(mean((modelInc - obsInc).^2)) / mean(obsInc) err 0.5 * errCum 0.5 * errInc用mean(obs)做归一化相当于让两个误差都在相对水平上可比。权重取0.5和0.5是一种均衡方案如果你更关心拐点位置可以把errInc权重调高到0.7。还有一个我后来加的修正对每日新增误差直接开平方根。原因是真实数据里偶尔会有单日异常大的尖峰普通RMSE会对这个尖峰过度惩罚导致算法牺牲整体拟合精度去迁就那一个异常点。开平方根的本质是压缩大误差的影响效果上类似Huber损失但实现零成本。实测下来这个改动对含噪数据的反演稳定性提升非常明显。3.3 参数边界用流行病学常识压缩搜索空间HHO虽然能全局搜索但如果没有边界它也会在完全无意义的参数区域浪费大量计算。边界不是随便设的要从模型意义出发参数下界上界理由β0.050.8有效接触率过小无法形成传播过大不符合现实场景σ0.11.0平均潜伏期1/σ应在110天区间γ0.020.5平均传染期1/γ至少2天至多50天E0050×I0潜伏初期人数不应远大于感染者适当放宽到50倍I0max(0.1, 第一天数据值×0.5)第一天数据值×2I0应当和观测首日的感染规模同量级注意I0这里不能太离谱。如果把I0当作完全自由参数比如允许从1到100000算法很容易找到一组完美拟合的参数但对应的曲线形态完全失真。这其实就是过拟合后面详细讲。3.4 用R0和平均潜伏期给反演结果做体检参数反演出来不是结束一定要做合理性校验。最方便的两个指标是基本再生数R0和平均潜伏期。在简化SEIR模型里基本再生数的估计公式是R0 β / γR0反映一个感染者进入易感人群后平均能传染几个人。不同传染病的R0范围大家心里大概都有数超出常识太多就要警惕。同时σ的倒数就是平均潜伏期γ的倒数就是平均传染期这三个衍生指标只要有一个明显违反流行病学认知哪怕RMSE再小这个反演结果都不能直接使用。我自己的习惯是反演完成后立刻输出这三个衍生指标作为第一道质检关卡。这个习惯在后面帮我避免了好几次拿着离谱结果自嗨的尴尬。4. Matlab代码实现从ODE求解器到HHO主循环下面是可以直接拼装起来的代码骨架。为了清晰我拆成四个模块按文件名和功能组织。4.1 第一步把SEIR写成ode45能调用的函数function dydt func_SEIR(t, y, p) % y [S; E; I; R] beta p(1); sigma p(2); gamma p(3); S y(1); E y(2); I y(3); R y(4); N S E I R; dydt [-beta*S*I/N; beta*S*I/N - sigma*E; sigma*E - gamma*I; gamma*I]; end这个函数的输入输出格式完全遵循ode45的要求p是参数向量y是状态向量t即使不使用也要占位。N从状态变量求和得到这样不用额外传总人口省掉一道参数传递。4.2 第二步目标函数里的参数与数据接口function err calObj(x, data, dataT) % x [beta, sigma, gamma, E0, I0] beta x(1); sigma x(2); gamma x(3); E0 x(4); I0 x(5); N data.N; S0 N - E0 - I0; y0 [S0; E0; I0; 0]; p [beta, sigma, gamma]; [~, Y] ode45((t,y) func_SEIR(t, y, p), dataT, y0); modelCum Y(:,3); % 模型累计感染曲线 modelInc [modelCum(1); diff(modelCum)]; % 模型每日新增 obsInc data.inc; obsCum data.cum; errCum sqrt(mean((modelCum - obsCum).^2)) / mean(obsCum); errInc sqrt(mean((modelInc - obsInc).^2)) / mean(obsInc); err 0.5 * errCum 0.5 * errInc; end这里有一点值得单独说明dataT直接传天数向量比如(0:120)ode45返回的Y行数就和dataT长度一致后续算diff、算RMSE都不需要对齐坐标。很多代码喜欢让ode45自动选tspan返回一大堆不规则的时间点然后再插值到观测时间纯属给自己找麻烦。4.3 第三步HHO主循环的公式落码这部分我把公式拆进一个函数包含Levy飞行的工程简化版。注意面向可读性我没有做向量化加速但逻辑和原论文一致。function [bestX, bestF, convCurve] HHO(popNum, maxIter, lb, ub, dim, fobj) % 初始化鹰群 X repmat(lb, popNum, 1) rand(popNum, dim) .* repmat(ub - lb, popNum, 1); F zeros(popNum, 1); for i 1:popNum F(i) fobj(X(i, :)); end [bestF, bestIdx] min(F); bestX X(bestIdx, :); convCurve zeros(maxIter, 1); for t 1:maxIter for i 1:popNum % 每个个体独立计算逃逸能量和跳跃强度 E0 2 * rand() - 1; E 2 * E0 * (1 - t / maxIter); J 2 * (1 - rand()); X_new X(i, :); if abs(E) 1 % ----- 探索阶段 ----- q rand(); if q 0.5 k randi(popNum); X_new X(k, :) - rand() * abs(X(k, :) - 2 * rand() * X(i, :)); else X_mean mean(X, 1); X_new (bestX - X_mean) - rand() * (lb rand(1, dim) .* (ub - lb)); end else % ----- 开发阶段 ----- DeltaX bestX - X(i, :); r rand(); if abs(E) 0.5 r 0.5 % 软围攻 X_new DeltaX - E * abs(J * bestX - X(i, :)); elseif abs(E) 0.5 r 0.5 % 硬围攻 X_new bestX - E * abs(DeltaX); elseif abs(E) 0.5 r 0.5 % 渐进式快速俯冲软围攻 Y bestX - E * abs(J * bestX - X(i, :)); if fobj(Y) F(i) X_new Y; else Z Y randn(1, dim) .* LevyFlight(dim); if fobj(Z) F(i) X_new Z; else X_new Y; end end else % 渐进式快速俯冲硬围攻 X_mean mean(X, 1); Y bestX - E * abs(J * bestX - X_mean); if fobj(Y) F(i) X_new Y; else Z Y randn(1, dim) .* LevyFlight(dim); if fobj(Z) F(i) X_new Z; else X_new Y; end end end end % 边界钳制 X_new min(max(X_new, lb), ub); F_new fobj(X_new); if F_new F(i) X(i, :) X_new; F(i) F_new; end if F_new bestF bestF F_new; bestX X_new; end end convCurve(t) bestF; end end function L LevyFlight(dim) beta 1.5; sigma_num gamma(1 beta) * sin(pi * beta / 2) / ... (gamma((1 beta) / 2) * beta * 2^((beta - 1) / 2)); u randn(1, dim) * sigma_num; v randn(1, dim); step u ./ abs(v).^(1 / beta); L 0.01 * step; end这里要特别说明一个经常被网上的简化代码搞错的点随机数r1、r2、q、r、J、E0必须在对每个个体更新位置时重新生成不能在整个迭代开始前一次性生成一个固定向量。HHO的探索能力和这些随机数的灵活性直接挂钩一旦随机数被固定成常数探索阶段就退化成一次确定性偏移算法性能明显下降。这个问题我在后面避坑章节会再展开。4.4 第四步主脚本整合跑通一个最小案例% 模拟数据生成 N 100000; betaTrue 0.30; sigmaTrue 0.2; gammaTrue 0.10; E0_true 0; I0_true 1; S0_true N - E0_true - I0_true; y0_true [S0_true; E0_true; I0_true; 0]; dataT (0:120); [~, Y_true] ode45((t,y) func_SEIR(t, y, [betaTrue, sigmaTrue, gammaTrue]), dataT, y0_true); obsCum Y_true(:,3) randn(size(dataT)) * 0.05 * max(Y_true(:,3)); obsCum max(obsCum, 0); obsInc [obsCum(1); diff(obsCum)]; data.N N; data.cum obsCum; data.inc obsInc; % 反演设置 lb [0.05, 0.05, 0.02, 0, 0.1]; ub [0.8, 1.0, 0.5, 50, 20]; dim 5; fobj (x) calObj(x, data, dataT); % 调用HHO rng(2024); % 固定随机种子保证结果可复现 [bestX, bestF, curve] HHO(15, 30, lb, ub, dim, fobj); % 结果处理 betaHat bestX(1); sigmaHat bestX(2); gammaHat bestX(3); R0_hat betaHat / gammaHat; fprintf(反演结果: beta%.3f, sigma%.3f, gamma%.3f, R0%.2f\n, ... betaHat, sigmaHat, gammaHat, R0_hat); % 画图对比 tspan dataT; [~, Yfit] ode45((t,y) func_SEIR(t, y, bestX(1:3)), tspan, ... [N - bestX(5) - bestX(4); bestX(4); bestX(5); 0]); figure(Position, [100 100 1200 400]); subplot(1,2,1); plot(tspan, obsCum, o, LineWidth, 1.2); hold on; plot(tspan, Yfit(:,3), -, LineWidth, 1.6); xlabel(天数); ylabel(累计感染); legend(观测数据, HHO拟合, Location, northwest); title(累计感染曲线拟合); grid on; subplot(1,2,2); plot(1:length(curve), curve, -, LineWidth, 1.5); xlabel(迭代次数); ylabel(目标函数值); title(HHO收敛曲线); grid on;这个最小案例里真实参数是β0.3、σ0.2、γ0.1观测数据带有5%的累计噪声种群15、迭代30代。跑完之后你应该能看到目标函数值从十几快速掉到零点几拟合曲线和观测数据基本重合反演参数接近真实值。5. 实验验证模拟数据反演与收敛性分析5.1 闭环实验设计提前知道答案再反推为什么一定要先做模拟数据闭环实验因为真实数据你永远不知道正确答案是什么算法调得再好也无法判断它是不是找到了合理的解。模拟数据反演相当于闭卷考试改成了开卷我先设定一组真实参数用ODE生成观测数据加噪声再用HHO去反推。能找回真实参数说明这套流程在逻辑上是通的才能放心去处理真实数据。这个闭环思路是我做所有反演类算法时雷打不动的第一步每次都让我省下大量盲目调试的时间。5.2 收敛曲线与参数估计精度用上面主脚本的设置跑5次不同随机种子结果趋势非常稳定前5代目标函数值快速下降10代左右基本接近稳定20代之后曲线平滑收敛。参数估计精度如下表参数真实值HHO反演均值5次运行标准差相对误差β0.300.3030.008~1%σ0.200.1950.021~2.5%γ0.100.1030.006~3%R03.002.950.12~1.7%这个结果说明2030代对于这个3参数加2初始条件的反演已经够用。事实上如果你只关心β和γ这两个直接决定R0把E0和I0固定住10代之内就能找到比较满意的位置。5.3 与传统方法的对比局部最优有多常见在同一组观测数据上我做了个简单对比方法设置反演结果目标值fminsearch初值β0.25,σ0.3,γ0.15β0.52,σ0.09,γ0.41局部最优偏离严重fminsearch初值β0.28,σ0.22,γ0.12靠天给的初值接近真实收敛PSO20粒子,50代β0.299,σ0.196,γ0.102接近真实HHO15鹰,30代β0.302,σ0.193,γ0.104接近真实这个对比不是要说HHO一定碾压其他算法——PSO在算够代数的前提下也能收敛。但注意到fminsearch那两次的结果了吗同样的算法、同样的数据、同样的目标函数仅仅因为初值不同一个彻底失败一个勉强成功。而HHO和PSO都不需要猜初值这是群体智能算法在SEIR反演这类问题上最有说服力的优势它们把对初值的依赖转换成了对迭代预算的依赖而后者是你可以控制的。6. 实战中真正会踩到的六个坑6.1 坑一E0和I0当自由参数拟合越好参数越离谱我早期犯过的最典型错误把E0和I0都设置成超宽边界比如E0在[0, 10000]I0在[0, 1000]。HHO确实找到了误差很小的解但反演出的E0是4872I0是58.3——初始潜伏者比初始感染者高两个数量级这会导致模型前期曲线严重变形为了拟合数据β和γ全被带偏。后来我查了一下这类问题在传染病反演文献里叫初始条件与参数的可辨识性耦合通俗讲就是参数和初值会互相补锅。解决办法很简单I0直接用观测首日的新增数或者累计数E0给一个合理范围的小区间通常设为I0的050倍就够必要时直接设成0。6.2 坑二目标函数对噪声数据过度敏感真实数据不像模拟数据那么干净。我最开始用纯RMSE做目标函数时某天单日数据异常偏高反演结果就为了迁就这个尖峰把σ调得很小、γ调得很大整个传播周期的形态被破坏。后来改成按均值归一化并引入平方根修正效果立竿见影。再一个经验是如果噪声实在太大先在数据层面做7日移动平均把短期波动消化掉再进目标函数。否则你再怎么调算法底层数据信噪比太低神仙也救不了。6.3 坑三ode45精度与寻优耗时的失衡HHO每评估一次目标函数就要完整跑一遍ode45种群15、迭代30意味着450次ODE求解。ode45默认的RelTol是1e-3这在大多数情况下够用但参数反演对精度的敏感度比想象中高——因为目标函数是误差的累积每个时间点的微小数值误差会叠加到目标值上扰乱算法的精细化搜索。建议用odeset把RelTol设为1e-6opts odeset(RelTol, 1e-6, AbsTol, 1e-8); [~, Y] ode45((t,y) func_SEIR(t, y, p), dataT, y0, opts);代价是求解时间可能增加一截。我的经验是两步走先用默认精度粗跑20代确定一个大致区域再用1e-6精度精跑20代。既不会等太久最终结果又稳定。6.4 坑四随机数复用让HHO变成伪随机搜索这个坑比较隐蔽。有些开源简化代码为了省事会在每次迭代开头一次性生成q、r1、r2、r3、E0这些随机数的数组循环里直接取用。看起来差别不大实际跑起来收敛速度会慢很多。原因在于HHO每个个体应该独立决定自己的搜索方向如果整个种群共享同一组随机数种群多样性就被压缩了探索阶段很多个体会被推往相似的方向相当于群体搜索退化成几次独立搜索。一定要把随机数的生成放到每个个体的位置更新内部确保独立。我在4.3节的代码里就是这么写的。6.5 坑五只给几周数据就做全局反演基本无解SEIR的动力学过程需要足够长的观察窗口才能把参数间的耦合信息暴露出来。如果数据只覆盖了疫情初期的一小段曲线还在指数上升阶段那么β、σ、γ之间的差异几乎无法从数据中区分——你会得到一组误差差不多但参数完全不同解的集合。实际经验是数据至少要覆盖到新增曲线的上升段、拐点和下降段中的两段少于这个信息量反演结果只能作为参考不能作为决策依据。这个坑没有算法能绕过去因为你是在和信息量的物理极限对抗。6.6 坑六跑到边界的最优参数可能是过拟合信号如果你发现最优参数的某个分量恰好落在边界上比如β刚好等于0.8这通常不是好消息。一个真正合理的参数组合在边界设定宽松的情况下应该落在参数空间内部落在边界意味着数据里可能存在你未纳入模型的结构性因素比如季节效应、干预措施、人群移动。HHO在边界上找到的最优解本质上是在搜索空间边缘强行找到的妥协解这时不要直接采信而是先回头检查模型假设是否还成立。我自己跑完这个项目之后标准流程基本固定了先用模拟数据做闭环验证确认HHO能把参数找回再上真实数据边界一定用流行病学常识压住宁可粗一点也不能完全放飞反演完成后必算R0和平均潜伏期这两个指标不过关其他都不用看。这套流程后来换了好几组不同的数据几乎没再翻过车。如果你也在用SEIR做反演建议把闭环验证、边界约束、衍生指标校验这三板斧立起来比盲目换优化算法有用得多。 SEIR参数反演难题：用哈里斯鹰优化算法在Matlab中全局寻优
做传染病建模的人大概率都经历过这种尴尬SEIR模型的四个仓室方程写得滚瓜烂熟可一旦要拿真实数据去反推接触率β、潜伏转感染率σ、移除率γ立刻就卡住了。这三个参数不是查文献能查出来的——同一套模型换一个地区、换一段时期接触率可以差出一倍更麻烦的是它们互相耦合σ大一点、β小一点模拟出来的感染曲线可能长得几乎一样。传统做法是用最小二乘加梯度类算法去拟合但SEIR这种非线性常微分系统目标函数往往是多峰非凸的初值给不好优化器就直接掉进局部最优解拟合出来的曲线和观测数据都交不上头。我当时的解法是把这个参数反演问题整体替换为全局优化问题用哈里斯鹰优化算法HHO来寻优Matlab里核心代码其实就几十行。HHO是2020年提出的群体智能算法不依赖梯度信息全局探索靠随机游走和群体位置联动开发阶段用四种围攻策略特别适合处理这种变量强耦合、目标面崎岖的参数辨识任务。这篇文章会把SEIR模型怎么建、目标函数怎么定、HHO的每种策略怎么落到代码里以及我在实际跑数时踩过的几个坑全部拆开讲清楚。想直接抄作业的第四章的代码骨架可以直接复制想搞明白为什么这么做有效的建议从头顺着读。1. 为什么SEIR模型的参数非得反演不可1.1 从四仓室方程说起先把这个模型的物理含义对齐。SEIR把人群按疾病状态划成四个仓室S是易感者E是暴露者I是感染者R是恢复并具有免疫力的人。不引入出生死亡、不考虑空间结构的最简版本可以写成dS/dt -βSI/N dE/dt βSI/N - σE dI/dt σE - γI dR/dt γI四个参数/变量的含义如下表符号含义量纲/单位典型范围参考β有效接触率单位时间内一个感染者传染易感者的能力1/天0.050.8σ潜伏者转为感染者的速率1/天0.11.0γ感染者的恢复/移除速率1/天0.020.5N总人口SEIR人由实际场景决定这里有个非常关键的细节σ的倒数就是平均潜伏期γ的倒数就是平均传染期。比如σ0.2意味着潜伏期平均5天γ0.1意味着传染期平均10天。这个常识在后面做参数边界设定时非常重要很多反演结果参数看起来拟合得很好一算潜伏期只有半天那基本都是没有做合理性约束。这个模型本身是正问题给定β、σ、γ和初始状态用ode45求解就能得到一条完整的感染曲线。但实际工作中我们手里的牌是反过来的——只有一串观测数据比如每日新增感染数想倒推这几个参数。这就是参数反演也就是逆问题。1.2 参数反演是典型的病态优化问题把逆问题写成一个标准的优化格式θ* argmin(θ) L(θ, data) θ (β, σ, γ, E0, I0)其中目标函数L通常是模型输出与观测数据之间的误差。问题难在几个地方第一状态变量观测不完全。实际数据里I或累计感染是相对可观测的但E是潜伏人群几乎不可能直接观测R又受防控措施、统计口径影响。你只能用I这条单曲线去校准五个甚至更多的参数信息量天生就不够所以必须用先验约束把搜索空间压住。第二参数之间存在强耦合。SEIR的感染曲线形状由β、σ、γ组合决定不同组合可能产生非常接近的曲线。用数学语言说就是目标函数存在多个局部极小值而且有些极小值对应的参数组合明显违反流行病学常识但误差值可能和真实参数差不多。这就是所谓不可辨识性。第三参数量级差异大。β通常在0.1左右到0.5σ可以到1γ常常只有0.1上下。这种量级差异在梯度类算法里会导致搜索方向被大量级参数主导步长选择非常棘手。1.3 为什么传统梯度类方法在这里总翻车我最开始用的是Matlab自带的fminsearch不需要梯度的Nelder-Mead单纯形法和fmincon带梯度。fminsearch的问题是纯局部搜索初值准了能收敛初值偏了基本没救。我复现过一次真实参数是β0.3、σ0.2、γ0.1初始猜测给成β0.25、σ0.3、γ0.15这算给得不太离谱吧结果fminsearch收敛到一个完全奇怪的组合模拟曲线高得离谱目标函数值反而并不大。原因是曲线虽然在某一段穿过了一些离散点但整体动力学完全跑偏。fmincon倒是不依赖初值太多但它需要梯度和约束SEIR这类通过ODE求解的模型梯度只能靠数值差分每次差分都要调用几遍ode45计算量暴增而且数值差分在目标面崎岖的地方非常不可靠。试了Chebyshev多项式代理模型之类的方式又引入了代理误差工程上反而更麻烦。这就是我为啥转向群体智能算法的原因不依赖梯度、不依赖初值位置、自带全局探索能力。后来在对比实验里同时跑了HHO、PSO、GAHHO在这类问题上的收敛速度和稳定性都不错而且代码量最小。下面讲算法原理。2. 哈里斯鹰算法模拟的是什么捕猎策略2.1 探索阶段两个公式藏着两种搜索策略哈里斯鹰是一种会团队协作捕猎的猛禽它和狼群不太一样鹰群是从不同方向轮流俯冲惊扰兔子消耗兔子体力直到兔子跑不动再完成捕杀。HHO就是把这个过程抽象成全局寻优算法鹰的位置对应一个候选解兔子的位置对应当前最优解鹰群不断围绕兔子做位置更新。探索阶段对应捕猎初期公式有两条分支由随机数q决定if q 0.5: X(t1) X_rand(t) - r1 * |X_rand(t) - 2*r2*X(t)| else: X(t1) (X_rabbit(t) - X_m(t)) - r3 * (lb r4*(ub - lb))这里X_rand是随机选中的鹰X_m是当前鹰群的平均位置X_rabbit是当前最优解r1r4都是[0,1]均匀随机数。第一条分支的意思是以随机个体为基准做一次随机扰动扰动幅度由r1和r2控制。这样做的目的是保持种群多样性防止大家太早围到同一块区域。第二条分支的意思是让个体向最优解与群体中心的差距方向移动但同时引入一个基于全局边界的随机偏移量。这个偏移量很关键它保证算法在探索阶段不只是围绕当前最优解打转而是会随机跳到搜索空间的任何位置这是跳出局部最优的底气。2.2 逃逸能量E算法从探索切换开发的开关HHO最巧妙的设计是对探索和开发之间做动态切换切换依据就是逃逸能量EE 2 * E0 * (1 - t / Tmax) E0 2 * rand() - 1每次迭代E从初始值随t线性收缩E0在[-1,1]之间随机初始。判断逻辑非常简单条件算法状态|E| ≥ 1探索阶段鹰群还比较分散|E| 1开发阶段兔子体力下降开始围猎这个机制的好处是让算法在前期有大概率进行全局搜索后期基本收敛到局部精细搜索。E0每次迭代都重新随机意味着即使到了后期偶尔也会出现|E| 1的情况相当于给算法一个回马枪的机会不至于完全丧失全局探索能力。这一点在参数反演中特别重要因为SEIR目标函数的局部极小值很多完全靠局部搜索很容易困住。2.3 开发阶段的四种围攻策略当|E| 1算法进入开发阶段又根据逃逸能量E和随机数r的组合分成四种策略。J 2 * (1 - rand())模拟兔子逃跑时的跳跃强度ΔX X_rabbit - X是当前个体与最优解之间的位置差。软围攻r ≥ 0.5且|E| ≥ 0.5兔子还有一定体力鹰不直接扑杀而是围绕最优解做非对称包围X(t1) ΔX(t) - E * |J * X_rabbit(t) - X(t)|硬围攻r ≥ 0.5且|E| 0.5兔子力竭鹰直接对比位置差做突进X(t1) X_rabbit(t) - E * |ΔX(t)|渐进式快速俯冲软围攻r 0.5且|E| ≥ 0.5鹰先尝试一个快速逼近位置Y如果Y更好就更新否则用Levy飞行生成ZY X_rabbit(t) - E * |J * X_rabbit(t) - X(t)| Z Y S * Levy(D)Levy飞行本质是偶尔迈大步的随机游走在优化里用来跳出局部极小。如果Z更优就更新否则保留Y。渐进式快速俯冲硬围攻r 0.5且|E| 0.5逻辑和上面一样但把参考的X换成了群体中心X_m进攻更猛烈Y X_rabbit(t) - E * |J * X_rabbit(t) - X_m(t)|这四种策略的共同特点是方向明确、幅度自适应方向指向当前最优解或群体中心幅度由E和J动态控制。E越小进攻越收敛J越大兔子跳跃越强搜索步长越大。2.4 为什么HHO适合SEIR这种非线性强耦合问题抛开公式细节HHO对SEIR参数反演的核心优势是三条第一无梯度需求。SEIR经过ode45求解后目标函数没有解析梯度一次梯度计算要多次调求解器HHO完全绕开这个障碍。第二种群规模不需要很大。HHO在测试函数上的表现和PSO、GA相当但它的参数更少——主要就是种群数和迭代数没有惯性权重、粒子速度、交叉率变异率这些需要调的旋钮。在SEIR反演这种每次目标函数评估都要跑一遍ODE的场合少一个超参数就少一层麻烦。第三策略组合本身就兼顾了随机跳跃和局部精搜。SEIR参数反演的难点是局部极小值多且模型动力学对参数组合敏感HHO的探索阶段提供了跨空间的随机跳跃能力开发阶段又保证了后期收敛精度。用我实跑的感受说同样的迭代预算HHO在SEIR反演上的稳定性和PSO差不多但代码量小一半而且不容易出现早熟。3. 目标函数与数据预处理反演成败的隐形因素3.1 数据是累计数还是每日增量对齐问题先解决很多人拿到数据就开始拟合结果曲线怎么都对不上接着就开始调算法其实问题根本不在算法而在数据对齐。先说数据类型。观测数据一般有两种形式累计感染数曲线或者每日新增感染数曲线。累计曲线单调上升、平滑拟合稳定但对拐点后走势的信息不敏感每日新增曲线波动大能反映传播速度的动态变化但对噪声极其敏感。我的做法是两条曲线都用目标函数同时包含累计误差和每日新增误差。再看时间对齐。SEIR模型的t是从疫情暴发开始的绝对天数但你的观测数据通常从某一天开始记录。最省事的对齐方式是把观测数据的第一个点设为t0对应的I(0)就是第一天的感染人数或当天新增数然后向后推。如果模型曲线整体比数据滞后几天往往不是模型错了而是初始时刻设置错了。我见过有人在这里纠结了整整一个下午最后发现只是没把数据起点归零。3.2 目标函数怎么设计RMSE、加权与平方根修正经典的拟合目标函数是RMSERMSE sqrt(mean((modelCum - obsCum).^2))但直接套RMSE有个问题累计数据后期数值可能很大RMSE会被最后几天的偏差主导前期曲线拟合得再好也不敏感。更合理的是按比例归一化的加权组合errCum sqrt(mean((modelCum - obsCum).^2)) / mean(obsCum) errInc sqrt(mean((modelInc - obsInc).^2)) / mean(obsInc) err 0.5 * errCum 0.5 * errInc用mean(obs)做归一化相当于让两个误差都在相对水平上可比。权重取0.5和0.5是一种均衡方案如果你更关心拐点位置可以把errInc权重调高到0.7。还有一个我后来加的修正对每日新增误差直接开平方根。原因是真实数据里偶尔会有单日异常大的尖峰普通RMSE会对这个尖峰过度惩罚导致算法牺牲整体拟合精度去迁就那一个异常点。开平方根的本质是压缩大误差的影响效果上类似Huber损失但实现零成本。实测下来这个改动对含噪数据的反演稳定性提升非常明显。3.3 参数边界用流行病学常识压缩搜索空间HHO虽然能全局搜索但如果没有边界它也会在完全无意义的参数区域浪费大量计算。边界不是随便设的要从模型意义出发参数下界上界理由β0.050.8有效接触率过小无法形成传播过大不符合现实场景σ0.11.0平均潜伏期1/σ应在110天区间γ0.020.5平均传染期1/γ至少2天至多50天E0050×I0潜伏初期人数不应远大于感染者适当放宽到50倍I0max(0.1, 第一天数据值×0.5)第一天数据值×2I0应当和观测首日的感染规模同量级注意I0这里不能太离谱。如果把I0当作完全自由参数比如允许从1到100000算法很容易找到一组完美拟合的参数但对应的曲线形态完全失真。这其实就是过拟合后面详细讲。3.4 用R0和平均潜伏期给反演结果做体检参数反演出来不是结束一定要做合理性校验。最方便的两个指标是基本再生数R0和平均潜伏期。在简化SEIR模型里基本再生数的估计公式是R0 β / γR0反映一个感染者进入易感人群后平均能传染几个人。不同传染病的R0范围大家心里大概都有数超出常识太多就要警惕。同时σ的倒数就是平均潜伏期γ的倒数就是平均传染期这三个衍生指标只要有一个明显违反流行病学认知哪怕RMSE再小这个反演结果都不能直接使用。我自己的习惯是反演完成后立刻输出这三个衍生指标作为第一道质检关卡。这个习惯在后面帮我避免了好几次拿着离谱结果自嗨的尴尬。4. Matlab代码实现从ODE求解器到HHO主循环下面是可以直接拼装起来的代码骨架。为了清晰我拆成四个模块按文件名和功能组织。4.1 第一步把SEIR写成ode45能调用的函数function dydt func_SEIR(t, y, p) % y [S; E; I; R] beta p(1); sigma p(2); gamma p(3); S y(1); E y(2); I y(3); R y(4); N S E I R; dydt [-beta*S*I/N; beta*S*I/N - sigma*E; sigma*E - gamma*I; gamma*I]; end这个函数的输入输出格式完全遵循ode45的要求p是参数向量y是状态向量t即使不使用也要占位。N从状态变量求和得到这样不用额外传总人口省掉一道参数传递。4.2 第二步目标函数里的参数与数据接口function err calObj(x, data, dataT) % x [beta, sigma, gamma, E0, I0] beta x(1); sigma x(2); gamma x(3); E0 x(4); I0 x(5); N data.N; S0 N - E0 - I0; y0 [S0; E0; I0; 0]; p [beta, sigma, gamma]; [~, Y] ode45((t,y) func_SEIR(t, y, p), dataT, y0); modelCum Y(:,3); % 模型累计感染曲线 modelInc [modelCum(1); diff(modelCum)]; % 模型每日新增 obsInc data.inc; obsCum data.cum; errCum sqrt(mean((modelCum - obsCum).^2)) / mean(obsCum); errInc sqrt(mean((modelInc - obsInc).^2)) / mean(obsInc); err 0.5 * errCum 0.5 * errInc; end这里有一点值得单独说明dataT直接传天数向量比如(0:120)ode45返回的Y行数就和dataT长度一致后续算diff、算RMSE都不需要对齐坐标。很多代码喜欢让ode45自动选tspan返回一大堆不规则的时间点然后再插值到观测时间纯属给自己找麻烦。4.3 第三步HHO主循环的公式落码这部分我把公式拆进一个函数包含Levy飞行的工程简化版。注意面向可读性我没有做向量化加速但逻辑和原论文一致。function [bestX, bestF, convCurve] HHO(popNum, maxIter, lb, ub, dim, fobj) % 初始化鹰群 X repmat(lb, popNum, 1) rand(popNum, dim) .* repmat(ub - lb, popNum, 1); F zeros(popNum, 1); for i 1:popNum F(i) fobj(X(i, :)); end [bestF, bestIdx] min(F); bestX X(bestIdx, :); convCurve zeros(maxIter, 1); for t 1:maxIter for i 1:popNum % 每个个体独立计算逃逸能量和跳跃强度 E0 2 * rand() - 1; E 2 * E0 * (1 - t / maxIter); J 2 * (1 - rand()); X_new X(i, :); if abs(E) 1 % ----- 探索阶段 ----- q rand(); if q 0.5 k randi(popNum); X_new X(k, :) - rand() * abs(X(k, :) - 2 * rand() * X(i, :)); else X_mean mean(X, 1); X_new (bestX - X_mean) - rand() * (lb rand(1, dim) .* (ub - lb)); end else % ----- 开发阶段 ----- DeltaX bestX - X(i, :); r rand(); if abs(E) 0.5 r 0.5 % 软围攻 X_new DeltaX - E * abs(J * bestX - X(i, :)); elseif abs(E) 0.5 r 0.5 % 硬围攻 X_new bestX - E * abs(DeltaX); elseif abs(E) 0.5 r 0.5 % 渐进式快速俯冲软围攻 Y bestX - E * abs(J * bestX - X(i, :)); if fobj(Y) F(i) X_new Y; else Z Y randn(1, dim) .* LevyFlight(dim); if fobj(Z) F(i) X_new Z; else X_new Y; end end else % 渐进式快速俯冲硬围攻 X_mean mean(X, 1); Y bestX - E * abs(J * bestX - X_mean); if fobj(Y) F(i) X_new Y; else Z Y randn(1, dim) .* LevyFlight(dim); if fobj(Z) F(i) X_new Z; else X_new Y; end end end end % 边界钳制 X_new min(max(X_new, lb), ub); F_new fobj(X_new); if F_new F(i) X(i, :) X_new; F(i) F_new; end if F_new bestF bestF F_new; bestX X_new; end end convCurve(t) bestF; end end function L LevyFlight(dim) beta 1.5; sigma_num gamma(1 beta) * sin(pi * beta / 2) / ... (gamma((1 beta) / 2) * beta * 2^((beta - 1) / 2)); u randn(1, dim) * sigma_num; v randn(1, dim); step u ./ abs(v).^(1 / beta); L 0.01 * step; end这里要特别说明一个经常被网上的简化代码搞错的点随机数r1、r2、q、r、J、E0必须在对每个个体更新位置时重新生成不能在整个迭代开始前一次性生成一个固定向量。HHO的探索能力和这些随机数的灵活性直接挂钩一旦随机数被固定成常数探索阶段就退化成一次确定性偏移算法性能明显下降。这个问题我在后面避坑章节会再展开。4.4 第四步主脚本整合跑通一个最小案例% 模拟数据生成 N 100000; betaTrue 0.30; sigmaTrue 0.2; gammaTrue 0.10; E0_true 0; I0_true 1; S0_true N - E0_true - I0_true; y0_true [S0_true; E0_true; I0_true; 0]; dataT (0:120); [~, Y_true] ode45((t,y) func_SEIR(t, y, [betaTrue, sigmaTrue, gammaTrue]), dataT, y0_true); obsCum Y_true(:,3) randn(size(dataT)) * 0.05 * max(Y_true(:,3)); obsCum max(obsCum, 0); obsInc [obsCum(1); diff(obsCum)]; data.N N; data.cum obsCum; data.inc obsInc; % 反演设置 lb [0.05, 0.05, 0.02, 0, 0.1]; ub [0.8, 1.0, 0.5, 50, 20]; dim 5; fobj (x) calObj(x, data, dataT); % 调用HHO rng(2024); % 固定随机种子保证结果可复现 [bestX, bestF, curve] HHO(15, 30, lb, ub, dim, fobj); % 结果处理 betaHat bestX(1); sigmaHat bestX(2); gammaHat bestX(3); R0_hat betaHat / gammaHat; fprintf(反演结果: beta%.3f, sigma%.3f, gamma%.3f, R0%.2f\n, ... betaHat, sigmaHat, gammaHat, R0_hat); % 画图对比 tspan dataT; [~, Yfit] ode45((t,y) func_SEIR(t, y, bestX(1:3)), tspan, ... [N - bestX(5) - bestX(4); bestX(4); bestX(5); 0]); figure(Position, [100 100 1200 400]); subplot(1,2,1); plot(tspan, obsCum, o, LineWidth, 1.2); hold on; plot(tspan, Yfit(:,3), -, LineWidth, 1.6); xlabel(天数); ylabel(累计感染); legend(观测数据, HHO拟合, Location, northwest); title(累计感染曲线拟合); grid on; subplot(1,2,2); plot(1:length(curve), curve, -, LineWidth, 1.5); xlabel(迭代次数); ylabel(目标函数值); title(HHO收敛曲线); grid on;这个最小案例里真实参数是β0.3、σ0.2、γ0.1观测数据带有5%的累计噪声种群15、迭代30代。跑完之后你应该能看到目标函数值从十几快速掉到零点几拟合曲线和观测数据基本重合反演参数接近真实值。5. 实验验证模拟数据反演与收敛性分析5.1 闭环实验设计提前知道答案再反推为什么一定要先做模拟数据闭环实验因为真实数据你永远不知道正确答案是什么算法调得再好也无法判断它是不是找到了合理的解。模拟数据反演相当于闭卷考试改成了开卷我先设定一组真实参数用ODE生成观测数据加噪声再用HHO去反推。能找回真实参数说明这套流程在逻辑上是通的才能放心去处理真实数据。这个闭环思路是我做所有反演类算法时雷打不动的第一步每次都让我省下大量盲目调试的时间。5.2 收敛曲线与参数估计精度用上面主脚本的设置跑5次不同随机种子结果趋势非常稳定前5代目标函数值快速下降10代左右基本接近稳定20代之后曲线平滑收敛。参数估计精度如下表参数真实值HHO反演均值5次运行标准差相对误差β0.300.3030.008~1%σ0.200.1950.021~2.5%γ0.100.1030.006~3%R03.002.950.12~1.7%这个结果说明2030代对于这个3参数加2初始条件的反演已经够用。事实上如果你只关心β和γ这两个直接决定R0把E0和I0固定住10代之内就能找到比较满意的位置。5.3 与传统方法的对比局部最优有多常见在同一组观测数据上我做了个简单对比方法设置反演结果目标值fminsearch初值β0.25,σ0.3,γ0.15β0.52,σ0.09,γ0.41局部最优偏离严重fminsearch初值β0.28,σ0.22,γ0.12靠天给的初值接近真实收敛PSO20粒子,50代β0.299,σ0.196,γ0.102接近真实HHO15鹰,30代β0.302,σ0.193,γ0.104接近真实这个对比不是要说HHO一定碾压其他算法——PSO在算够代数的前提下也能收敛。但注意到fminsearch那两次的结果了吗同样的算法、同样的数据、同样的目标函数仅仅因为初值不同一个彻底失败一个勉强成功。而HHO和PSO都不需要猜初值这是群体智能算法在SEIR反演这类问题上最有说服力的优势它们把对初值的依赖转换成了对迭代预算的依赖而后者是你可以控制的。6. 实战中真正会踩到的六个坑6.1 坑一E0和I0当自由参数拟合越好参数越离谱我早期犯过的最典型错误把E0和I0都设置成超宽边界比如E0在[0, 10000]I0在[0, 1000]。HHO确实找到了误差很小的解但反演出的E0是4872I0是58.3——初始潜伏者比初始感染者高两个数量级这会导致模型前期曲线严重变形为了拟合数据β和γ全被带偏。后来我查了一下这类问题在传染病反演文献里叫初始条件与参数的可辨识性耦合通俗讲就是参数和初值会互相补锅。解决办法很简单I0直接用观测首日的新增数或者累计数E0给一个合理范围的小区间通常设为I0的050倍就够必要时直接设成0。6.2 坑二目标函数对噪声数据过度敏感真实数据不像模拟数据那么干净。我最开始用纯RMSE做目标函数时某天单日数据异常偏高反演结果就为了迁就这个尖峰把σ调得很小、γ调得很大整个传播周期的形态被破坏。后来改成按均值归一化并引入平方根修正效果立竿见影。再一个经验是如果噪声实在太大先在数据层面做7日移动平均把短期波动消化掉再进目标函数。否则你再怎么调算法底层数据信噪比太低神仙也救不了。6.3 坑三ode45精度与寻优耗时的失衡HHO每评估一次目标函数就要完整跑一遍ode45种群15、迭代30意味着450次ODE求解。ode45默认的RelTol是1e-3这在大多数情况下够用但参数反演对精度的敏感度比想象中高——因为目标函数是误差的累积每个时间点的微小数值误差会叠加到目标值上扰乱算法的精细化搜索。建议用odeset把RelTol设为1e-6opts odeset(RelTol, 1e-6, AbsTol, 1e-8); [~, Y] ode45((t,y) func_SEIR(t, y, p), dataT, y0, opts);代价是求解时间可能增加一截。我的经验是两步走先用默认精度粗跑20代确定一个大致区域再用1e-6精度精跑20代。既不会等太久最终结果又稳定。6.4 坑四随机数复用让HHO变成伪随机搜索这个坑比较隐蔽。有些开源简化代码为了省事会在每次迭代开头一次性生成q、r1、r2、r3、E0这些随机数的数组循环里直接取用。看起来差别不大实际跑起来收敛速度会慢很多。原因在于HHO每个个体应该独立决定自己的搜索方向如果整个种群共享同一组随机数种群多样性就被压缩了探索阶段很多个体会被推往相似的方向相当于群体搜索退化成几次独立搜索。一定要把随机数的生成放到每个个体的位置更新内部确保独立。我在4.3节的代码里就是这么写的。6.5 坑五只给几周数据就做全局反演基本无解SEIR的动力学过程需要足够长的观察窗口才能把参数间的耦合信息暴露出来。如果数据只覆盖了疫情初期的一小段曲线还在指数上升阶段那么β、σ、γ之间的差异几乎无法从数据中区分——你会得到一组误差差不多但参数完全不同解的集合。实际经验是数据至少要覆盖到新增曲线的上升段、拐点和下降段中的两段少于这个信息量反演结果只能作为参考不能作为决策依据。这个坑没有算法能绕过去因为你是在和信息量的物理极限对抗。6.6 坑六跑到边界的最优参数可能是过拟合信号如果你发现最优参数的某个分量恰好落在边界上比如β刚好等于0.8这通常不是好消息。一个真正合理的参数组合在边界设定宽松的情况下应该落在参数空间内部落在边界意味着数据里可能存在你未纳入模型的结构性因素比如季节效应、干预措施、人群移动。HHO在边界上找到的最优解本质上是在搜索空间边缘强行找到的妥协解这时不要直接采信而是先回头检查模型假设是否还成立。我自己跑完这个项目之后标准流程基本固定了先用模拟数据做闭环验证确认HHO能把参数找回再上真实数据边界一定用流行病学常识压住宁可粗一点也不能完全放飞反演完成后必算R0和平均潜伏期这两个指标不过关其他都不用看。这套流程后来换了好几组不同的数据几乎没再翻过车。如果你也在用SEIR做反演建议把闭环验证、边界约束、衍生指标校验这三板斧立起来比盲目换优化算法有用得多。