#Suplementary Material 1 import pandas as pd import numpy as np import csv try: import docplex.mp except: raise Exception('Please install docplex. See https://pypi.org/project/docplex/') #load datasets #load excel file of stations with including ID and Number of ambulance vehicles def read_candidate_sites(candidates_csv_path): df = pd.read_csv(candidates_csv_path) df["NumberOfAmbulances"] = df["NumberOfAmbulances"].fillna(value=0) return df #load excel file of demand points including ID, Call frequancy, and location (urban/rural) def read_demands(demand_csv_path): df = pd.read_csv(demand_csv_path) df["CallFreq"] = df["CallFreq"].fillna(value=0) return df #load cost matrix including FacilityID, DemandID, DriveTime def create_OD_matrix(OD_csv_path, time_threshold): df = pd.read_csv(OD_csv_path) # determine the cover relationship based on the travel time and the service standard df["Covered"] = np.where(df["DriveTime"] < time_threshold, 1, 0) # create pivot table for OD matrix pivot = df.pivot("StationID", "DemandID", "Covered") return pivot # Scenario 0: Current distribution def calculate_objective_with_current_car_distribution(demand_csv, candidates_csv, ODMatrix_csv, time_threshold, unit_car_capacity): # read input data for the model demands = read_demands(demand_csv) candidates = read_candidate_sites(candidates_csv) coverage_matrix = create_OD_matrix(ODMatrix_csv, time_threshold) from docplex.mp.environment import Environment env = Environment() env.print_information() from docplex.mp.model import Model mdl = Model("EMS vehicles") # Define the decision variables # percentage of demand i covered by facility j y_i_j_vars = mdl.continuous_var_matrix(demands["DemandID"], candidates["StationID"], ub=1, name="y") # add constraints # ct1: the allocated demand should not exceed the capacity of the facility for j in candidates["StationID"]: mdl.add_constraint(mdl.scal_prod([y_i_j_vars[i, j] for i in demands["DemandID"]], demands["CallFreq"]) <= unit_car_capacity * candidates.loc[ candidates["StationID"] == j, "NumberOfAmbulances"].item()) # ct2: The allocated demand at i should not exceed 100% for i in demands["DemandID"]: mdl.add_constraint(mdl.sum(y_i_j_vars[i, j] for j in candidates["StationID"]) == 1) # express the objective total_covered_demand = mdl.sum(y_i_j_vars[i, j] * demands.loc[demands["DemandID"] == i, "CallFreq"].item() * coverage_matrix.loc[j, i] for i in demands["DemandID"] for j in candidates["StationID"]) mdl.maximize(total_covered_demand) mdl.print_information() # solve the model mdl.solve() # print the solution print("Total covered demand = %g" % mdl.objective_value) #Scenaro 1 and 2: Relocatoin and Allocation model def mcmclp_additive_model(demand_csv, candidates_csv, ODMatrix_csv, time_threshold, unit_car_capacity, maximal_cars_per_site,Output_csv, total_added_cars, additive_mode=0): # read input data for the model demands = read_demands(demand_csv) candidates = read_candidate_sites(candidates_csv) coverage_matrix = create_OD_matrix(ODMatrix_csv, time_threshold) max_matrix_rural = create_OD_matrix(ODMatrix_csv, max_time_rural) #Create an upper bound coverage matrix for rural demands max_matrix_urban = create_OD_matrix(ODMatrix_csv, max_time_urban) #Create an upper bound coverage matrix for urban demands from docplex.mp.environment import Environment env = Environment() # create a model from docplex.mp.model import Model mdl = Model("EMS vehicles") # define the decision variables # dv1: the percentage of demand i covered by facility j y_i_j_vars = mdl.continuous_var_matrix(demands["DemandID"], candidates["StationID"], ub=1, name="y") # dv2: the number of cars added at a station x_j_vars = mdl.integer_var_dict(candidates["StationID"], name="x") # add constraints # ct1: the allocated demand should not exceed the capacity of the facility # ct2: The total number of cars at each site should not exceed maximal_cars_per_site for j in candidates["StationID"]: num_existing_cars = candidates.loc[candidates["StationID"] == j, "NumberOfAmbulances"].item() if additive_mode == 0: num_existing_cars = 0 num_total_cars = num_existing_cars + x_j_vars[j] mdl.add_constraint(num_total_cars <= maximal_cars_per_site) mdl.add_constraint(mdl.scal_prod([y_i_j_vars[i, j] for i in demands["DemandID"]], demands["CallFreq"]) <= unit_car_capacity * num_total_cars) # ct3: the total number of EMS vehicles added should be equal to total_added_cars mdl.add_constraint(mdl.sum(x_j_vars[j] for j in candidates["StationID"]) == total_added_cars) # ct4: The allocated demand at i should not exceed 100% for i in demands["DemandID"]: mdl.add_constraint(mdl.sum(y_i_j_vars[i, j] for j in candidates["StationID"]) == 1) # ct5: Upper bound for urban and rural demands #NEW for i in demands.loc[demands["Urban"] == 1,"DemandID"]: mdl.add_constraint(mdl.sum(max_matrix_urban.loc[j,i]* x_j_vars[j] for j in candidates["StationID"]) >= 1) for i in demands.loc[demands["Urban"] == 0,"DemandID"]: mdl.add_constraint(mdl.sum(max_matrix_rural.loc[j,i]* x_j_vars[j] for j in candidates["StationID"]) >= 1) # express the objective total_covered_demand = mdl.sum(y_i_j_vars[i, j] * demands.loc[demands["DemandID"] == i, "CallFreq"].item() * coverage_matrix.loc[j, i] for i in demands["DemandID"] for j in candidates["StationID"]) mdl.maximize(total_covered_demand) mdl.print_information() # solve the model mdl.solve() # print the solution print("Total covered demand = %g" % mdl.objective_value) total_cars = 0 for j in candidates["StationID"]: num_added_cars = round(x_j_vars[j].solution_value) total_cars += num_added_cars if (num_added_cars>0): print("{0} #vehicles added: {1!s}".format(j, num_added_cars)) # print([y_i_j_vars[i, j].solution_value for i in demands["DemandID"]]) print("total cars added: {0!s}".format(total_cars)) ''' Run the models ''' # Model parameters input_data_folder = r"Folder_Path" demand_csv = input_data_folder + "\Demand.csv" candidates_csv = input_data_folder + "\Stations.csv" ODMatrix_csv = input_data_folder + "\ODMatrix.csv" Output_csv= input_data_folder + "\output.csv" time_threshold = 5 # minutes unit_car_capacity = 2387 # 224355/94 maximal_cars_per_site = 3 max_time_rural= 16 #minutes max_time_urban= 48 #minutes # use scenario variable to control which model will run # 0 - assess the covered demand based on the current distribution of EMS vehicles # 1 - scenario 1 where all cars can be relocated for optimization # 2 - scenario 2 where existing cars remain and added cars are optimally located scenario = 2 # run the model if scenario == 1: # run the model for scenario 1 and 2 mcmclp_additive_model(demand_csv, candidates_csv, ODMatrix_csv, time_threshold, unit_car_capacity, maximal_cars_per_site,Output_csv, total_added_cars=94, additive_mode=0) elif scenario == 2: # run the model for scenario 2 for i in range (0,11): #adding zero to 10 new vehicles mcmclp_additive_model(demand_csv, candidates_csv, ODMatrix_csv, time_threshold, unit_car_capacity, maximal_cars_per_site, Output_csv, total_added_cars=i, additive_mode=1) elif scenario == 0: # assess the covered demand based on the current distribution of EMS vehicles calculate_objective_with_current_car_distribution(demand_csv, candidates_csv, ODMatrix_csv, time_threshold, unit_car_capacity)