Where the network’s structure once provided the most efficient solution, it now serves as the intuitive model for the more general method of linear programming.
New Julia packages used
LightOSM downloads OpenStreetMap data, the streets of an area or its buildings, and builds a network from them, with each street’s tags, such as its class, name and speed limit, kept with its links. It can also compute shortest paths on the network it builds.
NearestNeighbors finds the points nearest to a given point among many, quickly, by first organizing the points into a tree, so that a search need look at only a few of them. It works in any number of dimensions and with several measures of distance, and other packages build on it wherever many points must be searched for the nearest one.
New Logjam functions used
mapbbox: Calculates the bounding box for a set of geographic coordinates, with optional expansion along the x and y axes.
osm_roads: Download OpenStreetMap drivable roads for a bounding box and return Logjam-compatible DataFrames. Requires LightOSM and NearestNeighbors to be loaded (using LightOSM, NearestNeighbors).
When should an actual road network be used, rather than the great circle distance scaled up by a circuity factor (see Lecture 2.2, Sec. 7)? The following is a good rule of thumb: below about 50 miles, a circuity factor is probably not what should be used, because the trip is within a metro area, and metro areas can be very roundabout: rarely can a trip go straight from here to there, and it typically has to go this way to get on that road and then around. So in a metro area, great circle distance is not recommended, and the actual street map has to be used. The good news is that a metro area’s street map is small, and OpenStreetMap supplies it. All of North Carolina’s streets would be a very large download with no easy way of filtering it, but the streets of Raleigh alone fit on a computer and can be worked with, so in practice the size is not an issue. For all of North Carolina, or all of the United States, the circuity factor of Lecture 2.2, Sec. 7 remains a good shortcut.
For a few trips, a distance could be looked up by hand in Google Maps, but with hundreds of facilities that is not a solution. Google Maps has a programming interface that returns the distance between two places automatically, but where it was once free up to a generous daily limit, it now charges for its use and has grown more restrictive. The good option is OpenStreetMap, which years ago was not a real alternative but has become a good one, and which can be used from Julia, so road distances can be computed within a script rather than looked up.
Logjam’s osm_roads downloads every drivable street inside a box of longitude and latitude from OpenStreetMap’s public server, Overpass, and returns the same two tables, nodes and links, as the FAF5 networks of Lecture 4.1, so every function of Lecture 4.1, Sec. 4 applies unchanged. Because every street is returned, the number of nodes and links can grow quite large. A bounding box around a portion of Raleigh’s streets about 6 miles on a side results in around 82,000 links requiring 6 MB and can be downloaded in a few seconds. A bounding box covering all of Raleigh would require around a quarter to half a million links and about 20 to 40 MB, and take several minutes to download. The entire state of North Carolina would represent around 15 million links and over a gigabyte, but for distances across North Carolina, great circle distances and the circuity factor are a reasonable way to estimate actual road distances.
The streets of Ex. 1 are found the way they would be for any problem: the places that matter are located first, a box is put around them, and the streets inside the box are downloaded. Logjam’s loc2lonlat can be used to look up an address or other place with OpenStreetMap’s online geocoder, Nominatim. It returns the exact longitude and latitude of an address, while for aggregate demand points like ZIP codes, cities, etc., it returns a central point of the area (for a county or a state, its center of population), along with an estimate of uncertainty, which corresponds to the average distance of a point within the area representing the aggregate demand point from the centroid.
The two NC State buildings of Ex. 1 are looked up by their street addresses, with four queries that give less than an address:
# Code block 1: the two buildings of Ex. 1, and four coarser queriesbldg = ["Fitts-Woolard Hall", "D. H. Hill Library"]qry = ["915 Partners Way, Raleigh, NC 27606", # the two buildings"2 Broughton Dr, Raleigh, NC 27695","Raleigh, NC", "27695", "Raliegh, NC", "Raleigh"]loc =loc2lonlat(qry; cache_dir ="data")prt(select(loc, Not(:LON, :LAT))) # how each query was answeredx =Float64.(loc.LON[1:2]) # the two buildings, for Ex. 1, asy =Float64.(loc.LAT[1:2]) # numbers (the column holds a missing)@show x@show y;
1
Given a vector of queries, loc2lonlat returns a table with a row for each: the query (INPUT), the LON and LAT found, the tier the answer came from (GC_SOURCE), its uncertainty in miles (GC_UNCERT) and its status (GC_STATUS). The table is printed without LON and LAT, which would make it wider than the page, and the two buildings’ are shown below it. Only an answer at the ADDRESS tier is saved, in data/geocode_cache.csv, so the two addresses are read from the copy saved in data/; the other four are answered from Logjam’s Census tables, and none of the six needs the Internet.
INPUT GC_SOURCE GC_UNCERT GC_STATUS
───────────────────────────────────────────────────────────────────────
915 Partners Way, Raleigh, NC 27606 ADDRESS 0.00 OK
2 Broughton Dr, Raleigh, NC 27695 ADDRESS 0.00 OK
Raleigh, NC PLACE 5.46 OK
27695 POSTALCODE 0.23 OK
Raliegh, NC STATE 99.23 PARTIAL
Raleigh missing FAIL
x = [-78.6758457, -78.6695716]
y = [35.7706306, 35.7875615]
Logjam’s mapbbox puts a box around the two buildings, widened on each side by a multiple of its width and of its height. The buildings lie nearly north and south of each other, so the box is widened more across than along, enough to take in the faster roads around them. osm_roads downloads the streets in the box, and cropnetwork keeps those inside it. A copy of the streets, fetched on 2026-10-14 and saved in the data folder, is read through the same call that would download it:1
# Code block 2: a box around the two buildings, and the streets in it# osm_roads needs LightOSM and NearestNeighbors loaded, even when, as# here, it reads the copy of the streets already saved in data/usingLightOSM, NearestNeighborsbb =mapbbox(x, y; xexpand =2, yexpand =0.5) # (lon, lon), (lat, lat)dfN, dfL =osm_roads((bb[1]..., bb[2]...); cache_dir ="data")dfN, dfL =cropnetwork(dfN, dfL, x, y; xexpand =2, yexpand =0.5)@showround.(bb[1], digits =3)@showround.(bb[2], digits =3)@shownrow(dfN)@shownrow(dfL);
1
osm_roads takes the box as one tuple of four numbers, the least and greatest longitude and then latitude, so mapbbox’s two pairs are joined with .... It downloads whole cells of a grid 0.05° of longitude and latitude on a side, every cell the box touches, so that nearby boxes share a download; here one cell, about 2.8 by 3.5 miles, holds the box. With cache_dir, it reads the cells from that folder if they are saved there, and downloads and saves them only if not; the copy read here was saved on 2026-10-14.
2
cropnetwork keeps the nodes inside the box mapbbox gives for the same points and expansions, and the links between them, so the network is the box’s and no more.
NoteA county’s bounding box from its Census block groups
The box around a county or a state can be computed from Logjam’s Census block groups; here, Wake County’s:
# Code block 3: a county's box, from the Census block groups in itco =only(filter(r -> r.ST ==:NC && r.NAME =="Wake", uscounty()))B =filter(r -> r.STFIP == co.STFIP && r.COFIP == co.COFIP,uscenblkgrp()) # Wake's block groupsmpd =dgc((0, 0), (0, 1)) # miles in a degree of latituderad =sqrt.(B.ALAND ./π) # radius as a disk, midlon = rad ./ (mpd .*cosd.(B.LAT)) # that radius in degrees of londlat = rad ./ mpd # and of latbbW = (minimum(B.LON .- dlon), maximum(B.LON .+ dlon),minimum(B.LAT .- dlat), maximum(B.LAT .+ dlat))wide = (bbW[2] - bbW[1]) * mpd *cosd((bbW[3] + bbW[4]) /2)high = (bbW[4] - bbW[3]) * mpd@shownrow(B)@showround.(bbW, digits =3)@showround(Int, wide) # mi@showround(Int, high); # mi
1
A block group is a point, its center of population, with its land area ALAND in square miles, so the box around the points alone falls short of the county’s edge. Taking each as a disk of its land area, of radius \sqrt{\text{ALAND}/\pi}, and widening the box by that radius around each point reaches the edge. A degree of latitude is mpd miles, the great circle distance of Lecture 2.2 between two points one degree apart on a meridian, and a degree of longitude is shorter by the cosine of the latitude.
NoteThe links by functional class, and their speeds
The tables have the columns of Lecture 4.1’s, and what is in them differs. The links can be counted by the functional class osm_roads gives each, from 1, Interstates, to 7, local streets:
FCLASS Class Links Miles mph
──────────────────────────────────────────────────
3 Principal arterial 863 10.10 30, 45
4 Minor arterial 1,195 13.50 25, 35
5 Major collector 1,321 14.50 20, 30
6 Minor collector 754 8.00 25
7 Local 3,867 54.40 25
Three of the five classes in the box show two speeds because osm_roads gives a speed to each OpenStreetMap road type, not to each class, and those classes hold two types each, a road and the ramps joining it, which OpenStreetMap tags as the road’s _link type (secondary_link for a secondary), the ramps at the lower speed. The larger box of Ex. 3 holds Interstates and other freeways as well, and in class 7 living streets, at 15 mph, beside residential streets, at 25. A living street is, in OpenStreetMap’s words, a “road with very low speed limits and other pedestrian friendly traffic rules,” where pedestrians have at least the same right of way as other road users.2 Classes 6 and 7 in this box are all at 25 mph.
# Code block 5: residential streets, one-way links and link lengthusingStatisticsarcs =Set(zip(dfL.SRC, dfL.DST)) # a one-way link has no reverseoneway =count(((s, d),) -> (d, s) ∉ arcs, zip(dfL.SRC, dfL.DST))pctlocal =round(100*count(==(7), dfL.FCLASS) /nrow(dfL), digits =1)pctoneway =round(100* oneway /nrow(dfL), digits =1)medft =round(Int, 5280*median(dfL.DIST))@show pctlocal # % of the links on local streets@show pctoneway # % of the links one-way@show medft # median link length, ft@showall(==(1), dfL.DIR); # every link has DIR 1
1
DIR is the direction column of Lecture 4.1’s links: 0 for a two-way road, which can be traveled from SRC to DST and back, 1 for a one-way link traveled from SRC to DST, and -1 for one traveled the other way. osm_roads gives a two-way street one link in each direction, so every link has DIR 1, which this line confirms (true); a one-way street is then simply a link with no reverse, which is how oneway counts them above.
Table 1: FAF5 and OpenStreetMap as road networks, the second measured on Code block 2’s streets.
FAF5 (Lecture 4.1)
OpenStreetMap
Source
federal freight network, released periodically
volunteer-edited map, changed daily
Access
part of Logjam: faf5nodes() and faf5links() return it, ready to use
on OpenStreetMap’s server, not part of Logjam: Logjam’s osm_roads downloads it a box at a time
Roads
freight roads, arterials and up
every street, 48.3% of the links local
Links
road segment, DIR 0 if two-way
one per direction of travel, DIR 1 on every link, 23.8% one-way
Link length
fraction of a mile
median 39 ft, every bend a node
Speed
posted speed of each link
one speed for each road type
License
public domain
ODbL; copyright OpenStreetMap contributors, credited on each map
Each difference in Table 1 shows in the examples that follow: a trip and its return can take different streets, a trip can begin and end at an address in a neighborhood, and a travel time is an estimate from speeds by class, which the Landmark check of Ex. 1 measures.
Example 1: Fitts-Woolard to D. H. Hill Library
Determine the shortest route by distance and the fastest route by time from Fitts-Woolard Hall, on NC State’s Centennial Campus, to D. H. Hill Library, on its main campus, and back, on the OpenStreetMap streets of Raleigh, and the streets each route follows.
The two graphs, on miles and on minutes, are built as Lecture 4.1, Code block 53 builds them, a minute at each link’s speed and a connector at 20 mph:
# Code block 6: the streets around the two buildings, as two graphsusingGraphs, SimpleWeightedGraphsNc, Lc =addconnectors(dfN, dfL, x, y; add_nf_nf =false)con =isequal.(Lc.SOURCE, "CONNECTOR")tmin = [c ? 60d /20:60d / v # minutes, a connector for (c, d, v) inzip(con, Lc.DIST, Lc.SPEED)] # 20 mphfw, bw = Lc.DIR .!=-1, Lc.DIR .!=1I, J = [Lc.SRC[fw]; Lc.DST[bw]], [Lc.DST[fw]; Lc.SRC[bw]]gD =SimpleWeightedDiGraph(I, J, [Lc.DIST[fw]; Lc.DIST[bw]]; combine = min) # milesgT =SimpleWeightedDiGraph(I, J, [tmin[fw]; tmin[bw]]; combine = min) # minutesD, PD =shortestpaths(gD, 2) # the buildings are nodes 1 and 2T, PT =shortestpaths(gT, 2)@shownrow(Nc)@shownrow(Lc);
1
addconnectors joins the two buildings to the streets of Code block 2’s box: each becomes a node, the buildings nodes 1 and 2, linked to the street nodes around it by connectors, links whose SOURCE is CONNECTOR.
2
c ? a : b is a when c is true and b otherwise. A connector has no speed of its own (addconnectors leaves its SPEED missing), so it is timed at 20 mph, and a street at its class’s speed; 60d / v is the minutes to go d miles at v mph.
3
As in Lecture 4.1, Code block 45, each link gives an arc in each direction it can be traveled. On OpenStreetMap every street link has DIR 1 (Code block 5), so fw is every link and bw only the connectors, whose DIR is 0 because they run both ways; of two arcs between the same nodes, the quicker is kept (combine = min).
nrow(Nc) = 4728
nrow(Lc) = 8006
# Code block 7: each route's miles and minutes, out and backWd, Wt =weights(gD), weights(gT)R =DataFrame(Trip =String[], Route =String[], mi =Float64[], min =Float64[])paths =Vector{Int}[]for (i, j, trip) in ((1, 2, "Out"), (2, 1, "Back")), (route, P) in (("Shortest", PD), ("Fastest", PT)) p =tracepath(P[i], i, j)push!(paths, p)push!(R, (trip, route, sum(segcost(p, Wd)), sum(segcost(p, Wt))))endR.mi, R.min =round.(R.mi, digits =3), round.(R.min, digits =2)prt(R)
Trip Route mi min
────────────────────────────
Out Shortest 1.91 4.08
Out Fastest 2.12 3.87
Back Shortest 2.00 3.90
Back Fastest 2.01 3.89
# Code block 8: the streets each route follows# The streets along path p, numbered in order, each with its first and# last step: step t is the link from p[t] to p[t+1], traveled either# way; a link with no name, such as a connector, is passed overfunctionstreets(p, L) k =Dict(zip(zip(L.SRC, L.DST), 1:nrow(L))) # link of (SRC, DST)lk(u, v) =get(k, (u, v), get(k, (v, u), 0)) # either direction nm = [L.NAME[lk(u, v)] for (u, v) inzip(p, p[2:end])] t =findall(!ismissing, nm) # the steps with a name a = [i for i ineachindex(t) if i ==1|| nm[t[i]] != nm[t[i-1]]] b = [a[2:end] .-1; length(t)] # each run's first, lastreturnDataFrame(Seg =eachindex(a), Street = nm[t[a]], Steps =string.(t[a], "-", t[b]))endfor (r, p) inzip(eachrow(R), paths)println(r.Trip, ", ", lowercase(r.Route), ":")prt(streets(p, Lc))end
Out, shortest:
Seg Street Steps
───────────────────────────────────────
1 Partners Way 2-20
2 Main Campus Drive 21-30
3 Oval Drive 31-37
4 Bilyeu Street 38-52
5 Pullen Road 66-127
6 Katharine Stinson Drive 141-181
7 West Broughton Drive 182-184
Out, fastest:
Seg Street Steps
───────────────────────────────────
1 Partners Way 2-7
2 Varsity Drive 8-26
3 Avent Ferry Road 27-56
4 Western Boulevard 57-68
5 Dan Allen Drive 69-117
6 Hillsborough Street 118-138
7 Gardner Street 139-141
8 Founders Drive 142-145
Back, shortest:
Seg Street Steps
───────────────────────────────────
1 Founders Drive 2-7
2 Lampe Drive 8-10
3 Hillsborough Street 11-35
4 Pullen Road 44-134
5 Bilyeu Street 137-151
6 Oval Drive 152-159
7 Main Campus Drive 160-168
8 Partners Way 169-187
Back, fastest:
Seg Street Steps
───────────────────────────────────
1 Founders Drive 2-7
2 Lampe Drive 8-10
3 Hillsborough Street 11-41
4 Pullen Road 43-135
5 Bilyeu Street 138-152
6 Oval Drive 153-160
7 Main Campus Drive 161-169
8 Partners Way 170-188
Out (Fig. 1), the shortest route, 1.906 mi, keeps to the campus streets, which OpenStreetMap classes with residential and local streets at 25 mph; the fastest, 2.125 mi but 3.87 min against 4.08, goes around by Varsity Drive, Western Boulevard, Dan Allen Drive and Hillsborough Street, faster roads for a longer way. Back, the routes differ from the routes out, because some of the streets are one-way, and here the shortest and the fastest take the same streets.
Show the code that draws this map
# Code block 9: the two routes out, on the streets of the boxusingGeoMakie# makemap draws with itsupc = colorant"#d21f26"# shortest routedcc = colorant"#1e8a3c"# fastest routebluec = colorant"#1f6fb5"# the two buildingsink = colorant"#252525"# The streets of a set of links as a graph x2ln can draw (Lecture 4.1)street(k) =SimpleWeightedGraph(Lc.SRC[k], Lc.DST[k], Lc.DIST[k]; combine = min)lo =street(.!con .& (Lc.FCLASS .>=6)) # local and minor collectorhi =street(.!con .& (Lc.FCLASS .<=5)) # the other classespS, pF = paths[1], paths[2] # out: shortest, fastestfig, ax =makemap(bb[1], bb[2]; doRoadbkgd =false, # the whole box xexpand =0, yexpand =0)resize!(fig, 540, 690)# Both street layers faint: a muted base under the routes and nameslines!(ax, x2ln(lo, Nc.LON), x2ln(lo, Nc.LAT); color = (ink, 0.12), linewidth =0.6)lines!(ax, x2ln(hi, Nc.LON), x2ln(hi, Nc.LAT); color = (ink, 0.3), linewidth =1.2)lS =lines!(ax, Nc.LON[pS], Nc.LAT[pS]; color = (supc, 0.85), linewidth =3.5)lF =lines!(ax, Nc.LON[pF], Nc.LAT[pF]; color = (dcc, 0.85), linewidth =3.5)scatter!(ax, x, y; color = bluec, markersize =14, strokecolor =:white, strokewidth =1)# The library's name on two lines, clear of the fastest routetext!(ax, x, y; text =replace.(bldg, "Hill "=>"Hill\n"), fontsize =14, color = bluec, font =:bold, align = [(:left, :top), (:right, :top)], offset = [(8, -6), (-8, -6)])# Name street nm along its direction, at fraction f of its extent in# the box xl, yl, gap px off it; side -1 puts the name below the streetusingLinearAlgebrafunctionnamestreet!(ax, N, L, nm, xl, yl; f =0.5, side =1, gap =5) k =findall(r ->isequal(r.NAME, nm) && xl[1] < N.LON[r.SRC] < xl[2] && yl[1] < N.LAT[r.SRC] < yl[2], eachrow(L))isempty(k) &&returnnothing lon, lat = N.LON[L.SRC[k]], N.LAT[L.SRC[k]] u = [lon .*cosd(mean(lat)) lat] # a degree of longitude is short u = u .-mean(u, dims =1) v =eigen(Symmetric(u'* u)).vectors[:, 2] # the street's direction v[1] <0&& (v =-v) i =sortperm(u * v)[clamp(round(Int, f *length(k)), 1, length(k))] θ =atan(v[2], v[1])text!(ax, lon[i], lat[i]; text = nm, fontsize =12, color = ink, rotation = θ, align = (:center, side >0 ? :bottom ::top), offset = gap * side .* (-sin(θ), cos(θ)))returnnothingendxl, yl = bb # the streets of the whole boxfor (nm, f, side) in (("Hillsborough Street", 0.35, 1), ("Western Boulevard", 0.49, -1), ("Avent Ferry Road", 0.12, 1), ("Pullen Road", 0.4, -1), ("Varsity Drive", 0.5, -1))namestreet!(ax, Nc, Lc, nm, xl, yl; f, side)endaxislegend(ax, [lS, lF], ["Shortest", "Fastest"]; position =:rb, labelsize =13)fig
The Open Source Routing Machine (OSRM), a router that runs on the same OpenStreetMap data, gives a landmark to check the routes against. It answers in meters and seconds, which the call converts at 1,609.344 m to the mile and 60 s to the minute before any number is compared, the Units check of Lecture 1.2. The call is shown and not run, as in Lecture 4.1, Code block 55, and its answer of 2026-10-12 is Table 2:3
# Code block 10: a landmark, OSRM's routes on the same OpenStreetMap datausingHTTP, JSON3a, b = loc[1, :], loc[2, :] # the two buildingsfor (s, t) in ((a, b), (b, a)) # out, then back url ="https://router.project-osrm.org/route/v1/driving/"*"$(s.LON),$(s.LAT);$(t.LON),$(t.LAT)"*"?alternatives=true&overview=false" r = JSON3.read(HTTP.get(url).body) # the router's answerprintln(DataFrame(DIST = [rt.distance for rt in r.routes] ./1609.344, TIME = [rt.duration for rt in r.routes] ./60))end
Table 2: OSRM’s routes between Fitts-Woolard Hall and D. H. Hill Library, fastest first, as the request of Code block 10 returned them on 2026-10-12.
Trip
Route
Distance (mi)
Time (min)
Out
1
2.25
7.5
Out
2
3.27
8.2
Back
1
1.90
6.1
The router’s routes differ from those of Code block 7 in three ways, each from a difference in how the two use the same OpenStreetMap data. Its second route out goes around by Gorman Street, 1.14 mi farther than the fastest route of Code block 7, because OpenStreetMap tags the stretch of Dan Allen Drive between Western Boulevard and Hillsborough Street closed to vehicles and open to buses: the gates NC State put at its railroad crossing restrict weekday traffic there from 9 a.m. to 5 p.m. to buses and to emergency and university service vehicles.4 The router’s car profile reads the tag and osm_roads does not, so the fastest route of Code block 7 is one a Wolfline bus can take at any hour, and a car only outside those hours. The other two, in where its fastest route out ends and in how long its routes take, come from how the router ends a route and the speed it assumes on a street, set out in the note below. Each difference traces to something the router reads and osm_roads or Code block 6 leaves out, a closed road, the direction of a one-way street, a posted speed, and finding them is what the Landmark check of Lecture 1.2 is for.
NoteWhere the router’s route ends, and the speeds it assumes
The router’s fastest route out takes the campus streets of the shortest route of Code block 7 as far as Katharine Stinson Drive and then circles by Lampe Drive, Hillsborough Street, Pogue Street, Vanderbilt Avenue and Gardner Street back to Founders Drive, 0.34 mi more, because the router ends a route at the nearest point of a drivable street, for the library a point on one-way Founders Drive 31.5 m away, which it must reach in the street’s direction; Code block 6 instead joins each building to the nodes around it by connectors that run both ways, so its route can end on West Broughton Drive. Back, the router’s route takes the streets of the route back.
The router’s times run longer, because its car profile drives a street at 80% of its posted speed limit, 20 mph on the campus streets OpenStreetMap posts at 25 mph, and adds up to 7.5 s for a turn, by the turn’s angle,5 where osm_roads gives each road type one speed, 25 mph for a residential or unclassified street, and no time for a turn: at 20 mph, the 2.25 mi of its fastest route out alone take 6.8 min of its 7.5. The difference in times is also what Table 4 in Lecture 4.1 found.
Out, the shortest route is 1.906 mi on campus streets and the fastest 3.87 min around by Western Boulevard and through the Dan Allen Drive gates, which on weekdays from 9 a.m. to 5 p.m. only buses and service vehicles may pass; back, the shortest route is 2.0 mi on other streets
2. Linear assignment problem
The linear assignment problem (LAP) is a special case of the transportation problem. Special procedures more efficient than linear programming (LP) were developed to solve the LAP and transportation problems, e.g., the Hungarian method. Now it is usually easier to transform them into an LP since solvers are so good, and special, very efficient procedures are only used for the shortest path problem (Dijkstra’s procedure, Lecture 4.1).
Many problems can be modeled as networks, where the nodes of the network represent entities and the arcs represent either flows or relationships, for example. All of the network models of Secs. 2 to 4 can be solved as linear programs. In the past, modeling a problem as a network allowed network-specific procedures to be used that were more efficient than solving a general linear programming formulation of the problem, but now even large-scale linear programs can be solved efficiently, and the use of network-specific procedures is not necessary for many problems.
NoteWhy these are called network problems
The assignment and transportation problems belong to Networks even though little in them concerns the network itself. The reason is that nowadays they are solved by formulating them as linear programs: they are really math programming problems, and calling them network problems is more a matter of name. The name came about because both can be looked at as networks, and when linear programs were slow to solve on a computer, special-purpose procedures that used the network structure, the Hungarian method the best known of them, could solve fairly large problems even on slow computers, much larger than could be solved as linear programs. Since the 2000s, linear programming solvers have been good enough that the network-specific procedures are no longer needed. The one problem that still calls for a specialized procedure is the shortest path problem of Lecture 4.1, Dijkstra’s procedure with its specialized data structures, which on a large road network is far faster than a linear program.
The linear assignment problem can be modeled as a complete bipartite network, which is a network consisting of two sets of nodes with all arcs from one set to the other (also called the matching problem on a weighted bipartite graph), as in Fig. 2, repeated from Lecture 4.1. The first set of “supply” nodes represent an entity that can be assigned to one of the second set of “demand” nodes. Each arc represents the cost of making the assignment, and the problem is solved by determining the least-cost assignment. That is, the problem is to find the least cost set of arcs that connect each destination node with a source node.
Show the code that draws this figure
# Code block 11: a network of two groups, suppliers 1-3 and DCs 4-5supc = colorant"#d21f26"# suppliersdcc = colorant"#1e8a3c"# DCsink = colorant"#252525"r =0.24# node radius, data unitsXY = [0.02.0; 0.01.0; 0.00.0; 2.41.6; 2.40.4] # node positions# One arc from node i to node j, outline to outline: kind :one (an arrow),# :two (a two-way arc, no arrowhead) or :absent (dashed); off shifts it# sideways, for a pair of parallel arcsfunctionarc!(ax, i, j, lab; kind =:one, off =0.0, lpos =0.5, lside =1) p, q = XY[i, :], XY[j, :] u = (q .- p) ./hypot((q .- p)...) nrm = [-u[2], u[1]] # unit normal s =asin(clamp(off / r, -1, 1)) # where a shifted arc meets t = p .+ r *cos(s) .* u .+ off .* nrm # the circle h = q .- (r *cos(s) +0.015) .* u .+ off .* nrmif kind ==:absentlines!(ax, [t[1], h[1]], [t[2], h[2]]; color = ink, linewidth =1, linestyle =:dash)elseif kind ==:twolines!(ax, [t[1], h[1] +0.015u[1]], [t[2], h[2] +0.015u[2]]; color = ink, linewidth =1.5)elsearrows2d!(ax, [t[1]], [t[2]], [h[1] - t[1]], [h[2] - t[2]]; color = ink, shaftwidth =1.5, tipwidth =9, tiplength =9)endif !isempty(lab) m = t .+ lpos .* (h .- t)text!(ax, m[1], m[2]; text = lab, color = ink, fontsize =14, align = (:center, :center), offset = (12lside * nrm[1], 12lside * nrm[2]))endreturnnothingendfunctiontwogroup(arcs; w =340, h =300) fig =Figure(size = (w, h)) ax =Axis(fig[1, 1]; aspect =DataAspect())hidedecorations!(ax); hidespines!(ax)for (i, j, lab, kw) in arcsarc!(ax, i, j, lab; kw...)endfor k in1:size(XY, 1) # the nodes, over the arcs x, y = XY[k, :] c = k <=3 ? supc : dccpoly!(ax, Circle(Point2f(x, y), r); color =:white, strokecolor = c, strokewidth =1.8)text!(ax, x, y; text =string(k), color = c, fontsize =15, align = (:center, :center))endxlims!(ax, -0.4, 2.8); ylims!(ax, -0.4, 2.4)return figendtwogroup([(1, 4, "6", (lpos =0.25,)), (1, 5, "10", (lpos =0.2, lside =-1)), (2, 4, "0", (lpos =0.3,)), (2, 5, "13", (lpos =0.3, lside =-1)), (3, 4, "9", (lpos =0.2,)), (3, 5, "16", (lpos =0.3, lside =-1))])
Figure 2: A complete bipartite digraph: three suppliers (red) each ship to both of two DCs (green), one arc per pair, each labeled with its weight.
Although the Hungarian method is a specialized procedure for the LAP, modeling the LAP as a network flow problem and then solving it as an LP is also effective and easier to implement. The Quadratic Assignment Problem (QAP) is when the cost of an assignment depends on the other assignments; as a result, it is more difficult to solve than the LAP.
The linear assignment problem is a basic, simple problem, also called the matching problem in graph theory. Its network is bipartite because it has two sets of nodes, the sources and the destinations, and the problem is to choose the set of arcs, each with a weight that could be a cost or a distance, that together have the minimum total weight while reaching every destination. In the network of Fig. 2, three sources and two destinations, only two arcs are needed. The arc of weight 0 from 2 to 4 is certainly used, and it covers node 4; the least-cost arc to node 5 is then the one from 1, for a total cost of 10, with source 3 left unused. With only two sources and three destinations the problem would not be feasible, since each destination needs a source of its own: there can be more sources than destinations, but not fewer. This is the typical assignment problem, and it is easily solved as a linear program. The quadratic assignment problem, used a lot in facility layout and circuit board layout, differs in that the cost of an assignment depends on what other assignments are made: the cost of assigning a source to node 5 depends on what is assigned to node 4. That allows more interesting problems to be posed, but they are much harder to solve; no procedure is known that solves the quadratic assignment problem to optimality in reasonable time as it grows, and approximations are used for it. The linear assignment problem is the easiest of them all.
The linear assignment problem is stated as a model in Model 1. It is the first of three models, each generalizing the one before it: the transportation problem of Sec. 3 generalizes it, and the minimum cost network flow problem of Sec. 4 generalizes the transportation problem. In these models S and D are the sets of supply and demand nodes, and s_i and d_j a node’s supply and demand, not the labeled set, origin and distance labels that the same letters name in Dijkstra’s procedure of Lecture 4.1.
minimize: total cost of assigning supply nodes to demand nodes
solve for:
(a) supply node assigned to each demand node.
subject to:
(a) supply: each supply node provides at most one unit;
(b) demand: each demand node requires exactly one unit.
assumptions:
(a) cost of assigning each supply node to each demand node is given;
(b) cost of an assignment does not depend on the other assignments;
(c) there are at least as many supply nodes as demand nodes.
=\begin{cases} 1, & \text{if supply node } i \text{ is assigned to demand node } j \\ 0, & \text{otherwise} \end{cases}
c_{ij}
=cost of assigning supply node i to demand node j (cost)
S
=supply nodes
D
=demand nodes.
Model 1 states two constraints, and both are rows. Row (a), supply, lets each supply node be assigned to at most one demand node, so a supply node can be left unassigned; row (b), demand, assigns each demand node exactly one supply node. The domain row asks only that each x_{ij} be nonnegative, not that it be 0 or 1, and the linear program’s solution nevertheless assigns each demand node a single supply node, as Ex. 2 shows.
Model 1 implementation: Linear assignment
# Model: linear assignmentusingJuMP, HiGHSfunctionLAP(c) S, D =1:size(c, 1), 1:size(c, 2)iflength(S) <length(D) # assumption (c)error("fewer supply nodes than demand nodes")end model =Model(HiGHS.Optimizer)@variable(model, 0≤ x[S, D])@objective(model, Min, sum(c[i, j] * x[i, j] for i ∈ S, j ∈ D))@constraint(model, supply[i ∈ S], sum(x[i, j] for j ∈ D) ≤1) # (a)@constraint(model, demand[j ∈ D], sum(x[i, j] for i ∈ S) ==1) # (b)set_silent(model)optimize!(model)iftermination_status(model) == JuMP.OPTIMALreturn (TC =objective_value(model), x =snapvals(value.(x)))enderror(termination_status(model))end
1
error stops the function there, with its message, and nothing after it runs: a call of LAP on a matrix with fewer rows (supply nodes) than columns (demand nodes) ends at this line, as in Ex. 2(c), before any model is built.
2
Only an optimal solution is returned; any other status stops with an error naming it, such as INFEASIBLE. snapvals, as in Lecture 2.7, sets a value close to an integer to that integer, so an assignment of 0.9999999997 reads as 1.
LAP (generic function with 1 method)
The implementation checks assumption (c) of Model 1 before it builds the linear program, and stops with an error if there are fewer supply nodes than demand nodes. The check is not in the formulation, Eq. 1, but it is in the concept model, and it has to be in the implementation: it is an example of what a formulation typically leaves out but an implementation needs. This is error checking: an assumption tested at the start of a function, and an error raised when it fails, as Lecture 2.2’s d2 raises one for two points of different lengths. Without it, the case would reach the solver and come back only as infeasible, which does not say which assumption was broken.
An example of the linear assignment problem would be its use by a ride-hailing service to assign available cars to passengers wanting a ride. One way of doing the assignment is to make the assignment based on which car is closest to any of the passengers and then continue until all the assignments are complete. Unfortunately, this simple greedy procedure does not, in general, provide an optimal assignment; instead, an optimal assignment can be found with the problem modeled as a network flow problem, where one unit of supply provided by a car can be used to satisfy one unit of demand required for each passenger, where the arc connecting each car to each passenger represents the cost (distance) of the car from the passenger.
Example 2: Unter’s passenger-to-car assignment
Unter, a ride-hailing service, has three cars available and three passengers wanting a ride, the distance of each car from each passenger given by the distance matrix
D = \begin{bmatrix} 4 & 8 & 5 \\ 2 & 1 & 9 \\ 6 & 3 & 7 \end{bmatrix},
one row for each car and one column for each passenger. Determine the assignment of cars to passengers, first by the greedy procedure and then as Model 1, and compare their total distances.
Unter is a competitor to Uber. A ride-hailing service has a set of cars that are available and a set of passenger requests, so it is essentially solving a linear assignment problem in determining which driver to assign to which passenger. Here there are three cars and three passengers, and the cost of an assignment is the distance of the car from the passenger, so the solution wanted is the one that minimizes the total distance with every passenger assigned a car. In terms of Model 1, the cars are the supply nodes, each providing one unit, the passengers are the demand nodes, each requiring one, and c_{ij} is d_{ij}, the distance of car i from passenger j, so the cost matrix is D.
Fig. 3 draws the problem as a complete bipartite network, as Fig. 2 does for its suppliers and distribution centers: the cars are its supply nodes, the passengers its demand nodes, and each arc carries the car’s distance from the passenger.
Show the code that draws this figure
# Code block 12: Unter's cars and passengers as a bipartite networklet D = [485; 219; 637], r =0.2# distances, node radius supc, demc = colorant"#d21f26", colorant"#1e8a3c" ink = colorant"#252525" XYc = [(0.0, 2.0), (0.0, 1.0), (0.0, 0.0)] # cars XYp = [(2.4, 2.0), (2.4, 1.0), (2.4, 0.0)] # passengers# Each distance near its car, on the side of its arc inside that# car's fan (1 above it, -1 below), so it reads with its own arc sd = [111; 11-1; -1-1-1] fig =Figure(size = (400, 330)) ax =Axis(fig[1, 1]; aspect =DataAspect())hidedecorations!(ax); hidespines!(ax)for (i, p) inenumerate(XYc), (j, q) inenumerate(XYp) u = (q .- p) ./hypot((q .- p)...) t, h = p .+ r .* u, q .- (r +0.015) .* uarrows2d!(ax, [t[1]], [t[2]], [h[1] - t[1]], [h[2] - t[2]]; color = ink, shaftwidth =1.2, tipwidth =8, tiplength =8) m = t .+0.18.* (h .- t) .+0.1* sd[i, j] .* (-u[2], u[1])text!(ax, m[1], m[2]; text =string(D[i, j]), color = ink, fontsize =13, align = (:center, :center))endfor (k, (x, y)) inenumerate([XYc; XYp]) # the nodes, over the arcs c = k <=3 ? supc : demcpoly!(ax, Circle(Point2f(x, y), r); color =:white, strokecolor = c, strokewidth =1.8)text!(ax, x, y; text =string((k -1) %3+1), color = c, fontsize =15, align = (:center, :center))endtext!(ax, [0.0, 2.4], [2.42, 2.42]; text = ["cars", "passengers"], color = [supc, demc], fontsize =14, align = (:center, :bottom))xlims!(ax, -0.5, 2.9); ylims!(ax, -0.35, 2.7) figend
Figure 3: Unter’s problem of Ex. 2 as a complete bipartite network: each of the three cars (red), the supply nodes, is joined to each of the three passengers (green), the demand nodes, and each arc is labeled with the car’s distance from the passenger, d_{ij}.
The greedy procedure comes first, in Julia, even though the linear program is the better approach. Developing a greedy procedure is not something this course asks for, but in research, small problems often come up that need some answer, and a greedy procedure can be the quickest way to get one.
Example 2(a): Greedy step by step
Assign the cars to the passengers by repeatedly making the assignment of least distance among those still possible.
argmin returns a Cartesian index, the row and the column together rather than a single number, here (2, 2): idx[1] is the row, idx[2] the column, and D[idx], in the next block, the value there.
CartesianIndex(2, 2)
D[idx] # Code block 15: its distance
1
The greedy procedure assigns the closest car and passenger first, then the closest of those left, and so on. There might be a better way of assigning the cars than this, but when an answer is needed, a greedy procedure at least gives a feasible solution, if not an optimal one, which makes it handy to have. The first step is to find the minimum value in the matrix, here in the second row and second column, where the smallest value is 1, so car 2 is assigned to passenger 2. Once an assignment is made, the car’s supply is used up, since in the linear assignment problem every supply and every demand is one; the transportation problem of Sec. 3 is different, because a supply can be more than one and the amount used is subtracted from it. So here, as soon as an assignment is made, the car’s row and the passenger’s column can be eliminated. One way to do that would be to delete the row and the column from D.
The 3 × 3 matrix would then become 2 × 2. That would work, but manipulating an array that way is usually slower and more complicated: every assignment would change the array’s size, and any programming language then has to keep rearranging it in memory, which can be really slow. In many cases it is a lot better to keep the matrix its size and make the assigned row and column impossible to choose again, which Code block 16 does by giving them a value larger than any in the matrix.
# Code block 16: assign large value so that assignment not made againD[idx[1],:] .=maximum(D)+1D[:,idx[2]] .=maximum(D)+1D
1
.= gives every element of row idx[1] the maximum, 9, plus one: 10.
2
The column gets 11, not 10, because the line before has already changed D, and its maximum is now 10: every step changes D itself.
3×3 Matrix{Int64}:
4 11 5
10 11 10
6 11 7
The value is the maximum in the matrix plus one, first for the row and then for the column. It does what is wanted, since a row or column at the maximum plus one is never picked again: there is always a smaller value left in the matrix.
idx =argmin(D) # Code block 17: find next closest passenger
CartesianIndex(1, 1)
D[idx] # Code block 18: its distance
4
# Code block 19: assign large values againD[idx[1],:] .=maximum(D)+1D[:,idx[2]] .=maximum(D)+1D
3×3 Matrix{Int64}:
13 12 12
13 11 10
13 11 7
idx =argmin(D) # Code block 20: the last passenger
CartesianIndex(3, 3)
D[idx] # Code block 21: its distance
7
TC =1+4+7# Code block 22: final cost of assignment
12
Applied again, the procedure finds the minimum at row 1, column 1, assigning car 1 to passenger 1, and a third time at row 3, column 3, assigning car 3 to passenger 3. The final cost is the sum of the original cost values of the three assignments, 1, 4 and 7, a total cost of 12.
Example 2(b): Greedy as a function
Make the greedy procedure a function, and use it to see what happens when the numbers of cars and passengers differ.
This procedure can be made into a function to allow it to be easily reused:
# Code block 23: the greedy procedure as a function# Need to re-define matrix since it was changed aboveD = [485; 219; 637]functiongreedyassignment(C) C =copy(C) # Copy-by-value instead of default copy-by-reference X =zeros(size(C))# Large value so that assignment not made again BigM =maximum(C) +1for i =1:size(C,2) idx =argmin(C) X[idx] =1 C[idx[1],:] .= BigM C[:,idx[2]] .= BigMendreturn XendX =greedyassignment(D)
1
Computed once, before the loop, where Code blocks 16 and 19 took the maximum again at each step; Big M is the name math programming gives a value larger than any in the problem.
Before the function runs, D has to be defined again, since the steps of Ex. 2(a) changed it. The function does the same thing those steps did, with one difference that matters. In MATLAB, a matrix passed into a function is copied automatically, so whatever the function does to it affects nothing outside. Julia, like Python, passes by reference instead, which is one of the reasons it is fast: what is passed is not the matrix but a pointer to where it is in memory, so a change made to it inside the function changes it outside as well, and any later code that uses D would then be working with a different matrix without knowing it. Forgetting to make a real copy is an easy error to make and a slow one to find. So the function’s first step is copy, after which it can change its own C freely without affecting the original. The rest goes through the columns making the same row and column assignments as before. The function returns X, which starts as all zeros and ends with a one wherever an assignment was made, an indicator matrix, so the total cost is the sum of the elements of D times X, again 12.
What if the number of cars does not equal the number of passengers?
# Code block 25: not feasible to have one car serve multiple passengersgreedyassignment(D[1:2,:])
2×3 Matrix{Float64}:
1.0 0.0 0.0
0.0 1.0 0.0
# Code block 26: feasible to have extra car not assignedgreedyassignment(D[:,1:2])
3×2 Matrix{Float64}:
1.0 0.0
0.0 1.0
0.0 0.0
With fewer cars than passengers, the first two rows of D only, the greedy function still returns a solution, but it is not feasible: the third passenger has not been assigned a car. The code runs without complaint and says nothing, so knowing that the answer is infeasible is left to whoever runs it. With fewer passengers than cars, the first two columns only, each passenger is still assigned a car and the third car is simply not used, which is feasible. Solved as a linear program, as in Ex. 2(c), the first case is reported as an error, because LAP checks assumption (c) of Model 1, that there are at least as many cars as passengers, before it solves anything. A greedy procedure that is going to be used could check the same thing itself, that supply equals or exceeds demand, and stop with an error message when it does not. On one small instance that seems a minor point, but code is often reused months later, by which time the limitation is forgotten, and an answer that is wrong without generating an error is the most dangerous kind: the problem can surface several steps further into the work, far from its cause.
Example 2(c): As a linear program
Solve the assignment as Model 1, and solve it again with fewer cars than passengers, and with fewer passengers than cars.
# Code block 27: Unter's assignment as a linear programTCᵒ, Xᵒ =LAP(D)@show TCᵒ; Int.(Xᵒ)
TCᵒ = 10.0
3×3 Matrix{Int64}:
0 0 1
1 0 0
0 1 0
# Code block 28: not feasible, one car shorttryLAP(D[1:2,:])catcheshowerror(stdout, e)end
1
try runs LAP, and when it stops with an error, catch prints the error’s message and the script runs on; without it, the companion script would stop at this line.
fewer supply nodes than demand nodes
# Code block 29: feasible, extra car not assignedTC₂, X₂ =LAP(D[:,1:2])@show TC₂; Int.(X₂)
TC₂ = 5.0
3×2 Matrix{Int64}:
0 0
1 0
0 1
The linear assignment problem is a really nice, simple formulation, the costs c_{ij} with the indicator variables x_{ij}: the sum over each car is at most one, each car assigned to up to one passenger, and the sum over each passenger is exactly one, each passenger assigned a car. It translates pretty much directly into a math program, Code block 27, which solves it to optimality, and the optimal total is 10. The greedy procedure’s 12 is not optimal: even on so simple a 3 × 3 instance, greedy is 20% above the optimum, which is pretty bad performance.
Before either procedure ran, the matrix alone set an expectation: each passenger’s nearest car is at distance 2, 1 and 5, so no assignment can total less than 8, and the cheapest entries, 1 and 2, sit in the same row, so not both can be used. An answer a little above 8 was the estimate to beat, which greedy’s 12 misses and the linear program’s 10 meets. Committing to an estimate before computing is the Prior check of Lecture 1.2.
Greedy assigns car 2 to passenger 2, car 1 to passenger 1 and car 3 to passenger 3, for a total distance of 12; the linear program assigns car 1 to passenger 3, car 2 to passenger 1 and car 3 to passenger 2, for a total of 10
Example 3: Unter on Raleigh’s streets
Twelve drivers and ten passengers wanting a ride are at Census block groups within 3 miles of downtown Raleigh. Assign the drivers to the passengers, by the greedy procedure and as Model 1, so that the total driving time to the pickups is least, and determine, by Little’s Law, how many drivers are on the way to a pickup at any time when 20 requests arrive an hour.
Code block 30 below generates random locations for drivers and passengers by drawing without replacement a random set of the Census block groups within 3 miles of the State Capitol, downtown, and using the population centroid of each group as its location. The seed and the two counts are set in the code, so a different seed, or more drivers or passengers, gives another instance:
# Code block 30: synthetic driver and passenger creationusingRandomnd, np =12, 10# drivers, passengerscap =loc2lonlat(["1 E Edenton St, Raleigh, NC 27601"]; cache_dir ="data")c0 = (Float64(cap.LON[1]), Float64(cap.LAT[1])) # the State Capitolbg =filter(r ->dgc(c0, (r.LON, r.LAT)) <=3, # block groups withinuscenblkgrp()) # 3 mi of itsort!(bg, [:TRFIP, :BGFIP]) # a fixed order for the drawRandom.seed!(754) # change it for another drawk =randperm(nrow(bg))[1:nd+np] # distinct block groupsU =DataFrame(ROLE = [fill("driver", nd); fill("passenger", np)], ID = [string.("D", 1:nd); string.("P", 1:np)], LON = bg.LON[k], LAT = bg.LAT[k])prt(U)
1
randperm(n) is the numbers 1 to n in a random order, so its first nd + np are that many different block groups, drawn without replacement: no two drivers or passengers share a center.
As in Ex. 1, the points come first and the streets after: Code block 31 puts a box around the 22 points with mapbbox, downloads its streets with osm_roads and keeps those inside it with cropnetwork. The cost of assigning a driver to a passenger is the driving time from the driver to the passenger, read in that direction, since with one-way streets the trip back would differ:
# Code block 31: the 22 points' box, its streets, and driving minutesdrv =findall(==("driver"), U.ROLE) # rows of U: 12 driverspas =findall(==("passenger"), U.ROLE) # and 10 passengersbbU =mapbbox(U.LON, U.LAT; xexpand =0.1, yexpand =0.1) # their boxNU, LU =osm_roads((bbU[1]..., bbU[2]...); cache_dir ="data")NU, LU =cropnetwork(NU, LU, U.LON, U.LAT; xexpand =0.1, yexpand =0.1)NU, LU =addconnectors(NU, LU, U.LON, U.LAT; add_nf_nf =false)cU =isequal.(LU.SOURCE, "CONNECTOR")tU = [c ? 60d /20:60d / v # minutes, as in Code block 6 for (c, d, v) inzip(cU, LU.DIST, LU.SPEED)]fw, bw = LU.DIR .!=-1, LU.DIR .!=1gU =SimpleWeightedDiGraph([LU.SRC[fw]; LU.DST[bw]], [LU.DST[fw]; LU.SRC[bw]], [tU[fw]; tU[bw]]; combine = min)TU, PU =shortestpaths(gU, nrow(U)) # among all 22 pointsC = TU[drv, pas] # driver i to passenger jround.(C, digits =1)
[ Info: Loading cached OSM data: data\osm_w-78.75_e-78.55_s35.7_n35.85_nodes.csv
The assignment is a snapshot: it assumes that every driver is free at the moment the ten requests are assigned together, the assumption to state whenever its result is read, as the Assumptions check of Lecture 1.2 asks. With more drivers than passengers, two drivers are left without a passenger, which row (a) of Eq. 1 allows. The greedy procedure of Ex. 2(b) and the linear program of Model 1 run unchanged on the 12 × 10 matrix:
# Code block 32: greedy against the linear program, in minutesX =greedyassignment(C)TC =sum(C .* X)TCᵒ, Xᵒ =LAP(C)@showround(TC, digits =2)@showround(TCᵒ, digits =2)# The passengers the two assign differently, and their driversd =findall(j -> X[:, j] != Xᵒ[:, j], 1:length(pas))car(Z, j) =findfirst(==(1), Z[:, j]) # the driver Z gives passenger jprt(DataFrame("Passenger"=> U.ID[pas[d]],"Greedy"=> [U.ID[drv[car(X, j)]] for j in d],"Greedy min"=> [round(C[car(X, j), j], digits =2) for j in d],"LP"=> [U.ID[drv[car(Xᵒ, j)]] for j in d],"LP min"=> [round(C[car(Xᵒ, j), j], digits =2) for j in d]))
Unter’s 3 × 3 matrix of Ex. 2 was chosen to make greedy look especially bad. On small instances drawn at random, with three cars and three passengers at random points, greedy and the linear program tend to come out much the same, and the difference is rarely as large as Ex. 2’s.
On Raleigh’s streets the same holds (Fig. 4): greedy is 1.7% above the optimum, and the two differ on only 3 of the ten passengers, where greedy’s first choices, each the least time still available, leave later passengers with longer pickups.
What the assignment saves can be put in the units a ride-hailing service manages, drivers. By Little’s Law of Lecture 1.3, Sec. 3, WIP = TH \times CT: with requests arriving at a throughput of TH an hour and each pickup taking a cycle time of CT, on average, the number of drivers on the way to a pickup at any time, the work in process, is WIP. The request rate below is an estimate, 20 an hour, set once so that it can be changed. A ride’s time is taken from the data: the average driving time between two of the 22 drawn points, as if each passenger rode to another of them. The drivers busy, picking up or carrying a passenger, are then TH times the pickup and ride together:
# Code block 33: Little's Law, drivers busy at TH requests an hourTH =20# requests per hour (an estimate; change it)CT = TC /length(pas) # average pickup minutes, greedyCTᵒ = TCᵒ /length(pas) # and by the linear programnU =nrow(U) # the 22 drawn pointsride =sum(TU[1:nU, 1:nU]) / (nU * (nU -1)) # minutes between two, avg# CT in minutes; WIP = TH*CT/60 drivers on the way to a pickup; and with# a ride, TH*(CT + ride)/60 drivers busy, of the 12@showround(ride, digits =2)ct = [CT, CTᵒ]prt(DataFrame("Assignment"=> ["Greedy", "LP"],"CT"=>round.(ct, digits =2),"WIP"=>round.(TH .* ct ./60, digits =2),"Busy"=>round.(TH .* (ct .+ ride) ./60, digits =2)))
At 20 requests an hour, about 1.2 driver is on the way to a pickup at any time, and the linear program’s shorter pickups free only 0.02 of a driver. The rides take more of the fleet’s time than the pickups, 6.5 min against 3.6: with them, about 3.4 of the 12 drivers are busy whichever assignment is used, so the fleet is sized by the rides, and a better assignment matters more as requests rise and pickups lengthen.
Show the code that draws this map
# Code block 34: the assignment on the streetsorc = colorant"#e08214"# greedy, where it differsstreetU(k) =SimpleWeightedGraph(LU.SRC[k], LU.DST[k], LU.DIST[k]; combine = min)lo =streetU(.!cU .& (LU.FCLASS .>=6)) # as in Code block 9hi =streetU(.!cU .& (LU.FCLASS .<=5))route(i, j) =tracepath(PU[drv[i]], drv[i], pas[j]) # driver i to jrts = [[route(car(Xᵒ, j), j) for j in1:length(pas)]; # LP, then [route(car(X, j), j) for j in d]] # greedy'sfig, ax =makemap(bbU[1], bbU[2]; doRoadbkgd =false, # the whole box xexpand =0, yexpand =0)resize!(fig, 687, 640)lines!(ax, x2ln(lo, NU.LON), x2ln(lo, NU.LAT); color = (ink, 0.12), linewidth =0.5)lines!(ax, x2ln(hi, NU.LON), x2ln(hi, NU.LAT); color = (ink, 0.3), linewidth =1)np =length(pas)for p in rts[1:np] # the linear program's routeslines!(ax, NU.LON[p], NU.LAT[p]; color = (dcc, 0.85), linewidth =3)endfor p in rts[np+1:end] # greedy's, where they differlines!(ax, NU.LON[p], NU.LAT[p]; color = (orc, 0.9), linewidth =2.5, linestyle =:dash)endxl, yl = bbU # the streets of the whole boxfor (nm, f, side) in (("Hillsborough Street", 0.72, -1), ("Wade Avenue", 0.5, 1), ("Western Boulevard", 0.85, -1))namestreet!(ax, NU, LU, nm, xl, yl; f, side) # Code block 9'send# Two places: NC State, below its two buildings of Ex. 1, clear of the# routes, and downtown, at the State Capitol of Code block 30text!(ax, [mean(x), c0[1]], [mean(y), c0[2]]; text = ["NC State", "Downtown"], fontsize =14, font =:bold, color = (ink, 0.6), align = (:center, :center), offset = [(-30, -95), (0, 0)])sD =scatter!(ax, U.LON[drv], U.LAT[drv]; marker =:rect, color = bluec, markersize =13, strokecolor =:white, strokewidth =1)sP =scatter!(ax, U.LON[pas], U.LAT[pas]; color = supc, markersize =13, strokecolor =:white, strokewidth =1)# Each ID above right of its point, unless a route leaves that waypl =Dict("D1"=> ((:right, :center), (-9, 0)),"D3"=> ((:right, :top), (-6, -4)),"D10"=> ((:right, :bottom), (-6, 4)),"D11"=> ((:right, :bottom), (-6, 4)),"P3"=> ((:right, :bottom), (-6, 4)),"P9"=> ((:right, :top), (-6, -4)))ap = [get(pl, id, ((:left, :bottom), (6, 4))) for id in U.ID]text!(ax, U.LON, U.LAT; text = U.ID, fontsize =11, color = ink, align =first.(ap), offset =last.(ap))lL =LineElement(color = dcc, linewidth =3)lG =LineElement(color = orc, linewidth =2.5, linestyle =:dash)axislegend(ax, [sD, sP, lL, lG], ["Driver", "Passenger", "LP","Greedy, where it differs"]; position =:rt, labelsize =12)fig
The linear program’s pickups total 35.5 min against greedy’s 36.2; at 20 requests an hour, about 1.2 driver is on the way to a pickup at any time
3. Transportation problem
The transportation problem is a generalization of the linear assignment problem (LAP). As in the LAP, nodes are separated into supply and demand nodes with arcs representing the associated cost per unit of flow, connecting each supply node to all of the demand nodes. Unlike the LAP, in which each supply node can provide only one unit and each demand node can require only one unit, the transportation problem supply nodes can provide multiple units, and the demand nodes can require multiple units.
The aim is to satisfy each node’s demand from the supply nodes. The transportation problem can be used for allocation in the alternating location–allocation (ALA) procedure of Lecture 2.3 when the new facilities (NFs) have capacity constraints, and the min cost/distance allocation of Lecture 2.3 is the transportation problem with infinite supply at each node.
It may seem presumptuous to call this the transportation problem, as if there were no other ways of transporting things, but that has always been its name. It goes back to the early days of linear programming, and for a long time afterward textbooks of logistics spent a good part of a chapter on solving these problems by hand, a tedious business of large tables of numbers. Unlike the linear assignment problem, in which each supply node has one unit of supply and each demand node demands one unit, a supply in the transportation problem can be any number, as can a demand. So the requirement for a feasible problem is not that there be more supply nodes than demand nodes: it is that the total supply equal or exceed the total demand, however few the supply nodes. If it does not, the problem is infeasible, and solved as a linear program it returns a message saying so. Because of their simple structure, these problems have traditionally been laid out as a tableau: a cost matrix with a row for each supply node and a column for each demand node, the arc costs in its cells, the supplies down the right side and the demands along the bottom. Feasibility can be checked from the tableau by adding up the supply column and the demand row. The transportation problem can be modeled as a network flow problem and then solved as an LP, as in Model 2.
minimize: total cost of shipping from the supply nodes to the demand nodes
solve for:
(a) amount shipped from each supply node to each demand node.
subject to:
(a) supply: each supply node provides at most one unitits supply;
(b) demand: each demand node requires exactly one unitits demand.
assumptions:
(a) cost per unit shipped from each supply node to each demand node is constant;
(b) total supply is at least total demand.
=amount shipped from supply node i to demand node j (units)
c_{ij}
=cost per unit shipped from i to j (cost/unit)
s_i
=supply at node i (units)
d_j
=demand at node j (units)
S
=supply nodes
D
=demand nodes.
Model 2 states two constraints, and both are rows, the rows of Eq. 1 with their right-hand sides 1 replaced by the supplies s_i and the demands d_j. Row (a), supply, lets a supply node ship less than its supply, so total supply can exceed total demand and the excess is simply not shipped; row (b), demand, meets each demand exactly. Assumption (b) is what makes the rows feasible together, and it is the throughput feasibility condition of Lecture 1.3, Sec. 5, in network form: no plan can ship more than the supply nodes provide.
Model 2 implementation: Transportation
# Model: transportationfunctiontrans(c, s, d) S, D =1:length(s), 1:length(d)ifsum(s) <sum(d) # assumption (b)error("total supply is less than total demand")end model =Model(HiGHS.Optimizer)@variable(model, 0≤ x[S,D])@objective(model, Min, sum(c[i,j]*x[i,j] for i ∈ S, j ∈ D))@constraint(model, supply[i ∈ S], sum(x[i,j] for j ∈ D) ≤ s[i]) # (a)@constraint(model, demand[j ∈ D], sum(x[i,j] for i ∈ S) == d[j]) # (b)set_silent(model)optimize!(model)iftermination_status(model) == JuMP.OPTIMALreturn (TC =objective_value(model), x =snapvals(value.(x)))enderror(termination_status(model))end
1
As in LAP, the model’s assumption is checked before the linear program is built: assumption (b) of Model 2, that total supply is at least total demand. When it fails, error stops the function with its message, as in Ex. 4(c), rather than leaving the solver to report the case as infeasible.
trans (generic function with 1 method)
Example 4: Transportation with excess supply
Supply nodes 1 to 3 of Fig. 5 ship to demand nodes 4 to 7 at the costs per unit of Table 3. Determine the shipments that meet every demand at least total cost, first by the greedy procedure and then as Model 2, and determine what happens when the total supply is less than the total demand.
Show the code that draws this figure
# Code block 35: the transportation networksupc, demc = colorant"#d21f26", colorant"#1e8a3c"ink = colorant"#252525"r =0.2# node radius, data unitsXYs = [(0.0, 2.2), (0.0, 1.1), (0.0, 0.0)]XYd = [(2.4, 2.4), (2.4, 1.6), (2.4, 0.8), (2.4, 0.0)]sup, dem = [55, 50, 40], [45, 20, 30, 30]cst = [86109; 912137; 149165] # the tableau's costs per unit# Each cost near its arc's tail: the fraction of the way along it, and# the side of it (1 to the left of the arc's direction) clear of othersat = [0.120.120.150.12; 0.120.120.150.12; 0.210.150.120.12]sd = [111-1; 1-1-1-1; -1-1-1-1]fig =Figure(size = (440, 390))ax =Axis(fig[1, 1]; aspect =DataAspect())hidedecorations!(ax); hidespines!(ax)# Every arc, outline to outline, and its costfor (a, p) inenumerate(XYs), (b, q) inenumerate(XYd) u = (q .- p) ./hypot((q .- p)...) t, h = p .+ r .* u, q .- (r +0.015) .* uarrows2d!(ax, [t[1]], [t[2]], [h[1] - t[1]], [h[2] - t[2]]; color = ink, shaftwidth =1.2, tipwidth =8, tiplength =8) m = t .+ at[a, b] .* (h .- t) .+0.08* sd[a, b] .* (-u[2], u[1])text!(ax, m[1], m[2]; text =string(cst[a, b]), color = ink, fontsize =10, align = (:center, :center))endfor (k, (x, y)) inenumerate([XYs; XYd]) # the nodes, over the arcs c = k <=3 ? supc : demcpoly!(ax, Circle(Point2f(x, y), r); color =:white, strokecolor = c, strokewidth =1.8)text!(ax, x, y; text =string(k), color = c, fontsize =15, align = (:center, :center)) amt, dx = k <=3 ? (sup[k], -0.45) : (dem[k -3], 0.45)text!(ax, x + dx, y; text =string(amt), color = c, fontsize =14, align = (:center, :center))end# The supply and the demand columns, titledtext!(ax, [-0.45, 2.85], [2.55, 2.75]; text = ["Supply", "Demand"], color = [supc, demc], fontsize =13, align = (:center, :bottom))xlims!(ax, -0.9, 3.4); ylims!(ax, -0.4, 3.1)fig
Figure 5: The transportation network of Ex. 4: each supply node (red) is joined to every demand node (green), its supply or demand beside it under the column’s title, and each arc is labeled with its cost per unit, near its supply node.
Table 3: Ex. 4’s tableau: the cost per unit from each supply node (row) to each demand node (column), with each node’s supply and demand.
Trans
4
5
6
7
Supply
1
8
6
10
9
55
2
9
12
13
7
50
3
14
9
16
5
40
Demand
45
20
30
30
Example 4(a): Greedy by hand
Ship from the supply nodes to the demand nodes by repeatedly making the shipment of least cost per unit still possible, as much as it can carry, and determine the total cost.
As was done for the LAP, a greedy procedure can be used to solve the transportation problem. The procedure is simply to
Greedy Transportation Procedure: Continue to select the lowest cost supply until all demand is satisfied.
The greedy procedure is fast but not always optimal for the transportation problem. By contrast, Dijkstra’s shortest path procedure of Lecture 4.1 and the simplex method for LP are optimal greedy procedures.
The greedy approach picks the smallest cost and assigns as much as it can on that arc, then looks for the next least cost and assigns as much as it can of that, and keeps going until all of the demand has been satisfied from the supply. It is called greedy because it always takes the next least cost. It is straightforward, but it does not guarantee an optimal solution; its answer is sometimes pretty good, but it is not optimal. Were the optimum hard to find, as it is for the uncapacitated facility location problem of Lecture 2.4, a heuristic would be welcome, and there it was: Lecture 2.7 found the heuristic within a small fraction of one percent of the optimum of the mixed-integer linear program. Here, even on a small problem, the greedy approach is a lot worse than the optimum, and the optimum is easy to get because the problem is a linear program, which today’s solvers handle at sizes far beyond any example here, and so there is no reason to accept greedy’s answer.
# Code block 36: the costs, demands and supplies of the tableauC = [86109; 912137; 149165]d = [45, 20, 30, 30]s = [55, 50, 40]@showsum(s)@showsum(d)sum(s) -sum(d)
sum(s) = 145
sum(d) = 125
20
By hand, on Table 3, with each step drawn in Fig. 6, the procedure starts with the minimum cost, 5 from node 3 to node 7, and ships all 30 that node 7 demands, which node 3 can do because it has 40, leaving it 10. The next smallest cost, 7 from node 2 to node 7, is no longer available, since node 7’s demand has been satisfied and its whole column is eliminated. So the next is 6, from node 1 to node 5, which satisfies all 20 of node 5’s demand and leaves node 1 with 35; the columns of nodes 5 and 7 are now eliminated. Of the two columns left, the next smallest cost is 8, from node 1 to node 4, but node 1 cannot ship node 4’s 45 because it has only 35, so it ships 35, its supply goes to zero and its row is eliminated, leaving node 4 ten units short. The 9 from node 2 to node 4 finishes node 4’s demand, and the 13 from node 2 to node 6 then ships all 30 of node 6’s demand. That is the greedy solution, and its total cost is the sum of each cost times the units shipped at it.
The procedure is called greedy because it only ever does what is immediately best: it always takes the smallest cost still available, and it does not look ahead to what a choice leaves for later. On this kind of problem that is not the best thing to do, and looking one or two steps ahead would have done better. There is one greedy procedure that does give the optimal solution to both the linear assignment and the transportation problems, and that is linear programming itself: the simplex method is greedy in that it always makes the most immediately favorable move. So the claim is not that greedy never works, only that this greedy procedure does not.
# Code block 37: the greedy procedure for the transportation problemfunctiongtrans(C, s, d) steps =DataFrame(From =Int[], To =Int[], Cost =eltype(C)[], Units =eltype(s)[]) s, d =float(copy(s)), float(copy(d)) F =zeros(size(C)) TC =0.0whileany(s .>0) &&any(d .>0) idxs, idxd =findall(s .>0), findall(d .>0) k =argmin(C[idxs, idxd]) # least cost still available i, j = idxs[k[1]], idxd[k[2]] f =min(s[i], d[j]) # as much as it can carry F[i, j] = f TC += C[i, j] * f s[i] -= f d[j] -= fpush!(steps, (i, length(s) + j, C[i, j], f))endreturn F, TC, stepsendFg, TCg, steps =gtrans(C, s, d)prt(steps)
1
Copies, as greedyassignment copies C in Code block 23: the loop uses up s and d, and the caller’s are left as they were for Code block 40, which solves the same s and d next.
2
The least cost is sought only among the rows with supply left and the columns with demand left, so k is a position in that smaller matrix, and idxs[k[1]] and idxd[k[2]] turn it back into a row and a column of C.
From To Cost Units
────────────────────────
1 3 7 5 30
2 1 5 6 20
3 1 4 8 35
4 2 4 9 10
5 2 6 13 30
TCg # Code block 38: total cost of greedy
1030.0
Show the code that draws this figure
# Code block 39: the steps of greedy, each frame from its row of steps# The supply and demand left after steps 1 to k, and the units shippedfunctionreplay(C, s, d, steps, k) m =size(C, 1) F, sl, dl =zeros(Int, size(C)), copy(s), copy(d)for t in1:k # steps 1 to k, replayed a, b, f = steps.From[t], steps.To[t] - m, steps.Units[t] F[a, b] += f sl[a] -= f dl[b] -= fendreturn F, sl, dlend# Value v at (x, y), struck through: the value before the step. The# stroke is a thin filled band rather than a line, because the figure# label check counts any line across a label as a defect, and here# crossing the label is the pointfunctionstruck!(ax, x, y, v, c)text!(ax, x, y; text =string(v), color = (c, 0.55), fontsize =12, align = (:center, :center)) w, h =0.1*length(string(v)) +0.06, 0.018# half width, thickness band = Point2f[(x - w, y -0.1- h), (x + w, y +0.1- h), (x + w, y +0.1+ h), (x - w, y -0.1+ h)]poly!(ax, band; color = c)returnnothingend# The tableau after step k (k = 0, as it starts): the costs, the units# shipped, the supply and demand left, what is used up shaded, step k's# cell in blue and the supply and demand it changed struck throughfunctiongreedystep!(ax, C, s, d, steps, k) m, n =size(C) F, sl, dl =replay(C, s, d, steps, k) _, sp, dp =replay(C, s, d, steps, max(k -1, 0)) # before step k ak, bk = k ==0 ? (0, 0) : (steps.From[k], steps.To[k] - m)for a in1:m, b in1:n # cell (a,b) from (b-1, m-a) here = (a, b) == (ak, bk) used = sl[a] ==0|| dl[b] ==0 fc = here ? (bluec, 0.18) : used ? (ink, 0.1) : (:white, 0.0)poly!(ax, Rect2f(b -1, m - a, 1, 1); color = fc)text!(ax, b -0.1, m - a +0.9; text =string(C[a, b]), color = used && !here ? (ink, 0.55) : ink, fontsize =11, align = (:right, :top)) F[a, b] >0||continuetext!(ax, b -0.55, m - a +0.35; text =string(F[a, b]), color = here ? bluec : ink, fontsize =13, font =:bold, align = (:center, :center))endlinesegments!(ax, [Point2f(x, y) for x in0:n for y in (0, m)]; color = ink, linewidth =0.8)linesegments!(ax, [Point2f(x, y) for y in0:m for x in (0, n)]; color = ink, linewidth =0.8) k ==0||poly!(ax, Rect2f(bk -1, m - ak, 1, 1); strokewidth =2, color = (:white, 0.0), strokecolor = bluec)for a in1:m # supply nodes and supply lefttext!(ax, -0.35, m - a +0.5; text =string(a), color = supc, fontsize =13, align = (:center, :center)) y = m - a +0.5if a == ak # struck, and the new valuestruck!(ax, n +0.4, y, sp[a], supc)text!(ax, n +1.05, y; text =string(sl[a]), color = supc, fontsize =12, font =:bold, align = (:center, :center))elsetext!(ax, n +0.4, y; text =string(sl[a]), color = supc, fontsize =12, align = (:center, :center))endendfor b in1:n # demand nodes and demand lefttext!(ax, b -0.5, m +0.35; text =string(m + b), color = demc, fontsize =13, align = (:center, :center))if b == bk # struck, the new value belowstruck!(ax, b -0.5, -0.35, dp[b], demc)text!(ax, b -0.5, -0.85; text =string(dl[b]), color = demc, fontsize =12, font =:bold, align = (:center, :center))elsetext!(ax, b -0.5, -0.35; text =string(dl[b]), color = demc, fontsize =12, align = (:center, :center))endendxlims!(ax, -0.7, n +1.45); ylims!(ax, -1.15, m +0.7)returnnothingendfig =Figure(size = (687, 470))for k in0:nrow(steps) # the start, then each step title = k ==0 ? "Start: the costs, supply and demand":"Step $k: $(steps.Cost[k]) from $(steps.From[k]) to "*"$(steps.To[k]), $(steps.Units[k]) units" axk =Axis(fig[k ÷3+1, k %3+1]; aspect =DataAspect(), title, titlesize =12, titlefont =:regular)hidedecorations!(axk); hidespines!(axk)greedystep!(axk, C, s, d, steps, k)end# Below the frames: each step's cost times its units, and their sumLabel(fig[3, 1:3], "Total cost: "*join(["$(st.Cost) × $(st.Units)" for st ineachrow(steps)], " + ") *" = "*commas(Int(TCg)); fontsize =13, color = ink)colgap!(fig.layout, 6); rowgap!(fig.layout, 4)fig
Figure 6: The greedy procedure on Ex. 4’s tableau: the tableau as it starts, then one frame per step, each drawn from the step gtrans records. Each cell holds its cost per unit at its corner and, once units are shipped on it, the units in bold. Each supply node (red) heads a row that ends in the supply it has left, and each demand node (green) heads a column that ends in the demand it has left; a row or column used up is shaded. The cell chosen at the step is blue, and the supply and demand it changed are struck through, with the value left beside or below them in bold. Below the frames, the total cost: each cost times its units, summed.
Example 4(b): As a linear program
Determine the shipments of least total cost as Model 2, and compare its total cost with greedy’s.
# Code block 40: the tableau solved as a linear programTCᵒ, Xᵒ =trans(C, s, d)@show TCᵒ; Xᵒ
prt labels a matrix’s rows and columns with rows and cols, and its corner with row_title; a NaN prints as a blank, for the corner cells that hold no total.
A total cost in Eq. 3, greedy or optimal, is calculated the same way: each cell of the tableau that ships anything contributes its cost per unit times the units it ships, and the total is the sum of those products, over the five cells greedy uses and the five the linear program uses.
The shipments of Code block 41 balance: each supply node ships no more than its supply (in its row, Shipped against Supply), and each demand node receives exactly its demand (in its column, Shipped against Demand); and the 145 units of supply less the 125 shipped leaves 20 unshipped, 10 at node 2 and 10 at node 3. That the parts add up to the whole is the Balance check of Lecture 1.2.
Both totals can also be bracketed. Greedy’s shipments are feasible, so no optimum can cost more than its 1,030. And no plan can cost less than one in which every demand node takes all of its demand at its cheapest arc, as if every supply were infinite, which is the minimum cost allocation:
# Code block 42: every demand at its cheapest arc, supply unlimitedLB =sum(minimum(C, dims =1)' .* d)
930
The optimum, 970, lies between 930 and 1,030, as it must: the Bounds check of Lecture 1.2.
Greedy’s 1,030 is significantly worse than the linear program’s optimal 970 (Eq. 3), about 6% more, and for such a small problem. That is not the case of the ADD procedure for the uncapacitated facility location problem, where the heuristic really is pretty good; here the difference is significant enough that the linear program is definitely wanted. The good news is that the linear program is very effective and easy to implement, so using it is no burden.
Example 4(c): Too little supply
Determine what the linear program returns when node 1’s supply falls from 55 to 25.
What if the total supply is less than the total demand?
sum(s) - sum(d) = -10
total supply is less than total demand
The total supply is now 115 against a demand of 125, 10 short, and trans stops with its error before any linear program is built: assumption (b) of Model 2 fails. Without the check, the linear program would be built and solved and come back INFEASIBLE, since no shipments can satisfy rows (a) and (b) of Eq. 2 together. The assumption is the throughput feasibility condition of Lecture 1.3, Sec. 5: the supply nodes together cannot provide the throughput the demand nodes require. Greedy, given the same data, would ship everything it has and stop, leaving a demand unmet and reporting nothing, the silent failure of Ex. 2(b).
Greedy ships at a total cost of 1,030, about 6% more than the linear program’s 970; with 10 units too little supply there is no feasible plan, and trans says so
The allocation step of the ALA procedure in Lecture 2.3, each customer to its nearest NF, can be considered uncapacitated allocation, since no NF is limited in what it can serve; capacitated allocation is its alternative. When NFs have capacity constraints, the allocation for the ALA procedure can be determined by solving the transportation problem. In Lecture 2.3, Ex. 2(a), the population of NC 100K cities was allocated to Raleigh and Charlotte based on closeness. All of the population of each city was allocated to either Raleigh or Charlotte. If the total population allocated was constrained, the population of some of the cities would have to be split between Raleigh and Charlotte. In the following, the capacity of Raleigh and Charlotte is 1.5 million people.
From a practical logistics point of view, the transportation problem comes up in the ALA procedure. Imagine a company with an existing distribution center (DC) that was built to a certain capacity and already runs three shifts a day, so that there is no way to add a shift to handle any increase: it is maxed out. If the company considers adding a second DC, or closing the existing one and adding one or two more, the allocation step of ALA as Lecture 2.3 ran it, every customer to its nearest DC, cannot respect that capacity. Solved as a transportation problem instead, the allocation can give the existing DC its real capacity, and each new facility a capacity equal to the total demand, which is effectively infinite, so that a new facility can be built to whatever size is wanted. The existing DC’s fixed location is handled as before, by always returning its current location from the location step; with the transportation problem as the allocation step, it can have a fixed capacity as well, and the rest of the ALA procedure is unchanged. The transportation problem is then a black box inside ALA: given the cost matrix, the supplies and the demands, it is translated into a linear program and solved, and as a linear program its answer is optimal, within the limits of how accurately a linear program models the problem.
Example 5: Capacitated allocation
The North Carolina cities over 100,000 people are served by DCs at Charlotte and Raleigh, as in Ex. 2(a) in Lecture 2.3. Determine the population-weighted total distance when each DC can serve at most 1.5 million people, which cities are then split between the two DCs, and what the capacity costs.
Example 5(a): Uncapacitated as in Lecture 2.3
Allocate each city to its nearer DC, and determine the population each DC serves.
# Code block 44: the NC cities over 100,000, and the two DCsdf =filter(r -> r.STFIP ==st2fips(:NC) && r.POP >100_000, usplace())select!(df, :NAME, :LAT, :LON, :POP)P =hcat(df.LON, df.LAT) # one row per cityw = df.POPname2lonlat(nm) =collect(filter(r -> r.NAME == nm, df)[1, [:LON, :LAT]])X =vcat(name2lonlat("Charlotte")', name2lonlat("Raleigh")')D =dists(X, P, :mi) # 2 x m great circle distancesprt(df)
NAME LAT LON POP
───────────────────────────────────────
Cary 35.78 -78.82 174,721
Charlotte 35.21 -80.83 874,579
Concord 35.39 -80.64 105,240
Durham 35.98 -78.90 283,506
Fayetteville 35.08 -78.97 208,501
Greensboro 36.10 -79.83 299,035
High Point 35.99 -79.99 114,059
Raleigh 35.83 -78.64 467,665
Wilmington 34.21 -77.89 115,451
Winston-Salem 36.10 -80.26 249,545
# Code block 45: each city to its nearer DC, uncapacitatedusingSparseArraysn, m =size(D)alloc(X) =sparse([argmin(c) for c ineachcol(dists(X, P, :mi))],1:m, w, n, m)W =alloc(X)'@showsum(W'.* D) # person-mi, as in Lecture 2.3sum(W, dims=1) # people served by each DC
sum(W' .* D) = 8.004697649736136e7
1×2 SparseMatrixCSC{Int64, Int64} with 2 stored entries:
1343423 1548879
Nearest allocation asks more of Raleigh than 1.5 million people, so with the capacity in place some of its population has to be served from Charlotte.
Example 5(b): Each DC at 1.5 million
Allocate the cities’ populations to the two DCs, each limited to 1.5 million, as Model 2, and compare with greedy.
# Code block 46: capacity constrained, as a transportation problems = [1_500_000, 1_500_000]d = walloc(X) =trans(dists(X, P, :mi), s, d)TCc, Wc =alloc(X)@show TCc@show TCc -sum(W'.* D) # what the capacity costsprt(DataFrame(NAME = df.NAME, POP = w, Charlotte = Wc[1, :], Raleigh = Wc[2, :]))
TCc = 8.07528287304478e7
TCc - sum(W' .* D) = 705852.233086437
NAME POP Charlotte Raleigh
────────────────────────────────────────────
Cary 174,721 0 174,721
Charlotte 874,579 874,579 0
Concord 105,240 105,240 0
Durham 283,506 0 283,506
Fayetteville 208,501 0 208,501
Greensboro 299,035 48,879 250,156
High Point 114,059 114,059 0
Raleigh 467,665 0 467,665
Wilmington 115,451 0 115,451
Winston-Salem 249,545 249,545 0
Show the code that draws this map
# Code block 47: the capacitated allocation, mappedsp = [k for k in1:m if count(>(0), Wc[:, k]) >1] # the split cityfig, ax =makemap(df.LON, df.LAT; xexpand =0.2)resize!(fig, 687, 500)for (k, col) in ((1, dcc), (2, supc)) # Charlotte, Raleigh Lx, Ly =alloclines(Wc[k:k, :], X[k:k, :], P)for i ineachindex(Lx)lines!(ax, Lx[i], Ly[i]; color = col, linewidth =1)endendscatter!(ax, P[:, 1], P[:, 2]; color = ink, markersize =6)scatter!(ax, X[:, 1], X[:, 2]; marker =:dtriangle, markersize =14, color = ink)lx, ly = [X[:, 1]; P[sp, 1]], [X[:, 2]; P[sp, 2]]text!(ax, lx, ly; text = ["Charlotte"; "Raleigh"; df.NAME[sp]], fontsize =12, # each clear of its lines and the border align = [(:right, :bottom), (:left, :center), (:center, :bottom)], offset = [(-6, 4), (9, 0), (0, 5)])fig
Figure 7: The allocation of Ex. 5(b), each DC serving a population of at most 1.5 million: each city is joined to the DC that serves its population, Charlotte’s in green and Raleigh’s in red, and Greensboro, whose population is split between the two DCs, is joined to both.
# Code block 48: the same allocation by the greedy procedure of Ex. 4Fg5, TCg5, _ =gtrans(D, s, d)@show TCg5@show TCg5 - TCc # what greedy costs over the optimumsplitg = df.NAME[[k for k in1:m if count(>(0), Fg5[:, k]) >1]]prt(DataFrame(NAME = df.NAME, POP = w, Charlotte = Fg5[1, :], Raleigh = Fg5[2, :]))
TCg5 = 8.303312161115827e7
TCg5 - TCc = 2.2802928807104677e6
NAME POP Charlotte Raleigh
────────────────────────────────────────────
Cary 174,721 0 174,721
Charlotte 874,579 874,579 0
Concord 105,240 105,240 0
Durham 283,506 0 283,506
Fayetteville 208,501 0 208,501
Greensboro 299,035 0 299,035
High Point 114,059 114,059 0
Raleigh 467,665 0 467,665
Wilmington 115,451 48,879 66,572
Winston-Salem 249,545 249,545 0
# Code block 49: the cost as each DC's capacity movescaps = [1_450_000, 1_500_000, 1_600_000]prt(DataFrame(Capacity = caps, PersonMiles = [trans(D, [c, c], d).TC for c in caps]))
The linear program splits one city, Greensboro, whose people are almost as near Charlotte as Raleigh (Fig. 7), and greedy, which fills each DC from its cheapest arcs first, splits Wilmington instead and costs more. As each DC’s capacity rises, the cost falls, and once a capacity covers what nearest allocation asks of it, the capacity costs nothing: the answer moves the way it should, the Nudge check of Lecture 1.2. The uncapacitated total is Lecture 2.3’s exactly, the same Census table for the same cities, the Source check of Lecture 1.2.
At 1.5 million people each, the allocation costs 80,752,829 person-mi against 80,046,976 uncapacitated, 0.9% more, with one city split between the two DCs; greedy costs 2,280,293 person-mi more than the linear program
4. Minimum cost network flow
Minimum cost network flow (MCNF) is the most general network model. It can be used to represent a variety of network models that can be solved as an LP. Special procedures more efficient than LP were developed to solve MCNF and Transportation problems, e.g.:
Network simplex algorithm (MCNF)
Hungarian method (Transportation and Transshipment)
It is usually easier to transform into an LP since current LP solvers are so good, with MCNF just aiding in the formulation of the problem. Special, very efficient procedures are only used for the shortest-path problem (e.g., Dijkstra).
Each row i of the arc-incidence matrix and the corresponding net supply represent the coefficients of the flow-balance equation for node i. The LP formulation of the MCNF problem is Model 3.
If all that were wanted were the transportation problem, there would be little reason for network flow problems, but they get a lot more interesting with additional capabilities. The more general problem, of which the transportation problem is just a special case, is the minimum cost network flow problem, and the good news is that it is still a linear program, so the ability to solve it effectively is not lost, while the representation handles a much bigger variety of scenarios. In the network of Fig. 8, supply nodes can be anywhere in the network: nodes 1 and 2 supply 6 and 2 units, and node 5 demands 8, so 8 units go into the network and 8 come out, and everything balances. Nodes 3 and 4 are transshipment nodes: whatever comes into one goes out. In a logistics setting a transshipment node would be like a distribution center, where the pallets that go in eventually leave, while the supply nodes produce pallets and the demand nodes consume them. Each arc carries a cost.
The difference from the transportation problem is that there, one level is clearly supply and the other demand, while in this general formulation nothing distinguishes one kind of node from another. So, by convention, each node is given a single number, its net supply: positive at a supply node, negative at a demand node, and zero at a transshipment node. Capacities can be added, both on the arcs and at the nodes; without them, as here, the only information is the cost. The costs form a vector, the objective of the linear program, one element per arc, and the arcs are the columns of the incidence matrix, whose rows are the nodes: the first arc goes from node 1 to node 2, so its column has +1 in row 1, where it starts, and -1 in row 2, where it ends. With the costs, the incidence matrix and the net supplies, the problem is ready to be solved as a linear program, one that minimizes total cost.
Show the code that draws this figure
# Code block 50: drawing a network with its net supplies# Nodes at XY; arc k from I[k] to J[k] labeled lab[k] in size ls, at# fraction at[k] along it, on side side[k], the arc shifted sideways by# off[k] (for parallel arcs) and drawn two-way if two[k]; net supply# s[i] beside node i, away from its arcs or at angle dirs[i], with its# role above it if roles (or the name roles[i]), and a zero shown if# showzero; notes, k => text, on the other side of arc k from its label;# and a legend at (x, y) reading an arc's label: given over defaultfunctiondrawnet(XY, I, J, lab, s; off =zeros(length(I)), r =0.17, w =380, h =260, lim =nothing, at =fill(0.5, length(I)), side = [o >=0 ? 1:-1 for o in off], ls =13, roles =false, showzero = roles !=false, two =falses(length(I)), rotate =false, notes =Dict{Int,String}(), legend =nothing, dirs =nothing) fig =Figure(size = (w, h)) ax =Axis(fig[1, 1]; aspect =DataAspect())hidedecorations!(ax); hidespines!(ax)for k ineachindex(I) p, q = XY[I[k]], XY[J[k]] u = (q .- p) ./hypot((q .- p)...) nrm = (-u[2], u[1]) t = p .+ r .* u .+ off[k] .* nrm hd = q .- (r +0.015) .* u .+ off[k] .* nrmarrows2d!(ax, [t[1]], [t[2]], [hd[1] - t[1]], [hd[2] - t[2]]; color = ink, shaftwidth =1.3, tipwidth =8, tiplength =8)if two[k] # two-way: a head at each end tl = p .+ (r +0.015) .* u .+ off[k] .* nrmarrows2d!(ax, [hd[1]], [hd[2]], [tl[1] - hd[1]], [tl[2] - hd[2]]; color = ink, shaftwidth =1.3, tipwidth =8, tiplength =8)end θ =atan(u[2], u[1]) # the arc's angle, kept readable θ <-π/2+0.01&& (θ +=π) θ >π/2+0.01&& (θ -=π) m = t .+ at[k] .* (hd .- t) .+0.01ls * side[k] .* nrmtext!(ax, m[1], m[2]; text = lab[k], color = ink, fontsize = ls, rotation = rotate ? θ :0.0, align = (:center, :center))ifhaskey(notes, k) # a note on the arc's other side n = t .+ at[k] .* (hd .- t) .-0.01ls * side[k] .* nrmtext!(ax, n[1], n[2]; text = notes[k], color = ink, fontsize = ls -1, rotation = θ, align = (:center, :center))endendfor (i, (x, y)) inenumerate(XY) # the nodes, over the arcs c = s[i] >0 ? supc : s[i] <0 ? demc : inkpoly!(ax, Circle(Point2f(x, y), r); color =:white, strokecolor = c, strokewidth =1.8)text!(ax, x, y; text =string(i), color = c, fontsize =14, align = (:center, :center)) s[i] ==0&& !showzero &&continue a =sum(collect(XY[j] .- (x, y)) ./hypot((XY[j] .- (x, y))...)for j in [J[I .== i]; I[J .== i]]) # toward its arcs a = dirs ===nothing ? a ./hypot(a...) :# or away from dirs[i] (-cosd(dirs[i]), -sind(dirs[i])) txt =replace(string(s[i]), "-"=>"−") # minus sign roles ==true&& (txt = (s[i] >0 ? "Supply": s[i] <0 ?"Demand":"Transshipment") *"\n"* txt) roles isa AbstractVector&& (txt = roles[i] *"\n"* txt) d = r + (roles !=false ? 0.27:0.16)text!(ax, x - d * a[1], y - d * a[2]; text = txt, color = c, fontsize =13, align = (:center, :center))endif legend !==nothing# how an arc's label reads x0, y0 = legendfor (x, t) in ((x0, "i"), (x0 +0.8, "j"))poly!(ax, Circle(Point2f(x, y0), 0.12); color =:white, strokecolor = ink, strokewidth =1.2)text!(ax, x, y0; text = t, color = ink, fontsize =12, font =:italic, align = (:center, :center))endarrows2d!(ax, [x0 +0.12], [y0], [0.8-0.255], [0.0]; color = ink, shaftwidth =1.2, tipwidth =7, tiplength =7)text!(ax, x0 +0.4, y0 +0.08; text ="c, u, l", color = ink, fontsize =12, align = (:center, :bottom))text!(ax, x0 +0.4, y0 -0.08; text ="c, ∞, 0", color = (ink, 0.6), fontsize =12, align = (:center, :top))end lim ===nothing|| (xlims!(ax, lim[1]...); ylims!(ax, lim[2]...))return figendXY5 = [(0.0, 1.0), (1.2, 1.8), (1.2, 0.2), (2.4, 1.8), (2.4, 0.2)]drawnet(XY5, [1, 1, 2, 2, 3, 4], [2, 3, 4, 5, 5, 5],string.([2, 3, 4, 5, 1, 3]), [6, 2, 0, 0, -8]; roles =true, notes =Dict(1=>"arc cost"), dirs = [180, 90, 270, 90, 270], w =440, h =320, lim = ((-0.8, 3.0), (-0.65, 2.55)))
Figure 8: A five-node network, each node labeled with its role and its net supply: nodes 1 and 2 supply 6 and 2 units (positive net supply), node 5 demands 8 (negative), and nodes 3 and 4 are transshipment nodes (zero); each arc is labeled with its cost per unit of flow.
minimize: total cost of the flow on the arcs of a network
solve for:
(a) flow on each arc.
subject to:
(a) flow balance: at each node, flow out minus flow in is its net supply, and no more than it at a supply node; (b) arc bounds: the flow on each arc is between its lower and upper bound.
assumptions:
(a) cost per unit of flow on each arc is constant;
(b) total supply is at least total demandeach node’s net supply is given: positive at a supply node, negative at a demand node, and zero at a node flow only passes through.
=\begin{cases} 1, & \text{if arc } k \text{ starts at node } i \\ -1, & \text{if arc } k \text{ ends at node } i \\ 0, & \text{otherwise} \end{cases}
s_i
=net supply of node i (units)
=\begin{cases} > 0, & \text{supply node } i \\ = 0, & \text{transshipment node } i \\ < 0, & \text{demand node } i \end{cases}
l_k
=lower bound on the flow on arc k, by default 0 (units)
u_k
=upper bound on the flow on arc k, by default \infty (units)
V
=nodes
E
=arcs.
Both lines labeled (a) in Eq. 4 are row (a), flow balance, on purpose: the one constraint is written as an inequality at the supply nodes and as an equation at every other node.
Model 3 states two constraints, and both are rows. Row (a), flow balance, is one row per node, written in two lines: the sum \sum_k a_{ik} x_k is the flow out of node i less the flow into it, which at a supply node can be no more than its supply, and elsewhere must equal the net supply, the negative of a demand at a demand node and zero at a transshipment node. At the supply and demand nodes of a transportation network, row (a) is exactly rows (a) and (b) of Eq. 2, which is why Model 2 is a special case of Model 3. Row (b), arc bounds, is new; with the defaults, 0 and \infty, it is just the nonnegativity of the flows.
Model 3 implementation: Minimum cost network flow
# Model: minimum cost network flow# Incidence matrix from an arc list: arc k from i[k] to j[k], +1 where it# starts and -1 where it ends, one column per arc in the list's orderlist2incid(i, j, n) =sparse([i; j], [1:length(i); 1:length(i)], [ones(length(i)); -ones(length(j))], n, length(i))functionmcnf(A, c, s; u =fill(Inf, size(A, 2)), l =zeros(size(A, 2))) V, E =1:size(A, 1), 1:size(A, 2) model =Model(HiGHS.Optimizer)@variable(model, l[k] ≤ x[k ∈ E] ≤ u[k]) # (b)@objective(model, Min, sum(c[k] * x[k] for k ∈ E))@constraint(model, supply[i ∈ V; s[i] >0],sum(A[i, k] * x[k] for k ∈ E) ≤ s[i]) # (a)@constraint(model, balance[i ∈ V; s[i] ≤0],sum(A[i, k] * x[k] for k ∈ E) == s[i]) # (a)set_silent(model)optimize!(model)iftermination_status(model) == JuMP.OPTIMALreturn (TC =objective_value(model), x =snapvals(value.(x)))enderror(termination_status(model))end
1
sparse(I, J, V, n, m), as in Lecture 2.3, puts V[k] at row I[k] and column J[k] of an n \times m matrix: here the rows are the arcs’ start nodes and then their end nodes, the columns the arc numbers twice, and the values +1 and then -1. Built from the arc list itself, every arc keeps its own column in the order listed, including the parallel arcs of Ex. 6(d), which a graph built from the same list would merge.
2
Row (b) is written as bounds on the variables rather than as constraints: each x[k] lies between its own l[k] and u[k], and an upper bound of Inf is no bound at all.
3
The condition after the semicolon keeps only the nodes that meet it: one row (a) for each supply node here, and one for every other node in the balance rows below.
mcnf (generic function with 1 method)
For the network of Fig. 8, the arc list gives the costs and the incidence matrix, one column per arc in the order listed, and the net supplies complete the data:
Labeled, the same data are Table 4, in which each column is an arc, with its cost above it, and each row is a node, with its net supply at the right:
Table 4: The data of Fig. 8: each column an arc, with its cost c_k, and each row a node, with its entries of the incidence matrix and its net supply s_i.
arc
(1,2)
(1,3)
(2,4)
(2,5)
(3,5)
(4,5)
c_k
2
3
4
5
1
3
s_i
node
1
1
1
0
0
0
0
6
2
-1
0
1
1
0
0
2
3
0
-1
0
0
1
0
0
4
0
0
-1
0
0
1
0
5
0
0
0
-1
-1
-1
-8
Example 6: Small networks as flows
Solve four small networks as Model 3: the transportation problem of Ex. 4, the shortest path problem of Ex. 1 in Lecture 4.1, the network of Fig. 8 with bounds on its arcs, and the multigraph of Lecture 4.1, Sec. 1, whose cheaper mode between two nodes is capacitated.
Example 6(a): Ex. 4’s transportation problem
Solve Ex. 4’s transportation problem as a minimum cost network flow, and compare the total cost with Ex. 4(b)’s.
In Ex. 4 (Fig. 5), the LP model was created directly from the interlevel matrix C, the matrix of arc costs from one level of nodes to the next, a row for each supply node and a column for each demand node. Since the transportation problem is a special case of the MCNF problem, C can be converted into an incidence matrix, A, and then solved as an MCNF instance.
# Code block 52: Ex. 4's tableau as an arc list, solved as a network flowC = [86109; 912137; 149165]m, n =size(C)i =repeat(1:m, n) # supply node of each arcj = m .+repeat(1:n, inner = m) # demand node of each arcc =vec(C) # cost of each arc, column by columns = [55, 50, 40, -45, -20, -30, -30]A =list2incid(i, j, m + n)TCᵒ, xᵒ =mcnf(A, c, s)@show TCᵒ; reshape(xᵒ, m, n)
1
vec(C) lists C column by column, so the arcs are listed the same way: supply nodes 1, 2, 3, 1, 2, 3, … and demand nodes 4, 4, 4, 5, 5, 5, …, numbered after the supply nodes. Listed this way, reshape(xᵒ, m, n) returns the flows in the tableau’s layout.
The minimum cost network flow returns the same total, 970, as the transportation model did in Ex. 4(b): two different formulations reach one answer, the Triangulate check of Lecture 1.2.
As a minimum cost network flow, Ex. 4’s tableau costs 970, the total of its transportation problem
Example 6(b): Shortest path
Determine the shortest path from node 1 to node 6 of Lecture 4.1’s six-node network, its arcs directed (Fig. 7 in Lecture 4.1, repeated here as Fig. 9 with what makes it a flow), as a minimum cost network flow, and then the shortest such path that uses the arc from node 2 to node 3.
The shortest path problem from Lecture 4.1, Ex. 1(c) can be formulated as an MCNF problem by adding a supply/demand of one to starting/terminal nodes and a net supply of 0 for the other (transshipment) nodes in the network:
# Code block 54: Ex. 6(b)'s network as a flow, drawn with Code block 50# The nodes where Lecture 4.1 draws them, and each cost on the side of# its arc where Lecture 4.1 puts the arc's lengthXY6 = [(0.0, 0.0), (1.0, 0.58), (1.0, -0.58), (2.0, 0.58), (2.0, -0.58), (3.0, 0.0)]fig =drawnet(XY6, i, j, string.(c), s; showzero =true, w =440, h =230, lim = ((-0.6, 3.6), (-0.95, 0.95)), side = [1, -1, -1, 1, 1, 1, -1, 1, -1])# The lower bound on arc (2,3), on the side of it clear of its costtext!(content(fig[1, 1]), 1.07, 0.24; text = L"l = 1", color = ink, fontsize =13, align = (:left, :center))fig
Figure 9: Lecture 4.1’s six-node network, its arcs directed, as a minimum cost network flow, each node labeled with its net supply. Node 1 (red) has a net supply of 1 and node 6 (green) a net supply of -1, so one unit flows from node 1 to node 6, and the other four nodes (black) are transshipment nodes, of zero net supply. Each arc is labeled with its cost per unit, the length Lecture 4.1 gives it, and l = 1 beside arc (2,3) is the lower bound on its flow, which forces the unit through it.
# Code block 55: the shortest path as a minimum cost network flowTCᵒwoBnds, xᵒwoBnds =mcnf(A, c, s)@show TCᵒwoBnds; xᵒwoBnds'
Note that although the solution gives the minimum cost flow and the sequence of edges, it would take additional processing to determine what the actual sequence of nodes is along the shortest path.
Additional constraints: The only reason to use the MCNF formulation for the shortest path problem would be if there needed to be additional constraints like arc bounds that need to be included in the model of the problem. For example, to find the shortest path in the directed network that always contains the arc from node 2 to node 3, solve as MCNF with a lower bound of arc (2,3) = 1 to force its use:
# Code block 56: a lower bound of 1 on arc (2,3)idx =findfirst(==((2, 3)), collect(zip(i, j))) # Arc (2,3)'s indexl =zeros(length(i))l[idx] =1u =fill(Inf, length(i))[u l]
Without bounds, the network flow sends its one unit along 1 → 2 → 4 → 5 → 6 at a total length of 14, the length Dijkstra’s procedure found for the directed network in Ex. 1(c) in Lecture 4.1: two procedures, one answer, the Triangulate check of Lecture 1.2. Forced through arc (2,3), the unit goes 1 → 2 → 3 and then on to node 6 at a length of 18, which Dijkstra’s procedure alone could not have been asked to find.
The shortest path from node 1 to node 6 has length 14; the shortest through arc (2,3), 18
Example 6(c): Arc bounds
Solve the network of Fig. 8, then solve it again with an upper bound of 3 on arcs (2,4) and (3,5), an upper bound of 9 on arc (2,5), and a lower bound of 1 on arc (4,5), and compare the two total costs.
Upper and lower arc bounds that differ from the default values of \infty and 0 can be added to the MCNF LP formulation.
Arc (4,5) is the only arc with a lower-bound >0, and arcs (2,4), (2,5), and (3,5) are the only arcs with upper-bounds <\infty. In order to provide flexibility with respect to inputting the arc information, the arc list is input as a comma-separated string, where any default upper and lower bounds are not input. This can then be used to create a DataFrame and then the missing default upper and lower bounds can be replaced with their default values.
The arc list, which so far has three columns, the start and end of each arc and its cost, takes the arc bounds as a fourth and a fifth column, the upper bound and then the lower bound. By default an upper bound is infinity and a lower bound is zero, and the columns are needed only for values other than the defaults: an upper bound less than infinity goes in the fourth column and a lower bound greater than zero in the fifth. The upper bound comes fourth and the lower bound last because the lower bound’s default of zero is more often what a problem wants than the upper bound’s of infinity, so a list of four columns leaves every lower bound at zero.
# Code block 58: the arc list, with its bounds, as a stringusingDataFrames, CSVarclist_str ="""i,j,c,u,l1,2,2,,1,3,3,,2,5,5,9,2,4,4,3,3,5,1,3,4,5,3,,1"""println(arclist_str)arclist_str
# Code block 59: the string read into a DataFramedf = CSV.read(IOBuffer(arclist_str), DataFrame)prt(string.(df))
1
IOBuffer makes the string read like a file, so CSV.read parses it into a DataFrame as it would a .csv file on disk.
2
Each empty field of the string is read as missing, which Code block 60 fills. prt cannot yet print a missing in a table of numbers, so string. turns every entry into text first; df itself keeps its numbers.
i j c u l
───────────────────────────
1 2 2 missing missing
1 3 3 missing missing
2 5 5 9 missing
2 4 4 3 missing
3 5 1 3 missing
4 5 3 missing 1
Code block 60 fills the empty bounds with their defaults using coalesce, used in Lecture 2.5, Sec. 7 for imputing missing data: coalesce.(x, v) replaces each missing in x with v and leaves every other value as it is, so an upper bound left empty becomes Inf and a lower bound left empty becomes 0.
# Code block 60: the default bounds filled indf.u =coalesce.(df.u, Inf) # Replace missing values in 'u' with infinitydf.l =coalesce.(df.l, 0) # Replace missing values in 'l' with zeroprt(df)
Fig. 10 draws the network with these bounds, each arc labeled with its cost and then with any bound that differs from its default:
Show the code that draws this figure
# Code block 61: the network with its bounds, drawn with Code block 50# An arc's label: its cost, then its upper and lower bound where either# differs from its default, an upper bound of Inf written as ∞arclab(c, u, l) = l >0 ? "$c, $(isinf(u) ? "∞":Int(u)), $l":isinf(u) ? "$c":"$c, $(Int(u))"drawnet(XY5, df.i, df.j, arclab.(df.c, df.u, df.l), [6, 2, 0, 0, -8]; showzero =true, rotate =true, dirs = [180, 90, 270, 90, 270], legend = (3.0, 1.0), w =520, h =300, lim = ((-0.6, 4.1), (-0.45, 2.3)))
Figure 10: The network of Fig. 8 with the bounds of Ex. 6(c), each node labeled with its net supply. Each arc is labeled c, u, l: its cost, its upper bound and its lower bound, a default left off the end, as the key at the right reads. Arc (4,5)’s lower bound of 1 is the only lower bound above zero, and arcs (2,4), (2,5) and (3,5) carry the only upper bounds below infinity.
# Code block 62: the network without its bounds, and with thems = [6, 2, 0, 0, -8]A =list2incid(df.i, df.j, length(s))TCᵒwoBnds, xᵒwoBnds =mcnf(A, df.c, s)TCᵒ, xᵒ =mcnf(A, df.c, s; u = df.u, l = df.l)@show TCᵒ@show TCᵒwoBndsdf.xᵒ, df.xᵒwoBnds = xᵒ, xᵒwoBndsprt(df)
Without bounds, all six of node 1’s units take the path 1 → 3 → 5 at 4 per unit, and node 2’s two units take arc (2,5). The upper bound of 3 on arc (3,5) halves the cheap path’s flow, so node 1’s other three units go by way of node 2; the lower bound of 1 on arc (4,5) forces one unit onto the costlier 2 → 4 → 5; and the total rises from 34 to 45. Bounds can only add constraints, so the total with them can be no less than without them: the Bounds check of Lecture 1.2, passed.
With its bounds, the five-node network’s least total cost is 45, against 34 without them
Example 6(d): A capacitated mode
In the multigraph of Fig. 4 in Lecture 4.1, nodes 1, 2 and 3 supply 30, 10 and 20 units, and nodes 4 and 5 demand 40 and 20. The cheaper of the two modes between nodes 1 and 4, at a cost of 6 per unit, can carry at most 20 units; the other, at 18, is uncapacitated. Determine the least-cost flows, with and without the capacity.
Of the two arcs between nodes 1 and 4, the 18, at the higher cost, might represent a more costly mode of transport, and the 6 a less costly one with limited capacity, so that only so much flow can use the cheaper mode, and beyond that quantity the higher-cost arc is needed. Representing that takes upper and lower bounds on the arcs. So there are contexts where a multigraph is wanted, and when it is, the arc list has to be used, since the adjacency matrix cannot represent it. That is one of the reasons for having several representations of a network: some cannot represent everything, and in other cases one is simply more convenient than another, as the adjacency matrix often is, even though it cannot represent everything.
Fig. 11 draws the multigraph of Lecture 4.1 with the net supplies and the capacity of this example:
Show the code that draws this figure
# Code block 63: Ex. 6(d)'s multigraph, drawn with Code block 50# The nodes where Lecture 4.1 draws them; each two-way pair of arcs 1-4# drawn as one arc with a head at each endXYm = [(0.0, 2.0), (0.0, 1.0), (0.0, 0.0), (2.4, 1.6), (2.4, 0.4)]drawnet(XYm, [1, 1, 1, 2, 2, 3, 3], [4, 4, 5, 4, 5, 5, 5], ["6, 20", "18", "10", "0", "13", "16", "3"], [30, 10, 20, -40, -20]; off = [0.07, -0.07, 0, 0, 0, 0.07, -0.07], two = [true, true, false, false, false, false, false], at = [0.5, 0.5, 0.2, 0.3, 0.3, 0.5, 0.5], side = [1, -1, -1, 1, -1, 1, -1], w =400, h =330, lim = ((-0.7, 3.1), (-0.4, 2.4)))
Figure 11: Lecture 4.1’s multigraph as a minimum cost network flow, each node labeled with its net supply: nodes 1, 2 and 3 (red) supply 30, 10 and 20 units, and nodes 4 and 5 (green) demand 40 and 20. Each arc is labeled with its cost per unit. The two two-way arcs between nodes 1 and 4 are the two modes, and the cheaper’s label, 6, 20, gives its cost per unit and the 20 units it can carry from node 1 to node 4.
The arc list is the one Lecture 4.1 draws as Fig. 4 in Lecture 4.1, as it stands, the parallel arcs 1 → 4 at 6 and at 18 each a row of their own, and the arc of cost 0 from node 2 to node 4 kept as it is. Uncapacitated, all 30 of node 1’s units take the cheaper mode, for a total of 240. Capacitated, the cheaper mode carries its 20 and the costlier one the other 10, for 360: the costlier mode is used only for the flow the cheaper one cannot carry. An adjacency matrix could hold only one of the two arcs from 1 to 4, and so could not pose this problem at all.
With the cheaper mode capacitated at 20 units, the least total cost is 360, against 240 uncapacitated
Example 7: Earthquake emergency response
Determine the least-cost flows, as Model 3, that deliver the 700 pallets the Field Hospital needs through the relief network of Table 6; if the network as it stands cannot deliver them, determine which arc’s capacity to raise, and the least total cost with it.
A major earthquake has struck a region, and emergency medical supplies must be transported to a central field hospital. Two warehouses have supplies, but the routes are compromised, with varying costs, travel risks, and transportation capacities. The goal is to meet the hospital’s demand at the lowest possible cost. The relief effort’s logistics network is composed of six key locations, as shown in Table 5. The supply originates from two source nodes: the Regional Warehouse (RW), which holds 500 pallets of general supplies, and a Specialty Pharmacy (SP) with 200 pallets of critical medicines. From these points, supplies travel through three transshipment nodes: a Local Airfield (LA) for consolidation, a Central Distribution Hub (CH) for sorting, and a Staging Area (SA) for final preparation. The ultimate goal is to deliver the combined 700 pallets to the single demand node, the Field Hospital (FH), to make it fully operational.
Table 5: The relief network’s nodes, each with its net supply: positive at a supply node, zero at a transshipment node and negative at the demand node.
Node
Name
Type
Net supply
1
Regional Warehouse (RW)
Supply
500
2
Specialty Pharmacy (SP)
Supply
200
3
Local Airfield (LA)
Transshipment
0
4
Central Distribution Hub (CH)
Transshipment
0
5
Staging Area (SA)
Transshipment
0
6
Field Hospital (FH)
Demand
-700
Table 6: The relief network’s arcs, each with its cost per pallet and the least and the most it can carry.
From Node
To Node
Cost ($/pallet)
Lower Bound (pallets)
Upper Bound (pallets)
Notes
1 (RW)
3 (LA)
4
0
400
Main road, limited trucks
1 (RW)
4 (CH)
6
0
300
Secondary road, some damage
2 (SP)
4 (CH)
3
0
200
Direct, efficient route
3 (LA)
5 (SA)
10
0
150
Expensive helicopter transport
4 (CH)
5 (SA)
5
100
600
Main convoy, needs scarce escort detail
5 (SA)
6 (FH)
2
0
800
Secure final-mile shuttle
The arcs are entered as Ex. 6(c)’s were, an arc list in a string, the empty lower bounds read as missing and set to zero, and the net supplies are those of Table 5:
# Code block 65: the relief network's arcs, with their bounds, as a flowarclist_str ="""i,j,c,u,l1,3,4,400,1,4,6,300,2,4,3,200,3,5,10,150,4,5,5,600,1005,6,2,800,"""df = CSV.read(IOBuffer(arclist_str), DataFrame)df.l =coalesce.(df.l, 0)s = [500, 200, 0, 0, 0, -700]A =list2incid(df.i, df.j, length(s))prt(df)
Since this problem has restrictive upper bounds, there is a significant probability that it may not be feasible to satisfy the required demand. As a result, the call to mcnf, the implementation of Model 3, is wrapped so that its error, which names the solver’s status, is printed rather than stopping the script.
# Code block 66: the network as it standstrymcnf(A, df.c, s; u = df.u, l = df.l)catcheshowerror(stdout, e)end
INFEASIBLE
Since the model is infeasible, likely due to the upper bounds on flows being too restrictive, one way to look at this is to examine the network structure, where the label on each arc gives its cost and its bounds (Fig. 12).
Figure 12: The relief network of Ex. 7, each node labeled with its name and its net supply: the Regional Warehouse (RW) and the Specialty Pharmacy (SP) supply 500 and 200 pallets, the Field Hospital (FH) needs 700, and the airfield (LA), the hub (CH) and the staging area (SA) are transshipment nodes. Each arc is labeled c, u, l: its cost per pallet, its upper bound and its lower bound, a default left off the end, as the key at the top right reads.
In looking at the graph, at most 500 units can flow into Node 4 and then on to Node 5. Since only 150 units can travel from 3 to 5, the total of 650 is less than the required demand of 700 at Node 6. One way to fix this is to “somehow try to” increase the capacity between 3 and 5 to 200.
# Code block 68: the helicopter's capacity, arc (3,5), raised to 200df2 =deepcopy(df) # deepcopy used so that `df2` changes don't change `df`df2.u[4] =200prt(df2)
Feasible now. But this was achieved by increasing the expensive helicopter capacity. Another possibility is to “somehow try to” increase the capacity between 1 and 4 to 350. At a cost of $6 per unit, increasing this by 50 is cheaper than increasing the expensive $10 per unit helicopter capacity by 50.
Raising road (1,4) instead delivers all 700 pallets for 8,950, 150 less than raising the helicopter’s: the 50 pallets the path 1 → 3 → 5 carried at 4 + 10 = 14 a pallet go by road, 1 → 4 → 5, at 6 + 5 = 11, a saving of 3 × 50 = 150, the Triangulate check of Lecture 1.2. The count before Code block 68, at most 650 pallets into node 5 against the 700 needed, is the Bounds check: it explains the infeasibility without the linear program.
As the network stands, no flow delivers the 700 pallets; raising road (1,4)’s capacity to 350 delivers them at a least total cost of 8,950, against 9,100 by raising the helicopter’s to 200
NC State University, “Work to Begin on Traffic Gates,” news.ncsu.edu/2012/12/gates-2, December 6, 2012. OpenStreetMap ways 236483130 and 1474912879, tagged vehicle=no, bus=yes and vehicle:conditional=yes @ (17:00-09:00), with Wolfline routes 30, 51 and 52 mapped over them, accessed 2026-10-14.↩︎
OSRM car profile, profiles/car.lua and profiles/lib/way_handlers.lua in github.com/Project-OSRM/osrm-backend, accessed 2026-10-14: a way’s maxspeed times speed_reduction = 0.8, and turn_penalty = 7.5 (s) scaled by the angle of the turn; vehicle is in its access_tags_hierarchy, and no in its access_tag_blacklist.↩︎
