from dataclasses import dataclass
@dataclass
class Facilities:
name: str
lengths: list
names: list
traffic: list
@property
def num_facilities(self):
return len(self.lengths)
@property
def span(self):
return sum(self.lengths)
@property
def total_trips(self):
return sum(w for _, _, w in self.traffic_pairs)
@property
def traffic_pairs(self):
n = self.num_facilities
pairs = ((i, j) for i in range(n) for j in range(i + 1, n))
return [(i, j, self.traffic[i][j]) for i, j in pairs if self.traffic[i][j] > 0]
@dataclass
class Layout:
cost: float
order: list
def listed_order(facilities):
return list(range(facilities.num_facilities))
def order_from_left_counts(n, is_left):
counts = [(sum(j != i and is_left(i, j) for j in range(n)), i) for i in range(n)]
counts.sort(reverse=True)
return [i for _, i in counts]Single-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.
import random
import string
from plot import draw_cost_distribution, draw_layout
from report import report_comparison, report_order
def 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],
],
)
def load_random(num_facilities=40, seed=4):
draw = random.Random(seed)
alphabet = string.ascii_uppercase + string.digits
indices = range(num_facilities)
lengths = [draw.randint(10, 30) for _ in indices]
names = ["".join(draw.choice(alphabet) for _ in range(6)) for _ in indices]
traffic = [[draw.randint(1, 15) if j > i else 0 for j in indices] for i in indices]
return Facilities(f"{num_facilities} random facilities", lengths, names, traffic)A 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.
import knitrodef minimize_cost(facilities, *, multistart=True, max_nodes=2, num_threads=1):
n = facilities.num_facilities
lengths, traffic = facilities.lengths, facilities.traffic
pairs = [(i, j) for i in range(n) for j in range(n) if i != j]
ordered = [(i, j) for i in range(n) for j in range(i + 1, n)]
prob = knitro.Problem()
a = {pair: prob.add_variable(vtype=knitro.KN_VARTYPE_BINARY) for pair in pairs}
def between(k, i, j):
return a[i, k] * a[k, j] + a[j, k] * a[k, i]
def distance(i, j):
others = [k for k in range(n) if k != i and k != j]
spanned = prob.nsum(lengths[k] * between(k, i, j) for k in others)
return spanned + (lengths[i] + lengths[j]) / 2
prob.add_objective(prob.nsum(traffic[i][j] * distance(i, j) for i, j in ordered))
for i, j in ordered:
prob.add_constraint(a[i, j] + a[j, i] == 1)
for i, j in ordered:
for k in range(j + 1, n):
prob.add_constraint(a[i, j] + a[j, k] + a[k, i] <= 2)
prob.add_constraint(a[i, k] + a[k, j] + a[j, i] <= 2)
prob.set_param(knitro.KN_PARAM_MIP_MULTISTART, int(multistart))
prob.set_param(knitro.KN_PARAM_MIP_MAXNODES, max_nodes)
prob.set_param(knitro.KN_PARAM_NUMTHREADS, num_threads)
prob.solve()
order = order_from_left_counts(n, lambda i, j: a[i, j].value >= 0.99)
cost = prob.get_attr(knitro.KN_ATTR_OBJ_VALUE)
return Layout(cost, order)Solving 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.
concurrent_evals 0
feastol 1e-06
feastol_abs 1e-06
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 expressions: 425 | 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.134
1 1644.75 0.00000e+00 7.60894 0.135
2 1420.53 2.22045e-16 13.0091 0.136
3 1216.91 2.22045e-16 0.985170 0.136
4 1034.70 2.22045e-16 0.270487 0.137
5 1022.60 1.11022e-16 1.06398e-03 0.137
6 1022.50 2.22045e-16 1.13268e-06 0.137
7 1022.50 4.18332e-13 6.06552e-09 0.138
8 1022.50 4.18332e-13 6.06552e-09 0.308
Tree search
-----------
Nodes Best solution Best bound Gap Time
Expl | Unexpl value value (secs)
--------------- ------------- ---------- --- ------
1 1 1022.50 LEAF -inf 0.309
2 1 1022.50 -inf 0.313
EXIT: Node limit reached. Integer feasible point found.
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.295s)
# 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.295s)
Total program time (secs) = 0.31347 (0.025 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.
concurrent_evals 0
feastol 1e-06
feastol_abs 1e-06
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 expressions: 162758 | 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.33856e+06 0.166769 8.60951 0.354
1 2.27532e+06 0.101940 20.9303 0.704
2 2.21369e+06 6.80025e-02 16.6109 1.075
3 2.15064e+06 3.52685e-02 19.8184 1.415
4 2.06546e+06 1.49496e-02 21.7320 1.780
5 1.97448e+06 3.87244e-03 28.0237 2.146
6 1.91181e+06 3.51283e-04 28.8892 2.496
7 1.82405e+06 2.22045e-16 13.1467 2.914
8 1.75933e+06 1.97620e-14 11.5862 3.136
9 1.64065e+06 1.35856e-07 14.4619 3.345
10 1.54201e+06 1.99510e-07 12.1321 3.558
11 1.44033e+06 1.21514e-07 9.90629 3.792
12 1.40425e+06 7.87751e-08 8.13230 4.000
13 1.39581e+06 7.16980e-08 6.71209 4.196
14 1.39300e+06 5.97955e-08 5.89443 4.323
15 1.38581e+06 6.21385e-08 5.97658 4.508
16 1.38329e+06 9.85575e-08 4.38112 4.701
17 1.38164e+06 9.29367e-08 4.95949 4.876
18 1.37950e+06 6.96680e-08 5.66136 5.068
19 1.37860e+06 5.74787e-08 4.65371 5.244
20 1.37762e+06 8.17737e-08 5.35671 5.437
30 1.37118e+06 5.25729e-08 3.70539 7.308
40 1.36277e+06 2.37942e-08 7.18900 11.397
50 1.35912e+06 2.56839e-08 4.12542 15.369
60 1.35481e+06 2.57526e-08 5.58718 19.036
70 1.35388e+06 4.44089e-16 6.54581 22.044
80 1.35094e+06 6.66134e-16 6.42211 24.554
90 1.34710e+06 8.88178e-16 5.28164 28.087
100 1.34379e+06 8.42589e-08 2.77818 30.727
110 1.34278e+06 3.11503e-08 4.34165 34.839
120 1.34140e+06 5.99691e-08 6.65781 37.722
130 1.34039e+06 3.32582e-08 2.14572 40.447
140 1.33988e+06 1.27143e-07 2.61922 42.975
150 1.33923e+06 8.17746e-08 4.47694 46.309
160 1.33874e+06 1.35420e-10 2.56859e-09 48.907
Tree search
-----------
Nodes Best solution Best bound Gap Time
Expl | Unexpl value value (secs)
--------------- ------------- ---------- --- ------
1 2 1.33874e+06 FCRD -inf 50.339
2 3 1.33843e+06 FCRD -inf 70.440
EXIT: Node limit reached. Integer feasible point found.
Final Statistics for MIP
------------------------
Final objective value = 1.33843050000000000e+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 (57.350s)
# 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 (66.447s)
Total program time (secs) = 70.44053 (69.392 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.012s
Rounding heuristic = 2 / 2 / 0.047s
MPEC heuristic = 1 / 0 / 9.103s
Local search heuristic = 6 / 0 / 2.282s
===========================================================================
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.