Combinatorial optimization by metaheuristic
Workshop

Tutor guide

Version	Date	Designer	Proofreader	English translator	Comments
2.1	03/02/2023	Benjamin COHEN BOULAKIA			Translation of the code into English, in view of the translation of the entire Notebook
2.1	21/09/2022	Benjamin COHEN BOULAKIA			Reorganisation of the Workshop as discussion points to be addressed with students A few more detailed explanations
2.0	15/03/2022	Benjamin COHEN BOULAKIA			Code optimization : Replacement of lists by more suitable structures for solutions (tuples), and the tabu list (dictionaries) Memoization of functions (neighbourhood, cost and value of a solution) Added questions on implementation optimization Reorganization of the sequence at the level of the basic functions A few more detailed explanations Addition of theoretical and technical questions Relocation of the Quality of the solution part to the next Workshop, dedicated to the experimental study and optimization of the metaheuristic.
1.2	22/06/2020	Benjamin COHEN BOULAKIA	Pierre HALFTERMEYER		The Quality part of the solution introduces the lower bound principle in more detail. Some text edits.
1.1	27/05/2019	Benjamin COHEN BOULAKIA			Removal of the presentation of KSP, incorporated into the previous Workshop.
1.0	03/05/2019	Benjamin COHEN BOULAKIA			Original version.

Workshop skills

• Testing one or more scientific computing libraries
• Producing optimized code generating the neighbourhood of a solution for a problem
• Developing a metaheuristic by exploiting a neighbourhood generation method

A Python script including the final version of the multi-start tabu of this Workshop is available, it’s the file Multistart Tabu Knapsack.py. It is possible to extract the simple version of the tabu search (without multistart) from it.

Goals

This Workshop addresses a problem known as $\mathcal{NP}$-Hard, the one-dimensional knapsack problem in 0-1. Students will implement a tabu search method using Python 3, determining the most suitable data structures and processing tasks in terms of runtime performance. Students must edit this Notebook by answering the questions asked, and by completing the code provided. They can also then export the Python code for execution outside of Jupyter.

It is fundamental that the students manage to implement an optimized version of their metaheuristic. On the one hand by using the right data structures (list, tuple, dictionary), on the other hand by optimizing processing (slicing, Comprehension, memoization). The following Workshop will look at the statistical behaviour of this algorithm, and unoptimized code will not allow them to obtain satisfactory results. The most delayed students can skip the multi-start version of the heuristic.

Philosophy and recommendations

The goal here is not simply to code an algorithm, but also to raise awareness of the multiple implementation possibilities, and their respective performance, as well as the multiple variants of any metaheuristic (not only the tabu method). It is therefore important not to ignore the questions posed in the text, and to lead to a debate on these questions.

In addition to the answers to the questions, elements of explanation are provided. It is imperative to present these explanations to the students, and to animate if necessary a discussion (or even some questions / answers) around these points.

In addition, many links are scattered throughout the Workshop, students should be encouraged to visit them systematically, especially if they are unfamiliar with the subject concerned. If necessary, they should be encouraged to navigate the technical documentation as quickly as possible without getting lost in it.

We meet up with our student, whose Smartphone is bursting at the seams. He has not given up taking his music library with him, but his attempt to rely on the knapsack problem was not very conclusive. The approach by linear programming did not work, due to the integer variables, which make its problem-solving by the simplex impossible. It's a shame, because this algorithmic problem corresponds to many challenges, especially from logistics. Which ones can you imagine ?

Solution In this case, we consider the problem of the one-dimensional knapsack (the objects have a one-dimensional size, and the bag is therefore defined by a single capacity). This variant is the simplest, the forms of the problem that are most encountered in the industry are more complex. More specifically, there are versions with dimensions and a given weight (weights are simply added up, but it’s more complicated for dimensions, because there may be empty spaces depending on how the objects are placed). We can even consider several containers (multiple knapsack, or MKP). All these variants have many applications :

• the most literal application : the loading of cargo (planes, trucks, boats…)
• tank storage of liquids that cannot be mixed (typical case of MKP)
• financial systems, where the idea is : given a certain amount of investment in projects, which projects to choose so that the whole thing brings in the most money possible
• the cutting of materials, in order to minimize the losses due to scraps
• etc.

The simple version considered here still has some practical applications. The best known one, which students are unlikely to imagine if they have not come across it in their research work, is in cryptography, where it is at the origin of the first asymmetric encryption algorithm in 1976, designed at MIT (and which has since been broken).

1 Basic modelling and implementation

We will start by setting up a simple generic approach to solve this knapsack problem.

Explanations Since the problem is $\mathcal{NP}$-Hard, it’s only natural that it lean towards a heuristic approach. We could build an ad hoc heuristic (designed specifically for this problem), but we will favour a more generic approach here. We will first implement a discrete local search.

1.1 Heuristic modelling

The Hill Climbing algorithm (that we should sometimes rather call Valley Descending) is a fairly simple solution to implement. What is this algorithm ?

Solution

1. We consider an initial solution.
2. A loop generates all the neighbours of the current solution, and chooses the best one as the new current solution.
3. We leave the loop when it no longer allows us to find a better solution than the current solution.

