[Python] Docplex and callback problems

I’m trying to implement a custom callback that solves a subproblem during the optimization process. I want to:

1. 
   Get the values of certain variables from the main model's solution.
   
   
2. 
   Solve a subproblem (a separate optimization problem) based on these values.
   
   
3. 
   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] * d * 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 * 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 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] * d * 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()

Continue reading...