How to make change in the algorithm to start and end at a same point?

Hey Everyone,

Please help me to make a logical change in this algorithm. I want to run this algorithm to find the best possible route such that the route will start and end at city 1. While there are 27 cities and the distance matrix represents the distances between them.

distMatrix = [0 25 4.2 4 23 6 23 1.4 9.5 2.3 27 26 28 2.5 8 6.3 22 5.4 22 31 20 1.8 25 7.8 5.7 28 9.8;
     25 0 24 23 12 21 13 25 17 24 6.8 18 16 23 18 20 10 25 11 5.7 5.4 25 18 18 7.6 4.3 17;
     4.2 24 0 0.2 20 4.7 17 3.1 4.8 3.9 25 30 39 2 4.3 2.6 20 8.4 20 29 18 2.7 23 4 3.1 25 4.9;
     4 23 0.2 0 20 4.8 17 2.7 4.9 3.5 25 30 39 1.6 4.4 2.8 20 7.7 20 29 18 2.3 23 4.1 3.2 25 5;
     23 12 20 20 0 20 24 24 16 23 9.1 30 28 23 17 21 1.9 24 1.7 17 7.7 25 6 16 9.1 11 16;
     6 21 4.7 4.8 20 0 18 5.1 5.2 3.9 22 38 37 3.4 3.2 1.6 17 5 21 26 15 5.5 21 3.7 1 23 5;
     23 13 17 17 24 18 0 22 12 21 19 18 17 19 14 17 23 22 23 7 18 22 30 14 0.4 17 13;
     1.4 25 3.1 2.7 24 5.1 22 0 9.7 0.5 26 27 29 1.5 7 5.4 21 4.7 25 30 19 0.8 24 7.6 4.8 27 9.2;
     9.5 17 4.8 4.9 16 5.2 12 9.7 0 8.2 19 35 33 6.4 1.9 4.6 14 9.5 14 23 12 9.7 17 1.1 3.8 19 0.5;
     2.3 24 3.9 3.5 23 3.9 21 0.5 8.2 0 25 27 40 2 6.3 4.6 20 4.7 25 29 18 1.3 24 6.8 4.1 26 8.1;
     27 6.8 25 25 9.1 22 19 26 19 25 0 25 23 25 20 22 7.5 27 7.9 12 10 27 15 19 12 2.4 19;
     26 18 30 30 30 38 18 27 35 27 25 0 2 28 36 38 28 29 28 12 23 27 36 36 18 22 35;
     28 16 39 39 28 37 17 29 33 40 23 2 0 30 35 37 27 31 27 11 22 29 35 35 17 21 34;
     2.5 23 2 1.6 23 3.4 19 1.5 6.4 2 25 28 30 0 5.8 4 20 6.1 24 29 18 0.7 23 5.3 3.5 25 6.3 ;
     8 18 4.3 4.4 17 3.2 14  7 1.9 6.3 20 36 35 5.8 0 2.7 15 7.4 16 24 13 7.8 18 0.9 2.3 21 1.8;
     6.3 20 2.6 2.8 21 1.6 17 5.4 4.6 4.6 22 38 37 4 2.7 0 17 5.7 18 26 15 4.7 20 3.5 0.4 23 4.4;
     22 10 20 20 1.9 17 23 21 14 20 7.5 28 27 20 15 17 0 31 1.2 16 6.1 22 7.8 15 7.6 9.8 14;
     5.4 25 8.4 7.7 24 5 22 4.7 9.5 4.7 27 29 31 6.1 7.4 5.7 31 0 24 31 20 5.4 29 7.9 1.7 27 9.5;
     22 11 20 20 1.7 21 23 25 14 25 7.9 28 27 24 16 18 1.2 24 0 16 6.5  25 7.3 15 7.9 11 14;
     31 5.7 29 29 17 26 7 30 23 29 12 12 11 29 24 26 16 31 16 0 11 31 23 23 7 10 23 ;
     20 5.4 18 18 7.7 15 18 19 12 18 10 23 22 18 13 15 6.1 20 6.5 11 0 20 14 13 2.3 7.6 12;
     1.8 25 2.7 2.3 25 5.5 22 0.8 9.7 1.3 27 27 29 0.7 7.8 4.7 22 5.4 25 31 20 0 25 6 5.6 27 7;
     25 18 23 23 6 21 30 24 17 24 15 36 35 23 18 20 7.8 29 7.3 23 14 25 0 18 15 17 17;
     7.8 18 4 4.1 16 3.7 14 7.6 1.1 6.8 19 36 35 5.3 0.9 3.5 15 7.9 15 23 13 6 18 0 2.8 20 1.2;
     5.7 7.6 3.1 3.2 9.1 1 0.4 4.8 3.8 4.1 12 18 17 3.5 2.3 0.4 7.6 1.7 7.9 7 2.3 5.6 15 2.8 0 9.9 4.3;
     28 4.3 25 25 11 23 17 27 19 26 2.4 22 21 25 21 23 9.8 27 11 10 7.6 27 17 20 9.9 0 19;
     9.8 17 4.9 5 16 5 13 9.2 0.5 8.1 19 35 34 6.3 1.8 4.4 14 9.5 14 23 12 7 17 1.2 4.3 19 0
   ];
