2.6: Logistics Network Design

ISE 754: Logistics Engineering, Fall 2026

A firm cannot easily tell you where its trucks went. It can tell you what they all cost, and that is enough.

No new Julia packages used.

No new Logjam functions used.

Companion script

2-loc-6.jl (PopcoData.csv)

1. What makes network design hard

A network is a set of nodes, the connections between them, and something that moves along those connections.1 Little else is common to all of them, which is why the idea reaches so far: a road network’s nodes are intersections and its connections streets, carrying vehicles; a power grid’s are substations and lines, carrying electricity. A node is a place where something is held, transformed or redirected; a connection is a route between two nodes. Both carry characteristics, and those distinguish one network from another. A node has a location, a capacity, and a cost of being there at all; a connection has a length, a capacity, and a cost per unit of what travels along it.

A logistics network gives them their logistics meaning. The nodes are facilities, plants and warehouses and distribution centers, and what matters most is where a facility is, because its location is usually what is being chosen. The connections are flows of materials, and what matters is how much moves and how far.

Network design is choosing that structure rather than taking it as given: how many nodes, where they go, which are connected, and the characteristics of each. Logistics network design is choosing how many facilities a system should have, where, and which demand each serves. The flows follow: once the facilities are placed and the customers assigned, what moves between them is determined. As introduced in Lecture 1.1, these types of decisions are often strategic because they involve, for example, committing capital to construct buildings in place for decades; they are tactical when a design must determine how production and inventory will be controlled throughout the network; and they are almost never operational, deciding, for example, this week’s deliveries.

Show the code that draws this figure
# The network half of lecture 2.2's weighted-flow figure, at the same
# geometry and in the same palette: node color means what kind of
# facility a circle is, and nothing else is colored.
supc = colorant"#1f6fc4"          # suppliers
nfc  = colorant"#1e8a3c"          # the new facility
cusc = colorant"#d21f26"          # customers
ink  = colorant"#252525"
rn   = 0.42                       # node radius, data units

y1, y2 = 3.55, 1.45               # suppliers
y3, y4, y5 = 4.05, 2.50, 0.95     # customers
ynf = 2.50
xsup, xnf, xcus = 11.1, 13.0, 14.9

fig = Figure(size = (520, 300); backgroundcolor = :white)
ax = Axis(fig[1, 1]; backgroundcolor = :white, aspect = DataAspect())
hidedecorations!(ax); hidespines!(ax)
limits!(ax, xsup - 1.0, xcus + 1.0, 0.3, 4.7)

# Each arrow runs edge to edge, so the star reads as connected.
function flow!(p, q)
    d = q .- p
    u = d ./ sqrt(sum(abs2, d))
    a, b = p .+ u .* rn, q .- u .* rn
    arrows2d!(ax, [a[1]], [a[2]], [b[1] - a[1]], [b[2] - a[2]];
              color = ink, shaftwidth = 2, tipwidth = 10, tiplength = 11)
end

for p in ([xsup, y1], [xsup, y2])
    flow!(p, [xnf, ynf])
end
for q in ([xcus, y3], [xcus, y4], [xcus, y5])
    flow!([xnf, ynf], q)
end

function node!(x, y, label, color; fs = 17)
    poly!(ax, Circle(Point2f(x, y), rn);
          color = :white, strokecolor = color, strokewidth = 2.5)
    text!(ax, x, y; text = label, color = color, fontsize = fs,
          align = (:center, :center))
end

node!(xsup, y1, "1", supc); node!(xsup, y2, "2", supc)
node!(xnf, ynf, "NF", nfc; fs = 14)
node!(xcus, y3, "3", cusc); node!(xcus, y4, "4", cusc)
node!(xcus, y5, "5", cusc)
fig
Figure 1: The network of lecture 2.2’s procurement-and-distribution example. Five nodes are fixed: two suppliers in blue, three customers in red. One node is to be determined, the new facility in green. The flow on each arc follows from the customers’ demand and from the bill-of-materials relationships that supply the plant.

Logistics networks have been used in several of the examples considered in previous lectures. Fig. 1 from Ex. 1 in Lecture 2.2 is a good example. Six nodes, five fixed and one, the new facility, determined by the problem. Five arcs. And the flow on each was not given: it followed from the demand at each customer and from the bill-of-materials relationships setting how much raw material the plant draws from each supplier to meet it. The logistics network was designed and used to determine the monetary weights needed to locate the new facility.

Uncertainty is not what makes logistics network design hard. It matters, and a proper logistics engineer accounts for it, but it does not drive the gross, macroscopic shape of a network. What drives that shape is lumpiness: the resources a network is built out of do not come in arbitrary sizes. And the lumpiest resource of all, the one that drives everything else, is people.

What the modeling has to carry is that lumpiness drives economies of scale, so that more of one thing is cheaper per unit, and economies of scope, so that making many things together is cheaper than making them apart, because the same lumpy resources serve all of them. Between them those reduce the cost of a network to two components: a fixed production cost at each facility, and a variable transport cost that depends on where the facilities are. That is the pair the uncapacitated facility location model takes as input, so the UFL determines the number and location of the new facilities.

The fixed production cost never has to be identified directly, and that is what makes the data obtainable. Before, when the UFL heuristics were the subject, a fixed cost and a variable cost were simply put in to illustrate the mechanics. For a real problem, separating the part of a plant’s cost that genuinely does not depend on output, the people and the leases, is difficult and somewhat arbitrary, and it is not where the economies mostly come from: those come from operating a larger machine at a lower cost per unit, and from running larger batches. So the parts are not separated at all. Total production cost is added up, a linear regression is fitted, and the intercept stands in for the fixed component, as in lecture 2.4 Sec. 8.

Lumpiness has a third consequence alongside those two: a facility has a minimum effective size, below which it is not worth having at all.

2. Bottom-up vs. top-down analysis

Almost every analysis of a logistics network, at the level where the number and location of facilities is decided, is one of two kinds. The difference between them is not the model, which is the same model, and not the objective, which is the same objective. It is where the transport rate comes from.

In a bottom-up analysis the rate is given. A carrier quotes a figure, or the firm knows what it pays per loaded mile, and total cost is built upward from that figure: rate times distance times volume, summed over the customers. Every location example so far in this course has been of this kind.

A firm that has been operating for years knows what it spent last year on outbound transportation, because it wrote the checks. It knows where its customers are and how much product went to each of them. What it very often does not know is the dollars per mile behind that spending, because no single number was ever paid: the total is an aggregate over many carriers, many lanes and several kinds of vehicle. Those three facts are enough. Divide what was spent by the demand-weighted distance it bought, and the result is a nominal transport rate, in dollars per ton-mile.

Which of the two applies is determined by what exists, not by preference. Bottom-up is the route when there is no existing system. A network being designed from the ground up has no spending history to reverse, so as long as a transport cost can be estimated, a solution is built upward from it and the differential analysis follows.

Top-down requires an operating system, because what it reverses is that system’s own spending. It is the more strategic of the two, and it extends naturally to a competitor: a rival’s network can be read the same way, and the question of how that rival would respond to a new facility becomes answerable in the same terms. A third case sits between the two. A rate calibrated on one region can be carried to another where the characteristics are expected to be similar, which is a top-down calibration serving a bottom-up design.

