How to implement lane-change rules for Nagel Schreckenberg model

Question

I am trying to extend the Nagel-Schreckenberg model for traffic flow to include two lanes of traffic.

I have a function that generates information about a random initial road configuration. A road in this case is thought of as two lines of cells side-by-side (representing two lanes) and each cell can be occupied by a car with an integer velocity less than some specified maximum velocity. The road is closed, ie once a car passes the last cell in a lane, it returns to the starting cell in that lane. The model is iterative. At each iteration, each car first decides whether to change lane (depending on lane-changing rules), and then the situation is updated according to the same rules as for the single-lane situation. It is specifically the lane-change step im struggling with.

The number of cars n, the length of the road l and the maximum velocity vmax are arguments for my road-generating function:

twoLanes[n_, l_, vmax_] := Module[{a = {}, b, lengtha},
Do[AppendTo[a, {b = RandomInteger[], If[b > 0, RandomInteger[{0, vmax}], 0], 
If[b > 0, i, 0], i, RandomInteger[{1, 2}]}]; 
If[Length[Cases[Transpose[a][[1]], 1]] == n, Break[]], {i, l}];
lengtha = Length[a];
Do[AppendTo[a, {0, 0, 0, i}], {i, lengtha + 1, l}];
Map[Rest, Cases[a, {x_ /; x == 1, _, _, _, _}]]]

For example

twoLanes[10, 500, 5]

gives the list

{{2, 2, 2, 1}, {5, 4, 4, 2}, {4, 5, 5, 2}, {0, 7, 7, 2}, {5, 11, 11, 
1}, {0, 13, 13, 2}, {2, 15, 15, 1}, {2, 16, 16, 1}, {0, 18, 18, 
2}, {5, 19, 19, 2}, {0, 20, 20, 2}, {5, 21, 21, 2}, {1, 22, 22, 
2}, {3, 24, 24, 2}, {4, 25, 25, 1}, {1, 26, 26, 2}, {2, 28, 28, 
2}, {5, 31, 31, 1}, {5, 32, 32, 2}, {2, 38, 38, 1}}

each element in this list is of the form {car velocity, car label, car position along road, lane number}. The label is for later use to track a specific car's journey.

How would I implement the following lane-changing procedure (here "gap" refers to the number of cells between two cars):

• the gap ahead in the same lane is less than v+1,
• the gap ahead in the other lane is greater than v+1,
• the gap behind in the other lane is greater than vmax.

I have tried a number of different ways of representing the two lane scenario, and my twoLanes function above is my latest attempt.

One lane code

One-lane road generating function:

ll3[n_, l_, vmax_] := Module[{a = {}, b, lengtha},
Do[AppendTo[
a, {b = RandomInteger[], If[b > 0, RandomInteger[{0, vmax}], 0], 
 If[b > 0, i, 0], i}]; 
If[Length[Cases[Transpose[a][[1]], 1]] == n, Break[]], {i, l}];
lengtha = Length[a];
Do[AppendTo[a, {0, 0, 0, i}], {i, lengtha + 1, l}];
Map[Rest, Cases[a, {x_ /; x == 1, _, _, _}]]]

one-lane update rules:

• If the velocity v of the car is lower than vmax , and the distance to the next car ahead is larger than v + 1, the speed is increased by one.
• If a driver at site i sees the next vehicle at site i+j, with j <= v, he reduces speed to j −1.
• The velocity of each vehicle (if greater than zero) is decreased by one with probability p (‘dawdling’).
• Each vehicle is advanced by v sites.

One-lane update functions :

