1. 这个项目到底在解决什么问题1.1 参数才是传染病模型的命门做传染病建模的人大概率都经历过这种尴尬SEIR模型的方程背得滚瓜烂熟示意图画得漂漂亮亮但一碰到真实数据就哑火——因为 beta、sigma、gamma 这些参数根本不知道。哈里斯鹰算法HHO配合SEIR模型做参数优化解决的就是这个“数据到参数”的反问题。把HHO当作一个全局搜索器把SEIR当作一个评价器两者一组合就能从每日新增或累计感染数据里反推出模型参数全程只用Matlab就能跑通。这条技术路线其实适用面很宽。数模竞赛里常常有“根据历史数据估计传染病发展趋势”的题目科研工作中做疫情反演、传播动力学分析也需要估计参数还有不少本科生、研究生的课程设计直接就拿这个当题目。市面上的常规做法是先手动调参凑出一条“看起来还行”的曲线然后开始分析。问题在于SEIR是一个四变量耦联的非线性常微分方程组参数之间互相牵制手动调参不仅慢而且很容易被局部最优骗过去。HHO这类群体智能算法的价值就是把这些不可直接观测的参数从“人工瞎猜”变成“机器搜索”。1.2 从正向模拟到反向反演如果把SEIR当成一台前向机器输入参数输出感染曲线那参数优化做的就是把这台机器倒过来用。给定一条可观测的累计感染曲线反推出是哪组参数生成了它。这就是典型的反问题。反问题最大的难点有两个一是目标函数是非凸的传统梯度下降容易卡在局部极值二是变量之间有强耦合比如beta大一点、gamma也大一点最终曲线可能和另一组参数很接近。HHO的优势在于它根本不关心目标函数是否光滑、是否可导它只需要反复采样、比较适应度就能在参数空间里逐步逼近全球最优区域。用生活里的话说这就像校准一台台上的弹道。你不关心风速、湿度、弹道系数到底是多少只要不断试射、对比落点、调整瞄准最终找到一组“打得很准”的参数就行。HHO干的活就是自动做这个试射和调整过程而且它不是单发试射是一群“鹰”同时在不同方向试射然后共享情报越试越准。1.3 这套框架能迁移到哪里这套HHO-SEIR框架并不仅仅能用来做SEIR。你把seir_rhs函数替换成SIR、SEIRD、SIER或者带时变接触率的扩展模型优化器那一部分完全不用改。甚至不限于传染病模型只要是“给定一个黑箱模拟器想通过历史观测数据反推内部参数”的问题HHO都能套进去。比如小范围的舆情传播模型、商品扩散模型、甚至简单的机械系统参数辨识思路都是一样。这也是我特别推荐大家把HHO代码吃透的原因你学到的不是一套孤立代码而是一个可复用的“模拟器-优化器”耦合范式。项目里带着完整Matlab代码看懂一遍后面改模型就是改个方程的事。2. SEIR模型与HHO算法的核心原理2.1 SEIR的四舱室结构和R0SEIR模型把人按疾病状态分成四类易感者S、暴露者E、感染者I、恢复者R。暴露者是那些已经被传染、但还处于潜伏期、没有表现出症状的人。标准形式的微分方程如下dS/dt -beta * S * I / N dE/dt beta * S * I / N - sigma * E dI/dt sigma * E - gamma * I dR/dt gamma * I其中beta是有效接触率也就是一个感染者每天能传染多少易感者的比例sigma是潜伏期向感染期转化的速率通常写成 1/平均潜伏期gamma是恢复速率通常写成 1/平均传染期。这组方程里N SEIR是总人口在封闭人群假设下保持不变。需要特别说明的是我后面代码里变量的名字都避开了Matlab内置函数gamma不然后面调用伽马函数时会被同名变量覆盖这种低级错误在实操中特别常见。基本再生数R0是传染病模型最常用的一个输出指标。在标准的SEIR框架下R0 beta / gamma含义是一个感染者进入完全易感人群之后平均能传染几个人。这个数字实际上可以当作目标函数的一部分来验证结果比如你反推出来的R0是否落在医学常识范围内。如果R0大于10或者小于0.5大概率是优化走到了错误的参数区域该检查边界设置或者数据口径了。2.2 目标函数设计的三个层次参数优化的质量完全取决于目标函数怎么定义。最简单粗暴的形式是均方误差MSE mean( (C_model(t) - C_data(t))^2 )其中C_model是SEIR输出的累计感染数通常用IR表示C_data是实际观测到的累计感染数。但这只是第一层。我在实际项目中更推荐对累计感染先做开方再算误差也就是MSE mean( ( sqrt(C_model) - sqrt(C_data) )^2 )为什么要开方因为累计感染曲线是单调递增且往往跨度很大的。如果直接用原始值后期的大数值会在误差里占据压倒性权重前期拟合得再好也被忽略。开方相当于做了一次压缩变换让前期和后期对误差的贡献更均衡。第三层考虑是数据口径问题。有些数据是“累计确诊”有些是“每日新增”有些是“现存阳性”。SEIR输出直接对应的是IR如果你手里的是每日新增那就应该用新增 sigma * E去对齐而不是把IR拿去对比。口径一旦错了优化的参数再漂亮也是错的。这个细节我在给研究生改代码时几乎每次都要强调。2.3 哈里斯鹰算法的探索与开发机制HHO算法的灵感来自哈里斯鹰捕猎兔子的行为。鹰群先在全场不同位置搜索猎物发现目标之后根据兔子的逃跑能量决定是继续包围还是发起俯冲。算法的位置更新分三个大阶段探索阶段、探索向开发转换阶段、开发阶段。在探索阶段个体位置通过两种随机策略更新。一种是完全随机跳到种群中另一个个体的附近另一种是围绕当前全局最优位置和种群平均位置生成新个体。这个阶段的核心是保证搜索范围足够广不容易陷入局部最优。随后算法计算“逃跑能量”E公式是E 2 * E0 * (1 - t / T)E0是每次迭代在-1到1之间重新取值的随机数。当|E|大于等于1说明兔子体力充足鹰群继续大范围探索当|E|小于1鹰群转入开发也就是围绕最优解精细搜刮。开发阶段又根据两个随机数r和E分成四种围攻策略软围攻、硬围攻、带渐进式快速俯冲的软围攻、带渐进式快速俯冲的硬围攻。其中Levy飞行用于模拟鹰的俯冲路径它能生成偶尔跳得很远的随机步长让算法在局部搜索时仍然保留一部分跳出能力。HHO算法对使用者的友好之处在于核心超参数只有种群规模和迭代次数没有交叉概率、惯性权重之类需要反复调的东西。只要设置好参数边界和适应度函数剩下的交给算法自己跑就行。这一点对比遗传算法和粒子群算法来说确实省心不少。3. Matlab代码实现与细节拆解3.1 工程文件结构与数据准备整个项目我拆成五个部分主脚本main_HHO_SEIR.m、优化器HHO.m、Levy飞行函数levy_flight.m、SEIR求解器seir_rk4.m和SEIR右端函数seir_rhs.m。主脚本负责读数据、设置边界、调用优化器、画图优化器是通用的换任何目标函数都能直接用SEIR这部分就是前向模拟器。如果你把这个框架迁移到其他模型只需要改seir_rhs.m这一个文件HHO完全不用动这是这个结构最值钱的地方。数据准备阶段需要确定总人口N、初始状态S0、E0、I0、R0和观测周期Tmax。在大多数课程设计和论文场景里初始感染人数是未知的但通常只占人口极小比例设为个位数到几十人都可以。如果你想连初始感染人数一起优化把I0加入theta向量目标函数维度就从3变成4。代码层面要做一个小处理在目标函数内部把I0取整并限制不小于1否则一个分数感染者在模型里看起来会很别扭。3.2 用RK4求解SEIR而不是直接用ODE45很多人在Matlab里解微分方程第一反应是ode45。但在这个项目里我更推荐自己写固定步长的RK4。原因很实际HHO每次迭代要评估种群内所有个体几十个个体乘几百代迭代意味着SEIR求解器要被调用上千次。ode45虽然有自适应步长控制但每次调用都有额外开销而且如果参数在探索期跑飞了ode45可能为了满足容差把步长压得很小导致一次评估异常地慢。固定步长RK4只要步长足够小精度对这类流行病模型完全够用而且耗时稳定可控。下面是一个干净可用的SEIR求解函数function [C, S, E, I, R] seir_rk4(theta, Tmax, N, S0, E0, I0, R0) beta_r theta(1); sigma theta(2); gamma_r theta(3); h 1; % 步长为1天如需更高精度可改为0.1 n round(Tmax / h); y zeros(4, n1); y(:,1) [S0; E0; I0; R0]; for k 1:n f1 seir_rhs(y(:,k), beta_r, sigma, gamma_r, N); f2 seir_rhs(y(:,k) h/2*f1, beta_r, sigma, gamma_r, N); f3 seir_rhs(y(:,k) h/2*f2, beta_r, sigma, gamma_r, N); f4 seir_rhs(y(:,k) h*f3, beta_r, sigma, gamma_r, N); y(:,k1) y(:,k) h/6 * (f1 2*f2 2*f3 f4); end S y(1,:); E y(2,:); I y(3,:); R y(4,:); C I R; % 累计感染口径现症感染 已恢复 end function dydt seir_rhs(y, beta_r, sigma, gamma_r, N) S max(y(1), 0); % 负值保护防止数值溢出 E max(y(2), 0); I max(y(3), 0); dS -beta_r * S * I / N; dE beta_r * S * I / N - sigma * E; dI sigma * E - gamma_r * I; dR gamma_r * I; dydt [dS; dE; dI; dR]; end负值保护这一行一定要写。HHO在探索期什么参数都可能试出来一旦beta过大或者S0很小数值求解时S可能变成负值负的S进到betaSI/N里会引发连锁错误最终结果就是NaN。有了max截断即使参数不合实际模型也能返回一个合理的有限数值只不过误差很大优化器会自然地避开这些区域。3.3 HHO优化器核心循环这段代码是整个项目的发动机我把它做了精简保留了每个策略的核心更新式function [rabbit, rabbit_f, convergence] HHO(pop_size, Tmax_iter, dim, lb, ub, fun) X repmat(lb, pop_size, 1) rand(pop_size, dim) .* repmat(ub-lb, pop_size, 1); rabbit X(1,:); rabbit_f inf; for i 1:pop_size fit_i fun(X(i,:)); if fit_i rabbit_f rabbit_f fit_i; rabbit X(i,:); end end convergence zeros(Tmax_iter, 1); for t 1:Tmax_iter E0 2 * rand - 1; E 2 * E0 * (1 - t / Tmax_iter); for i 1:pop_size if abs(E) 1 % 探索阶段 q rand; if q 0.5 X_rand X(randi(pop_size), :); X(i,:) X_rand - rand * abs(X_rand - 2*rand*X(i,:)); else X_m mean(X, 1); X(i,:) (rabbit - X_m) - rand*(lb rand*(ub - lb)); end else % 开发阶段 r rand; J 2 * (1 - rand); deltaX rabbit - X(i,:); if r 0.5 abs(E) 0.5 % 软围攻 X(i,:) deltaX - E * abs(J*rabbit - X(i,:)); elseif r 0.5 abs(E) 0.5 % 硬围攻 X(i,:) rabbit - E * abs(deltaX); elseif r 0.5 abs(E) 0.5 % 软围攻渐进俯冲 Y rabbit - E * abs(J*rabbit - X(i,:)); Z Y rand(1,dim) .* levy_flight(dim); if fun(Y) fun(X(i,:)), X(i,:) Y; end if fun(Z) fun(X(i,:)), X(i,:) Z; end else % 硬围攻渐进俯冲 X_m mean(X, 1); Y rabbit - E * abs(J*rabbit - X_m); Z Y rand(1,dim) .* levy_flight(dim); if fun(Y) fun(X(i,:)), X(i,:) Y; end if fun(Z) fun(X(i,:)), X(i,:) Z; end end end % 边界截断 X(i,:) max(min(X(i,:), ub), lb); fit_new fun(X(i,:)); if fit_new rabbit_f rabbit_f fit_new; rabbit X(i,:); end end convergence(t) rabbit_f; end endLevy飞行的实现里有一个需要特别注意的地方Matlab里gamma函数名和SEIR参数名撞车的问题我在前面提醒过这里再提醒一次。另外在计算步长时一定要防止除零function L levy_flight(dim) beta 1.5; sigma_num gamma(1beta) * sin(pi*beta/2); sigma_den gamma((1beta)/2) * beta * 2^((beta-1)/2); sigma (sigma_num / sigma_den)^(1/beta); u randn(1, dim) * sigma; v randn(1, dim); L 0.01 * u ./ (abs(v).^(1/beta) eps); end这里加eps就是为了防止v恰好为0时产生NaN。虽然概率很低但在上千次目标函数评估里一次NaN就可能导致整个优化结果报废所以这种防御性写法值得养成习惯。3.4 参数边界与种群规模的选择参数边界直接决定搜索空间大小给得太宽浪费时间给得太窄可能错过真解。我常用的经验值如下参数含义常见合理范围说明beta有效接触率[0.05, 1]取决于接触频率传染性强的场景取0.3-0.6sigma潜伏期转化速率[0.05, 0.5]对应平均潜伏期2-20天gamma恢复速率[0.02, 0.3]对应平均传染期约3-50天如果你的观测数据尺度完全不同这些边界需要重新审视。种群规模我建议30左右就够。太小容易早熟太大每代计算成本线性上涨收益却有限。迭代次数视目标函数复杂度而定SEIR这种3维参数问题200代基本够用如果发现收敛曲线还在下降就加到500。跑完一次如果结果不好不要急着改代码先固定随机种子多跑几遍然后选最优结果。4. 案例测试与结果分析4.1 构造一份带噪声的模拟数据为了验证整个链路是否正确我先用一组已知参数造一份“实测数据”。这样做的好处是心里有答案可以检查优化器有没有找回真实的参数值。假设N100万E010I05真实参数取beta0.35sigma1/6gamma1/12对应潜伏期6天、传染期12天、R0约4.2。在这个参数下跑60天取前30天作为历史数据然后叠加一点高斯噪声模拟真实上报数据的随机波动。我在主脚本里写的就是这个过程。噪音幅度我习惯控制在累计感染数的1%-2%量级太小就失去测试意义太大则任何算法都难以恢复真值。用带噪数据反演出来的参数不可能和真实值完全一致但只要优化后的曲线能贴合数据、参数落在合理区间就说明这套框架是能用的。4.2 读取收敛曲线并判断优化质量HHO跑完之后我会先看两个东西。第一个是收敛曲线也就是每次迭代的最优适应度。正常情况下曲线应该快速下降然后趋于平缓说明种群逐步锁定了最优区域。如果曲线中途出现断崖式下降通常意味着某个个体发生了Levy跳变越过了一个狭窄的“山谷”这是HHO的正常现象不用紧张。如果曲线一开始就平得像一潭死水那大概率是初始化出了问题所有个体挤在边界附近或者目标函数返回了同一个超大值。第二个是反推参数的合理性检验。假设优化结果是beta0.36sigma0.168gamma0.082和真实值0.35、0.1667、0.0833非常接近说明管道是通的。在实际数据场景下没有“真实值”可以参考这时就要靠R0 beta/gamma是否在合理范围、潜伏期1/sigma和传染期1/gamma是否符合医学常识来判断。如果反推出平均潜伏期只有半天那要么数据有误要么模型结构不对参数再好看也不能直接采用。4.3 与其他优化算法的对比体会为了确认HHO不是“唯一解”而是“合适解”我把同一份数据和同一个目标函数分别丢给遗传算法GA和粒子群算法PSO。GA用Matlab自带的ga函数PSO自己写了三十行。三种算法都能找到相近的适应度值但体验差别明显。GA的问题是参数多需要设置交叉比例、变异比例、精英保留数调起来繁琐PSO的问题是对惯性权重敏感跑几次结果波动大。HHO在低维问题上的优势主要是稳定且实现简单。表格对比一下算法需要调节的关键超参数代码行数我的使用感受HHO种群数、迭代数约80行默认参数就好用稳定GA交叉率、变异率、种群数自带/约100行也能收敛但调参费时PSO惯性权重、c1、c2约60行轻量但容易早熟需多次运行我最终的结论很明确三维参数反演问题上HHO是一个性价比很高的默认选择。如果将来优化维度升到10维以上HHO也需要和其他算法做对比不能盲目自信。5. 实操中踩过的坑与排查方法5.1 目标函数返回NaN整个优化直接崩溃这是新手最容易撞上的问题。现象是HHO跑几步之后收敛曲线突然变成NaN或者优化结果全部是NaN。根因通常出在SEIR数值求解上参数太极端S被算成负数然后E中出现NaN最终累计感染C全是NaN。应对方式有两层第一层在seir_rhs里对S、E、I做max截断保证状态量非负第二层在目标函数里做最终兜底function mse obj_fun(theta, data, N, S0, E0, I0, R0, Tmax) [C, ~, ~, ~, ~] seir_rk4(theta, Tmax, N, S0, E0, I0, R0); C C(1:length(data)); mse mean((sqrt(C) - sqrt(data(:))).^2); if ~isfinite(mse) mse 1e10; end end1e10这个惩罚值要足够大大到优化器不可能选择它但也不能是inf因为inf在某些比较逻辑里会有奇怪的传播。实测下来这个兜底策略非常有效。5.2 收敛到不合理的参数区域有时候HHO收敛很快适应度很低但beta和gamma同时很大R0高达20明显不符合常理。这其实不是HHO的错而是目标函数没有提供足够的约束信息。累计感染曲线只约束了beta和gamma的比值关系很难独立约束两者。解决思路有两个方向。第一个方向是缩边界。比如根据疾病常识把gamma限制在[0.02, 0.3]以内sigma限制在[0.05, 0.5]以内医学先验就能直接把不合理区域排除掉。第二个方向是改数据口径。如果手头除了累计感染还有每日新增数据把新增量也放进目标函数相当于多了一条观测通道beta和gamma的“分工”会更明确。5.3 参数不可辨识与补偿效应这是传染病反演里最难处理的一个问题。我举个例子一组参数beta0.3、gamma1/10另一组beta0.45、gamma1/15理论上两者的R0都是3生成的累计感染曲线在前期可能非常接近。这种“不同参数产生相似曲线”的现象就叫不可辨识性。HHO可能在这两组参数之间来回跳最终收敛到哪一组取决于噪声扰动和初始化位置。应对措施是加先验正则化。在目标函数里增加一个惩罚项比如如果潜伏期偏离5-7天就加上一个额外误差或者如果R0超出2-6就惩罚。这个惩罚项的权重不必很大它的作用只是“打破对称性”让优化器偏向更合理的解。这在实际项目中非常有用尤其是模型要用于政策评估时一组结构合理的参数比一组拟合误差略小但违背常识的参数重要得多。5.4 优化速度太慢怎么办当观测周期很长、RK4步长又设得很小时一次目标函数评估可能要算几百步一小时跑不完。我常用的加速手段有三个。第一把RK4步长从0.1天改成1天。对日粒度的累计感染数据来说1天的步长精度足够计算量减少到十分之一。第二缩短观测周期进行预跑。先用前15天数据快速跑一遍看HHO是否正常工作、参数是否大致合理再放到全周期。第三如果机器有多核把HHO内部循环改成parfor。因为每个个体的评估完全独立不存在数据竞争可以直接并行四核机器大约能快三倍。还有一个容易被忽略的小技巧提前把真实数据转换成行向量。Matlab里行向量和列向量的区别经常让人在mean和plot时踩坑目标函数里用data(:)强制转成行向量可以省掉很多莫名其妙的维度报错。6. 我的几点私人心得最后说几个我自己的习惯。第一任何群智能算法项目第一行就要写rng(固定种子)。HHO带随机性不固定种子每次跑出来的结果都不同写报告、做对比都会头疼。种子一旦固定结果就完全可复现评审和答辩时底气完全不一样。第二我强烈建议对优化变量做对数变换。直接优化theta的原始值时参数分布跨越几个数量级算法搜索效率不高改成优化x log(theta)目标函数里用exp(x)转回真实参数搜索空间更均匀HHO的边界问题也少很多。第三不要只跑一次HHO就下结论。我的做法是连续跑20次记录每次的最优参数和适应度然后看参数的分布区间。如果20次结果高度集中说明这个问题可辨识性较好如果结果满天飞说明目标函数约束不足应该回去检查模型结构而不是继续调HHO参数。这套HHO-SEIR框架最值得保存的地方就是“优化器与模拟器解耦”的写法。今天你拿它反演SEIR的参数明天换一个SIR模型或者SEIRD模型只需要改一个右端函数后天换一个完全不同的领域任何需要黑箱参数反演的问题这段HHO代码依然能用。多跑几组数据你就知道“参数反演”这四个字背后有多少细节值得琢磨了。