This algorithm uses the same notion of neighbourhood as some metaheuristics, which will therefore have to be defined first. One can wonder if it is a heuristic, or a metaheuristic. It's a question of vocabulary:

• On the one hand, we can consider that not, because we do not find the step of diversification typical of metaheuristics, which sometimes makes it possible to extract oneself from a local optimum. Here, it is strictly a greedy algorithm, which performs a local optimization, a bit like a gradient descent, but on a discrete structure.
• On the other hand, we can consider that yes, since it is a local search algorithm, exploiting the notion of neighbourhood. It is therefore a generic method (hence the term meta) that can be used to find an approximate solution to a difficult problem (hence the term heuristic). At the extreme, it can be considered as a metaheuristic doing only intensification, with zero diversification.

The literature tends to consider that the answer is 'no', so we stick to this vision.

Moreover, students have very recently seen an optimization algorithm that works according to a very similar principle, and which, apart from that, generates the optimal solution (i.e. it addresses problems using a notion of neighbourhood in which a local optimum is also a global optimum). This is the Simplex. In the case of linear programs in real numbers (but not in integers), we can see the simplex as an algorithm that considers only a subset of feasible solutions, those corresponding to a vertex of the polytope. The neighbourhood of solution corresponds to the fact that two vertices are linked by an edge on the polytope. The simplex starts from any vertex of the polytope (often the origin, of zero variables), and chooses, at each iteration, the neighbouring vertex that minimizes the fitness function, until reaching the optimum. Continuous linear programs are therefore a category of problems in which a local optimum is also a global optimum.

A small detail : in fact, the Simplex has a complexity that is exponential in the worst case (but in most cases, we remain polynomial).

The first (and main) step in modelling a problem for a metaheuristic is therefore the definition of the neighbourhood structure (at least in the case of a trajectory method). What simple modification of a solution can we consider in our case of knapsack problem ?

Solution As a general principle, we consider a solution close to the current solution as a valid solution which is equal to it within one atomic modification.

For our case, two solutions are neighbours if we can go from one to the other by adding or removing an object from the bag.

This means that the neighbourhood of a solution is the set of neighbouring solutions of this solution. Using this definition of the neighbourhood of a solution, it will be possible to implement various metaheuristics to deal with this problem.

What’s left to address is the computational representation of our data. Knowing that a solution represents, for each object, the fact that we take this object or not, how will we computationally model such a solution ?

Solution Be careful, we do not yet ask the question of Python structures, just their numerical representation).

A solution is a word of $n$ bits with $n$, the number of available objects. The value of the $i^{th}$ Boolean states whether the $i^{th}$ object is present in the solution.

Moreover, once we consider this representation, we can re-express the notion of neighbourhood according to this representation : The neighbours of a solution are all binary words of $n$ bits having a Hamming distance of 1 with the solution.

From there, representing the weights and values of each object is quite natural. What structures will be adapted ?

Solution Two arrays of $n$ values, with $n$ the number of available objects, one array for the weights, and one for the values.

1.2 Basic structures and functions

Explanations We determined above the data structures that we will use, we can now implement, write the associated basic functions (calculation of the value of a solution, of its weight), then implement a function to generate the neighbours of a solution. Once these tools are created, we can implement several metaheuristics, test them, compare them etc.

You already know the very general principle of the algorithms that we are going to implement. It's about generating valid solutions, assessing their value, and keeping the best one. Which Python data structures are we going to use to prioritize runtime performance ? Propose a structure for :

• an instance of the problem
• a solution of one instance

Solution

• Instance : An instance consists of :
  • the number of objects
  • the bag’s capacity
  • all objects, each defined by its weight and value
  The first two do not pose a problem, since they are integers. From the point of view of code performance, it is especially on the last two that you have to choose the right solution. In this case, dictionaries are appropriate. In fact, when calculating the weight or the value of a solution, we will not access all the objects, but only those present in the solution. And there, the ability that dictionaries have to access an element whose key is known in $O(1)$ (i.e. in constant time) thanks to a hashing method will be a valuable asset.
  
  
• Solution : an instance is represented by a binary word of size of the number of objects.
  The most suitable structure is a tuple of Booleans. In our case, we will have to intensively generate many solutions (whether one after the other in a trajectory approach, or at the same time in a population method), which, once created, will only be assessed (in terms of value and weight), and possibly retained. The performance of tuples are much more relevant to us than the use flexibility of lists, which we do not really need.

We are not going to complicate our lives, these structures will be stored in global variables. This is a bad habit for an application project, but here we are experimenting with algorithms, nothing will go into production.

Let's start with instances. If we want to test our functions, we must be able to randomly generate an instance (given a weight and a maximum value for any object). For the number of objects and the bag’s capacity, two integers will do :

nb_objects = 10
capacity = 20

With Python’s standard random library, our random instance generation function will be very simple to write. It will be able to return a tuple of dictionaries (weight, value), which it will have generated with Comprehension.

import random

