【优化算法】VMD分解算法的16种优化，对K和alpha参数寻优，附MATLAB代码

在上一篇文章中，我们介绍了优化算法的基本原理和一些常见的生物启发式算法。另外我们封装了一个16合1的寻优函数。

不过在上一篇中，我们举了一个简单的数值模型作为适应度函数的演示案例，然而在实际的研究中，适应度函数往往要复杂得多。本篇我们就以VMD算法的优化为例，讲一讲这种较为复杂的算法的寻优该怎么做，此外我还提供了修改好的VMD寻优代码，供大家参考使用！

一、为什么要对VMD做参数寻优

在一众的类EMD算法中，对于EMD、EEMD、CEEMD等等这些无法指定分解模态数的算法，经常有同学问我该怎样指定模态数量；当我告诉大家需要使用VMD才能制定模态数量K时，更多的是问我K值该怎么取。人类就是这么的纠结。

今天这篇文章就可以解答这个疑问了。话说回来，K的选取确实比较重要，K 的取值直接决定了分解后的模态分量的数目，若 K 设置得过大，会造成欠分解产生虚假模态；若 K 设置得偏小，会造成欠分解不能充分提取时间序列隐含特征。

此外惩罚因子alpha对分解结果也有重要影响，alpha的取值影响了模态分量的带宽，与带宽成反比，若alpha设置较小，则会出现模态混叠影响特征提取；若alpha设置较大，虽然能避免模态混叠但又造成了局部信息丢失的问题。

在有些研究中，对K和alpha的取值问题，会采用多次实验选取预测/分类效果最优值，但是这十分依赖研究人员的主观判断，缺乏客观的评判标准。因此,对VMD的参数进行优化是十分必要的。

二、VMD寻优的参数和适应度函数

1.待优化参数

就像上边所说，分解模态数K和惩罚因子alpha是比较关键的两个参数。其实还有第三个常用参数tol，不过因为tol对分解结果的影响相对较小，所以在本文中，不将tol作为寻优对象参数。

2.适应度函数

对VMD参数选择优化的关键是适应度函数的设置，用来对分解效果制定量化指标来衡量参数选择的优劣。

纵观众多论文，VMD基本都是使用的信息熵作为适应度函数[1]。

在之前的文章中讲过，信息熵越大，代表不确定性越大，信号中包含的信息量越少。

所以为了最大程度地提取出有效信息，就需要让分解出的各个模态的平均信息熵最小。

至于“信息熵”的具体类型选择就比较多了，我们之前讲过的功率谱熵、奇异谱熵、能量熵、近似熵、样本熵、模糊熵、排列熵，都可以选择。

三、应用kOptimizationAlgorithm函数实现16种优化算法寻优

有了上述定义,我们就可以利用kOptimizationAlgorithm函数来实现对VMD参数的自动寻优了。具体步骤如下：

3.1 定义适应度函数

下边是改好的适应度函数，其中tol值指定为1*10-6。适应度函数主要就三步：

第一步进行VMD分解，这里调用了之前封装过的kVMD函数，当然这里直接用MATLAB的vmd函数也可以，但是需要注意转置，因为MATLAB官方vmd函数得到的imf是按列方向排列的。

第二步是提取各个imf分量的排列熵特征，这里也使用了之前封装的熵特征提取函数，只需要指定一下熵类型名称以及必要的参数设置就行，大家可以更换这一步的熵类型，但是需要注意option也需要对应调整，具体设置方式可以参考这里：（功率谱熵、奇异谱熵、能量熵、近似熵、样本熵、模糊熵、排列熵）

第三步是求各个imf分量的排列熵的均值，并将其作为适应度函数值。

%% 适应度函数，此代码需要根据应用场景做适应性修改
function [fitness] = funFitness(x,dataforFitness)
% 输入：
% x是待优化参数，为数组。x(1)是第一个待优化参数，x(2)是第二个，以此类推。
% dataforFitness为结构体，表示导入到适应度函数的参数，按照数据实际情况设置。可以没有
% 输出：
% fitness是适应度函数
%%
% 待寻优参数读取，注意此处可能有些参数必须为整型，需要取整
alpha = round(x(1)); %alpha 
K = round(x(2)); %K
% 导入数据读取
FsOrT = dataforFitness.fs;
data  = dataforFitness.data;
try
    % 定义目标函数
    tol = 1e-6;
    [imf,CenFs]  = kVMD(data,FsOrT, alpha, K, tol);


    % 设置排列熵参数，Pedim为排列熵模式维度，Pet为排列熵的时间延迟，如果不提取排列熵特征可以删除以下两行
    option.Pedim = 6;
    option.Pet   = 1;
    featureNamesCell = {'PeEn'}; %要进行特征提取的特征名称，可以根据实际需要进行删减，留下的特征注意拼写正确
    fea = genFeatureEn(imf,featureNamesCell,option);  %调用genFeature函数，完成特征提取，算出的特征值会保存在fea变量里
                                             
    fitness = mean(fea);

    
catch ME 
    fprintf('发生错误: %s\n', ME.message);
end
end

3.2 设置寻优相关参数

以下代码是写在主程序里的了，必要的参数我都将其作为可调参数抽取出来了，大家可以根据需要进行调整，每个参数的含义在注释里标明了，这里不再赘述。

%% 1.设置寻优相关参数
dim = 2;             % 优化参数的个数，这里设定为2，表示有两个需要寻优的参数
lb = [100,2];        % 定义每个优化参数的取值下限，这里的写法代表：alpha下限设置为100，K下限设置为2
ub = [10000,10];     % 定义每个优化参数的取值上限，这里的写法代表：alpha上限设置为10000，K上限设置为10
OAmethod = 'SSA';    % 选择使用的优化算法，这里使用'SSA'即麻雀搜索算法
pop = 3;             % 设置优化算法的种群大小
Max_iteration = 40;  % 设置优化算法的最大迭代次数

3.3 调用kOptimizationAlgorithm函数实现寻优

以下代码也是写在主程序里，dataforFitness作为结构体用于传递必要的数值，这样做大大提升了代码的灵活性和简洁性，采样频率要根据你的信号情况设置（长周期信号设置为1）；dataforFitness.data则是要将你待分解信号的变量名赋值给他。最后一行直接调用即可实现寻优了。

%% 2.调用优化算法函数，传入上述定义的参数
dataforFitness.fs = 1000;  %设置采样频率
dataforFitness.data = x;   %设置待分解信号
% 返回最优位置Best_pos，最优适应度值Best_score，迭代过程的适应度曲线curve
[Best_pos,Best_score,curve] = kOptimizationAlgorithm(OAmethod, dataforFitness, lb, ub, dim, pop, Max_iteration);