populationSize = 10;  % Population size for GA
maxGenerations = 100

; % Number of generations for GA
initialTemperature = 100;  % Initial temperature for SA
coolingRate = 0.9;      % Cooling rate for SA
mutationRate = 2;    % Maximum number of iterations for SA

[bestTour, bestCost] = hybridTSP(distMatrix, populationSize, maxGenerations, initialTemperature, coolingRate, mutationRate);
disp('Best Tour:');

Best Tour:

disp(bestTour);

     1     8    10    22    14     4     3    16    25     7    20     2    26    11    17    19     5    23    21     9    27    24    15     6    18    12    13

disp(['Best Cost: ', num2str(bestCost)]);

Best Cost: 140.7

function [bestTour, bestCost] = hybridTSP(distMatrix, populationSize, maxGenerations, initialTemperature, coolingRate, mutationRate)
    numCities = size(distMatrix, 1);
    
    % Generate initial tour using the nearest neighbor heuristic
    initialTour = nearestNeighbor(distMatrix);
    currentTour = initialTour;
    bestTour = initialTour;
    currentCost = calculateTourCost(currentTour, distMatrix);
    bestCost = currentCost;
    
    for generation = 1:maxGenerations
        temperature = initialTemperature / (1 + coolingRate * generation);
        
        % Apply Genetic Algorithm
        population = generatePopulation(populationSize, numCities);
        for i = 1:populationSize
            population(i, :) = mutate(population(i, :), mutationRate);
        end
        [population, costs] = evaluatePopulation(population, distMatrix);
        
        % Apply Simulated Annealing
        for i = 1:populationSize
            newTour = simulatedAnnealing(population(i, :), currentCost, temperature, distMatrix);
            newCost = calculateTourCost(newTour, distMatrix);
            
            if newCost < currentCost || rand() < exp((currentCost - newCost) / temperature)
                currentTour = newTour;
                currentCost = newCost;
                
                if newCost < bestCost
                    bestTour = currentTour;
                    bestCost = newCost;
                end
            end
        end
    end
end

function tour = nearestNeighbor(distMatrix)
    numCities = size(distMatrix, 1);
    tour = zeros(1, numCities);
    unvisited = 1:numCities;
    currentCity = randi(numCities);
    tour(1) = currentCity;
    unvisited(unvisited == currentCity) = [];
    
    for i = 2:numCities
        nearestCity = unvisited(1);
        minDistance = distMatrix(currentCity, nearestCity);
        
        for j = 2:length(unvisited)
            city = unvisited(j);
            distance = distMatrix(currentCity, city);
            if distance < minDistance
                nearestCity = city;
                minDistance = distance;
            end
        end
        
        tour(i) = nearestCity;
        currentCity = nearestCity;
        unvisited(unvisited == nearestCity) = [];
    end
end

function cost = calculateTourCost(tour, distMatrix)
    numCities = length(tour);
    cost = 0;
    for i = 1:numCities-1
        cost = cost + distMatrix(tour(i), tour(i+1));
    end
    cost = cost + distMatrix(tour(end), tour(1)); % Return to the starting city
end

function population = generatePopulation(populationSize, numCities)
    population = zeros(populationSize, numCities);
    for i = 1:populationSize
        population(i, :) = randperm(numCities);
    end
end

