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;