feat(python): implement python bindings for all components and solvers
This commit is contained in:
Sepehr
@@ -104,6 +104,14 @@ pub enum ComponentOutput {
|
||||
component_id: String,
|
||||
},
|
||||
|
||||
/// Capacity (W).
|
||||
///
|
||||
/// Cooling or heating capacity of a component.
|
||||
Capacity {
|
||||
/// Component identifier
|
||||
component_id: String,
|
||||
},
|
||||
|
||||
/// Mass flow rate (kg/s).
|
||||
///
|
||||
/// Mass flow through a component.
|
||||
@@ -133,6 +141,7 @@ impl ComponentOutput {
|
||||
ComponentOutput::Superheat { component_id } => component_id,
|
||||
ComponentOutput::Subcooling { component_id } => component_id,
|
||||
ComponentOutput::HeatTransferRate { component_id } => component_id,
|
||||
ComponentOutput::Capacity { component_id } => component_id,
|
||||
ComponentOutput::MassFlowRate { component_id } => component_id,
|
||||
ComponentOutput::Pressure { component_id } => component_id,
|
||||
ComponentOutput::Temperature { component_id } => component_id,
|
||||
|
||||
@@ -175,7 +175,7 @@ impl ControlMapping {
|
||||
///
|
||||
/// Manages constraint-to-control-variable mappings for embedding constraints
|
||||
/// into the residual system.
|
||||
#[derive(Debug, Clone, Default)]
|
||||
#[derive(Debug, Clone)]
|
||||
pub struct InverseControlConfig {
|
||||
/// Mapping from constraint ID to control variable ID.
|
||||
constraint_to_control: HashMap<ConstraintId, BoundedVariableId>,
|
||||
@@ -183,15 +183,28 @@ pub struct InverseControlConfig {
|
||||
control_to_constraint: HashMap<BoundedVariableId, ConstraintId>,
|
||||
/// Whether inverse control is enabled globally.
|
||||
enabled: bool,
|
||||
/// Finite difference epsilon for numerical Jacobian computation.
|
||||
/// Default is 1e-6, which balances numerical precision against floating-point rounding errors.
|
||||
finite_diff_epsilon: f64,
|
||||
}
|
||||
|
||||
impl Default for InverseControlConfig {
|
||||
fn default() -> Self {
|
||||
Self::new()
|
||||
}
|
||||
}
|
||||
|
||||
impl InverseControlConfig {
|
||||
/// Default finite difference epsilon for numerical Jacobian computation.
|
||||
pub const DEFAULT_FINITE_DIFF_EPSILON: f64 = 1e-6;
|
||||
|
||||
/// Creates a new empty inverse control configuration.
|
||||
pub fn new() -> Self {
|
||||
InverseControlConfig {
|
||||
constraint_to_control: HashMap::new(),
|
||||
control_to_constraint: HashMap::new(),
|
||||
enabled: true,
|
||||
finite_diff_epsilon: Self::DEFAULT_FINITE_DIFF_EPSILON,
|
||||
}
|
||||
}
|
||||
|
||||
@@ -201,9 +214,25 @@ impl InverseControlConfig {
|
||||
constraint_to_control: HashMap::new(),
|
||||
control_to_constraint: HashMap::new(),
|
||||
enabled: false,
|
||||
finite_diff_epsilon: Self::DEFAULT_FINITE_DIFF_EPSILON,
|
||||
}
|
||||
}
|
||||
|
||||
/// Returns the finite difference epsilon used for numerical Jacobian computation.
|
||||
pub fn finite_diff_epsilon(&self) -> f64 {
|
||||
self.finite_diff_epsilon
|
||||
}
|
||||
|
||||
/// Sets the finite difference epsilon for numerical Jacobian computation.
|
||||
///
|
||||
/// # Panics
|
||||
///
|
||||
/// Panics if epsilon is non-positive.
|
||||
pub fn set_finite_diff_epsilon(&mut self, epsilon: f64) {
|
||||
assert!(epsilon > 0.0, "Finite difference epsilon must be positive");
|
||||
self.finite_diff_epsilon = epsilon;
|
||||
}
|
||||
|
||||
/// Returns whether inverse control is enabled.
|
||||
pub fn is_enabled(&self) -> bool {
|
||||
self.enabled
|
||||
|
||||
@@ -370,14 +370,14 @@ impl JacobianMatrix {
|
||||
// This optimizes the check from O(N^2 * C) to O(N^2)
|
||||
let mut row_block_cols = vec![None; nrows];
|
||||
for &(rs, re, cs, ce) in &blocks {
|
||||
for r in rs..re {
|
||||
row_block_cols[r] = Some((cs, ce));
|
||||
for block in &mut row_block_cols[rs..re] {
|
||||
*block = Some((cs, ce));
|
||||
}
|
||||
}
|
||||
|
||||
for row in 0..nrows {
|
||||
for (row, block) in row_block_cols.iter().enumerate().take(nrows) {
|
||||
for col in 0..ncols {
|
||||
let in_block = match row_block_cols[row] {
|
||||
let in_block = match *block {
|
||||
Some((cs, ce)) => col >= cs && col < ce,
|
||||
None => false,
|
||||
};
|
||||
|
||||
@@ -438,6 +438,13 @@ pub struct NewtonConfig {
|
||||
/// This is useful for HIL scenarios where the last known-good state should be used.
|
||||
pub previous_state: Option<Vec<f64>>,
|
||||
|
||||
/// Residual norm associated with `previous_state` for ZOH fallback (Story 4.5).
|
||||
///
|
||||
/// When using ZOH fallback, this residual is returned instead of `best_residual`,
|
||||
/// ensuring the returned state and residual are consistent.
|
||||
/// Should be set alongside `previous_state` by the HIL controller.
|
||||
pub previous_residual: Option<f64>,
|
||||
|
||||
/// Smart initial state for cold-start solving (Story 4.6).
|
||||
///
|
||||
/// When `Some`, the solver starts from this state instead of the zero vector.
|
||||
@@ -478,6 +485,7 @@ impl Default for NewtonConfig {
|
||||
divergence_threshold: 1e10,
|
||||
timeout_config: TimeoutConfig::default(),
|
||||
previous_state: None,
|
||||
previous_residual: None,
|
||||
initial_state: None,
|
||||
convergence_criteria: None,
|
||||
jacobian_freezing: None,
|
||||
@@ -530,7 +538,7 @@ impl NewtonConfig {
|
||||
/// - Previous state (ZOH) if `zoh_fallback` is true and previous state available
|
||||
fn handle_timeout(
|
||||
&self,
|
||||
best_state: Vec<f64>,
|
||||
best_state: &[f64],
|
||||
best_residual: f64,
|
||||
iterations: usize,
|
||||
timeout: Duration,
|
||||
@@ -545,15 +553,16 @@ impl NewtonConfig {
|
||||
// If ZOH fallback is enabled and previous state is available
|
||||
if self.timeout_config.zoh_fallback {
|
||||
if let Some(ref prev_state) = self.previous_state {
|
||||
let residual = self.previous_residual.unwrap_or(best_residual);
|
||||
tracing::info!(
|
||||
iterations = iterations,
|
||||
best_residual = best_residual,
|
||||
residual = residual,
|
||||
"Returning previous state (ZOH fallback) on timeout"
|
||||
);
|
||||
return Ok(ConvergedState::new(
|
||||
prev_state.clone(),
|
||||
iterations,
|
||||
best_residual,
|
||||
residual,
|
||||
ConvergenceStatus::TimedOutWithBestState,
|
||||
));
|
||||
}
|
||||
@@ -566,7 +575,7 @@ impl NewtonConfig {
|
||||
"Returning best state on timeout"
|
||||
);
|
||||
Ok(ConvergedState::new(
|
||||
best_state,
|
||||
best_state.to_vec(),
|
||||
iterations,
|
||||
best_residual,
|
||||
ConvergenceStatus::TimedOutWithBestState,
|
||||
@@ -623,6 +632,7 @@ impl NewtonConfig {
|
||||
///
|
||||
/// This method requires pre-allocated buffers to avoid heap allocation in the
|
||||
/// hot path. `state_copy` and `new_residuals` must have appropriate lengths.
|
||||
#[allow(clippy::too_many_arguments)]
|
||||
fn line_search(
|
||||
&self,
|
||||
system: &System,
|
||||
@@ -630,8 +640,9 @@ impl NewtonConfig {
|
||||
delta: &[f64],
|
||||
_residuals: &[f64],
|
||||
current_norm: f64,
|
||||
state_copy: &mut Vec<f64>,
|
||||
state_copy: &mut [f64],
|
||||
new_residuals: &mut Vec<f64>,
|
||||
clipping_mask: &[Option<(f64, f64)>],
|
||||
) -> Option<f64> {
|
||||
let mut alpha: f64 = 1.0;
|
||||
state_copy.copy_from_slice(state);
|
||||
@@ -641,9 +652,7 @@ impl NewtonConfig {
|
||||
|
||||
for _backtrack in 0..self.line_search_max_backtracks {
|
||||
// Apply step: x = x + alpha * delta
|
||||
for (s, &d) in state.iter_mut().zip(delta.iter()) {
|
||||
*s = *s + alpha * d;
|
||||
}
|
||||
apply_newton_step(state, delta, clipping_mask, alpha);
|
||||
|
||||
// Compute new residuals (uses pre-allocated buffer)
|
||||
if system.compute_residuals(state, new_residuals).is_err() {
|
||||
@@ -680,6 +689,24 @@ impl NewtonConfig {
|
||||
}
|
||||
}
|
||||
|
||||
/// Applies a Newton step to the state vector, clamping bounded variables.
|
||||
///
|
||||
/// Update formula: x_new = clamp(x_old + alpha * delta)
|
||||
fn apply_newton_step(
|
||||
state: &mut [f64],
|
||||
delta: &[f64],
|
||||
clipping_mask: &[Option<(f64, f64)>],
|
||||
alpha: f64,
|
||||
) {
|
||||
for (i, s) in state.iter_mut().enumerate() {
|
||||
let proposed = *s + alpha * delta[i];
|
||||
*s = match &clipping_mask[i] {
|
||||
Some((min, max)) => proposed.clamp(*min, *max),
|
||||
None => proposed,
|
||||
};
|
||||
}
|
||||
}
|
||||
|
||||
impl Solver for NewtonConfig {
|
||||
fn solve(&mut self, system: &mut System) -> Result<ConvergedState, SolverError> {
|
||||
let start_time = Instant::now();
|
||||
@@ -750,6 +777,11 @@ impl Solver for NewtonConfig {
|
||||
let mut frozen_count: usize = 0;
|
||||
let mut force_recompute: bool = true; // Always compute on the very first iteration
|
||||
|
||||
// Pre-compute clipping mask (Story 5.6)
|
||||
let clipping_mask: Vec<Option<(f64, f64)>> = (0..n_state)
|
||||
.map(|i| system.get_bounds_for_state_index(i))
|
||||
.collect();
|
||||
|
||||
// Initial residual computation
|
||||
system
|
||||
.compute_residuals(&state, &mut residuals)
|
||||
@@ -783,7 +815,11 @@ impl Solver for NewtonConfig {
|
||||
"System already converged at initial state (criteria)"
|
||||
);
|
||||
return Ok(ConvergedState::with_report(
|
||||
state, 0, current_norm, status, report,
|
||||
state,
|
||||
0,
|
||||
current_norm,
|
||||
status,
|
||||
report,
|
||||
));
|
||||
}
|
||||
} else {
|
||||
@@ -792,9 +828,7 @@ impl Solver for NewtonConfig {
|
||||
final_residual = current_norm,
|
||||
"System already converged at initial state"
|
||||
);
|
||||
return Ok(ConvergedState::new(
|
||||
state, 0, current_norm, status,
|
||||
));
|
||||
return Ok(ConvergedState::new(state, 0, current_norm, status));
|
||||
}
|
||||
}
|
||||
|
||||
@@ -815,7 +849,7 @@ impl Solver for NewtonConfig {
|
||||
);
|
||||
|
||||
// Story 4.5 - AC: #2, #6: Return best state or error based on config
|
||||
return self.handle_timeout(best_state, best_residual, iteration - 1, timeout);
|
||||
return self.handle_timeout(&best_state, best_residual, iteration - 1, timeout);
|
||||
}
|
||||
}
|
||||
|
||||
@@ -905,6 +939,7 @@ impl Solver for NewtonConfig {
|
||||
current_norm,
|
||||
&mut state_copy,
|
||||
&mut new_residuals,
|
||||
&clipping_mask,
|
||||
) {
|
||||
Some(a) => a,
|
||||
None => {
|
||||
@@ -915,9 +950,7 @@ impl Solver for NewtonConfig {
|
||||
}
|
||||
} else {
|
||||
// Full Newton step: x = x + delta (delta already includes negative sign)
|
||||
for (s, &d) in state.iter_mut().zip(delta.iter()) {
|
||||
*s = *s + d;
|
||||
}
|
||||
apply_newton_step(&mut state, &delta, &clipping_mask, 1.0);
|
||||
1.0
|
||||
};
|
||||
|
||||
@@ -988,7 +1021,11 @@ impl Solver for NewtonConfig {
|
||||
"Newton-Raphson converged (criteria)"
|
||||
);
|
||||
return Ok(ConvergedState::with_report(
|
||||
state, iteration, current_norm, status, report,
|
||||
state,
|
||||
iteration,
|
||||
current_norm,
|
||||
status,
|
||||
report,
|
||||
));
|
||||
}
|
||||
false
|
||||
@@ -1007,9 +1044,7 @@ impl Solver for NewtonConfig {
|
||||
final_residual = current_norm,
|
||||
"Newton-Raphson converged"
|
||||
);
|
||||
return Ok(ConvergedState::new(
|
||||
state, iteration, current_norm, status,
|
||||
));
|
||||
return Ok(ConvergedState::new(state, iteration, current_norm, status));
|
||||
}
|
||||
|
||||
// Check divergence (AC: #5)
|
||||
@@ -1099,6 +1134,13 @@ pub struct PicardConfig {
|
||||
/// This is useful for HIL scenarios where the last known-good state should be used.
|
||||
pub previous_state: Option<Vec<f64>>,
|
||||
|
||||
/// Residual norm associated with `previous_state` for ZOH fallback (Story 4.5).
|
||||
///
|
||||
/// When using ZOH fallback, this residual is returned instead of `best_residual`,
|
||||
/// ensuring the returned state and residual are consistent.
|
||||
/// Should be set alongside `previous_state` by the HIL controller.
|
||||
pub previous_residual: Option<f64>,
|
||||
|
||||
/// Smart initial state for cold-start solving (Story 4.6).
|
||||
///
|
||||
/// When `Some`, the solver starts from this state instead of the zero vector.
|
||||
@@ -1128,6 +1170,7 @@ impl Default for PicardConfig {
|
||||
divergence_patience: 5,
|
||||
timeout_config: TimeoutConfig::default(),
|
||||
previous_state: None,
|
||||
previous_residual: None,
|
||||
initial_state: None,
|
||||
convergence_criteria: None,
|
||||
}
|
||||
@@ -1167,7 +1210,7 @@ impl PicardConfig {
|
||||
/// - Previous state (ZOH) if `zoh_fallback` is true and previous state available
|
||||
fn handle_timeout(
|
||||
&self,
|
||||
best_state: Vec<f64>,
|
||||
best_state: &[f64],
|
||||
best_residual: f64,
|
||||
iterations: usize,
|
||||
timeout: Duration,
|
||||
@@ -1182,15 +1225,16 @@ impl PicardConfig {
|
||||
// If ZOH fallback is enabled and previous state is available
|
||||
if self.timeout_config.zoh_fallback {
|
||||
if let Some(ref prev_state) = self.previous_state {
|
||||
let residual = self.previous_residual.unwrap_or(best_residual);
|
||||
tracing::info!(
|
||||
iterations = iterations,
|
||||
best_residual = best_residual,
|
||||
residual = residual,
|
||||
"Returning previous state (ZOH fallback) on timeout"
|
||||
);
|
||||
return Ok(ConvergedState::new(
|
||||
prev_state.clone(),
|
||||
iterations,
|
||||
best_residual,
|
||||
residual,
|
||||
ConvergenceStatus::TimedOutWithBestState,
|
||||
));
|
||||
}
|
||||
@@ -1203,7 +1247,7 @@ impl PicardConfig {
|
||||
"Returning best state on timeout"
|
||||
);
|
||||
Ok(ConvergedState::new(
|
||||
best_state,
|
||||
best_state.to_vec(),
|
||||
iterations,
|
||||
best_residual,
|
||||
ConvergenceStatus::TimedOutWithBestState,
|
||||
@@ -1257,7 +1301,7 @@ impl PicardConfig {
|
||||
/// This is the standard Picard iteration: x_{k+1} = x_k - ω·F(x_k)
|
||||
fn apply_relaxation(state: &mut [f64], residuals: &[f64], omega: f64) {
|
||||
for (x, &r) in state.iter_mut().zip(residuals.iter()) {
|
||||
*x = *x - omega * r;
|
||||
*x -= omega * r;
|
||||
}
|
||||
}
|
||||
}
|
||||
@@ -1375,7 +1419,7 @@ impl Solver for PicardConfig {
|
||||
);
|
||||
|
||||
// Story 4.5 - AC: #2, #6: Return best state or error based on config
|
||||
return self.handle_timeout(best_state, best_residual, iteration - 1, timeout);
|
||||
return self.handle_timeout(&best_state, best_residual, iteration - 1, timeout);
|
||||
}
|
||||
}
|
||||
|
||||
@@ -2117,6 +2161,7 @@ mod tests {
|
||||
divergence_threshold: 1e10,
|
||||
timeout_config: TimeoutConfig::default(),
|
||||
previous_state: None,
|
||||
previous_residual: None,
|
||||
initial_state: None,
|
||||
convergence_criteria: None,
|
||||
jacobian_freezing: None,
|
||||
@@ -2427,6 +2472,7 @@ mod tests {
|
||||
divergence_patience: 7,
|
||||
timeout_config: TimeoutConfig::default(),
|
||||
previous_state: None,
|
||||
previous_residual: None,
|
||||
initial_state: None,
|
||||
convergence_criteria: None,
|
||||
}
|
||||
@@ -2712,4 +2758,63 @@ mod tests {
|
||||
"should not allow excessive switches"
|
||||
);
|
||||
}
|
||||
|
||||
// ─────────────────────────────────────────────────────────────────────────────
|
||||
// Story 5.6: Control Variable Step Clipping Tests
|
||||
// ─────────────────────────────────────────────────────────────────────────────
|
||||
|
||||
#[test]
|
||||
fn test_bounded_variable_clipped_at_max() {
|
||||
let mut state = vec![0.5];
|
||||
let delta = vec![2.0]; // Proposed step: 0.5 + 2.0 = 2.5
|
||||
let mask = vec![Some((0.0, 1.0))];
|
||||
super::apply_newton_step(&mut state, &delta, &mask, 1.0);
|
||||
assert_eq!(state[0], 1.0, "Should be clipped to max bound");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_bounded_variable_clipped_at_min() {
|
||||
let mut state = vec![0.5];
|
||||
let delta = vec![-2.0]; // Proposed step: 0.5 - 2.0 = -1.5
|
||||
let mask = vec![Some((0.0, 1.0))];
|
||||
super::apply_newton_step(&mut state, &delta, &mask, 1.0);
|
||||
assert_eq!(state[0], 0.0, "Should be clipped to min bound");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_edge_states_not_clipped() {
|
||||
let mut state = vec![0.5, 10.0];
|
||||
let delta = vec![-2.0, 50.0];
|
||||
// Only first variable is bounded
|
||||
let mask = vec![Some((0.0, 1.0)), None];
|
||||
super::apply_newton_step(&mut state, &delta, &mask, 1.0);
|
||||
assert_eq!(state[0], 0.0, "Bounded variable should be clipped");
|
||||
assert_eq!(state[1], 60.0, "Unbounded variable should NOT be clipped");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_saturation_detected_after_convergence() {
|
||||
use crate::inverse::{BoundedVariable, BoundedVariableId, SaturationType};
|
||||
|
||||
let mut sys = System::new();
|
||||
// A saturated variable (value = max bound)
|
||||
sys.add_bounded_variable(
|
||||
BoundedVariable::new(BoundedVariableId::new("v1"), 1.0, 0.0, 1.0).unwrap(),
|
||||
)
|
||||
.unwrap();
|
||||
// An unsaturated variable
|
||||
sys.add_bounded_variable(
|
||||
BoundedVariable::new(BoundedVariableId::new("v2"), 0.5, 0.0, 1.0).unwrap(),
|
||||
)
|
||||
.unwrap();
|
||||
|
||||
let saturated = sys.saturated_variables();
|
||||
assert_eq!(saturated.len(), 1, "Should detect 1 saturated variable");
|
||||
assert_eq!(
|
||||
saturated[0].saturation_type,
|
||||
SaturationType::UpperBound,
|
||||
"Variable v1 should be saturated at max"
|
||||
);
|
||||
assert_eq!(saturated[0].variable_id.as_str(), "v1");
|
||||
}
|
||||
}
|
||||
|
||||
@@ -353,8 +353,12 @@ impl System {
|
||||
let mut current_offset = 2 * self.graph.edge_count();
|
||||
|
||||
// Gather (node_idx, offset, incident_edge_indices) before mutating nodes.
|
||||
let mut node_context: Vec<(petgraph::graph::NodeIndex, usize, Vec<(usize, usize)>)> =
|
||||
Vec::new();
|
||||
#[allow(clippy::type_complexity)]
|
||||
let mut node_context: Vec<(
|
||||
petgraph::graph::NodeIndex,
|
||||
usize,
|
||||
Vec<(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)> = Vec::new();
|
||||
@@ -380,15 +384,46 @@ impl System {
|
||||
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 id_str.ends_with("f_m") || id_str == "f_m" {
|
||||
indices.f_m = Some(state_idx);
|
||||
} else if id_str.ends_with("f_dp") || id_str == "f_dp" {
|
||||
indices.f_dp = Some(state_idx);
|
||||
} else if id_str.ends_with("f_ua") || id_str == "f_ua" {
|
||||
indices.f_ua = Some(state_idx);
|
||||
} else if id_str.ends_with("f_power") || id_str == "f_power" {
|
||||
indices.f_power = Some(state_idx);
|
||||
} else if id_str.ends_with("f_etav") || id_str == "f_etav" {
|
||||
indices.f_etav = Some(state_idx);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Now mutate each node weight (component) with the gathered context.
|
||||
for (node_idx, offset, incident) in node_context {
|
||||
if let Some(component) = self.graph.node_weight_mut(node_idx) {
|
||||
component.set_system_context(offset, &incident);
|
||||
|
||||
// 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);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
self.total_state_len = current_offset;
|
||||
|
||||
if !self.constraints.is_empty() {
|
||||
match self.validate_inverse_control_dof() {
|
||||
Ok(()) => {
|
||||
@@ -484,18 +519,17 @@ impl System {
|
||||
"[P_edge0, h_edge0, P_edge1, h_edge1, ...] — 2 per edge (pressure Pa, enthalpy J/kg)"
|
||||
}
|
||||
|
||||
/// Returns the length of the state vector: `2 * edge_count`.
|
||||
/// Returns the length of the state vector: `2 * edge_count + internal_components_length`.
|
||||
///
|
||||
/// Note: This returns only the edge state length. For the full state vector
|
||||
/// including internal component states and control variables, use
|
||||
/// [`full_state_vector_len`](Self::full_state_vector_len).
|
||||
/// 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()");
|
||||
2 * self.graph.edge_count()
|
||||
self.total_state_len
|
||||
}
|
||||
|
||||
/// Returns the state indices (P, h) for the given edge.
|
||||
@@ -814,13 +848,13 @@ impl System {
|
||||
///
|
||||
/// ```rust,ignore
|
||||
/// let mut residuals = ResidualVector::new();
|
||||
/// let measured = system.extract_constraint_values(&state);
|
||||
/// 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 ResidualVector,
|
||||
residuals: &mut [f64],
|
||||
measured_values: &HashMap<ConstraintId, f64>,
|
||||
) -> usize {
|
||||
if self.constraints.is_empty() {
|
||||
@@ -840,42 +874,147 @@ impl System {
|
||||
constraint.target_value()
|
||||
});
|
||||
let residual = constraint.compute_residual(measured);
|
||||
residuals.push(residual);
|
||||
if count < residuals.len() {
|
||||
residuals[count] = residual;
|
||||
}
|
||||
count += 1;
|
||||
}
|
||||
count
|
||||
}
|
||||
|
||||
/// Extracts constraint output values from component state.
|
||||
/// Extracts measured values for all constraints, incorporating control variable effects.
|
||||
///
|
||||
/// This method attempts to extract measurable output values for all constraints
|
||||
/// from the current system state. For complex outputs (superheat, subcooling),
|
||||
/// additional thermodynamic calculations may be needed.
|
||||
/// 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)
|
||||
/// * `state` - Current system state (edge pressures and enthalpies)
|
||||
/// * `control_values` - Current values of control variables
|
||||
///
|
||||
/// # Returns
|
||||
///
|
||||
/// A map from constraint IDs to their measured values. Constraints whose
|
||||
/// outputs cannot be extracted will not appear in the map.
|
||||
/// A map from constraint ID to measured output value.
|
||||
///
|
||||
/// # Note
|
||||
/// # Cross-Coupling for MIMO Systems
|
||||
///
|
||||
/// Full implementation requires integration with ThermoState (Story 2.8) and
|
||||
/// component-specific output extraction. This MVP version returns an empty map
|
||||
/// and should be enhanced with actual component state extraction.
|
||||
pub fn extract_constraint_values(&self, _state: &StateSlice) -> HashMap<ConstraintId, f64> {
|
||||
/// 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.
|
||||
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 HashMap::new();
|
||||
return measured;
|
||||
}
|
||||
|
||||
tracing::debug!(
|
||||
constraint_count = self.constraints.len(),
|
||||
"Constraint value extraction called - MVP returns empty map"
|
||||
);
|
||||
HashMap::new()
|
||||
// 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) {
|
||||
// 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(&(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;
|
||||
}
|
||||
}
|
||||
|
||||
measured.insert(constraint.id().clone(), value);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
measured
|
||||
}
|
||||
|
||||
/// Computes the Jacobian entries for inverse control constraints.
|
||||
@@ -886,9 +1025,9 @@ impl System {
|
||||
///
|
||||
/// # Arguments
|
||||
///
|
||||
/// * `_state` - Current system state
|
||||
/// * `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)
|
||||
/// * `control_values` - Current values of control variables (for finite difference)
|
||||
///
|
||||
/// # Returns
|
||||
///
|
||||
@@ -898,11 +1037,16 @@ impl System {
|
||||
///
|
||||
/// 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,
|
||||
state: &StateSlice,
|
||||
row_offset: usize,
|
||||
_control_values: &[f64],
|
||||
control_values: &[f64],
|
||||
) -> Vec<(usize, usize, f64)> {
|
||||
let mut entries = Vec::new();
|
||||
|
||||
@@ -910,18 +1054,118 @@ impl System {
|
||||
return entries;
|
||||
}
|
||||
|
||||
for (i, (_constraint_id, bounded_var_id)) in self.inverse_control.mappings().enumerate() {
|
||||
let col = self.control_variable_state_index(bounded_var_id);
|
||||
if let Some(col_idx) = col {
|
||||
// Use configurable epsilon from InverseControlConfig
|
||||
let eps = self.inverse_control.finite_diff_epsilon();
|
||||
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
|
||||
for edge in self
|
||||
.graph
|
||||
.edges_directed(node_idx, petgraph::Direction::Incoming)
|
||||
{
|
||||
if let Some(&(p_idx, h_idx)) = self.edge_to_state.get(&edge.id()) {
|
||||
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(&(p_idx, h_idx)) = self.edge_to_state.get(&edge.id()) {
|
||||
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;
|
||||
entries.push((row, col_idx, 1.0));
|
||||
tracing::trace!(
|
||||
constraint = _constraint_id.as_str(),
|
||||
control = bounded_var_id.as_str(),
|
||||
row,
|
||||
col = col_idx,
|
||||
"Inverse control Jacobian entry (placeholder derivative = 1.0)"
|
||||
);
|
||||
|
||||
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"
|
||||
);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
@@ -1131,6 +1375,20 @@ impl System {
|
||||
self.inverse_control.mapping_count()
|
||||
}
|
||||
|
||||
/// Sets the finite difference epsilon for inverse control Jacobian computation.
|
||||
///
|
||||
/// # Panics
|
||||
///
|
||||
/// Panics if epsilon is non-positive.
|
||||
pub fn set_inverse_control_epsilon(&mut self, epsilon: f64) {
|
||||
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()
|
||||
@@ -1224,16 +1482,36 @@ impl System {
|
||||
}
|
||||
|
||||
let base = self.total_state_len;
|
||||
let mut index = 0;
|
||||
for linked_id in self.inverse_control.linked_controls() {
|
||||
for (index, linked_id) in self.inverse_control.linked_controls().enumerate() {
|
||||
if linked_id == id {
|
||||
return Some(base + index);
|
||||
}
|
||||
index += 1;
|
||||
}
|
||||
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))
|
||||
}
|
||||
|
||||
/// 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 the total state vector length including control variables.
|
||||
///
|
||||
/// ```text
|
||||
@@ -1399,7 +1677,7 @@ impl System {
|
||||
.map(|(_, c, _)| c.n_equations())
|
||||
.sum();
|
||||
total_eqs += self.constraints.len() + self.coupling_residual_count();
|
||||
|
||||
|
||||
if residuals.len() < total_eqs {
|
||||
return Err(ComponentError::InvalidResidualDimensions {
|
||||
expected: total_eqs,
|
||||
@@ -1419,13 +1697,15 @@ impl System {
|
||||
}
|
||||
|
||||
// Add constraints
|
||||
let measured = self.extract_constraint_values(state);
|
||||
let mut constraint_res = vec![];
|
||||
let n_constraints = self.compute_constraint_residuals(state, &mut constraint_res, &measured);
|
||||
if n_constraints > 0 {
|
||||
residuals[eq_offset..eq_offset + n_constraints].copy_from_slice(&constraint_res[0..n_constraints]);
|
||||
eq_offset += n_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);
|
||||
eq_offset += n_constraints;
|
||||
|
||||
// Add couplings
|
||||
let n_couplings = self.coupling_residual_count();
|
||||
@@ -1464,11 +1744,13 @@ impl System {
|
||||
}
|
||||
|
||||
// Add constraints jacobian
|
||||
let control_values: Vec<f64> = self.control_variable_indices()
|
||||
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);
|
||||
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);
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user