Fixed-Charge Network Flow Problem
AnsweredDear Sir/Madam,
In the process of solving my model using Benders' decomposition, I applied the callback function to the master problem. However, I am unsure what condition to use for where == ...
. Initially, I used MIPSOL
, but the model stops earlier than expected and does not reach the desired gap between the upper and lower bounds.
What should I do?
-
Hi Saeed,
Using MIPSOL is correct.
When using the callback approach you do not need to monitor lower and upper bounds (unlike the iterative approach). If the callback is implemented properly then you just let the solve run until the gap reaches the threshold defined by your MIP Gap parameter.
- Riley
0 -
Hi Riley,
Thanks for your response, I was dubious about that, since the Master Problem does not solve to optimal point. And therefore the lower bound seems not correct.
Also, you mentioned that the MIP Gap. do you mean that
Master.MIPGAP?
please take a look over the callback portion of my code and instruct me.
def Benders_dual_callback(model,where):
global lower_bound, upper_bound
global optim, feas
if where == GRB.Callback.MIPSOL:
x_vals = model.cbGetSolution(y)
si_val = model.cbGetSolution(si)
obj_val, uu = solve_dual_subproblem(x_values=x_vals)
dual.update()
if obj_val is None:
# Feasibility Cut
feas += 1
FC_1 = 0
print('Feasibility added')
Constraints_u_rhs = {m.getConstrByName(re.sub(r"^[^_]+_", "", uvar.VarName)).ConstrName: (math.ceil(uvar.UnbdRay),m.getConstrByName(re.sub(r"^[^_]+_", "", uvar.VarName)).RHS) for uvar in uu }
for key, values in Constraints_u_rhs.items():
if "BigM" not in key:
FC_1 += values[0]*values[1]
else:
cnstr = re.findall(r'\d+',key)
FC_1 += (m.getCoeff(m.getConstrByName(key),x[int(cnstr[2]),int(cnstr[1]),"F" + cnstr[0]])*y[int(cnstr[2]),int(cnstr[1]),"F" + cnstr[0]]) * values[0]
Master.cbLazy(FC_1 <= 0)
cuts_log.write(f"Feasibility Cut {feas}: {FC_1} <= 0\n")
else:
#Optimality Cut
OC_1 = 0
optim +=1
print('Optimality added')
Constraints_u_rhs = {m.getConstrByName(re.sub(r"^[^_]+_", "", uvar.VarName)).ConstrName: (uvar.X,m.getConstrByName(re.sub(r"^[^_]+_", "", uvar.VarName)).RHS) for uvar in uu }
for key, values in Constraints_u_rhs.items():
if "BigM" not in key:
OC_1 += values[0]*values[1]
else:
cnstr = re.findall(r'\d+',key)
OC_1 += -1*(m.getCoeff(m.getConstrByName(key),x[int(cnstr[2]),int(cnstr[1]),"F" + cnstr[0]])*y[int(cnstr[2]),int(cnstr[1]),"F" + cnstr[0]]) * values[0]
fy = gp.quicksum(y[p,t,f]*Setupcoff[p,t,f] for f in Group_nodes if f[0] == "F" for p in range(1,Produces+1) for t in range(1,T+1))
model.cbLazy(si >= OC_1 )
cuts_log.write(f"Optimality Cut {optim}: {si} >= 0 { OC_1 + fy }")
upper_bound = min(upper_bound,obj_val)
lower_bound = max(lower_bound,model.cbGet(GRB.Callback.MIPSOL_OBJ))
gap = abs(upper_bound-lower_bound)/abs(upper_bound) if upper_bound != 0 else 0
print(f"Gap: {gap:4f} | UB: {upper_bound} | LB: {lower_bound}")
if gap < gap_cutoff:
print('Overal Opt Sol found')
model.terminate()
Master.optimize(Benders_dual_callback)How should I modify this section ?0 -
Hi Saeed,
You should remove this:
upper_bound = min(upper_bound,obj_val) lower_bound = max(lower_bound,model.cbGet(GRB.Callback.MIPSOL_OBJ)) gap = abs(upper_bound-lower_bound)/abs(upper_bound) if upper_bound != 0 else 0 print(f"Gap: {gap:4f} | UB: {upper_bound} | LB: {lower_bound}") if gap < gap_cutoff: print('Overal Opt Sol found') model.terminate()
In the callback approach you do not calculate bounds yourself, or terminate the solve yourself. You only do this for the iterative approach where you loop, and solve master then subproblem, once each per iteration of the loop.
You instead need to let the callback add lazy constraints to the master and the solve will terminate when the MIP gap reaches the threshold set by the MIPGap parameter.
- Riley
0 -
Hi Riley;
Thanks for immediate response.
I am a bit confused. So how should I prove the optimality condition? You meant MIPGAP. so the MIP Gap must be defined on the Master Model?
Bests
Saeed
0 -
Dear Riely
I want to take your attention to the report given by Gurobi: Could you please interpret that for me. (This is the last iteration of the report, I am dubious what is the objective function indeed. Does it reached out to optimality?
Optimality added
Gap: 0.097940 | UB: 3343845.5114316954 | LB: 3016349.6841955036
0 0 3016349.68 0 - - 3016349.68 - - 272s
Set parameter InfUnbdInfo to value 1
Gurobi Optimizer version 12.0.0 build v12.0.0rc1 (win64 - Windows 11.0 (22631.2))
CPU model: 11th Gen Intel(R) Core(TM) i5-1135G7 @ 2.40GHz, instruction set [SSE2|AVX|AVX2|AVX512]
Thread count: 4 physical cores, 8 logical processors, using up to 8 threads
Non-default parameters:
InfUnbdInfo 1
Optimize a model with 8202 rows, 3828 columns and 20056 nonzeros
Model fingerprint: 0xd4424b28
Coefficient statistics:
Matrix range [1e-03, 9e+02]
Objective range [6e+01, 9e+03]
Bounds range [0e+00, 0e+00]
RHS range [2e+01, 3e+03]
Presolve removed 4823 rows and 1054 columns
Presolve time: 0.02s
Presolved: 3379 rows, 2774 columns, 17754 nonzeros
Iteration Objective Primal Inf. Dual Inf. Time
0 handle free variables 0s
1416 4.4859396e+06 0.000000e+00 0.000000e+00 0s
Solved in 1416 iterations and 0.07 seconds (0.06 work units)
Optimal objective 4.485939551e+06
Optimality added
Gap: 0.097464 | UB: 3343845.5114316954 | LB: 3017939.5508621708
* 0 0 0 3017939.5509 3017939.55 0.00% - 286s
Cutting planes:
Lazy constraints: 9
Explored 1 nodes (32 simplex iterations) in 286.60 seconds (0.00 work units)
Thread count was 8 (of 8 available processors)
Solution count 1: 3.01794e+06
Optimal solution found (tolerance 1.00e-04)
Best objective 3.017939550862e+06, best bound 3.017939550862e+06, gap 0.0000%
User-callback calls 121, time in user-callback 286.59 sec0 -
Hi Saeed,
Yes, the MIP gap is referring to the Master problem - I would guess your subproblem is continuous only and the concept of MIP gap for the subproblem doesn't exist.
In your log the gap is reported as 0:
gap 0.0000%
If your subproblem is generating cuts properly then this solution is optimal. You can see from your log you generated 9 benders cuts:
Cutting planes: Lazy constraints: 9
You don't need to prove optimality, in the same way that you don't need to prove optimality for solving a simple MIP model with Gurobi. Your callback just needs to provide feasibility and optimality cuts (using lazy constraints), and Gurobi does the rest.At the end, any solution that exists in your Solution Pool is feasible (if it wasn't you would have cut it off). For a MIP gap of 0 (or sufficiently close to it) you will have a lower bound, which has been shaped by your lazy constraints and Gurobi's cuts, equal to the best solution found. This is how you know it is optimal.
- Riley
0 -
Please Also note that Gurobi without Benders obtained the Incumbent as below:
This bound differs from what Benders found!
9173 7319 3096575.02 57 73 3099230.55 3090773.70 0.27% 144 121s
0 -
Dear Riley
The Process Time for reporting is 286.59 seconds?
Thanks for your reply
V/R
Saeed
0 -
Hi Saeed,
This bound differs from what Benders found!
If my read of your logfile is correct, you're maximizing your objective and it would imply the solution from the Benders approach is better than the one returned by a full model, i.e. non-Benders approach, (and should therefore be infeasible). Sounds like a benders cut is missing or miscalculated. If you save the solution from the benders approach you can fix the values of variables in the full model and solve to verify that it is in fact infeasible. You could also use this solution to help debug your subproblem and cut generation, however I would do this exercise on a small model, not one with thousands of constraints and variables.
The Process Time for reporting is 286.59 seconds?
Yes, this is the total run time for solving. Please don't expect it to be anywhere as fast as running Gurobi on the full model.
- Riley
0 -
Greetings Riley
Accordig to the action being called via Master problem, it can not be regarded as Benders Decomposition. Since the Master problem does not solve for each iteration until optimality. The Benders like cuts would be added while the Master model wolves not iteratively. Subsequently, this is a new way where cuts aded while optimization proceeds. I am still thinking When these cuts are added by gurobi. I am askig about theoretical actions under hood.
Is there an article discussing this call back function in details?
Saeed
0 -
Hi Saeed,
Accordig to the action being called via Master problem, it can not be regarded as Benders Decomposition.
This is not correct. Bender's decomposition simply reformulates the problem into a master and a subproblem (or multiple subproblems) and is independent of the approach used to solve it. Benders decomposition exists even if there are no methods to solve it, in the same way that a linear program still exists even if we did not know how to solve it.
There are two classes of approach, the iterative approach and callback approach. The following blog post should help explain: https://orinanobworld.blogspot.com/2011/10/benders-decomposition-then-and-now.html
Since the Master problem does not solve for each iteration until optimality.
If you're using the callback approach the concept of iterations does not apply. You also can apply the cut generation to any solution in the master, it does not need to be optimal.
The Benders like cuts would be added while the Master model wolves not iteratively.
Correct.
Subsequently, this is a new way where cuts aded while optimization proceeds.
It's new relative to the iterative approach, which was introduced in the 60s, but the approach with callbacks goes back at least 16 years, so not new in absolute terms. See https://pubsonline.informs.org/doi/abs/10.1287/opre.1090.0694 for example.
I am still thinking When these cuts are added by gurobi.
When an integer feasible solution to the master is found by Gurobi, a callback with the MIPSOL callback code is used to extract the solution, pass it to the subproblem(s), which generates a cut (or multiple cuts) which you then must add to Gurobi as a lazy constraint. It is immediately applied.
As an aside, when you solve the LP subproblem(s) you need an extreme point or ray (as described in the blog post). It will be a lot easier if you instruct the subproblems to solve with primal and/or dual simplex, as barrier requires crossover to run to provide a basic solution, and data from InfUnbdInfo will only be available when a basis is available (and infeasibility/unboundedness is detected during barrier you won't have basis information).
- Riley
0 -
Awesome response
Thanks a huge barrel of Vodka!
0
Please sign in to leave a comment.
Comments
12 comments