企业官网建设流程全解析

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的四舱室结构和R0

SEIR模型把人按疾病状态分成四类：易感者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 = S+E+I+R是总人口，在封闭人群假设下保持不变。需要特别说明的是，我后面代码里变量的名字都避开了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输出的累计感染数，通常用I+R表示；C_data是实际观测到的累计感染数。但这只是第一层。我在实际项目中更推荐对累计感染先做开方再算误差，也就是：

MSE = mean( ( sqrt(C_model) - sqrt(C_data) )^2 )

为什么要开方？因为累计感染曲线是单调递增且往往跨度很大的。如果直接用原始值，后期的大数值会在误差里占据压倒性权重，前期拟合得再好也被忽略。开方相当于做了一次压缩变换，让前期和后期对误差的贡献更均衡。第三层考虑是数据口径问题。有些数据是“累计确诊”，有些是“每日新增”，有些是“现存阳性”。SEIR输出直接对应的是I+R，如果你手里的是每日新增，那就应该用新增 = sigma * E去对齐，而不是把I+R拿去对比。口径一旦错了，优化的参数再漂亮也是错的。这个细节我在给研究生改代码时几乎每次都要强调。

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, n+1); 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(:,k+1) = 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 end

Levy飞行的实现里有一个需要特别注意的地方，Matlab里gamma函数名和SEIR参数名撞车的问题我在前面提醒过，这里再提醒一次。另外在计算步长时一定要防止除零：

function L = levy_flight(dim) beta = 1.5; sigma_num = gamma(1+beta) * sin(pi*beta/2); sigma_den = gamma((1+beta)/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 参数边界与种群规模的选择

参数边界直接决定搜索空间大小，给得太宽浪费时间，给得太窄可能错过真解。我常用的经验值如下：

如果你的观测数据尺度完全不同，这些边界需要重新审视。种群规模我建议30左右就够。太小容易早熟，太大每代计算成本线性上涨，收益却有限。迭代次数视目标函数复杂度而定，SEIR这种3维参数问题200代基本够用，如果发现收敛曲线还在下降就加到500。跑完一次如果结果不好，不要急着改代码，先固定随机种子多跑几遍，然后选最优结果。

4. 案例测试与结果分析

4.1 构造一份带噪声的模拟数据

为了验证整个链路是否正确，我先用一组已知参数造一份“实测数据”。这样做的好处是心里有答案，可以检查优化器有没有找回真实的参数值。假设N=100万，E0=10，I0=5，真实参数取beta=0.35，sigma=1/6，gamma=1/12，对应潜伏期6天、传染期12天、R0约4.2。在这个参数下跑60天，取前30天作为历史数据，然后叠加一点高斯噪声，模拟真实上报数据的随机波动。

我在主脚本里写的就是这个过程。噪音幅度我习惯控制在累计感染数的1%-2%量级，太小就失去测试意义，太大则任何算法都难以恢复真值。用带噪数据反演出来的参数不可能和真实值完全一致，但只要优化后的曲线能贴合数据、参数落在合理区间，就说明这套框架是能用的。

4.2 读取收敛曲线并判断优化质量

HHO跑完之后，我会先看两个东西。第一个是收敛曲线，也就是每次迭代的最优适应度。正常情况下曲线应该快速下降然后趋于平缓，说明种群逐步锁定了最优区域。如果曲线中途出现断崖式下降，通常意味着某个个体发生了Levy跳变，越过了一个狭窄的“山谷”，这是HHO的正常现象，不用紧张。如果曲线一开始就平得像一潭死水，那大概率是初始化出了问题，所有个体挤在边界附近，或者目标函数返回了同一个超大值。

第二个是反推参数的合理性检验。假设优化结果是beta=0.36，sigma=0.168，gamma=0.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是一个性价比很高的默认选择。如果将来优化维度升到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 end

1e10这个惩罚值要足够大，大到优化器不可能选择它，但也不能是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 参数不可辨识与补偿效应

这是传染病反演里最难处理的一个问题。我举个例子：一组参数beta=0.3、gamma=1/10，另一组beta=0.45、gamma=1/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代码依然能用。多跑几组数据，你就知道“参数反演”这四个字背后有多少细节值得琢磨了。