13
$\begingroup$

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]] &]]
$\endgroup$
  • $\begingroup$ By the way this is likely related: (46631). And a question on a similar subject: (135829) $\endgroup$ – Mr.Wizard Feb 19 '17 at 11:26
  • $\begingroup$ yes sorry i left out a few details to prevent the question being too long. Each car is advanced v sites at each iteration $\endgroup$ – jlrawden Feb 19 '17 at 11:35
  • $\begingroup$ Is it critical that each car preserve a unique label throughout? I think that might complicate application, which is why I ask. $\endgroup$ – Mr.Wizard Feb 19 '17 at 11:37
  • $\begingroup$ Your update says: It is specifically the lane-change step im struggling with. Does this mean you have working code apart from this rule? It may be easier to modify that than to start over. In any case several ideas come to mind but none of them are without some effort so I'd like to make sure my effort is not wasted. $\endgroup$ – Mr.Wizard Feb 19 '17 at 11:46
  • $\begingroup$ As long as there is a way for me to visualise the situation over time, the unique label is not needed. I didn't have a clue how to visualise the results of the one-lane case, so i just decided to use labels for each car and plot a graph of one car's position over time. $\endgroup$ – jlrawden Feb 19 '17 at 11:49
8
$\begingroup$

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

enter image description here

sometimes you can get a traffic jam

enter image description here

You can also do something like this to visualize trajectories:

Function[x, Cases[list, {_, x, _, _}, {2}][[All, 3]] //ListPlot[#,Joined -> True] &, Listable][indices]
$\endgroup$
  • $\begingroup$ @jlrawden please see the code and run it several times to note for any discrepancy. If you find that certain condition is missing then I can update the code here. $\endgroup$ – Ali Hashmi Feb 20 '17 at 6:38
  • $\begingroup$ Starting with twolanes = twoLanes[10, 500, 5]; I just observed one car going backward and ending up outside the black frame. I don't think that should happen? $\endgroup$ – Mr.Wizard Feb 20 '17 at 7:14
  • $\begingroup$ @Mr.Wizard oops ! this needs a fix. could you identify it in the code? $\endgroup$ – Ali Hashmi Feb 20 '17 at 7:15
  • $\begingroup$ Not now. I'm just waking up and checking messages. I'll come back to this problem later. $\endgroup$ – Mr.Wizard Feb 20 '17 at 7:16
  • 1
    $\begingroup$ @Mr.Wizard the code has now been fixed $\endgroup$ – Ali Hashmi Feb 21 '17 at 22:26
7
$\begingroup$

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]

enter image description here

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

enter image description here

The simulation itself:

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

enter image description here

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
  ] &;
$\endgroup$
  • $\begingroup$ I need to look into your code. looks remarkable ! $\endgroup$ – Ali Hashmi Feb 20 '17 at 6:07
  • 1
    $\begingroup$ @AliHashmi For an example that actually works well please see (46640). I ran into trouble here because these patterns must be quite explicit; using Condition etc. causes problems. This makes handling the v+1 rules difficult. Assuming I continue to work on this I may need to scrap CellularAutomaton entirely. But it looks like you've already produced a much more complete implementation so if anything I may attempt to improve that instead. $\endgroup$ – Mr.Wizard Feb 20 '17 at 6:58
  • $\begingroup$ I am upvoting this ! looks like it is getting closer. btw did you manage to see where the velocity could be getting negative in my code. or any condition that i may have missed. $\endgroup$ – Ali Hashmi Feb 21 '17 at 14:15
  • $\begingroup$ @AliHashmi Thank you. Sorry, but I did not; I chose to work on my own code instead. If I have time I shall try to look more closely through yours, but I often become occupied with a new Question and never return to such things. $\endgroup$ – Mr.Wizard Feb 21 '17 at 17:02
  • $\begingroup$ @jlrawden You're welcome. :-) Please let me know if you need anything explained; this code seems rather opaque to me and I didn't include any comments. $\endgroup$ – Mr.Wizard Feb 22 '17 at 18:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.