PuLP Model with Gurobi solver is very slow
Hi!
I have a model with about 21k rows, 18k columns and 179k nonzeros (I have seen that if I expand my model to 136k rows, 125k columns and 1M nonzeros it may run for many hours). When I use prob.solve(GUROBI()) it can take almost an hour be I get the optimal solution. Can anyone tell me why or what I can do with my model to improve the time solving my problem.
My Python file:
from pulp import *
import gurobipy
import numpy as np
import pandas as pd
import itertools
from itertools import combinations
###SETS###
OPERATION = ['depot', 'Ballasting', 'Mating', 'Towing', 'Anchoring', 'Chains', 'Hook-Up', 'Cable-Laying'] #indexed by i
VESSELTYPES = ['TUG','AHTS','HLV','CLV'] #indexed by v
VESSELNUMBER = np.arange(0, 5, 1) #indexed by k
TURBINENUMBER = np.arange(0, 4, 1)
#print(TURBINENUMBER)
#indexed by w
FUNCTIONS = ['BollardPull', 'LiftingCapacity', 'AnchorHandling', 'CableLaying'] #indexed by f
OPERATIONS = np.arange(0, len(OPERATION)+(len(OPERATION)-1)*(len(TURBINENUMBER)-1), 1)
#print(len(OPERATIONS))
#Oil price $25 to euros
oil_price = 22.77
#Penalty cost in euros/hour
penaltycost = 5000
#BIG-M
Big_M = 100000
###DICTIONARIES###
#Dictionary of duration for each operation in hours
dur = np.array([0, 42, 24, 30, 72, 72, 24, 96])
#Scaling the array to yield for all turbines
duration = np.concatenate((dur, np.tile(dur[1:8], len(TURBINENUMBER) - 1)))
#print(duration)
#3 suction anchors per turbine and 1 day per suction anchor
#3 chains per turbine and 1 day per chains
#1 day on each end and 2 days for the static part
#Dictionary of fixed costs for each vesseltype in 1000 euro
fixed_cost = {0 : 10,
1 : 40,
2 : 1000,
3 : 80}
#Dictionary of variable costs for each vesseltype in 1000 euro
variable_cost = {0 : np.round(30/24),
1 : np.round(150/24),
2 : np.round(500/24),
3 : np.round(300/24)}
#Fuelcost for each vesseltype in 1000 euro per hour. (from standard cubic meter/barrel of oil = 0.1590)
fuel_cost = {0 : np.round((15/24*0.1590*oil_price)/1000),
1 : np.round((34/24*0.1590*oil_price)/1000),
2 : np.round((50/24*0.1590*oil_price)/1000),
3 : np.round((30/24*0.1590*oil_price)/1000)}
#Dictionary of speed for each vesseltype in knots
speed = {0 : 15,
1 : 12,
2 : 9.5,
3 : 12}
#####PARAMETERS#####
#MATRICES
#Distance Matrix in nautical miles TUG = 15kn, AHTS = 12kn and HLV = 9.5kn
Matrix_Dist = np.array([[0, 0, 0, 0, 0, 86, 86, 86],
[0, 0, 0, 0, 0, 86, 86, 86],#Anchoring starts with loading of anchors in port
[0, 0, 0, 0, 0, 86, 86, 86],
[86, 86, 86, 86, 86, 0, 0, 0],
[86, 86, 86, 86, 86, 0, 0, 0],
[86, 86, 86, 86, 86, 0, 0, 0],
[86, 86, 86, 86, 86, 0, 0, 0],
[86, 86, 86, 86, 86, 0, 0, 0]])
#Scaling the matrix to yield for all turbines
Matrix_D = np.concatenate((Matrix_Dist, np.tile(Matrix_Dist[:,1:8], len(TURBINENUMBER) - 1)), axis=1)
Matrix_Distance = np.concatenate((Matrix_D, np.tile(Matrix_D[1:8,:], [len(TURBINENUMBER) - 1,1])), axis=0)
#Function matrix; Bollard Pull and Lifting capacity in tonnes, anchor handling equipment and cable laying equipment for
# each vesseltype
Matrix_Func = np.array([[100, 0, 0, 0],
[350, 0, 1, 0],
[0, 14000, 0, 0],
[0, 0, 0, 1]])
#Requirement matrix; Bollard Pull, Lifting capacity in tonnes, anchor handling equipment and cable laying equipment
Matrix_Requirements = np.array([[0, 0, 0, 0],
[0, 0, 0, 0],
[0, 500, 0, 0],
[300, 0, 0, 0],
[300, 0, 1, 0],
[300, 0, 1, 0],
[0, 0, 1, 0],
[0, 0, 0, 1]])
#Scaling the matrix to yield for all turbines
Matrix_Req = np.concatenate((Matrix_Requirements, np.tile(Matrix_Requirements[1:8,:], [len(TURBINENUMBER) - 1,1])), axis=0)
#Restrictions matrix; Significant waveheight [m] and wind speed [m/s] for each operation
Matrix_Restrictions = np.array([[np.inf, np.inf],
[np.inf, np.inf],
[0.5, 12],
[3, 20],
[2, np.inf],
[3.5, np.inf],
[3.5, np.inf],
[2.5, np.inf]])
#Scaling the matrix to yield for all turbines
Matrix_Restr = np.concatenate((Matrix_Restrictions, np.tile(Matrix_Restrictions[1:8,:], [len(TURBINENUMBER) - 1,1])), axis=0)
#Matrix showing which operations that are depending on each other
Matrix_Dependency = np.array([[0, 1, 1, 1, 1, 1, 1, 1],
[0, 0, 1, 1, 0, 0, 1, 0],
[0, 0, 0, 1, 0, 0, 1, 0],
[0, 0, 0, 0, 0, 0, 1, 0],
[0, 0, 0, 0, 0, 1, 1, 0],
[0, 0, 0, 0, 0, 0, 1, 0],
[0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 1, 0]])
#Scaling the matrix to yield for all turbines
Matrix_Dep = np.zeros([len(Matrix_Dependency)+(len(Matrix_Dependency)-1)*(len(TURBINENUMBER)-1),
len(Matrix_Dependency)+(len(Matrix_Dependency)-1)*(len(TURBINENUMBER)-1)])
Matrix_Dep[0:len(Matrix_Dependency),0:len(Matrix_Dependency)] = Matrix_Dependency
for i in range(1, len(TURBINENUMBER), 1):
Matrix_Dep[0:1, len(Matrix_Dependency):len(Matrix_Dependency)+i*(len(Matrix_Dependency)-1)] = np.tile(Matrix_Dependency[0:1,1:8],[1,i])
Matrix_Dep[len(Matrix_Dependency)+(i-1)*(len(Matrix_Dependency)-1):len(Matrix_Dependency)+i*(len(Matrix_Dependency)-1),len(Matrix_Dependency)+(i-1)*(len(Matrix_Dependency)-1):len(Matrix_Dependency)+i*(len(Matrix_Dependency)-1)] = Matrix_Dependency[1:8,1:8]
#####DECISION VARIABLES#####
#BINARY VARIABLES
#1 if vessel k is used for operation i, 0 else
vessop_vars = LpVariable.dicts("VesselUsedOperation",
[(i,v,k)
for i in range(len(OPERATIONS))
for v in range(len(VESSELTYPES))
for k in range(len(VESSELNUMBER))],0,1,LpBinary
)
#1 if vessel k of vesseltype v sails from operation i to operation j, 0 else
route_vars = LpVariable.dicts("VesselRoute",
[(i,j,v,k)
for i in range(len(OPERATIONS))
for j in range(len(OPERATIONS))
for v in range(len(VESSELTYPES))
for k in range(len(VESSELNUMBER))],0,1,LpBinary
)
#CONTINOUS VARIABLES
#Time which vessel k of vessel type v starts operation i for turbine number w (all activities preceding node i have been completed
starttime_vars = LpVariable.dicts("TimeOperationOccursForVessel",
[(i,v,k)
for i in range(len(OPERATIONS))
for v in range(len(VESSELTYPES))
for k in range(len(VESSELNUMBER))],0
)
#Time which operation i occurs for turbine number w (all activities preceding node i have been completed
startop_vars = LpVariable.dicts("TimeOperationOccurs",
[i for i in range(len(OPERATIONS))],0
)
#Time which vessel k of vesseltype v ends operation i for turbine number w
end_vars = LpVariable("TimeVesselEnd",0
)
#Total time of the project
totaltime_vars = LpVariable.dicts("TotalTime",
[(v,k)
for v in range(len(VESSELTYPES))
for k in range(len(VESSELNUMBER))],0
)
#SET PROBLEM VARIABLE
prob = LpProblem("FSM", LpMinimize)
#####OBJECTIVE FUNCTION#####
prob += lpSum(fixed_cost[v]*route_vars[(0,j,v,k)]
for j in range(len(OPERATIONS))
for v in range(len(VESSELTYPES))
for k in range(len(VESSELNUMBER))) \
+ lpSum(variable_cost[v]*totaltime_vars[(v,k)]
for v in range(len(VESSELTYPES))
for k in range(len(VESSELNUMBER))) \
+ lpSum(fuel_cost[v]*duration[i]*vessop_vars[(i,v,k)]
for i in range(len(OPERATIONS))
for v in range(len(VESSELTYPES))
for k in range(len(VESSELNUMBER))
) + end_vars*penaltycost
#####CONSTRAINTS#####
#A vessel can leave the depot maximum one time. This will happen if it is used for an operation.
for v in range(len(VESSELTYPES)):
for k in range(len(VESSELNUMBER)):
prob += lpSum(route_vars[(0,j,v,k)]
for j in range(1, len(OPERATIONS), 1)) <= 1
#All operations must be done
for i in range(1, len(OPERATIONS), 1):
prob += lpSum(route_vars[(i,j,v,k)]
for j in range(len(OPERATIONS))
for v in range(len(VESSELTYPES))
for k in range(len(VESSELNUMBER))) >= 1
#An operation can not be done more than once
for i in range(1, len(OPERATIONS), 1):
for v in range(len(VESSELTYPES)):
for k in range(len(VESSELNUMBER)):
prob += lpSum(route_vars[(i,j,v,k)]
for j in range(len(OPERATIONS))) <= 1
#Flow conservation: Every operation that is visited must also be left at some point
for i in range(len(OPERATIONS)):
for v in range(len(VESSELTYPES)):
for k in range(len(VESSELNUMBER)):
prob += lpSum(route_vars[(i,j,v,k)]
for j in range(len(OPERATIONS))) \
== lpSum(route_vars[(j,i,v,k)]
for j in range(len(OPERATIONS)))
#Functional requirement constraints
#The sum of the vessels used must satisfy the functional constraints for every operation
for j in range(len(OPERATIONS)):
for f in range(len(FUNCTIONS)):
prob += lpSum(route_vars[(i,j,v,k)]*Matrix_Func[v][f]
for i in range(len(OPERATIONS))
for v in range(len(VESSELTYPES))
for k in range(len(VESSELNUMBER))) >= Matrix_Req[j][f]
#Towing requires at least three vessels to keep the turbine stable
for j in range(3,len(OPERATIONS),7):
prob += lpSum(route_vars[(i,j,v,k)]
for i in range(len(OPERATIONS))
for v in range(len(VESSELTYPES))
for k in range(len(VESSELNUMBER))) >= 3
#Towing requires at least one AHTS
for j in range(3, len(OPERATIONS), 7):
prob += lpSum(route_vars[(i,j,1,k)]
for i in range(len(OPERATIONS))
for k in range(len(VESSELNUMBER))) >= 1
#Vessels used
#If the sum of routes to operation j for vessel k of vesseltype v and tubine number w are more than 0, it means that operation j is done with vessel k of vesseltype v for turbine number w
for j in range(len(OPERATIONS)):
for v in range(len(VESSELTYPES)):
for k in range(len(VESSELNUMBER)):
prob += lpSum(route_vars[(i,j,v,k)]
for i in range(len(OPERATIONS))) == vessop_vars[(j,v,k)]
#Time-constraints
#Starttime of depot must be before any other operations
for j in range(1, len(OPERATIONS),1):
for v in range(len(VESSELTYPES)):
for k in range(len(VESSELNUMBER)):
prob += starttime_vars[(0,v,k)] + np.round(Matrix_Distance[0][j]/speed[v])*route_vars[(0,j,v,k)] \
<= starttime_vars[(j,v,k)]
##If a vessel travels from operation 𝑖 and 𝑗 operation j must start after operation i.
for i in range(1, len(OPERATIONS), 1):
for j in range(1, len(OPERATIONS), 1):
for v in range(len(VESSELTYPES)):
for k in range(len(VESSELNUMBER)):
prob += starttime_vars[(i,v,k)] + duration[i] + np.round(Matrix_Distance[i][j]/speed[v]) \
<= starttime_vars[(j,v,k)] + (1-route_vars[(i,j,v,k)])*Big_M
#Time operation starts
for i in range(1, len(OPERATIONS), 1):
for v in range(len(VESSELTYPES)):
for k in range(len(VESSELNUMBER)):
prob += starttime_vars[(i,v,k)] - (1-vessop_vars[(i,v,k)])*Big_M <= startop_vars[i]
#If an operation j is dependent on operation i being done first, start time of operation i must be before operation j
for i in range(1, len(OPERATIONS), 1):
for j in range(1, len(OPERATIONS), 1):
prob += (startop_vars[i] + duration[i])*Matrix_Dep[i][j] <= startop_vars[j]
#The starttime for an operation must be the same for all vessels operating that operation
for i in range(1, len(OPERATIONS), 1):
for v in range(len(VESSELTYPES)):
for k in range(len(VESSELNUMBER)):
prob += starttime_vars[(i,v,k)] - (1-vessop_vars[(i,v,k)])*Big_M <= startop_vars[i]
for i in range(1, len(OPERATIONS), 1):
for v in range(len(VESSELTYPES)):
for k in range(len(VESSELNUMBER)):
prob += starttime_vars[(i,v,k)] + (1-vessop_vars[(i,v,k)])*Big_M >= startop_vars[i]
#Totaltime
for v in range(len(VESSELTYPES)):
for k in range(len(VESSELNUMBER)):
prob += lpSum((np.round(Matrix_Distance[i][j]/speed[v]) + duration[j]) * route_vars[(i,j,v,k)]
for i in range(len(OPERATIONS))
for j in range(len(OPERATIONS))) <= totaltime_vars[(v,k)]
for i in range(1, len(OPERATIONS), 1):
for v in range(len(VESSELTYPES)):
for k in range(len(VESSELNUMBER)):
prob += starttime_vars[(i,v,k)] + (duration[i] + np.round(Matrix_Distance[i][0]/speed[v]))*route_vars[(i,0,v,k)] - starttime_vars[(0,v,k)] <= totaltime_vars[(v,k)]
#Endtime of project
for v in range(len(VESSELTYPES)):
for k in range(len(VESSELNUMBER)):
prob += totaltime_vars[(v, k)] <= end_vars
#SOLUTION
prob.writeMPS("LPModel.mps")
prob.solve(GUROBI())
print("Status:",LpStatus[prob.status])
#PRINT DECISION VARIABLES
#o = [{'Vessel':"{}{}".format(v,k),
#'quant':vess_vars[(v,k)].varValue}
#for v in range(len(VESSELTYPES)) for k in range(len(VESSELNUMBER))]
#print(pd.DataFrame(o))
#df = pd.DataFrame(o)
#df.to_csv(r'Data.csv', sep='\t', index = False, header=True)
#####PRINT OPTIMAL SOLUTION#####
print("The cost of vessel fleet size and mix in 1000 euros=",value(prob.objective))
print()
for i in range(len(OPERATIONS)):
for v in range(len(VESSELTYPES)):
for k in range(len(VESSELNUMBER)):
if vessop_vars[(i,v,k)].varValue >= 1:
print(vessop_vars[(i,v,k)])
print(vessop_vars[(i,v,k)].varValue)
print()
for i in range(len(OPERATIONS)):
for v in range(len(VESSELTYPES)):
for k in range(len(VESSELNUMBER)):
if vessop_vars[(i,v,k)].varValue >= 1:
print(starttime_vars[(i,v,k)])
print(starttime_vars[(i,v,k)].varValue)
print()
for v in range(len(VESSELTYPES)):
for k in range(len(VESSELNUMBER)):
print(totaltime_vars[(v, k)])
print(totaltime_vars[(v, k)].varValue)
print()
for i in range(len(OPERATIONS)):
for j in range(len(OPERATIONS)):
for v in range(len(VESSELTYPES)):
for k in range(len(VESSELNUMBER)):
if route_vars[(i,j,v,k)].varValue >= 1:
print(route_vars[(i,j,v,k)])
print(route_vars[(i,j,v,k)].varValue)
print(starttime_vars[(13,2,1)].varValue)
print(totaltime_vars[(2,1)].varValue)
########################################
Then what the solver gives me:
Academic license - for non-commercial use only
Gurobi Optimizer version 9.0.2 build v9.0.2rc0 (mac64)
Optimize a model with 21196 rows, 18029 columns and 178929 nonzeros
Model fingerprint: 0x3c012e24
Variable types: 629 continuous, 17400 integer (0 binary)
Coefficient statistics:
Matrix range [1e+00, 1e+05]
Objective range [1e+00, 5e+03]
Bounds range [1e+00, 1e+00]
RHS range [1e+00, 1e+05]
Presolve removed 2524 rows and 600 columns
Presolve time: 0.43s
Presolved: 18672 rows, 17429 columns, 98020 nonzeros
Variable types: 629 continuous, 16800 integer (16800 binary)
Root relaxation: objective 8.601680e+05, 3361 iterations, 0.25 seconds
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
0 0 860168.000 0 108 - 860168.000 - - 2s
0 0 860168.000 0 166 - 860168.000 - - 5s
0 0 860235.746 0 175 - 860235.746 - - 65s
0 0 860235.746 0 175 - 860235.746 - - 65s
0 0 860332.962 0 97 - 860332.962 - - 68s
0 0 860332.962 0 170 - 860332.962 - - 72s
0 0 860332.962 0 173 - 860332.962 - - 73s
0 0 860332.962 0 169 - 860332.962 - - 74s
0 0 860332.962 0 133 - 860332.962 - - 78s
0 0 860332.962 0 182 - 860332.962 - - 79s
0 0 860332.962 0 134 - 860332.962 - - 82s
0 0 860332.962 0 175 - 860332.962 - - 83s
0 0 860332.962 0 128 - 860332.962 - - 87s
0 0 860332.962 0 180 - 860332.962 - - 88s
0 0 860332.962 0 133 - 860332.962 - - 90s
H 0 0 2051884.0000 860332.962 58.1% - 90s
0 0 860332.962 0 133 2051884.00 860332.962 58.1% - 92s
0 2 860342.962 0 109 2051884.00 860342.962 58.1% - 114s
3 8 1592168.00 2 129 2051884.00 860342.962 58.1% 2137 115s
27 27 1592383.67 8 115 2051884.00 860342.962 58.1% 478 137s
158 157 1592744.09 34 59 2051884.00 860342.962 58.1% 156 140s
459 416 1596790.00 98 40 2051884.00 860342.962 58.1% 103 145s
791 677 887649.714 21 99 2051884.00 860342.962 58.1% 75.1 150s
1090 835 923448.102 31 133 2051884.00 860342.962 58.1% 77.9 156s
1092 836 1592882.00 45 110 2051884.00 860342.962 58.1% 77.8 163s
1093 837 888426.000 40 93 2051884.00 860416.000 58.1% 77.7 221s
1095 838 917953.739 22 126 2051884.00 860556.512 58.1% 77.6 227s
1101 842 1159450.02 97 141 2051884.00 860673.788 58.1% 77.1 234s
1103 844 1756015.01 150 160 2051884.00 860767.776 58.0% 77.0 235s
1110 848 918776.320 46 129 2051884.00 860864.529 58.0% 76.5 241s
1115 852 1592704.00 19 127 2051884.00 860910.562 58.0% 76.2 252s
1117 853 1644034.00 58 135 2051884.00 860910.562 58.0% 76.0 259s
1118 854 1596790.00 93 158 2051884.00 860910.563 58.0% 76.0 261s
1122 858 1470038.00 144 133 2051884.00 860910.563 58.0% 175 265s
-
Official comment
-
Gurobi Staff
This post is more than three years old. Some information may not be up to date. For current information, please check the Gurobi Documentation or Knowledge Base. If you need more help, please create a new post in the community forum. Or why not try our AI Gurobot?. -
-
It seems like it is all the "Time Constraints" (or especially the second one) that makes the solver more time consuming. Because when I comment them out the solver finds a solution after 30 sec.
0 -
-
Can it be because I have (4x7)! = 3*10ˆ(29) routes and the solver has to check this for every vessel?
0 -
-
You wrote out the MPS file. How about you try running gurobi_cl.exe on the MPS file and see if you have similar performance?
0 -
Post is closed for comments.
Comments
4 comments