Download All

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

ctmBATCH() 
function ctmBATCH()
% CTMBATCH - equivalent of CTMGUI for batch mode, called from CTMSIM.
%
% Last modified:   03/18/2007.

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

global g_ctmGUI;

%%% INITIALIZATION %%%

if ~isstruct(g_ctmGUI) | ~isfield(g_ctmGUI, 'configFile')
  g_ctmGUI.configFile = 'ctm';  % default: 'ctm.mat'
end

g_ctmGUI.numColors = 1024;
g_ctmGUI.isStopped = 0;
g_ctmGUI.isSaved   = 1;

load(g_ctmGUI.configFile);

g_ctmGUI.cellData = celldata;  % MUST (!) be present in the config file

if ~exist('TS', 'var')
  TS = 1;
end
minl        = min(get_cell_lengths(g_ctmGUI.cellData));
maxv        = max(get_ff_speeds(g_ctmGUI.cellData));
g_ctmGUI.TS = min([TS ((minl / maxv) - eps)]);  % sampling period

if ~exist('plotTS', 'var')
  plotTS = 1/60;  % plote once a minute
end
g_ctmGUI.plotPeriod = max([round(plotTS/g_ctmGUI.TS) 1]);

if ~exist('dataFile', 'var')
  dataFile = [g_ctmGUI.configFile '_data'];
end
g_ctmGUI.dataFile = dataFile;

if ~exist('taxislim', 'var')
  taxislim = 100;  % time axis limit
end
g_ctmGUI.taxislim = taxislim;

if ~exist('timeout', 'var')
  timeout = 2;  % in seconds
end
g_ctmGUI.timeout = timeout;

if ~exist('yoColorRatio', 'var')
  yoColorRatio = [0.95 0.95];
end
g_ctmGUI.yoColorRatio = yoColorRatio;

if ~exist('densityCMF', 'var')
  densityCMF = jet(g_ctmGUI.numColors);
end
g_ctmGUI.densityCMF = densityCMF;

if ~exist('flowCMF', 'var')
  flowCMF = jet(g_ctmGUI.numColors);
end
g_ctmGUI.flowCMF = flowCMF;

if ~exist('orflowCMF', 'var')
  orflowCMF = jet(g_ctmGUI.numColors);
end
g_ctmGUI.orflowCMF = orflowCMF;

if ~exist('orqueueCMF', 'var')
  orqueueCMF = jet(g_ctmGUI.numColors);
end
g_ctmGUI.orqueueCMF = orqueueCMF;

if ~exist('frflowCMF', 'var')
  frflowCMF = jet(g_ctmGUI.numColors);
end
g_ctmGUI.frflowCMF = frflowCMF;

if ~exist('frbetaCMF', 'var')
  frbetaCMF = jet(g_ctmGUI.numColors);
end
g_ctmGUI.frbetaCMF = frbetaCMF;

if ~exist('speedCMF', 'var')
  speedCMF = jet(g_ctmGUI.numColors);
end
g_ctmGUI.speedCMF = speedCMF;

if ~exist('vhtCMF', 'var')
  vhtCMF = jet(g_ctmGUI.numColors);
end
g_ctmGUI.vhtCMF = vhtCMF;

if ~exist('vmtCMF', 'var')
  vmtCMF = jet(g_ctmGUI.numColors);
end
g_ctmGUI.vmtCMF = vmtCMF;

if ~exist('initialDensities', 'var')
  initialDensities = get_critical_densities(g_ctmGUI.cellData);
end
g_ctmGUI.densities = initialDensities;

if ~exist('inflow', 'var')
  inflow = g_ctmGUI.cellData(1).FDfmax;
end
g_ctmGUI.inflow = inflow;

if ~exist('ifKnob', 'var')
  ifKnob = 1;
end
g_ctmGUI.ifKnob = ifKnob;

if ~exist('outflow', 'var')
  outflow = g_ctmGUI.cellData(end).FDfmax;
end
g_ctmGUI.outflow = outflow;

if ~exist('ofKnob', 'var')
  ofKnob = 1;
end
g_ctmGUI.ofKnob = ofKnob;

if ~exist('demandProfile', 'var')
  demandProfile = [];
end
g_ctmGUI.demandProfile = demandProfile;

if ~exist('demandTS', 'var')
  demandTS = 1/12;  % 5 minutes
end
g_ctmGUI.demandPeriod = round(demandTS/TS);
g_ctmGUI.demandIndex  = 0;
g_ctmGUI.demandCount  = 0;

if ~exist('betaProfile', 'var')
  betaProfile = [];
end
g_ctmGUI.betaProfile = betaProfile;

if ~exist('betaTS', 'var')
  betaTS = 1/12;  % 5 minutes
