ISE 754: Logistics Engineering, Fall 2026
You can know where before you know what it costs.
Lecture 2.1 sorted location problems by their objective and solved the linear minisum objective by hand, in one dimension and in two dimensions under rectilinear distance. It left two questions open. First, where do the weights w_i in the minisum objective come from? Second, how is the location found when the problem is not one dimensional and the distances are not rectilinear, so that the by-hand median procedure no longer applies? Secs. 1 and 2 describe how the weights are determined and Sec. 3 describes a general-purpose numerical procedure that can be used to determine locations.
Fig. 1 shows a typical supply chain for a single plant. In the figure, items D and E provided by two different tier-two suppliers are shipped to a tier-one supplier where they are used to produce item B, units of which are then shipped to the plant. Similarly, items F and G are used to produce item C, which, along with B, are the two raw materials used to produce the finished good A. Units of A are shipped to distribution centers (DCs), from which they are delivered to the final customer. The number of units shown along each lane (or arc) in the figure is meant to indicate the typical size of each shipment; e.g., one of the reasons that a DC might be used is that it is cheaper to transport larger loads from the plant to the DC and then transport small loads to each customer as needed.
As shown in Fig. 2, a production system can be considered as a node (or facility) in a logistics network that converts raw materials procured from suppliers into finished goods that are distributed to customers. For most production systems, the material input to the system equals the material output from the system. Raw materials are those inputs that are transported to the production facility; ubiquitous inputs (e.g., water) are those available at any location, so that they do not need to be transported. Finished products are those outputs transported from the facility; while scrap is the output that is disposed of locally (although some outputs termed “scrap” are sometimes transported long distances from the facility for disposal or rework).
The bottom portion of Fig. 2 illustrates two different shipping terms that describe when the transfer of title occurs when goods are transported from the seller to the buyer: FOB Origin and FOB Destination, where FOB stands for free on board.1 In most cases, the cost of transporting the goods is paid for by whoever is the owner of the goods during the transport. Taking the production system as the buyer, the supplier (the seller) would pay for the transport of goods from the supplier’s location (the origin) to the production facility (the destination) if the shipping terms were FOB Destination.
The choice of FOB terms does not usually affect the location decision at all; it only fixes who pays for the freight. If a purchaser pays for the transport (FOB Origin), then they would typically pay less than they would if the supplier paid for the transport (FOB Destination). As a result, the sum of the transport cost and the purchase cost is unaffected by the FOB terms, and therefore does not impact the location decision. What typically determines who pays for the transport is who can provide the transport at lower cost. A large retailer with an extensive logistics network and its own fleet of trucks will typically be able to provide the transport at a lesser cost compared to a small supplier.
Assuming that local input costs are either the same at every location or are insignificant as compared to transport costs, the minisum transport-oriented single-facility location problem is to locate a new facility (NF) to minimize the sum of weighted distances between NF and m existing facilities EF_i, i = 1, \ldots, m:
\min_{X}\ TC(X) = \sum_{i=1}^{m} w_i\, d(X, P_i) = \sum_{i=1}^{m} q_i\, r_i\, d(X, P_i) \tag{1}
where
The quantity that drives the location is the monetary weight, and it is worth splitting into two parts. The physical weight is the easy one, tons per year. Multiplying it by the transport rate, in dollars per ton-mile, cancels the tons and leaves a rate in dollars per mile per year: the monetary weight. In practice people drop the qualifier and simply say “the weight,” always meaning this monetary one.
If the physical weight q_i is used, then TC in Eq. 1 is in units of $; if, instead, the physical weight rate f_i in units of tons per year is used, then TC is in units of $/year.
An activity is weight-gaining when its outbound monetary weight exceeds its inbound, and weight-losing in the reverse. The distribution center is the standard weight-gaining case: cases of a product arrive on a tractor-trailer and leave on a smaller delivery truck whose cost per ton-mile might be twice as high, so the transport rate, and with it the monetary weight, roughly doubles even though the product itself has not changed.
When the same vehicles carry the inbound and outbound flows, one transport rate is common to every term. Taking it as one turns the objective into minimizing total distance, so the minimum-cost problem becomes a minimum-distance problem, solved by the same median procedure as in Model 8 in Lecture 2.1. The transport rate has to be included explicitly only when the flows travel at different costs.
A plant draws raw materials from its suppliers and ships finished goods to its customers, so its location is pulled from both sides at once: procurement, the inbound flows, and distribution, the outbound flows. Each existing facility, supplier or customer, enters the minisum objective of Eq. 1 through its monetary weight, and the location balances all of them together.
The simplest case, one supplier and one customer, already shows how the weights decide the location. Two products illustrate the two extremes (Fig. 3 and Fig. 4), with the same transport rate on both sides, so only the physical weights matter.
Product A is physically weight-losing: two tons move in for every one ton that moves out, so every mile the plant sits closer to Asheville saves a net ton-mile of transport, and the plant belongs at the supplier. Product B is weight-gaining: its ubiquitous inputs, such as water or air, are available everywhere and incur no transport, so only one ton moves in for every three tons that move out, and the plant belongs at the customer. If the two products were made in a single shared plant, neither answer stands alone: the location would depend on how much of each product is made, so each product’s annual demand f_i enters, and the problem becomes exactly the weighted minisum of Eq. 1.
The worked example combines both sides at full scale: a plant on the I-40 corridor (Fig. 5) with two suppliers, three customers, and different inbound and outbound transport rates.
Example 1: Plant location with procurement and distribution costs
A product is to be produced in a plant located along I-40. Two tons of raw materials from a supplier in Asheville and a half ton of a raw material from a supplier in Durham are used to produce each ton of finished product that is shipped to customers in Statesville, Winston-Salem, and Wilmington. The demand of these customers is 10, 20, and 30 tons per year, respectively. Raw materials cost $0.33 per ton-mile to ship to the plant; finished goods cost $1.00 per ton-mile to ship from the plant. All scrap is disposed of locally. Determine where the plant should be located so that procurement and distribution costs (i.e., the transportation costs to and from the plant) are minimized, and whether the plant is weight-gaining or weight-losing.
Example 1(a): From the bill of material to the monetary weights
Determine each existing facility’s monetary weight, outbound from the demands and inbound through the bill of material.
The outbound side comes first, because the customers’ demands are given directly: each customer’s monetary weight is its demand times the outbound rate. The inbound side follows from the bill of material (Fig. 6): the customers together take \sum_i f_i = 60 tons of finished product per year, so the plant must pull in two tons of Asheville material and half a ton of Durham material for each of those tons, and each supplier’s monetary weight is that physical flow times the inbound rate.
fout = [10, 20, 30] # ton/yr to the three customers
rout = 1.00 # $/ton-mi outbound
wout = fout .* rout # customer monetary weights
BOM = [2.0, 0.5] # ton raw per ton out: Ashvl, Durham
fin = BOM .* sum(fout) # inbound physical flows (ton/yr)
win = fin ./ 3 # supplier wts; r_in = $1/3 ($0.33)
@show wout fin win;wout = [10.0, 20.0, 30.0]
fin = [120.0, 30.0]
win = [40.0, 10.0]The five monetary weights, in dollars per mile of plant movement, are 40 (Asheville), 10 (Statesville), 20 (Winston-Salem), 10 (Durham), and 30 (Wilmington). The suppliers’ weights land on the same footing as the customers’, so from here the problem no longer cares which facilities ship in and which ship out.
Example 1(b): Locating the plant
Order the facilities along I-40 and apply the median procedure of Model 8 in Lecture 2.1.
With every weight known, the corridor is one-dimensional and the median procedure applies unchanged: sweep west to east, accumulate weight, and stop at the first facility where the running total reaches half the grand total.
city = ["Asheville", "Statesville", "Winston-Salem",
"Durham", "Wilmington"] # west-to-east order
w = [win[1], wout[1], wout[2], win[2], wout[3]] # $/mi, W to E
W = sum(w)
@show W W/2 cumsum(w)
city[findfirst(>=(W / 2), cumsum(w))]W = 110.0
W / 2 = 55.0
cumsum(w) = [40.0, 50.0, 70.0, 80.0, 110.0]"Winston-Salem"The running total reaches 40, then 50, then 70 > 55 at Winston-Salem, and sweeping east to west confirms the same stop (30, 40, then 60 > 55). Fig. 7 shows both sweeps on the corridor.
The mile markers never enter the calculation. As in the one-dimensional median, only the west-to-east ordering of the facilities matters, not the distances between them.
The plant locates at Winston-Salem.
Example 1(c): Weight-gaining or weight-losing
Classify the plant by comparing its inbound and outbound weights, physical and monetary.
@show sum(fin) sum(fout) sum(win) sum(wout);sum(fin) = 150.0
sum(fout) = 60
sum(win) = 50.0
sum(wout) = 60.0The plant is (monetary) weight-gaining since \Sigma w_\text{in} = 50 < \Sigma w_\text{out} = 60, but the plant is physically weight-losing since \Sigma f_\text{in} = 150 > \Sigma f_\text{out} = 60; the difference is due to the higher outbound transport rate, finished goods costing three times as much per ton-mile to move as raw materials. Physically, 150 tons per year arrive to make 60 tons of product, with the difference discarded as scrap. The monetary weights are the ones the location responds to, which is why the plant lands at a customer, pulled toward its market even while shedding mass.
A plant can lose physical weight and gain monetary weight at the same time. The transport rates, not the tons alone, decide which way the location moves.
Setting the supply-chain framing aside, the minisum location problem is at heart a pure optimization problem. Take three existing facilities, as before. For a plain distance objective the goal is to minimize the total distance from the new facility to them; for a location problem each distance is first multiplied by its weight and the total weighted distance is minimized. The optimization logic is the same either way, and from an arbitrary starting point the total is generally not yet at its minimum.
Because that logic does not depend on the distance, one procedure serves every metric the course uses, a straight-line distance in the plane here and the great-circle distance of Sec. 5 on the globe; only the distance function changes. By hand the weighted objective yields only in special cases: the Majority Theorem of Model 8 in Lecture 2.1, where a single facility holding at least half the weight takes the location, and the two-dimensional rectilinear problem, solved by the median along each axis. Otherwise the location is found numerically.
This problem was solved in 1629 by Fermat. For three points, the location that minimizes the total distance to them is a fourth point, the Steiner point, characterized by the three lines from it to the points meeting at equal angles of 120° (Fig. 8), the angle soap bubbles meet at as they settle into a natural minimum-energy state. Seen another way, if roads run from a first city to a second and from the second to a third, adding one intersection point and running three road segments to it covers less total distance than the two direct segments, which is why the major intersections of an interstate network typically sit at these minimizing points.
Before computers the weighted minimum was found mechanically, on what is called the Varignon frame (Fig. 9). Holes are drilled in a board at the existing facilities; a physical weight proportional to each facility’s monetary weight hangs from a string through its hole, the strings meet at a common ring, and the ring settles where it comes to rest, at the point that minimizes the total weighted distance. It is more accurate than the rectilinear approximation but far more work: photographs survive of two or three engineers in the 1930s spending three or four days on a single location this way, whereas the rectilinear-ordering shortcut of Lecture 2.1 is a twenty-minute exercise.
The single-facility minisum is stated as a model carried from words to runnable code (Model 1). Its formulation is a math program in a generic distance; its implementation is a function given the weights, the facilities, and a distance function, returning the location, so the one model serves every metric the course uses.
minimize: the total weighted distance from the new facility to the m existing facilities
solve for:
(a) the location X of the new facility, any point in the plane.
subject to: none
return: the location X^\star of the new facility
assumptions:
(a) distance is a metric, for example the straight-line distance used below or the great-circle distance of Sec. 5;
(b) the cost of each flow is proportional to distance, so each monetary weight w_i is a constant.
\min_{X}\; TC(X) \;=\; \sum_{i=1}^{m} w_i\, d(X, P_i)
where:
# Model: single-facility minisum
function minisum(w, P, dist)
x0 = vec(sum(P; dims = 1)) ./ size(P, 1) # centroid start
TC(X) = sum(w[i] * dist(X, P[i, :]) for i in axes(P, 1))
return optimize(TC, x0).minimizer
endminisum (generic function with 1 method)The implementation names no method. It builds the objective and hands it to Optim.optimize, the general-purpose optimizer of the Julia ecosystem, which searches numerically; how that search works, and why the result can be trusted, is taken up after the example.
Example 2
Determine the minisum location for three facilities of equal weight at (1, 1), (6, 1), and (6, 5), under straight-line distance.
Binding the data and a straight-line distance to Model 1 and calling its implementation gives the location directly.
# d2 from 2-Loc-2.ipynb cell 3 (straight-line = Logjam.d2)
d2(x1, x2) = length(x1) == length(x2) ?
sqrt(sum((x1 .- x2).^2)) :
error("Inputs not same length.")
P = [1.0 1.0; 6.0 1.0; 6.0 5.0] # three EFs, one per row
w = [1.0, 1.0, 1.0] # equal weights
Xᵒ = minisum(w, P, d2) # the minisum location2-element Vector{Float64}:
5.089566905043874
1.9662907621468446X^\star \approx (5.09,\ 1.97), the Steiner point of Fig. 8.
The search reached that point without any care about where it started, and the reason is the shape of the objective. When the objective is convex, a single bowl (Fig. 10), a downhill search reaches the bottom from wherever it starts, so any starting point gives the same optimum. The weighted-minisum objective is convex, so choosing a start is not a concern here. A nonconvex objective is different: plotted, it shows a series of separate minima, and the starting point then has a large effect on which one the search settles in.
# From 2-Loc-2.ipynb cells 20–26 (tuples → matrix P)
using CairoMakie
xrng = 0:0.1:7
yrng = 0:0.1:6
Z = [sum(d2([x, y], P[i, :]) for i in axes(P, 1))
for x in xrng, y in yrng]
fig = Figure(size = (760, 320))
a1, _ = contour(fig[1, 1], xrng, yrng, Z; levels = 18)
scatter!(a1, P[:, 1], P[:, 2]; color = :firebrick)
surface(fig[1, 2], xrng, yrng, Z;
axis = (type = Axis3, azimuth = pi / 4))
fig
The optimizer steps downhill by one of two kinds of method. A gradient method, such as a quasi-Newton method, uses the objective’s derivatives; but since even simple location objectives rarely have an analytical gradient, it estimates one by finite differencing, which works well only on smooth objectives. A direct method, chiefly Nelder-Mead, never computes derivatives at all, reading only values of the objective, so it copes with the kinked or discontinuous objectives a gradient method cannot, less efficiently but far more robustly. Nelder-Mead is the basic method behind a spreadsheet’s Solver and the default in most optimization packages, and it is the default used here (Fig. 11).
# From 2-Loc-2.ipynb cell 12 (fargplot NM path; tuples → matrix P)
TC(X) = sum(d2(X, P[i, :]) for i in axes(P, 1))
fargplot(f, x, kw) = (scatter!(x...; kw...); f(x))
scatter(P[:, 1], P[:, 2]; color = :firebrick, markersize = 11)
kw = (; color = (:steelblue, 0.6), markersize = 5)
Xᵒ2 = optimize(x -> fargplot(TC, x, kw), [0.0, 0.0]).minimizer
scatter!([Xᵒ2[1]], [Xᵒ2[2]]; color = :red, marker = :star5,
markersize = 14)
current_figure()
Nelder-Mead is also called the amoeba method, for the way it moves. It works a triangle of points across the plane, evaluating the objective at each: from the current triangle it tries a new point, and when that point improves on the worst corner it moves the triangle that way, creeping downhill by function evaluations alone. The step rules that decide each move are collected in the callout below; the working point is that the method needs nothing but the objective’s values.
Each iteration orders the simplex’s corners from best to worst and tries to replace the worst. It first reflects that corner through the opposite edge: the reflected point is accepted if it is merely middling, stretched further by an expansion if it is the new best, and, if it is worse than the second-worst corner, pulled back by a contraction, taken outside the edge when the reflection still beats the current worst corner and inside when it does not. When even the contraction fails to improve, the whole simplex shrinks toward the best corner and the iteration restarts. The search reads only values of the objective, never its derivatives, which is what lets it handle kinked or discontinuous objectives a gradient method cannot.2 Fig. 12 shows why the four moves are fewer decisions than they look: reflection, expansion and both contractions all step along the same ray through the worst corner and the centroid, differing only in how far, and shrink is the only one that leaves it.
algorithm nelder-mead;
{ input: objective f; starting point x₀ ∈ ℝⁿ }
{ output: a vertex minimizing f }
begin
form a simplex of n+1 vertices about x₀;
repeat
order so that f(v₁) ≤ f(v₂) ≤ ⋯ ≤ f(vₙ₊₁);
c := centroid of the n best, v₁ … vₙ; { drop vₙ₊₁ }
xᵣ := c + (c − vₙ₊₁); { reflection }
if f(v₁) ≤ f(xᵣ) < f(vₙ) then vₙ₊₁ := xᵣ
else if f(xᵣ) < f(v₁) then
begin
xₑ := c + 2(c − vₙ₊₁); { expansion }
if f(xₑ) < f(xᵣ) then vₙ₊₁ := xₑ else vₙ₊₁ := xᵣ;
end
else if f(xᵣ) < f(vₙ₊₁) then
begin
x_c := c + ½(xᵣ − c); { outside contraction }
if f(x_c) ≤ f(xᵣ) then vₙ₊₁ := x_c else shrink;
end
else
begin
x_c := c + ½(vₙ₊₁ − c); { inside contraction }
if f(x_c) < f(vₙ₊₁) then vₙ₊₁ := x_c else shrink;
end
until the simplex is small enough;
return v₁;
end;
{ shrink: replace each vᵢ, i = 2 … n+1, by v₁ + ½(vᵢ − v₁) }Optim and the plotting used here are general-purpose Julia. The logistics-specific layer that the rest of the topic relies on, place names and distances over real geography, is introduced next.
The single-facility minisum of Sec. 3 was solved with optimize from Optim, a general-purpose optimizer taken from the Julia ecosystem. Logistics problems, though, are rarely posed in the language a solver expects: their data are place names, distances over real roads, and freight rates, not coordinate vectors and objective functions. Logjam is the Julia package written for this course to supply that missing layer. It does not reimplement what the ecosystem already does well, such as the optimization in Optim or the plotting in CairoMakie; it adds the logistics domain on top, turning place names into coordinates and measuring distance over real geography.
The distance and mapping functions the course needs were not available in any existing package, so they were gathered into a package of their own, Logjam. Packages in the general Julia registry are added by name; Logjam is still an early version and is not registered, so it is added instead by giving the URL of its GitHub repository, with registration by name left as a possibility for later.
Logjam collects tools across the logistics topics of the course: facility location, including discrete p-median and uncapacitated facility-location heuristics and continuous location by alternating location-allocation (ala) and the weighted center of gravity (wcentroid); transportation costing; network analysis and vehicle routing over the FAF5 highway network; a built-in U.S. gazetteer of cities, counties, and ZIP codes; forward and reverse geocoding; and the distance metrics used in this lecture.
The package cell in Lecture 1.2 adds Logjam by URL:
using Pkg
Pkg.add(url="https://github.com/mgkay/Logjam")Its features activate through package extensions, so its capabilities grow as more packages are loaded: using Logjam alone provides the core location, costing, distance, and geocoding tools; adding CairoMakie enables plotting; adding GeoMakie as well enables the mapping functions, such as makemap and aligntext, used in Sec. 5. The facility-location and routing maps introduced in Lecture 1.1 were drawn with these same mapping tools, and the package’s README is the reference for its full function set. The remainder of this lecture puts Logjam to work computing the distances that the minisum objective needs.
The minisum objective needs a distance d(X, P_i) between the facility and each existing point. In the plane, that distance belongs to a family of metrics selected by a single parameter; over the curved surface of the earth, it is the great-circle distance. Logjam supplies both, together with the geocoding that turns place names into the coordinates they require. The planar metrics come first.
Given two points P_1 = (x_1, y_1) and P_2 = (x_2, y_2), a general metric distance function termed the l_p distance or norm can be used to represent several of the most common distances used in practice. The function uses different values of the parameter p to represent these special cases:
d_p(P_1, P_2) = \left[\,|x_1 - x_2|^p + |y_1 - y_2|^p\,\right]^{1/p}, \quad p \ge 1 \tag{2}
corresponding to the rectilinear (p = 1), Euclidean (p = 2), and Chebychev (p \to \infty) distances:
d_1(P_1, P_2) = |x_1 - x_2| + |y_1 - y_2| \tag{3}
d_2(P_1, P_2) = \sqrt{(x_1 - x_2)^2 + (y_1 - y_2)^2} \tag{4}
d_\infty(P_1, P_2) = \max\{\,|x_1 - x_2|,\; |y_1 - y_2|\,\} \tag{5}
The above distances are defined for two-dimensional points, but they can be extended for points of any dimension. For 0 < p < 1, the l_p distance is not defined since the triangle inequality does not hold; for 1 < p < 2, the l_p distance lies between the rectilinear and the Euclidean distances; and for 2 < p < \infty, the l_p distance lies below the Euclidean distance and decreases as p increases.
Each value of p has a characteristic contour of points at equal distance from the center (Fig. 13): a diamond for the rectilinear distance, the familiar circle for the Euclidean, and a square for the Chebychev, with the intermediate values filling the range between them.
Logjam computes each of these directly. d1 and d2 take two points; the unified dists takes point sets and a metric parameter, so a single call covers the whole family, down to the Chebychev limit at Inf.
P1 = [0, 0] # two points in the plane
P2 = [4, 3]
d1(P1, P2), d2(P1, P2) # rectilinear, Euclidean(7.0, 5.0)X1 = [0 0] # point sets, one row per point
X2 = [4 3]
[dists(X1, X2, p)[1] for p in (1, 2, Inf)] # l₁, l₂, l∞3-element Vector{Float64}:
7.0
5.0
4.0Without loss of generality, let P_1 = (x, y), for x, y \ge 0, and P_2 = (0, 0). Then d_\infty(P_1, P_2) = \max\{x, y\} and d_p(P_1, P_2) = (x^p + y^p)^{1/p}.
If x = y, then \lim_{p \to \infty}(x^p + y^p)^{1/p} = \lim_{p \to \infty}(2x^p)^{1/p} = \lim_{p \to \infty} 2^{1/p}\, x = x .
If x < y, then \lim_{p \to \infty}(x^p + y^p)^{1/p} = \lim_{p \to \infty}\!\left[\left(\tfrac{x}{y}\right)^{p} + 1\right]^{1/p} y = 1 \cdot y = y .
A similar argument can be made if x > y. ∎
The Chebychev distance is the natural one when two motions run at once, so the longer of the two sets the total: a lift truck drives down the aisle while raising its forks, and an overhead crane bridges and trolleys at the same time. For a move with |\Delta x| = 6 and |\Delta y| = 3 (Fig. 14), the rectilinear distance adds the two (9) and the Euclidean takes the hypotenuse (\approx 6.71), but the time is set by the longer move alone, \max\{6, 3\} = 6.
The great circle distance is an exercise in spherical trigonometry. A compact formula gives the separation of two points as an angle in radians, as though on a sphere of unit radius; that angle becomes a distance by multiplying it by the radius of the earth.
Great circle, or geodesic, distances on the surface of a sphere, such as the earth, are the shortest distance between two points on the surface along the circle formed by the intersection of the surface and a plane passing through the center of the sphere (a “great circle”). The elevations of the points on the surface are usually ignored.
The great circle distance d_{GC} between points 1 and 2 on the surface of the earth, specified by their longitude (lon) and latitude (lat) angles in radians, is d_{GC} = R\,d_{rad}, where R is the radius of the earth and the central angle d_{rad} (in radians) is
d_{rad} = 2 \sin^{-1} \min\left\{1,\ \sqrt{\sin^2\!\left(\tfrac{y_1 - y_2}{2}\right) + \cos(y_1)\cos(y_2)\,\sin^2\!\left(\tfrac{x_1 - x_2}{2}\right)}\right\} \tag{6}
with x and y the longitude and latitude. To convert between decimal degrees and radians:
x_{rad} = \frac{x_{deg}}{180}\,\pi \qquad \text{and} \qquad x_{deg} = \frac{x_{rad} \cdot 180}{\pi}. \tag{7}
In Fig. 15, a, b, and c are the sides of a spherical triangle and A, B, and C are the corresponding angles. The great circle distance between points 1 and 2 in radians (d_{rad}) corresponds to side c of the triangle.
Given
a = 90^\circ - y_2, \quad b = 90^\circ - y_1, \quad C = x_1 - x_2
the formula for great circle distances can be derived using the spherical law of cosines for sides:
\begin{aligned} \cos(c) &= \cos(a)\cos(b) + \sin(a)\sin(b)\cos(C)\\ &= \cos(90^\circ - y_2)\cos(90^\circ - y_1) + \sin(90^\circ - y_2)\sin(90^\circ - y_1)\cos(x_1 - x_2)\\ &= \sin(y_2)\sin(y_1) + \cos(y_2)\cos(y_1)\cos(x_1 - x_2) \end{aligned}
Solving for c:
d_{rad} = c = \cos^{-1}\!\left[\sin(y_2)\sin(y_1) + \cos(y_2)\cos(y_1)\cos(x_1 - x_2)\right]
The above formula can result in round off error if the two points are located at exactly opposite sides of a sphere. Instead, the Haversine formula can be used (Eq. 6). The Haversine formula does not need to be used if the great circle distance is being calculated by hand and the two points are known to not be on opposite sides of the sphere.
The formula applied to a few points shows how far apart the distances can be. The chapter’s worked example measures three cities from Raleigh, each given as a (longitude, latitude) pair in decimal degrees; dists with the :mi unit (or :rad) returns the great-circle distance to each. Fig. 16 draws the same three distances on the globe. Each arc is the shortest path there is between its two points, and none of them is the straight line a flat map would draw between the same pair.
Raleigh = [-78.659 35.8219] # one point, as a row
XY = [-82.336 29.6742 # Gainesville, FL
44.3667 33.2333 # Baghdad
-43.2 -22.95] # Rio de Janeiro
prt(DataFrame(
City = ["Gainesville, FL", "Baghdad", "Rio de Janeiro"],
Radians = vec(dists(Raleigh, XY, :rad)),
Miles = vec(dists(Raleigh, XY, :mi)))) City Radians Miles
────────────────────────────────────
Gainesville, FL 0.12 475.37
Baghdad 1.62 6,412.00
Rio de Janeiro 1.18 4,675.44
A great-circle distance needs coordinates, and Logjam carries a built-in gazetteer of U.S. places to supply them. usplace() returns every city, town, and census-designated place as a table of names, coordinates, and populations; filter with st2fips, the state’s numeric code, narrows it to the places of interest, here the ten North Carolina cities above one hundred thousand people.
# NC cities over 100k (2-Loc-3 cells 17-18, Logjam gazetteer)
df = filter(r -> r.STFIP == st2fips(:NC) &&
r.POP > 100_000, usplace())
select!(df, :NAME, :LON, :LAT, :POP)
prt(df) NAME LON LAT POP
───────────────────────────────────────
Cary -78.82 35.78 174,721
Charlotte -80.83 35.21 874,579
Concord -80.64 35.39 105,240
Durham -78.90 35.98 283,506
Fayetteville -78.97 35.08 208,501
Greensboro -79.83 36.10 299,035
High Point -79.99 35.99 114,059
Raleigh -78.64 35.83 467,665
Wilmington -77.89 34.21 115,451
Winston-Salem -80.26 36.10 249,545Each row’s longitude and latitude are one point P_i, and dgc returns the great-circle distance between two of them. Raleigh to Charlotte is about one hundred thirty miles, as the crow flies.
ll(nm) = collect(filter(r -> r.NAME == nm,
df)[1, [:LON, :LAT]])
xyR = ll("Raleigh") # (lon, lat) in degrees
dgc(xyR, ll("Charlotte")) # great-circle miles130.3903764072013Broadcasting dgc over the whole set gives Raleigh’s distance to every city at once, and the standard index queries answer the natural questions: the farthest city, and the three nearest.
pt = collect.(zip(df.LON, df.LAT)) # one [lon,lat] per city
d = dgc.([xyR], pt) # Raleigh to each city
df[argmax(d), :NAME] # farthest city"Charlotte"df[sortperm(d)[2:4], :NAME] # three nearest (skip Raleigh)3-element Vector{String}:
"Cary"
"Durham"
"Fayetteville"Plotting longitude and latitude directly on flat axes distorts the map, so the geography is shown through a projection. The most convenient here is the Mercator projection, which has the useful property that a straight line between two points is a constant compass bearing, the line a navigator would steer.
The projection leaves longitude unchanged and stretches the latitude axis toward the poles, placing a point at latitude \varphi at ordinate \sinh^{-1}(\tan\varphi), so equal steps of latitude spread apart with distance from the equator. Logjam’s makemap applies the projection through GeoMakie, so the maps in this section are already Mercator-projected; Fig. 17 shows the effect on a world coastline.
These distances feed straight back into the single-facility minisum. The model of Model 1 never depended on the metric, so the same objective, now with the great-circle distance in place of the straight-line one, locates a facility over real geography. Taking each place’s population as its weight w_i turns the objective into total population-weighted travel, and its minimizer into the location that is, on average, closest to everyone.
Example 3
Locate a single facility to minimize the total population-weighted great-circle distance to every place in North Carolina.
Every place enters, not just the large cities: usplace filtered to the state, each place weighted by its population, with the weighted centroid wcentroid as the starting point for the search.
# every NC place, population-weighted (2-Loc-3 cells 40-41)
df = filter(r -> r.STFIP == st2fips(:NC), usplace())
select!(df, :NAME, :LON, :LAT, :POP)
pt = collect.(zip(df.LON, df.LAT)) # [lon,lat] per place
TC(xy) = sum(df.POP .* dgc.([xy], pt)) # pop-weighted miles
w = wcentroid(df.LON, df.LAT, df.POP) # weighted-centroid start
xyᵒ = optimize(TC, [w.LON, w.LAT]).minimizer2-element Vector{Float64}:
-79.70921230976549
35.64959801955728The nearest place names the answer. Fig. 18 sets that optimum against every place in the state, which is what makes the result readable: the ten largest cities are labeled, and the population-weighted optimum sits near none of them.
d = dgc.([xyᵒ], pt) # optimum to each place
df[argmin(d), :NAME], round(minimum(d)) # nearest place, miles("Franklinville", 6.0)The optimum lies about six miles from Franklinville, near the geographic center of the state.
The proximity of the population-weighted minisum location to the location chosen for the North Carolina Zoo near Asheboro seems most likely to just be coincidental. The minisum location is the one that, on average, minimizes the distance everyone in the state must travel, which would be a natural criterion for a facility the whole public visits. A zoo is such a facility. Accounts of the site’s selection do mention its central position in the state, but none gives evidence that centrality was a principal factor in the decision.3 Given what computing could routinely support in the early 1970s, it is also unlikely that a decision of that period rested on an analysis of the kind carried out here.
The population-weighted minisum location is a continuous answer, free to fall anywhere. A real single-facility decision usually chooses among a handful of actual sites, and one such choice played out in North Carolina. When FedEx was siting its Mid-Atlantic ground hub, the state weighed offering incentives to draw it to the North Carolina Global TransPark at Kinston, in the east. Kinston is a weak ground-hub site: it lies east of I-95, the easternmost interstate, and is served by no interstate of its own. Greensboro, by contrast, sits at the crossing of I-40 and I-85, central to Charlotte, the Triangle, and the Triad (Fig. 19).
An analysis carried out for the state, without access to FedEx’s own figures, used the guesstimation methods of Lecture 1.1 to put the annual cost penalty of the Kinston site over Greensboro at roughly $12 million a year, the size of the subsidy it would take to offset it. Years later a FedEx analyst who had run the company’s own internal study put the figure within a million or two of that estimate, a striking check on guesstimation for a quantity no outsider could measure directly. FedEx built its hub at Piedmont Triad International in Greensboro; the TransPark, long criticized for its distance from major highways, did not win it.4
The population-weighted center and the FedEx decision point the same way, to the central Piedmont, and together they close the loop opened in Lecture 2.1: where a facility sits is a long-lived, costly commitment, and the tools of this lecture, distances over real geography and a minisum solved numerically, are what turn that commitment into a computation. Choosing among discrete candidate sites, rather than a free point in the plane, is the subject the next lectures take up.
The following sources are recommended for further study:
Ahuja, R.K., Magnanti, T.L., and Orlin, J.B., 1993, Network Flows: Theory, Algorithms, and Applications, Englewood Cliffs, NJ: Prentice-Hall.
Ballou, R.H., 2004, Business Logistics/Supply Chain Management, 5th Ed., Englewood Cliffs, NJ: Pearson.
Daskin, M., 1995, Network and Discrete Location: Models, Algorithms, and Applications, New York: Wiley.
Francis, R.L., McGinnis Jr., L.F., and White, J.A., 1992, Facility Layout and Location: An Analytical Approach, 2nd Ed., Englewood Cliffs, NJ: Prentice-Hall.
Larson, R.C., and Odoni, A.R., 1981, Urban Operations Research, Englewood Cliffs, NJ: Prentice-Hall.
Mirchandani, P.B., and Francis, R.L., Eds., 1990, Discrete Location Theory, New York: Wiley.
Simchi-Levi, D., Kaminsky, P., and Simchi-Levi, E., 2008, Designing and Managing the Supply Chain: Concepts, Strategies, and Case Studies, 3rd Ed., Boston: McGraw-Hill.
“What are FOB shipping terms?” simplestudies.com, http://simplestudies.com/what-are-fob-shipping-terms.html (accessed January 2012); carried from the Facility Location chapter.↩︎
J. C. Lagarias, J. A. Reeds, M. H. Wright, and P. E. Wright, “Convergence Properties of the Nelder-Mead Simplex Method in Low Dimensions,” SIAM Journal on Optimization 9, no. 1 (1998): 112–147.↩︎
The zoo’s site was chosen in 1971 from six candidates, Albemarle, Concord, Statesville, Butner, Chatham County, and Asheboro, after a two-year search by a state site-selection committee; the conditions set for it were a parcel of at least 1,000 acres and a pledge of at least $150,000 from the host community. “The 1970s: Room to Roam at the North Carolina Zoo,” Our State, https://www.ourstate.com/the-1970s-room-to-roam-at-the-north-carolina-zoo/ (accessed August 2026). The chosen parcel, near Purgatory Mountain in the Uwharrie range of Randolph County, is described as lying “near the geographic center of the state”: William S. Powell, “North Carolina Zoological Park,” NCpedia, 2006, https://www.ncpedia.org/north-carolina-zoological-park (accessed August 2026). Neither account gives distance to the state’s population as a selection criterion, and the center computed in this example is population-weighted rather than geographic.↩︎
FedEx sited its Mid-Atlantic hub at Piedmont Triad International Airport in Greensboro; the North Carolina Global TransPark at Kinston was the promoted alternative, frequently criticized for its location away from major highways. “Global TransPark Seeks FedEx Hub,” WFMY News 2 (2003), https://www.wfmynews2.com/article/news/local/global-transpark-seeks-fedex-hub/83-404033858.↩︎