add congestion externality estimation method

Author: toruseo

tests/test_other_functions.py

@@ -1188,4 +1188,76 @@ def test_construct_time_space_network():
     construct_time_space_network(W)
 
     assert W.TSN_paths["4", 0]["7", "end"][-2] == ('7', 340)
-    assert equal_tolerance(W.TSN_costs["4", 0]["7", "end"], W.TSN_paths["4", 0]["7", "end"][-2][1])
+    assert equal_tolerance(W.TSN_costs["4", 0]["7", "end"], W.TSN_paths["4", 0]["7", "end"][-2][1])
+
+@pytest.mark.flaky(reruns=10)
+def test_estimate_congestion_externality_1link():
+    def create_world(seed):
+        W = World(
+            name="aaa",    # Scenario name
+            deltan=10,   # Simulation aggregation unit delta n
+            tmax=2400,  # Total simulation time (s)
+            print_mode=0, save_mode=0, show_mode=0,    # Various options
+            random_seed=seed    # Set the random seed
+        )
+
+        # Define the scenario
+        ## Create nodes
+        W.addNode(name="orig1", x=0, y=0)
+        W.addNode("mid", 1, 1)
+        W.addNode("dest", 2, 1)
+        ## Create links between nodes
+        link1 = W.addLink(name="link1", start_node="orig1", end_node="mid", length=2000, free_flow_speed=20, number_of_lanes=1, capacity_out=0.8)
+        link2 = W.addLink("link3", "mid", "dest", length=1000, free_flow_speed=20, number_of_lanes=1, capacity_in=0.4)
+
+        R1 = uxsim.Route(W, [link1, link2])
+        ## Create OD traffic demand between nodes
+        W.adddemand(orig="orig1", dest="dest", t_start=0, t_end=1200, flow=0.2)
+        W.adddemand(orig="orig1", dest="dest", t_start=300, t_end=600, flow=0.4)
+
+        return W
+
+    dep_ts = []
+    acts = []
+    ests = []
+    for _ in range(50):
+        seed = random.randint(0,999999)
+
+        W = create_world(seed)
+        W.exec_simulation()
+        W.analyzer.print_simple_stats()
+        #W.analyzer.time_space_diagram_traj_links(["link1", "link3"])
+
+        key = random.choice(list(W.VEHICLES.keys()))
+        veh = W.VEHICLES[key]
+        veh_tt = veh.travel_time*W.DELTAN
+        route = veh.traveled_route()[0]
+        dep_t = veh.log_t_link[1][0]
+        ext = estimate_congestion_externality_route(W, route, dep_t)
+        ttt1 = W.analyzer.total_travel_time
+
+
+        W2 = create_world(seed)
+        W2.VEHICLES[key].state="abort"
+        W2.exec_simulation()
+        W2.analyzer.print_simple_stats()
+        #W2.analyzer.time_space_diagram_traj_links(["link1", "link3"])
+        ttt2 = W2.analyzer.total_travel_time
+
+        #print(f"selected vehicle - dep time {veh.departure_time_in_second}")
+        #print("actual ext:", ttt1-ttt2-veh_tt, f"{ttt1}-{ttt2}-{veh_tt}")
+        #print("esti. ext: ", ext)
+
+        dep_ts.append(veh.departure_time_in_second)
+        acts.append(ttt1-ttt2-veh_tt)
+        ests.append(ext)
+
+    # W.analyzer.time_space_diagram_traj_links(["link1", "link3"])
+
+    # figure()
+    # plot(dep_ts, acts, "o", label="actual")
+    # plot(dep_ts, ests, "x", label="estimated")
+    # legend()
+    # grid()
+
+    assert uxsim.equal_tolerance(average(acts), average(ests), rel_tol=0.333)

uxsim/DTAsolvers/ALNS.py

@@ -2,10 +2,8 @@
 Adaptive Large Neighborhood Search (ALNS) for DSO
 
 - Initially coded by GPT-5, refined by human and Claude Code. 
-
-- 破壊/修復オペレータを辞書で差し替え可能な、本格ALNSの最小クリーン実装。
-- 主要API: init_alns / auto_step / run_auto / finalize
-- 拡張: _DESTROY_IMPLS / _REPAIR_IMPLS に関数を追加/置換するだけで新オペレータが自動採用される
+- Customizable minimum implementation
+- Includes domain-knowledge-based heuristics
 """
 
 from __future__ import annotations

uxsim/Utilities/Utilities.py

@@ -9,6 +9,7 @@
 import tqdm
 import random
 from collections import defaultdict
+import warnings
 
 def generate_grid_network(W, imax, jmax, **kwargs):
     """
