2. Loc 7: Logistics Network Design¶

ISE 754, Fall 2024

Package Used: No new packages used.

1. Popco Bottling Company Example¶

Input Data¶

In [1]:
using Logjam.DataTools, DataFrames, CSV

DC = DataFrame(CSV.File("PopcoData.csv"))
Out[1]:
42×5 DataFrame
17 rows omitted
RowLATLONDEMANDPROD_COSTDIST_COST
Float64Float64Int64Int64Int64
144.06-103.187534891864808011884204
245.7731-108.5241092058664711198643
345.9744-112.513466071525932511326980
445.6685-110.976516261707841313293327
546.905-114.04841968125809034487319
639.7747-104.9735262087312886332131128
738.8827-104.81663570156959842297982
838.2341-104.61341145142328169397847
939.0685-108.588563792134649015721042
1044.0189-107.958916991982061016644210
1142.8556-106.36658141071221019698
1241.6606-109.215601381812444710351078
1346.4217-117.002567892019666214598824
⋮⋮⋮⋮⋮⋮
3137.6122-122.0833291603925881816453611
3237.625-120.9921316643788287112624983
3338.4637-122.743118467765734762986
3438.4832-121.399225518383202058576069
3539.141-121.612578831534457413271637
3644.9671-122.99757603203287001666417
3745.0766-122.79840593136969664957374
3844.0199-123.03112759233102519885323
3942.369-122.88382690264119386740091
4044.0743-121.30149374193145931391628
4147.5827-122.2992029663702304311416920
4247.233-122.47390519261239244134356
In [2]:
rename!(DC, :DEMAND => :fDC, :PROD_COST => :TPC, :DIST_COST => :TDC)
first(DC, 5)
Out[2]:
5×5 DataFrame
RowLATLONfDCTPCTDC
Float64Float64Int64Int64Int64
144.06-103.187534891864808011884204
245.7731-108.5241092058664711198643
345.9744-112.513466071525932511326980
445.6685-110.976516261707841313293327
546.905-114.04841968125809034487319
In [3]:
using GeoMakie, Logjam.MapTools

fig, ax = makemap(DC.LON, DC.LAT)

scatter!(ax, DC.LON, DC.LAT; marker='.', markersize=36, color=:red)
fig
Out[3]:
No description has been provided for this image
In [4]:
using CairoMakie
dcf() = display(current_figure())

scatter(DC.fDC, DC.TPC)
dcf();
No description has been provided for this image

Step 1: Estimate fixed cost¶

In [5]:
using Optim

ŷ(p, x) = p[1] .+ p[2]*x
loss(p, x, y) = sum((y .- ŷ(p, x)).^2)
k, cp = optimize(p -> loss(p, DC.fDC, DC.TPC), [0., 1.]).minimizer
Out[5]:
2-element Vector{Float64}:
   1.0603652532813296e7
 126.15062860391504
In [6]:
lines!([0, maximum(DC.fDC)], [k, k + cp*maximum(DC.fDC)])
dcf();
No description has been provided for this image

Step 2: Allocate customer demand to DCs¶

In [7]:
# Get aggregate demand points
z3 = uszcta3()
z3.IDX .= 0     # Add index of allocated DC (0 if not allocated)
first(z3, 5)
Out[7]:
5×8 DataFrame
RowZCTA3LATLONPOPALANDAWATERISCUSIDX
Int64Float64Float64Int64Float64Float64BoolInt64
11042.2125-72.57324728391276.1338.492true0
21142.1064-72.550417176040.9461.716true0
31242.4562-73.2286128186902.05416.827true0
41342.589-72.488282872786.42827.535true0
51442.5801-71.7824224307474.26814.954true0
In [8]:
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

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

Daa(X, Xa, a) = [max(dgc(i, j), (2/3)*sqrt(k)/pi) 
    for i in eachrow(X), (j, k) in zip(eachrow(Xa), a)]
Out[8]:
Daa (generic function with 1 method)
In [9]:
# Allocate customer demand to DCs
for r in eachrow(z3)
    d = Dgc([r.LON r.LAT], hcat(DC.LON, DC.LAT)) * 1.2   # Only need for est road distance
    idx = argmin(d)[2]
    if d[idx] <= 200                                     # Max service distance
        r.IDX = idx
    end
end
filter!(r -> r.IDX != 0, z3)                             # Remove unallocated demand points
first(z3, 5)
Out[9]:
5×8 DataFrame
RowZCTA3LATLONPOPALANDAWATERISCUSIDX
Int64Float64Float64Int64Float64Float64BoolInt64
157543.9175-100.4565213214664.4307.285true1
257645.3991-101.1152031011570.9271.706true1
357744.1001-103.29220643617529.358.493true1
459045.7035-108.9418741025896.4105.174true2
559145.7593-108.511143646813.6693.984true2

Step 3: Allocate DC demand to customers based on population¶

