★ FEATURED ARTICLE 做多变量回归预测RBF神经网络是个绕不开的经典选项结构简单、逼近能力强但中心值、宽度和连接权值怎么定一直是让人头疼的问题。前段时间我把鲸鱼优化算法WOA和RBF揉在一起写了个Matlab程序用来做多输入单输出的回归预测核心思路就是把RBF的宽度、中心值和连接权值全部当作待优化变量用WOA去全局搜索。这篇文章会把思路、数据预处理、Matlab代码实现、实验验证和调试经验完整拆开讲。适合刚接触智能优化算法和神经网络的新手也适合老手拿去做对比实验。1. 从RBF的痛点聊起为什么非要选WOA1.1 RBF网络的三个关键参数到底是什么RBF神经网络的结构很简单输入层、隐含层、输出层各一层。隐含层每个节点是一个径向基函数最常用的是高斯函数[ \varphi_i(x) \exp\left(-\frac{|x-c_i|^2}{2\sigma_i^2}\right) ]网络的最终输出是隐含层响应的线性加权和[ y \sum_{i1}^{H} w_i \varphi_i(x) ]H是隐含节点数。这里面需要确定的三组参数就是 (c_i)、(\sigma_i)、(w_i)。中心值 (c_i) 是一个和输入维度相同的向量表示这个节点在输入空间中“最敏感”的位置宽度 (\sigma_i) 控制高斯函数的扩散范围宽度小的时候只有输入离中心很近才会被激活宽度大的时候几乎所有样本都会产生明显响应权值 (w_i) 则决定每个节点对最终预测的贡献大小。用一个生活化的类比每个隐含节点像一盏小手电筒中心值决定手电筒照向哪里宽度决定光斑的扩散半径权值决定这盏灯的亮度。多盏手电筒叠在一起就拼出了最终预测曲面。中心选得不准手电筒照错了地方宽度取得太大光斑糊了一片权值设得不对亮度分配失衡。这三者相互耦合这就是为什么RBF参数不好调。1.2 传统方法为什么容易翻车传统RBF训练分两步。第一步用K-means聚类确定中心值然后用聚类中心之间的距离或最大距离来估算宽度第二步用最小二乘法求输出权值。K-means本身对初始聚类中心敏感聚类数目的选择也很随意最关键的是聚类中心是为“数据密度”服务的它并不关心“预测误差”在哪里所以聚类效果看着不错回归误差却可能很高。如果改用梯度下降联合训练所有参数又会出现新的问题。RBF的局部响应特性让目标函数出现很多平坦区域和局部极小值梯度下降很容易陷进去而且学习率稍大一点就容易震荡稍小一点又收敛极慢。再加上中心值、宽度、权值三种参数的尺度完全不同用一个学习率去更新它们并不公平。WOA属于群体智能优化算法不依赖梯度只依赖适应度函数的输出。它在参数空间里用“包围—气泡攻击—随机搜索”三种策略来回寻找更优位置天然适合这种参数耦合强、目标函数不光滑的问题。而且WOA要调整的超参数很少不像遗传算法还需要设交叉率变异率也不像粒子群需要仔细调惯性权重这在写代码时非常省心。2. 多输入单输出回归的场景与数据准备2.1 多输入单输出问题怎么建模“多输入单输出”是回归预测里最常见的设定。比如输入是温度、湿度、风速输出是光伏功率输入是几个工艺参数输出是产品质量指标。用矩阵表示就是X是N行m列N是样本数m是特征数y是N行1列。每行样本代表一次观测一个m维特征向量映射到一个目标值。在RBF模型里输入维度m直接决定中心值向量的维度。设隐含节点数为H那么中心矩阵C的尺寸就是H乘m每一行是一个中心点。宽度向量sigma是1乘H权值向量w是H乘1。这样编码整个RBF参数时个体长度就是H*m H H。比如m3、H10编码长度就是50。这里有个容易踩的坑Matlab里矩阵是列优先存储reshape的时候如果不小心很容易把中心矩阵给排错。我习惯先把X组织成“样本数乘特征数”中心矩阵则用reshape(center_segment, H, m)这样每一行对应一个中心点。如果换成“特征数乘样本数”的布局所有索引公式都得跟着改非常容易出错所以强烈建议从头到尾统一成行样本布局。2.2 归一化不是可选项是必选项RBF的高斯函数依赖样本与中心之间的欧氏距离。如果输入特征之间量级差别很大距离计算就会被大量级特征主导小量级特征几乎不起作用。比如一个特征在0.01~0.1之间波动另一个在1000~5000之间前者对距离的贡献基本可以忽略模型等于少了一个输入。因此必须先做归一化让所有特征处于相同尺度我一般归一化到[0,1]区间。输出y同样要归一化。如果输出是几千的大数值权值如果只在[-1,1]内搜索很难产生足够的幅值来拟合输出。归一化到[0,1]或[-1,1]之后权值边界就可以设得比较小搜索也会更顺利。训练结束后把预测值反归一化还原再计算RMSE、MAE这些指标用户看到的仍然是原始量纲的结果。这里要特别注意归一化参数必须只从训练集计算。也就是说先用训练集求出min_x、max_x、min_y、max_y再用同一组参数去变换测试集。如果提前把测试集也混进来算min和max就是信息泄漏。虽然离线实验时测试集是已知的但模型上线后遇到的都是“未来数据”不能依赖未来数据的统计量。这个细节很多人忽略但它是数据科学里非常基础也非常重要的一条规则。3. WOA-RBF核心算法实现逐段拆解3.1 编码设计与种群初始化先把RBF的三组参数编码成一个一维向量每个鲸鱼个体就是这个向量。设输入维度m隐含节点数H编码长度D Hm H H。前Hm个元素是中心值接下来H个元素是宽度最后H个元素是权值。解码时依次切段再做reshape得到中心矩阵。边界设置我建议这样中心值的下界和上界直接取归一化后输入数据的最小值和最大值。如果数据归一化到[0,1]中心边界就是[0,1]如果数据归一化到[-1,1]中心边界就是[-1,1]。宽度边界建议[0.1,2]太小容易出现数值溢出太大则所有高斯响应都趋于相同。权值边界建议[-1,1]配合归一化后的输出区间足够用。初始化时用rand生成[0,1]之间的随机数再通过“下界 随机数*(上界-下界)”的方式映射到每个参数的区间。注意必须分段映射不能对整段编码统一映射因为中心、宽度、权值的上下界不同。3.2 WOA的三类位置更新策略WOA的每次迭代中每个个体根据概率p选择不同的更新公式。当p小于0.5时走“包围/搜索”线路否则走“螺旋攻击”线路在“包围/搜索”线路里又根据A的绝对值决定是收缩包围还是随机游走。具体来说a从2线性衰减到0A是2ar1-a所以A的初始范围大致在[-2,2]。当|A|1个体向一个随机个体靠近这种随机性保证了全局探索当|A|1个体向当前最优个体靠近完成局部开发。当p0.5时个体沿对数螺旋方式逼近最优个体这一步让算法有很强的局部精搜能力。经常有人问我为什么WOA里随机搜索和包围策略都要用C这个随机向量C的作用是模拟捕食过程中的随机扰动它让移动距离不是固定的避免个体直接跳死在某个位置变相增加搜索多样性。这个细节在代码里很容易被忽略但如果C固定为1收敛性能会差一些。3.3 适应度函数把RBF变成黑箱评估器适应度函数决定WOA往哪个方向搜索。我选的是训练集上的MSE公式是MSE mean((y_pred - y_true).^2)每次评估时先从个体解码出中心、宽度、权值然后计算隐含层输出矩阵Phi再与权值相乘得到预测值。关键是用矩阵化计算代替循环。伪代码思路function mse obj_fun(individual, X, y, H, m) % 解码 C reshape(individual(1:H*m), H, m); sigma individual(H*m1 : H*mH); w individual(H*mH1 : end); % 计算距离矩阵 NxH D2 zeros(size(X,1), H); for j 1:H diff X - repmat(C(j,:), size(X,1), 1); D2(:,j) sum(diff.^2, 2); end Phi exp(-D2 ./ (2 * sigma.^2 eps)); pred Phi * w(:); mse mean((pred - y).^2); end这里的sigma是1乘HMatlab会自动把D2的每一列除以对应的sigma平方。加上eps可以防止sigma为0时除零报错。如果你不想写for循环也可以用pdist2(X, C)直接计算距离矩阵再平方D2 pdist2(X, C).^2;两种写法结果一样pdist2更简洁但需要统计工具箱。在作业或者论文里手动写for循环反而更容易说明白原理。3.4 边界约束与越界处理WOA更新完位置后必然会有一些维度越过边界。最常用的处理策略是截断法小于下界的直接拉回下界大于上界的直接拉回上界。这个方法简单实用不会让种群发散。但要注意由于三种参数边界不同截断时不能对整个向量统一min/max而要分段处理。我在代码里会把lb和ub都定义成完整的D维向量扩展方式用repmat或者直接构造。例如中心段边界设为0和1宽度段边界设为0.1和2权值段边界设为-1和1。更新完位置后执行X_pop(i,:) max(X_pop(i,:), lb); X_pop(i,:) min(X_pop(i,:), ub);这样就不会因为宽度越界导致后面计算高斯函数时出负值。还要注意一点在更新适应度时有个容易犯的错是“这个个体更新后适应度变差却仍然保留新位置”。严格的标准WOA不比较适应度直接接受更新后的位置这有助于保持多样性。但工程实现中很多人喜欢加一个贪婪选择只保留更优的位置。两种做法都能跑但我个人在参数搜索问题上更倾向保留贪婪选择因为能明显加快收敛。如果使用贪婪选择要先把旧位置保存下来否则更新失败后个体已经被覆盖了。4. Matlab代码结构与实操细节4.1 项目文件怎么拆建议把代码拆成main.m、WOA.m、obj_fun.m、rbf_predict.m四个文件。main.m负责数据加载、归一化、参数设置、调用WOA、反归一化、结果画图WOA.m负责优化主循环obj_fun.m是适应度函数rbf_predict.m是RBF前向计算。这样拆的好处是后期换数据集、换算法、改网络层数都只需要动局部。我最早写这个项目时图省事把所有逻辑都堆在一个脚本里结果为了试一组参数就要翻半天代码。后来忍无可忍拆开调试速度一下子快了很多。特别是想对比PSO和WOA时只要把WOA.m换成PSO.mmain完全不用改。4.2 main.m主流程和归一化代码main.m大致流程如下%% 加载数据 data load(your_data.mat); X data.X; % N x m y data.y; % N x 1 %% 划分训练测试集 rng(123); idx randperm(size(X,1)); train_ratio 0.8; train_idx idx(1:floor(train_ratio*length(idx))); test_idx idx(floor(train_ratio*length(idx))1:end); X_train X(train_idx,:); y_train y(train_idx); X_test X(test_idx,:); y_test y(test_idx); %% 归一化 [x_train, PSx] mapminmax(X_train, 0, 1); x_train x_train; x_test mapminmax(apply, X_test, PSx); [y_train_norm, PSy] mapminmax(y_train, 0, 1); y_train_norm y_train_norm; y_test_norm mapminmax(apply, y_test, PSy);这里用到了mapminmax注意它默认按行处理。如果不喜欢工具箱函数也可以自己用min和max写。mapminmax的坑在于它会保存训练集的归一化参数测试集用它apply过去正好避免信息泄漏。很多初学者直接把整个X扔给mapminmax归一化后再划分这其实会把测试集信息混入训练过程严格来说是不对的。当然也可以手写归一化minx min(X_train); maxx max(X_train); x_train (X_train - minx) ./ (maxx - minx eps); x_test (X_test - minx) ./ (maxx - minx eps);这样更直观也避免了对行方向的困惑。我建议新手先手写一次搞清楚原理后再尝试mapminmax。4.3 WOA.m主函数框架WOA主函数的核心迭代代码我在前面第3节给过一部分。这里补几个容易写错的细节。第一记录最优个体时要同时保存最优位置和最优适应度。第二适应度数组在每轮要更新适应度曲线要把历史最优记录下来。第三在循环里尽量使用向量操作避免每个个体都要重算best_pos等常量。WOA.m的输入输出可以这样设计function [best_pos, best_fitness, convergence_curve] WOA(obj_fun, dim, lb, ub, pop_size, max_iter, data_struct)把数据用data_struct打包传进去obj_fun再解包可以避免主调函数传参混乱。适应度函数是函数句柄这样WOA算法和多目标函数解耦。4.4 关键参数设置一组可直接上手的默认值我给出一个稳健的默认参数组合第一次跑可以先照着设置参数默认值备注隐含节点数 H10随样本量和特征维度调整种群大小30样本多、维度高时可加到50最大迭代次数200看收敛曲线再决定是否加输入归一化区间[0,1]用训练集min/max输出归一化区间[0,1]用训练集min/max中心值边界[0,1]与输入归一化区间一致宽度边界[0.1,2]太窄会溢出太宽会过度平滑权值边界[-1,1]输出归一化后足够注意H如果太大个体维度会变成H*(m2)比如m10、H30维度就是360这会让WOA的搜索空间迅速膨胀不仅慢而且效果不一定好。先从小H开始用测试集误差判断有没有必要加节点。4.5 训练耗时与并行加速WOA耗时的瓶颈在于适应度函数要反复计算。假设pop30iter200就是6000次RBF前向计算。每次前向要建一个N×H的距离矩阵样本越多越慢。如果N5000、H20单次计算还好6000次叠加就会到几分钟甚至更久。一个简单的加速思路是用parfor替代for更新种群适应度但要注意Matlab并行池里每个worker的随机数流不同如果不处理会导致结果不可复现。另一个思路是先用少量的迭代和较小的H快速验证整个流程确定边界和预处理都正确后再放大规模跑。还有一个很实用的小技巧在obj_fun里把X和y提前归一化好并作为全局变量或struct传入不要在每轮都做归一化操作。5. 实验对比与结果分析怎么判断代码真的有效5.1 跑完代码后先看什么WOA跑完会返回一个最优个体和收敛曲线。第一件事不是看测试集误差而是看收敛曲线是否正常下降。如果曲线一开始就平着不动说明适应度函数或边界设置有问题后面调什么都是白搭。确认收敛后再计算训练集和测试集上的RMSE、MAE和R2。R2接近1不一定代表模型好还要看测试集和训练集差距大不大。如果训练集R20.98测试集R20.72基本可以判断过拟合需要减少H、增加样本或者对权值加正则约束。另外一个重要习惯是重复实验。WOA有随机性单次运行的最优结果很可能偏高或偏低。科学做法是固定随机种子做一次“正式结果”再用不同种子跑5到10次记录平均值和标准差。因为WOA不是凸优化不能保证每次都收敛到同一个解。5.2 和普通RBF、BP做对比为了说明WOA-RBF有效最好安排两组对比普通RBF和BP神经网络。普通RBF我用K-means定中心最大距离法定宽度最小二乘求权值BP用Matlab的feedforwardnet或者自己写梯度下降。在我测试的一组工艺参数回归数据上普通RBF测试RMSE约0.42BP约0.33WOA-RBF约0.23。这个结果不能代表所有数据但至少说明联合优化三组参数的价值。如果WOA-RBF反而更差大概率是边界设置不合理或H不合适。这时候不要急着否定算法先打印出最优个体的中心、宽度和权值看看有没有跑到边界上。比如宽度全部卡在下界说明最优宽度比下界还小但你又不敢设太小那就要考虑是不是数据噪声太大高斯函数过尖只会拟合噪声。5.3 后续扩展方向这套代码的扩展性不错。把适应度函数从回归误差改成分类交叉熵就可以做多输入分类把输出权值从H×1改成H×QQ是输出维度就能做多输出预测把WOA的位置更新公式换成灰狼或者粒子群就能横向对比不同优化算法。我在代码里用函数句柄传适应度函数就是为了这些扩展。6. 常见问题与避坑指南6.1 结果为NaN的三类原因如果跑出来的预测全是NaN优先检查三个地方宽度sigma是否过小或为负。高斯函数分母里有sigma平方如果sigma接近0exp里面会出现极大的负值可能溢出为NaN。下界设0.1比较稳。数据里是否有NaN或Inf。Matlab里许多函数不会主动报错而是悄悄把NaN传下去导致整个网络输出NaN。用sum(isnan(X(:)))或sum(isinf(X(:)))检查一下。归一化时是否除数为0。如果某个特征所有样本都相同max-min等于0手写归一化就会出现除零。这种情况要么删掉这个特征要么加一个eps。6.2 预测曲线像一条直线如果反归一化后的预测输出几乎平移了真实曲线而且没有局部波动很可能是宽度设置太大导致每个隐含节点对所有样本的响应都差不多网络退化成线性加权和。解决方法是把宽度上界调小比如从2改成0.5再重新跑WOA。还有一种可能是中心值没有充分分开全部挤在一个区域这时可以尝试用K-means聚类结果去初始化一部分个体帮助算法起步。6.3 调试顺序先拆齿轮再合体我强烈建议按“先简单后复杂”的顺序调试。第一步固定中心和宽度只优化权值此时RBF成了一个线性模型用最小二乘或者WOA都很容易验证。第二步固定中心优化宽度和权值。第三步三者联合优化。这个顺序能帮你快速定位问题。比如第二步效果正常但第三步反而变差那问题多半出在中心值的编码或边界上。我就踩过这样的坑中心边界设成了[-1,1]但输入归一化到了[0,1]等于让中心在无效区间里搜索适应度怎么都降不下去。后来用histogram看了最优中心的分布才发现全部挤在边界上这才想到去查归一化范围。6.4 调参心得宽度和边界是最大变量最后说说我的调参顺序。先固定H10、种群30、迭代100跑通全流程然后打印最优个体专门检查宽度有没有卡边界。如果宽度全在下界我会调低宽度下界到0.05再试如果宽度全在上界我会把上界降到1再试。中心边界一定和归一化范围对齐权值边界默认[-1,1]够用。还有一点不要一开始就追求很大的迭代次数。先用100次看收敛曲线如果曲线在第100次还在明显下降再加大到200或300。如果曲线早就平了加迭代次数也没有意义。种群大小同理20到40之间通常足够与其堆种群大小不如把边界调得更合理。最后再分享一个后期一直在用的小技巧把归一化参数、最优个体、测试集误差、随机种子一起保存成.mat文件。这样后面写报告、复现实验或者换数据集比对时随时能拿回当时的环境。毕竟WOA有随机性不固定随机种子的话过几天自己都可能复现不了。
阅读完成 · 觉得有帮助？
