how to write ant colony optimization code

Question

i have a function

function t=strength(x)
  
    t(1) = -0.804-0.001.*x(1)-0.001.*x(2)+0.0.*x(3)+0.0.*x(4)-0.031.*x(5)+0.005.*x(6)-0.002.*x(7)+0.001.*x(8)-0.735.*x(9)+0.097.*x(10)+0.012.*x(11)+0.078.*x(12)-0.220.*x(13);
end

with 13 variables and their lower and upper bounds are

lb = [490 0 0 19 0 170 700 990 0.35 1.5 16.8 18.5 1] 
ub = [500 0 0 22 0 180 710 1000 0.35 1.5 16.8 18.5 1]

i need to perform ACO to minimize the function and get t(1). this is my code so far. it may be complete wrong or half wrong. I have just started learning matlab and its not been easy for me. it would be very helpful if you guys could help me completing the code.

clear all
close all
% problem definition
obj = @(x)strength(x)
nVar = 13; % number of parameter
min = [490 0 0 19 0 170 700 990 0.35 1.5 16.8 18.5 1] 
max = [500 0 0 22 0 180 710 1000 0.35 1.5 16.8 18.5 1]
%parameters of ACO
MaxIt=500;      % Maximum Number of Iterations
nAnt=50;        % Number of Ants (Population Size)
Q=1;
tau0=10;        % Initial Phromone
alpha=0.3;      % Phromone Exponential Weight
rho=0.1; 
tau = tau0 * ones(nVar); % Phromone matirx 
Bestfitness=zeros(MaxIt,1);
% Empty Ant
empty_ant.Tour=[];
empty_ant.fitness=[];
% Ant Colony Matrix
ant=repmat(empty_ant,nAnt,1);
BestSol.fitness=inf;
% Main loop of ACO 
for it=1:MaxIt
    
    % Move Ants
    for k=1:nAnt
        
        ant(k).Tour=[];
        
        for l=1:nVar
            
            P=tau(:,l).^alpha;
            
            P(ant(k).Tour)=0;
            
            P=P/sum(P);
            
            j=RouletteWheelSelection(P);
            
            ant(k).Tour=[ant(k).Tour j];
            
        end
        
        ant(k).fitness=fitnessFunction(ant(k).Tour);
        
        if ant(k).fitness<BestSol.fitness
            BestSol=ant(k);
        end
        
    end
    
    % Update Phromones
    for k=1:nAnt
        
        tour=ant(k).Tour;
        
        for l=1:nVar
            
            tau(tour(l),l)=tau(tour(l),l)+Q/ant(k).fitness;
            
        end
        
    end
      % Evaporation
    tau=(1-rho)*tau;
    
    % Store Best Cost
    Bestfitness(it)=BestSol.fitness;
    
    % Show Iteration Information
    disp(['Iteration ' num2str(it) ': Best fitness = ' num2str(Bestfitness(it))]);
    % Plot Solution
    figure(1);
    PlotSolution(BestSol.Tour,model);
    pause(0.01);
    
% Results
figure;
plot(BestCost,'LineWidth',2);
xlabel('Iteration');
ylabel('Best fitness');
grid on;
end