%% A Model of Traffic Flow
%
% Everyone has had the experience of sitting in a traffic jam, or of seeing cars bunch up 
% on a road for no apparent good reason.  MATLAB and Simulink are good tools for studying 
% models of such behavior.  Our analysis here will be based on so-called ``follow-the-leader'' 
% theories of traffic flow, about which you can read more in 
% _Kinetic Theory of Vehicular Traffic_, by Ilya Prigogine and Robert Herman, Elsevier, 
% New York, 1971, or in _The Theory of Road Traffic Flow_, by Winifred Ashton, Methuen, London, 
% 1966.  We will analyze here an extremely simple model that already exhibits quite 
% complicated behavior.  We consider a one-lane, one-way, circular road with a number of cars
% on it (a very primitive model of, say, the Outer Loop of the Capital Beltway around 
% Washington, DC, since, in very dense traffic, it is hard to change lanes and each lane 
% behaves like a one-lane road).  Each driver slows down or speeds up on the basis of his 
% own speed, the speed of the car ahead of him, and the distance to the car ahead of him.  
% But human drivers have a finite reaction time.  In other words, it takes them a certain 
% amount of time (usually about a second) to observe what is going on around them and to 
% press the gas pedal or the brake, as appropriate.  The standard ``follow-the-leader'' theory 
% supposes that
% 
% $$\ddot u_n(t+T) = \lambda\bigl(\dot u_{n-1}(t) - \dot u_n(t)\bigr),\qquad\qquad(\dag)$$
% 
% where $t$ is time, $T$ is the reaction time, $u_n$ is the position of the $n$th car, and 
% the ``sensitivity coefficient'' $\lambda$ may depend on  
% $u_{n-1}(t) - u_n(t),$
% the spacing between cars, and/or 
% $\dot u_n(t),$
% the speed of the $n$th car.  The idea behind this equation is this.  A driver will 
% tend to decelerate if he is going faster than the car in front of him, or if he is 
% close to the car in front of him, and will tend to accelerate if he is going slower 
% than the car in front of him.  In addition, a driver (especially in light traffic) may 
% tend to speed up or slow down depending on whether he is going slower or faster 
% (respectively) than a ``reasonable'' speed for the road (often, but not always, equal 
% to the posted speed limit).  Since our road is circular, in this equation $u_0$ is 
% interpreted as $u_N$, where $N$ is the total number of cars.
%
% The simplest version of the model is the one in which the ``sensitivity coefficient'' 
% $\lambda$ is a (positive) constant.  Then we have a homogeneous linear 
% differential/difference equation with constant coefficients for the velocities
% $\dot u_n(t)$.
% Obviously, there is a ``steady-state'' solution when all the velocities are equal and 
% constant (i.e., traffic is flowing at a uniform speed), but what we are interested 
% in is the stability of the flow, or the question of what effect is produced by small 
% differences in the velocities of the cars.  The solution of ($\dag$) will be a superposition 
% of exponential solutions of the form 
% 
% $$u_n(t) = \exp(\alpha t)v_n,$$
% 
% where the $v_n$'s and $\alpha$ are (complex) constants, and the system will be unstable 
% if the velocities are unbounded, i.e., there are any solutions where the real part of 
% $\alpha$ is positive.  Using vector notation, we have
% 
% $$ \dot {\mathbf u}(t)= \exp(\alpha t) {\mathbf v},
% \qquad \ddot {\mathbf u}(t+T)= \alpha \exp(\alpha T) \exp(\alpha t) {\mathbf v}.$$
% 
% Substituting back in ($\dag$), we get the equation
% 
% $$\alpha  \exp(\alpha  T) \exp(\alpha  t) {\mathbf v} = 
% \lambda ( S - I ) \exp(\alpha  t) {\mathbf v},$$
%
% where 
% 
% $$S = \begin{pmatrix}0&0&\cdots&0&1\\ 1&0&\cdots&0&0\\ 0&1&\cdots&0&0\\
% \vdots&\vdots&\ddots&\vdots&\vdots\\ 0&0&\cdots&1&0\\\end{pmatrix}$$
% 
% is the ``shift'' matrix that, when it multiplies a vector on the left, cyclically 
% permutes the entries of the vector.  We can cancel the $\exp(\alpha t)$ on each side to get
% 
% $$\alpha \exp(\alpha T) {\mathbf v} = \lambda( S - I ) {\mathbf v}, $$
%
% or
%
% $$\left( S - \left(1 + \frac\alpha\lambda \exp(\alpha T) \right)I \right) {\mathbf v} = 0,  
% \qquad\qquad(**) $$
% 
% which says that $v$ is an eigenvector for $S$ with eigenvalue
% 
% $$1 + \frac\alpha\lambda \exp(\alpha T).$$
% 
% Since the eigenvalues of $S$ are the $N$th roots of unity, which are evenly spaced 
% around the unit circle in the 
% complex plane, and closely spaced together for large $N$, there is potential instability 
% whenever   
% 
% $$1 + \frac\alpha\lambda \exp(\alpha T)$$
% 
% has absolute value 1 for some $\alpha$ with positive real part; that is,
% whenever
% 
% $$\left( \frac{\alpha T}{\lambda T} \right)e^{\alpha T}$$
% 
% can be of the form $\exp(i\theta) - 1$ for some $\alpha T$ with positive real part.  Whether instability 
% occurs depends on the value of the product $\lambda T$.  We can see this by plotting 
% values of $z \exp(z)$ for $z = \alpha T = i y$ a complex number on the critical line Re $z = 0$, 
% and comparing with plots of  $\lambda T ( e^{i\theta} - 1 )$ for various values  
% of the parameter $\lambda T$:

