clc
clear
t1=cputime;
%declare the parameters of the optimization
max_iterations  = 30;
no_of_particles = 25;
dimensions = 2;

x_min=0.5;
x_max=100;


%initialise the particles
particle_position=unifrnd(x_min,x_max,no_of_particles,dimensions);

p_best=particle_position;%初始化pbest
%initialize the p_best_fitness array
for count = 1:no_of_particles
    p_best_fitness(count) = 10;%
end
%main QPSO
for count = 1:max_iterations
    %find the fitness of each particle
    %change fitness function as per equation requiresd and dimensions
    for count_x = 1:no_of_particles
       %for count_y=1:dimensions
       parameter=[particle_position(count_x,1);particle_position(count_x,2)]; 
       [a1,a2,a3]=mysvm1(parameter);
       current_fitness(count_x) =a1 ; %計(jì)算粒子的適應(yīng)度；                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                               
    end
    
    
    %decide on p_best etc for each particle
    for count_x = 1:no_of_particles
        if current_fitness(count_x)< p_best_fitness(count_x)
            p_best_fitness(count_x) = current_fitness(count_x);
            p_best(count_x,1:dimensions) = particle_position(count_x,1:dimensions);
        end
    end
    m_best=sum(p_best)/no_of_particles;
    
    %decide on the global best among all the particles
    [g_best_val,g_best_index] = min(current_fitness);%最小適應(yīng)度值，位置
    bestfitness(count)=g_best_val;%每次迭代的最優(yōu)適應(yīng)值
    %g_best contains the position of teh global best
    g_best=particle_position(g_best_index,1:dimensions);      
    %G_temp(count)=g_best;%每次迭代的GBEST;

    %update the position 
    W(count)=0.5-0.5*(max_iterations-count)/max_iterations;
    for count_x = 1:no_of_particles
        for count_y = 1:dimensions
            % 更新位置
            f=rand;
            pp=f*p_best(count_x,count_y)+(1-f)*g_best(1,count_y);
            u=rand;
            v=-log(u);
            if(u>0.5)
                temp_position = pp-W(count)*abs(m_best(1,count_y)-particle_position(count_x,count_y))*v;
            else
                temp_position = pp+W(count)*abs(m_best(1,count_y)-particle_position(count_x,count_y))*v;
            end
            if(temp_position >x_max)
                particle_position(count_x,count_y)=x_max;
            elseif(temp_position <x_min)
                particle_position(count_x,count_y)=x_min;
            else
                particle_position(count_x,count_y)=temp_position;
            end
        end
    end
end
current_fitness(g_best_index)
%---------------------------------------------------
t2=cputime-t1
% 結(jié)果作圖 
figure(1)                               
t=1:30 ;
plot(t,bestfitness(t)-0.2)    
xlabel('Iterations ');
ylabel('Fitness')
[Error,Yt,Yd]=mysvm1(parameter);
e=Yd-Yt;
e=e.^2;
rsme=sqrt(mean(e))
figure(2)
plot(1:length(Yt),Yt,'rd-',1:length(Yd),Yd,'b.:')
%axis([0,50,0.3,1]);
%title('+為真實(shí)值，.為預(yù)測(cè)值')
legend('實(shí)際測(cè)量值','QPSO-LS-SVM輸出值');
%figure(3)
%plot(1:length(Yt),Error,'b.:')
%title('誤差曲線(xiàn)')
%end