AMMM - COURSE PROJECT
Authors: Toni Ivanković, Alessandro Scannavini
The task was to implement an algorithm which would find a set of natural numbers in which
no two pairs of numbers share the same difference. The problem is known as a Golomb
ruler: such a set can be used in acoustics, where a period of a sinusoidal wave can be found
by measuring the pressure at the times specified by the numbers in the set. Additional
applications include radio antenna placement and also conference room wall configuration.
FORMAL STATEMENT
Variables:
𝑥𝑖 𝑠. 𝑡. 𝑖 = 1... 𝑁 - i-th integer in the set of distinct numbers in ascending order
Objective function:
𝑚𝑖𝑛𝑖𝑚𝑖𝑧𝑒 𝑥𝑁 − 𝑥1
*As described, the goal is to minimize the difference between the largest and the
smallest number in the set. For easier formulation, the numbers are in ascending
order, so the largest number is 𝑥𝑁 and the smallest is 𝑥1.
Constraints:
Ascending order:
∀𝑖 𝑠. 𝑡. 1 ≤ 𝑖 ≤ 𝑁 − 1, 𝑥𝑖 < 𝑥𝑖+1
Unique differences:
∀𝑖, 𝑗, 𝑘, 𝑙 𝑠. 𝑡. 1 ≤ 𝑖 ≤ 𝑁, 1 ≤ 𝑘 ≤ 𝑁, 1 ≤ 𝑗 < 𝑖, 1 ≤ 𝑙 < 𝑘,
¬(𝑘 = 𝑖 ∧ 𝑙 = 𝑗) ⇒ 𝑥𝑖 − 𝑥𝑗 ≠ 𝑥𝑘 − 𝑥𝑙
ILP MODEL
In the ILP model, only linear constraints with inequalities may be used. The model is, hence,
a bit more complex than the formal statement.
Variables:
𝑥𝑖 𝑠. 𝑡. 𝑖 = 1... 𝑁
- (integer variable) i-th integer in the set of distinct numbers in ascending order
𝑝𝑜𝑠𝑖𝑡𝑖𝑣𝑒𝑖,𝑗,𝑘,𝑙 𝑠. 𝑡. 𝑖, 𝑗, 𝑘, 𝑙 = 1... 𝑁
- (boolean variable) defined as true iff 𝑥𝑖 − 𝑥𝑗 − ( 𝑥𝑘 − 𝑥𝑙) > 0
𝑛𝑒𝑔𝑎𝑡𝑖𝑣𝑒𝑖,𝑗,𝑘,𝑙 𝑠. 𝑡. 𝑖, 𝑗, 𝑘, 𝑙 = 1... 𝑁
- (boolean variable) defined as true iff 𝑥𝑖 − 𝑥𝑗 − ( 𝑥𝑘 − 𝑥𝑙) < 0
Objective function:
𝑚𝑖𝑛𝑖𝑚𝑖𝑧𝑒 𝑥𝑁 − 𝑥1
Constraints:
Ascending order:
∀𝑖 𝑠. 𝑡. 1 ≤ 𝑖 ≤ 𝑁 − 1, 𝑥𝑖 + 1 ≤ 𝑥𝑖+1
Unique differences:
∀𝑖, 𝑗, 𝑘, 𝑙 𝑠. 𝑡. 1 ≤ 𝑖, 𝑘 ≤ 𝑁, 1 ≤ 𝑗 < 𝑖, 1 ≤ 𝑙 < 𝑘,
¬(𝑘 = 𝑖 ∧ 𝑙 = 𝑗) ⇒
[ 𝑝𝑜𝑠𝑖𝑡𝑖𝑣𝑒𝑖,𝑗,𝑘,𝑙 + 𝑛𝑒𝑔𝑎𝑡𝑖𝑣𝑒𝑖,𝑗,𝑘,𝑙 ≥ 1 ∧
𝑥𝑖 − 𝑥𝑗 − ( 𝑥𝑘 − 𝑥𝑙) − 1 ≥ − 10000 · (1 − 𝑝𝑜𝑠𝑖𝑡𝑖𝑣𝑒𝑖,𝑗,𝑘,𝑙) ∧
𝑥𝑖 − 𝑥𝑗 − ( 𝑥𝑘 − 𝑥𝑙) ≤ 10000 · 𝑝𝑜𝑠𝑖𝑡𝑖𝑣𝑒𝑖,𝑗,𝑘,𝑙 ∧
− (𝑥𝑖 − 𝑥𝑗 − ( 𝑥𝑘 − 𝑥𝑙)) − 1 ≥ − 10000 · (1 − 𝑛𝑒𝑔𝑎𝑡𝑖𝑣𝑒𝑖,𝑗,𝑘,𝑙) ∧
− (𝑥𝑖 − 𝑥𝑗 − ( 𝑥𝑘 − 𝑥𝑙)) ≤ 10000 · 𝑛𝑒𝑔𝑎𝑡𝑖𝑣𝑒𝑖,𝑗,𝑘,𝑙 ]
*In order to more easily formulate the constraints, the variables were constrained to
be in ascending order, meaning that each variable must be smaller by at least 1 than
its successor.
As for the unique differences, if two different pairs are chosen (it's not the case that
both their indices match), their difference has to be either positive or negative (1st
clause of the consequent).
Other clauses of the consequent are there to give the intended meaning to the
variables 𝑝𝑜𝑠𝑖𝑡𝑖𝑣𝑒 and 𝑛𝑒𝑔𝑎𝑡𝑖𝑣𝑒 - the 2nd clause of the consequent ensures that if
𝑝𝑜𝑠𝑖𝑡𝑖𝑣𝑒 = 1, the difference 𝑥𝑖 − 𝑥𝑗 − ( 𝑥𝑘 − 𝑥𝑙) is positive, the 3rd clause ensures
that if the difference is positive, then 𝑝𝑜𝑠𝑖𝑡𝑖𝑣𝑒 = 1, and the 4th and the 5th clause
achieve an analogous effect, but with 𝑛𝑒𝑔𝑎𝑡𝑖𝑣𝑒.
METAHEURISTICS
Constructive algorithm:
The function greedyAssign recursively appends the first number it can find that, when added,
does not produce any duplicate difference. The optimal solution is not guaranteed, but a
feasible solution is found quickly. The maximum number used is 10000, since the application
becomes too slow for the n-value big enough to come close to the ceiling. If needed, the
ceiling can easily be increased, or a while loop can be used to infinitely add to the i variable.
function greedyAssign(assignment, n):
if assignment.length == n and noCollisions(assignment) then
return assignment
end
for i = 1 to 10000 do
if i not in assignment and noCollisions(assignment +
[i]) then
return greedyAssign(assignment + [i], n)
end
end
The greedy cost function is embedded in the algorithm implicitly, but could be formalized as
𝑞(𝑐, 𝑆) = 𝑐 − 𝑚𝑎𝑥(𝑆). The new 𝑐 will by construction be larger than all of the elements in
𝑆, so the new objective function value will be the difference between it and the existing
largest element. Hence, the new objective value will be larger than the last objective function
value by 𝑐 − 𝑚𝑎𝑥(𝑆) .
The function noCollisions takes an assignment (a list of numbers) and checks whether any of
the differences between distinct pairs of numbers match. If that is the case, false is returned.
Otherwise, if all differences are checked and no collision has been found, true is returned.
function noCollisions(assignment):
for i = 0 to assignment.length do
for j = i+1 to assignment.length do
for k = 0 to assignment.length do
for l = k+1 to assignment.length do
if ¬(i==k ∧ j==l) and |assignment[i] -
assignment[j]| == |assignment[k] -
assignment[l]| then
return false
end
end
end
end
end
return true
Local search:
Neighboring states of an assignment are defined as all assignments in which all the numbers
differ by at most 1 from their counterparts in the original assignment. More formally:
𝑁(𝑆) = {𝑥 : ∃𝑥' ∈ 𝑆 𝑠. 𝑡. |𝑥' − 𝑥| ≤ 1} and all elements from S have a corresponding
element in N(S) which differs from it by at most 1.
𝑛
This is achieved by taking the set S, sorting it, and adding to it a vector from {− 1, 0, 1} .
All the neighbors of the initially constructed solution are searched and the best found is
reported.
The local search function itself does not iteratively update the solution, but just searches the
neighborhood of the given assignment. Searching until a local optimum is achieved by
repeatedly calling the function until the solution does not improve anymore.
It is worth noting that using this approach, the new set could contain 1 instead of 0. Even
though that would also be a perfectly feasible solution, for the sake of comparing solutions,
any solution found by the local search is afterwards shifted until its smallest element is 0.
This can be done because subtracting the same value from all elements does not tamper
with the validity of any constraints.
function localSearch(assignment, n, visited):
localBest := assignment
localBestValue := evaluateAssignment(assignment)
𝑛
for addition in {− 1, 0, 1} do
if addition[0] ≤ 0 or addition[n-1] ≥ 0 then
continue
assignmentCopy := assignment.add(addition)
if assignmentCopy not in visited and
feasibleAssignment(assignmentCopy, n) then
visited.add(assignmentCopy)
newValue := evaluate(assignmentCopy)
if newValue < localBestValue then
localBest := assignmentCopy
localBestValue := newValue
end
return localBest
function evaluate(assignment):
return max(assignment) - min(assignment)
function feasibleAssignment(assignment, n):
if assignment.length == n then
return noCollisions(assignment)
else
return false
GRASP:
The GRASP algorithm accepts the hyperparameters in the beginning, constructs a greedy
solution, and then iteratively tries to improve it.
In each iteration, an initial assignment is constructed and its neighborhood is searched.
The construction of the initial assignment in each iteration is done by repeatedly constructing
the RCL set of numbers, and taking the new value from a subset of the RCL set. The size of
the subset is determined by the alpha parameter. The formal definition of the feasible value
and RCL sets are:
𝐹𝑉(𝑆) = {𝑥: 𝑥 ∉ 𝑆 ∧ 𝑛𝑜𝐶𝑜𝑙𝑙𝑖𝑠𝑖𝑜𝑛𝑠(𝑆 ∪ {𝑥}) ∧ 𝑥 ≤ 𝑢𝑝𝑝𝑒𝑟𝐵𝑜𝑢𝑛𝑑}
𝑅𝐶𝐿(𝑆) = {𝑥: 𝑥 ∈ 𝐹𝑉(𝑆) ∧ 𝑥 ≤ [𝑚𝑎𝑥(𝐹𝑉(𝑆)) − α · (𝑚𝑎𝑥(𝐹𝑉(𝑆)) − 𝑚𝑖𝑛(𝐹𝑉(𝑆)))]}
Setting the alpha parameter to 0 produces a purely random search, and setting it to 1
produces a purely greedy search. The values in between produce partial randomness in
choosing each value in the assignment, which enables the algorithm not to get stuck in the
local optimum of the greedy solution.
To improve the performance, the RCL set was constrained not to include x values greater
than the current best solution. Such a modification does not eliminate any better solutions,
and it improves performance because any constructed assignment is either equal to or better
than the current best assignment by its objective function.
procedure GRASP(n, alpha, MAX_ITER):
greedySolution := greedyAssign([0], n)
bestSolution := greedySolution
bestSolutionValue := evaluateAssignment(greedySolution)
upperBound := greedySolution[n-1]
for i=0 to MAX_ITER do
assignment := semiGreedyAssign([0], alpha, upperBound)
if assignment exists then
localBest := localSearch(assignment,n, set())
localBestValue := evaluateAssignment(localBest)
if localBestValue < bestSolutionValue then
bestSolution := localBest
bestSolutionValue := localBestValue
upperBound := bestSolution[n-1]
done
return bestSolution
function semiGreedyAssign(assignment, alpha, upperBound):
assignmentCopy = copy(assignment)
while assignmentCopy.length < n do
RCL := getFeasibleValues(assignmentCopy, upperBound)
if RCL.length == 0 then
return None
qmin := RCL[0]
qmax := RCL[-1]
borderValue := qmax - alpha*(qmax-qmin)
RCL_reduced := list()
for x in RCL:
if x <= borderValue then
RCL_reduced.append(x)
done
RCL := RCL_reduced
index := random_int(0, RCL.length - 1)
value := RCL[index]
assignmentCopy.append(value)
done
if checkNoCollisions(assignmentCopy) then
return sort(assignmentCopy)
else
return None
funtcion getFeasibleValues(assignment, upperBound):
feasibleValues = list()
for i:=0 to upperBound+1:
if i not in assignment and noCollisions(assignment +
[i])):
feasibleValues.append(i)
return feasibleValues
RESULTS
Since setting the alpha parameter of GRASP to 1 gives a full greedy approach, there was no
need to perform the greedy search by itself, so only the GRASP experiments are reported.
Furthermore, different values of the hyperparameter MAX_ITER (the number of iterations of
the GRASP algorithm) were explored. However, setting a bigger MAX_ITER will, on average,
produce only equally good or better solutions. In that sense, there is no tradeoff with solution
quality if the parameter MAX_ITER is increased more than necessary.
Generally, bigger values of MAX_ITER were needed to reach better solutions as n and alpha
grew larger. For 𝑛 ≤ 10 the optimal solution could most often be found in less than 1000
iterations, while for larger n values, oftentimes more iterations were needed.
The biggest challenge as the n increases is the increase in time complexity. The time
increase is shown in the following graph, on a logarithmic scale. A nearly straight line in the
logarithmic scale shows an exponential trend as the n increases. The time is expressed in
seconds to obtain the solution.
It is worth noting that the exponential time needed by GRASP is much lower than the one
needed by CPLEX. CPLEX needed 50 minutes on average to compute the value for n=10,
while GRASP needed around 1 minute. Even though the asymptotic complexity should be
the same, this difference can be significant in practical applications.
Furthermore, different alpha values were explored to see how well they perform in reaching
the optimal value. The best solution over all the tested MAX_ITER was used in these plots.
To get a better sense of the quality of the solution, the ratio between the obtained solutions
and the optimal CPLEX solution was plotted as well.
It is clearly visible that alpha=1 (completely greedy approach) gives by far the worst solution,
but from the plots it is not clearly visible which alpha value would universally be the best. The
interest was to find the smallest number of iterations that produces the best solution. With
that goal, heatmaps showing average error (Table 1) and execution time by alpha and max
iterations (Table 2) were produced.
In the first heatmap the average error by number of iterations is shown with an additional
column describing the average error over each alpha value.
𝐻𝐸𝐴𝑇𝑀𝐴𝑃 𝐴𝑉𝐺 𝐸𝑅𝑅𝑂𝑅
NUMBER OF ITERATIONS
100 200 500 1000 2000 5000 AVERAGE
0 5,1 4,3 4,2 2,4 2,1 1,7 3,3
0,1 6,8 6,2 4,1 3,9 3,1 2,4 4,42
0,2 7,5 5,6 5,1 3,8 3,9 1,9 4,63
A 0,3 6,7 5,2 4 4,8 2,4 1,5 4,10
L
P 0,4 5,8 4,6 4,8 3,5 2,6 0,8 3,68
H
0,5 4,9 5,5 3,5 2,8 3,1 1,7 3,58
A
0,6 4,6 4,6 4,2 2 2,1 1,9 3,23
0,7 5,6 3,9 2,6 2 2,1 2,1 3,05
0,8 5 3,8 2,1 2,9 2,7 2,2 3,12
0,9 4,7 4,1 3,9 3,9 3,3 3,3 3,87
1 11,9 11,9 11,9 11,9 11,9 11,9 11,9
(Table 1)
In the second heatmap, the average time multiplied by the average error squared is shown
instead (the elevation is done in order to put more importance on the error, which was not
that influential without it). As in the previous heatmap, an additional column describing the
average values over each alpha value it’s shown.
2
𝐴𝑉𝐺 𝑇𝐼𝑀𝐸 𝑥 𝐴𝑉𝐺 𝐸𝑅𝑅𝑂𝑅
NUMBER OF ITERATIONS
100 200 500 1000 2000 5000 AVG
0 484,84 304,32 582,53 246,67 276,84 335,05 371,71
0,1 410,97 429,11 397,44 456,53 680,55 751,46 521,01
0,2 803,52 359,93 368,84 1095,86 724,65 365,77 619,76
0,3 354,34 216,29 217,72 719,50 251,74 209,63 328,20
A
0,4 393,69 393,34 598,18 423,88 359,15 68,19 372,74
L
P 0,5 237,13 325,64 248,88 274,78 505,98 284,35 312,79
H
A 0,6 356,26 788,60 319,01 80,14 336,65 331,12 368,63
0,7 272,15 257,41 168,34 181,99 580,35 384,36 307,43
0,8 284,54 165,14 192,06 324,77 811,74 531,03 384,88
0,9 153,74 424,68 362,17 472,09 1259,09 1233,77 650,92
1 433,36 865,76 2160,20 4130,72 8325,60 20522,31 6072,99
(Table 2)
From the last heatmap, alpha=0.7 was identified as the optimal hyperparameter value on
average. Another good combination was alpha=0.6 with 1000 iterations and alpha=0.4 with
5000 iterations. However, on average, the performance of those alpha values was worse.
The values obtained by alpha=0.7 with 500 iterations were compared to the CPLEX solution
(for n>10 a solution from the Internet was used) and plotted below:
It is visible that, at worst, the GRASP solution with these hyperparameters was around 15 %
higher than the optimal solution, which could be considered as an acceptable tradeoff for a
decrease in execution time.
