5.1 Facility Location Models: Center of Gravity & Minisum / Minimax Distance

Facility Location Models: Center of Gravity & Minisum / Minimax Distance

Facility location decisions represent some of the most critical long-term capital commitments in industrial engineering and supply chain systems. Siting a new plant, regional distribution center (DC), cross-dock terminal, or emergency response hub locks in fixed operating economics, freight transportation costs, labor accessibility, and customer service delivery times for decades. A sub-optimal layout within an existing building can be re-engineered with modest capital, but an improperly located facility imposes recurring transportation penalties that cannot be overcome by operational improvements.

Industrial engineers classify facility location problems across several fundamental dimensions:

• Single-Facility vs. Multi-Facility Models: Single-facility models determine the optimal spatial coordinates $(x, y)$ for one new facility interacting with $m$ existing fixed nodes (suppliers, customer clusters, ports). Multi-facility models determine the locations of $n$ new facilities simultaneously ($n \ge 2$), accounting not only for flows between the new facilities and existing demand points, but also for significant inter-facility transfers among the new facilities themselves.
• Continuous Space (Planar) vs. Discrete Network Models: Continuous location models allow the facility to be positioned anywhere on a continuous two-dimensional Euclidean or rectilinear coordinate plane $\mathbb{R}^2$. Discrete models restrict candidate locations to a finite set of pre-screened geographical parcels or graph network nodes, formulated using mixed-integer linear programming (MILP) such as the Capacitated Facility Location Problem (CFLP).

+-----------------------------------------------------------------------------------+
|                         FACILITY LOCATION PROBLEM TAXONOMY                        |
+-----------------------------------------------------------------------------------+
|  Single-Facility Continuous     |  Multi-Facility Continuous   |  Discrete Network |
|  - Center of Gravity (Centroid) |  - Multi-Weber Problem       |  - P-Median       |
|  - Rectilinear Minisum (Median) |  - Coupled Iterative Methods |  - P-Center       |
|  - Euclidean Minisum (Weiszfeld)|  - Allocation-Location Heur. |  - Capacitated FLP|
|  - Minimax (Emergency Siting)   |  - Transportation-Location   |  - Set Covering   |
+-----------------------------------------------------------------------------------+

---

The Center of Gravity (Centroid) Method

The Center of Gravity (Centroid) Method is the most widely utilized continuous single-facility location heuristic. It determines the spatial coordinates of a central facility (such as a regional distribution warehouse) that balances the inbound shipment volumes from suppliers and the outbound shipment volumes to demand nodes.

Mathematical Formulation

Let $m$ be the number of existing facilities (supply points or customer market zones). Each existing facility $i$ (for $i = 1, 2, \dots, m$) is defined by its known Cartesian coordinates $(x_i, y_i)$ and an associated weight $W_i$. The weight $W_i$ represents the flow intensity, typically expressed as annual tonnage, truckload volume, piece shipments, or the product of freight volume and ton-mile transportation freight rates ($W_i = F_i \times c_i$).

The coordinates of the Center of Gravity $(\bar{x}, \bar{y})$ are computed as the weighted average of the coordinates:

$\bar{x} = \frac{\sum_{i=1}^m W_i x_i}{\sum_{i=1}^m W_i}, \qquad \bar{y} = \frac{\sum_{i=1}^m W_i y_i}{\sum_{i=1}^m W_i}$

Physical and Mathematical Derivation

The physical analogy of the Center of Gravity method is finding the balance point (center of mass) of a flat, weightless rigid plate supporting $m$ concentrated vertical point masses $W_i$ positioned at points $(x_i, y_i)$. At the center of mass, the sum of the gravitational moments about any horizontal axis equals zero:

$\sum_{i=1}^m W_i (x_i - \bar{x}) = 0 \implies \bar{x} \sum_{i=1}^m W_i = \sum_{i=1}^m W_i x_i \implies \bar{x} = \frac{\sum_{i=1}^m W_i x_i}{\sum_{i=1}^m W_i}$

