From 58ed0fd1a4130fb8adb29ae125f26170b31dbaa9 Mon Sep 17 00:00:00 2001
From: Janek Laudan <laudan@tu-berlin.de>
Date: Thu, 31 Oct 2024 14:32:14 +0100
Subject: [PATCH] Add tests for storage capacity under queue and kinematicWaves
 traffic dynamics (#3533)

Co-authored-by: janekdererste <janek.laudan@mailbox.org>
---
 .../qsim/qnetsimengine/QueueWithBuffer.java   |  55 ++---
 .../qnetsimengine/QueueWithBufferTest.java    | 199 ++++++++++++++++++
 2 files changed, 228 insertions(+), 26 deletions(-)
 create mode 100644 matsim/src/test/java/org/matsim/core/mobsim/qsim/qnetsimengine/QueueWithBufferTest.java

diff --git a/matsim/src/main/java/org/matsim/core/mobsim/qsim/qnetsimengine/QueueWithBuffer.java b/matsim/src/main/java/org/matsim/core/mobsim/qsim/qnetsimengine/QueueWithBuffer.java
index 741763e4c87..947b2faa2a6 100644
--- a/matsim/src/main/java/org/matsim/core/mobsim/qsim/qnetsimengine/QueueWithBuffer.java
+++ b/matsim/src/main/java/org/matsim/core/mobsim/qsim/qnetsimengine/QueueWithBuffer.java
@@ -40,7 +40,6 @@
 import org.matsim.core.config.groups.QSimConfigGroup.LinkDynamics;
 import org.matsim.core.config.groups.QSimConfigGroup.TrafficDynamics;
 import org.matsim.core.gbl.Gbl;
-import org.matsim.core.gbl.MatsimRandom;
 import org.matsim.core.mobsim.framework.MobsimDriverAgent;
 import org.matsim.core.mobsim.qsim.interfaces.MobsimVehicle;
 import org.matsim.core.mobsim.qsim.interfaces.SignalGroupState;
@@ -63,8 +62,8 @@
  * Separating out the "lane" functionality from the "link" functionality.
  * <p></p>
  * Design thoughts:<ul>
- * <li> In fast capacity update, the flows are not accumulated in every time step, 
- * rather updated only if an agent wants to enter the link or an agent is added to buffer. 
+ * <li> In fast capacity update, the flows are not accumulated in every time step,
+ * rather updated only if an agent wants to enter the link or an agent is added to buffer.
  * Improvement of 15-20% in the computational performance is observed. amit feb'16
  * (I seem to recall that in the end that statement was not consistently correct.  kai, feb'18)</li>
  * <li>Currently (feb'18), the design is such that (possibly time-dep) flowCap and nEffectiveLanes are "pushed" into the
@@ -90,7 +89,7 @@ public final void addFromWait(final QVehicle veh) {
 
         addToBuffer(veh);
 	}
-	
+
 	/**
 	 * Stores the accumulated fractional parts of the flow capacity. See also
 	 * flowCapFraction.
@@ -189,7 +188,7 @@ private void addValue(double value1, double now) {
 	private double effectiveNumberOfLanesUsedInQsim = Double.POSITIVE_INFINITY ;
 
 	private double accumulatedInflowCap = 1. ;
-	
+
 	private final FlowEfficiencyCalculator flowEfficiencyCalculator;
 
 	private QueueWithBuffer(AbstractQLink.QLinkInternalInterface qlink, final VehicleQ<QVehicle> vehicleQueue, Id<Lane> laneId,
@@ -198,7 +197,7 @@ private QueueWithBuffer(AbstractQLink.QLinkInternalInterface qlink, final Vehicl
 		// the general idea is to give this object no longer access to "everything".  Objects get back pointers (here qlink), but they
 		// do not present the back pointer to the outside.  In consequence, this object can go up to qlink, but not any further. kai, mar'16
 		// Now I am even trying to get rid of the full qLink back pointer (since it allows, e.g., going back to Link). kai, feb'18
-		
+
 //		log.setLevel(Level.DEBUG);
 
 		this.flowEfficiencyCalculator = flowEfficiencyCalculator;
@@ -328,7 +327,7 @@ private void updateFastFlowAccumulation(){
 		double now = context.getSimTimer().getTimeOfDay() ;
 
 		double remainingFlowCapThisTimeStep = subtractConsumptionOfVehiclesThatAreAlreadyInTheBuffer();
-		
+
 		if( this.flowcap_accumulate.getTimeStep() < now
 				&& this.flowcap_accumulate.getValue() < remainingFlowCapThisTimeStep){
 
@@ -344,7 +343,7 @@ private void updateFastFlowAccumulation(){
 
 	private void updateSlowFlowAccumulation(){
 		double remainingFlowCapThisTimeStep = subtractConsumptionOfVehiclesThatAreAlreadyInTheBuffer();
-		
+
 		if (this.thisTimeStepGreen
 				&& this.flowcap_accumulate.getValue() < remainingFlowCapThisTimeStep){
 			double newFlowCap = Math.min(flowcap_accumulate.getValue() + flowCapacityPerTimeStep,
@@ -372,15 +371,15 @@ public final void initBeforeSimStep() {
 	private void calculateFlowCapacity() {
 		// the following is not looking at time because it simply assumes that the lookups are "now". kai, feb'18
 		// I am currently not sure if this statement is correct. kai, feb'18
-		
+
 		// we need the flow capacity per sim-tick and multiplied with flowCapFactor
 		flowCapacityPerTimeStep = unscaledFlowCapacity_s * context.qsimConfig.getTimeStepSize() * context.qsimConfig.getFlowCapFactor() ;
 		inverseFlowCapacityPerTimeStep = 1.0 / flowCapacityPerTimeStep;
-		
+
 		// start with the base assumption, might be adjusted below depending on the traffic dynamics
 		this.effectiveNumberOfLanesUsedInQsim = this.effectiveNumberOfLanes;
 		this.maxInflowUsedInQsim = this.flowCapacityPerTimeStep;
-		
+
 		switch (context.qsimConfig.getTrafficDynamics()) {
 			case queue:
 			case withHoles:
@@ -391,10 +390,10 @@ private void calculateFlowCapacity() {
 				// equal: rho * (vmax + vhole) = vhole * rhojam
 				// rho(qmax) = vhole * rhojam / (vmax + vhole)
 				// qmax = vmax * rho(qmax) = rhojam / (1/vhole + 1/vmax) ;
-				
+
 				// yyyyyy this should possibly be getFreespeed(now). But if that's the case, then maxFlowFromFdiag would
 				// also have to be re-computed with each freespeed change. kai, feb'18
-				
+
 				final double maxFlowFromFdiag = (this.effectiveNumberOfLanes/context.effectiveCellSize) / ( 1./(HOLE_SPEED_KM_H/3.6) + 1/this.qLinkInternalInterface.getFreespeed() ) ;
 				final double minimumNumberOfLanesFromFdiag = this.flowCapacityPerTimeStep * context.effectiveCellSize * ( 1./(HOLE_SPEED_KM_H/3.6) + 1/this.qLinkInternalInterface.getFreespeed() );
 
@@ -455,19 +454,19 @@ private void calculateFlowCapacity() {
 					}
 				}
 				break;
-			
+
 			default: throw new RuntimeException("The traffic dynamics "+context.qsimConfig.getTrafficDynamics()+" is not implemented yet.");
 		}
 //		log.debug( "linkId=" + this.qLink.getLink().getId() + "; flowCapPerTimeStep=" + flowCapacityPerTimeStep +
 //						   "; invFlowCapPerTimeStep=" + inverseFlowCapacityPerTimeStep + "; maxFlowFromFdiag=" + maxFlowFromFdiag ) ;
-		
+
 	}
 
 	private void calculateStorageCapacity() {
 		// The following is not adjusted for time-dependence!! kai, apr'16
 		// No, I think that it simply assumes that the lookups are "now". kai, feb'18
 //		double now = context.getSimTimer().getTimeOfDay() ;
-			
+
 		// first guess at storageCapacity:
 		storageCapacity = this.length * this.effectiveNumberOfLanesUsedInQsim / context.effectiveCellSize * context.qsimConfig.getStorageCapFactor() ;
 //		storageCapacity = this.length * this.qLink.getLink().getNumberOfLanes(now) / context.effectiveCellSize * context.qsimConfig.getStorageCapFactor() ;
@@ -489,7 +488,7 @@ private void calculateStorageCapacity() {
 		if (Double.isNaN(freespeedTravelTime)) {
 			throw new IllegalStateException("Double.NaN is not a valid freespeed travel time for a link. Please check the attributes length and freespeed!");
 		}
-		
+
 		//this assumes that vehicles have the flowEfficiencyFactor of 1.0; the actual flow can be different
 		double tempStorageCapacity = freespeedTravelTime * unscaledFlowCapacity_s * context.qsimConfig.getFlowCapFactor();
 		// yy note: freespeedTravelTime may be Inf.  In this case, storageCapacity will also be set to Inf.  This can still be
@@ -505,7 +504,7 @@ private void calculateStorageCapacity() {
 				QueueWithBuffer.spaceCapWarningCount++;
 			}
 			storageCapacity = tempStorageCapacity;
-			
+
 			// write out the modified qsim behavior as link attribute
 			qLinkInternalInterface.getLink().getAttributes().putAttribute("storageCapacityUsedInQsim", storageCapacity );
 		}
@@ -514,7 +513,7 @@ private void calculateStorageCapacity() {
 		 * () uncongested branch is q(rho) = rho * v_max
 		 * () congested branch is q(rho) = (rho - rho_jam) * v_holes
 		 * () rho_maxflow is where these two meet, resulting in rho_maxflow = v_holes * rho_jam / ( v_holes + v_max )
-		 * () max flow is q(rho_maxflow), resulting in v_max * v_holes * rho_jam / ( v_holes + v_max ) 
+		 * () max flow is q(rho_maxflow), resulting in v_max * v_holes * rho_jam / ( v_holes + v_max )
 		 * () Since everything else is given, rho_jam needs to be large enough so that q(rho_maxflow) can reach capacity, resulting in
 		 *    rho_jam >= capacity * (v_holes + v_max) / (v_max * v_holes) ;
 		 * () In consequence, storage capacity needs to be larger than curved_length * rho_jam .
@@ -554,7 +553,7 @@ private void calculateStorageCapacity() {
 	}
 
 	private double getBufferStorageCapacity() {
-		return flowCapacityPerTimeStep;//this assumes that vehicles have the flowEfficiencyFactor of 1.0 
+		return flowCapacityPerTimeStep;//this assumes that vehicles have the flowEfficiencyFactor of 1.0
 	}
 
 	@Override
@@ -664,7 +663,6 @@ private void removeVehicleFromQueue(final QVehicle veh2Remove) {
 			case withHoles:
 			case kinematicWaves:
 				QueueWithBuffer.Hole hole = new QueueWithBuffer.Hole() ;
-				double ttimeOfHoles = length*3600./HOLE_SPEED_KM_H/1000. ;
 
 				//			double offset = this.storageCapacity/this.flowCapacityPerTimeStep ;
 			/* NOTE: Start with completely full link, i.e. N_storageCap cells filled.  Now make light at end of link green, discharge with
@@ -680,7 +678,12 @@ private void removeVehicleFromQueue(final QVehicle veh2Remove) {
 				//			double nLanes = 2. * flowCapacityPerTimeStep ; // pseudo-lanes
 				//			double ttimeOfHoles = 0.1 * this.storageCapacity/this.flowCapacityPerTimeStep/nLanes ;
 
-				hole.setEarliestLinkExitTime( now + 1.0*ttimeOfHoles + 0.0*MatsimRandom.getRandom().nextDouble()*ttimeOfHoles ) ;
+				// The calculation of the earliest exit time looked like the formula below. It looks like someone tried to include some randomness,
+				// but the random part was multiplied with zero, therefore I removed it. Janek oct' 24
+				// now + 1.0*ttimeOfHoles + 0.0*MatsimRandom.getRandom().nextDouble()*ttimeOfHoles
+				var holeTravelTime = length * 3.6 / HOLE_SPEED_KM_H;
+				var earliestExitTime = now + holeTravelTime;
+				hole.setEarliestLinkExitTime(earliestExitTime) ;
 				hole.setSizeInEquivalents(veh2Remove.getSizeInEquivalents());
 				holes.add( hole ) ;
 				break;
@@ -747,9 +750,9 @@ public void recalcTimeVariantAttributes() {
 		// not speed, since that is looked up anyways.
 		// yy might also make flow and storage self-detecting changes (not really that
 		// much more expensive). kai, feb'18
-		
+
 //		log.debug("just entered recalcTimeVariantAttributes; now=" + this.context.getSimTimer().getTimeOfDay() ) ;
-		
+
 		calculateFlowCapacity();
 		calculateStorageCapacity();
 		flowcap_accumulate.setValue(flowCapacityPerTimeStep);
@@ -916,7 +919,7 @@ public final double getLastMovementTimeOfFirstVehicle() {
 	public final void addTransitSlightlyUpstreamOfStop( final QVehicle veh) {
 		this.vehQueue.addFirst(veh) ;
 	}
-	
+
 	@Override
 	public final void setSignalized( final boolean isSignalized) {
 		qSignalizedItem  = new DefaultSignalizeableItem( qLinkInternalInterface.getToNode().getOutLinks().keySet());
@@ -1020,7 +1023,7 @@ void setVisInfo(Coord upstreamCoord, Coord downstreamCoord) {
 			this.downstreamCoord = downstreamCoord;
 		}
 	}
-	
+
 	private int noOfSeepModeBringFwd = 0;
 
 	private QVehicle peekFromVehQueue(){
diff --git a/matsim/src/test/java/org/matsim/core/mobsim/qsim/qnetsimengine/QueueWithBufferTest.java b/matsim/src/test/java/org/matsim/core/mobsim/qsim/qnetsimengine/QueueWithBufferTest.java
new file mode 100644
index 00000000000..244f51b1246
--- /dev/null
+++ b/matsim/src/test/java/org/matsim/core/mobsim/qsim/qnetsimengine/QueueWithBufferTest.java
@@ -0,0 +1,199 @@
+package org.matsim.core.mobsim.qsim.qnetsimengine;
+
+import org.junit.jupiter.api.Test;
+import org.matsim.api.core.v01.Coord;
+import org.matsim.api.core.v01.Id;
+import org.matsim.core.api.experimental.events.EventsManager;
+import org.matsim.core.config.Config;
+import org.matsim.core.config.ConfigUtils;
+import org.matsim.core.config.groups.QSimConfigGroup;
+import org.matsim.core.mobsim.framework.MobsimDriverAgent;
+import org.matsim.core.mobsim.framework.MobsimTimer;
+import org.matsim.core.mobsim.qsim.interfaces.AgentCounter;
+import org.matsim.core.network.NetworkUtils;
+import org.matsim.vehicles.VehicleType;
+import org.matsim.vehicles.VehicleUtils;
+import org.matsim.vis.snapshotwriters.SnapshotLinkWidthCalculator;
+
+import static org.junit.jupiter.api.Assertions.*;
+import static org.mockito.Mockito.mock;
+
+/**
+ * This suite of unit tests is not exhaustive! I wrote these tests, because I wanted to understand what is happening when the QueueWithBuffer is operated
+ * with config::qsim.trafficDynamics = 'kinematicWaves'.
+ */
+class QueueWithBufferTest {
+
+	@Test
+	void initLinkQueue() {
+
+		var config = ConfigUtils.createConfig();
+		var context = createNetsimeEngineContext(config, new MobsimTimer());
+		var link = createQLink(10, 2, 5400, context, config);
+
+		assertEquals(1, link.getOfferingQLanes().size());
+		assertEquals(link.getAcceptingQLane(), link.getOfferingQLanes().getFirst());
+		QueueWithBuffer qwb = (QueueWithBuffer) link.getAcceptingQLane();
+		assertEquals(1.5, qwb.getSimulatedFlowCapacityPerTimeStep());
+		assertEquals(2 * 2, qwb.getStorageCapacity());
+	}
+
+	@Test
+	void initLinkKinematicWaves() {
+
+		var config = ConfigUtils.createConfig();
+		config.qsim().setTrafficDynamics(QSimConfigGroup.TrafficDynamics.kinematicWaves);
+		var context = createNetsimeEngineContext(config, new MobsimTimer());
+		var link = createQLink(10, 2, 5400, context, config);
+
+		assertEquals(1, link.getOfferingQLanes().size());
+		assertEquals(link.getAcceptingQLane(), link.getOfferingQLanes().getFirst());
+		QueueWithBuffer qwb = (QueueWithBuffer) link.getAcceptingQLane();
+		assertEquals(1.5, qwb.getSimulatedFlowCapacityPerTimeStep());
+
+		var vHole = 15 / 3.6;
+		// according to https://doi.org/10.1080/23249935.2017.1364802 equation 7
+		var expectedStorageCap = link.getSimulatedFlowCapacityPerTimeStep() * link.getLink().getLength() * ( 1 / link.getLink().getFreespeed() + 1 / vHole);
+		assertEquals(expectedStorageCap, qwb.getStorageCapacity(), 0.001);
+	}
+
+	@Test
+	void freeStorageQueue() {
+		var config = ConfigUtils.createConfig();
+		var timer = new MobsimTimer();
+		var context = createNetsimeEngineContext(config, timer);
+		var link = createQLink(10, 1, 1800, context, config);
+		var driver = mock(MobsimDriverAgent.class);
+		var vehicle1 = createVehicle("vehicle-1", driver, 10,1);
+		var vehicle2 = createVehicle("vehicle-2", driver, 10,1);
+
+		// the link should accept vehicles until the storage capacity is exhausted
+		assertTrue(link.getAcceptingQLane().isAcceptingFromUpstream());
+		link.getAcceptingQLane().addFromUpstream(vehicle1);
+		assertTrue(link.getAcceptingQLane().isAcceptingFromUpstream());
+		link.getAcceptingQLane().addFromUpstream(vehicle2);
+		assertFalse(link.getAcceptingQLane().isAcceptingFromUpstream());
+
+		// this should move one vehicle into the buffer, and free one pcu in the queue immediately
+		timer.setTime(1);
+		link.doSimStep();
+
+		// test that capacity is available at the upstream end of the queue
+		assertTrue(link.getAcceptingQLane().isAcceptingFromUpstream());
+
+		// test that only one vehicle was moved into the buffer. Remove the first vehicle, and then the link should not offer more vehicles
+		assertEquals(vehicle1.getId(), link.getOfferingQLanes().getFirst().popFirstVehicle().getId());
+		// this double negation is so terrible!
+		assertTrue(link.getOfferingQLanes().getFirst().isNotOfferingVehicle());
+	}
+
+	@Test
+	void freeStorageKinematicWaves() {
+		var config = ConfigUtils.createConfig();
+		config.qsim().setTrafficDynamics(QSimConfigGroup.TrafficDynamics.kinematicWaves);
+		var timer = new MobsimTimer();
+		var context = createNetsimeEngineContext(config, timer);
+		var link = createQLink(15, 1, 1800, context, config);
+		var driver = mock(MobsimDriverAgent.class);
+		var vehicle1 = createVehicle("vehicle-1", driver, 10,2);
+		var vehicle2 = createVehicle("vehicle-2", driver, 1, 1);
+		link.doSimStep();
+
+		// the link should accept vehicles according to its max inflow capacity
+		// vehicle consumes 2pcu. Max inflow should be: 1/cellSize / (1/vHole + 1/vMax) = 0.588...
+		// the queue with buffer seems to increase the accumulated inflow per 'doSimStep' regardless of the time between invocations. This works, because
+		// in the simulation 'doSimStep' is invoked every timestep. We need 3 'doSimStep's to increase the accumulated inflow to a value over 1.
+		assertTrue(link.getAcceptingQLane().isAcceptingFromUpstream());
+		link.getAcceptingQLane().addFromUpstream(vehicle1);
+		// acc inflow is -1.422
+		assertFalse(link.getAcceptingQLane().isAcceptingFromUpstream());
+		link.doSimStep();
+		// acc inflow is -0.842
+		assertFalse(link.getAcceptingQLane().isAcceptingFromUpstream());
+		link.doSimStep();
+		// acc inflow is -0.236
+		assertFalse(link.getAcceptingQLane().isAcceptingFromUpstream());
+		link.doSimStep();
+		// acc inflow is +0.352 > 0
+		assertTrue(link.getAcceptingQLane().isAcceptingFromUpstream());
+		link.getAcceptingQLane().addFromUpstream(vehicle2);
+		assertFalse(link.getAcceptingQLane().isAcceptingFromUpstream());
+
+		// this should move one vehicle into the buffer, and start a backwards travelling hole
+		timer.setTime(1);
+		link.doSimStep();
+
+		// now, one vehicle should be in the buffer and one vehicle should be in the queue.
+		// The link has free storage capacity, but it is not freed yet, because the leaving vehicle has sent a backwards travelling hole on its way.
+		// the earliest exit time of that hole should be: now + length / vHole -> 1 + 15 * 3.6 / 15km/h = 4.6s
+		assertFalse(link.getAcceptingQLane().isAcceptingFromUpstream());
+		// remove the first vehicle from the buffer and asser that it was the only vehicle in the buffer
+		assertEquals(vehicle1.getId(), link.getOfferingQLanes().getFirst().popFirstVehicle().getId());
+		assertTrue(link.getOfferingQLanes().getFirst().isNotOfferingVehicle());
+
+		// pretend we are doing 5 sim steps. I think we need to do this, as the inflow capacity accumulates per 'doSimStep' and does not keep track
+		// of the last update time.
+		timer.setTime(2);
+		link.doSimStep();
+		assertFalse(link.getAcceptingQLane().isAcceptingFromUpstream());
+		assertTrue(link.getOfferingQLanes().getFirst().isNotOfferingVehicle());
+		timer.setTime(3);
+		link.doSimStep();
+		assertFalse(link.getAcceptingQLane().isAcceptingFromUpstream());
+		assertTrue(link.getOfferingQLanes().getFirst().isNotOfferingVehicle());
+		timer.setTime(4);
+		link.doSimStep();
+		assertFalse(link.getAcceptingQLane().isAcceptingFromUpstream());
+		assertTrue(link.getOfferingQLanes().getFirst().isNotOfferingVehicle());
+		timer.setTime(5); // 5 > 4.6: 4.6 is the arrival time of the backwards travelling hole.
+		link.doSimStep();
+		assertTrue(link.getAcceptingQLane().isAcceptingFromUpstream());
+		assertTrue(link.getOfferingQLanes().getFirst().isNotOfferingVehicle());
+	}
+
+	private static NetsimEngineContext createNetsimeEngineContext(Config config, MobsimTimer timer) {
+		return new NetsimEngineContext(
+			mock(EventsManager.class),
+			5,
+			mock(AgentCounter.class),
+			mock(AbstractAgentSnapshotInfoBuilder.class),
+			config.qsim(),
+			timer,
+			mock(SnapshotLinkWidthCalculator.class)
+		);
+	}
+
+	QLinkImpl createQLink(double length, double lanes, double cap, NetsimEngineContext context, Config config) {
+		var net = NetworkUtils.createNetwork();
+		var n1 = net.getFactory().createNode(Id.createNodeId("n1"), new Coord(0, 0));
+		var n2 = net.getFactory().createNode(Id.createNodeId("n2"), new Coord(0, 100));
+		net.addNode(n1);
+		net.addNode(n2);
+		var link = net.getFactory().createLink(Id.createLinkId("test"), n1, n2);
+		link.setCapacity(cap);
+		link.setFreespeed(10);
+		link.setLength(length);
+		link.setNumberOfLanes(lanes);
+
+		var internalInterface = mock(QNetsimEngineI.NetsimInternalInterface.class);
+		QNodeImpl qNode = new QNodeImpl.Builder(internalInterface, context, config.qsim()).build(n2);
+		qNode.setNetElementActivationRegistry(mock(NetElementActivationRegistry.class));
+		var b = new QLinkImpl.Builder(context, internalInterface);
+		b.setLinkSpeedCalculator(new DefaultLinkSpeedCalculator());
+		var l =  b.build(link, qNode);
+		l.setNetElementActivationRegistry(mock(NetElementActivationRegistry.class));
+		return l;
+	}
+
+	QVehicle createVehicle(String id, MobsimDriverAgent driver, double maxV, double pcu) {
+
+		var type = VehicleUtils.createVehicleType(Id.create("type", VehicleType.class));
+		type.setMaximumVelocity(10);
+		type.setPcuEquivalents(pcu);
+		type.setMaximumVelocity(maxV);
+		var vehicle = VehicleUtils.createVehicle(Id.createVehicleId(id), type);
+		var result = new QVehicleImpl(vehicle);
+		result.setDriver(driver);
+		return result;
+	}
+ }