Dynamic Programming and Travelling Salesman Problem. C Code and Octave Script

Dynamic Programming and Travelling Salesman Problem. C Code and Octave Script

.

1. Abbreviation

Permutations of a set of n elements:

Note:

• 0! = 1! = 1

A combination is a selection of k elements from a set containing n elements:

Notes:

• The (2) are also called binomial coefficients

• See also Pascal’s triangle

.

2. Introduction

Dynamic programming was developed from Richard Bellman (see the [1], [2]) in the fifties of the twentieth century and occupies an important place for solving optimization problems, as well as problems of the calculus of variations.

Dynamic programming is the search for the optimal solution to a given problem, which can be divided into many steps and the optimal solution can be found only for the current step, which greatly simplifies calculations. Sometimes we can talk about a purposeful search of possible options with an assessment of their optimality, while the complexity of the task decreases compared to a complete search of all options.

Suppose there is a process of n steps. At each step a benefit function (or cost function) is set:

The calculation is performed using the formula:

The function f(x) is called the Bellman function.

In the first step the Bellman function f(x) is calculated:

The maximum in (3) is taken if d() is a benefit function and it must be maximized, and the minimum if d( ) is a cost function and it must be minimized. The function f(x) determines the optimal solution for the current step. Thus the dynamic programming method reduces the necessary search for options by dividing the task into steps. See below.

.

3. Travelling Salesman Problem

The traveling salesman problem is formulated as follows. There are n cities. When leaving one of the cities, the traveling salesman must visit each city once, and then return to the first city back. At the same time the total distance traveled should be minimal. If the starting city is fixed, then a complete search of all the options is (n-1)!. If the choice of the initial city is arbitrary, then a complete search of all options is n!

Note: The traveling salesman problem has already been discussed in the article: “Greedy Algorithm and Travelling Salesman Problem Using Octave GNU Tool”

Mathematically the problem can be described as follows. Let’s take for example 4 cities, the distance between which is shown in Figure 3-1:

.

Figure 3-1: Four Cities and the Distances between them

.

The distance between cities can be set as a matrix:

For example, the matrix will look like in Figure 3-1:

The matrix is symmetric and there are zeros along the diagonal since:

There are two methods for solving the traveling salesman problem using the dynamic method. The first method was proposed by M. Held and R. Karp, see the [3]. The second method was proposed by R. Bellman, see the [4], [5].

.

4. M. Held and R. Karp Method for Travelling Salesman Problem

There are n cities: 1, 2, 3, … , n. Let the traveling salesman start from city 1. Then we need to find the order of traversal of the remaining n-1 cities so that the distance traveled is minimal. We will consider various sets S made up of cities:

Let’s define the Bellman function f(S,last), which denotes the minimum path (or minimum cost) starting in city 1 and passing through all cities from the set S and ending in city last. Next we will denote the set S with the index k:

The number of sets S with index k is calculated using the formula:

The Bellman function is calculated as follows:

1. In the first step k=1, sets S with index k=1 are used. These sets contain only one city and the number of such sets is n-1:

   

   At the first stage the Bellman function is calculated using the formula:

   

2. Steps k = 2 … n-1. Sets S with index k are used here. These sets contain k cities and the number of such sets is calculated using formula (11). In this case the Bellman function has already been calculated for the previous step k-1.

   

   Next two cities are selected from this set, which we will call prev (previous) and last:

   

   Next the Bellman function is calculated for k-th step using the already calculated value of the Bellman function for (k-1)-th step using the formula:

   

   Using all combinations of prev and last from (14), the Bellman function is calculated for the k-th step. It is easy to understand that the number of combinations for (14) is calculated using the formula:

   

   And taking into account formula (11), we finally obtain:

   

   Then the calculations for all k = 2 … n-1:

   

   The proof of formula (18) is given in Appendix.

   For each set S with index k, it is necessary to calculate the Bellman function for each city from this set, taken as the last. This means that it will require k memory locations. Taking into account the number of such sets S with index k, we get the required memory size:

   

   Then the memory size for k = 1 … n-1 is calculated using the formula:

   

   This is the memory required for items 1. and 2. he proof of formula (20) is also given in Appendix.

3. In the last step, the traveling salesman returns to city 1. The Bellman function has already been calculated for the set S with index n-1:

   

   The optimal value is determined by the formula:

   

   for all possible combinations of (n-1) values of the last = 2, 3, … , n

   Next using the optimal last value and going from the end to the beginning, we find the optimal path. In this case an additional previous array is used. See below.

Let’s solve problem (6) using the described method. There are four cities 1, 2, 3, 4, and we will start from city 1:

1. Calculate the Bellman function in the first step k=1:

   

