diff --git a/pyomo/solvers/plugins/solvers/xpress_direct.py b/pyomo/solvers/plugins/solvers/xpress_direct.py
index 8e7975f784a..1af1084ce51 100644
--- a/pyomo/solvers/plugins/solvers/xpress_direct.py
+++ b/pyomo/solvers/plugins/solvers/xpress_direct.py
@@ -32,6 +32,7 @@
 from pyomo.opt.results.solver import TerminationCondition, SolverStatus
 from pyomo.opt.base import SolverFactory
 from pyomo.core.base.suffix import Suffix
+from pyomo.core.base.constraint import Constraint
 import pyomo.core.base.var
 
 
@@ -110,10 +111,476 @@ def __call__(self):
     callback=_finalize_xpress_import,
 )
 
+class CallbackContext(object):
+    """Base class for the contexts that are passed to Xpress callback functions.
+
+    Any callback function supported in Pyomo receives a single argument:
+    an instance of a subclass of `CallbackContext`. This instance carries
+    callback-specific information and also allows callback-specific actions.
+    """
+    def __init__(self, problem, solver, var2idx):
+        super(CallbackContext, self).__init__()
+        self._problem = problem
+        self._solver = solver
+        self._lpsol = None
+        self._mipsol = None
+        self._var2idx = var2idx
+    @property
+    def solver(self):
+        """The Pyomo solver object."""
+        return self._solver
+    @property
+    def problem(self):
+        """The callback local Xpress problem instance."""
+        return self._problem
+    @property
+    def attributes(self):
+        """The attributes of the callback local Xpress problem instance."""
+        return self._problem.attributes
+    @property
+    def controls(self):
+        """The controls of the callback local Xpress problem instance."""
+        return self._problem.controls
+    def _make_value_map(self, values, comp2idx, component_map):
+        """Create a `ComponentMap` that maps the keys in `component_map` to
+        data from `values`.
+        `values`:        Array of values the indices of which correspond to
+                         Xpress variable indices.
+        `comp2idx`:      Maps Xpress variables to their respective index.
+        `component_map`: Maps Pyomo variables (or other stuff) to
+                         Xpress variables (or other stuff)
+        """
+        m = ComponentMap()
+        for p in component_map:
+            m[p] = values[comp2idx[component_map[p]]]
+        return m
+    def getlpsol(self):
+        """Get LP solution for callback.
+
+        The meaning of the returned vector depends on the actual callback
+        context.
+        Returns a `ComponentMap` that maps Pyomo variables to values.
+        """
+        if self._lpsol is None:
+            x = [0.0] * self._problem.attributes.originalcols
+            self._problem.getlpsol(x)
+            self._lpsol = self._make_value_map(x, self._var2idx,
+                                               self._solver._pyomo_var_to_solver_var_map)
+        return self._lpsol
+
+    def getmipsol(self):
+        """Get MIP solution for callback.
+
+        The meaning of the returned vector depends on the actual callback
+        context.
+        Returns a `ComponentMap` that maps Pyomo variables to values.
+        """
+        if self._mipsol is None:
+            x = [0.0] * self._problem.attributes.originalcols
+            self._problem.getmipsol(x)
+            self._mipsol = self._make_value_map(x, self._var2idx,
+                                                self._solver._pyomo_var_to_solver_var_map)
+        return self._mipsol
+
+    def addmipsol(self, sol):
+        """Inject a feasible solution.
+
+        - `sol`: A map from Pyomo variables to solution values.
+                 May be incomplete. If multiple values are specified for the
+                 same variable then it is unspecified which of them will be
+                 used.
+        """
+        var2idx = self._var2idx
+        component_map = self._solver._pyomo_var_to_solver_var_map
+        data = dict()
+        for x in sol:
+            data[var2idx[component_map[x]]] = sol[x]
+        self._problem.addmipsol([data[x] for x in sorted(data)], sorted(data))
+
+    def addcut(self, cut, cuttype=0):
+        """Add a (local) cut to the current node.
+
+        - `cut`:     The cut to add. Must be a linear non-range constraint.
+        - `cuttype`: An integer chosen by the user. The Xpress optimizer does
+                     not use the actual value.
+        """
+
+        if not isinstance(cut, Constraint):
+            cut = Constraint(expr=cut)
+            cut.construct()
+
+        if cut.has_lb():
+            if cut.has_ub():
+                if cut.lower != cut.upper:
+                    raise ValueError("Cannot add ranged cut {0} ".format(cut))
+                sense = 'E'
+                orig_rhs = cut.lower
+            else:
+                sense = 'G'
+                orig_rhs = cut.lower
+        elif cut.has_ub():
+            sense = 'L'
+            orig_rhs = cut.upper
+        else:
+            raise ValueError("Cannot add free cut {0} ".format(cut))
+                
+        if cut._linear_canonical_form:
+            lhs = cut.canonical_form()
+        else:
+            lhs = generate_standard_repn(cut.body, quadratic=False)
+
+        degree = lhs.polynomial_degree()
+        if (degree is None) or (degree > 1):
+            raise ValueError("Cannot add non-linear cut {0} ".format(cut))
+
+        
+        orig_row_coeff = lhs.linear_coefs
+        orig_col_ind = [self._var2idx[self._solver._pyomo_var_to_solver_var_map[x]] for x in lhs.linear_vars]
+        max_coefs = self._problem.attributes.cols
+        col_ind = []
+        row_coef = []
+        red_rhs, status = self._problem.presolverow(sense, orig_col_ind,
+                                                    orig_row_coeff, orig_rhs,
+                                                    max_coefs, col_ind, row_coef)
+        if status != 0:
+            raise RuntimeError('Cut cannot be presolved: %d' % status)
+        self._problem.addcuts([cuttype], [sense], [red_rhs], [0, len(col_ind)],
+                              col_ind, row_coef)
+
+class MessageCallbackContext(CallbackContext):
+    """Data passed to message callbacks.
+
+    In addition to the super class, this class has the following properties:
+    - `msg`     [read]: The message that is sent (may be `None`).
+    - `msgtype` [read]: The type of message sent (info, warning, error).
+    """
+    def __init__(self, problem, solver, var2idx, msg, msgtype):
+        super(MessageCallbackContext, self).__init__(problem, solver, var2idx)
+        self._msg = msg
+        self._msgtype = msgtype
+    @property
+    def msg(self):
+        """Message sent by this callback."""
+        return self._msg
+    @property
+    def msgtype(self):
+        """Type of message sent."""
+        return self._msgtype
+    @property
+    def is_info(self):
+        """`True` if this is an info message."""
+        return self._msgtype == 1
+    @property
+    def is_warning(self):
+        """`True` if this is a warning message."""
+        return self._msgtype == 3
+    @property
+    def is_error(self):
+        """`True` if this is an error message."""
+        return self._msgtype == 4
+    @property
+    def is_flush(self):
+        """`True` if this message is a request to flush output."""
+        return self._msgtype < 0
+
+class OptNodeCallbackContext(CallbackContext):
+    """Data passed to optnode callbacks.
+
+    In addition to the super class, this class has one property:
+    - `infeas` [write]: Whether the node is considered infeasible.
+    """
+    def __init__(self, problem, solver, var2idx):
+        super(OptNodeCallbackContext, self).__init__(problem, solver, var2idx)
+        self._infeas = False
+    @property
+    def infeas(self):
+        """Whether the node is considered infeasible."""
+        return self._infeas
+    @infeas.setter
+    def infeas(self, value):
+        self._infeas = value
+
+class PreIntSolCallbackContext(CallbackContext):
+    """Data passed to preintsol callbacks.
+
+    In addition to the super class, this class has the following properties:
+    - `soltype`       [read]: Solution type (heuristic, optimal node, ...)
+    - `cutoff`  [read/write]: New cutoff value in case solution is not rejected
+    - `reject`       [write]: Whether solution should be rejected
+    - `candidate_sol` [read]: The candidate solution.
+    - `candidate_obj` [read]: The objective value for `candidate_sol`
+    """
+    def __init__(self, problem, solver, var2idx, soltype, cutoff):
+        super(PreIntSolCallbackContext, self).__init__(problem, solver, var2idx)
+        self._soltype = soltype
+        self._cutoff = cutoff
+        self._reject = False
+    @property
+    def soltype(self):
+        """Solution type."""
+        return self._soltype
+    @property
+    def cutoff(self):
+        """Cutoff if solution if accepted."""
+        return self._cutoff
+    @cutoff.setter
+    def cutoff(self, value):
+        """Cutoff if solution if accepted."""
+        self._cutoff = value
+    @property
+    def reject(self):
+        """Should solution be rejected?"""
+        return self._reject
+    @reject.setter
+    def reject(self, value):
+        """Should solution be rejected?"""
+        self._reject = value
+    @property
+    def candidate_sol(self):
+        """Returns the candidate solution as a map from Pyomo variables to
+        values.
+        """
+        return self.getlpsol()
+    @property
+    def candidate_obj(self):
+        """Get the objective value for `candidate_sol`."""
+        return self.attributes.lpobjval
+
+class IntSolCallbackContext(CallbackContext):
+    """Data passed to intsol callbacks.
+
+    In addition to the super class, this class has the following properties:
+    - `solution` [read]: The solution being reported.
+    - `objval`   [read]: The objective value of `solution`.
+    """
+    def __init__(self, problem, solver, var2idx):
+        super(IntSolCallbackContext, self).__init__(problem, solver, var2idx)
+    @property
+    def solution(self):
+        """Returns the solution as a map from Pyomo variables to values."""
+        return self.getmipsol()
+    @property
+    def objval(self):
+        """Get the objective value for `solution`."""
+        return self.attributes.mipobjval
+
+class ChangeBranchObjectCallbackContext(CallbackContext):
+    """Data passed to chgbranchobject callbacks.
+
+    In addition to the super class, this class has the following properties:
+    - `orig_branchobject`  [read]: The branching that Xpress suggests
+    - `branchobject' [read/write]: The branching that will be executed.
+    - `new_object()`:              Create a new branching object.
+    """
+    def __init__(self, problem, solver, var2idx, orig_branch):
+        super(ChangeBranchObjectCallbackContext, self).__init__(problem, solver, var2idx)
+        self._orig_branch = orig_branch
+        self._branch = orig_branch
+    @property
+    def orig_branchobject(self):
+        """The branching suggested by Xpress."""
+        return self._orig_branch
+    @property
+    def branchobject(self):
+        """The branching that will be carried out.
+
+        If you set this to `None` then Xpress will carry out the branching
+        it originally intended to do.
+        """
+        return self._branch
+    @branchobject.setter
+    def branchobject(self, value):
+        self._branch = value
+    def new_object(self, nbranch):
+        """Create a new branch object that branches in the original space."""
+        return xpress.branchobj(self.problem, nbranch, isoriginal=True)
+    def map_pyomo_var(self, var):
+        """Map a Pyomo variable to an Xpress variable."""
+        return self._solver._pyomo_var_to_solver_var_map[var]
+
+class CallbackInterface(object):
+    """This is an internal class that is used for callback handling.
+
+    Callbacks are stored internally and are only applied to an Xpress
+    problem instance when `_apply` is called. The usage pattern is this:
+      with callback_interface._apply(prob):
+         prob.solve()
+    """
+
+    class EmptyManager:
+        """Context manager used in `_apply` in case there are not callbacks."""
+        def __enter__(self):
+            return self
+        def __exit__(self, type, value, traceback):
+            pass
+    
+    def __init__(self):
+        super(CallbackInterface, self).__init__()
+        self._callbacks = {
+            'message': [],
+            'optnode': [],
+            'preintsol': [],
+            'intsol': [],
+            'chgbranchobject': []
+        }
+        self._numcbs = 0
+        self._nocb = CallbackInterface.EmptyManager()
+    def _add_callback(self, name, cb, priority):
+        self._callbacks[name].append((cb, priority))
+        self._numcbs += 1
+    def _del_callback(self, name, cb):
+        self._numcbs -= len(self._callbacks[name])
+        if cb is None:
+            # Delete all
+            self._callbacks[name] = []
+        else:
+            try:
+                self._callbacks[name].remove(cb)
+            except ValueError:
+                # It is ok if the callback was not even registered
+                pass
+        self._numcbs += len(self._callbacks[name])
+    def add_message(self, messagecb, priority = 0):
+        """Add a message callback.
+        `messagecb` must be a callable that takes exactly one argument.
+        It will be invoked with an instance of `MessageCallbackContext`
+        as argument.
+        """
+        self._add_callback('message', messagecb, priority)
+    def del_message(self, messagecb = None):
+        """Remove a specific or all (`messagecb` is `None`) message callbacks."""
+        self._del_callback('message', messagecb)
+    def _message_wrapper(self, solver, var2idx):
+        return lambda prob, cb, msg, msgtype: cb(MessageCallbackContext(prob, solver, var2idx, msg, msgtype))
+
+    def add_optnode(self, optnodecb, priority = 0):
+        """Add an optnode callback.
+        `optnodecb` must be a callable that takes exactly one argument.
+        It will be invoked with an instance of `OptNodeCallbackContext`
+        as argument.
+        """
+        self._add_callback('optnode', optnodecb, priority)
+    def del_optnode(self, optnodecb = None):
+        """Remove a specific or all (`optnodecb` is `None`) optnode callbacks."""
+        self._del_callback('optnode', optnodecb)
+    def _optnode_wrapper(self, solver, var2idx):
+        def wrapper(prob, cb):
+            data = OptNodeCallbackContext(prob, solver, var2idx)
+            cb(data)
+            return data.infeas
+        return wrapper
+
+    def add_preintsol(self, preintsolcb, priority = 0):
+        """Add a preintsol callback.
+        `preintsolcb` must be a callable that takes exactly one argument.
+        It will be invoked with an instance of `PreIntSolCallbackContext`
+        as argument.
+        """
+        self._add_callback('preintsol', preintsolcb, priority)
+    def del_preintsol(self, preintsolcb = None):
+        """Remove a specific or all (`preintsolcb` is `None`) preintsol callbacks."""
+        self._del_callback('preintsol', preintsolcb)
+    def _preintsol_wrapper(self, solver, var2idx):
+        def wrapper(prob, cb, soltype, cutoff):
+            data = PreIntSolCallbackContext(prob, solver, var2idx, soltype, cutoff)
+            cb(data)
+            return (data.reject, data.cutoff)
+        return wrapper
+
+    def add_intsol(self, intsolcb, priority = 0):
+        """Add an intsol callback.
+        `intsolcb` must be a callable that takes exactly one argument.
+        It will be invoked with an instance of `IntSolCallbackContext`
+        as argument.
+        """
+        self._add_callback('intsol', intsolcb, priority)
+    def del_intsol(self, intsolcb = None):
+        """Remove a specific or all (`intsolcb` is `None`) intsol callbacks."""
+        self._del_callback('intsol', preintsolcb)
+    def _intsol_wrapper(self, solver, var2idx):
+        return lambda prob, cb: cb(IntSolCallbackContext(prob, solver, var2idx))
+
+    def add_chgbranchobject(self, chgbranchobjectcb, priority = 0):
+        """Add an chgbranchobject callback."""
+        self._add_callback('chgbranchobject', chgbranchobjectcb, priority)
+    def del_chgbranchobject(self, chgbranchobjectcb = None):
+        """Remove a specific or all (`chgbranchobjectcb` is `None`) chgbranchobject callbacks."""
+        self._del_callback('chgbranchobject', prechgbranchobjectcb)
+    def _chgbranchobject_wrapper(self, solver, var2idx):
+        def wrapper(prob, cb, orig_branch):
+            data = ChangeBranchObjectCallbackContext(prob, solver, var2idx, orig_branch)
+            cb(data)
+            return data.branchobject
+        return wrapper
+
+    def _apply(self, prob, solver):
+        """Apply callback settings.
+
+        Installs all registered callbacks into `prob`.
+        `solver` is the solver object that will be available as "solver"
+        property in the objects that are passed to the callback functions.
+        Returns a context manager that will uninstall callbacks in `__exit__`.
+        """
+        if self._numcbs == 0:
+            return self._nocb
+
+        # This is required in the callbacks to map Pyomo variables to the
+        # respective variable indices in the solver.
+        var2idx = { v: i for i, v in enumerate(prob.getVariable()) }
+
+        class Manager(object):
+            def __init__(self, prob):
+                self._prob = prob
+                self._cbs = []
+            def __enter__(self):
+                return self
+            def __exit__(self, type, value, traceback):
+                self._prob.removecbmessage(None, None)
+                self._prob.removecboptnode(None, None)
+                self._prob.removecbpreintsol(None, None)
+                self._prob.removecbintsol(None, None)
+                self._prob.removecbchgbranchobject(None, None)
+        m = Manager(prob)
+        try:
+            for cb, prio in self._callbacks['message']:
+                prob.addcbmessage(self._message_wrapper(solver, var2idx), cb, prio)
+            for cb, prio in self._callbacks['optnode']:
+                prob.addcboptnode(self._optnode_wrapper(solver, var2idx), cb, prio)
+            for cb, prio in self._callbacks['preintsol']:
+                prob.addcbpreintsol(self._preintsol_wrapper(solver, var2idx), cb, prio)
+            for cb, prio in self._callbacks['intsol']:
+                prob.addcbintsol(self._intsol_wrapper(solver, var2idx), cb, prio)
+            for cb, prio in self._callbacks['chgbranchobject']:
+                prob.addcbchgbranchobject(self._chgbranchobject_wrapper(solver, var2idx), cb, prio)
+            return m
+        except Exception:
+            t, v, b = sys.exc_info()
+            m.__exit__(t, v, b)
+            raise
 
 @SolverFactory.register('xpress_direct', doc='Direct python interface to XPRESS')
 class XpressDirect(DirectSolver):