def random_objects(max_weight, max_val):

 objects_weight =\
 {i: random.randint(1, max_weight) for i in range(nb_objects)} #SOLUTION
 objects_value =\
 {i: random.randint(1, max_val) for i in range(nb_objects)} #SOLUTION
    
 return objects_weight, objects_value

random.seed(a=3) # using a seed explicitly will reproduce the initial conditions
 # and compare the behaviour of different algorithms on the same instance

objects_weight, objects_value = random_objects(10, 10)
print(str(objects_weight) + "\n" + str(objects_value))

{0: 4, 1: 10, 2: 9, 3: 3, 4: 6, 5: 10, 6: 8, 7: 10, 8: 2, 9: 10}
{0: 1, 1: 8, 2: 5, 3: 9, 4: 4, 5: 4, 6: 8, 7: 9, 8: 9, 9: 8}

Now that we can generate a random instance, we will be able to implement and test the basic functions that any metaheuristic needs. Which ones ?

Solution We need two functions :

• A function to calculate the weight of the solution, i.e. the sum of the weights of the objects in the bag, so that the metaheuristic checks if a solution is valid (if its weight complies with the bag’s capacity) :
  $\displaystyle\sum_{i=1}^{n}{w_ix_i}$
  
  with (notations taken from the previous Workshop)
  • $w_k\in\mathbb{Z}^{⁺}$ the weight of object $k$
  • $x_k\in\{0,1\}$ the decision variable of object $k$
  
  
• A function calculating the value of a solution, so that the metaheuristic can compare two solutions (regardless of the metaheuristic considered). The principle is the same as for the previous function, but using the value of the object rather than its weight.

It's simple to implement, especially with the function sum. But rather than literally using the formula we just talked about, we’ll use Python’s capabilities : we’ll build a subset of the values of object_weight that contains only the values of the objects in the bag. And there is a way to produce these values, it is the Comprehension (whether on lists or dictionaries).

Therefore, what must be done is to successfully write a Comprehension that returns the indices of the objects present in the bag. For that, we'll add a conditional List Comprehension, which will enumerate the content of the bag, and only filter the True values. Let's try to write this Comprehension.

Explanations As mentioned above, students can propose to implement the calculation of weight or value with a solution like : sum(objects_value[i]*val for i, val in enumerate(sack))

