From 45fa3e5f739124d69d792a2be649fa70af5f3915 Mon Sep 17 00:00:00 2001
From: GOELLER Adrien <adrien.goeller@rte-france.com>
Date: Mon, 22 Jun 2026 17:19:23 +0200
Subject: [PATCH 1/7] first prototype

Signed-off-by: GOELLER Adrien <adrien.goeller@rte-france.com>
---
 src/mechaphlowers/api/section_study.py        |  17 ++
 .../core/models/estimation/__init__.py        |  23 ++
 .../core/models/estimation/engine.py          | 283 ++++++++++++++++++
 .../core/models/estimation/methods.py         | 222 ++++++++++++++
 .../core/models/estimation/result.py          |  34 +++
 test/core/models/test_estimation.py           | 163 ++++++++++
 6 files changed, 742 insertions(+)
 create mode 100644 src/mechaphlowers/core/models/estimation/__init__.py
 create mode 100644 src/mechaphlowers/core/models/estimation/engine.py
 create mode 100644 src/mechaphlowers/core/models/estimation/methods.py
 create mode 100644 src/mechaphlowers/core/models/estimation/result.py
 create mode 100644 test/core/models/test_estimation.py

diff --git a/src/mechaphlowers/api/section_study.py b/src/mechaphlowers/api/section_study.py
index 6b80acb4..651b144b 100644
--- a/src/mechaphlowers/api/section_study.py
+++ b/src/mechaphlowers/api/section_study.py
@@ -31,6 +31,7 @@
 
 if TYPE_CHECKING:
     from mechaphlowers.core.models.cable.thermal import ThermalEngine
+    from mechaphlowers.core.models.estimation import EstimationEngine
     from mechaphlowers.core.models.guying import Guying
     from mechaphlowers.plotting.plot import PlotEngine
 
