Repository: qiskit-community/qiskit-optimization (language: python)

Issue:
[DNM] Add ScipyMilpOptimizer
<!--
⚠️ If you do not respect this template, your pull request will be closed.
⚠️ Your pull request title should be short detailed and understandable for all.
⚠️ Also, please add it in the CHANGELOG file under Unreleased section.
⚠️ If your pull request fixes an open issue, please link to the issue.

✅ I have added the tests to cover my changes.
✅ I have updated the documentation accordingly.
✅ I have read the CONTRIBUTING document.
-->

### Summary

Add optimizer method for linear problem based on `scipy.optimize.milp` solver. This is a [QAMP-fall-2022](https://github.com/qiskit-advocate/qamp-fall-22/issues/2) project.






Discussion / hints:

Fix patch:
diff --git a/.pylintdict b/.pylintdict
index 07a9e754b..23743a822 100644
--- a/.pylintdict
+++ b/.pylintdict
@@ -114,6 +114,7 @@ maxcut
 maxfun
 maxiter
 mdl
+milp
 minimizer
 minimumeigenoptimizer
 multiset
diff --git a/qiskit_optimization/algorithms/__init__.py b/qiskit_optimization/algorithms/__init__.py
index 445b3b6ef..c4ab6b0e8 100644
--- a/qiskit_optimization/algorithms/__init__.py
+++ b/qiskit_optimization/algorithms/__init__.py
@@ -55,6 +55,7 @@
    OptimizationResultStatus
    RecursiveMinimumEigenOptimizationResult
    RecursiveMinimumEigenOptimizer
+   ScipyMilpOptimizer
    SlsqpOptimizationResult
    SlsqpOptimizer
    SolutionSample
@@ -93,6 +94,7 @@
     RecursiveMinimumEigenOptimizationResult,
     IntermediateResult,
 )