2. Now we calculate the steps k=2, 3:

   

   We get three combinations of two cities:

   

   The Bellman function by formula (16):

   

   

3. In the last step the penultimate city will be found of the optimal path:

   
   Optmal path 1:
   
   Optmal path 2:
   

Notes:

• In the C code (see the TravellingSalesmanProblemHeldKarpWithVsCopilot.c/h), two arrays are created that correspond to the above f(S,last) and previous(S,last) arrays:

  

  Cities are coded: 0, 1, 2, …, n-1, which corresponds to cities 1, 2, 3, …, n.

  The traveling salesman starts from city 1 and arrives at the end in the same city 1 (coded as 0).

  The first index [numSubsets] is encoded bitwise:

  bit0 – city 2
  bit1 – city 3
  …
  bit(n-2) – city n

  For example:

  

  The second index is [n-1]:

  0, 1, 2, …, n-2 which corresponds to cities 2, 3, 4 …, n

  For example:

  

• In the Octave code (see the TravelingSalesmanProblemHeldKap.m), two arrays are created that correspond to the above f(S,last) and previous(S,last):

  

  Cities are encoded in the usual way: 1, 2, 3, …, n, which corresponds to cities 1, 2, 3, …, n.

  The traveling salesman starts from city 1 and arrives at the end in the same city 1.

  The first index (city Subsets) is encoded bitwise:

  bit0 – city 2
  bit1 – city 3
  …
  bit(n-2) – city n

  For example:
  
  

  The second index (cities_nmb-1):

  1, 2, …, n-1 which corresponds to cities 2, 3, …, n

  For example:

  

• It was previously proved that the array for f (…) and previous(…) must have the size determined by formula (20). We use arrays twice as large. This is due to the fact that the software implementation in this case is much simpler, although not all cells of the array are used in this case

• If we repeat all these calculations for different starting cities, we get the following optimal solutions:

  
  Table 4-1: Optimal Tracks for the Matrix (6)
  

.

5. R. Bellman Method for Travelling Salesman Problem

Bellman, the author of dynamic programming, proposed another solution to the problem. See the [4], as well as [5]. There are n cities and the traveling salesman starts from city i and then passes through the remaining cities and then returns to the i city:

Bellman introduced the optimal path function for such a transition:

I will change the original notation of function (24) from [4] to another notation:

This notation (25) is more appropriate to the previous notation used in M. Held and R. Karp.

In this case the recurrence relation for finding the optimal function f at the k step is:

In this case the function f at the last (or first from below) step is calculated using the formula:

In this case, functions (25) and (26) are not equivalent, since in function (25) the start takes place in city i and ends in city i, and in (26) and (27) the start takes place in the city j and ends in the city i.

The Bellman method is shown for four cities in the Figure 5-1.

Let’s solve problem (6) using the described method. There are four cities 1, 2, 3, 4. We will start from the city i = 1. Thus the following designations can be adopted:

The Bellman algorithm is easiest to implement using recursive calls. The algorithm starts from the top down (See the Figure 5-1), goes down, and then rises back up with the result. Here in the description I start from the bottom up so as not to complicate the analysis of the solution. The solution uses the auxiliary array next, which is then used to determine the optimal path.

.

Figure 5-1: Bellman Method Four Cities

.

Step 1

Step 2

Step 3

Step 4

Optimal Path 1:

Optimal Path 2:

Notes:

• In the C code (see the TravellingSalesmanProblemBellmanWithVsCopilot.c/h), two arrays are created:

  

  Cities are coded: 0, 1, 2, …, n-1, which corresponds to cities 1, 2, 3, …, n.

  The traveling salesman starts from city i and returns to the same city i.

  The first index is [n]:

  0, 1, 2, …, n-1 which corresponds to cities 1, 2, 3, …, n

  The second index [numSubsets] is encoded bitwise:

  bit0 – city 1
  bit1 – city 2
  …
  bit(n-1) – city n

  Note. This encoding is redundant compared to the encoding in the M. Held and R. Karp method. A more optimal encoding can be successfully used in the Bellman method too, which will lead to a reduction in working arrays.

  For example:

  

• In the software Octave code (see the TravelingSalesmanProblemBellman.m), two arrays are created:

  

  Cities are encoded in the usual way: 1, 2, 3, …, n, which corresponds to cities 1, 2, 3, …, n.

  The traveling salesman starts from city i and arrives at the end in city i.

  The first index (cities_nmb):

  1, 2, …, n which corresponds to cities 1, 2, …, n

  The second index (city Subsets) is encoded bitwise:

  bit0 – city 1
  bit1 – city 2
  …
  bit(n-1) – city n

  For example:

  

