diff --git a/README.md b/README.md index 62bc8ee..1566c24 100644 --- a/README.md +++ b/README.md @@ -22,8 +22,10 @@ pwr-hybrid-3-demo/ figures/ Shared figures for thesis + talks fret-pipeline/ FRET → ltlsynt → AIGER → state machine plant-model/ PWR point kinetics + thermal-hydraulics - reachability/ Continuous-mode verification (TBD) + reachability/ Continuous-mode verification (linear-model tube + Lyapunov barrier attempt; see README) + julia-port/ Parallel plant-model port + ReachabilityAnalysis.jl scaffold hardware/ Ovation HIL artifacts (TBD) + claude_memory/ Session notes by AI agents (distilled up into CLAUDE.md over time) thesis/ [submodule] PhD proposal presentations/ 2026DICE/ [submodule] DICE 2026 abstract @@ -48,12 +50,25 @@ python3 scripts/trace_aiger.py circuits/PWR_HYBRID_3_DRC.aag diagrams dot -Tpng diagrams/PWR_HYBRID_3_DRC_states.dot -o diagrams/PWR_HYBRID_3_DRC_states.png ``` -Run the plant model (MATLAB or GNU Octave in `plant-model/`): +Run the plant model (MATLAB in `plant-model/` — Octave compatibility not tested since the LQR pieces landed): ```matlab -main +main % original single-scenario demo (null vs operation) +main_mode_sweep % all five DRC modes back-to-back, writes to ../docs/figures/ +test_linearize % Jacobian sanity check, saves linearization for reach ``` +Run the reach artifacts (`reachability/`): + +```matlab +reach_operation % linear reach tube for operation-mode LQR +barrier_lyapunov % Lyapunov-ellipsoid barrier cert attempt (sweeps weights) +``` + +**Soundness note:** the current reach tube is the LINEAR model's tube; +it is not yet a sound over-approximation of the nonlinear plant. See +`reachability/README.md` § Soundness status. + ## Prerequisites - Python 3.10+ 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 index 8ce983f..fba6529 100644 --- a/claude_memory/2026-04-17-controllers-linearization-and-first-reach.md +++ b/claude_memory/2026-04-17-controllers-linearization-and-first-reach.md @@ -110,6 +110,28 @@ 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.** +## Soundness, α-drift, saturation — the three "this isn't actually rigorous yet" flags + +Dane caught these during walkthrough. All three now documented in the +relevant file headers so they don't live only in the transcript: + +- **Soundness**: `reachability/reach_operation.m` and + `reachability/README.md` now prominently state the current reach + tube is over-approximation of the LINEAR model, not the nonlinear + plant. Upgrade paths: nonlinear reach (CORA/JuliaReach + nonlinearSys) or linear reach + Taylor-remainder inflation. Either + is thesis-blocking for a real safety claim. +- **α-drift**: `ctrl_heatup.m` header and `plant-model/CLAUDE.md` now + note that the feedback-linearization cancellation assumes exact + α_f, α_c. Real reactors have α drifting (burnup ~20%, boron ~10x, + xenon). Robust treatment = α as interval, parametric reach. +- **Saturation semantics**: `ctrl_heatup.m` header now notes the + sat() is piecewise affine and must be either proven dormant or + modeled as a hybrid sub-mode in reach. Operation-mode LQR is + dormant (trivially); heatup will need explicit treatment. + +Documented in `README.md` top-level under the reach commands too. + ## Loose ends for the next session - Julia reach needs state rescaling. `S = diag(...)` chosen so all diff --git a/plant-model/CLAUDE.md b/plant-model/CLAUDE.md index b7fdf30..a00945f 100644 --- a/plant-model/CLAUDE.md +++ b/plant-model/CLAUDE.md @@ -172,3 +172,25 @@ running example. power-spike overshoot in the simulation. Physical resolution: the reactor should be taken critical in shutdown mode (at ~0.1% power) before DRC transitions to heatup. `main_mode_sweep.m` uses this IC. + +## Robustness caveats (idealized in current artifacts) + +- **α_f, α_c are treated as known exactly.** In reality α_f drifts + ~20% over burnup; α_c spans ~10x across soluble-boron dilution over + a cycle; xenon adds 2-3 $ reactivity on its own timescale. The + feedback-linearization in `ctrl_heatup.m` assumes the controller's + α matches the plant's; if not, the clean `rho_total = Kp*e` + property degrades to `Kp*e + delta*alpha*dT`, and the P term must + absorb the residual. Stabilization still holds but reach analysis + should eventually treat α as a bounded parametric uncertainty. +- **Saturation is a hybrid sub-mode.** The `sat(u, u_min, u_max)` in + `ctrl_heatup.m` is formally piecewise affine. Current reach + treats it as dormant (true for operation/LQR, near-true for the + demo heatup trajectory). A rigorous heatup reach has to model + the saturation regions explicitly. +- **Linear-model reach is not sound for the nonlinear plant.** The + reach artifacts in `../reachability/` use the linearization; the + result is a sound over-approximation of the LINEAR model's reach, + not of the plant's. To upgrade: nonlinear reach directly, or + linear reach + Taylor-remainder inflation. See + `../reachability/README.md` § Soundness status. diff --git a/plant-model/controllers/ctrl_heatup.m b/plant-model/controllers/ctrl_heatup.m index eca79a8..74371a3 100644 --- a/plant-model/controllers/ctrl_heatup.m +++ b/plant-model/controllers/ctrl_heatup.m @@ -13,6 +13,19 @@ function u = ctrl_heatup(t, x, plant, ref) % +0.5*beta guarantees rho_total < beta (below prompt), which in % turn bounds the neutron-kinetics excursion rate for reachability. % +% IMPORTANT for reach analysis: +% The sat() is a 3-mode piecewise-affine system (sat-low / linear / +% sat-high). Under linear reach assumptions it must either be +% (a) proven dormant (u_unsat stays in [u_min, u_max] across the +% reach set — trivial to check, expensive to over-approximate +% tightly), or +% (b) handled as an explicit hybrid automaton nested under the DRC +% mode, with transitions when u_unsat crosses the saturation +% bounds. +% The current reach_operation.m assumes (a) implicitly. Heatup +% reach will need option (b) because u_unsat is near the +0.5*beta +% bound during early ramp. +% % Why no integrator: % Ramp tracking has a structural lag proportional to ramp_rate / Kp_eff. % Acceptable because the DRC exits heatup on a predicate window @@ -20,6 +33,18 @@ function u = ctrl_heatup(t, x, plant, ref) % Adding PI would double-count the intrinsic plant integrator % (thermal mass) and make anti-windup a hybrid transition. % +% IMPORTANT caveat on the cancellation u_ff: +% The feedback linearization works only if the controller's values +% of alpha_f, alpha_c match the plant's. In reality alpha drifts: +% alpha_f ~20% over burnup, alpha_c by ~10x across boron dilution, +% plus xenon. With alpha_true = alpha_nom*(1+delta), the +% cancellation leaves residual reactivity delta*alpha*dT that the +% P term cleans up — stabilization still holds, but the clean +% "rho_total = Kp*e" property is gone. A robust deployment would +% treat alpha as an interval and propagate parametric uncertainty +% through reach (zonotope with parameter generators), or add +% adaptive alpha estimation. Idealized here. +% % Inputs: % t - time [s] % x - state vector (10 x 1) diff --git a/reachability/README.md b/reachability/README.md index 607104f..520391f 100644 --- a/reachability/README.md +++ b/reachability/README.md @@ -2,7 +2,47 @@ Continuous-mode verification for the PWR_HYBRID_3 hybrid controller. -## Status +## Soundness status: APPROXIMATE + +The current `reach_operation.m` result is **not a sound reach tube for +the physical plant**. It is a sound over-approximation of the +*linearized* closed-loop system (A_cl = A - BK around x_op) under +bounded disturbance. The linear model is itself an approximation of +the nonlinear plant (`../plant-model/pke_th_rhs.m`), and that +approximation error is not currently bounded or inflated into the tube. + +Two paths to upgrade to a sound result: + +1. **Nonlinear reach directly** — CORA `nonlinearSys`, JuliaReach + `BlackBoxContinuousSystem`, or equivalent. More expensive but the + honest answer. +2. **Linear reach + Taylor-remainder inflation** — compute an upper + bound on `||f_nl(x, u) - (A x + B u)||` over the reach set (via + Hessian norm estimate on each component of `f_nl`) and inflate the + linear tube by that bound. Less expensive, still rigorous. + +Both are thesis-blocking for any safety claim. Deferred only until +the per-mode plumbing is solid; it is not a "nice to have". + +The current 5-orders-of-margin buffer (reach envelope ~0.03 K against +a 5 K safety band) means linearization error would have to be huge to +invalidate the conclusion, but that is vibes, not a proof. + +## Related open issues + +- **Saturation semantics.** `ctrl_heatup.m` uses `sat(u, u_min, u_max)`. + Saturation is formally a 3-mode piecewise-affine system. For + heatup reach this has to be handled as (a) hybrid locations, or + (b) proven dormant via reach on `u_unsat`. Not modeled in the + current artifacts (operation-mode LQR saturation is dormant in + practice but the proof is implicit). +- **Parametric uncertainty in α_f, α_c.** Real plants have α drift + with burnup (~20%), boron (α_c ranges 10×), xenon. The + feedback-linearization in `ctrl_heatup.m` assumes exact α; a robust + treatment would make α an interval and propagate parametric reach. + Currently idealized — flag in the chapter. + +## What's here **Per-mode only.** Following the compositionality argument in the thesis: verify each continuous mode separately, let the DRC handle discrete @@ -48,8 +88,11 @@ options are CORA (MATLAB) or JuliaReach (port the plant to Julia). ## What this does NOT do yet +- Any sound reach tube (see top of this file). - Nonlinear reach for the original P controller on operation. -- Heatup reach (the ramped reference makes x* time-varying — needs - trajectory-LQR or a different formulation). +- Heatup reach (ramped reference makes x* time-varying — needs + trajectory-LQR or a different formulation, and the saturation + semantics need to be made explicit). - Shutdown, scram, initialization reach. - Hybrid-system level verification (mode switching validity). +- Parametric robustness to α_f, α_c drift. diff --git a/reachability/reach_operation.m b/reachability/reach_operation.m index 11106af..183d046 100644 --- a/reachability/reach_operation.m +++ b/reachability/reach_operation.m @@ -1,9 +1,30 @@ %% reach_operation.m — linear reach set for operation mode (LQR closed-loop) % -% Compute a sound over-approximation of the reach set starting from a -% box around x_op, under LQR feedback, with Q_sg in a specified -% interval. Check that T_avg stays inside the t_avg_in_range predicate -% for all t in [0, T_final]. +% *** SOUNDNESS STATUS: APPROXIMATE, NOT SOUND. *** +% +% This file computes a reach tube for the *linearized* closed-loop +% system (A_cl = A - BK around x_op) under bounded Q_sg. The tube +% itself is a sound over-approximation of the LINEAR model's reach +% set — it uses conservative box hulls and elementwise-absolute-value +% matrix propagation. But the LINEAR MODEL is only an approximation +% of the real nonlinear plant (pke_th_rhs.m), so the result is not a +% sound reach tube for the actual plant. +% +% To upgrade to a sound result, pick one: +% (a) Nonlinear reach directly (CORA's nonlinearSys, JuliaReach's +% BlackBoxContinuousSystem). Expensive, honest. +% (b) Linear reach + Taylor-remainder inflation: compute an upper +% bound on ||f_nonlinear(x,u) - (A x + B u)|| over the reach +% set and inflate the linear tube by that bound. Cheaper, +% requires a Hessian-norm estimate. +% Tracked as a thesis-blocking todo. For now, the 5-orders-of-margin +% buffer (|dT_c| ~ 0.03 K vs safety band 5 K) gives us a lot of room +% to absorb linearization error, but that's not a proof. +% +% Compute a reach-tube over-approximation starting from a box around +% x_op, under LQR feedback, with Q_sg in a specified interval. Check +% that T_avg stays inside the t_avg_in_range predicate for all t in +% [0, T_final]. % % This is the *continuous-mode obligation* for q_operation: % X_entry := { x : |x - x_op| <= delta_entry }