tic % 计时
%% 清空环境导入数据
clear
clc
close all
load wndspd % 示例数据为风速(时间序列)数据,共144个样本
%% PSO-SVR
% 训练/测试数据准备(用前3天预测后一天),用前100天做测试数据
train_input(1,:)=wndspd(1:97);
train_input(2,:)=wndspd(2:98);
train_input(3,:)=wndspd(3:99);
train_output=[wndspd(4:100)]';
test_input(1,:)=wndspd(101:end-3);
test_input(2,:)=wndspd(102:end-2);
test_input(3,:)=wndspd(103:end-1);
test_output=[wndspd(104:end)]';
[input_train,rule1]=mapminmax(train_input);
[output_train,rule2]=mapminmax(train_output);
input_test=mapminmax('apply',test_input,rule1);
output_test=mapminmax('apply',test_output,rule2);
%% 参数设置
% 粒子位置和速度更新参数
c1=1.5;
c2=1.5;
maxgen=50; % 迭代次数
sizepop=20; % 种群规模(粒子数)
% 目标函数信息
dim=2; % 待优化参数个数
xub=[100,100]; % 参数取值上界
xlb=[0.01,0.01]; % 参数取值下界
vub=[5,5]; % 速度上界
vlb=[-5,-5]; % 速度下界
%% 初始化
% 粒子位置初始化
xRange=repmat((xub-xlb),[sizepop,1]);
xLower=repmat(xlb,[sizepop,1]);
pop=rand(sizepop,dim).*xRange+xLower;
% 粒子速度初始化
vRange=repmat((vub-vlb),[sizepop,1]);
vLower=repmat(vlb,[sizepop,1]);
V=rand(sizepop,dim).*vRange+vLower;
% 计算初始化的目标函数值(适应度函数值)
fitness=ones(sizepop,1);
for k=1:sizepop
fitness(k)=fobj(pop(k,:),input_train,output_train,input_test,output_test);
end
%% 寻找初始极值
[bestfitness,bestindex]=min(fitness);
zbest=pop(bestindex,:); % 全局最优解的位置
gbest=pop; % 个体位置
fitnessgbest=fitness; % 个体适应度值
fitnesszbest=bestfitness; % 全局最优适应度值
%% 迭代寻优
yy=ones(sizepop,1); % 用于保存每次迭代的最优目标函数值
for k=1:maxgen
% 粒子位置和速度更新
for m=1:sizepop
% 粒子速度更新
sol=V(m,:)+c1*rand*(gbest(m,:)-pop(m,:))+c2*rand*(zbest-pop(m,:));
% 确保粒子速度取值范围不越界
ind=find(sol<vlb);
sol(ind)=vlb(ind);
ind=find(sol>vub);
sol(ind)=vub(ind);
V(m,:)=sol;
% 粒子位置更新
sol=pop(m,:)+0.5*V(m,:);
% 确保粒子位置取值范围不越界
ind=find(sol<xlb);
sol(ind)=xlb(ind);
ind=find(sol>xub);
sol(ind)=xub(ind);
pop(m,:)=sol;
% 更新粒子适应度值
fitness(m)=fobj(pop(m,:),input_train,output_train,input_test,output_test);
end
% 个体极值及位置和群体极值及位置更新
for m=1:sizepop
% 个体极值及其位置更新
if fitness(m)<fitnessgbest(m)
gbest(m,:)=pop(m,:);
fitnessgbest(m)=fitness(m);
end
% 群体极值及其位置更新
if fitness(m)<fitnesszbest
zbest=pop(m,:);
fitnesszbest=fitness(m);
end
end
% 记录每代最优值
yy(k)=fitnesszbest;
end
%% 打印参数选择结果
bestobjfun=fitnesszbest;
bestc=zbest(1);
bestg=zbest(2);
disp('打印参数选择结果');
str=sprintf('Best c = %g,Best g = %g',bestc,bestg);
disp(str)
%% 利用回归预测分析最佳的参数进行SVM网络训练
cmd_cs_svr=['-s 3 -t 2',' -c ',num2str(bestc),' -g ',num2str(bestg)];
model_cs_svr=svmtrain(output_train',input_train',cmd_cs_svr); % SVM模型训练
%% SVM网络回归预测
[output_test_pre,acc,~]=svmpredict(output_test',input_test',model_cs_svr); % SVM模型预测及其精度
test_pre=mapminmax('reverse',output_test_pre',rule2);
test_pre = test_pre';
err_pre=wndspd(104:end)-test_pre;
figure('Name','测试数据残差图')
set(gcf,'unit','centimeters','position',[0.5,5,30,5])
plot(err_pre,'*-');
figure('Name','原始-预测图')
plot(test_pre,'*r-');hold on;plot(wndspd(104:end),'bo-');
legend('预测','原始')
set(gcf,'unit','centimeters','position',[0.5,13,30,5])
result=[wndspd(104:end),test_pre];
error=wndspd(104:end)-test_pre;
MAE=norm(error)
MSE=mse(error)
MAPE=mean(abs(error)./wndspd(104:end))
%% 显示程序运行时间
toc
评论2