GWO 灰狼优化 VMD的两个参数，惩罚系数和分解层数。并做包络谱，频谱等分析。

% 灰狼优化算法（求函数极值） clc; clear; close all; tic %% 目标函数 %% 初始化参数 N=30; % 灰狼个数 dim=2; % 维度 Iter=9; % 最大迭代次数 a=2; % 收敛因子 lb=[50 2]; % 最小值限制 ub=[5000 9]; % 最大值限制 % 初始化alpha,beta,delta Alpha_pos=zeros(1,dim); Alpha_score=inf; %求最大值改为-inf Beta_pos=zeros(1,dim); Beta_score=inf; Delta_pos=zeros(1,dim); Delta_score=inf; Positions=rand(N,dim).*(ub-lb)+lb; % 初始化个体位置 Convergence_curve=zeros(1,Iter); % 收敛曲线 l=0; %循环次数记录 %% 迭代求解 while l<Iter toc for i=1:size(Positions,1) % 超出边界处理 Flag4ub=Positions(i,:)>ub; Flag4lb=Positions(i,:)<lb; Positions(i,:)=(Positions(i,:).*(~(Flag4ub+Flag4lb)))+ub.*Flag4ub+lb.*Flag4lb; % 计算个体适应度函数!!!!重点修改适应度函数f_xy fitness=f_xy(Positions(i,1), Positions(i,2)); % 更新 Alpha, Beta, and Delta if fitness< Alpha_score Alpha_score=fitness; Alpha_pos=Positions(i,:); end if fitness>Alpha_score && fitness<Beta_score Beta_score=fitness; Beta_pos=Positions(i,:); end if fitness>Alpha_score && fitness>Beta_score && fitness<Delta_score Delta_score=fitness; Delta_pos=Positions(i,:); end end a=2-l*((2)/Iter); % 收敛因子从2线性递减到0 % 更新灰狼个体的位置 for i=1:size(Positions,1) for j=1:size(Positions,2) r1=rand(); % r1 is a random number in [0,1] r2=rand(); % r2 is a random number in [0,1] A1=2*a*r1-a; C1=2*r2; D_alpha=abs(C1*Alpha_pos(j)-Positions(i,j)); X1=Alpha_pos(j)-A1*D_alpha; r1=rand(); r2=rand(); A2=2*a*r1-a; C2=2*r2; D_beta=abs(C2*Beta_pos(j)-Positions(i,j)); X2=Beta_pos(j)-A2*D_beta; r1=rand(); r2=rand(); A3=2*a*r1-a; C3=2*r2; D_delta=abs(C3*Delta_pos(j)-Positions(i,j)); X3=Delta_pos(j)-A3*D_delta; Positions(i,j)=(X1+X2+X3)/3;% Equation (3.7) end end l=l+1; Convergence_curve(l)=Alpha_score; end figure(2); plot(Convergence_curve); title('收敛过程'); xlabel('进化代数'); ylabel('适应度值'); display(['The best solution obtained by GWO is : ', num2str(Alpha_pos)]); display(['The best optimal value of the objective funciton found by GWO is : ', num2str(Alpha_score)]); %% 最优参数的VMD %% Best_alpha =floor( Alpha_pos(1)); Best_K =floor( Alpha_pos(2)); load FZ_VMD.mat data =xx; data = data - mean(data); %采样频率 fs = 5000; %信号长度 len=5000; N =len; %% x=data(1:len); t = (0:len-1)/fs; % 采样时间 N= 5000; f = data; L = len; Fs = fs; spectrum = fft(f)'; spectrum_abs=abs(spectrum)/Fs*2; nfft=length(spectrum_abs); freqs = (0:nfft/2-1)/nfft*Fs; %--------------- Run actual VMD code if(1) %[u, u_hat, omega] = VMD(f, alpha, tau, K, DC, init, tol); tau = 0; DC = 0; init = 1; tol = 1e-6; [u, u_hat, omega] = VMD(f, Best_alpha, tau, Best_K, DC, init, tol); % f= f/0.00316;%%信号单位换算 % u= u/0.00316; figure (3) subplot(size(u,1)+1,2,1); plot(t,f,'k');grid on; xlabel('t (s)'),ylabel('原信号'), spectrum_abs_0 = zeros(7,N/2); title('VMD分解'); subplot(size(u,1)+1,2,2); spectrum_abs= abs(fft(f)); spectrum_abs_0(1,:) = spectrum_abs(1:nfft/2); plot(freqs,spectrum_abs(1:nfft/2),'k');grid on; title('对应频谱'); kur =0; kur_index =1; for i = 2:size(u,1)+1 subplot(size(u,1)+1,2,i*2-1); plot(t,u(i-1,:),'k');grid on; sprintName = [sprintf('IMF%d',i-1)]; xlabel('t (s)'),ylabel(sprintName), subplot(size(u,1)+1,2,i*2); spectrum_abs = abs( fft(u(i-1,:)) ); spectrum_abs_0(i,:) = spectrum_abs(1:nfft/2); plot(freqs,spectrum_abs(1:nfft/2),'k');grid on; kur_temp0(i-1)= skewness(u(i-1,:)); % wai度 % kur_temp(i-1) = kurtosis(u(i-1,:)); % 峭度 % % kur_baoluo(i-1) = puqiandu(u(i-1,:),5000); % 谱峭度 % b=corrcoef(u(i-1,:),f'); %相关系数 % xiangguanxishu(1,i-1)=b(1,2); %相关系数 end X = zeros(1,5000); for k = 1:Best_K % X = X + u(k,:); % E(k) = baoluoshang(u(k,:),5000); % P_kur(k) = Power_spectrum_kur(u(k,:),5000); % E_P_kur(k) = E(k) + 1./P_kur(k); E_P_kur(k)= baoluoshang(u(k,: ),5000);%+ 1./Power_spectrum_kur(u(k,: ),Fs,N); end %% [kur , E_P_index] = min(E_P_kur); end uop = u(E_P_index,:);%绘制VMD最佳分量包络谱 figure (4) y_hht=hilbert(uop);%希尔伯特变换 y_an=abs(y_hht);%包络信号 y_an_fft = abs(fft(y_an)); plot(y_an_fft(2:300)) title('VMD最佳分量包络谱'); toc