%  Termite life cycle optimizer source code (Developed in MATLAB)
%
%  programming: Hoang-Le Minh
%  Soete Laboratory, Department of Electrical Energy, Metals, Mechanical Constructions, and Systems, 
%  Faculty of Engineering and Architecture, Ghent University, 9000 Gent, Belgium
%
%  paper: Termite life cycle optimizer
%    
%  Hoang-Le Minh,Thanh Sang-To,Guy Theraulazc, Magd Abdel Wahab, Thanh Cuong-Le
%
%  Email: hoang.lm@ou.edu.vn, minhhoang.le@ugent.be, cuong.lt@ou.edu.vn,magd.abdelwahab@ugent.be
%
%  Expert Systems with Applications
clc
clear all
close all
format short e
%% ---- PROGRAM FOR 23 BENCHMARK FUNCTIONS ---
fun = 8; % Select function from 1 to 23
[lb,ub,dim,fobj]=Get_Functions_details(fun);
%% INITIAL DEFINATION
nvar = dim;                  % Dimension  
Lower_B = lb;                % Lower boundary
Uper_B= ub;                  % Uper boundary
lb = Lower_B.*ones(1,nvar);  % Lower boundary (vector)
ub = Uper_B.*ones(1,nvar);   % Uper boundary (vector)
NP = 30;                     % Population size
cycle = 1000;                % Iteration
trial_worker = zeros(NP,1);  % Define according to Algorithm 2
trial_soldier = zeros(NP,1); % Define according to Algorithm 2
%% TERMITE LIFE CYCLE OPTIMIZER ------------------------------------------
empt.pos =[];          
empt.cost = [];         
particle_termine = repmat(empt,NP,1); % Define termite swarms
Int_Population = repmat(empt,NP,1);   % Define initial population
%% Creare initial poppulation
for i = 1:NP
    Int_Population(i).pos  = lb + rand(1,nvar).*(ub-lb);
    Int_Population(i).cost = fobj(Int_Population(i).pos);
end
% Assign the termite caste
 particle_soldier = Int_Population; 
 particle_worker  = Int_Population;
 particle_reproductive = Int_Population;
% Find the best solution
[value, index] = min([Int_Population.cost]);
Gbest = Int_Population(index);
%% Main iterations
best_TLCO = zeros(cycle,1);
for k = 1:cycle   
    % Calculate the proposed beta parameter for modification step length S
      beta = (0.5/cycle)*k +1.5; % Equation 6   
    % Calculate the timing of reproductive termites’ emergence 
      sigma = 1/(0.1*cycle); % Equation 17
      lamda_soldier = 1/(1+exp(-sigma*(k-0.5*cycle))); % Equation 16
      lamda_worker = 1- 1/(1+exp(-sigma*(k-0.5*cycle))); % Equation 15
    for i = 1:round(NP)
%% The phase of Worker termite 
        if i <= round(0.7*NP)
            % Update the new position using Eq.12
             particle_worker(i).pos =  particle_worker(i).pos + ...
                                       (-1+2*rand(1)).*(rand(1,nvar) + levy_fun_TLCO(1,nvar,beta)).*abs(Gbest.pos - particle_worker(i).pos);
             particle_worker(i).pos = boundary_condition(particle_worker(i).pos,lb,ub,dim);
             particle_worker(i).cost = fobj( particle_worker(i).pos);
            % The process of self-improvement using Algorithm 1
             if particle_worker(i).cost < Int_Population(i).cost;
                 Int_Population(i) = particle_worker(i);
            % End of Algorithm 1
               if  Int_Population(i).cost < Gbest.cost 
                 Gbest = Int_Population(i); % record the best solution
               end
             else % Recording the wasteful movement using Algorithm 2
                trial_worker(i) =  trial_worker(i)+1;
             end