Fig. 2 is the difference in one picture. Bottom-up starts from a rate and evaluates each site being compared. Top-down starts from a cost that was actually incurred, turns it into a rate, and only then can evaluate anywhere else.

Figure 2: Bottom-up begins with a known rate and fans out to evaluate each candidate site. Top-down is a chain: an incurred cost becomes a nominal rate, and only then can a new site be evaluated.

Read the arrows and the asymmetry is plain. Bottom-up has one input and fans out to as many sites as are worth evaluating. Top-down is a chain: nothing can be evaluated until the rate exists, and the rate cannot exist until something has already been spent.

Three-city total logistics cost

The two analyses are easiest to tell apart on one instance rather than two, and lecture 2.2 already built the instance: a company whose owners are in Cary, and three customers it ships to. What follows evaluates that firm twice, changing only what is known about it.

Example 1: Cary firm

Determine the increase in annual transport cost from staying in Cary rather than at the transport optimum, first from a given rate and then from what the firm spent last year, and determine what the two answers have in common.

Example 1(a): The rate is given (bottom-up)

Determine the optimum, the cost at Cary, and the difference between them, for three customers receiving 40, 25 and 35 truckloads a year at $2.00 per loaded mile.

# Code block 1: three cities, road against great-circle
Detroit     = [-83.1022, 42.3830]   # (lon, lat), degrees
Gainesville = [-82.3492, 29.6807]
Memphis     = [-89.9666, 35.1090]

gc = [dgc(Detroit, Gainesville),    # great circle, mi
      dgc(Detroit, Memphis),
      dgc(Gainesville, Memphis)]

road = [1024.2, 710.8, 632.8]       # road network, mi

g = road ./ gc

ḡ = sum(g) / length(g)              # the region's factor
1.1314763826552756
# Code block 2: truckload weights and the annual cost
P = permutedims(hcat(Detroit, Gainesville, Memphis))  # lon/lat, blk 1
wTL = [40.0, 25.0, 35.0]           # truckloads per year
rate = 2.00                        # $ per loaded mile

TCyr(xy) = rate * ḡ *
    sum(wTL[i] * dgc(xy, P[i, :]) for i in 1:3)
TCyr (generic function with 1 method)

The optimum is found the way 2.2 found one, a downhill search from the weighted centroid.

# Code block 3: the optimum site
c0 = wcentroid(P[:, 1], P[:, 2], wTL)
xᵒ = Optim.minimizer(optimize(TCyr, [c0.LON, c0.LAT]))
lonlat2loc(xᵒ, usplace()).desc
"3.8 mi E of Charlotte, TN"

Evaluating the two sites is then one call each.

# Code block 4: what an alternative site would cost
Cary = [-78.8190, 35.7814]
Δ = TCyr(Cary) - TCyr(xᵒ)
prt(DataFrame(
    Site = ["optimum", "Cary, NC", "increase"],
    Cost = round.([TCyr(xᵒ), TCyr(Cary), Δ], digits = 0)))
      Site     Cost
───────────────────
   optimum   87,149
  Cary, NC  122,536
  increase   35,386

Staying in Cary costs about $35.4k per year more than the transport optimum, an increase of 40.6%.

Example 1(b): The rate is backed out (top-down)

Determine the same increase for the same firm, given not the rate but what the firm spent last year, and determine what the two answers have in common.

Nothing about the firm changes. What changes is what is known about it: the $2.00 per loaded mile is withdrawn, and in its place is a single figure for what the Cary operation spent on outbound transport last year, together with the tonnage behind it at ten tons per truckload.

The circuity factor does not appear anywhere below, and its absence is not an approximation. A nominal rate is calibrated against the same distances the model will later evaluate, so multiplying every distance by a constant inflates the distances and deflates the rate by exactly the same factor. It cancels. Applying it consistently would not be wrong, merely extra work for an unchanged answer, and great-circle distance is used throughout.

# Code block 5: what the firm knows, and the rate it implies
tonperTL = 10.0                    # ton per truckload
f = wTL .* tonperTL                # ton/yr to each customer
TC₀ = TCyr(Cary)                   # $/yr spent from Cary last year

TD₀ = sum(f[i] * dgc(Cary, P[i, :]) for i in eachindex(f))
rnom = TC₀ / TD₀                   # $/ton-mi
0.2262952765310551

The objective is rebuilt from the nominal rate, and re-solved exactly as before.

# Code block 6: the optimum, priced at the nominal rate
TCnom(xy) = rnom * sum(f[i] * dgc(xy, P[i, :]) for i in eachindex(f))
xᵒn = Optim.minimizer(optimize(TCnom, [c0.LON, c0.LAT]))
Δn = TCnom(Cary) - TCnom(xᵒn)
prt(DataFrame(
    Site = ["optimum", "Cary, NC", "increase"],
    Cost = round.([TCnom(xᵒn), TCnom(Cary), Δn], digits = 0)))
      Site     Cost
───────────────────
   optimum   87,149
  Cary, NC  122,536
  increase   35,386

Three things in that table are worth checking rather than reading past, and each is one of the nine checks the course names.

# Code block 7: three checks on the top-down answer
prt(DataFrame(
    Check = ["Landmark: Cary re-priced equals what was spent",
             "Triangulate: the optimum is the same site",
             "Units: r_nom is \$/ton-mi"],
    Result = [round(TCnom(Cary) - TC₀, digits = 6),
              round(maximum(abs.(xᵒ .- xᵒn)), digits = 6),
              round(rnom, digits = 6)]))
                                           Check  Result
────────────────────────────────────────────────────────
  Landmark: Cary re-priced equals what was spent  0.0000
       Triangulate: the optimum is the same site  0.0000
                        Units: r_nom is $/ton-mi  0.2263

The first is an identity by construction: the rate was chosen so that the model reproduces the one number the firm already knows, which is what calibration means. The second is the substantive one. The two analyses put the facility in the same place, because the rate multiplies every term of the objective equally and so cannot move its minimum.

The nominal rate is $0.226 per ton-mile, and staying in Cary costs $35.4k per year more than the optimum: the same site and the same figure as part (a), from data that never mentioned a rate.

→ The rate sets the size of the answer, never its location. Bottom-up and top-down differ in where that one number comes from, and in nothing else.

Kudde Pillows

Kudde, a Swedish pillow manufacturer, is planning to build one or more plants in the U.S. to manufacture pillows and then distribute them to several major retailers located throughout the continental U.S. The shells for the pillows will be imported from Vietnam to each plant, 500,000 shells per 2-TEU container, where they will be stuffed with a locally sourced polyester filler and then shipped in 48 ft3 corrugated boxes, that each weigh 35 lb and contain 24 pillows, to the retailers’ DCs. Based on their experience in Europe, the total procurement and production cost to produce one million pillows is approximately $2 million, and, due to economies of scale, it has been estimated that it would cost around 2f^{0.85} million dollars to produce f million pillows at a new plant. Assuming that the annual demand will be one pillow per 100 people, how many plants should be constructed and where should they be located?

Kudde has no U.S. network. There is no spending to reverse, no incumbent plant whose cost can be re-evaluated, and no customer allocation on record. Top-down is not merely the weaker choice here; there is nothing for it to work on. So the rate is given, and everything else is built upward from the statement.