-
+    """
+    Interface to the Xpress solver.
+
+    Note that this solver also supports interaction with the solution process
+    via callbacks. The callbacks allow notification about new solutions,
+    give opportunities to reject solution, dynamically inject cuts and
+    constraints, etc.
+
+    Callbacks are handled by the `callbacks` property which is an instance of
+    `CallbackInterface`.
+    Callbacks are registered via
+      solver.callbacks.add_preintsol(callback_function)
+      solver.callbacks.add_intsol(callback_function)
+      solver.callbacks.add_optnode(callback_function)
+    etc.
+
+    See the documentation of `CallbackInterface` for further information and
+    the Xpress specific test cases in Pyomo for some examples how to use those.
+    Also see the general Xpress documentation to understand in which contexts
+    callbacks are invoked and what you can do with them.
+    """
     _name = None
     _version = None
     XpressException = RuntimeError
@@ -130,6 +597,7 @@ def __init__(self, **kwds):
         self._range_constraints = set()
 
         self._python_api_exists = xpress_available
+        self._callbacks = CallbackInterface()
 
         # TODO: this isn't a limit of XPRESS, which implements an SLP
         #       method for NLPs. But it is a limit of *this* interface
@@ -163,6 +631,9 @@ def available(self, exception_flag=True):
                 "No Python bindings available for %s solver plugin"
                 % (type(self),))
         return bool(xpress_available)
