diff --git a/pyprima/src/pyprima/__init__.py b/pyprima/src/pyprima/__init__.py
index 577e1d2661..1c7f54160c 100644
--- a/pyprima/src/pyprima/__init__.py
+++ b/pyprima/src/pyprima/__init__.py
@@ -151,6 +151,21 @@ def fixed_fun(x):
             newx[~_fixed_idx] = x
             return original_fun(newx, *args)
         fun = fixed_fun
+        # Adjust linear_constraint
+        if linear_constraint:
+            new_lb = linear_constraint.lb - linear_constraint.A[:, _fixed_idx] @ _fixed_values
+            new_ub = linear_constraint.ub - linear_constraint.A[:, _fixed_idx] @ _fixed_values
+            new_A = linear_constraint.A[:, ~_fixed_idx]
+            linear_constraint = LinearConstraint(new_A, new_lb, new_ub)
+        # Adjust nonlinear constraints
+        if nonlinear_constraint_function:
+            original_nlc_fun = nonlinear_constraint_function
+            def fixed_nlc_fun(x):
+                newx = np.zeros(lenx0)
+                newx[_fixed_idx] = _fixed_values
+                newx[~_fixed_idx] = x
+                return original_nlc_fun(newx, *args)
+            nonlinear_constraint_function = fixed_nlc_fun
 
 
     # Project x0 onto the feasible set
diff --git a/pyprima/src/pyprima/common/_bounds.py b/pyprima/src/pyprima/common/_bounds.py
index 89afe9529f..fdc585a94a 100644
--- a/pyprima/src/pyprima/common/_bounds.py
+++ b/pyprima/src/pyprima/common/_bounds.py
@@ -30,5 +30,9 @@ def process_bounds(bounds, lenx0):
     # If there were fewer bounds than variables, pad the rest with -/+ infinity
     lb = np.concatenate((lb, -np.inf*np.ones(lenx0 - len(lb))))
     ub = np.concatenate((ub, np.inf*np.ones(lenx0 - len(ub))))
+
+    # Check the infeasibility of the bounds.
+    infeasible = np.isposinf(lb) | np.isneginf(ub) | (lb > ub)
+    assert not np.any(infeasible), f"Some of the provided bounds are infeasible. {infeasible=} {lb[infeasible]=}, {ub[infeasible]=}"
     
     return lb, ub
diff --git a/pyprima/tests/test_bounds.py b/pyprima/tests/test_bounds.py
index 72d74db7b7..94505d6bcb 100644
--- a/pyprima/tests/test_bounds.py
+++ b/pyprima/tests/test_bounds.py
@@ -1,5 +1,6 @@
-from pyprima import minimize
+from pyprima import minimize, LinearConstraint as LC, NonlinearConstraint as NLC
 import numpy as np
+import pytest
 
 def test_eliminate_fixed_bounds():
     # Test the logic for detecting and eliminating fixed bounds
@@ -13,3 +14,60 @@ def f(x):
     res = minimize(f, x0=np.array([1, 2, 3, 4, 5]), bounds=bounds)
     assert np.allclose(res.x, np.array([-0.5, -0.5, 1, 0, -0.5]), atol=1e-3)
     assert np.allclose(res.f, 1.75, atol=1e-3)