This approach is valid, but will take a bit longer to run (less than we might imagine, though), because we replace the test of a Boolean value by a multiplication, and this for perfectible readability (even if it's closer to the linear program they saw during the previous Workshop).

sack = (True, False, True)

# we test the generator with a List Comprehension
objects =[\
 i for i, val in enumerate(sack) if val] #SOLUTION

print(objects)

[0, 2]

We obtain the list of indices of the objects we are interested in. Now, we are going to modify this generator to create the list of weight values of the objects present in the bag.

Explanations A hint : it is no longer i that we want to sum in the list, but weight_objects[i]. And this list, we will pass it to the function sum.

We will directly use this method to implement a function for calculating the weight of the bag. And by the way, we also implement the calculation of the value of the bag, on the same principle.

def sack_weight(sack):
 """
 This function returns the sum of the weights of the objects in the bag
 """
 # sum of the weights of the objects in the bag
 return\
 sum(objects_weight[i] for i, val in enumerate(sack) if val) #SOLUTION

def sack_value(sack):
 """
 This function returns the sum of the values of the objects in the bag
 """
 # sum of the values of the objects present in the bag
 return\
 sum(objects_value[i] for i, val in enumerate(sack) if val) #SOLUTION

Everything works correctly ? To check this, the easiest way is to generate a random solution, and to test our two functions on this solution.

import random

random.seed(a=3)
nb_objects = 5
objects_weight, objects_value = random_objects(10, 10)
print(objects_weight)
print(objects_value)

test = tuple(random.choice([False, True]) for _ in range(nb_objects))
print(test)
for o in range(nb_objects):
 print("%d:(%d,%d) "%(o, objects_weight[o], objects_value[o]), end='')

print("")
print([i for i, val in enumerate(test) if val]) # displays the indices of the True values of the solution
print("poids=%d, valeur=%d" %(sack_weight(test), sack_value(test)))

{0: 4, 1: 10, 2: 9, 3: 3, 4: 6}
{0: 10, 1: 8, 2: 10, 3: 2, 4: 10}
(False, True, True, False, False)
0:(4,10) 1:(10,8) 2:(9,10) 3:(3,2) 4:(6,10) 
[1, 2]
poids=19, valeur=18

test

(False, True, True, False, False)

Apparently, it works !

Explanations There is an implementation to filter the indices of the present objects that is less readable than Comprehension, but much faster. This implementation uses the compress function from the itertools library. This function does not accept a dictionary, so we must pass the values directly to it, for example by extracting them from the dictionary with the values method :

sum(itertools.compress(poids_objets.values(), test))

We can also directly replace dictionaries with tuples to represent values and weights, but the gain is insignificant.

We thus lose the access performance of the dictionary, but as the library is implemented in C, we still gain significantly (the calculation is 2 to 3 times faster, and on the set of the execution of the algorithm, we observe an average gain of 30 %). If the students are comfortable with the code, you can suggest that they look in this library for an alternative solution and test its performance.

2 simple heuristics

We can now focus on the implementation of a Hill Climbing. To begin with, it is necessary to implement the function for calculating the neighbourhood of a solution mentioned above.

2.1 Generation of neighbourhood

We therefore want to generate the set of all the neighbours of a solution (so, of a tuple of Booleans). To be sure, write this algorithm in more detail. You can consider that a function of cloning a solution is at your disposal.

Solution

It is necessary to implement a function that will take an input solution as a parameter, and which will sequentially generate all the valid neighbours of this solution. To do this, we only have to iterate over the number of objects in the instance and to build, at each iteration $i$, a tuple clone of the solution considered, whose presence in the cloned solution is reversed compared to the solution passed as a parameter.

function neighbourhood(bag)

Create the empty list of neighbours

For each object of the problem

Clone the bag
If the object is in the bag

Remove the object from the cloned bag

Otherwise:

Add object to the cloned bag

If cloned bag is a valid solution

Add the cloned bag in the list of neighbours

Return the list of neighbours

Now, what it left is to implement this algorithm.

Explanations We will implement this function of neighbourhood as a generator (it will make the implementation of the algorithm more memory efficient).

On the other hand, be careful, we have chosen tuples to represent our solutions. However, the tuples are immutable, which means that we cannot, as we wrote in the algorithm, implement a generator that clones and then modifies the current solution, but that we must generate only once a copy where a Boolean has been reversed. To do this, we will concatenate sliced tuples (pay attention to the bounds of slicing, in the notation Python, some are open, others closed).

And only valid solutions should be returned !

Explanations There is another approach than the one by slicing, which uses a Comprehension with ternary operator, which is passed to the function tuple :

voisin = tuple(v if i != k else not v for i, v in enumerate(sac))

Unless students are very familiar with Python’s advanced features, they are unlikely to have the idea for this approach. But if so, they can test both methods and compare their performance (using the same method as in Workshop 1 or the exercise series of sequence 0). They should realize that this implementation takes almost twice as long to run as the version with slicing, probably due to the enumerate(sac). But it's always a good approach to compare the running times of two possible implementations when doing scientific computing in Python.

def neighbourhood(sack):
 """
 This function is a generator of all valid neighbours of a solution
 """
 for k in range(len(sack)):
 # a hint: we want the states of the first k-1 objects and the last |sac|-k-1 objects
 # and the inverse state of object k (pay attention to the slicing bounds)
     neighbour = sack[:k] + (not(sack[k]),) + sack[k+1:] #SOLUTION
     if (sack_weight(neighbour) <= capacity): #SOLUTION
         yield(neighbour) #SOLUTION

# we test
test = (True, False, True)
print (test, "\nneighbours :")
for neighbour in neighbourhood(test):
 print(neighbour)

(True, False, True) 
neighbours :
(False, False, True)
(True, False, False)

Does that sound good to you ? Why aren't there more solutions ?

Because the other combinations of Booleans correspond to unfeasible solutions, as they exceed the bag’s capacity.

2.2 Hill Climbing

We can finally get down to implementing our heuristic Hill Climbing. If we were to write the algorithm in more detail, what could we propose ?

Solution

function HillClimbing (initial bag)

current bag ← initial bag
Do

Generate the neighbourhood of the current bag

For each bag in the neighbourhood

If the bag is of greater value than the current bag

The new bag becomes the current bag

As long as we find a new bag

Return the current bag

What is wrong with this algorithm ? How could it be fixed ?

Solution

It does not take into account the existence of a plateau, i.e. an area where the current solution stabilizes without degrading, and then eventually improves again.

We have to be careful about this problem, because we can’t simply replace the comparison of values with a test $\geq$ instead of $>$. In fact, in the case of several local optimums neighbouring each other, the algorithm would go back and forth between these solutions indefinitely.

A solution would be to do the comparison with $\geq$, but only allowing a certain number of iterations without improving the current solution.

And if students suggest setting up a memory to avoid backtracking, this ultimately comes down to implementing a Tabu list search, and this is the continuation of the workshop. Therefore, this problem is not very important, since it will be solved later with the Tabu list method.

Now let's get started with its implementation !

def hill_climbing(initial_solution):
 """
 1. We start from an element of our search set that we state current element
 2. We consider the neighbourhood of the current element and we choose the best of them
 as the new current element
 3. We loop until convergence on a local optimum
 """

 current_solution = initial_solution
 found_new_solution = True
 nb_iter = 0 # only used for display
    
 while (found_new_solution):
     nb_iter += 1
     # we go through all the neighbours of the current solution to keep the best one
     best_solution = current_solution #SOLUTION
     value_best = sack_value(current_solution) #SOLUTION
     #SOLUTION
     for neighbour in neighbourhood(current_solution): #SOLUTION
         neighbour_value = sack_value(neighbour) #SOLUTION
         if neighbour_value > value_best: #SOLUTION
             value_best = neighbour_value #SOLUTION
             best_solution = neighbour #SOLUTION
 #SOLUTION
     found_new_solution = (best_solution != current_solution)
     current_solution = best_solution
    
 return current_solution

2.3 Running the heuristic

Let's test this algorithm by taking an empty bag as an initial solution, and display the proposed solution :

random.seed(a=3)
capacity = 20
nb_objects = 100
objects_weight, objects_value = random_objects(10, 10)

sack = (False,)*nb_objects # a useful feature in Python

print("local optimization")
sol = hill_climbing(sack)
print("final value = " + str(sack_value(sol)) + ", capacite=" + str(sack_weight(sol)) + "/" + str(capacity))
print([i for i, val in enumerate(sol) if val]) # list of objects in the bag

local optimization
final value = 30, capacite=20/20
[5, 11, 27]

Three objects is not much ! That said, maybe we can do something to improve that. Think about it : this algorithm is deterministic. When considering an optimization algorithm like tabu search, what does that imply in terms of solution quality for a given instance ?

Solution

A deterministic algorithm will always produce the same solution for a given instance and a given starting solution (this is related to the deterministic character of the underlying Turing machine). The quality of the solution it generates depends solely on the starting solution.

Doesn't that give you an idea to improve the result produced by the algorithm ? How could one implement such an approach ?

Solution

We run the algorithm several times, starting from several different initial solutions, and we keep the best solution produced

It's simple to set up, especially in Python. Which function should you implement ? What exactly will she do ? A hint, Python’s random generation functions, e.g. random_choice should be helpful.

Solution

We only have to implement a function for generating a random valid solution. To do this, we can generate a random sequence of Booleans of size $n$ (the number of available objects). This solution is potentially unfeasible (too many objects in the bag), we will then randomly remove objects until the constraint of capacity of the bag is respected.

Ideally, such generator produces solutions according to a uniform distribution for the diversification step to work well. Here, the proposal is clearly not uniform, solutions with many objects will be over-represented. But uniformly generating a valid solution for the knapsack problem is a complex problem, which relies on advanced notions of probability (Markov chains, Monte Carlo method, etc.), which are beyond the students’ scope.

Note : In the same way as for the generation of the neighbourhood, it is possible to remove the selected object by Comprehension rather than by using the slicing, but there also it will be slower.

Let’s go !

Explanations This part makes it possible to approach the multi-start which will then be used on the metaheuristic. If the students are late, you can skip this part to go directly to optimizing the speed of execution of the metaheuristic. In fact, the following Workshop will not use the multi-start.

def random_solution():
 sack = tuple(\
     random.choice([False, True]) for _ in range(nb_objects)) #SOLUTION

 while (sack_weight(sack) > capacity):
      # list of `True` value indices in the bag
     selected_objects = tuple(
         i for i, val in enumerate(sack) if val) #SOLUTION
      # random choice of an object to delete among those present in the bag
     removed_object =\
         random.choice(selected_objects) #SOLUTION
        

 sack =\
 sack[:removed_object] + (False,) + sack[removed_object+1:] #SOLUTION
 # implementation: #SOLUTION
 # bag = tuple(val if i != objet_supprime else False for i, val in enumerate(sac)) #SOLUTION
 # also works but is a bit slower #SOLUTION

 return sack

Now we just have to test all this. To start, let's try with 10 iterations. At each iteration, display the search result (like the final display l.14 of the code to be completed, but without displaying the detail of the solution), so that we can see a little what is happening pass.

best_val = 0
for _ in range(5):
 sack = random_solution() #SOLUTION
 current_sol = hill_climbing(sack) #SOLUTION
 current_val = sack_value(current_sol) #SOLUTION
 print("value = %d, capacity=%d/%d" #SOLUTION
 %(sack_value(current_sol), #SOLUTION
 sack_weight(current_sol), #SOLUTION
 capacity)) #SOLUTION
 if (current_val > best_val): #SOLUTION
 best_val = current_val #SOLUTION
 best_sol = current_sol #SOLUTION

print("\nfinal value = %d, capacity=%d/%d\n" % (best_val, sack_weight(best_sol), capacity))
print([i for i, val in enumerate(best_sol) if val]) 

value = 40, capacity=20/20
value = 28, capacity=20/20
value = 33, capacity=20/20
value = 40, capacity=20/20
value = 44, capacity=20/20

final value = 44, capacity=20/20

[1, 27, 41, 62, 76, 81]

We have significantly improved the result ! And we can see that the local optimum found depends on the initial solution. By testing several initial solutions, one improves in certain cases the solution found. What is this technique called ?

Solution

The principle of this approach is called multi-start. Unlike simple Hill-Climbing, a Hill-Climbing with multi-start can be considered as a metaheuristic, since we have a intensification and a step of diversification. The intensification corresponds to the step of local optimization of the current solution and the diversification to the restart of a new starting solution.

Explanations Another existing metaheuristic exploits this idea of testing several starting solutions, but in a more sophisticated way. This is the GRASP (Greedy Randomized Adaptive Search Procedure) method, whose general principle is as follows :

function GRASP()

Repeat

Generate a valid solution $X$ per greedy randomized
Perform a local search starting from $X$
If Valeur(Recherche($X$)) > Valeur($meilleur$)

$meilleur \leftarrow X$

Up to the stop condition
Return $meilleur$

The difference between a simple multi-start and a GRASP is that the GRASP method uses a more elaborate approach to construct a random solution. The general idea is to iteratively build the starting solution, by randomly adding an element among the best $r$ according to a heuristic quality criterion (it should be emphasised that this criterion can only be heuristic, since the problem is $\mathcal{NP}$-Complete). This list of the best candidates, called RCL (Restricted Candidate List), is regenerated at each iteration :

function CréerRCL()

$Solution \leftarrow\{\emptyset\}$
Repeat

$RCL \leftarrow \{Candidats\}$
Sort $RCL$ by descending order of heuristic quality
$Candidat\leftarrow$ random element in $RCL[0..r]$
Add $Candidat$ to $Solution$

Until $Solution$ is complete
Return $Solution$

Here, one could consider as quality criterion heuristic the value/weight ratio of the objects. The choice of this ranking criterion is one of the criteria impacting the performance of a GRASP. The other parameter is the value of $r$. In the extreme cases, $r=1$ give a totally deterministic heuristic solution, $r=n$ (with $n$ the number of objects) gives a totally random solution, which is equivalent to a simple multi-start.

3 Metaheuristics

Finally, we have already implemented a first metaheuristic, but we are not going to stop there !

3.1 Search with a tabu list

A metaheuristic a little smarter than a Hill Climbing should give best results. Among all the metaheuristics in the literature, the tabu list search method (often called tabu method by shortcuts) is the one that could be adapted most easily based on the Hill Climbing code. Why ?

Solution

The general structure of the algorithm is very similar. We only have to adapt the search in the neighbourhood so it includes the tabu list management, and slightly modify the stop condition.

What algorithmic data structure do we need to handle this ?

Solution

A FIFO (First In, First Out) structure, i.e. a queue, with a maximum size.

To implement this in Python, which structure do you think is suitable from a performance point of view ?

Solution

The deque structure is exactly what we need. It performs addition operations at the beginning and at the end in O(1)$ instead of the $O(n)$ of other structures that could be used (typically lists), which represents a significant gain.

If students are comfortable with programming, there is a more efficient solution to running on large size lists, but a bit more subtle to implement. Adding in a deque is indeed very fast, but the defect of deque is that the search for a value is done in $O(n)$. For large sizes of list (which is what we will strive for), it is not ideal, it is better to be fast on the searches than on the addition.

The two structures that implement a search in $O(1)$ are set and dict. The first one seems ideal, but it has a huge flaw : the order of elements is not preserved, which disqualifies it from implementing a tabu list. This is however the case for dictionaries since Python version 3.6. One solution is therefore to create a dictionary whose keys are the tested solutions, and with None values. The only small difficulty is that you have to manage the maximum size of the list manually, by deleting the oldest element (hence the need to keep the order) if this maximum size is exceeded. With a tabu list of size 500, the gain varies between 10 % and 30 %, depending on the instance.

# we update the tabu list
tabu_list[current_solution] = None
if len(tabu_list) > tabu_size:
 tabu_list.pop(next(iter(tabu_list)), None)

Once again, this technique is a bit tricky to apply, to be reserved for people who are comfortable with code.

This time, we have everything we need. So let's get started with this implementation :

from collections import deque

def tabu_search(initial_solution, tabu_size, iter_max):
 """
 1. We start from an element of our search set that we state current element
 2. We consider the neighbourhood of the current element, we choose the best of them
 as a new current element, among those missing from the tabu list, and we add it
 to the tabu list
 3. We run a loop until the output condition is reached.
 """
 nb_iter = 0
 tabu_list =\
 deque((), maxlen = tabu_size) #SOLUTION

 # variables solutions for finding the neighbour optimal not tabu
 current_solution = initial_solution
 best = initial_solution
 global_best = initial_solution

 # variables values for finding neighbour optimal not tabu
 best_value = sack_value(initial_solution)
 global_best_value = best_value
    
 while (nb_iter < iter_max):
 best_value = -1
        
 # we go through all the neighbours of the current solution
 for neighbour in neighbourhood(current_solution):
 value_neighbour = sack_value(neighbour) #SOLUTION
 #SOLUTION
 # Update best non-tabu solution found #SOLUTION
 if value_neighbour > best_value and neighbour not in tabu_list: #SOLUTION
 best_value = value_neighbour #SOLUTION
 best = neighbour #SOLUTION

 # we update the best solution encountered since the beginning
 if best_value > global_best_value:
 global_best = best
 global_best_value = best_value
 nb_iter = 0
 else:
 nb_iter += 1

 # we go to the best non-tabu neighbour found
 current_solution = best

 # we update the tabu list
 tabu_list.append(current_solution) #SOLUTION

 return global_best

A test with a tabu list of size 5 on the same instance as before will allow us to compare the result with that obtained by local optimization :

nb_objects = 100
capacity = 20
random.seed(a=3)
objects_weight, objects_value = random_objects(10, 10)
sack = (False,)*nb_objects
tabu_size=5
iter_max=30

print("size 5 tabu list")
sol=tabu_search(sack, tabu_size, iter_max)
print("final value = " + str(sack_value(sol)) + ", capacity="+str(sack_weight(sol)) + "/" + str(capacity))

size 5 tabu list
final value = 62, capacity=20/20

We have significantly improved the quality compared to the solution obtained simply by local optimization ! That is to say if this first approach was ineffective… We are even better than with the multi-start !

3.2 Tabu with multi-start

Moreover, we can reuse this principle of multi-start. Implement it by strengthening diversification (500 iterations, for example) :

random.seed(a=3)

best_sol = None
best_val = 0 #SOLUTION
for _ in range (500): #SOLUTION
 sack = random_solution() #SOLUTION
 current_sol = tabu_search(sack, tabu_size, iter_max) #SOLUTION
 current_val = sack_value(current_sol) #SOLUTION
 if (current_val > best_val): #SOLUTION
 best_val = current_val #SOLUTION
 best_sol = current_sol #SOLUTION
print("final value = " + str(sack_value(best_sol)))

final value = 89

The solution is further improved. That said, the overhead in computation time is substantial.

Explanations This is problematic:

• if the algorithm must run on a very large amount of data
• if we run this algorithm a large number of times, to statistically assess its performance (which is the subject of the next Workshop)

However, there is a way to speed up the calculations without touching much. Some of the functions we have implemented are pure. And that will allow us to set up an acceleration of the execution of the code. Do you have any idea how to do it ?

Solution

A hint : we will have to modify the space–time trade-off that we made on the neighbourhood calculation.

The functions for calculating the weight and the value of a solution are pure, so they can be memoized. For the function of generation of the neighbourhood, it’s also feasible (and desirable), but not in the state, because this function is an iterable, so it needs to be reset with each new use. The function needs to be modified to generate a list that it will return at once. The code of the metaheuristic does not need to be modified.

This is a typical case where we change the space–time trade-off : we increase the memory used to reduce the computation time. This explains why the implementation of a generator must be abandoned here, it corresponds to a solution that is too economical in memory, and not fast enough.

Note : whether these functions are pure is actually debatable, since they depend on the global variables that describe the instance (objects_value, objects_weight, nb_objects, and capacity). If the students notice this, it can be pointed out that these values do not change, so it is not a problem. On the other hand, in the following Workshop, we will vary them to obtain statistical results. At this point, it will cause problems if these functions are memoized, because they will not produce the right values. This point will be dealt with in due course.

Moreover, fortunately we chose a tuple to represent a solution, because we could not have implemented this technique with lists (which are not hashable).

First of all, let's restart the calculation by counting the time CPU, to get an idea of the gain that will be obtained afterwards.

Explanations Another possibility that students can think of to speed up calculations is to parallelize the multi-start. Attention, in this case, the use of Threads in Python will not be interesting (apart from the fact that the generators, used for the neighbourhood of a solution, are not thread-safe) : the Global Interpreter Lock prevents the simultaneous execution of several Threads (a big flaw of Python). Only I/O can be parallelized with processing.

A real parallelization of calculations requires going through heavy process, with all the problems that this implies (including the slowdown due to context switching context switch, which makes an acceleration proportional to the number of cores seems out of reach).

And in this case, the benefit of memoization is clearly reduced (since each process will have its own cache). It is not rare that a sequential multi-start with memoization is more powerful than a parallel multi-start by process with memoization (here, it is the case).

There remains one last solution, which consists of implement a cache sharing of memoization between process via a shared memory, but this begins to become quite heavy in terms of development.

import time
random.seed(a=3)

best_sol = None
best_val = 0

start = time.process_time()

# multi-start of 500 iterations #SOLUTION
for _ in range (500): #SOLUTION
 sack = random_solution() #SOLUTION
 current_sol = tabu_search(sack, tabu_size, iter_max) #SOLUTION
 current_val = sack_value(current_sol) #SOLUTION
 if (current_val > best_val): #SOLUTION
 best_val = current_val #SOLUTION
 best_sol = current_sol #SOLUTION

stop = time.process_time()

print("final value = " + str(sack_value(best_sol)))
print("computed in", stop-start, 's')

final value = 89
computed in 11.25 s

Now let's test this technique, by redefining the functions involved.

from functools import lru_cache

@lru_cache(maxsize=None) #SOLUTION
def neighbourhood(sack): #SOLUTION
 neighbours = [] #SOLUTION
 for k in range(len(sack)): #SOLUTION
 neighbour = sack[:k] + (not(sack[k]),) + sack[k+1:] #SOLUTION
 if (sack_weight(neighbour) <= capacity): #SOLUTION
 neighbours.append(neighbour) #SOLUTION
 return neighbours #SOLUTION
 #SOLUTION
@lru_cache(maxsize=None) #SOLUTION
def sack_weight(sack): #SOLUTION
 return sum(objects_weight[i] for i, val in enumerate(sack) if val) #SOLUTION
 #SOLUTION
@lru_cache(maxsize=None) #SOLUTION
def sack_value(sack): #SOLUTION
 return sum(objects_value[i] for i, val in enumerate(sack) if val) #SOLUTION 

And let's test this by rerunning the previous calculation :

random.seed(a=3)

best_sol = None
best_val = 0

start = time.process_time()

# multi-start of 500 iterations #SOLUTION
for _ in range (500): #SOLUTION
 sack = random_solution() #SOLUTION
 current_sol = tabu_search(sack, tabu_size, iter_max) #SOLUTION
 current_val = sack_value(current_sol) #SOLUTION
 if (current_val > best_val): #SOLUTION
 best_val = current_val #SOLUTION
 best_sol = current_sol #SOLUTION

stop = time.process_time()

print("final value = " + str(sack_value(best_sol)))
print("computed in", stop-start, 's')

final value = 89
computed in 2.578125 s

The acceleration is impressive, the computation time has been almost divided by 4 ! And the more calculations we make on this instance, the more we will accelerate (in exchange for a higher memory consumption). To convince you of this, just run the same calculation again, but with a different seed :

random.seed(a=5)

best_sol = None
best_val = 0

start = time.process_time()

# multi-start of 500 iterations #SOLUTION
for _ in range (500): #SOLUTION
 sack = random_solution() #SOLUTION
 current_sol = tabu_search(sack, tabu_size, iter_max) #SOLUTION
 current_val = sack_value(current_sol) #SOLUTION
 if (current_val > best_val): #SOLUTION
 best_val = current_val #SOLUTION
 best_sol = current_sol #SOLUTION

stop = time.process_time()

print("final value = " + str(sack_value(best_sol)))
print("computed in", stop-start, 's')

final value = 89
computed in 2.40625 s

There too, we went faster, while the execution scenario is not the same, which means that the initial solutions considered are different.

Explanations This clearly shows that the memoization is more and more efficient as the functions are used. This makes sense, since each time we encounter a solution that we had already visited, we do not need to recalculate its weight, its value, and its neighbourhood.

And if you're wondering why we still get the same solution value, we'll have the answer at the next Workshop !

4. Conclusion

This Workshop is now over, and we have finally managed to do something of this knapsack problem, since you have implemented, not one, but three metaheuristics, a Hill-Climbing with multi-start, a simple tabu search, and one with multi-start.

Explanations Finally, if the theoretical presentation of a metaheuristic generally gives the impression of a complicated method, the implementation is often quite simple, especially with a language like Python. And, as long as we use the language in the right way, we keep quite satisfactory performances.

In this regard, many of the functions that we have implemented can be implemented in several different ways. The use of timeit and [cProfile] (https://www.machinelearningplus.com/python/cprofile-how-to-profile-your-python-code/) is more than recommended (even though reading the cProfile statistics in a Notebook is a bit tedious). Only with these tools will you be able to identify bottlenecks that require optimization work (after all, what's the point of optimizing code that is only responsible for a small portion of the total execution time? execution ?).

Even if you start with a different implementation for your project, you can see that it is not that complicated. Take the case of a genetic algorithm. In the case of a knapsack, what would a gene represent ?

Solution

A gene is the Boolean value representing the presence of an object in the bag.

If we had to implement, we would already have a lot of reusable functions (random solution generation, assessment of the feasibility of a solution, of its value). It would just be necessary to implement a crossing function and a mutation function. What kind of crossover could be considered ?

Solution

The simple solutions for crossover are generally to take the $i$ first genes of one individual, and the $n-i+1$ last genes of the other individual (with $i$ chosen randomly) :

Parent A													Parent B
0	1	1	0	0	1	0	0	1	1	0	0	+	1	0	1	0	0	1	0	1	0	1	1	0
⇩

Child
0	1	1	0	0	1	0	1	0	1	1	0

A slightly more advanced solution consists in taking the first $i$ genes and the last $n-j+1$ genes of the first individual, and the $i$ to $j$ genes of the second (with random $i$ and $j$, and $i<j$) :

Child
0	1	1	0	0	1	0	1	0	1	0	0

What about the mutation ?

Solution

Since the genes here are Booleans, mutating a gene consists in reversing its value.

Once these additional functions are implemented, the body of the algorithm will not be much more complicated than that of the tabu algorithm that we have just written.

Explanations The real challenge is not really in the choice of a metaheuristic, but especially in its tuning, to improve either the computation time or the quality of the solution. In a counter-intuitive way, these will strongly depend on the problem on which we are working and on its parameters (no free lunch theorem).

But the question that arises is : are the solutions we build of good quality ? Ultimately, we don't know.

Explanations If we are already very close to the optimal, it will not be of much interest to spend time implement optimizations that will only save us very little, especially if they induce a significant additional cost to the execution (or even to the design).

And even if we are far from the optimal, what to do ? There are techniques to improve the step of intensification, others which are interested in that of diversification, but in our case, which ones are we interested in ?

It will be necessary to analyse the behaviour of the algorithm according to certain parameters, to be able to determine if it is effective or not, and if necessary to optimize its operation. And that’s the subject of the next Workshop !