In [10]:
gdf = groupby(z3, :IDX)
Out[10]:

GroupedDataFrame with 42 groups based on key: IDX

First Group (4 rows): IDX = 1
RowZCTA3LATLONPOPALANDAWATERISCUSIDX
Int64Float64Float64Int64Float64Float64BoolInt64
157543.9175-100.4565213214664.4307.285true1
257645.3991-101.1152031011570.9271.706true1
357744.1001-103.29220643617529.358.493true1
469342.1465-103.2986722612281.290.06true1

⋮

Last Group (4 rows): IDX = 42
RowZCTA3LATLONPOPALANDAWATERISCUSIDX
Int64Float64Float64Int64Float64Float64BoolInt64
198347.4039-122.5568181345243.77282.374true42
298447.2011-122.475450094160.54812.763true42
398546.9764-123.0235341225358.61355.761true42
498946.6216-120.4493014725330.9755.173true42
In [11]:
length(gdf) == nrow(DC)    # Test for DCs without demand
Out[11]:
true
In [12]:
z3 = transform(gdf, :POP => sum => :DCPOP)
first(z3, 5)
Out[12]:
5×9 DataFrame
RowZCTA3LATLONPOPALANDAWATERISCUSIDXDCPOP
Int64Float64Float64Int64Float64Float64BoolInt64Int64
157543.9175-100.4565213214664.4307.285true1346104
257645.3991-101.1152031011570.9271.706true1346104
357744.1001-103.29220643617529.358.493true1346104
459045.7035-108.9418741025896.4105.174true2265435
559145.7593-108.511143646813.6693.984true2265435
In [13]:
DC.IDX .= 1:nrow(DC)
first(DC, 5)
Out[13]:
5×6 DataFrame
RowLATLONfDCTPCTDCIDX
Float64Float64Int64Int64Int64Int64
144.06-103.1875348918648080118842041
245.7731-108.52410920586647111986432
345.9744-112.5134660715259325113269803
445.6685-110.9765162617078413132933274
546.905-114.048419681258090344873195
In [14]:
z3 = leftjoin(z3, DC[!, [:IDX, :fDC]], on=:IDX)
Out[14]:
162×10 DataFrame
137 rows omitted
RowZCTA3LATLONPOPALANDAWATERISCUSIDXDCPOPfDC
Int64Float64Float64Int64Float64Float64BoolInt64Int64Int64?
157543.9175-100.4565213214664.4307.285true134610453489
257645.3991-101.1152031011570.9271.706true134610453489
357744.1001-103.29220643617529.358.493true134610453489
459045.7035-108.9418741025896.4105.174true226543510920
559145.7593-108.511143646813.6693.984true226543510920
659346.4994-105.4763437922237.8219.858true226543510920
759447.6418-111.37913327121677.6175.043true321918346607
859646.5932-111.989859126238.378.791true321918346607
959745.7681-111.6519274914989.685.468true443494051626
1059846.8957-114.12820906912823.1148.578true533786041968
1159948.2871-114.4751287918066.36252.885true533786041968
1269342.1465-103.2986722612281.290.06true134610453489
1379831.2768-105.703732949658.69.861true21128266337586
⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮
15198248.2761-122.2759840035441.9369.439true413834335202966
15298347.4039-122.5568181345243.77282.374true42210382290519
15398447.2011-122.475450094160.54812.763true42210382290519
15498546.9764-123.0235341225358.61355.761true42210382290519
15598645.793-122.5826653353831.0195.88true37281754340593
15698847.5845-119.88924195210324.0186.32true413834335202966
15798946.6216-120.4493014725330.9755.173true42210382290519
15899047.6799-117.3731686912223.1226.286true13161723756789
15999147.6773-117.70912101411778.0232.138true13161723756789
16099247.6847-117.374385114336.7641.33true13161723756789
16199346.2929-119.0444046367537.07155.382true13161723756789
16299446.3593-117.29522324664.8714.416true13161723756789
In [15]:
z3.f = map(r -> r.fDC * r.POP / r.DCPOP, eachrow(z3))
first(z3, 5)
Out[15]:
5×11 DataFrame
RowZCTA3LATLONPOPALANDAWATERISCUSIDXDCPOPfDCf
Int64Float64Float64Int64Float64Float64BoolInt64Int64Int64?Float64
157543.9175-100.4565213214664.4307.285true1346104534898056.79
257645.3991-101.1152031011570.9271.706true1346104534893138.83
357744.1001-103.29220643617529.358.493true13461045348931903.9
459045.7035-108.9418741025896.4105.174true2265435109203596.05
559145.7593-108.511143646813.6693.984true2265435109205909.6

Step 4: Nominal transport rate¶

In [16]:
dDC2r(r) = Daa([DC[r.IDX,:LON] DC[r.IDX, :LAT]],[r.LON r.LAT], r.ALAND)[1]
rnom = sum(DC.TDC) / sum(map(r -> r.f * dDC2r(r), eachrow(z3)))
Out[16]:
2.8625376133920297

