Add diagram workbench UI with Modelica DoF coaching and ISO glyphs.

Ship the Next.js cycle editor with CAD chrome, technical HX symbols, Fixed/Free boundary guidance, and secondary water/air pressure drop support in the solver stack.

Co-authored-by: Cursor <cursoragent@cursor.com>

crates/solver/src/coupling.rs

@@ -25,10 +25,28 @@ use petgraph::graph::{DiGraph, NodeIndex};
 use serde::{Deserialize, Serialize};
 use std::collections::HashMap;
 
+fn default_duty_scale() -> f64 {
+    1.0
+}
+
 /// Thermal coupling between two circuits via a heat exchanger.
 ///
 /// Heat flows from `hot_circuit` to `cold_circuit` proportional to the
 /// temperature difference and thermal conductance (UA value).
+///
+/// ## Physical (duty-transfer) mode
+///
+/// When [`hot_component`](Self::hot_component) and
+/// [`cold_component`](Self::cold_component) are set, the coupling is **physical**:
+/// the per-coupling state unknown is the transferred heat `Q` [W], closed by the
+/// residual `r = Q − η·duty(hot_component)` where the duty is *measured* from the
+/// solved state via the hot component's `measure_output(Capacity)` (e.g. the real
+/// ε-NTU condenser duty). The cold-side component (a
+/// [`ThermalLoad`](entropyk_components::ThermalLoad)) consumes `Q` in its energy
+/// balance `ṁ·(h_out − h_in) = Q`, so the cold circuit genuinely warms up.
+///
+/// Without the component references the legacy MVP stub applies (the residual
+/// simply pins `Q = 0`), preserved for backward compatibility.
 #[derive(Debug, Clone, Serialize, Deserialize, PartialEq)]
 #[serde(rename_all = "camelCase")]
 pub struct ThermalCoupling {
@@ -42,6 +60,19 @@ pub struct ThermalCoupling {
     pub ua: ThermalConductance,
     /// Efficiency factor (0.0 to 1.0). Default is 1.0 (no losses).
     pub efficiency: f64,
+    /// Multiplier applied to the measured duty before it is injected into the
+    /// receiver. Use `-1.0` when the receiver must be cooled by the source duty
+    /// (for example, chilled water across an evaporator).
+    #[serde(default = "default_duty_scale", alias = "duty_scale")]
+    pub duty_scale: f64,
+    /// Name of the registered component in the hot circuit whose *measured duty*
+    /// (`measure_output(Capacity)`) is transferred (e.g. an ε-NTU `Condenser`).
+    #[serde(default, alias = "hot_component")]
+    pub hot_component: Option<String>,
+    /// Name of the registered `ThermalLoad` component in the cold circuit that
+    /// receives `Q` in its energy balance.
+    #[serde(default, alias = "cold_component")]
+    pub cold_component: Option<String>,
 }
 
 impl ThermalCoupling {
@@ -71,6 +102,9 @@ impl ThermalCoupling {
             cold_circuit,
             ua,
             efficiency: 1.0,
+            duty_scale: 1.0,
+            hot_component: None,
+            cold_component: None,
         }
     }
 
@@ -82,6 +116,35 @@ impl ThermalCoupling {
         self.efficiency = efficiency.clamp(0.0, 1.0);
         self
     }
+
+    /// Sets the signed duty multiplier applied before injecting heat into the
+    /// receiver component.
+    pub fn with_duty_scale(mut self, duty_scale: f64) -> Self {
+        self.duty_scale = if duty_scale.is_finite() {
+            duty_scale
+        } else {
+            1.0
+        };
+        self
+    }
+
+    /// Enables the **physical duty-transfer mode**: the measured duty of
+    /// `hot_component` (via `measure_output(Capacity)`) is transferred, scaled
+    /// by `efficiency`, into `cold_component`'s energy balance (a `ThermalLoad`).
+    pub fn with_interface_components(
+        mut self,
+        hot_component: impl Into<String>,
+        cold_component: impl Into<String>,
+    ) -> Self {
+        self.hot_component = Some(hot_component.into());
+        self.cold_component = Some(cold_component.into());
+        self
+    }
+
+    /// Returns `true` when the coupling is in physical duty-transfer mode.
+    pub fn is_physical(&self) -> bool {
+        self.hot_component.is_some() && self.cold_component.is_some()
+    }
 }
 
 /// Computes heat transfer for a thermal coupling.

crates/solver/src/dof.rs

@@ -0,0 +1,349 @@
+//! System-wide degrees-of-freedom (DoF) bookkeeping.
+//!
+//! A thermodynamic cycle is a square nonlinear system:
+//!
+//! ```text
+//! n_equations == n_unknowns
+//! ```
+//!
+//! Every time a quantity is **fixed** (Dirichlet residual, outlet closure,
+//! inverse constraint), either another quantity must be **freed** (actuator,
+//! emergent pressure, free boundary) or an existing residual must be dropped.
+//!
+//! This module provides:
+//! - re-export of component-level [`EquationRole`];
+//! - unknown labels and a full system ledger;
+//! - hard validation (`SystemDofBalance`) used as a pre-solve gate.
+//!
+//! # Design rules
+//!
+//! 1. Parameters outside the state vector are **not** free unknowns.
+//! 2. Scalar secondary streams (`T_sec`, `C_sec`) are **rating-mode inputs**, not
+//!    system-mode substitutes for a live secondary loop.
+//! 3. Outlet closures (superheat / subcooling / quality / level) consume one DoF
+//!    and must be paired with a free actuator or an intentional residual drop.
+//! 4. Imbalance is an error, not a warning.
+
+use std::fmt;
+
+use thiserror::Error;
+
+pub use entropyk_components::{unspecified_roles, EquationRole};
+
+// ─────────────────────────────────────────────────────────────────────────────
+// Unknown labels
+// ─────────────────────────────────────────────────────────────────────────────
+
+/// Kind of Newton unknown in the full state vector.
+#[derive(Debug, Clone, PartialEq, Eq)]
+pub enum UnknownKind {
+    /// Shared branch mass-flow slot.
+    BranchMassFlow {
+        /// Branch id from topology presolve.
+        branch_id: usize,
+    },
+    /// Edge pressure.
+    EdgePressure {
+        /// Edge ordinal in finalize order.
+        edge_ordinal: usize,
+    },
+    /// Edge enthalpy.
+    EdgeEnthalpy {
+        /// Edge ordinal in finalize order.
+        edge_ordinal: usize,
+    },
+    /// Hard inverse-control variable.
+    InverseControl {
+        /// Bounded-variable id string.
+        id: String,
+    },
+    /// Thermal-coupling heat unknown Q [W].
+    CouplingHeat {
+        /// Coupling index.
+        index: usize,
+    },
+    /// Saturated-PI actuator `u`.
+    SaturatedActuator {
+        /// Controller index.
+        index: usize,
+    },
+    /// Saturated-PI internal state `x`.
+    SaturatedIntegrator {
+        /// Controller index.
+        index: usize,
+    },
+    /// Physical free actuator (EXV opening, fan speed, …).
+    FreeActuator {
+        /// Bounded-variable id string.
+        id: String,
+    },
+    /// Internal component state slot (macro-components / reserved).
+    Internal {
+        /// Slot index within the physical block.
+        index: usize,
+    },
+}
+
+impl fmt::Display for UnknownKind {
+    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
+        match self {
+            Self::BranchMassFlow { branch_id } => write!(f, "m_branch[{branch_id}]"),
+            Self::EdgePressure { edge_ordinal } => write!(f, "P_edge[{edge_ordinal}]"),
+            Self::EdgeEnthalpy { edge_ordinal } => write!(f, "h_edge[{edge_ordinal}]"),
+            Self::InverseControl { id } => write!(f, "ctrl[{id}]"),
+            Self::CouplingHeat { index } => write!(f, "Q_coupling[{index}]"),
+            Self::SaturatedActuator { index } => write!(f, "u_sat[{index}]"),
+            Self::SaturatedIntegrator { index } => write!(f, "x_sat[{index}]"),
+            Self::FreeActuator { id } => write!(f, "actuator[{id}]"),
+            Self::Internal { index } => write!(f, "internal[{index}]"),
+        }
+    }
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// Ledger entries
+// ─────────────────────────────────────────────────────────────────────────────
+
+/// One residual block contributed by a graph node.
+#[derive(Debug, Clone, PartialEq, Eq)]
+pub struct ComponentEquationBlock {
+    /// Human-readable component name (or node index fallback).
+    pub component_name: String,
+    /// Graph node index.
+    pub node_index: usize,
+    /// Declared `n_equations()`.
+    pub n_equations: usize,
+    /// Semantic roles (length should equal `n_equations` when fully migrated).
+    pub roles: Vec<EquationRole>,
+}
+
+/// Outcome of comparing equation count to unknown count.
+#[derive(Debug, Clone, Copy, PartialEq, Eq)]
+pub enum SystemDofBalance {
+    /// Square system — well posed at the dimension level.
+    Balanced,
+    /// More equations than unknowns.
+    OverConstrained {
+        /// `n_equations − n_unknowns`.
+        excess_equations: usize,
+    },
+    /// Fewer equations than unknowns.
+    UnderConstrained {
+        /// `n_unknowns − n_equations`.
+        free_dofs: usize,
+    },
+}
+
+impl SystemDofBalance {
+    /// Builds the balance enum from raw counts.
+    pub fn from_counts(n_equations: usize, n_unknowns: usize) -> Self {
+        match n_equations.cmp(&n_unknowns) {
+            std::cmp::Ordering::Equal => Self::Balanced,
+            std::cmp::Ordering::Greater => Self::OverConstrained {
+                excess_equations: n_equations - n_unknowns,
+            },
+            std::cmp::Ordering::Less => Self::UnderConstrained {
+                free_dofs: n_unknowns - n_equations,
+            },
+        }
+    }
+
+    /// Returns `true` when the system is square.
+    pub fn is_balanced(self) -> bool {
+        matches!(self, Self::Balanced)
+    }
+}
+
+impl fmt::Display for SystemDofBalance {
+    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
+        match self {
+            Self::Balanced => write!(f, "balanced"),
+            Self::OverConstrained { excess_equations } => {
+                write!(f, "over-constrained by {excess_equations}")
+            }
+            Self::UnderConstrained { free_dofs } => {
+                write!(f, "under-constrained by {free_dofs}")
+            }
+        }
+    }
+}
+
+/// Full DoF report for a finalized [`crate::System`].
+#[derive(Debug, Clone, PartialEq)]
+pub struct DofReport {
+    /// Total residual equations assembled by the solver.
+    pub n_equations: usize,
+    /// Total Newton unknowns (`full_state_vector_len`).
+    pub n_unknowns: usize,
+    /// Balance classification.
+    pub balance: SystemDofBalance,
+    /// Per-component residual blocks.
+    pub components: Vec<ComponentEquationBlock>,
+    /// System-level residual roles (constraints, couplings, saturated).
+    pub system_equations: Vec<EquationRole>,
+    /// Catalog of unknowns (may be compact when edge map is large).
+    pub unknowns: Vec<UnknownKind>,
+    /// Human diagnostics (pairing warnings, role-length mismatches, …).
+    pub diagnostics: Vec<String>,
+}
+
+impl DofReport {
+    /// Renders a concise multi-line summary suitable for logs and CLI.
+    pub fn summary(&self) -> String {
+        let mut lines = Vec::new();
+        lines.push(format!(
+            "DoF: equations={} unknowns={} → {}",
+            self.n_equations, self.n_unknowns, self.balance
+        ));
+        for block in &self.components {
+            let roles = if block.roles.is_empty() {
+                "(no roles declared)".to_string()
+            } else {
+                block
+                    .roles
+                    .iter()
+                    .map(|r| r.to_string())
+                    .collect::<Vec<_>>()
+                    .join(", ")
+            };
+            lines.push(format!(
+                "  node {} `{}`: {} eqs — {}",
+                block.node_index, block.component_name, block.n_equations, roles
+            ));
+        }
+        if !self.system_equations.is_empty() {
+            lines.push(format!(
+                "  system-level: {}",
+                self.system_equations
+                    .iter()
+                    .map(|r| r.to_string())
+                    .collect::<Vec<_>>()
+                    .join(", ")
+            ));
+        }
+        for d in &self.diagnostics {
+            lines.push(format!("  ! {d}"));
+        }
+        lines.join("\n")
+    }
+}
+
+/// Errors raised by the hard DoF gate.
+#[derive(Error, Debug, Clone, PartialEq)]
+pub enum SystemDofError {
+    /// Square-system check failed.
+    #[error(
+        "System DoF imbalance: {n_equations} equations vs {n_unknowns} unknowns ({balance}).\n{summary}"
+    )]
+    Imbalance {
+        /// Equation count.
+        n_equations: usize,
+        /// Unknown count.
+        n_unknowns: usize,
+        /// Balance tag.
+        balance: SystemDofBalance,
+        /// Full summary text.
+        summary: String,
+    },
+    /// Component declared roles that disagree with `n_equations()`.
+    #[error(
+        "Component `{component}` equation_roles length {roles_len} != n_equations {n_equations}"
+    )]
+    RoleCountMismatch {
+        /// Component name.
+        component: String,
+        /// Roles length.
+        roles_len: usize,
+        /// Declared equation count.
+        n_equations: usize,
+    },
+}
+
+/// Pads or trims roles to `n`, recording a diagnostic on mismatch.
+pub fn align_roles(
+    component: &str,
+    n_equations: usize,
+    mut roles: Vec<EquationRole>,
+    diagnostics: &mut Vec<String>,
+) -> Vec<EquationRole> {
+    if roles.len() == n_equations {
+        return roles;
+    }
+    if roles.is_empty() {
+        return unspecified_roles(n_equations);
+    }
+    diagnostics.push(format!(
+        "component `{component}`: equation_roles len {} != n_equations {n_equations} (aligned)",
+        roles.len()
+    ));
+    if roles.len() < n_equations {
+        let start = roles.len();
+        roles.extend(unspecified_roles(n_equations - start));
+    } else {
+        roles.truncate(n_equations);
+    }
+    roles
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn balance_from_counts() {
+        assert_eq!(
+            SystemDofBalance::from_counts(9, 9),
+            SystemDofBalance::Balanced
+        );
+        assert_eq!(
+            SystemDofBalance::from_counts(10, 9),
+            SystemDofBalance::OverConstrained {
+                excess_equations: 1
+            }
+        );
+        assert_eq!(
+            SystemDofBalance::from_counts(8, 9),
+            SystemDofBalance::UnderConstrained { free_dofs: 1 }
+        );
+    }
+
+    #[test]
+    fn align_roles_empty_becomes_unspecified() {
+        let mut diag = Vec::new();
+        let roles = align_roles("x", 3, Vec::new(), &mut diag);
+        assert_eq!(roles.len(), 3);
+        assert!(diag.is_empty());
+        assert!(matches!(
+            roles[2],
+            EquationRole::Unspecified { local_index: 2 }
+        ));
+    }
+
+    #[test]
+    fn summary_mentions_balance() {
+        let report = DofReport {
+            n_equations: 2,
+            n_unknowns: 1,
+            balance: SystemDofBalance::OverConstrained {
+                excess_equations: 1,
+            },
+            components: vec![ComponentEquationBlock {
+                component_name: "c".into(),
+                node_index: 0,
+                n_equations: 2,
+                roles: vec![
+                    EquationRole::EnergyBalance {
+                        stream: "refrigerant",
+                    },
+                    EquationRole::OutletClosure { kind: "quality" },
+                ],
+            }],
+            system_equations: vec![],
+            unknowns: vec![],
+            diagnostics: vec!["quality residual without free actuator".into()],
+        };
+        let s = report.summary();
+        assert!(s.contains("over-constrained"));
+        assert!(s.contains("quality"));
+    }
+}

crates/solver/src/error.rs

@@ -54,6 +54,16 @@ pub enum TopologyError {
         /// The circuit ID that was referenced but doesn't exist
         circuit_id: u16,
     },
+
+    /// System equation/unknown count is not square (DoF imbalance).
+    ///
+    /// A real-machine simulation requires `n_equations == n_unknowns`. Fixing a
+    /// quantity without freeing another (or dropping a residual) is rejected here.
+    #[error("System DoF imbalance: {message}")]
+    DofImbalance {
+        /// Human-readable ledger summary.
+        message: String,
+    },
 }
 
 /// Error when adding an edge with port validation.
@@ -95,7 +105,9 @@ pub enum ThermoError {
     },
 
     /// Required fluid backend is not available.
-    #[error("Fluid backend '{backend_name}' is not available. Required version: {required_version}")]
+    #[error(
+        "Fluid backend '{backend_name}' is not available. Required version: {required_version}"
+    )]
     BackendUnavailable {
         /// Name of the missing backend
         backend_name: String,

crates/solver/src/initializer.rs

@@ -10,14 +10,14 @@
 //! 1. Estimate evaporator pressure: `P_evap = P_sat(T_source - ΔT_approach)`
 //! 2. Estimate condenser pressure: `P_cond = P_sat(T_sink + ΔT_approach)`
 //! 3. Clamp `P_evap` to `0.5 * P_critical` if it exceeds the critical pressure
-//! 4. Fill the state vector with `[P, h_default]` per edge, using circuit topology
+//! 4. Fill the state vector with `[ṁ, P, h_default]` per edge, using circuit topology
 //!
 //! # Supported Fluids
 //!
 //! Built-in Antoine coefficients are provided for:
 //! - R134a, R410A, R32, R744 (CO2), R290 (Propane)
 //!
-//! Unknown fluids fall back to sensible defaults (5 bar / 20 bar) with a warning.
+//! Unknown fluids return an explicit error; no pressure guesses are invented.
 //!
 //! # No-Allocation Guarantee
 //!
@@ -29,6 +29,7 @@ use entropyk_core::{Enthalpy, Pressure, Temperature};
 use thiserror::Error;
 
 use crate::system::System;
+use serde::{Deserialize, Serialize};
 
 // ─────────────────────────────────────────────────────────────────────────────
 // Error types
@@ -57,6 +58,13 @@ pub enum InitializerError {
         /// Actual length of the provided slice.
         actual: usize,
     },
+
+    /// No Antoine coefficients are available for the configured fluid.
+    #[error("No Antoine saturation-pressure coefficients are available for {fluid}")]
+    UnsupportedFluid {
+        /// Fluid identifier string.
+        fluid: String,
+    },
 }
 
 // ─────────────────────────────────────────────────────────────────────────────
@@ -197,6 +205,80 @@ pub struct InitializerConfig {
     pub dt_approach: f64,
 }
 