From an optimization perspective, the Center of Gravity method is the exact analytical solution to the Squared Euclidean Distance Objective Function. Suppose the objective function minimizes the sum of weighted squared straight-line distances from the unknown location $(x, y)$ to all $m$ existing facilities:

$f(x, y) = \sum_{i=1}^m W_i d_i^2(x, y) = \sum_{i=1}^m W_i \left[ (x - x_i)^2 + (y - y_i)^2 \right]$

Because $f(x, y)$ is strictly convex and continuously differentiable, taking partial derivatives with respect to $x$ and $y$ and setting them to zero yields:

$\frac{\partial f}{\partial x} = 2 \sum_{i=1}^m W_i (x - x_i) = 0 \implies 2x \sum_{i=1}^m W_i - 2 \sum_{i=1}^m W_i x_i = 0 \implies x^* = \frac{\sum_{i=1}^m W_i x_i}{\sum_{i=1}^m W_i}$

$\frac{\partial f}{\partial y} = 2 \sum_{i=1}^m W_i (y - y_i) = 0 \implies 2y \sum_{i=1}^m W_i - 2 \sum_{i=1}^m W_i y_i = 0 \implies y^* = \frac{\sum_{i=1}^m W_i y_i}{\sum_{i=1}^m W_i}$

Critical Engineering Limitation

While the Center of Gravity method provides an outstanding initial heuristic, industrial engineers must recognize its primary mathematical limitation: actual freight transportation costs vary linearly with distance, not with the square of distance ($C \propto d$, not $C \propto d^2$). Squaring the distance penalizes distant outlier nodes disproportionately. For example, a customer located 100 miles away receives 10,000 times the penalty weight of a customer 1 mile away, rather than 100 times. Consequently, the centroid is pulled excessively toward remote, isolated demand nodes.

Step-by-Step Numerical Example: Center of Gravity Calculation

Problem Statement: A heavy machinery manufacturer plans to site a centralized regional parts distribution center to serve four assembly plants across the Midwest. The Cartesian coordinates (in miles, mapped onto a state plane coordinate system) and projected annual parts demand (in full trailer truckloads per year) are summarized below:

Step 1: Compute Total Demand Weight ($\sum W_i$) $\sum_{i=1}^4 W_i = 200 + 100 + 300 + 400 = 1,000\text{ truckloads/year}$

Step 2: Calculate Weighted Coordinates for the X-Axis ($\sum W_i x_i$)

$$\begin{aligned} \sum_{i=1}^4 W_i x_i &= (200 \times 20) + (100 \times 80) + (300 \times 50) + (400 \times 10) \\ &= 4,000 + 8,000 + 15,000 + 4,000 \\ &= 31,000\text{ truckload-miles} \end{aligned}$$

Step 3: Calculate Weighted Coordinates for the Y-Axis ($\sum W_i y_i$)

$$\begin{aligned} \sum_{i=1}^4 W_i y_i &= (200 \times 30) + (100 \times 20) + (300 \times 90) + (400 \times 60) \\ &= 6,000 + 2,000 + 27,000 + 24,000 \\ &= 59,000\text{ truckload-miles} \end{aligned}$$

Step 4: Solve for Centroid Coordinates ($\bar{x}, \bar{y}$) $\bar{x} = \frac{\sum W_i x_i}{\sum W_i} = \frac{31,000}{1,000} = 31.0\text{ miles}$ $\bar{y} = \frac{\sum W_i y_i}{\sum W_i} = \frac{59,000}{1,000} = 59.0\text{ miles}$

The optimal Center of Gravity location for the central parts warehouse is $(31.0, 59.0)$.

---

Spatial Distance Metrics in Industrial Logistics

The choice of distance metric $d(X, P_i)$ determines how physical transit across geographic terrain is modeled. Different transport modes and civil infrastructures follow different physical geometries:

       RECTILINEAR (L1)                    EUCLIDEAN (L2)                    CHEBYSHEV (L_inf)
      Orthogonal Grid Path             Direct Line-of-Sight                Simultaneous Dual-Axis
         +---------*                          /                                +---------*
         |         |                         /                                 |       / |
         |         |                        /                                  |      /  |
         *---------+                       *                                   *-----+---+
   d = |x1-x2| + |y1-y2|              d = sqrt((x1-x2)^2 + (y1-y2)^2)           d = max(|x1-x2|, |y1-y2|)

