implement departure time choice DUE

Author: toruseo

uxsim/DTAsolvers/DTAsolvers.py

@@ -888,6 +888,364 @@ def plot_link_stats(s):
 
         plt.show()
 
+    def plot_vehicle_stats(s, orig=None, dest=None):
+        ave_TT = []
+        std_TT = []
+        depature_time = []
+        for vehid in s.route_log[0].keys():
+            if (s.W_sol.VEHICLES[vehid].orig.name == orig or orig == None) and (s.W_sol.VEHICLES[vehid].dest.name == dest or dest == None):
+                length = len(s.route_log)
+                ts = [s.cost_log[day][vehid] for day in range(int(length/2), length)]
+                ave_TT.append(np.average(ts))
+                std_TT.append(np.std(ts))
+                depature_time.append(s.W_sol.VEHICLES[vehid].departure_time_in_second)
+
+        plt.figure()
+        orig_ = orig
+        if orig == None:
+            orig_ = "any"
+        dest_ = dest
+        if dest == None:
+            dest_ = "any"
+        plt.title(f"orig: {orig_}, dest: {dest_}")
+        plt.errorbar(x=depature_time, y=ave_TT, yerr=std_TT, 
+                fmt='bx', ecolor="#aaaaff", capsize=0, label="travel time (mean $\pm$ std)")
+        plt.xlabel("departure time of vehicle")
+        plt.ylabel("travel time")
+        plt.legend()
+        plt.show()
+
+
+class SolverDUE_departure_time_choice:
+    def __init__(s, func_World, desired_arrival_time, time_windows, alpha=1, beta=0.3, gamma=0.9, insensitive_ratio=0.333):
+        """
+        Solve quasi Dynamic User Equilibrium (DUE) problem with departure time choice using day-to-day dynamics. WIP
+
+        Parameters
+        ----------
+        func_World : function
+            function that returns a World object with nodes, links, and demand specifications
+        desired_arrival_time : float
+            desired arrival time to the destination in second
+        time_windows : list of float
+            Time windows for departure time. List of the start time of each time window in second. The final element denotes the end of the windows.
+            For example, `time_windows=[0,600,1200,1800]' means there are 3 windows, namely, 0-600 s, 600-1200 s, and 1200-1800 s.
+        alpha : float
+            travel time penalty coefficient per second
+        beta : float
+            early arrival penalty coefficient per second
+        gamma : float
+            late arrival penalty coefficient per second
+        insensitive_ratio : float
+            travelers become insensitive to the cost difference smaller than this ratio (i.e, bounded rationality, epsilon-Nash equilibirum)
+            
+        Notes
+        -----
+            This function computes a near dynamic user equilibirum state as a steady state of day-to-day dynamical routing game.
+
+            This considers departure time choice as well as route choice.
+
+            Specifically, on day `i`, vehicles choose their route based on actual travel time on day `i-1` with the same departure time.
+            If there are shorter travel time route, they will change with probability `swap_prob`.
+            This process is repeated until `max_iter` day.
+            It is expected that this process eventually reach a steady state.
+            Due to the problem complexity, it does not necessarily reach or converge to Nash equilibrium or any other stationary points.
+            However, in the literature, it is argued that the steady state can be considered as a reasonable proxy for Nash equilibrium or dynamic equilibrium state.
+            There are some theoretical background for it; but intuitively speaking, the steady state can be considered as a realistic state that people's rational behavior will reach.
+
+            This method is based on the following literature:
+            Route choice equilibrium:
+            Ishihara, M., & Iryo, T. (2015). Dynamic Traffic Assignment by Markov Chain. Journal of Japan Society of Civil Engineers, Ser. D3 (Infrastructure Planning and Management), 71(5), I_503-I_509. (in Japanese). https://doi.org/10.2208/jscejipm.71.I_503
+            Iryo, T., Urata, J., & Kawase, R. (2024). Traffic Flow Simulator and Travel Demand Simulators for Assessing Congestion on Roads After a Major Earthquake. In APPLICATION OF HIGH-PERFORMANCE COMPUTING TO EARTHQUAKE-RELATED PROBLEMS (pp. 413-447). https://doi.org/10.1142/9781800614635_0007
+            Iryo, T., Watling, D., & Hazelton, M. (2024). Estimating Markov Chain Mixing Times: Convergence Rate Towards Equilibrium of a Stochastic Process Traffic Assignment Model. Transportation Science. https://doi.org/10.1287/trsc.2024.0523
+
+            Departure time choice equilibrium:
+            Satsukawa, K., Wada, K., & Iryo, T. (2024). Stability analysis of a departure time choice problem with atomic vehicle models. Transportation Research Part B: Methodological, 189, 103039. https://doi.org/10.1016/j.trb.2024.103039
+
+        """
+        s.func_World = func_World
+
+        s.desired_arrival_time = desired_arrival_time
+        s.time_windows = time_windows
+        s.alpha = alpha
+        s.beta = beta
+        s.gamma = gamma
+        s.insensitive_ratio = insensitive_ratio
+
+        s.W_sol = None  #final solution
+        s.W_intermid_solution = None    #latest solution in the iterative process. Can be used when an user terminate the solution algorithm
+        s.dfs_link = []
+
+        warnings.warn("DTA solver is experimental and may not work as expected. It is functional but unstable.")
+    
+    def solve(s, max_iter, n_routes_per_od=10, swap_prob=0.05, print_progress=True, callback=None):
+        """
+        Solve quasi Dynamic User Equilibrium (DUE) problem using day-to-day dynamics. WIP
+
+        Parameters
+        ----------
+        max_iter : int
+            maximum number of iterations
+        n_routes_per_od : int
+            number of routes to enumerate for each OD pair
+        swap_prob : float
+            probability of route swap
+        print_progress : bool
+            whether to print the information
+
+        Returns
+        -------
+        W : World
+            World object with quasi DUE solution (if properly converged)
+        
+        Notes
+        -----
+        `self.W_sol` is the final solution. 
+        `self.W_intermid_solution` is a latest solution in the iterative process. Can be used when an user terminate the solution algorithm.
+        """
+        s.start_time = time.time()
+
+        W_orig = s.func_World()
+        if print_progress:
+            W_orig.print_scenario_stats()
+
+        # enumerate routes for each OD pair
+        n_routes_per_od = n_routes_per_od
+
+        dict_od_to_vehid = defaultdict(lambda: [])
+        for key, veh in W_orig.VEHICLES.items():
+            o = veh.orig.name
+            d = veh.dest.name
+            dict_od_to_vehid[o,d].append(key)
+
+        # dict_od_to_routes = {}
+        # for o,d in dict_od_to_vehid.keys():
+        #     routes = enumerate_k_shortest_routes(W_orig, o, d, k=n_routes_per_od)
+        #     dict_od_to_routes[o,d] = routes
+
+        if W_orig.finalized == False:
+            W_orig.finalize_scenario()
+        dict_od_to_routes = enumerate_k_random_routes(W_orig, k=n_routes_per_od)
+
+        if print_progress:
+            print(f"number of OD pairs: {len(dict_od_to_routes.keys())}, number of routes: {sum([len(val) for val in dict_od_to_routes.values()])}, number of departure time windows: {len(s.time_windows)}")
+
+        # day-to-day dynamics
+        s.ttts = []
+        s.n_swaps = []
+        s.potential_swaps = []
+        s.cost_gaps = []
+        s.route_log = []
+        s.dep_t_log = []
+        s.cost_log = []
+        swap_prob = swap_prob
+        max_iter = max_iter
+
+        print("solving DUE...")
+        for i in range(max_iter):
+            W = s.func_World()
+            # if i != max_iter-1:
+            #     W.vehicle_logging_timestep_interval = -1
+
+            if i != 0:
+                for key in W.VEHICLES:
+                    if key in routes_specified:
+                        W.VEHICLES[key].enforce_route(routes_specified[key])
+                    if key in dep_time_specified:
+                        W.VEHICLES[key].departure_time = int(dep_time_specified[key]/W.DELTAT)
+                        W.VEHICLES[key].departure_time_in_second = dep_time_specified[key]
+                        W.VEHICLES[key].log_t_link[0][0] = int(dep_time_specified[key]/W.DELTAT)
+            
+            route_set = defaultdict(lambda: []) #routes[o,d] = [Route, Route, ...]
+            for o,d in dict_od_to_vehid.keys():
+                for r in dict_od_to_routes[o,d]:
+                    route_set[o,d].append(W.defRoute(r))
+
+            # simulation
+            W.exec_simulation()
+
+            # results
+            W.analyzer.print_simple_stats()
+            #W.analyzer.network_average()
+
+            # attach route choice set to W object for later re-use at different solvers like DSO-GA
+            W.dict_od_to_routes = dict_od_to_routes
+            
+            s.W_intermid_solution = W
+
+            s.dfs_link.append(W.analyzer.link_to_pandas())
+
+            # route swap
+            routes_specified = {}
+            dep_time_specified = {}
+            route_actual = {}
+            dep_time_actual = {}
+            cost_actual = {}
+            n_swap = 0
+            total_cost_gap = 0
+            potential_n_swap = 0
+            for key,veh in W.VEHICLES.items():
+                o = veh.orig.name
+                d = veh.dest.name
+                r, ts = veh.traveled_route()
+                dep_time = veh.departure_time*W.DELTAT
+                travel_time = ts[-1]-ts[0]
+                
+                def schedule_cost(arr_t):
+                    if arr_t > 0:   #negative is null (not arrived)
+                        if arr_t < s.desired_arrival_time:
+                            arr_time_penalty = s.beta * (s.desired_arrival_time-arr_t)
+                        else:
+                            arr_time_penalty = s.gamma * (arr_t-s.desired_arrival_time)
+                    else:
+                        arr_time_penalty = (s.beta+s.gamma)*W.TMAX + W.TMAX*s.alpha
+                    return arr_time_penalty
+
+
+                route_actual[key] = [rr.name for rr in r]
+                dep_time_actual[key] = veh.departure_time*W.DELTAT
+                cost_actual[key] = s.alpha*travel_time + schedule_cost(veh.arrival_time*W.DELTAT)
+
+                veh.departure_time_choice_cost = cost_actual[key]
+                veh.departure_time_choice_cost_tt = s.alpha*travel_time
+                veh.departure_time_choice_cost_sc = schedule_cost(veh.arrival_time*W.DELTAT)
+
+                if veh.state != "end":
+                    continue
+
+                flag_changed = False
+                flag_potential_change = False
+                route_changed = None
+                dep_t_changed = None
+
+                #stabilization technique: portion change
+                change_possibility = random.random() < swap_prob
+                #stabilization technique: travelers with earlier departure time is unlikely to change
+                # dep_t_ratio = (dep_time_actual[key]-s.time_windows[0])/(s.time_windows[-1]-s.time_windows[0])
+                # if random.random() > dep_t_ratio+0.1:
+                #     change_possibility = False
+                
+                actual_cost = s.alpha*r.actual_travel_time(ts[0]) + schedule_cost(dep_time_actual[key]+r.actual_travel_time(ts[0]))
+                #print("time:", r.actual_travel_time(ts[0]), "+", dep_time_actual[key], "=", dep_time_actual[key]+r.actual_travel_time(ts[0]), "; cost:", actual_cost)
+                chosen_cost = actual_cost
+                best_cost = actual_cost #現在の選択と最適選択のコスト差cost_gap算出用．均衡に近いとコスト差が小さい
+                
+                for j,t1 in enumerate(s.time_windows[:-1]):
+                    t2 = s.time_windows[j+1] #time window [t1, t2]
+                    alt_t_dep = W.rng.random()*(t2-t1) + t1
+                    
+                    #stabilization technique: departure time change is relatively rare. choose only 1 time slot on average
+                    change_possibility_time = False
+                    if random.random() < 1/len(s.time_windows[:-1]):
+                        change_possibility_time = True
+
+                    for alt_route in route_set[o,d]:
+                        alt_cost = s.alpha*alt_route.actual_travel_time(alt_t_dep) + schedule_cost(alt_t_dep+alt_route.actual_travel_time(alt_t_dep))
+                        
+                        #print("time:", s.alpha*alt_route.actual_travel_time(ts[0]), "+", alt_t_dep, "=", alt_t_dep+alt_route.actual_travel_time(ts[0]), "; cost:", alt_cost)
+                        
+                        if alt_cost < chosen_cost:
+                            if alt_cost < best_cost:
+                                best_cost = alt_cost
+
+                                if actual_cost - alt_cost > actual_cost*s.insensitive_ratio:
+                                    #stabilization technique: bounded rational behavior. insensitve to difference smaller than x%
+                                    flag_potential_change = 1
+                                    
+                                    if change_possibility and change_possibility_time:
+                                        flag_changed = True
+                                        route_changed = alt_route
+                                        dep_t_changed = alt_t_dep
+                                        chosen_cost = alt_cost
+
+                total_cost_gap += actual_cost-best_cost
+                routes_specified[key] = r
+                dep_time_specified[key] = dep_time
+                if flag_potential_change:
+                    potential_n_swap += W.DELTAN
+                if flag_changed:
+                    n_swap += W.DELTAN
+                    routes_specified[key] = route_changed
+                    dep_time_specified[key] = dep_t_changed
+            
+                veh.departure_time_choice_changed = dep_time_specified[key]
+
+            cost_gap_per_vehicle = total_cost_gap/len(W.VEHICLES)
+            if print_progress:
+                print(f' iter {i}: cost gap: {cost_gap_per_vehicle:.1f}, potential change: {potential_n_swap}, change: {n_swap}, total travel time: {W.analyzer.total_travel_time: .1f}, delay ratio: {W.analyzer.average_delay/W.analyzer.average_travel_time: .3f}')
+            
+            if callback:
+                callback(i, W)
+
+            s.route_log.append(route_actual)
+            s.dep_t_log.append(dep_time_actual)
+            s.cost_log.append(cost_actual)
+
+            s.ttts.append(int(W.analyzer.total_travel_time))
+            s.n_swaps.append(n_swap)
+            s.potential_swaps.append(potential_n_swap)
+            s.cost_gaps.append(cost_gap_per_vehicle)
+
+        
+        s.end_time = time.time()
+
+        print("DUE summary:")
+        last_iters = int(max_iter/4)
+        print(f" total travel time: initial {s.ttts[0]:.1f} -> average of last {last_iters} iters {np.average(s.ttts[-last_iters:]):.1f}")
+        print(f" number of potential changes: initial {s.potential_swaps[0]:.1f} -> average of last {last_iters} iters {np.average(s.potential_swaps[-last_iters:]):.1f}")
+        print(f" cost gap: initial {s.cost_gaps[0]:.1f} -> average of last {last_iters} iters {np.average(s.cost_gaps[-last_iters:]):.1f}")
+        print(f" computation time: {s.end_time - s.start_time:.1f} seconds")
+
+        s.W_sol = W
+        return s.W_sol
+
+    def plot_convergence(s):
+        # iteration plot
+        plt.figure(figsize=(6,2))
+        plt.title("total travel time")
+        plt.plot(s.ttts)
+        plt.xlabel("iter")
+
+        plt.figure(figsize=(6,2))
+        plt.title("number of route change")
+        plt.plot(s.n_swaps)
+        plt.ylim(0,None)
+        plt.xlabel("iter")
+
+        plt.figure(figsize=(6,2))
+        plt.title("travel time difference between choosen route and minimum cost route")
+        plt.plot(s.t_gaps)
+        plt.ylim(0,None)
+        plt.xlabel("iter")
+
+        plt.show()
+
+        # plt.figure()
+        # plot_multiple_y(ys=[s.ttts, s.n_swaps, s.potential_swaps, s.total_t_gaps, np.array(s.total_t_gaps)/np.array(s.potential_swaps)], labels=["total travel time", "number of route change", "number of potential route change", "time gap for potential route change", "time gap per potential route change"])
+        # plt.xlabel("iter")
+        # plt.show()
+
+    def plot_link_stats(s):
+        plt.figure()
+        plt.title("traffic volume")
+        for i in range(len(s.dfs_link[0])):
+            vols = [df["traffic_volume"][i] for df in s.dfs_link]
+            plt.plot(vols, label=s.dfs_link[0]["link"][i])
+        plt.xlabel("iteration")
+        plt.ylabel("volume (veh)")
+        plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
+
+        plt.figure()
+        plt.title("average travel time")
+        for i in range(len(s.dfs_link[0])):
+            vols = [df["average_travel_time"][i] for df in s.dfs_link]
+            plt.plot(vols, label=s.dfs_link[0]["link"][i])
+        plt.xlabel("iteration")
+        plt.ylabel("time (s)")
+        plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
+
+        plt.show()
+    
     def plot_vehicle_stats(s, orig=None, dest=None):
         ave_TT = []
         std_TT = []