+
+
+def test_eliminate_fixed_bounds_with_linear_constraints():
+    # Ensure that the logic for fixed bounds also modifies linear constraints
+    # appropriately
+
+    def f(x):
+        return np.sum(x**2)
+
+    lb = [-1, None, None]
+    ub = [-1, None, None]
+    bounds = [(a, b) for a, b in zip(lb, ub)]
+    # Internally, the algorithm should modify lc to be a 1x2 matrix instead of 1x3,
+    # since it will modify x0 and the objective function to eliminate the first
+    # variable. If it fails to modify lc, we will get shape mismatch errors.
+    lc = LC(np.array([1, 1, 1]), lb=9, ub=15)
+    res = minimize(f, x0=np.array([1, 1, 1]), constraints=lc, bounds=bounds)
+    assert np.isclose(res.x[0], -1, atol=1e-6, rtol=1e-6)
+    assert np.isclose(res.x[1], 5, atol=1e-6, rtol=1e-6)
+    assert np.isclose(res.x[2], 5, atol=1e-6, rtol=1e-6)
+    assert np.isclose(res.f, 51, atol=1e-6, rtol=1e-6)
+
+
+def test_eliminate_fixed_bounds_with_nonlinear_constraints():
+    # Ensure that the logic for fixed bounds also modifies the nonlinear constraint
+    # function appropriately
+
+    def f(x):
+        return np.sum(x**2)
+
+    lb = [-1, None, None]
+    ub = [-1, None, None]
+    bounds = [(a, b) for a, b in zip(lb, ub)]
+    x0 = np.array([1, 1, 1])
+    # Have the nonlinear constraint function operate on the last element of x, but be
+    # explicit about the length of x. This ensures that the test is still valid if the
+    # fixed bound is removed. If we simply used x[-1] this test would pass but it
+    # wouldn't actually test if we had properly modified the nonlinear constraint
+    # function after removing the fixed bounds
+    nlc = NLC(lambda x: x[len(x0)-1]**2, lb=9, ub=15)
+    res = minimize(f, x0=x0, constraints=nlc, bounds=bounds)
+    assert np.isclose(res.x[0], -1, atol=1e-6, rtol=1e-6)
+    assert np.isclose(res.x[1], 0, atol=1e-6, rtol=1e-6)
+    assert np.isclose(res.x[2], 3, atol=1e-6, rtol=1e-6)
+    assert np.isclose(res.f, 10, atol=1e-6, rtol=1e-6)
+
+
+def test_infeasible_bounds():
+    def f(x):
+        return np.sum(x**2)
+
+    lb = [1, None, None]
+    ub = [-1, None, None]
+    bounds = [(a, b) for a, b in zip(lb, ub)]
+    with pytest.raises(AssertionError) as excinfo:
+        minimize(f, x0=np.array([1, 2, 3]), bounds=bounds)
+    assert str(excinfo.value) == "Some of the provided bounds are infeasible. infeasible=array([ True, False, False]) lb[infeasible]=array([1.]), ub[infeasible]=array([-1.])"
diff --git a/python/prima/__init__.py b/python/prima/__init__.py
index 33cbead755..0134755f02 100644
--- a/python/prima/__init__.py
+++ b/python/prima/__init__.py
@@ -10,6 +10,7 @@
 from enum import Enum
 import numpy as np
 from ._common import _project
+from ._common import get_arrays_tol
 
 
 class ConstraintType(Enum):