1. Rectilinear Distance ($L_1$ Norm / Manhattan Metric)

$d_1(X, P_i) = |x - a_i| + |y - b_i|$

• Physical Basis: Travel is constrained to a grid of orthogonal pathways parallel to coordinate axes (north-south and east-west).
• Applications: Urban street networks (e.g., Manhattan grid), rectangular warehouse aisles served by forklifts and automated guided vehicles (AGVs), and piping/wiring runs on factory utility trusses.
• Iso-Distance Contours: Diamond shapes oriented at $45^\circ$ angles to the coordinate axes.

2. Euclidean Distance ($L_2$ Norm / Straight-Line Metric)

$d_2(X, P_i) = \sqrt{(x - a_i)^2 + (y - b_i)^2}$

• Physical Basis: Direct, unobstructed line-of-sight transit between two points.
• Applications: Airfreight logistics, offshore pipeline networks, marine shipping channels, overhead traveling bridge cranes in heavy fabrication bays, and rural highway systems without rigid grid alignments.
• Iso-Distance Contours: Concentric circles centered at $(a_i, b_i)$.

3. Chebyshev Distance ($L_\infty$ Norm / Chessboard Metric)

$d_\infty(X, P_i) = \max\left( |x - a_i|, |y - b_i| \right)$

• Physical Basis: Travel occurs along both orthogonal axes simultaneously and independently at equal speeds. The time required to reach the destination is governed strictly by the axis requiring the maximum travel distance.
• Applications: Automated Storage and Retrieval Systems (AS/RS). In an AS/RS high-bay aisle, the storage/retrieval machine (SRM) travels horizontally down the track while its hoist carriage elevates vertically up the mast at the same time. If horizontal travel speed $v_x$ and vertical hoist speed $v_y$ are normalized, the transit cycle time is dictated by $t = \max(d_x/v_x, d_y/v_y)$.
• Iso-Distance Contours: Axis-aligned squares centered at $(a_i, b_i)$.

4. General Minkowski $L_p$ Metric

$d_p(X, P_i) = \left( |x - a_i|^p + |y - b_i|^p \right)^{1/p}$ For actual highway networks, actual driving distances are greater than Euclidean ($L_2$) but rarely as rigid as pure rectilinear ($L_1$). Transportation researchers fit empirical road travel using $p \approx 1.5$ to $1.8$, or apply a road circuity factor $k_c$ such that $d_{\text{road}} = k_c \times d_2$, where $k_c \approx 1.15$ to $1.30$ in developed road networks.

Mathematical Comparison Table

Minisum Location Models (The Weber Problem)

The Minisum Location Problem (historically termed the Weber Problem) seeks the coordinates $(x, y)$ of a single facility that minimizes the total weighted transportation cost, where cost is directly proportional to actual linear distance ($L_1$ or $L_2$):

$\min_{x, y} f(x, y) = \sum_{i=1}^m w_i d(X, P_i)$

Where $w_i$ represents the cost-weight product ($w_i = c_i F_i$) moving between the new facility and demand point $P_i = (a_i, b_i)$.

The Rectilinear Minisum Problem: The Median Method

When travel occurs along an orthogonal grid ($L_1$ metric), the Weber objective function exhibits a remarkable mathematical property: it completely decouples into two independent, one-dimensional optimization problems:

$f(x, y) = \sum_{i=1}^m w_i \left( |x - a_i| + |y - b_i| \right) = \sum_{i=1}^m w_i |x - a_i| + \sum_{i=1}^m w_i |y - b_i| = f_1(x) + f_2(y)$

Because $f_1(x)$ depends only on $x$ and $f_2(y)$ depends only on $y$, the optimal coordinate $x^*$ can be determined independently of $y^*$.

