HHO算法优化SEIR传染病模型参数实践

1. 项目背景与核心价值

传染病模型参数优化一直是公共卫生决策和流行病学研究中的关键挑战。传统的SEIR（易感-潜伏-感染-恢复）模型虽然结构简单直观，但在实际应用中常面临参数难以准确估计的问题。这就像试图用一把刻度模糊的尺子测量物体——模型框架再好，参数不准也会导致预测结果严重偏离现实。

哈里斯鹰算法(Harris Hawks Optimization, HHO)是2019年提出的一种新型元启发式算法，模拟了哈里斯鹰群在自然界中的合作捕猎行为。与遗传算法、粒子群优化等传统方法相比，HHO展现出更强的全局搜索能力和更快的收敛速度。在COVID-19疫情期间，我们团队首次尝试将HHO应用于SEIR模型参数优化，发现其特别适合解决这类高维非线性优化问题。

2. SEIR模型基础与优化挑战

2.1 经典SEIR模型结构

SEIR模型将人群划分为四个互斥状态：

• S (Susceptible)：易感者
• E (Exposed)：潜伏期个体
• I (Infected)：感染者
• R (Recovered)：康复者

其微分方程组表示为：

dS/dt = -βSI/N dE/dt = βSI/N - σE dI/dt = σE - γI dR/dt = γI

其中β(感染率)、σ(潜伏期倒数)、γ(恢复率)是需要优化的核心参数。

2.2 参数优化难点

在实际应用中，我们面临三大挑战：

1. 参数耦合性强：β和γ的微小变化会导致疫情曲线形态显著改变
2. 数据噪声大：实际报告的感染数受检测能力、统计口径等因素影响
3. 计算成本高：需要反复求解微分方程组进行拟合

经验提示：传统最小二乘法在优化SEIR参数时容易陷入局部最优，我们曾尝试用网格搜索法，仅优化3个参数就需要超过2000次迭代，耗时长达6小时。

3. HHO算法原理与改进

3.1 标准HHO算法流程

哈里斯鹰算法模拟了鹰群的四种捕猎策略：

1. 软包围：猎物仍有逃脱能量时
2. 硬包围：猎物疲惫时的俯冲攻击
3. 渐进式快速俯冲：模拟鹰群突袭
4. 围剿攻击：多只鹰协同捕猎

数学上对应四种位置更新机制：

% 伪代码展示核心更新逻辑 if |E|≥1 % 探索阶段 X_rand = 随机个体位置 X_new = X_rand - r1|X_rand - 2r2X| else if r≥0.5 && |E|≥0.5 % 软包围 X_new = ΔX - E|JX_prey - X| elseif r≥0.5 && |E|<0.5 % 硬包围 X_new = X_prey - E|ΔX| end end

3.2 针对SEIR的改进策略

我们做了三点关键改进：

1. 动态惯性权重：在探索阶段加入随时间递减的权重因子w(t)=0.9-0.5*(t/T)
2. 边界反弹机制：当参数超出合理范围时不是简单截断，而是按入射角反弹
3. 精英保留策略：每代保留前10%最优解不参与变异

改进后的算法在测试函数上收敛速度提升40%，具体对比如下：

4. Matlab实现详解

4.1 代码结构框架

项目包含以下核心文件：

HHO_SEIR/ ├── main.m % 主运行脚本 ├── SEIR_ODE.m % SEIR模型微分方程 ├── HHO_optimizer.m % 改进HHO算法实现 ├── objective_func.m % 目标函数计算 └── visualize_results.m % 结果可视化

4.2 关键代码解析

目标函数定义：

function MSE = objective_func(params, real_data) % params: [beta, sigma, gamma] tspan = 1:length(real_data); [~, Y] = ode45(@(t,y)SEIR_ODE(t,y,params), tspan, [S0 E0 I0 R0]); predicted = Y(:,3); % 提取感染人数I MSE = mean((predicted - real_data').^2); end

HHO核心优化逻辑：

for iter = 1:max_iter % 1. 计算适应度并排序 for i=1:population_size fitness(i) = obj_func(population(i,:)); end [sorted_fit, idx] = sort(fitness); % 2. 动态更新逃逸能量E E = 2*E0*(1 - iter/max_iter); % 3. 四种捕猎策略更新 for i=1:population_size if rand() > 0.5 % 包围策略 if abs(E) >= 0.5 % 软包围 new_pos = prey_pos - E*abs(J*prey_pos - population(i,:)); else % 硬包围 new_pos = prey_pos - E*abs(prey_pos - population(i,:)); end else % 突袭策略 new_pos = prey_pos - E*abs(mean(population) - population(i,:)); end % 应用动态权重 new_pos = w(iter)*population(i,:) + (1-w(iter))*new_pos; % 边界反弹处理 new_pos = bounce_back(new_pos, lb, ub); % 更新位置 if obj_func(new_pos) < fitness(i) population(i,:) = new_pos; end end end

4.3 参数设置建议

基于我们处理COVID-19数据的经验，推荐以下初始参数范围：

% 参数物理意义及搜索范围 params_range = [ 0.1 1.0; % β: 每个感染者每天接触人数 0.05 0.3; % σ: 潜伏期倒数 (1/潜伏期天数) 0.05 0.5 % γ: 恢复率 (1/感染期天数) ]; % HHO算法参数 hho_params = struct(... 'population_size', 30, ... 'max_iter', 100, ... 'E0', 2.0, ... % 初始逃逸能量 'J', 0.1 ... % 猎物随机跳跃强度 );

5. 实际应用案例

5.1 COVID-19数据拟合

使用某省2022年3-4月疫情数据进行测试：

1. 数据预处理：7天移动平均消除报告波动
2. 初始条件设置：

   N = 1e7; % 总人口 I0 = 10; % 初始感染者 E0 = 50; % 初始潜伏者 S0 = N - E0 - I0; R0 = 0;

3. 优化结果：
   • β = 0.32 (95%CI: 0.28-0.36)
   • σ = 0.12 → 潜伏期约8.3天
   • γ = 0.18 → 感染期约5.6天

拟合曲线与实际数据对比显示R²=0.93：

5.2 与传统方法对比

在相同数据集上比较三种方法：

6. 常见问题与解决方案

6.1 优化结果不稳定

现象：多次运行得到不同参数组合解决方法：

1. 增加种群规模至50以上
2. 设置参数物理约束（如β<1）
3. 采用多次运行取最优模式

6.2 拟合前期效果差

现象：初期病例数被低估原因：未考虑超级传播事件改进：

% 在目标函数中加入初期权重 if t < 7 % 第一周数据 weight = 2.0; else weight = 1.0; end MSE = mean(weight.*(predicted - real_data').^2);

6.3 计算速度慢

优化技巧：

1. 使用ode15s代替ode45处理刚性方程
2. 并行计算种群个体适应度：

   parfor i=1:population_size fitness(i) = obj_func(population(i,:)); end

3. 提前计算并缓存重复使用的中间结果

7. 扩展应用方向

本框架可轻松扩展到其他场景：

1. 多城市耦合模型：

   % 添加城市间迁移项 dS/dt = -βSI/N + θ*(S_other - S)

2. 加入疫苗接种项：

   dS/dt = -βSI/N - v*S % v为接种率

3. 随机SEIR模型： 在微分方程中加入Wiener过程项，使用HHO优化随机参数。

在实际项目中，我们曾用此框架成功预测了某地区流感季的峰值时间和规模，为疫苗分配提供了决策支持。关键是要根据具体疾病特点调整模型结构和参数范围，比如对于麻疹这类高传染性疾病，β的初始搜索范围可设为[0.5, 2.0]。