function mutatedTour = mutate(tour, mutationRate)
    numCities = length(tour);
    if rand() < mutationRate
        % Swap two random cities
        r1 = randi(numCities);
        r2 = randi(numCities);
        temp = tour(r1);
        tour(r1) = tour(r2);
        tour(r2) = temp;
    end
    mutatedTour = tour;
end

function [population, costs] = evaluatePopulation(population, distMatrix)
    populationSize = size(population, 1);
    costs = zeros(1, populationSize);
    for i = 1:populationSize
        costs(i) = calculateTourCost(population(i, :), distMatrix);
    end
    [~, idx] = sort(costs);
    population = population(idx, :);
    costs = costs(idx);
end

function newTour = simulatedAnnealing(tour, currentCost, temperature, distMatrix)
    numCities = length(tour);
    newTour = tour;
    
    % Perform 2-opt swap
    i = randi(numCities);
    j = randi(numCities);
    if i > j
        temp = i;
        i = j;
        j = temp;
    end
    
    while i < j
        temp = newTour(i);
        newTour(i) = newTour(j);
        newTour(j) = temp;
        i = i + 1;
        j = j - 1;
    end
    
    newCost = calculateTourCost(newTour, distMatrix);
    
    if newCost < currentCost || rand() < exp((currentCost - newCost) / temperature)
        currentCost = newCost;
    else
        % Revert the change
        while i < j
            temp = newTour(i);
            newTour(i) = newTour(j);
            newTour(j) = temp;
            i = i + 1;
            j = j - 1;
        end
    end
end

Answer 1

[Torsten]

Is this an attempt to code the Travelling Salesman Problem ?

If yes: it doesn't matter where the tour starts and ends since it is circular.

