Use bounds on delta rather than additional constraints

pyproject.toml

@@ -10,7 +10,7 @@ readme = "README.md"
 license = { file = "LICENSE" }
 authors = [{ name = "Timothy Nunn", email = "timothy.nunn@ukaea.uk" }]
 requires-python = ">=3.10"
-dependencies = ["numpy>=1.24", "cvxpy>=1.5.2"]
+dependencies = ["numpy>=1.24", "cvxpy>=1.6"]
 classifiers = [
     "Intended Audience :: Science/Research",
     "License :: OSI Approved :: MIT License",

src/pyvmcon/problem.py

@@ -93,20 +93,20 @@ def __init__(
         self,
         f: _ScalarReturnFunctionAlias,
         df: _VectorReturnFunctionAlias,
-        equality_constraints: list[_ScalarReturnFunctionAlias],
-        inequality_constraints: list[_ScalarReturnFunctionAlias],
-        dequality_constraints: list[_VectorReturnFunctionAlias],
-        dinequality_constraints: list[_VectorReturnFunctionAlias],
+        equality_constraints: list[_ScalarReturnFunctionAlias] | None = None,
+        inequality_constraints: list[_ScalarReturnFunctionAlias] | None = None,
+        dequality_constraints: list[_VectorReturnFunctionAlias] | None = None,
+        dinequality_constraints: list[_VectorReturnFunctionAlias] | None = None,
     ) -> None:
         """Construct the problem."""
         super().__init__()
 
         self._f = f
         self._df = df
-        self._equality_constraints = equality_constraints
-        self._inequality_constraints = inequality_constraints
-        self._dequality_constraints = dequality_constraints
-        self._dinequality_constraints = dinequality_constraints
+        self._equality_constraints = equality_constraints or []
+        self._inequality_constraints = inequality_constraints or []
+        self._dequality_constraints = dequality_constraints or []
+        self._dinequality_constraints = dinequality_constraints or []
 
     def __call__(self, x: VectorType) -> Result:
         """Evaluate the problem at input point x."""

src/pyvmcon/vmcon.py

@@ -306,23 +306,20 @@ def solve_qsp(
     different `solver` in the `options` dictionary.
 
     """
-    delta = cp.Variable(x.shape)
+    delta = cp.Variable(
+        x.shape,
+        bounds=[
+            lbs - x if lbs is not None else None,
+            ubs - x if ubs is not None else None,
+        ],
+    )
     problem_statement = cp.Minimize(
         result.f + (0.5 * cp.quad_form(delta, B)) + (delta.T @ result.df),
     )
 
-    equality_index = 0
-
     constraints = []
     if problem.has_inequality:
-        equality_index += 1
         constraints.append((result.die @ delta) + result.ie >= 0)
-    if lbs is not None:
-        equality_index += 1
-        constraints.append(x + delta >= lbs)
-    if ubs is not None:
-        equality_index += 1
-        constraints.append(x + delta <= ubs)
     if problem.has_equality:
         constraints.append((result.deq @ delta) + result.eq == 0)
 
@@ -338,13 +335,13 @@ def solve_qsp(
 
     if problem.has_inequality and problem.has_equality:
         lamda_inequality = qsp.constraints[0].dual_value
-        lamda_equality = -qsp.constraints[equality_index].dual_value
+        lamda_equality = -qsp.constraints[1].dual_value
 
     elif problem.has_inequality and not problem.has_equality:
         lamda_inequality = qsp.constraints[0].dual_value
 
     elif not problem.has_inequality and problem.has_equality:
-        lamda_equality = -qsp.constraints[equality_index].dual_value
+        lamda_equality = -qsp.constraints[0].dual_value
 
     return delta.value, lamda_equality, lamda_inequality
 

tests/test_vmcon_bounds.py

@@ -0,0 +1,51 @@
+"""Test VMCON's use of bounds in simple linear functions."""
+
+import numpy as np
+import pytest
+
+from pyvmcon.problem import Problem
+from pyvmcon.vmcon import solve
+
+
+@pytest.mark.parametrize(
+    ("problem", "expected"),
+    [
+        (Problem(f=lambda x: x[0], df=lambda _: np.array([1])), -10.0),
+        (Problem(f=lambda x: -x[0], df=lambda _: np.array([-1])), 10.0),
+    ],
+)
+def test_vmcon_1d_10bounds(problem, expected):
+    x, _, _, _ = solve(
+        problem=problem,
+        x=np.array([0.0]),
+        lbs=np.array([-10]),
+        ubs=np.array([10]),
+        max_iter=100,
+    )
+
+    assert x.item() == expected
+
+
+@pytest.mark.parametrize(
+    ("problem", "expected"),
+    [
+        (
+            Problem(f=lambda x: x[0] + x[1], df=lambda _: np.array([1, 1])),
+            [-10.0, -20.0],
+        ),
+        (
+            Problem(f=lambda x: -x[0] - x[1], df=lambda _: np.array([-1, -1])),
+            [20.0, 10.0],
+        ),
+    ],
+)
+def test_vmcon_2d_1020bounds(problem, expected):
+    x, _, _, _ = solve(
+        problem=problem,
+        x=np.array([0.0, 0.0]),
+        lbs=np.array([-10, -20]),
+        ubs=np.array([20, 10]),
+        max_iter=100,
+    )
+
+    assert (x == expected).all()

tests/test_vmcon_paper.py

@@ -140,6 +140,12 @@ def test_vmcon_paper_feasible_examples(vmcon_example: VMCONTestAsset):
     assert lamda_equality == pytest.approx(vmcon_example.expected_lamda_equality)
     assert lamda_inequality == pytest.approx(vmcon_example.expected_lamda_inequality)
 
+    if vmcon_example.lbs is not None:
+        assert (x >= vmcon_example.lbs).all()
+
+    if vmcon_example.ubs is not None:
+        assert (x <= vmcon_example.ubs).all()
+
 
 @pytest.mark.parametrize(
     "vmcon_example",