diff --git a/pybobyqa/solver.py b/pybobyqa/solver.py
index 3097131..fa76eb2 100644
--- a/pybobyqa/solver.py
+++ b/pybobyqa/solver.py
@@ -768,6 +768,10 @@ def solve(objfun, x0, args=(), bounds=None, projections=None, npt=None, rhobeg=N
     if exit_info is None and projections is not None and type(projections) != list:
         exit_info = ExitInformation(EXIT_INPUT_ERROR, "projections must be a list of functions")
 
+    if projections is not None and len(projections) == 0:
+        # empty list given
+        projections = None
+
     if maxfun <= npt:
         warnings.warn("maxfun <= npt: Are you sure your budget is large enough?", RuntimeWarning)
 
diff --git a/pybobyqa/tests/test_model.py b/pybobyqa/tests/test_model.py
index 1eaa988..0989ea3 100644
--- a/pybobyqa/tests/test_model.py
+++ b/pybobyqa/tests/test_model.py
@@ -55,7 +55,7 @@ def runTest(self):
         x0 = np.array([-1.2, 1.0])
         xl = -1e20 * np.ones((n,))
         xu = 1e20 * np.ones((n,))
-        model = Model(npt, x0, rosenbrock(x0), xl, xu, 1)
+        model = Model(npt, x0, rosenbrock(x0), xl, xu, [], 1)
         self.assertEqual(model.npt(), npt, 'Wrong npt after initialisation')
         self.assertTrue(array_compare(model.xopt(abs_coordinates=True), x0), 'Wrong xopt after initialisation')
         self.assertTrue(array_compare(model.fopt(), rosenbrock(x0)), 'Wrong fopt after initialisation')
@@ -103,7 +103,7 @@ def runTest(self):
         x0 = np.array([-1.2, 1.0])
         xl = -1e20 * np.ones((n,))
         xu = 1e20 * np.ones((n,))
-        model = Model(npt, x0, rosenbrock(x0), xl, xu, 1)
+        model = Model(npt, x0, rosenbrock(x0), xl, xu, [], 1)
         # Now add better point
         x1 = np.array([1.0, 0.9])
         f1 = rosenbrock(x1)
@@ -131,7 +131,7 @@ def runTest(self):
         x0 = np.array([-1.2, 1.0])
         xl = -1e2 * np.ones((n,))
         xu = 1e2 * np.ones((n,))
-        model = Model(npt, x0, rosenbrock(x0), xl, xu, 1)
+        model = Model(npt, x0, rosenbrock(x0), xl, xu, [], 1)
         self.assertTrue(array_compare(model.sl, xl - x0), 'Wrong sl after initialisation')
         self.assertTrue(array_compare(model.su, xu - x0), 'Wrong su after initialisation')
         x1 = np.array([1.0, 0.9])
@@ -204,7 +204,7 @@ def runTest(self):
         x0 = np.array([-1.2, 1.0])
         xl = -1e2 * np.ones((n,))
         xu = 1e2 * np.ones((n,))
-        model = Model(npt, x0, rosenbrock(x0), xl, xu, 1)
+        model = Model(npt, x0, rosenbrock(x0), xl, xu, [], 1)
         x1 = np.array([1.0, 0.9])
         model.change_point(1, x1 - model.xbase, rosenbrock(x1))
         x2 = np.array([1.0, 1.0])
@@ -224,13 +224,13 @@ def runTest(self):
         x0 = np.array([-1.2, 1.0])
         xl = -1e2 * np.ones((n,))
         xu = 1e2 * np.ones((n,))
-        model = Model(npt, x0, rosenbrock(x0), xl, xu, 1)
+        model = Model(npt, x0, rosenbrock(x0), xl, xu, [], 1)
         x1 = np.array([1.0, 0.9])
         model.change_point(1, x1 - model.xbase, rosenbrock(x1))
         x2 = np.array([2.0, 0.9])
         model.change_point(2, x2 - model.xbase, rosenbrock(x2))
         self.assertAlmostEqual(model.min_objective_value(), -1e20, msg='Wrong min obj value')
-        model = Model(npt, x0, rosenbrock(x0), xl, xu, 1, abs_tol=1.0)
+        model = Model(npt, x0, rosenbrock(x0), xl, xu, [], 1, abs_tol=1.0)
         self.assertAlmostEqual(model.min_objective_value(), 1.0, msg='Wrong min obj value 3')
 
 
@@ -241,7 +241,7 @@ def runTest(self):
         x0 = np.array([-1.2, 1.0])
         xl = -1e2 * np.ones((n,))
         xu = 1e2 * np.ones((n,))
-        model = Model(npt, x0, rosenbrock(x0), xl, xu, 1, precondition=False)
+        model = Model(npt, x0, rosenbrock(x0), xl, xu, [], 1, precondition=False)
         x1 = np.array([1.0, 0.9])
         model.change_point(1, x1 - model.xbase, rosenbrock(x1))
         x2 = np.array([2.0, 0.9])
@@ -280,7 +280,7 @@ def runTest(self):
         x0 = np.array([1.0, 1.0])
         xl = -1e2 * np.ones((n,))
         xu = 1e2 * np.ones((n,))
-        model = Model(npt, x0, objfun(x0), xl, xu, 1, precondition=False)
+        model = Model(npt, x0, objfun(x0), xl, xu, [], 1, precondition=False)
         x1 = x0 + np.array([1.0, 0.0])
         model.change_point(1, x1 - model.xbase, objfun(x1))
         x2 = x0 + np.array([0.1, 0.9])
@@ -323,7 +323,7 @@ def runTest(self):
         self.assertTrue(np.allclose(hess, model.model_hess), 'Bad Hessian')
 
         # Build a new model
-        model2 = Model(npt, x0, objfun(x0), xl, xu, 1, precondition=False)
+        model2 = Model(npt, x0, objfun(x0), xl, xu, [], 1, precondition=False)
         model2.change_point(1, x1 - model.xbase, objfun(x1))
         model2.change_point(2, x2 - model.xbase, objfun(x2))
         model2.change_point(3, x3 - model.xbase, objfun(x3))
@@ -350,7 +350,7 @@ def runTest(self):
         # print(model2.model_hess)
 
         # Build a new model
-        model3 = Model(npt, x0, objfun(x0), xl, xu, 1, precondition=False)
+        model3 = Model(npt, x0, objfun(x0), xl, xu, [], 1, precondition=False)
         model3.change_point(1, x1 - model.xbase, objfun(x1))
         model3.change_point(2, x2 - model.xbase, objfun(x2))
         model3.change_point(3, x3 - model.xbase, objfun(x3))
@@ -376,7 +376,7 @@ def runTest(self):
         x0 = np.array([1.0, 1.0])
         xl = -1e2 * np.ones((n,))
         xu = 1e2 * np.ones((n,))
-        model = Model(npt, x0, objfun(x0), xl, xu, 1, precondition=False)
+        model = Model(npt, x0, objfun(x0), xl, xu, [], 1, precondition=False)
         x1 = x0 + np.array([1.0, 0.0])
         model.change_point(1, x1 - model.xbase, objfun(x1))
         x2 = x0 + np.array([0.1, 0.9])
@@ -420,7 +420,7 @@ def runTest(self):
         self.assertTrue(np.allclose(hess, model.model_hess), 'Bad Hessian')
 
         # Build a new model
-        model2 = Model(npt, x0, objfun(x0), xl, xu, 1, precondition=False)
+        model2 = Model(npt, x0, objfun(x0), xl, xu, [], 1, precondition=False)
         model2.change_point(1, x1 - model.xbase, objfun(x1))
         model2.change_point(2, x2 - model.xbase, objfun(x2))
         model2.change_point(3, x3 - model.xbase, objfun(x3))
@@ -448,7 +448,7 @@ def runTest(self):
         # print(model2.model_hess)
 
         # Build a new model
-        model3 = Model(npt, x0, objfun(x0), xl, xu, 1, precondition=False)
+        model3 = Model(npt, x0, objfun(x0), xl, xu, [], 1, precondition=False)
         model3.change_point(1, x1 - model.xbase, objfun(x1))
         model3.change_point(2, x2 - model.xbase, objfun(x2))
         model3.change_point(3, x3 - model.xbase, objfun(x3))
@@ -475,7 +475,7 @@ def runTest(self):
         x0 = np.array([1.0, 1.0])
         xl = -1e2 * np.ones((n,))
         xu = 1e2 * np.ones((n,))
-        model = Model(npt, x0, objfun(x0), xl, xu, 1)
+        model = Model(npt, x0, objfun(x0), xl, xu, [], 1)
         x1 = x0 + np.array([1.0, 0.0])
         model.change_point(1, x1 - model.xbase, objfun(x1))
         x2 = x0 + np.array([0.1, 0.9])
@@ -511,7 +511,7 @@ def runTest(self):
         x0 = np.array([-1.2, 1.0])
         xl = -1e2 * np.ones((n,))
         xu = 1e2 * np.ones((n,))
-        model = Model(npt, x0, rosenbrock(x0), xl, xu, 1)
+        model = Model(npt, x0, rosenbrock(x0), xl, xu, [], 1)
         x1 = np.array([1.0, 0.9])
         model.change_point(1, x1 - model.xbase, rosenbrock(x1))
         x2 = np.array([2.0, 0.9])
@@ -534,7 +534,7 @@ def runTest(self):
         x0 = np.array([1.0, 1.0])
         xl = -1e2 * np.ones((n,))
         xu = 1e2 * np.ones((n,))
-        model = Model(npt, x0, objfun(x0), xl, xu, 1)
+        model = Model(npt, x0, objfun(x0), xl, xu, [], 1)
         x1 = x0 + np.array([1.0, 0.0])
         model.change_point(1, x1 - model.xbase, objfun(x1))
         x2 = x0 + np.array([0.1, 0.9])
@@ -559,7 +559,7 @@ def runTest(self):
         x0 = np.array([1.0, 1.0])
         xl = -1e2 * np.ones((n,))
         xu = 1e2 * np.ones((n,))
-        model = Model(npt, x0, objfun(x0), xl, xu, 1)
+        model = Model(npt, x0, objfun(x0), xl, xu, [], 1)
         x1 = x0 + np.array([1.0, 0.0])
         model.change_point(1, x1 - model.xbase, objfun(x1))
         x2 = x0 + np.array([0.1, 0.9])
@@ -586,7 +586,7 @@ def runTest(self):
         x0 = np.array([1.0, 1.0])
         xl = -1e2 * np.ones((n,))
         xu = 1e2 * np.ones((n,))
-        model = Model(npt, x0, objfun(x0), xl, xu, 1)
+        model = Model(npt, x0, objfun(x0), xl, xu, [], 1)
         x1 = x0 + np.array([1.0, 0.0])
         model.change_point(1, x1 - model.xbase, objfun(x1))
         x2 = x0 + np.array([0.1, 0.9])
@@ -616,7 +616,7 @@ def runTest(self):
         delta = 0.5
         xl = -1e2 * np.ones((n,))
         xu = 1e2 * np.ones((n,))
-        model = Model(npt, x0, rosenbrock(x0), xl, xu, 1)
+        model = Model(npt, x0, rosenbrock(x0), xl, xu, [], 1)
         model.add_new_sample(0, rosenbrock(x0))
         x1 = x0 + delta * np.array([1.0, 0.0])
         model.change_point(1, x1 - model.xbase, rosenbrock(x1))
@@ -636,7 +636,7 @@ def runTest(self):
         x0 = np.array([0.5, 0.5])
         xl = -1e2 * np.ones((n,))
         xu = 1e2 * np.ones((n,))
-        model = Model(npt, x0, objfun(x0), xl, xu, 1)
+        model = Model(npt, x0, objfun(x0), xl, xu, [], 1)
         x1 = np.array([0.05, 0.1])
         model.change_point(1, x1 - model.xbase, objfun(x1))
         x2 = np.array([0.1, 0.05])
@@ -660,7 +660,7 @@ def runTest(self):
         x0 = np.array([0.5, 0.5])
         xl = -1e2 * np.ones((n,))
         xu = 1e2 * np.ones((n,))
-        model = Model(npt, x0, objfun(x0), xl, xu, 1)
+        model = Model(npt, x0, objfun(x0), xl, xu, [], 1)
         x1 = np.array([0.524, 0.0006])
         model.change_point(1, x1 - model.xbase, objfun(x1))
         x2 = np.array([0.032, 0.323])
@@ -681,7 +681,7 @@ def runTest(self):
         x0 = np.array([-1.2, 1.0])
         xl = -1e2 * np.ones((n,))
         xu = 1e2 * np.ones((n,))
-        model = Model(npt, x0, rosenbrock(x0), xl, xu, 1)
+        model = Model(npt, x0, rosenbrock(x0), xl, xu, [], 1)
         x1 = np.array([1.0, 0.9])
         model.change_point(1, x1 - model.xbase, rosenbrock(x1))
         x2 = np.array([2.0, 0.9])
diff --git a/pybobyqa/tests/test_trust_region.py b/pybobyqa/tests/test_trust_region.py
index 282e482..b2ea10a 100644
--- a/pybobyqa/tests/test_trust_region.py
+++ b/pybobyqa/tests/test_trust_region.py
@@ -27,7 +27,7 @@
 import numpy as np
 import unittest
 
-from pybobyqa.trust_region import trsbox, trsbox_geometry
+from pybobyqa.trust_region import trsbox, trsbox_geometry, ctrsbox, ctrsbox_geometry
 from pybobyqa.util import model_value
 
 
@@ -70,6 +70,11 @@ def cauchy_pt_box(g, hess, delta, lower, upper):
         crvmin = -1.0
     return s, red, crvmin
 
+def p_box(x, l, u):
+    return np.minimum(np.maximum(x,l), u)
+
+def p_ball(x, c, r):
+    return c + (r/np.max([np.linalg.norm(x-c),r]))*(x-c)
 
 class TestUncInternal(unittest.TestCase):
     def runTest(self):
@@ -382,3 +387,198 @@ def runTest(self):
         # print(x)
         # print(xtrue)
         self.assertTrue(np.allclose(x, xtrue, atol=1e-3), "Wrong step")
+
+
+class TestGeom3CDFO(unittest.TestCase):
+    def runTest(self):
+        n = 2
+        xbase = np.zeros((n,))
+        xopt = np.ones((n,))
+        g = np.array([1.0, -1.0])
+        a = np.array([-2.0, -2.0])
+        b = np.array([1.0, 2.0])
+        H = np.zeros((n, n))
+        proj = lambda x: p_box(x, a, b)
+        delta = 5.0
+        c = 3.0  # may want to max instead
+        x = ctrsbox_geometry(xbase, xopt, c, g, H, a, b, [proj], delta)
+        xtrue = np.array([1.0, -2.0])
+        print("x = ", x)
+        print("xtrue = ", xtrue)
+        self.assertTrue(np.max(np.abs(x - xtrue)) < 1e-10, 'Wrong step')
+
+
+class TestGeomWithHessianCDFO(unittest.TestCase):
+    def runTest(self):
+        n = 3
+        c = -1.0
+        g = np.array([1.0, 0.0, 1.0])
+        H = np.array([[-2.0, 0.0, 0.0], [0.0, -1.0, 0.0], [0.0, 0.0, -1.0]])
+        Delta = 5.0 / 12.0
+        xbase = np.zeros((n,))
+        xopt = np.zeros((n,))
+        sl = -1e20 * np.ones((n,))
+        su = 1e20 * np.ones((n,))
+        proj = lambda x: p_box(x, sl, su)
+        x = ctrsbox_geometry(xbase, xopt, c, g, H, sl, su, [proj], Delta)
+        xtrue = np.array([-1.0 / 3.0, 0.0, -0.25])
+        # print(x)
+        # print(xtrue)
+        self.assertTrue(np.allclose(x, xtrue, atol=5e-2), "Wrong step")
+
+
+class TestGeomWithHessian2CDFO(unittest.TestCase):
+    def runTest(self):
+        n = 3
+        c = 1.0  # changed this value to force max rather than min
+        g = np.array([1.0, 0.0, 1.0])
+        H = np.array([[-2.0, 0.0, 0.0], [0.0, -1.0, 0.0], [0.0, 0.0, -1.0]])
+        Delta = 5.0 / 12.0
+        xbase = np.zeros((n,))
+        xopt = np.zeros((n,))
+        sl = -1e20 * np.ones((n,))
+        su = 1e20 * np.ones((n,))
+        proj = lambda x: p_box(x, sl, su)
+        x = ctrsbox_geometry(xbase, xopt, c, g, H, sl, su, [proj], Delta)
+        xtrue = np.array([0.25, 0.0, 1.0 / 3.0])  # max solution
+        # print(x)
+        # print(xtrue)
+        self.assertTrue(np.allclose(x, xtrue, atol=1e-3), "Wrong step")
+
+
+class TestUncInternalCDFO(unittest.TestCase):
+    def runTest(self):
+        n = 3
+        xbase = np.zeros((n,))
+        g = np.array([1.0, 0.0, 1.0])
+        H = np.array([[1.0, 0.0, 0.0], [0.0, 2.0, 0.0], [0.0, 0.0, 2.0]])
+        Delta = 2.0
+        xopt = np.ones((n,))  # trying nonzero (since bounds inactive)
+        sl = xopt-3*Delta * np.ones((n,))
+        su = xopt+3*Delta * np.ones((n,))
+        proj = lambda x: p_box(x, sl, su)
+        d, gnew, _crvmin = ctrsbox(xbase, xopt, g, H, sl, su, [proj], Delta)
+        true_d = np.array([-1.0, 0.0, -0.5])
+        est_min = model_value(g, H, d)
+        true_min = model_value(g, H, true_d)
+        s_cauchy, red_cauchy, _crvmin_cauchy = cauchy_pt(g, H, Delta)
+        self.assertTrue(est_min <= red_cauchy, 'Cauchy reduction not achieved')
+        self.assertTrue(np.all(gnew == g + H.dot(d)), 'Wrong gnew')
+
+
+class TestUncBdryCDFO(unittest.TestCase):
+    def runTest(self):
+        n = 3
+        xbase = np.zeros((n,))
+        g = np.array([1.0, 0.0, 1.0])
+        H = np.array([[1.0, 0.0, 0.0], [0.0, 2.0, 0.0], [0.0, 0.0, 2.0]])
+        Delta = 5.0 / 12.0
+        xopt = np.zeros((n,))
+        sl = xopt - 3 * Delta * np.ones((n,))
+        su = xopt + 3 * Delta * np.ones((n,))
+        proj = lambda x: p_box(x, sl, su)
+        d, gnew, _crvmin = ctrsbox(xbase, xopt, g, H, sl, su, [proj], Delta)
+        true_d = np.array([-1.0 / 3.0, 0.0, -0.25])
+        est_min = model_value(g, H, d)
+        true_min = model_value(g, H, true_d)
+        s_cauchy, red_cauchy, _crvmin_cauchy = cauchy_pt(g, H, Delta)
+        self.assertTrue(est_min <= red_cauchy, 'Cauchy reduction not achieved')
+        self.assertTrue(np.allclose(gnew, g + H.dot(d)), 'Wrong gnew')
+
+
+class TestUncHardCDFO(unittest.TestCase):
+    def runTest(self):
+        n = 3
+        xbase = np.zeros((n,))
+        g = np.array([0.0, 0.0, 1.0])
+        H = np.array([[-2.0, 0.0, 0.0], [0.0, -1.0, 0.0], [0.0, 0.0, -1.0]])
+        Delta = sqrt(2.0)
+        xopt = np.zeros((n,))
+        sl = xopt - 3 * Delta * np.ones((n,))
+        su = xopt + 3 * Delta * np.ones((n,))
+        proj = lambda x: p_box(x, sl, su)
+        d, gnew, _crvmin = ctrsbox(xbase, xopt, g, H, sl, su, [proj], Delta)
+        true_d = np.array([1.0, 0.0, -1.0])  # non-unique solution
+        est_min = model_value(g, H, d)
+        true_min = model_value(g, H, true_d)
+        s_cauchy, red_cauchy, _crvmin_cauchy = cauchy_pt(g, H, Delta)
+        self.assertTrue(est_min <= red_cauchy, 'Cauchy reduction not achieved')
+        self.assertTrue(np.allclose(gnew, g + H.dot(d)), 'Wrong gnew')
+
+
+class TestConInternalCDFO(unittest.TestCase):
+    def runTest(self):
+        n = 3
+        xbase = np.zeros((n,))
+        g = np.array([1.0, 0.0, 1.0])
+        H = np.array([[1.0, 0.0, 0.0], [0.0, 2.0, 0.0], [0.0, 0.0, 2.0]])
+        Delta = 2.0
+        xopt = np.ones((n,))  # trying nonzero (since bounds inactive)
+        sl = xopt + np.array([-0.5, -10.0, -10.0])
+        su = xopt + np.array([10.0, 10.0, 10.0])
+        proj = lambda x: p_box(x, sl, su)
+        d, gnew, _crvmin = ctrsbox(xbase, xopt, g, H, sl, su, [proj], Delta)
+        true_d = np.array([-1.0, 0.0, -0.5])
+        est_min = model_value(g, H, d)
+        true_min = model_value(g, H, true_d)
+        s_cauchy, red_cauchy, _crvmin_cauchy = cauchy_pt_box(g, H, Delta, sl-xopt, su-xopt)
+        self.assertTrue(est_min <= red_cauchy, 'Cauchy reduction not achieved')
+        self.assertTrue(np.all(gnew == g + H.dot(d)), 'Wrong gnew')
+
+
+class TestConBdryCDFO(unittest.TestCase):
+    def runTest(self):
+        n = 3
+        xbase = np.ones((n,))
+        g = np.array([1.0, 0.0, 1.0])
+        H = np.array([[1.0, 0.0, 0.0], [0.0, 2.0, 0.0], [0.0, 0.0, 2.0]])
+        Delta = 5.0 / 12.0
+        xopt = np.zeros((n,))
+        sl = xopt + np.array([-0.3, -0.01, -0.1])
+        su = xopt + np.array([10.0, 1.0, 10.0])
+        proj = lambda x: p_box(x, xbase + sl, xbase + su)  # xbase != 0, adjust projections accordingly
+        d, gnew, _crvmin = ctrsbox(xbase, xopt, g, H, sl, su, [proj], Delta)  # proj=duplicate of sl, su
+        true_d = np.array([-1.0 / 3.0, 0.0, -0.25])
+        est_min = model_value(g, H, d)
+        true_min = model_value(g, H, true_d)
+        s_cauchy, red_cauchy, _crvmin_cauchy = cauchy_pt_box(g, H, Delta, sl - xopt, su - xopt)
+        self.assertTrue(est_min <= red_cauchy, 'Cauchy reduction not achieved')
+        self.assertTrue(np.max(np.abs(gnew - g - H.dot(d))) < 1e-10, 'Wrong gnew')
+
+
+class TestBoxBallInternalCDFO(unittest.TestCase):
+    def runTest(self):
+        n = 3
+        xbase = np.zeros((n,))
+        g = np.array([1.0, 0.0, 1.0])
+        H = np.array([[1.0, 0.0, 0.0], [0.0, 2.0, 0.0], [0.0, 0.0, 2.0]])
+        Delta = 2.0
+        xopt = np.ones((n,))  # trying nonzero (since bounds inactive)
+        sl = xopt + np.array([-0.5, -10.0, -10.0])
+        su = xopt + np.array([10.0, 10.0, 10.0])
+        boxproj = lambda x: p_box(x, sl, su)
+        ballproj = lambda x: p_ball(x, xopt, 5)
+        d, gnew, _crvmin = ctrsbox(xbase, xopt, g, H, sl, su, [ballproj], Delta)
+        true_d = np.array([-0.5, 0.0, -0.5])
+        est_min = model_value(g, H, d)
+        true_min = model_value(g, H, true_d)
+        self.assertTrue(est_min <= true_min + 1e-3, 'Sufficient decrease not achieved')
+
+
+class TestBoxBallBdryCDFO(unittest.TestCase):
+    def runTest(self):
+        n = 3
+        g = np.array([1.0, 0.0, 1.0])
+        H = np.array([[1.0, 0.0, 0.0], [0.0, 2.0, 0.0], [0.0, 0.0, 2.0]])
+        Delta = 5.0 / 12.0
+        xbase = np.zeros((n,))
+        xopt = np.zeros((n,))
+        sl = xopt + np.array([-0.3, -0.01, -0.1])
+        su = xopt + np.array([10.0, 1.0, 10.0])
+        boxproj = lambda x: p_box(x, sl, su)
+        ballproj = lambda x: p_ball(x, xopt, 0.25)
+        d, gnew, _crvmin = ctrsbox(xbase, xopt, g, H, sl, su, [ballproj], Delta)
+        true_d = np.array([-0.22913085, 0.0, -0.09999527])
+        est_min = model_value(g, H, d)
+        true_min = model_value(g, H, true_d)
+        self.assertTrue(est_min <= true_min + 1e-3, 'Sufficient decrease not achieved')
\ No newline at end of file
diff --git a/pybobyqa/trust_region.py b/pybobyqa/trust_region.py
index e97fb1c..16fc078 100644
--- a/pybobyqa/trust_region.py
+++ b/pybobyqa/trust_region.py
@@ -417,6 +417,9 @@ def trsbox_geometry(xbase, c, g, H, lower, upper, Delta, use_fortran=USE_FORTRAN
 
 
 def ctrsbox(xbase, xopt, g, H, sl, su, projections, delta, dykstra_max_iters=100, dykstra_tol=1e-10, gtol=1e-8):
+    if projections is None or len(projections) == 0:
+        return trsbox(xopt, g, H, sl, su, delta)
+
     n = xopt.size
     assert xopt.shape == (n,), "xopt has wrong shape (should be vector)"
     assert xbase.shape == (n,), "xbase and xopt has wrong shape (should be vector)"
@@ -484,7 +487,7 @@ def proj(d0):
 
 
 def ctrsbox_geometry(xbase, xopt, c, g, H, lower, upper, projections, Delta, dykstra_max_iters=100, dykstra_tol=1e-10, gtol=1e-8):
-    # Given a Lagrange polynomial defined by: L(x) = c + g' * (x - xbase) + 0.5*(x-xbase)*H*(x-xbase)
+    # Given a Lagrange polynomial defined by: L(x) = c + g' * (x - xopt) + 0.5*(x-xopt)*H*(x-xopt)
     # Maximise |L(x)| in the feasible region + trust region - that is, solve:
     #   max_x  abs(c + g' * (x - xopt) + 0.5*(x-xopt)*H*(x-xopt))
     #    s.t.  lower <= x <= upper
@@ -495,6 +498,9 @@ def ctrsbox_geometry(xbase, xopt, c, g, H, lower, upper, projections, Delta, dyk
     #   s.t.   lower <= xopt + s <= upper
     #          ||s|| <= Delta
     #          xbase + xopt + s in the intersection of all convex sets C with proj_{C} in the list 'projections'
+    if projections is None or len(projections) == 0:
+        return trsbox_geometry(xopt, g, H, lower, upper, Delta)
+
     smin, gmin, crvmin = ctrsbox(xbase, xopt, g, H, lower, upper, projections, Delta,
                                  dykstra_max_iters=dykstra_max_iters, dykstra_tol=dykstra_tol, gtol=gtol)  # minimise L(x)
     smax, gmax, crvmax = ctrsbox(xbase, xopt, -g, -H, lower, upper, projections, Delta,