05_assignment.Rmd

```{r setup, include=FALSE, cache=FALSE}
knitr::opts_chunk$set(echo = TRUE, cache = TRUE)
library(tidyverse)
library(readxl)
```

# Network Assignment and Validation {#chap-assignment}

The purpose of network assignment is to estimate the traffic flows that will 
occur on each highway link, given the highway network and the trips flowing 
from all origins to all destinations (determined by the other three 
steps of the four-step model).

On average, we start from the assumption that people will take the shortest path
available to them. The travel time on a particular road is a function of the
road's capacity as well as its volume. Because the volume is not known when we
start a traffic assignment process, we will have to find the solution iteratively.

## Volume - Delay Functions

A volume-delay function (VDF) calculates the increase of travel time on a
roadway based on the ratio of volume $V$ to capacity $C$. A popular equation in
travel modeling is a function developed by the Buruea of Public Roads (the
predecessor to the US Department of Transportation). The BPR VDF is

\begin{equation}
  t = t_0[1 + \alpha * (V/C)^\beta] 
  (\#eq:bpr-vdf)
\end{equation}

Where $t$ is the travel time on the link, $t_0$ is the base travel time, 
and $\alpha, \beta$ are calibrated parameters. Figure \@ref(fig:bpr-coeffs) 
shows average values for BPR functions obtained from a sample of MPO travel models
of different sizes. As roads become more heavily loaded, the travel time increases
and other routes become more attractive.

```{r bpr-coeffs, fig.cap="Average BPR VDF curves in a sample of MPO models.", echo = FALSE, message = FALSE}
bpr <- function(vc, alpha, beta) { (1 + alpha * (vc)^beta) }
bpr_coefs <- read_csv("data/bpr_coeffs.csv") %>%
  left_join(tibble(Facility = rep(c("Freeway", "Arterial"), each = 131),
                   vc = rep(seq(0, 1.30, by = 0.01), 2))) %>%
  mutate( t1 = bpr(vc, alpha, beta) )

ggplot(bpr_coefs, aes(x = vc, y = t1, color = Facility)) + 
  facet_wrap(~`MPO Size`) +
  geom_line() + theme_bw()
  
```


## Assignment Algorithms

Consider that we have the network below, with two routes between nodes $A$ and
$B$. The bypass is longer initially, but its travel time will grow less quickly
with added volume.

```{r, echo=FALSE, engine='tikz', out.width='90%', fig.ext=if (knitr:::is_latex_output()) 'pdf' else 'png', engine.opts = list(template = "latex/tikz2pdf.tex"), fig.align="center"}
\begin{tikzpicture}
		\node [] (0) at (-9, 2) {};
		\node [] (1) at (0, 2) {};
		\node [] (2) at (-4.75, 5) {};
		\node [] (3) at (-4.5, 1.25) {$t_t = 10 + 0.02 V_t$};
		\node [] (4) at (-4.25, 5) {$t_b = 15 + 0.005 V_b$};
		\node [] (5) at (-10, 1.75) {A};
		\node [] (6) at (1, 1.75) {B};
		\draw (0.center) to (1.center);
		\draw [bend left=45, looseness=1.25] (0.center) to (1.center);
\end{tikzpicture}
```


In general, the operating theory of network assignment is called 
*static user equilibrium*,

> A network is in static user equilibrium if a person cannot find a shorter path 
between their origin and destination. That is, all paths that are used have the
same travel cost, and all longer paths are unused.

In a small and simple network, we could just generate a system of equations that
represent the SUE conditions, and solve for the values that will give us that
loading. In our example, we can write the system of equations as

\begin{align*}
    &   t_b &- 0.005 V_b  &           &= 15 \\
t_t &       &             &- 0.02 V_t &= 10 \\
t_t & - t_b &             &           &= 0\\
    &       &       V_b   & + V_t     &= 1000\\
\end{align*}

We can solve this using our matrix calculation skills from linear algebra. The
SUE assignment is reached when 600 vehicles take the bypass and 400 vehicles take
the through road, because when that happens both the routes have an equal
travel time of 18 minutes.

```{r sue_example}
(A <- matrix(c(0, 1, -0.005, 0, 
              1, 0, 0, -0.02, 
              1, -1, 0, 0,
              0, 0, 1, 1), byrow = TRUE, ncol = 4))
b <- c(15, 10, 0, 1000)
# Ax = b -> x = A^-1 b
solve(A) %*% b
```

It is not going to be feasible to construct this UE matrix for more complex networks.
So engineers have developed heuristic algorithms that iterate to find
a solution that replicates the UE conditions.

### All-or-Nothing
The most basic way to assign trips is with an "all-or-nothing" (AON) assignment.
This simply puts all the trips between $i$ and $j$ on the shortest route. This
is obviously not great in a lot of ways, because it will overload some
roads while leaving other roads completely empty. So if we assign 1000 trips to
this network, the volumes and travel times become

```{r, echo=FALSE, engine='tikz', out.width='90%', fig.ext=if (knitr:::is_latex_output()) 'pdf' else 'png', engine.opts = list(template = "latex/tikz2pdf.tex"), fig.align="center"}
\begin{tikzpicture}
		\node [] (0) at (-9, 2) {};
		\node [] (1) at (0, 2) {};
		\node [] (2) at (-4.75, 5) {};
		\node [] (3) at (-4.5, 1.25) {$t_t = 10 + 0.02 (1000) = 30$};
		\node [] (4) at (-4.25, 5) {$t_b = 15 + 0.005 (0) = 15$};
		\node [] (5) at (-10, 1.75) {A};
		\node [] (6) at (1, 1.75) {B};
		\draw (0.center) to (1.center);
		\draw [bend left=45, looseness=1.25] (0.center) to (1.center);
\end{tikzpicture}
```

We could repeat this process many times, assigning new AON loads to the updated
travel times. This won't converge to anything, but we could take the average of
all the different AON loadings and run with that. It's not perfect, but it's
easy.


### Incremental Assignment

Instead of assigning the flows all at once, we might be able to get a more
realistic loading by loading the flows in increments. Here's how this algorithm
works:

  1. Select increments $p_n$ that sum to 1 (e.g.,  0.4, 0.3, 0.2, and 0.1)
  1. Calculate the travel times on all links
  1. Assign $V_n * p_n$ trips to the network via All-or-Nothing
  1. Return to step 2 with the next increment
  
Assigning the 1000 trips to our two-route network using the increment values
results in the following successive assignments.
  
| Iteration | Increment | Vb  | tb   | Vt  | tt |
|-----------|-----------|-----|------|-----|----|
| 0         |           | 0   | 15   | 0   | 10 |
| 1         | 0.4       | 0   | 15   | 400 | 18 |
| 2         | 0.3       | 300 | 16.5 | 400 | 18 |
| 3         | 0.2       | 500 | 17.5 | 400 | 18 |
| 4         | 0.1       | 600 | 18   | 400 | 18 | 

### Successive Averages (FHWA) Assignment

A big problem with the AON assignment (and with incremental assignment) is the
large jump in travel times between iterations. It also is not guaranteed to
converge to any particular solution, and the outcome is determined by the
assumptions of the number of increments applied.

This can be improved with a method developed by FHWA that is designed to
repeatedly load the network and update the travel times by a diminishing rate.
In this method, the volume on any particular link after iteration  is given by

\begin{equation}
  V_{n} = (1-\phi) * V_{n-1} + \phi F
  (\#eq:fhwa)
\end{equation}

Where $\phi = 1 / n$ and F is the load of an All-or-Nothing assignment. 
As $n$ increases, the relative amount of weight given to the previous assignment
increases relative to the new AON assignment. The following table shows ten iterations
of this algorithm.

| Iterations | Load | phi        | Vt         | tt         | Vb         | tb         |
|------------|------|------------|------------|------------|------------|------------|
| 1          | F    | 1          | 1000       |            |            |            |
|            | V    |            | 1000       | 30         |            | 15         |
| 2          | F    | 0.5        |            |            | 1000       |            |
|            | V    |            | 500        | 20         | 500        | 17.5       |
| 3          | F    | 0.33333333 |            |            | 1000       |            |
|            | V    |            | 333.33     | 16.67      | 666.67     | 18.3333333 |
| 4          | F    | 0.25       | 1000       |            |            |            |
|            | V    |            | 500        | 20         | 500        | 17.5       |
| 5          | F    | 0.2        |            |            | 1000       |            |
|            | V    |            | 400        | 18         | 600        | 18         |
| 6          | F    | 0.16666667 | 1000       |            | 0          |            |
|            | V    |            | 500        | 20         | 500        | 17.5       |
| 7          | F    | 0.14285714 |            |            | 1000       |            |
|            | V    |            | 428.571429 | 18.5714286 | 571.428571 | 17.8571429 |
| 8          | F    | 0.125      |            | 10         | 1000       | 20         |
|            | V    |            | 375        | 17.5       | 625        | 18.125     |
| 9          | F    | 0.11111111 | 1000       |            |            |            |
|            | V    |            | 444.44     | 18.89      | 555.56     | 17.78      |
| 10         | F    | 0.1        |            |            | 1000       |            |
|            | V    |            | 400.00     | 18.00      | 600.00     | 18.00      |

### Frank-Wolfe 

We can represent SUE traffic assignment as a nonlinear optimization problem. 
We want to find the loading of links that minimizes the total delay in the 
system, subject to the UE constraints. We can represent the total vehicle travel
time on a single link as an integral of the travel time function:

\begin{equation}
\int_0^{v} S(x) dx
  (\#eq:optim-delay)
\end{equation}

where $v$ is the volume assigned to the link and $S(x)$ is the volume-delay
function for the link. As the volume increases, the total delay experienced by
each vehicle increases as well. The R function `integrate` computes a numerical 
integral for a function on a finite interval. We can compute the total vehicle
travel time for each link at an equal loading of 500 vehicles as 

```{r integrate-v}
# delay on bypass
vdb <- function(v){15 + 0.005*v} # volume-delay for bypass
integrate(vdb, 0, 500)$value

# delay on town
vdt <- function(v){10 + 0.02*v} # volume-delay for town
integrate(vdt, 0, 500)$value
```

If we plot these two integrals on the same graph (opposing each other), it is 
pretty easy to see that the 600, 400 loading condition minimizes the area under
each curve, and therefore the total travel time on the system.

```{r integrate-flow}
# reverse vdt function
rvdt <- function(x){30 - 0.02 * x}
ggplot() + 
  stat_function(fun = vdb, xlim = c(0,1000), color = "#ECCBAE") + 
  stat_function(fun = rvdt, xlim = c(0,1000), color = "#046C9A") + 
  stat_function(fun = vdb, xlim = c(0, 600), fill = "#ECCBAE", alpha = 0.2, geom = "area") + 
  stat_function(fun = rvdt, xlim = c(600, 1000), fill = "#046C9A", alpha = 0.2, geom = "area") + 
  xlab("Volume on Bypass (Town = 1000 - Bypass)") + ylab("Travel Time") + 
  theme_bw()
```


Keep in mind that though the travel *time* is the same on both links for the 
last vehicle, the total *delay* is different on each link because more people
use the bypass than the town route. We can now set up the rest of the
optimization. Let

  - $v_a =$ vehicles assigned to link $a$.
  - $S_a(v_a) =$ the travel cost on link $a$ as a function of its volume (VDF function)
  - $X_{ij}^r =$ the total number of vehicles traveling from $i$ to $j$ on
  the sum of links that represent route $r$.
  
We want to minimize the total travel cost for all users

\begin{equation}
  \sum_a \int_0^{v_a} S_a(x) dx
  (\#eq:fw-objective)
\end{equation}

subject to the constraints

\begin{align*}
  v_a &= \sum_i \sum_j \sum_r \delta_{ij}^{ar}X_{ij}{r}\\
  \sum_r X_{ij}^r &= T_{ij}\\
  X_{ij}^r &\geq 0 
\end{align*}  
 
In text, the constraints are as follow: the volume on a link is a
sum of the volume on all routes that use that link ($\delta$ is indicator), the 
total of all routes has to equal the total number of trips assigned, and the
paths on a route are not allowed to be negative. This is a *nonlinear programming*
problem. A number of libraries exist that will find the optimal solutions 
to these problems. The `Rsolnp` library finds our solution in 2 iterations.

```{r assign-nlm}
library(Rsolnp)  # nonlinear programming solver library

# total travel time objective function
# x is a vector of volumes on each link
total_time <- function(x){
  # integrate from 0 to estimated volume
  sum(integrate(vdb, 0, x[1])$value, 
      integrate(vdt, 0, x[2])$value)
}

opt <- solnp(
  c(500, 500), # starting values
  total_time, # objective function
  eqfun = sum, eqB = 1000, #the sum of volumes must be 1000
  LB = c(0, 0), # flows must be positive
  UB = c(1000, 1000) #flows cannot exceed total volume
)
```

Various algorithms can be used to find the values of $v_a$ that minimize this
objective function subject to these constraints. A popular algorithm is the 
Frank-Wolfe algorithm, though other algorithms have been developed that converge 
more quickly under different scenarios.

With these algorithms, it is essential to allow the algorithm to converge
appropriately. A measure of the convergence is a statistic called the "relative
gap", or the difference between the assignment at that iteration and an AON
assignment made with the calculated travel times. As this gap becomes smaller,
it means that the difference between travel times on the routes are becoming
closer to each other. The figure below shows the value of the relative gap 
after several thousand iterations in the Washington, D.C. travel model.

```{r relative-gap, fig.cap = "Relative gap after several thousand iterations.", echo = FALSE}
knitr::include_graphics("images/05_relativegap.png")
```

Large networks may take many hours to reach convergence that is acceptable for
policy analysis. There is a large incentive to "cut corners" by shrinking the
maximum number of iterations that are run, but this can lead to strange behavior.

## Homework {-#hw-assignment}


![Network](images/05_network.png)

The figure above^[This is an adaptation of a homework assignment from
Dr.\ John Ivan at the University of Connecticut.] represents a simple four-node
network where 7000 vehicles travel from A to D, and 5000 travel from B to D
(there are no additional trips from C to D). Link travel times for the network
are given by the functions below.

\begin{align*}
t_{AD} =& 20 + 0.01 q_{AD}\\
t_{AC} =& 10 + 0.005 q_{AC}\\
t_{CD} =& 12 + 0.005 q_{CD}\\
t_{BC} =& 7.25 + 0.005 q_{BC}\\
t_{BD} =& 20 + 0.01 q_{BD}
\end{align*}

Question 1: Solve for the user equilibrium (UE) link flows and travel times by
solving a set of simultaneous equations that explicitly define the UE
conditions.  Demonstrate that your solution is the user equilibrium by showing
through example that all UE conditions are satisfied.

```{r assign-ue, eval = FALSE, echo = FALSE}
system <- matrix(c(
  1, 0, 0, 0, 0,-0.01,      0,     0,     0,     0,
  0, 1, 0, 0, 0,    0, -0.005,     0,     0,     0,
  0, 0, 1, 0, 0,    0,      0,-0.005,     0,     0,
  0, 0, 0, 1, 0,    0,      0,     0,-0.005,     0,
  0, 0, 0, 0, 1,    0,      0,     0,     0, -0.01,
 -1, 1, 1, 0, 0,    0,      0,     0,     0,     0,
  0, 0, 1, 1,-1,    0,      0,     0,     0,     0,
  0, 0, 0, 0, 0,    0,      1,    -1,     1,     0,
  0, 0, 0, 0, 0,    1,      1,     0,     0,     0,
  0, 0, 0, 0, 0,    0,      0,     0,     1,     1), 
 nrow=10,ncol=10, byrow=TRUE)
# right-hand vector b
constraints <- matrix(c(20, 10, 12, 7.25, 20, 0, 0,0,7000,5000),
                      nrow=10, ncol=1)
# compute solution x, Ax=b
solve(system, constraints)
```


Question 2: Perform four iterations of All Or Nothing (AON) assignment on the network
and O/D volumes.  Show the link flows and travel times at the end of each
iteration and compute the average link loads and travel times.

Question 3: Perform four iterations of an incremental assignment assignment 
using the increment values 0.4, 0.3, 0.2, and 0.1.  Show the link flows and
travel times at the end of each iteration.

Question 4: Assign trips using the successive averages (FHWA) heuristic.  Show
the link flows and travel times for five successive assignments, and the final
assignment. 

Question 5: Compare these three traffic assignment heuristic approaches to the UE
assignment and to each other.  How do the resulting flow patterns differ (cite
specific differences)?  Which one comes closest to the UE flows? 

<!-- NEW QUESTION 
Construct a nonlinear programming solver to find the optimal loadings for this network.
You will need to construct an objective function that represents the total travel
time on all links, and a set of constraints that defines how the newtork is loaded.
Note that the solver might not succeed from every starting value; try a couple.
-->

```{r assign-nlp, eval = FALSE, echo = FALSE}
library(Rsolnp)  # nonlinear programming solver library

# volume-delay functions for links
tAD <- function(x) { 20   + 0.01  *x}
tAC <- function(x) { 10   + 0.005 *x}
tCD <- function(x) { 12   + 0.005 *x}
tBC <- function(x) { 7.25 + 0.005 *x}
tBD <- function(x) { 20   + 0.01  *x}

# total travel time objective function
# x is a vector of volumes on each link
total_time <- function(x){
  # integrate from 0 to estimated volume
  sum(
    integrate(tAD, 0, x[1])$value, 
    integrate(tAC, 0, x[2])$value, 
    integrate(tCD, 0, x[3])$value, 
    integrate(tBC, 0, x[4])$value, 
    integrate(tBD, 0, x[5])$value
  )
}

# constraints
constr <- function(x){
  z1 <- x[1] + x[2] # origins at A must equal a certain number
  z2 <- x[4] + x[5] # origins at B must equal a certain number
  z3 <- x[3] - (x[2] + x[4]) # flow on CD must equal sum of AC and BC
  return(c(z1, z2, z3))
}

opt <- solnp(
  c(2000, 2000, 2000, 2000, 2000), # starting values
  total_time, # objective function
  eqfun = constr, eqB = c(7000, 5000, 0), # 7k at A, 5k at B, CD - (AC + BC) = 0
  LB = c(0, 0, 0, 0, 0), # flows must be positive
  UB = c(7000, 7000, 13000, 5000, 5000) #flows cannot exceed total volume
)
```


Question 6: You are considering a road widening project in a suburb of a large
metropolitan area (indicated with the blue circle). The difference in loaded
volumes between your base scenario (no-build) and the widening is given in the figure
below. What is a likely explanation for the patterns shown in the figure?

```{r, fig.cap = "Difference in assigned volumes when adding a lane in area with blue circle.", echo = FALSE}
knitr::include_graphics("images/05_convergence.png")
```


## Laboratory Assignment {-}

The highway volume-to-capacity curves in the Roanoke Model have already been
largely calibrated^[To be specific, VDOT has values that they assert for all of
their models]. They use the Bureau of Public Roads (BPR) format,

$$T_c = T_0 * (1 + \alpha (V / C)^\beta)$$
Create a plot showing the values of these curves for varying VOC ratios and 
discuss the implications of the different curves on different facility types in
your report. Note that there are 5 facility types in the BPR table but 11
facility types in the model network. The assignment script files have comments
that build a crosswalk between the two facility type definitions.

For this lab, you will create a model validation report where you 
examine the following:

  - Root mean squared error (RMSE) by facility type, area type, volume group,
  and by screenline. Are there certain classes that are outperforming others? 
  - Observed vs Modeled link volume scatterplots: an X-Y fit
  line by facility type as well as a maximum desirable deviation plot defined in
  NCHRP 765.
  - Geographic distribution of link error. 
  
Comment on the Roanoke model's calibration.