import math import sys import weakref from collections import defaultdict from dataclasses import dataclass from functools import partial from threading import Timer from typing import Callable import pandas as pd from ortools.constraint_solver import pywrapcp, routing_enums_pb2 from pandas import DataFrame from core.types.Logger import Logger VEHICLE_COST = 16 * 3600 # Two working days. VEHICLE_PRIORITY_COST = 0 # Vehicle with priority has zero cost. VEHICLE_DUPLICATE_COST = 100_000_000 VEHICLE_DUPLICATE_FACTOR = 2 """ id name route_type ... cost max_time range 0 0 Kolo z pomožnim motorjem kpm ... 57600 3600 60000 1 1 Motorno kolo mk ... 0 3600 120000 2 2 Kolo z motorjem km ... 57600 3600 120000 3 3 Kolo bike ... 57600 3600 30000 4 4 Elektricni tro/štiri kolesnik ev ... 57600 3600 120000 5 5 Pes foot ... 57600 3600 6000 6 6 Avtomobil car ... 57600 3600 150000 7 0 Kolo z pomožnim motorjem kpm ... 100000000 3600 60000 8 1 Motorno kolo mk ... 100000000 3600 120000 9 2 Kolo z motorjem km ... 100000000 3600 120000 10 3 Kolo bike ... 100000000 3600 30000 11 4 Elektricni tro/štiri kolesnik ev ... 100000000 3600 120000 12 5 Pes foot ... 100000000 3600 6000 13 6 Avtomobil car ... 100000000 3600 150000 """ @dataclass class VrpInstance: """ Main "Instance" of the data to optimize """ vehicles: pd.DataFrame nodes: pd.DataFrame dist: dict time: dict initial_routes: list[list[int]] district_percentage: float def read_solution( manager: pywrapcp.RoutingIndexManager, routing: pywrapcp.RoutingModel, instance: VrpInstance, distance_evaluators: dict[callable], time_evaluators: dict[callable], ): routes = [] for vehicle_id, route_type in enumerate(instance.vehicles["route_type"]): distance_evaluator = distance_evaluators[route_type] time_evaluator = time_evaluators[route_type] points = [] route_distance = 0 route_time = 0 route_cost = 0 index = routing.Start(vehicle_id) while not routing.IsEnd(index): previous_index = index index = routing.NextVar(index).Value() route_distance += distance_evaluator(previous_index, index) route_time += time_evaluator(previous_index, index) route_cost += routing.GetArcCostForVehicle(previous_index, index, vehicle_id) node = manager.IndexToNode(index) point = instance.nodes.base_point.iloc[node] points.append(point) routes.append( { "vehicle": vehicle_id, "type": instance.vehicles.iloc[vehicle_id]["route_type"], "route": points, "total_distance": route_distance, "total_time": route_time, "total_cost": route_cost, "num_points": len(points), } ) routes = pd.DataFrame(routes) return routes class RepeatTimer(Timer): def run(self): while not self.finished.wait(self.interval): self.function() class SolutionCallback: def __init__( self, manager: pywrapcp.RoutingIndexManager, model: pywrapcp.RoutingModel, instance: VrpInstance, distance_evaluators: dict[callable], time_evaluators: dict[callable], solution_callback_fn: Callable[[int, pd.DataFrame], None], stop_callback_fn: callable ): self._routing_manager_ref = weakref.ref(manager) self._routing_model_ref = weakref.ref(model) self.objectives = [] self.instance = instance self.distance_evaluators = distance_evaluators self.time_evaluators = time_evaluators self.best_routes = None self.solution_callback_fn = solution_callback_fn self.stop_callback_fn = stop_callback_fn self._timer = RepeatTimer(10, self._check_terminated) self._timer.start() def __call__(self): # current objective value objective = int(self._routing_model_ref().CostVar().Value()) if not self.objectives or objective < self.objectives[-1]: self.objectives.append(objective) self.best_routes = read_solution( self._routing_manager_ref(), self._routing_model_ref(), self.instance, self.distance_evaluators, self.time_evaluators ) tmp = self.best_routes tmp = tmp[tmp["num_points"] > 2] vpd = defaultdict(set) districts = self.instance.nodes['district'].values for _, row in tmp.iterrows(): for p in row['route']: vpd[districts[p]].add(row['vehicle']) self.solution_callback_fn(objective, self.best_routes) # Num. clean districts: {sum(len(s) == 1 for s in vpd.values())} / {len(vpd.keys())} ") # log.info(f"Objective: {objective} Num. vehicles: {len(tmp)}") # self._routing_model_ref().solver().FinishCurrentSearch() def _check_terminated(self): """ if self.stop_callback_fn(None): ^^^^^^^^^^^^^^^^^^^^^^^^^^^ TypeError: SolvesallOptimizationService.vrpOptimization..stop_callback_fn() takes 0 positional arguments but 1 was given """ if self.stop_callback_fn(): self._timer.cancel() self._routing_model_ref().solver().FinishCurrentSearch() def solve(instance: VrpInstance, config, time_limit_sec, solution_callback_fn: Callable[[int, pd.DataFrame], None], stop_callback_fn, log: Logger, log_search=False): # with open(f"solve_args_{datetime.now().isoformat()}.pkl", "wb") as f: # pickle.dump((instance, config), f) sys.stdout.flush() assert config.objective in ['distance', 'time'] assert instance.nodes.iloc[0]["type"] == "depot", "Depot is expected to be at 0" manager = pywrapcp.RoutingIndexManager( len(instance.nodes), len(instance.vehicles), 0 ) routing = pywrapcp.RoutingModel(manager) def create_distance_evaluator(route_type, instance): dist_mat = instance.dist[route_type] base_point = instance.nodes["base_point"].values freq = instance.nodes['freq'].values def distance_evaluator(from_node, to_node): dst_node = manager.IndexToNode(to_node) src = base_point[manager.IndexToNode(from_node)] dst = base_point[manager.IndexToNode(to_node)] return round(dist_mat[src, dst]) return distance_evaluator distance_evaluators, distance_evaluators_index = {}, {} for route_type in instance.vehicles["route_type"].unique(): distance_evaluators[route_type] = create_distance_evaluator(route_type, instance) distance_evaluators_index[route_type] = routing.RegisterTransitCallback( distance_evaluators[route_type] ) def create_time_evaluator(route_type, instance): dist_mat = instance.dist[route_type] time_mat = instance.time[route_type] base_point = instance.nodes["base_point"].values service_time = instance.nodes["service_time"].values freq = instance.nodes['freq'].values hisa_ids = instance.nodes['hisa_id'].values def time_evaluator(from_node, to_node): src_node = manager.IndexToNode(from_node) dst_node = manager.IndexToNode(to_node) src = base_point[manager.IndexToNode(from_node)] dst = base_point[manager.IndexToNode(to_node)] src_hisa_id = hisa_ids[src] dst_hisa_id = hisa_ids[dst] # THIS MUST BE IN SYNC WITH Run_optimization_job.save WHERE OPTIMIZATION ROUTE IS CALCULATED!!! time = round(time_mat[src, dst] + freq[src_node] * service_time[src_node]) # log.info(f"({src} -> {dst} [{src_hisa_id} -> {dst_hisa_id}] [distance={dist_mat[src, dst]} time={time_mat[src, dst]} freq={freq[src_node]} service_time={service_time[src_node]}] = {time}") return time return time_evaluator time_evaluators, time_evaluators_index = {}, {} for route_type in instance.vehicles["route_type"].unique(): time_evaluators[route_type] = create_time_evaluator(route_type, instance) time_evaluators_index[route_type] = routing.RegisterTransitCallback( time_evaluators[route_type] ) def create_demand_evaluator(instance): demands = instance.nodes["demand"].values def demand_evaluator(from_node): return int(demands[manager.IndexToNode(from_node)]) return demand_evaluator demand_evaluator = create_demand_evaluator(instance) demand_evaluator_index = routing.RegisterUnaryTransitCallback(demand_evaluator) routing.AddDimensionWithVehicleTransitAndCapacity( [ distance_evaluators_index[route_type] for route_type in instance.vehicles["route_type"] ], 0, [1000000] * len(instance.vehicles), # [int(x) for x in instance.vehicles["range"]] if not config.set_initial else [1000000] * len(instance.vehicles), True, "Distance", ) """ With initial solution we must be aware that is in the feasable space. If it is not in the feasable space the solver can fail because it does not find an initial solution. That's why we will increase the vehicle time constraint to 10 hours, and create a soft penalty. On initial routes max_time constraint on vehicle is overacheived. """ routing.AddDimensionWithVehicleTransitAndCapacity( [ time_evaluators_index[route_type] for route_type in instance.vehicles["route_type"] ], 0, [int(x) for x in instance.vehicles["max_time"]] if not config.set_initial else [1000 * 3600] * len(instance.vehicles), True, "Time", ) routing.AddConstantDimension(1, len(instance.nodes), True, "Count") count_dimension = routing.GetDimensionOrDie("Count") for vehicle_id in range(len(instance.vehicles)): if instance.vehicles.iloc[vehicle_id]['cost'] == 0: index_end = routing.End(vehicle_id) count_dimension.SetCumulVarSoftLowerBound(index_end, 3, 1_000_000_000) routing.SetVehicleUsedWhenEmpty(True, vehicle_id) if config.set_initial: time_dimension = routing.GetDimensionOrDie('Time') for vehicle_id in range(len(instance.vehicles)): index = routing.End(vehicle_id) max_time = int(instance.vehicles.iloc[vehicle_id]['max_time']) time_dimension.SetCumulVarSoftUpperBound(index, max_time, 1_000) routing.AddDimensionWithVehicleCapacity( demand_evaluator_index, 0, [1000000] * len(instance.vehicles), # [int(x) for x in instance.vehicles["capacity"]], True, "Capacity", ) # District matching if config.set_initial: log.info("District matching ..") node_to_vehicle = {} district_size = {} for v, route in enumerate(instance.initial_routes): for n in route: node_to_vehicle[n] = v district_size[v] = len(route) def district_added_callback(vehicle_id, from_index): from_node = manager.IndexToNode(from_index) if from_node == 0: # If node == 0, then it is depo. return 1 # Check if node not belongs to vehicle's initial district return 1 if vehicle_id != node_to_vehicle[from_node] else 0 def district_required_callback(vehicle_id, from_index): from_node = manager.IndexToNode(from_index) if from_node == 0: # If node == 0, then it is depo. return 1 # Check if node belongs to vehicle's initial district return 1 if vehicle_id == node_to_vehicle[from_node] else 0 routing.AddDimensionWithVehicleTransitAndCapacity( [routing.RegisterUnaryTransitCallback(partial(district_added_callback, vehicle_id)) for vehicle_id in range(len(instance.vehicles)) ], 0, [len(instance.nodes)] * len(instance.vehicles), True, "District_added", ) routing.AddDimensionWithVehicleTransitAndCapacity( [routing.RegisterUnaryTransitCallback(partial(district_required_callback, vehicle_id)) for vehicle_id in range(len(instance.vehicles)) ], 0, [len(instance.nodes)] * len(instance.vehicles), True, "District_required", ) district_added_dimension = routing.GetDimensionOrDie('District_added') district_required_dimension = routing.GetDimensionOrDie('District_required') # Add soft lower bound for each vehicle for vehicle_id in range(len(instance.vehicles)): if vehicle_id not in district_size: continue # len(IR) * (1 - 0.8 (GASPER)) added_visits = int(district_size[vehicle_id] * (1 - instance.district_percentage)) # 80 % of district size index = routing.End(vehicle_id) district_added_dimension.SetCumulVarSoftUpperBound(index, added_visits, 10_000) district_required_dimension.SetCumulVarSoftLowerBound(index, 3, 10_000) # District must contains 3 initial points # One vehicle per street (or district) # if config.district_mode == 'single' and config.district_penalty > 0: # for _, ids in instance.nodes.groupby('district')['id']: # ids = [manager.NodeToIndex(x) for x in ids.values] # assert 0 not in ids, "Depot can't have an assigned district." # routing.AddSoftSameVehicleConstraint(ids, config.district_penalty) # elif config.district_mode == 'subsets' and config.district_penalty > 0: # for _, ids in instance.nodes.groupby('district')['id']: # ids = [manager.NodeToIndex(x) for x in ids.values] # assert 0 not in ids, "Depot can't have an assigned district." # log.info("Building pairwise constraints ...", end="") ## sys.stdout.flush() # combs = list(itertools.combinations(ids, 2))[:40] # combs.append(ids) # for subset in combs: # routing.AddSoftSameVehicleConstraint(subset, config.district_penalty) # log.info("finished") # elif config.district_mode == 'hard': # solver = routing.solver() # for _, ids in instance.nodes.groupby('district')['id']: # ids = [manager.NodeToIndex(x) for x in ids.values] # # v0 = routing.VehicleVar(ids[0]) # for i in ids[1:]: # solver.Add(v0 == routing.VehicleVar(i)) def create_objective_evaluator(route_type, instance): dist_mat = instance.dist[route_type] time_mat = instance.time[route_type] base_point = instance.nodes["base_point"].values service_time = instance.nodes["service_time"].values freq = instance.nodes['freq'].values hisa_ids = instance.nodes['hisa_id'].values def objective_evaluator(from_node, to_node): src_node = manager.IndexToNode(from_node) dst_node = manager.IndexToNode(to_node) src = base_point[manager.IndexToNode(from_node)] dst = base_point[manager.IndexToNode(to_node)] src_hisa_id = hisa_ids[src] dst_hisa_id = hisa_ids[dst] # THIS MUST BE IN SYNC WITH Run_optimization_job.save WHERE OPTIMIZATION ROUTE IS CALCULATED!!! if dist_mat[src, dst] > 3000: penalty = dist_mat[src, dst] else: distance = dist_mat[src, dst] max_distance_sqrt = math.sqrt(3000) penalty = (distance / max_distance_sqrt) ** 2 if config.useDistrictCentrality: total_cost = round(time_mat[src, dst] + freq[src_node] * service_time[src_node] + penalty) else: total_cost = round(time_mat[src, dst] + freq[src_node] * service_time[src_node]) # log.info(f"({src} -> {dst} [{src_hisa_id} -> {dst_hisa_id}] [distance={dist_mat[src, dst]} time={time_mat[src, dst]} freq={freq[src_node]} service_time={service_time[src_node]}] = {time}") return total_cost return objective_evaluator objective_evaluators, objective_evaluators_index = {}, {} for route_type in instance.vehicles["route_type"].unique(): objective_evaluators[route_type] = create_objective_evaluator(route_type, instance) objective_evaluators_index[route_type] = routing.RegisterTransitCallback( objective_evaluators[route_type] ) # Objective if config.objective == 'distance': obj_evaluators_index = distance_evaluators_index obj_dimension = routing.GetDimensionOrDie('Distance') elif config.objective == 'time': obj_evaluators_index = time_evaluators_index obj_dimension = routing.GetDimensionOrDie('Time') obj_evaluators_index = objective_evaluators_index # sum of distances (or travel times) for i, route_type in enumerate(instance.vehicles["route_type"]): routing.SetArcCostEvaluatorOfVehicle(obj_evaluators_index[route_type], i) # diff between max and min distance (or travel time) # obj_dimension.SetGlobalSpanCostCoefficient(100) # cost per each used vehicle for i, cost in enumerate(instance.vehicles["cost"]): routing.SetFixedCostOfVehicle(int(cost), i) solution_callback = SolutionCallback(manager, routing, instance, distance_evaluators, time_evaluators, solution_callback_fn, stop_callback_fn) routing.AddAtSolutionCallback(solution_callback) search_parameters = pywrapcp.DefaultRoutingSearchParameters() search_parameters.first_solution_strategy = ( routing_enums_pb2.FirstSolutionStrategy.LOCAL_CHEAPEST_COST_INSERTION ) search_parameters.local_search_metaheuristic = ( routing_enums_pb2.LocalSearchMetaheuristic.GUIDED_LOCAL_SEARCH ) search_parameters.time_limit.FromSeconds(time_limit_sec) search_parameters.log_search = log_search if config.set_initial: log.info("Initial solution added.") routing.CloseModelWithParameters(search_parameters) initial_solution = routing.ReadAssignmentFromRoutes(instance.initial_routes, True) assert initial_solution is not None, "Initial solution is not feasible." log.info("Initial solution found!") solution = routing.SolveFromAssignmentWithParameters( initial_solution, search_parameters ) else: solution = routing.SolveWithParameters(search_parameters) # Stop callback timer sice we dont need it anymore solution_callback._timer.cancel() assert solution, "No solution found." if log_search: debug_solution(instance.vehicles, instance.nodes, manager, routing, solution, log) obj, sol = solution.ObjectiveValue(), solution_callback.best_routes if config.set_initial: debug_solution_overrlapping(instance.initial_routes, sol, log) return obj, sol def debug_solution(vehicles, points, manager, routing, solution, log: Logger): objectiveValue: float = solution.ObjectiveValue() distanceDimension = routing.GetMutableDimension("Distance") timeDimension = routing.GetMutableDimension("Time") log.info(f"Objective value: {objectiveValue}") total_time = 0 total_distance = 0 total_cost = 0 for vehicle_idx in range(len(vehicles)): # add first node index = routing.Start(vehicle_idx) node = manager.IndexToNode(index) point = points.iloc[node].to_dict() log.info(f"Route for vehicle {vehicle_idx} = {vehicles.iloc[vehicle_idx].to_dict()}:") route_time = 0 route_distance = 0 route_cost = 0 start = True while not routing.IsEnd(index): # log.info(f"\t{node} = {point}") # Previous info ctime = solution.Value(timeDimension.CumulVar(index)) cdistance = solution.Value(distanceDimension.CumulVar(index)) # Next index previous_index = index index = solution.Value(routing.NextVar(index)) # Next info ntime = solution.Value(timeDimension.CumulVar(index)) ndistance = solution.Value(distanceDimension.CumulVar(index)) time = ntime - ctime distance = ndistance - cdistance cost = routing.GetArcCostForVehicle(previous_index, index, vehicle_idx) if start: log.info(f"STARTING COST: {cost}") start = False # log.info(f"\tCurrent time: {round(time / 3600, 3)}h") # log.info(f"\tCurrent distance: {round(distance, 3)}m") # log.info(f"\tCurrent cost: {round(cost / 3600, 3)}\n") route_time += time route_distance += distance route_cost += cost node = manager.IndexToNode(index) point = points.iloc[node].to_dict() # log.info(f"\t{node} = {point}") log.info(f"Route time: {round(route_time / 3600, 3)}h") log.info(f"Route distance: {round(route_distance, 3)}m") log.info(f"Route cost: {round(route_cost, 3)}\n") total_time += route_time total_distance += route_distance total_cost += route_cost log.info(f"\nAll routes time: {round(total_time / 3600, 3)}h") log.info(f"All routes distance: {round(total_distance, 3)}m") log.info(f"All routes cost: {round(total_cost, 3)}") def debug_solution_overrlapping(initial_routes: list[list[int]], solution: DataFrame, log: Logger): for id, vehicle, type, route, total_distance, total_time, total_cost, num_points in solution.to_records(): if len(initial_routes) == id: break initial_route = set(initial_routes[id]) route = set(route) crosSection = initial_route.intersection(route) if len(initial_route) > 0: log.info(f"Vehicle {id}. overlappings: {round(100 * len(crosSection) / len(initial_route), 1)}%")