diff --git a/.github/workflows/unopy-test-example.yml b/.github/workflows/unopy-test-example.yml
index 8e12a1e5dd..759d020b1f 100644
--- a/.github/workflows/unopy-test-example.yml
+++ b/.github/workflows/unopy-test-example.yml
@@ -43,8 +43,7 @@ jobs:
           python-version: ${{ matrix.python-version }}
           architecture: ${{ matrix.architecture }}
           
-      - name: Install pybind11 on Linux
-        if: runner.os == 'Linux'
+      - name: Install pybind11
         run: sudo apt install -y python3-pybind11 pybind11-dev
 
       - name: Download dependencies
@@ -55,4 +54,4 @@ jobs:
 
       - name: Run example
         working-directory: ${{github.workspace}}/interfaces/Python/example
-        run: python example_hs015.py
\ No newline at end of file
+        run: python example_hs015.py
diff --git a/.github/workflows/unopy-tests-gemseo.yml b/.github/workflows/unopy-tests-gemseo.yml
new file mode 100644
index 0000000000..9dfaa70591
--- /dev/null
+++ b/.github/workflows/unopy-tests-gemseo.yml
@@ -0,0 +1,59 @@
+name: Test unopy via GEMSEO interface
+
+on:
+  push:
+    branches: [ "main" ]
+    paths-ignore:
+      - '*.md'
+      - 'LICENSE'
+      - '*.cff'
+      - '*.yml'
+      - '*.yaml'
+      - 'docs/**'
+  pull_request:
+    branches: [ "main" ]
+    paths-ignore:
+      - '*.md'
+      - 'LICENSE'
+      - '*.cff'
+      - '*.yml'
+      - '*.yaml'
+      - 'docs/**'
+
+env:
+  BUILD_TYPE: Debug
+
+jobs:
+  build:
+    runs-on: ${{ matrix.os }}
+    strategy:
+      matrix:
+        os: [ubuntu-latest]
+        architecture: [x64]
+        python-version: ["3.12"]
+
+    steps:
+      - uses: actions/checkout@v4
+      
+      - name: Set up Python
+        # This is the version of the action for setting up Python, not the Python version.
+        uses: actions/setup-python@v5
+        with:
+          # Semantic version range syntax or exact version of a Python version
+          python-version: ${{ matrix.python-version }}
+          architecture: ${{ matrix.architecture }}
+          
+      - name: Install pybind11
+        run: sudo apt install -y python3-pybind11 pybind11-dev
+
+      - name: Install pytest, gemseo, and scipy
+        run: pip install pytest gemseo scipy
+
+      - name: Download dependencies
+        run: bash dependencies/scripts/download_dependencies.sh
+
+      - name: Compile and install unopy
+        run: CMAKE_BUILD_PARALLEL_LEVEL=$(nproc) pip install . -v
+
+      - name: Run tests
+        run: python -m pytest -s --import-mode=importlib ${{github.workspace}}/interfaces/Python/gemseo_uno/test_gemseo_uno.py
diff --git a/.github/workflows/unopy-tests-scipy.yml b/.github/workflows/unopy-tests-scipy.yml
new file mode 100644
index 0000000000..9b946cc700
--- /dev/null
+++ b/.github/workflows/unopy-tests-scipy.yml
@@ -0,0 +1,59 @@
+name: Test unopy via Scipy interface
+
+on:
+  push:
+    branches: [ "main" ]
+    paths-ignore:
+      - '*.md'
+      - 'LICENSE'
+      - '*.cff'
+      - '*.yml'
+      - '*.yaml'
+      - 'docs/**'
+  pull_request:
+    branches: [ "main" ]
+    paths-ignore:
+      - '*.md'
+      - 'LICENSE'
+      - '*.cff'
+      - '*.yml'
+      - '*.yaml'
+      - 'docs/**'
+
+env:
+  BUILD_TYPE: Debug
+
+jobs:
+  build:
+    runs-on: ${{ matrix.os }}
+    strategy:
+      matrix:
+        os: [ubuntu-latest]
+        architecture: [x64]
+        python-version: ["3.12"]
+
+    steps:
+      - uses: actions/checkout@v4
+      
+      - name: Set up Python
+        # This is the version of the action for setting up Python, not the Python version.
+        uses: actions/setup-python@v5
+        with:
+          # Semantic version range syntax or exact version of a Python version
+          python-version: ${{ matrix.python-version }}
+          architecture: ${{ matrix.architecture }}
+          
+      - name: Install pybind11
+        run: sudo apt install -y python3-pybind11 pybind11-dev
+        
+      - name: Install pytest and scipy
+        run: pip install pytest scipy
+
+      - name: Download dependencies
+        run: bash dependencies/scripts/download_dependencies.sh
+
+      - name: Compile and install unopy
+        run: CMAKE_BUILD_PARALLEL_LEVEL=$(nproc) pip install . -v
+
+      - name: Run tests
+        run: python -m pytest --import-mode=importlib ${{github.workspace}}/interfaces/Python/scipy_interface/test_scipy_interface.py
diff --git a/CMakeLists.txt b/CMakeLists.txt
index f2762db7b8..d7c656ead6 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -425,6 +425,9 @@ if(SKBUILD)
       target_compile_definitions(unopy PRIVATE PYBIND11_DETAILED_ERROR_MESSAGES)
    endif()
 
+   install(DIRECTORY interfaces/Python/unopy/ DESTINATION unopy)
+   install(DIRECTORY interfaces/Python/scipy_interface/ DESTINATION scipy_interface)
+   install(DIRECTORY interfaces/Python/gemseo_uno/ DESTINATION gemseo_uno)
    install(TARGETS unopy
            LIBRARY DESTINATION unopy    # Linux/macOS: unopy/unopy.so
            RUNTIME DESTINATION unopy    # Windows: unopy/unopy.pyd
diff --git a/interfaces/AMPL/AMPLModel.cpp b/interfaces/AMPL/AMPLModel.cpp
index 3dc9ed9cfa..6b59554ba8 100644
--- a/interfaces/AMPL/AMPLModel.cpp
+++ b/interfaces/AMPL/AMPLModel.cpp
@@ -200,7 +200,7 @@ namespace uno {
       fint error_flag = 0;
       this->asl->p.Jacval(this->asl, const_cast<double*>(x.data()), jacobian_values, &error_flag);
       if (0 < error_flag) {
-         throw GradientEvaluationError();
+         throw JacobianEvaluationError();
       }
       ++this->number_model_evaluations.jacobian;
    }
diff --git a/interfaces/C/Uno_C_API.cpp b/interfaces/C/Uno_C_API.cpp
index 8edfde7817..f81a5bac15 100644
--- a/interfaces/C/Uno_C_API.cpp
+++ b/interfaces/C/Uno_C_API.cpp
@@ -152,7 +152,7 @@ class UnoModel: public Model {
          const uno_int return_code = this->user_model.jacobian(this->user_model.number_variables,
             this->user_model.number_jacobian_nonzeros, x.data(), jacobian_values, this->user_model.user_data);
          if (0 < return_code) {
-            throw GradientEvaluationError();
+            throw JacobianEvaluationError();
          }
          ++this->number_model_evaluations.jacobian;
       }
@@ -189,7 +189,7 @@ class UnoModel: public Model {
          const uno_int return_code = this->user_model.jacobian_operator(this->user_model.number_variables,
             this->user_model.number_constraints, x, true, vector, result, this->user_model.user_data);
          if (0 < return_code) {
-            throw GradientEvaluationError();
+            throw JacobianEvaluationError();
          }
       }
       else {
@@ -202,7 +202,7 @@ class UnoModel: public Model {
          const uno_int return_code = this->user_model.jacobian_transposed_operator(this->user_model.number_variables,
             this->user_model.number_constraints, x, true, vector, result, this->user_model.user_data);
          if (0 < return_code) {
-            throw GradientEvaluationError();
+            throw JacobianEvaluationError();
          }
       }
       else {
@@ -1265,4 +1265,4 @@ void uno_destroy_solver(void* solver) {
       }
       delete uno_solver;
    }
-}
\ No newline at end of file
+}
diff --git a/interfaces/Python/cpp_classes/PythonModel.cpp b/interfaces/Python/cpp_classes/PythonModel.cpp
index 8b6ac969d0..11a67e4c27 100644
--- a/interfaces/Python/cpp_classes/PythonModel.cpp
+++ b/interfaces/Python/cpp_classes/PythonModel.cpp
@@ -142,7 +142,7 @@ namespace uno {
             ++this->number_model_evaluations.jacobian;
          }
          catch (const std::exception&) {
-            throw GradientEvaluationError();
+            throw JacobianEvaluationError();
          }
       }
    }
@@ -192,7 +192,7 @@ namespace uno {
             (*this->user_model.jacobian_operator)(x_py, true, vector_py, result_py);
          }
          catch (const std::exception&) {
-            throw GradientEvaluationError();
+            throw JacobianEvaluationError();
          }
       }
       else {
@@ -211,7 +211,7 @@ namespace uno {
             (*this->user_model.jacobian_transposed_operator)(x_py, true, vector_py, result_py);
          }
          catch (const std::exception&) {
-            throw GradientEvaluationError();
+            throw JacobianEvaluationError();
          }
       }
       else {
diff --git a/interfaces/Python/example/example_hs015.py b/interfaces/Python/example/example_hs015.py
index 5792ce387c..10b65f8800 100644
--- a/interfaces/Python/example/example_hs015.py
+++ b/interfaces/Python/example/example_hs015.py
@@ -2,6 +2,7 @@
 # Licensed under the MIT license. See LICENSE file in the project directory for details.
 
 import unopy
+
 Inf = float("inf")
 
 # hs015.mod
diff --git a/interfaces/Python/gemseo_uno/__init__.py b/interfaces/Python/gemseo_uno/__init__.py
new file mode 100644
index 0000000000..e69de29bb2
diff --git a/interfaces/Python/gemseo_uno/gemseo_uno.py b/interfaces/Python/gemseo_uno/gemseo_uno.py
new file mode 100644
index 0000000000..c34db65563
--- /dev/null
+++ b/interfaces/Python/gemseo_uno/gemseo_uno.py
@@ -0,0 +1,147 @@
+# Copyright 2021 IRT Saint Exupéry, https://www.irt-saintexupery.com
+#
+# This program is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License version 3 as published by the Free Software Foundation.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program; if not, write to the Free Software Foundation,
+# Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+# Contributors:
+#    INITIAL AUTHORS - initial API and implementation and/or initial
+#                           documentation
+#        :author: François Gallard
+#    OTHER AUTHORS   - MACROSCOPIC CHANGES
+"""The library of Uno constrained gradient-based optimization algorithms."""
+
+from __future__ import annotations
+
+from dataclasses import dataclass
+from typing import Any
+from typing import ClassVar
+from typing import TYPE_CHECKING
+
+from gemseo.algos.design_space_utils import get_value_and_bounds
+from gemseo.algos.opt.base_optimization_library import (
+    OptimizationAlgorithmDescription,
+    BaseOptimizationLibrary,
+)
+from numpy import isfinite
+from numpy import real
+
+from gemseo_uno.settings.base_uno_settings import UNO_Settings
+from scipy_interface.scipy_uno import minimize
+
+if TYPE_CHECKING:
+    from gemseo.algos.optimization_problem import OptimizationProblem
+
+
+@dataclass
+class UnoAlgorithmDescription(OptimizationAlgorithmDescription):
+    """The description of the Uno constrained grdient based optimization library."""
+
+    library_name: str = "Uno"
+    """The library name."""
+
+    handle_equality_constraints: bool = True
+    """Whether the optimization algorithm handles equality constraints."""
+
+    handle_inequality_constraints: bool = True
+    """Whether the optimization algorithm handles inequality constraints."""
+
+    positive_constraints: bool = True
+    """Whether the optimization algorithm requires positive constraints."""
+
+    require_gradient: bool = True
+    """Whether the optimization algorithm requires the gradient."""
+
+    Settings: type[UNO_Settings] = UNO_Settings
+    """The option validation model for Uno optimization library."""
+
+    website: str = "https://unosolver.readthedocs.io/en/latest/"
+    """The website of the wrapped library or algorithm."""
+
+
+class UnoOpt(BaseOptimizationLibrary[UNO_Settings]):
+    """The library of Uno optimization algorithms."""
+
+    ALGORITHM_INFOS: ClassVar[dict[str, UnoAlgorithmDescription]] = {
+        "UNO_Filter_SQP": UnoAlgorithmDescription(
+            algorithm_name="UNO_Filter_SQP",
+            description=(
+                "Sequential Quadratic Programming (SQP) "
+                "implemented in the Uno library"
+            ),
+            internal_algorithm_name="filtersqp",
+            Settings=UNO_Settings,
+        ),
+        "UNO_Filter_SLP": UnoAlgorithmDescription(
+            algorithm_name="UNO_Filter_SLP",
+            description=(
+                "Sequential Linear Programming (SLP) " "implemented in the Uno library"
+            ),
+            internal_algorithm_name="filterslp",
+            Settings=UNO_Settings,
+        ),
+        "UNO_Funnel_SQP": UnoAlgorithmDescription(
+            algorithm_name="UNO_Funnel_SQP",
+            description=(
+                "Funnel Sequential Quadratic Programming (SQP)"
+                "implemented in the Uno library"
+            ),
+            internal_algorithm_name="funnelsqp",
+            Settings=UNO_Settings,
+        ),
+        "UNO_IPOPT": UnoAlgorithmDescription(
+            algorithm_name="UNO_IPOPT",
+            description=(
+                "Interior Point Optimization (IPOPT)" "implemented in the Uno library"
+            ),
+            internal_algorithm_name="ipopt",
+            Settings=UNO_Settings,
+        ),
+    }
+
+    def _run(self, problem: OptimizationProblem) -> tuple[str, Any]:
+        # Get the normalized bounds:
+        x_0, l_b, u_b = get_value_and_bounds(
+            problem.design_space, self._settings.normalize_design_space
+        )
+        # Replace infinite values with None:
+        l_b = [val if isfinite(val) else None for val in l_b]
+        u_b = [val if isfinite(val) else None for val in u_b]
+        bounds = list(zip(l_b, u_b, strict=False))
+
+        # Get constraint in SciPy format
+        scipy_constraints = [
+            {
+                "type": constraint.f_type,
+                "fun": constraint.evaluate,
+                "jac": constraint.jac,
+            }
+            for constraint in self._get_right_sign_constraints(problem)
+        ]
+
+        # Filter settings to get only the uno ones
+        settings_ = self._filter_settings(self._settings.model_dump(), UNO_Settings)
+
+        # Deactivate stopping criteria which are handled by GEMSEO
+        tolerance = 0.0
+
+        opt_result = minimize(
+            fun=lambda x: real(problem.objective.evaluate(x)),
+            jac=problem.objective.jac,
+            x0=x_0,
+            method=self.ALGORITHM_INFOS[self._algo_name].internal_algorithm_name,
+            bounds=bounds,
+            constraints=scipy_constraints,
+            options=settings_,
+            tol=tolerance,
+        )
+
+        return opt_result.message, opt_result.status
diff --git a/interfaces/Python/gemseo_uno/settings/__init__.py b/interfaces/Python/gemseo_uno/settings/__init__.py
new file mode 100644
index 0000000000..e69de29bb2
diff --git a/interfaces/Python/gemseo_uno/settings/base_uno_settings.py b/interfaces/Python/gemseo_uno/settings/base_uno_settings.py
new file mode 100644
index 0000000000..3d25ac158b
--- /dev/null
+++ b/interfaces/Python/gemseo_uno/settings/base_uno_settings.py
@@ -0,0 +1,117 @@
+# Copyright 2021 IRT Saint Exupéry, https://www.irt-saintexupery.com
+#
+# This program is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public
+# License version 3 as published by the Free Software Foundation.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program; if not, write to the Free Software Foundation,
+# Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+"""Settings for the SciPy COBYQA algorithm."""
+
+from __future__ import annotations
+
+from gemseo.algos.opt.base_gradient_based_algorithm_settings import (
+    BaseGradientBasedAlgorithmSettings,
+)
+from gemseo.algos.opt.base_optimizer_settings import BaseOptimizerSettings
+from pydantic import PositiveFloat  # noqa:TC002
+
+
+class UNO_Settings(
+    BaseOptimizerSettings, BaseGradientBasedAlgorithmSettings
+):  # noqa: N801
+    """Settings for the UNO algorithm."""
+
+    _TARGET_CLASS_NAME = "UNO"
+    initial_tr_radius: PositiveFloat | None = None
+    primal_tolerance: float | None = None
+    dual_tolerance: float | None = None
+    loose_primal_tolerance: float | None = None
+    loose_dual_tolerance: float | None = None
+    loose_tolerance_consecutive_iteration_threshold: int | None = None
+    time_limit: float | None = None
+    print_solution: bool | None = None
+    unbounded_objective_threshold: float | None = None
+    enforce_linear_constraints: bool | None = None
+    logger: str | None = None
+    constraint_relaxation_strategy: str | None = None
+    inequality_handling_method: str | None = None
+    globalization_mechanism: str | None = None
+    globalization_strategy: str | None = None
+    hessian_model: str | None = None
+    inertia_correction_strategy: str | None = None
+    scale_functions: bool | None = None
+    function_scaling_threshold: float | None = None
+    function_scaling_factor: float | None = None
+    scale_residuals: bool | None = None
+    progress_norm: str | None = None
+    residual_norm: str | None = None
+    residual_scaling_threshold: float | None = None
+    protect_actual_reduction_against_roundoff: bool | None = None
+    print_subproblem: bool | None = None
+    armijo_decrease_fraction: float | None = None
+    armijo_tolerance: float | None = None
+    switching_delta: float | None = None
+    switching_infeasibility_exponent: float | None = None
+    filter_type: str | None = None
+    filter_beta: float | None = None
+    filter_gamma: float | None = None
+    filter_ubd: float | None = None
+    filter_fact: float | None = None
+    filter_capacity: int | None = None
+    filter_sufficient_infeasibility_decrease_factor: float | None = None
+    nonmonotone_filter_number_dominated_entries: int | None = None
+    funnel_kappa: float | None = None
+    funnel_beta: float | None = None
+    funnel_gamma: float | None = None
+    funnel_ubd: float | None = None
+    funnel_fact: float | None = None
+    funnel_update_strategy: int | None = None
+    funnel_require_acceptance_wrt_current_iterate: bool | None = None
+    LS_backtracking_ratio: float | None = None
+    LS_min_step_length: float | None = None
+    LS_scale_duals_with_step_length: bool | None = None
+    regularization_failure_threshold: float | None = None
+    regularization_initial_value: float | None = None
+    regularization_increase_factor: float | None = None
+    primal_regularization_initial_factor: float | None = None
+    dual_regularization_fraction: float | None = None
+    primal_regularization_lb: float | None = None
+    primal_regularization_decrease_factor: float | None = None
+    primal_regularization_fast_increase_factor: float | None = None
+    primal_regularization_slow_increase_factor: float | None = None
+    threshold_unsuccessful_attempts: int | None = None
+    TR_radius: float | None = None
+    TR_increase_factor: float | None = None
+    TR_decrease_factor: float | None = None
+    TR_aggressive_decrease_factor: float | None = None
+    TR_activity_tolerance: float | None = None
+    TR_min_radius: float | None = None
+    TR_radius_reset_threshold: float | None = None
+    switch_to_optimality_requires_linearized_feasibility: bool | None = None
+    l1_constraint_violation_coefficient: float | None = None
+    barrier_initial_parameter: float | None = None
+    barrier_default_multiplier: float | None = None
+    barrier_tau_min: float | None = None
+    barrier_k_sigma: float | None = None
+    barrier_smax: float | None = None
+    barrier_k_mu: float | None = None
+    barrier_theta_mu: float | None = None
+    barrier_k_epsilon: float | None = None
+    barrier_update_fraction: float | None = None
+    barrier_regularization_exponent: float | None = None
+    barrier_small_direction_factor: float | None = None
+    barrier_push_variable_to_interior_k1: float | None = None
+    barrier_push_variable_to_interior_k2: float | None = None
+    barrier_damping_factor: float | None = None
+    least_square_multiplier_max_norm: float | None = None
+    BQPD_kmax_heuristic: str | None = None
+    QP_solver: str | None = None
+    LP_solver: str | None = None
+    linear_solver: str | None = None
diff --git a/interfaces/Python/gemseo_uno/test_gemseo_uno.py b/interfaces/Python/gemseo_uno/test_gemseo_uno.py
new file mode 100644
index 0000000000..84813119d4
--- /dev/null
+++ b/interfaces/Python/gemseo_uno/test_gemseo_uno.py
@@ -0,0 +1,22 @@
+import pytest
+from gemseo.problems.optimization.power_2 import Power2
+from numpy import allclose
+
+from gemseo_uno.gemseo_uno import UnoOpt
+from gemseo_uno.settings.base_uno_settings import UNO_Settings
+
+
+@pytest.mark.parametrize(
+    "method", ["UNO_Filter_SQP", "UNO_Funnel_SQP"]
+)
+def test_power2(method):
+    problem = Power2()
+    problem.preprocess_functions()
+    res = UnoOpt(method).execute(
+        problem, settings_model=UNO_Settings(xtol_rel=1e-4, max_iter=50)
+    )
+    assert res.is_feasible
+    assert allclose(
+        res.x_opt, [0.5 ** (1 / 3), 0.5 ** (1 / 3), 0.9 ** (1 / 3)], rtol=1e-3
+    )
+    assert 3 <= len(problem.database) <= 50
diff --git a/interfaces/Python/scipy_interface/__init__.py b/interfaces/Python/scipy_interface/__init__.py
new file mode 100644
index 0000000000..e69de29bb2
diff --git a/interfaces/Python/scipy_interface/scipy_uno.py b/interfaces/Python/scipy_interface/scipy_uno.py
new file mode 100644
index 0000000000..b9157fa8c9
--- /dev/null
+++ b/interfaces/Python/scipy_interface/scipy_uno.py
@@ -0,0 +1,438 @@
+# Copyright (c) 2026 Francois Gallard, Jean-Christophe Giret
+# Licensed under the MIT license. See LICENSE file in the project directory for details.
+import numpy as np
+import unopy
+from numbers import Number
+from scipy.optimize import LinearConstraint, NonlinearConstraint, OptimizeResult
+from scipy.optimize._minimize import standardize_bounds, _validate_bounds, Bounds
+from typing import Iterable, Sequence, Any
+
+Inf = float("inf")
+AVAILABLE_METHODS = ["filterslp", "filtersqp", "funnelsqp", "ipopt"]
+
+
+def minimize(
+    fun: callable,
+    x0: np.ndarray,
+    args: tuple = (),
+    method: str = "filtersqp",
+    jac: callable = None,
+    hess: callable = None,
+    bounds: Iterable | None | Bounds = None,
+    constraints: Iterable[dict[str:Any] | NonlinearConstraint | LinearConstraint] = (),
+    tol: float | None = None,
+    options: dict | None = None,
+) -> OptimizeResult:
+    """A scipy.optimize.minimize like interface for Uno.
+
+    Minimization of scalar function of one or more variables under bounds and general constraints.
+
+    Parameters
+    ----------
+    fun : callable
+        The objective function to be minimized::
+
+            fun(x, *args) -> float
+
+        where ""x"" is a 1-D array with shape (n,) and ""args""
+        is a tuple of the fixed parameters needed to completely
+        specify the function.
+
+        Suppose the callable has signature ""f0(x, *my_args, **my_kwargs)"", where
+        ""my_args"" and ""my_kwargs"" are required positional and keyword arguments.
+        Rather than passing ""f0"" as the callable, wrap it to accept
+        only ""x""; e.g., pass ""fun=lambda x: f0(x, *my_args, **my_kwargs)"" as the
+        callable, where ""my_args"" (tuple) and ""my_kwargs"" (dict) have been
+        gathered before invoking this function.
+    x0 : ndarray, shape (n,)
+        Initial guess. Array of real elements of size (n,),
+        where ""n"" is the number of independent variables.
+    args : tuple, optional
+        Extra arguments passed to the objective function and its
+        derivatives ("fun", "jac" functions).
+    method : str
+        The UNO preset.  Should be one of
+
+        - 'filtersqp'
+        - 'filterslp
+        - "funnelsqp"
+        - "ipopt"
+
+    jac : {callable}, optional
+        Method for computing the gradient vector.
+        It should be a function that returns the gradient
+        vector::
+
+            jac(x, *args) -> array_like, shape (n,)
+
+        where ""x"" is an array with shape (n,) and ""args"" is a tuple with
+        the fixed parameters.
+
+    hess : callable, optional
+        Hessian of the objective function. Signature: hess(x, *args) -> ndarray (n, n).
+        If provided, constraint Hessians are also read from NonlinearConstraint.hess
+        attributes. The constraint Hessian callable has signature:
+        hess(x, v) -> ndarray (n, n), where v is the vector of constraint multipliers.
+
+    bounds : sequence or "Bounds", optional
+        Bounds on variables. There are two ways to specify the bounds:
+
+        1. Instance of "Bounds" class.
+        2. Sequence of ""(min, max)"" pairs for each element in "x". None
+           is used to specify no bound.
+
+    constraints : {Constraint, dict} or List of {Constraint, dict}, optional
+        Constraints definition.
+
+        Available constraints are:
+
+        - "LinearConstraint"
+        - "NonlinearConstraint"
+        - list of dictionaries.
+
+        For lists of dictionaries, each dictionary with fields:
+
+        type : str
+            Constraint type: 'eq' for equality, 'ineq' for inequality.
+        fun : callable
+            The function defining the constraint.
+        jac : callable, optional
+            The Jacobian of "fun" (only for SLSQP).
+        args : sequence, optional
+            Extra arguments to be passed to the function and Jacobian.
+
+        Equality constraint means that the constraint function result is to
+        be zero whereas inequality means that it is to be non-negative.
+
+    tol : float, optional
+        Tolerance for termination. When "tol" is specified, the selected
+        minimization algorithm sets some relevant solver-specific tolerance(s)
+        equal to "tol".
+
+    options : dict, optional
+        A dictionary of solver options.
+
+        The available options and their types are:
+
+        "primal_tolerance" : float
+        "dual_tolerance" : float
+        "loose_primal_tolerance" : float
+        "loose_dual_tolerance" : float
+        "loose_tolerance_consecutive_iteration_threshold" : int
+        "max_iterations" or "max_iter: int
+        "time_limit" : float
+        "print_solution" : bool
+        "unbounded_objective_threshold" : float
+        "enforce_linear_constraints" : bool
+        "logger" : str,
+        "constraint_relaxation_strategy" : str,
+        "inequality_handling_method" : str,
+        "globalization_mechanism" : str
+        "globalization_strategy" : str,
+        "hessian_model" : str,
+        "inertia_correction_strategy" : str,
+        "scale_functions" : bool
+        "function_scaling_threshold" : float
+        "function_scaling_factor" : float
+        "scale_residuals" : bool
+        "progress_norm" : str,
+        "residual_norm" : str,
+        "residual_scaling_threshold" : float
+        "protect_actual_reduction_against_roundoff" : bool
+        "print_subproblem" : bool
+        "armijo_decrease_fraction" : float
+        "armijo_tolerance" : float
+        "switching_delta" : float
+        "switching_infeasibility_exponent" : float
+        "filter_type" : str,
+        "filter_beta" : float
+        "filter_gamma" : float
+        "filter_ubd" : float
+        "filter_fact" : float
+        "filter_capacity" : int
+        "filter_sufficient_infeasibility_decrease_factor" : float
+        "nonmonotone_filter_number_dominated_entries" : int
+        "funnel_kappa" : float
+        "funnel_beta" : float
+        "funnel_gamma" : float
+        "funnel_ubd" : float
+        "funnel_fact" : float
+        "funnel_update_strategy" : int
+        "funnel_require_acceptance_wrt_current_iterate" : bool
+        "LS_backtracking_ratio" : float
+        "LS_min_step_length" : float
+        "LS_scale_duals_with_step_length" : bool
+        "regularization_failure_threshold" : float
+        "regularization_initial_value" : float
+        "regularization_increase_factor" : float
+        "primal_regularization_initial_factor" : float
+        "dual_regularization_fraction" : float
+        "primal_regularization_lb" : float
+        "primal_regularization_decrease_factor" : float
+        "primal_regularization_fast_increase_factor" : float
+        "primal_regularization_slow_increase_factor" : float
+        "threshold_unsuccessful_attempts" : int
+        "TR_radius" : float
+        "TR_increase_factor" : float
+        "TR_decrease_factor" : float
+        "TR_aggressive_decrease_factor" : float
+        "TR_activity_tolerance" : float
+        "TR_min_radius" : float
+        "TR_radius_reset_threshold" : float
+        "switch_to_optimality_requires_linearized_feasibility" : bool
+        "l1_constraint_violation_coefficient" : float
+        "barrier_initial_parameter" : float
+        "barrier_default_multiplier" : float
+        "barrier_tau_min" : float
+        "barrier_k_sigma" : float
+        "barrier_smax" : float
+        "barrier_k_mu" : float
+        "barrier_theta_mu" : float
+        "barrier_k_epsilon" : float
+        "barrier_update_fraction" : float
+        "barrier_regularization_exponent" : float
+        "barrier_small_direction_factor" : float
+        "barrier_push_variable_to_interior_k1" : float
+        "barrier_push_variable_to_interior_k2" : float
+        "barrier_damping_factor" : float
+        "least_square_multiplier_max_norm" : float
+        "BQPD_kmax_heuristic" : str
+        "QP_solver" : str
+        "LP_solver" : str
+        "linear_solver" : str
+
+
+    Returns
+    -------
+    res : OptimizeResult
+        The optimization result represented as a ""OptimizeResult"" object.
+        Important attributes are: ""x"" the solution array, ""success"" a
+        Boolean flag indicating if the optimizer exited successfully and
+        ""message"" which describes the cause of the termination. See
+        "OptimizeResult" for a description of other attributes.
+
+    Examples
+    --------
+    Let us consider the problem of minimizing the Rosenbrock function under constraints:
+    0.1 <= x**2 <= 0.8.
+
+    >>> res = minimize(rosen, x0 = np.array([1.3, 0.7, 0.8]), jac=rosen_der, method="filtersqp", tol=1e-3,
+    >>>                constraints=[NonlinearConstraint(fun=lambda x: x**2,  jac=lambda x: 2 * np.diag(x),
+    >>>                                                 lb=np.full(3, 0.1), ub=np.full(3, 0.8))],
+    >>>                options={"max_iterations": 10000})
+    >>> res.x
+
+    array([ 0.894427, 0.801979, 0.643170])
+
+    To define linear and non-linear constraints, please read the scipy.optimize. minimize documentation:
+    https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html
+
+
+    """
+    # Step 1: Input validation
+    x0 = np.atleast_1d(np.asarray(x0, dtype=float))
+    n = len(x0)
+    method_lower = method.lower()
+    if method_lower not in AVAILABLE_METHODS:
+        raise ValueError(
+            f"Unknown method '{method}'. Must be one of {AVAILABLE_METHODS}"
+        )
+    if options is None:
+        options = {}
+
+    # Step 2: Process bounds
+    if bounds is not None:
+        bounds = standardize_bounds(bounds, x0, "new")
+        _validate_bounds(bounds, x0, method)
+        var_lb = bounds.lb
+        var_ub = bounds.ub
+    else:
+        var_lb = [-Inf] * n
+        var_ub = [Inf] * n
+
+    # Step 3: Build objective callbacks
+    def objective_callback(x):
+        return float(fun(x, *args))
+
+    if jac is not None:
+        def gradient_callback(x, gradient):
+            gradient[:] = np.asarray(jac(x, *args), dtype=float).ravel()
+    else:
+        raise ValueError("Objective gradient is not defined.")
+
+    # Step 4: Normalize and merge constraints
+    # 4a: Normalize input to a list
+    if isinstance(constraints, dict) or isinstance(
+        constraints, (NonlinearConstraint, LinearConstraint)
+    ):
+        constraints = [constraints]
+    else:
+        constraints = list(constraints)
+
+    # 4b: Convert each constraint to (c_fun, c_jac, c_hess, c_lb, c_ub, m_i)
+    normalized = []
+    for i, con in enumerate(constraints):
+        if isinstance(con, NonlinearConstraint):
+            c_fun = con.fun
+            c_jac = con.jac
+            c_hess = getattr(con, "hess", None)
+            c_lb = np.atleast_1d(np.asarray(con.lb, dtype=float))
+            c_ub = np.atleast_1d(np.asarray(con.ub, dtype=float))
+            m_i = len(np.atleast_1d(c_fun(x0)))
+            normalized.append((c_fun, c_jac, c_hess, c_lb, c_ub, m_i))
+        elif isinstance(con, LinearConstraint):
+            A = np.atleast_2d(np.asarray(con.A, dtype=float))
+            m_i = A.shape[0]
+            c_lb = np.broadcast_to(
+                np.atleast_1d(np.asarray(con.lb, dtype=float)), (m_i,)
+            ).copy()
+            c_ub = np.broadcast_to(
+                np.atleast_1d(np.asarray(con.ub, dtype=float)), (m_i,)
+            ).copy()
+            # Use default-arg binding to capture A
+            c_fun = lambda x, _A=A: _A @ x
+            c_jac = lambda x, _A=A: _A
+            c_hess = None  # Linear constraints have zero Hessian
+            normalized.append((c_fun, c_jac, c_hess, c_lb, c_ub, m_i))
+        elif isinstance(con, dict):
+            c_fun_raw = con["fun"]
+            c_jac_raw = con.get("jac")
+            c_args = con.get("args", ())
+            c_type = con["type"]
+            # Evaluate to get dimension
+            val = np.atleast_1d(np.asarray(c_fun_raw(x0, *c_args), dtype=float))
+            m_i = len(val)
+            if c_type == "eq":
+                c_lb = np.zeros(m_i)
+                c_ub = np.zeros(m_i)
+            elif c_type == "ineq":
+                c_lb = np.zeros(m_i)
+                c_ub = np.full(m_i, Inf)
+            else:
+                raise ValueError(
+                    f"Unknown constraint type '{c_type}'. Must be 'eq' or 'ineq'"
+                )
+            # Use default-arg binding
+            c_fun = lambda x, _f=c_fun_raw, _a=c_args: np.atleast_1d(
+                np.asarray(_f(x, *_a), dtype=float)
+            )
+            if c_jac_raw is not None:
+                c_jac = lambda x, _j=c_jac_raw, _a=c_args: np.atleast_2d(
+                    np.asarray(_j(x, *_a), dtype=float)
+                )
+            else:
+                c_jac = None
+            c_hess = None  # Dict constraints don't support Hessians
+            normalized.append((c_fun, c_jac, c_hess, c_lb, c_ub, m_i))
+        else:
+            raise TypeError(f"Unsupported constraint type: {type(con)}")
+
+    # 4c: Merge
+    total_constraints = sum(item[5] for item in normalized)
+    if total_constraints > 0:
+        all_lb = np.concatenate([item[3] for item in normalized])
+        all_ub = np.concatenate([item[4] for item in normalized])
+    else:
+        all_lb = np.array([])
+        all_ub = np.array([])
+
+    # 4d: Build merged callbacks
+    def constraint_callback(x, constraint_values):
+        offset = 0
+        for c_fun_i, _, _, _, _, m_i in normalized:
+            constraint_values[offset:offset + m_i] = np.asarray(c_fun_i(x), dtype=float).ravel()
+            offset += m_i
+
+    def constraint_jacobian_callback(x, jacobian_values):
+        offset = 0
+        for c_fun_i, c_jac_i, _, _, _, m_i in normalized:
+            if c_jac_i is not None:
+                jac_block = np.atleast_2d(np.asarray(c_jac_i(x), dtype=float))
+            else:
+                raise ValueError("Constraint Jacobian is not provided.")
+            # Column-major (Fortran) ordering to match unopy convention
+            flat = jac_block.ravel(order="F")
+            for j in range(len(flat)):
+                jacobian_values[offset + j] = float(flat[j])
+            offset += len(flat)
+
+    # 4e: Dense sparsity pattern (column-major ordering)
+    if total_constraints > 0:
+        nnz_jac = total_constraints * n
+        # Column-major: iterate columns first, then rows within each column
+        col_indices = np.repeat(np.arange(n), total_constraints).tolist()
+        row_indices = np.tile(np.arange(total_constraints), n).tolist()
+    else:
+        nnz_jac = 0
+        row_indices = []
+        col_indices = []
+
+    # Step 5: Build unopy.Model
+    model = unopy.Model(
+        unopy.PROBLEM_NONLINEAR, n, var_lb, var_ub, unopy.ZERO_BASED_INDEXING
+    )
+    model.set_objective(unopy.MINIMIZE, objective_callback, gradient_callback)
+    if total_constraints > 0:
+        model.set_constraints(
+            total_constraints,
+            constraint_callback,
+            all_lb,
+            all_ub,
+            nnz_jac,
+            row_indices,
+            col_indices,
+            constraint_jacobian_callback,
+        )
+    model.set_initial_primal_iterate(list(x0))
+
+    # Step 5b: Lagrangian Hessian operator (if Hessian provided)
+    if hess is not None:
+        def hessian_operator_callback(x, evaluate_at_x, obj_mult,
+                                      multipliers, vector, result):
+            # Objective Hessian contribution
+            H = float(obj_mult) * np.atleast_2d(np.asarray(hess(x, *args), dtype=float))
+
+            # Constraint Hessian contributions
+            offset = 0
+            for _, _, c_hess_i, _, _, m_i in normalized:
+                if c_hess_i is not None:
+                    mult_slice = np.asarray(multipliers[offset:offset + m_i], dtype=float)
+                    H += np.atleast_2d(np.asarray(c_hess_i(x, mult_slice), dtype=float))
+                offset += m_i
+
+            # Compute H @ vector
+            result[:] = H @ vector
+
+        model.set_lagrangian_hessian_operator(hessian_operator_callback)
+        model.set_lagrangian_sign_convention(unopy.MULTIPLIER_POSITIVE)
+
+
+    # Step 6: Configure and run solver
+    uno_solver = unopy.UnoSolver()
+    uno_solver.set_preset(method_lower)
+    uno_solver.set_option("logger", "DEBUG3")
+    uno_solver.set_option("print_subproblem", True)
+    uno_solver.set_option("print_solution", True)
+    if tol is not None:
+        uno_solver.set_option("primal_tolerance", tol)
+        uno_solver.set_option("dual_tolerance", tol)
+    for key, value in options.items():
+        if key == "maxiter":
+            key = "max_iterations"
+        uno_solver.set_option(key, value)
+
+    result = uno_solver.optimize(model)
+
+    # Step 7: Return OptimizeResult
+    return OptimizeResult(
+        x=result.primal_solution,
+        success=str(result.optimization_status) == "OptimizationStatus.SUCCESS",
+        status=result.solution_status,
+        message=result.optimization_status,
+        fun=result.solution_objective,
+        nfev=result.number_objective_evaluations,
+        njev=result.number_jacobian_evaluations,
+        nhev=result.number_hessian_evaluations,
+        nit=result.number_iterations,
+        maxcv=result.solution_primal_feasibility,
+    )
diff --git a/interfaces/Python/scipy_interface/test_scipy_interface.py b/interfaces/Python/scipy_interface/test_scipy_interface.py
new file mode 100644
index 0000000000..b3c282a68f
--- /dev/null
+++ b/interfaces/Python/scipy_interface/test_scipy_interface.py
@@ -0,0 +1,101 @@
+# Copyright (c) 2026 Francois Gallard, Jean-Christophe Giret
+# Licensed under the MIT license. See LICENSE file in the project directory for details.
+
+import numpy as np
+import pytest
+from scipy.optimize import rosen, rosen_der, rosen_hess, NonlinearConstraint, LinearConstraint
+
+from scipy_interface.scipy_uno import minimize
+
+
+@pytest.mark.parametrize("method", ["filtersqp", "funnelsqp", "ipopt"])
+def test_rosen(method):
+    res = minimize(
+        rosen,
+        jac=rosen_der,
+        x0=np.zeros(2),
+        method=method,
+        options={"max_iterations": 10000},
+        tol=1e-4,
+    )
+    assert res.fun < 1e-4
+    assert res.success
+
+
+@pytest.mark.parametrize("method", ["filtersqp", "funnelsqp", "ipopt"])
+def test_rosen_bnds(method):
+    res = minimize(
+        rosen,
+        jac=rosen_der,
+        x0=np.zeros(2),
+        method=method,
+        bounds=[[-1.0, 0.5], [-1, 1]],
+        options={"max_iterations": 10000},
+    )
+    assert res.success
+
+
+@pytest.mark.parametrize(
+    "c_type", ["lambda", "NonLinearConstraint", "LinearConstraint"]
+)
+def test_rosen_constr(c_type):
+    if c_type == "NonLinearConstraint":
+        constr = NonlinearConstraint(
+            fun=lambda x: 1.0 - 2 * x,
+            jac=lambda x: -2 * np.eye(2),
+            lb=np.zeros(2),
+            ub=np.full(2, float("Inf")),
+        )
+    elif c_type == "LinearConstraint":
+        constr = LinearConstraint(A=np.eye(2), ub=0.5 * np.zeros(2), lb=-float("Inf"))
+    else:
+        constr = {
+            "type": "ineq",
+            "fun": lambda x: 1.0 - 2 * x,
+            "jac": lambda x: -2 * np.eye(2),
+        }
+    res = minimize(rosen, jac=rosen_der, x0=np.zeros(2), constraints=(constr,))
+
+    assert res.success
+
+@pytest.mark.parametrize("method", ["filtersqp", "funnelsqp", "ipopt"])
+def test_rosen_constr2(method):
+    res = minimize(
+        rosen,
+        np.array([1.3, 0.7, 0.8, 0.3]),
+        jac=rosen_der,
+        method=method,
+        tol=1e-3,
+        constraints=[
+            NonlinearConstraint(
+                fun=lambda x: x**2,
+                jac=lambda x: 2 * np.diag(x),
+                lb=np.full(4, 0.1),
+                ub=np.full(4, 0.8),
+            ),
+        ],
+        options={"max_iterations": 10000},
+    )
+    assert res.success
+
+@pytest.mark.parametrize("method", ["filtersqp", "funnelsqp"])
+def test_rosen_constr2_hess(method):
+    res = minimize(
+        rosen,
+        np.array([1.3, 0.7, 0.8]),
+        jac=rosen_der,
+        hess=rosen_hess,
+        method=method,
+        tol=1e-3,
+        constraints=[
+            NonlinearConstraint(
+                fun=lambda x: x**2,
+                jac=lambda x: 2 * np.diag(x),
+                hess=lambda x, v: 2 * np.diag(v),
+                lb=np.full(3, 0.1),
+                ub=np.full(3, 0.8),
+            ),
+        ],
+        options={"max_iterations": 10000},
+    )
+    assert res.success
diff --git a/pyproject.toml b/pyproject.toml
index 97d2647675..e76c08b0b8 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -31,7 +31,6 @@ Issues = "https://github.com/cvanaret/Uno/issues"
 cmake.source-dir = "."
 build.verbose = true
 wheel.cmake = true
-wheel.packages = ["interfaces/Python/unopy"]  # point to a dir containing __init__.py
 
 # support for passing paths to dependencies from environment to CMake
 [tool.scikit-build.cmake.define]
diff --git a/uno/optimization/EvaluationErrors.hpp b/uno/optimization/EvaluationErrors.hpp
index ee91d4bfa7..ca42c78f71 100644
--- a/uno/optimization/EvaluationErrors.hpp
+++ b/uno/optimization/EvaluationErrors.hpp
@@ -20,6 +20,12 @@ namespace uno {
          return "A numerical error was encountered while evaluating a gradient\n";
       }
    };
+   
+   struct JacobianEvaluationError : EvaluationError {
+      [[nodiscard]] const char* what() const noexcept override {
+         return "A numerical error was encountered while evaluating a Jacobian\n";
+      }
+   };
 
    struct HessianEvaluationError : EvaluationError {
       [[nodiscard]] const char* what() const noexcept override {
diff --git a/uno/optimization/Evaluations.cpp b/uno/optimization/Evaluations.cpp
index 6ea21567f8..1fbb951112 100644
--- a/uno/optimization/Evaluations.cpp
+++ b/uno/optimization/Evaluations.cpp
@@ -60,7 +60,7 @@ namespace uno {
          model.evaluate_jacobian(primals, this->jacobian_values.data());
          // check finiteness
          if (std::any_of(this->jacobian_values.begin(), this->jacobian_values.end(), invalid_value)) {
-            throw GradientEvaluationError();
+            throw JacobianEvaluationError();
          }
          this->is_jacobian_computed = true;
       }
@@ -102,4 +102,4 @@ namespace uno {
       this->is_objective_gradient_computed = false;
       this->is_jacobian_computed = false;
    }
-} // namespace
\ No newline at end of file
+} // namespace