diff --git a/src/python/miplearn/solvers/cplex.py b/src/python/miplearn/solvers/cplex.py
index 50d0a06..7041e9b 100644
--- a/src/python/miplearn/solvers/cplex.py
+++ b/src/python/miplearn/solvers/cplex.py
@@ -1,15 +1,14 @@
 #  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
 #  Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
 #  Released under the modified BSD license. See COPYING.md for more details.
-import re
 
 import pyomo.environ as pe
 from scipy.stats import randint
 
-from .internal import InternalSolver
+from .pyomo import PyomoSolver
 
 
-class CPLEXSolver(InternalSolver):
+class CPLEXSolver(PyomoSolver):
     def __init__(self, options=None):
         """
         Creates a new CPLEXSolver.
diff --git a/src/python/miplearn/solvers/gurobi.py b/src/python/miplearn/solvers/gurobi.py
index 139b767..0164d7a 100644
--- a/src/python/miplearn/solvers/gurobi.py
+++ b/src/python/miplearn/solvers/gurobi.py
@@ -8,13 +8,14 @@ from io import StringIO
 
 import pyomo.environ as pe
 from miplearn.solvers import RedirectOutput
-from miplearn.solvers.internal import InternalSolver
 from scipy.stats import randint
 
+from .pyomo import PyomoSolver
+
 logger = logging.getLogger(__name__)
 
 
-class GurobiSolver(InternalSolver):
+class GurobiSolver(PyomoSolver):
     def __init__(self,
                  use_lazy_callbacks=True,
                  options=None):
diff --git a/src/python/miplearn/solvers/internal.py b/src/python/miplearn/solvers/internal.py
index 77deb6f..dd54dec 100644
--- a/src/python/miplearn/solvers/internal.py
+++ b/src/python/miplearn/solvers/internal.py
@@ -20,29 +20,14 @@ logger = logging.getLogger(__name__)
 
 class InternalSolver(ABC):
     """
-    The MIP solver used internaly by LearningSolver.
-
-    Attributes
-    ----------
-    instance: miplearn.Instance
-        The MIPLearn instance currently loaded to the solver
-    model: pyomo.core.ConcreteModel
-        The Pyomo model currently loaded on the solver
+    Abstract class representing the MIP solver used internally by LearningSolver.
     """
 
-    def __init__(self):
-        self.instance = None
-        self.model = None
-        self._all_vars = None
-        self._bin_vars = None
-        self._is_warm_start_available = False
-        self._pyomo_solver = None
-        self._obj_sense = None
-        self._varname_to_var = {}
-
+    @abstractmethod
     def solve_lp(self, tee=False):
         """
-        Solves the LP relaxation of the currently loaded instance.
+        Solves the LP relaxation of the currently loaded instance. After this
+        method finishes, the solution can be retrieved by calling `get_solution`.
 
         Parameters
         ----------
@@ -55,20 +40,9 @@ class InternalSolver(ABC):
             A dictionary of solver statistics containing the following keys:
             "Optimal value".
         """
-        for var in self._bin_vars:
-            lb, ub = var.bounds
-            var.setlb(lb)
-            var.setub(ub)
-            var.domain = pyomo.core.base.set_types.Reals
-            self._pyomo_solver.update_var(var)
-        results = self._pyomo_solver.solve(tee=tee)
-        for var in self._bin_vars:
-            var.domain = pyomo.core.base.set_types.Binary
-            self._pyomo_solver.update_var(var)
-        return {
-            "Optimal value": results["Problem"][0]["Lower bound"],
-        }
+        pass
 
+    @abstractmethod
     def get_solution(self):
         """
         Returns current solution found by the solver.
@@ -80,45 +54,28 @@ class InternalSolver(ABC):
         The solution is a dictionary `sol`, where the optimal value of `var[idx]`
         is given by `sol[var][idx]`.
         """
-        solution = {}
-        for var in self.model.component_objects(Var):
-            solution[str(var)] = {}
-            for index in var:
-                solution[str(var)][index] = var[index].value
-        return solution
+        pass
 
+    @abstractmethod
     def set_warm_start(self, solution):
         """
         Sets the warm start to be used by the solver.
 
         The solution should be a dictionary following the same format as the
