代碼
先對外層代碼的揭露，包括：順序而下
1、
function s = Bounds( s, Lb, Ub)
% Apply the lower bound vector
temp = s;
I = temp < Lb;
temp(I) = Lb(I);
% Apply the upper bound vector
J = temp > Ub;
temp(J) = Ub(J);
% Update this new move
s = temp;
2、
function [in,out]=data_process(data,num)
% 采用1-num的各種值為輸入 第num+1的功率作為輸出
n=length(data)-num;
for i=1:n
x(i,:,:)=data(i:i+num,:);
end
in=x(1:end-1,:);
out=x(2:end,end);%功率是最后一列
3、
function y=fitness(x,train_x,train_y,test_x,test_y,typeID,kernelnamescell)
gam=x(1);
sig2=x(2);
model=initlssvm(train_x,train_y,typeID,gam,sig2,kernelnamescell);
model=trainlssvm(model);
y_pre_test=simlssvm(model,test_x);
y=mse(test_y-y_pre_test);
4、
clc;clear;close all
%%
lssvm=load('lssvm.mat');
ssa_lssvm=load('ssa_lssvm.mat');
vmd_lssvm=load('vmd_lssvm.mat');
vmd_ssa_lssvm=load('vmd_ssa_lssvm.mat');
disp('結(jié)果分析-lssvm')
result(lssvm.true_value,lssvm.predict_value)
fprintf('\n')
disp('結(jié)果分析-ssa-lssvm')
result(ssa_lssvm.true_value,ssa_lssvm.predict_value)
fprintf('\n')
disp('結(jié)果分析-vmd-lssvm')
result(vmd_lssvm.true_value,vmd_lssvm.predict_value)
fprintf('\n')
disp('結(jié)果分析-vmd-ssa-lssvm')
result(vmd_ssa_lssvm.true_value,vmd_ssa_lssvm.predict_value)
figure
plot(lssvm.true_value)
hold on;grid on
plot(lssvm.predict_value)
plot(ssa_lssvm.predict_value)
plot(vmd_lssvm.predict_value)
plot(vmd_ssa_lssvm.predict_value)
legend('真實(shí)值','lssvm預(yù)測','ssa-lssvm預(yù)測','vmd-lssvm預(yù)測','vmd-ssa-lssvm預(yù)測')
title('各算法預(yù)測效果對比')
ylabel('功率值')
xlabel('測試集樣本')
5、
%% 其他數(shù)據(jù)與功率數(shù)據(jù)一起作為輸入，參看data_process.m
clc;clear;close all
addpath(genpath('LSSVMlabv1_8'));
%%
data=xlsread('預(yù)測數(shù)據(jù).xls','B2:K1000');
[x,y]=data_process(data,24);%前24個時刻 預(yù)測下一個時刻
%歸一化
[xs,mappingx]=mapminmax(x',0,1);x=xs';
[ys,mappingy]=mapminmax(y',0,1);y=ys';
%劃分?jǐn)?shù)據(jù)
n=size(x,1);
m=round(n*0.7);%前70%訓(xùn)練，對最后30%進(jìn)行預(yù)測
train_x=x(1:m,:);
test_x=x(m+1:end,:);
train_y=y(1:m,:);
test_y=y(m+1:end,:);
%% lssvm 參數(shù)
gam=10;
sig2=1000;
typeID='function estimation';
kernelnamescell='RBF_kernel';
model=initlssvm(train_x,train_y,typeID,gam,sig2,kernelnamescell);
model=trainlssvm(model);
y_pre_test=simlssvm(model,test_x);
% 反歸一化
predict_value=mapminmax('reverse',y_pre_test',mappingy);
true_value=mapminmax('reverse',test_y',mappingy);
save lssvm predict_value true_value
disp('結(jié)果分析')
rmse=sqrt(mean((true_value-predict_value).^2));
disp(['根均方差(RMSE)：',num2str(rmse)])
mae=mean(abs(true_value-predict_value));
disp(['平均絕對誤差（MAE）：',num2str(mae)])
mape=mean(abs(true_value-predict_value)/true_value);
disp(['平均相對百分誤差（MAPE）：',num2str(mape*100),'%'])
fprintf('\n')
figure
plot(true_value,'-*','linewidth',3)
hold on
plot(predict_value,'-s','linewidth',3)
legend('實(shí)際值','預(yù)測值')
grid on
title('LSSVM')
6、
clc;clear;close all;format compact
addpath(genpath('LSSVMlabv1_8'));
%%
data=xlsread('預(yù)測數(shù)據(jù).xls','B2:K1000');
[x,y]=data_process(data,24);%前24個時刻 預(yù)測下一個時刻
%歸一化
[xs,mappingx]=mapminmax(x',0,1);x=xs';
[ys,mappingy]=mapminmax(y',0,1);y=ys';
%劃分?jǐn)?shù)據(jù)
n=size(x,1);
m=round(n*0.7);%前70%訓(xùn)練，對最后30%進(jìn)行預(yù)測
train_x=x(1:m,:);
test_x=x(m+1:end,:);
train_y=y(1:m,:);
test_y=y(m+1:end,:);
%% ssa優(yōu)化lssvm 參數(shù)
typeID='function estimation';
kernelnamescell='RBF_kernel';
[x,trace]=ssa_lssvm(typeID,kernelnamescell,train_x,train_y,test_x,test_y);
figure;plot(trace);title('適應(yīng)度曲線/mse')
%% 利用尋優(yōu)得到的參數(shù)重新訓(xùn)練lssvm
disp('尋優(yōu)得到的參數(shù)分別是：')
gam=x(1)
sig2=x(2)
model=initlssvm(train_x,train_y,typeID,gam,sig2,kernelnamescell);
model=trainlssvm(model);
y_pre_test=simlssvm(model,test_x);
% 反歸一化
predict_value=mapminmax('reverse',y_pre_test',mappingy);
true_value=mapminmax('reverse',test_y',mappingy);
save ssa_lssvm predict_value true_value
disp('結(jié)果分析')
rmse=sqrt(mean((true_value-predict_value).^2));
disp(['根均方差(RMSE)：',num2str(rmse)])
mae=mean(abs(true_value-predict_value));
disp(['平均絕對誤差（MAE）：',num2str(mae)])
mape=mean(abs(true_value-predict_value)/true_value);
disp(['平均相對百分誤差（MAPE）：',num2str(mape*100),'%'])
fprintf('\n')
figure
plot(true_value,'-*','linewidth',3)
hold on
plot(predict_value,'-s','linewidth',3)
legend('實(shí)際值','預(yù)測值')
grid on
title('SSA-LSSVM')
7、
clc;clear;format compact;close all;tic
rng('default')
%% 數(shù)據(jù)預(yù)處理
data=xlsread('預(yù)測數(shù)據(jù).xls','B2:K1000');
figure;plot(data(:,end),'-*','linewidth',3);title('原始功率曲線');grid on;ylabel('功率值')
alpha = 2000; % moderate bandwidth constraint
tau = 0; % noise-tolerance (no strict fidelity enforcement)
K = 5; % 3 modes
DC = 0; % no DC part imposed
init = 1; % initialize omegas uniformly
tol = 1e-7;
[imf, u_hat, omega] = VMD(data(:,end), alpha, tau, K, DC, init, tol);
figure
for i =1:size(imf,1)
subplot(size(imf,1),1,i)
plot(imf(i,:))
ylabel(['imf',num2str(i)])
end
ylabel('殘余')
suptitle('VMD')
c=size(imf,1);
pre_result=[];
true_result=[];
%% 對每個分量建模
for i=1:c
disp(['對第',num2str(i),'個分量建模'])
[x,y]=data_process([data(:,1:end-1) imf(i,:)'],24);%每個序列和原始的幾列數(shù)據(jù)合并 然后劃分
%歸一化
[xs,mappingx]=mapminmax(x',0,1);x=xs';
[ys,mappingy]=mapminmax(y',0,1);y=ys';
%劃分?jǐn)?shù)據(jù)
n=size(x,1);
m=round(n*0.7);%前70%訓(xùn)練，對最后30%進(jìn)行預(yù)測
train_x=x(1:m,:);
test_x=x(m+1:end,:);
train_y=y(1:m,:);
test_y=y(m+1:end,:);
%%
gam=100;
sig2=1000;
typeID='function estimation';
kernelnamescell='RBF_kernel';
model=initlssvm(train_x,train_y,typeID,gam,sig2,kernelnamescell);
model=trainlssvm(model);
pred=simlssvm(model,test_x);
pre_value=mapminmax('reverse',pred',mappingy);
true_value=mapminmax('reverse',test_y',mappingy);
pre_result=[pre_result;pre_value];
true_result=[true_result;true_value];
end
%% 各分量預(yù)測的結(jié)果相加
true_value=sum(true_result);
predict_value=sum(pre_result);
save vmd_lssvm predict_value true_value
%%
load vmd_lssvm
disp('結(jié)果分析')
rmse=sqrt(mean((true_value-predict_value).^2));
disp(['根均方差(RMSE)：',num2str(rmse)])
mae=mean(abs(true_value-predict_value));
disp(['平均絕對誤差（MAE）：',num2str(mae)])
mape=mean(abs(true_value-predict_value)/true_value);
disp(['平均相對百分誤差（MAPE）：',num2str(mape*100),'%'])
fprintf('\n')
figure
plot(true_value,'-*','linewidth',3)
hold on
plot(predict_value,'-s','linewidth',3)
legend('實(shí)際值','預(yù)測值')
grid on
title('VMD-LSSVM')
8、
clc;clear;format compact;close all;tic
rng('default')
%% 數(shù)據(jù)預(yù)處理
data=xlsread('預(yù)測數(shù)據(jù).xls','B2:K1000');
figure;plot(data(:,end),'-*','linewidth',3);title('原始功率曲線');grid on;ylabel('功率值')
alpha = 2000; % moderate bandwidth constraint
tau = 0; % noise-tolerance (no strict fidelity enforcement)
K = 5; % 3 modes
DC = 0; % no DC part imposed
init = 1; % initialize omegas uniformly
tol = 1e-7;
[imf, u_hat, omega] = VMD(data(:,end), alpha, tau, K, DC, init, tol);
figure
for i =1:size(imf,1)
subplot(size(imf,1),1,i)
plot(imf(i,:))
ylabel(['imf',num2str(i)])
end
ylabel('殘余')
suptitle('VMD')
c=size(imf,1);
pre_result=[];
true_result=[];
%% 對每個分量建模
for i=1:c
disp(['對第',num2str(i),'個分量建模'])
[x,y]=data_process([data(:,1:end-1) imf(i,:)'],24);
%歸一化
[xs,mappingx]=mapminmax(x',0,1);x=xs';
[ys,mappingy]=mapminmax(y',0,1);y=ys';
%劃分?jǐn)?shù)據(jù)
n=size(x,1);
m=round(n*0.7);%前70%訓(xùn)練，對最后30%進(jìn)行預(yù)測
train_x=x(1:m,:);
test_x=x(m+1:end,:);
train_y=y(1:m,:);
test_y=y(m+1:end,:);
%%
typeID='function estimation';
kernelnamescell='RBF_kernel';
[x,trace]=ssa_lssvm(typeID,kernelnamescell,train_x,train_y,test_x,test_y);
gam=x(1);
sig2=x(2);
model=initlssvm(train_x,train_y,typeID,gam,sig2,kernelnamescell);
model=trainlssvm(model);
pred=simlssvm(model,test_x);
pre_value=mapminmax('reverse',pred',mappingy);
true_value=mapminmax('reverse',test_y',mappingy);
pre_result=[pre_result;pre_value];
true_result=[true_result;true_value];
end
%% 各分量預(yù)測的結(jié)果相加
true_value=sum(true_result);
predict_value=sum(pre_result);
save vmd_ssa_lssvm predict_value true_value
%%
load vmd_ssa_lssvm
disp('結(jié)果分析')
rmse=sqrt(mean((true_value-predict_value).^2));
disp(['根均方差(RMSE)：',num2str(rmse)])
mae=mean(abs(true_value-predict_value));
disp(['平均絕對誤差（MAE）：',num2str(mae)])
mape=mean(abs(true_value-predict_value)/true_value);
disp(['平均相對百分誤差（MAPE）：',num2str(mape*100),'%'])
fprintf('\n')
figure
plot(true_value,'-*','linewidth',3)
hold on
plot(predict_value,'-s','linewidth',3)
legend('實(shí)際值','預(yù)測值')
grid on
title('VMD-SSA-LSSVM')
9、
clear;
close all;
clc;
format compact
addpath('vmd-verify')
data=xlsread('預(yù)測數(shù)據(jù).xls','B2:K1000');
f=data(:,end);
% some sample parameters for VMD
alpha = 2000; % moderate bandwidth constraint
tau = 0; % noise-tolerance (no strict fidelity enforcement)
K = 5; % 3 modes
DC = 0; % no DC part imposed
init = 1; % initialize omegas uniformly
tol = 1e-7;
%--------------- Run actual VMD code
[u, u_hat, omega] = VMD(f, alpha, tau, K, DC, init, tol);
figure
plot(f)
title('原始')
figure
for i=1:K
subplot(K,1,i)
plot(u(i,:))
ylabel(['IMF_',num2str(i)])
end
ylabel('res')
10、
function result(true_value,predict_value)
rmse=sqrt(mean((true_value-predict_value).^2));
disp(['根均方差(RMSE)：',num2str(rmse)])
mae=mean(abs(true_value-predict_value));
disp(['平均絕對誤差（MAE）：',num2str(mae)])
mape=mean(abs(true_value-predict_value)/true_value);
disp(['平均相對百分誤差（MAPE）：',num2str(mape*100),'%'])
11、
function [bestX,Convergence_curve]=ssa_lssvm(typeID,Kernel_type,inputn_train,label_train,inputn_test,label_test)
%% 麻雀優(yōu)化
pop=10; % 麻雀數(shù)
M=10; % Maximum numbef of iterations
c=1;
d=10000;
dim=2;
P_percent = 0.2; % The population size of producers accounts for "P_percent" percent of the total population size
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
pNum = round( pop * P_percent ); % The population size of the producers
lb= c.*ones( 1,dim ); % Lower limit/bounds/ a vector
ub= d.*ones( 1,dim ); % Upper limit/bounds/ a vector
%Initialization
for i = 1 : pop
x( i, : ) = lb + (ub - lb) .* rand( 1, dim );
fit( i )=fitness(x(i,:),inputn_train,label_train,inputn_test,label_test,typeID,Kernel_type);
end
pFit = fit;
pX = x; % The individual's best position corresponding to the pFit
[ fMin, bestI ] = min( fit ); % fMin denotes the global optimum fitness value
bestX = x( bestI, : ); % bestX denotes the global optimum position corresponding to fMin
for t = 1 : M
[ ans, sortIndex ] = sort( pFit );% Sort.
[fmax,B]=max( pFit );
worse= x(B,:);
r2=rand(1);
%%%%%%%%%%%%%5%%%%%%這一部位為發(fā)現(xiàn)者（探索者）的位置更新%%%%%%%%%%%%%%%%%%%%%%%%%
if(r2<0.8)%預(yù)警值較小，說明沒有捕食者出現(xiàn)
for i = 1 : pNum %r2小于0.8的發(fā)現(xiàn)者的改變（1-20） % Equation (3)
r1=rand(1);
x( sortIndex( i ), : ) = pX( sortIndex( i ), : )*exp(-(i)/(r1*M));%對自變量做一個隨機(jī)變換
x( sortIndex( i ), : ) = Bounds( x( sortIndex( i ), : ), lb, ub );%對超過邊界的變量進(jìn)行去除
fit( sortIndex( i ) )=fitness(x(sortIndex( i ),:),inputn_train,label_train,inputn_test,label_test,typeID,Kernel_type);
end
else %預(yù)警值較大，說明有捕食者出現(xiàn)威脅到了種群的安全，需要去其它地方覓食
for i = 1 : pNum %r2大于0.8的發(fā)現(xiàn)者的改變
x( sortIndex( i ), : ) = pX( sortIndex( i ), : )+randn(1)*ones(1,dim);
x( sortIndex( i ), : ) = Bounds( x( sortIndex( i ), : ), lb, ub );
fit( sortIndex( i ) )=fitness(x(sortIndex( i ),:),inputn_train,label_train,inputn_test,label_test,typeID,Kernel_type);
end
end
[ fMMin, bestII ] = min( fit );
bestXX = x( bestII, : );
%%%%%%%%%%%%%5%%%%%%這一部位為加入者（追隨者）的位置更新%%%%%%%%%%%%%%%%%%%%%%%%%
for i = ( pNum + 1 ) : pop %剩下20-100的個體的變換 % Equation (4)
A=floor(rand(1,dim)*2)*2-1;
if( i>(pop/2))%這個代表這部分麻雀處于十分饑餓的狀態(tài)（因?yàn)樗鼈兊哪芰亢艿?，也是是適應(yīng)度值很差），需要到其它地方覓食
x( sortIndex(i ), : )=randn(1)*exp((worse-pX( sortIndex( i ), : ))/(i)^2);
else%這一部分追隨者是圍繞最好的發(fā)現(xiàn)者周圍進(jìn)行覓食，其間也有可能發(fā)生食物的爭奪，使其自己變成生產(chǎn)者
x( sortIndex( i ), : )=bestXX+(abs(( pX( sortIndex( i ), : )-bestXX)))*(A'*(A*A')^(-1))*ones(1,dim);
end
x( sortIndex( i ), : ) = Bounds( x( sortIndex( i ), : ), lb, ub );%判斷邊界是否超出
fit( sortIndex( i ) )=fitness(x(sortIndex( i ),:),inputn_train,label_train,inputn_test,label_test,typeID,Kernel_type);
end
%%%%%%%%%%%%%5%%%%%%這一部位為意識到危險（注意這里只是意識到了危險，不代表出現(xiàn)了真正的捕食者）的麻雀的位置更新%%%%%%%%%%%%%%%%%%%%%%%%%
c=randperm(numel(sortIndex));%%%%%%%%%這個的作用是在種群中隨機(jī)產(chǎn)生其位置（也就是這部分的麻雀位置一開始是隨機(jī)的，意識到危險了要進(jìn)行位置移動，
%處于種群外圍的麻雀向安全區(qū)域靠攏，處在種群中心的麻雀則隨機(jī)行走以靠近別的麻雀）
b=sortIndex(c(1:10));
for j = 1 : length(b) % Equation (5)
if( pFit( sortIndex( b(j) ) )>(fMin) ) %處于種群外圍的麻雀的位置改變
x( sortIndex( b(j) ), : )=bestX+(randn(1,dim)).*(abs(( pX( sortIndex( b(j) ), : ) -bestX)));
else %處于種群中心的麻雀的位置改變
x( sortIndex( b(j) ), : ) =pX( sortIndex( b(j) ), : )+(2*rand(1)-1)*(abs(pX( sortIndex( b(j) ), : )-worse))/ ( pFit( sortIndex( b(j) ) )-fmax+1e-50);
end
x( sortIndex(b(j) ), : ) = Bounds( x( sortIndex(b(j) ), : ), lb, ub );
fit( sortIndex( b(j) ) )=fitness(x(sortIndex(b(j) ),:),inputn_train,label_train,inputn_test,label_test,typeID,Kernel_type);
end
for i = 1 : pop
if ( fit( i ) < pFit( i ) )
pFit( i ) = fit( i );
pX( i, : ) = x( i, : );
end
if( pFit( i ) < fMin )
fMin= pFit( i );
bestX = pX( i, : );
end
end
Convergence_curve(t,:)=[fMin mean(pFit)];
end
接下來是內(nèi)嵌代碼，就是下面兩個文件夾的代碼，實(shí)在是太多，想要的留言吧！
數(shù)據(jù)
數(shù)據(jù)是由時間、風(fēng)速、風(fēng)向等等因素組成的文件。
結(jié)果
文章來源：http://www.zghlxwxcb.cn/news/detail-517827.html
結(jié)果圖太多，就先給出這么多，如有需要代碼和數(shù)據(jù)的同學(xué)請在評論區(qū)發(fā)郵箱，一般一天之內(nèi)會回復(fù)，請點(diǎn)贊＋關(guān)注謝謝??！文章來源地址http://www.zghlxwxcb.cn/news/detail-517827.html
到了這里，關(guān)于機(jī)器學(xué)習(xí)之MATLAB代碼--基于VMD與SSA優(yōu)化lssvm的功率預(yù)測(多變量）（七）的文章就介紹完了。如果您還想了解更多內(nèi)容，請在右上角搜索TOY模板網(wǎng)以前的文章或繼續(xù)瀏覽下面的相關(guān)文章，希望大家以后多多支持TOY模板網(wǎng)！