The 1-Dimensional Median Theorem

The one-dimensional objective function $f_1(x) = \sum_{i=1}^m w_i |x - a_i|$ is a piecewise linear, convex function. Its derivative is discontinuous at the points $a_i$. The minimum of $f_1(x)$ occurs at the weighted median of the coordinates $a_i$.

Algorithm for the Rectilinear Median Method:

1. Sort along the X-axis: Order all demand points in ascending order of their $X$-coordinates: $a_{(1)} \le a_{(2)} \le \dots \le a_{(m)}$, keeping each point's associated weight $w_{(i)}$.
2. Compute cumulative weights: Calculate the total weight $W = \sum_{i=1}^m w_i$. Determine the half-weight threshold $W / 2$.
3. Identify the X-median: Accumulate the weights $\sum_{j=1}^k w_{(j)}$ starting from the lowest $X$-coordinate. The optimal coordinate $x^*$ is the coordinate $a_{(k)}$ where the cumulative weight first equals or exceeds $W / 2$: $\sum_{j=1}^{k-1} w_{(j)} < \frac{W}{2} \quad \text{and} \quad \sum_{j=1}^k w_{(j)} \ge \frac{W}{2}$ (If the cumulative weight at $a_{(k)}$ exactly equals $W / 2$, any coordinate in the interval $[a_{(k)}, a_{(k+1)}]$ is an alternative optimum.)
4. Repeat for the Y-axis: Independently sort all demand points in ascending order of their $Y$-coordinates: $b_{(1)} \le b_{(2)} \le \dots \le b_{(m)}$. The optimal coordinate $y^*$ is the coordinate $b_{(r)}$ where cumulative weight first equals or exceeds $W / 2$.

Worked Example: Rectilinear Median vs. Centroid

Using the four assembly plants from the previous example:

• Plant 1: $(20, 30)$, $w_1 = 200$
• Plant 2: $(80, 20)$, $w_2 = 100$
• Plant 3: $(50, 90)$, $w_3 = 300$
• Plant 4: $(10, 60)$, $w_4 = 400$

Total weight $W = 1,000$. Half-weight threshold $W / 2 = 500$.

Finding Optimal $x^*$: Sort by $X$-coordinate:

1. Plant 4: $x = 10$, $w = 400$, Cumulative Weight $= 400$ ($< 500$)
2. Plant 1: $x = 20$, $w = 200$, Cumulative Weight $= 600$ ($\ge 500$)

The threshold of 500 is crossed at $x = 20$. Therefore, $x^* = 20$.

Finding Optimal $y^*$: Sort by $Y$-coordinate:

1. Plant 2: $y = 20$, $w = 100$, Cumulative Weight $= 100$ ($< 500$)
2. Plant 1: $y = 30$, $w = 200$, Cumulative Weight $= 300$ ($< 500$)
3. Plant 4: $y = 60$, $w = 400$, Cumulative Weight $= 700$ ($\ge 500$)

The threshold of 500 is crossed at $y = 60$. Therefore, $y^* = 60$.

The optimal rectilinear Minisum location is $(20, 60)$.

Notice the crucial difference: the Center of Gravity yielded $(31.0, 59.0)$, whereas the rectilinear median yielded $(20, 60)$. The centroid was pulled toward Plant 2 ($x=80$) by the squared distance penalty, whereas the median method correctly identified that placing the facility directly at $x=20$ minimizes total linear travel miles.

The Euclidean Minisum Problem: The Weiszfeld Algorithm

When transportation occurs over straight-line distances ($L_2$ Euclidean metric), the objective function is:

$f(x, y) = \sum_{i=1}^m w_i \sqrt{(x - a_i)^2 + (y - b_i)^2}$

Unlike the rectilinear case, the Euclidean minisum problem cannot be decoupled into independent $x$ and $y$ problems, and it possesses no closed-form algebraic solution. The gradient components are coupled through the distance denominator $d_i(x, y) = \sqrt{(x - a_i)^2 + (y - b_i)^2}$:

