%  Kmean Optimizer source code (Developed in MATLAB)
%
%
%  paper:A new metaheuristic optimization based on K-Means clustering algorithm and its application for...
...structural damage identification in a complex 3D concrete structure
%  https://doi.org/10.1016/j.knosys.2022.109189 
 
%  Hoang-Le Minh,Thanh Sang-To, Magd Abdel Wahab, Thanh Cuong-Le
%
%  Email: hoang.lm@ou.edu.vn, minhhoang.le@ugent.be, magd.abdelwahab@ugent.be, cuong.lt@ou.edu.vn
%  Programmer: Hoang-Le Minh
%   
%  Knowledge-Based Systems
%% ---- PROGRAM FOR 23 BENCHMARK FUNCTIONS ---
clc
clear all
close all
format short e
%% SELECT FUNTION OF 23 CLASSICAL BENCHMARK FUNCTIONS---
fun = 2; % select function from 1 to 23
disp(['Function ' num2str(fun)])
[Lb,Ub,dim,fobj]=Get_Functions_details(fun);
%% INITIAL DEFINATION
nvar = dim;                 % Dimension  
Lb = Lb.*ones(1,nvar);      % Lower boundary
Ub = Ub.*ones(1,nvar);      % Uper boundary
number_centroid = 3;        % Number of Centroid
NP = 30;                    % Population size
cycle = 1000;               % Iteration
%% -------------------------KO_OPTIMZER ------------------------------------ 
%% INITIAL SOLUTION 
for q =1:NP  
    particle_kmean(q).pos = Lb + rand(1,nvar).*(Ub-Lb);
    particle_kmean(q).cost = fobj(particle_kmean(q).pos); 
end
SA =  particle_kmean; % create initial solution for self-improvement using Algorithm 3
[value, index_min] = min([particle_kmean.cost]);
[value, index_max] = max([particle_kmean.cost]);
G_best_kmean = particle_kmean(index_min); % find best solution
G_worst_kmean = particle_kmean(index_max); % find worst solution
centroid_matrix = rand(NP,nvar); % creare random initial centroids matrix
[~,Centroid_matrix] = kmeans(centroid_matrix,3);% create random centroids position
%% MAIN LOOP
best_kmean = zeros(cycle,1);
for k = 1:cycle
        number_sort = round(((NP-4)*k/(1-cycle)) + (NP- (NP-4)/(1-cycle))); % Reduction of population size (Eq.6)
        a =1 -(1/cycle)*k; % define value for calculating (Eq.14)
        Mean_centroid = (Centroid_matrix(1,:) + Centroid_matrix(2,:) + Centroid_matrix(3,:))/3;% Mean centroid (Eq.12)
        for i = 1:NP
         distance(i) = norm(SA(i).pos - G_best_kmean.pos);
         far_distance = norm(G_worst_kmean.pos - G_best_kmean.pos);
         D(i) = distance(i)./(far_distance+eps); % relative distance (Eq.9)
         CL(i) = exp(-D(i)); % characteristic length (Eq.10)
         PDP(i) = CL(i).*sin(k*pi/cycle); % position density probability (Eq.11)
         trend = randi(3); % select random trend for each seach agent 
         if  PDP(i)> 0.5 % move following the ability of exploitation using Algorithm 4
          switch trend
              case 1 % The first trend for exploitation (Eq.13)
                [d1 d2] =drift_control(cycle,a); % drift control in Case 1 for each search agent (Eq.14), (Eq.15), (Eq.16)
                SA(i).pos = 2*G_best_kmean.pos.*rand + d1.*abs(-G_best_kmean.pos.*(rand(1,nvar))+ SA(i).pos)+...
                                        d2.*abs(-Mean_centroid.*(rand(1,nvar)) + SA(i).pos);   
              case 2 % The second trend for exploitation (Eq.17)
                SA(i).pos =  2*SA(i).pos.*rand +  (rand(1,nvar)).*(G_best_kmean.pos - (rand(1,nvar)).*SA(i).pos) + ...
                                                (rand(1,nvar)).*(Centroid_matrix(1,:) - (rand(1,nvar)).*SA(i).pos)+...
                                                (rand(1,nvar)).*(Centroid_matrix(2,:) - (rand(1,nvar)).*SA(i).pos)+...
                                                (rand(1,nvar)).*(Centroid_matrix(3,:) - (rand(1,nvar)).*SA(i).pos);
               case 3 % The third trend for exploitation (Eq.18)
               SA(i).pos = 2*Mean_centroid.*rand  + (rand(1,nvar)).*(G_best_kmean.pos-(rand(1,nvar)).*Mean_centroid);
          end % end switch  
               SA(i).pos =  boundary_condition_KO(SA(i).pos,Lb,Ub,dim);                 
               SA(i).cost =fobj(SA(i).pos); 
            % Self-improvement of the initial solution using Algorithm 3
              if  SA(i).cost <  particle_kmean(i).cost
               particle_kmean(i) = SA(i);
              end    
              % update best solution using Algorithm 3
              if particle_kmean(i).cost < G_best_kmean.cost
                G_best_kmean = particle_kmean(i);
              end
              % update worst solution using Algorithm 3
              if particle_kmean(i).cost > G_worst_kmean.cost
                G_worst_kmean = particle_kmean(i);
              end              
          else  % move following the ability of exploration using Algorithm 4
             if rand>0.5
              SA(i).pos =  Mean_centroid + levy_fun_KO(1,nvar,(1+rand)); % the first trend of exploration (Eq.24)
             else
             SA(i).pos = Lb  + (Ub-Lb).*rand(1,nvar); % the second trend of exploration (Eq.25)  
             end
               SA(i).pos =  boundary_condition_KO(SA(i).pos,Lb,Ub,dim);               
               SA(i).cost =fobj(SA(i).pos);               
              % Self-improvement of the initial solution using Algorithm 3
               if  SA(i).cost <  particle_kmean(i).cost
                particle_kmean(i) = SA(i);
               end  
                 %  update best solution using Algorithm 3
               if particle_kmean(i).cost < G_best_kmean.cost
                G_best_kmean = particle_kmean(i);
               end
                % update worst solution using Algorithm 3
               if particle_kmean(i).cost > G_worst_kmean.cost  
                G_worst_kmean = particle_kmean(i);
               end             
         end   
        end 
          %% UPDATE CENTROIND VECTOR POSITION FOR THE NEXT ITERATION USING ALGORITHM 2
          [value_sort, index_sort] = sort([particle_kmean.cost]); % machenistic sort (Eq.4)
           matrix_pos_number_sort = zeros(number_sort,nvar);
          for j = 1:number_sort 
              matrix_pos_number_sort(j,:) = particle_kmean(index_sort(j)).pos; % Matrix assembly (Eq.7)
          end
          [~,Centroid_matrix] = kmeans( matrix_pos_number_sort,3); % Update centroid vector position(Eq.8)
          best_kmean(k,:) = G_best_kmean.cost; % write output
end
disp([' the best solution =' num2str(G_best_kmean.pos)])
disp([' the best fitness =' num2str(G_best_kmean.cost)])
figure (1)
plot(best_kmean(1:cycle),'-r','lineWidth',2);
xlabel('Iteration','FontName','Times New Roman','FontSize',16)
ylabel('Best fitness','FontName','Times New Roman','FontSize',16)
set(gca, 'FontName', 'Times New Roman','FontSize', 16)
legend('KO optimizer','FontName','Times New Roman','FontSize',16)
 