Example 2: Kudde’s U.S. plants

Determine the number and location of plants Kudde should build, given a transport rate of $2.5523 per ton-mile and a production cost subject to economies of scale.

Two quantities convert the statement into model data. A box of 24 pillows weighs 35 lb, which fixes the weight of a pillow; and one pillow per hundred people a year turns a census population into a tonnage.

# Code block 8: from the statement to demand, by ZIP code
ubox = 24                       # pillows per box
uwt  = 35 / ubox                # lb per pillow
udem = 1 / 100                  # pillows per person per year
rton = 2.5523                   # $/ton-mi, given
g    = 1.2                      # circuity

z = filter(r -> r.ISCUS && r.POP > 0, uszcta3())
P = hcat(z.LON, z.LAT)
f = float.(z.POP) .* (udem * uwt / 2000)     # ton/yr
(points = nrow(z), tons = round(sum(f), digits = 1))
(points = 882, tons = 2400.8)

The cost matrix is built once and never changes: only the fixed cost moves as the analysis proceeds.

# Code block 9: the variable cost of serving every ZIP from every site
D = g .* dgca(P, P, z.ALAND)
C = rton .* (f' .* D)
size(C)
(882, 882)

Lecture 2.4 obtained the fixed cost by fitting a line to the production costs of plants that already existed, and taking its intercept. Kudde has no plants to fit. What it has instead is a statement of how cost varies with size, and a line can be fitted to that curve just as well as to a scatter of observations: the intercept of the fit is the fixed component either way.

But a line fitted to a curve depends on the range it is fitted over,2 and the relevant range is the demand a single plant would see. That is precisely what the location problem is trying to determine. So the fixed cost and the plant count are solved alternately: assume a range, fit, solve, then re-fit over the demand the solution actually gives its largest plant.

# Code block 10: the scale curve, linearized over a demand range
f0   = 1e6 * uwt / 2000         # ton: what one million pillows weigh
TPC0 = 2e6                      # $/yr to produce that million
β    = 0.85                     # scale exponent, from the statement
TPC(x) = TPC0 .* (x ./ f0) .^ β

function fitk(fmax, fmin)
    x = collect(range(fmin, fmax; length = 2000))
    y = TPC(x)
    p = optimize(p -> sum((y .- (p[1] .+ p[2] .* x)) .^ 2), [0.0, 1.0])
    return Optim.minimizer(p)[1]   # the intercept is the fixed cost
end
fitk (generic function with 1 method)
# Code block 11: alternate between the fixed cost and the UFL
function alternate(C, f)
    tr = DataFrame(iter = Int[], fmax = Float64[], k = Float64[],
                   plants = Int[], TC = Float64[])
    fmin, fmax = minimum(f), sum(f)
    best = (TC = Inf, y = Int[], k = NaN)
    for i in 1:12
        k = fitk(fmax, fmin)
        y, TC, W = ufl(fill(k, size(C, 1)), C; verbose = false)
        push!(tr, (i, fmax, k, length(y), TC))
        TC < best.TC || break
        best = (TC = TC, y = y, k = k)
        fmax = maximum(W * f)   # the largest demand any plant now sees
    end
    return tr, best
end

tr, best = alternate(C, f)
prt(tr)
   iter      fmax           k  plants            TC
───────────────────────────────────────────────────
1     1  2,400.81  313,106.25       5  3,426,269.64
2     2    549.66   89,427.60       9  2,078,919.51
3     3    542.52   88,438.80       9  2,070,020.33
4     4    542.52   88,438.80       9  2,070,020.33

The first row is the analysis a reader would do without the alternation: assume one plant serves the country, and the implied fixed cost is large enough that only five are opened. Each subsequent row re-fits over a demand no plant actually exceeds, the fixed cost falls, and more plants become worth opening.

The trace is printed rather than the answer alone because it carries a Nudge check for free. Down the table the fixed cost falls and the plant count rises, and that is the one direction the model is allowed to move: a cheaper facility can only make more of them worth opening. A row where both fell would mean the fixed cost was not reaching the objective, and a single printed answer would have hidden it.

# Code block 12: where the plants go
prt(DataFrame(Plant = 1:length(best.y),
              Site = [lonlat2loc(P[i, :], usplace()).desc
                      for i in best.y]))
  Plant                                  Site
─────────────────────────────────────────────
      1             2.8 mi SE of Fircrest, WA
      2         1.7 mi SE of Altadena CDP, CA
      3                 in New Vernon CDP, NJ
      4           6.7 mi N of Gainesville, GA
      5  2.3 mi E of Bedford Park village, IL
      6           1.2 mi W of Wahneta CDP, FL
      7                      in Alamo CDP, CA
      8              8.5 mi SW of Elkhart, TX
      9  3.6 mi S of North Washington CDP, CO

9 plants, at a fixed cost of $88.4k a year each and a total cost of $2.07 million a year.

→ The fixed cost is not data here. It is an output of the same problem it is an input to, which is why the two are solved alternately rather than in sequence.

3. Popco Bottling Company

Popco is a century-old bottling company on the West Coast. It began with a single plant and grew by buying others, and a hundred years of that left it with 42 bottling plants, currently spread across the western U.S.. The question the company brings is whether 42 is the right number: whether profitability could be improved by reducing or adding plants, that is, whether some plants should be closed, or new ones opened, if the saving would be significant.

One thing about the units is worth fixing before the steps begin, because it is the most common place to lose the thread. The units of demand are tons. Population enters at exactly one step, to split a facility’s tonnage among the zones it serves in proportion to the people in them, and it appears nowhere else: not in the rate, not in the cost matrix, not in the objective. Every quantity downstream of that split is a tonnage.

Example 3: Popco’s 42 bottling plants

Determine the number and location of bottling plants that minimize Popco’s production and distribution cost, and the change against the 42 it operates.

Popco does its own distribution, so every plant knows what it spent last year delivering to its customers, and a century of accounts records what each plant produced and what that production cost. What the company will not hand over is customer order data. That combination is what makes this a top-down problem, and it is the ordinary case rather than the exception.

The following representative information is available for each current plant i \in N, four quantities that are the whole of the input, where

xy_i
= location of plant i
f^{DC}_i
= aggregate annual production, in tons
TPC_i
= total production and procurement cost, in dollars per year
TDC_i
= total distribution cost, in dollars per year.

The customers are not households. They are convenience stores, grocery chains and their distribution centers, and the quantities are carried in tons even though a bottler thinks in volume, because tons keep the units of the rate simple.

Because Popco’s plants are bottling plants, they are monetarily weight gaining: syrup arrives concentrated and leaves as filled bottles, so a truckload in becomes many truckloads out. The inbound procurement cost is therefore small beside the outbound distribution cost, and it barely depends on where the plant stands. The UFL can ignore it, which is what lets TPC_i and TDC_i be treated separately above.

That is an Assumptions check, and it is made here rather than at the end because it is what justifies the model that follows. Two assumptions carry the whole analysis. Weight gain is the first: it is what justifies dropping inbound procurement from the objective, and a plant that shipped out less than it took in would need the inbound leg evaluated. The second is that the fixed cost k is the same at every candidate site, which is what lets one number stand for every plant and every candidate location. The second is false in any real engagement: land and labor do not cost the same across the West. It is adopted because the alternative is a site-by-site cost estimate nobody has, and the place to say so is here rather than in a footnote to the answer.

# Code block 13: the 42 plants, as the company keeps them
DC = DataFrame(CSV.File("data/PopcoData.csv"))
rename!(DC, :DEMAND => :fDC, :PROD_COST => :TPC, :DIST_COST => :TDC)
nDC = nrow(DC)
prt(first(DC, 3))
     LAT      LON     fDC         TPC         TDC
─────────────────────────────────────────────────
1  44.06  -103.19  53,489  18,648,080  11,884,204
2  45.77  -108.52  10,920   5,866,471   1,198,643
3  45.97  -112.51  46,607  15,259,325  11,326,980

Fig. 3 is what the company has, and it is the only picture available before any analysis is done.

Show the code that draws this map
# Two colors carry the whole of Sec. 3: red is what Popco has.
incumbent = colorant"#c4342b"      # what Popco has today
chosen    = colorant"#1a7a3e"      # what the analysis recommends

function popcomap(; w = 640, h = 470)
    fig, ax = makemap(DC.LON, DC.LAT; xexpand = 0.04, yexpand = 0.04)
    resize!(fig, w, h)
    return fig, ax
end

fig, ax = popcomap()
scatter!(ax, DC.LON, DC.LAT; markersize = 10, color = :transparent,
         strokecolor = incumbent, strokewidth = 1.5)
ax.title = "$nDC bottling plants"
ax.titlesize = 15
fig
Figure 3: Popco’s 42 bottling plants. A century of acquisitions, and nobody has looked at the whole of it at once.

Before starting the procedure, here at a high level are the steps it takes.

  1. Regress total production cost on output across the plants, and keep the intercept.
  2. Allocate every ZIP code to its nearest plant, and drop the ones no plant reaches within a day’s round trip.
  3. Divide each plant’s tonnage among the ZIP codes it serves, in proportion to their populations.
  4. Divide what the network spent on distribution by the ton-miles it delivered, to determine the nominal transportation rate.
  5. Evaluate every candidate site against every customer at that rate, with the incumbent plants among the candidates.
  6. Run the UFL on the fixed cost and the cost matrix, and compare the result with the original network.

Step 1, the fixed cost. With 42 plants it is the regression of lecture 2.4, and the intercept is what is kept.

# Code block 14: regress production cost on output, keep the intercept
ŷ(p, x) = p[1] .+ p[2] .* x
loss(p) = sum((DC.TPC .- ŷ(p, DC.fDC)) .^ 2)
k, cp = Optim.minimizer(optimize(loss, [0.0, 1.0]))
(k = round(k), cp = round(cp, digits = 2))
(k = 1.0603653e7, cp = 126.15)

Fig. 4 is the fit, and the quantity the lecture wants is where the line meets the axis rather than the line itself.

Show the code that draws this figure
# The fit is neutral: the data and the intercept carry the color.
fitline = colorant"#4d565f"        # the least-squares fit
fig = Figure(size = (640, 400))
ax = Axis(fig[1, 1]; xlabel = "Annual production (ton/yr)",
          ylabel = "Production and procurement cost (\$/yr)",
          title = "The intercept is the fixed cost", titlesize = 15)
scatter!(ax, DC.fDC, DC.TPC; markersize = 9, color = incumbent)
xs = [0.0, maximum(DC.fDC)]
lines!(ax, xs, ŷ([k, cp], xs); color = fitline, linewidth = 1.8)
scatter!(ax, [0.0], [k]; markersize = 13, color = chosen)
text!(ax, 0.0, k; text = "  k", align = (:left, :bottom),
      fontsize = 16, color = chosen)
# The intercept IS what this figure exists to show, so the axis starts at
# zero: with Makie's default padding the y-axis sits left of x = 0 and k
# floats inside the plot instead of sitting on the axis.
xlims!(ax, 0, 1.02 * maximum(DC.fDC))
ylims!(ax, 0, 1.05 * maximum(DC.TPC))
fig
Figure 4: Total production and procurement cost against annual output, one point per plant, with the least-squares line. The intercept is the fixed cost the UFL uses; the slope is discarded, because a constant cost per ton is the same wherever the plant stands.

It is worth asking why the customer list is not simply requested, since Popco plainly has one. It could be, and the answer would come back in about a month, because each plant holds its own customer data and headquarters holds none of it. Some plants would answer and some would not, so the set might arrive covering thirty-five of the forty-two. What did arrive would be in whatever units each plant keeps: cases at one, tons at another, liters at a third. That is another week of cleaning before any of it can be used.

So the request is deliberately small: what each plant produced, and what that production and its distribution cost. Nothing else. Everything downstream of that is supplied by the analyst, and the substitution that makes it possible is that these plants ultimately serve people, so population is a usable proxy for demand.

Step 2, the market. Every three-digit ZIP goes to its nearest plant, and a plant serves nothing beyond 200 road miles. That number is not a round figure chosen for convenience. Popco supplies its own trucks, and those trucks must return to the plant at the end of the day, so 200 miles out is a 400-mile round trip, which is about the maximum a truck can do in a day. Beyond that a ZIP is not feasible to serve at all, even one holding a customer as large as a Walmart distribution center. A firm using common carriage would have no such screen.

Tracing it back that far is a Source check. The 200 miles is not a rule of the method and not a convention of the course: it is a fact about Popco’s own fleet, and the check is asking where the number came from before letting it decide which customers exist. A parameter that survives that question can be defended to a client and moved when the fleet changes; one that does not is a constant somebody once typed.

Eq. 1 is the market of plant i: every ZIP code whose nearest plant is i, and which lies within the screen:

M_i = \bigl\{\, j : \arg\min_h d_{hj} = i \ \text{ and } \ d_{ij} \le d_{\max} \,\bigr\}, \qquad d_{\max} = 200 \text{ mi} \tag{1}

The customers the analysis works on are M = \bigcup_{i \in N} M_i, which is smaller than the set of ZIP codes the country has. Running the screen over the whole country rather than the western states costs nothing and saves a judgment: the 200 miles removes the east automatically.

The screen also answers a question Sec. 2 left open. Circuity cancels out of a rate that is calibrated against the same distances it will later evaluate, which is why the top-down rate needs no circuity factor. It does not cancel out of a constraint. Two hundred miles is a distance a truck actually drives, so the comparison has to be made in road miles, so the circuity factor appears in the screen below even though it is absent from Eq. 3.

# Code block 15: allocate ZIP codes to plants, within range
z3 = filter(r -> r.ISCUS && r.POP > 0, uszcta3())
XDC = hcat(DC.LON, DC.LAT)
Dpz = dists(XDC, hcat(z3.LON, z3.LAT), :mi)
z3.IDX = [(i = argmin(@view Dpz[:, j]);
           Dpz[i, j] * 1.2 <= 200 ? i : 0) for j in 1:nrow(z3)]
nall = nrow(z3)
zout = filter(r -> r.IDX == 0, z3)   # beyond every plant's range
filter!(r -> r.IDX != 0, z3)
(candidates = nall, served = nrow(z3), unreachable = nrow(zout))
(candidates = 882, served = 162, unreachable = 720)

Fig. 5 is the market that screen defines, and the crosses are as much of the answer as the dots. The map is framed on the plants, so most of what the screen rejects lies off it entirely: the whole country east of these states is unreachable and takes no part in the analysis.

Show the code that draws this map
# A cross rather than a fainter dot: out of range is a rejection, not
# a weaker version of being served. Red stays what Popco has, here as
# in the first map, so the ZIP codes take the neutral and the marker
# shape carries the in-range/out-of-range distinction on its own.
served = colorant"#4d565f"         # a plant reaches these
beyond = colorant"#a9b0b6"         # no plant reaches these
fig, ax = popcomap()
scatter!(ax, zout.LON, zout.LAT; marker = :xcross, markersize = 6,
         color = beyond)
scatter!(ax, z3.LON, z3.LAT; markersize = 4.5, color = served)
scatter!(ax, DC.LON, DC.LAT; markersize = 10, color = :transparent,
         strokecolor = incumbent, strokewidth = 1.5)
ax.title = "$(nrow(z3)) of $nall ZIP codes are within range"
ax.titlesize = 14
fig
Figure 5: The market the 200-mile screen defines. Dots are ZIP codes a plant can serve; crosses are ZIP codes no plant reaches within a day’s round trip, and they take no part in the analysis.

Step 3, the split. Each plant’s tonnage is divided among the ZIP codes it serves in proportion to their populations:

f_{j \in M_i} = f^{DC}_i \, \frac{q_j}{\sum_{h \in M_i} q_h} \tag{2}

where

q_j
= population of existing facility j.

So if a plant’s market holds half a million people and one of its ZIP codes holds a hundred thousand, that ZIP code is allocated a fifth of the plant’s tonnage. This is the one step in which population appears at all.

# Code block 16: population-proportional demand, and the Balance check
z3 = transform(groupby(z3, :IDX), :POP => sum => :DCPOP)
DC.IDX = 1:nDC
z3 = leftjoin(z3, DC[!, [:IDX, :fDC]], on = :IDX)
z3.f = z3.fDC .* z3.POP ./ z3.DCPOP
(tons_at_plants = sum(DC.fDC), tons_at_zips = sum(z3.f),
 difference = sum(z3.f) - sum(DC.fDC))
(tons_at_plants = 5218228, tons_at_zips = 5.218228e6, difference = 0.0)

The third figure is a Balance check: splitting a total cannot change it, so the parts must sum to the whole. It is worth printing rather than assuming, because the split divides and division does not always return exactly what it started with.

Step 4, the nominal rate. What the network spent, divided by the ton-miles it delivered. This is the keystone of the whole lecture:

r_\text{nom} = \frac{\sum_{i \in N} TDC_i} {\sum_{i \in N} \sum_{j \in M_i} f_j \, d^a_{ij}} \tag{3}

where

d^a_{ij}
= area-adjusted distance from plant i to ZIP code j, in miles.

The adjustment matters here and nowhere more: a ZIP code served by a plant standing inside it is not served at zero cost, and without the floor of lecture 2.5 the denominator would be too small and the rate too high.

Reading Eq. 3 as a Units check is the cheapest guard the lecture has. Dollars per year over ton-miles per year gives dollars per ton-mile, which is what a transport rate is. The failure the check catches is a denominator that is only miles, because the demand weight was left out of the sum: the quotient is then dollars per mile, a different quantity of a different size, and it multiplies a ton-mile cost matrix downstream without complaint from anything but the units.

# Code block 17: back out the rate, and re-price the incumbent network
XZ = hcat(z3.LON, z3.LAT)
dsv = [dgca(XDC[z3.IDX[j]:z3.IDX[j], :], XZ[j:j, :], [z3.ALAND[j]])[1]
       for j in 1:nrow(z3)]
tonmi = sum(z3.f .* dsv)
rnomP = sum(DC.TDC) / tonmi
(rnom = round(rnomP, digits = 4), spent = sum(DC.TDC),
 repriced = round(rnomP * tonmi))
(rnom = 2.5597, spent = 449366635, repriced = 4.49366635e8)

The last two are a Landmark check, and this one is worth dwelling on. Evaluating the network Popco already runs, at the rate just derived from it, has to return the amount it actually spent. It does so by construction, which is what makes it useful: if the customers used in step 2 are not the customers the spending covered, the two will not agree, and nothing else in the analysis would have revealed it.

Step 5, the cost matrix. The candidate sites are every demand point and every plant Popco already has, so keeping a plant, moving it slightly and opening somewhere new are all available to the solver:

\mathbf{C} = \bigl[\, c_{ij} \,\bigr] = \bigl[\, r_\text{nom} f_j d^a_{ij} \,\bigr], \qquad i \in M \cup N, \ \ j \in M \tag{4}

The row index runs over M \cup N and the column index over M alone, which is the asymmetry that lets an incumbent plant be kept: a plant is a candidate site without being a customer.

# Code block 18: candidate sites are ZIP centroids plus incumbent plants
NF = vcat(XZ, XDC)
CP = rnomP .* (z3.f' .* dgca(NF, XZ, z3.ALAND))
size(CP)
(204, 162)

Step 6, the UFL. Before reading the answer, write down a guess: how many of the 42 plants survive, and does distribution cost rise or fall. That is a Prior, and it is worth spending here because the answer is genuinely hard to anticipate. Most people guess a mild consolidation and a saving in freight. One of those is wrong, and committing to the guess first is what makes the surprise informative rather than invisible.

# Code block 19: solve, and read the result against what Popco has
y, TC, W = ufl(k, CP)
TDCnew = TC - k * length(y)
TCorig = k * nDC + sum(DC.TDC)
prt(DataFrame(Case = ["existing", "recommended"],
              Plants = [nDC, length(y)],
              TDC = round.([sum(DC.TDC), TDCnew]),
              TC = round.([TCorig, TC])))
  Add: 7.512802738280858e8
 Xchg: 7.278099217481401e8
  Add: 7.278099217481401e8
 Drop: 7.260353173335743e8
 Xchg: 7.260353173335743e8
         Case  Plants          TDC           TC
───────────────────────────────────────────────
     existing      42  449,366,635  894,720,041
  recommended      24  471,547,657  726,035,317

Distribution cost goes up. That is the first thing a client will see, and it is worth stopping on before reading the second column: the analysis was commissioned to reduce cost and one of the costs has risen. What would be said to them?

The rise is not a defect, it is arithmetic. Serving the same demand from roughly half as many plants means shipping further, so distribution cost had to increase. What was bought with it is fixed cost: the variable production cost c_p is the same per ton wherever a plant stands, so it never entered the comparison, and what remains is one intercept k per plant. Halving the plants halves that, and the saving more than offsets the extra freight. The result is a rise of 4.9% in distribution against a fall of 18.9% in the total.

So the saving is not in the trucks. It is in the plants that no longer have to be paid for, and a reader who takes the first column alone will get the sign of the whole analysis wrong.

Two checks belong here, where a number has just been produced by a heuristic that carries no guarantee. The first is a Bounds check, and the bound is free: keeping every plant Popco already has is a feasible solution to the same UFL, so its cost is a ceiling the recommendation cannot exceed. A heuristic that returns more than the incumbent network costs has failed, and the comparison costs one line.

# Code block 20: Bounds, against the network Popco already runs
inc = nrow(z3) .+ (1:nDC)        # the incumbents, as candidate sites
TCceil = k * nDC + sum(minimum(CP[inc, :], dims = 1))
prt(DataFrame(Case = ["ceiling: keep all $nDC", "recommended"],
              TC = round.([TCceil, TC])))
                  Case           TC
───────────────────────────────────
  ceiling: keep all 42  894,720,041
           recommended  726,035,317

The ceiling lands on the existing network’s cost to the dollar, which is worth noticing rather than passing over. It was computed a different way: the table above added one intercept per plant to what Popco actually spent, while this one lets the solver allocate every ZIP code to whichever of the 42 is cheapest. The two agree because evaluating by cost and evaluating by distance rank the same sites in the same order, so the bound is tight and the Landmark check of step 4 holds at the level of the whole objective, not only the rate.

The second is a Nudge, and the direction is predictable before it runs, which is what makes it a check rather than an experiment. A plant that costs more to have should mean fewer plants, never more.

# Code block 21: Nudge, on the one parameter the answer turns on
nopen(kmul) = length(first(ufl(kmul * k, CP; verbose = false)))
prt(DataFrame(Fixed_cost = ["k/2", "k", "2k"],
              Plants = [nopen(0.5), nopen(1.0), nopen(2.0)]))
  Fixed_cost  Plants
────────────────────
         k/2      39
           k      24
          2k      16

Halving the fixed cost opens more plants and doubling it opens fewer, which is the monotone response the model must have. A run that moved the other way would mean the fixed cost had been wired into the objective with the wrong sign, and no amount of staring at a map would have shown it.

# Code block 22: what kind of site each recommendation is
kept = count(i -> i > nrow(z3), y)
prt(DataFrame(
    Kind = ["an incumbent plant location", "a new site"],
    Sites = [kept, length(y) - kept]))
                         Kind  Sites
────────────────────────────────────
  an incumbent plant location      9
                   a new site     15

Fig. 6 puts the recommendation against the network it replaces, and three categories are visible in it. A green ring sitting on a red dot is an existing plant the solution keeps. A ring standing alone is a genuinely new site, in country the current network does not reach well.

The third category is the one worth arguing with: a ring a few miles off a red dot. That is not a plant being kept. It is a ZIP centroid that happens to lie near an existing plant, so the recommendation is to close a working plant and build a new one a short distance away, which no client would accept and no analyst should propose. The model has no way of knowing the difference, because every candidate site is evaluated only on distance and fixed cost.

The fix is a second screen in the other direction. The 200 miles is an upper bound on how far a plant may serve; add a lower bound on how close a new site may be to an existing plant, perhaps 50 miles, and drop those candidates. The solver then chooses between keeping a plant and building somewhere genuinely different, which is the choice actually on offer.

The results are read off Fig. 6. The original 42 plants, in red, are reduced to the 24 green circles; total distribution cost stays about the same, while total cost decreases.

Show the code that draws this map
# Same red and green rings as the first map, so the two read together.
fig, ax = popcomap()
scatter!(ax, DC.LON, DC.LAT; markersize = 8, color = incumbent)
scatter!(ax, NF[y, 1], NF[y, 2]; markersize = 12, color = :transparent,
         strokecolor = chosen, strokewidth = 1.7)
ax.title = "$nDC existing plants, $(length(y)) recommended sites"
ax.titlesize = 15
fig
Figure 6: The recommendation against the incumbent network. Red dots are Popco’s existing plants; green rings are the sites the UFL selects. Some rings sit on a dot, some sit just beside one, and some stand alone.

42 plants become 24, distribution cost rises 4.9%, and total cost falls 18.9%.

→ A week’s work, on data the client already had, and no order history. The number it produces is not the decision; it is the one input the company cannot produce for itself.

Capacity and expansion

The recommendation is not the decision, and two of the things it does not determine are the ones a client asks first:

  • Question: If an original plant is kept, will it have sufficient capacity (since fewer plants)

  • Question: If Popco was planning to expand to serve the entire continental U.S., how could the current results be utilized

The first is a real gap in the model rather than a caution about it. The uncapacitated facility location problem is uncapacitated by construction: it places facilities as though each could make whatever is allocated to it, and nothing in Eq. 4 or in the solve says otherwise. Closing 42 plants down to 24 over the same total tonnage therefore hands each surviving site more work, and the model has no opinion about whether it can be done. Some of these plants would be asked to double their production rate, and a plant already running two or three shifts cannot. That is a capacity constraint, and the UFL carries no constraints at all, which is what the mixed-integer formulation of lecture 2.7 changes.

# Code block 23: what each site would have to produce
load = W * z3.f                       # ton/yr allocated to each open site
kept = [i for i in y if i > nrow(z3)]
prt(DataFrame(
    Network = ["existing, $nDC plants",
               "recommended, $(length(y)) sites"],
    Mean = round.([sum(DC.fDC) / nDC, sum(load[y]) / length(y)]),
    Largest = round.([maximum(DC.fDC), maximum(load[y])])))
                Network     Mean  Largest
─────────────────────────────────────────
    existing, 42 plants  124,244  526,208
  recommended, 24 sites  217,426  692,487

The mean rises by a factor of 1.75, which is the arithmetic and is not the interesting part. Two things in that table matter more. The largest site would have to produce more than any plant Popco runs today, so the recommendation quietly asks the company to operate at a scale it has never operated at. And the mean hides the spread: the plants that survive do not grow uniformly.

# Code block 24: how much each kept plant would grow
prt(DataFrame(
    Plant = [i - nrow(z3) for i in kept],
    Now = round.([DC.fDC[i - nrow(z3)] for i in kept]),
    Planned = round.([load[i] for i in kept]),
    Ratio = round.([load[i] / DC.fDC[i - nrow(z3)] for i in kept],
                   digits = 2)))
   Plant      Now  Planned  Ratio
─────────────────────────────────
1     19  197,496  197,496   1.00
2     38  112,759  162,133   1.44
3      1   53,489   57,763   1.08
4     12   60,138   45,911   0.76
5     17  375,863  456,469   1.21
6     37   40,593   98,196   2.42
7      6  526,208  540,301   1.03
8     31  329,160  462,679   1.41
9     10   91,699   90,785   0.99

The kept plants do not move together at all: the ratio of planned to current output runs from 0.76 to 2.42. So the question is asked per site and not of the average, and the sites it flags are not the ones the average would suggest.

The first evidence that it may be answerable is already on the page. Fig. 4 scatters the plants across a wide span of annual output, and a firm whose plants already run at very different scales is a firm whose plant capacity is flexible. Extra lines and second shifts are how that flexibility is usually realized. Had every plant sat at the same output, the span would have said the opposite, and a fixed capacity would have to be modeled rather than assumed away.

Lecture 1.3 Sec. 5 turns that from an impression into a test. A plant’s nameplate rate is optimistic: the throughput-feasible minimum inflates the arrival rate by yield loss and the process time by availability, so a site that looks adequate against its rated output can still be infeasible. Running it against the table above produces something more useful than a yes or a no. It produces pointed questions: whether the plant asked to more than double its current output can be taken there at all, and what it would cost. Those are questions a client can answer in a phone call, which a request for all their customer data is not.

Where the constraint genuinely cannot be met, the answer is not to re-run the UFL and hope. It is to say what capacity each site may have and let the model respect it, which is a mixed-integer formulation and is lecture 2.7.

The second question has an answer that costs surprisingly little. A national Popco needs no new method: take every three-digit ZIP in the continental United States as a candidate site, keep the western plants among the candidates, let each candidate’s market be the ZIPs within 200 miles of it, carry the same fixed cost, and solve the same UFL. The rate transfers on the assumption that a firm whose fleet, product and service pattern are unchanged buys miles at about the same rate in a new region, which is the transfer case of Sec. 2.

Two cautions come with it, and both are the kind that should be said aloud rather than discovered later. The eastern market may simply not behave like the western one, so this is a first cut and not a plan. And a national run raises a question the western one never had to face: whether every ZIP code should be served at all, since some are sparse enough that covering them cannot pay. The model as written assumes they all must be.

The whole of this analysis takes less than a week for a real company, and that speed is not a boast about efficiency. It is what makes the result useful. Presenting a number is what causes a client to say that the number is wrong, and to explain why: some constraint nobody thought to mention, a plant that cannot be closed for reasons that have nothing to do with cost, a customer who must be served from a particular site. That information is not withheld. It is latent, and a concrete recommendation is what brings it out.

So the first analysis is deliberately the simplest one that produces a number, and its value is partly in being wrong in ways the client can correct. An analysis that takes a month to produce arrives after the funding has moved and after the people who could have corrected it have stopped thinking about the question.

Eight of the nine checks have been used on the way here, each beside the number it had something to say about: Landmark twice, on the two calibrated rates; Assumptions at the setup, where the model is justified; Source at the 200-mile screen; Balance at the population split; Units on Eq. 3; Prior before the solve; Bounds and Nudge on the result. Which checks apply is a property of the problem, and a network design supports a different set from a distance calculation or a rate lookup.

Triangulate is the one this problem does not support, and saying so is part of the answer. There is no independent route to the number: a second method would need a second rate, and the rate is what the analysis exists to produce. Where a check has nothing to bite on, the honest report says so rather than manufacturing agreement.

4. Two plants on one road

The same six steps run at the scale a real engagement has are easiest to check at a size that fits on paper. They are run again here on an instance small enough to do by hand, where every quantity can be verified against the one before it. The zones used in this example are meant to be one-dimensional versions of an actual two-dimensional ZIP code area.

Example 4: Two plants on one road

Determine the change in total cost from closing both of a firm’s plants and opening a single plant at zone 2.

A firm has two plants, A and B, serving four zones that lie along one road. Everything is read off the road: plant A sits at mile 0, the zones at miles 30, 90, 160 and 420, and plant B at mile 200.

# Code block 25: the instance, read off one road
pos = (A = 0.0, Z1 = 30.0, Z2 = 90.0, Z3 = 160.0, B = 200.0, Z4 = 420.0)
zone, plant = [:Z1, :Z2, :Z3, :Z4], [:A, :B]

output = [300.0, 200.0]                    # ton/yr produced at A and B
prodcost = [460_000.0, 340_000.0]       # $/yr production + procurement
distcost = [162_000.0, 80_000.0]        # $/yr outbound distribution
popn = [120_000.0, 80_000.0, 100_000.0, 60_000.0]    # people
area = [1_200.0, 900.0, 1_500.0, 800.0]              # mi^2

d(i, j) = abs(pos[i] - pos[j])
Dpz = [d(p, z) for p in plant, z in zone]
prt(DataFrame(Plant = string.(plant), Z1 = Dpz[:, 1], Z2 = Dpz[:, 2],
              Z3 = Dpz[:, 3], Z4 = Dpz[:, 4]))
  Plant   Z1   Z2   Z3   Z4
───────────────────────────
      A   30   90  160  420
      B  170  110   40  220

Step 1. With two plants the regression of Sec. 3 is a line through two points, so the intercept can be read off directly and no fitting is needed.

# Code block 26: the fixed cost, from two points
cp = (prodcost[1] - prodcost[2]) / (output[1] - output[2])
k = prodcost[2] - cp * output[2]
(cp = cp, k = k)
(cp = 1200.0, k = 100000.0)

Step 2. Each zone goes to its nearer plant, and anything beyond 200 miles is dropped.

# Code block 27: the market, and what falls outside it
nearest = [argmin(@view Dpz[:, j]) for j in eachindex(zone)]
served = [Dpz[nearest[j], j] <= 200 for j in eachindex(zone)]
prt(DataFrame(Zone = string.(zone), Plant = string.(plant[nearest]),
              Miles = [Dpz[nearest[j], j] for j in eachindex(zone)],
              Served = [s ? "yes" : "no" for s in served]))
  Zone  Plant  Miles  Served
────────────────────────────
    Z1      A     30     yes
    Z2      A     90     yes
    Z3      B     40     yes
    Z4      B    220      no

Step 3. Each plant’s tonnage is split among its zones in proportion to population, by Eq. 2.

# Code block 28: population-proportional demand
f = zeros(length(zone))
for i in eachindex(plant)
    mkt = findall(j -> served[j] && nearest[j] == i, eachindex(zone))
    f[mkt] .= output[i] .* popn[mkt] ./ sum(popn[mkt])
end
(tons_at_plants = sum(output), tons_at_zones = sum(f))
(tons_at_plants = 500.0, tons_at_zones = 500.0)

The second figure is there to be compared with the first. Splitting a total among zones cannot change it, so printing both is a Balance check on the one step where tonnage changes hands. It costs a line, and it guards the step most likely to lose some: a zone dropped by the screen must have had its tonnage reassigned to the zones that remain, never simply discarded.

Step 4. The nominal rate is what was spent over the ton-miles delivered, by Eq. 3, with the area floor of lecture 2.5 standing in wherever a zone is served from inside itself.

# Code block 29: the nominal rate
floorz = (2 / 3) .* sqrt.(area ./ π)
da(dist, j) = max(dist, floorz[j])
tonmi = sum(f[j] * da(Dpz[nearest[j], j], j)
            for j in eachindex(zone) if served[j])
rnom = sum(distcost) / tonmi
(ton_miles = tonmi, spent = sum(distcost), rnom = rnom)
(ton_miles = 24200.0, spent = 242000.0, rnom = 10.0)

All three quantities are printed so the division can be read as a Units check: dollars a year over ton-miles a year is dollars per ton-mile, and a rate that came out in dollars per mile would mean the demand weight had been left out of the denominator. The rate itself is high for freight, which is the second reason to print it: these are short hauls carrying small tonnages, and a number that looks wrong is worth explaining rather than passing over.

Steps 5 and 6. With four zones there is no need for a solver: the question asks what one specific configuration costs, so it is evaluated directly against the network the firm runs today.

# Code block 30: closing both plants for one at zone 2
Dzz = [d(a, b) for a in zone, b in zone]
served_idx = findall(served)
TCnew = k + rnom * sum(f[j] * da(Dzz[2, j], j) for j in served_idx)
TCorig = length(plant) * k + sum(distcost)
prt(DataFrame(Case = ["two plants, as today", "one plant at zone 2"],
              TC = round.([TCorig, TCnew])))
                  Case       TC
───────────────────────────────
  two plants, as today  442,000
   one plant at zone 2  361,541

Everything this example does happens along one road, so the whole of it fits in Fig. 7: the instance, the screen that decides which zones exist at all, the tonnage each surviving zone is given, and the single site that replaces both plants.

Show the code that draws this figure
# Four rows on one shared mile axis. The road is one-dimensional, so
# no map is needed and none is used: Makie primitives on a common x.
# Red is what the firm has and green what the analysis picks, the same as
# every other figure in this lecture.
zonec = colorant"#4d565f"          # a zone inside the screen
outc  = colorant"#a9b0b6"          # a zone the screen drops
XMAX  = 520                        # the reach labels sit past Z4

roadfig = Figure(size = (660, 470))
titles = ["the road", "allocation, and the 200-mile screen",
          "tonnage given to each zone", "one site in place of two"]
axs = [Axis(roadfig[r, 1]; title = titles[r], titlesize = 13,
            titlealign = :left) for r in 1:4]
for (r, ax) in enumerate(axs)
    hideydecorations!(ax)
    hidespines!(ax)
    if r == 4
        ax.xlabel = "miles along the road"
        ax.xlabelsize = 13
    else
        hidexdecorations!(ax; grid = false)
    end
    xlims!(ax, -22, XMAX)   # room for the mile-0 and reach labels
    lines!(ax, [0, XMAX], [0, 0]; color = (:black, 0.2), linewidth = 1)
end

function tag!(ax, x, y, s, c; al = (:center, :bottom))
    text!(ax, x, y; text = s, color = c, fontsize = 12, align = al)
end

# Row 1: where everything is, and nothing else.
for p in plant
    scatter!(axs[1], [pos[p]], [0]; markersize = 13, color = incumbent)
    tag!(axs[1], pos[p], 0.15, "$p ($(Int(pos[p])))", incumbent)
end
for z in zone
    scatter!(axs[1], [pos[z]], [0]; markersize = 9, color = :transparent,
             strokecolor = zonec, strokewidth = 1.5)
    tag!(axs[1], pos[z], -0.55, "$z ($(Int(pos[z])))", zonec)
end
ylims!(axs[1], -1.0, 1.0)

# Row 2: the screen is a reach, so it is drawn as one. The bands sit
# at their own heights: drawn together they overlap from 0 to 200 and
# read as one bar, which hides the only thing the row is for.
for (i, p) in enumerate(plant)
    y = -0.55 - 0.45 * (i - 1)
    lo, hi = pos[p] - 200, pos[p] + 200
    lines!(axs[2], [max(lo, 0), hi], [y, y]; color = (incumbent, 0.30),
           linewidth = 8)
    lines!(axs[2], [hi, hi], [y - 0.14, y + 0.14]; color = incumbent,
           linewidth = 1.5)
    tag!(axs[2], hi + 6, y - 0.1, "$p reaches $(Int(hi))", incumbent;
         al = (:left, :bottom))
    scatter!(axs[2], [pos[p]], [0]; markersize = 13, color = incumbent)
end
for (j, z) in enumerate(zone)
    c = served[j] ? zonec : outc
    scatter!(axs[2], [pos[z]], [0]; markersize = 9, color = c,
             marker = served[j] ? :circle : :xcross)
    # Alternating heights: at 30 and 90 miles the labels collide.
    near = plant[nearest[j]]
    tag!(axs[2], pos[z], isodd(j) ? 0.16 : 0.46,
         served[j] ? "$z: $(Int(Dpz[nearest[j], j])) mi to $near" :
                     "$z: dropped", c)
end
ylims!(axs[2], -1.5, 1.25)

# Row 3: a bar per zone, so the split reads as a quantity.
for (j, z) in enumerate(zone)
    served[j] || continue
    lines!(axs[3], [pos[z], pos[z]], [0, f[j]];
           color = zonec, linewidth = 8)
    tag!(axs[3], pos[z], f[j] + 10, "$(Int(round(f[j]))) ton", zonec)
end
ylims!(axs[3], -25, 275)

# Row 4: the candidate, and what it would have to reach.
zc = pos[zone[2]]
for p in plant
    scatter!(axs[4], [pos[p]], [0]; markersize = 13, color = :transparent,
             strokecolor = (incumbent, 0.45), strokewidth = 1.5)
end
scatter!(axs[4], [zc], [0]; markersize = 15, color = chosen)
tag!(axs[4], zc, 0.18, "one plant at $(zone[2])", chosen)
reach = 0
for (j, z) in enumerate(zone)
    (served[j] && j != 2) || continue
    reach += 1
    y = -0.4 - 0.35 * (reach - 1)      # one line each, so none merge
    lines!(axs[4], [zc, pos[z]], [y, y];
           color = (chosen, 0.5), linewidth = 2)
    # Left-aligned at the far end, so the label runs away from the
    # site marker rather than across it.
    tag!(axs[4], min(zc, pos[z]) + 4, y + 0.05,
         "to $z, $(Int(abs(zc - pos[z]))) mi", chosen;
         al = (:left, :bottom))
end
ylims!(axs[4], -1.4, 1.0)
roadfig
Figure 7: The whole example on one road. Row 1 is the instance: plants at mile 0 and 200, zones at 30, 90, 160 and 420. Row 2 is the allocation, each zone joined to the nearer plant, with the band under each plant showing how far it can serve; Z4 lies outside both bands and leaves the problem. Row 3 is the tonnage each surviving zone receives, split by population. Row 4 is the single site that replaces both plants, with the distance it must cover to each zone.

One plant at zone 2 costs $80.5k a year less than the two the firm runs, a reduction of 18.2%.

Two things in that run are worth keeping. Zone 4 was dropped: it lies 220 miles from its nearer plant and so takes no part in the analysis at all, which is the screen doing its work. And the floor bound at zone 2, where a plant standing inside the zone would otherwise have distributed to it for nothing.

→ Consolidating and relocating are different questions. The UFL run over all six candidate sites opens two plants, at zones 1 and 3, for less than the single plant at zone 2 that the question asks about.

Endnotes

  1. R. K. Ahuja, T. L. Magnanti and J. B. Orlin, Network Flows: Theory, Algorithms, and Applications, Englewood Cliffs, NJ: Prentice-Hall, 1993.↩︎

  2. The intercept is an extrapolation to zero output, and zero output lies far outside any range of operating plants, so a slight bend inside the fitted range swings the line a long way by the time it reaches the axis. The effect is not sampling error and more observations do not remove it: it is the cost curve being concave, which is what economies of scale mean. It is starkest when the curve is given rather than observed, as it is here, since a stated scale curve of the form TPC(f) = TPC_0 (f/f_0)^\beta passes through the origin and so has no fixed cost at all. Whatever intercept a line through it reports is a property of where the line was fitted. That is why the fit range is stated whenever the curve is given, and why Sec. 3, fitting to plants that exist, takes the intercept once and does not iterate.↩︎