遗传算法优化VMD参数（GA-VMD）matlab代码

使用GA-VMD对一维时间序列进行分解，得到n个IMF分量，分别计算每个分量的包络熵，适应度函数为平均包络熵（mean envelope entropy)，优化VMD的两个参数（分量个数k和惩罚因子alpha）

部分matlab代码如下

主程序

%%%%%%%%%%%%%%%%%%%%GA-VMD%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%初始化%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc
clear 
fs=12000;%采样频率
Ts=1/fs;%采样周期
L=1500;%采样点数
t=(0:L-1)*Ts;%时间序列
STA=1; %采样起始位置
%----------------导入内圈故障的数据-----------------------------------------
load 105.mat
data = X105_DE_time(1:L)'; %这里可以选取DE(驱动端加速度)、FE(风扇端加速度)、BA(基座加速度)，直接更改变量名，挑选一种即可。
D=2;                                %基因数目    
NP=100;                              %染色体数目
Xs=[10,2000];                               %上限          
Xx=[3,500];                              %下限
G=1000;                               %最大遗传代数
f=zeros(D,NP);                       %初始种群赋空间
nf=zeros(D,NP);                      %子种群赋空间
Pc=0.8;                              %交叉概率
Pm=0.1;                              %变异概率
f1=rand(1,NP)*(Xs(1)-Xx(1))+Xx(1);             %随机获得初始种群
f2 = rand(1,NP)*(Xs(2)-Xx(2)) + Xx(2);
f = [f1;f2];
%%%%%%%%%%%%%%%%%%%%%%按适应度升序排列%%%%%%%%%%%%%%%%%%%%%%%
for np=1:NP
    FIT(np)=mean_envelope_entropy(data,floor(f(:,np)));
end
[SortFIT,Index]=sort(FIT);                            
Sortf=f(:,Index);
%%%%%%%%%%%%%%%%%%%%%%%遗传算法循环%%%%%%%%%%%%%%%%%%%%%%%%%%
for gen=1:G
    %%%%%%%%%%%%%%采用君主方案进行选择交叉操作%%%%%%%%%%%%%%%%
    Emper=Sortf(:,1);                      %君主染色体
    NoPoint=round(D*Pc);                   %每次交叉点的个数
    PoPoint=randi([1 D],NoPoint,NP/2);     %交叉基因的位置
    nf=Sortf;
    for i=1:NP/2
        nf(:,2*i-1)=Emper;
        nf(:,2*i)=Sortf(:,2*i);
        for k=1:NoPoint
            nf(PoPoint(k,i),2*i-1)=nf(PoPoint(k,i),2*i);
            nf(PoPoint(k,i),2*i)=Emper(PoPoint(k,i));
        end
    end
    %%%%%%%%%%%%%%%%%%%%%%%%%%变异操作%%%%%%%%%%%%%%%%%%%%%%%%%
    for m=1:NP
        for n=1:D
            r=rand(1,1);
            if r<Pm
                nf(n,m)=rand(1,1)*(Xs(n)-Xx(n))+Xx(n);
            end
        end
    end
    %%%%%%%%%%%%%%%%%%%%%子种群按适应度升序排列%%%%%%%%%%%%%%%%%%
    for np=1:NP 
          NFIT(np)=mean_envelope_entropy(data,floor(f(:,np)));   
    end
    [NSortFIT,Index]=sort(NFIT);           
    NSortf=nf(:,Index);
    %%%%%%%%%%%%%%%%%%%%%%%%%产生新种群%%%%%%%%%%%%%%%%%%%%%%%%%%
    f1=[Sortf,NSortf];                %子代和父代合并
    FIT1=[SortFIT,NSortFIT];       %子代和父代的适应度值合并
    [SortFIT1,Index]=sort(FIT1);    %适应度按升序排列
    Sortf1=f1(:,Index);               %按适应度排列个体
    SortFIT=SortFIT1(1:NP);         %取前NP个适应度值
    Sortf=Sortf1(:,1:NP);             %取前NP个个体
    trace(gen)=SortFIT(1);           %历代最优适应度值
    disp(['迭代次数:',num2str(gen)]);
end
Bestf=Sortf(:,1);                     %最优个体 
trace(end)                            %最优值
figure
plot(trace)
xlabel('迭代次数')
ylabel('目标函数值')
title('适应度进化曲线')