• If we repeat all these calculations for different starting cities, we get the solution already given in Table 4-1

.

6. Notes and Conclusions

1. The computational complexity of the Bellman method (see the [5]) is approximately the same as in the method of M. Held and R. Karp. But the implementation of the Bellman method uses recursion, which requires additional RAM on the stack and control over its possible overflow.

   There are applications where the use of recursion is undesirable due to high software reliability requirements.

   Therefore I prefer a software implementation by M. Held and R. Karp

2. In the traveling salesman problem (6) in the previous section a symmetric matrix was used. The elements of this matrix are the distance between cities i and j. It is clear that the distance between cities i and j is equal to the distance between cities j and i. The result is a symmetric matrix. Another assignment of this matrix is possible. For example, we will set the fuel consumption required to get from the city i to j. In this case the fuel consumption for the trip from the city of j in i will be different. This is because the road from the city of i in j, for example, goes uphill, and when traveling back from j in i, the road will go downhill and fuel consumption will be significantly lower. This will result in an unbalanced matrix. Next I will use an asymmetric matrix for the traveling salesman problem from [6]:

   

   Solving this problem, using the methods described above, we get:

   

   Table 6-1: Optimal Tracks for the Matrix (28)

   .

3. A direct search of all the options in the traveling salesman problem gives the complexity of the algorithm:

   

   It follows from (18) that the complexity of the algorithm by dynamic programming is:

   

   Dynamic programming provides a simpler option, but we still have exponential complexity, which corresponds to a Horrible algorithm. Thus the solution of the traveling salesman problem works only for small values of n.

   Approximate methods of solving the problem are used for large n

.

7. C Language TravellingSalesmanProblemHeldKarpWithVsCopilot.c/h

The code was generated by artificial intelligence (AI): Copilot, Claude 3.5, Sonnet model. Then the bugs were fixed and the code was optimized. The software solves the traveling salesman problem for a fixed starting city of 1.

/*
The number of units in a word
Input
uint32 n – word
Return
uint8 – the number of units in a word
*/
uint8 countBits(uint32 n)

/*
Solve the Travelling Salesman Problem (TSP)
using Dynamic Programming
Method of Michael Held and Richard Karp
Cities coding in C: 0, 1, 2, …, n-1 =>
=> according to 1, 2, …, n cities
Start and end in city 0
Input
float** dist – distance matrix
uint8 n – number of cities
Return
TSPResult_t – optimal path(s), number of solutions, and min. distance
*/
TSPResult_t solveTSP(float** dist, uint8 n)

.

8. Octave GNU file TravellingSalesmanProblemHeldKarp.m

The script solves the traveling salesman problem for a fixed and arbitrary starting city:

% Number of Ones in the Word
% Input:
% word
% Output parameters:
% ones_nmb – ones number in word
function ones_nmb = OnesNumberInWord(word)

% Travelling Salesman Problem, all optimal paths with Fixed Start Point
% Input:
% distanceMatrix – distance matrix between points
% startPoint – start point
% citiesBoundary – max. cities number
% Output parameters:
% all_tracks – all optimal point tracks
% minDistance – optimal distance
function [all_tracks minDistance] = TravellingSalesmanPrblmWithFixedStart(distanceMatrix, startPoint, citiesBoundary)

% Travelling Salesman Problem, all optimal paths
% Input:
% distanceMatrix – distance matrix between points
% citiesBoundary – max. cities number
% Output parameters:
% track – all optimal point tracks
% minDistance – optimal distance
function [all_tracks minDistance] = TravellingSalesmanPrblm(distanceMatrix,citiesBoundary)

.

9. C Language TravellingSalesmanProblemBellmanWithVsCopilot.c/h

The code was generated by artificial intelligence (AI): Copilot, Claude 3.5, Sonnet model. Then the bugs were fixed and the code was optimized. The software solves the traveling salesman problem for a fixed starting city

/*
To create a new set of cities
Input
void
Return
CitySet_t – the new set
*/
CitySet_t createSet(void)

/*
Add city to the set
Input
CitySet_t* set – pointer to the set
uint8 city – the city to add
Return
CitySet_t – the new set
*/
void addCity(CitySet_t* set, uint8 city)

/*
Convert set to BitSubset for memoization
Input
CitySet_t set – the set
Return
uint32 – subset
*/
uint32 setToBitSubset(CitySet_t set)

/*
Create set from BitSubset
Input
uint32 subset
uint8 n – number of cities
Return
CitySet_t – the set
*/
CitySet_t setFromBitSubset(uint32 subset, uint8 n)

