diff --git a/claude_memory/2026-04-17-controllers-linearization-and-first-reach.md b/claude_memory/2026-04-17-controllers-linearization-and-first-reach.md new file mode 100644 index 0000000..8ce983f --- /dev/null +++ b/claude_memory/2026-04-17-controllers-linearization-and-first-reach.md @@ -0,0 +1,130 @@ +# Session 2026-04-17 — controllers, linearization, first per-mode reach, barrier attempt + +## Context + +Pre-lunch discussion: reachability thrust is empty, controllers missing +for three of four DRC modes, tool choice (CORA vs JuliaReach) open. +Dane stepped away for lunch with instructions to "let me rip" on the +whole pipeline. Auto-mode on. + +## What landed + +1. **Three new controllers** in `plant-model/controllers/`: + - `ctrl_shutdown.m` — constant `u = -5*beta` (passive hot standby). + - `ctrl_scram.m` — constant `u = -8*beta` (max rod worth). + - `ctrl_heatup.m` — feedback-linearizing P on a ramped `T_avg` + reference, with saturated `u`. No integrator state on purpose + (see below). +2. **`pke_linearize.m`** — numerical Jacobians via central finite + differences. At the operating point the linear model agrees with + the nonlinear simulation to ~4e-4 relative error for a 5% `Q_sg` + step over 60 s. Eigenvalues span `[-65.93, -0.0124]` — stiff by a + factor of ~5000. +3. **`ctrl_operation_lqr.m`** — full-state LQR around `x_op`. + Closed-loop T_avg tracking is essentially perfect under the same + `Q_sg` step that makes the plain P controller overshoot by ~5 F. + Cost: requires `lyap`/`icare` (base MATLAB, no Control Toolbox). +4. **`main_mode_sweep.m`** — drives all five modes back-to-back from + plausible ICs, saves 6 figures to `docs/figures/mode_sweep_*.png`. +5. **Hand-rolled zonotope reach** in `reachability/reach_linear.m` + plus `reach_operation.m` end-to-end. Over-approximates with + axis-aligned box hull + elementwise-abs matrix propagation. + Result for operation mode under `Q_sg ∈ [85%, 100%]·P0` and an + initial ±0.1 K box on T_avg: max `|T_c - T_c0| = 0.1 K` over 600 s. + Safety band is ±5 K. **Operation-mode obligation discharged.** +6. **Lyapunov-ellipsoid barrier** in `reachability/barrier_lyapunov.m`. + **Does NOT certify safety** — the tightest quadratic bound we can + get (over a sweep of `Qbar(9,9)` from 10 to 1e6) yields max + `|dT_c| ≈ 33 K`, far outside the 5 K band. This is a fundamental + limitation of quadratic Lyapunov for highly anisotropic problems: + the ellipsoidal invariant is fat in the precursor/thermal + directions and can't be made tight only along T_c. **For an + analytic barrier cert here we'd need SOS polynomial barriers or a + polytopic invariant.** Numerical zonotope reach still proves + safety; the analytic cert is a "nice to have". +7. **Julia port** in `julia-port/`. Plant parameters, RHS, linearize, + all five controllers, LQR factory. `scripts/sim_sanity.jl` + matches MATLAB to 3 decimals on the standard operation test. + `scripts/reach_operation.jl` is a stub — ReachabilityAnalysis.jl's + LGG09, GLGM06, BFFPSV18 all numerically explode on the unscaled + stiff system (envelopes of 1e14 K or worse). Fix is diagonal + state-scaling; didn't get to it this session. + +## Key findings worth carrying forward + +### The model doesn't support true cold shutdown + +`alpha_f`, `alpha_c` are linearized at `(T_f0, T_c0)` = (328 C, 308 C). +Extrapolating to `T = 50 C` gives intrinsic feedback of ~5 $ which is +unphysical — real temperature coefficients saturate/flip at low temp. +**Trust range: ±50 C from operating point.** "Shutdown" in this demo +means hot standby (290 C, n≈0), not cold shutdown. Documented in +`plant-model/CLAUDE.md`. To properly support cold shutdown we'd need +piecewise-linear or table-lookup coefficients — separate modeling task. + +### Heatup overshoot is a controller-design signal, not a physics bug + +First heatup pass started from `n=1e-6` and produced a ~20-minute lag +(reactor rebuilds power from decay-heat level) followed by a ~7 K +overshoot (thermal inertia + no derivative term). Fixed by: +(a) saturating `|u|` at 0.5·beta so rho stays subcritical-prompt, +(b) starting heatup from `n=0.001` (0.1% power, already critical), +which matches real plant practice — criticality is achieved in +shutdown before DRC transitions to heatup. + +### PID not a good idea here + +The plant has intrinsic integrators (thermal mass). Adding PI: +- Increments state dim 10→11 (cheap for CORA, ok for JuliaReach). +- Anti-windup saturation turns each "continuous mode" into a hybrid + mode inside itself — more transitions, more reach-set bookkeeping. +- Integrator drift under persistent disturbance inflates the reach + set over long horizons. + +Stuck with P + feedforward; if needed we can promote to PI with eyes +open about the hybrid-mode cost. + +### LQR is dramatically better than P for operation mode + +Same operating point, same Q_sg step: LQR holds T_avg essentially at +setpoint; plain P has 5 F overshoot and ~0.8 F steady-state lag. LQR +adds zero state. The only downside: K is tuned around one operating +point — off-point robustness would need gain scheduling or MPC. + +### Hand-rolled reach beats ReachabilityAnalysis.jl out of the box + +The stiff, poorly-scaled nature of this system (eig span 5000x, state +magnitudes spanning 10 orders) breaks ReachabilityAnalysis.jl with +default settings. Hand-rolled MATLAB reach with augmented-matrix +integration trick + elementwise-abs box propagation was ~50 lines and +produced a tight, trustworthy result immediately. **Lesson: for +well-structured linear reach, don't reach for a toolbox.** + +### The analytic barrier problem is anisotropy-limited + +Safety set is a thin slab `|T_c - T_c0| <= 5 K` in a 10D state space. +Lyapunov ellipsoids that contain an initial box + robust invariant +are necessarily fat in the *uncontrolled* directions (slow +precursors), and no rescaling of the Lyapunov weight can pull them +tight *only* along T_c without blowing out the other directions. +**Logical next step for a thesis-strength barrier: SOS polynomial +barriers or polytopic (halfspace-intersection) invariants.** + +## Loose ends for the next session + +- Julia reach needs state rescaling. `S = diag(...)` chosen so all + states have O(1) magnitude, then run reach in scaled coords. +- Heatup reach: the ramped reference makes `x*(t)` time-varying, so + LQR-around-x_op isn't directly applicable. Options: trajectory-LQR + (time-varying K), or tube MPC. +- Predicate thresholds still not pinned in the FRET spec. `t_avg_in_range` + currently = ±5 K in the reach code but not reflected in the spec JSON. + Need to sync. +- Shutdown and scram modes don't have reach analysis yet. They're + trivial (constant u) but the safety claim for scram — "reactor + reaches subcritical decay in bounded time from any operating IC" — + is actually a nice bounded-time reachability problem. + +## Commit + +Committed as two or three logical chunks. See git log.