Matlog: Logistics Engineering Matlab Toolbox
Version 10 12-Jun-2008
| ala | - Alternate location-allocation procedure |
| alaplot | - Plot ALA iterations |
| minisumloc | - Continuous minisum facility location |
| randX | - Random generation of new points |
| sumcost | - Determine total cost of new facilities |
| ufl | - Discrete uncapacitated facility location |
| aggshmt | - Aggregate multiple shipments |
| maxpayld | - Determine maximum payload |
| plotshmt | - Plot shipments |
| rateLTL | - Determine estimated LTL rate |
| tlcLTL | - Total logistics cost for LTL shipments |
| tlcmin | - Minimum total logistics cost comparing TL and LTL |
| tlcTL | - Total logistics cost for TL shipments |
| add1sh2rte | - Add independent shipments to routes |
| checkrte | - Check if valid route vector |
| loccost | - Calculate location sequence cost |
| mincostinsert | - Min cost insertion of route i into route j |
| pairwisesavings | - Calculate pairwise savings |
| rte2idx | - Convert route to route index vector |
| rte2ld | - Convert route vector to load cell array |
| rte2loc | - Convert route to location vector |
| rtenorm | - Normalize route vector |
| rteTC | - Total cost of route |
| savings | - Savings procedure for route construction |
| sh2rte | - Create routes for shipments not in existing routes |
| twoopt | - 2-optimal exchange procedure for route improvement |
| combineloc | - Combine non-overlapping location sequences |
| locTC | - Calculate total cost of location sequences |
| tspchinsert | - Convex hull insertion algorithm for TSP construction |
| tspnneighbor | - Nearest neighbor algorithm for TSP construction |
| tspspfillcur | - Spacefilling curve algorithm for TSP construction |
| tsp2opt | - 2-optimal exchange procedure for TSP improvement |
| vrpcrossover | - Crossover procedure for VRP improvement |
| vrpexchange | - Two-vertex exchange procedure for VRP improvement |
| vrpinsert | - Insertion algorithm for VRP construction |
| vrpsavings | - Clark-Wright savings algorithm for VRP construction |
| vrpsweep | - Gillett-Miller sweep algorithm for VRP construction |
| vrptwout | - Generate output VRP with time windows |
| addconnector | - Add connector from new location to transport. network |
| adj2incid | - Convert weighted adjacency matrix to incidence matric |
| adj2lev | - Convert weighted adjacency to weighted interlevel matrix |
| adj2list | - Convert weighted adjacency matrix to arc list |
| dijk | - Shortest paths using Dijkstra algorithm |
| lev2adj | - Convert weighted interlevel to weighted adjacency matrix |
| incid2list | - Convert node-arc incidence matrix to arc list |
| list2adj | - Convert arc list to weighted adjacency matrix |
| list2incid | - Convert arc list to incidence matrix |
| mcnf | - Minimum cost network flow problem |
| minspan | - Minimum weight spanning tree using Kruskal algorithm |
| pred2path | - Convert predecessor indices to shortest path |
| loc2listidx | - Find indices of loc seq segments in arc list |
| subgraph | - Create subgraph from graph |
| thin | - Thin degree-two nodes from graph |
| trans | - Transportation and assignment problems |
| tri2list | - Convert triangle indices to arc list representation |
| trineighbors | - Find neighbors of a triangle |
| city2lonlat | - Select longitude and latitude of city |
| invproj | - Inverse Mercator projection |
| lonlat2city | - Determine nearest city given a location |
| makemap | - Create projection plot of World or US |
| normlonlat | - Convert longitude and latitude to normal form |
| proj | - Mercator projection |
| zip2lonlat | - Determine longitude and latitude of ZIP code |
| binaryorder | - Binary ordering of machine-part matrix |
| loc2W | - Converts routings to weight matrix |
| sdpi | - Steepest descent pairwise interchange heuristic for QAP |
| mheselect | - Material handling equipment selection MILP model |
| dists | - Metric distances between vectors of points |
| gantt | - Gantt chart |
| lprsm | - Linear programming using the revised simplex method |
| milp | - Mixed-integer linear programming branch and bound |
| pplot | - Projection plot |
| arcang | - Arc angles (in degrees) from xy to XY |
| argmax | - Indices of maximum component |
| argmin | - Indices of minimum component |
| argsort | - Indices of sort in ascending order |
| boundrect | - Bounding rectangle of XY points |
| cell2padmat | - Convert cell array of vectors to NaN-padded matrix |
| deletelntag | - Delete last tagged line objects from plot |
| editcellstr | - Edit or create a cell array of strings |
| file2cellstr | - Convert file to cell array |
| findXY | - Find XY points in rectangle |
| fixseq | - Fixed sequence proportional to percentage of demand |
| getfile | - Download file from Internet and save local copy |
| grect | - Get rectangular region in current axes |
| idx2is | - Convert index vector to n-element logical vector |
| inputevaldlg | - Input dialog with evaluation of values |
| invperm | - Inverse permutation |
| isinrect | - Are XY points in rectangle |
| isint | - True for integer elements (within tolerance) |
| is0 | - True for zero elements (within tolerance) |
| loaddatafile | - Common code for data file functions |
| matlogupdate | - Automatically download updates to Matlog files |
| mat2vec | - Convert columns of matrix to vectors |
| mdisp | - Display matrix |
| num2cellstr | - Create cell array of strings from numerical vector |
| optstructchk | - Validate 'opt' structure |
| padmat2cell | - Convert rows of NaN-padded matrix to cell array |
| pauseplot | - Drawnow and then pause after plotting |
| projtick | - Project tick marks on projected axes |
| randreorder | - Random re-ordering of an array |
| sdisp | - Display structure array with all scalar field values |
| subsetstruct | - Extract subset of each field of a structure |
| sfcpos | - Compute position of points along 2-D spacefilling curve |
| sub | - Create subarray |
| varnames | - Convert input to cell arrays of variable name and value |
| vec2struct | - Add each element in vector to field in structure array |
| vdisp | - Display vectors |
| wtbinselect | - Weighted binary selection |
| wtrandperm | - Weighted random permutation |
| wtrouselect | - Weighted roulette selection |
| mapdata | - Data for maps of the World, US, and North Carolina |
| nccity | - North Carolina cities data |
| uscity | - US cities data |
| uscity10k | - US cities with populations of at least 10,000 data |
| uscity50k | - US cities with populations of at least 50,000 data |
| usrdlink | - US highway network links data |
| usrdnode | - US highway network nodes data |
| uszip3 | - US 3-digit ZIP code data |
| uszip5 | - US 5-digit ZIP code data |
| vrpnc1 | - Christofiles' VRP problem 1 data |
| vrpsolrc101 | - Soloman's VRP with Time Windows problem RC101 data |
Add connector from new location to transportation network
[IJC11,IJC12,IJC22] = addconnector(XY1,XY2,IJC2,p,cir,thresh)
XY1 = m1 x 2 matrix of m1 new location longitude-latitude pairs
(in decimal degrees)
XY2 = m2 x 2 matrix of longitude-latitude pairs of original network
nodes (in decimal degrees)
IJC2 = n2 x 3 matrix arc list of original network
p = distance parameter for DISTS, should be in same units as the
distances in IJC2
= 'mi', default, great circle distance in statute miles
cir = circuity factor (between 1 and 2) for connectors that represents
the expected increase in arc distance compared to DISTS distance
= 1.50, default
thresh = threshold (between 0 and 1) to only add arc to closest node in
original network (if ratio of distance from new location to
closest node and to second closest node is less than threshold)
= 0.10, default
= 0 => don't consider threshold
= 1 => only add arc to closest node
IJC11 = 3-column matrix arc list of connectors between new locations
IJC12 = 3-column matrix arc list of connectors from new location to
original nodes
IJC22 = 3-column matrix arc list of modified original network, where arc
[i j] becomes [i+m1 j+m1]
Examples:
% 3 new nodes connected to 4 original nodes (Euclidean arc distances)
XY1 = [2 1; 5 -1];
XY2 = [0 0; 4 3; 4 -3; 8 0];
IJC2 = [1 2 5; 1 3 5; 2 4 5; 3 4 5];
[IJC11,IJC12,IJC22] = addconnector(XY1,XY2,IJC2,2)
pplot(IJC11,XY1,'r-')
pplot(IJC12,[XY1; XY2],'g-')
pplot(IJC22,[XY1; XY2],'b-')
% Connect cites around Raleigh, NC to road network
xy1xy2 = [-79 35.5; -78 36];
[XY,IJD] = subgraph(usrdnode('XY'),isinrect(usrdnode('XY'),xy1xy2),...
usrdlink('IJD'));
[Name,XYcity] = uscity10k('Name','XY',isinrect(uscity10k('XY'),xy1xy2));
[IJD11,IJD12,IJD22] = addconnector(XYcity,XY,IJD);
makemap(XY)
pplot(IJD11,XYcity,'r-')
pplot(IJD12,[XYcity; XY],'g-')
pplot(IJD22,[XYcity; XY],'b-')
pplot(XYcity,'k.')
pplot(XYcity,Name)
Note: (1) Arc between a new location and the node in the original network
is added to the combined network when the great circle distance times
circuity factor is less than shortest distance between the locations in
remainder of network. Nodes considered are those defining the location's
triangle (using DELAUNAY) and, if not a triangle node, the closest node
to the location (using DISTS).
(2) Single connector arc from new location to closest node in original
network if distance < thresh x distance to second closest node;
connectors to other new locations added only their distance x thresh <=
distance of shortest connector.
(3) Arc between two new locations in the same or adjacent triangles is
added to combined network when the great circle distance times circuity
factor is less than shortest distance between the locations in remainder
of network.
(4) Arc between new locations outside the convex hull of XY2 added to
any new location that is one of the three closest nodes to the location.
Node-node weighted adjacency matrix to node-arc incidence matrix
[I,c] = adj2incid(A)
A = m x m node-node weighted adjacency matrix of arc lengths
I = m x n node-arc incidence matrix
c = n-element vector of arc weights (nonzero elements of A)
Note: A(i,j) = 0 => Arc (i,j) does not exist
A(i,j) = NaN => Arc (i,j) exists with 0 weight
All A(i,j) = A(j,i) => A is symmetric => all undirected arcs in I
Wrapper for [I,c] = LIST2INCID(ADJ2LIST(A))
I is column major relative to A for speed purposes
Examples:
A = [0 1 0
0 0 0
3 2 0]
[I,c] = adj2incid(A) % I = -1 1 0 c = 3
0 -1 -1 1
1 0 1 2
A = [0 1 0 % Arc (3,1) has 0 weight
0 0 0
NaN 2 0]
[I,c] = adj2incid(A) % I = -1 1 0 c = 0
0 -1 -1 1
1 0 1 2
A = [0 1 3 % A is symmetric
1 0 2
3 2 0]
[I,c] = adj2incid(A) % I = 1 1 0 c = 1
1 0 1 3
0 1 1 2
See also LIST2INCID, LIST2ADJ, and ADJ2LIST
Weighted adjacency to weighted interlevel matrix representation
C = adj2lev(A,m)
[C12,C23,...] = adj2lev(A,m)
A = SUM(m) x SUM(m) node-node wt. adjacency matrix of arc lengths
= (SUM(m) x t) x (SUM(m) x t) matrix, t > 1, if multiple periods
m = [m1 m2 m3 ...], vector of number of nodes at each level, where
C = {C12,C23,...}, cell array of weighted interlevel matrices
Cij = mi x mj matrix from level i to level j
= mi x mj x t matrix, t > 1, if multiple periods
mi = number of nodes in level i
t = number of periods
Examples:
A = [0 3 4
0 0 0
0 0 0]
C = adj2lev(A,[1 2]) % (1 period C = 3 4
% 2 levels)
A = blkdiag(A,2*A)
C = adj2lev(A,[1 2]) % (2 periods C(:,:,1) = 3 4
% 2 levels) C(:,:,2) = 6 8
A = [0 3 4 0
0 0 0 5
0 0 0 6
0 0 0 0
[C12,C23] = adj2lev(A,[1 2 1]) % (1 period C12 = 3 4
% 3 levels) C23 = 5
6
See also LEV2ADJ
Node-node weighted adjacency matrix to arc list representation
IJC = adj2list(A)
[i,j,c] = adj2list(A)
A = m x m node-node weighted adjacency matrix of arc lengths
IJC = n x 2-3 matrix arc list [i j c], where
i = n-element vector of arc tails nodes
j = n-element vector of arc head nodes
c = n-element vector of arc weights
Note: All A(i,j) = A(j,i) => [i -j c] (symmetric A)
A(i,j) = 0 => Arc (i,j) does not exist
A(i,j) = NaN => Arc (i,j) exists with 0 weight
Wrapper for [i,j,c] = FIND(C); c(ISNAN(c)) = 0)
Examples:
A = [0 1 0
0 0 0
3 2 0]
IJC = adj2list(A) % IJC = 3 1 3
% 1 2 1
% 3 2 2
A = [0 1 0 % Arc (3,1) has 0 weight
0 0 0
NaN 2 0]
IJC = adj2list(A) % IJC = 3 1 0
% 1 2 1
% 3 2 2
A = [0 1 3 % A is symmetric
1 0 2
3 2 0]
IJC = adj2list(A) % IJC = 1 -2 1
% 1 -3 3
% 2 -3 2
See also LIST2INCID, LIST2ADJ, and ADJ2INCID
Aggregate multiple shipments
shagg = aggshmt(sh)
sh = structure array with fields:
.f = annual demand (tons/yr)
.q = single shipment weight (tons)
.s = shipment density (lb/ft^3)
.v = shipment value ($/ton)
shagg = single aggregate structure
Note: 1. Aggregation based on "f" values unless "f" field is missing or
empty, in which case it is based on "q" values.
2. Aggregate "s" and "v" determined only if respective fields are
present
3. All other fields are passed to aggregate shipment based on
their values in the first input shipment structure sh(1)
Alternate location-allocation procedure
[X,TC,W] = ala(X0,w,P,p) % Use default allocate and locate
= ala(X0,alloc_h,P,p) % Use MINISUMLOC(P,W,p) to locate NFs
= ala(X0,w,P,p,loc_h) % Use MIN(DISTS(X,P,p)) to allocate NFs
= ala(X0,alloc_h,loc_h) % User-defined allocate and locate
X0 = n-row matrix of n new facility (NF) starting locations
(use X0 = randX(P,n) to generate n random locations)
P = m-row matrix of m existing facility (EF) locations
w = m-element vector of weights
p = parameter for DISTS (default 'mi')
alloc_h = handle to anonymous function to allocate NFs
= @(X) myalloc(...,X,...), where MYALLOC must return W and TC
[W,TC] = myalloc(...,X,...)
loc_h = handle to anonymous function to locate NFs
= @(W) myloc(...,W,...), where MYLOC must return X
X = myloc(...,W,...)
X = n-row matrix of n NF locations
TC = total cost of allocation
W = n x m matrix of weights, where W(i,j) represents the weight
allocated to/from NF(i) from EF(j)
Run PPLOT(P,'r.') or MAKEMAP(P) before running ALA, if you want to plot
intermediate results, and use "Set Pauseplot Time" on the "Matlog"
menu in the figure to to adjust the frequency.
Run ALAPLOT(X,W,P) to plot only the final results.
%Example 1: Locate three NFs to serve customers in North Carolina cities
[P,w] = nccity('XY','Pop'); p = 'mi';
makemap(P), pplot(P,'r.')
ala(randX(P,3),w,P,p)
%Example 2: Take the best of 'nruns' runs of ALA
nruns = 5; TC = Inf;
for i=1:nruns, [X1,TC1,W1] = ala(randX(P,3),w,P,p); ...
fprintf('%d %e\n',i,TC1); if TC1 < TC, TC=TC1; X=X1; W=W1; end, end
%Example 3: Constrained NFs (each NF can handle a third of total demand)
alloc_h = @(XY) trans(dists(XY,P,p),ones(1,3)*sum(w)/3,w);
[XYc,TCc,Wc] = ala(randX(P,3),alloc_h,P,p); % Constrained NFs
TCc, pctdemC = full(sum(Wc,2))/sum(w) % Percentage of total demand
%Example 4: Compare constrained to unconstrained NFs
[XYu,TCu,Wu] = ala(XYc,w,P,p);
TCu, pctdemU = full(sum(Wu,2))/sum(w)
Plot ALA iterations
alaplot(X,W,P) Example: Plot final results from ALA P = [0 0;2 0;2 3], w = [1 2 1] [X,TC,W] = ala(2,w,P,2) pplot alaplot(X,W,P)
Arc angles (in degrees) from xy to XY
ang = arcang(xy,XY)
Indices of maximum component
[i,j,y] = argmax(X), where y = X(i,j) = MAX(MAX(X))
i = argmax(X), where [y,i] = MAX(X)
i = argmax(X,DIM), where [y,i] = MAX(X,[],dim)
Indices of minimum component
[i,j,y] = argmin(X), where y = X(i,j) = MIN(MIN(X))
i = argmin(X), where [y,i] = MIN(X)
i = argmin(X,DIM), where [y,i] = MIN(X,[],dim)
Indices of sort in ascending order
[i,j,y] = argsort(X), where [i,j] = IND2SUB(SIZE(X),ARGSORT(X(:)))
I = argsort(X), where [Y,II] = SORT(X);
I = argsort(X,DIM), where [Y,II] = SORT(X,DIM);
Binary ordering of machine-part matrix
[A,idxi,idxj] = binaryorder(Ain)
Ain = m x n 0-1 machine-part matrix
A = machine-part matrix in binary order, where A = Ain(idxi,idxj)
idxi = m-element index vector of machine orderings
idxj = n-element index vector of part orderings
Bounding rectangle of XY points
xy1xy2 = boundrect(XY,expand)
expand = nonnegative expansion factor for bounding rectangle, where
expansion of "expand" times max X,Y extent is added to all
sides of rectangle
= 0, default
= 0.1 => 10% expansion
xy1xy2 = bounding rectangle
= [min(X) min(Y); max(X) max(Y)]
Convert cell array of vectors to NaN-padded matrix
M = cell2padmat(C,align)
align = 1, left alignment (default)
= 2, right alignment
Check if valid route vector
checkrte(rte) % Check route and throw error if not correct
checkrte(rte,sh) % Also check that rte values in sh range
% and that sh has all scalar field values
checkrte([],sh) % Just check sh has all scalar field values
checkrte(rte,[],isvec) % Require single vector rte
checkrte([],[],[],iscklim) % Limit checking
rte = route vector
= m-element cell array of m route vectors
sh = structure array of shipments
isvec = rte must be only a single vector, not a cell array
iscklim = limit checking to inputs from calling function in base workspace
in order to eliminate duplicate checking
= true, default
Select longitude and latitude of city
XY = city2lonlat(name,st)
= city2lonlat(name,st,city)
s = city2lonlat(...), select city in dialog if multiple 'name' or
'st' matches and return stucture 's' of city
data if multiple cities selected in dialog
name = city name string or cell array of multiple strings
= [], all names
st = 2-char state abbreviation string or, if multiple values, cell
array of strings
= [], all states
city = (optional) stucture with fields:
.Name = m-element cell array of m city name strings
.ST = m-element cell array of m 2-char state abbreviations
.XY = m x 2 matrix of city lon-lat (in decimal deg)
.Pop = m-element vector of city population estimates (2000)
= USCITY10K, default (which contains all cities in US with
population of at least 10,000)
XY = 1 x 2 row vector of city 'name','st' lon-lat (in decimal,
degrees) if single state and multiple city names, or equal
number of citynames and states
s = stucture of city data, if multiple cities selected in dialog
Examples:
xy = city2lonlat('Raleigh','NC') % xy = -78.6446 35.8188
XY = city2lonlat({'Raleigh','Gainesville'},{'NC','FL'})
% XY = -78.6446 35.8188
% -82.3361 29.6652
XY = city2lonlat([],'NC') % Dialog to select cities in North
% Carolina with > 10,000 population
XY = city2lonlat([],'NC',uscity50k) % Dialog to select cities in North
% Carolina with > 50,000 population
See also LONLAT2CITY
Combine non-overlapping loc seqs
[cmbloc,st] = combineloc(out)
out = m-element struct array of timing output from LOCTC
cmbloc = combine non-overlapping loc seqs, where
cmbloc{i} = [j k] represents the i-th combined loc seq consisting
of seq j followed by seq k
Uses greedy combination heuristic
Example:
[TC,out] = locTC(loc,C,cap,twin)
[cmbloc,st] = combineloc(out);
[TC,out] = locTC(cmbloc,C,cap,twin); % Re-calc "out"
Delete last tagged line objects from plot
deletelntag()
Executes commands:
s = get(gca,'UserData');
delete(findobj(s.LastLineTag))
projtick
Shortest paths from nodes 's' to nodes 't' using Dijkstra algorithm
[D,P] = dijk(A,s,t)
A = n x n node-node weighted adjacency matrix of arc lengths
(Note: A(i,j) = 0 => Arc (i,j) does not exist;
A(i,j) = NaN => Arc (i,j) exists with 0 weight)
s = FROM node indices
= [] (default), paths from all nodes
t = TO node indices
= [] (default), paths to all nodes
D = |s| x |t| matrix of shortest path distances from 's' to 't'
= [D(i,j)], where D(i,j) = distance from node 'i' to node 'j'
= Inf, if no path from 'i' to 'j'
P = |s| x n matrix of predecessor indices, where P(i,j) is the
index of the predecessor to node 'j' on the path from 's(i)' to
'j',where P(i,i) = 0 and P(i,j) = NaN is 'j' not on path to
's(i)' (use PRED2PATH to convert P to paths)
= path from 's' to 't', if |s| = |t| = 1
Example:
A = [0 1 3 0
0 0 0 2
0 0 0 4
0 0 0 0]
[d,p] = dijk(A,1,4) % (Single path) d = 3
% p = 1 2 4
[D,P] = dijk(A) % (All paths) D = 0 1 3 3
% Inf 0 Inf 2
% Inf Inf 0 4
% Inf Inf Inf 0
% P = 0 1 1 2
% NaN 0 NaN 2
% NaN NaN 0 3
% NaN NaN NaN 0
p = pred2path(P,1,4) % p = 1 2 4
(If A is a triangular matrix, then computationally intensive node
selection step not needed since graph is acyclic (triangularity is a
sufficient, but not a necessary, condition for a graph to be acyclic)
and A can have non-negative elements)
(If |s| >> |t|, then DIJK is faster if DIJK(A',t,s) used, where D is now
transposed and P now represents successor indices)
(Based on Fig. 4.6 in Ahuja, Magnanti, and Orlin, Network Flows,
Prentice-Hall, 1993, p. 109.)
Metric distances between vectors of points
D = dists(X1,X2,p,e)
X1 = n x d matrix of n d-dimensional points
X2 = m x d matrix of m d-dimensional points
D = n x m matrix of distances
p = 2, Euclidean (default): D(i,j) = sqrt(sum((X1(i,:) - X2(j,:))^2))
= 1, rectilinear: D(i,j) = sum(abs(X1(i,:) - X2(j,:))
= Inf, Chebychev dist: D(i,j) = max(abs(X1(i,:) - X2(j,:))
= (1 Inf), lp norm: D(i,j) = sum(abs(X1(i,:) - X2(j,:))^p)^(1/p)
= 'rad', great circle distance in radians of a sphere
(where X1 and X2 are decimal degree longitudes and latitudes)
= 'mi' or 'sm', great circle distance in statute miles on the earth
= 'km', great circle distance in kilometers on the earth
e = epsilon for hyperboloid approximation gradient estimation
= 0 (default); no error checking if any non-empty 'e' input
~= 0 => general lp used for rect., Cheb., and p outside [1,2]
Examples:
x1 = [1 1], x2 = [2 3]
d = dists(x1,x2,1) % d = 3
X2 = [0 0;2 0;2 3]
d = dists(x1,X2,1) % d = 2 2 3
D = dists(X2,X2,1) % D = 0 2 5
% 2 0 3
% 5 3 0
d=dists(city2lonlat('Raleigh','NC'),city2lonlat('Gainesville','FL'),'mi')
% d = 475.8047
Great circle distances are calculated using the Haversine Formula (R.W.
Sinnott, "Virtues of the Haversine", Sky and Telescope, vol 68, no 2,
1984, p 159, reported in "http://www.census.gov/cgi-bin/geo/gisfaq?Q5.1")
Edit or create a cell array of strings
c = editcellstr(c) % Edit an existing cell array of strings "c"
= editcellstr(n) % Create an "n"-element array
= editcellstr % Create up to 10-element array, with empty trailing
% elements removed
Examples:
c = {'Boston','New York','Chicago'};
c = editcellstr(c)
c = editcellstr(3)
c = editcellstr
Convert file to cell array
c = file2cellstr(fname)
fname = name of file
= [], open dialog to get file name
c = cell array of strings, where each string is a line of the file
Find XY points inside rectangle
idx = findXY(XY,xy1xy2)
= findXY(XY) Graphical selection of rectangle
xy1xy2 = rectangle
= [min(X) min(Y); max(X) max(Y)]
Graphical selection of rectangle:
(1) Press or select point outside axes to toggle zoom
(2) Press key to erase last points selected
(3) Press to stop selecting points
(4) Use 'Erase Iter Plot' on Matlog figure menu to erase green circles
around selected points (or 'delete(findobj('Tag','iterplot'))')
Note: Wrapper for idx = FIND(ISINRECT(XY,xy1xy2))
Fixed sequence proportional to percentage of demand
s = fixseq(a,tol)
a = n-element vector of demands
tol = tolerance of rational approximation to demand percentage
= 1e-2, default
s = sequence of indices i = 1 to n, where the frequency of i's in "s"
is proportional to the percentage of a(i) in the total demand,
i.e., |sum(s == i)/length(s) - a(i)/sum(a)| < tol
Example:
If a = [0.1193 0.0504] then
a/sum(a) is [0.7030 0.2970]
fixseq(a) is [1 1 2 1 1 1 2 1 1 2]
fixseq(a,1e-1) is [1 1 2]
Gantt chart
h = gantt(Bars) % Plot activity bars
= gantt(Bars,'PropName',PropValue,...) % Bar properties
= gantt(Bars,Labels) % Label each bar
= gantt(Bars,Labels,'PropName',PropValue,...) % Label properties
Bars = m-element cell array, where Bars{i} = n x 2 matrix of n activity
bars for resource i and Bars{i}(j,1) and Bars{i}(j,2) are the
beginning and ending times for activity j
= m x 2 matrix, if each resource has only one activity
Labels = m-element cell array of bar labels
h = m-element cell array of handles to bar rectangles or labels
Example:
Bars = {[1 4]; [2 4; 5 6]; [1 2; 2 3]}
Labels = {'bar1';{'bar2','bar3'};{'bar4','bar5'}}
hbars = gantt(Bars)
hlabels = gantt(Bars,Labels)
Download file from Internet and save local copy
XFlg = getfile(baseurl,fnin,fnout)
baseurl = base URL address
= 'http://www.ie.ncsu.edu/kay/matlog/', default
fnin = name of file to read from from web
fnout = name of file to write to local directory
= 'fnin', default
XFlg = 1, able to get file
= 0, not able
Get rectangular region in current axes
xy1xy2 = grect
xy1xy2 = [min(X) min(Y); max(X) max(Y)]
In order to support the features of FINDXY:
if key pressed instead of mouse button, then decimal current character
is returned;
if first point selected is outside axes, then empty return and zoom
state is toggled;
if second point selected is outside axes, then empty return
Convert index vector to n-element logical vector
is = idx2is(idx,n)
idx = index vector
n = length of logical vector
is = logical vector, is(idx) = 1
Inverse of FIND(is) = idx, used to convert logical to index vectors.
Node-arc incidence matrix to arc list representation
IJC = incid2list(I,c)
[i,j,c] = incid2list(I,c)
I = m x n node-arc incidence matrix
c = (optional) n-element vector of arc weights
IJC = n x 2-3 matrix arc list [i j c], where
i = n-element vector of arc tails nodes
j = n-element vector of arc head nodes
Transforms: I[i(k),k] = 1 -> i(k)
I[j(k),k] = -sign(j(k)) -> j(k), where
I[j(k),k] = -1 => j(k) > 0, for directed arc (i(k),j(k))
= 1 => j(k) < 0, for undirected arc (i(k),j(k))
= I[j(k),k] = 1 => i(k) = j(k) (self-loop)
Examples:
I = [1 -1 0
-1 0 -1
0 1 1]
c = [1 3 2]
IJC = incid2list(I,c) % IJC = 1 2 1
% 3 1 3
% 3 2 2
I = [1 1 0 % Symmetric arcs
1 0 1
0 1 1]
IJC = incid2list(I,c) % IJC = 1 -2 1
% 3 -1 3
% 2 -3 2
See also LIST2INCID
Input dialog with evaluation of values
[val1,val2,...] = inputevaldlg(prompt,title,defval) % Cell array inputs
val = inputevaldlg(s,title) % Structure input
prompt = cell array of prompt strings for each input value
s = single structure, where the name of each field is prompt string
and default values are the field values (empty values are
ignored)
title = title for the dialog
= name of structure, default for structure input
defval = default values for inputs (optional)
(new dialog is created for any single structure default values)
[val1,val2,...] = output values, if cell array inputs (cell array of
output values, if single output argument)
val = structure of values, if structure input
= [], if dialog cancelled
Examples:
% Calling dialog directly
x = 0;
y = false;
z = 'Hello';
prompt = {'x = ','y = ','z = '};
title = 'Example Dialog';
defval = {x,y,z};
[x,y,z] = inputevaldlg(prompt,title,defval)
% Calling from within a function
function [x,y,z] = myfun(x,y,z)
%MYFUN My function that calls INPUTEVALDLG
if nargin < 1 | isempty(x), x = 0; end
if nargin < 2 | isempty(y), y = false; end
if nargin < 3 | isempty(z), z = 'Hello'; end
if nargin < 1 % Use dialog when no input arguments
prompt = {'x = ','y = ','z = '};
title = 'MYFUN My function that calls INPUTEVALDLG';
defval = {x,y,z};
[x,y,z] = inputevaldlg(prompt,title,defval);
if isempty(x), [x,y,z] = deal([]); return, end % Cancelled dialog
end
% Creating an option structure for MCNF
opt = mcnf('defaults')
opt = inputevaldlg(opt)
(When called within a function, the previous output values are used as
default values for inputs if function has previousely called the dialog.
Same as INPUTDLG except that the input values are evaluated instead of
being returned as strings.)
See also INPUTDLG
Inverse permutation
inva = invperm(a)
a = permutation vector of the integers 1 to length(a)
See also RANDPERM
Inverse Mercator projection
XY = invproj(XY)
y = invproj(y)
XY = two-column matrix of longitude-latitude pairs (in decimal degrees)
y = column vector of latitudes (in decimal degrees)
(Only latitudes are modified in the inverse Mercator projection,
longitudes are unmodified and are not required.)
(Based on Eric Weisstein's Mathworld website
"http://mathworld.wolfram.com/MapProjection.html)
True for zero elements (within tolerance)
y = is0(x,Tol)
= abs(x) < Tol
Tol = tolerance
= [0.01*sqrt(eps)], default
Are XY points in rectangle
isin = isinrect(XY,xy1xy2)
= isinrect(XY) % Graphical input of rectangle using GRECT
xy1xy2 = rectangle
= [min(X) min(Y); max(X) max(Y)]
Note: isin = XY(:,1) >= min(X) & XY(:,1) <= max(X) & ...
XY(:,2) >= min(Y) & XY(:,2) <= max(Y)
True for integer elements (within tolerance)
y = isint(x,TolInt)
= abs(x-round(x)) < TolInt
TolInt = integer tolerance
= [0.01*sqrt(eps)], default
Weighted interlevel to weighted adjacency matrix representation
A = lev2adj(C)
A = lev2adj(C12,C23,...)
C = {C12,C23,...}, cell array of weighted interlevel matrices
Cij = mi x mj matrix from level i to level j
mi = number of nodes in level i
m = m1 + m2 + ..., total number of nodes
A = m x m node-node weighted adjacency matrix of arc lengths
(returns sparse if less than 10% nonzero values)
Note: Predecessor Cij and successor Cjk must be (mi x mj) and (mj x mk)
matrices, respectively.
Zero values in interlevel matrices indicate no arc (must use NaN to
indicate arc with 0 weight).
Examples:
C = [3 4]
A = lev2adj(C) % (1 level) A = 0 3 4
% 0 0 0
% 0 0 0
C12 = C, C23 = [5 6]
A = lev2adj(C12,C23) % (2 levels) A = 0 3 4 0
% 0 0 0 5
% 0 0 0 6
% 0 0 0 0
See also ADJ2LEV, LIST2INCID, ADJ2LIST, and ADJ2INCID
Weighted interlevel to arc list representation
IJC = lev2list(C)
IJC = lev2adj(C12,C23,...)
C = {C12,C23,...}, cell array of weighted interlevel matrices
Cij = mi x mj matrix from level i to level j
mi = number of nodes in level i
m = m1 + m2 + ..., total number of nodes
IJC = n x 2-3 matrix arc list [i j c], where
i = n-element vector of arc tails nodes
j = n-element vector of arc head nodes
c = n-element vector of arc weights
Examples:
C = [3 4]
IJC = lev2list(C) % (1 level) IJC = 1 2 3
% 1 3 4
C12 = C, C23 = [5 6]'
IJC = lev2list(C12,C23) % (2 levels) IJC = 1 2 3
% 1 3 4
% 2 4 5
% 3 4 6
Wrapper for IJC = ADJ2LIST(LEV2ADJ(C))
See also LIST2INCID and ADJ2INCID
Arc list to node-node weighted adjacency matrix representation
A = list2adj(IJC,m,spA)
IJC = n x 2-5 matrix arc list [i j c u l], where
i = n-element vector of arc tails nodes
j = n-element vector of arc head nodes
c = (optional) n-element vector of arc costs, where n = number arcs
= (default) ONES(n,1)
u = (optional) ignored
l = (optional) ignored
m = (optional) scalar size of A if greater than
max{max(i),max(abs(j))}
spA = (optional) make A sparse matrix if n <= spA x m x m
= 1, always make A sparse
= 0.1 (default), A sparse if 10% arc density
= 0, always make A full matrix
A = m x m node-node weighted adjacency matrix
Transforms: If j(k) > 0, then [i(k) j(k) c(k)] -> A[i(k),j(k)] = c(k)
If j(k) < 0, then [i(k) j(k) c(k)] -> A[i(k),-j(k)] = c(k)
A[-j(k),i(k)] = c(k)
Note: Weights of any duplicate arcs added together in A
c(k) = 0 => A(i(k),j(k)) = NaN
Wrapper for c(c==0) = NaN; A = SPARSE(i,j,c,m,m);
Examples:
IJC = [1 2 1
3 1 3
3 2 2]
A = list2adj(IJC) % A = 0 1 0
0 0 0
3 2 0
IJC(3,1) = 0 % Arc (3,1) has 0 weight
A = list2adj(IJC) % A = 0 1 0
0 0 0
NaN 2 0
IJC = [1 -2 1 % Symmetric arcs
3 -1 3
3 -2 2]
A = list2adj(IJC) % A = 0 1 3
1 0 2
3 2 0
See also LIST2INCID, ADJ2LIST, and ADJ2INCID
Arc list to node-arc incidence matrix representation
[I,c,u,l] = list2incid(IJCUL,spI)
IJCUL = n x 2-5 matrix arc list [i j c u l], with i and j required and
c, u, and l optional (just passed through)
i = n-element vector of arc tails nodes
j = n-element vector of arc head nodes
c = n-element vector of arc costs, where n = number of arcs
u = n-element vector of arc upper bounds
l = n-element vector of arc lower bounds
spI = (optional) make I sparse matrix if no. of nodes (m) >= 1/spI
= 1, always make I sparse
= 0.05 (default), make I sparse if m >= 20
= 0, always make I full matrix
I = m x n node-arc incidence matrix
Transforms: i(k) -> I[i(k),k] = 1
j(k) -> I[j(k),k] = -sign(j(k)), where
j(k) > 0 => I[j(k),k] = -1, for directed arc (i(k),j(k))
j(k) < 0 => I[j(k),k] = 1, for undirected arc (i(k),j(k))
i(k) = j(k) (self-loop) => I[i(k),k] = I[j(k),k] = 1
Wrapper for I = SPARSE([1:n 1:n]',[i;abs(j)],[ones(n,1);-sign(j)],n,m)';
Examples:
IJC = [1 2 1
3 1 3
3 2 2]
[I,c] = list2incid(IJC) % I = 1 -1 0 c = 1
-1 0 -1 3
0 1 1 2
IJC(2,3) = 0 % Arc (3,1) has 0 weight
[I,c] = adj2incid(IJC) % I = 1 -1 0 c = 1
-1 0 -1 0
0 1 1 2
IJC = [1 -2 1 % Symmetric arcs
3 -1 3
3 -2 2]
[I,c] = adj2incid(IJC) % I = 1 1 0 c = 1
1 0 1 0
0 1 1 2
See also LIST2ADJ, ADJ2LIST, and ADJ2INCID
Common code for data file functions
varargout = loaddatafile(varargin0,nargout0,vn,fn)
varargin0 = 'varargin' of calling function
nargout0 = 'nargout' of calling function
vn = 'varnames' of calling function
fn = MFILENAME of calling function
Converts routings to weight matrix
W = loc2W(loc,f,h)
loc = routings; e.g., loc = {[1 2 3];[4 1]}, for routing 1-2-3 and 4-1
f = flow
h = handling cost (or equivalence factor)
Example:
loc = {[1 2 3 4],[2 4 1 2 3],[3 4 1 2 4]};
f = [8 5 12];
h = [3 2 1];
W = loc2W(loc,f,h)
Find indices of loc seq segments in arc list
idx = loc2listidx(loc,IJC) % Single arc list
[idx1,idx2] = loc2listidx(loc,IJC1,IJC2) % Multipart arc lists
[idx1,idx2,idx3] = loc2listidx(loc,IJC1,IJC2,IJC3)
loc = vector of single-seq vertices
IJC = n-row matrix arc list
= [IJC1; IJC2; IJC3], arc list composed of up to 3 separate parts
idx = index vector of rows in IJC, where idx(1) is arc IJC(idx(1),:)
corresponding to first segment [loc(1) loc(2)] of loc seq
idx1 = index vector of rows in IJC1
idx2 = index vector of rows in IJC2 offset by size(IJC1,1)
idx3 = index vector of rows in IJC3 offset by size([IJC1; IJC2],1)
Use of a multipart arc list allows reference to be made to the
original arc lists (e.g., in the label roads example, the "LinkTag" is
refinded with respect to the original list IJD, not IJD2).
Examples:
% 4-node graph
XY = [0 0; 1 1; 1 -1; 2 0];
IJC = [1 2 12; 1 3 13; 2 4 24; 3 4 34];
loc = [1 2 4];
idx = loc2listidx(loc,IJC) % idx = 1 3
% Label roads along shortest loc seq from Fayetteville to Raleigh
xy1xy2 = [-79 35; -78 36];
[XY,IJD,isXY,isIJD] = subgraph(usrdnode('XY'),...
isinrect(usrdnode('XY'),xy1xy2),usrdlink('IJD'));
s = usrdlink(isIJD);
[Name,XYcity] = uscity50k('Name','XY',isinrect(uscity50k('XY'),xy1xy2));
[IJD11,IJD12,IJD22] = addconnector(XYcity,XY,IJD);
makemap(XY)
pplot([IJD11; IJD12; IJD22],[XYcity; XY],'r-')
pplot(XYcity,'b.')
pplot(XYcity,Name)
idxFay = strmatch('Fayetteville',Name);
idxRal = strmatch('Raleigh',Name);
[d,loc] = dijk(list2adj([IJD11; IJD12; IJD22]),idxFay,idxRal);
[idx11,idx12,idx22] = loc2listidx(loc,IJD11,IJD12,IJD22);
pplot({loc},[XYcity; XY],'y-','Linewidth',2)
isEndTag = any(diff(double(s.LinkTag(idx22,:)),1,1)')'; % Don't use
isEndTag = [[isEndTag; 0] | [0; isEndTag]]; % repeated tags
pplot(IJD(idx22(isEndTag),:),cellstr(s.LinkTag(idx22(isEndTag),:)),XY)
Calculate total cost of loc seq with time windows
TC = locTC(loc,C)
[TC,XFlg,out] = locTC(loc,C,cap,twin,locfeas)
= locTC(loc), use arguments C,twin,... from previous call
and do not check input for errors
loc = vector of single-seq vertices
= m-element cell array of m loc seqs
(to get sum(TC) as output, use vector loc = [loc{:}] as input)
C = n x n matrix of costs between n vertices
TC = m-element vector of loc seq costs, where
TC(i) = Inf if loc seq i is infeasible
Optional input and output arguments used to determine loc seq feasibility:
cap = {q,Q} = cell array of capacity arguments, where
q = n-element vector of vertex demands, with depot q(1) = 0
Q = maximum loc seq load
twin = {ld,TW,st} = cell array of time window arguments, where
ld = n or (n+1)-element vector of loading/unloading
timespans, where
ld(loc(1)) = load at depot
ld(n+1) = unload at depot, if loc(1) == loc(end)
= scalar of constant values "ld" for loc(2) ... loc(end)
and 0 for loc(1); or loc(2) ... loc(end-1) and 0 for
loc(end), if loc(1) == loc(end)
= 0, default
TW = n or (n+1) x 2 matrix of time windows, where
TW(i,1) = start of time window for vertex i
TW(i,2) = end of time window
TW(loc(1),:) = start time window at depot
TW(n+1,:) = finish time window at depot,
if loc(1) = loc(end)
= (n+1)-element cell array, if multiple windows, where
TW{i} = (2 x p)-element vector of p window
(start,end) pairs
st = (optional) m-element vector of starting times at depot
= TW(1,1) or min(TW{1}), default (earliest starting time)
locfeas = {'locfeasfun',P1,P2,...} = cell array specifying user-defined
function to test the feasibility of a single loc seq (in addition
to time windows, capacity, and maximum cost), where LOCTC
argument out(i) along with user-specified arguments P1,P2,...
are passed to function and a logical value should be returned:
isfeas = LOCFEASFUN(out(i),P1,P2,...)
= {'maxTCfeas',maxTC} is a predefined loc seq feasibility function
to test if the total cost of a loc seq (including loading/
unloading times "ld") exceeds the maximum waiting time "maxTC"
(see below for code)
XFlg(i) = exitflag
= 1, if loc seq is feasible
= -1, if infeasible due to capacity
= -2, if infeasible due to time windows
= -3, if infeasible due to user-defined feasibility function
out = m-element struct array of outputs
out(i) = output structure with fields:
.Loc = loc seq indices, loc{i}
.Cost = cost from vertex j-1 to j,
Cost(j) = C(r{i}(j-1),r{i}(j)) and Cost(1) = 0
= drive timespan from vertex j-1 to j
.Demand = demands of vertices on loc seq, q(loc{i})
.Arrive = time of arrival
.Wait = wait timespan if arrival prior to beginning of window
.Start = time loading/unloading started (starting time for
loc seq is "Start(1)")
.LD = loading/unloading timespan, ld(loc{i})
.Depart = time of departure (finishing time is "Depart(end)")
.Total = total timespan from departing vtx j-1 to depart. vtx j
(= drive + wait + loading/unloading timespan)
.EarlySF = earliest starting and finishing times (default starting
time is "st" and default finish. time is "EarlySF(2)")
.LateSF = latest starting and finishing times
For each seq loc{i}, feasibility is determined in the following order:
1. Capacity feasibility: SUM(q(loc{i})) <= Q if feasible
2. Time window feasiblity: [TCi,ignore,outi] = LOCTC(loc{i},C,twin);
TCi < Inf if feasible
3. User defined feasibility: isfeas = LOCFEASFUN(outi,P1,P2,...);
isfeas == true if feasible
Code for example loc seq feasibility function:
function isfeas = maxTCfeas(outi,maxTC)
%MAXTCFEAS Maximum total cost loc seq feasibility function.
% isfeas = maxwaitfeas(outi,maxTC)
% outi = struct array of outputs from LOCTC for single seq i
% (automatically passed to function)
% maxTC = scalar max. total cost (including un/loading times) of loc seq
%
% Loc Seq is feasible if sum(outi.Total) <= maxTC
%
% This function can be used as a template for developing other
% loc seq feasibility functions.
% Input error check
if ~isnumeric(maxTC) || length(maxTC(:)) ~= 1 || maxTC < 0
error('"maxTC" must be a nonnegative scalar.')
end
% Feasibility test
if sum(outi.Total) <= maxTC
isfeas = true;
else
isfeas = false;
end
Examples:
% 4-vertex graph
IJD = [1 -2 1; 2 -3 2; 3 -4 1; 4 -1 1];
C = dijk(list2adj(IJD)) % C = 0 1 2 1
% 1 0 2 2
% 2 2 0 1
% 1 2 1 0
loc = [1 2 3 4 1];
TC = locTC(loc,C) % TC = 5
% Different loc seq, same C from previous call
TC = locTC([1 2 3 4]) % TC = 4
% Time-windows
ld = 0
TW = [ 6 18 % Start time window at depot
8 11
12 14
15 18
18 24] % Finish time window at depot
[TC,XFlg,out] = locTC([1 2 3 4 1],C,[],{ld,TW})
% TC = 8
% XFlg = 1
% out =
% Loc: [1 2 3 4 1]
% Cost: [0 1 2 1 1]
% Demand: []
% Arrive: [0 11 13 14 16]
% Wait: [0 0 0 1 2]
% Start: [10 11 13 15 18]
% LD: [0 0 0 0 0]
% Depart: [10 11 13 15 18]
% Total: [0 1 2 2 3]
% EarlySF: [10 18]
% LateSF: [10 18]
% Not feasible if total cost exceeds 4
[TC,XFlg] = locTC([1 2 3 4 1],C,[],[],{'maxTCfeas',4})
% TC = Inf
% XFlg = -3
[TC,XFlg] = locTC([1 2 3 4 1],C,[],[],{'maxTCfeas',4})
% TC = Inf
% XFlg = -3
Calculate location sequence cost
c = loccost(loc,C)
loc = location vector
C = n x n matrix of costs between n locations
Example:
loc = [1 2 4 3];
C = triu(magic(4),1); C = C + C'
% C = 0 2 3 13
% 2 0 10 8
% 3 10 0 12
% 13 8 12 0
c = loccost(loc,C)
% c = 2
8
12
Determine nearest city given longitude-latitude of location
[idx,dist,drt,dstr] = lonlat2city(XY,city)
XY = n x 2 matrix of lon-lat (in decimal degrees) of n locations
Optional input argument:
city = stucture with fields:
.Name = m-element cell array of m city name strings
.ST = m-element cell array of m 2-char state abbreviations
.XY = m x 2 matrix of city lon-lat (in decimal deg)
= USCITY50K, default (which contains all cities in US with
population of at least 50,000)
Output arguments: 'dstr' is displayed if no output arguments specified
idx = index of nearest city
dist = distance to nearest city (in miles)
drt = direction to nearest city, reported from 0 to 2PI
radians clockwise from north
dstr = display of nearest city to lat-lon (in city, if dist < 4 mi)
= string, if n == 1
= n-element cell array, if n > 1 (use disp(dstr{i}) to display
ith display string)
Example:
xy = zip2lonlat(28711);
lonlat2city(xy) % xy is 13.54 mi E of Asheville, NC
lonlat2city(xy,uscity) % xy is in Black Mountain, NC
See also CITY2LONLAT
Linear programming using the revised simplex method
[x,TC,XFlg,out] = lprsm(c,Alt,blt,Aeq,beq,LB,UB,idxB0,opt) solves
min TC = c'*x s.t. Alt*x <= blt, Aeq*beq == beq, LB <= x <= UB
x
c = n-element vector of variable costs
Alt = mlt x n inequality constraint matrix
blt = mlt-element RHS vector
Aeq = meq x n equality constraint matrix
beq = meq-element RHS vector
LB = n-element lower bound vector
= [0], default
UB = n-element upper bound vector
= [Inf], default
idxB0 = meq-element vector of initial basis indices
(problem must be in standard form: all Aeq, LB = 0, UB = Inf)
= [] (default) => Use Phase I to find initial basis
Options:
opt = options structure, with fields:
.Display = 2, show steps of evaluation
= 1, warning messages (default)
= -1, no warning messages
.Pivot = method used to select pivot variable
= 1, (default) Steepest edge
= 2, Bland's Rule (prevents cycling)
= 3, Dantzig's Rule
.Tol = tolerance for variable at zero
= 1e-8, default
.MaxCond = maximum initial basis condition estimate
= 1e8, default
.MaxIter = maximum number of simplex iterations before switching
to Bland's Rule ('Pivot' = 2)
= n + 1, default
= 'Field1',value1,'Field2',value2, ..., where multiple input
arguments or single cell array can be used instead of the full
options structure to change only selected fields
opt = LPRSM('defaults'), to get defaults values for option structure
XFlg = 1, solution found
= 0, maximum number of iterations reached
= -1, initial basis is ill conditioned or infeasible
= -2, problem is infeasible
= -3, problem is unbounded
= -4, redundant row in constraints
out = output structure, with fields:
.idxB = final basis indices
.r = reduced costs
.w = dual variables
.iter = number of simplex iterations
(Based on revised simplex method reported in Sec. 3.7 of S.C. Fang and
S.Puthenpura, Linear Optimization and Extensions: Theory and Algorithms,
Pretice-Hall: Englewood Cliffs, NJ, 1993.)
(Jeffery A. Joines developed the initial version of the RSM subroutine
used in this procedure.)
Create projection plot of World or US
h = makemap(region)
= makemap(XY,expand)
region = 'World', (default) world coastline and international borders
= 'US', United States coastline and state and international borders
= 'NC', North Carolina state border
XY = longitude-latitude pairs (in decimal degrees) used to fix axes
limits by calling BOUNDRECT(XY,expand) to get bounding rectangle
expand = nonnegative expansion factor for bounding rectangle (if XY input)
= 0.10, default
h = handle of line objects created
h(1) = coastlines
h(2) = international borders
h(3) = (World) arctic (66.5) and antarctic (-66.5) circles and tropics
of Cancer (23.5) and Capricorn (-23.5) parallels
= (US) State borders
Calls MAPDATA to get map data
Data for maps of the World, US, and North Carolina
s = mapdata Output all variables as structure 's'
[x1,x2] = mapdata Output only first 1, 2, etc., variables
s = mapdata('x',...) Output only variables 'x', ... as structure 's'
[x,...] = mapdata('x',...) Output variables 'x', ... as variables
= mapdata(...,is) Output subset 'x(is)' using SUBSETSTUCT(s,is)
where 'is' is vector of elements to extract
Loads data file "mapdata.mat" that contain the following variables:
World = structure of World data with fields:
.XYB = international borders
.XYC = world coastline
US = structure of United States data with fields:
.XYB = US international borders
.XYC = US coastline
.XYS = US state borders
NC = structure of North Carolina data with fields:
.XYS = NC state border
Data source: World from World Coast Line (1:5,000,000) and US from World
Data Bank II (1:2,000,000; N 53.61 S 14.9 W -125 E -65) of the Coastline
Extractor created by Rich Signell of US Geological Survey and hosted by
NOAA/National Geophysical Data Center, Marine Geology & Geophysics
Division (http://rimmer.ngdc.noaa.gov/coast/). Data thined and sorted
using REDUCEM from Mapping Toolbox and JOIN_CST from Coastline Extractor.
Convert columns of matrix to vectors
[X(:,1),X(:,2),...] = mat2vec(X) (Additional output vectors assigned as empty)
Automatically download updates to Matlog files
nfilesrepl = matlogupdate()
nfilesrepl = number of files updated (command line prompt suppressed if
this output argument is requested)
Note: Do not need to run this if connected to the NCSU network.
If off-campus, can add MATLOGUPDATE to your "startup.m" file.
Uses Matlog function GETFILE to copy files.
(Files that need updating are determined by checking their modification
dates to see if they are earlier than the dates of the same files in the
current Matlog directory at NCSU (http://www.ise.ncsu.edu/kay/matlog/
matlogdir.mat). Looks for the file "matlogupdate.m" to determine the
path to the Matlog directory on the local computer.)
Determine maximum payload
qmax = maxpayld(s,tr)
= maxpayld(sh,tr)
s = vector of shipment density (lb/ft^3)
sh = structure array with field:
.s = shipment density (lb/ft^3)
tr = structure with fields:
.Kwt = weight capacity of trailer (tons)
= Inf, default
.Kcu = cube capacity of trailer (ft^3)
= Inf, default
qmax = maximum payload (tons)
Minimum cost network flow problem
[f,TC,XFlg,nf,out] = mcnf(IJCUL,SCUL,opt)
IJCUL = [i j c u l], arc data
SCUL = [s nc nu nl], node data
i = n-element vector of arc tails nodes, where n = number of arcs
j = n-element vector of arc head nodes
c = n-element vector of arc costs
u = n-element vector of arc upper bounds
= [Inf], default
l = n-element vector of arc lower bounds
= [0], default
s = m-element vector of net node supply, where m = number of nodes
s(i) > 0 => node i is supply node (outflow > inflow)
s(i) = 0 => node i is transshipment node (outflow = inflow)
s(i) < 0 => node i is demand node (outflow < inflow)
nc = m-element vector of node costs
nu = m-element vector of upper bounds on total node flow
= [Inf], default
nl = m-element vector of lower bounds on total node flow
= [-Inf], default
f = n-element vector of arc flows
nf = m-element vector of node flows
Options:
opt = options structure, with fields:
.LP = linear programming solver
= 1, only use LPRSM
= 2, only use LINPROG, if available
= 3, (default) use LPRSM for n <= LrgProb, else LINPROG
.LrgProb = number of arcs in large problem instance
= 1000, default
.TolInt = tolerance for rounding to integer using ISINT
= [1e-6], default
.LPargOut = return c,Alt,blt,Aeq,beq,l,u struct w/o running LP
= 1, do
= 0, do not (default)
.Lprsm = option structure for LPRSM
.Linprog = option structure for LINPROG, with MCNF default:
.Display = 'none'
= 'Field1',value1,'Field2',value2, ..., where multiple input
arguments or single cell array can be used instead of the full
options structure to change only selected fields
opt = MCNF('defaults'), to get defaults values for option structure
TC = total cost
XFlg = exitflag from LP solver
out = output structure, with fields:
.LP = LP solver used
.Lmb = set of Lagrangian multipliers from LINPROG
.Outrsm = output structure from LPRSM
Matrix display
str = mdisp(X,row,col,rowtitle,fracdig,allfrac,isnosep)
X = m x n array of real numbers
row = m-element cell array of row heading strings
= m-element numeric array of row numbers
= numbers 1 to m, default
col = n-element cell array of column heading strings
= n-element numeric array of column numbers
= numbers 1 to n, default
rowtitle = scalar or string title for row headings
= name of first input argument, default
fracdig = digits for fractional portion of number
= 2, default
allfrac = digits for fractional portion of all fractional columns
= 4, default
isnosep = copy output to clipboard without horizontal and vertical
separators (to allow pasting into Excel as fixed width text)
= true, default
str = output as a string (string is also copied to clipboard)
(Adds comma separators to numbers in array.)
Examples:
X = rand(5,4);
mdisp(X)
mdisp(X < .5)
mdisp(X * 10000)
mdisp(fliplr(X),{'R1','row 2','R 3','RowFour','Row 5'},size(X,2):-1:1)
Material handling equipment selection MILP model
[X,Y,TC] = mheselect(c,C,t,T,display,LPargOut)
c = n x m matrix of variable costs of n equipment types and m moves
C = n-vector of equipment fixed costs; C(i) = 0 => fixed path equip
t = n x m matrix of equipment utilizations per move
T = n-vector of equipment availabilities
X = n x m matrix of equipment-move assignments
Y = n x 1 vector of equipment units requirements
TC = total move costs per period (optimal obj. function value)
display = 1, to show MILP steps
= 0, don't show (default)
LPargOut = 1, to only return MILP inputs
= 0, don't return (default)
Mixed-integer linear programming
[x,TC,nevals,XFlg] = milp(c,A,b,vlb,vub,N,R,opt)
N = First N constraints equality constraints (=0, default)
R = Last R variables not integer constrained (=0, default)
Options:
opt = options structure, with fields:
.Display =-1, no warning messages
= 0, initial LP message (default)
= 1, show steps of B&B evalutation
= 2, also show TolFun changes
.Martin = 0, don't use Martin cuts
= 1, use Martin cuts if possible (default)
.SOS = Special Ordered Sets of type one (SOS1)
= 0, don't use (default)
= 1, use if possible (and don't use Martin cuts)
.RCthresh = reduced cost threshold percent of obj for Martin cuts
= 0.05 (default)
.TolInt = integer tolerance for ISINT
= [1e-8], default
.ExTolFun = opt.LPopt.TolFun expansion factor for LINPROG feas.
= 10 (default)
.LPopt = option structure for LINPROG, with MCNF default:
.Display = 'none'
= 'Field1',value1,'Field2',value2, ..., where multiple input
arguments or single cell array can be used instead of the full
options structure to change only selected fields
opt = MILP('defaults'), to get defaults values for option structure
nevals = number of LP evaluations
XFlg = 2, initial LP solution feasible
= 1, solution found
=-1, no feasible integer solution found
=-2, initial LINPROG exitflag <= 0
Min cost insertion of route i into route j
[rte,TC] = mincostinsert(rtei,rtej,rteTC_h,sh,doNaN)
rtei = route i vector
rtej = route j vector
rteTC_h = handle to route total cost function, rteTC_h(rte)
sh = structure array with fields:
.b = beginning location of shipment
.e = ending location of shipment
doNaN = return NaN for rte if cost increase
= false, default
rte = combined route vector
= NaN if cost increase
TC = total cost of combined route
TCij = original total cost of routes i and j
Continuous minisum facility location
[X,TC,XFlg,X0] = minisumloc(P,W,p,V,X0,opt)
Determines locations X that minimizes TC = SUMCOST(X,P,W,p,V).
Uses best solution from gradient-based and Nelder-Mead searches.
P = m x d matrix of m d-dimensional points
W = n x m matrix of X-P weights
p = distance parameter for DISTS in SUMCOST
V = n x n matrix of X-X weights for multifacility problem
(all elements of V used)
= [] (default), solve n single-facility problems
X0 = n x d matrix of starting points
= RANDX(P,n) (default)
opt = options structure, with fields:
.SearchTech = search techniques to use
= 1 => only FMINUNC (gradient-based search)
= 2 => only FMINSEARCH (Nelder-Mead direct search)
= 3 => best of both (default)
.MajorityChk = check for majority single-facility solution
(i.e., X = P(FIND(W >= SUM(W)),:))
= 1, do check (default)
= 0, do not check
.e = epsilon parameter for DISTS passed through SUMCOST
= 1e-4 (default)
.IterPlot = plot each X evaluated by SUMCOST
(use 'delete(findobj('Tag','iterplot'))' to erase)
= h, handle of axes for plotting
= [] (defualt), no plotting
.NormLonLat = normalize lon-lat (i.e., X = NORMLONLAT(X))
= 1, normalize (default)
= 0, do not normalize
.Uopt = option structure for FMINUNC, with MINISUMLOC defaults:
.Display = 'off',.GradObj = 'on',.LargeScale = 'off',
.LineSearchType = 'cubicpoly'
.Sopt = option structure for FMINSEARCH, with MINISUMLOC defaults:
.Display = 'off',.TolFun = 1e-6,.TolX = 1e-6
= 'Field1',value1,'Field2',value2, ..., where multiple input
arguments or single cell array can be used instead of the full
options structure to change only selected fields
opt = MINISUMLOC('defaults'), to get defaults values for option struct.
XFlg = 3 => majority solution found
= 2 => best solution found using FMINSEARCH
= 1 => best solution found using FMINUNC
= 0 => maximum number of iterations reached, if best solution found
using FMINSEARCH; or maximum number of function evaluations
reached, if best solution found using FMINUNC
= -1 => did not converge to a solution, if only FMINUNC used
= -2 => all 0 row in W matrix (returns X(i,:) = X0(i,:); TC(i) = 0)
Examples:
P = [1 1;2 0;2 3], w = [4 3 5] % Single-facility location
[X,TC,XFlg,X0] = minisumloc(P,w,2) % X = 1.2790 1.1717
W = [4 3 5;2 3 0], V = [0 1;0 0], % Multifacility location
[X,TC,XFlg,X0] = minisumloc(P,W,2,V) % X = 1.3197 1.1093
% 2.0000 0.0000
V = [], % Two single-facility locations
[X,TC,XFlg,X0] = minisumloc(P,W,2) % X = 1.2790 1.1717
% 2.0000 0
Minimum weight spanning tree using Kruskal algorithm
[t,nk] = minspan(IJC)
IJC = n x 3 matrix arc list [i j c] of arc heads, tails, and costs
t = n-element logical vector, where
t(i) = 1, if IJC(i,:) arc in spanning tree
t(i) = k, if IJC(i,:) arc in component k of forest
nk = number of components
North Carolina cities with populations of at least 10,000 data
s = nccity Output all variables as structure 's'
[x1,x2] = nccity Output only first 1, 2, etc., variables
s = nccity('x',...) Output only variables 'x', ... as structure 's'
[x,...] = nccity('x',...) Output variables 'x', ... as variables
= nccity(...,is) Output subset 'x(is)' using SUBSETSTUCT(s,is)
where 'is' is vector of elements to extract
Loads data file "nccity.mat" that contain the following variables:
Name = m-element cell array of m city name strings
XY = m x 2 matrix of city lon-lat (in decimal deg)
Pop = m-element vector of total population estimates (2000)
Example: Extract name of all cities in North Carolina
NCcity = nccity('Name')
% NCcity = 'Albemarle'
% ...
% 'Winston-Salem'
(Subset of USCITY10K)
Derived from Source:
http://www.census.gov/ftp/pub/tiger/tms/gazetteer/places2k.txt
Convert longitude and latitude to normal form
XY = normlonlat(XY)
XY = n x 2 matrix of longitude-latitude pairs (in decimal degrees)
Normal form: -180 < X <= 180 and -90 <= Y <= 90
Create cell array of strings from numerical vector
c = num2cellstr(v)
v = n-element numerical vector
c = n-element cell array of strings
Wrapper for cellstr(strjust(num2str(v(:)),'left'))
Validate 'opt' structure
opt = optstructchk(optdef,varargin)
optdef = default opt structure
= MFUN('defaults'), to get default for function MFUN
optin = input opt cell from varargin
opt = structure, if valid
= string, if error
Note: Option values are not vailidated, just its structure
Used in MINISUMLOC and MCNF
Convert rows of NaN-padded matrix to cell array
C = padmat2cell(M)
Calculate pairwise savings
[IJS,S] = pairwisesavings(rteTC_h,sh,TC1)
rteTC_h = handle to route total cost function, rteTC_h(rte)
sh = structure array with fields:
.b = beginning location of shipment
.e = ending location of shipment
TC1 = (optional) user-supplied independent shipment total cost
= rteTC_h([i -i]) for shipment i, default
IJS = savings list in nonincreasing order, where for row IJS(i,j,s)
s is savings associated with adding shipments i and j to route
= [], default, savings
S = savings matrix
Drawnow and then pause after plotting
pauseplot(pt), drawnow and pause for 'pt' seconds (if current
figure UserData is empty)
pauseplot, drawnow and pause using 'pt' from current
figure's UserData (set 'pt' to 0 if empty)
pauseplot('set'), prompt user to input 'pt' into UserData
without drawnow and pause
pauseplot('set',pt), input 'pt' into UserData without drawnow and
pause
pt = pauseplot('get'), get 'pt' from UserData without drawnow & pause
pt = pause time (in seconds)
= Inf, to pause until any key is pressed
pt = [], returned if input canceled or invalid 'pt' value in UserData
Calls DRAWNOW followed by PAUSE(pt)
Plot shipments
h = plotshmt(sh,XY,rte,tr,doproj,maxlab)
sh = structure array with fields:
.b = beginning location of shipment
.e = ending location of shipment
.isLTL = (optional) true if LTL used for independent shipment
XY = 2-column matrix of location lon-lat (in decimal deg)
rte = (optional) route vector
= m-element cell array of m route vectors
tr = (optional) structure with fields:
.b = beginning location of truck = sh(rte(1)).b, default
.e = ending location of truck = sh(rte(end)).e, default
doproj = (optional) do projection using MAKEMAP
= true, default
maxlab = maximum characters on each L/U side of arc label
= 6, default
h = [hTL hLTL hRte] = handles
Projection plot
h = pplot('proj',...) % Current and subsequent plots w/ projection
= pplot(XY) % Plot lines and points
= pplot(XY,LineSpec)
= pplot(XY,LineSpec,'PropName',PropValue,...)
= pplot(XY,Labels,'PropName',PropValue,...) % Label points XY
= pplot(XY,'NumNode',...) % Number nodes
= pplot(IJ,XY,...) % Plot graph
= pplot(IJ,Labels,XY,...) % Label arcs IJ
= pplot(loc,XY,...) % Plot loc seq
= pplot(loc,Labels,XY,...) % Label loc seq 'loc'
= pplot([IJ | loc],XY,Labels,...) % Best-fit labels
= pplot(border,...) % Border offset percentage for XLim and YLim
= pplot('proj') % Create new figure & axes with projection
= pplot % Create without projection
'proj' = initialize current axes so that current and subsequent plots
are projected using PROJ; if axes not initialized, then PPLOT
plots data unprojected (axes initialized by setting Tag
property to 'proj')
XY = m x 2 matrix of m 2-D points, which, if projected, should be
longitude-latitude pairs (in decimal degrees)
IJ = n-row matrix arc list [i j] of arc heads and tails
loc = cell array of loc seq vectors, loc = {[loc1],[loc2],...}
(label arcs of single-seq and centroids of multiple loc seqs)
LineSpec = character string to display line object (see PLOT for options)
Labels = cell array of strings
= num2cell(1:m), if numbering XY
= num2cell(IJD(:,3)), if labeling arc costs in n x 3 matrix IJD
= cellstr(chararray), if using rows of character array as labels
PropName = any line or text object property name
PropValue = corresponding line object property value
(note: a structure or cell arrays can also be used for
property name-value pairs (see SET for details))
border = nonnegative scalar offset percentage for point/arc border,
where offset is percentage of max X,Y extent (only increases
XLim,YLim to MaxXYLim, which can be set as the field of a
structure stored in the "UserData" of the current axes (see
MAKEMAP))
= 0 (default, with projection), no offset
= 0.2 (default, no projection), 20% offset
h = handle of line object created, or, if labels, vector of text
object handles created, or, if PPLOT('proj'), handle of axes
Examples:
XY = [2 2; 3 1; 4 3; 1 4]; % Points
pplot(XY,'r.')
pplot(XY,num2cell(1:4))
% Names of North Carolina cities
[Name,XY] = nccity;
makemap(XY)
pplot(XY,'r.')
pplot(XY,Name)
IJ = [1 2; 1 3; 1 4]; % Arc List
h1 = pplot(IJ,XY,'g-');
IJlab = {'(1,2)','(1,3)','(1,4)'};
h2 = pplot(IJ,IJlab,XY);
delete([h1; h2])
loc = {[1 4 3 2 1]}; % Loc Seq Vector
h3 = pplot(loc,XY,'g-');
loclab = {'(1,4)','(4,3)','(3,2)','(2,1)'};
h4 = pplot(loc,loclab,XY);
set(h3,'color','b')
If current figure does not exist, then new figure created with
DOUBLEBUFFER ON and MATLOGMENU called. If current axes does not exist,
then new axes created with HOLD ON, AXIS EQUAL, and BOX ON.
If IJ or loc not specified for labeling points, then the Delaunay
triangulation of the points are used to fit the labels. Edges greater
than two times the shortest edge incident on each point are not used.
Convert predecessor indices to shortest paths from 's' to 't'
loc = pred2path(P,s,t)
P = |s| x n matrix of predecessor indices (from DIJK)
s = FROM node indices
= [] (default), paths from all nodes
t = TO node indices
= [] (default), paths to all nodes
loc = |s| x |t| cell array of paths (or loc seqs) from 's' to 't', where
loc{i,j} = path from s(i) to t(j)
= [], if no path exists from s(i) to t(j)
(Used with output of DIJK)
Mercator projection
XY = proj(XY)
y = proj(y)
XY = two-column matrix of longitude-latitude pairs (in decimal deg.)
y = column vector of latitudes (in decimal degrees)
(Only latitudes are modified in the Mercator projection, longitudes are
unmodified and are not required.)
(Based on Eric Weisstein's Mathworld website
"http://mathworld.wolfram.com/MapProjection.html)
Project tick marks on projected axes
out = projtick(h)
h = axes handle
= current axes, default
out = 1, ticks projected, if axes created using PPLOT('proj')
= 0, otherwise
Random generation of n points X that are within P's boundaries
X = randX(P,n)
Generates n x d matrix X of n random d-dimensional NF locations, where
each dimension is within the min and max EF locations of P (if m = 1,
then each X(i,:) = P).
P = m x d matrix of m d-dimensional EF locations
n = number of NF locations to generate
= 1, default
Example: Two X points within x = 0..2 and y = 0..3
P = [0 0;2 0;2 3]
X = randX(P,2)
Random re-ordering of an array
X = rankreorder(X,r)
X = array
r = scalar between 0 and 1
Initially, X = X(1:n,:); for i = 1:n, if rand < r, then
X(i,:) placed at the end of X.
Determine estimated LTL rate
rLTL = rateLTL(q,s,d,ppi)
q = weight of shipment (tons)
s = density of shipment (lb/ft^3)
d = distance of shipment (mi)
ppi = LTL Producer Price Index
= 104.2, default, 2004 LTL PPI value
rLTL = rate ($/ton-mi)
Convert route to route index vector
idx = rte2idx(rte)
rte = route vector
= m-element cell array of m route vectors
idx = route index vector, such that idx = rte(rte>0)
= m-element cell array of m route index vectors
Example:
rte = [23 15 6 -23 27 17 24 -27 -15 -17 -6 -24];
idx = rte2idx(rte) % idx = 23 15 6 27 17 24
Convert route vector to load cell array
ld = rte2ld(rte)
rte = route vector
ld = cell array of load vectors
Example:
rte = [1 2 -1 3 -3 -2];
ld = rte2ld(rte); ld{:} % ans = 1 2
% ans = 2 3
Convert route to location vector
loc = rte2loc(rte,sh)
= rte2loc(rte,sh,tr) % Include beginning/ending truck locations
rte = route vector
= m-element cell array of m route vectors
sh = structure array with fields:
.b = beginning location of shipment
.e = ending location of shipment
tr = (optional) structure with fields:
.b = beginning location of truck
= sh(rte(1)).b, default
.e = ending location of truck
= sh(rte(end)).e, default
loc = location vector
= m-element cell array of m location vectors
= NaN, degenerate location vector, which occurs if bloc = eloc and
truck returns to eloc before end of route (=> > one route)
Example:
rte = [1 2 -2 -1];
sh = vect2struct('b',[1 2],'e',[3 4]);
loc = rte2loc(rte,sh) % loc = 1 2 4 3
Total cost of route
[TC,Xflg] = rteTC(rte,sh,C)
= rteTC(rte,sh,C,tr)
rte = route vector
= m-element cell array of m route vectors
sh = structure array with fields:
.b = beginning location of shipment
.e = ending location of shipment
= include to check capacity feasibility of route:
.q = shipment weight (tons)
.s = shipment density (lb/ft^3)
= include to add shipment-related loading/unloading timespan to
route time/cost (location-related time should be added to C):
.tL = loading timespan = 0, default
.tU = unloading timespan = 0, default
= include to check time window feasibility of route:
.tbmin = earliest begin time = -Inf, default
.tbmax = latest begin time = Inf, default
.temin = earliest end time = -Inf, default
.temax = latest end time = Inf, default
C = n x n matrix of costs between n locations
tr = (optional) structure with fields:
.b = beginning location of truck = loc(1), default
.e = ending location of truck = loc(end), default
= include to check capacity feasibility of route:
.Kwt = weight capacity of trailer (tons) = Inf, default
.Kcu = cube capacity of trailer (ft^3) = Inf, default
= include to check max TC feasibility of route:
.maxTC = maximum route TC
= include to check time window feasibility of truck:
.tbmin = earliest begin time = -Inf, default
.tbmax = latest begin time = Inf, default
.temin = earliest end time = -Inf, default
.temax = latest end time = Inf, default
TC(i) = total cost of route i = sum of C and, if specified, tL and tU
= Inf if route i is infeasible
XFlg(i) = exitflag for route i
= 1, if route is feasible
= -1, if degenerate location vector for route (see RTE2LOC)
= -2, if infeasible due to excess weight
= -3, if infeasible due to excess cube
= -4, if infeasible due to exceeding max TC (TC(i) > tr.maxTC)
= -5, if infeasible due to time window violation
out = m-element struct array of outputs
out(i) = output structure with fields:
.Rte = route, with 0 at begin/end if truck locations provided
.Loc = location sequence
.Cost = cost (drive timespan) from location j-1 to j,
Cost(j) = C(j-1,j) and Cost(1) = 0
.Arrive = time of arrival
.Wait = wait timespan if arrival prior to beginning of window
.TWmin = beginning of time window (earliest time)
.Start = time loading/unloading started (starting time for
route is Start(1))
.LU = loading/unloading timespan
.Depart = time of departure (finishing time is Depart(end))
.TWmax = end of time window (latest time)
.Total = total timespan from departing loc j-1 to depart. loc j
(= drive + wait + loading/unloading timespan)
Note, costs C(i,i) ignored for same-location portions of route
Example 1:
rte = [1 2 -2 -1];
sh = struct('b',{1,2},'e',{3,4});
C = triu(magic(4),1); C = C + C'
% C = 0 2 3 13
% 2 0 10 8
% 3 10 0 12
% 13 8 12 0
TC = rteTC(rte,sh,C)
% TC = 22
Example 2: Time Windows
sh = vec2struct('b',[2:4],'e',1,'tbmin',[8 12 15],'tbmax',[11 14 18]);
tr = struct('b',1,'e',1,'tbmin',6,'tbmax',18,'temin',18,'temax',24);
C = [0 1 0 0; 0 0 2 0; 0 0 0 1; 1 0 0 0];
rte = [1:3 -(1:3)];
[TC,Xflg,out] = rteTC(rte,sh,C,tr)
% TC = 8
% Xflg = 1
% out = Rte: [0 1 2 3 -1 -2 -3 0]
% Loc: [1 2 3 4 1 1 1 1]
% Cost: [0 1 2 1 1 0 0 0]
% Arrive: [0 11 13 14 16 18 18 18]
% Wait: [0 0 0 1 2 0 0 0]
% TWmin: [6 8 12 15 18 18 18 18]
% Start: [10 11 13 15 18 18 18 18]
% LU: [0 0 0 0 0 0 0 0]
% Depart: [10 11 13 15 18 18 18 18]
% TWmax: [18 11 14 18 24 24 24 24]
% Total: [0 1 2 2 3 0 0 0]
Normalize route vector
nrte = rtenorm(rte)
rte = route vector
= m-element cell array of m route vectors
nrte = normalized route vector such that its indices range from 1 to
length(rte(rte>0))
Example:
rte = [100 22 33 -22 45 500 -33 66 -500 -100 -45 -66];
rtenorm(rte) % ans = 1 2 3 -2 4 5 -3 6 -5 -1 -4 -6
Savings procedure for route construction
[rte,TC] = savings(rteTC_h,sh,IJS,dodisp)
= savings(rteTC_h,sh,IJS,prte_h)
rteTC_h = handle to route total cost function, rteTC_h(rte)
sh = structure array with fields:
.b = beginning location of shipment
.e = ending location of shipment
IJS = 3-column savings list
= pairwisesavings(rteTC_h)
dodisp = display intermediate results = false, default
prte_h = handle to route plotting function, prte_h(rte)
(dodisp = true when handle input)
rte = route vector
= m-element cell array of m route vectors
TC(i) = total cost of route i
Display structure array with all scalar field values
sdisp(s)
M = sdisp(s) % Return matrix of values
s = structure array
Example:
s = vec2struct('a',[1 2],'b',[3 4]);
sdisp(s) % s: a b
% -:------
% 1: 1 3
% 2: 2 4
Steepest Descent Pairwise Interchange Heuristic for QAP
[TC,a] = sdpi(W,D,a0,C,XY)
m = number of machines and sites
W = m x m machine-machine weight matrix
D = m x m site-site distance matrix
a0 = m-element initial assignment vector (optional)
= [], generate random initial assignment using RANDPERM(m)
= if 'a0' is scalar, then perform 'a0' different runs and report
best solution found
C = m x m machine-site fixed cost matrix (optional), where
C(i,j) is the fixed cost of assigning Machine i to Site j
XY = m x 2 matrix of 2-D site locations to generate plots (optional)
Compute position of points along a 2-D spacefilling curve
theta = sfcpos(X,Y,k)
X,Y = single point in the unit square, or
= multiple 2-D points that will be scaled to the unit square
(single points must be in the unit square and are not scaled)
k = (optional) number of binary digits of X and Y to take into
account when computing theta to 2k binary digits (= 8, default)
(Based on algorithm in J.J. Bartholdi and L.K. Platzman, Management Sci.,
34(3):291-305, 1988)
Create routes for shipments not in existing routes
[rte,idx1,TC] = sh2rte(idx,rte,rteTC_h)
= sh2rte(sh,rte,rteTC_h)
idx = index vector of shipments
sh = structure array of shipments
(used just to create idx = 1:length(sh))
rte = (optional) existing routes = [], default
= m-element cell array of m route vectors
rteTC_h = (optional) handle to route total cost function used to display
results to command window if shipments added
idx1 = index vector of single-shipment routes created
TC = route total cost (only provided if rteTC_h inout)
SUB Create subarray
sub(x,is) % Output as 'ans'
sub(x1,x2,...,is) % Output same name subarray in base workspace
[sx1,...] = sub(x1,...,is) % Output subarrays
sx = sub(x,'# < .5') % Match # to x and evaluate result
= sub(x,'#(:,1)') % Extract first column of x
is = logical array used to create subarrays
Create subgraph from graph
[sXY,sIJC,sisXY,sisIJC] = subgraph(XY,isXY,IJC,isIJC)
psisIJC = subgraph(XY,isXY,IJC) % Partial arcs
XY = m x 2 matrix of node coordinates
isXY = m-element logical vector of node subgraph elements
IJC = n-row matrix arc list
isIJC = n-element logical vector of arc subgraph elements
sXY = nodes of subgraph
sIJC = arcs of subgraph (with endpoints renumbered to correspond to sXY)
sisXY = m-element logical vector of node subgraph elements
(returns index vector if isXY is an index vector)
sisIJC = n-element logical vector of arc subgraph elements
(returns index vector if isIJC is an index vector)
pisIJC = logical vector of partial arc elements, if isIJC empty
("partial arcs" are arcs with only one node in XY)
Examples:
% 7-node graph
XY = [0 0
1 1
1 -1
2 0
3 1
3 -1
4 0]
IJC = [1 -2 12
1 -3 13
2 -4 24
2 -5 25
3 -4 34
3 -6 36
4 -5 45
4 -6 46
5 -7 57
6 -7 67]
isXY = logical([0 1 1 1 1 1 0])' % Nodes 2 through 6
isIJC = logical([1 1 0 1 0 1 0 0 1 1])' % Exterior arcs not connected to
% node 4
[sXY,sIJC,sisXY,sisIJC] = subgraph(XY,isXY,IJC,isIJC)
% sXY = 1 1
% 1 -1
% 3 1
% 3 -1
% sIJC = 1 -3 25
% 2 -4 36
% sisXY' = 0 1 1 0 1 1 0
% sisIJC' = 0 0 0 1 0 1 0 0 0 0
pisIJC = subgraph(XY,isXY,IJC)
% pisIJC' = 1 1 0 0 0 0 0 0 1 1 % Partial arcs
% North Carolina (FIPS code 37) Interstate highways (Type 'I')
[XY,IJD,isXY,isIJD]=subgraph(usrdnode('XY'),usrdnode('NodeFIPS')==37,...
usrdlink('IJD'),usrdlink('Type')=='I');
s = usrdlink(isIJD); % Link data structure
makemap(XY)
pplot(IJD,XY,'r-')
Extract subset of each field of a structure
s2 = subsetstruct(s1,is)
s1 = full structure
is = n-element logical vector of field elements to extract, where n is
size of dimension to use for extraction (max of 2 dimensions, and
extracts first dimension if size of both dim are same--use
function twice in succession in this case)
= if index vector, converted to n-element logical vector assuming
maximum dimension of first field as size n
(use logical vector if this is not the case)
s2 = sub structure
Determine total cost TC = W*DISTS(X,P,p) + V*DISTS(X,X,p)
[TC,dF] = sumcost(X,P,W,p,V,e,h)
X = n x d matrix of n d-dimensional points
P = m x d matrix of m d-dimensional points
W = n x m matrix of X-P weights
= m-element vector of weights applied to each X(i,:), if V == []
= ones(1,m) (default), if V == []
p = distance parameter for DISTS
V = n x n matrix of X-X weights (all elements of V used in TC)
= [] (default)
e = epsilon parameter for DISTS
h = handle of axes to plot each X (see opt IterPlot in MINISUMLOC)
TC = scalar, if W is matrix
= n x 1 vector, if W is vector
dF = n x d gradient matrix
Examples:
x = [1 1], P = [0 0;2 0;2 3], w = [1 2 1]
TC = sumcost(x,P,w,1) % TC = 9
X = [1 1;2 1]
TC = sumcost(X,P,w,2) % TC = 9
7
W = [0 2 1;3 2 0], V = [0 1;0 0]
TC = sumcost(X,P,W,1,V) % TC = 19
Thin degree-two nodes from graph
[tIJC,idxIJC] = thin(IJC)
IJC = n-row matrix arc list of original graph
tIJC = tn-row matrix arc list of thinned graph
idxIJC = tn-element cell array, where idxIJC{i} is an index array of arcs
in IJC that comprise arc [tIJC(i,1) tIJC(i,2)]
Note: Any directed arcs in IJC are converted to undirected arcs.
Returns empty tIJC if all of the nodes are thinned.
Examples:
% 9-node graph
XY = [3 5; 2 4; 5 4; 4 3; 1 3; 0 2; 2 2; 4 2; 3 1];
IJC = [1 -2 12; 1 -3 13; 2 -4 24; 3 -4 34; 2 -5 25; 4 -5 45; 2 -6 26;
5 -7 57; 6 -7 67; 4 -8 48; 5 -8 58; 7 -9 79; 8 -9 89]
[tIJC,idxIJC] = thin(IJC);
tIJC, idxIJC{:}
pplot
pplot(tIJC,XY,'y-','LineWidth',4)
pplot(XY,'b.')
pplot(IJC,XY,'b-')
pplot([IJC; tIJC],XY,num2cell(1:size(XY,1)))
[tXY,stIJC] = subgraph(XY,[],tIJC) % Determine thinned XY and
% re-number nodes in tIJC
% Thin North Carolina Interstate highways
[XY,IJD]=subgraph(usrdnode('XY'),usrdnode('NodeFIPS')==37,...
usrdlink('IJD'),usrdlink('Type')=='I');
[tIJD,idxIJD] = thin(IJD);
makemap(XY)
pplot(tIJD,XY,'y-','LineWidth',3)
pplot(IJD,XY,'r-')
[tXY,stIJD] = subgraph(XY,[],tIJD) % Determine thinned XY and
% re-number nodes in tIJD
Total logistics cost for LTL shipments
[TLC,q,TC,IC] = tlcLTL(sh,ppi,q)
sh = structure array with fields:
.f = annual demand (tons/yr)
.d = shipment distance (mi)
.s = shipment density (lb/ft^3)
.a = inventory 0-D fraction
.v = shipment value ($/ton)
.h = annual holding cost fraction (1/yr)
ppi = LTL Producer Price Index (Series ID=PCU484122484122)
= 104.2, default, 2004 LTL PPI value
q = single shipment weight (tons)
(optional, if not specified and sh.a,v,h specified, then output q
minimizes TLC; required, if sh.a,v,h not specified)
TLC = total logistics cost ($/yr)
q = if input q not specified, optimal single shipment weight (tons)
TC = transport cost ($/yr)
IC = cycle inventory cost ($/yr)
= 0, if sh.a,v,h not specified
Total logistics cost for TL shipments
[TLC,q,TC,IC] = tlcTL(sh,tr,q)
sh = structure array with fields:
.f = annual demand (tons/yr)
.d = shipment distance (mi)
.s = shipment density (lb/ft^3)
.a = inventory 0-D fraction
.v = shipment value ($/ton)
.h = annual holding cost fraction (1/yr)
tr = structure with fields:
.r = transport rate ($/mi)
.Kwt = weight capacity of trailer (tons)
= Inf, default
.Kcu = cube capacity of trailer (ft^3)
= Inf, default
q = single shipment weight (tons)
(optional, if not specified and sh.a,v,h specified, then output q
minimizes TLC; required, if sh.a,v,h not specified)
TLC = total logistics cost ($/yr)
q = if input q not specified, optimal single shipment weight (tons)
TC = transport cost ($/yr)
= nq*r*d*f/q, where nq = max(1,ceil(q/qmax)) trucks per shipment
IC = cycle inventory cost ($/yr)
= 0, if sh.a,v,h not specified
Minimum total logistics cost comparing TL and LTL
[TLC,q,isLTL] = tlcmin(sh,tr,ppi,D,rte)
= tlcmin(sh,tr,ppi,D,q)
sh = structure array with fields:
.f = annual demand (tons/yr)
.s = shipment density (lb/ft^3)
.a = inventory 0-D fraction
.v = shipment value ($/ton)
.h = annual holding cost fraction (1/yr)
= required field if D not provided
.d = shipment distance (mi)
= required fields if D provided
.b = beginning location of shipment
.e = ending location of shipment
= optional fields to use to treat single shipment routes as
independent shipements
.TLC1 = minimum TLC for shipment ($/yr)
.q1 = optimal weight (tons)
.isLTL = true if LTL minimizes TLC for shipment
tr = structure with fields:
.r = transport rate ($/mi)
.Kwt = weight capacity of trailer (tons) = Inf, default
.Kcu = cube capacity of trailer (ft^3) = Inf, default
ppi = LTL Producer Price Index = 104.2, default
D = matrix of inter-location distances
rte = route vector of shipments (if specified, sh.TLC1,q1,isLTL are
used for single shipment routes)
= m-element cell array of m route vectors
q = required if sh.a,v,h not specified,
where TLC = min TL or LTL transport cost at q
TLC = minimum total logistics cost ($/yr)
= min TLC(i) for sh(i), if rte not provided
= min TLC for shipments in route, if rte provided
= m-element vector, min TLC(i) for rte{i}
q = optimal single shipment weight (tons)
= q(i) for sh(i), if rte not provided
= q(i) for sh(i) in route, if rte provided
isLTL = logical vector, true if LTL minimizes TLC for shipment
Transportation and assignment problems
[F,TC] = trans(C,s,d), transportation problem
= trans(C), assignment problem
C = m x n matrix of costs
s = m-element vector of supplies
= [1], default
d = n-element vector of demands
= [1], default
F = m x n matrix of flows
TC = total cost
Wrapper for: C(C == 0) = NaN; MCNF(LEV2LIST(C),[s; -d])
Triangle indices to adjacency matrix representation
A = tri2adj(T)
T = n x 3 matrix, where each row defines indices for one of n
triangles
A = (sparse) m x m adjacency matrix, m is the maximum index value
Converts set of triangles defined by X and Y vector indices into A.
Useful for ploting results from DELAUNAY triangulation using GPLOT.
See also CONVHULL and DSEARCH
Convert triangle indices to arc list representation
IJ = tri2list(T)
T = n x 3 matrix, where each row defines indices for one of n
triangles
IJ = 2-column arc list
Converts set of triangles defined by X and Y vector indices into IJ.
Useful for ploting results from DELAUNAY triangulation using PPLOT.
See also CONVHULL and DSEARCH
Find neighbors of a triangle
N = trineighbors(t,T)
t = m-element vector of triangle indices
T = n x 3 matrix of node indices, where each row defines node indices
for one of n triangles (the result from a DELAUNAY triangulation)
N = m x 3 matrix of triangle indices, where each row N(i,:) defines
indices of the neighboring triangles of triangle T(i,:) and
N(i,j) = NaN, if edge [T(i,j) T(i,j+1)] not shared with any other
triangle (e.g., if it is part of the convex hull)
2-optimal exchange procedure for TSP loc seq improvement
[loc,TC] = tsp2opt(loc,C,cap,twin,locfeas,h)
loc = vector of single-seq vertices, where loc(1) = loc(end)
(loc seq can be infeasible)
= m-element cell array of m loc seqs
C = n x n matrix of costs between n vertices
cap = {q,Q} = cell array of capacity arguments, where
q = n-element vector of vertex demands, with depot q(1) = 0
Q = maximum loc seq load
twin = {ld,TW} = cell array of time window arguments
ld = 1, n, or (n+1)-element vector of load/unload timespans
TW = n or (n+1) x 2 matrix of time windows
= (n+1)-element cell array, if multiple windows
locfeas = {'locfeasfun',P1,P2,...} = cell array specifying user-defined
function to test the feasibility of a single loc seq
h = (optional) handle to vertex plot, e.g., h = PPLOT(XY,'.');
use 'Set Pause Time' on figure's Matlog menu to pause plot
TC = m-element vector of loc seq costs
See LOCTC for information about the input parameters
Note: Also checks reverse of loc seq for improvement in addition to two-
edge exchanges
Example:
vrpnc1
C = dists(XY,XY,2);
h = pplot(XY,'r.');
rand('state',100)
loc = [1 randperm(size(XY,1)-1)+1 1];
[loc,TC] = tsp2opt(loc,C,[],[],[],h);
pplot({loc},XY,num2cellstr(1:size(XY,1)))
Convex hull insertion algorithm for TSP loc seq construction
loc = tspchinsert(XY,h)
XY = n x 2 matrix of vertex cooridinates
h = (optional) handle to vertex plot, e.g., h = PPLOT(XY,'.');
use 'Set Pause Time' on figure's Matlog menu to pause plot
loc = (n + 1)-vector of vertices
(a.k.a. the Vacuum Bag Algorithm)
Example:
vrpnc1
h = pplot(XY,'r.');
pplot(XY,num2cell(1:size(XY,1)))
loc = tspchinsert(XY,h);
TD = locTC(loc,dists(XY,XY,2))
Nearest neighbor algorithm for TSP loc seq construction
[loc,TC,bestvtx] = tspnneighbor(C,initvtx,h)
C = n x n matrix of costs between n vertices
initvtx = [] (default) determine NN loc seq for all vertices and report
best
= initial vertex (or vertices) of loc seq
h = (optional) handle to vertex plot, e.g., h = PPLOT(XY,'.');
use 'Set Pause Time' on figure's Matlog menu to pause plot
loc = (n + 1)-vector of vertices, best loc seq if initvtx = []
TC = total cost of loc seq, including C(n,1)
bestvtx = vertex corresponding to minimum TC found
Examples:
vrpnc1
C = dists(XY,XY,2);
h = pplot(XY,'r.');
pplot(XY,num2cell(1:size(XY,1)))
[loc,TC] = tspnneighbor(C,1,h); TC % Initial Velocx 1
[loc,TC,bestvtx] = tspnneighbor(C,[],h); TC, bestvtx % All vertices
Spacefilling curve algorithm for TSP loc seq construction
loc = tspspfillcur(XY,k)
XY = n x 2 matrix of vertex cooridinates
k = (optional) no. of binary digits to use in SFCPOS (= 8, default)
loc = (n + 1)-vector of vertices
(Based on algorithm in J.J. Bartholdi and L.K. Platzman, Management Sci.,
34(3):291-305, 1988)
Example:
vrpnc1
loc = tspspfillcur(XY);
TD = locTC(loc,dists(XY,XY,2))
h = pplot(XY,'r.');
pplot(XY,num2cell(1:size(XY,1)))
pplot({loc},XY,'m')
2-optimal exchange procedure for route improvement
[rte,TC] = twoopt(rtein,rteTC_h,dodisp)
= twoopt(rtein,rteTC_h,prte_h)
rte = route vector
= m-element cell array of m route vectors
rteTC_h = handle to route total cost function, rteTC_h(rte)
dodisp = display intermediate results = false, default
prte_h = handle to route plotting function, prte_h(rte)
(dodisp = true when handle input)
TC(i) = total cost of route i
Uncapacitated facility location
[y,x,z] = ufl(C,c,form)
M = number of EFs
P = number of sites for NFs
C = 1 x P fixed costs, where C(j) cost of locating NF at site j
c = M x P variable costs, where c(i,j) cost for EF i from site j
form = 0, strong formulation (default)
= 1, weak formulation
y = 1 x P NF assignments, where y(j) = 1 if NF at site j
x = M x P EF allocations, where x(i,j) fraction of EF i from site j
z = total cost
US cities data
s = uscity Output all variables as structure 's'
[x1,x2] = uscity Output only first 1, 2, etc., variables
s = uscity('x',...) Output only variables 'x', ... as struct. 's'
[x,...] = uscity('x',...) Output variables 'x', ... as variables
= uscity(...,is) Output subset 'x(is)' using SUBSETSTUCT(s,is)
where 'is' is vector of elements to extract
Loads data file "uscity.mat" that contain the following variables:
Name = m-element cell array of m city name strings
ST = m-element cell array of m 2-char state abbreviations
XY = m x 2 matrix of city lon-lat (in decimal deg)
Pop = m-element vector of total population estimates (2000)
House = m-element vector of total housing units (2000)
LandArea = m-element vector of land area (square miles)
WaterArea = m-element vector of water area (square miles)
StateFIPS = m-element vector of state FIPS code
PlaceFIPS = m-element vector of place FIPS Code
Example: Extract name of all cities in North Carolina
NCcity = uscity('Name',strcmp('NC',uscity('ST')))
% NCcity = 'Aberdeen'
% ...
% 'Zebulon'
Source:
http://www.census.gov/ftp/pub/tiger/tms/gazetteer/places2k.txt
The place file contains data for all Incorporated and Census Designated
places in the 50 states, the District of Columbia and Puerto Rico as of
the January 1, 2000.
US cities with populations of at least 10,000 data
s = uscity10k Output all variables as structure 's'
[x1,x2] = uscity10k Output only first 1, 2, etc., variables
s = uscity10k('x',...) Output only variables 'x', ... as struct. 's'
[x,...] = uscity10k('x',...) Output variables 'x', ... as variables
= uscity10k(...,is) Output subset 'x(is)' using SUBSETSTUCT(s,is)
where 'is' is vector of elements to extract
Loads data file "uscity10k.mat" that contain the following variables:
Name = m-element cell array of m city name strings
ST = m-element cell array of m 2-char state abbreviations
XY = m x 2 matrix of city lon-lat (in decimal deg)
Pop = m-element vector of total population estimates (2000)
Example: Extract name of all cities in North Carolina
NCcity = uscity10k('Name',strcmp('NC',uscity10k('ST')))
% NCcity = 'Albemarle'
% ...
% 'Winston-Salem'
(Subset of USCITY)
Derived from Source:
http://www.census.gov/ftp/pub/tiger/tms/gazetteer/places2k.txt
The place file contains data for all Incorporated and Census Designated
places in the 50 states, the District of Columbia and Puerto Rico as of
the January 1, 2000.
US cities with populations of at least 50,000 data
s = uscity50k Output all variables as structure 's'
[x1,x2] = uscity50k Output only first 1, 2, etc., variables
s = uscity50k('x',...) Output only variables 'x', ... as struct. 's'
[x,...] = uscity50k('x',...) Output variables 'x', ... as variables
= uscity50k(...,is) Output subset 'x(is)' using SUBSETSTUCT(s,is)
where 'is' is vector of elements to extract
Loads data file "uscity50k.mat" that contain the following variables:
Name = m-element cell array of m city name strings
ST = m-element cell array of m 2-char state abbreviations
XY = m x 2 matrix of city lon-lat (in decimal deg)
Pop = m-element vector of total population estimates (2000)
Example: Extract name of all cities in North Carolina
NCcity = uscity50k('Name',strcmp('NC',uscity50k('ST')))
% NCcity = 'Asheville'
% ...
% 'Winston-Salem'
(Subset of USCITY)
Derived from Source:
http://www.census.gov/ftp/pub/tiger/tms/gazetteer/places2k.txt
The place file contains data for all Incorporated and Census Designated
places in the 50 states, the District of Columbia and Puerto Rico as of
the January 1, 2000.
US highway network links data
s = usrdlink Output all variables as structure 's'
[x1,x2] = usrdlink Output only first 1, 2, etc., variables
s = usrdlink('x',...) Output only variables 'x', ... as struct 's'
[x,...] = usrdlink('x',...) Output variables 'x', ... as variables
= usrdlink(...,is) Output subset 'x(is)' using SUBSETSTUCT(s,is)
where 'is' is vector of elements to extract
Based on Oak Ridge National Highway Network, which contains approximately
500,000 miles of roadway in US, Canada, and Mexico, including virtually
all rural arterials and urban principal arterials in the US. It includes
a large attribute set relevant to routing. Version used (hp15) last
updated May 2001.
Loads data file "usrdlink.mat" that contain the following variables:
IJD = n x 3 matrix arc list [i j c] of link heads, tails,
and distances (in miles)
LinkFIPS = n-element vector of link state FIPS code, where no
individual link crosses a state boundary (see USRDNODE
for FIPS codes)
Type = n-element char vector, where
I - Interstate
U - US route
S - State route
T - Secondary state route
J - Interstate related route (business route, ramps)
C - County route number or name
L - Local road name
P - Parkway
R - Inventory route number (not normally posted)
V - Federal reservation road Used in national parks
and forests, and Indian and military reservation
"-" - (dash) Continuation of local road name from
previous
LinkTag = n x 5 char array of route numbers, with qualifiers
N, S, E, W - Directional qualifiers used to
distinguish
between distinct mainline routes, not to
indicate direction of travel on divided highways
or couplets
A - Alternate
BR - Business route or business loop
BY - Bypass
SP - Spur
LP - Loop
RM - Ramp
FR - Frontage road
PR - Proposed
T - Temporary
TR - Truck route
Y - Wye
Heading = n-element char vector of travel directions N, S, E,
or W
Urban = n-element char vector of flags indicating subjective
judgement about degree of urban congestion, where
U - Urban
V - Urban bypasses (usually circumferential routes)
S - Small urban or towns
T - Partially urban; that is, part of the link is
subject to congestion effects limiting speed
(links having 0.5 miles subject to urban speed
reduction)
X - In urban area, but attributes and topology
unchecked
Median = n-element char vector, where
M - Divided highway (with median)
C - Undivided (ie, "centerline")
F - Ferry
Access_control = n-element char vector, where
U - Uncontrolled access
G - Partial control access (some at-grade
intersections)
I - Fully controlled access (all intersections are
grade separated interchanges)
F - Ferry
Number_lanes = n-element vector, usually either 2 or 4, representing
all multilane roads (reported as 0 for ferries)
Traffic_restrict = n-element char vector, where (in order of priority)
Z - No vehicles (generally a passenger ferry)
P - Closed to public use
C - Commercial traffic prohibited
H - Hazardous materials prohibited
T - Large trucks prohibited
Q - Occasionally closed to public
L - Local commercial traffic only
W - Normally closed in winter
Toll = n-element char vector, where
T - Toll road
B - Link contains a toll bridge (or tunnel)
F - Free passage (roads not flagged can be assumed
toll free, but ferries are uncertain unless
marked)
Truck_route = n-element char vector, where
A - State designated route for STAA-dimensioned
vehicles (the "Staggers act" National Network)
B - Federally designated National Network routes for
large commercial vehicles (subset of State
routes)
C - While nominally Federally designated route, may
be construction-related restrictions for long
period
Major_hwy = n-element vector, integer represents 4-digit bit field
1 - Major highways subsystem
2 - Major highways exclusion
3 - 16-ft route flag (not in public version)
4 - Tunnel flag
Pavement = n-element char vector, where
P - Paved
Q - Secondary paved (through traffic discouraged)
G - Gravel, or otherwise below paved quality
D - Dirt or unimproved
F - Ferry
Admin_class = n-element char vector, where
I - Federal-Aid Interstate
P - Federal-Aid Primary
S - Federal-Aid Secondary
U - Federal-Aid Urban
T - Combination of "S" and "U"
N - Not on a Federal-Aid system
F - Direct Federal system
J - Ramps, connecting roadways, collector/
distributors administratively related to the
Interstate system, but not part of Interstate
mainline mileage
Q - Auxiliary roadways related to non-Interstate
primary routes
Function_class = n x 2 char array of functional class identifier, where
Rural Urban
01 - Interstate 11 - Interstate
02 - Principal arterial 12 - Other expressway
13 -
14 - Principal arterial
05 - Major arterial 15 -
06 - Minor arterial 16 - Minor arterial
07 - Major collector 17 - Collector
08 - Minor collector
09 - Local 19 - Local
Combination
_1 - 01 & 11
_2 - 02 & 12
_3 - 02 & 14
_4 - 06 & 12 05 & 14