The Original Frank-Wolfe Algorithm in Octave (MATLAB)
I recently took a course in network optimization where we learned a variety of ways to solve the static traffic assignment problem (TAP). For one of our problem sets, we had to use the Frank-Wolfe algorithm in order to define the user equilibrium (UE) in a standard TAP. So why is this interesting? Well, if you like transportation policy it probably is not, although its use in network optimization has influenced transportation system design for many years now. Most software packages have improved upon the original Frank-Wolfe algorithm.
Before I start, if anybody is struggling to understand the concepts I talk about, I highly recommend Sheffi's textbook which is available for free online (legally!). While the book is a little dated nowadays, it is very straightforward (for a book about the static TAP) and the concepts lay the groundwork for understanding modern-day research into traffic assignment.
Note: Equations require javascript, since they're created with Mathjax.
Background and Main Concepts
Thankfully this is not an academic paper, which means Wikipedia is a perfectly acceptable source for information. Moving on to the basic problem statement, if \(f(x)\) is differential convex real function and \(D\) is a compact convex set, then:
$$\min\ f(x),$$
$$\mathrm{s.t.}\ x\in D.$$
Convexity is incredibly important to many optimization techniques. While the Wikipedia article I linked above is far more rigorous in its definition, it's easy to think about this convexity requirement in two dimensions. Consider the decaying sine wave \(y(t)=e^{-t}\cdot\cos(2\pi t)\). A figure of this function is below.
The minimum point of this function in this range is about -0.6 from the graph. However, imagine your starting point is at 1.2. Moving along the descent direction would get you to the minimum point at about 1.5, which corresponds to a y-value of about -0.21. Without convexity, it is easily possible to end up in a local minimum instead of a global one. Moving on, the biggest advantage of this algorithm is that you only have to solve a small linear sub-problem instead of a full projection step. Essentially, for certain types of problems, Frank-Wolfe is much faster to solve. However, a crucial step in the algorithm is solving the shortest path problem. This algorithm needs to solve this problem hundreds (or thousands) of times, so it needs to be fairly efficient. But how do you actually find the shortest path anyway?
Bellman-Ford-Moore Shortest Path Code
You start by giving your origin a 'potential' (the time it takes to get to that node) of 0, and every other node a potential of some large number (MATLAB supports infinity, so I used that). Please note this is not vectorized code, so you probably don't want to use it when working on any networks with more than 4 or 5 nodes. Our professor gave us his code that was vectorized, but I have not obtained his permission to republish it here, so you are stuck with my inefficient version. MATLAB doesn't have any queue data structure (my friend who works there laughed at me when I asked, MATLAB isn't for this kind of thing), so I used CQueue from the MATLAB Workshop. CQueue is slightly faster (according to random internet forum posters...) than importing Java's queuing structure but I don't have Java installed so that wasn't an option anyway. Once again though, vectorized code will be much faster.
function linkPred = BellmanFord(sN,eN,nn,nl,frstOut,lstOut,bNode,t)
% Function takes in the start and end nodes along with the network
% characteristics pulled from a forward star network
% Inputs:
% sN - Start node
% eN - End node
% nn - Number of nodes
% nl - Number of links
% frstOut - first (by number) link out of a node
% lstOut - last (by number) link out of a node
% bNode - the end node of a link
% t - the travel times on each link
% Outputs:
% linkPred - an array of the precending link for a given node
nodePot = ones(1,nn).*Inf;
nodePot(sN) = 0; % Origin has a potential of zeros
SPLinks = zeros(1,nl);
nodePred = zeros(1,nn);
linkPred = zeros(1,nn);
q = CQueue(sN); % using CQueue package from zhang@zhiqiang.org
A = zeros(1,nn); % creating array A to check which nodes have been added to q
while q.isempty() ~= 1
i = q.pop(); % taking front node and removing from queue
if i == eN
break;
end
% Checking all links going out of node i
for j = frstOut(i):lstOut(i)
endNode = bNode(j); % finding the end node of link j
% Checking node j's potential, and adjusts it if needed
if nodePot(i) + t(j) < nodePot(endNode)
nodePot(endNode) = nodePot(i)+t(j);
nodePred(endNode) = i; % updating predecessors
linkPred(endNode) = j; % updating link predecessors
end
% Adds node j to end of the queue, if it's not already in it
if A(endNode) ~= 1
A(endNode) = 1;
q.push(endNode);
end
end
end
end
Frank-Wolfe Algorithm with Relative Gap Code
Once again, I didn't really bother to vectorize this code. This isn't nearly as big a problem as with the shortest path code though, as the number of iterations for the networks we had to deal with was always less than 250. For a long time, I used the absolute gap instead of the relative gap, and for the larger networks we had to do (400+ nodes), the algorithm wouldn't converge even at 10,000 iterations. For the same problem, it only took 226 iterations to converge. Also, bfsmp.m is my professor's shortest path code, so you would need to change that line if you want to experiment with my code.
function [LBB, OBJB, x] = FrankWolfeEquivGap(fname)
% This function solves the static Traffic Assignment Problem of a forward star
% network using a travel demand file
% Input:
% fname - Name of file, without the extension. For a network called NU, the
% forward star file would bu NU.1, and the travel demand NU.2. Then
% fname would just be 'NU'
% Output:
% LBB - Lower bound on objective function when finished
% OBJB - Objective function value when finished
% x - Array of flows on each link in the network
% Reading in Network Demand data
[no,orgid,startod,nod, dest, od_demand] = read2(strcat(fname,'.2'));
% Reading in FORT data about network
[nn,frstout,lstout,na,anode,bnode,sat,lngth,vmax] = read1(strcat(fname,'.1'));
% Necessary variables
k = 1; % iteration counter
epsilon0 = 0.001; % Acceptable gap amount
epsilon1 = 1*(10^-4); % gap criterion
gap = Inf; % initial gap size
gap_rel = Inf; % initial relative gap size
x = zeros(na,1); % no flows on any of the 24 links to start with
K0 = 10000; % max number of iterations
t0 = lngth./vmax; % free flow travel time
cost = Problem5TT(x, t0, sat); % Initial link costs
% Setting up buckets for later
OBJB = zeros(300,1);
LBB = zeros(300,1);
% Adding another entry to startod so that the for loop can go to the end of
% destinations
startod(no+1)=length(dest)+1;
% Going through all origins
for i = 1:no
% Finding the shortest path from the origin to all destinations
[sP, dist] = bfmsp(nn,frstout,lstout,bnode,i,cost);
for j = startod(i):(startod(i+1)-1)
% Finding the demand for the OD pair i-j
q_ij = od_demand(j);
% Assigning flow to this all links on the shortest path
indx = dest(j);
while indx != i
x(sP(indx)) += q_ij;
indx = anode(sP(indx));
end
end
end
% Actual algorithm
while (gap_rel > epsilon0 && k < K0)
% Updating y, where y is the solution to the linear subproblem
y = zeros(na,1);
cost = Problem5TT(x, t0, sat);
for i = 1:no
% Finding the shortest path from the origin to all destinations
[sP, dist] = bfmsp(nn,frstout,lstout,bnode,i,cost);
% Going through all destinations for origin i
for j = startod(i):(startod(i+1)-1)
% Finding the demand for the OD pair i-j
q_ij = od_demand(j);
% Assigning flow to this all links on the shortest path
indx = dest(j);
while indx != i
y(sP(indx)) += q_ij;
indx = anode(sP(indx));
end
end
end
% Computing the gap
gap = sum(cost.*(y.-x));
% Finding the lower bound,
low_bound = Problem5UEObj(x, t0, sat) + gap;
LBB(k) = low_bound/1000;
OBJB(k) = Problem5UEObj(x, t0, sat)/1000;
gap_rel = -gap/abs(low_bound);
% Line Search Time!
% Initializing necessary values
I = 0; a=0; b=1; alpha=0.5*(a+b);
d = y-x; % Descent direction
% Updating x and t, moving along descent direction with step size alpha
x_alpha = x + alpha*d;
t_alpha = Problem5TT(x_alpha, t0, sat);
s_alpha = sum(t_alpha.*d);
% bisection search
while (((I < K0) && ((b-a) > epsilon1) && (abs(s_alpha) > epsilon1)) || (s_alpha > 0))
% Updating bounds a and b
if s_alpha < 0
a = alpha;
else
b = alpha;
end
% Updating alpha with new boundaries
alpha = 0.5*(a+b);
% Updating x, t, and s_alpha for the next iteration
x_alpha = x + alpha*d;
t_alpha = Problem5TT(x_alpha, t0, sat);
s_alpha = sum(t_alpha.*d);
I++; % updating iteration counter
end
% Progressing!
x = x + alpha*(y-x);
% Updating t
t = Problem5TT(x, t0, sat);
k++; % updating iteration count
end
% Outputting important values
disp 'The objective function value is'
obj = Problem5UEObj(x, t0, sat)
disp 'The relative gap is'
gap_rel
disp 'The number of iterations is'
k
end
So is this useful?
In short, it used to be. Most modelling software nowadays uses more advanced algorithms than the original Frank-Wolfe Algorithm. More importantly though, they do not implement their code in MATLAB/Octave, which is far slower than almost any compiled language. Despite advancements in recent years, the Frank-Wolfe algorithm is still a great introduction to the static traffic assignment problem. While my code is certainly not efficient, hopefully it can still be informative to anyone who is struggling to figure it out.