$\frac{\partial f}{\partial x} = \sum_{i=1}^m \frac{w_i (x - a_i)}{\sqrt{(x - a_i)^2 + (y - b_i)^2}} = 0, \qquad \frac{\partial f}{\partial y} = \sum_{i=1}^m \frac{w_i (y - b_i)}{\sqrt{(x - a_i)^2 + (y - b_i)^2}} = 0$

Rearranging these equations yields the foundation of the Weiszfeld Iterative Algorithm:

$x = \frac{\sum_{i=1}^m \frac{w_i a_i}{d_i(x, y)}}{\sum_{i=1}^m \frac{w_i}{d_i(x, y)}}, \qquad y = \frac{\sum_{i=1}^m \frac{w_i b_i}{d_i(x, y)}}{\sum_{i=1}^m \frac{w_i}{d_i(x, y)}}$

The Weiszfeld Iteration Equations

Starting from an initial iteration $k = 0$ (commonly the Center of Gravity $(\bar{x}, \bar{y})$), the updated coordinates for iteration $k + 1$ are generated by:

$x^{(k+1)} = \frac{\sum_{i=1}^m \frac{w_i a_i}{d_i(x^{(k)}, y^{(k)})}}{\sum_{i=1}^m \frac{w_i}{d_i(x^{(k)}, y^{(k)})}}, \qquad y^{(k+1)} = \frac{\sum_{i=1}^m \frac{w_i b_i}{d_i(x^{(k)}, y^{(k)})}}{\sum_{i=1}^m \frac{w_i}{d_i(x^{(k)}, y^{(k)})}}$

Where: $d_i(x^{(k)}, y^{(k)}) = \sqrt{(x^{(k)} - a_i)^2 + (y^{(k)} - b_i)^2}$

The algorithm iterates until the Euclidean distance between successive estimates falls within an allowable convergence tolerance $\epsilon$ (e.g., $\sqrt{(x^{(k+1)} - x^{(k)})^2 + (y^{(k+1)} - y^{(k)})^2} < 0.001$).

Singularity Condition: If any trial point $(x^{(k)}, y^{(k)})$ lands exactly on an existing demand point $(a_i, b_i)$, the denominator $d_i$ becomes zero. In computer implementations, a small perturbation $\delta = 10^{-6}$ is added to avoid division by zero.

---

Minimax Location Models: Emergency & Critical Service Siting

While the Minisum objective focuses on economic efficiency (minimizing total or average transportation cost), the Minimax Location Problem focuses on equity and worst-case responsiveness. It minimizes the maximum distance or response time from the new facility to any demand point:

$\min_{x, y} g(x, y) = \max_{1 \le i \le m} \left\{ w_i d(X, P_i) \right\}$

Where $w_i$ represents a service priority weight (or the reciprocal of travel speed along the route to point $i$).

+-----------------------------------------------------------------------------------+
|                         MINISUM vs. MINIMAX COMPARISON                            |
+-----------------------------------------------------------------------------------+
| Attribute           | Minisum Location                 | Minimax Location         |
|---------------------|----------------------------------|--------------------------|
| Primary Objective   | Economic Efficiency              | Worst-case Equity / Life |
| Cost Metric         | Total Ton-Miles / Freight Cost   | Maximum Response Time    |
| Typical Systems     | Warehouses, Factories, DCs       | Fire, EMS, Hazmat, Police|
| Penalty Profile     | Linear with cumulative volume    | Governed by furthest node|
| Mathematical Basis  | Decoupled Medians / Weiszfeld    | Enclosing Bounding Boxes |
+-----------------------------------------------------------------------------------+

Rectilinear Minimax: The Coordinate Transformation Method

For unweighted rectilinear distance ($w_i = 1$), the minimax objective seeks:

$\min_{x, y} \max_{1 \le i \le m} \left( |x - a_i| + |y - b_i| \right)$

This non-differentiable problem is solved geometrically by rotating the coordinate axes by $45^\circ$ using the coordinate transformations $u = x + y$ and $v = x - y$. In the rotated coordinate plane, rectilinear distance transforms into the Chebyshev distance ($L_\infty$):