syms y; expand(i*y*(cos(y) + i*sin(y)))

%%

ezplot(-y*sin(y), y*cos(y), [-2*pi, 2*pi]); hold on
theta = 0:0.05*pi:2*pi;
plot((1/2)*(cos(theta) - 1), (1/2)*sin(theta), '-');
plot(cos(theta) - 1, sin(theta), ':')
plot(2*(cos(theta) - 1), 2*sin(theta), '--');
title('iye^{iy} and circles \lambda T(e^{i\theta}-1)')
hold off  

%%
% Here the small solid circle corresponds to $\lambda T = 1/2$, and we are just at the 
% limit of stability, since this circle does not cross the spiral produced by $z \exp(z)$ for 
% $z$ a complex number on the critical line Re $z = 0$, though it ``hugs'' the spiral closely.  
% The dotted and dashed circles, corresponding to $\lambda T = 1$ or $2$, do cross the spiral, 
% so they correspond to unstable traffic flow.
%
% We can check these theoretical predictions with a simulation using Simulink. 
% We'll give a picture of the Simulink model and then explain it:
 
open_system traf

%%
% Here the Subsystem, which corresponds to multiplication by $S - I$, looks like this:

open_system 'traf/Subsystem: computes velocity differences'
 
%%
% Most of the model is like 
% the example in Chapter 8, except that our unknown function (called $u$), representing the car 
% positions, is vector-valued, not scalar-valued.  The major exceptions are these.
% 
% * We need to incorporate the reaction-time delay, so we've inserted a
% *Transport Delay* block from the *Continuous* block library, with the
% ``Time delay'' parameter $T$ set to $0.5$.
% * The parameter $\lambda$ shows up as the value of the gain in the ``sensitivity
% parameter'' *Gain* block in the upper right.
% * Plotting car positions by themselves is not terribly useful, since only
% the relative positions matter.  So before outputting the car positions to the 
% *Scope* block labeled ``relative car positions,'' we've subtracted a constant linear 
% function (corresponding to uniform motion at the average car speed) created by 
% the *Ramp* block from the *Sources* block library.
% * We've made use of the option in the *Integrator* blocks to input the initial conditions, 
% instead of having them built into the block.  This makes the logical structure a little clearer.
% * We've used the subsystem feature of Simulink.  If you enclose a bunch of blocks with 
% the mouse and then click on ``Create subsystem'' in the model's *Edit* menu, Simulink will 
% package them as a subsystem.  This is helpful if your model is large or if there is some 
% combination of blocks that you expect to use more than once.  Our subsystem sends a vector 
% $v$ to $( S - I ) v = S v - v$.  A *Sum* block (with one of the signs changed to a minus) is used for 
% vector subtraction.  To model the action of $S$, we've used the *Demux* and *Mux* blocks from the 
% *Signal Routing* block library.  The *Demux* block, with the ``number of outputs'' parameter set 
% to *|[4, 1]|*, splits a five-dimensional vector into a pair consisting of a four-dimensional vector 
% and a scalar (corresponding to the last car).  Then we reverse the order of these and put 
% them back together with the *Mux* block, with the ``number of inputs'' parameter set to *|[1, 4]|*.
%
% Once the model has been assembled, it can be run with various inputs.  
% You can see the results yourself in the two *Scope* windows,
% but here we've run the simulation from the command line and
% plotted the results with the |*simplot*| command,
% that does almost the same thing as a *Scope* but in a 
% regular MATLAB figure window.  The following pictures
% are produced with $\lambda = 0.8$, corresponding to stable flow (though,
% to be honest, we've let two cars cross through each other briefly!):

set_param('traf/sensitivity parameter', 'Gain', '0.8');
trafout = sim('traf', 'SaveTime','on','TimeSaveName','t', ...
   'SaveState','on','StateSaveName','x');
t = trafout.find('t'); x = trafout.find('x');

%%
% The variable |*t*| stores the time parameter, the variable
% |*x*| stores car velocities in its first five columns and car
% positions in the second five columns.  In this example, the average
% velocity is $3.15$.  First, we plot the relative positions,
% then we plot the velocities.

relpos = x(:,6:10) - 3.15*t*ones(1,5);
simplot(t, relpos), title('Relative Positions')
axis([0, 20, 0, 1]), axis normal
set(gcf, 'InvertHardcopy', 'off')
%%

simplot(t, x(:,1:5)), title('Car Velocities')
axis([0, 20, 3, 3.3])
set(gcf, 'InvertHardcopy', 'off')
%%
% As you can see, the speeds fluctuate but eventually converge to a single value, and the 
% separations between cars eventually stabilize.
% On the other hand, if $\lambda$ is increased by changing the ``sensitivity parameter'' in the 
% *Gain* block in the upper right, say from $0.8$ to $2.0$, we get the following output, 
% which is typical of instability:

set_param('traf/sensitivity parameter', 'Gain', '2.0');
[t, x] = sim('traf');
relpos = x(:,6:10) - 3.15*t*ones(1,5);
simplot(t, relpos), title('Relative Positions')
axis([0, 20, -10, 10])
set(gcf, 'InvertHardcopy', 'off')
%%

simplot(t, x(:,1:5)), title('Car Velocities')
axis([0, 20, -10, 10])
set(gcf, 'InvertHardcopy', 'off')
%%
% We encourage you to go back and tinker with the model (for instance using a 
% sensitivity parameter that is also inversely proportional to the spacing between 
% cars) and study the results.
%
% Finally, you can create a movie with the following code:

%%
clf reset
set(gcf, 'Color', 'White')
clear M
theta = (0:0.025:2)*pi;
for j = 1:length(t)
    plot(cos(x(j, 1:5)), sin(x(j, 1:5)), 'o'); 
    axis([-1, 1, -1, 1]);
    hold on; plot(cos(theta), sin(theta), 'r'); hold off
    axis equal;
    M(j) = getframe;
end

%%
% The idea here is that we have taken the circular road to have radius 1 
% (in suitable units), so that the command |*plot(cos(theta),sin(theta),'r')*| 
% draws a red circle (representing the road) in each frame of the movie, and on
% top of that the cars are shown with moving little circles.
% The graph above is the last frame of the movie; you can view the entire
% movie by typing *|movie(M)|* or *|movieview(M)|*.  Try it!

%%
% We should mention here one fine point needed to create a realistic movie.  
% Namely, we need the values of |*t*| to be equally spaced -- otherwise the 
% cars will appear to be moving faster when the time steps are large, and will 
% appear to be moving slower when the time steps are small.  In its default 
% mode of operation, Simulink uses a variable-step differential-equation solver 
% based on MATLAB's command |*ode45*|, so the entries of |*t*| will not be equally
% spaced.  We have fixed this by opening the *Configuration Parameters* dialog box using the 
% *Simulation* menu in the model window, and, under the
% *Data Import/Export* item, changing
% the *Output options* box to read |Produce specified output only|, 
% with *Output times*
% chosen to be *|[0:0.5:20]|*.  Then the model will output the car 
% positions only at times that are multiples of $0.5$, and the MATLAB
% program above will produce a $41$-frame movie.