@@ -92,6 +93,7 @@ def __init__(
         self._plot_engine: PlotEngine | None = None
         self._thermal_engine: ThermalEngine | None = None
         self._guying: Guying | None = None
+        self._estimation_engine: EstimationEngine | None = None
         self._intermediate_memento: BalanceEngineMemento | None = None
 
     # ── Sub-engine properties ─────────────────────────────────────────────
@@ -130,6 +132,21 @@ def guying(self) -> Guying:
             self._guying = _G(self._balance_engine)
         return self._guying
 
+    @property
+    def estimation_engine(self) -> EstimationEngine:
+        """Lazy-loaded inverse estimation engine.
+
+        Returns an [`EstimationEngine`][mechaphlowers.core.models.estimation.EstimationEngine]
+        bound to this study, using Brent's method by default.
+        """
+        if self._estimation_engine is None:
+            from mechaphlowers.core.models.estimation import (
+                EstimationEngine as _EE,
+            )
+
+            self._estimation_engine = _EE(self)
+        return self._estimation_engine
+
     @property
     def intermediate_memento(self) -> BalanceEngineMemento | None:
         """The memento captured after the intermediate warm-start solve, if any."""
diff --git a/src/mechaphlowers/core/models/estimation/__init__.py b/src/mechaphlowers/core/models/estimation/__init__.py
new file mode 100644
index 00000000..a89e4e25
--- /dev/null
+++ b/src/mechaphlowers/core/models/estimation/__init__.py
@@ -0,0 +1,23 @@
+# Copyright (c) 2026, RTE (http://www.rte-france.com)
+# This Source Code Form is subject to the terms of the Mozilla Public
+# License, v. 2.0. If a copy of the MPL was not distributed with this
+# file, You can obtain one at http://mozilla.org/MPL/2.0/.
+# SPDX-License-Identifier: MPL-2.0
+
+from mechaphlowers.core.models.estimation.engine import EstimationEngine
+from mechaphlowers.core.models.estimation.methods import (
+    BisectionMethod,
+    BrentMethod,
+    NewtonMethod,
+    OptimizationMethod,
+)
+from mechaphlowers.core.models.estimation.result import EstimationResult
+
+__all__ = [
+    "EstimationEngine",
+    "EstimationResult",
+    "OptimizationMethod",
+    "BisectionMethod",
+    "BrentMethod",
+    "NewtonMethod",
+]
diff --git a/src/mechaphlowers/core/models/estimation/engine.py b/src/mechaphlowers/core/models/estimation/engine.py
new file mode 100644
index 00000000..9d2107f6
--- /dev/null
+++ b/src/mechaphlowers/core/models/estimation/engine.py
@@ -0,0 +1,283 @@
+# Copyright (c) 2026, RTE (http://www.rte-france.com)
+# This Source Code Form is subject to the terms of the Mozilla Public
+# License, v. 2.0. If a copy of the MPL was not distributed with this
+# file, You can obtain one at http://mozilla.org/MPL/2.0/.
+# SPDX-License-Identifier: MPL-2.0
+
+from __future__ import annotations
+
+import logging
+from typing import TYPE_CHECKING, Callable
+
+import numpy as np
+
+from mechaphlowers.core.models.estimation.methods import (
+    BrentMethod,
+    OptimizationMethod,
+)
+from mechaphlowers.core.models.estimation.result import EstimationResult
+
+if TYPE_CHECKING:
+    from mechaphlowers.api.section_study import SectionStudy
+
+logger = logging.getLogger(__name__)
+
+
+class EstimationEngine:
+    """Generic inverse-problem solver built on top of SectionStudy.
+
+    Wraps a `SectionStudy` and an `OptimizationMethod` to find the value of
+    a physical variable (temperature, wind, load) that produces a target
+    distance to an obstacle.
+
+    The engine saves/restores the balance-engine state around every objective
+    evaluation so that the study is left unchanged after estimation.
+
+    Args:
+        study: The `SectionStudy` instance (must have been solved via
+            `solve_adjustment` beforehand).
+        method: The optimization algorithm to use. Defaults to `BrentMethod`.
+
+    Examples:
+        >>> engine = EstimationEngine(study, method=BrentMethod(tol=0.01))
+        >>> result = engine.estimate_temperature(
+        ...     span_index=0,
+        ...     obstacle_point=np.array([150.0, 0.0, 5.0]),
+        ...     target_distance=8.0,
+        ...     bounds=(0.0, 200.0),
+        ... )
+        >>> print(result.value, result.converged)
+    """
+
+    def __init__(
+        self,
+        study: SectionStudy,
+        method: OptimizationMethod | None = None,
+    ) -> None:
+        self._study = study
+        self._method: OptimizationMethod = method or BrentMethod()
+
+    @property
+    def method(self) -> OptimizationMethod:
+        return self._method
+
+    @method.setter
+    def method(self, value: OptimizationMethod) -> None:
+        self._method = value
+
+    def estimate(
+        self,
+        objective: Callable[[float], float],
+        bounds: tuple[float, float],
+    ) -> EstimationResult:
+        """Run the optimization on a generic objective function.
+
+        The objective must be a callable ``f(x) -> float`` where the root
+        ``f(x) = 0`` corresponds to the desired solution. State management
+        (save/restore) is the caller's responsibility when using this method
+        directly.
+
+        Args:
+            objective: Function to zero. Signature: ``(x: float) -> float``.
+            bounds: ``(lower, upper)`` search interval.
+
+        Returns:
+            EstimationResult with the solution.
+        """
+        return self._method.solve(objective, bounds)
+
+    def estimate_temperature(
+        self,
+        span_index: int,
+        obstacle_point: np.ndarray,
+        target_distance: float,
+        bounds: tuple[float, float] = (0.0, 200.0),
+        wind_pressure: float | None = None,
+        ice_thickness: float | None = None,
+    ) -> EstimationResult:
+        """Find the cable temperature that yields a target distance to an obstacle.
+
+        Args:
+            span_index: Index of the span where the obstacle is located.
+            obstacle_point: 3D coordinates of the obstacle point (shape ``(3,)``).
+            target_distance: Desired distance between cable and obstacle (meters).
+            bounds: Search interval for temperature in °C.
+            wind_pressure: Fixed wind pressure in Pa (optional).
+            ice_thickness: Fixed ice thickness in m (optional).
+
+        Returns:
+            EstimationResult with the temperature value.
+        """
+
+        def objective(temperature: float) -> float:
+            return self._distance_residual(
+                span_index=span_index,
+                obstacle_point=obstacle_point,
+                target_distance=target_distance,
+                new_temperature=temperature,
+                wind_pressure=wind_pressure,
+                ice_thickness=ice_thickness,
+            )
+
+        logger.info(
+            "Estimating temperature for target distance %.3f m on span %d",
+            target_distance,
+            span_index,
+        )
+        return self._method.solve(objective, bounds)
+
+    def estimate_wind(
+        self,
+        span_index: int,
+        obstacle_point: np.ndarray,
+        target_distance: float,
+        bounds: tuple[float, float] = (0.0, 2000.0),
+        new_temperature: float | None = None,
+        ice_thickness: float | None = None,
+    ) -> EstimationResult:
+        """Find the wind pressure that yields a target distance to an obstacle.
+
+        Args:
+            span_index: Index of the span where the obstacle is located.
+            obstacle_point: 3D coordinates of the obstacle point (shape ``(3,)``).
+            target_distance: Desired distance between cable and obstacle (meters).
+            bounds: Search interval for wind pressure in Pa.
+            new_temperature: Fixed temperature in °C (optional).
+            ice_thickness: Fixed ice thickness in m (optional).
+
+        Returns:
+            EstimationResult with the wind pressure value.
+        """
+
+        def objective(wind: float) -> float:
+            return self._distance_residual(
+                span_index=span_index,
+                obstacle_point=obstacle_point,
+                target_distance=target_distance,
+                wind_pressure=wind,
+                new_temperature=new_temperature,
+                ice_thickness=ice_thickness,
+            )
+
+        logger.info(
+            "Estimating wind pressure for target distance %.3f m on span %d",
+            target_distance,
+            span_index,
+        )
+        return self._method.solve(objective, bounds)
+
+    def estimate_load(
+        self,
+        span_index: int,
+        obstacle_point: np.ndarray,
+        target_distance: float,
+        load_position_distance: float,
+        bounds: tuple[float, float] = (0.0, 100.0),
+        new_temperature: float | None = None,
+        wind_pressure: float | None = None,
+        ice_thickness: float | None = None,
+    ) -> EstimationResult:
+        """Find the load mass that yields a target distance to an obstacle.
+
+        Args:
+            span_index: Index of the span where the obstacle is located.
+            obstacle_point: 3D coordinates of the obstacle point (shape ``(3,)``).
+            target_distance: Desired distance between cable and obstacle (meters).
+            load_position_distance: Position of the load along the span (meters).
+            bounds: Search interval for load mass in kg.
+            new_temperature: Fixed temperature in °C (optional).
+            wind_pressure: Fixed wind pressure in Pa (optional).
+            ice_thickness: Fixed ice thickness in m (optional).
+
+        Returns:
+            EstimationResult with the load mass value.
+        """
+
+        def objective(load_mass: float) -> float:
+            return self._load_distance_residual(
+                span_index=span_index,
+                obstacle_point=obstacle_point,
+                target_distance=target_distance,
+                load_position_distance=load_position_distance,
+                load_mass=load_mass,
+                new_temperature=new_temperature,
+                wind_pressure=wind_pressure,
+                ice_thickness=ice_thickness,
+            )
+
+        logger.info(
+            "Estimating load mass for target distance %.3f m on span %d",
+            target_distance,
+            span_index,
+        )
+        return self._method.solve(objective, bounds)
+
+    # ── Private helpers ───────────────────────────────────────────────────
+
+    def _distance_residual(
+        self,
+        span_index: int,
+        obstacle_point: np.ndarray,
+        target_distance: float,
+        wind_pressure: float | None = None,
+        ice_thickness: float | None = None,
+        new_temperature: float | None = None,
+    ) -> float:
+        """Compute ``distance(x) - target`` with state save/restore.
+
+        Solves change-state with the given parameters, computes the distance
+        to the obstacle, then restores the engine to its original state.
+        """
+        memento = self._study.save_state()
+        try:
+            self._study.solve_change_state(
+                wind_pressure=wind_pressure,
+                ice_thickness=ice_thickness,
+                new_temperature=new_temperature,
+            )
+            distance_result = self._study.position_engine.point_distance(
+                span_index, obstacle_point
+            )
+            distance = distance_result.distance_3d
+        finally:
+            self._study.restore_state(memento)
+
+        return distance - target_distance
+
+    def _load_distance_residual(
+        self,
+        span_index: int,
+        obstacle_point: np.ndarray,
+        target_distance: float,
+        load_position_distance: float,
+        load_mass: float,
+        new_temperature: float | None = None,
+        wind_pressure: float | None = None,
+        ice_thickness: float | None = None,
+    ) -> float:
+        """Compute distance residual for a given load mass with state save/restore."""
+        memento = self._study.save_state()
+        try:
+            # Build load arrays for single point load
+            n_spans = len(
+                self._study.balance_engine.section_array.data.span_length
+            )
+            load_positions = [[] for _ in range(n_spans)]
+            load_masses = [[] for _ in range(n_spans)]
+            load_positions[span_index] = [load_position_distance]
+            load_masses[span_index] = [load_mass]
+
+            self._study.add_loads(load_positions, load_masses)
+            self._study.solve_change_state(
+                wind_pressure=wind_pressure,
+                ice_thickness=ice_thickness,
+                new_temperature=new_temperature,
+            )
+            distance_result = self._study.position_engine.point_distance(
+                span_index, obstacle_point
+            )
+            distance = distance_result.distance_3d
+        finally:
+            self._study.restore_state(memento)
+
+        return distance - target_distance
diff --git a/src/mechaphlowers/core/models/estimation/methods.py b/src/mechaphlowers/core/models/estimation/methods.py
new file mode 100644
index 00000000..3328156a
--- /dev/null
+++ b/src/mechaphlowers/core/models/estimation/methods.py
@@ -0,0 +1,222 @@
+# Copyright (c) 2026, RTE (http://www.rte-france.com)
+# This Source Code Form is subject to the terms of the Mozilla Public
+# License, v. 2.0. If a copy of the MPL was not distributed with this
+# file, You can obtain one at http://mozilla.org/MPL/2.0/.
+# SPDX-License-Identifier: MPL-2.0
+
+from __future__ import annotations
+
+import logging
+from typing import Callable, Protocol
+
+from mechaphlowers.core.models.estimation.result import EstimationResult
+
+logger = logging.getLogger(__name__)
+
+
+class OptimizationMethod(Protocol):
+    """Protocol for root-finding / optimization algorithms.
+
+    Implementations solve ``objective(x) = 0`` within given bounds.
+    """
+
+    def solve(
+        self,
+        objective: Callable[[float], float],
+        bounds: tuple[float, float],
+    ) -> EstimationResult: ...
+
+
+class BisectionMethod:
+    """Simple bisection root-finding.
+
+    Requires that ``objective`` changes sign over ``bounds``.
+
+    Args:
+        tol: Absolute tolerance on the root. Defaults to 1e-3.
+        maxiter: Maximum number of iterations. Defaults to 50.
+    """
+
+    def __init__(self, tol: float = 1e-3, maxiter: int = 50) -> None:
+        self.tol = tol
+        self.maxiter = maxiter
+
+    def solve(
+        self,
+        objective: Callable[[float], float],
+        bounds: tuple[float, float],
+    ) -> EstimationResult:
+        a, b = bounds
+        fa = objective(a)
+        fb = objective(b)
+        iterations = 2
+
+        if fa * fb > 0:
+            logger.warning(
+                "Bisection: objective does not change sign over bounds "
+                f"[{a}, {b}] (f(a)={fa:.3e}, f(b)={fb:.3e}). "
+                "Returning best bound."
+            )
+            best = a if abs(fa) < abs(fb) else b
+            residual = min(abs(fa), abs(fb))
+            return EstimationResult(
+                value=best,
+                residual=residual,
+                iterations=iterations,
+                converged=residual <= self.tol,
+            )
+
+        for _ in range(self.maxiter):
+            mid = (a + b) / 2.0
+            fmid = objective(mid)
+            iterations += 1
+
+            if abs(fmid) <= self.tol or (b - a) / 2.0 <= self.tol:
+                return EstimationResult(
+                    value=mid,
+                    residual=abs(fmid),
+                    iterations=iterations,
+                    converged=True,
+                )
+
+            if fa * fmid < 0:
+                b = mid
+            else:
+                a = mid
+                fa = fmid
+
+        mid = (a + b) / 2.0
+        return EstimationResult(
+            value=mid,
+            residual=abs(objective(mid)),
+            iterations=iterations + 1,
+            converged=False,
+        )
+
+
+class BrentMethod:
+    """Brent's method for root-finding (via scipy.optimize.brentq).
+
+    Requires that ``objective`` changes sign over ``bounds``.
+    Falls back to bisection if scipy is not available.
+
+    Args:
+        tol: Absolute tolerance on the root. Defaults to 1e-3.
+        maxiter: Maximum number of iterations. Defaults to 50.
+    """
+
+    def __init__(self, tol: float = 1e-3, maxiter: int = 50) -> None:
+        self.tol = tol
+        self.maxiter = maxiter
+
+    def solve(
+        self,
+        objective: Callable[[float], float],
+        bounds: tuple[float, float],
+    ) -> EstimationResult:
+        try:
+            from scipy.optimize import brentq
+        except ImportError:
+            logger.warning(
+                "scipy not available, falling back to BisectionMethod."
+            )
+            fallback = BisectionMethod(tol=self.tol, maxiter=self.maxiter)
+            return fallback.solve(objective, bounds)
+
+        a, b = bounds
+        iterations = 0
+
+        def counted_objective(x: float) -> float:
+            nonlocal iterations
+            iterations += 1
+            return objective(x)
+
+        try:
+            root, result_info = brentq(
+                counted_objective,
+                a,
+                b,
+                xtol=self.tol,
+                maxiter=self.maxiter,
+                full_output=True,
+            )
+            return EstimationResult(
+                value=root,
+                residual=abs(result_info.function_calls and objective(root) or 0.0),
+                iterations=result_info.iterations,
+                converged=result_info.converged,
+            )
+        except ValueError as e:
+            logger.warning(f"Brent's method failed: {e}. Trying bisection.")
+            fallback = BisectionMethod(tol=self.tol, maxiter=self.maxiter)
+            return fallback.solve(objective, bounds)
+
+
+class NewtonMethod:
+    """Newton-Raphson with finite-difference derivative.
+
+    Does not require bounds for the algorithm itself, but bounds are used
+    to clip iterates and provide the initial guess (midpoint).
+
+    Args:
+        tol: Absolute tolerance on the residual. Defaults to 1e-3.
+        maxiter: Maximum number of iterations. Defaults to 20.
+        dx: Step size for finite-difference derivative. Defaults to 1.0.
+    """
+
+    def __init__(
+        self, tol: float = 1e-3, maxiter: int = 20, dx: float = 1.0
+    ) -> None:
+        self.tol = tol
+        self.maxiter = maxiter
+        self.dx = dx
+
+    def solve(
+        self,
+        objective: Callable[[float], float],
+        bounds: tuple[float, float],
+    ) -> EstimationResult:
+        a, b = bounds
+        x = (a + b) / 2.0
+        iterations = 0
+
+        for _ in range(self.maxiter):
+            fx = objective(x)
+            iterations += 1
+
+            if abs(fx) <= self.tol:
+                return EstimationResult(
+                    value=x,
+                    residual=abs(fx),
+                    iterations=iterations,
+                    converged=True,
+                )
+
+            fx_dx = objective(x + self.dx)
+            iterations += 1
+            derivative = (fx_dx - fx) / self.dx
+
+            if abs(derivative) < 1e-12:
+                logger.warning(
+                    "Newton: near-zero derivative at x=%.6g, stopping.", x
+                )
+                return EstimationResult(
+                    value=x,
+                    residual=abs(fx),
+                    iterations=iterations,
+                    converged=False,
+                )
+
+            x_new = x - fx / derivative
+            # Clip to bounds
+            x_new = max(a, min(b, x_new))
+            x = x_new
+
+        fx = objective(x)
+        iterations += 1
+        return EstimationResult(
+            value=x,
+            residual=abs(fx),
+            iterations=iterations,
+            converged=abs(fx) <= self.tol,
+        )
diff --git a/src/mechaphlowers/core/models/estimation/result.py b/src/mechaphlowers/core/models/estimation/result.py
new file mode 100644
index 00000000..4f2b6686
--- /dev/null
+++ b/src/mechaphlowers/core/models/estimation/result.py
@@ -0,0 +1,34 @@
+# Copyright (c) 2026, RTE (http://www.rte-france.com)
+# This Source Code Form is subject to the terms of the Mozilla Public
+# License, v. 2.0. If a copy of the MPL was not distributed with this
+# file, You can obtain one at http://mozilla.org/MPL/2.0/.
+# SPDX-License-Identifier: MPL-2.0
+
+from __future__ import annotations
+
+from dataclasses import dataclass
+
+
+@dataclass(frozen=True)
+class EstimationResult:
+    """Result of an inverse estimation solve.
+
+    Attributes:
+        value: The estimated variable value that satisfies the target.
+        residual: Final residual ``|f(value) - target|``.
+        iterations: Number of objective function evaluations performed.
+        converged: Whether the algorithm converged within tolerance.
+    """
+
+    value: float
+    residual: float
+    iterations: int
+    converged: bool
+
+    def __repr__(self) -> str:
+        status = "converged" if self.converged else "NOT converged"
+        return (
+            f"EstimationResult(value={self.value:.6g}, "
+            f"residual={self.residual:.3e}, "
+            f"iterations={self.iterations}, {status})"
+        )
diff --git a/test/core/models/test_estimation.py b/test/core/models/test_estimation.py
new file mode 100644
index 00000000..0eb63861
--- /dev/null
+++ b/test/core/models/test_estimation.py
@@ -0,0 +1,163 @@
+# Copyright (c) 2026, RTE (http://www.rte-france.com)
+# This Source Code Form is subject to the terms of the Mozilla Public
+# License, v. 2.0. If a copy of the MPL was not distributed with this
+# file, You can obtain one at http://mozilla.org/MPL/2.0/.
+# SPDX-License-Identifier: MPL-2.0
+
+import numpy as np
+import pandas as pd
+import pytest
+
+from mechaphlowers.api.section_study import SectionStudy
+from mechaphlowers.core.models.balance.engine import BalanceEngine
+from mechaphlowers.core.models.estimation import (
+    BisectionMethod,
+    BrentMethod,
+    EstimationEngine,
+    EstimationResult,
+    NewtonMethod,
+)
+from mechaphlowers.entities.arrays import CableArray, SectionArray
+
+
+@pytest.fixture
+def solved_study(balance_engine_base_test: BalanceEngine) -> SectionStudy:
+    """A SectionStudy that has already been solved via solve_adjustment."""
+    study = SectionStudy(
+        cable_array=balance_engine_base_test.cable_array,
+        section_array=balance_engine_base_test.section_array,
+    )
+    study.solve_adjustment()
+    return study
+
+
+class TestEstimationResult:
+    def test_repr_converged(self):
+        r = EstimationResult(value=42.5, residual=1e-5, iterations=10, converged=True)
+        assert "converged" in repr(r)
+        assert "42.5" in repr(r)
+
+    def test_repr_not_converged(self):
+        r = EstimationResult(value=10.0, residual=5.0, iterations=50, converged=False)
+        assert "NOT converged" in repr(r)
+
+
+class TestOptimizationMethods:
+    """Test the optimization methods on a simple analytical function."""
+
+    def _simple_objective(self, x: float) -> float:
+        """f(x) = x^2 - 4, root at x=2."""
+        return x**2 - 4
+
+    def test_bisection_converges(self):
+        method = BisectionMethod(tol=1e-6, maxiter=100)
+        result = method.solve(self._simple_objective, bounds=(0.0, 5.0))
+        assert result.converged
+        assert abs(result.value - 2.0) < 1e-5
+
+    def test_brent_converges(self):
+        method = BrentMethod(tol=1e-6, maxiter=100)
+        result = method.solve(self._simple_objective, bounds=(0.0, 5.0))
+        assert result.converged
+        assert abs(result.value - 2.0) < 1e-5
+
+    def test_newton_converges(self):
+        method = NewtonMethod(tol=1e-6, maxiter=50, dx=0.01)
+        result = method.solve(self._simple_objective, bounds=(0.0, 5.0))
+        assert result.converged
+        assert abs(result.value - 2.0) < 1e-4
+
+    def test_bisection_no_sign_change(self):
+        """When bounds don't bracket a root, bisection returns best bound."""
+        method = BisectionMethod(tol=1e-6, maxiter=50)
+        # f(3)=5, f(5)=21, both positive
+        result = method.solve(self._simple_objective, bounds=(3.0, 5.0))
+        assert not result.converged or result.residual > 1e-6
+
+
+class TestEstimationEngine:
+    def test_lazy_property_access(self, solved_study: SectionStudy):
+        """estimation_engine property creates the engine lazily."""
+        engine = solved_study.estimation_engine
+        assert isinstance(engine, EstimationEngine)
+        # Second access returns same instance
+        assert solved_study.estimation_engine is engine
+
+    def test_estimate_generic(self, solved_study: SectionStudy):
+        """Generic estimate() works with a simple function."""
+        engine = EstimationEngine(solved_study, method=BisectionMethod(tol=1e-4))
+        result = engine.estimate(
+            objective=lambda x: x**2 - 9,
+            bounds=(0.0, 10.0),
+        )
+        assert result.converged
+        assert abs(result.value - 3.0) < 1e-3
+
+    def test_method_setter(self, solved_study: SectionStudy):
+        engine = EstimationEngine(solved_study)
+        assert isinstance(engine.method, BrentMethod)
+        engine.method = NewtonMethod(tol=1e-3)
+        assert isinstance(engine.method, NewtonMethod)
+
+    def test_estimate_temperature(self, solved_study: SectionStudy):
+        """Estimate temperature: solve at known temp, use that distance as target."""
+        # First, solve at a known temperature to get a reference distance
+        study = solved_study
+        study.solve_change_state(new_temperature=60.0)
+        ref_distance = study.position_engine.point_distance(
+            span_index=0, point=np.array([250.0, 10.0, 20.0])
+        ).distance_3d
+
+        # Restore state
+        study.solve_adjustment()
+
+        # Now estimate what temperature gives that distance
+        engine = EstimationEngine(study, method=BrentMethod(tol=0.1))
+        result = engine.estimate_temperature(
+            span_index=0,
+            obstacle_point=np.array([250.0, 10.0, 20.0]),
+            target_distance=ref_distance,
+            bounds=(0.0, 120.0),
+        )
+        assert result.converged
+        assert abs(result.value - 60.0) < 1.0
+
+    def test_estimate_wind(self, solved_study: SectionStudy):
+        """Estimate wind pressure using a known reference."""
+        study = solved_study
+        study.solve_change_state(wind_pressure=400.0)
+        ref_distance = study.position_engine.point_distance(
+            span_index=0, point=np.array([250.0, 10.0, 20.0])
+        ).distance_3d
+
+        # Restore state
+        study.solve_adjustment()
+
+        engine = EstimationEngine(study, method=BrentMethod(tol=0.5))
+        result = engine.estimate_wind(
+            span_index=0,
+            obstacle_point=np.array([250.0, 10.0, 20.0]),
+            target_distance=ref_distance,
+            bounds=(0.0, 1000.0),
+        )
+        assert result.converged
+        assert abs(result.value - 400.0) < 5.0
+
+    def test_state_preserved_after_estimation(self, solved_study: SectionStudy):
+        """Engine state is unchanged after estimation."""
+        study = solved_study
+        memento_before = study.save_state()
+
+        engine = EstimationEngine(study, method=BisectionMethod(tol=1.0, maxiter=5))
+        engine.estimate_temperature(
+            span_index=0,
+            obstacle_point=np.array([250.0, 10.0, 20.0]),
+            target_distance=50.0,
+            bounds=(0.0, 100.0),
+        )
+
+        # State should be unchanged
+        np.testing.assert_array_almost_equal(
+            study.balance_engine.balance_model.nodes.dxdydz,
+            memento_before.dxdydz,
+        )

From e8a322a165b997ff7f0f19dcfc558ef9c7c51dd9 Mon Sep 17 00:00:00 2001
From: lou-qui <184963772+lou-qui@users.noreply.github.com>
Date: Thu, 2 Jul 2026 16:54:33 +0200
Subject: [PATCH 2/7] Small improvements

Signed-off-by: lou-qui <184963772+lou-qui@users.noreply.github.com>
---
 .../core/models/estimation/engine.py          |  4 +--
 .../core/models/estimation/methods.py         | 20 +++++------
 test/core/models/test_estimation.py           | 35 ++++++++++++-------
 3 files changed, 34 insertions(+), 25 deletions(-)

diff --git a/src/mechaphlowers/core/models/estimation/engine.py b/src/mechaphlowers/core/models/estimation/engine.py
index 9d2107f6..b117b879 100644
--- a/src/mechaphlowers/core/models/estimation/engine.py
+++ b/src/mechaphlowers/core/models/estimation/engine.py
@@ -262,8 +262,8 @@ def _load_distance_residual(
             n_spans = len(
                 self._study.balance_engine.section_array.data.span_length
             )
-            load_positions = [[] for _ in range(n_spans)]
-            load_masses = [[] for _ in range(n_spans)]
+            load_positions: list[list[float]] = [[] for _ in range(n_spans)]
+            load_masses: list[list[float]] = [[] for _ in range(n_spans)]
             load_positions[span_index] = [load_position_distance]
             load_masses[span_index] = [load_mass]
 
diff --git a/src/mechaphlowers/core/models/estimation/methods.py b/src/mechaphlowers/core/models/estimation/methods.py
index 3328156a..348f347d 100644
--- a/src/mechaphlowers/core/models/estimation/methods.py
+++ b/src/mechaphlowers/core/models/estimation/methods.py
@@ -83,7 +83,6 @@ def solve(
                 b = mid
             else:
                 a = mid
-                fa = fmid
 
         mid = (a + b) / 2.0
         return EstimationResult(
@@ -115,7 +114,7 @@ def solve(
         bounds: tuple[float, float],
     ) -> EstimationResult:
         try:
-            from scipy.optimize import brentq
+            from scipy.optimize import brentq  # type: ignore
         except ImportError:
             logger.warning(
                 "scipy not available, falling back to BisectionMethod."
@@ -124,16 +123,10 @@ def solve(
             return fallback.solve(objective, bounds)
 
         a, b = bounds
-        iterations = 0
-
-        def counted_objective(x: float) -> float:
-            nonlocal iterations
-            iterations += 1
-            return objective(x)
 
         try:
             root, result_info = brentq(
-                counted_objective,
+                objective,
                 a,
                 b,
                 xtol=self.tol,
@@ -142,7 +135,9 @@ def counted_objective(x: float) -> float:
             )
             return EstimationResult(
                 value=root,
-                residual=abs(result_info.function_calls and objective(root) or 0.0),
+                residual=abs(
+                    result_info.function_calls and abs(objective(root)) or 0.0
+                ),
                 iterations=result_info.iterations,
                 converged=result_info.converged,
             )
@@ -165,7 +160,10 @@ class NewtonMethod:
     """
 
     def __init__(
-        self, tol: float = 1e-3, maxiter: int = 20, dx: float = 1.0
+        self,
+        tol: float = 1e-3,
+        maxiter: int = 20,
+        dx: float = 1.0,
     ) -> None:
         self.tol = tol
         self.maxiter = maxiter
diff --git a/test/core/models/test_estimation.py b/test/core/models/test_estimation.py
index 0eb63861..ecf14502 100644
--- a/test/core/models/test_estimation.py
+++ b/test/core/models/test_estimation.py
@@ -5,7 +5,6 @@
 # SPDX-License-Identifier: MPL-2.0
 
 import numpy as np
-import pandas as pd
 import pytest
 
 from mechaphlowers.api.section_study import SectionStudy
@@ -17,7 +16,6 @@
     EstimationResult,
     NewtonMethod,
 )
-from mechaphlowers.entities.arrays import CableArray, SectionArray
 
 
 @pytest.fixture
@@ -33,12 +31,16 @@ def solved_study(balance_engine_base_test: BalanceEngine) -> SectionStudy:
 
 class TestEstimationResult:
     def test_repr_converged(self):
-        r = EstimationResult(value=42.5, residual=1e-5, iterations=10, converged=True)
+        r = EstimationResult(
+            value=42.5, residual=1e-5, iterations=10, converged=True
+        )
         assert "converged" in repr(r)
         assert "42.5" in repr(r)
 
     def test_repr_not_converged(self):
-        r = EstimationResult(value=10.0, residual=5.0, iterations=50, converged=False)
+        r = EstimationResult(
+            value=10.0, residual=5.0, iterations=50, converged=False
+        )
         assert "NOT converged" in repr(r)
 
 
@@ -72,7 +74,8 @@ def test_bisection_no_sign_change(self):
         method = BisectionMethod(tol=1e-6, maxiter=50)
         # f(3)=5, f(5)=21, both positive
         result = method.solve(self._simple_objective, bounds=(3.0, 5.0))
-        assert not result.converged or result.residual > 1e-6
+        assert not result.converged
+        assert result.value == 3.0
 
 
 class TestEstimationEngine:
@@ -85,7 +88,9 @@ def test_lazy_property_access(self, solved_study: SectionStudy):
 
     def test_estimate_generic(self, solved_study: SectionStudy):
         """Generic estimate() works with a simple function."""
-        engine = EstimationEngine(solved_study, method=BisectionMethod(tol=1e-4))
+        engine = EstimationEngine(
+            solved_study, method=BisectionMethod(tol=1e-4)
+        )
         result = engine.estimate(
             objective=lambda x: x**2 - 9,
             bounds=(0.0, 10.0),
@@ -103,7 +108,8 @@ def test_estimate_temperature(self, solved_study: SectionStudy):
         """Estimate temperature: solve at known temp, use that distance as target."""
         # First, solve at a known temperature to get a reference distance
         study = solved_study
-        study.solve_change_state(new_temperature=60.0)
+        temperature = 60.0
+        study.solve_change_state(new_temperature=temperature)
         ref_distance = study.position_engine.point_distance(
             span_index=0, point=np.array([250.0, 10.0, 20.0])
         ).distance_3d
@@ -120,12 +126,13 @@ def test_estimate_temperature(self, solved_study: SectionStudy):
             bounds=(0.0, 120.0),
         )
         assert result.converged
-        assert abs(result.value - 60.0) < 1.0
+        assert abs(result.value - temperature) < 0.1
 
     def test_estimate_wind(self, solved_study: SectionStudy):
         """Estimate wind pressure using a known reference."""
         study = solved_study
-        study.solve_change_state(wind_pressure=400.0)
+        wind_pressure = 400.0
+        study.solve_change_state(wind_pressure=wind_pressure)
         ref_distance = study.position_engine.point_distance(
             span_index=0, point=np.array([250.0, 10.0, 20.0])
         ).distance_3d
@@ -141,14 +148,18 @@ def test_estimate_wind(self, solved_study: SectionStudy):
             bounds=(0.0, 1000.0),
         )
         assert result.converged
-        assert abs(result.value - 400.0) < 5.0
+        assert abs(result.value - wind_pressure) < 0.1
 
-    def test_state_preserved_after_estimation(self, solved_study: SectionStudy):
+    def test_state_preserved_after_estimation(
+        self, solved_study: SectionStudy
+    ):
         """Engine state is unchanged after estimation."""
         study = solved_study
         memento_before = study.save_state()
 
-        engine = EstimationEngine(study, method=BisectionMethod(tol=1.0, maxiter=5))
+        engine = EstimationEngine(
+            study, method=BisectionMethod(tol=1.0, maxiter=5)
+        )
         engine.estimate_temperature(
             span_index=0,
             obstacle_point=np.array([250.0, 10.0, 20.0]),

From ab7bb4fc2fcc342c01d34109bb37c7efcd7e899d Mon Sep 17 00:00:00 2001
From: lou-qui <184963772+lou-qui@users.noreply.github.com>
Date: Thu, 2 Jul 2026 17:26:37 +0200
Subject: [PATCH 3/7] tests: Add tests on section study

Signed-off-by: lou-qui <184963772+lou-qui@users.noreply.github.com>
---
 test/api/test_section_study.py | 27 +++++++++++++++++++--------
 1 file changed, 19 insertions(+), 8 deletions(-)

diff --git a/test/api/test_section_study.py b/test/api/test_section_study.py
index a4a905bd..402a6623 100644
--- a/test/api/test_section_study.py
+++ b/test/api/test_section_study.py
@@ -166,45 +166,56 @@ def test_restore_notifies_position_engine(
 
 
 class TestSectionStudyLazyEngines:
-    def test_plot_engine_not_created_eagerly(
+    def test_engines_not_created_eagerly(
         self, balance_engine_base_test: BalanceEngine
-    ):
+    ) -> None:
         study = SectionStudy(
             cable_array=balance_engine_base_test.cable_array,
             section_array=balance_engine_base_test.section_array,
         )
         assert study._plot_engine is None
+        assert study._thermal_engine is None
+        assert study._guying is None
 
     def test_plot_engine_created_on_access(
         self, balance_engine_base_test: BalanceEngine
-    ):
+    ) -> None:
         study = SectionStudy(
             cable_array=balance_engine_base_test.cable_array,
             section_array=balance_engine_base_test.section_array,
         )
-        _ = study.plot_engine
+        engine = study.plot_engine
         assert study._plot_engine is not None
 
+        # Second access returns same instance
+        assert study.plot_engine is engine
+
     def test_thermal_engine_created_on_access(
         self, balance_engine_base_test: BalanceEngine
-    ):
+    ) -> None:
         study = SectionStudy(
             cable_array=balance_engine_base_test.cable_array,
             section_array=balance_engine_base_test.section_array,
         )
-        _ = study.thermal_engine
+        engine = study.thermal_engine
         assert study._thermal_engine is not None
 
+        # Second access returns same instance
+        assert study.thermal_engine is engine
+
     def test_guying_created_on_access(
         self, balance_engine_base_test: BalanceEngine
-    ):
+    ) -> None:
         study = SectionStudy(
             cable_array=balance_engine_base_test.cable_array,
             section_array=balance_engine_base_test.section_array,
         )
-        _ = study.guying
+        guying = study.guying
         assert study._guying is not None
 
+        # Second access returns same instance
+        assert study.guying is guying
+
 
 class TestSectionStudySubEngines:
     def test_balance_engine_property(

From be793ebb7aa369b2ceffad9f82104f3a6b470019 Mon Sep 17 00:00:00 2001
From: lou-qui <184963772+lou-qui@users.noreply.github.com>
Date: Thu, 2 Jul 2026 17:32:32 +0200
Subject: [PATCH 4/7] tests: Move test on estimation engine lazy loading in
 section study

Signed-off-by: lou-qui <184963772+lou-qui@users.noreply.github.com>
---
 test/api/test_section_study.py      | 14 ++++++++++++++
 test/core/models/test_estimation.py |  7 -------
 2 files changed, 14 insertions(+), 7 deletions(-)

diff --git a/test/api/test_section_study.py b/test/api/test_section_study.py
index 402a6623..fcd5aac5 100644
--- a/test/api/test_section_study.py
+++ b/test/api/test_section_study.py
@@ -176,6 +176,7 @@ def test_engines_not_created_eagerly(
         assert study._plot_engine is None
         assert study._thermal_engine is None
         assert study._guying is None
+        assert study._estimation_engine is None
 
     def test_plot_engine_created_on_access(
         self, balance_engine_base_test: BalanceEngine
@@ -216,6 +217,19 @@ def test_guying_created_on_access(
         # Second access returns same instance
         assert study.guying is guying
 
+    def test_estimation_engine_created_on_access(
+        self, balance_engine_base_test: BalanceEngine
+    ) -> None:
+        study = SectionStudy(
+            cable_array=balance_engine_base_test.cable_array,
+            section_array=balance_engine_base_test.section_array,
+        )
+        estimation_engine = study.estimation_engine
+        assert study._estimation_engine is not None
+
+        # Second access returns same instance
+        assert study.estimation_engine is estimation_engine
+
 
 class TestSectionStudySubEngines:
     def test_balance_engine_property(
diff --git a/test/core/models/test_estimation.py b/test/core/models/test_estimation.py
index ecf14502..ecac0e5b 100644
--- a/test/core/models/test_estimation.py
+++ b/test/core/models/test_estimation.py
@@ -79,13 +79,6 @@ def test_bisection_no_sign_change(self):
 
 
 class TestEstimationEngine:
-    def test_lazy_property_access(self, solved_study: SectionStudy):
-        """estimation_engine property creates the engine lazily."""
-        engine = solved_study.estimation_engine
-        assert isinstance(engine, EstimationEngine)
-        # Second access returns same instance
-        assert solved_study.estimation_engine is engine
-
     def test_estimate_generic(self, solved_study: SectionStudy):
         """Generic estimate() works with a simple function."""
         engine = EstimationEngine(

From c37a299141b811cdae0934dccc21d72d7a4d2231 Mon Sep 17 00:00:00 2001
From: lou-qui <184963772+lou-qui@users.noreply.github.com>
Date: Thu, 2 Jul 2026 17:52:51 +0200
Subject: [PATCH 5/7] feat: Add estimation engine setter on section study

Signed-off-by: lou-qui <184963772+lou-qui@users.noreply.github.com>
---
 src/mechaphlowers/api/section_study.py | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/src/mechaphlowers/api/section_study.py b/src/mechaphlowers/api/section_study.py
index da6e604d..78ffc9f1 100644
--- a/src/mechaphlowers/api/section_study.py
+++ b/src/mechaphlowers/api/section_study.py
@@ -149,6 +149,10 @@ def estimation_engine(self) -> EstimationEngine:
             self._estimation_engine = _EE(self)
         return self._estimation_engine
 
+    @estimation_engine.setter
+    def estimation_engine(self, estimation_engine: EstimationEngine):
+        self._estimation_engine = estimation_engine
+
     @property
     def intermediate_memento(self) -> BalanceEngineMemento | None:
         """The memento captured after the intermediate warm-start solve, if any."""

From 968db5ee7bb47a57bb4b32aecde8eedd82739e3c Mon Sep 17 00:00:00 2001
From: lou-qui <184963772+lou-qui@users.noreply.github.com>
Date: Mon, 6 Jul 2026 15:29:38 +0200
Subject: [PATCH 6/7] feat(#120): Rename private methods

Signed-off-by: lou-qui <184963772+lou-qui@users.noreply.github.com>
---
 src/mechaphlowers/core/models/estimation/engine.py | 14 ++++++++------
 1 file changed, 8 insertions(+), 6 deletions(-)

diff --git a/src/mechaphlowers/core/models/estimation/engine.py b/src/mechaphlowers/core/models/estimation/engine.py
index b117b879..904e68a9 100644
--- a/src/mechaphlowers/core/models/estimation/engine.py
+++ b/src/mechaphlowers/core/models/estimation/engine.py
@@ -110,7 +110,7 @@ def estimate_temperature(
         """
 
         def objective(temperature: float) -> float:
-            return self._distance_residual(
+            return self._distance_difference(
                 span_index=span_index,
                 obstacle_point=obstacle_point,
                 target_distance=target_distance,
@@ -150,7 +150,7 @@ def estimate_wind(
         """
 
         def objective(wind: float) -> float:
-            return self._distance_residual(
+            return self._distance_difference(
                 span_index=span_index,
                 obstacle_point=obstacle_point,
                 target_distance=target_distance,
@@ -194,7 +194,7 @@ def estimate_load(
         """
 
         def objective(load_mass: float) -> float:
-            return self._load_distance_residual(
+            return self._distance_difference_with_load(
                 span_index=span_index,
                 obstacle_point=obstacle_point,
                 target_distance=target_distance,
@@ -214,7 +214,7 @@ def objective(load_mass: float) -> float:
 
     # ── Private helpers ───────────────────────────────────────────────────
 
-    def _distance_residual(
+    def _distance_difference(
         self,
         span_index: int,
         obstacle_point: np.ndarray,
@@ -244,7 +244,7 @@ def _distance_residual(
 
         return distance - target_distance
 
-    def _load_distance_residual(
+    def _distance_difference_with_load(
         self,
         span_index: int,
         obstacle_point: np.ndarray,
@@ -255,7 +255,9 @@ def _load_distance_residual(
         wind_pressure: float | None = None,
         ice_thickness: float | None = None,
     ) -> float:
-        """Compute distance residual for a given load mass with state save/restore."""
+        """Compute distance difference between target distance and distance with a given load mass.
+
+        Saves state before computiong and restores it afterwards."""
         memento = self._study.save_state()
         try:
             # Build load arrays for single point load

From 249aae92923bd999aee8d80d7b075097586d9cf3 Mon Sep 17 00:00:00 2001
From: lou-qui <184963772+lou-qui@users.noreply.github.com>
Date: Mon, 6 Jul 2026 16:15:37 +0200
Subject: [PATCH 7/7] Fixes after copilot review

Co-Authored-By: Github Copilot
Signed-off-by: lou-qui <184963772+lou-qui@users.noreply.github.com>
---
 src/mechaphlowers/core/models/estimation/engine.py | 10 +++++-----
 src/mechaphlowers/core/models/estimation/result.py |  2 +-
 2 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/src/mechaphlowers/core/models/estimation/engine.py b/src/mechaphlowers/core/models/estimation/engine.py
index 904e68a9..f44fdd85 100644
--- a/src/mechaphlowers/core/models/estimation/engine.py
+++ b/src/mechaphlowers/core/models/estimation/engine.py
@@ -264,12 +264,12 @@ def _distance_difference_with_load(
             n_spans = len(
                 self._study.balance_engine.section_array.data.span_length
             )
-            load_positions: list[list[float]] = [[] for _ in range(n_spans)]
-            load_masses: list[list[float]] = [[] for _ in range(n_spans)]
-            load_positions[span_index] = [load_position_distance]
-            load_masses[span_index] = [load_mass]
+            load_positions = np.zeros(n_spans)
+            load_masses = np.zeros(n_spans)
+            load_positions[span_index] = load_position_distance
+            load_masses[span_index] = load_mass
 
-            self._study.add_loads(load_positions, load_masses)
+            self._study.set_loads(load_positions, load_masses)
             self._study.solve_change_state(
                 wind_pressure=wind_pressure,
                 ice_thickness=ice_thickness,
diff --git a/src/mechaphlowers/core/models/estimation/result.py b/src/mechaphlowers/core/models/estimation/result.py
index 4f2b6686..14580bcf 100644
--- a/src/mechaphlowers/core/models/estimation/result.py
+++ b/src/mechaphlowers/core/models/estimation/result.py
@@ -16,7 +16,7 @@ class EstimationResult:
     Attributes:
         value: The estimated variable value that satisfies the target.
         residual: Final residual ``|f(value) - target|``.
-        iterations: Number of objective function evaluations performed.
+        iterations: Number of algorithm iterations.
         converged: Whether the algorithm converged within tolerance.
     """