|
- import pandas as pd
- import numpy as np
- from matplotlib import pyplot as plt
- from collections import namedtuple
- from ortools.constraint_solver import pywrapcp
- from ortools.constraint_solver import routing_enums_pb2
- from datetime import datetime, timedelta
-
- import datetime
-
-
- class Customers():
- # num_stops, demand, log,lat, min_tw, max_tw
- def __init__(self,
- num_stops,
- demands,
- lats,
- lons,
- tw_opens,
- tw_closes):
- self.number = num_stops #: The number of customers and depots
- # The 'name' of the stop, indexed from 0 to num_stops-1
- stops = np.array(range(0, num_stops))
-
- self.time_horizon = 24 * 60**2 # A 24 hour period.
-
- # The customers demand min_tw to max_tw hour time window for each
- # Make random timedeltas, nominally from the start of the day.
- start_times = []
- stop_times = []
- for idx in range(self.number):
- start_times.append(timedelta(seconds=tw_opens[idx]))
- stop_times.append(timedelta(seconds=tw_closes[idx]))
- # A named tuple for the customer
- Customer = namedtuple(
- 'Customer',
- [
- 'index', # the index of the stop
- 'demand', # the demand for the stop
- 'lat', # the latitude of the stop
- 'lon', # the longitude of the stop
- 'tw_open', # timedelta window open
- 'tw_close'
- ]) # timedelta window cls
-
- self.customers = [
- Customer(idx, dem, lat, lon, tw_open, tw_close)
- for idx, dem, lat, lon, tw_open, tw_close in zip(
- stops, demands, lats, lons, start_times, stop_times)
- ]
-
- # unload speed: 30 cubic meters per hour.
- # The number of seconds needed to 'unload' 1 cubic meters (2 min).
- self.service_time_per_dem = 120 / 1000 # seconds
-
- def set_manager(self, manager):
- self.manager = manager
-
- def central_start_node(self, invert=False):
- """
- Return a random starting node, with probability weighted by distance
- from the centre of the extents, so that a central starting node is
- likely.
-
- Args: invert (Optional bool): When True, a peripheral starting node is
- most likely.
-
- Returns:
- int: a node index.
-
- Examples:
- >>> customers.central_start_node(invert=True)
- 42
- """
- return 0
-
- def make_distance_mat(self, method='haversine'):
- """
- Return a distance matrix and make it a member of Customer, using the
- method given in the call. Currently only Haversine (GC distance) is
- implemented, but Manhattan, or using a maps API could be added here.
- Raises an AssertionError for all other methods.
-
- Args: method (Optional[str]): method of distance calculation to use. The
- Haversine formula is the only method implemented.
-
- Returns:
- Numpy array of node to node distances.
-
- Examples:
- >>> dist_mat = customers.make_distance_mat(method='haversine')
- >>> dist_mat = customers.make_distance_mat(method='manhattan')
- AssertionError
- """
- self.distmat = np.zeros((self.number, self.number))
- methods = {'haversine': self._haversine}
- assert (method in methods)
- for frm_idx in range(self.number):
- for to_idx in range(self.number):
- if frm_idx != to_idx:
- frm_c = self.customers[frm_idx]
- to_c = self.customers[to_idx]
- self.distmat[frm_idx, to_idx] = self._haversine(
- frm_c.lon, frm_c.lat, to_c.lon, to_c.lat)
- return (self.distmat)
-
- def _haversine(self, lon1, lat1, lon2, lat2):
- """
- Calculate the great circle distance between two points
- on the earth specified in decimal degrees of latitude and longitude.
- https://en.wikipedia.org/wiki/Haversine_formula
-
- Args:
- lon1: longitude of pt 1,
- lat1: latitude of pt 1,
- lon2: longitude of pt 2,
- lat2: latitude of pt 2
-
- Returns:
- the distace in km between pt1 and pt2
- """
- # convert decimal degrees to radians
- lon1, lat1, lon2, lat2 = map(np.radians, [lon1, lat1, lon2, lat2])
-
- # haversine formula
- dlon = lon2 - lon1
- dlat = lat2 - lat1
- a = (np.sin(dlat / 2)**2 +
- np.cos(lat1) * np.cos(lat2) * np.sin(dlon / 2)**2)
- c = 2 * np.arcsin(np.sqrt(a))
-
- # 6367 km is the radius of the Earth
- km = 6367 * c
- return km
-
- def get_total_demand(self):
- """
- Return the total demand of all customers.
- """
- return (sum([c.demand for c in self.customers]))
-
- def return_dist_callback(self, **kwargs):
- """
- Return a callback function for the distance matrix.
-
- Args: **kwargs: Arbitrary keyword arguments passed on to
- make_distance_mat()
-
- Returns:
- function: dist_return(a,b) A function that takes the 'from' node
- index and the 'to' node index and returns the distance in km.
- """
- self.make_distance_mat(**kwargs)
-
- def dist_return(from_index, to_index):
- # Convert from routing variable Index to distance matrix NodeIndex.
- from_node = self.manager.IndexToNode(from_index)
- to_node = self.manager.IndexToNode(to_index)
- return (self.distmat[from_node][to_node])
-
- return dist_return
-
- def return_dem_callback(self):
- """
- Return a callback function that gives the demands.
-
- Returns:
- function: dem_return(a) A function that takes the 'from' node
- index and returns the distance in km.
- """
-
- def dem_return(from_index):
- # Convert from routing variable Index to distance matrix NodeIndex.
- from_node = self.manager.IndexToNode(from_index)
- return (self.customers[from_node].demand)
-
- return dem_return
-
- def zero_depot_demands(self, depot):
- """
- Zero out the demands and time windows of depot. The Depots do not have
- demands or time windows so this function clears them.
-
- Args: depot (int): index of the stop to modify into a depot.
- Examples: >>> customers.zero_depot_demands(5) >>>
- customers.customers[5].demand == 0 True
- """
- start_depot = self.customers[depot]
- self.customers[depot] = start_depot._replace(
- demand=0, tw_open=None, tw_close=None)
-
- def make_service_time_call_callback(self):
- """
- Return a callback function that provides the time spent servicing the
- customer. Here is it proportional to the demand given by
- self.service_time_per_dem, default 300 seconds per unit demand.
-
- Returns:
- function [dem_return(a, b)]: A function that takes the from/a node
- index and the to/b node index and returns the service time at a
-
- """
-
- def service_time_return(a, b):
- return (self.customers[a].demand * self.service_time_per_dem)
-
- return service_time_return
-
- def make_transit_time_callback(self, speed_kmph=30):
- """
- Creates a callback function for transit time. Assuming an average
- speed of speed_kmph
- Args:
- speed_kmph: the average speed in km/h
-
- Returns:
- function [transit_time_return(a, b)]: A function that takes the
- from/a node index and the to/b node index and returns the
- transit time from a to b.
- """
-
- def transit_time_return(a, b):
- return (self.distmat[a][b] / (speed_kmph * 1.0 / 60**2))
-
- return transit_time_return
-
-
- class Vehicles():
- """
- A Class to create and hold vehicle information.
-
- The Vehicles in a CVRPTW problem service the customers and belong to a
- depot. The class Vehicles creates a list of named tuples describing the
- Vehicles. The main characteristics are the vehicle capacity, fixed cost,
- and cost per km. The fixed cost of using a certain type of vehicles can be
- higher or lower than others. If a vehicle is used, i.e. this vehicle serves
- at least one node, then this cost is added to the objective function.
-
- Note:
- If numpy arrays are given for capacity and cost, then they must be of
- the same length, and the number of vehicles are inferred from them.
- If scalars are given, the fleet is homogeneous, and the number of
- vehicles is determined by number.
-
- Args: capacity (scalar or numpy array): The integer capacity of demand
- units. cost (scalar or numpy array): The fixed cost of the vehicle. number
- (Optional [int]): The number of vehicles in a homogeneous fleet.
- """
-
- def __init__(self, capacity=12000, cost=1, number=None):
-
- Vehicle = namedtuple('Vehicle', ['index', 'capacity', 'cost'])
-
- if number is None:
- self.number = np.size(capacity)
- else:
- self.number = number
- idxs = np.array(range(0, self.number))
-
- if np.isscalar(capacity):
- capacities = capacity * np.ones_like(idxs)
- elif np.size(capacity) != self.number:
- print('capacity is neither scalar, nor the same size as num!')
- else:
- capacities = capacity
-
- if np.isscalar(cost):
- costs = cost * np.ones_like(idxs)
- elif np.size(cost) != self.number:
- print(np.size(cost))
- print('cost is neither scalar, nor the same size as num!')
- else:
- costs = cost
-
- self.vehicles = [
- Vehicle(idx, capacity, cost)
- for idx, capacity, cost in zip(idxs, capacities, costs)
- ]
-
- def get_total_capacity(self):
- return (sum([c.capacity for c in self.vehicles]))
-
- def return_starting_callback(self, customers, sameStartFinish=False):
- # create a different starting and finishing depot for each vehicle
- self.starts = [0 for o in range(self.number)]
- self.ends = self.starts
- # the depots will not have demands, so zero them.
- for depot in self.starts:
- customers.zero_depot_demands(depot)
- for depot in self.ends:
- customers.zero_depot_demands(depot)
-
- def start_return(v):
- return (self.starts[v])
-
- return start_return
-
-
- def discrete_cmap(N, base_cmap=None):
- """
- Create an N-bin discrete colormap from the specified input map
- """
- # Note that if base_cmap is a string or None, you can simply do
- # return plt.cm.get_cmap(base_cmap, N)
- # The following works for string, None, or a colormap instance:
-
- base = plt.cm.get_cmap(base_cmap)
- color_list = base(np.linspace(0, 1, N))
- cmap_name = base.name + str(N)
- return base.from_list(cmap_name, color_list, N)
-
-
- def vehicle_output_string(manager, routing, plan):
- """
- Return a string displaying the output of the routing instance and
- assignment (plan).
-
- Args: routing (ortools.constraint_solver.pywrapcp.RoutingModel): routing.
- plan (ortools.constraint_solver.pywrapcp.Assignment): the assignment.
-
- Returns:
- (string) plan_output: describing each vehicle's plan.
-
- (List) dropped: list of dropped orders.
-
- """
- dropped = []
- for order in range(routing.Size()):
- if (plan.Value(routing.NextVar(order)) == order):
- dropped.append(str(order))
-
- capacity_dimension = routing.GetDimensionOrDie('Capacity')
- time_dimension = routing.GetDimensionOrDie('Time')
- plan_output = ''
-
- for route_number in range(routing.vehicles()):
- order = routing.Start(route_number)
- plan_output += 'Route {0}:'.format(route_number)
- if routing.IsEnd(plan.Value(routing.NextVar(order))):
- plan_output += ' Empty \n'
- else:
- while True:
- load_var = capacity_dimension.CumulVar(order)
- time_var = time_dimension.CumulVar(order)
- node = manager.IndexToNode(order)
- plan_output += \
- ' {node} Load({load}) Time({tmin}, {tmax}) -> '.format(
- node=node,
- load=plan.Value(load_var),
- tmin=str(timedelta(seconds=plan.Min(time_var))),
- tmax=str(timedelta(seconds=plan.Max(time_var))))
-
- if routing.IsEnd(order):
- print(tlast) # new
- plan_output += ' EndRoute {0}. \n'.format(route_number)
- break
- else:
- tlast = str(timedelta(seconds=plan.Min(time_var))) # new
- order = plan.Value(routing.NextVar(order))
- plan_output += '\n'
-
- return (plan_output, dropped)
-
-
- def build_vehicle_route(manager, routing, plan, customers, veh_number, num_stops):
- """
- Build a route for a vehicle by starting at the strat node and
- continuing to the end node.
-
- Args: routing (ortools.constraint_solver.pywrapcp.RoutingModel): routing.
- plan (ortools.constraint_solver.pywrapcp.Assignment): the assignment.
- customers (Customers): the customers instance. veh_number (int): index of
- the vehicle
-
- Returns:
- (List) route: indexes of the customers for vehicle veh_number
- """
- veh_used = routing.IsVehicleUsed(plan, veh_number)
- print('Vehicle {0} is used {1}'.format(veh_number, veh_used))
- if veh_used:
- route = []
- vload = 0 # new
- node = routing.Start(veh_number) # Get the starting node index
- route.append(customers.customers[manager.IndexToNode(node)])
- while not routing.IsEnd(node):
- route.append(customers.customers[manager.IndexToNode(node)])
- if node < num_stops: # new
- vload += customers.customers[node].demand # new
- node = plan.Value(routing.NextVar(node))
- print('veh '+str(veh_number)+'\'s load:', vload)
- route.append(customers.customers[manager.IndexToNode(node)])
- #print("vehicle:",node,"load:",vload) # new
- return route
- else:
- return None
-
-
- def plot_vehicle_routes(veh_route, ax1, customers, vehicles):
- """
- Plot the vehicle routes on matplotlib axis ax1.
-
- Args: veh_route (dict): a dictionary of routes keyed by vehicle idx. ax1
- (matplotlib.axes._subplots.AxesSubplot): Matplotlib axes customers
- (Customers): the customers instance. vehicles (Vehicles): the vehicles
- instance.
- """
- veh_used = [v for v in veh_route if veh_route[v] is not None]
-
- cmap = discrete_cmap(vehicles.number + 2, 'nipy_spectral')
-
- for veh_number in veh_used:
-
- lats, lons = zip(*[(c.lat, c.lon) for c in veh_route[veh_number]])
- lats = np.array(lats)
- lons = np.array(lons)
- s_dep = customers.customers[vehicles.starts[veh_number]]
- s_fin = customers.customers[vehicles.ends[veh_number]]
- ax1.annotate(
- 'v({veh}) S @ {node}'.format(
- veh=veh_number, node=vehicles.starts[veh_number]),
- xy=(s_dep.lon, s_dep.lat),
- xytext=(10, 10),
- xycoords='data',
- textcoords='offset points',
- arrowprops=dict(
- arrowstyle='->',
- connectionstyle='angle3,angleA=90,angleB=0',
- shrinkA=0.05),
- )
- ax1.annotate(
- 'v({veh}) F @ {node}'.format(
- veh=veh_number, node=vehicles.ends[veh_number]),
- xy=(s_fin.lon, s_fin.lat),
- xytext=(10, -20),
- xycoords='data',
- textcoords='offset points',
- arrowprops=dict(
- arrowstyle='->',
- connectionstyle='angle3,angleA=-90,angleB=0',
- shrinkA=0.05),
- )
- ax1.plot(lons, lats, 'o', mfc=cmap(veh_number + 1))
- ax1.quiver(
- lons[:-1],
- lats[:-1],
- lons[1:] - lons[:-1],
- lats[1:] - lats[:-1],
- scale_units='xy',
- angles='xy',
- scale=1,
- color=cmap(veh_number + 1))
-
-
- def main():
- starttime = datetime.datetime.now()
- print(starttime)
- # 读取csv至字典
- df = pd.read_csv("./coldchain_20210730141118.csv",encoding='gbk')
- # 删去无效字段和缺失数据
- del df['sender_name']
- df = df.dropna(axis=0,how='any')
- # 统一体积单位为立方米
- df["schedule_total_volume"] = df["schedule_total_volume"] / 10**6
- df["goods_volume"] = df["goods_volume"] / 10**6
-
- # 历史派送地址集合,生成订单时从中随机有放回地抽取。
- # 当中有重复的,主要和自营还是外单有关,频率作为概率。
- locs = [(lon, lat) for lon, lat in zip(df["end_address_lng"], df["end_address_lat"])]
-
- # 100辆车
- num_vehicle = 100
- total_load = 12000 * num_vehicle * 0.65 # 12方容量车型,装载率65%
-
- num_stops = 1
- sum_load = 0
- demands = [0]
- lons = [116.596164]
- lats = [39.706433]
- tw_opens = [0]
- tw_closes = [0]
- while sum_load < total_load:
- num_stops += 1
- # 随机生成地址
- lon, lat = locs[np.random.randint(2766)]
- lats.append(lat)
- lons.append(lon)
- # 随机生成货量,80%概率属于head,20%属于long-tail,分别用两个指数分布来进行刻画
- if np.random.uniform(low = 0, high = 1) < 0.8:
- load = np.random.exponential(scale = 1/80) * 1000
- else:
- load = (np.random.exponential(scale = 0.8) + 0.1) * 1000
- demands.append(load)
- sum_load += load
- tw_opens.append(9 * 60**2)
- tw_closes.append(21 * 60**2)
-
- print("num_stops:",num_stops)
- # cust = Customers(num_stops, demands, lats, lons, tw_open, tw_close)
- # 设定目标
- aim_rate = 0.8
- aim_num_veh = int(np.ceil(total_load / (aim_rate * 12000) ))
-
- customers = Customers(num_stops, demands, lats, lons, tw_opens, tw_closes)
- capacity = 12000
- cost = 1
- vehicles = Vehicles(capacity=capacity, cost=cost, number = aim_num_veh)
-
- # Set the starting nodes, and create a callback fn for the starting node.
- start_fn = vehicles.return_starting_callback(
- customers, sameStartFinish=False)
-
- # Create the routing index manager.
- manager = pywrapcp.RoutingIndexManager(
- customers.number, # int number
- vehicles.number, # int number
- vehicles.starts, # List of int start depot
- vehicles.ends) # List of int end depot
- customers.set_manager(manager)
-
- # Set model parameters
- model_parameters = pywrapcp.DefaultRoutingModelParameters()
-
- # The solver parameters can be accessed from the model parameters. For example :
- # model_parameters.solver_parameters.CopyFrom(
- # pywrapcp.Solver.DefaultSolverParameters())
- # model_parameters.solver_parameters.trace_propagation = True
-
- # Make the routing model instance.
- routing = pywrapcp.RoutingModel(manager, model_parameters)
-
- parameters = pywrapcp.DefaultRoutingSearchParameters()
- # Setting first solution heuristic (cheapest addition).
- parameters.first_solution_strategy = (
- routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC)
- # Routing: forbids use of TSPOpt neighborhood, (this is the default behaviour)
- parameters.local_search_operators.use_tsp_opt = pywrapcp.BOOL_TRUE
- # Disabling Large Neighborhood Search, (this is the default behaviour)
- parameters.local_search_operators.use_path_lns = pywrapcp.BOOL_TRUE
- parameters.local_search_operators.use_inactive_lns = pywrapcp.BOOL_TRUE
-
- parameters.solution_limit = 10
- parameters.use_full_propagation = True
- #parameters.log_search = True
-
- # Create callback fns for distances, demands, service and transit-times.
- dist_fn = customers.return_dist_callback()
- dist_fn_index = routing.RegisterTransitCallback(dist_fn)
-
- dem_fn = customers.return_dem_callback()
- dem_fn_index = routing.RegisterUnaryTransitCallback(dem_fn)
-
- # Create and register a transit callback.
- serv_time_fn = customers.make_service_time_call_callback()
- transit_time_fn = customers.make_transit_time_callback()
- def tot_time_fn(from_index, to_index):
- """
- The time function we want is both transit time and service time.
- """
- # Convert from routing variable Index to distance matrix NodeIndex.
- from_node = manager.IndexToNode(from_index)
- to_node = manager.IndexToNode(to_index)
- return serv_time_fn(from_node, to_node) + transit_time_fn(from_node, to_node)
-
- tot_time_fn_index = routing.RegisterTransitCallback(tot_time_fn)
-
- # Set the cost function (distance callback) for each arc, homogeneous for
- # all vehicles.
- routing.SetArcCostEvaluatorOfAllVehicles(dist_fn_index)
-
- # Set vehicle costs for each vehicle, not homogeneous.
- for veh in vehicles.vehicles:
- routing.SetFixedCostOfVehicle(int(veh.cost), int(veh.index))
-
- # Add a dimension for vehicle capacities
- null_capacity_slack = 0
- capacities = [capacity for i in range(0,aim_num_veh)]
- routing.AddDimensionWithVehicleCapacity(
- dem_fn_index, # demand callback
- null_capacity_slack,
- capacities, # capacity array
- True,
- 'Capacity')
- # Add a dimension for time and a limit on the total time_horizon
- routing.AddDimension(
- tot_time_fn_index, # total time function callback
- customers.time_horizon,
- customers.time_horizon,
- True,
- 'Time')
-
- time_dimension = routing.GetDimensionOrDie('Time')
- for cust in customers.customers:
- if cust.tw_open is not None:
- time_dimension.CumulVar(manager.NodeToIndex(cust.index)).SetRange(
- cust.tw_open.seconds, cust.tw_close.seconds)
-
- """
- To allow the dropping of orders, we add disjunctions to all the customer
- nodes. Each disjunction is a list of 1 index, which allows that customer to
- be active or not, with a penalty if not. The penalty should be larger
- than the cost of servicing that customer, or it will always be dropped!
- """
- # To add disjunctions just to the customers, make a list of non-depots.
- non_depot = set(range(customers.number))
- non_depot.difference_update(vehicles.starts)
- non_depot.difference_update(vehicles.ends)
- # penalty = 400000 # The cost for dropping a node from the plan.
- # nodes = [routing.AddDisjunction([manager.NodeToIndex(c)], penalty) for c in non_depot]
-
- # This is how you would implement partial routes if you already knew part
- # of a feasible solution for example:
- # partial = np.random.choice(list(non_depot), size=(4,5), replace=False)
-
- # routing.CloseModel()
- # partial_list = [partial[0,:].tolist(),
- # partial[1,:].tolist(),
- # partial[2,:].tolist(),
- # partial[3,:].tolist(),
- # [],[],[],[]]
- # print(routing.ApplyLocksToAllVehicles(partial_list, False))
-
- # Solve the problem !
- assignment = routing.SolveWithParameters(parameters)
-
- # The rest is all optional for saving, printing or plotting the solution.
- if assignment:
- ## save the assignment, (Google Protobuf format)
- #save_file_base = os.path.realpath(__file__).split('.')[0]
- #if routing.WriteAssignment(save_file_base + '_assignment.ass'):
- # print('succesfully wrote assignment to file ' + save_file_base +
- # '_assignment.ass')
-
- print('The Objective Value is {0}'.format(assignment.ObjectiveValue()))
-
- plan_output, dropped = vehicle_output_string(manager, routing, assignment)
- print(plan_output)
- print('dropped nodes: ' + ', '.join(dropped))
-
- # you could print debug information like this:
- # print(routing.DebugOutputAssignment(assignment, 'Capacity'))
-
- vehicle_routes = {}
- for veh in range(vehicles.number):
- vehicle_routes[veh] = build_vehicle_route(manager, routing, assignment,
- customers, veh, num_stops)
- '''
- for j in delivered:
- if j in aim:
- aim.remove(j)
- else:
- print(str(j)+"is not in the aimed stops list")
- '''
-
- # Plotting of the routes in matplotlib.
- fig = plt.figure()
- ax = fig.add_subplot(111)
- # Plot all the nodes as black dots.
- clon, clat = zip(*[(c.lon, c.lat) for c in customers.customers])
- ax.plot(clon, clat, 'k.')
- # plot the routes as arrows
- plot_vehicle_routes(vehicle_routes, ax, customers, vehicles)
-
- endtime = datetime.datetime.now()
- print("程序运行时间:")
- print ((endtime - starttime))
- plt.show()
-
- #输出写入文件
- with open("data.txt","w") as f:
- f.write("Orders:"+'\n')
- for customer in customers.customers:
- f.write(str(customer)+'\n')
- f.write('\n'+"Plan_output:"+'\n')
- f.writelines(plan_output)
- f.write('\n'+"Dropped nodes:"+'\n')
- f.writelines(dropped)
-
-
- else:
- print('No assignment')
-
-
-
- if __name__ == '__main__':
- main()
-
|