diff --git a/uam_system_model/Pricing.py b/uam_system_model/Pricing.py index 2ea9da3..5e84fad 100644 --- a/uam_system_model/Pricing.py +++ b/uam_system_model/Pricing.py @@ -27,17 +27,27 @@ def optimize( uam_distance_matrix, optimality_gap, value_of_time, + utility_type, + beta_time=-0.0192, + beta_cost=-0.0353, uam_transition_time=10, time_limit=1800, CASM=0.79, verbose=True, ): if isinstance(value_of_time, int) or isinstance(value_of_time, float): - value_of_time = [value_of_time for _ in range(len(self.network.vertiport_dict))] - + value_of_time = [ + value_of_time for _ in range(len(self.network.vertiport_dict)) + ] + pax_arr = self.network.pax_arrival_times.copy() self.time_resolution = time_resolution self.num_vehicles = num_vehicles + self.beta_time = [ + beta_cost * value_of_time[i] / 60 + for i in range(len(self.network.vertiport_dict)) + ] + self.beta_cost = [beta_cost for _ in range(len(self.network.vertiport_dict))] pax_arr["passenger_arrival_time_slot"] = np.ceil( pax_arr["passenger_arrival_time_s"] / self.time_resolution / 60 @@ -93,9 +103,10 @@ def optimize( t_i_uam = [] first_last_mile_cost = [] flight_cost_uam = [] + beta_time_i = [] + beta_cost_i = [] p_i_bar = [] - for idx, row in self.pax_arr_grouped.iterrows(): origin = row["origin_vertiport_id"] destination = row["destination_vertiport_id"] @@ -126,9 +137,12 @@ def optimize( od_idx = max(origin, destination) uber_fare_i.append(uber_fare[od_idx]) elif len(uber_fare.shape) == 3: + od_idx = max(origin, destination) uber_fare_i.append(uber_fare[origin, destination, time]) p_i_bar.append(1 / value_of_time[od_idx]) + beta_time_i.append(self.beta_time[od_idx]) + beta_cost_i.append(self.beta_cost[od_idx]) distance = uam_distance_matrix[origin, destination] flight_cost_uam.append(CASM * 4 * distance) @@ -141,25 +155,30 @@ def optimize( non_zero_indices = [i for i, value in enumerate(self.di_bar) if value != 0] di_bar_selected_x = [self.di_bar[i] for i in non_zero_indices] - - if isinstance(value_of_time, float): - p_i_bar = [ - 1 / value_of_time for _ in non_zero_indices - ] # 32.63 is the VOT in dollars per minute - - else: - p_i_bar = [] - for idx, row in self.pax_arr_grouped.iterrows(): - origin = row["origin_vertiport_id"] - destination = row["destination_vertiport_id"] - od_idx = max(origin, destination) - p_i_bar.append(1 / (value_of_time[od_idx] / 60)) - v_i_bar_uber = -uber_travel_time_i - p_i_bar * uber_fare_i + if utility_type == "vot": + if isinstance(value_of_time, float): + p_i_bar = [ + 1 / value_of_time for _ in non_zero_indices + ] # 32.63 is the VOT in dollars per minute + + else: + p_i_bar = [] + for idx, row in self.pax_arr_grouped.iterrows(): + origin = row["origin_vertiport_id"] + destination = row["destination_vertiport_id"] + od_idx = max(origin, destination) + p_i_bar.append(1 / (value_of_time[od_idx] / 60)) + + v_i_bar_uber = -uber_travel_time_i - p_i_bar * uber_fare_i + elif utility_type == "betas": + beta_time_i = np.array(beta_time_i) + beta_cost_i = np.array(beta_cost_i) + v_i_bar_uber = beta_time_i * uber_travel_time_i + beta_cost_i * uber_fare_i bins = 20 max_flights = num_vehicles - eps = 0.01 + eps = 0.05 m = Model("Pricing Problem") m.Params.NonConvex = 2 @@ -215,7 +234,7 @@ def optimize( ) # Flow conservation constraints - for n in tqdm(self.nodes): + for n in self.nodes: m.addConstr( sum(m._x_vars[i, j] for i, j in self.edges if j == n) - sum(m._x_vars[i, j] for i, j in self.edges if i == n) @@ -243,36 +262,82 @@ def optimize( m._x_vars[self.edges[i]] * flight_cost_uam[i_non_zero] for i_non_zero, i in zip(range(len(di_bar_selected_x)), non_zero_indices) ) - cost_level_of_service = quicksum( - di_bar_selected_x[i_non_zero] - / p_i_bar[i_non_zero] - * m._x_inverse_vars[self.edges[i]] - * self.time_resolution - / 2 - for i_non_zero, i in zip(range(len(di_bar_selected_x)), non_zero_indices) - ) + if utility_type == "vot": + cost_level_of_service = quicksum( + di_bar_selected_x[i_non_zero] + / p_i_bar[i_non_zero] + * m._x_inverse_vars[self.edges[i]] + * self.time_resolution + / 2 + for i_non_zero, i in zip( + range(len(di_bar_selected_x)), non_zero_indices + ) + ) + theta_terms = quicksum( + di_bar_selected_x[i_non_zero] + / p_i_bar[i_non_zero] + * ( + m._theta_ln_theta_uam[i_non_zero] + - m._theta_ln_1_minus_theta_uam[i_non_zero] + ) + for i_non_zero in range(len(di_bar_selected_x)) + ) - theta_terms = quicksum( - di_bar_selected_x[i_non_zero] - / p_i_bar[i_non_zero] - * ( - m._theta_ln_theta_uam[i_non_zero] - - m._theta_ln_1_minus_theta_uam[i_non_zero] + other_terms = quicksum( + di_bar_selected_x[i_non_zero] + / p_i_bar[i_non_zero] + * m._theta_uam[i_non_zero] + * (v_i_bar_uber[i_non_zero] + t_i_uam[i_non_zero]) + for i_non_zero in range(len(di_bar_selected_x)) ) - for i_non_zero in range(len(di_bar_selected_x)) - ) - other_terms = quicksum( - di_bar_selected_x[i_non_zero] - / p_i_bar[i_non_zero] - * m._theta_uam[i_non_zero] - * (v_i_bar_uber[i_non_zero] + t_i_uam[i_non_zero]) - for i_non_zero in range(len(di_bar_selected_x)) - ) + objective = ( + theta_terms + other_terms + operating_cost + cost_level_of_service + ) + + elif utility_type == "betas": + theta_terms = quicksum( + -di_bar_selected_x[i_non_zero] + / beta_cost_i[i_non_zero] + * ( + m._theta_ln_theta_uam[i_non_zero] + - m._theta_ln_1_minus_theta_uam[i_non_zero] + ) + for i_non_zero in range(len(di_bar_selected_x)) + ) + + other_terms = quicksum( + -di_bar_selected_x[i_non_zero] + / beta_cost_i[i_non_zero] + * m._theta_uam[i_non_zero] + * ( + v_i_bar_uber[i_non_zero] + - beta_time_i[i_non_zero] * t_i_uam[i_non_zero] + ) + for i_non_zero in range(len(di_bar_selected_x)) + ) - objective = theta_terms + other_terms + operating_cost + cost_level_of_service + # objective = theta_terms + other_terms + operating_cost + cost_level_of_service + objective = theta_terms + other_terms + operating_cost m.setObjective(objective, GRB.MINIMIZE) + + if utility_type == "betas": + max_v = max_flights + for i, idx in enumerate(non_zero_indices): + current_los_costs = [] + factor = ( + di_bar_selected_x[i] * beta_time_i[i] * self.time_resolution + ) / (2 * beta_cost_i[i]) + for k in range(max_v + 1): + if k == 0: + current_los_costs.append(1e6) # Penalty for no service + else: + current_los_costs.append(factor * (1.0 / k)) + + m.setPWLObj( + m._x_vars[self.edges[idx]], range(max_v + 1), current_los_costs + ) if not verbose: m.setParam("OutputFlag", 0) m.update() @@ -308,17 +373,31 @@ def optimize( self.t_i_uam = t_i_uam self.value_of_time = value_of_time - output_merged["fare"] = output_merged.apply( - lambda row: self.calc_fare( - row, - 1/p_i_bar[int(row["flight_index"])], - v_i_bar_uber, - t_i_uam, - first_last_mile_cost, - self.time_resolution, - ), - axis=1, - ) + if utility_type == "vot": + output_merged["fare"] = output_merged.apply( + lambda row: self.calc_fare_vot( + row, + 1 / p_i_bar[int(row["flight_index"])], + v_i_bar_uber, + t_i_uam, + first_last_mile_cost, + self.time_resolution, + ), + axis=1, + ) + elif utility_type == "betas": + output_merged["fare"] = output_merged.apply( + lambda row: self.calc_fare_betas( + row, + beta_time_i, + beta_cost_i, + v_i_bar_uber, + t_i_uam, + first_last_mile_cost, + self.time_resolution, + ), + axis=1, + ) pax_arr_grouped_to_join = self.pax_arr_grouped.copy() pax_arr_grouped_to_join["flight_index"] = pax_arr_grouped_to_join.index @@ -355,7 +434,7 @@ def optimize( return df @staticmethod - def calc_fare( + def calc_fare_vot( row, value_of_time, v_i_bar_uber, t_i_uam, first_last_mile_cost, time_resolution ): num_flights = row["num_flights"] @@ -373,6 +452,41 @@ def calc_fare( return fare - first_last_cost + @staticmethod + def calc_fare_betas( + row, + beta_time, + beta_cost, + v_i_bar_uber, + t_i_uam, + first_last_mile_cost, + time_resolution, + ): + num_flights = row["num_flights"] + theta = row["percentage_uam"] + index = int(row["flight_index"]) + first_last_cost = first_last_mile_cost[index] + beta_time = beta_time[index] + beta_cost = beta_cost[index] + + if theta >= 1: + theta = 0.9999 + if theta <= 0: + theta = 0.0001 + + fare = ( + 1 + / beta_cost + * ( + math.log(theta) + - math.log(1 - theta) + + v_i_bar_uber[index] + - beta_time * (t_i_uam[index] + time_resolution / 2 / num_flights) + ) + ) + + return fare - first_last_cost + class FlightTask: def __init__( @@ -405,7 +519,10 @@ def next_task(self, next_task): + reposition_time ) - if ready_time <= next_task_start_time: + if ( + next_task_start_time - ready_time >= 0 + and next_task_start_time - ready_time <= 12 + ): return True else: return False diff --git a/uam_system_model/ScheduleGeneratorJFK.py b/uam_system_model/ScheduleGeneratorJFK.py index 37aed74..3cf2846 100644 --- a/uam_system_model/ScheduleGeneratorJFK.py +++ b/uam_system_model/ScheduleGeneratorJFK.py @@ -1,9 +1,7 @@ -import pandas as pd - import numpy as np +import pandas as pd -from .utils.schedule_utils import (build_schedules, - auto_regressive_poisson) +from .utils.schedule_utils import auto_regressive_poisson, build_schedules class ScheduleGenerator: @@ -20,46 +18,72 @@ def get_one_day( flight_distance_matrix, fare, ): - list_of_ods = arrival_rate_by_vertiport['od'].unique() + list_of_ods = arrival_rate_by_vertiport["od"].unique() - schedule = pd.DataFrame(columns=['schedule', 'od', 'num_pax', 'revenue']) - paxArrivalDf = pd.DataFrame(columns=['passenger_arrival_time_s', 'origin_vertiport_id', 'destination_vertiport_id']) + schedule = pd.DataFrame(columns=["schedule", "od", "num_pax", "revenue"]) + paxArrivalDf = pd.DataFrame( + columns=[ + "passenger_arrival_time_s", + "origin_vertiport_id", + "destination_vertiport_id", + ] + ) for od in list_of_ods: - od_data = arrival_rate_by_vertiport[arrival_rate_by_vertiport['od']==od] + od_data = arrival_rate_by_vertiport[arrival_rate_by_vertiport["od"] == od] realized_demand = auto_regressive_poisson( - rate=od_data['total_trips'].values / 365, - alpha=auto_regressive_alpha,) + rate=od_data["total_trips"].values / 365, + alpha=auto_regressive_alpha, + ) arrival_time = [] for hour in range(24): for _ in range(realized_demand[hour]): arrival_time.append(hour * 60 + np.random.randint(0, 60)) - departure_time, num_pax = build_schedules(arrival_time, max_waiting_time=max_waiting_time, occupancy=occupancy) + departure_time, num_pax = build_schedules( + arrival_time, max_waiting_time=max_waiting_time, occupancy=occupancy + ) - schedule = pd.concat([ - schedule, - pd.DataFrame({ - 'schedule': departure_time, - 'od': od, - 'num_pax': num_pax, - 'revenue': num_pax * flight_distance_matrix[vertiport_dict[od.split('_')[0]], vertiport_dict[od.split('_')[1]]] * fare # assume flat fare of $50 - }) - ], ignore_index=True) + schedule = pd.concat( + [ + schedule, + pd.DataFrame( + { + "schedule": departure_time, + "od": od, + "num_pax": num_pax, + "revenue": num_pax + * flight_distance_matrix[ + vertiport_dict[od.split("_")[0]], + vertiport_dict[od.split("_")[1]], + ] + * fare, # assume flat fare of $50 + } + ), + ], + ignore_index=True, + ) arrival_time = np.array(arrival_time) - paxArrivalDf = pd.concat([ - paxArrivalDf, - pd.DataFrame({ - 'passenger_arrival_time_s': arrival_time * 60, - 'origin_vertiport_id': vertiport_dict[od.split('_')[0]], - 'destination_vertiport_id': vertiport_dict[od.split('_')[1]], - }) - ], ignore_index=True) + paxArrivalDf = pd.concat( + [ + paxArrivalDf, + pd.DataFrame( + { + "passenger_arrival_time_s": arrival_time * 60, + "origin_vertiport_id": vertiport_dict[od.split("_")[0]], + "destination_vertiport_id": vertiport_dict[ + od.split("_")[1] + ], + } + ), + ], + ignore_index=True, + ) - schedule = schedule.sort_values(by=['schedule', "od"]).reset_index(drop=True) + schedule = schedule.sort_values(by=["schedule", "od"]).reset_index(drop=True) paxArrivalDf = paxArrivalDf.reset_index(drop=False) - paxArrivalDf = paxArrivalDf.rename(columns={'index':'passenger_id'}) + paxArrivalDf = paxArrivalDf.rename(columns={"index": "passenger_id"}) - return schedule, paxArrivalDf \ No newline at end of file + return schedule, paxArrivalDf diff --git a/uam_system_model/StarNetworkBJ.py b/uam_system_model/StarNetworkBJ.py index ce25441..b600265 100644 --- a/uam_system_model/StarNetworkBJ.py +++ b/uam_system_model/StarNetworkBJ.py @@ -158,7 +158,6 @@ def load_demand( return self.schedule, self.pax_arrival_times def plot_flight(self, ylim=(0, 25)): - # 假定：前两个是 hub，后面都是 spoke hubs = self.vertiports[:2] spokes = self.vertiports[2:] @@ -175,8 +174,7 @@ def plot_flight(self, ylim=(0, 25)): for hub_name in hubs: # 只保留与当前 hub 有关的航班（hub <-> 任意点） sub = schedule[ - (schedule["origin"] == hub_name) - | (schedule["destination"] == hub_name) + (schedule["origin"] == hub_name) | (schedule["destination"] == hub_name) ].copy() # 构造“局部视角”： @@ -208,7 +206,6 @@ def plot_flight(self, ylim=(0, 25)): input_to_viz[dummy_idx, :, :] = np.nan input_to_viz[:, dummy_idx, :] = np.nan - # 调用底层原始的 plot_travel_time（不改它） fig, ax = plot_travel_time( input_to_viz, @@ -218,10 +215,10 @@ def plot_flight(self, ylim=(0, 25)): ) # --- 覆盖 plot_travel_time 的标题 --- - hub_short = hub_name # 例如 "PEK", "PKX" + hub_short = hub_name # 例如 "PEK", "PKX" - ax[0].set_title(f"{hub_short}-Spokes") # 左图：hub → spokes - ax[1].set_title(f"Spokes-{hub_short}") # 右图：spokes → hub + ax[0].set_title(f"{hub_short}-Spokes") # 左图：hub → spokes + ax[1].set_title(f"Spokes-{hub_short}") # 右图：spokes → hub figs.append(fig) axes_list.append(ax) diff --git a/uam_system_model/StarNetworkJFK.py b/uam_system_model/StarNetworkJFK.py index 23abf7a..67c9990 100644 --- a/uam_system_model/StarNetworkJFK.py +++ b/uam_system_model/StarNetworkJFK.py @@ -27,13 +27,15 @@ def __init__( self.demand_generator = ScheduleGenerator() - def load_demand(self, - arrival_rate_by_vertiport, - auto_regressive_alpha=0, - max_waiting_time=5, - occupancy=4, - fare=3, - seed=9): + def load_demand( + self, + arrival_rate_by_vertiport, + auto_regressive_alpha=0, + max_waiting_time=5, + occupancy=4, + fare=3, + seed=9, + ): np.random.seed(seed) ( schedule, diff --git a/uam_system_model/utils/PricingOpSummary.py b/uam_system_model/utils/PricingOpSummary.py index de0ec87..8d1d0ec 100644 --- a/uam_system_model/utils/PricingOpSummary.py +++ b/uam_system_model/utils/PricingOpSummary.py @@ -56,11 +56,16 @@ def get_summary_statistics(self): return summary_stats - def plot_average_rasm(self, dpi=300): + def plot_average_rasm(self, ax=None, dpi=300): + if ax is None: + fig, ax = plt.subplots(figsize=(6, 4), dpi=dpi) + else: + fig = ax.figure + df_grouped = self.policy.groupby(["passenger_arrival_time_slot"])[ "rev_per_mile" ].mean() - fig, ax = plt.subplots(figsize=(6, 4), dpi=dpi) + sns.lineplot( x=np.arange(len(df_grouped)), y=df_grouped.values, @@ -87,6 +92,7 @@ def plot_rasm_by_od(self, dpi=300): data=self.policy, x="passenger_arrival_time_slot", hue="markets", + hue_order=self.network.vertiports[1:], y="rev_per_mile", marker="o", err_style=None, @@ -107,6 +113,34 @@ def plot_rasm_by_od(self, dpi=300): return fig, ax + def plot_fare_by_od(self, dpi=300): + fig, ax = plt.subplots(figsize=(6, 4), dpi=dpi) + sns.lineplot( + data=self.policy, + x="passenger_arrival_time_slot", + y="fare", + hue="markets", + hue_order=self.network.vertiports[1:], + err_style=None, + legend=True, + marker="o", + palette=custom_colors, + ax=ax, + ) + ax.set( + xticks=np.arange(0, 49, 12), + xticklabels=[str(i) + ":00" for i in range(0, 26, 6)], + xlabel="", + ylabel="Fare ($)", + xlim=(0, 48), + title="Average Fare by ODs", + ) + ax.xaxis.set_major_locator(MultipleLocator(2)) + ax.grid(True, linestyle="--", alpha=0.5) + ax.legend(bbox_to_anchor=(1.02, 0.85), loc="upper left", borderaxespad=0.0) + + return fig, ax + def plot_revenue_by_od(self, dpi=300): fig, ax = plt.subplots(figsize=(6, 4), dpi=dpi) sns.lineplot( @@ -114,6 +148,7 @@ def plot_revenue_by_od(self, dpi=300): x="passenger_arrival_time_slot", y="total_revenue", hue="markets", + hue_order=self.network.vertiports[1:], err_style=None, legend=True, marker="o", @@ -134,13 +169,17 @@ def plot_revenue_by_od(self, dpi=300): return fig, ax - def plot_uam_share_by_od(self, dpi=300): - fig, ax = plt.subplots(figsize=(6, 4), dpi=dpi) + def plot_uam_share_by_od(self, ax=None, dpi=300): + if ax is None: + fig, ax = plt.subplots(figsize=(6, 4), dpi=dpi) + else: + fig = ax.figure sns.lineplot( data=self.policy, x="passenger_arrival_time_slot", y="percentage_uam", hue="markets", + hue_order=self.network.vertiports[1:], err_style=None, legend=True, marker="o", diff --git a/uam_system_model/utils/visualize.py b/uam_system_model/utils/visualize.py index ed4425a..73da623 100644 --- a/uam_system_model/utils/visualize.py +++ b/uam_system_model/utils/visualize.py @@ -25,7 +25,7 @@ def plot_travel_time( travel_time_avg, vertiports, ylim=(0, 100), ylabel="TNC Trip Time (min)" ): fig, ax = plt.subplots(figsize=(12, 4), ncols=2, dpi=300) - for i in range(1, 9): + for i in range(1, len(vertiports)): ax[0].plot( np.arange(24), travel_time_avg[0, i, :],