Code covered by the BSD License  

Highlights from
CTMSIM - an interactive freeway traffic macrosimulator

image thumbnail
from CTMSIM - an interactive freeway traffic macrosimulator by Alex Kurzhanskiy
Freeway traffic simulation based on Asymmetric Cell Transmission Model

simulation_step(celldata, ...
function [simdata, tt] = simulation_step(celldata, ...
                                         densities, ...
                                         orflows, ...
                                         orqueues, ...
                                         flows, ...
                                         frflows, ...
                                         inflow, ...
                                         outflow, ...
                                         ts, ...
                                         sstep);
% SIMULATION_STEP - computes CTM values (densities; mainline, offramp flows;
%                   speeds; VHT; VMT; travel time; productivity loss) using
%                   initial data or data from previous simulation step.
%
% Call:   simdata = simulation_step(celldata, densities, orflows, orques, flows, frflows, inflow, outflow, ts, sstep)
%
% Parameters:
%             celldata  - array of freeway cells;
%             densities - vector of densities whose size is determined by the number of cells;
%             orflows   - vector of onramp flows from the previous step whose size is determined by the number of cells;
%             orqueues  - vector of onramp queues whose size is determined by the number of cells;
%             flows     - vector of mainline flows whose size is determined by the number of cells;
%             frflows   - vector of offramp flows whose size is determined by the number of cells;
%             inflow    - input flow into the most upstream cell
%             outflow   - mainline flow allowed to leave from the most downstream cell;
%             ts        - sampling time period; 
%             sstep     - simulation step.
%
% Returns:
%          simdata - matrix whose columns are:
%                       * vector of densities,
%                       * vector of mainline flows,
%                       * vector of offramp flows,
%                       * vector of speeds,
%                       * vector of VHT values,
%                       * vector of VMT values,
%                       * vector of productivity loss values;
%           tt     - travel time.
%
% Last modified:   02/24/2008.

%
% Alex Kurzhanskiy   <akurzhan@eecs.berkeley.edu>
%

N = size(celldata, 2);  % number of cells

if size(densities, 1) ~= N
  error('SIMULATION_STEP: number of densities does not match number of cells.');
end
if size(orqueues, 1) ~= N
  error('SIMULATION_STEP: number of onramp queues does not match number of cells.');
end
if size(flows, 1) ~= N
  error('SIMULATION_STEP: number of mainline flows does not match number of cells.');
end
if size(frflows, 1) ~= N
  error('SIMULATION_STEP: number of offramp flows does not match number of cells.');
end

if sstep < 0
  densities_out = densities;
else
  densities_out = [];
end
flows_out   = [];
frflows_out = [];
speeds_out  = [];
vht_out     = [];
vmt_out     = [];
ploss_out   = [];
tt          = 0;


for i = 1:N
  l = abs(celldata(i).PMend - celldata(i).PMstart);  % i-th cell length

  %%% DENSITY
  if sstep >= 0
     qi = flows(i, 1) + orflows(i, 1);  % total flow entering i-th cell
     if i == N  % last cell
       qo = get_rbflow(densities(N, 1), ...
                       orflows(N, 1), ...
                       flows(N, 1), ...
                       celldata(N), ...
                       outflow, ...
                       ts);
     else
       qo = flows(i+1, 1) + frflows(i, 1);
     end
     d = densities(i, 1) + (ts/l)*(qi - qo);  % new density
     densities_out = [densities_out; max([0 min([d celldata(i).FDrhojam])])];
  end

  %%% MAINLINE FLOW
  % number of vehicles that creates jam in i-th cell
  nJ   = celldata(i).FDrhojam * l;
  % actual number of vehicles in i-th cell
  n    = densities_out(i, 1) * l;
  % number of vehicles coming from on-ramp into i-th cell
  r    = celldata(i).ORflow * ts;
  % on-ramp flow blending coefficient in i-th cell
  g    = celldata(i).ORgamma;
  % congestion wave speed in i-th cell
  w    = celldata(i).FDfmax / (celldata(i).FDrhojam - celldata(i).FDrhocrit);
  % max flow that i-th cell allows to enter when congested
  val2 = max([0 ((w/l) * (nJ - (n + g*r)))]);
  val2 = min([val2 celldata(i).FDfmax]);
  if i == 1
    % max flow that can enter i-th cell 
    val1 = celldata(i).FDfmax;
    % demand on the left boundary of freeway section
    val3 = inflow;
  else
    % max flow that can enter i-th cell 
    val1 = celldata(i-1).FDfmax;
    % (i-1)th cell length
    ll   = abs(celldata(i-1).PMend - celldata(i-1).PMstart);
    % actual number of vehicles in (i-1)th cell
    n    = densities_out(i-1, 1) * ll;
    % number of vehicles coming from on-ramp into (i-1)th cell
    r    = celldata(i-1).ORflow * ts;
    % on-ramp flow blending coefficient in (i-1)th cell
    g    = celldata(i-1).ORgamma;
    % free-flow speed in (i-1)th cell
    v    = celldata(i-1).FDfmax / celldata(i-1).FDrhocrit;
    % off-ramp split ratio in (i-1)th cell
    b    = celldata(i-1).FRbeta;
    % max flow that can be generated by (i-1)th cell
    val3 = (1 - b) * (v/ll) * (n + g*r);
  end
  fl = min([val1 val2 val3]);
  % consistency check: main line flow cannot exceed ((1-b)/b)*FRfmax
  if i > 1
    b = celldata(i-1).FRbeta;
    if (b > 0) & (b < 1)
      if ((b/(1-b)) * fl) > celldata(i-1).FRfmax
        fl = ((1-b)/b) * celldata(i-1).FRfmax; 
      end
    elseif b > 0
      fl = 0;
    end
  end
  flows_out = [flows_out; fl];

  %%% OFFRAMP FLOW
  if i > 1
    if i == N
      if celldata(i-1).FRbeta < 1
        frf = (celldata(i-1).FRbeta / (1 - celldata(i-1).FRbeta)) * flows_out(i, 1);
      else
        v   = celldata(i-1).FDfmax / celldata(i-1).FDrhocrit;  % free-flow speed
        n   = densities_out(i-1, 1) * ll;  % number of vehicles
        r   = celldata(i-1).ORgamma * celldata(i-1).ORflow * ts;  % number of vehicles from on-ramp
        frf = (v/ll) * (n + r);
        frf = min([frf celldata(i-1).FRfmax]);
      end
      frflows_out = [frflows_out; frf];
      qo          = get_rbflow(densities_out(N, 1), ...
                               celldata(N).ORflow, ...
                               flows_out(N, 1), ...
                               celldata(N), ...
                               outflow, ...
                               ts);
      frf         = celldata(N).FRbeta * qo;
      frf         = min([frf celldata(N).FRfmax]);
    else
      if celldata(i-1).FRbeta < 1
        frf = (celldata(i-1).FRbeta / (1 - celldata(i-1).FRbeta)) * flows_out(i, 1);
      else
        v   = celldata(i-1).FDfmax / celldata(i-1).FDrhocrit;  % free-flow speed
        n   = densities_out(i-1, 1) * ll;  % number of vehicles
        r   = celldata(i-1).ORgamma * celldata(i-1).ORflow * ts;  % number of vehicles from on-ramp
        frf = (v/ll) * (n + r);
        frf = min([frf celldata(i-1).FRfmax]);
      end
    end
    frflows_out = [frflows_out; frf];

    %%% SPEED
    if i == N
      if densities_out(i-1, 1) <= 1e-7
        v = celldata(i-1).FDfmax / celldata(i).FDrhocrit;
      else
        v = (flows_out(i, 1) + frflows_out(i-1, 1)) / densities_out(i-1, 1);
      end
      speeds_out = [speeds_out; min([v (celldata(i-1).FDfmax/celldata(i-1).FDrhocrit)])];
      if densities_out(i, 1) <= 1e-7
        v = celldata(i).FDfmax / celldata(i).FDrhocrit;
      else
        v = qo / densities_out(i, 1);
      end
    else
      if densities_out(i-1, 1) <= 1e-7
        v = celldata(i-1).FDfmax / celldata(i).FDrhocrit;
      else
        v = (flows_out(i, 1) + frflows_out(i-1, 1)) / densities_out(i-1, 1);
      end
      v = min([v (celldata(i-1).FDfmax/celldata(i-1).FDrhocrit)]);
    end
    speeds_out = [speeds_out; v];

    %%% VMT & TRAVEL TIME
    if i == N
      vmt_out = [vmt_out; (densities_out(i-1, 1) * ll * ts * speeds_out(i-1, 1))];
      vmt_out = [vmt_out; (densities_out(i, 1) * l * ts * speeds_out(i, 1))];
      tt = tt  +  60 * (ll / speeds_out(i-1, 1));
      tt = tt  +  60 * (l / speeds_out(i, 1));
    else
      vmt_out = [vmt_out; (densities_out(i-1, 1) * ll * ts * speeds_out(i-1, 1))];
      tt = tt  +  60 * (ll / speeds_out(i-1, 1));
    end
  end

  %%% VHT
  vht_out = [vht_out; ((densities_out(i, 1) * l  +  orqueues(i, 1)) * ts)];

  %%% PRODUCTIVITY LOSS
  if densities_out(i, 1) > celldata(i).FDrhocrit
    pls = l  *  (1  -  (flows_out(i, 1) / celldata(i).FDfmax)) * ts;
  else
    pls = 0;
  end
  ploss_out = [ploss_out; pls];

end


simdata = [densities_out flows_out frflows_out speeds_out vht_out vmt_out ploss_out];

return;

Contact us at files@mathworks.com