/*
Remove city from the set
Input
CitySet_t set
uint8 city
Return
CitySet_t – the set
*/
CitySet_t removeCity(CitySet_t set, uint8 city)

/*
Bellman’s recursive function funcBellman(i, {j1,j2,…,jk}, startCity)
Input
uint8 i
CitySet_t remainingCities
uint8 startCity
Return
float – the cost of the optimal path starting from city i
and visiting all cities in remainingCities, returning to startCity
*/
float funcBellman(uint8 i, CitySet_t remainingCities, uint8 startCity)

/*
Efficient path reconstruction using next array
Input
uint8 i
CitySet_t remainingCities
uint8 start
uint8 n
TSPResult_t* result
Return
void
*/
void reconstructPath(uint8 start, uint8 n, TSPResult_t* result)

/*
Main function to solve TSP
Input
uint8 n
CitySet_t remainingCities
uint8 start
uint8 startCity
Return
TSPResult_t – the result containing the optimal path and cost
*/

TSPResult_t solveTSP(uint8 n, uint8 startCity)

.

10. Octave GNU file TravellingSalesmanProblemBellman.m

The script solves the traveling salesman problem for a fixed and arbitrary starting city:

% Convert set to BitSubset for memoization
% Input:
% remainingCities – structure of the remaining cities
% Output parameters:
% subset
function subset = setToBitSubset(remainingCities)

% Remove city from the set
% Input:
% prev_set – structure of the cities
% city – for remove
% Output parameters:
% new_set – structure of the cities without city
function new_set = removeCity(prev_set, city)

% Recursive Bellman function
% Input:
% i – start city for the remaining cities
% remainingCities – structure of the remaining cities
% startPoint – start city for all task
% distanceMatrix – distance matrix between points
% f_all_subsets – cost array
% next_all_subsets – next cities array
% solutionCounter – tracks support
% secondCity – tracks support
% Output parameters:
% cost – optimal distance
% f_all_subsets – cost array after update
% next_all_subsets – next cities array after update
% solutionCounter – tracks support after update
% secondCity – tracks support after update
function [cost, f_all_subsets, next_all_subsets, solutionCounter, secondCity] = funcBellman(i, remainingCities, startPoint, distanceMatrix, f_all_subsets, next_all_subsets, solutionCounter, secondCity)

% Travelling Salesman Problem, all optimal paths with Fixed Start Point
% Input:
% distanceMatrix – distance matrix between points
% startPoint – start point
% citiesBoundary – max. cities number
% Output parameters:
% all_tracks – all optimal point tracks
% minDistance – optimal distance
function [all_tracks minDistance] = TravellingSalesmanPrblmWithFixedStart(distanceMatrix, startPoint, citiesBoundary)

% Travelling Salesman Problem, all optimal paths
% Input:
% distanceMatrix – distance matrix between points
% citiesBoundary – max. cities number
% Output parameters:
% track – all optimal point tracks
% minDistance – optimal distance
function [all_tracks minDistance] = TravellingSalesmanPrblm(distanceMatrix,citiesBoundary)

.

11. Download the TravellingSalesmanProblem*.m/c/h

You can download the files:

TravellingSalesmanProblemHeldKarpWithVsCopilot.c/h

TravellingSalesmanProblemHeldKarp.m

TravellingSalesmanProblemBellmanWithVsCopilot.c/h

TravellingSalesmanProblemBellman.m

with the button:

.

12. Literature / References

[1] R. Bellman “Dynamic Programming“, Princeton University Press, New Jercy, 1957

[2] R. Bellman, S. Dreyfus “Applied Dynamic Programming“, Princeton University Press, New Jercy, 1962

[3] M. Held, R. Karp “A Dynamic Programming Approach to Sequencing Problems“, Journal of the Society for Industrial and Applied Mathematics (SIAM), Volume 10, No. 1, March 1962

[4] R. Bellman, “Dynamic Programming Treatment of the Travelling Salesman Problem“, Journal of the Association for Computing Machinery (ACM), Volume 9, Issue 1, 1962

[5] S.A. Kanzedal, M.V. Kostikova “Dynamic programming for the traveling salesman problem“, Magazine: “Automated control systems and automation devices“, Issue 166, Kharkiv National University of Radio Electronics, 2014

[6] V. I. Mudrov “The Traveling salesman problem“, Publishing House “Znanie“, Moscow, 1969

 

.

Appendix. Proof of several formulas

Let’s prove formulas (A-1) and (A-2):

Consider the function:

Expand the brackets in (A-3):

Let’s take the derivative of x of the function (A-4):

Let’s differentiate (A-5) again:

Substituting the value x=1 in (A-5) and (A-6), we prove the identities (A-1) and (A-2)

.