svm.rar_ABC_ABC-SVM_ABC_SVM_支持向量机_改进ABC

tic % 计时 %% 清空环境，准备数据 clear clc close all load wndspd % 示例数据为风速（时间序列）数据，共144个样本 % 训练/测试数据准备（用前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); %% %%%%%%%%%%%%%用ABC算法优化SVR中的参数c和g开始%%%%%%%%%%%%%%%%%%%% %% 参数初始化 NP=20; % 蜂群规模 FoodNumber=NP/2; % 蜜源（解）数量 limit=100; % 当有蜜源连续没被更新的次数超过limit时，该蜜源将被重新初始化 maxCycle=10; % 最大迭代次数 % 待优化参数信息 D=2; % 待优化参数个数，次数为c和g两个 ub=ones(1,D)*100; % 参数取值上界，此处将c和g的上界设为100 lb=ones(1,D)*(0.01); % 参数取值下界，此处将c和g的下界设为0.01 runtime=2; % 可用于设置多次运行（让ABC算法运行runtime次）以考察程序的稳健性 BestGlobalMins=ones(1,runtime); % 全局最小值初始化，这里的优化目标为SVR预测结果中的平均平方误差（MSE），初始化为最差值1 BestGlobalParams=zeros(runtime,D); % 用于存放ABC算法优化得到的最优参数 for r=1:runtime % 运行ABC算法runtime次 % 初始化蜜源 Range = repmat((ub-lb),[FoodNumber 1]); Lower = repmat(lb, [FoodNumber 1]); Foods = rand(FoodNumber,D) .* Range + Lower; % 计算每个蜜源（解）得目标函数值，fobj为计算SVR预测的平均平方误差（MSE）的函数，根据自己的实际问题变异目标函数即可 ObjVal=ones(1,FoodNumber); for k = 1:FoodNumber ObjVal(k) = fobj(Foods(k,:),input_train,output_train,input_test,output_test); end Fitness=calculateFitness(ObjVal); % 计算适应度函数值 trial=zeros(1,FoodNumber); % 用于记录第i个蜜源有连续trail(i)次没被更新过 % 标记最优蜜源（解） BestInd=find(ObjVal==min(ObjVal)); BestInd=BestInd(end); GlobalMin=ObjVal(BestInd); % 更新全局最优目标函数值 GlobalParams=Foods(BestInd,:); % 更新全局最优参数为最优蜜源 iter=1; % 迭代开始 while ((iter <= maxCycle)) % 循环条件 %%%%%%%%%%%%%%%%%%%%%引领蜂搜索解的过程%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for i=1:(FoodNumber) % 遍历每个蜜源（解） Param2Change=fix(rand*D)+1; % 随机选择需要变异的参数 neighbour=fix(rand*(FoodNumber))+1; % 随机选择相邻蜜源（解）以准备变异 % 需要保证选择的相邻蜜源不是当前蜜源（i） while(neighbour==i) neighbour=fix(rand*(FoodNumber))+1; end sol=Foods(i,:); % 提取当前蜜源（解）对应的的参数 % 参数变异得到新的蜜源：v_{ij}=x_{ij}+\phi_{ij}*(x_{kj}-x_{ij}) sol(Param2Change)=Foods(i,Param2Change)+(Foods(i,Param2Change)-Foods(neighbour,Param2Change))*(rand-0.5)*2; % 确保参数取值范围不越界 ind=find(sol<lb); sol(ind)=lb(ind); ind=find(sol>ub); sol(ind)=ub(ind); % 计算变异后蜜源的目标函数值和适应度函数值 ObjValSol=fobj(sol,input_train,output_train,input_test,output_test); FitnessSol=calculateFitness(ObjValSol); % 更新当前蜜源的相关信息 if (FitnessSol>Fitness(i)) Foods(i,:)=sol; Fitness(i)=FitnessSol; ObjVal(i)=ObjValSol; trial(i)=0; % 如果当前蜜源被更新了，则对应的trial归零 else trial(i)=trial(i)+1; % 如果当前蜜源没有被更新，则trial(i)加1 end end %%%%%%%%%%%%%%%%%%%%%%%% 跟随蜂搜索解的过程 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 计算解（蜜源）的概率 prob=(0.9.*Fitness./max(Fitness))+0.1; % 循环初始化 i=1; t=0; while(t<FoodNumber) % 循环条件 if(rand<prob(i)) % 若随机概率小于当前解（蜜源）的概率 t=t+1; % 循环计数器加1 Param2Change=fix(rand*D)+1; % 随机确定需要变异的参数 neighbour=fix(rand*(FoodNumber))+1; % 随机选择相邻蜜源（解） % 需要保证选择的相邻蜜源不是当前蜜源（i） while(neighbour==i) neighbour=fix(rand*(FoodNumber))+1; end sol=Foods(i,:); % 提取当前蜜源i（解）对应的的参数 % 参数变异得到新的蜜源：v_{ij}=x_{ij}+\phi_{ij}*(x_{kj}-x_{ij}) sol(Param2Change)=Foods(i,Param2Change)+(Foods(i,Param2Change)-Foods(neighbour,Param2Change))*(rand-0.5)*2; % 防止参数越界 ind=find(sol<lb); sol(ind)=lb(ind); ind=find(sol>ub); sol(ind)=ub(ind); % 计算变异后蜜源的目标函数值和适应度函数值 ObjValSol=fobj(sol,input_train,output_train,input_test,output_test); FitnessSol=calculateFitness(ObjValSol); % 更新当前蜜源的相关信息 if (FitnessSol>Fitness(i)) Foods(i,:)=sol; Fitness(i)=FitnessSol; ObjVal(i)=ObjValSol; trial(i)=0; % 如果当前蜜源被更新了，则对应的trial归零 else trial(i)=trial(i)+1; % 如果当前蜜源没有被更新，则trial(i)加1 end end i=i+1; % 更新i if (i==(FoodNumber)+1) % 若值超过蜜源数量，则i重新初始化 i=1; end end % 记住最优蜜源 ind=find(ObjVal==min(ObjVal)); ind=ind(end); if (ObjVal(ind)<GlobalMin) GlobalMin=ObjVal(ind); GlobalParams=Foods(ind,:); end %%%%%%%%%%%% 侦查蜂搜索解的过程 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 找出连续最多次都没有被更新的蜜源 ind=find(trial==max(trial)); ind=ind(end); % 如果连续没有更新的次数大于限定次数，则由侦查蜂重新初始化该蜜源 if (trial(ind)>limit) Bas(ind)=0; sol=(ub-lb).*rand(1,D)+lb; ObjValSol=fobj(sol,input_train,output_train,input_test,output_test); FitnessSol=calculateFitness(ObjValSol); Foods(ind,:)=sol; Fitness(ind)=FitnessSol; ObjVal(ind)=ObjValSol; end iter=iter+1; end % 一次ABC算法完结 BestGlobalMins(r)=GlobalMin; % 记录本次ABC算法的最优目标函数值 BestGlobalParams(r,:)=GlobalParams; % 记录本次ABC算法的最优参数 end % end of runs %% %%%%%%%%%%%%%用ABC算法优化SVR中的参数c和g结束%%%%%%%%%%%%%%%%%%%% %% 打印参数选择结果，这里输出的是最后一次ABC算法寻优得到的参数 bestc=GlobalParams(1); bestg=GlobalParams(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,'*-');legend('预测残差：实际-预测',0); 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]; MAE=mymae(wndspd(104:end),test_pre) MSE=mymse(wndspd(104:end),test_pre) MAPE=mymape(wndspd(104:end),test_pre) %% 显示程序运行时间 toc