% The phase of reproductive termites using algorithm 3
                    if (trial_worker(i))/cycle >  lamda_worker
                            trial_worker(i) = 0;
                            A=randperm(NP);
                            A = [1:i-1, i+1:NP];
                            r1 = randi(numel(A));  
                            A(r1)=[];
                            r2 = randi(numel(A));
                 % Update the new position using Eq.20
                            if Int_Population(r2).cost<Int_Population(r1).cost
                                particle_reproductive(i).pos  = Int_Population(r1).pos + rand(1,nvar).*(Int_Population(r2).pos - Int_Population(r1).pos);
                            else
                                particle_reproductive(i).pos  =Int_Population(r1).pos - rand(1,nvar).*(Int_Population(r2).pos - Int_Population(r1).pos);
                            end  
                              particle_reproductive(i).pos = boundary_condition(particle_reproductive(i).pos,lb,ub,dim);
                             particle_reproductive(i).cost = fobj( particle_reproductive(i).pos);
                 % The process of self-improvement using Algorithm 1                             
                             if particle_reproductive(i).cost < Int_Population(i).cost
                                 Int_Population(i) = particle_reproductive(i);
                 % End of Algorithm 1                                 
                                   if  Int_Population(i).cost < Gbest.cost 
                                      Gbest = Int_Population(i); % record the best solution
                                   end
                             else
                 % Update the new position using Eq.21               
                            particle_reproductive(i).pos = lb + rand(1,nvar).*(ub-lb);
                             end
                    particle_worker(i) = particle_reproductive (i);
% End of Algorithm 3   
                    end          
        end % End of the termite worker phase
%% -------------------------------------------------------------------       
%% The phase of Soldier termite 
        if  i >round(0.7*NP) && i <= round(NP)
            % Update the new position using Eq.14
             particle_soldier(i).pos =2*rand(1).*Gbest.pos +...
                                      (-1+2*rand(1)).*abs( particle_soldier(i).pos - (levy_fun_TLCO(1,nvar,beta).*Gbest.pos)) ;  
             particle_soldier(i).pos = boundary_condition(particle_soldier(i).pos,lb,ub,dim);
             particle_soldier(i).cost = fobj( particle_soldier(i).pos); 
             % The process of self-improvement using Algorithm 1            
            if  particle_soldier(i).cost < Int_Population(i).cost;
                 Int_Population(i) = particle_soldier(i); 
            % End of Algorithm 1                 
            if Int_Population(i).cost < Gbest.cost 
                Gbest =  Int_Population(i);
            end
            else % Recording the wasteful movement using Algorithm 2
              trial_soldier(i) =  trial_soldier(i)+1 ;
            end
% The phase of reproductive termites using algorithm 3           
                    if (trial_soldier(i))/cycle> lamda_soldier
                             trial_soldier(i)=0;
                             A = randperm(NP,0.7*NP);
                             A = [0.7*NP:i-1, i+1:NP];
                             r1 = randi(numel(A)); % 
                             A(r1)=[];
                             r2 = randi(numel(A));
                           % Update the new position using Eq.20
                            if Int_Population(r2).cost<Int_Population(r1).cost
                                particle_reproductive(i).pos  = Int_Population(r1).pos + rand(1,nvar).*(Int_Population(r2).pos - Int_Population(r1).pos);
                            else
                                particle_reproductive(i).pos  = Int_Population(r1).pos - rand(1,nvar).*(Int_Population(r2).pos - Int_Population(r1).pos);
                            end  
                             particle_reproductive(i).pos = boundary_condition(particle_reproductive(i).pos,lb,ub,dim);
                             particle_reproductive(i).cost = fobj( particle_reproductive(i).pos);
                         % The process of self-improvement using Algorithm 1                             
                             if particle_reproductive(i).cost < Int_Population(i).pos
                                 Int_Population(i) = particle_reproductive(i);
                         % End of Algorithm 1 
                                   if  Int_Population(i).cost < Gbest.cost 
                                      Gbest = Int_Population(i);% record the best solution
                                   end
                             else
                           % Update the new position using Eq.21       
                                       particle_reproductive(i).pos = lb + rand(1,nvar).*(ub-lb);
                             end
                    particle_soldier(i) = particle_reproductive(i);
                    end    
% End of Algorithm 3 
        end % End of termite soldier phase
end 
    best_TLCO(k) =  Gbest.cost;
end % end of cycle
disp ('==================================================================')
disp([ ' Best_Int_Population_TLCO = ' num2str(Gbest.pos)])
disp([ ' Best_fitness_TLCO  = ' num2str(Gbest.cost)])
figure(1)
semilogy(best_TLCO(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('TLCO','FontName','Times New Roman','FontSize',16)

 