update2[lane_, length_, vmax_, p_] := Module[{newlane},
newlane = lane;
Do[If[(newlane[[i, 1]] < 
   vmax) && (newlane[[i + 1, 3]] - 
    newlane[[i, 3]]) > (newlane[[i, 1]] + 1), 
newlane[[i, 1]] = newlane[[i, 1]] + 1, 
newlane[[i, 1]] = newlane[[i, 1]]], {i, 1, Length[newlane] - 1}];
If[(newlane[[-1, 1]] < 
  vmax) && (newlane[[1, 3]] - newlane[[-1, 3]] + 
   length) > (newlane[[-1, 1]] + 1), 
newlane[[-1, 1]] = newlane[[-1, 1]] + 1, 
newlane[[-1, 1]] = newlane[[-1, 1]]];
Do[
If[(newlane[[i + 1, 3]] - newlane[[i, 3]]) <= newlane[[i, 1]], 
newlane[[i, 1]] = (newlane[[i + 1, 3]] - newlane[[i, 3]]) - 1, 
newlane[[i, 1]] = newlane[[i, 1]]], {i, 1, Length[newlane] - 1}];
If[(newlane[[1, 3]] - newlane[[-1, 3]] + length) < newlane[[-1, 1]],
newlane[[-1, 
 1]] = (newlane[[1, 3]] - newlane[[-1, 3]] + length) - 1, 
newlane[[-1, 1]] = newlane[[-1, 1]]];
Do[
If[newlane[[i, 1]] > 0 && RandomReal[] < p, 
newlane[[i, 1]] = newlane[[i, 1]] - 1, 
newlane[[i, 1]] = newlane[[i, 1]]], {i, 1, Length[newlane] - 1}];
If[newlane[[-1, 1]] > 0 && RandomReal[] < p, 
newlane[[-1, 1]] = newlane[[-1, 1]] - 1, 
newlane[[-1, 1]] = newlane[[-1, 1]]];
Do[
If[(newlane[[i, 3]] + newlane[[i, 1]]) <= length, 
newlane[[i, 3]] = newlane[[i, 3]] + newlane[[i, 1]], 
newlane[[i, 3]] = newlane[[i, 3]] + newlane[[i, 1]] - length], {i,
 1, Length[lane]}];
Sort[newlane, #1[[3]] < #2[[3]] &]]

Answer 1

The code below has been fixed and is now working. It is based on RNSL model (a modified version of Nagel Schreckenberg model that works for the case of two lanes).

Clear@func;
With[{maxbound = 100, vmax = 5, prob = 0.1},
func[list_] := 
Module[{indices, switchLanes, singleLaneRules, templist, 
randomize, velocitiestoRand, indicessorted, vel, pos, templistSR},
(*functions to obtain position/velocity*)
vel[arg_] := arg[[1]];
pos[arg_] := arg[[3]];

(*switch lanes*)
switchLanes[lis_, car_] := 
Module[{temp = lis, current, currentposinList, indexlanecurrent, 
  indexlaneopposite, carAheadSameLane, carAheadDiffLane, 
  carBehindDiffLane,moveToNextLaneCheck, cond, j, checkcondNextLane},

 (*get information about the particular car*)
 current = Flatten@Cases[lis, {_, car, _, _}];
 currentposinList = First @@ Position[lis, current];
 indexlanecurrent = Last@current;
 indexlaneopposite = First@Complement[{1, 2}, {indexlanecurrent}];

 (*check to determine if there is an urgency to shift lane*)
 moveToNextLaneCheck := (pos[carAheadSameLane] - pos[current]) < (vel[current] + 1);

 (*check for determining if lane change can happen *)
 checkcondNextLane[curr_, ahead_, {}] := (pos[ahead] - pos[curr]) > (vel[curr] + 1); 
 checkcondNextLane[curr_, ahead_, behind_] := (pos[ahead] - pos[curr] > vel[curr] + 1) 
&& (pos[curr] - pos[behind] > vmax);

 (* obtain cars ahead in the same lane, 
 cars ahead and behind in the next lane*)
 carAheadSameLane = 
  Cases[lis, p : {___, current, ___, x : {_, _, _, indexlanecurrent}, ___} :> 
     x, {0}] /. {p__Integer} :> p;
 carAheadDiffLane = Cases[lis, {___, current, ___, x : {_, _, _, indexlaneopposite}, ___}
:> x, {0}] /. {p__Integer} :> p;
 carBehindDiffLane = Cases[lis, {___, x : {_, _, _, indexlaneopposite}, ___, 
      current, ___} :> x, {0}] /. {p__Integer} :> p;

 cond := checkcondNextLane[current, carAheadDiffLane, carBehindDiffLane];

 If[carAheadSameLane =!= {} && carAheadDiffLane =!= {},
  If[moveToNextLaneCheck~And~cond,
    temp = ReplacePart[temp, currentposinList -> 
       Join[current[[1 ;; 2]], {pos[current], indexlaneopposite}]]
    ];
  ];
 temp
 ];

(* Single lane rules applied after switching lanes *)
singleLaneRules[temp_, car_] := 
Module[{carAheadSameLane, lis = temp, current, indexlanecurrent, 
  currentposinList, j},
 current = Flatten@Cases[lis, {_, car, _, _}];
 indexlanecurrent = Last@current;
 currentposinList = First @@ Position[lis, current];

 carAheadSameLane = Cases[lis, p : {PatternSequence[___, current, ___, 
        x : {_, _, _, indexlanecurrent}, ___]} :> x, {0}] /. {p__Integer} :> p;

 If[carAheadSameLane =!= {},
  Which[(vel[current] < vmax) ~And~((pos[carAheadSameLane] - 
         pos[current]) > (vel[current] + 1)), 
    lis = ReplacePart[temp, currentposinList -> {vel[current] + 1, car, pos[current], 
indexlanecurrent}],
    (j = pos[carAheadSameLane] - pos[current]) <= vel[current],
    lis = 
     ReplacePart[temp, currentposinList -> {j - 1, car, pos[current],indexlanecurrent}]
    ];
  ];

 If[vel[current] > 0 && (RandomReal[] < prob),
  lis = ReplacePart[temp,currentposinList -> {vel[current] - 1, car, pos[current], 
      indexlanecurrent}]];
 lis
 ];

indices = list[[All, 2]];
templist = Fold[switchLanes, list, indices];

(* fix for out of bounds*)
templist = templist /. {p : PatternSequence[_, _], posi_?(# >= maxbound &),lane_} :>
Join[{p}, twolanes[[lane, 3 ;;]]];
templist = Sort[templist, #1[[3]] <= #2[[3]] &];

(* apply single lane rules to each lane *)
indicessorted = templist[[All, 2]];
templistSR = Fold[singleLaneRules, templist, indicessorted];

(*update positions*)
Apply[{#1, #2, #3 + #1, #4} &, templistSR , {1}]
]
]

To run or initialize the simulation

Block[{g,show}, 
Monitor[Nest[(show = func@#;
g = Graphics[{Thick, Line[{{2, 0}, {105, 0}, {105, 3}, {2, 3}, {2, 0}}],
{Red,Disk[#, 0.2] & /@ show[[All, {3, 4}]]}}, ImageSize -> Full, Background -> LightBlue];
Pause[0.1];
show) &, twolanes, 500], g]
]

You can use the labels together with NestList rather than Nest to get the history of cars if the need arises

Two simulation results are shown below:

the following shows the traffic flow

sometimes you can get a traffic jam

You can also do something like this to visualize trajectories:

Function[x, Cases[list, {_, x, _, _}, {2}][[All, 3]] //ListPlot[#,Joined -> True] &, Listable][indices]

Answer 2

Okay, I set aside enough time to make progress on this. I scrapped CellularAutomaton and went with a more manual approach.

Now as a self-contained function without global variables.

We might start by representing your road as an array. I leave out labels to simplify things.

Let's suppose your input from twoLanes is lanes.

lanes =
 {{4, 4, 4, 2}, {2, 7, 7, 2}, {0, 11, 11, 1}, {0, 15, 15, 1}, {5, 20, 20, 2},
 {1, 22, 22, 2}, {0, 27, 27, 2}, {5, 30, 30, 2}, {2, 36, 36, 1}, {3, 40, 40, 2}}

toSA =
  SparseArray[#, Automatic, -1] & @*
    Cases[{v_, name_, p_, lane_} :> {lane, p} -> v];

sa = toSA[lanes];

MatrixPlot[sa]

In the operation of the function I shall create distance tables for the gaps fore and aft of every car on the road, e.g.

gaps = toGaps[sa];   (* load code below *)

MatrixPlot[#, ImageSize -> 400] & /@ gaps // Column

The simulation itself:

MatrixPlot /@ NestList[cycle, sa, 99] // ListAnimate

Required code:

toGaps[a_?MatrixQ] := 
  Table[
    d[L]
      // Join[Take[#, -Max@a], #] &
      // FoldList[+## #2 &, #] &
      // RotateRight
      // d @ Drop[#, Max@a] &
    , {d, {Identity, Reverse}}
    , {L, 1 - UnitStep@a}
  ]

Attributes[advance] = HoldFirst;

advance[sa_, gaps_, n_][v_ /; v > 0, {L_, C_}] :=
  gaps[[All, {L, 3 - L}, C]] /. {{sb_, ob_}, {sf_, of_}} :> 
    Which[
      sf <= v && of > v && ob >= Max@sa,
        (sa[[L, C]] = -1; sa[[3 - L, Mod[C + v, n, 1]]] = v),
      sf >= v, 
        sa[[L, {C, Mod[C + v, n, 1]}]] = {-1, v},
      sf < v, 
        sa[[L, {C, Mod[C + sf, n, 1]}]] = {-1, sf}
    ]

cycle =
  Module[{sa = #},
    MapIndexed[advance[sa, toGaps @ #, Length @ First @ #], #, {2}];
    sa
  ] &;