end
g_ctmGUI.betaPeriod = round(betaTS/TS);
g_ctmGUI.betaIndex  = 0;
g_ctmGUI.betaCount  = 0;

if ~exist('frflowProfile', 'var')
  frflowProfile = [];
end
g_ctmGUI.frflowProfile = frflowProfile;

if ~exist('frflowTS', 'var')
  frflowTS = 1/12;  % 5 minutes
end
g_ctmGUI.frflowPeriod = round(frflowTS/TS);
g_ctmGUI.frflowIndex  = 0;
g_ctmGUI.frflowCount  = 0;

if ~exist('maxSimStep', 'var')
  maxSimStep = Inf;
end
g_ctmGUI.maxSimStep = maxSimStep;

if ~exist('maxSimTime', 'var')
  maxSimTime = 24;  % hours
end
g_ctmGUI.maxSimTime = maxSimTime;

g_ctmGUI.densityData    = [];
g_ctmGUI.flowData       = [];
g_ctmGUI.speedData      = [];
g_ctmGUI.orflowData     = [];
g_ctmGUI.orqueueData    = [];
g_ctmGUI.frflowData     = [];
g_ctmGUI.frbetaData     = [];
g_ctmGUI.vhtData        = [];
g_ctmGUI.vmtData        = [];
g_ctmGUI.delayData      = [];
g_ctmGUI.plossData      = [];
g_ctmGUI.traveltimeData = [];
g_ctmGUI.tvhData        = [];
g_ctmGUI.tvmData        = [];
g_ctmGUI.cdelayData     = [];
g_ctmGUI.cplossData      = [];
g_ctmGUI.timeMinutes    = [];  % data needed for the TIME axis
g_ctmGUI.timeStep       = 0;
g_ctmGUI.linIndex       = 0;
g_ctmGUI.loutIndex      = 1;

[g_ctmGUI.maxDensity, g_ctmGUI.maxFlow] = get_max_df(g_ctmGUI.cellData);
if exist('maxFlow', 'var')
  g_ctmGUI.maxFlow = maxFlow;
end
if exist('maxDensity', 'var')
  g_ctmGUI.maxDensity = maxDensity;
end
if exist('minFlow', 'var')
  g_ctmGUI.minFlow = minFlow;
else
  g_ctmGUI.minFlow = 0;
end
if exist('minDensity', 'var')
  g_ctmGUI.minDensity = minDensity;
else
  g_ctmGUI.minDensity = 0;
end

if ~exist('freeway', 'var')
  freeway = 'Freeway XYZ';
end
g_ctmGUI.freeway = freeway;

g_ctmGUI.pmData  = get_pm_data(g_ctmGUI.cellData);
g_ctmGUI.frbetas = get_fr_splitratios(g_ctmGUI.cellData)';
g_ctmGUI.orflows = get_or_flows(g_ctmGUI.cellData)';

xmin = g_ctmGUI.pmData(1);
xmax = g_ctmGUI.pmData(end);
if xmin > xmax
  xlim = [xmax xmin];
else
  xlim = [xmin xmax];
end
g_ctmGUI.xLims = xlim;





%%% SIMULATION %%%
NC = size(g_ctmGUI.cellData, 2);  % number of cells

if isempty(g_ctmGUI.timeMinutes)  % start new
  densities                 = g_ctmGUI.densities;
  orflows                   = get_or_fmax(g_ctmGUI.cellData);
  orflows                   = orflows';
  orqueues                  = zeros(NC, 1);
  g_ctmGUI.controllerCounts = zeros(NC, 1);
  g_ctmGUI.timeStep         = -1;
  flows                     = zeros(NC, 1);
  frflows                   = zeros(NC, 1);
else  % continue where we left off
  densities = g_ctmGUI.densityData(1:NC, end);
  flows     = g_ctmGUI.flowData(1:NC, end);
  orflows   = g_ctmGUI.orflowData(1:NC, end);
  orqueues  = g_ctmGUI.orqueueData(1:NC, end);
  frbetas   = g_ctmGUI.frbetaData(1:NC, end);
  frflows   = g_ctmGUI.frflowData(1:NC, end);
  speeds    = g_ctmGUI.speedData(1:NC, end);
  vht       = g_ctmGUI.vhtData(1:NC, end);
  vmt       = g_ctmGUI.vmtData(1:NC, end);
  delay     = g_ctmGUI.delayData(1:NC, end);
  ploss     = g_ctmGUI.plossData(1:NC, end);
end

g_ctmGUI.controllerPeriods = get_controller_periods(g_ctmGUI.cellData, g_ctmGUI.TS);
g_ctmGUI.isStopped         = 0;

orflows_prev = orflows;
onerow       = ones(1, NC);
timeStep     = g_ctmGUI.timeStep;
count        = -1;

