2. Loc 4: Allocation and ALA

ISE 754, Fall 2024

Package Used: No new packages used.

1. Allocation

Allocate NC 100K city population to Raleigh and Charlotte based on their closeness:

function dgc(xy₁, xy₂; unit=:mi)
    length(xy₁) == length(xy₂) == 2 || error("Inputs must have length 2.")
    unit in [:mi, :km] || error("Unit must be :mi or :km")

    Δx, Δy = xy₂[1] - xy₁[1], xy₂[2] - xy₁[2]
    a = sind(Δy / 2)^2 + cosd(xy₁[2]) * cosd(xy₂[2]) * sind(Δx / 2)^2
    2 * asin(min(sqrt(a), 1.0)) * (unit == :mi ? 3958.75 : 6371.00)
end

function name2lonlat(name, df)
    idx = findfirst(r -> r.NAME == name, eachrow(df))   # Find first match
    if idx === nothing
        error("'$name' not found in $df")
    end
    return collect(df[idx, [:LON, :LAT]])
end

name2lonlat (generic function with 1 method)

using Logjam.DataTools, DataFrames

df = filter(r -> (r.STFIP == st2fips(:NC)) && (r.POP > 100_000), usplace())
select!(df, :NAME, :LAT, :LON, :POP)

10×4 DataFrame

Row	NAME	LAT	LON	POP
	String	Float64	Float64	Int64
1	Cary	35.7814	-78.819	174721
2	Charlotte	35.209	-80.831	874579
3	Concord	35.3921	-80.6355	105240
4	Durham	35.98	-78.901	283506
5	Fayetteville	35.0828	-78.9735	208501
6	Greensboro	36.0951	-79.827	299035
7	High Point	35.9908	-79.9927	114059
8	Raleigh	35.8302	-78.6415	467665
9	Wilmington	34.2092	-77.8858	115451
10	Winston-Salem	36.1029	-80.2608	249545

P = hcat(df.LON, df.LAT)

10×2 Matrix{Float64}:
 -78.819   35.7814
 -80.831   35.209
 -80.6355  35.3921
 -78.901   35.98
 -78.9735  35.0828
 -79.827   36.0951
 -79.9927  35.9908
 -78.6415  35.8302
 -77.8858  34.2092
 -80.2608  36.1029

xyC = name2lonlat("Charlotte", df)
xyR = name2lonlat("Raleigh", df)
X = vcat(xyC', xyR')

2×2 Matrix{Float64}:
 -80.831   35.209
 -78.6415  35.8302

Dgc(X₁, X₂) = [dgc(i, j) for i in eachrow(X₁), j in eachrow(X₂)]

D = Dgc(X, P)

2×10 Matrix{Float64}:
 119.888     0.0    16.7769  120.81    …  71.661   130.39  180.976  69.5588
  10.5021  130.39  116.024    17.8347     76.4247    0.0   119.883  92.493

α = [argmin(c) for c in eachcol(D)]

10-element Vector{Int64}:
 2
 1
 1
 2
 2
 2
 1
 2
 2
 1

w = df.POP

10-element Vector{Int64}:
 174721
 874579
 105240
 283506
 208501
 299035
 114059
 467665
 115451
 249545

using SparseArrays

n, m = size(D)
W = sparse(α, 1:m, w, n, m)

2×10 SparseMatrixCSC{Int64, Int64} with 10 stored entries:
      ⋅  874579  105240       ⋅       ⋅  …  114059       ⋅       ⋅  249545
 174721       ⋅       ⋅  283506  208501          ⋅  467665  115451       ⋅

TC2 = sum(W .* D)

8.004697649736136e7

Add Greensboro: What would be the impact of adding Greensboro as potential allocation point?

