Stabilize Modelica Fixed/Free calibration and stop flaky Newton embeds.

Sort constraint/control assembly for deterministic Jacobians, keep live z_ua on legacy HX, and add generic Fixed/Free UI assist without removing Z factors.

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

crates/solver/src/inverse/bounded.rs

@@ -42,7 +42,7 @@ use thiserror::Error;
 /// Type-safe identifier for a bounded control variable.
 ///
 /// Uses a string internally but provides type safety and clear intent.
-#[derive(Debug, Clone, PartialEq, Eq, Hash)]
+#[derive(Debug, Clone, PartialEq, Eq, PartialOrd, Ord, Hash)]
 pub struct BoundedVariableId(String);
 
 impl BoundedVariableId {

crates/solver/src/inverse/constraint.rs

@@ -16,7 +16,7 @@ use thiserror::Error;
 /// Type-safe identifier for a constraint.
 ///
 /// Uses a string internally but provides type safety and clear intent.
-#[derive(Debug, Clone, PartialEq, Eq, Hash)]
+#[derive(Debug, Clone, PartialEq, Eq, PartialOrd, Ord, Hash)]
 pub struct ConstraintId(String);
 
 impl ConstraintId {

crates/solver/src/inverse/embedding.rs

@@ -334,19 +334,24 @@ impl InverseControlConfig {
         self.control_to_constraint.get(bounded_variable_id)
     }
 
-    /// Returns an iterator over all constraint-to-control mappings.
+    /// Returns constraint→control mappings in **stable** order (sorted by constraint id).
+    ///
+    /// HashMap iteration is randomized per process; residual rows, Jacobian rows/columns,
+    /// and control state indices must share one deterministic order or Newton flakes.
     pub fn mappings(&self) -> impl Iterator<Item = (&ConstraintId, &BoundedVariableId)> {
-        self.constraint_to_control.iter()
+        let mut items: Vec<_> = self.constraint_to_control.iter().collect();
+        items.sort_by(|a, b| a.0.cmp(b.0));
+        items.into_iter()
     }
 
-    /// Returns an iterator over linked constraint IDs.
+    /// Linked constraint IDs in the same order as [`Self::mappings`].
     pub fn linked_constraints(&self) -> impl Iterator<Item = &ConstraintId> {
-        self.constraint_to_control.keys()
+        self.mappings().map(|(c, _)| c)
     }
 
-    /// Returns an iterator over linked control variable IDs.
+    /// Linked control IDs in the same order as [`Self::mappings`] (column / state layout).
     pub fn linked_controls(&self) -> impl Iterator<Item = &BoundedVariableId> {
-        self.control_to_constraint.keys()
+        self.mappings().map(|(_, c)| c)
     }
 
     /// Checks if a constraint is linked.
@@ -583,15 +588,18 @@ mod tests {
     fn test_inverse_control_config_mappings_iterator() {
         let mut config = InverseControlConfig::new();
 
+        // Insert out of alphabetical order — iteration must still be sorted.
         config
-            .link(make_constraint_id("c1"), make_bounded_var_id("v1"))
+            .link(make_constraint_id("z_probe"), make_bounded_var_id("z_ua_cond"))
             .unwrap();
         config
-            .link(make_constraint_id("c2"), make_bounded_var_id("v2"))
+            .link(make_constraint_id("a_probe"), make_bounded_var_id("z_ua_evap"))
             .unwrap();
 
-        let mappings: Vec<_> = config.mappings().collect();
-        assert_eq!(mappings.len(), 2);
+        let ids: Vec<&str> = config.mappings().map(|(c, _)| c.as_str()).collect();
+        assert_eq!(ids, vec!["a_probe", "z_probe"]);
+        let controls: Vec<&str> = config.linked_controls().map(|c| c.as_str()).collect();
+        assert_eq!(controls, vec!["z_ua_evap", "z_ua_cond"]);
     }
 
     #[test]

crates/solver/src/inverse/expr.rs

@@ -0,0 +1,401 @@
+//! Minimal Modelica-style algebraic expression AST for model equations.
+//!
+//! Supports: `+ - * /`, parentheses, unary `-`, `min`, `max`, `exp`, `ln`/`log`,
+//! numeric literals, and references `component.field` or `"Component Name".field`.
+
+use std::fmt;
+
+/// Parsed expression node.
+#[derive(Debug, Clone, PartialEq)]
+pub enum Expr {
+    Const(f64),
+    /// `component` + `field` (field is temperature, z_ua, …).
+    Ref { component: String, field: String },
+    Neg(Box<Expr>),
+    Add(Box<Expr>, Box<Expr>),
+    Sub(Box<Expr>, Box<Expr>),
+    Mul(Box<Expr>, Box<Expr>),
+    Div(Box<Expr>, Box<Expr>),
+    Min(Box<Expr>, Box<Expr>),
+    Max(Box<Expr>, Box<Expr>),
+    Exp(Box<Expr>),
+    Ln(Box<Expr>),
+}
+
+impl Expr {
+    /// Evaluate with a resolver for `component.field` references.
+    pub fn eval<F>(&self, resolve: &mut F) -> Result<f64, String>
+    where
+        F: FnMut(&str, &str) -> Result<f64, String>,
+    {
+        match self {
+            Expr::Const(v) => Ok(*v),
+            Expr::Ref { component, field } => resolve(component, field),
+            Expr::Neg(a) => Ok(-a.eval(resolve)?),
+            Expr::Add(a, b) => Ok(a.eval(resolve)? + b.eval(resolve)?),
+            Expr::Sub(a, b) => Ok(a.eval(resolve)? - b.eval(resolve)?),
+            Expr::Mul(a, b) => Ok(a.eval(resolve)? * b.eval(resolve)?),
+            Expr::Div(a, b) => {
+                let den = b.eval(resolve)?;
+                if den.abs() < 1e-30 {
+                    return Err("division by zero in model equation".into());
+                }
+                Ok(a.eval(resolve)? / den)
+            }
+            Expr::Min(a, b) => Ok(a.eval(resolve)?.min(b.eval(resolve)?)),
+            Expr::Max(a, b) => Ok(a.eval(resolve)?.max(b.eval(resolve)?)),
+            Expr::Exp(a) => {
+                let x = a.eval(resolve)?;
+                let y = x.exp();
+                if y.is_finite() {
+                    Ok(y)
+                } else {
+                    Err(format!("exp({x}) overflow"))
+                }
+            }
+            Expr::Ln(a) => {
+                let x = a.eval(resolve)?;
+                if x <= 0.0 {
+                    return Err(format!("ln({x}) domain error"));
+                }
+                Ok(x.ln())
+            }
+        }
+    }
+
+    /// Collect all `component.field` references.
+    pub fn collect_refs(&self, out: &mut Vec<(String, String)>) {
+        match self {
+            Expr::Const(_) => {}
+            Expr::Ref { component, field } => out.push((component.clone(), field.clone())),
+            Expr::Neg(a) | Expr::Exp(a) | Expr::Ln(a) => a.collect_refs(out),
+            Expr::Add(a, b)
+            | Expr::Sub(a, b)
+            | Expr::Mul(a, b)
+            | Expr::Div(a, b)
+            | Expr::Min(a, b)
+            | Expr::Max(a, b) => {
+                a.collect_refs(out);
+                b.collect_refs(out);
+            }
+        }
+    }
+}
+
+/// Parse an expression string into an [`Expr`].
+pub fn parse_expr(input: &str) -> Result<Expr, String> {
+    let mut p = Parser::new(input);
+    let e = p.parse_expr()?;
+    p.skip_ws();
+    if p.pos < p.chars.len() {
+        return Err(format!(
+            "trailing junk at '{}'",
+            p.chars[p.pos..].iter().collect::<String>()
+        ));
+    }
+    Ok(e)
+}
+
+struct Parser {
+    chars: Vec<char>,
+    pos: usize,
+}
+
+impl Parser {
+    fn new(s: &str) -> Self {
+        Self {
+            chars: s.chars().collect(),
+            pos: 0,
+        }
+    }
+
+    fn skip_ws(&mut self) {
+        while self.pos < self.chars.len() && self.chars[self.pos].is_whitespace() {
+            self.pos += 1;
+        }
+    }
+
+    fn peek(&self) -> Option<char> {
+        self.chars.get(self.pos).copied()
+    }
+
+    fn bump(&mut self) -> Option<char> {
+        let c = self.peek()?;
+        self.pos += 1;
+        Some(c)
+    }
+
+    fn parse_expr(&mut self) -> Result<Expr, String> {
+        let mut left = self.parse_term()?;
+        loop {
+            self.skip_ws();
+            match self.peek() {
+                Some('+') => {
+                    self.bump();
+                    let right = self.parse_term()?;
+                    left = Expr::Add(Box::new(left), Box::new(right));
+                }
+                Some('-') => {
+                    self.bump();
+                    let right = self.parse_term()?;
+                    left = Expr::Sub(Box::new(left), Box::new(right));
+                }
+                _ => break,
+            }
+        }
+        Ok(left)
+    }
+
+    fn parse_term(&mut self) -> Result<Expr, String> {
+        let mut left = self.parse_unary()?;
+        loop {
+            self.skip_ws();
+            match self.peek() {
+                Some('*') => {
+                    self.bump();
+                    let right = self.parse_unary()?;
+                    left = Expr::Mul(Box::new(left), Box::new(right));
+                }
+                Some('/') => {
+                    self.bump();
+                    let right = self.parse_unary()?;
+                    left = Expr::Div(Box::new(left), Box::new(right));
+                }
+                _ => break,
+            }
+        }
+        Ok(left)
+    }
+
+    fn parse_unary(&mut self) -> Result<Expr, String> {
+        self.skip_ws();
+        if self.peek() == Some('-') {
+            self.bump();
+            let e = self.parse_unary()?;
+            return Ok(Expr::Neg(Box::new(e)));
+        }
+        if self.peek() == Some('+') {
+            self.bump();
+            return self.parse_unary();
+        }
+        self.parse_primary()
+    }
+
+    fn parse_primary(&mut self) -> Result<Expr, String> {
+        self.skip_ws();
+        match self.peek() {
+            Some('(') => {
+                self.bump();
+                let e = self.parse_expr()?;
+                self.skip_ws();
+                if self.bump() != Some(')') {
+                    return Err("expected ')'".into());
+                }
+                Ok(e)
+            }
+            Some('"') => self.parse_quoted_ref(),
+            Some(c) if c.is_ascii_digit() || c == '.' => self.parse_number(),
+            Some(c) if is_ident_start(c) => self.parse_ident_or_call(),
+            Some(c) => Err(format!("unexpected '{c}'")),
+            None => Err("unexpected end of expression".into()),
+        }
+    }
+
+    fn parse_number(&mut self) -> Result<Expr, String> {
+        let start = self.pos;
+        while let Some(c) = self.peek() {
+            if c.is_ascii_digit() || c == '.' || c == 'e' || c == 'E' {
+                self.bump();
+                if (c == 'e' || c == 'E') && matches!(self.peek(), Some('+') | Some('-')) {
+                    self.bump();
+                }
+            } else {
+                break;
+            }
+        }
+        let s: String = self.chars[start..self.pos].iter().collect();
+        s.parse::<f64>()
+            .map(Expr::Const)
+            .map_err(|_| format!("invalid number '{s}'"))
+    }
+
+    fn parse_quoted_ref(&mut self) -> Result<Expr, String> {
+        self.bump(); // "
+        let start = self.pos;
+        while let Some(c) = self.peek() {
+            if c == '"' {
+                break;
+            }
+            self.bump();
+        }
+        let name: String = self.chars[start..self.pos].iter().collect();
+        if self.bump() != Some('"') {
+            return Err("unterminated quoted component name".into());
+        }
+        self.skip_ws();
+        if self.bump() != Some('.') {
+            return Err("expected '.' after quoted component name".into());
+        }
+        let field = self.parse_ident_string()?;
+        Ok(Expr::Ref {
+            component: name,
+            field,
+        })
+    }
+
+    fn parse_ident_string(&mut self) -> Result<String, String> {
+        self.skip_ws();
+        let start = self.pos;
+        if !self.peek().is_some_and(is_ident_start) {
+            return Err("expected identifier".into());
+        }
+        self.bump();
+        while self.peek().is_some_and(is_ident_cont) {
+            self.bump();
+        }
+        Ok(self.chars[start..self.pos].iter().collect())
+    }
+
+    fn parse_ident_or_call(&mut self) -> Result<Expr, String> {
+        let ident = self.parse_ident_string()?;
+        self.skip_ws();
+        if self.peek() == Some('(') {
+            return self.parse_call(&ident);
+        }
+        // component.field[.more] — last segment is field, rest is component
+        let mut parts = vec![ident];
+        while self.peek() == Some('.') {
+            self.bump();
+            parts.push(self.parse_ident_string()?);
+            self.skip_ws();
+        }
+        if parts.len() < 2 {
+            return Err(format!(
+                "bare identifier '{}' — expected component.field or function call",
+                parts[0]
+            ));
+        }
+        let field = parts.pop().unwrap();
+        let component = parts.join(".");
+        Ok(Expr::Ref { component, field })
+    }
+
+    fn parse_call(&mut self, name: &str) -> Result<Expr, String> {
+        self.bump(); // (
+        self.skip_ws();
+        let lower = name.to_ascii_lowercase();
+        match lower.as_str() {
+            "min" | "max" => {
+                let a = self.parse_expr()?;
+                self.skip_ws();
+                if self.bump() != Some(',') {
+                    return Err(format!("expected ',' in {lower}(...)"));
+                }
+                let b = self.parse_expr()?;
+                self.skip_ws();
+                if self.bump() != Some(')') {
+                    return Err(format!("expected ')' after {lower}"));
+                }
+                Ok(if lower == "min" {
+                    Expr::Min(Box::new(a), Box::new(b))
+                } else {
+                    Expr::Max(Box::new(a), Box::new(b))
+                })
+            }
+            "exp" => {
+                let a = self.parse_expr()?;
+                self.skip_ws();
+                if self.bump() != Some(')') {
+                    return Err("expected ')' after exp".into());
+                }
+                Ok(Expr::Exp(Box::new(a)))
+            }
+            "ln" | "log" => {
+                let a = self.parse_expr()?;
+                self.skip_ws();
+                if self.bump() != Some(')') {
+                    return Err("expected ')' after ln".into());
+                }
+                Ok(Expr::Ln(Box::new(a)))
+            }
+            other => Err(format!(
+                "unknown function '{other}' (supported: min, max, exp, ln)"
+            )),
+        }
+    }
+}
+
+fn is_ident_start(c: char) -> bool {
+    c.is_ascii_alphabetic() || c == '_'
+}
+
+fn is_ident_cont(c: char) -> bool {
+    c.is_ascii_alphanumeric() || c == '_'
+}
+
+impl fmt::Display for Expr {
+    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
+        match self {
+            Expr::Const(v) => write!(f, "{v}"),
+            Expr::Ref { component, field } => write!(f, "{component}.{field}"),
+            Expr::Neg(a) => write!(f, "-({a})"),
+            Expr::Add(a, b) => write!(f, "({a} + {b})"),
+            Expr::Sub(a, b) => write!(f, "({a} - {b})"),
+            Expr::Mul(a, b) => write!(f, "({a} * {b})"),
+            Expr::Div(a, b) => write!(f, "({a} / {b})"),
+            Expr::Min(a, b) => write!(f, "min({a}, {b})"),
+            Expr::Max(a, b) => write!(f, "max({a}, {b})"),
+            Expr::Exp(a) => write!(f, "exp({a})"),
+            Expr::Ln(a) => write!(f, "ln({a})"),
+        }
+    }
+}
+
+/// Algebraic model equation: `lhs - rhs = 0`.
+#[derive(Debug, Clone)]
+pub struct AlgebraicEquation {
+    pub id: String,
+    pub lhs: Expr,
+    pub rhs: Expr,
+}
+
+impl AlgebraicEquation {
+    pub fn residual<F>(&self, mut resolve: F) -> Result<f64, String>
+    where
+        F: FnMut(&str, &str) -> Result<f64, String>,
+    {
+        Ok(self.lhs.eval(&mut resolve)? - self.rhs.eval(&mut resolve)?)
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    fn eval_const(e: &Expr) -> f64 {
+        e.eval(&mut |_, _| Err("no refs".into())).unwrap()
+    }
+
+    #[test]
+    fn arithmetic_and_funcs() {
+        assert!((eval_const(&parse_expr("1 + 2 * 3").unwrap()) - 7.0).abs() < 1e-12);
+        assert!((eval_const(&parse_expr("min(2, 5)").unwrap()) - 2.0).abs() < 1e-12);
+        assert!((eval_const(&parse_expr("max(2, 5)").unwrap()) - 5.0).abs() < 1e-12);
+        assert!((eval_const(&parse_expr("exp(0)").unwrap()) - 1.0).abs() < 1e-12);
+        assert!((eval_const(&parse_expr("ln(exp(1))").unwrap()) - 1.0).abs() < 1e-9);
+        assert!((eval_const(&parse_expr("-(1+2)").unwrap()) + 3.0).abs() < 1e-12);
+    }
+
+    #[test]
+    fn refs_and_quoted() {
+        let e = parse_expr(r#"evap.z_ua + "SST probe".T"#).unwrap();
+        let v = e
+            .eval(&mut |c, f| match (c, f) {
+                ("evap", "z_ua") => Ok(0.1),
+                ("SST probe", "T") => Ok(280.0),
+                _ => Err("bad".into()),
+            })
+            .unwrap();
+        assert!((v - 280.1).abs() < 1e-12);
+    }
+}

crates/solver/src/inverse/mod.rs

@@ -45,6 +45,7 @@ pub mod bounded;
 pub mod calibration;
 pub mod constraint;
 pub mod embedding;
+pub mod expr;
 pub mod override_network;
 pub mod saturated_control;
 
@@ -58,5 +59,6 @@ pub use calibration::{
 };
 pub use constraint::{ComponentOutput, Constraint, ConstraintError, ConstraintId};
 pub use embedding::{ControlMapping, DoFError, InverseControlConfig};
+pub use expr::{parse_expr, AlgebraicEquation, Expr};
 pub use override_network::{eval_error_signal, eval_error_weights, Combine, Objective};
 pub use saturated_control::{SaturatedControlError, SaturatedController, Saturation};

crates/solver/src/system.rs

@@ -16,7 +16,7 @@ use petgraph::algo;
 use petgraph::graph::{EdgeIndex, Graph, NodeIndex};
 use petgraph::visit::EdgeRef;
 use petgraph::Directed;
-use std::collections::HashMap;
+use std::collections::{BTreeSet, HashMap};
 
 use crate::coupling::{has_circular_dependencies, ThermalCoupling};
 use crate::dof::{
@@ -208,6 +208,8 @@ pub struct System {
     thermal_couplings: Vec<ThermalCoupling>,
     /// Constraints for inverse control (output - target = 0)
     constraints: HashMap<ConstraintId, Constraint>,
+    /// Modelica algebraic equations (`lhs - rhs = 0`), e.g. `evap.z_ua = f(...)`.
+    algebraic_equations: Vec<crate::inverse::AlgebraicEquation>,
     /// Bounded control variables for inverse control (with box constraints)
     bounded_variables: HashMap<BoundedVariableId, BoundedVariable>,
     /// Inverse control configuration (constraint → control variable mappings)
@@ -246,6 +248,7 @@ impl System {
             node_to_circuit: HashMap::new(),
             thermal_couplings: Vec::new(),
             constraints: HashMap::new(),
+            algebraic_equations: Vec::new(),
             bounded_variables: HashMap::new(),
             inverse_control: InverseControlConfig::new(),
             saturated_controllers: Vec::new(),
@@ -735,18 +738,20 @@ impl System {
             }
         }
 
-        if !self.constraints.is_empty() {
+        if !self.constraints.is_empty() || !self.algebraic_equations.is_empty() {
             match self.validate_inverse_control_dof() {
                 Ok(()) => {
                     tracing::debug!(
-                        constraint_count = self.constraints.len(),
+                        constraint_count =
+                            self.constraints.len() + self.algebraic_equations.len(),
                         control_count = self.inverse_control.mapping_count(),
                         "Inverse control DoF validation passed"
                     );
                 }
                 Err(DoFError::UnderConstrainedSystem { .. }) => {
                     tracing::warn!(
-                        constraint_count = self.constraints.len(),
+                        constraint_count =
+                            self.constraints.len() + self.algebraic_equations.len(),
                         control_count = self.inverse_control.mapping_count(),
                         "Under-constrained inverse control system - solver may still converge"
                     );
@@ -1240,6 +1245,21 @@ impl System {
         Ok(())
     }
 
+    /// Adds a Modelica algebraic equation residual (`lhs − rhs = 0`).
+    pub fn add_algebraic_equation(
+        &mut self,
+        equation: crate::inverse::AlgebraicEquation,
+    ) -> Result<(), ConstraintError> {
+        let id = ConstraintId::new(equation.id.clone());
+        if self.constraints.contains_key(&id)
+            || self.algebraic_equations.iter().any(|e| e.id == equation.id)
+        {
+            return Err(ConstraintError::DuplicateId { id });
+        }
+        self.algebraic_equations.push(equation);
+        Ok(())
+    }
+
     /// Removes a constraint by ID.
     ///
     /// # Arguments
@@ -1259,9 +1279,16 @@ impl System {
         self.constraints.len()
     }
 
-    /// Returns a reference to all constraints.
+    /// Number of Modelica algebraic equations (`lhs − rhs = 0`).
+    pub fn algebraic_equation_count(&self) -> usize {
+        self.algebraic_equations.len()
+    }
+
+    /// Returns all constraints in **stable** order (sorted by constraint id).
     pub fn constraints(&self) -> impl Iterator<Item = &Constraint> {
-        self.constraints.values()
+        let mut items: Vec<_> = self.constraints.values().collect();
+        items.sort_by(|a, b| a.id().cmp(b.id()));
+        items.into_iter()
     }
 
     /// Returns a reference to a specific constraint by ID.
@@ -1317,8 +1344,12 @@ impl System {
             return Ok(0);
         }
 
+        // Must match Jacobian row order (`inverse_control.mappings()` / sorted ids).
+        let mut sorted: Vec<_> = self.constraints.values().collect();
+        sorted.sort_by(|a, b| a.id().cmp(b.id()));
+
         let mut count = 0;
-        for constraint in self.constraints.values() {
+        for constraint in sorted {
             let measured = match measured_values.get(constraint.id()).copied() {
                 Some(v) => v,
                 None => {
@@ -2076,7 +2107,12 @@ impl System {
         constraint_id: &ConstraintId,
         bounded_variable_id: &BoundedVariableId,
     ) -> Result<(), DoFError> {
-        if !self.constraints.contains_key(constraint_id) {
+        let has_constraint = self.constraints.contains_key(constraint_id);
+        let has_algebraic = self
+            .algebraic_equations
+            .iter()
+            .any(|e| e.id == constraint_id.as_str());
+        if !has_constraint && !has_algebraic {
             return Err(DoFError::ConstraintNotFound {
                 constraint_id: constraint_id.clone(),
             });
@@ -2191,7 +2227,7 @@ impl System {
     pub fn validate_inverse_control_dof(&self) -> Result<(), DoFError> {
         let n_edge_unknowns = self.total_state_len;
         let n_controls = self.inverse_control.mapping_count();
-        let n_constraints = self.constraints.len();
+        let n_constraints = self.constraints.len() + self.algebraic_equations.len();
         let n_unknowns = n_edge_unknowns + n_controls;
 
         let n_edge_eqs: usize = self
@@ -2242,6 +2278,7 @@ impl System {
             .sum();
         n_comp
             + self.constraints.len()
+            + self.algebraic_equations.len()
             + self.coupling_residual_count()
             + 2 * self.saturated_controllers.len()
     }
@@ -2319,6 +2356,11 @@ impl System {
                 name: id.to_string(),
             });
         }
+        for eq in &self.algebraic_equations {
+            system_equations.push(EquationRole::ControlTracking {
+                name: eq.id.clone(),
+            });
+        }
         for _ in 0..self.coupling_residual_count() {
             system_equations.push(EquationRole::CouplingDuty);
         }
@@ -2701,6 +2743,7 @@ impl System {
             .map(|(_, c, _)| c.n_equations())
             .sum();
         total_eqs += self.constraints.len()
+            + self.algebraic_equations.len()
             + self.coupling_residual_count()
             + 2 * self.saturated_controllers.len();
 
@@ -2734,6 +2777,15 @@ impl System {
             .map_err(|e| ComponentError::CalculationFailed(e.to_string()))?;
         eq_offset += n_constraints;
 
+        // Modelica algebraic equations: lhs − rhs = 0
+        for (i, eq) in self.algebraic_equations.iter().enumerate() {
+            let r = eq
+                .residual(|comp, field| self.resolve_expr_ref(comp, field, state))
+                .map_err(ComponentError::CalculationFailed)?;
+            residuals[eq_offset + i] = r;
+        }
+        eq_offset += self.algebraic_equations.len();
+
         // Add couplings. Each coupling owns one unknown Q = state[coupling_state_index(i)].
         //
         // * Physical mode (`hot_component`/`cold_component` set): the residual
@@ -2836,12 +2888,19 @@ impl System {
             jacobian.add_entry(r, c, v);
         }
 
+        let algebraic_row_offset = row_offset + self.constraints.len();
+        let algebraic_jac = self.compute_algebraic_equation_jacobian(state, algebraic_row_offset);
+        for (r, c, v) in algebraic_jac {
+            jacobian.add_entry(r, c, v);
+        }
+
         // Thermal coupling rows: r_i = Q_i − η·duty_hot(state).
         // ∂r/∂Q = 1 always (also fixes the legacy stub's singular row); in
         // physical mode the plant coupling ∂r/∂col = −η·∂duty/∂col is formed by
         // central finite differences over the hot component's incident states
         // (same pattern as the saturated-controller measurement Jacobian).
-        let coupling_row_offset = row_offset + self.constraints.len();
+        let coupling_row_offset =
+            row_offset + self.constraints.len() + self.algebraic_equations.len();
         if !self.thermal_couplings.is_empty() {
             let eps = self.inverse_control.finite_diff_epsilon();
             let mut state_mut = state.to_vec();
@@ -2882,8 +2941,10 @@ impl System {
             }
         }
 
-        let saturated_row_offset =
-            row_offset + self.constraints.len() + self.coupling_residual_count();
+        let saturated_row_offset = row_offset
+            + self.constraints.len()
+            + self.algebraic_equations.len()
+            + self.coupling_residual_count();
         let saturated_jac = self.compute_saturated_control_jacobian(state, saturated_row_offset);
         for (r, c, v) in saturated_jac {
             jacobian.add_entry(r, c, v);
@@ -2891,6 +2952,147 @@ impl System {
         Ok(())
     }
 
+    /// Resolve `component.field` for algebraic model equations.
+    ///
+    /// Fields may be Z-factors (`z_ua`, …) or measured outputs (`T`, `Tsat`, …).
+    pub fn resolve_expr_ref(
+        &self,
+        component: &str,
+        field: &str,
+        state: &StateSlice,
+    ) -> Result<f64, String> {
+        if let Some(canon) = entropyk_core::normalize_factor_name(field) {
+            let indices = self.calib_indices_by_name.get(component).ok_or_else(|| {
+                format!("unknown component '{component}' for Z-factor '{field}'")
+            })?;
+            let idx = match canon {
+                entropyk_core::Z_FLOW => indices.z_flow,
+                entropyk_core::Z_FLOW_ECO => indices.z_flow_eco,
+                entropyk_core::Z_DP => indices.z_dp,
+                entropyk_core::Z_UA => indices.z_ua,
+                entropyk_core::Z_POWER => indices.z_power,
+                entropyk_core::Z_ETAV => indices.z_etav,
+                entropyk_core::F_W => indices.f_w,
+                _ => None,
+            }
+            .ok_or_else(|| {
+                format!("Z-factor '{field}' is not a free unknown on component '{component}'")
+            })?;
+            return state.get(idx).copied().ok_or_else(|| {
+                format!("state index {idx} out of range for {component}.{field}")
+            });
+        }
+
+        let kind = Self::measured_output_from_field(field).ok_or_else(|| {
+            format!(
+                "unknown field '{field}' on '{component}' \
+                 (expected z_*, T/temperature, Tsat/SST/SDT, P/pressure, …)"
+            )
+        })?;
+        let &node_idx = self.component_names.get(component).ok_or_else(|| {
+            format!("unknown component '{component}'")
+        })?;
+        let component_ref = self.graph.node_weight(node_idx).ok_or_else(|| {
+            format!("component '{component}' missing from graph")
+        })?;
+        component_ref
+            .measure_output(kind, state)
+            .filter(|v| v.is_finite())
+            .ok_or_else(|| {
+                format!("measure '{field}' unavailable on component '{component}'")
+            })
+    }
+
+    fn measured_output_from_field(field: &str) -> Option<entropyk_components::MeasuredOutput> {
+        use entropyk_components::MeasuredOutput as MO;
+        let n = field.trim().to_ascii_lowercase().replace(['_', '-'], "");
+        Some(match n.as_str() {
+            "t" | "temperature" => MO::Temperature,
+            "tsat" | "sst" | "sdt" | "saturationtemperature" => MO::SaturationTemperature,
+            "p" | "pressure" => MO::Pressure,
+            "m" | "mdot" | "massflow" | "massflowrate" => MO::MassFlowRate,
+            "q" | "capacity" => MO::Capacity,
+            "duty" | "heattransferrate" | "heatrate" => MO::HeatTransferRate,
+            "superheat" | "tsh" | "sh" => MO::Superheat,
+            "subcooling" | "tsc" | "sc" => MO::Subcooling,
+            _ => return None,
+        })
+    }
+
+    fn compute_algebraic_equation_jacobian(
+        &self,
+        state: &StateSlice,
+        row_offset: usize,
+    ) -> Vec<(usize, usize, f64)> {
+        let mut entries = Vec::new();
+        if self.algebraic_equations.is_empty() {
+            return entries;
+        }
+        let eps = self.inverse_control.finite_diff_epsilon();
+        if state.len() < self.full_state_vector_len() {
+            tracing::error!(
+                state_len = state.len(),
+                required = self.full_state_vector_len(),
+                "compute_algebraic_equation_jacobian: state slice too short"
+            );
+            return entries;
+        }
+        let mut state_mut = state.to_vec();
+
+        for (i, eq) in self.algebraic_equations.iter().enumerate() {
+            let row = row_offset + i;
+            let mut refs = Vec::new();
+            eq.lhs.collect_refs(&mut refs);
+            eq.rhs.collect_refs(&mut refs);
+
+            let mut cols = BTreeSet::new();
+            for (comp, field) in &refs {
+                if let Some(canon) = entropyk_core::normalize_factor_name(field) {
+                    if let Some(indices) = self.calib_indices_by_name.get(comp) {
+                        let idx = match canon {
+                            entropyk_core::Z_FLOW => indices.z_flow,
+                            entropyk_core::Z_FLOW_ECO => indices.z_flow_eco,
+                            entropyk_core::Z_DP => indices.z_dp,
+                            entropyk_core::Z_UA => indices.z_ua,
+                            entropyk_core::Z_POWER => indices.z_power,
+                            entropyk_core::Z_ETAV => indices.z_etav,
+                            entropyk_core::F_W => indices.f_w,
+                            _ => None,
+                        };
+                        if let Some(idx) = idx {
+                            cols.insert(idx);
+                        }
+                    }
+                }
+                if let Some(&node_idx) = self.component_names.get(comp) {
+                    for col in self.incident_state_indices_for_component(node_idx) {
+                        cols.insert(col);
+                    }
+                }
+            }
+
+            for col in cols {
+                let orig = state_mut[col];
+                state_mut[col] = orig + eps;
+                let r_plus = eq
+                    .residual(|c, f| self.resolve_expr_ref(c, f, &state_mut))
+                    .unwrap_or(f64::NAN);
+                state_mut[col] = orig - eps;
+                let r_minus = eq
+                    .residual(|c, f| self.resolve_expr_ref(c, f, &state_mut))
+                    .unwrap_or(f64::NAN);
+                state_mut[col] = orig;
+                if r_plus.is_finite() && r_minus.is_finite() {
+                    let dr = (r_plus - r_minus) / (2.0 * eps);
+                    if dr.abs() > 1e-10 {
+                        entries.push((row, col, dr));
+                    }
+                }
+            }
+        }
+        entries
+    }
+
     /// Tolerance for mass balance validation [kg/s].
     ///
     /// This value (1e-9 kg/s) is tight enough to catch numerical issues