+    @property
+    def callbacks(self):
+        return self._callbacks
 
     def _apply_solver(self):
         StaleFlagManager.mark_all_as_stale()
@@ -473,24 +944,25 @@ def _get_nlp_results(self, results, soln):
             return have_soln
 
     def _solve_model(self):
-        xprob = self._solver_model
-
-        is_mip = (xprob.attributes.mipents > 0) or (xprob.attributes.sets > 0)
-        # Check for quadratic objective or quadratic constraints. If there are
-        # any then we call nlpoptimize since that can handle non-convex
-        # quadratics as well. In case of convex quadratics it will call
-        # mipoptimize under the hood.
-        if (xprob.attributes.qelems > 0) or (xprob.attributes.qcelems > 0):
-            xprob.nlpoptimize("g" if is_mip else "")
-            self._get_results = self._get_nlp_results
-        elif is_mip:
-            xprob.mipoptimize()
-            self._get_results = self._get_mip_results
-        else:
-            xprob.lpoptimize()
-            self._get_results = self._get_lp_results
+        with self._callbacks._apply(self._solver_model, self):
+            xprob = self._solver_model
+
+            is_mip = (xprob.attributes.mipents > 0) or (xprob.attributes.sets > 0)
+            # Check for quadratic objective or quadratic constraints. If there are
+            # any then we call nlpoptimize since that can handle non-convex
+            # quadratics as well. In case of convex quadratics it will call
+            # mipoptimize under the hood.
+            if (xprob.attributes.qelems > 0) or (xprob.attributes.qcelems > 0):
+                xprob.nlpoptimize("g" if is_mip else "")
+                self._get_results = self._get_nlp_results
+            elif is_mip:
+                xprob.mipoptimize()
+                self._get_results = self._get_mip_results
+            else:
+                xprob.lpoptimize()
+                self._get_results = self._get_lp_results
         
