How to change it be to multi-objective pso (for FOPID)? Can anyone help me?

The program below is a single-objective pso (for fractional order pid). ( attach file)

How to change it be to multi-objective pso (for FOPID)?

%%%%%% Tunning of FOPID controller using Particle Swarm Optimization 

%%% %% Initialization
clc
clear

n = 15;           % Size of the swarm " no of birds "
bird_setp = 3;   % Maximum number of "birds steps"
dim = 15;          % Dimension of the problem

c2 = 2.1;          % PSO parameter C1 
c1 = 1.2;        % PSO parameter C2 
w =0.9;           % pso momentum or inertia  
fitness=zeros(n,bird_setp);

                                       %-----------------------------%
                                       %    initialize the parameter %
                                       %-----------------------------%
                                       
R1 = rand(dim, n);
R2 = rand(dim, n);
current_fitness =zeros(n,1);

                                 %------------------------------------------------%
                                 % Initializing swarm and velocities and position %
                                 %------------------------------------------------%
current_position = rand(dim, n);
velocity = .3*rand(dim, n) ;
local_best_position  = current_position ;


                                 %-------------------------------------------%
                                 %     Evaluate initial population           %           
                                 %-------------------------------------------%

for i = 1:n

         Kp1 = current_position(1,i);  Ki1 = current_position(2,i); Kd1 = current_position(3,i);
         Kp2 = current_position(4,i);  Ki2 = current_position(5,i); Kd2 = current_position(6,i);
         Kp3 = current_position(7,i);  Ki3 = current_position(8,i); Kd3 = current_position(9,i);
         lam1 = current_position(10,i); mu1 = current_position(11,i);
         lam2 = current_position(12,i); mu2 = current_position(13,i);
         lam3 = current_position(14,i); mu3 = current_position(15,i);

         simopt = simset('solver','ode45');  % Initialize sim options
         [tout,xout,yout] = sim('Model',[0 100],simopt);
          F = sum(e2);
          current_fitness(i) = F;
end


local_best_fitness  = current_fitness ;
[global_best_fitness,g] = min(local_best_fitness) ;

for i=1:n
    globl_best_position(:,i) = local_best_position(:,g) ;
end
                                               %-------------------%
                                               %  VELOCITY UPDATE  %
                                               %-------------------%

velocity = w *velocity + c1*(R1.*(local_best_position-current_position)) + c2*(R2.*(globl_best_position-current_position));

                                               %------------------%
                                               %   SWARMUPDATE    %
                                               %------------------%
                                               
            
current_position = current_position + velocity ;
    for i = 1:n
        for j = 1:n
            if current_position(i,j)>1
                current_position(i,j) = 1;
            elseif current_position(i,j)<0
                current_position(i,j) = 0;
            end
        end
    end
                                               %------------------------%
                                               %  evaluate anew swarm   %
                                               %------------------------%
                                               

%% Main Loop
iter = 0 ;        % Iterations’counter
while  ( iter < bird_setp )
iter = iter + 1;

for i = 1:n
         Kp1 = current_position(1,i);  Ki1 = current_position(2,i); Kd1 = current_position(3,i);
         Kp2 = current_position(1,i);  Ki2 = current_position(2,i); Kd2 = current_position(3,i);
         Kp3 = current_position(1,i);  Ki3 = current_position(2,i); Kd3 = current_position(3,i);
         lam1 = current_position(10,i); mu1 = current_position(11,i);
         lam2 = current_position(12,i); mu2 = current_position(13,i);
         lam3 = current_position(14,i); mu3 = current_position(15,i);
          simopt = simset('solver','ode45');  % Initialize sim options
         [tout,xout,yout] = sim('Model',[0 100],simopt);
         F = sum(e2);
       current_fitness(i) = F;
end


for i = 1 : n
        if current_fitness(i) < local_best_fitness(i)
           local_best_fitness(i)  = current_fitness(i);  
           local_best_position(:,i) = current_position(:,i)   ;
        end   
 end

  
 [current_global_best_fitness,g] = min(local_best_fitness);
  
    
if current_global_best_fitness < global_best_fitness
   global_best_fitness = current_global_best_fitness;
   
    for i=1:n
        globl_best_position(:,i) = local_best_position(:,g);
    end
   
end


 velocity = w *velocity + c1*(R1.*(local_best_position-current_position)) + c2*(R2.*(globl_best_position-current_position));
 current_position = current_position + velocity; 
     for i = 1:n
        for j = 1:n
            if current_position(i,j)>1
                current_position(i,j) = 1;
            elseif current_position(i,j)<0
                current_position(i,j) = 0;
            end
        end
    end
 sprintf('The value of interation iter %3.0f ', iter )


end % end of while loop its mean the end of all step that the birds move it 
                      

            xx=fitness(:,3);
            [Y,I] = min(xx);
            current_position(:,I)

               


    

%