Stock level optimisation
Awaiting user inputI am trying to solve the following problem:
Given the dataframe with columns:
- location
- product
- current_stock
- demand
- price
I want to optimize the unmet demand given additional constraints:
- the number of shipments from one location cannot exceed a given threshold like 50
- the value of each package sent cannot be lower then $250
I use the code below. I think it works quite well if the number of products/locations is not that big like 50 or 100. But when I run it on a full problem with more than 100 locations and couple thousands of products the only solution it finds is 0.
I don't get it, why there is a non zero solution which improves the satisfied demand for 50 products and it cannot improve for a full set. I expected that in worst case it should find the solution for 50 products leaving all other 0.
def optimize_demand_gurobi_with_more_constraints(min_stock, current_stock, demand, price, locations, products,max_time):
log_system_status("Opt started")
m = Model("Unsatisfied_Demand_Optimization")
m.setParam('OutputFlag', 1) # Enable output for more verbose logging
# m.params.Method = 1
# m.params.SoftMemLimit = 170
m.params.Presolve = 2
m.params.PreSparsify = 1
# m.params.NodefileStart = 0.5
# m.params.Threads =
m.params.LogFile = add_timestamp_to_filename('logs/GRB_opt.txt')
m.params.DegenMoves = 0
m.params.TimeLimit = max_time
log_system_status("x")
# Defining variables: Only create shipment variables x where conditions are met
x = m.addVars(
[(loc_a, loc_b, p) for loc_a in tqdm(locations) for loc_b in locations for p in products
if loc_a != loc_b and current_stock[loc_a, p] > demand[loc_a, p] and current_stock[loc_b, p] < demand[loc_b, p]],
vtype=GRB.INTEGER, name="Ship"
)
log_system_status("unmet")
unmet = m.addVars([(loc_b, p) for loc_b in tqdm(locations) for p in products],
vtype=GRB.INTEGER, name="Unmet")
log_system_status("y")
y = m.addVars([(loc_a, loc_b, p) for loc_a in tqdm(locations) for loc_b in locations for p in products if loc_a != loc_b],
vtype=GRB.BINARY, name="NoCrossShip")
# log_system_status("obj")
# Objective: Minimize the total unmet demand
# m.setObjective(quicksum(unmet[loc_b, p] for loc_b in tqdm(locations) for p in products), GRB.MINIMIZE)
log_system_status("Hints")
for k, v in tqdm(x.items()):
x[k].VarHintVal = shipments.get(k,0)
x[k].VarHintPri = 10
for k, v in tqdm(y.items()):
if shipments.get(k,0) == 0:
y[k].VarHintVal = 0
y[k].VarHintPri = 5
else:
y[k].VarHintVal = 1
y[k].VarHintPri = 5
log_system_status("stock & demand")
# Constraints to ensure shipments do not exceed stock available or demand
for loc_a in tqdm(locations):
for p in products:
if any((loc_a, loc_b, p) in x for loc_b in locations if loc_a != loc_b):
m.addConstr(quicksum(x[loc_a, loc_b, p] for loc_b in locations if loc_a != loc_b and (loc_a, loc_b, p) in x) <= current_stock[loc_a, p])
for loc_b in tqdm(locations):
for p in products:
if any((loc_a, loc_b, p) in x for loc_a in locations if loc_a != loc_b):
m.addConstr(quicksum(x[loc_a, loc_b, p] for loc_a in locations if loc_a != loc_b and (loc_a, loc_b, p) in x) <= demand[loc_b, p])
m.addConstr(unmet[loc_b, p] >= price[p]*(demand[loc_b, p]-current_stock[loc_b, p]- quicksum(x[loc_a, loc_b, p] for loc_a in locations if loc_a != loc_b and (loc_a, loc_b, p) in x)))
m.addConstr(unmet[loc_b, p] >=0 )
log_system_status("Crosshipment")
# Constraint to prevent shipment in both directions for the same product
for loc_a in tqdm(locations):
for loc_b in locations:
if loc_a != loc_b:
for p in products:
if (loc_a, loc_b, p) in y and (loc_b, loc_a, p) in y:
m.addConstr(y[loc_a, loc_b, p] + y[loc_b, loc_a, p] <= 1)
if (loc_a, loc_b, p) in x:
m.addConstr(x[loc_a, loc_b, p] <= y[loc_a, loc_b, p] * max(current_stock[loc_a, p], demand[loc_b, p]))
if (loc_b, loc_a, p) in x:
m.addConstr(x[loc_b, loc_a, p] <= y[loc_b, loc_a, p] * max(current_stock[loc_b, p], demand[loc_a, p]))
count_parcels = m.addVars(
locations,
vtype=GRB.INTEGER, name="Ship"
)
total_value = m.addVars([(loc_a, loc_b) for loc_a in tqdm(locations) for loc_b in locations],
vtype=GRB.INTEGER, name="toal_value")
shipment = m.addVars([(loc_a, loc_b) for loc_a in tqdm(locations) for loc_b in locations],
vtype=GRB.BINARY, name="shipment")
for loc_a in tqdm(locations):
for loc_b in locations:
if loc_a != loc_b:
m.addConstr(shipment[loc_a, loc_b] == max_(y[loc_a, loc_b, p] for p in products))
else:
m.addConstr(shipment[loc_a, loc_b] == 0)
for loc_a in tqdm(locations):
m.addConstr(count_parcels[loc_a] == quicksum(shipment[loc_a, loc_b] for loc_b in locations ))
m.addConstr(count_parcels[loc_a] <= 50)
for loc_a in tqdm(locations):
for loc_b in locations:
if loc_a != loc_b:
m.addConstr(total_value[loc_a, loc_b] == quicksum(x.get((loc_a, loc_b, p),0)*price[p] for p in products))
m.addConstr(total_value[loc_a, loc_b] >= 250*shipment[loc_a, loc_b])
log_system_status("Vars and const defined")
log_system_status("obj")
# Objective: Minimize the total unmet demand
# m.setObjective(-unmet.sum(), GRB.MAXIMIZE)
log_system_status("obj")
# Objective: Minimize the total unmet demand
m.setObjective(
quicksum(x[src, dest, prod] for src in locations for dest in locations if src != dest for prod in products if (src, dest, prod) in x),
GRB.MAXIMIZE
)
log_system_status("Optimizing")
m._solution_count = 0
m.optimize()
if m.status in [GRB.OPTIMAL, GRB.TIME_LIMIT, GRB.ITERATION_LIMIT]:
log_system_status("GRB is optimal")
shipment_plan = {k: v.X for k, v in x.items() if v.X > 0}
m.write(add_timestamp_to_filename('model/with_value_constraint.mps'))
m.write(add_timestamp_to_filename('model/model.sol'))
return shipment_plan, m
else:
log_system_status("Optimization was not successful. Status code:"+str(m.status))
if m.status in [GRB.INFEASIBLE, GRB.INF_OR_UNBD]:
m.computeIIS()
m.write("model.ilp")
return {}
-
Hi Krzysztof,
it would be much easier to help you, if you could also include a mathematical formulation of your problem. Also, two examples of input data, for which you get satisfactory and unsatisfactory results would be helpful.
On the very high level, I would guess that the full problem you consider is by a few orders of magnitude larger than the problem for which you found a solution. It may be the case that you will need to turn to some more sophisticated modeling techniques to make your full problem solvable.
Hope this helps.
Best regards
Jonasz0
Please sign in to leave a comment.
Comments
1 comment