-        self._solver_model.postsolve()
+            xprob.postsolve()
 
     def _get_expr_from_pyomo_repn(self, repn, max_degree=2):
         referenced_vars = ComponentSet()
diff --git a/pyomo/solvers/tests/checks/test_xpress_persistent.py b/pyomo/solvers/tests/checks/test_xpress_persistent.py
index 0acb9c09777..611d06c9362 100644
--- a/pyomo/solvers/tests/checks/test_xpress_persistent.py
+++ b/pyomo/solvers/tests/checks/test_xpress_persistent.py
@@ -4,6 +4,209 @@
 from pyomo.solvers.plugins.solvers.xpress_direct import xpress_available
 from pyomo.opt.results.solver import TerminationCondition, SolverStatus
 
+from pyomo.environ import value
+from pyomo.opt import SolverFactory
+from pyomo.common.collections import ComponentMap
+import random
+import math
+
+# This serves as example as well as test case. It illustrates how to use
+# callbacks with Pyomo and the Xpress solver. Of course, this is very
+# solver specific since callbacks are inherently solver specific
+class TSP:
+    """
+    Solve a MIP using cuts/constraints that are lazily separated.
+
+    We take a random instance of the symmetric TSP and solve that using
+    lazily separated constraints.
+
+    The model is based on a graph G = (V,E).
+    We have one binary variable x[e] for each edge e in E. That variable
+    is set to 1 if edge e is selected in the optimal tour and 0 otherwise.
+
+    The model contains only these explicit constraint:
+        for each v in V: sum(u in V : u != v) x[uv] == 1
+        for each v in V: sum(u in V : u != v) x[vu] == 1
+
+    This states that each node must be entered and exited exactly once in the
+    optimal tour.
+
+    The above constraints ensures that the selected edges form tours. However,
+    it allows multiple tours, also known as subtours. So we need a constraint
+    that requires that there is only one tour (which then necessarily hits
+    all the nodes). For this we just create no-good constraints for subtours:
+    If S is a set of edges that forms a subtour, then we add constraint
+        sum(e in S) x[e] <= |S|-1
+
+    Since there are exponentially many subtours in a graph, this constraint
+    is not stated explicitly. Instead we check for any solution that the
+    optimizer finds, whether it satisfies the subtour elimination constraint.
+    If it does then we accept the solution. Otherwise we reject the solution
+    and augment the model by the violated subtour elimination constraint.
+
+    This lazy addition of constraints is implemented using two callbacks:
+    - a preintsol callback that rejects any solution that violates a
+      subtour elimination constraint,
+    - an optnode callback that injects any violated subtour elimination
+      constraints.
+
+    An important thing to note about this strategy is that dual reductions
+    have to be disabled. Since the optimizer does not see the whole model
+    (subtour elimination constraints are only generated on the fly), dual
+    reductions may cut off the optimal solution.
+    """
+    def __init__(self, nodes, seed=0):
+        """Construct a new random instance with the given seed."""
+        self._nodes = nodes  # Number of nodes/cities in the instance.
+        self._nodex = [0.0] * nodes # X coordinate of nodes.
+        self._nodey = [0.0] * nodes # Y coordinate of nodes.
+
+        random.seed(seed)
+        for i in range(nodes):
+            self._nodex[i] = 4.0 * random.random()
+            self._nodey[i] = 4.0 * random.random()
+
+    def distance(self, u, v):
+        """Get the distance between two nodes."""
+        return math.sqrt((self._nodex[u] - self._nodex[v]) ** 2 +
+                         (self._nodey[u] - self._nodey[v]) ** 2)
+
+    def find_tour(self, sol, prev=None):
+        """Find the tour rooted at the first city in a solution.
+        """
+        if prev is None:
+            prev = dict()
+        tour = set()
+        x = self._model.x
+
+        u = 1
+        used = 0
+        print('1', end='')
+        while True:
+            for v in self._model.cities:
+                if u == v:
+                    continue # no self-loops
+                elif v in prev:
+                    continue # node already on tour
+                elif sol[x[u, v]] < 0.5:
+                    continue # edge not selected in solution
+                elif (u, v) in tour:
+                    continue # edge already on tour
+                else:
+                    print(' -> %d' % v, end='')
+                    tour.add((u, v))
+                    prev[v] = u
+                    used += 1;
+                    u = v;
+                    break
+            if u == 1:
+                break
+        print()
+        return used
+
+    def preintsol(self, data):
+        """Integer solution check callback."""
+        print("Checking feasible solution ...")
+
+        # Get current solution and check whether it is feasible
+        used = self.find_tour(data.candidate_sol)
+        print("Solution is ", end='')
+        if used < len(self._model.cities):
+            print('infeasible (%d edges)' % used)
+            data.reject = True
+        else:
+            print('feasible with length %f' % data.candidate_obj)
+
+    def optnode(self, data):
+        """Optimal node callback.
+        This callback is invoked after the LP relaxation of a node is solved.
+        This is where we can inject additional constraints as cuts.
+        """
+        # Only separate constraints on nodes that are integer feasible.
+        if data.attributes.mipinfeas != 0:
+            return
+
+        # Get the current solution
+        sol = data.getlpsol()
+
+        # Get the tour starting at the first city and check whether it covers
+        # all nodes. If it does not then it is infeasible and we must
+        # generate a subtour elimination constraint.
+        prev = dict()
+        used = self.find_tour(sol, prev)
+        if used < len(self._model.cities):
+            # The tour is too short. Get the edges on the tour and add a
+            # subtour elimination constraint
+            lhs = sum(self._model.x[u, prev[u]] for u in self._model.cities if u in prev)
+            data.addcut(lhs <= (used - 1))
+            data.infeas = True
+
+    def create_initial_tour(self):
+        """Create a feasible tour and add this as initial MIP solution."""
+        for u in self._model.cities:
+            for v in self._model.cities:
+                # A starting solution in Pyomo is set by assigning the values
+                # to the variables and then passing `warmstart=True` to the
+                # `solve()` call.
+                if v == u + 1 or (u == self._nodes and v == 1):
+                    self._model.x[u, v] = 1.0
+                else:
+                    self._model.x[u, v] = 0.0
+
+    def solve(self):
+        """Solve the TSP represented by this instance."""
+        self._model = pe.ConcreteModel()
+        self._model.cities = pe.RangeSet(self._nodes)
+        # Create variables. We create one variable for each edge in
+        # the complete directed graph. x[u,v] is set to 1 if the tour goes
+        # from u to v, otherwise it is set to 0.
+        # All variables are binary.
+        self._model.x = pe.Var(self._model.cities * self._model.cities,
+                               within=pe.Binary)
+        self._model.cons = pe.ConstraintList()
+
+        # Do not allow self loops.
+        # We could have skipped creating the variables but fixing them to
+        # 0 here is slightly easier.
+        for u in self._model.cities:
+            self._model.cons.add(self._model.x[u, u] <= 0.0)
+
+        # Objective function.
+        obj = sum(self._model.x[u, v] * self.distance(u-1, v-1) for u in self._model.cities for v in self._model.cities)
+        self._model.obj = pe.Objective(expr=obj)
+
+        # Constraint: Each node must be exited and entered exactly once.
+        for u in self._model.cities:
+            self._model.cons.add(sum(self._model.x[u, v] for v in self._model.cities if v != u) == 1)
+            self._model.cons.add(sum(self._model.x[v, u] for v in self._model.cities if v != u) == 1)
+
+        # Create a starting solution.
+        # This is optional but having a feasible solution available right
+        # from the beginning can improve optimizer performance.
+        self.create_initial_tour()
+
+        # We don't have all constraints explicitly in the matrix, hence
+        # we must disable dual reductions. Otherwise MIP presolve may
+        # cut off the optimal solution.
+        opt = SolverFactory('xpress_direct')
+        opt.options['mipdualreductions'] = 0
+
+        # Add a callback that rejects solutions that do not satisfy
+        # the subtour constraints.
+        opt.callbacks.add_preintsol(self.preintsol)
+
+        # Add a callback that separates subtour elimination constraints
+        opt.callbacks.add_optnode(self.optnode)
+
+        opt.solve(self._model, tee=True, warmstart=True)
+
+        # Print the optimal tour.
+        print("Tour with length %f:" % value(self._model.obj))
+        self.find_tour(ComponentMap([(self._model.x[e], self._model.x[e].value) for e in self._model.x]))
+
+        self._model = None # cleanup
+
+
 class TestXpressPersistent(unittest.TestCase):
     @unittest.skipIf(not xpress_available, "xpress is not available")
     def test_basics(self):
