%%==================清空环境=============
clc
clear
%%导入训练数据和测试数据
load data1.mat
%其中train_data test_data分别为训练集和测试集
%train_result test_result为训练结果和测试结果
%%数据归一化
%%归一化方法1（利用libsvm工具箱函数归一化）
[train_x,test_x]=scaleForSVM(train_data,test_data,0,1)
[train_y,test_y,pstrain1]=scaleForSVM(train_result,test_result,0,1)
%%归一化方法2（利用svm自带mapminmax函数归一化）
%[train_data ,pstrain0] = mapminmax(train',0,1);
%[test_data] = mapminmax('apply',test',pstrain0);
%[train_result,pstrain1] = mapminmax(train_out,0,1);
%[test_result] = mapminmax('apply',test_out,pstrain1);
%train_data = train_data'
%train_result=train_result'
%test_data = test_data'
%% 参数初始化
%粒子群算法中的两个参数
c1 = 1.5; % c1 belongs to [0,2] c1:初始为1.5,pso参数局部搜索能力
c2 = 1.5; % c2 belongs to [0,2] c2:初始为1.7,pso参数全局搜索能力
maxgen=300; % 进化次数
sizepop=30; % 种群规模
popcmax=10^(3); % popcmax:初始为1000,SVM 参数c的变化的最大值.
popcmin=10^(-1); % popcmin:初始为0.1,SVM 参数c的变化的最小值.
popgmax=10^(2); % popgmax:初始为1000,SVM 参数g的变化的最大值
popgmin=10^(-1); % popgmin:初始为0.01,SVM 参数c的变化的最小值.
k = 0.5; % k belongs to [0.1,1.0];
Vcmax = k*popcmax;%参数 c 迭代速度最大值
Vcmin = -Vcmax ;
Vgmax = k*popgmax;%参数 g 迭代速度最大值
Vgmin = -Vgmax ;
eps = 10^(-7);
%%定义lssvm相关参数
type = 'function estimation';
kernel = 'RBF_kernel';
proprecess='original'
%% 产生初始粒子和速度
for i=1:sizepop
% 随机产生种群
pop(i,1) = (popcmax-popcmin)*rand(1,1)+popcmin ; % 初始种群
pop(i,2) = (popgmax-popgmin)*rand(1,1)+popgmin;
V(i,1)=Vcmax*rands(1,1); % 初始化速度
V(i,2)=Vgmax*rands(1,1);
% 计算初始适应度
gam = pop(i,1)
sig2 = pop(i,2)
model=initlssvm(train_x,train_y,type,gam,sig2,kernel,proprecess)
model=trainlssvm(model) % 训练svm模型
%求出训练集和测试集的预测值
[ptrain,zt,model]=simlssvm(model,train_x);
%预测数据反归一化
% train_predict=mapminmax('reverse',train_predict_y,pstrain1)%训练集预测值
%计算均方差
trainmse=sum((ptrain-train_y).^2)/length(train_y)
fitness(i)=trainmse%以训练集的预测值计算的均方差为适应度值
end
% 找极值和极值点
[global_fitness bestindex]=min(fitness) % 全局极值
local_fitness=fitness % 个体极值初始化
global_x=pop(bestindex,:) % 全局极值点
local_x=pop % 个体极值点初始化
% 每一代种群的平均适应度
avgfitness_gen = zeros(1,maxgen);
tic
%% 迭代寻优
for i=1:maxgen
for j=1:sizepop
%速度更新
wV = 1; % wV best belongs to [0.8,1.2]为速率更新公式中速度前面的弹性系数
V(j,:) = wV*V(j,:) + c1*rand*(local_x(j,:) - pop(j,:)) + c2*rand*(global_x - pop(j,:));
if V(j,1) > Vcmax %以下几个不等式是为了限定速度在最大最小之间
V(j,1) = Vcmax;
end
if V(j,1) < Vcmin
V(j,1) = Vcmin;
end
if V(j,2) > Vgmax
V(j,2) = Vgmax;
end
if V(j,2) < Vgmin
V(j,2) = Vgmin; %以上几个不等式是为了限定速度在最大最小之间
end
%种群更新
wP = 1; % wP:初始为1,种群更新公式中速度前面的弹性系数
pop(j,:)=pop(j,:)+wP*V(j,:);
if pop(j,1) > popcmax %以下几个不等式是为了限定 c 在最大最小之间
pop(j,1) = popcmax;
end
if pop(j,1) < popcmin
pop(j,1) = popcmin;
end
if pop(j,2) > popgmax %以下几个不等式是为了限定 g 在最大最小之间
pop(j,2) = popgmax;
end
if pop(j,2) < popgmin
pop(j,2) = popgmin;
end
% 自适应粒子变异
if rand>0.5
k=ceil(2*rand);%ceil 是向离它最近的大整数圆整
if k == 1
pop(j,k) = (20-1)*rand+1;
end
if k == 2
pop(j,k) = (popgmax-popgmin)*rand+popgmin;
end
%%新粒子适应度值
gam=pop(j,1)
sig2=pop(j,2)
model=initlssvm(train_x,train_y,type,gam,sig2,kernel,proprecess)
model=trainlssvm(model) % 训练svm模型
%求出训练集和测试集的预测值
[ptrain,zt,model]=simlssvm(model,train_x)
%计算均方差
trainmse=sum((ptrain-train_y).^2)/length(train_y)
fitness(j)=trainmse%以训练集的预测值计算的均方差为适应度值
end
%个体最优更新
if fitness(j) < local_fitness(j)
local_x(j,:) = pop(j,:);
local_fitness(j) = fitness(j);
end
if fitness(j) == local_fitness(j) && pop(j,1) < local_x(j,1)
local_x(j,:) = pop(j,:);
local_fitness(j) = fitness(j);
end
%群体最优更新
if fitness(j) < global_fitness
global_x = pop(j,:);
global_fitness = fitness(j);
end
if abs( fitness(j)-global_fitness )<=eps && pop(j,1) < global_x(1)
global_x = pop(j,:);
global_fitness = fitness(j);
end
end
fit_gen(i)=global_fitness;
avgfitness_gen(i) = sum(fitness)/sizepop;
end
toc
%% 结果分析
plot(fit_gen,'LineWidth',5);
title(['适应度曲线','(参数c1=',num2str(c1),',c2=',num2str(c2),',终止代数=',num2str(maxgen),')'],'FontSize',13);
xlabel('进化代数');ylabel('适应度');
bestc = global_x(1)
bestg = global_x(2)
gam=bestc
sig2=bestg
model=initlssvm(train_data,train_result,type,gam,sig2,kernel,proprecess)
model=trainlssvm(model)
%求出训练集和测试集的预测值
[train_predict,zt,model]=simlssvm(model,train_data)
[test_predict,zt,model]=simlssvm(model,test_data)
%预测数据反归一化
%train_predict=mapminmax('reverse',train_predict_y,pstrain1)
%test_predict=mapminmax('reverse',test_predict_y,pstrain1)
%计算均方差
trainmse=sum((train_predict-train_result).^2)/length(train_result)
testmse=sum((test_predict-test_result).^2)/length(test_result)
test00=sum(abs(test_predict-test_result)./test_result)/length(test_predict)
test0=(test_predict-test_result)./test_result
- 1
- 2
- 3
- 4
- 5
- 6
前往页