VMD函数如下 

function [u, u_hat, omega] = VMD(signal, alpha, tau, K, DC, init, tol)
% Variational Mode Decomposition
 
% Input and Parameters:
% ---------------------
% signal  - the time domain signal (1D) to be decomposed
% alpha   - the balancing parameter of the data-fidelity constraint 惩罚因子
% tau     - time-step of the dual ascent ( pick 0 for noise-slack )
% K       - the number of modes to be recovered 模态分量
% DC      - true if the first mode is put and kept at DC (0-freq)
% init    - 0 = all omegas start at 0
%                    1 = all omegas start uniformly distributed
%                    2 = all omegas initialized randomly
% tol     - tolerance of convergence criterion; typically around 1e-6
%
% Output:
% -------
% u       - the collection of decomposed modes
% u_hat   - spectra of the modes
% omega   - estimated mode center-frequencies
 
 
 
%---------- Preparations
% Period and sampling frequency of input signal
save_T = length(signal);
fs = 1/save_T;
 
% extend the signal by mirroring
T = save_T;
f_mirror(1:T/2) = signal(T/2:-1:1);
f_mirror(T/2+1:3*T/2) = signal;
f_mirror(3*T/2+1:2*T) = signal(T:-1:T/2+1);
f = f_mirror;
 
% Time Domain 0 to T (of mirrored signal)
T = length(f);
t = (1:T)/T;
 
% Spectral Domain discretization
freqs = t-0.5-1/T;
 