@@ -263,6 +466,177 @@ def test_add_column_exceptions(self):
         # var already in solver model
         self.assertRaises(RuntimeError, opt.add_column, m, m.y, -2, [m.c], [1])
 
+    def _markshare(self):
+        """Create a model that is non-trivial to solve.
+        The returned model has two variables: `x` and `s`. It also has an
+        objective function that is stored in `obj`.
+        """
+        model = pe.ConcreteModel()
+        model.X = pe.RangeSet(50)
+        model.S = pe.RangeSet(6)
+        model.x = pe.Var(model.X, within=pe.Binary)
+        x = model.x
+        model.s = pe.Var(model.S, bounds = (0, None))
+        s = model.s
+        model.obj = pe.Objective(expr=s[1] + s[2] + s[3] + s[4] + s[5] + s[6])
+        model.cons = pe.ConstraintList()
+        
+        model.cons.add(s[1] + 25*x[1] + 35*x[2] + 14*x[3] + 76*x[4] + 58*x[5] + 10*x[6] + 20*x[7]
+                       + 51*x[8] + 58*x[9] + x[10] + 35*x[11] + 40*x[12] + 65*x[13] + 59*x[14] + 24*x[15]
+                       + 44*x[16] + x[17] + 93*x[18] + 24*x[19] + 68*x[20] + 38*x[21] + 64*x[22] + 93*x[23]
+                       + 14*x[24] + 83*x[25] + 6*x[26] + 58*x[27] + 14*x[28] + 71*x[29] + 17*x[30]
+                       + 18*x[31] + 8*x[32] + 57*x[33] + 48*x[34] + 35*x[35] + 13*x[36] + 47*x[37]
+                       + 46*x[38] + 8*x[39] + 82*x[40] + 51*x[41] + 49*x[42] + 85*x[43] + 66*x[44]
+                       + 45*x[45] + 99*x[46] + 21*x[47] + 75*x[48] + 78*x[49] + 43*x[50] == 1116)
+        model.cons.add(s[2] + 97*x[1] + 64*x[2] + 24*x[3] + 63*x[4] + 58*x[5] + 45*x[6] + 20*x[7]
+                       + 71*x[8] + 32*x[9] + 7*x[10] + 28*x[11] + 77*x[12] + 95*x[13] + 96*x[14]
+                       + 70*x[15] + 22*x[16] + 93*x[17] + 32*x[18] + 17*x[19] + 56*x[20] + 74*x[21]
+                       + 62*x[22] + 94*x[23] + 9*x[24] + 92*x[25] + 90*x[26] + 40*x[27] + 45*x[28]
+                       + 84*x[29] + 62*x[30] + 62*x[31] + 34*x[32] + 21*x[33] + 2*x[34] + 75*x[35]
+                       + 42*x[36] + 75*x[37] + 29*x[38] + 4*x[39] + 64*x[40] + 80*x[41] + 17*x[42]
+                       + 55*x[43] + 73*x[44] + 23*x[45] + 13*x[46] + 91*x[47] + 70*x[48] + 73*x[49]
+                       + 28*x[50] == 1325)
+        model.cons.add(s[3] + 95*x[1] + 71*x[2] + 19*x[3] + 15*x[4] + 66*x[5] + 76*x[6] + 4*x[7]
+                       + 50*x[8] + 50*x[9] + 97*x[10] + 83*x[11] + 14*x[12] + 27*x[13] + 14*x[14]
+                       + 34*x[15] + 9*x[16] + 99*x[17] + 62*x[18] + 92*x[19] + 39*x[20] + 56*x[21]
+                       + 53*x[22] + 91*x[23] + 81*x[24] + 46*x[25] + 94*x[26] + 76*x[27] + 53*x[28]
+                       + 58*x[29] + 23*x[30] + 15*x[31] + 63*x[32] + 2*x[33] + 31*x[34] + 55*x[35]
+                       + 71*x[36] + 97*x[37] + 71*x[38] + 55*x[39] + 8*x[40] + 57*x[41] + 14*x[42]
+                       + 76*x[43] + x[44] + 46*x[45] + 87*x[46] + 22*x[47] + 97*x[48] + 99*x[49] + 92*x[50]
+                       == 1353)
+        model.cons.add(s[4] + x[1] + 27*x[2] + 46*x[3] + 48*x[4] + 66*x[5] + 58*x[6] + 52*x[7] + 6*x[8]
+                       + 14*x[9] + 26*x[10] + 55*x[11] + 61*x[12] + 60*x[13] + 3*x[14] + 33*x[15]
+                       + 99*x[16] + 36*x[17] + 55*x[18] + 70*x[19] + 73*x[20] + 70*x[21] + 38*x[22]
+                       + 66*x[23] + 39*x[24] + 43*x[25] + 63*x[26] + 88*x[27] + 47*x[28] + 18*x[29]
+                       + 73*x[30] + 40*x[31] + 91*x[32] + 96*x[33] + 49*x[34] + 13*x[35] + 27*x[36]
+                       + 22*x[37] + 71*x[38] + 99*x[39] + 66*x[40] + 57*x[41] + x[42] + 54*x[43] + 35*x[44]
+                       + 52*x[45] + 66*x[46] + 26*x[47] + x[48] + 26*x[49] + 12*x[50] == 1169)
+        model.cons.add(s[5] + 3*x[1] + 94*x[2] + 51*x[3] + 4*x[4] + 25*x[5] + 46*x[6] + 30*x[7]
+                       + 2*x[8] + 89*x[9] + 65*x[10] + 28*x[11] + 46*x[12] + 36*x[13] + 53*x[14]
+                       + 30*x[15] + 73*x[16] + 37*x[17] + 60*x[18] + 21*x[19] + 41*x[20] + 2*x[21]
+                       + 21*x[22] + 93*x[23] + 82*x[24] + 16*x[25] + 97*x[26] + 75*x[27] + 50*x[28]
+                       + 13*x[29] + 43*x[30] + 45*x[31] + 64*x[32] + 78*x[33] + 78*x[34] + 6*x[35]
+                       + 35*x[36] + 72*x[37] + 31*x[38] + 28*x[39] + 56*x[40] + 60*x[41] + 23*x[42]
+                       + 70*x[43] + 46*x[44] + 88*x[45] + 20*x[46] + 69*x[47] + 13*x[48] + 40*x[49]
+                       + 73*x[50] == 1160)
+        model.cons.add(s[6] + 69*x[1] + 72*x[2] + 94*x[3] + 56*x[4] + 90*x[5] + 20*x[6] + 56*x[7]
+                       + 50*x[8] + 79*x[9] + 59*x[10] + 36*x[11] + 24*x[12] + 42*x[13] + 9*x[14]
+                       + 29*x[15] + 68*x[16] + 10*x[17] + x[18] + 44*x[19] + 74*x[20] + 61*x[21] + 37*x[22]
+                       + 71*x[23] + 63*x[24] + 44*x[25] + 77*x[26] + 57*x[27] + 46*x[28] + 51*x[29]
+                       + 43*x[30] + 4*x[31] + 85*x[32] + 59*x[33] + 7*x[34] + 25*x[35] + 46*x[36] + 25*x[37]
+                       + 70*x[38] + 78*x[39] + 88*x[40] + 20*x[41] + 40*x[42] + 40*x[43] + 16*x[44]
+                       + 3*x[45] + 3*x[46] + 5*x[47] + 77*x[48] + 88*x[49] + 16*x[50] == 1163)
+        
+        return model
+        
+
+    @unittest.skipIf(not xpress_available, "xpress is not available")
+    def test_callbacks_01(self):
+        """Simple callback test.
+
+        Tests that optnode, preintsol, intsol callbacks are invoked.
+        Also tests that information between preintsol and intsol callbacks
+        is consistent.
+        """
+        model = self._markshare()
+        opt = pe.SolverFactory('xpress_direct')
+        opt.options['MAXNODE'] = 5
+        opt.options['THREADS'] = 1 # for interaction between preintsol and intsol
+        test = self
+
+        lastnode = [0]
+        noptnode = [0]
+        def optnode(data):
+            try:
+                noptnode[0] += 1
+                node = data.attributes.nodes
+                test.assertGreaterEqual(node, lastnode[0])
+                lastnode[0] = node
+            except Exception as ex:
+                print('optnode:', ex)
+                raise ex
+        opt.callbacks.add_optnode(optnode)
+
+        announced = [None]
+        npreintsol = [0]
+        def preintsol(data):
+            try:
+                npreintsol[0] += 1
+                test.assertIsNone(announced[0])
+                test.assertGreater(data.attributes.mipobjval,
+                                   data.candidate_obj)
+                announced[0] = (data.candidate_obj, data.candidate_sol)
+            except Exception as ex:
+                print('preintsol:', ex)
+                raise ex
+        opt.callbacks.add_preintsol(preintsol)
+
+        nintsol = [0]
+        def intsol(data):
+            try:
+                nintsol[0] += 1
+                self.assertIsNotNone(announced[0])
+                obj, x = announced[0]
+                announced[0] = None
+                self.assertEqual(data.objval, obj)
+                sol = data.solution
+                for p in sol:
+                    self.assertEqual(sol[p], x[p])
+            except Exception as ex:
+                print('intsol:', ex)
+                raise ex
+        opt.callbacks.add_intsol(intsol)
+
+        opt.solve(model, tee=True)
+
+        self.assertGreater(noptnode[0], 0)
+        self.assertGreater(npreintsol[0], 0)
+        self.assertGreater(nintsol[0], 0)
+
+    @unittest.skipIf(not xpress_available, "xpress is not available")
+    def test_callbacks_02(self):
+        """Test branching callback by doing most fractional branching on markshare."""
+
+        model = self._markshare()
+        opt = pe.SolverFactory('xpress_direct')
+        opt.options['MAXNODE'] = 5
+        test = self
+        
+        called = [0]
+        def chgbranchobject(data):
+            try:
+                test.assertIsNotNone(data.branchobject)
+                test.assertEqual(data.branchobject, data.orig_branchobject)
+                called[0] += 1
+                sol = data.getlpsol()
+                maxfrac = 0.0
+                maxvar = None
+                for var in sol:
+                    if var.domain == pe.Binary:
+                        frac = abs(round(sol[var]) - sol.var)
+                        if frac > maxfrac:
+                            maxfrac = frac
+                            maxvar = var
+                test.assertIsNotNone(maxvar)
+                if maxvar is not None:
+                    b = data.new_object(2)
+                    b.addbounds(0, ['U'], [data.map_pyomo_var(maxvar)], [0])
+                    b.addbounds(1, ['L'], [data.map_pyomo_var(maxvar)], [1])
+                    data.branchobject = b
+            except Exception as ex:
+                print('chgbranchobject:', ex)
+                raise ex
+        opt.callbacks.add_chgbranchobject(chgbranchobject)
+
+        opt.solve(model, tee=True)
+
+        self.assertGreater(called[0], 0)
+
+    @unittest.skipIf(not xpress_available, "xpress is not available")
+    def test_callbacks_03(self):
+        """Test the TSP example."""
+        TSP(10).solve()
+
     @unittest.skipIf(not xpress_available, "xpress is not available")
     def test_nonconvexqp_locally_optimal(self):
         """Test non-convex QP for which xpress_direct should find a locally