Step 5: Cost matrix for UFL¶

In [17]:
# Include demand + DC locations for NF sites
NF = hcat(vcat(z3.LON, DC.LON), vcat(z3.LAT, DC.LAT))
D = Daa(NF, hcat(z3.LON, z3.LAT), z3.ALAND)
C = rnom*z3.f'.*D
Out[17]:
204×162 Matrix{Float64}:
 5.92658e5  9.64745e5  1.292e7    4.46354e6  …  3.73774e7       1.87629e6
 2.47632e6  2.05098e5  1.27387e7  3.90231e6     3.51108e7       1.74961e6
 3.26274e6  1.25328e6  2.56588e6  3.06537e6     3.17408e7       1.56637e6
 1.00004e7  3.40615e6  2.71957e7  3.51524e5     1.97778e7       9.04456e5
 9.55826e6  3.22017e6  2.55531e7  3.51524e5     2.05911e7       9.49078e5
 6.9755e6   2.00204e6  1.79764e7  1.80016e6  …  2.62702e7       1.262e6
 1.34957e7  4.59604e6  4.19432e7  1.82073e6     1.51729e7       6.56167e5
 1.36062e7  4.74373e6  4.16259e7  1.62926e6     1.36833e7       5.67583e5
 1.2975e7   4.57954e6  3.87712e7  1.34531e6     1.44956e7       6.13991e5
 1.5991e7   5.66665e6  5.10211e7  2.68552e6     9.64018e6       3.47211e5
 1.695e7    5.94058e6  5.54607e7  3.25191e6  …  1.03528e7       4.20822e5
 4.34985e6  2.24362e6  1.23272e7  3.8397e6      3.37139e7       1.68301e6
 2.1192e7   9.04313e6  8.17951e7  1.04155e7     5.11748e7       2.71901e6
 ⋮                                           ⋱  ⋮          
 2.7867e7   1.08369e7  9.82877e7  9.02859e6     2.52029e7       1.46342e6
 2.66438e7  1.03851e7  9.35689e7  8.5855e6      2.46943e7       1.41917e6
 2.80007e7  1.08332e7  9.85572e7  8.90098e6     2.32962e7       1.37301e6
 2.64742e7  1.0266e7   9.26283e7  8.32624e6  …  2.2482e7        1.30789e6
 2.62754e7  1.01495e7  9.16694e7  8.11076e6     2.07827e7       1.22181e6
 2.56158e7  9.54883e6  8.85837e7  7.03751e6     8.61394e6       6.54261e5
 2.53801e7  9.45183e6  8.76597e7  6.92664e6     8.12283e6       6.27847e5
 2.58083e7  9.67618e6  8.92877e7  7.19273e6     1.01561e7       7.23769e5
 2.61116e7  9.88787e6  9.05093e7  7.50102e6  …  1.34562e7  875048.0
 2.38326e7  8.91583e6  8.1461e7   6.3279e6      7.67121e6  562839.0
 2.48895e7  9.12562e6  8.61517e7  6.649e6       7.22077e6       5.6225e5
 2.50296e7  9.19869e6  8.66208e7  6.7083e6      7.10896e6       5.65906e5

Step 6: Solve UFL¶

In [18]:
include("ufl_heuristics.jl")

y, TC = ufl(fill(k, size(C, 1)), C)

  Add: 7.645948563275278e8
 Xchg: 7.35417248086411e8
  Add: 7.35417248086411e8
 Drop: 7.351590430439293e8
 Xchg: 7.34021255286948e8
  Add: 7.34021255286948e8
 Drop: 7.34021255286948e8
Out[18]:
([181, 17, 143, 82, 150, 163, 101, 24, 198, 128  …  108, 14, 20, 146, 54, 158, 36, 133, 8, 44], 7.34021255286948e8)
In [19]:
TCorig = k*nrow(DC) + sum(DC.TDC)
println("% reduction is TC = ", 100*(TCorig - TC)/TCorig, "\n")
DataFrame(_=["No. of NFs", "TDC","TC"],
    Orig=[nrow(DC), sum(DC.TDC), TCorig],
    New=[length(y), TC - k*length(y), TC])

% reduction is TC = 17.96078981797282
Out[19]:
3×3 DataFrame
Row_OrigNew
StringFloat64Float64
1No. of NFs42.027.0
2TDC4.49367e84.47723e8
3TC8.9472e87.34021e8
In [20]:
scatter!(ax, NF[y, 1], NF[y, 2]; marker=:circle, markersize=9, color=:transparent, 
    strokecolor=:green, strokewidth=1)
fig
Out[20]:
No description has been provided for this image

Save C matrix as CSV and re-print k¶

In [21]:
CSV.write("PopcoCmatrix.csv", DataFrame(C, :auto))
println("k = ", k)

k = 1.0603652532813296e7