% Maximum number of iterations (if not converged yet, then it won't anyway)
N = 500;
% For future generalizations: individual alpha for each mode
Alpha = alpha*ones(1,K);
% Construct and center f_hat
f_hat = fftshift((fft(f)));
f_hat_plus = f_hat;
f_hat_plus(1:T/2) = 0;
% matrix keeping track of every iterant // could be discarded for mem
u_hat_plus = zeros(N, length(freqs), K);
% Initialization of omega_k
omega_plus = zeros(N, K);
switch init
    case 1
        for i = 1:K
            omega_plus(1,i) = (0.5/K)*(i-1);
        end
    case 2
        omega_plus(1,:) = sort(exp(log(fs) + (log(0.5)-log(fs))*rand(1,K)));
    otherwise
        omega_plus(1,:) = 0;
end
% if DC mode imposed, set its omega to 0
if DC
    omega_plus(1,1) = 0;
end
% start with empty dual variables
lambda_hat = zeros(N, length(freqs));
% other inits
uDiff = tol+eps; % update step
n = 1; % loop counter
sum_uk = 0; % accumulator
% ----------- Main loop for iterative updates
while ( uDiff > tol &&  n < N ) % not converged and below iterations limit
    
    % update first mode accumulator
    k = 1;
    sum_uk = u_hat_plus(n,:,K) + sum_uk - u_hat_plus(n,:,1);
    
    % update spectrum of first mode through Wiener filter of residuals
    u_hat_plus(n+1,:,k) = (f_hat_plus - sum_uk - lambda_hat(n,:)/2)./(1+Alpha(1,k)*(freqs - omega_plus(n,k)).^2);
    
    % update first omega if not held at 0
    if ~DC
        omega_plus(n+1,k) = (freqs(T/2+1:T)*(abs(u_hat_plus(n+1, T/2+1:T, k)).^2)')/sum(abs(u_hat_plus(n+1,T/2+1:T,k)).^2);
    end
    
    % update of any other mode
    for k=2:K
        
        % accumulator
        sum_uk = u_hat_plus(n+1,:,k-1) + sum_uk - u_hat_plus(n,:,k);
        
        % mode spectrum
        u_hat_plus(n+1,:,k) = (f_hat_plus - sum_uk - lambda_hat(n,:)/2)./(1+Alpha(1,k)*(freqs - omega_plus(n,k)).^2);
        
        % center frequencies
        omega_plus(n+1,k) = (freqs(T/2+1:T)*(abs(u_hat_plus(n+1, T/2+1:T, k)).^2)')/sum(abs(u_hat_plus(n+1,T/2+1:T,k)).^2);
        
    end
    
    % Dual ascent
    lambda_hat(n+1,:) = lambda_hat(n,:) + tau*(sum(u_hat_plus(n+1,:,:),3) - f_hat_plus);
    
    % loop counter
    n = n+1;
    
    % converged yet?
    uDiff = eps;
    for i=1:K
        uDiff = uDiff + 1/T*(u_hat_plus(n,:,i)-u_hat_plus(n-1,:,i))*conj((u_hat_plus(n,:,i)-u_hat_plus(n-1,:,i)))';
    end
    uDiff = abs(uDiff);
    
end
 
 
%------ Postprocessing and cleanup
 
 
% discard empty space if converged early
N = min(N,n);
omega = omega_plus(1:N,:);
 
% Signal reconstruction
u_hat = zeros(T, K);
u_hat((T/2+1):T,:) = squeeze(u_hat_plus(N,(T/2+1):T,:));
u_hat((T/2+1):-1:2,:) = squeeze(conj(u_hat_plus(N,(T/2+1):T,:)));
u_hat(1,:) = conj(u_hat(end,:));
 
u = zeros(K,length(t));
 
for k = 1:K
    u(k,:)=real(ifft(ifftshift(u_hat(:,k))));
end
 
% remove mirror part
u = u(:,T/4+1:3*T/4);
 
% recompute spectrum
clear u_hat;
for k = 1:K
    u_hat(:,k)=fftshift(fft(u(k,:)))';
end
end

适应度函数（平均包络熵）如下 

 

function fitness = mean_envelope_entropy(data,x)
%--------- some sample parameters forVMD：对于VMD样品参数进行设置---------------
alpha = x(2);       % moderate bandwidth constraint：适度的带宽约束/惩罚因子
tau = 0;          % noise-tolerance (no strict fidelity enforcement)：噪声容限（没有严格的保真度执行）
K = x(1);              % modes：分解的模态数，可以自行设置，这里以8为例。
DC = 0;             % no DC part imposed：无直流部分
init = 1;           % initialize omegas uniformly  ：omegas的均匀初始化
tol = 1e-7;        
%--------------- Run actual VMD code:数据进行vmd分解---------------------------
[u, u_hat, omega] = VMD(data, alpha, tau, K, DC, init, tol); %其中u为分解得到的IMF分量

% 假设 u 是VMD分解得到的IMFs集合  
% u 的每一行是一个IMF  
  
% 初始化包络熵数组  
envelopeEntropy = zeros(size(u, 1), 1);  
  
% 对每个IMF计算包络和熵  
for i = 1:size(u, 1)  
    % 计算IMF的包络  
    % 使用Hilbert变换获取解析信号，然后取其绝对值作为包络  
    analyticSignal = hilbert(u(i, :));  
    envelope = abs(analyticSignal);  
      
    % 归一化包络以便计算熵  
    envelope = envelope / sum(envelope);  
      
    % 计算包络的熵  
    % 香农熵计算公式  
    envelopeEntropy(i) = -sum(envelope .* log2(envelope + eps));  
end  
  fitness = mean(envelopeEntropy);
% 显示每个IMF的包络熵  
disp('Envelope Entropy for each IMF:');  
disp(envelopeEntropy);
disp(['平均包络熵:',num2str(fitness)]);
end

VMD代码来自如下网站

https://blog.csdn.net/woaipythonmeme/article/details/128702229?ops_request_misc=%257B%2522request%255Fid%2522%253A%2522171214981916800213025060%2522%252C%2522scm%2522%253A%252220140713.130102334..%2522%257D&request_id=171214981916800213025060&biz_id=0&utm_medium=distribute.pc_search_result.none-task-blog-2~all~top_positive~default-1-128702229-null-null.142^v100^pc_search_result_base1&utm_term=VMD%20matlab&spm=1018.2226.3001.4187https://blog.csdn.net/woaipythonmeme/article/details/128702229?ops_request_misc=%257B%2522request%255Fid%2522%253A%2522171214981916800213025060%2522%252C%2522scm%2522%253A%252220140713.130102334..%2522%257D&request_id=171214981916800213025060&biz_id=0&utm_medium=distribute.pc_search_result.none-task-blog-2~all~top_positive~default-1-128702229-null-null.142^v100^pc_search_result_base1&utm_term=VMD%20matlab&spm=1018.2226.3001.4187