Entropyk/_bmad-output/implementation-artifacts/5-6-control-variable-step-clipping-in-solver.md

Story 5.6: Control Variable Step Clipping in Solver

Status: review

Story

As a control engineer, I want bounded control variables to be clipped at each Newton iteration, so that the solver never proposes physically impossible values (e.g., valve > 100%, frequency < min).

Acceptance Criteria

1. Step Clipping During Newton Iteration

   • Given a bounded variable with bounds [min, max]
   • When the solver computes a Newton step Δx
   • Then the new value is clipped: x_new = clamp(x + Δx, min, max)
   • And the variable never goes outside bounds during ANY iteration

2. State Vector Integration

   • Given control variables in the state vector at indices [2*edge_count, ...]
   • When the solver updates the state vector
   • Then bounded variables are clipped
   • And regular edge states (P, h) are NOT clipped

3. Saturation Detection After Convergence

   • Given a converged solution
   • When one or more bounded variables are at bounds
   • Then ConvergenceStatus::ControlSaturation is returned
   • And saturated_variables() returns the list of saturated variables

4. Performance

   • Given the clipping logic
   • When solving a system
   • Then no measurable performance degradation (< 1% overhead)

5. Backward Compatibility

   • Given existing code that doesn't use bounded variables
   • When solving
   • Then behavior is unchanged (no clipping applied)

Tasks / Subtasks

• Add get_bounded_variable_by_state_index() method to System
  • Efficient lookup from state index to bounded variable
  • Return Option<&BoundedVariable>
• Modify Newton-Raphson step application in solver.rs
  • Identify which state indices are bounded variables
  • Apply clip_step() only to bounded variable indices
  • Regular edge states pass through unchanged
• Update line search to respect bounds
  • When testing step lengths, ensure bounded vars stay in bounds
  • Or skip line search for bounded variables
• Add unit tests
  • Test clipping at max bound
  • Test clipping at min bound
  • Test multiple bounded variables
  • Test that edge states are NOT clipped
  • Test saturation detection after convergence
• Update documentation
  • Document clipping behavior in solver.rs
  • Update inverse control module docs

Dev Notes

Problem Statement

Currently in solver.rs (line ~921), the Newton step is applied without clipping:

// Current code (BUG):
for (s, &d) in state.iter_mut().zip(delta.iter()) {
    *s += d;  // No clipping for bounded variables!
}

This means:

• Valve position could become 1.5 (> 100%)
• VFD frequency could become 0.1 (< 30% minimum)
• The solver may diverge or converge to impossible solutions

Solution

Add clipping logic that:

1. Checks if state index corresponds to a bounded variable
2. If yes, applies clamp(min, max)
3. If no, applies the step normally

// Fixed code:
for (i, s) in state.iter_mut().enumerate() {
    let delta_i = delta[i];
    
    // Check if this index is a bounded variable
    if let Some((min, max)) = system.get_bounds_for_state_index(i) {
        *s = (*s + delta_i).clamp(min, max);
    } else {
        *s = *s + delta_i;
    }
}

State Vector Layout

State Vector Layout:
┌──────────────────────────────────────────────────────────────────┐
│ Index 0..2*N_edges-1  │ Edge states (P, h) - NOT clipped        │
├──────────────────────────────────────────────────────────────────┤
│ Index 2*N_edges..     │ Control variables - CLIPPED to bounds    │
├──────────────────────────────────────────────────────────────────┤
│ After controls        │ Thermal coupling temps - NOT clipped     │
└──────────────────────────────────────────────────────────────────┘

Implementation Strategy

Option A: Check on every iteration (simple)

for (i, s) in state.iter_mut().enumerate() {
    let proposed = *s + delta[i];
    *s = system.clip_state_index(i, proposed);
}

Option B: Pre-compute clipping mask (faster)

// Before Newton loop:
let clipping_mask: Vec<Option<(f64, f64)>> = (0..state.len())
    .map(|i| system.get_bounds_for_state_index(i))
    .collect();

// In Newton loop:
for (i, s) in state.iter_mut().enumerate() {
    let proposed = *s + delta[i];
    *s = match &clipping_mask[i] {
        Some((min, max)) => proposed.clamp(*min, *max),
        None => proposed,
    };
}