@@ -115,6 +116,59 @@ def minimize(fun, x0, args=(), method=None, bounds=None, constraints=(), callbac
 
     lb, ub = process_bounds(bounds, lenx0)
 
+    # Check which variables are fixed and eliminate them from the problem.
+    # Save the indices and values so that we can call the original function with
+    # an array of the appropriate size, and so that we can add the fixed values to the
+    # result when COBYLA returns.
+    tol = get_arrays_tol(lb, ub)
+    _fixed_idx = (
+        (lb <= ub)
+        & (np.abs(lb - ub) < tol)
+    )
+    if any(_fixed_idx):
+        _fixed_values = 0.5 * (
+            lb[_fixed_idx] + ub[_fixed_idx]
+        )
+        _fixed_values = np.clip(
+            _fixed_values,
+            lb[_fixed_idx],
+            ub[_fixed_idx],
+        )
+        if isinstance(x0, np.ndarray):
+            x0 = np.array(x0)[~_fixed_idx]
+        elif x0_is_scalar:
+            raise Exception("You have provided a scalar x with fixed bounds, there is" \
+            "no optimization to be done.")
+        else:
+            # In this case x is presumably a list, so we turn it into a numpy array
+            # for the convenience of indexing and then turn it back into a list
+            x0 = np.array(x0)[~_fixed_idx].tolist()
+        lb = lb[~_fixed_idx]
+        ub = ub[~_fixed_idx]
+        original_fun = fun
+        def fixed_fun(x):
+            newx = np.zeros(lenx0)
+            newx[_fixed_idx] = _fixed_values
+            newx[~_fixed_idx] = x
+            return original_fun(newx, *args)
+        fun = fixed_fun
+        # Adjust linear_constraint
+        if linear_constraint:
+            new_lb = linear_constraint.lb - linear_constraint.A[:, _fixed_idx] @ _fixed_values
+            new_ub = linear_constraint.ub - linear_constraint.A[:, _fixed_idx] @ _fixed_values
+            new_A = linear_constraint.A[:, ~_fixed_idx]
+            linear_constraint = LinearConstraint(new_A, new_lb, new_ub)
+        # Adjust nonlinear constraints
+        if nonlinear_constraint_function:
+            original_nlc_fun = nonlinear_constraint_function
+            def fixed_nlc_fun(x):
+                newx = np.zeros(lenx0)
+                newx[_fixed_idx] = _fixed_values
+                newx[~_fixed_idx] = x
+                return original_nlc_fun(newx, *args)
+            nonlinear_constraint_function = fixed_nlc_fun
+
+
     # Project x0 onto the feasible set
     if nonlinear_constraint_function is None:
         result = _project(x0, lb, ub, {"linear": linear_constraint, "nonlinear": None})
@@ -142,7 +196,7 @@ def minimize(fun, x0, args=(), method=None, bounds=None, constraints=(), callbac
         options['nlconstr0'] = nlconstr0
         options['m_nlcon'] = len(nlconstr0)
 
-    return _minimize(
+    result = _minimize(
         fun,
         x0,
         args,
@@ -157,3 +211,10 @@ def minimize(fun, x0, args=(), method=None, bounds=None, constraints=(), callbac
         callback,
         options
     )
+
+    if any(_fixed_idx):
+        newx = np.zeros(lenx0)
+        newx[_fixed_idx] = _fixed_values
+        newx[~_fixed_idx] = result.x
+        result.x = newx
+    return result
diff --git a/python/prima/_bounds.py b/python/prima/_bounds.py
index 89afe9529f..6ca465a811 100644
--- a/python/prima/_bounds.py
+++ b/python/prima/_bounds.py
@@ -30,5 +30,9 @@ def process_bounds(bounds, lenx0):
     # If there were fewer bounds than variables, pad the rest with -/+ infinity
     lb = np.concatenate((lb, -np.inf*np.ones(lenx0 - len(lb))))
     ub = np.concatenate((ub, np.inf*np.ones(lenx0 - len(ub))))
-    
+
+    # Check the infeasibility of the bounds.
+    infeasible = np.isposinf(lb) | np.isneginf(ub) | (lb > ub)
+    assert not np.any(infeasible), f"Some of the provided bounds are infeasible. {infeasible=} {lb[infeasible]=}, {ub[infeasible]=}"
+
     return lb, ub
diff --git a/python/prima/_common.py b/python/prima/_common.py
index 950c247bf1..d94f1e8d0d 100644
--- a/python/prima/_common.py
+++ b/python/prima/_common.py
@@ -165,3 +165,32 @@ def _project(x0, lb, ub, constraints):
             return OptimizeResult(x=x0_c)
 
     return OptimizeResult(x=x0_c)
+
+
+def get_arrays_tol(*arrays):
+    """
+    Get a relative tolerance for a set of arrays. Borrowed from COBYQA
+
+    Parameters
+    ----------
+    *arrays: tuple
+        Set of `numpy.ndarray` to get the tolerance for.
+
+    Returns
+    -------
+    float
+        Relative tolerance for the set of arrays.
+
+    Raises
+    ------
+    ValueError
+        If no array is provided.
+    """
+    if len(arrays) == 0:
+        raise ValueError("At least one array must be provided.")
+    size = max(array.size for array in arrays)
+    weight = max(
+        np.max(np.abs(array[np.isfinite(array)]), initial=1.0)
+        for array in arrays
+    )
+    return 10.0 * eps * max(size, 1.0) * weight
diff --git a/python/tests/test_bounds.py b/python/tests/test_bounds.py
new file mode 100644
index 0000000000..6beeba0087
--- /dev/null
+++ b/python/tests/test_bounds.py
@@ -0,0 +1,73 @@
+from prima import minimize, LinearConstraint as LC, NonlinearConstraint as NLC
+import numpy as np
+import pytest
+
+def test_eliminate_fixed_bounds():
+    # Test the logic for detecting and eliminating fixed bounds
+
+    def f(x):
+        return np.sum(x**2)
+
+    lb = [-1, None, 1, None, -0.5]
+    ub = [-0.5, -0.5, None, None, -0.5]
+    bounds = [(a, b) for a, b in zip(lb, ub)]    
+    res = minimize(f, x0=np.array([1, 2, 3, 4, 5]), bounds=bounds)
+    assert np.allclose(res.x, np.array([-0.5, -0.5, 1, 0, -0.5]), atol=1e-3)
+    assert np.allclose(res.fun, 1.75, atol=1e-3)
+
+
+def test_eliminate_fixed_bounds_with_linear_constraints():
+    # Ensure that the logic for fixed bounds also modifies linear constraints
+    # appropriately
+
+    def f(x):
+        return np.sum(x**2)
+
+    lb = [-1, None, None]
+    ub = [-1, None, None]
+    bounds = [(a, b) for a, b in zip(lb, ub)]
+    # Internally, the algorithm should modify lc to be a 1x2 matrix instead of 1x3,
+    # since it will modify x0 and the objective function to eliminate the first
+    # variable. If it fails to modify lc, we will get shape mismatch errors.
+    lc = LC(np.array([1, 1, 1]), lb=9, ub=15)
+    res = minimize(f, x0=np.array([1, 1, 1]), constraints=lc, bounds=bounds)
+    assert np.isclose(res.x[0], -1, atol=1e-6, rtol=1e-6)
+    assert np.isclose(res.x[1], 5, atol=1e-6, rtol=1e-6)
+    assert np.isclose(res.x[2], 5, atol=1e-6, rtol=1e-6)
+    assert np.isclose(res.fun, 51, atol=1e-6, rtol=1e-6)
+
+
+def test_eliminate_fixed_bounds_with_nonlinear_constraints():
+    # Ensure that the logic for fixed bounds also modifies the nonlinear constraint
+    # function appropriately
+
+    def f(x):
+        return np.sum(x**2)
+
+    lb = [-1, None, None]
+    ub = [-1, None, None]
+    bounds = [(a, b) for a, b in zip(lb, ub)]
+    x0 = np.array([1, 1, 1])
+    # Have the nonlinear constraint function operate on the last element of x, but be
+    # explicit about the length of x. This ensures that the test is still valid if the
+    # fixed bound is removed. If we simply used x[-1] this test would pass but it
+    # wouldn't actually test if we had properly modified the nonlinear constraint
+    # function after removing the fixed bounds
+    nlc = NLC(lambda x: x[len(x0)-1]**2, lb=9, ub=15)
+    res = minimize(f, x0=x0, constraints=nlc, bounds=bounds)
+    assert np.isclose(res.x[0], -1, atol=1e-6, rtol=1e-6)
+    assert np.isclose(res.x[1], 0, atol=1e-6, rtol=1e-6)
+    assert np.isclose(res.x[2], 3, atol=1e-6, rtol=1e-6)
+    assert np.isclose(res.fun, 10, atol=1e-6, rtol=1e-6)
+
+
+def test_infeasible_bounds():
+    def f(x):
+        return np.sum(x**2)
+
+    lb = [1, None, None]
+    ub = [-1, None, None]
+    bounds = [(a, b) for a, b in zip(lb, ub)]
+    with pytest.raises(AssertionError) as excinfo:
+        minimize(f, x0=np.array([1, 2, 3]), bounds=bounds)
+    assert str(excinfo.value) == "Some of the provided bounds are infeasible. infeasible=array([ True, False, False]) lb[infeasible]=array([1.]), ub[infeasible]=array([-1.])"