CPLEX Python: Current subproblem model in branch and bound

Question

I have an MILP problem and use CPLEX (Python interface). I am working on user heuristics for branching in the branch-and-bound procedure. With HSCallback I managed to get the information about the current subproblem (all integer variables are binary), in particular, what the variables that have been previously fixed/branched on are and how they've been fixed.

Now, I would like to do something with the current subproblem. Is there a way to obtain it as a model? Do I have to construct it each time in the callback from the original model by adding the constraints corresponding to the fixed variables? Or is there a more efficient way of doing it?

Below is a working example for a problem encoding Boolean satisfiability of CNF with 5 variables that has exactly one assignment (corresponding to the commented out line, uncommenting it makes the problem infeasible). I disable all possible heauristics and cuts to make sure that branch and bound has to do some work. Here I construct the subproblem by creating a copy of the original model and adding the constraints corresponding to the current branch.

import datetime
import docplex.mp.model as cpx
import cplex.callbacks as CPX_CB
from docplex.mp.relax_linear import LinearRelaxer


class MyHSCallback(CPX_CB.HSCallback):
    def __call__(self, *args, **kwargs):
        print()
        print(datetime.datetime.now(), "MyHSCallback working.............")
        print("node ID", self.get_node_ID())

        fixed = [(ind, int(self.get_lower_bounds(ind)))
                 for ind in range(self.var_count) if self.get_lower_bounds(ind) == self.get_upper_bounds(ind)]
        print("Fixed variables", ", ".join([str(ind) + ":" + str(value) for ind, value in fixed]))
        print("Feasible subproblem", self.solve())

        # creating a copy of the subproblem
        model_copy = self.original_model.copy()
        for ind, value in fixed:
            model_copy.add_constraint(model_copy.get_var_by_index(ind) == value)

        # doing something with the current subproblem
        lp_model = LinearRelaxer().linear_relaxation(mdl=model_copy)
        sol = lp_model.solve()
        if sol is not None:
            print("Relaxation Values", sol.get_values([lp_model.get_var_by_index(ind) for ind in range(self.var_count)]))


model = cpx.Model(name="cnf")
vars = model.binary_var_list(5, name="var")
x = vars[0]
y = vars[1]
z = vars[2]
w = vars[3]
v = vars[4]

model.add_constraint(w + x + y >= 1)
model.add_constraint(v + w - y >= 0)
model.add_constraint(w - x + y + z >= 0)
model.add_constraint(v - x + y - z >= -1)
model.add_constraint(v - w + y + z >= 0)
model.add_constraint(v - w + x - z >= -1)
model.add_constraint(v - w + x - y + z >= -1)
model.add_constraint(- w - x - y >= -2)
model.add_constraint(-v + w - y + z >= -1)
# model.add_constraint(-v + w + x - y - z >= -2)
model.add_constraint(-v + w - x - z >= -2)
model.add_constraint(-v - w + x >= -1)
model.add_constraint(-v - w - x + y >= -2)

hs_cb = model.register_callback(MyHSCallback)
hs_cb.var_count = len(vars)
hs_cb.original_model = model

model.parameters.preprocessing.symmetry = 0
model.parameters.preprocessing.presolve = 0
model.parameters.preprocessing.boundstrength = 0
model.parameters.preprocessing.coeffreduce = 0
model.parameters.mip.strategy.search = 1 #traditional
model.parameters.mip.strategy.probe = -1
model.parameters.mip.strategy.heuristicfreq = -1
model.parameters.mip.strategy.fpheur = -1
model.parameters.mip.limits.cutpasses = -1
model.parameters.mip.cuts.bqp = -1
model.parameters.mip.cuts.cliques = -1
model.parameters.mip.cuts.covers = -1
model.parameters.mip.cuts.disjunctive = -1
model.parameters.mip.cuts.flowcovers = -1
model.parameters.mip.cuts.gomory = -1
model.parameters.mip.cuts.gubcovers = -1
model.parameters.mip.cuts.implied = -1
model.parameters.mip.cuts.liftproj = -1
model.parameters.mip.cuts.localimplied = -1
model.parameters.mip.cuts.mcfcut = -1
model.parameters.mip.cuts.mircut = -1
model.parameters.mip.cuts.pathcut = -1
model.parameters.mip.cuts.rlt = -1
model.parameters.mip.cuts.zerohalfcut = -1

sol = model.solve()
print()
if sol is None:
    print("Infeasible")
else:
    model.print_solution()
    values = sol.get_values(vars)
    print(values)

```