Recommendation: Option B - Pre-compute the mask once after finalize() for minimal overhead.

Code Locations

File	Location	Change
system.rs	New method	get_bounds_for_state_index()
solver.rs	Line ~920-925	Add clipping in Newton step
solver.rs	Line ~900-918	Consider bounds in line search

Existing Code to Leverage

From bounded.rs:

pub fn clip_step(current: f64, delta: f64, min: f64, max: f64) -> f64 {
    let proposed = current + delta;
    proposed.clamp(min, max)
}

From system.rs:

pub fn control_variable_state_index(&self, id: &BoundedVariableId) -> Option<usize>
pub fn saturated_variables(&self) -> Vec<SaturationInfo>

Anti-Patterns to Avoid

• DON'T clip edge states (P, h) - they can be any physical value
• DON'T allocate in the Newton loop - pre-compute the mask
• DON'T use unwrap() - the mask contains Option<(f64, f64)>
• DON'T forget line search - bounds must be respected there too

Testing Strategy

#[test]
fn test_bounded_variable_clipped_at_max() {
    // Valve at 95%, Newton wants +10% → should clip to 100%
    let mut system = System::new();
    // ... setup ...
    
    let valve = BoundedVariable::new(
        BoundedVariableId::new("valve"),
        0.95, 0.0, 1.0
    ).unwrap();
    system.add_bounded_variable(valve);
    
    let result = solver.solve(&system).unwrap();
    
    let final_valve = system.get_bounded_variable(&BoundedVariableId::new("valve")).unwrap();
    assert!(final_valve.value() <= 1.0);
}

#[test]
fn test_bounded_variable_clipped_at_min() {
    // VFD at 35%, Newton wants -10% → should clip to 30%
    let mut system = System::new();
    // ... setup ...
    
    let vfd = BoundedVariable::new(
        BoundedVariableId::new("vfd"),
        0.35, 0.30, 1.0  // min = 30%
    ).unwrap();
    system.add_bounded_variable(vfd);
    
    let result = solver.solve(&system).unwrap();
    
    let final_vfd = system.get_bounded_variable(&BoundedVariableId::new("vfd")).unwrap();
    assert!(final_vfd.value() >= 0.30);
}

#[test]
fn test_edge_states_not_clipped() {
    // Edge pressures can go negative in Newton step (temporary)
    // They should NOT be clipped
    // ... verify edge states can exceed any bounds ...
}

#[test]
fn test_saturation_detected_after_convergence() {
    // If solution requires valve > 100%, detect saturation
    // ... verify ConvergenceStatus::ControlSaturation ...
}

References

• [Source: bounded.rs:502-517] clip_step() function
• [Source: solver.rs:920-925] Newton step application (current bug)
• [Source: system.rs:1305-1320] control_variable_state_index()
• [Source: Story 5.2] BoundedVariable implementation
• [Source: FR23] "solving respecting Bounded Constraints"
• [Source: Story 5.2 AC] "variable never outside bounds during iterations"

Related Stories

• Story 5.2: Bounded Control Variables (DONE) - provides clip_step(), BoundedVariable
• Story 5.3: Residual Embedding (DONE) - adds controls to state vector
• Story 5.7: Prioritized Constraints (FUTURE) - protection vs performance
• Story 5.8: Coupling Functions (FUTURE) - SDT_max = f(SST)

Dev Agent Record

Agent Model Used

Antigravity

Debug Log References

• Iteration loop uses apply_newton_step to enforce bounds cleanly, avoiding heap allocations during Newton loop.
• Line search is updated to test step lengths while still honoring the bound limits.

Completion Notes List

• Implemented get_bounded_variable_by_state_index and get_bounds_for_state_index in System.
• Precomputed clipping_mask in Newton solver before loops to cache bound lookup without allocations.
• Overhauled Newton step updates via apply_newton_step to properly clamp variables defined in the clipping_mask.
• Ensured Armijo line search evaluates bounded states correctly without crashing.
• Verified physical edge variables (P, h) continue unhindered by bounds.
• Added comprehensive tests for behavior including saturation detection.

File List

• crates/solver/src/system.rs
• crates/solver/src/solver.rs
• crates/solver/tests/newton_raphson.rs