Traffic_flow.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Traffic flow: the Lighthill-Whitham-Richards model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this chapter we investigate a conservation law that models the flow of traffic.  This model is sometimes referred to as the Lighthill-Whitham-Richards (or LWR) traffic model (see <cite data-cite=\"lighthill1955kinematic\"><a href=\"riemann.html#lighthill1955kinematic\">(Lighthill, 1955)</a></cite> and <cite data-cite=\"richards1956shock\"><a href=\"riemann.html#richards1956shock\">(Richards, 1956)</a></cite>).  This model and the corresponding Riemann problem are discussed in many places; the discussion here is most closely related to that in Chapter 11 of <cite data-cite=\"fvmhp\"><a href=\"riemann.html#fvmhp\">(LeVeque, 2002)</a></cite>.\n",
    "\n",
    "This nonlinear scalar problem is similar to Burgers' equation that we already discussed in [Burgers](Burgers.ipynb) in many ways, since both involve a quadratic (and hence convex) flux function.  In this notebook we repeat some of the discussion from [Burgers](Burgers.ipynb) in order to reinforce essential concepts that will be important throughout the remainder of the book."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you wish to examine the Python code for this chapter, please see:\n",
    "\n",
    " - [exact_solvers/traffic_LWR.py](exact_solvers/traffic_LWR.py) ...\n",
    "   [on github,](https://github.com/clawpack/riemann_book/blob/FA16/exact_solvers/traffic_LWR.py)\n",
    " - [exact_solvers/traffic_demos.py](exact_solvers/traffic_demos.py) ...\n",
    "   [on github.](https://github.com/clawpack/riemann_book/blob/FA16/exact_solvers/traffic_demos.py)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide"
    ]
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide"
    ]
   },
   "outputs": [],
   "source": [
    "%config InlineBackend.figure_format = 'svg'\n",
    "from ipywidgets import interact\n",
    "from ipywidgets import FloatSlider, fixed\n",
    "from exact_solvers import traffic_LWR\n",
    "from exact_solvers import traffic_demos\n",
    "from IPython.display import Image"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The LWR model\n",
    "Recall the continuity equation for any density that is advected with a flow: \n",
    "$$\\rho_t + (u\\rho)_x = 0.$$  \n",
    "In this chapter, $\\rho$ represents the density of cars on a road, traveling with velocity $u$.  Note that we're not keeping track of the individual cars, but just of the average number of cars per unit length of road.  Thus $\\rho=0$ represents an empty stretch of road, and we can choose the units so that $\\rho=1$ represents bumper-to-bumper traffic.\n",
    "\n",
    "We'll also choose units so that the speed limit is $u_\\text{max}=1$, and assume that drivers never go faster than this (yeah, right!).  If we assume that drivers always travel at a single uniform velocity, we obtain once again the advection equation that we studied in [Advection](Advection.ipynb).  But we all know that's not accurate in practice -- cars go faster in light traffic and slower when there is congestion.  The simplest way to incorporate this effect is to make the velocity a linearly decreasing function of the density:\n",
    "$$u(\\rho) = 1 - \\rho.$$\n",
    "Notice that $u$ goes to zero as $\\rho$ approaches the maximum density of 1, while $u$ goes to the maximum value of 1 as traffic density goes to zero.  Obviously, both $\\rho$ and $u$ should always stay in the interval $[0,1]$.\n",
    "\n",
    "Here is a plot of this velocity function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Image('figures/LWR-Velocity.png', width=350)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Combining the two equations above, our conservation law says\n",
    "$$\\rho_t + (\\rho (1-\\rho))_x = 0$$\n",
    "with the flux function \n",
    "$$f(\\rho) = \\rho(1-\\rho)$$\n",
    "giving the rate of flow of cars.  Notice how the flux is zero when there are no cars ($\\rho=0$) and also when the road is completely full ($\\rho=1$).  The maximum flow of traffic actually occurs when the road is half full, as the plot below shows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = lambda rho: rho*(1-rho)\n",
    "traffic_demos.plot_flux(f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Like the flux in Bugers' equation, the LWR flux is **nonlinear**, so we again expect to see shock waves and rarefaction waves in the solution.  We can superficially make this equation look like the advection equation by using the chain rule to write it in quasilinear form:  \n",
    "$$f(\\rho)_x = f'(\\rho) \\rho_x = (1-2\\rho)\\rho_x.$$  \n",
    "Then we have  \n",
    "$$\\rho_t + (1-2\\rho)\\rho_x = 0.$$  \n",
    "This is like the advection equation, but with a velocity $1-2\\rho$ that depends on the density of cars.  The value $f'(\\rho)=1-2\\rho$ is referred to as the *characteristic speed*.  This characteristic speed is not the speed at which cars move (notice that it can even be negative, whereas cars only drive to the right in our model)  Rather, it is the speed at which *information* is transmitted along the road.  Notice that the LWR flux is not convex but concave; because of this, the characteristic speed is a decreasing function of the density."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example: Traffic jam\n",
    "\n",
    "What does our model predict when traffic approaches a totally congested ($\\rho=1$) area?  This might be due to construction, an accident or a red light somewhere to the right; upstream of the obstruction, cars will be bumper-to-bumper, so we set $\\rho=1$ for $x>0$ (supposing that traffic has backed up to that point).  For $x<0$ we'll assume a lower density $\\rho_\\ell<1$.  This is another example of a Riemann problem: two constant states separated by a discontinuity.  We have  \n",
    "$$\n",
    "\\rho(x,t=0) = \\begin{cases} \\rho_\\ell & x<0 \\\\\n",
    "                            1 & x>0. \\end{cases}\n",
    "$$  \n",
    "What will happen as time goes forward?  Intuitively, we expect traffic to continue backing up to the left, so the region with $\\rho=1$ will extend further and further to the left.  This corresponds to the discontinuity (or shock wave) moving to the left.  How quickly will it move?  The example below shows the solution (on the left) and individual vehicle trajectories in the $x-t$ plane (on the right)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "interact(traffic_demos.jam, \n",
    "         rho_l=FloatSlider(min=0.1,max=0.9,value=0.5,\n",
    "                           description=r'$\\rho_l$'),\n",
    "         t=FloatSlider(min=0.,max=1.,value=0.5), fig=fixed(0));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the vehicle trajectory plot above shows the motion of cars moving at velocity $u_\\ell = 1 - \\rho_\\ell>0$ as they approach the traffic jam, and at speed $u_r = 1-\\rho_r = 0$ to the right of the shock, where the cars are stationary.\n",
    "\n",
    "Unlike the case of linear advection, the *characteristic speeds*  are different, with $f'(\\rho_\\ell) = 1 - 2\\rho_\\ell$, which could be either negative or positive depending on $\\rho_\\ell$, and $f'(\\rho_r) = -1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Speed of a shock wave: the Rankine-Hugoniot condition\n",
    "\n",
    "In the plot above, we see a shock wave (i.e., a discontinuity) that moves to the left as more and more cars pile up behind the traffic jam.  How quickly does this discontinuity move to the left?\n",
    "\n",
    "We can figure it out by putting an imaginary line at the location of the shock, as shown in the next figure."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let $\\rho_\\ell$ be the density of cars just to the left of the line, and let $\\rho_r$ be the density of cars just to the right.  Imagine for a moment that the line is stationary.  Then the rate of cars reaching the line from the left is $f(\\rho_\\ell)$ and the rate of cars departing from the line to the right is $f(\\rho_r)$.  If the line really were stationary, we would need to have $f(\\rho_\\ell)-f(\\rho_r)=0$ to avoid cars accumulating at the line."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Image('figures/shock_diagram_traffic_a.png', width=350)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, the shock is not stationary, so the line is moving.  Let $s$ be the speed of the shock.  Then as the line moves to the left, some cars that were to the left are now to the right of the line.  The rate of cars removed from the left is $s \\rho_\\ell$ and the rate of cars added on the right is $s \\rho_r$, as shown in this figure:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Image('figures/shock_diagram_traffic_b.png', width=350)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So in order to avoid an infinite density of cars at the shock, these two effects need to be balanced:  \n",
    "\n",
    "$$f(\\rho_\\ell) - f(\\rho_r) = s(\\rho_\\ell - \\rho_r).$$  \n",
    "\n",
    "This same condition was used for Burgers' equation in [Burgers](Burgers.ipynb), and is known as the **Rankine-Hugoniot condition**. It holds for any shock wave in the solution of any hyperbolic PDE (even systems of equations, where the corresponding vector version gives even more information about the structure of allowable shock waves).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Returning to our traffic jam scenario, we set $\\rho_r=1$.  Then we find that the Rankine-Hugoniot condition gives the shock speed\n",
    "$$s = \\frac{f(\\rho_\\ell)-f(\\rho_r)}{\\rho_\\ell-\\rho_r} = \\frac{f(\\rho_\\ell)}{\\rho_\\ell-1} =  -\\rho_\\ell.$$ \n",
    "This makes sense: the traffic jam propagates back along the road, and it does so more quickly if there is a greater density of approaching cars."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example: green light\n",
    "\n",
    "What about when a traffic light turns green?  At $t=0$, when the light changes, there will be a discontinuity, with\n",
    "traffic backed up behind the light but little or no traffic after the light.  With the light at $x=0$, this takes the form of another Riemann problem:  \n",
    "$$\n",
    "\\rho(x,t=0) = \\begin{cases} 1 & x<0, \\\\\n",
    "                            \\rho_r & x>0, \\end{cases}\n",
    "$$  \n",
    "with $\\rho_r = 0$, for example.\n",
    "In this case we don't expect the discontinuity in density to propagate.  Physically, the reason is clear: after the light turns green, the cars in front accelerate and spread out; then the cars behind them accelerate, and so forth.  This kind of expansion wave is referred to as a *rarefaction wave* because the drivers experience a decrease in density (a rarefaction) as they pass through this wave.   Initially, the solution is discontinuous, but after time zero it becomes continuous."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Similarity solutions\n",
    "The exact form of the solution at a green light can be determined by assuming that the solution $\\rho(x,t)$ depends only on $x/t$.  A solution with this property is referred to as a *similarity solution* because it remains the same if we rescale both $x$ and $t$ by the same factor.  The solution of any Riemann problem is, in fact, a similarity solution.  Writing $\\rho(x,t) = \\tilde{\\rho}(x/t)$ we have (with $\\xi = x/t$):  \n",
    "\\begin{align*}\n",
    "    \\rho_t & = -\\frac{x}{t^2}\\tilde{\\rho}'(\\xi) & f(\\rho)_x & = \\frac{1}{t}\\tilde{\\rho}'(\\xi) f'(\\tilde{\\rho}(\\xi)).\n",
    "\\end{align*}  \n",
    "Thus  \n",
    "\\begin{align}\n",
    "    \\rho_t + f(\\rho)_x = -\\frac{x}{t^2}\\tilde{\\rho}'(\\xi) + \\frac{1}{t}\\tilde{\\rho}'(\\xi) f'(\\tilde{\\rho}(\\xi)) = 0.\n",
    "\\end{align}  \n",
    "This can be solved to find\n",
    "\\begin{align}\n",
    "    f'(\\tilde{\\rho}(\\xi)) & = \\frac{x}{t}\n",
    "\\end{align}  \n",
    "or, since $f'(\\tilde{\\rho}) = 1-2\\tilde{\\rho}$,  \n",
    "\\begin{align}\n",
    "    \\tilde{\\rho}(\\xi) & = \\frac{1}{2}\\left(1 - \\frac{x}{t}\\right).\n",
    "\\end{align}  \n",
    "We know that the solution far enough to the left is just $\\rho_\\ell=1$, and far enough to the right it is $\\rho_r$.  The formula above gives the solution in the region between these constant states.  For instance, if $\\rho_r=0$ (i.e., the road beyond the light is empty at time zero), then  \n",
    "\\begin{align}\n",
    "\\rho(x,t) & = \\begin{cases}\n",
    "                1 & x/t \\le -1 \\\\\n",
    "                \\frac{1}{2}\\left(1 - x/t\\right) & -1 < x/t < 1 \\\\\n",
    "                0 & 1 \\le x/t.\n",
    "            \\end{cases}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The plot below shows the solution density and vehicle trajectories for a green light at $x=0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "interact(traffic_demos.green_light,\n",
    "         rho_r=FloatSlider(min=0.,max=0.9,value=0.3,\n",
    "                           description=r'$\\rho_r$'),\n",
    "         t=FloatSlider(min=0.,max=1.), fig=fixed(0));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the $x$-$t$ plane plot above, the black curves show vehicle trajectories while the blue rays are characteristics corresponding to values of $\\rho$ between $\\rho_\\ell = 1$ and $\\rho_r$, with the left-most characteristic following $x=f'(q_\\ell)t$ and the right-most characteristic following $x= f'(q_r)t$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Entropy condition\n",
    "\n",
    "How can we determine whether an initial discontinuity will lead to a shock or a rarefaction? We have already addressed this for scalar equations in [Burgers](Burgers.ipynb).  Recall the Lax Entropy Condition introduced there:\n",
    "\n",
    "- Shocks appear in regions where characteristics converge, as in the traffic jam example above.  \n",
    "- Rarefactions appear in regions where characteristics are spreading out, as in the green light example."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "More precisely, if the solution is a shock wave with left state $\\rho_\\ell$ and right state $\\rho_r$, then it must be that $f'(\\rho_\\ell)>f'(\\rho_r)$.  In fact the shock speed must lie between these characteristic speeds:\n",
    "\n",
    "$$f'(\\rho_\\ell) > s > f'(\\rho_r).$$\n",
    "\n",
    "We say that the characteristics *impinge* on the shock."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is an example showing the characteristics near a shock:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "traffic_LWR.plot_riemann_traffic(0.2,0.4,t=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On the other hand, if $f'(\\rho_\\ell)< f'(\\rho_r)$, then a rarefaction wave results and the initial discontinuity immediately spreads out.\n",
    "\n",
    "Here is an example showing the characteristics in a rarefaction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "traffic_LWR.plot_riemann_traffic(0.4,0.2,t=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": [
     "hide"
    ]
   },
   "source": [
    "## Interactive Riemann solution\n",
    "\n",
    "In the live notebook, the interactive module below shows the solution\n",
    "of the Riemann problem for any inputs $(\\rho_\\ell,\\rho_r)$ with slider bars to adjust these. The characteristics and vehicle trajectories are also plotted."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide"
    ]
   },
   "outputs": [],
   "source": [
    "traffic_LWR.riemann_solution_interact()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Other resources\n",
    "Many other traffic flow models exist.  Some of them are *continuum models*, like the one presented here, that model traffic density and velocity as an aggregate.  Others are *particle* or *agent* models, that simulate individual vehicles.  You can see simulations of the latter kind on [this webpage](http://www.traffic-simulation.de/routing.html).  Shock waves and rarefaction waves naturally appear in such agent-based models too, and viewing some of these simulations may give you better intuition for the development of shocks and rarefactions in the context of traffic jams."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  },
  "toc": {
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}