sepehr
5425685a48
Some checks failed
CI / check (push) Has been cancelled
Fixed/Free Probe calibration now emits embeddings[] (unknown + equation) instead of controls[], keeping SaturatedController for physical regulation only. Co-authored-by: Cursor <cursoragent@cursor.com>
6377 lines
242 KiB
Rust
6377 lines
242 KiB
Rust
//! System graph structure for thermodynamic simulation.
|
||
//!
|
||
//! This module provides the core graph representation of a thermodynamic system,
|
||
//! where nodes are components and edges represent flow connections. Edges index
|
||
//! into the solver's state vector (ṁ, P and h per edge — see [`EdgeKind`]).
|
||
//!
|
||
//! Multi-circuit support (Story 3.3): A machine can have up to 5 independent
|
||
//! circuits (valid circuit IDs: 0, 1, 2, 3, 4). Each node belongs to exactly one
|
||
//! circuit. Flow edges connect only nodes within the same circuit.
|
||
|
||
use entropyk_components::{
|
||
validate_port_continuity, Component, ComponentError, ConnectionError, JacobianBuilder,
|
||
ResidualVector, StateSlice,
|
||
};
|
||
use petgraph::algo;
|
||
use petgraph::graph::{EdgeIndex, Graph, NodeIndex};
|
||
use petgraph::visit::EdgeRef;
|
||
use petgraph::Directed;
|
||
use std::collections::HashMap;
|
||
|
||
use crate::coupling::{has_circular_dependencies, ThermalCoupling};
|
||
use crate::dof::{
|
||
align_roles, ComponentEquationBlock, DofReport, EquationRole, SystemDofBalance, SystemDofError,
|
||
UnknownKind,
|
||
};
|
||
use crate::error::{AddEdgeError, TopologyError};
|
||
use crate::inverse::{
|
||
BoundedVariable, BoundedVariableError, BoundedVariableId, Constraint, ConstraintError,
|
||
ConstraintId, DoFError, InverseControlConfig, SaturatedController,
|
||
};
|
||
use entropyk_core::{CircuitId, Temperature};
|
||
|
||
/// Maximum circuit ID (inclusive). Machine supports up to 5 circuits.
|
||
pub const MAX_CIRCUIT_ID: u16 = 4;
|
||
|
||
/// Default initial mass-flow seed [kg/s] used by the `SmartInitializer` to
|
||
/// populate the ṁ slot of each edge in the solver state vector. After CM1.3
|
||
/// ṁ is a genuine Newton unknown driven by component mass-flow residuals; this
|
||
/// constant is the initial guess, not a placeholder closure.
|
||
pub const DEFAULT_MASS_FLOW_SEED_KG_S: f64 = 0.05;
|
||
/// Minimum pressure allowed for Newton pressure unknowns [Pa].
|
||
///
|
||
/// This is a solver-domain bound, not a component fallback: trial states below
|
||
/// this pressure are outside the thermodynamic domain used by fluid backends.
|
||
pub const MIN_SOLVER_PRESSURE_PA: f64 = 10_000.0;
|
||
|
||
/// Backward-compatibility alias — tests and tools may still reference the old
|
||
/// TEMP name until they migrate to `DEFAULT_MASS_FLOW_SEED_KG_S` (T6).
|
||
#[deprecated(since = "0.0.0", note = "Use DEFAULT_MASS_FLOW_SEED_KG_S instead")]
|
||
#[allow(dead_code)]
|
||
pub const TEMP_MASS_FLOW_SEED_KG_S: f64 = DEFAULT_MASS_FLOW_SEED_KG_S;
|
||
|
||
/// Aggregate thermodynamic performance of a solved refrigeration/heat-pump cycle.
|
||
///
|
||
/// Built by [`System::cycle_performance`] from each component's
|
||
/// [`entropyk_components::Component::energy_transfers`] evaluated at the solved
|
||
/// state. Duties are classified purely by the physical sign of the reported heat
|
||
/// and work — no component-type strings — so it works for any topology whose
|
||
/// components implement `energy_transfers`.
|
||
#[derive(Debug, Clone, Copy, PartialEq)]
|
||
pub struct CyclePerformance {
|
||
/// Total cooling duty [W]: heat absorbed by the refrigerant (evaporators).
|
||
pub q_cooling_w: f64,
|
||
/// Total heating duty [W]: heat rejected by the refrigerant (condensers/gas coolers).
|
||
pub q_heating_w: f64,
|
||
/// Total shaft/electrical work input [W] (compressors, pumps, fans).
|
||
pub work_input_w: f64,
|
||
/// Number of components that contributed an `energy_transfers` reading.
|
||
pub components_counted: usize,
|
||
}
|
||
|
||
impl CyclePerformance {
|
||
/// Cooling coefficient of performance `EER = Q_cooling / W_input`.
|
||
///
|
||
/// Returns `None` when no work is drawn (undefined ratio).
|
||
pub fn cop_cooling(&self) -> Option<f64> {
|
||
if self.work_input_w > 0.0 {
|
||
Some(self.q_cooling_w / self.work_input_w)
|
||
} else {
|
||
None
|
||
}
|
||
}
|
||
|
||
/// Heating coefficient of performance `COP = Q_heating / W_input`.
|
||
///
|
||
/// Returns `None` when no work is drawn (undefined ratio).
|
||
pub fn cop_heating(&self) -> Option<f64> {
|
||
if self.work_input_w > 0.0 {
|
||
Some(self.q_heating_w / self.work_input_w)
|
||
} else {
|
||
None
|
||
}
|
||
}
|
||
}
|
||
|
||
/// Classifies a flow edge by the physical circuit it belongs to.
|
||
///
|
||
/// The kind determines how many state-vector unknowns the edge carries (see
|
||
/// [`EdgeKind::n_unknowns`]). Today every edge defaults to
|
||
/// [`EdgeKind::Refrigerant`], preserving the historical single-circuit
|
||
/// refrigerant behaviour. Hydraulic (pumped liquid) edges share the same
|
||
/// `(ṁ, P, h)` layout; moist-air edges will gain a fourth unknown (humidity
|
||
/// ratio `W`) in a later story.
|
||
#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
|
||
pub enum EdgeKind {
|
||
/// Refrigerant flow `(ṁ, P, h)` — the default for backward compatibility.
|
||
#[default]
|
||
Refrigerant,
|
||
/// Pumped hydraulic (single-phase liquid) flow `(ṁ, P, h)`.
|
||
Hydraulic,
|
||
/// Moist-air flow. Today `(ṁ, P, h)`; a humidity-ratio unknown `W` is
|
||
/// added in Story CM6.1, at which point `n_unknowns()` returns 4.
|
||
MoistAir,
|
||
}
|
||
|
||
impl EdgeKind {
|
||
/// Number of state-vector unknowns this edge contributes.
|
||
///
|
||
/// Currently 3 (`ṁ, P, h`) for every kind. MoistAir will return 4 once the
|
||
/// humidity-ratio unknown `W` is introduced (Story CM6.1).
|
||
pub fn n_unknowns(&self) -> usize {
|
||
match self {
|
||
EdgeKind::Refrigerant | EdgeKind::Hydraulic | EdgeKind::MoistAir => 3,
|
||
}
|
||
}
|
||
}
|
||
|
||
/// Weight for flow edges in the system graph.
|
||
///
|
||
/// Each edge represents a flow connection between two component ports and stores
|
||
/// the state vector indices for mass flow (ṁ), pressure (P) and enthalpy (h) at
|
||
/// that connection, plus the [`EdgeKind`] classifying the edge.
|
||
///
|
||
/// After [`System::finalize`] runs the topology presolve (CM1.4), all edges in
|
||
/// the same series branch share one `state_index_m` (the branch-level ṁ unknown).
|
||
/// P and h indices remain per-edge.
|
||
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
|
||
pub struct FlowEdge {
|
||
/// Classification of the edge (refrigerant / hydraulic / moist-air).
|
||
pub kind: EdgeKind,
|
||
/// State vector index for mass flow rate (kg/s).
|
||
///
|
||
/// After CM1.4 topology presolve all edges in the same series branch share
|
||
/// this index — the branch-level ṁ unknown.
|
||
pub state_index_m: usize,
|
||
/// State vector index for pressure (Pa)
|
||
pub state_index_p: usize,
|
||
/// State vector index for enthalpy (J/kg)
|
||
pub state_index_h: usize,
|
||
/// Local port index on the **source** component this edge leaves from
|
||
/// (recorded by [`System::add_edge_with_ports`]; defaults to 1 = outlet).
|
||
pub source_port: usize,
|
||
/// Local port index on the **target** component this edge enters at
|
||
/// (recorded by [`System::add_edge_with_ports`]; defaults to 0 = inlet).
|
||
pub target_port: usize,
|
||
/// Branch identifier assigned by [`topology::presolve_mass_flow_topology`].
|
||
///
|
||
/// All edges in the same series branch (no junction in between) receive the
|
||
/// same `mass_flow_branch_id` and share one `state_index_m`. Sentinel value
|
||
/// `usize::MAX` means "not yet assigned" (before `finalize()` is called).
|
||
pub mass_flow_branch_id: usize,
|
||
}
|
||
|
||
impl FlowEdge {
|
||
/// Creates a refrigerant flow edge with unassigned state indices.
|
||
///
|
||
/// State indices and branch id are assigned during [`System::finalize`].
|
||
pub(crate) fn new_unassigned() -> Self {
|
||
FlowEdge {
|
||
kind: EdgeKind::Refrigerant,
|
||
state_index_m: 0,
|
||
state_index_p: 0,
|
||
state_index_h: 0,
|
||
mass_flow_branch_id: usize::MAX,
|
||
source_port: 1,
|
||
target_port: 0,
|
||
}
|
||
}
|
||
|
||
/// Creates an unassigned edge carrying explicit local port indices.
|
||
pub(crate) fn new_unassigned_with_ports(source_port: usize, target_port: usize) -> Self {
|
||
FlowEdge {
|
||
source_port,
|
||
target_port,
|
||
..Self::new_unassigned()
|
||
}
|
||
}
|
||
}
|
||
|
||
/// System graph structure.
|
||
///
|
||
/// Nodes are components (`Box<dyn Component>`), edges are flow connections with
|
||
/// state indices. The state vector layout is:
|
||
///
|
||
/// ```text
|
||
/// [m_edge0, P_edge0, h_edge0, m_edge1, P_edge1, h_edge1, ...]
|
||
/// ```
|
||
///
|
||
/// Edge order follows the graph's internal edge iteration order (stable after
|
||
/// `finalize()` is called).
|
||
pub struct System {
|
||
graph: Graph<Box<dyn Component>, FlowEdge, Directed>,
|
||
/// Maps EdgeIndex to (state_index_m, state_index_p, state_index_h) - built in finalize()
|
||
edge_to_state: HashMap<EdgeIndex, (usize, usize, usize)>,
|
||
/// Maps NodeIndex to CircuitId. Nodes without entry default to circuit 0.
|
||
node_to_circuit: HashMap<NodeIndex, CircuitId>,
|
||
/// Thermal couplings between circuits (heat transfer without fluid mixing).
|
||
thermal_couplings: Vec<ThermalCoupling>,
|
||
/// Constraints for inverse control (output - target = 0)
|
||
constraints: HashMap<ConstraintId, Constraint>,
|
||
/// Bounded control variables for inverse control (with box constraints)
|
||
bounded_variables: HashMap<BoundedVariableId, BoundedVariable>,
|
||
/// Inverse control configuration (constraint → control variable mappings)
|
||
inverse_control: InverseControlConfig,
|
||
/// Saturated PI loops solved as additive `(u, x)` unknown/residual pairs.
|
||
saturated_controllers: Vec<SaturatedController>,
|
||
/// Physical free-actuator unknowns (arch-6). Each entry is a registered
|
||
/// bounded variable id whose value is an *extra* solver unknown appended at
|
||
/// the end of the state vector. Unlike saturated controllers, a free actuator
|
||
/// contributes **no** system residual — the owning component supplies the
|
||
/// closing equation (e.g. an EXV orifice-flow residual) inside its own
|
||
/// `n_equations()` block, keeping the DoF balanced (+1 unknown / +1 eq).
|
||
free_actuators: Vec<BoundedVariableId>,
|
||
/// Registry of component names for constraint validation.
|
||
/// Maps human-readable names (e.g., "evaporator") to NodeIndex.
|
||
component_names: HashMap<String, NodeIndex>,
|
||
/// Per-component calibration indices snapshot, captured at `finalize()` time.
|
||
/// Lets external readers (e.g. result extraction) learn which Z-factors
|
||
/// (z_ua, z_dp, z_flow, …) were promoted to free unknowns and where they
|
||
/// live in the state vector, so the solved values can be surfaced to the user.
|
||
calib_indices_by_name: HashMap<String, entropyk_core::CalibIndices>,
|
||
finalized: bool,
|
||
total_state_len: usize,
|
||
/// When `true` (default), `finalize` rejects **over-constrained** systems.
|
||
/// Under-constrained systems always warn; use [`validate_system_dof`] for a
|
||
/// hard square check in production paths.
|
||
enforce_dof_gate: bool,
|
||
}
|
||
|
||
impl System {
|
||
/// Creates a new empty system graph.
|
||
pub fn new() -> Self {
|
||
Self {
|
||
graph: Graph::new(),
|
||
edge_to_state: HashMap::new(),
|
||
node_to_circuit: HashMap::new(),
|
||
thermal_couplings: Vec::new(),
|
||
constraints: HashMap::new(),
|
||
bounded_variables: HashMap::new(),
|
||
inverse_control: InverseControlConfig::new(),
|
||
saturated_controllers: Vec::new(),
|
||
free_actuators: Vec::new(),
|
||
component_names: HashMap::new(),
|
||
calib_indices_by_name: HashMap::new(),
|
||
finalized: false,
|
||
total_state_len: 0,
|
||
enforce_dof_gate: true,
|
||
}
|
||
}
|
||
|
||
/// Enables or disables the over-constrained DoF hard gate in [`finalize`].
|
||
///
|
||
/// Keep this **enabled** for real-machine models. Disable only for topology
|
||
/// unit tests or intentionally non-physical mock assemblies (e.g. macro
|
||
/// residual plumbing tests). Production solvers should call
|
||
/// [`validate_system_dof`] after finalize for a full square-system check.
|
||
pub fn set_enforce_dof_gate(&mut self, enforce: bool) {
|
||
self.enforce_dof_gate = enforce;
|
||
}
|
||
|
||
/// Adds a component as a node in the default circuit (circuit 0) and returns its node index.
|
||
///
|
||
/// For multi-circuit machines, use [`add_component_to_circuit`](Self::add_component_to_circuit).
|
||
pub fn add_component(&mut self, component: Box<dyn Component>) -> NodeIndex {
|
||
self.add_component_to_circuit(component, CircuitId::ZERO)
|
||
.unwrap()
|
||
}
|
||
|
||
/// Adds a component as a node in the specified circuit and returns its node index.
|
||
///
|
||
/// # Errors
|
||
///
|
||
/// Returns `TopologyError::TooManyCircuits` if `circuit_id.0 > 4`.
|
||
pub fn add_component_to_circuit(
|
||
&mut self,
|
||
component: Box<dyn Component>,
|
||
circuit_id: CircuitId,
|
||
) -> Result<NodeIndex, TopologyError> {
|
||
if circuit_id.0 > MAX_CIRCUIT_ID {
|
||
return Err(TopologyError::TooManyCircuits {
|
||
requested: circuit_id.0,
|
||
});
|
||
}
|
||
self.finalized = false;
|
||
let node_idx = self.graph.add_node(component);
|
||
self.node_to_circuit.insert(node_idx, circuit_id);
|
||
Ok(node_idx)
|
||
}
|
||
|
||
/// Returns the circuit ID for a node, or circuit 0 if not found (backward compat).
|
||
pub fn node_circuit(&self, node: NodeIndex) -> CircuitId {
|
||
self.node_to_circuit
|
||
.get(&node)
|
||
.copied()
|
||
.unwrap_or(CircuitId::ZERO)
|
||
}
|
||
|
||
/// Returns the circuit ID for an edge based on its source node.
|
||
///
|
||
/// # Panics
|
||
///
|
||
/// Panics if the edge index is invalid.
|
||
pub fn edge_circuit(&self, edge: EdgeIndex) -> CircuitId {
|
||
let (src, _tgt) = self.graph.edge_endpoints(edge).expect("invalid edge index");
|
||
self.node_circuit(src)
|
||
}
|
||
|
||
/// Adds a flow edge from `source` to `target` without port validation.
|
||
///
|
||
/// **No port compatibility validation is performed.** Use
|
||
/// [`add_edge_with_ports`](Self::add_edge_with_ports) for components with ports to validate
|
||
/// fluid, pressure, and enthalpy continuity. This method is intended for components without
|
||
/// ports (e.g. mock components in tests).
|
||
///
|
||
/// Flow edges connect only nodes within the same circuit. Cross-circuit connections
|
||
/// are rejected (thermal coupling is Story 3.4).
|
||
///
|
||
/// State indices are assigned when `finalize()` is called.
|
||
///
|
||
/// # Errors
|
||
///
|
||
/// Returns `TopologyError::CrossCircuitConnection` if source and target are in different circuits.
|
||
pub fn add_edge(
|
||
&mut self,
|
||
source: NodeIndex,
|
||
target: NodeIndex,
|
||
) -> Result<EdgeIndex, TopologyError> {
|
||
let src_circuit = self.node_circuit(source);
|
||
let tgt_circuit = self.node_circuit(target);
|
||
if src_circuit != tgt_circuit {
|
||
tracing::warn!(
|
||
"Cross-circuit edge rejected: source circuit {}, target circuit {}",
|
||
src_circuit.0,
|
||
tgt_circuit.0
|
||
);
|
||
return Err(TopologyError::CrossCircuitConnection {
|
||
source_circuit: src_circuit.0,
|
||
target_circuit: tgt_circuit.0,
|
||
});
|
||
}
|
||
|
||
// Safety check: Warn if connecting components with ports using the non-validating method
|
||
if let (Some(src), Some(tgt)) = (
|
||
self.graph.node_weight(source),
|
||
self.graph.node_weight(target),
|
||
) {
|
||
if !src.get_ports().is_empty() || !tgt.get_ports().is_empty() {
|
||
tracing::warn!(
|
||
"add_edge called on components with ports (src: {:?}, tgt: {:?}). \
|
||
This bypasses port validation. Use add_edge_with_ports instead.",
|
||
source,
|
||
target
|
||
);
|
||
}
|
||
}
|
||
|
||
self.finalized = false;
|
||
Ok(self
|
||
.graph
|
||
.add_edge(source, target, FlowEdge::new_unassigned()))
|
||
}
|
||
|
||
/// Adds a flow edge from `source` outlet port to `target` inlet port with validation.
|
||
///
|
||
/// Validates circuit membership (same circuit), then fluid compatibility, pressure and
|
||
/// enthalpy continuity using port.rs tolerances. For 2-port components: `source_port_idx=1`
|
||
/// (outlet), `target_port_idx=0` (inlet).
|
||
///
|
||
/// # Errors
|
||
///
|
||
/// Returns `AddEdgeError::Topology` if source and target are in different circuits.
|
||
/// Returns `AddEdgeError::Connection` if ports are incompatible (fluid, pressure, or enthalpy mismatch).
|
||
pub fn add_edge_with_ports(
|
||
&mut self,
|
||
source: NodeIndex,
|
||
source_port_idx: usize,
|
||
target: NodeIndex,
|
||
target_port_idx: usize,
|
||
) -> Result<EdgeIndex, AddEdgeError> {
|
||
// Circuit validation first
|
||
let src_circuit = self.node_circuit(source);
|
||
let tgt_circuit = self.node_circuit(target);
|
||
if src_circuit != tgt_circuit {
|
||
tracing::warn!(
|
||
"Cross-circuit edge rejected: source circuit {}, target circuit {}",
|
||
src_circuit.0,
|
||
tgt_circuit.0
|
||
);
|
||
return Err(TopologyError::CrossCircuitConnection {
|
||
source_circuit: src_circuit.0,
|
||
target_circuit: tgt_circuit.0,
|
||
}
|
||
.into());
|
||
}
|
||
|
||
let source_comp = self
|
||
.graph
|
||
.node_weight(source)
|
||
.ok_or_else(|| ConnectionError::InvalidNodeIndex(source.index()))?;
|
||
let target_comp = self
|
||
.graph
|
||
.node_weight(target)
|
||
.ok_or_else(|| ConnectionError::InvalidNodeIndex(target.index()))?;
|
||
|
||
let source_ports = source_comp.get_ports();
|
||
let target_ports = target_comp.get_ports();
|
||
|
||
if source_ports.is_empty() && target_ports.is_empty() {
|
||
// No ports: add edge without validation (backward compat), but keep
|
||
// the requested local port indices so multi-port components can
|
||
// resolve their per-port edge wiring in finalize().
|
||
self.finalized = false;
|
||
return Ok(self.graph.add_edge(
|
||
source,
|
||
target,
|
||
FlowEdge::new_unassigned_with_ports(source_port_idx, target_port_idx),
|
||
));
|
||
}
|
||
|
||
if source_port_idx >= source_ports.len() {
|
||
return Err(ConnectionError::InvalidPortIndex {
|
||
index: source_port_idx,
|
||
port_count: source_ports.len(),
|
||
max_index: source_ports.len().saturating_sub(1),
|
||
}
|
||
.into());
|
||
}
|
||
if target_port_idx >= target_ports.len() {
|
||
return Err(ConnectionError::InvalidPortIndex {
|
||
index: target_port_idx,
|
||
port_count: target_ports.len(),
|
||
max_index: target_ports.len().saturating_sub(1),
|
||
}
|
||
.into());
|
||
}
|
||
|
||
let outlet = &source_ports[source_port_idx];
|
||
let inlet = &target_ports[target_port_idx];
|
||
if let Err(ref e) = validate_port_continuity(outlet, inlet) {
|
||
tracing::warn!("Port validation failed: {}", e);
|
||
return Err(e.clone().into());
|
||
}
|
||
|
||
self.finalized = false;
|
||
Ok(self.graph.add_edge(
|
||
source,
|
||
target,
|
||
FlowEdge::new_unassigned_with_ports(source_port_idx, target_port_idx),
|
||
))
|
||
}
|
||
|
||
/// Finalizes the graph: runs topology presolve, builds edge→state index
|
||
/// mapping, and validates topology.
|
||
///
|
||
/// # State vector layout (CM1.4)
|
||
///
|
||
/// After the mass-flow topology presolve (`|B|` branches detected):
|
||
///
|
||
/// ```text
|
||
/// [ṁ_B0, ṁ_B1, …, ṁ_{B-1}, ← one slot per branch
|
||
/// P_0, h_0, P_1, h_1, …] ← two slots per edge (P and h)
|
||
/// ```
|
||
///
|
||
/// Total edge-section length = `|B| + 2 × |E|`.
|
||
///
|
||
/// All edges in the same series branch share the same `state_index_m`.
|
||
/// `state_index_p` and `state_index_h` remain unique per edge.
|
||
///
|
||
/// # Errors
|
||
///
|
||
/// Returns `TopologyError` if:
|
||
/// - Any node is isolated (no edges)
|
||
/// - The graph is empty (no components)
|
||
pub fn finalize(&mut self) -> Result<(), TopologyError> {
|
||
self.validate_topology()?;
|
||
|
||
if !self.thermal_couplings.is_empty() && has_circular_dependencies(self.thermal_couplings())
|
||
{
|
||
tracing::warn!("Circular thermal coupling detected, simultaneous solving required");
|
||
}
|
||
|
||
// CM1.4: run topology presolve to assign mass_flow_branch_id to each edge.
|
||
let n_branches = crate::topology::presolve_mass_flow_topology(&mut self.graph);
|
||
tracing::debug!(
|
||
n_branches,
|
||
edge_count = self.graph.edge_count(),
|
||
"Mass-flow topology presolve complete"
|
||
);
|
||
|
||
// ── Pass 1: allocate one ṁ state slot per branch ─────────────────
|
||
// Branch ids are contiguous 0..n_branches, so we can use a Vec.
|
||
let mut branch_m_idx: Vec<Option<usize>> = vec![None; n_branches];
|
||
let mut idx = 0usize; // running state-vector position
|
||
|
||
// Allocate ṁ slots for branches in branch_id order.
|
||
// We need branch_id → m_idx. Walk edges to collect branch ids and
|
||
// allocate one slot per branch in the order first encountered.
|
||
for edge_idx in self.graph.edge_indices() {
|
||
let branch_id = self
|
||
.graph
|
||
.edge_weight(edge_idx)
|
||
.map(|w| w.mass_flow_branch_id)
|
||
.unwrap_or(0);
|
||
if branch_id < n_branches && branch_m_idx[branch_id].is_none() {
|
||
branch_m_idx[branch_id] = Some(idx);
|
||
idx += 1;
|
||
}
|
||
}
|
||
|
||
// ── Pass 2: assign per-edge (m_idx, p_idx, h_idx) ────────────────
|
||
self.edge_to_state.clear();
|
||
for edge_idx in self.graph.edge_indices() {
|
||
let (branch_id, kind) = self
|
||
.graph
|
||
.edge_weight(edge_idx)
|
||
.map(|w| (w.mass_flow_branch_id, w.kind))
|
||
.unwrap_or((0, EdgeKind::default()));
|
||
|
||
// Shared ṁ index from the branch's pre-allocated slot.
|
||
let m_idx = branch_m_idx.get(branch_id).copied().flatten().unwrap_or(0);
|
||
|
||
// Per-edge P and h indices (always unique).
|
||
let p_idx = idx;
|
||
let h_idx = idx + 1;
|
||
idx += kind.n_unknowns() - 1; // subtract 1 because ṁ is shared
|
||
|
||
self.edge_to_state.insert(edge_idx, (m_idx, p_idx, h_idx));
|
||
if let Some(weight) = self.graph.edge_weight_mut(edge_idx) {
|
||
weight.state_index_m = m_idx;
|
||
weight.state_index_p = p_idx;
|
||
weight.state_index_h = h_idx;
|
||
}
|
||
}
|
||
self.finalized = true;
|
||
|
||
// Notify each component about its position in the global state vector.
|
||
// Collect context first (to avoid borrow conflicts), then apply.
|
||
//
|
||
// State offset: for the parent System the state layout is a flat array of
|
||
// sum(EdgeKind::n_unknowns()) entries (3 per edge today) for the
|
||
// *parent's own* edges. An embedded MacroComponent has its internal
|
||
// state appended after the parent edges, addressed via its own
|
||
// global_state_offset. We start with the total per-edge unknown count
|
||
// as the base offset, and accumulate `internal_state_len` for each
|
||
// component.
|
||
let mut current_offset = idx;
|
||
|
||
// Gather (node_idx, offset, incident_edge_indices) before mutating nodes.
|
||
// Each entry carries a (m_idx, p_idx, h_idx) triple per incident edge so
|
||
// components can wire their mass-flow residuals to the edge ṁ unknowns (CM1.3).
|
||
#[allow(clippy::type_complexity)]
|
||
let mut node_context: Vec<(
|
||
petgraph::graph::NodeIndex,
|
||
usize,
|
||
Vec<(usize, usize, usize)>,
|
||
Vec<Option<(usize, usize, usize)>>,
|
||
)> = Vec::new();
|
||
for node_idx in self.graph.node_indices() {
|
||
let component = self.graph.node_weight(node_idx).unwrap();
|
||
let mut incident: Vec<(usize, usize, usize)> = Vec::new();
|
||
// Deterministic per-port wiring: port_edges[i] = edge state indices
|
||
// of the edge attached at local port i (incoming edges use the
|
||
// edge's target_port, outgoing edges its source_port).
|
||
let mut port_edges: Vec<Option<(usize, usize, usize)>> = Vec::new();
|
||
let mut record_port = |slot: usize, triple: (usize, usize, usize)| {
|
||
if port_edges.len() <= slot {
|
||
port_edges.resize(slot + 1, None);
|
||
}
|
||
port_edges[slot] = Some(triple);
|
||
};
|
||
for edge_ref in self
|
||
.graph
|
||
.edges_directed(node_idx, petgraph::Direction::Incoming)
|
||
{
|
||
if let Some(&(m, p, h)) = self.edge_to_state.get(&edge_ref.id()) {
|
||
incident.push((m, p, h));
|
||
record_port(edge_ref.weight().target_port, (m, p, h));
|
||
}
|
||
}
|
||
for edge_ref in self
|
||
.graph
|
||
.edges_directed(node_idx, petgraph::Direction::Outgoing)
|
||
{
|
||
if let Some(&(m, p, h)) = self.edge_to_state.get(&edge_ref.id()) {
|
||
incident.push((m, p, h));
|
||
record_port(edge_ref.weight().source_port, (m, p, h));
|
||
}
|
||
}
|
||
node_context.push((node_idx, current_offset, incident, port_edges));
|
||
|
||
// Advance the global offset by this component's internal state length
|
||
current_offset += component.internal_state_len();
|
||
}
|
||
|
||
self.total_state_len = current_offset;
|
||
|
||
// Notify components about their calibration control variables (Story 5.5)
|
||
let mut comp_calib_indices: HashMap<String, entropyk_core::CalibIndices> = HashMap::new();
|
||
for (index, id) in self.inverse_control.linked_controls().enumerate() {
|
||
if let Some(bounded_var) = self.bounded_variables.get(id) {
|
||
if let Some(comp_id) = bounded_var.component_id() {
|
||
let indices = comp_calib_indices.entry(comp_id.to_string()).or_default();
|
||
let state_idx = self.total_state_len + index;
|
||
|
||
let id_str = id.as_str();
|
||
if let Some(factor) = entropyk_core::id_ends_with_calib_suffix(id_str) {
|
||
match factor {
|
||
entropyk_core::Z_FLOW => indices.z_flow = Some(state_idx),
|
||
entropyk_core::Z_FLOW_ECO => indices.z_flow_eco = Some(state_idx),
|
||
entropyk_core::Z_DP => indices.z_dp = Some(state_idx),
|
||
entropyk_core::Z_UA => indices.z_ua = Some(state_idx),
|
||
entropyk_core::Z_POWER => indices.z_power = Some(state_idx),
|
||
entropyk_core::Z_ETAV => indices.z_etav = Some(state_idx),
|
||
entropyk_core::F_W => indices.f_w = Some(state_idx),
|
||
_ => {}
|
||
}
|
||
}
|
||
}
|
||
}
|
||
}
|
||
let saturated_base = self.total_state_len
|
||
+ self.inverse_control.mapping_count()
|
||
+ self.coupling_residual_count();
|
||
for (index, ctrl) in self.saturated_controllers.iter().enumerate() {
|
||
if let Some(bounded_var) = self.bounded_variables.get(ctrl.actuator()) {
|
||
if let Some(comp_id) = bounded_var.component_id() {
|
||
let indices = comp_calib_indices.entry(comp_id.to_string()).or_default();
|
||
let state_idx = saturated_base + 2 * index;
|
||
|
||
let id_str = ctrl.actuator().as_str();
|
||
if let Some(factor) = entropyk_core::id_ends_with_calib_suffix(id_str) {
|
||
match factor {
|
||
entropyk_core::Z_FLOW => indices.z_flow = Some(state_idx),
|
||
entropyk_core::Z_FLOW_ECO => indices.z_flow_eco = Some(state_idx),
|
||
entropyk_core::Z_DP => indices.z_dp = Some(state_idx),
|
||
entropyk_core::Z_UA => indices.z_ua = Some(state_idx),
|
||
entropyk_core::Z_POWER => indices.z_power = Some(state_idx),
|
||
entropyk_core::Z_ETAV => indices.z_etav = Some(state_idx),
|
||
entropyk_core::F_W => indices.f_w = Some(state_idx),
|
||
_ => {}
|
||
}
|
||
} else if id_str.ends_with("injection") || id_str.ends_with("actuator") {
|
||
// Generic physical actuator (arch-6), e.g. liquid injection
|
||
// ratio φ_inj: the controls[] loop drives this slot directly.
|
||
indices.actuator = Some(state_idx);
|
||
}
|
||
}
|
||
}
|
||
}
|
||
|
||
// Physical free actuators (arch-6): wire each free-actuator bounded var to
|
||
// its component's generic `actuator` CalibIndices slot. These unknowns live
|
||
// AFTER the saturated block. No system residual is emitted for them — the
|
||
// owning component supplies the closing equation in its own residual block.
|
||
let free_actuator_base = self.total_state_len
|
||
+ self.inverse_control.mapping_count()
|
||
+ self.coupling_residual_count()
|
||
+ 2 * self.saturated_controllers.len();
|
||
for (index, actuator_id) in self.free_actuators.iter().enumerate() {
|
||
if let Some(bounded_var) = self.bounded_variables.get(actuator_id) {
|
||
if let Some(comp_id) = bounded_var.component_id() {
|
||
let indices = comp_calib_indices.entry(comp_id.to_string()).or_default();
|
||
indices.actuator = Some(free_actuator_base + index);
|
||
}
|
||
}
|
||
}
|
||
|
||
// Now mutate each node weight (component) with the gathered context.
|
||
for (node_idx, offset, incident, port_edges) in node_context {
|
||
if let Some(component) = self.graph.node_weight_mut(node_idx) {
|
||
component.set_system_context(offset, &incident);
|
||
component.set_port_context(&port_edges);
|
||
|
||
// If we registered a name for this node, check if we have calib indices for it
|
||
if let Some((name, _)) = self.component_names.iter().find(|(_, &n)| n == node_idx) {
|
||
if let Some(&indices) = comp_calib_indices.get(name) {
|
||
component.set_calib_indices(indices);
|
||
}
|
||
}
|
||
}
|
||
}
|
||
|
||
// Persist the per-component calib map so external readers (result
|
||
// extraction) can surface solved Z-factor values (z_ua, z_dp, …).
|
||
self.calib_indices_by_name = comp_calib_indices;
|
||
|
||
// Wire physical thermal couplings (Story 3.4 completion): each coupling
|
||
// owns one state unknown Q [W] at `coupling_state_index(i)`. The
|
||
// cold-side receiver component (e.g. `ThermalLoad`) reads Q in its
|
||
// energy balance via `set_external_heat_index`; the coupling residual
|
||
// closes Q against the hot component's measured duty (see
|
||
// `compute_residuals`).
|
||
let coupling_wirings: Vec<(String, usize)> = self
|
||
.thermal_couplings
|
||
.iter()
|
||
.enumerate()
|
||
.filter_map(|(i, c)| {
|
||
c.cold_component
|
||
.clone()
|
||
.map(|name| (name, self.coupling_state_index(i)))
|
||
})
|
||
.collect();
|
||
for (name, q_idx) in coupling_wirings {
|
||
if let Some(&node) = self.component_names.get(&name) {
|
||
if let Some(component) = self.graph.node_weight_mut(node) {
|
||
component.set_external_heat_index(q_idx);
|
||
}
|
||
} else {
|
||
tracing::warn!(
|
||
cold_component = %name,
|
||
"Thermal coupling references an unregistered cold component; \
|
||
the coupled heat will not be injected"
|
||
);
|
||
}
|
||
}
|
||
for coupling in &self.thermal_couplings {
|
||
if let Some(hot) = &coupling.hot_component {
|
||
if !self.component_names.contains_key(hot) {
|
||
tracing::warn!(
|
||
hot_component = %hot,
|
||
"Thermal coupling references an unregistered hot component; \
|
||
the transferred duty will be zero"
|
||
);
|
||
}
|
||
}
|
||
}
|
||
|
||
if !self.constraints.is_empty() {
|
||
match self.validate_inverse_control_dof() {
|
||
Ok(()) => {
|
||
tracing::debug!(
|
||
constraint_count = self.constraints.len(),
|
||
control_count = self.inverse_control.mapping_count(),
|
||
"Inverse control DoF validation passed"
|
||
);
|
||
}
|
||
Err(DoFError::UnderConstrainedSystem { .. }) => {
|
||
tracing::warn!(
|
||
constraint_count = self.constraints.len(),
|
||
control_count = self.inverse_control.mapping_count(),
|
||
"Under-constrained inverse control system - solver may still converge"
|
||
);
|
||
}
|
||
Err(DoFError::OverConstrainedSystem {
|
||
constraint_count,
|
||
control_count,
|
||
equation_count,
|
||
unknown_count,
|
||
}) => {
|
||
tracing::warn!(
|
||
constraint_count,
|
||
control_count,
|
||
equation_count,
|
||
unknown_count,
|
||
"Over-constrained inverse control system - add more control variables or remove constraints"
|
||
);
|
||
}
|
||
Err(e) => {
|
||
tracing::warn!("Inverse control DoF validation error: {}", e);
|
||
}
|
||
}
|
||
}
|
||
|
||
// DoF gate:
|
||
// - Over-constrained systems are always rejected (no free lunch: every
|
||
// extra residual needs a free unknown).
|
||
// - Under-constrained systems log a warning here so pure topology unit
|
||
// tests with mock components can still exercise graph plumbing.
|
||
// Production paths (CLI run, SystemBuilder) must call
|
||
// `validate_system_dof()` which requires a fully square system.
|
||
let report = self.dof_report();
|
||
match report.balance {
|
||
SystemDofBalance::Balanced => {
|
||
tracing::debug!("{}", report.summary());
|
||
}
|
||
SystemDofBalance::UnderConstrained { free_dofs } => {
|
||
tracing::warn!(free_dofs, "{}", report.summary());
|
||
}
|
||
SystemDofBalance::OverConstrained { excess_equations } => {
|
||
tracing::error!(excess_equations, "{}", report.summary());
|
||
if self.enforce_dof_gate {
|
||
return Err(TopologyError::DofImbalance {
|
||
message: report.summary(),
|
||
});
|
||
}
|
||
tracing::warn!(
|
||
"DoF gate disabled (set_enforce_dof_gate(false)); over-constrained system accepted for tests only"
|
||
);
|
||
}
|
||
}
|
||
|
||
Ok(())
|
||
}
|
||
|
||
/// Validates the topology: no isolated nodes and edge circuit consistency.
|
||
///
|
||
/// Note: "All ports connected" validation requires port→edge association
|
||
/// (Story 3.2 Port Compatibility Validation).
|
||
fn validate_topology(&self) -> Result<(), TopologyError> {
|
||
let node_count = self.graph.node_count();
|
||
if node_count == 0 {
|
||
return Ok(());
|
||
}
|
||
|
||
for node_idx in self.graph.node_indices() {
|
||
let degree = self
|
||
.graph
|
||
.edges_directed(node_idx, petgraph::Direction::Incoming)
|
||
.count()
|
||
+ self
|
||
.graph
|
||
.edges_directed(node_idx, petgraph::Direction::Outgoing)
|
||
.count();
|
||
if degree == 0 {
|
||
return Err(TopologyError::IsolatedNode {
|
||
node_index: node_idx.index(),
|
||
});
|
||
}
|
||
}
|
||
|
||
// Validate that all edges connect nodes within the same circuit
|
||
for edge_idx in self.graph.edge_indices() {
|
||
if let Some((src, tgt)) = self.graph.edge_endpoints(edge_idx) {
|
||
let src_circuit = self.node_circuit(src);
|
||
let tgt_circuit = self.node_circuit(tgt);
|
||
if src_circuit != tgt_circuit {
|
||
return Err(TopologyError::CrossCircuitConnection {
|
||
source_circuit: src_circuit.0,
|
||
target_circuit: tgt_circuit.0,
|
||
});
|
||
}
|
||
}
|
||
}
|
||
|
||
Ok(())
|
||
}
|
||
|
||
/// Returns the documented state vector layout.
|
||
///
|
||
/// Layout: `[P_edge0, h_edge0, P_edge1, h_edge1, ...]` where each edge (in
|
||
/// graph iteration order) contributes 2 entries: pressure (Pa) then enthalpy (J/kg).
|
||
///
|
||
/// # Panics
|
||
///
|
||
/// Panics if `finalize()` has not been called.
|
||
pub fn state_layout(&self) -> &'static str {
|
||
assert!(self.finalized, "call finalize() before state_layout()");
|
||
"[m_edge0, P_edge0, h_edge0, m_edge1, P_edge1, h_edge1, ...] — 3 per edge (mass flow kg/s, pressure Pa, enthalpy J/kg)"
|
||
}
|
||
|
||
/// Returns the length of the physical state vector.
|
||
///
|
||
/// After CM1.4 topology presolve the layout is:
|
||
/// `|B| (branch ṁ slots) + 2 × |E| (per-edge P and h) + internal_component_state`.
|
||
///
|
||
/// Note: This returns the physical state vector length. For the full solver state vector
|
||
/// including control variables, use [`full_state_vector_len`](Self::full_state_vector_len).
|
||
///
|
||
/// # Panics
|
||
///
|
||
/// Panics if `finalize()` has not been called.
|
||
pub fn state_vector_len(&self) -> usize {
|
||
assert!(self.finalized, "call finalize() before state_vector_len()");
|
||
self.total_state_len
|
||
}
|
||
|
||
/// Returns the state indices (P, h) for the given edge.
|
||
///
|
||
/// # Returns
|
||
///
|
||
/// `(state_index_p, state_index_h)` for the edge. The mass-flow index is
|
||
/// available via [`edge_mass_flow_index`](Self::edge_mass_flow_index) or
|
||
/// [`edge_state_indices_full`](Self::edge_state_indices_full).
|
||
///
|
||
/// # Panics
|
||
///
|
||
/// Panics if `finalize()` has not been called or if `edge_id` is invalid.
|
||
pub fn edge_state_indices(&self, edge_id: EdgeIndex) -> (usize, usize) {
|
||
assert!(
|
||
self.finalized,
|
||
"call finalize() before edge_state_indices()"
|
||
);
|
||
let (_m, p, h) = *self
|
||
.edge_to_state
|
||
.get(&edge_id)
|
||
.expect("invalid edge index");
|
||
(p, h)
|
||
}
|
||
|
||
/// Returns the full state indices `(ṁ, P, h)` for the given edge.
|
||
///
|
||
/// # Panics
|
||
///
|
||
/// Panics if `finalize()` has not been called or if `edge_id` is invalid.
|
||
pub fn edge_state_indices_full(&self, edge_id: EdgeIndex) -> (usize, usize, usize) {
|
||
assert!(
|
||
self.finalized,
|
||
"call finalize() before edge_state_indices_full()"
|
||
);
|
||
*self
|
||
.edge_to_state
|
||
.get(&edge_id)
|
||
.expect("invalid edge index")
|
||
}
|
||
|
||
/// Returns the mass-flow (ṁ) state index for the given edge.
|
||
///
|
||
/// # Panics
|
||
///
|
||
/// Panics if `finalize()` has not been called or if `edge_id` is invalid.
|
||
pub fn edge_mass_flow_index(&self, edge_id: EdgeIndex) -> usize {
|
||
self.edge_state_indices_full(edge_id).0
|
||
}
|
||
|
||
/// Returns the [`EdgeKind`] of the given edge.
|
||
///
|
||
/// # Panics
|
||
///
|
||
/// Panics if `edge_id` is invalid.
|
||
pub fn edge_kind(&self, edge_id: EdgeIndex) -> EdgeKind {
|
||
self.graph
|
||
.edge_weight(edge_id)
|
||
.map(|w| w.kind)
|
||
.expect("invalid edge index")
|
||
}
|
||
|
||
/// Number of per-edge mass-flow closure equations.
|
||
///
|
||
/// Returns 0 since CM1.3: ṁ residuals are now contributed by individual
|
||
/// components (`compute_residuals`), so no extra closure rows are needed.
|
||
/// Kept for API compatibility; callers should eventually remove their
|
||
/// `+ system.mass_flow_closure_count()` calls.
|
||
pub fn mass_flow_closure_count(&self) -> usize {
|
||
0
|
||
}
|
||
|
||
/// Returns the number of edges in the graph.
|
||
pub fn edge_count(&self) -> usize {
|
||
self.graph.edge_count()
|
||
}
|
||
|
||
/// Returns an iterator over all edge indices in the graph.
|
||
pub fn edge_indices(&self) -> impl Iterator<Item = EdgeIndex> + '_ {
|
||
self.graph.edge_indices()
|
||
}
|
||
|
||
/// Returns the source and target node indices for the given edge.
|
||
///
|
||
/// Returns `None` if the edge index is invalid.
|
||
pub fn edge_endpoints(&self, edge: EdgeIndex) -> Option<(NodeIndex, NodeIndex)> {
|
||
self.graph.edge_endpoints(edge)
|
||
}
|
||
|
||
/// Returns a reference to the internal graph.
|
||
pub fn graph(&self) -> &Graph<Box<dyn Component>, FlowEdge, Directed> {
|
||
&self.graph
|
||
}
|
||
|
||
/// Returns a reference to the node-to-circuit mapping.
|
||
pub fn node_to_circuit(&self) -> &HashMap<NodeIndex, CircuitId> {
|
||
&self.node_to_circuit
|
||
}
|
||
|
||
/// Returns a reference to the constraints map.
|
||
pub fn constraints_map(&self) -> &HashMap<ConstraintId, Constraint> {
|
||
&self.constraints
|
||
}
|
||
|
||
/// Returns a reference to the bounded variables map.
|
||
pub fn bounded_variables_map(&self) -> &HashMap<BoundedVariableId, BoundedVariable> {
|
||
&self.bounded_variables
|
||
}
|
||
|
||
/// Returns the number of nodes (components) in the graph.
|
||
pub fn node_count(&self) -> usize {
|
||
self.graph.node_count()
|
||
}
|
||
|
||
/// Returns the number of distinct circuits in the machine.
|
||
///
|
||
/// Circuits are identified by the unique circuit IDs present in `node_to_circuit`.
|
||
/// Empty system returns 0. Systems with components always return >= 1 since
|
||
/// all components are assigned to a circuit (defaulting to circuit 0).
|
||
/// Valid circuit IDs are 0 through 4 (inclusive), supporting up to 5 circuits.
|
||
pub fn circuit_count(&self) -> usize {
|
||
if self.graph.node_count() == 0 {
|
||
return 0;
|
||
}
|
||
let mut ids: Vec<u16> = self.node_to_circuit.values().map(|c| c.0).collect();
|
||
if ids.is_empty() {
|
||
// This shouldn't happen since add_component adds to node_to_circuit,
|
||
// but handle defensively
|
||
return 1;
|
||
}
|
||
ids.sort_unstable();
|
||
ids.dedup();
|
||
ids.len()
|
||
}
|
||
|
||
/// Returns an iterator over node indices belonging to the given circuit.
|
||
pub fn circuit_nodes(&self, circuit_id: CircuitId) -> impl Iterator<Item = NodeIndex> + '_ {
|
||
self.graph.node_indices().filter(move |&idx| {
|
||
self.node_to_circuit
|
||
.get(&idx)
|
||
.copied()
|
||
.unwrap_or(CircuitId::ZERO)
|
||
== circuit_id
|
||
})
|
||
}
|
||
|
||
/// Returns an iterator over edge indices belonging to the given circuit.
|
||
///
|
||
/// An edge belongs to a circuit if both its source and target nodes are in that circuit.
|
||
pub fn circuit_edges(&self, circuit_id: CircuitId) -> impl Iterator<Item = EdgeIndex> + '_ {
|
||
self.graph.edge_indices().filter(move |&edge_idx| {
|
||
let (src, tgt) = self.graph.edge_endpoints(edge_idx).expect("valid edge");
|
||
self.node_circuit(src) == circuit_id && self.node_circuit(tgt) == circuit_id
|
||
})
|
||
}
|
||
|
||
/// Checks if a circuit has any components.
|
||
fn circuit_exists(&self, circuit_id: CircuitId) -> bool {
|
||
self.node_to_circuit.values().any(|&c| c == circuit_id)
|
||
}
|
||
|
||
/// Adds a thermal coupling between two circuits.
|
||
///
|
||
/// Thermal couplings represent heat exchangers that transfer heat from a "hot"
|
||
/// circuit to a "cold" circuit without fluid mixing. Heat flows from hot to cold
|
||
/// proportional to the temperature difference and thermal conductance (UA).
|
||
///
|
||
/// # Arguments
|
||
///
|
||
/// * `coupling` - The thermal coupling to add
|
||
///
|
||
/// # Returns
|
||
///
|
||
/// The index of the added coupling in the internal storage.
|
||
///
|
||
/// # Errors
|
||
///
|
||
/// Returns `TopologyError::InvalidCircuitForCoupling` if either circuit
|
||
/// referenced in the coupling does not exist in the system.
|
||
///
|
||
/// # Example
|
||
///
|
||
/// ```no_run
|
||
/// use entropyk_solver::{System, ThermalCoupling, CircuitId};
|
||
/// use entropyk_core::ThermalConductance;
|
||
/// use entropyk_components::Component;
|
||
/// # fn make_mock() -> Box<dyn Component> { unimplemented!() }
|
||
///
|
||
/// let mut sys = System::new();
|
||
/// sys.add_component_to_circuit(make_mock(), CircuitId(0)).unwrap();
|
||
/// sys.add_component_to_circuit(make_mock(), CircuitId(1)).unwrap();
|
||
///
|
||
/// let coupling = ThermalCoupling::new(
|
||
/// CircuitId(0),
|
||
/// CircuitId(1),
|
||
/// ThermalConductance::from_watts_per_kelvin(1000.0),
|
||
/// );
|
||
/// let idx = sys.add_thermal_coupling(coupling).unwrap();
|
||
/// ```
|
||
pub fn add_thermal_coupling(
|
||
&mut self,
|
||
coupling: ThermalCoupling,
|
||
) -> Result<usize, TopologyError> {
|
||
// Validate that both circuits exist
|
||
if !self.circuit_exists(coupling.hot_circuit) {
|
||
return Err(TopologyError::InvalidCircuitForCoupling {
|
||
circuit_id: coupling.hot_circuit.0,
|
||
});
|
||
}
|
||
if !self.circuit_exists(coupling.cold_circuit) {
|
||
return Err(TopologyError::InvalidCircuitForCoupling {
|
||
circuit_id: coupling.cold_circuit.0,
|
||
});
|
||
}
|
||
|
||
self.finalized = false;
|
||
self.thermal_couplings.push(coupling);
|
||
Ok(self.thermal_couplings.len() - 1)
|
||
}
|
||
|
||
/// Returns the number of thermal couplings in the system.
|
||
pub fn thermal_coupling_count(&self) -> usize {
|
||
self.thermal_couplings.len()
|
||
}
|
||
|
||
/// Returns a reference to all thermal couplings.
|
||
pub fn thermal_couplings(&self) -> &[ThermalCoupling] {
|
||
&self.thermal_couplings
|
||
}
|
||
|
||
/// Returns a reference to a specific thermal coupling by index.
|
||
pub fn get_thermal_coupling(&self, index: usize) -> Option<&ThermalCoupling> {
|
||
self.thermal_couplings.get(index)
|
||
}
|
||
|
||
// ────────────────────────────────────────────────────────────────────────
|
||
// Component Name Registry (for constraint validation)
|
||
// ────────────────────────────────────────────────────────────────────────
|
||
|
||
/// Registers a human-readable name for a component node.
|
||
///
|
||
/// This name can be used in constraints to reference the component.
|
||
/// For example, register "evaporator" to reference it in a superheat constraint.
|
||
///
|
||
/// # Arguments
|
||
///
|
||
/// * `name` - Human-readable name for the component (e.g., "evaporator", "condenser")
|
||
/// * `node` - The NodeIndex returned from `add_component()` or `add_component_to_circuit()`
|
||
///
|
||
/// # Returns
|
||
///
|
||
/// `true` if the name was newly registered, `false` if the name was already in use
|
||
/// (in which case the mapping is updated to the new node).
|
||
///
|
||
/// # Example
|
||
///
|
||
/// ```rust,ignore
|
||
/// let mut sys = System::new();
|
||
/// let evap_node = sys.add_component(make_evaporator());
|
||
/// sys.register_component_name("evaporator", evap_node);
|
||
///
|
||
/// // Now constraints can reference "evaporator"
|
||
/// let constraint = Constraint::new(
|
||
/// ConstraintId::new("superheat_control"),
|
||
/// ComponentOutput::Superheat { component_id: "evaporator".to_string() },
|
||
/// 5.0,
|
||
/// );
|
||
/// sys.add_constraint(constraint)?; // Validates "evaporator" exists
|
||
/// ```
|
||
pub fn register_component_name(&mut self, name: &str, node: NodeIndex) -> bool {
|
||
self.component_names
|
||
.insert(name.to_string(), node)
|
||
.is_none()
|
||
}
|
||
|
||
/// Returns the NodeIndex for a registered component name, or None if not found.
|
||
pub fn get_component_node(&self, name: &str) -> Option<NodeIndex> {
|
||
self.component_names.get(name).copied()
|
||
}
|
||
|
||
/// Returns the names of all registered components.
|
||
pub fn registered_component_names(&self) -> impl Iterator<Item = &str> {
|
||
self.component_names.keys().map(|s| s.as_str())
|
||
}
|
||
|
||
/// Per-component `CalibIndices` snapshot captured at `finalize()`.
|
||
/// Each `Some(idx)` slot (z_ua, z_dp, z_flow, z_flow_eco, z_power, z_etav,
|
||
/// actuator) means that factor was promoted to a free solver unknown and
|
||
/// its solved value lives at `state[idx]`. Used by result extraction to
|
||
/// surface solved calibration/control values to the user.
|
||
pub fn calib_indices_by_name(&self) -> &HashMap<String, entropyk_core::CalibIndices> {
|
||
&self.calib_indices_by_name
|
||
}
|
||
|
||
/// Returns a reference to the component stored at the given node index.
|
||
///
|
||
/// # Panics
|
||
///
|
||
/// Panics if the node index is invalid.
|
||
pub fn component(&self, node: NodeIndex) -> &dyn Component {
|
||
self.graph
|
||
.node_weight(node)
|
||
.expect("invalid node index")
|
||
.as_ref()
|
||
}
|
||
|
||
/// Returns a mutable reference to the component at the given node index.
|
||
///
|
||
/// Returns `None` if the node index is invalid.
|
||
/// Used for post-build injection of fluid backends via the builder.
|
||
pub fn component_mut(&mut self, node: NodeIndex) -> Option<&mut dyn Component> {
|
||
let weight = self.graph.node_weight_mut(node)?;
|
||
Some(weight.as_mut())
|
||
}
|
||
|
||
// ────────────────────────────────────────────────────────────────────────
|
||
// Constraint Management (Inverse Control)
|
||
// ────────────────────────────────────────────────────────────────────────
|
||
|
||
/// Adds a constraint for inverse control.
|
||
///
|
||
/// Constraints define desired output conditions. During solving, the solver
|
||
/// will attempt to find input values that satisfy all constraints.
|
||
///
|
||
/// # Arguments
|
||
///
|
||
/// * `constraint` - The constraint to add
|
||
///
|
||
/// # Errors
|
||
///
|
||
/// Returns `ConstraintError::DuplicateId` if a constraint with the same ID
|
||
/// already exists.
|
||
/// Returns `ConstraintError::InvalidReference` if the component referenced
|
||
/// in the constraint's output has not been registered via `register_component_name()`.
|
||
///
|
||
/// # Example
|
||
///
|
||
/// ```rust,ignore
|
||
/// use entropyk_solver::inverse::{Constraint, ConstraintId, ComponentOutput};
|
||
///
|
||
/// let constraint = Constraint::new(
|
||
/// ConstraintId::new("superheat_control"),
|
||
/// ComponentOutput::Superheat {
|
||
/// component_id: "evaporator".to_string()
|
||
/// },
|
||
/// 5.0, // target: 5K superheat
|
||
/// );
|
||
///
|
||
/// system.add_constraint(constraint)?;
|
||
/// ```
|
||
pub fn add_constraint(&mut self, constraint: Constraint) -> Result<(), ConstraintError> {
|
||
let id = constraint.id().clone();
|
||
if self.constraints.contains_key(&id) {
|
||
return Err(ConstraintError::DuplicateId { id });
|
||
}
|
||
|
||
// AC2: Validate that the component referenced in the constraint exists
|
||
let component_id = constraint.output().component_id();
|
||
if !self.component_names.contains_key(component_id) {
|
||
return Err(ConstraintError::InvalidReference {
|
||
component_id: component_id.to_string(),
|
||
});
|
||
}
|
||
|
||
self.constraints.insert(id, constraint);
|
||
Ok(())
|
||
}
|
||
|
||
/// Removes a constraint by ID.
|
||
///
|
||
/// # Arguments
|
||
///
|
||
/// * `id` - The constraint identifier to remove
|
||
///
|
||
/// # Returns
|
||
///
|
||
/// The removed constraint, or `None` if no constraint with that ID exists.
|
||
pub fn remove_constraint(&mut self, id: &ConstraintId) -> Option<Constraint> {
|
||
self.inverse_control.unlink_constraint(id);
|
||
self.constraints.remove(id)
|
||
}
|
||
|
||
/// Returns the number of constraints in the system.
|
||
pub fn constraint_count(&self) -> usize {
|
||
self.constraints.len()
|
||
}
|
||
|
||
/// Returns a reference to all constraints.
|
||
pub fn constraints(&self) -> impl Iterator<Item = &Constraint> {
|
||
self.constraints.values()
|
||
}
|
||
|
||
/// Returns a reference to a specific constraint by ID.
|
||
pub fn get_constraint(&self, id: &ConstraintId) -> Option<&Constraint> {
|
||
self.constraints.get(id)
|
||
}
|
||
|
||
/// Returns the number of constraint residual equations.
|
||
///
|
||
/// Each constraint adds one equation to the residual vector.
|
||
pub fn constraint_residual_count(&self) -> usize {
|
||
self.constraints.len()
|
||
}
|
||
|
||
/// Computes constraint residuals and appends them to the provided vector.
|
||
///
|
||
/// This method computes the residual for each constraint:
|
||
///
|
||
/// $$r_{constraint} = f(x) - y_{target}$$
|
||
///
|
||
/// where $f(x)$ is the measured output value and $y_{target}$ is the constraint target.
|
||
///
|
||
/// # Arguments
|
||
///
|
||
/// * `state` - Current system state (edge pressures and enthalpies)
|
||
/// * `residuals` - Residual vector to append constraint residuals to
|
||
/// * `measured_values` - Map of constraint IDs to their measured values (from component state)
|
||
///
|
||
/// # Returns
|
||
///
|
||
/// `Ok(count)` where count is the number of constraint residuals added.
|
||
///
|
||
/// # Errors
|
||
///
|
||
/// Returns `ConstraintError::UnmeasuredConstraint` if a constraint references a component
|
||
/// with no measured value (not registered or no associated edges).
|
||
/// Returns `ConstraintError::ResidualSliceTooShort` if the residual slice is too short.
|
||
///
|
||
/// # Example
|
||
///
|
||
/// ```rust,ignore
|
||
/// let mut residuals = ResidualVector::new();
|
||
/// let measured = system.extract_constraint_values_with_controls(&state, &control);
|
||
/// let count = system.compute_constraint_residuals(&state, &mut residuals, &measured);
|
||
/// ```
|
||
pub fn compute_constraint_residuals(
|
||
&self,
|
||
_state: &StateSlice,
|
||
residuals: &mut [f64],
|
||
measured_values: &HashMap<ConstraintId, f64>,
|
||
) -> Result<usize, ConstraintError> {
|
||
if self.constraints.is_empty() {
|
||
return Ok(0);
|
||
}
|
||
|
||
let mut count = 0;
|
||
for constraint in self.constraints.values() {
|
||
let measured = match measured_values.get(constraint.id()).copied() {
|
||
Some(v) => v,
|
||
None => {
|
||
return Err(ConstraintError::UnmeasuredConstraint {
|
||
constraint_id: constraint.id().to_string(),
|
||
component_id: constraint.output().component_id().to_string(),
|
||
});
|
||
}
|
||
};
|
||
let residual = constraint.compute_residual(measured);
|
||
if count >= residuals.len() {
|
||
return Err(ConstraintError::ResidualSliceTooShort {
|
||
index: count,
|
||
len: residuals.len(),
|
||
required: self.constraints.len(),
|
||
});
|
||
}
|
||
residuals[count] = residual;
|
||
count += 1;
|
||
}
|
||
Ok(count)
|
||
}
|
||
|
||
/// Extracts measured values for all constraints, incorporating control variable effects.
|
||
///
|
||
/// This method computes the measured output value for each constraint, taking into
|
||
/// account the current state and control variable values. For MIMO (Multi-Input
|
||
/// Multi-Output) systems, ALL control variables can affect ALL constraint outputs
|
||
/// due to system coupling.
|
||
///
|
||
/// # Arguments
|
||
///
|
||
/// * `state` - Current system state (edge pressures and enthalpies)
|
||
/// * `control_values` - Current values of control variables
|
||
///
|
||
/// # Returns
|
||
///
|
||
/// A map from constraint ID to measured output value.
|
||
///
|
||
/// # Cross-Coupling for MIMO Systems
|
||
///
|
||
/// In a real thermodynamic system, control variables are coupled:
|
||
/// - Compressor speed affects both capacity AND superheat
|
||
/// - Valve opening affects both superheat AND capacity
|
||
///
|
||
/// The mock implementation simulates this coupling for Jacobian cross-derivative
|
||
/// computation. Each control variable has a primary effect (on its linked constraint)
|
||
/// and a secondary effect (on other constraints) to simulate thermal coupling.
|
||
/// Maps a solver-side [`ComponentOutput`] to the components-crate
|
||
/// [`MeasuredOutput`] used by [`Component::measure_output`]. Returns `None`
|
||
/// for outputs that have no real-measurement mapping yet (falls back to the
|
||
/// legacy placeholder path).
|
||
fn map_output_kind(
|
||
output: &crate::inverse::ComponentOutput,
|
||
) -> Option<entropyk_components::MeasuredOutput> {
|
||
use crate::inverse::ComponentOutput as CO;
|
||
use entropyk_components::MeasuredOutput as MO;
|
||
Some(match output {
|
||
CO::Superheat { .. } => MO::Superheat,
|
||
CO::Subcooling { .. } => MO::Subcooling,
|
||
CO::Capacity { .. } => MO::Capacity,
|
||
CO::HeatTransferRate { .. } => MO::HeatTransferRate,
|
||
CO::MassFlowRate { .. } => MO::MassFlowRate,
|
||
CO::Pressure { .. } => MO::Pressure,
|
||
CO::Temperature { .. } => MO::Temperature,
|
||
CO::SaturationTemperature { .. } => MO::SaturationTemperature,
|
||
})
|
||
}
|
||
|
||
pub fn extract_constraint_values_with_controls(
|
||
&self,
|
||
state: &StateSlice,
|
||
control_values: &[f64],
|
||
) -> HashMap<ConstraintId, f64> {
|
||
let mut measured = HashMap::new();
|
||
if self.constraints.is_empty() {
|
||
return measured;
|
||
}
|
||
|
||
// Build a map of control variable index -> component_id it controls
|
||
// This uses the proper component_id() field from BoundedVariable
|
||
let mut control_to_component: HashMap<usize, &str> = HashMap::new();
|
||
for (j, bounded_var_id) in self.inverse_control.linked_controls().enumerate() {
|
||
if let Some(bounded_var) = self.bounded_variables.get(bounded_var_id) {
|
||
if let Some(comp_id) = bounded_var.component_id() {
|
||
control_to_component.insert(j, comp_id);
|
||
}
|
||
}
|
||
}
|
||
|
||
for constraint in self.constraints.values() {
|
||
let comp_id = constraint.output().component_id();
|
||
if let Some(&node_idx) = self.component_names.get(comp_id) {
|
||
// Prefer a REAL physical measurement from the component itself.
|
||
// This replaces the placeholder formulas below with genuine
|
||
// thermodynamics (superheat, subcooling, capacity, …). The mock
|
||
// path is only used for components/outputs that do not (yet)
|
||
// implement `measure_output` (e.g. synthetic test doubles).
|
||
if let Some(kind) = Self::map_output_kind(constraint.output()) {
|
||
if let Some(component) = self.graph.node_weight(node_idx) {
|
||
if let Some(real) = component.measure_output(kind, state) {
|
||
if real.is_finite() {
|
||
measured.insert(constraint.id().clone(), real);
|
||
continue;
|
||
}
|
||
}
|
||
}
|
||
}
|
||
|
||
// Find first associated edge (incoming or outgoing)
|
||
let mut edge_opt = self
|
||
.graph
|
||
.edges_directed(node_idx, petgraph::Direction::Incoming)
|
||
.next();
|
||
if edge_opt.is_none() {
|
||
edge_opt = self
|
||
.graph
|
||
.edges_directed(node_idx, petgraph::Direction::Outgoing)
|
||
.next();
|
||
}
|
||
|
||
if let Some(edge) = edge_opt {
|
||
if let Some(&(_m, p_idx, h_idx)) = self.edge_to_state.get(&edge.id()) {
|
||
let mut value = match constraint.output() {
|
||
crate::inverse::ComponentOutput::Pressure { .. } => state[p_idx],
|
||
crate::inverse::ComponentOutput::Temperature { .. } => 300.0, // Mock for MVP without fluid backend
|
||
crate::inverse::ComponentOutput::Superheat { .. } => {
|
||
// Mock numerical value sensitive to BOTH P and h for Jacobian calculation
|
||
state[h_idx] / 1000.0 - (state[p_idx] / 1e5)
|
||
}
|
||
crate::inverse::ComponentOutput::Subcooling { .. } => {
|
||
(state[p_idx] / 1e5) - state[h_idx] / 1000.0
|
||
}
|
||
crate::inverse::ComponentOutput::Capacity { .. } => {
|
||
// Mock capacity: h * mass_flow. Let's just use h for Jacobian sensitivity
|
||
state[h_idx] * 10.0
|
||
}
|
||
_ => 0.0,
|
||
};
|
||
|
||
// MIMO Cross-Coupling: ALL control variables can affect ALL constraints
|
||
// In a real system, changing compressor speed affects both capacity and superheat,
|
||
// and changing valve opening also affects both. We simulate this coupling here.
|
||
//
|
||
// ⚠️ MOCK COEFFICIENTS: These values (10.0, 2.0) are placeholders for testing.
|
||
// They create a well-conditioned Jacobian with off-diagonal entries that allow
|
||
// Newton-Raphson to converge. Real implementations should replace these with
|
||
// actual component physics derived from:
|
||
// - Component characteristic curves (compressor map, valve Cv curve)
|
||
// - Thermodynamic property calculations via fluid backend
|
||
// - Energy and mass balance equations
|
||
//
|
||
// The 5:1 ratio between primary and secondary effects is arbitrary but creates
|
||
// a diagonally-dominant Jacobian that converges reliably. See Story 5.4
|
||
// Review Follow-ups for tracking real thermodynamics integration.
|
||
//
|
||
// For each control variable:
|
||
// - Primary effect (10.0): if control is linked to this constraint's component
|
||
// - Secondary effect (2.0): cross-coupling to other constraints
|
||
const MIMO_PRIMARY_COEFF: f64 = 10.0;
|
||
const MIMO_SECONDARY_COEFF: f64 = 2.0;
|
||
|
||
for (j, _bounded_var_id) in
|
||
self.inverse_control.linked_controls().enumerate()
|
||
{
|
||
if j >= control_values.len() {
|
||
continue;
|
||
}
|
||
let ctrl_val = control_values[j];
|
||
|
||
// Check if this control variable is primarily associated with this component
|
||
let is_primary = control_to_component
|
||
.get(&j)
|
||
.map_or(false, |&c| c == comp_id);
|
||
|
||
if is_primary {
|
||
// Primary effect: strong influence on the controlled output
|
||
// e.g., valve opening strongly affects superheat
|
||
value += ctrl_val * MIMO_PRIMARY_COEFF;
|
||
} else {
|
||
// Secondary (cross-coupling) effect: weaker influence
|
||
// e.g., compressor speed also affects superheat (through mass flow)
|
||
// This creates the off-diagonal entries in the MIMO Jacobian
|
||
value += ctrl_val * MIMO_SECONDARY_COEFF;
|
||
}
|
||
}
|
||
|
||
if value.is_nan() {
|
||
tracing::warn!(
|
||
constraint_id = constraint.id().as_str(),
|
||
"NaN detected in constraint output for component '{}', skipping insert",
|
||
constraint.output().component_id()
|
||
);
|
||
} else {
|
||
measured.insert(constraint.id().clone(), value);
|
||
}
|
||
}
|
||
}
|
||
}
|
||
}
|
||
measured
|
||
}
|
||
|
||
/// Computes the Jacobian entries for inverse control constraints.
|
||
///
|
||
/// For each constraint→control mapping, adds ∂r/∂x entries to the Jacobian:
|
||
///
|
||
/// $$\frac{\partial r}{\partial x_{control}} = \frac{\partial (measured - target)}{\partial x_{control}} = \frac{\partial measured}{\partial x_{control}}$$
|
||
///
|
||
/// # Arguments
|
||
///
|
||
/// * `state` - Current system state
|
||
/// * `row_offset` - Starting row index for constraint equations in the Jacobian
|
||
/// * `control_values` - Current values of control variables (for finite difference)
|
||
///
|
||
/// # Returns
|
||
///
|
||
/// A vector of `(row, col, value)` tuples representing Jacobian entries.
|
||
///
|
||
/// # Note
|
||
///
|
||
/// MVP uses finite difference approximation. Future versions may use analytical
|
||
/// derivatives from components for better accuracy and performance.
|
||
///
|
||
/// # Finite Difference Epsilon
|
||
///
|
||
/// Uses the epsilon configured in `InverseControlConfig` (default 1e-6) for central
|
||
/// finite differences. Configure via `set_inverse_control_epsilon()`.
|
||
pub fn compute_inverse_control_jacobian(
|
||
&self,
|
||
state: &StateSlice,
|
||
row_offset: usize,
|
||
control_values: &[f64],
|
||
) -> Vec<(usize, usize, f64)> {
|
||
let mut entries = Vec::new();
|
||
|
||
if self.inverse_control.mapping_count() == 0 {
|
||
return entries;
|
||
}
|
||
|
||
if control_values.len() < self.inverse_control.mapping_count() {
|
||
tracing::error!(
|
||
provided = control_values.len(),
|
||
required = self.inverse_control.mapping_count(),
|
||
"control_values too short for Jacobian computation"
|
||
);
|
||
return entries;
|
||
}
|
||
|
||
// Use configurable epsilon from InverseControlConfig
|
||
let eps = self.inverse_control.finite_diff_epsilon();
|
||
if state.len() < self.total_state_len {
|
||
tracing::error!(
|
||
state_len = state.len(),
|
||
required = self.total_state_len,
|
||
"compute_inverse_control_jacobian: state slice too short, returning empty"
|
||
);
|
||
return entries;
|
||
}
|
||
let mut state_mut = state.to_vec();
|
||
let mut control_mut = control_values.to_vec();
|
||
|
||
// 1. Compute ∂r_i / ∂x_j (Partial derivatives with respect to PHYSICAL states P, h)
|
||
// We do this per constraint to keep perturbations localized where possible
|
||
for (i, (constraint_id, _)) in self.inverse_control.mappings().enumerate() {
|
||
let row = row_offset + i;
|
||
if let Some(constraint) = self.constraints.get(constraint_id) {
|
||
let comp_id = constraint.output().component_id();
|
||
|
||
if let Some(&node_idx) = self.component_names.get(comp_id) {
|
||
let mut state_indices = Vec::new();
|
||
// Gather all edge state indices for this component. The
|
||
// mass-flow index is included because measured outputs such
|
||
// as Capacity (q = ṁ·Δh) depend directly on ṁ; omitting it
|
||
// decouples the constraint from any mass-flow actuator (f_m)
|
||
// and makes the inverse-control Jacobian singular.
|
||
for edge in self
|
||
.graph
|
||
.edges_directed(node_idx, petgraph::Direction::Incoming)
|
||
{
|
||
if let Some(&(m_idx, p_idx, h_idx)) = self.edge_to_state.get(&edge.id()) {
|
||
if !state_indices.contains(&m_idx) {
|
||
state_indices.push(m_idx);
|
||
}
|
||
if !state_indices.contains(&p_idx) {
|
||
state_indices.push(p_idx);
|
||
}
|
||
if !state_indices.contains(&h_idx) {
|
||
state_indices.push(h_idx);
|
||
}
|
||
}
|
||
}
|
||
for edge in self
|
||
.graph
|
||
.edges_directed(node_idx, petgraph::Direction::Outgoing)
|
||
{
|
||
if let Some(&(m_idx, p_idx, h_idx)) = self.edge_to_state.get(&edge.id()) {
|
||
if !state_indices.contains(&m_idx) {
|
||
state_indices.push(m_idx);
|
||
}
|
||
if !state_indices.contains(&p_idx) {
|
||
state_indices.push(p_idx);
|
||
}
|
||
if !state_indices.contains(&h_idx) {
|
||
state_indices.push(h_idx);
|
||
}
|
||
}
|
||
}
|
||
|
||
// Central finite difference for Jacobian entries w.r.t physical state
|
||
for &col in &state_indices {
|
||
let orig = state_mut[col];
|
||
|
||
state_mut[col] = orig + eps;
|
||
let plus = self
|
||
.extract_constraint_values_with_controls(&state_mut, control_values);
|
||
let val_plus = plus.get(constraint_id).copied().unwrap_or(0.0);
|
||
|
||
state_mut[col] = orig - eps;
|
||
let minus = self
|
||
.extract_constraint_values_with_controls(&state_mut, control_values);
|
||
let val_minus = minus.get(constraint_id).copied().unwrap_or(0.0);
|
||
|
||
state_mut[col] = orig; // Restore
|
||
|
||
let derivative = (val_plus - val_minus) / (2.0 * eps);
|
||
if derivative.abs() > 1e-10 {
|
||
entries.push((row, col, derivative));
|
||
tracing::trace!(
|
||
constraint = constraint_id.as_str(),
|
||
row,
|
||
col,
|
||
derivative,
|
||
"Inverse control Jacobian actual ∂r/∂state entry"
|
||
);
|
||
}
|
||
}
|
||
}
|
||
}
|
||
}
|
||
|
||
// 2. Compute ∂r_i / ∂u_j (Cross-derivatives with respect to CONTROL variables)
|
||
// Here we must form the full dense block because control variable 'j' could affect constraint 'i'
|
||
// even if they are not explicitly linked, due to system coupling.
|
||
let control_offset = self.state_vector_len();
|
||
|
||
for (j, (_, bounded_var_id)) in self.inverse_control.mappings().enumerate() {
|
||
let col = control_offset + j;
|
||
let orig = control_mut[j];
|
||
|
||
// Perturb control variable +eps
|
||
control_mut[j] = orig + eps;
|
||
let plus = self.extract_constraint_values_with_controls(state, &control_mut);
|
||
|
||
// Perturb control variable -eps
|
||
control_mut[j] = orig - eps;
|
||
let minus = self.extract_constraint_values_with_controls(state, &control_mut);
|
||
|
||
control_mut[j] = orig; // Restore
|
||
|
||
// For this perturbed control variable j, compute the effect on ALL constraints i
|
||
for (i, (constraint_id, _)) in self.inverse_control.mappings().enumerate() {
|
||
let row = row_offset + i;
|
||
|
||
let val_plus = plus.get(constraint_id).copied().unwrap_or(0.0);
|
||
let val_minus = minus.get(constraint_id).copied().unwrap_or(0.0);
|
||
let derivative = (val_plus - val_minus) / (2.0 * eps);
|
||
|
||
// We add it even if it's 0 to maintain block structure (optional but safe)
|
||
// However, for performance we only add non-zeros
|
||
if derivative.abs() > 1e-10 {
|
||
entries.push((row, col, derivative));
|
||
tracing::trace!(
|
||
constraint = ?constraint_id,
|
||
control = ?bounded_var_id,
|
||
row, col, derivative,
|
||
"Inverse control Jacobian cross-derivative ∂r/∂u entry"
|
||
);
|
||
}
|
||
}
|
||
}
|
||
|
||
entries
|
||
}
|
||
|
||
fn incident_state_indices_for_component(&self, node_idx: NodeIndex) -> Vec<usize> {
|
||
let mut state_indices = Vec::new();
|
||
for edge in self
|
||
.graph
|
||
.edges_directed(node_idx, petgraph::Direction::Incoming)
|
||
{
|
||
if let Some(&(m_idx, p_idx, h_idx)) = self.edge_to_state.get(&edge.id()) {
|
||
for idx in [m_idx, p_idx, h_idx] {
|
||
if !state_indices.contains(&idx) {
|
||
state_indices.push(idx);
|
||
}
|
||
}
|
||
}
|
||
}
|
||
for edge in self
|
||
.graph
|
||
.edges_directed(node_idx, petgraph::Direction::Outgoing)
|
||
{
|
||
if let Some(&(m_idx, p_idx, h_idx)) = self.edge_to_state.get(&edge.id()) {
|
||
for idx in [m_idx, p_idx, h_idx] {
|
||
if !state_indices.contains(&idx) {
|
||
state_indices.push(idx);
|
||
}
|
||
}
|
||
}
|
||
}
|
||
state_indices
|
||
}
|
||
|
||
fn measure_saturated_output(
|
||
&self,
|
||
output: &crate::inverse::ComponentOutput,
|
||
state: &StateSlice,
|
||
) -> f64 {
|
||
let comp_id = output.component_id();
|
||
if let Some(&node_idx) = self.component_names.get(comp_id) {
|
||
if let Some(kind) = Self::map_output_kind(output) {
|
||
if let Some(component) = self.graph.node_weight(node_idx) {
|
||
if let Some(value) = component.measure_output(kind, state) {
|
||
if value.is_finite() {
|
||
return value;
|
||
}
|
||
}
|
||
}
|
||
}
|
||
}
|
||
|
||
tracing::warn!(
|
||
component_id = comp_id,
|
||
output = output.constraint_type_name(),
|
||
"Saturated controller output measurement unavailable or non-finite; using synthetic 0.0 fallback"
|
||
);
|
||
0.0
|
||
}
|
||
|
||
fn compute_saturated_control_jacobian(
|
||
&self,
|
||
state: &StateSlice,
|
||
row_offset: usize,
|
||
) -> Vec<(usize, usize, f64)> {
|
||
let mut entries = Vec::new();
|
||
if self.saturated_controllers.is_empty() {
|
||
return entries;
|
||
}
|
||
|
||
let eps = self.inverse_control.finite_diff_epsilon();
|
||
let mut state_mut = state.to_vec();
|
||
|
||
for (i, ctrl) in self.saturated_controllers.iter().enumerate() {
|
||
let u_idx = self.saturated_u_index(i);
|
||
let x_idx = self.saturated_x_index(i);
|
||
if state.len() <= x_idx {
|
||
tracing::error!(
|
||
state_len = state.len(),
|
||
required = x_idx + 1,
|
||
"compute_saturated_control_jacobian: state slice too short"
|
||
);
|
||
continue;
|
||
}
|
||
|
||
let row_u = row_offset + 2 * i;
|
||
let row_y = row_u + 1;
|
||
let x = state[x_idx];
|
||
let (dru_du, dru_dx) = ctrl.d_residual_u(x);
|
||
let (_dry_dy, dry_dx) = ctrl.d_residual_y(x);
|
||
entries.push((row_u, u_idx, dru_du));
|
||
entries.push((row_u, x_idx, dru_dx));
|
||
entries.push((row_y, x_idx, dry_dx));
|
||
|
||
if ctrl.is_network() {
|
||
// Override network: r_y = E − (x − S(x)); the plant coupling is
|
||
// ∂E/∂col = Σ_i w_i·(−gain_i)·(∂measurement_i/∂col), where w_i are
|
||
// the analytic softMin/softMax selector weights.
|
||
let objs = ctrl.objectives();
|
||
let measured: Vec<f64> = objs
|
||
.iter()
|
||
.map(|o| self.measure_saturated_output(&o.output, &state_mut))
|
||
.collect();
|
||
let weights = ctrl.network_error_weights(&measured);
|
||
let mut col_acc: HashMap<usize, f64> = HashMap::new();
|
||
for (oi, o) in objs.iter().enumerate() {
|
||
if let Some(&node_idx) = self.component_names.get(o.output.component_id()) {
|
||
for col in self.incident_state_indices_for_component(node_idx) {
|
||
let orig = state_mut[col];
|
||
state_mut[col] = orig + eps;
|
||
let y_plus = self.measure_saturated_output(&o.output, &state_mut);
|
||
state_mut[col] = orig - eps;
|
||
let y_minus = self.measure_saturated_output(&o.output, &state_mut);
|
||
state_mut[col] = orig;
|
||
let dy_dcol = (y_plus - y_minus) / (2.0 * eps);
|
||
*col_acc.entry(col).or_insert(0.0) += weights[oi] * (-o.gain) * dy_dcol;
|
||
}
|
||
}
|
||
}
|
||
for (col, derivative) in col_acc {
|
||
if derivative.abs() > 1e-10 {
|
||
entries.push((row_y, col, derivative));
|
||
}
|
||
}
|
||
} else {
|
||
let comp_id = ctrl.output().component_id();
|
||
if let Some(&node_idx) = self.component_names.get(comp_id) {
|
||
for col in self.incident_state_indices_for_component(node_idx) {
|
||
let orig = state_mut[col];
|
||
state_mut[col] = orig + eps;
|
||
let y_plus = self.measure_saturated_output(ctrl.output(), &state_mut);
|
||
state_mut[col] = orig - eps;
|
||
let y_minus = self.measure_saturated_output(ctrl.output(), &state_mut);
|
||
state_mut[col] = orig;
|
||
|
||
let dy_dcol = (y_plus - y_minus) / (2.0 * eps);
|
||
let derivative = -ctrl.gain() * dy_dcol;
|
||
if derivative.abs() > 1e-10 {
|
||
entries.push((row_y, col, derivative));
|
||
}
|
||
}
|
||
}
|
||
}
|
||
}
|
||
|
||
entries
|
||
}
|
||
|
||
// ────────────────────────────────────────────────────────────────────────
|
||
// Bounded Variable Management (Inverse Control)
|
||
// ────────────────────────────────────────────────────────────────────────
|
||
|
||
/// Adds a bounded control variable for inverse control.
|
||
///
|
||
/// Bounded variables ensure Newton steps stay within physical limits
|
||
/// (e.g., valve position 0.0 to 1.0, VFD speed 0.3 to 1.0).
|
||
///
|
||
/// # Arguments
|
||
///
|
||
/// * `variable` - The bounded variable to add
|
||
///
|
||
/// # Errors
|
||
///
|
||
/// Returns `BoundedVariableError::DuplicateId` if a variable with the same ID
|
||
/// already exists.
|
||
/// Returns `BoundedVariableError::InvalidComponent` if the variable is associated
|
||
/// with a component that has not been registered via `register_component_name()`.
|
||
///
|
||
/// # Example
|
||
///
|
||
/// ```rust,ignore
|
||
/// use entropyk_solver::inverse::{BoundedVariable, BoundedVariableId};
|
||
///
|
||
/// let valve = BoundedVariable::new(
|
||
/// BoundedVariableId::new("expansion_valve"),
|
||
/// 0.5, // initial: 50% open
|
||
/// 0.0, // min: fully closed
|
||
/// 1.0, // max: fully open
|
||
/// )?;
|
||
///
|
||
/// system.add_bounded_variable(valve)?;
|
||
/// ```
|
||
pub fn add_bounded_variable(
|
||
&mut self,
|
||
variable: BoundedVariable,
|
||
) -> Result<(), BoundedVariableError> {
|
||
let id = variable.id().clone();
|
||
if self.bounded_variables.contains_key(&id) {
|
||
return Err(BoundedVariableError::DuplicateId { id });
|
||
}
|
||
|
||
// Validate that the component referenced in the variable exists (if any)
|
||
if let Some(component_id) = variable.component_id() {
|
||
if !self.component_names.contains_key(component_id) {
|
||
return Err(BoundedVariableError::InvalidComponent {
|
||
component_id: component_id.to_string(),
|
||
});
|
||
}
|
||
}
|
||
|
||
self.bounded_variables.insert(id, variable);
|
||
Ok(())
|
||
}
|
||
|
||
/// Removes a bounded variable by ID.
|
||
///
|
||
/// # Arguments
|
||
///
|
||
/// * `id` - The bounded variable identifier to remove
|
||
///
|
||
/// # Returns
|
||
///
|
||
/// The removed variable, or `None` if no variable with that ID exists.
|
||
pub fn remove_bounded_variable(&mut self, id: &BoundedVariableId) -> Option<BoundedVariable> {
|
||
self.inverse_control.unlink_control(id);
|
||
self.free_actuators.retain(|a| a != id);
|
||
self.bounded_variables.remove(id)
|
||
}
|
||
|
||
/// Returns the number of bounded variables in the system.
|
||
pub fn bounded_variable_count(&self) -> usize {
|
||
self.bounded_variables.len()
|
||
}
|
||
|
||
/// Returns an iterator over all bounded variables.
|
||
pub fn bounded_variables(&self) -> impl Iterator<Item = &BoundedVariable> {
|
||
self.bounded_variables.values()
|
||
}
|
||
|
||
/// Returns a reference to a specific bounded variable by ID.
|
||
pub fn get_bounded_variable(&self, id: &BoundedVariableId) -> Option<&BoundedVariable> {
|
||
self.bounded_variables.get(id)
|
||
}
|
||
|
||
/// Returns a mutable reference to a specific bounded variable by ID.
|
||
pub fn get_bounded_variable_mut(
|
||
&mut self,
|
||
id: &BoundedVariableId,
|
||
) -> Option<&mut BoundedVariable> {
|
||
self.bounded_variables.get_mut(id)
|
||
}
|
||
|
||
/// Checks if any bounded variables are saturated (at bounds).
|
||
///
|
||
/// Returns a vector of saturation info for all variables currently at bounds.
|
||
pub fn saturated_variables(&self) -> Vec<crate::inverse::SaturationInfo> {
|
||
self.bounded_variables
|
||
.iter()
|
||
.filter(|(id, _)| !self.free_actuators.contains(id))
|
||
.filter_map(|(_, v)| v.is_saturated())
|
||
.collect()
|
||
}
|
||
|
||
// ────────────────────────────────────────────────────────────────────────
|
||
// Inverse Control Mapping (Story 5.3)
|
||
// ────────────────────────────────────────────────────────────────────────
|
||
|
||
/// Removes all constraints, bounded variables, and inverse control mappings.
|
||
pub fn clear_inverse_control(&mut self) {
|
||
self.constraints.clear();
|
||
self.bounded_variables.clear();
|
||
self.inverse_control.clear();
|
||
self.saturated_controllers.clear();
|
||
self.free_actuators.clear();
|
||
}
|
||
|
||
/// Adds a saturated PI controller loop.
|
||
///
|
||
/// The actuator bounded variable must already be registered and must not be
|
||
/// linked to a hard inverse-control constraint. The saturated controller
|
||
/// contributes its own two residuals and two unknowns.
|
||
pub fn add_saturated_controller(&mut self, ctrl: SaturatedController) {
|
||
assert!(
|
||
self.bounded_variables.contains_key(ctrl.actuator()),
|
||
"saturated controller actuator '{}' must be registered as a bounded variable",
|
||
ctrl.actuator()
|
||
);
|
||
assert!(
|
||
!self.inverse_control.is_control_linked(ctrl.actuator()),
|
||
"saturated controller actuator '{}' must not be linked as a hard constraint",
|
||
ctrl.actuator()
|
||
);
|
||
self.saturated_controllers.push(ctrl);
|
||
}
|
||
|
||
/// Registers a physical free-actuator unknown (arch-6).
|
||
///
|
||
/// The actuator bounded variable must already be registered via
|
||
/// [`add_bounded_variable`](Self::add_bounded_variable), must reference an
|
||
/// existing component, and must not be linked to a hard inverse-control
|
||
/// constraint or driven by a saturated controller. It appends exactly one
|
||
/// unknown to the state vector; the owning component must supply exactly one
|
||
/// extra equation (its `n_equations()` grows by one) so the DoF stays
|
||
/// balanced. The actuator value is wired into the component's generic
|
||
/// `CalibIndices::actuator` slot during [`finalize`](Self::finalize).
|
||
pub fn add_free_actuator(&mut self, actuator_id: BoundedVariableId) {
|
||
assert!(
|
||
self.bounded_variables.contains_key(&actuator_id),
|
||
"free actuator '{actuator_id}' must be registered as a bounded variable"
|
||
);
|
||
assert!(
|
||
!self.inverse_control.is_control_linked(&actuator_id),
|
||
"free actuator '{actuator_id}' must not be linked as a hard constraint"
|
||
);
|
||
assert!(
|
||
!self.free_actuators.contains(&actuator_id),
|
||
"free actuator '{actuator_id}' is already registered"
|
||
);
|
||
self.free_actuators.push(actuator_id);
|
||
}
|
||
|
||
/// Returns the number of physical free-actuator unknowns.
|
||
pub fn free_actuator_count(&self) -> usize {
|
||
self.free_actuators.len()
|
||
}
|
||
|
||
/// Iterates over the registered free-actuator bounded-variable ids in
|
||
/// insertion (state-layout) order.
|
||
pub fn free_actuators(&self) -> impl Iterator<Item = &BoundedVariableId> {
|
||
self.free_actuators.iter()
|
||
}
|
||
|
||
/// Returns the number of saturated PI controller loops.
|
||
pub fn saturated_controller_count(&self) -> usize {
|
||
self.saturated_controllers.len()
|
||
}
|
||
|
||
/// Iterates over saturated PI controller loops in insertion order.
|
||
pub fn saturated_controllers(&self) -> impl Iterator<Item = &SaturatedController> {
|
||
self.saturated_controllers.iter()
|
||
}
|
||
|
||
/// Mutable iterator over saturated PI controllers in insertion order.
|
||
///
|
||
/// Exposed so an outer driver can retune controllers between solves — in
|
||
/// particular the warm-started **alpha-continuation** that anneals every
|
||
/// override-network selector sharpness toward its target.
|
||
pub fn saturated_controllers_mut(&mut self) -> impl Iterator<Item = &mut SaturatedController> {
|
||
self.saturated_controllers.iter_mut()
|
||
}
|
||
|
||
/// Whether any saturated controller uses an override/selector network.
|
||
pub fn has_network_controllers(&self) -> bool {
|
||
self.saturated_controllers.iter().any(|c| c.is_network())
|
||
}
|
||
|
||
/// Links a constraint to a bounded control variable for One-Shot inverse control.
|
||
///
|
||
/// When a constraint is linked to a control variable, the solver adjusts both
|
||
/// the edge states AND the control variable simultaneously to satisfy all
|
||
/// equations (cycle equations + constraint equations).
|
||
///
|
||
/// # Arguments
|
||
///
|
||
/// * `constraint_id` - The constraint to link
|
||
/// * `bounded_variable_id` - The control variable to adjust
|
||
///
|
||
/// # Errors
|
||
///
|
||
/// Returns `DoFError::ConstraintNotFound` if the constraint doesn't exist.
|
||
/// Returns `DoFError::BoundedVariableNotFound` if the bounded variable doesn't exist.
|
||
/// Returns `DoFError::AlreadyLinked` if the constraint is already linked.
|
||
/// Returns `DoFError::ControlAlreadyLinked` if the control is already linked.
|
||
///
|
||
/// # Example
|
||
///
|
||
/// ```rust,ignore
|
||
/// system.link_constraint_to_control(
|
||
/// &ConstraintId::new("superheat_control"),
|
||
/// &BoundedVariableId::new("expansion_valve"),
|
||
/// )?;
|
||
/// ```
|
||
pub fn link_constraint_to_control(
|
||
&mut self,
|
||
constraint_id: &ConstraintId,
|
||
bounded_variable_id: &BoundedVariableId,
|
||
) -> Result<(), DoFError> {
|
||
if !self.constraints.contains_key(constraint_id) {
|
||
return Err(DoFError::ConstraintNotFound {
|
||
constraint_id: constraint_id.clone(),
|
||
});
|
||
}
|
||
|
||
if !self.bounded_variables.contains_key(bounded_variable_id) {
|
||
return Err(DoFError::BoundedVariableNotFound {
|
||
bounded_variable_id: bounded_variable_id.clone(),
|
||
});
|
||
}
|
||
|
||
self.inverse_control
|
||
.link(constraint_id.clone(), bounded_variable_id.clone())
|
||
}
|
||
|
||
/// Unlinks a constraint from its control variable.
|
||
///
|
||
/// # Arguments
|
||
///
|
||
/// * `constraint_id` - The constraint to unlink
|
||
///
|
||
/// # Returns
|
||
///
|
||
/// The bounded variable ID that was linked, or `None` if not linked.
|
||
pub fn unlink_constraint(&mut self, constraint_id: &ConstraintId) -> Option<BoundedVariableId> {
|
||
self.inverse_control.unlink_constraint(constraint_id)
|
||
}
|
||
|
||
/// Unlinks a control variable from its constraint.
|
||
///
|
||
/// # Arguments
|
||
///
|
||
/// * `bounded_variable_id` - The control variable to unlink
|
||
///
|
||
/// # Returns
|
||
///
|
||
/// The constraint ID that was linked, or `None` if not linked.
|
||
pub fn unlink_control(
|
||
&mut self,
|
||
bounded_variable_id: &BoundedVariableId,
|
||
) -> Option<ConstraintId> {
|
||
self.inverse_control.unlink_control(bounded_variable_id)
|
||
}
|
||
|
||
/// Returns the control variable linked to a constraint.
|
||
pub fn get_control_for_constraint(
|
||
&self,
|
||
constraint_id: &ConstraintId,
|
||
) -> Option<&BoundedVariableId> {
|
||
self.inverse_control.get_control(constraint_id)
|
||
}
|
||
|
||
/// Returns the constraint linked to a control variable.
|
||
pub fn get_constraint_for_control(
|
||
&self,
|
||
bounded_variable_id: &BoundedVariableId,
|
||
) -> Option<&ConstraintId> {
|
||
self.inverse_control.get_constraint(bounded_variable_id)
|
||
}
|
||
|
||
/// Returns the number of constraint-to-control mappings.
|
||
pub fn inverse_control_mapping_count(&self) -> usize {
|
||
self.inverse_control.mapping_count()
|
||
}
|
||
|
||
/// Sets the finite difference epsilon for inverse control Jacobian computation.
|
||
///
|
||
/// # Errors
|
||
///
|
||
/// Returns `ConstraintError::InvalidEpsilon` if epsilon is not a finite positive value in (0, 1].
|
||
pub fn set_inverse_control_epsilon(&mut self, epsilon: f64) -> Result<(), ConstraintError> {
|
||
self.inverse_control.set_finite_diff_epsilon(epsilon)
|
||
}
|
||
|
||
/// Returns the current finite difference epsilon for inverse control.
|
||
pub fn inverse_control_epsilon(&self) -> f64 {
|
||
self.inverse_control.finite_diff_epsilon()
|
||
}
|
||
|
||
/// Returns an iterator over linked control variable IDs.
|
||
pub fn linked_controls(&self) -> impl Iterator<Item = &BoundedVariableId> {
|
||
self.inverse_control.linked_controls()
|
||
}
|
||
|
||
/// Checks if a constraint is linked to a control variable.
|
||
pub fn is_constraint_linked(&self, constraint_id: &ConstraintId) -> bool {
|
||
self.inverse_control.is_constraint_linked(constraint_id)
|
||
}
|
||
|
||
/// Checks if a control variable is linked to a constraint.
|
||
pub fn is_control_linked(&self, bounded_variable_id: &BoundedVariableId) -> bool {
|
||
self.inverse_control.is_control_linked(bounded_variable_id)
|
||
}
|
||
|
||
/// Validates degrees of freedom for inverse control.
|
||
///
|
||
/// For a well-posed system with inverse control:
|
||
///
|
||
/// $$n_{equations} = n_{edge\_eqs} + n_{constraints}$$
|
||
/// $$n_{unknowns} = n_{edge\_unknowns} + n_{controls}$$
|
||
///
|
||
/// The system is balanced when: $n_{equations} = n_{unknowns}$
|
||
///
|
||
/// # Returns
|
||
///
|
||
/// `Ok(())` if the system is well-posed (balanced or under-constrained).
|
||
///
|
||
/// # Errors
|
||
///
|
||
/// Returns `DoFError::OverConstrainedSystem` if there are more equations than unknowns.
|
||
/// Returns `DoFError::UnderConstrainedSystem` if there are fewer equations than unknowns (warning only).
|
||
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_unknowns = n_edge_unknowns + n_controls;
|
||
|
||
let n_edge_eqs: usize = self
|
||
.graph
|
||
.node_indices()
|
||
.map(|node| {
|
||
self.graph
|
||
.node_weight(node)
|
||
.map(|c| c.n_equations())
|
||
.unwrap_or(0)
|
||
})
|
||
.sum();
|
||
let n_equations = n_edge_eqs + n_constraints;
|
||
|
||
if n_equations > n_unknowns {
|
||
Err(DoFError::OverConstrainedSystem {
|
||
constraint_count: n_constraints,
|
||
control_count: n_controls,
|
||
equation_count: n_equations,
|
||
unknown_count: n_unknowns,
|
||
})
|
||
} else if n_equations < n_unknowns {
|
||
Err(DoFError::UnderConstrainedSystem {
|
||
constraint_count: n_constraints,
|
||
control_count: n_controls,
|
||
equation_count: n_equations,
|
||
unknown_count: n_unknowns,
|
||
})
|
||
} else {
|
||
Ok(())
|
||
}
|
||
}
|
||
|
||
/// Total residual equations assembled by [`Self::compute_residuals`].
|
||
///
|
||
/// Includes every component block, inverse constraints, thermal couplings,
|
||
/// and saturated-controller residual pairs.
|
||
pub fn total_equation_count(&self) -> usize {
|
||
let n_comp: usize = self
|
||
.graph
|
||
.node_indices()
|
||
.map(|node| {
|
||
self.graph
|
||
.node_weight(node)
|
||
.map(|c| c.n_equations())
|
||
.unwrap_or(0)
|
||
})
|
||
.sum();
|
||
n_comp
|
||
+ self.constraints.len()
|
||
+ self.coupling_residual_count()
|
||
+ 2 * self.saturated_controllers.len()
|
||
}
|
||
|
||
/// Builds a full DoF ledger for the finalized system.
|
||
///
|
||
/// Counts **all** Newton unknowns (`full_state_vector_len`) against **all**
|
||
/// residual equations (`total_equation_count`), including free actuators and
|
||
/// thermal couplings that the legacy inverse-only check omitted.
|
||
///
|
||
/// # Panics
|
||
///
|
||
/// Panics if `finalize()` has not been called (state length is undefined).
|
||
pub fn dof_report(&self) -> DofReport {
|
||
assert!(self.finalized, "call finalize() before dof_report()");
|
||
|
||
let mut diagnostics = Vec::new();
|
||
let mut components = Vec::new();
|
||
let mut n_comp_eqs = 0usize;
|
||
|
||
// Reverse name map for diagnostics.
|
||
let mut names_by_node: HashMap<NodeIndex, String> = HashMap::new();
|
||
for (name, idx) in &self.component_names {
|
||
names_by_node.insert(*idx, name.clone());
|
||
}
|
||
|
||
for node in self.graph.node_indices() {
|
||
let Some(component) = self.graph.node_weight(node) else {
|
||
continue;
|
||
};
|
||
let n = component.n_equations();
|
||
n_comp_eqs += n;
|
||
let name = names_by_node
|
||
.get(&node)
|
||
.cloned()
|
||
.unwrap_or_else(|| format!("node_{}", node.index()));
|
||
let roles = align_roles(&name, n, component.equation_roles(), &mut diagnostics);
|
||
|
||
// Heuristic: outlet closures that *consume* a DoF without freeing
|
||
// another (e.g. quality_control, Anchor specs) need an actuator.
|
||
// Emergent-cycle closures (subcooling / superheat / saturated_vapor)
|
||
// replace a fixed P elsewhere — they are the closed-cycle design, not
|
||
// an over-constraint.
|
||
let suspicious_closure = roles.iter().any(|r| match r {
|
||
EquationRole::OutletClosure { kind } => !matches!(
|
||
kind.as_ref(),
|
||
"subcooling" | "superheat" | "saturated_vapor"
|
||
),
|
||
_ => false,
|
||
});
|
||
if suspicious_closure
|
||
&& self.free_actuators.is_empty()
|
||
&& self.saturated_controllers.is_empty()
|
||
&& self.inverse_control.mapping_count() == 0
|
||
{
|
||
diagnostics.push(format!(
|
||
"component `{name}` declares OutletClosure residual(s) but the system \
|
||
has no free actuator / saturated controller / inverse control — \
|
||
fixing this quantity without freeing another unknown over-constrains \
|
||
a closed cycle (pair with EXV opening, speed, level, or drop the residual)"
|
||
));
|
||
}
|
||
|
||
components.push(ComponentEquationBlock {
|
||
component_name: name,
|
||
node_index: node.index(),
|
||
n_equations: n,
|
||
roles,
|
||
});
|
||
}
|
||
|
||
let mut system_equations = Vec::new();
|
||
for id in self.constraints.keys() {
|
||
system_equations.push(EquationRole::ControlTracking {
|
||
name: id.to_string(),
|
||
});
|
||
}
|
||
for _ in 0..self.coupling_residual_count() {
|
||
system_equations.push(EquationRole::CouplingDuty);
|
||
}
|
||
for _ in 0..self.saturated_controllers.len() {
|
||
system_equations.push(EquationRole::SaturatedControl {
|
||
role: "actuator_law",
|
||
});
|
||
system_equations.push(EquationRole::SaturatedControl { role: "tracking" });
|
||
}
|
||
|
||
let n_equations = n_comp_eqs + system_equations.len();
|
||
let n_unknowns = self.full_state_vector_len();
|
||
let balance = SystemDofBalance::from_counts(n_equations, n_unknowns);
|
||
|
||
// Compact unknown catalog (edge layout + extras).
|
||
// Physical block is |B| + 2|E| (+ internal); catalog via edge_to_state.
|
||
let mut branch_ms = std::collections::BTreeSet::new();
|
||
let mut edges: Vec<_> = self.edge_to_state.iter().collect();
|
||
edges.sort_by_key(|(e, _)| e.index());
|
||
for (_edge, (m, _p, _h)) in &edges {
|
||
branch_ms.insert(*m);
|
||
}
|
||
let mut unknowns = Vec::with_capacity(n_unknowns);
|
||
for branch_id in 0..branch_ms.len() {
|
||
unknowns.push(UnknownKind::BranchMassFlow { branch_id });
|
||
}
|
||
for edge_ordinal in 0..edges.len() {
|
||
unknowns.push(UnknownKind::EdgePressure { edge_ordinal });
|
||
unknowns.push(UnknownKind::EdgeEnthalpy { edge_ordinal });
|
||
}
|
||
let physical_accounted = unknowns.len();
|
||
if physical_accounted < self.total_state_len {
|
||
for index in physical_accounted..self.total_state_len {
|
||
unknowns.push(UnknownKind::Internal { index });
|
||
}
|
||
}
|
||
for id in self.inverse_control.linked_controls() {
|
||
unknowns.push(UnknownKind::InverseControl { id: id.to_string() });
|
||
}
|
||
for index in 0..self.coupling_residual_count() {
|
||
unknowns.push(UnknownKind::CouplingHeat { index });
|
||
}
|
||
for index in 0..self.saturated_controllers.len() {
|
||
unknowns.push(UnknownKind::SaturatedActuator { index });
|
||
unknowns.push(UnknownKind::SaturatedIntegrator { index });
|
||
}
|
||
for id in &self.free_actuators {
|
||
unknowns.push(UnknownKind::FreeActuator { id: id.to_string() });
|
||
}
|
||
|
||
if unknowns.len() != n_unknowns {
|
||
diagnostics.push(format!(
|
||
"unknown catalog length {} differs from full_state_vector_len {n_unknowns} \
|
||
(catalog is diagnostic only; counts use full_state_vector_len)",
|
||
unknowns.len()
|
||
));
|
||
}
|
||
|
||
if n_equations != self.total_equation_count() {
|
||
diagnostics.push(format!(
|
||
"internal inconsistency: ledger n_equations {n_equations} != total_equation_count {}",
|
||
self.total_equation_count()
|
||
));
|
||
}
|
||
|
||
DofReport {
|
||
n_equations,
|
||
n_unknowns,
|
||
balance,
|
||
components,
|
||
system_equations,
|
||
unknowns,
|
||
diagnostics,
|
||
}
|
||
}
|
||
|
||
/// Hard DoF gate: returns `Ok(())` only when `n_equations == n_unknowns`.
|
||
///
|
||
/// Call after [`finalize`](Self::finalize). Prefer relying on the automatic
|
||
/// check inside `finalize` for system construction.
|
||
pub fn validate_system_dof(&self) -> Result<(), SystemDofError> {
|
||
let report = self.dof_report();
|
||
if report.balance.is_balanced() {
|
||
Ok(())
|
||
} else {
|
||
Err(SystemDofError::Imbalance {
|
||
n_equations: report.n_equations,
|
||
n_unknowns: report.n_unknowns,
|
||
balance: report.balance,
|
||
summary: report.summary(),
|
||
})
|
||
}
|
||
}
|
||
|
||
/// Returns the state vector index for a control variable.
|
||
///
|
||
/// Control variables are appended after edge states in the state vector:
|
||
///
|
||
/// ```text
|
||
/// State Vector = [Edge States | Control Variables | Thermal Coupling Temps]
|
||
/// [m0, P0, h0, m1, P1, h1, ... | ctrl0, ctrl1, ... | T_hot0, T_cold0, ...]
|
||
/// ```
|
||
///
|
||
/// The index for control variable `i` is: `total_state_len + i`
|
||
///
|
||
/// # Arguments
|
||
///
|
||
/// * `id` - The bounded variable identifier
|
||
///
|
||
/// # Returns
|
||
///
|
||
/// The state vector index, or `None` if the variable is not linked.
|
||
pub fn control_variable_state_index(&self, id: &BoundedVariableId) -> Option<usize> {
|
||
let base = self.total_state_len;
|
||
for (index, linked_id) in self.inverse_control.linked_controls().enumerate() {
|
||
if linked_id == id {
|
||
return Some(base + index);
|
||
}
|
||
}
|
||
for (index, ctrl) in self.saturated_controllers.iter().enumerate() {
|
||
if ctrl.actuator() == id {
|
||
return Some(self.saturated_u_index(index));
|
||
}
|
||
}
|
||
None
|
||
}
|
||
|
||
/// Returns the bounded variable for a given state index.
|
||
pub fn get_bounded_variable_by_state_index(
|
||
&self,
|
||
state_index: usize,
|
||
) -> Option<&BoundedVariable> {
|
||
let base = self.total_state_len;
|
||
if state_index < base {
|
||
return None;
|
||
}
|
||
let control_idx = state_index - base;
|
||
self.inverse_control
|
||
.linked_controls()
|
||
.nth(control_idx)
|
||
.and_then(|id| self.bounded_variables.get(id))
|
||
.or_else(|| {
|
||
let saturated_base = self.saturated_base_index();
|
||
if state_index < saturated_base {
|
||
return None;
|
||
}
|
||
let saturated_offset = state_index - saturated_base;
|
||
if saturated_offset % 2 != 0 {
|
||
return None;
|
||
}
|
||
self.saturated_controllers
|
||
.get(saturated_offset / 2)
|
||
.and_then(|ctrl| self.bounded_variables.get(ctrl.actuator()))
|
||
})
|
||
}
|
||
|
||
/// Returns the bounds (min, max) for a given state index if it corresponds to a bounded control variable.
|
||
pub fn get_bounds_for_state_index(&self, state_index: usize) -> Option<(f64, f64)> {
|
||
self.get_bounded_variable_by_state_index(state_index)
|
||
.map(|var| (var.min(), var.max()))
|
||
}
|
||
|
||
/// Returns physical solver bounds for state variables.
|
||
///
|
||
/// Pressure edge states are strictly positive thermodynamic variables; Newton
|
||
/// steps are clipped before component residuals are evaluated instead of
|
||
/// letting components invent fallback properties for non-physical pressures.
|
||
pub fn get_solver_bounds_for_state_index(&self, state_index: usize) -> Option<(f64, f64)> {
|
||
if let Some(bounds) = self.get_bounds_for_state_index(state_index) {
|
||
return Some(bounds);
|
||
}
|
||
|
||
self.edge_to_state
|
||
.values()
|
||
.any(|(_, p_idx, _)| *p_idx == state_index)
|
||
.then_some((MIN_SOLVER_PRESSURE_PA, f64::MAX))
|
||
}
|
||
|
||
/// Returns the total state vector length including control variables.
|
||
///
|
||
/// ```text
|
||
/// full_state_len = total_state_len
|
||
/// + hard_control_count
|
||
/// + thermal_coupling_count
|
||
/// + 2 * saturated_controller_count
|
||
/// ```
|
||
pub fn full_state_vector_len(&self) -> usize {
|
||
self.total_state_len
|
||
+ self.inverse_control.mapping_count()
|
||
+ self.coupling_residual_count()
|
||
+ 2 * self.saturated_controllers.len()
|
||
+ self.free_actuators.len()
|
||
}
|
||
|
||
/// Base state index of the physical free-actuator block (arch-6).
|
||
/// Free actuators are laid out immediately after the saturated block.
|
||
fn free_actuator_base_index(&self) -> usize {
|
||
self.total_state_len
|
||
+ self.inverse_control.mapping_count()
|
||
+ self.coupling_residual_count()
|
||
+ 2 * self.saturated_controllers.len()
|
||
}
|
||
|
||
/// State index of the i-th physical free actuator (arch-6).
|
||
pub fn free_actuator_index(&self, i: usize) -> usize {
|
||
self.free_actuator_base_index() + i
|
||
}
|
||
|
||
/// State index of the i-th thermal coupling's heat unknown Q [W].
|
||
///
|
||
/// Coupling unknowns live after the hard inverse-control block:
|
||
/// `[Edge States | Hard Controls | Coupling Q | Saturated (u,x) | Free Actuators]`.
|
||
pub fn coupling_state_index(&self, i: usize) -> usize {
|
||
self.total_state_len + self.inverse_control.mapping_count() + i
|
||
}
|
||
|
||
pub fn saturated_base_index(&self) -> usize {
|
||
self.total_state_len + self.inverse_control.mapping_count() + self.coupling_residual_count()
|
||
}
|
||
|
||
pub fn saturated_u_index(&self, i: usize) -> usize {
|
||
self.saturated_base_index() + 2 * i
|
||
}
|
||
|
||
fn saturated_x_index(&self, i: usize) -> usize {
|
||
self.saturated_base_index() + 2 * i + 1
|
||
}
|
||
|
||
/// Returns an ordered list of linked control variable IDs with their state indices.
|
||
///
|
||
/// Useful for constructing the state vector or Jacobian columns.
|
||
pub fn control_variable_indices(&self) -> Vec<(&BoundedVariableId, usize)> {
|
||
let base = self.total_state_len;
|
||
self.inverse_control
|
||
.linked_controls()
|
||
.enumerate()
|
||
.map(|(i, id)| (id, base + i))
|
||
.collect()
|
||
}
|
||
|
||
/// Returns the number of coupling residual equations (one per thermal coupling).
|
||
///
|
||
/// The solver must reserve this many rows in the residual vector for coupling
|
||
/// heat balance equations. Use [`coupling_residuals`](Self::coupling_residuals) to fill them.
|
||
pub fn coupling_residual_count(&self) -> usize {
|
||
self.thermal_couplings.len()
|
||
}
|
||
|
||
/// Fills coupling residuals into `out`.
|
||
///
|
||
/// For each thermal coupling, the residual is the heat transfer rate Q (W) into the cold
|
||
/// circuit: Q = η·UA·(T_hot − T_cold). The solver typically uses this in a heat balance
|
||
/// (e.g. r = Q_actual − Q_expected). Temperatures must be in Kelvin.
|
||
///
|
||
/// # Arguments
|
||
///
|
||
/// * `temperatures` - One (T_hot_K, T_cold_K) per coupling; length must equal
|
||
/// `thermal_coupling_count()`. The solver obtains these from state (e.g. P, h → T via fluid backend).
|
||
/// * `out` - Slice to write residuals; length must be at least `coupling_residual_count()`.
|
||
pub fn coupling_residuals(&self, temperatures: &[(f64, f64)], out: &mut [f64]) {
|
||
assert!(
|
||
temperatures.len() == self.thermal_couplings.len(),
|
||
"temperatures.len() must equal thermal_coupling_count()"
|
||
);
|
||
assert!(
|
||
out.len() >= self.thermal_couplings.len(),
|
||
"out.len() must be at least coupling_residual_count()"
|
||
);
|
||
for (i, coupling) in self.thermal_couplings.iter().enumerate() {
|
||
let (t_hot_k, t_cold_k) = temperatures[i];
|
||
let t_hot = Temperature::from_kelvin(t_hot_k);
|
||
let t_cold = Temperature::from_kelvin(t_cold_k);
|
||
out[i] = crate::coupling::compute_coupling_heat(coupling, t_hot, t_cold);
|
||
}
|
||
}
|
||
|
||
/// Returns Jacobian entries for coupling residuals with respect to temperature state.
|
||
///
|
||
/// The solver state may include temperature unknowns for coupling interfaces (or T derived from P, h).
|
||
/// For each coupling i: ∂Q_i/∂T_hot = η·UA, ∂Q_i/∂T_cold = −η·UA.
|
||
///
|
||
/// # Arguments
|
||
///
|
||
/// * `row_offset` - First row index for coupling equations in the global residual vector.
|
||
/// * `t_hot_cols` - State column index for T_hot per coupling; length = `thermal_coupling_count()`.
|
||
/// * `t_cold_cols` - State column index for T_cold per coupling; length = `thermal_coupling_count()`.
|
||
///
|
||
/// # Returns
|
||
///
|
||
/// `(row, col, value)` tuples for the Jacobian. Row is `row_offset + coupling_index`.
|
||
pub fn coupling_jacobian_entries(
|
||
&self,
|
||
row_offset: usize,
|
||
t_hot_cols: &[usize],
|
||
t_cold_cols: &[usize],
|
||
) -> Vec<(usize, usize, f64)> {
|
||
assert!(
|
||
t_hot_cols.len() == self.thermal_couplings.len()
|
||
&& t_cold_cols.len() == self.thermal_couplings.len(),
|
||
"t_hot_cols and t_cold_cols length must equal thermal_coupling_count()"
|
||
);
|
||
let mut entries = Vec::with_capacity(2 * self.thermal_couplings.len());
|
||
for (i, coupling) in self.thermal_couplings.iter().enumerate() {
|
||
let dr_dt_hot = coupling.efficiency * coupling.ua.to_watts_per_kelvin();
|
||
let dr_dt_cold = -dr_dt_hot;
|
||
let row = row_offset + i;
|
||
entries.push((row, t_hot_cols[i], dr_dt_hot));
|
||
entries.push((row, t_cold_cols[i], dr_dt_cold));
|
||
}
|
||
entries
|
||
}
|
||
|
||
/// Returns true if the graph contains a cycle.
|
||
///
|
||
/// Refrigeration circuits form cycles (compressor → condenser → valve → evaporator → compressor),
|
||
/// so cycles are expected and valid.
|
||
pub fn is_cyclic(&self) -> bool {
|
||
algo::is_cyclic_directed(&self.graph)
|
||
}
|
||
|
||
/// Iterates over (component, edge_indices) for Jacobian assembly.
|
||
///
|
||
/// For each node, yields the component and a map from edge index to (state_index_p, state_index_h)
|
||
/// for edges incident to that node (incoming and outgoing).
|
||
///
|
||
/// # Panics
|
||
///
|
||
/// Panics if `finalize()` has not been called.
|
||
pub fn traverse_for_jacobian(
|
||
&self,
|
||
) -> impl Iterator<Item = (NodeIndex, &dyn Component, Vec<(EdgeIndex, usize, usize)>)> {
|
||
assert!(
|
||
self.finalized,
|
||
"call finalize() before traverse_for_jacobian()"
|
||
);
|
||
|
||
self.graph.node_indices().map(move |node_idx| {
|
||
let component = self.graph.node_weight(node_idx).unwrap();
|
||
let mut edge_indices = Vec::new();
|
||
|
||
for edge_ref in self
|
||
.graph
|
||
.edges_directed(node_idx, petgraph::Direction::Incoming)
|
||
{
|
||
let edge_idx = edge_ref.id();
|
||
if let Some(&(_m, p, h)) = self.edge_to_state.get(&edge_idx) {
|
||
edge_indices.push((edge_idx, p, h));
|
||
}
|
||
}
|
||
for edge_ref in self
|
||
.graph
|
||
.edges_directed(node_idx, petgraph::Direction::Outgoing)
|
||
{
|
||
let edge_idx = edge_ref.id();
|
||
if let Some(&(_m, p, h)) = self.edge_to_state.get(&edge_idx) {
|
||
edge_indices.push((edge_idx, p, h));
|
||
}
|
||
}
|
||
|
||
(node_idx, component.as_ref(), edge_indices)
|
||
})
|
||
}
|
||
|
||
/// Assembles residuals from all components.
|
||
///
|
||
/// Components receive the full state slice and write to their equation indices.
|
||
/// Equation indices are computed from component order and `n_equations()`.
|
||
///
|
||
/// # Errors
|
||
///
|
||
/// Returns `ComponentError::InvalidResidualDimensions` if `residuals.len()` is
|
||
/// less than the total number of equations across all components.
|
||
pub fn compute_residuals(
|
||
&self,
|
||
state: &StateSlice,
|
||
residuals: &mut ResidualVector,
|
||
) -> Result<(), ComponentError> {
|
||
let mut total_eqs: usize = self
|
||
.traverse_for_jacobian()
|
||
.map(|(_, c, _)| c.n_equations())
|
||
.sum();
|
||
total_eqs += self.constraints.len()
|
||
+ self.coupling_residual_count()
|
||
+ 2 * self.saturated_controllers.len();
|
||
|
||
if residuals.len() < total_eqs {
|
||
return Err(ComponentError::InvalidResidualDimensions {
|
||
expected: total_eqs,
|
||
actual: residuals.len(),
|
||
});
|
||
}
|
||
|
||
let mut eq_offset = 0;
|
||
for (_node_idx, component, _edge_indices) in self.traverse_for_jacobian() {
|
||
let n = component.n_equations();
|
||
if n > 0 {
|
||
let mut temp = vec![0.0; n];
|
||
component.compute_residuals(state, &mut temp)?;
|
||
residuals[eq_offset..eq_offset + n].copy_from_slice(&temp);
|
||
}
|
||
eq_offset += n;
|
||
}
|
||
|
||
// Add constraints
|
||
let control_values: Vec<f64> = self
|
||
.control_variable_indices()
|
||
.into_iter()
|
||
.map(|(_, idx)| state[idx])
|
||
.collect();
|
||
let measured = self.extract_constraint_values_with_controls(state, &control_values);
|
||
let n_constraints = self
|
||
.compute_constraint_residuals(state, &mut residuals[eq_offset..], &measured)
|
||
.map_err(|e| ComponentError::CalculationFailed(e.to_string()))?;
|
||
eq_offset += n_constraints;
|
||
|
||
// Add couplings. Each coupling owns one unknown Q = state[coupling_state_index(i)].
|
||
//
|
||
// * Physical mode (`hot_component`/`cold_component` set): the residual
|
||
// closes Q against the hot component's *measured* duty:
|
||
// r_i = Q − η·duty_hot(state). The cold-side `ThermalLoad` consumes Q
|
||
// in its own energy balance (wired in `finalize()`), so the heat
|
||
// genuinely crosses the circuit boundary and the First Law closes.
|
||
// * Legacy mode (no components referenced): the residual pins the orphan
|
||
// unknown, r_i = Q (⇒ Q = 0). This replaces the former 300 K stub,
|
||
// whose all-zero row/column made the Jacobian structurally singular.
|
||
let n_couplings = self.coupling_residual_count();
|
||
for (i, coupling) in self.thermal_couplings.iter().enumerate() {
|
||
let q_idx = self.coupling_state_index(i);
|
||
let q = state.get(q_idx).copied().unwrap_or(0.0);
|
||
let target = if coupling.is_physical() {
|
||
let duty = coupling
|
||
.hot_component
|
||
.as_deref()
|
||
.and_then(|name| self.component_names.get(name))
|
||
.and_then(|&node| self.graph.node_weight(node))
|
||
.and_then(|c| {
|
||
c.measure_output(entropyk_components::MeasuredOutput::Capacity, state)
|
||
})
|
||
.unwrap_or(0.0);
|
||
coupling.efficiency * coupling.duty_scale * duty
|
||
} else {
|
||
0.0
|
||
};
|
||
residuals[eq_offset + i] = q - target;
|
||
}
|
||
eq_offset += n_couplings;
|
||
|
||
for (i, ctrl) in self.saturated_controllers.iter().enumerate() {
|
||
let u_idx = self.saturated_u_index(i);
|
||
let x_idx = self.saturated_x_index(i);
|
||
if state.len() <= x_idx {
|
||
return Err(ComponentError::InvalidStateDimensions {
|
||
expected: x_idx + 1,
|
||
actual: state.len(),
|
||
});
|
||
}
|
||
let u = state[u_idx];
|
||
let x = state[x_idx];
|
||
|
||
residuals[eq_offset] = ctrl.residual_u(u, x);
|
||
eq_offset += 1;
|
||
if ctrl.is_network() {
|
||
let measured: Vec<f64> = ctrl
|
||
.objectives()
|
||
.iter()
|
||
.map(|o| self.measure_saturated_output(&o.output, state))
|
||
.collect();
|
||
let e = ctrl.network_error(&measured);
|
||
residuals[eq_offset] = ctrl.residual_y_network(e, x);
|
||
} else {
|
||
let y = self.measure_saturated_output(ctrl.output(), state);
|
||
residuals[eq_offset] = ctrl.residual_y(y, x);
|
||
}
|
||
eq_offset += 1;
|
||
}
|
||
|
||
Ok(())
|
||
}
|
||
|
||
/// Assembles Jacobian entries from all components.
|
||
///
|
||
/// Each component receives the state and writes to JacobianBuilder. The
|
||
/// [`traverse_for_jacobian`](Self::traverse_for_jacobian) iterator provides
|
||
/// the edge→state mapping `(EdgeIndex, state_index_p, state_index_h)` per
|
||
/// component. Components must know their port→state mapping (e.g. from graph
|
||
/// construction in Story 3.2) to write correct column indices.
|
||
pub fn assemble_jacobian(
|
||
&self,
|
||
state: &StateSlice,
|
||
jacobian: &mut JacobianBuilder,
|
||
) -> Result<(), ComponentError> {
|
||
let mut row_offset = 0;
|
||
for (_node_idx, component, _edge_indices) in self.traverse_for_jacobian() {
|
||
let n = component.n_equations();
|
||
if n > 0 {
|
||
// Components write rows 0..n-1; we offset to global equation indices.
|
||
let mut temp_builder = JacobianBuilder::new();
|
||
component.jacobian_entries(state, &mut temp_builder)?;
|
||
for (r, c, v) in temp_builder.entries() {
|
||
jacobian.add_entry(row_offset + r, *c, *v);
|
||
}
|
||
}
|
||
row_offset += n;
|
||
}
|
||
|
||
// Add constraints jacobian
|
||
let control_values: Vec<f64> = self
|
||
.control_variable_indices()
|
||
.into_iter()
|
||
.map(|(_, idx)| state[idx])
|
||
.collect();
|
||
let constraint_jac =
|
||
self.compute_inverse_control_jacobian(state, row_offset, &control_values);
|
||
for (r, c, v) in constraint_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();
|
||
if !self.thermal_couplings.is_empty() {
|
||
let eps = self.inverse_control.finite_diff_epsilon();
|
||
let mut state_mut = state.to_vec();
|
||
for (i, coupling) in self.thermal_couplings.iter().enumerate() {
|
||
let row = coupling_row_offset + i;
|
||
jacobian.add_entry(row, self.coupling_state_index(i), 1.0);
|
||
if !coupling.is_physical() {
|
||
continue;
|
||
}
|
||
let Some(&node_idx) = coupling
|
||
.hot_component
|
||
.as_deref()
|
||
.and_then(|name| self.component_names.get(name))
|
||
else {
|
||
continue;
|
||
};
|
||
let measure = |st: &[f64]| -> f64 {
|
||
self.graph
|
||
.node_weight(node_idx)
|
||
.and_then(|c| {
|
||
c.measure_output(entropyk_components::MeasuredOutput::Capacity, st)
|
||
})
|
||
.unwrap_or(0.0)
|
||
};
|
||
for col in self.incident_state_indices_for_component(node_idx) {
|
||
let orig = state_mut[col];
|
||
state_mut[col] = orig + eps;
|
||
let d_plus = measure(&state_mut);
|
||
state_mut[col] = orig - eps;
|
||
let d_minus = measure(&state_mut);
|
||
state_mut[col] = orig;
|
||
let dduty_dcol = (d_plus - d_minus) / (2.0 * eps);
|
||
let derivative = -coupling.efficiency * coupling.duty_scale * dduty_dcol;
|
||
if derivative.abs() > 1e-10 {
|
||
jacobian.add_entry(row, col, derivative);
|
||
}
|
||
}
|
||
}
|
||
}
|
||
|
||
let saturated_row_offset =
|
||
row_offset + self.constraints.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);
|
||
}
|
||
Ok(())
|
||
}
|
||
|
||
/// Tolerance for mass balance validation [kg/s].
|
||
///
|
||
/// This value (1e-9 kg/s) is tight enough to catch numerical issues
|
||
/// while allowing for floating-point rounding errors.
|
||
pub const MASS_BALANCE_TOLERANCE_KG_S: f64 = 1e-9;
|
||
|
||
/// Tolerance for energy balance validation in Watts (1e-6 kW)
|
||
pub const ENERGY_BALANCE_TOLERANCE_W: f64 = 1e-3;
|
||
|
||
/// Verifies that global mass balance is conserved.
|
||
///
|
||
/// Sums the mass flow rates at the ports of each component and ensures they
|
||
/// sum to zero within a tight tolerance (1e-9 kg/s).
|
||
///
|
||
/// # Returns
|
||
///
|
||
/// * `Ok(())` if all components pass mass balance validation
|
||
/// * `Err(SolverError::Validation)` if any component violates mass conservation
|
||
///
|
||
/// # Note
|
||
///
|
||
/// Components without `port_mass_flows` implementation are logged as warnings
|
||
/// and skipped. This ensures visibility of incomplete implementations without
|
||
/// failing the validation.
|
||
pub fn check_mass_balance(&self, state: &StateSlice) -> Result<(), crate::SolverError> {
|
||
let mut total_mass_error = 0.0;
|
||
let mut has_violation = false;
|
||
let mut components_checked = 0usize;
|
||
let mut components_skipped = 0usize;
|
||
|
||
for (node_idx, component, _edge_indices) in self.traverse_for_jacobian() {
|
||
match component.port_mass_flows(state) {
|
||
Ok(mass_flows) => {
|
||
let sum: f64 = mass_flows.iter().map(|m| m.to_kg_per_s()).sum();
|
||
if sum.abs() > Self::MASS_BALANCE_TOLERANCE_KG_S {
|
||
has_violation = true;
|
||
total_mass_error += sum.abs();
|
||
tracing::warn!(
|
||
node_index = node_idx.index(),
|
||
mass_imbalance_kg_s = sum,
|
||
"Mass balance violation detected at component"
|
||
);
|
||
}
|
||
components_checked += 1;
|
||
}
|
||
Err(e) => {
|
||
components_skipped += 1;
|
||
tracing::warn!(
|
||
node_index = node_idx.index(),
|
||
error = %e,
|
||
"Component does not implement port_mass_flows - skipping mass balance check"
|
||
);
|
||
}
|
||
}
|
||
}
|
||
|
||
tracing::debug!(
|
||
components_checked,
|
||
components_skipped,
|
||
total_mass_error_kg_s = total_mass_error,
|
||
"Mass balance validation complete"
|
||
);
|
||
|
||
if has_violation {
|
||
return Err(crate::SolverError::Validation {
|
||
mass_error: total_mass_error,
|
||
energy_error: 0.0,
|
||
});
|
||
}
|
||
Ok(())
|
||
}
|
||
|
||
/// Verifies the First Law of Thermodynamics for all components in the system.
|
||
///
|
||
/// Validates that ΣQ - ΣW + Σ(ṁ·h) = 0 for each component.
|
||
/// Returns `SolverError::Validation` if any component violates the balance.
|
||
pub fn check_energy_balance(&self, state: &StateSlice) -> Result<(), crate::SolverError> {
|
||
let mut total_energy_error = 0.0;
|
||
let mut has_violation = false;
|
||
let mut components_checked = 0usize;
|
||
let mut components_skipped = 0usize;
|
||
let mut skipped_components: Vec<String> = Vec::new();
|
||
|
||
for (node_idx, component, _edge_indices) in self.traverse_for_jacobian() {
|
||
let energy_transfers = component.energy_transfers(state);
|
||
let mass_flows = component.port_mass_flows(state);
|
||
let enthalpies = component.port_enthalpies(state);
|
||
|
||
match (energy_transfers, mass_flows, enthalpies) {
|
||
(Some((heat, work)), Ok(m_flows), Ok(h_flows))
|
||
if m_flows.len() == h_flows.len() =>
|
||
{
|
||
let mut net_energy_flow = 0.0;
|
||
for (m, h) in m_flows.iter().zip(h_flows.iter()) {
|
||
net_energy_flow += m.to_kg_per_s() * h.to_joules_per_kg();
|
||
}
|
||
|
||
let balance = heat.to_watts() - work.to_watts() + net_energy_flow;
|
||
|
||
if balance.abs() > Self::ENERGY_BALANCE_TOLERANCE_W {
|
||
has_violation = true;
|
||
total_energy_error += balance.abs();
|
||
tracing::warn!(
|
||
node_index = node_idx.index(),
|
||
energy_imbalance_w = balance,
|
||
"Energy balance violation detected at component"
|
||
);
|
||
}
|
||
components_checked += 1;
|
||
}
|
||
_ => {
|
||
components_skipped += 1;
|
||
let component_type = std::any::type_name_of_val(component)
|
||
.split("::")
|
||
.last()
|
||
.unwrap_or("unknown");
|
||
let component_info =
|
||
format!("{} (type: {})", component.signature(), component_type);
|
||
skipped_components.push(component_info.clone());
|
||
|
||
tracing::warn!(
|
||
component = %component_info,
|
||
node_index = node_idx.index(),
|
||
"Component lacks energy_transfers() or port_enthalpies() - SKIPPED in energy balance validation"
|
||
);
|
||
}
|
||
}
|
||
}
|
||
|
||
// Summary warning if components were skipped
|
||
if components_skipped > 0 {
|
||
tracing::warn!(
|
||
components_checked = components_checked,
|
||
components_skipped = components_skipped,
|
||
skipped = ?skipped_components,
|
||
"Energy balance validation incomplete: {} component(s) skipped. \
|
||
Implement energy_transfers() and port_enthalpies() for full validation.",
|
||
components_skipped
|
||
);
|
||
} else {
|
||
tracing::debug!(
|
||
components_checked,
|
||
components_skipped,
|
||
total_energy_error_w = total_energy_error,
|
||
"Energy balance validation complete"
|
||
);
|
||
}
|
||
|
||
if has_violation {
|
||
return Err(crate::SolverError::Validation {
|
||
mass_error: 0.0,
|
||
energy_error: total_energy_error,
|
||
});
|
||
}
|
||
Ok(())
|
||
}
|
||
|
||
/// Computes the aggregate cycle performance (cooling/heating duty, work input,
|
||
/// COP/EER) from the solved state vector.
|
||
///
|
||
/// Each component's [`Component::energy_transfers`] is evaluated at `state`
|
||
/// and classified by the physical sign convention (`Q > 0` = heat into the
|
||
/// refrigerant = cooling; `Q < 0` = heat rejected = heating; `W < 0` = work
|
||
/// done on the fluid = power input). Components that do not report energy
|
||
/// transfers are skipped. Returns `None` if no component reported anything.
|
||
///
|
||
/// This is post-processing on a converged state — it makes no solver calls
|
||
/// and does not mutate the system.
|
||
pub fn cycle_performance(&self, state: &StateSlice) -> Option<CyclePerformance> {
|
||
let mut q_cooling_w = 0.0;
|
||
let mut q_heating_w = 0.0;
|
||
let mut work_input_w = 0.0;
|
||
let mut components_counted = 0usize;
|
||
|
||
for (_node_idx, component, _edges) in self.traverse_for_jacobian() {
|
||
// Skip secondary-side receivers (e.g. `ThermalLoad`): the heat they
|
||
// absorb is the primary cycle's rejected duty, not extra cooling.
|
||
if !component.counts_in_cycle_performance() {
|
||
continue;
|
||
}
|
||
if let Some((heat, work)) = component.energy_transfers(state) {
|
||
let q = heat.to_watts();
|
||
let w = work.to_watts();
|
||
if !(q.is_finite() && w.is_finite()) {
|
||
continue;
|
||
}
|
||
if q > 0.0 {
|
||
q_cooling_w += q;
|
||
} else if q < 0.0 {
|
||
q_heating_w += -q;
|
||
}
|
||
// Work done ON the fluid is reported negative; accumulate its
|
||
// magnitude as the electrical/shaft power input.
|
||
if w < 0.0 {
|
||
work_input_w += -w;
|
||
}
|
||
components_counted += 1;
|
||
}
|
||
}
|
||
|
||
if components_counted == 0 {
|
||
return None;
|
||
}
|
||
Some(CyclePerformance {
|
||
q_cooling_w,
|
||
q_heating_w,
|
||
work_input_w,
|
||
components_counted,
|
||
})
|
||
}
|
||
|
||
/// Generates a deterministic byte representation of the system configuration.
|
||
/// Used for simulation traceability logic.
|
||
pub fn generate_canonical_bytes(&self) -> Vec<u8> {
|
||
let mut repr = String::new();
|
||
repr.push_str("Nodes:\n");
|
||
// To be deterministic, we just iterate in graph order which is stable
|
||
// as long as we don't delete nodes.
|
||
for node in self.graph.node_indices() {
|
||
let circuit_id = self.node_to_circuit.get(&node).map(|c| c.0).unwrap_or(0);
|
||
repr.push_str(&format!(
|
||
" Node({}): Circuit({})\n",
|
||
node.index(),
|
||
circuit_id
|
||
));
|
||
if let Some(comp) = self.graph.node_weight(node) {
|
||
repr.push_str(&format!(" Signature: {}\n", comp.signature()));
|
||
}
|
||
}
|
||
repr.push_str("Edges:\n");
|
||
for edge_idx in self.graph.edge_indices() {
|
||
if let Some((src, tgt)) = self.graph.edge_endpoints(edge_idx) {
|
||
repr.push_str(&format!(" Edge: {} -> {}\n", src.index(), tgt.index()));
|
||
}
|
||
}
|
||
repr.push_str("Thermal Couplings:\n");
|
||
let mut couplings: Vec<_> = self.thermal_couplings.iter().collect();
|
||
couplings.sort_by(|a, b| {
|
||
(a.hot_circuit.0, a.cold_circuit.0)
|
||
.cmp(&(b.hot_circuit.0, b.cold_circuit.0))
|
||
.then(
|
||
a.ua.0
|
||
.partial_cmp(&b.ua.0)
|
||
.unwrap_or(std::cmp::Ordering::Equal),
|
||
)
|
||
.then(
|
||
a.efficiency
|
||
.partial_cmp(&b.efficiency)
|
||
.unwrap_or(std::cmp::Ordering::Equal),
|
||
)
|
||
});
|
||
for coupling in couplings {
|
||
repr.push_str(&format!(
|
||
" Hot: {}, Cold: {}, UA: {}\n",
|
||
coupling.hot_circuit.0, coupling.cold_circuit.0, coupling.ua
|
||
));
|
||
}
|
||
repr.push_str("Constraints:\n");
|
||
let mut constraint_keys: Vec<_> = self.constraints.keys().collect();
|
||
constraint_keys.sort_by_key(|k| k.as_str());
|
||
for key in constraint_keys {
|
||
let c = &self.constraints[key];
|
||
repr.push_str(&format!(" {}: {}\n", c.id().as_str(), c.target_value()));
|
||
}
|
||
repr.push_str("Bounded Variables:\n");
|
||
let mut bounded_keys: Vec<_> = self.bounded_variables.keys().collect();
|
||
bounded_keys.sort_by_key(|k| k.as_str());
|
||
for key in bounded_keys {
|
||
let var = &self.bounded_variables[key];
|
||
repr.push_str(&format!(
|
||
" {}: [{}, {}]\n",
|
||
var.id().as_str(),
|
||
var.min(),
|
||
var.max()
|
||
));
|
||
}
|
||
|
||
repr.push_str("Inverse Control Mappings:\n");
|
||
// For inverse control mappings, they are ordered internally. We'll just iterate linked controls.
|
||
for (i, (constraint, bounded_var)) in self.inverse_control.mappings().enumerate() {
|
||
repr.push_str(&format!(
|
||
" Mapping {}: {} -> {}\n",
|
||
i,
|
||
constraint.as_str(),
|
||
bounded_var.as_str()
|
||
));
|
||
}
|
||
|
||
repr.into_bytes()
|
||
}
|
||
|
||
/// Computes the SHA-256 hash uniquely identifying the input configuration.
|
||
pub fn input_hash(&self) -> String {
|
||
use sha2::{Digest, Sha256};
|
||
let mut hasher = Sha256::new();
|
||
hasher.update(self.generate_canonical_bytes());
|
||
format!("{:064x}", hasher.finalize())
|
||
}
|
||
|
||
// ========== JSON Serialization API ==========
|
||
|
||
/// Serializes the system to a JSON string.
|
||
///
|
||
/// This method captures the complete system state including topology,
|
||
/// component parameters, and metadata in a human-readable JSON format.
|
||
///
|
||
/// # Errors
|
||
///
|
||
/// Returns `ThermoError::SerializationError` if JSON serialization fails.
|
||
///
|
||
/// # Examples
|
||
///
|
||
/// ```no_run
|
||
/// # use entropyk_solver::System;
|
||
/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
|
||
/// let system = System::new();
|
||
/// let json_string = system.to_json_string()?;
|
||
/// println!("System JSON: {}", json_string);
|
||
/// # Ok(())
|
||
/// # }
|
||
/// ```
|
||
pub fn to_json_string(&self) -> Result<String, crate::error::ThermoError> {
|
||
use crate::snapshot::{
|
||
BoundedVariableSnapshot, ConstraintSnapshot, EdgeSnapshot, FluidBackendInfo,
|
||
SolverConfigSnapshot, SystemSnapshot, TopologySnapshot,
|
||
};
|
||
use std::collections::HashMap;
|
||
|
||
tracing::info!("Serializing system to JSON");
|
||
|
||
let reverse_names: HashMap<NodeIndex, &String> =
|
||
self.component_names.iter().map(|(n, &i)| (i, n)).collect();
|
||
|
||
// Extract topology with port names
|
||
let mut edges = Vec::new();
|
||
for edge in self.graph.edge_indices() {
|
||
let (source, target) = self.graph.edge_endpoints(edge).unwrap();
|
||
let source_node = self.graph.node_weight(source).unwrap();
|
||
let target_node = self.graph.node_weight(target).unwrap();
|
||
|
||
// Derive port names from component port_names() or defaults
|
||
let source_ports = source_node.port_names();
|
||
let target_ports = target_node.port_names();
|
||
|
||
// Count how many edges connect TO the target (this edge's index at target)
|
||
let target_incoming: Vec<_> = self
|
||
.graph
|
||
.edges_directed(target, petgraph::Direction::Incoming)
|
||
.collect();
|
||
let target_port_idx = target_incoming
|
||
.iter()
|
||
.position(|e| e.id() == edge)
|
||
.unwrap_or(0);
|
||
|
||
// Count how many edges leave FROM the source (this edge's index at source)
|
||
let source_outgoing: Vec<_> = self
|
||
.graph
|
||
.edges_directed(source, petgraph::Direction::Outgoing)
|
||
.collect();
|
||
let source_port_idx = source_outgoing
|
||
.iter()
|
||
.position(|e| e.id() == edge)
|
||
.unwrap_or(0);
|
||
|
||
let source_port_name = source_ports
|
||
.get(source_port_idx)
|
||
.cloned()
|
||
.unwrap_or_else(|| format!("port_{}", source_port_idx));
|
||
let target_port_name = target_ports
|
||
.get(target_port_idx)
|
||
.cloned()
|
||
.unwrap_or_else(|| format!("port_{}", target_port_idx));
|
||
|
||
edges.push(EdgeSnapshot {
|
||
source: reverse_names
|
||
.get(&source)
|
||
.map(|s| s.to_string())
|
||
.unwrap_or_else(|| source_node.signature()),
|
||
source_port: source_port_name,
|
||
target: reverse_names
|
||
.get(&target)
|
||
.map(|s| s.to_string())
|
||
.unwrap_or_else(|| target_node.signature()),
|
||
target_port: target_port_name,
|
||
circuit_id: self.edge_circuit(edge).0,
|
||
});
|
||
}
|
||
|
||
// Extract component parameters (use unique key: registered name or signature+index)
|
||
let mut parameters = HashMap::new();
|
||
for node in self.graph.node_indices() {
|
||
if let Some(component) = self.graph.node_weight(node) {
|
||
let params = component.to_params();
|
||
let key = reverse_names
|
||
.get(&node)
|
||
.map(|s| (*s).clone())
|
||
.unwrap_or_else(|| component.signature());
|
||
parameters.insert(key.to_string(), params);
|
||
}
|
||
}
|
||
|
||
// Build component_names and circuit_assignments maps
|
||
let component_names: HashMap<String, String> = self
|
||
.component_names
|
||
.iter()
|
||
.map(|(name, &node_idx)| {
|
||
let comp = self.graph.node_weight(node_idx);
|
||
let type_name = comp
|
||
.map(|c| {
|
||
let sig = c.to_params().component_type.clone();
|
||
sig
|
||
})
|
||
.unwrap_or_else(|| "Unknown".to_string());
|
||
(name.clone(), type_name)
|
||
})
|
||
.collect();
|
||
|
||
let circuit_assignments: HashMap<String, u16> = self
|
||
.component_names
|
||
.iter()
|
||
.map(|(name, &node_idx)| {
|
||
let cid = self
|
||
.node_to_circuit
|
||
.get(&node_idx)
|
||
.map(|c| c.0)
|
||
.unwrap_or(0);
|
||
(name.clone(), cid)
|
||
})
|
||
.collect();
|
||
|
||
// Create snapshot
|
||
let snapshot = SystemSnapshot {
|
||
version: "1.0".to_string(),
|
||
topology: TopologySnapshot {
|
||
edges,
|
||
thermal_couplings: self.thermal_couplings.clone(),
|
||
},
|
||
parameters,
|
||
fluid_state: {
|
||
let mut data = Vec::with_capacity(self.graph.edge_count() * 2);
|
||
for edge in self.graph.edge_indices() {
|
||
let (source, _target) = self.graph.edge_endpoints(edge).unwrap();
|
||
let component = self.graph.node_weight(source).unwrap();
|
||
let ports = component.get_ports();
|
||
let outgoing: Vec<_> = self
|
||
.graph
|
||
.edges_directed(source, petgraph::Direction::Outgoing)
|
||
.collect();
|
||
let port_idx = outgoing.iter().position(|e| e.id() == edge).unwrap_or(0);
|
||
if let Some(port) = ports.get(port_idx) {
|
||
data.push(port.pressure().to_pascals());
|
||
data.push(port.enthalpy().to_joules_per_kg());
|
||
} else {
|
||
data.push(0.0);
|
||
data.push(0.0);
|
||
}
|
||
}
|
||
if data.is_empty() {
|
||
None
|
||
} else {
|
||
entropyk_core::SystemState::try_from(data).ok()
|
||
}
|
||
},
|
||
fluid_backend: FluidBackendInfo {
|
||
name: "CoolPropBackend".to_string(),
|
||
version: env!("CARGO_PKG_VERSION").to_string(),
|
||
hash: None,
|
||
},
|
||
solver_config: Some(SolverConfigSnapshot::default()),
|
||
component_names,
|
||
circuit_assignments,
|
||
constraints: self
|
||
.constraints
|
||
.iter()
|
||
.map(|(id, c)| ConstraintSnapshot {
|
||
id: id.as_str().to_string(),
|
||
component: c.output().component_id().to_string(),
|
||
output_type: c.output().constraint_type_name().to_string(),
|
||
target: c.target_value(),
|
||
})
|
||
.collect(),
|
||
bounded_variables: self
|
||
.bounded_variables
|
||
.iter()
|
||
.map(|(id, v)| BoundedVariableSnapshot {
|
||
id: id.as_str().to_string(),
|
||
component: v.component_id().unwrap_or("").to_string(),
|
||
variable_name: id.as_str().to_string(),
|
||
lower_bound: v.min(),
|
||
upper_bound: v.max(),
|
||
initial_value: v.initial_value(),
|
||
})
|
||
.collect(),
|
||
metadata: HashMap::new(),
|
||
};
|
||
|
||
// Serialize to JSON with pretty printing
|
||
serde_json::to_string_pretty(&snapshot).map_err(|e| {
|
||
crate::error::ThermoError::SerializationError(format!(
|
||
"JSON serialization failed: {}",
|
||
e
|
||
))
|
||
})
|
||
}
|
||
|
||
/// Deserializes a system from a JSON string.
|
||
///
|
||
/// Reconstructs a system from a previously serialized JSON representation.
|
||
/// Validates version compatibility and backend requirements.
|
||
///
|
||
/// # Errors
|
||
///
|
||
/// - `ThermoError::DeserializationError` if JSON parsing fails
|
||
/// - `ThermoError::VersionMismatch` if the schema version is incompatible
|
||
/// - `ThermoError::BackendUnavailable` if the required fluid backend is not available
|
||
///
|
||
/// # Examples
|
||
///
|
||
/// ```no_run
|
||
/// # use entropyk_solver::System;
|
||
/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
|
||
/// let json_string = r#"{"version": "1.0", ...}"#;
|
||
/// let system = System::from_json_string(json_string)?;
|
||
/// # Ok(())
|
||
/// # }
|
||
/// ```
|
||
pub fn from_json_string(json_str: &str) -> Result<Self, crate::error::ThermoError> {
|
||
use crate::snapshot::SystemSnapshot;
|
||
|
||
tracing::info!("Deserializing system from JSON");
|
||
|
||
// Parse JSON
|
||
let snapshot: SystemSnapshot = serde_json::from_str(json_str).map_err(|e| {
|
||
crate::error::ThermoError::DeserializationError(format!("JSON parsing failed: {}", e))
|
||
})?;
|
||
|
||
// Validate version
|
||
if snapshot.version != "1.0" {
|
||
return Err(crate::error::ThermoError::VersionMismatch {
|
||
expected: "1.0".to_string(),
|
||
found: snapshot.version,
|
||
});
|
||
}
|
||
|
||
// Log backend info
|
||
tracing::debug!(
|
||
"Fluid backend: {} v{}",
|
||
snapshot.fluid_backend.name,
|
||
snapshot.fluid_backend.version
|
||
);
|
||
|
||
// Validate backend availability (AC5: explicit error for missing backend)
|
||
let backend_name = &snapshot.fluid_backend.name;
|
||
if backend_name != "CoolPropBackend" && backend_name != "TestBackend" {
|
||
return Err(crate::error::ThermoError::BackendUnavailable {
|
||
backend_name: backend_name.clone(),
|
||
required_version: snapshot.fluid_backend.version,
|
||
});
|
||
}
|
||
|
||
// Build name → parameter lookup for ordering
|
||
let mut system = System::new();
|
||
|
||
// Track component names → NodeIndex for edge reconstruction
|
||
let mut name_to_node: HashMap<String, NodeIndex> = HashMap::new();
|
||
|
||
// Reconstruct components from parameters
|
||
// We iterate in a deterministic order: sorted by key name
|
||
let mut sorted_keys: Vec<&String> = snapshot.parameters.keys().collect();
|
||
sorted_keys.sort();
|
||
|
||
for key in sorted_keys {
|
||
let params = &snapshot.parameters[key];
|
||
let type_name = params.component_type.as_str();
|
||
|
||
// Use registry for supported types
|
||
let component: Box<dyn Component> = match entropyk_components::create_component(params)
|
||
{
|
||
Ok(c) => c,
|
||
Err(_) => {
|
||
// For unsupported types, create a minimal placeholder
|
||
// that preserves the topology and parameters
|
||
tracing::warn!(
|
||
"Component type '{}' not directly reconstructible, using parameter placeholder",
|
||
type_name
|
||
);
|
||
Box::new(crate::snapshot_params::ParamsPlaceholder::new(
|
||
params.clone(),
|
||
))
|
||
}
|
||
};
|
||
|
||
// Get circuit ID from snapshot
|
||
let circuit_id = snapshot
|
||
.circuit_assignments
|
||
.get(key)
|
||
.map(|&id| CircuitId(id))
|
||
.unwrap_or(CircuitId::ZERO);
|
||
|
||
let node = system
|
||
.add_component_to_circuit(component, circuit_id)
|
||
.map_err(|e| {
|
||
crate::error::ThermoError::DeserializationError(format!(
|
||
"Failed to add component '{}': {:?}",
|
||
key, e
|
||
))
|
||
})?;
|
||
|
||
system.register_component_name(key, node);
|
||
name_to_node.insert(key.clone(), node);
|
||
}
|
||
|
||
// Reconstruct edges
|
||
for edge in &snapshot.topology.edges {
|
||
let source_node = name_to_node.get(&edge.source).ok_or_else(|| {
|
||
crate::error::ThermoError::DeserializationError(format!(
|
||
"Edge source '{}' not found in parameters",
|
||
edge.source
|
||
))
|
||
})?;
|
||
let target_node = name_to_node.get(&edge.target).ok_or_else(|| {
|
||
crate::error::ThermoError::DeserializationError(format!(
|
||
"Edge target '{}' not found in parameters",
|
||
edge.target
|
||
))
|
||
})?;
|
||
|
||
system.add_edge(*source_node, *target_node).map_err(|e| {
|
||
crate::error::ThermoError::DeserializationError(format!(
|
||
"Failed to add edge {} → {}: {:?}",
|
||
edge.source, edge.target, e
|
||
))
|
||
})?;
|
||
}
|
||
|
||
// Restore thermal couplings
|
||
for coupling in &snapshot.topology.thermal_couplings {
|
||
system.add_thermal_coupling(coupling.clone()).map_err(|e| {
|
||
crate::error::ThermoError::DeserializationError(format!(
|
||
"Failed to restore thermal coupling ({:?} → {:?}): {}",
|
||
coupling.hot_circuit, coupling.cold_circuit, e
|
||
))
|
||
})?;
|
||
}
|
||
|
||
// Restore constraints
|
||
for cs in &snapshot.constraints {
|
||
use crate::inverse::{ComponentOutput, Constraint, ConstraintId};
|
||
let output = match cs.output_type.as_str() {
|
||
"superheat" => ComponentOutput::superheat_for(&cs.component),
|
||
"subcooling" => ComponentOutput::subcooling_for(&cs.component),
|
||
"capacity" => ComponentOutput::capacity_for(&cs.component),
|
||
"heatTransferRate" => ComponentOutput::HeatTransferRate {
|
||
component_id: cs.component.clone(),
|
||
},
|
||
"massFlowRate" => ComponentOutput::MassFlowRate {
|
||
component_id: cs.component.clone(),
|
||
},
|
||
"pressure" => ComponentOutput::Pressure {
|
||
component_id: cs.component.clone(),
|
||
},
|
||
"temperature" => ComponentOutput::Temperature {
|
||
component_id: cs.component.clone(),
|
||
},
|
||
"saturationTemperature" => ComponentOutput::SaturationTemperature {
|
||
component_id: cs.component.clone(),
|
||
},
|
||
other => {
|
||
return Err(crate::error::ThermoError::DeserializationError(format!(
|
||
"Unknown constraint output type '{}' for component '{}'",
|
||
other, cs.component
|
||
)));
|
||
}
|
||
};
|
||
let id = ConstraintId::new(&cs.id);
|
||
let constraint = Constraint::new(id, output, cs.target);
|
||
system.add_constraint(constraint).map_err(|e| {
|
||
crate::error::ThermoError::DeserializationError(format!(
|
||
"Could not restore constraint '{}': {:?}",
|
||
cs.id, e
|
||
))
|
||
})?;
|
||
}
|
||
|
||
// Restore bounded variables
|
||
for bv in &snapshot.bounded_variables {
|
||
use crate::inverse::{BoundedVariable, BoundedVariableId};
|
||
let var = BoundedVariable::with_component(
|
||
BoundedVariableId::new(&bv.id),
|
||
&bv.component,
|
||
bv.initial_value,
|
||
bv.lower_bound,
|
||
bv.upper_bound,
|
||
)
|
||
.map_err(|e| {
|
||
crate::error::ThermoError::DeserializationError(format!(
|
||
"Failed to restore bounded variable '{}': {:?}",
|
||
bv.id, e
|
||
))
|
||
})?;
|
||
system.add_bounded_variable(var).map_err(|e| {
|
||
crate::error::ThermoError::DeserializationError(format!(
|
||
"Failed to add bounded variable '{}': {:?}",
|
||
bv.id, e
|
||
))
|
||
})?;
|
||
}
|
||
|
||
// Restore fluid state if present
|
||
if let Some(ref fluid_state) = snapshot.fluid_state {
|
||
tracing::debug!("Restoring fluid state: {} edges", fluid_state.edge_count());
|
||
// Fluid state is stored for hot-start scenarios.
|
||
// Apply to the system's internal state vector during solve initialization.
|
||
}
|
||
|
||
system.finalize().map_err(|e| {
|
||
crate::error::ThermoError::DeserializationError(format!(
|
||
"Failed to finalize reconstructed system: {:?}",
|
||
e
|
||
))
|
||
})?;
|
||
|
||
Ok(system)
|
||
}
|
||
|
||
/// Saves the system to a JSON file.
|
||
///
|
||
/// # Errors
|
||
///
|
||
/// Returns `ThermoError::IoError` if file writing fails.
|
||
///
|
||
/// # Examples
|
||
///
|
||
/// ```no_run
|
||
/// # use entropyk_solver::System;
|
||
/// # use std::path::Path;
|
||
/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
|
||
/// let system = System::new();
|
||
/// system.save_json(Path::new("system.json"))?;
|
||
/// # Ok(())
|
||
/// # }
|
||
/// ```
|
||
pub fn save_json<P: AsRef<std::path::Path>>(
|
||
&self,
|
||
path: P,
|
||
) -> Result<(), crate::error::ThermoError> {
|
||
use std::io::Write;
|
||
|
||
let json_str = self.to_json_string()?;
|
||
let path_ref = path.as_ref();
|
||
|
||
let mut file = std::fs::File::create(path_ref).map_err(|e| {
|
||
crate::error::ThermoError::IoError(format!("Failed to create file: {}", e))
|
||
})?;
|
||
|
||
file.write_all(json_str.as_bytes()).map_err(|e| {
|
||
crate::error::ThermoError::IoError(format!("Failed to write to file: {}", e))
|
||
})?;
|
||
|
||
tracing::info!("System saved to JSON file: {}", path_ref.display());
|
||
|
||
Ok(())
|
||
}
|
||
|
||
/// Loads a system from a JSON file.
|
||
///
|
||
/// # Errors
|
||
///
|
||
/// - `ThermoError::IoError` if file reading fails
|
||
/// - `ThermoError::DeserializationError` if JSON parsing fails
|
||
/// - See `from_json_string` for additional error conditions
|
||
///
|
||
/// # Examples
|
||
///
|
||
/// ```no_run
|
||
/// # use entropyk_solver::System;
|
||
/// # use std::path::Path;
|
||
/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
|
||
/// let system = System::load_json(Path::new("system.json"))?;
|
||
/// # Ok(())
|
||
/// # }
|
||
/// ```
|
||
pub fn load_json<P: AsRef<std::path::Path>>(
|
||
path: P,
|
||
) -> Result<Self, crate::error::ThermoError> {
|
||
use std::io::Read;
|
||
|
||
let path_ref = path.as_ref();
|
||
|
||
let mut file = std::fs::File::open(path_ref).map_err(|e| {
|
||
crate::error::ThermoError::IoError(format!("Failed to open file: {}", e))
|
||
})?;
|
||
|
||
let mut json_str = String::new();
|
||
file.read_to_string(&mut json_str).map_err(|e| {
|
||
crate::error::ThermoError::IoError(format!("Failed to read file: {}", e))
|
||
})?;
|
||
|
||
tracing::info!("System loaded from JSON file: {}", path_ref.display());
|
||
|
||
Self::from_json_string(&json_str)
|
||
}
|
||
}
|
||
|
||
impl Default for System {
|
||
fn default() -> Self {
|
||
Self::new()
|
||
}
|
||
}
|
||
|
||
#[cfg(test)]
|
||
mod tests {
|
||
use super::*;
|
||
use approx::assert_relative_eq;
|
||
use entropyk_components::port::{FluidId, Port};
|
||
use entropyk_components::{ConnectedPort, StateSlice};
|
||
use entropyk_core::{Enthalpy, Pressure};
|
||
|
||
/// Minimal mock component for testing.
|
||
struct MockComponent {
|
||
n_equations: usize,
|
||
}
|
||
|
||
impl Component for MockComponent {
|
||
fn compute_residuals(
|
||
&self,
|
||
_state: &StateSlice,
|
||
residuals: &mut entropyk_components::ResidualVector,
|
||
) -> Result<(), ComponentError> {
|
||
for r in residuals.iter_mut().take(self.n_equations) {
|
||
*r = 0.0;
|
||
}
|
||
Ok(())
|
||
}
|
||
|
||
fn jacobian_entries(
|
||
&self,
|
||
_state: &StateSlice,
|
||
jacobian: &mut JacobianBuilder,
|
||
) -> Result<(), ComponentError> {
|
||
for i in 0..self.n_equations {
|
||
jacobian.add_entry(i, i, 1.0);
|
||
}
|
||
Ok(())
|
||
}
|
||
|
||
fn n_equations(&self) -> usize {
|
||
self.n_equations
|
||
}
|
||
|
||
fn get_ports(&self) -> &[ConnectedPort] {
|
||
&[]
|
||
}
|
||
}
|
||
|
||
fn make_mock(n: usize) -> Box<dyn Component> {
|
||
Box::new(MockComponent { n_equations: n })
|
||
}
|
||
|
||
struct SaturatedControlMockComponent {
|
||
z_flow_idx: Option<usize>,
|
||
}
|
||
|
||
impl Component for SaturatedControlMockComponent {
|
||
fn compute_residuals(
|
||
&self,
|
||
state: &StateSlice,
|
||
residuals: &mut entropyk_components::ResidualVector,
|
||
) -> Result<(), ComponentError> {
|
||
let f_m = self.z_flow_idx.map(|idx| state[idx]).unwrap_or(1.0);
|
||
residuals[0] = state[0] - 2.0 * f_m;
|
||
Ok(())
|
||
}
|
||
|
||
fn jacobian_entries(
|
||
&self,
|
||
_state: &StateSlice,
|
||
jacobian: &mut JacobianBuilder,
|
||
) -> Result<(), ComponentError> {
|
||
jacobian.add_entry(0, 0, 1.0);
|
||
if let Some(idx) = self.z_flow_idx {
|
||
jacobian.add_entry(0, idx, -2.0);
|
||
}
|
||
Ok(())
|
||
}
|
||
|
||
fn n_equations(&self) -> usize {
|
||
1
|
||
}
|
||
|
||
fn get_ports(&self) -> &[ConnectedPort] {
|
||
&[]
|
||
}
|
||
|
||
fn set_calib_indices(&mut self, indices: entropyk_core::CalibIndices) {
|
||
self.z_flow_idx = indices.z_flow;
|
||
}
|
||
|
||
fn measure_output(
|
||
&self,
|
||
kind: entropyk_components::MeasuredOutput,
|
||
state: &StateSlice,
|
||
) -> Option<f64> {
|
||
match kind {
|
||
entropyk_components::MeasuredOutput::Capacity => {
|
||
Some(state[0] + 2.0 * state[1] + 3.0 * state[2])
|
||
}
|
||
// Extra linear output so override-network tests can exercise
|
||
// several distinct measurements on the same plant.
|
||
entropyk_components::MeasuredOutput::Temperature => Some(state[1]),
|
||
_ => None,
|
||
}
|
||
}
|
||
}
|
||
|
||
/// Mock component that reports fixed energy transfers, for cycle-performance tests.
|
||
struct EnergyMock {
|
||
heat_w: f64,
|
||
work_w: f64,
|
||
}
|
||
|
||
impl Component for EnergyMock {
|
||
fn compute_residuals(
|
||
&self,
|
||
_state: &StateSlice,
|
||
_residuals: &mut entropyk_components::ResidualVector,
|
||
) -> Result<(), ComponentError> {
|
||
Ok(())
|
||
}
|
||
|
||
fn jacobian_entries(
|
||
&self,
|
||
_state: &StateSlice,
|
||
_jacobian: &mut JacobianBuilder,
|
||
) -> Result<(), ComponentError> {
|
||
Ok(())
|
||
}
|
||
|
||
fn n_equations(&self) -> usize {
|
||
0
|
||
}
|
||
|
||
fn get_ports(&self) -> &[ConnectedPort] {
|
||
&[]
|
||
}
|
||
|
||
fn energy_transfers(
|
||
&self,
|
||
_state: &StateSlice,
|
||
) -> Option<(entropyk_core::Power, entropyk_core::Power)> {
|
||
Some((
|
||
entropyk_core::Power::from_watts(self.heat_w),
|
||
entropyk_core::Power::from_watts(self.work_w),
|
||
))
|
||
}
|
||
}
|
||
|
||
/// Mock component with 2 ports (inlet=0, outlet=1) for port validation tests.
|
||
fn make_ported_mock(fluid: &str, pressure_pa: f64, enthalpy_jkg: f64) -> Box<dyn Component> {
|
||
let inlet = Port::new(
|
||
FluidId::new(fluid),
|
||
Pressure::from_pascals(pressure_pa),
|
||
Enthalpy::from_joules_per_kg(enthalpy_jkg),
|
||
);
|
||
let outlet = Port::new(
|
||
FluidId::new(fluid),
|
||
Pressure::from_pascals(pressure_pa),
|
||
Enthalpy::from_joules_per_kg(enthalpy_jkg),
|
||
);
|
||
let (connected_inlet, connected_outlet) = inlet.connect(outlet).unwrap();
|
||
let ports: Vec<ConnectedPort> = vec![connected_inlet, connected_outlet];
|
||
Box::new(PortedMockComponent { ports })
|
||
}
|
||
|
||
struct PortedMockComponent {
|
||
ports: Vec<ConnectedPort>,
|
||
}
|
||
|
||
impl Component for PortedMockComponent {
|
||
fn compute_residuals(
|
||
&self,
|
||
_state: &StateSlice,
|
||
residuals: &mut entropyk_components::ResidualVector,
|
||
) -> Result<(), ComponentError> {
|
||
for r in residuals.iter_mut() {
|
||
*r = 0.0;
|
||
}
|
||
Ok(())
|
||
}
|
||
|
||
fn jacobian_entries(
|
||
&self,
|
||
_state: &StateSlice,
|
||
_jacobian: &mut JacobianBuilder,
|
||
) -> Result<(), ComponentError> {
|
||
Ok(())
|
||
}
|
||
|
||
fn n_equations(&self) -> usize {
|
||
0
|
||
}
|
||
|
||
fn get_ports(&self) -> &[ConnectedPort] {
|
||
&self.ports
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
fn test_simple_cycle_builds() {
|
||
let mut sys = System::new();
|
||
let n0 = sys.add_component(make_mock(0));
|
||
let n1 = sys.add_component(make_mock(0));
|
||
let n2 = sys.add_component(make_mock(0));
|
||
let n3 = sys.add_component(make_mock(0));
|
||
|
||
sys.add_edge(n0, n1).unwrap();
|
||
sys.add_edge(n1, n2).unwrap();
|
||
sys.add_edge(n2, n3).unwrap();
|
||
sys.add_edge(n3, n0).unwrap();
|
||
|
||
assert_eq!(sys.node_count(), 4);
|
||
assert_eq!(sys.edge_count(), 4);
|
||
|
||
let result = sys.finalize();
|
||
assert!(
|
||
result.is_ok(),
|
||
"finalize should succeed: {:?}",
|
||
result.err()
|
||
);
|
||
}
|
||
|
||
#[test]
|
||
fn test_state_vector_length() {
|
||
let mut sys = System::new();
|
||
let n0 = sys.add_component(make_mock(0));
|
||
let n1 = sys.add_component(make_mock(0));
|
||
sys.add_edge(n0, n1).unwrap();
|
||
sys.add_edge(n1, n0).unwrap();
|
||
|
||
sys.finalize().unwrap();
|
||
// CM1.4: 2-edge series cycle → 1 branch + 2×2 (P,h) = 5
|
||
assert_eq!(sys.state_vector_len(), 5); // 1 ṁ_branch + 4 (P,h per edge) = 5
|
||
}
|
||
|
||
#[test]
|
||
fn test_edge_kind_defaults_to_refrigerant() {
|
||
// AC#1: edges added via the existing API default to EdgeKind::Refrigerant.
|
||
let mut sys = System::new();
|
||
let n0 = sys.add_component(make_mock(0));
|
||
let n1 = sys.add_component(make_mock(0));
|
||
sys.add_edge(n0, n1).unwrap();
|
||
sys.add_edge(n1, n0).unwrap();
|
||
sys.finalize().unwrap();
|
||
|
||
for edge_idx in sys.edge_indices() {
|
||
assert_eq!(sys.edge_kind(edge_idx), EdgeKind::Refrigerant);
|
||
}
|
||
// EdgeKind::default() is Refrigerant.
|
||
assert_eq!(EdgeKind::default(), EdgeKind::Refrigerant);
|
||
}
|
||
|
||
#[test]
|
||
fn test_state_vector_len_is_branches_plus_2_edges() {
|
||
// CM1.4: state_vector_len() == n_branches + 2 * edge_count for a simple cycle
|
||
// with no control variables / thermal couplings / internal component state.
|
||
// A 3-edge series cycle forms 1 branch → 1 + 2×3 = 7.
|
||
let mut sys = System::new();
|
||
let n0 = sys.add_component(make_mock(0));
|
||
let n1 = sys.add_component(make_mock(0));
|
||
let n2 = sys.add_component(make_mock(0));
|
||
sys.add_edge(n0, n1).unwrap();
|
||
sys.add_edge(n1, n2).unwrap();
|
||
sys.add_edge(n2, n0).unwrap();
|
||
sys.finalize().unwrap();
|
||
|
||
assert_eq!(sys.edge_count(), 3);
|
||
// 1 branch (all 3 edges in series) + 2 × 3 (P, h per edge) = 7
|
||
assert_eq!(sys.state_vector_len(), 7);
|
||
}
|
||
|
||
#[test]
|
||
fn test_edge_state_mapping_distinct_in_range_ordered() {
|
||
// CM1.4: edge→state mapping: ṁ is SHARED per branch, P and h are per-edge.
|
||
// In a 3-edge series cycle: all edges share m_idx=0.
|
||
// P and h indices are unique per edge, packed consecutively starting at 1.
|
||
let mut sys = System::new();
|
||
let n0 = sys.add_component(make_mock(0));
|
||
let n1 = sys.add_component(make_mock(0));
|
||
let n2 = sys.add_component(make_mock(0));
|
||
sys.add_edge(n0, n1).unwrap();
|
||
sys.add_edge(n1, n2).unwrap();
|
||
sys.add_edge(n2, n0).unwrap();
|
||
sys.finalize().unwrap();
|
||
|
||
let len = sys.state_vector_len();
|
||
let mut seen_p = std::collections::HashSet::new();
|
||
let mut seen_h = std::collections::HashSet::new();
|
||
let mut seen_m = std::collections::HashSet::new();
|
||
|
||
for edge_idx in sys.edge_indices() {
|
||
let (m, p, h) = sys.edge_state_indices_full(edge_idx);
|
||
// ṁ index matches the dedicated accessor
|
||
assert_eq!(m, sys.edge_mass_flow_index(edge_idx));
|
||
// CM1.4: all edges in series share m_idx (same branch)
|
||
seen_m.insert(m);
|
||
// P < h (contiguous per edge)
|
||
assert!(p < h, "expected P < h, got ({p}, {h})");
|
||
assert_eq!(h, p + 1);
|
||
// In range
|
||
assert!(m < len, "m {m} out of range (len {len})");
|
||
assert!(p < len, "p {p} out of range (len {len})");
|
||
assert!(h < len, "h {h} out of range (len {len})");
|
||
// P and h are unique across edges
|
||
assert!(seen_p.insert(p), "P index {p} duplicated across edges");
|
||
assert!(seen_h.insert(h), "h index {h} duplicated across edges");
|
||
// Backwards-compatible (P, h) accessor agrees.
|
||
assert_eq!(sys.edge_state_indices(edge_idx), (p, h));
|
||
}
|
||
// A 3-edge series cycle has exactly 1 unique ṁ index (shared branch).
|
||
assert_eq!(seen_m.len(), 1, "3-edge series cycle must share 1 ṁ index");
|
||
}
|
||
|
||
#[test]
|
||
fn test_mass_flow_closure_count_is_zero() {
|
||
// Since CM1.3: mass_flow_closure_count() always returns 0 — closures removed.
|
||
let mut sys = System::new();
|
||
let n0 = sys.add_component(make_mock(0));
|
||
let n1 = sys.add_component(make_mock(0));
|
||
sys.add_edge(n0, n1).unwrap();
|
||
sys.add_edge(n1, n0).unwrap();
|
||
sys.finalize().unwrap();
|
||
|
||
assert_eq!(sys.mass_flow_closure_count(), 0);
|
||
}
|
||
|
||
#[test]
|
||
fn test_edge_indices_contiguous() {
|
||
// CM1.4: for a 3-edge series cycle (1 branch), state layout is:
|
||
// [0: ṁ_branch, 1: P_0, 2: h_0, 3: P_1, 4: h_1, 5: P_2, 6: h_2]
|
||
// All edges share m_idx=0. P and h are packed consecutively starting at 1.
|
||
let mut sys = System::new();
|
||
let n0 = sys.add_component(make_mock(0));
|
||
let n1 = sys.add_component(make_mock(0));
|
||
let n2 = sys.add_component(make_mock(0));
|
||
|
||
sys.add_edge(n0, n1).unwrap();
|
||
sys.add_edge(n1, n2).unwrap();
|
||
sys.add_edge(n2, n0).unwrap();
|
||
|
||
sys.finalize().unwrap();
|
||
|
||
let n = sys.edge_count(); // 3
|
||
// state_len = 1 branch + 2×3 edges = 7
|
||
assert_eq!(sys.state_vector_len(), 7);
|
||
|
||
// All edges share branch m_idx=0.
|
||
for edge_idx in sys.edge_indices() {
|
||
let (m, _p, _h) = sys.edge_state_indices_full(edge_idx);
|
||
assert_eq!(m, 0, "all edges in a series cycle share m_idx=0");
|
||
}
|
||
|
||
// P and h indices for the 3 edges should span slots 1..7 without gaps.
|
||
let mut ph_slots: Vec<usize> = sys
|
||
.edge_indices()
|
||
.flat_map(|e| {
|
||
let (_m, p, h) = sys.edge_state_indices_full(e);
|
||
[p, h]
|
||
})
|
||
.collect();
|
||
ph_slots.sort_unstable();
|
||
let expected_ph: Vec<usize> = (1..1 + 2 * n).collect();
|
||
assert_eq!(
|
||
ph_slots,
|
||
expected_ph,
|
||
"P,h slots must cover 1..{}",
|
||
1 + 2 * n
|
||
);
|
||
}
|
||
|
||
#[test]
|
||
fn test_cycle_detected() {
|
||
let mut sys = System::new();
|
||
let n0 = sys.add_component(make_mock(0));
|
||
let n1 = sys.add_component(make_mock(0));
|
||
let n2 = sys.add_component(make_mock(0));
|
||
|
||
sys.add_edge(n0, n1).unwrap();
|
||
sys.add_edge(n1, n2).unwrap();
|
||
sys.add_edge(n2, n0).unwrap();
|
||
|
||
sys.finalize().unwrap();
|
||
assert!(
|
||
sys.is_cyclic(),
|
||
"refrigeration cycle should be detected as cyclic"
|
||
);
|
||
}
|
||
|
||
#[test]
|
||
fn test_dangling_node_error() {
|
||
let mut sys = System::new();
|
||
let n0 = sys.add_component(make_mock(0));
|
||
let n1 = sys.add_component(make_mock(0));
|
||
let n2 = sys.add_component(make_mock(0)); // isolated
|
||
|
||
sys.add_edge(n0, n1).unwrap();
|
||
// n2 has no edges
|
||
|
||
let result = sys.finalize();
|
||
assert!(result.is_err());
|
||
match &result {
|
||
Err(TopologyError::IsolatedNode { node_index }) => {
|
||
assert!(
|
||
*node_index < sys.node_count(),
|
||
"isolated node index {} must be < node_count {}",
|
||
node_index,
|
||
sys.node_count()
|
||
);
|
||
assert_eq!(*node_index, n2.index(), "isolated node should be n2");
|
||
}
|
||
other => panic!("expected IsolatedNode error, got {:?}", other),
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
fn test_traverse_components() {
|
||
let mut sys = System::new();
|
||
let n0 = sys.add_component(make_mock(1));
|
||
let n1 = sys.add_component(make_mock(1));
|
||
sys.add_edge(n0, n1).unwrap();
|
||
sys.add_edge(n1, n0).unwrap();
|
||
|
||
sys.finalize().unwrap();
|
||
|
||
let mut count = 0;
|
||
for (_node_idx, component, edge_indices) in sys.traverse_for_jacobian() {
|
||
count += 1;
|
||
assert_eq!(component.n_equations(), 1);
|
||
assert_eq!(
|
||
edge_indices.len(),
|
||
2,
|
||
"each node has 2 incident edges in 2-node cycle"
|
||
);
|
||
for (_edge_idx, p, h) in &edge_indices {
|
||
assert!(p < &sys.state_vector_len());
|
||
assert!(h < &sys.state_vector_len());
|
||
assert_eq!(h, &(p + 1));
|
||
}
|
||
}
|
||
assert_eq!(count, 2);
|
||
}
|
||
|
||
#[test]
|
||
fn test_empty_graph_finalize_ok() {
|
||
let mut sys = System::new();
|
||
let result = sys.finalize();
|
||
assert!(result.is_ok());
|
||
}
|
||
|
||
#[test]
|
||
fn test_state_layout_integration() {
|
||
let mut sys = System::new();
|
||
let n0 = sys.add_component(make_mock(1));
|
||
let n1 = sys.add_component(make_mock(1));
|
||
sys.add_edge(n0, n1).unwrap();
|
||
sys.add_edge(n1, n0).unwrap();
|
||
|
||
sys.finalize().unwrap();
|
||
|
||
let layout = sys.state_layout();
|
||
assert!(layout.contains("P_edge"));
|
||
assert!(layout.contains("h_edge"));
|
||
assert!(layout.contains("m_edge"));
|
||
|
||
let state_len = sys.state_vector_len();
|
||
// CM1.4: 2-edge series cycle → 1 branch + 2×2 (P,h) = 5
|
||
assert_eq!(state_len, 5); // 1 ṁ_branch + 4 (P,h per edge)
|
||
|
||
// Layout: [0: ṁ_branch, 1: P_0, 2: h_0, 3: P_1, 4: h_1]
|
||
let mut state = vec![0.0; state_len];
|
||
state[0] = 0.05; // ṁ_branch
|
||
state[1] = 1e5; // P_edge0
|
||
state[2] = 250000.0; // h_edge0
|
||
state[3] = 5e5; // P_edge1
|
||
state[4] = 300000.0; // h_edge1
|
||
|
||
// 2 component equations (each mock has 1 equation)
|
||
let mut residuals = vec![0.0; 2];
|
||
let result = sys.compute_residuals(&state, &mut residuals);
|
||
assert!(result.is_ok());
|
||
}
|
||
|
||
#[test]
|
||
fn test_valid_connection_same_fluid() {
|
||
let p = 100_000.0;
|
||
let h = 400_000.0;
|
||
let mut sys = System::new();
|
||
let n0 = sys.add_component(make_ported_mock("R134a", p, h));
|
||
let n1 = sys.add_component(make_ported_mock("R134a", p, h));
|
||
|
||
let result = sys.add_edge_with_ports(n0, 1, n1, 0);
|
||
assert!(
|
||
result.is_ok(),
|
||
"R134a to R134a should succeed: {:?}",
|
||
result.err()
|
||
);
|
||
sys.add_edge(n1, n0).unwrap(); // backward edge, no ports
|
||
sys.finalize().unwrap();
|
||
}
|
||
|
||
#[test]
|
||
fn test_incompatible_fluid_rejected() {
|
||
let p = 100_000.0;
|
||
let h = 400_000.0;
|
||
let mut sys = System::new();
|
||
let n0 = sys.add_component(make_ported_mock("R134a", p, h));
|
||
let n1 = sys.add_component(make_ported_mock("Water", p, h));
|
||
|
||
let result = sys.add_edge_with_ports(n0, 1, n1, 0);
|
||
assert!(result.is_err());
|
||
match result {
|
||
Err(AddEdgeError::Connection(ConnectionError::IncompatibleFluid { from, to })) => {
|
||
assert_eq!(from, "R134a");
|
||
assert_eq!(to, "Water");
|
||
}
|
||
other => panic!(
|
||
"expected AddEdgeError::Connection(IncompatibleFluid), got {:?}",
|
||
other
|
||
),
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
fn test_pressure_mismatch_rejected() {
|
||
let h = 400_000.0;
|
||
let mut sys = System::new();
|
||
let n0 = sys.add_component(make_ported_mock("R134a", 100_000.0, h));
|
||
let n1 = sys.add_component(make_ported_mock("R134a", 200_000.0, h));
|
||
|
||
let result = sys.add_edge_with_ports(n0, 1, n1, 0);
|
||
assert!(result.is_err());
|
||
match result {
|
||
Err(AddEdgeError::Connection(ConnectionError::PressureMismatch {
|
||
from_pressure,
|
||
to_pressure,
|
||
tolerance,
|
||
})) => {
|
||
assert_relative_eq!(from_pressure, 100_000.0, epsilon = 1.0);
|
||
assert_relative_eq!(to_pressure, 200_000.0, epsilon = 1.0);
|
||
assert!(tolerance >= 1.0, "tolerance should be at least 1 Pa");
|
||
}
|
||
other => panic!(
|
||
"expected AddEdgeError::Connection(PressureMismatch), got {:?}",
|
||
other
|
||
),
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
fn test_enthalpy_mismatch_rejected() {
|
||
let p = 100_000.0;
|
||
let mut sys = System::new();
|
||
let n0 = sys.add_component(make_ported_mock("R134a", p, 400_000.0));
|
||
let n1 = sys.add_component(make_ported_mock("R134a", p, 500_000.0));
|
||
|
||
let result = sys.add_edge_with_ports(n0, 1, n1, 0);
|
||
assert!(result.is_err());
|
||
match result {
|
||
Err(AddEdgeError::Connection(ConnectionError::EnthalpyMismatch {
|
||
from_enthalpy,
|
||
to_enthalpy,
|
||
tolerance,
|
||
})) => {
|
||
assert_relative_eq!(from_enthalpy, 400_000.0, epsilon = 1.0);
|
||
assert_relative_eq!(to_enthalpy, 500_000.0, epsilon = 1.0);
|
||
assert_relative_eq!(tolerance, 100.0, epsilon = 1.0); // ENTHALPY_TOLERANCE_J_KG
|
||
}
|
||
other => panic!(
|
||
"expected AddEdgeError::Connection(EnthalpyMismatch), got {:?}",
|
||
other
|
||
),
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
fn test_pressure_tolerance_boundary() {
|
||
let h = 400_000.0;
|
||
let base_pressure = 100_000.0; // 100 kPa
|
||
let tolerance: f64 = (base_pressure * 1e-4f64).max(1.0f64); // 10 Pa for 100 kPa
|
||
|
||
let mut sys = System::new();
|
||
let n0 = sys.add_component(make_ported_mock("R134a", base_pressure, h));
|
||
|
||
// Exactly at tolerance - should succeed
|
||
let n1 = sys.add_component(make_ported_mock("R134a", base_pressure + tolerance, h));
|
||
let result = sys.add_edge_with_ports(n0, 1, n1, 0);
|
||
assert!(
|
||
result.is_ok(),
|
||
"Connection at exact tolerance ({:.1} Pa diff) should succeed",
|
||
tolerance
|
||
);
|
||
|
||
// Just outside tolerance - should fail
|
||
let n2 = sys.add_component(make_ported_mock(
|
||
"R134a",
|
||
base_pressure + tolerance + 1.0,
|
||
h,
|
||
));
|
||
let result = sys.add_edge_with_ports(n0, 1, n2, 0);
|
||
assert!(
|
||
result.is_err(),
|
||
"Connection just outside tolerance ({:.1} Pa diff) should fail",
|
||
tolerance + 1.0
|
||
);
|
||
match result {
|
||
Err(AddEdgeError::Connection(ConnectionError::PressureMismatch {
|
||
from_pressure,
|
||
to_pressure,
|
||
tolerance: tol,
|
||
})) => {
|
||
assert_relative_eq!(from_pressure, base_pressure, epsilon = 0.1);
|
||
assert_relative_eq!(to_pressure, base_pressure + tolerance + 1.0, epsilon = 0.1);
|
||
assert_relative_eq!(tol, tolerance, epsilon = 0.1);
|
||
}
|
||
other => panic!("expected PressureMismatch at boundary, got {:?}", other),
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
fn test_invalid_port_index_rejected() {
|
||
let p = 100_000.0;
|
||
let h = 400_000.0;
|
||
let mut sys = System::new();
|
||
let n0 = sys.add_component(make_ported_mock("R134a", p, h));
|
||
let n1 = sys.add_component(make_ported_mock("R134a", p, h));
|
||
|
||
// Port index 2 out of bounds for 2-port component
|
||
let result = sys.add_edge_with_ports(n0, 2, n1, 0);
|
||
assert!(result.is_err());
|
||
match result {
|
||
Err(AddEdgeError::Connection(ConnectionError::InvalidPortIndex {
|
||
index,
|
||
port_count,
|
||
max_index,
|
||
})) => {
|
||
assert_eq!(index, 2);
|
||
assert_eq!(port_count, 2);
|
||
assert_eq!(max_index, 1);
|
||
}
|
||
other => panic!("expected InvalidPortIndex for source, got {:?}", other),
|
||
}
|
||
|
||
// Target port index out of bounds
|
||
let result = sys.add_edge_with_ports(n0, 1, n1, 5);
|
||
assert!(result.is_err());
|
||
match result {
|
||
Err(AddEdgeError::Connection(ConnectionError::InvalidPortIndex {
|
||
index,
|
||
port_count,
|
||
max_index,
|
||
})) => {
|
||
assert_eq!(index, 5);
|
||
assert_eq!(port_count, 2);
|
||
assert_eq!(max_index, 1);
|
||
}
|
||
other => panic!("expected InvalidPortIndex for target, got {:?}", other),
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
fn test_simple_cycle_port_validation() {
|
||
let p = 100_000.0;
|
||
let h = 400_000.0;
|
||
let mut sys = System::new();
|
||
let n0 = sys.add_component(make_ported_mock("R134a", p, h));
|
||
let n1 = sys.add_component(make_ported_mock("R134a", p, h));
|
||
let n2 = sys.add_component(make_ported_mock("R134a", p, h));
|
||
let n3 = sys.add_component(make_ported_mock("R134a", p, h));
|
||
|
||
sys.add_edge_with_ports(n0, 1, n1, 0).unwrap();
|
||
sys.add_edge_with_ports(n1, 1, n2, 0).unwrap();
|
||
sys.add_edge_with_ports(n2, 1, n3, 0).unwrap();
|
||
sys.add_edge_with_ports(n3, 1, n0, 0).unwrap();
|
||
|
||
assert_eq!(sys.edge_count(), 4);
|
||
sys.finalize().unwrap();
|
||
}
|
||
|
||
#[test]
|
||
fn test_compute_residuals_bounds_check() {
|
||
let mut sys = System::new();
|
||
let n0 = sys.add_component(make_mock(2));
|
||
let n1 = sys.add_component(make_mock(2));
|
||
sys.add_edge(n0, n1).unwrap();
|
||
sys.add_edge(n1, n0).unwrap();
|
||
sys.finalize().unwrap();
|
||
|
||
let state = vec![0.0; sys.state_vector_len()];
|
||
let mut residuals = vec![0.0; 1]; // Too small: need 4 (2 components × 2 equations, no closures)
|
||
let result = sys.compute_residuals(&state, &mut residuals);
|
||
assert!(result.is_err());
|
||
match result {
|
||
Err(ComponentError::InvalidResidualDimensions { expected, actual }) => {
|
||
assert_eq!(expected, 4);
|
||
assert_eq!(actual, 1);
|
||
}
|
||
other => panic!("expected InvalidResidualDimensions, got {:?}", other),
|
||
}
|
||
}
|
||
|
||
// --- Story 3.3: Multi-Circuit Machine Definition tests ---
|
||
|
||
#[test]
|
||
fn test_two_circuit_machine() {
|
||
let mut sys = System::new();
|
||
let c0 = CircuitId::ZERO;
|
||
let c1 = CircuitId(1);
|
||
|
||
let n0 = sys.add_component_to_circuit(make_mock(0), c0).unwrap();
|
||
let n1 = sys.add_component_to_circuit(make_mock(0), c0).unwrap();
|
||
let n2 = sys.add_component_to_circuit(make_mock(0), c1).unwrap();
|
||
let n3 = sys.add_component_to_circuit(make_mock(0), c1).unwrap();
|
||
|
||
sys.add_edge(n0, n1).unwrap();
|
||
sys.add_edge(n1, n0).unwrap();
|
||
sys.add_edge(n2, n3).unwrap();
|
||
sys.add_edge(n3, n2).unwrap();
|
||
|
||
assert_eq!(sys.circuit_count(), 2);
|
||
assert_eq!(sys.circuit_nodes(c0).count(), 2);
|
||
assert_eq!(sys.circuit_nodes(c1).count(), 2);
|
||
assert_eq!(sys.circuit_edges(c0).count(), 2);
|
||
assert_eq!(sys.circuit_edges(c1).count(), 2);
|
||
|
||
sys.finalize().unwrap();
|
||
}
|
||
|
||
#[test]
|
||
fn test_cross_circuit_edge_rejected() {
|
||
let mut sys = System::new();
|
||
let c0 = CircuitId::ZERO;
|
||
let c1 = CircuitId(1);
|
||
|
||
let n0 = sys.add_component_to_circuit(make_mock(0), c0).unwrap();
|
||
let n1 = sys.add_component_to_circuit(make_mock(0), c1).unwrap();
|
||
|
||
let result = sys.add_edge(n0, n1);
|
||
assert!(result.is_err());
|
||
match result {
|
||
Err(TopologyError::CrossCircuitConnection {
|
||
source_circuit,
|
||
target_circuit,
|
||
}) => {
|
||
assert_eq!(source_circuit, 0);
|
||
assert_eq!(target_circuit, 1);
|
||
}
|
||
other => panic!("expected CrossCircuitConnection, got {:?}", other),
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
fn test_circuit_count_and_accessors() {
|
||
let mut sys = System::new();
|
||
let c0 = CircuitId::ZERO;
|
||
let c1 = CircuitId(1);
|
||
let c2 = CircuitId(2);
|
||
|
||
let _n0 = sys.add_component_to_circuit(make_mock(0), c0).unwrap();
|
||
let _n1 = sys.add_component_to_circuit(make_mock(0), c0).unwrap();
|
||
let _n2 = sys.add_component_to_circuit(make_mock(0), c1).unwrap();
|
||
let _n3 = sys.add_component_to_circuit(make_mock(0), c2).unwrap();
|
||
|
||
assert_eq!(sys.circuit_count(), 3);
|
||
assert_eq!(sys.circuit_nodes(c0).count(), 2);
|
||
assert_eq!(sys.circuit_nodes(c1).count(), 1);
|
||
assert_eq!(sys.circuit_nodes(c2).count(), 1);
|
||
}
|
||
|
||
#[test]
|
||
fn test_max_five_circuits() {
|
||
// Test: 5 circuits accepted (0, 1, 2, 3, 4), 6th circuit (5) rejected
|
||
let mut sys = System::new();
|
||
for i in 0..=4 {
|
||
let cid = CircuitId(i);
|
||
let result = sys.add_component_to_circuit(make_mock(0), cid);
|
||
assert!(
|
||
result.is_ok(),
|
||
"circuit {} should be accepted (max 5 circuits: 0-4)",
|
||
i
|
||
);
|
||
}
|
||
assert_eq!(sys.circuit_count(), 5, "should have exactly 5 circuits");
|
||
|
||
// 6th circuit should be rejected
|
||
let result = sys.add_component_to_circuit(make_mock(0), CircuitId(5));
|
||
assert!(
|
||
result.is_err(),
|
||
"circuit 5 should be rejected (exceeds max of 4)"
|
||
);
|
||
match result {
|
||
Err(TopologyError::TooManyCircuits { requested }) => assert_eq!(requested, 5),
|
||
other => panic!("expected TooManyCircuits, got {:?}", other),
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
fn test_single_circuit_backward_compat() {
|
||
let mut sys = System::new();
|
||
let n0 = sys.add_component(make_mock(0));
|
||
let n1 = sys.add_component(make_mock(0));
|
||
let edge0 = sys.add_edge(n0, n1).unwrap();
|
||
let edge1 = sys.add_edge(n1, n0).unwrap();
|
||
|
||
assert_eq!(sys.circuit_count(), 1);
|
||
assert_eq!(sys.circuit_nodes(CircuitId::ZERO).count(), 2);
|
||
|
||
// Verify edge circuit membership
|
||
assert_eq!(sys.edge_circuit(edge0), CircuitId::ZERO);
|
||
assert_eq!(sys.edge_circuit(edge1), CircuitId::ZERO);
|
||
|
||
// Verify circuit_edges returns correct edges
|
||
let circuit_0_edges: Vec<_> = sys.circuit_edges(CircuitId::ZERO).collect();
|
||
assert_eq!(circuit_0_edges.len(), 2);
|
||
assert!(circuit_0_edges.contains(&edge0));
|
||
assert!(circuit_0_edges.contains(&edge1));
|
||
|
||
sys.finalize().unwrap();
|
||
}
|
||
|
||
#[test]
|
||
fn test_cross_circuit_add_edge_with_ports_rejected() {
|
||
let p = 100_000.0;
|
||
let h = 400_000.0;
|
||
let mut sys = System::new();
|
||
let n0 = sys
|
||
.add_component_to_circuit(make_ported_mock("R134a", p, h), CircuitId::ZERO)
|
||
.unwrap();
|
||
let n1 = sys
|
||
.add_component_to_circuit(make_ported_mock("R134a", p, h), CircuitId(1))
|
||
.unwrap();
|
||
|
||
let result = sys.add_edge_with_ports(n0, 1, n1, 0);
|
||
assert!(result.is_err());
|
||
match result {
|
||
Err(AddEdgeError::Topology(TopologyError::CrossCircuitConnection {
|
||
source_circuit,
|
||
target_circuit,
|
||
})) => {
|
||
assert_eq!(source_circuit, 0);
|
||
assert_eq!(target_circuit, 1);
|
||
}
|
||
other => panic!(
|
||
"expected AddEdgeError::Topology(CrossCircuitConnection), got {:?}",
|
||
other
|
||
),
|
||
}
|
||
}
|
||
|
||
// --- Story 3.4: Thermal Coupling Between Circuits tests ---
|
||
|
||
#[test]
|
||
fn test_add_thermal_coupling_valid() {
|
||
use entropyk_core::ThermalConductance;
|
||
|
||
let mut sys = System::new();
|
||
let _n0 = sys
|
||
.add_component_to_circuit(make_mock(0), CircuitId(0))
|
||
.unwrap();
|
||
let _n1 = sys
|
||
.add_component_to_circuit(make_mock(0), CircuitId(1))
|
||
.unwrap();
|
||
|
||
let coupling = ThermalCoupling::new(
|
||
CircuitId(0),
|
||
CircuitId(1),
|
||
ThermalConductance::from_watts_per_kelvin(1000.0),
|
||
);
|
||
|
||
let idx = sys.add_thermal_coupling(coupling).unwrap();
|
||
assert_eq!(idx, 0);
|
||
assert_eq!(sys.thermal_coupling_count(), 1);
|
||
|
||
let retrieved = sys.get_thermal_coupling(0).unwrap();
|
||
assert_eq!(retrieved.hot_circuit, CircuitId(0));
|
||
assert_eq!(retrieved.cold_circuit, CircuitId(1));
|
||
}
|
||
|
||
#[test]
|
||
fn test_add_thermal_coupling_invalid_circuit() {
|
||
use entropyk_core::ThermalConductance;
|
||
|
||
let mut sys = System::new();
|
||
let _n0 = sys
|
||
.add_component_to_circuit(make_mock(0), CircuitId(0))
|
||
.unwrap();
|
||
// Circuit 1 has no components
|
||
|
||
let coupling = ThermalCoupling::new(
|
||
CircuitId(0),
|
||
CircuitId(1), // This circuit doesn't exist
|
||
ThermalConductance::from_watts_per_kelvin(1000.0),
|
||
);
|
||
|
||
let result = sys.add_thermal_coupling(coupling);
|
||
assert!(result.is_err());
|
||
match result {
|
||
Err(TopologyError::InvalidCircuitForCoupling { circuit_id }) => {
|
||
assert_eq!(circuit_id, 1);
|
||
}
|
||
other => panic!("expected InvalidCircuitForCoupling, got {:?}", other),
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
fn test_add_thermal_coupling_hot_circuit_invalid() {
|
||
use entropyk_core::ThermalConductance;
|
||
|
||
let mut sys = System::new();
|
||
let _n0 = sys
|
||
.add_component_to_circuit(make_mock(0), CircuitId(0))
|
||
.unwrap();
|
||
|
||
let coupling = ThermalCoupling::new(
|
||
CircuitId(99), // This circuit doesn't exist
|
||
CircuitId(0),
|
||
ThermalConductance::from_watts_per_kelvin(1000.0),
|
||
);
|
||
|
||
let result = sys.add_thermal_coupling(coupling);
|
||
assert!(result.is_err());
|
||
match result {
|
||
Err(TopologyError::InvalidCircuitForCoupling { circuit_id }) => {
|
||
assert_eq!(circuit_id, 99);
|
||
}
|
||
other => panic!("expected InvalidCircuitForCoupling, got {:?}", other),
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
fn test_multiple_thermal_couplings() {
|
||
use entropyk_core::ThermalConductance;
|
||
|
||
let mut sys = System::new();
|
||
let _n0 = sys
|
||
.add_component_to_circuit(make_mock(0), CircuitId(0))
|
||
.unwrap();
|
||
let _n1 = sys
|
||
.add_component_to_circuit(make_mock(0), CircuitId(1))
|
||
.unwrap();
|
||
let _n2 = sys
|
||
.add_component_to_circuit(make_mock(0), CircuitId(2))
|
||
.unwrap();
|
||
|
||
let coupling1 = ThermalCoupling::new(
|
||
CircuitId(0),
|
||
CircuitId(1),
|
||
ThermalConductance::from_watts_per_kelvin(1000.0),
|
||
);
|
||
let coupling2 = ThermalCoupling::new(
|
||
CircuitId(1),
|
||
CircuitId(2),
|
||
ThermalConductance::from_watts_per_kelvin(500.0),
|
||
);
|
||
|
||
let idx1 = sys.add_thermal_coupling(coupling1).unwrap();
|
||
let idx2 = sys.add_thermal_coupling(coupling2).unwrap();
|
||
|
||
assert_eq!(idx1, 0);
|
||
assert_eq!(idx2, 1);
|
||
assert_eq!(sys.thermal_coupling_count(), 2);
|
||
|
||
let all_couplings = sys.thermal_couplings();
|
||
assert_eq!(all_couplings.len(), 2);
|
||
}
|
||
|
||
#[test]
|
||
fn test_thermal_coupling_same_circuit() {
|
||
// It's valid to couple a circuit to itself (internal heat exchanger / economizer)
|
||
use entropyk_core::ThermalConductance;
|
||
|
||
let mut sys = System::new();
|
||
let _n0 = sys
|
||
.add_component_to_circuit(make_mock(0), CircuitId(0))
|
||
.unwrap();
|
||
let _n1 = sys
|
||
.add_component_to_circuit(make_mock(0), CircuitId(0))
|
||
.unwrap();
|
||
|
||
let coupling = ThermalCoupling::new(
|
||
CircuitId(0),
|
||
CircuitId(0), // Same circuit
|
||
ThermalConductance::from_watts_per_kelvin(1000.0),
|
||
);
|
||
|
||
let result = sys.add_thermal_coupling(coupling);
|
||
assert!(result.is_ok(), "Same-circuit coupling should be allowed");
|
||
}
|
||
|
||
// --- Story 3.5: Zero-Flow Branch Handling ---
|
||
|
||
/// Mock component that behaves like an Off branch: residual = state[0] (ṁ - 0), Jacobian ∂r/∂state[0] = 1.
|
||
/// At state[0] = 0 (zero flow) residuals and Jacobian remain finite (no division by zero).
|
||
struct ZeroFlowMock;
|
||
|
||
impl Component for ZeroFlowMock {
|
||
fn compute_residuals(
|
||
&self,
|
||
state: &StateSlice,
|
||
residuals: &mut entropyk_components::ResidualVector,
|
||
) -> Result<(), ComponentError> {
|
||
if !state.is_empty() {
|
||
residuals[0] = state[0];
|
||
}
|
||
Ok(())
|
||
}
|
||
|
||
fn jacobian_entries(
|
||
&self,
|
||
_state: &StateSlice,
|
||
jacobian: &mut JacobianBuilder,
|
||
) -> Result<(), ComponentError> {
|
||
jacobian.add_entry(0, 0, 1.0);
|
||
Ok(())
|
||
}
|
||
|
||
fn n_equations(&self) -> usize {
|
||
1
|
||
}
|
||
|
||
fn get_ports(&self) -> &[ConnectedPort] {
|
||
&[]
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
fn test_zero_flow_branch_residuals_and_jacobian_finite() {
|
||
// Story 3.5: System with one branch at zero flow must not produce NaN/Inf in residuals or Jacobian.
|
||
let mut sys = System::new();
|
||
let n0 = sys.add_component(Box::new(ZeroFlowMock));
|
||
let n1 = sys.add_component(make_mock(1));
|
||
sys.add_edge(n0, n1).unwrap();
|
||
sys.add_edge(n1, n0).unwrap();
|
||
sys.finalize().unwrap();
|
||
|
||
let state_len = sys.state_vector_len();
|
||
let mut state = vec![0.0; state_len];
|
||
state[0] = 0.0;
|
||
|
||
let total_eqs: usize = 1 + 1; // ZeroFlowMock (1 eq) + make_mock(1) (1 eq)
|
||
let mut residuals = vec![0.0; total_eqs];
|
||
let result = sys.compute_residuals(&state, &mut residuals);
|
||
assert!(result.is_ok());
|
||
for (i, &r) in residuals.iter().enumerate() {
|
||
assert!(r.is_finite(), "residual[{}] must be finite, got {}", i, r);
|
||
}
|
||
|
||
let mut jacobian = JacobianBuilder::new();
|
||
let result = sys.assemble_jacobian(&state, &mut jacobian);
|
||
assert!(result.is_ok());
|
||
|
||
let entries = jacobian.entries();
|
||
for (row, col, value) in entries {
|
||
assert!(
|
||
value.is_finite(),
|
||
"Jacobian ({}, {}) must be finite, got {}",
|
||
row,
|
||
col,
|
||
value
|
||
);
|
||
}
|
||
|
||
// Check for zero rows: each equation should have at least one non-zero derivative
|
||
let mut row_has_nonzero = vec![false; total_eqs];
|
||
for (row, _col, value) in entries {
|
||
if value.abs() > 1e-15 {
|
||
row_has_nonzero[*row] = true;
|
||
}
|
||
}
|
||
for (row, &has_nonzero) in row_has_nonzero.iter().enumerate() {
|
||
assert!(
|
||
has_nonzero,
|
||
"Jacobian row {} is all zeros (degenerate equation)",
|
||
row
|
||
);
|
||
}
|
||
}
|
||
|
||
// ────────────────────────────────────────────────────────────────────────
|
||
// Constraint Management Tests
|
||
// ────────────────────────────────────────────────────────────────────────
|
||
|
||
#[test]
|
||
fn test_add_constraint() {
|
||
use crate::inverse::{ComponentOutput, Constraint, ConstraintId};
|
||
|
||
let mut system = System::new();
|
||
|
||
// Register a component first
|
||
let node = system.add_component(make_mock(0));
|
||
system.register_component_name("evaporator", node);
|
||
|
||
let constraint = Constraint::new(
|
||
ConstraintId::new("superheat_control"),
|
||
ComponentOutput::Superheat {
|
||
component_id: "evaporator".to_string(),
|
||
},
|
||
5.0,
|
||
);
|
||
|
||
assert!(system.add_constraint(constraint).is_ok());
|
||
assert_eq!(system.constraint_count(), 1);
|
||
}
|
||
|
||
#[test]
|
||
fn test_add_constraint_unregistered_component() {
|
||
use crate::inverse::{ComponentOutput, Constraint, ConstraintId};
|
||
|
||
let mut system = System::new();
|
||
// No component registered with name "evaporator"
|
||
|
||
let constraint = Constraint::new(
|
||
ConstraintId::new("superheat_control"),
|
||
ComponentOutput::Superheat {
|
||
component_id: "evaporator".to_string(),
|
||
},
|
||
5.0,
|
||
);
|
||
|
||
let result = system.add_constraint(constraint);
|
||
assert!(matches!(
|
||
result,
|
||
Err(ConstraintError::InvalidReference { .. })
|
||
));
|
||
assert_eq!(system.constraint_count(), 0);
|
||
}
|
||
|
||
#[test]
|
||
fn test_register_component_name() {
|
||
let mut system = System::new();
|
||
let n0 = system.add_component(make_mock(0));
|
||
let n1 = system.add_component(make_mock(0));
|
||
|
||
// Register first name - should return true (newly registered)
|
||
assert!(system.register_component_name("evaporator", n0));
|
||
|
||
// Register second name - should return true
|
||
assert!(system.register_component_name("condenser", n1));
|
||
|
||
// Re-register same name with different node - should return false
|
||
assert!(!system.register_component_name("evaporator", n1));
|
||
|
||
// Verify lookup works
|
||
assert_eq!(system.get_component_node("evaporator"), Some(n1)); // Updated to n1
|
||
assert_eq!(system.get_component_node("condenser"), Some(n1));
|
||
assert_eq!(system.get_component_node("nonexistent"), None);
|
||
|
||
// Verify iteration
|
||
let names: Vec<&str> = system.registered_component_names().collect();
|
||
assert_eq!(names.len(), 2);
|
||
}
|
||
|
||
#[test]
|
||
fn test_add_duplicate_constraint() {
|
||
use crate::inverse::{ComponentOutput, Constraint, ConstraintId};
|
||
|
||
let mut system = System::new();
|
||
|
||
// Register components
|
||
let n0 = system.add_component(make_mock(0));
|
||
let n1 = system.add_component(make_mock(0));
|
||
system.register_component_name("evaporator", n0);
|
||
system.register_component_name("condenser", n1);
|
||
|
||
let constraint1 = Constraint::new(
|
||
ConstraintId::new("superheat"),
|
||
ComponentOutput::Superheat {
|
||
component_id: "evaporator".to_string(),
|
||
},
|
||
5.0,
|
||
);
|
||
let constraint2 = Constraint::new(
|
||
ConstraintId::new("superheat"),
|
||
ComponentOutput::Superheat {
|
||
component_id: "condenser".to_string(),
|
||
},
|
||
3.0,
|
||
);
|
||
|
||
assert!(system.add_constraint(constraint1).is_ok());
|
||
let result = system.add_constraint(constraint2);
|
||
assert!(matches!(result, Err(ConstraintError::DuplicateId { .. })));
|
||
assert_eq!(system.constraint_count(), 1);
|
||
}
|
||
|
||
#[test]
|
||
fn test_remove_constraint() {
|
||
use crate::inverse::{ComponentOutput, Constraint, ConstraintId};
|
||
|
||
let mut system = System::new();
|
||
|
||
// Register component
|
||
let node = system.add_component(make_mock(0));
|
||
system.register_component_name("evaporator", node);
|
||
|
||
let id = ConstraintId::new("superheat");
|
||
let constraint = Constraint::new(
|
||
id.clone(),
|
||
ComponentOutput::Superheat {
|
||
component_id: "evaporator".to_string(),
|
||
},
|
||
5.0,
|
||
);
|
||
|
||
system.add_constraint(constraint).unwrap();
|
||
assert_eq!(system.constraint_count(), 1);
|
||
|
||
let removed = system.remove_constraint(&id);
|
||
assert!(removed.is_some());
|
||
assert_eq!(system.constraint_count(), 0);
|
||
|
||
// Removing non-existent constraint returns None
|
||
let removed_again = system.remove_constraint(&id);
|
||
assert!(removed_again.is_none());
|
||
}
|
||
|
||
#[test]
|
||
fn test_get_constraint() {
|
||
use crate::inverse::{ComponentOutput, Constraint, ConstraintId};
|
||
|
||
let mut system = System::new();
|
||
|
||
// Register component
|
||
let node = system.add_component(make_mock(0));
|
||
system.register_component_name("compressor", node);
|
||
|
||
let id = ConstraintId::new("pressure_control");
|
||
let constraint = Constraint::new(
|
||
id.clone(),
|
||
ComponentOutput::Pressure {
|
||
component_id: "compressor".to_string(),
|
||
},
|
||
300000.0,
|
||
);
|
||
|
||
system.add_constraint(constraint).unwrap();
|
||
|
||
let retrieved = system.get_constraint(&id);
|
||
assert!(retrieved.is_some());
|
||
assert_relative_eq!(retrieved.unwrap().target_value(), 300000.0);
|
||
|
||
let missing = system.get_constraint(&ConstraintId::new("nonexistent"));
|
||
assert!(missing.is_none());
|
||
}
|
||
|
||
#[test]
|
||
fn test_multiple_constraints() {
|
||
use crate::inverse::{ComponentOutput, Constraint, ConstraintId};
|
||
|
||
let mut system = System::new();
|
||
|
||
// Register components
|
||
let n0 = system.add_component(make_mock(0));
|
||
let n1 = system.add_component(make_mock(0));
|
||
system.register_component_name("evaporator", n0);
|
||
system.register_component_name("condenser", n1);
|
||
|
||
let c1 = Constraint::new(
|
||
ConstraintId::new("superheat"),
|
||
ComponentOutput::Superheat {
|
||
component_id: "evaporator".to_string(),
|
||
},
|
||
5.0,
|
||
);
|
||
let c2 = Constraint::new(
|
||
ConstraintId::new("subcooling"),
|
||
ComponentOutput::Subcooling {
|
||
component_id: "condenser".to_string(),
|
||
},
|
||
3.0,
|
||
);
|
||
let c3 = Constraint::new(
|
||
ConstraintId::new("capacity"),
|
||
ComponentOutput::HeatTransferRate {
|
||
component_id: "evaporator".to_string(),
|
||
},
|
||
10000.0,
|
||
);
|
||
|
||
assert!(system.add_constraint(c1).is_ok());
|
||
assert!(system.add_constraint(c2).is_ok());
|
||
assert!(system.add_constraint(c3).is_ok());
|
||
|
||
assert_eq!(system.constraint_count(), 3);
|
||
assert_eq!(system.constraint_residual_count(), 3);
|
||
|
||
// Iterate over all constraints
|
||
let count = system.constraints().count();
|
||
assert_eq!(count, 3);
|
||
}
|
||
|
||
// ────────────────────────────────────────────────────────────────────────
|
||
// Bounded Variable Tests
|
||
// ────────────────────────────────────────────────────────────────────────
|
||
|
||
#[test]
|
||
fn test_add_bounded_variable() {
|
||
use crate::inverse::{BoundedVariable, BoundedVariableId};
|
||
|
||
let mut system = System::new();
|
||
|
||
let valve =
|
||
BoundedVariable::new(BoundedVariableId::new("expansion_valve"), 0.5, 0.0, 1.0).unwrap();
|
||
|
||
assert!(system.add_bounded_variable(valve).is_ok());
|
||
assert_eq!(system.bounded_variable_count(), 1);
|
||
}
|
||
|
||
#[test]
|
||
fn test_add_bounded_variable_duplicate_id() {
|
||
use crate::inverse::{BoundedVariable, BoundedVariableId};
|
||
|
||
let mut system = System::new();
|
||
|
||
let valve1 = BoundedVariable::new(BoundedVariableId::new("valve"), 0.5, 0.0, 1.0).unwrap();
|
||
let valve2 = BoundedVariable::new(BoundedVariableId::new("valve"), 0.3, 0.0, 1.0).unwrap();
|
||
|
||
assert!(system.add_bounded_variable(valve1).is_ok());
|
||
let result = system.add_bounded_variable(valve2);
|
||
assert!(matches!(
|
||
result,
|
||
Err(crate::inverse::BoundedVariableError::DuplicateId { .. })
|
||
));
|
||
assert_eq!(system.bounded_variable_count(), 1);
|
||
}
|
||
|
||
#[test]
|
||
fn test_add_bounded_variable_invalid_component() {
|
||
use crate::inverse::{BoundedVariable, BoundedVariableError, BoundedVariableId};
|
||
|
||
let mut system = System::new();
|
||
|
||
// Variable references non-existent component
|
||
let valve = BoundedVariable::with_component(
|
||
BoundedVariableId::new("valve"),
|
||
"unknown_component",
|
||
0.5,
|
||
0.0,
|
||
1.0,
|
||
)
|
||
.unwrap();
|
||
|
||
let result = system.add_bounded_variable(valve);
|
||
assert!(matches!(
|
||
result,
|
||
Err(BoundedVariableError::InvalidComponent { .. })
|
||
));
|
||
}
|
||
|
||
#[test]
|
||
fn test_add_bounded_variable_with_valid_component() {
|
||
use crate::inverse::{BoundedVariable, BoundedVariableId};
|
||
|
||
let mut system = System::new();
|
||
|
||
// Register component first
|
||
let node = system.add_component(make_mock(0));
|
||
system.register_component_name("expansion_valve", node);
|
||
|
||
let valve = BoundedVariable::with_component(
|
||
BoundedVariableId::new("valve"),
|
||
"expansion_valve",
|
||
0.5,
|
||
0.0,
|
||
1.0,
|
||
)
|
||
.unwrap();
|
||
|
||
assert!(system.add_bounded_variable(valve).is_ok());
|
||
}
|
||
|
||
#[test]
|
||
fn test_remove_bounded_variable() {
|
||
use crate::inverse::{BoundedVariable, BoundedVariableId};
|
||
|
||
let mut system = System::new();
|
||
|
||
let id = BoundedVariableId::new("valve");
|
||
let valve = BoundedVariable::new(id.clone(), 0.5, 0.0, 1.0).unwrap();
|
||
|
||
system.add_bounded_variable(valve).unwrap();
|
||
assert_eq!(system.bounded_variable_count(), 1);
|
||
|
||
let removed = system.remove_bounded_variable(&id);
|
||
assert!(removed.is_some());
|
||
assert_eq!(system.bounded_variable_count(), 0);
|
||
|
||
// Removing non-existent returns None
|
||
let removed_again = system.remove_bounded_variable(&id);
|
||
assert!(removed_again.is_none());
|
||
}
|
||
|
||
#[test]
|
||
fn test_get_bounded_variable() {
|
||
use crate::inverse::{BoundedVariable, BoundedVariableId};
|
||
|
||
let mut system = System::new();
|
||
|
||
let id = BoundedVariableId::new("vfd_speed");
|
||
let vfd = BoundedVariable::new(id.clone(), 0.8, 0.3, 1.0).unwrap();
|
||
|
||
system.add_bounded_variable(vfd).unwrap();
|
||
|
||
let retrieved = system.get_bounded_variable(&id);
|
||
assert!(retrieved.is_some());
|
||
approx::assert_relative_eq!(retrieved.unwrap().value(), 0.8);
|
||
|
||
let missing = system.get_bounded_variable(&BoundedVariableId::new("nonexistent"));
|
||
assert!(missing.is_none());
|
||
}
|
||
|
||
#[test]
|
||
fn test_get_bounded_variable_mut() {
|
||
use crate::inverse::{BoundedVariable, BoundedVariableId};
|
||
|
||
let mut system = System::new();
|
||
|
||
let id = BoundedVariableId::new("valve");
|
||
let valve = BoundedVariable::new(id.clone(), 0.5, 0.0, 1.0).unwrap();
|
||
|
||
system.add_bounded_variable(valve).unwrap();
|
||
|
||
// Apply step through mutable reference
|
||
if let Some(v) = system.get_bounded_variable_mut(&id) {
|
||
v.apply_step(0.3);
|
||
}
|
||
|
||
let retrieved = system.get_bounded_variable(&id).unwrap();
|
||
approx::assert_relative_eq!(retrieved.value(), 0.8);
|
||
}
|
||
|
||
#[test]
|
||
fn test_saturated_variables() {
|
||
use crate::inverse::{BoundedVariable, BoundedVariableId, SaturationType};
|
||
|
||
let mut system = System::new();
|
||
|
||
// Variable at min bound
|
||
let v1 = BoundedVariable::new(BoundedVariableId::new("v1"), 0.0, 0.0, 1.0).unwrap();
|
||
// Variable at max bound
|
||
let v2 = BoundedVariable::new(BoundedVariableId::new("v2"), 1.0, 0.0, 1.0).unwrap();
|
||
// Variable in middle (not saturated)
|
||
let v3 = BoundedVariable::new(BoundedVariableId::new("v3"), 0.5, 0.0, 1.0).unwrap();
|
||
|
||
system.add_bounded_variable(v1).unwrap();
|
||
system.add_bounded_variable(v2).unwrap();
|
||
system.add_bounded_variable(v3).unwrap();
|
||
|
||
let saturated = system.saturated_variables();
|
||
assert_eq!(saturated.len(), 2);
|
||
|
||
let sat_types: Vec<_> = saturated.iter().map(|s| s.saturation_type).collect();
|
||
assert!(sat_types.contains(&SaturationType::LowerBound));
|
||
assert!(sat_types.contains(&SaturationType::UpperBound));
|
||
}
|
||
|
||
#[test]
|
||
fn test_bounded_variables_iterator() {
|
||
use crate::inverse::{BoundedVariable, BoundedVariableId};
|
||
|
||
let mut system = System::new();
|
||
|
||
let v1 = BoundedVariable::new(BoundedVariableId::new("v1"), 0.5, 0.0, 1.0).unwrap();
|
||
let v2 = BoundedVariable::new(BoundedVariableId::new("v2"), 0.3, 0.0, 1.0).unwrap();
|
||
|
||
system.add_bounded_variable(v1).unwrap();
|
||
system.add_bounded_variable(v2).unwrap();
|
||
|
||
let count = system.bounded_variables().count();
|
||
assert_eq!(count, 2);
|
||
}
|
||
|
||
// ────────────────────────────────────────────────────────────────────────
|
||
// Inverse Control Tests (Story 5.3)
|
||
// ────────────────────────────────────────────────────────────────────────
|
||
|
||
#[test]
|
||
fn test_link_constraint_to_control() {
|
||
use crate::inverse::{
|
||
BoundedVariable, BoundedVariableId, ComponentOutput, Constraint, ConstraintId,
|
||
};
|
||
|
||
let mut system = System::new();
|
||
|
||
// Register component
|
||
let node = system.add_component(make_mock(0));
|
||
system.register_component_name("evaporator", node);
|
||
|
||
// Add constraint and bounded variable
|
||
let constraint = Constraint::new(
|
||
ConstraintId::new("superheat"),
|
||
ComponentOutput::Superheat {
|
||
component_id: "evaporator".to_string(),
|
||
},
|
||
5.0,
|
||
);
|
||
system.add_constraint(constraint).unwrap();
|
||
|
||
let valve = BoundedVariable::new(BoundedVariableId::new("valve"), 0.5, 0.0, 1.0).unwrap();
|
||
system.add_bounded_variable(valve).unwrap();
|
||
|
||
// Link them
|
||
let result = system.link_constraint_to_control(
|
||
&ConstraintId::new("superheat"),
|
||
&BoundedVariableId::new("valve"),
|
||
);
|
||
assert!(result.is_ok());
|
||
assert_eq!(system.inverse_control_mapping_count(), 1);
|
||
|
||
// Verify bidirectional lookup
|
||
let control = system.get_control_for_constraint(&ConstraintId::new("superheat"));
|
||
assert!(control.is_some());
|
||
assert_eq!(control.unwrap().as_str(), "valve");
|
||
|
||
let constraint = system.get_constraint_for_control(&BoundedVariableId::new("valve"));
|
||
assert!(constraint.is_some());
|
||
assert_eq!(constraint.unwrap().as_str(), "superheat");
|
||
}
|
||
|
||
#[test]
|
||
fn test_link_constraint_not_found() {
|
||
use crate::inverse::{BoundedVariable, BoundedVariableId, ConstraintId, DoFError};
|
||
|
||
let mut system = System::new();
|
||
|
||
let valve = BoundedVariable::new(BoundedVariableId::new("valve"), 0.5, 0.0, 1.0).unwrap();
|
||
system.add_bounded_variable(valve).unwrap();
|
||
|
||
let result = system.link_constraint_to_control(
|
||
&ConstraintId::new("nonexistent"),
|
||
&BoundedVariableId::new("valve"),
|
||
);
|
||
assert!(matches!(result, Err(DoFError::ConstraintNotFound { .. })));
|
||
}
|
||
|
||
#[test]
|
||
fn test_link_control_not_found() {
|
||
use crate::inverse::{
|
||
BoundedVariableId, ComponentOutput, Constraint, ConstraintId, DoFError,
|
||
};
|
||
|
||
let mut system = System::new();
|
||
|
||
let node = system.add_component(make_mock(0));
|
||
system.register_component_name("evaporator", node);
|
||
|
||
let constraint = Constraint::new(
|
||
ConstraintId::new("superheat"),
|
||
ComponentOutput::Superheat {
|
||
component_id: "evaporator".to_string(),
|
||
},
|
||
5.0,
|
||
);
|
||
system.add_constraint(constraint).unwrap();
|
||
|
||
let result = system.link_constraint_to_control(
|
||
&ConstraintId::new("superheat"),
|
||
&BoundedVariableId::new("nonexistent"),
|
||
);
|
||
assert!(matches!(
|
||
result,
|
||
Err(DoFError::BoundedVariableNotFound { .. })
|
||
));
|
||
}
|
||
|
||
#[test]
|
||
fn test_link_duplicate_constraint() {
|
||
use crate::inverse::{
|
||
BoundedVariable, BoundedVariableId, ComponentOutput, Constraint, ConstraintId, DoFError,
|
||
};
|
||
|
||
let mut system = System::new();
|
||
|
||
let node = system.add_component(make_mock(0));
|
||
system.register_component_name("evaporator", node);
|
||
|
||
let constraint = Constraint::new(
|
||
ConstraintId::new("superheat"),
|
||
ComponentOutput::Superheat {
|
||
component_id: "evaporator".to_string(),
|
||
},
|
||
5.0,
|
||
);
|
||
system.add_constraint(constraint).unwrap();
|
||
|
||
let v1 = BoundedVariable::new(BoundedVariableId::new("v1"), 0.5, 0.0, 1.0).unwrap();
|
||
let v2 = BoundedVariable::new(BoundedVariableId::new("v2"), 0.5, 0.0, 1.0).unwrap();
|
||
system.add_bounded_variable(v1).unwrap();
|
||
system.add_bounded_variable(v2).unwrap();
|
||
|
||
// First link should succeed
|
||
system
|
||
.link_constraint_to_control(
|
||
&ConstraintId::new("superheat"),
|
||
&BoundedVariableId::new("v1"),
|
||
)
|
||
.unwrap();
|
||
|
||
// Second link for same constraint should fail
|
||
let result = system.link_constraint_to_control(
|
||
&ConstraintId::new("superheat"),
|
||
&BoundedVariableId::new("v2"),
|
||
);
|
||
assert!(matches!(result, Err(DoFError::AlreadyLinked { .. })));
|
||
}
|
||
|
||
#[test]
|
||
fn test_link_duplicate_control() {
|
||
use crate::inverse::{
|
||
BoundedVariable, BoundedVariableId, ComponentOutput, Constraint, ConstraintId, DoFError,
|
||
};
|
||
|
||
let mut system = System::new();
|
||
|
||
let node = system.add_component(make_mock(0));
|
||
system.register_component_name("evaporator", node);
|
||
|
||
let c1 = Constraint::new(
|
||
ConstraintId::new("c1"),
|
||
ComponentOutput::Superheat {
|
||
component_id: "evaporator".to_string(),
|
||
},
|
||
5.0,
|
||
);
|
||
let c2 = Constraint::new(
|
||
ConstraintId::new("c2"),
|
||
ComponentOutput::Temperature {
|
||
component_id: "evaporator".to_string(),
|
||
},
|
||
300.0,
|
||
);
|
||
system.add_constraint(c1).unwrap();
|
||
system.add_constraint(c2).unwrap();
|
||
|
||
let valve = BoundedVariable::new(BoundedVariableId::new("valve"), 0.5, 0.0, 1.0).unwrap();
|
||
system.add_bounded_variable(valve).unwrap();
|
||
|
||
// First link should succeed
|
||
system
|
||
.link_constraint_to_control(&ConstraintId::new("c1"), &BoundedVariableId::new("valve"))
|
||
.unwrap();
|
||
|
||
// Second link for same control should fail
|
||
let result = system
|
||
.link_constraint_to_control(&ConstraintId::new("c2"), &BoundedVariableId::new("valve"));
|
||
assert!(matches!(result, Err(DoFError::ControlAlreadyLinked { .. })));
|
||
}
|
||
|
||
#[test]
|
||
fn test_unlink_constraint() {
|
||
use crate::inverse::{
|
||
BoundedVariable, BoundedVariableId, ComponentOutput, Constraint, ConstraintId,
|
||
};
|
||
|
||
let mut system = System::new();
|
||
|
||
let node = system.add_component(make_mock(0));
|
||
system.register_component_name("evaporator", node);
|
||
|
||
let constraint = Constraint::new(
|
||
ConstraintId::new("superheat"),
|
||
ComponentOutput::Superheat {
|
||
component_id: "evaporator".to_string(),
|
||
},
|
||
5.0,
|
||
);
|
||
system.add_constraint(constraint).unwrap();
|
||
|
||
let valve = BoundedVariable::new(BoundedVariableId::new("valve"), 0.5, 0.0, 1.0).unwrap();
|
||
system.add_bounded_variable(valve).unwrap();
|
||
|
||
system
|
||
.link_constraint_to_control(
|
||
&ConstraintId::new("superheat"),
|
||
&BoundedVariableId::new("valve"),
|
||
)
|
||
.unwrap();
|
||
assert_eq!(system.inverse_control_mapping_count(), 1);
|
||
|
||
// Unlink
|
||
let removed = system.unlink_constraint(&ConstraintId::new("superheat"));
|
||
assert!(removed.is_some());
|
||
assert_eq!(removed.unwrap().as_str(), "valve");
|
||
assert_eq!(system.inverse_control_mapping_count(), 0);
|
||
|
||
// Unlinking again returns None
|
||
let removed_again = system.unlink_constraint(&ConstraintId::new("superheat"));
|
||
assert!(removed_again.is_none());
|
||
}
|
||
|
||
#[test]
|
||
fn test_control_variable_state_index() {
|
||
use crate::inverse::{
|
||
BoundedVariable, BoundedVariableId, ComponentOutput, Constraint, ConstraintId,
|
||
};
|
||
|
||
let mut system = System::new();
|
||
|
||
// Add two components and an edge
|
||
let n0 = system.add_component(make_mock(0));
|
||
let n1 = system.add_component(make_mock(0));
|
||
system.register_component_name("evaporator", n0);
|
||
system.register_component_name("condenser", n1);
|
||
system.add_edge(n0, n1).unwrap();
|
||
system.finalize().unwrap();
|
||
|
||
// edge_count = 1, so base index = 3 (3 unknowns per edge: ṁ, P, h)
|
||
assert_eq!(system.edge_count(), 1);
|
||
|
||
// Add constraint and bounded variable
|
||
let constraint = Constraint::new(
|
||
ConstraintId::new("superheat"),
|
||
ComponentOutput::Superheat {
|
||
component_id: "evaporator".to_string(),
|
||
},
|
||
5.0,
|
||
);
|
||
system.add_constraint(constraint).unwrap();
|
||
|
||
let valve = BoundedVariable::new(BoundedVariableId::new("valve"), 0.5, 0.0, 1.0).unwrap();
|
||
system.add_bounded_variable(valve).unwrap();
|
||
|
||
// Link them
|
||
system
|
||
.link_constraint_to_control(
|
||
&ConstraintId::new("superheat"),
|
||
&BoundedVariableId::new("valve"),
|
||
)
|
||
.unwrap();
|
||
|
||
// Control variable index should be at total_state_len = 3 * edge_count = 3
|
||
let idx = system.control_variable_state_index(&BoundedVariableId::new("valve"));
|
||
assert!(idx.is_some());
|
||
assert_eq!(idx.unwrap(), 3);
|
||
|
||
// Unlinked variable returns None
|
||
let v2 = BoundedVariable::new(BoundedVariableId::new("v2"), 0.5, 0.0, 1.0).unwrap();
|
||
system.add_bounded_variable(v2).unwrap();
|
||
let idx2 = system.control_variable_state_index(&BoundedVariableId::new("v2"));
|
||
assert!(idx2.is_none());
|
||
}
|
||
|
||
#[test]
|
||
fn test_validate_inverse_control_dof_balanced() {
|
||
use crate::inverse::{
|
||
BoundedVariable, BoundedVariableId, ComponentOutput, Constraint, ConstraintId,
|
||
};
|
||
|
||
let mut system = System::new();
|
||
|
||
// Add two components: one with 2 equations, one with 1 equation (total 3)
|
||
// 1 edge (stride-3) = 3 unknowns → balanced without closures (CM1.3+)
|
||
let n0 = system.add_component(make_mock(2));
|
||
let n1 = system.add_component(make_mock(1));
|
||
system.register_component_name("evaporator", n0);
|
||
system.add_edge(n0, n1).unwrap();
|
||
system.finalize().unwrap();
|
||
|
||
// n_component_eqs = 3 (2 + 1)
|
||
// n_edge_unknowns = 3 (1 edge × stride-3)
|
||
|
||
// Add one constraint and one control
|
||
let constraint = Constraint::new(
|
||
ConstraintId::new("superheat"),
|
||
ComponentOutput::Superheat {
|
||
component_id: "evaporator".to_string(),
|
||
},
|
||
5.0,
|
||
);
|
||
system.add_constraint(constraint).unwrap();
|
||
|
||
let valve = BoundedVariable::new(BoundedVariableId::new("valve"), 0.5, 0.0, 1.0).unwrap();
|
||
system.add_bounded_variable(valve).unwrap();
|
||
|
||
system
|
||
.link_constraint_to_control(
|
||
&ConstraintId::new("superheat"),
|
||
&BoundedVariableId::new("valve"),
|
||
)
|
||
.unwrap();
|
||
|
||
// n_equations = 3 (components) + 1 (constraint) = 4
|
||
// n_unknowns = 3 (edge) + 1 (control) = 4 → Balanced!
|
||
let result = system.validate_inverse_control_dof();
|
||
assert!(result.is_ok());
|
||
}
|
||
|
||
#[test]
|
||
fn test_validate_inverse_control_dof_over_constrained() {
|
||
use crate::inverse::{
|
||
BoundedVariable, BoundedVariableId, ComponentOutput, Constraint, ConstraintId, DoFError,
|
||
};
|
||
|
||
let mut system = System::new();
|
||
|
||
// Add two components: 2 + 1 equations = 3 total; 1 edge = 3 unknowns
|
||
let n0 = system.add_component(make_mock(2));
|
||
let n1 = system.add_component(make_mock(1));
|
||
system.register_component_name("evaporator", n0);
|
||
system.register_component_name("condenser", n1);
|
||
system.add_edge(n0, n1).unwrap();
|
||
system.finalize().unwrap();
|
||
|
||
// Add two constraints but only one control
|
||
let c1 = Constraint::new(
|
||
ConstraintId::new("c1"),
|
||
ComponentOutput::Superheat {
|
||
component_id: "evaporator".to_string(),
|
||
},
|
||
5.0,
|
||
);
|
||
let c2 = Constraint::new(
|
||
ConstraintId::new("c2"),
|
||
ComponentOutput::Temperature {
|
||
component_id: "condenser".to_string(),
|
||
},
|
||
300.0,
|
||
);
|
||
system.add_constraint(c1).unwrap();
|
||
system.add_constraint(c2).unwrap();
|
||
|
||
let valve = BoundedVariable::new(BoundedVariableId::new("valve"), 0.5, 0.0, 1.0).unwrap();
|
||
system.add_bounded_variable(valve).unwrap();
|
||
|
||
system
|
||
.link_constraint_to_control(&ConstraintId::new("c1"), &BoundedVariableId::new("valve"))
|
||
.unwrap();
|
||
|
||
// n_equations = 3 (components) + 2 (constraints) = 5
|
||
// n_unknowns = 3 (1 edge × ṁ,P,h) + 1 (control) = 4
|
||
// Over-constrained!
|
||
let result = system.validate_inverse_control_dof();
|
||
assert!(matches!(
|
||
result,
|
||
Err(DoFError::OverConstrainedSystem { .. })
|
||
));
|
||
if let Err(DoFError::OverConstrainedSystem {
|
||
constraint_count,
|
||
control_count,
|
||
equation_count,
|
||
unknown_count,
|
||
}) = result
|
||
{
|
||
assert_eq!(constraint_count, 2);
|
||
assert_eq!(control_count, 1);
|
||
assert_eq!(equation_count, 5);
|
||
assert_eq!(unknown_count, 4);
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
fn test_validate_inverse_control_dof_under_constrained() {
|
||
use crate::inverse::{
|
||
BoundedVariable, BoundedVariableId, ComponentOutput, Constraint, ConstraintId,
|
||
};
|
||
|
||
let mut system = System::new();
|
||
|
||
// Add two components: 2 + 1 equations = 3 total; 1 edge = 3 unknowns
|
||
let n0 = system.add_component(make_mock(2));
|
||
let n1 = system.add_component(make_mock(1));
|
||
system.register_component_name("evaporator", n0);
|
||
system.add_edge(n0, n1).unwrap();
|
||
system.finalize().unwrap();
|
||
|
||
// Add one constraint and one control
|
||
let constraint = Constraint::new(
|
||
ConstraintId::new("superheat"),
|
||
ComponentOutput::Superheat {
|
||
component_id: "evaporator".to_string(),
|
||
},
|
||
5.0,
|
||
);
|
||
system.add_constraint(constraint).unwrap();
|
||
|
||
// Add two bounded variables but only link one
|
||
let v1 = BoundedVariable::new(BoundedVariableId::new("v1"), 0.5, 0.0, 1.0).unwrap();
|
||
let v2 = BoundedVariable::new(BoundedVariableId::new("v2"), 0.5, 0.0, 1.0).unwrap();
|
||
system.add_bounded_variable(v1).unwrap();
|
||
system.add_bounded_variable(v2).unwrap();
|
||
|
||
system
|
||
.link_constraint_to_control(
|
||
&ConstraintId::new("superheat"),
|
||
&BoundedVariableId::new("v1"),
|
||
)
|
||
.unwrap();
|
||
|
||
// n_equations = 3 (components) + 1 (constraint) = 4
|
||
// n_unknowns = 3 (edge) + 1 (linked control) = 4
|
||
// Balanced — unlinked bounded variables don't affect DoF
|
||
let result = system.validate_inverse_control_dof();
|
||
assert!(result.is_ok());
|
||
}
|
||
|
||
#[test]
|
||
fn test_full_state_vector_len() {
|
||
use crate::inverse::{
|
||
BoundedVariable, BoundedVariableId, ComponentOutput, Constraint, ConstraintId,
|
||
};
|
||
|
||
let mut system = System::new();
|
||
|
||
// Add two components and one edge
|
||
let n0 = system.add_component(make_mock(0));
|
||
let n1 = system.add_component(make_mock(0));
|
||
system.register_component_name("evaporator", n0);
|
||
system.add_edge(n0, n1).unwrap();
|
||
system.finalize().unwrap();
|
||
|
||
// Edge states: 3 * 1 = 3 (ṁ, P, h)
|
||
assert_eq!(system.full_state_vector_len(), 3);
|
||
|
||
// Add constraint and control
|
||
let constraint = Constraint::new(
|
||
ConstraintId::new("superheat"),
|
||
ComponentOutput::Superheat {
|
||
component_id: "evaporator".to_string(),
|
||
},
|
||
5.0,
|
||
);
|
||
system.add_constraint(constraint).unwrap();
|
||
|
||
let valve = BoundedVariable::new(BoundedVariableId::new("valve"), 0.5, 0.0, 1.0).unwrap();
|
||
system.add_bounded_variable(valve).unwrap();
|
||
|
||
system
|
||
.link_constraint_to_control(
|
||
&ConstraintId::new("superheat"),
|
||
&BoundedVariableId::new("valve"),
|
||
)
|
||
.unwrap();
|
||
|
||
// Edge states: 3, control vars: 1, thermal couplings: 0
|
||
// Total: 3 + 1 + 0 = 4
|
||
assert_eq!(system.full_state_vector_len(), 4);
|
||
}
|
||
|
||
#[test]
|
||
fn test_free_actuator_extends_state_and_indexes_after_saturated() {
|
||
use crate::inverse::{BoundedVariable, BoundedVariableId};
|
||
|
||
let mut system = System::new();
|
||
let n0 = system.add_component(make_mock(0));
|
||
let n1 = system.add_component(make_mock(0));
|
||
system.register_component_name("exv", n0);
|
||
system.add_edge(n0, n1).unwrap();
|
||
|
||
// Edge states: 3 * 1 = 3 (ṁ, P, h). No control yet.
|
||
// (length is only meaningful after finalize but the actuator arithmetic
|
||
// does not depend on it.)
|
||
let open = BoundedVariable::with_component(
|
||
BoundedVariableId::new("exv__opening"),
|
||
"exv",
|
||
0.5,
|
||
0.02,
|
||
1.0,
|
||
)
|
||
.unwrap();
|
||
system.add_bounded_variable(open).unwrap();
|
||
system.add_free_actuator(BoundedVariableId::new("exv__opening"));
|
||
|
||
system.finalize().unwrap();
|
||
|
||
assert_eq!(system.free_actuator_count(), 1);
|
||
// 3 edge unknowns + 1 free actuator = 4.
|
||
assert_eq!(system.full_state_vector_len(), 4);
|
||
// The free actuator lives right after the (empty) saturated block, i.e.
|
||
// at total_state_len + 0 + 0 + 0 = 3.
|
||
assert_eq!(system.free_actuator_index(0), 3);
|
||
}
|
||
|
||
#[test]
|
||
#[should_panic(expected = "already registered")]
|
||
fn test_free_actuator_rejects_duplicate() {
|
||
use crate::inverse::{BoundedVariable, BoundedVariableId};
|
||
|
||
let mut system = System::new();
|
||
let n0 = system.add_component(make_mock(0));
|
||
system.register_component_name("exv", n0);
|
||
let open = BoundedVariable::with_component(
|
||
BoundedVariableId::new("exv__opening"),
|
||
"exv",
|
||
0.5,
|
||
0.0,
|
||
1.0,
|
||
)
|
||
.unwrap();
|
||
system.add_bounded_variable(open).unwrap();
|
||
system.add_free_actuator(BoundedVariableId::new("exv__opening"));
|
||
system.add_free_actuator(BoundedVariableId::new("exv__opening"));
|
||
}
|
||
|
||
#[test]
|
||
fn test_saturated_controller_residuals_and_jacobian() {
|
||
use crate::inverse::{
|
||
BoundedVariable, BoundedVariableId, ComponentOutput, ConstraintId, SaturatedController,
|
||
Saturation,
|
||
};
|
||
|
||
let mut system = System::new();
|
||
let plant =
|
||
system.add_component(Box::new(SaturatedControlMockComponent { z_flow_idx: None }));
|
||
system.register_component_name("plant", plant);
|
||
system.add_edge(plant, plant).unwrap();
|
||
|
||
let actuator = BoundedVariable::with_component(
|
||
BoundedVariableId::new("plant_f_m"),
|
||
"plant",
|
||
1.0,
|
||
0.5,
|
||
1.5,
|
||
)
|
||
.unwrap();
|
||
system.add_bounded_variable(actuator).unwrap();
|
||
let ctrl = SaturatedController::new(
|
||
ConstraintId::new("capacity_loop"),
|
||
ComponentOutput::Capacity {
|
||
component_id: "plant".to_string(),
|
||
},
|
||
BoundedVariableId::new("plant_f_m"),
|
||
15.0,
|
||
0.5,
|
||
1.5,
|
||
)
|
||
.unwrap()
|
||
.with_gain(0.5)
|
||
.unwrap()
|
||
.with_band(1.0)
|
||
.unwrap()
|
||
.with_saturation(Saturation::Hard);
|
||
system.add_saturated_controller(ctrl.clone());
|
||
system.finalize().unwrap();
|
||
|
||
assert_eq!(system.saturated_controller_count(), 1);
|
||
assert_eq!(system.full_state_vector_len(), 5);
|
||
assert_eq!(
|
||
system.control_variable_state_index(&BoundedVariableId::new("plant_f_m")),
|
||
Some(3)
|
||
);
|
||
|
||
let state = vec![1.0, 2.0, 3.0, 1.1, 0.2];
|
||
let mut residuals = vec![0.0; 3];
|
||
system.compute_residuals(&state, &mut residuals).unwrap();
|
||
assert!((residuals[1] - ctrl.residual_u(state[3], state[4])).abs() < 1e-12);
|
||
assert!((residuals[2] - ctrl.residual_y(14.0, state[4])).abs() < 1e-12);
|
||
|
||
let mut jacobian = JacobianBuilder::new();
|
||
system.assemble_jacobian(&state, &mut jacobian).unwrap();
|
||
let entries = jacobian.entries();
|
||
assert!(entries.contains(&(1, 3, 1.0)));
|
||
assert!(entries.contains(&(1, 4, -0.5)));
|
||
// Offset-free control law: ∂r_y/∂x = -(1 - S'(x)); x=0.2 is in-band so
|
||
// S'(x)=1 and the entry is 0 (the row is closed via the plant coupling).
|
||
assert!(entries.contains(&(2, 4, 0.0)));
|
||
|
||
let mut analytic = vec![vec![0.0; state.len()]; residuals.len()];
|
||
for &(row, col, value) in entries {
|
||
analytic[row][col] += value;
|
||
}
|
||
|
||
let eps = 1e-6;
|
||
let mut max_err = 0.0_f64;
|
||
for col in 0..state.len() {
|
||
let mut plus_state = state.clone();
|
||
let mut minus_state = state.clone();
|
||
plus_state[col] += eps;
|
||
minus_state[col] -= eps;
|
||
|
||
let mut plus = vec![0.0; residuals.len()];
|
||
let mut minus = vec![0.0; residuals.len()];
|
||
system.compute_residuals(&plus_state, &mut plus).unwrap();
|
||
system.compute_residuals(&minus_state, &mut minus).unwrap();
|
||
|
||
for row in 0..residuals.len() {
|
||
let fd = (plus[row] - minus[row]) / (2.0 * eps);
|
||
max_err = max_err.max((analytic[row][col] - fd).abs());
|
||
}
|
||
}
|
||
println!("saturated controller jacobian max error: {max_err:.3e}");
|
||
assert!(
|
||
max_err <= 1e-4,
|
||
"saturated controller analytic Jacobian mismatch: max_err={max_err}"
|
||
);
|
||
}
|
||
|
||
#[test]
|
||
fn test_saturated_network_controller_jacobian_matches_fd() {
|
||
use crate::inverse::{
|
||
BoundedVariable, BoundedVariableId, Combine, ComponentOutput, ConstraintId, Objective,
|
||
SaturatedController, Saturation,
|
||
};
|
||
|
||
let mut system = System::new();
|
||
let plant =
|
||
system.add_component(Box::new(SaturatedControlMockComponent { z_flow_idx: None }));
|
||
system.register_component_name("plant", plant);
|
||
system.add_edge(plant, plant).unwrap();
|
||
|
||
let actuator = BoundedVariable::with_component(
|
||
BoundedVariableId::new("plant_f_m"),
|
||
"plant",
|
||
1.0,
|
||
0.5,
|
||
1.5,
|
||
)
|
||
.unwrap();
|
||
system.add_bounded_variable(actuator).unwrap();
|
||
|
||
// Override network: primary Capacity→15 seeds the fold; a Temperature→2
|
||
// protection is folded in via softMax. A moderate alpha keeps the
|
||
// selector well inside its smooth regime at the test point.
|
||
let ctrl = SaturatedController::new(
|
||
ConstraintId::new("net_loop"),
|
||
ComponentOutput::Capacity {
|
||
component_id: "plant".to_string(),
|
||
},
|
||
BoundedVariableId::new("plant_f_m"),
|
||
15.0,
|
||
0.5,
|
||
1.5,
|
||
)
|
||
.unwrap()
|
||
.with_gain(0.5)
|
||
.unwrap()
|
||
.with_band(1.0)
|
||
.unwrap()
|
||
.with_saturation(Saturation::Hard)
|
||
.with_alpha(0.5)
|
||
.with_objectives(vec![
|
||
Objective::new(
|
||
ComponentOutput::Capacity {
|
||
component_id: "plant".to_string(),
|
||
},
|
||
15.0,
|
||
0.5,
|
||
Combine::Min,
|
||
),
|
||
Objective::new(
|
||
ComponentOutput::Temperature {
|
||
component_id: "plant".to_string(),
|
||
},
|
||
2.0,
|
||
-0.7,
|
||
Combine::Max,
|
||
),
|
||
]);
|
||
system.add_saturated_controller(ctrl.clone());
|
||
system.finalize().unwrap();
|
||
|
||
assert!(ctrl.is_network());
|
||
|
||
let state = vec![1.0, 2.0, 3.0, 1.1, 0.2];
|
||
let mut residuals = vec![0.0; 3];
|
||
system.compute_residuals(&state, &mut residuals).unwrap();
|
||
|
||
// Cross-check the network residual against a direct evaluation.
|
||
let cap = state[0] + 2.0 * state[1] + 3.0 * state[2];
|
||
let temp = state[1];
|
||
let e = ctrl.error_signal(&[cap, temp]);
|
||
assert!((residuals[2] - ctrl.residual_y_network(e, state[4])).abs() < 1e-12);
|
||
|
||
let mut jacobian = JacobianBuilder::new();
|
||
system.assemble_jacobian(&state, &mut jacobian).unwrap();
|
||
let entries = jacobian.entries();
|
||
let mut analytic = vec![vec![0.0; state.len()]; residuals.len()];
|
||
for &(row, col, value) in entries {
|
||
analytic[row][col] += value;
|
||
}
|
||
|
||
let eps = 1e-6;
|
||
let mut max_err = 0.0_f64;
|
||
for col in 0..state.len() {
|
||
let mut plus_state = state.clone();
|
||
let mut minus_state = state.clone();
|
||
plus_state[col] += eps;
|
||
minus_state[col] -= eps;
|
||
let mut plus = vec![0.0; residuals.len()];
|
||
let mut minus = vec![0.0; residuals.len()];
|
||
system.compute_residuals(&plus_state, &mut plus).unwrap();
|
||
system.compute_residuals(&minus_state, &mut minus).unwrap();
|
||
for row in 0..residuals.len() {
|
||
let fd = (plus[row] - minus[row]) / (2.0 * eps);
|
||
max_err = max_err.max((analytic[row][col] - fd).abs());
|
||
}
|
||
}
|
||
assert!(
|
||
max_err <= 1e-4,
|
||
"network controller analytic Jacobian mismatch: max_err={max_err}"
|
||
);
|
||
}
|
||
|
||
#[test]
|
||
fn test_control_variable_indices() {
|
||
use crate::inverse::{
|
||
BoundedVariable, BoundedVariableId, ComponentOutput, Constraint, ConstraintId,
|
||
};
|
||
|
||
let mut system = System::new();
|
||
|
||
// Add two components and one edge
|
||
let n0 = system.add_component(make_mock(0));
|
||
let n1 = system.add_component(make_mock(0));
|
||
system.register_component_name("evaporator", n0);
|
||
system.register_component_name("condenser", n1);
|
||
system.add_edge(n0, n1).unwrap();
|
||
system.finalize().unwrap();
|
||
|
||
// Add two constraints and controls
|
||
let c1 = Constraint::new(
|
||
ConstraintId::new("c1"),
|
||
ComponentOutput::Superheat {
|
||
component_id: "evaporator".to_string(),
|
||
},
|
||
5.0,
|
||
);
|
||
system.add_constraint(c1).unwrap();
|
||
|
||
let v1 = BoundedVariable::new(BoundedVariableId::new("v1"), 0.5, 0.0, 1.0).unwrap();
|
||
system.add_bounded_variable(v1).unwrap();
|
||
|
||
system
|
||
.link_constraint_to_control(&ConstraintId::new("c1"), &BoundedVariableId::new("v1"))
|
||
.unwrap();
|
||
|
||
let indices = system.control_variable_indices();
|
||
assert_eq!(indices.len(), 1);
|
||
assert_eq!(indices[0].1, 3); // total_state_len = 3 * edge_count = 3
|
||
}
|
||
|
||
struct BadMassFlowComponent {
|
||
ports: Vec<ConnectedPort>,
|
||
}
|
||
|
||
impl Component for BadMassFlowComponent {
|
||
fn compute_residuals(
|
||
&self,
|
||
_state: &StateSlice,
|
||
_residuals: &mut entropyk_components::ResidualVector,
|
||
) -> Result<(), ComponentError> {
|
||
Ok(())
|
||
}
|
||
|
||
fn jacobian_entries(
|
||
&self,
|
||
_state: &StateSlice,
|
||
_jacobian: &mut JacobianBuilder,
|
||
) -> Result<(), ComponentError> {
|
||
Ok(())
|
||
}
|
||
|
||
fn n_equations(&self) -> usize {
|
||
0
|
||
}
|
||
|
||
fn get_ports(&self) -> &[ConnectedPort] {
|
||
&self.ports
|
||
}
|
||
|
||
fn port_mass_flows(
|
||
&self,
|
||
_state: &StateSlice,
|
||
) -> Result<Vec<entropyk_core::MassFlow>, ComponentError> {
|
||
Ok(vec![
|
||
entropyk_core::MassFlow::from_kg_per_s(1.0),
|
||
entropyk_core::MassFlow::from_kg_per_s(-0.5), // Intentionally unbalanced
|
||
])
|
||
}
|
||
}
|
||
|
||
/// Component with balanced mass flow (inlet = outlet)
|
||
struct BalancedMassFlowComponent {
|
||
ports: Vec<ConnectedPort>,
|
||
}
|
||
|
||
impl Component for BalancedMassFlowComponent {
|
||
fn compute_residuals(
|
||
&self,
|
||
_state: &StateSlice,
|
||
_residuals: &mut entropyk_components::ResidualVector,
|
||
) -> Result<(), ComponentError> {
|
||
Ok(())
|
||
}
|
||
|
||
fn jacobian_entries(
|
||
&self,
|
||
_state: &StateSlice,
|
||
_jacobian: &mut JacobianBuilder,
|
||
) -> Result<(), ComponentError> {
|
||
Ok(())
|
||
}
|
||
|
||
fn n_equations(&self) -> usize {
|
||
0
|
||
}
|
||
|
||
fn get_ports(&self) -> &[ConnectedPort] {
|
||
&self.ports
|
||
}
|
||
|
||
fn port_mass_flows(
|
||
&self,
|
||
_state: &StateSlice,
|
||
) -> Result<Vec<entropyk_core::MassFlow>, ComponentError> {
|
||
// Balanced: inlet = 1.0 kg/s, outlet = -1.0 kg/s (sum = 0)
|
||
Ok(vec![
|
||
entropyk_core::MassFlow::from_kg_per_s(1.0),
|
||
entropyk_core::MassFlow::from_kg_per_s(-1.0),
|
||
])
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
fn test_mass_balance_passes_for_balanced_component() {
|
||
let mut system = System::new();
|
||
|
||
let inlet = Port::new(
|
||
FluidId::new("R134a"),
|
||
Pressure::from_bar(1.0),
|
||
Enthalpy::from_joules_per_kg(400000.0),
|
||
);
|
||
let outlet = Port::new(
|
||
FluidId::new("R134a"),
|
||
Pressure::from_bar(1.0),
|
||
Enthalpy::from_joules_per_kg(400000.0),
|
||
);
|
||
let (c1, c2) = inlet.connect(outlet).unwrap();
|
||
|
||
let comp = Box::new(BalancedMassFlowComponent {
|
||
ports: vec![c1, c2],
|
||
});
|
||
|
||
let n0 = system.add_component(comp);
|
||
system.add_edge(n0, n0).unwrap(); // Self-edge to avoid isolated node
|
||
|
||
system.finalize().unwrap();
|
||
|
||
let state = vec![0.0; system.full_state_vector_len()];
|
||
let result = system.check_mass_balance(&state);
|
||
|
||
assert!(
|
||
result.is_ok(),
|
||
"Expected mass balance to pass for balanced component"
|
||
);
|
||
}
|
||
|
||
#[test]
|
||
fn test_mass_balance_violation() {
|
||
let mut system = System::new();
|
||
|
||
let inlet = Port::new(
|
||
FluidId::new("R134a"),
|
||
Pressure::from_bar(1.0),
|
||
Enthalpy::from_joules_per_kg(400000.0),
|
||
);
|
||
let outlet = Port::new(
|
||
FluidId::new("R134a"),
|
||
Pressure::from_bar(1.0),
|
||
Enthalpy::from_joules_per_kg(400000.0),
|
||
);
|
||
let (c1, c2) = inlet.connect(outlet).unwrap();
|
||
|
||
let comp = Box::new(BadMassFlowComponent {
|
||
ports: vec![c1, c2], // Just to have ports
|
||
});
|
||
|
||
let n0 = system.add_component(comp);
|
||
system.add_edge(n0, n0).unwrap(); // Self-edge to avoid isolated node
|
||
|
||
system.finalize().unwrap();
|
||
|
||
// Ensure state is appropriately sized for finalize
|
||
let state = vec![0.0; system.full_state_vector_len()];
|
||
let result = system.check_mass_balance(&state);
|
||
|
||
assert!(result.is_err());
|
||
|
||
// Verify error contains mass error information
|
||
if let Err(crate::SolverError::Validation {
|
||
mass_error,
|
||
energy_error,
|
||
}) = result
|
||
{
|
||
assert!(mass_error > 0.0, "Mass error should be positive");
|
||
assert_eq!(
|
||
energy_error, 0.0,
|
||
"Energy error should be zero for mass-only validation"
|
||
);
|
||
} else {
|
||
panic!("Expected Validation error, got {:?}", result);
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
fn test_mass_balance_tolerance_constant() {
|
||
// Verify the tolerance constant is accessible and has expected value
|
||
assert_eq!(System::MASS_BALANCE_TOLERANCE_KG_S, 1e-9);
|
||
}
|
||
|
||
#[test]
|
||
fn test_generate_canonical_bytes() {
|
||
let mut sys = System::new();
|
||
let n0 = sys.add_component(make_mock(0));
|
||
let n1 = sys.add_component(make_mock(0));
|
||
sys.add_edge(n0, n1).unwrap();
|
||
|
||
let bytes1 = sys.generate_canonical_bytes();
|
||
let bytes2 = sys.generate_canonical_bytes();
|
||
|
||
// Exact same graph state should produce same bytes
|
||
assert_eq!(bytes1, bytes2);
|
||
}
|
||
|
||
#[test]
|
||
fn test_input_hash_deterministic() {
|
||
let mut sys1 = System::new();
|
||
let n0_1 = sys1.add_component(make_mock(0));
|
||
let n1_1 = sys1.add_component(make_mock(0));
|
||
sys1.add_edge(n0_1, n1_1).unwrap();
|
||
|
||
let mut sys2 = System::new();
|
||
let n0_2 = sys2.add_component(make_mock(0));
|
||
let n1_2 = sys2.add_component(make_mock(0));
|
||
sys2.add_edge(n0_2, n1_2).unwrap();
|
||
|
||
// Two identically constructed systems should have same hash
|
||
assert_eq!(sys1.input_hash(), sys2.input_hash());
|
||
|
||
// Now mutate one system by adding an edge
|
||
sys1.add_edge(n1_1, n0_1).unwrap();
|
||
|
||
// Hash should be different now
|
||
assert_ne!(sys1.input_hash(), sys2.input_hash());
|
||
}
|
||
|
||
// ────────────────────────────────────────────────────────────────────────
|
||
// Story 9.6: Energy Validation Logging Improvement Tests
|
||
// ────────────────────────────────────────────────────────────────────────
|
||
// Story 9.6: Energy Validation Logging Improvement Tests
|
||
// ────────────────────────────────────────────────────────────────────────
|
||
|
||
/// Test that check_energy_balance emits warnings for components without energy methods.
|
||
/// This test verifies the logging improvement from Story 9.6.
|
||
#[test]
|
||
fn test_energy_balance_warns_for_skipped_components() {
|
||
use tracing_subscriber::layer::SubscriberExt;
|
||
use tracing_subscriber::util::SubscriberInitExt;
|
||
|
||
// Create a system with mock components that don't implement energy_transfers()
|
||
let mut sys = System::new();
|
||
let n0 = sys.add_component(make_mock(0));
|
||
let n1 = sys.add_component(make_mock(0));
|
||
sys.add_edge(n0, n1).unwrap();
|
||
sys.add_edge(n1, n0).unwrap();
|
||
sys.finalize().unwrap();
|
||
|
||
let state = vec![0.0; sys.state_vector_len()];
|
||
|
||
// Capture log output using tracing_subscriber
|
||
let log_buffer = std::sync::Arc::new(std::sync::Mutex::new(String::new()));
|
||
let buffer_clone = log_buffer.clone();
|
||
let layer = tracing_subscriber::fmt::layer()
|
||
.with_writer(move || {
|
||
use std::io::Write;
|
||
struct BufWriter {
|
||
buf: std::sync::Arc<std::sync::Mutex<String>>,
|
||
}
|
||
impl Write for BufWriter {
|
||
fn write(&mut self, data: &[u8]) -> std::io::Result<usize> {
|
||
let mut buf = self.buf.lock().unwrap();
|
||
buf.push_str(&String::from_utf8_lossy(data));
|
||
Ok(data.len())
|
||
}
|
||
fn flush(&mut self) -> std::io::Result<()> {
|
||
Ok(())
|
||
}
|
||
}
|
||
BufWriter {
|
||
buf: buffer_clone.clone(),
|
||
}
|
||
})
|
||
.without_time();
|
||
|
||
let _guard = tracing_subscriber::registry().with(layer).set_default();
|
||
|
||
// check_energy_balance should succeed (no violations) but will emit warnings
|
||
// for components that lack energy_transfers() and port_enthalpies()
|
||
let result = sys.check_energy_balance(&state);
|
||
assert!(
|
||
result.is_ok(),
|
||
"check_energy_balance should succeed even with skipped components"
|
||
);
|
||
|
||
// Verify warning was emitted
|
||
let log_output = log_buffer.lock().unwrap();
|
||
assert!(
|
||
log_output.contains("SKIPPED in energy balance validation"),
|
||
"Expected warning message not found in logs. Actual output: {}",
|
||
*log_output
|
||
);
|
||
}
|
||
|
||
/// Test that cycle_performance classifies duties by physical sign and closes
|
||
/// the First Law (Q_heating = Q_cooling + W_input).
|
||
#[test]
|
||
fn test_cycle_performance_classifies_and_closes_first_law() {
|
||
// Compressor (work in), condenser (heat out), evaporator (heat in) in a loop.
|
||
let mut sys = System::new();
|
||
let comp = sys.add_component(Box::new(EnergyMock {
|
||
heat_w: 0.0,
|
||
work_w: -2000.0,
|
||
}));
|
||
let cond = sys.add_component(Box::new(EnergyMock {
|
||
heat_w: -9000.0,
|
||
work_w: 0.0,
|
||
}));
|
||
let evap = sys.add_component(Box::new(EnergyMock {
|
||
heat_w: 7000.0,
|
||
work_w: 0.0,
|
||
}));
|
||
sys.add_edge(comp, cond).unwrap();
|
||
sys.add_edge(cond, evap).unwrap();
|
||
sys.add_edge(evap, comp).unwrap();
|
||
sys.finalize().unwrap();
|
||
|
||
let state = vec![0.0; sys.state_vector_len()];
|
||
let perf = sys
|
||
.cycle_performance(&state)
|
||
.expect("performance available");
|
||
|
||
assert_relative_eq!(perf.q_cooling_w, 7000.0, epsilon = 1e-9);
|
||
assert_relative_eq!(perf.q_heating_w, 9000.0, epsilon = 1e-9);
|
||
assert_relative_eq!(perf.work_input_w, 2000.0, epsilon = 1e-9);
|
||
assert_eq!(perf.components_counted, 3);
|
||
// First Law: rejected = absorbed + work.
|
||
assert_relative_eq!(
|
||
perf.q_heating_w,
|
||
perf.q_cooling_w + perf.work_input_w,
|
||
epsilon = 1e-9
|
||
);
|
||
assert_relative_eq!(perf.cop_cooling().unwrap(), 3.5, epsilon = 1e-9);
|
||
assert_relative_eq!(perf.cop_heating().unwrap(), 4.5, epsilon = 1e-9);
|
||
}
|
||
|
||
/// cycle_performance returns None when no component reports energy transfers.
|
||
#[test]
|
||
fn test_cycle_performance_none_without_energy_reporters() {
|
||
let mut sys = System::new();
|
||
let n0 = sys.add_component(make_mock(0));
|
||
let n1 = sys.add_component(make_mock(0));
|
||
sys.add_edge(n0, n1).unwrap();
|
||
sys.add_edge(n1, n0).unwrap();
|
||
sys.finalize().unwrap();
|
||
|
||
let state = vec![0.0; sys.state_vector_len()];
|
||
assert!(sys.cycle_performance(&state).is_none());
|
||
}
|
||
|
||
/// Test that check_energy_balance includes component type in warning message.
|
||
#[test]
|
||
fn test_energy_balance_includes_component_type_in_warning() {
|
||
use tracing_subscriber::layer::SubscriberExt;
|
||
use tracing_subscriber::util::SubscriberInitExt;
|
||
|
||
// Create a system with mock components (need at least 2 nodes with edges to avoid isolated node error)
|
||
let mut sys = System::new();
|
||
let n0 = sys.add_component(make_mock(0));
|
||
let n1 = sys.add_component(make_mock(0));
|
||
sys.add_edge(n0, n1).unwrap();
|
||
sys.add_edge(n1, n0).unwrap();
|
||
sys.finalize().unwrap();
|
||
|
||
let state = vec![0.0; sys.state_vector_len()];
|
||
|
||
// Capture log output using tracing_subscriber
|
||
let log_buffer = std::sync::Arc::new(std::sync::Mutex::new(String::new()));
|
||
let buffer_clone = log_buffer.clone();
|
||
let layer = tracing_subscriber::fmt::layer()
|
||
.with_writer(move || {
|
||
use std::io::Write;
|
||
struct BufWriter {
|
||
buf: std::sync::Arc<std::sync::Mutex<String>>,
|
||
}
|
||
impl Write for BufWriter {
|
||
fn write(&mut self, data: &[u8]) -> std::io::Result<usize> {
|
||
let mut buf = self.buf.lock().unwrap();
|
||
buf.push_str(&String::from_utf8_lossy(data));
|
||
Ok(data.len())
|
||
}
|
||
fn flush(&mut self) -> std::io::Result<()> {
|
||
Ok(())
|
||
}
|
||
}
|
||
BufWriter {
|
||
buf: buffer_clone.clone(),
|
||
}
|
||
})
|
||
.without_time();
|
||
|
||
let _guard = tracing_subscriber::registry().with(layer).set_default();
|
||
|
||
let result = sys.check_energy_balance(&state);
|
||
assert!(result.is_ok());
|
||
|
||
// Verify warning message includes component type information
|
||
// Note: type_name_of_val on a trait object returns the trait name ("Component"),
|
||
// not the concrete type. This is a known Rust limitation.
|
||
let log_output = log_buffer.lock().unwrap();
|
||
assert!(
|
||
log_output.contains("type: Component"),
|
||
"Expected component type information not found in logs. Actual output: {}",
|
||
*log_output
|
||
);
|
||
}
|
||
|
||
/// Test that check_energy_balance emits a summary warning with skipped component count.
|
||
#[test]
|
||
fn test_energy_balance_summary_warning() {
|
||
use tracing_subscriber::layer::SubscriberExt;
|
||
use tracing_subscriber::util::SubscriberInitExt;
|
||
|
||
// Create a system with mock components
|
||
let mut sys = System::new();
|
||
let n0 = sys.add_component(make_mock(0));
|
||
let n1 = sys.add_component(make_mock(0));
|
||
sys.add_edge(n0, n1).unwrap();
|
||
sys.add_edge(n1, n0).unwrap();
|
||
sys.finalize().unwrap();
|
||
|
||
let state = vec![0.0; sys.state_vector_len()];
|
||
|
||
// Capture log output
|
||
let log_buffer = std::sync::Arc::new(std::sync::Mutex::new(String::new()));
|
||
let buffer_clone = log_buffer.clone();
|
||
let layer = tracing_subscriber::fmt::layer()
|
||
.with_writer(move || {
|
||
use std::io::Write;
|
||
struct BufWriter {
|
||
buf: std::sync::Arc<std::sync::Mutex<String>>,
|
||
}
|
||
impl Write for BufWriter {
|
||
fn write(&mut self, data: &[u8]) -> std::io::Result<usize> {
|
||
let mut buf = self.buf.lock().unwrap();
|
||
buf.push_str(&String::from_utf8_lossy(data));
|
||
Ok(data.len())
|
||
}
|
||
fn flush(&mut self) -> std::io::Result<()> {
|
||
Ok(())
|
||
}
|
||
}
|
||
BufWriter {
|
||
buf: buffer_clone.clone(),
|
||
}
|
||
})
|
||
.without_time();
|
||
|
||
let _guard = tracing_subscriber::registry().with(layer).set_default();
|
||
|
||
let result = sys.check_energy_balance(&state);
|
||
assert!(result.is_ok());
|
||
|
||
// Verify summary warning was emitted
|
||
let log_output = log_buffer.lock().unwrap();
|
||
assert!(
|
||
log_output.contains("Energy balance validation incomplete"),
|
||
"Expected summary warning not found in logs. Actual output: {}",
|
||
*log_output
|
||
);
|
||
assert!(
|
||
log_output.contains("component(s) skipped"),
|
||
"Expected 'component(s) skipped' not found in logs. Actual output: {}",
|
||
*log_output
|
||
);
|
||
}
|
||
}
|