+from .scipy_milp_optimizer import ScipyMilpOptimizer
 from .slsqp_optimizer import SlsqpOptimizer, SlsqpOptimizationResult
 from .warm_start_qaoa_optimizer import (
     BaseAggregator,
@@ -124,6 +126,7 @@
     "RecursiveMinimumEigenOptimizer",
     "RecursiveMinimumEigenOptimizationResult",
     "IntermediateResult",
+    "ScipyMilpOptimizer",
     "SlsqpOptimizer",
     "SlsqpOptimizationResult",
     "SolutionSample",
diff --git a/qiskit_optimization/algorithms/scipy_milp_optimizer.py b/qiskit_optimization/algorithms/scipy_milp_optimizer.py
new file mode 100644
index 000000000..ec86c0af8
--- /dev/null
+++ b/qiskit_optimization/algorithms/scipy_milp_optimizer.py
@@ -0,0 +1,187 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2023.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+"""The SciPy MILP optimizer wrapped to be used within Qiskit's optimization module."""
+
+from warnings import warn
+
+import numpy as np
+
+import qiskit_optimization.optionals as _optionals
+from qiskit_optimization import INFINITY
+from qiskit_optimization.algorithms.optimization_algorithm import (
+    OptimizationAlgorithm,
+    OptimizationResult,
+    OptimizationResultStatus,
+)
+from qiskit_optimization.problems.quadratic_program import (
+    ConstraintSense,
+    QuadraticProgram,
+    VarType,
+)
+
+
+def _conv_inf(val):
+    # Note: qiskit-optimization treats INFINITY as infinity
+    # while scipy treats np.inf as infinity
+    if val <= -INFINITY:
+        return -np.inf
+    elif val >= INFINITY:
+        return np.inf
+    else:
+        return val
+
+
+@_optionals.HAS_SCIPY_MILP.require_in_instance
+class ScipyMilpOptimizer(OptimizationAlgorithm):
+    """The MILP optimizer from Scipy wrapped as a Qiskit :class:`OptimizationAlgorithm`.
+
+    This class provides a wrapper for ``scipy.milp``
+    (https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.milp.html)
+    to be used within the optimization module.
+    """
+
+    def __init__(self, disp: bool = False) -> None:
+        """Initializes the ScipyMILPOptimizer.
+
+        Args:
+            disp: Whether to print MILP output or not.
+        """
+        self._disp = disp
+
+    @property
+    def disp(self) -> bool:
+        """Returns the display setting.
+
+        Returns:
+            Whether to print scipy.milp information or not.
+        """
+        return self._disp
+
+    @disp.setter
+    def disp(self, disp: bool):
+        """Set the display setting.
+
+        Args:
+            disp: The display setting.
+        """
+        self._disp = disp
+
+    # pylint:disable=unused-argument
+    def get_compatibility_msg(self, problem: QuadraticProgram) -> str:
+        """Checks whether a given problem can be solved with this optimizer.
+
+        Checks if the problem has only linear objective function and linear constraints.
+        The ``scipy.milp`` supports only linear objective function and linear constraints.
+
+        Args:
+            problem: The optimization problem to check compatibility.
+
+        Returns:
+            An empty string (if compatible) or a string describing the incompatibility.
+        """
+        msg = []
+        if problem.objective.quadratic.to_dict():
+            msg.append("scipy.milp supports only linear objective function")
+        if problem.quadratic_constraints:
+            msg.append("scipy.milp supports only linear constraints")
+        return "; ".join(msg)
+
+    def _generate_problem(self, problem: QuadraticProgram):
+        from scipy.optimize import Bounds, LinearConstraint
+        from scipy.sparse import lil_array  # pylint: disable=no-name-in-module
+
+        sense = problem.objective.sense.value
+        objective = problem.objective.linear.to_array() * sense
+
+        integrality = []
+        for variable in problem.variables:
+            if variable.vartype == VarType.CONTINUOUS:
+                integrality.append(0)
+            else:
+                integrality.append(1)
+
+        lower_bounds = [_conv_inf(variable.lowerbound) for variable in problem.variables]
+        upper_bounds = [_conv_inf(variable.upperbound) for variable in problem.variables]
+        bounds = Bounds(lb=lower_bounds, ub=upper_bounds)
+
+        lower_bounds = []
+        upper_bounds = []
+        mat = lil_array((problem.get_num_linear_constraints(), problem.get_num_vars()))
+        for i, constraint in enumerate(problem.linear_constraints):
+            for variable_id, val in constraint.linear.to_dict().items():
+                mat[i, variable_id] = val
+
+            rhs_val = _conv_inf(constraint.rhs)
+            if constraint.sense == ConstraintSense.GE:
+                lower_bounds.append(rhs_val)
+                upper_bounds.append(np.inf)
+            elif constraint.sense == ConstraintSense.LE:
+                lower_bounds.append(-np.inf)
+                upper_bounds.append(rhs_val)
+            else:
+                # ConstraintSense.EQ
+                lower_bounds.append(rhs_val)
+                upper_bounds.append(rhs_val)
+
+        constraints = LinearConstraint(mat, lower_bounds, upper_bounds)
+
+        return objective, integrality, bounds, constraints
+
+    def solve(self, problem: QuadraticProgram) -> OptimizationResult:
+        """Tries to solve the given problem using the optimizer.
+
+        Runs the optimizer to try to solve the optimization problem. If problem is not convex,
+        this optimizer may raise an exception due to incompatibility, depending on the settings.
+
+        Args:
+            problem: The problem to be solved.
+
+        Returns:
+            The result of the optimizer applied to the problem.
+
+        Raises:
+            QiskitOptimizationError: If the problem is incompatible with the optimizer.
+        """
+        from scipy.optimize import milp  # pylint: disable=no-name-in-module
+
+        self._verify_compatibility(problem)
+
+        objective, integrality, bounds, constraints = self._generate_problem(problem)
+        raw_result = milp(
+            c=objective,
+            integrality=integrality,
+            bounds=bounds,
+            constraints=constraints,
+            options={"disp": self._disp},
+        )
+
+        if raw_result.x is None:
+            warn("scipy.milp cannot solve the model. See `raw_results` for details.")
+            x = [0.0] * problem.get_num_vars()
+            status = OptimizationResultStatus.FAILURE
+        else:
+            x = []
+            for i, ele in enumerate(raw_result.x):
+                if integrality[i]:
+                    x.append(round(ele))
+                else:
+                    x.append(ele)
+            status = self._get_feasibility_status(problem, x)
+
+        return OptimizationResult(
+            x=x,
+            fval=problem.objective.evaluate(x),
+            variables=problem.variables,
+            status=status,
+            raw_results=raw_result,
+        )
diff --git a/qiskit_optimization/optionals.py b/qiskit_optimization/optionals.py
index 4293f5144..093e13ac8 100644
--- a/qiskit_optimization/optionals.py
+++ b/qiskit_optimization/optionals.py
@@ -1,6 +1,6 @@
 # This code is part of Qiskit.
 #
-# (C) Copyright IBM 2022.
+# (C) Copyright IBM 2022, 2023.
 #
 # This code is licensed under the Apache License, Version 2.0. You may
 # obtain a copy of this license in the LICENSE.txt file in the root directory
@@ -47,3 +47,11 @@
     name="Matplotlib",
     install="pip install 'qiskit-optimization[matplotlib]'",
 )