$|x - a_i| + |y - b_i| = \max\left( |u - u_i|, |v - v_i| \right)$

Where $u_i = a_i + b_i$ and $v_i = a_i - b_i$. The optimal solution is determined by finding the midpoints of the enclosing bounding intervals:

$c_1 = \min_i (a_i + b_i), \qquad c_2 = \max_i (a_i + b_i)$ $c_3 = \min_i (-a_i + b_i), \qquad c_4 = \max_i (-a_i + b_i)$ $c_5 = \max(c_2 - c_1, c_4 - c_3)$

The minimax distance is $r^* = c_5 / 2$. Every point on the line segment joining the two points below is optimal:

$(x_1^*, y_1^*) = \frac{1}{2}\left(c_1 - c_3,\; c_1 + c_3 + c_5\right), \qquad (x_2^*, y_2^*) = \frac{1}{2}\left(c_2 - c_4,\; c_2 + c_4 - c_5\right)$

Check: for two demand points at $(0, 0)$ and $(10, 0)$, $c_1 = 0$, $c_2 = 10$, $c_3 = -10$, $c_4 = 0$, and $c_5 = 10$. Both formulas give $(5, 0)$, which is 5 miles from each point, so $r^* = 5$.

This ensures that the furthest emergency caller experiences the minimum possible response delay.

Qualitative Factor Rating Method

While quantitative transportation models optimize logistics costs, site selection in practice involves numerous non-quantifiable, regulatory, and environmental factors that cannot be translated directly into ton-mile freight rates. The Qualitative Factor Rating Method (multi-criteria scoring model) provides a structured, defensible engineering decision framework for synthesizing qualitative and quantitative attributes.

Step-by-Step Engineering Procedure

1. Develop a Comprehensive Factor List: Identify all critical decision criteria impacting facility success, including:
   • Labor availability, technical skill level, prevailing wage rates, and unionization history.
   • Transportation infrastructure (proximity to interstate highways, Class I rail spurs, intermodal ramps, deepwater ports, cargo airports).
   • State and municipal tax climates, corporate tax abatements, capital investment tax credits, and enterprise zone incentives.
   • Utility infrastructure (electrical grid reliability, dual-feed substation availability, water and sewer capacity, natural gas lines, industrial fiber internet).
   • Environmental regulations, air permit attainment status, seismic zone classification, and 100-year flood plain boundaries.
   • Quality of life, cost of living index, public school ratings, and executive housing availability.
2. Assign Normalized Weights ($w_k$): Assign an importance weight to each criterion reflecting corporate strategic objectives. The weights must sum to unity: $\sum_{k=1}^K w_k = 1.00 \quad (\text{or } 100\%)$
3. Establish a Common Scoring Scale: Rate each candidate site on each criterion using a standardized scale (e.g., 1 to 10 or 1 to 100, where 100 represents exceptional capability and 1 represents unacceptable deficiency).
4. Calculate Composite Weighted Scores: For each candidate site $j$, compute the weighted score $S_j$: $S_j = \sum_{k=1}^K w_k \cdot s_{kj}$
5. Perform Sensitivity Analysis: Evaluate how rankings shift if key subjective weights (e.g., labor availability vs. tax incentives) vary by $\pm 20\%$.

Worked Engineering Comparison Table

An industrial engineering site selection team is evaluating three candidate locations (Site Alpha, Site Beta, and Site Gamma) for a new $85M aerospace composites manufacturing facility:

Analysis: Site Beta achieves the highest composite score (80.00), driven heavily by its exceptional rating in the highest-weighted category (Skilled Labor Pool at 95). However, Site Alpha (79.50) and Site Gamma (79.75) are extremely competitive. If the company negotiates enhanced tax abatements at Site Beta or if power reliability at Site Alpha is upgraded via redundant utility feeds, the ranking could invert. This demonstrates why qualitative factor scoring must be coupled with sensitivity analysis before final board approval.