-        one produced by `get_solution`. Only one warm start is currently
-        supported. Calling this function when a warm start already exists will
+        one produced by `get_solution`. Only one warm start is supported.
+        Calling this function when a warm start already exists will
         remove the previous warm start.
         """
-        self.clear_warm_start()
-        count_total, count_fixed = 0, 0
-        for var_name in solution:
-            var = self._varname_to_var[var_name]
-            for index in solution[var_name]:
-                count_total += 1
-                var[index].value = solution[var_name][index]
-                if solution[var_name][index] is not None:
-                    count_fixed += 1
-        if count_fixed > 0:
-            self._is_warm_start_available = True
-        logger.info("Setting start values for %d variables (out of %d)" %
-                    (count_fixed, count_total))
+        pass
 
+    @abstractmethod
     def clear_warm_start(self):
         """
         Removes any existing warm start from the solver.
         """
-        for var in self._all_vars:
-            if not var.fixed:
-                var.value = None
-        self._is_warm_start_available = False
+        pass
 
+    @abstractmethod
     def set_instance(self, instance, model=None):
         """
         Loads the given instance into the solver.
@@ -127,34 +84,14 @@ class InternalSolver(ABC):
         ----------
         instance: miplearn.Instance
             The instance to be loaded.
-        model: pyomo.core.ConcreteModel
-            The corresponding Pyomo model. If not provided, it will be
-            generated by calling `instance.to_model()`.
+        model:
+            The concrete optimization model corresponding to this instance
+            (e.g. JuMP.Model or pyomo.core.ConcreteModel). If not provided,
+            it will be generated by calling `instance.to_model()`.
         """
-        if model is None:
-            model = instance.to_model()
-        assert isinstance(instance, Instance)
-        assert isinstance(model, pe.ConcreteModel)
-        self.instance = instance
-        self.model = model
-        self._pyomo_solver.set_instance(model)
-
-        # Update objective sense
-        self._obj_sense = "max"
-        if self._pyomo_solver._objective.sense == pyomo.core.kernel.objective.minimize:
-            self._obj_sense = "min"
-
-        # Update variables
-        self._all_vars = []
-        self._bin_vars = []
-        self._varname_to_var = {}
-        for var in model.component_objects(Var):
-            self._varname_to_var[var.name] = var
-            for idx in var:
-                self._all_vars += [var[idx]]
-                if var[idx].domain == pyomo.core.base.set_types.Binary:
-                    self._bin_vars += [var[idx]]
+        pass
 
+    @abstractmethod
     def fix(self, solution):
         """
         Fixes the values of a subset of decision variables.
@@ -163,28 +100,20 @@ class InternalSolver(ABC):
         `get_solution`. Missing values in the solution indicate variables
         that should be left free.
         """
-        count_total, count_fixed = 0, 0
-        for varname in solution:
-            for index in solution[varname]:
-                var = self._varname_to_var[varname]
-                count_total += 1
-                if solution[varname][index] is None:
-                    continue
-                count_fixed += 1
-                var[index].fix(solution[varname][index])
-                self._pyomo_solver.update_var(var[index])
-        logger.info("Fixing values for %d variables (out of %d)" %
-                    (count_fixed, count_total))
+        pass
 
+    @abstractmethod
     def add_constraint(self, constraint):
         """
         Adds a single constraint to the model.
         """
-        self._pyomo_solver.add_constraint(constraint)
+        pass
 
+    @abstractmethod
     def solve(self, tee=False):
         """
-        Solves the currently loaded instance.
+        Solves the currently loaded instance. After this method finishes,
+        the best solution found can be retrieved by calling `get_solution`.
 
         Parameters
         ----------
