模拟退火粒子群优化算法控制程序

%%%%%基于模拟退火粒子群优化函数——轮盘赌策略迭代
%%%%惯性权重二次缩减

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('适应度');