diff --git a/extras.py b/extras.py new file mode 100644 index 0000000000000000000000000000000000000000..ed5360cad0b8fc65e4c367551f78ed34eb6810e4 --- /dev/null +++ b/extras.py @@ -0,0 +1,128 @@ +import os +import json +import pandas as pd +import argparse +import random +import utilities, solve_ilp, solve_bap, prepare_instance + + +""" +COMPARE INSTANCES +""" +folder = 'instances/bug' +instances_path = [os.path.join(folder, f) for f in os.listdir(folder) if os.path.isfile(os.path.join(folder, f)) and f.endswith('.json')] + +names, ilp_obj, bap_obj = [], [], [] +for instance in instances_path: + print(instance) + # _, solution_ilp = solve_ilp.main(instance) + _, solution_bap, _ = solve_bap.main(instance, utilities.parse_arguments_from_file()) + names.append(instance) + # ilp_obj.append(round(solution_ilp.objective_value,2)) + bap_obj.append(round(solution_bap.objective_value,2)) + +# for i in range(len(names)): +# print(f"{names[i]}\t\t{ilp_obj[i]}\t\t{bap_obj[i]}\t\t{'NO' if ilp_obj[i] != bap_obj[i] else ''}") + + +""" +COMPUTE LOAD FACTOR +""" +# folder = 'instances/experiments_quadratic_2/bap_before' +# subfolders = [f.path for f in os.scandir(folder) if f.is_dir()] +# +# for subfolder_path in subfolders: +# instances = ['2_p60.json', '2_p100.json', '1_p100.json', '2_p80.json', '0_p100.json', '0_p80.json', '0_p60.json', +# '1_p80.json', '1_p60.json'] +# instances_path = [os.path.join(subfolder_path, i) for i in instances] +# print(subfolder_path) +# for instance_path in instances_path: +# instance = prepare_instance.main(instance_path) +# +# # load factor = total surgery time / available time +# load_factor = sum([patient.duration for patient in instance.patients]) / (instance.n_blocks * instance.blocks[0].capacity) +# print(round(load_factor, 2)) + + +""" +EDIT JSON (add feature) +""" +# folder = 'instances/ilp_after_overtime_undertime_sequence' +# subfolders = [f.path for f in os.scandir(folder) if f.is_dir()] +# json_files = [] +# for subfolder in subfolders: +# json_files.extend([f.path for f in os.scandir(subfolder) if f.path.endswith('.json')]) +# +# for path in json_files: +# with open(path, 'r') as f: +# data = json.load(f) +# data["sequence_duration"] = [ +# [0, 20, 25, 30, 35, 40], +# [20, 0, 20, 25, 30, 35], +# [25, 20, 0, 20, 25, 30], +# [30, 25, 20, 0, 20, 25], +# [35, 30, 25, 20, 0, 20], +# [40, 35, 30, 25, 20, 0] +# ] +# for i in data['patients']: +# i['group'] = random.randint(0, 5) +# with open(path, 'w') as f: +# json.dump(data, f) + +""" +EDIT CSV +""" +# path = 'instances/record_instances_2_ranking/instances_r1_s8_d8/subproblem_features_recording.csv' +# days = 8 +# batch_size = 2 * days +# +# # open CSV +# df = pd.read_csv(path) +# +# # prepare column vector with batches +# batches = [] +# print(len(df)/batch_size) +# for b in range(int(len(df)/batch_size)): +# for i in range(batch_size): +# batches.append(b) +# +# # add batches as new column +# df['batch index'] = batches +# +# # save CSV +# df.to_csv(path, index=False) + + +""" +MERGE CSV +""" +# # define folder to search in by parsing input arguments +# parser = argparse.ArgumentParser() +# parser.add_argument('--folder', '-f', help='path to the folder where the generated instances should be stored') +# args = parser.parse_args() +# folder = args.folder +# +# # 1. define final CSV file +# csv_name = 'subproblem_features_recording.csv' +# final_csv_path = os.path.join(folder, csv_name) +# +# # 2. list all subfolder's CSV files +# subfolders = [f.path for f in os.scandir(folder) if f.is_dir()] +# subfolder_csv_files = [os.path.join(subfolder_path, csv_name) for subfolder_path in subfolders if os.path.exists(os.path.join(subfolder_path, csv_name))] +# +# # 3. merge them together and save final CSV +# subfolder_dfs = [] +# batch_increase = 0 +# for f in subfolder_csv_files: +# df = pd.read_csv(f) +# +# # in merging datasets, get unique value of batch index, if the column is present +# if 'batch index' in df.columns: +# batch_increase_new = max(df['batch index']) + 1 +# df['batch index'] = df['batch index'] + batch_increase +# batch_increase += batch_increase_new +# +# subfolder_dfs.append(df) +# +# df_concat = pd.concat(subfolder_dfs, ignore_index=True) +# df_concat.to_csv(final_csv_path, index=False) diff --git a/generate_instances.py b/generate_instances.py index 432f892b6c06cf87af2ea893cbe2c5e77be69c11..57821044ea3c367f4821e55363322ff63abe64ec 100644 --- a/generate_instances.py +++ b/generate_instances.py @@ -63,12 +63,12 @@ def sample_surgeons(number_of_surgeons, number_of_blocks): def sample_patients(number_of_days, number_of_surgeons, number_of_patients): """ Samples patient object with its assignment to surgeon. - > fixed parameters: surgery duration range, patient priority + > fixed parameters: surgery duration range, patient priority, patient group :param number_of_days: number of days to randomly sample from :param number_of_surgeons: number of surgeons to randomly sample from :param number_of_patients: number of patients to randomly sample from - :return: sampled surgeon (with release day, surgery duration, patient priority and surgeon) + :return: sampled surgeon (with release day, surgery duration, patient priority, surgeon and patient's group) """ patients = [] @@ -88,7 +88,7 @@ def sample_instance(number_of_days, number_of_rooms, number_of_surgeons, number_ """ Samples instance object with all the sampled blocks, surgeons and patients. > fixed parameters: maximum blocks assigned to a surgeon per day, - weights m1, m5, m6, m7, rooms' overtime + weights m1, m5, m6, m7, rooms' overtime, sequence duration :param number_of_days: number of days to randomly sample from :param number_of_rooms: number of rooms to randomly sample from @@ -115,7 +115,7 @@ def sample_instance(number_of_days, number_of_rooms, number_of_surgeons, number_ "m6": 100, "m7": 1 }, - "rooms": [{"overtime": 120} for _ in range(number_of_rooms)], + "rooms": [{"overtime": 90} for _ in range(number_of_rooms)], "blocks": blocks, "surgeons": surgeons, "patients": patients, @@ -146,10 +146,10 @@ def main(path, number_of_instances): os.mkdir(path) # prepare constants - number_of_days_array = [10, 20] + number_of_days_array = [2,4,6] number_of_rooms_array = [1] number_of_surgeons_array = [4, 8] - number_of_patients_array = [60, 80, 100] + number_of_patients_array = [10,15,20] for number_of_days in number_of_days_array: for number_of_rooms in number_of_rooms_array: diff --git a/merge_csv.py b/merge_csv.py deleted file mode 100644 index 94bb5d74579f9446be0f0597e9f2c6d853019cb0..0000000000000000000000000000000000000000 --- a/merge_csv.py +++ /dev/null @@ -1,34 +0,0 @@ -import os -import pandas as pd -import argparse - -# define folder to search in by parsing input arguments -parser = argparse.ArgumentParser() -parser.add_argument('--folder', '-f', help='path to the folder where the generated instances should be stored') -args = parser.parse_args() -folder = args.folder - -# 1. define final CSV file -csv_name = 'subproblem_features_recording.csv' -final_csv_path = os.path.join(folder, csv_name) - -# 2. list all subfolder's CSV files -subfolders = [f.path for f in os.scandir(folder) if f.is_dir()] -subfolder_csv_files = [os.path.join(subfolder_path, csv_name) for subfolder_path in subfolders if os.path.exists(os.path.join(subfolder_path, csv_name))] - -# 3. merge them together and save final CSV -subfolder_dfs = [] -batch_increase = 0 -for f in subfolder_csv_files: - df = pd.read_csv(f) - - # in merging datasets, get unique value of batch index, if the column is present - if 'batch index' in df.columns: - batch_increase_new = max(df['batch index']) + 1 - df['batch index'] = df['batch index'] + batch_increase - batch_increase += batch_increase_new - - subfolder_dfs.append(df) - -df_concat = pd.concat(subfolder_dfs, ignore_index=True) -df_concat.to_csv(final_csv_path, index=False) diff --git a/parameters.yaml b/parameters.yaml index fda86ac74f19c2706eb4c730059a998ce2f83776..61d1dfa0b48ce33fc6ae3aa983909a4fbdfb1137 100644 --- a/parameters.yaml +++ b/parameters.yaml @@ -1,4 +1,4 @@ -threshold_depth: 20 +threshold_depth: 40 subproblem_strategy: 'return k best' sorting_strategy: 'random' model_path: '' diff --git a/prepare_bap.py b/prepare_bap.py index f9bd00691ac09dc522a05719b1d017638e8e0baa..75708054eaffa7564f6cc4396f94356354ecb840 100644 --- a/prepare_bap.py +++ b/prepare_bap.py @@ -88,12 +88,14 @@ class MasterProblemPrimal: o = self.model.addVars(self.instance.n_blocks, lb=0, ub=g.GRB.INFINITY, vtype=g.GRB.INTEGER, name="o") z = self.model.addVars(self.instance.n_blocks, lb=0, ub=g.GRB.INFINITY, vtype=g.GRB.INTEGER, name="z") n = self.model.addVars(self.instance.n_surgeons, self.instance.n_days, vtype=g.GRB.BINARY, name="n") + y = self.model.addVars(self.instance.n_patients + 1, self.instance.n_patients + 1, self.instance.n_blocks, vtype=g.GRB.BINARY, name="y") + e = self.model.addVars(self.instance.n_blocks, vtype=g.GRB.BINARY, name="e") # add constraints for patient in self.instance.patients: self.model.addConstr(g.LinExpr((1, x[patient.id, block.id]) for block in self.instance.blocks) <= 1, name=f"constraint_1_p{patient.id}") for block in self.instance.blocks: - self.model.addConstr(g.LinExpr((patient.duration, x[patient.id, block.id]) for patient in self.instance.patients) - block.capacity == o[block.id] - z[block.id], name=f"constraint_2_b{block.id}") + self.model.addConstr(g.LinExpr((patient.duration, x[patient.id, block.id]) for patient in self.instance.patients) + g.LinExpr((self.instance.sequence_duration[patient_1.group][patient_2.group], y[patient_1.id, patient_2.id, block.id]) for patient_1 in self.instance.patients for patient_2 in self.instance.patients if patient_2.id != patient_1.id) - block.capacity == o[block.id] - z[block.id], name=f"constraint_2_b{block.id}") self.model.addConstrs((o[block.id] <= block.overtime for block in self.instance.blocks), name="constraint_3") for room in self.instance.rooms: for day in self.instance.days: @@ -104,6 +106,16 @@ class MasterProblemPrimal: self.model.addConstrs((x[patient.id, block.id] == 0 for block in self.instance.blocks for patient in self.instance.patients if sum([1 if patient.surgeon == surgeon.id and block.id in surgeon.blocks else 0 for surgeon in self.instance.surgeons]) == 0), name="constraint_6") self.model.addConstrs((x[patient.id, block.id] == 0 for patient in self.instance.patients for block in self.instance.blocks if patient.release > block.day), name="constraint_7") + for patient in self.instance.patients: + for block in self.instance.blocks: + self.model.addConstr(x[patient.id, block.id] == g.LinExpr((1, y[patient.id, patient_next.id, block.id]) for patient_next in self.instance.patients if patient_next.id != patient.id) + y[patient.id, self.instance.n_patients, block.id]) + self.model.addConstr(x[patient.id, block.id] == g.LinExpr((1, y[patient_prev.id, patient.id, block.id]) for patient_prev in self.instance.patients if patient_prev.id != patient.id) + y[self.instance.n_patients, patient.id, block.id]) + self.model.addConstr(x[patient.id, block.id] <= e[block.id]) + for block in self.instance.blocks: + self.model.addConstr(g.LinExpr((1, y[patient.id, self.instance.n_patients, block.id]) for patient in self.instance.patients) == e[block.id]) + self.model.addConstr(g.LinExpr((1, y[self.instance.n_patients, patient.id, block.id]) for patient in self.instance.patients) == e[block.id]) + self.model.addConstr(e[block.id] <= g.LinExpr((1, x[patient.id, block.id]) for patient in self.instance.patients)) + # add surgeon-day branching constraints for branching_constraint in self.branching_constraints['surgeon']: branching_surgeon_id, branching_day_id, indication = branching_constraint[0], branching_constraint[1], branching_constraint[2] @@ -129,6 +141,7 @@ class MasterProblemPrimal: # optimize the model self.model.optimize() + print(f'BaP ILP model status: {self.model.status}', flush=True) # if the model is infeasible or cut off, exit the function if (self.model.status != g.GRB.OPTIMAL and self.model.status != g.GRB.TIME_LIMIT) or self.model.SolCount == 0: @@ -138,13 +151,13 @@ class MasterProblemPrimal: self.optimality_gap = self.model.MIPGap # store objective value - self.objective_value = round(float(self.model.objVal), 5) + self.objective_value = round(float(self.model.objVal), 1) # store all the model variables self.primal_variables = { 'x': utilities.convert_dictionary_to_array(self.model.getAttr('X', x), np.zeros((self.instance.n_patients, self.instance.n_blocks))), 'o': self.model.getAttr('X', o.values()), - 'z':self.model.getAttr('X', z.values()), + 'z': self.model.getAttr('X', z.values()), 'n': utilities.convert_dictionary_to_array(self.model.getAttr('X', n), np.zeros((self.instance.n_surgeons, self.instance.n_days))) } @@ -165,6 +178,11 @@ class MasterProblemPrimal: self.model.setParam(g.GRB.Param.Seed, 1) self.model.setParam(g.GRB.Param.TimeLimit, timelimit) + # TODO: tune tolerance parameters + # self.model.setParam(g.GRB.Param.Presolve, 0) + # self.model.setParam(g.GRB.Param.MIPGap, 1e-6) + # self.model.setParam(g.GRB.Param.IntFeasTol, 1e-6) + # create variables if integer: theta = self.model.addVars(len(self.instance.patterns), vtype=g.GRB.BINARY, name="theta") @@ -211,13 +229,14 @@ class MasterProblemPrimal: # optimize the model self.model.optimize() + print(f'BaP model status: {self.model.status}', flush=True) # if the model is infeasible or cut off, exit the function if (self.model.status != g.GRB.OPTIMAL and self.model.status != g.GRB.TIME_LIMIT) or self.model.SolCount == 0: return # store objective value - objective_value = round(float(self.model.objVal), 5) + objective_value = round(float(self.model.objVal), 1) # store optimality gap if integer: @@ -225,7 +244,7 @@ class MasterProblemPrimal: # store primal variables to dictionary primal_variables = { - 'theta': {idx: round(theta, 5) for idx, theta in zip([p.id for p in self.instance.patterns], self.model.getAttr('X', theta.values()))}, + 'theta': {idx: theta for idx, theta in zip([p.id for p in self.instance.patterns], self.model.getAttr('X', theta.values()))}, 'n': utilities.convert_dictionary_to_array(self.model.getAttr('X', n), np.zeros((self.instance.n_surgeons, self.instance.n_days))) } @@ -266,7 +285,7 @@ class MasterProblemPrimal: for surgeon in self.instance.surgeons: for day in self.instance.days: reduced_cost = self.model.getVarByName(f'n[{surgeon.id},{day}]').getAttr('RC') - variable_value = round(self.model.getVarByName(f'n[{surgeon.id},{day}]').X, 5) + variable_value = round(self.model.getVarByName(f'n[{surgeon.id},{day}]').X, 2) if z_lp + reduced_cost > z_ilp_best and variable_value == 0: fixing_combinations.append((surgeon.id, day)) @@ -280,7 +299,7 @@ class MasterProblemPrimal: """ is_fractional_patient, is_fractional_surgeon = False, False - epsilon = 1e-5 + epsilon = 1e-3 # check for a fractionality of patient for key, value in self.primal_variables['theta'].items(): @@ -322,7 +341,7 @@ class MasterProblemPrimal: for block in self.instance.blocks: patterns = [p for p in self.instance.patterns if p.block == block] for patient in self.instance.patients: - value = round(sum([(1 if patient.id in pattern.patient_id else 0) * self.primal_variables['theta'][pattern.id] for pattern in patterns]), 3) + value = round(sum([(1 if patient.id in pattern.patient_id else 0) * self.primal_variables['theta'][pattern.id] for pattern in patterns]), 2) if not value.is_integer(): variables.append([patient.id, block.id]) values.append(value) @@ -331,7 +350,7 @@ class MasterProblemPrimal: elif variable_type == 'surgeon': for surgeon in self.instance.surgeons: for day in self.instance.days: - value = round(self.primal_variables['n'][surgeon.id, day], 3) + value = round(self.primal_variables['n'][surgeon.id, day], 2) if not value.is_integer(): variables.append([surgeon.id, day]) values.append(value) @@ -387,23 +406,35 @@ class Subproblem: self.model = g.Model() self.model.setParam(g.GRB.Param.OutputFlag, 0) self.model.setParam(g.GRB.Param.Seed, 1) + self.model.setParam(g.GRB.Param.MIPGap, 1e-6) + self.model.setParam(g.GRB.Param.IntFeasTol, 1e-6) # create variables self.x = self.model.addVars(self.instance.n_patients, vtype=g.GRB.BINARY, name="x") self.o = self.model.addVar(lb=0, ub=g.GRB.INFINITY, vtype=g.GRB.INTEGER, name="o") self.z = self.model.addVar(lb=0, ub=g.GRB.INFINITY, vtype=g.GRB.INTEGER, name="z") self.n = self.model.addVars(self.instance.n_surgeons, vtype=g.GRB.BINARY, name="n") + self.y = self.model.addVars(self.instance.n_patients + 1, self.instance.n_patients + 1, vtype=g.GRB.BINARY, name="y") + self.e = self.model.addVar(vtype=g.GRB.BINARY, name="e") # add initial constraints - self.model.addConstr(g.LinExpr((patient.duration, self.x[patient.id]) for patient in self.instance.patients) - block.capacity == self.o - self.z, name="constraint_1") + self.model.addConstr(g.LinExpr((patient.duration, self.x[patient.id]) for patient in self.instance.patients) + g.LinExpr((self.instance.sequence_duration[patient_1.group][patient_2.group], self.y[patient_1.id, patient_2.id]) for patient_1 in self.instance.patients for patient_2 in self.instance.patients if patient_2.id != patient_1.id) - block.capacity == self.o - self.z, name=f"constraint_1") self.model.addConstr(self.o <= block.overtime, name="constraint_2") self.model.addConstrs((self.x[patient.id] <= self.n[surgeon.id] for surgeon in self.instance.surgeons if block.id in surgeon.blocks for patient in self.instance.patients if patient.surgeon == surgeon.id), name="constraint_3") self.model.addConstrs((self.x[patient.id] == 0 for patient in self.instance.patients if sum([1 if patient.surgeon == surgeon.id and block.id in surgeon.blocks else 0 for surgeon in self.instance.surgeons]) == 0), name="constraint_4") self.model.addConstrs((self.x[patient.id] == 0 for patient in self.instance.patients if patient.release > day), name="constraint_5") self.model.addConstrs((self.n[surgeon.id] == 0 for surgeon in self.instance.surgeons if block.id not in surgeon.blocks), name="constraint_6") + for patient in self.instance.patients: + self.model.addConstr(self.x[patient.id] == g.LinExpr((1, self.y[patient.id, patient_next.id]) for patient_next in self.instance.patients if patient_next.id != patient.id) + self.y[patient.id, self.instance.n_patients]) + self.model.addConstr(self.x[patient.id] == g.LinExpr((1, self.y[patient_prev.id, patient.id]) for patient_prev in self.instance.patients if patient_prev.id != patient.id) + self.y[self.instance.n_patients, patient.id]) + self.model.addConstr(self.x[patient.id] <= self.e) + self.model.addConstr(g.LinExpr((1, self.y[patient.id, self.instance.n_patients]) for patient in self.instance.patients) == self.e) + self.model.addConstr(g.LinExpr((1, self.y[self.instance.n_patients, patient.id]) for patient in self.instance.patients) == self.e) + self.model.addConstr(self.e <= g.LinExpr((1, self.x[patient.id]) for patient in self.instance.patients)) + # add patient-block branching constraints - for branching_constraint in branching_constraints: + for branching_constraint in branching_constraints['patient']: branching_patient_id, branching_block_id, indication = branching_constraint[0], branching_constraint[1], branching_constraint[2] if block.id == branching_block_id: # operating with block b if indication == 0: # 'patient p not assigned to block b' @@ -414,6 +445,12 @@ class Subproblem: if indication == 1: # 'patient p assigned to block b' self.model.addConstr((self.x[branching_patient_id] == 0), name=f"patient_branching_constraint_{branching_patient_id}_{branching_block_id}") + # add surgeon-day branching constraints (only n_sd = 1 as we cannot assign a surgeon to specific block, only day) + for branching_constraint in branching_constraints['surgeon']: + branching_surgeon_id, branching_day_id, indication = branching_constraint[0], branching_constraint[1], branching_constraint[2] + if indication == 0 and branching_day_id == day: + self.model.addConstr((self.n[branching_surgeon_id] == 0), name=f"surgeon_branching_constraint_{branching_surgeon_id}_{branching_day_id}") + # set objective part1 = g.LinExpr((master_problem_variables['lambda'][patient.id], self.x[patient.id]) for patient in self.instance.patients) \ + master_problem_variables['mu'][block.id] \ @@ -437,7 +474,7 @@ class Subproblem: return # store objective value - self.objective_value = round(float(self.model.objVal), 5) + self.objective_value = round(float(self.model.objVal), 1) # store optimized variables to dictionary self.variables = { diff --git a/prepare_instance.py b/prepare_instance.py index 225daefbdfaa71df43aab0c4d9e16a8729f87e7c..3562949b6eb3b37738662832a1461ce105397476 100644 --- a/prepare_instance.py +++ b/prepare_instance.py @@ -58,35 +58,39 @@ class Patient: """ A class used to represent one patient from the given set of patients. """ - def __init__(self, id, release, duration, priority, surgeon): + def __init__(self, id, release, duration, priority, surgeon, group): """ :param id: ID of the patient :param release: release time/date for the surgery of the patient :param duration: expected surgery duration of the patient :param priority: clinical priority coefficient of the patient :param surgeon: surgeon by which the patient needs to be treated + :param group: patient group to which the patient belongs to """ self.id = id self.release = release self.duration = duration self.priority = priority self.surgeon = surgeon + self.group = group class Instance: """ A class used to represent the instance on which the OR scheduling task should be performed. """ - def __init__(self, max_blocks_surgeon_day, weights): + def __init__(self, max_blocks_surgeon_day, weights, sequence_duration): """ :param max_blocks_surgeon_day: maximum number of blocks assigned to a surgeon in any day :param weights: dictionary of weights m_1,m_5,m_6,m_7 of term m_j in the objective function + :param sequence_duration: matrix defining sequence dependent setup time """ self.rooms, self.blocks, self.surgeons, self.patients, self.days = [], [], [], [], [] self.n_rooms, self.n_blocks, self.n_surgeons, self.n_patients, self.n_days = 0, 0, 0, 0, 0 self.max_blocks_surgeon_day = max_blocks_surgeon_day self.max_patients_surgeon_day = {} self.weights = weights + self.sequence_duration = sequence_duration self.patterns = [] def count_items(self): @@ -161,7 +165,7 @@ class Instance: # print patients and their surgeries for patient in self.patients: - print(f"PATIENT: {patient.id} (release: {patient.release}, duration: {patient.duration}, priority: {patient.priority}, surgeon: {patient.surgeon})") + print(f"PATIENT: {patient.id} (release: {patient.release}, duration: {patient.duration}, priority: {patient.priority}, surgeon: {patient.surgeon}), group: {patient.group}") print(f"\n{'-' * 150}\n") @@ -181,7 +185,7 @@ def convert_json_to_instance(json): """ # prepare Instance object - instance = Instance(json['max_blocks_surgeon_day'], json['weights']) + instance = Instance(json['max_blocks_surgeon_day'], json['weights'], json['sequence_duration']) # add all elements of the JSON file to the Instance object for key, item in enumerate(json['rooms']): @@ -191,7 +195,7 @@ def convert_json_to_instance(json): for key, item in enumerate(json['surgeons']): instance.surgeons.append(Surgeon(key, item['blocks'])) for key, item in enumerate(json['patients']): - instance.patients.append(Patient(key, item['release'], item['duration'], item['priority'], item['surgeon'])) + instance.patients.append(Patient(key, item['release'], item['duration'], item['priority'], item['surgeon'], item['group'])) for day in range(json['number_of_days']): instance.days.append(day) diff --git a/scripts/run.sh b/scripts/run.sh new file mode 100644 index 0000000000000000000000000000000000000000..e171417b197138259c017288befb99421f4ca12b --- /dev/null +++ b/scripts/run.sh @@ -0,0 +1,17 @@ +#!/bin/bash + +#SBATCH --cpus-per-task=4 +#SBATCH --mem-per-cpu=16G +#SBATCH --time=4-00:00:00 +#SBATCH --partition=compute,bigmem + +module purge +module load Python/3.7.2-GCCcore-8.2.0 +module load Gurobi/9.0.0 +module load SciPy-bundle/2019.10-fosscuda-2019b-Python-3.7.4 +pip3 install lightgbm --user +pip3 install PyYAML --user + +cd /home/koutepa2/advance-or-scheduling-with-ml/ +python validate_instances.py -f $1 + diff --git a/scripts/run_experiment.sh b/scripts/run_experiment.sh new file mode 100644 index 0000000000000000000000000000000000000000..d240e251cc9d521f1c014826ca50fccc8171b2a5 --- /dev/null +++ b/scripts/run_experiment.sh @@ -0,0 +1,15 @@ +#!/bin/bash + +# define folder +NAME="bap_thresh40_after_overtime_undertime_sequence" +FOLDER="/home/koutepa2/advance-or-scheduling-with-ml/instances/experiments_quadratic/$NAME" + +# search through all the subfolders in the defined folder +for SUBDIR in $(find $FOLDER -mindepth 1 -maxdepth 1 -type d); do + + echo "Running validation for $SUBDIR" + + # run job for each of the subfolders + sbatch --job-name="$NAME" --output="$SUBDIR/output.out" --error="$SUBDIR/output.err" scripts/run.sh "$SUBDIR" + +done diff --git a/scripts/run_single.sh b/scripts/run_single.sh new file mode 100644 index 0000000000000000000000000000000000000000..2afba03ad7984281ee257a470a0cc3b746ad1393 --- /dev/null +++ b/scripts/run_single.sh @@ -0,0 +1,22 @@ +#!/bin/bash + +#SBATCH --job-name=single +#SBATCH --output="/home/koutepa2/advance-or-scheduling-with-ml/instances/bug/output.out" +#SBATCH --error="/home/koutepa2/advance-or-scheduling-with-ml/instances/bug/output.err" +#SBATCH --cpus-per-task=4 +#SBATCH --mem-per-cpu=16G +#SBATCH --time=4-00:00:00 +#SBATCH --partition=compute,bigmem + +module purge +module load Python/3.7.2-GCCcore-8.2.0 +module load Gurobi/9.0.0 +module load SciPy-bundle/2019.10-fosscuda-2019b-Python-3.7.4 +pip3 install lightgbm --user +pip3 install PyYAML --user + +FOLDER="/home/koutepa2/advance-or-scheduling-with-ml/instances/bug/" +cd /home/koutepa2/advance-or-scheduling-with-ml +python validate_instances.py -f "$FOLDER" + + diff --git a/scripts/template.sh b/scripts/template.sh new file mode 100644 index 0000000000000000000000000000000000000000..e378420aa8269dee438cfbfb56b2685be0300e80 --- /dev/null +++ b/scripts/template.sh @@ -0,0 +1,21 @@ +#!/bin/bash +#SBATCH --job-name=t1400_30h_throughput +#SBATCH --output=/home/suchap/dev/Prevedig/results/out/t1400_30h_throughput.out +#SBATCH --error=/home/suchap/dev/Prevedig/results/out/t1400_30h_throughput.err +#SBATCH --cpus-per-task=16 +#SBATCH --mem=80G +#SBATCH --time=2-00:00:00 +#SBATCH --mail-user=suchap@cvut.cz +#SBATCH --partition=compute +#SBATCH --mail-type=ALL + +#module load Anaconda3/5.3.0 +module load Python/3.7.2-GCCcore-8.2.0 +module load Gurobi/9.0.0 + +pip3 install numpy --user + +#. /opt/apps/software/Anaconda3/5.3.0/etc/profile.d/conda.sh + +cd /home/suchap/dev/Prevedig/python_scripts/prevedig/ +python ILP_biochem_throughput_maxim_fixed_transport.py 1400 108000 diff --git a/scripts/template_detailed.sh b/scripts/template_detailed.sh new file mode 100644 index 0000000000000000000000000000000000000000..7d361cc74faa5638613cf09913b121c1f709ad6a --- /dev/null +++ b/scripts/template_detailed.sh @@ -0,0 +1,50 @@ +#!/bin/bash -l + +# Give your job a name, so you can recognize it in the queue overview +#SBATCH --job-name=test_or + +# Define, how many nodes you need. +#SBATCH --nodes=2 +#SBATCH --ntasks=4 # ask for 4 cpus + +# Define, how long the job will run in real time. +# d-hh:mm:ss +#SBATCH --time=2-00:00:00 + +# Define the partition on which the job shall run. May be omitted. +#SBATCH --partition gpu + +# How much memory you need. Choose one of those. +#SBATCH --mem-per-cpu=5000MB # defines memory per CPU/core +##SBATCH --mem=5GB # defines memory per node + +# Turn on mail notification. +#SBATCH --mail-type=END,FAIL + +# You may not place any commands before the last SBATCH directive + +# Define and create a unique scratch directory for this job +# /lscratch is local ssd disk on particular node which is faster than your network home dir +SCRATCH_DIRECTORY=/lscratch/${USER}/${SLURM_JOBID}.stallo-adm.uit.no +mkdir -p ${SCRATCH_DIRECTORY} +cd ${SCRATCH_DIRECTORY} + +# You can copy everything you need to the scratch directory +# ${SLURM_SUBMIT_DIR} points to the path where this script was submitted from +cp ${SLURM_SUBMIT_DIR}/myfiles*.txt ${SCRATCH_DIRECTORY} + +# This is where the actual work is done. In this case, the script only waits. +# The time command is optional, but it may give you a hint on how long the command worked. +time sleep 150 +#sleep 150 + +# After the job is done we copy our output back to $SLURM_SUBMIT_DIR +cp ${SCRATCH_DIRECTORY}/my_output ${SLURM_SUBMIT_DIR} + +# After everything is saved to the home directory, delete the work directory to +# save space on /lscratch +cd ${SLURM_SUBMIT_DIR} +rm -rf ${SCRATCH_DIRECTORY} + +# Finish the script +exit 0 diff --git a/solve_bap.py b/solve_bap.py index 3a79e80e22a619020d86abfb27339db731c9edb6..8ec7e25d0b5a381fdda3688ef8eb3e17b2684ec2 100644 --- a/solve_bap.py +++ b/solve_bap.py @@ -1,5 +1,5 @@ """ -In this file, branch and price algorithm is solved using only primal formulation. It consist of preparing the instance using +In this file, branch and price algorithm is solved using only primal formulation. It consists of preparing the instance using prepare_instance.py script, constructing initial set of patterns and running column generation and branching in a loop till the incumbent solution is found. """ @@ -18,7 +18,6 @@ import prepare_instance, prepare_bap, utilities params = utilities.Parameters() TIMELIMIT = 4 * 60 * 60 - def generate_initial_patterns(instance): """ This function generates initial set of feasible patterns for the column generation algorithm. @@ -55,8 +54,10 @@ def generate_initial_patterns(instance): # if we can still fit the patient into the block, do it for block in patient_blocks: - if block_vocab[block.id]['leftover_capacity'] >= patient.duration: - block_vocab[block.id]['leftover_capacity'] -= patient.duration + sequence_duration = 0 if block_vocab[block.id]['patients'] == [] else instance.sequence_duration[block_vocab[block.id]['patients'][-1].group][patient.group] + operation_duration = patient.duration + sequence_duration + if block_vocab[block.id]['leftover_capacity'] >= operation_duration: + block_vocab[block.id]['leftover_capacity'] -= operation_duration block_vocab[block.id]['patients'].append(patient) block_vocab[block.id]['surgeons'].add(patient_surgeon) break @@ -78,7 +79,7 @@ def generate_initial_patterns(instance): return patterns -def generate_improving_pattern(instance, primal_variables, dual_variables, branching_constraints): +def generate_improving_pattern(master_problem): """ This function generates improving pattern/s (= pattern with positive reduced cost) for the column generation algorithm. If the objective value of the solved subproblem is positive, we have found @@ -87,10 +88,7 @@ def generate_improving_pattern(instance, primal_variables, dual_variables, branc subproblems) at the beginning as well. Subproblems are then searched in the order predicted by the selected strategy. - :param instance: instance on which the OR scheduling task is performed - :param branching_constraints: branching constraints present in the current node - :param primal_variables: primal variables outputted by the master problem - :param dual_variables: dual variables outputted by the master problem + :param master_problem: current state of master problem :return: improving pattern/s found by the selected strategy """ @@ -98,26 +96,24 @@ def generate_improving_pattern(instance, primal_variables, dual_variables, branc # sort the block in order in which they should be searched if params.sorting_strategy == 'random': - sorted_blocks, sorted_predictions = random.sample(instance.blocks, len(instance.blocks)), np.zeros(len(instance.blocks)) + sorted_blocks, sorted_predictions = random.sample(master_problem.instance.blocks, len(master_problem.instance.blocks)), np.zeros(len(master_problem.instance.blocks)) if params.sorting_strategy == 'ranked': - sorted_blocks, sorted_predictions = rank_subproblems(instance, branching_constraints, primal_variables, dual_variables) - if params.sorting_strategy == 'dummy ranked': - sorted_blocks, sorted_predictions = rank_subproblems_dummy(instance, branching_constraints, dual_variables) + sorted_blocks, sorted_predictions = rank_subproblems(master_problem.instance, master_problem.branching_constraints, master_problem.primal_variables, master_problem.dual_variables) # if recording is performed, write down subproblem features if params.type_of_validation == 'record': - record_subproblems(instance, branching_constraints, primal_variables, dual_variables) + record_subproblems(master_problem.instance, master_problem.branching_constraints, master_problem.primal_variables, master_problem.dual_variables) # solve subproblem for (block, room, day)... for block, prediction in zip(sorted_blocks, sorted_predictions): - room, day = instance.get_room(block.room), block.day + room, day = master_problem.instance.get_room(block.room), block.day # update number of subproblems value for analysis params.analysis.number_of_subproblems += 1 # solve the subproblem for the pattern - subproblem = prepare_bap.Subproblem(instance) - subproblem.optimize_model(max(TIMELIMIT - (time.time() - params.analysis.time_start), 0), dual_variables, branching_constraints['patient'], block, day, room) + subproblem = prepare_bap.Subproblem(master_problem.instance) + subproblem.optimize_model(max(TIMELIMIT - (time.time() - params.analysis.time_start), 0), master_problem.dual_variables, master_problem.branching_constraints, block, day, room) # if the SP is not solved to optimality, continue with different subproblem if subproblem.model.status != g.GRB.OPTIMAL: @@ -131,11 +127,11 @@ def generate_improving_pattern(instance, primal_variables, dual_variables, branc if subproblem.objective_value > 0: # crete Pattern corresponding to solved subproblem - improving_item_patients = [instance.get_patient(idx) for idx, item in enumerate(subproblem.variables['x']) if round(item, 5) == 1] - improving_item_surgeons = [instance.get_surgeon(idx) for idx, item in enumerate(subproblem.variables['n']) if round(item, 5) == 1] - improving_pattern = prepare_bap.Pattern(max([pattern.id for pattern in instance.patterns + improving_patterns]) + 1 if instance.patterns + improving_patterns else 0, + improving_item_patients = [master_problem.instance.get_patient(idx) for idx, item in enumerate(subproblem.variables['x']) if round(item, 2) == 1] + improving_item_surgeons = [master_problem.instance.get_surgeon(idx) for idx, item in enumerate(subproblem.variables['n']) if round(item, 2) == 1] + improving_pattern = prepare_bap.Pattern(max([pattern.id for pattern in master_problem.instance.patterns + improving_patterns]) + 1 if master_problem.instance.patterns + improving_patterns else 0, block, improving_item_patients, improving_item_surgeons, subproblem.variables['o'], subproblem.variables['z'], - instance.n_blocks, instance.n_patients, instance.n_surgeons) + master_problem.instance.n_blocks, master_problem.instance.n_patients, master_problem.instance.n_surgeons) improving_patterns.append(improving_pattern) improving_values.append(subproblem.objective_value) @@ -184,15 +180,15 @@ def record_subproblem(instance, branching_constraints, primal_variables, dual_va # -- patients and patterns that are fully involved in the block (integral) and surgeons that are fully involved in the day sumproducts = {patient.id: sum([(1 if patient.id in pattern.patient_id else 0) * primal_variables['theta'][pattern.id] for pattern in possible_patterns]) for patient in instance.patients} - fully_involved_patients = [patient for patient in instance.patients if round(sumproducts[patient.id], 5) == 1] - # fully_involved_patterns = [pattern for pattern in possible_patterns if round(primal_variables['theta'][pattern.id], 5) >= 1] - # fully_involved_surgeons = [surgeon for surgeon in instance.surgeons if round(primal_variables['n'][surgeon.id, day], 5) >= 1] + fully_involved_patients = [patient for patient in instance.patients if round(sumproducts[patient.id], 2) == 1] + # fully_involved_patterns = [pattern for pattern in possible_patterns if round(primal_variables['theta'][pattern.id], 2) >= 1] + # fully_involved_surgeons = [surgeon for surgeon in instance.surgeons if round(primal_variables['n'][surgeon.id, day], 2) >= 1] # # -- is the schedule w.r.t. the block integral or not? # possible_patterns_theta = [primal_variables['theta'][pattern.id] for pattern in possible_patterns] # is_block_integral = 1 # for p_value in possible_patterns_theta: - # if not round(p_value, 5).is_integer(): + # if not round(p_value, 2).is_integer(): # is_block_integral = 0 # break @@ -254,7 +250,7 @@ def record_subproblems(instance, branching_constraints, primal_variables, dual_v # solve the subproblem for the pattern subproblem = prepare_bap.Subproblem(instance) - subproblem.optimize_model(max(TIMELIMIT - (time.time() - params.analysis.time_start), 0), dual_variables, branching_constraints['patient'], block, day, room) + subproblem.optimize_model(max(TIMELIMIT - (time.time() - params.analysis.time_start), 0), dual_variables, branching_constraints, block, day, room) # write down subproblem features subproblem_features = record_subproblem(instance, branching_constraints, primal_variables, dual_variables, block, day, subproblem.objective_value) @@ -282,46 +278,6 @@ def record_subproblems(instance, branching_constraints, primal_variables, dual_v params.analysis.subproblem_features.append(f) -def rank_subproblems_dummy(instance, branching_constraints, dual_variables): - - blocks, predictions, randomized_predictions = [], [], [] - - # obtain randomized predictions for all the blocks - for block in instance.blocks: - room, day = instance.get_room(block.room), block.day - - # solve the subproblem for the pattern - subproblem = prepare_bap.Subproblem(instance) - subproblem.optimize_model(max(TIMELIMIT - (time.time() - params.analysis.time_start), 0), dual_variables, branching_constraints['patient'], block, day, room) - - # if the SP is not solved to optimality, continue with different subproblem - if subproblem.model.status != g.GRB.OPTIMAL: - continue - - # write down the prediction - prediction = 1 if subproblem.objective_value > 0 else 0 - blocks.append(block.id) - predictions.append(prediction) - - return predictions - - # threshold = 1 - # random_number = random.uniform(0, 1) - # randomized_prediction = prediction if random_number < threshold else 1 - prediction - # randomized_predictions.append(randomized_prediction) - # - # # sort the blocks - # sorted_blocks = np.argsort(randomized_predictions)[::-1] - # - # # retrieve objects for the sorted blocks - # sorted_blocks_objects, sorted_blocks_predictions = [], [] - # for block_id in sorted_blocks: - # sorted_blocks_objects.append(instance.get_block(block_id)) - # sorted_blocks_predictions.append(randomized_predictions[block_id]) - # - # return sorted_blocks_objects, sorted_blocks_predictions - - def rank_subproblems(instance, branching_constraints, primal_variables, dual_variables): """ This function extracts features of a subproblem for each block in the current state (branching constrains, primal and @@ -390,6 +346,10 @@ def column_generation(master_problem): mp_time_start = time.time() + # if the timelimit is reached, kill the CG procedure + if time.time() - params.analysis.time_start > TIMELIMIT: + break + # 1. pass improving patterns from last iteration subproblem to restricted MP, define it and solve it to optimize dual variables print("Optimizing master problem...") master_problem.optimize_model(timelimit=max(TIMELIMIT - (time.time() - params.analysis.time_start), 0), integer=False) @@ -398,16 +358,22 @@ def column_generation(master_problem): params.analysis.mp_cpu_time += time.time() - mp_time_start params.analysis.number_of_cg_iterations += 1 + # if the model is infeasible, cut off or time limited do not look for improving patterns + if (master_problem.model.status != g.GRB.OPTIMAL and master_problem.model.status != g.GRB.TIME_LIMIT) or master_problem.model.SolCount == 0: + print("Warning: MP model in CG is infeasible or cut off.") + break + # 2. pass variables from restricted MP to subproblem, define it and solve it to find improving patterns print("Generating improving patterns...") sp_time_start = time.time() sp_before = params.analysis.number_of_subproblems - improving_patterns = generate_improving_pattern(master_problem.instance, master_problem.primal_variables, master_problem.dual_variables, master_problem.branching_constraints) + improving_patterns = generate_improving_pattern(master_problem) params.analysis.sp_cpu_time += time.time() - sp_time_start number_of_subproblems_array.append(params.analysis.number_of_subproblems - sp_before) if not improving_patterns: + print("No improving patterns, ending the column generation.") break # 3. if improving pattern (pattern with positive reduced cost) is found, add it to set of MP patterns and return to step 1, otherwise stop the process @@ -498,7 +464,7 @@ def branching(master_problem, master_problem_index, subject, horizon, indication # --> dummy pattern should be substituted with minimal one elif p_indication == 1: - # remove current dummy pattern for block b and substitue it with minimal one + # remove current dummy pattern for block b and substitute it with minimal one dummy_pattern = branched_master_problem.instance.get_pattern(block) patient_object = branched_master_problem.instance.get_patient(patient) extended_patients, extended_surgeons = dummy_pattern.patients + [patient_object], dummy_pattern.surgeons + [branched_master_problem.instance.get_surgeon(patient_object.surgeon)] @@ -577,7 +543,10 @@ def branch_and_price(instance): stack = [] # generate initial set of feasible patterns according to provided instance and add them to the instance + print("Generating initial feasible patterns...") instance.patterns = generate_initial_patterns(instance) + for p in instance.patterns: + p.print_pattern() # initialize root master problem and add it to the queue master_problem = prepare_bap.MasterProblemPrimal(instance, {'patient': [], 'surgeon': [], 'aggregated_surgeon': [], 'extended_aggregated_surgeon': []}, 0) @@ -618,11 +587,12 @@ def branch_and_price(instance): # make a note that ILP was performed in this node params.graph.update_node(str(master_problem.idx), {'method': 'ILP'}) - # if the ILP is infeasible or cut off because of the objective value, solution cannot be better, continue + # if the ILP is infeasible or cut off because of the objective value, solution cannot be better, continue with next node if (master_problem.model.status != g.GRB.OPTIMAL and master_problem.model.status != g.GRB.TIME_LIMIT) or master_problem.model.SolCount == 0: + print("ILP model is infeasible or cut off, solution cannot be better, fathoming current node and moving to the next one.",flush=True) continue - # check if the solution is better than optimal solution - if so, update it and continue + # check if the solution is better than optimal solution - if so, update it and continue with next node if master_problem.objective_value < params.solution.objective_value: print(f"Solution of ILP ({master_problem.objective_value}) is better than the incumbent solution ({params.solution.objective_value}).",flush=True) print("No fractional variable in the primal master problem revealed, solution is integer, updating the incumbent solution, fathoming current node and moving to the next one.",flush=True) @@ -635,25 +605,28 @@ def branch_and_price(instance): 'total_patterns': str(number_of_patterns_before), 'added_patterns': str(len(master_problem.instance.patterns) - number_of_patterns_before) }) - params.solution.update('original', master_problem.idx, master_problem.objective_value, master_problem.optimality_gap, master_problem.primal_variables, master_problem.instance.patterns, time.time() - params.analysis.time_start) - continue + ('original', master_problem.idx, master_problem.objective_value, master_problem.optimality_gap, master_problem.primal_variables, master_problem.instance.patterns, time.time() - params.analysis.time_start) + # otherwise, don't update the optimal solution and directly move to next node + else: + print(f"Solution of ILP ({master_problem.objective_value}) is not better than the incumbent solution ({params.solution.objective_value}), fathoming current node and moving to the next one.",flush=True) + continue # otherwise, run normally CG else: column_generation(master_problem) - # update the node - params.graph.update_node(str(master_problem.idx), { - 'obj_val': str(master_problem.objective_value), - 'incumbent_solution': params.solution.objective_value, - 'cpu_time': str(round(time.time() - start_time, 3)), - 'total_patterns': str(number_of_patterns_before), - 'added_patterns': str(len(master_problem.instance.patterns) - number_of_patterns_before) - }) + # if the column generation is not optimal, fathom the current node and move to the next one + if (master_problem.model.status != g.GRB.OPTIMAL and master_problem.model.status != g.GRB.TIME_LIMIT) or master_problem.model.SolCount == 0: + print(f"Solution of master problem is infeasible or cut off. Fathoming current node and moving to the next one.",flush=True) + continue - # if the column generation is not optimal, fathom the current node and move to the next one - if (master_problem.model.status != g.GRB.OPTIMAL and master_problem.model.status != g.GRB.TIME_LIMIT) or master_problem.model.SolCount == 0: - print(f"Solution of master problem is not optimal. Fathoming current node and moving to the next one.",flush=True) - continue + # update the node + params.graph.update_node(str(master_problem.idx), { + 'obj_val': str(master_problem.objective_value), + 'incumbent_solution': params.solution.objective_value, + 'cpu_time': str(round(time.time() - start_time, 3)), + 'total_patterns': str(number_of_patterns_before), + 'added_patterns': str(len(master_problem.instance.patterns) - number_of_patterns_before) + }) # >>>>>>> 2. a) YES, move to 3. if master_problem.objective_value < params.solution.objective_value: @@ -689,7 +662,7 @@ def branch_and_price(instance): # >>>>>> 4. b) NO, move to 5 if (master_problem.model.status != g.GRB.OPTIMAL and master_problem.model.status != g.GRB.TIME_LIMIT) or master_problem.model.SolCount == 0: - print("Solution of master heuristic is infeasible. Moving to the branching.",flush=True) + print("Solution of master heuristic is infeasible or cut off. Moving to the branching.",flush=True) # >>>>>> 4. b) YES, if the solution is better than the solution of incumbent solution, update it and move to 5 else: @@ -746,7 +719,7 @@ def branch_and_price(instance): # colour the node where optimal solution was found params.graph.update_node(str(params.solution.mp_index), {'color': 'yellow'}) - print("All nodes fathomed, returning the best solution.",flush=True) + print("Returning the best solution...",flush=True) def main(path, arguments): diff --git a/solve_ilp.py b/solve_ilp.py index 04dd0bb0c66442af70f6053b15195bb2ee6b89d4..f361f9915c83b936608834bef444e6a2a5e61c24 100644 --- a/solve_ilp.py +++ b/solve_ilp.py @@ -36,12 +36,14 @@ def ilp(instance): o = model.addVars(instance.n_blocks, lb=0, ub=g.GRB.INFINITY, vtype=g.GRB.INTEGER, name="o") z = model.addVars(instance.n_blocks, lb=0, ub=g.GRB.INFINITY, vtype=g.GRB.INTEGER, name="z") n = model.addVars(instance.n_surgeons, instance.n_days, vtype=g.GRB.BINARY, name="n") + y = model.addVars(instance.n_patients+1, instance.n_patients+1, instance.n_blocks, vtype=g.GRB.BINARY, name="y") + e = model.addVars(instance.n_blocks, vtype=g.GRB.BINARY, name="e") # add constraints for patient in instance.patients: model.addConstr(g.LinExpr((1, x[patient.id, block.id]) for block in instance.blocks) <= 1, name=f"constraint_1_p{patient.id}") for block in instance.blocks: - model.addConstr(g.LinExpr((patient.duration, x[patient.id, block.id]) for patient in instance.patients) - block.capacity == o[block.id] - z[block.id], name=f"constraint_2_b{block.id}") + model.addConstr(g.LinExpr((patient.duration, x[patient.id, block.id]) for patient in instance.patients) + g.LinExpr((instance.sequence_duration[patient_1.group][patient_2.group], y[patient_1.id, patient_2.id, block.id]) for patient_1 in instance.patients for patient_2 in instance.patients if patient_2.id != patient_1.id) - block.capacity == o[block.id] - z[block.id], name=f"constraint_2_b{block.id}") model.addConstrs((o[block.id] <= block.overtime for block in instance.blocks), name="constraint_3") for room in instance.rooms: for day in instance.days: @@ -52,6 +54,16 @@ def ilp(instance): model.addConstrs((x[patient.id, block.id] == 0 for block in instance.blocks for patient in instance.patients if sum([1 if patient.surgeon == surgeon.id and block.id in surgeon.blocks else 0 for surgeon in instance.surgeons]) == 0), name="constraint_6") model.addConstrs((x[patient.id, block.id] == 0 for patient in instance.patients for block in instance.blocks if patient.release > block.day), name="constraint_7") + for patient in instance.patients: + for block in instance.blocks: + model.addConstr(x[patient.id, block.id] == g.LinExpr((1, y[patient.id, patient_next.id, block.id]) for patient_next in instance.patients if patient_next.id != patient.id) + y[patient.id, instance.n_patients, block.id]) + model.addConstr(x[patient.id, block.id] == g.LinExpr((1, y[patient_prev.id, patient.id, block.id]) for patient_prev in instance.patients if patient_prev.id != patient.id) + y[instance.n_patients, patient.id, block.id]) + model.addConstr(x[patient.id, block.id] <= e[block.id]) + for block in instance.blocks: + model.addConstr(g.LinExpr((1, y[patient.id, instance.n_patients, block.id]) for patient in instance.patients) == e[block.id]) + model.addConstr(g.LinExpr((1, y[instance.n_patients, patient.id, block.id]) for patient in instance.patients) == e[block.id]) + model.addConstr(e[block.id] <= g.LinExpr((1, x[patient.id, block.id]) for patient in instance.patients)) + # set objective obj1 = g.LinExpr((((day - patient.release)) * patient.priority, x[patient.id, block.id]) for patient in instance.patients for day in instance.days for block in instance.blocks if block.day == day) obj5 = g.LinExpr((1, x[patient.id, block.id]) for patient in instance.patients for block in instance.blocks) @@ -103,8 +115,7 @@ def main(path): print('# SOLUTION') print(f"{'#' * 50}" + utilities.Color.END + '\n') for b in instance.blocks: - patient_vec = variables_ilp['x'][:, b.id] - patients = [instance.get_patient(idx) for idx, item in enumerate(patient_vec) if item == 1] + patients = [instance.get_patient(idx) for idx, item in enumerate(variables_ilp['x'][:, b.id]) if item == 1] print(f"| BLOCK: {b.id} | PATIENTS: {[patient.id for patient in patients]} | SURGEONS: {list(set([patient.surgeon for patient in patients]))} | OVERTIME: {variables_ilp['o'][b.id]} | UNDERTIME: {variables_ilp['z'][b.id]}") print(f"OBJECTIVE VALUE: {model_ilp.objVal}") print(f"OPTIMALITY GAP: {model_ilp.MIPGap}") diff --git a/utilities.py b/utilities.py index 424f7861cfdf444c4d8d485cf998f4e30d6e6b80..2cbe1cb875d82a727f75c2869f8589debbdffb54 100644 --- a/utilities.py +++ b/utilities.py @@ -169,19 +169,13 @@ class Solution: if self.type == 'original': for b in instance.blocks: - patient_vec = self.variables['x'][:, b.id] - patients = [instance.get_patient(idx) for idx, item in enumerate(patient_vec) if item == 1] + patients = [instance.get_patient(idx) for idx, item in enumerate(self.variables['x'][:, b.id]) if item == 1] print( f"| BLOCK: {b.id} | PATIENTS: {[patient.id for patient in patients]} | SURGEONS: {list(set([patient.surgeon for patient in patients]))} | OVERTIME: {self.variables['o'][b.id]} | UNDERTIME: {self.variables['z'][b.id]}", flush=True) else: for block in instance.blocks: - patterns = [p for p in self.patterns if p.block.id == block.id] - sumproducts = [sum([(1 if patient.id in pattern.patient_id else 0) * self.variables['theta'][pattern.id] for pattern in patterns]) for patient in instance.patients] - patients = [patient for idx, patient in enumerate(instance.patients) if sumproducts[idx] == 1] - surgeons = list(set([patient.surgeon for patient in patients])) - overtime = max(0, sum([patient.duration for patient in patients]) - block.capacity) - undertime = max(0, block.capacity - sum([patient.duration for patient in patients])) - print(f'| BLOCK: {block.id} | PATIENTS: {[patient.id for patient in patients]} | SURGEONS: {surgeons} | OVERTIME: {overtime} | UNDERTIME: {undertime}', flush=True) + pattern = [p for p in self.patterns if p.block.id == block.id and round(self.variables['theta'][p.id], 2) == 1][0] + pattern.print_pattern() print(f"OBJECTIVE VALUE: {self.objective_value}", flush=True) print(f"OPTIMALITY GAP: {self.optimality_gap}", flush=True)