@@ -198,107 +127,22 @@ class InternalSolver(ABC):
             "Lower bound", "Upper bound", "Wallclock time", "Nodes", "Sense",
              "Log" and "Warm start value".
         """
-        total_wallclock_time = 0
-        streams = [StringIO()]
-        if tee:
-            streams += [sys.stdout]
-        self.instance.found_violations = []
-        while True:
-            logger.debug("Solving MIP...")
-            with RedirectOutput(streams):
-                results = self._pyomo_solver.solve(tee=True,
-                                                   warmstart=self._is_warm_start_available)
-            total_wallclock_time += results["Solver"][0]["Wallclock time"]
-            if not hasattr(self.instance, "find_violations"):
-                break
-            logger.debug("Finding violated constraints...")
-            violations = self.instance.find_violations(self.model)
-            if len(violations) == 0:
-                break
-            self.instance.found_violations += violations
-            logger.debug("    %d violations found" % len(violations))
-            for v in violations:
-                cut = self.instance.build_lazy_constraint(self.model, v)
-                self.add_constraint(cut)
-
-        log = streams[0].getvalue()
-        return {
-            "Lower bound": results["Problem"][0]["Lower bound"],
-            "Upper bound": results["Problem"][0]["Upper bound"],
-            "Wallclock time": total_wallclock_time,
-            "Nodes": self._extract_node_count(log),
-            "Sense": self._obj_sense,
-            "Log": log,
-            "Warm start value": self._extract_warm_start_value(log),
-        }
-
-    @staticmethod
-    def __extract(log, regexp, default=None):
-        value = default
-        for line in log.splitlines():
-            matches = re.findall(regexp, line)
-            if len(matches) == 0:
-                continue
-            value = matches[0]
-        return value
-
-    def _extract_warm_start_value(self, log):
-        """
-        Extracts and returns the objective value of the user-provided MIP start
-        from the provided solver log. If more than one value is found, returns
-        the last one. If no value is present in the logs, returns None.
-        """
-        value = self.__extract(log, self._get_warm_start_regexp())
-        if value is not None:
-            value = float(value)
-        return value
-
-    def _extract_node_count(self, log):
-        """
-        Extracts and returns the number of explored branch-and-bound nodes.
-        """
-        return int(self.__extract(log,
-                                  self._get_node_count_regexp(),
-                                  default=1))
-
-    def set_threads(self, threads):
-        key = self._get_threads_option_name()
-        self._pyomo_solver.options[key] = threads
-
-    def set_time_limit(self, time_limit):
-        key = self._get_time_limit_option_name()
-        self._pyomo_solver.options[key] = time_limit
-
-    def set_node_limit(self, node_limit):
-        key = self._get_node_limit_option_name()
-        self._pyomo_solver.options[key] = node_limit
-
-    def set_gap_tolerance(self, gap_tolerance):
-        key = self._get_gap_tolerance_option_name()
-        self._pyomo_solver.options[key] = gap_tolerance
-
-    @abstractmethod
-    def _get_warm_start_regexp(self):
         pass
 
     @abstractmethod
-    def _get_node_count_regexp(self):
-        pass
-
-    @abstractmethod
-    def _get_threads_option_name(self):
+    def set_threads(self, threads):
         pass
 
     @abstractmethod
-    def _get_time_limit_option_name(self):
+    def set_time_limit(self, time_limit):
         pass
 
     @abstractmethod
-    def _get_node_limit_option_name(self):
+    def set_node_limit(self, node_limit):
         pass
 
     @abstractmethod
-    def _get_gap_tolerance_option_name(self):
+    def set_gap_tolerance(self, gap_tolerance):
         pass
 
 
diff --git a/src/python/miplearn/solvers/learning.py b/src/python/miplearn/solvers/learning.py
index ecfd0ad..2b4f858 100644
--- a/src/python/miplearn/solvers/learning.py
+++ b/src/python/miplearn/solvers/learning.py
@@ -139,7 +139,7 @@ class LearningSolver:
 
         self.tee = tee
         self.internal_solver = self._create_internal_solver()
-        self.internal_solver.set_instance(instance, model=model)
+        self.internal_solver.set_instance(instance, model)
 
         logger.debug("Solving LP relaxation...")
         results = self.internal_solver.solve_lp(tee=tee)
diff --git a/src/python/miplearn/solvers/pyomo.py b/src/python/miplearn/solvers/pyomo.py
new file mode 100644
index 0000000..a8bcb99
--- /dev/null
+++ b/src/python/miplearn/solvers/pyomo.py
@@ -0,0 +1,306 @@
+#  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
+#  Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
+#  Released under the modified BSD license. See COPYING.md for more details.
+
+import logging
+import re
+import sys
+from abc import ABC, abstractmethod
+from io import StringIO
+
+import pyomo
+import pyomo.environ as pe
+from pyomo.core import Var
+
+from . import RedirectOutput
+from .internal import InternalSolver
+from ..instance import Instance
+
+logger = logging.getLogger(__name__)
+
+
+class PyomoSolver(InternalSolver):
+    """
+    Base class for all Pyomo-based InternalSolvers.
+
+    Attributes
+    ----------
+    instance: miplearn.Instance
+        The MIPLearn instance currently loaded to the solver
+    model: pyomo.core.ConcreteModel
+        The Pyomo model currently loaded on the solver
+    """
+
+    def __init__(self):
+        self.instance = None
+        self.model = None
+        self._all_vars = None
+        self._bin_vars = None
+        self._is_warm_start_available = False
+        self._pyomo_solver = None
+        self._obj_sense = None
+        self._varname_to_var = {}
+
+    def solve_lp(self, tee=False):
+        """
+        Solves the LP relaxation of the currently loaded instance.
+
+        Parameters
+        ----------
+        tee: bool
+            If true, prints the solver log to the screen.
+
+        Returns
+        -------
+        dict
+            A dictionary of solver statistics containing the following keys:
+            "Optimal value".
+        """
+        for var in self._bin_vars:
+            lb, ub = var.bounds
+            var.setlb(lb)
+            var.setub(ub)
+            var.domain = pyomo.core.base.set_types.Reals
+            self._pyomo_solver.update_var(var)
+        results = self._pyomo_solver.solve(tee=tee)
+        for var in self._bin_vars:
+            var.domain = pyomo.core.base.set_types.Binary
+            self._pyomo_solver.update_var(var)
+        return {
+            "Optimal value": results["Problem"][0]["Lower bound"],
+        }
+
+    def get_solution(self):
+        """
+        Returns current solution found by the solver.
+
+        If called after `solve`, returns the best primal solution found during
+        the search. If called after `solve_lp`, returns the optimal solution
+        to the LP relaxation.
+
+        The solution is a dictionary `sol`, where the optimal value of `var[idx]`
+        is given by `sol[var][idx]`.
+        """
+        solution = {}
+        for var in self.model.component_objects(Var):
+            solution[str(var)] = {}
+            for index in var:
+                solution[str(var)][index] = var[index].value
+        return solution
+
+    def set_warm_start(self, solution):
+        """
+        Sets the warm start to be used by the solver.
+
+        The solution should be a dictionary following the same format as the
+        one produced by `get_solution`. Only one warm start is currently
+        supported. Calling this function when a warm start already exists will
+        remove the previous warm start.
+        """
+        self.clear_warm_start()
+        count_total, count_fixed = 0, 0
+        for var_name in solution:
+            var = self._varname_to_var[var_name]
+            for index in solution[var_name]:
+                count_total += 1
+                var[index].value = solution[var_name][index]
+                if solution[var_name][index] is not None:
+                    count_fixed += 1
+        if count_fixed > 0:
+            self._is_warm_start_available = True
+        logger.info("Setting start values for %d variables (out of %d)" %
+                    (count_fixed, count_total))
+
+    def clear_warm_start(self):
+        """
+        Removes any existing warm start from the solver.
+        """
+        for var in self._all_vars:
+            if not var.fixed:
+                var.value = None
+        self._is_warm_start_available = False
+
+    def set_instance(self, instance, model=None):
+        """
+        Loads the given instance into the solver.
+
+        Parameters
+        ----------
+        instance: miplearn.Instance
+            The instance to be loaded.
+        model: pyomo.core.ConcreteModel
+            The corresponding Pyomo model. If not provided, it will be
+            generated by calling `instance.to_model()`.
+        """
+        if model is None:
+            model = instance.to_model()
+        assert isinstance(instance, Instance)
+        assert isinstance(model, pe.ConcreteModel)
+        self.instance = instance
+        self.model = model
+        self._pyomo_solver.set_instance(model)
+
+        # Update objective sense
+        self._obj_sense = "max"
+        if self._pyomo_solver._objective.sense == pyomo.core.kernel.objective.minimize:
+            self._obj_sense = "min"
+
+        # Update variables
+        self._all_vars = []
+        self._bin_vars = []
+        self._varname_to_var = {}
+        for var in model.component_objects(Var):
+            self._varname_to_var[var.name] = var
+            for idx in var:
+                self._all_vars += [var[idx]]
+                if var[idx].domain == pyomo.core.base.set_types.Binary:
+                    self._bin_vars += [var[idx]]
+
+    def fix(self, solution):
+        """
+        Fixes the values of a subset of decision variables.
+
+        The values should be provided in the dictionary format generated by
+        `get_solution`. Missing values in the solution indicate variables
+        that should be left free.
+        """
+        count_total, count_fixed = 0, 0
+        for varname in solution:
+            for index in solution[varname]:
+                var = self._varname_to_var[varname]
+                count_total += 1
+                if solution[varname][index] is None:
+                    continue
+                count_fixed += 1
+                var[index].fix(solution[varname][index])
+                self._pyomo_solver.update_var(var[index])
+        logger.info("Fixing values for %d variables (out of %d)" %
+                    (count_fixed, count_total))
+
+    def add_constraint(self, constraint):
+        """
+        Adds a single constraint to the model.
+        """
+        self._pyomo_solver.add_constraint(constraint)
+
+    def solve(self, tee=False):
+        """
+        Solves the currently loaded instance.
+
+        Parameters
+        ----------
+        tee: bool
+            If true, prints the solver log to the screen.
+
+        Returns
+        -------
+        dict
+            A dictionary of solver statistics containing the following keys:
+            "Lower bound", "Upper bound", "Wallclock time", "Nodes", "Sense",
+             "Log" and "Warm start value".
+        """
+        total_wallclock_time = 0
+        streams = [StringIO()]
+        if tee:
+            streams += [sys.stdout]
+        self.instance.found_violations = []
+        while True:
+            logger.debug("Solving MIP...")
+            with RedirectOutput(streams):
+                results = self._pyomo_solver.solve(tee=True,
+                                                   warmstart=self._is_warm_start_available)
+            total_wallclock_time += results["Solver"][0]["Wallclock time"]
+            if not hasattr(self.instance, "find_violations"):
+                break
+            logger.debug("Finding violated constraints...")
+            violations = self.instance.find_violations(self.model)
+            if len(violations) == 0:
+                break
+            self.instance.found_violations += violations
+            logger.debug("    %d violations found" % len(violations))
+            for v in violations:
+                cut = self.instance.build_lazy_constraint(self.model, v)
+                self.add_constraint(cut)
+
+        log = streams[0].getvalue()
+        return {
+            "Lower bound": results["Problem"][0]["Lower bound"],
+            "Upper bound": results["Problem"][0]["Upper bound"],
+            "Wallclock time": total_wallclock_time,
+            "Nodes": self._extract_node_count(log),
+            "Sense": self._obj_sense,
+            "Log": log,
+            "Warm start value": self._extract_warm_start_value(log),
+        }
+
+    @staticmethod
+    def __extract(log, regexp, default=None):
+        value = default
+        for line in log.splitlines():
+            matches = re.findall(regexp, line)
+            if len(matches) == 0:
+                continue
+            value = matches[0]
+        return value
+
+    def _extract_warm_start_value(self, log):
+        """
+        Extracts and returns the objective value of the user-provided MIP start
+        from the provided solver log. If more than one value is found, returns
+        the last one. If no value is present in the logs, returns None.
+        """
+        value = self.__extract(log, self._get_warm_start_regexp())
+        if value is not None:
+            value = float(value)
+        return value
+
+    def _extract_node_count(self, log):
+        """
+        Extracts and returns the number of explored branch-and-bound nodes.
+        """
+        return int(self.__extract(log,
+                                  self._get_node_count_regexp(),
+                                  default=1))
+
+    def set_threads(self, threads):
+        key = self._get_threads_option_name()
+        self._pyomo_solver.options[key] = threads
+
+    def set_time_limit(self, time_limit):
+        key = self._get_time_limit_option_name()
+        self._pyomo_solver.options[key] = time_limit
+
+    def set_node_limit(self, node_limit):
+        key = self._get_node_limit_option_name()
+        self._pyomo_solver.options[key] = node_limit
+
+    def set_gap_tolerance(self, gap_tolerance):
+        key = self._get_gap_tolerance_option_name()
+        self._pyomo_solver.options[key] = gap_tolerance
+
+    @abstractmethod
+    def _get_warm_start_regexp(self):
+        pass
+
+    @abstractmethod
+    def _get_node_count_regexp(self):
+        pass
+
+    @abstractmethod
+    def _get_threads_option_name(self):
+        pass
+
+    @abstractmethod
+    def _get_time_limit_option_name(self):
+        pass
+
+    @abstractmethod
+    def _get_node_limit_option_name(self):
+        pass
+
+    @abstractmethod
+    def _get_gap_tolerance_option_name(self):
+        pass
+
+
+