+/// Optional start values for one solver edge or auxiliary unknown group.
+///
+/// These values are numerical guesses only. They must never be interpreted as
+/// imposed boundary conditions or component equations.
+#[derive(Debug, Clone, Default, PartialEq, Serialize, Deserialize)]
+pub struct StartValues {
+    /// Pressure start value [Pa].
+    pub pressure_pa: Option<f64>,
+    /// Enthalpy start value [J/kg].
+    pub enthalpy_j_kg: Option<f64>,
+    /// Mass-flow start value [kg/s].
+    pub mass_flow_kg_s: Option<f64>,
+    /// Temperature start value [K], useful for diagnostics and backend conversion.
+    pub temperature_k: Option<f64>,
+    /// Vapour quality start value [-], when the intended regime is two-phase.
+    pub vapor_quality: Option<f64>,
+}
+
+/// Regime label used to explain why a start value was assigned.
+#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize, Deserialize)]
+#[serde(rename_all = "snake_case")]
+pub enum InitializationRegime {
+    /// High-pressure superheated vapour, typically compressor discharge.
+    HighPressureVapor,
+    /// High-pressure liquid, typically condenser outlet.
+    HighPressureLiquid,
+    /// Low-pressure two-phase mixture, typically EXV outlet.
+    LowPressureTwoPhase,
+    /// Low-pressure superheated vapour, typically compressor suction.
+    LowPressureVapor,
+    /// Secondary water/brine/air branch.
+    Secondary,
+    /// Generic fallback seed.
+    Generic,
+    /// Boundary condition seed from a source/sink component.
+    Boundary,
+    /// Control or actuator unknown seed.
+    Control,
+}
+
+/// One initialization diagnostic entry.
+#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
+pub struct InitializationSeed {
+    /// Human-readable edge/control label.
+    pub label: String,
+    /// Assigned regime.
+    pub regime: InitializationRegime,
+    /// Values written to the state vector.
+    pub values: StartValues,
+}
+
+/// Diagnostics emitted by an initialization pass.
+#[derive(Debug, Clone, Default, PartialEq, Serialize, Deserialize)]
+pub struct InitializationDiagnostics {
+    /// Ordered seed records for edges and control variables.
+    pub seeds: Vec<InitializationSeed>,
+}
+
+impl InitializationDiagnostics {
+    /// Appends a diagnostic seed record.
+    pub fn push(
+        &mut self,
+        label: impl Into<String>,
+        regime: InitializationRegime,
+        values: StartValues,
+    ) {
+        self.seeds.push(InitializationSeed {
+            label: label.into(),
+            regime,
+            values,
+        });
+    }
+}
+
 impl Default for InitializerConfig {
     fn default() -> Self {
         Self {
@@ -248,13 +330,12 @@ impl SmartInitializer {
     /// - `P_evap = P_sat(T_source - ΔT_approach)`, clamped to `0.5 * P_critical`
     /// - `P_cond = P_sat(T_sink + ΔT_approach)`
     ///
-    /// For unknown fluids, returns sensible defaults (5 bar / 20 bar) with a
-    /// `tracing::warn!` log entry.
-    ///
     /// # Errors
     ///
     /// Returns [`InitializerError::TemperatureAboveCritical`] if the adjusted
     /// source temperature exceeds the critical temperature for a known fluid.
+    /// Returns [`InitializerError::UnsupportedFluid`] if no Antoine coefficients
+    /// are available for the configured fluid.
     pub fn estimate_pressures(
         &self,
         t_source: Temperature,
@@ -263,15 +344,7 @@ impl SmartInitializer {
         let fluid_str = self.config.fluid.to_string();
 
         match AntoineCoefficients::for_fluid(&fluid_str) {
-            None => {
-                // Unknown fluid: emit warning and return sensible defaults
-                tracing::warn!(
-                    fluid = %fluid_str,
-                    "Unknown fluid for Antoine estimation — using fallback pressures \
-                     (P_evap = 5 bar, P_cond = 20 bar)"
-                );
-                Ok((Pressure::from_bar(5.0), Pressure::from_bar(20.0)))
-            }
+            None => Err(InitializerError::UnsupportedFluid { fluid: fluid_str }),
             Some(coeffs) => {
                 let t_source_c = t_source.to_celsius();
                 let t_sink_c = t_sink.to_celsius();
@@ -332,9 +405,10 @@ impl SmartInitializer {
     /// Fill a pre-allocated state vector with smart initial guesses.
     ///
     /// No heap allocation is performed. The `state` slice must have length equal
-    /// to `system.state_vector_len()` (i.e., `2 * edge_count`).
+    /// to `system.state_vector_len()` (i.e., `3 * edge_count` for a system of
+    /// refrigerant/hydraulic edges).
     ///
-    /// State layout per edge: `[P_edge_i, h_edge_i]`
+    /// State layout per edge: `[ṁ_edge_i, P_edge_i, h_edge_i]`
     ///
     /// Pressure assignment follows circuit topology:
     /// - Edges in circuit 0 → `p_evap`
@@ -344,7 +418,8 @@ impl SmartInitializer {
     /// # Errors
     ///
     /// Returns [`InitializerError::StateLengthMismatch`] if `state.len()` does
-    /// not match `system.state_vector_len()`.
+    /// not match `system.full_state_vector_len()` (edges plus any inverse-control
+    /// and coupling auxiliary unknowns).
     pub fn populate_state(
         &self,
         system: &System,
@@ -353,7 +428,9 @@ impl SmartInitializer {
         h_default: Enthalpy,
         state: &mut [f64],
     ) -> Result<(), InitializerError> {
-        let expected = system.state_vector_len();
+        // Size against the FULL state vector (the length Newton/Picard expect):
+        // base edge unknowns + inverse-control mappings + coupling residual slots.
+        let expected = system.full_state_vector_len();
         if state.len() != expected {
             return Err(InitializerError::StateLengthMismatch {
                 expected,
@@ -365,11 +442,28 @@ impl SmartInitializer {
         let p_cond_pa = p_cond.to_pascals();
         let h_jkg = h_default.to_joules_per_kg();
 
-        for (i, edge_idx) in system.edge_indices().enumerate() {
+        for edge_idx in system.edge_indices() {
             let circuit = system.edge_circuit(edge_idx);
             let p = if circuit.0 == 0 { p_evap_pa } else { p_cond_pa };
-            state[2 * i] = p;
-            state[2 * i + 1] = h_jkg;
+            let (m_idx, p_idx, h_idx) = system.edge_state_indices_full(edge_idx);
+            // CM1.4: m_idx is BRANCH-shared — multiple edges in the same series
+            // branch point to the same slot. Writing the same seed value multiple
+            // times is idempotent and stays within state bounds (m_idx < state_len).
+            state[m_idx] = crate::system::DEFAULT_MASS_FLOW_SEED_KG_S;
+            state[p_idx] = p;
+            state[h_idx] = h_jkg;
+        }
+
+        // Seed inverse-control unknowns (fan speed, opening, frequency, …) to the
+        // midpoint of their bounds so the cold start sits inside the feasible box
+        // instead of at zero (often an out-of-bounds, non-physical control value).
+        // Coupling auxiliary slots keep their 0.0 default.
+        for (_, idx) in system.control_variable_indices() {
+            if let Some((min, max)) = system.get_bounds_for_state_index(idx) {
+                if min.is_finite() && max.is_finite() && min <= max {
+                    state[idx] = 0.5 * (min + max);
+                }
+            }
         }
 
         Ok(())
@@ -385,6 +479,30 @@ mod tests {
     use super::*;
     use approx::assert_relative_eq;
 
+    #[test]
+    fn test_initialization_diagnostics_records_start_values() {
+        let mut diagnostics = InitializationDiagnostics::default();
+        diagnostics.push(
+            "comp->cond",
+            InitializationRegime::HighPressureVapor,
+            StartValues {
+                pressure_pa: Some(1.2e6),
+                enthalpy_j_kg: Some(430_000.0),
+                mass_flow_kg_s: Some(0.05),
+                temperature_k: None,
+                vapor_quality: None,
+            },
+        );
+
+        assert_eq!(diagnostics.seeds.len(), 1);
+        assert_eq!(diagnostics.seeds[0].label, "comp->cond");
+        assert_eq!(
+            diagnostics.seeds[0].regime,
+            InitializationRegime::HighPressureVapor
+        );
+        assert_eq!(diagnostics.seeds[0].values.pressure_pa, Some(1.2e6));
+    }
+
     // ── Antoine equation unit tests ──────────────────────────────────────────
 
     /// AC: #1, #5 — R134a at 0°C: P_sat ≈ 2.93 bar (293,000 Pa), within 5%
@@ -465,9 +583,9 @@ mod tests {
         assert_relative_eq!(p_cond.to_pascals(), expected_pa, max_relative = 1e-9);
     }
 
-    /// AC: #6 — Unknown fluid returns fallback (5 bar / 20 bar) without panic
+    /// AC: #6 — Unknown fluid returns an explicit error instead of invented pressures.
     #[test]
-    fn test_unknown_fluid_fallback() {
+    fn test_unknown_fluid_returns_error() {
         let init = SmartInitializer::new(InitializerConfig {
             fluid: FluidId::new("R999-Unknown"),
             dt_approach: 5.0,
@@ -476,10 +594,12 @@ mod tests {
             Temperature::from_celsius(5.0),
             Temperature::from_celsius(40.0),
         );
-        assert!(result.is_ok(), "Unknown fluid should not return Err");
-        let (p_evap, p_cond) = result.unwrap();
-        assert_relative_eq!(p_evap.to_bar(), 5.0, max_relative = 1e-9);
-        assert_relative_eq!(p_cond.to_bar(), 20.0, max_relative = 1e-9);
+        assert_eq!(
+            result,
+            Err(InitializerError::UnsupportedFluid {
+                fluid: "R999-Unknown".to_string()
+            })
+        );
     }
 
     /// AC: #1 — Verify evaporator pressure uses T_source - ΔT_approach
@@ -559,12 +679,21 @@ mod tests {
         init.populate_state(&sys, p_evap, p_cond, h_default, &mut state)
             .unwrap();
 
-        // All edges in circuit 0 (single-circuit) → p_evap
-        assert_eq!(state.len(), 4); // 2 edges × 2 entries
-        assert_relative_eq!(state[0], p_evap.to_pascals(), max_relative = 1e-9);
-        assert_relative_eq!(state[1], h_default.to_joules_per_kg(), max_relative = 1e-9);
-        assert_relative_eq!(state[2], p_evap.to_pascals(), max_relative = 1e-9);
-        assert_relative_eq!(state[3], h_default.to_joules_per_kg(), max_relative = 1e-9);
+        // CM1.4: 2-edge linear chain → 1 branch → state_len = 1 + 2×2 = 5
+        // Layout: [0:ṁ_branch, 1:P_e0, 2:h_e0, 3:P_e1, 4:h_e1]
+        assert_eq!(state.len(), 5);
+        // Branch ṁ seeded at DEFAULT_MASS_FLOW_SEED_KG_S
+        assert_relative_eq!(
+            state[0],
+            crate::system::DEFAULT_MASS_FLOW_SEED_KG_S,
+            max_relative = 1e-9
+        );
+        // P and h for edge 0 (circuit 0 → p_evap)
+        assert_relative_eq!(state[1], p_evap.to_pascals(), max_relative = 1e-9);
+        assert_relative_eq!(state[2], h_default.to_joules_per_kg(), max_relative = 1e-9);
+        // P and h for edge 1 (circuit 0 → p_evap)
+        assert_relative_eq!(state[3], p_evap.to_pascals(), max_relative = 1e-9);
+        assert_relative_eq!(state[4], h_default.to_joules_per_kg(), max_relative = 1e-9);
     }
 
     /// AC: #4 — populate_state uses P_cond for circuit 1 edges in multi-circuit system.
@@ -633,13 +762,33 @@ mod tests {
         init.populate_state(&sys, p_evap, p_cond, h_default, &mut state)
             .unwrap();
 
-        assert_eq!(state.len(), 4); // 2 edges × 2 entries
-                                    // Edge 0 (circuit 0) → p_evap
-        assert_relative_eq!(state[0], p_evap.to_pascals(), max_relative = 1e-9);
-        assert_relative_eq!(state[1], h_default.to_joules_per_kg(), max_relative = 1e-9);
-        // Edge 1 (circuit 1) → p_cond
-        assert_relative_eq!(state[2], p_cond.to_pascals(), max_relative = 1e-9);
-        assert_relative_eq!(state[3], h_default.to_joules_per_kg(), max_relative = 1e-9);
+        // CM1.4: 2 isolated 1-edge chains → 2 branches → state_len = 2 + 2×2 = 6
+        // Layout: [0:ṁ_B0, 1:ṁ_B1, 2:P_e0, 3:h_e0, 4:P_e1, 5:h_e1]
+        assert_eq!(state.len(), 6);
+        // Branch ṁ slots seeded at DEFAULT_MASS_FLOW_SEED_KG_S
+        assert_relative_eq!(
+            state[0],
+            crate::system::DEFAULT_MASS_FLOW_SEED_KG_S,
+            max_relative = 1e-9
+        );
+        assert_relative_eq!(
+            state[1],
+            crate::system::DEFAULT_MASS_FLOW_SEED_KG_S,
+            max_relative = 1e-9
+        );
+        // Verify P and h values are seeded correctly for each edge using actual state indices.
+        // Edge 0 is circuit 0 (p_evap), edge 1 is circuit 1 (p_cond).
+        for edge_idx in sys.edge_indices() {
+            let circuit = sys.edge_circuit(edge_idx);
+            let (_m, p, h) = sys.edge_state_indices_full(edge_idx);
+            let expected_p = if circuit.0 == 0 {
+                p_evap.to_pascals()
+            } else {
+                p_cond.to_pascals()
+            };
+            assert_relative_eq!(state[p], expected_p, max_relative = 1e-9);
+            assert_relative_eq!(state[h], h_default.to_joules_per_kg(), max_relative = 1e-9);
+        }
     }
 
     /// AC: #7 — populate_state returns error on length mismatch (no panic).
@@ -689,13 +838,13 @@ mod tests {
         let p_cond = Pressure::from_bar(15.0);
         let h_default = Enthalpy::from_joules_per_kg(400_000.0);
 
-        // Wrong length: system has 2 state entries (1 edge × 2), we provide 5
+        // Wrong length: system has 3 state entries (1 edge × 3), we provide 5
         let mut state = vec![0.0f64; 5];
         let result = init.populate_state(&sys, p_evap, p_cond, h_default, &mut state);
         assert!(matches!(
             result,
             Err(InitializerError::StateLengthMismatch {
-                expected: 2,
+                expected: 3,
                 actual: 5
             })
         ));

crates/solver/src/inverse/calibration.rs

@@ -2,12 +2,12 @@
 //!
 //! This module provides a higher-level orchestration layer on top of the existing
 //! one-shot calibration infrastructure (Story 5.5). It automates the process of
-//! estimating Calib parameters (f_m, f_ua, f_power, etc.) from measured data.
+//! estimating Calib Z-factors (z_flow, z_ua, z_power, etc.) from measured data.
 //!
 //! # Modes
 //!
 //! - **Sequential** (default): calibrates one factor at a time in the recommended order
-//!   f_m → f_dp → f_ua → f_power → f_etav. More stable for nonlinear interactions.
+//!   z_flow → z_flow_eco → z_dp → z_ua → z_power → z_etav. More stable for nonlinear interactions.
 //! - **Simultaneous**: swaps all factors at once for a single One-Shot solve. Faster
 //!   but less robust.
 //!
@@ -21,12 +21,12 @@
 //! let problem = CalibrationProblem::new()
 //!     .with_mode(CalibrationMode::Sequential)
 //!     .add_request(CalibRequest::new(
-//!         CalibFactor::FUa, "evaporator", (0.1, 10.0), 1.0,
+//!         CalibFactor::ZUa, "evaporator", (0.1, 10.0), 1.0,
 //!     ))
 //!     .add_target(CalibrationTarget::capacity("evaporator", 4015.0));
 //!
 //! let result = problem.calibrate(&mut sys, &mut solver)?;
-//! println!("f_ua = {}", result.estimated_factor("evaporator.f_ua"));
+//! println!("f_ua = {}", result.estimated_factor("evaporator.z_ua"));
 //! ```
 
 use std::collections::{HashMap, HashSet};
@@ -36,8 +36,7 @@ use serde::{Deserialize, Serialize};
 use thiserror::Error;
 
 use super::{
-    BoundedVariable, BoundedVariableId, ComponentOutput, Constraint, ConstraintId,
-    ConstraintError,
+    BoundedVariable, BoundedVariableId, ComponentOutput, Constraint, ConstraintError, ConstraintId,
 };
 use crate::solver::{Solver, SolverError};
 use crate::strategies::NewtonConfig;
@@ -52,49 +51,54 @@ use crate::system::System;
 /// Each variant maps to a field in [`entropyk_core::Calib`].
 #[derive(Debug, Clone, Copy, PartialEq, Eq, Hash, Serialize, Deserialize)]
 pub enum CalibFactor {
-    /// f_m: mass flow multiplier (Compressor, Expansion Valve)
-    FM,
-    /// f_dp: pressure drop multiplier (Pipe, Heat Exchanger)
-    FDp,
-    /// f_ua: UA multiplier (Evaporator, Condenser)
-    FUa,
-    /// f_power: power multiplier (Compressor)
-    FPower,
-    /// f_etav: volumetric efficiency multiplier (Compressor)
-    FEtav,
+    /// z_flow: mass flow multiplier (BOLT Z_flow_suc; Compressor, Expansion Valve)
+    ZFlow,
+    /// z_flow_eco: economizer injection mass flow multiplier (BOLT Z_flow_eco)
+    ZFlowEco,
+    /// z_dp: pressure drop multiplier (BOLT Z_dpc; Pipe, Heat Exchanger)
+    ZDp,
+    /// z_ua: UA multiplier (BOLT Z_UA; Evaporator, Condenser)
+    ZUa,
+    /// z_power: power multiplier (BOLT Z_power; Compressor)
+    ZPower,
+    /// z_etav: volumetric efficiency multiplier (Compressor)
+    ZEtav,
 }
 
 impl CalibFactor {
-    /// Returns the canonical short name (e.g. "f_m").
+    /// Returns the canonical Z-factor name (e.g. `"z_flow"`).
     pub fn as_str(&self) -> &'static str {
         match self {
-            CalibFactor::FM => "f_m",
-            CalibFactor::FDp => "f_dp",
-            CalibFactor::FUa => "f_ua",
-            CalibFactor::FPower => "f_power",
-            CalibFactor::FEtav => "f_etav",
+            CalibFactor::ZFlow => entropyk_core::Z_FLOW,
+            CalibFactor::ZFlowEco => entropyk_core::Z_FLOW_ECO,
+            CalibFactor::ZDp => entropyk_core::Z_DP,
+            CalibFactor::ZUa => entropyk_core::Z_UA,
+            CalibFactor::ZPower => entropyk_core::Z_POWER,
+            CalibFactor::ZEtav => entropyk_core::Z_ETAV,
         }
     }
 
-    /// Recommended calibration order: f_m → f_dp → f_ua → f_power → f_etav.
+    /// Recommended calibration order: z_flow → z_flow_eco → z_dp → z_ua → z_power → z_etav.
     pub fn calibration_order() -> &'static [CalibFactor] {
         &[
-            CalibFactor::FM,
-            CalibFactor::FDp,
-            CalibFactor::FUa,
-            CalibFactor::FPower,
-            CalibFactor::FEtav,
+            CalibFactor::ZFlow,
+            CalibFactor::ZFlowEco,
+            CalibFactor::ZDp,
+            CalibFactor::ZUa,
+            CalibFactor::ZPower,
+            CalibFactor::ZEtav,
         ]
     }
 
     /// Returns the default bounds for this factor type.
     pub fn default_bounds(&self) -> (f64, f64) {
         match self {
-            CalibFactor::FM => (0.5, 2.0),
-            CalibFactor::FDp => (0.5, 2.0),
-            CalibFactor::FUa => (0.1, 10.0),
-            CalibFactor::FPower => (0.5, 2.0),
-            CalibFactor::FEtav => (0.5, 2.0),
+            CalibFactor::ZFlow => (0.5, 2.0),
+            CalibFactor::ZFlowEco => (0.5, 2.0),
+            CalibFactor::ZDp => (0.5, 2.0),
+            CalibFactor::ZUa => (0.1, 10.0),
+            CalibFactor::ZPower => (0.5, 2.0),
+            CalibFactor::ZEtav => (0.5, 2.0),
         }
     }
 }
@@ -352,10 +356,7 @@ pub enum CalibrationError {
         "DoF mismatch: {n_targets} targets for {n_requests} calibration requests \
          (must be equal)"
     )]
-    DoFMismatch {
-        n_targets: usize,
-        n_requests: usize,
-    },
+    DoFMismatch { n_targets: usize, n_requests: usize },
 
     /// A referenced component does not exist in the system.
     #[error("Component '{component_id}' not registered in the system")]
@@ -522,7 +523,10 @@ impl CalibrationProblem {
         let mut pairs: Vec<(&CalibRequest, &CalibrationTarget)> =
             self.requests.iter().zip(self.targets.iter()).collect();
         pairs.sort_by_key(|(req, _)| {
-            order.iter().position(|f| *f == req.factor).unwrap_or(usize::MAX)
+            order
+                .iter()
+                .position(|f| *f == req.factor)
+                .unwrap_or(usize::MAX)
         });
 
         for (req, target) in pairs {
@@ -547,14 +551,12 @@ impl CalibrationProblem {
                     reason: format!("Bounded variable error for {}: {e}", req.key()),
                 })
             })?;
-            system
-                .add_bounded_variable(bv)
-                .map_err(|e| {
-                    let _ = system.remove_constraint(&constraint_id);
-                    CalibrationError::ConstraintError(ConstraintError::InvalidConfiguration {
-                        reason: format!("Failed to add bounded variable: {e}"),
-                    })
-                })?;
+            system.add_bounded_variable(bv).map_err(|e| {
+                let _ = system.remove_constraint(&constraint_id);
+                CalibrationError::ConstraintError(ConstraintError::InvalidConfiguration {
+                    reason: format!("Failed to add bounded variable: {e}"),
+                })
+            })?;
 
             system
                 .link_constraint_to_control(&constraint_id, &var_id)
@@ -567,7 +569,9 @@ impl CalibrationProblem {
                 })?;
 
             // Re-finalize to update state vector layout
-            system.finalize().map_err(|_| CalibrationError::SystemNotFinalized)?;
+            system
+                .finalize()
+                .map_err(|_| CalibrationError::SystemNotFinalized)?;
 
             // Create solver with correct initial state
             let initial_state = vec![0.0; system.full_state_vector_len()];
@@ -621,10 +625,8 @@ impl CalibrationProblem {
                         &converged.state[..base_len],
                         &control_values,
                     );
-                    let computed_output = computed_outputs
-                        .get(&constraint_id)
-                        .copied()
-                        .unwrap_or(0.0);
+                    let computed_output =
+                        computed_outputs.get(&constraint_id).copied().unwrap_or(0.0);
                     let residual = target.measured_value - computed_output;
 
                     // P-6: Populate residuals HashMap
@@ -645,7 +647,10 @@ impl CalibrationProblem {
                         result.saturated_factors.push(req.key());
                     }
                 }
-                Err(SolverError::NonConvergence { iterations, final_residual }) => {
+                Err(SolverError::NonConvergence {
+                    iterations,
+                    final_residual,
+                }) => {
                     let _ = system.remove_constraint(&constraint_id);
                     let _ = system.remove_bounded_variable(&var_id);
                     let _ = system.finalize();
@@ -728,13 +733,11 @@ impl CalibrationProblem {
                     reason: format!("Bounded variable error for {}: {e}", req.key()),
                 })
             })?;
-            system
-                .add_bounded_variable(bv)
-                .map_err(|e| {
-                    CalibrationError::ConstraintError(ConstraintError::InvalidConfiguration {
-                        reason: format!("Failed to add bounded variable: {e}"),
-                    })
-                })?;
+            system.add_bounded_variable(bv).map_err(|e| {
+                CalibrationError::ConstraintError(ConstraintError::InvalidConfiguration {
+                    reason: format!("Failed to add bounded variable: {e}"),
+                })
+            })?;
 
             system
                 .link_constraint_to_control(&constraint_id, &var_id)
@@ -802,15 +805,16 @@ impl CalibrationProblem {
                 );
                 for (i, req) in self.requests.iter().enumerate() {
                     let constraint_id = ConstraintId::new(format!("calib_{}", req.key()));
-                    let computed_output = computed_outputs
-                        .get(&constraint_id)
-                        .copied()
-                        .unwrap_or(0.0);
+                    let computed_output =
+                        computed_outputs.get(&constraint_id).copied().unwrap_or(0.0);
                     let residual = self.targets[i].measured_value - computed_output;
                     result.residuals.insert(req.key(), residual);
                 }
             }
-            Err(SolverError::NonConvergence { iterations, final_residual }) => {
+            Err(SolverError::NonConvergence {
+                iterations,
+                final_residual,
+            }) => {
                 for cid in &constraint_ids {
                     system.remove_constraint(cid);
                 }
@@ -899,40 +903,46 @@ mod tests {
     #[test]
     fn test_calib_factor_order() {
         let order = CalibFactor::calibration_order();
-        assert_eq!(order.len(), 5);
-        assert_eq!(order[0], CalibFactor::FM);
-        assert_eq!(order[1], CalibFactor::FDp);
-        assert_eq!(order[2], CalibFactor::FUa);
-        assert_eq!(order[3], CalibFactor::FPower);
-        assert_eq!(order[4], CalibFactor::FEtav);
+        assert_eq!(order.len(), 6);
+        assert_eq!(order[0], CalibFactor::ZFlow);
+        assert_eq!(order[1], CalibFactor::ZFlowEco);
+        assert_eq!(order[2], CalibFactor::ZDp);
+        assert_eq!(order[3], CalibFactor::ZUa);
+        assert_eq!(order[4], CalibFactor::ZPower);
+        assert_eq!(order[5], CalibFactor::ZEtav);
     }
 
     #[test]
     fn test_calib_factor_default_bounds() {
-        let (lo, hi) = CalibFactor::FUa.default_bounds();
+        let (lo, hi) = CalibFactor::ZUa.default_bounds();
         assert_eq!(lo, 0.1);
         assert_eq!(hi, 10.0);
 
-        let (lo, hi) = CalibFactor::FM.default_bounds();
+        let (lo, hi) = CalibFactor::ZFlow.default_bounds();
+        assert_eq!(lo, 0.5);
+        assert_eq!(hi, 2.0);
+
+        let (lo, hi) = CalibFactor::ZFlowEco.default_bounds();
         assert_eq!(lo, 0.5);
         assert_eq!(hi, 2.0);
     }
 
     #[test]
     fn test_calib_factor_display() {
-        assert_eq!(format!("{}", CalibFactor::FM), "f_m");
-        assert_eq!(format!("{}", CalibFactor::FUa), "f_ua");
+        assert_eq!(format!("{}", CalibFactor::ZFlow), "z_flow");
+        assert_eq!(format!("{}", CalibFactor::ZFlowEco), "z_flow_eco");
+        assert_eq!(format!("{}", CalibFactor::ZUa), "z_ua");
     }
 
     #[test]
     fn test_calib_request_key() {
-        let req = CalibRequest::new(CalibFactor::FUa, "evaporator", (0.1, 10.0), 1.0);
-        assert_eq!(req.key(), "evaporator.f_ua");
+        let req = CalibRequest::new(CalibFactor::ZUa, "evaporator", (0.1, 10.0), 1.0);
+        assert_eq!(req.key(), "evaporator.z_ua");
     }
 
     #[test]
     fn test_calib_request_default_bounds() {
-        let req = CalibRequest::with_default_bounds(CalibFactor::FUa, "evaporator", 1.0);
+        let req = CalibRequest::with_default_bounds(CalibFactor::ZUa, "evaporator", 1.0);
         assert_eq!(req.bounds, (0.1, 10.0));
     }
 
@@ -980,7 +990,7 @@ mod tests {
         let p = CalibrationProblem::new()
             .with_mode(CalibrationMode::Simultaneous)
             .add_request(CalibRequest::new(
-                CalibFactor::FUa,
+                CalibFactor::ZUa,
                 "evaporator",
                 (0.1, 10.0),
                 1.0,
@@ -994,15 +1004,14 @@ mod tests {
 
     #[test]
     fn test_calibration_error_dof_mismatch() {
-        let p = CalibrationProblem::new()
-            .add_request(CalibRequest::new(
-                CalibFactor::FUa,
-                "evaporator",
-                (0.1, 10.0),
-                1.0,
-            ));
+        let p = CalibrationProblem::new().add_request(CalibRequest::new(
+            CalibFactor::ZUa,
+            "evaporator",
+            (0.1, 10.0),
+            1.0,
+        ));
 
-        let mut sys = System::new();
+        let sys = System::new();
         let err = p.validate(&sys).unwrap_err();
         assert!(matches!(err, CalibrationError::DoFMismatch { .. }));
     }
@@ -1011,14 +1020,14 @@ mod tests {
     fn test_calibration_error_component_not_found() {
         let p = CalibrationProblem::new()
             .add_request(CalibRequest::new(
-                CalibFactor::FUa,
+                CalibFactor::ZUa,
                 "nonexistent",
                 (0.1, 10.0),
                 1.0,
             ))
             .add_target(CalibrationTarget::capacity("nonexistent", 4015.0));
 
-        let mut sys = System::new();
+        let sys = System::new();
         let err = p.validate(&sys).unwrap_err();
         assert!(matches!(err, CalibrationError::ComponentNotFound { .. }));
     }
@@ -1028,16 +1037,18 @@ mod tests {
         let mut result = CalibrationResult::new();
         result
             .estimated_factors
-            .insert("evaporator.f_ua".to_string(), 1.15);
+            .insert("evaporator.z_ua".to_string(), 1.15);
         result
             .estimated_factors
-            .insert("compressor.f_m".to_string(), 0.95);
-        result.residuals.insert("evaporator.f_ua".to_string(), 0.02);
+            .insert("compressor.z_flow".to_string(), 0.95);
+        result.residuals.insert("evaporator.z_ua".to_string(), 0.02);
         result.mape = 1.5;
         result.max_abs_error = 0.05;
         result.iterations = 42;
         result.converged = true;
-        result.saturated_factors.push("compressor.f_m".to_string());
+        result
+            .saturated_factors
+            .push("compressor.z_flow".to_string());
 
         let json = serde_json::to_string(&result).unwrap();
         let result2: CalibrationResult = serde_json::from_str(&json).unwrap();
@@ -1046,7 +1057,7 @@ mod tests {
 
     #[test]
     fn test_calib_factor_serialize_roundtrip() {
-        let factor = CalibFactor::FUa;
+        let factor = CalibFactor::ZUa;
         let json = serde_json::to_string(&factor).unwrap();
         let factor2: CalibFactor = serde_json::from_str(&json).unwrap();
         assert_eq!(factor, factor2);

crates/solver/src/inverse/constraint.rs

@@ -164,17 +164,23 @@ impl ComponentOutput {
 
     /// Creates a Superheat output for the given component.
     pub fn superheat_for(component_id: &str) -> Self {
-        ComponentOutput::Superheat { component_id: component_id.to_string() }
+        ComponentOutput::Superheat {
+            component_id: component_id.to_string(),
+        }
     }
 
     /// Creates a Subcooling output for the given component.
     pub fn subcooling_for(component_id: &str) -> Self {
-        ComponentOutput::Subcooling { component_id: component_id.to_string() }
+        ComponentOutput::Subcooling {
+            component_id: component_id.to_string(),
+        }
     }
 
     /// Creates a Capacity output for the given component.
     pub fn capacity_for(component_id: &str) -> Self {
-        ComponentOutput::Capacity { component_id: component_id.to_string() }
+        ComponentOutput::Capacity {
+            component_id: component_id.to_string(),
+        }
     }
 }
 

crates/solver/src/inverse/mod.rs

@@ -45,6 +45,8 @@ pub mod bounded;
 pub mod calibration;
 pub mod constraint;
 pub mod embedding;
+pub mod override_network;
+pub mod saturated_control;
 
 pub use bounded::{
     clip_step, BoundedVariable, BoundedVariableError, BoundedVariableId, SaturationInfo,
@@ -56,3 +58,5 @@ pub use calibration::{
 };
 pub use constraint::{ComponentOutput, Constraint, ConstraintError, ConstraintId};
 pub use embedding::{ControlMapping, DoFError, InverseControlConfig};
+pub use override_network::{eval_error_signal, eval_error_weights, Combine, Objective};
+pub use saturated_control::{SaturatedControlError, SaturatedController, Saturation};

crates/solver/src/inverse/override_network.rs

@@ -0,0 +1,242 @@
+//! Steady-state **override / selector control** network.
+//!
+//! Real supervisory controllers drive a *single* actuator from *several*
+//! competing objectives: a primary setpoint (e.g. capacity, superheat) plus a
+//! set of operating-envelope protections (SST low, SDT high, DGT high,
+//! min/max frequency, …). Only one objective is "in authority" at a time; the
+//! others act as overrides that take over when a limit is about to be crossed.
+//!
+//! This mirrors the `BOLT.Control.SteadyState.SetpointControl` library used in
+//! the reference Modelica chillers (61WH / 61AQ / NG-Screw), where the pattern
+//! is `ErrorCalculation` blocks feeding a tree of `Min` / `Max` selectors into a
+//! single `SetpointController`. See also the ALES/UTC report *Supervisory
+//! Control Formulation: Centrifugal System* (Mancuso & Morari, 2016).
+//!
+//! # Formulation
+//!
+//! Each objective `i` computes a **normalized** error
+//!
+//! ```text
+//!   e_i = gain_i · (setpoint_i − measurement_i)
+//! ```
+//!
+//! The `gain_i` normalizes every objective to a comparable scale (e.g.
+//! `1/(freq_max − freq_min)`, `−1/(T_dgt_max − T_dgt_min)`), so that the
+//! selector compares apples to apples — this is the "same-gain" principle from
+//! the reference: after normalization a *single* unit controller integrates the
+//! selected error.
+//!
+//! Errors are folded left-to-right into a single selected error `E`:
+//!
+//! ```text
+//!   acc_0 = e_0
+//!   acc_i = combine_i(acc_{i-1}, e_i)     with combine_i ∈ {Min, Max}
+//!   E     = acc_{n-1}
+//! ```
+//!
+//! The fold order encodes **priority**: place higher-priority protections later
+//! in the chain (this reproduces the linear `min/max/min/…` selector chains of
+//! `CompressorControl` / `EXVControl`).
+//!
+//! # Smoothing (convergence)
+//!
+//! `Min` / `Max` are replaced by the C^∞ `softMin` / `softMax`
+//! (`entropyk_core::smoothing`) with sharpness `alpha`. Using a smooth selector
+//! with an **exact analytic Jacobian** (rather than a non-smooth `min`/`max`
+//! with a semismooth Newton step) is the "Jacobian-smoothing" approach that the
+//! nonlinear-complementarity literature reports as markedly more robust and
+//! faster to converge (fewer Newton iterations, no chattering at the selector
+//! kinks). `alpha` can be annealed toward zero by an outer continuation loop for
+//! a sharp final solution.
+
+use entropyk_core::smoothing::{smooth_max, smooth_min};
+
+use super::constraint::ComponentOutput;
+
+/// How an objective combines with the running selected error.
+#[derive(Debug, Clone, Copy, PartialEq, Eq)]
+pub enum Combine {
+    /// Take the (smooth) minimum of the accumulator and this objective's error.
+    Min,
+    /// Take the (smooth) maximum of the accumulator and this objective's error.
+    Max,
+}
+
+/// A single control objective feeding an override network.
+///
+/// The normalized error is `gain · (setpoint − measurement)`, where
+/// `measurement` is the current value of [`Objective::output`].
+#[derive(Debug, Clone)]
+pub struct Objective {
+    /// The measured plant output for this objective.
+    pub output: ComponentOutput,
+    /// Target value for the measured output (SI units).
+    pub setpoint: f64,
+    /// Normalization/sign gain for this objective's error.
+    pub gain: f64,
+    /// Selector applied between the running accumulator and this objective.
+    /// Ignored for the first objective (which seeds the accumulator).
+    pub combine: Combine,
+}
+
+impl Objective {
+    /// Builds an objective with the given output, setpoint, gain and combinator.
+    pub fn new(output: ComponentOutput, setpoint: f64, gain: f64, combine: Combine) -> Self {
+        Self {
+            output,
+            setpoint,
+            gain,
+            combine,
+        }
+    }
+
+    /// Normalized error `e = gain · (setpoint − measurement)`.
+    #[inline]
+    pub fn error(&self, measurement: f64) -> f64 {
+        self.gain * (self.setpoint - measurement)
+    }
+}
+
+/// `softMin` value and partials `(value, ∂/∂a, ∂/∂b)`.
+#[inline]
+fn soft_min_partials(a: f64, b: f64, k: f64) -> (f64, f64, f64) {
+    let d = ((a - b) * (a - b) + k * k).sqrt();
+    let s = if d > 0.0 { (a - b) / d } else { 0.0 };
+    (smooth_min(a, b, k), 0.5 * (1.0 - s), 0.5 * (1.0 + s))
+}
+
+/// `softMax` value and partials `(value, ∂/∂a, ∂/∂b)`.
+#[inline]
+fn soft_max_partials(a: f64, b: f64, k: f64) -> (f64, f64, f64) {
+    let d = ((a - b) * (a - b) + k * k).sqrt();
+    let s = if d > 0.0 { (a - b) / d } else { 0.0 };
+    (smooth_max(a, b, k), 0.5 * (1.0 + s), 0.5 * (1.0 - s))
+}
+
+/// Evaluates the selected error `E` for the given objectives and their measured
+/// values (`measured[i]` corresponds to `objectives[i]`).
+///
+/// Panics in debug builds if the slice lengths differ. Returns `0.0` for an
+/// empty objective list.
+pub fn eval_error_signal(objectives: &[Objective], measured: &[f64], alpha: f64) -> f64 {
+    debug_assert_eq!(objectives.len(), measured.len());
+    if objectives.is_empty() {
+        return 0.0;
+    }
+    let mut acc = objectives[0].error(measured[0]);
+    for i in 1..objectives.len() {
+        let e = objectives[i].error(measured[i]);
+        acc = match objectives[i].combine {
+            Combine::Min => smooth_min(acc, e, alpha),
+            Combine::Max => smooth_max(acc, e, alpha),
+        };
+    }
+    acc
+}
+
+/// Computes the selector weights `w_i = ∂E/∂e_i` for each objective via a
+/// forward/backward sweep over the fold. These let the caller assemble the
+/// exact plant-coupling Jacobian: `∂E/∂measurement_i = w_i · (−gain_i)`.
+pub fn eval_error_weights(objectives: &[Objective], measured: &[f64], alpha: f64) -> Vec<f64> {
+    debug_assert_eq!(objectives.len(), measured.len());
+    let n = objectives.len();
+    let mut weights = vec![0.0; n];
+    if n == 0 {
+        return weights;
+    }
+    if n == 1 {
+        weights[0] = 1.0;
+        return weights;
+    }
+
+    // Forward: accumulate value and store per-step partials.
+    let mut pa = vec![0.0; n]; // ∂acc_i/∂acc_{i-1}
+    let mut pb = vec![0.0; n]; // ∂acc_i/∂e_i
+    let mut acc = objectives[0].error(measured[0]);
+    for i in 1..n {
+        let e = objectives[i].error(measured[i]);
+        let (val, da, db) = match objectives[i].combine {
+            Combine::Min => soft_min_partials(acc, e, alpha),
+            Combine::Max => soft_max_partials(acc, e, alpha),
+        };
+        pa[i] = da;
+        pb[i] = db;
+        acc = val;
+    }
+
+    // Backward: propagate ∂E/∂acc back to each e_i.
+    let mut g = 1.0;
+    for i in (1..n).rev() {
+        weights[i] = g * pb[i];
+        g *= pa[i];
+    }
+    weights[0] = g;
+    weights
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    fn obj(setpoint: f64, gain: f64, combine: Combine) -> Objective {
+        Objective::new(
+            ComponentOutput::Temperature {
+                component_id: "c".to_string(),
+            },
+            setpoint,
+            gain,
+            combine,
+        )
+    }
+
+    #[test]
+    fn single_objective_is_plain_error() {
+        let objs = vec![obj(5.0, -0.5, Combine::Min)];
+        let e = eval_error_signal(&objs, &[7.0], 1e-3);
+        assert!((e - (-0.5 * (5.0 - 7.0))).abs() < 1e-12);
+        let w = eval_error_weights(&objs, &[7.0], 1e-3);
+        assert_eq!(w, vec![1.0]);
+    }
+
+    #[test]
+    fn min_selects_smaller_error_and_routes_weight() {
+        // Two objectives; e_0 large, e_1 small → Min picks ~e_1, so weight ~1 on
+        // objective 1 and ~0 on objective 0.
+        let objs = vec![obj(10.0, 1.0, Combine::Min), obj(0.0, 1.0, Combine::Min)];
+        // measured: obj0 at 5 → e0 = 5; obj1 at 5 → e1 = -5. min → ~-5.
+        let e = eval_error_signal(&objs, &[5.0, 5.0], 1e-4);
+        assert!((e - (-5.0)).abs() < 1e-2, "E={e}");
+        let w = eval_error_weights(&objs, &[5.0, 5.0], 1e-4);
+        assert!(w[1] > 0.98 && w[0] < 0.02, "weights={w:?}");
+        // Weights of a smooth selector sum to 1 (convex combination).
+        assert!((w[0] + w[1] - 1.0).abs() < 1e-9);
+    }
+
+    #[test]
+    fn weights_match_finite_difference() {
+        let objs = vec![
+            obj(8.0, 0.7, Combine::Min),
+            obj(2.0, -1.3, Combine::Max),
+            obj(-1.0, 0.9, Combine::Min),
+        ];
+        let measured = [6.0, 3.0, 0.5];
+        let alpha = 0.05;
+        let w = eval_error_weights(&objs, &measured, alpha);
+        let h = 1e-6;
+        for i in 0..objs.len() {
+            // dE/de_i via FD on the measurement, then convert: dE/dm_i = -gain_i·w_i.
+            let mut mp = measured;
+            let mut mm = measured;
+            mp[i] += h;
+            mm[i] -= h;
+            let de_dm = (eval_error_signal(&objs, &mp, alpha)
+                - eval_error_signal(&objs, &mm, alpha))
+                / (2.0 * h);
+            let expected = -objs[i].gain * w[i];
+            assert!(
+                (de_dm - expected).abs() < 1e-4,
+                "objective {i}: FD {de_dm} vs analytic {expected}"
+            );
+        }
+    }
+}

crates/solver/src/inverse/saturated_control.rs

@@ -0,0 +1,631 @@
+//! Steady-state saturated PI supervisory control, expressed as algebraic
+//! residual pairs solved jointly with the plant model.
+//!
+//! # Motivation
+//!
+//! Real HVAC machines are run by supervisory controllers: a compressor tracks a
+//! leaving-water-temperature (LWT) setpoint, an expansion valve tracks a suction
+//! saturated temperature (SST) / superheat setpoint, a hot-gas-bypass valve
+//! keeps the compressor off its surge line, etc. To *qualify a machine with its
+//! controls* at steady state, each control loop must be represented as equations
+//! that are solved together with the thermodynamic residuals — not applied as an
+//! outer optimisation loop.
+//!
+//! A naive embedding (`measure(x) − setpoint = 0` plus a hard-clipped actuator)
+//! suffers from **integrator wind-up**: when the actuator saturates, the
+//! setpoint can no longer be met, yet the hard constraint still demands it,
+//! leaving the system either infeasible or stuck against a bound with no smooth
+//! release.
+//!
+//! # Formulation
+//!
+//! This module implements the steady-state saturated-PI formulation: for a loop
+//! with actuator `u ∈ [u_min, u_max]`, controlled output `y`, and setpoint
+//! `y_ref`, we introduce **one internal variable `x`** and **two residuals**:
+//!
+//! ```text
+//!   r_u = u − ( Z · S(x) + Y )                         (actuator law)
+//!   r_y = K · ( y_ref − y ) − ( x − S(x) )             (control law)
+//! ```
+//!
+//! with
+//!
+//! ```text
+//!   S(x) = ( |x + Q| − |x − Q| ) / 2   = clamp(x, −Q, Q)   (saturation)
+//!   Z    = ( u_max − u_min ) / 2
+//!   Y    = ( u_max + u_min ) / 2
+//! ```
+//!
+//! * `K` carries the sign of the proportional gain (`+1` when raising `u` raises
+//!   `y`, `−1` otherwise). With the offset-free control law below its **magnitude
+//!   no longer changes the tracked value** (only convergence scaling), so tuning
+//!   reduces to picking the correct sign.
+//! * `Q > 0` sets the width of the linear (unsaturated) band of the internal
+//!   variable.
+//!
+//! # Offset-free (integral-equivalent) formulation
+//!
+//! This is the canonical steady-state saturated-PI formulation of Mancuso &
+//! Morari (*Supervisory Control Formulation: Centrifugal System*, UTC, 2016,
+//! eq. A.1): `K·e + S(x) − x = 0` with `e = y_ref − y`, i.e.
+//! `r_y = K·(y_ref − y) − (x − S(x))`. The key term is **`x − S(x)`**, which is
+//! **exactly zero inside the band** (`S(x) = x`): the control-law residual then
+//! collapses to `K·(y_ref − y) = 0 ⇒ y = y_ref` — **perfect steady-state
+//! tracking, independent of `K` and of where the actuator sits in its range**.
+//! (An earlier `x + S(x)` form was a droop/proportional controller with a
+//! residual steady-state offset `≈ 2x/K`; it is replaced by this offset-free
+//! form to remove the fragile gain tuning it required.)
+//!
+//! **Behaviour.** When the internal variable `x` sits inside `(−Q, Q)`,
+//! `S(x) = x`, so `x − S(x) = 0 ⇒ y = y_ref` (perfect tracking) and `u` is free
+//! between its bounds. When `x` leaves the band, `S(x)` saturates to `±Q`,
+//! `x − S(x) = ±(x − Q) ≠ 0`, pinning `u` to `u_max` / `u_min` while the
+//! tracking error is *released* (no wind-up): the solver naturally finds the
+//! best achievable `y` at the saturated actuator. This reproduces the
+//! anti-wind-up behaviour of a real supervisory PI controller at steady state.
+//!
+//! # Differentiability
+//!
+//! The exact `S(x)` is only C⁰ (kinks at `x = ±Q`), which hurts Newton
+//! convergence. [`Saturation::Smooth`] replaces `|·|` with the analytic
+//! [`smooth_abs`](entropyk_core::smoothing::smooth_abs) so the whole loop has a
+//! continuous Jacobian, at the cost of a small, tunable rounding of the corners.
+
+use entropyk_core::smoothing::{smooth_abs, smooth_abs_derivative};
+
+use super::bounded::BoundedVariableId;
+use super::constraint::ComponentOutput;
+use super::constraint::ConstraintId;
+use super::override_network::{eval_error_signal, eval_error_weights, Objective};
+
+/// Saturation-function variant used inside a [`SaturatedController`].
+#[derive(Debug, Clone, Copy, PartialEq)]
+pub enum Saturation {
+    /// Exact `S(x) = clamp(x, −Q, Q)`. Continuous but not differentiable at the
+    /// two corners `x = ±Q`.
+    Hard,
+    /// C¹ approximation using [`smooth_abs`] with rounding parameter `eps`
+    /// (larger `eps` ⇒ softer corners, smaller ⇒ closer to [`Saturation::Hard`]).
+    Smooth { eps: f64 },
+}
+
+impl Default for Saturation {
+    fn default() -> Self {
+        // A mild rounding that keeps the Jacobian continuous without materially
+        // shifting the saturation corners for typical Q≈1 loops.
+        Saturation::Smooth { eps: 1e-3 }
+    }
+}
+
+/// Errors raised while constructing a [`SaturatedController`].
+#[derive(Debug, Clone, PartialEq)]
+pub enum SaturatedControlError {
+    /// `u_max` was not strictly greater than `u_min`.
+    InvalidBounds { u_min: f64, u_max: f64 },
+    /// `Q` was not strictly positive.
+    NonPositiveQ(f64),
+    /// The gain `K` was zero (the loop would be inert).
+    ZeroGain,
+    /// A supplied parameter was not finite.
+    NonFinite(&'static str),
+}
+
+impl std::fmt::Display for SaturatedControlError {
+    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
+        match self {
+            SaturatedControlError::InvalidBounds { u_min, u_max } => write!(
+                f,
+                "actuator bounds invalid: u_max ({u_max}) must exceed u_min ({u_min})"
+            ),
+            SaturatedControlError::NonPositiveQ(q) => {
+                write!(f, "saturation band Q ({q}) must be strictly positive")
+            }
+            SaturatedControlError::ZeroGain => write!(f, "loop gain K must be non-zero"),
+            SaturatedControlError::NonFinite(p) => write!(f, "parameter `{p}` must be finite"),
+        }
+    }
+}
+
+impl std::error::Error for SaturatedControlError {}
+
+/// A single steady-state saturated-PI control loop.
+///
+/// It links a measurable plant output (`output`, the controlled variable `y`)
+/// to a bounded actuator (`actuator`, the manipulated variable `u`) so that the
+/// solver drives `y` to `setpoint` while respecting `[u_min, u_max]` with
+/// anti-wind-up. Each loop contributes two residuals and two unknowns (`u` and
+/// the internal variable `x`) to the global system.
+#[derive(Debug, Clone)]
+pub struct SaturatedController {
+    id: ConstraintId,
+    output: ComponentOutput,
+    actuator: BoundedVariableId,
+    setpoint: f64,
+    u_min: f64,
+    u_max: f64,
+    q: f64,
+    k: f64,
+    saturation: Saturation,
+    /// Optional override/selector network. When non-empty, the control law is
+    /// driven by the selected error `E` over these objectives instead of the
+    /// single `(output, setpoint, k)` triple above. See [`super::override_network`].
+    objectives: Vec<Objective>,
+    /// Selector smoothing sharpness for the override network (`softMin`/`softMax`).
+    alpha: f64,
+    /// Override **activation** homotopy parameter `λ ∈ [0, 1]`. The effective
+    /// selected error blends the primary objective with the full network:
+    /// `E_eff = (1−λ)·e_primary + λ·E`. At `λ = 0` only the primary objective is
+    /// active (a feasible baseline that always solves like the single-loop
+    /// controller); at `λ = 1` the full override network is in force. Ramping `λ`
+    /// from 0 → 1 with warm starts engages protections gradually and keeps every
+    /// intermediate solve near-feasible — the robust cure for hard cold-start
+    /// override switching. Default `1.0` (fully active).
+    activation: f64,
+}
+
+impl SaturatedController {
+    /// Builds a controller loop.
+    ///
+    /// * `output` — the controlled variable `y` (a measurable plant output).
+    /// * `actuator` — the manipulated variable `u` (a bounded control variable).
+    /// * `setpoint` — the target value `y_ref` for `y`.
+    /// * `u_min`, `u_max` — the actuator saturation bounds.
+    ///
+    /// Defaults: unit gain (`K = +1`), unit band (`Q = 1`), and the default
+    /// [`Saturation`]. Tune with [`Self::with_gain`], [`Self::with_band`] and
+    /// [`Self::with_saturation`].
+    pub fn new(
+        id: ConstraintId,
+        output: ComponentOutput,
+        actuator: BoundedVariableId,
+        setpoint: f64,
+        u_min: f64,
+        u_max: f64,
+    ) -> Result<Self, SaturatedControlError> {
+        if !setpoint.is_finite() {
+            return Err(SaturatedControlError::NonFinite("setpoint"));
+        }
+        if !u_min.is_finite() {
+            return Err(SaturatedControlError::NonFinite("u_min"));
+        }
+        if !u_max.is_finite() {
+            return Err(SaturatedControlError::NonFinite("u_max"));
+        }
+        if u_max <= u_min {
+            return Err(SaturatedControlError::InvalidBounds { u_min, u_max });
+        }
+        Ok(Self {
+            id,
+            output,
+            actuator,
+            setpoint,
+            u_min,
+            u_max,
+            q: 1.0,
+            k: 1.0,
+            saturation: Saturation::default(),
+            objectives: Vec::new(),
+            alpha: 1e-3,
+            activation: 1.0,
+        })
+    }
+
+    /// Sets the loop gain `K` (its sign must match the actuator→output
+    /// sensitivity; magnitude scales the internal variable).
+    pub fn with_gain(mut self, k: f64) -> Result<Self, SaturatedControlError> {
+        if !k.is_finite() {
+            return Err(SaturatedControlError::NonFinite("k"));
+        }
+        if k == 0.0 {
+            return Err(SaturatedControlError::ZeroGain);
+        }
+        self.k = k;
+        Ok(self)
+    }
+
+    /// Sets the saturation band half-width `Q > 0`.
+    pub fn with_band(mut self, q: f64) -> Result<Self, SaturatedControlError> {
+        if !q.is_finite() {
+            return Err(SaturatedControlError::NonFinite("q"));
+        }
+        if q <= 0.0 {
+            return Err(SaturatedControlError::NonPositiveQ(q));
+        }
+        self.q = q;
+        Ok(self)
+    }
+
+    /// Selects the saturation-function variant (hard vs. C¹ smooth).
+    pub fn with_saturation(mut self, saturation: Saturation) -> Self {
+        self.saturation = saturation;
+        self
+    }
+
+    /// Attaches an override/selector network of [`Objective`]s. When set, the
+    /// control-law residual is driven by the selected error `E` (a `softMin`/
+    /// `softMax` fold over the objectives) instead of the single `(output,
+    /// setpoint, gain)` triple. The primary setpoint should be the first
+    /// objective; higher-priority protections come later in the chain.
+    ///
+    /// The per-objective `gain`s carry the normalization and sign, so the loop
+    /// gain `K` is implicitly `1` in network mode.
+    pub fn with_objectives(mut self, objectives: Vec<Objective>) -> Self {
+        self.objectives = objectives;
+        self
+    }
+
+    /// Sets the selector smoothing sharpness `alpha > 0` for the override
+    /// network (smaller ⇒ sharper `min`/`max`, larger ⇒ smoother/more robust).
+    pub fn with_alpha(mut self, alpha: f64) -> Self {
+        if alpha.is_finite() && alpha > 0.0 {
+            self.alpha = alpha;
+        }
+        self
+    }
+
+    /// In-place setter for the selector smoothing sharpness `alpha`. Used by the
+    /// warm-started **alpha-continuation** (homotopy on the selector sharpness):
+    /// solve first with a large `alpha` (very smooth, well-conditioned), then
+    /// anneal toward the target while warm-starting from the previous solution.
+    /// Ignores non-finite or non-positive values.
+    pub fn set_alpha(&mut self, alpha: f64) {
+        if alpha.is_finite() && alpha > 0.0 {
+            self.alpha = alpha;
+        }
+    }
+
+    /// Whether this controller uses an override/selector network.
+    pub fn is_network(&self) -> bool {
+        !self.objectives.is_empty()
+    }
+
+    /// The override objectives (empty in single-objective mode).
+    pub fn objectives(&self) -> &[Objective] {
+        &self.objectives
+    }
+
+    /// Selector smoothing sharpness for the override network.
+    pub fn alpha(&self) -> f64 {
+        self.alpha
+    }
+
+    /// Raw selected error `E` over the override network for the given measured
+    /// objective values (`measured[i]` ↔ `objectives()[i]`), ignoring activation.
+    pub fn error_signal(&self, measured: &[f64]) -> f64 {
+        eval_error_signal(&self.objectives, measured, self.alpha)
+    }
+
+    /// Raw selector weights `w_i = ∂E/∂e_i` for the override network, ignoring
+    /// activation.
+    pub fn error_weights(&self, measured: &[f64]) -> Vec<f64> {
+        eval_error_weights(&self.objectives, measured, self.alpha)
+    }
+
+    /// Override activation homotopy parameter `λ ∈ [0, 1]`.
+    pub fn activation(&self) -> f64 {
+        self.activation
+    }
+
+    /// In-place setter for the override activation `λ` (clamped to `[0, 1]`),
+    /// used by the warm-started activation continuation.
+    pub fn set_activation(&mut self, activation: f64) {
+        if activation.is_finite() {
+            self.activation = activation.clamp(0.0, 1.0);
+        }
+    }
+
+    /// **Activation-blended** selected error
+    /// `E_eff = (1−λ)·e_primary + λ·E`. Equal to the raw network error at
+    /// `λ = 1`, and to the primary objective's error alone at `λ = 0`.
+    pub fn network_error(&self, measured: &[f64]) -> f64 {
+        let full = eval_error_signal(&self.objectives, measured, self.alpha);
+        if self.activation >= 1.0 || self.objectives.is_empty() {
+            return full;
+        }
+        let primary = self.objectives[0].error(measured[0]);
+        (1.0 - self.activation) * primary + self.activation * full
+    }
+
+    /// Activation-blended selector weights `∂E_eff/∂e_i`.
+    pub fn network_error_weights(&self, measured: &[f64]) -> Vec<f64> {
+        let mut w = eval_error_weights(&self.objectives, measured, self.alpha);
+        if self.activation >= 1.0 || w.is_empty() {
+            return w;
+        }
+        for wi in w.iter_mut() {
+            *wi *= self.activation;
+        }
+        w[0] += 1.0 - self.activation;
+        w
+    }
+
+    /// Control-law residual in override-network mode:
+    /// `r_y = E − (x − S(x))`, where `E` is the selected error.
+    pub fn residual_y_network(&self, e: f64, x: f64) -> f64 {
+        e - (x - self.saturation_s(x))
+    }
+
+    /// The controller's identifier.
+    pub fn id(&self) -> &ConstraintId {
+        &self.id
+    }
+
+    /// The controlled plant output `y`.
+    pub fn output(&self) -> &ComponentOutput {
+        &self.output
+    }
+
+    /// The manipulated actuator `u`.
+    pub fn actuator(&self) -> &BoundedVariableId {
+        &self.actuator
+    }
+
+    /// The output setpoint `y_ref`.
+    pub fn setpoint(&self) -> f64 {
+        self.setpoint
+    }
+
+    /// Actuator lower bound `u_min`.
+    pub fn u_min(&self) -> f64 {
+        self.u_min
+    }
+
+    /// Actuator upper bound `u_max`.
+    pub fn u_max(&self) -> f64 {
+        self.u_max
+    }
+
+    /// `Z = (u_max − u_min) / 2` — half the actuator span.
+    pub fn z(&self) -> f64 {
+        0.5 * (self.u_max - self.u_min)
+    }
+
+    /// `Y = (u_max + u_min) / 2` — actuator mid-point.
+    pub fn y_offset(&self) -> f64 {
+        0.5 * (self.u_max + self.u_min)
+    }
+
+    /// Loop gain `K` used by the control-law residual.
+    pub fn gain(&self) -> f64 {
+        self.k
+    }
+
+    /// Saturation function `S(x)` (exact or smoothed per [`Saturation`]).
+    ///
+    /// Exact form: `S(x) = (|x + Q| − |x − Q|)/2 = clamp(x, −Q, Q)`.
+    pub fn saturation_s(&self, x: f64) -> f64 {
+        match self.saturation {
+            Saturation::Hard => x.clamp(-self.q, self.q),
+            Saturation::Smooth { eps } => {
+                0.5 * (smooth_abs(x + self.q, eps) - smooth_abs(x - self.q, eps))
+            }
+        }
+    }
+
+    /// Derivative `S'(x)` of [`Self::saturation_s`].
+    pub fn saturation_ds(&self, x: f64) -> f64 {
+        match self.saturation {
+            Saturation::Hard => {
+                if x > -self.q && x < self.q {
+                    1.0
+                } else {
+                    0.0
+                }
+            }
+            Saturation::Smooth { eps } => {
+                0.5 * (smooth_abs_derivative(x + self.q, eps)
+                    - smooth_abs_derivative(x - self.q, eps))
+            }
+        }
+    }
+
+    /// Actuator-law residual `r_u = u − (Z·S(x) + Y)`.
+    ///
+    /// Zero when the actuator equals the value dictated by the internal
+    /// variable through the saturation map.
+    pub fn residual_u(&self, u: f64, x: f64) -> f64 {
+        u - (self.z() * self.saturation_s(x) + self.y_offset())
+    }
+
+    /// Control-law residual `r_y = K·(y_ref − y) − (x − S(x))`
+    /// (offset-free / integral-equivalent form, Mancuso & Morari eq. A.1).
+    ///
+    /// Inside the band `x − S(x) = 0`, so the residual is `K·(y_ref − y)` and
+    /// vanishes exactly at `y = y_ref` (perfect steady-state tracking). Once `x`
+    /// (hence `u`) saturates, `x − S(x) ≠ 0` releases the tracking error
+    /// smoothly (anti-wind-up).
+    pub fn residual_y(&self, y: f64, x: f64) -> f64 {
+        self.k * (self.setpoint - y) - (x - self.saturation_s(x))
+    }
+
+    /// Jacobian partials of the actuator-law residual `r_u`.
+    ///
+    /// Returns `(∂r_u/∂u, ∂r_u/∂x)`. The dependence on the plant state enters
+    /// only through `u` and `x`, both system unknowns.
+    pub fn d_residual_u(&self, x: f64) -> (f64, f64) {
+        (1.0, -self.z() * self.saturation_ds(x))
+    }
+
+    /// Jacobian partials of the control-law residual `r_y`.
+    ///
+    /// Returns `(∂r_y/∂y, ∂r_y/∂x)`. `∂r_y/∂y = −K` is the coupling into the
+    /// plant (via whatever state `y` is measured from); `∂r_y/∂x = −(1 − S'(x))`.
+    ///
+    /// In the unsaturated band `S'(x) = 1`, so `∂r_y/∂x = 0`: the control-law
+    /// row constrains only `y` (`y = y_ref`) and is closed through the plant
+    /// coupling `∂r_y/∂y = −K` (well-posed whenever the loop is controllable,
+    /// `∂y/∂u ≠ 0`). Under saturation `S'(x) → 0`, restoring `∂r_y/∂x → −1`.
+    pub fn d_residual_y(&self, x: f64) -> (f64, f64) {
+        (-self.k, -(1.0 - self.saturation_ds(x)))
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    fn ctrl(saturation: Saturation) -> SaturatedController {
+        SaturatedController::new(
+            ConstraintId::new("lwt"),
+            ComponentOutput::Temperature {
+                component_id: "evaporator".to_string(),
+            },
+            BoundedVariableId::new("compressor_f_m"),
+            280.0, // y_ref
+            0.5,   // u_min
+            2.0,   // u_max
+        )
+        .unwrap()
+        .with_saturation(saturation)
+    }
+
+    #[test]
+    fn test_saturation_equals_hard_clamp() {
+        let c = ctrl(Saturation::Hard).with_band(1.0).unwrap();
+        assert_eq!(c.saturation_s(0.4), 0.4);
+        assert_eq!(c.saturation_s(2.0), 1.0);
+        assert_eq!(c.saturation_s(-3.0), -1.0);
+        // Exact analytic identity S(x) = (|x+Q|-|x-Q|)/2.
+        for &x in &[-2.0f64, -0.7, 0.0, 0.3, 5.0] {
+            let exact = 0.5 * ((x + 1.0).abs() - (x - 1.0).abs());
+            assert!((c.saturation_s(x) - exact).abs() < 1e-12);
+        }
+    }
+
+    #[test]
+    fn test_z_and_y_offset() {
+        let c = ctrl(Saturation::Hard);
+        assert!((c.z() - 0.75).abs() < 1e-12); // (2.0-0.5)/2
+        assert!((c.y_offset() - 1.25).abs() < 1e-12); // (2.0+0.5)/2
+    }
+
+    #[test]
+    fn test_unsaturated_loop_tracks_setpoint() {
+        // Offset-free: in the linear band S(x)=x ⇒ x−S(x)=0, so r_y=0 collapses
+        // to K(y_ref−y)=0 ⇒ y=y_ref for ANY x in (−Q,Q), and r_u=0 ⇒ u=Z·x+Y.
+        let c = ctrl(Saturation::Hard).with_band(1.0).unwrap();
+        let x = 0.3;
+        // Perfect tracking, independent of the internal variable x.
+        let y = c.setpoint();
+        let u = c.z() * x + c.y_offset();
+        assert!(c.residual_y(y, x).abs() < 1e-12);
+        assert!(c.residual_u(u, x).abs() < 1e-12);
+        // Actuator lands strictly inside its bounds (not saturated).
+        assert!(u > c.u_min() && u < c.u_max());
+    }
+
+    #[test]
+    fn test_saturated_actuator_releases_tracking_error() {
+        // Drive the internal variable well past the band: the actuator law must
+        // pin u to u_max and the control law must NOT force y = y_ref (the error
+        // is released — anti-wind-up).
+        let c = ctrl(Saturation::Hard).with_band(1.0).unwrap();
+        let x = 10.0; // deep in saturation, S(x)=Q=1
+        let u = c.z() * c.saturation_s(x) + c.y_offset();
+        assert!((u - c.u_max()).abs() < 1e-12, "u must pin to u_max: {u}");
+        // At r_y = 0: K(y_ref − y) = x − S(x) = 9 ⇒ y = y_ref − 9 ≠ y_ref.
+        let y = c.setpoint() - (x - c.saturation_s(x));
+        assert!(c.residual_y(y, x).abs() < 1e-12);
+        assert!(
+            (y - c.setpoint()).abs() > 1.0,
+            "tracking error must be released under saturation"
+        );
+    }
+
+    #[test]
+    fn test_smooth_saturation_is_c1_and_close_to_hard() {
+        let hard = ctrl(Saturation::Hard).with_band(1.0).unwrap();
+        let soft = ctrl(Saturation::Smooth { eps: 1e-3 })
+            .with_band(1.0)
+            .unwrap();
+        // Away from the corners the smooth and hard forms nearly coincide.
+        for &x in &[-3.0, -0.5, 0.0, 0.5, 3.0] {
+            assert!(
+                (soft.saturation_s(x) - hard.saturation_s(x)).abs() < 5e-3,
+                "smooth S deviates too far at x={x}"
+            );
+        }
+        // Finite-difference check of the analytic smooth derivative.
+        let x = 0.9;
+        let h = 1e-6;
+        let fd = (soft.saturation_s(x + h) - soft.saturation_s(x - h)) / (2.0 * h);
+        assert!(
+            (soft.saturation_ds(x) - fd).abs() < 1e-4,
+            "S'(x) mismatch: analytic {} vs FD {}",
+            soft.saturation_ds(x),
+            fd
+        );
+    }
+
+    #[test]
+    fn test_residual_jacobians_match_finite_difference() {
+        let c = ctrl(Saturation::Smooth { eps: 1e-3 })
+            .with_band(1.0)
+            .unwrap()
+            .with_gain(-1.5)
+            .unwrap();
+        let (x, u, y) = (0.6, 1.1, 279.0);
+        let h = 1e-6;
+
+        // ∂r_u/∂u and ∂r_u/∂x
+        let (dru_du, dru_dx) = c.d_residual_u(x);
+        let fd_dru_du = (c.residual_u(u + h, x) - c.residual_u(u - h, x)) / (2.0 * h);
+        let fd_dru_dx = (c.residual_u(u, x + h) - c.residual_u(u, x - h)) / (2.0 * h);
+        assert!((dru_du - fd_dru_du).abs() < 1e-5);
+        assert!((dru_dx - fd_dru_dx).abs() < 1e-4);
+
+        // ∂r_y/∂y and ∂r_y/∂x
+        let (dry_dy, dry_dx) = c.d_residual_y(x);
+        let fd_dry_dy = (c.residual_y(y + h, x) - c.residual_y(y - h, x)) / (2.0 * h);
+        let fd_dry_dx = (c.residual_y(y, x + h) - c.residual_y(y, x - h)) / (2.0 * h);
+        assert!((dry_dy - fd_dry_dy).abs() < 1e-5);
+        assert!((dry_dx - fd_dry_dx).abs() < 1e-4);
+    }
+
+    #[test]
+    fn test_construction_validates_parameters() {
+        let out = ComponentOutput::Temperature {
+            component_id: "e".to_string(),
+        };
+        // u_max <= u_min rejected.
+        assert!(matches!(
+            SaturatedController::new(
+                ConstraintId::new("c"),
+                out.clone(),
+                BoundedVariableId::new("a"),
+                1.0,
+                2.0,
+                2.0
+            ),
+            Err(SaturatedControlError::InvalidBounds { .. })
+        ));
+        // Zero gain rejected.
+        let c = SaturatedController::new(
+            ConstraintId::new("c"),
+            out.clone(),
+            BoundedVariableId::new("a"),
+            1.0,
+            0.0,
+            1.0,
+        )
+        .unwrap();
+        assert!(matches!(
+            c.with_gain(0.0),
+            Err(SaturatedControlError::ZeroGain)
+        ));
+        // Non-positive Q rejected.
+        let c2 = SaturatedController::new(
+            ConstraintId::new("c"),
+            out,
+            BoundedVariableId::new("a"),
+            1.0,
+            0.0,
+            1.0,
+        )
+        .unwrap();
+        assert!(matches!(
+            c2.with_band(0.0),
+            Err(SaturatedControlError::NonPositiveQ(_))
+        ));
+    }
+}

crates/solver/src/jacobian.rs

@@ -103,8 +103,17 @@ impl JacobianMatrix {
 
     /// Solves the linear system J·Δx = -r and returns Δx.
     ///
-    /// Uses LU decomposition with partial pivoting. Returns `None` if the
-    /// matrix is singular (no unique solution exists).
+    /// Uses **Ruiz equilibration** (iterative row/column scaling) followed by LU
+    /// decomposition with partial pivoting. Equilibration rescales the rows and
+    /// columns so their ∞-norms approach 1, which dramatically lowers the
+    /// condition number of badly-scaled Jacobians — exactly the situation that
+    /// arises when a thermodynamic system mixes pressures (~1e6), enthalpies
+    /// (~1e5) and dimensionless controls (~1) in the same matrix. The scaling is
+    /// solution-preserving (it is undone on the returned step), so the result is
+    /// mathematically identical to an unscaled solve but numerically far more
+    /// robust on stiff, >50-variable systems.
+    ///
+    /// Returns `None` if the matrix is singular (no unique solution exists).
     ///
     /// # Arguments
     ///
@@ -138,16 +147,45 @@ impl JacobianMatrix {
             return None;
         }
 
-        // For square systems, use LU decomposition
+        // For square systems, use Ruiz-equilibrated LU decomposition.
         if self.0.nrows() == self.0.ncols() {
-            let lu = self.0.clone().lu();
+            let n = self.0.nrows();
 
-            // Solve J·Δx = -r
-            let r_vec = DVector::from_row_slice(residuals);
-            let neg_r = -r_vec;
+            // Diagonal scalings D_r, D_c such that the scaled matrix
+            // Ĵ = diag(d_r) · J · diag(d_c) has near-unit row/column ∞-norms.
+            let (d_r, d_c) = crate::scaling::equilibrate(&self.0);
 
-            match lu.solve(&neg_r) {
-                Some(delta) => Some(delta.iter().copied().collect()),
+            // Build the scaled matrix in place on a single clone (same number of
+            // allocations as the previous unscaled path).
+            let mut scaled = self.0.clone();
+            for i in 0..n {
+                for j in 0..n {
+                    scaled[(i, j)] *= d_r[i] * d_c[j];
+                }
+            }
+
+            let lu = scaled.lu();
+
+            // Scaled right-hand side: D_r · (-r).
+            let neg_r_scaled = DVector::from_iterator(n, (0..n).map(|i| -residuals[i] * d_r[i]));
+
+            match lu.solve(&neg_r_scaled) {
+                // Undo the column scaling to recover the true step: Δx = D_c · y.
+                Some(y) => {
+                    let delta = crate::scaling::unscale_dx(y.as_slice(), &d_c);
+                    // A NaN/Inf entry anywhere in the Jacobian (or RHS) silently
+                    // slips through LU as a non-finite step. Reject it here so the
+                    // caller treats the iteration as a clean failure instead of
+                    // propagating NaN into the state for another iteration.
+                    if delta.iter().all(|v| v.is_finite()) {
+                        Some(delta)
+                    } else {
+                        tracing::warn!(
+                            "LU solve produced a non-finite step - Jacobian may contain NaN/Inf"
+                        );
+                        None
+                    }
+                }
                 None => {
                     tracing::warn!("LU solve failed - Jacobian may be singular");
                     None
@@ -168,7 +206,17 @@ impl JacobianMatrix {
             // Use SVD for robust least-squares solution
             let svd = self.0.clone().svd(true, true);
             match svd.solve(&neg_r, 1e-10) {
-                Ok(delta) => Some(delta.iter().copied().collect()),
+                Ok(delta) => {
+                    let v: Vec<f64> = delta.iter().copied().collect();
+                    if v.iter().all(|x| x.is_finite()) {
+                        Some(v)
+                    } else {
+                        tracing::warn!(
+                            "SVD solve produced a non-finite step - Jacobian may contain NaN/Inf"
+                        );
+                        None
+                    }
+                }
                 Err(e) => {
                     tracing::warn!("SVD solve failed - Jacobian may be rank-deficient: {}", e);
                     None
@@ -177,6 +225,29 @@ impl JacobianMatrix {
         }
     }
 
+    /// Returns the Ruiz row/column equilibration factors `(d_r, d_c)`.
+    ///
+    /// The scaled matrix `diag(d_r) · J · diag(d_c)` has row and column
+    /// ∞-norms close to 1. This is the same scaling applied internally by
+    /// [`solve`](Self::solve); it is exposed for diagnostics (e.g. inspecting how
+    /// ill-scaled a Jacobian is, or estimating the conditioning improvement).
+    ///
+    /// # Example
+    ///
+    /// ```rust
+    /// use entropyk_solver::jacobian::JacobianMatrix;
+    ///
+    /// // Badly scaled diagonal: entries span 12 orders of magnitude.
+    /// let entries = vec![(0, 0, 1e6), (1, 1, 1e-6)];
+    /// let j = JacobianMatrix::from_builder(&entries, 2, 2);
+    /// let (d_r, d_c) = j.ruiz_scaling_factors();
+    /// assert_eq!(d_r.len(), 2);
+    /// assert_eq!(d_c.len(), 2);
+    /// ```
+    pub fn ruiz_scaling_factors(&self) -> (Vec<f64>, Vec<f64>) {
+        crate::scaling::equilibrate(&self.0)
+    }
+
     /// Estimates the condition number of the Jacobian matrix.
     ///
     /// The condition number κ = σ_max / σ_min indicates how ill-conditioned
@@ -225,7 +296,11 @@ impl JacobianMatrix {
         }
 
         let sigma_max = singular_values.max();
-        let sigma_min = singular_values.iter().filter(|&&s| s > 0.0).min_by(|a, b| a.partial_cmp(b).unwrap()).copied();
+        let sigma_min = singular_values
+            .iter()
+            .filter(|&&s| s > 0.0)
+            .min_by(|a, b| a.partial_cmp(b).unwrap())
+            .copied();
 
         match sigma_min {
             Some(min) => Some(sigma_max / min),
@@ -471,6 +546,15 @@ mod tests {
     use super::*;
     use approx::assert_relative_eq;
 
+    /// A NaN entry in the Jacobian must yield `None` (clean failure) rather than
+    /// a `Some(..)` step containing NaN that would poison the next iteration.
+    #[test]
+    fn test_solve_with_nan_entry_returns_none() {
+        let entries = vec![(0, 0, f64::NAN), (1, 1, 1.0)];
+        let j = JacobianMatrix::from_builder(&entries, 2, 2);
+        assert!(j.solve(&[1.0, 1.0]).is_none());
+    }
+
     #[test]
     fn test_from_builder_simple() {
         let entries = vec![(0, 0, 1.0), (0, 1, 2.0), (1, 0, 3.0), (1, 1, 4.0)];
@@ -694,4 +778,82 @@ mod tests {
         assert_relative_eq!(j_num.get(1, 0).unwrap(), j10, epsilon = 1e-5);
         assert_relative_eq!(j_num.get(1, 1).unwrap(), j11, epsilon = 1e-5);
     }
+
+    #[test]
+    fn test_ruiz_equilibration_unit_norms() {
+        // After equilibration, scaled row/column ∞-norms should be ≈ 1.
+        let entries = vec![(0, 0, 1e6), (0, 1, 1e3), (1, 0, 1e-3), (1, 1, 1e-6)];
+        let j = JacobianMatrix::from_builder(&entries, 2, 2);
+        let (d_r, d_c) = j.ruiz_scaling_factors();
+
+        let m = j.as_matrix();
+        for i in 0..2 {
+            let row_max = (0..2)
+                .map(|jj| (d_r[i] * m[(i, jj)] * d_c[jj]).abs())
+                .fold(0.0_f64, f64::max);
+            assert!(
+                (row_max - 1.0).abs() < 0.05,
+                "row {} ∞-norm not unit: {}",
+                i,
+                row_max
+            );
+        }
+        for jj in 0..2 {
+            let col_max = (0..2)
+                .map(|i| (d_r[i] * m[(i, jj)] * d_c[jj]).abs())
+                .fold(0.0_f64, f64::max);
+            assert!(
+                (col_max - 1.0).abs() < 0.05,
+                "col {} ∞-norm not unit: {}",
+                jj,
+                col_max
+            );
+        }
+    }
+
+    #[test]
+    fn test_ruiz_reduces_condition_number() {
+        // Badly-scaled but SVD-resolvable matrix (dynamic range ~1e6 keeps the
+        // small singular value above underflow, so the condition number is
+        // measured reliably). Row 0 lives at ~1e6, row 1 at ~1.
+        let entries = vec![(0, 0, 1.0e6), (0, 1, 2.0e6), (1, 0, 3.0), (1, 1, 1.0)];
+        let j = JacobianMatrix::from_builder(&entries, 2, 2);
+
+        let cond_before = j.estimate_condition_number().unwrap();
+
+        // Apply the Ruiz scaling and recompute the condition number.
+        let (d_r, d_c) = j.ruiz_scaling_factors();
+        let m = j.as_matrix();
+        let mut scaled_entries = Vec::new();
+        for i in 0..2 {
+            for jj in 0..2 {
+                scaled_entries.push((i, jj, d_r[i] * m[(i, jj)] * d_c[jj]));
+            }
+        }
+        let scaled = JacobianMatrix::from_builder(&scaled_entries, 2, 2);
+        let cond_after = scaled.estimate_condition_number().unwrap();
+
+        assert!(
+            cond_after < cond_before / 100.0 && cond_after < 10.0,
+            "Ruiz should slash the condition number: before={:.3e}, after={:.3e}",
+            cond_before,
+            cond_after
+        );
+    }
+
+    #[test]
+    fn test_solve_illconditioned_matches_known_solution() {
+        // Badly-scaled system with a known solution. J·x = b where
+        // J = [[1e6, 2e6], [3e-6, 4e-6]], x = [1, -1]  =>  b = [-1e6, -1e-6].
+        // solve() returns Δx for J·Δx = -r, so set r = -b to get Δx = x.
+        let entries = vec![(0, 0, 1e6), (0, 1, 2e6), (1, 0, 3e-6), (1, 1, 4e-6)];
+        let j = JacobianMatrix::from_builder(&entries, 2, 2);
+
+        let b = [-1e6, -1e-6];
+        let r = [-b[0], -b[1]];
+        let delta = j.solve(&r).expect("non-singular");
+
+        assert_relative_eq!(delta[0], 1.0, epsilon = 1e-9);
+        assert_relative_eq!(delta[1], -1.0, epsilon = 1e-9);
+    }
 }

crates/solver/src/lib.rs

@@ -8,6 +8,7 @@
 
 pub mod coupling;
 pub mod criteria;
+pub mod dof;
 pub mod error;
 pub mod graph;
 pub mod initializer;
@@ -15,36 +16,44 @@ pub mod inverse;
 pub mod jacobian;
 pub mod macro_component;
 pub mod metadata;
+pub mod scaling;
 pub mod snapshot;
 pub mod snapshot_params;
 pub mod solver;
 pub mod strategies;
 pub mod system;
+pub mod topology;
 
 pub use coupling::{
     compute_coupling_heat, coupling_groups, has_circular_dependencies, ThermalCoupling,
 };
+pub use dof::{
+    align_roles, unspecified_roles, ComponentEquationBlock, DofReport, EquationRole,
+    SystemDofBalance, SystemDofError, UnknownKind,
+};
 pub use criteria::{CircuitConvergence, ConvergenceCriteria, ConvergenceReport};
 pub use entropyk_components::ConnectionError;
 pub use entropyk_core::CircuitId;
 pub use error::{AddEdgeError, ThermoError, TopologyError};
 pub use initializer::{
-    antoine_pressure, AntoineCoefficients, InitializerConfig, InitializerError, SmartInitializer,
+    antoine_pressure, AntoineCoefficients, InitializationDiagnostics, InitializationRegime,
+    InitializationSeed, InitializerConfig, InitializerError, SmartInitializer, StartValues,
 };
 pub use inverse::{ComponentOutput, Constraint, ConstraintError, ConstraintId};
 pub use jacobian::JacobianMatrix;
 pub use macro_component::{MacroComponent, MacroComponentSnapshot, PortMapping};
 pub use metadata::SimulationMetadata;
+pub use scaling::{equilibrate, unscale_dx};
 pub use snapshot::{
     BoundedVariableSnapshot, ConstraintSnapshot, EdgeSnapshot, FluidBackendInfo,
     SolverConfigSnapshot, SystemSnapshot, TopologySnapshot,
 };
 pub use solver::{
-    ConvergedState, ConvergenceStatus, ConvergenceDiagnostics, IterationDiagnostics,
-    JacobianFreezingConfig, Solver, SolverError, SolverSwitchEvent, SolverType, SwitchReason,
-    TimeoutConfig, VerboseConfig, VerboseOutputFormat,
+    dominant_residual, ConvergedState, ConvergenceDiagnostics, ConvergenceStatus,
+    IterationDiagnostics, JacobianFreezingConfig, Solver, SolverError, SolverSwitchEvent,
+    SolverType, SwitchReason, TimeoutConfig, VerboseConfig, VerboseOutputFormat,
 };
 pub use strategies::{
-    FallbackConfig, FallbackSolver, NewtonConfig, PicardConfig, SolverStrategy,
+    FallbackConfig, FallbackSolver, HomotopyConfig, NewtonConfig, PicardConfig, SolverStrategy,
 };
-pub use system::{FlowEdge, System, MAX_CIRCUIT_ID};
+pub use system::{CyclePerformance, FlowEdge, System, MAX_CIRCUIT_ID};

crates/solver/src/macro_component.rs

@@ -26,13 +26,17 @@
 //!
 //! When the `MacroComponent` is connected to external edges in the parent graph,
 //! coupling residuals enforce continuity between those external edges and the
-//! corresponding exposed internal edges:
+//! corresponding exposed internal edges.  With the stride-3 `(ṁ, P, h)` layout
+//! introduced in Story CM1.2, P is at `offset + 3 * pos + 1` and h at
+//! `offset + 3 * pos + 2` (resolved symbolically via `edge_state_indices()`):
 //!
 //! ```text
-//! r_P = state[p_ext] − state[offset + 2 * internal_edge_pos]     = 0
-//! r_h = state[h_ext] − state[offset + 2 * internal_edge_pos + 1] = 0
+//! r_P = state[p_ext] − state[offset + p_local] = 0
+//! r_h = state[h_ext] − state[offset + h_local] = 0
 //! ```
 //!
+//! Note: ṁ continuity at ports is not yet enforced (deferred to Story CM1.3).
+//!
 //! ## Serialization (AC #4)
 //!
 //! The full component graph inside a `MacroComponent` cannot be trivially
@@ -82,7 +86,7 @@ pub struct PortMapping {
 /// ```json
 /// {
 ///   "label": "chiller_1",
-///   "internal_edge_states": [1.5e5, 4.2e5, 8.0e4, 3.8e5],
+///   "internal_edge_states": [0.05, 1.5e5, 4.2e5, 0.05, 8.0e4, 3.8e5],
 ///   "port_names": ["evap_in", "evap_out"]
 /// }
 /// ```
@@ -90,7 +94,9 @@ pub struct PortMapping {
 pub struct MacroComponentSnapshot {
     /// Optional human-readable label for the subsystem.
     pub label: Option<String>,
-    /// Flat state vector for the internal edges: `[P_e0, h_e0, P_e1, h_e1, ...]`.
+    /// Flat state vector for the internal edges: `[ṁ_e0, P_e0, h_e0, ṁ_e1, P_e1, h_e1, ...]`.
+    ///
+    /// Per-edge layout is `(ṁ, P, h)` (stride 3, introduced in Story CM1.2).
     pub internal_edge_states: Vec<f64>,
     /// Names of exposed ports, in the same order as `port_mappings`.
     pub port_names: Vec<String>,
@@ -113,9 +119,16 @@ pub struct MacroComponentSnapshot {
 /// `compute_residuals` then appends 2 coupling residuals per exposed port:
 ///
 /// ```text
-/// r_P[i] = state[p_ext_i] − state[offset + 2·internal_edge_pos_i]     = 0
-/// r_h[i] = state[h_ext_i] − state[offset + 2·internal_edge_pos_i + 1] = 0
+/// r_P[i] = state[p_ext_i] − state[offset + p_local_i] = 0
+/// r_h[i] = state[h_ext_i] − state[offset + h_local_i] = 0
 /// ```
+/// (where `p_local_i`, `h_local_i` are the internal-local state indices of the
+/// mapped internal edge — stride 3 in the `(ṁ, P, h)` layout)
+///
+/// TEMP (CM1.2): ṁ continuity at ports is not yet coupled (`r_m = ṁ_ext − ṁ_int`
+/// is missing). Each side has its own independent mass-flow closure (`ṁ = seed`).
+/// Story CM1.3 must add the ṁ coupling residual + Jacobian entries here, and
+/// update `n_equations()` from `2 * n_ports` to `3 * n_ports`.
 ///
 /// # Usage
 ///
@@ -151,10 +164,10 @@ pub struct MacroComponent {
     /// internal edge. Set automatically via `set_system_context` during parent
     /// `System::finalize()`.  Defaults to 0.
     global_state_offset: usize,
-    /// State indices `(p_idx, h_idx)` of every parent-graph edge incident to
+    /// State indices `(m_idx, p_idx, h_idx)` of every parent-graph edge incident to
     /// this node (incoming and outgoing), in traversal order.
     /// Populated by `set_system_context`; empty until finalization.
-    external_edge_state_indices: Vec<(usize, usize)>,
+    external_edge_state_indices: Vec<(usize, usize, usize)>,
 }
 
 impl MacroComponent {
@@ -239,12 +252,12 @@ impl MacroComponent {
         &self.port_mappings
     }
 
-    /// Number of internal edges (each contributes 2 state variables: P, h).
+    /// Number of internal edges (each contributes 3 state variables: ṁ, P, h).
     pub fn internal_edge_count(&self) -> usize {
         self.internal.edge_count()
     }
 
-    /// Total number of internal state variables (2 per edge).
+    /// Total number of internal state variables (3 per edge: ṁ, P, h).
     pub fn internal_state_len(&self) -> usize {
         self.internal.state_vector_len()
     }
@@ -259,6 +272,24 @@ impl MacroComponent {
             .sum()
     }
 
+    /// Number of internal residuals produced by `internal.compute_residuals`:
+    /// component equations plus constraint and coupling rows.
+    fn n_internal_residuals(&self) -> usize {
+        self.n_internal_equations()
+            + self.internal.constraint_residual_count()
+            + self.internal.coupling_residual_count()
+    }
+
+    /// Internal-local `(p_idx, h_idx)` of the internal edge at graph-order
+    /// position `pos`, accounting for the per-edge stride (CM1.2: `(ṁ, P, h)`).
+    /// Returns `None` if no such internal edge exists.
+    fn internal_edge_ph_local(&self, pos: usize) -> Option<(usize, usize)> {
+        self.internal
+            .edge_indices()
+            .nth(pos)
+            .map(|e| self.internal.edge_state_indices(e))
+    }
+
     // ─── snapshot ─────────────────────────────────────────────────────────────
 
     /// Captures the current internal state as a serializable snapshot.
@@ -292,14 +323,14 @@ impl Component for MacroComponent {
     /// Called by `System::finalize()` to inject the parent-level state offset
     /// and the external edge state indices for this MacroComponent node.
     ///
-    /// `external_edge_state_indices` contains one `(p_idx, h_idx)` pair per
+    /// `external_edge_state_indices` contains one `(m_idx, p_idx, h_idx)` triple per
     /// parent edge incident to this node (in traversal order: incoming, then
     /// outgoing).  The *i*-th entry is matched to `port_mappings[i]` when
     /// emitting coupling residuals.
     fn set_system_context(
         &mut self,
         state_offset: usize,
-        external_edge_state_indices: &[(usize, usize)],
+        external_edge_state_indices: &[(usize, usize, usize)],
     ) {
         self.global_state_offset = state_offset;
         self.external_edge_state_indices = external_edge_state_indices.to_vec();
@@ -326,9 +357,9 @@ impl Component for MacroComponent {
             });
         }
 
-        let n_int_eqs = self.n_internal_equations();
+        let n_int_res = self.n_internal_residuals();
         let n_coupling = 2 * self.port_mappings.len();
-        let n_total = n_int_eqs + n_coupling;
+        let n_total = n_int_res + n_coupling;
 
         if residuals.len() < n_total {
             return Err(ComponentError::InvalidResidualDimensions {
@@ -339,22 +370,28 @@ impl Component for MacroComponent {
 
         // --- 1. Delegate internal residuals ----------------------------------
         let internal_state: Vec<f64> = state[start..end].to_vec();
-        let mut internal_residuals = vec![0.0; n_int_eqs];
+        let mut internal_residuals = vec![0.0; n_int_res];
         self.internal
             .compute_residuals(&internal_state, &mut internal_residuals)?;
-        residuals[..n_int_eqs].copy_from_slice(&internal_residuals);
+        residuals[..n_int_res].copy_from_slice(&internal_residuals);
 
         // --- 2. Port-coupling residuals --------------------------------------
         // For each exposed port mapping we append two residuals that enforce
         // continuity between the parent-graph external edge and the
         // corresponding internal edge:
         //
-        //   r_P = state[p_ext] − state[offset + 2·internal_edge_pos]     = 0
-        //   r_h = state[h_ext] − state[offset + 2·internal_edge_pos + 1] = 0
+        //   r_P = state[p_ext] − state[offset + p_local]  = 0
+        //   r_h = state[h_ext] − state[offset + h_local]  = 0
         for (i, mapping) in self.port_mappings.iter().enumerate() {
-            if let Some(&(p_ext, h_ext)) = self.external_edge_state_indices.get(i) {
-                let int_p = self.global_state_offset + 2 * mapping.internal_edge_pos;
-                let int_h = int_p + 1;
+            if let Some(&(_, p_ext, h_ext)) = self.external_edge_state_indices.get(i) {
+                let (p_local, h_local) = self
+                    .internal_edge_ph_local(mapping.internal_edge_pos)
+                    .ok_or(ComponentError::InvalidStateDimensions {
+                        expected: mapping.internal_edge_pos,
+                        actual: self.internal.edge_count(),
+                    })?;
+                let int_p = self.global_state_offset + p_local;
+                let int_h = self.global_state_offset + h_local;
 
                 if state.len() <= int_h || state.len() <= p_ext || state.len() <= h_ext {
                     return Err(ComponentError::InvalidStateDimensions {
@@ -363,8 +400,8 @@ impl Component for MacroComponent {
                     });
                 }
 
-                residuals[n_int_eqs + 2 * i] = state[p_ext] - state[int_p];
-                residuals[n_int_eqs + 2 * i + 1] = state[h_ext] - state[int_h];
+                residuals[n_int_res + 2 * i] = state[p_ext] - state[int_p];
+                residuals[n_int_res + 2 * i + 1] = state[h_ext] - state[int_h];
             }
         }
 
@@ -387,7 +424,7 @@ impl Component for MacroComponent {
             });
         }
 
-        let n_int_eqs = self.n_internal_equations();
+        let n_int_res = self.n_internal_residuals();
 
         // --- 1. Internal Jacobian entries ------------------------------------
         let internal_state: Vec<f64> = state[start..end].to_vec();
@@ -410,10 +447,15 @@ impl Component for MacroComponent {
         //   ∂r_h/∂state[h_ext] = +1
         //   ∂r_h/∂state[int_h] = −1
         for (i, mapping) in self.port_mappings.iter().enumerate() {
-            if let Some(&(p_ext, h_ext)) = self.external_edge_state_indices.get(i) {
-                let int_p = self.global_state_offset + 2 * mapping.internal_edge_pos;
-                let int_h = int_p + 1;
-                let row_p = n_int_eqs + 2 * i;
+            if let Some(&(_, p_ext, h_ext)) = self.external_edge_state_indices.get(i) {
+                let (p_local, h_local) =
+                    match self.internal_edge_ph_local(mapping.internal_edge_pos) {
+                        Some(v) => v,
+                        None => continue,
+                    };
+                let int_p = self.global_state_offset + p_local;
+                let int_h = self.global_state_offset + h_local;
+                let row_p = n_int_res + 2 * i;
                 let row_h = row_p + 1;
 
                 jacobian.add_entry(row_p, p_ext, 1.0);
@@ -427,8 +469,10 @@ impl Component for MacroComponent {
     }
 
     fn n_equations(&self) -> usize {
-        // Internal equations + 2 coupling equations per exposed port.
-        self.n_internal_equations() + 2 * self.port_mappings.len()
+        // Internal residuals (component eqs + mass-flow closures) + 2 coupling
+        // equations per exposed port (P and h only).
+        // TEMP (CM1.2): becomes `3 * n_ports` in CM1.3 when ṁ coupling is added.
+        self.n_internal_residuals() + 2 * self.port_mappings.len()
     }
 
     fn get_ports(&self) -> &[ConnectedPort] {
@@ -505,8 +549,12 @@ mod tests {
     }
 
     /// Build a simple linear subsystem: A → B → C (2 edges, 3 components).
+    ///
+    /// Intentionally not square (6 eqs / 5 unknowns) — used only to exercise
+    /// macro residual plumbing, not physical DoF. DoF gate disabled for that reason.
     fn build_simple_internal_system() -> System {
         let mut sys = System::new();
+        sys.set_enforce_dof_gate(false);
         let a = sys.add_component(make_mock(2));
         let b = sys.add_component(make_mock(2));
         let c = sys.add_component(make_mock(2));
@@ -521,11 +569,10 @@ mod tests {
         let sys = build_simple_internal_system();
         let mc = MacroComponent::new(sys);
 
-        // 3 components × 2 equations each = 6 equations (no ports exposed yet,
-        // so no coupling equations).
+        // 3 components × 2 equations each = 6 internal residuals (no ports exposed yet).
         assert_eq!(mc.n_equations(), 6);
-        // 2 edges → 4 state variables
-        assert_eq!(mc.internal_state_len(), 4);
+        // CM1.4: 2-edge linear chain → 1 branch → state_len = 1 + 2×2 = 5
+        assert_eq!(mc.internal_state_len(), 5);
         // No ports exposed yet
         assert!(mc.get_ports().is_empty());
     }
@@ -538,7 +585,7 @@ mod tests {
         let port = make_connected_port("R134a", 100_000.0, 400_000.0);
         mc.expose_port(0, "inlet", port.clone());
 
-        // 6 internal + 2 coupling = 8
+        // 6 internal residuals + 2 coupling = 8
         assert_eq!(mc.n_equations(), 8);
         assert_eq!(mc.get_ports().len(), 1);
         assert_eq!(mc.port_mappings()[0].name, "inlet");
@@ -556,7 +603,7 @@ mod tests {
         mc.expose_port(0, "inlet", port_in);
         mc.expose_port(1, "outlet", port_out);
 
-        // 6 internal + 4 coupling = 10
+        // 6 internal residuals + 4 coupling = 10
         assert_eq!(mc.n_equations(), 10);
         assert_eq!(mc.get_ports().len(), 2);
         assert_eq!(mc.port_mappings()[0].name, "inlet");
@@ -577,13 +624,13 @@ mod tests {
         let sys = build_simple_internal_system();
         let mc = MacroComponent::new(sys);
 
-        // 4 state variables for 2 internal edges (no external coupling)
-        let state = vec![1.0, 2.0, 3.0, 4.0];
+        // 6 state variables for 2 internal edges (ṁ,P,h each; no external coupling)
+        let state = vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0];
         let mut residuals = vec![0.0; mc.n_equations()];
 
         mc.compute_residuals(&state, &mut residuals).unwrap();
 
-        // 6 equations (no coupling ports)
+        // 6 internal residuals (3 components × 2 equations, no coupling ports)
         assert_eq!(residuals.len(), 6);
     }
 
@@ -593,8 +640,8 @@ mod tests {
         let mut mc = MacroComponent::new(sys);
         mc.set_global_state_offset(4);
 
-        // State vector: 4 padding + 4 internal = 8 total
-        let state = vec![0.0, 0.0, 0.0, 0.0, 1.0, 2.0, 3.0, 4.0];
+        // State vector: 4 padding + 6 internal = 10 total
+        let state = vec![0.0, 0.0, 0.0, 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0];
         let mut residuals = vec![0.0; mc.n_equations()];
 
         mc.compute_residuals(&state, &mut residuals).unwrap();
@@ -607,7 +654,7 @@ mod tests {
         let mut mc = MacroComponent::new(sys);
         mc.set_global_state_offset(4);
 
-        let state = vec![0.0; 5]; // Needs at least 8 (offset 4 + 4 internal vars)
+        let state = vec![0.0; 5]; // Needs at least 10 (offset 4 + 6 internal vars)
         let mut residuals = vec![0.0; mc.n_equations()];
 
         let result = mc.compute_residuals(&state, &mut residuals);
@@ -647,7 +694,7 @@ mod tests {
 
     #[test]
     fn test_coupling_residuals_and_jacobian() {
-        // 2-edge internal system: edge0 = (P0, h0), edge1 = (P1, h1)
+        // 2-edge internal system: edge0 = (ṁ0, P0, h0), edge1 = (ṁ1, P1, h1)
         let sys = build_simple_internal_system();
         let mut mc = MacroComponent::new(sys);
 
@@ -655,36 +702,28 @@ mod tests {
         let port = make_connected_port("R134a", 100_000.0, 400_000.0);
         mc.expose_port(0, "inlet", port);
 
-        // Simulate finalization: inject external edge state index (p=6, h=7)
-        // and global offset = 4 (4 parent edges before the macro's internal block).
+        // Concrete global layout with internal block at [4..10] (stride 3):
+        //   index 4 = ṁ_int_e0, 5 = P_int_e0, 6 = h_int_e0
+        //   index 7 = ṁ_int_e1, 8 = P_int_e1, 9 = h_int_e1
+        //   index 10 = ṁ_ext, 11 = P_ext, 12 = h_ext (external edge)
         mc.set_global_state_offset(4);
-        mc.set_system_context(4, &[(6, 7)]);
+        mc.set_system_context(4, &[(10, 11, 12)]);
 
-        // Global state: [*0, 1, 2, 3*, 4, 5, 6, 7, P_ext=1e5, h_ext=4e5]
-        //                                 ^--- internal block at [4..8]
-        //                                                  ^--- ext edge at (6,7)... wait,
-        // Let's use a concrete layout:
-        // indices 0..3: some other parent edges
-        // indices 4..7: internal block (2 edges * 2 vars)
-        //   4=P_int_e0, 5=h_int_e0, 6=P_int_e1, 7=h_int_e1
-        // indices 8,9: external edge (p_ext=8, h_ext=9)
-        mc.set_system_context(4, &[(8, 9)]);
+        let mut state = vec![0.0; 13];
+        state[5] = 1.5e5; // P_int_e0
+        state[6] = 3.9e5; // h_int_e0
+        state[11] = 2.0e5; // P_ext
+        state[12] = 4.1e5; // h_ext
 
-        let mut state = vec![0.0; 10];
-        state[4] = 1.5e5; // P_int_e0
-        state[5] = 3.9e5; // h_int_e0
-        state[8] = 2.0e5; // P_ext
-        state[9] = 4.1e5; // h_ext
-
-        let n_eqs = mc.n_equations(); // 6 internal + 2 coupling = 8
+        let n_eqs = mc.n_equations(); // 6 internal residuals + 2 coupling = 8
         assert_eq!(n_eqs, 8);
 
         let mut residuals = vec![0.0; n_eqs];
         mc.compute_residuals(&state, &mut residuals).unwrap();
 
-        // Coupling residuals:
-        //   r[6] = state[8] - state[4] = 2e5 - 1.5e5 =  0.5e5
-        //   r[7] = state[9] - state[5] = 4.1e5 - 3.9e5 = 0.2e5
+        // Coupling residuals occupy the last 2 rows (indices 6, 7):
+        //   r[6] = state[11] - state[5] = 2e5 - 1.5e5 = 0.5e5
+        //   r[7] = state[12] - state[6] = 4.1e5 - 3.9e5 = 0.2e5
         assert!(
             (residuals[6] - 0.5e5).abs() < 1.0,
             "r_P mismatch: {}",
@@ -701,28 +740,29 @@ mod tests {
         mc.jacobian_entries(&state, &mut jac).unwrap();
 
         let entries = jac.entries();
-        // Check that we have entries for p_ext=8 → +1, int_p=4 → -1
         let find = |row: usize, col: usize| -> Option<f64> {
             entries
                 .iter()
                 .find(|&&(r, c, _)| r == row && c == col)
                 .map(|&(_, _, v)| v)
         };
-        assert_eq!(find(6, 8), Some(1.0), "expect ∂r_P/∂p_ext = +1");
-        assert_eq!(find(6, 4), Some(-1.0), "expect ∂r_P/∂int_p = -1");
-        assert_eq!(find(7, 9), Some(1.0), "expect ∂r_h/∂h_ext = +1");
-        assert_eq!(find(7, 5), Some(-1.0), "expect ∂r_h/∂int_h = -1");
+        assert_eq!(find(6, 11), Some(1.0), "expect ∂r_P/∂p_ext = +1");
+        assert_eq!(find(6, 5), Some(-1.0), "expect ∂r_P/∂int_p = -1");
+        assert_eq!(find(7, 12), Some(1.0), "expect ∂r_h/∂h_ext = +1");
+        assert_eq!(find(7, 6), Some(-1.0), "expect ∂r_h/∂int_h = -1");
     }
 
     #[test]
     fn test_n_equations_empty_system() {
         let mut sys = System::new();
+        sys.set_enforce_dof_gate(false);
         let a = sys.add_component(make_mock(0));
         let b = sys.add_component(make_mock(0));
         sys.add_edge(a, b).unwrap();
         sys.finalize().unwrap();
 
         let mc = MacroComponent::new(sys);
+        // 2 components × 0 equations = 0 (no mass-flow closures since CM1.3).
         assert_eq!(mc.n_equations(), 0);
     }
 
@@ -744,6 +784,7 @@ mod tests {
 
         // Place it in a parent system alongside another component
         let mut parent = System::new();
+        parent.set_enforce_dof_gate(false);
         let mc_node = parent.add_component(Box::new(mc));
         let other = parent.add_component(make_mock(2));
         parent.add_edge(mc_node, other).unwrap();
@@ -765,14 +806,15 @@ mod tests {
         let port = make_connected_port("R134a", 1e5, 4e5);
         mc.expose_port(0, "inlet", port);
 
-        // Fake global state with known values in the internal block [0..4]
-        let global_state = vec![1.5e5, 3.9e5, 8.0e4, 4.2e5];
+        // CM1.4: internal block has 5 slots (1 ṁ_branch + 2×2 P,h for 2 edges)
+        // Layout: [0:ṁ_branch, 1:P_e0, 2:h_e0, 3:P_e1, 4:h_e1]
+        let global_state = vec![0.05, 1.5e5, 3.9e5, 8.0e4, 4.2e5];
         let snap = mc
             .to_snapshot(&global_state, Some("chiller_1".into()))
             .unwrap();
 
         assert_eq!(snap.label.as_deref(), Some("chiller_1"));
-        assert_eq!(snap.internal_edge_states.len(), 4);
+        assert_eq!(snap.internal_edge_states.len(), 5);
         assert_eq!(snap.port_names, vec!["inlet"]);
 
         // Round-trip through JSON

crates/solver/src/scaling.rs

@@ -0,0 +1,254 @@
+//! Jacobian row/column equilibration (Ruiz-style scaling).
+//!
+//! The Newton state vector mixes quantities of vastly different magnitudes:
+//! mass flow (ṁ ≈ 1 kg/s), pressure (P ≈ 1e5–4e6 Pa), enthalpy (h ≈ 1e5–5e5
+//! J/kg) and — once moist-air edges arrive — humidity ratio (W ≈ 0.005–0.025).
+//! Without rescaling, the Jacobian condition number can reach 1e10+, so the LU
+//! solve loses significant digits and Newton stalls or diverges.
+//!
+//! This module applies a **solution-preserving** diagonal scaling. Given the
+//! row/column factors `(d_r, d_c)` from [`equilibrate`], the scaled matrix
+//! `Ĵ = diag(d_r) · J · diag(d_c)` has row and column ∞-norms ≈ 1. The Newton
+//! system is then solved as
+//!
+//! ```text
+//!   Ĵ · y = diag(d_r) · (−r)
+//!   Δx    = diag(d_c) · y          (see `unscale_dx`)
+//! ```
+//!
+//! which is mathematically identical to the unscaled solve but numerically far
+//! more robust on stiff, mixed-unit, >50-variable systems.
+//!
+//! Row scaling uses the inverse row ∞-norm; column scaling the inverse column
+//! ∞-norm. Both are applied iteratively (Ruiz equilibration) until the norms are
+//! within tolerance of 1. This data-driven scheme subsumes the nominal
+//! per-variable column scaling described in the design doc (§3.5.2) without
+//! needing hand-tuned reference magnitudes.
+//!
+//! **Ref:** IDAES scaling theory,
+//! `idaes-pse.readthedocs.io/en/stable/explanations/scaling_toolbox/scaling_theory.html`
+//! and the Ruiz (2001) iterative equilibration algorithm. Design context:
+//! `_bmad-output/planning-artifacts/complex-machine-design.md#3.5`.
+
+use nalgebra::DMatrix;
+
+/// Maximum Ruiz sweeps before giving up on the unit-norm target.
+const MAX_ITERS: usize = 20;
+/// Convergence tolerance on the row/column ∞-norm deviation from 1.
+const TOL: f64 = 1e-3;
+
+/// Computes Ruiz row/column equilibration factors for a matrix.
+///
+/// Iteratively rescales rows and columns by the square root of their current
+/// ∞-norm until the norms are within tolerance of 1 (or [`MAX_ITERS`] is
+/// reached). Returns `(d_r, d_c)` such that `diag(d_r) · m · diag(d_c)` is
+/// equilibrated (row/column ∞-norms ≈ 1).
+///
+/// Zero rows/columns are left untouched (factor `1.0`), so the rank — and
+/// therefore singularity detection in a subsequent LU — is preserved exactly.
+///
+/// # Example
+///
+/// ```rust
+/// use entropyk_solver::scaling::equilibrate;
+/// use nalgebra::DMatrix;
+///
+/// // Badly scaled diagonal spanning 12 orders of magnitude.
+/// let m = DMatrix::from_row_slice(2, 2, &[1e6, 0.0, 0.0, 1e-6]);
+/// let (d_r, d_c) = equilibrate(&m);
+/// assert_eq!(d_r.len(), 2);
+/// assert_eq!(d_c.len(), 2);
+/// ```
+pub fn equilibrate(m: &DMatrix<f64>) -> (Vec<f64>, Vec<f64>) {
+    let n_rows = m.nrows();
+    let n_cols = m.ncols();
+    let mut d_r = vec![1.0_f64; n_rows];
+    let mut d_c = vec![1.0_f64; n_cols];
+
+    if n_rows == 0 || n_cols == 0 {
+        return (d_r, d_c);
+    }
+
+    for _ in 0..MAX_ITERS {
+        let mut max_deviation = 0.0_f64;
+
+        // Row scaling: divide each row factor by sqrt of its current ∞-norm.
+        for i in 0..n_rows {
+            let mut row_max = 0.0_f64;
+            for j in 0..n_cols {
+                let v = (d_r[i] * m[(i, j)] * d_c[j]).abs();
+                if v > row_max {
+                    row_max = v;
+                }
+            }
+            if row_max > 0.0 {
+                d_r[i] /= row_max.sqrt();
+                max_deviation = max_deviation.max((1.0 - row_max).abs());
+            }
+        }
+
+        // Column scaling: divide each column factor by sqrt of its current ∞-norm.
+        for j in 0..n_cols {
+            let mut col_max = 0.0_f64;
+            for i in 0..n_rows {
+                let v = (d_r[i] * m[(i, j)] * d_c[j]).abs();
+                if v > col_max {
+                    col_max = v;
+                }
+            }
+            if col_max > 0.0 {
+                d_c[j] /= col_max.sqrt();
+                max_deviation = max_deviation.max((1.0 - col_max).abs());
+            }
+        }
+
+        if max_deviation < TOL {
+            break;
+        }
+    }
+
+    (d_r, d_c)
+}
+
+/// Un-scales a solved step back to physical units.
+///
+/// After solving the equilibrated system `diag(d_r) · J · diag(d_c) · y =
+/// diag(d_r) · (−r)`, the true Newton step is recovered by reapplying the column
+/// scaling: `Δx[j] = d_c[j] · y[j]`.
+///
+/// If `y` and `d_c` differ in length, the shorter length governs the result
+/// (defensive — callers always pass matching lengths).
+///
+/// # Example
+///
+/// ```rust
+/// use entropyk_solver::scaling::unscale_dx;
+///
+/// let y = [2.0, -3.0];
+/// let d_c = [0.5, 10.0];
+/// let dx = unscale_dx(&y, &d_c);
+/// assert_eq!(dx, vec![1.0, -30.0]);
+/// ```
+pub fn unscale_dx(y: &[f64], d_c: &[f64]) -> Vec<f64> {
+    y.iter().zip(d_c.iter()).map(|(&yj, &cj)| cj * yj).collect()
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// Tests
+// ─────────────────────────────────────────────────────────────────────────────
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use approx::assert_relative_eq;
+
+    /// After equilibration, scaled row/column ∞-norms should be ≈ 1.
+    #[test]
+    fn test_equilibrate_unit_norms() {
+        let m = DMatrix::from_row_slice(2, 2, &[1e6, 1e3, 1e-3, 1e-6]);
+        let (d_r, d_c) = equilibrate(&m);
+
+        for i in 0..2 {
+            let row_max = (0..2)
+                .map(|j| (d_r[i] * m[(i, j)] * d_c[j]).abs())
+                .fold(0.0_f64, f64::max);
+            assert!(
+                (row_max - 1.0).abs() < 0.05,
+                "row {} ∞-norm not unit: {}",
+                i,
+                row_max
+            );
+        }
+        for j in 0..2 {
+            let col_max = (0..2)
+                .map(|i| (d_r[i] * m[(i, j)] * d_c[j]).abs())
+                .fold(0.0_f64, f64::max);
+            assert!(
+                (col_max - 1.0).abs() < 0.05,
+                "col {} ∞-norm not unit: {}",
+                j,
+                col_max
+            );
+        }
+    }
+
+    /// `unscale_dx` must be the exact inverse of column scaling: given a step
+    /// expressed in scaled coordinates `y = D_c^{-1} · x`, reapplying `d_c`
+    /// recovers `x` exactly.
+    #[test]
+    fn test_unscale_dx_inverts_column_scaling() {
+        let x = [3.0_f64, -7.0, 0.25];
+        let d_c = [0.5_f64, 10.0, 2.0];
+        // y = D_c^{-1} · x
+        let y: Vec<f64> = x.iter().zip(d_c.iter()).map(|(&xj, &cj)| xj / cj).collect();
+        let recovered = unscale_dx(&y, &d_c);
+        for (got, want) in recovered.iter().zip(x.iter()) {
+            assert_relative_eq!(got, want, epsilon = 1e-12);
+        }
+    }
+
+    /// Zero rows and columns must keep factor 1.0 so the rank is preserved and
+    /// singularity detection in the downstream LU is unchanged.
+    #[test]
+    fn test_zero_row_and_column_keep_unit_factor() {
+        // Row 1 is all zero; column 1 is all zero.
+        let m = DMatrix::from_row_slice(2, 2, &[5.0, 0.0, 0.0, 0.0]);
+        let (d_r, d_c) = equilibrate(&m);
+        assert_relative_eq!(d_r[1], 1.0, epsilon = 1e-12);
+        assert_relative_eq!(d_c[1], 1.0, epsilon = 1e-12);
+    }
+
+    /// Equilibration must slash the condition number of a badly-scaled matrix.
+    #[test]
+    fn test_equilibrate_reduces_condition_number() {
+        // Row 0 lives at ~1e6, row 1 at ~1; dynamic range keeps the small
+        // singular value well above underflow so κ is measured reliably.
+        let m = DMatrix::from_row_slice(2, 2, &[1.0e6, 2.0e6, 3.0, 1.0]);
+
+        let cond_before = condition_number(&m).unwrap();
+
+        let (d_r, d_c) = equilibrate(&m);
+        let mut scaled = m.clone();
+        for i in 0..2 {
+            for j in 0..2 {
+                scaled[(i, j)] *= d_r[i] * d_c[j];
+            }
+        }
+        let cond_after = condition_number(&scaled).unwrap();
+
+        assert!(
+            cond_after < cond_before / 100.0 && cond_after < 10.0,
+            "Ruiz should slash κ: before={:.3e}, after={:.3e}",
+            cond_before,
+            cond_after
+        );
+    }
+
+    /// Empty matrices return empty/unit factors without panicking.
+    #[test]
+    fn test_empty_matrix() {
+        let m: DMatrix<f64> = DMatrix::zeros(0, 0);
+        let (d_r, d_c) = equilibrate(&m);
+        assert!(d_r.is_empty());
+        assert!(d_c.is_empty());
+    }
+
+    /// Local κ helper (SVD σ_max/σ_min) for the test above only.
+    fn condition_number(m: &DMatrix<f64>) -> Option<f64> {
+        if m.nrows() == 0 || m.ncols() == 0 {
+            return None;
+        }
+        let svd = m.clone().svd(false, false);
+        let sv = svd.singular_values;
+        if sv.len() == 0 {
+            return None;
+        }
+        let sigma_max = sv.max();
+        let sigma_min = sv
+            .iter()
+            .filter(|&&s| s > 0.0)
+            .min_by(|a, b| a.partial_cmp(b).unwrap())
+            .copied();
+        sigma_min.map(|min| sigma_max / min)
+    }
+}

crates/solver/src/snapshot_params.rs

@@ -51,7 +51,7 @@ impl ParamsPlaceholder {
             "FloodedEvaporator" => 2,
             "Node" => 2,
             "Drum" => 8,
-            "ScrewEconomizerCompressor" => 5,
+            "ScrewEconomizerCompressor" | "ScrewCompressor" => 6,
             "RefrigerantSource" | "RefrigerantSink" => 2,
             "AirSource" | "AirSink" => 2,
             "BrineSource" | "BrineSink" => 2,
@@ -59,14 +59,22 @@ impl ParamsPlaceholder {
         }
     }
 
-    fn infer_ports(_type_name: &str) -> usize {
-        2 // Most components have 2 ports
+    fn infer_ports(type_name: &str) -> usize {
+        match type_name {
+            "ScrewEconomizerCompressor" | "ScrewCompressor" => 3,
+            _ => 2, // Most components have 2 ports
+        }
     }
 
     /// Returns the stored parameters.
     pub fn params(&self) -> &ComponentParams {
         &self.params
     }
+
+    /// Returns the inferred port count preserved for topology sizing.
+    pub fn n_ports(&self) -> usize {
+        self.n_ports
+    }
 }
 
 impl Component for ParamsPlaceholder {

crates/solver/src/solver.rs

@@ -80,6 +80,65 @@ pub enum SolverError {
         /// Energy balance error in W
         energy_error: f64,
     },
+
+    /// Solver failure with attached post-mortem convergence diagnostics.
+    ///
+    /// This wrapper keeps the existing concrete error variants unchanged, so
+    /// callers that do not need diagnostics can keep matching those variants.
+    #[error("{error}")]
+    WithDiagnostics {
+        /// Original solver error.
+        error: Box<SolverError>,
+        /// Diagnostics captured before the failure was returned.
+        diagnostics: Box<ConvergenceDiagnostics>,
+    },
+}
+
+impl SolverError {
+    /// Attach diagnostics to an error without changing its display text.
+    pub fn with_diagnostics(self, diagnostics: ConvergenceDiagnostics) -> Self {
+        if diagnostics.iteration_history.is_empty() {
+            return self;
+        }
+
+        let base_error = self.without_diagnostics();
+        Self::WithDiagnostics {
+            error: Box::new(base_error),
+            diagnostics: Box::new(diagnostics),
+        }
+    }
+
+    /// Attach diagnostics only when they are available.
+    pub fn with_optional_diagnostics(self, diagnostics: Option<ConvergenceDiagnostics>) -> Self {
+        match diagnostics {
+            Some(diagnostics) => self.with_diagnostics(diagnostics),
+            None => self,
+        }
+    }
+
+    /// Returns diagnostics captured for this failure, if any.
+    pub fn diagnostics(&self) -> Option<&ConvergenceDiagnostics> {
+        match self {
+            Self::WithDiagnostics { diagnostics, .. } => Some(diagnostics),
+            _ => None,
+        }
+    }
+
+    /// Returns the original concrete error variant, ignoring diagnostics.
+    pub fn base_error(&self) -> &SolverError {
+        match self {
+            Self::WithDiagnostics { error, .. } => error.base_error(),
+            _ => self,
+        }
+    }
+
+    /// Removes any diagnostics wrapper and returns the original error.
+    pub fn without_diagnostics(self) -> SolverError {
+        match self {
+            Self::WithDiagnostics { error, .. } => error.without_diagnostics(),
+            _ => self,
+        }
+    }
 }
 
 // ─────────────────────────────────────────────────────────────────────────────
@@ -509,7 +568,7 @@ impl VerboseConfig {
 ///
 /// Records the state of the solver at each iteration for debugging
 /// and post-mortem analysis.
-#[derive(Debug, Clone, Serialize, Deserialize)]
+#[derive(Debug, Clone, Default, PartialEq, Serialize, Deserialize)]
 pub struct IterationDiagnostics {
     /// Iteration number (0-indexed).
     pub iteration: usize,
@@ -532,6 +591,41 @@ pub struct IterationDiagnostics {
     ///
     /// Only populated when `log_jacobian_condition` is enabled.
     pub jacobian_condition: Option<f64>,
+
+    /// Index of the equation with the largest absolute residual this iteration.
+    ///
+    /// Inspired by IPM's NLPD toolchain: pinpoints which residual (and hence
+    /// which component / state slot) is dominating convergence so a slow or
+    /// stalled solve can be diagnosed. `None` if the residual vector is empty.
+    #[serde(default)]
+    pub max_residual_index: Option<usize>,
+
+    /// Magnitude of the largest absolute residual this iteration
+    /// ($\max_i |r_i|$, the $\ell_\infty$ residual norm).
+    #[serde(default)]
+    pub max_residual: f64,
+}
+
+/// Identifies the dominant (largest-magnitude) entry of a residual vector.
+///
+/// Returns the `(index, magnitude)` of the equation with the largest absolute
+/// residual — the $\ell_\infty$ norm together with its location. Returns
+/// `(None, 0.0)` for an empty vector. NaN entries are treated as dominant so
+/// they are surfaced rather than hidden.
+pub fn dominant_residual(residuals: &[f64]) -> (Option<usize>, f64) {
+    let mut index = None;
+    let mut max = 0.0_f64;
+    for (i, &r) in residuals.iter().enumerate() {
+        let magnitude = r.abs();
+        if index.is_none() || magnitude.is_nan() || magnitude > max {
+            index = Some(i);
+            max = magnitude;
+            if magnitude.is_nan() {
+                break;
+            }
+        }
+    }
+    (index, max)
 }
 
 /// Type of solver being used.
@@ -541,6 +635,8 @@ pub enum SolverType {
     NewtonRaphson,
     /// Sequential Substitution (Picard) solver.
     SequentialSubstitution,
+    /// Newton-homotopy continuation solver.
+    Homotopy,
 }
 
 impl std::fmt::Display for SolverType {
@@ -548,6 +644,7 @@ impl std::fmt::Display for SolverType {
         match self {
             SolverType::NewtonRaphson => write!(f, "Newton-Raphson"),
             SolverType::SequentialSubstitution => write!(f, "Sequential Substitution"),
+            SolverType::Homotopy => write!(f, "Newton-Homotopy"),
         }
     }
 }
@@ -579,7 +676,7 @@ impl std::fmt::Display for SwitchReason {
 /// Event record for solver switches in fallback strategy.
 ///
 /// Captures when and why the solver switched between strategies.
-#[derive(Debug, Clone, Serialize, Deserialize)]
+#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
 pub struct SolverSwitchEvent {
     /// Solver being switched from.
     pub from_solver: SolverType,
@@ -601,7 +698,7 @@ pub struct SolverSwitchEvent {
 ///
 /// Contains all diagnostic information collected during solving,
 /// suitable for JSON serialization and post-mortem analysis.
-#[derive(Debug, Clone, Default, Serialize, Deserialize)]
+#[derive(Debug, Clone, Default, PartialEq, Serialize, Deserialize)]
 pub struct ConvergenceDiagnostics {
     /// Total iterations performed.
     pub iterations: usize,
@@ -658,6 +755,18 @@ impl ConvergenceDiagnostics {
         self.solver_switches.push(event);
     }
 
+    /// Returns the location and magnitude of the dominant residual at the final
+    /// recorded iteration: `(equation_index, magnitude)`.
+    ///
+    /// This is the NLPD-style "bottleneck equation" — the residual that was
+    /// hardest to drive to zero when the solve stopped. Returns `None` if no
+    /// iterations were recorded.
+    pub fn final_dominant_residual(&self) -> Option<(usize, f64)> {
+        self.iteration_history
+            .last()
+            .and_then(|it| it.max_residual_index.map(|i| (i, it.max_residual)))
+    }
+
     /// Returns a human-readable summary of the diagnostics.
     pub fn summary(&self) -> String {
         let converged_str = if self.converged { "YES" } else { "NO" };
@@ -691,6 +800,10 @@ impl ConvergenceDiagnostics {
             summary.push_str(&format!("\nFinal Solver: {}", solver));
         }
 
+        if let Some((idx, mag)) = self.final_dominant_residual() {
+            summary.push_str(&format!("\nDominant Residual: eq[{}] = {:.3e}", idx, mag));
+        }
+
         summary
     }
 
@@ -726,8 +839,13 @@ pub(crate) fn apply_newton_step(
     for (i, s) in state.iter_mut().enumerate() {
         let proposed = *s + alpha * delta[i];
         *s = match &clipping_mask[i] {
-            Some((min, max)) => proposed.clamp(*min, *max),
-            None => proposed,
+            // `f64::clamp` panics if min > max or either bound is NaN. Guard the
+            // bounds explicitly to honour the zero-panic policy: an invalid bound
+            // pair is treated as "no constraint" rather than crashing the solver.
+            Some((min, max)) if min.is_finite() && max.is_finite() && min <= max => {
+                proposed.clamp(*min, *max)
+            }
+            _ => proposed,
         };
     }
 }
@@ -748,6 +866,25 @@ mod tests {
         accepts_dyn_solver(&mut newton);
     }
 
+    /// Inverted or NaN bound pairs must not panic `apply_newton_step` (the
+    /// `f64::clamp` zero-panic guard); they are treated as "no constraint".
+    #[test]
+    fn test_apply_newton_step_invalid_bounds_do_not_panic() {
+        let mut state = vec![0.0, 0.0, 0.0];
+        let delta = vec![1.0, 1.0, 1.0];
+        let mask = vec![
+            Some((10.0, -10.0)),   // inverted: min > max
+            Some((f64::NAN, 1.0)), // NaN bound
+            Some((-5.0, 5.0)),     // valid: should clamp normally
+        ];
+        apply_newton_step(&mut state, &delta, &mask, 1.0);
+        // Invalid bounds → unconstrained update applied.
+        assert_eq!(state[0], 1.0);
+        assert_eq!(state[1], 1.0);
+        // Valid bounds → normal clamp (1.0 is within [-5, 5]).
+        assert_eq!(state[2], 1.0);
+    }
+
     /// Verify that `Box<dyn Solver>` can be constructed from concrete types.
     #[test]
     fn test_box_dyn_solver_compiles() {
@@ -830,4 +967,93 @@ mod tests {
         assert_eq!(state[0], 0.0, "Bounded variable should be clipped");
         assert_eq!(state[1], 60.0, "Unbounded variable should NOT be clipped");
     }
+
+    // ── NLPD-style dominant-residual diagnostics ──────────────────────────────
+
+    #[test]
+    fn test_dominant_residual_picks_largest_magnitude() {
+        let residuals = vec![0.1, -5.0, 2.0, -0.3];
+        let (idx, mag) = dominant_residual(&residuals);
+        assert_eq!(idx, Some(1), "index 1 has the largest |residual|");
+        assert!((mag - 5.0).abs() < 1e-15);
+    }
+
+    #[test]
+    fn test_dominant_residual_empty_vector() {
+        let (idx, mag) = dominant_residual(&[]);
+        assert_eq!(idx, None);
+        assert_eq!(mag, 0.0);
+    }
+
+    #[test]
+    fn test_dominant_residual_surfaces_nan() {
+        let residuals = vec![1.0, f64::NAN, 0.5];
+        let (idx, _mag) = dominant_residual(&residuals);
+        assert_eq!(idx, Some(1), "NaN residual must be surfaced, not hidden");
+    }
+
+    #[test]
+    fn test_final_dominant_residual_reports_last_iteration() {
+        let mut diag = ConvergenceDiagnostics::new();
+        assert_eq!(diag.final_dominant_residual(), None);
+
+        diag.push_iteration(IterationDiagnostics {
+            iteration: 0,
+            residual_norm: 10.0,
+            max_residual_index: Some(2),
+            max_residual: 8.0,
+            ..Default::default()
+        });
+        diag.push_iteration(IterationDiagnostics {
+            iteration: 1,
+            residual_norm: 1.0,
+            max_residual_index: Some(5),
+            max_residual: 0.9,
+            ..Default::default()
+        });
+
+        assert_eq!(diag.final_dominant_residual(), Some((5, 0.9)));
+    }
+
+    #[test]
+    fn test_diagnostics_dominant_residual_serializes_to_json() {
+        let mut diag = ConvergenceDiagnostics::new();
+        diag.push_iteration(IterationDiagnostics {
+            iteration: 0,
+            max_residual_index: Some(3),
+            max_residual: 4.2,
+            ..Default::default()
+        });
+        let json = diag.dump_diagnostics(VerboseOutputFormat::Json);
+        assert!(json.contains("max_residual_index"));
+        assert!(json.contains("max_residual"));
+    }
+
+    #[test]
+    fn test_solver_error_exposes_failure_diagnostics() {
+        let mut diag = ConvergenceDiagnostics::new();
+        diag.iterations = 2;
+        diag.final_residual = 12.0;
+        diag.push_iteration(IterationDiagnostics {
+            iteration: 2,
+            residual_norm: 12.0,
+            max_residual_index: Some(7),
+            max_residual: 11.5,
+            ..Default::default()
+        });
+
+        let err = SolverError::NonConvergence {
+            iterations: 2,
+            final_residual: 12.0,
+        }
+        .with_diagnostics(diag);
+
+        assert!(matches!(
+            err.base_error(),
+            SolverError::NonConvergence { .. }
+        ));
+        let diagnostics = err.diagnostics().expect("diagnostics should be attached");
+        assert_eq!(diagnostics.final_residual, 12.0);
+        assert_eq!(diagnostics.final_dominant_residual(), Some((7, 11.5)));
+    }
 }

crates/solver/src/strategies/fallback.rs

@@ -31,7 +31,7 @@ use crate::solver::{
 };
 use crate::system::System;
 
-use super::{NewtonConfig, PicardConfig};
+use super::{HomotopyConfig, NewtonConfig, PicardConfig};
 
 /// Configuration for the intelligent fallback solver.
 ///
@@ -143,7 +143,7 @@ impl FallbackState {
             self.best_residual = Some(residual);
         }
     }
-    
+
     /// Record a solver switch event (Story 7.4)
     fn record_switch(
         &mut self,
@@ -188,15 +188,29 @@ pub struct FallbackSolver {
     pub newton_config: NewtonConfig,
     /// Sequential Substitution (Picard) configuration.
     pub picard_config: PicardConfig,
+    /// Optional Newton-homotopy continuation used as a last-resort recovery.
+    ///
+    /// When set, the homotopy solver is invoked if both Newton and Picard fail
+    /// (divergence or non-convergence) from the cold start. This mirrors IPM
+    /// BOLT's escalating initialization cascade (`GLBL` → `iGLBL` → `PRVS` →
+    /// `iPRVS`): cheap solvers first, robust continuation only when needed.
+    /// `None` (default) preserves the original Newton↔Picard-only behaviour.
+    pub homotopy_config: Option<HomotopyConfig>,
 }
 
 impl FallbackSolver {
     /// Creates a new fallback solver with the given configuration.
+    ///
+    /// The Picard fallback stage is configured with Anderson acceleration
+    /// (depth 3) by default, which converges the fixed-point iteration
+    /// super-linearly and typically halves the fallback iteration count without
+    /// affecting robustness. Override via [`Self::with_picard_config`].
     pub fn new(config: FallbackConfig) -> Self {
         Self {
             config,
             newton_config: NewtonConfig::default(),
-            picard_config: PicardConfig::default(),
+            picard_config: PicardConfig::default().with_anderson(3),
+            homotopy_config: None,
         }
     }
 
@@ -217,13 +231,27 @@ impl FallbackSolver {
         self
     }
 
+    /// Enables Newton-homotopy continuation as a last-resort recovery stage.
+    ///
+    /// If both Newton and Picard fail from the cold start, the homotopy solver
+    /// is invoked to walk in from `λ = 0` to `λ = 1` (see [`HomotopyConfig`]).
+    /// When the supplied config has no `initial_state`, the fallback solver's
+    /// Newton initial state is reused so all stages share the same cold start.
+    pub fn with_homotopy(mut self, config: HomotopyConfig) -> Self {
+        self.homotopy_config = Some(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.
+    /// Delegates to `newton_config`, `picard_config`, and the optional homotopy
+    /// stage so the initial state is used regardless of which solver runs.
     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.picard_config.initial_state = Some(state.clone());
+        if let Some(ref mut h) = self.homotopy_config {
+            h.initial_state = Some(state);
+        }
         self
     }
 
@@ -251,10 +279,10 @@ impl FallbackSolver {
         timeout: Option<Duration>,
     ) -> Result<ConvergedState, SolverError> {
         let mut state = FallbackState::new();
-        
+
         // Verbose mode setup
-        let verbose_enabled = self.config.verbose_config.enabled 
-            && self.config.verbose_config.is_any_enabled();
+        let verbose_enabled =
+            self.config.verbose_config.enabled && self.config.verbose_config.is_any_enabled();
         let mut diagnostics = if verbose_enabled {
             Some(ConvergenceDiagnostics::with_capacity(100))
         } else {
@@ -264,7 +292,7 @@ impl FallbackSolver {
         // Pre-configure solver configs once
         let mut newton_cfg = self.newton_config.clone();
         let mut picard_cfg = self.picard_config.clone();
-        
+
         // Propagate verbose config to child solvers
         newton_cfg.verbose_config = self.config.verbose_config.clone();
         picard_cfg.verbose_config = self.config.verbose_config.clone();
@@ -301,17 +329,18 @@ impl FallbackSolver {
                     if let Some(ref mut diag) = diagnostics {
                         diag.iterations = state.total_iterations;
                         diag.final_residual = converged.final_residual;
-                        diag.best_residual = state.best_residual.unwrap_or(converged.final_residual);
+                        diag.best_residual =
+                            state.best_residual.unwrap_or(converged.final_residual);
                         diag.converged = true;
                         diag.timing_ms = start_time.elapsed().as_millis() as u64;
                         diag.final_solver = Some(state.current_solver.into());
                         diag.solver_switches = state.switch_events.clone();
-                        
+
                         // Merge iteration history from child solver if available
                         if let Some(ref child_diag) = converged.diagnostics {
                             diag.iteration_history = child_diag.iteration_history.clone();
                         }
-                        
+
                         if self.config.verbose_config.log_residuals {
                             tracing::info!("{}", diag.summary());
                         }
@@ -327,287 +356,385 @@ impl FallbackSolver {
                         switch_count = state.switch_count,
                         "Fallback solver converged"
                     );
-                    
+
                     // Return with diagnostics if verbose mode enabled
                     return Ok(if let Some(d) = diagnostics {
-                        ConvergedState { diagnostics: Some(d), ..converged }
-                    } else { converged });
+                        ConvergedState {
+                            diagnostics: Some(d),
+                            ..converged
+                        }
+                    } else {
+                        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,
-                            state.total_iterations,
-                            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 => {
-                            // Get residual from error context (use best known)
-                            let residual_at_switch = state.best_residual.unwrap_or(f64::MAX);
-                            
-                            // 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;
-                                let prev_solver = state.current_solver;
-                                state.current_solver = CurrentSolver::Picard;
-                                
-                                // Record switch event
-                                state.record_switch(
-                                    prev_solver,
-                                    state.current_solver,
-                                    SwitchReason::Divergence,
-                                    residual_at_switch,
-                                );
-                                
-                                // Verbose logging
-                                if verbose_enabled && self.config.verbose_config.log_solver_switches {
-                                    tracing::info!(
-                                        from = "NewtonRaphson",
-                                        to = "Picard",
-                                        reason = "divergence",
-                                        switch_count = state.switch_count,
-                                        residual = residual_at_switch,
-                                        "Solver switch (max switches reached)"
-                                    );
-                                }
-                                
+                Err(err) => {
+                    let child_diagnostics = err.diagnostics().cloned();
+                    match err.without_diagnostics() {
+                        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,
+                                    state.total_iterations,
+                                    best_residual,
+                                    ConvergenceStatus::TimedOutWithBestState,
+                                    SimulationMetadata::new(system.input_hash()),
+                                ));
+                            }
+                            return Err(SolverError::Timeout { timeout_ms }
+                                .with_optional_diagnostics(child_diagnostics));
+                        }
+                        SolverError::Divergence { 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 Err(SolverError::Divergence { reason }
+                                    .with_optional_diagnostics(child_diagnostics));
+                            }
+
+                            match state.current_solver {
+                                CurrentSolver::Newton => {
+                                    // Get residual from error context (use best known)
+                                    let residual_at_switch =
+                                        state.best_residual.unwrap_or(f64::MAX);
+
+                                    // 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;
+                                        let prev_solver = state.current_solver;
+                                        state.current_solver = CurrentSolver::Picard;
+
+                                        // Record switch event
+                                        state.record_switch(
+                                            prev_solver,
+                                            state.current_solver,
+                                            SwitchReason::Divergence,
+                                            residual_at_switch,
+                                        );
+
+                                        // Verbose logging
+                                        if verbose_enabled
+                                            && self.config.verbose_config.log_solver_switches
+                                        {
+                                            tracing::info!(
+                                                from = "NewtonRaphson",
+                                                to = "Picard",
+                                                reason = "divergence",
+                                                switch_count = state.switch_count,
+                                                residual = residual_at_switch,
+                                                "Solver switch (max switches reached)"
+                                            );
+                                        }
+
+                                        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;
-                                let prev_solver = state.current_solver;
-                                state.current_solver = CurrentSolver::Picard;
-                                
-                                // Record switch event
-                                state.record_switch(
-                                    prev_solver,
-                                    state.current_solver,
-                                    SwitchReason::Divergence,
-                                    residual_at_switch,
-                                );
-                                
-                                // Verbose logging
-                                if verbose_enabled && self.config.verbose_config.log_solver_switches {
-                                    tracing::info!(
-                                        from = "NewtonRaphson",
-                                        to = "Picard",
-                                        reason = "divergence",
-                                        switch_count = state.switch_count,
-                                        residual = residual_at_switch,
-                                        "Solver switch"
-                                    );
+                                    } else {
+                                        // Switch to Picard
+                                        state.switch_count += 1;
+                                        let prev_solver = state.current_solver;
+                                        state.current_solver = CurrentSolver::Picard;
+
+                                        // Record switch event
+                                        state.record_switch(
+                                            prev_solver,
+                                            state.current_solver,
+                                            SwitchReason::Divergence,
+                                            residual_at_switch,
+                                        );
+
+                                        // Verbose logging
+                                        if verbose_enabled
+                                            && self.config.verbose_config.log_solver_switches
+                                        {
+                                            tracing::info!(
+                                                from = "NewtonRaphson",
+                                                to = "Picard",
+                                                reason = "divergence",
+                                                switch_count = state.switch_count,
+                                                residual = residual_at_switch,
+                                                "Solver switch"
+                                            );
+                                        }
+
+                                        tracing::warn!(
+                                            switch_count = state.switch_count,
+                                            reason = reason,
+                                            "Newton diverged, switching to Picard"
+                                        );
+                                    }
+                                    // Continue loop with 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!(
+                                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;
+                                        // 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 Err(SolverError::Divergence { reason }
+                                            .with_optional_diagnostics(child_diagnostics));
+                                    }
+                                }
                             }
                         }
-                    }
-                }
-                Err(SolverError::NonConvergence {
-                    iterations,
-                    final_residual,
-                }) => {
-                    state.total_iterations += iterations;
-                    
-                    // Non-convergence: check if we should try the other solver
-                    if !self.config.fallback_enabled {
-                        return Err(SolverError::NonConvergence {
+                        SolverError::NonConvergence {
                             iterations,
                             final_residual,
-                        });
-                    }
+                        } => {
+                            state.total_iterations += iterations;
 
-                    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;
-                                let prev_solver = state.current_solver;
-                                state.current_solver = CurrentSolver::Picard;
-                                
-                                // Record switch event
-                                state.record_switch(
-                                    prev_solver,
-                                    state.current_solver,
-                                    SwitchReason::SlowConvergence,
+                            // Non-convergence: check if we should try the other solver
+                            if !self.config.fallback_enabled {
+                                return Err(SolverError::NonConvergence {
+                                    iterations,
                                     final_residual,
-                                );
-                                
-                                // Verbose logging
-                                if verbose_enabled && self.config.verbose_config.log_solver_switches {
-                                    tracing::info!(
-                                        from = "NewtonRaphson",
-                                        to = "Picard",
-                                        reason = "slow_convergence",
-                                        switch_count = state.switch_count,
-                                        residual = final_residual,
-                                        "Solver switch (max switches reached)"
-                                    );
                                 }
-                                
-                                tracing::info!(
+                                .with_optional_diagnostics(child_diagnostics));
+                            }
+
+                            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;
+                                        let prev_solver = state.current_solver;
+                                        state.current_solver = CurrentSolver::Picard;
+
+                                        // Record switch event
+                                        state.record_switch(
+                                            prev_solver,
+                                            state.current_solver,
+                                            SwitchReason::SlowConvergence,
+                                            final_residual,
+                                        );
+
+                                        // Verbose logging
+                                        if verbose_enabled
+                                            && self.config.verbose_config.log_solver_switches
+                                        {
+                                            tracing::info!(
+                                                from = "NewtonRaphson",
+                                                to = "Picard",
+                                                reason = "slow_convergence",
+                                                switch_count = state.switch_count,
+                                                residual = final_residual,
+                                                "Solver switch (max switches reached)"
+                                            );
+                                        }
+
+                                        tracing::info!(
                                     switch_count = state.switch_count,
                                     "Max switches reached, committing to Picard permanently"
                                 );
-                            } else {
-                                state.switch_count += 1;
-                                let prev_solver = state.current_solver;
-                                state.current_solver = CurrentSolver::Picard;
-                                
-                                // Record switch event
-                                state.record_switch(
-                                    prev_solver,
-                                    state.current_solver,
-                                    SwitchReason::SlowConvergence,
-                                    final_residual,
-                                );
-                                
-                                // Verbose logging
-                                if verbose_enabled && self.config.verbose_config.log_solver_switches {
-                                    tracing::info!(
-                                        from = "NewtonRaphson",
-                                        to = "Picard",
-                                        reason = "slow_convergence",
-                                        switch_count = state.switch_count,
-                                        residual = final_residual,
-                                        "Solver switch"
-                                    );
-                                }
-                                
-                                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,
-                                });
-                            }
+                                    } else {
+                                        state.switch_count += 1;
+                                        let prev_solver = state.current_solver;
+                                        state.current_solver = CurrentSolver::Picard;
 
-                            // Check if residual is low enough to try Newton
-                            if final_residual < self.config.return_to_newton_threshold {
-                                state.switch_count += 1;
-                                let prev_solver = state.current_solver;
-                                state.current_solver = CurrentSolver::Newton;
-                                
-                                // Record switch event
-                                state.record_switch(
-                                    prev_solver,
-                                    state.current_solver,
-                                    SwitchReason::ReturnToNewton,
-                                    final_residual,
-                                );
-                                
-                                // Verbose logging
-                                if verbose_enabled && self.config.verbose_config.log_solver_switches {
-                                    tracing::info!(
-                                        from = "Picard",
-                                        to = "NewtonRaphson",
-                                        reason = "return_to_newton",
-                                        switch_count = state.switch_count,
-                                        residual = final_residual,
-                                        threshold = self.config.return_to_newton_threshold,
-                                        "Solver switch (Picard stabilized)"
-                                    );
+                                        // Record switch event
+                                        state.record_switch(
+                                            prev_solver,
+                                            state.current_solver,
+                                            SwitchReason::SlowConvergence,
+                                            final_residual,
+                                        );
+
+                                        // Verbose logging
+                                        if verbose_enabled
+                                            && self.config.verbose_config.log_solver_switches
+                                        {
+                                            tracing::info!(
+                                                from = "NewtonRaphson",
+                                                to = "Picard",
+                                                reason = "slow_convergence",
+                                                switch_count = state.switch_count,
+                                                residual = final_residual,
+                                                "Solver switch"
+                                            );
+                                        }
+
+                                        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,
+                                        }
+                                        .with_optional_diagnostics(child_diagnostics));
+                                    }
+
+                                    // Check if residual is low enough to try Newton
+                                    if final_residual < self.config.return_to_newton_threshold {
+                                        state.switch_count += 1;
+                                        let prev_solver = state.current_solver;
+                                        state.current_solver = CurrentSolver::Newton;
+
+                                        // Record switch event
+                                        state.record_switch(
+                                            prev_solver,
+                                            state.current_solver,
+                                            SwitchReason::ReturnToNewton,
+                                            final_residual,
+                                        );
+
+                                        // Verbose logging
+                                        if verbose_enabled
+                                            && self.config.verbose_config.log_solver_switches
+                                        {
+                                            tracing::info!(
+                                                from = "Picard",
+                                                to = "NewtonRaphson",
+                                                reason = "return_to_newton",
+                                                switch_count = state.switch_count,
+                                                residual = final_residual,
+                                                threshold = self.config.return_to_newton_threshold,
+                                                "Solver switch (Picard stabilized)"
+                                            );
+                                        }
+
+                                        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,
+                                        }
+                                        .with_optional_diagnostics(child_diagnostics));
+                                    }
                                 }
-                                
-                                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,
-                                });
                             }
                         }
+                        other => {
+                            // InvalidSystem or other errors - propagate immediately
+                            return Err(other.with_optional_diagnostics(child_diagnostics));
+                        }
                     }
                 }
-                Err(other) => {
-                    // InvalidSystem or other errors - propagate immediately
-                    return Err(other);
-                }
+            }
+        }
+    }
+
+    /// Attempts Newton-homotopy continuation as a last-resort recovery after the
+    /// primary Newton/Picard stages have failed.
+    ///
+    /// Returns the homotopy result on success; otherwise returns the *primary*
+    /// error (the homotopy failure is logged but not surfaced, since the primary
+    /// error is the more actionable diagnostic). Structural `InvalidSystem`
+    /// errors are never retried — they indicate a malformed model, not a hard
+    /// cold start.
+    fn try_homotopy_recovery(
+        &self,
+        system: &mut System,
+        primary_err: SolverError,
+        remaining: Option<Duration>,
+    ) -> Result<ConvergedState, SolverError> {
+        if matches!(primary_err.base_error(), SolverError::InvalidSystem { .. }) {
+            return Err(primary_err);
+        }
+
+        let Some(mut homotopy) = self.homotopy_config.clone() else {
+            return Err(primary_err);
+        };
+
+        // Share the cold start with the primary solvers unless explicitly set.
+        if homotopy.initial_state.is_none() {
+            homotopy.initial_state = self.newton_config.initial_state.clone();
+        }
+        // Inherit the remaining global time budget if the stage has none.
+        if homotopy.timeout.is_none() {
+            homotopy.timeout = remaining;
+        }
+
+        tracing::info!(
+            error = %primary_err,
+            "Primary solvers failed; attempting Newton-homotopy continuation as last resort"
+        );
+
+        match homotopy.solve(system) {
+            Ok(converged) => {
+                tracing::info!(
+                    iterations = converged.iterations,
+                    final_residual = converged.final_residual,
+                    "Homotopy continuation recovered convergence"
+                );
+                Ok(converged)
+            }
+            Err(homotopy_err) => {
+                let primary_diagnostics = primary_err.diagnostics().cloned();
+                let homotopy_diagnostics = homotopy_err.diagnostics().cloned();
+                let diagnostics = match (primary_diagnostics, homotopy_diagnostics) {
+                    (Some(primary), Some(homotopy)) => {
+                        if homotopy.iterations >= primary.iterations {
+                            Some(homotopy)
+                        } else {
+                            Some(primary)
+                        }
+                    }
+                    (Some(primary), None) => Some(primary),
+                    (None, Some(homotopy)) => Some(homotopy),
+                    (None, None) => None,
+                };
+                tracing::warn!(
+                    error = %homotopy_err,
+                    "Homotopy continuation also failed; returning primary error"
+                );
+                Err(primary_err
+                    .without_diagnostics()
+                    .with_optional_diagnostics(diagnostics))
             }
         }
     }
@@ -622,20 +749,32 @@ impl Solver for FallbackSolver {
             fallback_enabled = self.config.fallback_enabled,
             return_to_newton_threshold = self.config.return_to_newton_threshold,
             max_fallback_switches = self.config.max_fallback_switches,
+            homotopy_recovery = self.homotopy_config.is_some(),
             "Fallback solver starting"
         );
 
-        if self.config.fallback_enabled {
+        let primary = if self.config.fallback_enabled {
             self.solve_with_fallback(system, start_time, timeout)
         } else {
             // Fallback disabled - run pure Newton
             self.newton_config.solve(system)
+        };
+
+        match primary {
+            Ok(converged) => Ok(converged),
+            Err(primary_err) => {
+                let remaining = timeout.map(|t| t.saturating_sub(start_time.elapsed()));
+                self.try_homotopy_recovery(system, primary_err, remaining)
+            }
         }
     }
 
     fn with_timeout(mut self, timeout: Duration) -> Self {
         self.newton_config.timeout = Some(timeout);
         self.picard_config.timeout = Some(timeout);
+        if let Some(ref mut h) = self.homotopy_config {
+            h.timeout = Some(timeout);
+        }
         self
     }
 }
@@ -684,4 +823,60 @@ mod tests {
         system.finalize().unwrap();
         assert!(boxed.solve(&mut system).is_err());
     }
+
+    // ── Homotopy last-resort recovery wiring ──────────────────────────────────
+
+    #[test]
+    fn test_fallback_homotopy_disabled_by_default() {
+        let solver = FallbackSolver::default_solver();
+        assert!(solver.homotopy_config.is_none());
+    }
+
+    #[test]
+    fn test_fallback_with_homotopy_sets_config() {
+        let solver = FallbackSolver::default_solver()
+            .with_homotopy(HomotopyConfig::default().with_initial_steps(20));
+        let h = solver
+            .homotopy_config
+            .expect("homotopy should be configured");
+        assert_eq!(h.initial_steps, 20);
+    }
+
+    #[test]
+    fn test_with_initial_state_propagates_to_homotopy() {
+        let solver = FallbackSolver::default_solver()
+            .with_homotopy(HomotopyConfig::default())
+            .with_initial_state(vec![1.0, 2.0, 3.0]);
+        assert_eq!(
+            solver.newton_config.initial_state,
+            Some(vec![1.0, 2.0, 3.0])
+        );
+        assert_eq!(
+            solver.picard_config.initial_state,
+            Some(vec![1.0, 2.0, 3.0])
+        );
+        assert_eq!(
+            solver.homotopy_config.unwrap().initial_state,
+            Some(vec![1.0, 2.0, 3.0])
+        );
+    }
+
+    #[test]
+    fn test_with_timeout_propagates_to_homotopy() {
+        let timeout = Duration::from_millis(750);
+        let solver = FallbackSolver::default_solver()
+            .with_homotopy(HomotopyConfig::default())
+            .with_timeout(timeout);
+        assert_eq!(solver.homotopy_config.unwrap().timeout, Some(timeout));
+    }
+
+    #[test]
+    fn test_invalid_system_not_retried_by_homotopy() {
+        // An empty (degenerate) system yields InvalidSystem; the homotopy stage
+        // must NOT retry it — a malformed model is not a hard cold start.
+        let mut solver = FallbackSolver::default_solver().with_homotopy(HomotopyConfig::default());
+        let mut system = System::new();
+        system.finalize().unwrap();
+        assert!(solver.solve(&mut system).is_err());
+    }
 }

crates/solver/src/strategies/homotopy.rs

@@ -0,0 +1,496 @@
+//! Newton-homotopy continuation solver for robust cold starts.
+//!
+//! This module provides [`HomotopyConfig`], a globally-convergent continuation
+//! solver that improves on a naive cold start without requiring a database of
+//! previous solutions (the IPM BOLT `GLBL`/`iPRVS` approach) or any manual
+//! tuning.
+//!
+//! # The Newton homotopy
+//!
+//! Given the target system `F(x) = 0` and an arbitrary initial guess `x₀`, define
+//! the homotopy
+//!
+//! ```text
+//! H(x, λ) = F(x) − (1 − λ) · F(x₀),   λ ∈ [0, 1]
+//! ```
+//!
+//! At `λ = 0`, `H(x₀, 0) = F(x₀) − F(x₀) = 0`, so the initial guess is an
+//! **exact** solution of the deformed system. At `λ = 1`, `H(x, 1) = F(x)`, the
+//! real system. The solver walks `λ` from 0 to 1, solving `H(·, λ) = 0` with an
+//! inner Newton iteration at each step and using the previous converged point as
+//! the next initial guess.
+//!
+//! # Why it reuses the analytic Jacobian unchanged
+//!
+//! The subtracted term `(1 − λ)·F(x₀)` is **constant in `x`**, so
+//!
+//! ```text
+//! ∂H/∂x = ∂F/∂x = J(x)
+//! ```
+//!
+//! The inner Newton step therefore uses the exact, component-wise analytic
+//! Jacobian assembled by [`System::assemble_jacobian`] — no finite differences
+//! and no changes to any component are required. This keeps Entropyk's
+//! structural advantage over finite-difference solvers (IPM eKINSOL) while adding
+//! cold-start robustness. A `use_numerical_jacobian` flag is provided for parity
+//! with [`crate::strategies::NewtonConfig`] when a component's analytic Jacobian
+//! is unavailable.
+//!
+//! # Adaptive step control
+//!
+//! The `λ` increment starts at `1 / initial_steps` and is **halved** whenever an
+//! inner Newton solve fails to converge, retrying from the last good `λ`. On
+//! success the increment is gently grown again. This predictor–corrector scheme
+//! automatically takes small steps through difficult regions (phase boundaries,
+//! stiff correlations) and large steps through easy ones.
+
+use std::time::{Duration, Instant};
+
+use crate::jacobian::JacobianMatrix;
+use crate::metadata::SimulationMetadata;
+use crate::solver::{
+    apply_newton_step, dominant_residual, ConvergedState, ConvergenceDiagnostics,
+    ConvergenceStatus, IterationDiagnostics, Solver, SolverError, SolverType,
+};
+use crate::system::System;
+use entropyk_components::JacobianBuilder;
+
+/// Configuration for the Newton-homotopy continuation solver.
+///
+/// Solves `F(x) = 0` by continuation on the homotopy
+/// `H(x, λ) = F(x) − (1 − λ)·F(x₀)` from `λ = 0` (where `x₀` is exact) to
+/// `λ = 1` (the real system).
+#[derive(Debug, Clone, PartialEq)]
+pub struct HomotopyConfig {
+    /// Initial number of `λ` subdivisions. The starting step is `1 / initial_steps`.
+    /// Default: 10.
+    pub initial_steps: usize,
+    /// Maximum Newton iterations allowed for each inner `λ` solve. Default: 50.
+    pub inner_max_iterations: usize,
+    /// Convergence tolerance (L2 residual norm) for each inner `λ` solve.
+    /// Default: 1e-8.
+    pub inner_tolerance: f64,
+    /// Final convergence tolerance (L2 residual norm) checked at `λ = 1`.
+    /// Default: 1e-6.
+    pub tolerance: f64,
+    /// Smallest allowed `λ` increment before the solver gives up. Default: 1e-4.
+    pub min_lambda_step: f64,
+    /// Divergence guard: inner residual norm above this aborts the inner solve.
+    /// Default: 1e12.
+    pub divergence_threshold: f64,
+    /// Use a finite-difference Jacobian instead of the analytic one. Default: false.
+    pub use_numerical_jacobian: bool,
+    /// Relative step for the finite-difference Jacobian. Default: 1e-5.
+    pub numerical_epsilon: f64,
+    /// Optional overall time budget.
+    pub timeout: Option<Duration>,
+    /// Initial guess `x₀`. When `None`, a zero vector is used.
+    pub initial_state: Option<Vec<f64>>,
+}
+
+impl Default for HomotopyConfig {
+    fn default() -> Self {
+        Self {
+            initial_steps: 10,
+            inner_max_iterations: 50,
+            inner_tolerance: 1e-8,
+            tolerance: 1e-6,
+            min_lambda_step: 1e-4,
+            divergence_threshold: 1e12,
+            use_numerical_jacobian: false,
+            numerical_epsilon: 1e-5,
+            timeout: None,
+            initial_state: None,
+        }
+    }
+}
+
+impl HomotopyConfig {
+    /// Sets the initial guess `x₀` for the continuation.
+    pub fn with_initial_state(mut self, state: Vec<f64>) -> Self {
+        self.initial_state = Some(state);
+        self
+    }
+
+    /// Sets the initial number of `λ` subdivisions.
+    pub fn with_initial_steps(mut self, steps: usize) -> Self {
+        self.initial_steps = steps.max(1);
+        self
+    }
+
+    /// Selects the finite-difference Jacobian (analytic is the default).
+    pub fn with_numerical_jacobian(mut self, enabled: bool) -> Self {
+        self.use_numerical_jacobian = enabled;
+        self
+    }
+
+    /// L2 norm of a residual vector.
+    fn residual_norm(residuals: &[f64]) -> f64 {
+        residuals.iter().map(|r| r * r).sum::<f64>().sqrt()
+    }
+
+    fn failure_diagnostics(
+        &self,
+        iterations: usize,
+        final_residual: f64,
+        residuals: &[f64],
+        elapsed_ms: u64,
+    ) -> Option<ConvergenceDiagnostics> {
+        if iterations == 0 {
+            return None;
+        }
+
+        let (max_residual_index, max_residual) = dominant_residual(residuals);
+        let mut diagnostics = ConvergenceDiagnostics::with_capacity(1);
+        diagnostics.iterations = iterations;
+        diagnostics.final_residual = final_residual;
+        diagnostics.best_residual = final_residual;
+        diagnostics.converged = false;
+        diagnostics.timing_ms = elapsed_ms;
+        diagnostics.final_solver = Some(SolverType::Homotopy);
+        diagnostics.push_iteration(IterationDiagnostics {
+            iteration: iterations,
+            residual_norm: final_residual,
+            delta_norm: 0.0,
+            alpha: Some(1.0),
+            jacobian_frozen: false,
+            jacobian_condition: None,
+            max_residual_index,
+            max_residual,
+        });
+        Some(diagnostics)
+    }
+
+    /// Runs the inner Newton iteration for `H(x, λ) = F(x) − (1 − λ)·r0 = 0`.
+    ///
+    /// Mutates `state` in place. Returns `Ok(iterations)` if the inner system
+    /// converged below `inner_tolerance`, or `Err(())` if it diverged, the
+    /// Jacobian was singular, or the iteration budget was exhausted. On failure
+    /// the caller restores the previous good state.
+    #[allow(clippy::too_many_arguments)]
+    fn inner_newton(
+        &self,
+        system: &mut System,
+        state: &mut [f64],
+        r0: &[f64],
+        lambda: f64,
+        clipping_mask: &[Option<(f64, f64)>],
+        residuals: &mut Vec<f64>,
+        residuals_h: &mut Vec<f64>,
+        jacobian: &mut JacobianMatrix,
+        jacobian_builder: &mut JacobianBuilder,
+    ) -> Result<usize, ()> {
+        let offset = 1.0 - lambda;
+
+        for k in 0..self.inner_max_iterations {
+            // Evaluate F(x) and form the homotopy residual H = F − (1 − λ)·r0.
+            if system.compute_residuals(state, residuals).is_err() {
+                return Err(());
+            }
+            for i in 0..residuals.len() {
+                residuals_h[i] = residuals[i] - offset * r0[i];
+            }
+
+            let norm = Self::residual_norm(residuals_h.as_slice());
+            if norm < self.inner_tolerance {
+                return Ok(k);
+            }
+            if !norm.is_finite() || norm > self.divergence_threshold {
+                return Err(());
+            }
+
+            // ∂H/∂x = ∂F/∂x, so the Jacobian of F is used unchanged.
+            if self.use_numerical_jacobian {
+                let eps = self.numerical_epsilon;
+                let compute = |s: &[f64], r: &mut [f64]| {
+                    let s_vec = s.to_vec();
+                    let mut r_vec = vec![0.0; r.len()];
+                    let res = system.compute_residuals(&s_vec, &mut r_vec);
+                    r.copy_from_slice(&r_vec);
+                    res.map(|_| ()).map_err(|e| format!("{:?}", e))
+                };
+                match JacobianMatrix::numerical(compute, state, residuals.as_slice(), eps) {
+                    Ok(jm) => jacobian.as_matrix_mut().copy_from(jm.as_matrix()),
+                    Err(_) => return Err(()),
+                }
+            } else {
+                jacobian_builder.clear();
+                if system.assemble_jacobian(state, jacobian_builder).is_err() {
+                    return Err(());
+                }
+                jacobian.update_from_builder(jacobian_builder.entries());
+            }
+
+            // Solve J·Δx = −H (the solve routine negates the supplied residual).
+            let delta = match jacobian.solve(residuals_h.as_slice()) {
+                Some(d) => d,
+                None => return Err(()),
+            };
+
+            apply_newton_step(state, &delta, clipping_mask, 1.0);
+        }
+
+        // Final convergence check after the last step.
+        if system.compute_residuals(state, residuals).is_err() {
+            return Err(());
+        }
+        for i in 0..residuals.len() {
+            residuals_h[i] = residuals[i] - offset * r0[i];
+        }
+        if Self::residual_norm(residuals_h.as_slice()) < self.inner_tolerance {
+            Ok(self.inner_max_iterations)
+        } else {
+            Err(())
+        }
+    }
+}
+
+impl Solver for HomotopyConfig {
+    fn solve(&mut self, system: &mut System) -> Result<ConvergedState, SolverError> {
+        let start_time = Instant::now();
+
+        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()
+            + 2 * system.saturated_controller_count()
+            + system.mass_flow_closure_count();
+
+        if n_state == 0 || n_equations == 0 {
+            return Err(SolverError::InvalidSystem {
+                message: "Empty system has no state variables or equations".to_string(),
+            });
+        }
+
+        // Working buffers (allocated once, reused across every λ step).
+        // A caller-supplied initial guess MUST match the system size: silently
+        // substituting zeros would hide a caller bug behind an opaque later failure.
+        let mut state: Vec<f64> = match self.initial_state.as_ref() {
+            Some(s) if s.len() == n_state => s.clone(),
+            Some(s) => {
+                return Err(SolverError::InvalidSystem {
+                    message: format!(
+                        "initial_state length {} does not match system state length {}",
+                        s.len(),
+                        n_state
+                    ),
+                });
+            }
+            None => vec![0.0; n_state],
+        };
+        let mut residuals = vec![0.0; n_equations];
+        let mut residuals_h = vec![0.0; n_equations];
+        let mut jacobian = JacobianMatrix::zeros(n_equations, n_state);
+        let mut jacobian_builder = JacobianBuilder::new();
+        let mut state_saved = vec![0.0; n_state];
+
+        let clipping_mask: Vec<Option<(f64, f64)>> = (0..n_state)
+            .map(|i| system.get_bounds_for_state_index(i))
+            .collect();
+
+        // r0 = F(x0). By construction H(x0, 0) = 0, so x0 is exact at λ = 0.
+        system
+            .compute_residuals(&state, &mut residuals)
+            .map_err(|e| SolverError::InvalidSystem {
+                message: format!("Failed to compute initial residuals: {:?}", e),
+            })?;
+        let r0 = residuals.clone();
+        let initial_norm = Self::residual_norm(&r0);
+
+        // F(x0) must be finite for the deformation H(x,λ)=F(x)-(1-λ)F(x0) to be
+        // well defined. A non-finite r0 (e.g. a zero (P,h) cold start hitting the
+        // fluid backend) would make every continuation step doomed; fail early
+        // with an actionable message instead of running the whole loop.
+        if !initial_norm.is_finite() {
+            return Err(SolverError::InvalidSystem {
+                message: "Initial residual F(x0) is non-finite; the initial guess is infeasible \
+                          for the fluid backend (provide a physical initial_state)"
+                    .to_string(),
+            });
+        }
+
+        // Already solved? Skip the continuation entirely.
+        if initial_norm < self.tolerance {
+            return Ok(ConvergedState::new(
+                state,
+                0,
+                initial_norm,
+                ConvergenceStatus::Converged,
+                SimulationMetadata::new(system.input_hash()),
+            ));
+        }
+
+        let max_step = 4.0 / self.initial_steps.max(1) as f64;
+        // Guard against a non-positive min step (e.g. struct-literal misconfig),
+        // which would otherwise let dlambda shrink forever and hang the solver.
+        let min_lambda_step = self.min_lambda_step.max(1e-12);
+        let mut lambda = 0.0_f64;
+        let mut dlambda = 1.0 / self.initial_steps.max(1) as f64;
+        let mut total_iterations = 0usize;
+
+        while lambda < 1.0 {
+            if let Some(timeout) = self.timeout {
+                if start_time.elapsed() > timeout {
+                    let compute_ok = system.compute_residuals(&state, &mut residuals).is_ok();
+                    let final_residual = if compute_ok {
+                        Self::residual_norm(&residuals)
+                    } else {
+                        f64::INFINITY
+                    };
+                    let diagnostics = self.failure_diagnostics(
+                        total_iterations,
+                        final_residual,
+                        if compute_ok { &residuals } else { &[] },
+                        start_time.elapsed().as_millis() as u64,
+                    );
+                    return Err(SolverError::Timeout {
+                        timeout_ms: timeout.as_millis() as u64,
+                    }
+                    .with_optional_diagnostics(diagnostics));
+                }
+            }
+
+            let target = (lambda + dlambda).min(1.0);
+            state_saved.copy_from_slice(&state);
+
+            match self.inner_newton(
+                system,
+                &mut state,
+                &r0,
+                target,
+                &clipping_mask,
+                &mut residuals,
+                &mut residuals_h,
+                &mut jacobian,
+                &mut jacobian_builder,
+            ) {
+                Ok(iters) => {
+                    total_iterations += iters;
+                    lambda = target;
+                    // Step succeeded: gently grow the increment for the next step.
+                    dlambda = (dlambda * 1.5).min(max_step);
+                }
+                Err(()) => {
+                    // Step failed: restore and halve the increment, then retry.
+                    state.copy_from_slice(&state_saved);
+                    dlambda *= 0.5;
+                    if dlambda < min_lambda_step {
+                        // Report the residual at the restored (last-good) state so
+                        // final_residual matches the state we actually return from.
+                        let compute_ok = system.compute_residuals(&state, &mut residuals).is_ok();
+                        let final_residual = if compute_ok {
+                            Self::residual_norm(&residuals)
+                        } else {
+                            f64::INFINITY
+                        };
+                        let diagnostics = self.failure_diagnostics(
+                            total_iterations,
+                            final_residual,
+                            if compute_ok { &residuals } else { &[] },
+                            start_time.elapsed().as_millis() as u64,
+                        );
+                        return Err(SolverError::NonConvergence {
+                            iterations: total_iterations,
+                            final_residual,
+                        }
+                        .with_optional_diagnostics(diagnostics));
+                    }
+                }
+            }
+        }
+
+        // At λ = 1, H == F: verify the real system is actually solved.
+        system
+            .compute_residuals(&state, &mut residuals)
+            .map_err(|e| SolverError::InvalidSystem {
+                message: format!("Failed to compute final residuals: {:?}", e),
+            })?;
+        let final_norm = Self::residual_norm(&residuals);
+
+        if final_norm < self.tolerance {
+            let status = if !system.saturated_variables().is_empty() {
+                ConvergenceStatus::ControlSaturation
+            } else {
+                ConvergenceStatus::Converged
+            };
+            Ok(ConvergedState::new(
+                state,
+                total_iterations,
+                final_norm,
+                status,
+                SimulationMetadata::new(system.input_hash()),
+            ))
+        } else {
+            let diagnostics = self.failure_diagnostics(
+                total_iterations,
+                final_norm,
+                &residuals,
+                start_time.elapsed().as_millis() as u64,
+            );
+            Err(SolverError::NonConvergence {
+                iterations: total_iterations,
+                final_residual: final_norm,
+            }
+            .with_optional_diagnostics(diagnostics))
+        }
+    }
+
+    fn with_timeout(self, timeout: Duration) -> Self {
+        Self {
+            timeout: Some(timeout),
+            ..self
+        }
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn test_homotopy_default_config() {
+        let cfg = HomotopyConfig::default();
+        assert_eq!(cfg.initial_steps, 10);
+        assert!(!cfg.use_numerical_jacobian);
+        assert!((cfg.tolerance - 1e-6).abs() < 1e-15);
+    }
+
+    #[test]
+    fn test_homotopy_builders() {
+        let cfg = HomotopyConfig::default()
+            .with_initial_steps(20)
+            .with_numerical_jacobian(true)
+            .with_initial_state(vec![1.0, 2.0]);
+        assert_eq!(cfg.initial_steps, 20);
+        assert!(cfg.use_numerical_jacobian);
+        assert_eq!(cfg.initial_state, Some(vec![1.0, 2.0]));
+    }
+
+    #[test]
+    fn test_homotopy_initial_steps_floor_is_one() {
+        let cfg = HomotopyConfig::default().with_initial_steps(0);
+        assert_eq!(cfg.initial_steps, 1);
+    }
+
+    #[test]
+    fn test_homotopy_residual_norm() {
+        assert!((HomotopyConfig::residual_norm(&[3.0, 4.0]) - 5.0).abs() < 1e-12);
+    }
+
+    #[test]
+    fn test_homotopy_empty_system_errors() {
+        let mut system = System::new();
+        system.finalize().unwrap();
+        let mut solver = HomotopyConfig::default();
+        assert!(solver.solve(&mut system).is_err());
+    }
+
+    #[test]
+    fn test_homotopy_with_timeout_sets_field() {
+        let cfg = HomotopyConfig::default().with_timeout(Duration::from_millis(250));
+        assert_eq!(cfg.timeout, Some(Duration::from_millis(250)));
+    }
+}

crates/solver/src/strategies/mod.rs

@@ -21,10 +21,12 @@
 //! ```
 
 mod fallback;
+mod homotopy;
 mod newton_raphson;
 mod sequential_substitution;
 
 pub use fallback::{FallbackConfig, FallbackSolver};
+pub use homotopy::HomotopyConfig;
 pub use newton_raphson::NewtonConfig;
 pub use sequential_substitution::PicardConfig;
 
@@ -83,11 +85,12 @@ impl Solver for SolverStrategy {
 
         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)?;
+                // Post-solve validation checks. Components index the state slice by
+                // global index, so pass the raw (ṁ, P, h)-strided vector directly
+                // rather than through the stride-2 SystemState conversion (CM1.2).
+                let state_slice: &[f64] = &state.state;
+                system.check_mass_balance(state_slice)?;
+                system.check_energy_balance(state_slice)?;
             }
         }
 

crates/solver/src/strategies/newton_raphson.rs

@@ -9,9 +9,9 @@ use crate::criteria::ConvergenceCriteria;
 use crate::jacobian::JacobianMatrix;
 use crate::metadata::SimulationMetadata;
 use crate::solver::{
-    apply_newton_step, ConvergedState, ConvergenceDiagnostics, ConvergenceStatus,
-    IterationDiagnostics, JacobianFreezingConfig, Solver, SolverError, SolverType,
-    TimeoutConfig, VerboseConfig,
+    apply_newton_step, dominant_residual, ConvergedState, ConvergenceDiagnostics,
+    ConvergenceStatus, IterationDiagnostics, JacobianFreezingConfig, Solver, SolverError,
+    SolverType, TimeoutConfig, VerboseConfig,
 };
 use crate::system::System;
 use entropyk_components::JacobianBuilder;
@@ -154,7 +154,10 @@ impl NewtonConfig {
     ) -> Option<SolverError> {
         if current_norm > self.divergence_threshold {
             return Some(SolverError::Divergence {
-                reason: format!("Residual {} exceeds threshold {}", current_norm, self.divergence_threshold),
+                reason: format!(
+                    "Residual {} exceeds threshold {}",
+                    current_norm, self.divergence_threshold
+                ),
             });
         }
 
@@ -162,7 +165,10 @@ impl NewtonConfig {
             *divergence_count += 1;
             if *divergence_count >= 3 {
                 return Some(SolverError::Divergence {
-                    reason: format!("Residual increased 3x: {:.6e} → {:.6e}", previous_norm, current_norm),
+                    reason: format!(
+                        "Residual increased 3x: {:.6e} → {:.6e}",
+                        previous_norm, current_norm
+                    ),
                 });
             }
         } else {
@@ -201,7 +207,12 @@ impl NewtonConfig {
 
             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");
+                tracing::debug!(
+                    alpha,
+                    old_norm = current_norm,
+                    new_norm,
+                    "Line search accepted"
+                );
                 return Some(alpha);
             }
 
@@ -209,9 +220,45 @@ impl NewtonConfig {
             alpha *= 0.5;
         }
 
-        tracing::warn!("Line search failed after {} backtracks", self.line_search_max_backtracks);
+        tracing::warn!(
+            "Line search failed after {} backtracks",
+            self.line_search_max_backtracks
+        );
         None
     }
+
+    fn finalize_failure_diagnostics(
+        &self,
+        mut diagnostics: Option<ConvergenceDiagnostics>,
+        iterations: usize,
+        final_residual: f64,
+        best_residual: f64,
+        elapsed_ms: u64,
+        jacobian_condition_final: Option<f64>,
+        final_state: Option<Vec<f64>>,
+    ) -> Option<ConvergenceDiagnostics> {
+        if let Some(ref mut diag) = diagnostics {
+            diag.iterations = iterations;
+            diag.final_residual = final_residual;
+            diag.best_residual = best_residual;
+            diag.converged = false;
+            diag.timing_ms = elapsed_ms;
+            diag.jacobian_condition_final = jacobian_condition_final;
+            diag.final_solver = Some(SolverType::NewtonRaphson);
+
+            if self.verbose_config.dump_final_state {
+                diag.final_state = final_state;
+                let json_output = diag.dump_diagnostics(self.verbose_config.output_format);
+                tracing::warn!(
+                    iterations,
+                    final_residual,
+                    "Non-convergence diagnostics:\n{}",
+                    json_output
+                );
+            }
+        }
+        diagnostics
+    }
 }
 
 impl Solver for NewtonConfig {
@@ -240,7 +287,9 @@ impl Solver for NewtonConfig {
             .map(|(_, c, _)| c.n_equations())
             .sum::<usize>()
             + system.constraints().count()
-            + system.coupling_residual_count();
+            + system.coupling_residual_count()
+            + 2 * system.saturated_controller_count()
+            + system.mass_flow_closure_count();
 
         if n_state == 0 || n_equations == 0 {
             return Err(SolverError::InvalidSystem {
@@ -248,15 +297,22 @@ impl Solver for NewtonConfig {
             });
         }
 
-        // 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]);
+        // Pre-allocate all buffers. A caller-supplied initial state MUST match
+        // the full state length: a debug_assert would abort (violating zero-panic)
+        // and a silent zeros fallback would solve a different problem. Fail cleanly.
+        let mut state: Vec<f64> = match self.initial_state.as_ref() {
+            Some(s) if s.len() == n_state => s.clone(),
+            Some(s) => {
+                return Err(SolverError::InvalidSystem {
+                    message: format!(
+                        "initial_state length {} does not match system state length {}",
+                        s.len(),
+                        n_state
+                    ),
+                });
+            }
+            None => 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;
@@ -273,13 +329,13 @@ impl Solver for NewtonConfig {
         let mut jacobian_matrix = JacobianMatrix::zeros(n_equations, n_state);
         let mut frozen_count: usize = 0;
         let mut force_recompute: bool = true;
-        
+
         // Cached condition number (for verbose mode when Jacobian frozen)
         let mut cached_condition: Option<f64> = None;
 
         // Pre-compute clipping mask
         let clipping_mask: Vec<Option<(f64, f64)>> = (0..n_state)
-            .map(|i| system.get_bounds_for_state_index(i))
+            .map(|i| system.get_solver_bounds_for_state_index(i))
             .collect();
 
         // Initial residual computation
@@ -306,15 +362,32 @@ impl Solver for NewtonConfig {
             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)");
+                    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()),
+                        state,
+                        0,
+                        current_norm,
+                        status,
+                        report,
+                        SimulationMetadata::new(system.input_hash()),
                     ));
                 }
             } else {
-                tracing::info!(iterations = 0, final_residual = current_norm, "Converged at initial state");
+                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()),
+                    state,
+                    0,
+                    current_norm,
+                    status,
+                    SimulationMetadata::new(system.input_hash()),
                 ));
             }
         }
@@ -327,7 +400,18 @@ impl Solver for NewtonConfig {
             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);
+                    let failure_diagnostics = self.finalize_failure_diagnostics(
+                        diagnostics.take(),
+                        iteration - 1,
+                        current_norm,
+                        best_residual,
+                        start_time.elapsed().as_millis() as u64,
+                        cached_condition,
+                        Some(state.clone()),
+                    );
+                    return self
+                        .handle_timeout(&best_state, best_residual, iteration - 1, timeout, system)
+                        .map_err(|err| err.with_optional_diagnostics(failure_diagnostics));
                 }
             }
 
@@ -346,7 +430,7 @@ impl Solver for NewtonConfig {
             };
 
             let jacobian_frozen_this_iter = !should_recompute;
-            
+
             if should_recompute {
                 // Fresh Jacobian assembly (in-place update)
                 jacobian_builder.clear();
@@ -359,13 +443,15 @@ impl Solver for NewtonConfig {
                         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 {
+                    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)
+                    system
+                        .assemble_jacobian(&state, &mut jacobian_builder)
                         .map_err(|e| SolverError::InvalidSystem {
                             message: format!("Failed to assemble Jacobian: {:?}", e),
                         })?;
@@ -374,19 +460,27 @@ impl Solver for NewtonConfig {
 
                 frozen_count = 0;
                 force_recompute = false;
-                
+
                 // Compute and cache condition number if verbose mode enabled
                 if verbose_enabled && self.verbose_config.log_jacobian_condition {
                     let cond = jacobian_matrix.estimate_condition_number();
                     cached_condition = cond;
                     if let Some(c) = cond {
-                        tracing::info!(iteration, condition_number = c, "Jacobian condition number");
+                        tracing::info!(
+                            iteration,
+                            condition_number = c,
+                            "Jacobian condition number"
+                        );
                         if c > 1e10 {
-                            tracing::warn!(iteration, condition_number = c, "Ill-conditioned Jacobian detected (κ > 1e10)");
+                            tracing::warn!(
+                                iteration,
+                                condition_number = c,
+                                "Ill-conditioned Jacobian detected (κ > 1e10)"
+                            );
                         }
                     }
                 }
-                
+
                 tracing::debug!(iteration, "Fresh Jacobian computed");
             } else {
                 frozen_count += 1;
@@ -397,23 +491,49 @@ impl Solver for NewtonConfig {
             let delta = match jacobian_matrix.solve(&residuals) {
                 Some(d) => d,
                 None => {
+                    let failure_diagnostics = self.finalize_failure_diagnostics(
+                        diagnostics.take(),
+                        iteration,
+                        current_norm,
+                        best_residual,
+                        start_time.elapsed().as_millis() as u64,
+                        cached_condition,
+                        Some(state.clone()),
+                    );
                     return Err(SolverError::Divergence {
                         reason: "Jacobian is singular".to_string(),
-                    });
+                    }
+                    .with_optional_diagnostics(failure_diagnostics));
                 }
             };
 
             // 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,
+                    system,
+                    &mut state,
+                    &delta,
+                    &residuals,
+                    current_norm,
+                    &mut state_copy,
+                    &mut new_residuals,
+                    &clipping_mask,
                 ) {
                     Some(a) => a,
                     None => {
+                        let failure_diagnostics = self.finalize_failure_diagnostics(
+                            diagnostics.take(),
+                            iteration,
+                            current_norm,
+                            best_residual,
+                            start_time.elapsed().as_millis() as u64,
+                            cached_condition,
+                            Some(state.clone()),
+                        );
                         return Err(SolverError::Divergence {
                             reason: "Line search failed".to_string(),
-                        });
+                        }
+                        .with_optional_diagnostics(failure_diagnostics));
                     }
                 }
             } else {
@@ -421,16 +541,18 @@ impl Solver for NewtonConfig {
                 1.0
             };
 
-            system.compute_residuals(&state, &mut 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);
-            
+
             // Compute delta norm for diagnostics
-            let delta_norm: f64 = state.iter()
+            let delta_norm: f64 = state
+                .iter()
                 .zip(prev_iteration_state.iter())
                 .map(|(s, p)| (s - p).powi(2))
                 .sum::<f64>()
@@ -444,9 +566,16 @@ impl Solver for NewtonConfig {
 
             // 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 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");
+                        tracing::debug!(
+                            iteration,
+                            current_norm,
+                            previous_norm,
+                            "Unfreezing Jacobian"
+                        );
                     }
                     force_recompute = true;
                     frozen_count = 0;
@@ -464,9 +593,10 @@ impl Solver for NewtonConfig {
                     "Newton iteration"
                 );
             }
-            
+
             // Collect iteration diagnostics
             if let Some(ref mut diag) = diagnostics {
+                let (max_residual_index, max_residual) = dominant_residual(&residuals);
                 diag.push_iteration(IterationDiagnostics {
                     iteration,
                     residual_norm: current_norm,
@@ -474,21 +604,29 @@ impl Solver for NewtonConfig {
                     alpha: Some(alpha),
                     jacobian_frozen: jacobian_frozen_this_iter,
                     jacobian_condition: cached_condition,
+                    max_residual_index,
+                    max_residual,
                 });
             }
 
-            tracing::debug!(iteration, residual_norm = current_norm, alpha, "Newton iteration complete");
+            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);
+                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
                     };
-                    
+
                     // Finalize diagnostics
                     if let Some(ref mut diag) = diagnostics {
                         diag.iterations = iteration;
@@ -498,19 +636,33 @@ impl Solver for NewtonConfig {
                         diag.timing_ms = start_time.elapsed().as_millis() as u64;
                         diag.jacobian_condition_final = cached_condition;
                         diag.final_solver = Some(SolverType::NewtonRaphson);
-                        
+
                         if self.verbose_config.log_residuals {
                             tracing::info!("{}", diag.summary());
                         }
                     }
-                    
-                    tracing::info!(iterations = iteration, final_residual = current_norm, "Converged (criteria)");
+
+                    tracing::info!(
+                        iterations = iteration,
+                        final_residual = current_norm,
+                        "Converged (criteria)"
+                    );
                     let result = ConvergedState::with_report(
-                        state, iteration, current_norm, status, report, SimulationMetadata::new(system.input_hash()),
+                        state,
+                        iteration,
+                        current_norm,
+                        status,
+                        report,
+                        SimulationMetadata::new(system.input_hash()),
                     );
                     return Ok(if let Some(d) = diagnostics {
-                        ConvergedState { diagnostics: Some(d), ..result }
-                    } else { result });
+                        ConvergedState {
+                            diagnostics: Some(d),
+                            ..result
+                        }
+                    } else {
+                        result
+                    });
                 }
                 false
             } else {
@@ -523,7 +675,7 @@ impl Solver for NewtonConfig {
                 } else {
                     ConvergenceStatus::Converged
                 };
-                
+
                 // Finalize diagnostics
                 if let Some(ref mut diag) = diagnostics {
                     diag.iterations = iteration;
@@ -533,54 +685,76 @@ impl Solver for NewtonConfig {
                     diag.timing_ms = start_time.elapsed().as_millis() as u64;
                     diag.jacobian_condition_final = cached_condition;
                     diag.final_solver = Some(SolverType::NewtonRaphson);
-                    
+
                     if self.verbose_config.log_residuals {
                         tracing::info!("{}", diag.summary());
                     }
                 }
-                
-                tracing::info!(iterations = iteration, final_residual = current_norm, "Converged");
+
+                tracing::info!(
+                    iterations = iteration,
+                    final_residual = current_norm,
+                    "Converged"
+                );
                 let result = ConvergedState::new(
-                    state, iteration, current_norm, status, SimulationMetadata::new(system.input_hash()),
+                    state,
+                    iteration,
+                    current_norm,
+                    status,
+                    SimulationMetadata::new(system.input_hash()),
                 );
                 return Ok(if let Some(d) = diagnostics {
-                    ConvergedState { diagnostics: Some(d), ..result }
-                } else { result });
+                    ConvergedState {
+                        diagnostics: Some(d),
+                        ..result
+                    }
+                } else {
+                    result
+                });
             }
 
-            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);
+            if let Some(err) =
+                self.check_divergence(current_norm, previous_norm, &mut divergence_count)
+            {
+                tracing::warn!(
+                    iteration,
+                    residual_norm = current_norm,
+                    "Divergence detected"
+                );
+                let failure_diagnostics = self.finalize_failure_diagnostics(
+                    diagnostics.take(),
+                    iteration,
+                    current_norm,
+                    best_residual,
+                    start_time.elapsed().as_millis() as u64,
+                    cached_condition,
+                    Some(state.clone()),
+                );
+                return Err(err.with_optional_diagnostics(failure_diagnostics));
             }
         }
 
         // Non-convergence: dump diagnostics if enabled
-        if let Some(ref mut diag) = diagnostics {
-            diag.iterations = self.max_iterations;
-            diag.final_residual = current_norm;
-            diag.best_residual = best_residual;
-            diag.converged = false;
-            diag.timing_ms = start_time.elapsed().as_millis() as u64;
-            diag.jacobian_condition_final = cached_condition;
-            diag.final_solver = Some(SolverType::NewtonRaphson);
-            
-            if self.verbose_config.dump_final_state {
-                diag.final_state = Some(state.clone());
-                let json_output = diag.dump_diagnostics(self.verbose_config.output_format);
-                tracing::warn!(
-                    iterations = self.max_iterations,
-                    final_residual = current_norm,
-                    "Non-convergence diagnostics:\n{}",
-                    json_output
-                );
-            }
-        }
+        let failure_diagnostics = self.finalize_failure_diagnostics(
+            diagnostics.take(),
+            self.max_iterations,
+            current_norm,
+            best_residual,
+            start_time.elapsed().as_millis() as u64,
+            cached_condition,
+            Some(state.clone()),
+        );
 
-        tracing::warn!(max_iterations = self.max_iterations, final_residual = current_norm, "Did not converge");
+        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,
-        })
+        }
+        .with_optional_diagnostics(failure_diagnostics))
     }
 
     fn with_timeout(mut self, timeout: Duration) -> Self {

crates/solver/src/strategies/sequential_substitution.rs

@@ -3,13 +3,16 @@
 //! Provides [`PicardConfig`] which implements Picard iteration for solving
 //! systems of non-linear equations. Slower than Newton-Raphson but more robust.
 
+use std::collections::VecDeque;
 use std::time::{Duration, Instant};
 
+use nalgebra::{DMatrix, DVector};
+
 use crate::criteria::ConvergenceCriteria;
 use crate::metadata::SimulationMetadata;
 use crate::solver::{
-    ConvergedState, ConvergenceDiagnostics, ConvergenceStatus, IterationDiagnostics, Solver,
-    SolverError, SolverType, TimeoutConfig, VerboseConfig,
+    dominant_residual, ConvergedState, ConvergenceDiagnostics, ConvergenceStatus,
+    IterationDiagnostics, Solver, SolverError, SolverType, TimeoutConfig, VerboseConfig,
 };
 use crate::system::System;
 
@@ -43,6 +46,13 @@ pub struct PicardConfig {
     pub convergence_criteria: Option<ConvergenceCriteria>,
     /// Verbose mode configuration for diagnostics.
     pub verbose_config: VerboseConfig,
+    /// Anderson acceleration depth `m` (history window). `0` disables acceleration
+    /// and the solver behaves as plain relaxed Picard (default). Typical useful
+    /// values are 3–5. See [`PicardConfig::with_anderson`].
+    pub anderson_depth: usize,
+    /// Tikhonov regularization added to the Anderson least-squares normal matrix
+    /// for numerical stability. Default: 1e-10. Only used when `anderson_depth > 0`.
+    pub anderson_regularization: f64,
 }
 
 impl Default for PicardConfig {
@@ -60,6 +70,8 @@ impl Default for PicardConfig {
             initial_state: None,
             convergence_criteria: None,
             verbose_config: VerboseConfig::default(),
+            anderson_depth: 0,
+            anderson_regularization: 1e-10,
         }
     }
 }
@@ -90,6 +102,23 @@ impl PicardConfig {
         self
     }
 
+    /// Enables Anderson acceleration with history depth `m` (Story: solver speed).
+    ///
+    /// Anderson acceleration (Walker & Ni, 2011) turns the linearly-convergent
+    /// relaxed Picard fixed-point iteration into a super-linearly convergent one by
+    /// extrapolating from the last `m` residual/map-value pairs via a small
+    /// least-squares problem. `m = 0` disables it (plain relaxed Picard). Values of
+    /// 3–5 typically cut the iteration count by 2–3× on stiff refrigeration cycles
+    /// while adding only an `O(m² · n)` least-squares solve per iteration.
+    ///
+    /// # Reference
+    /// Walker, H.F., Ni, P. (2011). "Anderson acceleration for fixed-point
+    /// iterations." *SIAM J. Numerical Analysis*, 49(4):1715–1735.
+    pub fn with_anderson(mut self, depth: usize) -> Self {
+        self.anderson_depth = depth;
+        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()
@@ -200,6 +229,37 @@ impl PicardConfig {
             *x -= omega * r;
         }
     }
+
+    fn finalize_failure_diagnostics(
+        &self,
+        mut diagnostics: Option<ConvergenceDiagnostics>,
+        iterations: usize,
+        final_residual: f64,
+        best_residual: f64,
+        elapsed_ms: u64,
+        final_state: Option<Vec<f64>>,
+    ) -> Option<ConvergenceDiagnostics> {
+        if let Some(ref mut diag) = diagnostics {
+            diag.iterations = iterations;
+            diag.final_residual = final_residual;
+            diag.best_residual = best_residual;
+            diag.converged = false;
+            diag.timing_ms = elapsed_ms;
+            diag.final_solver = Some(SolverType::SequentialSubstitution);
+
+            if self.verbose_config.dump_final_state {
+                diag.final_state = final_state;
+                let json_output = diag.dump_diagnostics(self.verbose_config.output_format);
+                tracing::warn!(
+                    iterations,
+                    final_residual,
+                    "Non-convergence diagnostics:\n{}",
+                    json_output
+                );
+            }
+        }
+        diagnostics
+    }
 }
 
 impl Solver for PicardConfig {
@@ -231,7 +291,9 @@ impl Solver for PicardConfig {
             .map(|(_, c, _)| c.n_equations())
             .sum::<usize>()
             + system.constraints().count()
-            + system.coupling_residual_count();
+            + system.coupling_residual_count()
+            + 2 * system.saturated_controller_count()
+            + system.mass_flow_closure_count();
 
         // Validate system
         if n_state == 0 || n_equations == 0 {
@@ -251,25 +313,22 @@ impl Solver for PicardConfig {
         }
 
         // 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]);
+        // Story 4.6 - AC: #8: Use initial_state if provided, else start from zeros.
+        // A mismatched length is a hard error (zero-panic; no silent zeros fallback
+        // that would solve a different problem) — consistent with Newton/Homotopy.
+        let mut state: Vec<f64> = match self.initial_state.as_ref() {
+            Some(s) if s.len() == n_state => s.clone(),
+            Some(s) => {
+                return Err(SolverError::InvalidSystem {
+                    message: format!(
+                        "initial_state length {} does not match system state length {}",
+                        s.len(),
+                        n_state
+                    ),
+                });
+            }
+            None => 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;
@@ -310,6 +369,16 @@ impl Solver for PicardConfig {
             ));
         }
 
+        // Optional Anderson accelerator (disabled when depth == 0).
+        let mut anderson = if self.anderson_depth > 0 {
+            Some(AndersonAccelerator::new(
+                self.anderson_depth,
+                self.anderson_regularization,
+            ))
+        } else {
+            None
+        };
+
         // Main Picard iteration loop
         for iteration in 1..=self.max_iterations {
             // Save state before step for convergence criteria delta checks
@@ -327,18 +396,28 @@ impl Solver for PicardConfig {
                     );
 
                     // Story 4.5 - AC: #2, #6: Return best state or error based on config
-                    return self.handle_timeout(
-                        &best_state,
-                        best_residual,
+                    let failure_diagnostics = self.finalize_failure_diagnostics(
+                        diagnostics.take(),
                         iteration - 1,
-                        timeout,
-                        system,
+                        current_norm,
+                        best_residual,
+                        start_time.elapsed().as_millis() as u64,
+                        Some(state.clone()),
                     );
+                    return self
+                        .handle_timeout(&best_state, best_residual, iteration - 1, timeout, system)
+                        .map_err(|err| err.with_optional_diagnostics(failure_diagnostics));
                 }
             }
 
-            // Apply relaxed update: x_new = x_old - omega * residual (AC: #2, #3)
-            Self::apply_relaxation(&mut state, &residuals, self.relaxation_factor);
+            // Apply update. With Anderson acceleration enabled, extrapolate from the
+            // residual/map-value history; otherwise use plain relaxed Picard.
+            // Both share the same underlying fixed-point map G(x) = x - ω·F(x).
+            if let Some(acc) = anderson.as_mut() {
+                acc.next_state_into(&mut state, &residuals, self.relaxation_factor);
+            } else {
+                Self::apply_relaxation(&mut state, &residuals, self.relaxation_factor);
+            }
 
             // Compute new residuals
             system
@@ -349,9 +428,10 @@ impl Solver for PicardConfig {
 
             previous_norm = current_norm;
             current_norm = Self::residual_norm(&residuals);
-            
+
             // Compute delta norm for diagnostics
-            let delta_norm: f64 = state.iter()
+            let delta_norm: f64 = state
+                .iter()
                 .zip(prev_iteration_state.iter())
                 .map(|(s, p)| (s - p).powi(2))
                 .sum::<f64>()
@@ -378,16 +458,19 @@ impl Solver for PicardConfig {
                     "Picard iteration"
                 );
             }
-            
+
             // Collect iteration diagnostics
             if let Some(ref mut diag) = diagnostics {
+                let (max_residual_index, max_residual) = dominant_residual(&residuals);
                 diag.push_iteration(IterationDiagnostics {
                     iteration,
                     residual_norm: current_norm,
                     delta_norm,
-                    alpha: None, // Picard doesn't use line search
-                    jacobian_frozen: false, // Picard doesn't use Jacobian
+                    alpha: None,              // Picard doesn't use line search
+                    jacobian_frozen: false,   // Picard doesn't use Jacobian
                     jacobian_condition: None, // No Jacobian in Picard
+                    max_residual_index,
+                    max_residual,
                 });
             }
 
@@ -411,12 +494,12 @@ impl Solver for PicardConfig {
                         diag.converged = true;
                         diag.timing_ms = start_time.elapsed().as_millis() as u64;
                         diag.final_solver = Some(SolverType::SequentialSubstitution);
-                        
+
                         if self.verbose_config.log_residuals {
                             tracing::info!("{}", diag.summary());
                         }
                     }
-                    
+
                     tracing::info!(
                         iterations = iteration,
                         final_residual = current_norm,
@@ -432,8 +515,13 @@ impl Solver for PicardConfig {
                         SimulationMetadata::new(system.input_hash()),
                     );
                     return Ok(if let Some(d) = diagnostics {
-                        ConvergedState { diagnostics: Some(d), ..result }
-                    } else { result });
+                        ConvergedState {
+                            diagnostics: Some(d),
+                            ..result
+                        }
+                    } else {
+                        result
+                    });
                 }
                 false
             } else {
@@ -449,12 +537,12 @@ impl Solver for PicardConfig {
                     diag.converged = true;
                     diag.timing_ms = start_time.elapsed().as_millis() as u64;
                     diag.final_solver = Some(SolverType::SequentialSubstitution);
-                    
+
                     if self.verbose_config.log_residuals {
                         tracing::info!("{}", diag.summary());
                     }
                 }
-                
+
                 tracing::info!(
                     iterations = iteration,
                     final_residual = current_norm,
@@ -469,8 +557,13 @@ impl Solver for PicardConfig {
                     SimulationMetadata::new(system.input_hash()),
                 );
                 return Ok(if let Some(d) = diagnostics {
-                    ConvergedState { diagnostics: Some(d), ..result }
-                } else { result });
+                    ConvergedState {
+                        diagnostics: Some(d),
+                        ..result
+                    }
+                } else {
+                    result
+                });
             }
 
             // Check divergence (AC: #5)
@@ -482,30 +575,27 @@ impl Solver for PicardConfig {
                     residual_norm = current_norm,
                     "Divergence detected"
                 );
-                return Err(err);
+                let failure_diagnostics = self.finalize_failure_diagnostics(
+                    diagnostics.take(),
+                    iteration,
+                    current_norm,
+                    best_residual,
+                    start_time.elapsed().as_millis() as u64,
+                    Some(state.clone()),
+                );
+                return Err(err.with_optional_diagnostics(failure_diagnostics));
             }
         }
 
         // Non-convergence: dump diagnostics if enabled
-        if let Some(ref mut diag) = diagnostics {
-            diag.iterations = self.max_iterations;
-            diag.final_residual = current_norm;
-            diag.best_residual = best_residual;
-            diag.converged = false;
-            diag.timing_ms = start_time.elapsed().as_millis() as u64;
-            diag.final_solver = Some(SolverType::SequentialSubstitution);
-            
-            if self.verbose_config.dump_final_state {
-                diag.final_state = Some(state.clone());
-                let json_output = diag.dump_diagnostics(self.verbose_config.output_format);
-                tracing::warn!(
-                    iterations = self.max_iterations,
-                    final_residual = current_norm,
-                    "Non-convergence diagnostics:\n{}",
-                    json_output
-                );
-            }
-        }
+        let failure_diagnostics = self.finalize_failure_diagnostics(
+            diagnostics.take(),
+            self.max_iterations,
+            current_norm,
+            best_residual,
+            start_time.elapsed().as_millis() as u64,
+            Some(state.clone()),
+        );
 
         // Max iterations exceeded
         tracing::warn!(
@@ -516,7 +606,8 @@ impl Solver for PicardConfig {
         Err(SolverError::NonConvergence {
             iterations: self.max_iterations,
             final_residual: current_norm,
-        })
+        }
+        .with_optional_diagnostics(failure_diagnostics))
     }
 
     fn with_timeout(mut self, timeout: Duration) -> Self {
@@ -525,6 +616,110 @@ impl Solver for PicardConfig {
     }
 }
 
+/// Anderson acceleration state for the relaxed Picard fixed-point iteration.
+///
+/// The underlying fixed-point map is `G(x) = x - ω·F(x)` where `F` is the residual
+/// vector and `ω` the relaxation factor. Define the map residual `f(x) = G(x) - x =
+/// -ω·F(x)`. Anderson acceleration maintains the last `m` differences of `f` and `G`
+/// and, each iteration, solves the small least-squares problem
+/// `min_γ ‖f_k - ΔF·γ‖` then sets `x_{k+1} = G_k - ΔG·γ` (Walker & Ni, 2011,
+/// following H. Walker's reference `anderson.m`). With `m = 0` (empty history) it
+/// reduces exactly to the plain step `x_{k+1} = G_k`.
+struct AndersonAccelerator {
+    depth: usize,
+    regularization: f64,
+    /// Previous map-residual f = G(x) - x.
+    f_prev: Option<Vec<f64>>,
+    /// Previous map value G(x).
+    g_prev: Option<Vec<f64>>,
+    /// History of Δf columns (most-recent at back), capped at `depth`.
+    df: VecDeque<Vec<f64>>,
+    /// History of ΔG columns (most-recent at back), capped at `depth`.
+    dg: VecDeque<Vec<f64>>,
+}
+
+impl AndersonAccelerator {
+    fn new(depth: usize, regularization: f64) -> Self {
+        Self {
+            depth,
+            regularization,
+            f_prev: None,
+            g_prev: None,
+            df: VecDeque::with_capacity(depth),
+            dg: VecDeque::with_capacity(depth),
+        }
+    }
+
+    /// Advances `state` in place from `x_k` to the accelerated `x_{k+1}`, given the
+    /// current residual vector `F(x_k)` and relaxation factor `ω`.
+    fn next_state_into(&mut self, state: &mut [f64], residual: &[f64], omega: f64) {
+        let n = state.len();
+        // Map residual f = -ω·F and fixed-point map value G = x + f.
+        let fval: Vec<f64> = residual.iter().map(|r| -omega * r).collect();
+        let gval: Vec<f64> = state.iter().zip(&fval).map(|(x, f)| x + f).collect();
+
+        // Push newest history differences.
+        if let (Some(fp), Some(gp)) = (self.f_prev.as_ref(), self.g_prev.as_ref()) {
+            let df_col: Vec<f64> = fval.iter().zip(fp).map(|(a, b)| a - b).collect();
+            let dg_col: Vec<f64> = gval.iter().zip(gp).map(|(a, b)| a - b).collect();
+            self.df.push_back(df_col);
+            self.dg.push_back(dg_col);
+            while self.df.len() > self.depth {
+                self.df.pop_front();
+                self.dg.pop_front();
+            }
+        }
+        self.f_prev = Some(fval.clone());
+        self.g_prev = Some(gval.clone());
+
+        let m = self.df.len();
+        if m == 0 {
+            // No history yet — plain relaxed step.
+            state.copy_from_slice(&gval);
+            return;
+        }
+
+        // Solve the small least-squares problem for γ via regularized normal
+        // equations: (ΔFᵀΔF + λI)·γ = ΔFᵀ·f_k. `m` is at most `depth` (small).
+        let mut ata = DMatrix::<f64>::zeros(m, m);
+        let mut atb = DVector::<f64>::zeros(m);
+        for i in 0..m {
+            for j in i..m {
+                let mut s = 0.0;
+                for k in 0..n {
+                    s += self.df[i][k] * self.df[j][k];
+                }
+                ata[(i, j)] = s;
+                ata[(j, i)] = s;
+            }
+            ata[(i, i)] += self.regularization;
+            let mut s = 0.0;
+            for k in 0..n {
+                s += self.df[i][k] * fval[k];
+            }
+            atb[i] = s;
+        }
+
+        let gamma = match ata.clone().lu().solve(&atb) {
+            Some(g) => g,
+            None => {
+                // Singular even with regularization — fall back to plain step.
+                state.copy_from_slice(&gval);
+                return;
+            }
+        };
+
+        // x_{k+1} = G_k - ΔG·γ.
+        for k in 0..n {
+            let mut acc = gval[k];
+            for (i, g) in gamma.iter().enumerate() {
+                acc -= g * self.dg[i][k];
+            }
+            state[k] = acc;
+        }
+    }
+}
+
 #[cfg(test)]
 mod tests {
     use super::*;
@@ -570,4 +765,99 @@ mod tests {
         system.finalize().unwrap();
         assert!(boxed.solve(&mut system).is_err());
     }
+
+    // ── Anderson acceleration ────────────────────────────────────────────────
+
+    /// Reference linear residual F(x) = A·x - b. Its unique root is x* = A⁻¹·b.
+    /// The relaxed Picard map is x_{k+1} = x_k - ω·(A·x_k - b).
+    fn linear_residual(a: &[[f64; 2]; 2], b: &[f64; 2], x: &[f64]) -> Vec<f64> {
+        vec![
+            a[0][0] * x[0] + a[0][1] * x[1] - b[0],
+            a[1][0] * x[0] + a[1][1] * x[1] - b[1],
+        ]
+    }
+
+    fn residual_norm2(r: &[f64]) -> f64 {
+        r.iter().map(|v| v * v).sum::<f64>().sqrt()
+    }
+
+    #[test]
+    fn test_anderson_depth_zero_matches_plain_relaxation() {
+        // With no history, next_state_into must equal x - ω·F(x).
+        let mut acc = AndersonAccelerator::new(0, 1e-10);
+        let mut state = vec![10.0, 20.0];
+        let residuals = vec![1.0, 2.0];
+        acc.next_state_into(&mut state, &residuals, 0.5);
+        assert!((state[0] - 9.5).abs() < 1e-15);
+        assert!((state[1] - 19.0).abs() < 1e-15);
+    }
+
+    #[test]
+    fn test_anderson_converges_faster_than_plain_picard() {
+        // Stiff-ish SPD system where plain relaxed Picard converges slowly.
+        let a = [[8.0, 1.0], [1.0, 3.0]];
+        let b = [9.0, 4.0]; // exact root x* = [1, 1]
+        let omega = 0.12; // deliberately small → slow plain Picard
+        let tol = 1e-9;
+        let max_iter = 2000;
+
+        let count_iters = |depth: usize| -> (usize, Vec<f64>) {
+            let mut state = vec![0.0, 0.0];
+            let mut acc = AndersonAccelerator::new(depth, 1e-12);
+            for it in 1..=max_iter {
+                let r = linear_residual(&a, &b, &state);
+                if residual_norm2(&r) < tol {
+                    return (it - 1, state);
+                }
+                acc.next_state_into(&mut state, &r, omega);
+            }
+            (max_iter, state)
+        };
+
+        let (plain_iters, _) = count_iters(0);
+        let (anderson_iters, sol) = count_iters(3);
+
+        // Anderson must converge, land on the true root, and use far fewer steps.
+        assert!(anderson_iters < max_iter, "Anderson did not converge");
+        assert!((sol[0] - 1.0).abs() < 1e-6 && (sol[1] - 1.0).abs() < 1e-6);
+        assert!(
+            anderson_iters * 3 < plain_iters,
+            "Anderson ({}) should be much faster than plain Picard ({})",
+            anderson_iters,
+            plain_iters
+        );
+    }
+
+    #[test]
+    fn test_anderson_solves_where_plain_diverges_marginally() {
+        // Anderson should still hit the exact root of a well-posed linear system.
+        let a = [[4.0, 1.0], [2.0, 5.0]];
+        let b = [6.0, 9.0];
+        // exact root: solve → x=[1, 1.4? ] compute: 4x+y=6, 2x+5y=9
+        // From first: y = 6-4x; sub: 2x+5(6-4x)=9 → 2x+30-20x=9 → -18x=-21 → x=7/6
+        // y = 6-4*7/6 = 6-28/6 = 8/6 = 4/3
+        let omega = 0.15;
+        let mut state = vec![0.0, 0.0];
+        let mut acc = AndersonAccelerator::new(4, 1e-12);
+        let mut converged = false;
+        for _ in 0..5000 {
+            let r = linear_residual(&a, &b, &state);
+            if residual_norm2(&r) < 1e-9 {
+                converged = true;
+                break;
+            }
+            acc.next_state_into(&mut state, &r, omega);
+        }
+        assert!(converged);
+        assert!((state[0] - 7.0 / 6.0).abs() < 1e-6);
+        assert!((state[1] - 4.0 / 3.0).abs() < 1e-6);
+    }
+
+    #[test]
+    fn test_with_anderson_builder_sets_depth() {
+        let cfg = PicardConfig::default().with_anderson(5);
+        assert_eq!(cfg.anderson_depth, 5);
+        // Default remains disabled.
+        assert_eq!(PicardConfig::default().anderson_depth, 0);
+    }
 }

crates/solver/src/topology.rs

@@ -0,0 +1,438 @@
+//! Mass-flow topology presolve for the system graph (CM1.4).
+//!
+//! This module identifies **series branches** — maximal sequences of flow edges
+//! where every intermediate node has exactly one incoming and one outgoing edge
+//! (no junction). All edges in such a branch share a single mass-flow Newton
+//! unknown, reducing the state-vector size from `3|E|` to `|B| + 2|E|`.
+//!
+//! ## Algorithm
+//!
+//! 1. Iterate over all graph edges in petgraph iteration order.
+//! 2. For each unvisited edge, trace the maximal series branch by walking
+//!    forward (following the target node's single outgoing edge) and backward
+//!    (following the source node's single incoming edge) while nodes have
+//!    in-degree == 1 and out-degree == 1.
+//! 3. Assign the same `mass_flow_branch_id` to every edge in the branch.
+//! 4. Return the total branch count (`|B|`).
+//!
+//! ## Design reference
+//!
+//! Based on TESPy `presolve_massflow_topology()` in
+//! `src/tespy/networks/network.py` (lines 983–1177).
+//! See `complex-machine-design.md §1.6` and `epics-complex-machine.md` Story 1.4.
+
+use entropyk_components::Component;
+use petgraph::graph::{EdgeIndex, Graph};
+use petgraph::visit::EdgeRef;
+use petgraph::Directed;
+use std::collections::HashSet;
+
+use crate::system::FlowEdge;
+
+/// Assigns `mass_flow_branch_id` to every edge in `graph` by grouping edges
+/// into series branches (maximal paths with no junction node).
+///
+/// Returns the number of independent branches `|B|`. Every edge in the same
+/// branch receives the same `mass_flow_branch_id`; edges in different branches
+/// receive distinct ids.
+///
+/// # Series branch definition
+///
+/// Two adjacent edges belong to the same branch when their shared node has
+/// **in-degree == 1 and out-degree == 1** in the directed graph. Nodes with
+/// more than one incoming or outgoing edge are junction boundaries that start
+/// a new branch.
+///
+/// For a pure series directed cycle (every node in-degree=1, out-degree=1), all
+/// edges form one branch. The walk terminates when it revisits a branch-member
+/// edge (cycle guard).
+pub fn presolve_mass_flow_topology(
+    graph: &mut Graph<Box<dyn Component>, FlowEdge, Directed>,
+) -> usize {
+    let mut visited: HashSet<EdgeIndex> = HashSet::new();
+    let mut branch_id: usize = 0;
+
+    // Collect edge indices first to avoid borrow conflict on the mutable graph.
+    let edge_indices: Vec<EdgeIndex> = graph.edge_indices().collect();
+
+    for start_edge in edge_indices {
+        if visited.contains(&start_edge) {
+            continue;
+        }
+
+        // Trace the maximal series branch containing start_edge.
+        let branch = trace_series_branch(graph, start_edge, &visited);
+
+        // Assign branch_id and mark all branch edges as visited.
+        for &eid in &branch {
+            if let Some(weight) = graph.edge_weight_mut(eid) {
+                weight.mass_flow_branch_id = branch_id;
+            }
+            visited.insert(eid);
+        }
+
+        branch_id += 1;
+    }
+
+    branch_id
+}
+
+/// Traces the maximal series branch containing `start_edge`.
+///
+/// Walks forward from the target node and backward from the source node of
+/// `start_edge`, collecting all edges reachable through 1-in-1-out nodes **or
+/// through declared internal flow paths** ([`Component::flow_paths`]) of
+/// multi-port components: a Modelica-style 4-port heat exchanger declares
+/// `[(0, 1), (2, 3)]`, so the refrigerant inlet/outlet edges (ports 0→1) and
+/// the secondary inlet/outlet edges (ports 2→3) each form one continuous
+/// series branch, exactly like `port_a.m_flow + port_b.m_flow = 0`.
+///
+/// Stops when a genuine junction node (splitter/merger/drum without a matching
+/// flow path) or a previously-visited/branch edge is encountered.
+fn trace_series_branch(
+    graph: &Graph<Box<dyn Component>, FlowEdge, Directed>,
+    start_edge: EdgeIndex,
+    visited: &HashSet<EdgeIndex>,
+) -> Vec<EdgeIndex> {
+    let mut branch: Vec<EdgeIndex> = vec![start_edge];
+
+    let (src, tgt) = graph
+        .edge_endpoints(start_edge)
+        .expect("start_edge must be a valid edge");
+
+    // ── Walk forward from the target node ──────────────────────────────────
+    let mut cur = tgt;
+    let mut cur_in = start_edge;
+    loop {
+        let next_eid = match forward_continuation(graph, cur, cur_in) {
+            Some(e) => e,
+            None => break,
+        };
+
+        // Cycle guard and already-visited guard.
+        if branch.contains(&next_eid) || visited.contains(&next_eid) {
+            break;
+        }
+
+        branch.push(next_eid);
+        cur = graph
+            .edge_endpoints(next_eid)
+            .expect("next_eid must be valid")
+            .1;
+        cur_in = next_eid;
+    }
+
+    // ── Walk backward from the source node ────────────────────────────────
+    let mut cur = src;
+    let mut cur_out = start_edge;
+    loop {
+        let prev_eid = match backward_continuation(graph, cur, cur_out) {
+            Some(e) => e,
+            None => break,
+        };
+
+        if branch.contains(&prev_eid) || visited.contains(&prev_eid) {
+            break;
+        }
+
+        branch.push(prev_eid);
+        cur = graph
+            .edge_endpoints(prev_eid)
+            .expect("prev_eid must be valid")
+            .0;
+        cur_out = prev_eid;
+    }
+
+    branch
+}
+
+/// Returns the unique outgoing edge continuing the series branch through
+/// `node`, entered via `in_edge` — either the plain 1-in-1-out rule or a
+/// declared internal flow path matching the entry port.
+fn forward_continuation(
+    graph: &Graph<Box<dyn Component>, FlowEdge, Directed>,
+    node: petgraph::graph::NodeIndex,
+    in_edge: EdgeIndex,
+) -> Option<EdgeIndex> {
+    let in_edges: Vec<_> = graph
+        .edges_directed(node, petgraph::Direction::Incoming)
+        .collect();
+    let out_edges: Vec<_> = graph
+        .edges_directed(node, petgraph::Direction::Outgoing)
+        .collect();
+
+    // Plain series node.
+    if in_edges.len() == 1 && out_edges.len() == 1 {
+        return Some(out_edges[0].id());
+    }
+
+    // Multi-port component with declared internal flow paths (Modelica-style):
+    // continue through the path whose inlet port matches the entry port.
+    let entry_port = graph.edge_weight(in_edge).map(|w| w.target_port)?;
+    let paths = graph.node_weight(node)?.flow_paths();
+    let out_port = paths
+        .iter()
+        .find(|(p_in, _)| *p_in == entry_port)
+        .map(|(_, p_out)| *p_out)?;
+    let matching: Vec<EdgeIndex> = out_edges
+        .iter()
+        .filter(|e| e.weight().source_port == out_port)
+        .map(|e| e.id())
+        .collect();
+    match matching.as_slice() {
+        [only] => Some(*only),
+        _ => None,
+    }
+}
+
+/// Backward analogue of [`forward_continuation`]: unique incoming edge whose
+/// target port feeds the flow path that exits through `out_edge`.
+fn backward_continuation(
+    graph: &Graph<Box<dyn Component>, FlowEdge, Directed>,
+    node: petgraph::graph::NodeIndex,
+    out_edge: EdgeIndex,
+) -> Option<EdgeIndex> {
+    let in_edges: Vec<_> = graph
+        .edges_directed(node, petgraph::Direction::Incoming)
+        .collect();
+    let out_edges: Vec<_> = graph
+        .edges_directed(node, petgraph::Direction::Outgoing)
+        .collect();
+
+    if in_edges.len() == 1 && out_edges.len() == 1 {
+        return Some(in_edges[0].id());
+    }
+
+    let exit_port = graph.edge_weight(out_edge).map(|w| w.source_port)?;
+    let paths = graph.node_weight(node)?.flow_paths();
+    let in_port = paths
+        .iter()
+        .find(|(_, p_out)| *p_out == exit_port)
+        .map(|(p_in, _)| *p_in)?;
+    let matching: Vec<EdgeIndex> = in_edges
+        .iter()
+        .filter(|e| e.weight().target_port == in_port)
+        .map(|e| e.id())
+        .collect();
+    match matching.as_slice() {
+        [only] => Some(*only),
+        _ => None,
+    }
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// Unit tests
+// ─────────────────────────────────────────────────────────────────────────────
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use entropyk_components::{ComponentError, JacobianBuilder, ResidualVector, StateSlice};
+    use petgraph::graph::Graph;
+
+    // ── Minimal mock component for topology tests ─────────────────────────
+
+    struct MockComponent;
+
+    impl entropyk_components::Component for MockComponent {
+        fn n_equations(&self) -> usize {
+            2
+        }
+
+        fn compute_residuals(
+            &self,
+            _state: &StateSlice,
+            _residuals: &mut ResidualVector,
+        ) -> Result<(), ComponentError> {
+            Ok(())
+        }
+
+        fn jacobian_entries(
+            &self,
+            _state: &StateSlice,
+            _jacobian: &mut JacobianBuilder,
+        ) -> Result<(), ComponentError> {
+            Ok(())
+        }
+
+        fn get_ports(&self) -> &[entropyk_components::ConnectedPort] {
+            &[]
+        }
+
+        fn signature(&self) -> String {
+            "mock".to_string()
+        }
+    }
+
+    /// Build a directed graph with `n_nodes` arranged as a directed cycle:
+    /// node 0 → node 1 → ... → node n-1 → node 0.
+    /// Returns (graph, edge_indices_in_order).
+    fn build_series_cycle(
+        n: usize,
+    ) -> (
+        Graph<Box<dyn Component>, FlowEdge, Directed>,
+        Vec<EdgeIndex>,
+    ) {
+        let mut g: Graph<Box<dyn Component>, FlowEdge, Directed> = Graph::new();
+        let nodes: Vec<_> = (0..n)
+            .map(|_| g.add_node(Box::new(MockComponent) as Box<dyn Component>))
+            .collect();
+        let mut edges = Vec::new();
+        for i in 0..n {
+            let src = nodes[i];
+            let tgt = nodes[(i + 1) % n];
+            edges.push(g.add_edge(src, tgt, FlowEdge::new_unassigned()));
+        }
+        (g, edges)
+    }
+
+    /// Build a graph that represents a main loop (4 edges) plus a bypass
+    /// branch: node 0 → (splitter) node 1 → two paths → (merger) node 4 → node 0.
+    /// Topology: 0→1→2→4→0 (main) and 1→3→4 (bypass), so node 1 has out-degree 2
+    /// and node 4 has in-degree 2.
+    fn build_splitter_merger_topology() -> (
+        Graph<Box<dyn Component>, FlowEdge, Directed>,
+        Vec<EdgeIndex>,
+    ) {
+        // Nodes: 0 (source/sink), 1 (splitter), 2 (main path), 3 (bypass path), 4 (merger)
+        let mut g: Graph<Box<dyn Component>, FlowEdge, Directed> = Graph::new();
+        let n: Vec<_> = (0..5)
+            .map(|_| g.add_node(Box::new(MockComponent) as Box<dyn Component>))
+            .collect();
+
+        // Main series edges: 0→1, 2→4, 4→0 (these are outside the junction nodes)
+        let e01 = g.add_edge(n[0], n[1], FlowEdge::new_unassigned()); // branch A
+        let e12 = g.add_edge(n[1], n[2], FlowEdge::new_unassigned()); // branch B (main)
+        let e24 = g.add_edge(n[2], n[4], FlowEdge::new_unassigned()); // branch B (main)
+        let e13 = g.add_edge(n[1], n[3], FlowEdge::new_unassigned()); // branch C (bypass)
+        let e34 = g.add_edge(n[3], n[4], FlowEdge::new_unassigned()); // branch C (bypass)
+        let e40 = g.add_edge(n[4], n[0], FlowEdge::new_unassigned()); // branch A
+
+        (g, vec![e01, e12, e24, e13, e34, e40])
+    }
+
+    // ── Test: 4-edge series cycle → 1 branch ─────────────────────────────
+
+    #[test]
+    fn test_series_cycle_4_edges_one_branch() {
+        let (mut g, edges) = build_series_cycle(4);
+
+        let n_branches = presolve_mass_flow_topology(&mut g);
+
+        assert_eq!(
+            n_branches, 1,
+            "pure series cycle must yield exactly 1 branch"
+        );
+
+        // All 4 edges must share branch_id == 0.
+        let branch_ids: Vec<usize> = edges
+            .iter()
+            .map(|&e| g.edge_weight(e).unwrap().mass_flow_branch_id)
+            .collect();
+        assert!(
+            branch_ids.iter().all(|&id| id == 0),
+            "all edges must be in branch 0, got {:?}",
+            branch_ids
+        );
+    }
+
+    // ── Test: 2-edge series → 1 branch ───────────────────────────────────
+
+    #[test]
+    fn test_series_cycle_2_edges_one_branch() {
+        let (mut g, edges) = build_series_cycle(2);
+        let n_branches = presolve_mass_flow_topology(&mut g);
+        assert_eq!(n_branches, 1);
+        let ids: Vec<usize> = edges
+            .iter()
+            .map(|&e| g.edge_weight(e).unwrap().mass_flow_branch_id)
+            .collect();
+        assert!(ids.iter().all(|&id| id == 0));
+    }
+
+    // ── Test: single-node self-loop (degenerate) — 1 branch ──────────────
+
+    #[test]
+    fn test_single_edge_one_branch() {
+        let (mut g, edges) = build_series_cycle(1);
+        let n_branches = presolve_mass_flow_topology(&mut g);
+        assert_eq!(n_branches, 1);
+        assert_eq!(g.edge_weight(edges[0]).unwrap().mass_flow_branch_id, 0);
+    }
+
+    // ── Test: splitter-merger topology → 3 branches ───────────────────────
+
+    #[test]
+    fn test_splitter_merger_three_branches() {
+        // Topology:
+        //   node 0 →[e01]→ node 1 (splitter, out-degree 2)
+        //   node 1 →[e12]→ node 2 →[e24]→ node 4 (merger, in-degree 2)
+        //   node 1 →[e13]→ node 3 →[e34]→ node 4
+        //   node 4 →[e40]→ node 0
+        //
+        // Expected branches:
+        //   Branch A: e01, e40  (series path from node 0 through junction node 1 back via node 4)
+        //   Wait — node 0 has in-degree 1 (from e40) and out-degree 1 (to e01), so it IS a 1-in-1-out node.
+        //   node 1 has out-degree 2 → junction boundary.
+        //   node 4 has in-degree 2 → junction boundary.
+        //
+        //   So tracing from e01:
+        //     start: e01 (0→1)
+        //     forward from node 1: out-degree=2 → STOP
+        //     backward from node 0: in-edges = [e40], out-edges = [e01] → 1-in-1-out → add e40, move to node 4
+        //     node 4: in-degree=2 → STOP
+        //   Branch A: {e01, e40}
+        //
+        //   From e12 (unvisited):
+        //     forward from node 2: in-degree=1, out-degree=1 → add e24, move to node 4; node 4 in-degree=2 → STOP
+        //     backward from node 1: out-degree=2 → STOP
+        //   Branch B: {e12, e24}
+        //
+        //   From e13 (unvisited):
+        //     forward from node 3: in-degree=1, out-degree=1 → add e34, move to node 4; node 4 in-degree=2 → STOP
+        //     backward from node 1: out-degree=2 → STOP
+        //   Branch C: {e13, e34}
+        //
+        //   Total: 3 branches ✓
+
+        let (mut g, edges) = build_splitter_merger_topology();
+        let [e01, e12, e24, e13, e34, e40] = edges.as_slice() else {
+            panic!()
+        };
+
+        let n_branches = presolve_mass_flow_topology(&mut g);
+
+        assert_eq!(n_branches, 3, "splitter-merger must yield 3 branches");
+
+        let id = |e: &EdgeIndex| g.edge_weight(*e).unwrap().mass_flow_branch_id;
+
+        // e01 and e40 must share a branch.
+        assert_eq!(id(e01), id(e40), "e01 and e40 must share branch");
+
+        // e12 and e24 must share a branch (main path through node 2).
+        assert_eq!(id(e12), id(e24), "e12 and e24 must share branch");
+
+        // e13 and e34 must share a branch (bypass through node 3).
+        assert_eq!(id(e13), id(e34), "e13 and e34 must share branch");
+
+        // All three pairs must be in DISTINCT branches.
+        let a = id(e01);
+        let b = id(e12);
+        let c = id(e13);
+        assert_ne!(a, b, "branch A and B must differ");
+        assert_ne!(a, c, "branch A and C must differ");
+        assert_ne!(b, c, "branch B and C must differ");
+    }
+
+    // ── Test: no double-assignment on a 6-edge series ─────────────────────
+
+    #[test]
+    fn test_no_double_assignment_large_cycle() {
+        let (mut g, edges) = build_series_cycle(6);
+        let n_branches = presolve_mass_flow_topology(&mut g);
+        assert_eq!(n_branches, 1, "6-edge series cycle is still 1 branch");
+        for &e in &edges {
+            assert_eq!(g.edge_weight(e).unwrap().mass_flow_branch_id, 0);
+        }
+    }
+}