% main loop
while (timeStep <= g_ctmGUI.maxSimStep) & ...
        (timeStep * g_ctmGUI.TS <= g_ctmGUI.maxSimTime) & ...
        (g_ctmGUI.isStopped == 0)

  if (count + 1) == g_ctmGUI.plotPeriod
    count = -1;
  end
  count = count + 1;

  % determine actual on-ramp demands, flows, queues
  if ~isempty(g_ctmGUI.demandProfile)
    if g_ctmGUI.demandCount == 0
      g_ctmGUI.demandIndex = g_ctmGUI.demandIndex + 1;
      if g_ctmGUI.demandIndex > size(g_ctmGUI.demandProfile, 1)
        g_ctmGUI.demandIndex = 1;
        g_ctmGUI.isStopped   = 1;
      end
    end
    if (g_ctmGUI.demandCount + 1) == g_ctmGUI.demandPeriod
      g_ctmGUI.demandCount = -1;
    end
    g_ctmGUI.demandCount = g_ctmGUI.demandCount + 1;
    demands              = g_ctmGUI.demandProfile(g_ctmGUI.demandIndex, :);
    g_ctmGUI.inflow      = g_ctmGUI.ifKnob * demands(1);
    g_ctmGUI.outflow     = g_ctmGUI.ofKnob * demands(end);
    demands              = adjust_or_demands(demands(2:(NC+1))', g_ctmGUI.cellData);
    
    if 1
      [orflows, g_ctmGUI.controllerCounts] = ...
         auto_control(densities, ...
                      [demands orflows orqueues], ...
                      g_ctmGUI.controllerCounts, ...
                      g_ctmGUI.controllerPeriods, ...
                      g_ctmGUI.cellData, ...
                      g_ctmGUI.TS, ...
                      1);
    else
      orflows = adjust_or_flow_limits(densities, g_ctmGUI.cellData, g_ctmGUI.TS);
    end
  else
    demands  = get_or_flows(g_ctmGUI.cellData);
    demands  = demands';
    orflows  = demands;
  end

  % make sure on-ramp flows do nor exceed allowed limits
  orflows = adjust_or_flows(densities, orflows, g_ctmGUI.cellData, g_ctmGUI.TS);

  % adjust on-ramp flow and queue values
  [orflows, orqueues] = adjust_or_data(demands, ...
                                       orflows, ...
                                       orqueues, ...
                                       g_ctmGUI.TS);

  % update on-ramp flows in cell structures
  g_ctmGUI.cellData   = set_or_flows(g_ctmGUI.cellData, orflows);
  if timeStep == -1
    orflows_prev = orflows;
  end

  % determine off-ramp split ratio values
  if ~isempty(g_ctmGUI.betaProfile)
    if g_ctmGUI.betaCount == 0
      g_ctmGUI.betaIndex = g_ctmGUI.betaIndex + 1;
      if g_ctmGUI.betaIndex > size(g_ctmGUI.betaProfile, 1)
        g_ctmGUI.betaIndex = 1;
        g_ctmGUI.isStopped = 1;
      end
    end
    if (g_ctmGUI.betaCount + 1) == g_ctmGUI.betaPeriod
      g_ctmGUI.betaCount = -1;
    end
    g_ctmGUI.betaCount = g_ctmGUI.betaCount + 1;
    frbetas            = g_ctmGUI.betaProfile(g_ctmGUI.betaIndex, :);
    frbetas            = adjust_fr_split_ratios(frbetas', g_ctmGUI.cellData);
    g_ctmGUI.cellData  = set_fr_split_ratios(g_ctmGUI.cellData, frbetas);
  end


  if ~isempty(g_ctmGUI.frflowProfile)
    if g_ctmGUI.frflowCount == 0
      g_ctmGUI.frflowIndex = g_ctmGUI.frflowIndex + 1;
      if g_ctmGUI.frflowIndex > size(g_ctmGUI.frflowProfile, 1)
        g_ctmGUI.frflowIndex = 1;
        g_ctmGUI.isStopped = 1;
      end
    end
    if (g_ctmGUI.frflowCount + 1) == g_ctmGUI.frflowPeriod
      g_ctmGUI.frflowCount = -1;
    end
    g_ctmGUI.frflowCount = g_ctmGUI.frflowCount + 1;
    frflows              = g_ctmGUI.frflowProfile(g_ctmGUI.frflowIndex, :);
    frflows              = adjust_fr_flows(densities, ...
                                           frflows', ...
                                           g_ctmGUI.cellData, ...
                                           g_ctmGUI.TS);

    [simdata, traveltime] = simulation_step_2(g_ctmGUI.cellData, ...
                                              densities, ...
                                              orflows_prev, ...
                                              orqueues, ...
                                              flows, ...
                                              frflows, ...
                                              g_ctmGUI.inflow, ...
                                              g_ctmGUI.outflow, ...
                                              g_ctmGUI.TS, ...
                                              timeStep);
    frbetas               = simdata(:, 3);

    % update beta values in the cell structures
    g_ctmGUI.cellData     = set_fr_split_ratios(g_ctmGUI.cellData, frbetas);
  else
    [simdata, traveltime] = simulation_step(g_ctmGUI.cellData, ...
                                            densities, ...
                                            orflows_prev, ...
                                            orqueues, ...
                                            flows, ...
                                            frflows, ...
                                            g_ctmGUI.inflow, ...
                                            g_ctmGUI.outflow, ...
                                            g_ctmGUI.TS, ...
                                            timeStep);
    frflows               = simdata(:, 3);

    % get off-ramp split ratios
    frbetas               = get_fr_splitratios(g_ctmGUI.cellData);
    frbetas               = frbetas';
  end

  densities    = simdata(:, 1);
  flows        = simdata(:, 2);
  speeds       = simdata(:, 4);
  orflows_prev = orflows;

  % compute Vehicle Hours Traveled (VHT), Vehicle Miles Traveled (VMT)
  % and Productivity Loss
  if timeStep >= 0
    vht   = vht + simdata(:, 5);
    vmt   = vmt + simdata(:, 6);
    ploss = ploss + simdata(:, 7);
  else
    vht   = g_ctmGUI.plotPeriod * simdata(:, 5);
    vmt   = g_ctmGUI.plotPeriod * simdata(:, 6);
    ploss = g_ctmGUI.plotPeriod * simdata(:, 7);
  end

  % compute total Productivity Loss
  cploss    = sum(ploss);

  % compute new Total Vehicle Hours
  tvh        = sum(vht);

  % compute new Total Vehicle Miles
  tvm        = sum(vmt);

  % compute delay
  delay      = vht  -  (vmt ./ get_ff_speeds(g_ctmGUI.cellData));
  delay      = max([delay'; zeros(1, NC)])';
  cdelay     = sum(delay);

  % next time step
  timeStep   = timeStep + 1;

  if count > 0
    continue;
  end

  % update data arrays
  g_ctmGUI.timeStep       = timeStep;
  tm                      = g_ctmGUI.timeStep * g_ctmGUI.TS * 60;
  g_ctmGUI.timeMinutes    = [g_ctmGUI.timeMinutes tm];
  g_ctmGUI.densityData    = [g_ctmGUI.densityData [densities; densities(end, 1)]];
  g_ctmGUI.flowData       = [g_ctmGUI.flowData [flows; flows(end, 1)]];
  g_ctmGUI.orflowData     = [g_ctmGUI.orflowData [orflows; orflows(end, 1)]];
  g_ctmGUI.orqueueData    = [g_ctmGUI.orqueueData [orqueues; 200]];
  g_ctmGUI.frflowData     = [g_ctmGUI.frflowData [frflows; frflows(end, 1)]];
  g_ctmGUI.frbetaData     = [g_ctmGUI.frbetaData [frbetas; frbetas(end, 1)]];
  g_ctmGUI.speedData      = [g_ctmGUI.speedData [speeds; speeds(end, 1)]];
  g_ctmGUI.vhtData        = [g_ctmGUI.vhtData [vht; vht(end, 1)]];
  g_ctmGUI.vmtData        = [g_ctmGUI.vmtData [vmt; vmt(end, 1)]];
  g_ctmGUI.delayData      = [g_ctmGUI.delayData [delay; delay(end, 1)]];
  g_ctmGUI.plossData      = [g_ctmGUI.plossData [ploss; ploss(end, 1)]];
  g_ctmGUI.traveltimeData = [g_ctmGUI.traveltimeData traveltime];
  g_ctmGUI.tvhData        = [g_ctmGUI.tvhData tvh];
  g_ctmGUI.tvmData        = [g_ctmGUI.tvmData tvm];
  g_ctmGUI.cdelayData     = [g_ctmGUI.cdelayData cdelay];
  g_ctmGUI.cplossData     = [g_ctmGUI.cplossData cploss];
  vht                     = zeros(NC, 1);
  vmt                     = zeros(NC, 1);
  ploss                   = zeros(NC, 1);

end

l_save_simulation(g_ctmGUI);

return;




% save simulation
function l_save_simulation(ctmGUI_data)
if str2num(version('-release')) >= 14
  feval(@save, '-v6', ctmGUI_data.dataFile, 'ctmGUI_data');
else
  feval(@save, ctmGUI_data.dataFile, 'ctmGUI_data');
end

return;

Contact us at files@mathworks.com