%% Travelling Salesman Problem
% Adapted from the TSP Example, Matlab Optimization Toolbox (https://mathworks.com/help/optim/ug/travelling-salesman-problem.html)
% by Santhanakrishnan Narayanan ([e-mail])

% Use this code to solve both symmetrical and asymmetrical TSPs based on binary integer programming.
% Required inputs: Distance matrix file
% Place the file in the same folder as the script
% Enter the file name along with extension like .csv/.xls
% The matrix file should be a square matrix
% Distance between (i,i) should be zero. Also the values for non-existing routes should be zero.
% Copyright 2014 The MathWorks, Inc.

%reading distanceMatrix from csv file
%distanceMatrixFile = input('Enter name of the file containing distance matrix (include extension like .csv/.xls): ','s');
%distanceMatrix = readmatrix(distanceMatrixFile);
%distanceMatrix = distanceMatrix(:,2:end);
distanceMatrix = [0 25 4.2 4 23 6 23 1.4 9.5 2.3 27 26 28 2.5 8 6.3 22 5.4 22 31 20 1.8 25 7.8 5.7 28 9.8;
     25 0 24 23 12 21 13 25 17 24 6.8 18 16 23 18 20 10 25 11 5.7 5.4 25 18 18 7.6 4.3 17;
     4.2 24 0 0.2 20 4.7 17 3.1 4.8 3.9 25 30 39 2 4.3 2.6 20 8.4 20 29 18 2.7 23 4 3.1 25 4.9;
     4 23 0.2 0 20 4.8 17 2.7 4.9 3.5 25 30 39 1.6 4.4 2.8 20 7.7 20 29 18 2.3 23 4.1 3.2 25 5;
     23 12 20 20 0 20 24 24 16 23 9.1 30 28 23 17 21 1.9 24 1.7 17 7.7 25 6 16 9.1 11 16;
     6 21 4.7 4.8 20 0 18 5.1 5.2 3.9 22 38 37 3.4 3.2 1.6 17 5 21 26 15 5.5 21 3.7 1 23 5;
     23 13 17 17 24 18 0 22 12 21 19 18 17 19 14 17 23 22 23 7 18 22 30 14 0.4 17 13;
     1.4 25 3.1 2.7 24 5.1 22 0 9.7 0.5 26 27 29 1.5 7 5.4 21 4.7 25 30 19 0.8 24 7.6 4.8 27 9.2;
     9.5 17 4.8 4.9 16 5.2 12 9.7 0 8.2 19 35 33 6.4 1.9 4.6 14 9.5 14 23 12 9.7 17 1.1 3.8 19 0.5;
     2.3 24 3.9 3.5 23 3.9 21 0.5 8.2 0 25 27 40 2 6.3 4.6 20 4.7 25 29 18 1.3 24 6.8 4.1 26 8.1;
     27 6.8 25 25 9.1 22 19 26 19 25 0 25 23 25 20 22 7.5 27 7.9 12 10 27 15 19 12 2.4 19;
     26 18 30 30 30 38 18 27 35 27 25 0 2 28 36 38 28 29 28 12 23 27 36 36 18 22 35;
     28 16 39 39 28 37 17 29 33 40 23 2 0 30 35 37 27 31 27 11 22 29 35 35 17 21 34;
     2.5 23 2 1.6 23 3.4 19 1.5 6.4 2 25 28 30 0 5.8 4 20 6.1 24 29 18 0.7 23 5.3 3.5 25 6.3 ;
     8 18 4.3 4.4 17 3.2 14  7 1.9 6.3 20 36 35 5.8 0 2.7 15 7.4 16 24 13 7.8 18 0.9 2.3 21 1.8;
     6.3 20 2.6 2.8 21 1.6 17 5.4 4.6 4.6 22 38 37 4 2.7 0 17 5.7 18 26 15 4.7 20 3.5 0.4 23 4.4;
     22 10 20 20 1.9 17 23 21 14 20 7.5 28 27 20 15 17 0 31 1.2 16 6.1 22 7.8 15 7.6 9.8 14;
     5.4 25 8.4 7.7 24 5 22 4.7 9.5 4.7 27 29 31 6.1 7.4 5.7 31 0 24 31 20 5.4 29 7.9 1.7 27 9.5;
     22 11 20 20 1.7 21 23 25 14 25 7.9 28 27 24 16 18 1.2 24 0 16 6.5  25 7.3 15 7.9 11 14;
     31 5.7 29 29 17 26 7 30 23 29 12 12 11 29 24 26 16 31 16 0 11 31 23 23 7 10 23 ;
     20 5.4 18 18 7.7 15 18 19 12 18 10 23 22 18 13 15 6.1 20 6.5 11 0 20 14 13 2.3 7.6 12;
     1.8 25 2.7 2.3 25 5.5 22 0.8 9.7 1.3 27 27 29 0.7 7.8 4.7 22 5.4 25 31 20 0 25 6 5.6 27 7;
     25 18 23 23 6 21 30 24 17 24 15 36 35 23 18 20 7.8 29 7.3 23 14 25 0 18 15 17 17;
     7.8 18 4 4.1 16 3.7 14 7.6 1.1 6.8 19 36 35 5.3 0.9 3.5 15 7.9 15 23 13 6 18 0 2.8 20 1.2;
     5.7 7.6 3.1 3.2 9.1 1 0.4 4.8 3.8 4.1 12 18 17 3.5 2.3 0.4 7.6 1.7 7.9 7 2.3 5.6 15 2.8 0 9.9 4.3;
     28 4.3 25 25 11 23 17 27 19 26 2.4 22 21 25 21 23 9.8 27 11 10 7.6 27 17 20 9.9 0 19;
     9.8 17 4.9 5 16 5 13 9.2 0.5 8.1 19 35 34 6.3 1.8 4.4 14 9.5 14 23 12 7 17 1.2 4.3 19 0
   ];
%creating city pairs and converting distance square matrix to distance
%column vector
fprintf('Creating city pairs\n');

Creating city pairs

numberOfCities = size(distanceMatrix,1); %number of cities
c=1;
for count = 1:numberOfCities:(numberOfCities*numberOfCities)
    cityPairs(count:numberOfCities*c, 1) = c;
    cityPairs(count:numberOfCities*c, 2) = 1:numberOfCities;
    distanceVector(count:numberOfCities*c, 1) = distanceMatrix(c,:)';
    c=c+1;
end
lengthDistanceVector = length(distanceVector);



%% Equality Constraints
fprintf('Creating equality constraints\n');

Creating equality constraints

%Number of trips = number of cityPairs
Aeq = spones(1:length(cityPairs));
beq = numberOfCities;

%Number of trips to a city = 1 and from a city = 1
Aeq = [Aeq;spalloc(2*numberOfCities,length(cityPairs),2*numberOfCities*(numberOfCities+numberOfCities-1))]; %allocate a sparse matrix to preallocate memory for the equality constraints;
c=1;
for count = 1:2:((2*numberOfCities)-1)
    columnSum = sparse(cityPairs(:,2)==c);
    Aeq(count+1,:) = columnSum'; % include in the constraint matrix
    rowSum = cityPairs(:,1)==c;
    Aeq(count+2,:) = rowSum';
    c=c+1;
end
beq = [beq; ones(2*numberOfCities,1)];

%Non-existing routes
nonExists = sparse(distanceVector == 0);
Aeq(2*c,:) = nonExists';
beq = [beq; 0];


%% Binary Bounds
%Setting the decision variables as binary variables
intcon = 1:lengthDistanceVector;
lb = zeros(lengthDistanceVector,1);
ub = ones(lengthDistanceVector,1);

%% Optimize Using intlinprog
fprintf('Solving the problem\n');

Solving the problem

opts = optimoptions('intlinprog','CutGeneration','Advanced','NodeSelection','mininfeas','Display','off');
[decisionVariables,optimumCost,exitflag,output] = intlinprog(distanceVector,intcon,[],[],Aeq,beq,lb,ub,opts);

%% Subtour Detection
tours = detectSubtours(decisionVariables,cityPairs);
numberOfTours = length(tours);
fprintf('Number of subtours: %d\n',numberOfTours);

Number of subtours: 13

%% Subtour Constraints
A = spalloc(0,lengthDistanceVector,0); % creating sparse inequality constraint matrix
b = [];
while numberOfTours > 1 % repeat until there is just one subtour
    b = [b;zeros(numberOfTours,1)]; % entering inequality constraints RHS
    A = [A;spalloc(numberOfTours,lengthDistanceVector,numberOfCities)]; % entering inequality constraints LHS
    for count = 1:numberOfTours
        inequalityConstraintNumber = size(A,1)+1;
        subTourId = tours{count}; % Extracting subtour one by one

        % adding subtour constraints (inequality constraints)
        subTourPairs = nchoosek(1:length(subTourId),2);
        for jj = 1:size(subTourPairs,1) % Finding variables associated with the current sub tour
            subTourVariable = (sum(cityPairs==subTourId(subTourPairs(jj,1)),2)) & ...
                       (sum(cityPairs==subTourId(subTourPairs(jj,2)),2));
            A(inequalityConstraintNumber,subTourVariable) = 1;
        end
        b(inequalityConstraintNumber) = length(subTourId)-1; % reducing number of trips allowed by One Ex., A-B-A: 2 -> 1
    end

    % Optimize again
    fprintf('\nsolving the problem again eliminating subtours\n');
    [decisionVariables,optimumCost,exitflag,output] = intlinprog(distanceVector,intcon,A,b,Aeq,beq,lb,ub,opts);

    % Check for subtours again
    fprintf('Checking again for subtours\n');
    tours = detectSubtours(decisionVariables,cityPairs);
    numberOfTours = length(tours);
    fprintf('Number of subtours: %d\n',numberOfTours);
end

solving the problem again eliminating subtours

Checking again for subtours

Number of subtours: 7

solving the problem again eliminating subtours

Checking again for subtours

Number of subtours: 7

solving the problem again eliminating subtours

Checking again for subtours

Number of subtours: 2

solving the problem again eliminating subtours

Checking again for subtours

Number of subtours: 3

solving the problem again eliminating subtours

Checking again for subtours

Number of subtours: 1

%% Solution Quality
%smaller the value better the solution
fprintf('\nSolution Quality: %f (lesser the better)\n',output.absolutegap);

Solution Quality: 0.000000 (lesser the better)

fprintf('Optimized tour route:');

Optimized tour route:

celldisp(tours);

 
tours{1} =
 
     8    10    18    25     7    13    12    20     2    26    11    19    23     5    17    21    27     9    24    15     6    16     3     4    14    22     1

 

fprintf('Note: The numbers correspond to order of cities in the input file\n');

Note: The numbers correspond to order of cities in the input file

fprintf('Total distance of the optimal route: %d\n', optimumCost);

Total distance of the optimal route: 1.075000e+02

function subTours = detectSubtours(x,idxs)
% Returns a cell array of subtours. The first subtour is the first row of x, etc.
x(x<0.0001)=0;
r = find(x); % indices of the trips that exist in the solution
substuff = idxs(r,:); % the collection of node pairs in the solution
unvisitedSubToursSubTours = ones(length(r),1); % keep track of places not yet visitedSubTours
curr = 1; % subtour we are evaluating
startour = find(unvisitedSubToursSubTours,1); % first unvisitedSubToursSubTours trip
    while ~isempty(startour)
        home = substuff(startour,1); % starting point of subtour
        nextpt = substuff(startour,2); % next point of tour
        visitedSubTours = nextpt; unvisitedSubToursSubTours(startour) = 0; % update unvisitedSubToursSubTours points
        while nextpt ~= home
            % Find the other trips that starts at nextpt
            [srow,scol] = find(substuff == nextpt);
            % Find just the new trip
            trow = srow(srow ~= startour);
            scol = 3-scol(trow == srow); % turn 1 into 2 and 2 into 1
            startour = trow; % the new place on the subtour
            nextpt = substuff(startour,scol); % the point not where we came from
            visitedSubTours = [visitedSubTours,nextpt]; % update nodes on the subtour
            unvisitedSubToursSubTours(startour) = 0; % update unvisitedSubToursSubTours
        end
        subTours{curr} = visitedSubTours; % store in cell array
        curr = curr + 1; % next subtour
        startour = find(unvisitedSubToursSubTours,1); % first unvisitedSubToursSubTours trip
    end
end

[Amna Habib]

I appreciate your effort dear @Torsten. Yes, you are right. It is TSP but I have to sove a case study in which node 1 is distribution center. I have to start and end at Node 1.

Moreover, my algorithm is a hybrid algorithm. I need to apply this hybrid algorithm rather than any other TSP algorithm

[Torsten]

But how can a node like node 1 in your case have a special meaning for the solution of the TSP ? The solution is a circular path - so each node can be considered to be the "distribution center". Or is it a speciality of the solution algorithm you use that it needs a certain node as anchor point ?

[Amna Habib]

Yes, node 1 has special meaning in this case study. You are right, TSP is circular but i need this special one.

I have made this hybrid algorithm. But now I have to solve a case study in which Node 1 is origin. We have to start and end at Node 1.

[Torsten]

As already said, a solution of a TSP starts and ends in every node.

So if you have a tour

5 3 1 2 4 5

starting and ending in 5, you can equivalently write it as

1 2 4 5 3 1

2 4 5 3 1 2

3 1 2 4 5 3

4 5 3 1 2 4

for tours starting and ending in 1, 2, 3 or 4.

[Amna Habib]

Got it! Thank you so much @Torsten

I was missing this point. Your contribution is highly appreciated.

[Amna Habib]

Dear @Torsten, please see my code, the resultant output do not start and end on th same node. While this should happen. Can you please mark the change, I have to make???

[Torsten]

It's not usual for the output that a node is repeated because it's implicitly assumed that the circle is closed by connecting the last node to the first node.

If you wish, you can add the first node to the end of the array "bestTour" before printing it :

[bestTour, bestCost] = hybridTSP(distMatrix, populationSize, maxGenerations, initialTemperature, coolingRate, mutationRate);
bestTour(end+1) = bestTour(1);

But "BestCost" must be 107.5 - you get something around 140. Does your method only give an approximate solution or should it end up optimal ?

[Amna Habib]

Thanks a lot @Torsten! I have now obtained the required result.

My method gives optimal soluation. I could not understand the meanig 'But "BestCost" must be 107.5 - you get something around 140.'

I think the best cost should be 140.1.

[Torsten]

Did you look at the result from "my" code and the cost for the tour

8 10 18 25 7 13 12 20 2 26 11 19 23 5 17 21 27 9 24 15 6 16 3 4 14 22 1 8

?

The cost of this tour is 107.5 .

[Amna Habib]

Thanks @Torsten. Got it!

Your expertise in this fielf is really appreciateable!

Answer 2

[John D'Errico]

Your problem is, you want to force the code to return a tour that starts at point 1, and ends there. And that is where you are WRONG.

All you need to do is find the order of points 2 through 27 on that tour. You already know the start and end points are point number 1.

[Amna Habib]

Thanks @John D'Errico for your time. The order of points 2 through 27 will not enough for this purpose.