-
Notifications
You must be signed in to change notification settings - Fork 1
/
model_PAMIP.py
140 lines (115 loc) · 6.5 KB
/
model_PAMIP.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
"""File Name: model_PAMIP.py
MIP Model formulation of PA-MIP
Purpose of this function: Build optimization model and ensure Model-data separation
Input: Weight vectors of decision variables in the objective function, sets and indices of the variables, and related input data
Output: Decision variables of the model and their attributes
Usage: This function will be called inside main.ipynb to invoke the necessary optimization model using data in that Python code.
Developer: Tanmoy Das
Date: Aug 2023
"""
# %% Model
"""### Model Deployment
We now determine the MIP model for the facility location problem, by defining the decision variables, constraints, and objective function. Next, we start the optimization process and Gurobi finds the plan to build facilities that minimizes total costs.
"""
# import libraries
import gurobipy as gp
from gurobipy import GRB
import pandas as pd
import custom_functions as func
def solve(W, coordinates_st, coordinates_spill, pairings, SizeSpill, Sensitivity, TimeR, NumberStMax, m, spill_data):
# Organizing data
w1 = W[0];
w2 = W[1];
w3 = W[2]
num_facilities = len(coordinates_st);
num_customers = len(coordinates_spill) # rename customer
# MIP model formulation
model = gp.Model('model PA-MIP')
# Decision variables
cover = model.addVars(pairings, vtype=GRB.BINARY, name='Cover')
select = model.addVars(num_facilities, vtype=GRB.BINARY, name='Select')
amount = model.addVars(pairings, vtype=GRB.CONTINUOUS, lb=0, name='Amount')
# Objective function
if m == 'm1':
objective_m1 = gp.quicksum(
(SizeSpill[c] + Sensitivity[c] - TimeR[c, f]) * cover[c, f]
for c, f in pairings)
print('m:', m)
model.setObjective(objective_m1, GRB.MAXIMIZE)
elif m == 'm2':
objective_m1 = gp.quicksum(
(w1 * SizeSpill[c] + w2 * Sensitivity[c] - w3 * TimeR[c, f]) * cover[c, f]
for c, f in pairings)
print('m:', m)
model.setObjective(objective_m1, GRB.MAXIMIZE)
# Constraints
# C2
C_cover_by_open_facility = model.addConstrs((cover[c, f] <= select[f] for c, f in pairings), name='C_cover_by_open_facility')
# C3
C_max_facility = model.addConstr((gp.quicksum(select[f] for f in range(num_facilities)) == NumberStMax), name='C_max_facility')
# C4
C_one_facility = model.addConstrs((cover.sum(c, '*') <= 1 for c in range(num_customers)), name='C_one_facility')
# C5: deploy less than demand
C_demand_allocation = model.addConstrs((deploy[o, s, r] <= Demand[o, r]
for o, s, r in osr_pair), name='C_demand_allocation')
# C6: resource capacity constaint & deploy only when facility is open
C_capacity_constraint = model.addConstrs((deploy.sum('*', s, r) <= BigM * Availability[s, r] * select[s] #
for s, r in sr_pair), name='C_capacity_constraint')
# C7
C_Arctic_remoteness = model.addConstr((gp.quicksum(select[f] for f in range(num_facilities_uN)) <= n_u), name='C_Arctic_remoteness')
# Write and run the model
model.write('Outputs/model PA-MIP.lp')
model.optimize()
# %% Analyse model
# Obtain model results & carry them outside of the model scope
mvars = model.getVars() # these values are NOT accessible outside of the model scope
names = model.getAttr('VarName', mvars)
values = model.getAttr('X', mvars)
# def extract_ones_DV(model, cover, select, amount, spill_data):
cover_series = pd.Series(model.getAttr('X', cover))
cover_1s = cover_series[cover_series > 0.5]
select_series = pd.Series(model.getAttr('X', select))
select_1s = select_series[select_series > 0.5]
cover_1s_df = cover_1s.reset_index()
cover_1s_df.rename(columns={'level_0': 'Spill #', 'level_1': 'Station no.'}, inplace=True)
# Demand_df = pd.DataFrame(spill_data[['Resource needed']]) #DS1
Demand_df = pd.DataFrame(SizeSpill)
Demand_df.columns=['Resource needed'] # dont need this for DS1
Demand_df['Spill #'] = Demand_df.index
amountSt = pd.merge(cover_1s_df, Demand_df)
amountSt_groupby = amountSt.groupby(by="Station no.").sum()
# normalize amount in each station
amountSt_groupby['amountSt_display'] = 800 * amountSt_groupby['Resource needed'] / (
amountSt_groupby['Resource needed'].max())
# Calculate Coverage
coverage_percentage = int(100 * len(cover_1s) / len(SizeSpill))
# Extract assignment variables
sol_y = pd.Series(model.getAttr('X', cover))
sol_y.name = 'Assignments'
sol_y.index.names = ['Spill #', 'Station no.']
assignment = sol_y[sol_y > 0.5].to_frame()
assignment_name = assignment.reset_index()
# organize data
spill_df = pd.DataFrame(coordinates_spill, columns=['Spill_Latitude', 'Spill_Longitude']).reset_index() # only 3 columns
# spill_df['Resource needed'] = spill_data['Resource needed'] # DS1
spill_df['Resource needed'] = pd.DataFrame(SizeSpill)
station_df = pd.DataFrame(coordinates_st, columns=['St_Latitude', 'St_Longitude']).reset_index() # only 3 columns
spill_df.rename(columns={'index': 'Spill #'}, inplace=True)
station_df.rename(columns={'index': 'Station no.'}, inplace=True)
assignment2 = pd.merge(assignment_name.reset_index()[['Spill #', 'Station no.']],
station_df[['Station no.', 'St_Latitude', 'St_Longitude']])
#df_temp = spill_df[['Spill #', 'St_Latitude', 'St_Longitude']]
#df_temp.rename(columns={'Extracted_1': 'Spill_Latitude', 'Extracted_2': 'Spill_Longitude'}, inplace=True)
assignment3 = pd.merge(assignment2, spill_df[['Spill #', 'Spill_Latitude', 'Spill_Longitude']])
# Calculate total distance travelled
DistanceTravelled = []
for i in range(len(assignment3)):
st_coord = (assignment3[['St_Latitude', 'St_Longitude']]).iloc[i, :].values
sp_coord = (assignment3[['Spill_Latitude', 'Spill_Longitude']]).iloc[i, :].values
aaa = DistanceTravelled.append(func.compute_distance(st_coord, sp_coord))
DistanceTravelled = sum(DistanceTravelled) # Section 4.1 response time equation ; paragraph inline
ResponseTimeT = int(DistanceTravelled / 60)
return model, cover, select, amount, mvars, names, values, \
cover_1s, select_1s, amountSt_groupby, coverage_percentage, \
ResponseTimeT, assignment3, spill_df, station_df, \
sol_y, assignment, assignment2, assignment_name