chore: sync project state and current artifacts

crates/solver/src/strategies/fallback.rs

@@ -0,0 +1,490 @@
+//! Intelligent fallback solver implementation.
+//!
+//! This module provides the [`FallbackSolver`] which implements an intelligent
+//! fallback strategy between Newton-Raphson and Sequential Substitution (Picard).
+//!
+//! # Strategy
+//!
+//! The fallback solver implements the following algorithm:
+//!
+//! 1. Start with Newton-Raphson (quadratic convergence)
+//! 2. If Newton diverges, switch to Picard (more robust)
+//! 3. If Picard stabilizes (residual < threshold), try returning to Newton
+//! 4. If max switches reached, stay on current solver permanently
+//!
+//! # Features
+//!
+//! - Automatic fallback from Newton to Picard on divergence
+//! - Return to Newton when Picard stabilizes the solution
+//! - Maximum switch limit to prevent infinite oscillation
+//! - Time-budgeted solving with graceful degradation (Story 4.5)
+//! - Smart initialization support (Story 4.6)
+//! - Multi-circuit convergence criteria (Story 4.7)
+
+use std::time::{Duration, Instant};
+
+use crate::criteria::ConvergenceCriteria;
+use crate::metadata::SimulationMetadata;
+use crate::solver::{ConvergedState, ConvergenceStatus, Solver, SolverError};
+use crate::system::System;
+
+use super::{NewtonConfig, PicardConfig};
+
+/// Configuration for the intelligent fallback solver.
+///
+/// The fallback solver starts with Newton-Raphson (quadratic convergence) and
+/// automatically switches to Sequential Substitution (Picard) if Newton diverges.
+/// It can return to Newton when Picard stabilizes the solution.
+///
+/// # Example
+///
+/// ```rust
+/// use entropyk_solver::solver::{FallbackConfig, FallbackSolver, Solver};
+/// use std::time::Duration;
+///
+/// let config = FallbackConfig {
+///     fallback_enabled: true,
+///     return_to_newton_threshold: 1e-3,
+///     max_fallback_switches: 2,
+/// };
+///
+/// let solver = FallbackSolver::new(config)
+///     .with_timeout(Duration::from_secs(1));
+/// ```
+#[derive(Debug, Clone, PartialEq)]
+pub struct FallbackConfig {
+    /// Enable automatic fallback from Newton to Picard on divergence.
+    ///
+    /// When `true` (default), the solver switches to Picard if Newton diverges.
+    /// When `false`, the solver runs pure Newton or Picard without fallback.
+    pub fallback_enabled: bool,
+
+    /// Residual norm threshold for returning to Newton from Picard.
+    ///
+    /// When Picard reduces the residual below this threshold, the solver
+    /// attempts to return to Newton for faster convergence.
+    /// Default: $10^{-3}$.
+    pub return_to_newton_threshold: f64,
+
+    /// Maximum number of solver switches before staying on current solver.
+    ///
+    /// Prevents infinite oscillation between Newton and Picard.
+    /// Default: 2.
+    pub max_fallback_switches: usize,
+}
+
+impl Default for FallbackConfig {
+    fn default() -> Self {
+        Self {
+            fallback_enabled: true,
+            return_to_newton_threshold: 1e-3,
+            max_fallback_switches: 2,
+        }
+    }
+}
+
+/// Tracks which solver is currently active in the fallback loop.
+#[derive(Debug, Clone, Copy, PartialEq, Eq)]
+enum CurrentSolver {
+    Newton,
+    Picard,
+}
+
+/// Internal state for the fallback solver.
+struct FallbackState {
+    current_solver: CurrentSolver,
+    switch_count: usize,
+    /// Whether we've permanently committed to Picard (after max switches or Newton re-divergence)
+    committed_to_picard: bool,
+    /// Best state encountered across all solver invocations (Story 4.5 - AC: #4)
+    best_state: Option<Vec<f64>>,
+    /// Best residual norm across all solver invocations (Story 4.5 - AC: #4)
+    best_residual: Option<f64>,
+}
+
+impl FallbackState {
+    fn new() -> Self {
+        Self {
+            current_solver: CurrentSolver::Newton,
+            switch_count: 0,
+            committed_to_picard: false,
+            best_state: None,
+            best_residual: None,
+        }
+    }
+
+    /// Update best state if the given residual is better (Story 4.5 - AC: #4).
+    fn update_best_state(&mut self, state: &[f64], residual: f64) {
+        if self.best_residual.is_none() || residual < self.best_residual.unwrap() {
+            self.best_state = Some(state.to_vec());
+            self.best_residual = Some(residual);
+        }
+    }
+}
+
+/// Intelligent fallback solver that switches between Newton-Raphson and Picard.
+///
+/// The fallback solver implements the following algorithm:
+///
+/// 1. Start with Newton-Raphson (quadratic convergence)
+/// 2. If Newton diverges, switch to Picard (more robust)
+/// 3. If Picard stabilizes (residual < threshold), try returning to Newton
+/// 4. If max switches reached, stay on current solver permanently
+///
+/// # Timeout Handling
+///
+/// The timeout applies to the total solving time across all solver switches.
+/// Each solver inherits the remaining time budget.
+///
+/// # Pre-Allocated Buffers
+///
+/// All buffers are pre-allocated once before the fallback loop to avoid
+/// heap allocation during solver switches (NFR4).
+#[derive(Debug, Clone)]
+pub struct FallbackSolver {
+    /// Fallback behavior configuration.
+    pub config: FallbackConfig,
+    /// Newton-Raphson configuration.
+    pub newton_config: NewtonConfig,
+    /// Sequential Substitution (Picard) configuration.
+    pub picard_config: PicardConfig,
+}
+
+impl FallbackSolver {
+    /// Creates a new fallback solver with the given configuration.
+    pub fn new(config: FallbackConfig) -> Self {
+        Self {
+            config,
+            newton_config: NewtonConfig::default(),
+            picard_config: PicardConfig::default(),
+        }
+    }
+
+    /// Creates a fallback solver with default configuration.
+    pub fn default_solver() -> Self {
+        Self::new(FallbackConfig::default())
+    }
+
+    /// Sets custom Newton-Raphson configuration.
+    pub fn with_newton_config(mut self, config: NewtonConfig) -> Self {
+        self.newton_config = config;
+        self
+    }
+
+    /// Sets custom Picard configuration.
+    pub fn with_picard_config(mut self, config: PicardConfig) -> Self {
+        self.picard_config = config;
+        self
+    }
+
+    /// Sets the initial state for cold-start solving (Story 4.6 — builder pattern).
+    ///
+    /// Delegates to both `newton_config` and `picard_config` so the initial state
+    /// is used regardless of which solver is active in the fallback loop.
+    pub fn with_initial_state(mut self, state: Vec<f64>) -> Self {
+        self.newton_config.initial_state = Some(state.clone());
+        self.picard_config.initial_state = Some(state);
+        self
+    }
+
+    /// Sets multi-circuit convergence criteria (Story 4.7 — builder pattern).
+    ///
+    /// Delegates to both `newton_config` and `picard_config` so criteria are
+    /// applied regardless of which solver is active in the fallback loop.
+    pub fn with_convergence_criteria(mut self, criteria: ConvergenceCriteria) -> Self {
+        self.newton_config.convergence_criteria = Some(criteria.clone());
+        self.picard_config.convergence_criteria = Some(criteria);
+        self
+    }
+
+    /// Main fallback solving loop.
+    ///
+    /// Implements the intelligent fallback algorithm:
+    /// - Start with Newton-Raphson
+    /// - Switch to Picard on Newton divergence
+    /// - Return to Newton when Picard stabilizes (if under switch limit and residual below threshold)
+    /// - Stay on Picard permanently after max switches or if Newton re-diverges
+    fn solve_with_fallback(
+        &mut self,
+        system: &mut System,
+        start_time: Instant,
+        timeout: Option<Duration>,
+    ) -> Result<ConvergedState, SolverError> {
+        let mut state = FallbackState::new();
+
+        // Pre-configure solver configs once
+        let mut newton_cfg = self.newton_config.clone();
+        let mut picard_cfg = self.picard_config.clone();
+
+        loop {
+            // Check remaining time budget
+            let remaining = timeout.map(|t| t.saturating_sub(start_time.elapsed()));
+
+            // Check for timeout before running solver
+            if let Some(remaining_time) = remaining {
+                if remaining_time.is_zero() {
+                    return Err(SolverError::Timeout {
+                        timeout_ms: timeout.unwrap().as_millis() as u64,
+                    });
+                }
+            }
+
+            // Run current solver with remaining time
+            newton_cfg.timeout = remaining;
+            picard_cfg.timeout = remaining;
+
+            let result = match state.current_solver {
+                CurrentSolver::Newton => newton_cfg.solve(system),
+                CurrentSolver::Picard => picard_cfg.solve(system),
+            };
+
+            match result {
+                Ok(converged) => {
+                    // Update best state tracking (Story 4.5 - AC: #4)
+                    state.update_best_state(&converged.state, converged.final_residual);
+
+                    tracing::info!(
+                        solver = match state.current_solver {
+                            CurrentSolver::Newton => "NewtonRaphson",
+                            CurrentSolver::Picard => "Picard",
+                        },
+                        iterations = converged.iterations,
+                        final_residual = converged.final_residual,
+                        switch_count = state.switch_count,
+                        "Fallback solver converged"
+                    );
+                    return Ok(converged);
+                }
+                Err(SolverError::Timeout { timeout_ms }) => {
+                    // Story 4.5 - AC: #4: Return best state on timeout if available
+                    if let (Some(best_state), Some(best_residual)) =
+                        (state.best_state.clone(), state.best_residual)
+                    {
+                        tracing::info!(
+                            best_residual = best_residual,
+                            "Returning best state across all solver invocations on timeout"
+                        );
+                        return Ok(ConvergedState::new(
+                            best_state,
+                            0, // iterations not tracked across switches
+                            best_residual,
+                            ConvergenceStatus::TimedOutWithBestState,
+                            SimulationMetadata::new(system.input_hash()),
+                        ));
+                    }
+                    return Err(SolverError::Timeout { timeout_ms });
+                }
+                Err(SolverError::Divergence { ref reason }) => {
+                    // Handle divergence based on current solver and state
+                    if !self.config.fallback_enabled {
+                        tracing::info!(
+                            solver = match state.current_solver {
+                                CurrentSolver::Newton => "NewtonRaphson",
+                                CurrentSolver::Picard => "Picard",
+                            },
+                            reason = reason,
+                            "Divergence detected, fallback disabled"
+                        );
+                        return result;
+                    }
+
+                    match state.current_solver {
+                        CurrentSolver::Newton => {
+                            // Newton diverged - switch to Picard (stay there permanently after max switches)
+                            if state.switch_count >= self.config.max_fallback_switches {
+                                // Max switches reached - commit to Picard permanently
+                                state.committed_to_picard = true;
+                                state.current_solver = CurrentSolver::Picard;
+                                tracing::info!(
+                                    switch_count = state.switch_count,
+                                    max_switches = self.config.max_fallback_switches,
+                                    "Max switches reached, committing to Picard permanently"
+                                );
+                            } else {
+                                // Switch to Picard
+                                state.switch_count += 1;
+                                state.current_solver = CurrentSolver::Picard;
+                                tracing::warn!(
+                                    switch_count = state.switch_count,
+                                    reason = reason,
+                                    "Newton diverged, switching to Picard"
+                                );
+                            }
+                            // Continue loop with Picard
+                        }
+                        CurrentSolver::Picard => {
+                            // Picard diverged - if we were trying Newton again, commit to Picard permanently
+                            if state.switch_count > 0 && !state.committed_to_picard {
+                                state.committed_to_picard = true;
+                                tracing::info!(
+                                    switch_count = state.switch_count,
+                                    reason = reason,
+                                    "Newton re-diverged after return from Picard, staying on Picard permanently"
+                                );
+                                // Stay on Picard and try again
+                            } else {
+                                // Picard diverged with no return attempt - no more fallbacks available
+                                tracing::warn!(
+                                    reason = reason,
+                                    "Picard diverged, no more fallbacks available"
+                                );
+                                return result;
+                            }
+                        }
+                    }
+                }
+                Err(SolverError::NonConvergence {
+                    iterations,
+                    final_residual,
+                }) => {
+                    // Non-convergence: check if we should try the other solver
+                    if !self.config.fallback_enabled {
+                        return Err(SolverError::NonConvergence {
+                            iterations,
+                            final_residual,
+                        });
+                    }
+
+                    match state.current_solver {
+                        CurrentSolver::Newton => {
+                            // Newton didn't converge - try Picard
+                            if state.switch_count >= self.config.max_fallback_switches {
+                                // Max switches reached - commit to Picard permanently
+                                state.committed_to_picard = true;
+                                state.current_solver = CurrentSolver::Picard;
+                                tracing::info!(
+                                    switch_count = state.switch_count,
+                                    "Max switches reached, committing to Picard permanently"
+                                );
+                            } else {
+                                state.switch_count += 1;
+                                state.current_solver = CurrentSolver::Picard;
+                                tracing::info!(
+                                    switch_count = state.switch_count,
+                                    iterations = iterations,
+                                    final_residual = final_residual,
+                                    "Newton did not converge, switching to Picard"
+                                );
+                            }
+                            // Continue loop with Picard
+                        }
+                        CurrentSolver::Picard => {
+                            // Picard didn't converge - check if we should try Newton
+                            if state.committed_to_picard
+                                || state.switch_count >= self.config.max_fallback_switches
+                            {
+                                tracing::info!(
+                                    iterations = iterations,
+                                    final_residual = final_residual,
+                                    "Picard did not converge, no more fallbacks"
+                                );
+                                return Err(SolverError::NonConvergence {
+                                    iterations,
+                                    final_residual,
+                                });
+                            }
+
+                            // Check if residual is low enough to try Newton
+                            if final_residual < self.config.return_to_newton_threshold {
+                                state.switch_count += 1;
+                                state.current_solver = CurrentSolver::Newton;
+                                tracing::info!(
+                                    switch_count = state.switch_count,
+                                    final_residual = final_residual,
+                                    threshold = self.config.return_to_newton_threshold,
+                                    "Picard stabilized, attempting Newton return"
+                                );
+                                // Continue loop with Newton
+                            } else {
+                                // Stay on Picard and keep trying
+                                tracing::debug!(
+                                    final_residual = final_residual,
+                                    threshold = self.config.return_to_newton_threshold,
+                                    "Picard not yet stabilized, aborting"
+                                );
+                                return Err(SolverError::NonConvergence {
+                                    iterations,
+                                    final_residual,
+                                });
+                            }
+                        }
+                    }
+                }
+                Err(other) => {
+                    // InvalidSystem or other errors - propagate immediately
+                    return Err(other);
+                }
+            }
+        }
+    }
+}
+
+impl Solver for FallbackSolver {
+    fn solve(&mut self, system: &mut System) -> Result<ConvergedState, SolverError> {
+        let start_time = Instant::now();
+        let timeout = self.newton_config.timeout.or(self.picard_config.timeout);
+
+        tracing::info!(
+            fallback_enabled = self.config.fallback_enabled,
+            return_to_newton_threshold = self.config.return_to_newton_threshold,
+            max_fallback_switches = self.config.max_fallback_switches,
+            "Fallback solver starting"
+        );
+
+        if self.config.fallback_enabled {
+            self.solve_with_fallback(system, start_time, timeout)
+        } else {
+            // Fallback disabled - run pure Newton
+            self.newton_config.solve(system)
+        }
+    }
+
+    fn with_timeout(mut self, timeout: Duration) -> Self {
+        self.newton_config.timeout = Some(timeout);
+        self.picard_config.timeout = Some(timeout);
+        self
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::solver::Solver;
+    use crate::system::System;
+    use std::time::Duration;
+
+    #[test]
+    fn test_fallback_config_defaults() {
+        let cfg = FallbackConfig::default();
+        assert!(cfg.fallback_enabled);
+        assert!((cfg.return_to_newton_threshold - 1e-3).abs() < 1e-15);
+        assert_eq!(cfg.max_fallback_switches, 2);
+    }
+
+    #[test]
+    fn test_fallback_solver_new() {
+        let config = FallbackConfig {
+            fallback_enabled: false,
+            return_to_newton_threshold: 5e-4,
+            max_fallback_switches: 3,
+        };
+        let solver = FallbackSolver::new(config.clone());
+        assert_eq!(solver.config, config);
+    }
+
+    #[test]
+    fn test_fallback_solver_with_timeout() {
+        let timeout = Duration::from_millis(500);
+        let solver = FallbackSolver::default_solver().with_timeout(timeout);
+        assert_eq!(solver.newton_config.timeout, Some(timeout));
+        assert_eq!(solver.picard_config.timeout, Some(timeout));
+    }
+
+    #[test]
+    fn test_fallback_solver_trait_object() {
+        let mut boxed: Box<dyn Solver> = Box::new(FallbackSolver::default_solver());
+        let mut system = System::new();
+        system.finalize().unwrap();
+        assert!(boxed.solve(&mut system).is_err());
+    }
+}

crates/solver/src/strategies/mod.rs

@@ -0,0 +1,232 @@
+//! Solver strategy implementations for thermodynamic system solving.
+//!
+//! This module provides the concrete solver implementations that can be used
+//! via the [`Solver`] trait or the [`SolverStrategy`] enum for zero-cost
+//! static dispatch.
+//!
+//! # Available Strategies
+//!
+//! - [`NewtonRaphson`] — Newton-Raphson solver with quadratic convergence
+//! - [`SequentialSubstitution`] — Picard iteration solver, more robust for non-linear systems
+//! - [`FallbackSolver`] — Intelligent fallback between Newton and Picard
+//!
+//! # Example
+//!
+//! ```rust
+//! use entropyk_solver::solver::{Solver, SolverStrategy};
+//! use std::time::Duration;
+//!
+//! let solver = SolverStrategy::default()
+//!     .with_timeout(Duration::from_millis(500));
+//! ```
+
+mod fallback;
+mod newton_raphson;
+mod sequential_substitution;
+
+pub use fallback::{FallbackConfig, FallbackSolver};
+pub use newton_raphson::NewtonConfig;
+pub use sequential_substitution::PicardConfig;
+
+use crate::solver::{ConvergedState, Solver, SolverError};
+use crate::system::System;
+use std::time::Duration;
+
+/// Enum-based solver strategy dispatcher.
+///
+/// Provides zero-cost static dispatch to the selected solver strategy via
+/// `match` (monomorphization), avoiding vtable overhead while still allowing
+/// runtime strategy selection.
+///
+/// # Default
+///
+/// `SolverStrategy::default()` returns `NewtonRaphson(NewtonConfig::default())`.
+///
+/// # Example
+///
+/// ```rust
+/// use entropyk_solver::solver::{Solver, SolverStrategy, PicardConfig};
+/// use std::time::Duration;
+///
+/// let strategy = SolverStrategy::SequentialSubstitution(
+///     PicardConfig { relaxation_factor: 0.3, ..Default::default() }
+/// ).with_timeout(Duration::from_secs(1));
+/// ```
+#[derive(Debug, Clone, PartialEq)]
+pub enum SolverStrategy {
+    /// Newton-Raphson solver (quadratic convergence, requires Jacobian).
+    NewtonRaphson(NewtonConfig),
+    /// Sequential Substitution / Picard iteration (robust, no Jacobian needed).
+    SequentialSubstitution(PicardConfig),
+}
+
+impl Default for SolverStrategy {
+    /// Returns `SolverStrategy::NewtonRaphson(NewtonConfig::default())`.
+    fn default() -> Self {
+        SolverStrategy::NewtonRaphson(NewtonConfig::default())
+    }
+}
+
+impl Solver for SolverStrategy {
+    fn solve(&mut self, system: &mut System) -> Result<ConvergedState, SolverError> {
+        tracing::info!(
+            strategy = match self {
+                SolverStrategy::NewtonRaphson(_) => "NewtonRaphson",
+                SolverStrategy::SequentialSubstitution(_) => "SequentialSubstitution",
+            },
+            "SolverStrategy::solve dispatching"
+        );
+        let result = match self {
+            SolverStrategy::NewtonRaphson(cfg) => cfg.solve(system),
+            SolverStrategy::SequentialSubstitution(cfg) => cfg.solve(system),
+        };
+
+        if let Ok(state) = &result {
+            if state.is_converged() {
+                // Post-solve validation checks
+                // Convert Vec<f64> to SystemState for validation methods
+                let system_state: entropyk_components::SystemState = state.state.clone().into();
+                system.check_mass_balance(&system_state)?;
+                system.check_energy_balance(&system_state)?;
+            }
+        }
+
+        result
+    }
+
+    fn with_timeout(self, timeout: Duration) -> Self {
+        match self {
+            SolverStrategy::NewtonRaphson(cfg) => {
+                SolverStrategy::NewtonRaphson(cfg.with_timeout(timeout))
+            }
+            SolverStrategy::SequentialSubstitution(cfg) => {
+                SolverStrategy::SequentialSubstitution(cfg.with_timeout(timeout))
+            }
+        }
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::system::System;
+    use std::time::Duration;
+
+    /// Verify that `SolverStrategy::default()` returns Newton-Raphson.
+    #[test]
+    fn test_solver_strategy_default_is_newton_raphson() {
+        let strategy = SolverStrategy::default();
+        assert!(
+            matches!(strategy, SolverStrategy::NewtonRaphson(_)),
+            "Default strategy must be NewtonRaphson, got {:?}",
+            strategy
+        );
+    }
+
+    /// Verify that the Newton-Raphson variant wraps a `NewtonConfig`.
+    #[test]
+    fn test_solver_strategy_newton_raphson_variant() {
+        let strategy = SolverStrategy::NewtonRaphson(NewtonConfig::default());
+        match strategy {
+            SolverStrategy::NewtonRaphson(cfg) => {
+                assert_eq!(cfg.max_iterations, 100);
+                assert!((cfg.tolerance - 1e-6).abs() < 1e-15);
+                assert!(!cfg.line_search);
+                assert!(cfg.timeout.is_none());
+            }
+            other => panic!("Expected NewtonRaphson, got {:?}", other),
+        }
+    }
+
+    /// Verify that the Sequential Substitution variant wraps a `PicardConfig`.
+    #[test]
+    fn test_solver_strategy_sequential_substitution_variant() {
+        let strategy = SolverStrategy::SequentialSubstitution(PicardConfig::default());
+        match strategy {
+            SolverStrategy::SequentialSubstitution(cfg) => {
+                assert_eq!(cfg.max_iterations, 100);
+                assert!((cfg.tolerance - 1e-6).abs() < 1e-15);
+                assert!((cfg.relaxation_factor - 0.5).abs() < 1e-15);
+                assert!(cfg.timeout.is_none());
+            }
+            other => panic!("Expected SequentialSubstitution, got {:?}", other),
+        }
+    }
+
+    /// Verify that `with_timeout` on `SolverStrategy::NewtonRaphson` propagates to inner config.
+    #[test]
+    fn test_solver_strategy_newton_with_timeout() {
+        let timeout = Duration::from_millis(500);
+        let strategy = SolverStrategy::default().with_timeout(timeout);
+        match strategy {
+            SolverStrategy::NewtonRaphson(cfg) => {
+                assert_eq!(cfg.timeout, Some(timeout));
+            }
+            other => panic!("Expected NewtonRaphson after with_timeout, got {:?}", other),
+        }
+    }
+
+    /// Verify that `with_timeout` on `SolverStrategy::SequentialSubstitution` propagates.
+    #[test]
+    fn test_solver_strategy_picard_with_timeout() {
+        let timeout = Duration::from_secs(1);
+        let strategy =
+            SolverStrategy::SequentialSubstitution(PicardConfig::default()).with_timeout(timeout);
+        match strategy {
+            SolverStrategy::SequentialSubstitution(cfg) => {
+                assert_eq!(cfg.timeout, Some(timeout));
+            }
+            other => panic!(
+                "Expected SequentialSubstitution after with_timeout, got {:?}",
+                other
+            ),
+        }
+    }
+
+    /// Verify that `SolverStrategy::NewtonRaphson` dispatches to the Newton implementation.
+    #[test]
+    fn test_solver_strategy_newton_dispatch_reaches_stub() {
+        let mut strategy = SolverStrategy::default(); // NewtonRaphson
+        let mut system = System::new();
+        system.finalize().unwrap();
+        let result = strategy.solve(&mut system);
+        // Empty system should return InvalidSystem
+        assert!(
+            result.is_err(),
+            "Newton solver must return Err for empty system"
+        );
+        match result {
+            Err(SolverError::InvalidSystem { ref message }) => {
+                assert!(
+                    message.contains("Empty") || message.contains("no state"),
+                    "Newton dispatch must detect empty system, got: {}",
+                    message
+                );
+            }
+            other => panic!("Expected InvalidSystem from Newton solver, got {:?}", other),
+        }
+    }
+
+    /// Verify that `SolverStrategy::SequentialSubstitution` dispatches to the Picard implementation.
+    #[test]
+    fn test_solver_strategy_picard_dispatch_reaches_implementation() {
+        let mut strategy = SolverStrategy::SequentialSubstitution(PicardConfig::default());
+        let mut system = System::new();
+        system.finalize().unwrap();
+        let result = strategy.solve(&mut system);
+        assert!(
+            result.is_err(),
+            "Picard solver must return Err for empty system"
+        );
+        match result {
+            Err(SolverError::InvalidSystem { ref message }) => {
+                assert!(
+                    message.contains("Empty") || message.contains("no state"),
+                    "Picard dispatch must detect empty system, got: {}",
+                    message
+                );
+            }
+            other => panic!("Expected InvalidSystem from Picard solver, got {:?}", other),
+        }
+    }
+}

crates/solver/src/strategies/newton_raphson.rs

@@ -0,0 +1,491 @@
+//! Newton-Raphson solver implementation.
+//!
+//! Provides [`NewtonConfig`] which implements the Newton-Raphson method for
+//! solving systems of non-linear equations with quadratic convergence.
+
+use std::time::{Duration, Instant};
+
+use crate::criteria::ConvergenceCriteria;
+use crate::jacobian::JacobianMatrix;
+use crate::metadata::SimulationMetadata;
+use crate::solver::{
+    apply_newton_step, ConvergedState, ConvergenceStatus, JacobianFreezingConfig, Solver,
+    SolverError, TimeoutConfig,
+};
+use crate::system::System;
+use entropyk_components::JacobianBuilder;
+
+/// Configuration for the Newton-Raphson solver.
+///
+/// Solves F(x) = 0 by iterating: x_{k+1} = x_k - α·J^{-1}·r(x_k)
+/// where J is the Jacobian matrix and α is the step length.
+#[derive(Debug, Clone, PartialEq)]
+pub struct NewtonConfig {
+    /// Maximum iterations before declaring non-convergence. Default: 100.
+    pub max_iterations: usize,
+    /// Convergence tolerance (L2 norm). Default: 1e-6.
+    pub tolerance: f64,
+    /// Enable Armijo line-search. Default: false.
+    pub line_search: bool,
+    /// Optional time budget.
+    pub timeout: Option<Duration>,
+    /// Use numerical Jacobian (finite differences). Default: false.
+    pub use_numerical_jacobian: bool,
+    /// Armijo condition constant. Default: 1e-4.
+    pub line_search_armijo_c: f64,
+    /// Max backtracking iterations. Default: 20.
+    pub line_search_max_backtracks: usize,
+    /// Divergence threshold. Default: 1e10.
+    pub divergence_threshold: f64,
+    /// Timeout behavior configuration.
+    pub timeout_config: TimeoutConfig,
+    /// Previous state for ZOH fallback.
+    pub previous_state: Option<Vec<f64>>,
+    /// Residual for previous_state.
+    pub previous_residual: Option<f64>,
+    /// Smart initial state for cold-start.
+    pub initial_state: Option<Vec<f64>>,
+    /// Multi-circuit convergence criteria.
+    pub convergence_criteria: Option<ConvergenceCriteria>,
+    /// Jacobian-freezing optimization.
+    pub jacobian_freezing: Option<JacobianFreezingConfig>,
+}
+
+impl Default for NewtonConfig {
+    fn default() -> Self {
+        Self {
+            max_iterations: 100,
+            tolerance: 1e-6,
+            line_search: false,
+            timeout: None,
+            use_numerical_jacobian: false,
+            line_search_armijo_c: 1e-4,
+            line_search_max_backtracks: 20,
+            divergence_threshold: 1e10,
+            timeout_config: TimeoutConfig::default(),
+            previous_state: None,
+            previous_residual: None,
+            initial_state: None,
+            convergence_criteria: None,
+            jacobian_freezing: None,
+        }
+    }
+}
+
+impl NewtonConfig {
+    /// Sets the initial state for cold-start solving.
+    pub fn with_initial_state(mut self, state: Vec<f64>) -> Self {
+        self.initial_state = Some(state);
+        self
+    }
+
+    /// Sets multi-circuit convergence criteria.
+    pub fn with_convergence_criteria(mut self, criteria: ConvergenceCriteria) -> Self {
+        self.convergence_criteria = Some(criteria);
+        self
+    }
+
+    /// Enables Jacobian-freezing optimization.
+    pub fn with_jacobian_freezing(mut self, config: JacobianFreezingConfig) -> Self {
+        self.jacobian_freezing = Some(config);
+        self
+    }
+
+    /// Computes the L2 norm of the residual vector.
+    fn residual_norm(residuals: &[f64]) -> f64 {
+        residuals.iter().map(|r| r * r).sum::<f64>().sqrt()
+    }
+
+    /// Handles timeout based on configuration.
+    fn handle_timeout(
+        &self,
+        best_state: &[f64],
+        best_residual: f64,
+        iterations: usize,
+        timeout: Duration,
+        system: &System,
+    ) -> Result<ConvergedState, SolverError> {
+        if !self.timeout_config.return_best_state_on_timeout {
+            return Err(SolverError::Timeout {
+                timeout_ms: timeout.as_millis() as u64,
+            });
+        }
+
+        if self.timeout_config.zoh_fallback {
+            if let Some(ref prev_state) = self.previous_state {
+                let residual = self.previous_residual.unwrap_or(best_residual);
+                tracing::info!(iterations, residual, "ZOH fallback");
+                return Ok(ConvergedState::new(
+                    prev_state.clone(),
+                    iterations,
+                    residual,
+                    ConvergenceStatus::TimedOutWithBestState,
+                    SimulationMetadata::new(system.input_hash()),
+                ));
+            }
+        }
+
+        tracing::info!(iterations, best_residual, "Returning best state on timeout");
+        Ok(ConvergedState::new(
+            best_state.to_vec(),
+            iterations,
+            best_residual,
+            ConvergenceStatus::TimedOutWithBestState,
+            SimulationMetadata::new(system.input_hash()),
+        ))
+    }
+
+    /// Checks for divergence based on residual growth.
+    fn check_divergence(
+        &self,
+        current_norm: f64,
+        previous_norm: f64,
+        divergence_count: &mut usize,
+    ) -> Option<SolverError> {
+        if current_norm > self.divergence_threshold {
+            return Some(SolverError::Divergence {
+                reason: format!("Residual {} exceeds threshold {}", current_norm, self.divergence_threshold),
+            });
+        }
+
+        if current_norm > previous_norm {
+            *divergence_count += 1;
+            if *divergence_count >= 3 {
+                return Some(SolverError::Divergence {
+                    reason: format!("Residual increased 3x: {:.6e} → {:.6e}", previous_norm, current_norm),
+                });
+            }
+        } else {
+            *divergence_count = 0;
+        }
+
+        None
+    }
+
+    /// Performs Armijo line search. Returns Some(alpha) if valid step found.
+    /// hot path. `state_copy` and `new_residuals` must have appropriate lengths.
+    #[allow(clippy::too_many_arguments)]
+    fn line_search(
+        &self,
+        system: &System,
+        state: &mut Vec<f64>,
+        delta: &[f64],
+        _residuals: &[f64],
+        current_norm: f64,
+        state_copy: &mut [f64],
+        new_residuals: &mut Vec<f64>,
+        clipping_mask: &[Option<(f64, f64)>],
+    ) -> Option<f64> {
+        let mut alpha: f64 = 1.0;
+        state_copy.copy_from_slice(state);
+        let gradient_dot_delta = -current_norm;
+
+        for _backtrack in 0..self.line_search_max_backtracks {
+            apply_newton_step(state, delta, clipping_mask, alpha);
+
+            if system.compute_residuals(state, new_residuals).is_err() {
+                state.copy_from_slice(state_copy);
+                alpha *= 0.5;
+                continue;
+            }
+
+            let new_norm = Self::residual_norm(new_residuals);
+            if new_norm <= current_norm + self.line_search_armijo_c * alpha * gradient_dot_delta {
+                tracing::debug!(alpha, old_norm = current_norm, new_norm, "Line search accepted");
+                return Some(alpha);
+            }
+
+            state.copy_from_slice(state_copy);
+            alpha *= 0.5;
+        }
+
+        tracing::warn!("Line search failed after {} backtracks", self.line_search_max_backtracks);
+        None
+    }
+}
+
+impl Solver for NewtonConfig {
+    fn solve(&mut self, system: &mut System) -> Result<ConvergedState, SolverError> {
+        let start_time = Instant::now();
+
+        tracing::info!(
+            max_iterations = self.max_iterations,
+            tolerance = self.tolerance,
+            line_search = self.line_search,
+            "Newton-Raphson solver starting"
+        );
+
+        let n_state = system.full_state_vector_len();
+        let n_equations: usize = system
+            .traverse_for_jacobian()
+            .map(|(_, c, _)| c.n_equations())
+            .sum::<usize>()
+            + system.constraints().count()
+            + system.coupling_residual_count();
+
+        if n_state == 0 || n_equations == 0 {
+            return Err(SolverError::InvalidSystem {
+                message: "Empty system has no state variables or equations".to_string(),
+            });
+        }
+
+        // Pre-allocate all buffers
+        let mut state: Vec<f64> = self
+            .initial_state
+            .as_ref()
+            .map(|s| {
+                debug_assert_eq!(s.len(), n_state, "initial_state length mismatch");
+                if s.len() == n_state { s.clone() } else { vec![0.0; n_state] }
+            })
+            .unwrap_or_else(|| vec![0.0; n_state]);
+        let mut residuals: Vec<f64> = vec![0.0; n_equations];
+        let mut jacobian_builder = JacobianBuilder::new();
+        let mut divergence_count: usize = 0;
+        let mut previous_norm: f64;
+        let mut state_copy: Vec<f64> = vec![0.0; n_state]; // Pre-allocated for line search
+        let mut new_residuals: Vec<f64> = vec![0.0; n_equations]; // Pre-allocated for line search
+        let mut prev_iteration_state: Vec<f64> = vec![0.0; n_state]; // For convergence delta check
+
+        // Pre-allocate best-state tracking buffer (Story 4.5 - AC: #5)
+        let mut best_state: Vec<f64> = vec![0.0; n_state];
+        let mut best_residual: f64;
+
+        // Jacobian-freezing tracking state
+        let mut jacobian_matrix = JacobianMatrix::zeros(n_equations, n_state);
+        let mut frozen_count: usize = 0;
+        let mut force_recompute: bool = true;
+
+        // Pre-compute clipping mask
+        let clipping_mask: Vec<Option<(f64, f64)>> = (0..n_state)
+            .map(|i| system.get_bounds_for_state_index(i))
+            .collect();
+
+        // Initial residual computation
+        system
+            .compute_residuals(&state, &mut residuals)
+            .map_err(|e| SolverError::InvalidSystem {
+                message: format!("Failed to compute initial residuals: {:?}", e),
+            })?;
+
+        let mut current_norm = Self::residual_norm(&residuals);
+        best_state.copy_from_slice(&state);
+        best_residual = current_norm;
+
+        tracing::debug!(iteration = 0, residual_norm = current_norm, "Initial state");
+
+        // Check if already converged
+        if current_norm < self.tolerance {
+            let status = if !system.saturated_variables().is_empty() {
+                ConvergenceStatus::ControlSaturation
+            } else {
+                ConvergenceStatus::Converged
+            };
+
+            if let Some(ref criteria) = self.convergence_criteria {
+                let report = criteria.check(&state, None, &residuals, system);
+                if report.is_globally_converged() {
+                    tracing::info!(iterations = 0, final_residual = current_norm, "Converged at initial state (criteria)");
+                    return Ok(ConvergedState::with_report(
+                        state, 0, current_norm, status, report, SimulationMetadata::new(system.input_hash()),
+                    ));
+                }
+            } else {
+                tracing::info!(iterations = 0, final_residual = current_norm, "Converged at initial state");
+                return Ok(ConvergedState::new(
+                    state, 0, current_norm, status, SimulationMetadata::new(system.input_hash()),
+                ));
+            }
+        }
+
+        // Main Newton-Raphson iteration loop
+        for iteration in 1..=self.max_iterations {
+            prev_iteration_state.copy_from_slice(&state);
+
+            // Check timeout
+            if let Some(timeout) = self.timeout {
+                if start_time.elapsed() > timeout {
+                    tracing::info!(iteration, elapsed_ms = ?start_time.elapsed(), best_residual, "Solver timed out");
+                    return self.handle_timeout(&best_state, best_residual, iteration - 1, timeout, system);
+                }
+            }
+
+            // Jacobian Assembly / Freeze Decision
+            let should_recompute = if let Some(ref freeze_cfg) = self.jacobian_freezing {
+                if force_recompute {
+                    true
+                } else if frozen_count >= freeze_cfg.max_frozen_iters {
+                    tracing::debug!(iteration, frozen_count, "Jacobian freeze limit reached");
+                    true
+                } else {
+                    false
+                }
+            } else {
+                true
+            };
+
+            if should_recompute {
+                // Fresh Jacobian assembly (in-place update)
+                jacobian_builder.clear();
+                if self.use_numerical_jacobian {
+                    // Numerical Jacobian via finite differences
+                    let compute_residuals_fn = |s: &[f64], r: &mut [f64]| {
+                        let s_vec = s.to_vec();
+                        let mut r_vec = vec![0.0; r.len()];
+                        let result = system.compute_residuals(&s_vec, &mut r_vec);
+                        r.copy_from_slice(&r_vec);
+                        result.map(|_| ()).map_err(|e| format!("{:?}", e))
+                    };
+                    let jm = JacobianMatrix::numerical(compute_residuals_fn, &state, &residuals, 1e-5)
+                        .map_err(|e| SolverError::InvalidSystem {
+                            message: format!("Failed to compute numerical Jacobian: {}", e),
+                        })?;
+                    jacobian_matrix.as_matrix_mut().copy_from(jm.as_matrix());
+                } else {
+                    system.assemble_jacobian(&state, &mut jacobian_builder)
+                        .map_err(|e| SolverError::InvalidSystem {
+                            message: format!("Failed to assemble Jacobian: {:?}", e),
+                        })?;
+                    jacobian_matrix.update_from_builder(jacobian_builder.entries());
+                };
+
+                frozen_count = 0;
+                force_recompute = false;
+                tracing::debug!(iteration, "Fresh Jacobian computed");
+            } else {
+                frozen_count += 1;
+                tracing::debug!(iteration, frozen_count, "Reusing frozen Jacobian");
+            }
+
+            // Solve J·Δx = -r
+            let delta = match jacobian_matrix.solve(&residuals) {
+                Some(d) => d,
+                None => {
+                    return Err(SolverError::Divergence {
+                        reason: "Jacobian is singular".to_string(),
+                    });
+                }
+            };
+
+            // Apply step with optional line search
+            let alpha = if self.line_search {
+                match self.line_search(
+                    system, &mut state, &delta, &residuals, current_norm,
+                    &mut state_copy, &mut new_residuals, &clipping_mask,
+                ) {
+                    Some(a) => a,
+                    None => {
+                        return Err(SolverError::Divergence {
+                            reason: "Line search failed".to_string(),
+                        });
+                    }
+                }
+            } else {
+                apply_newton_step(&mut state, &delta, &clipping_mask, 1.0);
+                1.0
+            };
+
+            system.compute_residuals(&state, &mut residuals)
+                .map_err(|e| SolverError::InvalidSystem {
+                    message: format!("Failed to compute residuals: {:?}", e),
+                })?;
+
+            previous_norm = current_norm;
+            current_norm = Self::residual_norm(&residuals);
+
+            if current_norm < best_residual {
+                best_state.copy_from_slice(&state);
+                best_residual = current_norm;
+                tracing::debug!(iteration, best_residual, "Best state updated");
+            }
+
+            // Jacobian-freeze feedback
+            if let Some(ref freeze_cfg) = self.jacobian_freezing {
+                if previous_norm > 0.0 && current_norm / previous_norm >= (1.0 - freeze_cfg.threshold) {
+                    if frozen_count > 0 || !force_recompute {
+                        tracing::debug!(iteration, current_norm, previous_norm, "Unfreezing Jacobian");
+                    }
+                    force_recompute = true;
+                    frozen_count = 0;
+                }
+            }
+
+            tracing::debug!(iteration, residual_norm = current_norm, alpha, "Newton iteration complete");
+
+            // Check convergence
+            let converged = if let Some(ref criteria) = self.convergence_criteria {
+                let report = criteria.check(&state, Some(&prev_iteration_state), &residuals, system);
+                if report.is_globally_converged() {
+                    let status = if !system.saturated_variables().is_empty() {
+                        ConvergenceStatus::ControlSaturation
+                    } else {
+                        ConvergenceStatus::Converged
+                    };
+                    tracing::info!(iterations = iteration, final_residual = current_norm, "Converged (criteria)");
+                    return Ok(ConvergedState::with_report(
+                        state, iteration, current_norm, status, report, SimulationMetadata::new(system.input_hash()),
+                    ));
+                }
+                false
+            } else {
+                current_norm < self.tolerance
+            };
+
+            if converged {
+                let status = if !system.saturated_variables().is_empty() {
+                    ConvergenceStatus::ControlSaturation
+                } else {
+                    ConvergenceStatus::Converged
+                };
+                tracing::info!(iterations = iteration, final_residual = current_norm, "Converged");
+                return Ok(ConvergedState::new(
+                    state, iteration, current_norm, status, SimulationMetadata::new(system.input_hash()),
+                ));
+            }
+
+            if let Some(err) = self.check_divergence(current_norm, previous_norm, &mut divergence_count) {
+                tracing::warn!(iteration, residual_norm = current_norm, "Divergence detected");
+                return Err(err);
+            }
+        }
+
+        tracing::warn!(max_iterations = self.max_iterations, final_residual = current_norm, "Did not converge");
+        Err(SolverError::NonConvergence {
+            iterations: self.max_iterations,
+            final_residual: current_norm,
+        })
+    }
+
+    fn with_timeout(mut self, timeout: Duration) -> Self {
+        self.timeout = Some(timeout);
+        self
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::solver::Solver;
+    use crate::system::System;
+    use std::time::Duration;
+
+    #[test]
+    fn test_newton_config_with_timeout() {
+        let cfg = NewtonConfig::default().with_timeout(Duration::from_millis(100));
+        assert_eq!(cfg.timeout, Some(Duration::from_millis(100)));
+    }
+
+    #[test]
+    fn test_newton_config_default() {
+        let cfg = NewtonConfig::default();
+        assert_eq!(cfg.max_iterations, 100);
+        assert!(cfg.tolerance > 0.0 && cfg.tolerance < 1e-3);
+    }
+
+    #[test]
+    fn test_newton_solver_trait_object() {
+        let mut boxed: Box<dyn Solver> = Box::new(NewtonConfig::default());
+        let mut system = System::new();
+        system.finalize().unwrap();
+        assert!(boxed.solve(&mut system).is_err());
+    }
+}

crates/solver/src/strategies/sequential_substitution.rs

@@ -0,0 +1,467 @@
+//! Sequential Substitution (Picard iteration) solver implementation.
+//!
+//! Provides [`PicardConfig`] which implements Picard iteration for solving
+//! systems of non-linear equations. Slower than Newton-Raphson but more robust.
+
+use std::time::{Duration, Instant};
+
+use crate::criteria::ConvergenceCriteria;
+use crate::metadata::SimulationMetadata;
+use crate::solver::{ConvergedState, ConvergenceStatus, Solver, SolverError, TimeoutConfig};
+use crate::system::System;
+
+/// Configuration for the Sequential Substitution (Picard iteration) solver.
+///
+/// Solves x = G(x) by iterating: x_{k+1} = (1-ω)·x_k + ω·G(x_k)
+/// where ω ∈ (0,1] is the relaxation factor.
+#[derive(Debug, Clone, PartialEq)]
+pub struct PicardConfig {
+    /// Maximum iterations. Default: 100.
+    pub max_iterations: usize,
+    /// Convergence tolerance (L2 norm). Default: 1e-6.
+    pub tolerance: f64,
+    /// Relaxation factor ω ∈ (0,1]. Default: 0.5.
+    pub relaxation_factor: f64,
+    /// Optional time budget.
+    pub timeout: Option<Duration>,
+    /// Divergence threshold. Default: 1e10.
+    pub divergence_threshold: f64,
+    /// Consecutive increases before divergence. Default: 5.
+    pub divergence_patience: usize,
+    /// Timeout behavior configuration.
+    pub timeout_config: TimeoutConfig,
+    /// Previous state for ZOH fallback.
+    pub previous_state: Option<Vec<f64>>,
+    /// Residual for previous_state.
+    pub previous_residual: Option<f64>,
+    /// Smart initial state for cold-start.
+    pub initial_state: Option<Vec<f64>>,
+    /// Multi-circuit convergence criteria.
+    pub convergence_criteria: Option<ConvergenceCriteria>,
+}
+
+impl Default for PicardConfig {
+    fn default() -> Self {
+        Self {
+            max_iterations: 100,
+            tolerance: 1e-6,
+            relaxation_factor: 0.5,
+            timeout: None,
+            divergence_threshold: 1e10,
+            divergence_patience: 5,
+            timeout_config: TimeoutConfig::default(),
+            previous_state: None,
+            previous_residual: None,
+            initial_state: None,
+            convergence_criteria: None,
+        }
+    }
+}
+
+impl PicardConfig {
+    /// Sets the initial state for cold-start solving (Story 4.6 — builder pattern).
+    ///
+    /// The solver will start from `state` instead of the zero vector.
+    /// Use [`SmartInitializer::populate_state`] to generate a physically reasonable
+    /// initial guess.
+    pub fn with_initial_state(mut self, state: Vec<f64>) -> Self {
+        self.initial_state = Some(state);
+        self
+    }
+
+    /// Sets multi-circuit convergence criteria (Story 4.7 — builder pattern).
+    ///
+    /// When set, the solver uses [`ConvergenceCriteria::check()`] instead of the
+    /// raw L2-norm `tolerance` check.
+    pub fn with_convergence_criteria(mut self, criteria: ConvergenceCriteria) -> Self {
+        self.convergence_criteria = Some(criteria);
+        self
+    }
+
+    /// Computes the residual norm (L2 norm of the residual vector).
+    fn residual_norm(residuals: &[f64]) -> f64 {
+        residuals.iter().map(|r| r * r).sum::<f64>().sqrt()
+    }
+
+    /// Handles timeout based on configuration (Story 4.5).
+    ///
+    /// Returns either:
+    /// - `Ok(ConvergedState)` with `TimedOutWithBestState` status (default)
+    /// - `Err(SolverError::Timeout)` if `return_best_state_on_timeout` is false
+    /// - Previous state (ZOH) if `zoh_fallback` is true and previous state available
+    fn handle_timeout(
+        &self,
+        best_state: &[f64],
+        best_residual: f64,
+        iterations: usize,
+        timeout: Duration,
+        system: &System,
+    ) -> Result<ConvergedState, SolverError> {
+        // If configured to return error on timeout
+        if !self.timeout_config.return_best_state_on_timeout {
+            return Err(SolverError::Timeout {
+                timeout_ms: timeout.as_millis() as u64,
+            });
+        }
+
+        // If ZOH fallback is enabled and previous state is available
+        if self.timeout_config.zoh_fallback {
+            if let Some(ref prev_state) = self.previous_state {
+                let residual = self.previous_residual.unwrap_or(best_residual);
+                tracing::info!(
+                    iterations = iterations,
+                    residual = residual,
+                    "Returning previous state (ZOH fallback) on timeout"
+                );
+                return Ok(ConvergedState::new(
+                    prev_state.clone(),
+                    iterations,
+                    residual,
+                    ConvergenceStatus::TimedOutWithBestState,
+                    SimulationMetadata::new(system.input_hash()),
+                ));
+            }
+        }
+
+        // Default: return best state encountered during iteration
+        tracing::info!(
+            iterations = iterations,
+            best_residual = best_residual,
+            "Returning best state on timeout"
+        );
+        Ok(ConvergedState::new(
+            best_state.to_vec(),
+            iterations,
+            best_residual,
+            ConvergenceStatus::TimedOutWithBestState,
+            SimulationMetadata::new(system.input_hash()),
+        ))
+    }
+
+    /// Checks for divergence based on residual growth pattern.
+    ///
+    /// Returns `Some(SolverError::Divergence)` if:
+    /// - Residual norm exceeds `divergence_threshold`, or
+    /// - Residual has increased for `divergence_patience`+ consecutive iterations
+    fn check_divergence(
+        &self,
+        current_norm: f64,
+        previous_norm: f64,
+        divergence_count: &mut usize,
+    ) -> Option<SolverError> {
+        // Check absolute threshold
+        if current_norm > self.divergence_threshold {
+            return Some(SolverError::Divergence {
+                reason: format!(
+                    "Residual norm {} exceeds threshold {}",
+                    current_norm, self.divergence_threshold
+                ),
+            });
+        }
+
+        // Check consecutive increases
+        if current_norm > previous_norm {
+            *divergence_count += 1;
+            if *divergence_count >= self.divergence_patience {
+                return Some(SolverError::Divergence {
+                    reason: format!(
+                        "Residual increased for {} consecutive iterations: {:.6e} → {:.6e}",
+                        self.divergence_patience, previous_norm, current_norm
+                    ),
+                });
+            }
+        } else {
+            *divergence_count = 0;
+        }
+
+        None
+    }
+
+    /// Applies relaxation to the state update.
+    ///
+    /// Update formula: x_new = x_old - omega * residual
+    /// where residual = F(x_k) represents the equation residuals.
+    ///
+    /// This is the standard Picard iteration: x_{k+1} = x_k - ω·F(x_k)
+    fn apply_relaxation(state: &mut [f64], residuals: &[f64], omega: f64) {
+        for (x, &r) in state.iter_mut().zip(residuals.iter()) {
+            *x -= omega * r;
+        }
+    }
+}
+
+impl Solver for PicardConfig {
+    fn solve(&mut self, system: &mut System) -> Result<ConvergedState, SolverError> {
+        let start_time = Instant::now();
+
+        tracing::info!(
+            max_iterations = self.max_iterations,
+            tolerance = self.tolerance,
+            relaxation_factor = self.relaxation_factor,
+            divergence_threshold = self.divergence_threshold,
+            divergence_patience = self.divergence_patience,
+            "Sequential Substitution (Picard) solver starting"
+        );
+
+        // Get system dimensions
+        let n_state = system.full_state_vector_len();
+        let n_equations: usize = system
+            .traverse_for_jacobian()
+            .map(|(_, c, _)| c.n_equations())
+            .sum::<usize>()
+            + system.constraints().count()
+            + system.coupling_residual_count();
+
+        // Validate system
+        if n_state == 0 || n_equations == 0 {
+            return Err(SolverError::InvalidSystem {
+                message: "Empty system has no state variables or equations".to_string(),
+            });
+        }
+
+        // Validate state/equation dimensions
+        if n_state != n_equations {
+            return Err(SolverError::InvalidSystem {
+                message: format!(
+                    "State dimension ({}) does not match equation count ({})",
+                    n_state, n_equations
+                ),
+            });
+        }
+
+        // Pre-allocate all buffers (AC: #6 - no heap allocation in iteration loop)
+        // Story 4.6 - AC: #8: Use initial_state if provided, else start from zeros
+        let mut state: Vec<f64> = self
+            .initial_state
+            .as_ref()
+            .map(|s| {
+                debug_assert_eq!(
+                    s.len(),
+                    n_state,
+                    "initial_state length mismatch: expected {}, got {}",
+                    n_state,
+                    s.len()
+                );
+                if s.len() == n_state {
+                    s.clone()
+                } else {
+                    vec![0.0; n_state]
+                }
+            })
+            .unwrap_or_else(|| vec![0.0; n_state]);
+        let mut prev_iteration_state: Vec<f64> = vec![0.0; n_state]; // For convergence delta check
+        let mut residuals: Vec<f64> = vec![0.0; n_equations];
+        let mut divergence_count: usize = 0;
+        let mut previous_norm: f64;
+
+        // Pre-allocate best-state tracking buffer (Story 4.5 - AC: #5)
+        let mut best_state: Vec<f64> = vec![0.0; n_state];
+        let mut best_residual: f64;
+
+        // Initial residual computation
+        system
+            .compute_residuals(&state, &mut residuals)
+            .map_err(|e| SolverError::InvalidSystem {
+                message: format!("Failed to compute initial residuals: {:?}", e),
+            })?;
+
+        let mut current_norm = Self::residual_norm(&residuals);
+
+        // Initialize best state tracking with initial state
+        best_state.copy_from_slice(&state);
+        best_residual = current_norm;
+
+        tracing::debug!(iteration = 0, residual_norm = current_norm, "Initial state");
+
+        // Check if already converged
+        if current_norm < self.tolerance {
+            tracing::info!(
+                iterations = 0,
+                final_residual = current_norm,
+                "System already converged at initial state"
+            );
+            return Ok(ConvergedState::new(
+                state,
+                0,
+                current_norm,
+                ConvergenceStatus::Converged,
+                SimulationMetadata::new(system.input_hash()),
+            ));
+        }
+
+        // Main Picard iteration loop
+        for iteration in 1..=self.max_iterations {
+            // Save state before step for convergence criteria delta checks
+            prev_iteration_state.copy_from_slice(&state);
+
+            // Check timeout at iteration start (Story 4.5 - AC: #1)
+            if let Some(timeout) = self.timeout {
+                if start_time.elapsed() > timeout {
+                    tracing::info!(
+                        iteration = iteration,
+                        elapsed_ms = start_time.elapsed().as_millis(),
+                        timeout_ms = timeout.as_millis(),
+                        best_residual = best_residual,
+                        "Solver timed out"
+                    );
+
+                    // Story 4.5 - AC: #2, #6: Return best state or error based on config
+                    return self.handle_timeout(
+                        &best_state,
+                        best_residual,
+                        iteration - 1,
+                        timeout,
+                        system,
+                    );
+                }
+            }
+
+            // Apply relaxed update: x_new = x_old - omega * residual (AC: #2, #3)
+            Self::apply_relaxation(&mut state, &residuals, self.relaxation_factor);
+
+            // Compute new residuals
+            system
+                .compute_residuals(&state, &mut residuals)
+                .map_err(|e| SolverError::InvalidSystem {
+                    message: format!("Failed to compute residuals: {:?}", e),
+                })?;
+
+            previous_norm = current_norm;
+            current_norm = Self::residual_norm(&residuals);
+
+            // Update best state if residual improved (Story 4.5 - AC: #2)
+            if current_norm < best_residual {
+                best_state.copy_from_slice(&state);
+                best_residual = current_norm;
+                tracing::debug!(
+                    iteration = iteration,
+                    best_residual = best_residual,
+                    "Best state updated"
+                );
+            }
+
+            tracing::debug!(
+                iteration = iteration,
+                residual_norm = current_norm,
+                relaxation_factor = self.relaxation_factor,
+                "Picard iteration complete"
+            );
+
+            // Check convergence (AC: #1, Story 4.7 — criteria-aware)
+            let converged = if let Some(ref criteria) = self.convergence_criteria {
+                let report =
+                    criteria.check(&state, Some(&prev_iteration_state), &residuals, system);
+                if report.is_globally_converged() {
+                    tracing::info!(
+                        iterations = iteration,
+                        final_residual = current_norm,
+                        relaxation_factor = self.relaxation_factor,
+                        "Sequential Substitution converged (criteria)"
+                    );
+                    return Ok(ConvergedState::with_report(
+                        state,
+                        iteration,
+                        current_norm,
+                        ConvergenceStatus::Converged,
+                        report,
+                        SimulationMetadata::new(system.input_hash()),
+                    ));
+                }
+                false
+            } else {
+                current_norm < self.tolerance
+            };
+
+            if converged {
+                tracing::info!(
+                    iterations = iteration,
+                    final_residual = current_norm,
+                    relaxation_factor = self.relaxation_factor,
+                    "Sequential Substitution converged"
+                );
+                return Ok(ConvergedState::new(
+                    state,
+                    iteration,
+                    current_norm,
+                    ConvergenceStatus::Converged,
+                    SimulationMetadata::new(system.input_hash()),
+                ));
+            }
+
+            // Check divergence (AC: #5)
+            if let Some(err) =
+                self.check_divergence(current_norm, previous_norm, &mut divergence_count)
+            {
+                tracing::warn!(
+                    iteration = iteration,
+                    residual_norm = current_norm,
+                    "Divergence detected"
+                );
+                return Err(err);
+            }
+        }
+
+        // Max iterations exceeded
+        tracing::warn!(
+            max_iterations = self.max_iterations,
+            final_residual = current_norm,
+            "Sequential Substitution did not converge"
+        );
+        Err(SolverError::NonConvergence {
+            iterations: self.max_iterations,
+            final_residual: current_norm,
+        })
+    }
+
+    fn with_timeout(mut self, timeout: Duration) -> Self {
+        self.timeout = Some(timeout);
+        self
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::solver::Solver;
+    use crate::system::System;
+    use std::time::Duration;
+
+    #[test]
+    fn test_picard_config_with_timeout() {
+        let timeout = Duration::from_millis(250);
+        let cfg = PicardConfig::default().with_timeout(timeout);
+        assert_eq!(cfg.timeout, Some(timeout));
+    }
+
+    #[test]
+    fn test_picard_config_default_sensible() {
+        let cfg = PicardConfig::default();
+        assert_eq!(cfg.max_iterations, 100);
+        assert!(cfg.tolerance > 0.0 && cfg.tolerance < 1e-3);
+        assert!(cfg.relaxation_factor > 0.0 && cfg.relaxation_factor <= 1.0);
+    }
+
+    #[test]
+    fn test_picard_apply_relaxation_formula() {
+        let mut state = vec![10.0, 20.0];
+        let residuals = vec![1.0, 2.0];
+        PicardConfig::apply_relaxation(&mut state, &residuals, 0.5);
+        assert!((state[0] - 9.5).abs() < 1e-15);
+        assert!((state[1] - 19.0).abs() < 1e-15);
+    }
+
+    #[test]
+    fn test_picard_residual_norm() {
+        let residuals = vec![3.0, 4.0];
+        let norm = PicardConfig::residual_norm(&residuals);
+        assert!((norm - 5.0).abs() < 1e-15);
+    }
+
+    #[test]
+    fn test_picard_solver_trait_object() {
+        let mut boxed: Box<dyn Solver> = Box::new(PicardConfig::default());
+        let mut system = System::new();
+        system.finalize().unwrap();
+        assert!(boxed.solve(&mut system).is_err());
+    }
+}