+
+HAS_SCIPY_MILP = LazyImportTester(
+    {
+        "scipy.optimize": ("milp",),
+    },
+    name="SciPy MILP",
+    install="pip install scipy>=1.9.0'",
+)
diff --git a/releasenotes/notes/add-scipy-milp-7d78a5ae8c893cc5.yaml b/releasenotes/notes/add-scipy-milp-7d78a5ae8c893cc5.yaml
new file mode 100644
index 000000000..223a73e71
--- /dev/null
+++ b/releasenotes/notes/add-scipy-milp-7d78a5ae8c893cc5.yaml
@@ -0,0 +1,9 @@
+---
+features:
+  - |
+    Adds :class:`~qiskit_optimization.algorithms.ScipyMilpOptimizer`
+    as a classical solver based on
+    `scipy.optimize.milp <https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.milp.html>`_,
+    which can solve mixed-integer linear program (MILP) problems.
+    Note that ``scipy.optimize.milp`` has been introduced by
+    `SciPy 1.9.0 <https://docs.scipy.org/doc/scipy/release.1.9.0.html>`_.


Test patch:
diff --git a/test/algorithms/test_scipy_milp_optimizer.py b/test/algorithms/test_scipy_milp_optimizer.py
new file mode 100644
index 000000000..babea8d2f
--- /dev/null
+++ b/test/algorithms/test_scipy_milp_optimizer.py
@@ -0,0 +1,270 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2023.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+""" Test ScipyMilP Optimizer """
+
+import unittest
+from test.optimization_test_case import QiskitOptimizationTestCase
+
+import numpy as np
+from ddt import ddt
+
+import qiskit_optimization.optionals as _optionals
+from qiskit_optimization import INFINITY
+from qiskit_optimization.algorithms import OptimizationResultStatus, ScipyMilpOptimizer
+from qiskit_optimization.exceptions import QiskitOptimizationError
+from qiskit_optimization.problems import QuadraticProgram
+
+
+@ddt
+class TestScipyMilpOptimizer(QiskitOptimizationTestCase):
+    """ScipyMilp Optimizer Tests."""
+
+    @unittest.skipIf(not _optionals.HAS_SCIPY_MILP, "Scipy MILP solver not available.")
+    def test_scipy_milp_optimizer(self):
+        """ScipyMilp Optimizer Test"""
+        problem = QuadraticProgram()
+        problem.continuous_var(upperbound=4)
+        problem.continuous_var(upperbound=4)
+        problem.linear_constraint(linear=[1, 1], sense="=", rhs=2)
+        problem.minimize(linear=[2, -2])
+
+        optimizer = ScipyMilpOptimizer()
+        result = optimizer.solve(problem)
+
+        self.assertAlmostEqual(result.fval, -4)
+
+    @unittest.skipIf(not _optionals.HAS_SCIPY_MILP, "Scipy MILP solver not available.")
+    def test_disp(self):
+        """Test disp"""
+        with self.subTest("default"):
+            optimizer = ScipyMilpOptimizer()
+            self.assertFalse(optimizer.disp)
+        with self.subTest("init True"):
+            optimizer = ScipyMilpOptimizer(disp=True)
+            self.assertTrue(optimizer.disp)
+        with self.subTest("setter / True"):
+            optimizer = ScipyMilpOptimizer(disp=False)
+            optimizer.disp = True
+            self.assertTrue(optimizer.disp)
+        with self.subTest("setter / False"):
+            optimizer = ScipyMilpOptimizer(disp=True)
+            optimizer.disp = False
+            self.assertFalse(optimizer.disp)
+
+    @unittest.skipIf(not _optionals.HAS_SCIPY_MILP, "Scipy MILP solver not available.")
+    def test_compatibility(self):
+        """Test compatibility"""
+        with self.subTest("quadratic obj min"):
+            problem = QuadraticProgram()
+            _ = problem.binary_var("x")
+            problem.minimize(quadratic={(0, 0): 1})
+            optimizer = ScipyMilpOptimizer()
+            msg = optimizer.get_compatibility_msg(problem)
+            self.assertEqual(msg, "scipy.milp supports only linear objective function")
+        with self.subTest("quadratic obj max"):
+            problem = QuadraticProgram()
+            _ = problem.binary_var("x")
+            problem.maximize(quadratic={(0, 0): 1})
+            optimizer = ScipyMilpOptimizer()
+            msg = optimizer.get_compatibility_msg(problem)
+            self.assertEqual(msg, "scipy.milp supports only linear objective function")
+        with self.subTest("quadratic constraint"):
+            problem = QuadraticProgram()
+            _ = problem.binary_var("x")
+            problem.quadratic_constraint(quadratic={(0, 0): 1})
+            optimizer = ScipyMilpOptimizer()
+            msg = optimizer.get_compatibility_msg(problem)
+            self.assertEqual(msg, "scipy.milp supports only linear constraints")
+        with self.subTest("quadratic obj / quadratic constraint"):
+            problem = QuadraticProgram()
+            _ = problem.binary_var("x")
+            problem.minimize(quadratic={(0, 0): 1})
+            problem.quadratic_constraint(quadratic={(0, 0): 1})
+            optimizer = ScipyMilpOptimizer()
+            msg = optimizer.get_compatibility_msg(problem)
+            self.assertEqual(
+                msg,
+                "scipy.milp supports only linear objective function; "
+                "scipy.milp supports only linear constraints",
+            )
+
+    @unittest.skipIf(not _optionals.HAS_SCIPY_MILP, "Scipy MILP solver not available.")
+    def test_error(self):
+        """Test errors"""
+        with self.subTest("quadratic obj min"):
+            problem = QuadraticProgram()
+            _ = problem.binary_var("x")
+            problem.minimize(quadratic={(0, 0): 1})
+            optimizer = ScipyMilpOptimizer()
+            with self.assertRaises(QiskitOptimizationError):
+                _ = optimizer.solve(problem)
+        with self.subTest("quadratic obj max"):
+            problem = QuadraticProgram()
+            _ = problem.binary_var("x")
+            problem.maximize(quadratic={(0, 0): 1})
+            optimizer = ScipyMilpOptimizer()
+            with self.assertRaises(QiskitOptimizationError):
+                _ = optimizer.solve(problem)
+        with self.subTest("quadratic constraint"):
+            problem = QuadraticProgram()
+            _ = problem.binary_var("x")
+            problem.quadratic_constraint(quadratic={(0, 0): 1})
+            optimizer = ScipyMilpOptimizer()
+            with self.assertRaises(QiskitOptimizationError):
+                _ = optimizer.solve(problem)
+
+    @unittest.skipIf(not _optionals.HAS_SCIPY_MILP, "Scipy MILP solver not available.")
+    def test_senses(self):
+        """Test constraint senses"""
+        with self.subTest("min / <="):
+            problem = QuadraticProgram()
+            _ = problem.continuous_var(0, 1, "x")
+            _ = problem.binary_var("y")
+            problem.minimize(linear={"x": 1, "y": -1})
+            problem.linear_constraint({"x": 2, "y": 1}, "<=", 1)
+            optimizer = ScipyMilpOptimizer()
+            result = optimizer.solve(problem)
+            self.assertAlmostEqual(result.fval, -1)
+            np.testing.assert_allclose(result.x, [0, 1])
+            self.assertEqual(result.status, OptimizationResultStatus.SUCCESS)
+
+        with self.subTest("max / <="):
+            problem = QuadraticProgram()
+            _ = problem.continuous_var(0, 1, "x")
+            _ = problem.binary_var("y")
+            problem.maximize(linear={"x": 1, "y": -1})
+            problem.linear_constraint({"x": 2, "y": 1}, "<=", 1)
+            optimizer = ScipyMilpOptimizer()
+            result = optimizer.solve(problem)
+            self.assertAlmostEqual(result.fval, 0.5)
+            np.testing.assert_allclose(result.x, [0.5, 0])
+            self.assertEqual(result.status, OptimizationResultStatus.SUCCESS)
+
+        with self.subTest("min / >="):
+            problem = QuadraticProgram()
+            _ = problem.continuous_var(0, 1, "x")
+            _ = problem.binary_var("y")
+            problem.minimize(linear={"x": 1, "y": -1})
+            problem.linear_constraint({"x": 2, "y": 1}, ">=", 1)
+            optimizer = ScipyMilpOptimizer()
+            result = optimizer.solve(problem)
+            self.assertAlmostEqual(result.fval, -1.0)
+            np.testing.assert_allclose(result.x, [0, 1.0])
+            self.assertEqual(result.status, OptimizationResultStatus.SUCCESS)
+
+        with self.subTest("max / >="):
+            problem = QuadraticProgram()
+            _ = problem.continuous_var(0, 1, "x")
+            _ = problem.binary_var("y")
+            problem.maximize(linear={"x": 1, "y": -1})
+            problem.linear_constraint({"x": 2, "y": 1}, ">=", 1)
+            optimizer = ScipyMilpOptimizer()
+            result = optimizer.solve(problem)
+            self.assertAlmostEqual(result.fval, 1)
+            np.testing.assert_allclose(result.x, [1, 0])
+            self.assertEqual(result.status, OptimizationResultStatus.SUCCESS)
+
+        with self.subTest("min / =="):
+            problem = QuadraticProgram()
+            _ = problem.continuous_var(0, 1, "x")
+            _ = problem.binary_var("y")
+            problem.minimize(linear={"x": 1, "y": -1})
+            problem.linear_constraint({"x": 2, "y": 1}, "==", 2)
+            optimizer = ScipyMilpOptimizer()
+            result = optimizer.solve(problem)
+            self.assertAlmostEqual(result.fval, -0.5)
+            np.testing.assert_allclose(result.x, [0.5, 1.0])
+            self.assertEqual(result.status, OptimizationResultStatus.SUCCESS)
+
+        with self.subTest("max / =="):
+            problem = QuadraticProgram()
+            _ = problem.continuous_var(0, 1, "x")
+            _ = problem.binary_var("y")
+            problem.maximize(linear={"x": 1, "y": -1})
+            problem.linear_constraint({"x": 2, "y": 1}, "==", 2)
+            optimizer = ScipyMilpOptimizer()
+            result = optimizer.solve(problem)
+            self.assertAlmostEqual(result.fval, 1)
+            np.testing.assert_allclose(result.x, [1, 0])
+            self.assertEqual(result.status, OptimizationResultStatus.SUCCESS)
+
+    @unittest.skipIf(not _optionals.HAS_SCIPY_MILP, "Scipy MILP solver not available.")
+    def test_infinity(self):
+        """Test infinity"""
+        with self.subTest("ub = infinity default"):
+            problem = QuadraticProgram()
+            _ = problem.continuous_var(name="x")
+            problem.maximize(linear={0: 1})
+            optimizer = ScipyMilpOptimizer()
+            with self.assertWarns(UserWarning):
+                result = optimizer.solve(problem)
+            self.assertEqual(result.status, OptimizationResultStatus.FAILURE)
+
+        with self.subTest("ub = infinity manual"):
+            problem = QuadraticProgram()
+            _ = problem.continuous_var(upperbound=INFINITY, name="x")
+            problem.maximize(linear={0: 1})
+            optimizer = ScipyMilpOptimizer()
+            with self.assertWarns(UserWarning):
+                result = optimizer.solve(problem)
+            self.assertEqual(result.status, OptimizationResultStatus.FAILURE)
+
+        with self.subTest("ub > infinity"):
+            problem = QuadraticProgram()
+            _ = problem.continuous_var(upperbound=1e100, name="x")
+            problem.maximize(linear={0: 1})
+            optimizer = ScipyMilpOptimizer()
+            with self.assertWarns(UserWarning):
+                result = optimizer.solve(problem)
+            self.assertEqual(result.status, OptimizationResultStatus.FAILURE)
+
+        with self.subTest("ub < infinity"):
+            problem = QuadraticProgram()
+            _ = problem.continuous_var(upperbound=1e10, name="x")
+            problem.maximize(linear={0: 1})
+            optimizer = ScipyMilpOptimizer()
+            result = optimizer.solve(problem)
+            np.testing.assert_allclose(result.x, [1e10])
+            self.assertEqual(result.status, OptimizationResultStatus.SUCCESS)
+
+        with self.subTest("lb = -infinity manual"):
+            problem = QuadraticProgram()
+            _ = problem.continuous_var(lowerbound=-INFINITY, name="x")
+            problem.minimize(linear={0: 1})
+            optimizer = ScipyMilpOptimizer()
+            with self.assertWarns(UserWarning):
+                result = optimizer.solve(problem)
+            self.assertEqual(result.status, OptimizationResultStatus.FAILURE)
+
+        with self.subTest("lb < -infinity"):
+            problem = QuadraticProgram()
+            _ = problem.continuous_var(lowerbound=-1e100, name="x")
+            problem.minimize(linear={0: 1})
+            optimizer = ScipyMilpOptimizer()
+            with self.assertWarns(UserWarning):
+                result = optimizer.solve(problem)
+                print(result)
+            self.assertEqual(result.status, OptimizationResultStatus.FAILURE)
+
+        with self.subTest("lb > -infinity"):
+            problem = QuadraticProgram()
+            _ = problem.continuous_var(lowerbound=-1e10, name="x")
+            problem.minimize(linear={0: 1})
+            optimizer = ScipyMilpOptimizer()
+            result = optimizer.solve(problem)
+            np.testing.assert_allclose(result.x, [-1e10])
+            self.assertEqual(result.status, OptimizationResultStatus.SUCCESS)
+
+
+if __name__ == "__main__":
+    unittest.main()