Packages Used: No new packages are used in this notebook.
The effectiveness of a heuristic can be measured by comparing the result obtained from the heuristic to an equal number of brute force evaluations. Using the QAP considered in Comb Opt 2, a larger instance (n = 25) is randomly generated. Both the W and D generated are not non-symmetric and have zeros along their main diagonals.
using Random, LinearAlgebra
Random.seed!(543453)
n = 25
W = rand(n,n)
W = triu(W,1) + tril(W,-1)
D = rand(n,n)
D = triu(D,1) + tril(D,-1)
25×25 Matrix{Float64}:
0.0 0.734633 0.295956 … 0.111347 0.157435 0.949252
0.774889 0.0 0.0606261 0.981961 0.209453 0.632379
0.196373 0.514708 0.0 0.876814 0.676663 0.0913014
0.100732 0.00135112 0.477068 0.49046 0.611562 0.78657
0.434407 0.34739 0.887897 0.119165 0.543803 0.404344
0.875129 0.0942881 0.775923 … 0.234693 0.432229 0.150349
0.718132 0.999431 0.858756 0.5132 0.904724 0.934413
0.93238 0.57109 0.711163 0.147829 0.43591 0.378114
0.916641 0.693598 0.23702 0.61198 0.527348 0.491715
0.904217 0.494604 0.234618 0.989766 0.617034 0.255759
0.0119593 0.0222515 0.867342 … 0.818101 0.445842 0.953445
0.259317 0.3544 0.786769 0.007798 0.558257 0.403668
0.481527 0.275997 0.159614 0.045784 0.945741 0.451291
0.431843 0.461648 0.373341 0.247006 0.130917 0.750905
0.818006 0.370336 0.00819065 0.723119 0.13314 0.814841
0.101392 0.913905 0.707019 … 0.0684296 0.313747 0.436001
0.814685 0.570696 0.809299 0.648238 0.583333 0.845044
0.554121 0.115248 0.0851838 0.56172 0.729648 0.964829
0.744955 0.949927 0.0194064 0.942244 0.00838548 0.411264
0.34982 0.73143 0.0306511 0.241418 0.228093 0.325457
0.910521 0.930545 0.157131 … 0.350923 0.702914 0.596449
0.757109 0.655005 0.0556803 0.949813 0.473079 0.91058
0.755232 0.354019 0.424286 0.0 0.818444 0.235044
0.136098 0.164871 0.816721 0.67091 0.0 0.752949
0.554567 0.611774 0.0880611 0.107 0.453536 0.0
Since each element of W and D ranges from 0 to 1, an upper bound on the solution would correspond to the case where all of the elements of each matrix are at 1:
TC_UB = n^2
625
The TC of the expected random solution would correspond to each element of W and D being equal to 0.5:
TCexp = (0.5)*(0.5)*n^2
156.25
function sdpi(αᵒ,W,D)
TCᵒ, n, cnt = sum(W[αᵒ,αᵒ].*D), length(αᵒ), 1
done = false
while !done
done = true # Assume no improvement will be found
α₀ = copy(αᵒ) # Copy as current in case improvement found
for i = 1:n-1, j = i+1:n
cnt += 1
α = copy(α₀) # Work on copy of current assignment
α[[j,i]] .= α[[i,j]] # Pairwise interchange of i and j
TC = sum(W[α,α].*D) # Calc. TC of interchange
if TC < TCᵒ # Improvement found
done = false
TCᵒ, αᵒ = copy(TC), copy(α)
end
end
end
return (TC = TCᵒ, α = αᵒ, cnt = cnt) # Named tuple used to pass multiple outputs
end
function multisdpi(W,D)
TCᵒ, αᵒ, n, cnt = Inf, [], size(W,1), 0
for i = 1:n
res = sdpi(randperm(n),W,D)
cnt += res.cnt
if res.TC < TCᵒ
(TCᵒ,αᵒ) = res
end
end
return (TC = TCᵒ, α = αᵒ, cnt = cnt)
end
Random.seed!(3242)
res = multisdpi(W,D)
(TC = 136.88790675024478, α = [20, 19, 2, 5, 24, 15, 1, 21, 8, 22 … 3, 13, 4, 10, 7, 11, 17, 18, 12, 16], cnt = 108325)
In order to make a brute-force solution commensurate with the heuristic solution, the number of times the cost of a solution was evaluated in sdpi (i.e., TC was calculated) was output as cnt. This number of random solutions was then evaluated:
factorial(big(n))
15511210043330985984000000
res.cnt
108325
TCmin,TCmax,TCest = Inf,-Inf,0
for i = 1:res.cnt
α = randperm(n)
TC = sum(W[α,α].*D)
TCest += TC
if TC < TCmin
TCmin = TC
elseif TC > TCmax
TCmax = TC
end
end
TCest /= res.cnt
res.TC,TCmin,TCexp,TCest,TCmax,TC_UB
(136.88790675024478, 145.56531871064269, 156.25, 154.75799502592483, 164.3423276331147, 625)
TCmin/res.TC
1.063390639585351
Scheduling decisions involve both (1) the allocation of tasks to resources and (2) determining, for each resource, the sequence in which its allocated tasks should occur. Unlike bin packing, the number of resources (bins) is predetermined, and unlike QAP, a resource (cf. site) can be allocated (cf. assigned) more than one task. If there is only a single resource, then only the sequencing decision needs to be determined; for multiple resources, if all of the tasks are independent of each other and do not have due dates or priorities, then only the allocation decision needs to be made.
Given three identical machines (resources) and eight independent jobs (tasks) 1 to 8 that each take 1 to 8 units of time to complete, respectively, determine how the jobs should be allocated to the machines so that all processing is completed as soon as possible (i.e., minimize the makespan). Once processing starts on a machine, a job cannot be preempted.
m = 3 # Number of resources
n = 8 # Number of tasks
t = [1:n;] # Time of each task
8-element Vector{Int64}:
1
2
3
4
5
6
7
8
With preemption, it is easy to determine an optimal allocation by just filling each resource in sequence, and when a task exceeds the lower bound, stop processing at that point, move the preempted task to the beginning of processing on the current resource, and then allocate the remaining processing later on a different resource. Without preemption, it may not be possible to reach the LB.
LB = sum(t)/m
12.0
For this simple problem instance, an optimal solution can be generated by just pairing tasks on each machine so that they total the LB:
T = [[8,4], [7,5], [1,2,3,6]]
3-element Vector{Vector{Int64}}:
[8, 4]
[7, 5]
[1, 2, 3, 6]
An array comprehension can be used to determine the total time on each resource, and then taking the maximum gives the makespan:
[sum(i) for i in T]
3-element Vector{Int64}:
12
12
12
maximum([sum(i) for i in T]) # Makespan
12
The number of ways that n distinguishable balls can be placed into m indistinguishable cells can be determined using the Stirling number of the second kind:
$ \begin{eqnarray} \quad \left\{{n \atop m}\right\} = \frac {1}{m!} \sum _{i=0}^{m}(-1)^{i}{\binom {m}{i}}(m-i)^{n} \implies \left\{{8 \atop 3}\right\} = 966. \end{eqnarray}$
using Combinatorics
stirlings2(n,m)
966
Even for moderately size problem instances, the number of possible solutions makes any type of brute force approach infeasible:
stirlings2(36,6)
5204262846887190720
A simple approach to randomly generating an allocation array is to randomly generate integers between 1 and m:
using Random
Random.seed!(234215)
a = 0
done = false
while !done
a = rand(1:m,n)
if maximum(a) == m
done = true
end
end
a
8-element Vector{Int64}:
1
2
1
1
2
3
1
1
In the above, if the maximum index in the allocation array does not include m, then the array is regenerated. This is to ensure that only the array is needed to create m resources. As a result, only a is needed to create an m-element array T listing the times of the tasks allocated to each resource (any unallocated resources less than m are not a problem since they would just result in an empty task-index array being created in T):
T = [[] for i = 1:maximum(a)]
for i = 1:length(a)
push!(T[a[i]],t[i])
end
T
3-element Vector{Vector{Any}}:
[1, 3, 4, 7, 8]
[2, 5]
[6]
The total time on each resource can then be determined. Since the makespan will be used in the improvement heuristics, it can be defined as the one-line function makespan that determines the makespan directly from the allocation vector without the creation of the m-element array T. In makespan, the maximum element in the allocation array a is used as the upper limit to search for matches:
[sum(i) for i in T]
3-element Vector{Int64}:
23
7
6
makespan(a,t) = maximum([sum(t[a.==i]) for i=1:maximum(a)])
makespan(a,t)
23
The following function is used to display the total time concatenated with the task times allocated to each resource:
function displaytasktimes(a,t)
T = [Int64[] for i = 1:maximum(a)]
for i = 1:length(a)
push!(T[a[i]],t[i])
end
for i in T
println(sum(i),": ",i)
end
end
displaytasktimes(a,t)
23: [1, 3, 4, 7, 8] 7: [2, 5] 6: [6]
A simple heuristic that can be applied to improve the solution is to relocate a task from one resource to another if it results in a decrease in makespan (e.g., moving any task from resource 1 to resource 3 reduces the makespan from 16 to 15); continuing to consider relocations until there is no improvement and makespan. This is the idea of the following relocate improvement heuristics:
function relocate(aᵒ,t)
cnt = 1
n = length(aᵒ)
m = maximum(aᵒ)
TCᵒ = makespan(aᵒ,t)
done = false
while !done
done = true # Assume no improvement will be found
a₀ = copy(aᵒ) # Copy as current in case improvement found
for i = 1:n, j = 1:m
if a₀[i] != j # Don't relocate if task already allocated to resource j
cnt += 1
a = copy(a₀) # Work on copy of current assignment
a[i] = j # Relocate task i to resource j
TC = makespan(a,t) # Calc. TC of relocation
if TC < TCᵒ # Improvement found
done = false
TCᵒ, aᵒ = copy(TC), copy(a)
end
end
end
end
return (TC = TCᵒ, a = aᵒ, cnt = cnt)
end
relocate (generic function with 1 method)
(TC,a) = relocate(a,t)
(TC = 15, a = [1, 2, 1, 1, 2, 3, 1, 2], cnt = 33)
displaytasktimes(a,t)
15: [1, 3, 4, 7] 15: [2, 5, 8] 6: [6]
Another simple improvement procedure is to consider pairwise interchanges, similar to what was done for the QAP:
function sdpisched(aᵒ,t)
cnt = 1
n = length(aᵒ)
TCᵒ = makespan(aᵒ,t)
done = false
while !done
done = true # Assume no improvement will be found
a₀ = copy(aᵒ) # Copy as current in case improvement found
for i = 1:n-1, j = i+1:n
if a₀[i] != a₀[j]
cnt += 1
a = copy(a₀) # Work on copy of current assignment
a[[j,i]] .= a[[i,j]] # Pairwise interchange of i and j
TC = makespan(a,t) # Calc. TC of interchange
if TC < TCᵒ # Improvement found
done = false
TCᵒ, aᵒ = copy(TC), copy(a)
end
end
end
end
return (TC = TCᵒ, a = aᵒ, cnt = cnt)
end
sdpisched (generic function with 1 method)
(TC,a) = sdpisched(a,t)
(TC = 15, a = [1, 2, 1, 1, 2, 3, 1, 2], cnt = 20)
Not unsurprisingly, in this case, there was no improvement with pairwise interchanges. This was due to their being two resources with total times equal to the makespan; in this case, there is no way of reducing makespan because an interchange between these two resources would, at best, not change and, more likely, increase the makespan. Interchanges between only one of these two resources and any other resource could reduce its total time but would not reduce the makespan because the total time on the other resource would not change.
Constructing a different random solution and then applying the relocate heuristic finds, by chance, the global optimal solution:
Random.seed!(34534)
a = rand(1:m,n)
displaytasktimes(a,t)
7: [1, 6] 10: [3, 7] 19: [2, 4, 5, 8]
(TC,a) = relocate(a,t)
(TC = 12, a = [1, 2, 2, 3, 1, 1, 2, 3], cnt = 49)
displaytasktimes(a,t)
12: [1, 5, 6] 12: [2, 3, 7] 12: [4, 8]
Another construction procedure is to try to allocate an equal number of tasks to each resource in order to avoid, as was seen in just randomly allocating tasks, having one resource with significantly larger numbers of tasks allocated to it. First, the tasks are randomly reordered:
using Random
Random.seed!(48329)
idx = randperm(n)
8-element Vector{Int64}:
3
7
2
6
8
5
1
4
Then, a sequence of m task is then allocated to each resource:
a = zeros(n)
step = Int(ceil(n/m)) # Tasks per resource
for i = 1:m
a[idx[(i-1)*step+1:min(i*step,n)]] .= i
end
a = Int.(a)
8-element Vector{Int64}:
3
1
1
3
2
2
1
2
In order to see how the indices of each sequence is determined, the following code just outputs the indices (note: min(i*step,n) is used so that the final ending index does not exceed n):
for i = 1:m
println([(i-1)*step+1:min(i*step,n)])
end
UnitRange{Int64}[1:3]
UnitRange{Int64}[4:6]
UnitRange{Int64}[7:8]
displaytasktimes(a,t)
12: [2, 3, 7] 19: [5, 6, 8] 5: [1, 4]
The relocation and pairwise improvement heuristics are then applied to the constructed solution:
(TC,a) = relocate(a,t)
(TC = 13, a = [3, 1, 1, 3, 2, 3, 1, 2], cnt = 33)
displaytasktimes(a,t)
12: [2, 3, 7] 13: [5, 8] 11: [1, 4, 6]
(TC,a) = sdpisched(a,t)
(TC = 12, a = [3, 1, 1, 2, 3, 3, 1, 2], cnt = 43)
A repair shop has four identical stations available. Today, 10, 7, 4, 2, 9, 8, 4, and 10 repairs will be performed for each of eight different types of repair, where each type takes, on average, 53, 20, 53, 15, 30, 55, 32, and 21 minutes, respectively. Determine a lower bound on the maximum time that it will take to complete all of the repairs, assuming that it takes five minutes to switch between each repair at the station no matter the type of repair or the type of the preceding repair.
# Create data
f = [10, 7, 4, 2, 9, 8, 4, 10] # Number of each repair type
toper = [53, 20, 53, 15, 30, 55, 32, 21] # Time (min) of each repair type
t0 = vcat([repeat([i],j) for (i,j) in zip(toper,f)]...) # Time (min) of each repair
54-element Vector{Int64}:
53
53
53
53
53
53
53
53
53
53
20
20
20
⋮
32
32
21
21
21
21
21
21
21
21
21
21
t = t0 .+ 5 # Repair and switch time (min)
m = 4 # Number of identical stations
LB = ceil(sum(t)/m) # Round up since all times integer
558.0
MILP formulation of schedule independent tasks on parallel identical resources without preemption to minimize makespan problem:
$ \begin{eqnarray} \quad \mbox{Minimize} \quad y \\ \quad \mbox{subject to} \quad y &\ge& \sum_{i \in N} t_i x_{ij}, &\quad j \in M \\ \sum_{j \in M} x_{ij} &=& 1, &\quad i \in N \\ y &\ge& 0 \\ x_{ij} &\in& \bigl\{ 0, 1 \bigr\}, \end{eqnarray} $
where
$ \begin{eqnarray} \quad y &=& \mbox{makespan} \\ x_{ij} &=& \begin{cases} 1, & \mbox{if task } i \mbox{ is allocated to resource } j \\ 0, & \mbox{otherwise} \end{cases} \\ M &=& \bigl\{ 1, \dots, m \bigr\}, \quad \mbox{set of } m \mbox{ resources} \\ N &=& \bigl\{ 1, \dots, n \bigr\}, \quad \mbox{set of } n \mbox{ tasks} \\ \end{eqnarray} $