I’m trying to implement a custom callback that solves a subproblem during the optimization process. I want to:
Get the values of certain variables from the main model's solution.
Solve a subproblem (a separate optimization problem) based on these values.
If the subproblem solution suggests a cut, add it to the main model as a new constraint.
I'll share my code here, i really wanted to do it in python but im so lost right now, ive been trying also xpress and similar but documentation is useless, any help would be greatly appreciated.
from docplex.mp.model import Model
from cplex.callbacks import UserCutCallback
from docplex.mp.callbacks.cb_mixin import ConstraintCallbackMixin
import random
class CustomCutCallback(ConstraintCallbackMixin, UserCutCallback):
def __init__(self, env):
UserCutCallback.__init__(self, env)
ConstraintCallbackMixin.__init__(self)
self.eps = 1e-6
self.nb_cuts = 0
self.cts = []
def add_cut_constraint(self, cut):
self.register_constraint(cut)
def __call__(self, context):
"""
This method is invoked by CPLEX during optimization.
It receives the 'context' to interact with the model and variables.
"""
print("Callback called with context")
m = context.model # The main model from the context
m_sol = m.solution
x_values = {var: m_sol.get_value(var) for var in m.variables if var.name.startswith('x')}
w_values = {var: m_sol.get_value(var) for var in m.variables if var.name.startswith('w')}
m2_cuts = self.solve_subproblem(x_values, w_values, m, context)
def solve_subproblem(self, x_values, w_values, m, context):
"""
Solves the subproblem and generates cuts if necessary.
"""
m2 = Model(name='subproblem')
client_range = range(len(x_values))
deposit_range = range(len(w_values))
plant_range = range(len(w_values[0]))
alpha = m2.continuous_var_matrix(client_range, deposit_range, name='alpha', lb=0)
beta = m2.continuous_var_matrix(deposit_range, plant_range, name='beta', lb=0)
m2.add_constraints(alpha[i, j] + (x_values.get((i, j), 0) * beta[j, k]) <= x_values.get((i, j), 0) * w_values.get((j, k), 0)
for i in client_range for j in deposit_range for k in plant_range)
m2.maximize(m2.sum(alpha[i, j] * x_values.get((i, j), 0) for i in client_range for j in deposit_range) +
m2.sum(beta[j, k] * w_values.get((j, k), 0) for j in deposit_range for k in plant_range))
m2.solve()
print(m2.solution)
# Here, you perform an evaluation of the cut values
for i in client_range:
for j in deposit_range:
for k in plant_range:
if sum(m2.solution.get_value(alpha[i, j]) * x_values.get((i, j), 0) for i in client_range for j in deposit_range) + \
sum(m2.solution.get_value(beta[j, k]) * w_values.get((j, k), 0) for j in deposit_range for k in plant_range) > \
sum(transport_cost_deposit_client[j][i] * d[i] * x[i, j] for i in client_range for j in deposit_range) + \
sum(transport_cost_plant_deposit[j][k] * w[j, k] for j in deposit_range for k in plant_range):
m.add_constraint(sum(m2.solution.get_value(alpha[i, j]) * x[i, j] for i in client_range for j in deposit_range) + \
sum(m2.solution.get_value(beta[j, k]) * w[j, k] for j in deposit_range for k in plant_range))
# Main model
def build_location_model(transport_cost_deposit_client, transport_cost_plant_deposit, p, q, d, **kwargs):
m = Model(name='location', **kwargs)
num_deposits = len(transport_cost_deposit_client)
num_plants = len(transport_cost_plant_deposit[0])
num_clients = len(d)
deposit_range = range(num_deposits)
plant_range = range(num_plants)
client_range = range(num_clients)
x = m.binary_var_matrix(client_range, deposit_range, name='x')
w = m.integer_var_matrix(deposit_range, plant_range, name='w', lb=0)
y = m.binary_var_list(deposit_range, name='y')
h = m.binary_var_list(plant_range, name='h')
m.add_constraints(m.sum(x[i, j] for j in deposit_range) == 1 for i in client_range)
m.add_constraints(m.sum(w[j, k] for k in plant_range) == y[j] for j in deposit_range)
m.add_constraints(x[i, j] <= y[j] for i in client_range for j in deposit_range)
m.add_constraint(m.sum(h[k] for k in plant_range) == q)
m.add_constraint(m.sum(y[j] for j in deposit_range) == p)
m.add_constraints(m.sum(d[i] * x[i,j] for i in client_range) == m.sum(w[j, k] for k in plant_range) for j in deposit_range)
m.add_constraints(w[j, k] <= (m.sum(d[i] for i in client_range) * h[k]) for j in deposit_range for k in plant_range)
transport_cost = m.sum(transport_cost_deposit_client[j][i] * d[i] * x[i, j] for i in client_range for j in deposit_range) + \
m.sum(transport_cost_plant_deposit[j][k] * w[j, k] for j in deposit_range for k in plant_range)
m.minimize(transport_cost)
m.parameters.preprocessing.presolve = 0
# Register the callback
cut_cb = m.register_callback(CustomCutCallback)
# Configure CPLEX parameters for cuts
params = m.parameters
params.mip.cuts.mircut = -1
m.solve()
return m
# Test function
def solve_model():
num_deposits = 10
num_plants = 4
num_clients = 20
TRANSPORT_COST_DEPOSITS_CLIENTS = [
[random.randint(20, 100) for _ in range(num_clients)] for _ in range(num_deposits)
]
TRANSPORT_COST_PLANTS_DEPOSITS = [
[random.randint(30, 80) for _ in range(num_plants)] for _ in range(num_deposits)
]
p = 5
q = 3
d = [random.randint(5, 20) for _ in range(num_clients)]
m = build_location_model(TRANSPORT_COST_DEPOSITS_CLIENTS, TRANSPORT_COST_PLANTS_DEPOSITS, p, q, d)
if m.solution is None:
print("No valid solution found.")
return None
print(m.solution)
return m
if __name__ == "__main__":
solve_model()
Your basic model ("location") is infeasible. I ran your code with output on solve(log_output=True)
and it stops at once as infeasible. Hence the callback is never called. I suggest you modify the code to start from a feasible model, and then add the callbacks.
Here is the output I of the solve:
C:\python\anaconda202210\envs\docplex_dev38\python.exe C:\OPTIM\PYLAB\stackov\cutcb.py
Model: location
- number of variables: 254
- binary=214, integer=40, continuous=0
- number of constraints: 282
- linear=282
- parameters:
parameters.mip.cuts.mircut = -1
parameters.preprocessing.presolve = 0
- objective: minimize
- problem type is: MILP
Version identifier: 20.1.0.0 | 2020-11-10 | 9bedb6d68
CPXPARAM_Preprocessing_Presolve 0
CPXPARAM_Read_DataCheck 1
CPXPARAM_MIP_Cuts_MIRCut -1
Legacy callback UD
Warning: Control callbacks may disable some MIP features.
Clique table members: 211.
MIP emphasis: balance optimality and feasibility.
MIP search method: traditional branch-and-cut.
Parallel mode: none, using 1 thread.
Root relaxation solution time = 0.00 sec. (0.23 ticks)
Nodes Cuts/
Node Left Objective IInf Best Integer Best Bound ItCnt Gap Variable B NodeID Parent Depth
0 0 infeasible 2
Root node processing (before b&c):
Real time = 0.00 sec. (1.57 ticks)
Sequential b&c:
Real time = 0.00 sec. (0.00 ticks)
------------
Total (root+branch&cut) = 0.00 sec. (1.57 ticks)
No valid solution found.
Process finished with exit code 0