Algorithm 2 The multi-state heuristic float x = 1.3; float pr2[i in 1..12] = (x^i)*(x-1.0)/(x^(13.0)-x); Closure step = closure {
if(state == 1){
selectPr(i in 1..12)(pr2[i])
select(d in ((12-i)*10)..((13-i)*10): d>0)
selectMin(f in RWF, ad = S.getAssignDelta(deltaT[f],d):
tabu[f] <= it && S.violations(deltaT[f]) > 0 && deltaT[f] != d && ad < 0)(ad){
deltaT[f] := d;
}
}
else if(state == 2){
selectMax(f in RWF: tabu[f] <= it && S.violations(deltaT[f]) > 0)
(S.violations(deltaT[f])){
selectMin(d in 0..g, ad = S.getAssignDelta(deltaT[f],d): deltaT[f] != d && ad < 0)(ad,d){
deltaT[f] := d;
}
}
}
else if(state == 3){
selectMin(f in RWF, d in 0..g, ad = S.getAssignDelta(deltaT[f],d): tabu[f] <= it && S.violations(deltaT[f]) > 0 && deltaT[f] != d && ad < 0)(ad,d){
deltaT[f] := d;
}
}
};
We can define the required exponential probability distribution function by dividing the geometric series by the sum of its term:

f(x, y) = x
$$v)={\frac{x^{y}\cdot(x-1)}{x^{n+1}-x^{m}}},$$
$$f(x,y)=1$$
$$\leq y\leq n$$
, with m ≤ y ≤ n &
n
∑
y=m
f(x, y) = 1 (4)
where m and n are the powers of the first and last terms in the geometric series, respectively; these values set the *scale* of the probability distribution. Variable x is a real number called the *ratio* of the series; it affects the speed of growth of the function. Variable y is an integer ranging from m to n. The delay domain [0*,...,*g] being potentially large, we divide it into ten-minute-steps by making y range from m = 1 to n = 12 when g = 120; the ratio was experimentally set to x = 1.3. Algorithm 2 depicts the Comet implementation of the heuristic.

In the first state, when the violation of the constraints falls below an experimentally determined threshold of 300, the search procedure switches to the second state, which selects the most violated