%%%%%基于模拟退火粒子群优化函数——轮盘赌策略迭代
%%%%惯性权重二次缩减
clc
clear
%%%%
%导入训练数据和测试数据
train_input=xlsread('Sample Data2',['AA2:AG111']);
train_out=xlsread('Sample Data2',['Z2:Z111']);
test_input=xlsread('Sample Data2',['AA112:AH141']);
test_out=xlsread('Sample Data2',['Z112:Z141']);
%%归一化
[train_data ,pstrain0] = mapminmax(train_input',0,1);
[test_data] = mapminmax('apply',test_input',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';
%%参数初始化
N=80; %%%种群规模
c1=1.5;
c2=2.0;
c1f=0.5;c1i=2.5; %c1f、c1i表示对应的最终和初始值
c2f=2.5;c2i=0.5; %%%%认知因子
wmax=1.2;
wmin=0.8; %%%%惯性权重的取值范围
xmax=[1000 100];
xmin=[0.1 0.01];%%%%%%位置范围
M=400; %%%%%迭代次数
%MaxC=10; %%%%%混沌搜索步数
D=2; %%%%%问题维数
vmax=[500 50];
vmin=[-500 -50]; %%%%速度范围
lamda=0.8; %%%%退火常数
eps=10^(-8); %机器的浮点运算误差限，当某个数小于它时则记为0
%%定义LSSVM相关参数
type='function estimation'; %回归
%type=’c’，表示classification 分类，type=’f’，表示回归。
kernel='RBF_kernel'; %LSSVM的核函数选取为径向基内核
proprecess='proprecess';
%proprecess表明数据已归一化处理，original表明数据没有进行归一化处理
format long;
%%
%------初始化种群的个体------------
for i=1:N
for j=1:D
x(i,j)=(xmax(1,j)-xmin(1,j))*rand(1,1)+xmin(1,j); %随机初始化位置
v(i,j)=vmax(1,j)*rands(1,1); %随机初始化速度
end
end
%------先计算各个粒子的适应度，并初始化个体适应度P(i)和群体最优适应度Pg----------------------
for i=1:N
%%%%%计算各个粒子的适应度p(i)
gam=x(i,1);sig2=x(i,2);
model=initlssvm(train_data,train_result,type,gam,sig2,kernel,proprecess);
model=trainlssvm(model);
%求出训练集和测试集的预测值
[train_predict_y,zt,model]=simlssvm(model,train_data);
[test_predict_y,zt,model]=simlssvm(model,test_data);
%训练集和测试集的预测数据反归一化
train_predict=mapminmax('reverse',train_predict_y',pstrain1);%训练集预测值
test_predict=mapminmax('reverse',test_predict_y',pstrain1); %测试集预测值
%%%%预测数据四舍五入保留一位有效数字
train_predict=roundn(train_predict,-1);
test_predict=roundn(test_predict,-1);
%计算均方差
result=[train_out;test_out];
predict=[train_predict';test_predict'];
allmse=(sum(((result-predict)./result).^2))/length(result);
fitness(i)=allmse; %以所有数据的预测值计算的均方差为适应度值
p(i)=fitness(i);
y(i,:)=x(i,:);
end
pg = x(N,:); %Pg为全局最优
for i=1:(N-1)
%%%%%计算群体粒子的最优适应度pg
gam=pg(1);sig2=pg(2);
model=initlssvm(train_data,train_result,type,gam,sig2,kernel,proprecess);
model=trainlssvm(model);
%求出训练集和测试集的预测值
[train_predict_y,zt,model]=simlssvm(model,train_data);
[test_predict_y,zt,model]=simlssvm(model,test_data);
%训练集和测试集的预测数据反归一化
train_predict=mapminmax('reverse',train_predict_y',pstrain1);%训练集预测值
test_predict=mapminmax('reverse',test_predict_y',pstrain1); %测试集预测值
%%%%预测数据四舍五入保留一位有效数字
train_predict=roundn(train_predict,-1);
test_predict=roundn(test_predict,-1);
%计算均方差
result=[train_out;test_out];
predict=[train_predict';test_predict'];
allmse=(sum(((result-predict)./result).^2))/length(result);
global_fitness=allmse; %以所有数据的预测值计算的均方差为适应度值
if p(i)
pg=x(i,:);
end
end
%%
%------进入主要循环，按照公式依次迭代------------
T=global_fitness/log(5); %%%设置初始温度
for t=1:M
groupFit=global_fitness %%%%群体最优
%%%根据轮盘赌策略公式确定当前温度下各粒子的适配值
for i=1:N
Tfit(i) = exp( - (p(i) - groupFit)/T);
end
SumTfit = sum(Tfit);
Tfit = Tfit/SumTfit;
pBet = rand();
for i=1:N
ComFit(i)=sum(Tfit(1:i));
if pBet<=ComFit(i)
pg_plus=x(i,:);
break;
end
end
c1f=c1i+(c1f-c1i)*t/M;%%c1f递减
c2f=c2i+(c2f-c2i)*t/M; %%c2f递增
C=c1+c2;
ksi=2/abs(2-C-sqrt(C^2-4*C));
%%%%%%%更新粒子的速度与位置
for i=1:N
v(i,:)=ksi*(v(i,:)+c1*rand*(y(i,:)-x(i,:))+c2*rand*(pg_plus-x(i,:)));
x(i,:)=x(i,:)+v(i,:);
%%%限制粒子的速度与位置的大小在取值范围之内
for j=1:D
if v(i,j)>vmax(1,j)
v(i,j)=vmax(1,j);
end
if v(i,j)
v(i,j)=vmin(1,j);
end
if x(i,j)>xmax(1,j)
x(i,j)=xmax(1,j);
end
if x(i,j)
x(i,j)=xmin(1,j);
end
end
%%%%%%%计算更新后粒子的适应度
gam=x(i,1);sig2=x(i,2);
model=initlssvm(train_data,train_result,type,gam,sig2,kernel,proprecess);
model=trainlssvm(model);
%求出训练集和测试集的预测值
[train_predict_y,zt,model]=simlssvm(model,train_data);
[test_predict_y,zt,model]=simlssvm(model,test_data);
%训练集和测试集的预测数据反归一化
train_predict=mapminmax('reverse',train_predict_y',pstrain1);%训练集预测值
test_predict=mapminmax('reverse',test_predict_y',pstrain1); %测试集预测值
%%%%预测数据四舍五入保留一位有效数字
train_predict=roundn(train_predict,-1);
test_predict=roundn(test_predict,-1);
%计算均方差
result=[train_out;test_out];
predict=[train_predict';test_predict'];
allmse=(sum(((result-predict)./result).^2))/length(result);
fitness(i)=allmse; %以所有数据的预测值计算的均方差为适应度值
local_fit(i)=fitness(i);
if local_fit(i)
p(i)=local_fit(i);
y(i,:)=x(i,:);
end
if p(i) pg=y(i,:);
end
end
T=T*lamda;%%%%%%降温
yy(t)=groupFit;
end
xm=pg'
fv=groupFit
hold on
plot(yy,'g')
hold off
%title(['适应度曲线 ' '终止代数＝' num2str(M)]);
%xlabel('进化代数');ylabel('适应度');