X = vcat(X, name2lonlat("Greensboro", df)')

3×2 Matrix{Float64}:
 -80.831   35.209
 -78.6415  35.8302
 -79.827   36.0951

D = Dgc(X, P)

3×10 Matrix{Float64}:
 119.888     0.0      16.7769  …  71.661   130.39    180.976  69.5588
  10.5021  130.39    116.024      76.4247    0.0     119.883  92.493
  60.4097   83.2192   66.448      11.7341   68.7784  170.301  24.2256

α = [argmin(c) for c in eachcol(D)]
n, m = size(D)
TC3 = sum(sparse(α, 1:m, w, n, m) .* D)

4.133195732115901e7

println("Deduction in person-miles = ", 100*(TC2 - TC3)/TC2, "%")

Deduction in person-miles = 48.36537352223233%

2. Facility Location–Allocation Problem

Ex: I-40 Cities (1-D location problem)

d1(x1, x2) = length(x1) == length(x2) ? sum(abs.(x1 .- x2)) : 
    error("Inputs not same length.")
D1(X₁, X₂) = [d1(i, j) for i in eachrow(X₁), j in eachrow(X₂)]

P = [50 150 190 220 270 295 420]'
m = size(P, 1)
w = collect(1:m)
X = [100 300]'
n = size(X, 1)
D1(X, P)

2×7 Matrix{Int64}:
  50   50   90  120  170  195  320
 250  150  110   80   30    5  120

# First Approach: Integrated Formulation
using SparseArrays

function TCint(X)
    D = D1(X, P)
    α = [argmin(c) for c in eachcol(D)]
    n, m = size(D)
    return sum(sparse(α, 1:m, w, n, m) .* D)
end

TCint (generic function with 1 method)

using CairoMakie

xrng = 0:500
Z = [TCint([x1; x2]) for x1 in xrng, x2 in xrng]
contour(xrng, xrng, Z; levels = 100)

using Optim

Xᵒ = optimize(TCint, [0.; 200.]).minimizer

2-element Vector{Float64}:
 189.99999999849905
 294.9999999987741

Xᵒ = optimize(TCint, [200.; 500.]).minimizer

2-element Vector{Float64}:
 269.9999999550589
 419.99999999652727

# Second Approach: Alternating Location-Allocation (ALA) Procedure

TC(W, X) = sum(W .* D1(X, P))
alloc(X) = sparse([argmin(c) for c in eachcol(D1(X, P))], 1:m, w, n, m)
loc(W, X₀) = optimize(X -> TC(W, X), X₀).minimizer

function ala(Xᵒ)            # Input initial location
    TCᵒ, done = Inf, false
    while !done
        W = alloc(Xᵒ)       # Allocate
        X′ = loc(W, Xᵒ)     # Locate
        TC′ = TC(W, X′)
        println(TC′, X′)
        if TC′ < TCᵒ
            TCᵒ, Xᵒ = TC′, X′
        else
            done = true
        end
    end
end

ala (generic function with 1 method)

ala([0.; 200.])

1880.0000000155678[50.00000000602211, 270.0000000095456]
1715.0000000164055[149.99999999875615, 294.9999999848385]
1340.0000000094158[190.00000000182195, 295.00000000144297]
1340.0000000094158[190.00000000182195, 295.00000000144297]

ala([200.; 500.])

1050.000003411217[269.9999969104532, 419.9999999540471]
1050.0000000605023[269.99999996976754, 419.99999999567575]
1050.0000000122438[269.9999999911894, 419.99999999950955]
1050.0000000122438[269.9999999911894, 419.99999999950955]

Ex: NC & SC 10K+ Cities (2-D location problem)

using Logjam.DataTools, Logjam.MapTools
using GeoMakie, DataFrames

df = filter!(r -> any(i -> r.STFIP == i, st2fips.([:NC,:SC])) && 
    r.POP > 10_000, usplace())
select!(df, :NAME, :ST, :LAT, :LON, :POP)

fig, ax = makemap(df.LON, df.LAT; xexpand=0.1)
scatter!(ax, df.LON, df.LAT; marker='.', markersize=24, color=:red)
fig

P = hcat(df.LON, df.LAT)

155×2 Matrix{Float64}:
 -80.1921  35.36
 -78.9583  35.2638
 -82.5527  35.571
 -81.6652  36.2138
 -79.4683  36.0761
 -78.819   35.7814
 -79.0383  35.9259
 -80.831   35.209
 -80.6355  35.3921
 -78.901   35.98
 -76.2354  36.294
 -78.9735  35.0828
 -81.1864  35.2488
   ⋮       
 -82.4976  34.78
 -81.2339  33.9295
 -79.0151  33.778
 -81.1056  34.0511
 -81.1409  34.0462
 -82.2564  34.7285
 -79.0089  33.6874
 -80.1796  33.0021
 -82.3139  34.9143
 -81.0167  35.0346
 -82.3317  34.8827
 -81.0994  33.9901

using Random

function randX(P, n)
    (xmn, xmx), (ymn, ymx) = extrema(P, dims=1)
    return [xmn ymn] .+ rand(n, 2) .* [(xmx - xmn) (ymx - ymn)]
end

Random.seed!(8345)

X₀ = randX(P, 3)

3×2 Matrix{Float64}:
 -79.0625  35.8591
 -76.3702  36.299
 -79.2144  34.6258

n = 2
m = size(P, 1)
w = df.POP

TC(W, X) = sum(W .* Dgc(X, P))
alloc(X) = sparse([argmin(c) for c in eachcol(Dgc(X, P))], 1:m, w, n, m)
loc(W, X₀) = optimize(X -> TC(W, X), X₀).minimizer

function ala(Xᵒ)            # Input initial location
    TCᵒ, done = Inf, false
    while !done
        W = alloc(Xᵒ)       # Allocate
        X′ = loc(W, Xᵒ)     # Locate
        TC′ = TC(W, X′)
        println(TC′)
        if TC′ < TCᵒ
            TCᵒ, Xᵒ = TC′, X′
        else
            done = true
        end
    end
    return Xᵒ, TCᵒ 
end

@show X₀ = randX(P, n)
Xᵒ, TCᵒ = ala(X₀)
Xᵒ

X₀ = randX(P, n) = [-77.8709055268364 34.55741060286532; -76.98603711094563 33.71346377992808]
6.119252565942339e8
6.119252565942339e8

2×2 Matrix{Float64}:
 -80.2176  35.3065
 -98.1239  28.4809

@show X₀ = randX(P, n)
Xᵒ, TCᵒ = ala(X₀)
Xᵒ

X₀ = randX(P, n) = [-79.34104766328115 33.43515740590162; -79.41217429380765 33.28446258172177]
5.0895689005692464e8
4.9188933126056933e8
4.908051153938427e8
4.9076019385930514e8
4.90760193859305e8
4.90760193859305e8

2×2 Matrix{Float64}:
 -80.0565  35.6086
 -80.0569  33.0971

@show X₀ = randX(P, n)
Xᵒ, TCᵒ = ala(X₀)
Xᵒ

X₀ = randX(P, n) = [-81.14941461895812 35.32133853473314; -79.45311349760652 35.07362594110691]
4.2413052084111845e8
3.9351789995158535e8
3.935178998918117e8
3.935178998918115e8
3.935178998918115e8

2×2 Matrix{Float64}:
 -80.831   35.209
 -78.8084  35.775

scatter!(ax, Xᵒ[:, 1], Xᵒ[:, 2]; marker=:dtriangle, markersize=12, color=:black)
fig

function lonlat2name(X, D, name, st)
    for i = 1:size(D, 1)
        d = D[i, :]
        println("NF", i, ": ", round(minimum(d)), " miles from ", name[argmin(d)], ", ", st[argmin(d)])
    end
end

lonlat2name(Xᵒ, Dgc(Xᵒ, P), df.NAME, df.ST)

NF1: 0.0 miles from Charlotte, NC
NF2: 1.0 miles from Cary, NC

# Take the best of 'nruns' runs of ALA
nruns = 5
TCᵒ, Xᵒ = Inf, nothing
for i = 1:nruns
    println("\nRun ", i)
    X′, TC′ = ala(randX(P, n))
    if TC′ < TCᵒ
        Xᵒ, TCᵒ = X′, TC′
    end
end
println("\nTC = ", TCᵒ, " person-miles")
lonlat2name(Xᵒ, Dgc(Xᵒ, P), df.NAME, df.ST)

Run 1
4.914255798427955e8
4.215160072931946e8
3.9387212718630224e8
3.935178998510067e8
3.935178998510067e8

Run 2
4.228849214247457e8
3.974554709664992e8
3.9351693654424417e8
3.935169365442441e8
3.9351693654424393e8
3.9351693654424393e8

Run 3
4.9323565770943356e8
4.920951595691688e8
4.9072388041677475e8
4.877664908572087e8
4.8657635131324047e8
4.863563996193143e8
4.424586829287798e8
3.985090982174633e8
3.935290585734499e8
3.935169365442902e8
3.9351693654429007e8
3.9351693654429007e8

Run 4
5.313595083725199e8
4.93081054245647e8
4.396751987433205e8
3.985090982174644e8
3.935290585734522e8
3.9351693654424703e8
3.935169365442469e8
3.9351693654424685e8
3.9351693654424685e8

Run 5
6.119252565942339e8
6.119252565942338e8
6.119252565942338e8

TC = 3.9351693654424393e8 person-miles
NF1: 0.0 miles from Charlotte, NC
NF2: 0.0 miles from Cary, NC

Allocation constraint: Due to a preexising contract. the first three EFs must be allocated to NF1, and the next three EFs must be allocated to NF2.

# Re-define `alloc`:
function alloc(X)
    α = [argmin(c) for c in eachcol(Dgc(X, P))]
    α[1:3] .= 1
    α[4:6] .= 2
    return sparse(α, 1:m, w, n, m)
end

# Need to re-init `ala` so that reads the new `alloc`
function ala(Xᵒ)            # Input initial location
    TCᵒ, done = Inf, false
    while !done
        W = alloc(Xᵒ)
        X′ = loc(W, Xᵒ)
        TC′ = TC(W, X′)
        println(TC′)
        if TC′ < TCᵒ
            TCᵒ, Xᵒ = TC′, X′
        else
            done = true
        end
    end
    return Xᵒ, TCᵒ 
end

Xᵒ, TCᵒ = ala(randX(P, n))
lonlat2name(Xᵒ, Dgc(Xᵒ, P), df.NAME, df.ST)
W = alloc(Xᵒ)
W[:, 1:10]

5.422553719422339e8
5.232959819069187e8
5.228999972556481e8
5.217395515582672e8
4.8770484544541514e8
4.523056879540197e8
4.286632288615675e8
4.2850005115017635e8
4.2850005115017575e8
4.2850005115017575e8
NF1: 3.0 miles from Cary, NC
NF2: 0.0 miles from Charlotte, NC

2×10 SparseMatrixCSC{Int64, Int64} with 10 stored entries:
 16432  13636  94589      ⋅      ⋅       ⋅  61960       ⋅       ⋅  283506
     ⋅      ⋅      ⋅  19092  57303  174721      ⋅  874579  105240       ⋅

---

Question 2.4.1

After adding the allocation constraint, still found Cary and Charlotte locations for the two NFs, but the TC increased from 3.9e8 to 4.3e8. Why?

(a) Adding any constraint to a minimization problem usually increases TC even though the same final NF locations were found.

(b) Since ALA was started from random initial locations, the local minimum TCs found vary.

(c) The ALA with the constraint only ran once; the best of multiple runs would have found the same result.

(d) Adding any constraint to a minimization problem usually increases TC only when different final NF locations are found.

---

Ex: Best Retail Warehouse Locations in U.S.

using Logjam.DataTools, Logjam.MapTools
using GeoMakie, DataFrames

df = filter!(r -> r.ISCUS == true && r.POP > 0, uszcta3())
select!(df, :ZCTA3, :LAT, :LON, :POP)

fig, ax, h = makemap(region=:CUS)
h[1].alpha = 0.5
scatter!(ax, df.LON, df.LAT; marker='.', markersize=12, color=:red)
fig

Random.seed!(4161)

P = hcat(df.LON, df.LAT)
w = df.POP

TC(X) = sum(w .* Dgc(X, P))
Xᵒ = optimize(TC, randX(P, 1)).minimizer

1×2 Matrix{Float64}:
 -87.4162  39.2563

df2 = filter!(r -> r.ISCUS == true && r.POP > 100_000, usplace())
lonlat2name(Xᵒ, Dgc(Xᵒ, hcat(df2.LON, df2.LAT)), df2.NAME, df2.ST)

NF1: 77.0 miles from Indianapolis (balance), IN

TC(W, X) = sum(W .* Dgc(X, P))
alloc(X) = sparse([argmin(c) for c in eachcol(Dgc(X, P))],  # Using P and X to calc m and n
    1:size(P, 1), w, size(X, 1), size(P, 1))
loc(W, X₀) = optimize(X -> TC(W, X), X₀).minimizer

function ala(Xᵒ)
    TCᵒ, done = Inf, false
    while !done
        W = alloc(Xᵒ)
        X′ = loc(W, Xᵒ)
        TC′ = TC(W, X′)
        println(TC′)
        if TC′ < TCᵒ
            TCᵒ, Xᵒ = TC′, X′
        else
            done = true
        end
    end
    return Xᵒ, TCᵒ 
end

Xᵒ = ala(randX(P, 2))[1]
lonlat2name(Xᵒ, Dgc(Xᵒ, hcat(df2.LON, df2.LAT)), df2.NAME, df2.ST)

1.584087205166322e11
1.5840781896108435e11
1.5840781896108435e11
NF1: 106.0 miles from Visalia, CA
NF2: 55.0 miles from Lexington-Fayette urban county, KY

Xᵒ = ala(randX(P, 3))[1]
lonlat2name(Xᵒ, Dgc(Xᵒ, hcat(df2.LON, df2.LAT)), df2.NAME, df2.ST)

1.6769439557681995e11
1.5051806714205005e11
1.391806394177039e11
1.2636438638319363e11
1.2441268829488632e11
1.2414860506145445e11
1.2403524390594615e11
1.2385531450674121e11
1.2372610832284201e11
1.2299806129175824e11
1.2250583409057196e11
1.223157341239782e11
1.2224607368276477e11
1.2223110787475739e11
1.2221027277468652e11
1.2219654053297597e11
1.2219601276871376e11
1.2219601276871376e11
NF1: 95.0 miles from Visalia, CA
NF2: 14.0 miles from Allentown, PA
NF3: 86.0 miles from Clarksville, TN

scatter!(ax, Xᵒ[:, 1], Xᵒ[:, 2]; marker=:dtriangle, markersize=12, color=:black)
fig

X₀ = randX(P, 9)
Xᵒ = ala(X₀)[1]
lonlat2name(Xᵒ, Dgc(Xᵒ, hcat(df2.LON, df2.LAT)), df2.NAME, df2.ST)

7.44567408497061e10
6.436790841140675e10
6.165268657216493e10
5.9510791738966606e10
5.849734291500484e10
5.8399751612027115e10
5.838690504341328e10
5.838028904374029e10
5.83760496196261e10
5.837450190632683e10
5.837431496183205e10
5.8374295417651276e10
5.837428392252685e10
5.8374255341685455e10
5.837424822464474e10
5.837424822464474e10
NF1: 6.0 miles from Denver, CO
NF2: 11.0 miles from Elizabeth, NJ
NF3: 5.0 miles from Mesa, AZ
NF4: 61.0 miles from Waco, TX
NF5: 2.0 miles from Chicago, IL
NF6: 53.0 miles from Tacoma, WA
NF7: 54.0 miles from Athens-Clarke County unified government (balance), GA
NF8: 9.0 miles from Santa Clarita, CA
NF9: 19.0 miles from Lakeland, FL

# Integrated Formulation
using SparseArrays

function TCint(X)
    D = Dgc(X, P)
    α = [argmin(c) for c in eachcol(D)]
    n, m = size(D)
    return sum(sparse(α, 1:m, w, n, m) .* D)
end

Xᵒ = optimize(TCint, X₀).minimizer
lonlat2name(Xᵒ, Dgc(Xᵒ, hcat(df2.LON, df2.LAT)), df2.NAME, df2.ST)
TCint(Xᵒ)

NF1: 53.0 miles from Salt Lake City, UT
NF2: 11.0 miles from Newark, NJ
NF3: 5.0 miles from Centennial, CO
NF4: 43.0 miles from Mesquite, TX
NF5: 10.0 miles from Chicago, IL
NF6: 18.0 miles from Tacoma, WA
NF7: 54.0 miles from Athens-Clarke County unified government (balance), GA
NF8: 14.0 miles from Pasadena, CA
NF9: 8.0 miles from Orlando, FL

5.9777689716620674e10

Xᵒ = optimize(TCint, randX(P, 9)).minimizer
lonlat2name(Xᵒ, Dgc(Xᵒ, hcat(df2.LON, df2.LAT)), df2.NAME, df2.ST)
TCint(Xᵒ)

NF1: 54.0 miles from Tacoma, WA
NF2: 11.0 miles from Pasadena, CA
NF3: 17.0 miles from Des Moines, IA
NF4: 40.0 miles from Athens-Clarke County unified government (balance), GA
NF5: 22.0 miles from Lakewood, CO
NF6: 3.0 miles from Newark, NJ
NF7: 12.0 miles from Fort Wayne, IN
NF8: 18.0 miles from Lakeland, FL
NF9: 54.0 miles from Waco, TX

5.7143602374936134e10

To add a row to the table: