diff --git a/CHANGELOG.md b/CHANGELOG.md
index 03ccc2a9..5c09ca90 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -13,12 +13,19 @@ Note: This is the main Changelog file for the Rust solver. The Changelog file fo
      --------------------- -->
 ## [v0.12.0] - Unreleased
 
+### Added
+
+- Richer Rust-side solver errors with human-readable reasons for projection failures, non-finite computations, linear algebra failures, and invalid internal solver states
+- Fallible constraint projections via `Constraint::project(...) -> FunctionCallResult`, with error propagation through FBS, PANOC, and ALM
+- Checked affine-space construction through `AffineSpace::try_new(...)` and `AffineSpaceError`
 
 ### Changed
 
 - Rust solver supports generic float types
-- Expanded Rust constraint test coverage with constructor validation, boundary/idempotence checks, and additional `BallP` / epigraph projection cases
-- Swap the cross-platform timer dependency to web-time, remove instant-specific wasm feature wiring, update optimizer timing call sites to use `web_time::Instant`, keep existing native and wasm timing behavior without stdweb risk
+- Expanded Rust constraint and solver test coverage with constructor validation, boundary/idempotence checks, additional `BallP` / epigraph projection cases, and broader `f32`/`f64` coverage
+- Swapped the cross-platform timer dependency to `web-time`, removed the old `instant`-specific wasm feature wiring, and updated optimizer timing call sites to use `web_time::Instant`
+- Improved Rust-side error handling across constraints and core solvers so projection failures and invalid numerical states are reported explicitly instead of being silently flattened
+- Refined `BallP`, `EpigraphSquaredNorm`, and related constraint implementations and docs for stronger numerical robustness and clearer behavior
 
 <!-- ---------------------
       v0.11.1
diff --git a/docs/openrust-basic.md b/docs/openrust-basic.md
index edaaaf13..4539fd63 100644
--- a/docs/openrust-basic.md
+++ b/docs/openrust-basic.md
@@ -208,7 +208,8 @@ fn main() {
 			.with_max_iter(max_iters);
 
 	// Invoke the solver
-	let status = panoc.solve(&mut u);
+	let status = panoc.solve(&mut u).unwrap();
+	assert!(status.has_converged());
 }
 ```
 
diff --git a/docs/python-tcp-ip.md b/docs/python-tcp-ip.md
index a8ee9d48..fdb5d3c3 100644
--- a/docs/python-tcp-ip.md
+++ b/docs/python-tcp-ip.md
@@ -131,24 +131,56 @@ $  echo '{ "Kill" : 1 }' | nc localhost 4598
 
 ### Error reporting
 
-In case a request cannot be processed, e.g., because the provided JSON is malformed, the provided vectors have incompatible dimensions, the TCP server will return to the client an error report. This is a JSON with three attributes: (i) a key-value pair `"type": "Error"`, to allow the client to tell that an error has occurred, (ii) a `code`, which can be used to uniquely identify the type of error and (iii) a `message`, which offers some human-readable details.
+If a request cannot be processed, for example because the JSON payload is
+malformed, the request body is not valid UTF-8, or one of the provided vectors
+has the wrong dimension, the TCP server returns an error report instead of a
+solution.
 
-For example, if the client provides an incompatible number of parameters, that is, if vector `parameter` is of the wrong length, then the server will return the following error:
+The error report is a JSON document with three fields:
+
+- `"type": "Error"` so the client can distinguish errors from successful solver
+  responses
+- `code`, a machine-readable integer code
+- `message`, a human-readable string with more context
+
+In particular, the `message` field is now intended to be descriptive. For
+dimension-related errors it includes the provided and expected lengths, and for
+solver failures it includes the underlying solver-side reason whenever
+available.
+
+For example, if the client provides an incompatible number of parameters, that
+is, if vector `parameter` has the wrong length, the server returns an error like
+the following:
+
+```json
+{
+  "type": "Error",
+  "code": 3003,
+  "message": "wrong number of parameters: provided 1, expected 2"
+}
+```
+
+Likewise, if the solver itself fails, the server returns code `2000` together
+with the propagated solver reason, for example:
 
 ```json
 {
-	"type": "Error", 
-	"code": 3003, 
-	"message": "wrong number of parameters"
+  "type": "Error",
+  "code": 2000,
+  "message": "problem solution failed: non-finite computation: gradient evaluation returned a non-finite value during an FBS step"
 }
 ```
 
-The following errors may be returned to the client
+The following error codes may be returned to the client:
 
 | Code      | Explanation                                 |
 |-----------|---------------------------------------------|
-| 1000      | Invalid request: Malformed or invalid JSON  |
+| 1000      | Invalid request: malformed JSON or invalid UTF-8 payload |
 | 1600      | Initial guess has incompatible dimensions   |
-| 1700      | Wrong dimension of Langrange multipliers    |
-| 2000      | Problem solution failed (solver error)      |
+| 1700      | Wrong dimension of Lagrange multipliers     |
+| 2000      | Problem solution failed; the message contains the solver-side reason |
 | 3003      | Vector `parameter` has wrong length         |
+
+When using the Python TCP client, these responses are surfaced as
+`opengen.tcp.solver_error.SolverError`, whose `code` and `message` properties
+mirror the JSON payload returned by the TCP server.
diff --git a/examples/panoc_ex1.rs b/examples/panoc_ex1.rs
index 01f3e497..c50d3301 100644
--- a/examples/panoc_ex1.rs
+++ b/examples/panoc_ex1.rs
@@ -29,7 +29,9 @@ fn main() {
     // define the cost function and its gradient
     let df = |u: &[f64], grad: &mut [f64]| -> Result<(), SolverError> {
         if a < 0.0 || b < 0.0 {
-            Err(SolverError::Cost)
+            Err(SolverError::Cost(
+                "Rosenbrock parameters must be nonnegative",
+            ))
         } else {
             rosenbrock_grad(a, b, u, grad);
             Ok(())
@@ -38,7 +40,9 @@ fn main() {
 
     let f = |u: &[f64], c: &mut f64| -> Result<(), SolverError> {
         if a < 0.0 || b < 0.0 {
-            Err(SolverError::Cost)
+            Err(SolverError::Cost(
+                "Rosenbrock parameters must be nonnegative",
+            ))
         } else {
             *c = rosenbrock_cost(a, b, u);
             Ok(())
diff --git a/open-codegen/CHANGELOG.md b/open-codegen/CHANGELOG.md
index e6940a39..87bcd1c5 100644
--- a/open-codegen/CHANGELOG.md
+++ b/open-codegen/CHANGELOG.md
@@ -14,10 +14,12 @@ Note: This is the Changelog file of `opengen` - the Python interface of OpEn
 
 - ROS2 package generation support via `BuildConfiguration.with_ros2(...)`, including auto-generated ROS2 templates, launcher, messages, and package wrapper code
 - Dedicated ROS2 tests covering package generation, build configuration behavior, rendered custom package settings, and end-to-end execution of a generated ROS2 node
+- More informative TCP solver error payloads, including clearer dimension/parameter validation failures and propagated solver-side failure messages
 
 ### Changed
 
 - Extended `RosConfiguration` so it can be used for both ROS and ROS2 package generation
+- Updated generated TCP server and C interface templates to work with the richer Rust solver error model and expose better failure information to clients
 
 
 ## [0.10.1] - 2026-03-25
diff --git a/open-codegen/opengen/tcp/solver_error.py b/open-codegen/opengen/tcp/solver_error.py
index 61ab7fc8..e2d2d72c 100644
--- a/open-codegen/opengen/tcp/solver_error.py
+++ b/open-codegen/opengen/tcp/solver_error.py
@@ -19,10 +19,10 @@ def code(self):
 
         Possible error codes are:
 
-        - **1000**: Invalid request: Malformed or invalid JSON
+        - **1000**: Invalid request: malformed JSON or invalid UTF-8 payload
         - **1600**: Initial guess has incomplete dimensions
         - **1700**: Wrong dimension of Lagrange multipliers
-        - **2000**: Problem solution failed (solver error)
+        - **2000**: Problem solution failed (message may include the solver reason)
         - **3003**: Parameter vector has wrong length
 
         :return: Error code
diff --git a/open-codegen/opengen/templates/c/optimizer_cinterface.rs.jinja b/open-codegen/opengen/templates/c/optimizer_cinterface.rs.jinja
index 039dfe0f..c27c356e 100644
--- a/open-codegen/opengen/templates/c/optimizer_cinterface.rs.jinja
+++ b/open-codegen/opengen/templates/c/optimizer_cinterface.rs.jinja
@@ -179,8 +179,11 @@ pub unsafe extern "C" fn {{meta.optimizer_name|lower}}_solve(
         },
         Err(e) => {{meta.optimizer_name}}SolverStatus {
             exit_status: match e {
-                SolverError::Cost => {{meta.optimizer_name}}ExitStatus::{{meta.optimizer_name}}NotConvergedCost,
-                SolverError::NotFiniteComputation => {{meta.optimizer_name}}ExitStatus::{{meta.optimizer_name}}NotConvergedNotFiniteComputation,
+                SolverError::Cost(_)
+                | SolverError::ProjectionFailed(_)
+                | SolverError::LinearAlgebraFailure(_)
+                | SolverError::InvalidProblemState(_) => {{meta.optimizer_name}}ExitStatus::{{meta.optimizer_name}}NotConvergedCost,
+                SolverError::NotFiniteComputation(_) => {{meta.optimizer_name}}ExitStatus::{{meta.optimizer_name}}NotConvergedNotFiniteComputation,
             },
             num_outer_iterations: u64::MAX as c_ulong,
             num_inner_iterations: u64::MAX as c_ulong,
@@ -209,4 +212,4 @@ pub unsafe extern "C" fn {{meta.optimizer_name|lower}}_free(instance: *mut {{met
     assert!(!instance.is_null());
     drop(Box::from_raw(instance));
 }
-{% endif %}
\ No newline at end of file
+{% endif %}
diff --git a/open-codegen/opengen/templates/tcp/tcp_server.rs b/open-codegen/opengen/templates/tcp/tcp_server.rs
index 042e176f..c676a7f5 100644
--- a/open-codegen/opengen/templates/tcp/tcp_server.rs
+++ b/open-codegen/opengen/templates/tcp/tcp_server.rs
@@ -9,7 +9,7 @@ extern crate clap;
 
 use std::{
     io::{prelude::Read, Write},
-    net::TcpListener,
+    net::{TcpListener, TcpStream},
 };
 
 use clap::{Arg, App};
@@ -80,27 +80,54 @@ struct OptimizerSolution<'a> {
     cost: f64,
 }
 
-fn pong(stream: &mut std::net::TcpStream, code: i32) {
-    let error_message = format!(
+fn write_bytes_to_stream(stream: &mut TcpStream, payload: &[u8], context: &str) -> bool {
+    if let Err(err) = stream.write_all(payload) {
+        warn!("{}: {}", context, err);
+        return false;
+    }
+    true
+}
+
+fn write_json_to_stream<T: Serialize>(
+    stream: &mut TcpStream,
+    payload: &T,
+    context: &str,
+) -> bool {
+    let payload_json = match serde_json::to_vec_pretty(payload) {
+        Ok(payload_json) => payload_json,
+        Err(err) => {
+            error!("{}: {}", context, err);
+            return false;
+        }
+    };
+    write_bytes_to_stream(stream, &payload_json, context)
+}
+
+fn pong(stream: &mut TcpStream, code: i32) {
+    let pong_message = format!(
         {% raw %}"{{\n\t\"Pong\" : {}\n}}\n"{% endraw %},
         code
     );
-    stream
-        .write_all(error_message.as_bytes())
-        .expect("cannot write to stream");
+    write_bytes_to_stream(stream, pong_message.as_bytes(), "could not write pong to stream");
 }
 
 /// Writes an error to the communication stream
-fn write_error_message(stream: &mut std::net::TcpStream, code: i32, error_msg: &str) {
-    let error_message = format!(
-        {% raw %}"{{\n\t\"type\" : \"Error\", \n\t\"code\" : {}, \n\t\"message\" : \"{}\"\n}}\n"{% endraw %},
+#[derive(Serialize)]
+struct ErrorResponse<'a> {
+    #[serde(rename = "type")]
+    response_type: &'a str,
+    code: i32,
+    message: &'a str,
+}
+
+fn write_error_message(stream: &mut TcpStream, code: i32, error_msg: &str) {
+    let error_response = ErrorResponse {
+        response_type: "Error",
         code,
-        error_msg
-    );
-    warn!("Invalid request {:?}", code);
-    stream
-        .write_all(error_message.as_bytes())
-        .expect("cannot write to stream");
+        message: error_msg,
+    };
+    warn!("TCP error {}: {}", code, error_msg);
+    write_json_to_stream(stream, &error_response, "could not write error response to stream");
 }
 
 /// Serializes the solution and solution status and returns it
@@ -108,7 +135,7 @@ fn write_error_message(stream: &mut std::net::TcpStream, code: i32, error_msg: &
 fn return_solution_to_client(
     status: AlmOptimizerStatus,
     solution: &[f64],
-    stream: &mut std::net::TcpStream,
+    stream: &mut TcpStream,
 ) {
     let empty_vec : [f64; 0] = Default::default();
     let solution: OptimizerSolution = OptimizerSolution {
@@ -125,10 +152,7 @@ fn return_solution_to_client(
         cost: status.cost(),
 
     };
-    let solution_json = serde_json::to_string_pretty(&solution).unwrap();
-    stream
-        .write_all(solution_json.as_bytes())
-        .expect("cannot write to stream");
+    write_json_to_stream(stream, &solution, "could not write optimizer solution to stream");
 }
 
 /// Handles an execution request
@@ -137,7 +161,7 @@ fn execution_handler(
     execution_parameter: &ExecutionParameter,
     u: &mut [f64],
     p: &mut [f64],
-    stream: &mut std::net::TcpStream,
+    stream: &mut TcpStream,
 ) {
     // ----------------------------------------------------
     // Set initial value
@@ -150,7 +174,12 @@ fn execution_handler(
         Some(u0) => {
             if u0.len() != {{meta.optimizer_name|upper}}_NUM_DECISION_VARIABLES {
                 warn!("initial guess has incompatible dimensions");
-                write_error_message(stream, 1600, "Initial guess has incompatible dimensions");
+                let error_message = format!(
+                    "initial guess has incompatible dimensions: provided {}, expected {}",
+                    u0.len(),
+                    {{meta.optimizer_name|upper}}_NUM_DECISION_VARIABLES
+                );
+                write_error_message(stream, 1600, &error_message);
                 return;
             }
             u.copy_from_slice(u0);
@@ -162,7 +191,12 @@ fn execution_handler(
     // ----------------------------------------------------
     if let Some(y0) = &execution_parameter.initial_lagrange_multipliers {
         if y0.len() != {{meta.optimizer_name|upper}}_N1 {
-            write_error_message(stream, 1700, "wrong dimension of Langrange multipliers");
+            let error_message = format!(
+                "wrong dimension of Langrange multipliers: provided {}, expected {}",
+                y0.len(),
+                {{meta.optimizer_name|upper}}_N1
+            );
+            write_error_message(stream, 1700, &error_message);
             return;
         }
     }
@@ -172,7 +206,12 @@ fn execution_handler(
     // ----------------------------------------------------
     let parameter = &execution_parameter.parameter;
     if parameter.len() != {{meta.optimizer_name|upper}}_NUM_PARAMETERS {
-        write_error_message(stream, 3003, "wrong number of parameters");
+        let error_message = format!(
+            "wrong number of parameters: provided {}, expected {}",
+            parameter.len(),
+            {{meta.optimizer_name|upper}}_NUM_PARAMETERS
+        );
+        write_error_message(stream, 3003, &error_message);
         return;
     }
     p.copy_from_slice(parameter);
@@ -185,8 +224,9 @@ fn execution_handler(
         Ok(ok_status) => {
             return_solution_to_client(ok_status, u, stream);
         }
-        Err(_) => {
-            write_error_message(stream, 2000, "Problem solution failed (solver error)");
+        Err(err) => {
+            let error_message = format!("problem solution failed: {}", err);
+            write_error_message(stream, 2000, &error_message);
         }
     }
 }
@@ -196,22 +236,52 @@ fn run_server(tcp_config: &TcpServerConfiguration) {
     let mut p = [0.0; {{meta.optimizer_name|upper}}_NUM_PARAMETERS];
     let mut cache = initialize_solver();
     info!("Done");
-    let listener = TcpListener::bind(format!("{}:{}", tcp_config.ip, tcp_config.port)).unwrap();
+    let listener = match TcpListener::bind(format!("{}:{}", tcp_config.ip, tcp_config.port)) {
+        Ok(listener) => listener,
+        Err(err) => {
+            error!(
+                "failed to bind TCP server at {}:{}: {}",
+                tcp_config.ip,
+                tcp_config.port,
+                err
+            );
+            return;
+        }
+    };
     let mut u = [0.0; {{meta.optimizer_name|upper}}_NUM_DECISION_VARIABLES];
     info!("listening started, ready to accept connections at {}:{}", tcp_config.ip, tcp_config.port);
-    for stream in listener.incoming() {
-        let mut stream = stream.unwrap();
+    'incoming: for stream in listener.incoming() {
+        let mut stream = match stream {
+            Ok(stream) => stream,
+            Err(err) => {
+                warn!("failed to accept incoming TCP connection: {}", err);
+                continue;
+            }
+        };
 
         //The following is more robust compared to `read_to_string`
         let mut bytes_buffer = vec![0u8; READ_BUFFER_SIZE];
         let mut read_data_length = 1;
         let mut buffer = String::new();
         while read_data_length != 0 {
-            read_data_length = stream
-                .read(&mut bytes_buffer)
-                .expect("could not read stream");
-            let new_string = String::from_utf8(bytes_buffer[0..read_data_length].to_vec())
-                .expect("sent data is not UFT-8");
+            read_data_length = match stream.read(&mut bytes_buffer) {
+                Ok(read_data_length) => read_data_length,
+                Err(err) => {
+                    warn!("could not read stream: {}", err);
+                    continue 'incoming;
+                }
+            };
+            let new_string = match String::from_utf8(bytes_buffer[0..read_data_length].to_vec()) {
+                Ok(new_string) => new_string,
+                Err(err) => {
+                    let error_message = format!(
+                        "invalid request: request body is not valid UTF-8 ({})",
+                        err.utf8_error()
+                    );
+                    write_error_message(&mut stream, 1000, &error_message);
+                    continue 'incoming;
+                }
+            };
             buffer.push_str(&new_string);
         }
 
@@ -236,8 +306,9 @@ fn run_server(tcp_config: &TcpServerConfiguration) {
                     pong(&mut stream, ping_code);
                 }
             },
-            Err(_) => {
-                write_error_message(&mut stream, 1000, "Invalid request");
+            Err(err) => {
+                let error_message = format!("invalid request: {}", err);
+                write_error_message(&mut stream, 1000, &error_message);
             }
         }
     }
diff --git a/open-codegen/test/test.py b/open-codegen/test/test.py
index d56f3075..dcebf68b 100644
--- a/open-codegen/test/test.py
+++ b/open-codegen/test/test.py
@@ -1,5 +1,7 @@
 import os
 import unittest
+import json
+import socket
 import casadi.casadi as cs
 import opengen as og
 import subprocess
@@ -227,6 +229,53 @@ def setUpHalfspace(cls):
                                                   cls.solverConfig())
         builder.build()
 
+    @classmethod
+    def setUpSolverError(cls):
+        u = cs.MX.sym("u", 1)
+        p = cs.MX.sym("p", 1)
+        phi = cs.dot(u, u)
+        bounds = og.constraints.Rectangle(xmin=[-1.0], xmax=[1.0])
+        tcp_config = og.config.TcpServerConfiguration(bind_port=3310)
+        meta = og.config.OptimizerMeta() \
+            .with_optimizer_name("solver_error")
+        problem = og.builder.Problem(u, p, phi) \
+            .with_constraints(bounds)
+        build_config = og.config.BuildConfiguration() \
+            .with_open_version(local_path=RustBuildTestCase.get_open_local_absolute_path()) \
+            .with_build_directory(RustBuildTestCase.TEST_DIR) \
+            .with_build_mode(og.config.BuildConfiguration.DEBUG_MODE) \
+            .with_tcp_interface_config(tcp_interface_config=tcp_config)
+        og.builder.OpEnOptimizerBuilder(problem,
+                                        metadata=meta,
+                                        build_configuration=build_config,
+                                        solver_configuration=cls.solverConfig()) \
+            .build()
+
+        target_lib = os.path.join(
+            RustBuildTestCase.TEST_DIR, "solver_error", "src", "lib.rs")
+        with open(target_lib, "r", encoding="utf-8") as fh:
+            solver_lib = fh.read()
+
+        # Look for this excerpt inside lib.rs (in the auto-generated solver)...
+        anchor = (
+            '    assert_eq!(u.len(), SOLVER_ERROR_NUM_DECISION_VARIABLES, '
+            '"Wrong number of decision variables (u)");\n'
+        )
+        # Replace the anchor with this so that if p[0] < 0, the function `solve`
+        # will reutrn an error of type SolverError::Cost
+        injected_guard = (
+            anchor +
+            '\n'
+            '    if p[0] < 0.0 {\n'
+            '        return Err(SolverError::Cost("forced solver error for TCP test"));\n'
+            '    }\n'
+        )
+        if anchor not in solver_lib:
+            raise RuntimeError("Could not inject deterministic solver error")
+
+        with open(target_lib, "w", encoding="utf-8") as fh:
+            fh.write(solver_lib.replace(anchor, injected_guard, 1))
+
     @classmethod
     def setUpClass(cls):
         cls.setUpPythonBindings()
@@ -237,6 +286,38 @@ def setUpClass(cls):
         cls.setUpPlain()
         cls.setUpOnlyParametricF2()
         cls.setUpHalfspace()
+        cls.setUpSolverError()
+
+    @staticmethod
+    def raw_tcp_request(ip, port, payload, buffer_size=4096):
+        with socket.socket(socket.AF_INET, socket.SOCK_STREAM, 0) as conn_socket:
+            conn_socket.connect((ip, port))
+            if isinstance(payload, str):
+                payload = payload.encode()
+            conn_socket.sendall(payload)
+            conn_socket.shutdown(socket.SHUT_WR)
+
+            data = b''
+            while True:
+                data_chunk = conn_socket.recv(buffer_size)
+                if not data_chunk:
+                    break
+                data += data_chunk
+
+        return json.loads(data.decode())
+
+    @staticmethod
+    def send_partial_tcp_payload_and_close(ip, port, payload):
+        with socket.socket(socket.AF_INET, socket.SOCK_STREAM, 0) as conn_socket:
+            conn_socket.connect((ip, port))
+            if isinstance(payload, str):
+                payload = payload.encode()
+            conn_socket.sendall(payload)
+
+    def start_tcp_manager(self, manager):
+        manager.start()
+        self.addCleanup(manager.kill) # at the end, kill the TCP mngr
+        return manager
 
     def test_python_bindings(self):
         import sys
@@ -314,22 +395,18 @@ def test_start_multiple_servers(self):
 
         # Start all servers
         for m in all_managers:
-            m.start()
+            self.start_tcp_manager(m)
 
         # Ping all
         for m in all_managers:
             m.ping()
 
-        # Kill all
-        for m in all_managers:
-            m.kill()
-
     def test_rust_build_only_f1(self):
         # Start the server using a custom bind IP and port
-        mng = og.tcp.OptimizerTcpManager(RustBuildTestCase.TEST_DIR + '/only_f1',
-                                         ip='0.0.0.0',
-                                         port=13757)
-        mng.start()
+        mng = self.start_tcp_manager(
+            og.tcp.OptimizerTcpManager(RustBuildTestCase.TEST_DIR + '/only_f1',
+                                       ip='0.0.0.0',
+                                       port=13757))
         pong = mng.ping()  # check if the server is alive
         self.assertEqual(1, pong["Pong"])
 
@@ -351,79 +428,89 @@ def test_rust_build_only_f1(self):
         self.assertFalse(response.is_ok())
         self.assertEqual(True, isinstance(status, og.tcp.SolverError))
         self.assertEqual(3003, status.code)
+        self.assertEqual(
+            "wrong number of parameters: provided 3, expected 2",
+            status.message)
 
         response = mng.call(p=[2.0, 10.0], initial_guess=[0.1, 0.2])
         self.assertFalse(response.is_ok())
         status = response.get()
         self.assertEqual(True, isinstance(status, og.tcp.SolverError))
         self.assertEqual(1600, status.code)
+        self.assertEqual(
+            "initial guess has incompatible dimensions: provided 2, expected 5",
+            status.message)
 
         response = mng.call(p=[2.0, 10.0], initial_y=[0.1])
         status = response.get()
         self.assertFalse(response.is_ok())
         self.assertEqual(True, isinstance(status, og.tcp.SolverError))
         self.assertEqual(1700, status.code)
-
-        mng.kill()
+        self.assertEqual(
+            "wrong dimension of Langrange multipliers: provided 1, expected 2",
+            status.message)
 
     def test_rust_build_only_f2_preconditioned(self):
-        mng1 = og.tcp.OptimizerTcpManager(
-            RustBuildTestCase.TEST_DIR + '/only_f2')
-        mng2 = og.tcp.OptimizerTcpManager(
-            RustBuildTestCase.TEST_DIR + '/only_f2_precond')
-        mng1.start()
-        mng2.start()
-
-        try:
-            response1 = mng1.call(p=[0.5, 8.5], initial_guess=[
-                                  1, 2, 3, 4, 0]).get()
-            response2 = mng2.call(p=[0.5, 8.5], initial_guess=[
-                                  1, 2, 3, 4, 0]).get()
-
-            self.assertEqual("Converged", response1.exit_status)
-            self.assertEqual("Converged", response2.exit_status)
-
-            # Further testing
-            slv_cfg = RustBuildTestCase.solverConfig()
-            # check that the solution is (near-) feasible
-            self.assertTrue(response1.f2_norm < slv_cfg.constraints_tolerance)
-            self.assertTrue(response2.f2_norm < slv_cfg.constraints_tolerance)
-            # check the nrom of the FPR
-            self.assertTrue(response1.last_problem_norm_fpr <
-                            slv_cfg.tolerance)
-            self.assertTrue(response2.last_problem_norm_fpr <
-                            slv_cfg.tolerance)
-            # compare the costs
-            self.assertAlmostEqual(response1.cost, response2.cost, 4)
-
-            x1, x2 = response1.solution, response2.solution
-            for i in range(len(x1)):
-                self.assertAlmostEqual(x1[i], x2[i], delta=5e-4)
-
-            response = mng1.call(p=[2.0, 10.0, 50.0])
-            self.assertFalse(response.is_ok())
-            status = response.get()
-            self.assertEqual(True, isinstance(status, og.tcp.SolverError))
-            self.assertEqual(3003, status.code)
-
-            response = mng1.call(p=[2.0, 10.0], initial_guess=[0.1, 0.2])
-            self.assertFalse(response.is_ok())
-            status = response.get()
-            self.assertEqual(True, isinstance(status, og.tcp.SolverError))
-            self.assertEqual(1600, status.code)
-
-            response = mng1.call(p=[2.0, 10.0], initial_y=[0.1])
-            self.assertFalse(response.is_ok())
-            status = response.get()
-            self.assertEqual(True, isinstance(status, og.tcp.SolverError))
-            self.assertEqual(1700, status.code)
-        finally:
-            mng1.kill()
-            mng2.kill()
+        mng1 = self.start_tcp_manager(og.tcp.OptimizerTcpManager(
+            RustBuildTestCase.TEST_DIR + '/only_f2'))
+        mng2 = self.start_tcp_manager(og.tcp.OptimizerTcpManager(
+            RustBuildTestCase.TEST_DIR + '/only_f2_precond'))
+
+        response1 = mng1.call(p=[0.5, 8.5], initial_guess=[
+                              1, 2, 3, 4, 0]).get()
+        response2 = mng2.call(p=[0.5, 8.5], initial_guess=[
+                              1, 2, 3, 4, 0]).get()
+
+        self.assertEqual("Converged", response1.exit_status)
+        self.assertEqual("Converged", response2.exit_status)
+
+        # Further testing
+        slv_cfg = RustBuildTestCase.solverConfig()
+        # check that the solution is (near-) feasible
+        self.assertTrue(response1.f2_norm < slv_cfg.constraints_tolerance)
+        self.assertTrue(response2.f2_norm < slv_cfg.constraints_tolerance)
+        # check the nrom of the FPR
+        self.assertTrue(response1.last_problem_norm_fpr <
+                        slv_cfg.tolerance)
+        self.assertTrue(response2.last_problem_norm_fpr <
+                        slv_cfg.tolerance)
+        # compare the costs
+        self.assertAlmostEqual(response1.cost, response2.cost, 4)
+
+        x1, x2 = response1.solution, response2.solution
+        for i in range(len(x1)):
+            self.assertAlmostEqual(x1[i], x2[i], delta=5e-4)
+
+        response = mng1.call(p=[2.0, 10.0, 50.0])
+        self.assertFalse(response.is_ok())
+        status = response.get()
+        self.assertEqual(True, isinstance(status, og.tcp.SolverError))
+        self.assertEqual(3003, status.code)
+        self.assertEqual(
+            "wrong number of parameters: provided 3, expected 2",
+            status.message)
+
+        response = mng1.call(p=[2.0, 10.0], initial_guess=[0.1, 0.2])
+        self.assertFalse(response.is_ok())
+        status = response.get()
+        self.assertEqual(True, isinstance(status, og.tcp.SolverError))
+        self.assertEqual(1600, status.code)
+        self.assertEqual(
+            "initial guess has incompatible dimensions: provided 2, expected 5",
+            status.message)
+
+        response = mng1.call(p=[2.0, 10.0], initial_y=[0.1])
+        self.assertFalse(response.is_ok())
+        status = response.get()
+        self.assertEqual(True, isinstance(status, og.tcp.SolverError))
+        self.assertEqual(1700, status.code)
+        self.assertEqual(
+            "wrong dimension of Langrange multipliers: provided 1, expected 0",
+            status.message)
 
     def test_rust_build_plain(self):
-        mng = og.tcp.OptimizerTcpManager(RustBuildTestCase.TEST_DIR + '/plain')
-        mng.start()
+        mng = self.start_tcp_manager(
+            og.tcp.OptimizerTcpManager(RustBuildTestCase.TEST_DIR + '/plain'))
         pong = mng.ping()  # check if the server is alive
         self.assertEqual(1, pong["Pong"])
 
@@ -433,13 +520,84 @@ def test_rust_build_plain(self):
         status = response.get()
         self.assertEqual("Converged", status.exit_status)
 
-        mng.kill()
+    def test_rust_build_plain_invalid_request_details(self):
+        self.start_tcp_manager(og.tcp.OptimizerTcpManager(
+            RustBuildTestCase.TEST_DIR + '/plain',
+            ip='127.0.0.1',
+            port=13758))
+
+        malformed_response = og.tcp.SolverResponse(
+            RustBuildTestCase.raw_tcp_request('127.0.0.1', 13758, '{"Run":'))
+        self.assertFalse(malformed_response.is_ok())
+        malformed_status = malformed_response.get()
+        self.assertEqual(1000, malformed_status.code)
+        self.assertTrue(
+            malformed_status.message.startswith("invalid request:"))
+        self.assertIn("line 1 column", malformed_status.message)
+
+        utf8_response = og.tcp.SolverResponse(
+            RustBuildTestCase.raw_tcp_request('127.0.0.1', 13758, b'\xff\xfe'))
+        self.assertFalse(utf8_response.is_ok())
+        utf8_status = utf8_response.get()
+        self.assertEqual(1000, utf8_status.code)
+        self.assertTrue(
+            utf8_status.message.startswith(
+                "invalid request: request body is not valid UTF-8"))
+
+    def test_rust_build_plain_survives_invalid_request(self):
+        mng = self.start_tcp_manager(og.tcp.OptimizerTcpManager(
+            RustBuildTestCase.TEST_DIR + '/plain',
+            ip='127.0.0.1',
+            port=13759))
+
+        malformed_response = og.tcp.SolverResponse(
+            RustBuildTestCase.raw_tcp_request('127.0.0.1', 13759, '{"Run":'))
+        self.assertFalse(malformed_response.is_ok())
+        malformed_status = malformed_response.get()
+        self.assertEqual(1000, malformed_status.code)
+
+        pong = mng.ping()
+        self.assertEqual(1, pong["Pong"])
+
+        response = mng.call(p=[2.0, 10.0])
+        self.assertTrue(response.is_ok())
+        self.assertEqual("Converged", response.get().exit_status)
+
+    def test_rust_build_plain_survives_disconnect_mid_request(self):
+        mng = self.start_tcp_manager(og.tcp.OptimizerTcpManager(
+            RustBuildTestCase.TEST_DIR + '/plain',
+            ip='127.0.0.1',
+            port=13760))
+
+        RustBuildTestCase.send_partial_tcp_payload_and_close(
+            '127.0.0.1',
+            13760,
+            '{"Run":{"parameter":[2.0'
+        )
+
+        pong = mng.ping()
+        self.assertEqual(1, pong["Pong"])
+
+        response = mng.call(p=[2.0, 10.0])
+        self.assertTrue(response.is_ok())
+        self.assertEqual("Converged", response.get().exit_status)
+
+    def test_rust_build_solver_error_details(self):
+        mng = self.start_tcp_manager(og.tcp.OptimizerTcpManager(
+            RustBuildTestCase.TEST_DIR + '/solver_error'))
+
+        response = mng.call(p=[-1.0])
+        self.assertFalse(response.is_ok())
+        status = response.get()
+        self.assertEqual(2000, status.code)
+        self.assertEqual(
+            "problem solution failed: cost or gradient evaluation failed: forced solver error for TCP test",
+            status.message)
 
     def test_rust_build_parametric_f2(self):
         # introduced to tackle issue #123
-        mng = og.tcp.OptimizerTcpManager(
-            RustBuildTestCase.TEST_DIR + '/parametric_f2')
-        mng.start()
+        mng = self.start_tcp_manager(og.tcp.OptimizerTcpManager(
+            RustBuildTestCase.TEST_DIR + '/parametric_f2'))
         pong = mng.ping()  # check if the server is alive
         self.assertEqual(1, pong["Pong"])
 
@@ -449,12 +607,10 @@ def test_rust_build_parametric_f2(self):
         status = response.get()
         self.assertEqual("Converged", status.exit_status)
         self.assertTrue(status.f2_norm < 1e-4)
-        mng.kill()
 
     def test_rust_build_parametric_halfspace(self):
-        mng = og.tcp.OptimizerTcpManager(
-            RustBuildTestCase.TEST_DIR + '/halfspace_optimizer')
-        mng.start()
+        mng = self.start_tcp_manager(og.tcp.OptimizerTcpManager(
+            RustBuildTestCase.TEST_DIR + '/halfspace_optimizer'))
         pong = mng.ping()  # check if the server is alive
         self.assertEqual(1, pong["Pong"])
 
@@ -470,8 +626,6 @@ def test_rust_build_parametric_halfspace(self):
         self.assertTrue(sum([u[i] * c[i] for i in range(5)]) - b <= eps)
         self.assertTrue(-sum([u[i] * c[i] for i in range(5)]) + b <= eps)
 
-        mng.kill()
-
     @staticmethod
     def c_bindings_helper(optimizer_name):
         p = subprocess.Popen(["/usr/bin/gcc",
@@ -515,6 +669,41 @@ def test_c_bindings(self):
         self.assertEqual(0, rc1)
         self.assertEqual(0, rc2)
 
+    def test_tcp_generated_server_builds(self):
+        tcp_iface_dir = os.path.join(
+            RustBuildTestCase.TEST_DIR, "plain", "tcp_iface_plain")
+        process = subprocess.Popen(
+            ["cargo", "build", "--quiet"],
+            cwd=tcp_iface_dir,
+            stdout=subprocess.PIPE,
+            stderr=subprocess.PIPE,
+        )
+        _stdout, stderr = process.communicate()
+
+        self.assertEqual(
+            0,
+            process.returncode,
+            msg=stderr.decode()
+        )
+
+    def test_tcp_manager_start_fails_cleanly_when_port_is_in_use(self):
+        mng1 = self.start_tcp_manager(og.tcp.OptimizerTcpManager(
+            RustBuildTestCase.TEST_DIR + '/plain',
+            ip='127.0.0.1',
+            port=13761))
+
+        mng2 = og.tcp.OptimizerTcpManager(
+            RustBuildTestCase.TEST_DIR + '/only_f1',
+            ip='127.0.0.1',
+            port=13761)
+        with self.assertRaises(Exception) as context:
+            mng2.start()
+
+        self.assertIn("Port 13761 not available", str(context.exception))
+
+        pong = mng1.ping()
+        self.assertEqual(1, pong["Pong"])
+
     def test_tcp_manager_remote_cannot_start(self):
         remote_tcp_manager = og.tcp.OptimizerTcpManager(
             ip='10.8.0.1', port=3345)
diff --git a/src/alm/alm_factory.rs b/src/alm/alm_factory.rs
index ffc55a4c..81987c8d 100644
--- a/src/alm/alm_factory.rs
+++ b/src/alm/alm_factory.rs
@@ -285,7 +285,7 @@ where
                 .zip(y_lagrange_mult.iter())
                 .for_each(|(ti, yi)| *ti += *yi / penalty_scale);
             s.copy_from_slice(&f1_u_plus_y_over_c);
-            set_c.project(&mut s);
+            set_c.project(&mut s)?;
             let dist_sq: T = matrix_operations::norm2_squared_diff(&f1_u_plus_y_over_c, &s);
             let scaling: T = cast::<T>(0.5) * penalty_parameter;
             *cost += scaling * dist_sq;
@@ -403,7 +403,7 @@ where
                 .zip(y_lagrange_mult.iter())
                 .for_each(|(ti, yi)| *ti += *yi / c_penalty_parameter);
             s_aux_var.copy_from_slice(&f1_u_plus_y_over_c); // s = t
-            set_c.project(&mut s_aux_var); // s = Proj_C(F1(u) + y/c)
+            set_c.project(&mut s_aux_var)?; // s = Proj_C(F1(u) + y/c)
 
             // t = F1(u) + y/c - Proj_C(F1(u) + y/c)
             f1_u_plus_y_over_c
@@ -519,7 +519,9 @@ mod tests {
         let f2 = mapping_f2;
         let jac_f2_tr =
             |_u: &[f64], _d: &[f64], _res: &mut [f64]| -> Result<(), crate::SolverError> {
-                Err(SolverError::NotFiniteComputation)
+                Err(SolverError::NotFiniteComputation(
+                    "mock Jacobian-transpose product returned a non-finite result",
+                ))
             };
         let factory = AlmFactory::new(
             mocks::f0,
diff --git a/src/alm/alm_optimizer.rs b/src/alm/alm_optimizer.rs
index 92c7dffa..e55ce7b3 100644
--- a/src/alm/alm_optimizer.rs
+++ b/src/alm/alm_optimizer.rs
@@ -674,7 +674,7 @@ where
                 .for_each(|((y_plus_i, y_i), w_alm_aux_i)| *y_plus_i = *w_alm_aux_i + *y_i / c);
 
             // Step #3: y_plus := Proj_C(y_plus)
-            alm_set_c.project(y_plus);
+            alm_set_c.project(y_plus)?;
 
             // Step #4
             y_plus
@@ -691,7 +691,7 @@ where
     }
 
     /// Project y on set Y
-    fn project_on_set_y(&mut self) {
+    fn project_on_set_y(&mut self) -> FunctionCallResult {
         let problem = &self.alm_problem;
         if let Some(y_set) = &problem.alm_set_y {
             // NOTE: as_mut() converts from &mut Option<T> to Option<&mut T>
@@ -700,9 +700,10 @@ where
             // *  which can be treated as  Option<&mut [T]>
             // * y_vec is                  &mut [T]
             if let Some(xi_vec) = self.alm_cache.xi.as_mut() {
-                y_set.project(&mut xi_vec[1..]);
+                y_set.project(&mut xi_vec[1..])?;
             }
         }
+        Ok(())
     }
 
     /// Solve inner problem
@@ -764,7 +765,7 @@ where
         inner_solver.solve(u)
     }
 
-    fn is_exit_criterion_satisfied(&self) -> bool {
+    fn is_exit_criterion_satisfied(&self) -> Result<bool, SolverError> {
         let cache = &self.alm_cache;
         let problem = &self.alm_problem;
         // Criterion 1: ||Delta y|| <= c * delta
@@ -788,8 +789,14 @@ where
         //              This should never happen because we set the AKKT tolerance
         //              in the constructor and can never become `None` again
         let criterion_3 =
-            cache.panoc_cache.akkt_tolerance.unwrap() <= self.epsilon_tolerance + T::epsilon();
-        criterion_1 && criterion_2 && criterion_3
+            cache
+                .panoc_cache
+                .akkt_tolerance
+                .ok_or(SolverError::InvalidProblemState(
+                    "missing inner AKKT tolerance while checking the exit criterion",
+                ))?
+                <= self.epsilon_tolerance + T::epsilon();
+        Ok(criterion_1 && criterion_2 && criterion_3)
     }
 
     /// Whether the penalty parameter should not be updated
@@ -826,10 +833,17 @@ where
         }
     }
 
-    fn update_inner_akkt_tolerance(&mut self) {
+    fn update_inner_akkt_tolerance(&mut self) -> FunctionCallResult {
         let cache = &mut self.alm_cache;
         // epsilon_{nu+1} := max(epsilon, beta*epsilon_nu)
-        let next_tolerance = cache.panoc_cache.akkt_tolerance.unwrap() * self.epsilon_update_factor;
+        let akkt_tolerance =
+            cache
+                .panoc_cache
+                .akkt_tolerance
+                .ok_or(SolverError::InvalidProblemState(
+                    "missing inner AKKT tolerance while updating it",
+                ))?;
+        let next_tolerance = akkt_tolerance * self.epsilon_update_factor;
         cache
             .panoc_cache
             .set_akkt_tolerance(if next_tolerance > self.epsilon_tolerance {
@@ -837,6 +851,7 @@ where
             } else {
                 self.epsilon_tolerance
             });
+        Ok(())
     }
 
     fn final_cache_update(&mut self) {
@@ -871,7 +886,7 @@ where
         let mut inner_exit_status: ExitStatus = ExitStatus::Converged;
 
         // Project y on Y
-        self.project_on_set_y();
+        self.project_on_set_y()?;
 
         // If the inner problem fails miserably, the failure should be propagated
         // upstream (using `?`). If the inner problem has not converged, that is fine,
@@ -895,7 +910,7 @@ where
         self.compute_alm_infeasibility()?; // ALM: ||y_plus - y||
 
         // Check exit criterion
-        if self.is_exit_criterion_satisfied() {
+        if self.is_exit_criterion_satisfied()? {
             // Do not continue the outer iteration
             // An (epsilon, delta)-AKKT point has been found
             return Ok(InnerProblemStatus::new(false, inner_exit_status));
@@ -904,7 +919,7 @@ where
         }
 
         // Update inner problem tolerance
-        self.update_inner_akkt_tolerance();
+        self.update_inner_akkt_tolerance()?;
 
         // conclusive step: updated iteration count, resets PANOC cache,
         // sets f2_norm = f2_norm_plus etc
@@ -1010,12 +1025,11 @@ where
             .with_penalty(c)
             .with_cost(cost);
         if self.alm_problem.n1 > 0 {
-            let status = status.with_lagrange_multipliers(
-                self.alm_cache
-                    .y_plus
-                    .as_ref()
-                    .expect("Although n1 > 0, there is no vector y (Lagrange multipliers)"),
-            );
+            let status = status.with_lagrange_multipliers(self.alm_cache.y_plus.as_ref().ok_or(
+                SolverError::InvalidProblemState(
+                    "missing Lagrange multipliers at the ALM solution",
+                ),
+            )?);
             Ok(status)
         } else {
             Ok(status)
@@ -1163,7 +1177,7 @@ mod tests {
             .with_initial_penalty(25.0)
             .with_initial_lagrange_multipliers(&[2., 3., 4., 10.]);
 
-        alm_optimizer.project_on_set_y();
+        alm_optimizer.project_on_set_y().unwrap();
         if let Some(xi_after_proj) = &alm_optimizer.alm_cache.xi {
             println!("xi = {:#?}", xi_after_proj);
             let y_projected_correct = [
@@ -1316,7 +1330,7 @@ mod tests {
             .with_initial_inner_tolerance(1e-1)
             .with_inner_tolerance_update_factor(0.2);
 
-        alm_optimizer.update_inner_akkt_tolerance();
+        alm_optimizer.update_inner_akkt_tolerance().unwrap();
 
         unit_test_utils::assert_nearly_equal(
             0.1,
@@ -1339,7 +1353,7 @@ mod tests {
         );
 
         for _i in 1..=5 {
-            alm_optimizer.update_inner_akkt_tolerance();
+            alm_optimizer.update_inner_akkt_tolerance().unwrap();
         }
         unit_test_utils::assert_nearly_equal(
             2e-5,
@@ -1445,20 +1459,20 @@ mod tests {
 
         // should not exit yet...
         assert!(
-            !alm_optimizer.is_exit_criterion_satisfied(),
+            !alm_optimizer.is_exit_criterion_satisfied().unwrap(),
             "exists right away"
         );
 
         let alm_optimizer = alm_optimizer
             .with_initial_inner_tolerance(1e-3)
             .with_epsilon_tolerance(1e-3);
-        assert!(!alm_optimizer.is_exit_criterion_satisfied());
+        assert!(!alm_optimizer.is_exit_criterion_satisfied().unwrap());
 
         alm_optimizer.alm_cache.delta_y_norm_plus = 1e-3;
-        assert!(!alm_optimizer.is_exit_criterion_satisfied());
+        assert!(!alm_optimizer.is_exit_criterion_satisfied().unwrap());
 
         alm_optimizer.alm_cache.f2_norm_plus = 1e-3;
-        assert!(alm_optimizer.is_exit_criterion_satisfied());
+        assert!(alm_optimizer.is_exit_criterion_satisfied().unwrap());
     }
 
     #[test]
diff --git a/src/alm/tests.rs b/src/alm/tests.rs
index 09b6e909..08e5ba8f 100644
--- a/src/alm/tests.rs
+++ b/src/alm/tests.rs
@@ -153,7 +153,9 @@ fn t_create_alm_optimizer() {
         .with_initial_lagrange_multipliers(&vec![5.0; n1]);
 
     let mut u = vec![0.0; nx];
-    println!("result = {:?}", alm_optimizer.solve(&mut u));
+    let result = alm_optimizer.solve(&mut u);
+    println!("result = {:?}", result);
+    assert!(result.is_ok());
 }
 
 #[test]
@@ -454,7 +456,7 @@ fn t_alm_numeric_test_1_f32() {
     assert!(mocks::mapping_f1_affine(&u, &mut f1u).is_ok());
     let mut projection = f1u.clone();
     let set_c_check = Ball2::new(None, 1.0_f32);
-    set_c_check.project(&mut projection);
+    set_c_check.project(&mut projection).unwrap();
     assert!((f1u[0] - projection[0]).abs() < 2e-3_f32);
     assert!((f1u[1] - projection[1]).abs() < 2e-3_f32);
 
@@ -618,7 +620,7 @@ fn t_alm_numeric_test_2_f32() {
     let mut f1u = vec![0.0_f32; n1];
     assert!(mocks::mapping_f1_affine(&u, &mut f1u).is_ok());
     let mut f1_proj = f1u.clone();
-    Ball2::new(None, 1.0_f32).project(&mut f1_proj);
+    Ball2::new(None, 1.0_f32).project(&mut f1_proj).unwrap();
     assert!((f1u[0] - f1_proj[0]).abs() < 2e-3_f32);
     assert!((f1u[1] - f1_proj[1]).abs() < 2e-3_f32);
 
diff --git a/src/constraints/affine_space.rs b/src/constraints/affine_space.rs
index bfa8c304..2bb2fec8 100644
--- a/src/constraints/affine_space.rs
+++ b/src/constraints/affine_space.rs
@@ -1,6 +1,7 @@
 use super::Constraint;
-use crate::matrix_operations;
-use crate::{CholeskyError, CholeskyFactorizer};
+use crate::{
+    matrix_operations, CholeskyError, CholeskyFactorizer, FunctionCallResult, SolverError,
+};
 
 use ndarray::{ArrayView1, ArrayView2, LinalgScalar};
 use num::Float;
@@ -82,7 +83,7 @@ where
     /// let affine_space = AffineSpace::try_new(a, b).unwrap();
     ///
     /// let mut x = [2.0, 2.0];
-    /// affine_space.project(&mut x);
+    /// affine_space.project(&mut x).unwrap();
     /// ```
     pub fn try_new(a: Vec<T>, b: Vec<T>) -> Result<Self, AffineSpaceError> {
         let n_rows = b.len();
@@ -144,30 +145,34 @@ where
     /// let b = vec![1., 2., -0.5];
     /// let affine_set = AffineSpace::new(a, b);
     /// let mut x = [1., -2., -0.3, 0.5];
-    /// affine_set.project(&mut x);
+    /// affine_set.project(&mut x).unwrap();
     /// ```
     ///
     /// The result is stored in `x` and it can be verified that $Ax = b$.
-    fn project(&self, x: &mut [T]) {
+    fn project(&self, x: &mut [T]) -> FunctionCallResult {
         let n = self.n_cols;
         assert!(x.len() == n, "x has wrong dimension");
 
         // Step 1: Compute e = Ax - b
-        let a = ArrayView2::from_shape((self.n_rows, self.n_cols), &self.a_mat)
-            .expect("invalid A shape");
+        let a = ArrayView2::from_shape((self.n_rows, self.n_cols), &self.a_mat).map_err(|_| {
+            SolverError::InvalidProblemState("failed to construct the affine-space matrix view")
+        })?;
         let x_view = ArrayView1::from(&x[..]);
         let b = ArrayView1::from(&self.b_vec[..]);
         let e = a.dot(&x_view) - b;
-        let e_slice: &[T] = e.as_slice().unwrap();
+        let e_slice: &[T] = e.as_slice().ok_or(SolverError::InvalidProblemState(
+            "affine-space residual vector is not stored contiguously",
+        ))?;
 
         // Step 2: Solve AA' z = e and compute z
-        let z = self.factorizer.solve(e_slice).unwrap();
+        let z = self.factorizer.solve(e_slice)?;
 
         // Step 3: Compute x = x - A'z
         let at_z = a.t().dot(&ArrayView1::from(&z[..]));
         for (xi, corr) in x.iter_mut().zip(at_z.iter()) {
             *xi = *xi - *corr;
         }
+        Ok(())
     }
 
     /// Affine sets are convex.
diff --git a/src/constraints/ball1.rs b/src/constraints/ball1.rs
index de2ee721..0039b268 100644
--- a/src/constraints/ball1.rs
+++ b/src/constraints/ball1.rs
@@ -1,5 +1,6 @@
 use super::Constraint;
 use super::Simplex;
+use crate::FunctionCallResult;
 use num::Float;
 use std::iter::Sum;
 
@@ -23,7 +24,7 @@ impl<'a, T: Float> Ball1<'a, T> {
     ///
     /// let ball = Ball1::new(None, 1.0);
     /// let mut x = [2.0, -0.5];
-    /// ball.project(&mut x);
+    /// ball.project(&mut x).unwrap();
     /// ```
     pub fn new(center: Option<&'a [T]>, radius: T) -> Self {
         assert!(radius > T::zero());
@@ -35,7 +36,7 @@ impl<'a, T: Float> Ball1<'a, T> {
         }
     }
 
-    fn project_on_ball1_centered_at_origin(&self, x: &mut [T])
+    fn project_on_ball1_centered_at_origin(&self, x: &mut [T]) -> FunctionCallResult
     where
         T: Sum<T>,
     {
@@ -46,11 +47,12 @@ impl<'a, T: Float> Ball1<'a, T> {
                 .zip(x.iter())
                 .for_each(|(ui, &xi)| *ui = xi.abs());
             // u = P_simplex(u)
-            self.simplex.project(&mut u);
+            self.simplex.project(&mut u)?;
             x.iter_mut()
                 .zip(u.iter())
                 .for_each(|(xi, &ui)| *xi = xi.signum() * ui);
         }
+        Ok(())
     }
 }
 
@@ -58,7 +60,7 @@ impl<'a, T> Constraint<T> for Ball1<'a, T>
 where
     T: Float + Sum<T>,
 {
-    fn project(&self, x: &mut [T]) {
+    fn project(&self, x: &mut [T]) -> FunctionCallResult {
         if let Some(center) = &self.center {
             assert_eq!(
                 x.len(),
@@ -68,13 +70,14 @@ where
             x.iter_mut()
                 .zip(center.iter())
                 .for_each(|(xi, &ci)| *xi = *xi - ci);
-            self.project_on_ball1_centered_at_origin(x);
+            self.project_on_ball1_centered_at_origin(x)?;
             x.iter_mut()
                 .zip(center.iter())
                 .for_each(|(xi, &ci)| *xi = *xi + ci);
         } else {
-            self.project_on_ball1_centered_at_origin(x);
+            self.project_on_ball1_centered_at_origin(x)?;
         }
+        Ok(())
     }
 
     fn is_convex(&self) -> bool {
diff --git a/src/constraints/ball2.rs b/src/constraints/ball2.rs
index 347654ad..1c0672c2 100644
--- a/src/constraints/ball2.rs
+++ b/src/constraints/ball2.rs
@@ -1,4 +1,5 @@
 use super::Constraint;
+use crate::FunctionCallResult;
 use num::Float;
 use std::iter::Sum;
 
@@ -21,7 +22,7 @@ impl<'a, T: Float> Ball2<'a, T> {
     ///
     /// let ball = Ball2::new(None, 1.0);
     /// let mut x = [2.0, 0.0];
-    /// ball.project(&mut x);
+    /// ball.project(&mut x).unwrap();
     /// ```
     pub fn new(center: Option<&'a [T]>, radius: T) -> Self {
         assert!(radius > T::zero());
@@ -34,7 +35,7 @@ impl<'a, T> Constraint<T> for Ball2<'a, T>
 where
     T: Float + Sum<T>,
 {
-    fn project(&self, x: &mut [T]) {
+    fn project(&self, x: &mut [T]) -> FunctionCallResult {
         if let Some(center) = &self.center {
             assert_eq!(
                 x.len(),
@@ -61,6 +62,7 @@ where
                 x.iter_mut().for_each(|x_| *x_ = *x_ / norm_over_radius);
             }
         }
+        Ok(())
     }
 
     fn is_convex(&self) -> bool {
diff --git a/src/constraints/ballinf.rs b/src/constraints/ballinf.rs
index 50b8fa37..6c4eb1f3 100644
--- a/src/constraints/ballinf.rs
+++ b/src/constraints/ballinf.rs
@@ -1,4 +1,5 @@
 use super::Constraint;
+use crate::FunctionCallResult;
 use num::Float;
 
 #[derive(Copy, Clone)]
@@ -22,7 +23,7 @@ impl<'a, T: Float> BallInf<'a, T> {
     ///
     /// let ball = BallInf::new(None, 1.0);
     /// let mut x = [2.0, -0.2, -3.0];
-    /// ball.project(&mut x);
+    /// ball.project(&mut x).unwrap();
     /// ```
     pub fn new(center: Option<&'a [T]>, radius: T) -> Self {
         assert!(radius > T::zero());
@@ -52,7 +53,7 @@ impl<'a, T: Float> Constraint<T> for BallInf<'a, T> {
     ///
     /// for all $i=1,\ldots, n$.
     ///
-    fn project(&self, x: &mut [T]) {
+    fn project(&self, x: &mut [T]) -> FunctionCallResult {
         if let Some(center) = &self.center {
             assert_eq!(
                 x.len(),
@@ -68,6 +69,7 @@ impl<'a, T: Float> Constraint<T> for BallInf<'a, T> {
                 .filter(|xi| xi.abs() > self.radius)
                 .for_each(|xi| *xi = xi.signum() * self.radius);
         }
+        Ok(())
     }
 
     fn is_convex(&self) -> bool {
diff --git a/src/constraints/ballp.rs b/src/constraints/ballp.rs
index 1d2930b4..5bd375c2 100644
--- a/src/constraints/ballp.rs
+++ b/src/constraints/ballp.rs
@@ -1,6 +1,7 @@
 use crate::numeric::cast;
 
 use super::Constraint;
+use crate::{FunctionCallResult, SolverError};
 use num::Float;
 
 #[derive(Copy, Clone)]
@@ -67,7 +68,7 @@ use num::Float;
 ///
 /// let ball = BallP::new(None, 1.0, 1.5, 1e-10, 100);
 /// let mut x = vec![3.0, -1.0, 2.0];
-/// ball.project(&mut x);
+/// ball.project(&mut x).unwrap();
 /// ```
 ///
 /// Project onto a translated \(\ell_p\)-ball:
@@ -78,7 +79,7 @@ use num::Float;
 /// let center = vec![1.0, 1.0, 1.0];
 /// let ball = BallP::new(Some(&center), 2.0, 3.0, 1e-10, 100);
 /// let mut x = vec![4.0, -1.0, 2.0];
-/// ball.project(&mut x);
+/// ball.project(&mut x).unwrap();
 /// ```
 ///
 /// # Notes
@@ -154,7 +155,7 @@ impl<'a, T: Float> BallP<'a, T> {
             .powf(T::one() / self.p)
     }
 
-    fn project_lp_ball(&self, x: &mut [T]) {
+    fn project_lp_ball(&self, x: &mut [T]) -> FunctionCallResult {
         let p = self.p;
         let r = self.radius;
         let tol = self.tolerance;
@@ -162,7 +163,7 @@ impl<'a, T: Float> BallP<'a, T> {
 
         let current_norm = self.lp_norm(x);
         if current_norm <= r {
-            return;
+            return Ok(());
         }
 
         let abs_x: Vec<T> = x.iter().map(|xi| xi.abs()).collect();
@@ -185,7 +186,9 @@ impl<'a, T: Float> BallP<'a, T> {
         while radius_error(lambda_hi) > T::zero() {
             lambda_hi = lambda_hi * cast::<T>(2.0);
             if lambda_hi > cast::<T>(1e20) {
-                panic!("Failed to bracket the Lagrange multiplier");
+                return Err(SolverError::ProjectionFailed(
+                    "failed to bracket the Lagrange multiplier",
+                ));
             }
         }
 
@@ -212,6 +215,7 @@ impl<'a, T: Float> BallP<'a, T> {
             let u = Self::solve_coordinate_newton(a, lambda_star, p, tol, max_iter);
             *xi = xi.signum() * u;
         });
+        Ok(())
     }
 
     /// Solve for u >= 0 the equation u + lambda * p * u^(p-1) = a
@@ -268,7 +272,7 @@ impl<'a, T: Float> BallP<'a, T> {
 }
 
 impl<'a, T: Float> Constraint<T> for BallP<'a, T> {
-    fn project(&self, x: &mut [T]) {
+    fn project(&self, x: &mut [T]) -> FunctionCallResult {
         if let Some(center) = &self.center {
             assert_eq!(
                 x.len(),
@@ -282,14 +286,15 @@ impl<'a, T: Float> Constraint<T> for BallP<'a, T> {
                 .zip(x.iter().zip(center.iter()))
                 .for_each(|(s, (xi, ci))| *s = *xi - *ci);
 
-            self.project_lp_ball(&mut shifted);
+            self.project_lp_ball(&mut shifted)?;
 
             x.iter_mut()
                 .zip(shifted.iter().zip(center.iter()))
                 .for_each(|(xi, (si, ci))| *xi = *ci + *si);
         } else {
-            self.project_lp_ball(x);
+            self.project_lp_ball(x)?;
         }
+        Ok(())
     }
 
     fn is_convex(&self) -> bool {
diff --git a/src/constraints/cartesian_product.rs b/src/constraints/cartesian_product.rs
index 2e7d12b0..df11e8ef 100644
--- a/src/constraints/cartesian_product.rs
+++ b/src/constraints/cartesian_product.rs
@@ -1,4 +1,5 @@
 use super::Constraint;
+use crate::FunctionCallResult;
 
 /// Cartesian product of constraints.
 ///
@@ -48,7 +49,7 @@ impl<'a, T> CartesianProduct<'a, T> {
     ///     .add_constraint(4, Ball2::new(None, 1.0));
     ///
     /// let mut x = [3.0, -2.0, 2.0, 0.0];
-    /// cartesian.project(&mut x);
+    /// cartesian.project(&mut x).unwrap();
     /// ```
     ///
     pub fn new() -> Self {
@@ -159,16 +160,14 @@ impl<'a, T> Constraint<T> for CartesianProduct<'a, T> {
     ///
     /// The method will panic if the dimension of `x` is not equal to the
     /// dimension of the Cartesian product (see `dimension()`)
-    fn project(&self, x: &mut [T]) {
+    fn project(&self, x: &mut [T]) -> FunctionCallResult {
         assert!(x.len() == self.dimension(), "x has wrong size");
         let mut j = 0;
-        self.idx
-            .iter()
-            .zip(self.constraints.iter())
-            .for_each(|(&i, c)| {
-                c.project(&mut x[j..i]);
-                j = i;
-            });
+        for (&i, c) in self.idx.iter().zip(self.constraints.iter()) {
+            c.project(&mut x[j..i])?;
+            j = i;
+        }
+        Ok(())
     }
 
     fn is_convex(&self) -> bool {
diff --git a/src/constraints/epigraph_squared_norm.rs b/src/constraints/epigraph_squared_norm.rs
index d4ff435b..9242cb34 100644
--- a/src/constraints/epigraph_squared_norm.rs
+++ b/src/constraints/epigraph_squared_norm.rs
@@ -1,4 +1,5 @@
 use crate::{matrix_operations, numeric::cast};
+use crate::{FunctionCallResult, SolverError};
 
 use super::Constraint;
 use num::Float;
@@ -27,7 +28,7 @@ impl EpigraphSquaredNorm {
     ///
     /// let epi = EpigraphSquaredNorm::new();
     /// let mut x = [1.0, 2.0, 1.0];
-    /// epi.project(&mut x);
+    /// epi.project(&mut x).unwrap();
     /// ```
     #[must_use]
     pub fn new() -> Self {
@@ -62,7 +63,10 @@ where
     ///
     /// Panics if:
     ///
-    /// - `x.len() < 2`,
+    /// - `x.len() < 2`.
+    ///
+    /// Returns an error if:
+    ///
     /// - no admissible real root is found,
     /// - the Newton derivative becomes too small,
     /// - the final scaling factor is numerically singular.
@@ -77,9 +81,9 @@ where
     /// // Here, z = [1., 2., 3.] and t = 4.
     /// let mut x = [1., 2., 3., 4.];
     ///
-    /// epi.project(&mut x);
+    /// epi.project(&mut x).unwrap();
     /// ```
-    fn project(&self, x: &mut [T]) {
+    fn project(&self, x: &mut [T]) -> FunctionCallResult {
         assert!(
             x.len() >= 2,
             "EpigraphSquaredNorm::project requires x.len() >= 2"
@@ -92,7 +96,7 @@ where
 
         // Already feasible
         if norm_z_sq <= t {
-            return;
+            return Ok(());
         }
 
         // Cubic:
@@ -122,8 +126,9 @@ where
             }
         }
 
-        let mut zsol =
-            right_root.expect("EpigraphSquaredNorm::project: no admissible real root found");
+        let mut zsol = right_root.ok_or(SolverError::ProjectionFailed(
+            "no admissible real root found for the cubic projection equation",
+        ))?;
 
         // Newton refinement
         let newton_max_iters: usize = 5;
@@ -139,10 +144,11 @@ where
             }
 
             let dp_z = cast::<T>(3.0) * a3 * zsol_sq + cast::<T>(2.0) * a2 * zsol + a1;
-            assert!(
-                num::Float::abs(dp_z) > cast::<T>(10.0) * T::epsilon(),
-                "EpigraphSquaredNorm::project: Newton derivative too small"
-            );
+            if num::Float::abs(dp_z) <= cast::<T>(1e-15) {
+                return Err(SolverError::ProjectionFailed(
+                    "Newton refinement derivative is too small",
+                ));
+            }
 
             zsol = zsol - p_z / dp_z;
         }
@@ -150,16 +156,18 @@ where
         let right_root = zsol;
         let scaling = cast::<T>(1.0) + cast::<T>(2.0) * (right_root - t);
 
-        assert!(
-            num::Float::abs(scaling) > cast::<T>(10.0) * T::epsilon(),
-            "EpigraphSquaredNorm::project: scaling factor too small"
-        );
+        if num::Float::abs(scaling) <= cast::<T>(1e-15) {
+            return Err(SolverError::ProjectionFailed(
+                "projection scaling factor is numerically singular",
+            ));
+        }
 
         // Projection
         for xi in x.iter_mut().take(nx) {
             *xi = *xi / scaling;
         }
         x[nx] = right_root;
+        Ok(())
     }
 
     /// This is a convex set, so this function returns `true`.
diff --git a/src/constraints/finite.rs b/src/constraints/finite.rs
index 2c76ede8..30ab9c48 100644
--- a/src/constraints/finite.rs
+++ b/src/constraints/finite.rs
@@ -1,4 +1,5 @@
 use super::Constraint;
+use crate::FunctionCallResult;
 use num::Float;
 use std::iter::Sum;
 
@@ -92,7 +93,7 @@ where
     /// ];
     /// let finite_set = FiniteSet::new(data);
     /// let mut x = [0.7, 0.6];
-    /// finite_set.project(&mut x); // compute projection
+    /// finite_set.project(&mut x).unwrap(); // compute projection
     /// ```
     ///
     /// # Panics
@@ -100,7 +101,7 @@ where
     /// This method panics if the dimension of `x` is not equal to the
     /// dimension of the points in the finite set.
     ///
-    fn project(&self, x: &mut [T]) {
+    fn project(&self, x: &mut [T]) -> FunctionCallResult {
         assert_eq!(x.len(), self.data[0].len(), "x has incompatible dimension");
         let mut idx: usize = 0;
         let mut best_distance = T::infinity();
@@ -112,6 +113,7 @@ where
             }
         }
         x.copy_from_slice(self.data[idx]);
+        Ok(())
     }
 
     fn is_convex(&self) -> bool {
diff --git a/src/constraints/halfspace.rs b/src/constraints/halfspace.rs
index 40ed4104..cd039128 100644
--- a/src/constraints/halfspace.rs
+++ b/src/constraints/halfspace.rs
@@ -1,5 +1,6 @@
 use super::Constraint;
 use crate::matrix_operations;
+use crate::FunctionCallResult;
 use num::Float;
 use std::iter::Sum;
 
@@ -47,7 +48,7 @@ where
     /// let offset = 1.0;
     /// let halfspace = Halfspace::new(&normal_vector, offset);
     /// let mut x = [-1., 3.];
-    /// halfspace.project(&mut x);
+    /// halfspace.project(&mut x).unwrap();
     /// ```
     ///
     pub fn new(normal_vector: &'a [T], offset: T) -> Self {
@@ -87,7 +88,7 @@ where
     /// This method panics if the length of `x` is not equal to the dimension
     /// of the halfspace.
     ///
-    fn project(&self, x: &mut [T]) {
+    fn project(&self, x: &mut [T]) -> FunctionCallResult {
         let inner_product = matrix_operations::inner_product(x, self.normal_vector);
         if inner_product > self.offset {
             let factor = (inner_product - self.offset) / self.normal_vector_squared_norm;
@@ -95,6 +96,7 @@ where
                 .zip(self.normal_vector.iter())
                 .for_each(|(x, normal_vector_i)| *x = *x - factor * *normal_vector_i);
         }
+        Ok(())
     }
 
     /// Halfspaces are convex sets
diff --git a/src/constraints/hyperplane.rs b/src/constraints/hyperplane.rs
index 1400fbad..958223c6 100644
--- a/src/constraints/hyperplane.rs
+++ b/src/constraints/hyperplane.rs
@@ -1,5 +1,6 @@
 use super::Constraint;
 use crate::matrix_operations;
+use crate::FunctionCallResult;
 use num::Float;
 use std::iter::Sum;
 
@@ -46,7 +47,7 @@ where
     /// let offset = 1.0;
     /// let hyperplane = Hyperplane::new(&normal_vector, offset);
     /// let mut x = [-1., 3.];
-    /// hyperplane.project(&mut x);
+    /// hyperplane.project(&mut x).unwrap();
     /// ```
     ///
     pub fn new(normal_vector: &'a [T], offset: T) -> Self {
@@ -87,13 +88,14 @@ where
     /// This method panics if the length of `x` is not equal to the dimension
     /// of the hyperplane.
     ///
-    fn project(&self, x: &mut [T]) {
+    fn project(&self, x: &mut [T]) -> FunctionCallResult {
         assert_eq!(x.len(), self.normal_vector.len(), "x has wrong dimension");
         let inner_product = matrix_operations::inner_product(x, self.normal_vector);
         let factor = (inner_product - self.offset) / self.normal_vector_squared_norm;
         x.iter_mut()
             .zip(self.normal_vector.iter())
             .for_each(|(x, nrm_vct)| *x = *x - factor * *nrm_vct);
+        Ok(())
     }
 
     /// Hyperplanes are convex sets
diff --git a/src/constraints/mod.rs b/src/constraints/mod.rs
index aed5a22e..478d54e0 100644
--- a/src/constraints/mod.rs
+++ b/src/constraints/mod.rs
@@ -8,6 +8,8 @@
 //!
 //! [`Constraint`]: trait.Constraint.html
 
+use crate::FunctionCallResult;
+
 mod affine_space;
 mod ball1;
 mod ball2;
@@ -57,7 +59,17 @@ pub trait Constraint<T = f64> {
     ///
     /// - `x`: The given vector $x$ is updated with the projection on the set
     ///
-    fn project(&self, x: &mut [T]);
+    /// ## Errors
+    ///
+    /// Implementations should return `Err(...)` when the projection cannot be
+    /// computed reliably because of a numerical or internal failure.
+    ///
+    /// ## Panics
+    ///
+    /// Implementations may still panic on clear API misuse, such as calling the
+    /// projection with a slice of incompatible dimension.
+    ///
+    fn project(&self, x: &mut [T]) -> FunctionCallResult;
 
     /// Returns true if and only if the set is convex
     fn is_convex(&self) -> bool;
diff --git a/src/constraints/no_constraints.rs b/src/constraints/no_constraints.rs
index 1d156db4..c04447b3 100644
--- a/src/constraints/no_constraints.rs
+++ b/src/constraints/no_constraints.rs
@@ -1,4 +1,5 @@
 use super::Constraint;
+use crate::FunctionCallResult;
 
 /// The whole space, no constraints
 #[derive(Default, Clone, Copy)]
@@ -14,7 +15,7 @@ impl NoConstraints {
     ///
     /// let no_constraints = NoConstraints::new();
     /// let mut x = [1.0, -2.0, 3.0];
-    /// no_constraints.project(&mut x);
+    /// no_constraints.project(&mut x).unwrap();
     /// ```
     ///
     #[must_use]
@@ -24,7 +25,9 @@ impl NoConstraints {
 }
 
 impl<T> Constraint<T> for NoConstraints {
-    fn project(&self, _x: &mut [T]) {}
+    fn project(&self, _x: &mut [T]) -> FunctionCallResult {
+        Ok(())
+    }
 
     fn is_convex(&self) -> bool {
         true
diff --git a/src/constraints/rectangle.rs b/src/constraints/rectangle.rs
index d4df786f..e81caabf 100644
--- a/src/constraints/rectangle.rs
+++ b/src/constraints/rectangle.rs
@@ -1,4 +1,5 @@
 use super::Constraint;
+use crate::FunctionCallResult;
 use num::Float;
 
 #[derive(Clone, Copy)]
@@ -44,15 +45,18 @@ impl<'a, T: Float> Rectangle<'a, T> {
     /// let xmax = [1.0, 2.0];
     /// let rectangle = Rectangle::new(Some(&xmin), Some(&xmax));
     /// let mut x = [3.0, -4.0];
-    /// rectangle.project(&mut x);
+    /// rectangle.project(&mut x).unwrap();
     /// ```
     ///
     pub fn new(xmin: Option<&'a [T]>, xmax: Option<&'a [T]>) -> Self {
         assert!(xmin.is_some() || xmax.is_some()); // xmin or xmax must be Some
-        assert!(
-            xmin.is_none() || xmax.is_none() || xmin.unwrap().len() == xmax.unwrap().len(),
-            "incompatible dimensions of xmin and xmax"
-        );
+        if let (Some(xmin), Some(xmax)) = (xmin, xmax) {
+            assert_eq!(
+                xmin.len(),
+                xmax.len(),
+                "incompatible dimensions of xmin and xmax"
+            );
+        }
         if let (Some(xmin), Some(xmax)) = (xmin, xmax) {
             assert!(
                 xmin.iter()
@@ -67,7 +71,7 @@ impl<'a, T: Float> Rectangle<'a, T> {
 }
 
 impl<'a, T: Float> Constraint<T> for Rectangle<'a, T> {
-    fn project(&self, x: &mut [T]) {
+    fn project(&self, x: &mut [T]) -> FunctionCallResult {
         if let Some(xmin) = &self.xmin {
             assert_eq!(
                 x.len(),
@@ -93,6 +97,7 @@ impl<'a, T: Float> Constraint<T> for Rectangle<'a, T> {
                 };
             });
         }
+        Ok(())
     }
 
     fn is_convex(&self) -> bool {
diff --git a/src/constraints/simplex.rs b/src/constraints/simplex.rs
index f0e34cb5..da1bf1bf 100644
--- a/src/constraints/simplex.rs
+++ b/src/constraints/simplex.rs
@@ -1,6 +1,7 @@
 use crate::numeric::cast;
 
 use super::Constraint;
+use crate::FunctionCallResult;
 use num::Float;
 
 #[derive(Copy, Clone)]
@@ -23,7 +24,7 @@ impl<T: Float> Simplex<T> {
     ///
     /// let simplex = Simplex::new(1.0);
     /// let mut x = [0.5, -0.5, 2.0];
-    /// simplex.project(&mut x);
+    /// simplex.project(&mut x).unwrap();
     /// ```
     pub fn new(alpha: T) -> Self {
         assert!(alpha > T::zero(), "alpha is nonpositive");
@@ -37,7 +38,7 @@ impl<T: Float> Constraint<T> for Simplex<T> {
     /// See: Laurent Condat. Fast Projection onto the Simplex and the $\ell_1$ Ball.
     /// <em>Mathematical Programming, Series A,</em> Springer, 2016, 158 (1), pp.575-585.
     /// ⟨<a href="https://dx.doi.org/10.1007/s10107-015-0946-6">10.1007/s10107-015-0946-6</a>⟩.
-    fn project(&self, x: &mut [T]) {
+    fn project(&self, x: &mut [T]) -> FunctionCallResult {
         assert!(!x.is_empty(), "x must be nonempty");
         let a = &self.alpha;
 
@@ -95,6 +96,7 @@ impl<T: Float> Constraint<T> for Simplex<T> {
         // ---- step 6
         let zero = T::zero();
         x.iter_mut().for_each(|x_n| *x_n = zero.max(*x_n - rho));
+        Ok(())
     }
 
     fn is_convex(&self) -> bool {
diff --git a/src/constraints/soc.rs b/src/constraints/soc.rs
index 0f12b40b..7b0de112 100644
--- a/src/constraints/soc.rs
+++ b/src/constraints/soc.rs
@@ -1,5 +1,6 @@
 use super::Constraint;
 use crate::matrix_operations;
+use crate::FunctionCallResult;
 use num::Float;
 use std::iter::Sum;
 
@@ -48,7 +49,7 @@ impl<T: Float> SecondOrderCone<T> {
     ///
     /// let cone = SecondOrderCone::new(1.0);
     /// let mut x = [2.0, 0.0, 0.5];
-    /// cone.project(&mut x);
+    /// cone.project(&mut x).unwrap();
     /// ```
     pub fn new(alpha: T) -> SecondOrderCone<T> {
         assert!(alpha > T::zero()); // alpha must be positive
@@ -71,7 +72,7 @@ where
     ///
     /// This method panics if the length of `x` is less than 2.
     ///
-    fn project(&self, x: &mut [T]) {
+    fn project(&self, x: &mut [T]) -> FunctionCallResult {
         // x = (z, r)
         let n = x.len();
         assert!(n >= 2, "x must be of dimension at least 2");
@@ -87,6 +88,7 @@ where
                 .for_each(|v| *v = *v * self.alpha * beta / norm_z);
             x[n - 1] = beta;
         }
+        Ok(())
     }
 
     fn is_convex(&self) -> bool {
diff --git a/src/constraints/sphere2.rs b/src/constraints/sphere2.rs
index 11a7f57f..9a702872 100644
--- a/src/constraints/sphere2.rs
+++ b/src/constraints/sphere2.rs
@@ -1,6 +1,7 @@
 use crate::numeric::cast;
 
 use super::Constraint;
+use crate::FunctionCallResult;
 use num::Float;
 use std::iter::Sum;
 
@@ -30,7 +31,7 @@ impl<'a, T: Float> Sphere2<'a, T> {
     ///
     /// let sphere = Sphere2::new(None, 1.0);
     /// let mut x = [3.0, 4.0];
-    /// sphere.project(&mut x);
+    /// sphere.project(&mut x).unwrap();
     /// ```
     pub fn new(center: Option<&'a [T]>, radius: T) -> Self {
         assert!(radius > T::zero());
@@ -62,7 +63,7 @@ where
     /// Panics if `x` is empty or, when a center is provided, if `x` and
     /// `center` have incompatible dimensions.
     ///
-    fn project(&self, x: &mut [T]) {
+    fn project(&self, x: &mut [T]) -> FunctionCallResult {
         let epsilon = cast::<T>(1e-12);
         assert!(!x.is_empty(), "x must be nonempty");
         if let Some(center) = &self.center {
@@ -75,7 +76,7 @@ where
             if norm_difference <= epsilon {
                 x.copy_from_slice(center);
                 x[0] = x[0] + self.radius;
-                return;
+                return Ok(());
             }
             x.iter_mut().zip(center.iter()).for_each(|(x, c)| {
                 *x = *c + self.radius * (*x - *c) / norm_difference;
@@ -84,11 +85,12 @@ where
             let norm_x = crate::matrix_operations::norm2(x);
             if norm_x <= epsilon {
                 x[0] = x[0] + self.radius;
-                return;
+                return Ok(());
             }
             let norm_over_radius = self.radius / norm_x;
             x.iter_mut().for_each(|x_| *x_ = *x_ * norm_over_radius);
         }
+        Ok(())
     }
 
     /// Returns false (the sphere is not a convex set)
diff --git a/src/constraints/tests.rs b/src/constraints/tests.rs
index 55c97bee..328c4b52 100644
--- a/src/constraints/tests.rs
+++ b/src/constraints/tests.rs
@@ -11,7 +11,7 @@ fn t_zero_set() {
     let zero = Zero::new();
     let mut x = [1.0, 2.0, 3.0];
     let x_projection = [0.0; 3];
-    zero.project(&mut x);
+    zero.project(&mut x).unwrap();
     unit_test_utils::assert_nearly_equal_array(
         &x_projection,
         &x,
@@ -25,7 +25,7 @@ fn t_zero_set() {
 fn t_zero_set_f32() {
     let zero = Zero::new();
     let mut x = [1.0_f32, -2.0, 3.5];
-    zero.project(&mut x);
+    zero.project(&mut x).unwrap();
     assert_eq!([0.0_f32, 0.0, 0.0], x);
 }
 
@@ -40,7 +40,7 @@ fn t_hyperplane() {
         0.285_714_285_714_286,
         0.928_571_428_571_429,
     ];
-    hyperplane.project(&mut x);
+    hyperplane.project(&mut x).unwrap();
     unit_test_utils::assert_nearly_equal_array(
         &x,
         &x_proj_expected,
@@ -63,7 +63,7 @@ fn t_hyperplane_wrong_dimension() {
     let normal_vector = [1.0, 2.0, 3.0];
     let hyperplane = Hyperplane::new(&normal_vector, 1.0);
     let mut x = [1.0, 2.0];
-    hyperplane.project(&mut x);
+    hyperplane.project(&mut x).unwrap();
 }
 
 #[test]
@@ -73,7 +73,7 @@ fn t_halfspace_project_inside() {
     let halfspace = Halfspace::new(&normal_vector, offset);
     let mut x = [-1., 3.];
     let x_expected = [-1., 3.];
-    halfspace.project(&mut x);
+    halfspace.project(&mut x).unwrap();
     unit_test_utils::assert_nearly_equal_array(
         &x,
         &x_expected,
@@ -90,7 +90,7 @@ fn t_halfspace_project_outside() {
     let halfspace = Halfspace::new(&normal_vector, offset);
     let mut x = [-1., 3.];
     let x_expected = [-1.8, 1.4];
-    halfspace.project(&mut x);
+    halfspace.project(&mut x).unwrap();
     unit_test_utils::assert_nearly_equal_array(
         &x,
         &x_expected,
@@ -121,7 +121,7 @@ fn t_finite_set() {
     let data: &[&[f64]] = &[&[0.0, 0.0], &[1.0, 1.0], &[0.0, 1.0], &[1.0, 0.0]];
     let finite_set = FiniteSet::new(data);
     let mut x = [0.7, 0.6];
-    finite_set.project(&mut x);
+    finite_set.project(&mut x).unwrap();
     unit_test_utils::assert_nearly_equal_array(
         &[1.0, 1.0],
         &x,
@@ -130,7 +130,7 @@ fn t_finite_set() {
         "projection is wrong (should be [1,1])",
     );
     x = [-0.1, 0.2];
-    finite_set.project(&mut x);
+    finite_set.project(&mut x).unwrap();
     unit_test_utils::assert_nearly_equal_array(
         &[0.0, 0.0],
         &x,
@@ -139,7 +139,7 @@ fn t_finite_set() {
         "projection is wrong (should be [0,0])",
     );
     x = [0.48, 0.501];
-    finite_set.project(&mut x);
+    finite_set.project(&mut x).unwrap();
     unit_test_utils::assert_nearly_equal_array(
         &[0.0, 1.0],
         &x,
@@ -148,7 +148,7 @@ fn t_finite_set() {
         "projection is wrong (should be [0,1])",
     );
     x = [0.7, 0.2];
-    finite_set.project(&mut x);
+    finite_set.project(&mut x).unwrap();
     unit_test_utils::assert_nearly_equal_array(
         &[1.0, 0.0],
         &x,
@@ -163,7 +163,7 @@ fn t_finite_set_f32() {
     let data: &[&[f32]] = &[&[0.0_f32, 0.0], &[1.0, 1.0], &[0.0, 1.0], &[1.0, 0.0]];
     let finite_set = FiniteSet::new(data);
     let mut x = [0.7_f32, 0.2];
-    finite_set.project(&mut x);
+    finite_set.project(&mut x).unwrap();
     assert_eq!([1.0_f32, 0.0], x);
 }
 
@@ -173,7 +173,7 @@ fn t_finite_set_project_wrong_dimension() {
     let data: &[&[f64]] = &[&[0.0, 0.0], &[1.0, 1.0]];
     let finite_set = FiniteSet::new(data);
     let mut x = [0.5, 0.5, 0.5];
-    finite_set.project(&mut x);
+    finite_set.project(&mut x).unwrap();
 }
 
 #[test]
@@ -183,7 +183,7 @@ fn t_rectangle_bounded() {
     let rectangle = Rectangle::new(Some(&xmin[..]), Some(&xmax[..]));
     let mut x = [1.0, 2.0, 3.0, 4.0, 5.0];
 
-    rectangle.project(&mut x);
+    rectangle.project(&mut x).unwrap();
 
     unit_test_utils::assert_nearly_equal_array(
         &[2.0, 2.0, 3.0, 4.0, 4.5],
@@ -201,7 +201,7 @@ fn t_rectangle_bounded_f32() {
     let rectangle = Rectangle::new(Some(&xmin[..]), Some(&xmax[..]));
     let mut x = [1.0_f32, 3.0, 5.0];
 
-    rectangle.project(&mut x);
+    rectangle.project(&mut x).unwrap();
 
     assert_eq!([2.0_f32, 3.0, 4.5], x);
 }
@@ -213,7 +213,7 @@ fn t_rectangle_infinite_bounds() {
     let rectangle = Rectangle::new(Some(&xmin[..]), Some(&xmax[..]));
     let mut x = [-2.0, 3.0, 1.0];
 
-    rectangle.project(&mut x);
+    rectangle.project(&mut x).unwrap();
 
     unit_test_utils::assert_nearly_equal_array(
         &[-1.0, 3.0, 1.0],
@@ -247,7 +247,7 @@ fn t_rectangle_bounded_negative_entries() {
     let rectangle = Rectangle::new(Some(&xmin[..]), Some(&xmax[..]));
     let mut x = [-6.0, -3.0, 0.0, 3.0, -5.0, 1.0, 2.0, 3.0, -1.0, 0.0, 0.0];
 
-    rectangle.project(&mut x);
+    rectangle.project(&mut x).unwrap();
 
     unit_test_utils::assert_nearly_equal_array(
         &[-5.0, -3.0, -1.0, 2.0, -1.0, 0.0, 2.0, 3.0, 3.0, 4.0, 5.0],
@@ -264,7 +264,7 @@ fn t_rectangle_only_xmin() {
     let rectangle = Rectangle::new(Some(&xmin[..]), None);
     let mut x = [1.0, 2.0, 3.0, 4.0, 5.0];
 
-    rectangle.project(&mut x);
+    rectangle.project(&mut x).unwrap();
 
     unit_test_utils::assert_nearly_equal_array(
         &[2.0, 2.0, 3.0, 4.0, 5.0],
@@ -281,7 +281,7 @@ fn t_rectangle_only_xmax() {
     let rectangle = Rectangle::new(None, Some(&xmax[..]));
     let mut x = [-10.0, -20.0, 0.0, 5.0, 3.0];
 
-    rectangle.project(&mut x);
+    rectangle.project(&mut x).unwrap();
 
     unit_test_utils::assert_nearly_equal_array(
         &[-10.0, -20.0, -3.0, -3.0, -3.0],
@@ -298,7 +298,7 @@ fn t_ball2_at_origin() {
     let mut x = [1.0, 1.0];
     let ball = Ball2::new(None, radius);
 
-    ball.project(&mut x);
+    ball.project(&mut x).unwrap();
 
     unit_test_utils::assert_nearly_equal_array(
         &[
@@ -318,7 +318,7 @@ fn t_ball2_at_origin_f32() {
     let mut x = [1.0_f32, 1.0];
     let ball = Ball2::new(None, radius);
 
-    ball.project(&mut x);
+    ball.project(&mut x).unwrap();
 
     let expected = std::f32::consts::FRAC_1_SQRT_2;
     assert!((x[0] - expected).abs() < 1e-6);
@@ -330,7 +330,7 @@ fn t_ball2_at_origin_different_radius_outside() {
     let radius = 0.8;
     let mut x = [1.0, 1.0];
     let ball = Ball2::new(None, radius);
-    ball.project(&mut x);
+    ball.project(&mut x).unwrap();
     let norm_proj_x = crate::matrix_operations::norm2(&x);
     unit_test_utils::assert_nearly_equal(radius, norm_proj_x, 1e-10, 1e-12, "wrong norm");
 }
@@ -340,7 +340,7 @@ fn t_ball2_at_origin_different_radius_inside() {
     let radius = 0.8;
     let mut x = [-0.2, 0.15];
     let ball = Ball2::new(None, radius);
-    ball.project(&mut x);
+    ball.project(&mut x).unwrap();
     unit_test_utils::assert_nearly_equal_array(&x, &[-0.2, 0.15], 1e-10, 1e-12, "wrong");
 }
 
@@ -350,9 +350,8 @@ fn t_ball2_at_center_different_radius_outside() {
     let mut x = [1.0, 1.0];
     let center = [-0.8, -1.1];
     let ball = Ball2::new(Some(&center), radius);
-    ball.project(&mut x);
-    let norm_sq_x_minus_c: f64 = crate::matrix_operations::norm2_squared_diff(&x, &center);
-    let norm_x_minus_c = norm_sq_x_minus_c.sqrt();
+    ball.project(&mut x).unwrap();
+    let norm_x_minus_c = crate::matrix_operations::norm2_squared_diff(&x, &center).sqrt();
     unit_test_utils::assert_nearly_equal(radius, norm_x_minus_c, 1e-10, 1e-12, "wrong norm");
 }
 
@@ -362,7 +361,7 @@ fn t_ball2_at_center_different_radius_inside() {
     let mut x = [-0.9, -0.85];
     let center = [-0.8, -1.1];
     let ball = Ball2::new(Some(&center), radius);
-    ball.project(&mut x);
+    ball.project(&mut x).unwrap();
     unit_test_utils::assert_nearly_equal_array(&[-0.9, -0.85], &x, 1e-10, 1e-12, "wrong result");
 }
 
@@ -373,7 +372,7 @@ fn t_ball2_elsewhere() {
     let mut x = [2.0, 2.0];
     let ball = Ball2::new(Some(&center[..]), radius);
 
-    ball.project(&mut x);
+    ball.project(&mut x).unwrap();
 
     let expected_proj_element = std::f64::consts::FRAC_1_SQRT_2 + 1.;
     unit_test_utils::assert_nearly_equal_array(
@@ -391,7 +390,7 @@ fn t_ball2_boundary_no_change() {
     let mut x = [0.0, 2.0];
     let x_expected = x;
     let ball = Ball2::new(None, radius);
-    ball.project(&mut x);
+    ball.project(&mut x).unwrap();
     unit_test_utils::assert_nearly_equal_array(&x_expected, &x, 1e-12, 1e-12, "wrong result");
 }
 
@@ -400,7 +399,7 @@ fn t_no_constraints() {
     let mut x = [1.0, 2.0, 3.0];
     let whole_space = NoConstraints::new();
 
-    whole_space.project(&mut x);
+    whole_space.project(&mut x).unwrap();
 
     unit_test_utils::assert_nearly_equal_array(&[1., 2., 3.], &x, 1e-10, 1e-15, "x is wrong");
 }
@@ -410,7 +409,7 @@ fn t_no_constraints_f32() {
     let mut x = [1.0_f32, 2.0, 3.0];
     let whole_space = NoConstraints::new();
 
-    whole_space.project(&mut x);
+    whole_space.project(&mut x).unwrap();
 
     assert_eq!([1.0_f32, 2.0, 3.0], x);
 }
@@ -434,7 +433,7 @@ fn t_cartesian_product_constraints_wrong_vector_dim() {
         .add_constraint(3, ball1)
         .add_constraint(10, ball2);
     let mut x = [0.0; 30];
-    cart_prod.project(&mut x);
+    cart_prod.project(&mut x).unwrap();
 }
 
 #[test]
@@ -449,7 +448,7 @@ fn t_cartesian_product_constraints() {
         .add_constraint(idx1, ball1)
         .add_constraint(idx2, ball2);
     let mut x = [3.0, 4.0, 5.0, 2.0, 1.0];
-    cart_prod.project(&mut x);
+    cart_prod.project(&mut x).unwrap();
     let r1 = crate::matrix_operations::norm2(&x[0..idx1]);
     let r2 = crate::matrix_operations::norm2(&x[idx1..idx2]);
     unit_test_utils::assert_nearly_equal(r1, radius1, 1e-8, 1e-12, "r1 is wrong");
@@ -480,7 +479,7 @@ fn t_cartesian_product_ball_and_rectangle() {
 
     /* Projection */
     let mut x = [-10.0, 0.5, 10.0, 0.01, -0.01, 0.1, 10.0, -1.0, 1.0];
-    cart_prod.project(&mut x);
+    cart_prod.project(&mut x).unwrap();
 
     unit_test_utils::assert_nearly_equal_array(
         &x[0..3],
@@ -516,7 +515,7 @@ fn t_cartesian_product_ball_and_rectangle_f32() {
         .add_constraint(5, ball);
 
     let mut x = [-4.0_f32, 0.25_f32, 2.0_f32, -1.0_f32, 2.0_f32];
-    cart_prod.project(&mut x);
+    cart_prod.project(&mut x).unwrap();
 
     assert_eq!([-1.0_f32, 0.25_f32], x[..2]);
     let ball_norm = crate::matrix_operations::norm2(&x[2..5]);
@@ -528,7 +527,7 @@ fn t_second_order_cone_case_i() {
     let soc = SecondOrderCone::new(1.0);
     let mut x = vec![1.0, 1.0, 1.42];
     let x_copy = x.clone();
-    soc.project(&mut x);
+    soc.project(&mut x).unwrap();
     unit_test_utils::assert_nearly_equal_array(&x, &x_copy, 1e-10, 1e-12, "x has been modified");
 }
 
@@ -537,7 +536,7 @@ fn t_second_order_cone_case_ii() {
     let alpha = 0.5;
     let soc = SecondOrderCone::new(alpha);
     let mut x = vec![1.0, 1.0, -0.71];
-    soc.project(&mut x);
+    soc.project(&mut x).unwrap();
     let expected = vec![0.0; 3];
     unit_test_utils::assert_nearly_equal_array(
         &x,
@@ -553,7 +552,7 @@ fn t_second_order_cone_case_iii() {
     let alpha = 1.5;
     let soc = SecondOrderCone::new(alpha);
     let mut x = vec![1.0, 1.0, 0.1];
-    soc.project(&mut x);
+    soc.project(&mut x).unwrap();
     // make sure the new `x` is in the cone
     let norm_z: f64 = crate::matrix_operations::norm2(&x[..=1]);
     assert!(norm_z <= alpha * x[2]);
@@ -566,7 +565,7 @@ fn t_second_order_cone_case_iii_f32() {
     let alpha = 1.5_f32;
     let soc = SecondOrderCone::new(alpha);
     let mut x = vec![1.0_f32, 1.0_f32, 0.1_f32];
-    soc.project(&mut x);
+    soc.project(&mut x).unwrap();
     let norm_z = crate::matrix_operations::norm2(&x[..=1]);
     assert!(norm_z <= alpha * x[2] + 1e-5_f32);
     assert!((norm_z - alpha * x[2]).abs() <= 1e-4_f32);
@@ -592,7 +591,7 @@ fn t_second_order_cone_short_vector() {
     let alpha = 1.0;
     let soc = SecondOrderCone::new(alpha);
     let mut _x = vec![1.0];
-    soc.project(&mut _x);
+    soc.project(&mut _x).unwrap();
 }
 
 #[test]
@@ -612,7 +611,7 @@ fn t_cartesian_product_dimension() {
     // let's do a projection to make sure this works
     // Note: we've used the same set (finite_set), twice
     let mut x = [-0.5, 1.1, 0.45, 0.55, 10.0, 10.0, -500.0, 1.0, 1.0, 1.0];
-    cartesian.project(&mut x);
+    cartesian.project(&mut x).unwrap();
     println!("X = {:#?}", x);
     let sqrt_3_over_3 = 3.0_f64.sqrt() / 3.;
     unit_test_utils::assert_nearly_equal_array(
@@ -648,7 +647,7 @@ fn t_cartesian_ball_no_constraint() {
         .add_constraint(9, no_constraints);
     assert_eq!(9, cartesian.dimension());
     let mut x = [100., -200., 0.5, 1.5, 3.5, 1000., 5., -500., 2_000_000.];
-    cartesian.project(&mut x);
+    cartesian.project(&mut x).unwrap();
     let x_proj_ball = [0.869811089019176, 0.390566732942472, 0.911322376865767];
     unit_test_utils::assert_nearly_equal_array(
         &x[0..=1],
@@ -672,7 +671,7 @@ fn t_ball_inf_origin() {
     let ball_inf = BallInf::new(None, 1.0);
     let mut x = [0.0, -0.5, 0.5, 1.5, 3.5, 0.8, 1.1, -5.0, -10.0];
     let x_correct = [0.0, -0.5, 0.5, 1.0, 1.0, 0.8, 1.0, -1.0, -1.0];
-    ball_inf.project(&mut x);
+    ball_inf.project(&mut x).unwrap();
     unit_test_utils::assert_nearly_equal_array(
         &x_correct,
         &x,
@@ -688,31 +687,31 @@ fn t_ball_inf_center() {
     let xc = [5.0, -6.0];
     let ball_inf = BallInf::new(Some(&xc), 1.5);
     let mut x = [11.0, -0.5];
-    ball_inf.project(&mut x);
+    ball_inf.project(&mut x).unwrap();
     unit_test_utils::assert_nearly_equal_array(&[6.5, -4.5], &x, 1e-10, 1e-12, "upper right");
 
     let mut x = [3.0, -7.0];
-    ball_inf.project(&mut x);
+    ball_inf.project(&mut x).unwrap();
     unit_test_utils::assert_nearly_equal_array(&[3.5, -7.0], &x, 1e-10, 1e-12, "left");
 
     let mut x = [800.0, -5.0];
-    ball_inf.project(&mut x);
+    ball_inf.project(&mut x).unwrap();
     unit_test_utils::assert_nearly_equal_array(&[6.5, -5.0], &x, 1e-10, 1e-12, "right");
 
     let mut x = [9.0, -10.0];
-    ball_inf.project(&mut x);
+    ball_inf.project(&mut x).unwrap();
     unit_test_utils::assert_nearly_equal_array(&[6.5, -7.5], &x, 1e-10, 1e-12, "down right");
 
     let mut x = [3.0, 0.0];
-    ball_inf.project(&mut x);
+    ball_inf.project(&mut x).unwrap();
     unit_test_utils::assert_nearly_equal_array(&[3.5, -4.5], &x, 1e-10, 1e-12, "top left");
 
     let mut x = [6.0, -5.0];
-    ball_inf.project(&mut x);
+    ball_inf.project(&mut x).unwrap();
     unit_test_utils::assert_nearly_equal_array(&[6.0, -5.0], &x, 1e-10, 1e-12, "inside");
 
     let mut x = [5.0, -6.0];
-    ball_inf.project(&mut x);
+    ball_inf.project(&mut x).unwrap();
     unit_test_utils::assert_nearly_equal_array(&[5.0, -6.0], &x, 1e-10, 1e-12, "centre");
 }
 
@@ -721,7 +720,7 @@ fn t_ball_inf_boundary_no_change() {
     let ball_inf = BallInf::new(None, 1.0);
     let mut x = [-1.0, 0.2, 1.0];
     let x_expected = x;
-    ball_inf.project(&mut x);
+    ball_inf.project(&mut x).unwrap();
     unit_test_utils::assert_nearly_equal_array(&x_expected, &x, 1e-12, 1e-12, "wrong result");
 }
 
@@ -795,7 +794,7 @@ fn t_simplex_projection() {
     let mut x = [1.0, 2.0, 3.0];
     let alpha = 3.0;
     let my_simplex = Simplex::new(alpha);
-    my_simplex.project(&mut x);
+    my_simplex.project(&mut x).unwrap();
     unit_test_utils::assert_nearly_equal(
         crate::matrix_operations::sum(&x),
         alpha,
@@ -810,7 +809,7 @@ fn t_simplex_projection_f32() {
     let mut x = [1.0_f32, 2.0, 3.0];
     let alpha = 3.0_f32;
     let simplex = Simplex::new(alpha);
-    simplex.project(&mut x);
+    simplex.project(&mut x).unwrap();
 
     let sum = x[0] + x[1] + x[2];
     assert!((sum - alpha).abs() < 1e-5);
@@ -824,7 +823,7 @@ fn t_halfspace_boundary_no_change() {
     let halfspace = Halfspace::new(&normal_vector, offset);
     let mut x = [1.0, 2.0];
     let x_expected = x;
-    halfspace.project(&mut x);
+    halfspace.project(&mut x).unwrap();
     unit_test_utils::assert_nearly_equal_array(&x_expected, &x, 1e-12, 1e-12, "wrong result");
 }
 
@@ -840,8 +839,7 @@ fn t_simplex_projection_random_spam() {
         let alpha_scale = 20.;
         let alpha = alpha_scale * rand::random::<f64>();
         let simplex = Simplex::new(alpha);
-        simplex.project(&mut x);
-        println!("x = {:?}", x);
+        simplex.project(&mut x).unwrap();
         assert!(x.iter().all(|&xi| xi >= -1e-12));
         unit_test_utils::assert_nearly_equal(
             crate::matrix_operations::sum(&x),
@@ -865,7 +863,7 @@ fn t_simplex_projection_random_optimality() {
             let alpha = alpha_scale * rand::random::<f64>();
             let simplex = Simplex::new(alpha);
             let y = z.clone();
-            simplex.project(&mut z);
+            simplex.project(&mut z).unwrap();
             for j in 0..n {
                 let w = alpha * (y[j] - z[j]) - crate::matrix_operations::inner_product(&z, &y)
                     + crate::matrix_operations::norm2_squared(&z);
@@ -897,7 +895,7 @@ fn t_simplex_alpha_negative() {
 fn t_simplex_empty_vector() {
     let simplex = Simplex::new(1.0);
     let mut x = [];
-    simplex.project(&mut x);
+    simplex.project(&mut x).unwrap();
 }
 
 #[test]
@@ -914,7 +912,7 @@ fn t_ball1_random_optimality_conditions() {
             x.copy_from_slice(&x_star);
             let radius = 5. * rand::random::<f64>();
             let ball1 = Ball1::new(None, radius);
-            ball1.project(&mut x_star);
+            ball1.project(&mut x_star).unwrap();
             // make sure |x|_1 <= radius
             assert!(
                 crate::matrix_operations::norm1(&x_star) <= radius * (1. + 1e-9),
@@ -952,7 +950,7 @@ fn t_ball1_random_optimality_conditions() {
 fn t_ball1_projection_f32() {
     let ball1 = Ball1::new(None, 1.0_f32);
     let mut x = [2.0_f32, -1.0_f32, 0.0_f32];
-    ball1.project(&mut x);
+    ball1.project(&mut x).unwrap();
     assert!((x[0] - 1.0_f32).abs() < 1e-6_f32);
     assert!(x[1].abs() < 1e-6_f32);
     assert!(x[2].abs() < 1e-6_f32);
@@ -974,13 +972,13 @@ fn t_ball1_random_optimality_conditions_centered() {
                 .for_each(|xi| *xi = scale_xc * (2. * rand::random::<f64>() - 1.));
             let radius = 5. * rand::random::<f64>();
             let ball1 = Ball1::new(Some(&xc), radius);
-            ball1.project(&mut x);
+            ball1.project(&mut x).unwrap();
             // x = x - xc
             x.iter_mut()
                 .zip(xc.iter())
                 .for_each(|(xi, &xci)| *xi -= xci);
             assert!(
-                crate::matrix_operations::norm1(&x) <= radius * (1. + 1e-9),
+                crate::matrix_operations::norm1(&x) <= radius * (1. + 1e-9) + 1e-12,
                 "norm(x - xc, 1) > radius"
             );
         }
@@ -994,7 +992,7 @@ fn t_ball1_wrong_dimensions() {
     let mut x = vec![3.0, 4.0, 5.0];
     let radius = 1.0;
     let ball1 = Ball1::new(Some(&xc), radius);
-    ball1.project(&mut x);
+    ball1.project(&mut x).unwrap();
 }
 
 #[test]
@@ -1003,8 +1001,8 @@ fn t_sphere2_no_center() {
     let mut x_out = [1.0, 1.0];
     let mut x_in = [-0.3, -0.2];
     let unit_sphere = Sphere2::new(None, radius);
-    unit_sphere.project(&mut x_out);
-    unit_sphere.project(&mut x_in);
+    unit_sphere.project(&mut x_out).unwrap();
+    unit_sphere.project(&mut x_in).unwrap();
     let norm_out = crate::matrix_operations::norm2(&x_out);
     let norm_in = crate::matrix_operations::norm2(&x_in);
     unit_test_utils::assert_nearly_equal(radius, norm_out, 1e-10, 1e-12, "norm_out is not 1.0");
@@ -1016,7 +1014,7 @@ fn t_sphere2_no_center_projection_of_zero() {
     let radius = 0.9;
     let mut x = [0.0, 0.0];
     let unit_sphere = Sphere2::new(None, radius);
-    unit_sphere.project(&mut x);
+    unit_sphere.project(&mut x).unwrap();
     let norm_result = crate::matrix_operations::norm2(&x);
     unit_test_utils::assert_nearly_equal(radius, norm_result, 1e-10, 1e-12, "norm_out is not 1.0");
 }
@@ -1028,7 +1026,7 @@ fn t_sphere2_center() {
     let mut x = [1.0, 1.0];
     let unit_sphere = Sphere2::new(Some(&center), radius);
 
-    unit_sphere.project(&mut x);
+    unit_sphere.project(&mut x).unwrap();
     let mut x_minus_c = [0.0; 2];
     x.iter()
         .zip(center.iter())
@@ -1048,7 +1046,7 @@ fn t_sphere2_center_projection_of_center() {
     let mut x = [-3.0, 5.0];
     let unit_sphere = Sphere2::new(Some(&center), radius);
 
-    unit_sphere.project(&mut x);
+    unit_sphere.project(&mut x).unwrap();
     let mut x_minus_c = [0.0; 2];
     x.iter()
         .zip(center.iter())
@@ -1067,7 +1065,7 @@ fn t_sphere2_empty_vector() {
     let radius = 1.0;
     let unit_sphere = Sphere2::new(None, radius);
     let mut x = [];
-    unit_sphere.project(&mut x);
+    unit_sphere.project(&mut x).unwrap();
 }
 
 #[test]
@@ -1077,7 +1075,7 @@ fn t_sphere2_center_wrong_dimension() {
     let center = [1.0, 2.0, 3.0];
     let unit_sphere = Sphere2::new(Some(&center), radius);
     let mut x = [1.0, 2.0];
-    unit_sphere.project(&mut x);
+    unit_sphere.project(&mut x).unwrap();
 }
 
 #[test]
@@ -1090,8 +1088,8 @@ fn t_ball1_alpha_negative() {
 fn t_epigraph_squared_norm_inside() {
     let epi = EpigraphSquaredNorm::new();
     let mut x = [1., 2., 10.];
+    epi.project(&mut x).unwrap();
     let x_correct = x;
-    epi.project(&mut x);
     unit_test_utils::assert_nearly_equal_array(
         &x_correct,
         &x,
@@ -1106,7 +1104,7 @@ fn t_epigraph_squared_norm_boundary_no_change() {
     let epi = EpigraphSquaredNorm::new();
     let mut x = [1.0, 2.0, 5.0];
     let x_expected = x;
-    epi.project(&mut x);
+    epi.project(&mut x).unwrap();
     unit_test_utils::assert_nearly_equal_array(&x_expected, &x, 1e-12, 1e-12, "wrong result");
 }
 
@@ -1116,7 +1114,7 @@ fn t_epigraph_squared_norm() {
     for i in 0..100 {
         let t = 0.01 * i as f64;
         let mut x = [1., 2., 3., t];
-        epi.project(&mut x);
+        epi.project(&mut x).unwrap();
         let err = (matrix_operations::norm2_squared(&x[..3]) - x[3]).abs();
         assert!(err < 1e-10, "wrong projection on epigraph of squared norm");
     }
@@ -1132,7 +1130,7 @@ fn t_epigraph_squared_norm_correctness() {
         1.680426686710711,
         4.392630432414829,
     ];
-    epi.project(&mut x);
+    epi.project(&mut x).unwrap();
     unit_test_utils::assert_nearly_equal_array(
         &x_correct,
         &x,
@@ -1146,7 +1144,7 @@ fn t_epigraph_squared_norm_correctness() {
 fn t_epigraph_squared_norm_f32() {
     let epi = EpigraphSquaredNorm::new();
     let mut x = [1.0_f32, 0.0, 0.0];
-    epi.project(&mut x);
+    epi.project(&mut x).unwrap();
     let err = (matrix_operations::norm2_squared(&x[..2]) - x[2]).abs();
     assert!(err < 1e-4);
 }
@@ -1159,7 +1157,7 @@ fn t_affine_space() {
     let b = vec![1., 2., -0.5];
     let affine_set = AffineSpace::new(a, b);
     let mut x = [1., -2., -0.3, 0.5];
-    affine_set.project(&mut x);
+    affine_set.project(&mut x).unwrap();
     let x_correct = [
         1.888564346697095,
         5.629857182200888,
@@ -1183,7 +1181,7 @@ fn t_affine_space_f32() {
     let b = vec![1.0_f32, 2.0, -0.5];
     let affine_set = AffineSpace::new(a.clone(), b.clone());
     let mut x = [1.0_f32, -2.0, -0.3, 0.5];
-    affine_set.project(&mut x);
+    affine_set.project(&mut x).unwrap();
 
     let x_correct = [
         1.888_564_3_f32,
@@ -1213,7 +1211,7 @@ fn t_affine_space_projection_feasibility() {
     let b = vec![1., 2., -0.5];
     let affine_set = AffineSpace::new(a.clone(), b.clone());
     let mut x = [1., -2., -0.3, 0.5];
-    affine_set.project(&mut x);
+    affine_set.project(&mut x).unwrap();
     let residual = [
         a[0] * x[0] + a[1] * x[1] + a[2] * x[2] + a[3] * x[3] - b[0],
         a[4] * x[0] + a[5] * x[1] + a[6] * x[2] + a[7] * x[3] - b[1],
@@ -1233,7 +1231,7 @@ fn t_affine_space_larger() {
     let b = vec![1., -2., 3., 4.];
     let affine_set = AffineSpace::new(a, b);
     let mut x = [10., 11., -9., 4., 5.];
-    affine_set.project(&mut x);
+    affine_set.project(&mut x).unwrap();
     let x_correct = [
         9.238095238095237,
         -0.714285714285714,
@@ -1256,7 +1254,7 @@ fn t_affine_space_single_row() {
     let b = vec![1.];
     let affine_set = AffineSpace::new(a, b);
     let mut x = [5., 6., 10., 25.];
-    affine_set.project(&mut x);
+    affine_set.project(&mut x).unwrap();
     let s = x.iter().sum();
     unit_test_utils::assert_nearly_equal(1., s, 1e-12, 1e-14, "wrong sum");
 }
@@ -1533,7 +1531,7 @@ where
             let p = cast::<T>(p_f64);
             let radius = cast::<T>(radius_f64);
             let ball = BallP::new(center.as_deref(), radius, p, solver_tol, 300);
-            ball.project(&mut x);
+            ball.project(&mut x).unwrap();
 
             let shifted_projection: Vec<T> = if let Some(center) = center.as_ref() {
                 x.iter()
@@ -1550,7 +1548,7 @@ where
             );
 
             let mut reproj = x.clone();
-            ball.project(&mut reproj);
+            ball.project(&mut reproj).unwrap();
             let max_reproj_diff = reproj
                 .iter()
                 .zip(x.iter())
@@ -1620,7 +1618,7 @@ where
             x.push(cast::<T>(rng.random_range(-2.0..4.0)));
             let x_before = as_f64_vec(&x);
 
-            epi.project(&mut x);
+            epi.project(&mut x).unwrap();
 
             let z = &x[..dim];
             let t = x[dim];
@@ -1631,7 +1629,7 @@ where
             );
 
             let mut reproj = x.clone();
-            epi.project(&mut reproj);
+            epi.project(&mut reproj).unwrap();
             let max_reproj_diff = reproj
                 .iter()
                 .zip(x.iter())
@@ -1674,9 +1672,9 @@ where
 fn assert_projection_idempotent<C: Constraint>(constraint: &C, x0: &[f64], message: &'static str) {
     let mut once = x0.to_vec();
     let mut twice = x0.to_vec();
-    constraint.project(&mut once);
-    constraint.project(&mut twice);
-    constraint.project(&mut twice);
+    constraint.project(&mut once).unwrap();
+    constraint.project(&mut twice).unwrap();
+    constraint.project(&mut twice).unwrap();
     unit_test_utils::assert_nearly_equal_array(&once, &twice, 1e-10, 1e-12, message);
 }
 
@@ -1689,7 +1687,7 @@ fn t_ballp_at_origin_projection() {
     let tol = 1e-16;
     let max_iters: usize = 200;
     let ball = BallP::new(None, radius, p, tol, max_iters);
-    ball.project(&mut x);
+    ball.project(&mut x).unwrap();
     assert!(is_norm_p_projection(&x0, &x, p, radius, 10_000));
 }
 
@@ -1699,7 +1697,7 @@ fn t_ballp_at_origin_projection_preserves_signs() {
     let mut x = [1.0, -3.0, 2.5, -0.7];
     let x0 = x;
     let ball = BallP::new(None, radius, 3.0, 1e-14, 200);
-    ball.project(&mut x);
+    ball.project(&mut x).unwrap();
     for (proj, original) in x.iter().zip(x0.iter()) {
         assert!(proj.abs() <= original.abs() + 1e-12);
         if *original != 0.0 {
@@ -1715,7 +1713,7 @@ fn t_ballp_zero_coordinates_branch() {
     let mut x = [0.0, -2.0, 0.0, 1.5];
     let x0 = x;
     let ball = BallP::new(None, radius, p, 1e-14, 300);
-    ball.project(&mut x);
+    ball.project(&mut x).unwrap();
     assert_eq!(x[0], 0.0);
     assert_eq!(x[2], 0.0);
     assert!(is_norm_p_projection(&x0, &x, p, radius, 10_000));
@@ -1734,7 +1732,7 @@ fn t_ballp_outside_projection_lands_on_boundary_for_multiple_p() {
     for (p, x_init) in test_cases {
         let mut x = x_init;
         let ball = BallP::new(None, radius, p, 1e-14, 400);
-        ball.project(&mut x);
+        ball.project(&mut x).unwrap();
         let norm_p = x
             .iter()
             .map(|xi| xi.abs().powf(p))
@@ -1757,7 +1755,7 @@ fn t_ballp_boundary_no_change() {
     let mut x = [1.0, 0.0];
     let x_expected = x;
     let ball = BallP::new(None, radius, p, 1e-14, 200);
-    ball.project(&mut x);
+    ball.project(&mut x).unwrap();
     unit_test_utils::assert_nearly_equal_array(&x_expected, &x, 1e-12, 1e-12, "wrong result");
 }
 
@@ -1775,7 +1773,7 @@ fn t_ballp_translated_projection_multiple_p_values() {
     for (p, x_init) in cases {
         let mut x = x_init;
         let ball = BallP::new(Some(&center), radius, p, 1e-14, 400);
-        ball.project(&mut x);
+        ball.project(&mut x).unwrap();
         let norm_p = x
             .iter()
             .zip(center.iter())
@@ -1868,7 +1866,7 @@ fn t_ballp_at_origin_x_already_inside() {
     let tol = 1e-16;
     let max_iters: usize = 1200;
     let ball = BallP::new(None, radius, p, tol, max_iters);
-    ball.project(&mut x);
+    ball.project(&mut x).unwrap();
     unit_test_utils::assert_nearly_equal_array(
         &x0,
         &x,
@@ -1887,7 +1885,7 @@ fn t_ballp_at_xc_projection() {
     let tol = 1e-16;
     let max_iters: usize = 200;
     let ball = BallP::new(Some(&x_center), radius, p, tol, max_iters);
-    ball.project(&mut x);
+    ball.project(&mut x).unwrap();
 
     let nrm: f64 = (x.iter().zip(x_center.iter()).fold(0.0_f64, |s, (x, y)| {
         let diff: f64 = *x - *y;
@@ -1915,7 +1913,7 @@ fn t_ballp_at_xc_projection_f32() {
     let tol = 1e-6_f32;
     let max_iters: usize = 200;
     let ball = BallP::new(Some(&x_center), radius, p, tol, max_iters);
-    ball.project(&mut x);
+    ball.project(&mut x).unwrap();
 
     let nrm = x
         .iter()
@@ -1961,7 +1959,7 @@ fn t_rectangle_only_xmin_wrong_dimension() {
     let xmin = [1.0, 2.0, 3.0];
     let rectangle = Rectangle::new(Some(&xmin), None);
     let mut x = [0.0, 1.0];
-    rectangle.project(&mut x);
+    rectangle.project(&mut x).unwrap();
 }
 
 #[test]
@@ -1970,7 +1968,7 @@ fn t_rectangle_only_xmax_wrong_dimension() {
     let xmax = [1.0, 2.0, 3.0];
     let rectangle = Rectangle::new(None, Some(&xmax));
     let mut x = [0.0, 1.0];
-    rectangle.project(&mut x);
+    rectangle.project(&mut x).unwrap();
 }
 
 #[test]
@@ -1979,7 +1977,7 @@ fn t_halfspace_wrong_dimension() {
     let normal_vector = [1.0, 2.0, 3.0];
     let halfspace = Halfspace::new(&normal_vector, 1.0);
     let mut x = [1.0, 2.0];
-    halfspace.project(&mut x);
+    halfspace.project(&mut x).unwrap();
 }
 
 #[test]
@@ -1988,7 +1986,7 @@ fn t_ball2_wrong_dimensions() {
     let center = [1.0, 2.0];
     let ball = Ball2::new(Some(&center), 1.0);
     let mut x = [1.0, 2.0, 3.0];
-    ball.project(&mut x);
+    ball.project(&mut x).unwrap();
 }
 
 #[test]
@@ -2003,7 +2001,7 @@ fn t_ball_inf_wrong_dimensions() {
     let center = [1.0, 2.0];
     let ball_inf = BallInf::new(Some(&center), 1.0);
     let mut x = [1.0, 2.0, 3.0];
-    ball_inf.project(&mut x);
+    ball_inf.project(&mut x).unwrap();
 }
 
 #[test]
@@ -2017,7 +2015,7 @@ fn t_ball_inf_nonpositive_radius() {
 fn t_epigraph_squared_norm_short_vector() {
     let epi = EpigraphSquaredNorm::new();
     let mut x = [1.0];
-    epi.project(&mut x);
+    epi.project(&mut x).unwrap();
 }
 
 #[test]
@@ -2033,7 +2031,7 @@ fn t_affine_space_project_wrong_dimension() {
     let b = vec![0.0, 0.0];
     let affine_set = AffineSpace::new(a, b);
     let mut x = [1.0];
-    affine_set.project(&mut x);
+    affine_set.project(&mut x).unwrap();
 }
 
 #[test]
@@ -2132,5 +2130,5 @@ fn t_ballp_wrong_dimensions() {
     let center = [1.0, 2.0];
     let ballp = BallP::new(Some(&center), 1.0, 3.0, 1e-12, 100);
     let mut x = [1.0, 2.0, 3.0];
-    ballp.project(&mut x);
+    ballp.project(&mut x).unwrap();
 }
diff --git a/src/constraints/zero.rs b/src/constraints/zero.rs
index e2bb1541..c9a7e503 100644
--- a/src/constraints/zero.rs
+++ b/src/constraints/zero.rs
@@ -1,4 +1,5 @@
 use super::Constraint;
+use crate::FunctionCallResult;
 use num::Float;
 
 #[derive(Clone, Copy, Default)]
@@ -15,7 +16,7 @@ impl Zero {
     ///
     /// let zero = Zero::new();
     /// let mut x = [1.0, -2.0, 3.0];
-    /// zero.project(&mut x);
+    /// zero.project(&mut x).unwrap();
     /// ```
     #[must_use]
     pub fn new() -> Self {
@@ -26,8 +27,9 @@ impl Zero {
 impl<T: Float> Constraint<T> for Zero {
     /// Computes the projection on $\\{0\\}$, that is, $\Pi_{\\{0\\}}(x) = 0$
     /// for all $x$
-    fn project(&self, x: &mut [T]) {
+    fn project(&self, x: &mut [T]) -> FunctionCallResult {
         x.iter_mut().for_each(|xi| *xi = T::zero());
+        Ok(())
     }
 
     fn is_convex(&self) -> bool {
diff --git a/src/core/fbs/fbs_engine.rs b/src/core/fbs/fbs_engine.rs
index c13ae870..a7a38f32 100644
--- a/src/core/fbs/fbs_engine.rs
+++ b/src/core/fbs/fbs_engine.rs
@@ -45,22 +45,24 @@ where
         FBSEngine { problem, cache }
     }
 
-    fn gradient_step(&mut self, u_current: &mut [T]) {
-        assert_eq!(
-            Ok(()),
-            (self.problem.gradf)(u_current, &mut self.cache.work_gradient_u),
-            "The computation of the gradient of the cost failed miserably"
-        );
+    fn gradient_step(&mut self, u_current: &mut [T]) -> FunctionCallResult {
+        (self.problem.gradf)(u_current, &mut self.cache.work_gradient_u)?;
+        if !crate::matrix_operations::is_finite(&self.cache.work_gradient_u) {
+            return Err(SolverError::NotFiniteComputation(
+                "gradient evaluation returned a non-finite value during an FBS step",
+            ));
+        }
 
         // take a gradient step: u_currect -= gamma * gradient
         u_current
             .iter_mut()
             .zip(self.cache.work_gradient_u.iter())
             .for_each(|(u, w)| *u = *u - self.cache.gamma * *w);
+        Ok(())
     }
 
-    fn projection_step(&mut self, u_current: &mut [T]) {
-        self.problem.constraints.project(u_current);
+    fn projection_step(&mut self, u_current: &mut [T]) -> FunctionCallResult {
+        self.problem.constraints.project(u_current)
     }
 }
 
@@ -89,8 +91,13 @@ where
     /// or the cost function panics.
     fn step(&mut self, u_current: &mut [T]) -> Result<bool, SolverError> {
         self.cache.work_u_previous.copy_from_slice(u_current); // cache the previous step
-        self.gradient_step(u_current); // compute the gradient
-        self.projection_step(u_current); // project
+        self.gradient_step(u_current)?; // compute the gradient
+        self.projection_step(u_current)?; // project
+        if !crate::matrix_operations::is_finite(u_current) {
+            return Err(SolverError::NotFiniteComputation(
+                "projected iterate contains a non-finite value during an FBS step",
+            ));
+        }
         self.cache.norm_fpr =
             matrix_operations::norm_inf_diff(u_current, &self.cache.work_u_previous);
 
diff --git a/src/core/fbs/fbs_optimizer.rs b/src/core/fbs/fbs_optimizer.rs
index f45d0a13..076a750a 100644
--- a/src/core/fbs/fbs_optimizer.rs
+++ b/src/core/fbs/fbs_optimizer.rs
@@ -122,7 +122,9 @@ where
         (self.fbs_engine.problem.cost)(u, &mut cost_value)?;
 
         if !matrix_operations::is_finite(u) || !cost_value.is_finite() {
-            return Err(SolverError::NotFiniteComputation);
+            return Err(SolverError::NotFiniteComputation(
+                "final FBS iterate or cost is non-finite",
+            ));
         }
 
         Ok(SolverStatus::new(
diff --git a/src/core/fbs/tests.rs b/src/core/fbs/tests.rs
index e3ed7a31..765caff7 100644
--- a/src/core/fbs/tests.rs
+++ b/src/core/fbs/tests.rs
@@ -55,7 +55,7 @@ fn t_solve_fbs_hard_failure_nan() {
     let mut u = [-12., -160., 55.];
     let mut optimizer = FBSOptimizer::new(problem, &mut fbs_cache).with_max_iter(10000);
     let status = optimizer.solve(&mut u);
-    assert_eq!(Err(SolverError::NotFiniteComputation), status);
+    assert!(matches!(status, Err(SolverError::NotFiniteComputation(_))));
 }
 
 #[test]
diff --git a/src/core/mod.rs b/src/core/mod.rs
index f0ea08a8..a329bdaa 100644
--- a/src/core/mod.rs
+++ b/src/core/mod.rs
@@ -35,9 +35,14 @@ pub trait Optimizer<T = f64>
 where
     T: Float,
 {
-    /// solves a given problem and updates the initial estimate `u` with the solution
+    /// Solves a given problem and updates the initial estimate `u` with the solution.
     ///
-    /// Returns the solver status
+    /// Returns the solver status on success.
+    ///
+    /// If the algorithm cannot proceed because a user callback fails, a
+    /// projection fails, non-finite values are encountered, or an internal
+    /// numerical/kernel inconsistency is detected, this method returns
+    /// `Err(SolverError)`.
     ///
     fn solve(&mut self, u: &mut [T]) -> Result<SolverStatus<T>, SolverError>;
 }
@@ -52,9 +57,16 @@ where
 /// the specified termination criterion is satisfied
 ///
 pub trait AlgorithmEngine<T = f64> {
-    /// Take a step of the algorithm and return `Ok(true)` only if the iterations should continue
+    /// Take a step of the algorithm and return `Ok(true)` only if the iterations should continue.
+    ///
+    /// Returns `Err(SolverError)` if a callback or projection fails, if a
+    /// non-finite value is produced, or if the engine detects an invalid
+    /// numerical state.
     fn step(&mut self, u: &mut [T]) -> Result<bool, SolverError>;
 
-    /// Initializes the algorithm
+    /// Initializes the algorithm.
+    ///
+    /// Returns `Err(SolverError)` if initialization requires evaluating a
+    /// callback/projection and that operation fails.
     fn init(&mut self, u: &mut [T]) -> FunctionCallResult;
 }
diff --git a/src/core/panoc/panoc_engine.rs b/src/core/panoc/panoc_engine.rs
index f021117c..c558ceb4 100644
--- a/src/core/panoc/panoc_engine.rs
+++ b/src/core/panoc/panoc_engine.rs
@@ -155,13 +155,14 @@ where
     }
 
     /// Computes a projection on `gradient_step`
-    fn half_step(&mut self) {
+    fn half_step(&mut self) -> FunctionCallResult {
         let cache = &mut self.cache;
         // u_half_step ← projection(gradient_step)
         cache.u_half_step.copy_from_slice(&cache.gradient_step);
-        self.problem.constraints.project(&mut cache.u_half_step);
+        self.problem.constraints.project(&mut cache.u_half_step)?;
         cache.gradient_step_u_half_step_diff_norm_sq =
             norm2_squared_diff(&cache.gradient_step, &cache.u_half_step);
+        Ok(())
     }
 
     /// Computes an LBFGS direction; updates `cache.direction_lbfgs`
@@ -201,7 +202,11 @@ where
 
         // Compute the cost at the half step
         (self.problem.cost)(&self.cache.u_half_step, &mut cost_u_half_step)?;
-        debug_assert!(matrix_operations::is_finite(&[self.cache.cost_value]));
+        if !matrix_operations::is_finite(&[self.cache.cost_value, cost_u_half_step]) {
+            return Err(SolverError::NotFiniteComputation(
+                "cost evaluation returned a non-finite value during Lipschitz estimation",
+            ));
+        }
 
         let mut it_lipschitz_search = 0;
 
@@ -217,11 +222,16 @@ where
 
             // recompute the half step...
             self.gradient_step(u_current); // updates self.cache.gradient_step
-            self.half_step(); // updates self.cache.u_half_step
+            self.half_step()?; // updates self.cache.u_half_step
 
             // recompute the cost at the half step
             // update `cost_u_half_step`
             (self.problem.cost)(&self.cache.u_half_step, &mut cost_u_half_step)?;
+            if !cost_u_half_step.is_finite() {
+                return Err(SolverError::NotFiniteComputation(
+                    "half-step cost became non-finite during Lipschitz backtracking",
+                ));
+            }
 
             // recompute the FPR and the square of its norm
             self.compute_fpr(u_current);
@@ -278,10 +288,17 @@ where
         // point `u_plus`
         (self.problem.cost)(&self.cache.u_plus, &mut self.cache.cost_value)?;
         (self.problem.gradf)(&self.cache.u_plus, &mut self.cache.gradient_u)?;
+        if !self.cache.cost_value.is_finite()
+            || !matrix_operations::is_finite(&self.cache.gradient_u)
+        {
+            return Err(SolverError::NotFiniteComputation(
+                "line-search candidate produced a non-finite cost or gradient",
+            ));
+        }
         self.cache_gradient_norm();
 
         self.gradient_step_uplus(); // gradient_step ← u_plus - gamma * gradient_u
-        self.half_step(); // u_half_step ← project(gradient_step)
+        self.half_step()?; // u_half_step ← project(gradient_step)
 
         // Update the LHS of the line search condition
         self.cache.lhs_ls = self.cache.cost_value - half * gamma * self.cache.gradient_u_norm_sq
@@ -295,9 +312,16 @@ where
         u_current.copy_from_slice(&self.cache.u_half_step); // set u_current ← u_half_step
         (self.problem.cost)(u_current, &mut self.cache.cost_value)?; // cost value
         (self.problem.gradf)(u_current, &mut self.cache.gradient_u)?; // compute gradient
+        if !self.cache.cost_value.is_finite()
+            || !matrix_operations::is_finite(&self.cache.gradient_u)
+        {
+            return Err(SolverError::NotFiniteComputation(
+                "first PANOC iterate produced a non-finite cost or gradient",
+            ));
+        }
         self.cache_gradient_norm();
         self.gradient_step(u_current); // updated self.cache.gradient_step
-        self.half_step(); // updates self.cache.u_half_step
+        self.half_step()?; // updates self.cache.u_half_step
 
         Ok(())
     }
@@ -326,6 +350,11 @@ where
     pub(crate) fn cost_value_at_best_half_step(&mut self) -> Result<T, SolverError> {
         let mut cost = T::zero();
         (self.problem.cost)(&self.cache.best_u_half_step, &mut cost)?;
+        if !cost.is_finite() {
+            return Err(SolverError::NotFiniteComputation(
+                "best cached half-step cost is non-finite",
+            ));
+        }
         Ok(cost)
     }
 }
@@ -387,12 +416,19 @@ where
         self.cache.reset();
         (self.problem.cost)(u_current, &mut self.cache.cost_value)?; // cost value
         self.estimate_loc_lip(u_current)?; // computes the gradient as well! (self.cache.gradient_u)
+        if !self.cache.cost_value.is_finite()
+            || !matrix_operations::is_finite(&self.cache.gradient_u)
+        {
+            return Err(SolverError::NotFiniteComputation(
+                "initial PANOC cost or gradient is non-finite",
+            ));
+        }
         self.cache_gradient_norm();
         self.cache.gamma =
             gamma_l_coeff::<T>() / self.cache.lipschitz_constant.max(min_l_estimate());
         self.cache.sigma = (T::one() - gamma_l_coeff::<T>()) / (cast::<T>(4.0) * self.cache.gamma);
         self.gradient_step(u_current); // updated self.cache.gradient_step
-        self.half_step(); // updates self.cache.u_half_step
+        self.half_step()?; // updates self.cache.u_half_step
 
         Ok(())
     }
@@ -504,7 +540,7 @@ mod tests {
             .gradient_step
             .copy_from_slice(&[40., 50.]);
 
-        panoc_engine.half_step(); // u_half_step ← projection(gradient_step)
+        panoc_engine.half_step().unwrap(); // u_half_step ← projection(gradient_step)
 
         unit_test_utils::assert_nearly_equal_array(
             &[0.312_347_523_777_212, 0.390_434_404_721_515],
diff --git a/src/core/panoc/panoc_optimizer.rs b/src/core/panoc/panoc_optimizer.rs
index 19f16ce2..b8525ef2 100644
--- a/src/core/panoc/panoc_optimizer.rs
+++ b/src/core/panoc/panoc_optimizer.rs
@@ -165,7 +165,9 @@ where
 
         // check for possible NaN/inf
         if !matrix_operations::is_finite(u) {
-            return Err(SolverError::NotFiniteComputation);
+            return Err(SolverError::NotFiniteComputation(
+                "final PANOC iterate contains a non-finite value",
+            ));
         }
 
         // exit status
@@ -258,7 +260,7 @@ mod tests {
         /* CHECK FEASIBILITY */
         let mut u_project = [0.0; 2];
         u_project.copy_from_slice(&u_solution);
-        bounds.project(&mut u_project);
+        bounds.project(&mut u_project).unwrap();
         unit_test_utils::assert_nearly_equal_array(
             &u_solution,
             &u_project,
@@ -301,7 +303,7 @@ mod tests {
 
         let mut u_project = [0.0_f32; 2];
         u_project.copy_from_slice(&u_solution);
-        bounds.project(&mut u_project);
+        bounds.project(&mut u_project).unwrap();
         assert!((u_solution[0] - u_project[0]).abs() < 1e-5_f32);
         assert!((u_solution[1] - u_project[1]).abs() < 1e-5_f32);
     }
diff --git a/src/core/panoc/tests.rs b/src/core/panoc/tests.rs
index 0edbed55..e94200e3 100644
--- a/src/core/panoc/tests.rs
+++ b/src/core/panoc/tests.rs
@@ -167,7 +167,7 @@ fn t_zero_gamma_l() {
     let mut panoc_engine = PANOCOptimizer::new(problem, &mut panoc_cache).with_max_iter(100);
 
     // Invoke the solver.
-    let _status = panoc_engine.solve(u);
+    panoc_engine.solve(u).unwrap();
     println!("norm_gamma_fpr = {}", panoc_cache.norm_gamma_fpr);
     println!("u = {:?}", u);
     println!("iters = {}", panoc_cache.iteration);
@@ -210,7 +210,7 @@ fn t_zero_gamma_huber() {
     let mut panoc_engine = PANOCOptimizer::new(problem, &mut panoc_cache).with_max_iter(100);
 
     // Invoke the solver.
-    let _status = panoc_engine.solve(u);
+    panoc_engine.solve(u).unwrap();
     println!("norm_gamma_fpr = {}", panoc_cache.norm_gamma_fpr);
     println!("u = {:?}", u);
     println!("iters = {}", panoc_cache.iteration);
diff --git a/src/lib.rs b/src/lib.rs
index ece5f024..d122b87a 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -40,16 +40,76 @@
 
 extern crate num;
 
-/// Exceptions/Errors that may arise while solving a problem
+use std::fmt;
+
+/// Exceptions/errors that may arise while solving a problem.
+///
+/// In the Rust API, numerical failures are generally reported through
+/// `Result<_, SolverError>`:
+///
+/// - user-provided callbacks may return [`SolverError::Cost`],
+/// - set projections may return [`SolverError::ProjectionFailed`],
+/// - non-finite intermediate values are reported via
+///   [`SolverError::NotFiniteComputation`], and
+/// - internal numerical/kernel issues are reported through
+///   [`SolverError::LinearAlgebraFailure`] or
+///   [`SolverError::InvalidProblemState`].
+///
+/// By contrast, blatant API misuse such as inconsistent slice dimensions may
+/// still panic in some low-level routines.
 #[derive(Debug, Clone, Copy, PartialEq, Eq)]
 pub enum SolverError {
     /// If the gradient or cost function cannot be evaluated
-    Cost,
+    Cost(&'static str),
     /// Computation failed and NaN/Infinite value was obtained
-    NotFiniteComputation,
+    NotFiniteComputation(&'static str),
+    /// A projection could not be computed numerically
+    ProjectionFailed(&'static str),
+    /// A linear algebra operation failed
+    LinearAlgebraFailure(&'static str),
+    /// The solver reached an unexpected internal state
+    InvalidProblemState(&'static str),
+}
+
+impl fmt::Display for SolverError {
+    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
+        match self {
+            SolverError::Cost(reason) => {
+                write!(f, "cost or gradient evaluation failed: {}", reason)
+            }
+            SolverError::NotFiniteComputation(reason) => {
+                write!(f, "non-finite computation: {}", reason)
+            }
+            SolverError::ProjectionFailed(reason) => write!(f, "projection failed: {}", reason),
+            SolverError::LinearAlgebraFailure(reason) => {
+                write!(f, "linear algebra failure: {}", reason)
+            }
+            SolverError::InvalidProblemState(reason) => {
+                write!(f, "invalid internal problem state: {}", reason)
+            }
+        }
+    }
+}
+
+impl std::error::Error for SolverError {}
+
+impl From<crate::matrix_operations::MatrixError> for SolverError {
+    fn from(_: crate::matrix_operations::MatrixError) -> Self {
+        SolverError::LinearAlgebraFailure("matrix operation failed")
+    }
+}
+
+impl From<crate::cholesky_factorizer::CholeskyError> for SolverError {
+    fn from(_: crate::cholesky_factorizer::CholeskyError) -> Self {
+        SolverError::LinearAlgebraFailure("Cholesky factorization or solve failed")
+    }
 }
 
-/// Result of a function call (status)
+/// Standard result type used by user callbacks and internal projection/codegen routines.
+///
+/// A successful call returns `Ok(())`. Failures should be reported with a
+/// descriptive [`SolverError`] so optimizers can propagate the reason to the
+/// caller.
 pub type FunctionCallResult = Result<(), SolverError>;
 
 pub mod alm;