@@ -419,3 +420,82 @@ def construct_time_space_network(W, dt=None, from_origin_only=True):
     #             print(f"Shortest path from {source[0]} to {target[0]} on t = {source[1]}: cost = {costs[source][target]}: {paths[source][target]}")
 
 
+
+def estimate_congestion_externality_link(W, link, t):
+    """
+    Estimate externality caused by a hypothetical vehicle traveling the link in the world. Subfunction for `estimate_congestion_externality_route`.
+
+    Parameters
+    ----------
+    W : World
+        The World object.
+    link : str | Link
+        Link object or link name.
+    t : float
+        The departure time of the hypothetical vehicle.
+    """
+    try:
+        link = W.get_link(link)
+        ts = int(t/W.DELTAT)
+
+        #同じ待ち行列で自分より後ろで待っている車両台数の算出
+        ts_end = ts
+        for i in range(ts, len(link.traveltime_actual), 1):
+            if link.traveltime_actual[i] > link.length/link.free_flow_speed*1.01 and link.cum_arrival[i+int(link.traveltime_actual[i]/W.DELTAT)]-link.cum_departure[i+int(link.traveltime_actual[i]/W.DELTAT)] > 0:
+            #if link.traveltime_actual[i] > link.length/link.free_flow_speed*1.01:
+                ts_end = i
+            else:
+                break
+        
+        veh_count = link.cum_arrival[ts_end]-link.cum_arrival[ts]
+
+        #局所的容量の算出
+        #TODO:微妙に系統的にずれる
+        #TODO:スピルオーバー時の除外 
+        #TODO:信号待ちの一台目のときに破綻してしまう．頑健で賢い方法あるか？信号，このリンクの容量制約，下流側リンク・ノード容量制約を全部考えられるもの
+        ts_out = int(ts+int(link.traveltime_actual[ts]/W.DELTAN))+1
+        ts_out_leader = ts_out
+        for i in range(ts_out, 0, -1):
+            if link.cum_departure[i] < link.cum_departure[ts_out]:
+                ts_out_leader = i
+                break
+
+        headway_consumed = (ts_out-ts_out_leader)*W.DELTAT 
+
+        ext = headway_consumed*veh_count
+        
+        return ext
+    
+    except Exception as e:
+        #細かいインデックスエラーなどを一時的に回避 TODO: fix
+
+        warnings.warn(f"Error at `estimate_congestion_externality_link`: {e}", UserWarning)
+
+        return 0
+
+
+
+def estimate_congestion_externality_route(W, route, t):
+    """
+    Estimate externality caused by a hypothetical vehicle traveling the route in the world. 
+    
+    Specifically, it estimates the total travel time that would be achieved if the vehicle did not exist and then calculates the difference. This estimation is based on a simple queuing model for each link of the route that the vehicle passes through. Does not fully capture queue spilovers.
+
+    Parameters
+    ----------
+    W : World
+        The World object.  
+    route : Route
+        Route object.
+    t : float
+        The departure time of the hypothetical vehicle.
+    """
+
+    tts = route.actual_travel_time(t, return_details=True)[1]
+    exts = 0
+    for i,link in enumerate(route):
+        exts += estimate_congestion_externality_link(W, link, t)
+        t += tts[i]
+    
+    return exts
+

uxsim/analyzer.py

@@ -325,6 +325,7 @@ def print_simple_stats(s, force_print=False):
         s.W.print(" number of completed trips:\t", s.trip_completed, "/", s.trip_all)
         #s.W.print(" number of completed trips:\t", s.trip_completed, "/", len(s.W.VEHICLES)*s.W.DELTAN)
         if s.trip_completed > 0:
+            s.W.print(f" total travel time:\t\t {s.total_travel_time:.1f} s")
             s.W.print(f" average travel time of trips:\t {s.average_travel_time:.1f} s")
             s.W.print(f" average delay of trips:\t {s.average_delay:.1f} s")
             s.W.print(f" delay ratio:\t\t\t {s.average_delay/s.average_travel_time:.3f}")
@@ -336,6 +337,7 @@ def print_simple_stats(s, force_print=False):
             print(" number of completed trips:\t", s.trip_completed, "/", s.trip_all)
             #print(" number of completed trips:\t", s.trip_completed, "/", len(s.W.VEHICLES)*s.W.DELTAN)
             if s.trip_completed > 0:
+                print(f" total travel time:\t\t {s.total_travel_time:.1f} s")
                 print(f" average travel time of trips:\t {s.average_travel_time:.1f} s")
                 print(f" average delay of trips:\t {s.average_delay:.1f} s")
                 print(f" delay ratio:\t\t\t {s.average_delay/s.average_travel_time:.3f}")