struct Facilities
name::String
lengths::Vector{Int}
names::Vector{String}
traffic::Matrix{Int}
end
num_facilities(facilities::Facilities) = length(facilities.lengths)
span(facilities::Facilities) = sum(facilities.lengths)
function traffic_pairs(facilities::Facilities)
n = num_facilities(facilities)
out = Tuple{Int,Int,Int}[]
for i in 1:n, j in (i + 1):n
facilities.traffic[i, j] == 0 && continue
push!(out, (i, j, facilities.traffic[i, j]))
end
return out
end
total_trips(facilities::Facilities) = sum(w for (_, _, w) in traffic_pairs(facilities))
struct Layout
cost::Float64
order::Vector{Int}
end
listed_order(facilities::Facilities) = collect(1:num_facilities(facilities))
function order_from_left_counts(n, is_left)
counts = [(count(j -> j != i && is_left(i, j), 1:n), i) for i in 1:n]
sort!(counts; rev=true)
return [i for (_, i) in counts]
endSingle-Row Facility Layout
Order facilities of given lengths along a line to minimize the weighted traffic between them, a nonconvex binary quadratic program over the pairwise ordering.
Introduction
We consider the problem of finding an arrangement of elements of given lengths (facilities, rooms, machines…) on a line so as to minimize the weighted sum of the distances between all the pairs of elements, that is, minimize the traffic intensity.
Thus, the problem can be reduced to finding an ordering.
Here is an illustration of the 3 different (non-symmetric) solutions of an example problem with 3 facilities.
Each facility is represented by a colored rectangle, with the length of the rectangle corresponding to the facility’s size. The rectangles above them represent the traffic. The width of each one corresponds to the distance between the centers of the corresponding facilities, and its height corresponds to the traffic intensity. Thus, the objective value corresponds to the sum of their areas. That total is carried across on the right and repacked into a single box, so that the height it fills is the cost of that ordering.
In this example, the best solution is the third one. Scissor Lift has high traffic. Therefore, putting it in the middle minimizes the overall traffic.
Problem description
Input
- A set of \(N\) facilities. For each facility \(i = 1, \dots, N\), its length \(L_i \in \mathbb{R}^+\)
- The traffic density between each pair of facilities \(C_{i,j} \in \mathbb{R}^+\), \(i = 1, \dots, N\), \(j = i + 1, \dots, N\)
Problem: find a permutation of facilities
Objective: minimize the traffic intensity
\[ \sum_{i = 1}^N \sum_{j = i + 1}^N C_{i, j}\, l_{i, j} \]
where \(l_{i, j}\) is the distance between facility \(i\) and facility \(j\) in the solution.
Mixed-integer nonlinear program
There are various possible models for the problem. We propose here one relying on boolean variables indicating the relative positions between each pair of facilities.
Variables
- \(a_{i, j} \in \{0, 1\}\), \(i = 1, \dots, N\), \(j = 1, \dots, N\): \(a_{i, j} = 1\) if facility \(i\) is to the left of facility \(j\)
Objective
\[ \min \sum_{i=1}^{N-1} \sum_{j=i+1}^N C_{i, j} \left( \sum_{\substack{k=1 \\ k \neq i,\, k \neq j}}^{N} L_k \big(a_{i, k} a_{k, j} + a_{j, k} a_{k, i}\big) + \frac{1}{2}(L_i + L_j) \right) \]
To understand this formula, let’s note that \(a_{i, k} a_{k, j} + a_{j, k} a_{k, i} = 1\) if facility \(k\) is located between facility \(i\) and facility \(j\). In this case, a cost \(C_{i, j} L_k\) is paid. Moreover, the distance between 2 machines is from center to center, hence the term \(\frac{1}{2}(L_i + L_j)\).
Constraints
- Either facility \(i\) is to the left of facility \(j\), or facility \(j\) is to the left of facility \(i\):
\[ \forall i = 1, \dots, N-1, \quad \forall j = i+1, \dots, N \qquad a_{i, j} + a_{j, i} = 1 \]
- If facility \(i\) is located before facility \(j\) and facility \(j\) is located before facility \(k\), then facility \(i\) must be located before facility \(k\). Equivalently, no three facilities may form a cycle:
\[ \forall i, j, k \text{ distinct} \qquad a_{i, j} + a_{j, k} + a_{k, i} \leq 2 \]
The model has the following properties:
- Binary variables only
- Non-convex quadratic objective and linear constraints only
- Therefore, it is a non-convex BQP (Binary Quadratic Problem)
The number of variables is quadratic in the number of facilities. It is possible to linearize this model by introducing three-index variables \(a_{i, j, k}\). However, the number of variables then becomes cubic in the number of facilities and, even if it is linear, the model quickly becomes intractable when the number of facilities increases.
Input data
An instance gives a length for every facility and a traffic[i][j] count of the trips per period between facilities i and j. Only the entries above the diagonal are filled in, since a trip counts the same in either direction. A Layout holds the total cost and the facilities in left-to-right order; the model works in pairwise variables instead, so order_from_left_counts turns them back into that order.
Two instances are used below, a six-machine workshop written out in full and a larger one drawn at random. Each language draws with its own generator, so the Python and Julia pages lay out different facilities and their costs are not comparable side by side.
using Random
function load_workshop()
return Facilities(
"Workshop",
[6, 12, 2, 7, 13, 6],
["Lathe", "Drill", "Mill", "Punch press", "Centreless grinder", "Shaper"],
[
0 5 1 10 4 19
0 0 9 10 7 4
0 0 0 6 2 7
0 0 0 0 1 4
0 0 0 0 0 0
0 0 0 0 0 0
],
)
end
function load_random(; num_facilities=40, seed=4)
rng = MersenneTwister(seed)
alphabet = ['A':'Z'; '0':'9']
lengths = [rand(rng, 10:30) for _ in 1:num_facilities]
names = [String(rand(rng, alphabet, 6)) for _ in 1:num_facilities]
traffic = zeros(Int, num_facilities, num_facilities)
for i in 1:num_facilities, j in (i + 1):num_facilities
traffic[i, j] = rand(rng, 1:15)
end
return Facilities("$num_facilities random facilities", lengths, names, traffic)
endA layout figure carries one rectangle per pair of facilities, as wide as the distance between them and as tall as the trips they exchange, so its area is what that pair costs and the whole stack is the total. The box on the right repacks those same areas into one column, and that is what makes two layouts comparable at a glance.
The first instance is a workshop of six machines, drawn in the order the data lists them.
facilities = load_workshop()draw_layout(facilities, listed_order(facilities), "Machines as listed")Every pair has its own color, and the total traffic does not depend on the ordering, so the stack is always the same height whatever the layout. Only the widths change, so a good layout is one where the tall rectangles are narrow. The dashed rule on the repacked box marks the cost, and that box is drawn to a fixed size, so its fill can be read against any other layout of the same instance.
Model implementation
minimize_cost(facilities) builds the model and hands it to Knitro’s MIP multistart.
using JuMP
using KNITROfunction minimize_cost(facilities; multistart=true, max_nodes=2, num_threads=1)
n = num_facilities(facilities)
lengths, traffic = facilities.lengths, facilities.traffic
ordered = [(i, j) for i in 1:n for j in (i + 1):n]
model = Model(KNITRO.Optimizer)
set_attribute(model, "mip_multistart", Int(multistart))
set_attribute(model, "mip_maxnodes", max_nodes)
set_attribute(model, "numthreads", num_threads)
@variable(model, a[i in 1:n, j in 1:n; i != j], Bin)
between(k, i, j) = a[i, k] * a[k, j] + a[j, k] * a[k, i]
function distance(i, j)
spanned = sum(lengths[k] * between(k, i, j) for k in 1:n if k != i && k != j)
return spanned + (lengths[i] + lengths[j]) / 2
end
@objective(model, Min, sum(traffic[i, j] * distance(i, j) for (i, j) in ordered))
for (i, j) in ordered
@constraint(model, a[i, j] + a[j, i] == 1)
end
for (i, j) in ordered, k in (j + 1):n
@constraint(model, a[i, j] + a[j, k] + a[k, i] <= 2)
@constraint(model, a[i, k] + a[k, j] + a[j, i] <= 2)
end
optimize!(model)
order = order_from_left_counts(n, (i, j) -> value(a[i, j]) >= 0.99)
cost = objective_value(model)
return Layout(cost, order)
endSolving the workshop
The model reorders the same six machines. The traffic is unchanged, so the stack keeps its height; the whole gain is in the widths.
layout = minimize_cost(facilities);=======================================
Commercial License
Artelys Knitro 16.0.0
=======================================
Knitro changing mip_method from AUTO to 1.
No start point provided -- Knitro computing one.
Knitro presolve eliminated 0 variables (0%) and 0 constraints (0%) in 0.00s.
datacheck 0
feastol 1e-06
feastol_abs 1e-06
hessian_no_f 1
mip_maxnodes 2
mip_multistart 1
numthreads 1
opttol 1e-06
opttol_abs 0.001
Knitro changing mip_root_nlpalg from AUTO to 1.
Knitro changing mip_node_nlpalg from AUTO to 1.
Knitro changing mip_branchrule from AUTO to 2.
Knitro changing mip_selectrule from AUTO to 2.
Knitro changing mip_mir from AUTO to 2.
Knitro changing mip_clique from AUTO to 0.
Knitro changing mip_zerohalf from AUTO to 0.
Knitro changing mip_liftproject from AUTO to 0.
Knitro changing mip_knapsack from AUTO to 2.
Knitro changing mip_gomory from AUTO to 0.
Knitro changing mip_cut_flowcover from AUTO to 2.
Knitro changing mip_cut_probing from AUTO to 1.
Knitro changing mip_rounding from AUTO to 3.
Knitro changing mip_heuristic_strategy from AUTO to 1.
Knitro changing mip_heuristic_feaspump from AUTO to 1.
Knitro changing mip_heuristic_misqp from AUTO to 0.
Knitro changing mip_heuristic_mpec from AUTO to 1.
Knitro changing mip_heuristic_diving from AUTO to 1926.
Knitro changing mip_heuristic_fixpropagate from AUTO to 62.
Knitro changing mip_heuristic_lns from AUTO to 0.
Knitro changing mip_heuristic_localsearch from AUTO to 1.
Knitro changing mip_pseudoinit from AUTO to 1.
Problem Characteristics | Presolved
-----------------------
Problem type: MIQP
Objective: minimize / quadratic
Number of variables: 30 | 30
bounds: lower upper range | lower upper range
0 0 30 | 0 0 30
free fixed | free fixed
0 0 | 0 0
cont. binary integer | cont. binary integer
0 30 0 | 0 30 0
Number of constraints: 55 | 55
eq. ineq. range | eq. ineq. range
linear: 15 40 0 | 15 40 0
quadratic: 0 0 0 | 0 0 0
Number of nonzeros:
objective Jacobian Hessian | objective Jacobian Hessian
linear: 0 150 | 0 150
quadratic: 30 0 112 | 30 0 112
total: 30 150 112 | 30 150 112
Knitro using Branch and Bound method with 1 thread.
Initial points
--------------
No initial point provided for the root node relaxation.
No primal point provided for the MIP.
Coefficient range:
linear objective: [0e+00, 0e+00] | [0e+00, 0e+00]
linear constraints: [1e+00, 1e+00] | [1e+00, 1e+00]
quadratic objective: [2e+00, 2e+02] | [6e-01, 7e+01]
quadratic constraints: [0e+00, 0e+00] | [0e+00, 0e+00]
variable bounds: [1e+00, 1e+00] | [1e+00, 1e+00]
constraint bounds: [1e+00, 2e+00] | [1e+00, 2e+00]
Root node relaxation
--------------------
Iter Objective Feasibility Optimality Time
error error (secs)
---- --------- ----------- ---------- ------
0 2034.05 0.106794 15.0701 0.012
1 1644.75 0.00000e+00 7.60894 0.013
2 1420.53 2.22045e-16 13.0091 0.013
3 1216.91 2.22045e-16 0.985170 0.014
4 1034.70 2.22045e-16 0.270487 0.014
5 1022.60 1.11022e-16 1.06398e-03 0.014
6 1022.50 2.22045e-16 1.13268e-06 0.014
7 1022.50 4.18332e-13 6.06552e-09 0.014
8 1022.50 4.18332e-13 6.06552e-09 0.015
Tree search
-----------
Nodes Best solution Best bound Gap Time
Expl | Unexpl value value (secs)
--------------- ------------- ---------- --- ------
1 1 1022.50 LEAF -inf 0.015
2 1 1022.50 -inf 0.019
EXIT: Node limit reached. Integer feasible point found.
HINT: Knitro spent 1.5% of solution time (0.000283 secs) checking model
convexity. To skip the automatic convexity checker for QPs and QCQPs,
explicity set the user option convex=0 or convex=1.
Final Statistics for MIP
------------------------
Final objective value = 1.02250000054890518e+03
Final bound value = -inf
Final optimality gap (abs / rel) = inf / inf
# of root cutting plane rounds = 0
# of restarts = 0
# of nodes processed = 2 (0.012s)
# of strong branching evaluations = 0 (0.000s)
# of function evaluations = 0 (0.000s)
# of gradient evaluations = 0 (0.000s)
# of hessian evaluations = 0 (0.000s)
# of hessian-vector evaluations = 0
# of subproblems processed = 2 (0.012s)
Total program time (secs) = 0.01892 (0.016 CPU time)
Time spent in evaluations (secs) = 0.00000
Cuts statistics (gen / add)
---------------------------
Knapsack cuts = 0 / 0
Mixed-integer rounding cuts = 0 / 0
Flow-cover cuts = 0 / 0
Probing cuts = 0 / 0
Heuristics statistics (calls / successes / time)
------------------------------------------------
Feasibility pump = 0 / 0 / 0.000s
Rounding heuristic = 0 / 0 / 0.000s
MPEC heuristic = 0 / 0 / 0.000s
Local search heuristic = 0 / 0 / 0.000s
===========================================================================
Show the full outputHide the full output
draw_layout(facilities, layout.order, "Knitro solution")report_order(facilities, layout) # facility length centre
----------------------------------------
1 Centreless grinder 13 6.5
2 Drill 12 19.0
3 Mill 2 26.0
4 Punch press 7 30.5
5 Lathe 6 37.0
6 Shaper 6 43.0
Total cost 1,022One pair dominates the cost. Moving it from the ends of the row to the middle is most of the gain:
report_comparison(facilities, listed_order(facilities), layout.order)Busiest pair: Lathe and Shaper, 19 trips
apart pair cost share total
----------------------------------------------
start 40.0 760 40% 1,912
solved 6.0 114 11% 1,022
All 360 orderings checked: the cheapest costs 1,022This instance is small enough to enumerate every ordering, so the last line is a proof rather than a claim: Knitro returns the cheapest layout that exists.
A larger instance
The number of binary variables grows quadratically with the facilities, so a larger instance is a harder problem.
facilities_large = load_random(; num_facilities=40)layout_large = minimize_cost(facilities_large);=======================================
Commercial License
Artelys Knitro 16.0.0
=======================================
Knitro changing mip_method from AUTO to 1.
No start point provided -- Knitro computing one.
Knitro presolve eliminated 0 variables (0%) and 0 constraints (0%) in 0.03s.
datacheck 0
feastol 1e-06
feastol_abs 1e-06
hessian_no_f 1
mip_maxnodes 2
mip_multistart 1
numthreads 1
opttol 1e-06
opttol_abs 0.001
Knitro changing mip_root_nlpalg from AUTO to 1.
Knitro changing mip_node_nlpalg from AUTO to 1.
Knitro changing mip_branchrule from AUTO to 2.
Knitro changing mip_selectrule from AUTO to 2.
Knitro changing mip_mir from AUTO to 2.
Knitro changing mip_clique from AUTO to 0.
Knitro changing mip_zerohalf from AUTO to 0.
Knitro changing mip_liftproject from AUTO to 0.
Knitro changing mip_knapsack from AUTO to 2.
Knitro changing mip_gomory from AUTO to 0.
Knitro changing mip_cut_flowcover from AUTO to 2.
Knitro changing mip_cut_probing from AUTO to 1.
Knitro changing mip_rounding from AUTO to 3.
Knitro changing mip_heuristic_strategy from AUTO to 1.
Knitro changing mip_heuristic_feaspump from AUTO to 1.
Knitro changing mip_heuristic_misqp from AUTO to 0.
Knitro changing mip_heuristic_mpec from AUTO to 1.
Knitro changing mip_heuristic_diving from AUTO to 1926.
Knitro changing mip_heuristic_fixpropagate from AUTO to 62.
Knitro changing mip_heuristic_lns from AUTO to 0.
Knitro changing mip_heuristic_localsearch from AUTO to 1.
Knitro changing mip_pseudoinit from AUTO to 1.
Problem Characteristics | Presolved
-----------------------
Problem type: MIQP
Objective: minimize / quadratic
Number of variables: 1560 | 1560
bounds: lower upper range | lower upper range
0 0 1560 | 0 0 1560
free fixed | free fixed
0 0 | 0 0
cont. binary integer | cont. binary integer
0 1560 0 | 0 1560 0
Number of constraints: 20540 | 20540
eq. ineq. range | eq. ineq. range
linear: 780 19760 0 | 780 19760 0
quadratic: 0 0 0 | 0 0 0
Number of nonzeros:
objective Jacobian Hessian | objective Jacobian Hessian
linear: 0 60840 | 0 60840
quadratic: 1560 0 59280 | 1560 0 59280
total: 1560 60840 59280 | 1560 60840 59280
Knitro using Branch and Bound method with 1 thread.
Initial points
--------------
No initial point provided for the root node relaxation.
No primal point provided for the MIP.
Coefficient range:
linear objective: [0e+00, 0e+00] | [0e+00, 0e+00]
linear constraints: [1e+00, 1e+00] | [1e+00, 1e+00]
quadratic objective: [1e+01, 4e+02] | [1e-01, 5e+00]
quadratic constraints: [0e+00, 0e+00] | [0e+00, 0e+00]
variable bounds: [1e+00, 1e+00] | [1e+00, 1e+00]
constraint bounds: [1e+00, 2e+00] | [1e+00, 2e+00]
Root node relaxation
--------------------
Iter Objective Feasibility Optimality Time
error error (secs)
---- --------- ----------- ---------- ------
0 2.53658e+06 0.166769 5.65160 0.310
1 2.46239e+06 0.101940 22.7792 0.593
2 2.39141e+06 6.85958e-02 25.1507 0.852
3 2.30895e+06 4.41083e-02 28.1867 1.189
4 2.23220e+06 2.22195e-02 35.7614 1.538
5 2.14139e+06 1.10704e-02 24.9259 1.885
6 2.04018e+06 5.31888e-03 18.3201 2.251
7 1.93055e+06 2.40186e-03 11.6622 2.614
8 1.86032e+06 8.15117e-04 13.6948 2.983
9 1.68914e+06 1.28082e-05 10.8198 3.211
10 1.60070e+06 4.49224e-06 13.8748 3.326
11 1.54559e+06 1.33550e-06 10.2525 3.491
12 1.52955e+06 9.69071e-07 9.36728 3.596
13 1.51560e+06 5.72264e-07 7.79460 3.755
14 1.50641e+06 3.27611e-07 6.89730 3.925
15 1.49855e+06 1.57133e-07 5.27955 4.089
16 1.49500e+06 8.85535e-08 4.13449 4.254
17 1.49408e+06 8.54033e-08 4.16849 4.368
18 1.49255e+06 9.49949e-08 5.27102 4.559
19 1.49245e+06 9.63298e-08 5.25282 4.744
20 1.49204e+06 1.05469e-07 3.80944 4.928
30 1.48636e+06 2.68194e-08 0.814890 7.153
40 1.48480e+06 1.65480e-07 5.73967 12.254
50 1.48446e+06 1.10379e-02 1.11885 16.685
60 1.48260e+06 5.93381e-05 4.29542 20.688
70 1.48024e+06 4.79342e-06 9.62325 26.092
80 1.47960e+06 4.74954e-02 11.0461 26.957
90 1.47957e+06 4.26961e-06 3.04204e-04 27.972
Tree search
-----------
Nodes Best solution Best bound Gap Time
Expl | Unexpl value value (secs)
--------------- ------------- ---------- --- ------
1 2 1.47957e+06 FCRD -inf 28.365
Knitro deduced that the problem is non-convex.
2 3 1.47622e+06 FCRD -inf 43.823
EXIT: Node limit reached. Integer feasible point found.
Final Statistics for MIP
------------------------
Final objective value = 1.47622050000000000e+06
Final bound value = -inf
Final optimality gap (abs / rel) = inf / inf
# of root cutting plane rounds = 1
# of restarts = 0
# of nodes processed = 2 (34.798s)
# of strong branching evaluations = 0 (0.000s)
# of function evaluations = 0 (0.000s)
# of gradient evaluations = 0 (0.000s)
# of hessian evaluations = 0 (0.000s)
# of hessian-vector evaluations = 0
# of subproblems processed = 6 (40.765s)
Total program time (secs) = 43.82273 (43.354 CPU time)
Time spent in evaluations (secs) = 0.00000
Cuts statistics (gen / add)
---------------------------
Knapsack cuts = 0 / 0
Mixed-integer rounding cuts = 0 / 0
Flow-cover cuts = 0 / 0
Probing cuts = 1560 / 0
Heuristics statistics (calls / successes / time)
------------------------------------------------
Feasibility pump = 1 / 0 / 0.016s
Rounding heuristic = 2 / 2 / 0.039s
MPEC heuristic = 1 / 0 / 5.975s
Local search heuristic = 6 / 0 / 2.468s
===========================================================================
Show the full outputHide the full output
Every pair of machines carries traffic here, so one rectangle per pair would draw a solid block. The figure below compares the solved cost against orderings drawn at random instead.
draw_cost_distribution(facilities_large, layout_large)The node limit stops the search before it produces a bound, so the reported gap stays infinite and this layout is not proven optimal. It is still far cheaper than anything chance produced. The MIP multistart heuristic finds it in two nodes.